道路是中国重要的基础交通设施，国家在道路建设上不断加大投资力度，呈现出逐年增长状态。对在建道路项目的实时监测，以保证项目顺利进行，对国家交通管理、经济战略规划都具有十分重要的意义。随着中国高分系列卫星的发射，高分辨率遥感影像的获取更加便捷，高分辨率遥感影像数据包含丰富的光谱、纹理、空间结构特征，逐渐成为遥感技术应用的主要数据来源，为地物提取提供了很好的数据源。遥感具有观测范围大、实时性强的优点，利用高分辨率遥感影像提取规划区内的道路，检测道路变化，为监测道路项目的施工进度提供依据。

目前遥感在道路提取方面，特别是矢量数据和遥感影像相结合的道路提取研究比较多，但用于道路施工监测的方案较少。董明等^[1]用缓冲区方法减少道路搜索区域，用多尺度模板匹配法检测道路，用LSB-Snake (least square B-spline snake)模型进行半自动检测新增道路。张剑清等^[2]通过道路矢量信息获取初始点，并利用均值漂移算法迭代收敛到道路中心点，实现了小比例尺航空影像上道路的自动提取。刘朋飞^[3]在矢量道路附近通过二值模板匹配获取初始点，结合水平集方法迭代至道路轮廓；王双^[4]利用矢量数据构造的B样条曲线作为道路的初始中心线，并采用改进的分级演化方法使模型收敛到正确的位置，得到道路的两条边线；傅嘉政等^[5]使用小波变换进行图像降噪和道路边缘检测, 将图像分为道路区域和非道路区域，利用霍夫变换提取分割道路；丁磊等^[6]从矢量数据获取候选种子点，并通过提炼同质区域的形状特征剔除错误候选点，自动获取负样本点以进行朴素贝叶斯分类，采用邻域质心投票算法从分类影像提取道路中心线。

本文以规划矢量和高分辨率遥感影像为数据基础，提出一种基于图割算法的道路提取方法，监测道路施工进度。

1 研究方法 1.1 图像预处理

在高分辨率影像中，道路不仅仅是简单的带状地物，还包含丰富的细节信息，如车辆、标志线和绿化带等，为道路提取带来了一些困难。双边滤波对图像进行预处理，能有效地平滑道路上的标志线等噪声，同时又较好地保持道路边缘信息。

双边滤波是一种非线性滤波方法，将高斯滤波与邻域滤波结合，不仅考虑了中心像素点与其邻域像素点之间的灰度值的相似程度，而且考虑了两者之间几何距离，使得双边滤波在有效地平滑噪声的同时较好地保持图像边缘信息^[7]。双边滤波由式(1)定义为：

$ I^{\prime}(x, y)=\frac{1}{w_{p}} \sum\limits_{i, j \in \mathit{\Omega }} w_{d}(i, j) w_{r}(i, j) I(i, j) $

(1)

式中，I′(x, y)为位置(x, y)滤波后像素值；I(i, j)为噪声图像I在位置(i, j)上的像素值；Ω为像素(x, y)的空间邻域像素集合；w_d(i, j)度量了邻域内像素点(i, j)与中心像素(x, y)的空间邻近度；w_r(i, j)度量了邻域内像素(i, j)与中心像素(x, y)的光度相似性；w_p为归一化系数。定义分别为：

$ \left\{\begin{array}{l} w_{d}(i, j)=\exp \left(-\frac{(x-i)^{2}+(y-j)^{2}}{2 \sigma_{d}^{2}}\right) \\ w_{r}(i, j)=\exp \left(-\frac{\|I(x, y)-I(i, j)\|^{2}}{2 \sigma_{r}^{2}}\right) \\ w_{p}=\sum\limits_{i, j \in \mathit{\Omega }} w_{d}(i, j) w_{r}(i, j) \end{array}\right. $

(2)

当像素点与中心点的欧几里德距离增大时，w_d(i, j)会随之减小；而当像素点亮度值与中心点亮度值的差值增大时，w_r(i, j)会随之减小。图像在变化平缓的区域，邻域内像素亮度较为接近，双边滤波即为高斯低通滤波器；图像在变化剧烈的区域，滤波器则用边缘点附近亮度值相近的像素点亮度的平均值替代原来像素点的亮度值。

双边滤波算法主要在于选择合理的空间标准差σ_d和灰度标准差σ_r，根据遥感影像的分辨率和道路的特征，经过实验发现，1 m分辨率的影像，选取经验值σ_d=75，σ_r=25，效果较好，如图 1所示。

1.2 图像缓冲区建立

为了减少影像上其他地物的干扰，同时减少冗余计算，加快算法的处理速度，可以先将不感兴趣区域剔除，突出感兴趣区域。将规划矢量数据和遥感影像配准，矢量叠加影像如图 2(a)所示。根据道路规划的修建宽度，将线状道路中心线向外扩展一定距离的区域建立缓冲区。本文建立了两级缓冲区，一级缓冲区需包含所有道路区域，考虑到配准精度和施工偏差，本文采用的道路规划修建宽度的2倍；二级缓冲区是影像需要处理的部分，本文采用道路规划修建宽度的4倍，如图 2(b)所示，浅色标识为一级缓冲区，浅色区域与深色区域合起来为二级缓冲区。

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

1.3 获取道路区域

1) 图割概述。基于图割的交互式图像分割技术可应用在灰度图像分割问题上。图割使用的图比普通的无向图多了两个顶点，分别用符号S和T表示，如图 3所示。普通顶点之间的边，称为N链，普通顶点和终端顶点之间的边，称为T链。N链的权值反映了相邻像素点的特征差异，T链的权值反映了像素与背景、前景的相似程度。

设X表示图像, X= {x₁, …, x_n, …, x_N}，x_n为第n个像素的特征向量，N为图像像素总数。假设图像的分割S={s₁, …, s_n, …, s_N}，s_n∈{0, 1}，0为背景，1为前景，其能量函数可表示为以下等式：

$ E(S, X)=U(S, X)+V(S, X) $

(3)

式中, U(·)是区域项，表示像素被归类为目标或者背景的惩罚；V(·)是边界项，表示邻域像素之间不连续的惩罚。图割的目标就是把图的顶点划分为两个不相交的子集，此时能量函数的值达到最小。

2) 图的构造。为了利用Graph Cuts解算经典能量函数，从影像中分割出道路区域，必须在图像上构造满足式(3)的图，图的N链与T链分别对应式(3)中边界项、区域项。N链的构造公式表示为：

$ V(S, X)=\sum\limits_{(m, n) \in C} \frac{1}{\operatorname{dis}(n, m)} \exp \left(-\left\|\mathit{\boldsymbol{x}}_{n}-\mathit{\boldsymbol{x}}_{m}\right\|^{2}\right) $

(4)

式中，C表示相邻像素的集合; dis(·)是相邻像素的欧氏距离。当两个相邻像素特征值差别‖ x_n－ x_m‖很小，那么它属于同一个目标或者同一背景的概率就很大，能量越大；如果他们的差别很大，说明这两个像素很有可能处于目标和背景的边缘部分，则被分割开的概率比较大，能量越小。

本文利用高斯混合模型(Gaussian mixture model, GMM)构造N链，对目标和背景分别使用K个高斯分量的全协方差GMM来进行建模，引入向量k =(k₁, …, k_n, …, k_N)，其中k_n是第n个像素对应于的高斯分量，k_n∈{1, …, K}。每个像素只能属于目标GMM的某个高斯分量，或者背景GMM的某个高斯分量。高斯混合模型的概率密度函数为：

$ P(\boldsymbol{x})=\sum\limits_{k=1}^{K} \pi_{k} \cdot N\left(\boldsymbol{x} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right) $

(5)

式中，π_k表示每一个高斯分量的权重；N(·)表示高斯概率分布函数为：

$ \begin{array}{c} N\left(\boldsymbol{x} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)= \\ \frac{1}{\sqrt{(2 \pi)^{d}\left|\boldsymbol{\Sigma}_{k}\right|}} \exp \left[-\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{\mu}_{k}\right)^{\mathrm{T}} \boldsymbol{\Sigma}_{k}^{-1}\left(\boldsymbol{x}-\boldsymbol{\mu}_{k}\right)\right] \end{array} $

(6)

式中，μ_k和Σ_k分别表示高斯混合模型第k个高斯分量的均值向量和协方差矩阵。

T链的构造公式表示为：

$ U(S, X)=\sum\limits_{n=1}^{N}-\ln P\left(\boldsymbol{x}_{n} | s_{n}, \theta\right) $

(7)

式中，θ={π(s_n, k_n), μ (s_n, k_n), Σ (s_n, k_n)}，表示前景和背景GMM的参数，包括高斯分量的权重、均值向量和协方差矩阵。

计算T链的具体步骤为：①道路一级缓冲区作为“可能包含道路目标的区域”，初始化标签α=1。二级缓冲区减去一级缓冲区(图 2(b)中的深色区域)作为背景区域，初始化背景像素的标签α=0。②将α=1的像素作为前景像素，通过k-means算法来初始化前景GMM，使用α=0的背景像素，通过k-means算法来初始化前景GMM。③对待处理图像(二级缓冲区)中的每一个像素，通过初始化的标签来判断该像素属于背景像素还是前景像素，然后分配GMM中的高斯分量。④从待处理的图像数据X中，学习GMM的参数，得到前景高斯混合模型概率函数P_F(x)，背景高斯混合模型概率函数P_B(x)。⑤将“可能包含道路目标的区域”中的像素的特征向量带入前景和背景的高斯混合模型概率函数，得到每一像素点分别属于目标模型和背景模型的概率P_F和P_B。⑥f_S、t_T分别表示像素与S点和T点相连接的T链的权值, 则“可能包含道路目标的区域”中的像素的T链的权值为：f_S=－ln(P_F), t_T=－ln(P_B)；背景区域像素的T链的权值为：f_S=0, t_T=max，max为人工设置的极大值。

高斯混合模型中需要确定参数K。研究发现遥感影像每一类地物的光谱特征空间用3个分量左右的高斯混合模型，拟合精度较高^[8-9]。由于背景一般包含植被、裸地，最少两种地物类型，前景一般也有两种以上的地物，但高斯分量多容易过拟合，同时增加算法的计算量，经过试验得到，本文方案的前景和背景高斯模型分量选取5较为合适。

3) 能量最小化。通过构造图，能量函数最小化转化为求解图的最小割问题，网路的最大流等于最小割。因此，针对构造的图，使用最大流算法^[10]求解图的最小割，其对应的能量函数的值达到最小, 得到图分割结果，提取道路区域，如图 4(a)所示。

1.4 道路施工进度

道路作为一种人工地物，路面的主要材质为沥青或者水泥。路面已经铺好的道路，才算是已修建完成的道路。§1.3节中提取的道路区域是已经具备道路形状的区域，包含已修建完成的道路和正在施工的道路。首先，对§1.3节中道路提取结果，进行腐蚀操作，去除独立小区域和道路两边缓冲带，并使用中值滤波去除路面噪声；然后，使用对种子点不敏感的k-mean++聚类算法^[11]将提取出的道路聚为两类，同时得到两个类别的聚类中心。由于水泥在高分辨率遥感影像中亮度最高，土壤稍低，沥青最低。对比两个聚类中心，如果道路规划修建成沥青路面，则聚类中心亮度低的类别为沥青，即已修建道路，另一个聚类中心对应为正在修建道路。如果规划道路为水泥路面，则聚类中心亮度高的类别为水泥，即已修建道路，另一个类别对应为正在修建道路。将分类图合并独立小区域后，得到的道路分类结果如图 4(b)所示。

将道路分类结果与规划适量叠加分析，就可以得到规划图中已修建道路、正在修建道路和未修建道路区域，得到道路施工进度图，如图 4(c)所示。

图 4中，绿色为已修建好道路；红色为正在修建道路；黑色为未修建道路。

2 实验结果与分析

本文道路施工监测方案的具体技术方案如图 5所示。首先，将矢量与高分遥感影像配准，建立道路规划缓冲区；然后，使用双边滤波对遥感图像的进行预处理，在保持边缘信息的同时减少图像噪声，增强道路信息；第三，结合缓冲区，使用高斯混合模型对前景和背景建模，构建能量方程；第四，使用最小割算法获得能量最小化，对道路目标进行自动图像分割；第五，将分割结果进行形态学处理，去除孤立小区域后，使用k-mean++算法将道路分为正在修建的道路和已修建完成的道路，得到道路分类图；最后，和规划矢量图叠置分析得到已建道路、在建道路和未建道路，形成道路施工进度结果。

实验区域选取湖北省武汉市某开发区的道路项目施工区域。实验数据选用同一地区2景不同时间的1 m分辨率高分遥感影像，影像大小1 209行，1 872列，占地面积226.3万m²，规划修建道路11.74 km。

选取道路施工阶段的两个时相的影像进行实验，验证本文算法。时相一的影像拍摄于2012年9月14日，其中有部分道路已经修建完成，有部分道路正在修建中，也有还未开始修建的道路。时相二的影像拍摄于2013年7月12日，计划修建的道路基本已经修建完成，只有很小一段没有铺上路面材质。实验结果如图 6所示，从实验结果可以看出，整体提取效果较好。

为了衡量本文方法的提取精度，将实验结果和人工提取结果经过叠置分析得到误差矩阵，如表 1所示。

由表 1计算得到，时相一的总体准确率为95.3%，已修建道路的准确率为97.4%;时相二的总体准确率为92.1%，已修建道路的准确率为94.6%。不管是从实验结果的目视效果还是统计数据，本文方法都表现出一定的有效性。

3 结束语

为解决使用遥感技术动态监测道路施工进度，提出了一种规划区内道路自动提取识别方案。充分利用规划矢量数据的特征，达到自动化提取的效果。实验结果表明，利用该方法处理已修建道路和正在修建道路特征明显的影像，能够得到较理想的提取精度。当道路形状不明显和附近背景光谱接近，会出现漏提取，同时由于施工时期，邻近修建完成的道路可能被尘土覆盖，而被错误分为正在修建的道路。不过整体精度达到90%以上，可满足实际应用需求。道路施工场景较为复杂，实验结果中仍存在的道路漏提取和错误分类等问题，需在下一步研究中进行改善。