一种结合多种遥感指标的裸土变化识别方法与流程

本发明涉及时序遥感变化监测技术领域，特别是一种结合多种遥感指标的裸土变化识别方法。

背景技术：

水土流失是当今世界面临的一个严峻的生态环境问题，由其诱发的各种灾害也给人们的生产生活带来了威胁和损害，而其中裸土是水土流失区的主要地表景观。对此，快速并准确地定位裸土区域，进行地表裸土的时空变化分析，对于水土流失的治理工作尤其重要。近年来，遥感技术已成为水土流失治理和评价不可或缺的高新技术。对此，众多学者采用遥感技术对裸土变化进行研究，其中主要方法包括三类：(1)通过人工目视解译方法进行对研究区遥感影像进行裸土识别；(2)通过采取数据样本，采用监督分类方法对研究区进行裸土识别，结合多个时相分类结果对比分析裸土时空变化特征。(3)为了解决建筑用地和裸土混杂的情况，有学者提出一种双重指数法，该方法首先利用裸土指数并阈值法提取裸土，接受利用归一化建筑指数来去除裸土中混杂的建设用地。

当前众多学者对裸土变化识别方法进行研究，但依然存在以下不足：1，基于人工目视解译的方法存在工作量大，获取大范围裸土信息比较困难等问题。2，基于监督分类的方法，在训练样区选取上难度大，而非监督分类方法获取的分类结果精度低。3，由于存在类内异质性，该方法不适于大范围裸土识别；4，进行裸土时空演变分析时，采用分类后比较方法会带来较大误差。

技术实现要素：

有鉴于此，本发明的目的是提供一种结合多种遥感指标的裸土变化识别方法，结合植被丰度和裸土指数变化趋势，最终确定裸土增加或减少的区域；简单易操作且结果准确度高，不需要设置阈值，能很好地应用于大范围长时间序列裸土变化区域快速监测。

本发明采用以下方案实现：一种结合多种遥感指标的裸土变化识别方法，包括以下步骤：

步骤s1：建立植被近红外反射指数时序数据集；

步骤s2：建立裸土指数时序数据集；

步骤s3：建立特征指标时序数据集；

步骤s4：特征指标评估和选取；

步骤s5：计算特征指标的变化趋势；

步骤s6：建立裸土变化识别方案；

步骤s7：进行裸土变化识别，获得研究区裸土变化空间分布图：根据步骤s6所建立的裸土变化识别方案逐像元判断裸土是否发生变化，最终生成研究区裸土变化识别结果图。

进一步地，步骤s1的具体内容为：在研究区范围内，对采集的图像逐像元逐年建立研究时段内植被近红外反射指数的8天最大化合成多年时序数据集；利用whittakersmoother平滑方法，对原始8天最大化合成植被近红外反射指数多年时序数据平滑处理，从而逐像元获得多年平滑的逐日植被近红外反射指数时序数据集。

进一步地，所述步骤s2的具体内容为：

裸土指数(简称ribs)，其计算公式为：

其中，nndsi指归一化差分积雪指数，nbrightness是通过缨帽变换获取的亮度指标归一化得到，其计算公式分别为：

逐像元逐年逐期计算裸土指数，生成研究时段内研究区裸土指数时序数据集。

进一步地，所述步骤s3的具体内容为：对研究区，逐年逐像元提取包括植被近红外反射指数和裸土指数的年最大值、年均值和50th分位值即植被丰度特征指标，并建立研究时段内研究区各特征指标时序数据集。

进一步地，所述步骤s4的具体内容为：选取若干不同地物类型的稳定点位，具体地物类型包括：植被、不透水面和裸土，分别评估基于植被近红外反射率指数和裸土指数的年最大值、年均值和50th分位值即植被丰度特征指标在不同地物类型中的分离度，分离度计算公式：

其中，μi、μs分别代表地物类型a和地物类型b的均值，σi、σs分别代表地物类型a和地物类型b的标准差，sdi用于比较两种不同地物类型的分离度；sdi越大，分离度越大；如果值小于1，表明分离度差；最终选择拥有高分离度的特征指标，经过计算选取植被丰度和裸土指数年均值作为裸土变化识别的指标。

进一步地，所述步骤s5的具体内容为：利用一元线性拟合的方法，分别基于研究时段内特征指标时序数据集，逐像元依次计算研究时段内特征指标的变化趋势k；计算公式如下：

y＝kx+b

当k>0时，表示该时序曲线呈上升趋势；当k<0时，表示该时序曲线具有下降趋势；基于mann-kendall方法，进一步判断该时序曲线的变化趋势是否显著；依据显著性检验的结果，将变化趋势分为三种情况：显著的上升趋势、没有趋势和显著的下降趋势；

基于mann-kendall方法，进一步判断该时序曲线的变化趋势是否显著；假设有n个样本量x1……xn的时间序列，对于所有k,j≤n,且k≠j，xk和xj的分布是不同的，计算检验统计量s，公式如下:

其中，s为正态分布，均值为0，方差var(s)＝n(n-1)(2n+5)/18；当n>10时，标准的正态统计变量通过下式计算:

对于统计值z来说，大于0时，表示增加趋势；小于0时，表示减少趋势；z的绝对值在大于1.96时，表示通过了置信度95％的显著性检验；当z>1.96时，表示该时序曲线呈显著上升趋势；当z<-1.96时，表示该时序曲线具有显著下降趋势；当z在正负1.96之间时为无趋势。

进一步地，步骤s6中所述的建立裸土变化识别方案的具体内容为：根据特征指标的变化趋势进行组合，将裸土指数和植被丰度均呈显著下降趋势的像元，判断为裸土增加；将裸土指数和植被丰度均呈显著上升趋势的像元，判断为裸土减少。

与现有技术相比，本发明具有以下有益效果：

(1)本发明适用于大范围裸土变化监测，时空连续性强，不需要事先提取裸土，且能够快速获取长时间序列裸土变化空间分布图。

(2)本发明充分考虑到不同地物类型在不同遥感指标中呈现的差异性，即植被丰度和裸土指数均最低，而其他地物通常在植被丰度较低的情况下，具有较高的裸土指数，进而依据植被丰度、裸土指数判断是否发生裸土变化。

(3)本发明充分利用多指标多年变化趋势，判断是否发生裸土增加或减少情况，相较于现有技术，能够提高识别结果的精度。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的特征指标评估图；其中图2(a)为不同地物类型点位的植被指数年最大值箱线图；图2(b)为不同地物类型点位的植被指数年均值箱线图，图2(c)为不同地物类型点位的植被丰度(植被指数计算的)箱线图，图2(d)为不同地物类型点位的裸土指数年最大值箱线图，图2(e)为不同地物类型点位的裸土指数年均值箱线图，图2(f)为不同地物类型点位的植被丰度(裸土指数计算的)箱线图。

图3为本发明实施例的裸土减少点位各遥感指标时序信号图,其中，图3(a)为裸土减少点位裸土指数原始信号图，图3(b)为裸土指数年均值年际趋势图，图3(c)为植被指数原始信号图，图3(d)为植被丰度指标年际趋势图。

图4为本发明实施例的裸土增加点位各遥感指标时序信号图,其中，图4(a)为裸土增加点位裸土指数原始信号图，图4(b)为裸土指数年均值年际趋势图，图4(c)为植被指数原始信号图，图4(d)为植被丰度指标年际趋势图。

图5为本发明实施例的特征指标趋势组合图。

图6为本发明实施例的裸土变化识别流程图。

图7为本发明实施例的研究区裸土变化空间分布图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

如图1所示，本实施例提供一种结合多种遥感指标的裸土变化识别方法，

首先基于植被近红外反射指数和裸土指数构建了多种特征指标，并进一步评估和选取最能体现不同地物类型差异的特征指标，根据不同地物类型之间的分离度结果表明，基于植被近红外反射指数构建的植被丰度和裸土指数的年均值最能体现不同地物类型的差异。通过判断植被丰度、裸土指数是否存在显著的变化趋势，最终依据植被丰度、裸土指数的变化趋势，建立一种结合多种遥感指标的裸土变化识别方法。

包括以下步骤：

步骤s1：建立植被近红外反射指数时序数据集；

步骤s2：建立裸土指数时序数据集；

步骤s3：建立特征指标时序数据集；

步骤s4：特征指标评估和选取；

步骤s5：计算特征指标的变化趋势；

步骤s6：建立裸土变化识别方案；

步骤s7：对研究区裸土变化进行识别，获得研究区裸土变化空间

分布图：根据步骤s6建立的裸土变化识别方案逐像元判断裸土是否发生变化，最终生成研究区裸土变化识别结果图。

在本实施例中，步骤s1的具体内容为：在研究区范围内，对采集的图像逐像元逐年建立研究时段内植被近红外反射指数的8天最大化合成多年时序数据集；利用whittakersmoother平滑方法，对原始8天最大化合成植被近红外反射指数多年时序数据平滑处理，从而逐像元获得多年平滑的逐日植被近红外反射指数时序数据集。

在本实施例中，在遥感领域，大多时序变化监测方法一般时间跨度都会在10年以上，但这个也没有具体的规定，时间跨度越长，方法监测结果精度会更高。在本实施例中，研究时段在10年以上。

在本实施例中，所述步骤s2的具体内容为：

裸土指数(简称ribs)，其计算公式为：

其中，nndsi指归一化差分积雪指数，nbrightness是通过缨帽变换获取的亮度指标归一化得到，其计算公式分别为：

逐像元逐年逐期计算裸土指数，生成研究时段内研究区裸土指数时序数据集。

在本实施例中，所述步骤s3的具体内容为：对研究区，逐年逐像元提取包括植被近红外反射指数和裸土指数的年最大值、年均值和50th分位值即植被丰度特征指标，并建立研究时段内研究区各特征指标时序数据集。

在本实施例中，所述步骤s4的具体内容为：选取若干不同地物类型的稳定点位，具体地物类型包括：植被、不透水面和裸土。探索植被、不透水面和裸土在裸土指数和植被指数不同分位值上表现出的差异性；分别评估基于植被近红外反射率指数和裸土指数的年最大值、年均值和50th分位值即植被丰度特征指标在不同地物类型中的分离度如图2，并选择拥有最佳分离的特征指标；根据评估结果表明：基于植被近红外反射指数构建植被丰度指标和裸土指数年均值能更好区分三种地物类型，因此，最终选取植被丰度和裸土指数年均值作为裸土变化识别的指标。

分离度计算公式：

其中μi、μs分别代表地物类型a和地物类型b的均值，σi、σs分别代表地物类型a和地物类型b的标准差，sdi用于比较两种不同地物类型的分离度；sdi越大，分离度越大；如果值小于1，表明分离度差；

根据对比结果表明，裸土指数年均值计算的不同地物类型两两之间的sdi数值要高于年最大值和50th分位值；而植被指数50th分位值计算的不同地物类型两两之间的sdi数值要高于年最大值和年均值。因此。最终选取植被丰度和裸土指数年均值作为裸土变化识别的指标。

在本实施例中，所述步骤s5的具体内容为：利用一元线性拟合的方法，分别基于研究时段内特征指标时序数据集，逐像元依次计算研究时段内特征指标的变化趋势k；计算公式如下：

y＝kx+b

当k>0时，表示该时序曲线呈上升趋势；当k<0时，表示该时序曲线具有下降趋势；基于mann-kendall方法，进一步判断该时序曲线的变化趋势是否显著；依据显著性检验的结果，将变化趋势分为三种情况：显著的上升趋势、没有趋势和显著的下降趋势；

基于mann-kendall方法，进一步判断该时序曲线的变化趋势是否显著；假设有n个样本量x1……xn的时间序列，对于所有k,j≤n,且k≠j，xk和xj的分布是不同的，计算检验统计量s，公式如下:

其中，s为正态分布，均值为0，方差var(s)＝n(n-1)(2n+5)/18；当n>10时，标准的正态统计变量通过下式计算:

对于统计值z来说，大于0时，表示增加趋势；小于0时，表示减少趋势；z的绝对值在大于1.96时，表示通过了置信度95％的显著性检验；在本实施例选择了z的绝对值大于1.96，即置信度为95％。当z>1.96时，表示该时序曲线呈显著上升趋势；当z<-1.96时，表示该时序曲线具有显著下降趋势；当z在正负1.96之间时为无趋势。

在本实施例中，步骤s6中所述的建立裸土变化识别方案的具体内容为：根据特征指标的变化趋势进行组合，将裸土指数和植被丰度均呈显著下降趋势的像元，判断为裸土增加；将裸土指数和植被丰度均呈显著上升趋势的像元，判断为裸土减少。

对于植被、不透水面、裸地和水体四种不同地物类型而言，裸土的裸土指数最低，其次为植被，再次为不透水面，而水体的裸土指数最高；植被丰度也存在明显差异。具体来说，植被的植被丰度最高，其次为不透水面、裸地，水体的植被丰度最低；裸土变化识别方案主要依据特征指标的数值在不同地物类型中存在明显的差异，因此，假设当地物类型发生变化时比如裸土变成植被，特征指标在年际水平上存在明显上升或者下降的变化趋势；因此，选取若干已知变化类型的点位，验证上述假设的合理性，并随机选取两个点位进行说明，第一个为裸土减少点位，第二个为裸土增加点位。对于第一个点位，当裸土减少时，基于植被近红外反射指数构建的植被丰度指标和裸土指数年均值呈现显著上升趋势如图3；对于第二个点位，当裸土增加时，基于植被近红外反射指数构建的植被丰度指标和裸土指数年均值呈现显著下降趋势如图4。依据实验结果，假设成立。

由于水体不稳定，易受包括降雨或冰川融化因素影响，所以不考虑水体，因此需对其进行掩膜处理；根据特征指标趋势组合的结果(图5)，采用研究时段开始年份和结束年份均值大于-0.05掩模掉水体发生变化的区域和变成水体的区域；最后得到裸土变化识别的方案：将裸土指数和植被丰度均呈显著下降趋势的像元，判断为裸土增加；将裸土指数和植被丰度均呈显著上升趋势的像元，判断为裸土减少。裸土变化识别方法技术流程图见图6。

以上上述流程，逐像元判断裸土是否发生变化，最终生成研究区裸土变化识别结果图。依据本实施例中提供的方法，以我国云南省为例，所获得的研究区2001-2018年裸土变化结果分布图见图7。

在本实施例中，在特征指标评估和选取时本实施例没有随意使用某一个特征指标，而是经过评估后才挑选的。特征指标评估和选取的步骤如下：

(a)选取若干不同地物类型的稳定点位，具体地物类型包括：植被、不透水面和裸土。

(b)根据不同地物类型点位分布位置提取不同特征指标的数值。

(c)分别评估基于植被近红外反射率指数和裸土指数的年最大值、年均值、50th分位值(植被丰度)等特征指标在不同地物类型中的分离度(图2)，并选择拥有最佳分离的特征指标。

在本实施例中，在建立裸土变化识别方法时。主要步骤如下：(a)根据特征指标评估时所表现出来的特点，特征指标在不同地物类型中存在明显的差异，因此，提出假设，当地物类型发生变化时比如裸土变成植被，特征指标在年际水平时序信号上存在明显上升或者下降的变化趋势。

(b)选取若干已知变化类型的点位，验证上述假设的合理性。

(c)在假设成立的基础上，继续统计裸土增加和裸土减少点位所满足的特征指标趋势组合的情况，并把90％以上点位均满足的规则作为本实施例裸土变化识别的方法。

依据本实施例提供的方法，以我国云南省为例，所获得的研究区2001-2018年裸土变化结果分布图见图7。结果显示，2001-2018裸土增加的区域相对较为集中，主要分布在大理、曲靖和九乡市；裸土减少的区域则比较分散，分布在云南省各个地方，没有相对较集中的区域。为了验证该结果的精度，随机挑选几处裸土发生变化的区域，与谷歌地图高分影像进行目视解译对比验证。通过谷歌历史影像与挑选的出来的区域进行

对比发现，该区域确实发生明显变化，与变化识别的结果一致。以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。