于暢，路來君，于洪池，呂鐵鑫，魏美璇，陳卓

1.吉林大學 地球科學學院，長春 130061； 2.吉林省地震局 吉林地震臺，長春 130117； 3.吉林省地震局 合隆地震臺，長春 130200； 4.吉視傳媒股份有限公司 四平分公司，吉林 四平 136000； 5.吉林省地震局 信息中心，長春 130117

0 引言

地磁又稱“地球磁場”，是指地球周圍空間分布的磁場。隨著地球探測技術的發(fā)展，人們逐步展開對地球磁場的深入認識和本質(zhì)探討。在地震研究中，地磁學科的基本任務是基于地磁學的理論和方法技術，探索與孕震過程有關的地磁變化規(guī)律，為地震預報及孕震過程的理論研究服務。地磁場共有7個地磁要素，分別為F、H、D、X、Y、Z、I。其中，F(xiàn)為地磁場總強度，H為地磁場水平分量，D為磁偏角(磁北與地理北的夾角)，X為地磁場北向分量，Y為地磁場東向分量，Z為地磁場垂直分量，I為磁傾角(磁子午面內(nèi)F與H的夾角)[1]。

在中國用于地磁數(shù)據(jù)分析的方法有多種[2--5]，不同分析方法使用的地磁分量和計算原理不同。中國學者多運用地磁數(shù)據(jù)分析方法對距震中一定范圍內(nèi)多個臺站的地磁數(shù)據(jù)分量進行處理分析，研究震前、震中、震后地磁數(shù)據(jù)的異常變化，從地磁異常持續(xù)時間、變化幅度和異常分布范圍等方面預報地震。國外學者在地磁數(shù)據(jù)分析中也取得大量成果。日本學者Hayakawa et al.[6]針對1993年8月8日關島7.1級地震研究超低頻磁噪聲的變化特性，發(fā)現(xiàn)地磁垂直分量Z與水平分量H的比值在描述空間等離子體波地震源發(fā)射方面具有重要意義。Hayakawa et al.[7]認為短期地震預報對日本等地震活躍的國家來說是緊迫的重要課題，而實驗表明與地震相關的地磁異?，F(xiàn)象確有存在。

在國內(nèi)外的地磁學研究中，應用地磁異常進行地震預報的準確率不高，因此需不斷探索新的地磁數(shù)據(jù)分析方法，以提高運用地磁異常進行地震預報的可靠性。奇異譜分析能夠強化原始數(shù)據(jù)的真實特征，目前在地磁數(shù)據(jù)中應用奇異譜分析的相關研究較少，因而筆者以長春地磁臺采集的地磁數(shù)據(jù)為研究對象，采用奇異譜分析方法對時間序列進行轉換，對比原始序列與重構序列的趨勢及異常變化，通過對重構序列展開數(shù)學研究達到對原始序列的趨勢預測，實現(xiàn)地磁異常提取，為地磁異常與地震預報的準確對應提供新途徑。

1 奇異譜分析原理

奇異譜分析(Singular Spectrum Analysis, SSA)是一種廣義功率譜分析。據(jù)資料記載[8]，奇異譜分析理論最早于1978年被提出并應用于海洋學研究中，它是針對一維時間序列的主成分分析方法[9]，具有穩(wěn)定的識別和強化信號功能，優(yōu)點在于能夠過濾噪聲并提取蘊含在時間序列中的不規(guī)則波動和隨機性特征[10]。

奇異譜分析的核心思想是將一維時間序列轉換為多維空間序列的過程[11]，通過對轉換的軌跡矩陣構造、分解、重構進而提取出信號的變化趨勢、周期和噪聲等成分。現(xiàn)將主要算法敘述如下：

確定一維時間序列X=(x1，x2，… ，xN)，其中N為序列長度。構造軌跡矩陣D，式中L為窗口長度。

(1)

根據(jù)軌跡矩陣D求協(xié)方差矩陣C：

C=DDT

(2)

對協(xié)方差矩陣C進行特征值分解，求出特征值λ和特征向量v。則第k個特征值對應的主成分表示為：

(3)

(4)

計算主成分對原始序列的貢獻率，若前r個主成分可表現(xiàn)出原始序列的全部特征，則重構序列的表達式如下，原始序列的余項為噪聲信號。

(5)

2 數(shù)據(jù)預處理

筆者選用吉林省長春地磁臺2017—2019年FHDZ--M15相對記錄儀采集的Z、H、D、F分量原始分鐘值數(shù)據(jù)進行奇異譜分析。Z、H、D、F分量的原始數(shù)據(jù)曲線如圖1所示。

圖1 各分量原始分鐘值Fig.1 Original minute value for each component

采用每日一值的思想提取世界時16～20時Z、H、D、F原始分鐘平均值(夜均值)替代當日各分量數(shù)據(jù)。取此預處理方法是由于夜間地磁場受空間磁場影響和環(huán)境干擾相對較小，數(shù)據(jù)變化較其他時段更穩(wěn)定，表現(xiàn)出的地磁特征更準確，同時亦可避免計算量冗長[13]。預處理后Z、H、D、F各分量的每日一值曲線見圖2。

圖2 各分量預處理分鐘值Fig.2 Preprocessed minute value for each component

3 奇異譜分析在地磁數(shù)據(jù)中的應用

3.1 基于SSA的數(shù)據(jù)擬合

在奇異譜分析中窗口長度的選取對信號提取和信號分析具有較大影響，一般取周期的整數(shù)倍且滿足N/3≤L≤N/2[14]。通過多次試驗，選定窗口長度L=320，以便較好地提取時間序列的信號成分。

對Z、H、D、F各分量創(chuàng)建軌跡矩陣并求其主成分在總特征值中的占比(表1)得出，各分量的前10個特征值的貢獻率之和分別為99.996 0%、99.882 5%、99.832 4%、99.999 9%，這表明運用前10個特征值能夠充分表現(xiàn)原始序列的信號特征，其后的特征值貢獻率較弱，因而采用前10個特征值對各分量進行奇異譜分析。擬合數(shù)據(jù)與預處理數(shù)據(jù)對比見圖3。

表1 各分量前10個特征值貢獻率

由圖3可知，擬合后的曲線較平滑且能夠準確表現(xiàn)出原始數(shù)據(jù)的變化趨勢。其中Z分量、F分量呈逐年上升趨勢，H分量相對較緩、年變化趨勢略有下降，D分量呈逐年下降趨勢。

圖3 各分量夜均值與SSA擬合值Fig.3 Night mean and SSA fitting values for each component

3.2 擬合精度檢驗

為檢驗奇異譜分析方法的擬合精度，選取一元線性回歸分析(Unary Linear Regression，ULR)作為對比方法，對各分量預處理數(shù)據(jù)再次進行數(shù)據(jù)擬合。根據(jù)各分量夜均值與時間刻度建立一元線性回歸分析表達式為：

y=β0+β1x+ε

(6)

式中：x為時間刻度，y為夜均值擬合值，ε為誤差。經(jīng)計算不同分量的各項系數(shù)見表2。由此得出Z、H、D、F各分量夜均值、SSA擬合值和ULR擬合值的對比曲線見圖4。

表2 一元線性回歸分析系數(shù)

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

從對比曲線可以看出，奇異譜分析和一元線性回歸分析都較好地體現(xiàn)出數(shù)據(jù)的變化趨勢，但相比之下奇異譜分析能更精準地表達出數(shù)據(jù)變化的細節(jié)。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

經(jīng)奇異譜分析的擬合曲線與原始分量觀測值曲線相比更加連續(xù)穩(wěn)定，可作為趨勢變化的背景值。相比之下原始觀測數(shù)據(jù)則含有真實物理量變化、干擾噪聲和異常信號等信息。如果用原始觀測數(shù)據(jù)減掉背景數(shù)據(jù)所獲得的時間序列，則突出了背景噪聲和異常信號，這有助于分析背景噪聲并放大短臨異常信號。在本論文中原始數(shù)據(jù)與背景數(shù)據(jù)的差值序列為奇異譜分析的殘差序列，這突出了殘差序列在地磁異常研究中的重要性。

運用RMSE(均方根誤差)、MAE(平均絕對誤差)、MSPE(均方百分比誤差)3項誤差評價指標對奇異譜分析和一元線性回歸的擬合結果進行歸算[15]，數(shù)值越小表明擬合方法的精度越高。各方法的數(shù)學表達式為：

(7)

(8)

(9)

公式7～9中：m表示序列長度，h表示奇異譜分析的擬合值，y表示一元線性回歸分析的擬合值，所求結果見表4。從精度評價結果中可見奇異譜分析的數(shù)據(jù)擬合精度優(yōu)于一元線性回歸分析。

表4 各分量奇異譜分析與一元線性回歸擬合精度

3.3 數(shù)據(jù)預測

選用ARMA模型(自回歸滑動平均模型)分別對4個分量的夜均值序列進行數(shù)值預測。經(jīng)檢驗各分量夜均值一階差分序列為平穩(wěn)時間序列[16--17]。運用自相關及偏自相關系數(shù)初步定階，殘差檢驗結果顯示其不存在一階相關性，且各殘差序列接近正態(tài)分布，由此證明上述流程符合ARMA建模要求。由各分量的夜均值一階差分預測值歸算至夜均值預測值，并給出95%置信區(qū)間。Z、H、D、F各分量的原始數(shù)據(jù)夜均值后30 d的預測結果見圖6。

圖6 各分量ARMA預測值Fig.6 Prediction value of ARMA for each component

4 結論

(1)與一元線性回歸分析相比，奇異譜分析后的擬合誤差遠小于一元線性回歸分析，其殘差曲線更接近0值線，這表明奇異譜分析的擬合精度更高，擬合效果更加顯著。

(2)通過對奇異譜分析后的序列進行ARMA模型預測得出，預測值符合原始序列的趨勢變化，且預測值未超出置信區(qū)間范圍，表明短期內(nèi)將無地磁異常發(fā)生。