0 引言

1967年河间地震以来，中国大陆地电阻率观测记录到较多震例，对于地电数据分析基础、孕震过程、异常特征等在诸多文献中详实记载，如：1976年唐山7.8级地震(钱家栋等，1998)、1995年永登5.8级地震(杜学彬等，1997)、2008年汶川8.0级特大地震(张学民等，2009；肖武军等，2009；高曙德等，2010；杜学彬，2010)等，说明地电阻率观测与其他前兆观测相比，具有较显著的映震信号和可信度(冯锐等，2001)。

原始观测得到的地电阻率年变化一般认为是正常变化，这种变化涵盖降雨、气温、气压、潜水位、地下结构应力场等因素。对于介质的各向异性研究，钱复业等(1990)做了改进，地电阻率异常效果明显，后经毛桐恩、阮爱国等进一步研究，认为在地震孕育过程中应力场变化和介质各向异性在一定范围内是存在的(阮爱国，2002)。陆阳泉等(1980)、陈有发(1993)通过实验研究，认为同一测点不同测向地电阻率年变幅存在偏差，说明地电阻率变化具有方向性，地电阻率年变化幅度极小值方向指向当时当地应力场最大主压应力方向。陈有发科研组还利用全国98个地电台作为研究对象，计算1968—1979年地电阻率年变幅极小值，并给出中国大陆应力场3大分布区，结果与阚荣举等(1977)、邓起东等(1979)的研究具有一致性。陈有发(1993)结合1976年唐山7.8级地震前后地电阻率观测资料，分析震中距250 km范围内15个台站地电阻率年变化幅度极小值方向的变化，印证了构造应力场变化与地电阻率异常的必然联系。

随着经济建设的发展，地电阻率观测环境发生较大变化，如基础设施建设、地下水开采量、人为活动干扰、自然环境变化等，均对地电观测产生一定影响，其中降雨量和地下潜水位(或浅水位)变化至关重要，而多数观测台站无同台潜水位(或浅水位)观测，对年变消除产生一定影响，给一些具有物理意义的调查和分析工作带来困难。数学处理方法具有一定科学性和有效性，针对所发生的地震，在不同研究领域有较好应用(杜学彬，2010；刘君等，2013；史红军等，2014)。本文采用斜率合成方法，对中国大陆东部河北省的昌黎、阳原地震台及大陆西部甘肃省的嘉峪关、山丹地震台地电阻率年变幅进行分析，希望多角度认识地电阻率观测数据。

1 台站信息

选取河北省昌黎、阳原地震台及甘肃省嘉峪关、山丹地震台2002—2016年地电观测数据，对比分析相互垂直的2道(NS、EW或NW、NE)电阻率年变幅。各地电台站布极参数见表 1，台站空间分布见图 1。

河北省昌黎地电台位于张家口—北京—蓬莱活动构造带北侧，宁河—昌黎断裂西约5 km处，西侧约20 km处分布滦县—乐亭断裂、约65 km处分布唐山断裂。台站第四系覆盖层厚70 m，基岩风化层厚度8 m，下伏燕山期花岗岩，历史记载含水层厚度20—30 m。采用对称四极装置垂直浅埋方式，布设NS、EW测道，电极埋深2 m，供电极距AB = 1.0 km，测量极距MN_NS = 0.2 km，MN_EW = 0.2 km。年变受季节影响明显，冬高夏低型，平均年变幅约9 — 11 Ω·m。地下电性结构分4层，其中：第1层：h₁ = 6 m，ρ₁ = 45 Ω·m；第2层：h₂ = 24 m，ρ₂ = 90 Ω·m；第3层：h₃ = 85 m，ρ₃ = 55 Ω·m；第4层：h₄ = _∞，ρ₄ = 1 300 Ω·m。第1层至第3层属第四系覆盖层，第4层属基岩层。由于下面为高阻层，所以地电阻率观测受浅层电性层影响较大。

河北省阳原地电台附近主要活动断裂为NEE向延展的南山山前断裂和北山山前断裂，2条断裂均为倾向盆地中心的正断层。年变受季节影响明显，冬高夏低型，平均年变幅约1—3 Ω·m。采用对称四极装置垂直浅埋方式，布设N22°E、N68°W测道，电极埋深2.5 m，供电极距AB = 1.2 km，测量极距MN = 0.2 km。电测深曲线为QQ型。地下电性结构分3层，其中：第1层为第四纪黄土层夹砾石，ρ₁ = 125 — 131 Ω·m，h₁= 18.0— 19.8 Ω·m；第2层为泥层，ρ₂ = 11— 16 Ω·m，h₂ = 50— 100 Ω· m；第3层为黏土层，ρ₃ = 4 Ω·m，h₃= 113— 145 Ω·m；第3层以下为砂质灰岩，ρ₄ = 2— 3 Ω·m。

甘肃山丹地电台测区几乎无第四纪覆盖层，均为上第三系底层，岩性为粘土质沙砾岩。紧靠测区，以西为NNW走向隐伏断层，以南5 — 6 km处为瞭高山逆冲断层，呈NWW—SEE向延伸。年变受季节影响明显，夏高冬低型，平均年变幅约0.3 — 0.6 Ω·m。台址岩性表层10 m为砂砾黄土层，以下为砂岩层、砾层，500 m内红层。采用对称四极装置垂直浅埋方式，布设NS、EW、N45°W测道，电极埋深1.8 m，供电极距AB =1.0 km，测量极距MN_NS = 0.25 km，MN_EW = 0.226 km。

甘肃嘉峪关地电台位于嘉峪关断裂带上，介于阿拉善断块南缘超壳断裂带和祁连山北缘超壳断裂带之间，台址岩性系砾石沙土岩，砂质胶结，质地坚硬，下伏岩石为红粘土，地表为厚30 m的第四系全新统冲积层所覆盖。年变受季节影响明显，夏高冬低型，平均年变幅约0.3 — 0.6 Ω·m。采用对称四极装置垂直浅埋方式，布设N50°E、N45°W测道，电极埋深2.0 m，供电极距AB_NE = 1.120 km，AB_NW = 1.139 km，测量极距MN_NE = 0.37 km，MN_NW = 0.30 km。

2 计算方法 2.1 思路

采用傅里叶滑动消除年周期法、K—L最佳直线拟合消除线性变化及平滑滤波去除粗差方法，从地震前兆数据中提取“真实”变化，以获得有价值的观测值时间序列等(蒋骏等，2000)，也可采用动态距平、去倾去周期等方法提取异常(杜学彬，2010；史红军等，2014)。地电阻率监测受诸多因素影响，原始数据既包含季节变化引起的正常年变，又包含场地环境影响导致的线性趋势变化及其他非震兆成分。本文提出年变幅斜率合成方法，斜率值年动态在时间轴上如一个动点做周而复始的运动，若把该值作为一个动点看待，在“纯净”环境下，NS、EW向变化应具有同步性或一致性特征，尝试对相互垂直的2个测道的斜率值做合成计算，可消除年变，对异常进行放大，削弱某一测道局部异常或其他个别干扰，获得长趋势变化序列，仅包含地壳应变引起的电信号异常和长趋势干扰变化信息，更好地体现介质的物理意义。具体计算思路为：①计算地电阻率年变幅对时间轴的斜率；②计算斜率合成。

2.2 原理

(1) 斜率值。设观测物理量时序序列为y={y₁, y₂, L, y_n}，采用最小二乘法，进行去倾平滑处理，得到{y'}序列，计算该序列的变化范围W_f。利用公式计算曲线变化斜率，即

$ {K_i} = \frac{{\sum\limits_{j = 1}^N {{T_j}} \sum\limits_{j = 1}^N {{y_i} - N} \sum\limits_{j = 1}^N {{T_j}{y_j}} }}{{{{\left({\sum\limits_{j = 1}^N {{T_j}} } \right)}^2} - N\sum\limits_{j = 1}^N {T_j^2} }} $

(1)

自相关系数为

$ {R_i} = \frac{{\sum\limits_{j = 1}^N {{T_j}} {y_i} - \frac{1}{n}\left({\sum\limits_{j = 1}^N {{T_j}} \sum\limits_{j = 1}^N {{y_j}} } \right)}}{{\sqrt {\left[ {\sum\limits_{j = 1}^N {T_j^2} - \frac{1}{n}{{\left({\sum\limits_{j = 1}^N {{T_j}} } \right)}^2}} \right] \cdot \left[ {\sum\limits_{j = 1}^n {y_j^2} - \frac{1}{n}{{\left({\sum\limits_{j = 1}^n {{y_j}} } \right)}^2}} \right]} }} $

(2)

式中，n为滑动步长，N为资料长度，{y}为等间隔前兆数据时间序列，{T}为相应等间隔时间序列。

$ 斜率值序列\;\;\;\;\;\;\;\;\;\;\;{V_i} = {R_i}{K_i}\;\;\;\;\;\;\;\;\;\;\;\;\;\left({i = n + 1, n + 2, L, N} \right) $

(3)

(2) 斜率合成。斜率合成有平均和乘积2种合成方式，某一台站计算方法以体现震例异常明显为选择标准。

$ 平均斜率合成\;\;\;\;\;\;\;\;\;\bar V\left({{t_i}} \right) = \sum\limits_{j = 1}^n {|{V_{ji}}|} /n\;\;\;\;\;\;\;\;\left({i = 1, 2, L, m} \right) $

(4)

$ 乘积斜率合成\;\;\;\;\;\;\;\;\;\prod {V\left({{t_i}} \right)} = \prod\limits_{j = 1}^n {{V_{ji}}} = |{V_{1i}} \times {V_{2i}} \times L \times {V_{ni}}|\;\;\;\;\;\;\;\left({i = 1, 2, L, m} \right) $

(5)

式中，V_ji为斜率，j为各时间序列的编号，i为时间序列内部时元序号，n为时间序列总数，m为时元总数。

3 数据计算与分析

选取2002—2016年昌黎、阳原、嘉峪关、山丹地电台地电阻率观测数据，修正阶变、剔除错误数据，绘制各地电台垂直两测道原始年变幅曲线，见图 2。由图 2可见，昌黎台(属高阻介质台站)年变幅相对偏大，平均地电阻率约10 Ω·m，其余台站年变幅小于3 Ω·m。各地电台变化特征不一致，可能与降雨量、地下潜水位埋深、地应力结构变化等不同有关(《一九七六年唐山地震》编辑组，1982；马钦忠等，2013；刘君等，2013)。

3.1 震例选取

中国东、西部地区地震活动水平不同，为了合理选取地震，西部地区以甘肃地电台周边300 km范围内M ≥ 5.0地震、东部地区以河北地电台周边250 km范围内M ≥ 4.0地震为主，若半年之内同一区域(或地震条带)出现多个地震，则提取震级最大地震。根据以上原则，选取15个地震进行地电阻率年变幅值斜率及垂直2测道合成斜率映震分析。为合理体现地电前兆异常，昌黎、阳原地电台采用乘积合成计算，嘉峪关、山丹地电台采用平均合成计算。

3.2 斜率分析

选取的15个地震地电阻率年变幅值斜率及合成斜率映震曲线见图 3。由图 3可见：斜率及其合成计算对年变有一定抑制作用，且合成斜率优于单一斜率的计算结果；各震例发震前后，合成斜率年变趋势具有一定变化特征，统计结果见表 2。

由图 3(b)可见，9个震例出现在合成斜率1 — 2年下降期内，6个震例出现在合成斜率高值后6个月内，其中：①昌黎台：2006年下降，2011年下降—2012年上升，2014 — 2015年连续下降，分别对应2006年7月4日文安M 5.0地震、2012年5月28日唐山M 4.8地震、2015年9月14日昌黎M 4.2地震；②阳原台：2005年下降—2006年上升，2010年下降，分别对应：2006年7月4日文安M 5.0地震、2010年4月2日大同M 4.5地震、2010年6月5日阳曲M 4.6地震；③嘉峪关台：2005年下降、2009 — 2011年下降— 2012年上升、2015年略下降，分别对应2004年5月4日德令哈M 5.5地震、2009年10月2日肃北M 5.2地震、2009年12月21日德令哈M 5.0地震、2012年5月3日金塔M 5.4地震、2015年11月23日祁连M 5.2地震；④山丹台：2008年下降、2013 — 2016年持续低值，分别对应2008年3月30日肃南M 5.0地震、2013年9月20日门源M 5.3地震、2015年11月23日祁连M 5.2地震、2016年1月21日门源M 6.4地震。

4个地电台地电阻率合成斜率具有1年尺度的下降、2年尺度的下降转上升、3—5年低值背景的变化特征，漏震频次分别为：昌黎台0次、阳原台0次、嘉峪关台2次(2007年和2013年合成斜率下降，无震例对应)、山丹台2次(2006年和2009年合成斜率下降，无震例对应)。因此，采用地电阻率合成斜率分析，甘肃省地电台站周边300 km范围内M ≥ 5.0地震、河边省地电台站周边250 km范围内M ≥ 4.0地震的有效率在70 %以上(概值)；异常变化幅度大小因台而异，可能与地下电性结构有关，如趋势相对变化明显则异常的可能性较大。

4 结论

在地电阻率原始观测数据中，多数中等地震或中强震前后不存在显著异常现象，但并不意味着无电信号。强震前地电阻率异常持续时间可能较长，期间也会出现持续时间较短的短临异常(钱复业等，1990；杜学彬等, 1999, 2000)。在强震前1年内，在大部分地电台站中，与震源机制解最大主压应力方位正交的测道所测电阻率的月速率值出现下降、短临上升变化，与主压应力方位一致的测道所测视电阻率的月速率值则变化相反(李相平等，2016)。中等地震的出现是能量正常释放的表现，孕育的异常电信号可通过合理方法提取(朱宗德等，2001；冯志生等，2004)。由此可见，合理利用数学方法提取短临异常具有一定参考意义。

本文尝试对地电阻率观测中互相垂直的测道做年变幅斜率及斜率合成计算，理论上可以消除年变、放大异常，同时削弱某一测道局部异常或个别干扰，得到长趋势变化序列。对昌黎、阳原、嘉峪关、山丹地电台2003 — 2016年地电阻率观测数据进行计算，并与提取的15个震例做对比，发现计算结果对中、短临异常有一定显示，有效率在70 %以上。地电阻率合成斜率异常主要特征为：1年尺度的下降、2年尺度的下降转上升、3—5年低值背景；异常变化幅度大小因台而异，可能与地下电性结构有关，如趋势相对变化相对明显则异常的可能性较大。

综上所述，地电阻率观测数据合成斜率计算生成新的时间序列值，不仅包含地壳应变引起的电信号异常和长趋势干扰变化信息，而且更接近现实物理意义，可以更好地体现地电观测资料的使用价值，对于其他地震前兆观测物理量，也可通过合成计算得到较好效果。