一种面向高分辨率遥感影像的植被环岛自动化提取方法与流程

本发明属于视觉处理
技术领域：
，涉及一种从高分辨率多波段遥感影像中自动化提取植被环岛的方法。
背景技术：
：环岛是一种常见于交通压力较小且位置充足的地段的交通设施，由于环岛将车辆的汇合点转变成行驶点，其能够减缓车辆行驶速度、减少碰撞事故的发生，因而改善交通质量。环岛的半径一般在12至30米，随道路等级不同而变化，有一些环岛甚至能够成为地标或休闲广场。环岛中常种有植被，这些植被一般为常青植物、花朵或低矮灌木，应当有明显的边缘和不影响驾驶员驾驶的高度。对于植被环岛的提取可以了解交通状况、帮助市政规划，完善辅助人工或自动驾驶的电子地图，或衡量一片区域的繁荣程度，对于城市的交通、设计和人文研究都具有一定意义。对于从遥感影像进行环岛提取的研究较少，在2009年ravanbakhsh提出了地理空间数据库作为先验知识的环岛提取方法，从而发展了环岛提取的研究。然而，目前的环岛提取方法仍有很多不足之处，如提取准确率低，没有很好地利用遥感影像中的植被信息，或必须依赖其他辅助数据等。因此，提出一种从高分辨率多波段遥感影像中自动化提取植被环岛的方法是本领域的一项待攻克的难题。技术实现要素：为了解决上述技术问题，本发明提供了一种面向高分辨率遥感影像的植被环岛自动化提取方法。本发明的技术方案为：一种从高分辨率多波段遥感影像中自动化提取植被环岛的算法，包含以下步骤：步骤1：并行执行流程a和流程b；流程a：对原始影像进行影像分割，提取环岛和环岛块；流程b：对原始影像进行道路分类，提取道路空洞和空洞块；步骤2：检验环岛块和空洞的拓扑关系；步骤3：获取环岛提取结果。本发明中，图像分割有两种分水岭算法和两种融合模式可供选择，支持向量机分类有四种核函数可供选择，比较灵活；本发明对于处理的影像没有较高的要求，能够适应不同大小、不同剪裁、不同分辨率、不同波段配置的高分辨率多波段的遥感影像，具有很强的适应能力和鲁棒性；本发明的提取结果具有较高的准确率，和人工目视解译的结果有较高的重合率；本发明的数据处理时间较短，能够在短时间内处理大量的环岛提取；本发明充分利用了影像本身的信息，充分考虑了地物的植被信息，并没有借助其他辅助知识或数据，有十分强的独立性；本发明在数据情况不理想或算法本身导致的误差的情况下，仍能进行粗略的拓扑关系判断，具有一定的容错率，保证了提取结果的可靠性。本发明改良和创新了从遥感影像中提取环岛地物的算法，为环岛的自动识别和提取任务降低了成本、增加了稳定性和准确率，对遥感地物提取的相关研究具有积极的推动作用。附图说明图1为本发明实施例的流程简图；图2为本发明实施例的预处理流程简图；图3为本发明实施例的面向对象的图像分割流程简图；图4为本发明实施例的基于规则的提取流程简图；图5为本发明实施例的基于支持向量机的道路分类的流程简图；图6为本发明实施例的空洞检测的流程简图；具体实施方式为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。本发明是一种从高分辨率多波段遥感影像中自动化提取植被环岛的算法，包含以下四步：面向对象的图像分割和基于规则的提取，基于支持向量机的道路分类，道路空洞提取，检验环岛块和空洞的拓扑关系。在第一步之前，可以有适当的图像预处理；在第四步后，可以将提取结果与人工目视解译的结果进行对比，以对自动化算法的成果进行评价。第一步可以与第二步和第三步并行，第二步必须在第三步前，第四步必须在第一、二、三步完成后再进行，如图1所示。预处理，是指对原始的遥感影像进行几何校正、辐射定标、天地图配准、裁剪等处理，使得尽可能提高原始遥感影像的质量，希望从数据源上提高本发明的成果的准确率。顺序上，几何校正和辐射定标的没有指定的先后顺序，天地图配准必须待前两者完成后才可进行，剪裁必须待前三者完成后才可进行，如图2所示。几何校正是指消除或减少原始影像的几何变形；辐射定标是指将传感器记录的无量纲的dn值转换成具有实际物理意义的大气顶层辐射亮度或反射率；天地图配准是将几何校正、辐射定标的影像与标准的、正射的世界地图进行配准；裁剪是为了减少处理时间，将大幅的遥感影像裁为小幅的包含一个或多个完整环岛的过程。整个预处理并非必须，但适当的预处理能够在一定范围内保证本发明提取结果的准确率，并能够缩减数据的处理时间。实施例具体实施方案为：几何校正：在envi5.3sp软件中，打开高分二号原始影像数据及其元数据(其他高分辨率的多波段影像亦可)，选择geometriccorrection下orthorectification下的rpcorthorectificationworkflow工具，选择输出像素大小如“4米”，选择重采样方式为“三次项方程”，完成输出路径等设置后即可按设定参数实施几何校正。辐射定标：在envi5.3sp软件中，打开高分二号原始影像数据或几何校正后的数据，选择radiometriccorrection下的radiometriccalibration工具，定标类型设为辐射率数据radiance，由envi软件自动设置flaash大气校正工具需要的数据类型，存储顺序为bil或者bip，数据类型为float，设置辐射率数据单位调整系数为0.1，完成相关设置后即可按设定参数实施辐射定标。天地图配准，可以使用arcgis10.2软件中的arcmap10.2下，调用其georeferencing功能，针对正射的标准世界地图的影像拍摄到的区域，选择几对明显的人工地物通过点击(标准世界地图http://www.scgis.net.cn/imap/imapserver/defaultrest/services/newtianditudom/wms)进行配准。裁剪，可以使用envi5.3sp软件中的resize功能，指定裁剪后角点位置和剪裁大小后(如左上角坐标114.468,30.471，大小为400*400像素)，即可完成裁剪。本发明的第一个正式步骤：面向对象的图像分割，使用了分水岭算法针对整幅高分辨率多波段的遥感影像进行分割，分割的结果为影像块，针对众多影像块，将进入基于规则的提取进行处理。分水岭算法起源于水文上的概念，即随着水位的上升，不同高程区域的水体将会汇集，分水岭就建在来自不同高程区域的水体的汇集处。在数字图像处理领域，分水岭算法是常用的图像分割算法，理想的分割结果，应当与现实世界中的对象一致，即理想情况下，一个分割块应该对应一个现实中的实体，这是面向对象思想的体现。分水岭分割包括分割和融合两步，如图3所示。分割有两种融合模式可选：边缘或强度，边缘模式是指对原图进行sobel算子的边缘检测，再形成梯度图，进行分水岭算法的分割，强度模式是指针对原图，从中选择部分波段，取一定范围内的均值生成强度图，再进行分水岭算法的分割。融合是为了避免过度分割的问题，判断分割后的各影像块根据块间的相关程度后，进行融合。融合共有两种融合模式可选：全λ模式或快速λ模式；全λ模式公式如下，若两像素的合并代价ti,j小于阈值，则合并：其中：oi是影像中的区域i；|oi|是区域i的面积；μi是区域i的均值；||μi-μj||是区域i和区域j光谱值的欧氏距离；是oi和oj公共边界的长度；快速λ模式公式如下，lambda值越小，说明欧式颜色距离越小，公共边界长度越大，越应当合并，若lambda值小于阈值，则合并；其中：n1是区域1中的像素数；e是区域1和区域2的欧式颜色距离；l是区域1和区域2的公共边界长度。第一步的基于规则的提取，是对分割结果同时考虑了植被环岛的明显特征：大小、圆度和植被指数，将三者的综合考虑设置为规则，即：大小在影像所属国家和地区的交通部门指定的环岛半径为半径的圆面积为上下限的范围内；圆度(公式如下)在接近1的一定范围内；植被指数使用ndvi(normalizeddifferentialvegetationindex，归一化植被指数，公式如下)并且植被环岛所在的影像块的ndvi值应当不低于整幅影像ndvi值的前40％，如图4所示。规则的三种属性设置完毕后，参照规则对所有影像块进行筛选；筛选后的留下的影像块应当有且仅有一个，即环岛块。其中：r是区域的圆度；s是区域的面积；c是区域的周长；其中：ndvi是归一化植被指数；eps是一个足够小的且能避免分母不为0的常数；b1是红外波段；b2是近红外波段；实施例具体实施方案为：对预处理后的影像，选择envi中的rulebasedfeatureextractionworkflow工具，勾选normalizeddifference选项，并将影像的第三波段(红波段)设为band1，第四波段(红外波段)设为band2，以此为设置参与ndvi的计算。接下来进入图像分割，需要分别设置图像分割的分水岭算法模式和参数，以及融合的模式和参数，可以设为：“边缘模式，程度50，全λ模式，0”。接下来设定规则。由于本发明中的规则包括大小、圆度、植被指数三项，则先创建规则，再点击创建好的规则，创建三个属性，属性可以分别如下设置：“空间属性-大小：452-2826；空间属性-圆度：0.5-1.2；波段属性-归一化指数(ndvi)：60％-100％”。按照规则运行，进行筛选，查看筛选后的结果，主要关注两个方面：环岛在基于规则的筛选后是否留下了；筛选后是否有多余的非环岛的块。如果前者为是，后者为否，则完成了面向对象的图像分割和基于规则的提取。可将提取结果转为evf再转为shp文件，以便后续操作。值得说明的是，esri的shapefile数据格式(.shp文件)并非是本发明的唯一可行的数据格式，视具体操作软件和生产需求而不同。基于支持向量机的道路分类，是指在选取适当数量的道路样本点后，使用支持向量机算法(supportvectormachine,svm)进行影像分类。svm是模式识别常用的算法，其思想为：对于线性不可分的情况，通过使用非线性映射算法将低维输入空间线性不可分的样本，转化为高维特征空间(超平面)使其线性可分，从而使得高维特征空间采用线性算法对样本的非线性特征进行线性分析成为可能。道路样本点可以在envi软件中以人工方式，鼠标选中多边形的方式形成roi(regionsofinterest)获得道路像素的样本，也可以通过其他方式指定分类样本。分类的核函数有四种备选，分别为线性、多项式、径向基函数、s形(公式如下)，输出为包含带状实体的分类结果。建议如果无特殊要求时，可以选择径向基函数为核函数，因为其在大多数情况下效果较好。此外还要设置svm算法要求的三个参数：金字塔层数，惩罚参数，分类概率阈值(如图5所示)。金字塔层数是在支持向量机中应用的分级处理级别的数量，以避免分类或过度分类。乘法参数允许一定程度的误判，理论上允许一些训练点在错误的一边的超平面，从而使得支持向量机模型具有较强的容错性和灵活性。分类概率阈值表示在超平面中，最接近的部分点代表该类的置信度，因此较高的阈值会导致较少的最近点被分类和图像中的未分类部分的增加。简单地说，分类概率阈值可以控制分类的程度。在本发明中，一个令人满意的基于svm的道路分类结果是连续的带状的实体，即带状的道路，其中应当有一个接近圆形的大洞。由于这个大洞应当是由于环岛导致的，因此这个洞和环岛的空间位置应该基本重合。细微的误差难以避免，如车辆或行道树阴影导致的道路提取结果的细碎破洞和非平整边缘，以及洞和环岛的边缘存在一些空隙或重合的现象。本发明对这些难以避免的误差皆能够容忍，在之后的两步有所体现。线性：k(xi,xj)＝xitxj多项式：k(xi,xj)＝(gxitxj+r)d径向基函数：k(xi,xj)＝exp(-g||xi-xj||2)s形：k(xi,xj)＝tanh(gxitxj+r)其中：g是在除线性核函数外其他函数所需的gamma参数；d是多项式核函数所需的参数；r是多项式和s形核函数所需要的偏项参数。实施例具体实施方案为：打开待分类的影像，在envi中选择classification下的supervised的supportvectormachine工具，进行svm分类。可以进行如下设置：选择径向基函数(radialbasisfunction)，设置gamma为0.25,惩罚参数(penaltyparameter)为100，金字塔层数(pyramidlevel)为0，分类概率阈值(classificationprobabilitythreshold)为0.7。需要说明的是，这些设置的分类效果随不同的影像而不同，有时可能需要多次适当的调整才能达到理想的效果。道路空洞提取，只提取带状道路中的空洞，将最大的空洞视为环岛造成的空洞。这是因为正常情况下，道路的基于svm的分类效果难以产生半径超过12米的提取误差，即在面积上大于环岛的误差。空洞提取包括四步，每一步皆用gis(地理信息系统)工具完成，其中三步是转换；第一步，转换多边形为点，即转换带状道路多边形为中心点，中心点属性包括编号fid、形状shape、类型_名称class_name、类型_编号class_id、区块数parts、长度length、面积area、源编号orig_fid；第二步，转换多边形为线，即转换带状道路多边形为边界线，这些边界线只有少量属性，边界线属性包括编号fid、形状shape、类型_名称class_name、类型_编号class_id、区块数parts、长度length、面积area、源要素id；第三步，转换线为多边形，即根据边界线再生成多边形，多边形属性包括编号fid、形状shape、源要素id、类型_编号class_id、区块数parts、长度length、面积area、源编号orig_fid；根据每一多边形的id或其他标志，将第一步产生的点的属性加到第三步生成的多边形中，若源多边形不存在，则附加属性为空；第四步，类型_编号class_id、区块数parts、长度length、面积area、源要素id皆为0的多边形为空洞；最后将空洞按面积排序，选出最大的，即为环岛导致的道路空洞，称为空洞块(流程如图6所示)。也就是说，道路的空洞提取利用了在gis多种转换中，属性的遗失和保留现象，筛选出了空洞；由于一开始时空洞并不为一个单独的多边形实体，因此在三步转换后，由空洞的边界生成了和空洞形状一致的多边形，和所有转换前的多边形相比，自然存在属性上完整度和项目数的不同。实施例具体实施方案为：在arcmap10.2中，可逐次使用“featuretopoints”,“polygontolines”,“featuretopolygons”,完成转换多边形为点，转换多边形为线，转换线为多边形。需要强调的是，在第三步转换操作时，用第一步转成的点做labelfeatures，并选中preserveattributes选项，才能对第三步产生的多边形添加上第一步获取的属性，进而找到所有空洞。最后，选中所有第三步产生的多边形shapefile，利用基于属性的查询找到那些遗失了例如classname,area和originalfid属性的多边形，再点击排序按钮，按大小排序。选中面积最大的多边形，保存为shp格式，称为空洞块。检验环岛块和空洞的拓扑关系，目的是希望容忍了由于众多误差导致的环岛块和空洞的不完全重合现象，允许了环岛块和空洞边界的空隙或重合，即包容了道路和环岛块的关系并不为理想状态下的包含关系的现象，但与此同时，仍然粗略地检验了环岛块和道路的拓扑关系，为判断环岛块是否为环岛提供依据。具体做法是，通过分别计算环岛块和空洞的中心点，再测量两中心点的距离，将此距离与一定阈值进行比较，若距离小于阈值，则判定提取的环岛块是环岛，若距离大于阈值，则提取无效。实施例具体实施方案为：在arcmap10.2中，对前述步骤提取出的环岛块和空洞块shapefile文件使用“featuretopoints”,获得两者的重心，或称中心。重心的计算是取该实体所有边界点的二维坐标的均值。将两个中心以不同的方式显示，如红色三角与黑色圆点。再用arcmap的measure工具测量两个点的距离。参考遥感影像的尺寸和分辨率，结合道路的一般宽度，可以设置一个粗略的阈值。如原始遥感影像为7300*6908像素，裁剪后的遥感影像为400*400像素，每一像素代表4*4米的地面实际面积(地面分辨率为4m*4m)，待判断的环岛所在的地区为中国湖北省黄石，道路宽度约8-30米，综合这些因素，可设置阈值为5米。将距离与阈值(5m)相比，做出相应的判断，粗略地判断环岛块和道路空洞的大致拓扑关系是否为重合，即提取的环岛块是否可以认作是真正的环岛。阈值的具体数值没有太多的局限，只要能够体现出本发明的从两点距离判断环岛提取效果的特点，并能够大致评判环岛的提取效果，即可。测试本发明的环岛提取效果的最好的方法就是将多次多个提取出的环岛块与空洞块的距离进行统计和综合分析，以及将多个提取出的环岛块和人工目视解译的效果进行对比。表1所示显示了在中国湖北省武汉、黄石、鄂州的共10个环岛的实验中，距离的统计结果和与目视解译的对比结果。表1显示出，在地面分辨率4m*4m或3.7m*3.7m的情况下，两个中心点的距离均小于4米，即约一个像素的边长，粗略的拓扑判断能够发挥出比较好的作用，对不同的数据源都表现出了鲁棒性，提取效果可认为是比较精确的。表1还显示出，生产准确率(produceraccuracy，同时被自动化和目视解译标为环岛的像素数占目视解译标为环岛的像素数的百分比)在10次实验中都超过了85％，具有相当好的提取准确率，与已有的算法相比，具有明显的准确率优势。表2另外展示了10个环岛的位置、面积大小、圆度和ndvi，环岛的面积基本在预想范围内，符合交通部门的要求，圆度都基本在0.4-1.0间，符合几何规律和本发明的设想，ndvi基本都在该环岛所在的整幅或裁剪影像中拥有较高的数值，由此可以看出基于规则的筛选的思路是十分正确且有效的。这些结果可以证明，从高分辨率多波段遥感影像中自动化提取植被环岛的算法结果可用且准确，说明了从高分辨率多波段遥感影像中自动化提取植被环岛的算法的方案可行，具有高准确率、适应性强、自动化、不依赖其他辅助数据的优势。gtr(m)distance(m)produceraccuracy13.73.03691.67％23.73.27186.02％33.73.40985.71％43.72.30892.63％53.73.82187.37％63.73.28690.32％740.85895.47％842.46194.06％942.60496.67％1042.48694.74％表1十次实验的地面分辨率(gtr)、中心点距离(distance)和生产准确率(produceraccuracy)idcoordinatesrd.namesize(m2)roundnessndvi1114.368,30.571黄鹂路555.48030.5842010.0441862114.377,30.569沿湖路1353.0930.4739940.1224963114.255,30.626常青路897.31430.7107750.0723434115.059,30.255黄石路1623.7120.7374490.0941005114.378,30.580梨园路1285.9340.8167960.1173526115.097,30.211颐阳路1373.6330.7648840.0354407114.825,30.417四海路3926.9750.9081190.2295628114.920,30.398司徒路1600.0000.8109960.2397809114.882,30.391东吴路2432.0000.8555900.20769810114.882,30.367官柳路1104.0000.8206110.166888表2十次实验的环岛坐标(coordinates)、名字(rd.name)、大小(size)、圆度(roundness)和植被指数(ndvi)应当理解的是，本说明书未详细阐述的部分均属于现有技术。应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。当前第1页12