一种基于特征点过程的遥感影像道路中心线提取算法

孔令然，刘立志，吴 优，张秋良，魏向祎

(1.中水北方勘测设计研究有限责任公司，天津 300222；2.黑龙江工业学院，黑龙江 鸡西 158100；3.内蒙古农业大学 林学院，内蒙古 呼和浩特 010019；4.辽宁工程技术大学 测绘与地理科学学院，辽宁 阜新 123000)

0 引言

道路是最常见的基础设施之一，因此提取道路信息对交通运输规划、GPS导航监测和地图制图更新等领域有着极其重要的作用[1-3]。而且，随着遥感影像分辨率的提高，道路在影像中的轮廓愈加清晰，道路的提取可为其他地物目标识别和提取提供上下文信息[4-5]，因此道路提取已逐渐成为受关注的课题之一。到目前为止，已有大量用于道路提取的算法被提出和应用，其中具有代表性的算法主要有基于主动轮廓模型的提取算法[6-7]、基于数学形态学的算法[8-9]、基于神经网络的算法[10-11]和基于点过程的算法[12-13]等。基于主动轮廓模型的道路提取算法是一种半自动算法，主要根据人为给定的道路种子点和方向遍历整个影像以获取道路网，虽然可以较准确地获取道路网，但由于其需借助人力获得目标的初始位置和方向，降低了算法的提取效率。用于道路提取的数学形态学算法是一种基于像素的提取算法，利用道路光谱与其他地物光谱的异质性进行二值化处理，再根据形态学边缘检测等后处理获取道路，虽然该算法的提取速度较快，但算法容易受到噪声干扰而使结果的完整性和连接性大大降低。基于神经网络的道路提取算法的核心是从训练数据中大量采样获得训练集，通过训练集检验测试数据提取道路，虽然可以比较精确地提取出道路网，但极其依赖训练集，因此算法泛化性较低并且属于半自动化算法。基于点过程的道路提取算法是目前研究最通用的算法之一，比如基于标值点过程的提取模型研究，其主要思想是利用点过程中的随机点定位道路的位置，再根据道路的局部形态标值点的几何结构，对点过程和几何结构进行连接性约束和光谱约束获得目标提取模型，考虑了道路的几何特性、局部网络结构特性和光谱特性，但模型在一定程度上变得更复杂。由此可见，基于主动轮廓模型的提取算法侧重考虑了道路的局部几何结构而忽略了道路的光谱性，基于数学形态学的算法和基于神经网络的算法侧重考虑了道路的光谱性而忽略了道路的局部几何特性，而基于点过程的方法兼顾了道路的2种特性，更具有鲁棒性，但由于点过程中的点是随机分布的，从而算法求解过程中遍历了非必要的道路点，因此降低了算法的提取效率。因此，为了改善传统基于点过程的道路提取算法提取效率低的问题，本文提出了一种基于特征点过程的遥感影像道路中心线提取算法。

在遥感影像中，道路可看作是由线段连接构成的线状网络结构，并且道路内部具有光谱匀质性，道路内部与邻近两侧地物存在光谱异质性。因此，根据道路光谱性可构建特征点过程，根据道路网络结构特性，可构建线状网络结构约束模型，进而根据道路内部光谱一致性和与两侧地物光谱异质性建立学生-检验光谱测度约束模型，最后结合以上模型构建道路提取模型，并设计RJMCMC模拟算法[14-15]，以后验概率最大为准则[16-17]求解道路提取模型的待解参数集。采用不同的高分辨率遥感影像对提出算法和对比算法进行实验，以缓冲区评价方法进行定性和定量评价。结果表明，提出算法在效率和质量上均优于对比算法，验证了提出算法的有效性。

1 算法描述

对于一幅给定的道路遥感影像，设其所有的像素灰度值为矩阵Z={Zi=Z(xi,yi):i=1,2,…,m, (xi,yi)∈D}，其中i为像素索引，(xi,yi)为定义在图像域D(平面)上的像素格点坐标，m为总像素数，Zi为第i个像素的灰度值。定义在D上的离散随机场为z={zi=z(xi,yi):i=1,2,…,m, (xi,yi)∈D}，其中zi为第i个像素灰度值的随机变量，从空间统计学的观点来看，Z可以看作是离散随机场z的一个实现。

在遥感影像中，可将道路中心线视作由一条条线段连接而成，并且线段所在的道路内部具有光谱匀质性，所在的道路内部与两侧邻域具有光谱异质性。因此，根据道路的光谱特征、几何特征和网络的局部特征建立特征点过程、网络结构约束模型和光谱测度模型。

1.1 特征点过程

道路特征点过程如图1所示。在影像平面上的道路区域可生成n个满足道路光谱特性的随机点构成点集F={(xj,yj):j=1,2,…,n,(xj,yj)∈D}，其中(xj,yj)为点的位置，j为点的索引，设点的个数n服从强度为β的泊松分布(如式(1)所示)，其位置分布在满足道路光谱特性的区域(点之间的距离可根据道路的曲率进行约束)，如式(2)所示：

(1)

(2)

式中，μZ和σZ分别为道路区域的光谱测度均值和标准差。值得注意的是，考虑到道路区域可能含有其他地物(如阴影、车辆等)，为了解决当前具有复杂光谱的道路漏提问题，将含有其他地物的道路体看作是由多种地物共同构成的一个反射体，可利用高斯混合模型进行建模。设某一点(xj,yj)所在的道路体A(xj,yj)为一定长宽比的矩形区域，则该区域中的光谱测度Z(A(xj,yj))满足一定的高斯混合分布：

(3)

式中，O为道路区域中的混合地物个数；o为混合地物索引；ϖo为混合地物o的权值；O1和O2为根据道路区域中的混合地物情况预估的值，O在[O1,O2]满足均匀分布。值得说明的是，尽管道路区域混有其他地物，从整体上看道路区域的光谱仍具有一定的匀质性。综上，随机特征点过程可表示为Φ=(F,n)，其联合分布的概率为p(Φ)=p(F,n)=p(F|n)p(n)。

在遥感影像中，道路结构具有显著的连续性、光滑性和曲折性等，因此为约束提取结构满足以上结构特性，可在特征点过程的基础上建立线结构约束模型。首先，采用最大最小法则连接所有特征点构成初始道路结构，如图2所示。

图2 线网络连接Fig.2 Connection of line network

图2中任意2点(xj,yj)和(xj′,yj′)之间的距离ljj′若满足lmin≤ljj′≤lmax(lmin和lmax分别为给定的最小距离和最大距离)，则这2点可用线连接成道路段，否则不连接。因此初始道路结构可表达为：L={Ljj′:lmin≤ljj′≤lmax,j∈F,j′∈F,j≠j′}。

1.2 网络结构约束模型

由于实际的道路结构基本呈现的是一条条连续起伏的线状结构，因此可设定被2条线连接的特征点的先验概率最大，被一条线连接的特征点的先验概率次之，孤立的特征点的先验概率最小，具体约束如下:

(4)

式中，γ={γj:j∈Φ}代表点连接的线段条数；τ1,τ2,…,τ3为能量值，且满足τ2<τ1≤τ3≤τ0。此外，线与线之间的分布关系可按照道路沿向的特点进行约束，如图3所示。

图3 夹角与中心距离关系Fig.3 Relationship between included angles and center distances

设L,L′∈L是2条不同的线段，根据线段端点可计算出其长度lL,lL′，角度θL,θL′和中点间的距离dLL′，进而判断2条线段之间的走向是否满足道路结构特征。若dLL′≤max{lL,lL′}/2，则2线段为重叠走向，但如果角度差θLL′=min{|θL-θL′|，π-|θL-θL′|}大于走向阈值θ，则2线段的走向可接受，否则就应该被约束，其约束如下：

(5)

式中，θ={θLL′:L，L′∈L}代表线段角度差集合；d={dLL′:L，L′∈L}代表线段间中点距离集合。综上，道路结构约束模型可表示为p(L|Φ)=p(γ)p(θ，d)。

1.3 光谱测度约束模型

为了使提取的道路段与影像中的道路准确吻合，以学生-检验[18-19](Student’s-test)为依据，比较道路段所在区域的光谱匀质性和其与邻域之间的光谱差异性，进而建立符合数据特性的势能函数。首先，在已知影像分辨率的前提下，为验证提取线段符合数据特性，进行如下操作：

图4 道路段缓冲区模型Fig.4 Buffer model of road segment

(6)

(7)

(8)

为了使提取线段的光谱性同时满足内部匀质和内部与邻域异质的特点，将数据势能函数的变量设为式(8)与式(7)的比值TL(数据势能变量)来进行判断和约束，此外令TI(L)的最小值为1以避免线段所在的道路内部超匀质而导致比值过大影响正确判断。

(9)

显然TL∈[0, ∞)，如式(10)所示，通过数据势能函数δL将TL变换到[-1，2]，对该提取线段所处的状态进行奖惩。

(10)

式中，η1和η2为实验所得的经验值，具有鲁棒性，适用于具有相似光谱特性的地物目标提取。

最后，根据Gibbs能量函数的原理[20-21]，设τd为正的权重值(随模拟退火温度变化而变化)，由此得出所有提取线段与数据一致性的似然概率：

(11)

1.4 提取模型及其模拟算法

根据贝叶斯理论[22-23]可知，后验概率正比于先验概率与似然概率的乘积，因此结合几何特征、光谱特征和拓扑特征可计算出后验概率，则获得道路目标提取模型：

p(L,Φ|Z)∝p(Z|L,Φ)p(L,Φ)=p(Z|L)p(L|Φ)p(Φ)=

(12)

由于S={L,Φ}={L,F,n}的维度随着n的改变而改变，因此，为高效求解道路目标提取模型，采用RJMCMC和模拟退火优化算法实现在不同维度的空间中模拟采样，随着温度T不断降低，后验概率不断变大，最终获得全局最优解。RJMCMC是一种关于Green比[24]R的模拟算法，即计算从当前状态S(t)={L(t),F(t),n(t)}(t为迭代次数)跳跃转移至新状态S*={L*,F*,n*}的接受率α(S,S*)，此外，设计高效的转移操作是提高接受率的关键，加快算法收敛于最优分布。

(13)

式中，Q(S→S*)为状态S根据转移操作到状态S*的概率密度函数；Q(S*→S)为状态S*根据转移操作到状态S的概率密度函数；φ为引入向量以确保状态S和状态S*之间转移的维度平衡；p(φ)为向量φ的概率密度函数；|∂(S*)/∂(S(t),φ)|为雅克比行列式。可以看出，设计高效的转移操作是提高接受率的关键，从而加快算法收敛于最优分布。

1.4.1 特征点和道路段的简单变换转移操作

图5 特征点的简单变换Fig.5 Simple changes of feature point

αFM(S,S*)=min{1,RFM}=

(14)

式中，QFM(S→S*)为状态S根据特征点的变换转移操作到状态S*的概率密度函数；QFM*(S*→S)为该变换操作的反向操作的概率密度函数。

图6 道路段的位移、伸缩和旋转简单变换Fig.6 Move, expansion and rotation changes of road segment

(15)

式中，QLM(S→S*)为状态S根据道路段的某一简单变换转移操作到状态S*的概率密度函数；QLM*(S*→S)为该变换操作的反向操作的概率密度函数。

1.4.2 特征点和道路段的生灭转移操作

(1)特征点的生灭转移操作。在当前的状态S(t)中，为了使新生成的特征点在道路体内，根据道路的光谱特征可生成相关点。由于道路内部的光谱具有一致性，因此，特征点的生灭操作如图7所示，一个满足道路光谱特征的点可以按照以下步骤生成：

图7 特征点的生灭操作Fig.7 Birth and death of feature point

① 在影像域中找到满足道路光谱特征的点(x*,y*)，即满足，zmin

② 计算点(x*,y*)与邻近点的距离，若距离大于等于lmin且小于等于lmax，则将其与邻近点连接生成新的道路段L*，否则不连接；

③ 计算接受率。

此时，新状态S*={{L(t),L*}, {F(t), (x*,y*)},n(t)+1}，其接受率可计算如下：

(16)

(2)道路段的生灭转移操作。为了提高道路网络结构的连接性，如图8所示，在当前结构S中随机选择一点(xj,yj)，若存在另外的点(xj′,yj′)与其距离满足大于等于lmin且小于等于lmax，则将其连接生成新的道路段L*；否则以(xj,yj)为端点，以θ*∈[0, π]为方向，以l*∈[lmin,lmax]为长度，生成新的道路段L*和特征点(端点)(x*,y*)。

图8 道路段的生灭操作Fig.8 Birth and death of road segment

此时，新状态S*={{L(t),L*}, {F(t), (x*,y*)},n(t)+1}或S*={{L(t),L*}, {F(t),n(t)}}，其接受率可计算为：

(17)

当转移操作的接受率α低于1时，则在[0, 1]随机产生一个数，若α高于，则接受该操作，保留结构S*；若接受率低于，则拒绝该操作，保留结构S。同时使用合适的转移核，充分利用各自的优越性，以加快算法收敛于最优分布。

2 实验与结果

所有实验均在同一CPU配置下执行，表1列出了实验数据，其中，Ikonos卫星影像的分辨率为1 m，被选取的道路场景位于德国科隆市，主要有乡村道路和城镇道路；GeoEye-1卫星影像的分辨率为0.5 m，被选取的道路场景位于日本仙台市，主要为乡村道路；Gaofen-2卫星影像的分辨率为0.8 m，被选取的道路场景位于中国辽宁省沈阳市区，其中的道路区域的光谱特征复杂，含有其他非道路地物。

表1 实验数据Tab.1 Experimental data

为了有效评估提出算法的优越性，采用缓冲区评价法[13]对提取结果进行定量分析。缓冲区评价法的主要思想是以提取的道路中心线为基准向两侧做缓冲区域(像素格点)，则提取线与参考线的像素重合率可定量评价提取结果的精度，参考线所在的像素格点落入缓冲区域的个数百分比可定量评价提取结果的完整性，如表2所示(取第0，2，中间和边缘层作为评价结果，分别表示为表中的ρ0，ρ1，ρ2和ρ3)。

表2 提取结果定量评价Tab.2 Quantitative evaluation of extraction result

首先，利用Ikonos卫星影像对提出算法进行可行性验证。局部道路的提取结果如图9所示，其中，图9(a)～图9(d)分别是原始影像，大小均为256 pixel×256 pixel，分辨率为1 m，覆盖的区域中具有四岔路、三岔路、弯曲道路和阴影遮挡的道路等；图9(e)～图9(h)分别是原始影像对应的提取结果。从结果中可以看到，提出算法能有效地提取出交叉道路和弯曲道路，并且可提取出树冠阴影遮挡的道路段，一定程度上说明了提出方法的可行性和有效性。此外，如表2的缓冲区评价结果所示，4组提取结果与标准提取线(由多位专家老师绘制得到)的重叠率分别可达54%，52%，55%和57%，提取的完整率可达99%以上，提取的时间都在2 min以内，进一步说明了提出方法的有效性。

图9 局部道路提取结果Fig.9 Extraction result of local road

其次，为了进一步说明提出算法的优越性，采用文献[13]中的算法作为对比算法，对2种来自不同数据源的遥感影像进行实验。2种算法的提取结果如图10所示。其中，图10(a)和图10(b)分别为Ikonos卫星影像和GeoEye-1卫星影像(800 pixel×800 pixel，分辨率：0.5 m)，图10(a1)和图10(b1)分别为标准提取线。如图10(a2)和图10(b2)所示，对比算法的提取结果不仅在道路拐弯转折处、阴影处不连贯(如图中红色框区域)，而且在交叉口处形成了缝隙(如图中绿色框区域)，这是由于文献[13]算法中的局部结构模型仅仅考虑了2条道路段的端点连接，并未充分考虑多条线段的端点连接；然而，如图10(a3)和图10(b3)所示，提出算法在很大程度上改善了道路提取的连通性和完整性。此外，提取结果定量评价如表3所示。由表3可以看出，提出算法的提取结果精度和模拟运行速率远远高于对比算法，这是由于对比算法在对道路位置分布建模时，认为道路位置分布为随机分布，而提出算法充分考虑了道路的光谱特性，建立特征点过程可在一定程度上正确聚类道路的位置分布和分布趋势，从而建立一个完备的网络结构模型，使得算法的模拟运行速率也大大降低。2种算法平均能量收敛的对比如图11所示。可以看出，从初始状态到最后收敛状态，基于特征点过程的道路提取算法更加平稳，说明特征点过程在最初状态就能确定大量道路的分布位置，将提高求解速率。显然，对比算法在一点无法快速获得道路的准确位置，使得整个过程缓慢。综上，提出算法在精度和连续性等方面有了大大提高，证明了算法的有效性。

图10 提出算法与文献[13]算法的结果比较Fig.10 Comparison between results from reference [13] algorithm and the proposed algorithm

表3 提取结果定量评价Tab.3 Quantitative evaluation of extraction result

图11 平均能量收敛对比Fig.11 Comparison of average energy convergences

较大尺度的Ikonos卫星彩色合成影像(大小：1 951 pixel×1 207 pixel，分辨率：1 m)的道路提取结果如图12所示。从提取结果中可看到，提出算法可检测到并提取出肉眼可见的全部道路；从提取的细节图(如图中蓝色和绿色方框所示)中可以看到，提出方法可有效地提取出岔路结构和被大面积树冠遮挡的道路结构。根据缓冲区评价结果(如表4所示)可知，提取的完整性可达98%，充分说明了提出算法的有效性。值得一提的是，所有实验采用的模拟退火温度下降速率皆为0.999，且实验结果的质量较好，说明提出算法对温度下降速率不敏感，在同等温度下降速率的条件下，皆可提取全局最优的道路。

图12 Ikonos卫星影像提取结果Fig.12 Extraction result from Ikonos satellite image

表4 提取结果定量评价Tab.4 Quantitative evaluation of extraction result

最后，为了说明提出算法可适用于具有复杂光谱特征的道路提取，采用了Gaofen-2卫星影像进行了实验，该影像的的分辨率为0.8 m，实验场景大小为1 299 pixel× 1 104 pixel，场景中的道路区域存在2种情形，一是道路区域的光谱特征极其匀质(见黄色虚线框区域)，二是道路区域含有其他大量非道路地物，主要为车辆和交通标线(如绿色虚线框区域)，针对这一道路场景，将道路区域中混合地物个数设置在[O1,O2] = [1, 3]，进而可进行实验。Gaofen-2卫星影像的实验结果如图13所示，可看到提取结果的连通性较高，能够有效地提取出以上2种情形的道路。图13(d)为提取结果与参考线及其缓冲区(根据道路的平均宽度和分辨率设定，本文设定为24个pixel)的叠加结果，可以看出，所有提取结果均落入了缓冲区，说明提出算法能够有效地实现具有复杂光谱的道路提取。提取结果定量评价如表5所示，由缓冲区评价结果可知，提取的精度可达96.88%，说明提取结果的完整性较高。另外，提出算法难以提取出不清晰且宽度极窄的道路结构。

(a)原始数据

(b)参考线

(c)提取结果

(d)叠加结果图13 Gaofen-2卫星影像提取结果Fig.13 Extraction result from Gaofen-2 satellite image

表5 提取结果定量评价Tab.5 Quantitative evaluation of extraction result

综上所述，提出算法有效地解决了传统基于标识点过程道路提取方法存在的完整性差、速率低、交叉口提取不精准等问题，可有效地应用到具有不同光谱特征的道路提取，且提出算法适用于不同数据源的遥感数据，具有一定的普适性和可靠性。

3 结束语

为了提取出连续、完整且精确的道路网络结构，提出了一种改进标识点过程的道路提取方法。利用道路的光谱特性建立了随机特征点，进而根据最大最小距离建立一个网络结构模型。根据道路的结构连续性、道路内部的光谱匀质性和道路内部与邻近两侧地物的光谱异质性可建立道路提取模型。通过道路特征点建模网络结构解决了传统标识点过程模型的随机性，进而提高了提取效率；对网络结构进行连续性约束改善了传统局部结构约束模型，可提高整体网络的连接性。设计的RJMCMC模拟算法实现了算法的全局最优收敛，其中，特征点和道路段的生灭转移核提高了道路网络提取的连接性，也在一定程度上提高了算法求解的效率。显然，提出算法通过道路的光谱性建立特征点过程算法依赖于道路的光谱特征且算法的参数较多，提取效率仍然难以满足大规模的城市道路提取，未来将结合道路的网络拓扑特征和光谱特征建立提取算法，以实现城镇道路提取。