新疆北天山地區(qū)破年變異常識別及其效能檢驗①

顏 龍， 史勇軍， 閆 瑋， 雷 晴， 朱治國

(新疆維吾爾自治區(qū)地震局，新疆 烏魯木齊 830011)

地震孕育過程要經(jīng)歷彈性變形和非彈性變形兩個重要階段，其中非彈性變形可能引起巖石物理性質(zhì)的顯著變化以及伴隨地殼形變、地下水等非線性、非均勻的變化[1-3]。這些低頻非線性過程大都受氣壓負(fù)載、地下水負(fù)載、地表氣溫以及地殼基巖熱脹冷縮等物理效應(yīng)的共同作用，尤其以季節(jié)冷熱變化帶來的周期成分尤為顯著，并表現(xiàn)為年周期震蕩運動。已有研究表明，地震發(fā)生前不僅可能有短臨異常出現(xiàn)，當(dāng)應(yīng)力積累到一定程度，也可能會出現(xiàn)大量破年變異常，這些破年變異常在地震預(yù)測過程中可能具有中長期尺度的指示意義[4]。例如，米亞齊金等利用彈性空間孕震模型分析前兆異常在孕震空間尺度與震級之間的關(guān)系[5]，認(rèn)為地震發(fā)生前不僅會有短臨異常出現(xiàn)，當(dāng)應(yīng)力積累到一定程度，也會出現(xiàn)趨勢性轉(zhuǎn)折和破年變異常。劉耀煒等[6]提出用月均值計算月變差率方法來定量描述破年變異常,以此方法確定本年相對上年在同一月的變化程度。劉冬英等[7]針對大華北地區(qū)地下流體資料的轉(zhuǎn)折時間點和形態(tài)特征進行分析，并指出異常集中分布的區(qū)域可能是當(dāng)前或今后中強地震發(fā)生的主體區(qū)域。以往研究中，多是通過對數(shù)據(jù)曲線形態(tài)的經(jīng)驗辨別，將不同學(xué)科、不同手段的觀測資料進行聯(lián)合對比[8-10]。但不同地震前的異常測項不同，各測項的異常形態(tài)、幅度和持續(xù)時間也不盡相同，異常的判斷存在明顯的主觀性。本文中擬對新疆北天山地區(qū)具有年動態(tài)的定點形變資料和地下流體資料進行系統(tǒng)梳理，開展破年變異常定量識別研究，提取出異常信息，并在此基礎(chǔ)上分析新疆北天山地區(qū)破年變異常的時空演化特征，為客觀識別出地震異常信息提供參考依據(jù)。

1 研究區(qū)域概況

新疆地球物理觀測臺站主要沿南北天山地震帶布設(shè)，在昆侖山、阿爾泰地震帶上也分布有個別臺站，以北天山地震帶中段和南天山地震帶中、西段分布最為密集，北天山地震帶歷史上中強震多發(fā)，是新疆地震監(jiān)測的重中之重，該區(qū)域主要以形變和地下流體等多項監(jiān)測為主，廣泛分布于烏魯木齊—奎屯—溫泉線(圖1)。

圖1 新疆北天山地區(qū)地球物理觀測點分布圖Fig.1 Distribution map of precursor observation point in North Tianshan area，Xinjiang

為了進一步定量分析北天山地區(qū)地球物理觀測資料中具有年變異常的測項及其預(yù)報效能，本文中基于地球物理觀測資料預(yù)報效能評估結(jié)果，選取預(yù)報效能B級以上的資料，包括5個臺站的10個定點形變和10個臺站的10個地下流體測項。其中A類資料9項(水平擺3項、水位2項、水氡3項、流量1項)，B類資料11項(水管儀3項、洞體應(yīng)變3項、水平擺1項、水位2項、水氡1項，流量1項)(表1)。

表1 新疆北天山地區(qū)地球物理觀測測項預(yù)報效能評估表

2 研究方法

根據(jù)不同地球物理觀測資料變化特征選取適用的年變擬合方法，即制作標(biāo)準(zhǔn)年變模型，然后利用R值方法進行預(yù)報效能檢驗，提取異常閾值。首先采用破年變異常自動識別軟件ADDv1.exe，該軟件主要利用4種數(shù)學(xué)模型(矩平法、傅里葉滑動法、最小二乘法、分段最小二乘法)對具有年動態(tài)的觀測數(shù)據(jù)制作 “標(biāo)準(zhǔn)年變”模型，即年變擬合曲線，以此獲取觀測數(shù)據(jù)清晰可識別的年變形態(tài)，然后根據(jù)擬合曲線和觀測數(shù)據(jù)的殘差判斷異常，即殘差序列中高于閾值部分，再將破年變異常的測項進行R值檢驗，提取異常指標(biāo)。

R值評分方法由許紹燮提出[11]，該方法在時間序列中的研究較多，20世紀(jì)90年代以來，該方法多被用于對歷年年度預(yù)測效能的檢驗，之后又被多次完善[12-14]，在檢驗前兆異常和地震預(yù)報效能得到了普遍應(yīng)用和認(rèn)可。

(1)

式中，當(dāng)R>0表示預(yù)報有效，R值越大效果越佳，R=1表示預(yù)測完全正確，所有地震全部落入危險區(qū)，且不存在任何虛報。通過對一段時間不同目標(biāo)震級地震進行若干次R值動態(tài)檢驗，可以得到各測項異常的最佳閾值。

3 破年變異常檢驗示例

根據(jù)表1中新疆北天山地區(qū)地球物理觀測測項，本文中將分2個震級檔分別進行檢驗。測點震中距300 km范圍內(nèi)的MS5.0地震和震中距500 km范圍內(nèi)的MS≥6.0地震，各地球物理觀測數(shù)據(jù)均截止到2019年12月31日。

3.1 阜康洞體應(yīng)變EW分量(矩平法)

阜康洞體應(yīng)變自觀測以來積累近10年的觀測資料，其NS分量無顯著年變動態(tài)，而EW分量年變形態(tài)較顯著。圖2(a)對阜康洞體應(yīng)變EW分量進行去趨勢處理，黑色實線為EW分量預(yù)處理后的觀測曲線，紅色虛線為通過矩平法計算出的年變擬合曲線。而圖2(d)中，色標(biāo)差異則代表預(yù)報效能大小，黑色圓點代表對應(yīng)最佳閾值和預(yù)測時間的最優(yōu)R值，通過對測點周邊MS≥6.0地震進行檢驗，可以看出阜康洞體應(yīng)變EW分量異常出現(xiàn)后102 d為優(yōu)勢對應(yīng)時間，最佳閾值為1.2，對應(yīng)的預(yù)報效能R=0.418 3。圖2(c)為觀測值和擬合值的殘差絕對值曲線，紅色虛線代表殘差序列所帶入閾值，從全部篩選的MS≥6.0地震檢驗結(jié)果來看，檢驗效果不是很理想，僅報準(zhǔn)2016年12月8日呼圖壁MS6.2地震和2017年8月9日精河MS6.6地震[15-16]，其原因可能是阜康洞體應(yīng)變EW分量長時間受天池庫容變化引起的載荷變化影響，剔除明確的干擾時段后，在有效的預(yù)報時間內(nèi)僅發(fā)生了4次MS≥6.0地震，該資料積累時間較短也是預(yù)報效能較低的可能原因之一。

3.2 溫泉水平擺傾斜EW分量(傅里葉滑動法)

溫泉水平擺傾斜EW分量數(shù)據(jù)存在11個完整的年變，年變振幅大小雖有小幅變化，但周期形態(tài)較為穩(wěn)定[17]。通過傅里葉滑動方法擬合標(biāo)準(zhǔn)年變選取周期T=365 d，采用3T窗長滑動，擬合結(jié)果如圖3(a)，可以看出該擬合結(jié)果和原始曲線較貼合。利用R-TT曲線(效能-預(yù)測時間/閾值)分析發(fā)現(xiàn)溫泉水平擺傾斜EW分量對應(yīng)周邊300 km內(nèi)的MS5.0地震具有較好的映震效能，但無優(yōu)勢集中區(qū)域，地震分布如圖3(b)。圖3(d)的檢驗結(jié)果中，當(dāng)閾值取1.5時，預(yù)報優(yōu)勢時效為82 d，對應(yīng)的最佳預(yù)報效能R=0.802 5，通過了預(yù)報效能檢驗。圖3(c)中紅色虛線上的綠色方塊代表超過殘差閾值的異常對應(yīng)時段，可以看出溫泉水平擺傾斜EW分量破年變異?？偣渤霈F(xiàn)了9次異常時段，在此期間總共發(fā)生了9次MS5.0地震，報準(zhǔn)了其中的5次地震，說明溫泉水平擺傾斜EW分量在出現(xiàn)破年變異常后發(fā)生地震的可能性較大。

圖2 阜康洞體應(yīng)變EW分量年變識別結(jié)果(a) 年變曲線擬合圖 (b) 震例空間分布圖 (c) 異常自動識別結(jié)果 (d) R-TT曲線Fig.2 Annual change anomaly recognition results of EW component of cavern strain in Fukang

圖3 溫泉水平擺傾斜EW分量年變識別結(jié)果(a) 年變曲線擬合圖 (b) 震例空間分布圖 (c) 異常自動識別結(jié)果 (d) R-TT曲線Fig.3 Annual change anomaly recognition results of EW component of horizontal pendulum tilt in Wenquan

3.3 榆樹溝洞體應(yīng)變NS分量(最小二乘法)

榆樹溝洞體應(yīng)變NS分量數(shù)據(jù)累計超過12年，且年變幅周期和振幅均無明顯規(guī)律。由于該點NS分量數(shù)據(jù)年變無顯著的分段點，因此利用最小二乘法和分段最小二乘法年變擬合后效果一致。利用年變擬合殘差對該測項500 km內(nèi)的MS6.0地震進行檢驗如圖4(d)，R-TT曲線顯示，當(dāng)閾值取0.6時，預(yù)報優(yōu)勢時效為62 d內(nèi)，對應(yīng)的最佳預(yù)報效能R=0.4612。從圖4(c)中看出，資料的異常時間段較多，有效預(yù)報時間段內(nèi)發(fā)生28次MS6.0地震，但僅報準(zhǔn)4次，并且這4次報準(zhǔn)的地震并非發(fā)生在最大異常之后。從空間分布上來看，4次地震都在該測點西側(cè)且無優(yōu)勢區(qū)域，該測項的破年變異常的信度較低。

圖4 榆樹溝洞體應(yīng)變NS分量年變識別結(jié)果(a) 年變曲線擬合圖 (b) 震例空間分布圖 (c) 異常自動識別結(jié)果 (d) R-TT曲線Fig.4 Annual change anomaly recognition results of NS component of cavern strain in Yushugou

3.4 新10井水位(分段最小二乘法)

新10井水位自觀測以來無明顯年變規(guī)律。2000年以來的變化量均不相同，多處數(shù)據(jù)還存在相位改變的的情況[18]，但在2011年前后存在年變幅差別，2011年后整體年變的振幅明顯增大。根據(jù)新10井水位的變化特征，可采用分段曲線擬合，在每段區(qū)間上進行局部最小二乘擬合，然后通過F檢驗判斷分段點的顯著性，在顯著性水平α>0.2的水平時，即為分段點。圖5(a)可以看出，年變擬合曲線較好的適應(yīng)原始數(shù)據(jù)2011年前后的年變幅變化，分段處理可有效減少擬合曲線與原始數(shù)據(jù)的殘差序列。利用殘差序列對區(qū)域300 km范圍內(nèi)的MS5.0地震進行預(yù)報效能檢驗如圖5(d)所示，預(yù)報效能最高為R=0.592 7，閾值為2.2。有效預(yù)報時間段內(nèi)總共發(fā)生10次MS5.0地震，報準(zhǔn)7次，291 d內(nèi)的預(yù)測效果最好。綜合考慮整體預(yù)報效能和周邊地震分布特征，新10井水位破年變異常發(fā)生后對周邊MS5.0地震的預(yù)報效能較好。

圖5 新10井水位年變識別結(jié)果(a) 年變曲線擬合圖 (b) 震例空間分布圖 (c) 異常自動識別結(jié)果 (d) R-TT曲線Fig.5 Annual change anomaly recognition results of Xin 10 well

4 檢驗結(jié)果

參照上述示例，分別利用4種擬合方法對表1中的測項進行年變擬合，并利用殘差序列定量化提取異常指標(biāo)：對應(yīng)時間、預(yù)報效能、最佳閾值(表2)，主要包括北天山區(qū)域10個定點形變的20分量測項和10個地下流體測項，其中阜康洞體應(yīng)變NS分量由于年變不顯著，暫不做分析。

按照對應(yīng)的空間(天山以北、周邊)劃分，地震優(yōu)勢對應(yīng)空間為天山以北范圍的測項有23項(占研究測項總數(shù)的76.7%)，無優(yōu)勢對應(yīng)區(qū)域的7項(占研究測項總數(shù)的23.3%)。按照對應(yīng)的時間劃分，優(yōu)勢對應(yīng)時間段為短期的測項有22項(占研究測項總數(shù)的73.3%)，優(yōu)勢對應(yīng)時間段為中期的6項(占研究測項總數(shù)的20.0%)，優(yōu)勢對應(yīng)時間段為長期的2項(占研究測項總數(shù)的6.7%)。按照對應(yīng)的震級劃分，破年變異常優(yōu)勢對應(yīng)為MS≥6.0地震的測項有22項(占研究測項總數(shù)的73.3%)，優(yōu)勢對應(yīng)為MS5.0地震的8項(占研究測項總數(shù)的26.7%)。按照測項所適應(yīng)的擬合方法劃分，適應(yīng)矩平法擬合的測項有3項(占研究測項總數(shù)的10.0%)，適應(yīng)傅里葉滑動法擬合的17項(占研究測項總數(shù)的56.7%)，適應(yīng)最小二乘法擬合的4項(占研究測項總數(shù)的13.3%)，適應(yīng)分段最小二乘法擬合的6項(占研究測項總數(shù)的20.0%)。從整體研究結(jié)果看，多數(shù)地球物理觀測測項對北天山地區(qū)中強震預(yù)測效果較好，對應(yīng)時間以中短期居多。定點形變測項較適應(yīng)選擇最小二乘法和分段最小二乘法進行年變擬合，地下流體測項較適用傅里葉滑動法進行年變擬合。

表2 新疆北天山地區(qū)地球物理觀測測項破年變異常預(yù)報效能評估與預(yù)報指標(biāo)

5 結(jié)論與討論

破年變異常作為地震預(yù)測中常見的異常形態(tài)，如果只利用單學(xué)科單手段進行預(yù)報，往往可能出現(xiàn)虛報漏報，開展多學(xué)科破年變異常特征研究可以更好地為各前兆異常信度進行相互佐證，大大消除干擾，降低前兆異常判定的復(fù)雜度。通過對北天山區(qū)域的定點形變和地下流體中具有年動態(tài)的觀測資料進行破年變異常分析，選用4種不同的擬合方法進行年變周期擬合，得到正常年變時間擬合序列，從分級震例統(tǒng)計角度，得到最佳預(yù)測效能的參數(shù)組合，即異常最佳閾值和預(yù)測時間，并嘗試探討該區(qū)域地球物理觀測資料破年變異常的時空特征。

(1) 在識別破年變異常過程中，針對地球物理資料各自變化特征需要選取不同年變擬合方法。矩平法適合積累時間較短年變形態(tài)變化相對較少的資料，此方法對大多已積累10年以上年變形態(tài)復(fù)雜多樣的觀測資料來說，擬合度相對過低，得到殘差序列取值過大，且不利于異常的提取。最小二乘擬合法對年變周期比較規(guī)律的資料擬合度較高，但對于振幅多變的觀測資料擬合度就較差，若振幅變化存在明顯前后分段現(xiàn)象，則可使用分段最小二乘法尋找出顯著分段點，這樣得到的擬合結(jié)果更接近于原始觀測數(shù)據(jù)。傅里葉滑動法在觀測數(shù)據(jù)較長年變形態(tài)比較復(fù)雜且存在非規(guī)律的振幅改變時擬合效果最佳，一般選取1年為周期，3年窗長進行滑動，這樣可以較好避免某年年變幅異常帶來的擬合不穩(wěn)定。

(2) 從異常對應(yīng)時間來看，北天山地區(qū)地球物理定點觀測資料的破年變異常大多為中短期異常，這與以往分析人員利用數(shù)據(jù)曲線形態(tài)的經(jīng)驗辨別的破年變異常震例結(jié)果較為一致，即破年變異常對應(yīng)的地震以半年內(nèi)的居多。說明破年變異常一定程度上反映了區(qū)域構(gòu)造應(yīng)力狀態(tài)與異常所處區(qū)域整體動力學(xué)背景的調(diào)整和演化特征，此類異常對區(qū)域內(nèi)中強地震的對應(yīng)性較好。

(3) 從對應(yīng)空間上來看，主要優(yōu)勢對應(yīng)區(qū)域為北天山地區(qū)。其原因是所選地球物理觀測測項主要沿北天山地震構(gòu)造帶分布，該區(qū)域中強地震頻發(fā)。因此，對北天山區(qū)域地震檢驗效果較好，雖然部分南天山地震也在所選空間范圍內(nèi)，但從檢驗結(jié)果上來看，對南天山區(qū)域的預(yù)報效能較差。

北天山地區(qū)沿構(gòu)造具有明顯的分帶性，強震構(gòu)造帶之間力學(xué)性質(zhì)不盡相同，不同測項對不同構(gòu)造帶上的地震預(yù)報效能也各不相同。定量化破年變異常識別不僅能直觀的對各測項整體預(yù)報效能進行評價，還能克服人為因素影響，更為客觀和科學(xué)的顯示定量化結(jié)果，具有一定優(yōu)勢，但單從震級和震中距兩個方面作為約束條件，不同的測項和不同方向分量對周邊地震的反應(yīng)都存在其特異性，需要結(jié)合具體的地質(zhì)構(gòu)造背景等進一步深入研究。