CN104851087B - 多尺度森林动态变化监测方法 - Google Patents

Abstract

本发明公开了一种多尺度森林动态变化监测方法，包括：遥感影像几何校正与配准；利用一年内不同季节的1KM分辨率的MODIS时间序列NDVI数据，获得1KM分辨率的土地覆盖类型图；利用多年1KM分辨率的土地覆盖图，生成粗尺度土地覆盖变化图；利用粗尺度土地覆盖变化类型图建立粗尺度森林植被变化掩模文件；根据粗尺度森林植被变化掩模文件，在30m分辨率TM影像上构建森林植被特征指数；粗尺度土地覆盖变化图与时间序列的森林特征指数相结合的森林动态变化信息提取。本发明提供的方法利用不同空间分辨率的时间序列遥感数据，可对大面积区域实现从粗尺度土地覆盖到精尺度森林变化的逐级细化的森林动态监测，不仅能提供监测效率，而且还能提高监测准确性。

Description

多尺度森林动态变化监测方法

技术领域

本发明涉及一种多尺度森林动态变化监测方法，特别涉及多尺度时间序列遥感数据在森林资源动态变化监测方面的技术，属于遥感技术在林业中的应用领域。

背景技术

森林植被不仅为人类的生存发展提供了最基本的物质和能量来源，而且在维持生态过程和生态平衡中也具有重要的作用，特别是保护脆弱的生态系统、水域和淡水资源方面起着关键的作用。森林自然的更新演替、森林灾害(火灾，病虫害)以及人类活动的干扰会造成森林植被发生变化，森林植被的动态直接影响到陆地生态系统的平衡和健康。森林动态变化的研究，可以及时了解森林资源的数量和质量，掌握森林资源的消长变化规律和趋势，分析影响与制约森林生长的自然、经济、社会客观条件，建立或更新森林资源档案，对揭示生态系统环境变化及植被恢复和重建布局等具有重要意义。

遥感技术为森林动态监测提供了新的手段。与传统的森林资源清查相比，遥感能获得大量不同空间分辨率、多时相的数据，为在森林变化监测提供了丰富的数据源。在全球或区域范围的大尺度森林监测中，中分辨率成像光谱设备(MODIS)等低空间分辨率的影像被大量的用在森林的变化监测中，这些数据具有高时间分辨率的特性，可以获得更多时间点的地表覆盖信息。但是这些数据由于空间分辨率较低，主要被用来监测大面积的土地利用类型的变化，更多的反应了森林和非森林的变化情况，对变化相对较小的区域很难监测。高分辨率的遥感影像具有更高的空间分辨率，为更精细的森林变化检测提供了很好的数据，但是绝大多数的高分辨率遥感影像需要购买如QuickBird、SPOT，数据成本较高。中等分辨率的影像数据如LandSat卫星的TM/ETM+影像免费，同时也具有较好的光谱与空间分辨率可以很好的用于变化检测。但是由于云等天气的影响，能够获得的有效时间序列数据质量不高。

在动态变化监测算法方面，有许多学者研究了针对各种不同空间分辨率的遥感影像提出了许多土地覆盖和植被的变化检测算法。这些变化检测一般对两个不同时期的影像进行分析，提取变化信息，算法可大致分为两类：一类是采用非监督的方式，对两个不同时期的影像，利用差值、比值、典型相关变换等方法构建差异影像，然后设置变化阈值，确定变化和非变化区域。这种变化检测算法的关键在于如何确定变化的阈值。另外一类变化检测算法主要采用监督分类的方式，对两个不同时期的影像分别分类，然后比较分类结果确定变化区域及变化类型。这种变化检测方法需要大量的样本。这些变化检测算法在土地覆盖变化监测方面取得了较好的效果，但是在应用到森林的变化检测中存在一定的局限性。由于植被的生长具有季节性规律，不同植被类型在不同季节的影像上会表现出明显的差异，这种植被季相性的差异会影响到变化检测的结果。

综上所述，现有的变化检测算法大都只利用了两个时期的影像，且只在相同空间分辨率的影像上进行变化检测，监测的范围和效率相对较低。此外，现有技术还存在以下突出的缺陷和问题，包括：

(1)只利用了单一空间分辨率的影像，在数据处理的效率和空间范围上有一定的局限性。而一般森林覆盖区域面积大，区域范围广，需要能够及时快速的定位到发生变化区域。

(2)只利用了两个时期的影像，而植被本身的生长就具有季节性规律，两个单一时相的数据必须要求在同一生长季，限制了算法的应用范围。

(3)只两个时间点的影像数据只能反映在这时间段内的森林状态信息，无法反映森林随时间的动态变化的规律。

发明内容

针对现有技术存在的问题，本发明提供了一种基于多尺度时间序列遥感影像的森林资源动态监测方法。

为了解决上述技术问题，本发明采用如下的技术方案：

本发明提供一种多尺度森林动态变化监测方法，主要包括以下步骤：

步骤一、影像预处理；

步骤二、获取土地覆盖类型图；

步骤三、生成森林植被变化信息；

步骤四、云及阴影掩模处理；

步骤五、构建森林植被特征指数；

步骤六、提取动态变化信息。

具体包括以下步骤：

步骤一、获取目标对象的MODIS-NDVI时间序列数据和对应的TM影像时间序列数据，并对影像进行几何校正与配准；

步骤二、利用不同的MODIS时间序列NDVI数据，获得土地覆盖类型图；

步骤三、利用土地覆盖类型图，生成粗尺度森植被变化信息；

步骤四、对TM影像进行云及阴影掩模处理；

步骤五、根据粗尺度森林植被变化掩模文件，在TM影像上构建森林植被特征指数；

步骤六、粗尺度土地覆盖变化图与时间序列的森林特征指数相结合的森林动态变化信息提取。

优选的，上述步骤一的影像几何校正与配准包括：

步骤1.1 MODIS-NDVI时间序列的影像几何校正；

步骤1.2 TM影像时间序列的影像几何校正；

步骤1.3 MODIS-NDVI与TM影像间的相互配准；

所述TM影像间的匹配根据TM自带的UTM坐标进行匹配，如果匹配误差较大在选择控制点，利用多项式纠正方法纠正；

所述MODIS-NDVI数据与TM数据间直接利用数据自带的地理坐标确定位置。

优选的，上述步骤二具体为利用一年内不同季节的1KM分辨率的MODIS时间序列NDVI数据，获得1KM分辨率的土地覆盖类型图，在1KM分辨率粗尺度数据上，将地表覆盖类型分为森林、农田、城市建筑、水体四大基本类型，其NDVI值由高到低依次为森林、农田、城市建筑、水体，通过统计分析不同土地覆盖类别各月份的NDVI直方图，确定对应月份各类别NDVI阈值，具体处理过程包括以下步骤：

步骤2.1各地表覆盖类别阈值确定；选择各种覆盖类型的样本，计算各个月份中的均值和方差，各个类别阈值区间被定义为其中j是月份，是类别i所选取样本的NDVI平均值，S_i是类别i所选取样本的NDVI方差；

步骤2.2基于决策树的粗尺度土地覆盖分类；四大基本土地覆盖类型的NDVI值的高低依次为森林(赋值1)、农田(赋值2)、城市建筑(赋值3)、水体(赋值4)。

优选的，上述步骤三具体为利用多年1KM分辨率的土地覆盖图，提取粗尺度森林植被变化信息，包括以下步骤：

步骤3.1对多年的1KM分辨率土地覆盖图逐年进行比较，分别计算相邻时间的两幅分类影像差值，像元计算结果为0就代表该像元类别未发生改变，而计算结果非0则表示该像元的土地覆盖类型在监测期间发生了改变；

步骤3.2将差值结果中值为非0的像素，且该像素在任意土地覆盖图中为森林类型(像素值为1)的所有像素提取出来；

步骤3.3利用该信息建立森林植被变化掩模文件，值为1区域表示森林变化，非0表示森林没有变化。

优选的，上述步骤四具体为根据粗尺度森林植被变化掩模文件，选择30M空间分辨率的TM影像，并进行云及阴影掩模处理，包括以下步骤：

步骤4.1云的识别；

第一步，处理1，通过波段2-波段6这五个波段的影像来识别云光谱，判断暖云、冷云、模糊像元、雪以及实现无云的掩膜；

第二步，波段6的云标记，将第一步已经确定是云的像元标记在波段6中。

第三步，处理2，剩余的模糊像元类别通过波段6再进行一次分类，在模糊像元中能区分一部分冷云和暖云，剩余不能识别的还是留在模糊类别中；

第四步，云像元聚合，将处理1和处理2中识别出来的云聚合在一起就是得到的云的结果；

步骤4.2云阴影的识别；通过地理位置来判别和确定来进行云阴影的识别，

第一步，反推云的高度，温度在高程上以1℃/hm递减，波段6转换成温度以后结合DEM反推出云的高度；

第二步，推测云影的位置，结合太阳高度角和第一步中云的高度反推出云像元对应的该像元对应云影的位置，对云团对应云影的位置区域进行5个像元(合理阈值范围)的缓冲；

第三步，确定云影，推测云影的缓冲区像元在波段4和波段5的灰度值都很低的位置被标记为云影。

优选的，上述步骤五具体为根据粗尺度森林植被变化掩模文件，在30m分辨率TM影像上构建森林植被特征指数，包括以下步骤：

步骤5.1森林样本自动选择；

第一步，将整个影像区域按照300*300像素，划分为若干小区域；

第二步，在每个小区域内，计算NDVI(归一化植被指数)，剔除NDVI小于0.2的像素，可剔除非植被暗目标；

第三步，在小区域内统计掩膜红光波段(TM3)的直方图，统计第二步中剩余的像元在红光波段的反射和频数的关系，得到反射-频数分步曲线，对曲线进行中值滤波得到统计直方图；

第四步，确定直方图中峰的位置，以频数达到1％的灰度值作为DNred(min)，当直方图中只有一个峰时，峰值对应的灰度值为DNred(max)；当该直方图有多个峰时，而且后面出现峰离第一个峰很近的时候，也应该被识别为森林峰，DNred(max)是最后一个森林峰对应的灰度值；

第五步，得到森林样本掩膜文件。根据阈值[DNred(min)，DNred(max)]，掩膜红光波段(TM3)得到森林样本，其中森林样本像元值为1，其它像元值为0；

步骤5.2森林特征指数计算；

第一步，根据选择的森林样本，计算波段i(i＝3，5，7)上光谱值的平均值和标准差，分别记为M_i和S_i；

第二步，计算每个波段的特征指数FZ_i，对任意像元在该波段的值为b_i，那么这个像元的森林特征指数(FZ_i)可以用公式(1)计算：

其中，公式(1)中FZ_i的实质是描述了任意像元与森林的偏离程度，该指数越小，对应像元是森林的可能性就越大；

对于具有多个波段的TM影像而言，对所有波段的FZ_i值做如公式(2)的积分运算能得到该像元的综合性森林特征指数(IFZ)；

其中NB是所使用的波段总数。

优选的，上述步骤六具体为粗尺度森林植被变化图与时间序列的森林特征指数相结合的森林动态变化信息提取，以时间序列的IFZ值与粗尺度的森林植被变化图采用决策树方法提取森林变化信息，通过以下方法进行：

设定IFZ阈值，结合粗检测结果确定，确定每个像素的时间序列IFZ值与阈值的比较结果：

如果时间序列的IFZ值都低于给定的阈值，则该像元在监测期间一直是典型的森林的像元；

如果时间序列IFZ均值大于给定阈值，或者绝大多数的IFZ值大于给定的阈值(＞90％)，则该像元被确定为非森林植被；

如果时间序列的IFZ值小部分(＜15％)高于给定的阈值，且这种高于阈值的时间不能连续超过2年，则在监测期间是疏林地或者非典型森林像元；

如果在前面步骤中没有被分出，且属于粗尺度的植被变化区域的像元，则归于森林变化发生区域。

优选的，上述森林变化有三种情况：

a类别为不可逆转的毁林，即森林变为非森林，且在监测期内未变成森林；

b类别为造林，即表现为非森林变为森林；

c类别为森林干扰，包括采伐更新、自然灾害、病虫害和火灾；

其中，a类别表现为IFZ连续很低，突然增高到阈值以上，此后就不再降低或者在后来的某个年份或者或某几个不连续的年份由于云和阴影的偶然因素低于阈值。

b类别表现为IFZ连续很高，随后慢慢的降低到到阈值以下，此后的时间内就不再高于阈值或在后来的某个年份或者或某几个不连续的年份偶然高于阈值；

c类别表现，森林干扰虽然通过IFZ突发的增高易于识别，但森林恢复过程中IFZ不会迅速地从较高的IFZ值降到不受干扰水平的IFZ值，而是表现为逐渐减少。

相对于现有技术，本发明提供的多尺度森林动态变化监测方法有以下优点：

(1)本发明充分利用了不同空间分辨率的遥感影像对不同土地覆盖及植被的响应特点，构建了有粗到细的变化检测技术，提高了植被监测效率；

(2)本发明针对植被随时间生长变化规律的特点，构建了基于时间序列的植被变化提取方法，提高监测准确性；

(3)本发明利用时间序列的数据，不仅可以监测变化，还可以监测提取森林扰动信息。

附图说明

图1为本发明流程示意图；

图2为土地覆盖类型分类决策流程示意图；

图3为TM数据通过ACCA除云流程示意图；

图4为云影识别过程示意图；

图5为VCT算法识别森林变化的决策流程示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及具体实施方式对本发明作进一步的详细描述。

(S1)遥感影像几何校正与配准；

(S2)利用一年内不同季节的1KM分辨率的MODIS时间序列NDVI数据，获得1KM分辨率的土地覆盖类型图；

(S3)利用多年1KM分辨率的土地覆盖图，生成粗尺度森植被变化信息；

(S4)对30M空间分辨率的TM数据进行云及阴影掩模处理；

(S5)根据粗尺度森林植被变化掩模文件，在30m分辨率TM影像上构建森林植被特征指数；

(S6)粗尺度土地覆盖变化图与时间序列的森林特征指数相结合的森林动态变化信息提取。

具体的方法流程如图1所示，具体的步骤详细描述如下：

(S1)遥感影像几何校正与配准：包括MODIS-NDVI时间序列的影像几何校正；TM影像时间序列的影像几何校正；MODIS-NDVI与TM影像间的相互配准。由于MODIS-NDVI数据是正弦曲线投影需要转换成地理经纬度投影(可利用ENVI软件转换完成)。TM影像间的匹配可根据TM自带的UTM坐标进行匹配，如果匹配误差较大在选择控制点，利用多项式纠正方法纠正。MODIS-NDVI数据与TM数据间直接利用数据自带的地理坐标确定位置。

(S2)利用一年内不同季节的1KM分辨率的MODIS时间序列NDVI数据，获得1KM分辨率的土地覆盖类型图。

在1KM分辨率粗尺度数据上，将地表覆盖类型分为森林、农田、城市建筑、水体四大基本类型。其NDVI值由高到低依次为森林、农田、城市建筑、水体。通过统计分析不同土地覆盖类别各月份的NDVI直方图，确定对应月份各类别NDVI阈值。具体处理过程包括：(1)各地表覆盖类别阈值确定。选择各种覆盖类型的样本，计算各个月份中的均值和方差，其中j是月份，是类别i所选取样本的NDVI平均值，S_i是类别i所选取样本的NDVI方差，各个类别阈值区间被定义为(2)基于决策树的粗尺度土地覆盖分类。由于不同土地覆盖类型的阈值不同，植被的NDVI值通常较高，四大基本土地覆盖类型的NDVI值的高低依次为森林(赋值1)、农田(赋值2)、城市建筑(赋值3)、水体(赋值4)，基本的分类决策思想如图2所示。

(S3)利用多年1KM分辨率的土地覆盖图，提取粗尺度森林植被变化信息；

对多年的1KM分辨率土地覆盖图逐年进行比较，分别计算相邻时间的两幅分类影像差值，像元计算结果为0就代表该像元类别未发生改变，而计算结果非0则表示该像元的土地覆盖类型在监测期间发生了改变。将差值结果中值为非0的像素，且该像素在任意土地覆盖图中为森林类型(像素值为1)的所有像素提取出来，这些像素区域代表了在一段时间范围内粗尺度的森林植被变化信息。利用该信息建立森林植被变化掩模文件，值为1区域表示森林变化，非0表示森林没有变化。

(S4)根据粗尺度森林植被变化掩模文件，选择30M空间分辨率的TM影像，并进行云及阴影掩模处理；

根据粗尺度森林变化掩模文件，提取森林变化发生的地理位置，结合Landsat的Worldwide Reference System 2(WRS-2)计算TM影像的轨道号，并获取该区域位置的时间序列影像。对获取的影像进行云和阴影的掩膜处理。

云一般在反射带表现为高亮，在温度带表现为低温，可以通过这个特点来对云进行掩膜。阴影表现的比较暗，而且云形成的阴影通常在云像元的西北角，其面积和云的面积相近，通过空间分析能确定云形成的阴影。地形造成的阴影在森林指数计算过程中会被筛选并剔除。

云的识别选用最常用的Automated Cloud-Cover Assessment(ACCA)算法，该算法的流程(如图3所示)主要包括四个处理步骤，而每一个步骤中详细处理和参数设定在Characterization of the Landsat-7ETM Automated Cloud-Cover Assessment(ACCA)Algorithm中有详细介绍。

第一步，处理1。通过波段2-波段6这五个波段的影像来识别云光谱，判断暖云、冷云、模糊像元、雪以及实现无云的掩膜。

第二步，波段6的云标记。即使是暖云，其温度也会比一般地物低，将上面已经确定是云的像元标记在波段6中。

第三步，处理2。剩余的模糊像元类别通过波段6再进行一次分类，在模糊像元中能区分一部分冷云和暖云，剩余不能识别的还是留在模糊类别中。

第四步，云像元聚合。将处理1和处理2中识别出来的云聚合在一起就是最终识别的云的结果。

云阴影的识别通过地理位置来判别和确定，由于landsat是上午10点左右过境，因此阴影一般都在云的西北方向。温度随着海拔增高而降低，通过云的温度能反推出云高度的大概位置，然后结合太阳高度角能确定阴影所在的大概位置以寻找云阴影。图4是云阴影识别的详细流程，主要包括三个步骤：

第一步，反推云的高度。温度在高程上以1℃/hm递减，波段6转换成温度以后结合DEM可以反推出云的高度。

第二步，推测云影的位置。结合太阳高度角和第一步中云的高度能反推出云像元对应的该像元对应云影的位置。虽然对于单个云像元而言，这样推出来的结果不准确，但云以聚集的形式出现，云团推测出的结果可信度高。对云团对应云影的位置区域进行5个像元(合理阈值范围)的缓冲。

第三步，确定云影。推测云影的缓冲区像元在波段4和波段5的灰度值都很低的位置被标记为云影。使用b4和b5是因为植被在可见光波段和中红外波段(b7)的灰度值也会很低。由于云的厚度会影像云影的黑暗程度，即b4和b5的阈值需要合理设定或者多次设定。

(S5)根据粗尺度森林植被变化掩模文件，在30m分辨率TM影像上构建森林植被特征指数；

主要步骤包括：1)森林样本自动选择：2)森林特征指数计算。

森林样本自动选择的步骤包括：a)将整个影像区域按照300*300像素，划分为若干小区域；b)在每个小区域内，计算NDVI(Normalized Difference Vegetation Index，归一化植被指数，)。剔除NDVI小于0.2的像素，可剔除非植被暗目标；C)在小区域内统计TM3波段的直方图。统计b)中剩余的像元在红光波段的反射和频数的关系，得到反射-频数分步曲线(y轴为频数，x轴为灰度值)。对曲线进行中值滤波(IDL中median函数实现)得到统计直方图，该图最少有一个峰；d)确定直方图中峰的位置。以频数达到1％的灰度值作为DNred(min)。当直方图中只有一个峰时，峰值对应的灰度值为DNred(max)；当该图有多个峰时，而且后面出现峰离第一个峰很近(两个峰对应的灰度值之差小于10)的时候，也应该被识别为森林峰，DNred(max)是最后一个森林峰对应的灰度值；e)得到森林样本掩膜文件。根据阈值[DNred(min)，DNred(max)]，掩膜红光波段(TM3)得到森林样本(森林样本像元值为1，其它像元值为0)。

森林特征指数计算步骤主要包括：a)根据选择的森林样本，计算波段i(i＝3，5，7)上光谱值的平均值和标准差，分别记为M_i和S_i，b)计算每个波段的特征指数FZ_i，对任意像元在该波段的值为b_i，那么这个像元的森林特征指数(FZ_i)可以用公式(1)计算：

公式(1)中FZ_i的实质是描述了任意像元与森林的偏离程度，该指数越小，对应像元是森林的可能性就越大。

对于具有多个波段的TM影像而言，对所有波段的FZ_i值做如式(2)的积分运算能得到该像元的综合性森林特征指数(IFZ)。

其中NB是所使用的波段总数，本发明中使用的波段是3，5，7波段。

这样计算出来的森林指数能够反映出该像元是森林的可能性，且该指数与森林可能性成负相关，也就是IFZ越高代表该像元是森林的可能性越小，IFZ越低代表该像元是森林的可能性越大。

(S6)粗尺度森林植被变化图与时间序列的森林特征指数相结合的森林动态变化信息提取。

以时间序列的IFZ值与粗尺度的森林植被变化图采用决策树方法提取森林变化信息。具体的详细过程见如图5。

设定IFZ阈值，结合粗检测结果确定，确定每个像素的时间序列IFZ值与阈值的比较结果：

(1)如果时间序列的IFZ值都低于给定的阈值，则该像元在监测期间一直是典型的森林的像元。

(2)如果时间序列IFZ均值大于给定阈值，或者绝大多数的IFZ值大于给定的阈值(＞90％)，则该像元被确定为非森林植被。

(3)如果时间序列的IFZ值小部分(＜15％)高于给定的阈值，且这种高于阈值的时间不能连续超过2年，则在监测期间是疏林地或者非典型森林像元。因为在监测期间疏林地或者非典型森林可能不会一直出现低于阈值的IFZ值，某个年份或者或某几个不连续的年份中高于阈值，这样的像元被确定为疏林地或者非典型森林像元。

(4)如果在前面步骤中没有被分出，且属于粗尺度的植被变化区域的像元，则归于森林变化发生区域。森林的变化有三种情况：a不可逆转的毁林(森林变为非森林，且在监测期内未变成森林)；b造林(表现为非森林变为森林)；c森林干扰(包括采伐更新、自然灾害、病虫害、火灾等)。

其中，a表现为IFZ连续很低(至少连续3年)，突然增高到阈值以上，此后就不再降低或者在后来的某个年份或者或某几个不连续的年份(连续低的时间不超过2年)由于云、阴影等偶然因素低于阈值。

b表现为IFZ连续很高(至少连续3年)，随后慢慢的降低到到阈值以下(这个过程有一个变化趋势)，此后的时间内(时间阈值设定为3年)就不再高于阈值或在后来的某个年份或者或某几个不连续的年份(连续高的时间不超过2年)偶然高于阈值。因为林分的建立是一个渐进的过程，树木虽然在生长，但是要达到能显示出“森林外观”的光谱数据至少需要几年(阈值设定为3年)。云和阴影以及其他因素造成的干扰虽然存在，但同一个像元连续几年(阈值设定为2年)出现偶然干扰的概率非常小。

c类别非常复杂，森林干扰虽然通过IFZ突发的增高易于识别，但森林恢复过程中IFZ不会迅速地从较高的IFZ值降到不受干扰水平的IFZ值，而是表现为逐渐减少。与此同时干扰发生和恢复完成的时间点也不易确定，因此森林动态像元未被分到前两个类别的都属于干扰类。最典型的干扰是监测中期采伐及更新，根据时间序列IFZ可以确定获得森林动态变化信息，包括干扰年度和干扰程度。此时采伐的时间的IFZ表现为忽然涨高点，更新完成时间点IFZ表现为低于阈值，扰动时间段是两个时间点相隔的时间段。而干扰程度可以有干扰发生时间点IFZ大小除以干扰持续的时间长度，或者干扰发生时间点IFZ降低到到阈值整个过程中变化趋势来确定。

以上步骤能实现森林变化由粗到细的检测，其中最复杂的是完成森林变化监测以后的森林干扰特征提取。

本发明提供的多尺度森林动态变化监测方法，充分利用了不同空间分辨率的遥感影像对不同土地覆盖及植被的响应特点，构建了有粗到细的变化检测技术，提高了植被监测效率；针对植被随时间生长变化规律的特点，构建了基于时间序列的植被变化提取方法，提高监测准确性；利用时间序列的数据，不仅可以监测变化，还可以监测提取森林扰动信息。

以上所述，仅是用以说明本发明的具体实施案例而已，并非用以限定本发明的可实施范围，举凡本领域熟练技术人员在未脱离本发明所指示的精神与原理下所完成的一切等效改变或修饰，仍应由本发明权利要求的范围所覆盖。

Claims (7)

1.一种多尺度森林动态变化监测方法，其特征在于所述方法主要包括以下步骤：

步骤一、获取目标对象的MODIS-NDVI时间序列数据和对应的TM影像时间序列数据，并对影像进行几何校正与配准；

步骤二、利用不同的MODIS时间序列NDVI数据，获得土地覆盖类型图；

步骤三、利用土地覆盖类型图，生成粗尺度森植被变化信息；

步骤四、对TM影像进行云及阴影掩模处理；

步骤五、根据粗尺度森林植被变化掩模文件，在TM影像上构建森林植被特征指数，具体为：

根据粗尺度森林植被变化掩模文件，在30m分辨率TM影像上构建森林植被特征指数，包括以下步骤：

步骤5.1森林样本自动选择；

第一步，将整个影像区域按照300*300像素，划分为若干小区域；

第二步，在每个小区域内，计算NDVI，剔除NDVI小于0.2的像素，可剔除非植被暗目标；

第三步，在小区域内统计掩膜红光波段波段的直方图，统计第二步中剩余的像元在红光波段的反射和频数的关系，得到反射-频数分步曲线，对曲线进行中值滤波得到统计直方图；

第四步，确定直方图中峰的位置，以频数达到1％的灰度值作为DNred(min)，当直方图中只有一个峰时，峰值对应的灰度值为DNred(max)；当该图有多个峰时，而且后面出现峰离第一个峰很近的时候，也应该被识别为森林峰，DNred(max)是最后一个森林峰对应的灰度值；

第五步，得到森林样本掩膜文件，根据阈值[DNred(min)，DNred(max)]，掩膜红光波段得到森林样本，其中森林样本像元值为1，其它像元值为0；

步骤5.2森林特征指数计算；

第一步，根据选择的森林样本，计算波段i上光谱值的平均值和标准差，分别记为M_i和S_i，其中i＝3，5，7；

第二步，计算每个波段的特征指数FZ_i，对任意像元在该波段的值为b_i，那么这个像元的森林特征指数FZ_i可以用公式(1)计算：

<mrow> <msub> <mi>FZ</mi> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>b</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>M</mi> <mi>i</mi> </msub> <mo>|</mo> </mrow> <msub> <mi>S</mi> <mi>i</mi> </msub> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

其中，公式(1)中FZ_i的实质是描述了任意像元与森林的偏离程度，该指数越小，对应像元是森林的可能性就越大；

对于具有多个波段的TM影像而言，对所有波段的FZ_i值做如公式(2)的积分运算能得到该像元的综合性森林特征指数IFZ；

<mrow> <mi>I</mi> <mi>F</mi> <mi>Z</mi> <mo>=</mo> <msqrt> <mrow> <mfrac> <mn>1</mn> <mrow> <mi>N</mi> <mi>B</mi> </mrow> </mfrac> <munderover> <mi>&amp;Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mi>B</mi> </mrow> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>FZ</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

其中NB是所使用的波段总数；

步骤六、粗尺度土地覆盖变化图与时间序列的森林特征指数相结合的森林动态变化信息提取。

2.根据权利要求1所述的多尺度森林动态变化监测方法，其特征在于：所述步骤一的遥感影像几何校正与配准包括：

步骤1.1MODIS-NDVI时间序列的影像几何校正；

步骤1.2TM影像时间序列的影像几何校正；

步骤1.3MODIS-NDVI与TM影像间的相互配准；

所述TM影像间的匹配根据TM自带的UTM坐标进行匹配，如果匹配误差较大在选择控制点，利用多项式纠正方法纠正；

所述MODIS-NDVI数据与TM数据间直接利用数据自带的地理坐标确定位置。

3.根据权利要求1所述的多尺度森林动态变化监测方法，其特征在于：所述步骤二具体为利用一年内不同季节的1KM分辨率的MODIS时间序列NDVI数据，获得1KM分辨率的土地覆盖类型图，在1KM分辨率粗尺度数据上，将地表覆盖类型分为森林、农田、城市建筑、水体四大基本类型，其NDVI值由高到低依次为森林、农田、城市建筑、水体，通过统计分析不同土地覆盖类别各月份的NDVI直方图，确定对应月份各类别NDVI阈值，具体处理过程包括以下步骤：

步骤2.1各地表覆盖类别阈值确定；选择各种覆盖类型的样本，计算各个月份中的均值和方差，各个类别阈值区间被定义为其中j是月份，是类别i所选取样本的NDVI平均值，S_i是类别i所选取样本的NDVI方差；

步骤2.2基于决策树的粗尺度土地覆盖分类；四大基本土地覆盖类型的NDVI值的高低依次为森林、农田、城市建筑、水体。

4.根据权利要求1所述的多尺度森林动态变化监测方法，其特征在于：所述步骤三具体为利用多年1KM分辨率的土地覆盖图，提取粗尺度森林植被变化信息，包括以下步骤：

步骤3.1对多年的1KM分辨率土地覆盖图逐年进行比较，分别计算相邻时间的两幅分类影像差值，像元计算结果为0就代表该像元类别未发生改变，而计算结果非0则表示该像元的土地覆盖类型在监测期间发生了改变；

步骤3.2将差值结果中值为非0的像素，且该像素在任意土地覆盖图中为森林类型的所有像素提取出来；

步骤3.3利用该信息建立森林植被变化掩模文件，值为1区域表示森林变化，非0表示森林没有变化。

5.根据权利要求1所述的多尺度森林动态变化监测方法，其特征在于：所述步骤四具体为根据粗尺度森林植被变化掩模文件，选择30M空间分辨率的TM影像，并进行云及阴影掩模处理，包括以下步骤：

步骤4.1云的识别；

第一步，处理1，通过波段2-波段6这五个波段的影像来识别云光谱，判断暖云、冷云、模糊像元、雪以及实现无云的掩膜；

第二步，波段6的云标记，将第一步已经确定是云的像元标记在波段6中；

第三步，处理2，剩余的模糊像元类别通过波段6再进行一次分类，在模糊像元中能区分一部分冷云和暖云，剩余不能识别的还是留在模糊类别中；

第四步，云像元聚合，将处理1和处理2中识别出来的云聚合在一起就是得到的云的结果；

步骤4.2云阴影的识别；通过地理位置来判别和确定来进行云阴影的识别，

第一步，反推云的高度，温度在高程上以1℃/hm递减，波段6转换成温度以后结合DEM反推出云的高度；

第二步，推测云影的位置，结合太阳高度角和第一步中云的高度反推出云像元对应的该像元对应云影的位置，对云团对应云影的位置区域进行5个像元的缓冲；

第三步，确定云影，推测云影的缓冲区像元在波段4和波段5的灰度值都很低的位置被标记为云影。

6.根据权利要求1所述的多尺度森林动态变化监测方法，其特征在于：所述步骤六具体为粗尺度森林植被变化图与时间序列的森林特征指数相结合的森林动态变化信息提取，以时间序列的IFZ值与粗尺度的森林植被变化图采用决策树方法提取森林变化信息，通过以下方法进行：

设定IFZ阈值，结合粗检测结果确定，确定每个像素的时间序列IFZ值与阈值的比较结果：

如果时间序列的IFZ值都低于给定的阈值，则该像素的像元在监测期间一直是典型的森林的像元；

如果时间序列IFZ均值大于给定阈值，或者绝大多数的IFZ值大于给定的阈值，则该像元被确定为非森林植被；

如果时间序列的IFZ值小部分高于给定的阈值，且这种高于阈值的时间不能连续超过2年，则在监测期间是疏林地或者非典型森林像元；

如果在前面步骤中没有被分出，且属于粗尺度的植被变化区域的像元，则归于森林变化发生区域。

7.根据权利要求6所述的多尺度森林动态变化监测方法，其特征在于：所述森林变化有三种情况：

a类别为不可逆转的毁林，即森林变为非森林，且在监测期内未变成森林；

b类别为造林，即表现为非森林变为森林；

c类别为森林干扰，包括采伐更新、自然灾害、病虫害和火灾；

其中，a类别表现为IFZ连续很低，突然增高到阈值以上，此后就不再降低或者在后来的某个年份或者或某几个不连续的年份由于云和阴影的偶然因素低于阈值；

b类别表现为IFZ连续很高，随后慢慢的降低到到阈值以下，此后的时间内就不再高于阈值或在后来的某个年份或者或某几个不连续的年份偶然高于阈值；

c类别表现，森林干扰虽然通过IFZ突发的增高易于识别，但森林恢复过程中IFZ不会迅速地从较高的IFZ值降到不受干扰水平的IFZ值，而是表现为逐渐减少。