苑爭(zhēng)一，閆 偉，牛安福，趙 靜，高 歌

(1.中國(guó)地震臺(tái)網(wǎng)中心，北京 100045；2.新疆維吾爾自治區(qū)地震局，新疆 烏魯木齊 830011)

0 引言

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

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

SSA方法在對(duì)時(shí)間序列趨勢(shì)項(xiàng)和周期項(xiàng)的識(shí)別、提取等方面具有顯著的優(yōu)勢(shì)。因此，利用SSA方法進(jìn)行背景年變序列的擬合研究，既能克服傳統(tǒng)距平法年變擬合振幅不準(zhǔn)確的缺陷，又能克服滑動(dòng)傅里葉變換方法嚴(yán)重依賴于正余弦重構(gòu)的不穩(wěn)定性問(wèn)題。本文基于SSA方法，以新疆東風(fēng)煤礦鉆孔傾斜EW分量和巴里坤水平擺傾斜NS分量為例，在去除典型干擾及長(zhǎng)周期趨勢(shì)變化的基礎(chǔ)上，對(duì)觀測(cè)資料背景年變序列進(jìn)行擬合。

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

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

1.1 構(gòu)建時(shí)間滯后矩陣X

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

(1)

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

1.2 計(jì)算奇異值及特征向量

設(shè)矩陣B=XXT，其中XT為軌跡矩陣X的轉(zhuǎn)置。接下來(lái)計(jì)算矩陣B的所有特征值λi及其特征向量Ui，按照降序?qū)⑻卣髦蹬帕袨棣?，λ2，…，λL(λ1≥λ2≥…≥λL≥0)，其對(duì)應(yīng)的特征向量分別為U1，U2，…，UL。

(2)

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

X=X1+X2+…+Xd

(3)

1.3 分組

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

(4)

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

1.4 對(duì)角平均

通過(guò)上述計(jì)算得到的重構(gòu)矩陣X*，需要進(jìn)行對(duì)角平均計(jì)算(Zhigljavsky，2010)，從而得到一個(gè)重構(gòu)的時(shí)間序列y1，y2，…，yN，展開(kāi)后的計(jì)算方法為：

(5)

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

2 R值檢驗(yàn)基本原理及其計(jì)算方法

R值評(píng)分方法(許紹燮，1989；張國(guó)民等，2002；閆偉等，2019)是當(dāng)前普遍認(rèn)可和使用的地震效能評(píng)價(jià)方法，其計(jì)算方法如下式：

(6)

式中：c1為報(bào)準(zhǔn)地震的次數(shù)；c2為應(yīng)預(yù)測(cè)的地震總次數(shù)；t1為預(yù)測(cè)占用時(shí)間；t2為預(yù)測(cè)總時(shí)間。其基本原理如圖1所示。R>0表示預(yù)報(bào)有效，且R值越大效果越好，最大值為1；R=0表示預(yù)報(bào)無(wú)意義；R<0為預(yù)報(bào)無(wú)效。另外結(jié)合統(tǒng)計(jì)分布檢驗(yàn)，利用報(bào)準(zhǔn)、漏報(bào)地震的個(gè)數(shù)，可以計(jì)算出大于97.5%置信度的R臨界值R0，當(dāng)R>R0時(shí)，認(rèn)為其預(yù)測(cè)的有效性高于自然發(fā)生概率，該測(cè)項(xiàng)具有較高的預(yù)測(cè)效能。

3 技術(shù)路線

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

圖2 算法流程圖

(1)首先加載數(shù)據(jù)，并對(duì)原始觀測(cè)數(shù)據(jù)進(jìn)行預(yù)處理，包括去突跳、臺(tái)階及插值補(bǔ)缺等。

(2)對(duì)預(yù)處理數(shù)據(jù)進(jìn)行滑動(dòng)傅里葉周期顯著性檢驗(yàn)，若不顯著，先進(jìn)行高階小波去趨勢(shì)處理，進(jìn)而將具有顯著周期的觀測(cè)數(shù)據(jù)進(jìn)行奇異譜分析(SSA)。

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

(4)設(shè)置初始預(yù)測(cè)時(shí)間閾值T0，將超限的時(shí)間點(diǎn)加時(shí)間初始閾值T0得到預(yù)警時(shí)段，檢索震例并代入R值公式計(jì)算得到R00。

(5)增加預(yù)測(cè)時(shí)間步長(zhǎng)step1，循環(huán)進(jìn)行上述R值檢驗(yàn)，得到對(duì)應(yīng)于同一異常判定閾值、不同預(yù)測(cè)時(shí)間的R值序列[R00，R01，…，R0n-1]。

(6)增加異常判定閾值步長(zhǎng)step2，按照步驟(5)重新遍歷預(yù)測(cè)時(shí)間段集合并進(jìn)行R值檢驗(yàn)，得到[R10，R12，…，R1n-1]。

(7)以此遞推，求得異常判定閾值、預(yù)測(cè)時(shí)間滑動(dòng)下的二維R值序列數(shù)組，將[T1，T2，…，Tm]與[S1，S2，…，Sn]、R值序列之間的關(guān)系擬合成連續(xù)的二維曲面，取R值的最高點(diǎn)所對(duì)應(yīng)的閾值T、S作為該測(cè)項(xiàng)異常提取和地震預(yù)測(cè)的標(biāo)準(zhǔn)。

4 實(shí)驗(yàn)結(jié)果

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

4.1 數(shù)據(jù)預(yù)處理及殘差序列SSA

圖3原始觀測(cè)時(shí)間序列顯示，2個(gè)測(cè)項(xiàng)在年變

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

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

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

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

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

圖5 重構(gòu)年變成分的FFT變換結(jié)果及震例空間分布

序列做差，得到了去趨勢(shì)和年變(SSA重構(gòu)的正常年變)后的殘差時(shí)間序列，在此基礎(chǔ)上進(jìn)行動(dòng)態(tài)R值檢驗(yàn)。

4.2 震例回溯及動(dòng)態(tài)R值檢驗(yàn)

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

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

動(dòng)態(tài)R值檢驗(yàn)過(guò)程對(duì)應(yīng)于2個(gè)互相嵌套的循環(huán)，需要設(shè)置2個(gè)循環(huán)變量：

(1)異常判定閾值

通常取均值加減幾倍的均方差作為異常判定的標(biāo)準(zhǔn)，但不同觀測(cè)方法、不同臺(tái)站測(cè)項(xiàng)的信噪比差別較大，取單一的閾值會(huì)漏掉部分異常信息，或引入部分噪聲信息。因此本文將上述異常判定閾值擴(kuò)大到某一個(gè)動(dòng)態(tài)范圍，循環(huán)進(jìn)行R值檢驗(yàn)，根據(jù)檢驗(yàn)結(jié)果確定最佳閾值。本例中設(shè)定閾值范圍為1～3倍的均方差，循環(huán)步長(zhǎng)設(shè)為0.1倍的均方差。

(2)預(yù)測(cè)時(shí)間

預(yù)測(cè)時(shí)間即異常從開(kāi)始到結(jié)束所持續(xù)的時(shí)間，通過(guò)循環(huán)設(shè)定該變量并進(jìn)行R值檢驗(yàn)，可以得到某一觀測(cè)具有最佳預(yù)測(cè)效能的預(yù)測(cè)時(shí)間取值，本文設(shè)定該變量的取值范圍0～720 d，循環(huán)步長(zhǎng)為10 d。

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

圖7 東風(fēng)煤礦鉆孔傾斜EW動(dòng)態(tài)R值掃描(a)與異常自動(dòng)識(shí)別結(jié)果

圖8 巴里坤水平擺NS動(dòng)態(tài)R值掃描(a)與異常自動(dòng)識(shí)別結(jié)果(b)

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

5 結(jié)論與討論

本文通過(guò)對(duì)去除長(zhǎng)周期趨勢(shì)變化的平穩(wěn)時(shí)間序列進(jìn)行SSA分析，得到了基于SSA的正常年變時(shí)間序列的擬合；在此基礎(chǔ)上，對(duì)扣掉趨勢(shì)變化和年變成分的變化部分做動(dòng)態(tài)R值檢驗(yàn)，從震例回溯和統(tǒng)計(jì)的角度，得到了具有最佳映震效能的預(yù)測(cè)參數(shù)組合，即異常判定閾值和時(shí)間預(yù)測(cè)準(zhǔn)則。應(yīng)用該準(zhǔn)則，可以跟蹤該測(cè)項(xiàng)未來(lái)的異常變化情況，并結(jié)合R值評(píng)分和以往的震例對(duì)應(yīng)情況，給出相應(yīng)的預(yù)測(cè)意見(jiàn)，對(duì)異常判定的定量化和自動(dòng)化工作具有一定的參考價(jià)值。

基于以上研究，對(duì)未來(lái)工作提出如下的建議：

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

(2)對(duì)比不同研究方法對(duì)原始數(shù)據(jù)去趨勢(shì)處理的結(jié)果及可靠性，優(yōu)選出最佳的趨勢(shì)擬合方案，確保該處理過(guò)程的可靠性。

(3)目前震例的選取采用的是經(jīng)驗(yàn)做法，單從與臺(tái)站的距離和震級(jí)兩個(gè)方面作為約束，不同臺(tái)站對(duì)周邊地震的反映具有特異性，因此需要結(jié)合具體地質(zhì)構(gòu)造背景，利用數(shù)據(jù)挖掘等方法綜合確定針對(duì)某一特定臺(tái)站和測(cè)項(xiàng)的震例，輸入動(dòng)態(tài)R值檢驗(yàn)?zāi)K求解最佳預(yù)測(cè)參數(shù)組合。