一种基于DI的森林扰动变化遥感监测方法及装置与流程

一种基于di的森林扰动变化遥感监测方法及装置
技术领域
1.本发明属于遥感图像处理技术领域，具体涉及一种基于di的森林扰动变化遥感监测方法及装置。


背景技术：

2.森林作为最大的陆地生态系统，是地球生物圈的重要组成部分，也是人类赖以生存和发展的物质基础。森林在全球碳循环和气候调节扮演重要的角色，通过光合作用将大量co2固定起来，成为最大的陆地碳汇。人类活动、病虫害、雨雪灾害和火灾等人为和自然因素均会引起不同程度的森林扰动，进而影响森林碳汇。因此，监测森林扰动对于景观、区域甚至全球尺度的森林碳核算和气候变化至关重要。
3.传统的森林扰动及其变化监测主要依赖于野外调查，耗时费力，且难以在大区域及偏远山区开展作业。遥感凭借其宏观、快速、准确、实时等优势，森林扰动监测取得了大量成果，监测理论、方法及技术手段不断丰富和完善。过去20年使用遥感影像进行长时间序列分析主要基于modis和avhrr等有高时间分辨率、低空间分辨率和覆盖范围大的影像，但这类低空间分辨率影像开发的扰动产品无法准确监测小规模扰动，如森林采伐活动引起的扰动。
4.近年来，基于30m分辨率的landsat卫星遥感技术的发展为森林扰动的监测提供了数据支撑。healey等(2005)则采用30m分辨率landsat数据缨帽变换的亮度、绿度和湿度分量在扰动林和非扰动林中不同的特性构建了扰动指数(disturbance index,di)。该指数是一种基于缨帽变换的森林扰动监测方法，能够凸显非扰动森林(低亮度、高绿度和高湿度)和受扰动森林(高亮度、低绿度和低湿度)之间的差异，成功保留并增强森林扰动的相关信息。已有研究表明，di指数对于不同类型(如森林砍伐、火灾和虫害等)的扰动都具有良好的监测能力，模型简洁可靠且易于计算解释。huang等(2010)是根据高置信度森林像元的归一化红光波段和两个短波红外波段构建了landsat数据的综合森林指数(integrated forest z-score,ifz)。然而，目前大多数研究往往采用时间跨度大的两期或少数几期两两组合对同一区域影像进行对比分析。因此，如何通过快速有效的手段开展长时间序列的森林扰动变化分析是当前亟待解决的问题。


技术实现要素：

5.本发明的目的在于实现长时间序列大范围森林扰动变化遥感监测，结合变化矢量法、sen趋势度和mann-kendall方法，提供一种基于di的森林扰动变化遥感监测方法及装置。
6.为实现上述目的，本发明的技术方案是：一种基于di的森林扰动变化遥感监测方法，包括如下步骤：
7.s1、影像预处理
8.选取覆盖研究区的遥感图像，进行包括辐射校正和几何校正的影像预处理，获取
研究区地表反射率数据；
9.s2、基于遗传-随机森林特征筛选算法的林地信息提取
10.采用基于遗传-随机森林特征筛选算法进行研究区土地覆盖分类，根据分类结果掩膜非林地像元，获取林地信息；
11.s3、计算森林扰动指数di
12.对研究区涉及的包括landsat-5 tm和landsat-8 oli的影像进行缨帽变换，获取亮度、绿度和湿度分量，再采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，从而计算di指数；
13.s4、构建森林扰动变化趋势模型
14.结合变化矢量法、sen趋势度和mann-kendall方法，构建森林扰动变化趋势模型，以实现长时间序列研究区森林扰动变化趋势和特征分析。
15.在本发明一实施例中，基于遗传-随机森林特征筛选算法的林地信息提取具体包括如下步骤：
16.s21、制定研究区土地覆盖分类体系；
17.s22、基于遗传-随机森林特征筛选算法的研究区土地覆盖分类：首先，基于影像自带的光谱波段，引入包括水体指数、植被指数、建筑指数、裸土指数和dem数据的衍生波段共同进行土地覆盖分类；其次，通过遗传算法和随机森林相结合来进行特征筛选，获得目标地物之间最大区分性和最小冗余度的最优特征波段子集；最后，采用随机森林对影像进行分类；
18.s23、根据土地覆盖分类结果，掩膜非林地像元，获取林地信息；
19.s24、采用混淆矩阵进行林地分类结果精度评价。
20.在本发明一实施例中，所述步骤s22中，衍生波段包括：
21.(1)水体指数：归一化差值水体指数ndwi、改进型归一化差值水体指数mndwi、自动水体提取指数awei、水体指数wi2015；
22.(2)植被指数：归一化差异植被指数ndvi、土壤调节植被指数savi、修正型土壤调节植被指数msavi、归一化差值山地植被指数ndmvi、绿度指数greenness；
23.(3)建筑指数：归一化差值建筑指数ndbi、基于指数的建筑用地指数ibi；
24.(4)裸土指数：裸土指数bsi、归一化差值裸土指数ndbai；
25.(5)纹理特征波段：均值mea、方差var、均匀性hom、对比度con、非相似性di、熵ent、角二阶asm、相关度cor；
26.(6)dem数据。
27.在本发明一实施例中，所述步骤s3中，计算森林扰动指数di具体包括如下步骤：
28.s31、对研究区涉及的包括landsat-5 tm和landsat-8 oli的影像进行缨帽变换，获取亮度、绿度和湿度分量；缨帽变换的基本模型为：
29.u＝r
t
x r
30.其中，u为缨帽变换后增强的数据信息；r为缨帽变换系数；x为各波段的灰度值或反射率；r为偏移的常数；缨帽变换的第1、2和3分量分别为所需研究区的亮度、绿度和湿度分量；
31.s32、采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，消除影像
时相差异带来的误差；具体公式如下：
32.br＝(b-b
μ
)/b
σ
33.gr＝(g-g
μ
)/g
σ
34.wr＝(w-w
μ
)/w
σ
35.其中，b、g、w分别为缨帽变换的亮度、绿度和湿度值，br、gr和wr分别为进行标准化之后的亮度、绿度和湿度值；b
μ
、g
μ
和w
μ
分别为亮度、绿度和湿度的均值；b
σ
、g
σ
和w
σ
分别为亮度、绿度和湿度的标准差；
36.s33：在数据标准化的基础上，计算森林扰动指数di；具体公式如下：
37.di＝b
r-(gr wr)
38.其中，di值越高，代表森林受到扰动的程度越高；di值越低，表示森林越健康。
39.在本发明一实施例中，所述步骤s4中，构建森林扰动变化趋势模型具体包括如下步骤：
40.s41、采用变化矢量法计算di变化强度，通过先计算欧氏距离获得相邻两年间矢量的模，再累加矢量的模即可获得di变化强度；具体公式如下：
[0041][0042]
其中，p为di变化强度，n是序列年份数，x为森林扰动指数di，y为年份；
[0043]
s42、采用sen趋势度方法检验di的变化趋势；具体公式如下：
[0044][0045]
其中，x为森林扰动指数di，1《i《j《n，当β》0时表示扰动指数呈现上升趋势，即森林受到扰动的程度越来越高；反之呈下降趋势；
[0046]
s43：采用mann-kendall趋势检验法检验di变化趋势的显著性，可以考察时间序列的di数据是否存在持续上升或下降的趋势；具体公式如下：对于序列(x1,x2,
…
,xn)，先确定所有对偶值(xi,xj,j》i)中xi与xj的大小关系(设为s)；做如下假设：h0：序列中数据随机排列，即无显著趋势；h1：序列存在上升或下降单调趋势；检验统计量s由以下公式计算：
[0047][0048][0049]
根据时间序列长度n值大小的不同，显著性检验统计量的选取有所不同：
[0050]
当n≤10时，直接使用统计量s进行双边趋势检验；在给定显著性水平α下，若|s|≥s
α/2
，表明原序列存在显著趋势；且当s》0，则认为序列存在上升趋势；s＝0，无趋势；s《0，则认为序列存在下降趋势；
[0051]
当n》10时，统计量s近似服从标准正态分布，使用检验统计量zs进行趋势检验，其公式如下：
[0052][0053][0054]
其中，n是序列中数据个数，m是重复出现的数据组的个数，ti是第i组重复数据组中的重复数据个数；
[0055]
s44、采用自然断点法将di矢量变化强度实行图像分割，分别划分为无变化、低变化、中变化和高变化区域，定义低变化和中变化的分割阈值为k；结合矢量变化强度、sen趋势度和mann-kendall方法，提出了森林扰动变化趋势模型；具体公式如下：
[0056][0057]
其中，p为di变化强度，β为sen趋势度，s为mann-kendall趋势检验法中的统计量，s0为当研究期时间序列长度n和显著性水平α时，介于森林扰动趋势显著和不显著之间所对应的统计量临界值；平稳型表示在研究时期内森林未受到扰动的像元；持续变好型表示在研究时期内森林早期受到扰动，但后期重新生长并持续变好的像元；持续变差型表示在研究时期内森林持续受到扰动的像元；一般扰动型表示森林在某一时期受到干扰但扰动不会持续发生的像元。
[0058]
本发明还提供了一种基于di的森林扰动变化遥感监测装置，包括：数据输入输出模块、遥感影像预处理模块、林地信息提取模块、缨帽变换模块、森林扰动指数(di)模块、森林扰动变化趋势模块；所述数据输入输出模块连接遥感影像预处理模块，所述林地信息提取模块分别连接遥感影像预处理模块和缨帽变换模块，所述森林扰动指数(di)模块分别连接缨帽变换模块和森林扰动变化趋势模块；
[0059]
所述数据输入输出模块：用于读取研究区的遥感影像、研究区矢量边界，保存森林扰动变化趋势结果；
[0060]
所述遥感影像预处理模块：用于实现遥感影像的辐射校正、几何校正和影像裁剪等处理，消除影像辐射畸变和几何畸变等影响；
[0061]
所述林地信息提取模块：利用遗传-随机森林特征筛选算法实现土地覆盖分类，根据分类结果掩膜非林地像元，获取林地信息；
[0062]
所述缨帽变换模块：利用传感器对应的缨帽变换系数进行影像缨帽变换，以获取亮度、绿度和湿度分量；
[0063]
所述森林扰动指数(di)模块：采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，进而计算di指数，以实现不同年份的森林扰动程度反演；
[0064]
所述森林扰动变化趋势模块：结合变化矢量法、sen趋势度和mann-kendall方法，构建森林扰动变化趋势模型，实现长时间序列森林扰动变化趋势分析。
[0065]
在本发明一实施例中，在所述林地信息提取模块中，首先通过构建水体指数、植被指数、建筑指数、裸土指数和dem数据等22个衍生波段来与原始波段进行组合分类，以丰富初始光谱信息；再通过遗传算法和随机森林相结合来进行特征筛选，可删除部分冗余信息，降低数据集维数，获得目标地物之间最大区分性和最小冗余度的最优特征波段子集，提高分类处理效率；最后采用随机森林获得土地覆盖分类结果，通过掩膜非林地像元，从而获取准确的林地信息。
[0066]
在本发明一实施例中，所述森林扰动变化趋势模块中，采用变化矢量法计算di变化强度，利用自然断点法设置研究区di变化强度的分割阈值k，并结合sen趋势度和mann-kendall方法检验di变化趋势及其显著性，从而构建森林扰动变化趋势模型，得到准确的平稳型、持续变好型、持续变差型和一般扰动型森林分布图。
[0067]
相较于现有技术，本发明具有以下有益效果：
[0068]
(1)本发明采用变化矢量法、sen趋势度和mann-kendall方法相结合的方法，创新性地提出一种简易可行的森林扰动变化趋势模型，该模型能较好地实现研究区长时间序列的森林扰动变化趋势监测。
[0069]
(2)本发明通过引入一系列光谱特征指数、纹理特征和dem数据，建立了基于遗传-随机森林特征筛选算法的土地覆盖分类模型，改善了传统的随机森林分类存在地物混淆的问题，提高林地的提取精度，为后续森林扰动监测提供基础数据。
附图说明
[0070]
图1为本发明基于di的森林扰动变化遥感监测方法的流程示意图；
[0071]
图2为本发明基于di的森林扰动变化遥感监测装置的结构示意图；
[0072]
图3为本发明实施例研究区的区位示意图；
[0073]
图4为本发明实施例研究区各年份土地覆盖分类图；
[0074]
图5为本发明实施例研究区各年份林地提取结果；
[0075]
图6为本发明实施例研究区各年份森林扰动di指数空间分布图；
[0076]
图7为本发明实施例研究区长时间序列森林扰动变化趋势图；
[0077]
图8为本发明实施例研究区持续变差型(a)和一般扰动型(b)区域局部示意图。
具体实施方式
[0078]
下面结合附图，对本发明的技术方案进行具体说明。
[0079]
本发明一种基于di的森林扰动变化遥感监测方法，包括如下步骤：
[0080]
s1、影像预处理
[0081]
选取覆盖研究区的遥感图像，进行包括辐射校正和几何校正的影像预处理，获取研究区地表反射率数据；
[0082]
s2、基于遗传-随机森林特征筛选算法的林地信息提取
[0083]
采用基于遗传-随机森林特征筛选算法进行研究区土地覆盖分类，根据分类结果掩膜非林地像元，获取林地信息；
[0084]
s3、计算森林扰动指数di
[0085]
对研究区涉及的包括landsat-5 tm和landsat-8 oli的影像进行缨帽变换，获取亮度、绿度和湿度分量，再采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，从而计算di指数；
[0086]
s4、构建森林扰动变化趋势模型
[0087]
结合变化矢量法、sen趋势度和mann-kendall方法，构建森林扰动变化趋势模型，以实现长时间序列研究区森林扰动变化趋势和特征分析。
[0088]
本发明还提供了一种基于di的森林扰动变化遥感监测装置，包括：数据输入输出模块、遥感影像预处理模块、林地信息提取模块、缨帽变换模块、森林扰动指数(di)模块、森林扰动变化趋势模块；所述数据输入输出模块连接遥感影像预处理模块，所述林地信息提取模块分别连接遥感影像预处理模块和缨帽变换模块，所述森林扰动指数(di)模块分别连接缨帽变换模块和森林扰动变化趋势模块；
[0089]
所述数据输入输出模块：用于读取研究区的遥感影像、研究区矢量边界，保存森林扰动变化趋势结果；
[0090]
所述遥感影像预处理模块：用于实现遥感影像的辐射校正、几何校正和影像裁剪等处理，消除影像辐射畸变和几何畸变等影响；
[0091]
所述林地信息提取模块：利用遗传-随机森林特征筛选算法实现土地覆盖分类，根据分类结果掩膜非林地像元，获取林地信息；
[0092]
所述缨帽变换模块：利用传感器对应的缨帽变换系数进行影像缨帽变换，以获取亮度、绿度和湿度分量；
[0093]
所述森林扰动指数(di)模块：采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，进而计算di指数，以实现不同年份的森林扰动程度反演；
[0094]
所述森林扰动变化趋势模块：结合变化矢量法、sen趋势度和mann-kendall方法，构建森林扰动变化趋势模型，实现长时间序列森林扰动变化趋势分析。
[0095]
以下为本发明具体实现过程。
[0096]
请参照图1-2，本发明提供了一种基于di的森林扰动变化遥感监测方法及装置。其中，一种基于di的森林扰动变化遥感监测方法(图1)，包括以下步骤：
[0097]
本实施例的研究区为：敖江流域，敖江流域地处东经118
°
50
′
至119
°
40
′
，北纬26
°
08
′
至26
°
26
′
，位于福建省东北部、福州市与宁德市交界处，是福建省重点流域之一。流域面积2655km2，涉及福州市的晋安区、连江县、罗源县、闽侯县、闽清县和宁德市古田县，如图3所示。林地作为敖江流域最主要的土地覆盖类型，占流域总面积的78％以上，及时、准确地掌握敖江流域森林扰动的时空变化特征对于流域森林管理和生态保护意义深远。
[0098]
s1：影像预处理
[0099]
本实施例选取敖江流域2004、2009、2013、2015和2018年5个年份的landsat-5 tm、landsat-8 oli多光谱影像，进行辐射校正和几何校正等影像预处理，获取后续处理所需的研究区地表反射率数据，以消除地球大气条件、太阳高度角、地形以及遥感传感器本身误差所带来的几何变形和光谱失真等影响。
[0100]
对于lansat-5影像，通过引进chavez的cost模型进行大气校正，即可获得iacm模型，公式为：
[0101][0102]
l
λ
＝gain
λ
·qλ
bias
λ
[0103]
lh＝gain
λ
·hλ
bias
λ
[0104]
其中，ρ
λ-iacm
为波段λ经iacm模型辐射校正后的地表反射率；l
λ
为波段λ在传感器处光谱辐射值(w
·
m-2
·
sr-1
·
μm-1
)；lh为大气修正参数(w
·
m-2
·
sr-1
·
μm-1
)；d为日-地天文单位距离；esun
λ
为波段λ在大气顶部的平均太阳辐照度(w
·
m-2
·
μm-1
)；θz为太阳天顶角(
°
)；q
λ
为波段λ的像元亮度值；h
λ
为波段λ最暗像元的亮度值；gain
λ
和bias
λ
分别为波段λ的增益值和偏置值(w
·
m-2
·
sr-1
·
μm-1
)。
[0105]
对于lansat-8影像，同样采用cost模型对landsat-8进行大气校正，公式为：
[0106][0107]
其中，ρ
λ-iacm
为波段λ经iacm模型辐射校正后的地表反射率；q
λ
为波段λ的像元亮度值；m
ρ
和a
ρ
分别为波段反射率的调整因子和调整参数，可从影像头文件中获取；h
λ
为波段λ最暗像元的亮度值；θz为太阳天顶角(
°
)。
[0108]
此外，由于从usgs网站下载的landsat数据产品均是l1t级，已经过地面控制点的几何校正和系统辐射校正，且通过dem数据的地形校正处理，因此，可直接用于landsat系列影像之间的叠加分析。
[0109]
s2：基于遗传-随机森林特征筛选算法的林地信息提取
[0110]
在本实施例中，采用基于遗传-随机森林特征筛选算法进行敖江流域土地覆盖分类，根据分类结果掩膜非林地像元，获取所需的林地像元，为后续敖江流域森林扰动监测提供基础信息。
[0111]
具体包括如下步骤：
[0112]
s21：根据敖江流域的地物特征以及所用的遥感影像分辨率，制定研究区土地覆盖分类体系，将其分为林地、耕地、草地、水体、建筑用地、工矿用地、裸土和火烧迹地这8大土地覆盖类型。
[0113]
s22：基于遗传-随机森林特征筛选算法的敖江流域土地覆盖分类：首先，基于影像自带的光谱波段，引入一系列的水体指数、植被指数、建筑指数、裸土指数和dem数据等22个衍生波段共同进行土地覆盖分类；其次，通过遗传算法和随机森林相结合来进行特征筛选，获得目标地物之间最大区分性和最小冗余度的最优特征波段子集；最后，采用随机森林对影像进行分类，得到5个研究年份的土地覆盖分类图，如图4所示。
[0114]
s23：根据土地覆盖分类结果，掩膜非林地像元，获取敖江流域林地信息，如图5所示。
[0115]
s24：采用混淆矩阵进行林地分类结果精度评价，敖江流域各年份的精度验证结果如表1所示。
[0116]
从表1可以看出，基于遗传-随机森林特征筛选算法的土地覆盖分类模型的分类结果较好，其获得的5个年份的林地分类总精度均在95％以上，kappa系数均高于0.9，表明本发明所提取的林地信息较为准确，满足研究的需要。
[0117]
表1敖江流域5个年份林地精度验证结果
[0118][0119]
s3：计算森林扰动指数(di)
[0120]
在本实施例中，对敖江流域涉及的landsat-5 tm和landsat-8 oli等影像进行缨帽变换，获取亮度、绿度和湿度分量，再采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，从而计算di指数。具体包括如下步骤：
[0121]
s31：对敖江流域涉及的landsat-5 tm和landsat-8 oli等影像进行缨帽变换，获取亮度、绿度和湿度分量。
[0122]
缨帽变换的基本模型为：
[0123]
u＝r
t
x r
[0124]
其中，u为缨帽变换后增强的数据信息；r为缨帽变换系数；x为各波段的灰度值或反射率；r为偏移的常数。缨帽变换的第1、2和3分量分别为所需研究区的亮度、绿度和湿度分量。
[0125]
由于缨帽变换依赖于传感器本身的波段设置和特性，因此不同传感器的缨帽变换系数不同。在本实施例中，landsat-5 tm影像所对应的缨帽变换亮度brightness、绿度greenness和湿度wetness公式如下：
[0126]
brightness＝0.2043
×
ρ1 0.4158
×
ρ2 0.5524
×
ρ3 0.5741
×
ρ4 0.3124
×
ρ5 0.2303
×
ρ6[0127]
greenness＝-0.1603
×
ρ
1-0.2819
×
ρ
2-0.4934
×
ρ3 0.7940
×
ρ
4-0.0002
×
ρ
5-0.1446
×
ρ6[0128]
wetness＝0.0315
×
ρ1 0.2021
×
ρ2 0.3102
×
ρ3 0.1594
×
ρ
4-0.6806
×
ρ
5-0.6109
×
ρ6[0129]
其中，波段ρ1、ρ2、ρ3、ρ4、ρ5和ρ6分别为landsat-5影像第1、2、3、4、5、7波段(即：蓝光、绿光、红光、近红外、中红外1和中红外2波段)的表观反射率。
[0130]
landsat-8oli影像所对应的缨帽变换亮度brightness、绿度greenness和湿度wetness公式如下：
[0131]
brightness＝0.3029
×
ρ1 0.2786
×
ρ2 0.4733
×
ρ3 0.5599
×
ρ4 0.508
×
ρ5 0.1872
×
ρ6[0132]
greenness＝-0.2941
×
ρ
1-0.243
×
ρ
2-0.5424
×
ρ3 0.7276
×
ρ4 0.0713
×
ρ
5-0.1608
×
ρ6[0133]
wetness＝0.1511
×
ρ1 0.1973
×
ρ2 0.3283
×
ρ3 0.3407
×
ρ
4-0.7117
×
ρ
5-0.4559
×
ρ6[0134]
其中，波段ρ1、ρ2、ρ3、ρ4、ρ5和ρ6分别为landsat-8影像第2、3、4、5、6、7波段(即：蓝光、绿光、红光、近红外、中红外1和中红外2波段)的表观反射率。
[0135]
s32：采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，消除影像时相差异带来的误差，保证不同年份的反演结果具有可比性。具体公式如下：
[0136]br
＝(b-b
μ
)/b
σ
[0137]gr
＝(g-g
μ
)/g
σ
[0138]
wr＝(w-w
μ
)/w
σ
[0139]
其中，b、g、w分别为缨帽变换的亮度、绿度和湿度值，br、gr和wr分别为进行标准化之后的亮度、绿度和湿度值；b
μ
、g
μ
和w
μ
分别为亮度、绿度和湿度的均值；b
σ
、g
σ
和w
σ
分别为亮度、绿度和湿度的标准差。
[0140]
s33：在数据标准化的基础上，计算di指数。具体公式如下：
[0141]
di＝b
r-(gr wr)
[0142]
其中，di值越高，代表森林受到扰动的程度越高；di值越低，表示森林越健康。
[0143]
图6为敖江流域5个研究年份的森林扰动di指数空间分布图。图中浅灰色代表di值较低的区域，而黑色代表di值较高的区域。可以发现：森林受到扰动的程度越高的区域，如火烧迹地和建筑用地，它们的di值均较高；而未受到扰动的区域则di值均较低，说明所反演的di值可以客观地反映森林扰动的实际情况。
[0144]
s4：构建森林扰动变化趋势模型
[0145]
在本实施例中，结合变化矢量法、sen趋势度和mann-kendall方法，构建一种简易可行的敖江流域森林扰动变化趋势模型，以实现长时间序列敖江流域森林扰动变化趋势和特征分析。
[0146]
具体包括如下步骤：
[0147]
s41：采用变化矢量法计算di变化强度。本实施例以年为间隔，计算5个年份di影像，定义二维矢量空间，通过先计算欧氏距离获得相邻两年间矢量的模，再累加矢量的模即可获得di变化强度。具体公式如下：
[0148][0149]
其中，p为di变化强度，n是序列年份数，x为扰动指数di，y为年份。
[0150]
s42：采用sen趋势度方法检验di的变化趋势。具体公式如下：
[0151][0152]
其中，x为扰动指数di，1《i《j《n，当β》0时表示扰动指数呈现上升趋势，即森林受到扰动的程度越来越高；反之呈下降趋势。
[0153]
s43：采用mann-kendall趋势检验法检验di变化趋势的显著性，可以考察时间序列的di数据是否存在持续上升或下降的趋势。
[0154]
具体公式如下：对于序列(x1,x2,
…
,xn)，先确定所有对偶值(xi,xj,j》i)中xi与xj的大小关系(设为s)。做如下假设：h0：序列中数据随机排列，即无显著趋势；h1：序列存在上升或下降单调趋势。检验统计量s由以下公式计算：
[0155]
[0156][0157]
根据时间序列长度n值大小的不同，显著性检验统计量的选取有所不同：
[0158]
当n≤10时，直接使用统计量s进行双边趋势检验。在给定显著性水平α下，若|s|≥s
α/2
，表明原序列存在显著趋势。且当s》0，则认为序列存在上升趋势；s＝0，无趋势；s《0，则认为序列存在下降趋势。
[0159]
当n》10时，统计量s近似服从标准正态分布，使用检验统计量zs进行趋势检验，其公式如下：
[0160][0161][0162]
其中，n是序列中数据个数，m是重复出现的数据组的个数，ti是第i组重复数据组中的重复数据个数。
[0163]
在本实施例中，时间序列长度n为5，所以采用检验统计量s来进行趋势检验，检验中取显著性水平α＝0.05，对应的s＝8。当β》0且|s|＞8，表明di值呈显著上升趋势；当β》0且|s|≤8，表明di值呈上升但不显著趋势；同理，当β《0且|s|＞8，表明di值呈显著下降趋势；当β《0且|s|≤8，表明di值呈下降但不显著趋势。
[0164]
s44：采用自然断点法将di矢量变化强度实行图像分割，分别划分为无变化、低变化、中变化和高变化区域，定义低变化和中变化的分割阈值为k。结合矢量变化强度、sen趋势度和mann-kendall方法，提出了森林扰动变化趋势模型。具体公式如下：
[0165][0166]
其中，p为di变化强度，β为sen趋势度，s为mann-kendall趋势检验法中的统计量，s0为当研究期时间序列长度n和显著性水平α时，介于森林扰动趋势显著和不显著之间所对应的统计量临界值。平稳型表示在研究时期内森林未受到扰动的像元；持续变好型表示在研究时期内森林早期受到扰动，但后期重新生长并持续变好的像元；持续变差型表示在研究时期内森林持续受到扰动的像元；一般扰动型表示森林在某一时期受到干扰但扰动不会持续发生的像元。在本实施例中，分割阈值k＝17.75；当时间序列长度n为5、显著性水平α＝0.05时，s0＝8。
[0167]
图7为敖江流域森林扰动变化趋势图。图中中灰色区域代表平稳型，是指长时间序列未受到扰动的森林，而持续变好型(黑色区域)、持续变差型(深灰色区域)和一般扰动型(浅灰色区域)都是指长时间序列受到扰动的森林。其中持续变好型表示长时间序列森林早期受到扰动，但后期重新生长并持续变好的像元，呈现利好趋势，因此未来无需过多关注。持续变差型是指长时间序列森林持续受到扰动的像元，表明其未来发生扰动的趋势明显，应该重点关注。如图8(a)黑色框内所示，2004年该区域森林长势较好，但逐年衰退，到了2018年森林长势欠佳，在影像中表现为红中泛白(影像rgb:543)。一般扰动型是指某一时期森林受到干扰而扰动不会持续发生的像元，但未来再次发生扰动的可能性较大，表明其生态也较脆弱。如图8(b)黑色框内所示，该区域在2009年发生大规模火灾，但在2013年森林得以恢复生长；在2015年该区域的西北方向上再次发生小面积的森林火灾，随后在2018年又得以恢复，这类区域的森林扰动属于动态波动，也是容易受到干扰的脆弱区，未来也应该重点保护。因此，将持续变差型和一般扰动型这两类容易受到扰动的生态脆弱区作为本研究的重点关注对象，以期为未来流域森林管理和政策制定提供技术支撑。
[0168]
此外，本发明还提供了一种基于di的森林扰动变化遥感监测装置(图2)，该装置包括：数据输入输出模块、遥感影像预处理模块、林地信息提取模块、缨帽变换模块、森林扰动指数(di)模块、森林扰动变化趋势模块；所述数据输入输出模块连接遥感影像预处理模块，所述林地信息提取模块分别连接遥感影像预处理模块和缨帽变换模块，所述森林扰动指数(di)模块分别连接缨帽变换模块和森林扰动变化趋势模块。
[0169]
所述数据输入输出模块：用于读取研究区的遥感影像、研究区矢量边界，保存森林扰动变化趋势结果；在本实施例中，所述的遥感影像为敖江流域2004、2009、2013、2015和2018年5个年份的landsat-5 tm、landsat-8 oli多光谱影像。
[0170]
所述遥感影像预处理模块：用于实现遥感影像的辐射校正、几何校正和影像裁剪等处理，消除影像辐射畸变和几何畸变等影响；
[0171]
所述林地信息提取模块：利用遗传-随机森林特征筛选算法实现土地覆盖分类，根据分类结果掩膜非林地像元，获取林地信息；
[0172]
所述缨帽变换模块：利用传感器对应的缨帽变换系数进行影像缨帽变换，以获取亮度、绿度和湿度分量；在本实施例中，所述的传感器为landsat tm和oli传感器。
[0173]
所述森林扰动指数(di)模块：采用z-scores方法对所得到的亮度、绿度和湿度分量进行标准化，进而计算di指数，以实现不同年份的森林扰动程度反演；
[0174]
所述森林扰动变化趋势模块：结合变化矢量法、sen趋势度和mann-kendall方法，构建森林扰动变化趋势模型，实现长时间序列森林扰动变化趋势分析。
[0175]
在本实施例中的林地信息提取模块，首先通过构建水体指数、植被指数、建筑指数、裸土指数和dem数据等22个衍生波段来与原始波段进行组合分类，以丰富初始光谱信息；再通过遗传算法和随机森林相结合来进行特征筛选，可删除部分冗余信息，降低数据集维数，获得目标地物之间最大区分性和最小冗余度的最优特征波段子集，提高分类处理效率；最后采用随机森林获得土地覆盖分类结果，通过掩膜非林地像元，从而获取准确的林地信息。
[0176]
在本实施例中的森林扰动变化趋势模块，采用变化矢量法计算di变化强度，利用自然断点法设置敖江流域di变化强度的分割阈值k，并结合sen趋势度和mann-kendall方法
检验di变化趋势及其显著性，从而构建敖江流域森林扰动变化趋势模型，得到敖江流域准确的平稳型、持续变好型、持续变差型和一般扰动型森林分布图。
[0177]
通过以上实施例的描述，本领域的技术人员可以清楚地了解到本发明可以通过硬件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，本发明的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是cd-rom，u盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明实施例所述的方法。
[0178]
本领域技术人员可以理解附图只是一个优选实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。
[0179]
以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。