基于分层乘法模型和高分1号数据的森林覆盖度提取方法

1.本发明属于森林遥感技术领域，涉及一种基于分层乘法模型和高分1号数据的森林覆盖度提取方法。

背景技术：

2.生态环境保护是全球的共识，也是当今社会的热点问题。作为生态建设主体的林业在维护环境健康和生态平衡、改善全球碳循环和水质、提供休闲潜力等方面发挥着至关重要的作用。然而，经济的快速增长、人口的急速膨胀和自然灾害的频发等因素对林业的健康生长构成了严重威胁。森林覆盖率是反映森林资源的丰富程度和生态平衡状况的重要指标，因此，基于高分辨率遥感影像进行快速、高效的森林覆盖率调查已成为生态环境建设与保护的重要研究课题。
3.传统的森林覆盖率及其变化监测主要依赖国外遥感数据，例如landsat，modis，spot，quickbird等。但随着我国高分辨率对地观测系统重大专项(简称高分专项)的实施，基于国产高分卫星的森林研究成为林业遥感的重要内容。国产高分数据已经被广泛的用于估算森林落面积、调查的森林生产变化、估算森林生物量和碳储量等多个领域。
4.基于遥感影像提取森林覆盖度的常用方法有单波段反射率计算方法、基于代数计算的光谱指数、监督分类、非监督分类、人工神经网络、专家系统等。其中单个光谱带或带的总和对植被特性的变化很敏感，经常用于时间序列图像的创建。归一化差异湿度指数(ndmi)和归一化燃烧比(nbr)指数经常用于检测森林采伐和灾害。土壤调整植被指数(savi)和抗大气植被指数(arvi)用于分析土壤效应和大气噪声。监督和非监督分类方法精度较低，人工神经网络和专家系统较为复杂，可操作性差。这些方法无法消除混合像元的影响。
5.光谱混合分析(sma)是一种成熟且有效的解决波谱混合问题的技术。sma将混合光谱建模为光谱端元的线性或非线性组合从而获取亚像素植被信息，并且被广泛用于林业生态调查和监测。sma是从混合像素中提取纯净光谱成分的方法，但在运用多光谱影像如gaofen-1 wfv时存在两个挑战：1)端元的数量受限于遥感数据的波段数。根据混合像素分解理论，最大端元数等于光谱波段数加一，因此无法从只有四个波段的gf-1wfv数据中准确提取更多的地物类型；2)很难直接评估混合像素的解混合精度。由于通过地面测量难以确定图像中每个像素的真实组成，直接解混合精度评估成为一个难题。
6.因此如何提供一种利用有限波段的多光谱遥感数据快速、准确的提取森林覆盖度的方法是当前林业遥感亟待解决的问题。
7.为了解决这些问题，提出了一种基于像元分解的分层乘法模型(hierarchical multiplication model,hmm)。本研究的内容包括：(1)利用分层乘法模型(hmm)和gf-1wfv数据提取某区域森林覆盖度；2)基于高分辨率(2m)的gf-1pms数据和svm方法对试验区域土地覆盖类型进行精细分类，根据分类结果计算森林覆盖度，并以此覆盖度作为验证数据对hmm结果的精度进行分析。

技术实现要素：

8.本发明的目的就是提供一种基于混合像元分解的分层乘法模型(hmm)和高分1号数据(gf1数据)的森林覆盖度提取方法，利用分层乘法模型从多光谱遥感数据中准确、高效的提取森林覆盖度。在对遥感数据预处理的基础上，首先按照地物波谱特征差异将地物类型划分为三个主要层级。其次在每一层级中利用线性像元分解模型进行地物丰度计算，最后通过乘法模型计算林业植被在第一级地物中的丰度。
9.本发明具体包括以下步骤：
10.步骤一、影像预处理：
11.获取目标对象的gf1数据，借助gf1传感器的辐射定标参数对原始影像进行辐射定标、大气校正和几何校正。辐射定标将原始dn值转化为像元辐射亮度值，大气校正是利用flaash模型(fast line-of-sight atmospheric analysis of spectral hypercubes)将辐射亮度值转化为地表实际反射率，几何校正首先选取同名地物点，采用最邻近像元法进行重采样，校正误差控制在0.2个像元。
12.步骤二、土地覆盖类型层级划分：根据地物波谱差异的大小将地物覆盖类型分类三个层级，第一层级为水体-非水体层，第二层级为植被-非植被层，第三层级为森林和非森林层；
13.步骤三、各层级内进行像元分解，提取各层级对应元素的丰度：
14.像元分解是从多地物光谱混合的数据中提取各种地物(端元)以及各成分的比例(丰度)。线性光谱混合模型(linear spectral mixture analysis，lsma)定义像元在某一光谱波段的反射率是由构成像元的基本组分的反射率以其所占像元面积比例为权重系数的线性组合。lsma通过下式计算：
[0015][0016]
式中r
iλ
是λ波段上的i像元的波谱反射；n是地物数量；f
ki
是像元i内端元k所占的面积比例；r
kλ
是λ波段上端元k的波谱反射率；ε
iλ
是λ波段上像元i的误差。
[0017]
比例系数fk应满足如下两个约束条件：
[0018][0019]fki
≥1
ꢀꢀ
(3)；
[0020]
比例系数f
ki
通过最小二乘法计算，计算公式如下：
[0021][0022]
式中ε
iλ
,r
iλ
,r
iλ
分别是λ波段上像元i的误差值，观察波谱值和估算波谱值。
[0023]
分别按照步骤二所述三个层次选取端元，并利用lsma模型在各层级进行像元分解提取对应地物的丰度。
[0024]
步骤3.1、在第一层内提取水体和非水体：
[0025]
水体端元选取自河流、湖泊、水库等清澈水体，非水体主要取自建筑物、道路、植被、裸地等，再利用lsma模型和最小二乘法计算出水体和非水体的丰度。虽然纯净水体和含水量低的非水体的波谱差异明显，但高含水量的物质受水份影响较大，为了减少这种干扰，归一化水指数(ndwi)被计算以进一步修正水体和非水体的提取误差。ndwi计算公式为(绿波-近红外)/(绿波+近红外)，只用ndwi指数小于零的非水体像元才能在下一级中被进一步
分解。
[0026]
步骤3.2在第二层提取出植被和非植被:
[0027]
由于光线的影响，光照下的植被和阴影下的植被光谱有较大的差异，因此本层级的植被被进一步分为光照植被和阴影植被，植被端元分别由光照条件下的植被和阴影下的植被选取。非植被被进一步分为裸地、高光不透水层和普通不透水层三类，非植被端元分别从对应的裸地、高光建筑和普通不透水层选取。利用lsma模型分别计算出光照植被、阴影植被、裸地、高光不透水层和一般不透水层的丰度，进而叠加光照植被和阴影植被为植被，叠加其余子类为非植被。
[0028]
步骤3.3、在第三层提取出森林和非森林:
[0029]
植被主要包括森林、草地和农田，草地表面均匀，反射率高，森林表面粗糙反射率低，农田的反射率介于二者之间。相对于草地和农田，森林受光照影响明显，因此本层森林首先被分解为光照森林和阴影森林，与农田、草地一起利用lsma模型提取出各自丰度，然后光照森林和阴影森林叠加为森林类型，农田和草地叠加为非森林类。
[0030]
步骤四、利用丰度值和分层乘法模型(hmm)计算森林的覆盖度：
[0031]
分层乘法模型(hmm)是一种分步计算丰度问题的方法，在计算时有明显的先后顺序,先完成上一步再进行下一步，即多次计算各层级丰度求综合丰度。假设第一级包含m个元素，计算每种元素的丰度为f
11
,f
12
,
……f1m
,其中每种元素又包含若干次级元素；若第一级第i个元素包含n个次级元素，其对应的丰度为f
21
,f
22
,
…f2n
,则第二级第j个元素在第一级的丰度，即综合丰度计算为f＝f
1i
*f
2j
.。以此类推，任一个元素的综合丰度可以通过这个元素在本级的丰度乘以上一级的综合丰度计算获得。例如第k级第t个元素的综合丰度计算公式为：
[0032]
f＝f
1i
*f
2j
*
…fkt
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)；
[0033]
在森林覆盖度计算过程中，利用像元分解模型计算出来的森林覆盖度是第三层级即森林-非森林层级的森林丰度，其综合丰度是所属植被类型在第二层级丰度的函数，而植被的综合丰度又是第一级非水体丰度的函数。根据lsma模型计算结果，在第三层级森林丰度为f3，森林所属的植被类型在第二级的丰度为f2，植被所属的非水体类型在第一级的丰度为f1，根据分层乘法模型，森林在第一级的丰度即森林的综合覆盖度为f＝f1*f2*f3。
[0034]
与现有技术相比，本发明具有如下优点：
[0035]
(1)本发明创建了基于乘法模型的像元分解算法，提高了基于多光谱影像提取地物丰度值的适用性和准确性。
[0036]
(2)本发明解决了多光谱波段不足的限制，扩大了多光谱影像和像元分解算法在地物分类、识别中的应用。
[0037]
(3)本发明针对像元分解精度难以验证的问题，提出了利用高精度影像进行分类、并以分类结果计算大尺度上各端元的丰度以验证像元分解精确度，提高了像元分解验证的准确性。
附图说明
[0038]
图1为实施例中三个层级的波谱特征提取；
[0039]
图2为实施例中基于hmm方法获得的森林丰度；
[0040]
图3为实施例中基于svm和高分辨率影像的地物分类结果；
[0041]
图4为实施例中基于svm分类结果的森林覆盖度；
[0042]
图5为实施例中hmm方法的误差结果分布；
[0043]
图6为实施例中误差分布分析。
具体实施方式
[0044]
为了使本发明的上述目的、特征和优点能够更加明确易懂，下面结合具体案例对本发明的具体实施方式做详细说明。
[0045]
本实施例中研究区为浙江省陆地区域，浙江面积10.18万平方公里，70.4％的山地，23.2％的平原和盆地，6.4％的水体。大部分山地位于浙江中南部区域，平均海拔高度800m.平原主要包括北部的冲积平原，东北的海岸和中部的盆地。浙江整体呈现西南高、东北低。浙江具有明显的亚热带季风气候，四季分明、阳光充沛，年平均温度15到19℃，平均降水量1100至1900mm,优越的地理条件造就了浙江的“鱼米之乡”称号。
[0046]
基于分层乘法模型和高分1号数据的提取浙江省森林覆盖度的具体步骤如下：
[0047]
步骤一、从中国资源卫星应用中心下载国产高分1号宽幅遥感影像gf-1wfv，影像空间分辨率位16m,包含4个波段。利用gf1传感器的辐射定标参数对原始影像进行辐射定标、将原始dn值转化为像元辐射亮度值；大气校正是利用flaash模型将辐射亮度值转化为地表实际反射率，实现大气校正；选取同名地物点，采用最邻近像元法进行重采样进行几何校正，校正误差控制在0.2个像元。
[0048]
步骤二、根据波谱差异大小进行土地覆盖类型分层，三个层级分别是水体-非水体层，植被-非植被层，以及森林和非森林层；
[0049]
步骤2.1、从影像中选择河流、湖泊、人工池塘、森林、农田、草地、裸地、岩石、城市建筑、道路等各种类型，提取每种类型的光谱特征，通过聚类分析，确定第一级土地类型分类为水体-非水体。
[0050]
步骤2.2、选取一级分类中非水体的森林、农田、草地、裸地、岩石、城市建筑、道路等类型，提取每种类型的光谱特征，通过聚类方法分析确定第二级土地类型分类为植被、非植被。
[0051]
步骤2.3、选取第二级分类中属于植被的森林、草地、农田等，把它们分成第三级类型，即森林和非森林。
[0052]
步骤三、在各层级内进行像元分解，提取各层级对应元素的丰度。
[0053]
步骤3.1、在第一层级内分解出水体和非水体。从遥感影像选取水体和非水体的端元信息。水体端元选取自河流、湖泊、水库等清澈水体，非水体主要取自建筑物、道路、植被、裸地等,如图1a所示，水体和非水体端元的波谱特征曲线差别显著。
[0054]
利用线性像元分解方法提取水体和非水体的丰度。考虑浑浊水体和高含水量的物质有较高的相似度，为了减少这种干扰，归一化水指数(ndwi)被计算以掩膜水体。ndwi计算公式为(绿波-近红外)/(绿波+近红外)，只用ndwi指数小于零的非水体像元才能在下一级中被进一步分解。
[0055]
步骤3.2、在第二层级分解出植被和非植被。植被冠层表面粗糙，由于光线的影响，光照下的植被和阴影下的植被光谱有较大的差异，因此本层级的植被被进一步细分为光照
植被和阴影植被，其端元分别选自光照条件下的植被和阴影下的植被。非植被包括裸地、高光效的不透水层和不透水层三类，它们的端元分别从对应的裸地、高光建筑和普通不透水层选取，其端元的特征波谱曲线如图1b。
[0056]
基于以上五类端元，利用线性像元分解模型分别获得光照植被、阴影植被、裸地、高光不透水层和一般不透水层的丰度，最后合并光照植被和阴影植被为本层级植被丰度，合并其余子类丰度为非植被类型丰度。
[0057]
步骤3.3、在第三层级分解出森林和非森林。考虑光照条件对森林光谱的影响，本层森林首先被分为光照森林和阴影森林两个次级类，光照森林端元从冠层顶部和光照条件较好的区域选取，阴影森林从背光处冠层选取。非森林主要包括农田和草地，其端元分别从农田去和草地区域选取，其端元波谱特征曲线见图1c，
[0058]
利用以上端元波谱和线性像元分解方法提取出光照森林、阴影森林、农田和草地在本层级的丰度，最后归并光照森林和阴影森林为本层级森林丰度，农田和草地为本层级非森林丰度。
[0059]
步骤四、利用丰度值和乘法模型计算森林在的覆盖度。
[0060]
研究区森林的综合覆盖度是第三级森林丰度、第二级植被丰度和第一级非水体丰度的函数。在第三级森林-非森林层级中，通过lsma提取的森林丰度为f
forest
，森林在第二层级属于植被类型，植被类型丰度通过lsma模型计算为f
vegetation
，而植被在第一层级属于非水体类型，非水体丰度通过lsma模型计算为f
non-water
，通过乘法模型计算。因此，研究区的森林综合覆盖度通过乘法模型进行计算，计算公式为：f＝f
non-water
,*f
vegetation
*f
forest
，计算结果如图2所示。
[0061]
步骤五、对hmm方法的精度进行检验。精度检验主要通过两种方法实现：
[0062]
a)以官方统计数据作为验证数据计算hmm方法的精度：
[0063]
通过统计数据获取浙江省森林覆盖度，并以此数据为验证数据，利用hmm方法获取的研究区森林覆盖度与验证数据进行比较，分析hmm方法的整体精度和计算误差。浙江省林业局统计报告显示2019年浙江森林覆盖度为61.15％，而hmm方法计算的森林综合覆盖度为60.09％，计算误差仅为1.57％。
[0064]
b)以高空间分辨数据计算植被覆盖度作为验证数据分析hmm方法的精度。
[0065]
在研究区选取一块1.5km*1.2km、包含本研究所有土地覆盖类型的区域为试验区，下载覆盖试验区的2米分辨率的高分1号全波段/多光谱遥感数据(gf-1pms)，利用支持向量机对验证区的地表进行分类，并依据分类结果计算试验区森林的覆盖度作为验证数据，
[0066]
①
对gf-1pms影像进行辐射校正、大气校正、几何校正，并进行多光谱和全波段的融合，获取2米分辨率的影像。在envi中，启动toolbox/image sharpening/nndiffuse pan sharpening，分别选择正射校正好的多光谱和全图像，选择输出路径及文件名，融合后为2m分辨率多光谱影像。
[0067]
②
基于融合后的影像利用支持向v机进行监督分类，分类类型与hmm分类一致，包括水体、森林、裸地、高反射物体、普通建筑物、农田、草地等。在envi中打开pms融合影像，根据影像分类数量新建相应8类感兴趣区roi文件，左键双击新建的roi文件，弹出工具框，修改roi名称，选取代表颜，在影像上选取对应的区域，呼出不规则闭合多边形，生成水体、森林、裸地、高反射物、普通建筑物、农田、草地的roi样本，并保存。在软件工具栏toolbox中
选择“raster management-masking-building mask”，选择要做掩膜的pms融合影像，制作掩膜文件。在软件由此工具栏toolbox中，点击“classification-supervised classification”选择support vector machine classification”,在弹出的对话框中选择融合后的影像、掩膜文件和分类样本，并输出文件名完成监督分类，如图3所示，根据分类结果将所有地物分为森林和非森林两大类。
[0068]
③
基于监督分类结果，计算验证区的森林丰度。在试验区按照wfv的分辨率划分16*16米的网格，每个网格内包含8*8＝64个空间分辨率为2m的像元。图4所示，由于每个小像元的类型已经被识别为森林或非森林，因此被划分为森林类型的小像元个数占整个大像元的比例被计算出来，即为试验区的精细森林丰度(覆盖率)。
[0069]
④
利用融合影像计算出的森林丰度分析hmm方法的精度。对比基于高精度影像和svm结合的方法与hmm方法获取的森林丰度可以发现其相关度达到0.95。空间分布上，森林内部的相关性高，如图5所示，结果误差主要分布于森林边缘地带。
[0070]
如图6所示，经统计发现43.12％的像元的森林覆盖度完全一致，误差在0.1以内的像元占比为25.30％，误差在0.2以内像元为15.9，误差为0.3以内的占9.08％，误差在0.4以内的像元占比为1.46％，误差超过0.4的像元仅有2.61％。
[0071]
本发明提供的基于分层乘法模型和高分1号数据的森林覆盖度提取方法，充分利用了像元分解模型降低了森林覆盖度提取过程中混合像元的不利影响，同时利用乘法传输模型解决了像元分解过程中多光谱影像波段不足的问题，不仅极大的提高了森林覆盖度提取精度，而且扩大了多光谱影像在森林覆盖度提取中的应用。国产高分1号遥感数据的应用为大范围大尺度的森林覆盖度提取提供了有力的数据支撑。
[0072]
还能结合时间序列数据，准确计算林业覆盖度，监测森林扰动信息，分析影响森林生态变化的因素。

技术特征：

1.基于分层乘法模型和高分1号数据的森林覆盖度提取方法，其特征在于：具体包括如下步骤：步骤一、影像预处理：获取目标对象的gf1数据，借助gf1传感器的辐射定标参数对原始影像进行辐射定标、大气校正和几何校正；步骤二、土地覆盖类型分层：根据地物波谱差异的大小分为第一层水体-非水体层、第二层植被-非植被层、第三层森林和非森林层；步骤三、在各层级内进行像元分解，提取各层级对应元素的丰度：步骤3.1、在第一层内分解出水体和非水体：选取水体端元和非水体端元，利用线性光谱混合模型lsma模型和最小二乘法计算出水体和非水体的丰度；通过计算归一化水指数ndwi，减少高含水量物质的干扰；ndwi计算公式为：(绿波-近红外)/(绿波+近红外)，只有ndwi指数小于零的非水体像元才能在下一级中被进一步分解；步骤3.2在第二层内分解出植被和非植被：植被分为光照植被和阴影植被，植被端元分别由光照条件下的植被和阴影下的植被选取；非植被分为裸地、高光效的不透水层和不透水层三类，非植被端元分别从对应的裸地、高光建筑和普通不透水层选取；利用lsma模型分别获得光照植被、阴影植被、裸地、高光不透水层和一般不透水层的丰度，进而叠加光照植被和阴影植被为植被，叠加其余子类为非植被；步骤3.3、在第三层内分解出森林和非森林：森林端元分为光照森林和阴影森林，与非森林端元利用lsma模型提取出各自丰度，然后光照森林和阴影森林叠加为森林类型，其余叠加为非森林类；步骤四、利用丰度值和分层乘法模型hmm计算森林的覆盖度：森林覆盖度是森林在最低层级的丰度与其所属类型在上一层级的丰度的函数，并递推至最高层级；根据lsma模型计算结果，在第三层级森林丰度为f3，森林所属的植被类型在第二级的丰度为f2，植被所属的非水体类型在第一级的丰度为f1，根据分层乘法模型，森林在第一级的丰度即森林覆盖度的计算公式为：f＝f1*f2*f3。2.如权利要求1所述的基于分层乘法模型和高分1号数据的森林覆盖度提取方法，其特征在于：步骤一所述的辐射定标将原始dn值转化为像元辐射亮度值，大气校正是利用flaash模型将辐射亮度值转化为地表实际反射率，几何校正首先选取同名地物点，采用最邻近像元法进行重采样，校正误差控制在0.2个像元。3.如权利要求1所述的基于分层乘法模型和高分1号数据的森林覆盖度提取方法，其特征在于：步骤三所述的像元分解是从多地物光谱混合的数据中提取各种地物以及各成分的比例，线性光谱混合模型lsma定义像元在某一光谱波段的反射率是由构成像元的基本组分的反射率以其所占像元面积比例为权重系数的线性组合；lsma计算方法为：lsma计算方法为：f
ki
≥1
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)；
式中r
iλ
是λ波段上的i像元的波谱反射；n是地物数量；f
ki
是像元i内端元k所占的面积比例；r
kλ
是λ波段上端元k的波谱反射率；ε
iλ
是λ波段上像元i的误差；丰度f
ki
通过最小二乘法计算，计算公式如下：式中ε
iλ
,r
iλ
,r
iλ
分别是λ波段上像元i的误差值，观察波谱值和估算波谱值；分别按照步骤二所述三个层次选取端元，并利用lsma模型在各层级进行像元分解提取对应地物的丰度。

技术总结

本发明涉及基于分层乘法模型和高分1号数据的森林覆盖度提取方法。本发明包括利用分层乘法模型从多光谱遥感数据中准确、高效的提取森林覆盖度。在对遥感数据预处理的基础上，首先按照地物波谱特征差异将地物类型划分为三个主要层级。其次在每一层级中利用线性像元分解模型进行地物丰度计算，最后通过乘法模型计算林业植被在第一级地物中的丰度。本发明创建了基于乘法模型的像元分解算法，提高了基于多光谱影像提取地物丰度值的适用性和准确性。解决了多光谱波段不足的限制，针对像元分解精度难以验证的问题，提出了利用高精度影像进行分类、并以分类结果计算大尺度上各端元的丰度以验证像元分解精确度，提高了像元分解验证的准确性。确性。确性。