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

王志豪 李 剛 蔣 驍

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

1 引言

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

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

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

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

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

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

2 洪災區域檢測方法的發展

2.1 傳統的洪災區域檢測方法

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

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

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

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

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

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

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圖像洪災區域的檢測時間、降低洪災區域的錯誤檢測概率以及抑制SAR圖像中存在的大量隨機分布的相干斑噪聲對檢測結果的影響。根據河流連續體(River Continuum Concept,RCC)[21]的概念，在已知的河流系統內，自上游到下游的物理變量的梯度變化是連續的[22]。因此發生洪水災害之后，若地理環境沒有發生劇烈變化，洪水的形態仍然符合RCC。并且，根據典型地物的光譜特征曲線圖[23](見圖3)，水體對近紅外波段的電磁波吸收能力比大部分地物強。

表1 FCM算法Tab.1 FCM algorithm

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

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

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

圖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圖像洪災區域檢測的改進方法

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

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

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

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

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

3.2.2 檢測河流位置的區域生長算法規則

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

為了充分消除區塊噪聲的影響，并保留河流區域的完整性，步驟(3)的卷積操作分為兩步。首先，確定大小為α×α的全1卷積核Wα，將Pt與Wα卷積得到初步消除區塊噪聲的結果P1，如式(9)和式(10)；其次，將P1與Wβ卷積得到基本消除區塊噪聲的結果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 區域生長算法的參數選擇

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

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

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

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

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

3.2.4 利用區域生長算法獲得河流位置

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

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

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

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

將洪災前的河流位置Eb和洪災后的SAR圖像Ea融合。融合規則為將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

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

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

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

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

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

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

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

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

3.4 洪災區域檢測的空間約束

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 實驗結果與分析

4.1 檢測結果的評價指標

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

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

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

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

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

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

(4) 卡帕系數

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

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

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

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

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

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

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

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

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

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

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

圖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算法獲得的洪災區域檢測結果與真值圖(圖23)比較，得到4種算法的檢測結果與真實災區的比較結果(見圖26)。在圖26中，白色部分表示實際洪災區域被檢測為洪災區域；綠色部分表示實際洪災區域被檢測為非洪災區域；紅色區域表示實際非洪災區域被檢測為洪災區域；黑色部分表示實際非洪災區域被檢測為非洪災區域。本文后續實驗的檢測比較結果與圖26的表示方法一致。

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

4.3.2 實驗2

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

實驗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種方法對降噪處理后的圖片進行洪災區域檢測。定義新的檢測方法分別為RL_H-FCM算法、RL_WA算法、RL_Snake算法、RL_H-Kmeans算法，其檢測結果如圖30所示，檢測性能評價指標如表5所示。

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

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

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

5 結束語

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

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

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

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

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

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

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