一种利用激光点云计算三维森林冠层消光系数的方法

摘要

本发明提供了一种改进的利用计算几何算法来计算三维森林冠层任意给定入射光线下的消光系数的方法，属于森林冠层结构参数获取方法的研究领域。其步骤为：植被冠层的三维激光点云数据的获取及预处理；点云数据的三维网格化；基于体元数据结构的点云切片算法；点云切片的线采样分析；第k层切片的平均投影系数的计算；各层切片任意给定入射光线下的消光系数及整个冠层消光系数的计算。本发明与传统观测手段相比，工作量小，无需接触式观测，不破坏冠层结构和辐射特性，具有客观高效精确的特点；开发了从激光雷达数据中提取三维结构和生物物理多样性信息的方法，将叶片的水平和垂直分布变化规律特征化。

说明书

一、技术领域

本发明涉及一种利用地面三维激光扫描仪获取的点云数据计算森林冠层消光系数的方法， 具体地说，是指一种改进的利用计算几何算法来计算三维森林冠层任意给定入射光线下的消光 系数的方法(流程如附图1所示)。

二、背景技术

森林冠层结构在土壤-冠层-大气的相互作用中占有重要的地位，并与生物地球化学循环中 的物质能量交换密切相关。消光系数定义为单位叶面积在垂直于光线方向的平面上的平均投影 面积，是定量的理解森林冠层结构和辐射传输的重要因子。他是基于叶倾角的方向性分布函数 (包括倾角和方位角分布函数)，对于任意给定角度的入射光线来计算的。

间接、准确地获取植被冠层的叶倾角分布密度函数一直以来是一个具有挑战性的工作，传 统的方法主要是利用逐叶片直接观测的方法来获取给定单株树的叶倾角分布密度函数。比如： Norman等人利用量角器和指南针等设备制作了一种较为简单方便的测量叶片角度和空间分布 的仪器(附图2左)；Lang等人利用机械臂的灵活性和巧妙的几何传递关系制作了一种可以自 由移动和伸缩的用来测量叶片空间和角度分布的仪器(附图2右)。

但由于其工作量大，主观性强，而且对所测量的目标树木的高度有限制，所以不能够广泛 地应用于实际研究中。更为重要的一点是，直接观测往往会影响甚至破坏植被冠层的结构，进 而破坏冠层内部和以下的辐射分布情况。而且就森林而言，实际中很难利用直接观测仪器来进 行逐叶片测量。因此，需要一种新的快速、准确、间接的方法来获取冠层叶片的角度和空间分 布情况。近年来，也有人开始尝试进行间接地获取叶片倾角分布。例如，2007年Shilbayama 和Watanabe发表在期刊《Plant Production Science》第10卷的“Estimating the mean leaf  inclination angle of wheat canopies using reflected polarized light”一文中，提出利用反射极化光 来进行平均叶倾角的估算方法。.

目前，还有一种较为流行的方法是利用数学几何模型来近似地逼近真实树木的叶片实际分 布，其中被广泛接受的是Campbell在20世纪80年代提出的椭球体儿何模型。其核心思想是 椭球体表面各切面倾角的不同分布来近似表达实际叶片的角度分布，并通过调整椭球体的长、 短轴来调节不同倾角分布密度函数的变化。该方法可用以下公式表达：

$\xi \left(\alpha \right)=\frac{2{\chi }^3\sin \alpha }{\Lambda {({\cos }^2\alpha +{\chi }^2{\sin }^2\alpha )}^2}---\left(1\right)$

其中$\chi =b/a,\Lambda =\left(\begin{array}{cc}\chi +\left({\sin }^{-1}ϵ\right)/ϵ,& \chi <1,ϵ={(1-{\chi }^2)}^{1/2}\\ \chi +\frac{\ln [(1+{ϵ}_1)/(1-{ϵ}_1)]}{2{ϵ}_1\chi },& \chi >1,{ϵ}_1={(1-{\chi }^{-2})}^{1/2}\end{array}\right),0\le \alpha \le \pi /2,$

当χ＝1时，该分布将成为球面分布，此时Λ＝2，a、b分别是椭球体的半短、半长轴。但 这种方法只是在理论上近似的逼近叶片的实际分布，不足以描述变化多端的真实情况。

地面激光雷达扫描仪的出现使得我们可以直接从三维角度获取高精度(毫米级)的森林冠 层结构信息，进而计算冠层的消光系数。目前国内关于激光雷达在森林中的应用研究，特别是 针对三维冠层结构的研究还处在起步和探索阶段。已有的研究大部分集中在利用地面激光雷达 进行单株立木或森林样地水平的基本每木检尺的参数(包括胸径、树高和生物量等)提取方面。 例如：2011年马利群等人在《世界林业研究》第1卷的“激光雷达在森林垂直结构参数估算中 的应用”一文中介绍了激光雷达在提取森林垂直结构参数方面的应用；在激光雷达的森林应用 方面，生成高精度的数字地面模型是其中关键的一步，周淑芳等人在2007年《遥感技术与应 用》第22卷的“基于机载激光雷达数据的DEM获取及应用”一文中探讨了利用机载激光雷达 系统生成数字地面模型的方法和精确评价。除此之外，庞勇等人在2006-2009年间对利用激光 雷达系统提取森林单木和样方水平的树高估算进行了大量较深入的研究。但目前国内对利用激 光雷达所获取的三维点云数据提取冠层消光系数的研究理论和技术还不够深入，需要进一步加 强。

三、发明内容

本发明的目的是：

提供一套直接利用地面激光雷达系统生成的三维点云数据，结合计算几何的方法，从三维 角度来计算任意给定入射光的单株树木或森林样方消光系数的算法。并提供非接触式的快速准 确观测森林冠层结构特征参数的方法。

本发明的原理如下：

利用较新的遥感技术手段(地面三维激光雷达扫描系统)，结合计算几何的技术手段，直 接从三维的角度获取给定方向入射光下的单株树木和森林样方尺度叶片的消光系数，进而计算 植被冠层的消光系数。首先建立三维空间的数学模型来精确定位单个叶片的位置，在此基础上 采用计算几何方法，对三维点云数据进行网格化，利用点云切片算法对不同高度的点云进行分 类。然后进行线采样分析，获取叶片截获光线的相关信息，进而计算每个切片层的投影系数， 即消光系数，从而得到整个冠层的消光系数。完成给定入射光线下的消光系数及冠层整体的平 均投影系数的计算。

本发明的技术方案主要包括以下步骤：

(1)首先利用地面激光雷达扫描系统，获取植被冠层的三维点云数据。其中包含了扫描 目标点的空间几何和激光束弹回的能量信息，三维坐标直接提供了任意一点的空间位置坐标信 息，这也是传统光学理论中的数学模型的基础。获取的点云数据首先进行图像拼接，并手动去 除地面点云。

(2)三维网格化。首先在点云区域定义一个X，Y，Z三个坐标轴的笛卡尔坐标系，将所 有点云数据包含其中。并将点云数据分为有限的小区域，建立以体元(voxel)为基础的数据结 构，这个过程叫做三维网格化。每个体元由长(l)宽(w)、高(h)三个参数决定其大小， 所有的点云数据被分为m×n×p个体元，其中，m＝(X_max-X_min)/w，n＝(Y_max-Y_min)/l，P＝(Z_max-Z_min)/h。 如果，l＝w＝h，则体元为一个立方体。以这些点云数据的三维中心点，即[(X_max-X_min)/2，(Y_max-Y_min) /2，(Z_max-Z_min)/2]为新的原点，建立笛卡尔坐标系。坐标轴Z为树干的生长方向，垂直于水平面， X、Y轴位于与Z轴垂直的平面上，Z轴为纵向的轴，X、Y轴为横向的轴。如附图3(a，b，c) 所示。

(3)点云切片算法。三维网格化过程结束后，点云数据被分为纵横方向的众多数据层， 即切片，点云数据可视为纵横方向切片的叠加。如附图3(d)所示水平切片。所有体元设为(i， j，k)(i＝1，2，...，m；j＝1，2，...n；k＝1，2，...，p)。当k＝1，i和j为给定区域内的任意值，则可表示为第 一层或第一横切面的所有体元。

为了近似模拟入射光线与森林冠层的相对位置关系，发展了另外一种旋转点云的切片方 式。首先，将点云数据旋转至我们要模拟的相对位置，然后对其进行三维网格化。这样水平切 片与XY轴所在的平面平行，垂直的切片与Z轴所在的平面平行。我们假设太阳直射光从天顶 方向射入，与Z轴方向平行。当入射光方向改变时，可通过旋转点云数据来实现。

具体的说，三维点云数据的中心点为笛卡尔坐标系的原点，点云数据则根据点云与切片相 对位置的需要围绕X、Y、Z轴旋转。例如，将点云数据围绕Y轴旋转30°，相当于点云被切 片或太阳入射光以30°的倾斜角进行切片。因此，点云可在任意位置(水平0°-360°，垂直0°-90°) 被切片，称作方向性切片。任何点云数据都可通过方向性切片生成薄板以提取三维冠层信息。

(4)线采样分析。叶片在切片中的空间分布通过线采样方法进行分析。点云数据每个位 置的切片都是由许多样线组成的，这些样线垂直或平行于太阳入射光r(θ，β)。当研究沿r方向 叶片的分布时，是考虑与r平行的样线。如果将每个体元的三维参数设置为激光雷达采样距离 的1.5倍，我们可以保证每个体元最多只包含一个样点。那么，每条样线中非空的体元(N) 个数代表样线与叶片相交的次数，且N服从随机分布。如果大于采样距离的1.5倍，一个体元 内包含多个样点时，我们需要将样线上所有的点投影到与样线方向平行的面上，计算不同的N 值，作为与叶片相交的次数。

定义P_n为样线穿过整个点云数据时，与叶片相交n次的可能性。P₀表示样线穿越的都是空 隙，叠加指数为0。通过求解所有样线穿过植被时N的平均值，可以求得整个点云数据中所有 样线与叶片相交次数的平均值(m)和方差(σ²)。这个相交的次数并不是实际的叶片点云密度， 而是投影到与样线方向垂直的平面上的投影。每条样线的方差(σ²)/均值(m)的相对变化可 以反映叶片在真实空间的分布状况。单株树或林分的叶片分布可根据其相对变化(σ²/m＝1，＞1， 和＜1)相应的分为三种情况：规则分布、随机分布和聚集分布。

线采样分析是计算整个冠层和不同高度切片消光系数的关键，同时有利于点云密度、辐射 通量和辐射分布状况等参数的获取。

(5)太阳入射光倾斜角(θ)固定时，计算第k层切片的平均投影系数。当只研究r(θ，β) 一个方向的太阳入射光时，也就是说第k层切片的投影系数G_k(r)是指固定的r(θ，β)方向的投影 系数。当不考虑入射光方位角，而倾斜角为θ时的k层平均投影系数G_k(θ)，是通过方位角β 在[0，2π]上的积分计算的：

$G_k\left(\theta \right)={\int }_0^{2\pi }{G}_k\left(r\right)\mathrm{d\beta }---\left(2\right)$

(6)计算任意给定入射光线下的投影系数，即消光系数。首先考虑以固定方向r(θ，β) 进行切片，在点云数据完成三维网格化之后，所有的三维点云都被包围在立体边界内，生成水 平和垂直的切片及样线。如果所有的切片与方向r(θ，β)平行，那么所有的样线也平行于方向r(θ， β)，且纵轴的方向即代表此方向。每个体元可以构建一个垂直于这个主轴的平面，所有的样点 都可以投影到这个平面上。在三维网格框架中，这个投影面是两个垂直于样线方向的体元边界 之。设，每个非空体元内包含的样点个数为n：

a.若n＝1，每个二维的正方形以用户采样间距为边长，则样点代表这个正方形叶片区域的 中心。若采样间距为s，则叶片区域面积为s×s。投影系数为以取样间距为边长的正方形面积(即 s×s)，与体元的长度(l)×宽度(w)的面积之比。

b.若n＝2，可以构建一条样线来表示叶片区域，叶片投影面积为这条样线的投影长度乘以 正方形的单位面积。投影系数为投影后和投影前的样线长度比。

c.若n≥3，首先检测是否所有的点都在一个面或是一条线上。如果在同一样线上，则将这 些点分解到更小的尺度上(n＝2)。如果在同一个面上，则构建一个三角形来表示叶片区域，叶 片投影面积由三个顶点的投影坐标来计算。三个投影坐标设为(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，三角形的 投影面积A_p即叶片投影面积，由以下矩阵计算：

$A_p=\frac{1}{2}\times \left(\begin{array}{ccc}{x}_1& y_1& 1\\ x_2& y_2& 1\\ x_3& y_3& 1\end{array}\right)---\left(3\right)$

d.若体元内的所有点都不共面或在一条线上，则构建一个三维的凸包，凸包表面积的一 半即表示实际的叶片面积(A)。然后，将所有的点投影到与样线平行的平面上，用这些点建立 二维的凸包，其面积即表示投影的叶片面积(A_p)

综上所述，投影系数A_p/A是由非空体元计算。第k层垂直于入射光的投影系数G_k(r)可由 以下公式计算：

$G_k\left(r\right)=\frac{1}{N}\underset{\xi =1}{\overset{N}{\Sigma }}{A}_{\mathrm{p\xi }}/A_{\xi }---\left(4\right)$

其中N为第k层切片非空体元的总个数，k为切片的层数。

求出G_k(r)后，代入公式(2)，求出各切片层的平均投影系数，整个冠层的投影系数(即 消光系数)为各层切片投影系数的总和。

与现有技术相比，本发明优势体现在：激光雷达技术为从三维角度获取森林冠层结构参数 提供了有效的技术手段，结合计算几何方法进行间接观测，与传统观测手段相比，工作量小， 无需接触式观测，不破坏冠层结构和辐射特性，具有客观高效精确的特点；开发了从激光雷达 数据中提取三维结构和生物物理多样性信息的方法，将叶片的水平和垂直分布变化规律特征 化。

具体有益效果如下：

本发明利用激光扫描仪获取的三维点云数据，以体元为基础计算森林冠层的消光系数，提 供了一种间接的无需接触的观测森林冠层生物物理参数的方法。本发明通过旋转点云数据而保 持太阳光从树冠上方沿纵轴平行方向入射，来模拟现实情况中太阳光在不同角度照射树冠的情 景，既简便又高效，且不会对森林结构和辐射特性造成任何不良影响。本发明不仅可以计算冠 层整体的消光系数，还可以计算树冠不同高度切片的消光系数，这有利于进一步研究光照辐射 在冠层内部的截留和分布机制，以及其他生物物理参数的研究。本发明从三维空间对森林冠层 机构进行观测，其关键是利用几何投影技术将三维点云数据转换为二维栅格图像，从不同方向 计算冠层不同高度平面的消光系数，不仅极大提高了估算精度，而且为冠层内部的辐射传输机 制的研究提供了大量信息。

实际应用表明，本发明提供了间接观测森林冠层消光系数的有效方法，无论是单株树木还 是森林样方冠层的消光系数，均可直接使用本方法求得。该方法克服了传统方法需花费大量时 间和大量人力物力，影响森林冠层结构特性，且误差大的缺陷。提高了植被生物物理参数估算 的效率，增强了三维激光扫描技术应用的普及性和有效性。可以更好地服务于森林资源调查、 植被生态遥感等资源环境研究项目。

四、附图说明

图1为本发明的流程示意图；

图2为传统的叶片角度和空间分布观测仪器；

图3为体元数据结构示意图和水平切片示意图

a.两个单独的体元，以及决定其大小的长、宽、高三个参数；

b.一个2×2×2的三维栅格空间；

c.体元结构的三维空间网格化示意图；

d.基于体元结构的点云水平切片算法示意图；

图4为一株北美花旗松的三维点云数据；

图5为三维网格化和点云切片示意图

a.针对单株北美花旗松的基于体元结构的点云水平切片算法示意图；

b.单层切片的体元三维空间网格化(5×5×5个体元)示意图；

图6为单株立木的三维点云数据在三维网络化的基础上的线采样示意图

a.“线采样”(图中红色长方体)的提取和分析；

b.线采样提取样本的横断面；

c.线采样提取样本的直方图分析。

图7为单个体元内部叶片消光系数计算方法示意图

a.单个体元内部叶片三维点云示意图；

b.叶片三维凸包结果及面表达；

c.叶片三维凸包结果及点、线表达；

d.叶片在水平面上投影点集的二维凸包计算结果及面表达。

图8为14m-16m高度的点云切片层消光系数分布特征示意图

a.本切片层点云分布示意图；

b.本切片层三维消光系数空间分布特征示意图；

c.本切片层二维点云密度空间分布特征示意图；

d.本切片层二维消光系数空间分布特征示意图。

五、具体实施方式

以下通过实例对本发明作进一步解释：

以一株北美花旗松(Douglas-fir)为研究对象(树高约25m)，使用地面三维激光扫描仪Leica  ScanStation 2(其参数如表1所示)和高精度GPS，在树的一个侧面进行三维点云数据的采集， 采样间距为2cm。手动去除地面点云及其他噪声点云，得到单株北美花旗松的三维点云数据， 如附图4所示。

表1 三维激光扫描仪Leica ScanStation 2参数

在取得花旗松的点云数据并预处理后，如附图1所示，采用计算几何技术对三维点云数据 进行处理，第一步为三维网格化，第二步为基于体元的点云切片，第三步为线采样分析。

根据技术方案步骤(2)所述，对点云数据进行网格化，建立体元数据结构，每个体元大小 为0.5m×0.5m×1m。如附图5(b)中所示，为单层切片的体元三维空间网格化(5×5×5个体元) 示意图，图中的点云颜色代表不同高度，所示点云为附图5(a)中高度为17至18米间的点云。

对三维网格化后的点云数据进行切片，如附图5(b)为针对一株北美花旗松的基于体元结 构的水平点云切片算法示意图，灰色为水平的7个切片平面，点云的不同颜色代表不同的高度。 为了模拟太阳光从不同方向照射到树冠的情形，我们将太阳入射光固定在与Z轴平行的正上 方，通过旋转花旗松的点云数据，来模拟太阳入射光与树冠相对角度的变化，以计算任意角度 入射光下的投影系数。

根据技术方案步骤(4)所述，对三维网格化的点云数据进行线采样提取和分析。我们对 花旗松点云数据进行垂直方向的线采样提取分析，如附图6(a)所示，红色长方体为“线采样”， 附图6(b)为线采样提取样本的横断面，附图6(c)为线采样提取样本的直方图分析，可见 14m-18m间的点云数据最多。

为更好地描述单个体元内部叶片的消光系数计算方法，下面以大叶枫树的一个叶片为例进 行说明，如附图7所示，7(a)为枫树叶片的三维点云。根据技术方案步骤(6)中d所述， 对叶片点云构建一个三维凸包，如附图7(b)所示，为凸包的计算结果与面表达。附图7(c) 为凸包的计算结果及点、线表达。计算三维凸包的表面积，其表面积的一半即表示实际的叶片 面积(A)。然后将三维点云投影到与太阳入射光垂直的平面上(即与Z轴垂直、与XY轴平行 的平面)，用这些点建立二维的凸包，其面积即表示投影的叶片面积(A_p)，如附图7(d)所示。 A_p/A即该叶片的消光系数。

根据本发明提出的算法，对花旗松点云进行分析，依照技术方案步骤(5)和(6)所述， 获取不同高度切片中每个体元的点云密度空间分布和消光系数。以高度为14m-16m处的切片 为例，如附图8所示：图8(a)为本切片层点云分布示意图，图8(b)为本切片层三维消光 系数空间分布特征示意图，图8(c)和8(d)分别为本切片层二维点云密度和消光系数空间 分布特征示意图。

从以上结果可以看出，本发明还可以提供树冠叶片密度垂直变化及其特定高度水平分布的 三维信息，是其他方法难以比拟的。