基于Welch 谱分析的全球MS ≥7.0 地震活动应变能周期特征

王利兵，郭 骄，郭建芳*，于春颂

（1.河北省地震局红山基准台, 邢台 055350；2.河北省地震局唐山地震监测中心站, 秦皇岛 066300）

0 引言

很多学者对历史强震活动的研究表明，全球8 级地震活动具有规律性[1]，主要地震带的MS≥7.0 地震活动存在彼此交替现象，低纬度各地震带之间有共同的活动周期[2]。近二十多年来中国学者对地震活动周期进行了广泛深入的研究，张郢珍等[3]利用最大熵谱法对全球巨大地震能量释放谱进行计算，得到10.53 a 的最显著周期；宋治平等[4]利用Morlet 小波分析方法提取全球MS≥7.0 地震活动周期谱，显示存在多个显著稳定的活动周期，50 a、80 a周期置信度达到95%，其次还有置信度相对较低的10 a 周期存在。这些周期现象与天体运行的各种周期有密切关系，从地球自转速率、太阳活动的周期分析中可以得到佐证[5-9]。在地震活动周期与天体运动关系的研究文献中，小波分析法比较常见。

提取周期的方法很多，不同方法得到的结果不尽相同。本文利用Welch 平均法周期图谱提取全球MS≥7.0 地震的应变能周期。为验证结果的合理性和可靠性，结合太阳黑子活动周期做对比，太阳黑子活动周期的计算采用傅里叶去多周期方法，同时参考相关研究成果进行深入讨论。

本文使用地震目录的时间范围是1700—2019年。全球7 级以上地震目录源自中国地震台网中心EQ11 目录（http://www.csndmc.ac.cn），并依据宋治平等[5]2011 年编著的《全球地震目录》对1700—1899年历史地震目录进行校正，1900 年之后的目录较为完整，根据相关研究成果对其进行修正，其最小完备震级均为7 级[4]，因此本文分析资料震级最小为7 级。另外，本文还利用太阳黑子数进行周期统计分析与检验，探讨其与全球大震活动周期的关系，数 据 来 源NASA （ National Aeronautics and Space Administration），研究时段为1900—2019 年。

1 地震活动应变能

地震释放的弹性应变能E是与构造应力、地震过程联系最为密切的物理量，反映地震活动的强弱，在研究地震活动周期中较为常用[10]。取某一地区一定时间段内地震序列能量平方根组成时间序列f(t)。

其中：N 为一定时间段内的地震次数，

利 用MAPSEIS 软 件对1700—2019 年全球MS≥ 7.0 地震目录进行时间扫描，月频度见图1，分别进行MS7.0～7.9、MS≥8.0 地震应变能E年值周期谱分析，时间步长、窗长均为1 年（图2）。结果显示：1900 年前、后7 级水平地震释放的能量存在阶变特征，后期明显高于前期，而8 级水平地震相差并不是很明显。该现象说明实际7 级水平地震频次远高于8 级地震，可见1900 年以前的目录限于观测技术，其完备性仍有欠缺。

图1 1700—2019 年全球7 级以上地震月频度Fig.1 The monthly frequency of global earthquakes with magnitude 7 or above from 1700 to 2019

2 功率谱估计方法及原理概述

功率谱估计的目的是根据有限数据给出信号、随机过程的频率成分分布，是频域内提取淹没在噪声中有用信息的分析方法[11]。全球强、大震活动具有一定规律，但又不能按确定的数学规律进行判断，可认为具有一定随机性。

Welch 功率谱法是基于周期法改进而来。在实际应用中，对得到的数据样本x(n)(0 ≤n≤N-1)切分为K个分段，令

其中：w(n)是一个长度为L的窗函数；D是偏移量。第i段的周期图定义为

3 周期计算

3.1 谱分析窗函数的选择

上百年尺度的地震能量信号有多个频率周期，频谱特征较为复杂，谱分析更多关注频率点而非能量的大小，故谱估计中窗函数的选择比较重要。矩形窗的主瓣较窄，分辨率好，方差较大，噪声水平较高；汉明窗和布莱克曼窗主瓣较宽，分辨率较低，方差较小，噪声水平较低；汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减小泄露观点出发，它与矩形窗相比，泄露、波动都减小了，并且选择性也得到提高，但其分辨率较低[11-12]。以1024 个FFT 频率点数绘制各窗函数幅频特征（图3a）。

图3 Welch 功率谱估计的各窗函数特征及算例Fig.3 The characteristics and examples of each window function using Welch power spectrum estimation

以MS7.0～7.9 地震为例，利用Welch 法分别选取矩形窗、汉宁窗、汉明窗、布莱克曼窗，256 个点数进行计算，以对数幅频谱形式产出，即S=lg(R2+I2)，R为 实部，I为虚部。如图3（b）所示，汉宁窗具有比较清晰的主瓣和旁瓣，第一旁瓣相对主瓣衰减明显，幅值最高，周期最全，出现了39 a、23 a、19 a、16 a、11 a 等大于10 a 周期。其他三个窗函数不论周期和幅值都基本一致，但主瓣和旁瓣十分不清晰，大于10 a 的周期只有19 a、11 a、10 a，说明估计地震应变能周期谱的能力相对偏弱，故采用汉宁窗比较合理。

3.2 谱估计结果

采用汉宁窗函数，谱计算用256 个点数，对MS7.0～7.9、MS≥8.0 地震应变能年值分别计算。本文主要分析地震活动周期T，选择前8 阶中T＞8 a的显著周期进行介绍，小于8 a 的周期谱值相对较小且繁多，不再详述。

1）计算MS7.0～7.9 地震时，经多次尝试选择合适点数，汉宁窗段点数取值范围在63～100 之间有较理想周期出现，小于63 时周期泄露偏多，大于100 时周期相同，故总计抽取19 个段点，分别为63、65、68、70、72、75、77、78、80、82、85、88、90、92、95、97、98、99、100，曲线代表不同段点同一周期的连线（图4a）。这些段点中，提取到的周期组数相对最完整的为8 组；最少的为5 组，共计27 个周期成分（说明：由于数据量偏大，8.1～8.9 a均归为8 a 期，其他依此类推），由小到大为8 a、9 a、10 a、11 a、12 a、13 a、14 a、15 a、16 a、17 a、18 a、19 a、20 a、21 a、22 a、23 a、24 a、26 a、28 a、30 a、32 a、36 a、39 a、46 a、51 a、56 a、64 a。

图4 MS7.0～7.9 地震、MS≥8.0 地震周期及其贡献率、频数变化Fig.4 The period of earthquakes with MS 7.0～7.9 and earthquake with MS≥8.0, their contribution rates and frequency changes

假设8 个周期在不同段点中均应存在，不同段点未提取到的周期以零补充，此时以不同周期在各个段点中出现的频数排列成19 行27 列的矩阵，用数据降维方法[12]提取主要成分及其贡献率，该方法首先对原始数据进行标准化转换、计算样本相关系数矩阵，其次计算相关系数矩阵的特征值和特征向量，最后计算主成分贡献率，选择重要的主成分，主成分贡献率越大，说明它包含的原始信息量越大。主成分分析用较少的变量代替了原来较多的变量。

结果显示，27 个成分中前8 个相对突出，其贡献率依次为16.80 %、13.83 %、12.72 %、10.41 %、9.77 %、8.58 %、5.38 %、5.26 %，总贡献率82.75 %，较为显著。依据图4a 各周期出现的总频次高低，前8 个周期分别为8 a（14 次）、10 a（11 次）、51 a（7 次）、11 a（7 次）、17 a（7 次）、12 a（6 次）、19 a（6 次）、39 a（6 次），总占比约64.41 %，故推测该8 个周期为主要周期成分（图4b）。需要说明的是图4b、4d 贡献率代表27 个贡献率由大到小分布，频数代表27 个周期由小到大分布）。

2）计算MS≥8.0 地震时思路同上，经挑选，段点数设置范围在65～100 之间，共20 个段点数（分别为65、66、67、69、71、77、78、80、82、83、84、85、87、89、90、91、96、98、99、100），提取到的周期组数相对最完整的为5 组；最少的为3 组，共计17 个周期成分，由小到大为8 a、10 a、12 a、13 a、14 a、18 a、19 a、20 a、21 a、22 a、26 a、28 a、30 a、42 a、46 a、51 a、56 a（图4c），可见8 级水平地震的出现周期比7 级水平地震少很多。

以同样方法提取主要成分及其贡献率，有5 个成 分 较 显 著，贡 献 率 依 次 为：29.35 %、20.72 %、17.33 %、9.00 %、5.39 %，总贡献率81.79 %。依据图4c 各周期出现的总频次，相对突出的分别为10 a（16 次）、46 a（8 次）、22 a（7 次）、8 a（5 次）、20 a（5 次），总占比约64.56 %，故推测该5 个周期为主要周期成分。

3.3 太阳黑子活动周期

太阳黑子是太阳光球表面经常出现的阴暗斑点，是太阳活动的主体，黑子数量的多少表征太阳活动的强弱，太阳风暴的成因往往归因于太阳黑子的影响[13]。前人认为地球电磁学与太阳黑子等宇宙天体也存在密切联系[14-15]，在天文研究的基础上，地震活动规律与太阳黑子变化周期之间的关系得到印证[8,16-17]。屠泓为等[7]利用FFT 法计算了1700—2008 年太阳黑子活动周期，并与中国大陆西部（E108°以西）MS≥5.0 地震活动规律做对比，认为不同震级的地震活动周期均受一定程度的太阳黑子活动影响。

在前人研究基础上，本文采用傅里叶去多个周期法提取1749—2019 年太阳黑子活动周期。该方法首先将时序数据去倾和中心化，然后根据一定条件确定周期成分的振幅和周期，从中优选出最大振幅和其对应的周期，利用傅里叶滑动法排除原始数据中的周期成分，得到残差。重复计算直到残差序列中没有满足条件的周期成分，具体计算方法见文献 [18]。

图5 为该方法计算的周期数据和残差结果，图6为提取到的共计30 个周期及其振幅，按振幅大小，表1 计算了前10 个显著周期成分的显著性检验，其中6 个为特别显著周期（显著性大于90 %）分别为10.8 a、10.0 a、11.8 a、10.4 a、90.3 a、11.3 a；4 个为显著周期（显著性大于70 %），分别为67.8 a、8.5 a、9.3 a、54.2 a。综合看，显著性最高的周期约为10～11 a（平均10.9 a），其次约为67 a、8～9 a、54 a。

表1 太阳黑子10 个主要周期成分及显著性检验Table 1 The ten main periodic components of sunspots and their significance test

图5 太阳黑子傅里叶去多个周期分析Fig.5 The analysis of sunspots with Fourier eliminating multiple period method

图6 太阳黑子傅里叶周期分析提取到的周期及其振幅Fig.6 The periods and amplitudes of sunspots extracted from Fourier cycle analysis

3.4 太阳黑子活动周期与地震周期的关系

根据提取到的太阳黑子活动显著周期和地震活动显著周期，做二者线性相关分析及显著性检验。由于3.2、3.3 节使用了不同数学方法提取不同对象的周期特征，其结果存在差异，计算线性相关时，太阳黑子显著周期数量与MS7.0～7.9、MS≥8.0 地震显著周期数量并不一致。为此，以太阳黑子周期数量为参考（10 个），地震周期数量不足部分从主要成分计算结果中弥补，贡献率大的优先选取（表2），表中周期的排列顺序由小到大，只是为了计算同类时间尺度周期的相关性，不能反映各周期的显著性对应顺序。

表2 显著性活动周期统计表Table 2 The statistical table of significant activity cycle

回归计算结果见图7 及表3。太阳黑子活动周期分别与MS7.0～7.9、MS≥8.0 地震的线性相关度达到0.9873、0.9704，二者均显示高度显著。此外，还显示了各周期在8～22 a（前7 个周期）左右的另外一组线性关系，其相关系数为0.9554、0.9784，同样具有高度相关性，但在7 级水平地震中，此相关系数略低于10 个周期的系数约0.03，而在8 级以上地震中两个相关系数基本一致。同时该组周期斜率大、拟合点集中，反映了各参数存在的20 年以下周期，尤其是8 a、10～11 a 左右周期很显著，其中太阳黑子活动显著于地震活动。

表3 太阳黑子活动周期与地震活动周期线性回归分析结果统计Table 3 The statistical results of linear regression analysis between sunspot activity period and seismic activity period

图7 太阳黑子活动周期与MS7.0～7.9（a）、MS≥8.0（b）地震活动周期线性回归相关计算Fig.7 The linear regression correlation between sunspot activity period and seismic activity period of MS7.0～7.9 (a) and MS≥8.0 (b)

4 结论与讨论

本文得到的全球MS≥7.0 地震的显著周期有以下几个：MS7.0～7.9 地震约为8 a、10 a、51 a、11 a、17 a、12 a、19 a、39 a；MS≥8.0 地震约为10 a、46 a、22 a、8 a、20 a。利用傅里叶去多周期方法提取1749—2019 年太阳黑子周期，出现6 个特别显著周期成分，分别为10.8 a、10.0 a、11.8 a、10.4 a、90.3 a、11.3 a；4 个显著周期成分分别为67.8 a、8.5 a、9.3 a、54.2 a。综合分析，显著性最高的周期约为10～11 a（平均10.9 a），其次约为90 a、67 a、8～9 a、54 a。

10～11 a 周期在MS7.0～7.9 地震、MS≥8.0 地震和太阳黑子活动周期中均为优势周期，与前人使用最大熵谱法、Morlet 小波分解法获得的最显著周期是一致的[3，19]，说明太阳黑子11 周年变化与地球内部活动的关系是非常密切的，比如华北地区6 级以上强震活动有11 年左右的起伏特征[16]、华北北部地下水位存在10 年准周期变化，高值段是5 级以上地震从的出现时间 [20]等。

8 a 周期在MS7.0～7.9 地震周期中出现概率偏高，在MS≥8.0 地震周期和太阳黑子活动周期中虽然存在，但不显著。该周期在20 世纪环太平洋强震的周期研究中曾呈现最高峰值[3]，目前仍是显著周期，说明是天体和板块活动联合作用的结果，因此比10～11 a 太阳活动周期短。

稍长约50 a 左右的周期在MS7.0～7.9 地震、MS≥8.0 地震和太阳黑子活动中出现，但有几年偏差，分别是51 a、46 a、54 a。尹继尧等[21]利用Morlet小波分析法提取到：1700—2010 年MS≥7.0 地震置信度65% 的47.5 a 周期，MS≥7.5 地震置信度95%的48.6 a 周 期，MS≥8.0 地 震 置 信 度95 % 的50.6 a周期，MS≥8.3 地震置信度95% 的51.7 a 周期。本文结果与此大致相同，大致都在50 a 左右。

21～22 a 周期反映了太阳活动大周期事件，它可能引发全球地震年统计约150 的波动，且引起地球自转速率2 ms 的起伏[8]。但在全球各地震带中，该活动周期并不十分显著，只有少部分地震带有显著活动，如宋治平等利用Morlet 小波分析法提取到1900—2010 年中国大陆置信度95 %的23.8 a 周期，阿拉斯加湾—瓜达卢佩岛带置信度50 % 的23.8 a周期，开曼群岛—委内瑞拉带置信度70 %的23.8 a周期，1900—2010 年阿拉斯加湾—瓜达卢佩岛带置信度90 %的22.7 a 周期，开曼群岛—委内瑞拉带置信度50 %的23.8 a 周期，其他地震带类似周期置信度偏低[4]。

本文提取到的太阳黑子90 a、67 a 长周期，未在MS≥8.0 地震周期中显示，但在MS7.0～7.9 震级档段点数取99 时提取到的最长周期为64 a，该周期（64～67 a 左右）在以上提到的文献中也较少出现，可见其可信度低。约90 a 周期可能是一个弱周期，Welch 方法未提取到与其接近的周期变化，尹继尧等利用Morlet 小波分析法对全球巨大地震活动周期的分析认为，50 a 左右周期最显著，但还存在80～100 a 左右显著性不强的活动周期，其中对1700—2010 年全球MS≥8.0 地震周期分析中，87.8 a周期置信度为50%，而1900—2010 年的88.8 a 周期置信度仅为20%[21]。李文龙等利用Morlet 小波分析法提取到了1700—2017 年太阳黑子8～11 a、96～102 a 的周期[22]；白春华等利用小波变换方法提取到了1590—1990 年 地 磁 主 磁 场 存 在 稳 定 的30 a 周期、准50 a 周期和强度变化的110 a 周期，其中50 a周期成分主要由偶极子场的轴向分量贡献[23]，这与全球强震50 a 周期应该有一定相关性。

由此看出，地磁场与太阳黑子活动的相关性高于地震活动，说明了地震活动周期的复杂性。另外，长周期变化的结构非常复杂，周期分量的强度和大小会随研究时间尺度变化，这也是不同分析方法的研究结果容易出现偏差或遗漏的原因。