基于雙邊加權核圖割的SAR圖像變化檢測

丁 稀，董 張 玉*，楊 學 志

(1.合肥工業大學計算機與信息學院/工業安全與應急技術安徽省重點實驗室，安徽 合肥 230601；2.合肥工業大學軟件學院/智能互聯網系統安徽省實驗室，安徽 合肥 230601)

0 引言

SAR(Synthetic Aperture Radar)圖像具有全天時、全天候監測特性，已成為遙感圖像變化檢測[1]、分割等應用中的主要數據源之一[2]，在災害監測、農業、軍事等領域具有極高的研究和應用價值[3]。SAR圖像變化檢測的核心步驟是對兩時刻SAR圖像生成的差異圖進行分析，從而獲取最終的變化二值圖[4]。由于差異圖分析可視為二分類問題，因此學者們多利用影像分割算法進行差異圖分析，如Bazi等對差異圖進行廣義高斯模型建模，并基于KI閾值算法選取最優分割閾值對差異圖進行二分類[5]。然而，閾值法易造成分割絕對化，對相干斑噪聲嚴重的SAR圖像難以保證分割精度。近年來，部分學者引入聚類法進行差異圖分析，獲取變化類和非變化類的兩個聚類中心，如Gong等在模糊聚類算法中引入MRF方法，根據局部區域中心像素和鄰域像素的標簽關系修正能量函數，有效減少了FCM算法的耗時，并在一定程度上削弱了散斑噪聲的影響[6]。此外，由于圖割(graph cuts)模型具有良好的空間上下文利用能力和影像分割能力，逐漸成為SAR圖像變化檢測研究的熱點。例如：Moser等將圖割算法與MRF相結合對差異圖進行分析，基于每個像素點與其周圍像素點的標簽類型極可能相同的假設，對差異圖的分布模型進行參數估計，并設計新的MRF能量函數，最后利用圖割對該能量函數進行優化，從而獲取變化檢測結果[7]；Salah等將核理論和圖割算法相結合提出核圖割(KGC)圖像分割算法，在差異圖上構建圖結構，利用核函數將差異圖轉到特征空間，使傳統圖割的分段常數模型更契合SAR圖像，最后設計能量函數完成圖像分割，從而提高差異圖的可分性[8]；Gong等進而提出基于局部擬合搜索模型的核圖割變化檢測方法，該模型逼近直方圖，可為核圖割提供優異的初始化結果，提升了圖割算法精度[9]。

上述算法在變化檢測任務中表現良好，但對空間上下文信息的挖掘程度不夠，僅從像素級角度考慮鄰接頂點的相似程度，并未從區域級角度考慮以像素為中心的區域塊之間的關系，導致差異圖分析過程中鄰域信息丟失。鑒于此，本文提出雙邊加權核圖割算法，即在均值比差異圖上構建雙邊加權核圖割模型，利用雙邊權重將鄰接像素點的像素級信息和區域級信息結合，解決圖像塊鄰域信息難以利用的問題，進而依托該雙邊加權核圖割模型，提出新的圖割能量函數，并利用能量最小化算法對該能量函數進行求解，獲取變化二值圖。

1 研究方法

本文提出的基于雙邊加權核圖割的SAR圖像變化檢測算法流程(圖1)為：1)將經過配準的同一地區不同時刻的兩幅SAR圖像X1、X2通過高斯濾波器進行降噪平滑預處理，以減少相干斑噪聲的影響；2)使用均值比差異算子[10]生成差異圖，利用鄰域信息降低噪聲對像素點的干擾，為后續差異圖分析提供初始分類結果；3)根據得到的差異圖構建雙邊加權核圖割模型，圖像各像素點為圖割模型頂點，圖像的像素級信息和區域級信息為圖割模型的邊權重；4)依托雙邊加權核圖割模型設計相應的圖割能量函數；5)依據網絡流理論獲取圖割模型的最大流/最小割，從而完成SAR圖像變化檢測任務。

圖1 基于雙邊加權核圖割的SAR圖像變化檢測流程Fig.1 Flowchart of SAR image change detection based on bilateral weighted kernel graph cuts

1.1 雙邊加權核圖割模型

圖割算法本質是將圖像中的每個像素與有限標簽集合(變化類和非變化類[11]二元標簽集合)進行一一配對。該算法將圖像建模成一個帶權無向圖，圖中頂點與圖像的像素點或關鍵特征點一一對應，圖的邊對應像素間按傳統四鄰域、八鄰域或K近鄰原則聯系起來的鄰接邊，而邊上的權重則反映某頂點的連接頂點對其影響力的大小[11]。本文根據生成的差異圖構建圖割模型：定義P為差異圖內所有像素點的集合，p為集合內某像素點，N為圖像內所有八鄰域像素點對(pn，pk)構成的集合，無向加權圖G={V，E，ω}(其中V為圖的頂點集，包含圖像內的所有像素點以及源點S(變化類標簽)和匯點T(非變化類標簽)(式(1))；E為圖的邊集，用于描述頂點間的鄰接關系，根據圖割理論，邊集可分為兩部分：一是鄰接像素點對N的n-links邊，二是像素點與源點和匯點連接的t-links邊(式(2))；ω為像素點間的相似度權重)。一般情況下，距離越近的像素彼此間的影響力越強，僅考慮單個像素的變化不具備物理意義。因此本文提出假設：空間最近的像素點對其中心像素點的影響最大，其影響不僅與自身強度相關，還與其所在的區域信息有關。依據以上假設，本文利用像素相似度權重ωI和區域相似度權重ωp計算ω(式(3))，其中，像素相似度權重ωI主要考慮兩頂點間的強度差異，常利用指數函數計算(式(4))。

V=P∪{S，T}

(1)

E=N∪{{p，S}，{p，T}}

(2)

ω=ωI·ωp

(3)

ωI(pn，pk)=exp{-|log[I(pn)]-log[I(pk)]|}

(4)

式中：I(pn)表示像素點pn的灰度值；pk表示像素點pn的八鄰接像素點。

利用圖像塊的結構信息可有效降低圖像噪聲，增加算法的魯棒性，已被應用于圖像去噪和分割等領域[12]。針對SAR圖像變化檢測過程中未充分利用圖像塊的區域級空間上下文信息問題，本文受非局部均值去噪算法[13]啟發，提出新的區域相似度公式，即采用區域局部均值的歸一化比率衡量兩個成對圖像塊的強度相似度，以減少SAR圖像相干斑噪聲的影響；然而均值信息作為圖像塊的特征信息有一定局限性，當兩個圖像塊均值相差不大時，均勻圖像塊和包含邊界的圖像塊之間的相似性很難通過均值歸一化比率衡量。為解決圖像塊邊緣信息利用困難的問題，本文將圖像塊的邊緣強度映射(Edge Strength Maps,ESM)矩陣[14]作為邊緣細節信息引入區域相似度特征度量中，以增強邊界信息精準分割能力。最終得到的區域相似度權重ωp的計算公式為：

(5)

(6)

式中：ESD(pn,pk)為像素點pn和pk的邊緣相似度距離，本文利用高斯伽馬窗算法[14]獲取兩時刻SAR圖像的邊緣信息，得到ESM矩陣，該矩陣各元素與像素點一一對應，值越大，對應像素點位于邊緣的可能性越大，反之該像素點處于均勻區域的可能性越大；h為邊緣相似度距離參數；d為區域塊的半徑；μ(pk)、μ(pn)分別為以像素點pk、pn為中心的區域塊的平均值；H(pn)、H(pk)分別為圖像ESM矩陣中以像素點pn、pk為中心的圖像塊，當H(pn)、H(pk)包含不同的邊緣信息時，其ESM矩陣差異越大，邊緣相似性距離值越大，當二者具有相同的邊緣或均處于均勻區域時，其邊緣相似性距離為0。

1.2 圖割能量函數

在雙邊加權核圖割模型構建完成后，兩時刻SAR圖像上的八鄰域像素點可被連接，其對應的相似度信息可通過鄰接邊獲取。依托以上信息，本文提出雙邊加權核圖割能量函數F(L)，計算公式為：

F(L)=D(L)+αR(L)

(7)

(8)

(9)

(10)

式中：D(L)為圖割能量函數的一元區域項，主要描述像素點屬于某類標簽L={Lc,Ln}的概率：Lc為變化類標簽，Ln為非變化類標簽，本文采用經典的分段常數圖割模型[11]，假設μl為分段常數模型參數,即源點、匯點的實際值，在給定觀測數據I(pn)下，使用參數μl與像素點pn的灰度值差異進行計算；R(L)為圖割能量函數的二元邊界項，用于衡量鄰接像素點間的相似度，并對非鄰接的像素點對進行不連續懲罰；α為常系數，用于平衡一元區域項和二元邊界項的比重，一般取0<α<1；δ(·)為狄拉克函數。

由于SAR圖像具有非線性特性，分段常數模型難以直接表示變化類與非變化類的分布，因此采用核函數將圖像轉換到特征空間，以此提高數據的可分性。根據Mercer定理[15]，核函數本質是數據在高維空間中的內積，因此，可以利用核函數計算能量函數的一元項。假設φ(·)是從數據觀測空間到更高維(可能是無限維)特征/映射空間的非線性映射，同時引入高斯徑向基核函數(RBF)，即：K(i，j)=exp(-‖i-j‖/2σ2)，則在核誘導空間中能量函數計算公式為：

(11)

‖φ(I(pn))-φ(μ)‖2=K(I(pn),I(pn))+
K(μ,μ)-2K(I(pn),μ)

(12)

該能量函數的一元項由核函數隱式映射到特征空間，二元去斑鄰域項綜合衡量像素相似度和區域相似度，并在區域級信息引入邊緣相似度特征。至此，雙邊加權核圖割能量函數建立完畢，其最佳切割方法可依據α-β移動交換標簽優化算法[11]求得。

圖2展示了圖割模型結構和算法處理過程：不同顏色的圓圈表示不同灰度的像素點，點S和點T分別表示變化類標簽和非變化類標簽，每一頂點均與兩終端頂點S和T連接，形成t-links邊，表示為虛箭頭，該邊上的權重反映各像素點屬于兩類標簽的概率；鄰域內的頂點相連，形成n-links邊，表示為實線，邊上的權重反映像素點間的相似度。當圖割模型建立后，為完成圖像分割任務，需尋找一條割線(Cut，由虛線表示)將圖中頂點分隔為兩部分，而割線所包含的邊稱為割集；當割集所在邊的權重之和最小時，即為圖像的最終分割結果，該割線被稱為最小割[16]。

圖2 圖割算法示意Fig.2 Schematic diagram of graph cuts

2 實驗結果與分析

2.1 實驗數據與參數設置

實驗選取Ottawa和Mexico fire兩組真實SAR數據，且均已經過預處理并達到亞像素級的配準精度(表1)。算法參數中，高斯濾波參數如果過大，會加深濾波程度，導致檢測結果出現過平滑現象，如果過小，檢測結果零碎點過多，抗噪性減弱，故本文選取常見的高斯濾波參數0.5；鄰域加權斑塊尺寸過大會導致算法復雜度變高和出現過平滑現象，所以本文選擇較小的斑塊尺寸d=3，在提高算法運行效率的同時，能保證良好的平滑降噪性能；邊緣相似度參數h在合理范圍內使邊緣距離權重由邊緣相似性距離值決定，過大或過小都會使該權重失去意義，本文設置h=0.002。

表1 SAR圖像數據的先驗信息Table 1 Prior information of SAR image dataset

選取基于修正馬爾科夫隨機場的模糊C均值算法[6](MRFFCM)、基于核的圖割算法(Kernel-induced Graph Cuts,KGC)、PCA-Kmeans算法[17]和FLICM聚類算法[18]作為對比實驗組，采用變化檢測常用的指標[4]對各算法精度進行評價,包括：1)錯檢率(FPR)：變化類型的像素數占實際非變化類型像素總數的比例；2)漏檢率(FNR)：遺漏的變化類型像素數占實際變化類像素總數的比例；3)總體精度(OA)：檢測正確像素數占圖像像素總數的比例；4)Kappa系數：綜合反映變化檢測的二值結果圖與實際參考圖的趨近程度。FPR和FNR越小，OA和Kappa系數越接近1，表明檢測結果精度越高、越接近參考圖[19]。

2.2 實驗對比分析

2.2.1 Ottawa數據集變化檢測結果分析 第一組數據集呈現了Ottawa地區洪水前后的真實地物變化情況(圖3a、圖3b)，從真實變化參考圖(圖3c)可以看出，變化區域(圖中白色區域)多為線狀或帶狀，且塊狀變化區域內有很多碎點和不規則邊緣，因此變化檢測算法要具有良好的抗噪性和準確的邊緣定位能力。由5種算法檢測結果(圖3d-圖3h)對比可知，PCA-Kmeans算法檢測結果良好，但在中上部區域漏檢點較多；FLICM算法碎點較多，中上部區域也存在較多漏檢點，對變化區域和非變化區域表征均不理想；MRFFCM算法在非變化區域(圖中黑色區域)存在大量碎點，抗噪性能不佳；KGC算法的噪聲抑制能力較好，但變化區域和非變化區域邊界模糊，可分性較差；本文算法合理利用了像素級和區域級空間上下文信息，抗噪性能較好且變化區域邊緣定位準確。

圖3 Ottawa數據集變化檢測結果Fig.3 Change detection results of the Ottawa flood dataset

由5種算法在Ottawa數據集上的定量檢測結果(表2)可以看出，PCA-Kmeans算法的錯檢率較低，但漏檢率較高，影響了檢測精度；FLICM和MRFFCM算法的錯檢率較高，表明其抗噪能力不強；KGC算法錯檢率較低，但漏檢率較高，二值變化圖出現過平滑現象，導致最終的檢測結果不理想；本文算法OA和Kappa系數均最高，驗證了算法的有效性。

表2 不同算法Ottawa數據集變化檢測結果比較Table 2 Comparison of change detection results of the Ottawa flood dataset for various algorithms

2.2.2 Mexico fire數據集變化檢測結果分析 第二組數據集呈現了Mexico地區火災前后的真實地物變化情況(圖4a、圖4b)，與前文分析類似，PCA-Kmeans、FLICM算法檢測結果中變化區域與非變化區域邊界表現良好，但在較小變化區域表現不理想，漏檢率較高；MRFFCM和KGC算法背景區域表現良好，但變化區域邊界較平滑，且丟失部分較小變化區域，致使漏檢率較高；本文算法(圖4d)確保漏檢較少，同時對變化區域邊緣定位準確。由表3可以看出，本文算法的漏檢率最低，總體精度OA和Kappa系數最高，進而驗證了本文算法在檢測精度和邊緣檢測能力上的優異性。

表3 不同算法Mexico fire數據集變化檢測結果比較Table 3 Comparison of change detection results of the Mexico fire dataset for various algorithms

圖4 Mexico fire數據集變化檢測結果Fig.4 Change detection results of the Mexico fire dataset

3 結論

本文針對SAR圖像固有的相干斑噪聲問題，提出了雙邊加權核圖割算法，該算法特點為：1)圖割結構上，從像素級角度和區域級角度挖掘圖像的空間上下文信息，并將兩者結合作為圖割的雙邊權重二元項，提升了抗噪能力；2)能量函數上，在去斑鄰域項中引入邊緣相似度特征，通過邊緣強度映射矩陣獲取成對斑塊的邊緣相似距離，緩解圖割模型邊緣細節信息丟失嚴重問題。與其他4種算法對比結果表明，本文算法總體精度和Kappa系數更高，邊緣細節信息更豐富，且具有良好的抗噪能力。本文提出的SAR圖像變化檢測算法，主要利用像素的強度特征和邊緣特征，尚未研究其他特征信息，今后將進一步挖掘圖像多種特征，并將其引入雙邊核圖割模型中，探究多種特征對變化檢測的影響。