CN101976335A - 基于改进c-v模型的遥感图路网提取方法 - Google Patents

Abstract

本发明公开了一种基于改进C-V模型的遥感图路网提取方法，包括如下步骤：(1)对原始遥感图像预进行处理；(2)选定种子点的并进行首张路网子图的分割；(3)利用改进C-V模型的水平集方法进行路网区域提取；(4)利用形态学细化的方法进行路网中心线的提取；(5)利用子图位置判定规则分割下一幅路网子图、利用阈值分割以及形态学方法自动获取路网子图中的路网初始曲线。本发明可以准确高效地从遥感图中提取路网，能够实时更新城市道路网，有效地辅助决策者进行交通规划的制定。

Description

基于改进C-V模型的遥感图路网提取方法

技术领域

本发明涉及一种基于C-V模型的遥感图路网提取方法，属于高分辨率卫星遥感影像处理与信息提取应用领域。

背景技术

随着我国社会经济的快速发展，城市化进程不断加快，城市人口和机动车拥有量持续增加，交通运输需求和城市道路条件能力之间的矛盾日益突出。目前，我国很多城市交通拥堵现象以及由此衍生引发的交通事故与交通环境污染等社会问题已经极其严重。作为城市命脉的交通已成为城市发展的瓶颈，严重地制约了社会、经济的发展，影响了人民群众的日常工作和生活。因此，加强城市交通规划的研究，建立通畅良好的交通环境，对保证社会和经济的可持续发展具有重要意义。

城市交通规划的关键信息之一是地图。传统电子交通地图和交通网络的采集与获取主要依靠人力测绘纸质地图继而进行扫描矢量化得到的，由于地图测绘工作量大，绘制周期长，交通地图的精度和实时性都会与当前实际交通环境存在一定的差距。随着传感技术和遥感技术的发展，利用遥感影像进行地理数据的采集与获取，使得上述问题的解决成为可能。

正常的高分辨率卫星图像容量都以GB为单位，受限于电脑处理能力，直接处理大幅遥感图像不现实。

目前，高分辨率遥感图像的路网提取算法主要有动态规划、模板匹配、数学形态学、水平集、边缘检测、区域生长等方法，其中基于水平集算法的图像分割方法因其可以对任意复杂的形状进行模型化，且可以隐性解决拓扑形状的分裂、合并等拓扑变化的特点而备受关注。C-V模型就是一种常见的水平集分割模型。然而，传统的水平集C-V模型存在以下不足：一是需要进行迭代中的初始曲线修正，即重新初始化符号距离函数，造成了计算量增大，降低了分割速度，并且分割准确度也难以保证；二是对区域的划分仅仅考虑到图像的RGB空间灰度特征，损失了彩色图像丰富的颜色信息，由于异物同灰度值的现象普遍存在，势必会造成目标的误分割。同时，背景物体颜色的相似也严重地影响了路网的识别效果。因此，传统水平集C-V模型在图像分割速度、效率、准确性方面都存在不足之处，以致传统的基于C-V模型的遥感图路网提取方法在速度、效率、准确性等方面也相应地存在不足。

发明内容

本发明的目的是提供一种基于改进C-V模型的遥感图路网提取方法。

为实现上述目的，本发明所采取的技术方案是：其基于改进C-V模型的遥感图路网提取方法包括以下步骤：

(1)对原始遥感图像进行预处理以提高该遥感图像的清晰度和对比度，得到该遥感图像的预处理图；将所述遥感图像的预处理图置于任意建立的一个平面直角坐标系中；

(2)用户在所述遥感图像的预处理图的路网中选定种子点；根据所选定的种子点获取所述遥感图像的预处理图的路网初始曲线，根据所述路网初始曲线对所述遥感图像的预处理图进行分割获得第一幅路网子图，所述路网子图中含有所述路网；

(3)利用基于C-V模型的水平集方法对当前路网子图进行路网区域提取，获取路网区域子图，所述C-V模型的水平集方法的演化方程如式(1)所示：

$\frac{\∂\mathrm{\φ}}{\∂t}={\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)[v\▿\·\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+{\mathrm{\β}}_1({\mathrm{\λ}}_2{(u_{\mathrm{RGB}}-c_2)}^2-{\mathrm{\λ}}_1{(u_{\mathrm{RGB}}-c_1)}^2)$

$+{\mathrm{\β}}_2({{\mathrm{\λ}}_4(u_{\mathrm{HSI}}-c_4)}^2-{\mathrm{\λ}}_3{(u_{\mathrm{HSI}}-c_3)}^2)]+\mathrm{\μ}(\mathrm{\Δ\φ}-\▿\·\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right))---\left(1\right)$

式(1)中，φ为水平集函数，t为时间，

为φ的偏微分，δ_ε(φ)为Dirac函数，δ_ε(φ)＝H′_ε(φ)，H′_ε(φ)为正则化的Heaviside函数的导数，

为φ二维平面中的梯度，

为梯度范式，u_RGB为RGB空间灰度值，u_HSI为HIS空间混合通道值，c₁、c₂分别表示曲线内、外平均灰度值，v、λ₁、λ₂为控制演化曲线收缩或者膨胀的权重参数，c₃、c₄分别为HSI空间轮廓内、外的平均灰度值，λ₃、λ₄为HIS空间权重参数，β₁、β₂为RGB空间灰度信息与HSI空间混合通道信息的权重系数，μ为距离函数权重参数；

(4)利用形态学细化方法对步骤(3)所述的路网区域子图进行路网中心线提取，获取路网中心线子图；

(5)利用子图位置判定规则判断所述遥感图像的预处理图中是否存在下一幅路网子图：

若存在下一幅路网子图，则利用所述子图位置判定规则从所述遥感图像的预处理图中分割得到所述下一幅路网子图，再利用自适应灰度阈值、饱和度阈值分割以及形态学方法自动获取该下一幅路网子图中的路网初始曲线，然后对该下一幅路网子图依次执行步骤(3)和步骤(4)得到该下一幅路网子图相应的路网区域子图和路网中心线子图；

若不存在下一幅路网子图，则利用各路网子图的路网中心线子图生成路网图。

进一步地，本发明在步骤(5)中，利用子图位置判定规则判断所述遥感图像的预处理图中是否存在下一幅路网子图的方法为：在所述平面直角坐标系中，如果当前路网中心线子图中的路网中心线的末端点的横坐标和纵坐标对应地小于遥感图像的预处理图的像素点的最大横坐标和最大纵坐标，则存在下一幅路网子图；否则，不存在下一幅路网子图；

利用所述子图位置判定规则从所述遥感图像的预处理图中分割得到所述下一幅路网子图的方法为：

在所述平面直角坐标系中，用直线连接当前路网中心线子图中路网中心线的首端点和末端点，计算该直线的斜率以及该直线与横坐标的夹角θ，并以该当前路网中心线子图中路网中心线的末端点为圆心生成圆，根据夹角θ的大小和所述圆分割下一幅路网子图，该下一幅路网子图为正方形框图并与夹角θ和所述圆满足以下关系：

如果-22.5°≤θ＜22.5°，则所述圆内切于下一幅路网子图的左边，切点为该左边的中点；

如果22.5°≤θ＜67.5°，则所述圆同时内切于下一幅路网子图的左边和下边；

如果67.5°≤θ＜112.5°，则所述圆内切于下一幅路网子图的下边，切点为该下边的中点；

如果112.5°≤θ＜157.5°，则所述圆同时内切于下一幅路网子图的下边和右边；

如果157.5°≤θ＜180°或者-180°≤θ＜-157.5°，则所述圆内切于下一幅路网子图的右边，切点为该右边的中点；

如果-157.5°≤θ＜-112.5°，则所述圆内切于下一幅路网子图的上边和右边；

如果-112.5°≤θ＜-67.5°，则所述圆内切于下一幅路网子图的上边，切点为该上边的中点；

如果-67.5°≤θ＜-22.5°，则所述圆内切于下一幅路网子图的上边和左边。

与现有路网提取方法相比，本发明的有益效果是：

(1)本发明通过在传统的C-V模型水平集方法中引入惩罚函数避免了重新初始化，加快了曲线迭代过程，提高了路网提取的速度和效率；

(2)本发明在传统的C-V模型水平集方法中整合RGB空间和HSI空间的各通道特性，有效地抑制边缘噪声的干扰，提高了路网提取的准确性；

(3)本发明利用子图位置判定规则将遥感图分割成路网子图，使提取的目标路网集中在少量的路网子图中，明显地消除了无目标背景对路网提取的干扰，提高路网提取的准确性，同时大大降低了数据处理量，加快了路网提取过程；

(4)本发明通过自适应灰度阈值、饱和度阈值分割以及形态学方法实现了路网子图初始曲线的自动获取，提高路网识别的速度和效率；

(5)本发明的遥感图路网提取方法可以直接从遥感图中提取路网，相比地图测绘成本低、周期短、实时性高。

附图说明

图1为遥感图路网提取流程图；

图2为基于子图位置判定规则的路网子图分割示意图；

图3为子图位置判定规则示例图；

图4为本发明提取方法利用子图位置判定规则前、后的道路提取结果比较图：其中，(a)为原始路网图，图中方框为利用子图位置判定规则生成的路网子图，(b)为利用子图位置判定规则前的道路提取结果示意图，(c)为利用本发明的改进C-V模型和子图位置判定规则后的道路提取结果示意图；

图5为利用本发明改进C-V模型水平集方法实现高速公路初始曲线自动获取的过程示意图：其中，(a)为待分割原图；(b)为利用灰度阈值和饱和度阈值进行自适应双阈值分割的结果；(c)为阈值分割后进行形态学操作的结果；(d)为对形态学操作结果轮廓跟踪的结果；

图6为利用本发明提取方法实现城区路网初始曲线自动获取的过程示意图：其中，(a)为待分割原图；(b)为利用灰度阈值和饱和度阈值进行自适应双阈值分割的结果；(c)为阈值分割后进行形态学操作的结果；(d)为对形态学操作结果轮廓跟踪的结果；

图7为不同路网提取方法的道路提取结果对比图：其中，(a)为原始路网图，图中正方形框是初始化曲线传播过程中生成的路网子图，(b)为灰度图二值分割的结果，(c)为灰度图的区域生长结果，(d)为基于灰度图的水平集分割结果，(e)为基于HIS空间的C-V分割结果，(f)为本发明的提取结果。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步的详细描述。

以下以广东省中山市地区的遥感影像图(QuickBird公司提供、精度0.61m)为例具体说明本发明提取方法的各步骤，并以完整性和正确性指标对路网提取结果进行评价。

为了量化地评价路网提取算法，需要将从遥感影像图中提取的路网区域结果(Estimate)和实际测绘所得的路网区域结果(GT)进行比较，并引入完整性(c)和正确性(p)两个指标进行评价。式(2)分别给出了完整性(c)和正确性(p)的定义

$\left\{\begin{array}{c}c=\frac{|\mathrm{GT}\∩\mathrm{Estimate}|}{\left|\mathrm{GT}\right|}\×100\%\\ p=\frac{|\mathrm{GT}\∩\mathrm{Estimate}|}{\left|\mathrm{Estimate}\right|}\×100\%\end{array}\right.---\left(2\right)$

式(2)中，|·|表示区域像素点数。完整性是由遥感图提取得到的路网区域与实际测绘所得路网区域相匹配部分的面积和测绘所得路网区域总面积的百分比；正确性是实际由遥感图提取得到的路网区域与测绘所得路网区域结果相匹配部分的面积和实际遥感图提取得到的路网区域总面积的百分比。通过计算完整性和正确性，进行路网提取算法评价，为本发明路网提取方法的优劣提供了一个量化评价。

本发明基于改进C-V模型的遥感图路网提取方法的流程如图1所示，以下举例说明具体步骤：

(1)对原始遥感图像预进行处理

对原始遥感图像预进行处理，得到该遥感图像的预处理图。对原始遥感图预处理包括图像几何变换、图像增强、颜色变换、阈值分割以及绘图板等功能。通过预处理实现24位真彩图到256色位图的转换、RGB与HSI空间变换，同时可以对图像进行缩放和旋转，从而减少图像的噪声，提高图像清晰度，增强目标与背景的对比度，以取得更好的分割识别效果。

将以上得到的遥感图像的预处理图置于任意建立的一个平面直角坐标系中。

(2)种子点的选定及首张路网子图的分割

用户在遥感图像的预处理图的待提取的目标路面区域选定至少一个种子点。根据用户选定的种子点，生成以种子点为圆心、半径为15pixels的圆(路宽的1/2左右)，将该圆作为遥感图像的预处理图的路网初始曲线，根据所述路网初始曲线对所述遥感图像的预处理图进行分割获得第一幅路网子图，第一幅路网子图大小为256×256pixels，该路网子图中含有所述路网。

(3)基于C-V模型的水平集方法的路网区域提取

获得路网子图后，本发明利用改进C-V模型的水平集方法来实现路网子图的路网区域的提取。

1)基于改进C-V模型的子图路网区域提取

目前，遥感图路网提取以水平集方法为主，水平集方法实质上将曲线运动转化为曲面运动的过程，在图像分割应用中即使目标边界分裂或者合并，但在高一维的曲面上其拓扑结构并不改变。常见的水平集分割模型有Chan和Vese在Mumford-Shah模型基础上改进得到的C-V模型。该模型不依赖局部图像的梯度信息，能够达到全局最优分割，可以处理复杂多变的拓扑结构。其能量函数通常包括闭合曲线C的长度项，C内部形成的面积项以及两个拟合项，简化为

$E(C,c_1,c_2)=v\·\mathrm{length}\left(C\right)+{\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}{|u(x,y)-c_1|}^2\mathrm{dxdy}+{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|u(x,y)-c_2|}^2\mathrm{dxdy}---\left(3\right)$

式(3)中，C是一条封闭曲线(活动轮廓)，Ω表示图像区域，

{x，y}∈Ω，length(C)为闭合曲线C的长度项，u(x，y)为原图像，c₁和c₂分别表示曲线内外平均灰度值，权重参数v≥0，λ₁，λ₂＞0。通过对权重参数的调整可以实现演化曲线收缩或者膨胀，从而放宽了对初始曲线位置选取的限制。

在传统的C-V模型水平集方法中，利用零水平集函数替代

同时引入正则化的Heaviside函数

式(3)的能量泛函可表示为

$E(\mathrm{\φ},c_1,c_2)=v{\∫}_{\mathrm{\Ω}}{|\▿H_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)|}^2\mathrm{dxdy}+{\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}{|u(x,y)-c_1|}^2{H}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)\mathrm{dxdy}$

$+{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|u(x,y)-c_2|}^2(1-H_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right))\mathrm{dxdy}---\left(4\right)$

其中Dirac函数δ_ε(φ)＝H′_ε(φ)＝ε/π(ε²+φ²)。c₁和c₂可通过下式计算获得

$\left\{\begin{array}{cc}{c}_1=\frac{\&Integral;u{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}{\&Integral;H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)\mathrm{dxdy}}& \{\mathrm{\&phi;}\&GreaterEqual;0\}\\ c_2=\frac{\&Integral;u(1-H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}{\&Integral;(1-H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right))\mathrm{dxdy}}& \{\mathrm{\&phi;}<0\}\end{array}\right.---\left(5\right)$

传统的C-V模型水平集方法中，为了保证水平集函数φ演化稳定，需要进行重新初始化处理，否则演化过程会产生剧烈的震荡，导致图像分割效果不理想；在处理遥感图像时，直接表现为道路轮廓不光滑，边缘参差不齐，路面提取不完整甚至偏离路面。重新初始化作为一种数值修正，其代价是昂贵的。为了避免重新初始化的问题，同时又保证水平集演化过程的稳定，本发明引入距离函数惩罚项替代重新初始化过程。距离函数惩罚项是水平集函数与符号距离函数偏差的一个度量函数。在区域

该项能够使演化曲线自动保持为大概的距离函数，有效控制水平集与符号函数的偏差，用积分表示如下

$P\left(\mathrm{\&phi;}\right)={\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dx}---\left(6\right)$

增加距离函数惩罚项后可以得到改进的C-V模型的能量函数为

E(φ)＝E(φ，c₁，c₂)+μP(φ)(7)

其中μ为距离函数权重参数，于是改进后的能量函数表示为

$E\left(\mathrm{\&phi;}\right)={\mathrm{\&lambda;}}_1{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_1|}^2{H}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)\mathrm{dxdy}+{\mathrm{\&lambda;}}_2{\&Integral;}_{\mathrm{\&Omega;}}{|u(x,y)-c_2|}^2(1-H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right))\mathrm{dxdy}$

$+v{\&Integral;}_{\mathrm{\&Omega;}}{|\&dtri;H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)|}^2\mathrm{dxdy}+{\mathrm{\&mu;}\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^2\mathrm{dxdy}---\left(8\right)$

最小化能量函数E(φ)，得到相应的Euler-Lagrange方程，并由水平集方法得到

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)[v\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)-{\mathrm{\&lambda;}}_1{(u(x,y)-c_1)}^2+{\mathrm{\&lambda;}}_2{(u(x,y)-c_2)}^2]+\mathrm{\&mu;}(\mathrm{\&Delta;\&phi;}-\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right))---\left(9\right)$

式(9)中，t为时间，

为φ的偏微分，

为φ二维平面中的梯度，

为梯度范式。通过对输入图像u(x，y)迭代计算得到曲线演化的最终结果，并以此进行图像分割，提取路网区域。

本发明通过增加距离惩罚项改进传统的C-V模型水平集方法，解决了水平集函数在演化过程中需要重新初始化的问题，加快了曲线迭代过程，并能有效抑制区域噪声干扰，保持复杂拓扑结构。从而在进行遥感图路网提取时完全能够针对具有较弱边界的复杂图像进行目标分割，比如对于场景复杂的卫星图像，可以对其中的道路目标实现较为光滑的轮廓提取。所以，本发明通过在传统的C-V模型水平集方法中增加距离惩罚项可以提高路网识别的速度、效率和准确性。

通常，式(9)的输入图像u(x，y)是RGB空间的灰度信息u_RGB，RGB空间是最常用的颜色空间，各通道值可以直接从BMP位图获取得到。目前，常常将彩色图像R、G、B三个通道信息线性合并获取灰度信息，并以此作为进一步图像处理的基础，由于灰度图简单方便，得到了广泛的应用，许多图像处理方面的研究算法都是针对灰度图提出。然而简单的将RGB各通道信息线性相加，势必造成信息的损失，导致异物同灰度值，造成目标的误分割。为此，本发明为式(9)模型又引进了另一种常用的颜色空间HSI(hue，saturation，intensity)空间。H为色调，表示不同颜色种类；S为饱和度，表示颜色深浅；I为亮度，表示颜色的明亮。Gonzalez提供RGB到HSI空间的转换算法，如式(10)所示

$\left\{\begin{array}{c}H=\left\{\begin{array}{cc}\mathrm{\&theta;}& B\&le;G\\ 360-\mathrm{\&theta;}& B\&GreaterEqual;G\end{array}\right.\\ S=1-\frac{3}{(R+G+B)}\left[\min (R,G,B)\right]\\ I=\frac{1}{3}(R+G+B)\\ \mathrm{\&theta;}=\arccos \left\{\frac{1/2[(R-G)+(R-B)]}{{(R-G)}^2+(R-G)\&times;(G-B)}\right\}\end{array}\right.---\left(10\right)$

HSI颜色模型中H分量是确定颜色的关键因素，S分量越大(接近于1)，颜色越纯，S分量越小(接近于0)，颜色越接近纯灰色，I分量类似RGB空间的灰度信息。H和S与人感受颜色的方式紧密相关，彩色图像中每个均匀区域对应相对一致的色调和饱和度。I分量和彩色信息无关。HSI模型的三个基本属性符合人类对颜色的描述和解释方式，提供了人类视觉对事物色彩的直接认识。

由于H分量是确定颜色的主要因素，常被单独用作水平集方法的输入信息。当亮度I很大接近于白色或很小接近于黑色，色调H和饱和度S变化不大。鉴于H分量、S分量与亮度信息I关系不大，在不均匀光照的情况下，例如阴影、树荫，利用H分量、S分量进行分割特别具有优势。但是H分量在R＝G＝B时具有奇异点，造成无穷大的跳跃，颜色可能呈现不连续。如果遥感图像路面区域的饱和度较低，色调较高，绿色植被区域的饱和度较高，色调较低，若单独采用H分量或S分量，会丢失部分信息，为此，考虑同时采用S分量和H分量进行分割，可采用下式对H分量和S分量信息进行整合

式(11)中，k为可调参数，由式(11)可知，图像饱和度较大，色调较高时，u_HSI较大，反之则较小。式(11)将像素点的颜色色调和饱和度信息综合在一起得到混合通道值，突出了路面的颜色特征。将u_HSI作为水平集方法的另一个输入图像整合到式(9)中，得到基于多颜色空间无需重新初始化的Hamilton-Jacobi方程为

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)[v\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+{\mathrm{\&beta;}}_1({\mathrm{\&lambda;}}_2{(u_{\mathrm{RGB}}-c_2)}^2-{\mathrm{\&lambda;}}_1{(u_{\mathrm{RGB}}-c_1)}^2)$

$+{\mathrm{\&beta;}}_2({\mathrm{\&lambda;}}_4{(u_{\mathrm{HSI}}-c_4)}^2-{\mathrm{\&lambda;}}_3{(u_{\mathrm{HSI}}-c_3)}^2)]+\mathrm{\&mu;}(\mathrm{\&Delta;\&phi;}-\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right))---\left(12\right)$

式(12)即为本发明改进型的C-V模型水平集方法的演化方程。

式(12)满足初始条件φ(x，y，t)|_t＝0＝φ₀(x，y)。其中u_RGB为RGB空间灰度值，

u_HSI根据式(11)计算得到；c₃，c₄是类似c₁，c₂，为HSI空间轮廓内外的均值；β₁，β₂为RGB空间灰度信息与HSI空间混合通道信息的权重系数。式(12)利用惩罚函数避免了重新初始化，并综合RGB空间和HSI空间的各通道特性，通过调整β₁，β₂改变色调饱和度信息和灰度信息的权重，可以提高了演化速度，抑制边缘噪声的干扰，对边缘离散、背景复杂、色彩信息丰富的图像具有理想的分割效果，从而提高路网提取的准确性。

在式(12)所示的模型中，由于所涉及参数过多，限制了该模型的应用能力，同时也降低了该模型的执行效率。因此，根据实际遥感图的特点，设定该模型中参数λ₁＝λ₂＝λ₃＝λ₄＝1，μ＝3。这样该模型中只保留了参数v、β₁、β₂。参数简化后进行离散化得到

$\left\{\begin{array}{c}{\mathrm{\&phi;}}^{n+1}={\mathrm{\&phi;}}^n+\mathrm{F\&Delta;t}\\ {\mathrm{\&phi;}}^0:x_n^2+y_n^2=15\end{array}\right.---\left(13\right)$

式(13)中，Δt表示演化曲线的时间步长；F表示式(12)中等式右边项；其中φ⁰为初始曲线，本文选定初始曲线为半径15pixels的圆(路宽的1/2左右)。利用式(13)的基于C-V模型的水平集方法对当前路网子图进行迭代运算以提取路网区域，获取路网区域子图。

2)区域滤波

利用改进C-V模型的水平集方法对当前路网子图进行路网区域提取后，得到的路网区域子图会出现一些大小不等图斑。为了提高路网区域提取的准确性，在利用改进C-V模型的水平集方法对当前路网子图进行路网区域提取后，可以进行区域滤波，滤除非路网区域的图斑。对于小图斑，采用下面办法：标出子图中所有区域，计算各区域的面积S(像素点数)，并与设定面积阈值TH1比较，若小于TH1则滤除该区域，否则保留。对于大的图斑，利用道路的形状特征进行去除。道路的形状特征指数t可以用以下公式表示

$t=\frac{P^2}{4\mathrm{\&pi;S}}---\left(14\right)$

式(14)中，P为路网区域周长(区域轮廓像素点数)，S为路网区域面积(像素点数)。形状指数t，对圆为1，正方形为4/π。遥感图像中道路区域为带状，近似为一个具有长宽的矩形，形状指数t较大，将t与形状指数阈值TH2进行比较，若小于TH2则滤除该区域，否则保留。根据以上可知，只有同时满足面积阈值和形状指数阈值条件的区域才保留。

通过区域滤波可以滤除当前路网区域子图中非路网区域的图斑，从而使当前提取出来的路网区域子图更加准确。

(4)基于形态学细化的路网区域子图的路网中心线的提取

利用改进C-V模型提取出路网区域后，为了便于后续的交通流量预测，需要提取路网中心线。路网中心线的提取过程本质上是求取路网骨架的过程，反映了道路特性。本发明采用形态学细化方法来获取路网中心线。形态学细化方法的基本思想是：在给定一系列的具有一定形状的结构元素(模板)后，顺序地删除满足某种或者几种变换的像素。

形态学细化的关键在于寻找合适的结构元素，亦即结构元素决定了算法的目的和性能。结构元素保证了目标在每次迭代过程中结构的连通性和整个图形的构造不变。在细化过程中，始终要遵循以下原则：(1)内部点不能删除；(2)孤立点不能删除；(3)直线端点不能删除；(4)如果边界点删除后不影响连通分量，该边界点才能被删除。

(5)子图位置判定规则与路网初始曲线的自动获取

1)基于子图位置判定规则的遥感图像分割

正常的遥感图容量都以GB为单位，受限于电脑处理能力，直接处理大幅遥感图像不现实，本发明利用子图位置判定规则对遥感图像的预处理图分块，在大幅遥感图像的预处理图中自动分割出包含路网的小幅遥感图像，然后对各个路网子图进行路网提取。

首先利用子图位置判定规则判断遥感图像的预处理图中是否存在下一幅路网子图，方法如下：

在步骤(1)建立的平面直角坐标系中，如果当前路网中心线子图中的路网中心线的末端点的横坐标和纵坐标对应地小于遥感图像的预处理图的像素点的最大横坐标和最大纵坐标，则存在下一幅路网子图；否则不存在下一幅路网子图。

如果存在下一幅路网子图，则利用子图位置判定规则来判定下一幅路网子图的位置，并分割出下一幅路网子图；如果不存在下一幅路网子图，则利用各路网子图的路网中心线子图生成路网图。

以下具体说明利用子图位置判定规则判断遥感图像的预处理图中是否存在下一幅路网子图的方法。

本发明中子图位置判定规则是在路网中心线接力引导的模式下，根据每个路网中心线子图得到下一幅路网子图位置的过程，其示意图如图2所示。图2中，两条实线表示道路的边缘，虚线表示路网中心线，虚线框为生成的路网子图位置，虚线框A表示当前路网子图，虚线框B、C、D、E、F是利用子图位置判定规则分割得到的路网子图，圆a为以当前路网子图的路网中心线首端点为圆心、半径为15pixels的圆，圆b、c、d、e、f是以各个路网子图中路网中心线末端点为圆心、半径为15pixels的圆。

由图2可知，子图位置的获取的两个关键步骤：一是路网中心线的提取，二是如何根据路网中心线得到相应路网子图的位置，即子图位置判定规则。路网中心线的提取在上一步骤中已经详细说明，以下结合图3详细说明本发明的子图位置判定规则。图3中，I_n和I_n-1为相邻的两个路网子图，其中，I_n-1为当前路网子图，I_n为根据子图位置判定判定规则由I_n-1分割得到的下一幅路网子图。O_n-1(x_n-1，y_n-1)、O_n(x_n，y_n)分别为当前路网子图I_n-1的路网中心线的首端点和末端点，圆a_n-1为以点O_n-1(x_n-1，y_n-1)为圆心、半径为15pixels的圆，圆a_n为以点O_n(x_n，y_n)为圆心、半径为15pixels的圆，l为经过点O_n-1(x_n-1，y_n-1)与O_n(x_n，y_n)的直线，l₁为经过O_n-1(x_n-1，y_n-1)的水平轴方向参考线；θ_n为路网中心线l与水平参考线l₁的夹角。

路网子图位置获取的过程是从圆a_n-1接力传播到a_n，并依据夹角θ_n值大小范围来生成子图I_n的过程，θ_n、a_n与I_n的关系就是子图位置判定的规则。利用θ_n来判定路网子图生成位置的前提是，图像中l与路网中心线近似接近，l的斜率能够反映路网子图中路网中心线的趋势。这里，通过选择路网子图的规模大小来满足这个要求，即路网子图不宜过大或过小，保证在高分辨率卫星子图像中道路呈现近似矩形，路网子图中尽量只包含一条道路，从而l的斜率能够反映图像中道路趋势；另外，路网子图之间设置重叠区域，弥补道路切割可能造成图像边界的畸变和目标漏失。由于道路是相互连通的，势必有些路网子图中包含十字交叉路口。处理这类情况的方法是，以上一幅路网子图的斜率来判定路网走势，然后找到该路网子图中相应走势的路网中心线的末端点，以该端点和前一子图的θ_n-1生成下一子图的位置。

夹角θ_n可由计算得到。

路网子图的大小根据道路宽度进行设置。本实施例中，路网子图的大小可设置为256×256pixels。根据夹角θ_n的大小和圆a_n分割下一幅路网子图，该下一幅路网子图与夹角θ_n和圆a_n满足以下关系：

如果-22.5°≤θ_n＜22.5°，则圆a_n内切于下一幅路网子图I_n的左边，切点为该左边的中点；

如果22.5°≤θ_n＜67.5°，则圆a_n同时内切于下一幅路网子图I_n的左边和下边；

如果67.5°≤θ_n＜112.5°，则圆a_n内切于下一幅路网子图I_n的下边，切点为该下边的中点；

如果112.5°≤θ_n＜157.5°，则圆a_n同时内切于下一幅路网子图I_n的下边和右边；

如果157.5°≤θ_n＜180°或者-180°≤θ＜-157.5°，则圆a_n内切于下一幅路网子图I_n的右边，切点为该右边的中点；

如果-157.5°≤θ_n＜-112.5°，则圆a_n内切于下一幅路网子图I_n的上边和右边；

如果-112.5°≤θ_n＜-67.5°，则圆a_n内切于下一幅路网子图I_n的上边，切点为该上边的中点；

如果-67.5°≤θ_n＜-22.5°，则圆a_n内切于下一幅路网子图I_n的上边和左边。

图3中θ_n近似于20度，则圆a_n内切于下一幅路网子图I_n的左边，切点为该左边的中点。

基于子图位置判定规则的遥感图分割方法使待提取的目标路网集中在少量子图中，明显消除了无目标背景对路网提取的干扰，抑制了边缘噪声的影响，大大降低了数据处理量，提高了得到路网识别的准确性。

图4显示了是否基于子图位置判定规则进行道路提取的结果比较。图4(a)为利用子图位置判定规则实现图像分割结果示意图，图像大小为1024×1024pixels，图中位于左下方的圆A是由种子点生成的初始曲线，路网子图大小均为256×256pixels；图4(b)中为未用子图判定规则进行图像分割方法的水平集路网提取的结果；图4(c)为本发明提出的基于改进C-V模型水平集方法的遥感图路网提取结果。图4(c)在道路提取过程中，利用子图位置判定规则将图像分割成5张256×256pixels的路网子图分别进行路网提取。

计算图4(b)和图4(c)算法的评价指标得到，图4(b)所示方法的完整性和正确性分别为80.87％、62.25％，图4(c)所示方法的完整性和正确性为96.62、96.90％；并且图4(b)采用的算法耗时是图4(c)的4倍左右。根据以上可以得到，图4(c)采用的本发明提出基于子图位置判定规则的水平集路网提取方法可以明显消除无目标背景对路网提取的干扰，抑制了边缘噪声的影响，大大降低了数据处理量，加快了路网提取过程，并且基于子图位置判定规则的路网提取将整图路网提取结合到一起，获得光滑完整的道路轮廓。

2)路网初始曲线的自动获取

分割出路网子图后，需要获取路网子图的初始化曲线，然后再利用改进C-V模型水平集方法对当前路网子图进行路网区域提取。

目前，C-V模型水平集方法通常采用人工选定以及预先设置的办法得到初始曲线位置，在应用于图像分割、目标提取时，通常采用固定阈值分割办法获得演化初始轮廓。C-V模型水平集方法的初始化曲线的设置不当会明显影响模型的分割效率，甚至导致分割失败。高分辨率遥感图中图像背景复杂，灰度不均匀，路面灰度和一些建筑物的灰度近似，如果初始曲线的位置不当，会导致迭代运算时间过长，容易造成轮廓泄露，甚至无法分割出路网区域。为避免出现上述现象，合理选择初始曲线成为了路网区域提取首先需要解决的问题。

为实现初始曲线的自动获取，本发明采用自适应阈值算法动态计算遥感图像的灰度值阈值和饱和度阈值，利用灰度阈值和饱和度阈值对图像进行预分割，并对分割结果执行形态学操作，跟踪形态学处理结果得到轮廓，将该轮廓作为水平集演化的初始曲线。

彩色图像一般由R、G、B三个分量组成，一般图像处理时都是针对灰度图进行。通过考察遥感图像中道路三个分量的特点，可以发现RGB空间G分量比较清晰干净。为了充分地利用彩色图像各通道信息，进行RGB空间到HSI空间变换，发现HSI空间道路区域的饱和度S分量较低。因此采用G分量和饱和度S分量进行双阈值分割。阈值分割是实现图像二值化的一种手段，某点灰度值大于阈值则成为背景点，其余的则设置为目标点，加入目标区域。本发明阈值采用自适应阈值计算方法，通过迭代方法来计算灰度阈值和饱和度阈值。只有图像中像素点灰度值Z(i，j)同时小于等于灰度阈值T_G和饱和度阈值T_S，才将其加入目标区域。

自适应阈值是对一个邻域窗口内像素点进行的阈值计算，窗口大小根据实际图像进行调整，某一窗口内阈值计算步骤如下：

①求出图像中最小灰度值Z_min和最大灰度值Z_max，令

②根据阈值

(k＝0，1，2…)将图像分成目标和背景两部分，计算两部分的灰度均值：

$Z_O=\frac{\underset{Z(i,j)\&le;T_G^k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)\&le;T_G^k}{\mathrm{\&Sigma;}}N(i,j)}---\left(15\right)$

$Z_B=\frac{\underset{Z(i,j)>T_G^k}{\mathrm{\&Sigma;}}Z(i,j)\&times;N(i,j)}{\underset{Z(i,j)>T_G^k}{\mathrm{\&Sigma;}}N(i,j)}---\left(16\right)$

其中，Z(i，j)是图像中点(i，j)的灰度值，N(i，j)是点(i，j)的权重值，一般令N(i，j)＝1.0。

③更新阈值

$T_G^{k+1}=\frac{Z_O+Z_B}{2}---\left(17\right)$

如果

则为最终阈值，否则k++，并回到步骤②继续迭代计算。

利用上述方法计算灰度阈值T_G和饱和度阈值T_S，并根据Z(i，j)同时小于等于T_G和T_S对路网子图进行双阈值粗分割。粗分割得到的结果一般会有许多噪声和杂质，目标可能不连通，或者过度分割，因此需要进一步处理。由于对初始曲线位置要求不是很精准，在此运用形态学操作，采用米字形结构元素对分割结果进行腐蚀操作，去除小杂质噪声和欠分割区域，然后进行膨胀操作，去除目标路网过度分割造成的内部空洞，得到大概的能够连通目标路网的路网区域。最后采用8邻域轮廓跟踪得到接近路网区域的初始轮廓线，完成C-V模型水平集演化初始曲线的自动获取。

路网初始曲线自动获取实验中，待分割遥感图像有单独包含高速公路的简单图像和包含复杂场景的城区图像，其初始曲线自动获取过程分别如图5、图6所示。图5和图6中，(a)为待分割原图；(b)为利用灰度阈值和饱和度阈值进行自适应双阈值分割的结果；(c)为阈值分割后进行形态学操作的结果；(d)为对形态学操作结果轮廓跟踪的结果。

图5和图6显示了自适应双阈值方法得到的粗分割结果并不是很准确，但是基本实现了路网的定位；形态学操作可以消除小面积块的干扰，填补路网区域内部的空洞，得到连通的路面。尽管图6得到的轮廓曲线并不是道路精准的位置，但是比较靠近道路边缘，实现了路网子图初始曲线的自动获取，提高路网识别的速度和效率，为后续水平集演化提供了初始演化曲线。

3)遥感图像整图路网提取

获取路网子图初始曲线后：

首先利用步骤(3)的改进C-V模型的水平集方法来实现路网子图的路网区域的提取；

然后再利用步骤(4)的形态学细化方法实现路网区域子图的路网中心线的提取；

再利用步骤(5)中的子图位置判定规则判断是否存在下一幅路网子图，若存在下一幅路网子图，则利用步骤(5)中的子图位置判定规则分割出下一幅路网子图，再次利用步骤(5)中的自适应灰度阈值、饱和度阈值分割以及形态学方法自动获取下一幅路网子图中的路网初始曲线，最后再跳转至步骤(3)顺序重复操作；若不存在下一幅路网子图，则利用各路网子图的路网中心线子图生成路网图。

图7显示了不同路网提取方法在场景复杂情况下的道路提取结果比较，图像大小为1024×1024pixels，图像中道路受到车辆、白线、楼房等噪声影响。图7(a)为待提取路网的原图，图中正方形框是初始化曲线传播过程中生成的路网子图，大小均为256×256pixels；图7(b)为灰度图二值分割的结果；图7(c)为灰度图的区域生长结果；图7(d)为基于灰度图的水平集分割结果，图7(e)为基于HIS空间的C-V分割结果；图7(f)为本发明改进C-V模型水平集方法的提取结果。

图7(b)、7(c)、7(d)中路网区域与路边灰度值相似的建筑物连为一体，道路边缘不连续，拓扑结构受到破坏，无法获得道路准确的位置。图7(e)中得到了路网区域的基本形状，但是提取得到的道路边缘非常不光滑，有小突起存在。相比图7(e)，图7(f)获得的路网区域边缘比较完整光滑，受背景噪声影响较少。为了更好地对比各种不同路网提取方法，利用式(2)计算图7(a)、7(b)、7(c)、7(d)、7(e)、7(f)中应用的各类图像提取方法的评价指标，具体数值如表1所示。

表1各类路网提取方法的评价指标

图7	所用提取方法	完整性(c)	准确性(p)
				b	二值分割	70.43	27.61
c	区域生长	76.58	42.93
				d	灰度图的水平集分割	74.47	26.68
e	HIS空间C-V模型分割	89.08	90.25
				f	本发明路网提取方法	89.92	93.69

表1显示了图7(e)和图7(f)得到的路网路面的完整性和正确性高于其他两种方法，图7(f)方法的完整性和正确性高于图7(e)，这表明了本发明改进C-V模型水平集的路网提取方法在处理拓扑变化时较为优越，引入的距离惩罚项保证了水平集演化的稳定，得到更光滑的道路边缘，RGB空间和HIS空间各通道的融合信息能够比较准确全面地反应了目标路网的特征。综合图7和表1，可以看出，图7(f)采用本发明改进C-V模型水平集方法结合路网子图位置判定规则可以很好地抑制区域噪声和背景信息的干扰，能够处理高分辨率卫星图像中路网的拓扑变化，路面不连续，边缘模糊的情况，在受房屋、树荫、白线、车辆影响的复杂路段也能获得光滑完整的路网区域。

Claims (2)

1.一种基于改进C-V模型的遥感图路网提取方法，其特征是包括以下步骤：

(1)对原始遥感图像进行预处理以提高该遥感图像的清晰度和对比度，得到该遥感图像的预处理图；将所述遥感图像的预处理图置于任意建立的一个平面直角坐标系中；

(2)用户在所述遥感图像的预处理图的路网中选定种子点；根据所选定的种子点获取所述遥感图像的预处理图的路网初始曲线，根据所述路网初始曲线对所述遥感图像的预处理图进行分割获得第一幅路网子图，所述路网子图中含有所述路网；

(3)利用基于C-V模型的水平集方法对当前路网子图进行路网区域提取，获取路网区域子图，所述C-V模型的水平集方法的演化方程如式(1)所示：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}={\mathrm{\&delta;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)[v\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)+{\mathrm{\&beta;}}_1({\mathrm{\&lambda;}}_2{(u_{\mathrm{RGB}}-c_2)}^2-{\mathrm{\&lambda;}}_1{(u_{\mathrm{RGB}}-c_1)}^2)$

$+{\mathrm{\&beta;}}_2({{\mathrm{\&lambda;}}_4(u_{\mathrm{HSI}}-c_4)}^2-{\mathrm{\&lambda;}}_3{(u_{\mathrm{HSI}}-c_3)}^2)]+\mathrm{\&mu;}(\mathrm{\&Delta;\&phi;}-\&dtri;\&CenterDot;\left(\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right))---\left(1\right)$

式(1)中，φ为水平集函数，t为时间，

为φ的偏微分，δ_ε(φ)为Dirac函数，δ_ε(φ)＝H′_ε(φ)，H′_ε(φ)为正则化的Heaviside函数的导数，

为φ二维平面中的梯度，

为梯度范式，u_RGB为RGB空间灰度值，u_HSI为HIS空间混合通道值，c₁、c₂分别表示曲线内、外平均灰度值，v、λ₁、λ₂为控制演化曲线收缩或者膨胀的权重参数，c₃、c₄分别为HSI空间轮廓内、外的平均灰度值，λ₃、λ₄为HIS空间权重参数，β₁、β₂为RGB空间灰度信息与HSI空间混合通道信息的权重系数，μ为距离函数权重参数；

(4)利用形态学细化方法对步骤(3)所述的路网区域子图进行路网中心线提取，获取路网中心线子图；

(5)利用子图位置判定规则判断所述遥感图像的预处理图中是否存在下一幅路网子图：

若存在下一幅路网子图，则利用所述子图位置判定规则从所述遥感图像的预处理图中分割得到所述下一幅路网子图，再利用自适应灰度阈值、饱和度阈值分割以及形态学方法自动获取该下一幅路网子图中的路网初始曲线，然后对该下一幅路网子图依次执行步骤(3)和步骤(4)得到该下一幅路网子图相应的路网区域子图和路网中心线子图；

若不存在下一幅路网子图，则利用各路网子图的路网中心线子图生成路网图。

2.根据权利要求1所述的基于改进C-V模型的遥感图路网提取方法，其特征是：在步骤(5)中，

利用子图位置判定规则判断所述遥感图像的预处理图中是否存在下一幅路网子图的方法为：在所述平面直角坐标系中，如果当前路网中心线子图中的路网中心线的末端点的横坐标和纵坐标对应地小于遥感图像的预处理图的像素点的最大横坐标和最大纵坐标，则存在下一幅路网子图；否则，不存在下一幅路网子图；

利用所述子图位置判定规则从所述遥感图像的预处理图中分割得到所述下一幅路网子图的方法为：

在所述平面直角坐标系中，用直线连接当前路网中心线子图中路网中心线的首端点和末端点，计算该直线的斜率以及该直线与横坐标的夹角θ，并以该当前路网中心线子图中路网中心线的末端点为圆心生成圆，根据夹角θ的大小和所述圆分割下一幅路网子图，该下一幅路网子图为正方形框图并与夹角θ和所述圆满足以下关系：

如果-22.5°≤θ＜22.5°，则所述圆内切于下一幅路网子图的左边，切点为该左边的中点；

如果22.5°≤θ＜67.5°，则所述圆同时内切于下一幅路网子图的左边和下边；

如果67.5°≤θ＜112.5°，则所述圆内切于下一幅路网子图的下边，切点为该下边的中点；

如果112.5°≤θ＜157.5°，则所述圆同时内切于下一幅路网子图的下边和右边；

如果157.5°≤θ＜180°或者-180°≤θ＜-157.5°，则所述圆内切于下一幅路网子图的右边，切点为该右边的中点；

如果-157.5°≤θ＜-112.5°，则所述圆内切于下一幅路网子图的上边和右边；

如果-112.5°≤θ＜-67.5°，则所述圆内切于下一幅路网子图的上边，切点为该上边的中点；

如果-67.5°≤θ＜-22.5°，则所述圆内切于下一幅路网子图的上边和左边。

Images