基于多回波机载激光扫描数据的城市树木碳含量估计方法

摘要

本发明公开了一种基于多回波机载激光扫描数据的城市树木含碳量估计方法，包括：S1、使用多光谱机载激光扫描仪对城市树木进行扫描，获得多个航带的多回波ALS点云数据；S2、对多个航带的多回波ALS点云数据进行预处理，获得nDSM数据及各激光通道的强度图像；S3、对nDSM数据进行分类，提取其中的树木类别数据，获得树冠高度模型CHM；S4、由树冠高度模型CHM估计树木的形状参数,所述树木的形状参数包括树冠高度与树冠直径；S5、先由经验模型构建ALS‑DBH回归模型，再通过现场选取若干树木样本进行验证；S6、根据所有物种的DBH和树冠高度来估算碳含量。本发明基于ALS的树木形状参数和异速生长模型，可对城市树木碳含量进行快速无损的估计。

说明书

技术领域

本发明涉及一种基于多回波机载激光扫描数据的城市树木碳含量估计方法。

背景技术

随着人口的增长，城镇化成了不可阻挡的趋势。当城市扩大的同时，生态和地面附着物也随之改变。地面附着物从自然环境变成各种城市建筑物，比如：不透水的公路、高楼大厦以及公园等等。这些转变带来了些紧迫的问题，比如城市热岛效应。城市扩大和经济活动还增加能量消耗并导致温室气体释放。而城市里的植被对城市气候的改善有着重要的作用，例如，树荫和风屏蔽效应可以减弱直接太阳热量和空气渗透到房屋。同时城市植被存储着大量的碳，能很大程度上减少城市热岛效应。因此，估计和监测城市树木含碳量就尤为重要了。

估计城市树木碳含量需要大量的现场数据，如果采用人工测量，将需要巨大的劳动力和时间。此外，估计树木干重往往又是存在破坏性的，因此，有必要提出一种快捷、高效的方法对城市树木碳含量进行估计。

发明内容

本发明的目的在于克服上述现有技术的不足，提供一种基于多回波机载激光扫描数据的城市树木碳含量估计方法。

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

基于多回波机载激光扫描数据的城市树木含碳量估计方法，包括以下步骤：

S1、使用包括绿激光、近红外激光及短波红外激光的机载激光扫描仪对城市树木进行扫描，获得多个航带的多回波ALS点云数据；

S2、对多个航带的多回波ALS点云数据进行预处理，获得nDSM数据及各激光通道的强度图像；

S3、植被分类：对nDSM数据进行分类，提取其中的树木类别数据，获得树冠高度模型CHM；

S4、由树冠高度模型CHM估计树木的形状参数,所述树木的形状参数包括树冠高度与树冠直径；

S5、ALS-DBH回归模型的选择：先由经验模型构建ALS-DBH回归模型，再通过现场选取若干树木样本获得偶数组采样组，将其中一半采样组作为模型组用于模型开发，另一半采样组作为验证组，将所有的组合情况进行模型开发，并通过交叉验证过程拟合，收集每个组合产生的参数；通过测定系数R²检查回归模型的预测能力和性能，并且通过预测参数的RMSE来检查预测的准确度，比较所有组合情况中的R²和RMSE，选择在模型拟合中具有高R²和在验证中具有低RMSE的一个ALS-DBH回归模型来预测DBH；

S6、碳含量的预测：

根据所有物种的DBH和树冠高度来计算总生物量，存储在树中的碳存储量估计为总生物量的一半，将预测的碳储存与通过现场测量的DBH和高度估计的碳储存进行比较，并且通过RMSE和R²进行评估。

优选地，步骤S1中所述绿激光、近红外激光及短波红外激光的波长分别为532nm、1064nm及1550nm。

优选地，步骤S2具体包括：

S21、分别对每个航带的多回波ALS点云数据去噪；

首先，采用强度滤波器将强度在3σ以外的点去除；

其次，采用统计滤波器统计某一点与其邻近若干点之间的平均距离，若该平均距离超过全局平均距离与一个标准偏差的和，则去除该点；

S22、对去噪后的每个航带的多回波ALS点云数据的强度进行幅度归一化；

式中，I′是幅度归一化强度，I_raw是ALS点云数据原始强度，a是由目标几何结构确定的指数，R是飞机到目标的距离，R_ref为平均参考距离；

S23、生成各通道的强度图像：将归一化后的多个航带的多回波ALS点云数据融合在一起，每个航带带有三个波长的点云数据，对融合的点云数据进行栅格化处理，则各格子中的像素值为格子中所有点的强度平均值；对于内部没有点的格子，由其邻近格子的像素值进行线性插值得到；

S24、地面点和非地面点的分类；

将归一化后的多个航带的多回波ALS点云数据融合在一起，用基于斜率的滤波算法对融合后的点云数据进行地面点和非地面点分类；

S25、将非地面点数据进行栅格化处理，获得nDSM数据；

将非地面点云数据分到边长为L的栅格中，以格子内所有点的最大高度值为格子的像素值，对于内部没有点的格子的像素值由其邻近格子的像素值进行线性插值得到。

优选地，步骤S24获得的nDSM数据中，若某一格子的像素值与其八邻域间有2m以上的高度差，则用其八邻域的高度均值作为该格子的像素值。

优选地，步骤S3具体为：

S31、以绿激光强度、近红外激光强度、短波红外激光强度、nDSM、pNDWI、pNDVI为输入，在输入数据选择训练样本，并通过谷歌地球的航拍照片获取训练样本所属类别，将训练样本用于SVM分类器训练，使用训练好的SVM分类器对nDSM数据进行分类，提取出其中的树木类别；

其中，

C_Green，C_NIR和C_SWIR分别为绿激光强度、近红外激光强度、短波红外激光强度；

S32、采用3*3的众数滤波器对提取出的树木类别进行若干次的滤波；

S33、经验地测试滤波次数，产生最佳树冠，生成树冠高度模型CHM，所述最佳树冠指可以移除孤立的像素并且保留优势树的树冠形状。

优选地，步骤S4具体为：

S41、获得局部最大像素集合：对CHM采用3*3的局部最大值滤波器进行滤波，将像素值比其八邻域都高的像素点定义为局部最大值，获得局部最大像素集合；

S42、树梢的确定：使用ALS强度数据对局部最大像素集合进行进一步抽取，将在ALS强度图像中也是局部最大的像素点记为树梢点；

S43、将检测到的树梢作为标记，采用区域生长分割算法将CHM分割成独立树冠，分割后每个树梢点都有一个封闭的切块，用绝对精度表示树木分离效果，则有：

其中n_1，1是检测到的树冠切块数量，n_total是实际树冠数量；树高定义为在每个树冠切块中局部最大值的均值，树冠直径则定义为经过树冠最高点中心的最大树冠直径与其垂直方向直径的均值。

优选地，步骤S5具体为：

S51、ALS-DBH回归模型构建：

以树高度和树冠直径作为独立预测因子，根据经验模型构建ALS-DBH回归模型以预测DBH，则有，

DBH_Field＝a·CD_ALS+b·H_ALS+c；

式中，CD是树冠直径，以米为单位；H则是树冠高度，以米为单位；a,b,c是从经验模型导出的参数；

其中CD是树冠直径，以米为单位；H则是树冠高度，以米为单位；a,b,c指的是从经验模型导出的参数，及分别为DBH、CD及H的均值；

S52、ALS-DBH回归模型选择：

现场选取若干树木样本，将树木样本分成a组，a为偶数，选择其中a/2组数用于模型开发，a/2组用作验证，则存在种组合情况，将种组合情况通过交叉验证过程拟合，收集每个组合产生的参数，通过测定系数R²检查回归模型的预测能力和性能，并且通过预测参数的RMSE来检查预测的准确度；

其中y是树木形状参数，y_fitted是从ALS数据导出的树木形状参数，是田间实测树木形状参数，比较组ALS-DBH回归模型的R²和RMSE，选择在模型拟合中具有高R²和在验证中具有低RMSE的一个回归模型来预测DBH。

优选地，步骤S6中，总生物量为树叶、树枝、木材及树皮中的生物量的总和，则有，

y_total＝y_foliage+y_branches+y_wood+y_bark；

则，

其中y是树木组分干重，单位为Kg；D是DBH，单位为cm；β是模型参数，每种树木组分有一个相应的β；e则是误差补偿。

优选地，还包括步骤S7；

S7、创建城市碳存储地图：

S71、将ALS数据导出的树冠部分以及其中存储的碳含量在GIS中转换成矢量数据；

S72、将城市土地分为若干种覆盖类型，将每种土地覆盖类型中的碳存储总量除以该土地覆盖类型的总面积获得每种土地覆盖类型的单位面积碳存储量；

S73、根据土地覆盖类型对城市土地进行区域划分，将每个区域的总面积乘以其对应的土地覆盖类型的单位面积碳存储量，估计出每一区域的碳含量，最后创建城市碳存储地图。

采用上述技术方案后，本发明与背景技术相比，具有如下优点：

1、本发明采用的数据是通过多光谱ALS获得，根据不同附着物对不同波长的反射率不同，可以有效的分离出植被，植被分类准确率达到90.23％，比单光谱的准确率高了11％；

2、对CHM采用局部最大滤波后得到树梢，再次通过最大强度图像滤波，最终留下树梢，可以有效的提高树梢检测准确度；

3、通过ALS数据得出树木形状参数(树高、树冠直径和树干胸径)建立碳含量模型，能够无破坏的、快捷的计算出城市树木碳含量，并绘制成城市碳含量地图。

附图说明

图1为本发明总体流程框图；

图2为本发明细节流程框图；

图3为地面附着物分类示意图；

图4为本发明步骤S4中各步骤的结果示意图，其中，4(a)为树冠高度模型局部最大值滤波结果示意图；图4(b)为最大强度图像局部最大值滤波结果示意图；图4(c)为区域生长分割结果；

图5为本发明所创建的城市树木碳含量地图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

实施例

如图1及图2所示，基于多回波机载激光扫描数据的城市树木含碳量估计方法，包括以下步骤：

S1、获得多个航带的多回波ALS点云数据；

S2、获得nDSM数据及各激光通道的强度图像；

S3、植被分类，获得树冠高度模型CHM；

S4、获得树木的树冠高度与树冠直径；

S5、ALS-DBH回归模型的选择；

S6、碳含量的预测；

S7创建城市碳存储地图。

步骤S1中所述绿激光、近红外激光及短波红外激光的波长分别为532nm、1064nm及1550nm。

步骤S2具体包括：

S21、分别对每个航带的多回波ALS点云数据去噪：

首先，采用强度滤波器将强度在3σ以外的点去除，以去除具有高强度的镜面反射点；

其次，采用统计滤波器统计某一点与其邻近若干点之间的平均距离，若该平均距离超过全局平均距离与一个标准偏差的和，则去除该点。

S22、对去噪后的每个航带的多回波ALS点云数据的强度进行幅度归一化，归一化方程为：

式中，I′是幅度归一化强度，I_raw是ALS点云数据原始强度，a是由目标几何结构确定的指数，R是飞机到目标的距离，R_ref为平均参考距离。一般，对于有大表面的物体比如阔叶树，指数a取值为2最好，而对于点状的物体，比如针叶树，指数a取值为3最好，本发明中a的取值为2。

S23、生成各通道的强度图像：将归一化后的多个航带的多回波ALS点云数据融合在一起，每个航带带有三个波长的点云数据，然后将每个点云栅格化成具有1m地面分辨率的强度图像。这里像素的大小根据数据集的点间距进行选择，选择像素大小接近点距离，可以使得多数的像素都有至少一个点且点的垂直分布会显著突出。这样就将点分到了边长为1m的栅格里，该格子里的像素值就为格子内所有点的强度平均值；对于内部没有点的格子，有其邻近格子的像素值进行线性插值得到。分别对波长为532nm、1064nm和1550nm的通道生成各自的强度图像。

S24、地面点和非地面点的分类；

将归一化后的多个航带的多回波ALS点云数据融合在一起。用基于斜率的滤波算法对融合后的点云数据进行地面点和非地面点分类。

S25、将非地面点数据进行栅格化处理，获得nDSM数据：

将所有的地面点分到边长为L的栅格里，本实施例L＝1，并用线性插值出格内没点的格子，得到栅格化数据；

非地面点也用相似的方法进行栅格化得到nDSM，不同的是，格子上的像素值改为格内所有点的最大高度，这样就能最大的体现树高；

上面得到的nDSM中含有凹洞或者像素值为负数的点，这些凹洞往往是由于栅格化、地面点滤除和插值等操作产生。如果某一格子的像素值与其八邻域间有2m以上的高度差，就用其八邻域高度均值填充该格子。

步骤S3具体为：

S31、以绿激光强度、近红外激光强度、短波红外激光强度、nDSM、pNDWI、pNDVI为输入，在输入数据中截取训练样本，并通过谷歌地球的航拍照片获取样本所属类别，将训练样本用于SVM分类器训练，使用训练好的SVM分类器对nDSM数据进行分类，提取出其中的树木类别；

其中，

C_Green，G_NIR和C_SWIR分别为绿激光强度、近红外激光强度、短波红外激光强度，pNDWI及pNDVI指标可以很好地区分出人造物体。

由于研究区域是简单的住宅区，如图3所示，可以将目标分为水域、房子、公路、草地、树木和开放区域六类地面附着物，但由于水和水中的有机物使得对红光的反射强度变得不规则，因此将水域从点云中去除，只对五种地面附着物进行分类。在所有输入数据中截取训练样本，并通过谷歌地球的航拍照片获取样本所属类别。

对每种地面附着物至少选取600像素作为训练样本，而对于植被则有1000像素，占据了最大比例。将这些像素点用于SVM分类器训练。选择SVM分类器用于执行具有多光谱ALS衍生数据的土地覆盖分类是因为：由于SVM在典型的ALS相关分类研究中的普及，使得该研究的分类结果可以和以前的研究相比较。在遥感数据分类中常常采用的核为：

并将核宽Υ取值为0.167，正则化参数C取值100。

为了估计分类结果的质量，为了评估分类结果的质量，通过研究区域的分层随机抽样选择总共522个参考像素，并手动标记为土地覆盖类型。对于每个类别，至少标记60个参考像素，并且对于占据研究区域的大部分的植被覆盖物，包括超过100个参考像素，然后根据标记错误和遗漏错误生成错误矩阵，用于估计分类结果质量。如图2所示的是

S32、采用3*3的众数滤波器对提取出的树木类别进行若干次的滤波，能有效地减少椒盐噪声；

S33、经验的测试滤波次数，产生最佳树冠，生成树冠高度模型CHM，所述最佳树冠指可以移除孤立的像素并且保留优势树的树冠形状。

通过这种方法，可以填充树冠内的凹坑；同时过滤的树类点用于掩蔽nDSM以生成孤立树的nDSM。孤立树的nDSM也被称为树冠高度模型(CHM)，即通过树冠在自顶向下视图中显示树位置并以像素存储高度值的模型。

步骤S4具体为：

S41、获得局部最大像素集合：对CHM采用3*3的局部最大值滤波器进行滤波，将像素值比其八邻域都高的像素点定义为局部最大值，获得局部最大像素集合；如图4(a)所示的是树冠高度模型局部最大值滤波结果。

S42、树梢的确定：为了消除与局部最大值滤波相关的标记错误并检测真实的树梢像素点，使用ALS强度数据对局部最大像素集合进行进一步抽取，将在ALS强度图像中也是局部最大的像素点记为树梢点。

ALS强度值不仅取决于对象的反射率，而且取决于传感器和对象之间的距离。树梢像素点既有较高的强度值也有高的像素值。基于这点，对至少有15个像素点组成的局部最大像素集合再进一步抽取，只有那些在最大强度图像中也是局部最大的像素点才最终保留为树梢点。如图4(b)所示的是最大强度图像局部最大值滤波结果。

S43、将检测到的树梢作为标记，采用区域生长分割算法将CHM分割成独立树冠，分割后每个树梢点都有一个封闭的切块，用绝对精度表示树木分离效果则有：

其中n_1，1是检测到的树冠切块数量，n_total是实际树冠数量；树高定义为在每个树冠切块中局部最大值的均值，树冠直径则定义为经过树冠最高点中心的最大树冠直径与其垂直方向直径的均值。如图4(c)所示的是区域生长分割结果。

步骤S5具体为：

S51、ALS-DBH回归模型构建：

通过ALS导出树高度和树冠直径作为独立预测因子，根据经验模型构建ALS-DBH回归模型以预测DBH，则有，

DBH_Field＝a·CD_ALS+b·H_ALS+c；

式中，CD是树冠直径，以米为单位；H则是树冠高度，以米为单位；a,b,c是从经验模型导出的参数；

其中CD是树冠直径，以米为单位；H则是树冠高度，以米为单位；a,b,c指的是从经验模型导出的参数，及分别为DBH、CD及H的均值；

S52、ALS-DBH回归模型选择：

现场选取若干树木样本(本实施例为40棵)，将树木样本分成a组，a为偶数，选择其中a/2组数用于模型开发，a/2组用作验证，则存在种组合情况。本实施例a＝4，则存在6种组合情况，将6种情况均以0.05显著性水平开发模型，并通过交叉验证过程拟合，收集每个组合产生的参数，通过测定系数R²检查回归模型的预测能力和性能，并且通过预测参数的RMSE来检查预测的准确度。

其中，

式中，y是树木形状参数，y_fitted是从ALS数据导出的树木形状参数，是田间实测树木形状参数，比较组ALS-DBH回归模型的R²和RMSE，选择在模型拟合中具有高R²和在验证中具有低RMSE的一个回归模型来预测DBH。

步骤S6根据所有物种的DBH和树冠高度来计算总生物量，存储在树中的碳存储量估计为总生物量的一半，将预测的碳储存与通过现场测量的DBH和高度估计的碳储存进行比较，并且通过RMSE和R²进行评估。

选择基于所有物种的DBH和高度的方程组来计算生物量，因为在本实施例中没有属性或物种信息。因而，总生物量为树叶、树枝、木材及树皮中的生物量的总和，则有，

y_total＝y_foliage+y_branches+y_wood+y_bark；

则，

其中y是树木组分干重，单位为Kg；D是DBH，单位为cm；β是模型参数，每种树木组分有一个相应的β；e则是误差补偿。

S7、创建城市碳存储地图：

S71、将ALS数据导出的树冠部分以及其中存储的碳含量在GIS中转换成矢量数据；

S72、将城市土地分为若干种覆盖类型，将每种土地覆盖类型中的碳存储总量除以该土地覆盖类型的总面积获得每种土地覆盖类型的单位面积碳存储量；

S73、根据土地覆盖类型对城市土地进行区域划分，将每个区域的总面积乘以其对应的土地覆盖类型的单位面积碳存储量，估计出每一区域的碳含量，最后创建城市碳存储地图，如图5所示的是城市树木碳含量地图。

采用本发明所述的方法，在地面点和非地面点分类中正确分类的比例为96.4％，在地面附着物分类中，树木正确分类达到90.23％，而经典ALS数据只有79％的正确率。在树木形状参数测量中，树高的均方根误差RMSE为1.21m，树冠直径的均方根误差RMSE则为1.47m。ALS-DBH模型预测中DBH的相关系数R²为0.80，均方根误差RMSE为4.6cm。碳含量的RMSE则为142Kg。这结果表明，基于ALS的树木形状参数和异速生长模型可以产生一致的性能和准确的估计，充分揭示了多光谱ALS数据在估计个体树水平的碳储存和在城市环境中绘制植被的潜力。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。