0 引言

经过50多年的发展，地电阻率定点连续观测已成为地震地球物理观测的重要测项之一。在多次中强地震前，特别是在1976年唐山7.8级、2008年汶川8.0级地震前，均记录到显著的地电阻率异常。但是，目前地电阻率台站受到的干扰愈来愈多，且受干扰情况较复杂，特别是地铁运营、城际高铁、高压输电线路、管道输送等大型工程建设产生的干扰，严重影响了地电台站观测数据产出质量；另外，还有水渠、大型铁质用具、金属管道布设、铁丝网等金属管网及工业、农业游散电流干扰，这些干扰均会引起地电阻率测值的变化。

地电阻率观测数据用于地震预测时，首先要解决干扰问题。解滔等（2015）、刘素珍等（2017）基于影响系数理论，采用数值模拟方法对地震电阻率观测受局部介质电阻率变化的影响机理进行了研究；解滔等(2015, 2016)依据三维影响系数在地表的分布情况定性判断干扰源对地电阻率观测的影响形态，并结合不同测道测值排除或锁定某些干扰源；解滔等（2015）、王同利等（2017）采用有限元数值计算方法建立三维模型，分析了地电阻率测区位于地表的金属导线和局部异常体对观测产生的干扰形态和幅度随时间的变化特征。本文利用浚县地震台测区铁质异常体相关信息及地质剖面、岩层电性、电测深等资料，建立三维有限元模型，经模型检验后，结合实际观测定量地分析局部铁质异常体对观测的影响幅度和形态，以期为进一步跟踪浚县地震台地电阻率观测数据异常变化和地震前兆异常判定提供量化依据。

1 浚县地震台简介

浚县地震台位于位于河南省浚县白寺乡中国人民解放军73569部队农场东部。该台地处鲁西隆起、太行山隆起及开封断块之间，构造单元以块状结构为主，凹凸相间。该区域的濮阳断块为NE向菱形块体，其北界为磁县—大名断裂，南界为新乡—商丘断裂，西有太行山山前断裂与太行山隆起接壤，东为聊兰断裂与鲁西隆起分界。区域地层受青羊口断裂、汤西断裂及EW向的安阳南断裂控制（图 1），浚县地震台距活动构造汤东断裂仅25 km左右。台站海拔58 m，测区地貌平坦开阔，坡度高差极小，周围主要为农田，局部为林地，植被发育，环境干扰源较少。岩性表层360 m以上为粘土层，以下为新近系泥灰岩。地表大部分为全新统、上更新统地层所覆盖，而且第四纪覆盖层仅185 m，观测环境较有利。

浚县地震台地电阻率观测系统于2013年11月30日建成，2014年初投入运行，观测数据质量符合规范要求，年变化形态正常。台站地电阻率观测布极方式采用对称四极法（图 2），共布设NS、EW两个测道，观测室偏离布极区中心的距离为180 m。供电电极极距均为1 km，测量电极级距均为0.3 km，电极均采用900 mm×900 mm×8 mm的铅板电极，外线路采用地埋方式，电极埋设深度3.0 m，观测仪器为ZD8M型地电仪。

在日常数据监控中发现，自2016年5月13日起EW测道地电阻率数据出现明显的下降阶变，之后趋于稳定。经环境调查发现，2016年5月12日，浚县地震台地电阻率EW测线东供电极东端，新建1条三角铁轨，该铁轨NS向长约300 m，宽约0.05 m，距东供电极约9.2 m。由于异常出现时间与铁轨铺设时间较接近，因此，EW测道的数据变化是否与铁轨有关，以及铁轨对地电阻率观测的影响需要进一步验证。

2 有限元数值分析 2.1 稳恒电流场有限元法

地电阻率定点观测中采用对称四极装置，观测过程中，在供电电极A、B输入直流电流，在测量电极M、N测量电势差。此问题可视为稳恒电流场计算，满足Maxwell方程组和电荷守恒定律，因此稳恒电流场问题可表示为如下泊松方程

$ \nabla \cdot (\sigma \nabla V) = - I\delta (x, y, z) $

(1)

式中，V为由电流源产生的电位；σ为介质电导率；δ(x, y, z)为Dirac delta函数。有限介质空间的全部边界为Г。其中，一部分边界没有电流流出（如地表），满足Neumann边界条件，记为Г_S；其余边界记为Г_V，满足Dirichlet边界条件。因此，式（1）满足如下边界条件

$ \mathit{\Gamma } = {\mathit{\Gamma }_{\rm{S}}} + {\mathit{\Gamma }_{\rm{V}}} $

(2)

其中，${\left. {\frac{{\partial V}}{{\partial n}}} \right|_{{\mathit{\Gamma }_{\rm{S}}}}} = 0$，n为边界指向区域外的法线方向；${\left. V \right|_{{\mathit{\Gamma }_{\rm{V}}}}} = p$。

应用虚功原理可得到稳恒电流场泊松方程的有限元弱解形式

$ \int_\mathit{\Omega } \sigma \left[ {\frac{{\partial V}}{{\partial x}} \cdot \frac{{\partial \varphi }}{{\partial x}} + \frac{{\partial V}}{{\partial y}} \cdot \frac{{\partial \varphi }}{{\partial y}} + \frac{{\partial V}}{{\partial z}} \cdot \frac{{\partial \varphi }}{{\partial z}}} \right]{\rm{d}}\mathit{\Omega } = \int_\mathit{\Omega } I \delta (x, y, z)\varphi {\rm{d}}\mathit{\Omega } $

(3)

式中，Ω为计算区域；φ为任意的虚位移函数，在满足Dirichlet边界条件的边界上，虚位移函数φ = 0。

地电阻率观测在地表自然地满足Neumann边界条件，在水平方向和垂直方向（深度）可视为无穷远边界，可以施加Dirichlet边界条件（V = 0），也可以施加Neumann边界条件（Coggon，1971）。但是，建立的模型在水平、垂直方向上的尺度不可能是无限的，对于一固定尺寸的模型，在供电极距AB大于一定的值后，对无穷远边界施加Diriehlet边界条件时计算得到的地电阻率值将小于实际值；而对无穷远边界施加Neumann边界条件时计算得到的地电阻率值将大于实际值（Dey et al，1979；Li et al，2005）。对固定的供电极距AB，模型尺寸越大，边界对计算结果的影响就越小。但是模型越大，计算量也就越大，因此需要合理地选择模型水平方向的尺寸和最底层厚度。

模型经单元离散化和施加电流源、边界条件后可对单元节点上的自由度（电位）进行数值求解，求解出电位分布后可以获得测量电极间的电位差，进而依据对称四极装置系数计算地电阻率（Huang et al，2010）。

2.2 模型建立

根据浚县地震台测区内电测深和电性剖面的测量资料，该台测区场地电性结构横向分布较均匀，具有明显的水平层状结构的特点。该区域第四系沉积层主要由粉土和泥灰岩组成，粉土层电阻率为10—40 Ω·m，泥灰岩层电阻率为9—57 Ω·m，总体电阻率较低。

依据浚县地震台EW向电测深曲线，反演得到测区电测深曲线（图 3），将其作为三维有限元模型的基础数据，建立水平层状模型并假定浅层3层介质电阻率均匀，参考地电阻率观测资料得到测区地下电性结构，具体数值见表 1。一些研究者（解滔等，2013；王同利等，2017）通过对多个台站地电阻率观测资料进行有限元数值分析认为，模型水平尺寸越大，地层越厚，计算结果就越接近理论值，同时计算量也越大。当模型水平尺寸D＞6倍AB、模型厚度H＞2倍AB时，计算结果越接近理论值。为了保证计算结果与理论值间具有较好的一致性，建立一定水平尺寸、不同厚度的模型。经过多次试验，最终确定模型尺寸为长6 000 m，宽6 000 m，厚3 000 m，模型水平面中心始终位于台站观测装置布极中心，此时边界效应对计算的影响已小于仪器的观测精度，可满足模型计算的要求。

设测区内铁轨为长方形金属铁板，长为300 m，宽为0.05 m，计算铁板的面积，在模型中使用导线单元对铁板划分网格，从里到外依次采用10×10、20×20、50×50进行有限元网格划分，图 4分别为铁轨所在位置有限元模型水平剖面和铁轨所在水平面中心布极区网格划分剖面。

2.3 计算结果

表 2为铁轨铺设前后浚县地震台地电阻率EW向观测值与计算值。由表 2可见，铁轨铺设前后地电阻率观测值变化为-0.04 Ω·m，有限元模型计算变化值为-0.08 Ω·m，此结果与实际观测值的变化趋势基本相同，但变化量较实际观测值大。铁轨对观测值的影响趋势与实际观测值的下降同步，这主要是由局部低阻体的影响所致，与解滔等（2016）低阻干扰源影响所表现出的电阻率下降的研究结果相符。所以，2016年5月浚县地震台地电阻率EW测向测值下降变化与铺设的铁轨有关。

3 结论

采用三维有限元模型对浚县地震台地电阻率测区内铁轨对观测数据的影响进行了数值模拟分析，得到以下结论。

（1）三维有限元初始模型的建立较关键。在建立模型时发现，当对铁质异常体分别采用导线单元和铁板进行网格划分时，计算所得结果有所不同，铁板模拟比导线单元模拟计算结果明显偏低，更接近实际观测值，但不同模型的计算结果变化形态和趋势一致。计算表明，铁质异常体产生的干扰与其相对电极的位置（王同利等，2016）和自身长度间也有关系，当低阻体位于供电极附近时，干扰形态为下降阶变，并且距电极越近，干扰幅度越大；当干扰源长度不同时，产生干扰的幅度也不同。

（2）由有限元数值计算结果发现，浚县地震台东供电极附近的铁轨会导致地电阻率观测值的下降，数值分析结果比实际观测数据变化偏大，这可能是数值分析模型与测区地下实际电性结构之间存在一定差异，而干扰幅值和形态分析结果与实际观测结果趋势一致，则主要是由局部低阻异常体的影响所致。因此认为，浚县地震台地电阻率EW测向测值出现的下降阶变是由环境干扰引起的。