一种全波形激光雷达和高光谱数据联合反演森林结构参数的方法与流程

本发明属于森林资源监测和环境因子调查等技术领域，涉及一种全波形激光雷达和高光谱数据联合反演森林结构参数的方法。

背景技术

精确的森林结构参数提取，对于森林资源监测和环境因子调查具有重要意义。同时，这些信息也可以用于掌握森林植物与环境的关系，进一步掌握森林的生长、发育规律，对于森林的可持续经营管理、生态模型构建及全球碳循环研究具有重要意义。常规的森林结构参数提取主要依赖于野外调查及航片或卫片判读等，其精度往往不高，且难于在“面”上实用化推广。

近年来，基于激光雷达数据与高光谱数据进行森林结构参数提取的研究为：latifi等2012年在《remotesensingofenvironment》第121卷上发表的“foreststructuremodelingwithcombinedairbornehyperspectralandlidardata”，该研究采用航空飞机获取的高光谱和激光雷达数据，结合光谱特征、高度/强度特征和遗传方法估测了森林结构参数。kandare等2017年在《remotesensing》第9卷上发表的“predictionofspecies-specificvolumeusingdifferentinventoryapproachesbyfusingairbomelaserscanningandhyperspectraldata”，该研究采用航空飞机获取激光雷达和高光谱影像，并基于自动和半自动方法提取了树冠信息，并据此估测了复杂高山森林的森林结构参数。sankey等2017年在《remotesensingofenvironment》第195卷上发表的“uavlidarandhyperspectralfusionforforestmonitoringinthesouthwesternusa”，该研究采用无人机获取的激光雷达和高光谱影像，结合冠层结构特征和光谱特征估测了森林结构参数。然而，以上方法都是基于激光雷达点云数据，并没有利用能够更好记录冠层结构特性的全波形激光雷达数据。同时，更未见全面深入提取森林冠层激光雷达特征(波形特征和点云特征)、光谱特征及计算森林结构参数提取的方法。

技术实现要素：

发明目的：针对现有技术中存在的不足，本发明的目的是提供一种全波形激光雷达和高光谱数据联合反演森林结构参数的方法，通过联合全波形激光雷达数据和高光谱数据的方法，可以从波形和光谱数据中提取出更好描述森林结构的特征，将有助于提高森林结构参数的反演精度，并有效抑制森林覆盖度高、生物量高林分结构参数反演的“饱和”问题。

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

一种全波形激光雷达和高光谱数据联合反演森林结构参数的方法，先对机载全波形激光雷达数据进行去噪、平滑、强度校正、滤波，插值生成数字地形模型，点云及波形数据高度归一化处理；对高光谱影像进行辐射定标、大气校正、几何校正预处理；然后，分别基于归一化点云和波形数据、预处理高光谱数据分别提取特征变量；最后，结合地面实测数据和提取的特征变量分别构建多元回归模型以预测各森林结构参数。

所述的全波形激光雷达和高光谱数据联合反演森林结构参数的方法，步骤如下：

1)在地面设置样地，借助航空飞机采集全波形激光雷达数据和高光谱数据，并在样地中记录树种并计数，同时测量每木的胸径和树高；

2)去除激光雷达波形数据的背景噪声并对波形数据进行平滑，基于高斯分解方法和levenberg-marquardt方法对波形数据进行拟合以提取点云数据；同时基于体元的方法对波形数据进行校正以降低扫描角的影响；基于滤波方法去除非地面点，然后通过计算每个像元内激光点高度的平均值，生成数字地形模型；

3)借助传感器辐射定标数据对原始高光谱数据进行辐射定标，并利用经验线性模型结合地面实测标靶光谱数据进行大气校正；同时，利用地面实测控制点对高光谱影像进行几何精校正；

4)提取激光雷达特征变量和光谱特征变量；

5)通过相关性分析筛选特征变量，首先筛选特征变量之间相关性低于0.6的特征变量，然后进一步筛选特征变量与各林分结构参数相关性高于0.6的特征变量；

6)将地面实测森林结构参数作为因变量，各特征变量作为自变量，建立多元回归模型；运用逐步回归法选择进入模型的变量，即在预先给定的f水平下进行显著性检验，如果t检验达不到显著水平(p＞0.1)，则予以剔除；t检验达到显著水平(p＜0.05)则予以进入。

步骤1)中，蓄积量根据一元材积公式结合实测胸径进行估算，地上生物量通过异速生长方程结合胸径和树高进行计算。

步骤2)中，数字地形模型的空间分辨率为0.5m。

步骤4)中，激光雷达特征变量包括：冠层高度分布百分位数，冠层点云分布平均高度以上的覆盖度，冠层点云分布的变异系数，点云数量在各百分数高度以上的点占所有点云的百分比，冠层体积与剖面特征变量，冠层各结构类别体积占比，波形几何和辐射特性。

步骤4)中，剖面特征变量为weibull函数对冠层高度分布剖面进行拟合得到2个剖面特征量α，β。

步骤4)中，冠层各结构类别包括开放层，透光层，低光层和封闭层四个冠层结构类别，

步骤4)中，光谱特征变量包括：原始光谱波段，植被指数，主成分分析前十主成分，独立成分分析前十变量，最小噪声分离前十变量。

步骤6)中，通过主成分分析计算相关关系矩阵得到控制因子k，k小于30则模型进一步被选中。

所述的全波形激光雷达和高光谱数据联合反演森林结构参数的方法，采用决定系数r²、均方根误差rmse和相对均方根误差rrmse评价回归模型拟合的效果及估测精度，各自的计算公式如下：

式中，xi为某森林结构参数实测值；为某森林结构参数实测平均值；为某森林结构参数的模型估测值；n为样地的数量；i为某一个样地。

有益效果：与现有技术相比，本发明具有以下优势：

1)全波形激光雷达数据能够获得准确的森林冠层三维结构特征信息，高光谱数据具有光谱分辨率高的特点，能够记录森林冠层结构特性对不同波长电磁波的反射特性的影响。通过联合全波形激光雷达数据和高光谱数据的方法，可以从波形和光谱数据中提取出更好描述森林结构的特征，将有助于提高森林结构参数的反演精度，并有效抑制森林覆盖度高、生物量高林分结构参数反演的“饱和”问题。

2)本方法通过联合激光雷达和高光谱两类数据源，综合所提取的特征变量建立森林结构参数估算模型，激光雷达与高光谱数据分别从空间三维和光谱维的角度记录了森林结构信息，其数据之间互为补充，有效增强了森林结构参数反演的能力和精度；

3)本方法将记录森林冠层全部反射信号的全波形激光雷达技术与高光谱数据相结合，可增强数据的冠层空间结构刻画能力，进而提高森林结构参数的估算精度；

4)本方法全面深入地提取了多组森林冠层激光雷达特征(波形特征和点云特征)、光谱特征，并进行了特征变量优选，从而高质量地提取了森林结构参数。同时，该发明不仅利于挖掘特征变量与森林结构参量之间的联系，从机理上对特征变量与森林结构参数之间的关系进行解释，也易于进行方法移植(即在不同地区的不同森林类型中也可以进行应用)；

5)验证结果表明，通过本发明对亚热带天然次生林的森林结构参数进行提取，与使用其他相近遥感方法进行林分结构参数相比，其相对均方根误差提升了5％以上。

附图说明

图1是全波形激光雷达数据三维显示及样地照片、球面镜照片与合成波形、光谱反射率曲线图；

图2是全波形激光雷达和高光谱数据联合反演森林结构参数原理图(含数据获取、数据融合与特征提取)。

具体实施方式

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

以下实施例的试验区位于江苏省常熟市国营虞山林场(120.70°e，31.67°n)，面积约1422hm²，高程变化范围为2-261m。试验区所处区域为亚热带季风气候，年降水量1062.5mm。其森林类型属于亚热带次生混交林，可细分为针叶林、阔叶林和混交林。其中主要针叶和阔叶落叶树种包含马尾松(pinusmassoniana)、麻栎(quercusacutissima)、枫香树(liquidambarformosana)和板栗(castaneamollissima)等，同时伴生部分常绿阔叶树种。

实施例1

一种全波形激光雷达和高光谱数据联合反演森林结构参数的方法，如图1所示，步骤如下：

1)在试验区范围内设置67个方形样地(30×30m²)，借助航空飞机采集全波形激光雷达数据和高光谱影像数据，如图2所示。样地中心点坐标使用gps(trimblegeoxh6000)测定，gps通过接收广域差分信号定位，精度优于0.5m。并在样地中记录树种并计数，同时测量每木的胸径和树高。蓄积量根据一元材积公式结合实测胸径进行估算，地上生物量通过异速生长方程结合胸径和树高进行计算。根据单木调查数据汇总到样地尺度的地上生物量，结果见表1。

表1样地实测林分特征信息汇总表

2)数据预处理时，首先采用高通滤波方法去除激光雷达波形数据的背景噪声，并对波形数据进行平滑，基于高斯分解方法和levenberg-marquardt方法对波形数据进行拟合以提取点云数据。同时基于体元的方法对波形数据进行校正以降低扫描角的影响。基于滤波方法去除非地面点，然后通过计算每个像元内激光点高度的平均值，生成数字地形模型(dem)(空间分辨率为0.5m)。

3)借助传感器辐射定标数据对原始高光谱数据进行辐射定标，并利用经验线性模型结合地面实测标靶光谱数据进行大气校正(校正方式为首先采集多个标靶的影像光谱值与实测光谱值，然后建立这些标靶影像光谱值与实测光谱值的线性模型，最后利用该模型对整幅影像进行校正)。同时，利用地面实测控制点对高光谱影像进行几何精校正(校正方式为首先采集多个地面控制点的影像坐标值与实测投影坐标系值，然后建立这些地面控制点的影像坐标值与实测投影坐标系值的转换模型，最后利用该模型对整幅影像进行几何位置校正，并对像元值进行重采样)。

4)特征变量提取时，共提取两组特征变量，即激光雷达特征变量(含波形及点云特征)和光谱特征变量。

激光雷达特征变量包括：冠层高度分布百分位数(h25，h50，h75，h95)；冠层点云分布平均高度以上的覆盖度(ccmean)；冠层点云分布的变异系数(hcv)；点云数量在各百分数高度(30th，50th，70th，90th，即d3，d5，d7，d9)以上的点占所有点云的百分比；冠层体积与剖面特征变量：weibull函数对冠层高度分布剖面进行拟合得到2个剖面特征量α，β(即weibullα和weibullβ)；冠层各结构类别体积占比，包括开放层，透光层，低光层和封闭层四个冠层结构类别，每个冠层结构类别的体积所占百分比(即opengap，oligophotic，euphotic，closedgap)；波形几何和辐射特性(home，wd，vdr，np，rough，fs，rwe，int，fwhm)。

光谱特征变量包括：原始光谱波段(band1-band64)；植被指数(sr，ndvi，evi，gndvi，savi，arvi，rvsi，psri，vog1，vog2，rgri，pri，prr，wbi，cri1，cri2，ari1，ari2)；主成分分析前十主成分(pc1-pc10)；独立成分分析前十变量(ic1-ic10)；最小噪声分离前十变量(mnf1-mnf10)。

5)通过相关性分析(即求取各特征变量和各林分结构参数之间的pearson’s相关系数)筛选特征变量。即首先筛选特征变量之间相关性低于0.6的特征变量，然后进一步筛选特征变量与各林分结构参数相关性高于0.6的特征变量。

6)将地面实测森林结构参数作为因变量，各特征变量作为自变量，建立多元回归模型。运用逐步回归法选择进入模型的变量，即在预先给定的f水平下进行显著性检验，如果t检验达不到显著水平(p＞0.1)，则予以剔除；t检验达到显著水平(p＜0.05)则予以进入。为了降低自变量之间的相关性，本方法通过主成分分析计算相关关系矩阵得到控制因子k(即最大特征根的平方根与最小特征根的比值)，k小于30则模型进一步被选中。本方法采用决定系数(r²)、均方根误差(rmse)和相对均方根误差(rrmse)评价回归模型拟合的效果及估测精度：

式中，xi为某森林结构参数实测值；为某森林结构参数实测平均值；为某森林结构参数的模型估测值；n为样地的数量；i为某一个样地。

各林分结构参数估测模型及模型预测精度见表2，可见该方法对于各林分结构参量估测精度较高，因此通过全波形激光雷达和高光谱数据联合反演森林结构参数的方法是可靠且可行的。

表2各林分结构参数估测模型及模型预测精度

注：h50，h75，d1，e，o分别为冠层为25％，75％的高度分布百分位数，点云数量在百分数高度10th上的点占所有点云的百分比，透光层和低光层体积所占百分比；homeμ，vdrμ，npμ，rweμ，andwdσ分别为波形质心到地面的平均距离，波形起点到波形质心的平均距离除以波形起点到地面的距离，波形脉冲上被探测点的平均个数，平均接收总能量，波形起点到地面点距离标准差；ndvi，rvsi，arvi，cri1andpca2分别为归一化植被指数，红边胁迫指数，大气调节植被指数，胡萝卜素反射指数，主成分分析第二主成分。