给定嵌入长度L,即窗口长度,构建给定时间序列X_i=(x_i,…,x_i+L-1)^T[KG*2](1≤i≤k)的轨迹矩阵X_L×K,其中K=N-L+1,1<L≤N/2,计算方法见下式:

X=[X₁:…:X_k]=(x_ij)^L,K_i,j=1=[JB([][HL(4:1,*2/3; 2,*2/3; 3,*2/3; 4,*2/3]x₁[]x₂[]…[]x_N-L+1

x₂[]x₃[]…[]x_N-L+2

[][][]

x_L<[]x_L+1[]…[]x_N[HL)][JB)]]

SSA实现的过程需要确定2个重要参数——嵌入维度L和重构阶次P。嵌入维度L的选择直接影响SSA分解信号的效果,研究表明:若L较小,将会导致信号分解不彻底; 反之,将会使低频信号的分解产生谱裂现象(Vautard et al,1992)。对于待提取年变信号的观测数据,L的选取一般为年变周期的整数倍,这样分解效果较好。经过试验,本文设置L为整周年天数365 d。原始曲线如果存在显著的趋势,先利用高阶小波去除长周期趋势项,根据奇异谱曲线动态选取重构阶次P,确保重构序列位于年频段范围。

设矩阵[WTHX]B=XX^T,其中[WTHX]X^T为轨迹矩阵[WTHX]X的转置。接下来计算矩阵B的所有特征值λ_i及其特征向量U_i,按照降序将特征值排列为λ₁,λ₂,…,λ_L(λ₁≥λ₂≥…≥λ_L≥0),其对应的特征向量分别为[WTHX]U₁,U₂,…,U[WTBX]_L。

_i=[WTHX]X^T[WTHX]U_i/[KG-*3][KF(]λ_i[KF)](i=1,2,…,d)(2)

式中:[WTHX]U[WTBX]_i和[WTHX]V[WTBX]_i是轨迹矩阵[WTHX]X的左右特征向量; d=L^*; L^*=min{L,K}; 轨迹矩阵[WTHX]X可以由初等矩阵按照下式合成:

X=X₁+X₂+…+X_[WT]d(3)

式中:初等矩阵[WTHX]X[WTBX]_i=λ_i[KF)][WTHX]U[WTBX]_i V^T_i, rank(X_i)=1, 而[KG-*3] λ₁≥[KG-*3][KF(]λ₂[KF)]≥…≥[KG-*3][KF(]λ_L≥0为时间序列[WTHX]X的奇异谱(Broomhead,King,1986; Shen,Li,2012)。由最大特征值及其对应的特征向量计算得到的初等矩阵,对应于原始时间序列的最大变化方向,反之则视为噪声。

将初等矩阵[WTHX]X_i的下标(i=1,2,…,d)分成p个不相交的子集I₁,I₂,…,I_p, 设I={i₁,i₂,…,i_m}, 则:

X^*=X_I1+X_I2+…+X_Ip=(x^*_ij)^L,K_i,j=1(4)

选取集合I₁,I₂,…,I_p的过程即为分组。研究表明,奇异值接近的2个初等矩阵重构得到一个具有优势周期的成分。通过求解并观测奇异值的分布情况,对观测时间序列分组重构,进而通过快速傅里叶变化(FFT)对重构的年变序列进行周期检验,确保其优势周期位于年频段范围内。

通过上述计算得到的重构矩阵[WTHX]X^*,需要进行对角平均计算(Zhigljavsky,2010),从而得到一个重构的时间序列y₁,y₂,…,y_N,展开后的计算方法为:y_k=[JB({][SX(]1[]k[SX)]∑[DD(]k[]m=1[DD)]x^*_m,k-m+1 1≤k≤L^*

[SX(]1[]L^*[SX)]∑[DD(]L^*[]m=1[DD)]x^*_m,k-m+1 L^*≤k≤K^*

[SX(]1[]N-K+1[SX)] ∑[DD(]N-K^*+1[]m=K-K^*+1[DD)]x^*_m,k-m+1[KG*1]K^*≤k≤N[JB)]

式中:L^*=min{L,K}, K^*=max{L,K}。

图3原始观测时间序列显示,2个测项在年变的基础上叠加了显著的趋势项变化,其中东风煤矿钻孔倾斜EW存在单一的西倾趋势(图3a),而巴里坤水平摆倾斜NS存在多个线性变化趋势(图3b)。因此,首先利用db5小波进行高阶滤波,去除长周期的趋势变化,进而对变化相对平稳的去趋势时间序列做奇异谱分析。

去趋势时间序列的奇异谱曲线显示(图4a),2个测项的第1,2个奇异值位于噪声平台的上方,占比较高且存在显著的周期台阶,因此取重构阶次P为[1,2],重构结果见图4b,c中红色和蓝色虚线部分,灰色曲线为去趋势后的观测时间序列。从图4中可以看出,东风煤矿钻孔倾斜EW分量去掉趋势项后,部分年份如2013,2016年,明显偏离SSA重构得到的年周期变化规律,出现破年变信号; 巴里坤水平摆NS的多个长周期趋势

图3 原始观测时间序列
Fig.3 Original observation time series

图4 奇异谱曲线(a)及SSA年变重构结果(b)(c)
Fig.4 Singular spectral analysis curve(a)and results of annual variation reconstruction based on the SSA algorithm(b)(c)

项,使得趋势转折与破年变变化混在一起难以区分,通过去趋势并结合SSA重构的结果来看,2017,2018以及2019年均存在显著的破年变异常,与笔者从原始曲线观测到的破年变异常年份基本一致,只是幅度略有差异。

对SSA算法重构得到的年变时间序列做快速傅里叶变换(FFT)的结果显示(图5),2个测项通过SSA重构后的时间序列,其优势周期均为372 d,这表明重构时间序列成功分离出了原始曲线中年频段的成分。将去趋势的观测序列与重构时间

图5 重构年变成分的FFT变换结果及震例空间分布
Fig.5 The FFT transformation results of annual variation reconstruction and the spatial distribution of earthquake cases

序列做差,得到了去趋势和年变(SSA重构的正常年变)后的残差时间序列,在此基础上进行动态R值检验。

取台站周边200 km范围内5级以上、300 km范围内6级以上地震(Dobrovolsky et al,1979),作为输入震例进行R值检验,检索到的震例及台站的空间分布情况如图6所示,2个蓝色圆圈分别代表了以台站为中心的200,300 km的圆形震例检索区。

动态R值检验过程对应于2个互相嵌套的循环,需要设置2个循环变量:

(1)异常判定阈值

通常取均值加减几倍的均方差作为异常判定的标准,但不同观测方法、不同台站测项的信噪比差别较大,取单一的阈值会漏掉部分异常信息,或引入部分噪声信息。因此本文将上述异常判定阈值扩大到某一个动态范围,循环进行R值检验,根据检验结果确定最佳阈值。本例中设定阈值范围为1～3倍的均方差,循环步长设为0.1倍的均方差。

(2)预测时间

预测时间即异常从开始到结束所持续的时间,通过循[JP2]环设定该变量并进行R值检验,可以得到某一观测具有最佳预测效能的预测时间取值,本文设定该变量的取值范围0～720 d,循环步长为10 d。

循环计算后得到了“R值-异常判定阈值-预

图6 震例空间分布(此2例中均无漏报地震)
Fig.6 Spatial distribution of earthquake cases(There was no under reported earthquake in both cases)

测时间”的空间插值图像(图7a,8a),其色标代表了不同的异常判定阈值-预测时间组合所对应的预测效能检验结果,即R值评分,取值范围-1～1。从中可以看到不同组合的明显差异性,图中黑色圆点代表了R值最高点,将对应的异常判定阈值及预测时间代入残差时间序列,提取异常特征时段进行R值检验,得到了具有最佳映震效能的异常识别及震例回溯结果如图7b,8b所示。

图7 东风煤矿钻孔倾斜EW动态R值扫描(a)与异常自动识别结果
Fig.7 The dynamic R-value scanning(a)and automatic anomaly recognition results of EW component of borehole tilt(b)in Dongfeng coal mine station

图8 巴里坤水平摆NS动态R值扫描(a)与异常自动识别结果(b)
Fig.8 The dynamic R-value scanning(a)and automatic anomaly recognition results of NS component of horizonal pendulum tilt(b)in Balikun station

图7b,8b中蓝色虚线为具有最佳R值评分的异常判定阈值线,黑色实线为数据均值线。从动态检验的结果来看:①上述2个测项具有最佳预测效能的异常判定阈值分别为均值加减1.6倍、2倍的均方差; ②异常出现后最佳预测时间分别为91 d,61 d,目前2个测项均处在异常阶段; ③东风煤矿钻孔倾斜EW的破年变异常对其周边300 km内的5级以上地震具有较好的映震效能,震例主要集中在台站的ES方向; ④巴里坤水平摆倾斜NS破年异常主要对应200 km内的5级以上地震,震例较少且无优势集中方向; ⑤应用本研究中的方法提取的破年变异常,东风煤矿钻孔倾斜EW破年变异常的预测效能大于97.5%置信度的临界值R₀,预测结果较为可靠; 而巴里坤水平摆倾斜NS虚报率较高,小于临界值R₀,虽然目前处于异常时段,但信度较低,建议参考使用。