-
摘要: 针对星载合成孔径雷达 (Synthetic aperture radar, SAR) 图像信噪比低、建筑物目标几何变形大以及周围背景复杂的特点, 本文提出了一种基于能量最小化的星载SAR图像建筑物分割方法.基于星载SAR图像数据构造条件概率能量项, 推动变形曲线向建筑物目标边界演化; 在能量泛函模型中定义长度能量项以保证变形曲线的平滑; 在水平集方法获取的SAR图像初始分割结果的基础上, 以高分辨率光学遥感影像中建筑物目标的轮廓作为先验信息, 构造先验形状能量项约束曲线在第二阶段的演化, 最终实现SAR图像建筑物的分割.实验结果表明, 该方法显著提高了建筑物目标轮廓的分割精度.Abstract: To aim at the difficulties of low signal to noise ratio (SNR) of spaceborne synthetic aperture radar (SAR) image, high geometric deformation and the complex background of buildings, a variational method is proposed for segmenting the building of spaceborne SAR image based on energy minimization. The condition probability energy term is defined with the data of SAR image to drive the deformation curve to the boundary of the building. The length energy term is constructed with the length of the evolving contour to ensure the smoothness of the deformation curve. The prior shape of the building is acquired from the corresponding high-resolution optical remote sensing image and the prior shape energy term is introduced into the energy functional. On the basis of a preliminary segmentation result which is achieved by the level set method, the building segmentation is implemented eventually by imposing the constraint of the prior shape. Experimental results demonstrate that the proposed model can improve the segmentation accuracy of the building target significantly.
-
Key words:
- Spaceborne SAR images /
- building segmentation /
- energy minimization /
- prior shape
-
星载合成孔径雷达 (Synthetic aperture radar, SAR) 具有运行周期固定、地面覆盖范围广、分辨率高等特点, 因此被应用于地球表面各种资源的探测.星载SAR图像被广泛应用于制图、土地监测、灾害监测、地质、水文、海洋、冰川等许多方面[1]. SAR图像中建筑物轮廓的抽取是一个难点, 其分割精度直接影响着后续解译的质量.
目前, SAR图像的分割方法主要有聚类方法[2-4]、阈值方法[5-7]、基于马尔科夫随机场的方法[8-12]和水平集方法[13-23]等.马秀丽等[2]结合分水岭和谱聚类算法实现了对SAR图像的分割.邓晓政等[3]提出了一种基于非负矩阵分解的谱聚类SAR图像分割方法.徐海霞等[4]针对谱聚类方法应用SAR图像分割时Laplace矩阵的特征值和特征向量难以计算的问题, 结合SAR图像的多尺度统计信息, 提出了基于谱聚类和混合模型的SAR图像分割方法.薛景浩等[5]在Rayleigh分布假设下, 提出一种最小误差的阈值化分割算法.聚类方法和阈值方法难以有效地利用SAR图像的边缘信息和空间信息, 当目标外观有歧义时, 这类算法的鲁棒性有待提高.余航等[8]针对SAR图像的特点, 提出了一种基于上下文分析的无监督分层迭代算法, 该算法结合了聚类算法和区域增长算法的优点, 提高了SAR图像分割的准确率. Deng等[9]采用MRF方法实现了SAR图像的非监督分割.傅兴玉等[11]在传统的MRF方法基础上, 提出了一种基于像素Gabor纹理相似度的邻域势函数模型, 改进了高分辨率SAR图像建筑物的分割效果.然而, MRF方法计算量巨大, 得到的分割结果受到相干斑噪声影响严重.
与上述分割方法相比, 水平集方法更能够适应图像的拓扑关系变化, 具有较强的实用性. Chan-Vese (CV) 模型[13]以分片光滑的区域描述图像的目标和背景, 对光学图像取得了较好的分割结果.为有效分割复杂背景下的感兴趣区域, Bresson等[14]基于待分割目标轮廓的先验知识, 在模型中增加了先验形状约束. Chen等[15]通过深度Boltzman机提取先验形状的分层结构, 在变分框架下通过耦合优化Boltzman机参数和变形曲线实现感兴趣目标的分割.田昊等[16]将多个先验形状竞争模型引入水平集方法中, 在标记函数的指导下, 利用先验形状能量来约束曲线的演化, 在对图像进行分割的同时完成建筑物的检测和提取.然而, 由于SAR发射的是纯相干波, 相干波经过与地物的相干作用, 特别是地物的后向散射作用, 目标的回波信号产生了衰减, 在图像中表现为相干斑噪声, 为了提高分割精度, 需要在模型中考虑噪声对分割结果的影响.涂松等[17]综合分析了基于活动轮廓模型的SAR图像分割方法. Mao等[18]使用复小波变换获取SAR图像的边缘信息, 然后使用CV模型抽取图像中的水边线. Zhang等[19]基于支持向量机获取SAR图像目标的初始轮廓, 再使用几何活动轮廓模型 (Geometric active contour, GAC) 实现SAR图像的分割. Ben等[20]对SAR图像基于Gamma分布建模, 通过水平集方法分割出SAR图像中的目标. Shuai等[21]在此基础上提出了一种静态全局最小能量函数分割模型, 改善了分割结果. Sui等[22]针对CV模型的缺点, 将CV模型中的灰度均值项改为Gamma概率分布项, 有效地降低了相干斑噪声对分割结果的影响.
上述文献提出的分割方法基于SAR图像的低层特征, 当目标周围背景复杂或在目标被部分遮挡的情况下会出现一定的分割误差.苏娟等[24]提出了一种基于SAR和可见光图像融合的SAR图像建筑物分割算法, 该方法充分利用可见光图像解译性好的特点, 提高了分割精度.本文结合高分辨率光学遥感影像, 提出了一种基于能量最小化的星载SAR图像建筑物分割模型:基于Gamma分布拟合SAR图像数据以构造条件概率能量项推动演化曲线向目标边缘演化, 定义演化曲线长度能量项以保证曲线在变形过程中的平滑; 根据先验知识构造先验形状能量项, 演化曲线在该能量项的推动下不断变形, 有效克服了目标被遮挡或复杂背景对分割结果造成的影响.实验结果表明, 该方法能够有效地提高复杂背景下SAR图像中目标轮廓的分割精度.
1. 星载SAR图像分割方法
1.1 分割模型算法框架
本文提出的分割算法流程如图 1所示.首先, 基于Gamma分布拟合SAR图像数据, 定义几何活动轮廓模型的曲线演化速度, 基于水平集方法对SAR图像分割得到粗分割结果; 再利用高分辨率光学遥感图像提取建筑物目标的先验形状, 以粗分割结果作为初始的轮廓线, 在先验形状约束下优化能量泛函直至收敛.在每一步迭代过程中, 需要将当前的演化曲线与先验形状进行配准.
1.2 基于能量最小化的SAR图像分割模型
由于SAR图像相干斑噪声的存在, 传统的混合Gaussian模型难以描述SAR图像的灰度分布.本文采用Gamma分布拟合图像数据. Gamma分布的概率密度函数 (Probability density function, PDF) 为[25]
$ \begin{align} \label{eq1} P_{\Gamma}^L(I\vert R) = \frac{L^L}{\mu _R(L - 1)!}\left( {\frac{I}{\mu _R }} \right)^{L - 1}{\rm e}^{ - \frac{LI}{\mu _R }} \end{align} $
(1) 其中, $I$ 为SAR图像的灰度, $R$ 为待分割SAR图像的目标和背景区域, $L$ 为SAR图像的等效视数, 本文取等效视数 $L = 1$ , $\mu _R $ 为区域的平均灰度值.
根据SAR图像数据的概率分布构造条件概率能量项
$ \begin{align} \label{eq2} E_{\rm con\_prob} = \sum\limits_{i = 1}^2 {\int_{R_i } { - \log P_{\Gamma}^L(I\vert R_i ){\rm d}\Omega } } \end{align} $
(2) 式中, $R_1 $ 和 $R_2 $ 分别为SAR图像的目标及背景区域, $\Omega \in {\bf R}^2$ 为图像的定义域.条件能量项的大小由目标和背景区域中SAR数据的概率分布给出, 该能量项推动曲线朝着目标的边界演化.
为使曲线在变形过程中保持平滑定义长度能量项
$ \begin{align} \label{eq3} E_{\rm length} = \int_\gamma {\rm d}s \end{align} $
(3) 式中, $\gamma $ 为当前演化曲线, 长度能量项的大小为演化曲线的长度.
由于SAR图像的信噪比低, 加之目标背景复杂, 仅依据条件概率能量项和长度能量项难以驱动演化曲线获取清晰、完整的轮廓.为了提高分割精度, 本文从高分辨率光学遥感图像中根据建筑物目标的先验形状构造先验能量项约束曲线的演化, 使得抽取的建筑物轮廓与先验形状相似.将先验形状以水平集方式表示, 设为 $\phi _s $ , 定义先验形状能量项
$ \begin{align} \label{eq4} E_{\rm prior} = \int_\Omega {(H(\phi(x, y)) - H(\tilde {\phi }_s(x_s, y_s )))^2} {\rm d}\Omega \end{align} $
(4) 式中, $H$ 为Heaviside函数, 当 $\phi \ge 0$ 时, $H(\phi ) = 1$ , 当 $\phi < 0$ 时, $H(\phi ) = 0$ . $\phi(x, y)$ 为当前演化曲线的水平集函数, $\tilde {\phi }_s(x_s, y_s )$ 为经仿射变换后的先验形状的水平集函数, $\tilde {\phi }_s(x_s, y_s )$ 定义为
$ \begin{align} \label{eq5} \tilde {\phi }_s(x_s, y_s ) = \phi _s(T^{ - 1}(x, y)) \end{align} $
(5) 式中, $T^{ - 1}$ 为演化曲线到先验形状之间的仿射变换, 考虑形状变化的几种可能情形:平移、缩放、旋转, 本文定义仿射变换参数 $P = \left[{t_x t_y s_x s_y \theta } \right]^{\rm T}$ , 其中, $t_x $ 为水平位移, $t_y $ 为垂直位移, $s_x $ 为水平尺度变换, $s_y $ 为垂直尺度变换, $\theta $ 为旋转角度. $(x, y)$ 和 $(x_s, y_s )$ 满足下面关系:
$ \begin{align} \left[{\begin{array}{l} x_s \\ y_s \\ \end{array}} \right] = &\ T^{ - 1}\left[{\begin{array}{l} x \\ y \\ \end{array}} \right] =\notag\\[2mm] & \left[{\begin{array}{l} {\kern 1pt} \dfrac{1}{s_x }[(x-x_g-t_x )\cos (\theta ) +\\ \qquad\quad(y-y_g - t_y )\sin(\theta )] + x_g \\[1mm] \dfrac{1}{s_y }[-(x-x_g-t_x )\sin(\theta ) +\\ \qquad\quad(y -y_g - t_y )\cos (\theta )] + y_g \\ \end{array}} \right] \end{align} $
(6) 式中, $(x_g, y_g )$ 为先验形状的中心点坐标.
为获取仿射变换参数, 使用梯度下降流极小化 $E_{\rm prior}$ 得到下式:
$ \begin{align} \label{eq7}\left\{ \begin{gathered} \frac{{{\rm d}{t_x}}}{{{\rm d}t}} = \iint_\Omega {(H({{\tilde \phi }_s}) - H(\phi )) \times 2\delta({{\tilde \phi }_s})} \times \\ \qquad\quad \left(\frac{{{{\tilde \phi }_{sx}}}}{{{s_x}}}\cos (\theta ) - \frac{{{{\tilde \phi }_{sy}}}}{{{s_y}}}{\rm sin}(\theta )\right){\rm d}x{\rm d}y \\[2mm] \frac{{{\rm d}{t_y}}}{{{\rm d}t}} = \iint_\Omega {(H({{\tilde \phi }_s}) - H(\phi )) \times 2\delta({{\tilde \phi }_s})} \times \\ \qquad\quad \left(\frac{{{{\tilde \phi }_{sx}}}}{{{s_x}}}\sin(\theta ) + \frac{{{{\tilde \phi }_{sy}}}}{{{s_y}}}\cos(\theta )\right){\rm d}x{\rm d}y \\[2mm] \frac{{{\rm d}{s_x}}}{{{\rm d}t}} = \iint_\Omega {(H({{\tilde \phi }_s}) - H(\phi ))} \times \\ \qquad\quad 2\delta({{\tilde \phi }_s}) \times \frac{{{{\tilde \phi }_{sx}}}}{{{s_x}}} \times {x_s}{\rm d}x{\rm d}y \\[2mm] \frac{{{\rm d}{s_y}}}{{{\rm d}t}} = \iint_\Omega {(H({{\tilde \phi }_s}) - H(\phi ))} \times \\ \qquad\quad 2\delta({{\tilde \phi }_s}) \times \frac{{{{\tilde \phi }_{sy}}}}{{{s_y}}} \times {y_s}{\rm d}x{\rm d}y \\[2mm] \frac{{{\rm d}\theta }}{{{\rm d}t}} = \iint_\Omega {(H({{\tilde \phi }_s}) - H(\phi )) \times 2\delta({{\tilde \phi }_s})} \times \\ \qquad\quad \left(\frac{{{{\tilde \phi }_{sy}}{s_x}}}{{{s_y}}} \times {x_s} - \frac{{{{\tilde \phi }_{sx}}{s_y}}}{{{s_x}}}{y_s}\right){\rm d}x{\rm d}y \end{gathered} \right. \end{align} $
(7) 式中, ${\tilde \phi _{sx}} =\dfrac {{\partial {{\tilde \phi }_s}}}{{\partial x}}$ , ${\tilde \phi _{sy}} =\dfrac{{\partial {{\tilde \phi }_s}}}{{\partial y}}$ .
演化曲线在条件概率能量项、长度能量项和先验形状能量项的作用下逐渐向目标轮廓演化, 模型的能量泛函为
$ \begin{align} \label{eq8} E = {\mu _1 E_{\rm con\_pro} + \mu _2 E_{\rm length} + \mu _3 E_{\rm prior} } \end{align} $
(8) 式中, $\mu _1 $ , $\mu _2 $ 和 $\mu _3 $ 分别为条件概率能量项、长度能量项和先验能量项权重系数. $\mu _1 $ 越大, 演化曲线在条件概率能量的驱动下向目标边界演化的速度越快; 增大 $\mu _2 $ , 可保证曲线演化平滑; 当SAR图像的信噪比较低或建筑物背景复杂时, 可增大 $\mu _3 $ , 演化曲线在先验能量的作用下将向与先验形状经仿射变换后相似的轮廓逼近.
计算式 (8) 的一阶变分, 基于梯度下降流得到曲线的演化方程
$ \begin{align} \frac{\partial \phi }{\partial t} =&\ \mu _1V\left\| {\nabla \phi } \right\| + \mu _2 \kappa \left\| {\nabla \phi } \right\| + \notag\\[1mm] &\ 2\mu _3(H(\phi ) - H(\tilde {\phi }_s ))\delta(\phi ) \end{align} $
(9) 式中, $V = -( {\log \hat {\mu }_{R_1 } + {I}/{\hat {\mu }_{R_1 } } - \log \hat {\mu }_{R_2 } - {I}/{\hat {\mu }_{R_2 } }} )$ , $\hat {\mu }_{R_1 } $ 和 $\hat {\mu }_{R_2 } $ 分别为当前演化曲线所划分的内外区域的灰度均值, $\delta $ 为Dirac函数, $\delta(\phi ) = {H}'(\phi )$ , $\kappa $ 为演化曲线的曲率.
为提高数值计算的效率, 避免迭代过程中重新初始化水平集函数, 本文使用有符号距离约束[26]
$ \begin{align} \label{eq10} E_{\rm Disreg} = \int_\Omega {\frac{1}{2}}(\left\| {\nabla \phi } \right\| - 1)^2{\rm d}\Omega \end{align} $
(10) 则式 (9) 变为
$ \begin{align} \frac{\partial \phi }{\partial t} =&\ \mu _1 V\left\| {\nabla \phi } \right\| + \mu _2 \kappa \left\| {\nabla \phi } \right\| +\notag\\[1mm] &\ 2\mu _3(H(\phi ) - H(\tilde {\phi }_s ))\delta(\phi ) +\\[1mm] &\ \mu _4 \left[\Delta \phi-{\rm div}\left(\frac{\nabla \phi }{\left\| {\nabla \phi } \right\|}\right)\right] \end{align} $
(11) 在SAR图像分割过程中, 为提高算法的执行效率, 首先令 $\mu _3 = 0$ , 利用条件概率能量和长度能量驱动演化曲线进行初步分割, 得到待分割目标的初步轮廓, 之后加入先验形状能量项进行分割, 从而得到最终的分割结果.
1.3 基于能量最小化的SAR建筑物分割算法步骤
基于能量最小化的SAR建筑物分割算法主要步骤如下:
步骤1.基于高分辨率光学图像获取感兴趣目标的轮廓, 作为目标的先验形状 $\phi _s $ .
步骤2.初始化:演化曲线的水平集函数为 $\phi ^0$ , 算法第一阶段最大迭代次数为 $Max\_Iter1$ , $n = 0$ , 设置模型各项的权重系数 $\mu _1$ , $\mu _2 $ , $\mu _3 $ 和 $\mu _4 $ .
步骤3.曲线进行第一阶段演化, 令式 (11) 中 $\mu _3= 0$ , 迭代:
while $n \le Max\_Iter1$ do
$ \begin{align*} &{\phi ^{n + 1}} = {\phi ^n} + \Delta t \times \\[1mm] &\qquad \Bigg\{ {\mu _1}V\left\| {\nabla {\phi ^n}} \right\| + {\mu _2}\kappa \left\| {\nabla {\phi ^n}} \right\| + \\[1mm] &\qquad {\mu _4}\left[\Delta {\phi ^n}-{\rm div}\left(\frac{{\nabla {\phi ^n}}}{{\left\| {\nabla {\phi ^n}} \right\|}}\right)\right]\Bigg\} \\[1mm] & n = n + 1 \end{align*} $
零水平曲线若停止演化, 转步骤4, 否则继续循环;
end
将输出结果作为粗分割结果 $\phi '$ .
步骤4.取 $\phi ^0=\phi '$ , 令算法第二阶段最大迭代次数为 $ Max\_Iter2$ , $n = 0$ , 仿射变换最大迭代次数为 $M$ , $m = 0$ ;
while $n \le Max\_Iter2$ do
初始化仿射变换的参数: $t_x = 0$ , $t_y = 0$ , $s_x$ $=$ $1$ , $s_y= 1$ , $\theta = 0$ ;
while $m \le M$ do
$ \begin{align*} &t_x^{m + 1} = t_x^m + \Delta t \times\\ &\qquad \Bigg\{ (H(\tilde \phi _s^m) - H({\phi ^n})) \times2\delta(\tilde \phi _s^m) \times\\ &\qquad \left(\frac{{\tilde \phi _{sx}^m}}{{s_x^m}}\cos({\theta ^m}) - \frac{{\tilde \phi _{sy}^m}}{{s_y^m}}\sin({\theta ^m})\right)\Bigg\}\\[1mm] &t_y^{m + 1} = t_y^m + \Delta t \times \\ &\qquad \Bigg\{(H(\tilde \phi _s^m) - H({\phi ^n})) \times 2\delta(\tilde \phi _s^m) \times \\ &\qquad \left(\frac{{\tilde \phi _{sx}^m}}{{s_x^m}}\sin({\theta ^m}) + \frac{{\tilde \phi _{sy}^m}}{{s_y^m}}\cos({\theta ^m})\right)\Bigg\}\\[1mm] &s_x^{m + 1} =s_x^m + \Delta t \times\\ &\qquad \Bigg\{ \left(H(\tilde \phi _s^m) - H({\phi ^n})\right) \times \\ &\qquad 2\delta(\tilde \phi _s^m) \times \frac{{\tilde \phi _{sx}^m}}{{s_x^m}} \cdot x_s^m\Bigg\}\\[1mm] & s_y^{m + 1} = s_y^m + \Delta t \times\\ &\qquad \Bigg\{ \left(H(\tilde \phi _s^m) - H({\phi ^n})\right) \times \\ &\qquad 2\delta(\tilde \phi _s^m) \times \frac{{\tilde \phi _{sy}^m}}{{s_y^m}} \times y_s^m\Bigg\}\\[1mm] &{\theta ^{m + 1}} = {\theta ^m} + \Delta t \times\\ &\qquad \Bigg\{ \left(H(\tilde \phi_s^m) - H({\phi ^n})\right) \cdot 2\delta(\tilde \phi _s^m) \times \\ &\qquad \left(\frac{{\tilde \phi _{sy}^ms_x^m}}{{s_y^m}} \times x_s^m - \frac{{\tilde \phi _{sx}^ms_y^m}}{{s_x^m}} \times y_s^m\right)\Bigg\}\\[1mm] & m = m + 1 \end{align*} $
end
$ \begin{align*} &{\phi ^{n + 1}} = {\phi ^n} + \Delta t \times \Bigg\{ {\mu _1}V\left\| {\nabla {\phi ^n}} \right\| + \\ &\qquad 2{\mu _3}\left(H({\phi ^n}) - H({\tilde \phi _s})\right)\delta({\phi ^n}) + {\mu _2}\kappa \left\| {\nabla {\phi ^n}} \right\| +\\ &\qquad {\mu _4}\left[\Delta {\phi ^n}-{\rm div}(\frac{{\nabla {\phi ^n}}}{{\left\| {\nabla {\phi ^n}} \right\|}})\right]\Bigg\} \end{align*} $
$n = n + 1$
零水平曲线若停止演化, 转步骤5, 否则继续循环;
end.
步骤5.输出最终分割结果 $\phi $
1.4 数值计算
离散化式 (11) 得到
$ \begin{align} \phi _{ij}^{n + 1} =&\ \phi _{ij}^n + \Delta t \times \bigg\{ {\mu _1}(\max(V_{ij}^n, 0){\nabla ^ + } + \\[1mm] &\ \min(V_{ij}^n, 0){\nabla ^ - }) + {\mu _2}\kappa _{ij}^n({(D_{ij}^{0x})^2} + \notag\\[1mm] & \ {(D_{ij}^{0y})^2})^{\frac{1}{2}} +2{\mu _3}\left( {H(\phi _{ij}^n) - H(\tilde \phi _{{s_{ij}}}^n)} \right)\times\notag\\[1mm] & \ \delta(\phi _{ij}^n) +{\mu _4}\left( {\Delta \phi _{ij}^n - \kappa _{ij}^n} \right)\bigg\} \end{align} $
(12) 式中, $\Delta t$ 为时间步长, $n$ 为迭代次数, $i$ , $j$ 为图像坐标, $\nabla ^ + $ 和 $\nabla ^ - $ 定义如下:
$ \begin{align} {\nabla ^ + } =&(\max {(D_{ij}^{ - x}, 0)^2} + \min {(D_{ij}^{ + x}, 0)^2} + \\[1mm] &\ \max {(D_{ij}^{ - y}, 0)^2} + \min {(D_{ij}^{ + y}, 0)^2}{)^{\frac{1}{2}}} \notag\\[2mm] {\nabla ^ - } =&(\max {(D_{ij}^{ + x}, 0)^2} + \min {(D_{ij}^{ - x}, 0)^2} +\\[1mm] & \ \max {(D_{ij}^{ + y}, 0)^2} + \min {(D_{ij}^{ - y}, 0)^2}{)^{\frac{1}{2}}} \end{align} $
(13) 其中, $D_{ij}^{ + x} $ , $D_{ij}^{ - x} $ , $D_{ij}^{0x} $ 和 $D_{ij}^{ + y} $ , $D_{ij}^{ - y} $ , $D_{ij}^{0y} $ 分别为 $x$ 和 $y$ 的前向差分、后向差分和中心差分格式.
曲率 $\kappa $ 采用中心差分格式离散化:
$ \begin{align} \label{eq14} \kappa _{ij}^k = \frac{D_{ij}^{xx}(D_{ij}^{0y} )^2 + 2D_{ij}^{0x} D_{ij}^{0y} D_{ij}^{xy} - D_{ij}^{yy}(D_{ij}^{0x} )^2}{\left( {(D_{ij}^{0x} )^2 +(D_{ij}^{0y} )^2} \right)^{\frac{3}{2}}} \end{align} $
(14) 水平集函数 $\phi $ 的拉普拉斯算子 $\Delta \phi _{ij}^n $ 的差分格式为
$ \begin{align} \label{eq15} \Delta \phi _{ij}^n = \phi _{i - 1, j}^n + \phi _{i + 1, j}^n + \phi _{ij - 1}^n + \phi _{ij + 1}^n - 4\phi _{ij}^n \end{align} $
(15) 在数值计算中, Heaviside函数和Dirac函数采用下面的形式进行计算:
$ \begin{align} &{H_\varepsilon }(z) =\begin{cases} 1, &z > \varepsilon \\ 0, & z < - \varepsilon \\ \dfrac{1}{2}\left(1 + \dfrac{2}{\pi }\arctan \dfrac{z}{\varepsilon }\right), &\left| z \right| \leq \varepsilon \end{cases} \notag\\[2mm] & {\delta _\varepsilon }(z) = \begin{cases} 0, &\left| { z} \right| > \varepsilon \\[1mm] \dfrac{1}{{2\varepsilon }}\left(1 + \cos \dfrac{{\pi z}}{\varepsilon }\right), & \left| { z} \right| \leq \varepsilon \end{cases} \end{align} $
(16) 2. 实验结果及分析
为了验证本文方法的有效性, 分别在带斑点噪声的合成图像和真实的SAR图像上进行分割实验, 并与CV方法[13]、Ben方法[20](算法实现的代码通过作者主页提供的链接下载, 参数设置通过实验进行了优化) 和Sui方法[22]所得到的分割结果进行比较.采用Dice (Dice similarity coefficient, DSC) 系数[27]作为指标对分割结果进行定量评价, DSC定义如下:
$ \begin{align} \label{eq17} DSC = \frac{2\times \vert A \cap B\vert }{\vert A\vert + \vert B\vert }\times 100 \% \end{align} $
(17) 式中, $A$ 为自动分割结果, $B$ 为手工分割结果. $DSC$ 越大, 分割效果越好, 当 $DSC$ 为100%时, $A$ 和 $B$ 完全匹配, 达到最佳分割效果.
2.1 带斑点噪声的合成图像分割
图 2(a)为使用Matlab的imnoise函数生成的带斑点噪声的合成图像, 斑点噪声的均值为1, 方差为0.04, 图像大小为100像素 $\times$ 100像素, 正方形目标右下角被遮挡; 图 2(b)为目标的先验形状.
图 3(a)中蓝色的初始曲线在条件概率能量和长度能量的共同作用下逐渐向内收缩, 红色曲线为迭代250次之后演化曲线的位置.当条件概率能量和长度能量达到最小时, 分割曲线停止了演化. 图 3(a)中的绿色曲线为模型没有加入先验形状约束的分割结果.为获取目标的完整轮廓, 在模型中引入由先验形状构造的先验形状能量项, 权重参数取 $\mu _1$ $=$ $1$ , $\mu _2 = 1$ , $\mu _3 = 1.5$ , $\mu _4 = 0.1$ , 时间步长 $\Delta t$ $=$ $0.5$ , 由于受到先验能量的作用, 图 3(e)中右下角处的曲线逐渐向外扩张并向先验形状逼近, 绿色曲线为迭代45次的结果, 迭代190次后, 模型最终获得了目标的较为完整的轮廓, 如图 3(e)中的红色曲线所示. 图 3(f)为手工分割的结果, 与本文模型相比, 图 3(b)的CV方法、图 3(c)的Ben方法和图 3(d)的Sui方法的演化曲线向目标的边缘演化, 可以有效地将目标和背景分割出来, 但是由于出现了遮挡, 这些模型均未能获取完整的目标轮廓.
为了研究图像噪声对本文方法的影响, 选取斑点噪声方差分别为0.1, 0.2, 0.3和0.4的图像进行分割实验, 结果如图 4所示.
从图 4可以看出, 本文方法在不同方差斑点噪声的影响下仍然能够分割出较为完整目标轮廓.但随着图像噪声方差的逐渐增大, 条件概率能量项受到较大的影响, 演化曲线难以完整分割出目标右下角的缺失部分, 得到的目标轮廓精度降低.
2.2 真实SAR图像建筑物分割
本文选取三幅真实SAR图像进行分割实验, 实验数据来自TerraSAR-X于2008年10月采集的上海交通大学区域的单视斜距图像 (Single look slant range complex, SSC) 及IKONOS 2007年3月采集的上海交通大学区域图像.如图 5(a)所示, SAR图像大小为106像素 $\times$ 145像素, 从对应的分辨率为1 m的IKONOS全色波段图像 (图 5(b)) 中抽取" H "状楼的轮廓作为先验形状 (图 5(c)).由于在分割过程中需要将变形曲线与先验形状配准, 因此不需要预先配准SAR和IKONOS图像.
图 5 " H "状楼的SAR图像和全色波段图像及其先验形状Fig. 5 SAR image of the building with " H " contour, its corresponding panchromatic image and the prior shape图 6为不同方法的分割结果.由于SAR图像存在大量强散射斑块, 在同质区域内表现为部分像素灰度值较高, 部分像素灰度值较低; 建筑物内部纹理复杂且轮廓不完整. CV方法和Ben方法只能分割出图像中的一些高亮线条, 无法得到完整的轮廓, 分割结果中存在较多孤立像素. Gamma分布能够较好地反映SAR图像数据的统计特性, Sui模型结合了Gamma分布和CV模型, 因此Sui模型的分割结果优于前两种方法.图 6(e)为手工标注的分割结果, 本文模型引入了建筑物的先验形状信息, 模型参数设置为 $\mu _1 = 1$ , $\mu _2 = 1$ , $\mu _3 = 2$ , $\mu _4 = 1$ , 时间步长 $\Delta t = 0.1$ , 最终分割结果由SAR图像建筑物的灰度分布特性及先验形状确定, 与Sui模型相比, 抽取的轮廓更为完整, 分割精度更高.
图 6 不同方法对" H "字楼的分割结果Fig. 6 Segmentation results of the building with " H " contour by different methods第二组实验SAR图像如图 7(a)所示, 图像大小为136像素 $\times$ 73像素, 从对应的IKONOS全色波段图像 (图 7(b)) 中抽取目标建筑物的轮廓作为先验形状 (图 7(c)).
图 7 " L "状楼的SAR图像和全色波段图像及其先验形状Fig. 7 SAR image of the building with " L " contour, its corresponding panchromatic image and the prior shape图 8为不同方法的分割结果.由于SAR图像相干斑噪声严重, 目标建筑物在SAR图像中无明显轮廓, CV方法和Ben方法只能分割出图像中的小块的分块区域. Sui方法优于前两种方法, 可以得到较为完整的形状, 但未能有效地将目标与背景分割开来, 且存在着明显的误分割现象.本文模型的参数设置为 $\mu _1 = 1$ , $\mu _2 = 1$ , $\mu _3 = 2$ , $\mu _4 = 0.5$ , 时间步长 $\Delta t$ $=$ $0.1$ , 在先验形状的约束下分割得到的目标建筑物轮廓更加精确完整.
图 8 不同方法对" L "状楼的分割结果Fig. 8 Segmentation results of the building with " L " contour by different methods为进一步验证本文方法在复杂背景下的分割效果, 选择了体育馆作为待分割的目标如图 9(a)所示.图像的大小为221像素 $\times$ 201像素; 从对应的IKONOS全色波段图像 (图 9(b)) 中抽取建筑物的轮廓作为先验形状 (图 9(c)).
图 9 体育馆的SAR图像和全色波段图像及其先验形状Fig. 9 The SAR image, its corresponding panchromatic image and the prior shape of the gymnasium图 10为不同方法的分割结果.由于SAR图像中建筑物右上方边缘断裂, CV方法和Ben方法得到的轮廓在该处出现凹陷现象.加之建筑物内部纹理结构复杂, Ben方法得到的轮廓内部存在较多的孔洞. Sui方法虽然优于前两种方法, 但由于该建筑物周围背景较为复杂, Sui方法未能有效地将目标与背景分割开来, 存在大量误分割现象.本文模型的参数设置为 $\mu _1 = 1$ , $\mu _2 = 1$ , $\mu _3 = 2$ , $\mu _4 = 1$ , 时间步长 $\Delta t = 0.1$ , 在先验形状的约束下分割得到的建筑物边界更为清晰平滑.
分割实验的精度指标DSC及运行时间 (运行环境: Intel Core I7, 4 GB RAM; Matlab 2012 a, 操作系统为Windows 7) 如表 1所示.从表中可以看出, 在SAR仿真图像和真实图像分割实验中, 本文方法的分割精度接近或超过90%, 优于其他方法.引入待分割目标的先验信息, 显著提高了分割精度.在算法的时间效率上, 本文方法优于CV方法, 但低于Ben和Sui方法.综合考虑分割精度和运行时间, 本文的方法具有更强的分割性能.
表 1 实验结果比较Table 1 Comparison of segmentation results实验名称 CV方法 Ben方法 Sui方法 本文方法 精度 (%) 耗时 (s) 精度 (%) 耗时 (s) 精度 (%) 耗时 (s) 精度 (%) 耗时 (s) 合成图像 97.4 15.4 97.4 21.8 97.3 16.9 99.4 113.7 " H "状楼 28.7 133.9 57.3 82.9 84.0 52.9 89.6 116.1 " L "状楼 13.4 162.3 21.7 20.8 62.9 13.7 91.5 67.1 体育馆 82.0 326.2 76.1 88.0 86.3 101.9 94.1 206.1 平均值 55.4 159.5 63.1 53.3 82.6 46.4 93.7 125.8 3. 结论
根据星载SAR图像的特点, 本文提出了一种基于能量最小化的星载SAR图像建筑物分割方法, 能够有效地提取建筑物目标的完整轮廓.通过星载SAR图像的灰度特征和GAC模型, 获取初步的分割结果, 然后从高分辨率光学图像提取感兴趣建筑物的先验形状以构造先验能量, 约束变形曲线直至收敛到目标的边界.模型中的参数通过实验设置, 如何自适应地选取权重系数以及根据SAR图像的特点采用更合适的分布模型以进一步提高建筑物目标的分割精度是下一步研究的方向.
-
图 5 " H "状楼的SAR图像和全色波段图像及其先验形状
Fig. 5 SAR image of the building with " H " contour, its corresponding panchromatic image and the prior shape
图 6 不同方法对" H "字楼的分割结果
Fig. 6 Segmentation results of the building with " H " contour by different methods
图 7 " L "状楼的SAR图像和全色波段图像及其先验形状
Fig. 7 SAR image of the building with " L " contour, its corresponding panchromatic image and the prior shape
图 8 不同方法对" L "状楼的分割结果
Fig. 8 Segmentation results of the building with " L " contour by different methods
图 9 体育馆的SAR图像和全色波段图像及其先验形状
Fig. 9 The SAR image, its corresponding panchromatic image and the prior shape of the gymnasium
表 1 实验结果比较
Table 1 Comparison of segmentation results
实验名称 CV方法 Ben方法 Sui方法 本文方法 精度 (%) 耗时 (s) 精度 (%) 耗时 (s) 精度 (%) 耗时 (s) 精度 (%) 耗时 (s) 合成图像 97.4 15.4 97.4 21.8 97.3 16.9 99.4 113.7 " H "状楼 28.7 133.9 57.3 82.9 84.0 52.9 89.6 116.1 " L "状楼 13.4 162.3 21.7 20.8 62.9 13.7 91.5 67.1 体育馆 82.0 326.2 76.1 88.0 86.3 101.9 94.1 206.1 平均值 55.4 159.5 63.1 53.3 82.6 46.4 93.7 125.8 
下载: 导出CSV
-
[1] Bovolo F, Marin C, Bruzzone L. A hierarchical approach to change detection in very high resolution SAR images for surveillance applications. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(4):2042-2054 doi: 10.1109/TGRS.2012.2223219 [2] 马秀丽, 焦李成.基于分水岭--谱聚类的SAR图像分割.红外与毫米波学报, 2008, 27(6):452-456 http://www.cnki.com.cn/Article/CJFDTOTAL-HWYH200806013.htmMa Xiu-Li, Jiao Li-Cheng. SAR image segmentation based on watershed and spectral clustering. Journal of Infrared and Millimeter Waves, 2008, 27(6):452-456 http://www.cnki.com.cn/Article/CJFDTOTAL-HWYH200806013.htm [3] 邓晓政, 焦李成, 卢山.基于非负矩阵分解的谱聚类集成SAR图像分割.电子学报, 2011, 39(12):2905-2909 http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201112030.htmDeng Xiao-Zheng, Jiao Li-Cheng, Lu Shan. Spectral clustering ensemble applied to SAR image segmentation using nonnegative matrix factorization. Acta Electronica Sinica, 2011, 39(12):2905-2909 http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201112030.htm [4] 徐海霞, 田铮, 丁明涛.基于谱聚类与混合模型的SAR图像多尺度分割.中国图象图形学报, 2010, 15(3):450-454 http://www.cnki.com.cn/Article/CJFDTOTAL-ZGTB201003018.htmXu Hai-Xia, Tian Zheng, Ding Ming-Tao. Multiscale segmentation for SAR image based on spectral clustering and mixture model. Journal of Image and Graphics, 2010, 15(3):450-454 http://www.cnki.com.cn/Article/CJFDTOTAL-ZGTB201003018.htm [5] 薛景浩, 章毓晋, 林行刚. SAR图像基于Rayleigh分布假设的最小误差阂值化分割.电子科学学刊, 1999, 21(2):219-225 http://www.cnki.com.cn/Article/CJFDTOTAL-DZYX199902012.htmXue Jing-Hao, Zhang Yu-Jin, Lin Xing-Gang. Rayleigh-distribution based minimum error thresholding for SAR images. Journal of Electronics, 1999, 21(2):219-225 http://www.cnki.com.cn/Article/CJFDTOTAL-DZYX199902012.htm [6] Han C M, Guo H D, Shao Y, Liao J J. A method to segment SAR images based on histogram. In:Proceedings of the 2005 IEEE International Geoscience and Remote Sensing Symposium. Seoul, Korea (South):IEEE, 2005. 3694-3696 [7] 龙建武, 申铉京, 陈海鹏.自适应最小误差阈值分割算法.自动化学报, 2012, 38(7):1134-1144 doi: 10.3724/SP.J.1004.2012.01134Long Jian-Wu, Shen Xuan-Jing, Chen Hai-Peng. Adaptive minimum error thresholding algorithm. Acta Automatica Sinica, 2012, 38(7):1134-1144 doi: 10.3724/SP.J.1004.2012.01134 [8] 余航, 焦李成, 刘芳.基于上下文分析的无监督分层迭代算法用于SAR图像分割.自动化学报, 2014, 40(1):100-116 http://www.aas.net.cn/CN/abstract/abstract18271.shtmlYu Hang, Jiao Li-Cheng, Liu Fang. Context based unsupervised hierarchical iterative algorithm for SAR segmentation. Acta Automatica Sinica, 2014, 40(1):100-116 http://www.aas.net.cn/CN/abstract/abstract18271.shtml [9] Deng H W, Clausi D A. Unsupervised segmentation of synthetic aperture radar sea ice imagery using a novel Markov random field model. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(3):528-538 doi: 10.1109/TGRS.2004.839589 [10] 侯彪, 徐婧, 刘枫, 焦李成.基于第二代Bandelet域隐马尔科夫树模型的图像分割.自动化学报, 2009, 35(5):498-504 http://www.aas.net.cn/CN/abstract/abstract13437.shtmlHou Biao, Xu Jing, Liu Feng, Jiao Li-Cheng. Image segmentation using second generation Bandelet-domain hidden Markov tree models. Acta Automatica Sinica, 2009, 35(5):498-504 http://www.aas.net.cn/CN/abstract/abstract13437.shtml [11] 傅兴玉, 尤红健, 付琨.基于改进Markov随机场的高分辨率SAR图像建筑物分割算法.电子学报, 2012, 40(6):1141-1147 http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201206013.htmFu Xing-Yu, You Hong-Jian, Fu Kun. Building segmentation from high-resolution SAR images based on improved Markov random field. Chinese Journal of Electronics, 2012, 40(6):1141-1147 http://www.cnki.com.cn/Article/CJFDTOTAL-DZXU201206013.htm [12] 宋艳涛, 纪则轩, 孙权森.基于图像片马尔科夫随机场的脑MR图像分割算法.自动化学报, 2014, 40(8):1754-1763 http://www.aas.net.cn/CN/abstract/abstract18442.shtmlSong Yan-Tao, Ji Ze-Xuan, Sun Quan-Sen. Brain MR image segmentation algorithm based on Markov random field with image patch. Acta Automatica Sinica, 2014, 40(8):1754-1763 http://www.aas.net.cn/CN/abstract/abstract18442.shtml [13] Chan T F, Vese L A. Active contours without edges. IEEE Transactions on Image Processing, 2001, 10(2):266-277 doi: 10.1109/83.902291 [14] Bresson X, Vandergheynst P, Thiran J P. A variational model for object segmentation using boundary information and shape prior driven by the Mumford-Shah functional. International Journal of Computer Vision, 2006, 68(2):145-162 doi: 10.1007/s11263-006-6658-x [15] Chen F, Yu H M, Hu R, Zeng X X. Deep learning shape priors for object segmentation. In:Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition. Portland, OR:IEEE, 2013. 1870-1877 [16] 田昊, 杨剑, 汪彦明, 李国辉.基于先验形状约束水平集模型的建筑物提取方法.自动化学报, 2010, 36(11):1502-1511 doi: 10.3724/SP.J.1004.2010.01502Tian Hao, Yang Jian, Wang Yan-Ming, Li Guo-Hui. Towards automatic building extraction:variational level set model using prior shape knowledge. Acta Automatica Sinica, 2010, 36(11):1502-1511 doi: 10.3724/SP.J.1004.2010.01502 [17] 涂松, 李禹, 粟毅.基于主动轮廓模型的SAR图像分割方法综述.系统工程与电子技术, 2015, 37(8):1754-1766 http://www.cnki.com.cn/Article/CJFDTOTAL-XTYD201508008.htmTu Song, Li Yu, Su Yi. Overview of SAR image segmentation based on active contour model. Systems Engineering and Electronics, 2015, 37(8):1754-1766 http://www.cnki.com.cn/Article/CJFDTOTAL-XTYD201508008.htm [18] Mao C L, Wan S H. A water/land segmentation algorithm based on an improved Chan-Vese model with edge constraints of complex wavelet domain. Chinese Journal of Electronics, 2015, 24(2):361-365 doi: 10.1049/cje.2015.04.023 [19] Zhang S M, He F, Zhang Y L, Wang J M, Mei X, Feng T T. Geometric active contour based approach for segmentation of high-resolution spaceborne SAR images. Journal of Systems Engineering and Electronics, 2015, 26(1):69-76 doi: 10.1109/JSEE.2015.00010 [20] Ben A I, Mitiche A, Belhadj Z. Multiregion level-set partitioning of synthetic aperture radar images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2005, 27(5):793-800 doi: 10.1109/TPAMI.2005.106 [21] Shuai Y M, Sun H, Xu G. SAR image segmentation based on level set with stationary global minimum. IEEE Geoscience and Remote Sensing Letters, 2008, 5(4):644-648 doi: 10.1109/LGRS.2008.2001768 [22] Sui H G, Xu C, Liu J Y, Sun K M, Wen C F. A novel multi-scale level set method for SAR image segmentation based on a statistical model. International Journal of Remote Sensing, 2012, 33(17):5600-5614 doi: 10.1080/01431161.2012.666814 [23] 唐利明, 田学全, 黄大荣, 王晓峰.结合FCMS与变分水平集的图像分割模型.自动化学报, 2014, 40(6):1233-1248 http://www.aas.net.cn/CN/abstract/abstract18394.shtmlTang Li-Ming, Tian Xue-Quan, Huang Da-Rong, Wang Xiao-Feng. Image segmentation model combined with FCMS and variational level set. Acta Automatica Sinica, 2014, 40(6):1233-1248 http://www.aas.net.cn/CN/abstract/abstract18394.shtml [24] 苏娟, 鲜勇, 宋建社. SAR图像与可见光图像融合的建筑物提取算法.兵工学报, 2010, 31(11):1448-1454 http://www.cnki.com.cn/Article/CJFDTOTAL-BIGO201011008.htmSu Juan, Xian Yong, Song Jian-She. A building extraction algorithm based on the fusion of SAR and optical images. Acta Armamentarii, 2010, 31(11):1448-1454 http://www.cnki.com.cn/Article/CJFDTOTAL-BIGO201011008.htm [25] Goodman J W. Statistical properties of laser speckle patterns. Laser Speckle and Related Phenomena. Berlin Heidelberg:Spring-Verlag, 1975. 9-75 [26] Li C M, Xu C Y, Gui C F, Fox M D. Level set evolution without re-initialization:a new variational formulation. In:Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. San Diego, CA, USA:IEEE, 2005. 430-436 [27] Zou K H, Warfield S K, Bharatha A, Tempany C M C, Kaus M R, Haker S J, Wells W M Ⅲ, Jolesz F A, Kikinis R. Statistical validation of image segmentation quality based on a spatial overlap index. Academic Radiology, 2004, 11(2):178-189 doi: 10.1016/S1076-6332(03)00671-8 -
下载:
计量
- 文章访问数: 1725
- HTML全文浏览量: 289
- PDF下载量: 866
- 被引次数: 0
