謝 庭1a,1b,1c，陳 忠1a,1b,1c，李志平2，張寧新1a,1b,1c，郭莉莉1a,1b,1c

(1.華中科技大學 a.自動化學院；b.多譜信息處理技術國家級重點實驗室；c.圖像信息處理與智能控制教育部重點實驗室，武漢 430074；2.大港油田集團通信公司，天津 300280)

1 概述

到目前為止，地震預報，尤其是短臨地震預報，是一個世界性難題。地震預報必須同時包括時間、地點和強度。由于地震情況復雜，有些地震能預報，有些則無法預報，現(xiàn)在全球預報地震的準確率只有20%多。然而地震是極具破壞性的自然災害，給人類帶來了巨大的物質和精神損失，因此探索和發(fā)現(xiàn)預報地震的新方法、新途徑，在防震減災方面具有十分深遠的意義。

隨著高新技術的不斷發(fā)展，衛(wèi)星熱紅外遙感為地震預測預報提供了新的方法和手段。早在1971年，我國著名地震學家傅承義提出紅腫學說，為后來利用衛(wèi)星遙感、監(jiān)測地表與低空大氣的異常增溫預報地震奠定了理論基礎[1]。1988年，前蘇聯(lián)科學家Gorny、Salman和Tronin在研究中亞以及東地中海地區(qū)的地震時就發(fā)現(xiàn)在中強地震震前，震區(qū)及其周邊區(qū)域對應的衛(wèi)星紅外圖像上存在熱紅外異常[2]。意識到這一發(fā)現(xiàn)的重要性，我國一大批學者投入到了衛(wèi)星熱紅外異常預測地震的研究中。例如，強祖基、賃常恭等利用衛(wèi)星遙感熱紅外圖像，經(jīng)過多年對震前熱紅外異常進行分析和研究，取得了顯著的地震預報成果[3-5]；徐秀登、徐向民通過分析反生在中國及相鄰地區(qū)的40多次地震的熱紅外圖像，總結了衛(wèi)星紅外臨震異常的基本特征[6-7]；劉德富等研究了地球長波輻射時空變化，結果表明在地震發(fā)生前的震區(qū)呈現(xiàn)出比周圍區(qū)域更顯著的輻射增強變化特征[8]。還有其他很多學者也參與到了衛(wèi)星熱紅外異常預測地震的研究中，并取得了一定的成果。在上述傳統(tǒng)方法研究過程中，紅外異常均由人眼識別或手工計算得出，而當今時代人工智能發(fā)展日趨成熟，利用圖像識別技術提取紅外異常是這一研究領域發(fā)展的必然趨勢。

本文提出一種基于圖像識別熱紅外遙感異常數(shù)據(jù)的自動地震預測方法。不僅為基于熱紅外遙感數(shù)據(jù)的地震預測研究領域打開了人工智能的大門，而且實現(xiàn)國產(chǎn)風云2號氣象衛(wèi)星數(shù)據(jù)在地震預測中的應用。傳統(tǒng)方法使用的熱紅外數(shù)據(jù)取自國美NOAA衛(wèi)星或日本MTSAT衛(wèi)星，而本文方法使用國產(chǎn)風云2號氣象衛(wèi)星(FY_2C)1通道分區(qū)圖。該圖一般應用于天氣預報，也是首次將其應用于地震預測預報中。

2 震象云地震預測原理及圖像識別方法描述

2.1 震象云地震預測原理

通過對衛(wèi)星云圖的觀察，發(fā)現(xiàn)地震前震中上方的衛(wèi)星紅外云圖會出現(xiàn)增溫型特殊云團和降溫型特殊云團，一般接近發(fā)震時間時出現(xiàn)較頻繁，稱這種特殊云團為震象云。

震象云是地震前震區(qū)上空出現(xiàn)增溫和降溫異常的特殊云跡。圖1所示是在汶川地震震前成都上空出現(xiàn)的震象云。

圖1 汶川地震震前汶川附近出現(xiàn)的震象云

震象云與氣象中的云團是有差別的，震象云一般位于對流層頂，與周圍氣象云團有溫度差異，有一定的紋理特征，并且這些震象云移動速度要遠小于周圍氣象云團的移動速度，且接近發(fā)震時間時出現(xiàn)較為頻繁。因此，可以綜合震象云的顏色、紋理以及浮現(xiàn)頻率來預測地震。

2.2 圖像識別方法描述

從紅外云圖中觀察到增溫型震象云浮現(xiàn)頻率相對較為頻繁，顏色、紋理特征較降溫型震象云更為突出，因此，本文針對增溫型震象云進行自動檢測和識別。

紅外云圖中增溫型的震象云表現(xiàn)為在高層云中出現(xiàn)中、低層云，灰度級集中出現(xiàn)在110～145之間，云團大小在一定的范圍變化，形狀不確定，有一定模式的紋理特征，而且震象云是在地震地點上空附近浮動，可以根據(jù)以上特征對其進行目標識別和跟蹤。針對震象云的以上特點，對其采取了以下方案進行處理。首先對輸入數(shù)據(jù)進行預處理，拉伸震象云所在目標灰度區(qū)段的灰度級，壓縮非目標灰度區(qū)段的灰度級；利用灰度共生矩陣對熱紅外數(shù)據(jù)進行紋理特征提取，使用BP神經(jīng)網(wǎng)絡模型訓練目標神經(jīng)網(wǎng)絡，將紋理特征輸入目標神經(jīng)網(wǎng)絡進行識別，提取疑似目標；根據(jù)目標的大小范圍濾掉非目標；對剩下的疑似目標進行跟蹤，依據(jù)其浮現(xiàn)頻率精確定位目標，提取得到震象云。主要分為4個階段：(1)預處理階段；(2)目標識別階段；(3)目標過濾階段；(4)目標跟蹤階段。處理流程如圖2所示。

圖2 震象云提取整體流程

2.2.1 預處理

根據(jù)觀測，增溫型震象云主要位于是中層云區(qū)段，在紅外云圖中灰度級集中于110～145之間。為了在后續(xù)流程中能夠充分提取增溫型震象云的紋理特征，必須擴大中層云灰度級的浮動范圍。因此，要壓縮低層云和高層云，拉伸中層云。本文采用分段線性拉伸方法增強震象云區(qū)域對比度，采用的公式如下所示：

其中，r是輸入灰度級；s為輸出灰度級。

2.2.2 目標識別

由于BP網(wǎng)絡能學習和存儲大量的輸入-輸出模式映射關系，而無需事前揭示描述這種映射關系的數(shù)學方程，并且具有相當長的研究歷史，發(fā)展成熟。因此，本文利用BP神經(jīng)網(wǎng)絡模型訓練目標網(wǎng)絡，進行目標識別。

識別過程為：對某一像素點的n×n鄰域進行紋理特征提取，將特征向量輸入到目標神經(jīng)網(wǎng)絡中進行識別，輸出該像素點的特性判別值，若特性判別值落在接受域(ac1,ac2)之中，則斷定為疑似目標點(由于震象云的紋理特征不具備100%的識別力，因此用其特征進行區(qū)別時，會將一些非震象云區(qū)域的像素點也包含進來)，否則為非目標點。經(jīng)實驗表明，鄰域n×n中取n=51時較好，因為取樣時，發(fā)現(xiàn)在51× 51正方形區(qū)域內(nèi)可以較好地包含震象云，這樣可以盡量獲取震象云的紋理特征。

此過程的算法流程如下：

Step1特征提取。掃描以某一像素點為中心51× 51鄰域，統(tǒng)計目標灰度級區(qū)段(拉伸前是110～145，拉伸后是55～200)0°、45°、90°、135°等4 個方向上的灰度共生矩陣，然后分別計算每個矩陣的對比度、相關性、能量和同質性這4個特性的值。這樣將獲得16個紋理特征值組成的特征向量。主要計算公式如下：

上述公式中用到的頻度 p(i,j/d,θ) 、期望(ux,uy)、方差(σx,σy)等計算公式及各元素含義可參考文獻[9-10]。

Step2獲取目標神經(jīng)網(wǎng)絡：應用BP神經(jīng)網(wǎng)絡學習模型訓練目標網(wǎng)絡。輸入為震象云和非震象云的特征向量，特征向量按照Step1進行計算。期望輸出為特性判別值，震象云的為1.0，非震象云的為–1.0。按照BP神經(jīng)網(wǎng)絡模型理論[11-12]，設置訓練模型的相關參數(shù)，進行訓練，最終輸出目標神經(jīng)網(wǎng)絡。

Step3利用目標神經(jīng)網(wǎng)絡對輸入圖像進行目標識別和分類。掃描圖片，對于灰度級落在目標灰度區(qū)域內(nèi)的像素點統(tǒng)計其鄰域紋理特征，輸入到目標神經(jīng)網(wǎng)絡進行識別，并且將特性判別值落在接受域內(nèi)的標記為疑似目標點，否則標記為非目標點；對于落在目標灰度區(qū)域之外的像素點，直接標記為非目標點，不需要進行識別。

同一片震象云中的像素點雖然紋理是相似的，但并非一模一樣，因此，計算出的特征向量會有差異，通過神經(jīng)網(wǎng)絡映射后的特性判別值不一定會與正樣本的特性判別值相等，而是會在其值附近波動。上文實驗規(guī)定震象云樣本的特性判別值為1.0，實驗驗證目標神經(jīng)網(wǎng)絡輸出值接受域為(ac1,a c2)=(0.5,1.2)較好，這樣可以完整包含震象云區(qū)域的像素點。

2.2.3 目標過濾

對整張圖片的每個像素點進行處理之后，疑似目標點會聚積形成疑似目標區(qū)域。目標過濾就是過濾掉疑似目標區(qū)域中過大或過小的區(qū)域，其中過大的區(qū)域是背景的云層圖像，過小的區(qū)域是圖像中的噪聲或云團運動過程中拉伸出來的小空洞。具體實現(xiàn)為：統(tǒng)計每個8連通疑似目標區(qū)域的像素點個數(shù)Area(即面積)，若Area的取值落在接受域(a c3,a c4)之外，則標記為非目標塊，過濾掉；否則保留疑似目標。實驗證明，疑似目標塊面積Area接受域取值為(a c3,a c4)=(250,1 2 50)較為理想，既可以濾掉大的背景云團，又可濾掉小的噪聲或空洞。

2.2.4 目標跟蹤

目標跟蹤是對多張圖片的疑似目標塊進行處理的操作，分為目標標準化和跟蹤統(tǒng)計2個步驟：

(1)目標標準化：掃描每一個目標過濾剩下的疑似目標塊，計算其幾何中心(即點集在水平和豎直方向上的平均值)和半徑radious(公式為radious=)，即用圓形擬合跟蹤區(qū)域)。

(2)跟蹤統(tǒng)計：這個步驟需要建立疑似目標集，對每個出現(xiàn)出新的疑似目標進行判斷是否屬于已經(jīng)存在的目標集。判斷依據(jù)：如果疑似目標中心點和集合中心點的距離不超過前集合半徑和疑似目標半徑中的最大值，則該疑似目標屬于該集合。如果屬于某一個目標集，則將該疑似目標作為新元素加入到該集合中，調(diào)整這個疑似目標集的中心(即計算集合中所有元素中心的平均值)和半徑(取當前集合半徑和加入元素半徑中的最大值)；如果不屬于任何一個已存在的集合，則新建一個集合，將該疑似目標加入進去，將集合的中心和半徑初始化為疑似目標的中心和半徑。

當集合中的元素(疑似目標)個數(shù)超過一定的閾值5時，可以報告這個范圍出現(xiàn)的疑似目標為目標，根據(jù)目標出現(xiàn)的周期長短，預測地震的發(fā)生；否則，繼續(xù)跟蹤，若出現(xiàn)次數(shù)太少(如只有一兩次)且長期未出現(xiàn)，則可以刪除這個集合。這是因為目標(增溫型震象云)區(qū)域所對應的地面斷斷續(xù)續(xù)釋放能量，為其加熱，故其位置可以保留；而非目標區(qū)域為高層云運動時云團內(nèi)部速度不一致，導致撕裂偶然形成的，沒有能源支撐，故其不能保留。

3 實驗及結果分析

3.1 實驗

本文以汶川地震為例進行實驗。實驗對象為風云2C衛(wèi)星的長波紅外1通道云圖，云圖格式為.AWX，波長范圍為10.3～11.3。

以FY2C_ANI_IR1_R01_20071014_0800為例說明單張圖片處理過程，該圖像大小為1200×1200。首先根據(jù)云圖格式[13]讀入圖像數(shù)據(jù)，進行預處理，提取特征輸入到目標神經(jīng)網(wǎng)絡進行識別。

圖3為長波紅外云圖的原始數(shù)據(jù)所呈現(xiàn)的圖像。圖4為灰度經(jīng)過式(1)分段線性拉伸處理之后的圖像，從這幅圖像中可以明顯地看到在圖像中央的白色云團中有一個長條形的孤立灰黑色的云斑，這個云斑就是震象云，即圖1中所標識區(qū)域。圖5為經(jīng)過目標神經(jīng)網(wǎng)絡識別后留下的疑似目標所形成的圖像，從圖像中可以看出圖4中的震象云區(qū)域得到完整的保留，其中白色點集表示疑似目標。圖6是對圖5進行目標過濾之后剩下的疑似目標所形成的圖像。從這一系列圖中看出，經(jīng)過一系列處理之后，疑似目標區(qū)域只剩下少數(shù)幾個，為后續(xù)的跟蹤奠定了基礎。

圖3 長波紅外原始圖像

圖4 灰度分段線性拉伸

圖5 目標神經(jīng)網(wǎng)絡識別

圖6 目標過濾輸出

目標跟蹤：為了簡化跟蹤，先將目標標準化，再在標準化的目標區(qū)域進行跟蹤，統(tǒng)計同一區(qū)域疑似目標出現(xiàn)的次數(shù)，最后輸出結果。圖7為圖6目標標準化后的目標分布圖，圖8為從2007年10月到2008年4月底跟蹤一系列晚間云圖所得到的疑似云目標分布。

圖7 單張圖片的標準化目標

圖8 多張圖片的標準化目標

跟蹤晚間云圖的原因是晚間沒有陽光輻射且人類活動也相對白天較少，能夠較好地反映地殼釋放能量對大氣增溫的影響，有利于震象云跟蹤的準確性。圖9(a)是進行疑似目標跟蹤后，同一區(qū)域出現(xiàn)疑似目標的次數(shù)和分布的示意圖，同一矩形區(qū)域內(nèi)的點表示該區(qū)域出現(xiàn)過的疑似目標。當目標出現(xiàn)次數(shù)達到一定的值時，這些疑似目標就很有可能是同一個地震帶活動釋放能量所產(chǎn)生的震象云。圖9(b)為出現(xiàn)2次或2次以上疑似目標區(qū)域的分布示意圖，圖9(c)為出現(xiàn)3次或3次以上疑似目標區(qū)域的分布示意圖，圖9(d)為出現(xiàn)8次或8次以上疑似目標區(qū)域的分布示意圖。

圖9 跟蹤結果

在實驗中進行目標跟蹤后，得到同一區(qū)域疑似目標出現(xiàn)次數(shù)與區(qū)域個數(shù)的關系如下：只出現(xiàn)疑似目標1次的區(qū)域有84個，出現(xiàn)疑似目標2次的區(qū)域有9個，出現(xiàn)疑似目標8次的區(qū)域只有1個，沒有區(qū)域出現(xiàn)3次～7次疑似目標。本實驗最終結果是將出現(xiàn)8次的疑似目標斷定為目標，即將在同一區(qū)域出現(xiàn)的8次疑似目標確定為震象云。這些震象云的幾何中心均包含在圖片中以(688,690)為中心，半徑是19單位像素的圓形區(qū)域內(nèi)。圖10為本實驗最終跟蹤到的目標，按照時間序列由遠及近8次出現(xiàn)時的區(qū)域形態(tài)圖，即汶川地震震前的震象云所在區(qū)域圖。

圖10 汶川地震震前震象云出現(xiàn)區(qū)域

3.2 結果分析

根據(jù)各圖震象云出現(xiàn)的位置，結合當時風級、風向對震象云帶來的偏移效果，將圖像坐標轉化成地理坐標，預測震中坐標是北緯30.885°± 0.1°，東經(jīng)103.632°± 0.1°，而汶川地震震中是北緯30.986°，東經(jīng)103.364°，雖沒有精準地與實際震中重合，但偏移范圍較小。

從實驗結果看，汶川地震發(fā)生前，其上空至少出現(xiàn)了8次的增溫型震象云，第一次檢測到該區(qū)域的震象云是2007年10月14號(云圖名為FY2C_ANI_IR1_R01_20071014_0800)，最后一次是2008年4月13日(云圖名為FY2C_ANI_IR2_R01_20080413_1000)，距離地震發(fā)生時間5月12日是29天，屬于短期地震預測。

4 結束語

本文利用圖像識別技術提取震象云，實現(xiàn)了地震預測的自動化，因為基于震象云的地震預測本身具有震中位置預測較準確的特點，這就彌補了傳統(tǒng)熱紅外地震預測方法中的主要缺陷。

由于震象云形成過程非常復雜：圖像上的震象云云團體現(xiàn)的是對應位置對流層云頂?shù)臏囟确植迹贿@個溫度分布不僅與地殼活動強烈和釋放的各物質成分有關，還和當時周圍的大氣壞境有關。下一步的研究方向是單純從震象云的形態(tài)上準確預測地震震級。

[1]傅承義.地球十講[M].北京:科學出版社,1976.

[2]Gorny V I,Salman A G,Tronin AA,et al.The Earth Outgoing IR Radiation as an Indicator of Seismic Activity[J].Proceedings of the USSR Academy of Sciences,1988,30(1):67-69.

[3]強祖基,孔令昌,王戈平,等.地球放氣，熱紅外異常與地震活動[J].科學通報,1992,37(24):2259-2262.

[4]強祖基,孔令昌,郭滿紅,等.衛(wèi)星熱紅外增溫機制的實驗研究[J].地震學報,1997,19(2):197-201.

[5]強祖基,賃常恭,李玲芝,等.衛(wèi)星熱紅外圖像亮溫異常-短臨震兆[J].中國科學:D輯,1998,28(6):564-573.

[6]徐秀登,張行才,李貴達.張北地震與大氣增溫異常[J].西北地震學報,2000,22(1):44-47.

[7]徐秀登,徐向民.地震前紅外異常的基本特征與成因機制[J].西北地震學報,2001,23(3):310-312.

[8]劉德富,康春麗.地球長波輻射(OLR)遙感與重大自然災害預測[J].地學前緣,2003,10(2):427-435.

[9]侯海苗,冀小平.基于灰度共生矩陣的紋理特征[J].長治學院學報,2008,25(5):31-32.

[10]寧順剛,白萬民,喻 鈞.基于灰度共生矩陣的圖像分割方法研究[J].電子科技,2009,22(11):69-71.

[11]樊振宇.BP神經(jīng)網(wǎng)絡模型與學習算法[J].軟件導刊,2011,10(7):66-68.

[12]李秀珍,孔紀名,李朝鳳.基于Matlab的BP神經(jīng)網(wǎng)絡在泥石流危險性評價中的應用[J].工程勘察,2010,38(1):47-50.

[13]許健明,張文建,楊 軍,等.風云二號衛(wèi)星業(yè)務產(chǎn)品與衛(wèi)星數(shù)據(jù)格式實用手冊[M].北京:氣象出版社,2008:229-232.