CN112381133A - 一种森林地上生物量遥感估测中合理样本数的确定方法 - Google Patents

Abstract

本发明公开了一种森林地上生物量遥感估测中合理样本数的确定方法，以高山松、云冷杉2个树种为研究对象，结合Landsat 8影像和地面调查样地数据，基于地统计学半方差函数理论，以k‑NN遥感估测模型为案例，以均方根误差为评价指标，计算得到了高山松和云冷杉模型精度最佳时的样本数量大小，确定了一个合理的样本数量范围。本发明为了探讨定量遥感反演中样本数量的不确定性，基于地统计学半变异函数理论和k‑NN模型，分析了样本数量在森林地上生物量估测中的不确定性影响，为传统定量遥感反演中样本数量难以确定的问题提供了参考依据。

Description

一种森林地上生物量遥感估测中合理样本数的确定方法

技术领域

本发明涉及一种森林生物量遥感估测中样本数量的不确定性分析技术，具体为一种森林地上生物量遥感估测中合理样本数的确定方法。

背景技术

森林生物量是森林是评价森林生态系统生产力、陆地生态系统功能和可持续性的重要指标之一。随着遥感技术的飞速发展，利用多源遥感数据替代传统调查方法进行定量化森林地上生物量，不仅能获取森林资源的数量、空间分布和动态变化信息，而且可以结合各种模型和样地调查实现森林参数的定量反演(何兴元等，2018)，在满足不同尺度森林资源、生态过程监测分析的需求的同时节约调查成本。近年来，国内外对用光学遥感影像信息数据进行地上生物量反演都进行了大量的研究(汤旭光，2012；Beaudoin，2014；谢福明，2018)。其中，使用基于非参数的方法，如k-最近邻法、支持向量机和人工神经网络进行地上生物量反演成为研究热点之一(Wang，2009)。同时，基于遥感反演森林生态系统生物量的不确定性研究也成为学者们关注的重点(Wang，2011；Lu，2012)，其不确定性的主要来源包括研究区的典型性，遥感数据和地面样地数据获取中存在的误差，模型的类型、参数和变量的选择，样本数量大小以及样本的空间分布以及模型的评价指标(Phillips，2000)。生物量遥感估测不确定性分析目的在于，准确获取不确定的来源，然后从数据采集和模型方法优化等方面做出相应的改进，降低不确定性因素的影响，提高遥感估测反演产品的可靠性和实用性。其中，地面实测样本数据的不确性一项关键的不确定性因素，地面观测的数据的采集是一项费时费力又必不可少的工作，如何采集合理的样本数是一个值得探究的问题。

当前，国内外对自由样本数量研究尚少：傅煜等(2015)基于单木生物量建模样本量对估计区域生物量的不确定性分析认为，增加建模数据量能有效提高生物量模型的估计精度、准确度以及工作效率，并降低不确定性。在实际应用中，森林资源调查数据往往难以获取，在样本数量有限的情况下模型会发生“过学习”现象。而采用非参数化方法估测森林地上生物量能有效解决这一问题(韩宗涛，2018)。吴超凡基于区域森林生物量遥感估测与应用研究表明，样本数量的增加有利于提高建模精度，特别是对于支持向量算法；但对于PLS和KNN算法表现出的精度变化也反映出，并非样本越多，精度越高，不同的估算方法需要寻找最合适的样本条件(吴超凡，2015)。由此可见，样本容量越大，样本的代表性就越好，建模的精度也就越高，但耗费的人力、物力也就越多，甚至难以完成。相反，如果样本量太小，其建模精度不能够满足需求。本发明旨在分析样本数量在森林地上生物量估测中的不确定性影响，对解决传统定量遥感反演中样本数量难以确定的问题具有重要的科学价值。

发明内容

针对现有技术的不足，本发明的目的在于提供一种基于遥感数据的森林地上生物量估测中合理样本数的确定方法，实现生物量的科学准确估测。

为实现上述目的，本发明采用的技术方案是：一种森林地上生物量遥感估测中合理样本数的确定方法，包括如下步骤：

步骤1、获取美国陆地卫星Landsat 8OLI图像，经过辐射定标、大气校正、几何校正，得到地表反射率数据集；

步骤2、进行野外样地调查，对伐倒木进行取样称重，带回实验室烘干称重，采用不同的函数形式拟合得到单木的地上生物量估测模型；

在本发明中，单木生物量模型也可直接引用其他学者已经建立的单木生物量模型；

步骤3、根据森林资源二类调查角规控制样地数据，计算样地水平的地上生物量；

步骤4、以步骤1和步骤3得到的因子作为变量，构建模型数据集；

步骤5、构建森林地上生物量的非参数估测模型—k-NN估测模型，在MATLAB软件平台上编程实现k-NN地上生物量估测反演算法；

步骤6、基于k-NN模型和半变异函数理论分析探讨样本数量的不确定性分析；

步骤7、模型检验，利用均方根误差RMSE作为模型精度的评价指标，均方根误差的计算方法见公式(1)，

式中，y_i和

分别为第i个样本的实测值与模型预测值；N为样本个数

步骤8、分析建模结果数据，确定模型精度最佳时的样本数量，针对本案例确定合理的样本数范围。

进一步的：步骤1中，所述的遥感数据来自2013年2月11日发射的Landsat 8卫星，旨在获取农业、教育、商业、科学和政府中有价值的数据和图像，Landsat 8卫星包含OLI(Operational Land Imager陆地成像仪)和TIRS(Thermal Infrared Sensor热红外传感器)两种传感器；OLI包括了ETM+的所有波段，为了避免大气吸收部分特征，OLI对波段进行了重新调整，主要包括：Band 5(0.84～0.88μm)，排除了0.82μm处水汽吸收特征；全色波段Band 8波段范围较窄，这种方式可以在全色图像上更好区分植被和无植被特征；新增两个波段：海蓝波段(band 1(0.43～0.45μm)主要应用海岸带观测；短波红外波段，又称卷云波段(band 9(1.36～1.39μm)包含水汽强吸收特征，可用于云检测；近红外band 5和短波红外band 9与MODIS对应的波段更加接近。

进一步的：步骤2中，所述的单木地上生物量模型仅针对高山松树种，高山松样木数据记录了不同龄组下包括幼龄林、中龄林、近熟林、成熟林、过熟林116株高山松胸径(DBH)和树高(H)，并测定了树干、树皮、树叶、树枝、树冠的生物量，用于拟合高山松地上生物量计算模型，利用公式(2)和(3)计算得到高山松单木地上生物量和云冷杉单木地上生物量；生物量调查参照胥辉等生物量测定方法(胥辉等，2002)，而云冷杉单木地上生物量模型则引用其他学者的研究结果(赵菡，2017)；

AGB＝0.0955(DBH²H)^0.8329 (2)

AGB＝0.048(DBH²H)^0.880 (3)

进一步的：步骤3中，所述样地数据是2016年度香格里拉市森林资源二类调查的角规控制样地，每个角规控制检尺样地记录了小班地类、优势树种、郁闭度、样木胸径和树高等因子，同时计算了每个样地的平均树高、平均胸径、总断面积和蓄积量等信息。

进一步的：步骤5中，所述k-NN模型，其原理描述如下，在k-NN的专业术语中，将待测变量及其特征变量的观测值样本指定为参考集，将待测变量的预测集指定为目标集，特征变量定义的空间成为特征空间，对于诸如生物量或蓄积量等连续性变量M在像元p上的预测值mp的计算方法如下，

式(4)中：mi为变量M参考样地点i上的实测值；k为计算预测值mp时考虑的近邻个数；w_ip为像元权重值，其计算如下，

式(5)中：i是参考集样本；p是目标集像元；p_j是与参考集样本j对应的样本；

为距离分解因子；k，t为常量，一般通过实验反复测试选取最佳值；{i₁(p)，…，i_k(p)}是与待测像元p在特征空间上最相似的k个参考集样本，特征变量空间相似度由d_pi，p度量，其计算方法如下，

式(6)中：f_l，pj和f_l，p分别为参考集和目标集样本对应的遥感影像光谱波段及其派生因子等特征变量；n_f为特征变量个数；p为目标集像元；pi为参考集样本i对应的像元。

进一步的：步骤6中，所述半变异函数的拟合采用球状模型，具体描述如下，半方差函数的理论模型分为有基台值模型(球状、高斯、指数、线性有基台值和纯块金效应模型)、无基台值模型(幂函数、抛物线、线性无基台值模型)、孔穴效应模型，本次研究半方差函数的拟合使用球状模型，其计算公式见方程式(7)

式(7)中，C_o为块金常数，C为拱高，C_o+C为基台值，a为变程；基于球状模型，在进行最佳的样本数量求解时，h即为样本数量(h＝s+20)；γ(h)为k-NN模型估测精度；C₀为样本数h＝20(s＝0)时的k-NN模型估测精度；C是拱高，即模型精度变化率；C_o+C为样本数量达到最佳时所对应的模型估测精度；a为最佳的样本数量。在球状模型参数估计时，记γ(h)＝Y，h＝X₁，h³＝X₂，c₀＝B₀，

变换后的线性模型如方程式(8)，利用最小二乘法对模型公式(8)参数(B₀，B₁，B₂)进行求解，即可求得不同树种对应的最佳样本数

Y＝B₀+B₁X₁+B₂X₂ (8)

进一步的：步骤7中，所述的精度评价方法采用留一法交叉验证，即对于N个样本量，每次从N个样本中抽出1个样本作为测试集，利用剩余的N-1个样本作为参考集，重复N次循环，直至结束；本研究将N个样本的模型预测值

与对应样本的实测值(y_i)进行统计分析，利用均方根误差RMSE来检验模型的精度，

式中，y_i和

分别为第i个样本的实测值与模型预测值；

本项发明具有技术创新点：为了探讨定量遥感反演中样本数量的不确定性，本发明基于地统计学半变异函数理论和k-NN模型，分析了样本数量在森林地上生物量估测中的不确定性影响，为传统定量遥感反演中样本数量难以确定的问题提供了参考依据。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实验区和样地分布图；

图2是本发明200次实验结果的波动性(高山松，样本数量＝65)；

图3是本发明200次实验结果的波动性(云冷杉，样本数量＝70)；

图4是本发明高山松地上生物量估测精度随样本数量大小的变化；

图5是本发明云冷杉地上生物量估测精度随样本数量大小的变化；

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实例仅仅室本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

一种森林地上生物量遥感估测中合理样本数的确定方法，包括以下步骤：

1)研究区概况

研究区香格里拉市原名中甸县，藏语称“建塘”，隶属于云南省迪庆藏族自治州，位于云南省西北部、青藏高原横断山区腹地，纬度范围为26°52′11.44″～28°50′59.57″，经度范围为99°23′6.08″～100°18′29.15″(图1)，总面积114.2万hm²，东、南、西三面被金沙江环绕，是滇、川、藏三省区交界地，也是世界自然遗产“三江并流”景区所在地。研究区地势高耸，热量不足，气温偏低，海拔为1503～5545m，多年平均气温为5.5℃，历年平均降水量为618.4mm，平均降雪日为35.7d，年日照率为40％～50％，属山地寒温带季风气候。境内密集的金沙江水系支流、冰雪融水和高原湖泊等水资源以及以棕壤、红壤为主的森林土壤类型孕育了丰富的植物资源。森林植被面积大，覆盖率高，南北差异分布明显，主要分布有10种植被类型，常见的树种有云杉Picea asperata，冷杉Abies fabri，高山松Pinus densata，云南松Pinus yunnanensis，高山栎Quercus semicarpifolia等

2)遥感数据的采集与处理

遥感数据来自2013年2月11日发射的Landsat 8卫星，包含OLI(Operational LandImager陆地成像仪)和TIRS(Thermal Infrared Sensor热红外传感器)两种传感器。从地理空间数据云(http：//www.gscloud.cn/)里经过反复筛选，筛选出数据质量较高的三景影像(参数见表1)，利用软件ENVI 5.3对卫星影像进行辐射定标、大气校正(FLAASH)等预处理，获取Band1～Band7的信息作为本研究中的遥感特征变量。

表1 Landsat8 OLI影像参数

3)样地生物量数据的采集与处理

本研究中的地面样地数据是2016年度香格里拉市森林资源二类调查的角规控制样地，每个角规控制检尺样地记录了小班地类、优势树种、郁闭度、单木胸径和树高等因子，同时计算了每个样地的平均树高、平均胸径、总断面积和蓄积量等信息。样地地上生物量的计算分为两个环节：首先，利用样地的平均树高和平均胸径来计算样地的单木平均地上生物量；其次，依据公式9计算每公顷面积上的林木株数，结合平均单木地上生物量求得样地的地上生物量。

本研究中的样地数据共计191个样地，其中高山松91个，云冷杉100个。表2详细记录不同树种样地地上生物量的最小值、最大值、均值和标准差，最大值为224.12Mg/ha(云冷杉)，最小值为3.36Mg/ha(高山松)，样地在研究区范围均匀分布(图1)，是本研究中重要的地面参考数据。

式中，AGB表示地上生物量；DBH表示胸径；H表示树高；Z_j为第j径阶(设林分中林木共有K个径阶)的计数木株数，g_j为该径阶中值的断面积；F_g为断面积系数；N为求得的每公顷林木株数。

表2森林实测样地地上生物量描述统计分析表

4)模型方法的构建

本发明基于k-NN模型设计实验，分析不同数量样本对模型运行结果的影响。为了便于描述，记x为样本数量，Z(x)为k-NN模型估测精度与样本数量间的函数。设计实验：每次实验从总体样本中随机抽取x个样本，样本数量x(x≥20)从20开始计数，每增加一个样本，运行k-NN模型并记录模型的估测精度，重复200次实验，直到样本数量x等于总体的样本数，结束实验。设置200次的循环实验是为了避免单次实验的随机性，提高实验的科学性与可信度。最后，通过绘制k-NN模型估测精度Z(x)随样本数量x变化的散点图，分析样本数量对模型精度的影响效益，确定合理的样本数量范围。此外，除了定性的分析样本数量对模型估测结果的影响之外，还通过定量的方法分析估算了基于k-NN模型的最佳样本数量。

5)结果分析

本技术发明的主要目的是分析样本数量对模型精度的影响情况，其前提是保证样本数量的充足，而对k-NN模型的估测精度不作要求。所以，本次实验用高山松(91个样本)和云冷杉(100个样本)样地数据的实测地上生物量作为响应变量，选取OLI影像的7个原始波段光谱反射值作为特征变量，建立了基于k-NN模型的样本数量不确定性分析技术流程。采用留一法交叉验证精度评价方法，用均方根误差RMSE记录模型的预测精度。本发明中k-NN模型的参数的选择为：距离度量方式采用欧式距离，k值等于5，t值等于2。

①不同样本数量下模型精度的统计分析

本发明为避免单次实验的随机性误差，每次从总体样本中抽取确定的样本数量N，有放回地随机抽取200次，采用留一法交叉验证评价每一次抽样的模型精度，并求得200次实验结果得均值作为样本数量大小为N下的模型精度。图2描绘的是高山松样本数量为65时的200次实验的结果，估测精度指标RMSE最大值40.44t/ha，最小值30.35t/ha，平均值36.36t/ha；图3描绘的是云冷杉样本数量为70时的200次实验的结果，估测精度指标RMSE最大值53.59t/ha，最小值39.92t/ha，平均值47.78t/ha。从散点图数据分析显示，实验次数与精度指标RMSE之间的决定系数均R2小于0.00，这一数据充分展示了随机抽样中的随机性误差。

表3和表4是分别是不同样本数量下高山松和云冷杉的模型估测精度均方根误差RMSE，数据已排除抽样中随机性误差。

表3不同样本数量的模型估测精度RMSE(高山松)

表4不同样本数量的模型估测精度RMSE(云冷杉)

②样本数量的不确定性分析

图4是高山松地上生物量k-NN模型估测精度随样本数量的变化及其线性拟合残差分布图。数据分析显示，高山松样本数量＜60时的模型估测精度方差是样本数量≥60时的8.8倍，当样本数<60时，模型精度不稳定，波动大；当样本数＞60时，模型精度波动减小，并趋向于稳定。利用最小二乘法求解出地上生物量估测模型精度最佳时的样本数量为72。

图5分别是云冷杉地上生物量k-NN模型估测精度随样本数量的变化及其线性拟合残差分布图。数据分析显示，云冷杉样本数量＜60时的模型估测精度方差是样本数量≥60时的2.7倍，当样本数<60时，模型精度不稳定，波动大；当样本数＞60时，模型精度波动减小，并趋向于稳定。利用最小二乘法求解出地上生物量估测模型精度最佳时的样本数量为65。

综合上述的分析，总体上呈现出相同的规律，高山松和云冷杉地上生物量估测中样本数量的合理范围均介于60到70之间，即在基于k-NN模型的地上生物量遥感估测中，为了保证模型精度的准确性，要保证充足的样本数量，即样本数量在60个及以上可以保证模型估测结果的科学性。这一结论也将进一步论证本文采用的样本数的合理性。

本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

最后所应说明的是：以上实施例仅用以说明而非限制本发明的技术方案，尽管参照上述实施例对本发明进行了详细说明，本领域的普通技术人员应该理解：依然可以对本发明进行修改或者等同替换，而不脱离本发明的精神和范围的任何修改或局部替换，其均应涵盖在本发明的权利要求范围当中。

Claims (7)

1.一种森林地上生物量遥感估测中合理样本数的确定方法，其特征在于，以高山松、云冷杉2个树种为研究对象，步骤包括：

步骤1、获取美国陆地卫星Landsat 8OLI图像，经过辐射定标、大气校正、几何校正，得到地表反射率数据集；

步骤2、进行野外样地调查，对伐倒木进行取样称重，带回实验室烘干称重，采用不同的函数形式拟合得到单木的地上生物量估测模型；

步骤3、根据森林资源二类调查角规控制样地数据，计算样地水平的地上生物量；

步骤4、以步骤1和步骤3得到的因子作为变量，构建模型数据集；

步骤5、构建森林地上生物量的非参数估测模型k-NN估测模型，在MATLAB软件平台上编程实现k-NN地上生物量估测反演算法；

步骤6、基于k-NN估测模型和半变异函数理论分析探讨样本数量的不确定性分析；

步骤7、模型检验，采用交叉验证留一法，利用均方根误差RMSE作为模型精度的评价指标，均方根误差的计算方法见公式(1)，

式中，y_i和

分别为第i个样本的实测值与模型预测值；N为样本个数；

步骤8、分析建模结果数据，确定模型精度最佳时的样本数量，确定合理的样本数范围。

2.根据权利要求1所述的方法，其特征在于：步骤1中，所述遥感数据来自2013年2月11日发射的Landsat 8卫星，Landsat 8卫星包含OLI和TIRS两种传感器，OLI包括了ETM+的所有波段，为了避免大气吸收部分特征，OLI对波段进行了重新调整，主要包括：Band 50.84～0.88μm，排除了0.82μm处水汽吸收特征；全色波段Band 8波段范围较窄，这种方式可以在全色图像上更好区分植被和无植被特征；新增两个波段：海蓝波段band 1 0.43～0.45μm主要应用海岸带观测；短波红外波段，又称卷云波段band 9 1.36～1.39μm包含水汽强吸收特征，可用于云检测；近红外band 5和短波红外band 9与MODIS对应的波段更加接近。

3.根据权利要求1所述的方法，其特征在于：步骤2中，所述单木的地上生物量模型仅针对高山松树种，高山松样木数据记录了不同龄组下包括幼龄林、中龄林、近熟林、成熟林、过熟林116株高山松胸径(DBH)和树高(H)，并测定了树干、树皮、树叶、树枝、树冠的生物量，用于拟合高山松地上生物量计算模型，利用公式(2)和(3)计算得到高山松单木地上生物量和云冷杉单木地上生物量；

AGB＝0.0955(DBH²H)^0.8329 (2)

AGB＝0.048(DBH²H)^0.880 (3)

4.根据权利要求1所述的方法，其特征在于：步骤3中，所述样地数据是2016年度香格里拉市森林资源二类调查的角规控制样地，每个角规控制检尺样地记录了小班地类、优势树种、郁闭度、样木胸径和树高等因子，同时计算了每个样地的平均树高、平均胸径、总断面积和蓄积量信息。

5.根据权利要求1所述的方法，其特征在于：步骤5中，所述k-NN估测模型，其原理描述如下，在k-NN的专业术语中，将待测变量及其特征变量的观测值样本指定为参考集，将待测变量的预测集指定为目标集，特征变量定义的空间成为特征空间，对于诸如生物量或蓄积量等连续性变量M在像元p上的预测值mp的计算方法如下，

式(4)中：mi为变量M参考样地点i上的实测值；k为计算预测值mp时考虑的近邻个数；w_ip为像元权重值，其计算如下，

式(5)中：i是参考集样本；p是目标集像元；p_j是与参考集样本j对应的样本；

为距离分解因子；k，t为常量，一般通过实验反复测试选取最佳值；{i₁(p)，…，i_k(p)}是与待测像元p在特征空间上最相似的k个参考集样本，特征变量空间相似度由

度量，其计算方法如下，

式(6)中：

和f_l，p分别为参考集和目标集样本对应的遥感影像光谱波段及其派生因子等特征变量；n_f为特征变量个数；p为目标集像元；pi为参考集样本i对应的像元。

6.根据权利要求1所述的方法，其特征在于：步骤6中，所述半变异函数的拟合采用球状模型，具体描述如下，半方差函数的理论模型分为有基台值模型、无基台值模型和孔穴效应模型，本次研究半方差函数的拟合使用球状模型，其计算公式见方程式(7)，

式(7)中，C_o为块金常数，C为拱高，C_o+C为基台值，a为变程；基于球状模型，在进行最佳的样本数量求解时，h即为样本数量(h＝s+20)；γ(h)为k-NN模型估测精度；C₀为样本数h＝20(s＝0)时的k-NN模型估测精度；C是拱高，即模型精度变化率；C_o+C为样本数量达到最佳时所对应的模型估测精度；a为最佳的样本数量。在球状模型参数估计时，记γ(h)＝Y，h＝X₁，h³＝X₂，c₀＝B₀，

变换后的线性模型如方程式(8)，利用最小二乘法对模型公式(8)参数(B₀，B₁，B₂)进行求解，即可求得不同树种对应的最佳样本数；

Y＝B₀+B₁X₁+B₂X₂ (8)。

7.根据权利要求1所述的方法，其特征在于：步骤7中，所述精度评价方法采用留一法交叉验证，即对于N个样本量，每次从N个样本中抽出1个样本作为测试集，利用剩余的N-1个样本作为参考集，重复N次循环，直至结束，本研究将N个样本的模型预测值

与对应样本的实测值(y_i)进行统计分析，利用均方根误差RMSE来检验模型的精度，均方根误差的计算方法见公式(1)

式中，y_i和

分别为第i个样本的实测值与模型预测值。

Images