开展地下介质波速变化测量研究对于认识地震孕育机理和探索地震物理预报手段具有重要意义。随着背景噪声层析成像技术的发展，经验格林函数也被应用于地壳内部地震波速度的时空变化特征研究^[1]。背景噪声分布广且不依赖于特定震源，同时具有可重复性和经济环保等优势，非常适合跟踪地壳内部结构的物性变化^[1]。

背景噪声较天然震源具有更高的时间采样率，也比人工震源更加经济。通过背景噪声波速方法，地震学家观测到火山区波速的季节性变化^[2]，发现火山喷发前数周内波速具有下降特征^[3]。部分学者^[4-10]对同震和震后地震波速度的变化特征进行了研究，但大部分都基于波形资料来反演强震发生过程及震后的波速变化。

近年来，滇西北地区构造运动强烈，震情形势严峻。本文以滇西北作为实验地点，基于部分固定台站记录的宽频带连续波形资料，采用背景噪声互相关方法研究地壳介质内部地震波速随时间的变化特征，通过对台站对走时偏移时间序列进行深加工，提取适合于滇西北地区M≥5.0地震的短临异常识别指标。该工作的开展对深入研究滇西北地区地下介质的变化过程，探索测震资料及背景噪声技术在日常震情跟踪及地震预测中的应用，提高异常识别能力及推进指标体系建设等具有重要意义。

1 数据收集与整理

本文主要研究滇西北地区(25°~28°N，99°~101.5°E)波速变化，选取云南数字地震台网中5个固定台站的垂直分量连续波形资料，时间范围为2012-01~2020-11，该时间段内所有台站的数据资料连续性都很好，数据基本能够反映介质波速的变化。研究时段内，研究区共发生6次M≥5.0地震，图 1为5个地震台站、2个GNSS连续观测站及6次地震震中分布情况，相关地震参数见表 1。

2 背景噪声波速测量方法 2.1 单台站数据预处理

采用滇西北地区5个宽频带地震台记录的垂直分量连续波形数据计算地震波速度。首先将连续波形数据按天截取^[10]，采用1 Hz对波形数据进行抽样，以降低数据采样率；然后进行去除仪器响应、去均值、去趋势、带通滤波处理^{[1, 8-10]}；接着对波形数据进行时域平均，以去除地震事件的影响；最后进行频谱白化处理^{[8, 10]}。以上数据处理的目的是为了消除天然地震和仪器本身异常信号的影响，以获取高质量的地震背景噪声。

2.2 背景噪声互相关及数据叠加

单台数据预处理后，对由5个台站构成的10个台站对每天的背景噪声进行互相关处理，以获取各台站对单日的经验格林函数^[8]。Stehly等^[11]研究认为，对任意2个台站记录的背景噪声进行互相关处理，理论上得到的经验格林函数中有一正、一负2个分支，分别表示台站对路径上的因果和非因果信号。但背景噪声的来源及能量差异，会极大程度地影响经验格林函数的形态^[12]。当台站两侧的噪声源分布均匀时，因果信号和非因果信号的到时一致，振幅相同；而当噪声源分布不均匀时，2个方向信号的到时相同，但振幅不同，在噪声源能量较强一侧产生的信号振幅较大^[11-12](图 2)。

Bensen等^[13]研究认为，对因果和非因果信号进行叠加处理可以等效噪声源的均匀分布。虽然该处理会损失振幅信息，但可以在一定程度上提高信号的信噪比^[12-13]。因此，本文以0为界线，将经验格林函数的“负支”反号后与“正支”进行对称叠加，形成“叠加分量”(图 3)。

由于受噪声源的不均匀分布和噪声成分的不确定性等因素影响，同时台站间距离较远，导致台站对单日经验格林函数的信噪比普遍较低，难以分辨地震波形态^[10]。为提高信噪比，本文将各台站对单日前59 d(共60 d)的经验格林函数进行叠加，作为当天经验格林函数(图 4(b))。在此基础上，通过进一步滑动叠加，获取2012-02-29~2020-11-16当天经验格林函数(前59 d因叠加数据不足60 d，不作进一步分析)。

为获得台站对经验格林函数随时间变化的定量信息，需确定各台站对的参考经验格林函数^[8-10]。Stehly等^[11]研究发现，背景噪声源的季节性变化会引起由互相关方法重构的面波信号走时变化，因此建议尽量采用较长时间范围(至少1 a)的资料来确定参考经验格林函数。本文将各台站对2012-02-29~2020-11-16共3 184 d的单日经验格林函数进行叠加，作为各自的参考经验格林函数(图 4(a))。

2.3 走时偏移提取

在获取台站对当天经验格林函数和参考经验格林函数后，对当天经验格林函数和参考经验格林函数进行互相关计算^{[8, 10]}，获取相似波形窗口的互相关系数(图 4(c))，互相关系数最大时对应的时间延迟就是这两段相似波形窗口的走时偏移。对计算结果进行余弦插值以提取更高精度的走时偏移。以互相关系数为约束(互相关系数要求在0.9以上)，通过提取直达瑞利波中部分震相的走时偏移来计算台站对之间瑞利波的速度。图 5(a)为TUS-YOS台站对瑞利波走时偏移变化曲线。

从图 5(a)可以看出，TUS-YOS台站对瑞利波走时偏移变化曲线存在明显的年变成分。为消除年变成分对走时偏移的影响，采用傅立叶变换以获取长趋势变化背景下隐含的波速变化特征。对于去除年变成分后的走时偏移，划定±1.5倍标准差作为单条曲线的异常指标(以±1.5倍标准差作为异常阈值)。从图 5(c)可以看出，该台站对的走时偏移在洱源M5.5和M5.0、云龙M5.0、漾濞M5.1等地震前存在较为明显的超指标异常。

2.4 台站对走时偏移可靠性分析

由图 1可知，地震台站YUL和TUS附近分别有云龙和下关2个GNSS连续观测站，且TUS-YUL台站对的路径分布与云龙-下关GNSS基线基本保持平行。图 6为TUS-YUL台站对瑞利波走时偏移与云龙-下关GNSS基线相对变化的对应关系(GNSS数据获取时间截至2017-05-20)。从趋势上来看，TUS-YUL台站对瑞利波走时偏移与云龙-下关GNSS基线相对变化存在较好的耦合性，即在一段时间内，当TUS-YUL台站对瑞利波走时偏移缩短(波速加快)时，云龙-下关GNSS基线长度也在缩短(挤压增强)。由于本文提取的数据为台站对直达瑞利波走时偏移，反映的是中上地壳物性随时间的变化特征，而GNSS基线监测的是中上地壳的应力变化，两者具有较好的耦合性，表明本文提取的瑞利波走时偏移能够监测到中上地壳介质的物性变化，这也间接表明本文的数据结果具有可靠性。

3 短临异常指标提取及效能评价

由于需要获取短临异常信息，地震事件对应的时间间隔应不超过90 d。对10个台站对的瑞利波走时偏移时间序列进行傅立叶滑动去年变处理，并设定±1.5倍标准差作为异常阈值，结合研究区6次地震对各台站对走时偏移时间序列的映震能力进行效能评价。现阶段应用最为广泛的预报效能评价方法为R值评分法^[14]，其主要依据预测指标实际预报地震的有震报准率c和预报占时率(预报时间占有率)b^[15]，可反映预测方法或指标与地震的相关程度。R值计算公式^[15]为：

$ R=c-b=有震报准率-预报占时率=\\ \frac{\text { 预报有震命中次数 }}{\text { 地震总次数 }}-\frac{\text { 预报有震时间 }}{\text { 预报研究总时间 }} $

(1)

式中，R=1表示全报准；R=0表示预报未起作用。R值越大，预报效果越好。根据R值评分规则，在97.5%置信度下，当预报效能R＞R₀(预报效能检验的最低值)时^[16]，表明所评估的预测方法可通过统计检验，高于随机预测的预报效能。

表 2为10个台站对瑞利波走时偏移时间序列的预报效能检验结果。

由表 2可知，在10个台站对中，有7个台站对通过了检验，3个台站对未通过检验，但通过检验的7个台站对的对应率普遍偏低(各台站对均存在较高的虚报)。考虑到单个台站对的虚报率偏高，且可能存在不确定性等因素，需要将通过检验的7个台站对进行综合分析，以获取稳定可靠的地震短临异常识别指标(综合指标)。本文采用自适应加权综合预测方法获取综合概率指标，计算公式为^[17]：

$ P_{t}=\frac{\sum\limits_{t=1}^{n} K_{t i} P_{t i} \lg N_{t i}}{\sum\limits_{t=1}^{n} P_{t i} \lg N_{t i}} $

(2)

式中，n为预报因子总数；K_ti为到时间t为止的第i项因子动态值(当第i项因子出现异常后，其所对应的地震危险时段K=1，其余时段K=0)；P_ti为到时间t为止的第i项因子历史映震率；N_ti为到时间t为止的第i项因子历史映震次数；P_t为时间t的综合概率值^[17-18]。

从理论和实际预测中均可发现，在各项因子与地震的相关性及预测贡献中，最重要的表征值为映震率P_ti。另外，各因子的映震次数是因子与地震相关性置信度评价最主要的因素，通过置信度检验的相关因子，映震次数N_ti越大，可信度就越高^[18]。

图 7为采用自适应加权综合预测方法对通过检验的7个台站对进行综合分析后获取的地震短临异常识别指标(综合指标)。

在获取滇西北地区的综合指标后，设定2倍均值线作为异常阈值(此处为0.60)，当综合概率值P_t超过2倍均值线时，表示出现异常，认为滇西北地区未来90 d内存在发生M≥5.0地震的危险。

在提取滇西北地区的综合指标和划定2倍均值作为异常阈值后，需要结合对应的地震事件对综合指标进行效能评价(效能评价采用R值评分法)。预报效能检验结果表明，利用该综合指标对滇西北地区2012年以来发生的6次M≥5.0地震进行90 d短临预报时，异常指标共出现8次(2016年异常指标超过90 d，故划分为2次)，其中准确预报地震5次，漏报1次，虚报4次，预报效能评分R为0.692，R₀为0.475，R>R₀，效能检验通过。该综合指标的地震对应率为62.50%，概括率为83.33%，表 3为综合指标异常出现时间与地震发生时间的对应情况。

由表 3可知，大部分地震的发生与异常指标的出现时间间隔均不超过45 d，表明该综合指标在时间上具有较好的短临预报能力。结合图 1中台站对路径及地震震中分布可知，该综合指标对台站对覆盖范围内的洱源M5.5、洱源M5.0、云龙M5.0和漾濞M5.1地震都作出了准确的短临预报，不存在漏报现象，从而检验了该指标的预测能力。此外，该指标对台站对范围之外的昌宁M5.0地震也可作出一定预报，表明该指标对台站对覆盖范围及邻区的地震均有较好的预报能力。

4 讨论

通过背景噪声互相关计算可获取台站对的经验格林函数。研究表明，噪声源往往呈不均匀分布，且具有一定的季节性变化，这与理论假设存在差异^[12]。该差异使得通过互相关计算获取的经验格林函数并不完全收敛于台站间的理论格林函数，并存在季节性变化，从而导致获取的走时偏移(波速变化)受季节变化影响存在明显的年变特性。Wang等^[19]基于日本台站数据对地震波的季节性变化特征进行研究发现，地震波速度存在夏季低、冬季高的特点。根据波速和走时偏移的计算关系^[3]dv/v=-dt/t可知，本文计算得到的瑞利波走时偏移dt存在夏季高、冬季低(对应波速dv为夏季低、冬季高)的变化特征，与Wang等^[19]的结果一致。为消除季节性变化对走时偏移的影响，本文对走时偏移时间序列进行傅立叶滑动去年变处理，该方法理论基础严密，是一种线性、无偏、最优的周期估计方法，利用该方法去除年变成分的可行性已得到认可^[20]。傅立叶滑动去年变方法能最大限度地去除年变特性对走时偏移的影响，因此认为本文提取的短临异常识别指标(综合指标)具有可靠性。

本文基于6个震例，采用R值评分法对台站对走时偏移时间序列的映震能力进行检验发现，有7个台站对通过检验，3个未通过检验。究其原因，可能与强震发生区域的波速变化存在一定的地域性，且不同台站对之间的面波传播路径存在差异有关。具体原因还有待进一步研究。

2021-05-21在本文台站对覆盖范围内发生漾濞6.4级地震，为验证本文提取的短临异常识别指标(综合指标)是否对该地震的发震时刻具有指示意义，对综合指标进行数据更新，结果如图 8所示。

由图 8可知，该短临异常识别指标(综合指标)在2021-04-20(距漾濞6.4级地震发震时刻32 d)出现明显的超指标异常，因此漾濞6.4级地震的发生是对本文提取的短临异常识别指标(综合指标)的很好检验。

在上述分析的基础上，对2021-04-20异常指标出现时的异常台站对进行提取，以探讨异常台站对的分布与漾濞6.4级地震震中位置的对应关系，结果如图 9所示。从图中可以看出，2021-04-20综合指标出现异常时，震中附近DAY-TUS、DAY-YUL和YUL-ZOD台站对的走时偏移出现超指标异常，漾濞6.4级地震发生在3个异常台站对的附近区域。由此可见，基于综合指标提取的异常台站对分布对未来地震的发震地点具有一定的指示意义。

5 结语

本文基于地震台站记录的宽频带连续波形资料，采用背景噪声互相关及傅立叶变换等方法，获取台站对瑞利波走时偏移。通过对台站对走时偏移时间序列进行深加工处理，进而提取出适合于滇西北地区的地震短临异常识别指标(综合指标)。结果表明，该指标物理意义明确，具有对背景噪声波速测量结果高度综合的特点，可提高对异常的识别能力，对滇西北地区M≥5.0地震的发震时间具有较好的指示意义。但有几点需要注意：

1) 在对每个台站对瑞利波走时偏移时间序列进行效能评价时，本文选取±1.5倍标准差作为异常阈值，主要是为提高信噪比，阈值设置过低会带入部分噪声。

2) 在综合异常指标中，选取2倍均值线作为异常阈值，大部分指标在超过阈值指标线后较短时间内发生地震，因此在后续预测实践中，及时处理数据尤为重要。另外，该阈值可根据实际需要进行动态调整。

3) 本文研究是将背景噪声波速测量结果应用于地震短临预报的一次很好的尝试。研究发现，单个台站对的瑞利波走时偏移时间序列对整个云南地区M≥5.0地震的发震时刻均有较好的指示意义，因此本文方法可拓展至整个云南地区。基于云南地区固定台站的连续波形资料，提取适合于云南地区M≥5.0地震的短临异常识别指标，对整个云南地区的震情跟踪及预测预报工作都具有重要意义。

致谢: 本文研究所用连续波形资料来源于云南数字地震台网，GNSS基线数据来源于云南省地震局信息中心，背景噪声互相关计算程序由中国地震局地球物理研究所地球物理先导技术研究室王伟涛研究员提供，综合指标基于云南省保山市地震局李宗兴高级工程师开发的前兆处理软件计算得到，在此表示衷心感谢。