CN105046188A - Modis混合像元分解森林信息提取方法 - Google Patents

Abstract

一种MODIS混合像元分解森林信息提取方法，综合MOD09A1地表反射率产品和MOD13Q1植被指数产品标准影像，其特征是MODIS不同产品影像具有不同的森林信息，但由于其分辨率低，影像中存在大量混合像元，通过改进端元提纯方法，利用线性混合像元分解模型提取森林信息。本发明的有益效果是，可以在大范围快速提取森林覆盖信息的同时，提高森林制图精度以及森林类型识别精度。信息提取中仅采用MOD09A1地表反射率产品和MOD13Q1植被指数产品，效果理想，简单易行。

Description

MODIS混合像元分解森林信息提取方法

所属技术领域

本发明涉及一种MODIS遥感影像森林信息提取的技术，尤其是能利用混合像元分解提高森林覆盖制图与森林类型识别精度遥感信息提取方法。

背景技术

目前中高分辨率的遥感影像，虽然测量精度比较高，但数据的价格比较贵、获取不容易；低分辨率遥感影像覆盖范围比较广，且费用比较低廉，如MODIS影像，由于其最高空间分辨率为250m，且MODIS遥感影像具有很高的光谱辐***度，能够免费获取，但存在大量混合像元，在地物的分类时，常常会造成一定的误差，利用传统方法难以获得理想的分类精度。

发明内容

为了克服现有的MODIS影像森林信息提取方法分类精度不能满足生产需要和森林类型识别能力较差的不足,本发明提供一种森林信息提取方法，该方法既能有效利用MODIS数据光谱辐***度极佳的特点，又可改善MODIS空间分辨率的不足，且能方便快捷地利用MODIS影像提供的MOD09A1地表反射率产品和MOD13Q1植被指数产品有效提高森林覆盖制图和森林类型识别精度。

为了实现上述技术目的，本发明的技术方案是，一种MODIS混合像元分解森林信息提取方法，包括如下步骤：

步骤一：数据获取步骤，获取森林遥感影像；

步骤二：预处理步骤，对步骤一得到的遥感影像进行预处理，获取预处理后的遥感影像数据；

步骤三：森林类型物候差异特征分析步骤，找出能够明显区分各森林类型的植被指数时期；

步骤四：端元组分提纯步骤，获取终端像元组分；

步骤五：端元组分光谱反射率求取步骤，获取各地物光谱反射率；

步骤六：线性混合像元分解步骤，求取各地物分布图；

步骤七：森林信息提取步骤，获取各地物分布图。

所述的方法，所述的步骤二中遥感数据的预处理包括：对遥感数据进行投影转换，大气校正，图像拼接，消除黑色条带，几何校正，研究区提取。

所述的方法，所述的步骤三中包括的步骤为：

利用GPS获取典型样点的坐标信息，利用GIS软件获取典型样点在遥感影像上的灰度值，构建植被指数时间序列剖面曲线，然后通过物候分析找出的能够明显区分各森林类型的植被指数遥感数据与MOD09A1地表反射率产品相结合所得到的数据。

所述的方法，所述的典型样点包括针叶林、阔叶林、竹林、灌木林、水域、耕地和建设用地。

所述的方法，所述的通过物候分析找出的能够明显区分各森林类型的植被指数遥感数据与MOD09A1地表反射率产品相结合所得到的数据包括以下步骤：

首先构建端元提纯模型，利用决策树分类计算机图形学分类方法，根据各森林类型物候差异的特征，利用找出的能够明显区分各森林类型的植被指数时期遥感数据，并且利用能够区分各地物的植被指数灰度值作为决策树分类模型的阈值：即

当第8期NDVI的阈值小于0.6则为耕地、水域、建设用地，否则为植被覆盖高的森林；第9期的NDVI值小于0.25时，则为水域，否则为耕地和建设用地；当第2期的NDVI小于0.45，并且第9期的EVI小于0.32时，则为建设用地，否则为耕地；第16期NDVI小于0.45时，则为灌木林，否则为针叶林、阔叶林、竹林；当第18期的NDVI小于0.75，并且第16期的EVI大于0.35时，则为针叶林，否则为阔叶林和竹林；当第11期的EVI小于0.65，并且第15期的NDVI小于0.75时，则为竹林，反之则为阔叶林；最后选择NDVI的第2、8、9、15、16、18期以及EVI的第9、11、16期数据进行合成。

所述的方法，所述的步骤四中包括的步骤为：

对端元组分进行初步提纯，对步骤三中得到的数据进行最小噪声分离变换、纯净指数分析，利用N维散点图进行N维散度分析，初步确定终端像元，将在N维散点图中得到各自的初始终端像元，作为感兴趣区代入到决策树分类后的图像，利用决策树分类建立端元提纯模型对初始终端像元进一步提纯，对与决策树分类得到的分布图不一致的终端像元的散点进行剔除或修改，得到最终的终端像元。

所述的方法，所述的步骤五中包括的步骤为：

利用最终的终端像元求取各地物在各波段的光谱反射率，并结合GPS定位样地点，通过遥感影像获取的各波段光谱反射率，确定各地物在各波段的最终地物反射率值。

所述的方法，所述的步骤六中包括的步骤为：

利用各地物在各波段的最终地物反射率值，通过线性混合像元分解模型

$R=\underset{k=1}{\overset{m}{\Σ}}{F}_k{r}_k+{\mathrm{\ϵ}}_k,$

其中m为端元数目，r_k是第k个端元组的反射率，F_k像元中第k端元组分在像元中所占的面积比，ε_k为第k波段的误差，其中带约束性的线性混合像元的分解要满足各像元的组分比F_k之和为1,对遥感影像进行混合像元分解，得到各地物的丰度图。

所述的方法，所述的步骤七中包括的步骤为：

在每组组分图像中选择终端单元丰度值大于0.5的像元，再把每组组分图像组成一幅图像，利用最大似然法分类对该图像进行分类，得到各森林类型分布图，实现混合像元分解森林信息提取。

本发明解决其技术问题所采用的技术方案是：在对采用的MODIS标准产品影像进行投影转换、图像拼接、裁剪、去黑色条带等预处理基础上，首先，利用森林类型物候差异特征，采用MOD13Q1植被指数产品，通过对归一化植被指数NDVI以及增强型植被指数EVI进行处理，构建NDVI、EVI时间序列剖面曲线，通过决策树分类计算机图形学分类方法，构建端元提纯模型；其次，采用构建端元提纯模型的MOD13Q1植被指数产品与MOD09A1地表反射率产品相结合组成新的数据，并利用构建的端元提纯模型进行端元进一步提纯；再次，采用线性混合像元分解模型对组合成的新数据进行分解；最后，采用MOD09A1地表反射率产品和MOD13Q1植被指数产品标准影像产品，通过最大似然法分类，实现森林信息提取。本发明的有益效果是，可以在大范围快速提取森林覆盖信息的同时，提高森林制图精度以及森林类型识别精度。信息提取中仅采用MOD09A1地表反射率产品和MOD13Q1植被指数产品，效果理想，简单易行。

附图说明

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

图1是本发明的技术原理图；

图2是本发明实施例的森林信息提取结果图；

图3是本发明的端元提纯决策树分类模型图。

具体实施方式

本发明包括以下步骤：

步骤一：数据获取步骤，获取森林遥感影像；

步骤二：预处理步骤，对步骤一得到的遥感影像进行预处理，获取预处理后的遥感影像数据，预处理包括：对遥感数据进行投影转换，大气校正，图像拼接，消除黑色条带，几何校正，研究区提取；

步骤三：森林类型物候差异特征分析步骤，找出能够明显区分各森林类型的植被指数时期；具体来说，首先利用GPS获取(针叶林、阔叶林、竹林、灌木林、水域、耕地、建设用地)的坐标信息，利用GIS软件获取针叶林、阔叶林、竹林、灌木林、水域、耕地、建设用地野外典型样点在遥感影像上的灰度值，构建植被指数时间序列剖面曲线，通过物候分析找出的能够明显区分各森林类型的植被指数遥感数据与MOD09A1地表反射率产品相结合组成新的数据。当第8期NDVI的阈值小于0.6则为耕地、水域、建设用地，否则为植被覆盖高的森林；第9期的NDVI值小于0.25时，则为水域，否则为耕地和建设用地；当第2期的NDVI小于0.45，并且第9期的EVI小于0.32时，则为建设用地，否则为耕地；第16期NDVI小于0.45时，则为灌木林，否则为针叶林、阔叶林、竹林；当第18期的NDVI小于0.75，并且第16期的EVI大于0.35时，则为针叶林，否则为阔叶林和竹林；当第11期的EVI小于0.65，并且第15期的NDVI小于0.75时，则为竹林，反之则为阔叶林；因此选择了NDVI的第2、8、9、15、16、18期以及EVI的第9、11、16期数据进行合成。

步骤四：端元组分提纯步骤，获取终端像元组分；

首先构建端元提纯模型，利用决策树分类计算机图形学分类方法，根据各森林类型物候差异的特征，利用找出的能够明显区分各森林类型的植被指数时期遥感数据，并且利用能够区分各地物的植被指数灰度值作为决策树分类模型的阈值，得到端元提纯决策树分类模型如图3，模型中，“ZYL”表示针叶林，“KYL”表示阔叶林，“ZL”表示竹林，“GML”表示灌木林，“SY”表示水域，“GD”表示耕地，“JSYD”表示建设用地，B2、B8、B9、B15、B16、B18分别代表NDVI的第2、8、9、15、16、18期，b9、b11、b16、分别代表EVI的第9、11、16期。

然后对端元组分进行初步提纯，对步骤三中得到的数据进行最小噪声分离变换、纯净指数分析，利用N维散点图进行N维散度分析，初步确定终端像元，将在N维散点图中得到各自的初始终端像元，作为感兴趣区代入到决策树分类后的图像，利用决策树分类建立端元提纯模型对初始终端像元进一步提纯，对一些不符合规律及错误(即与决策树分类得到的分布图不一致)的终端像元的散点进行剔除或修改，得到最终的终端像元。

步骤五：端元组分光谱反射率求取步骤，获取各地物光谱反射率；即利用最终的终端像元求取各地物在各波段的光谱反射率，并结合GPS定位样地点，通过遥感影像获取的各波段光谱反射率，确定各地物在各波段的最终地物反射率值。

步骤六：线性混合像元分解步骤，求取各地物分布图；即利用各地物在各波段的最终地物反射率值，通过线性混合像元分解模型其中m为端元数目，r_k是第k个端元组的反射率，F_k像元中第k端元组分在像元中所占的面积比，ε_k为第k波段的误差。其中带约束性的线性混合像元的分解要满足各像元的组分比F_k之和为1,对遥感影像进行混合像元分解，得到各地物的丰度图。

步骤七：森林信息提取步骤，获取各地物分布图。即在每组组分图像中选择终端单元丰度值大于0.5的像元，再把每组组分图像组成一幅图像，利用最大似然法分类对该图像进行分类，得到各森林类型分布图，实现混合像元分解森林信息提取。

参见图1，在实施时，需要对MOD09A1地表反射率产品和MOD13Q1植被指数产品进行投影转换、图像拼接、裁剪、去黑色条带等预处理，通过对MOD13Q1归一化植被指数NDVI以及增强植被指数EVI，构建NDVI、EVI时间序列剖面曲线，分析得出构建决策树提纯模型的阈值。并把构建决策树模型的阈值数据产品与MOD09A1地表反射率产品相结合生产新的数据，利用线性混合像元分解模型进行分解，再利用最大似然法进行分类提取森林信息。

在图2所示实施例中，选取湖南省为对象，采用2009年的MOD09A1地表反射率产品和MOD13Q1植被指数产品进行投影转换、图像拼接、裁剪、去黑色条带等预处理。首先，利用森林类型物候差异特征，采用MOD13Q1植被指数产品，通过对归一化植被指数NDVI以及增强型植被指数EVI(全年各23期)进行处理，构建NDVI、EVI时间序列剖面曲线，筛选出最优决策树分类模型的数据，构建端元提纯模型；其次，采用构建端元提纯模型的MOD13Q1植被指数产品与MOD09A1地表反射率产品相结合组成新的数据，并进行最小噪音分离(MNF)，像元纯度指数计算(PPI),密度分割，N维散点图，端元组分初步提纯处理，并利用决策树提纯模型进行进一步提纯，结合影像建标各波段的光谱反射率，得出最终端元组分反射率；再次，采用线性混合像元分解模型对组合成的数据进行分解；最后，采用MOD09A1地表反射率产品和筛选出来的MOD13Q1植被指数产品组合成的数据，通过最大似然法分类，实现森林信息提取，处理输出湖南省森林分布图。

本发明提出的MODIS混合像元分解森林信息提取方法大幅度提高了森林制图与森林类型识别精度，分类结果理想，具有重要的实用推广意义。

Claims (9)

1.一种MODIS混合像元分解森林信息提取方法，其特征在于，包括如下步骤：

步骤一：数据获取步骤，获取森林遥感影像；

步骤二：预处理步骤，对步骤一得到的遥感影像进行预处理，获取预处理后的遥感影像数据；

步骤三：森林类型物候差异特征分析步骤，找出能够明显区分各森林类型的植被指数时期；

步骤四：端元组分提纯步骤，获取终端像元组分；

步骤五：端元组分光谱反射率求取步骤，获取各地物光谱反射率；

步骤六：线性混合像元分解步骤，求取各地物分布图；

步骤七：森林信息提取步骤，获取各地物分布图。

2.根据权利要求1所述的方法，其特征在于，所述的步骤二中遥感数据的预处理包括：对遥感数据进行投影转换，大气校正，图像拼接，消除黑色条带，几何校正，研究区提取。

3.根据权利要求2所述的方法，其特征在于，所述的步骤三中包括的步骤为：

利用GPS获取典型样点的坐标信息，利用GIS软件获取典型样点在遥感影像上的灰度值，构建植被指数时间序列剖面曲线，然后通过物候分析找出的能够明显区分各森林类型的植被指数遥感数据与MOD09A1地表反射率产品相结合所得到的数据。

4.根据权利要求3所述的方法，其特征在于，所述的典型样点包括针叶林、阔叶林、竹林、灌木林、水域、耕地和建设用地。

5.根据权利要求3所述的方法，其特征在于，所述的通过物候分析找出的能够明显区分各森林类型的植被指数遥感数据与MOD09A1地表反射率产品相结合所得到的数据包括以下步骤：

首先构建端元提纯模型，利用决策树分类计算机图形学分类方法，根据各森林类型物候差异的特征，利用找出的能够明显区分各森林类型的植被指数时期遥感数据，并且利用能够区分各地物的植被指数灰度值作为决策树分类模型的阈值：即

当第8期NDVI的阈值小于0.6则为耕地、水域、建设用地，否则为植被覆盖高的森林；第9期的NDVI值小于0.25时，则为水域，否则为耕地和建设用地；当第2期的NDVI小于0.45，并且第9期的EVI小于0.32时，则为建设用地，否则为耕地；第16期NDVI小于0.45时，则为灌木林，否则为针叶林、阔叶林、竹林；当第18期的NDVI小于0.75，并且第16期的EVI大于0.35时，则为针叶林，否则为阔叶林和竹林；当第11期的EVI小于0.65，并且第15期的NDVI小于0.75时，则为竹林，反之则为阔叶林；最后选择NDVI的第2、8、9、15、16、18期以及EVI的第9、11、16期数据进行合成。

6.根据权利要求3所述的方法，其特征在于，所述的步骤四中包括的步骤为：

对端元组分进行初步提纯，对步骤三中得到的数据进行最小噪声分离变换、纯净指数分析，利用N维散点图进行N维散度分析，初步确定终端像元，将在N维散点图中得到各自的初始终端像元，作为感兴趣区代入到决策树分类后的图像，利用决策树分类建立端元提纯模型对初始终端像元进一步提纯，对与决策树分类得到的分布图不一致的终端像元的散点进行剔除或修改，得到最终的终端像元。

7.根据权利要求6所述的方法，其特征在于，所述的步骤五中包括的步骤为：

利用最终的终端像元求取各地物在各波段的光谱反射率，并结合GPS定位样地点，通过遥感影像获取的各波段光谱反射率，确定各地物在各波段的最终地物反射率值。

8.根据权利要求7所述的方法，其特征在于，所述的步骤六中包括的步骤为：

利用各地物在各波段的最终地物反射率值，通过线性混合像元分解模型

$R=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{F}_k{r}_k+{\mathrm{\ϵ}}_k,$

其中m为端元数目，r_k是第k个端元组的反射率，F_k像元中第k端元组分在像元中所占的面积比，ε_k为第k波段的误差，其中带约束性的线性混合像元的分解要满足各像元的组分比F_k之和为1,对遥感影像进行混合像元分解，得到各地物的丰度图。

9.根据权利要求8所述的方法，其特征在于，所述的步骤七中包括的步骤为：

在每组组分图像中选择终端单元丰度值大于0.5的像元，再把每组组分图像组成一幅图像，利用最大似然法分类对该图像进行分类，得到各森林类型分布图，实现混合像元分解森林信息提取。