1 引言

合成孔径雷达(SAR)是一种主动式微波传感器,由于其具有全天时、全天候、穿透力强以及能连续观测运动目标等特点,SAR已经在民用及**领域得到日益广泛的应用^[1-2]。水体是SAR图像解译的重要内容,如在2008年5·12汶川地震抗震救灾中,堰塞湖及水库对灾区构成了极大威胁,而利用SAR对其进行持续的监控及评估是减少次生灾害的重要手段。在SAR图像解译过程中,往往需要对图像进行变化检测、信息融合及图像镶嵌等处理,这就需要对SAR图像进行配准,配准精度将直接影响后续图像处理的效果。目前,图像配准方法主要可分为基于灰度的方法和基于特征的方法。由于结构特征对图像的灰度变换、变形及遮挡具有更好的稳健性,因而基于特征的方法是SAR图像配准的主要方法。而精确地提取SAR图像特征是实现高精度图像配准的基础。目前常见的用于SAR图像配准的特征有点特征^[3]、线特征^[4]和区域特征^[5]等。由于在SAR成像的过程中,水体表面在针对微波波长成像时具有镜面反射特征而呈现为低强度区域,而与此相对的那些非水体通常具有粗糙表面,它们在成像过程中由于漫散射而呈现为亮的区域^[6]。水体和非水体的成像特征使得它们在连接处的灰度值具有明显的差异,不易出现SAR图像中常见的边缘模糊问题,比较适合用于提取精确的目标轮廓,因而轮廓是包含水体的SAR图像较为理想的配准特征。目前,目标轮廓的提取方法主要包含两类,一类是直接利用边缘检测方法提取目标的边缘轮廓,如Canny边缘检测算子、Sobel方法和Log算子等,这些方法对光学图像边缘检测具有较好的效果,但在检测SAR图像的边缘时则容易出现边缘断裂及虚假边缘等问题;另一类方法是先对SAR图像进行分割,再由分割图像来提取目标轮廓,这种方法比较适合用于含水体SAR图像目标轮廓的提取。

由于乘性斑点噪声的影响及地物回波的反射特性,SAR图像表现出很强的不确定性和模糊性。结合了多特征信息和上下文信息的马尔科夫随机场(MRF)模型是SAR图像分割的有效模型^[7-8]。 根据贝叶斯准则, 图像分割可以表示为已知观测图像下标记场的最大后验概率或等价的Gibbs场的能量最小化问题。能量函数可分为两个部分:数据模型因子和区域标记因子,数据模型因子包含了图像分割过程中图像强度等数据特征的约束;区域标记因子则包含了图像分割过程中标记场的光滑性约束。模型中的加权参数反映了两种约束之间的相互作用^[9-10]。通常对于图像中的均质区域需要减少数据约束以利于得到光滑的分割区域,而对于非均质区域则需要增加数据约束,以利于保持目标的边缘与细节。为此,本文将观测图像的局部均质性特征信息融入MRF分割模型,使模型的加权参数能根据观测图像的均质性自适应地改变,以引导其得到和观测图像特征相一致的分割结果。

常用的基于轮廓的图像配准方法可分为两类,一类方法是直接利用轮廓上的点或角点进行匹配^[11-12];另一类方法是利用轮廓的描述子进行匹配,如形状上下文和不变量等^[13]。近年来,基于高斯混合模型(GMM)的相干点漂移算法得到了人们广泛的研究。该算法将两个点集的对应问题转化为GMM概率密度函数的参数估计问题,可以对两个具有刚体变换及非刚体变换的点集实现有效的匹配^[14-17]。文献[ 18]给出了同时融合特征点位置信息和尺度不变特征转换(SIFT)特征的非均匀GMM,该模型能显著提高点集匹配的精度和稳健性,但该方法难以应用于轮廓的匹配。为此,本文提出了另一种非均匀GMM,该模型除了考虑轮廓点的位置信息外,还包含了以轮廓点为中心的窗口内图像灰度的统计相似性特征。

2 轮廓特征提取

2.1 图像分割的MRF模型

设S={s=(i,j)|1≤i≤H,1≤j≤W,i,j,H,W∈Z⁺}表示一个二维点集,每一个点表示图像中的一个网格位置,H和W分别表示图像的长度和宽度。X=(Xs)s∈S和x=(xs)s∈S分别表示不可观测的随机场及其实现。Y=(Ys)s∈S和y=(ys)s∈S表示观测图像对应的随机场及其实现。x_s∈{1,2,…,L}表示位置s所对应的观测图像像素所属类的标记。N={N_s,s∈S}为邻域系统,其中N_s,s∈S表示位置s的邻居。C={c|c⊂N_s}为此邻域系统的基团集合。图像分割就是要确定观测图像中每个像素类的标记。

在贝叶斯框架下,分割结果 x^通常表示为p(X|Y)的最大后验概率(MAP)解,也即:

x^=argmaxx[p(X=x|Y=y)]=argmaxx[p(X=x)p(Y=y|X=x)]。(1)

为了计算的可行性,通常假设MRF分割模型的观测场是关于标记场独立的,且标记场中的每一个特征值只和与其对应的标记场的值有关,即:

p(Y=y|X=x)=∏s∈Sp(Ys=ys|Xs=xs)=∏s∈Spxs(ys),(2)

与此对应的局部最优解可以表示为

x^s=argmaxxs∈Lpys|xs)p(xs|xNs)]。(3)

对多层逻辑(MLL)模型,由贝叶斯准则和Hamersley-Clifford 定理可知图像的分割问题可以表示为最大后验概率或Gibbs场能量最小化问题,即:

p(X=x|Y=y)=1Zexp[-U(x)+∑s∈Slogf(ys|Xs=xs)]。(4)

令

ER(x)=∑sUs[xs,(xt)t∈Ns]=∑s∑t∈Nsβsδ(xs,xt),(5)EY=∑s∈S-logf(ys|Xs=xs)(6)

式中δ(x_s,x_t)=0,xs=xt1,xs≠xt。上述标记场的最大后验概率等价于最小化的能量函数为

E=ER+αEY,(7)

式中α为加权参数。为叙述方便,称E_R为区域标记因子,E_Y为数据模型因子。

2.2 融合观测图像局部信息的MRF分割模型

经典的MRF模型为了应用的方便,通常假设每一观测图像只与该位置的标记相关。由此得到的模型只可以对标记图像中的上下文信息建模,没有有效地利用观测图像的上下文信息,适当地利用观测图像的局部信息显然能提高模型的表现力。在(7)式表示的模型中,α表示区域标记因子和数据模型因子在分割模型中的相对重要性。α越小,区域标记因子在分割模型中的作用越大,得到的分割图像也越光滑,但可能导致欠分割问题;α越大,将会保留更多的细节,但可能会导致过分割问题。由于分割图像在初始阶段具有较高的复杂性,在收敛过程中,分割图像会越来越光滑。为了和这个过程相一致,模型需要在初始阶段使用一个较大的α值,此后,该值需随着迭代的过程逐渐变小。为此,文献[ 9]采用了可变参数,表示为

αs(k)=80×0.9k+1,(8)

式中k表示迭代数。

由于图像在不同的区域具有不同的特征,如在灰度比较均匀的均质区域,需要取较小的α值,以便得到光滑的分割结果;而在非均质区域,则需要取较大的α值,以保持边缘和细节。这就需要模型能根据观测图像的局部特征自适应地改变模型参数。由于图像中均质区域的灰度值取值更集中,变化较小,具有较小的样本方差。相反地,当图像局部区域内包含了细节或边缘时,区域内的灰度取值将更分散,具有较大的样本方差。因而对于模型中的点s的加权参数α_s,可以以窗口N_s(窗口中心为s,半径为n)内灰度的样本标准差作为均质性的衡量标准,其表达式为

VstdNs=1(2n+1)2-1∑i=-nn∑ j=-nn(xij-x-)2,(9)

式中x_ij表示在图像位置(i,j)的灰度值, x-表示窗口N_s内像素灰度的均值。综合考虑文献[ 9]提及的因素及(8)式,给出的α值为

αs(k)=c×VstdNs×0.95k+1,(10)

式中k表示迭代数,c为某个正的常数。α_s(k)值能根据观测图像的内容而自适应地改变,对于均质区域,α_s(k)取较小的值;对于非均质区域,则取较大的值。

3 非均匀分布的相干点漂移算法

3.1 相干点漂移算法

相干点漂移算法将两个点集的对应问题转化为概率密度的估计问题。其核心思想是设模板点集为一GMM中心,另一个点集为GMM生成的数据点集。在优化过程中,两个点集的对应通过最大后验概率问题来获取,而GMM中心点集按照一定的形变模型朝数据点集进行连续地漂移运动,在此过程中点集拓扑结构保持不变。

给定两个点集X和Y。设Y=yj|yj∈R2,j=1,2,…,m为GMM的点集的m个中心,X=xi|xi∈R2,i=1,2,…,n为由GMM产生的数据点。对数据点x_i可以由GMM通过下式产生:

p(xi)=∑j=1mpijp(xi|j),(11)

式中p_ij(i=1,2,…,n;j=1,2,…,m)表示第i个数据点的第j个混合系数。x_i的第j个高斯分布的概率密度函数为

p(xi|j)=1(2πσ2)D2exp-‖xi-yj*‖22σ2,(12)

式中σ为标准差,D为点集的维数。

设C(Y)为Y中点的位移函数,即任意取y_j∈Y, yj*=y_j+C(y_j)为第j个高斯因子的新位置。为了解释异常点及噪声,需要在原有GMM中额外增加一个均匀分布,总的数据分布为

p(xi|Θ,w)=w∑j=1mpijp(xi|j)+(1-w)p(Voutlier),(13)

式中p(V_outlier)=1/n表示异常点和噪声点的概率模型。p_ij(i=1,2,…,n;j=1,2,…,m)表示第i个数据点的第j个混合因子。Θ=C(Y),σ表示待估参数。

假设数据点是独立同分布,则可以得到的似然函数为

L(Θ,w)=∏i=1np(xi|Θ,w)。(14)

定义的先验分布为

p(C)=exp-λ2CH2,(15)

式中‖C‖ H2表示C(Y)在再生核Hilbert空间的模。理想形式可以表示为核的线性组合C(z)=∑j=1mφ_jK(z,y_j)=KΦ,其中K为一个m×m核矩阵,K(y_i,y_j)=exp -‖yi-yj‖22β。于是可得后验分布为

p(Θ,w|X)∝L(Θ,w)p(C),(16)

式中参数Θ和w可以通过求解最大后验概率得到。它等价于求解的最小化的负对数后验概率为

lnL(Θ,w)=-∑i=1nlog[w∑j=1mpijp(xi|j)+(1-w)p(Voutlier)]+λ2CH2。(17)

参数Θ和w的最小化可以通过期望最大(EM)算法来实现。

E步: 计算混合因子后验概率为

p(j|xi,Θ,w)=w(pij)/[(2πσ2)D2]exp[-(‖xi-yj*‖)/(2σ2)]w∑k=1mpikp(xi|k)+(1-w)/n,(18)

异常点的后验概率为

p(Voutlier|xi,Θ,w)=1/nw1-w∑k=1mpikp(xi|k)+1。(19)

M步: 通过最小化(18)和(19)式来更新参数Θ和w,通过求参数的偏导并令其为0可得

σ2=∑i,j=1n,mpoldjxi,Θ,w)‖xi-yj*‖2D∑i,j=1n,mpoldjxi,Θ,w),(20)

式中w=[ ∑i,j=1n,mp^old(j|x_i,Θ,w)] /[∑i,j=1n,mpoldjxi,Θ,w)+∑i=1npold(Voutlier|xi,Θ,w)]。类似地,可得:

[diag(P1)K+λσ2]Φ=PX-diag(P1)Y,(21)

式中diag表示矩阵对角元素,1为所有元素都为1的列向量,P为一个m×n矩阵,其元素p(j,i)=p(j|x_i,Θ,w)。

3.2 非均匀GMM

在(13)式表示的模型中,GMM中数据点的权重p_ij采用了均匀分布,即无论这个点是否具有精确的对应点,其混合因子权重p_ij都采用的如下确定的权值:

pij=(1-ω)/m。(22)

这种均匀的混合因子将导致GMM只利用了特征点的空间位置信息,使得匹配概率矩阵中的概率由于异常点的增加而将显著变小,降低了匹配模型对异常点的稳健性和匹配的精度。为了使匹配概率更好地反映点的匹配情况,应让那些更有可能具有精确的匹配点对在GMM中的混合系数中占据主导位置。

如果在GMM中的中心点y_j和数据点x_i的特征空间越相似,那么它们在模型中对应的权重系数p_ij就应该越大,这样也将使得具有相似特征的点更有可能成为匹配点对。因而可以将轮廓上对应点的特征空间的相似性作为混合模型权重系数。对于两个相匹配的点对,包含该两点的对应的窗口应该具有正相关性。因而可以采用标准化互相关系数^[19]作为相似性度量来构造非均匀GMM中的权重,其形式为

pij=exp-αρij∑k=1mexp-αρik,(23)

式中α为控制参数,ρ_ij表示以第i个特征点和第j个特征点为中心的两个窗口的标准化互相关系数,可以表示为

ρij=∑m=-NN∑ n=-NNIi(m,n)Ij(m,n)∑m=-NN∑ n=-NNIi2(m,n)∑m=-NN∑ n=-NNIj2(m,n),(24)

式中I_i(m,n)表示以第i个特征点为中心的位置为(m,n)的像素的强度。由此得到的非均匀GMM将特征点的位置信息和相似性信息统一到一个匹配模型之中,从而使匹配概率p(j|x_i,Θ,w)能更正确地反映特征点之间的匹配程度。

在(18)式中,p(j|x_i,Θ,w)表示数据点x_i和中心点y_j之间的匹配概率,概率值越小则该匹配越不可靠。为了得到精确的匹配点对,可以选择阈值T,当匹配概率小于阈值T时认为该匹配是误匹配而将该匹配点对剔除,从而得到精确的匹配点对。

4 实验结果

分别对合成数据和真实SAR数据进行图像分割与轮廓匹配实验来对所提出的方法进行验证。首先利用合成图像和真实SAR图像验证所提出的分割模型的有效性,然后利用所提出的分割模型提取SAR图像中水体目标的轮廓,并利用所提出的非均匀混合模型对轮廓上的点进行匹配,进而实现含水体目标的SAR图像配准。

4.1 图像分割实验

由于在利用水体目标轮廓进行配准时主要考虑分割方法对目标细节保持的能力和目标边缘定位的精度,因而图像分割实验只考虑目标和背景的分割实验。图像分割实验中对比了分割的水平集方法^[6]、可变权MRF(V-MRF)模型^[9]、自适应权MRF(A-MRF)模型^[10]和所提出的结合观测图像局部统计信息的A-MRF分割模型。基于MRF模型的初始分割由K-means方法得到。图1(a)为像素强度均值分别为90和130的剪切板在加入方差为0.05的乘性斑点噪声后的图像。为了显示分割模型保持图像细节的能力,剪切板中加入了细线。图1(c)~(f)分别表示水平集方法方法、V-MRF模型、 A-MRF模型和所提出的方法的分割结果,其正确率分别为88.61%,93.35%,93.87%,94.54%。由加噪合成图像的分割实验结果可以看出,所提出的方法能较好地保持图像边缘和细节。

图 1. 加噪合成图像分割的实验结果。(a)加噪合成图像;(b) K-means方法得到的初始分割结果;(c)水平集方法分割结果;(d) V-MRF模型分割结果;(e) A-MRF模型分割结果;(f)所提方法分割结果

Fig. 1. Segmentation experiment results of synthetic image with noises. (a) Synthetic image with noises;(b) initial segmentation result with K-means method; (c) segmentation result with level set method; (d) segmentation result with V-MRF model; (e) segmentation result with A-MRF model; (f) segmentation result with the proposed model

下载图片 查看所有图片

图2显示了三种方法对真实SAR图像的分割结果。由于在SAR图像中水体目标和非水体目标边缘的灰度具有比较大的反差,因而基于V-MRF模型、A-MRF模型和所提出的方法都能较好地分割水体目标,但前两种方法包含更多的噪声点,而且在图像的左上角将部分非水体区域误分为水体目标。而水平集方法具有比较光滑的分割结果,但不能很好地保持目标细节和定位边缘位置,同时,在图像的左上角包含较大的误分区域。

图 2. 真实图像分割实验结果。 (a) SAR图像; (b) K-means方法得到的初始分割结果;(c)水平集方法分割结果;(d) V-MRF模型分割结果; (e) A-MRF模型分割结果;(f)所提方法分割结果

Fig. 2. Segmentation experiment results of real image. (a) SAR image; (b) initial segmentation result with K-means method; (c) segmentation result with level set method; (d) segmentation result with V-MRF model; (e) segmentation result with A-MRF model; (f) segmentation result with the proposed model

下载图片 查看所有图片

4.2 含水体SAR图像的轮廓提取及轮廓匹配实验

本节实验针对5·12汶川地震前后所获取的包含河流及水库区域的SAR图像进行配准。图像来自日本ALOS卫星的L波段单视HH极化SAR传感器,幅值图像的分辨率为6.5 m。 SAR图像的拍摄时间分别为地震前的2008年2月17日和地震后的2008年5月19日。图3(a)、 (b)和图4(a)、(b)分别为四川绵阳地区的金华水库和狮儿河水库。从图像可以看到,地面场景由于受到了地震的毁坏而出现了局部扭曲变形及灰度变化,这增加了精确提取非水体部分的同名特征点的难度。SAR图像中的水体也受到了地震的影响,如金华水库入水口上游出现了比较严重的堵塞情况,水库部分位置由于水位变化或堤岸边缘出现山体滑坡及山上滑落的土石堆积而发生变化。这使得提取的水体轮廓部分位置出现了变形,同时出现了大量的外点[见图3(c)、 (d)和图4(c)、 (d)]。但是由于水体的特征,目标边缘和表面形状能基本保持稳定,部分边缘位置并没有发生改变,这为利用水体目标轮廓对SAR图像进行配准提供了可能。

图 3. 金华水库含水体的SAR图像配准。 (a)待配准图像; (b)参考图像; (c)~(d)提取的水体目标轮廓;(e)轮廓匹配结果;(f)配准后的图像;(g)变化检测结果

Fig. 3. Registration of SAR image including body of water in Jinhua reservoir. (a) Sensed image; (b) referenced image; (c)-(d) target contours of extracted water body; (e) matching result of contours; (f) result after registration; (g) change detection result

下载图片 查看所有图片

图 4. 狮儿河水库含水体的SAR图像配准。 (a)待配准图像;(b)参考图像; (c)~(d)提取的水体目标轮廓; (e)轮廓匹配结果;(f)配准后的图像;(g)变化检测结果

Fig. 4. Registration of SAR image including body of water in Lion Creek reservoir. (a) Sensed image; (b) referenced image; (c)-(d) target contours of extracted water body; (e) matching result of contours; (f) result after registration; (g) change detection result

下载图片 查看所有图片

在利用相干点漂移算法进行匹配时,用于计算标准化互相关系数的窗口大小为5×5。(23)式中控制参数α=0.1,用于剔除误配的阈值取匹配概率矩阵的十分之一分位点。通过多次迭代,直到轮廓上匹配点对关于位置的均方误差(RMSE)小于1或匹配点数等于3。其中描述两个点集匹配精度的均方根误差公式为

VRMSE=1N∑i,j(xi,yj)-(x',y')ij‖2。(25)

图3(e)和图4(e)为利用所提方法得到的匹配结果。可以看到,通过剔除匹配矩阵中概率较小的点对,使得轮廓上的大量的点对被剔除,只保留了部分精确的匹配点对,从而得到了较好的配准图像[图3(f)和图4(f)]。 图3(g)和图4(g)表示变化检测结果。其中红色部分表示地震后水体目标减少的区域,即地震后被山体滑坡挤压的区域或山石滚落后掩埋的区域,蓝色部分为震后水面增加的部分。

5 结论

提出了针对含水体目标的SAR图像配准的方法。该方法主要在两个方面进行了改进:一方面,在轮廓特征提取时,提出了一个包含观测图像局部信息的MRF分割模型,该模型能较好地保持目标细节和精确地定位目标边缘;另一方面,在轮廓匹配阶段,给出了一种非均匀系数的GMM,该模型同时融合了轮廓上特征点的位置信息和以特征点为中心的窗口的统计相似性特征,能提高轮廓的匹配精度。利用模拟图像和真实图像进行轮廓的提取和匹配实验,实验结果验证了模型的有效性。