基于雙邊加權(quán)核圖割的SAR圖像變化檢測(cè)

丁 稀，董 張 玉*，楊 學(xué) 志

(1.合肥工業(yè)大學(xué)計(jì)算機(jī)與信息學(xué)院/工業(yè)安全與應(yīng)急技術(shù)安徽省重點(diǎn)實(shí)驗(yàn)室，安徽 合肥 230601；2.合肥工業(yè)大學(xué)軟件學(xué)院/智能互聯(lián)網(wǎng)系統(tǒng)安徽省實(shí)驗(yàn)室，安徽 合肥 230601)

0 引言

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

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

1 研究方法

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

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

1.1 雙邊加權(quán)核圖割模型

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

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

(5)

(6)

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

1.2 圖割能量函數(shù)

在雙邊加權(quán)核圖割模型構(gòu)建完成后，兩時(shí)刻SAR圖像上的八鄰域像素點(diǎn)可被連接，其對(duì)應(yīng)的相似度信息可通過(guò)鄰接邊獲取。依托以上信息，本文提出雙邊加權(quán)核圖割能量函數(shù)F(L)，計(jì)算公式為：

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

(7)

(8)

(9)

(10)

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

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

(11)

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

(12)

該能量函數(shù)的一元項(xiàng)由核函數(shù)隱式映射到特征空間，二元去斑鄰域項(xiàng)綜合衡量像素相似度和區(qū)域相似度，并在區(qū)域級(jí)信息引入邊緣相似度特征。至此，雙邊加權(quán)核圖割能量函數(shù)建立完畢，其最佳切割方法可依據(jù)α-β移動(dòng)交換標(biāo)簽優(yōu)化算法[11]求得。

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

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

2 實(shí)驗(yàn)結(jié)果與分析

2.1 實(shí)驗(yàn)數(shù)據(jù)與參數(shù)設(shè)置

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

表1 SAR圖像數(shù)據(jù)的先驗(yàn)信息Table 1 Prior information of SAR image dataset

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

2.2 實(shí)驗(yàn)對(duì)比分析

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

圖3 Ottawa數(shù)據(jù)集變化檢測(cè)結(jié)果Fig.3 Change detection results of the Ottawa flood dataset

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

表2 不同算法Ottawa數(shù)據(jù)集變化檢測(cè)結(jié)果比較Table 2 Comparison of change detection results of the Ottawa flood dataset for various algorithms

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

表3 不同算法Mexico fire數(shù)據(jù)集變化檢測(cè)結(jié)果比較Table 3 Comparison of change detection results of the Mexico fire dataset for various algorithms

圖4 Mexico fire數(shù)據(jù)集變化檢測(cè)結(jié)果Fig.4 Change detection results of the Mexico fire dataset

3 結(jié)論

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