结合双树复小波纹理特征和MRF模型的遥感图像分割

韦春桃，赵 平，2，肖博林，白 风，李小勇，杨晚芸，3

(1.重庆交通大学，重庆 400074；2.中国电建集团贵州电力设计研究院有限公司，贵州 贵阳 550002；3.中核咨询有限公司四川分公司，四川 乐山 614000)

随着航空航天技术的发展，高分辨率影像中各类地物细节变得清晰，纹理特征呈现多样性，但同时也增加了同质区域内部光谱的差异[1]。为弥补光谱特征应用于影像分类的不足，由光谱、纹理形成的混合特征被广泛应用于高分辨率遥感影像分割方法中且取得了一定的成效[2-4]。

Gabor纹理[5]、LBP纹理[6]、灰度共生矩阵[7]、小波纹理[8-9]等方法在遥感影像分割中都取得了很好的应用效果。文献[9]弥补了光谱特征的不足，运用小波纹理特征进行高分辨率遥感影像分类。由于二维小波变换(2D DWT)存在平移变动性和有限方向选择的缺陷，文献[10]提出了具有近似平移不变性和良好方向选择的双树复小波变换(dual-tree complex wavelet transform,DT-CWT)；文献[11]将双树复小波纹理应用于遥感影像土地覆盖分类，但该方法缺少对光谱特征的利用，影像信息表达不够全面。马尔可夫随机场(Markov random field,MRF)模型由于具有能引入影像上下文关系的优势在遥感图像分类中有很好的应用[12]。

针对高分辨率遥感影像分割中小波纹理特征表达能力不足，像元的空间关系和分布特征利用不充分的缺陷，本文提出结合复小波纹理和MRF模型的高分辨率遥感图像分割方法。首先，使用双树复小波变换提取影像纹理特征，并与光谱信息形成混合特征；然后,将混合特征进行高斯归一化处理，并用K-means算法对混合特征聚类，得到高分辨率遥感影像的初始分割；最后，由MRF模型引入上下文信息，基于贝叶斯最大后验概率准则(maximum a posteriori,MAP)优化初始分割结果，减弱分割图中同质区域的“椒盐噪声”，得到最终的分割结果。

1 双树复小波变换原理

双树复小波由2组并行小波树构成。如图1所示，在2组小波树Tree-a与Tree-b上进行实数小波变换，并将2组树的输出分别称为双树复小波变换的实部和虚部。一维双树复小波变换为

Ψ(t)=Ψa(t)+jΨb(t)

(1)

二维双树复小波函数定义为

Ψ(x,y)=[ψa(x)+jψb(x)]×[ψa(y)+jψb(y)]

(2)

如图2、图3所示，双树复小波实部和虚部均有6个方向信息，分别为+75°、+45°、+15°、-15°、-45°、-75°。

双树复小波分解如图4所示，其中图4(a)为原始影像，图4(b)、图4(c)表示获得的1个低频子影像和6个高频子影像，高频子图对应6个方向信息。

2 MRF图像分割

MRF描述像素点间的空间关系，在贝叶斯框架下做最优分割。设待分割影像大小为s=m×n，每个像素赋值标号Y={Y|y1，y2，…，yi}，i为标号类别。每个像素特征矢量X={X|x1，x2，…，xs}。图5为标记场与特征场的对应关系。根据贝叶斯准则将图像分割表示为

(3)

由最大后验概率准则(MAP)使后验概率取得P(Y=y|X=x)最大值[13]，影像的最优分割结果为

(4)

2.1 建立标号场模型

标号场用二阶邻域系统的MRF表示，图6所示为二阶领域及二阶领域内双点势团类型，根据Hammersley-Clifford理论可知，先验概率服从吉布斯分布[14]

(5)

式中，U(X=x)为能量函数；Z为配分函数；T为温度，通常设为1。

(6)

Vc为基团的势函数，本文使用Potts模型[15]定义二元势函数，其中β取值为0.5。

(7)

此时Potts局部概率模型可表示为

(8)

式中，Ni(yi)为二阶邻域系统中yi≠yi的个数。

2.2 建立特征场模型

假设影像中不同类别的观测特征相互独立，每个像元相互独立且服从高斯分布，则一幅影像的概率分布可表示为

P(X=x|Y=y)=

(9)

2.3 MAP估计与求解

根据式(4)、式(8)、式(9)，在给定特征场和标记场的条件下，基于最大后验概率准则(MAP)，图像分割的全局最优估计可表示为

(10)

本文采用条件迭代模型(ICM)算法[15]找到使能量最小的条件，迭代更新出最优标记场(Y)～，从而得到图像分割的最佳结果。

3 分割算法步骤及试验

3.1 分割算法步骤描述

(1)用中值滤波对待分割遥感影像预处理，去除由成像传感器性能产生的噪声点。在大小为L×L的滑动窗口内做双树复小波纹理变换，取得复小波纹理特征。本文窗口大小选择为15×15。窗口内的像元做一层双树复小波变换，得到一个低频子带系数，并产生6个高频子带系数hl；l表示高频子带的方向(l=-15°、-45°、-75°、+15°、+45°、+75°)。本文使用每个子带系数的均值与方差作为像元的纹理测度。

(11)

式中，ul、vl分别表示6个方向高频子带的均值与方差。

(2)把原始影像中每个像元的r、g、b数值作为光谱特征，并与步骤(1)获得的纹理特征结合形成混合特征向量

c=[ulvlrgb]

(12)

由于混合特征向量中的纹理特征分量和光谱特征分量具有不同的物理意义，为使特征向量中的元素处于同一量纲下，本文使用高斯归一化处理混合特征向量，使特征向量内的元素服从均值为0、方差为1的标准正态分布。

(13)

式中，U、V为特征向量c的均值与方差；E是维数与c向量维数相同的单位向量。使用K-means算法对高斯归一化后的混合特征向量C进行聚类，得到影像的初始分割结果Y。

(3)初始分割结果Y表示影像的标记场，综合原始影像各地物的分布特征和初始分割结果建立MRF模型，通过条件迭代算法(ICM)求得MAP准则下影像的最优分割结果。

综上，本文所提基于双树复小波纹理与上下文信息高分辨率遥感影像分割算法步骤如图7所示。

3.2 试验结果与分析

在2幅高分辨率遥感影像上验证本文算法的有效性，并将该算法与基于光谱特征、基于光谱特征+MRF、基于小波纹理特征的结果进行比较。试验算法通过Matlab 2014b软件平台编写，在Intel core i5处理器，内存4 GB，CPU 1.6 GHz环境下运行。

试验1使用的待分割图像为高分二号遥感影像，如图8(a)所示。原始影像大小为256×256，其中光谱空间分辨率为1 m，包含R、G、B 3个波段信息，图像中包含树林、草地、浅滩、水域等地物。图8(b)—(e)分别为4种方法的分割结果图。试验1中截取了原始影像中4种地物进行精度评定，结果见表1。4种地物在不同方法下得到的分割精度，由光谱特征得到的分割结果中在浅滩和绿地种类型地物上的分割精度不高。由于地物内存在其他次要的光谱信息，图8(b)相较于图8(c)、(d)、(e)在同质区域产生了很多的“椒盐噪声”，而其他3种分割方法考虑像元间的空间关系和分布特征获得了很好的区域一致性。本文方法在江面的分割精度上略低于其他方法。但浅滩和树林两种地物的分割精度分别为54.6%和84.9%，对比光谱+MRF方法、小波纹理方法树林的分割精度提高很大。绿地部分的分割精度为86.8%。仅次于小波纹理88.3%的分割精度。

为进一步验证算法的有效性，试验2使用Google Earth遥感影像作为待分割影像，如图9(a)所示。影像大小为256×256，其中光谱空间分辨率为0.5 m，包含R、G、B 3个波段信息，图像中有林地、水域、房屋、耕地、草地、道路等地物。试验2中增加了影像的地物种类，地物之间存在更加复杂的关系。图9(b)—(e)分别为4种方法的分割结果图。如图9(b)所示，分割出来的区域比较细碎，噪声多，房屋、耕地、道路3种地物的分割精度不高，试验精度(见表2)分别为37.8%、39.5%和8.5%。光谱+MRF方法虽然比仅仅使用光谱特征的分割方法获得了很好的区域一致性，但由于特征向量判断能力不足，在房屋和耕地区域的分割精度也不高。小波纹理的分割方法在视觉上获得了很好的区域一致性，但也存在特征向量判断力不足的问题，分割影像在房屋、耕地、水域误分割较多。本文方法增强了纹理特征的表达能力，同时考虑了分割影像像元间的空间关系和分布特征，虽然在林地和水域的分割精度稍低于其他方法，但在房屋、道路、耕地等复杂区域，本文分割方法得到了很好的视觉效果和精度。

表1 试验1精度统计

4 结 语

本文针对小波纹理特征在方向选择上的不足和平移变动性的缺点，使用双树复小波纹理，增强了纹理特征的表达能力。将光谱特征和纹理特征应用于图像分割。使用MRF引入上下文信息，建立初始分割图与原始影像间的MRF模型，在MAP准则下优化分割结果，分割结果图中同质区域的“椒盐噪点”得到有效的抑制，获得了更好的区域一致性，提高了影像分割的精度。试验表明本文算法与基于光谱特征的影像分割方法、基于光谱特征+MRF影像方法及基于小波纹理特征的影像分割方法相比，获得了更高的精度。

表2 试验2精度统计