CN106248003B - 一种三维激光点云提取植被冠层聚集度指数的方法 - Google Patents

Abstract

本发明属于激光雷达遥感技术领域，具体为一种三维激光点云提取植被冠层聚集度指数的方法。本发明利用遥感技术手段，获取植被样方冠层三维点云数据，通过构建三维体元模型、转换坐标系统、计算冠层间隙率等过程，建立基于地面高分辨率激光点云数据提取植被冠层的聚集度指数的方法。该方法可快速、准确的提取植被冠层聚集度指数，数据采集简便，且不受观测时光照条件的影响，在研究时也不用考虑卫星遥感数据或产品的像元尺度问题。并且不会对植被结构和辐射特性造成任何不良影响，同时还可以永久性的记录植被样地的三维结构特征。本发明简便，高效，对植被无副作用，且相对现有技术极大地减少了计算量。

Description

一种三维激光点云提取植被冠层聚集度指数的方法

技术领域

本发明属于激光雷达遥感技术领域，涉及一种利用地面三维激光扫描仪获取的点云数据来评估森林冠层聚集情况的方法，具体为一种三维激光点云提取植被冠层聚集度指数的方法。

背景技术

植被冠层是植被与外界环境相互作用最直接和最活跃的界面层，对生态系统物质、能量交换，生物多样性，气候变化等具有重要的影响。冠层结构是冠层研究的一个重要方面，对植被冠层结构的准确描述是理解植被生态系统格局、过程及其运作机制的重要基础。事实上，植被的冠层结构不是随机分布的，植被的冠层叶片会因为空间的限制而出现不同程度的聚集。聚集度指数(clumping index,Ω)是一个重要的植被冠层结构参数，表征了冠层的空间分布集聚特征。聚集度指数描述了有效叶面积指数与真实叶面积指数的偏离程度，是精确获取叶面积指数的重要参数。此外聚集度指数还可以区分冠层“光照叶”和“阴叶”，从而改进各种地表过程模型。激光雷达(Light Detection and Ranging，LiDAR)，是近年来国际上发展十分迅速的主动遥感技术，可以快速获取目标物的三维结构信息，在反演和提取各种生态物理参数的研究中取得了成功的应用。

高空间分辨率的半球图像可用来反演聚集度指数(Walter 2009)。利用半球图像来提取聚集度指数的研究有很多，Chen和Cihlar(1995,CCI)提出利用间隙率和间隙大小分布计算聚集度指数，这种方法最初被用到叶面积指数测量仪器TRAC中，后来被用于校正半球图像。Lang和Xiang(1986)提出一种基于对数间隙平均的方法(CLX)来计算聚集度指数。Leblanc等人(2005)提出一种新的计算聚集度指数的方法，这种方法是由CCI和CLX两者结合而得的，处理了先前方法中的许多限制。Walter等人(2003)也提出对CLX方法的一种分层校正方法。Pielou提出Pielou的空间分段系数(pielou,1962,PCS)计算聚集度指数。对于半球图像，以上的研究方法基本都能达到较好的预期效果。不过这类光学遥感技术存在诸多外界因素的影响，如拍摄时的光线条件等。基于多角度的卫星遥感数据或产品也可实现聚集度指数的计算，如POLDER，MODIS的BRDF产品等，利用热点和暗点之间的归一化差值NDHD进行计算(Lacaze和Roujean，2001；Lacaze等，2002；Chen等，2005；Simic等，2010；Pisek等，2011)。目前，没有像元尺度聚集效应的描述因子，这类研究通常只是将像元中主要植被类型聚集指数的经验值作为整个像元的聚集度指数(Plummer等，2005，2006)，这种估算方法没有考虑像元(尤其是混合像元)内的不均匀性，并且精度上也存在着很大的不确定性。另一方面，由于遥感资料的限制,基于多角度遥感的聚集度系数反演理论与技术还存在诸多困难。目前，已有部分研究利用激光雷达技术实现了聚集度指数的提取。在地基激光雷达方面，Moorthy等人(2011；2008)基于Chen和Cihlar(1995)提出的间隙大小分布理论计算聚集度指数。通过激光束的拦截信息来模拟研究区域的空间分布并比较树冠的间隙率。这类方法多是结合激光雷达技术和半球摄影技术将点云数据模拟成半球图像而实现的。此外，Zhao等人(2012)利用地基扫描仪ECHIDNA得到的波形数据估算了聚集度指数。该方法同样是基于Chen和Cihlar(1995)的间隙大小分布理论得到的。在机载激光雷达方面，Thomas等人(2011)基于包括平均、中间、标准误差的几个机载激光扫描尺度运用了一个新进的方法计算聚集度指数。

较之于光学等其他技术领域，激光雷达领域对于聚集度指数的研究有其它技术领域所不能比拟的优势，不过该技术领域对于聚集度指数的研究还不是很成熟，还有较大的进步空间。与机载激光雷达技术相比，地基激光雷达数据的获取较为简便。此外，目前存在的一些利用激光雷达技术计算聚集度指数的研究中，还没有很好的能直接从点云数据中提取聚集度指数的方法。因此在利用地基激光雷达技术反演聚集度指数这方面有较好的研究前景。

地基激光雷达，作为一种主动遥感技术，具有分辨率高、光斑小、携带便捷等特点，能够以非接触方式快速、高精度地从地面测量树冠层的内部结构，获取海量点云数据。利用地基激光雷达技术反演冠层聚集度指数在一定程度上克服了其它技术领域存在的一些缺点。本发明利用实验中所得的冠层的三维点云数据来研究一种三维激光点云提取植被冠层聚集度指数的方法。

发明内容

针对上述存在问题或不足，为解决光线条件、尺度选取、数据获取的不便等问题，本发明提供了一种三维激光点云提取植被冠层聚集度指数的方法。

具体技术方案如下：

步骤1、利用地面激光雷达扫描系统，获取植被冠层三维点云数据：

首先，在目标区样方中心点和样方外侧分别架设三维激光扫描仪，以样方中心点为观测点所得的坐标系统为标准，将得到的多站点云数据进行点云配准。

然后，在水平方向将样方的三维点云数据切割为以样方中心点为圆心，3≤r≤10m为半径的圆形区域；再将所有低于三维激光扫描仪高度的点剔除，得到植被冠层三维点云数据。

步骤2、三维体元模型构建：

依据步骤1得到的冠层三维点云数据求出笛卡尔坐标X、Y、Z的最小值(X_min，Y_min，Z_min)和最大值(X_max，Y_max，Z_max)，以X、Y、Z的最小值(X_min，Y_min，Z_min)为起始点，以体元大小为步长划分冠层三维点云数据，并确定点云在体元坐标系中对应的体元坐标值及体素值。体元大小由体元的长L、宽W、高H决定，整个数据区域被划分为N_L×N_W×N_H个体元，其中，N_L＝(X_max-X_min)/L，N_W＝(Y_max-Y_min)/W，N_H＝(Z_max-Z_min)/H。点云体元化后的坐标值由以下公式得到：

式中，int是取整符，直接取出小数前面的整数部分，(i，j，k)是点云数据笛卡尔坐标(X，Y，Z)对应的体元坐标，体元大小L×W×H与扫描采用的点间距一致。

体元的体素值通过判断体元内包含的激光点个数来确定，如果体元内激光点个数大于等于1，代表激光束被体元拦截，体元体素值赋为1，否则体素值赋为0。

步骤3、坐标系统的转换：

将体元化后的冠层三维点云数据从笛卡尔坐标系转为半径为1的球面坐标系统。去除重复的体元，即确保转换后的球面坐标系统中每个点只有一个体元。如果该方向有体元属性为1的，则保留为1的体元，即视作该方向的体元属性为1，否则视作该方向的体元属性为0。

步骤4、间隙率(gap fraction,P)的计算：

在天顶方向以5°为间隔将0°到90°的天顶角分成18个区域，并以每个区域的中间天顶角值代表该区域的天顶角。同时在方位角方向以45°为间隔将0°到360°的方位角分成8个区域，并以每个区域的中间天顶角值代表该区域的方位角，得到144个扇形区域。通过统计各个扇形区域的总体元个数和属性为0的体元个数，得到扇形区域的间隙率为属性为0的体元个数与总体元个数之比，公式如下：

式中，θ为天顶角，为方位角。

步骤5、聚集度指数(clumping index,Ω)的计算：

由步骤4得到各个扇形区域的间隙率假设每个扇形区域都有间隙，即可通过以下公式求出各天顶方向的聚集度指数(Lang和Xiang，1986)：

式中θ为天顶角，为方位角，为冠层平均间隙率，为间隙率的对数平均。

本发明利用遥感技术手段(地面三维激光雷达扫描系统)，获取植被样方冠三维点云数据，通过构建三维体元模型、转换坐标系统、计算冠层间隙率等过程，建立基于地面高分辨率激光点云数据提取植被冠层的聚集度指数的方法。本发明适用于所有植被冠层，但是较之于针叶林，阔叶林利用该方法得到的结果会更加精确。

该方法可快速、准确的提取植被冠层聚集度指数，数据采集简便，且不受观测时光照条件的影响，在研究时也不用考虑卫星遥感数据或产品的像元尺度问题。并且采用的激光雷达技术手段不会对植被结构和辐射特性造成任何不良影响，同时还可以永久性的记录植被样地的三维结构特征，这将有利于进一步研究其它生物物理参数。此外，本发明通过构建三维体元模型，极大地减少了计算量。

综上所述，本发明简便，高效，对植被无副作用，且相对现有技术极大地减少了计算量。

附图说明

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

图2为数据获取的示意图；a.三维激光扫描仪现场工作图；b.样方的冠层点云侧视图；c.样方的冠层点云仰视图；

图3为构建三维体元模型的示意图；

图4为同一样方真实的数字半球摄影照片；

图5为玉兰树样方的地面激光雷达技术与数字半球摄影技术计算聚集度指数结果的比较分析。

具体实施方式

以下通过实例并结合附图对本发明作进一步说明：

步骤1，以玉兰树样方为研究对象(面积10m*10m，平均树高约7m),使用地面三维激光扫描仪Leica ScanStation C10(其参数如表1所示)在样方的中心位置及个侧面进行多站扫描，扫描仪离地高度为1米，扫描分辨率为高分辨率。在进行数据配准后，手动去除地面点云及其他噪声点云，得到玉兰树样方的冠层三维点云数据，如附图2示。

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

步骤2，在取得玉兰树样方的冠层三维点云数据并预处理后，利用体元化方法，构建三维体元模型。将体元大小设置为0.1m*0.1m*0.1m，通过判断体元内包含的激光点个数来确定各体元的体素值为1还是0。如果体元内激光点个数大于等于1，则将体元体素值赋为1，否则体素值赋为0。

步骤3，将体元化后的冠层三维点云数据从笛卡尔坐标系转为半径为1的球面坐标系统。同时计算每个体元的天顶角和方位角。

步骤4，在天顶方向以5°为间隔将0°到90°的天顶角分成18个区域，并以每个区域的中间天顶角值代表该区域的天顶角。同时在方位角方向以45°为间隔将0°到360°的方位角分成8个区域，并以每个区域的中间天顶角值代表该区域的方位角。这样便得到了144个扇形区域。通过统计各个扇形区域的总体元个数和属性为0的体元个数，利用式(2)计算得到各扇区的间隙率

步骤5，得到了各个扇形区域的间隙率后，通过式(3)得到各天顶方向的聚集度指数(Lang和Xiang，1986)(参见图5)。

综上可见，该实例根据本发明提出的方法，对玉兰树样方的激光雷达点云数据进行分析，依照技术方案所述，获取了样方冠层的聚集度指数。同时，在相同位置和高度采集同一样方真实的数字半球摄影照片(参见图4)，利用数字半球摄影技术计算得到该样方的聚集度指数。将利用实施例所用的激光雷达技术(LIDAR-based)和数字半球摄影技术(DHP-based)所得到的相同样方的聚集度指数的结果进行比较分析(参见图5)，可以看出，在天顶角为0°到65°之间，两种方法得到的聚集度指数的值是相近的。在天顶角为65°到90°之间，由于截取的冠层三维点云数据范围限制，该区域也没有点云数据。故这个范围的聚集度指数不考虑。当然两种方法得到的聚集度指数的值在天顶角为0°到65°之间也不是完全相同的，具有一定的差异性，这是因为利用数字半球摄影技术计算得到的聚集度指数本身具有光学测量上的误差，同时本发明所用的激光雷达技术也存在一定的误差，其中包括植被冠层三维点云数据预处理时期的配准误差、体元模型构建时期所选体元大小带来的误差等。不过本发明的误差可以通过改善实验手段、选取体元大小来减小。综上所述，本发明的方法是可行且有效的。

Claims (1)

1.一种三维激光点云提取植被冠层聚集度指数的方法，具体包括以下步骤：

步骤1、利用地面激光雷达扫描系统，获取植被冠层三维点云数据：

首先，在目标区样方中心点和样方外侧分别架设三维激光扫描仪，以样方中心点为观测点所得的坐标系统为标准，将得到的多站点云数据进行点云配准；

然后，在水平方向将样方的三维点云数据切割为以样方中心点为圆心，3≤r≤10m为半径的圆形区域；再将所有低于三维激光扫描仪高度的点剔除，得到植被冠层三维点云数据；

步骤2、三维体元模型构建：

依据步骤1得到的冠层三维点云数据求出笛卡尔坐标X、Y、Z的最小值(X_min，Y_min，Z_min)和最大值(X_max，Y_max，Z_max)，以X、Y、Z的(X_min，Y_min，Z_min)为起始点，以体元大小为步长划分冠层三维点云数据，并确定点云在体元坐标系中对应的体元坐标值及体素值；体元大小由体元的长L、宽W、高H决定，整个数据区域被划分为N_L×N_W×N_H个体元，其中，N_L＝(X_max-X_min)/L，N_W＝(Y_max-Y_min)/W，N_H＝(Z_max-Z_min)/H；点云体元化后的坐标值由以下公式得到：

式中，int是取整符，直接取出小数前面的整数部分，(i，j，k)是点云数据笛卡尔坐标(X，Y，Z)对应的体元坐标，体元大小L×W×H与扫描采用的点间距一致；

体元的体素值通过判断体元内包含的激光点个数来确定，如果体元内激光点个数大于等于1，代表激光束被体元拦截，则该体元体素值赋为1，否则体素值赋为0；

步骤3、坐标系统的转换：

将体元化后的冠层三维点云数据从笛卡尔坐标系转为半径为1的球面坐标系统，去除重复的体元，即确保转换后的球面坐标系统中每个点只有一个体元，如果该方向有体元属性为1的，则保留为1的体元，即视作该方向的体元属性为1，否则视作该方向的体元属性为0；

步骤4、间隙率的计算：

在天顶方向以5°为间隔将0°到90°的天顶角分成18个区域，并以每个区域的中间天顶角值代表该区域的天顶角；同时在方位角方向以45°为间隔将0°到360°的方位角分成8个区域，并以每个区域的中间天顶角值代表该区域的方位角，得到144个扇形区域；通过统计各个扇形区域的总体元个数和属性为0的体元个数，得到扇形区域的间隙率为属性为0的体元个数与总体元个数之比，公式如下：

式中，θ为天顶角，为方位角；

步骤5、聚集度指数Ω(θ)的计算：

由步骤4得到各个扇形区域的间隙率假设每个扇形区域都有间隙，即可通过以下公式求出各天顶方向的聚集度指数：

式中θ为天顶角，为方位角，为冠层平均间隙率，为间隙率的对数平均。