1 引言

20世纪70年代，有学者开始研究利用计算机提取遥感影像中的道路，并且取得了一定的进展。文献^[1]基于局部灰度一致性并结合形状特征提取道路；文献^[2]给出了一种从高分辨率遥感影像中自动提取无清晰边缘线道路的方法；文献^[3]结合LiDAR数据和遥感影像两种信息提取城市道路并进行建模。此外，还有分别基于平行线对^{[4, 5]}、基于二值化知识^[6]、基于窗口模型特征^[7]、基于模糊聚类^[8]、基于灰度形态学^{[9, 21]}等理论的道路提取方法。

提取道路时，多数方法均会得到断裂的道路段这一中间结果，如何连接断裂的道路段是提取过程中的难点，而已有的研究大多针对其余阶段展开。对此，本文探索了一种连接初始道路段的策略：首先构建马尔可夫随机场下的经典能量函数，为了连接道路时能获得无粘连块、无孔洞的道路段，在经典能量函数中分别加入两种形状先验，构成新的能量函数，利用图割解算新能量函数后即可达到连接初始道路段的目的。

在连接初始道路段方法的基础上，给出了一种从高分辨率影像中提取道路段的方法：首先利用Gabor滤波器进行预处理操作，获得待分割特征矢量；基于图割理论解算经典能量函数，从待分割特征矢量中检测道路候选点；接着应用直线匹配法获取初始道路段；最后利用图割解算融合两种形状先验信息的能量函数，结合动态外推的思路连接初始道路段，获取完整道路网。

2 算法具体步骤

2.1 影像预处理

仅依靠光谱信息难以区分遥感影像中的道路与其他地物，利用Gabor滤波器获取影像的角度纹理特征^[10]，并与灰度特征共同组成待分割特征矢量。

2.2 获取道路候选点

2.2.1 经典能量函数

本文基于图割理论解算经典能量函数，从而从待分割特征矢量中检测出道路候选点。本文使用的经典能量函数为

点击浏览公式

式中，等号右边第1项为数据项，第2项为边缘项。为了利用图割解算经典能量函数，必须在2.1节特征矢量的基础上构造满足式(1)的图。

2.2.2图的构造

构造一幅图的实质是分别构造出N链与T链，N链与T链分别对应式(1)的边缘项、数据项。

本文利用 Boykov 在文献^[11]中的公式构造N链

点击浏览公式

大多数的文献一般都利用高斯密度函数刻画T链，高斯密度函数包括单高斯模型和高斯混合模型 (Gaussian mixture model，GMM)^{[11, 12, 13, 14]}。单高斯模型适合于简单数据的建模，高分辨率遥感影像中存在道路、房屋、树木、河流等多种地物，并且相同的地物可能会呈现出不同的形态，因此本文利用GMM 构造T链。

高斯混合模型的概率密度函数如公式(3)所示

点击浏览公式

式中

点击浏览公式

本文计算T链的具体步骤为：

(1) 在特征矢量图上选择若干有代表性的目标样本和背景样本。

(2) 目标区域高斯混合模型表示为p_Fx=∑Kk=1α_kNx|μ_k,∑_k，确定K值。利用 EM 算法估计模型中的所有参数，得到p_F(x)。同理得背景区域模型的概率密度函数p_B(x)。

(3) 将特征矢量值分别代入p_F(x)和p_B(x)，得到每一像素点分别属于目标模型和背景模型的概率密度值，分别表示为P_F和P_B。

(4) 令(P,S)=- log (P_B)，(P,T)= -log (P_F)，Value _min 为 (P,S) 、(P,T) 中的最小值，Source、Sink分别表示像素与源、汇相连接的T链值，则Source=(P,S)-Value _min ，Sink=(P,T)-Value _min 。

高斯混合模型中唯一需要确定的参数为K，一般通过目视观察的方法确定。一般情况下，K值越小，运行时间越短，但是精度将受到影响；K值越大，精度将得到保证，但运行速度可能会变慢。需要结合影像确定最佳K值。

选择一种最大流最小割算法解算图后，即从待分割特征矢量中检测出道路候选点。本文使用 Boykov 的基于增广路径的最大流算法^[12]获取最小割。

2.3 获取初始道路段

利用文献^[15]提出的直线匹配法处理道路候选点以得到初始道路段。

2.4 连接初始道路段

在2.3节中获得的初始道路段存在着断裂、缝隙和孔洞等问题，需要进一步处理。本文主要在式(1)基础上分别引入两种形状先验能量项，得到新的能量函数；接着基于图割理论解算新的能量函数，并结合动态外推的思路连接断裂的道路、填充道路孔洞和缝隙。

将物体的形状先验知识与图像分割相结合的思路，已经在水平集和一些其他的曲线演化方法中得到了较为广泛和深入的研究。文献[13,14,16—20]都在分割图像的过程中应用了物体的形状先验。

2.4.1 新的能量函数

为了保证在连接初始道路段时，仅得到待连接道路、排除其他道路及干扰物的影响，首先引入一种用无符号距离函数表达的形状先验。为无符号距离函数，(x)表示点x与轮廓间的最小欧氏距离。由于一条道路段可以用多个大小、方向各异的小矩形近似表示。因此，本文构造矩形形状的模板。用E_s1表示由无符号距离函数获得的形状模板的能量项

点击浏览公式

将式(5)与式(1)相结合，得到一个新的能量函数

点击浏览公式

式中，λ₂表示E_i与E_s1间的权重，取值在^{[0, 1]}之间。为了利用图割理论最小化式(6)表示的能量函数，同样需要根据公式(6)分别构建N链和T链。将E_s1表示为N链，并与式(1)的N链合并，公式(6)的T链则与公式(1)一致，从而得到公式(6)的图，利用 Boykov 算法解算即可。

图 1中的第1组试验影像中，道路与背景的光谱值差别较为明显，道路周围存在部分与道路相似的干扰物。首先根据公式(1)分割试验影像，获取了完整的道路，但同时也分割出了部分干扰物(如图(b)所示)。添加图(a)中的矩形模板后，再基于公式(6)分割影像。在矩形模板的作用下，道路保持完整，同时剔除了图(b)中的粘连块(如图(c)所示)。

图 1匹配的矩形模板在道路分割中的作用(第1列：试验影像 第2列：无矩形模板的分割结果 第3列：匹配的矩形模板的分割结果) Fig. 1 Influences of a matching rectangular template in segmenting roads

图选项

第2组试验影像中，道路与背景相似，不易区分。未添加矩形模板的分割结果如图(e)所示，分割结果中包括两条主干道和大量干扰物。解算含有如图(d)所示的矩形模板的能量函数后的结果如图(f)所示。可以看到，已经排除了矩形模板外的绝大部分道路与干扰物。但由于道路周围的屋顶与道路的灰度值非常相似，并且两者距离很近，因此即使添加了矩形模板，道路周围依然存在干扰物(图(f)中的椭圆处)，这是需要解决的一个问题。

图 2显示了分别经过旋转和平移的矩形模板对分割结果的影响。可看出，结合经旋转或者平移后的矩形模板，能够完整地提取出道路，同时去除了道路周围的干扰物。而由于平移、旋转后的矩形模板轮廓上像素的N链发生不可知的变化，导致分割出的道路内部出现了孔洞(椭圆所示)。这是需要避免的另一个问题。为了解决仅利用矩形模板出现的上述两个问题，下面将在公式(6)的基础上进一步引入另一种形状先验信息—星形形状先验^[16]。

图 2 旋转、平移后的矩形模板对分割结果的影响 Fig. 2 Influences of a rotated or translated rectangular template in segmenting roads

图选项

将星形形状作为图像分割中的先验信息，能保证仅分割出一个满足星形形状的物体，而满足星形形状意味着一个不含粘连块、无孔洞的连通成分，这正解决了上文出现的两种问题。

星形形状先验能量项表示为

点击浏览公式

将公式(7)与公式(6)结合，得到最终的能量函数

点击浏览公式

构造满足式(8)的图时，只需将式(7)表示为N链，并与式(6)的N链合并即可。

2.4.2 基于动态外推的思路连接初始道路段

在式(8)基础上，基于动态外推的思路连接初始道路段，图 5为连接原理图。连接初始道路段的步骤包括：

步骤1标记初始道路段。2.3节中，直线匹配法处理道路候选点后同时得到了方向图，方向图表示了每点所属的道路方向。根据这些道路点的方向值为道路段编号，得到一个标记数组。

步骤2计算道路段相关初始值。根据标记数组选择一条初始道路段，计算得到初始道路段主方向、道路段长、宽值及重心位置，这些值主要用于构建矩形模板与星形形状先验。

步骤3动态外推连接初始道路段。在步骤2选择的初始道路段基础上，围绕道路段两端中心点构建分割候选区，在候选区内构建满足能量函数式(8)的图，解算图后即获得新生长的道路段。重复本步骤，直到满足终止条件。图 3是道路连接原理的示意图。步骤3的终止条件包括：每次得到的新生长的道路段的像素个数小于阈值thre1，则退出运算；一条道路段经外推已经与另一条道路相连时，退出运算。

重复上述3个步骤，直至处理完所有的初始道路段，此时也已完成连接。

图 3连接原理示意图 Fig. 3 Diagram of connection

图选项

3 试 验

3.1 单条道路段动态分割试验

为了验证本文的连接方法是否可行，首先动态外推分割一条道路，如果能得到完整的道路段，说明思路合理且适用。通过人工给定相关初始值后解算公式(8)，以得到初始道路段。

试验1影像大小为1600像素×261像素，影像中道路笔直，与周围背景容易区分，初始道路段主方向86 ° ，宽24像素，模板长130像素，thre1为150像素，λ₁=5，λ₂=0.55。动态外推19次后，道路端中心点超出图像范围而退出运算。从图 4( f )可看出，已经较为准确和完整地分割出道路。

图 4 动态分割试验1 Fig. 4 Example one of dynamic segmentation

图选项

图 5 动态分割试验2 Fig. 5 Example two of dynamic segmentation

图选项

试验2中影像大小为1600像素×209像素，影像中道路呈线状，与周围建筑物光谱值近似，并且受到汽车、建筑物阴影的干扰。初始道路段主方向92 ° ，宽30像素，模板长设为100像素，thre1为150像素，λ₁=0.5，λ₂=0.55。动态外推11次后，受到建筑物阴影的影响，新生长道路段像素个数小于阈值而结束运算。

两组动态外推分割试验均得到了较为完整、没有干扰物影响的道路段，可知本文的连接思路可行、有效。

3.2 道路段提取试验

选取两幅多波段遥感影像进行道路段提取试验。图 6(a)显示了两幅不同地区的多波段遥感影像，大小均为500像素×500像素，图 6(b)为图割解算经典能量函数后获得的道路候选点，图 6(c)为经过直线匹配后获取的初始道路段图，去除了杂点、散点，但道路段明显不完整，存在多处断裂及孔洞，图 6(d)为通过解算式(8)并结合动态外推思路获得的连接结果，道路段断裂问题得到明显改善(椭圆处)。

图 6 道路段提取试验 Fig. 6 Experiments of road extraction

图选项

图 7为利用文献^[15]提取的道路段，该方法中，为了消除非道路干扰，需加严道路认定准则，由此会将部分道路误认为非道路，导致道路段不完整(图 7(a)椭圆处)；若为了保持道路段的完整性而放宽道路认定准则，则导致大量的粘连块出现(图 7(b)椭圆处)。由于本文在后处理中运用了两种形状先验，在保持道路段完整性的同时较好地去除了粘连块。

利用下面的公式进行定量评价

准确率=正确提取的道路段面积/人工解译的道路段面积

遗漏率=遗漏提取的道路段面积/人工解译的道路段面积

根据上述评价指标，对本文得到的道路段做了如下统计。

图 7 文献[15]的提取结果 Fig. 7 Extraction results of reference[15]

图选项

从表 1中可看出，提取出的道路段平均准确率可达88.49%，较为理想。

4 结 论

本文给出的从高分辨率遥感影像中提取道路段的方法的创新点在于将两种形状先验能量项引入经典能量函数，并结合动态外推的思路，得到一种连接初始道路段的新方法。

本文未验证该方法在大幅影像中的性能，没有考虑立交路、环形路等形状特别复杂的道路，这都是需要进一步研究的内容。