地电阻率是目前国内用于地震预报的定点前兆观测重要手段之一。随着地电仪的改进，台站数字化仪器稳定性有所提升，但地表的气象变化、人为施工等环境因素对地电阻率数据造成干扰，尤其是降雨量、地下水等干扰严重，常常在降雨后表现为地电阻率测值下降，且测值下降幅度随降雨量的增加而增大，随着地面降雨逐渐蒸发，地电阻率数值也逐渐恢复原来水平^[1-3]。因此为使地电阻率资料分析结果可靠且严谨，在跟踪分析各台站地电阻率数据时，有必要对地电阻率数据进行去降雨因素分析。去除降雨量对地电阻率影响的方法包括相关分析校正^[2]及褶积滤波法^{[1, 4-8]}，其中褶积滤波法应用效果最好且应用最为广泛。月速率法是描述地电阻率异常的常用因子之一，适用于中短期异常判定，在地电阻率异常分析中应用广泛^{[5, 9]}。

本文选取祁连山地震带地电阻率观测场地效果最佳的山丹台数据，利用褶积滤波去除降雨影响，依据地电阻率观测数据和台站附近中强地震震例确定归一化速率分析最优数据处理方式，并对比去除降雨前后归一化速率超限异常的映震效果。2022-01祁连山地震带中强地震活动进入空前活跃期，8月才逐渐趋于平静。结合此前中强地震震例总结得到的超限异常映震指标，分析8月之后山丹地电阻率的超限异常速率，对祁连山地震带地震活动预测有重要意义，同时为祁连山地震带地震危险性评价提供科学依据。

1 山丹观测站基本情况

山丹地震台位于山丹县城西8 km处，海拔1 748 m；台站观测场岩性为砂岩，测区地表为旷野型戈壁沙土堆，覆盖层较深，地表几乎没有植被。研究区发育河西走廊区域性主干大断裂，切割地壳较深，走向呈NWW向，活动性较强(图 1)。

山丹地震台地电阻率观测采用对称四极法(图 2)布极，布设有NS、EW、N45°W共3条测道。其中，NS测道和N45°W测道供电极极距AB=1.000 km，测量极极距MN=0.200 km，装置系数K=(π×AM×AN)÷MN=3.770 km；EW测道供电极极距AB=1.000 km，测量极极距MN=0.226 km，装置系数K=(π×AM×AN)÷MN=3.290 km。

2 褶积滤波去降雨分析

假设B是一个系统, 响应函数为H(t), t时刻输入X(t)，则输出Y(t)=H(t)×X(t)，二者的褶积为$\sum\limits_{\tau=0}^t H(\tau) X(t-\tau) $。本文将土层看作一个系统，降雨量C_R(t)作为输入信号，降雨量对地电阻率ρ_s的影响系统函数为R(t)，则降雨量对ρ_s的影响可表示为：

$ \begin{gathered} \rho_s(t)=R(t) \times C_R(t)= \\ \sum\limits_{\tau=0}^t R(\tau) C_R(t-\tau) \end{gathered} $

(1)

考虑到降雨对地电阻率的影响特点，R(t)为褶积算法的系统函数且满足式(2)特征，其中由于降雨下渗过程体现降雨的滞后效应，根据降雨渗透过程对地电阻率影响的研究结果^{[1, 6-7]}, 选择M=1 d，表示降雨即时影响时间为1 d，N=6 m表示降雨的滞后影响6个月：

$ \begin{gathered} &R(t)=\\ &\left\{\begin{array}{l} 0, t \leqslant 0 \\ A_0+A_1 t, 0<t \leqslant M \\ A_2+A_3 t+A_4 t^{-1}+A_5 t^{-2}, M<t \leqslant N \\ 0, t>N \end{array}\right. \end{gathered} $

(2)

联合式(1)和(2)可得ρ_s(t)= KX，其中，

$ \boldsymbol{X}=\left[A_0, A_1, A_2, A_3, A_4, A_5\right]^{\mathrm{T}} $

$ \boldsymbol{K}=\left[\begin{array}{cccccc} \sum\limits_{\tau=0}^M C_R\left(t_1-\tau\right) & \sum\limits_{\tau=0}^M \tau C_R\left(t_1-\tau\right) & \sum\limits_{\tau=M}^N C_R\left(t_1-\tau\right) & \sum\limits_{\tau=M}^N \tau C_R\left(t_1-\tau\right) & \sum\limits_{\tau=M}^N \tau^{-1} C_R\left(t_1-\tau\right) & \sum\limits_{\tau=M}^N \tau^{-2} C_R\left(t_1-\tau\right) \\ \sum\limits_{r=0}^M C_R\left(t_2-\tau\right) & \sum\limits_{\tau=0}^M \tau C_R\left(t_2-\tau\right) & \sum\limits_{\tau=M}^N C_R\left(t_2-\tau\right) & \sum\limits_{\tau=M}^N \tau C_R\left(t_2-\tau\right) & \sum\limits_{\tau=M}^N \tau^{-1} C_R\left(t_2-\tau\right) & \sum\limits_{\tau=M}^N \tau^{-2} C_R\left(t_2-\tau\right) \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ \sum\limits_{\tau=0}^M C_R\left(t_n-\tau\right) & \sum\limits_{\tau=0}^M \tau C_R\left(t_n-\tau\right) & \sum\limits_{\tau=M}^N C_R\left(t_n-\tau\right) & \sum\limits_{\tau=M}^N \tau C_R\left(t_n-\tau\right) & \sum\limits_{\tau=M}^N \tau^{-1} C_R\left(t_n-\tau\right) & \sum\limits_{\tau=M}^N \tau^{-2} C_R\left(t_n-\tau\right) \end{array}\right] $

选取多个时段的降雨数据和地电阻率整点值或日均值观测数据，给出ρ_s(t)的变化，得到降雨对地电阻率的影响特征曲线，经残差计算即可获得降雨对地电阻率的短期变化影响。再选用月均值通过褶积滤波去除降雨影响，便可分析降雨对地电阻率年变趋势的影响。选取观测站每年降雨量较为集中时段及集中降雨开始前1 d各观测站不同测道观测日值，利用褶积滤波去除短期集中降雨对地电阻率的影响，再分析降雨量对山丹电阻率年变和趋势性变化的响应。

2.1 降雨对地电阻率的短期变化影响

选取2014~2022-08山丹观测站同台观测的气象三要素观测资料中降雨量日均值数据及2014-06-11~07-22、2014-08-02~09-11、2015-05-27~07-22、2015-09-02~21、2016-04-15~05-01、2016-05-11~06-03、2016-06-20~07-21、2016-08-10~29、2017-04-07~24、2017-05-03~21、2017-06-03~07-31、2017-08-07~09-04、2017-10-17~31、2018-06-16~07-26、2018-08-04~09-17、2018-10-26~30、2019-04-08~05-09、2019-06-09~07-23、2019-08-05~09-19、2020-05-05~06-22、2020-07-17~09-21、2021-04-13~05-22、2021-06-06~07-21、2021-08-02~09-14、2022-05-28~06-20、2022-07-13~08-28共26个时段的降雨数据和降雨时段前1 d NS测道、EW测道及NW测道地电阻率数据，分别给出系数矩阵K和降雨影响地电阻率变化量ρ_s(t)，根据ρ_s(t)= KX得出参数矩阵X，再将X作用于山丹地电阻率历年降雨观测数据日均值，给出降雨对地电阻率观测影响量ρ_s(t)的短期变化，结果如图 3所示。图中，蓝色代表原始曲线，红色代表通过褶积滤波法去除降雨对各观测站地电阻率影响后的残差变化，可以看出，降雨多集中在夏季。山丹地电阻率EW测道和NW测道褶积拟合基本与实际观测曲线重合，NS测道拟合后残差数值高于实际观测值，说明短期降雨可能造成地电阻率降低，降雨结束后测值恢复。

2.2 降雨对地电阻率的季节性影响

选取2014~2022-08山丹观测站气象三要素中降雨量月均值数据，利用褶积滤波计算降雨对山丹地电阻率年变和趋势变化的响应，结果如图 4所示。由图可知，降雨对电阻率季节性的最大影响出现在每年最高降雨之后，山丹地电阻率在剔除月降雨量对低电阻率的影响值后，原始观测曲线“冬高夏低”的年变规律均被削弱，但不改变长趋势变化形态，因此地电阻率去除降雨季节性周期影响之后变化曲线较为平稳，年变形态相比原始观测曲线不再清晰，趋势变化依旧存在，说明降雨量对地电阻率的年变形态有一定的影响，但对趋势性变化影响不大。

3 归一化月速率分析

月速率分析的原理为：设定各个观测站地电阻率观测序列ρ={ρ₁, ρ₂, …, ρ_N}，以滑动步长计算归一化月速率的序列为$S_i=\frac{R_i \times K_i}{\sigma_{n-1}}(i=n, $n+1, …, N)，其中R_i为序列曲线斜率，K_i为自相关系数，分母为分子乘积关于时间序列的均方差，N为选取资料长度，滑动步长n视情况选定。若最终处理结果S_i绝对值超过2.4，则视为异常^{[5, 10-13]}。

采用GSEP软件^[14]进行归一化月速率分析时，数据选用月均值，推荐滑动步长为4~12，优势滑动步长为6~8^[15]。采用不同滑动步长分别对去干扰数据进行归一化分析发现，选用步长越大，归一化异常形态越不显著，滑动步长为6~8时，异常形态基本一致，但随着步长增大，异常超限部分变小，从异常识别和指标提取的角度，选取步长6更为合理。月相对变化速率进行地电阻率前兆异常提取的关键是对地电阻率资料进行去倾、去年变及排除周期变化，因此在分别对不同测道数据采用去年变、排除周期变化和去倾3种不同处理方式单独及组合处理数据后再进行月速率分析，并选择最优分析方法。

4 去降雨后地电阻率映震效能分析

选取2014~2021年山丹地电阻率月均值及祁连山地震带10次5级以上地震震例，包括2014年青海乌兰5.1级、2015年青海祁连5.2级、2016年青海门源6.4级、2017年内蒙古阿拉善左旗5级、2019年青海茫崖5级、甘州5级、2021年青海茫崖5.8级、阿克塞5.5级及2017年九寨沟7级、2021年青海玛多7.4级地震，对山丹地电阻率原始数据和去除降雨后数据分别进行归一化速率计算。月速率分析方法是异常分辨率较高的中短期预测方法^[14-15]，超限异常的预测有效时间自超限异常开始第2个月起，到异常结束后6个月止，根据异常预测时段和发震时间分析地震与超限异常的对应情况。表 1分别统计了异常结束后半年内有震次数(发震)、异常结束后无震次数(虚报)、震前半年内无异常次数(漏报)，结果显示，去除降雨后数据异常映震次数显著增多，异常的发震率增加，虚报率降低，多种组合方式获得的异常对应结果基本一致。

对比原始数据和去除降雨后数据归一化结果，NS测道和NW测道在去除降雨影响后，经归一化速率分析的超限异常对应地震效果有所下降，EW测道去除降雨后的归一化速率分析超限异常映震效果显著增强。对于原始数据，NS测道仅需去周期处理即可有较好的映震效能，而去除降雨后需同时去除年变和周期才能有较好的映震效能；EW测道去除降雨前后的数据在进行去年变+去周期的相同处理后，剔除降雨的映震效能更好，说明EW测道受降雨影响更为显著；NW测道剔除降雨前后的数据在进行归一化分析时的最优处理方式均为去年变+去周期，但明显在去除降雨后虚报率有所下降，说明对于山丹地电阻率而言，降雨对地电阻率的异常提取产生了一定影响。去降雨之后山丹地电阻率进行归一化月速率分析最优结果为：去除降雨后的山丹地电阻率NS测道归一化处理既可选择去倾+去年变，也可选择去倾+去周期；EW测道和NW测道归一化处理则选择去周期的方式。

将山丹地电阻率数据研究时段和震例分析时段均延长至2022-08，根据选取的归一化前最佳处理方式对山丹地电阻率进行计算，结果如图 5和6所示，超限时段和有效预测时间内发震情况统计结果如表 2所示。结果发现，NS测道去倾+去周期后共有4次超限异常，在有效预测期内对应3次地震，去倾+去年变后共有5次超限异常，有效预测期内对应地震1次；EW测道去降雨+去周期后共有8次超限异常，有效预测期内发震12次，其中2020年年底出现的超限异常时段内共发生3次地震，2021年底超限异常时段内共发生5次地震；NW测道去降雨+去周期后共出现超限异常3次，发生4次地震，其中2021年年初超限异常有效预测期内发生3次地震。除2014~2021年10次地震，2022年门源6.9级地震前，EW测道也出现超限异常。

山丹地电阻率3个测道中映震效果最好的是EW测道，NS测道和NW测道映震效果相当，测道间异常变化特征并不相同，表现出各向异性，与地电阻率观测场地的地层岩性及层状岩石含水量影响有关^[16]。因此归一化速率分析结果显示，并不是所有地震前3个测道均能同步出现异常特征，且由于各测道归一化前数据处理方式不同，映震效果差异也较大。地电阻率各向异性变化存在一个高值易发震变化区间，对台站周围500 km范围内5级以上中强地震具有明显的异常特征^[17]，而满足300 km范围内5.0≤M_S≤5.9、400 km范围内6.0≤M_S≤6.9及500 km内M_S≥7.0的地震视为与异常对应^[18]，由此看来，表 2中山丹地电阻率NS测道在茫崖5.0级地震前、EW测道在茫崖5.0级、5.8级及九寨沟7.0级地震前出现的异常属于虚报，并非前兆异常。

去除超出500 km范围的地震后，山丹地电阻率NS测道经去倾+去周期处理后映震效果较去倾+去年变好。山丹地电阻率NS测道和EW测道均在甘州5.0级、玛多7.4级地震前出现超限异常且与其对应，NW测道也在玛多7.4级地震前出现异常且对应，EW测道和NW测道则在阿克塞5.5级地震前同时出现超限异常且与地震对应。EW测道超限异常持续时间为4~6个月，且对应多次中强地震。山丹地电阻率超限异常对祁连山地震带的预测指标如表 3所示，其中NS测道速率超限异常持续1~2个月时，在异常开始至异常结束后半年内存在发生5级左右地震的可能；EW测道速率超限异常持续4~6个月时，异常开始至结束后半年内祁连山地震带存在发生6级左右地震的可能；NW测道速率超限异常持续1个月时，异常开始至结束后半年内祁连山地震带存在发生5~6级地震的可能。

5 结语

通过对2014~2022-08山丹地电阻率观测数据进行褶积滤波去除降雨影响后，利用归一化月速率分析结果的映震效果显著好于原始数据映震效果，根据去除降雨后的电阻率数据映震指标可知，山丹地电阻率3个测道的结果均可作为祁连山地震带中强地震前的短期异常指标，根据震例情况，对台站周边300 km范围内5~6级地震的指示意义更强。

由2022-08~12山丹地电阻率各测道速率超限异常可知，EW测道2022-08以来持续5个月出现超限异常，NW测道2022-12出现超限异常，说明2023年祁连山地震带中西段发生中强地震的危险性仍未得到缓解，需重点关注速率超限异常。后期应结合震例指标情况，根据速率超限异常结束时间来判定未来发震的优势时间，同时密切关注研究区地震活动性及应力演化特征，若出现地震活动显著增强及应力积累特征，则可进一步验证山丹地电阻率指标的可靠性及可行性。

致谢: 感谢河北省地震局张国苓高级工程师和陕西省地震局石富强副研究员为本文褶积滤波降雨分析给予指导和帮助。