技术领域本发明属于遥感影像处理与目标提取技术领域，尤其涉及一种基于OSM(OpenStreetMap，众源地理数据)的高分辨遥感影像道路变化检测方法。
背景技术：
道路是高分辨遥感影像中的主要目标，道路变化检测技术对GIS数据更新、目标检测和识别、影像匹配意义重大。国内外许多学者对道路变化检测技术进行了研究，特别是近几年，随着我国经济快速发展，我国道路建设取得长足发展。同时，我国城乡一体化建设使得许多区域发生了巨大变化，道路信息的快速变化对现有的地图数据更新速度提出了更高的要求。如何快速检测道路重要变化并实时更新，对我国的基础地理数据更新有着重要的意义。另一方面，由于近几年遥感技术的快速发展，遥感数据数量庞大。随着遥感卫星性能提高，高分辨遥感数据的获取已经变得容易，但是现有的高分辨道路影像的提取和变化检测技术、方法和理论远远不能满足现有的图像数据处理的需要，基于这种情况，加强高分辨遥感道路变化检测的研究势在必行。实验证明，遥感影像和矢量数据结合进行道路变化检测可以提高检测速度和精度。但是，由于矢量数据的更新速度滞后于遥感影像的获取速度，所以如何快速地结合矢量数据和遥感数据进行道路检测成为重要课题。而且，在道路变化检测过程中，矢量数据通常作为旧数据和新遥感影像作对比来检测道路变化情况。OpenStreetMap，又称为“面向众源开放街道地图”，简称为OSM，具有现势性高、数据丰富、获取方便等特点。将其与高分辨遥感影像结合，可以作为一种获取道路变化检测结果的重要技术和手段。变化检测技术是对不同时段的目标或现象状态发生的变化进行识别和分析，从而判断目标是否发生变化，并确定变化的空间分布和变化类别。由于高分辨遥感影像中道路信息的复杂性以及干扰性较多，采用常规的变化检测方法(如图像差值法、图像比值法)容易造成误判，这将无法准确获取道路变化信息。
技术实现要素：
针对现有技术存在的问题，本发明结合高分辨遥感影像和OSM数据，提供了一种高效准确的、基于OSM的高分辨遥感影像道路变化检测方法。本发明充分结合遥感影像数据和OSM矢量数据信息，利用OSM矢量数据数据量大、现势性好、信息丰富、成本低等特点和优势，将提取的矢量化道路数据和OSM矢量数据比对，以检测道路是否变化。为解决上述技术问题，本发明采用如下的技术方案：基于OSM的高分辨遥感影像道路变化检测方法，包括：采用长度完整性和定位精度评价OSM数据质量，将长度完整性大于第一阈值且定位精度大于第二阈值的OSM数据用于道路变化检测；第一阈值和第二阈值根据经验或采用专家投票法确定，长度完整性QL＝LOSM/LR，定位精度LOSM和LR分别为OSM道路和正规参考地图中道路的总长度，为缓冲区内OSM道路总长度；所述的道路变化检测包括步骤：S1将OSM数据和高分辨遥感影像进行配准；S2提取高分辨遥感影像的矢量道路，即道路中心线，本步骤依次包括子步骤：2.1对高分辨遥感影像进行傅立叶变换；2.2采用多方向的Log-Gabor小波滤波器对高分辨遥感影像像素点进行滤波；2.3基于相位一致性原理，获得高分辨遥感影像多方向的相位一致性灰度图和边缘信息；2.4基于形态学对边缘信息进行滤波，得道路中心线，具体为：人为选择一方向的相位一致性灰度图进行二值化，得二值化灰度图；分割二值化灰度图得分割图像；对分割图像进行连通区域分析，计算各连通区域的面积指数SArea和狭长指数SLFI，并保留满足条件SArea≥TAreaSLFI≥TLFI]]>的连通区域；TArea、TLFI分别为根据经验设定的面积指数阈值和狭长指数阈值；S3采用仿射变换纠正道路中心线，并通过叠加纠正后道路中心线和OSM道路的缓冲区识别道路变化；上述，道路中心线和OSM道路的缓冲区均按照道路实际宽度生成。作为优选，子步骤2.2中，Log-Gabor小波滤波器的方向数为6。作为优选，子步骤2.3中边缘信息可采用如下方法获得：对高分辨遥感影像中所有像素点，分别计算其在各方向的局部能量，各方向局部能量之和即该像素点的局部能量，局部能量峰值对应的像素点构成边缘信息。子步骤2.4中，面积指数SArea即连通区域内像素点数。子步骤2.4中，狭长指数其中，LMER和WMER分别表示连通区域最小外接矩形的长和宽；np表示连通区域面积。步骤S4中，采用加色法区识别道路变化。和现有技术相比，本发明具有如下优点和有益效果：(1)目前，大多数基于高分辨遥感影像的矢量道路均是基于亮度和梯度进行提取，提取质量很大程度上依赖原始影像的成像质量。而本发明基于高分辨遥感影像的矢量道路提取则不受原始影像亮度、梯度和对比度等成像质量的影响，且可省去遥感影像预处理过程。(2)其他的基于高分辨遥感影像的矢量道路提取易产生双边缘问题，本发明基于高分辨遥感影像的矢量道路提取不存在双边缘问题，可获得更完整准确的道路信息。(3)OSM数据来源丰富，现势性好，易获取且成本低廉，但准确性缺乏保证。本发明采用长度完整性和定位精度两个质量要素评价OSM数据质量，采用质量符合要求的OSM数据进行道路变化检测。(4)采用加色法显示道路变化，简洁明了。附图说明图1是本发明的具体流程图；图2是实施例中分辨率为1.5米的原始高分辨遥感影像；图3是图2的相位一致性灰度图；图4是图3的二值化灰度图；图5是OSM数据与遥感数据的配准；图6是纠正后的道路中心线图；图7是实施例的道路变化检测结果。具体实施方式具体实施方式中，根据相位一致性原理(phasecongruency，PC)，采用Log-Gabor滤波器，提取高分辨遥感影像中道路的边缘特征，矢量化获取道路中心线；结合OSM数据，通过仿射变换进行叠加纠正后，分别按照道路实际宽度做缓冲区，使用加色法进行道路变化检测。图1为本发明的具体流程图，下面将结合图1和具体实施方式进一步说明本发明步骤。步骤1，OSM数据的质量评价。OSM数据是指由大众采集并向大众提供的开放地理空间数据。与传统地理信息采集和更新方式相比，来自非专业大众的OSM数据具有数据量大、现势性好、信息丰富、成本低等特点和优势。本发明采用长度完整性和定位精度两个质量要素评价OSM数据的质量，将质量符合要求的OSM数据用于后续的道路变化检测。(1)长度完整性长度完整性是表征道路数据的覆盖程度和几何质量的特征，该特征对道路变化检测有着重要作用，只有OSM数据在长度完整性上达到一定标准，才能作为检测道路变化的参考数据，从而减少误差。长度完整性QL定义为OSM数据中道路(后文简记为“OSM道路”)总长度LOSM占正规参考地图中道路总长度LR的比例，见公式(1)：QL＝LOSM/LR(1)长度完整性QL是OSM数据覆盖情况最直接、最直观的表现，反映了OSM数据的可用性，是最为重要的质量评价要素之一。(2)定位精度定位精度是用来评价OSM数据几何精度的指标，和长度完整性一样，定位精度是评价OSM可用性的另一重要质量要素。本发明采用缓冲区分析法评价OSM数据的定位精度。各地区的道路都有一定的宽度标准，进行缓冲区分析时，首先按照实际道路的宽度标准做缓冲区，然后在缓冲区覆盖下评价OSM数据的定位精度QLP，见公式(2)：QLP=LOSMP/LOSM---(2)]]>式(2)中，LOSM为OSM道路总长度，为落在缓冲区中OSM道路总长度。当长度完整性大于预设的第一阈值且定位精度大于预设的第二阈值时，OSM数据才可用于后续的道路变化检测中。第一阈值和第二阈值可根据经验设定，也可通过多次试验结果调整第一阈值和第二阈值。第一阈值和第二阈值也可根据专家投票法设定。步骤2，OSM数据和高分辨遥感影像的配准。OSM数据和高分辨遥感影像的配准为本领域的常规技术。具体实施时，可将OSM数据和高分辨遥感影像导入Arcgis软件，使用工具EditorforOSM建立地理数据库，即可对导入的OSM数据和高分辨遥感影像进行自动配准。步骤3，基于相位一致性原理，提取高分辨遥感影像的矢量道路。本步骤依次包括以下子步骤：3.1高分辨遥感影像的傅立叶变换。本子步骤为本领域内常规技术，在此不作赘述。3.2采用多方向的Log-Gabor小波滤波器对高分辨遥感影像像素点进行滤波。由于Log-Gabor小波函数能真实地反映自然影像的频率响应，所以本具体实施方式中滤波器选择Log-Gabor小波函数。滤波器组的设计涉及到滤波器带宽与相邻滤波器间中心频率比例因子间的关系，目标是用较少滤波器得到宽度合理、统一覆盖的谱。从Gabor滤波器可得到的最大带宽约为1倍频，要得到统一的谱覆盖，相邻滤波器中心频率之间的比值不能大于1.5。这样，构造跨越4倍频的滤波器组就需要8个滤波器。Log-Gabor小波滤波器的大带宽可以为设计滤波器组提供更大的灵活性。例如，要得到4倍频的滤波器组，可以不用8个1倍频的滤波器，只用4个相邻中心频率之比为2.6的2倍频滤波器(可以用中心频率之比为3的滤波器，而且仍有平坦的谱覆盖)。使用较少的滤波器意味着减少计算量；且，Log-Gabor小波函数形状中的尖峰可以为空间域提供精确定位。3.3基于相位一致性原理，根据滤波后的高分辨遥感影像获得多方向的相位一致性灰度图和边缘信息。相位一致性的检测原理是：通过Fourier变换，在影像频率域中寻找不受亮度、对比度影响的像素点，即边缘像素点。相位一致性的检测过程可以完整检测灰度亮度不明显的图像特征，在待处理遥感图像未进行预处理时，也能保持提取的边缘信息的精度及连贯性。按Fourier级数展开高分辨遥感影像像素点x的灰度值I(x)，见式(3)：I(x)=ΣnAncos(nωx+φn)=ΣnAncos(φn(x))---(3)]]>像素点x的相位一致性值PC(x)见式(4)：式(3)～(4)中，n为Fourier谐波次数，由谐波频率同基波频率之比确定；An表示n次谐波分量的幅值；ω为常数，一般取2π；φn为n次谐波分量的相位偏移量；函数表示像素点x的Fourier分量的局部相位；为像素点x处所有Fourier项平均局部相位。相位一致性值PC(x)的极值点即局部相位变化最小的点，记为最大相位一致性点。直接利用式(4)计算相位一致性值比较复杂，为简化计算，本发明通过寻找像素点灰度值的局部能量峰值获得最大相位一致性点。灰度值的局部能量是按照像素点本身以及它的Hilbert变换进行定义的，局部能量E(x)定义如下：E(x)=F2(x)+H2(x)---(5)]]>式(5)中：F(x)表示像素的直流分量；H(x)表示F(x)的Hilbert变换，即由F(x)经相位移动90°得来。局部能量和相位一致性值PC(x)存在如下关系：E(x)=PC(x)ΣnAn---(6)]]>因此，局部能量峰值对应的像素点即最大相位一致性点。本具体实施中，采用式(5)计算像素点在各方向的局部能量E(x)，各方向的局部能量E(x)之和即该像素点的局部能量；局部能量峰值对应的像素点构成边缘信息。3.4对边缘信息进行形态学处理，得矢量道路，即道路中心线。从多方向的相位一致性灰度图中人为选择特征清晰、对比度强的一幅影像进行二值化，得二值化灰度图，利用形态学对二值化灰度图进行滤波，得到去除了非道路噪声信息的道路信息。道路在高空间分辨率遥感影像上具有以下特征：在一定方向是连续的；具有一定的曲率，除了环形道路外，曲率变化很小；道路宽度在一个局部区域内是一个定值；道路和背景的亮度比较显著，要么偏暗，要么明亮。本子步骤进一步包括：(1)人为选择一方向的相位一致性灰度图进行二值化，得二值化灰度图。二值化阈值根据经验设定，要保证既能检测出所有道路，又不能包含过多噪声。(2)分割二值化灰度图，并对分割图像进行连通区域分析(connectedcomponentsanalysis，CCA)；(3)计算各连通区域的面积指数SArea和狭长指数SLFI；(4)保留分割图像中满足公式(7)的连通区域，并剔除不满足公式(7)的连通区域。SArea≥TAreaSLFI≥TLFI---(7)]]>式(7)中，TArea、TLFI分别为根据经验设定的面积指数阈值和狭长指数阈值，这两个阈值可根据实际精度要求灵活设置。下面将提供具体实施方式中面积指数和狭长度指数的计算方法。(a)面积指数道路网在高分空间辨率遥感影像上呈现的是连续的面状区域目标，道路的连续面积不会很小，所以可使用面积指数区分道路与小面积噪声。连通区域的面积指数即该连通区域内像素点数。(b)狭长指数根据道路特征描述可知，道路是连续的狭长区域，所以可以利用线性度量指标长宽比指数来去掉非狭长区域，长宽比指数LFI即狭长指数：LFI=LMERWMER---(8)]]>式(8)中：LMER和WMER分别表示MER区域的长和宽，MER区域即连通区域的最小外接矩形。根据公式(9)定义新矩形，新矩形的长宽比指数即连通区域的狭长指数SLFI。LW=npL2=LMER2+WMER2---(9)]]>式(9)中，W和L分别代表新矩形的宽度和长度；np即新矩形面积，连通区域的狭长指数SLFI见式(10)：SLFI=LW=Lnp/L=L2np---(10)]]>步骤4，道路变化检测。4.1道路中心线的纠正。步骤3提取的矢量道路即道路中心线，而OSM数据则是用户的轨迹折线，道路中心线和OSM数据会出现不吻合现象，从而会对加色法生成缓冲区造成一定影响，所以需要对OSM数据和道路中心线进行纠正。本具体实施中采用仿射变换对道路中心线进行纠正。4.2叠加道路中心线和OSM道路的缓冲区识别道路变化。将纠正后的道路中心线与OSM道路叠加。为使结果更为直观，本具体实施方式采用加色法进行道路变化识别。将提取的道路中心线缓冲区设为绿色，OSM道路缓冲区设为红色。叠加道路中心线和OSM道路，通过颜色即可直观识别道路变化情况。实施例本实施例中采用的原始高分辨遥感影像见图2，分辨率为1.5米，OSM数据从网络获取。(1)OSM数据的质量评价。研究区域为辽宁省大连市城区，总面积约12574平方公里，包括西岗区、中山区、沙河口区、甘井子区、旅顺口区等4个行政区域，共1471条道路，道路总长约为3466km。本实施例中OSM数据来源于OpenStreetMap网站，为WGS-84经纬度坐标系，正规参考地图为北京四维图新2011年生产的GPS版导航地图，定位精度约4m。本实施例OSM数据的质量评价结果见表1。采取专家投票法，认为长度完整性和定位精度超过60％即可作为道路变化检测的参考数据。由表1可看出，OSM数据的长度完整性和定位精度都比较理想。表1OSM数据的质量评价结果(2)将高分辨遥感影像和OSM数据进行配准。高分辨遥感影像和OSM数据的配准是本领域内常规技术，可采用Arcgis软件工具进行。高分辨遥感影像和OSM数据的配准见图5，其中灰色部分为缓冲区。(3)提取高分辨遥感影像的矢量道路。本步骤进一步包括：3.1对高分辨遥感影像进行傅立叶变换。3.2对傅立叶变换后的高分辨遥感影像像素点进行滤波。计算相位一致性值前需要对高分辨遥感影像各像素点进行滤波，计算量较大，所以滤波器方向不宜选用过多。结合遥感影像上道路分布特性，本实施例选择6个方向，即0°、30°、60°、90°、120°、150°方向的Log-Gabor小波滤波器提取道路特征。方向数并不限于6。实验表明滤波器方向间隔为30°时，即可得到平坦的频谱覆盖，又能使用最少的方向数。使用更多的方向数并不会明显改变结果性能，方向过多计算量也将增大。在滤波器的尺度和频率选择上，选用4个频率尺度，其中最小尺度波长为3，尺度倍数为2.1，滤波器的方向角度标准差为1.5，频率扩散的阈值为0.4，噪声能量阈值为2.0，用于计算频率扩散相位一致性的调和函数的收敛权值为10。采用Log-Gabor小波滤波器对图2遥感影像进行滤波处理。Log-Gabor小波滤波器可选用Matlab语言编程实现。3.3对滤波后高分辨遥感影像进行相位一致性计算，得相位一致性灰度图和边缘信息。具体实施时，对滤波后高分辨遥感影像的所有像素点，分别计算其6个方向的局部能量，各方向的局部能量之和即当前像素点的局部能量，局部能量峰值对应的像素点构成边缘信息。图2的相位一致性灰度图见图3。从图3可以看出，道路不是单像素宽，是脉冲状边缘，其一阶导数有两个极大值，这样边缘检测时，会出现双边缘响应，对道路检测效果影响较大。3.4对相位一致性灰度图进行二值化，得二值化灰度图。本实施例中二值化阈值设为0.9。当像素点6个方向的局部能量之和大于0.9，该像素点的灰度值设为1，否则为0，图4是图3的二值化灰度图，其中白色区域代表道路特征。分割二值化灰度图，并采用形态学方法对分割图像进行滤波，获得道路中心线。(4)道路变化识别。提取道路中心线缓冲区纠正图，见图6。道路变化识别结果见图7，图中，椭圆圈出道路为错检道路，即本实施例中错检道路数位2条。表2本实施例道路变换检测漏检精度表总道路数错检数错检率3026.7％