摘 要：提出了地电阻率滑动自相关数据处理方法，利用随机时间序列对方法进行了数值检验，验证了方法的有效性。将滑动自相关方法应用于3个地电阻率台站实际观测资料，经分析研究得到了一些比较有意义的结果： 1)甘孜台和山丹台地电阻率滑动相关系数变化与地震事件有较好的对应关系，主要表现为，在震前半年到1年时间范围内，会出现相关系数阶段性增大现象，同时，这2个台站的相关系数异常也表现出比较强的各向异性特征； 2)成都台地电阻率观测虽受干扰较为严重，但相关系数异常现象与地震事件仍有比较好的对应关系，除了说明所使用方法的有效性之外，还可能与参加计算的地震事件震级比较大，震中距较小，以及发震时间比较集中有关，在该台站同样也存在异常的各向异性特征； 3)通过对比不同震级特征地震的相关研究结果得出，当选择的特征地震震级较小时，不同地震的震前异常幅度会有所不同，但当震级较大时(如MS≥5.0)，其对研究结果产生的影响非常有限。此外，还对地电阻率地震异常的各向异性，以及特征地震的选取问题进行了简单论述。
关键词：滑动自相关方法 地电阻率 地震异常 相关系数 各向异性
0 引言
将地电阻率观测资料应用于地震预测是中国的一大特色。地电阻率观测始于1966年邢台MS7.2地震，其观测装置为对称四极观测(杜学彬等，2006，2007； 钱家栋等，2009)，经历了近50a的观测积累和分析研究，得到了与大震和强震突出相关的多次异常现象。如： 对1976年7月28日唐山7.8级地震进行的地电阻率相关研究(赵玉林等，1978； 钱复业等，1990； 钱家栋等，1993； Lu et al.，1999)，以及对2008年5月12日汶川MS8.0和2013年4月20日芦山MS7.0地震地电阻率观测开展的相关分析研究(张学民等，2009； Huang，2011； 钱家栋等，2013； 朱涛，2013； 杜学彬等，2015)，都得到了非常好的结果。总体来看，大地震前的地电阻率变化特征一般表现为： 震前2～3a开始出现趋势下降异常现象； 震前地电阻率异常具有阶段性，个别大震前还有短临和临震前兆现象； 震前趋势异常表现出各向异性(钱家栋等，1993)等。
地电阻率地震异常特征在1993年就已经得出(钱家栋等，1993)，但2～3a的时间预测范围，所能起到的减灾实效比较有限。到目前为止，针对地电阻率观测资料在实际地震预测中进行应用的方法研究仍显乏力。杜学彬等2001年提出的归一化月速率方法，推进了地电阻率观测手段在预测预报中的实际应用，并开展了一些较好的分析工作(叶青等，2005)。
以往开展的地电阻率观测异常分析，很多是关于特定地震事件进行的回溯性研究。虽然也得到了较好的对应结论，但这些分析研究只是增加了典型震例的数量，对地电阻率观测在实际应用工作的推动有限。本文借鉴地震电磁异常滑动相关分析方法(Jiang et al.，2016)的研究思想，提出了滑动自相关地震预测方法，并初步应用于四川省甘孜、 成都、 甘肃省山丹地电阻率台的实际观测资料，进行一些比较有意义的探索性分析研究。
1 方法简介及数值检验
1.1 方法简介
依据经验性地震预测的特点，假设在观测台站一定震中距范围内发生的地震事件，所引起的地电阻率异常现象具有一定的相似性。根据这一假设条件，为突出地电阻率地震异常的中期和短期特点，文中选取某次地震(下文称 “特征地震”)前一段时间(180d)的观测资料作为特征时间序列，通过滑动特征时间序列与原始时间序列得到相关系数曲线。利用相关系数曲线与地震事件序列的对应性，分析研究震前相关系数的变化特点，从而对未来的地震事件进行预测。相关系数的计算公式：
(1)
式(1)中，X1和Y是2个长度相等的时间序列，RME 是这2个离散时间序列的相关系数，
1和
是相应的平均值，m是采样点数。
计算分析过程中采用地电阻率日均值。与地震事件进行相关性计算分析时，应考虑1d内距离台站一定范围内所有地震事件的总体响应，借鉴Jiang等(2016)的方法，进行等效地震计算。首先，利用式(2)计算出1d内所应考虑地震事件的总释放能量(Gutenberg et al.，1956)，再由式(3)计算出相应的等效震级(Jiang et al.，2016)。
(2)
(3)
式(2)、(3)中，E为同一天内地震释放的总能量，Mi为第i个地震的震级，s为1d内的总地震数量，A和B为常系数，分别为1.96和2.05，
为等效震级。
1.2 数值检验
为了验证方法的有效性，利用随机时间序列对该方法进行数值检验。首先，需要产生1段随机时间序列，选取其中1段时间序列作为滑动时间序列； 其次，将滑动时间序列，按不同的时间间隔叠加在原始时间序列上，从而形成新的时间序列，模拟地电阻率观测中震前异常具有相似性的特点； 最后，通过滑动所选取的特征时间序列，计算与新时间序列之间的滑动相关系数，通过比较不同时间相关系数的大小，达到验证该方法有效性的目的。具体过程如下：
图1 数值检验
Fig. 1 The demonstration of numerical examination.
纵坐标为数据大小； a 随机时间序列； b 随机时间序列中的1段数据； c 新合成的时间序列； d 相关系数图
信号X(t)(图1a)是随机信号，X1(t)(图1b)是从X(t)中随机选出的1段信号，长度为100个数据点。将X1(t)分别叠加在501～600，1,201～1,300，1,501～1,600和1,901～2,000区间上，从而形成新的时间序列Y(t)(图1c)，通过滑动X1(t)与Y(t)做相关性计算，从而得到相关系数R曲线(图1d)，由于开始滑动的位置为1，所以结束时的滑动位置为2,000，故相关系数只有2,000个。从图1d中可以看出，在101位置时相关系数为1，表示为完全相关，符合实际情况； 在其他4个插入滑动信号的地方，其相关系数都在0.6以上，呈现出很强的相关性，从而证实该方法的正确有效性。
2 实测数据分析
分析研究过程中使用地电阻率观测日均值，计算过程中需对观测资料进行一些预处理。首先，计算观测台站每天观测值的均方根误差； 其次，对>3倍均方根误差的观测数据进行删除，然后计算出相应的日均值。此外，在滑动计算时采用地电阻率日均值的一阶差分(前差)，可以降低大阶跃和其他干扰带来的影响，使计算结果更为理想。分析过程中所使用的地震目录从中国地震台网中心网站下载，目录截止日期为2014年11月30日。
2.1 甘孜台结果分析
甘孜地电观测台站位于四川省甘孜县城郊区东南方向斯俄乡境内的布绒郎山下，现有2个观测方向(N30°E和 N60°W)，本次分析研究所使用观测资料的时间段为： 2007年9月26日至2015年1月26日。参加滑移的特征地震事件为2013年1月18日发生的MS5.5地震(30.95°N，99.4°E)，考虑该方法对地电阻率短期和临震异常的实用性，参加滑动的特征长度取180d，参加等效地震计算的地震事件选取范围： MS≥6.5地震的震中距范围≤350km，MS 在5.0～6.5地震的震中距范围≤250km，MS<5.0地震的震中距范围≤100km。
图2 是四川省甘孜地电台的滑动自相关结果图。从图中可以看出，地震事件与震前的相关系数异常有较好的对应关系，结果初步显示，在震前半年到1年时间范围内，会出现相关系数增大异常现象，并且这种异常现象有一定的阶段性，并不是持续发生； 从这几次的对应情况来看，N60°W测向比N30°E效果好，这可能与地震震中相对于台站的位置有一定的关系，同时，也是地电阻率震前异常各向异性的1个具体体现。由于该台所选时间范围内，地震具有成簇发生的特点，导致多次震前异常与震后效应可能相互叠加，这个特点可能给震前异常的判断带来一定的困难。
图2 甘孜台滑动自相关结果
Fig. 2 The results of shifting self-correlation at Ganzi station.
虚线表示等效地震的发生时间位置； 相关系数达到1.0时，其所对应的地震事件为该台站选定的特征地震
2.2 成都台结果分析
成都地电观测台站位于四川省都江堰市与郫县交界的走石山，现有2个观测方向(N58°E和 N49°W)，本次分析研究所使用观测资料的时间段为： 2010年10月29日至2015年1月26日。参加滑移的特征地震事件为2011年9月4日四川发生的MS4.2地震(31.27°N，103.62°E)，参加滑动的特征长度为180d，参加等效地震计算的地震事件选取范围： MS≥6.5地震的震中距范围≤350km，MS 在5.0～6.5地震的震中距范围≤200km，MS<5.0地震的震中距范围≤100km。
图3 是成都地电阻率观测的处理结果图。台站受工业游散电流、 成青快铁等干扰较为严重，影响了传统处理方法的正常使用，但滑动自相关方法仍能得到较理想的结果。从结果图中可以看出，相关系数异常与地震事件有比较好的对应关系，利用滑动自相关方法处理得到的异常现象，出现在震前大概半年时间内，其异常时段的相关系数值与甘孜台相比要大许多，这种现象可能是由于，在分析研究的时间段内所选地震事件的震级比较大，而且相对比较集中发生的原因； 此外，还有震中距较小的原因，如芦山MS7.0地震距离成都台约100km。该台站N49°W测向对四川芦山MS7.0地震的异常反映比N58°E好，这种现象与王成虎等(2014)得到的该区域为逆冲断层应力状态，最大水平主压应力方向(N44°～64°W)相一致，说明地电阻率异常与主压应力之间可能存在一定的关系。
图3 成都台滑动自相关结果图
Fig. 3 The results of shifting self-correlation at Chengdu station.
2.3 山丹台结果分析
山丹地电观测台站位于甘肃省山丹县清泉镇南湾村，现有3个观测方向(SN、 EW和 N45°W)，本次分析研究所使用观测资料的时间段为： 2007年5月11日至2015年1月26日。参加滑移的特征地震事件为2008年3月30日青海发生的MS5.2地震(37.97°N，101.92°E)，参加滑动的特征长度为180d，参加等效地震计算的地震事件选取范围： MS≥6.5地震的震中距范围≤300km，MS5.0～6.5地震的震中距范围≤200km，MS<5.0地震的震中距范围≤100km。
图4 为甘肃山丹地电阻率观测的处理结果图。该台站EW测向映震效能较好，SN和 N45°W测向次之，相关系数异常现象多发生在震前半年至1年时间范围内。其中SN测向表现出的相关系数变化形态较为特殊，其原因可能是由于该测向对特征地震的反映比较弱，参加滑动的特征时间序列包含的异常信息非常少所导致。但这种现象在EW和 N45°W 测向上却表现得截然相反，这种现象从侧面反映出，山丹台存在明显的地电阻率各向异性现象。EW和 N45°W 测向在2015年初出现的这种相关系数异常现象，或许与该区域未来一定时期内某次MS5.0左右的地震相对应。
图4 山丹台滑动自相关结果图
Fig. 4 The results of shifting self-correlation at Shandan station.
3 相关讨论
图5 不同特征地震得到的相关系数图
Fig. 5 The correlation coefficient map of different characteristic earthquakes.
在分析研究的过程中，特征地震(Meq)的选取比较关键，不同震级的地震事件在震前引起的异常变化不会完全相同。针对这种可能情况，以甘肃省山丹地电阻率台为例，分析讨论了参加滑动的特征地震不同时会出现何种变化情况。图5 是山丹台地电阻率选取不同特征地震时得到的相关系数图，图5a是Meq震级为3.6的结果，其震中距约118km； 图5b是Meq震级为3.8的结果，其震中距约54km； 图5c是Meq震级为5.3的结果，其震中距约95km。初步研究发现，参与滑动的特征地震震级与等效地震震级相差不大时，其相关系数异常较为明显，如图5a和图5b中震前相关系数异常幅度均比图5c的大； 当两者相差较大时，仍有相关系数异常现象出现，但异常幅度较小(相关系数值较小)。针对不同的特征地震，所得到的相关系数图存在一定的相似性，反映了假设条件的正确性和所采用方法的有效性。
从文中3个台站观测资料的分析能够看出，不同测向对地电阻率地震异常的反映各不相同，除了有震级和震中距的因素以外，可能存在着地电阻率各向异性的原因，关于这个现象，杜学彬等(2007)研究了强地震附近视电阻率各向异性变化的原因，通过建立真、 视电阻率各向异性变化与裂隙率v、 骨架电阻率r、 饱和裂隙电阻率的本构关系，解释了视电阻率各向异性的原因是，在强地震孕震晚期阶段，在震源区及附近地壳近地表的较深部分，介质微裂隙发育，其走向沿最大主应力方位优势取向，低阻水填充微裂隙，导电通道连通，引起最大主压应力方向的真电阻率变化最为显著，导致电阻率的各向异性变化。文中各测向滑动自相关系数对地震事件的反映能力的不同，只是地电阻率各向异性的一个侧面体现。
4 结论
本文基于距离地电阻率观测台站一定范围内发生的地震，在震前会出现相关性较高的异常现象这一前提假设条件，进行随机时间序列的滑动自相关数值检验，其结果说明，通过相关系数的变化能够对特定的时间序列进行检出，从而为后续特征地震序列的检测提供可能。随后，将滑动自相关分析方法，应用于3个地电阻率台站实际观测资料的分析研究，得到了一些比较有意义的结论：
(1)甘孜台地电阻率分析结果显示，地震事件与震前相关系数异常有较好的对应关系。主要表现为，在震前半年到1年时间范围内，会出现相关系数增大的异常现象，并且这种异常现象有一定的阶段性，并不是持续发生。从几次地震的对应情况来看，该台 N60°W 测向比N30°E向效果好，除了与地震震中相对于台站的位置有一定的关系外，非常可能也是地电阻率震前异常各向异性的具体体现。
(2)成都地电阻率台站虽受工业游散电流、 成青快铁等干扰较为严重，但相关系数异常仍与地震事件有比较好的对应关系。但该台站异常时段的相关系数值与甘孜台相比要大许多，这种现象可能是由于，在分析研究的时间段内所选地震事件的震级较大，震中距较小，以及发震时间比较集中所导致。此外，该台站同样存在异常的各向异性特征。
(3)与甘孜台相类似，山丹台地电阻率相关系数的增大异常也发生在震前半年至1年时间范围内。此外，该台站也存在明显的地电阻率异常各向异性现象。
(4)参与计算的特征地震震级对所得到的结果有一定的影响，当选择的特征地震震级较小时，对震前异常的判断会有一定的影响，但当震级较大时(如MS≥5.0)，其对研究结果产生的影响非常有限。 MS≥5.0的地震正是我们非常关心的地震事件，更反映了本文所提方法的实用性。
从上述研究结果可看出，滑动自相关方法应用于地电阻率观测资料，能够得到比较好的结果，该方法为提高地电阻率观测资料在地震预测中的应用率提供了新的思路。文中所涉及3个台站都存在异常的各向异性现象，研究出现异常的方位与震中的位置关系，将是进一步工作的重点。未来可以对更多观测台站、 更长观测资料进行分析研究，还可以针对不同观测台站，从统计学的角度分析研究属于各个台站的特征地震异常和特征异常时间尺度，逐步提高地电阻率地震异常的判定能力。
致谢 本文在中国地震局地质研究所完成，对中国地震局人事教育司和甘肃省地震局人事教育处提供的相关帮助，以及审稿专家的有益建议表示衷心感谢!
