基于STL和GESD算法的于田兩次Ms7.3地震熱異常時空分布特征識別研究

米日阿依·買土地， 玉素甫江·如素力， 喀迪阿依·阿力木， 艾爾肯·圖爾孫

（1.新疆師范大學(xué)地理科學(xué)與旅游學(xué)院，烏魯木齊 830054；2.新疆干旱區(qū)湖泊環(huán)境與資源重點(diǎn)實(shí)驗(yàn)室，烏魯木齊 830054）

地震是地殼快速釋放能量過程中造成的振動. 地震給人類所帶來的損失很大. 為了盡可能地減少地震災(zāi)害造成的損失，相關(guān)領(lǐng)域的學(xué)者多年一直專注于地震前兆信息提取和地震預(yù)測研究［1-3］. 地震前發(fā)生的，所有與地震有關(guān)系的一切不正?，F(xiàn)象一般稱為地震前兆，分為宏觀前兆和微觀前兆. 常見的地震前兆現(xiàn)象有：地震活動異常，地震波速度變化，地下水不正常改變，重力異常，動物異常，地聲，地光，地表溫度（LST，Land Surface Temperate）和地表亮溫（Brightness Temperature）異常等［4-5］. 其中，震前亮度溫度，地表溫度（LST）異常是利用熱紅外遙感監(jiān)測地震前兆的主要參量.

最初對地球內(nèi)部的熱紅外異常與地震活動的對應(yīng)關(guān)系進(jìn)行探討的是蘇聯(lián)Gorny 等，隨著空間技術(shù)的發(fā)展，采用熱紅外遙感數(shù)據(jù)監(jiān)測震前LST異常成為可能［6-7］. 熱紅外遙感在地震監(jiān)測以及地震預(yù)報中的應(yīng)用越來越豐富. 地震所涉及的異常中存在很多復(fù)雜問題，怎樣用有效方法，從地震熱異常中獲取地震相關(guān)的異常信息也是個難點(diǎn). 最前期的熱紅外異?，F(xiàn)象主要是通過目視判讀發(fā)現(xiàn)的，現(xiàn)在跟著研究的需求，一個個又復(fù)雜又精細(xì)的研究思路，研究方法開始出現(xiàn)在地震熱紅外信息的提取研究當(dāng)中. 當(dāng)今，涉及遙感數(shù)據(jù)檢測地震熱異常的主要方法有背景場差值法［8］、小波分析法［9-11］、異常點(diǎn)比值時序法［12］、RST（Robust Satellite Techniques）等［13］. 目前為止，不同算法在不同地震活動異常信息監(jiān)測中的作用和效果都不同，提高異常信息提取中的不確定性. 因此，熱異常監(jiān)測的算法及其應(yīng)用還需深入研究.

然而，LST的變化過程不僅具有明顯的季節(jié)性變化，還有年際趨勢變化的特征，是由多種因素共同作用的結(jié)果，不能直接對其進(jìn)行監(jiān)測分析，而要先對其進(jìn)行分解. 相關(guān)研究表明，時間序列分解法（STL，Seasonal-Trend Decomposition using LOESS）［14］可將時間序列分解為趨勢項(xiàng)（trend）、季節(jié)項(xiàng)（seasonal）和殘余項(xiàng)（remainder）.如果STL分解法能用于地震區(qū)地表溫度的時空數(shù)據(jù)集，有效去除趨勢和季節(jié)性變化的背景場信息，其對殘余項(xiàng)用GESD檢驗(yàn)（Generlized Extreme Studentized Deviate Test）進(jìn)一步進(jìn)行異常檢測，取出地震熱紅外異常. 這對防震減災(zāi)有著重要的實(shí)際應(yīng)用價值和現(xiàn)實(shí)指導(dǎo)意義.

本文擬利用基于MODIS LST 數(shù)據(jù)，以于田縣2008年3月21日7.3級與2014年2月12日7.3級地震為例，采用STL分解法、GESD檢驗(yàn)，以震中周圍區(qū)異常面積時空變化作為主要參量，獲取與此次地震相關(guān)的地震熱紅外異常演化特征.

1 研究區(qū)概況

于田縣位于新疆維吾爾自治區(qū)西南部，東臨民豐縣，北鄰塔克拉瑪干沙漠與沙雅縣接壤，西與策勒縣相毗鄰，南與西藏自治區(qū)改則縣、日土縣相接. 由南部山區(qū)、中部平原、北部沙漠區(qū)三大部分組成，南部山區(qū)海拔高度平均在4000 m以上.

2008 年3 月21 日新疆于田發(fā)生Ms7.3 地震. 震中（35.60°N，81.60°E），震源深度33 km. 地震發(fā)生后發(fā)生345 次余震. 2014 年2 月12 日新疆于田發(fā)生Ms7.3 地震，據(jù)中國地震臺網(wǎng)中心測定，震中位于（36.1°N，82.5°E），震源深度為12 km. 地震的前1天曾發(fā)生Ms5.4地震. Ms7.3主震后余震活動頻繁，最大余震為2月14日17時24分發(fā)生的Ms5.7地震.

野外地震地質(zhì)調(diào)查表明，2008年于田地震震中位于西昆侖山中的阿什庫勒盆地，距阿爾金斷裂南緣約50 km. 處于西昆侖地震帶與阿爾金地震帶的交會區(qū)，即塔里木盆地的南部，靠近西昆侖山區(qū)的邊緣部分.造成了全長約30 km，走向南北—北東東，同時兼有左旋走滑和正斷分量的張剪地表破裂，最大左旋位移1.8 m，最大垂直位移約2.0 m，發(fā)震構(gòu)造屬于阿爾金斷裂西南尾端的張性構(gòu)造區(qū). 美國地質(zhì)調(diào)查局（USGS）與哈佛大學(xué)CMT震源機(jī)制解結(jié)果表明，此次地震為略帶走滑的正斷層性質(zhì)，斷層走向?yàn)榻媳毕? 而2014年于田地震也是發(fā)生在阿爾金斷裂西南段，這次地震因發(fā)生在高山無人地區(qū)并未造成人員傷亡，其發(fā)震斷裂為阿爾金斷裂帶西南段，阿爾金斷裂帶是亞洲大陸內(nèi)部一條長達(dá)2000 km三疊紀(jì)以來長期左旋走滑的巖石圈斷裂. 地貌上表現(xiàn)為醒目的線狀展布特征，呈東西至北東東向，是青藏高原北邊界斷裂. 但在較短的時間內(nèi)在同一地區(qū)發(fā)生2次7級以上地震，引起了地震學(xué)者的廣泛關(guān)注［15-17］.

2 數(shù)據(jù)及其處理

MODIS也被稱之為中分辨率成像光譜儀（Moderate-resolution Imaging Spectroradiometer），具有高輻射分辨率（36通道）、中等空間分辨率及高飽和度，覆蓋范圍廣、時間分辨率為每天1～2次，具有很高的地震熱異常監(jiān)測精度. 用于研究的LST 數(shù)據(jù)來自于美國航空航天局（NSIDC，http：//nsidc.org/NASA/MODIS）網(wǎng)站下載的MODIS/Aqua 8d合成的MYD21A2夜間LST產(chǎn)品，空間分辨率為1 km，時間范圍為2003—2019年（782幅圖像）.

MYD21A2的LST數(shù)據(jù)，首先利用MYD21A2數(shù)據(jù)投影及格式轉(zhuǎn)換工具（MODIS reprojection tool，MRT）進(jìn)行拼接和重投影，投影轉(zhuǎn)換為ALBERS，橢球體選為WGS84，圖像文件轉(zhuǎn)換為GeoTIFF格式；其次對于地震熱異常提取，云層是一種干擾. 有云時MODIS數(shù)據(jù)或?yàn)殛懕頊囟?，或云頂溫度，或混合溫度，需要去除云的干擾. 所以，本文中地表溫度數(shù)據(jù)經(jīng)過云識別并做去云處理. 對于云像元干擾而形成的空值區(qū)域，為了解決此問題，低質(zhì)量和缺失的值通常被排除在分析之外，或者被各種預(yù)測方法的預(yù)測結(jié)果所取代. 許多預(yù)測方法利用了數(shù)據(jù)的時間相關(guān)性，但迄今為止，只有很少人嘗試?yán)脭?shù)據(jù)的時空性質(zhì). 為了解決云干擾而形成的缺失值區(qū)域，本文依賴R包中的gapfill缺失值預(yù)測方法來替代地表溫度數(shù)據(jù)集中的缺失值. gapfill方法是利用空間連貫性的特征趨勢以及遙感數(shù)據(jù)集的時間季節(jié)性動態(tài)來預(yù)測缺失值. 此方法優(yōu)點(diǎn)包括：①能夠預(yù)測具有很大比例的缺失值（高達(dá)50%）的數(shù)據(jù)集中的所有缺失值. 該方法結(jié)果均方根相對誤差（RMSPE，Root Mean Square Percentage Error）和平均絕對百分比誤差（MAPE，Mean Absolute Percentage Error）都較低［18］. ②填補(bǔ)缺失值完成之后，利用STL分解（stats和raster包）、gesd法異常檢測（Anomalize包）等處理檢測遙感數(shù)據(jù)集中的異常信息.

3 研究方法說明

3.1 STL分解法

STL是非參數(shù)的統(tǒng)計方法，此方法把時間序列數(shù)據(jù)分成三項(xiàng)，分別是低頻率的趨勢項(xiàng)，高頻率的周期項(xiàng)（季節(jié)項(xiàng)），不規(guī)則變化的趨勢項(xiàng)：

圖1 基于STL的地震熱異常分解流程Fig.1 Workflow of earthquake thermal anomaly decomposition by STL

3.2 GESD檢驗(yàn)

大多數(shù)情況下數(shù)據(jù)集中異常值不止一個，因此Rosner提出了GESD（泛化版ESD）. GESD（Generlized Extreme Studentized Deviate Test）是Rosner 教授基于Grubb’s Test（或ESD）改進(jìn)的識別離散值的方法，此方法中比較統(tǒng)計量Rj與臨界值λj，若Rj＞λj，該點(diǎn)為異常點(diǎn).

式中：λj為第i 大的數(shù)據(jù)臨界值；n 為離散點(diǎn)數(shù)量；tn-i-1,p是(n-i-1)自由度的t 分布右側(cè)尾概率為p 時的值.

3.3 基于STL分解和GESD檢驗(yàn)的異常檢測算法

LST 的變化過程不僅具有明顯的季節(jié)性變化，還有年際趨勢變化的特征，是由多種因素共同作用的結(jié)果，不能直接對其進(jìn)行監(jiān)測分析. 因此首先利用STL分解法，取出地表溫度空間數(shù)據(jù)集中周期分量、趨勢分量，對剩余的殘余分量（沒有季節(jié)及趨勢特征的時間序列）采用GESD檢驗(yàn)進(jìn)行異常監(jiān)測. 此數(shù)據(jù)處理過程借助R語言所實(shí)現(xiàn). 基于STL分解和GESD檢驗(yàn)的異常檢測算法的流程圖如圖2所示.

圖2 異常檢測算法流程圖Fig.2 Flow chart of anomaly detection algorithm

4 震例結(jié)果與分析

4.1 平均LST時間序列的STL分解結(jié)果

圖3是于田縣Ms7.3地震震中附近2003—2019年LST的8 d平均時間序列數(shù)據(jù)用STL分解后的結(jié)果圖.分解結(jié)果包括3 個部分：數(shù)據(jù)（data）反映時間序列數(shù)據(jù)的統(tǒng)計信息，季節(jié)項(xiàng)（seasonal）反映該時間序列的周期信息，趨勢項(xiàng)（trend）為長期趨勢，即該時間序列的發(fā)展方向，殘余項(xiàng)（remainder）為去除周期信息以及長期趨勢后的剩余部分. 從LST 夜間數(shù)據(jù)的STL 分解結(jié)果（圖3）可以看出，LST 數(shù)據(jù)7 月或8 月份達(dá)到最大，最大值約-0.24 ℃. 每年的12月和次年的1月份達(dá)到最低值，最低值約為-31.59 ℃. 季節(jié)性變化是LST年內(nèi)變化的重要特征，其變化范圍為-11～12 ℃之間. 趨勢分量的變化范圍為-13～-15 ℃左右，研究區(qū)地表溫度趨勢的總體變化幅度為2 ℃左右，一般與全球氣候變化有關(guān). LST殘余項(xiàng)的波動范圍-6～5 ℃之間，可以看出，2003—2019年間LST殘余項(xiàng)一直處于不規(guī)則的波動，而且波動幅度較大.

圖3 震中附近LST時序STL分解結(jié)果圖Fig.3 STL decomposition of time series data of LST around the earthquake center

通過對地表溫度趨勢項(xiàng)的分析可以看出，2003—2019年地表溫度具有明顯的變化. 為了明確地表溫度趨勢項(xiàng)的變化點(diǎn)，對研究區(qū)地表溫度數(shù)據(jù)進(jìn)行差分，文中第一步對地表溫度趨勢序列做一階差分，初步得到其趨勢變化點(diǎn). 然后對一階差分序列再做一次一階差分，使數(shù)據(jù)保持在0值附近. 本文界限為二次差分序列的2倍標(biāo)準(zhǔn)差. 超過上界，下界就歸為趨勢發(fā)生變化. 看圖4可知，2008年1月前后至2008年3月21日地震發(fā)生前地表溫度具有上升趨勢，2014年2月前后研究區(qū)地表溫度偏增，呈上升趨勢. 這兩次地表溫度上升趨勢時間點(diǎn)臨近地震發(fā)震期. 初步認(rèn)為可能與研究區(qū)2008年與2014年7.3Ms地震有一定的關(guān)系.

圖4 地表溫度趨勢項(xiàng)的二次差分結(jié)果Fig.4 The result of the quadratic differences of the surface temperature trend terms

4.2 LST異常面積的時間變化特征

圖5為2003年1月至2019年12月發(fā)生異常像素數(shù)量的時間變化圖. 從圖5、6、7可以看出，2008年3月21日地震前，1月17日約2個月前出現(xiàn)異常的面積覆蓋率達(dá)到19.5%，異常最明顯. 之后1月25日的異常面積覆蓋率為9.1%，2月2日的為16.6%. 2月2日之后異常面積覆蓋率依次減少.

圖5 震中附近LST時序熱異常提取結(jié)果Fig.5 Results of thermal infrared anomalies extracted of LST time series around the earthquake center

同樣，對于2014年2月12日地震，分析震前2個月，2013年12月份至2014年3月份研究區(qū)異常結(jié)果. 結(jié)果可知，異常較明顯的是2013年12月27日研究區(qū)1.4%范圍出現(xiàn)了異常，之后1月17日，1月25日異常覆蓋率為1.8%，1.0%. 發(fā)震前4天（2月10日）異常面積覆蓋率達(dá)到最大為2.1%.

4.3 LST異常的空間變化特征

以圖6，2008年3月21日地震前1月至3月LST異常時空演化過程可見，空間異常范圍主要出現(xiàn)在震中西北部. 隨著發(fā)震期將近，空間異常范圍擴(kuò)展到震中東北部，異常區(qū)域面積1月25日達(dá)到最大，出現(xiàn)異常的面積覆蓋率19.5%，2月18日起異常覆蓋率逐漸減少.

圖6 2008年于田Ms7.3地震熱異常時空演化Fig.6 Space-time evolution processes of thermal anomaly before and after the Yutian Ms7.3 earthquake in 2008

圖7 2014年于田Ms7.3地震熱異常時空演化Fig.7 Space-time evolution processes of thermal anomaly before and after the Yutian Ms7.3 earthquake in 2014

2014年2月12日相隔6年于田縣發(fā)生Ms7.3地震. 2013年12月至2014年2月地震熱異常時空演化圖如圖7，此次地震發(fā)生在冬季，比較2008年地震出現(xiàn)異常面積不明顯，異常范圍比較零散，震中西北以及東北部出現(xiàn)小范圍異常. 2月10日震前2天較小的異常區(qū)域出現(xiàn)在研究區(qū)東南部，且異常最明顯，異常面積覆蓋率達(dá)到2.1%. 發(fā)震之后，異常范圍隨著時間逐漸減少. 圖8表示研究區(qū)2003—2019年地震熱紅外異常的空間分布，2003年6月10日最為明顯，異常區(qū)域達(dá)到最大，熱異常方向主要表現(xiàn)為東南部分.

圖8 2003年6月10日的LST與熱異常空間分布圖Fig.8 Spatial distributions of LST and thermal anomaly on June 10，2003

5 討論與結(jié)論

本文以2003—2019 年連續(xù)17 年的MODIS LST 數(shù)據(jù)，采用基于STL 分解和GESD 檢驗(yàn)的異常檢測法，對于田縣兩次Ms7.3地震，地震熱異常演化過程進(jìn)行了分析，得到如下結(jié)論：

1）以于田縣地震為例，本論文中采用基于STL和GESD的異常檢測. 首先利用STL分解法將時間序列進(jìn)行分解，除了年變趨勢及其季節(jié)因素周期性的影響后，用GESD 檢驗(yàn)對STL 分解后的殘余項(xiàng)進(jìn)行異常檢測.研究結(jié)果可見，震前研究區(qū)LST有明顯的增溫現(xiàn)象. 此方法對地震熱異常的監(jiān)測提供了參考依據(jù).

2）對于2008年3月21日地震，異常區(qū)域面積1月25日達(dá)到最大，出現(xiàn)異常的面積覆蓋率19.5%，2月18日起異常覆蓋率逐漸減少. 對于2014年2月12日地震2月10日異常最明顯異常面積覆蓋率達(dá)到2.1%，發(fā)震之后，異常范圍隨著時間逐漸減少.

3）兩次地震異常區(qū)域主要靠近震中西北及東北部，但2008年地震異常比2014年異常較明顯，此研究區(qū)內(nèi)異常范圍較大.

總之，本文基于STL分解和GESD檢驗(yàn)的異常檢測算法提取了研究區(qū)地震熱紅外異常，探討了此方法對地震熱異常的監(jiān)測有一定的參考價值，并總結(jié)了地震演化過程. 此外，由圖可知，2003年6月10日異常達(dá)到最大，但此日前后沒發(fā)生過地震. 熱紅外異常與地震的關(guān)聯(lián)性可提高熱異常對于地震的檢測能力，這仍需要進(jìn)一步研究.