利用增量式马尔科夫随机场分割提取高空间分辨率遥感影像道路

吕 野， 胡翔云

(武汉大学遥感信息工程学院,武汉 430072)

0 引言

随着信息产业的日益发展以及信息技术的不断提高，人类与地理信息技术联系越来越紧密。对地理信息准确、高效的获取成了众多信息学科工作者们共同的研究目标。空间数据的人工采集与人工目视判读等方式已经不适合当今的发展趋势，也难以满足人们的需求。如今，摄影测量与遥感技术和计算机视觉技术能够完成大批量地理空间数据的采集与处理分析，能够为当今的城市发展提供强有力的支持。从遥感影像数据中提取兴趣目标是遥感领域的重要任务，而道路提取则可为城镇发展建设需要提供服务。

道路在地理信息中占有重要地位，与交通、城市规划、应急响应等行业有密切联系。道路的识别与提取研究一直是近20 a来的一个重要研究课题，从城市地区中提取道路网更是尚未解决的难题。目前，在中低空间分辨率遥感影像中提取郊区道路已有不少可靠稳健的算法策略，但是，这些地区地物简单，场景简单，而且道路多呈线状； 而在高空间分辨率城镇影像中，地物种类多样，背景复杂，道路呈面状，路面还受道路线、车流、行人等因素的干扰，并且道路形状、种类多变，易受周围地物的干扰、遮挡，致使道路特征变化复杂，产生“同物异谱”、“异物同谱”的现象[1]。常见的道路提取难点包括以下几个方面： ①影像空间分辨率高，纹理特征复杂，影像噪声大； ②路面、周围高大地物阴影及车流遮挡； ③不同等级道路差异大，路面材质与道路宽度变化大； ④非道路条带状地物的干扰； ⑤道路拓扑关系复杂； ⑥道路连接方式多样，包括十字路口、丁字路口和高架桥等。

通常，道路提取工作依据如下5种主要特征[2]： ①几何特征。在高空间分辨率遥感影像中，道路表现为具有一定长度的带状目标，其边缘与中心线都具有明显的线状几何特征。道路的长度远大于其宽度，故长宽比非常大。通常道路的曲率变化很小，边界平滑。沿着道路方向，宽度变化也较小，一般道路的宽度不会突然剧烈变化，而在道路交叉处宽度变化则较大。②辐射特征。道路内部灰度变化均匀，且与其周边地区灰度反差较大，道路一般表现为黑、白、灰的颜色特征，但会受路面上的车辆、道路线和行人等因素的影响。③拓扑特征。道路存在的意义决定了它的拓扑特征。道路间一般是相互连通的，且不会突然中断，各部分道路会整体形成一个连通域，组成路网。④功能特征。道路与人类的生产生活密切相关，所以道路一般都有通达的目的地，会与村庄、城镇等居民地或人工设施相连接。⑤上下文特征。通常指与道路相关的地物所具有的特征信息，包括了局部上下文信息和全局上下文信息，如电线杆、建筑物、行道树、立交桥、高大地物的阴影、路面上的汽车、路中央的隔离带、交通管理线等。高空间分辨率影像中存在许多不同种类的上下文特征，可以被合理地用于提取道路。

利用上述特征，许多道路提取算法被相继提出，例如数学形态学方法[3]、模板匹配方法[4]、snake模型法[5]、动态规划与优化算法[6]、分割与分类方法[7]、边缘与直线检测法[8]、道路跟踪法[9]、知觉组织方法[10]和机器学习方法[11]等。为了能够更可靠地提取道路，从多种类型的数据中共同提取道路也成为热门的研究方向，例如激光雷达(light detection and ranging，LiDAR)点云数据与遥感影像的结合[12]，数字高程模型(digital elevation model，DEM)与遥感影像的结合等[13]。

本文将综合利用道路的辐射特征(道路的灰度信息)、拓扑特征(道路相互连接的特性)以及上下文特征(道路与背景的模型信息)进行有效的道路提取。

1 马尔科夫随机场影像分割

1.1 影像分割

本文将以马尔科夫随机场影像分割[14-15]方法为基础。影像分割问题可以根据如下描述定义。

影像是一个以灰度为值的数组z=(z1,…,zn…,zN)，n为编号；N为像素数。影像分割被表达为一个不透明度的数组α=(α1,…,αN)，与各像素相对应。通常，0≤αn≤1，但对于硬分割，αn∈{0,1}，0代表背景，1代表前景。参数θ用来描述前景和背景的像素灰度级别的分布，由灰度直方图构成，即

θ={h(z;α)},

(1)

一个用来描述前景，一个用来描述背景。灰度直方图直接通过相应的影像获得。灰度直方图被规范化，其和为1，即

(2)

影像分割则是通过z以及θ来估计α。

1.2 基于能量优化的影像分割

定义能量函数E，通过求能量函数的最小值来求取最优的影像分割，该能量函数由以下几个方面进行引导： 前景与背景灰度直方图，以及不同透明度的统一性，即同一个目标应该具有相同的不透明度。以上几点可通过Gibbs能量形式表达为

E(α,θ,z)=U(α,θ,z)+V(α,z),

(3)

式中： 数据项U为在给定灰度直方图θ的条件下，α符合z的程度；V为平滑项。U和V分别定义为

(4)

(5)

式中： 符号[Φ]为指示函数，对于任意变量Φ，函数值为0或者1；C为邻接的像素对； dis(.)为欧式距离函数。平滑项使相邻像素拥有相似的灰度。若η=0，则平滑项将使得所有区域更加平滑，而平滑程度取决于常数γ； 若η>0将会在高对比度的区域放松平滑的程度[14]。η的取值设置为

η=(2<(zm-zn)2>)-1,

(6)

式中<.>表示影像样本中的期望。通过对η设置不同的取值，可以调整影像整体高低对比度。

至此，能量模型就被完全定义了，影像分割可以通过求全局能量最小值获得，即

α=argminαE(α,θ)。

(7)

整体优化求解通过标准的最小割算法获得。

对于包含道路的影像，通过上述的硬分割方式，将道路与非道路划分开。由于直接通过灰度直方图估计影像前景和背景模型在效率和算法稳定性方面存在欠缺，本文将使用高斯混合模型来估计前景和背景区域的影像模型，然后利用此模型对各像素区域的数据项进行估计。

1.3 颜色建模

高斯混合模型是较为有效、稳定，且对颜色空间表征具有一定抽象性的模型[15]，对彩色影像建模具有较强的适应性。对于包含道路的影像，用一个高斯混合模型表征前景，即道路区域； 用另一个高斯混合模型表征背景，即非道路区域。每个高斯混合模型都由T个成分构成。在此，定义额外的向量K={K1,…,Kn,…,KN}，其中，Kn∈{1,…T}，对于每个像素，使用唯一一个高斯混合模型成分。该成分来自前景道路模型或者背景非道路区域模型，由αn=0或者αn=1决定。

由此，Gibbs能量函数变成

E(α,K,θ,z)=U(α,K,θ,z)+V(α,z),

(8)

式中K为选择的高斯混合模型的成分编号。U现在定义为

(9)

其中，

D(αn,Kn,θ,zn)=-lnp(zn|αn,Kn,θ)-lnп(αn,Kn)，

(10)

式中： 函数p(.)为高斯概率分布； п(.)为混合模型的权重系数。据此得

(11)

因此，模型的参数现在变为

θ={п(α,K),μ(α,K),Σ(α,K)}。

(12)

对于平滑项V，基本保持不变，仅将对比项改为颜色空间的欧氏距离，即

(13)

式中η依据影像信息获得，即

(14)

式中： mean为取平均值；m和n为邻接像素索引。

1.4 迭代求解能量最小值

根据上述得到的与道路局部区域分割等价的能量函数模型，利用最小割、最大流算法[16]迭代求解获取较优的解。首先，设置初始的前景、背景局部区域，标记部分像素隶属类别α，估计更为准确的道路前景和背景颜色模型θ，此步骤将降低整体能量E； 接着，依据获取的前景、背景颜色模型θ，再次更新影像中的各像素点类别α，此步骤同样将降低整体能量E。通过不断循环迭代，调整颜色模型θ与像素隶属类别α，最终获得较优的能量函数，将道路区域与非道路区域分割。整体流程如下：

1)用户给定局部初始前景、背景标记α，前景αn=1，背景αn=0。

2)初始化前景与背景高斯混合模型θ。

3)对所有待标记像素进行高斯混合模型的成分赋值，即

Kn: =argminKnDn(αn,Kn,θ,zn)。

(15)

4)根据影像数据z，更新高斯混合模型参数，即

θ: =argminθU(α,K,θ,z)。

(16)

5)估计分割，即更新隶属标记α，利用最小割算法求解，即

minαminKE(α,K,θ,z)。

(17)

6)重复步骤3)―步骤5)，直到能量收敛。

其中步骤1)—步骤2)为初始化，步骤3)—步骤6)为迭代优化。

2 增量式分割道路提取

2.1 局部有效区域中的道路提取

在高空间分辨率遥感影像中，道路要素的分布与草地、林区和城镇建筑具有明显区别。道路通常贯穿于整个遥感影像，并在影像中仅占有部分区域，其功能特征直接导致道路分布十分分散。对于马尔科夫随机场框架下的道路分割，只有在局部进行道路提取才能获得鲁棒的道路模型； 如果在影像全部区域中进行道路模型的训练与估计，则很有可能使道路分割出现明显错误。由于路面仅与邻近道路的背景具有相关性，离道路较远的背景区域对局部道路提取并无益处； 相反，由于高空间分辨率遥感影像中“异物同谱”、“同物异谱”现象的干扰，冗余的背景反而会对马尔科夫随机场模型中的数据项U的估计产生干扰，致使前景模型和背景模型混淆。如图1所示，绿色区域为马尔科夫随机场模型的建立范围，黄色区域为无效区域。增量式分割道路提取仅考虑局部有效区域。因此，本文通过增量式道路分割提取策略对所提取的有效区域进行划分，并从有效区域中提取道路区域。

图1 道路提取有效区域示意图Fig.1 Validate area of road extraction

2.2 增量式道路分割提取

道路相比其他地物具有独特的拓扑特征，道路路面彼此相连，相邻的道路间具有灰度、纹理、几何特征方面的正相关性； 同样，相邻的背景间也具有高度正相关的灰度、纹理特征。在连续道路空间中，利用邻接已知道路与非道路区域，提取道路特征，可有效地对相邻区域的像素类别隶属变量α进行估计。相比其他全局提取道路的方法，该策略更为快速、准确。本文通过对前景目标进行膨胀，不断更新有效区域，然后在有效区域中不断更新道路分割。

道路增量式提取流程如图2所示。

图2 道路增量提取原理示意图Fig.2 Sketch map of road incremental extraction principal

图2中，白色为待分割有效区域TU，蓝色为提取的道路区域TF，红色为提取的非道路区域TB。道路增量式提取流程依据ABCBCBC…的顺序进行。首先，如图A进行道路增量式提取初始化，提供局部初值道路区域TF与背景区域TB； 然后根据道路前景TF进行形态学膨胀运算，获取有效待分割的区域TU。图B中，利用马尔科夫随机场方法对图A中获得的TU中的像素进行学习估计，利用前景TF与背景TB像素估计前景和背景RGB三通道高斯混合模型θ； 然后进行判别，将TU中所有像素划分为道路区域TF与非道路区域TB。至此，相比图A，前景道路区域TF与背景区域TB均有所增长。在此过程中，前景道路的增长对有效区域的扩展起决定性作用。图C中，固定图B中已判别像素的前景和背景类别，直接对相应像素点数据项U赋值，对于所属类别赋值为λ，反之赋值为0； 然后，仅对前景道路TF进行膨胀运算，获得更大的待分割区域，即有效区域TU，再利用图B中的步骤对前景道路进行提取。循环B，C步骤，直到完整道路提取完毕，即前景TF不再增长，有效区域TU不再增长。整体分割流程如图3所示。

图3 道路增量式提取流程Fig.3 Road incremental extraction process

3 实验与分析

本文实验数据为高空间分辨率RGB波段遥感影像和全色遥感影像，道路宽度范围在8～30个像素之间，主要为无过多遮挡的城际道路与乡村道路。

实验参数设置如下： 马尔科夫随机场模型中，采用8邻接像素邻接关系； 对于TU中的像素学习与估计优化的迭代次数为10次； 公式(13)中γ设为10，λ设为9倍γ，高斯混合模型分量个数K为3，道路区域增长利用半径为20个像素的圆形要素进行前景膨胀。

3.1 设置初值

本实验由人工设置初值，初值设置非常简单，仅需要标记局部道路的部分像素与该像素附近的部分非道路区域即可。如图4所示，红色为前景道路标记区，蓝色为背景非道路标记区。若提供部分道路矢量，可利用已有矢量自动设置前景和背景初值。

图4 初值设置Fig.4 Initial setting

3.2 增量式马尔科夫随机场分割

增量式分割流程效果如图5和图6所示。每个场景从初始分割，经过若干次前景膨胀再分割的迭代操作，逐步扩张其分割区域，实现分割区域的不断扩大，从而提取出完整的道路范围。

图5增量式分割提取流程

Fig.5Incrementalsegmentationprocess

(a) 乡村道路

(b) 公路道路

3.3 实验结果评价与分析

在本文的实验中，分别对全色影像公路、多光谱影像公路以及全色影像乡村道路进行提取实验。实验中发现，不论对单条道路还是平行并排道路，本文算法均有效，克服了传统追踪方法只能对单一道路追踪的弱点。同时，该算法能适应道路变化的情况，对复杂道路情况的适应性好。像素级精度评价结果如图7所示。

(a) 全色影像公路提取结果 (b) 多光谱影像公路提取结果 (c) 全色影像乡村道路提取结果

(d) 全色影像公路参考结果 (e) 多光谱影像公路参考结果 (f) 全色影像乡村道路参考结果

(g) 全色影像公路提取结果划分(h) 多光谱影像公路提取结果划分(i) 全色影像乡村道路提取结果划分

图7实验结果评价

Fig.7Evaluationofexperimentalresults

图 7(a)—(c)为本文道路分割结果中的道路区域； 图7(d)—(f)为道路区域覆盖的真实区域； 图7(g)—(i)为道路提取的结果划分。红色表示误检像素，即非道路区域被判别为道路区域的像素； 蓝色表示漏检像素，即未被检测出的道路区域像素； 绿色表示正确检测的像素，即属于道路区域且被正确检测的像素。

采用F-Beta测度[17-18]作为量化评价精度测量指标。F-Beta测度是精确率与召回率的调和平均数，综合考虑提取算法的精确率和召回率以修正真值中存在的误差。F-Beta测度定义为

(18)

式中：Precision为精确率；Recall为召回率；β为参数，本文中取β2=0.3[19-20]。

对全色影像公路、多光谱影像公路与全色影像乡村道路分别另选20段区域进行像素级道路分割精度评估(表1)，以各影像精度指标平均值作为最终统计量。

图8展示了一个更加复杂场景的道路增量式分割结果。依照图3所示流程，逐步对各条道路设置初始前景和背景，并进行分割。值得注意的是，对于大范围的遥感影像，可将其划分为若干小影像区域进行处理，以避免不同区域影像前景和背景存在的较大差异。图8(a)中的红线与蓝线分别为各条道路初始前景和背景设置，绿色区域为分割获得的道路区域； 图8(b)为图8(a)分割出的绿色道路区域对应原始影像的裁剪图。

表1道路提取结果评估

Tab.1Evaluationofroadextractionresult

道路类型精确率召回率F-Beta测度全色影像公路0.920.740.87多光谱影像公路0.840.950.87全色影像乡村道路0.710.900.75

(a) 道路提取结果(b) 道路提取结果裁剪图

图8一个更加复杂场景的高空间分辨率遥感影像道路提取结果

Fig.8Roadextractionresultfromamorecomplexscene

4 结论

1)本文提出了一种基于增量式马尔科夫随机场分割的高空间分辨率遥感影像道路提取方法。利用基于RBG波段的高斯混合模型对道路影像前景和背景建模，在马尔科夫随机场模型框架中进行优化，完成对道路区域的分割； 并利用分割前景区域膨胀的方法增量式确定道路提取有效区域，仅在有效区域内完成道路提取，以剔除对于道路提取冗余无效的干扰区域的影响，实现道路区域的精确分割。

2)实验结果表明，所提方法对于前景和背景具有一定灰度差异地区的道路提取效果显著，并且彩色影像在颜色空间比全色影像具有更高的前景、背景模型区分度和更高的道路提取精度。

3)该算法对于平行道路和变宽道路提取适应性强，能够有效应对城际公路与乡村道路提取； 但对于边界模糊的影像区域，道路边界提取仍存在误差。后续研究可通过向分割模型中引入纹理信息并结合道路的几何结构特征，更精确地提取道路。