奇异谱分析在地磁数据研究中的应用

于畅，路来君，于洪池，吕铁鑫，魏美璇，陈卓

1.吉林大学 地球科学学院，长春 130061； 2.吉林省地震局 吉林地震台，长春 130117； 3.吉林省地震局 合隆地震台，长春 130200； 4.吉视传媒股份有限公司 四平分公司，吉林 四平 136000； 5.吉林省地震局 信息中心，长春 130117

0 引言

地磁又称“地球磁场”，是指地球周围空间分布的磁场。随着地球探测技术的发展，人们逐步展开对地球磁场的深入认识和本质探讨。在地震研究中，地磁学科的基本任务是基于地磁学的理论和方法技术，探索与孕震过程有关的地磁变化规律，为地震预报及孕震过程的理论研究服务。地磁场共有7个地磁要素，分别为F、H、D、X、Y、Z、I。其中，F为地磁场总强度，H为地磁场水平分量，D为磁偏角(磁北与地理北的夹角)，X为地磁场北向分量，Y为地磁场东向分量，Z为地磁场垂直分量，I为磁倾角(磁子午面内F与H的夹角)[1]。

在中国用于地磁数据分析的方法有多种[2--5]，不同分析方法使用的地磁分量和计算原理不同。中国学者多运用地磁数据分析方法对距震中一定范围内多个台站的地磁数据分量进行处理分析，研究震前、震中、震后地磁数据的异常变化，从地磁异常持续时间、变化幅度和异常分布范围等方面预报地震。国外学者在地磁数据分析中也取得大量成果。日本学者Hayakawa et al.[6]针对1993年8月8日关岛7.1级地震研究超低频磁噪声的变化特性，发现地磁垂直分量Z与水平分量H的比值在描述空间等离子体波地震源发射方面具有重要意义。Hayakawa et al.[7]认为短期地震预报对日本等地震活跃的国家来说是紧迫的重要课题，而实验表明与地震相关的地磁异常现象确有存在。

在国内外的地磁学研究中，应用地磁异常进行地震预报的准确率不高，因此需不断探索新的地磁数据分析方法，以提高运用地磁异常进行地震预报的可靠性。奇异谱分析能够强化原始数据的真实特征，目前在地磁数据中应用奇异谱分析的相关研究较少，因而笔者以长春地磁台采集的地磁数据为研究对象，采用奇异谱分析方法对时间序列进行转换，对比原始序列与重构序列的趋势及异常变化，通过对重构序列展开数学研究达到对原始序列的趋势预测，实现地磁异常提取，为地磁异常与地震预报的准确对应提供新途径。

1 奇异谱分析原理

奇异谱分析(Singular Spectrum Analysis, SSA)是一种广义功率谱分析。据资料记载[8]，奇异谱分析理论最早于1978年被提出并应用于海洋学研究中，它是针对一维时间序列的主成分分析方法[9]，具有稳定的识别和强化信号功能，优点在于能够过滤噪声并提取蕴含在时间序列中的不规则波动和随机性特征[10]。

奇异谱分析的核心思想是将一维时间序列转换为多维空间序列的过程[11]，通过对转换的轨迹矩阵构造、分解、重构进而提取出信号的变化趋势、周期和噪声等成分。现将主要算法叙述如下：

确定一维时间序列X=(x1，x2，… ，xN)，其中N为序列长度。构造轨迹矩阵D，式中L为窗口长度。

(1)

根据轨迹矩阵D求协方差矩阵C：

C=DDT

(2)

对协方差矩阵C进行特征值分解，求出特征值λ和特征向量v。则第k个特征值对应的主成分表示为：

(3)

(4)

计算主成分对原始序列的贡献率，若前r个主成分可表现出原始序列的全部特征，则重构序列的表达式如下，原始序列的余项为噪声信号。

(5)

2 数据预处理

笔者选用吉林省长春地磁台2017—2019年FHDZ--M15相对记录仪采集的Z、H、D、F分量原始分钟值数据进行奇异谱分析。Z、H、D、F分量的原始数据曲线如图1所示。

图1 各分量原始分钟值Fig.1 Original minute value for each component

采用每日一值的思想提取世界时16～20时Z、H、D、F原始分钟平均值(夜均值)替代当日各分量数据。取此预处理方法是由于夜间地磁场受空间磁场影响和环境干扰相对较小，数据变化较其他时段更稳定，表现出的地磁特征更准确，同时亦可避免计算量冗长[13]。预处理后Z、H、D、F各分量的每日一值曲线见图2。

图2 各分量预处理分钟值Fig.2 Preprocessed minute value for each component

3 奇异谱分析在地磁数据中的应用

3.1 基于SSA的数据拟合

在奇异谱分析中窗口长度的选取对信号提取和信号分析具有较大影响，一般取周期的整数倍且满足N/3≤L≤N/2[14]。通过多次试验，选定窗口长度L=320，以便较好地提取时间序列的信号成分。

对Z、H、D、F各分量创建轨迹矩阵并求其主成分在总特征值中的占比(表1)得出，各分量的前10个特征值的贡献率之和分别为99.996 0%、99.882 5%、99.832 4%、99.999 9%，这表明运用前10个特征值能够充分表现原始序列的信号特征，其后的特征值贡献率较弱，因而采用前10个特征值对各分量进行奇异谱分析。拟合数据与预处理数据对比见图3。

表1 各分量前10个特征值贡献率

由图3可知，拟合后的曲线较平滑且能够准确表现出原始数据的变化趋势。其中Z分量、F分量呈逐年上升趋势，H分量相对较缓、年变化趋势略有下降，D分量呈逐年下降趋势。

图3 各分量夜均值与SSA拟合值Fig.3 Night mean and SSA fitting values for each component

3.2 拟合精度检验

为检验奇异谱分析方法的拟合精度，选取一元线性回归分析(Unary Linear Regression，ULR)作为对比方法，对各分量预处理数据再次进行数据拟合。根据各分量夜均值与时间刻度建立一元线性回归分析表达式为：

y=β0+β1x+ε

(6)

式中：x为时间刻度，y为夜均值拟合值，ε为误差。经计算不同分量的各项系数见表2。由此得出Z、H、D、F各分量夜均值、SSA拟合值和ULR拟合值的对比曲线见图4。

表2 一元线性回归分析系数

图4 各分量夜均值、SSA拟合值和ULR拟合值对比Fig.4 Night mean, SSA fitting values and ULR fitting values for each component

从对比曲线可以看出，奇异谱分析和一元线性回归分析都较好地体现出数据的变化趋势，但相比之下奇异谱分析能更精准地表达出数据变化的细节。Z、H、D、F各分量的残差对比曲线见图5。计算各分量残差值距0值线的偏离程度，用以比较两种方法的拟合效果，所求结果见表3。由此可见4个分量的奇异谱分析残差量均小于一元线性回归分析，故奇异谱分析的拟合效果更显著。

表3 奇异谱分析和一元线性回归分析的残差偏离值

图5 各分量SSA和ULR残差值对比Fig.5 Comparision of residual error values of SSA and ULR for each component

经奇异谱分析的拟合曲线与原始分量观测值曲线相比更加连续稳定，可作为趋势变化的背景值。相比之下原始观测数据则含有真实物理量变化、干扰噪声和异常信号等信息。如果用原始观测数据减掉背景数据所获得的时间序列，则突出了背景噪声和异常信号，这有助于分析背景噪声并放大短临异常信号。在本论文中原始数据与背景数据的差值序列为奇异谱分析的残差序列，这突出了残差序列在地磁异常研究中的重要性。

运用RMSE(均方根误差)、MAE(平均绝对误差)、MSPE(均方百分比误差)3项误差评价指标对奇异谱分析和一元线性回归的拟合结果进行归算[15]，数值越小表明拟合方法的精度越高。各方法的数学表达式为：

(7)

(8)

(9)

公式7～9中：m表示序列长度，h表示奇异谱分析的拟合值，y表示一元线性回归分析的拟合值，所求结果见表4。从精度评价结果中可见奇异谱分析的数据拟合精度优于一元线性回归分析。

表4 各分量奇异谱分析与一元线性回归拟合精度

3.3 数据预测

选用ARMA模型(自回归滑动平均模型)分别对4个分量的夜均值序列进行数值预测。经检验各分量夜均值一阶差分序列为平稳时间序列[16--17]。运用自相关及偏自相关系数初步定阶，残差检验结果显示其不存在一阶相关性，且各残差序列接近正态分布，由此证明上述流程符合ARMA建模要求。由各分量的夜均值一阶差分预测值归算至夜均值预测值，并给出95%置信区间。Z、H、D、F各分量的原始数据夜均值后30 d的预测结果见图6。

图6 各分量ARMA预测值Fig.6 Prediction value of ARMA for each component

4 结论

(1)与一元线性回归分析相比，奇异谱分析后的拟合误差远小于一元线性回归分析，其残差曲线更接近0值线，这表明奇异谱分析的拟合精度更高，拟合效果更加显著。

(2)通过对奇异谱分析后的序列进行ARMA模型预测得出，预测值符合原始序列的趋势变化，且预测值未超出置信区间范围，表明短期内将无地磁异常发生。