一种面向遥感影像的多尺度道路自动提取方法

1.本发明涉及遥感地学应用技术领域，具体为一种面向遥感影像的多尺度道路自动提取方法。


背景技术：

2.随着我国加快实施建设交通强国任务，智能交通在国家战略中扮演重要角色，准确的道路信息是智能交通和智慧城市的基础，遥感技术因其具有多源、覆盖广泛、现势性强等优势，已成为道路监测的重要手段。
3.现有技术的不足：
4.遥感影像道路提取一般分为半自动和自动道路提取方法，半自动道路提取方法主要通过人机交互提取道路信息，工作量大、效率低。道路在高分辨遥感影像上因与其他不同地物相比，有着明显的几何特征、光谱特征和纹理特征；
5.自动道路提取方法在影像道路特征分析和理解的基础上，计算机自动准确识别影像中的道路信息。目前已有的自动道路提取算法，著如基于hough变换和数学形态学算法，受限于遥感影像结构复杂、信息丰富等特点，均得不到理想的道路提取结果。


技术实现要素：

6.本发明的目的在于提供一种面向遥感影像的多尺度道路自动提取方法，以解决上述背景技术中提出的问题。
7.一种面向遥感影像的多尺度道路自动提取方法，本提取方法包括以下步骤：
8.s1、对sentinel-2遥感影像进行大气校正、背景均一化等预处理；
9.s2、采用frangi算法对遥感影像进行多尺度道路系统目标增强；
10.s3、基于遥感影像道路几何特征，多次循环ransac算法提取道路；
11.s4、利用监督分类算法提取遥感影像大尺度道路信息；
12.s5、合并道路提取结果且评定精度，实现道路自动提取。
13.作为本发明的进一步改进，所述步骤s1中：
14.大气校正：对sentinel-2遥感影像近红外波段采用flassh算法进行校正；
15.背景均一化为：使用模糊c均值聚类算法即matlab的fcm函数实现该影像复杂异质背景均一化处理，得到其预处理后图像。
16.作为本发明的进一步改进，所述步骤s2包括：
17.a、构建尺度因子σ下二阶高斯核函数，通过卷积运算建立尺度不变的特征空间，得到预处理后图像对应的hessian矩阵；
18.b、根据对称矩阵的性质，计算出hessian矩阵的两个特征值λ1和λ2及两个特征向量e1和e2，hessian矩阵特征值和对应的特征向量分别代表该点处沿某一方向上图形曲率大小和方向；
19.c、以线性特征最优滤波器响应为依据，采用多尺度融合策略，实现道路系统目标
增强。
20.作为本发明的进一步改进，本方法中步骤a构建尺度因子σ下二阶高斯核函数其中高斯核函数g为：
[0021][0022]
二阶高斯核函数为：
[0023][0024][0025][0026]
作为本发明的进一步改进，所述步骤s1得到的图像由所述步骤s2构建的二阶高斯核函数，使用matlab的imfilter函数的卷积运算来构建尺度不变的特征空间：
[0027][0028]
得hessian矩阵为：
[0029][0030]
其中，为对应尺度因子σ下的二阶偏导图像，为卷积运算。
[0031]
作为本发明的进一步改进，本方法中步骤c遥感影像道路线性特征最优滤波器响应为：
[0032][0033][0034][0035]
其中，l
σ
为对应尺度因子σ下的frangi道路目标增强图像，rb能够用来区分块状要素和目标要素；s用于区分前景线状要素以及背景噪声，β、γ是加权因子控制rb和s的内部响应。
[0036]
作为本发明的进一步改进，本方法中步骤c采用多尺度融合策略，根据尺度空间理论，在空域内按照先验知识确定目标道路的尺度变程范围为σ
min
≤σ≤σ
max
，以最大的道路响应值作为道路目标增强结果：
[0037]
[0038]
作为本发明的进一步改进，所述步骤s3包括：
[0039]
随机选择两个道路预选点构成直线，计算所有的道路预选点到该直线的距离，当满足到直线距离小于2个像素且相互连通的道路像素数大于20时，则保留该结果为道路集；
[0040]
重复以上步骤，多次循环ransac算法直至选出所有的道路。
[0041]
作为本发明的进一步改进，所述步骤s4为：随机选取步骤s3中ransac算法得到的道路提取结果作为道路训练集，则随机选取非道路数据作为非道路训练集，实现训练集数据的自动选取。
[0042]
作为本发明的进一步改进，所述步骤s5合并s3和s4得到的道路提取结果，目视解译地面真值数据由解释人员使用b
red-b
green-b
blue
波段组合手动描绘，将其作为验证样本评定提取道路精度。
[0043]
与现有技术相比，本发明的有益效果是：
[0044]
本发明在对影像进行大气校正、背景均一化等预处理的基础上，采用基于hessian矩阵的frangi滤波算法增强道路系统目标特征，然后利用多次循环ransac算子实现中小尺度道路自动提取，对大尺度的道路，根据道路特征建立道路和非道路训练集，使用最大似然监督分类算法识别较宽道路，最后将中小尺度道路与大尺度道路提取结果合并，生成目标区道路提取结果并评估道路算法提取精度。
附图说明
[0045]
图1为本发明一种面向遥感影像的多尺度道路自动提取方法流程图；
[0046]
图2为上海地区sentinel-2影像近红外波段b
nir
原始影像；(b)为近红外波段b
nir
对应的大气校正及背景均一化等预处理后的影像；
[0047]
图3为多尺度道路系统目标增强，(a)-(c)分别为b
nir
影像对应frangi滤波尺度σ＝1,3,5；(d)b
nir
影像对应frangi滤波尺度σ∈[1,5]；
[0048]
图4为阈值分割道路提取结果；
[0049]
图5为ransac算法多次循环提取的道路结果，(i)为非长条形面状几何特征的道路示例；(ii)为大尺度道路提取缺失示例；
[0050]
图6为ransac算法在图3(i)非长条形面状几何特征的道路中直线构建；
[0051]
图7为最大似然监督分类算法大尺度道路提取结果；
[0052]
图8为道路提取最终结果。
具体实施方式
[0053]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0054]
实施例
[0055]
请参阅图1-8，本发明提供如下技术方案：一种面向遥感影像的多尺度道路自动提取方法：
[0056]
本方法采用的数据是2022年3月9号覆盖上海地区的sentinel-2数据。
[0057]
s1、对sentinel-2遥感影像进行大气校正、背景均一化等预处理：
[0058]
下载目标区域高质量的sentinel-2影像数据，利用envi软件采用flassh算法对该影像做大气校正，并随机选择该影像中的部分区域进行掩膜作为本发明的样本区域，如图2中a区域所示，为sentinel-2近红外波段b
nir
，使用matlab的fcm函数将影像通过模糊c均值聚类算法(fuzzy c-means algorithm)实现复杂异质背景均一化处理，得到预处理后图像，如图2中b区域所示；
[0059]
s2、采用frangi算法对遥感影像进行多尺度道路系统目标增强：
[0060]
在步骤s1的基础上，基于解构局部几何微分结构的核心思想，结合多尺度hessian矩阵分析，构建以尺度空间理论为基础，针对道路特征的多尺度线状目标描述子，通过连续变化的尺度参数来获得目标对象的特征信息，以挖掘影像的道路信息，具体为：
[0061]
a、构建尺度因子σ下二阶高斯核函数其中高斯核函数g为：
[0062][0063]
则二阶高斯核函数为：
[0064][0065][0066][0067]
b、对步骤s1得到的影像，由步骤a构建的二阶高斯核函数，使用matlab的imfilter函数的卷积运算来构建尺度不变的特征空间：
[0068][0069]
即得hessian矩阵为：
[0070][0071]
其中，为对应尺度因子σ下的二阶偏导图像，为卷积运算，高斯尺度空间具有平移、旋转和尺度不变性，并且微分几何属性在不同分辨率重采样条件下能够同一地展示图像结构，frangi算法依据图像的微分几何性质，hessian矩阵能够表示某像元点邻域范围内图像梯度的变化情况，图像的二阶微分几何性质可通过hessian对称矩阵分析。
[0072]
c、利用对称矩阵的性质，计算出hessian矩阵h
σ
(x,y)的两个特征值λ1和λ2及对应的两个特征向量e1和e2，hessian矩阵特征值和对应的特征向量分别代表该点处沿某一方向上图形曲率大小和方向。
[0073]
d、通常在影像中会出现暗背景下的亮特征目标和亮背景下的暗特征目标两种类型，特征值λ1(假定两个特征值中绝对值较小的一项)在两种类型中数值近似为0，而λ2(假定绝对值较大的一项)相对变化较大，其正负值可以区分这两种不同亮度目标，即λ2为正值时道路是暗特征目标，λ2为负值时道路为亮特征目标，如果道路内的某一个像元值满足λ1≈0
且λ2》0，则目标像元所在的道路被认为是暗特征，周围异质物为亮特征。
[0074]
由于道路在sentinel-2影像近红外波段b
nir
上显示为暗特征，故按照以下式计算frangi暗道路线状要素描述子为：
[0075][0076][0077][0078]
其中，l
σ
为对应尺度因子σ下的frangi道路目标增强图像，rb能够用来区分块状要素和目标要素；s用于区分前景线状要素以及背景噪声，β、γ是加权因子控制rb和s的内部响应，线性特征乘以来保证仅在两个标准下的最优滤波器响应。
[0079]
e、采用多尺度融合策略，根据尺度空间理论，在空域内按照先验知识确定目标道路的尺度变程范围为σ
min
≤σ≤σ
max
，设置σ
min
＝1和σ
max
＝5，尺度因子步长为2，其余参数分别设置为：描述子采用的窗口大小w＝6scale，β＝0.5，γ＝15，按照步骤a至d遍历各个尺度因子σ进行目标增强，将探测尺度空间的极值响应记录到最优尺度匹配线性特征指标中，线性特征宽度满足最优尺度因子σ，l
σ
就会达到最大值。在计算一系列尺度时，取最大的响应值尺度作为最终道路尺度，即：
[0080][0081]
不同尺度下道路提取结果如图3所示，采用阈值分割(灰度值》0.1)得到道路预选点提取二值化结果，如图4所示；
[0082]
s3、基于遥感影像道路几何特征，多次循环ransac算法提取道路：
[0083]
ransac算法是根据一组包含离群的被观测数据中估算出数据的数学模型参数，得到有效样本数据的算法；
[0084]
依据高分辨率遥感中道路具有长条形面状几何特征，在第二步提取道路点的基础上进行ransac算法运算：随机选择两个道路预选点构成直线，计算出这两个点所表示的直线方程：y＝ax b，并计算所有的道路预选点到该直线的距离，当满足到直线距离小于2个像素且相互连通的道路像素数大于20时，则保留该结果为道路集；
[0085]
重复以上步骤，多次循环ransac算法(在此设置循环次数为5000)直至选出所有的道路，如图5所示，需要提到的是，尽管除大部分长条形面状几何特征的道路外，还存在一些曲线道路，如附图5(i)所示，对于这种情况也可以通过以上ransac算法构建直线方程得到，附图6所示；
[0086]
s4、利用监督分类算法提取遥感影像大尺度道路信息：
[0087]
高分辨率遥感影像中存在大尺度的道路，如附图5(ii)所示，可以通过监督分类算法提取得到，从而缩小尺度空间的搜索范围；
[0088]
首先按照b
nir-b
green-b
blue
波段组合进行rgb假彩色图像合成，然后依据第三步道路信息制作道路和非道路两类文件作为训练集，通过训练集数据计算每一个类别的判别函数，经与k-近邻算法、决策树算法等其他监督分类算法对比，最大似然算法得到道路解译结果最优，其判别函数为：
[0089][0090]
其中，(n维列向量)
[0091][0092][0093][0094]
按照上式分别构建道路和非道路的判别函数，然后将每一个像素点的值带入进行计算，计算结果大的函数所对应的类别属于该像素点的类别，直至分类结束；
[0095]
对以上分类结果，可以使用matlab的函数bwareaopen(参数设置为500)，去除小斑块，最终实现影像大尺度道路提取结果，如附图7所述；
[0096]
s5、合并道路提取结果且评定精度，实现道路自动提取：
[0097]
合并步骤s3和s4得到的道路提取结果，如图8所示，目视解译地面真值数据由经验丰富的解释人员使用b
red-b
green-b
blue
波段组合手动描绘，将其作为验证样本评定提取道路精度，并参考表1计算kappa系数(kappa coefficient，k)、总体分类精度(overall accuracy，oa)、错分误差(commission error，ce)和漏分误差(omission error，oe)作为精度评估指标。
[0098]
表1定量评估系数及计算结果
[0099][0100]
其中p0和pc分别为和(tp)为正的真值，(tn)为负的真值，(fn)为假的负值和(fp)为假的正值；
[0101]
经精度评估后，得到的结果为表1所示，kappa系数和总体精度分别高达96％和98％，且错分误差和漏分误差则较低，分别为3％和5％，均达到了满意精度。
[0102]
本发明在对影像进行大气校正、背景均一化等预处理的基础上，采用基于hessian矩阵的frangi滤波算法增强道路系统目标特征，然后利用多次循环ransac算子实现中小尺度道路自动提取，对大尺度的道路，根据道路特征建立道路和非道路训练集，使用最大似然监督分类算法识别较宽道路。最后将中小尺度道路与大尺度道路合并，生成目标区道路提取结果并评估道路算法提取精度。
[0103]
需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性地包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。
[0104]
尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。