一种通过长时序卫星遥感评价湿地退化的方法

技术领域

本发明涉及湿地资源保护领域，尤其涉及一种通过长时序卫星遥感评价湿地退化的方法。

背景技术

湿地作为全球三大生态系统之一，在蓄水、调节河川径流、补给地下水以及为多种动植物提供栖息环境等方面发挥重要作用。目前湿地退化问题十分严重,不仅退化湿地面积大、类型多,而且受人类活动干扰强度大,并且退化过程和成因复杂。因此，明确湿地退化范围，掌握湿地退化机制，确定湿地退化程度,为退化湿地恢复提供前提保障，是当前湿地退化评估面临的重要问题。

目前湿地退化评价指标多侧重于小尺度湿地的生物、物理、化学指标，集中在水质、土壤和生物等单指标或者多指标的比较，研究范围小，不利于对整个生态系统退化进行系统性的综合评价。遥感技术由于其监测范围广、获取信息快、获取信息受限制条件少、获取信息量大、成本低的优势，已经广泛应用于湿地的提取和动态监测。但是当前遥感技术在研究湿地退化评估方面面临棘手挑战，大多方法利用退化前后两期或者多期土地利用数据建立转移矩阵来确定，土地利用数据间隔时期较长，无法反映湿地退化的演变序列和发展趋势。

因此针对上述问题，本发明提出一种通过湿地退化指数评价湿地退化程度的方法。该方法基于长时序卫星遥感数据集，根据每相邻两年(两期)的湿地水体频率信息提取湿地变化的时空特征，建立湿地退化指数，确定湿地退化程度。通过该方法，能够确定湿地退化范围，确定湿地退化程度从而为湿地资源的管理和保护提供依据，为其管理部门的决策提供一定的参考指导。

发明内容

本发明的目的为了利用长时序的水体数据获取水体的时空变化特征，从而通过湿地退化指数的计算确定湿地退化程度。

为了实现上述目的，本发明采用了如下技术方案：

一种通过长时序卫星遥感评价湿地退化的方法，包括以下步骤：

S1、获取覆盖研究区范围的中分辨率卫星遥感地表反射率数据；将地表反射率数据通过裁剪、镶嵌合成研究区范围的遥感影像；

S2、对地表反射率数据进行去云/云阴影、去冰/雪预处理，去除云/云阴影覆盖、冰、雪的影响；

S3、通过波段合成计算水体指数和植被指数，并通过多种指数阈值法以及逻辑判别法获取水体信息；首先根据绿、近红外、短波红外波段的地表反射率计算自动水体提取指数AWEI

上述式(1)一般适用于平原和城市区域，式(2)适用于山区；式中，

AWEI

AWEI

其中G、B、R、NIR、MIR、SWIR1和SWIR2分别表示绿波段(0.492～0.577μm)、蓝波段(0.44～0.475μm)、红波段(0.625～0.74μm)、近红外(0.7～2.5μm)、中红外(2.5～25μm)、短波红外1(1.56～1.66μm)和短波红外2(2.1～2.3μm)的地表反射率；

S4、利用坡度信息去除山体阴影对水体提取的影响；利用DEM数据生成坡度栅格数据，当上述S3判别为水体的像元所在坡度大于10°时，将该像元修改为非水体；

S5、精度验证，结合高分辨率影像进行精度验证；通过生成混淆矩阵来判断精度的高低，其中，混淆矩阵包括：整体精度、用户精度、产品精度和Kappa系数；

S6、逐年或逐期计算水体频率，2-3年为一期，生成水体频率栅格数据；水体频率反映的是像元在一年内(或一期内)被水体覆盖的频率，通过以下公式计算：

WF为水体频率，取值范围是0～1；N表示一年内(或一期内)的有效观测次数，即像元在卫星过境时不为冰、雪、云以及阴影，可以获取有效观测的次数；W表示像元被识别为水体的次数；水体频次越高，表明为水体的次数越多；

S7、对逐年(逐期)WF栅格数据进行重分类；按照WF的大小，逐像元进行阈值划分，分为非水体(WF≤0.25)、季节性水体(WF∈(0.25,0.75))永久性水体(WF≥0.75)；

S8、根据上述WF重分类数据的逐年(逐期)变化，提取水体时空动态变化特征，识别湿地退化范围，确定湿地退化类型；通过比较相邻两年(或两期)三种状态的转换，分别提取由水体转为非水体且持续为非水体，以及由永久性水体转为季节性水体且持续为季节性水体的像元；第一种转换类型的像元属于水体→非水体的退化类型；对于第二种转换类型的像元，进一步比较转换前后的水体频率，若像元满足如下计算条件时，视为退化，退化类型为永久性水体→季节性水体；

其中，∑WF

S9、对于上述S8识别的退化像元，确定退化时长，结合退化之前的频率，计算湿地退化指数，表征湿地退化程度；

对水体-非水体以及永久-季节两种退化类型，湿地退化指数的计算方式如下：

式中，WLI为湿地退化指数，∑WF

S10、计算归一化湿地退化指数；为方便进行计算和理解，对湿地退化指数在研究区范围内进行归一化，利用湿地退化指数和研究区WLI最大值的比值构建归一化湿地退化指数，计算公式如下：

WLINor为归一化湿地退化指数，WLIMax-指湿地退化指数的最大值。

优选的，所述S2-S10中的计算都可以利用软件：ENVI，Erdas、ArcGIS、Matlab、Python、C++、Google Earth Engine，通过湿地退化指数确定湿地退化程度。

与现有技术相比，本发明提供了一种通过长时序卫星遥感评价湿地退化的方法，具备以下有益效果：

1、本发明所提出的通过长时序卫星遥感评价湿地退化的方法，在数据上不同于以往采用数据时期间隔长且数据少的土地利用数据来识别退化湿地，而是充分利用卫星遥感影像的全部有效观测值生成长时序数据集，用该数据集来监测、识别退化湿地，评估湿地退化程度。长时序数据集能够长期监测湿地的时空动态变化情况，对于掌握湿地变化的规律能够提供更为详细的信息，有利于湿地退化的评价。

2、本发明所提出的方法，创建了一种新的评估湿地退化程度的湿地退化指数，该指数对湿地退化程度的评估不同于前人评估时只考虑退化湿地当前的状态，而是综合湿地长期的动态变化过程，创建相应的参数分别体现湿地退化前后的状态对湿地退化程度的影响，既能体现湿地退化前状态的差异性，又能表征湿地退化后状态的一致性，即湿地退化前的状态越好，退化后时长越长，退化程度越高，也就是能够更加客观地反映湿地的退化程度。

附图说明

以下附图将结合本法明的具体构思、结构以及实施例中的技术结果进行说明，以充分了解本发明的特征、目的和结果。

图1为本发明流程图；

图2为归一化湿地退化指数结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。

实施例：

一种通过长时序卫星遥感评价湿地退化的方法，包括以下步骤：

S1、获取覆盖研究区范围的中分辨率卫星遥感地表反射率数据。将地表反射率数据通过裁剪、镶嵌合成研究区范围的遥感影像。

S2、对地表反射率数据进行去云/云阴影、去冰/雪预处理，去除云/云阴影覆盖、冰、雪的影响。

S3、通过波段合成计算水体指数和植被指数，并通过多种指数阈值法以及逻辑判别法获取水体信息。首先根据绿、近红外、短波红外波段的地表反射率计算自动水体提取指数AWEI

上述式(1)一般适用于平原和城市区域，式(2)适用于山区。式中，

AWEI

AWEI

其中G、B、R、NIR、MIR、SWIR1和SWIR2分别表示绿波段(0.492～0.577μm)、蓝波段(0.44～0.475μm)、红波段(0.625～0.74μm)、近红外(0.7～2.5μm)、中红外(2.5～25μm)、短波红外1(1.56～1.66μm)和短波红外2(2.1～2.3μm)的地表反射率。

S4、利用坡度信息去除山体阴影对水体提取的影响。利用DEM数据生成坡度栅格数据，当上述S3判别为水体的像元所在坡度大于10°时，将该像元修改为非水体。

S5、精度验证，结合高分辨率影像进行精度验证。通过生成混淆矩阵来判断精度的高低，其中，混淆矩阵包括：整体精度、用户精度、产品精度和Kappa系数。

S6、逐年或逐期(2-3年为一期)计算水体频率，生成水体频率栅格数据。水体频率反映的是像元在一年内(或一期内)被水体覆盖的频率，通过以下公式计算：

WF为水体频率，取值范围是0～1；N表示一年内(或一期内)的有效观测次数，即像元在卫星过境时不为冰、雪、云以及阴影，可以获取有效观测的次数；W表示像元被识别为水体的次数。水体频次越高，表明为水体的次数越多。

S7、对逐年(逐期)WF栅格数据进行重分类。按照WF的大小，逐像元进行阈值划分，分为非水体(WF≤0.25)、季节性水体(WF∈(0.25,0.75))永久性水体(WF≥0.75)。

S8、根据上述WF重分类数据的逐年(逐期)变化，提取水体时空动态变化特征，识别湿地退化范围，确定湿地退化类型。通过比较相邻两年(或两期)三种状态的转换，分别提取由水体转为非水体且持续为非水体，以及由永久性水体转为季节性水体且持续为季节性水体的像元。第一种转换类型的像元属于水体→非水体的退化类型。对于第二种转换类型的像元，进一步比较转换前后的水体频率，若像元满足如下计算条件时，视为退化，退化类型为永久性水体→季节性水体。

其中，∑WF

S9、对于上述S8识别的退化像元，确定退化时长，结合退化之前的频率，计算湿地退化指数，表征湿地退化程度。

对水体-非水体以及永久-季节两种退化类型，湿地退化指数的计算方式如下：

式中，WLI为湿地退化指数，∑WF

S10、计算归一化湿地退化指数。为方便进行计算和理解，对湿地退化指数在研究区范围内进行归一化，利用湿地退化指数和研究区WLI最大值的比值构建归一化湿地退化指数，计算公式如下：

WLI

所述S2-S10中的计算都可以利用软件：ENVI，Erdas、ArcGIS、Matlab、Python、C++、Google Earth Engine等。

实施例2：

本方案以张家口坝上地区为研究区，利用本发明提出的一种通过长时序卫星遥感评价湿地退化的方法，评估张家口坝上地区的湿地退化程度，主要步骤如下：

(1)基于Google Earth Engine(GEE)平台，获取研究区在1985-2019年的Landsat5、7、8地表反射率数据，其中，共5422景地表反射率数据；并基于GEE平台裁剪、镶嵌地表反射率数据，合成研究区范围的遥感影像。

(2)对地表反射率数据进行去云/云阴影、去冰/雪预处理，去除云/云阴影覆盖、冰、雪的影响。

基于GEE云平台，利用Landsat质量评估波段QA波段，根据波段信息筛选并剔除含有云/云阴影、冰/雪的像元。

(3)通过波段合成计算水体指数和植被指数，并通过多种指数阈值法以及逻辑判别法获取水体信息。

直接利用GEE云平台对预处理之后的地表反射率数据进行计算，根据绿、蓝、红、近红外、中红外、短波红外波段的地表反射率计算自动水体提取指数AWEInsh和AWEIsh、改进的归一化差异水体指数MNDWI、归一化植被指数NDVI、以及增强型植被指数EVI。

基于GEE利用多种指数阈值法以及逻辑判别法提取出每景数据中的水体部分，根据判别式的适用范围，在本研究区采用的水体判别式如下：

Water

AWEI

AWEI

其中G、B、R、NIR、MIR、SWIR1和SWIR2分别表示蓝波段(0.44～0.475μm)、绿波段(0.492～0.577μm)、红波段(0.625～0.74μm)、近红外(0.7～2.5μm)、中红外(2.5～25μm)、短波红外1(1.56～1.66μm)和短波红外2(2.1～2.35μm)的地表反射率。在Landsat中对应的波段信息如下：

(4)利用坡度信息去除山体阴影对水体提取的影响。

首先利用ArcGIS的坡度工具，根据SRTM-3高程数据生成坡度栅格数据；然后将坡度栅格数据和水体提取结果数据叠加，当像元识别为水体且所在坡度大于10°时，利用ArcGIS的栅格计算器工具将该像元修改为非水体。

(5)对水体提取结果数据结合谷歌影像进行精度验证。

分别选取2000年、2006年、2010年和2016年谷歌地球高清影像，通过目视解译建立样本点，对同一时期基于Landsat数据的水体提取结果进行验证，生成混淆矩阵。

(6)逐期(2年为一期)计算水体频率，生成水体频率栅格数据。

直接利用GEE云平台，分别统计该像元在两年内的有效观测次数和被识别为水体的次数，逐像元计算水体频率值，计算公式如下：

WF为水体频率，取值范围是0～1；N表示一年内(或一期内)的有效观测次数，即像元在卫星过境时不为冰、雪、云以及阴影，可以获取有效观测的次数；W表示像元被识别为水体的次数。最终得到1985-2019年共18期水体频率WF栅格数据。

(7)水体频率WF栅格数据重分类

对逐期的WF栅格数据进行重分类。按照WF的大小，利用ArcGIS重分类工具逐像元进行阈值划分，分为非水体(WF<0.25)、季节性水体(WF∈(0.25,0.75))、永久性水体(WF>0.75)。

(8)分析WF重分类数据的逐期变化，提取水体时空动态变化特征，识别湿地退化范围，确定湿地退化类型。

退化类型的确定主要基于matlab计算实现，利用matlab根据相邻两期三种状态的转换，分别提取由水体转为非水体且持续为非水体，以及由永久性水体转为季节性水体且持续为季节性水体的像元。第一种转换类型的像元属于水体→非水体的退化类型。对于第二种转换类型的像元，进一步比较转换前后的水体频率，若像元满足如下计算条件时，视为退化，退化类型为永久性水体→季节性水体

∑WF

最终确定湿地退化类型及退化范围。

(9)根据退化像元，确定退化时长；结合退化之前的频率，计算湿地退化指数，表征湿地退化程度。

对水体-非水体以及永久-季节两种退化类型，分别利用ArcGIS的栅格计算器计算湿地退化指数，计算方式如下：

式中，WLI为湿地退化指数，∑WF

∑WF

最终得到湿地退化指数，确定湿地退化程度。

(10)计算归一化湿地退化指数。

利用ArcGIS合成两种退化类型的湿地退化指数的计算结果，并对湿地退化指数进行统计，得到研究区内湿地退化指数的最大值，并利用ArcGIS的栅格计算器计算每个像元与最大值的比值，得到归一化湿地退化指数。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。