CN105913016A - 一种通过条带LiDAR数据升尺度估算森林生物量的方法 - Google Patents

Abstract

本发明公开了一种通过条带LiDAR数据升尺度估算森林生物量的方法，以江苏南部丘陵地区的亚热带天然次生林为方法实施对象，首先从LiDAR条带数据中提取9个特征变量并与地面估算的生物量结合反演条带内的生物量连续分布信息；然后，在条带反演结果区域内采样并与覆盖整个研究区域的Landsat OLI影像特征变量结合，从而升尺度估算整个研究区域的生物量。从而在充分挖掘遥感数据特征信息的基础上，通过一个条带的LiDAR数据在降低林场级生物量估算成本的同时，也提升了遥感方法在该尺度上的生物量估算精度。

Description

一种通过条带LiDAR数据升尺度估算森林生物量的方法

技术领域

本发明属于林业调查、动态监测及生物多样性等技术领域，涉及一种通过条带LiDAR数据升尺度估算森林生物量的方法。

背景技术

森林生态系统是陆地生物圈的主体，其生物量约占整个陆地生态系统的85％，在减缓全球气候变化中具有不可替代的作用。森林在生长过程中通过同化作用吸收大气中的CO₂，并以生物量形式将其长期固定。亚热带树种丰富且森林生产力高，不仅对区域生态环境有很大影响，而且对维持全球碳平衡也具有重要作用。传统的生物量调查方法耗时费力，且只能获得有限的“点”上信息。而遥感技术却能够准确、快速地获取各个尺度的“面”上连续分布的森林信息，在森林资源动态监测及森林碳储量反演方面具有重要作用。

多光谱遥感数据可以获取连续、实时的森林生物物理特性的水平空间分布信息，且免费获取，但易“饱和”且难于获取森林的结构参数；而LiDAR发射的激光脉冲可以穿透植被冠层获得其三维结构和能量信息，但是覆盖整个研究区的LiDAR数据通常获取成本高，难以在生产上应用。Guo等在《Science China(Earth Sciences)》2010年第1期发表了“Estimating forest aboveground biomass using HJ–A satellite and ICESat GLASwaveform data”，利用环境卫星多光谱数据和ICESat/GLAS波形数据，以塔河林场为研究区反演了地上生物量，结果表明80％的验证点生物量误差在20t·hm^-2以内。但该要求LiDAR数据与多光谱数据覆盖整个研究区，成本较高，工作量较大，不适合在更大尺度上推广和应用。黄克标等在《遥感学报》2013年第1期发表了“基于ICESat GLAS的云南省森林地上生物量反演”，结合机载LiDAR数据对星载GLAS光斑范围内的森林地上生物量进行估测，并利用MODIS以及MERIS土地覆盖产品进行了云南省森林地上生物量升尺度估算，其构建的生物量估算模型的决定系数(R²)为0.52，均方根误差(RMSE)为31Mg·hm^-2。庞勇等在《资源科学》2011年第10期发表“基于遥感的湄公河次区域森林地上生物量分析”，通过非连续的机载LiDAR和星载LiDAR GLAS数据以及连续的光学遥感数据MERIS对大湄公河次区域的森林地上生物量进行连续分布制图。连续覆盖生物量估计的模型总体误差为34t·hm^-2，升尺度模型的相关系数为0.7。以上研究证明了非连续LiDAR数据进行森林参数升尺度反演和连续分布制图的可行性及有效性，但其升尺度模型自变量相对单一(如MERIS数据只提供了NDVI作为反演参数)且制图分辨率较低(300m)。同时，以上研究及应用都未在北亚热带地区进行实施。

发明内容

发明目的：针对现有技术中存在的不足，本发明的目的是提供一种通过条带LiDAR数据升尺度估算森林生物量的方法，有效提高精度，降低成本。

技术方案：为了实现上述发明目的，本发明采用的技术方案为：

一种通过条带LiDAR数据升尺度估算森林生物量的方法，包括以下步骤：

1)采用Riegl LMS-Q680i LiDAR传感器获取待检测区域内的LiDAR点云数据，传感器记录完整的激光脉冲返回波形信息，其时间采样间隔为1ns；

2)在待检测区域范围内设置多个方形样地，样地调查过程中，对于胸径大于5cm的树，逐一测定单木的树种、胸径、树高和枝下高以及冠幅；根据单木调查数据汇总样地尺度的相关森林参数，包括样地尺度上的单位面积地上生物量、地下生物量；生物量信息通过异速生长方程计算单木的生物量，并汇总得到每块样地单位面积地上生物量W_A以及地下生物量W_B；

3)根据LiDAR点云数据计算9个特征变量：冠层覆盖度变量c；5个高度百分位数变量：h₂₅、h₅₀、h₆₀、h₇₅、h₉₅；4个高度变量：最大高度h_max、最小高度h_min、高度变异系数h_cv以及平均高度h_mean；同时，提取出5组OLI多光谱特征变量：包括6个原始单波段特征、10个波段组合特征、10个信息增强组、18个植被指数特征以及9个纹理特征；其中，OLI多光谱特征变量已由像元尺度经均值计算转化为样地尺度；

4)利用样地GPS坐标点位置，基于9景包含LiDAR特征变量信息的影像和53景包含OLI特征变量信息的影像，分别提取样地中心位置的像元值，并借助pearson’s相关系数矩阵分析特征变量之间的相关性，去除相关性高的特征变量；pearson’s相关性分析：

$r=\frac{\underset{i=1}{\overset{n}{\Σ}}(x_i-{\stackrel{\‾}{x}}_i)(y_i-{\stackrel{\‾}{y}}_i)}{\sqrt{\underset{i=1}{\overset{n}{\Σ}}{(x_i-{\stackrel{\‾}{x}}_i)}^2\·\underset{i=1}{\overset{n}{\Σ}}{(y_i-{\stackrel{\‾}{y}}_i)}^2}}---\left(1\right)$

式中：x_i为地面实测的某林分特征；y_i为某特征变量；为x_i的平均值；为y_i的平均值；

将地面实测汇总的生物量信息作为因变量，遥感方法提取的特征变量作为自变量，建立多元回归模型；运用逐步进入法和检验决定系数R²的变化情况来选择进入模型的合适变量；如果有自变量使统计量F值过小并且T检验达不到显著水平，则予以剔除；F值较大且T检验达到显著水平则得以进入；通过控制因子κ进一步筛选模型；κ通过PCA分析计算相关关系矩阵得到，其值低于30则表明自变量之间相关性较低；采用决定系数R²、均方根误差RMSE和相对均方根误差rRMSE评价回归模型的精度，计算公式为：

$R^2=1-\frac{\underset{i=1}{\overset{n}{\Σ}}{(x_i-{\hat{x}}_i)}^2}{\underset{i=1}{\overset{n}{\Σ}}{(x_i-{\stackrel{\‾}{x}}_i)}^2}---\left(2\right)$

式中：x_i为地面实测的某林分特征；为x_i的平均值；为模型估算的某林分特征；n为样地数量；

$RMSE=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{(x_i-{\hat{x}}_i)}^2}---\left(3\right)$

$rRMSE=\frac{RMSE}{{\stackrel{\‾}{x}}_i}\×100\%---\left(4\right)$

森林生物量升尺度反演步骤如下：先以重叠区样地LiDAR特征变量为自变量，样地实测生物量为因变量构建LiDAR生物量反演模型，再利用LiDAR模型计算出的重叠区所有像元的生物量，以其作为因变量，与重叠区的OLI特征变量进行多元线性拟合，建立重叠区生物量反演模型。

Riegl LMS-Q680i LiDAR传感器的激光发射脉冲所在波段为近红外，波长为1055nm，脉冲发射频率为360kHz，扫描频率为112线·s^-1，扫描角为±30°；遥感平台的飞行高度为900m，飞行速度为55m·s^-1，旁向重叠度为60％；地面光斑直径约为25cm，平均地面光斑距离约为0.48m。

Riegl LMS-Q680i LiDAR传感器最终获取的数据格式为LAS 1.3标准格式，每个激光返回点内包含三维坐标值、强度值和返回类型信息。

在待检测区域范围内设置60个大小为30m×30m方形样地。样地西南角坐标使用差分GPS测定，通过接收JSCORS广域差分信号定位，精度优于0.5米。

本发明基于一个条带的LiDAR数据，在低成本数据获取的前提下，升尺度估算林场级的森林生物量。以江苏南部丘陵地区的亚热带天然次生林为方法实施对象，首先从LiDAR条带数据中提取9个特征变量并与地面估算的生物量结合反演条带内的生物量连续分布信息；然后，在条带反演结果区域内采样并与覆盖整个研究区域的Landsat OLI影像特征变量(5组)结合，从而升尺度估算整个研究区域的生物量。从而在充分挖掘遥感数据特征信息的基础上，通过一个条带的LiDAR数据在降低林场级生物量估算成本的同时，也提升了遥感方法在该尺度上的生物量估算精度。

有益效果：与现有技术相比，本发明通过提取9个LiDAR特征及5组光谱特征变量，从而更为充分地挖掘遥感数据所包含的森林信息特征；只用了一个条带，从而LiDAR数据的获取成本低，且获取了森林冠层的结构信息，故在亚热带地区反演生物量不易“饱和”；制图分辨率提升为30m，且为在北亚热带地区的首次实施。试验结果表明，通过本发明对北亚热带天然次生林的进行生物量估算，取得了较高的精度：模型对地上生物量的估算精度为：R²＝0.69，RMSE＝17.48t·hm^-2，rRMSE＝18％,交叉验证的差值的平均值为-12.92t·hm^-2(差异不显著)。对地下生物量的估算精度为：R²＝0.56，RMSE＝3.33t·hm^-2，rRMSE＝13％,交叉验证的差值的平均值为-0.15t·hm^-2(差异不显著)。

附图说明

图1是试验区和样地分布图；

图2是基于LiDAR-OLI模型的地上生物量模型预测值与重叠区样地实测值的对比散点图及1:1线图；

图3是基于LiDAR-OLI模型的地下生物量模型预测值与重叠区样地实测值的对比散点图及1:1线图；

图4是整个研究区地上生物量信息空间分布图；

图5是整个研究区地下生物量信息空间分布图。

具体实施方式

下面结合具体实例对本发明作进一步描述。

实施例1

一种通过条带LiDAR数据升尺度估算森林生物量的方法，包括以下步骤：

1)试验区概况

研究区位于江苏省常熟市国营虞山林场(120°42′9.4″E,31°40′4.1″N)，属亚热带季风气候，气候温和，年平均降水量1054毫米，面积约1103hm²，其海拔高度为20-261m。虞山林场属于北亚热带次生混交林，主要森林类型为针叶林，阔叶林和混交林，其中主要针叶树种有马尾松(Pinus massoniana)、杉木(Cunninghamia lanceolata)和湿地松(Pinuselliottii)等；主要阔叶树种有麻栎(Quercus acutissima)、枫香(Liquidambarformosana)及部分常绿阔叶树种，如壳斗科植物(Fagaceae)、樟科植物(Lauraceae)和山茶科植物(Theaceae)。试验区及样地分布如图1。

2)遥感数据获取及预处理

采用Riegl LMS-Q680i LiDAR传感器获取的数据(2013年8月17日)。传感器记录了完整的激光脉冲返回波形信息，其时间采样间隔为1ns(即每个采样单元间隔为15cm)。激光发射脉冲所在波段为近红外(波长为1055nm)，脉冲发射频率(PRF)为360kHz，扫描频率(SF)为112线·s^-1，扫描角为±30°；遥感平台的飞行高度为900m，飞行速度为55m·s^-1，旁向重叠度为60％；地面光斑直径约为25cm，平均地面光斑距离约为0.48m。最终获取的数据格式为LAS 1.3(美国摄影测量与遥感协会)标准格式。每个激光返回点内包含三维坐标值(x,y,z)、强度值(intensity)和返回类型(return type)等信息。本实施例使用LiDAR数据其中一个条带范围的数据。

3)汇总地面调查数据

根据虞山林场历史调查数据中的森林类型、年龄和立地指数等指标在试验区范围内设置60个方形样地(大小：30×30m，设置时间：2013年8月)。样地调查过程中，对于胸径大于5cm的树，逐一测定单木的树种、胸径(用围尺测量)、树高和枝下高(利用Vertex IV激光测高器测量)以及冠幅(即两个主方向上的投影距离，用皮尺测量)，对于胸径小于5cm的和枯死木进行计数，但不参与生物量的计算。样地西南角坐标使用差分GPS测定，通过接收JSCORS广域差分信号定位，精度优于0.5米。

根据单木调查数据汇总样地尺度的相关森林参数，包括样地尺度上的单位面积地上、地下生物量(t·hm^-2)。生物量信息通过异速生长方程(见表1)计算单木的生物量，并汇总得到每块样地的单位面积地上生物量(W_A)以及地下生物量(W_B)(见表2)。

表1 研究区各主要树种组分生物量的异速生长模型

注:模型的格式为W＝a(D²H)^b.H＝树高(m),D＝胸径(cm)；a,b为模型参数。

表2 地面估算样地生物量信息汇总

4)特征变量提取

根据LiDAR点云数据计算了9个特征变量：冠层覆盖度变量(c)；5个高度百分位数变量：h₂₅、h₅₀、h₆₀、h₇₅、h₉₅；4个高度变量：最大高度(h_max)、最小高度(h_min)、高度变异系数(h_cv)以及平均高度(h_mean)。9个特征变量的含义及计算公式见表3。

表3 LiDAR特征变量汇总

同时，提取出5组OLI多光谱特征变量：包括6个原始单波段特征、10个波段组合特征、10个信息增强组、18个植被指数特征以及9个纹理特征。其中，OLI多光谱特征变量已由像元尺度经均值计算转化为样地尺度。其中，OLI多光谱特征变量已由像元尺度经均值计算转化为样地尺度。53个光谱特征变量的含义及计算公式见表4。

表4 多光谱特征变量汇总表

注：表中B_i表示高光谱第i波段数据；其中D_i,j表示第i行和第j列位置上的像元亮度值，N表示计算纹理特征时窗口的大小或尺寸。

5)特征变量优化及统计建模

利用样地GPS坐标点位置，基于9景包含LiDAR特征变量信息的影像和53景包含OLI特征变量信息的影像，分别提取样地中心位置的像元值，并借助pearson’s相关系数矩阵分析特征变量之间的相关性，去除相关性高的特征变量。pearson’s相关性分析：

$r=\frac{\underset{i=1}{\overset{n}{\Σ}}(x_i-{\stackrel{\‾}{x}}_i)(y_i-{\stackrel{\‾}{y}}_i)}{\sqrt{\underset{i=1}{\overset{n}{\Σ}}{(x_i-{\stackrel{\‾}{x}}_i)}^2\·\underset{i=1}{\overset{n}{\Σ}}{(y_i-{\stackrel{\‾}{y}}_i)}^2}}---\left(1\right)$

式中：x_i为地面实测的某林分特征；y_i为某特征变量；为x_i的平均值；为y_i的平均值。本方法共建立了2个生物量反演模型(LiDAR模型、LiDAR-OLI模型)。

将地面实测汇总的生物量信息作为因变量，遥感方法提取的特征变量作为自变量，建立多元回归模型。运用逐步进入法(stepwise)和检验决定系数(R²)的变化情况来选择进入模型的合适变量。如果有自变量使统计量F值过小并且T检验达不到显著水平(P值>0.1)，则予以剔除；F值较大且T检验达到显著水平(P值<0.05)则得以进入。为保证自变量之间的低信息重叠，本方法通过控制因子κ(即最大特征根的平方根除以最小特征根所得到的值)进一步筛选模型。κ通过PCA分析计算相关关系矩阵得到，其值低于30则表明自变量之间相关性较低。采用决定系数(R²)、均方根误差(RMSE)和相对均方根误差(rRMSE)评价回归模型的精度，计算公式为：

$R^2=1-\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{(x_i-{\hat{x}}_i)}^2}{\underset{i=1}{\overset{n}{\&Sigma;}}{(x_i-{\stackrel{\&OverBar;}{x}}_i)}^2}---\left(2\right)$

式中：x_i为地面实测的某林分特征；为x_i的平均值；为模型估算的某林分特征；n为样地数量。

$RMSE=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{(x_i-{\hat{x}}_i)}^2}---\left(3\right)$

$rRMSE=\frac{RMSE}{{\stackrel{\&OverBar;}{x}}_i}\&times;100\%---\left(4\right)$

由式(3)可见，rRMSE为RMSE(均方根误差)与(实测值均值)的百分比，作为一个相对量，与评价量本身的数量级水平无关。

LiDAR模型及LiDAR-OLI模型的主要参数及精度评价见表5。LiDAR-OLI模型的生物量模型预测值与重叠区样地实测值的对比散点图及1:1线图见图2和图3。

表5 生物量反演模型主要参数及其精度评价

森林生物量升尺度反演步骤如下：先以重叠区(少量LiDAR数据与OLI影像的重叠区域)样地(31个)LiDAR特征变量为自变量，样地实测生物量为因变量构建LiDAR生物量反演模型(以下简称LiDAR模型)，再利用LiDAR模型计算出的重叠区所有像元的生物量，以其作为因变量，与重叠区的OLI特征变量进行多元线性拟合，建立重叠区生物量反演模型(以下称为LiDAR-OLI模型)。本研究采用交叉验证法评价LiDAR-OLI模型的精度。该方法是在模型参数已经确定的情况下，随机选取所有样地(n)中的一个样地作为验证样地，而其余n-1个样地进行建模，利用拟合得到的模型对随机选取的样地进行验证，循环往复，直至所有样地都验证完毕。方法中，建模组的所有样地参与交叉验证，逐个样地进行检验。LiDAR-OLI模型交叉验证均值差值分析统计结果见表6。最后，将该模型应用于整个试验区，绘制出试验区地上和地下生物量信息空间分布图。整个研究区生物量信息空间分布图见图4和图5。

表6 LiDAR-OLI模型交叉验证均值差值分析统计结果

注：NS表示差异不显著(P>0.05)。

6)方法运算结果

结果表明，通过本方法对北亚热带天然次生林的进行生物量估算，取得了较高的精度：模型对地上生物量的估算精度为：R²＝0.69，RMSE＝17.48t·hm^-2，rRMSE＝18％，交叉验证的差值的平均值为-12.92t·hm^-2(差异不显著)。对地下生物量的估算精度为：R²＝0.56，RMSE＝3.33t·hm^-2，rRMSE＝13％,交叉验证的差值的平均值为-0.15t·hm^-2(差异不显著)。

Claims (6)

1.一种通过条带LiDAR数据升尺度估算森林生物量的方法，其特征在于，包括以下步骤：

1)采用Riegl LMS-Q680i LiDAR传感器获取待检测区域内的LiDAR点云数据，传感器记录完整的激光脉冲返回波形信息，其时间采样间隔为1ns；

2)在待检测区域范围内设置多个方形样地，样地调查过程中，对于胸径大于5cm的树，逐一测定单木的树种、胸径、树高和枝下高以及冠幅；根据单木调查数据汇总样地尺度的相关森林参数，包括样地尺度上的单位面积地上生物量、地下生物量；生物量信息通过异速生长方程计算单木的生物量，并汇总得到每块样地单位面积地上生物量W_A以及地下生物量W_B；

3)根据LiDAR点云数据计算9个特征变量：冠层覆盖度变量c；5个高度百分位数变量：h₂₅、h₅₀、h₆₀、h₇₅、h₉₅；4个高度变量：最大高度h_max、最小高度h_min、高度变异系数h_cv以及平均高度h_mean；同时，提取出5组OLI多光谱特征变量：包括6个原始单波段特征、10个波段组合特征、10个信息增强组、18个植被指数特征以及9个纹理特征；其中，OLI多光谱特征变量已由像元尺度经均值计算转化为样地尺度；

4)利用样地GPS坐标点位置，基于9景包含LiDAR特征变量信息的影像和53景包含OLI特征变量信息的影像，分别提取样地中心位置的像元值，并借助pearson’s相关系数矩阵分析特征变量之间的相关性，去除相关性高的特征变量；pearson’s相关性分析：

$r=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}(x_i-{\stackrel{\&OverBar;}{x}}_i)(y_i-{\stackrel{\&OverBar;}{y}}_i)}{\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{(x_i-{\stackrel{\&OverBar;}{x}}_i)}^2\&CenterDot;\underset{i=1}{\overset{n}{\&Sigma;}}{(y_i-{\stackrel{\&OverBar;}{y}}_i)}^2}}---\left(1\right)$

式中：x_i为地面实测的某林分特征；y_i为某特征变量；为x_i的平均值；为y_i的平均值；

将地面实测汇总的生物量信息作为因变量，遥感方法提取的特征变量作为自变量，建立多元回归模型；运用逐步进入法和检验决定系数R²的变化情况来选择进入模型的合适变量；如果有自变量使统计量F值过小并且T检验达不到显著水平，则予以剔除；F值较大且T检验达到显著水平则得以进入；通过控制因子κ进一步筛选模型；κ通过PCA分析计算相关关系矩阵得到，其值低于30则表明自变量之间相关性较低；采用决定系数R²、均方根误差RMSE和相对均方根误差rRMSE评价回归模型的精度，计算公式为：

$R^2=1-\frac{\underset{i=1}{\overset{n}{\&Sigma;}}{(x_i-{\hat{x}}_i)}^2}{\underset{i=1}{\overset{n}{\&Sigma;}}{(x_i-{\stackrel{\&OverBar;}{x}}_i)}^2}---\left(2\right)$

式中：x_i为地面实测的某林分特征；为x_i的平均值；为模型估算的某林分特征；n为样地数量；

$RMSE=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\&Sigma;}}{(x_i-{\hat{x}}_i)}^2}---\left(3\right)$

$rRMSE=\frac{RMSE}{{\stackrel{\&OverBar;}{x}}_i}\&times;100\%---\left(4\right)$

2.根据权利要求1所述的通过条带LiDAR数据升尺度估算森林生物量的方法，其特征在于，森林生物量升尺度反演步骤如下：先以重叠区样地LiDAR特征变量为自变量，样地实测生物量为因变量构建LiDAR生物量反演模型，再利用LiDAR模型计算出的重叠区所有像元的生物量，以其作为因变量，与重叠区的OLI特征变量进行多元线性拟合，建立重叠区生物量反演模型。

3.根据权利要求1所述的通过条带LiDAR数据升尺度估算森林生物量的方法，其特征在于，Riegl LMS-Q680i LiDAR传感器的激光发射脉冲所在波段为近红外，波长为1055nm，脉冲发射频率为360kHz，扫描频率为112线·s^-1，扫描角为±30°；遥感平台的飞行高度为900m，飞行速度为55m·s^-1，旁向重叠度为60％；地面光斑直径约为25cm，平均地面光斑距离约为0.48m。

4.根据权利要求1所述的通过条带LiDAR数据升尺度估算森林生物量的方法，其特征在于，Riegl LMS-Q680i LiDAR传感器最终获取的数据格式为LAS 1.3标准格式，每个激光返回点内包含三维坐标值、强度值和返回类型信息。

5.根据权利要求1所述的通过条带LiDAR数据升尺度估算森林生物量的方法，其特征在于，在待检测区域范围内设置60个大小为30m×30m方形样地。

6.根据权利要求1所述的通过条带LiDAR数据升尺度估算森林生物量的方法，其特征在于，样地西南角坐标使用差分GPS测定，通过接收JSCORS广域差分信号定位，精度优于0.5米。