基于高分遥感图像和地基三维激光雷达数据的城市三维绿量快速测算*

刘炜 赵尔平 雒伟群 刘爱华

（西藏民族大学 西藏光信息处理与可视化技术重点实验室，陕西 咸阳 712082）

0 引言

三维绿量(Living Vegetation Volume,LVV)是指地表生长植物的茎、叶所占有的空间体积，单位为立方米(m3)。相较于以平方米(m2)为量纲的二维绿量指数，三维绿量从立体空间的角度评定植被空间结构差异和生态功能水平高低，因而能够更准确的反映城市绿色景观格局变化及其在缓解热岛效应、净化大气、释氧固碳、水土保持、维持生物多样性等方面功能作用［1］。

近几年，我国学者利用遥感技术或三维激光雷达扫描技术，在北京、上海、沈阳、武汉、福州等城市针对城市三维绿量测算开展了大量研究工作。其中，冯仲科［2-3］、刘常富［4-5］、姚崇怀［6］等的研究工作具有代表性。在国外，基于三维激光雷达扫描技术制作树冠三维模型提取测树因子，进而建立估算单木三维绿量的方法已被推广使用［7-8］。然而采用传统的人工测量方法，大面积测算三维绿量的工作包含多个环节，步骤复杂，精度有限。鉴于此，本文提出了结合高分遥感图像和地基三维激光雷达扫描数据，采用平面量推算立体量的方法测算城市三维绿量。

该方法首先在样地按树种类型不同进行分类，如分成针叶林地、阔叶林地和混交林地等，再基于实测的测树因子（冠幅、冠径、冠高、冠下高）得到各树种样地的实测三维绿量数值，将其作为因变量；然后对城市遥感图像进行景观分类并提取不同树种类型的空间覆盖面积和植被指数，如归一化差值植被指数（NDVI）等，将植被指数等因子作为解释变量；之后利用多元回归分析建立不同树种的三维绿量的遥感估算模型；最后结合各树种类型的分布面积，汇总得到研究区三维绿量总值。

由于三维激光雷达扫描获取技术能够更加快速、准确的采集测树因子，因此可以考虑在采用平面量推算立体量的方法时，采用三维激光雷达扫描的方法替代常规测量工作。使用三维激光雷达对单株树木进行不同方向的扫描，获取整株树木的完整点云数据，然后在计算机上将植株原型直接转换到三维数字模型，再通过相关软件提取冠幅、冠径、树高、胸径等测树因子指标。如此可大幅减少外业采集样本数据时间，并且测算结果受人为因素干扰小，时效性强［9-10］。

自2016年以来，我国很多城市都积极响应国家创建“森林城市”的号召，持续推进多期大尺度城市森林景观建设，拉萨、西宁、重庆、兰州陆续展开了打造“国家级环城生态公园”“西部特色山水城市”“森林进城、森林围城”等生态景观建设工程。采用本文提出的基于高分遥感图像和地基三维激光雷达扫描数据，利用平面量推算立体量的方法快速、准确、大面积测算城市三维绿量，对于准确评价绿色景观生态效益，科学提升城市生态系统服务价值，避免当前部分粗放式生态景观建设改造活动具有重要借鉴意义。

1 实验设计与方法

结合高分遥感图像数据和地基三维激光雷达扫描数据测算三维绿量，具体工作内容包括以下8个方面：（1）采集城市园林绿地系统规划建设图件、高分遥感图像、统计资料。（2）为配合地基三维激光雷达扫描系统作业，以及后续遥感图像处理，结合系统随机抽样和典型抽样的方法，在研究区布设包含主要绿化树种的阔叶林、针叶林、灌草丛的样地。（3）在样地对3个树种类别的主要树种每木检尺冠幅、冠高、胸径等测树因子；测算单株立木的三维绿量值。（4）利用地基三维激光雷达扫描系统对每个树种类别单株立木进行扫描，得到单木树冠外部点云数据，依据点云数据逼近树冠体积，得到各树种类别单株立木的三维绿量扫描值。（5）用各树种类别单株立木三维绿量的实测值校验、修正三维激光雷达扫描值。（6）结合高分遥感图像数据，抽取各树种类别单株立木的光谱特征、GIS特征。以三维绿量激光雷达扫描值作为因变量，以光谱特征、GIS特征（高程、坡度、坡向）作为自变量，采用逐步回归分析的方法，建立各树种类别单株立木的三维绿量多元线性回归方程。（7）结合地面实测三维绿量数据，对回归方程进行验证和精度评价。（8）依据每个树种类别的专题图谱，利用回归方程推算研究区各树种类别三维绿量，再统计出研究区全域三维绿量；绘制城市三维绿量等级分布图，分析评价城市三维绿量空间格局。

1.1 样地选取及样木选择

为确保样本点(单棵样木)具有代表性且在研究区均匀分布，配合地基三维激光雷达扫描系统作业，以及高分遥感图像几何精校正和分类精度验证，收集整理研究区园林绿地系统规划图、植物志、植物图鉴、园林植物配置表等图件、统计调查资料；然后结合系统随机抽样和典型抽样的方法，在研究区均匀布设包含主要树种的阔叶林、针叶林、灌草丛样地，并在样地中选择具有代表性的样木。具体实施过程分为以下8步：（1）结合实地调研并分析研究区城市园林绿地系统规划图，勾绘研究区生态网络线。（2）对于阔叶林、针叶林、灌草丛这3种类别，根据系统随机抽样理论，分别在生态网络线上随机布设各类别样地。（3）对于位于生态网络线上重要公园斑块、绿地斑块、林地斑块、农地斑块，以及沿江绿廊、交通沿线绿廊中布设包含主要绿化树种的样地。（4）在区园林绿地系统规划图上标定每一块样地具体的经纬度、高程、坡向、坡度；并勾绘斑块边界范围，斑块样地可以设置为正方形，样地面积为225m2；廊道样地可以设置为长方形，样地面积为500m2。（5）在每一块研究样地实地调查，观测样地是否包括各树种类别的主要绿化树种。选择生长状况良好，间隔稀疏，与周围立木树冠交叉重叠少、互相遮挡少的单株树作为样木。这有利于在业内处理点云数据和建立三维绿量多元回归方程。（6）确保样木在生态网络线上均匀分布且具有代表性。若不符合要求，记录位置、编号说明原因；之后返回上述流程第（2）步，在生态网络线上调整样地，重新选择样木。（7）对于阔叶林、针叶林、灌草丛这3个树种类别，每个类别至少选定150棵样木作为总样本数据集，然后按照训练样本集、测试样本集为70%:30% 的比例，从总样本数据集中随机抽取样本点（单棵样木）配置到训练样本集、测试样本集中。（8）每棵样木确定后，在遥感图像上标定使用激光雷达扫描仪扫描样木的测站位置（每棵样需3个测站，呈等边三角形分布）；检查每棵样木的编号、经纬度坐标及测站编号，确保正确无误后导入GIS数据库。

1.2 基于传统测量工具采集样木测树因子

在夏季的7～8月间在研究区按照预设路线在样地根据精准测树原理进行立地调查和每木检尺：使用全站仪、手持GPS、测距仪、3D 地面标识球、罗盘仪、胸径卷尺、围尺、测角规、数码摄像机等工具，记录数据采集样地的经纬度、海拔高程，样地编号；测量样地主要树种的胸径，冠高，东西、南北冠径，冠下高以及树高等测树因子的数值；全角度录像样木树冠形态；在内业将所测数据及相关内容作为样点属性导入GIS数据库；在内业依据所测数据计算主要树种类别单棵样木的三维绿量的实测值。

在上述工作中，每木检尺样木测树因子是一个重要环节，具体实施可以分为5步：（1）设定样地中心点，选择样地时需保证在外延8个方向上均为有林地，要求样地坡度小且视野开阔。记录样地主要树种、周侧景观、观测日期以及测量小组成员。（2）确定外延角规点和样地边界，使用罗盘仪、皮尺确认样地边界以及外延8个方向的角规点。并用长绳圈定样地的边界。（3）每木检尺，每个测量小组分配4～5名成员：其中2人配合使用全站仪测量树高和第一枝高，1人操作全站仪另1人则移动棱镜；另外2～3人配合移动棱镜测量每棵样木的东西冠径、南北冠径，以及测量样木胸径并设置标号。当林区不通视时，采用激光测高仪测量样木树高和树胸径。（4）测量样地参数，使用罗盘仪、皮尺测量样地坡度、坡向、坡位和植株密度、郁闭度。（5）内业数据处理，在外业测量工作完成后对每个样本点（单棵样木）的数据进行整理归档，然后作为样本点的属性数据导入至GIS数据库中。

1.3 基于LiDAR 点云数据采集样木测树因子

利用地基三维激光雷达扫描仪扫描样木，采集测树因子的实施过程可以分为5步：（1）已利用传统测量工具采集了全部样木的测树因子，在研究区高分遥感图像上标定每棵样木和3个测站的位置，检查确保每棵样木的编号、方位坐标以及测站编号正确无误，设定外业扫描样木的路线。（2）在晴朗无风，树木枝条不晃动的天气条件下，在样地架设测站，利用地基三维激光雷达扫描仪对每棵样木进行扫描。每单棵树木扫描需要3个测站，3个测站呈等边三角形分布，相互之间间隔夹角为120度（以待扫描树为参考）。（3）在内业对每棵样木的原始点云数据依次进行去噪、匹配拼接、抽稀处理。（4）参照在外业扫描样木时对样木的全视角录像，从点云数据中提取出样木的三维点云斑块，再通过旋转、缩放等操作，对样木三维点云斑块进行立体化判读检查，去除点云斑块中不属于树冠表面的数据点。（5）依据剩余的点云数据，提取单棵样木的胸径，冠高，东西冠径、南北冠径、冠下高以及树高等测树因子的数值。

在上述工作中，第（2）步对样地样木进行扫描是一个重要环节，需注意完成以下三点：第一,为配合后续处理内业数据，在扫描单棵样木前需用数码摄像机对样地周侧环境、目标树种进行全视角录像，然后画出目标样木树冠草图，并标注序号和说明文字；第二，为确保点云数据的完整性和密集性，在样木树冠下方四周外侧扫描时，让扫描仪避开密集的树枝遮挡，并以仰视的角度环绕植株扫描。可以架设多个测站以保证扫描数据的密集性；第三，为确保后期点云数据能够被正确拼接，在对单棵样木的扫描的过程中，必须保证相邻两个测站至少有3个共同标靶。

1.4 误差分析与样本筛选

基于样木点云数据，可以提取出单棵样木的胸径，冠高，东西冠径、南北冠径、冠下高的数值，藉此计算树冠体积，进而建立单木三维绿量的回归方程。因此，为保证三维绿量的回归方程的准确度，有必要将样木测树因子的实测值与基于点云数据的采集值进行对比分析；然后从三维激光雷达扫描仪使用、内业点云数据处理、测树因子提取这3个方面入手，查找出现偏差的原因，再做出相应改进，重新采集样木的各测树因子的数值。具体实施过程可以分为4步：（1）把利用传统测量工具所得样木各项测树因子数值作为实测值；把利用LiDAR 点云数据提取样木各项测树因子数值作为采集值。对每一棵样木，计算每个测树因子采集值与实测值之间的相对误差（采集值-实测值）/实测值×100%），以及各测树因子的平均相对误差。（2）标识某项测树因子相对误差超过5% 的样木；标识平均相对误差超过5% 的测树因子项。（3）从三维激光雷达扫描仪使用、内业点云数据处理、测树因子提取这3个方面入手，查找出现偏差的原因；做出针对性改进，对不合格数据的样木重新进行扫描，重新采集各项测树因子值。（4）再次计算测树因子采集值与实测值之间的相对误差，以及各测树因子的平均相对误差。若经过多次扫描量测后，某棵样木测树因子值的相对误差仍然超过5%，则剔除该样木数据，以确保后续计算树冠体积的精度。

1.5 单木三维绿量计算

采用AutoCAD 建模计算法测算样木树冠体积，算法流程分为4步：（1）截取三维激光雷达扫描仪获得的点云数据，进行预处理，获取单棵树木的树冠点云数据。去除第一主枝（或第一个分枝）以下的部分的点云数据。（2）对获取单棵树木的点云数据进行旋转、缩放等操作，采用人机交互的方法进行立体化判读检查，除去点云中不属于树冠表面的点。（3）将点云数据导入软件AutoCAD，再AutoCAD 在上加载软件Point Clouds，利用Point Cluods使点云数据生成Mesh网格，然后得到适配线自动拟合出整个单木的树冠三维模型。（4）在构建出单木的树冠三维模型之后，利用查询命令选择计算模型测算体积，则可获得单木的三维绿量值。

1.6 单木三维绿量遥感回归方程构建

从高分遥感图像上提取出阔叶林地、针叶林地、灌草丛的空间分布信息，然后利用掩膜运算等操作分别制作这三个类别的专题信息图。基于专题信息图为各类别样木提取光谱特征，用于构建单木三维绿量的遥感回归方程；在确立各类别单木三维绿量的回归方程后，利用专题信息图测算本类别全域三维绿量。具体实施过程可以分为6步：（1）对样地进行实地调查，利用全站仪测定每块样地测量点的经纬度坐标，将离散分在遥感图像上的测量点作为地面控制点，对高分遥感图像的所有波段进行几何精校正。（2）将高分遥感图像的全色波段和多光谱波段进行像素级融合，融合后的图像空间分辨率为2m。（3）对同期多景相邻融合后的图像进行拼接，得到研究区全域遥感多波段图像。（4）结合实地调查在遥感图像上为非植被地类（水体、不透水面、砂石地，其它非植被类）、植被地类(阔叶林地、针叶林地、灌草丛)建立解译标志；结合分层分类和基于支持向量机的监督分类方法识别以上各个地类；基于分类结果提取出各植被地类的空间分布信息图像，分别设置为一个图层。（5）采用人机交互目视判读的方法检查分类结果；制作误差混淆矩阵，计算分类图像的制图精度、用户精度、总体精度、Kappa系数进行定量评价，检验空间分布信息图像对各植被类型的识别精度。检验合格则统计各植被地类的面积；若不合格返回第（4）步，重新进行图像分类。（6）依据植被分类结果，对全域高分遥感图像进行掩膜运算、符号注记、配置图面要素等处理，生成3个树种类别空间分布信息图像，并作为3个特征波段图像加入融合后的遥感图像中。为下一步的为提取建模特征图像做准备。

1.7 单木三维绿量回归方程构建

对阔叶树、针叶树、灌草丛这3个类别，从高分遥感图像确定各自全部样本点（每棵样木）的光谱特征；获取研究区的数字地面模型（digital elevation model，DEM），从DEM图像上提取同名像点的GIS特征（高程、坡度、坡向）。对于每个树种类别，按70%:30% 的比例，从各类别总样本数据集中随机抽取样点配置到训练样本集和测试样本集；以训练样本集中样本点的光谱特征、GIS特征作为自变量；三维绿量实测值作为因变量，利用逐步回归的方法进行分析，拟合得到每个树种类别单木三维绿量的回归方程；通过相关系数与F 值检验评估模型精度；基于测试样本集评价回归方程的预测精度。具体实施过程可以分为7步：（1）对3个树种类别，在高分遥感图像上分别确定各类别每棵样木（样本点）的经纬度坐标，对样本点像元进行缓冲区处理（以坐标点为圆心,以样木冠径值的1/2为半径），然后分别提取样本点缓冲区内的遥感图像像元的近红外波段、红波段、绿波段的灰度值；再对近红外波段、红波段灰度值进行运算，进一步得到样本点缓冲区像元的归一化差值植被指数（NDVI-Normalized Difference Vegetation Index）、比值植被指数（RVI-Ratio Vegetation Index）、改进型土壤调节植被指数（MSAVI-Modified Soil Adjusted Vegetation Index）；从DEM图像上提取同名像点的高程、坡度、坡向。每个样本点的以上9项数值作为对应单棵样木的遥感/GIS特征。（2）在每个树种类别会议训练样本集内，将样本点以上9个光谱特征和GIS特征值设为各自类别的自变量；将已获取了每个类别样本点的三维绿量实测值，设为各自树种类别的因变量。（3）对3个树种类别，分别利用各自测试样本集，对因变量和与其对应的9个自变量进行相关分析，计算Pearson简单相关系数，根据相关系数值去除与因变量无强线性关系的自变量。（4）采对3个树种类别，都利用逐步回归分析的方法筛选并去除引起多重共线性的自变量，拟合本类别单木三维绿量的回归方程：对每个类别，首先用因变量对每一个自变量做线性回归；然后以对因变量贡献最大的自变量所对应的回归方程作为基础，逐次引入其余自变量。每引入一个新自变量进行F 检验，对已选入回归模型的原变量逐一进行t检验，去除经t检验认为不显著的原自变量（设定相伴概率值Sig0.05时，将自变量引入回归方程；Sig0.10时将自变量从回归方程中去除），以确保当前自变量集中每一个变量都为显著自变量；将迭代过程执行若干步后确保既无显著自变量被选入回归方程，也无不显著自变量从回归方程被去除。至此保留在回归模型中的自变量既重要又相互之间无严重多重共线性。由此方法得到3个树种类别单木三维绿量的回归方程。（5）通过相关系数与F 值检验评估模型精度：对每个树种类别，利用测试样本集数据和回归估算数据，计算回归方程的相关系数并在方差分析的基础上利用F 检验，对回归方程进行显著性检验。考查回归方程整体上是否能够通过P 值在0.05置信水平或P 值在0.01置信水平的显著性检验；评价回归标准残差直方图是否服从正态分布。（6）对3个类别的回归方程，进一步评价回归方程的预测精度：计算三维绿量实测值与回归方程三维绿量估算值之间的均方根误差；然后计算0.05置信水平和0.01置信水平下的绝对误差限和相对误差限；再基于相对误差限计算得到回归方程的预测精度。（7）对阔叶树、针叶树、灌草丛这3个类别，检查各类别回归方程的预测精度是否大于85%.若某一树种类别回归方程的预测精度小于85%，分析建模过程疏漏之处查找原因，然后返回第（1）步重新设置训练样本集、测试样本集构建回归方程。

1.8 制作单位面积三维绿量等级分布图

制作单位面积三维绿量等级分布图的实施过程可以分为3步：（1）已分别制作出了阔叶树、针叶树、灌草丛这3个类别的空间分布信息图像。基于各类别空间分布信息图，对参与各回归方程的每个自变量图像（如NDVI图像、高程图像）进行掩膜运算；对每个类别，将所有自变量图像作为特征波段，并与本类别遥感专题信息图组合在一起构成专题图谱。至此，阔叶树、针叶树、灌草丛这3个类别均形成了与回归方程对应的专题图谱。（2）对各树种类别，从专题图谱上提取每个数据点（单颗树木）的各个自变量的数值，将其输入回归方程，计算得到该点的三维绿量的测算值；根据测算结果制作本类别三维绿量专题信息图。然后加入到本类别专题图谱中。（3）综合3个树种类别三维绿量专题信息图，统计整个研究区三维绿量值。制作研究区单位面积三维绿量等级分布图。

2 结论

相较于以平方米为量纲的二维绿量指数，三维绿量从立体空间的角度评定植被空间结构差异和生态功能水平高低，因而能够更准确的反映城市绿色景观格局变化。然而，采用传统的人工测量方法，大面积测算三维绿量的工作环节多，步骤复杂，精度有限。为此，本文提出了基于高分遥感图像和地基三维激光雷达扫描数据，采用平面量推算立体量的方法测算城市三维绿量。该方法首先利用地基三维激光扫描点云数据建立样本点单木三维模型，再利用模型提取单木三维绿量；然后利用高分遥感图像和DEM图像获取样本点的光谱特征和GIS特征作为自变量；之后利用逐步回归分析的方法筛选并去除引起多重共线性的自变量，拟合单木三维绿量的回归方程；最后结合不同植被类别的空间分布信息遥感图像，推算出研究区全域的三维绿量总值。该方法对于快速测算城市三维绿量提供了一个新的思路，可为准确定量评价城市绿色景观生态效益提供理论依据和技术支持。