基于马尔科夫随机场的图像分割算法一直以来被国内外学者广泛研究, 并在各领域广泛应用(李旭超等, 2007; Li S Z, 2009; Li Y et al., 2011; 温航宇, 2013; 徐迪, 2013; 徐胜军等, 2013; Zheng et al., 2013; 李铁, 2014; Mignotte, 2014; 宋艳涛等, 2014; 朱瑶, 2014; 曹家梓等, 2015; 罗欣, 2015; 许妙忠等, 2015; 杨鹏, 2015)。 该模型具有以下几个鲜明的特点: 1)空间约束性强, 将像素的邻域系统与先验分布模型相结合, 能够有效地描述空间像素与像素之间的依赖关系; 2)既能反映图像的随机性, 又能反映图像的潜在结构, 可有效地描述图像的性质; 3)既从物理模型出发, 又与图像数据(灰度值或特征)拟合直接联系起来。 MRF模型运用Gibbs场和最大后验概率(Maximum A Posteriori, MAP)实现图像分割, 其过程可概述如下:

假定S是n个位置的像素集合, S={1, 2, …, n}, 其中{1, 2, …, n}代表位置索引, 不同的位置索引表示不同的像素或特征。 将一个大小为M× M的图像作为二维方形格网, 则S可表示为S={(i, j)|1≤ i≤ M, 1≤ j≤ M}。L和D为2个符号集, L={1, 2, …, l}, D={1, 2, …, d} , X和Y为2个随机场, 分别对应于图像的类别标识和灰度, 它们的状态空间分别是L和D, 即对于∀ i∈ S, X_i∈ L, Y_i∈ D。 在集合S中可用邻域系统来表示每个像素之间的关系, 二维网格S中邻域系统的定义是N={N_i|∀ i∈ S}, 其中N_i代表位置i处的邻域集合。

根据贝叶斯估计, x的后验概率可以表示为

P(x|y)=P(y|x)P(x)P(y)(1)

其中, y代表观测图像, P(y|x)代表在给定x的情况下y的条件概率, P(x)代表x的先验概率, P(y)代表观测密度。

令$\hat{x}$是图像真实类别标号的估计, 由最大后验概率(MAP)准则有:

$\hat{x}$=argmaxx∈M{P(y|x)P(x)}(2)

其中, x表示随机场X的一组配置, M是X所有可能配置的集合, y表示随机场Y的一组配置。要得到 $\hat{x}$, 首先需计算类别的先验概率P(x)和观测量的似然概率P(y|x)。

由Hammersley-Clifford定理可知, MRF能够等价描述为Gibbs分布; x可视为MRF的具体实例, 其先验概率近似表示为

P(x)∝exp(-U(x))(3)

能量函数 U(x)=-β∑c∈CVc(x),为基团C的势能V_c(x)之和。 引入Potts模型有V_c(x_i)=δ (x_i, x_j)-1, j∈ N_i, 则:

u(xi)=-βi∑c∈CVc(xi)=βi∑j∈Ni(δ(xi, xj)-1), 其中δ(xi, xj)=1xi=xj0xi≠xj(4)

式中, x_i代表当前分割结果, x_j代表真实分割结果。 对比当前的分割结果和真实分割结果, 若像素之间的标记不同, 则δ (x_i, x_j)取1, 反之取0。 MAP准则的目的是将图像中每个像素的错分概率降到最低。

而对于给定类别标号x_i=ℓ, 通常假定像素颜色深度值y_i服从高斯分布θ _i={μ _ℓ, σ _ℓ}, 即:

P(yi|xi=l)=(1/2πσl2)·exp(-(yi-μl)2/2σl2)(5)

常用的MRF-MAP求解方法有很多, 通常通过迭代搜索能量函数最小的解作为建筑物的分割结果。 本文选择条件迭代(Iteration Condition Model, ICM)算法进行迭代运算, 以求出最优解。

迭代条件模型是基于局部条件概率求解局部最小值的算法, 该模型通过对每个像素点进行逐点扫描来更新图像的标记, 经过数次迭代以寻找局部最优或全局最优的建筑物分割结果。 其公式为

$\hat{x}$=argmaxxi∈Lf(yi|xi)·f(xi|xNi), i∈S(6)

由于ICM算法的分割效果较为依赖初始分割, 因此选择模糊k-means聚类算法求取ICM迭代算法的初始分割。 ICM算法的具体流程为:

(1)影像读取, 根据影像中目标的特征设置需要分割的类别数K、 势函数β 和迭代次数。

(2)输入影像, 利用k-means分类方法得到影像初始分割结果。

(3)针对特征场数据, 选取高斯混合模型建模, 计算出特征场参数均值μ 和方差σ ²。

(4)对特征场的能量进行求解。

(5)根据标记场数据, 选取Potts模型进行建模, 并求出标记场的能量函数。

(6)计算全局能量的值, 通过迭代搜索求出能量函数最小的解, 得到新的分割结果。

(7)判断是否满足收敛条件。 若满足条件, 迭代结束, 得到最终的分割结果; 若不满足条件, 则转到步骤(3)。

随着遥感影像空间分辨率的提升, 图像信息也更加丰富。 建筑物的典型特征是具有明显的外形轮廓和规律的排列形式, 在遥感影像上表现为灰度及其结构具有规则性和有序性。 地震使建筑物发生不同程度的倒塌, 建筑物的规则结构也受到了不同程度的破坏, 主要体现在波谱特征和纹理结构2个方面。 在波谱特征方面, 完好的建筑物往往呈现出高亮度的特征, 与倒塌建筑物、 半倒塌建筑物及其它地物形成鲜明对比, 因而可凭借影像区域的明暗程度对受损程度不同的建筑物进行区分。 在纹理结构方面, 完好的建筑物与受损建筑物的纹理特征具有明显差异。 完好建筑物的几何结构较规整, 且表面往往呈现出均一的图像纹理; 而对于倒塌建筑物或者半倒塌建筑物, 由于破坏截面比较粗糙, 在图像上通常呈现出斑块状的低灰度值区域(柳稼航, 2003)。 因此, 可以依据影像中建筑物的纹理特征来区分震害影像中的完好建筑物和受损建筑物。

灰度共生矩阵(Gray Level Co-occurrence Matrix, GLCM)是描述影像灰度值变化的二阶统计量, 它体现了影像中每个像元的灰度在方向、 邻域间隔以及变化幅度上的特征, 是研究影像局部区域特征及排列属性的基础。 本文在大量实验后, 从中选取了3种纹理特征参数参与分类, 分别是对比度(Contrast)、 非相似性(Dissimilarity)及方差(Variance)。 研究发现, 完好建筑物的对比度较大, 非相似性与对比度线性正相关, 方差值较小; 而受损建筑物对比度较小, 非相似性较低, 方差较高。

对比度反映了像素与其周围像素灰的对比情况。 其计算公式为

C=∑i=0n-1∑j=0n-1p(i, j)(i-j)2(7)

其中, i和j分别代表灰度共生矩阵中行、 列的灰度值, p(i, j)代表第i行、 第j列的标准化值, n代表行或列的总数。

非相似性的计算公式为

D=∑i=0n-1∑j=0n-1p(i, j)|i-j|(8)

方差的计算公式为

V=12∑i=0n-1∑j=0n-1[p(i, j)(i-μ)2+p(i, j)(j-μ)2](9)

其中, μ 代表均值。

支持向量机(Support Vector Machine, SVM)是基于统计学习理论而发展起来的一种新的学习方法, 在遥感影像分类、 人脸检测与识别、 回归分析、 图像分类和模式识别等领域都得到了广泛的应用(Weston, 1999; 杜培军等, 2006; 谭琨等, 2008; 刘晓亮等, 2010; Niu et al., 2012; 高晓健, 2013; Vijayan et al., 2016; Dilmen et al., 2017; Ebrahimi et al., 2017)。 支持向量机具有众多优越性, 该模型基于结构风险最小的原则, 具有较高的泛化能力, 不仅学习速度快, 在高维数据的处理上更是有着独特的优势。

支持向量机是一种二分类模型, 其目的是寻找一个超平面对样本进行分割, 最终转化为一个凸优化问题, 以获得全局范围内的最优解。 该模型遵循以下原则: 对于线性可分的情况, 直接进行分类; 对于线性不可分的情况, 利用非线性映射核函数(如线性核函数、 多项式核函数、 径向基函数等)将输入向量映射到一个高维甚至无穷维的特征空间, 在特征空间中构造最优分类超平面, 使原本非线性不可分的问题转换为线性可分的问题, 即找到一个超平面, 在该平面内训练样本的样本点被分为2类, 且2类样本点与该超平面的距离之和最大。

如图 2所示(以二维情况为例), 图中粉色和蓝色的点分别代表2类样本, 超平面A和超平面B是2种分类超平面, 其中超平面B为最终的分类超平面, 因为超平面B不仅能将2类样本点正确区分开, 而且使得2类样本点与超平面的距离之和最大。

图 2

Fig. 2

	Figure Option ViewDownloadNew Window
	图 2 最优分类超平面Fig. 2 Optimal classification hyperplane.

根据支持向量机理论, 其最优分类超平面为

f(x)=sgn[∑i=1naiyiK(xi, x)+b](10)

其中, a_i> 0, 为拉格朗日乘值, 满足a_i> 0的样本点代表支持向量; y_i为结果类别标签, 通常取-1或1; K(x_i, x)是支持向量机中的核函数, 满足K(x_i, x)=< φ (x_i)· φ (x)> ; b为分类域值。

支持向量机对二分问题具有较好的分类效果, 对于多类分类问题, 一般采取多次二分的方法实现(Debnath et al., 2004; 郭显娥等, 2010; 陈中杰等, 2013; 刘绍毓等, 2015)。

采用分辨率为0.5m的ADS40数字航空影像作为实验数据, 以2008年5月12日发生M8.0地震的四川省汶川县作为研究区域。 汶川M8.0地震的震源深度为14km, 其将汶川县的大部分建筑物摧毁。 从2008年5月16日汶川县ADS40震后遥感影像中截取了5块子区域, 同时, 为了验证本文方法对于不同分辨率影像的受损建筑物的提取能力, 从同等设备获取的具有不同分辨率的影像中选取了1块子区域。 实验区的详细信息见表1。

表1

Table1

表1(Table1)

表1 实验子区域概况 Table1 Overview of subareas 实验区域 影像尺寸/像素 影像空间分辨率/m 实验区地貌景观 T1 696× 806 0.5 建筑物, 植被, 山体滑坡, 河流, 裸地 T2 488× 298 0.5 建筑物, 植被, 道路, 裸地 T3 480× 139 0.5 建筑物, 植被, 道路, 山体滑坡 T4 409× 385 0.5 建筑物, 植被, 裸地, 道路 T5 298× 319 0.5 建筑物, 植被, 道路, 裸地, 滑坡 T6 166× 206 1 建筑物, 植被, 裸地

表1 实验子区域概况 Table1 Overview of subareas

实验包括2个部分: 1)利用本文提出的算法对所选取的6个实验子区域进行实验与精度评定; 2)与面向对象分类方法所得的结果进行比较研究, 以证明本方法的优越性。

2.2.1 本文方法的实验结果与分析

运用本文方法对6个子区域中不同程度的受损建筑物进行提取, 结果如图3— 8所示。 此方法需设置的参数包括: 1)影像分类的类别数; 2)势函数的取值, 本方法利用Potts模型构建势函数; 3)算法的迭代次数。

从图3— 8可以看出, 本文所提出的方法能够较好地提取出完整建筑物, 在提取半倒塌建筑物及完全倒塌建筑物时, 本方法所得结果与真实结果相比面积偏大, 这是由于完全倒塌建筑物部分与裸地部分极为相似; 此外, 半倒塌建筑物与完全倒塌建筑物之间存在少许错分现象, 这是由于在半倒塌建筑物中, 受损的部分与倒塌建筑物的特征十分相似, 故本方法将一些半倒塌建筑物的受损部分识别为倒塌建筑物。

图 3

Fig. 3

	Figure Option ViewDownloadNew Window
	图 3 T1的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 3 Original image(a), classification results(b) and groundtruth(c) of T1.

图 4

Fig. 4

	Figure Option ViewDownloadNew Window
	图 4 T2的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 4 Original image(a), classification results(b) and groundtruth(c) of T2.

图 5

Fig. 5

	Figure Option ViewDownloadNew Window
	图 5 T3的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 5 Original image(a), classification results(b) and groundtruth(c) of T3.

图 6

Fig. 6

	Figure Option ViewDownloadNew Window
	图 6 T4的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 6 Original image(a), classification results(b) and groundtruth(c) of T4.

图 7

Fig. 7

	Figure Option ViewDownloadNew Window
	图 7 T5的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 7 Original image(a), classification results(b) and groundtruth(c) of T5.

图 8

Fig. 8

	Figure Option ViewDownloadNew Window
	图 8 T6的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 8 Original image(a), classification results(b) and groundtruth(c) of T6.

2.2.2 面向对象方法的实验结果与分析

eCognition是目前广泛应用的面向对象分类方法的智能化影像分析软件。 为了验证本文方法的有效性和优越性, 利用eCognition对研究区进行实验, 并将其结果(图9— 14)与本文所提出方法的结果进行对比。 eCognition软件提取震害信息主要分为3个步骤, 即影像分割、 特征提取和影像监督分类。 影像分割是指基于颜色、 形状和灰度等性质将影像分割成具有不同特征且互不相交的区域, 每个区域内部的灰度、 纹理等特征相近, 而相邻区域之间不具有一致性特征。 由于本文主要研究建筑物的损毁情况, 因此实验中采取层次分割的方法进行分割, 首先划分出建筑物区域和非建筑物区域, 再在建筑物区域内提取出完好建筑物和损毁建筑物。

震后遥感影像中建筑物的大小、 形状各异, 且不同影像、 不同区域中建筑物的密度、 类型也不尽相同, 在适宜的尺度上进行震害信息提取非常重要, 因此, 2层分割均采用多尺度分割方法。 影像分割完成后, 根据不同建筑物以及建筑物与周围背景的差异, 选取了亮度值、 长宽比、 形状指数和灰度共生矩阵等指标构建特征规则集。 面向对象分类的最后一步是影像监督分类, 进行分类之前, 在影像中随机选择一定数量的样本区域, 再从这些区域中选取一部分作为训练样本, 剩下的区域作为检验样本, 采用最邻近距离法进行分类。

从图9— 14可以看出, 利用面向对象分类法提取震害影像中的受损建筑物具有较好的效果。 对于完整建筑物和半倒塌建筑物, 面向对象的分类方法能够较好地对其边界轮廓进行描述, 保存建筑物的边界信息, 这是由于面向对象分类方法在分类前先将具有相同特征的像素聚为图像块, 从而使分类结果均以 “ 块状” 进行呈现, 很少出现细小的噪声点。 但是, 面向对象的分类方法对于像素灰度值具有一定的依赖性, 使得部分灰度值较低的完整建筑物被分为半倒塌建筑物, 如图中红色圆圈所示。 同时, 对于完全倒塌建筑物, 面向对象的分类方法不能将其与土地进行有效区分, 如T₁的影像所示。 除此之外, 利用面向对象分类法对震后遥感影像中的建筑物进行提取时需要设置大量的参数, 且分类结果极易受到操作者主观性的影响, 因而不能快速准确地获取震害信息。

图 9

Fig. 9

	Figure Option ViewDownloadNew Window
	图 9 T1的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 9 Original image(a), classification results(b) and groundtruth(c) of T1.

图 10

Fig. 10

	Figure Option ViewDownloadNew Window
	图 10 T2的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 10 Original image(a), classification results(b) and groundtruth(c) of T2.

图 11

Fig. 11

	Figure Option ViewDownloadNew Window
	图 11 T3的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 11 Original image(a), classification results(b) and groundtruth(c) of T3.

图 12

Fig. 12

	Figure Option ViewDownloadNew Window
	图 12 T4的原始影像(a)、 分类结果(b)及真实分类结果(c) 红色圆圈表示完整建筑物被错分为半倒塌建筑物的部分Fig. 12 Original image(a), classification results(b) and groundtruth(c) of T4.

图 13

Fig. 13

	Figure Option ViewDownloadNew Window
	图 13 T5的原始影像(a)、 分类结果(b)及真实分类结果(c) 红色圆圈表示完整建筑物被错分为半倒塌建筑物的部分Fig. 13 Original image(a), classification results(b) and groundtruth(c) of T5.

图 14

Fig. 14

	Figure Option ViewDownloadNew Window
	图 14 T6的原始影像(a)、 分类结果(b)及真实分类结果(c)Fig. 14 Original image (a), classification results (b) and groundtruth (c) of T6.

将本文所提出的方法与eCognition的实验结果(图3— 14)进行对比可知, 首先, 由于受像素灰度值影响较小, 本文的方法与面向对象分类方法相比, 能够更好地区分完整建筑物和半倒塌建筑物, 从而使最终结果在精度上有较大幅度的提升; 其次, 在建筑物的形状轮廓描述方面, 虽然本文方法略逊于面向对象分类方法, 但也能较为准确地保存建筑物的边界信息; 最后, 本文的方法对于倒塌建筑物和土地也有较好的区分能力。 地震发生后, 救援的黄金时间只有短短的72h, 快速准确地获取灾情信息尤为重要。 而本文方法最大的优点就是只需设置少量参数(影像分类类别数、 势函数取值和算法迭代次数), 且运行速度快, 缩短了获取震害灾情信息的时间。

对比2种方法的混淆矩阵, 通过APA、 AUA、 OA及Kappa系数4种指标对震害建筑物分类精度进行评价, 结果见表2。

表2

Table2

表2(Table2)

表2 利用本文所提方法与面向对象对比实验精度结果 Table2 Comparison of accuracy between the method proposed in this paper and the object-oriented experiment 评价指标 实验方法 实验区域 T1 T2 T3 T4 T5 T6 APA/% 本文方法 93.2 92.8 92.9 95.4 91.4 91.6 eCognition 81.3 81.7 90.6 90.8 86.1 80.6 AUA/% 本文方法 94.2 91.3 95.8 94.3 91.2 92.9 eCognition 82.7 82.6 85.5 92.4 86.4 74.5 OA/% 本文方法 93.9 92.3 92.9 94.6 92.0 92.4 eCognition 80.6 80.8 89.3 91.4 86.7 85.2 Kappa 本文方法 0.916 0.895 0.900 0.918 0.892 0.842 eCognition 0.740 0.745 0.847 0.882 0.820 0.642

表2 利用本文所提方法与面向对象对比实验精度结果 Table2 Comparison of accuracy between the method proposed in this paper and the object-oriented experiment

从表2中可以看出, 本文的方法在提取完整和倒塌建筑物时准确性较高, 无论在何种评价指标下, 识别精度均在90%以上; 总体精度最高可达94.6%。 与面向对象方法相比, 本文所提出的方法在T₁影像中的总体分类精度提高了13.3%。 此外, 本文方法比面向对象方法的速度更快、 参数更少且更易调节。 因此, 利用该方法可在地震后快速提取建筑物的损伤信息, 为救灾和减轻灾害损失赢得时间。