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

張澤均 水鵬朗

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

3 邊緣懲罰層次區域合并模型求解

本文的模型求解方法由初始分割和邊緣懲罰強度遞增層次區域合并構成。初始分割將原始圖像分割成均質的小區域，利用邊緣懲罰強度遞增的層次區域合并技術遞歸地合并這些小區域，實現模型式(4)最小化。同時，利用RAG表示分割結果，并用NNG加速區域合并過程。

3.1 初始分割

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

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

圖1 邊緣檢測濾波器配置

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

對RESM進行閾值化處理，得到閾值處理以后的RESM：

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

圖2中給出了強度格式SAR圖像的初始過分割結果。可以看出：(1)圖像中絕大部分邊緣被提取出來；(2)每個區域的邊界單像素寬度，區域 Ri的邊界像素集為?Ri。

3.2 方向邊緣強度懲罰遞增區域合并技術

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

圖2 初始分割結果

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

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

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

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

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

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

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

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

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

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

3.3 快速區域合并算法

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

更新WM-1為

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

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

圖4 RAG及其NNG示意圖

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

表2 遞增邊緣懲罰區域合并

3.4 算法特性

在區域合并算法中，需要預先設定權值λ和η與參數T，分析區域合并準則式(7)的競爭關系知，權值λ和η與參數T之間存在如下關系：

(1)權值λ的取值較小時，參數T的步長 Ts可以取得較大，保證分割質量的同時降低算法的時間復雜度(見式(13))；當權值λ較大時，參數T的步長 Ts應取得較小。

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

4 實驗結果與性能分析

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

4.1 實驗結果

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

4.2 性能分析

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

圖5 原始SAR圖像

圖6 分割結果1

圖7 分割結果2

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

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

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

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

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

5 結束語

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

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

[1] 李春升, 楊威, 王鵬波. 星載SAR成像處理算法綜述[J]. 雷達學報, 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技術的發展趨勢及應用淺析[J]. 雷達學報, 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] 賈亞飛, 趙鳳軍, 禹衛東, 等. 基于擴散方程和MRF的SAR圖像分割[J]. 電子與信息學報, 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.