奇異譜分析在地磁數(shù)據(jù)研究中的應(yīng)用

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

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

0 引言

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

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

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

1 奇異譜分析原理

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

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

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

(1)

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

C=DDT

(2)

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

(3)

(4)

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

(5)

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

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

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

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

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

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

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

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

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

表1 各分量前10個(gè)特征值貢獻(xiàn)率

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

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

3.2 擬合精度檢驗(yàn)

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

y=β0+β1x+ε

(6)

式中：x為時(shí)間刻度，y為夜均值擬合值，ε為誤差。經(jīng)計(jì)算不同分量的各項(xià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ù)的變化趨勢，但相比之下奇異譜分析能更精準(zhǔn)地表達(dá)出數(shù)據(jù)變化的細(xì)節(jié)。Z、H、D、F各分量的殘差對比曲線見圖5。計(jì)算各分量殘差值距0值線的偏離程度，用以比較兩種方法的擬合效果，所求結(jié)果見表3。由此可見4個(gè)分量的奇異譜分析殘差量均小于一元線性回歸分析，故奇異譜分析的擬合效果更顯著。

表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í)物理量變化、干擾噪聲和異常信號等信息。如果用原始觀測數(shù)據(jù)減掉背景數(shù)據(jù)所獲得的時(shí)間序列，則突出了背景噪聲和異常信號，這有助于分析背景噪聲并放大短臨異常信號。在本論文中原始數(shù)據(jù)與背景數(shù)據(jù)的差值序列為奇異譜分析的殘差序列，這突出了殘差序列在地磁異常研究中的重要性。

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

(7)

(8)

(9)

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

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

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

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

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

4 結(jié)論

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

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