基于区域不一致性评价自动优选遥感影像分割参数的方法

摘要

本发明公开了一种基于区域不一致性评价自动优选遥感影像分割参数的方法。基于5组等间隔的尺度参数和固定的形状因子和紧凑度因子依次构建分割数据集，然后与相对应的参考数据集基于PSE‑NSR‑ED2不一致性分割参数评价体系获得ED

说明书

技术领域

本发明属于地球科学研究领域，特别涉及遥感地学空间统计分析和模式识别等方面，具体涉及基于区域不一致性评价自动优选遥感影像分割参数的方法。

背景技术

在遥感领域，针对不同应用和目的，不同的影像分割算法不断被提出，而多分辨率分割(Multiresolution Segmentation)算法的出现，被认为是遥感影像分割的一个里程碑。该算法在影像分割时，综合了影像的光谱信息和空间信息，能产生内部同质性最高的影像对象，其主要参数有尺度、形状因子、紧凑度因子，这些参数的不同组合会产生不同的分割结果。选择质量最好的分割结果的过程被称为参数优选，而参数优选必须要解决的问题是对分割质量的具体评价。因此，如何获得最优分割参数组合来评价影像分割质量是OBIA中必须解决的一个问题。

不一致性分割结果评价方法基于参考多边形(Reference Polygon)和对应的匹配分割多边形(Corresponding Polygon)之间的不一致性(Discrepancy)来度量当前参数组合所产生的分割数据的质量。它是一种客观的经验性评价方法(Empirical Method)，它是用几何不一致性度量的是参考多边形与匹配多边形之间面积的差异，而算术不一致性度量的是两者多边形数量的差异。

在不一致性评价体系(Potential Segmentation Error,PSE-Number-of-Segments Ratio,NSR-Euclidean Distance 2,ED2)中，PSE是潜在分割误差面积比，NSR是分割多边形数量比，ED2是PSE与NSR的欧几里得距离。斜U型(Euclidean Distance2,ED2-Scale Patterns,SP)模式是基于对PSE-SP,NSR-SP,ED2-SP曲线的分析提出的。在给定形状、紧凑度参数的情况下，影像分割单元平均面积随尺度参数的增序近似呈幂函数的形式单调递增，分割单元的面积随尺度参数的递增近似呈幂函数。相应地，分割单元的数量随尺度参数递减近似呈幂函数。ED2作为PSE和NSR的组合形式会随着尺度参数的变化呈现倾斜的U型曲线形式，如图1所示。

2、现有技术方案

分割质量评价的方法有不一致性法和优度法。基于不一致性评价方法通过比较参考数据集和分割数据集，从几何不一致性和算术不一致性两个方面对分割质量进行综合评价。主要的不一致性评价指标有Clinton等、Weidner等提出的一系列分割质量评价指标QR(Quality Rate)、UR(Under-Segmentation Rate)、OR(Over-Segmentation Rate)和ED1；Liu等基于几何不一致性和算术不一致性提出了ED2(Euclidean Distance 2)指标；而Yang等通过分析ED1、ED2系列评价指标，进一步提出了ED3、ED3-Modified和SEI(Segmentation Evaluation Index)。同时，Zhang等提出了F-measure指标和MOA(Multiscale Object Accuracy)和BCA(Bidirectional Consistency Accuracy)。基于优度法的分割参数质量评价指标主要是通过局部方差来进行构建的；等(2010)将Kim提出的局部方差评价最优分割尺度的方法进行了自动化，构建了最优尺度参数选取的工具ESP；等(2014)对ESP工具进行了改进。

通过分析和总结最优尺度分割参数选择方法，提出尺度评定应该结合形状、纹理等信息，最终实现尺度分割参数的自动选取，然而目前自动化方法多为非监督分割，在一定程度上存在选择的分割参数尺度偏大，对地物类别的分辨针对性不强，欠分割现象明显等问题，这会给后续影像分类带来不利的影响。

发明内容

本发明的目的是克服现有技术的不足，提供了一种基于区域不一致性评价自动优选遥感影像分割参数组合的方法。

本发明的技术方案如下：

本发明公开了一种基于区域不一致性评价自动优选遥感影像分割参数的方法，包括如下步骤：

步骤一：输入待分割的遥感影像和参考数据集，初始化参数，所述的初始化参数包括给定初始尺度分割参数区间值s₁和s₅，且s₅>s₁，给定尺度分割参数间最小步距d_min，给定ED2_min最大值L、形状因子、紧凑度因子和ED2_min极小值ζ；

步骤二：如果s₅-s₁>4d_min，则在初始尺度分割参数值s₁和s₅的基础上设置5个尺度分割参数及其步距，否则需重新设置初始尺度分割参数值s₁和s₅，这里d>d_min为约束条件，s₁，s₂，s₃，s₄，s₅五个尺度分割参数具体计算如下：

s₂＝s₁+d

s₄＝s₅-d

步骤三：根据s₁，s₂，s₃，s₄，s₅尺度分割参数上各自的分割数据集和参考数据集，分别统计其参考多边形数量、分割多边形数量，重叠面积、过分割面积、欠分割面积、依次计算出每一个尺度分割参数对应的不一致性度量参数对应的PSE_i，NSR_i和ED2_i，其中i＝{2,3,4,5}。具体计算流程如下：

在每个尺度分割参数中，用R＝{r_i:i＝1,2,...,m}表示m个参考多边形的集合，S＝{s_j:j＝1,2,...,n}表示分割多边形的集合，|r_i∩s_j|表示参考多边形r_i和分割多边形s_j相交部分的面积，|r_i|和|s_j|分别为参考多边形r_i和分割多边形s_j的面积，S'＝{s_k:k＝1,2,...,v}表示与参考数据集相对应的分割数据集的集合。定义Sa和Sb为集合S的两个子集，且满足匹配准则：

也就是说匹配准则是参考多边形和分割多边形相交部分的面积至少是参考多边形或匹配多边形的面积一半，则参考多边形相匹配的分割多边形的集合S'就为Sa和Sb的并集。同时定义∑|R_i|为m个参考多边形的总面积，∑|S_k|为跟m个参考多边形相匹配的分割多边形总面积，∑|R_i∩S_k|为参考数据集和匹配的分割数据集重叠面积，|r_i-s_j|＝|r_i|-|r_i∩s_j|是在匹配多边形之外的那部分参考多边形的面积为过分割面积，|s_j-r_i|＝|s_j|-|r_i∩s_j|是在参考多边形之外的那部分匹配多边形的面积为欠分割面积。因此，PSE、NSR和ED2可以表示为：

最后得到计算5点间的ED2最小值和最大值：

ED2_min＝min{ED2₁,ED2₂,ED2₃,ED2₄,ED2₅}

ED2_max＝max{ED2₁,ED2₂,ED2₃,ED2₄,ED2₅}

步骤四:模式匹配过程

分析ED2随s₁～s₅尺度分割参数变化的趋势并动态调整五个尺度分割参数，不断地匹配PSE-NSR-ED2不一致评价模型Case>

步骤五：依次按照形状因子＝0.1,0.2,…,0.9，紧凑度因子＝0.1,0.2,…,0.9的组合方式进行步骤一到步骤四的迭代运算，得到81个参数组合上的最优尺度分割参数；然后对这81组最优尺度分割参数进行再排序，从中选出的最小ED2组对应的尺度分割参数，形状因子及紧凑度因子三个参数组合即为最优分割参数组合。

优选的，所述的步骤四ED2随尺度分割参数变化的趋势具有Case a～Case q17种模式：

Case a.ED2_i不随尺度分割参数发生变化，即ED2₁＝ED2₂＝ED2₃＝ED2₄＝ED2₅，或者其间的最大、最小值之差小于预设的极小值ζ：

|ED2_max-ED2_min|＜ζ

若如此且ED2_max≥L，则意味着预设尺度分割参数明显偏大至ED2值变化不稳定区域，需要将尺度分割参数范围向小值方向，将s₅左移到s₁的位置，s₁左移4个步距，s₃为动态调整后s₁和s₅的值一半，s₂和s₄也相应左移4个步距：

s₅←s₁

s₁←s₁-4·d

s₂←s₁+d

s₄←s₅-d

其中“←”表示将箭头右端的值赋值给箭头左端的变量，即为尺度分割参数动态调整的过程，五个尺度分割参数的调整是同步进行的，没有先后之分；箭头右边的s1～s5为调整前的尺度分割参数值；以s₅←s₁为例；即表示将原s₁的值赋值给新的s₅，即表示将s₅左移到原s₁的位置。

参数的调整具体是根据尺度分割参数变化的趋势对s₁～s₅的动态调整，完成相应尺度分割参数“左移”、“右移”、“放大”和“缩小”等，本发明中所有提到的“←”均代表这个意思。

如果s₁<0，则使s₁＝5，s₅左移s₄的位置，d为动态调整后s₁和s₅的值四分之一，s₃为动态调整后s₁和s₅的值二分之一，s₂和s₄左移1个新获的步距：

s₁←5

s₅←s₄

s₂←s₁+d

s₄←s₅-d

回到步骤三继续；

否则当满足条件ED2_max<L时扩大尺度分割参数搜索范围，相应的s₃不变，s₁减小原来的2倍步距，s₅扩大原来的2倍步距，s₂和s₄分别替换原来的s₁和s₅)

s₂←s₁

s₄←s₅

s₁←s₂-2·d

s₅←s₄+2·d

如果扩大出现s₁<0，则使s₁＝5，s₅左移s₄的位置，d为动态调整后s₁和s₅的值四分之一，s₃为动态调整后s₁和s₅的值二分之一，s₂和s₄左移1个新获的步距：

s₁←5

s₅←s₄

s₂←s₁+d

s₄←s₅-d

回到步骤三继续；

如果扩大尺度分割参数范围搜索，仍出现ED2_i不随尺度参数发生变化，或其间的最大、最小值之差小于预设的极小值ζ其中之一的情况：

|ED2_max-ED2_min|＜ζ

且ED2_max≤L，则最优尺度分割参数可能位居其一，s_opt∈[s₁,s₅]；需作进一步的判断：如果前后两次运算均出现上述状况，对比前后两次运算ED_2min，取ED2_min较小的一次作为最优尺度分割参数的范围，s_opt∈[s₁,s₅]，则s_opt＝s₅，运行结束并报告结果。

Case b.ED2_i随尺度分割参数的增加而递减，即ED2₁≥ED2₂≥ED2₃≥ED2₄≥ED2₅，则最优尺度分割参数大于s₅，调整尺度分割参数设置，将s₁右移到s₃的位置，s₂右移到s₄的位置，s₃右移到s₅的位置，s₂和s₄也相应右移2个步距：

s₁←s₃

s₂←s₄

s₃←s₅

s₄←s₃+2·d

s₅←s₄+2·d

回到步骤三继续；

Case c.ED2_i在s₄处呈现最小值拐点，即ED2₁≥ED2₂≥ED2₃≥ED2₄并且ED2₄≤ED2₅，则最优尺度分割参数可能处于s₃和s₅之间，将尺度分割参数范围向小值方向，则使s₁右移到s₂的位置，s₂右移到s₃的位置，s₃右移到s₄的位置，s₄右移到s₅的位置，s₅右移1个步距：

s₁←s₂

s₂←s₃

s₃←s₄

s₄←s₅

s₅←s₄+d

回到步骤三继续；

Case d.ED2_i在s₃处呈现最小值拐点，即ED2₁≥ED2₂≥ED2₃，并且ED2₃≤ED2₄≤ED2₅，则最优尺度分割参数处于s₂和s₄之间；进而如果ED2₃≥L，则意味着预设尺度分割参数明显偏大至ED2值变化不稳定区域，需要将尺度分割范围向小值方向，将s₅左移到s₁的位置，s₁左移4个步距，s₃为动态调整后s₅和s₁的二分之一，s₂和s₄也相应左移4个步距：

s₅←s₁

s₁←s₁-4·d

s₂←s₁+d

s₄←s₂-d

如果s1<0，则使s₁＝5，s₅左移s₄的位置，d为动态调整后s₁和s₅的值四分之一，s₃为动态调整后s₁和s₅的值二分之一，s₂和s₄也相应左移1个新获的步距：

s₁←5

s₅←s₄

s₂←s₁+d

s₄←s₅-d

回到步骤三继续；

否则ED2₃<L，并且d＞d_min，则需要在s₂和s₄之间加密搜索，相应的s₃不变，将s₁左移到s₂的位置，s₅左移到s₄的位置，s₂和s₄也相应左移二分之一个步距：

s₁←s₂

s₅←s₄

直到d≤d_min，且满足ED2₃<L，则s_opt＝s₃。

Case e.ED2_i在s₂处呈现最小值拐点，即ED2₁≥ED2₂，并且ED2₂≤ED2₃≤ED2₄≤ED2₄，则最优尺度分割参数可能处于s₁和s₃之间，若ED2₂<L，且d<d_min，则s₂为最优参数；否则如果ED2₂≥L，则需要将尺度分割参数范围向小值方向，则使s₅左移到s₄的位置，s₄左移到s₃的位置，s₃左移到s₂的位置，s₂左移到s₁的位置，s₁左移1个步距：

s₅←s₄

s₄←s₃

s₃←s₂

s₂←s₁

s₁←s₁-d

如果s₁<0，则使s₁＝5，s₅左移s₄的位置，d为动态调整后s₁和s₅的值四分之一，s₃为动态调整后s₁和s₅的值二分之一，s₂和s₄也相应左移1个新获的步距：

s₁←5

s₅←s₄

s₂←s₁+d

s₄←s₅-d

回到步骤三继续；

否则当ED2₂<L，且d>d_min，在s₁和s₃处加密，则需要在s₂和s₄之间加密搜索，相应的s₃不变，将s₁左移到s₂的位置，s₅左移到s₄的位置，s₂和s₄也相应左移二分之一个步距：

s₁←s₂

s₅←s₄

回到步骤三继续；

Case f.ED2_i随尺度分割参数递增，即ED2₁≤ED2₂≤ED2₃≤ED2₄≤ED2₂，则最优尺度分割参数小于s₁或者邻近s₁，则使s₅左移到s₃的位置，s₃左移到s₁的位置，s₄左移到s₂的位置，s₂左移2个步距，s₁左移2个步距：

s₅←s₃

s₃←s₁

s₄←s₂

s₂←s₃-d

s₁←s₂-d

如果s₁<0，则使s₁＝5，s₅左移s₄的位置，d为动态调整后s₁和s₅的值四分之一，s₃为动态调整后s₁和s₅的值二分之一，s₂和s₄也相应左移1个新获的步距：

s₁←5

s₅←s₄

s₂←s₁+d

s₄←s₅-d

回到步骤三继续；

Case g-Case q.随尺度分割参数递增，ED2_i随尺度分割参数的变化趋势趋于不稳定，ED2i值随尺度分割参数变化而至不稳定区域，且ED2_min≥L，则意味着预设尺度分割参数明显偏大至ED2_i值变化不稳定区域，需要将分割尺度参数范围向小值方向左移，将s₅左移到s₁的位置，s₁左移4个步距，s₃为动态调整后s₁和s₅的值一半，s₂和s₄也相应左移4个步距：

s₅←s₁

s₁←s₁-4·d

s₂←s₁+d

s₄←s₂-d

如果s₁<0，则使s₁＝5，s₅左移s₄的位置，d为动态调整后s₁和s₅的值四分之一，s₃为动态调整后s₁和s₅的值二分之一，s₂和s₄也相应左移1个新获步距：

s₁←5

s₅←s₄

s₂←s₁+d

s₄←s₅-d

回到步骤三继续:否则遇到ED2_max≤L时，将放大后计算获得的这组数据与上一组数据的ED2_min进行比较，如果出现小于上一组最小值的情况就继续放大，直到找到最小的ED2_min，放大结束，确定最小ED2_min的区间，然后在最小区间进行加密，加密到d≤d_min，加密结束，取其ED2_min对应的尺度分割参数为最优尺度分割参数。

优选的，所述的步骤一中给定尺度分割参数间最小步距d_min＝1，给定ED2_min最大值L＝1.0、形状因子＝0.1、紧凑度因子＝0.1、ED2_min极小值ζ＝0.0001。

本发明将地学空间分析和模式匹配用到具体问题上，确定的最优分割参数可用地学空间分析以及基于几何不一致性和算术不一致性原理明确地加以解释，具有更高的理论可信度。

本发明方法，不拘泥于遥感影像分割参数优选算法，可以广泛应用于任何服从斜U型曲线分布的模型求最优值的问题，该方法具有广泛的实用性和应用价值。

本发明方法，实质上是遥感影像分析之前分割参数的预测，避免了选择参数的盲目性，解决了试错法带来的不确定性和穷举法带来的耗时，兼顾了工作效率的同时，也提高了面向对象遥感影像处理与分析的精度和自动化程度。

附图说明

图1为给定形状因子和紧凑度因子下斜U型ED2-SP模式示意图；

图2为ED2随尺度分割参数变化的17种基本模式示意图；

图3为本发明实施例的基于区域不一致性评价自动优选遥感影像分割参数的方法流程流程图；

图4为本发明实施例与穷举法选择遥感影像最优分割参数组合自动搜索的轨迹对比图；

图5为本发明实施例应用在广东省东莞市Quick bird高分辨率遥感影像中的分割效果图；

图6为本发明实施例应用在宁夏回族自治州中卫市World view2高分辨率遥感影像中的分割效果图。

具体实施方式

如图3所示，为本发明实施例的基于区域不一致性评价自动优选遥感影像分割参数的方法流程流程图；在本发明的所有实施例中，步骤一给定尺度分割参数间最小步距d_min＝1，给定ED2_min最大值L＝1.0、形状因子＝0.1、紧凑度因子＝0.1、ED2_min极小值ζ＝0.0001。

如图2所示为ED2随尺度分割参数变化的17种基本模式示意图。

1.本发明的方法与穷举法选择最优分割参数组合的对比

将本发明实施例应用到广东省东莞市和宁夏回族自治州中卫市三种传感器(Quick bird、Alos、World view2)多光谱和融合的6景高分辨率遥感影像中3种不同地物(耕地FL、居民住宅RB、坑塘WB)最优分割参数组合的选择，并与穷举法获得的结果进行比较，实验中用相同配置的计算机进行测试，将初始尺度参数范围进行统一，确保两种方法具有可比性，获得实验结果，其中穷举法实验结果如表1所示，本发明的方法结果如表2所示。

对比表1和表2，可以发现本发明实施例在相同的初始尺度参数范围下，基本都能产生较小数据量级的分割数据集，所耗费时间也相对也较短，且可以将步长精确到1，获得更加精确的ED2的取值来确定遥感影像的最优分割参数组合，然而目前流行的穷举法需要机械地产生步长规定的分割数据集，步长一旦确定也不能动态调整，找到的初始尺度参数范围，也会因为初始尺度参数范围考虑不全面，没有列举到最优分割参数组合存在的范围，给寻找参数带来了具体的制约性。通过这个实验的对比，本发明实施例表现良好，不但可以通过设置合适的初始尺度参数范围快速的得到最优尺度分割参数，而且可以动态改变步长，可以精准预测PSE‐NSR‐ED2不一致性评价模型取得最优尺度分割参数的区间，节约时间的同时避免了选择最优分割参数组合的盲目性。

表1穷举法结果统计表

表2本发明方法结果统计表

从表1和表2获得大量的数据中抽取一组形状因子和紧凑度因子都为0.1时，选择最优尺度分割参数的轨迹进行分析，可以清楚的看到本发明实施例中可以在初始尺度分割参数的范围内快速预测到最优尺度分割参数的存在的区间，然后开始加密确定最优，而穷举法没有任何搜索轨迹，需要人工判断斜U型曲线最优的位置，如图4所示。

2.本发明的方法与其他自动化分割参数优选工具对比

将本发明实施例与等(2014)提出的自动优选遥感影像尺度参数的工具ESP进行对比，实验数据为广东省东莞市和宁夏回族自治州中卫市三种传感器(Quick bird、Alos、World view2)3景高分辨率遥感影像中3种不同地物(耕地、林地、坑塘)，由于ESP只能选择最优尺度分割参数，将本发明实施例中找到的形状因子和紧凑度因子作为标准，并设置相同的初始尺度参数和步长进行对比实验。根据分割精度评价指标中的规定，生产者精度和用户精度都接近1，分割效果最优的原则衡量表3中的结果，可以发现ESP工具的用户精度很不稳定，最低达到0.11，而且找到的最优尺度分割参数明显偏大，分割数据集量级很大，整体效率较差。

表3最优分割参数结果及精度评价对比统计表

3.图5为本发明实施例和穷举法获得的最优分割参数组合对广东省东莞市Quick bird高分辨率遥感影像对水体的分割效果，可以看出已经能很好的划分水体的区域，过分割和欠分割现象都解决良好；图6为本发明实施例和穷举法获得的最优分割参数组合对宁夏回族自治州中卫市World view2高分辨率遥感影像对耕地的分割效果,，明显可以看出本发明实施例已经可以预测出此类地物的最优尺度分割参数范围比较大，但是穷举法已经盲目认为初始尺度参数范围已经估计到了最大范围，这样人工的判读存在严重的主观性。