定點形變破年變異常自動識別應用研究*

苑爭一，閆 偉，牛安福，趙 靜，高 歌

(1.中國地震臺網中心，北京 100045；2.新疆維吾爾自治區地震局，新疆 烏魯木齊 830011)

0 引言

近年來，不同學者對多種類型的定點形變異常識別、判定及映震效能檢驗工作做了大量研究(牛安福等，2013；楊玲英等，2014；崔青發等，2014；劉琦等，2016；馬棟等，2019)，破年變類異常是其中的一個重要部分。目前該類異常的判定大多采用基于分析人員主觀經驗的曲線形態辨識法，不同的人對于同一個異常的認定具有不同的看法，導致異常的判定準則及預報效能缺乏定量化的表達。隨著前兆臺網的密集布設，時間序列數據的積累量已經達到TB級，為實現海量觀測數據的快速高效提取，迫切需要研究針對不同異常類型、切實可行的自動識別算法。

破年變異常的提取以正常年變的擬合為基礎，目前常用的方法有距平法和滑動傅里葉方法。距平法年變擬合通過不同年份內同一天內的觀測取平均值構建標準背景年變序列，其振幅是固定不變的；滑動傅里葉方法受正余弦波假定的約束，對于非正弦波及非平穩的時間序列周期成分的識別不穩定。奇異譜分析(singular spectrum analysis，SSA)是在Karhumen-Loeve分解理論的基礎上發展起來的(Vautardetal，1992)，是從時間序列的動力重構出發的自適應分析方法。該方法無需先驗信息(江志紅，丁裕國，1998)，不僅能對信號的不同頻率進行識別，還可以按照信號能量的大小進行嚴格排序，穩定識別非正弦波的周期信號，比傳統的頻譜分析更具優勢(趙佳佳等，2017)。該方法適應面廣、物理意義清晰，廣泛應用于氣象學、海洋學等領域(吳洪寶，1997；余錦華，丁裕國，2001；李亞偉等，2010)。近年來在大地測量學及地形變數據處理領域也有了初步應用：①數據預處理應用方面，王解先等(2013)利用SSA分析方法對中國地殼運動觀測網絡GPS數據服務提供的測站坐標時間(NEU，以BJFS站為例)序列進行了補缺，并根據降噪重構序列提取了坐標時間序列中的趨勢和周期成分。②干擾信息提取方面，趙佳佳(2017)將SSA應用于實際傾斜應變固體潮數據的處理中，結果表明SSA能有效提取數據中的長周期成分和固體潮各潮波，有效分離數據中的干擾成分，對干擾排除和異常提取具有一定的參考意義。③周期成分重構與分析方面，李世友等(2018)利用SSA分析算法提取了大壩變形時間序列的趨勢和周期分量，分析了各分量與時效、溫度和水位因子的關聯性；張旺等(2017)改進了傳統SSA的“相移現象缺點”，提出一種新的“SSA-PD”分析方法，并應用該方法對GNSS坐標時間序列中的趨勢項和季節項進行了分析；姚宜斌等(2016)利用SSA方法提取出電離層TEC時間序列中除去異常擾動以及觀測噪聲部分的日周期部分作為主要成分，在此基礎上對震前電離層異常進行了探測與分析。

SSA方法在對時間序列趨勢項和周期項的識別、提取等方面具有顯著的優勢。因此，利用SSA方法進行背景年變序列的擬合研究，既能克服傳統距平法年變擬合振幅不準確的缺陷，又能克服滑動傅里葉變換方法嚴重依賴于正余弦重構的不穩定性問題。本文基于SSA方法，以新疆東風煤礦鉆孔傾斜EW分量和巴里坤水平擺傾斜NS分量為例，在去除典型干擾及長周期趨勢變化的基礎上，對觀測資料背景年變序列進行擬合。

1 奇異譜分析(SSA)基本原理計算方法

SSA的實現過程主要分為構建時間滯后矩陣(軌跡矩陣)、奇異值分解(SVD)、分組重構、對角平均4大步驟(趙佳佳等，2017)，具體操作過程如下。

1.1 構建時間滯后矩陣X

給定嵌入長度L，即窗口長度，構建給定時間序列Xi=(xi，…，xi+L-1)T(1≤i≤k)的軌跡矩陣XL×K，其中K=N-L+1，1

(1)

SSA實現的過程需要確定2個重要參數——嵌入維度L和重構階次P。嵌入維度L的選擇直接影響SSA分解信號的效果，研究表明：若L較小，將會導致信號分解不徹底；反之，將會使低頻信號的分解產生譜裂現象(Vautardetal，1992)。對于待提取年變信號的觀測數據，L的選取一般為年變周期的整數倍，這樣分解效果較好。經過試驗，本文設置L為整周年天數365 d。原始曲線如果存在顯著的趨勢，先利用高階小波去除長周期趨勢項，根據奇異譜曲線動態選取重構階次P，確保重構序列位于年頻段范圍。

1.2 計算奇異值及特征向量

設矩陣B=XXT，其中XT為軌跡矩陣X的轉置。接下來計算矩陣B的所有特征值λi及其特征向量Ui，按照降序將特征值排列為λ1，λ2，…，λL(λ1≥λ2≥…≥λL≥0)，其對應的特征向量分別為U1，U2，…，UL。

(2)

式中：Ui和Vi是軌跡矩陣X的左右特征向量；d=L*；L*=min{L，K}；軌跡矩陣X可以由初等矩陣按照下式合成：

X=X1+X2+…+Xd

(3)

1.3 分組

將初等矩陣Xi的下標(i=1，2，…，d)分成p個不相交的子集I1，I2，…，Ip， 設I={i1，i2，…，im}， 則：

(4)

選取集合I1，I2，…，Ip的過程即為分組。研究表明，奇異值接近的2個初等矩陣重構得到一個具有優勢周期的成分。通過求解并觀測奇異值的分布情況，對觀測時間序列分組重構，進而通過快速傅里葉變化(FFT)對重構的年變序列進行周期檢驗，確保其優勢周期位于年頻段范圍內。

1.4 對角平均

通過上述計算得到的重構矩陣X*，需要進行對角平均計算(Zhigljavsky，2010)，從而得到一個重構的時間序列y1，y2，…，yN，展開后的計算方法為：

(5)

式中：L*=min{L，K}，K*=max{L，K}。

2 R值檢驗基本原理及其計算方法

R值評分方法(許紹燮，1989；張國民等，2002；閆偉等，2019)是當前普遍認可和使用的地震效能評價方法，其計算方法如下式：

(6)

式中：c1為報準地震的次數；c2為應預測的地震總次數；t1為預測占用時間；t2為預測總時間。其基本原理如圖1所示。R>0表示預報有效，且R值越大效果越好，最大值為1；R=0表示預報無意義；R<0為預報無效。另外結合統計分布檢驗，利用報準、漏報地震的個數，可以計算出大于97.5%置信度的R臨界值R0，當R>R0時，認為其預測的有效性高于自然發生概率，該測項具有較高的預測效能。

3 技術路線

本文的基本思路和算法流程如圖2所示：

圖2 算法流程圖

(1)首先加載數據，并對原始觀測數據進行預處理，包括去突跳、臺階及插值補缺等。

(2)對預處理數據進行滑動傅里葉周期顯著性檢驗，若不顯著，先進行高階小波去趨勢處理，進而將具有顯著周期的觀測數據進行奇異譜分析(SSA)。

(3)設置初始的異常判定閾值S0，按此閾值提取殘差時間序列(原始觀測值減去趨勢和SSA擬合年變成分)的超限部分。

(4)設置初始預測時間閾值T0，將超限的時間點加時間初始閾值T0得到預警時段，檢索震例并代入R值公式計算得到R00。

(5)增加預測時間步長step1，循環進行上述R值檢驗，得到對應于同一異常判定閾值、不同預測時間的R值序列[R00，R01，…，R0n-1]。

(6)增加異常判定閾值步長step2，按照步驟(5)重新遍歷預測時間段集合并進行R值檢驗，得到[R10，R12，…，R1n-1]。

(7)以此遞推，求得異常判定閾值、預測時間滑動下的二維R值序列數組，將[T1，T2，…，Tm]與[S1，S2，…，Sn]、R值序列之間的關系擬合成連續的二維曲面，取R值的最高點所對應的閾值T、S作為該測項異常提取和地震預測的標準。

4 實驗結果

奇異值的選取需根據奇異譜曲線的形態確定，因此本文設計了一套可以交互選取奇異值進行年變成分重構的計算程序。該程序主要分為“奇異譜分析SSA”和“R值遞歸檢驗”兩大模塊，其中SSA模塊提供了數據預處理(包括去除臺階、突跳和趨勢項)、預覽、奇異值選取、重構序列周期檢驗以及原始及處理結果可視化幾個功能，支持交互確定重構參數。以新疆東風煤礦鉆孔傾斜EW和目前存在顯著破年變異常的巴里坤水平擺傾斜NS分量為例，簡要介紹提取過程并分析結果。

4.1 數據預處理及殘差序列SSA

圖3原始觀測時間序列顯示，2個測項在年變

(a)東風煤礦鉆孔傾斜EW (b)巴里坤水平擺傾斜NS

的基礎上疊加了顯著的趨勢項變化，其中東風煤礦鉆孔傾斜EW存在單一的西傾趨勢(圖3a)，而巴里坤水平擺傾斜NS存在多個線性變化趨勢(圖3b)。因此，首先利用db5小波進行高階濾波，去除長周期的趨勢變化，進而對變化相對平穩的去趨勢時間序列做奇異譜分析。

去趨勢時間序列的奇異譜曲線顯示(圖4a)，2個測項的第1，2個奇異值位于噪聲平臺的上方，占比較高且存在顯著的周期臺階，因此取重構階次P為[1，2]，重構結果見圖4b，c中紅色和藍色虛線部分，灰色曲線為去趨勢后的觀測時間序列。從圖4中可以看出，東風煤礦鉆孔傾斜EW分量去掉趨勢項后，部分年份如2013，2016年，明顯偏離SSA重構得到的年周期變化規律，出現破年變信號；巴里坤水平擺NS的多個長周期趨勢項，使得趨勢轉折與破年變變化混在一起難以區分，通過去趨勢并結合SSA重構的結果來看，2017，2018以及2019年均存在顯著的破年變異常，與筆者從原始曲線觀測到的破年變異常年份基本一致，只是幅度略有差異。

圖4 奇異譜曲線(a)及SSA年變重構結果(b)(c)

對SSA算法重構得到的年變時間序列做快速傅里葉變換(FFT)的結果顯示(圖5)，2個測項通過SSA重構后的時間序列，其優勢周期均為372 d，這表明重構時間序列成功分離出了原始曲線中年頻段的成分。將去趨勢的觀測序列與重構時間

圖5 重構年變成分的FFT變換結果及震例空間分布

序列做差，得到了去趨勢和年變(SSA重構的正常年變)后的殘差時間序列，在此基礎上進行動態R值檢驗。

4.2 震例回溯及動態R值檢驗

取臺站周邊200 km范圍內5級以上、300 km范圍內6級以上地震(Dobrovolskyetal，1979)，作為輸入震例進行R值檢驗，檢索到的震例及臺站的空間分布情況如圖6所示，2個藍色圓圈分別代表了以臺站為中心的200，300 km的圓形震例檢索區。

(a)東風煤礦鉆孔傾斜EW (b)巴里坤水平擺傾斜N

動態R值檢驗過程對應于2個互相嵌套的循環，需要設置2個循環變量：

(1)異常判定閾值

通常取均值加減幾倍的均方差作為異常判定的標準，但不同觀測方法、不同臺站測項的信噪比差別較大，取單一的閾值會漏掉部分異常信息，或引入部分噪聲信息。因此本文將上述異常判定閾值擴大到某一個動態范圍，循環進行R值檢驗，根據檢驗結果確定最佳閾值。本例中設定閾值范圍為1～3倍的均方差，循環步長設為0.1倍的均方差。

(2)預測時間

預測時間即異常從開始到結束所持續的時間，通過循環設定該變量并進行R值檢驗，可以得到某一觀測具有最佳預測效能的預測時間取值，本文設定該變量的取值范圍0～720 d，循環步長為10 d。

循環計算后得到了“R值-異常判定閾值-預測時間”的空間插值圖像(圖7a，8a)，其色標代表了不同的異常判定閾值-預測時間組合所對應的預測效能檢驗結果，即R值評分，取值范圍-1～1。從中可以看到不同組合的明顯差異性，圖中黑色圓點代表了R值最高點，將對應的異常判定閾值及預測時間代入殘差時間序列，提取異常特征時段進行R值檢驗，得到了具有最佳映震效能的異常識別及震例回溯結果如圖7b，8b所示。

圖7 東風煤礦鉆孔傾斜EW動態R值掃描(a)與異常自動識別結果

圖8 巴里坤水平擺NS動態R值掃描(a)與異常自動識別結果(b)

圖7b，8b中藍色虛線為具有最佳R值評分的異常判定閾值線，黑色實線為數據均值線。從動態檢驗的結果來看：①上述2個測項具有最佳預測效能的異常判定閾值分別為均值加減1.6倍、2倍的均方差；②異常出現后最佳預測時間分別為91 d，61 d，目前2個測項均處在異常階段；③東風煤礦鉆孔傾斜EW的破年變異常對其周邊300 km內的5級以上地震具有較好的映震效能，震例主要集中在臺站的ES方向；④巴里坤水平擺傾斜NS破年異常主要對應200 km內的5級以上地震，震例較少且無優勢集中方向；⑤應用本研究中的方法提取的破年變異常，東風煤礦鉆孔傾斜EW破年變異常的預測效能大于97.5%置信度的臨界值R0，預測結果較為可靠；而巴里坤水平擺傾斜NS虛報率較高，小于臨界值R0，雖然目前處于異常時段，但信度較低，建議參考使用。

5 結論與討論

本文通過對去除長周期趨勢變化的平穩時間序列進行SSA分析，得到了基于SSA的正常年變時間序列的擬合；在此基礎上，對扣掉趨勢變化和年變成分的變化部分做動態R值檢驗，從震例回溯和統計的角度，得到了具有最佳映震效能的預測參數組合，即異常判定閾值和時間預測準則。應用該準則，可以跟蹤該測項未來的異常變化情況，并結合R值評分和以往的震例對應情況，給出相應的預測意見，對異常判定的定量化和自動化工作具有一定的參考價值。

基于以上研究，對未來工作提出如下的建議：

(1)進一步論證SSA算法中嵌入維度L和重構階次P對年變重構的影響，并給出更加合理的參數設置。

(2)對比不同研究方法對原始數據去趨勢處理的結果及可靠性，優選出最佳的趨勢擬合方案，確保該處理過程的可靠性。

(3)目前震例的選取采用的是經驗做法，單從與臺站的距離和震級兩個方面作為約束，不同臺站對周邊地震的反映具有特異性，因此需要結合具體地質構造背景，利用數據挖掘等方法綜合確定針對某一特定臺站和測項的震例，輸入動態R值檢驗模塊求解最佳預測參數組合。