自1966年邢台M_S7.2地震后我国开始进行大规模规范化的地电阻率观测， 在40多年的观测中报道了许多大震、 中强震前突出的地电阻率中期—短临阶段的异常(钱复业等， 1982， 1990； 桂燮泰等， 1989； Lu et al, 1999； 张学民等， 2009； 杜学彬， 2010)．在长期的观测中， 发现几乎所有的台站都存在季节性的年周期变化(简称年变化或年变)， 对年变的分析研究表明， 年变主要是由地电阻率测区表层介质随季节变化而产生的(金安忠， 1981； 王志贤， 1981； 钱家栋等， 1985； 刘允秀等， 1999； 张学民， 2004)．直观地看， 在降水多的季节， 表层介质含水率升高， 电阻率降低， 引起地电阻率观测值降低； 在降水少的季节， 表层介质电阻率上升， 地电阻率观测值升高. 这是我国大多数地电阻率台站的年变形态， 称之为“正常年变”．同时也发现有一部分台站所有测道地电阻率在降水多时观测值升高， 而降雨少时观测值降低， 与正常年变截然相反．结合台站的电测深资料， 采用水平层状均匀介质模型, 计算了各层介质电阻率变化对地电阻率的影响系数. 结果表明， 这些台站表层介质电阻率的影响系数为负(钱家栋等， 1985； Lu et al， 2004)， 即表层介质电阻率升高， 地电阻率降低. 反之亦然．

除上述两类变化外， 还存在这样一类年变化， 同一台站不同测道地电阻率年变形态相 反， 即“相反年变化”， 本文所要分析的新沂台就是其中之一．新沂台地电阻率EW和N45°E 测道地电阻率年变形态为表层介质电阻率降低时观测值降低， 表层介质电阻率升高时观测值升高; 而NS测道年变形态则在表层介质电阻率降低时观测值升高， 表层介质电阻率升高时观测值降低. 地电阻率台站在建台时一般沿相互垂直的两个方向进行电测深测量， 如果测区地下介质大致是水平的、 均匀的分层介质， 则两个方向的电测深曲线类型一般相近， 与水平层状均匀介质模型反演时所反映的电性剖面大体相同， 因而表层介质电阻率随季节变化时不同测道(极距相同)地电阻率年变形态总体上一致； 如果两个方向的电测深曲线存在显著差异， 说明测区地下介质在不同方向上的不均匀性显著， 不能用单一的层状介质模型来近似电性剖面， 因为采用单一水平层状均匀介质模型难以对同一台站不同测道地电阻率年变形态相反这一现象做出合理的解释． 本文结合新沂台电测深和测区地质剖面资料， 建立三维空间非均匀介质有限元模型， 计算3个测道地电阻率随表层介质电阻率变化的年变形态. 计算结果符合该台3个测道实测的年变形态. 其结果对分析、 研究地电阻率台站年变化有启示意义．

1.   新沂地电阻率台站简介

新沂地电阻率台位于苏鲁交界地区、 新沂市市区东郊的闻马村， 地处郯庐断裂中南段. 台站东侧约1 500 m处为山左口—王庄集断裂F₁， 台站西侧为近NS走向的大贺山—桥北镇断裂F₅， 紧邻F₅断裂东侧为近NS走向的f₅断裂， 且F₅和f₅断裂经过地电阻率测区(图1a)．F₅断裂是第四纪以来活动表现强烈的一条断裂. 该断裂由两条断裂组成， 新沂官庄以西为接合部位， 北支断裂由新沂何庄向北经马陵山西侧的山外岩大贤庄向北延伸， 南支断裂由神山向南经马陵山、 宿迁地区， 断续出露^①．台站测区地势平坦， 地下水潜水位约3 m. f₅断裂以东地表覆盖层较浅， 基岩埋深约2 m; f₅断裂以西过渡至F₅覆盖层厚度介于2—6 m之间. F₅以西基岩埋深迅速增至70—80 m左右(张秀霞等， 2009)．

江苏省地震局. 2011. 江苏省电磁观测台网地震预报效能评估报告．

图  1  新沂台布极区断层、 地电阻率布极(a)和地电阻率日均值观测曲线(b)图

Figure  1.  (a) Schematic diagram of faults and three Schlumberger monitoring arrays used at Xinyi station; (b) Daily average observations of georesistivity at Xinyi station

下载: 全尺寸图片 幻灯片

新沂台地电阻率观测于1978年9月完成建设并投入使用， 目前使用ZD8B数字地电仪．地电阻率观测布设EW、 N45°E和NS方向3条测道(图1a)， 采用对称四极观测装置， 3个测道沿布极中心点对称分布， 布极中心位于F₅与f₅断裂之间， 其中EW和N45°E测道跨越了F₅和f₅断裂， 3条测道供电电极极距AB均为1 000 m， 测量电极极距MN均为240 m．自观测以来台站周围发生了1995年山东苍山M_S5.2地震， 震中距为77 km， 震前新沂台地电阻率记录到了与之对应的地震异常(李飞， 陈有发， 1999； 冯志生等， 2004)， NS测道观测值上升， EW和N45°E测道则出现观测值下降异常．

2.   资料分析

图1b是新沂台3个测道2008年1月—2012年5月的地电阻率日均值观测曲线． 由图可见, EW测道年变形态清晰， 且年变幅度最大， 夏秋季降水量增加时观测值处于低值， 而冬季和春季降水量减少， 观测值处于高值． 新沂台地处苏北， 降水量集中在夏秋季， 因而年变形态中观测值处于高值状态的时间要多于观测值处于低值的时间．N45°E测道年变形态与EW测道在同步性上有差异， 但整体上仍是降水量增加时观测值处于低值， 降水量减少时观测值处于高值．NS测道地电阻率年变形态则与N45°E测道年变形态时间较同步. 在EW、 N45°E两测道观测值处于高值时， NS测道观测值则处于低值. 因此可以认为, NS测道的年变形态与EW、 N45°E测道年变形态相反． 2007年10月—2008年5月观测资料受到测区内蔬菜大棚的干扰， 2008年5月对蔬菜金属支架和铁丝做绝缘处理后基本消除了蔬菜大棚的影响^① . 同时该台2007年之前的观测资料中 NS测道年变形态也与EW、 N45°E 测道年变形态相反．

钱家栋． 2008． 江苏新沂和安徽蒙城地电阻率异常落实汇报．

新沂台地层结构沿NS方向基本一致， 沿EW方向F₅断裂两侧基岩埋深差异较大．图2a给出了EW测道F₅和f₅断裂附近的地层剖面， f₅以东第四纪覆盖层较浅， 基岩以砂岩为主．NS、 EW方向实测电测深曲线如图2b， c所示， 电测深测量采用对称四极装置．表1显示了以水平均匀层状模型对两条实测电测深曲线进行反演得到的电性分层结构. 表中EW向和NS向电性结构在第二层和第三层介质的厚度上存在很大的差异， 说明该台地电阻率测区不同方向的地下介质非均匀性显著．

图  2  (a) 新沂台地电阻率EW向F5和f5断裂附近地质剖面； (b) NS向 电测深曲线实测值与模型计算值； (c) EW向电测深曲线实测值和模型计算值； (d) 由NS向电测深曲线反演的各层介质的影响系数B， 虚线表示负值

Figure  2.  (a) Geologic section along EW observation line around faults F₅ and f₅; (b) Observed and modeled electric sounding data in NS direction; (c) Observed and modeled electric sounding data in EW direction; (d) Sensitivities B of each layer derived from electric sounding data along NS direction where dashed lines denote negative sensitivities

下载: 全尺寸图片 幻灯片

表  1  新沂台电测深曲线反演的电性结构

Table  1.  Electrical structure obtained by inverting electric sounding data at Xinyi station

分层	EW		NS
	ρ/Ω · m	h/m	ρ/Ω · m	h/m
1	17.5	2.1	18.5	2.7
2	59.4	15.1	71.5	32.7
3	104.7	407.4	120.8	108.2
4	68		55
注： ρ为电阻率， h为各层的厚度.

下载: 导出CSV 

| 显示表格

NS向电测深测点位于台站观测室以东300 m处. 该测点距F₅断裂较远， 可将其地下介质电性剖面近似为水平层状均匀介质电性剖面. 则由NS向AK型电测深曲线可知, 测点下介质可简化为4层地电断面(表1)， 并以此地电断面作为模型中F₅断裂以东的初始地层电性结构. 在水平层状模型下计算该地电断面各层介质电阻率变化对地电阻率观测的影 响系数(O′Neill， Merrick， 1984; 钱家栋等， 1985)， 如图2d所示. 在供电极距AB/2=500 m 时表层介质电阻率变化对观测的影响系数为负. 在表层介质电阻率降低时地电阻率观测值升高， 表层介质电阻率增大时地电阻率观测值降低， 这符合新沂台地电阻率NS测道的实测年变形态．EW向电测深测点在观测室以西300 m处， 位于F₅断裂东侧地电阻率观测装置布极中心点附近， 并靠近F₅断裂. EW向电测深曲线与NS向电测深曲线之间有一定的差异．由于F₅断裂EW向两侧基岩埋深差达70—80 m左右， EW向电测深观测跨越F₅断裂， 不宜将EW向电性剖面近似为单一的水平层状均匀介质模型， 因此将地电阻率测区简化为非均匀水平层状模型， 在模型的东盘与西盘采用不同的层参数．显然这样一个简化的非均匀层状模型与测区实际电性结构之间是有一定差别的， 仅是在一定程度上的近似．

3.   有限元数值分析

新沂台NS测道地电阻率年变形态与EW、 N45°E测道年变形态相反， 虽然依据台站以东300 m处的NS向电测深资料, 及水平层状均匀介质模型所反演的电性结构表层介质电阻率变化对地电阻率观测的影响系数为负， 并与NS测道实际观测的年变形态一致， 但是NS测道位于台站以西并靠近F₅断裂， F₅断裂以西较厚的低阻覆盖层对NS测道地电阻率观测的影响不容忽视．EW、 N45°E测道跨越F₅断裂， 不能用单一的水平层状均匀介质模型来计算表层介质对地电阻率的影响系数， 因此本文根据新沂台电测深资料和测区地质剖面资料建立了三维空间非均匀介质有限元模型， 计算3个测道地电阻率随表层介质电阻率变化时的年变形态， 计算方法参照阎照文(2006)文章．

3.1   稳恒电流场有限元方法

地电阻率定点台站观测采用对称四极装置， 观测时在供电电极A， B输入直流电流， 在测量电极M， N测量电势差. 此问题可视为稳恒电流场计算， 依据麦克斯韦方程组和电荷守恒定律， 在供电电极以外区域电位分布可表示为拉普拉斯方程

稳恒电流场满足狄利克雷-纽曼(Dirichlet-Neumann)边界条件

其中

应用变分原理可得到稳恒电流场拉普拉斯方程的有限元弱解形式

该弱解方程经模型单元离散化、 施加边界条件后就可以对单元节点上的自由度(电位)进行有限元数值求解, 求解出电位分布后可以获得测量电极间的电势差， 进而依据对称四极装置系数计算地电阻率.

3.2   模型建立与检验

新沂台地电阻率测区地势平坦， 地下浅水位维持在3 m左右. F₅断裂两侧基岩面埋深差异达70—80 m， 在建立模型时以F₅断裂为分界面， 东西两侧各自采用水平层状均匀介质模型．新沂台NS向电测深测点距F₅断裂较远， 且台站测区地层构造沿NS向基本一致， 所以依据NS向电测深资料以水平层状均匀介质模型反演其地层的电性参数， 并以此作为有限元模型中F₅断裂以东的初始电性剖面．EW向电测深测点位于F₅以东， 紧邻F₅断裂， 且电测深测线跨越了F₅断裂， 不宜反映F₅断层两侧的差异， 因此在确定模型F₅以西部分电性参数时结合了测区地质剖面．以地下浅水面以上部分作为模型第一层; 同时为了简化模型， 第四层埋深采用与F₅ 以东第四层相同的埋深; F₅西侧基岩埋深在70—80 m左右， 取第二层厚度为70 m; 余下部分为第三层. 结合F₅断层的走向建立初始有限元模型， 分别在EW、 NS向电测深测点以实测电测深曲线时采用的极距对计算模型的电测深曲线. 将计算电测深曲线与实测电测深曲线进行比较， 逐步修改模型， 直至二者大体一致， 从而得到最终的模型参数．模型EW向电性剖面如图3所示. 西盘地下电性层参数与表1中EW方向的反演结果有较大的差异， 主要是由EW向电测深测线跨越了F₅断裂所致. 电测深资料是对F₅断裂两侧不同电性结构的综合反映， 用水平层状均匀介质模型反演的EW向地层电性参数应有差别． 原电测深测点的计算电测深曲线示于图2b, c， 从图中可以看出, 模型在NS和EW向的计算电测深曲线与实测电测深曲线基本一致．

图  3  新沂台有限元模型EW向电性剖面

Figure  3.  Electric section of Xinyi station in EW direction for the finite element model

下载: 全尺寸图片 幻灯片

台站地电阻率观测在水平方向和垂直方向(深度)可视为无穷远边界， 因此有限元模型水平方向的尺寸和模型最底层的有限尺寸将影响计算结果．在建立有限元模型时， 水平尺寸越大, 最底层厚度越厚， 计算结果就越准确. 但模型尺寸越大， 计算量也越大， 因此需要合理地选择模型水平方向的尺寸和最底层厚度．新沂台3个测道的供电极距相同， 且都沿布极中心对称分布， 因而模型在EW方向和NS方向采用相同的尺寸．图4a显示了在图3所示模型F₅两侧最底层厚度为2 km时， NS测道地电阻率计算值随模型水平尺寸D的变化. 在D>5 km之后, NS测道地电阻率计算值已基本上不随D的增加而变化， 所以模型取D=6 km．图4b表示模型水平尺寸D=6 km时， NS测道地电阻率计算值随模型最底层厚度H的变化. 在H>2 km之后, NS测道地电阻率计算值基本上也不随H的增加而变化， 所以模型取H=2.5 km， 最终模型大小为6 km×6 km×2.644 km. 在建立模型时台站观测装置布极中心始终位于模型水平面中心．

图  4  (a) NS向地电阻率计算值随模型水平长度D的变化； (b) NS向地电阻率计算值随模型最底层厚度H的变化

Figure  4.  (a) Calculated georesistivity in NS direction versus D, the horizontal length of the model; (b) Calculated georesistivity in NS direction versus H, the thickness of the bottom layer of the model

下载: 全尺寸图片 幻灯片

3.3   计算结果

新沂台EW测道年变形态相较于其它两测道相对规则， 计算时假定EW测道地电阻率年变形态为正弦型曲线(EW测道基准值+正弦曲线)， 并在模型中计算与该年变曲线相对应的表层介质的真电阻率变化曲线. 以该表层介质真电阻率年变曲线计算了N45°E、 NS测道的地电阻率年变形态， 计算结果示于图5. 图5中所示曲线为3个测道地电阻率年变计算值扣除各自基准值后的年周期成分．计算时一年按月取12个数据点， 以图3所示参数作为计算起始点， 起始点年变幅度为0， 共计算了两年的年变形态． N45°E测道地电阻率年变形态与EW测道相同， 在表层介质电阻率上升和下降时两测道地电阻率也表现为上升和相应的下降变化. N45°E测道年变幅度小于EW测道．NS测道地 电阻率年变形态则与EW、 N45°E测道相反， 年变幅度也小于EW测道．新沂台EW、 N45°E 测道地电阻率年变形态与NS测道相反， 从模型计算结果来看基本符合新沂台3个测道的年变形态．

图  5  有限元模型计算的新沂台三测道地电阻率年变形态

Figure  5.  Annual variations of georesistivity along the three channels at Xinyi station based on finite element model

下载: 全尺寸图片 幻灯片

4.   讨论与结论

依据新沂地电阻率台站电测深资料和测区地质剖面资料， 建立三维空间非均匀层状介质有限元模型， 计算了新沂台EW、 N45°E和NS测道的年变形态. 计算结果表明, EW与 N45°E测道年变形态相同， 随表层介质真电阻率同向变化； NS 测道年变形态与EW、 N45°E 测道相反， 随表层介质真电阻率反向变化．模型的计算结果与新沂台NS测道实测年变形态相反， 与EW、 N45°E测道实测年变形态相符合， 说明用这样一个非均匀层状介质模型可以解释新沂台不同测道地电阻率年变形态相反这一观测现象．

台站所观测的地电阻率变化是地下各种岩石介质电阻率变化的综合反映， 地电阻率影响系数理论认为， 地电阻率相对变化可以表示为每种介质电阻率相对变化的加权和(钱家栋等， 1985； 钱家栋， 曹爱民， 1998)． 地电阻率年变规律完全受测区地下介质电性结构和观测装置控制， 在观测装置固定时， 测区各区域的影响系数则取决于地下电性结构. 将测区地层电性结构简化为水平层状均匀介质模型， 表层介质对地电阻率相对变化的影响系数为正时， 地电阻率年变形态表现为“夏低冬高”; 而当影响系数为负时则表现为“夏高冬低”. 同一台站不同测道年变形态相反， 本质上也是表层介质对各测道地电阻率相对变化的影响系数出现正负．对于确定的电性结构， 测区各区域对地电阻率的影响系数则随观测装置和测点的不同而异， 采用对称四极装置观测时水平层状介质各层的影响系数随极距的不同而变化． 多极距实验观测表明, 在一些电性结构中采用对称四极装置观测时， 长、 短极距所观测的地电阻率会出现年变相反的情况(赵和云， 钱家栋， 1987); 采用多极距装置进行地电阻率观测， 可以依据层状介质影响系数理论反演出测区各层介质真电阻率的变化(钱家栋， 赵和云， 1988)， 从而直观地分析表层介质真电阻率变化与各极距地电阻率年变化的关系．

新沂台NS测道与N45°E测道年变形态相反， 同步性较好， 变化显著的时段在年变中相对较短, 这与苏北地区降雨集中在夏季和秋初有关．从观测资料上可以看出, EW测道与NS、 N45°E测道的同步性相对差一些， NS测道地电阻率年变形态也并非与EW、 N45°E测道完全镜像相反， 这些都是本文所建立的模型无法解释的. 其原因可能是测区地下介质， 甚至是同一层介质不同区域的电阻率在时间上并非同步地发生变化所致. 因此合理地解释这些现象还需要用更多的资料进一步分析．

江苏省新沂地震台全体工作人员参与了异常落实与台站环境考察工作， 彭润禾台长提供了台站电测深资料， 并详细介绍了台站地电阻率测区的地质构造； 评审专家给出了中肯和合理的评审建议. 在此一并表示衷心的感谢．