引言

自1966年河北邢台M_S7.2地震后我国将物探中的视电阻率方法引入地震监测中(国家地震局预测预防司, 1998), 并开展了地电阻率定点连续观测，经过50余年的研究和发展，地电阻率已成为目前地震前兆观测的重要手段之一.目前全国80多个地电阻率台站分布在主要的地震活动断裂带和大城市周围，承担着常规的地震监测任务，形成了一定规模和规范化的地电阻率监测台网，对我国大陆地震监测预报起着重要的支撑作用.在地电阻率的长期观测中，这些台站曾多次记录到大陆及周边地区中强地震前的地电阻率异常变化, 如1976年唐山M_S7.8和松潘—平武M_S7.2、1998年张北M_S6.2、2003年大姚M_S6.2、2008年汶川M_S8.0和2016年民乐—山丹M_S6.1地震等(钱复业，赵玉林，1982；桂燮泰等, 1989；钱复业等，1990；高立新等，1999；杜学彬等，2015)，并对其中部分地震实施了一年尺度的预测(叶青等，2005；杜学彬等，2013).已有研究结果和观测事实表明，中强地震前在震中一定范围内, 地电阻率可能出现较为显著的异常现象.但是近年来许多台站受到了严重的环境干扰，削弱了地电阻率对地震监测预测的服务能力.对于地电阻率观测中受到的各类干扰，我国地电工作者作了许多理论研究和实验分析工作，例如: 金安忠等(1990)关于高压输电线路对地电阻率观测的影响分析表明，其干扰强度随输电线路与观测装置距离的增大而逐渐减小，且高压输电线路平行于布极方向时干扰最强，垂直于布极方向时干扰最弱；Lu等(2004)模拟分析了一圆盘介质异常体在二维模型中不同位置时地电阻率观测影响系数的分布，分析结果表明测量电极与供电电极之间的影响系数为负，测量电极之间和供电电极之外影响系数为正；针对金属导线类干扰，张磊等(2010)在二维均匀介质模型下计算了其对地电观测的干扰幅度和形态特征，结果显示干扰幅度由金属导线的长度、方位和位置所决定.这些地表金属管线和大面积灌溉等干扰源会影响布极区表层的电性结构，进而影响地电阻率的观测结果.为了取得优质的观测数据，我国地电工作者一直在对观测方法进行研究和改进，例如：解滔等(2013, 2014)分析了地电阻率观测中不同测道的相反年变现象，认为测区地下介质的不均匀性是造成该现象的主要原因；刘昌谋等(1994)和聂永安等(2010)发展了井下地电阻率观测方法以抑制地表的干扰、降低季节变化等.

延庆台定点地电阻率数字化观测始于1999年，观测数据质量良好，多年来一直受到地震地电阻率预测研究人员的重视，同时也是京津冀地区地电阻率映震能力较好的地震台站之一. 闫睿等(2015)利用归一化月速率方法对北京及邻区近30年内发生的7次中强地震进行了分析，讨论了该台地电阻率异常与地震的对应关系，得出延庆台地电阻率对于测区附近中等以上地震和远距离大震具有一定的映震能力.之后的延庆台地电阻观测表明，测区出现的干扰源对观测资料的质量产生了相当程度的影响，削弱了该台地电阻率资料在地震趋势预测中的作用.因此，本文拟采用三维有限元数值分析方法，分析近年来测区内地下铁质干扰源对地电阻率观测的影响幅度和形态，为今后延庆台地电阻率资料异常变化的进一步跟踪分析和判定提供参考.

1.   台站简介

延庆地震台位于北京市延庆县城关镇张庄村南(图 1)，是首都圈地区重点监测的台站之一.该台站地处延庆—怀来盆地的边缘，地质构造属于燕山褶皱带的一部分，在台站NW向大约3 km处为狼山—黄柏寺断裂带.地电阻率观测区域的第四系地层厚度为329 m，下伏侏罗系火成岩，但沉积层分布极不均匀，主要为黏土、粉黏土、砂土、砾石等互层分布.台站海拔高度为484 m，属于山前地带，地势平坦.测区位于农田内，无山川、河流，亦无大型厂矿、公路及高压线路.

图  1  延庆台及北京周边断层分布

Figure  1.  Fault distribution at Yanqing station and around Beijing

下载: 全尺寸图片 幻灯片

延庆地电阻率观测采用对称四极装置(图 2)，分别布设NS和EW两个测道，两测道供电极距AB均为1.5 km，测量极距MN均为0.5 km.台站观测装置、测量仪器、外线路和电极等均符合地电阻率台站观测的技术要求(中国地震局，2001)，观测仪器为DDC-2系列和ZD8系列地电仪.

图  2  延庆地电阻率观测线路布极图

Figure  2.  Electrodes arrangement for resistivity observation at Yanqing station

下载: 全尺寸图片 幻灯片

2003年以前观测装置以共中心点方式布设(图 2a)，由于受到当地环境干扰(养鸡场铁质地下管线)，观测数据年变形态发生变化，观测质量下降；2003年底进行了观测装置改造，将NS测道的观测线路向西平移(图 2b)，改造后观测资料质量得到了一定程度的提高；2008—2011年，测区内出现两次铁丝网接地干扰地电阻率观测的情况，因铁丝网均平行于NS方向排列，故NS测向观测数据变化较大(王同利等，2016)；2011年年底又对测区内台站周围的线路、观测仪器和电极进行了更新改造，至2012年3月，正式更换为新观测系统，观测质量得到进一步提高. 图 3为1995—2016年延庆地电阻率观测曲线，其清楚地反映了观测数据的真实变化情况，包括环境干扰、观测系统改造及装置系统故障对观测数据的影响，是地震前兆监测分析的基础.延庆台地电阻率对周边地震异常反应的灵敏度较高，自2000年以来，台站周边300 km内发生M_S4.0以上地震34次，M_S5.0以上地震4次.最近一次地震为2014年9月6日河北省涿鹿县M_S4.3地震，虽然震级不大，但影响范围较广，北京、天津、河北等地均有震感.延庆台EW测道观测到的地电阻率自2013年年底开始呈下降趋势，2015年7—8月延庆台地电阻率受到台站院内房屋装修使用的铁脚手架及EW测道东侧400 m处铺设楼地基的干扰，2016年NS测道观测的地电阻率也呈下降趋势(图 3).为了分析延庆地电阻率持续下降的原因，对上述干扰事件进行了数值模拟分析.

图  3  1995—2016年延庆台地电阻率观测数据月均值曲线

Figure  3.  The monthly mean value curves of georesistivity observation data at Yanqing station from 1995 to 2016

下载: 全尺寸图片 幻灯片

2011年进行台站改造时，测出了延庆地电阻率测区电测深曲线，结果如图 4所示.为了更好地利用三维有限元模拟计算干扰源的影响情况，研究人员于2014年对测区高密度电性剖面进行了测量，结果显示第四系沉积层的电阻率较低，一般在20—36 Ω·m范围内变化，局部地段略高些，最高可达45 Ω·m(图 5)，电阻率分层不明显, 大致可分为3层.

图  4  延庆观测场地电测深曲线

Figure  4.  Electrical sounding curves of Yanqing station

下载: 全尺寸图片 幻灯片

图  5  延庆观测场地高密度电法测量图

Figure  5.  High density electrical measurement at Yanqing station

下载: 全尺寸图片 幻灯片

2.   有限元数值分析

2.1   稳恒电流场有限元方法

有限元方法将求解域划分为一系列小单元，其最大优点是通过采用小单元最大程度地逼近几何模型.地电阻率定点观测采用对称四极装置，在观测过程中，供电电流I自供电电极A，B流入地下介质，在测量电极M，N两端测量电位差ΔU，由

计算出地电阻率值(K为装置系数).该观测系统可视为稳恒电流场进行计算，当直流电流固定时，电位满足泊松分布，即

式中，U为由直流电流源I产生的电位，σ为介质电导率，δ(x, y, z)为狄拉克δ函数(Ida, Bastos，1997).由于总电位U在电流源处存在奇异性，用数值计算方法求解式(2) 时在电流源附近得到的结果误差较大，常用的解决方法是将总电位分解为在均匀介质σ₀中产生的电位U₀和在非均匀介质σ_a中产生的电位U_a，且满足U=U₀+U_a，σ=σ₀+σ_a.因此等效变分问题的能量函数为

式中，F(U_a)为变分能量函数，v为计算区域，Γ_∞为模型除自由界面外的边界，Γ为全部边界，包含Γ_∞和地表自由界面，r为电流源的径向矢量，n为边界上的外法线矢量(Huang, Lin，2010).系统稳定的条件为能量函数的变分应满足

有限元方法是通过变分原理使能量函数F(U_a)达到最小，而非直接求解方程.因在式(3) 中无电流源，从而消除了式(2) 中的奇异点.直流电流源I在均匀介质中的电位分布U₀由解析解算出，在非均匀介质中的电位分布U_a通过有限元方法计算得到(徐世浙, 1994；Wu, 2003).地表的地电阻率观测满足纽曼(Neumann)边界条件，有限元模型的其余边界既可以施加狄利克雷(Dirichlet)边界条件，也可以施加纽曼边界条件(Coggon，1971).但是对于一个尺寸固定的模型，当观测供电极距大于某个定值后，施加狄利克雷边界条件计算得到的理论值小于实际观测值，而施加纽曼边界条件计算得到的理论值则大于实际观测值(Dey, Morrison，1979；Li, Spitzer，2005).有限元模型尺寸越大，边界条件对计算结果产生的影响越小，计算量则越大，因此需要合理地估计模型的尺寸.模型建立后经过单元离散化、施加电流源和边界条件，再对每个单元的电位进行有限元数值计算，得出测量电极之间的电势差，最后通过装置系数计算出地电阻率值(Huang, Lin，2010；解滔, 卢军，2015).

2.2   模型建立

针对延庆台地电阻率的观测情况建立稳恒电流源对称四极观测的三维有限元模型.首先根据延庆台测区内电测深和电性剖面的测量结果，将测区地下电性结构分为3层(图 5):上层ρ₁＜24 Ω·m，厚度为17 m；中层28 Ω·m＜ρ₂＜38 Ω·m，厚度为80 m；下层ρ₃＜28 Ω·m，厚度为3000 m.边界条件采用纽曼方式.为保证理论值与计算值具有较好的一致性，建立不同水平尺度和不同底层厚度的多个模型来对比分析计算值随模型尺寸的变化，最后选定模型的水平尺寸为7AB，底层厚度为2AB，如图 6所示.为了检验此有限元模型计算结果的准确性，将数值计算结果与理论解析解进行对比分析，结果如图 7所示，可以看出, 理论解析解与模型数值计算值有较好的一致性.

图  6  延庆台地电阻率有限元模型示意图

Figure  6.  The three-dimensional finite element model of georesistivity at Yanqing station

下载: 全尺寸图片 幻灯片

图  7  有限元模型计算值和解析解曲线

Figure  7.  The numerical results and analytic solution curve of the finite element model

下载: 全尺寸图片 幻灯片

测区内铁质干扰介质由铁管组成，因此，在模型中使用导线单元，计算中只需导线的横截面积和电阻率值，而导线单元的长度则在模型建立时设定为固定值.模型中其它部位为三维立体单元，单元网格划分时布极区较密，向外及深度方向逐渐变稀.经过对模型网格划分的多次调试比对，发现深度方向每层网格划分为5个以上单元时，理论值与计算值基本接近，因此，在模型网格划分时保证每层为5个以上单元，这样既保证了计算精度，又可节省计算时间(王同利等，2016).

地电阻率观测中存在年变化现象.已有研究表明，年变化是由浅层介质的电阻率随季节性降雨(浅层介质含水率变化)及温度变化引起的(张学民，2004).夏季降雨量增加，温度上升，浅层介质电阻率下降；冬季降雨量减少，温度降低，浅层介质电阻率上升，因而地电阻率观测值也表现出与季节性降雨和温度变化一致的年周期性变化.尽管延庆台的年变化幅值很小，但变化形式相同.因此假设模型第一层电阻率变化而厚度不变，其余层位电阻率和厚度均不变，即变量为第一层电阻率.

本文分析的铁质干扰介质的空间位置如图 6所示，其中：

1) 养鸡场地下铁管(暖气管，水管并排)始建于1999年4月，至2000年年底完全建成，铁管距离EW测线垂直距离约为280 m，EW向长为420 m，直径为0.01 m，埋深为0.2 m，壁厚为0.005 m.干扰数据出现于铁管始建时，之后持续变大，直至2003年12月，NS测向线路向西移动约300 m后干扰才消失，至2004年观测数据变化正常；

2) 2015年7月17日，延庆地震台站院内进行房屋装修，在水泥地上搭建铁脚架，形状为正方形，大小为25 m×25 m，脚手架位于EW测线附近，距离西端测量极M₃ 240 m，距离北端测量极N₁ 310 m；

3) 铁质地基建于2015年8月中旬，位于延庆地电阻率EW测线的东端，距离供电极B₃约400 m，地基大小约为200 m×200 m.

3.   计算结果

3.1   养鸡场铁质管线干扰

首先对养鸡场地下铁管进行数值模拟分析.模型第一层介质的电阻率获取方式为：采用1998年地电阻率观测数据的月均值，逐步修改模型第一层介质的地电阻率，直至模型计算值与观测值的误差小于仪器观测误差时，将该计算值作为该月均值对应的第一层介质的电阻率，并将由12个月均值计算得到的电阻率作为第一层介质电阻率的年变化值.严格地讲，2002年第一层介质电阻率的年变化与1998年不可能完全相同，但由于2002年受到养鸡场建设的干扰，因此用1998年第一层介质电阻率的年变化代替2002年的数据，据此计算有铁管干扰情况下地电阻率随第一层介质变化时的变化情况，便于分析延庆台地电阻率观测受金属管线干扰的动态变化特征. 图 8给出了2002年延庆地电阻率观测受观测地内养鸡场铁质管线影响的有限元分析和观测系统改造前后的数据变化情况，表 1列出了有限元模拟得到的NS和EW测向干扰源的年变化量.

图  8  延庆台地电阻率观测数据及有限元分析

Figure  8.  Georesistivity observation data and finite element analysis for Yanqing station

下载: 全尺寸图片 幻灯片

表  1  2002年延庆地电阻率有限元模拟计算得到的干扰源年变化量

Table  1.  Annual variations of interference sources calculated from FEM simulation for georesistivity of Yanqing station in 2002

月份	干扰源年变化量
	EW测向/(Ω·m)	NS测向/(Ω·m)
1	-0.236	0.111
2	-0.242	0.070
3	-0.227	0.165
4	-0.283	-0.144
5	-0.273	-0.104
6	-0.249	-0.054
7	-0.231	-0.025
8	-0.249	0.026
9	-0.246	0.049
10	-0.242	0.070
11	-0.227	0.165
12	-0.239	0.091

下载: 导出CSV 

| 显示表格

从图 8和表 1可以看出，1998年延庆台地电阻率观测数据的年变化幅值很小，EW测向年变幅度为0.1 Ω·m, NS测向年变幅度小于0.01 Ω·m, 几乎可以忽略；而2002年年变幅度增大，且数据变化不符合季节变化规律，可能由测区内铁管异性体干扰所致；2003年NS测线西移后，2004年的年变形态基本恢复为受干扰前(即1998年)的变化形态.同时，三维有限元数值模拟得到的2002年的数据变化与2002年实际观测的数据变化基本接近，无干扰下的计算值基本接近于测区无铁管干扰的实际观测数据变化，说明此模型真实可靠地反映了延庆台地电阻率观测区环境的变化情况.因此，此模型可以作为延庆台地电阻率观测区的环境模型来处理数据异常变化，分析铁质干扰导致的数据变化，为鉴定前兆异常和干扰异常奠定基础.

3.2   铁脚架干扰

2015年7月，延庆台地电阻率出现异常变化，尤其是EW测向变化幅值较大，干扰源的具体位置见图 6.经实地异常落实，发现7月17日延庆台站院内进行房屋装修时搭建了铁脚架.在进行三维有限元数值模拟分析时，考虑到异常变化时间短，对2015年6—9月的数据采用日均值进行计算，其它时间采用月均值.数值分析结果反映出2015年7月17日台站安装的铁脚架对EW测向的影响量为-0.02 Ω·m，低于正常观测值，但相对变化量很小，仅0.07%左右；对NS测向的影响为0.009 Ω·m，高于正常观测值，但变化量很小，考虑到误差的可能性，认为数值模拟曲线与实际观测曲线变化基本一致.该结果与安装铁脚架后EW测向短期内(7月中下旬)的变化基本一致.数值分析结果如图 9a所示(NS测向异常变化可忽略)，表 2列出了NS和EW测向的干扰年变化量.

图  9  铁脚架(a)和铁质地基(b)干扰的有限元分析

Figure  9.  Finite element analysis for three-dimensional structures of the iron frame (a) and the iron foundation (b)

下载: 全尺寸图片 幻灯片

表  2  铁脚架有限元分析的干扰量

Table  2.  Disturbances achieved from finite element analysis of iron frame

月份	干扰源年变化量
	EW测向/(Ω·m)	NS测向/(Ω·m)
2014-10	0.0017	-0.0004
11	0.0019	-0.0005
12	0.0023	-0.0007
2015-01	0.0027	-0.0007
02	0.0023	-0.0007
03	0.0019	-0.0005
2015-04	0.0018	-0.0005
05	0.0016	-0.0005
06	-0.0035	-0.0004
07	-0.0218	-0.0003
08	-0.0195	-0.0002
09	-0.0292	-0.0002

下载: 导出CSV 

| 显示表格

3.3   铁地基干扰

2015年8月中旬，延庆台地电阻率观测EW测向数据持续下降，且变化幅值相对增大，经调查落实，发现延庆台地电阻率观测EW测线东端，距离供电极B₃约400 m处在建楼施工，地基建造物为铁质，其具体位置见图 6.利用三维有限元数值模拟方法进行模拟计算，结果如图 9b所示(EW测向)，表 3列出了NS和EW测向的干扰年变化量.

表  3  铁地基有限元分析的干扰量

Table  3.  Disturbances achieved from finite element analysis of iron foundation

月份	干扰源年变化量
	EW测向/(Ω·m)	NS测向/(Ω·m)
2014-10	-0.0049	-0.0029
11	-0.0016	-0.0013
12	-0.0054	-0.0015
2015-01	0.0055	-0.0013
02	0.0013	-0.0093
03	-0.0019	-0.0069
2015-04	-0.0022	-0.0004
05	0.0044	-0.0012
06	-0.0052	-0.0006
07	-0.0332	-0.0023
08	-0.0540	-0.0008
09	-0.0425	-0.0033

下载: 导出CSV 

| 显示表格

数值模拟得到的铁质地基对延庆台地电阻率EW测向的影响量为-0.04 Ω·m，即相对变化量约为0.15%. NS测向的影响值小于0.0001 Ω·m，可认为铁质地基对NS测向未造成干扰. EW测向的变化趋势与2015年8—9月实际观测值的变化趋势基本相同，但变化量较实际观测变化量小.

4.   讨论与结论

在将三维有限元数值模拟方法运用于地电阻率分析方面, 前人已积累了很多实例，例如解滔和卢军(2015)总结了地电阻率观测中三维影响系数的分布情况，计算结果显示，当干扰介质位于测区不同位置时，其对观测数据的影响系数不同.本文采用三维有限元模型对延庆台地电阻率测区内铁管、铁脚架和铁质地基对观测数据的影响进行了数值模拟分析，计算结果与实际观测数据变化一致，反映出干扰源位置不同，即相对测量电极和供电电极的距离和位置不同，则对实际观测的影响幅度不同.同时，铁质干扰介质的排列方向与测线方向也密切相关：当铁质干扰介质与观测测线方向一致时，干扰变化较大，当铁质干扰介质与观测测线方向垂直时，干扰变化较小，即铁脚架和铁质地基对EW测向影响较大，而对NS测向影响很小.养鸡场地下铁管对两个测向均有影响，但对NS测向的影响程度远大于EW测向.

文中铁质地基干扰对延庆台地电阻率影响的数值分析结果比实际观测数据变化偏小，原因可能为数值分析模型与测区地下实际电性结构之间存在一定的差异，尤其是延庆台地电阻率年变较小，年变幅值仅为0.17 Ω·m左右，可见其观测场地地下电性结构层状不明显，分层分析存在一定的局限性.同时此分析结果也说明，对于延庆台地电阻率细小、异常的变化应进行分析研究，因其反映了地下构造的明显变化.计算结果显示2015年铁脚架干扰对延庆台地电阻率下降无影响，铁质地基的影响并不明显，说明延庆地电阻率的同期下降变化不是由干扰引起的.文中的干扰幅值和形态分析结果与实际观测资料基本吻合，从数值分析的角度确认了测区内已有干扰源对地电阻率观测的影响，为今后该台地电阻率资料的跟踪分析和异常性质的判定提供了一定的参考.