樂應(yīng)波， 陳福春， 陳桂林

1. 中國科學(xué)院上海技術(shù)物理研究所 紅外探測與成像技術(shù)重點(diǎn)實(shí)驗(yàn)室, 上海 200083;

2. 中國科學(xué)院大學(xué), 北京 100049

1 引 言

20 世紀(jì)80 年代，蘇聯(lián)學(xué)者分析熱紅外圖像發(fā)現(xiàn)了震前熱紅外輻射異常（Gornyi 等，1988）。國內(nèi)外學(xué)者對(duì)此進(jìn)行了大量研究，包括機(jī)理研究、特征分析和算法等（申旭輝 等，2018）。強(qiáng)祖基等（1997）和Qiang 等（1998）提出了震前溫度異常升高機(jī)理的初步模型：在高應(yīng)力作用下，沿地震斷裂帶釋放CH4、CO2和帶電粒子進(jìn)入大氣，在太陽輻射和電場激發(fā)的作用下大氣溫度升高。這就是震前紅外異常的理論模型。

隨著衛(wèi)星遙感技術(shù)和現(xiàn)代信息技術(shù)的發(fā)展，海量的熱紅外遙感數(shù)據(jù)和先進(jìn)的信號(hào)處理方法被應(yīng)用于提取震前紅外異常信號(hào)（Bellaoui 等，2017；Jiao 等，2018；荊鳳 等，2018；Ouzounov 等，2006；屈春燕 等，2006；張?jiān)?等，2004）。目前，常用于地震紅外遙感研究的數(shù)據(jù)主要來自美國發(fā)射的NOAA 系列衛(wèi)星 （Ouzounov 等，2007）、EOS 系列衛(wèi)星（Zoran 等，2014）和中國發(fā)射的風(fēng)云系列氣象衛(wèi)星（Yao 和Qiang，2012）。靜止軌道衛(wèi)星可以大范圍地連續(xù)觀測地表和大氣參數(shù)，在時(shí)間和空間上更具有可比性（Bormann 等，2003；Turk 等，2000）。Choudhury 等（2006）利用NOAA-AVHRR數(shù)據(jù)對(duì)伊朗的幾次地震進(jìn)行研究，發(fā)現(xiàn)該地區(qū)震前普遍存在大幅度的紅外增溫。

受到季節(jié)變化、天氣條件、地質(zhì)活動(dòng)和人類活動(dòng)的共同影響，無法直接觀測到地質(zhì)活動(dòng)引起的紅外增溫，可靠的異常識(shí)別和分析算法至關(guān)重要（Shen 等，2013；Rivera 等，2012）。常見的震前紅外異常提取算法包括RST （Robust Satellite Techniques） 算法（Tramutoli 等，2013）、四分位法、小波變換法、卡爾曼濾波法（Saradjian 和Akhoondzadeh，2011）、自回歸整合移動(dòng)平均法、人工神經(jīng)網(wǎng)絡(luò)和支持向量機(jī)（Akhoondzade，2013）和一些集成算法（Akhoondzade，2014）等。這些算法在個(gè)別典型震例中得到了驗(yàn)證，缺乏大量震例統(tǒng)計(jì)驗(yàn)證。在震前異常統(tǒng)計(jì)中，陽性預(yù)測值是異?，F(xiàn)象中地震前兆的比例，用于表征地震預(yù)測的準(zhǔn)確性；真陽性率是地震中有前兆的比例，用于表征震前異常的普遍性（Filizzola 等，2022；Zhang 等，2021；Fu 等，2020；Jiao和Shan，2021）。Lu 等（2016）利用FY-2E 衛(wèi)星數(shù)據(jù)對(duì)西藏地區(qū)多年地震熱紅外異常進(jìn)行了研究，結(jié)果發(fā)現(xiàn)5級(jí)以上地震震前普遍存在異常，但因未對(duì)所有異常信號(hào)進(jìn)行統(tǒng)計(jì)，尚無法預(yù)測準(zhǔn)確率；Zhang 和Meng（2019）對(duì)四川地區(qū)的熱異常進(jìn)行了統(tǒng)計(jì)分析，結(jié)果表明熱紅外異常與5 級(jí)以下的地震相關(guān)性不明顯，與5級(jí)以上的地震有較低的陽性預(yù)測值和真陽性率，預(yù)測潛力有限，反映異常信號(hào)識(shí)別算法效果不太理想。

功率譜相對(duì)變化法是震例研究中常用的紅外遙感地震信息提取算法（郭曉 等，2010；孟慶巖等，2017）。在典型震例的研究中，由于樣本數(shù)量少，尚不能從統(tǒng)計(jì)學(xué)的角度驗(yàn)證算法的準(zhǔn)確性和普適性；個(gè)別地區(qū)的統(tǒng)計(jì)分析中，存在研究區(qū)域小，無法觀測到大面積異常現(xiàn)象的整體空間分布情況的問題（Zhang 等，2019；Fu 等，2020）。為此，本文提出一種基于連通域的方法對(duì)中國及周邊地區(qū)的功率譜異常信號(hào)進(jìn)行統(tǒng)計(jì)，從統(tǒng)計(jì)學(xué)的角度評(píng)估算法的準(zhǔn)確性和普適性。該方法的核心思想就是將時(shí)間和空間上連續(xù)的異常體視為一個(gè)樣本，計(jì)算不同參數(shù)條件下異常信號(hào)的陽性預(yù)測值和地震的真陽性率，并分析異常信號(hào)的幅度、空間范圍特征和相關(guān)的地震信息。

2 數(shù)據(jù)和方法

2.1 功率譜相對(duì)變化法

本文所用的數(shù)據(jù)來自FY-2G 衛(wèi)星搭載的多通道掃描輻射計(jì)，所用波段為長波紅外波段（10.3—11.3 μm）。地表除了發(fā)出自身輻射外，還反射其他物體的輻射，為了避免太陽輻射的影響，選取午夜北京時(shí)間00：00—04：00 的亮度溫度進(jìn)行分析（張?jiān)?等，2011）。原始數(shù)據(jù)經(jīng)幾何校正和輻射定標(biāo)得到每個(gè)像素點(diǎn)對(duì)應(yīng)的經(jīng)緯度和亮度溫度數(shù)據(jù)。熱紅外輻射難以穿透云層，并且云頂?shù)臏囟冗h(yuǎn)低于地面，在云層覆蓋區(qū)域會(huì)形成孤立的低溫區(qū)（Kato 等，2011）。本文中剔除低于平均值且差值超過1.5 倍標(biāo)準(zhǔn)差的數(shù)據(jù)，然后求每日的平均值作為當(dāng)日值，時(shí)間序列中缺失的數(shù)據(jù)按近鄰法進(jìn)行插值，得到日亮度溫度序列。

在日亮度溫度數(shù)據(jù)中，地球溫度場的年變化是低頻信息，天氣變化和人類活動(dòng)引起的溫度變化是高頻信息（Ma 等， 2010）。小波分解可以分離信號(hào)中的高頻信息和低頻信息（Wang 等，2018；Wei 等，2013）。本文以7 階小波分解的低頻部分作為背景場，2 階小波分解的高頻部分可代表天氣變化信息等，低頻部分包含地震異常信息。用2 階小波分解的低頻信息減去7 階小波分解的低頻信息，可以在一定程度上消除季節(jié)和天氣的影響，又包含有可能與地震相關(guān)的異常信號(hào)（張?jiān)?等，2010）；最后，計(jì)算功率譜的相對(duì)變化值，從相對(duì)變化的時(shí)頻數(shù)據(jù)中選取極大值作為特征頻率，再尋找對(duì)應(yīng)的時(shí)間序列。從而獲得涉及到整個(gè)研究區(qū)域溫度信息功率譜的時(shí)空分布數(shù)據(jù)。

2.2 異常信號(hào)識(shí)別

異常信號(hào)就是功率譜相對(duì)變化的時(shí)空分布數(shù)據(jù)中時(shí)間上和空間上連續(xù)的高幅值點(diǎn)集（Zhang等，2017）。每個(gè)點(diǎn)集被視為一個(gè)異常信號(hào)樣本。

具體的算法流程如圖1所示，包括二值化、形態(tài)學(xué)處理和連通域識(shí)別4個(gè)步驟。

（1）二值化：設(shè)定幅值閾值為A0，第d天第m行第n列的功率譜相對(duì)變化的幅值為P（d，m，n）。當(dāng)P（d，m，n）≥A0，則認(rèn)為該點(diǎn)為高幅值點(diǎn)，即δ（d，m，n）=1；當(dāng)P（d，m，n）

（2）形態(tài)學(xué)處理：在數(shù)字圖像處理中，閉運(yùn)算是先膨脹后腐蝕的操作，用于填充空隙，連接相鄰區(qū)域，在本文中用于消除高幅值區(qū)域內(nèi)的小區(qū)域空隙，對(duì)比圖1（b）和圖1（c），閉運(yùn)算增強(qiáng)了異常點(diǎn)之間的空間連續(xù)性，避免將相同因素引起的異常判斷為多個(gè)異常。開運(yùn)算是先腐蝕后膨脹的操作，可以消除小的離散點(diǎn)（Sundararajan，2017），在纖細(xì)部位分離兩個(gè)區(qū)域，在本文中用于消除小區(qū)域的高幅值點(diǎn)，對(duì)比圖1（c）和圖1（d），開運(yùn)算消除了部分小面積異常區(qū)域，避免將其識(shí)別為異常信號(hào)。

圖1 功率譜異常信號(hào)識(shí)別算法Fig. 1 Recognition algorithm of abnormal signals in power spectrum

（3）連通域識(shí)別：對(duì)于一個(gè)區(qū)域D，若D中任意兩個(gè)點(diǎn)A1（d1，m1，n1）和A2（d2，m2，n2）都可以用完全屬于D的一條折線連接起來，則稱區(qū)域D是連通區(qū)域。本文通過遍歷和迭代的方法，識(shí)別所有高幅值點(diǎn)中的連通域。識(shí)別結(jié)果如圖1（e）所示，圖像中不同顏色表示不同的連通域，面積最大的連通域適合用于計(jì)算異常區(qū)域的位置和面積。

（4）樣本選擇：大氣氣團(tuán)的覆蓋面積為數(shù)十至數(shù)千公里，人類活動(dòng)的影響范圍約為數(shù)十公里（Song 等，2018）。Dobrovolsky 等（1979）依據(jù)理論模型得到孕震區(qū)域半徑R與地震震級(jí)M之間滿足R=100.43M。本文將覆蓋面積較大的連通域標(biāo)記為異常信號(hào)樣本，剔除人為因素造成的小范圍升溫。如圖1（f）所示，僅面積最大的連通域被標(biāo)記為異常信號(hào)樣本并用于下一步的統(tǒng)計(jì)分析。

2.3 統(tǒng)計(jì)分析

若在異常發(fā)生后t0天內(nèi)，有5 級(jí)及以上的地震發(fā)生，且震中與異常區(qū)域的距離小于d0，則該異常信號(hào)可視為震前征兆，屬于正樣本，否則屬于負(fù)樣本。其中，t0為時(shí)間閾值，d0為空間閾值，反映地震預(yù)測的空間精度和時(shí)間精度。陽性預(yù)測值是被判斷為正樣本的樣本中實(shí)際分類也為正樣本的比例，反映地震預(yù)測的準(zhǔn)確性；真陽性率指的是震前識(shí)別出異常的地震占所有地震的比例，反映算法預(yù)測地震的普適性（Zhang和Meng，2019）。陽性預(yù)測值（PPV）和真陽性率（TPR）如式（1）、（2）所示。其中，NA為異常信號(hào)樣本總數(shù)，NP為正樣本數(shù)量，NE為地震總數(shù)，N0為有震前異常的地震數(shù)量。

Molchan 圖表法是經(jīng)驗(yàn)性地震預(yù)測中常用的顯著性檢驗(yàn)方法（Zechar 和Jordan，2008）。圖表的橫坐標(biāo)為異常時(shí)空占有率τ，即預(yù)測有震區(qū)域在研究時(shí)空范圍內(nèi)的占比；縱坐標(biāo)為地震漏報(bào)率ν，即1 -TPR。其概率增益Gain如式（3）所示（蔣長勝 等，2011）。

為了研究震前異常的特征，本文對(duì)異常信號(hào)的空間分布參數(shù)進(jìn)行了統(tǒng)計(jì)。異常信號(hào)樣本中的峰值可以表征異常的幅度，幅度越大，則異常越明顯（馬未宇 等，2018）。由于異常區(qū)域并非規(guī)則的幾何圖形，本文用異常區(qū)域的最小外接矩形來描述異常的空間范圍，如圖2所示，最小外接矩形的4個(gè)頂角經(jīng)幾何校正后發(fā)生形變，用四邊形的最長邊長度來表示異常區(qū)域的長度。對(duì)正負(fù)樣本的異常幅度和長度進(jìn)行統(tǒng)計(jì)，可以直觀地看出震前異常與其它異常的差異。震前異常與震級(jí)相關(guān)，也受到周邊地形的影響（Jiao 和Shan，2022）。因此，本文對(duì)地震的時(shí)間、地點(diǎn)和震級(jí)進(jìn)行了統(tǒng)計(jì)分析。

圖2 異常區(qū)域長度Fig. 2 The length of abnormal region

3 研究區(qū)應(yīng)用實(shí)例與結(jié)果分析

3.1 研究區(qū)選擇與相關(guān)信息概述

選取2018 年功率譜相對(duì)變化的時(shí)空分布數(shù)據(jù)進(jìn)行研究，該數(shù)據(jù)中共包含365幅圖像，每幅圖像包含1192819個(gè)像素點(diǎn)。整個(gè)數(shù)據(jù)集中的大部分幅值集中分布于0—5。

本文研究的區(qū)域?yàn)?°—60°N，70°E—140°E，包括了整個(gè)中國及其周邊地區(qū)。在2018年2月1日至2019 年1 月1 日期間，該區(qū)域內(nèi)發(fā)生5 級(jí)及以上地震共47 次，無7 級(jí)以上地震。震級(jí)最大的是2018 年12 月29 日11：39：08 在棉蘭老島附近海域發(fā)生的6.9級(jí)地震，震中位置為5.85°N，126.89°E，震源深度為50 km，地震前的相對(duì)功率譜圖像如圖3 所 示。地震前8 天，即2018 年12 月21 日，震中周邊開始出現(xiàn)離散的高幅值點(diǎn)，范圍逐漸擴(kuò)大，強(qiáng)度逐漸增強(qiáng)；24 日到25 日期間，震中東北方向明顯形成孤立的高溫區(qū)域，然后逐漸消散；異常共持續(xù)了7天，在震前一天完全消失；震中距異常邊緣約231 km，距異常峰值約683 km。根據(jù)Dobrovolsky（1979）的理論模型，5級(jí)以上的地震，受影響區(qū)域的半徑應(yīng)大于141.25 km，7 級(jí)以上的地震影響區(qū)域的半徑應(yīng)大于1023.29 km，此次異常出現(xiàn)的位置，在地震的影響范圍內(nèi)。本文中所研究的地震，其影響區(qū)域半徑為數(shù)百公里。Tronin（2006）通過不同地區(qū)的衛(wèi)星數(shù)據(jù)研究，發(fā)現(xiàn)地震前6—24天出現(xiàn)熱異常，地震后持續(xù)約一周，且異常對(duì)大于4.5 級(jí)的地震敏感。因此，本文取幅值閾值A(chǔ)0分別為5、10、15、20、25，距離閾值d0分別為200 km、400 km、600 km、800 km、1000 km，時(shí)間閾值t0為30 d，陽性預(yù)測值和真陽性率如表1所示。A0越小，d0越大，則陽性預(yù)測值和真陽性率越大。約66%的地震在震前一個(gè)月內(nèi)，出現(xiàn)幅值大于5的異常信號(hào)，且震中與異常區(qū)域的最短距離不超過400 km。陽性預(yù)測值較低，表明有大量的異常信號(hào)與短期內(nèi)的地震無關(guān)。

圖3 震前相對(duì)功率譜圖像Fig. 3 Power spectrum image before the earthquake

本文通過Molchan 圖表法 （蔣長勝 等，2011）進(jìn)行統(tǒng)計(jì)檢驗(yàn)，A0取15、20、25 的真陽性率偏低。A0取5 和10的統(tǒng)計(jì)檢驗(yàn)結(jié)果如圖4所示，概率增益最大為1.76，對(duì)應(yīng)的閾值選擇為A0=5，d0=400 km，對(duì)應(yīng)的異常時(shí)空占有率為0.3753，漏報(bào)率為0.3403。通過表1 可知，該閾值條件對(duì)應(yīng)的陽性預(yù)測值為0.2037（20.37%），真陽性率為0.6596（65.96%）。

表1 陽性預(yù)測值和真陽性率統(tǒng)計(jì)結(jié)果Table 1 The statistical results of positive predictive value and true positive rate

圖4 Molchan圖表法進(jìn)行統(tǒng)計(jì)檢驗(yàn)Fig. 4 Statistic test using Molchan diagram method

3.2 異常信號(hào)特征分析

為達(dá)到較高的空間精度和真陽性率，本文取A0=5，d0=400 km，t0=30 d。正負(fù)樣本的峰值分布如圖5 所示，長度分布如圖6 所示，橫坐標(biāo)為異常信號(hào)的峰值和長度，縱坐標(biāo)為峰值或長度在給定范圍內(nèi)的樣本頻率。正樣本集中高幅度、大范圍的異常信號(hào)所占的比例高于負(fù)樣本集。異常樣本的峰值和長度如圖7 所示。峰值和長度小的樣本中，絕大部分是負(fù)樣本。峰值高于10，長度大于500 km 的異常信號(hào)樣本，可達(dá)到33.71%的陽性預(yù)測值。峰值高于20，長度大于2000 km的異常信號(hào)樣本，可達(dá)到80.00%的陽性預(yù)測值。

圖5 異常信號(hào)峰值分布Fig. 5 The frequency distribution of abnormal peaks

圖6 異常區(qū)域長度分布Fig. 6 The frequency distribution of the length of abnormal region

圖7 異常的峰值和長度Fig. 7 The peak and length of anomalies

3.3 地震信息統(tǒng)計(jì)

在本文統(tǒng)計(jì)的47 次地震中，有25 次地震震級(jí)不超過5.4 級(jí)，不同震級(jí)的震前異常統(tǒng)計(jì)如圖8 所示，橫坐標(biāo)為震級(jí)，縱坐標(biāo)為對(duì)應(yīng)震級(jí)中有震前異常和無震前異常的地震數(shù)量。5.4 級(jí)及以下的地震，有無異常的地震數(shù)量接近。22次大于5.4 級(jí)的地震中，僅有4 次地震在震前一個(gè)月內(nèi)未發(fā)現(xiàn)異常，真陽性率可以達(dá)到81.82%。圖9 為地震發(fā)生與異常消失的時(shí)間間隔分布圖，橫坐標(biāo)為地震與異常結(jié)束的日期間隔，正數(shù)表示地震發(fā)生在異常消失后，負(fù)數(shù)表示地震發(fā)生在異常消失前，研究表明，大部分地震發(fā)生在異常消失后。圖10 為震中的地理位置分布，日本、中國的臺(tái)灣省、菲律賓和印度尼西亞位于環(huán)太平洋地震帶上，這幾個(gè)國家和地區(qū)在研究區(qū)域和時(shí)段內(nèi)共發(fā)生25次地震，僅有4次地震在震前一個(gè)月內(nèi)未發(fā)現(xiàn)異常。新疆維吾爾自治區(qū)、西藏自治區(qū)、青海、四川和云南位于地中海-喜馬拉雅地震帶上，該地震帶上地震的真陽性率明顯低于環(huán)太平洋地震帶。西藏自治區(qū)和青海位于中國地殼厚度最厚的地區(qū)，在該地區(qū)的5次地震中僅有2次在震前一個(gè)月內(nèi)發(fā)現(xiàn)了異常。

圖8 不同震級(jí)的震前異常統(tǒng)計(jì)Fig. 8 Statistics of anomalies before earthquakes with different magnitude

圖9 地震與異常消失的時(shí)間間隔分布Fig. 9 Time intervals distribution between the earthquake and the disappearance of the anomaly

圖10 震中位置分布Fig. 10 Distribution of epicenters

4 結(jié) 論

本文通過功率譜相對(duì)變化方法提取紅外遙感圖像中的異常信號(hào)，并基于連通域方法對(duì)異常信號(hào)進(jìn)行識(shí)別和統(tǒng)計(jì)分析。研究空間范圍大，可以比較完整地觀測大面積紅外異?，F(xiàn)象。本文在不同參數(shù)條件下進(jìn)行統(tǒng)計(jì)分析和顯著性檢驗(yàn)，在A0=5，d0=400 km 的條件下在得到20.37%的陽性預(yù)測值和65.96%的真陽性率，其概率增益為1.76。本文的結(jié)果驗(yàn)證了功率譜相對(duì)變化法能在大部分地震前提取到紅外遙感圖像異常，但陽性預(yù)測值偏低，預(yù)測的準(zhǔn)確率有待提升。本文對(duì)研究區(qū)域異常信號(hào)的特征參數(shù)和地震信息進(jìn)行統(tǒng)計(jì)分析，得到以下結(jié)論：

（1）功率譜相對(duì)變化法提取出的高幅度、大范圍的異常信號(hào)（峰值大于20，長度大于2000 km），可以達(dá)到80.00%的陽性預(yù)測值，說明剔除小幅度小范圍的異常信號(hào)可以在一定程度上提高地震預(yù)測的陽性預(yù)測值；

（2）對(duì)于大于5.4 級(jí)的地震，真陽性率可以達(dá)到81.82%，說明本文的方法對(duì)大于5.4 級(jí)的地震敏感；

（3）真陽性率具有明顯的地區(qū)差異，環(huán)太平洋地震帶的地震真陽性率要高于地中海-喜馬拉雅地震帶。

本文提出的方法適用于大范圍、長期的震前異常統(tǒng)計(jì)和特征分析，解決了樣本數(shù)據(jù)少、研究區(qū)域小的問題。在此基礎(chǔ)上，可以進(jìn)行更深入的研究，包括但不限于以下幾個(gè)方面：

（1）增加數(shù)據(jù)樣本，分類統(tǒng)計(jì)，研究不同地震類型的震前異常情況。

（2）提取更多的異常信號(hào)特征，分析正負(fù)樣本的差異，提高陽性預(yù)測值。

（3）結(jié)合多種數(shù)據(jù)源和數(shù)據(jù)處理方法，對(duì)比不同算法的準(zhǔn)確性和普適性，探索更可靠的震前征兆。

志 謝本文所用部分靜止氣象衛(wèi)星數(shù)據(jù)資料來自中國氣象局國家衛(wèi)星氣象中心，感謝郭強(qiáng)研究員所提供的幫助和指導(dǎo)。