本发明涉及地表覆盖变化识别方法，具体涉及一种单期影像地表覆盖变化检测自动化方法。
背景技术：
遥感影像变化检测是全球变化研究的重要内容，已应用于诸多领域，如灾后应急与评估、环境变化检测和空间数据更新。
对变化检测方法总体可以分为两大类：一类是基于两期遥感影像的变化检测方法，另一类是使用矢量数据和遥感影像相结合的变化检测方法。由于影响地物光谱特性变化的因素比较复杂，因此，基于两期遥感影像的变化检测方法存在工作量大、分类误差累积现象严重且对数据条件要求苛刻等诸多不足。
现有技术中多通过影像对象间纹理特征的差异性进行变化检测，有效地提高了变化检测结果的精度和效率。但是，在变化检测中，需要以人工目视判别的方式从分割后的影像对象中抽取一定数量、具有相同类别属性的样本，工作量依然巨大，费时费力，且样本的抽取结果存在一定的主观性，会严重影像遥感影像变化检测的精度。
有鉴于此需要提供一种单期影像地表覆盖变化检测自动化方法。
技术实现要素：
本发明所要提供的是一种单期影像地表覆盖变化检测自动化方法，其能够利用一期已有的矢量数据结合单期新遥感影像实现地表覆盖变化目标检测的自动化，且检测精度高。
为实现以上发明目的，本发明提供一种单期影像地表覆盖变化检测自动化方法，包括如下步骤：a)利用已有矢量数据分割遥感影像，获取具有先验地表覆盖类型的初始影像对象；b)根据所述初始影像对象中包含的地表覆盖类型以及所述初始影像对象的空间分布抽取该研究区域内各类所述地表覆盖类型的影像样本；c)计算各类地表覆盖类型的所述影像样本的纹理特征值，以基于纹理特征贡献度得到各类地表覆盖类型的优选纹理特征，并剔除异常样本，保留正确样本，形成待检测影像对象；d)计算各所述待检测影像对象的优选纹理特征值，将所述待检测影像对象的优选纹理特征值与同类别影像样本中所述正确样本的纹理特征值对比，以获得变化影像对象；e)对所述变化影像对象重新分割并分类，以得到变化检测结果。
具体地，所述步骤b)中，所述初始影像对象的空间布设根据研究区域范围、研究区域内所述影像对象的分布特征和研究区域内的地形差异得出，包括如下步骤：
1)将所述研究区域划分为r×c个抽样格网；
2)根据所述地形差异将该研究区域划分为t个层级；
3)根据所述抽样格网内所述初始影像对象的个数和所述地形差异计算所述抽样格网内布设的影像样本个数。
进一步具体地，所述步骤c)中，各类所述地表覆盖类型的优选纹理特征由该类地表覆盖类型的所述影像样本的纹理特征贡献度得出；各类所述地表覆盖类型的影像样本的所述异常样本由该类地表覆盖类型的影像样本的优选纹理特征与纹理异常度指数得出。
进一步具体地，各类所述影像样本的纹理特征贡献度由所述纹理特征的信息增益率得出。
进一步具体地，各类所述影像样本的所述纹理异常度指数由该影像样本中样本对象的局部可达密度得出。
进一步具体地，所述局部可达密度由所述样本对象的第k邻域和该样本对象与另一样本对象间的第k可达距离得出。
进一步具体地，所述第k可达距离为所述样本对象的第k距离和该样本对象与所述另一样本对象在优选特征空间向量内的欧式距离中的最大值。
进一步具体地，所述优选特征空间向量由同一类地表覆盖类型的所述优选纹理特征根据所述纹理特征贡献度加权构建。
进一步具体地，所述步骤d)中，所述变化影像对象由所述待检测影像对象的所述优选纹理特征值与同类别地表覆盖类型的所述影像样本中所述正确样本的纹理特征值之间的纹理异常度指数，根据设定的异常度阈值筛选得出。
更具体地，所述步骤e)中，所述变化影像对象采用多尺度分割算法进行分割，计算分割后的所述变化影像对象在各类所述地表覆盖类型的所述影像样本中的纹理异常度指数，并据此划分该变化影像对象的分类类别。
本发明的单期影像地表覆盖变化检测自动化方法，先将具有先验地表覆盖类型的初始影像对象按照地表覆盖类型以及各地表覆盖类型在空间上的分布进行抽样，得到各类地表覆盖类型的影像样本，并将各地表覆盖类型的影像样本根据纹理特征贡献度加权构建成优选特征空间向量，同时根据影像样本的纹理异常度指数判断该影像样本对象是否为异常样本，若是则剔除该异常样本，并将保留的正确样本形成为待检测影像对象的样本对象，异常样本即是指与先验类别属性不一致的样本，检测并剔除这类样本能够确保自动抽取样本的正确性，从而能够提高变化检测结果的精度。随后计算待检测影像对象的优选纹理特征值与同类别影像样本中正确样本的纹理特征值之间的纹理异常度指数，并根据设定的异常度阈值筛选得出变化影像对象，最后采用多尺度分割算法对变化影像对象进行分割，计算分割后的变化影像对象在各类所述影像样本中的纹理异常度指数，并据此划分该变化影像对象的地物类别，得到各类地表覆盖的变化情况。本发明能够实现对已有一期地表覆盖矢量数据基础和单期新影像条件的地表覆盖变化的自动化检测，既降低了两期(多期)遥感影像变化检测对数据的苛刻要求；又实现了样本的自动抽取，避免了由人工目视判别造成的主观性，且由于剔除了自动抽取样本中的异常样本(包括错分或已变化样本)，并选用贡献度高的特征进行异常检测，能够提高变化检测结果的精度。
本发明实施例的其它特征和优点将在随后的具体实施方式部分予以详细说明。
附图说明
图1是本发明单期影像地表覆盖变化检测自动化方法的原理图；
图2是本发明单期影像地表覆盖变化检测自动化方法的样本自动提取流程图；
图3是样本对象在三维特征空间的可达距离计算原理图；
图4是本发明单期影像地表覆盖变化检测自动化方法的地表覆盖变化检测流程图；
图5是本发明单期影像地表覆盖变化检测自动化方法的变化影像对象分类流程图；
图6是本发明单期影像地表覆盖变化检测自动化方法中的耕地空间分布特征图，其中，图6-1是耕地类影像对象图；图6-2是海拔650～660m的耕地类影像对象图；图6-3是海拔660～670m的耕地类影像；图6-4是海拔670～680m的耕地类影像；图6-5是海拔680～690m的耕地类影像；图6-6是海拔690～700m的耕地类影像；图6-7是海拔700～710m的耕地类影像；图6-8是海拔710～720m的耕地类影像；
图7是本发明单期影像地表覆盖变化检测自动化方法中的林地空间分布特征图，其中，图7-1是林地类影像对象图；图7-2是海拔650～660m的林地类影像对象图；图7-3是海拔660～670m的林地类影像；图7-4是海拔670～680m的林地类影像；图7-5是海拔680～690m的林地类影像；图7-6是海拔690～700m的林地类影像；图7-7是海拔700～710m的林地类影像；图7-8是海拔710～720m的林地类影像；
图8是本发明单期影像地表覆盖变化检测自动化方法中的居民地空间分布特征图，其中，图8-1是居民地类影像对象图；图8-2是海拔650～660m的居民地类影像对象图；图8-3是海拔660～670m的居民地类影像；图8-4是海拔670～680m的居民地类影像；图8-5是海拔680～690m的居民地类影像；图8-6是海拔690～700m的居民地类影像；图8-7是海拔700～710m的居民地类影像；图8-8是海拔710～720m的居民地类影像；
图9是本发明单期影像地表覆盖变化检测自动化方法中影像样本的布设结果图，其中，图9-1是耕地类影像样本布设结果图；图9-2是林地类影像样本布设结果图；图9-3是居民地类影像样本布设结果图；
图10是本发明单期影像地表覆盖变化检测自动化方法中耕地类异常度指数检测频率分布图；
图11是本发明单期影像地表覆盖变化检测自动化方法中林地类异常度指数检测频率分布图；
图12是本发明单期影像地表覆盖变化检测自动化方法中居民地类异常度指数检测频率分布图；
图13是本发明单期影像地表覆盖变化检测自动化方法中各类影像样本的提取结果与异常样本图，其中，图13-1是耕地类影像样本提取结果；图13-2是林地类影像样本提取结果；图13-3是居民地类影像样本提取结果；图13-4是耕地异常样本；图13-5是林地异常样本；图13-6是居民地异常样本。
具体实施方式
以下结合附图对本发明实施例的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明实施例，并不用于限制本发明实施例。
如图1和图4所示，在本发明所提供的单期影像地表覆盖变化检测自动化方法的一个实例中，该方法包括如下步骤：
a)利用已有矢量数据分割遥感影像，获取具有先验地表覆盖类型的初始影像对象：
具体地，是利用已有的一期历史矢量数据对最新的遥感图像数据进行分割从而得到具有先验地表覆盖类型的初始影像对象，该方法能够充分利用历史矢量数据中图斑的大小、形状、类型信息，避免了现有基于影像特征的遥感影像分割效果不理想、分类精度不高等问题，提高了新影像分割与分类的精度。
b)根据初始影像对象中包含的地表覆盖类型以及初始影像对象的空间分布抽取研究区域内各类地表覆盖类型的影像样本；
具体地，初始影像对象的空间布设是根据该究区域的范围、研究区域内初始影像对象的分布特征和研究区域内的地形差异得出的，其具体的得出步骤为：
首先，将研究区域划分为r×c个抽样格网；其次将抽样区域按地形差异划分为t个层级，然后根据抽样格网内初始影像对象个数和地形起伏情况计算抽样格网内布设的影像样本个数，则第r行第c列抽样格网内影像样本布设个数snr×c的计算方法为：
式中，ion为抽样区域内初始影像对象的总数，stn为影像样本布设总数，
c)计算各类地表覆盖类型的影像样本的纹理特征值，基于纹理特征贡献度得到各类地表覆盖类型的优选纹理特征，并剔除异常样本，保留正确样本，形成待检测影像对象：
首先，优选纹理特征的获取过程如下：计算影像样本的纹理特征的信息增益，设gain(ci,fj)表示第j个纹理特征参数fj的信息增益，则有：
gain(ci,fj)＝h(ci)-h(ci/fi)
即gain(ci,fj)表示了第j个纹理特征fj的特征值已知时，地物类别ci信息量增减的程度，纹理特征的信息增益越大，表示纹理特征fj对地物类别ci分类结果的影响越大，因此，在进行纹理特征选择时，通常选择信息增益较大的纹理特征构建优选特征空间向量。但是在使用信息增益选择纹理特征时，选择结果往往会偏向于具有更多取值区间的纹理特征，导致纹理特征评价结果不准确，因此，在获取各类地表覆盖类型的优选纹理特征时可以进一步优选采用信息增益率作为选择纹理特征的依据，以效消除上述的不良影响。具体方法是，设gainrat(ci,fj)表示第j个纹理特征参数fj的信息增益率，则有：
gainrat(ci,fj)＝gain(ci,fj)/h(fj)
其中，
其中，maxf(·)表示纹理特征fj在不同地物类别中的信息增益率的最大值，maxc{·}表示识别地物类别为ci的不同纹理特征的信息增益率的最大值。采用纹理特征贡献度能够量化同一纹理特征对不同类别的地物以及不同纹理特征对同类地物的相对贡献大小。
随后，剔除某类地表覆盖类型的影像样本中与该影像样本所属的地表覆盖类型不一致的异常样本，以保留的正常样本形成待检测影像对象的样本对象，以能够对在利用已有的历史矢量数据对最新的遥感图像数据进行分割时对影像对象地表覆盖类型的判别是否产生错误进行检验，从而能够进一步地提高新影像分割与分类的精度，如图2所示，筛选异常样本的具体步骤为：
首先，建立优选特征空间向量，该优选特征空间向量由同一类地表覆盖类型的优选纹理特征根据纹理特征贡献度加权构建，例如，图3中样本对象obj的特征空间fsp(obj)可以表示为：
其中
其次，计算样本对象的局部可达密度，具体方法为：首先，计算该样本对象obj到另一样本对象obi之间第k可达距离(rdisk(obj,obi))：
rdisk(obj,obi)＝max(k-dis(obi),d(obj,obi))
其中，d(obj,obi)表示优选特征空间向量内样本对象obi到样本对象obj的欧氏距离，例如图3中的d′2即表示样本对象obj至样本对象ob2之间的欧氏距离；k-dis(obi)表示样本对象obi的第k距离，即样本对象obi与含有k个样本对象的邻域中与样本对象obi相距最远的样本对象之间的距离，例如图3中的d1即表示样本对象ob1与含有k个样本对象的邻域中和样本对象ob1相距最远的样本对象之间的距离。随后，根据样本对象obi到样本对象obj之间第k可达距离计算样本对象obj的局部可达密度(lrd(obj))：
其中，k为邻域参数，表示一个样本对象邻域内应包含最少的样本对象个数，nk(obj)表示对象obj的第k邻域。
最后根据样本对象obj的局部可达密度计算样本对象obj的纹理异常度指数(fsoi(obj))：
其中，d表示样本布设对象的集合，maxi∈d{lrd(obj)}表示样本对象集合中lrd(obj)的最大值。
将计算出的样本对象obj的纹理异常度指数fsoi(obj)与设定的阈值对比，若大于设定的阈值则该样本对象为异常样本，则将该样本对象剔除，以使得影像样本成为只保留有正确的样本对象。在异常样本筛选中，“样本对象obi的第k距离”中k的取值会影响样本对象的纹理异常度指数，一般按样本总数的1/5～1/3取值时，纹理异常度指数相对趋于稳定。纹理异常度指数阈值的设定不同也会使得异常样本筛选结果出现差异，较低的阈值可以获得异常样本较低的漏检率和较高的误检率，而较高的阈值会获得异常样本较高的漏检率和较低的误检率。在样本提取结果中，应当更侧重漏检率，可以以牺牲一定的误检率来换取0漏检率，一般情况下，将异常度阈值设定在80％或70％时可以获得0漏检率和较高的检测精度。
d)计算各待检测影像对象的优选纹理特征值，将该待检测影像对象的优选纹理特征值与同地表覆盖类别的影像样本中正确样本的纹理特征值进行对比，以获得变化影像对象；即是通过计算待检测影像的纹理异常度指数，并设定异常度阈值以区分变化和未变化影像对象。
e)对变化影像对象重新分割并分类，以得到变化检测结果，具体是根据变化影像本身的特征采用多尺度分割方法对变化影像对象进行重新分割，然后对重新分割后的变化影像对象进行分类。由于重新分割后变化影像对象的先验类别未知，因此，需要计算变化影像对象在所有类别的影像样本中的纹理异常度值，异常度值最小且小于一定阈值(例如20％)的类别即为该变化影像对象的分类类别，如图5所示，以影像样本包括耕地类影像样本、林地类影像样本、水体类影像样本和居民地类影像样本为例，则待分类的变化影像对象需要计算在上述四个类别的影像样本中的纹理异常度值，以判断该变化影像对象属于哪一类地表覆盖类型。
上述的技术方案首先能够有效实现研究区域遥感影像变化检测所需样本的自动抽样，避免了繁琐的人工抽样，减少了抽样的工作量，提高了遥感影像变化检测的效率；其次，该方法是通过矢量数据提取抽样图层的先验信息进行分层样本空间布设，能够针对特定的变化检测目标自动抽取样本，如在灾后应急救援中，一般更关注的是房屋和道路的破坏程度，则该方法能够提取特定的房屋和道路样本；最后该方法实现了遥感影像变化检测的自动化，特别是单期影像变化检测与分类的自动化，提高了变化检测精度，降低了遥感影像变化检测对遥感影像的获取传感器、时间、分辨率等苛刻条件。
下面将以一具体实施例对本发明的地表覆盖变化检测过程作进一步说明：
以耕地、林地和居民地三种地表覆盖类型为例，如图2所示，首先从基准期矢量地图中分别提取耕地、林地和居民地抽样图层，并将抽样区域依据制图比例尺划分为10cm*10cm的规则抽样格网，如制图比例尺为1:1000，则10cm*10cm的规则格网的实地水平距离100m×100m。然后用抽样图层分割影像数据，获取影像对象。根据dem数据(高程数据)表示的地形特征，根据研究区域高差，按照等高距原则将抽样区域划分为10个以内高程抽样层级，耕地、林地和居民地影像对象在不同抽样层级中的分布特征如图6至图8所示。
实验数据抽样区域内包含耕地影像对象总数1028个、林地影像对象总数625个、居民地影像对象总数159个，本试验中，在抽样区域内计划布设耕地类影像样本224个、林地类影像样本167个、居民地类影像样本92个。根据耕地、林地和居民地影像对象在各层级的分布特征，首先计算出各抽样格网内每一层级布设的样本个数，然后采用随机抽样方法，对每个格网进行样本布设，抽样区域样本布设结果如图9所示。
在影像对象异常样本筛选的过程中，首先构建影像对象的优选特征空间向量，根据初始影像对象先验类别属性，耕地由角二阶矩、对比度、逆差矩、熵、平方和、差熵和差方差7个特征参数构成优选特征空间向量，林地由角二阶矩、逆差矩、熵、均值、总方差和总平均6个特征参数构成优选特征空间向量，居民地由逆差矩、熵、对比度、差熵和角二阶矩5个特征参数构成优选特征空间向量。其次是计算样本对象的第k距离和第k可达距离，然后计算局部可达密度，并由局部可达密度计算异常度指数。k值的大小，会影响第k距离的大小，进而会影响影像对象的异常度指数。图10至图12表示了不同地表覆盖类别在不同k值下样本异常度大小的频率分布直方图。
从图10至12中可以直观看出，随着k值的增大，低异常度样本的个数逐渐增多，高异常度样本的个数逐渐减少。而且当k值达到一定取值区间时(如图10中k取值在35～90)，高异常度样本个数趋于稳定，而当k取值过大时(接近样本总数)，样本中几乎全是低异常度样本(如图12中k取120)，则会导致异常检测结果失败。影像对象异常检测就是要去除样本中高异常度的影像对象，因此，由图10至图12中的试验数据可知，k值按样本总数的1/5左右选择时，不会影响异常检测结果的准确性。
为了定量分析影像对象异常检测的总体效果，在不同k值和不同异常度阈值下，对异常检测中的误检、漏检和总体精度进行了分析，其中误检指原本为与先验类别一致的样本被误判为异常样本，漏检指原本为异常样本而没有被检测出来。耕地、林地和居民地异常检测结果精度分析见表1。
表1耕地、林地和居民地异常检测精度分析
从表1可知，在异常度阈值取80％，k值低于90时，耕地的漏检率为0，林地的漏检率为21.3，居民地的漏检率为32.8。在异常度阈值取80％，k值低于65时，林地的漏检率为0,而居民地只有当k值低于50时，漏检率才为0，其原因是林地的样本总数比耕地少，而居民地的样本总数最少，所以k值相对于样本总数的比例对异常检测结果精度有直接影响。而且，当耕地k值在区间50～90取值，林地k值在区间20～65取值，居民地k值在区间10～50取值，异常度阈值取80％时，误检率和漏检率均为零，总体检测精度为100％。在设定异常度阈值为70％的条件下，无论k取何值时，耕地、林地和居民地的漏检率均为0。而在样本影像对象异常检测中，只要能满足所抽取样本的验后类别与先验类别属性一致，则为成功抽样，因此，在异常检测中可以以牺牲一定的误检率来换取0漏检率。所以，为了获得0漏检率和较高的检测精度，一般情况下，k可以按布设样本总数的1/5～1/3取值，异常度阈值可设定为80％。本试验中，耕地取k＝50、林地取k＝35、居民地取k＝20，异常度阈值设定为80％，异常样本筛选结果，以及剔除了异常样本的地表覆盖提取样本的结果如图13所示。
以上结合附图详细描述了本发明实施例的可选实施方式，但是，本发明实施例并不限于上述实施方式中的具体细节，在本发明实施例的技术构思范围内，可以对本发明实施例的技术方案进行多种简单变型，这些简单变型均属于本发明实施例的保护范围。
另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明实施例对各种可能的组合方式不再另行说明。
此外，本发明实施例的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明实施例的思想，其同样应当视为本发明实施例所公开的内容。
