引言

以地震监测预报为目的的地电阻率试验观测始于1966年河北邢台地震，目前已在全国建成了由80多个台站组成的地电阻率台网，积累了丰富的观测资料和震例，在理论、观测技术、观测环境、异常特征、非震变化、异常机理、承载试验及数据处理等方面均进行了有益的探索(赵家骝，陈才军，1980；钱复业等，1982；陈大元等，1983；钱家栋等，1985；桂燮泰等，1989；陆阳泉等，1990；Ma，2002；Lu et al，2004；杜学彬，2010；Huang，2011；王兰炜等，2011；毛先进等，2014)，并基于台网观测资料对数次中等以上地震作了较准确的中短期(1年尺度)三要素预测(叶青等，2005)，然而在地震预测预报实践中，地震前兆异常的震前研判和确认仍然相当困难.这是由于前兆异常的时空强分布相当复杂，而且异常变化往往是“弱幅度”的，类似于1976年唐山M_S7.8地震中马家沟台(震中距20 km)和2008年汶川M_S8.0地震中成都台(35 km)记录到的那样突出的异常非常少，再说从叠加了年变化成分和背景噪声、观测噪声的时间序列数据上甄别数据异常本身就很困难. “六五”、“七五”期间，地电学科发展了滑动傅氏分析(以下简称“傅氏滑动”)、傅氏拟合分析、相关分析、回归分析、差分频次等地电阻率数据处理方法，引进了气象数据处理中常用的距平方法用于消除月均值数据中的年变化成分.自“九五”以来，^①改进了其中常用的傅氏滑动、距平等方法，发展了归一化变化速率、相对均方差(即变异系数)、差分能量等方法，将其用于处理地电阻率小时测值和各类均值时间序列数据，并研发了地电学方法地震分析预报软件GSEP(郝臻等，2000)，在实际应用中均显示出其效用.但是，随着这些被改进的方法和新方法的广泛使用，相关研究人员迫切需要详细了解其原理、使用条件、异常指标以及前期的数据准备、预处理要求等.鉴于此，本文针对上述需求对其中常用的8种方法的原理和使用条件等进行系统阐述，以期这些方法能更有效地用于地电阻率数据的处理、分析及异常甄别中.

① 杜学彬等杜学彬，卢军，王志贤，薛顺章，郝臻，钱复业，赵和云，赵玉林，张学民. 1999.地电中短期前兆识别、标志体系及预报方法研究.项目结题报告.

1.   消除年变化方法

经多年观测的震例总结，一般把地电阻率日、5日、旬均值和月均值等原始曲线或消除年变化曲线上的相对变化幅度不小于1%的持续变化判定为异常.绝大多数台站的地电阻率均存在年变化，年变化幅度为2%—3%，有时甚至高达20%—30%不等(钱家栋等，1985).年变化的幅度、形态、相位变化本身是一种异常，但是，年变化又为研判幅度为1%的持续变化异常带来了困难，异常往往被“淹没”在年变化中，所以需要消除年变化.目前，常用的消除年变化的方法主要有傅氏滑动、距平和动态距平等3种，下面将逐一介绍.

1.1   傅氏滑动方法

傅氏滑动方法是赵跃辰和刘小伟(1984)提出的消除年变化的数据处理方法，当时主要用于消除地电阻率月均值数据中的年变化，“九五”期间该方法发展成可处理小时值和各类均值数据^①.该方法原理为：设地电阻率观测数据的时间序列为{y}={y₁, y₂, …, y_N}，即{y_n}, n=1, 2, …, N，N为有限长度，年变化周期为T，即一年中有等时间间隔的T个年变化成分测值，序列{x}={x₁, x₂, …, x_N}为{y}中的年变化成分.使用三角级数公式由序列{y}中的T个数据组成的长度为T的子序列{y_n-T+1, y_n-T+2, y_n-T+3, …, y_n}拟合出年变化成分(即基波)的第n个值x_n，再由{y_n-T+2, y_n-T+3, …, y_n+1}拟合出x_n+1，依此类推.计算年变化的公式为

① 杜学彬，卢军，王志贤，薛顺章，郝臻，钱复业，赵和云，赵玉林，张学民. 1999.地电中短期前兆识别、标志体系及预报方法研究.项目结题报告.

其中

令式(1) 中j=n, 有

令

由式(3) 则可得到消除年变化的新序列{G}.对序列{y}中前T-1个数据{y₁, y₂, …, y_T-1}通过三角展开补充基波成分{x₁, x₂, …, x_T-1}，然后由式(3) 得{G₀，G₁，…, G_T-1}.分别令式(1)—(3) 中的T=12，73，365和8760，则可消除月、5日、日均值和小时测值数据中的年变化成分.需要指出的是，年变化周期按1年365天计，因年变化是连续渐变的，366天的年份来说误差可忽略.

一般而言，大部分台站连续数年的年变化形态连续、稳定，但也有一些台站的年变化形态、幅度逐年变化或某些年份年变化起伏时间时移.傅氏滑动方法的特点是三角系数a_n随年变化强弱而动态变化，有自适应消除年变化成分的效果, 适合处理上述年变化起伏时移的数据. 图 1a为四川成都台在汶川M_S8.0地震前的地电阻率日均值原始曲线，从该曲线仅能判定自2004年开始地电阻率长期下降、2007年年变化消失以及震前约两个月出现上升(恢复)变化，而难以确定震前下降异常究竟从何时起算；但从图 1b经傅氏滑动消除年变化的曲线上则可明确甄别出显著下降异常始于2006年8月，最大相对变化Δρ_s/ρ_s约为-7%，且自2006年8月、2007年4月和12月开始出现了3次地电阻率的加速下降.

图  1  2004年1月—2008年6月成都台N58°E测道地电阻率日均值曲线

(a)原始曲线；(b)傅氏滑动曲线

Figure  1.  Daily mean value curves of georesistivity measured by channel N58°E of station Chengdu from January of 2004 to June of 2008

(a) Raw curve; (b) Moving Fourier's curve

下载: 全尺寸图片 幻灯片

上述相对变化幅度Δρ_s/ρ_s估算如下：首先，选定正常变化时段地电阻率平均值ρ_s1，若曲线呈现如图 1b所示的下降趋势，则需估算正常变化趋势的平均值ρ_s1；然后，估算“异常”时段的ρ_s2，再通过

估算出异常的相对变化幅度.

一般来说，经傅氏滑动方法消除了年变化的曲线大致保持了与原始曲线除年变化之外的吻合，并可显示出近距离大震、强震前后地电阻率异常变化的持续时间、幅度和形态特点.由图 1可见：在汶川地震前成都台(距离主震的震中距35 km)的地电阻率从2006年8月下降，经3次加速下降后于2008年2月降至4.5年以来的最小值，最大相对变化Δρ_s/ρ_s约为-7%；之后开始上升(恢复)，直至发生汶川M_S8.0地震，呈现出“震前中短期阶段ρ_s下降—短临阶段ρ_s上升—发震”的过程，该过程与膨胀-扩散模式(图 2)所示的震源区电阻率变化的3个阶段(第Ⅱ阶段电阻率减小，第Ⅲ阶段加速减小，第Ⅵ阶段陡然增大)(Nur，1972；Scholz et al，1973；Mjachkin et al，1975)十分吻合.

图  2  膨胀-扩散模式所示的震源区岩石电阻率变化(引自Mjachkin et al, 1975)

Figure  2.  The resistivity change within the focal area predicted by the dilatancy-diffusion model (after Mjachkin et al, 1975)

下载: 全尺寸图片 幻灯片

1.2   距平方法

傅氏滑动方法简便，对多数的长程观测数据去年变化很有效，但也注意到，在曲线上与背景噪声接近时，要消除年变化，该方法的效果欠佳，特别是大幅度的突跳变化易造成虚假的振荡变化.为此，引进了气象领域广泛应用的距平方法.该方法前期主要用于消除地电阻率月均值数据中的年变化成分，“九五”期间发展成为消除地电阻率小时测值和日、5日、旬均值数据中的年变化成分^①.

① 杜学彬，卢军，王志贤，薛顺章，郝臻，钱复业，赵和云，赵玉林，张学民. 1999.地电中短期前兆识别、标志体系及预报方法研究.项目结题报告.

距平方法的原理为：设地电阻率观测数据的时间序列为{y}={y₁, y₂, …, y_N}，长度为N，其中的年变化成分为{x}={x₁, x₂, …, x_T}，年变化周期为T.若年变化形态、幅度在多年(如k年)观测中(基本)不变，即x_i=x_i+kT(不考虑其它因素引起的地电阻率变化)，则y_i-x_i=y_i+kT-x_i，即可求出所谓的“标准年变化”.对时间序列{y}去倾和中心化(中值为M)，得到{y′}.标准年变化的模型为

式中：k表示年份排序，k=0, 1，…，L-1；L为序列{y′}中y′_i, y′_i+T, …, y′_i+kT出现的频次，若序列前、后时段有非整年的数据，缺数时段不参加平均.序列{x}={x₁, x₂, …, x_T)则是所谓的“标准年变化”，然后由

得到计算消除年变化的序列{F}.同傅氏滑动方法，分别令式(5)—(6) 中T=12，73，365和8760则可得消除月、5日、日均值和小时测值曲线的年变化成分.

图 3a为云南元谋台地电阻率月均值原始曲线，在2009年楚雄M_S6.1地震(84 km)前后出现了年变化消失的中短期异常和两次短时下降异常.由图 3b可见，在经过距平方法消除年变化的曲线上中短期异常更突出，且可以判定在汶川M_S8.0地震(610 km)和附近的会理M_S6.1地震(67 km)前未记录到任何异常.需要说明的是，图 3中汶川地震、会理地震未记录到异常并不能否定楚雄地震前的异常，地电阻率异常与活动构造展布、震源机制有关，杜学彬(2010)就此问题已作过讨论.

图  3  2006年9月—2012年12月元谋台EW测道地电阻率月均值变化曲线

(a)原始曲线；(b)距平曲线；(c)动态距平曲线

Figure  3.  Monthly mean value curves of georesistivity measured by channel EW of station Yuanmou from September of 2006 to December of 2012

(a) Raw curve; (b) Monthly departure curve; (c) Dynamic monthly departure curve

下载: 全尺寸图片 幻灯片

1.3   动态距平方法

距平方法对消除年变化相位、年变化幅度差异不大的年变化成分是可行的，且简便、直观，不会因个别大幅度的短时突跳等变化造成虚假异常，但是，如果年变化存在逐年相位差、幅差或者在某年份无年变化的情况下，该方法消除年变化会产生畸变.例如，图 3a中2008年下半年至2009年上半年无年变化时段的数据也参与了“标准年变化”的估算，2009年楚雄地震前的原始曲线(图 3a)上本无上升变化，但距平去年变化曲线(图 3b)上对应的时段却出现了上升变化，这样看来显然是产生了畸变.为此“九五”期间又发展了动态距平方法，即先检验每年是否存在年变化成分，仅用有显著年变化成分的数据估算“标准年变化”，然后再去除年变化成分^①.

① 杜学彬，卢军，王志贤，薛顺章，郝臻，钱复业，赵和云，赵玉林，张学民. 1999.地电中短期前兆识别、标志体系及预报方法研究.项目结题报告.

动态距平方法的原理为：设地电阻率时间序列为y={y₁, y₂, …, y_N}，长度为N，某年的年变化成分为序列x={x₁, x₂, …, x_T}，先对{y}序列去倾和中心化得到{y′}序列.某年子序列中的年变化成分为

其中

令，以为判断年变化成分显著性的条件^①，对{y′}序列中满足该条件的年份的数据由式(5)—(6) 估算出“标准年变化”{x}序列(L为满足条件的年份数)，然后扣除其{y′}序列中的标准年变化成分，再得到年变化校正序列{F}；对不满足条件的年份数据则不予扣除.式(7) 中，对测量数据时间序列前、后端不满足1年时段的数据，向后或向前延续选取1年的数据，作为缺数年份的数据.

由元谋台地电阻率月均值距平曲线(图 3b)与动态距平曲线(图 3c)可看出，动态距平方法消除年变化较距平方法更为客观.例如，在2009年初距平曲线上出现了上升“畸变”，而动态距平曲线保留了原始曲线的变化形态.同理，该方法可消除日、5日、旬均值和小时值数据的年变化成分.

2.   无量纲方法

2.1   归一化变化速率方法

通常由各类地电阻率原始曲线及其消除年变化的曲线来研判异常，但震后能够识别出如1976年河北唐山M_S7.8地震前天津宝坻台(80 km)和河北昌黎台(70 km)、马家沟台(10 km)以及1976年四川松平M_S7.2地震前甘肃武都台(105 km)和2008年四川汶川M_S8.0地震前成都台(35 km)出现的显著中期异常(国家地震局《一九七六年唐山地震》编辑组，1982；钱家栋等，1985；桂燮泰等，1989；杜学彬，2010)的情况是极少的.甚至，在M_S7.0—8.0地震震中区的台站所观测到的地电阻率曲线上也不是均能在震后识别出异常，当然对于M_S5.0—6.0地震，异常识别则更难.因为在绝大多数情况下，原始曲线或消除年变化曲线上的异常幅度很小，其变化形态并不像图 2中膨胀-扩散模式中的Ⅱ—Ⅳ阶段、汶川地震前成都台中期下降异常那样可识别.例如，唐山地震发生后，钱复业等(1982)和王新华等(1984)总结出部分地电台记录到了地电阻率中长期、中短期和短临异常，但认为距离震中仅114 km的平谷台未记录到地电阻率异常.为此，杜学彬等(2001)在“九五”期间发展了归一化变化速率方法(normalized variation rate method，简写为NVRM)，该方法初期称为月速率方法.

归一化变化速率方法的原理为：先剔除地电阻率时间序列{ρ_s}中个别大幅度突跳的错误数据、去倾和去年变化，然后对ζ个数据组成的子序列线性回归，计算其斜率，且保持步长ζ逐一向后滑动计算初始速率时间序列. 序列为

对进行归一化处理，得到时间序列为

式中：K_i为{ρ_s}中第i个子序列的斜率；R_i为相关系数，用来压制子序列中奇异数据对K_i的影响；m为初始速率的平均值；σ_n-1为的均方根误差；δ为正小量.在前ζ-1个数据补“0”即可得到完整的归一化变化速率时间序列.

实际上，从时间序列数据曲线上甄别异常的变化形态、幅度等，其关键可归结为分析曲线变化速率的大小及其正(上升)、负(下降)变化的形态；归一化变化速率方法紧扣这个关键将人工分析的过程用“变化速率”表征，并统一异常指标(阈值).只要被分析的数据时段固定，即可避免人为研判正常与异常的不确定性，且实际应用证明该方法用于月均值数据分析1年尺度中短期异常的分辨力较高(杜学彬等，2015).

杜学彬等(2001, 2015)曾介绍了归一化变化速率方法的原理及其使用条件，但后来注意到相关研究人员在应用该方法时并未遵循其使用条件，还有误报该方法的原理以及质疑归一化变化速率异常真实性的情形出现.关于归一化变化速率异常的质疑涉及两方面的问题：其一，归一化变化速率曲线上出现类似图 4a中昆仑山口西M_S8.1地震前远距离(距甘肃山丹台979 km)的异常，甚至于在2004年苏门答腊M_S8.9地震和东日本M_S9.0地震这样远的大震前中国大陆相当多的台站(最大震中距约4000 km)出现异常；其二，假设观测数据是多年持续的一条水平直线(即均方差σ_n-1=0)，在此水平直线上加上大幅度的突跳，这样归一化变化速率曲线就出现了虚假异常，且突跳变化的持续时间、幅度在数据中所占比例越大，虚假异常越显著.关于前一个问题，作者认为此类远距离异常不是地震“源”过程引起的，而是属于地球自转年尺度加速所造成的地震与异常所表现出来的“同源异象”现象(杜学彬，谭大成，2000；杜学彬，2010).至于后一个问题，上述假设数据在正常的地电阻率观测中并不存在.我国绝大多数台站的供电极距AB=1000—1200 m，测量极距MN=AB/5-AB/3，在温度0℃—40℃范围内和供电电位差大于10 mV的条件下，测量仪器的电阻率最大容许误差为±(0.1%读数+0.02 Ω·m)，整个台网的地电阻率测值是处于5—750 Ω·m范围内，再加上测量线路、电极极化等因素以及测量数据产出为小数点后3位，实际测值不可能出现长期不变的一条水平直线，除非台站的观测系统(测量仪器、电极、室内外线路等)出现故障而造成错误数据.从台站观测系统和观测环境(自然和人为电磁噪声、地下介质应力扰动等)的影响来说，测量数据总有误差，式(9) 中的σ_n-1不可能为零，由正小量δ约束也不会得到的归一化变化速率异常.另一方面，在《地震及前兆数字观测技术规范》(中国地震局, 2001)中规定“连续4小时以内地电阻率小时值与正常值的相对变化|Δρ_s/ρ_s|＞5%的粗差数据应剔除”，同时在测量过程中剔除相对均方差大于3%的小时测值.据地电台网运行管理的要求，小时测值(日次数大于5) 超过前一天正常日均值1%的突跳和日均值出现大于上月均值1%并持续3天以上要及时检查和说明原因，已明确人为、观测系统等所造成的变化应在预处理时剔除.这样，通过上述的数据产出和质量监控过程可得到正常观测时的日均值，特别是月均值不可能出现瞬间的大幅度突跳，因此上述假设的测量数据属于错误数据.事实上，很难相信地下介质结构变化能引起地下介质特别是近地表浅层介质的视电阻率会出现瞬间的大幅度突跳.一般来说，正常观测数据经去倾、去年变化后，数据时间序列会稳定在基准值附近，再经ζ长度子序列回归、归一化处理，序列为似稳甚至平稳变化. 图 5给出了1986—1989年甘肃山丹台、定西台和天津宝坻台地电阻率月均值元素的分布检验，当时这些台站无环境干扰.可见，优势集中在-2.4—2.4之间(含附近发生M4.0以上地震时的异常数据)，似呈正态分布，说明正常情况下，序列中出现异常的几率很低.

图  4  1997—2004年山丹台EW测道地电阻率月均值曲线

(a)归一化变化速率曲线；(b)原始曲线；(c)傅氏滑动曲线；(d)距平曲线；(e)动态距平曲线

Figure  4.  Monthly mean value curves of georesistivity measured by channel EW of station Shandan in 1997-2004

(a) Normalized variation rate curve; (b) Raw curve; (c) Moving Fourier's curve; (d) Monthly departure curve; (e) Dynamic monthly departure curve

下载: 全尺寸图片 幻灯片

图  5  山丹台(a)、定西台(b)和宝坻台(c)的归一化变化速率值(NVRM)分布检验

图中N为各台个数，n为按速率值0.5间隔统计的| |个数

Figure  5.  Statistical distribution of normalized variation rate values for stations Shandan (a), Dingxi (b) and Baodi (c)

N is the number of , and n is the number of || in the value interval of 0.5

下载: 全尺寸图片 幻灯片

通过上述讨论可知，归一化变化速率方法的原理和异常物理机制均很清晰.该方法得到了较广泛的应用，并取得了明显的处理效果(张学民，武有文，2000；郭宝娥等，2002；张学民，翟彦忠，2002；王福才等，2004；胡哲等，2011；史红军等，2014；周剑青等，2014).由图 4b-e所示的甘肃山丹台地电阻率月均值的原始曲线及其傅氏滑动、距平、动态距平曲线，很难甄别出2003年甘肃民乐—山丹M_S6.1地震(43 km)前相对变化幅度不小于1%、变化形态清晰的中短期异常，但在归一化变化速率曲线(图 4a)上清晰地展示了震前1年尺度下降、临震前夕显著上升的异常，且先期下降—短临阶段上升—发震的变化过程与膨胀-扩散模式中震源区电阻率变化过程(图 2)类似.此外，唐山M_S7.8地震中平谷台(111 km)和1991年大同M_S5.8地震中阳原台(33 km)的归一化变化速率曲线也呈现出类似的变化过程(杜学彬等，2006；杜学彬，2010).不过，归一化变化速率方法也有其不足，杜学彬等(2015)已作详细讨论，不再赘述.另外，归一化变化速率方法也可用于处理其它观测物理量时间序列的异常.

2.2   年变化识别方法

在地电阻率时间序列曲线上识别年变化幅度、形态、相位变化也是探索地震前兆异常的常用方法之一.某些台站年变化幅度小，且有时迭加有复杂的短时变化，使得人工分析年变化较为困难，所以在“九五”期间发展了年变化识别方法^①.该方法原理简单，先对原始曲线去倾，然后由式(8)—(9) 计算变化速率即可.以汶川M_S8.0大震震中区附近的成都台为例，汶川大震后于2008年6月—2010年5月在该台布极区内安置灾民，造成了NE测道的严重干扰，但对NW测道影响小(杜学彬等，2015). 图 6给出了该台NW测道地电阻率月均值原始数据曲线和年变化形态曲线，可见：自2004年后难以从原始曲线(图 6a)上识别年变化，但从年变化形态曲线(图 6b)看，每年年底至次年年初地电阻率减小的年变化较原始曲线要清晰得多；从全程年变化曲线来看，年变化幅度增大前、后时段在台站附近及其周围发生了M_S≥5.0地震，且汶川地震前后(2007—2009年)出现了一年两次起伏的年变化形态“畸变”现象，年变化幅度减小，似乎是年变化“畸变”异常.该方法用于定性地分析年变化异常，可作为异常识别的参考方法.

① 杜学彬，卢军，王志贤，薛顺章，郝臻，钱复业，赵和云，赵玉林，张学民. 1999.地电中短期前兆识别、标志体系及预报方法研究.项目结题报告.

图  6  1998年5月—2010年6月成都台地电阻率月均值原始数据曲线(a)和年变化形态曲线(b)

Figure  6.  The monthly mean value curve (a) and annual variation curve (b) of georesistivity measured at the station Chengdu from May of 1998 to June of 2010

下载: 全尺寸图片 幻灯片

3.   相对均方差方法

在地电阻率观测中采用正、反向供电方式测量供电电位差ΔV.首先由电流电极A和B正向导入电流I测量ΔV₊(ΔV₊=V₊+V_sp+，其中V₊和V_sp+为正向供电电位差和自然电位差)，然后反向导入电流-I测量ΔV_-(ΔV_-=V_-+V_sp-，其中V_-和V_sp-是反向供电电位差和自然电位差)，再得到单次供电的电位差ΔV(ΔV=|ΔV₊-ΔV_-|/2).一般每单次供电时长约5 s，每小时的测值为5或10次单次供电ΔV的均值，这样多数台站测得的地电阻率小时值的相对均方差日均值可达0.1%，月均值可达10^-4.由于地下介质充放电(傅良魁，李金铭，1980)的影响，从开始供电到电场稳定后测量有短时的延时，再加上每小时5或10次单次供电也有时间先后，当V_sp叠加有电脉冲扰动等快变化时会使测量误差增大.例如，地电暴发生时一些台站出现地电阻率测值增大或减小、均方差增大的现象即源于此.地震发生和岩石破裂前已记录到了地震电信号、电脉冲扰动等快变化现象(郭自强等，1988; Varotsos et al，1993；钱书清等，2003；郝锦绮等，2004；马钦忠等，2013).大震强震前夕，震中附近及周围的地电日变化“波形畸变”的研究表明，引起畸变的快变化的功率谱密度呈数量级变化(李宁等，2007；范莹莹等，2010；刘君等，2011).此类地电场快变化也会造成地电阻率增大或减小、均方差增大的现象.因此，地电阻率相对均方差也是探索地震短临前兆的测项.

相对均方差法的原理为：地电阻率小时测值或日均值的相对均方差序列K_σn-1为

计算小时值相对均方差K_σn-1时，σ_n-1为用该小时5或10次地电阻率单次测量的均方差，ρ_s为对应的地电阻率小时值；计算日均值相对均方差K_σn-1时，均方差σ_n-1为该日每小时均方差的日均值，ρ_s为该日地电阻率均值.

图 7a给出了汶川M_S8.0地震前后2008年1—7月甘肃陇南台(297 km)N73°W测道地电阻率小时值的K_σn-1曲线，可见：① 主震前和当年7月下旬台站附近的强余震前约10—20天内，K_σn-1显著增大，属震前短临阶段的变化；② 大约从2008年3月开始K_σn-1整体上比以往增大，似乎远距离的新疆于田M_S7.3地震(2147 km)前K_σn-1也呈增大趋势.上述K_σn-1增大的时段(特别是主震前)与震中区及周围的成都台地电场、剑阁台电磁扰动和青藏高原东北缘中法国际合作的松山台、代乾台地电场功率谱密度(power spectrum density, 简写为PSD)增大以及法国地震电磁卫星Demeter磁场PSD增大的时段相吻合(范莹莹等，2010；刘君等，2011). 图 7b给出了同期四川江油台(距主震震中151 km)NE测道地电阻率日均值K_σn-1，可以看出大约从主震前两天开始K_σn-1明显增大，显示出突出的临震阶段异常变化.当时该台使用DDC-2A电子自动补偿仪测量地电阻率，每天定时数次人工观测，因在主震前未校检仪器，数据可靠性存在质疑，但在震前较长时段内测量仪器运行状态未改变，上述短临阶段的异常变化还是可信的.此外，还可见汶川主震后的强余震期间K_σn-1整体上较主震前增大，且表现出与主震前类似的K_σn-1逐步减小的变化趋势.从这个意义来说，震前K_σn-1增大的异常变化仍有可信性.然而，在该台站地电阻率日均值、月均值曲线上并无显著的短临变化.另一方面，据震后了解，观测人员在主震前测量时DDC-2A仪器指针快速地大幅摆动，读数困难.以往台站近距离发生地震之前，也曾注意到DDC-2系列仪器指针剧烈摆动、数字仪器均方差增大的现象.这种现象也反映了震前电磁扰动变化剧烈，与K_σn-1增大现象的物理机制相同.

图  7  2008年1—7月地电阻率相对均方差K_σn-1曲线

(a)陇南台N73°W测道小时值曲线(ZD8BI地电仪测量)；(b)江油台NE测道日均值曲线(DDC-2A电子自动补偿仪测量)

Figure  7.  The curves of relative RMSE of georesistivity from January to July, 2008

(a) Hourly-value curve measured by channel N73°W of station Longnan (ZD8BI instrument); (b) Daily mean value curve measured by channel NE of station Jiangyou (DDC-2A instrument)

下载: 全尺寸图片 幻灯片

4.   差分方法

4.1   差分频次方法

差分频次方法能够抑制长周期变化成分而突出扰动变化，常用于地震前兆异常分析，习惯上用一阶、二阶差分频次.方法原理如下：设长度为N的日均值或小时值时间序列{y}={y₁, y₂, …, y_N}，由{y}可得一阶差分序列{F}={F₁, F₂, …, F_N-1}或二阶差分序列{G}={G₁, G₂, …, G_N-2}，进而得到{F}或{G}的均方差σ_n-1，按月子序列(日均值差分)或按日子序列(小时值差分)统计{F}或{G}中元素绝对值|F_i|(或|G_i|)≥2σ_n-1的频次.该方法用于地电阻率及其(相对)均方差、地电场观测数据的小时测值、日均值时间序列, 能突出地震短临阶段的异常变化.例如，汶川M_S8.0地震前甘肃天水台(411 km)EW测道地电阻率小时值曲线出现突跳但不集中(图 8a)，而在一阶差分频次图(图 8b)上主震前约两个月至宁强M_S6.4余震(242 km)以及8月初青川M_S6.1余震前均出现密集的一阶差分高频次.

图  8  2008年1—7月天水台EW测道地电阻率小时测值(a)、一阶差分频次N(b)和一阶差分能量E (c)

Figure  8.  Hourly value (a), the first-order difference frequency N (b) and the first-order difference energy E (c) of georesistivity measured by channel EW of station Tianshui from January to July, 2008

下载: 全尺寸图片 幻灯片

4.2   差分能量方法

在某时段内地电阻率差分频次高，可能反映了在应力作用下地下介质中激发的低频电磁场扰动，在震前多属短临信息.按月或日统计一、二阶差分频次仅反映了差分值不小于阈值(例如{F}或{G}的2σ_n-1)的频次，而未考虑差分值的大小.例如：差分值1.0 Ω·m与10.0 Ω·m所反映的介质活动强烈程度完全不同；不同台站和测道的构造、介质、电性条件等不同，其对应力扰动的响应也有差别.因此，在不同台站或测道之间的差分频次图上很难确定统一的异常判断标准，只能根据该台站或测道以往对应地震的情况进行定性地对比分析.为此，“九五”期间杜学彬等^①发展了差分能量方法，其原理为：按月或日(也可用其它的统计时间间隔)统计序列{F}或{G}的各子序列中元素|F_i|(或|G_i|)≥2σ_n-1的能量E_k，得到一、二阶差分能量时间序列{E}，这样便兼顾了差分频次和差分值大小两个因素.统计公式为

① 杜学彬，卢军，王志贤，薛顺章，郝臻，钱复业，赵和云，赵玉林，张学民. 1999.地电中短期前兆识别、标志体系及预报方法研究.项目结题报告.

式中，|x_i|为子序列中|F_i|(或|G_i|)≥2σ_n-1的元素，k为子序列的序数，m为第k个子序列中|x_i|≥2σ_n-1的元素个数，i为第k个子序列中|x_i|≥2σ_n-1元素序数.由式(11) 可得一、二阶差分能量时间序列{E}.对不满1个月或1天的数据时段，按实际天数或小时值的个数统计.为了使序列{E}中异常指标具有可操作性，又定义了对数能量式

{E}序列中的元素满足式(12) 为异常，不满足的则赋予0.

比较图 8b与图 8c可见，在汶川M_S8.0主震和宁强M_S6.4强余震前天水台差分能量图中出现的异常少于差分频次图，后者异常较明显. 图 9为甘肃武都台(297 km，距宁强M_S6.4强余震震中区113 km)N73°W测道地电阻率小时值相对均方差的一阶差分频次(图 9a)和一阶差分能量图(图 9b)，可见主震前差分频次、差分能量均显示了突出的高频次以及能量增大的短临异常.武都台距7月24日宁强M_S5.7强余震震中区仅77 km，图中差分频次、差分能量异常比主震前还显著.

图  9  2008年1—7月武都台N73°W测道地电阻率小时值相对均方差的一阶差分频次N(a)和一阶差分能量E(b)

Figure  9.  The first-order difference frequency N (a) and the first-order differential energy E (b) of the relative RSME of hourly value of georesistivity measured by channel N73°W of station Longnan from January to July, 2008

下载: 全尺寸图片 幻灯片

另外，上述一、二阶差分频次和差分能量统计方法还可以用于地电场分钟值、小时均值和日均值曲线上的异常分析(刘君等，2011).

5.   讨论与结论

5.1   异常识别指标

本文分别介绍了傅氏滑动、距平、动态距平、归一化变化速率、年变化形态识别、相对均方差、差分频次和差分能量共8种方法的原理、数据处理过程、异常分辨能力、异常指标和异常的物理机制及其在数据处理中仍存在的不足.

傅氏滑动、距平和动态距平方法以及从原始曲线上识别异常的方法通常用于分析中期和短临异常.经多年的震例积累和分析，这些方法的异常指标为地电阻率相对变化幅度不小于1%(即Δρ_s/ρ_s≥1%)的连续上升或下降的持续性变化，包括异常幅度、持续时间和变化形态等信息，适用于在各类均值和小时值曲线上研判异常.近年来，随着震例积累和观测技术的发展，对台站观测环境加大了保护力度，有的台站出现0.6%左右的相对变化，已与正常背景变化明显不同，这种情况下即使变化幅度小也可识别为异常.这3种方法的异常识别标准是基本明确的，不足之处是人为研判正常背景变化和异常变化往往有不确定性，特别是不同研究人员研判中等地震或远距离地震对应的“弱变化”异常往往存在不确定的结果.其次，在不同均值曲线和是否消除了年变化的曲线上识别异常也有差异.例如，成都台N58°E测道在距离仅35 km的2008年汶川M_S8.0地震前记录到了十分突出的地电阻率中期下降异常，但对该异常的起始时间、变化幅度的认定仍有差异. 钱家栋等(2013)基于月均值原始曲线分析异常认为，从2006年2月开始趋势下降，最大幅度为-6.8%；而杜学彬等(2015)用傅氏滑动方法消除年变化的月均值曲线分析异常认为，2006年8月开始趋势下降，幅度为-5.9%，再用消除年变化的日均值曲线分析认为，异常起始时间同上，幅度为-7%(图 2).确定震前异常幅度、持续时间和形态是研判震级和发震时间的重要依据，所以采用原始曲线或消除年变化曲线以及用日均值、月均值等哪种均值曲线研判异常也很关键，应由多种方法综合分析，确认可信的异常参数.

归一化变化速率方法的异常指标为变化速率值不小于|±2.4|，对所有台站和测道的数据指标统一.用该方法识别的异常仍包含异常变化的幅度、形态和持续时间信息，只要选定数据时段，并经剔除错误数据等预处理，异常判定则无需人为干预，避免了在原始曲线、消除年变化曲线上识别异常的不确定性，特别对识别震前“弱幅度”异常有明显优势(张学民，武有文，2000；张学民，翟彦忠，2002).一般地，该方法处理地电阻率月均值数据识别1年尺度中短期异常的效果更佳.其次，用该方法判定出变化速率异常后再对比分析消除年变化曲线或原始曲线上的变化能较准确地判定异常.但是，该方法仍需依据地电阻率长趋势变化转折的时间选取拟处理数据的时段，且存在选定数据滑动步长ζ等不足(杜学彬等，2015)，该算法尚需进一步改进.

年变化形态识别方法适用于分析地电阻率中短期异常，相对均方差、差分频次方法和差分能量方法适用于研判短临异常.据目前的震例积累和分析，尚不能得到定量甚至于半定量、半定性的异常指标.目前是通过学习“正常”时段的变化来定性分析是否出现了异常，所以使用上述4种方法处理的异常仅有参考意义，异常指标尚需震例分析的积累.

5.2   数据异常与前兆异常

文中8种方法处理的地电阻率异常属于地电阻率时间序列出现的异常(简称为“数据异常”)，严格地讲，数据异常与地震“前兆异常”并不相同.数据异常仅指通过原始曲线或经数据处理得到的时间序列上满足异常指标的变化，某台站出现了异常并不预示在该台站可监测范围内一定发生显著地震.地震孕育、发生是地下介质的运动，地电阻率异常也是地下介质运动过程中微裂隙、导电流体活动等造成的介质(真)电阻率变化引起的视电阻率变化.视电阻率相对变化与介质微裂隙、骨架电阻率和裂隙水电阻率之间的本构关系(杜学彬等，2007)以及扩展的阿尔奇(Archies)公式(钱家栋等，2013)较清晰地解释了地电阻率异常变化的物理机制.由此，非震活动引起的台站所处地下介质物理环境改变也能造成地电阻率异常.另一方面，长期地电阻率观测证明，在数千千米远的远距离M_S8.0及以上大震前(例如，2004年苏门答腊M_S8.9、2011年东日本M_S9.0大震等)，在中国大陆的多数地电阻率台站均出现显著的中短期异常，此类异常与具体的震源孕育发生的过程并无直接联系，异常与大震发生属于“同源异象”现象(杜学彬，2010)，也属于数据异常.数据异常也不是台站观测系统、测量过程、场地环境变化等引起的干扰变化，而是地下介质运动自然产生的异常.而前兆异常属于数据异常中与震源孕育、地震发生有直接关系的异常，其时、空、强分布和变化形态符合经震例剖析而总结出的异常特征，主要集中在震中区及周围.据大量的震例统计，M_S6.0—7.8地震前中短期异常绝大多数分布在约400 km范围内(杜学彬，2010)，若台站与震源区有构造关联性，异常空间分布范围会更大.其中，大震、强震近震中区负异常集中的空间范围与震源断层破裂尺度对等、与震源特性直接关联的各向异性变化以及与震源区-台站的活动构造关联的异常(杜学彬，2010)属于典型的前兆异常.所以，不能将数据异常与前兆异常混为一谈，不要期望在地电阻率时间序列上出现异常就一定会在台站附近发生显著地震.

5.3   结论

在长期观测中，极少出现类似2008年汶川M_S8.0地震前成都台(35 km)、1976年唐山M_S7.8地震前昌黎台(72 km)等显著异常，甚至在M_S7.0—8.0地震的震中区也如此.相比正常背景变化和观测系统、场地环境等引起的观测噪声来说，绝大多数异常为“弱幅度”变化，需用有效的方法处理观测数据才能分析、识别出异常，若台站距离远、地震震级小更需如此.

文中介绍的8种数据处理方法的原理简明，异常物理机制较清晰，其中：归一化变化速率方法有量化、统一的异常识别指标；傅氏滑动、距平、动态距平方法有半定量、半定性的指标；相对均方差、差分频次、差分能量和年变化形态识别方法有定性的指标，用这些方法判定的异常只有参考意义.

归一化变化速率、傅氏滑动、距平、动态距平和年变化形态方法通常用于识别震前中期至短临阶段的异常，相对均方差、差分频次和差分能量方法则主要用于识别短临异常.

受台站构造介质环境和观测噪声等影响，地电阻率异常往往表现为多样性，原始观测数据是最真实的，通过数据处理方法得到的异常应保持与原始曲线上“弱幅度”变化的一致性或与之基本吻合，至少由方法的物理意义分析应一致或吻合.

地电阻率时间序列曲线上出现的异常仅为数据异常，而数据异常并不等同于地震前兆异常，不应期望出现数据异常的台站附近及周围一定发生显著地震.出现数据异常后，应根据前兆异常的时空强及形态特征分析确认异常的属性，进而研判未来震情.