基于融合差異圖的變化檢測方法及其在洪災中的應用

黃平平 段盈宏 譚維賢 徐 偉

(內蒙古工業大學信息工程學院 呼和浩特 010051)

(內蒙古自治區雷達技術與應用重點實驗室 呼和浩特 010051)

1 引言

由于惡劣的氣候和天氣影響自然災害頻發，其中，洪災嚴重破壞了生態環境，損害了私人和公共財產，對人類社會的生存發展有重大影響[1]。因此，對洪害前后進行變化檢測，確定洪災的淹沒范圍，對災情的評估和災后的管理規劃都有一定的指導作用。合成孔徑雷達(Synthetic Aperture Radar,SAR)可以全天時、全天候地獲取遙感數據，不受光照和天氣的影響，彌補了光學和紅外遙感的不足。隨著遙感技術的發展，不同時相、不同頻帶以及不同極化方式的多種SAR圖像資源為SAR圖像變化檢測提供了數據基礎，并廣泛應用在災區定位和災害評估等方面[2]。

SAR圖像的變化檢測方法一般可分為兩種[3]：監督和非監督。監督類變化檢測方法需要大量的先驗知識，而實際中往往缺乏真實的參考信息，適用范圍較窄。非監督類變化檢測方法不需要先驗信息，可以直接對獲得的兩時相圖像進行變化檢測，在實際應用中更為廣泛，且非監督類變化檢測技術多在差異圖生成和變化信息提取兩方面進行創新。傳統方法多使用不同的變化檢測算子構建差異圖，如經典的差值法和比值法[4,5]。相對熵，又被稱為KL散度(Kullback-Leibler divergence)或信息散度，是兩個概率分布間差異的非對稱性度量，以此構造的差異圖可衡量兩幅影像間的相似程度。考慮到散斑噪聲影響以及檢測結果的精度，Gong等人[6]利用鄰域像元灰度值和中心像元灰度值的相關性，提出一種基于鄰域比值算子的差異圖生成方法。單一類型的差異圖像有自身的局限性，如均值比差異圖對目標區域的群分性較好，但對紋理信息感知較弱，且存有大量噪聲，對數比差異圖在一定程度上削弱了背景部分的斑點噪聲，但殘留的噪聲仍對后續處理有很大的影響[7–9]。基于此，一些研究人員選擇使用圖像融合的方法，利用不同差異圖像的互補信息以及不同傳感器的并行信息，在城市化和河口區域建筑的變化檢測中展現了較好的結果[10–12]。在不同雷達圖像數據中，通過融合多時相雷達數據以及相干和非相干性分析，對氣旋、滑坡等災害區域的評估也有一定的效果[13–15]。小波變換相對于輪廓波等多尺度變換而言，具有更小的計算復雜度，同時，小波變換可以將圖像的時空信息轉換為頻域信息，進而可以方便地提取圖像的詳細信息，在圖像融合中有廣泛的應用[16–18]。在變化信息提取方面，Sumaiya等人[19]使用加權閾值分割對數均值比差異圖，獲得變化檢測結果，計算過程快速簡便。Zhang等人[20]利用廣義高斯分布對差異圖像的每一類進行建模，采用圖割算法提取空間先驗信息并通過模糊C均值算法很好地初始化模型的參數，在最大后驗的貝葉斯準則下提取變化信息，其精度與算法速度都有一定的提升。Li等人[21]提出了一種稱為SDPCANet的監督學習網絡模型，該模型將顯著性檢測與主成分分析網絡相結合，并在多時相SAR圖像中驗證了方法的有效性。Li等人[22]提出了一種基于卷積神經網絡(Convolutional Neural Networks,CNN)的SAR圖像變化檢測新方法，該方法通過CNN直接從原始的兩個SAR圖像中生成分類結果，無需進行任何預處理操作，消除了生成差異圖像的過程，從而減少了差異圖對源圖像的影響，并且在異質圖像的變化檢測中也有較好的效果。

洪災發生時往往伴隨著多云多雨，而光學傳感器無法穿透云層，無法獲取有效的地面信息，因此，SAR圖像的變化檢測更有效地應用在洪災區域中。由于洪災發生時降雨影響了地物散射特性[23]，在整幅圖像中許多沒有受洪災影響的區域也呈現為偽暗區。在傳統的差異圖中，偽暗區與背景信息差異較大，而與變化信息更為接近，一般的變化檢測算法更容易將偽暗區作為變化信息處理，這也導致了檢測結果中較高的虛警率和較低的檢測精度[24]。由于偽暗區的存在，對差異圖的直接聚類易導致較多的錯誤信息，故可利用皮爾遜相關系數對模糊局部信息C均值(Fuzzy Local Information C-Means,FLICM)聚類方法的初始聚類結果進行二次分類，減小對檢測結果的影響。一般而言，洪水淹沒區域的相關系數更低，而標準未變化區域的相關系數更高，所以通過比較待定區域的相關系數，待定區域的類別，更容易劃分為變化區域和未變化區域中兩者相關系數均值更小的類別[25–28]。

本文結構安排如下：本文第2部分提出了一種利用融合差異圖的SAR圖像變化檢測方法，然后在此差異圖下得到基于相關系數FLICM算法二次分類后的檢測結果，最后經迭代條件模型和馬爾科夫隨機場(Iterative Conditional Model-Markov Random Field,ICM-MRF)算法優化得到變化檢測最終結果；第3部分首先利用Bern地區的ERS-2實測SAR數據以及加拿大Ottawa地區的Radarsat實測SAR數據進行結果分析并驗證本文方法的有效性，之后利用該方法對中國鄱陽湖區域2020年6月和7月洪災前后的Sentinel-1-A實測SAR數據進行分析，進而對受災面積和其后的災害趨勢進行評估；第4部分對全文做出總結。

2 SAR圖像變化檢測算法

本文算法流程主要由融合差異圖構造和ICMMRF分割[29,30]兩部分組成，算法流程如圖1所示。

融合差異圖構造：首先得到兩幅SAR圖像的均值比差異圖和利用鄰域異質性信息的相對熵差異圖，然后經小波變換，通過加權平均的方式將兩幅圖像的低頻系數重構，與各自的高頻系數，經小波逆變換得到兩幅融合低頻系數后的差異圖，通過局部能量法融合為最終的差異圖。

ICM-MRF分割：在ICM-MRF算法中，使用FLICM算法將變化檢測結果分為3類，分別為變化區域、未變化區域和待定區域，對待定區域通過洪災前后圖像的皮爾遜相關系數(后簡稱為相關系數)進一步分為變化區域和未變化區域，并以此作為初始分割，得到最終的變化檢測結果。

2.1 融合差異圖構造

假設X1,X2為已經過配準預處理的兩時像SAR洪災前后變化區域圖像，則離散圖像數據X2對X1的相對熵可表示為

其中，N表示像元總個數，xn表示第n個像元的值，Ys表示X2對X1的相對熵。

圖1 基于融合差異圖的變化檢測算法流程圖Fig.1 Flow chart of change detection algorithm in this paper

在相對熵的計算中引入像元空間鄰域的異質性信息，突出圖像邊緣，增強連續性，同時抑制背景信息。鄰域異質性信息可表示為

?表示圖像像元間的異質度，其定義為像元鄰域的方差δ(x)和其鄰域的均值μ(x)的比值。則引入鄰域異質性信息后的相對熵計算公式可定義為

其中，Yd為X2對X1的鄰域異質性相對熵，?n1,?n2分別為圖像X1,X2同一位置像元對應的異質性系數歸一化后的結果，μ1,μ2為其分別對應的去除中心像元點的鄰域均值

式(4)為相對熵計算的對稱性處理，Yk為本文改進的相對熵差異圖的最終結果。使用鄰域異質信息來調節圖像相對熵的計算，熵值越小表示兩幅圖像同一位置對應像元的相似程度越高，反之則相似程度越低。

引入鄰域異質性信息的相對熵計算方法，相對于傳統的單一像元值的熵值計算方法，有較好的背景抑制和邊緣連續性效果，但目標區域的完整性和平滑性一般，故將目標區域完整性和平滑性較好的均值比算子差異圖與本文引入鄰域異質性信息的差異圖，通過小波變換進行融合。均值比算子可定義為

其中，μ1(n)和μ2(n)分別表示在圖像X1和X2第n個像元其鄰域(一般取3×3大小)內的平均值，Ym表示均值比差異圖。

本文圖像通過小波變換進行融合，小波變換一般先將圖像分解為更低分辨率水平上的低頻輪廓信息，以及在水平垂直和對角線方向的高頻細節信息，再經過小波系數的替換、選擇權值或者疊加運算進行融合。本文選擇的融合規則如下：(1)取兩圖像小波分解后的加權低頻系數均值作為新的低頻系數；(2)將新的低頻系數與各自的高頻系數經小波逆變換得到融合低頻系數后的差異圖；(3)將兩幅新的差異圖通過局部能量法得到目標差異圖。經過本文規則得到的差異圖中，具有均值比差異圖中較好的區域完整性和平滑性，也具有改進相對熵差異圖中背景抑制性較強、邊緣連續性較好的優點。

2.2 ICM-MRF分割

FLICM是一種無監督的模糊聚類算法，通過圖像的統計特性可以較好地完成分類，進而提取圖像的變化信息。使用FLICM算法將融合差異圖分為兩類的目標函數為

若以μn1和μn2分別表示為第n個像元對于兩類的隸屬度，通過比較兩類的隸屬度大小，則可以確定初始分割結果。

通過FLICM算法將差異圖分為3類，分別代表變化區域、待定區域以及未變化區域。計算洪災前后兩SAR圖像的相關系數，如式(10)所示

式中，Ic(i,j)表示為(i,j)像元點對應的相關系數值；X1(i,j)和X2(i,j)均表示圖像的(i,j)像元點對應的像元值；表示X(i,j)像元點周圍鄰域像元值的均值。計算各類區域相關系數均值，比較待定區域相關系數與變化區域和未變化區域相關系數的均值，進行二次分類。

融合差異圖通過以上算法二次分類后，其結果可作為ICM-MRF算法的初始分割。假設融合差異圖像為Y={yη=yij}，初始分割圖像X={xη=xij}為定義在Y上的馬爾科夫隨機場，是對圖像所做的標記；xη是點η的類別標記，xη ∈I(I=1,2,···,k),k為類別數，本文中k=2。設點η的鄰域值記為sη,s={sη,η ∈Y}是η的鄰域系統。求解最優圖像標記問題的關鍵在于求解先驗概率P(X)和概率密度函數P(Y |X)。基于貝葉斯決策理論以及馬爾科夫隨機場和吉布斯分布的等價性，選取2階鄰域系統和合適的模型，用β表示馬爾科夫參數，則局部條件先驗概率可表示為

Vh為使用Potts數學模型的基團勢能函數，H為所有基團的集合。由于要得到后驗概率P(X|Y)最大，即等效于求解U(X|Y)最小，采用迭代條件模型(Iterative Conditional Model,ICM)求解，它的核心思想是用逐像元點最大化條件概率，迭代地實現每個像元點的更新，如式(15)所示

在馬爾科夫隨機場(Markov Random Field,MRF)算法中，作為初始分割的標簽圖對MRF的分割結果有較大的影響，因此將利用圖像相關性對FLICM算法二次分類后的結果作為MRF算法的初始標簽圖，以期得到更好的初始分割結果，則基于融合差異圖的ICM-MRF算法如表1所示。

3 實驗結果與分析

本文選擇的驗證數據為瑞士Bern地區在1999年4月和5月的ERS-2數據，兩時相遙感圖像為圖2(a)、圖2(b)，大小為301×301，變化檢測參考圖為圖2(c)；加拿大Ottawa地區在1997年5月和8月的Radarsat遙感數據，兩時相遙感圖像為圖3(a)、圖3(b)，大小為290×350。兩數據集產生圖像變化的主要原因均為洪水淹沒陸地所致[2]。

表1 基于融合差異圖的ICM-MRF圖像分割計算過程Tab.1 ICM-MRF image segmentation calculation process based on fusion difference map

使用上述數據進行仿真試驗，通過對Bern地區數據以及Ottawa地區數據均值比、相對熵方法得到的差異圖，以及本文融合差異圖進行比較，得到了上述所有方法在兩數據集下不同差異圖的期望最大化(Expectation-Maximization,EM)算法、K-means算法、FLICM算法以及本文ICM-MRF算法的檢測結果，并分別在兩數據集下進行主觀和客觀的評價。其后通過本文方法對2020年中國鄱陽湖地區數據進行變化檢測，進而估計由于鄱陽湖湖水覆蓋而受災的區域面積和變化趨勢。

本文采用虛警率(False Positive,FP)、總錯誤率(Overall Error,OE)、準確率(Probability Correct Classification,PCC)和Kappa系數作為量化評估的指標。其中，PCC表示總體的分類正確性，Kappa系數綜合考慮了像元點的檢測狀況，是衡量檢測精度的關鍵指標，其值越大表示檢測精度越高。

3.1 差異圖2分類和3分類影響分析

Bern地區數據集的不同算法差異圖如圖4(a)—圖4(c)所示，并以A表示均值比差異圖，B表示相對熵差異圖，C表示本文差異圖，在此3種差異圖下，以EM算法、K-means算法和FLICM算法為例，分析了2分類和3分類后基于相關系數二次分類對結果的影響，如圖5、圖6和圖7所示。對比各算法的差異圖可得，均值比算法在背景中存在大量噪聲，相對而言，相對熵法雖大量抑制了噪聲，但模糊了目標區域的結構。本文構造的差異圖綜合了相對熵和均值比差異圖的優勢，邊緣連續性和區域完整性較相對較好，與背景的對比度較兩差異圖更為明顯，且噪聲點較均值比差異圖更少。從以下變化檢測結果圖中不難看出，當僅將差異圖劃分為2類時，噪聲點和一些區域也將直接定義為變化區域或者為變化區域，直接影響了結果的區域劃分和檢測精度。而將差異圖先劃分為變化區域，未變化區域以待定區域3類區域，再將其根據不同劃分區域的相關系數二次分類后，將有效地減小噪聲的影響，并加強了區域的完整性，進而提升檢測結果的精度。相對于均值比和異質相對熵差異圖，在不同算法下，本文差異圖也顯示了更好的分類效果。

圖2 Bern地區數據Fig.2 Bern region data

圖4 Bern數據差異圖Fig.4 Bern data difference map

圖5 基于均值比差異圖2分類和3分類后二次分類結果Fig.5 Based on the difference map in mean ratio 2 and 3 classification results after secondary classification

圖6 基于相對熵差異圖2分類和3分類后二次分類結果Fig.6 Based on the difference map in relative entropy 2 and 3 classification results after secondary classification

3.2 Bern和Ottawa地區實驗結果與分析

Bern地區數據集對應的每種差異圖的不同算法變化檢測結果的相關指標分析如表2所示。對于不同差異圖Bern數據集的本文算法變化檢測結果如圖8(a)—圖8(c)所示。檢測結果中白色點為變化區域的像元點。從變化檢測結果圖中也可看出，在本文ICM-MRF算法下，融合差異圖的檢測結果顯示了更好的變化區域邊界，與實際的變化情況更為接近。表2中也顯示了針對不同差異圖下同一種算法的變化檢測結果，本文差異圖相對于其他差異圖而言，具較低的虛警率與總錯誤率，更高的正確率和Kappa系數。針對本文差異圖而言，由于EM算法檢測結果中存在較多噪聲，故Kappa系數較低，但仍可看清變化區域的完整結構，并且本文ICM-MRF算法具有更高的Kappa系數和較低的錯誤率，相對于檢測效果較好的FLICM算法，Kappa系數提高至0.837。

圖7 基于融合差異圖2分類和3分類后二次分類結果Fig.7 Based on the fusion difference and the results of secondary classification after classification 2 and 3

表2 Bern數據集不同差異圖不同算法變化檢測指標分析Tab.2 Analysis of change detection indicators of different algorithms in Bern dataset

圖8 Bern數據本文ICM-MRF算法變化檢測結果Fig.8 Change detection results of this paper algorithm in Bern data

圖9 Ottawa數據差異圖及其變化檢測結果Fig.9 Ottawa data difference map change detection result

Ottawa地區數據集的不同算法差異圖如圖9(a)—圖9(c)所示，并以A表示均值比差異圖，B表示相對熵差異圖，C表示本文差異圖，其對應的每種差異圖的不同算法變化檢測結果的相關指標分析如表3所示。對于不同差異圖Ottawa數據集的本文ICM-MRF算法變化檢測結果如圖9(d)—圖9(f)所示，EM算法變化檢測結果如圖9(g)—圖9(i)所示，K-means算法變化檢測結果為圖9(j)—圖9(l)所示，FLICM算法變化檢測結果為圖9(m)—圖9(o)所示。檢測結果中白色點為變化區域的像元點。通過對比各算法的差異圖不難看出，均值比算法、相對熵算法與Bern區域結果分析一致，且本文的差異圖構造方法在Ottawa數據集內也體現出邊緣連續性，區域完整性以及背景對比度均相對較好的效果。不同差異圖的變化檢測結果中，構造的融合差異圖在本文ICM-MRF算法下變化檢測結果相對于均值比差異圖的檢測結果具有更好的細節表現，而相對于相對熵差異圖的檢測結果具有更好的區域完整性，且Kappa系數也有一定提升。表3中也顯示了與Bern數據集基本相似的指標分析結果，由此可見本文算法在兩數據集均具有較好的適用性。整體而言，本文算法可以提高變化檢測結果的精度，減小錯誤率。

3.3 鄱陽湖區域洪災前后變化檢測結果分析

由于連續多日暴雨的影響，中國鄱陽湖區域發生洪災，淹沒了周邊大片陸地區域，人工無法有效地實時實地測量受災區域而制訂善后處理計劃，利用SAR衛星遙感圖像，可以快速準確地提取受災區域，評估災區的空間分布和動態變化。本文選取中國鄱陽湖地區2020年6月和7月Sentinel-1-A數據進行實驗來評估此次洪災。由于完整的鄱陽湖區域分布在Sentinel-1-A數據的兩景圖像內，故本文使用對兩景數據分別檢測的方法進行實驗。具體實驗數據為2020年6月26日數據，拼合后完整鄱陽湖區域圖像如圖10(a)所示，7月8日數據如圖10(b)所示，7月20日數據如圖10(c)所示，3日包含鄱陽湖整體區域的Sentinel-1-A數據。圖10(a)、圖10(b)和圖10(c)大小均為999×651，像素大小為243×243 m2，產生圖像變化的主要原因是鄱陽湖水淹沒附近區域陸地所致。

通過EM算法、FLICM算法、K-means算法和本文ICM-MRF算法，在融合差異圖下對鄱陽湖Sentinel-1-A衛星遙感圖像進行變化檢測，EM算法變化檢測結果如圖11(a)、圖11(b)所示，FLICM算法變化檢測結果如圖11(c)、圖11(d)所示，K-means算法變化檢測結果如圖11(e)、圖11(f)所示，本文ICM-MRF算法變化檢測結果如圖11(g)、圖11(h)所示。通過對比6月26日與7月8日的變化檢測結果，明顯看到在本文差異圖下兩種算法均有較好的檢測結果，得到的變化區域也較為相似，但不同的是，在圖11(c)、圖11(d)的FLICM算法中，黃色方框區域為虛警區域，此區域在真實雷達圖中并未被洪水淹沒，圖11(a)、圖11(b)EM算法和圖11(e)、圖11(f)K-means算法中，黃色方框區域為漏警區域，此區域為被洪水淹沒的區域，而本文算法均未將兩區域錯誤劃分，整體來看，本文算法下，變化區域結構更為完整，檢測結果更為準確。

基于本文算法，圖11(g)—圖11(j)變化檢測結果為3幅鄱陽湖洪災前后變化檢測結果。其中，圖11(g)為6月26日與7月8日的變化檢測結果，圖11(h)為其對應在實際中的變化區域。圖11(i)為7月8日與7月20日的變化檢測結果，圖11(j)為其對應在實際中的變化區域。白色、紅色和藍色像元點代表變化區域。由此可見，本文算法較完整且清晰地檢測出鄱陽湖區域的洪災前后變化范圍。選取的包含完整鄱陽湖區域的SAR圖像實際面積約為38402.46 km2，由于受暴雨影響，鄱陽湖水位上升，7月8日—7月20日，鄱陽湖西南部區域受災嚴重，湖水淹沒了湖周邊大片陸地，包括附近水域淹沒的水中小島和陸地等。由表4鄱陽湖區域變化檢測結果量化分析可知，此期間變化像素數為12589，占總像素點的1.94%，總體受災面積大約為 743.37 km2。7月8日—7月20日，在表4中仍可看到5475個變化像素點，占總像素點的0.84%，這表明總體受災面積仍在提升，但變化像素點減少了大約56.5%，洪水淹沒區域較之前明顯減少，期間總受災面積大約為323.29 km2，洪災擴散勢頭逐漸減小。值得注意的是，圖11(j)中藍色區域雖為變化區域，但不是洪水淹沒區域，而是由水體轉為陸地的區域，這也說明部分區域的洪水被引流或者逐漸退去，總體形趨于穩定。

表3 Ottawa數據集不同差異圖不同算法變化檢測指標分析Tab.3 Analysis of change detection indicators of different algorithms in Ottawa dataset

圖10 鄱陽湖地區數據Fig.10 Poyang Lake area data

圖11 鄱陽湖區域EM,FLICM,K-means和本文算法洪災前后變化檢測結果Fig.11 Change detection results of EM,FLICM,K-means and this paper algorithm before and after the flood disaster in Poyang Lake area

表4 鄱陽湖區域變化檢測結果量化分析Tab.4 Quantitative analysis of the detection results of Poyang Lake area change

4 結論

針對洪災區域的變化檢測，本文提出一種基于融合差異圖的變化檢測方法，利用鄰域異質性信息保持相對熵差異圖中變化區域的連續性，抑制背景信息，并通過小波變換，構造與均值比差異圖融合的最終差異圖。通過FLICM聚類算法與圖像的相關系數結合而得到更接近于洪災變化區域的二次分類結果，然后將其經過ICM-MRF算法優化得到最終的變化檢測結果。在Bern地區ERS-2遙感數據以及Ottawa地區Radarsat遙感數據實驗結果均表明，本文方法對洪災區域的變化檢測結果有較好的檢測精度，變化區域結構更為完整準確。通過本文算法得到2020年6,7月份中國鄱陽湖地區的3景Sentinel-1-A遙感圖像變化檢測結果，并對其受災區域面積以及變化趨勢進行估計，具有一定的實用價值。今后將繼續研究洪災區域產生虛警、漏警的原因以及在變化檢測中自動區別洪水淹沒區域和水體轉為陸地區域的方法。