边缘惩罚层次区域合并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.