基于光學和SAR遙感圖像融合的洪災區(qū)域檢測方法

王志豪 李 剛 蔣 驍

(清華大學電子工程系 北京 100084)

1 引言

遙感技術為地理勘探[1]、建筑測繪[2]、植被監(jiān)測[3]、水域監(jiān)測[4]、災區(qū)探測[5]等提供了重要的技術支持。在洪災區(qū)域檢測中，光學遙感能夠提供豐富的光譜信息和較高的圖像分辨率，而雷達遙感則具備全天候的洪災區(qū)域檢測能力[6]。異質圖像融合已經在復Contourlet域得到了實現[7]，因此結合光學遙感數據和雷達遙感數據能實現洪災區(qū)域的全天候、高時效檢測[8]。

當前，學者們提出許多基于遙感圖像的洪災區(qū)域檢測方法。其中，應用人工神經網絡和自組織映射網絡[9]方法處理洪災前后的同質SAR圖像時能大致區(qū)分洪災區(qū)域，但是受到大量相干斑噪聲的影響，檢測結果仍然存在很高的虛警率。為了結合光學遙感圖像和SAR圖像的成像優(yōu)勢，一些學者提出了基于異質圖像的洪災區(qū)域檢測方法，例如像素變換法[10]和特征空間變換法[11]。上述方法對洪災前的SAR圖像和洪災后的光學圖像有良好的檢測效果，但是對洪災前的光學圖像和洪災后的SAR圖像的檢測效果較差。

針對SAR圖像中洪災區(qū)域的檢測難點，本文提出了H-FCM算法。該算法的主要貢獻如下：

(1) 提出SAR圖像的分級聚類模型；

(2) 將洪災前的近紅外波段光學遙感圖像中提取出的河流位置融合到洪災后的SAR圖像中，并作為洪災區(qū)域檢測的空間約束，進一步提升了檢測結果的準確率。

文章結構如下：第2節(jié)分析了近年來洪災區(qū)域檢測方法的發(fā)展和存在的問題。第3節(jié)介紹了分級C均值聚類算法Hierarchical Fuzzy C-Means (H-FCM)算法的原理、步驟以及參數設定對算法性能的影響。第4節(jié)應用包括H-FCM算法的4種算法對不同類型遙感數據進行洪災區(qū)域檢測，并分析了圖像降噪處理對4種算法檢測性能的影響。第5節(jié)系統(tǒng)總結了全文的主要工作和存在的不足。

2 洪災區(qū)域檢測方法的發(fā)展

2.1 傳統(tǒng)的洪災區(qū)域檢測方法

目前，學者們提出的基于雷達遙感圖像的水體檢測方法包括灰度閾值分割法、濾波法、機器學習法和結合輔助信息的提取方法等[12]；基于光學遙感圖像的提取水體方法包括基于像元分類的閾值法和基于目標分類的分類法等[12]。洪災是由于大雨、暴雨或者持續(xù)降雨造成的，此時無法通過實時獲取的光學遙感數據檢測地面的洪災區(qū)域。

以2019年8月9日對利奇馬臺風實時追蹤的GF4數據(見圖1)為例，光學圖像無法顯示地表的洪災區(qū)域。但是通過適當的衛(wèi)星調度方案可以短時間內獲取雷達遙感數據。以2019年8月8日對利奇馬臺風實時追蹤的GF3數據(見圖2)為例，SAR圖像可以清晰地顯示地面受災情況。

應用閾值分割法中的最大類間方差算法[13]對災后的SAR圖像進行洪災區(qū)域檢測時，由于受到相干斑噪聲的影響，檢測精度較低。為克服相干斑噪聲的影響，學者們提出了濾波法[14]、邊界追蹤法[15]、Markov分割法[16]、基于鄰域最小生成樹的半監(jiān)督分類法[17]、基于概率轉移卷積神經網絡的分類法[18]等。但是洪災區(qū)域檢測問題的核心在于保證高檢測準確率的同時，盡可能縮短檢測時間，而上述方法的計算復雜度普遍較高，并不能滿足洪災區(qū)域檢測的高時效性需求。

圖1 利奇馬臺風的GF4圖像(資源衛(wèi)星中心)Fig.1 The GF4 image of Typhoon Lekima(Resource Satellite Centre)

圖2 利奇馬臺風的GF3圖像(資源衛(wèi)星中心)Fig.2 The GF3 image of Typhoon Lekima(Resource Satellite Centre)

直接利用K-means算法[19]、C均值聚類方法(FCM)算法[20]等無監(jiān)督聚類算法對洪災后的SAR圖像進行洪災區(qū)域檢測可以縮短檢測時間，但是檢測結果的分類準確率很低。具體原因如下：根據灰度閾值分割法的原理，在SAR圖像中，水體的散射值較低，整體呈現為暗區(qū)。在沒有發(fā)生暴雨或者持續(xù)降雨的情況下，可以通過分析SAR圖像的整體像素值分布，應用FCM算法等聚類算法確定水體的分割閾值。根據分割閾值，將大于閾值的標記為非水體，小于閾值的標記為水體。但是洪水災害發(fā)生時伴隨著暴雨或者持續(xù)降雨，導致很多未受到洪災影響的區(qū)域在SAR圖像中也呈現為暗區(qū)。因此，直接通過閾值法獲得的檢測結果存在較高的虛警率。

2.2 FCM算法原理

FCM算法的步驟總結為表1。以洪災后的大小為n1×n2像素的SAR圖像為例，將其所有像素點的像素值記作二維矩陣X。將X分成c類，則可以得到c個聚類中心。將每一個像素點作為單個樣本xj，則xj對于第i類聚類中心的隸屬度可記作uij,FCM算法的目標函數、約束條件的定義式分別為式(1)和式(2)

其中，Ci表示第i類聚類中心；m表示隸屬度因子；n=n1×n2;i=1,2,…,c;j=1,2,…,n。

結合式(2)的約束條件，利用拉格朗日乘數法，將J分別對uij和Ci求導，并令其結果等于0，可得uij和Ci的迭代式分別為式(3)和式(4)

3 H-FCM算法原理

本文提出的H-FCM算法致力于縮短SAR圖像洪災區(qū)域的檢測時間、降低洪災區(qū)域的錯誤檢測概率以及抑制SAR圖像中存在的大量隨機分布的相干斑噪聲對檢測結果的影響。根據河流連續(xù)體(River Continuum Concept,RCC)[21]的概念，在已知的河流系統(tǒng)內，自上游到下游的物理變量的梯度變化是連續(xù)的[22]。因此發(fā)生洪水災害之后，若地理環(huán)境沒有發(fā)生劇烈變化，洪水的形態(tài)仍然符合RCC。并且，根據典型地物的光譜特征曲線圖[23](見圖3)，水體對近紅外波段的電磁波吸收能力比大部分地物強。

表1 FCM算法Tab.1 FCM algorithm

H-FCM算法的應用前提為可獲得洪災前的光學圖像以及洪災后的SAR圖像。該算法通過對遙感圖像的像素強度劃分進行分級聚類。以uint8格式的圖像數據為例，應用FCM算法獲得待測數據像素值的8個聚類中心，并且由低至高排列。根據圖3,SAR圖像中代表水體的像素點的像素值主要分布在0～64之間，即處于低像素值區(qū)間。但是由于發(fā)生洪災時，SAR圖像中存在大量隨機分布的相干斑噪聲以及低洼地區(qū)存在的積水，直接用閾值法進行洪災區(qū)域檢測的準確率較低。

H-FCM算法采用反向聚類的方法，考慮到在SAR圖像中，代表水體的像素點幾乎不可能處于高像素值區(qū)間。H-FCM算法在聚類過程中將處于最高級的像素值區(qū)間的像素點記為1，其余像素點記為0，對所有零像素點進行區(qū)域聚類，獲得初步的洪災區(qū)域檢測結果。又由RCC原理，洪災發(fā)生后，主要受災區(qū)域與洪災前的河流相連通，利用災害前的河流位置作為空間約束可以得到近似真實洪災區(qū)域的檢測結果。因此，應用H-FCM算法檢測洪災發(fā)生后的SAR圖像的洪災區(qū)域，可以降低相干斑噪聲對檢測結果的影響。

H-FCM算法主要分為4步。首先，應用FCM算法獲得洪災后SAR圖像的聚類中心，進而設置分類閾值，并確定分級聚類模型；其次，提取災害前近紅外波段光學遙感數據中的河流位置，并將其融合到洪災后的SAR圖像中；再次，以聚類模型為基礎進行分級聚類，獲得洪災區(qū)域的初步聚類結果；最后，通過空間約束矩陣對初步聚類結果進行約束，進一步提升洪災區(qū)域檢測的準確率。

圖3 典型地物的光譜特征曲線[23]Fig.3 Spectral characteristic curves of typical features[23]

3.1 建立分級聚類模型

利用FCM算法獲得洪災后SAR圖像的八分類的聚類中心矩陣C和閾值矩陣T，并由此獲得7個聚類模型

其中，一維矩陣C表示聚類中心矩陣，包含8個按照升序排列的聚類中心

其中，T0=0,q=1,2,…,7。

以1999年英國格洛斯特洪災后獲得的SAR圖像Ea[10]為例，其大小為n1×n2(見圖4)。遍歷Ea中的所有像素點，若其像素值大于Tq–1且小于Tq則在對應位置標記為1，否則標記為0，將所得到的矩陣記作聚類模型Gmq，如式(8)

其中，“(Gmq)m n∈Gmq”在本文中均代表(Gmq)mn是Gmq矩陣中第m行，第n列的元素，下文符號“∈”的含義與此處一致；(Ea)m n∈Ea;q=1,2,…,7。

3.2 災后SAR圖像融合災前河流空間信息

3.2.1 SAR圖像洪災區(qū)域檢測的改進方法

圖4 洪災后的SAR圖像Ea[10]Fig.4 SAR imageEa after the flood[10]

由于SAR圖像的成像模式一般為寬幅掃描模式，同時受到SAR圖像中相干斑噪聲的影響，在沒有災害前河流信息的情況下直接采用無監(jiān)督聚類的方法來檢測Ea中的洪災區(qū)域存在較高的虛警率，導致檢測結果的準確率較低。H-FCM算法將洪災前河流的空間信息融合到洪災后的SAR圖像中，通過局部空間約束，可以提高洪災區(qū)域檢測的準確率。但是在圖4中，Ea的河流位置有一部分被背景地物或者噪聲覆蓋，另一部分和洪災區(qū)域重合，導致大部分的河流位置已經無法識別。為了獲得洪災前的河流位置信息，H-FCM算法進一步檢測洪災前的光學圖像。

以1999年英國格洛斯特洪災前的SPOT圖像中提取的近紅外波段圖像Ei[10]為例(見圖5)。由于不同河段的河流泥沙含量不同，造成不同河段對近紅外波段電磁波的吸收能力不同[24]，通過閾值法無法直接獲得所有河流區(qū)域的空間位置?？紤]到洪災前的水體表面近似為平滑面，可以應用區(qū)域生長算法[25]提取洪災前近紅外波段遙感圖像的河流位置。本文提出的區(qū)域生長算法為表2，根據圖5中近紅外波段的光學遙感圖像可以利用該算法獲得洪災前的河流位置。

圖5 洪災前的近紅外波段影像Ei[10]Fig.5 Near-infrared imageEi before the flood[10]

表2 區(qū)分河流和道路的區(qū)域生長算法Tab.2 Region growing algorithm to distinguish rivers and roads

3.2.2 檢測河流位置的區(qū)域生長算法規(guī)則

由于步驟(2)獲得的第2類閾值圖Pt(見圖10)中存在大量的區(qū)塊噪聲的干擾，無法從中直接獲得災害前的具體河流的空間位置。因此，需要篩除大面積區(qū)塊性的非河流區(qū)域。像素點的面積定義為：將像素值為1的像素點的面積記作1，像素值為0的像素點面積記作0。連通區(qū)域的面積定義為：將連通區(qū)域內所有非零像素點的面積累加得到該區(qū)域的面積。4連通區(qū)域的定義為目標像素點(紅色方塊)的上、下、左、右4個像素點位置(見圖11)。

為了充分消除區(qū)塊噪聲的影響，并保留河流區(qū)域的完整性，步驟(3)的卷積操作分為兩步。首先，確定大小為α×α的全1卷積核Wα，將Pt與Wα卷積得到初步消除區(qū)塊噪聲的結果P1，如式(9)和式(10)；其次，將P1與Wβ卷積得到基本消除區(qū)塊噪聲的結果P2，如式(11)和式(12)。為保證卷積之后的輸出結果的尺寸與原圖像大小保持一致，本文所有的卷積模式均采取matlab中conv2函數的‘same’模式進行卷積

圖6 方向矩陣Jd1(向下)Fig.6 Direction matrix Jd1(down)

圖7 方向矩陣Jd2(向上)Fig.7 Direction matrix Jd2(up)

圖8 方向矩陣Jd3(向左)Fig.8 Direction matrix Jd3(left)

圖9 方向矩陣Jd4(向右)Fig.9 Direction matrix Jd4(right)

3.2.3 區(qū)域生長算法的參數選擇

參數α一般為固定的大小，當選擇的α較大或者較小時，在篩除區(qū)塊噪聲的同時，也篩除了河流。本文固定參數α=18的卷積核Wα初步篩除大面積聚點，獲得P1；其次，固定參數β=9的卷積核Wβ進一步篩除小面積聚點，獲得P2。利用式(13)至式(15)確定河流區(qū)域種子生長點的位置Po，在Po的每一個4連通區(qū)域中等間距選擇10個像素點作為初始生長點(見圖12)。

圖10 第2類閾值圖PtFig.10 Threshold graphPt

圖11 4連通區(qū)域的定義Fig.11 Definition of four-connected regions

其中，Objj表示P2中的第j個4連通區(qū)域；P2中4連通區(qū)域總數為na；每個連通區(qū)域的背景大小均為n1×n2;j=1,2,…,na

其中，Areaj表示P2中第j個4連通區(qū)域Objj的面積。

3.2.4 利用區(qū)域生長算法獲得河流位置

根據圖12中獲得初始生長點，按照步驟(5)進行4個方向的區(qū)域生長，最后將各個方向的生長結果累加得到初步生長結果(見圖13)。但是由于本例中道路和河流交錯，當道路和河流寬度相近且沒有先驗測繪信息時，直接通過機器識別提取并區(qū)分河流和道路的空間位置是非常困難的。所以在步驟(7)中利用霍夫變換[26]檢測初步生長結果中屬于公路的長距離直線段區(qū)域，進而識別道路所在區(qū)域(見圖14)。最后，將初步生長結果去除道路區(qū)域后，得到河流位置Eb(見圖15)。

3.2.5 災后SAR圖像與災前河流位置融合

圖12 初始生長點位置Fig.12 Initial growth point location

圖13 初步生長結果Fig.13 Preliminary growth results

將洪災前的河流位置Eb和洪災后的SAR圖像Ea融合。融合規(guī)則為將Ea中與Eb中災前河流位置對應的所有像素點的像素值記作T1，其余像素點的像素值保持不變，如式(16)，將融合后的災后SAR圖像記作Fu(見圖16)

其中，(Fu)ij∈Fu,(Eb)ij∈Eb,(Ea)ij∈Ea。

3.3 分級聚類

圖14 道路的識別結果Fig.14 Road recognition results

圖15 洪災前的河流位置EbFig.15 River locationEb before the flood

圖16 融合后的災后SAR圖像FuFig.16 Fusion post-disaster SAR image Fu

光學圖像的水體和非水體的區(qū)分可以通過計算歸一化水指數(Normalised Difference Water Index,NDWI)[27]檢測，但是SAR圖像中不含多通道圖像信息，不能直接計算其NDWI來區(qū)分水體和非水體。通過觀察融合后SAR圖像的高像素值區(qū)間Gm7(見圖17)，可以發(fā)現真實洪災區(qū)域的像素點在其中幾乎全部表現為零像素值。同時，根據RCC原理，真實洪災區(qū)域在整個圖像空間呈現較高的區(qū)塊性的稀疏度。即代表洪災區(qū)域的非零像素值區(qū)域是一個連通區(qū)域，其余零像素值區(qū)域均為非洪災區(qū)域，所以作為洪災區(qū)域檢測的聚類模型應當具備較高的稀疏度。Gm7是一個反向模型，洪災區(qū)域基本呈現為零像素區(qū)域，非洪災區(qū)域基本呈現為非零像素的散點分布。

在目標檢測中，文獻[28]定義圖像稀疏度為非零像素點個數K1，考慮到洪災區(qū)域檢測的目標地物覆蓋范圍較廣，本文進一步引入文獻[29]中圖像稀疏占比的定義。即模型的稀疏占比為待測圖像中非零像素點個數和所有像素點個數的比例，稀疏占比K2的定義式為公式(17)。聚類模型的稀疏度分析，以Gm7聚類模型為例，檢測結果如圖18。

圖17 聚類模型Gm7Fig.17 Clustering model Gm7

圖18 聚類模型的稀疏度分析Fig.18 Analysis of sparsity of clustering model

通過分析聚類模型Gm7的稀疏度，得到其稀疏占比為0.10905，可以確認該模型本身具有較高的稀疏度。并且該模型的非零像素點代表最高級像素值的分布，實際洪災區(qū)域在該模型中基本處于零像素值區(qū)域。對聚類模型Gm7進行卷積運算后，根據卷積結果得到每一個像素點周圍非零像素點的總數，并由此判斷該像素點是否為洪災區(qū)域。

其中，Wk表示大小為k×k的全1卷積核。

H-FCM算法對聚類系數φ和卷積系數k的選擇較為敏感。當選擇的φ過大時，聚類結果包含較多的非洪災區(qū)域；當選擇的φ較小時，聚類結果舍棄了較多真實洪災區(qū)域。相應地，當選擇的k過大或者過小時，聚類結果都將包含較多的非洪災區(qū)域。為使這兩類參數選擇達到最優(yōu)檢測性能，應當在選擇最佳聚類模型的同時，依據模型的稀疏度來確定相應的聚類系數和卷積系數。考慮到SAR圖中大量隨機分布的相干斑噪聲，為減弱噪聲對檢測性能的影響，φ,k的選擇要比實際的稀疏占比更小，φ,k與K2的關系定義為式(19)。

圖19 初步聚類結果Gr7Fig.19 Preliminary clustering results Gr7

3.4 洪災區(qū)域檢測的空間約束

H-FCM算法在對洪災區(qū)域的空間約束上提出了一種參數自適應的空間約束方法，有效地篩除疑似洪災區(qū)域，并提升檢測性能。由于洪災區(qū)域的水體滿足RCC原理，定義洪災區(qū)域主要范圍為Rt。Rt表示初步聚類結果Gr7的4連通區(qū)域面積最大的連通區(qū)域。

由于真實洪災區(qū)域在整個圖像空間存在較高的稀疏度，即非洪災區(qū)域處于零像素值區(qū)域，洪災區(qū)域呈現為非零像素值的連通區(qū)域。此時，等間距增大空間約束范圍會增加計算復雜度。引入約束矩陣CSβp作為初步聚類結果的空間約束，當約束系數βp不斷減少時，整體的空間約束范圍按照e指數級不斷增大。由于真實洪災區(qū)域在全局上具備較高的稀疏度，約束系數βp按照等間隔減小的過程，對應約束范圍呈e指數級的擴大。因此，式(21)的搜索規(guī)則符合洪災區(qū)域的區(qū)塊稀疏分布特點，既能保證獲得適應待檢測圖像的約束系數，又不增加額外的運算量。

約束系數βp的自適應選擇方法，主要是依據兩個因素：

(1) 根據RCC原理，洪災發(fā)生時，洪災區(qū)域在遙感圖像上滿足4連通區(qū)域的連通性；

(2) 遠離洪災前河流位置的低洼地區(qū)在洪災發(fā)生時與洪災區(qū)域相鄰接的概率較低。

其中，e為自然常數，p=1,2,…,100

本文采用曼哈頓距離來表示像素點之間的空間距離，可以在保證檢測性能的基礎上，降低算法的運算量。洪災區(qū)域的空間約束式exp(-dij/dmax)的取值范圍為(0.367,1](見圖20)。并且當dij越小時，exp(-dij/dmax)的值越接近1，當d ij越大時，exp(-dij/dmax)的值越接近0.367，確保距離災害前河流位置越近的像素點被檢測為洪災區(qū)域的概率越大。

圖20 洪災區(qū)域空間約束曲線Fig.20 Spatial constraint curve of flooded area

將洪災區(qū)域主要范圍Rt加上空間約束CSβp，可以進一步提升洪災區(qū)域檢測結果的準確率。由式(23)可以獲得約束系數為βp時的洪災區(qū)域檢測結果Outβp。

其中，“A·B”表示矩陣A與矩陣B的點乘運算，即矩陣A和矩陣B對應點相乘。

約束系數βp自適應選擇的目的是判定疑似災區(qū)S是否為真正災區(qū)，以獲得洪災區(qū)域的最終檢測結果。所以在βp的搜索過程中，若對應的疑似災區(qū)的面積增長比率遠小于聚類模型本身的稀疏度時，則判定疑似災區(qū)為未受災區(qū)域

其中，“–”表示矩陣對應點像素值作差

其中，“≈”近似的誤差為±0.001;p=1,2,…,100

其中，Sj表示S中的第j個4連通區(qū)域，其背景大小為n1×n2

其中，Aj為Sj的面積，α1,α2為判定系數

其中，lj表示重合邊界長度；(CHj)mq∈CHj; CHj為Sj與Outβa的重合邊界，其寬度為2個像素單元。

將得到的約束系數βa,βb帶入到式(24)，獲得待定的疑似災區(qū)S。根據RCC原理，若疑似的4連通區(qū)域為真實洪災區(qū)域，則其像素面積Aj與Sj和Outβa鄰接的邊長lj的比值應處于(α1,α2)區(qū)間內。當面積邊長比過大時，說明疑似災區(qū)與真實洪災區(qū)域相鄰接的邊界較短，可判定為非災區(qū)；當面積邊長比過小時，說明疑似災區(qū)面積較小，屬于區(qū)塊噪聲的可能性比較大。最后，根據每一個疑似區(qū)域的判定結果dtj確定該區(qū)域是否為洪災區(qū)域，并將判定結果與Outβa累加得到最終的洪災區(qū)域檢測結果FF。

4 實驗結果與分析

4.1 檢測結果的評價指標

(1) 正確檢測概率。Righta表示所有像素點被正確檢測為變化的概率與被正確檢測為未變化的概率之和

其中，Nr_c表示正確檢測為變化的像素點總數，Nr_uc表示正確檢測為未變化的像素點總數，Nt表示像素點總數。

(2) 遺漏檢測概率。Missa表示所有實際發(fā)生變化的像素點沒有被檢測為發(fā)生變化的概率

其中，Nc表示實際發(fā)生變化的像素點總數；Nm表示實際變化的像素點未被檢測為變化的像素點總數。

(3) 錯誤檢測概率。Wra表示所有像素點被錯誤檢測的概率，即實際未發(fā)生變化的像素點被錯誤檢測為發(fā)生變化的概率與實際發(fā)生變化的像素點被檢測為未變化的概率之和

其中，Nf_u表示實際發(fā)生變化的像素點被檢測為未變化的總數；Nun_c表示實際未發(fā)生變化的像素點總數，Nf_c表示實際未發(fā)生變化像素點被檢測為變化的總數。

(4) 卡帕系數

其中，po表示總體分類精度。

Nd_c表示檢測為發(fā)生變化的像素點總數，Nd_uc表示檢測為未變化的像素點總數。

(5) 計算時間。以計算時間衡量算法的時間復雜度，計算時間定義為Time，其單位為“s”(秒)。

(6) 迭代次數。以運算的迭代次數體現算法的運算復雜度。迭代次數定義為Iter，其單位為“次”，表示對應算法獲得聚類中心的迭代次數。

4.2 實驗數據的預處理和參數選擇

本文的實驗數據是在確定目標災害區(qū)域地理坐標的基礎上，獲得洪災前遙感圖像和洪災后遙感圖像的歷史數據。實驗數據的預處理過程主要分為4步。首先，利用ENVI5.3軟件的Orthorectification模塊對洪災前后的遙感數據分別進行正射矯正；其次，利用ENVI5.3軟件Registration模塊中的Image to Image方式，以洪災前的光學圖像為基準，選擇關鍵角點為控制點對同地區(qū)的其它遙感數據進行配準；再次，裁剪相同大小配準后遙感數據作為實驗數據；最后，利用ENVI軟件獲得洪災后光學圖像的NDVI圖像，通過人工標注的方法確定真實洪災區(qū)域，作為各類算法檢測結果性能指標評定的參照對象。

H-FCM算法的參數選擇，只需要確定判定系數α1,α2。卷積系數k、聚類系數φ可以通過分析聚類模型的稀疏度獲得，約束系數βa,βb可以通過自適應選擇獲得。判定系數的范圍是通過統(tǒng)計英國格洛斯特的歷史洪災數據的貝葉斯更新結果[30]獲得的最優(yōu)區(qū)間，在H-FCM算法中被設置為固定的區(qū)間。并且該區(qū)間定義在4.3節(jié)中實驗1和實驗2的結果中均獲得了較好的檢測性能。

為了驗證H-FCM算法在洪災檢測領域的普適性，本文通過收集案例中的河流歷史洪災信息，利用歷史實測數據來確定不同河段發(fā)生洪災的空間約束系數?；谠摵恿鞯臍v史災害數據進行貝葉斯更新[30]，可以發(fā)現某一河段附近地區(qū)受洪災影響的范圍與該河段的河流寬度、河床高度、河道轉角系數密切相關。參考英國格洛斯特的歷史災害數據[31]，可以得到洪災區(qū)域的約束系數βa,βb與河流平均寬度lm的近似對應關系

其中，lm通過計算圖15的平均河寬得到。

計算實際數據得到dmax=1338,lm=4.5，其單位均為像素單元。依據真實災害數據統(tǒng)計的近似對應關系得到的βa=0.8997,βb=0.7719。對比自適應選擇的約束系數，βa的選擇結果與洪災的統(tǒng)計歷史數據較為吻合，但是βb的選擇結果與洪災的統(tǒng)計歷史數據有一部分偏差。這是由于該流域的歷史洪災規(guī)模不同，本例所體現的洪災區(qū)域要大于歷史平均水平，所以本例的βb的自適應選擇結果要小于根據歷史統(tǒng)計得出的結果，對應實際更大的洪災區(qū)域。

4.3 4種檢測方法的比較

為充分論證H-FCM算法對洪災前光學圖像和洪災后SAR圖像的變化檢測性能，在對比實驗中加入了圖像分割算法—分水嶺算法(Watershed Algorithm,WA)[32]，邊緣檢測算法(Snake)[33],H-Kmeans算法等3種方法作為對比實驗方法。此外，將1999年英國格洛斯特洪災數據、2019年中國南昌附近洪災的洪災數據等兩個場景作為實驗數據。

4.3.1 實驗1

以1999年英國格洛斯特洪災前的近紅外波段圖像和洪災后的SAR圖像[10]為例，檢驗H-FCM算法的有效性。將融合后的SAR圖像(見圖16)記作Fu，其大小為2359×1318像素，采用的聚類模型為Gm7。設定判定系數α1=120,α2=160，計算獲得其卷積系數k=85，聚類系數φ=0.054。通過自適應選擇得到的約束系數βa=0.9199,βb=0.7199，最終的檢測結果為圖21。

首先，將實驗1中洪災后的NDVI圖[10](見圖22)進行人工標注獲得洪災區(qū)域的真值圖(見圖23)。WA算法的洪災區(qū)域檢測方法是根據整體像素值分布(見圖24)，將所有像素點進行二值化，然后求出零值與最近非零值的距離，畫出分水嶺，進而對圖像進行二分類。Snake算法獲得Ea的邊緣輪廓(見圖25)，按照式(18)和式(20)對邊緣輪廓檢測結果進行卷積運算得到最終洪災區(qū)域檢測結果。H-Kmeans算法是在Kmeans聚類(八分類)的基礎上，進行分級聚類，從而得到洪災區(qū)域的檢測結果。

圖21 最終檢測結果FFFig.21 Final experimental result FF

圖22 洪災后的NDVI圖像[10]Fig.22 NDVI images after the flood[10]

圖23 人工標注的真值圖Fig.23 Manually labeled truth map

圖24 像素值分布三維圖Fig.24 3-D map of the distribution of pixel values

圖25Ea的邊緣輪廓Fig.25 Edge profile ofEa

其次，將H-FCM算法、WA算法、Snake算法以及H-Kmeans算法獲得的洪災區(qū)域檢測結果與真值圖(圖23)比較，得到4種算法的檢測結果與真實災區(qū)的比較結果(見圖26)。在圖26中，白色部分表示實際洪災區(qū)域被檢測為洪災區(qū)域；綠色部分表示實際洪災區(qū)域被檢測為非洪災區(qū)域；紅色區(qū)域表示實際非洪災區(qū)域被檢測為洪災區(qū)域；黑色部分表示實際非洪災區(qū)域被檢測為非洪災區(qū)域。本文后續(xù)實驗的檢測比較結果與圖26的表示方法一致。

根據4.1節(jié)中對6項評價指標的規(guī)定，按照4種算法分別用matlab進行計算分析，得到每一種算法檢測結果對應的性能指標(見表3)。觀察表3的各項評價指標可以發(fā)現H-FCM算法有最高的正確檢測概率、最低的錯誤檢測概率以及最高的Kappa系數。H-Kmean算法雖然有著與H-FCM算法接近一致的檢測性能，但是其計算成本較高。綜合考慮洪災區(qū)域檢測的高準確檢測率和高時效性的要求，H-FCM算法具有相對最優(yōu)的檢測性能。

4.3.2 實驗2

實驗2的數據來源于歐洲航天局的開源數據，分別為哨兵1號和哨兵2號衛(wèi)星在2019年6月對中國南昌附近的觀測數據(見圖28)。按照H-FCM算法，首先，將圖28(a)中提取的的洪災前河流空間位置與圖28(d)中洪災后的SAR圖像融合；其次，對融合圖像進行分級聚類；最后，對聚類結果加以空間約束并獲得檢測結果。根據圖28(c)的洪災后的光學數據，人工標注真實洪災區(qū)域范圍(見圖27)。比較4種檢測方法對融合后SAR圖像的洪災區(qū)域檢測結果(見圖29)和性能指標(見表4)，可以得出H-FCM算法具有相對最優(yōu)的檢測性能。

實驗2的參數設定如下，將融合后的SAR圖像記作NCa，其大小為1025×1025像素，采用的聚類模型為Gm7；設定判定系數α1=120,α2=160；卷積系數k，聚類系數φ和約束系數βa,βb均基于NCa由式(19)和式(26)獲得。

4.4 圖像降噪處理對算法檢測性能的影響

Refined Lee filter(RL)[34]在濾除相干斑噪聲的同時能較好地保持SAR圖像的邊緣特征。為分析圖像降噪處理對上述算法檢測性能的影響，首先對實驗1中洪災后的SAR圖像進行RL濾波降噪，再利用上述4種方法對降噪處理后的圖片進行洪災區(qū)域檢測。定義新的檢測方法分別為RL_H-FCM算法、RL_WA算法、RL_Snake算法、RL_H-Kmeans算法，其檢測結果如圖30所示，檢測性能評價指標如表5所示。

表3 融合后SAR圖像洪災區(qū)域的檢測結果Tab.3 Detection results of flooded area in fusioned SAR image

圖27 人工標注的洪災區(qū)域Fig.27 Manually labeled flooded areas

比較表3和表5中上述4種算法對應的性能指標，可以發(fā)現RL_Snake算法的檢測性能提升最為明顯。這是由于Snake算法是先獲得圖像的邊緣特征，再對邊緣輪廓的內外進行區(qū)分。通過RL濾波可以減少相干斑噪聲的影響，使得圖像整體像素值的均值和方差更穩(wěn)定，有利于提取圖像的邊緣特征。但是通過RL濾波對H-FCM算法、WA算法和H-Kmeans算法的性能提升有限，反而極大地增加了算法的計算復雜度?？傮w而言，利用H-FCM算法直接對洪災前的光學圖像和洪災后的SAR圖像進行洪災區(qū)域檢測時，其檢測性能與降噪處理后的圖像檢測性能幾乎一致，而直接檢測的運算時間大幅度減少。所以應用H-FCM算法直接處理配準后的光學和SAR圖像，既可以縮短洪災區(qū)域的檢測時間又可以保證相對較高的檢測準確率。

5 結束語

圖28 洪災前后哨兵1，哨兵2數據Fig.28 Sentinel 1,2 data before and after flood

圖29 4種算法的洪災區(qū)域檢測結果Fig.29 Detection results of flooded area based on four algorithms

表4 融合后SAR圖像洪災區(qū)域的檢測結果Tab.4 Detection results of flooded area in the fusioned SAR image

圖30 RL濾波后4種算法的洪災區(qū)域檢測結果Fig.30 Detection results of flooded area based on four algorithms after RL filtering

表5 RL濾波后SAR圖像的洪災區(qū)域檢測結果Tab.5 Detection results of flooded area in the SAR image after RL filtering

本文提出的H-FCM算法綜合考慮洪災區(qū)域檢測的高準確率和高時效性的需求，既提升了SAR圖像的洪災區(qū)域檢測精度又縮短了檢測時間。H-FCM算法區(qū)別于傳統(tǒng)的無監(jiān)督聚類算法，沒有直接對待處理的圖像進行無監(jiān)督聚類，而是提出了一種反向聚類的方法。首先，H-FCM算法將洪災前光學圖像的河流位置檢測結果與洪災后的SAR圖像融合；其次，在融合圖像的最高級像素值區(qū)間對零像素值點進行分級聚類；最后，H-FCM算法在對融合圖像的洪災區(qū)域進行分級聚類的基礎上，將洪災前的河流位置作為空間約束，進一步篩除疑似洪災區(qū)域以提升檢測性能。

H-FCM算法將洪災前的光學圖像和洪災后的SAR圖像融合后進行洪災區(qū)域檢測，實現了全天侯、高準確率、高時效的洪災區(qū)域檢測要求。雖然H-FCM算法在一定程度上解決了SAR圖像相干斑噪聲給洪災區(qū)域檢測造成的問題，但是還有部分災區(qū)的細節(jié)沒有精確檢測到，這有待在未來的工作中解決。