邊緣懲罰層次區(qū)域合并SAR圖像分割算法

張澤均 水鵬朗

1 引言

最近，合成孔徑雷達(dá)(SAR)成像技術(shù)得到快速發(fā)展[1]和廣泛應(yīng)用[2]。SAR圖像分割是其處理和應(yīng)用中的關(guān)鍵步驟[3]，它將一幅SAR圖像分割成互不相交的同質(zhì)區(qū)域。圖像中的相干斑噪聲增大了它的分割難度[4]。兩大類SAR圖像分割算法是：基于區(qū)域特征和區(qū)域合并的算法[5,6]和基于優(yōu)化模型的算法[3,7 11]- 。

第1類算法以初始過分割為基礎(chǔ)，設(shè)計(jì)相鄰區(qū)域之間的相似性度量，利用區(qū)域合并技術(shù)獲得最終分割結(jié)果。文獻(xiàn)[6]利用超像素技術(shù)獲得圖像的初始分割，同時利用Gabor算子和Prewitt算子提取圖像的區(qū)域和邊緣特征，設(shè)計(jì)相鄰區(qū)域之間的相似性度量，該方法增大了SAR圖像中乘性斑點(diǎn)噪聲對分割結(jié)果的影響。

第2類算法是基于模型優(yōu)化的方法。首先建立圖像分割模型，然后對模型優(yōu)化求解實(shí)現(xiàn)圖像分割[3,711]-。文獻(xiàn)[7]基于最短描述長度(Minimum Description Length, MDL)準(zhǔn)則建立SAR圖像分割模型，遞歸地變形多邊形網(wǎng)格實(shí)現(xiàn)模型優(yōu)化求解。該方法的缺點(diǎn)是模型優(yōu)化的時間復(fù)雜度高，且算法性能和效率對初始分割敏感[7]。文獻(xiàn)[11]將 SAR 圖像的邊緣信息融入高斯馬爾科夫隨機(jī)場模型的懲罰項(xiàng)中，建立圖像分割模型，提出一種基于語義迭代區(qū)域生長(Iterative Region Growing using Semantics, IRGS)的模型優(yōu)化方法[11]。不足之處：(1)高斯分布不適合對強(qiáng)度格式 SAR圖像建模[12]；(2)基于梯度算子的邊緣強(qiáng)度提取方法不具備的恒虛警(Constant False-Alarm Rate, CFAR)特性[13]，梯度算子對乘性噪聲的敏感性使得獲得的邊緣具有明顯的鋸齒狀；(3)算法中用來控制欠分割的區(qū)域標(biāo)記算法的時間復(fù)雜度高。

本文利用 gamma分布[3,7]對強(qiáng)度 SAR圖像建模，結(jié)合方向邊緣強(qiáng)度信息，建立一種新的邊緣懲罰SAR圖像分割模型，提出一種最小化該模型的層次區(qū)域合并算法。該算法具有如下特點(diǎn)：(1)與文獻(xiàn)[11]不同，本文利用多方向比例邊緣檢測(Multi-Direction Ratio Edge Detector, MDRED)算子提取圖像的比例邊緣強(qiáng)度映射(Ratio Edge Strength Map, RESM)，該方法有兩個優(yōu)點(diǎn)：第一，MDRED算子的 CFAR特性避免了文獻(xiàn)[11]中的區(qū)域標(biāo)記運(yùn)算，降低了分割算法的時間復(fù)雜度；第二，降低了斑點(diǎn)噪聲對邊緣信息的影響。(2)分割模型中懲罰項(xiàng)中的邊緣方向信息增強(qiáng)了邊緣信息的抗噪能力，提高分割算法的性能。(3)分割模型中區(qū)域個數(shù)懲罰項(xiàng)，減少了噪聲引起的細(xì)碎區(qū)域。(4)利用區(qū)域鄰接圖(Region Adjacency Graph, RAG)表示分割結(jié)果和最近鄰圖(Nearest Neighbor Graph, NNG)加速模型求解。

2 邊緣懲罰SAR圖像分割模型

假設(shè)I:Ω →?,Ω={(x, y ) |x ∈{1,2,…,Nx},y∈ { 1,2,… , Ny}}為一幅強(qiáng)度SAR圖像，其中， Nx和 Ny分別表示圖像的行和列。圖像分割的目的是將圖像分割成 M 個互不相交的區(qū)域={ R1, R2,…,RM}。設(shè)區(qū)域Ri內(nèi)的像素值I（ x, y）服從視數(shù)為L和均值為μi的伽馬分布[7]。

式(1)假設(shè) SAR圖像中像素值滿足完全發(fā)展相干斑模型[12]，它適合描述大多數(shù)農(nóng)田場景的強(qiáng)度格式SAR圖像數(shù)據(jù)[3,7]。

假設(shè)圖像中任意兩個子區(qū)域 Ri與 Rj，像素值I( x1,y1)與 I ( x2,y2)之間相互獨(dú)立，SAR圖像的最終分割結(jié)果 ?*通過最大化似然函數(shù)（I）得到

通過對P?（I）取負(fù)對數(shù)，c）=- l n P?（I），式(2)等價于： ?*= a rg（?），

式中， Ni為區(qū)域 Ri的像素個數(shù)，μ=(1/ Ni)?I（ x, y ）。式(3)丟掉 - ln P?（I）中的常數(shù)項(xiàng)。

當(dāng)原始SAR圖像被分割成N,N = Nx× Ny個區(qū)域時，式(3)取最小值，然而，這樣的分割結(jié)果無意義。因此，需要在式(3)中加入對分割結(jié)果的懲罰項(xiàng)[3,7,11,14]。與文獻(xiàn)[10]中不同，本文利用方向邊緣強(qiáng)度信息，構(gòu)建一種新的邊緣懲罰項(xiàng)；同時，引入?yún)^(qū)域個數(shù)懲罰項(xiàng)，以防止分割結(jié)果中出現(xiàn)細(xì)碎的區(qū)域，得到本文提出的SAR圖像分割模型：

其中，iR?表示區(qū)域iR的邊緣像素集合，OESM(,,x y(,))x yθ是利用MDRED提取像素點(diǎn)(x,y)處（）,x yθ方向上的邊緣強(qiáng)度值，T是邊緣懲罰強(qiáng)度參數(shù)。式(4)等號右邊第1項(xiàng)為統(tǒng)計(jì)擬合度；第2項(xiàng)是區(qū)域邊緣懲罰項(xiàng)，其懲罰強(qiáng)弱與邊緣強(qiáng)度呈反比關(guān)系。權(quán)值λ起到均衡作用；第 3項(xiàng)是區(qū)域個數(shù)的懲罰項(xiàng)，它防止分割結(jié)果中出現(xiàn)細(xì)碎的小區(qū)域，η為權(quán)重參數(shù)。

3 邊緣懲罰層次區(qū)域合并模型求解

本文的模型求解方法由初始分割和邊緣懲罰強(qiáng)度遞增層次區(qū)域合并構(gòu)成。初始分割將原始圖像分割成均質(zhì)的小區(qū)域，利用邊緣懲罰強(qiáng)度遞增的層次區(qū)域合并技術(shù)遞歸地合并這些小區(qū)域，實(shí)現(xiàn)模型式(4)最小化。同時，利用RAG表示分割結(jié)果，并用NNG加速區(qū)域合并過程。

3.1 初始分割

利用 MDRED[13]提取 SAR圖像的 RESM，對閾值處理后的RESM(Thresholded RESM, TRESM)進(jìn)行分水嶺變換獲得初始過分割結(jié)果?。

MDRED是一種如圖1所示旋轉(zhuǎn)的雙邊平行矩形濾波器，其參數(shù)配置 Kf={l, w, d, θf}，分別為檢測器的長度、寬度、矩形之間的寬度和檢測器的方向。對于某一特定方向θf，先分別計(jì)算中心像素點(diǎn)(x,y)兩 邊 矩 形 區(qū) 域 內(nèi) 像 素 均 值m︿1( x, y, θf)和m︿2(x, y, θf)，然后計(jì)算(x,y)點(diǎn)沿 θf方向上的邊緣強(qiáng)度映射OESM(x, y, θf)：

圖1 邊緣檢測濾波器配置

OESM(x, y, θf) 具有如下特征：(1)在同質(zhì)區(qū)域內(nèi)部，任意方向上的OESM均很??；(2)在邊緣處，當(dāng)θf與邊緣方向一致時，OESM(x, y, θf)最大，當(dāng)θf與邊緣方向垂直時，OESM(x, y, θf) 最小。像素點(diǎn)(x,y)處的RESM為： E SM （ x, y ） =maxf{OESM(x, y, θf)}。

對RESM進(jìn)行閾值化處理，得到閾值處理以后的RESM：

其中，閾值 β 為 E SM（x, y ）的統(tǒng)計(jì)直方圖hist（ i）在α% 的左分位點(diǎn)處的值。在本文所有實(shí)驗(yàn)中，方向數(shù)K和分位數(shù)α分別取為16和0.35，邊緣檢測器的參數(shù)l, w和d分別為：9, 3和1。

圖2中給出了強(qiáng)度格式SAR圖像的初始過分割結(jié)果?？梢钥闯觯?1)圖像中絕大部分邊緣被提取出來；(2)每個區(qū)域的邊界單像素寬度，區(qū)域 Ri的邊界像素集為?Ri。

3.2 方向邊緣強(qiáng)度懲罰遞增區(qū)域合并技術(shù)

在區(qū)域合并技術(shù)中，RAG是表示分割結(jié)果的一種有效方法[15]。圖像分割的 RAG 被定義為一個無向圖G,G =(V, E, W )，其中，V為區(qū)域集，即v∈V,v ={(x, y ) |(x, y) ∈ Rv}; E為公共邊界集。若區(qū)域 Ru和Rv相鄰，則euv∈E,euv={(x, y ) |(x, y ) ∈?Ru∩ ? Rv}；相鄰區(qū)域 Ru和 Rv之間的權(quán)值 w ( u, v),w（ u, v）∈W，為合并相鄰區(qū)域Ru和Rv，式(4)的減少量為

圖2 初始分割結(jié)果

其中，|? Ru∩ ?Rv|表示公共邊界的長度，Ni和分別為區(qū)域 Ri( i = u , v)的均值和像素個數(shù)，為公共邊界像素的均值，和Nk分別為合并后區(qū)域Rk,Rk=Ru∪Rv∪（?Ru∩ ?Rv）的均值和像素個數(shù)。當(dāng)w（ u, v）＜ 0 時，合并相鄰區(qū)域 Ru和 Rv。

式(7)中邊界長度懲罰項(xiàng)的強(qiáng)弱由邊緣強(qiáng)度和參數(shù)T控制：邊緣懲罰強(qiáng)度與參數(shù)T值成正比關(guān)系，同時，與邊緣強(qiáng)度成反比關(guān)系。隨著參數(shù)T的增大，邊緣懲罰強(qiáng)度也增強(qiáng)。

式(7)等號右邊的3項(xiàng)之間是彼此競爭關(guān)系，第1項(xiàng)可以改寫成

其中， α =Nu/ Nk, α = Nv/ Nk，由算術(shù)與幾何平均之間的關(guān)系知，式(8)大于或者等于零，其作用是阻止區(qū)域合并，阻止力度隨著區(qū)域合并近似平方速度增大。它大于式(7)右邊促使區(qū)域合并的懲罰項(xiàng)的懲罰強(qiáng)度線性增強(qiáng)速度。隨著T增大，SAR圖像中較長弱邊緣得到保留。

利用如下兩個步驟計(jì)算式(7)中公共邊界長度懲罰項(xiàng)所需要的邊緣像素方向。

步驟 1 用分段線段近似邊緣。令 BAB表示點(diǎn)A= ( xA, yA)和B = ( xB, yB)之間的邊緣像素集，LAB表示連接點(diǎn)A和B之間的線段。表1給出了提取 BAB的分段線段 SAB的迭代二分法。

步驟2 對 SAB中的每條線段 LAB，計(jì)算方向：θAB= a rctan((yB- yA) /(xB- xA))。為簡化計(jì)算，邊緣 BAB上所有像素點(diǎn)的方向用初始分割中MDRED的均勻方向采樣θk來近似表示θAB，即

圖3顯示了利用表1中算法提取區(qū)域 Ri和 Rj之間的公共邊界的分段線段的示意圖。首先，用線段LAB近似公共輪廓 BAB(圖3(a))；然后，利用表1中算法遞歸地對線段 LAB進(jìn)行二分(圖3(b))，直到算法結(jié)束(圖3(d))。

表1 提取邊緣的線段連接近似表示的迭代二分法

圖3 遞歸提取公共邊界的線段(虛線表示邊緣，實(shí)線表示線段)

3.3 快速區(qū)域合并算法

假設(shè)由M個子區(qū)域構(gòu)成的分割結(jié)果的RAG為GM,GM=(VM, EM, WM)，當(dāng)WM中最小權(quán)值 wmin(u,v)滿足條件： wmin（u, v ）＜ 0 時，合并相鄰區(qū)域 Ru和Rv為新區(qū)域 Rk，更新RAG GM為 GM-1=(VM-1,EM-1, WM-1)，更新后的VM-1為

更新WM-1為

其中，利用式(7)重新計(jì)算與區(qū)域 Rk相鄰的區(qū)域之間的權(quán)值w（ k, j）。

為了加速區(qū)域合并過程，利用RAG的NNG加速最小權(quán)值 wmin（u, v ）的搜索。一個RAG G的NNG定義為有向圖 Gd,Gd=(Vd, Ed, Wd)，其中， Ed是有向 邊 集 ，如 果 w ( u, v) = min{w( u, k ) |k ∈N(u)},N (u )表示與節(jié)點(diǎn) u相鄰的節(jié)點(diǎn)集，存在有向邊＜u, v ＞∈Ed，而且，w（ k, j）∈Wd。圖4給出了RAG及其NNG的示意圖。一個RAG的NNG由V/ 2個子圖構(gòu)成，每個子圖僅有一個該子圖的最小權(quán)值環(huán)。因此，NNG中的最小權(quán)值環(huán)即為RAG中最小權(quán)值。

圖4 RAG及其NNG示意圖

表2給出本文提出的基于邊緣懲罰遞增的快速區(qū)域合并算法。

表2 遞增邊緣懲罰區(qū)域合并

3.4 算法特性

在區(qū)域合并算法中，需要預(yù)先設(shè)定權(quán)值λ和η與參數(shù)T，分析區(qū)域合并準(zhǔn)則式(7)的競爭關(guān)系知，權(quán)值λ和η與參數(shù)T之間存在如下關(guān)系：

(1)權(quán)值λ的取值較小時，參數(shù)T的步長 Ts可以取得較大，保證分割質(zhì)量的同時降低算法的時間復(fù)雜度(見式(13))；當(dāng)權(quán)值λ較大時，參數(shù)T的步長 Ts應(yīng)取得較小。

(2)權(quán)值λ和η的取值與 SAR 圖像的場景復(fù)雜程度有關(guān)。場景越復(fù)雜，其取值越小，以減少欠分割；對于簡單的場景，其取值較大，以降低過分割。因此，對于不同的SAR圖像，最優(yōu)權(quán)值λ和η與參數(shù)T的取值不同，折中考慮，本文實(shí)驗(yàn)中，λ=1.5,η= 4 , Ts= 0 .05。

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

為了驗(yàn)證本文方法的有效性，對圖5中4幅不同視數(shù)不同場景SAR圖像，將本文方法的分割結(jié)果與兩種基于模型優(yōu)化的分割方法(MDL方法[7]和IRGS方法[11])和一種基于特征的分割方法(Contextbased Hierarchical Unequal Merging, CHUM方法[5])進(jìn)行比較(圖6和圖7)。利用文獻(xiàn)[14]中的數(shù)值指標(biāo)來度量分割算法的性能，分析了本文算法的時間復(fù)雜度。

4.1 實(shí)驗(yàn)結(jié)果

對單視農(nóng)田場景(圖5(a))，MDL方法中存在大量的過分割和明顯欠分割現(xiàn)象。CHUM和IRGS方法中也同樣存在明顯的欠分割現(xiàn)象。同時，由于這兩種方法中梯度算子對斑點(diǎn)噪聲的敏感，使得區(qū)域輪廓出現(xiàn)不規(guī)則的鋸齒狀。對3視農(nóng)田場景(圖5(b)),MDL和 CHUM 方法中存在大量的漏檢邊緣。CHUM和IRGS方法中區(qū)域輪廓不光滑。而本文方法獲得光滑的區(qū)域輪廓，減少了欠分割和過分割。對于建筑物場景(圖 5(c))和城市區(qū)域(圖 5(d))的分割結(jié)果中。MDL方法中存在不同程度的過分割。雖然CHUM方法降低了過分割，但出現(xiàn)了不少漏檢輪廓，且輪廓不光滑。IRGS方法中，存在明顯的過分割現(xiàn)象，區(qū)域輪廓不光滑。本文方法中，不存在明顯的過分割和欠分割現(xiàn)象，提高了建筑物區(qū)域和城市區(qū)域的檢測能力。

4.2 性能分析

利用比值圖像的方差V和歸一化對數(shù)度量D[14]來度量不同SAR圖像分割算法的性能。方差V越小，算法性能越好；度量D的絕對值D越小，算法性能越好。

圖5 原始SAR圖像

圖6 分割結(jié)果1

圖7 分割結(jié)果2

表3給出了4種分割算法的數(shù)值指標(biāo)。本文方法的數(shù)值指標(biāo)大多優(yōu)于其他 3種方法。城市區(qū)域SAR圖像圖5(d)的分割結(jié)果中，MDL和IRGS方法的過分割使得它們的數(shù)值指標(biāo)均優(yōu)于本文方法。結(jié)合視覺和數(shù)值指標(biāo)比較，本文方法取得較好的分割結(jié)果。

表3 4種分割方法的性能比較

本文算法的時間復(fù)雜度主要由初始分割和區(qū)域合并的時間復(fù)雜度決定。初始分割的時間與SAR圖像的大小和場景復(fù)雜程度有關(guān)；表2中的區(qū)域合并算法的時間復(fù)雜度為

其中，NT為外循環(huán)次數(shù)，由步長因子 Ts決定，Ts越大，循環(huán)次數(shù)越少，時間復(fù)雜度越低；反之亦然；tini為算法初始化當(dāng)前RAG和NNG的時間，tup_RAG和tup_NNG分別為每次合并局部更新RAG和NNG的時間；M 為初始分割區(qū)域個數(shù)； tcmp為比較兩個權(quán)值的時間。從式(13)知，本文的區(qū)域合并算法的時間復(fù)雜度由外循環(huán)次數(shù)和初始分割區(qū)域個數(shù)決定。

MDL, CHUM和IRGS方法的時間復(fù)雜度分別由遞歸地變形多邊形網(wǎng)格、兩階段區(qū)域合并與區(qū)域分組和區(qū)域標(biāo)記與合并運(yùn)算決定。表4給出了4種分割算法的運(yùn)行時間，計(jì)算機(jī)平臺為：Pentium (R)Dual-Core, 2.93 GHz CPU, 2 GB 內(nèi)存，MATLAB 2010b。

5 結(jié)束語

本文提出一種新的SAR圖像分割算法。利用方向邊緣信息，構(gòu)造一種新的邊緣懲罰SAR圖像分割模型；利用MDRED和分水嶺變換獲得SAR圖像初始過分割結(jié)果；利用多邊形近似區(qū)域的輪廓提取邊緣的方向，結(jié)合MDRED提取圖像的方向邊緣信息；提出一種層次區(qū)域合并算法求解模型。利用RAG及其NNG表示圖像分割，提高區(qū)域合并的速度。實(shí)驗(yàn)結(jié)果表明：本文方法與其它3種方法相比在性能和效率方面都有優(yōu)勢，獲得更好的分割結(jié)果。

表4 4種分割方法的運(yùn)行時間(s)

[1] 李春升, 楊威, 王鵬波. 星載SAR成像處理算法綜述[J]. 雷達(dá)學(xué)報, 2013, 2(1): 111-122.Li Chun-sheng, Yang Wei, and Wang Peng-bo. A review of spaceborne SAR algorithm for image formation[J]. Journal of Radars, 2013, 2(1): 111-122.

[2] 鄧云凱, 趙鳳軍, 王宇. 星載SAR技術(shù)的發(fā)展趨勢及應(yīng)用淺析[J]. 雷達(dá)學(xué)報, 2012, 1(1): 1-10.Deng Yun-kai, Zhao Feng-jun, and Wang Yu. Brief analysis on the development and application of spaceborne SAR[J].Journal of Radars, 2012, 1(1): 1-10.

[3] Ayed I B, Mitiche A, and Belhadj Z. Multiregion level-set partitioning of synthetic aperture radar images[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence,2005, 27(5): 793-800.

[4] Gan Lu, Wu Yan, Wang Fan, et al.. Unsupervised SAR image segmentation based on triplet Markov fields with graph cuts[J]. IEEE Geoscience and Remote Sensing Letters, 2014,11(4): 853-857.

[5] Carvalho E A, Ushizima D M, Medeiros F N S, et al.. SAR imagery segmentation by statistical region growing and hierarchical merging[J]. Digital Signal Processing, 2009, 20(5):1365-1378.

[6] Yu Hang, Zhang Xiang-rong, Wang Shuang, et al.. Contextbased hierarchical unequal merging for SAR image segmentation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(2): 995-1009.

[7] Galland F, Bertaux N, and Refregier P. Minimum description length synthetic aperture radar image segmen-tation[J].IEEE Transactions on Image Processing, 2003, 12(9):995-1006.

[8] Kwon T J, Li J, and Wong A. ETVOS: an enhanced total variation optimization segmentation approach for SAR sea-ice image segmentation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(2): 925-934.

[9] Marques R C P, Medeiros F N, and Nobre J S. SAR image segmentation based on level set approach and G-zero-A model[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2012, 34(10): 2046-2057.

[10] 賈亞飛, 趙鳳軍, 禹衛(wèi)東, 等. 基于擴(kuò)散方程和MRF的SAR圖像分割[J]. 電子與信息學(xué)報, 2011, 33(2): 363-368.Jia Ya-fei, Zhao Feng-jun, Yu Wei-dong, et al.. SAR image segmentation based on diffusion equations and MRF[J].Journal of Electronics & Information Technology, 2011, 33(2):363-368.

[11] Yu Qiyao and Clausi D A. IRGS: image segmentation using edge penalties and region growing[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008, 30(12):2126-2139.

[12] Oliver C J. Information from SAR images[J]. Journal of Physics D: Applied Physics, 1991, 24(9): 1493-1514.

[13] Schou J, Skriver H, Nielsen A A, et al.. CFAR edge detector for polarimetric SAR images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(1): 20-32.

[14] Zhang Peng, Wu Yan, Li Ming, et al.. Unsupervised multiclass segmentation of SAR images using fuzzy triplet Markov fields model[J]. Pattern Recognition, 2012, 45(11): 4018-4033.[15] Peng B, Zhang L, and Zhang D. Automatic image segmentation by dynamic region merging[J]. IEEE Transactions on Image Processing, 2011, 20(12): 3592-3605.