利用均線差值振幅波動水平分析汶川、雅安蘆山、九寨溝地震前熱紅外異常＊

廖洪月 董 娜 王 剛

1) 西安市地震局，西安 710007

2) 西安市地震監測中心，西安 710007

引言

在地震監測預報中，衛星熱紅外觀測資料與其他地震前兆資料相比，其突出的優點是覆蓋面廣，平面分辨率高，時間連續好，且不受地理環境限制[1]。在多山的四川省，于2008年、2013年和2017年分別發生了汶川MS8.0地震（31°N，103.4°E）、雅安蘆山MS7.0地震（30.3°N，103°E）和九寨溝MS7.0地震（33.2°N，103.8°E）。地震發生后，國內學者相繼進行了與熱紅外相關的研究并取得了一定進展。目前，熱紅外異常提取方法主要包括4大類：①目視解譯；②基于差值分析的異常提取算法；③基于信號分析的異常提取算法；④基于背景場分析的異常提取算法[2]。本文主要運用第②類方法，較為成功的提取出了汶川、雅安蘆山、九寨溝地震前熱紅外異常的時空序列。

當前學術界主要研究方法為小波功率譜分析法、距平法等，主要側重于短時間內絕對亮溫升溫異常的研究。本文的差值振幅分析法則側重于中短期相對升溫異常的研究，且在不需要剔除錯誤數據和去云干擾影響等人工干預手段的情況下，提取的震前震后異常時空特征較為明顯，異常時間和異常區域的研究結果與有關文獻有相同之處（在后面具體震例分析中有描述）。

1 數據選取

本文熱紅外數據是由四川省地震局提供的TERRA衛星MODIS遙感傳感器亮溫數據。數據時間范圍為2004年1月1日—2020年12月31日，數據空間區域范圍為（30°N—45°N，100°E—130°E），數據分辨率為1 km，由450萬個像元組成，單個像元每天一組數據，數據覆蓋中國大陸西南、西北、華北、華東等部分地區。

在數據分析過程中，將分辨率為1 km數據分別反演成10 km、20 km、30 km、50 km和100 km分辨率數據后進行相關處理分析。對比發現，數據分辨率并非越高越好，雖然不同分辨率數據的時間曲線、空間圖形態總體相似，但是分辨率較低時，空間圖相關指數等值線平滑度更佳，且分辨率越低后續數據處理效率越高。但分辨率過低，存在影響異常空間定位精度和異常峰值時間滯后等缺點。本文以100 km分辨率熱紅外數據作為研究的基礎數據。

2 數據處理及研究方法

前人利用不同的衛星遙感熱紅外數據，采用不同方法研究大震前熱紅外異常的結果顯示，獲取的異常信息幾乎都是升溫異常[3]。本文研究也是以升溫異常為研究對象，不同之處在于側重研究升溫過程中均線差值振幅的變化水平異常，而非絕對的亮溫升溫異常。

2.1 數據處理

本文的數據處理流程如下：

（1）數據反演：將空間分辨率為1 km數據分別反演為10 km、20 km、30 km、50 km和100 km數據。以100 km分辨率為例，即將100 km×100 km區域內10000組1 km分辨率網格數據進行均值化處理，轉化為一組數據，將450萬個1 km網格單元轉換為450個100 km網格單元。

（2）滑動均線計算：以365天作為年周期對反演后的分辨率為100 km亮溫數據進行滑動均線計算。將數據區域內450個網格單元每天的亮溫數據轉化為年均線數據。

（3）差值分析：對每個單元網格年均線差值的振幅變化特征進行分析。

（4）產出異常區域單元網格時間曲線、區域空間異常疊加圖和剖面圖等。

實現流程（3）和流程（4）相關功能的計算公式為：

式（1）中，Ai和Ai?t分 別為第i天、第（i?t）天的亮溫年均線值。在研究中發現，數據分辨率越低，代表網格區域空間范圍越大，能觀察到區域內的亮溫異常現象就可能越多，觀察到異常的時間也可能越長。當用均線分析異常時，均線周期應大于熱紅外異常發展周期，這樣的均線值就可能包含整個異常發展周期的信息。為避免均線值受到季節變化的影響，本文特將均線周期設定為365天，默認為1年，在理想情況下，年均線任何時間點的均線值相等。當出現亮溫升溫異常時，年均線值會伴隨亮溫升溫發展過程不斷增大。△Bi是紅外亮溫均線值Ai的 差值，t為差值分析的時間間隔，可根據需要設置。文獻[4-6]表明，熱紅外異常通常只持續10天到幾十天之間。當前均值與t天前均線值進行差值對比時，假設t天前處于正常狀態。t值不宜過小，當t值過小時，則可能出現異常值與異常值間的差值比較，因異常對沖而導致提取的異常減弱或異常被屏蔽現象；t值也不宜過大，t值過大時可能導致本次周期異常與上一個周期異常進行差值對比，也可能出現異常屏蔽現象。本研究經過反復測試，認為t值約為熱紅外異常的發展周期為最佳，但判斷異常發展周期非常復雜。不同的t值，產生的差值系列不一致，最后提取的異常結果也有所不同（表1）。本文研究上述3次地震時，t值統一設定為20天。式（2）中，△Bmin為數據起始時間到當前時間段內的最小差值，可能隨時間變化而變化，是一個變量，△Ti是當前亮溫均線差值與最小均線差值的差值，也就是當前差值的振幅變化值，振幅變化值越小，表示熱紅外活躍性越低，振幅變化值越大，表示活躍性越高。式（3）中，是當前的平均值，n為數據起始時間到當前的總天數，隨著時間變化而變化，也是一個變量，值高意味著對應時段內熱紅外活躍性高。本研究認為，n的最大值也就是選擇的數據時段不宜過長，要避免數據時段內存在多個熱紅外異常時段，包括已經對應發生地震的異常和其他離當前時間較長的熱紅外異常，數據時段內只存在一個異常時段最為理想。在本文中，選擇連續數據時長一般不超過2年，不小于1年。式（4）中，用于計算均線差值振幅變化值與平均變化值的變化比，比值為增大倍數，將其作為亮溫異常參數，其基準值為零，大于零時表示處于升溫異常狀態，小于零時表示升溫活躍性水平下降，其活躍性低于歷史平均水平，也可理解為處于降溫異常狀態。在研究中，為突顯異常觀測效果，在不改變異常參數正負屬性條件下，對增大倍數進行3次冪運算放大，當存在異常時，放大后的冪值在時間曲線中有一目了然的視覺效果。空間剖面圖等值線Yi則為原始值。

表1 地震熱紅外異常特征相似點對比Table 1 Comparison of the thermal infrared anomaly similarity

2.2 異常分析原理

異常分析原理和異常判斷標準：①時間異常分析。理論上均線差值應該是一個穩定值，包括兩種情形：當均線為水平線時，差值應該為0；當均線為一定斜率的直線時，差值為一個穩定值，此時差值振幅變化值應為0。當均線不穩定時，差值數據系列的振幅也發生變化。振幅變化越大，出現異常的可能性越大。本文設定的異常判別條件是：當振幅變化大于平均值的2倍時，則判斷為增溫異常，值越大，表示時間異常越高。②空間異常分析。當空間平面某區域的異常值高于2且在整個空間中處于相對高位時，認為該區域存在空間異常。當時間和空間異常同時存在且相互對應時，將該現象判定為震前異常。

3 結果分析

3.1 汶川8.0級地震

以汶川地震為例，實際工作中快速異常分析分為3步：①通過時空掃描繪制2006年8月1日—2008年5月12日即汶川地震前這段時間的異常疊加圖，通過異常疊加圖找出異常區域；②提取異常核心區時間曲線，分析該區域震前熱紅外活動特征，判斷異常發生時間等；③以5天為步長，繪制連續的異常平面圖，分析震前與震后熱紅外活動的空間變化特征。

本文的異常疊加圖是指一段時間內產生的異常在同一平面上的投影，它能顯示異常發生的位置但不能顯示異常發生的時間。本文3個震例分析時選取空間的范圍為（30°N—45°N，100°E—115°E）。

由2006年8月1日—2008年5月12日時間段內時間掃描異常疊加圖發現，以（37°N，102°E）為核心的較大區域出現異常（圖1）。異常核心區網格單元的時間曲線顯示異常峰值發生于3月18日，距離汶川地震55天，異常放大值約為55，折算值≈3.8（圖2），即最大峰值倍數為3.8倍，表示當天差值振幅變化水平高于歷史平均水平3.8倍。從異常疊加圖可以看出，異常區域在震中的北部偏西，地理上異常主體區域位于柴達木塊體，該結果與路茜等[3]的研究有相同之處，異常發生時段的結果與Xie等[4]的研究有相同之處。

圖1 2006年8月1日—2008年5月12日異常疊加圖（紅圓為震中）Fig.1 Abnormal stacking diagram from August 1，2006 to May 12，2008（the red circle is the epicenter）

圖2 異常區均線差值振幅隨時間變化關系（縱軸為異常放大系數）Fig.2 The relation between the amplitude of mean difference in the abnormal area and time（the vertical axis is the abnormal amplification index）

對步長為5天的連續空間剖面分析發現：從2008年1月下旬開始，汶川地震中北部區域出現大范圍的低值異常，表明該區域熱紅外活躍性降低，該狀態持續約1個多月。從3月初開始，震中北部區域異常值開始緩慢上升，3月中旬高值異常范圍達到最大，之后開始減弱。至汶川地震時，異常值較低，震后約20天之后的整個6月份，整個數據區域內大范圍呈現低值異常狀態，說明這個時段熱紅外活躍性水平非常弱（圖3）。

圖3 汶川地震前后熱紅外活躍性水平時空變化過程Fig.3 Temporal and spatial variation of thermal infrared activity level before and after Wenchuan earthquake

3.2 雅安蘆山7.0級地震

雅安蘆山地震的分析方法步驟與汶川地震相同。2012年1月1月—2013年4月20日異常疊加圖顯示存在兩個高異常區域，一個異常位于（40°N，112°E）區域（圖4右上），異常區呈短條帶狀。該區域時間曲線顯示于2012年5月出現過一次高值異常，經查詢發現同年5月9日和6月18日在山西與陜西交界處各發生一次4.3級地震（數據來源于USGS），本文不作詳細描述。另一個異常位于雅安蘆山震中的西北方向（圖4左下），由于受熱紅外原始數據空間范圍局限，只能顯示部分異常區域而不能看到全貌。

圖4 2012年1月1月—2013年4月20日異常疊加圖（紅圓為震中）Fig.4 Abnormal stacking diagram from January 1，2012 to April 20，2013（the red circle is the epicenter）

異常核心區域時間曲線顯示，異常峰值發生于3月16日，距雅安蘆山地震35天，最大幅度變化高于歷史平均水平3.1倍（圖5）。

圖5 異常區均線差值振幅隨時間變化關系（縱軸為異常放大指數）Fig.5 The relation between the amplitude of mean difference in the abnormal area and time（the vertical axis is the abnormal amplification index）

步長為5天的連續異常圖（圖6）顯示：從2013年2月上旬至2月下旬，震中的偏西區域，熱紅外活躍性水平低值異常持續了20多天，之后活躍性強勁上升，3月中旬達到峰值后迅速下降。至4月20日地震時，震中及附近區域無明顯異常。震后約20天之后，震中附近區域熱紅外活躍性水平開始下降，低活躍性持續了約20天，該異常現象與汶川地震后低值異常相似。本次地震前活躍性高值異常區域的研究結果與魏樂軍等[5]的研究結果相近。

圖6 雅安蘆山地震前后熱紅外活躍性水平時空變化過程Fig.6 Temporal and spatial variation of thermal infrared activity level before and after Lushan earthquake in Yaan

3.3 九寨溝7.0級地震

分析方法與汶川地震相同。2016年1月1日—2017年8月8日時空異常疊加圖顯示，以（35°N，101°E）為核心的較大區域存在高值異常（圖7），以（35°N，101°E）為中心的異常區域時間曲線顯示，分別在1月、5月和7月出現過3次較大異常，最大異常峰值發生于7月17日，距離九寨溝地震22天，當日最大振幅變化高于歷史平均水平3.1倍（圖8）。

圖7 2016年1月1日—2017年8月8日異常疊加圖（紅圓為震中）Fig.7 Anomalous stacking diagram from January 1，2016 to August 8，2017（the red circle is the epicenter）

圖8 異常區均線差值振幅隨時間變化關系（縱軸為異常放大系數）Fig.8 The relation between the amplitude of mean difference in the abnormal area and time（the vertical axis is the abnormal amplification index）

異常剖面圖變化過程顯示，從2017年4月中旬至5月初，震中西部區域有一個明顯的異常指數下降過程，說明該區域熱紅外活動性水平在這段時間呈下降狀態，之后活動性呈現寬幅的振蕩變化，分別于5月20日和7月20日在該區域產生高值異常后迅速下降，第2次異常強度和異常區域面積均大于第1次，8月8日發震后震中南部區域活躍性指數有下降現象（圖9）。

圖9 九寨溝地震前后熱紅外活躍性水平時空變化過程Fig.9 Temporal and spatial variation of thermal infrared activity level before and after Jiuzhaigou earthquake

本次地震的震前異常時間與異常區域研究結果與戴勇等[6]、楊星等[7]研究結果有相同之處，但異常區域研究結果與張麗峰等[8]的研究結果不一致。

3.4 亮溫異常的時空變化特征

不同的時間間隔，均線差值振幅變化結果有所不同。汶川地震、雅安蘆山地震、九寨溝地震前與震后相關異常特征簡要對比如表1。

3.5 低分辨率反演與年均線差值分析的潛在優點

在處理亮溫數據時，通常需要人工干預：剔除錯值、云干擾值及不符合黑體輻射公式的高值，計算剩余數據的均值，從而得到扣除部分云影響的亮溫日值[9]，本文數據處理過程中省略人工剔除錯誤數據和去云干擾工作等環節，也達到幾乎相同異常提取的效果，也就是說，本文方法具有一定的容錯功能和自動去云干擾影響功能，其原理如下：

容錯原理：將高分辨率由1 km×1 km數據反演為100 km×100 km數據后，單個錯誤值的影響稀釋到原來的1/10000，將分辨率為100 km×100 km網格數據進行年均線轉換后，可再將影響再稀釋為現有的1/365，總影響量降為最初的1/3650000，此時錯誤數據負面影響非常有限，說明低分辨率反演和年均值轉換處理具有非常理想的容錯功能。

自動去云干擾影響原理：我們可以簡單地認為，任意一組年均值由受云干擾影響的亮溫值和正常陸地亮溫值組成，同時也可以簡單地認為，每組年均值受云干擾影響的量值是相同的。這樣，通過差值計算，直接將雙方均擁有的云干擾影響去除。本文中時間相差20天的年均線數據包含了相差無多的云干擾影響，在差值計算過程中云干擾影響被自動全部對沖或絕大部分被對沖掉，因此，無需人工處理。本研究認為，當均線周期較短時，計算差值、時間間隔中可能存在季節差，尤其是存在陰雨季節差別時，不能忽略云干擾影響，仍有必要做去云干擾工作。但當均線周期為年或年的整數倍時，均線中每組數據中包含相同的季節，滑動差值計算時不存在季節差，因此，可以忽略云干擾影響。

由于不需人工干預，且不同的數據處理流程均有統一的算法，因而數據處理的全流程均可實現計算機全自動分析完成。這樣，不但提高了數據處理效率，同時也降低了數據處理和分析的技術難度，從技術角度為研究成果轉化和應用提供了便利條件。

4 結論與討論

4.1 結論

通過以上的系統分析，得出主要結論如下：

（1）本文熱紅外異常參數Yi的實質是變化率，其基準值為零，其異常活躍性水平是針對升溫異常而言，其數值的大小僅代表升溫異常活躍性水平的高低，其內在的意義如下：異常值越趨近于零，表示熱紅外活動越穩定；當異常值大于零時，異常值越高表示熱紅外升溫異常活躍性水平越強；當異常值小于零時，異常值越低則表示熱紅外降溫異常活躍性水平越強。

（2）汶川、雅安蘆山和九寨溝地震前均出現熱紅外亮溫均線差值振幅變化正異常，也就是活躍性上升異常，且地震均在峰值后較長時間發生，該研究結果具有地震短臨預測的意義。

（3）上述3次地震中，汶川地震級別最高，其震前異常值大于2倍的異常區域也最大。3次地震震中均在異常區域外圍而不在異常優勢區域，震中均位于異常區東南方向。

（4）上述3次地震前伴隨有活躍性低值異常和活躍性高值異常。先有低活躍性異常，再演變為高活躍性異常，低值時低至歷史平均水平50%以下，且低值狀態存續時間較長，高值時高于平均水平3倍以上。震前熱紅外存在兩種性質相反的異常現象，在文獻中較為罕見。

（5）本方法在提取異常過程中雖然缺少了剔除錯誤數據和去云干擾影響兩個看似不可或缺的環節，但仍能提取到與前人相近或相同的研究結果，說明本方法可行性較高。

（6）通過低分辨率反演原始熱紅外數據和對反演后日均值數據實行年周期均值轉換，有效地稀釋了原始數據中可能存在的錯誤數據帶來的負面影響，使本方法具有一定的容錯功能。

（7）本方法巧妙利用年周期均線差值對沖功能，從而自動地實現了最大限度去云干擾影響。換句話說，通過均線年周期的設定，使本方法具有抗云干擾功能。

（8）本研究方法是對全時段、大區域的衛星原始數據進行系統的整體的分析，不對原始數據做任何的刪改，確保在百分之百真實的原始數據基礎上提取異常，所有數據參與處理的流程統一，分析方法統一，每組原始數據均同權參與全流程計算，最后得出的數據結果不含任何主觀因素。由于實現了高速高效的自動化數據處理流程，本文認為該方法適合大規模、大范圍數據甚至全球熱紅外數據同步或準同步處理，可用于全國性大區域熱紅外實時異常動態監測。

4.2 討論

在研究過程中，同時也發現了另一些現象，值得進一步討論和研究：

（1）與前人通常側重于絕對升溫異常不同，本文是以相對異常作為研究方向。亮溫升溫絕對高值異常與相對高值異常并不一定同步，絕對亮溫升溫高值異常時并不一定產生相對亮溫高值異常，同理，相對高值異常也不一定會同時出現絕對升溫高值異常，因此，本文關于異常區域和異常時間的研究結論與其他學者研究結論存在出入屬于正常現象。

（2）汶川、雅安蘆山、九寨溝3個地震前均存在相同的異常發展過程，先有活躍性水平大幅下降，后又轉為活躍性大幅上升，高值異常峰值過后發震。震前兩種性質相反的異常現象是否為大地震前共有現象以及其產生的原因，值得深入研究。

（3）對沖的作用與副作用。本文差值分析時應用對沖特性去云干擾影響，從而省略了去云干擾數據處理人工干預工作，但差值間隔時間設置不當可能會對沖掉部分真實異常，影響最終的異常提取效果。因此，差值間隔時間參數設置是一個值得深入研究的問題。

（4）3次地震震前異常值在3.1—3.8倍之間，以本文的2倍可作為判斷異常的參考指標是否恰當，需更多震例驗證和深入研究討論。

（5）汶川地震前，距震中500—600 km的大面積異常區域是不是該震前兆值得討論。本文認為，雖然異常區域距震中較遠，但從圖2的時間曲線分析，此區域在超過10年的時間內，僅在汶川地震前出現過一次這種罕見異常現象，因此，從時間和概率上對應分析，它是汶川地震前兆的可能性較大。根據馬瑾院士[10]的研究，一個地震的發生可對周邊不同構造區造成升溫或降溫影響，距離有遠有近，500—600 km在該文獻的震例中不算是遠距離。因此，這種現象在空間上認為是汶川地震前兆也是可能的。

（6）本研究方法的缺點，在特定條件下可能出現偽異常。假設熱紅外歷史數據非常平穩，則其年線差值振幅變化的歷史平均值趨近于零，此時如果熱紅外出現小幅異常變化，則可能出現數值特別高的異常值，因此，要慎重對待高值異常的計算結果。

本文采用與眾不同的熱紅外分析方法，發現了一些有趣的現象，這些現象是否可作為地震前兆還需進一步研究。由于本文所用到的數據范圍較小，震例僅3例，雖然本文研究結果顯示，震前異常較為明顯且異常峰值出現在地震發生前，但由于3個震例的震中均不在異常區內，因此，本研究方法用于實際地震預測還有很遠的距離，研究方法的科學性也有待更大范圍的數據累積和更多震例的驗證，希望獲得更多學者的批判指正。

致謝

作為一名從事行政工作的業余研究者，感謝西安市地震局領導對本人和本研究提供寬松的研究環境與相關支持，同時感謝同事田暉先生為本文圖片合成的辛苦付出。