基于移动激光扫描的橡胶林风害相关参数精准反演

引用本文

云挺, 张艳侠, 王佳敏, 胡春华, 陈邦乾, 薛联凤, 陈凡迪. 基于移动激光扫描的橡胶林风害相关参数精准反演[J]. 光谱学与光谱分析, 2018,38(11): 3452-3463.
YUN Ting, ZHANG Yan-xia, WANG Jia-min, HU Chun-hua, CHEN Bang-qian, XUE Lian-feng, CHEN Fan-di. Quantitative Inversion for Wind Injury Assessment of Rubber Trees by Using Mobile Laser Scanning[J]. Spectroscopy and Spectral Analysis, 2018,38(11): 3452-3463.
Doi:10.3964/j.issn.1000-0593(2018)11-3452-12 复制到剪切板

Permissions

《光谱学与光谱分析》期刊社所有

基于移动激光扫描的橡胶林风害相关参数精准反演

云挺^1,³, 张艳侠¹, 王佳敏¹, 胡春华¹, 陈邦乾², 薛联凤^1,^*, 陈凡迪¹

1. 南京林业大学信息科学技术学院, 江苏南京 210037

2. 中国热带农业科学院橡胶研究所/农业部儋州热带作物野外观测实验站, 海南儋州 571737

3. 南京林业大学生物与环境学院, 江苏南京 210037

*通讯联系人 e-mail: xuelianfeng@njfu.edu.cn;云挺, 张艳侠, 王佳敏: 并列第一作者

作者简介: 云挺, 1980年生, 南京林业大学信息科学技术学院副教授 e-mail: njyunting@qq.com;张艳侠, 1994年生, 南京林业大学信息科学技术学院本科生 e-mail: 2569593898@qq.com;王佳敏, 1997年生, 南京林业大学信息科学技术学院本科生 e-mail: 15298393123@163.com

收稿日期: 2018-02-13

基金: 国家重点研发计划(2017YFD0600900), 国家自然科学基金项目(31770591, 41701510), 中国博士后面上基金项目(2016M601823), 江苏高校优势学科建设工程项目资助

摘要

林木参数反演是森林资源管理与培育经营的关键环节。迅速发展的激光探测与测量技术突破了传统测量方法, 可以快速的获取林木的空间三维信息, 在林业普查中具备高效率、高精度的优势。结合计算机图形学与图像学方法, 以中国最大的橡胶生产基地海南省儋州市长期受台风侵害下的两个不同品种橡胶林段(林段1 PR107, 林段2 CATAS7-20-59)为研究对象, 设计了面向离散激光点云的单株林木参数提取方法, 自动获取橡胶林木风害后的生物量指标。首先, 通过人背负移动激光雷达获取林段点云数据并使用瑞利商求取台风造成的主枝干倾角, 以找寻每株橡胶树的树冠中心点。其次, 对点云进行垂直投影, 并采用分水岭与Meanshift算法实现株株分离。最后, 基于以上操作自动获取与实测值相近的林木参数, 例如冠幅、冠径、冠积、叶面积密度、叶面积分布以及主枝与分枝之间夹角等。计算表明, 林段1与林段2东西冠幅分别为3.95和3.73 m, 与实测相差1.74%~6.27%; 林段1与林段2南北向冠幅分别为6.47和6.51 m, 与实测相差2.54%~4.02%; 林段1与林段2平均胸径分别为5.20和4.73 cm, 与实测相差0.64%~2.44%; 林段1与林段2平均冠积分别为168.01和141.80 m³, 与实测相差0.67%~0.85%; 林段1与林段2主枝干倾斜角分别为18.80°和13.11°, 与实测相差5.53%~7.09%; 林段1与林段2二级分支与主枝干的夹角分别为40.21°~69.23°和10.63°~32.14°, 它们相差62.63%; 林段1在不同天顶角下的叶面积指数均大于林段2。通过对一定样本(150棵/每类林段)的分析结果与真实比对表明, 该方法对林木参数反演结果精度较高, 有效地评估了不同品种橡胶树在台风下的损伤度(如主枝干歪斜率、叶面积密度及衰减分布)。参数反演结果与实测值仅有8%的偏差, 此偏差主要由橡胶林分叶片稠密, 导致林分中叶面与枝干扫描数据获取缺失, 以及外界环境干扰如风力扰动、点云拼接误差、激光束发散率、激光扫描范围等原因造成。同时, 由于林段1(PR107)的橡胶树的主枝与分枝夹角、冠积以及叶面积指数均大于林段2(CATAS7-20-59)的橡胶树, 导致林段1的橡胶树比林段的2橡胶树在台风侵害下更脆弱。因此, 该研究可用于研究风力侵害对于不同森林地块的影响, 以及量化风害造成的生态系统紊乱的影响。同时, 该方法解决了人背负激光扫描数据进行单株提取与林木参数反演的问题, 为激光测绘在林业中的应用提供了新思路。

关键词: 激光探测与测量技术; 单株树冠提取; 林木参数反演; 风害

中图分类号:TP79 文献标志码:A

Quantitative Inversion for Wind Injury Assessment of Rubber Trees by Using Mobile Laser Scanning

YUN Ting^1,³, ZHANG Yan-xia¹, WANG Jia-min¹, HU Chun-hua¹, CHEN Bang-qian², XUE Lian-feng^1,^*, CHEN Fan-di¹

1. College of Information Science and Technology, Nanjing Forestry University, Nanjing 210037, China

2. Rubber Research Institute/Danzhou Investigation & Experiment Station of Tropical Crops of Ministry of Agriculture, CATAS, Danzhou 571737, China;

3. College of Biology and the Environment, Nanjing Forestry University, Nanjing 210037, China

Abstract

The light detection and ranging (LiDAR) technique, which has the advantages of high efficiency and high accuracy in forest survey and is superior to the traditionalinformation acquisitionmethods, can be used to quickly obtain high-resolution mapping of morphological structures of forest. In this paper, two rubber forest plots (forest plot 1, clonePR107; forest plot 2, clone CATAS7-20-59) are taken as the study subjects, which are under the long-term impact of wind-induced disturbance severity and located in Danzhou city, the largest rubber production base of Hainan Island, China. First, point cloud of the forest plots through man-loaded mobile LiDAR was collectedand Ruili entropy method was designed to process the scanned data for calculating the slope angle of tree trunk (typhoon-induced) in order to find the canopy centre of each tree. Second, after the vertical projection of scanned forest points, Watershed and Mean shift algorithm were adopted to realize individual tree canopy delineation. Finally, many tree parameters, such as crown breadth, Diameter at Breast Height (DBH), crown volume, leaf area density, leaf distribution and included angle between trunk and main branches, were deduced automatically for analyzing the impact of typhoon disturbance on the two forest plots. Overall parameters obtained from our methods were compared with manual field measurements. The calculated average crown diameter in west-east direction of rubber trees in forest plot 1 and plot 2 using our method were 3.95 and 3.73 m, respectively, with false rate of 1.74% for forest plot 1 and 6.27% for plot 2. The calculated average crown diameter in north-south direction of rubber trees in forest plot 1 and plot 2 using our method were 6.47 and 6.51 m, respectively, with false rate of 4.02% for forest plot 1 and 2.54% for plot 2. The calculated average diameter at breast height (DBH) for forest plot 1 and plot 2 using our method were 5.20 and 4.73 cm, respectively, with false rate of 2.44% for forest plot 1 and 0.64% for plot 2. The calculated average crown volume for forest plot 1 and plot 2 using our method were 168.01 and 141.80 m³, respectively, with false rate of 0.67% for forest plot 1 and 0.85% for plot 2. The calculated average inclination angle of rubber trunk for forest plot 1 and plot 2 using our method were 18.80° and 13.11°, respectively, with false rate of 5.53% for forest plot 1 and 7.09% for plot 2. The calculated average included angle between trunk and branch for forest plot 1 ranged from 45.21° to 69.23°, and the calculated average included angle between trunk and branch for forest plot 2 ranged from 10.63° to 32.14°. Thedifference in the included angles of two forest plots was nearly 62.63%. Meanwhile, the leaf area index (LAI) of forest plot 1 derived fromhemispherical photos of various zenith angles was generally higher than forest plot 2. Compared with the in-situ measurements, the forest parameters from the subsample (scanned data of 150 trees per forest plot) were accurately retrieved using our method with a deviation of less than 8%. A variety of disturbance, such as the perspective occlusion caused by closed forest canopies, the error produced by multi-scan registration, vegetative elements moved by wind during the scanning process, beam divergence and scanning range constraint of the scanner, hampers the accurate scanned data acquisition and generates computer errors in the algorithm. Meanwhile, the included angle between trunks and branches, canopy volume and leaf area index of rubber tree clone PR107 (in forest plot 1) were overall higher than the parameters of rubber tree clone CATAS7-20-59 (in forest plot 2), resulting in the existence of higher vulnerability of clone PR107 than clone CATAS7-20-59 when wind damage propagation occurred in the forest plots. Thus, our research can be used to study the effects of wind disturbance on different forest plots and to quantify ecological instability of the forest in response to wind excitation. Our method makes a contribution to solving the problem of tree canopy delineation and forest parameter retrieval using man-loaded laser scanning technique, showing promise for further exploration of utilizing mobile terrestrial LiDAR as an effective tool for the applications in forest survey.

Keyword: Light detection and ranging (LiDAR); Individual tree canopy delineation; Forest parameter retrieval; Wind damage

文章图片

引言

森林对二氧化碳下降、动物群落、调节水文流动和巩固土壤起着重要的作用。若森林遭受破坏, 则会产生长期和不可逆的影响。因此, 如何对森林进行实时监测, 准确而又高效地获取森林相关资源信息对研究森林的生长状况及其动态变化具有非常重要的意义^[1]。

正是由于森林的重要性, 近年来, 国内外许多的学者做了大量的面向森林监测的相关工作。按研究区域的尺度可分为以下三类:

(1) 在宏观尺度下, 运用星载影像数据获取大尺度林地影像并结合图像处理算法实现单株定位与树冠参数提取, 如采用QuickBird^[2]、 WorldView-2^[3]、 Ikonos卫星数据、 Landsat 8多波段成像数据^[4]、 ALOS PALSAR微波波段成像数据^[5]、 Sentinel-2多光谱遥感数据^[6]等。虽然宏观尺度下可获取大范围林木的波段影像信息, 但大尺度下精度低、高光谱敏感波段难找寻, 以及树种、树龄、环境等因素的影响, 都会造成该测量方法对不同种类的树冠识别与冠中心定位不准确^[7]。

(2) 在中观尺度下, 如机载(light detection and ranging, LiDAR)^[8]、航空影像、移动式激光扫描、地基激光雷达^[9]、高光谱小型机载光谱成像仪等, 都应用于林业样地调查与参数反演中。然而, 树木自身或相互间的遮挡、扫描仪参数设置和林木结构差异性等因素都对单株树木提取及林木监测有着一定的影响。

(3) 在微观尺度下, 主要有: 落叶法、角规测树法、半球摄影法, 叶面积测定法和LAI-2200植物冠层分析仪等方法。该类方法较为简单, 无需昂贵仪器, 但人力成本高、效率低、易受外界因素(如季节、天气、光照和色差等)和采样频率与精度的影响。

基于上述情况, 以海南儋州的两个橡胶林段作为研究对象, 采用自下而上的人背负激光雷达的扫描方式获取林木点云数据。并结合图形图像学方法, 设计了基于人背负激光雷达点云的橡胶林木株株分离与参数反演算法。同时考虑到橡胶林木长期受到台风胁迫出现歪倒现象, 以歪斜的主枝干作为主要特征进行单株定位, 并获取对应的林木参数, 以研究不同橡胶树品种在台风侵害下林木损伤度评估与成因。整个项目的实验流程如图1所示。

	Figure Option View Download New Window
	图1 实验主要步骤流程图Fig.1 Flowchart of experimental process

1 研究区概况与数据获取

1.1 研究区概况

研究区选在海南儋州的橡胶实验基地, 地处海南岛的西部工业走廊北段(19° 32'47.89″N, 109° 28'29.33″E)。儋州市是海南岛最大的橡胶生产基地, 橡胶林众多, 橡胶资源极其丰富。儋州市属于热带季风海洋性气候区, 夏无酷暑, 冬无严寒。阳光充足, 气候温和, 年均气温23.2 ℃, 极端最高气温为40.0 ℃。雨量充沛但分布不均, 冬春雨量较稀少, 夏秋雨水较充足, 年均降雨量1 815 mm, 5月至10月份为雨季, 11月至翌年4月份为旱季。强台风平均每4年一遇。受热带风暴及台风影响, 部分树木已出现弯折, 且橡胶树的树冠顶并未出现形态结构明显的冠形特征。因此如何精确有效地实现强风力载荷后橡胶林分精准参数的反演较为困难。研究区概况如图2所示。

	Figure Option View Download New Window
	图2 中国海南岛儋州市实验农场(a)和用于研究的橡胶树(b)Fig.2 Study area map (a) and rubber trees (b) at CATAS experimental farm, Danzhou city, Hainan Island, China

1.2 激光雷达数据的获取

使用的原始点云数据获取装置为Velodyne HDL-32E激光雷达传感器, 其创新的激光阵列可使系统观测到比其他激光雷达传感器更多的信息。相关参数: 重量小于2 kg, 操作温度: -10~+60℃, 输入电压为9~32V, 垂直视野为+10.67° ~-30.67° , 水平视野可达360° 。其发出的激光束波段为红外波段(896~910 nm), 为排除环境噪声干扰, 扫描器使用波长905 nm, 脉冲宽度为5 ns的激光进行扫描。人背负仪器, 按照“ 之” 字型在林段内来回走动进行测量。雷达以每秒10 Hz的频率进行水平旋转并发射波长905 nm的激光, 经物体散射返回, 被接收器探测后, 每秒可产生约70万个点, 可获得探测距离为80~100 m的实时数据。对于整个系统数据的拼接采用的是simultaneous localization and mapping(SLAM)算法^[10]。此外该系统还具有重量轻, 体积小, 结构坚固, 便于携带, 高分辨率、高精确, 且选用的工作波长不受阳光的影响, 有良好的穿透烟、雾等的特性优势。数据获取时间为2017年7月30日。实时作业的温度为+26~+34 ℃, 获取到两个林段点云数据, 其中林段1橡胶树种(PR107, 台风侵害下易伏倒), 林段2橡胶树种(CATAS7-20-59, 台风侵害下不易伏倒)。选取部分林段点云数据, 分别共有1 759 870(林段1)和1 326 925(林段2)个扫描点。

2 研究方法

2.1 消除地面高程不一致

获取到的激光雷达数据, 由于林地起伏, 其地面高程不一致, 会对树干检测造成误差。因此, 设置了平滑窗口, 对从林段1选取的N₁(N₁≈ 150)棵树、林段2选取的N₂(N₂≈ 150)棵树的原始点云数据进行平滑滤波处理, 以消除实验数据中的地面高程差异性。设一窗口下的点云数据为: p_i(x_i, y_i, z_i)。令另一窗口下点云数据中的任一点集为p_j(x_j, y_j, z_j), 为使这两个窗口中的点云数据处于同一地面水平上, 令diff=min(z_j)-min(z_i), 则z_j=z_j-diff。随着滑动窗口连续平移, 不断消除地形高度差异性, 进而使得整片林分处于相同的地形高度上。如图3(a)和(b)所示。

	Figure Option View Download New Window
	图3 (a)林段1(PR107)未消除地面高度差的原始扫描点云扫描图; (b)消除地面高度差后的扫描点云图Fig.3 (a) The initial scanned points (forest segment 1, clone PR107); (b) The processed scanned data after elimination of the altitude difference

2.2 林分枝叶分离及主枝干提取

激光扫描后, 获取到的仅是稠密的林木点云数据, 并非是橡胶林木的具体器官信息, 如单株枝干的空间信息、树冠的叶子点云等。因此, 采用文献[9]的方法对林木点云进行计算, 实现枝叶分离。

通过该算法, 获取每个激光点p_i的特征, 记为: $f_{p_{i}}$ ={e_ix, e_iy, e_iz, c_0,_i, c_1,_i, c_2,_i, l_0,_i, l_1,_i, l_2,_i}, 其中: {e_ix, e_iy, e_iz}为点p_i的法向量, {c_0,_i, c_1,_i, c_2,_i}为点p_i的结构张量特征, 而{l_0,_i, l_1,_i, l_2,_i}为点p_i的形状特征值。

根据获得的每个激光扫描点p_i的特征 $f_{p_{i}}$ , 结合高斯分类器实现对橡胶林段的枝叶分离操作。按照橡胶树自然生长规律, 以树木总高度的1/2以下区域作为主枝干数据(其余部分为分支干数据), 从而获取到两个林段每棵树的主枝干信息。如图4(a)和(b)所示。

2.3 冠中心定位与主干倾斜度分析

根据2.2节已获取到每株橡胶树的主枝干信息, 发现大部分的橡胶树的树干倾斜, 由此造成冠中心位置与主枝干天顶方向有偏差, 导致无法直接获取树冠的中心点。

因此对每棵树树冠中心进行如计算下: 选取林段中第k棵树, 以根部往上取相对于整棵树高度h^k的下部1/20区域的点云数据( $p_{1}^{k}$ , …, $p_{m}^{k}$ )为基点, 拟合出m条沿树干方向的最佳拟合直线, 其中m为整棵树高度h^k的下部1/20区域的所有点云数据的个数, 拟合直线的点向式表现为: $L_{i}^{k}$ ( ${\overset{⇀}{v}}_{i}^{k}$ ( $v_{i, x}^{k}$ , $v_{i, y}^{k}$ , $v_{i, z}^{k}$ ), $p_{i}^{k}$ ( $p_{i, x}^{k}$ , $p_{i, y}^{k}$ , $p_{i, z}^{k}$ )), 其中 ${\overset{⇀}{v}}_{i}^{k}$ 为直线的方向向量(未知), $p_{i}^{k}$ 为点云( $p_{1}^{k}$ , …, $p_{m}^{k}$ )中的一点, i=1, …, m。为求最佳拟合直线, 即求解 ${\overset{⇀}{v}}_{i}^{k}$ , 采用瑞利商的方法。由于最佳拟合直线即为树干的所有点云数据到直线 $L_{i}^{k}$ 的正交距离的平方和最小

$\begin{array}{l} argmin ({\overset{⇀}{v}}_{i}^{k} (v_{i, x}^{k}, v_{i, y}^{k}, v_{i, z}^{k})) = \\ \overset{S}{\sum_{j = 1}} \frac{| p_{i}^{k} - p_{j}^{k} |^{2} | {\overset{⇀}{v}}_{i}^{k} |^{2} - ((p_{i}^{k} - p_{j}^{k}) {\overset{⇀}{v}}_{i}^{k})^{2}}{| {\overset{⇀}{v}}_{i}^{k} |^{2}} (1) \end{array}$

其中S为整个主枝干0~h_trunk高度的点云个数(h_trunk为树干基底到树冠下部的高度[如图4(a)所示)], $p_{j}^{k}$ 为橡胶树整个主枝干部分S个点云数据中的任意一点。令w_ij= $p_{i}^{k}$ - $p_{j}^{k}$ , $c_{ij}^{2}$ =|w_ij|², 其中w_ij(w_{ij, x}, w_{ij, y}, w_{ij, z})= $p_{i}^{k}$ ( $p_{i, x}^{k}$ , $p_{i, y}^{k}$ , $p_{i, z}^{k}$ )- $p_{j}^{k}$ ( $p_{j, x}^{k}$ , $p_{j, y}^{k}$ , $p_{j, z}^{k}$ ), 所以式(1)可改写为

$\begin{array}{l} argmin ({\overset{⇀}{v}}_{i}^{k} (v_{i, x}^{k}, v_{i, y}^{k}, v_{i, z}^{k})) = \\ \frac{\overset{S}{\sum_{j = 1}} c_{ij}^{2} ({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k} - (v_{i, x}^{k} w_{ij, x} + v_{i, y}^{k} w_{ij, y} + v_{i, z}^{k} w_{ij, z})^{2}}{({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k}} \\ = \frac{α ({\overset{⇀}{v}}_{i}^{k}) - β ({\overset{⇀}{v}}_{i}^{k})}{({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k}} (2) \end{array}$

	Figure Option View Download New Window
	图4 枝叶分离结果 (a): 利用本方法对林段1(PR107)中典型橡胶树的扫描点自动分类为叶子(绿色)和枝干(红色); (b): 林段2(CATAS 7-20-59)中典型橡胶树的枝叶分离结果图Fig.4 Results of tree organ classification (a): Leaves (green) and branches (red) are separated for clone PR107; (b): Separated result for clone CATAS 7-20-59

其中:

$\begin{array}{l} α ({\overset{⇀}{v}}_{i}^{k}) = ({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k} \overset{S}{\sum_{j = 1}} c_{ij}^{2} \\ β ({\overset{⇀}{v}}_{i}^{k}) = (v_{i, x}^{k})^{2} \overset{S}{\sum_{j = 1}} w_{ij, x}^{2} + (v_{i, y}^{k})^{2} \overset{S}{\sum_{j = 1}} w_{ij, y}^{2} + \\ (v_{i, z}^{k})^{2} \overset{S}{\sum_{j = 1}} w_{ij, z}^{2} + 2 v_{i, x}^{k} v_{i, y}^{k} \overset{S}{\sum_{j = 1}} w_{ij, x} w_{ij, y} + \\ 2 v_{i, x}^{k} v_{i, z}^{k} \overset{S}{\sum_{j = 1}} w_{ij, x} w_{ij, z} + 2 v_{i, y}^{k} v_{i, z}^{k} \overset{S}{\sum_{j = 1}} w_{ij, y} w_{ij, z} \end{array}$

设矩阵B为点集( $p_{1}^{k}$ , …, $p_{m}^{k}$ )的协方差矩阵乘S(S为整个主枝干0~h_trunk高度的点云个数)

$B = U^{T} U U = (\begin{array}{l} {w^{T}}_{i 1} \\ ︙ \\ {w^{T}}_{iS} \end{array}) (3)$

其中U是S× 3的矩阵, B是3× 3的矩阵, w_ij是3× 1的矩阵。则α 和β 可以表示为

$β ({\overset{⇀}{v}}_{i}^{k}) = ({\overset{⇀}{v}}_{i}^{k})^{T} B {\overset{⇀}{v}}_{i}^{k}, α ({\overset{⇀}{v}}_{i}^{k}) = ({\overset{⇀}{v}}_{i}^{k})^{T} (el) {\overset{⇀}{v}}_{i}^{k}$

其中e= $\overset{S}{\sum_{j = 1}} c_{ij}^{2}$ , l是3× 3的单位矩阵。则式(2)则可以写为

$\begin{array}{l} argmin ({\overset{⇀}{v}}_{i}^{k} (v_{i, x}^{k}, v_{i, y}^{k}, v_{i, z}^{k})) = \frac{({\overset{⇀}{v}}_{i}^{k})^{T} (el) {\overset{⇀}{v}}_{i}^{k} - ({\overset{⇀}{v}}_{i}^{k})^{T} B {\overset{⇀}{v}}_{i}^{k}}{({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k}} \\ = \frac{({\overset{⇀}{v}}_{i}^{k})^{T} (el - B) {\overset{⇀}{v}}_{i}^{k}}{({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k}} = \frac{({\overset{⇀}{v}}_{i}^{k})^{T} M {\overset{⇀}{v}}_{i}^{k}}{({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k}} (4) \end{array}$

其中M=el-B。瑞利商 $\frac{({\overset{⇀}{v}}_{i}^{k})^{T} M {\overset{⇀}{v}}_{i}^{k}}{({\overset{⇀}{v}}_{i}^{k})^{T} {\overset{⇀}{v}}_{i}^{k}}$ 中, 当 ${\overset{⇀}{v}}_{i}^{k}$ 等于矩阵M的最小特征值对应的特征向量时, 其表达式的值达到最小值。根据特征值分解得到

$Mξ = (el - B) ξ = emξ - Bξ = eξ - λξ = (e - λ) ξ (5)$

其中λ 为矩阵B的任意特征值, ξ 为对应的特征向量。

由式(5)可得, 矩阵B的特征值都与矩阵M的特征值相关联, 当取到矩阵M的最小特征值对应的特征向量时, 即为矩阵B的最大特征值对应的特征向量。因此, 矩阵B的最大特征值相对应的特征向量 ${\overset{⇀}{v}}_{i}^{k}$ ( $v_{i, x}^{k}$ , $v_{i, y}^{k}$ , $v_{i, z}^{k}$ ), 作为整个主枝干0~h_trunk高度部分S个点到拟合直线正交距离最小所对应的拟合直线的方向向量。则第i个橡胶树主枝干点所对应的最佳拟合直线为

$\frac{x - p_{i, x}^{k}}{v_{i, x}^{k}} = \frac{y - p_{i, y}^{k}}{v_{i, y}^{k}} = \frac{z - p_{i, z}^{k}}{v_{i, z}^{k}} (6)$

反复选取不同的 $p_{i}^{k}$ , 迭代求解多个 ${\overset{⇀}{v}}_{i}^{k}$ ( $v_{i, x}^{k}$ , $v_{i, y}^{k}$ , $v_{i, z}^{k}$ ), 并获取平均值 $\bar{{\overset{⇀}{v}}^{k}}$ ( $\bar{v_{x}^{k}}$ , $\bar{v_{y}^{k}}$ , $\bar{v_{z}^{k}}$ ), 同时求解所选取所有的 $p_{i}^{k}$ 的平均点 $\bar{p^{k}}$ ( $\bar{x_{avg}^{k}}$ , $\bar{y_{avg}^{k}}$ , $\bar{{z^{k}}_{avg}}$ ), 由此得到每株树干的最佳拟合点向式直线方程 $\bar{L^{k}}$ , 如式(7)

$\frac{x - \bar{x_{avg}^{k}}}{\bar{v_{x}^{k}}} = \frac{y - \bar{y_{avg}^{k}}}{\bar{v_{y}^{k}}} = \frac{z - \bar{z_{avg}^{k}}}{\bar{v_{z}^{k}}} (7)$

同时根据经验推断, 得平均直线 $\bar{L^{k}}$ 在高度 $\frac{5}{3}$ h_trunk下的交点即为冠中心点, 其中h_trunk为树干基底到树冠下部的高度[如图5(a)所示], 带入直线方程, 由此可求第k棵数的冠中心点坐标为 $P_{cc}^{k}$ ( $x_{cc}^{k}$ , $y_{cc}^{k}$ , $z_{cc}^{k}$ )

$\begin{array}{l} (\frac{(\frac{5}{3} h_{trunk} - \bar{z_{avg}^{k}}) \bar{v_{x}^{k}}}{\bar{v_{z}^{k}}} + \bar{x_{avg}^{k}}, \frac{(\frac{5}{3} h_{trunk} - \bar{z_{avg}^{k}}) \bar{v_{y}^{k}}}{\bar{v_{z}^{k}}} + \\ \bar{y_{avg}^{k}}, \frac{5}{3} h_{trunk}) \end{array}$

Figure Option
View Download New Window

图5 枝干倾斜度与冠中心点计算示意图
(a): 树冠中心点获取以及结果示意图(左图为单株点云侧视图, 右图为冠中点计算示意图); (b): 立体图显示林段1(PR107)计算的树冠中心点及主枝干信息; (c): 立体图显示林段2(CATAS 7-20-59)计算的树冠中心点及主枝干信息Fig.5 Schematic diagrams of calculation for the oblique angle of the trunk and the central point of canopy
(a): Rubber tree (left) and calculation (right); (b): 3D image for clone PR107; (c): 3D image for clone CATAS-7-20-59

图5(a)即为冠中心点的获取过程示意图。

联合拟合直线的平均法向量 $\bar{{\overset{⇀}{v}}^{k}}$ 与天顶方向向量 $\vec{v_{2}}$ (0, 0, 1)(z轴方向), 得到每树的主枝干倾斜度, 其中k代表第几棵树。

$α^{k} = \arccos \frac{\bar{{\overset{⇀}{v}}^{k}} \cdot \vec{v_{2}}}{| \bar{{\overset{⇀}{v}}^{k}} \vec{v_{2}} |} = \arccos \frac{\bar{v_{z}^{k}}}{\sqrt[]{{\bar{v_{x}^{k}}}^{2} + {\bar{v_{y}^{k}}}^{2} + {\bar{v_{z}^{k}}}^{2}}} (8)$

图5(b)和(c)为针对不同林段数据样本处理得到的主干信息与倾斜度显示。

2.4 基于Meanshift和分水岭的株株分离

根据2.3节计算得到的单株橡胶树主干信息和树冠的中心点信息, 结合分水岭与Meanshift算法^{[11, 12]}, 对橡胶林段进行单株提取, 实现橡胶林的株株分离。

根据冠中心点坐标, 结合分水岭算法, 将该三维坐标投影到二维灰度图上, 得到点云数据的聚类图。如图6(a)和(b)中左图所示, 其中白色标记为每棵树的树冠中心。再结合分水岭与Meanshift算法, 实现单株树的空间分离。分类结果如图6(a)和(b)中右图所示, 其中林段1(PR107)共155棵橡胶树, 林段2(CATAS 7-20-59)共142棵橡胶树。各林段对应俯视图、侧视图如图6(c)和(d)所示。

Figure Option
View Download New Window

图6 单株树冠提取结果图
(a): 树冠中心点俯视图(林段1PR107)并结合分水岭与MeanShift算法实现株株分离; (b): 林段2(CATAS 7-20-59)的等价图; (c): 林段1(左图)、 2(右图)单株提取俯视图; (d)林段1(上图)、 2(下图)单株提取侧视图Fig.6 Delineation results of individual tree canopy
(a): Top view of the central points of canopies for clone PR107; (b): Top view of the central points of canopies for clone CATAS 7-20-59; (c): Top views of individul tree canopys of segment 1 (left) and segment 2(right); (d): Lateral views of segment 1 (upper) and segment 2 (lower)

3 结果与讨论

实现了对单株橡胶树的分离和提取, 同时也计算得到了主枝干最佳拟合直线, 对林段中每棵橡胶树分别取整棵树高度h的下部1/20区域点云的平均点 $\overset{̅}{p}$ (* 表示)以及最佳拟合直线在高度(2/3)h_trunk(■表示)、 (5/3)h_trunk(★表示)处交点[该参数的侧视图参见图5(a)], 整片林分的俯视图如图7所示, 其中三个标记点的重合度越高, 表示该橡胶树的倾斜程度越低。由此可以明确两个林段中每株橡胶树在台风侵害下的倾斜程度。

	Figure Option View Download New Window
	图7 不同种类橡胶树主枝干倾斜程度计算示意图 (a): 林段1 (PR107); (b): 林段2 (CATAS 7-20-59)Fig.7 Inclination level calculation of rubber tree trunles with different varieties (a): Forest segment 1 (clone PR107); (b): Forest segment 2 (clone CATAS 7-20-59)

3.1 林木参量计算值与实测值对比

从林段1和林段2中分别随机均匀选取十五棵树, 基于Convex Hull^[13]算法求得最小树冠凸包体积, 将凸包体积的平均值, 作为两个林段每株树的平均冠体积。另外, 对橡胶树的冠幅、倾斜度、平均胸径和平均冠体积等进行了参数分析, 两个林段的各参数准确度的定量评估如表1所示。并与人工实测数据相比较, 其相关度也罗列在表1中, 其中α 和β 所表示的含义见图5(a)所示。

表1 两个林分的算法参数估测值与人工测量值比对 Table 1 Accuracy of the retrieved parameters for two forest plots using our method compared with in-situ measurements

3.2 不同林段抗风性能分析

由于树木的冠幅、骨架拓扑结构、叶面空间分布情况的差异性都会对树木的抗风性能产生较大的影响, 因此本节将分析林段1(PR107)与林段2(CATAS 7-20-59)在这些参量上的差异性对其抗风性能的影响。

首先分析树木冠幅以及骨架拓扑结构差异性对于树木抗风性能的影响。本文选取两林段在东西和南北方向上的冠幅、主枝干倾斜度α 以及分支与主枝干间的夹角β [图5(a)]进行对比分析, 其对比箱体图如图8所示。表1表明, 林段1与林段2的冠幅相差不大, 林段1的树高明显低于林段2而平均冠积却大于林段2。同时, 林段1分支与主枝干间的夹角β 远大于林段2, 这些都会造成林段1的迎风面积大于林段2, 且林段1的冠内风力扰动高于林段2。

Figure Option
View Download New Window

图8 不同种类橡胶树冠幅以及分枝与主枝夹角对比图
(a): 两类橡胶树在不同方向上的冠幅比对箱线图; (b): 两类橡胶树(PR107和CATAS 7-20-59)的主枝干倾斜角α 和分支与主枝干间的夹角β 分布箱线图Fig.8 Comparison of canopy breadths, inclination angles (α ) of trunks, angles (β ) between trunks and branches
(a): Canopy breadths in different directions; (b): Angle α and angle β

其次分析树木叶密度对其抗风性能的影响, 本文从林段1和林段2中各选取30棵橡胶树, 分别求得其叶子点云数据分布的雷达图, 结果如图9所示, 其中每一个虚线标记的闭合区域代表对应单株树木的扫描点云在空间方向上的分布。由图发现, 两林段叶面点主要分布在东西方向上。当风力方向如图9所示时, 由于林段1的平均冠积大于林段2, 且叶片点云数量高于林段2, 使得林段1的受风面积增加。由此可知相同台风侵害下, 林段1受风力侵害程度较林段2严重, 即林段1(PR107)中橡胶树比林段2(CATAS 7-20-59)中橡胶树更易倒伏。

	Figure Option View Download New Window
	图9 两林段中叶子点云数据空间分布雷达图 (a): 林段1(PR107); (b): 林段2(CATAS 7-20-59)Fig.9 Radar charts showing the spatial distribution of leaf points of each rubber tree in two forest segments (a): Clone PR107; (b): Clone CATAS 7-20-59

最后, 为验证不同林段叶密度比对结果, 本文运用半球摄影获取两林段内的鱼目图像, 在不同天顶角下计算并获取叶面积指数, 结果如图10所示。图10(c)表明, 林段1的叶面积指数都高于林段2的叶面积指数。再结合其他树木指标, 可以验证林段1受风力侵害程度比林段2要大, 即林段1的橡胶树比林段2的橡胶树易倒伏。

	Figure Option View Download New Window
	图10 利用鱼目图像反演不同林段的叶面积指数 (a): 林段1(PR107); (b): 林段2(CATAS 7-20-59); (c): 不同天顶角获取的叶面积指数Fig.10 Using the hemispherical photography to retrieve leaf area index of different forest segments (a): PR107; (b): CATAS 7-20-59; (c): Retrieval leaf area index from different annulus sectors (representing different zenith angles)

4.3 林段内风力破坏定位

本文进一步对林段遭受风力破坏的区域进行定位, 通过对林段1与林段2的叶面积密度分析, 结合单株分离结果, 将叶面积密度相对稀疏的橡胶树体标记为红色(如图11所示)。

通过图示发现, 南北走向的风力遭受橡胶林阻遏后在红色区域进行突破, 穿越橡胶林分并带来区域林分损伤。

Figure Option
View Download New Window

图11 基于点云数据计算叶面积密度以确定林段中受风力损伤的区域
(a): 林段, 其中红色代表叶面积密度较低区域(台风导致); (b): 林段2Fig.11 Calculating the leaf area density from the scanned points to estimate the regions damaged by wind
(a): Top view of the analyzed result of forest segment 1, where red color represents the tree canopies with lower leaf area density (wind-induced); (b): Result for forest segment 2

5 结论

使用人背负移动LiDAR, 采用自下而上的移动扫描方式获取两个橡胶林段点云数据, 并设计图形图像学算法实现对橡胶林段的单株提取与参数反演。通过对林段1(PR107)和林段2(CATAS 7-20-59)中橡胶树参数的对比分析发现, 其林木中的生理参数的不同导致林段1受风力侵害程度要大于林段2, 也更易倒伏。

由于橡胶树树木高大、树冠郁闭度大、吸水能力强等特性, 使得其种植区域内的灌木与草本层植被生长受限, 林下缺乏中间的灌木层和地面植被。而这一生态特点也使得本文采用的地面移动激光扫描方式能够更加全面并无遮挡的获取到橡胶林木自下而上的点云数据。而对于林内中灌木层植被较为丰富的研究区域, 该扫描方式会受到遮挡影响并致使树木下部枝干区域扫描数据的缺失。同时, 本文后续工作将运用移动地基与机载LiDAR相结合的扫描方式, 开展更多品种橡胶树林分参数的精准反演, 为橡胶树的栽培种植与防风营造提供有效的理论依据。

The authors have declared that no competing interests exist.

参考文献

文献列表

[1]	XU Wei-heng, FENG Zhong-ke, SU Zhi-fang, et al(徐伟恒, 冯仲科, 苏志芳, 等). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2014, 34(2): 465. [本文引用:1]
[2]	Ardila J P, Bijker W, Tolpekin V A, et al. International Journal of Applied Earth Observation and Geoinformation, 2012, 15(4): 57. [本文引用:1]
[3]	Karlson M, Ostwald M, Reese H, et al. International Journal of Applied Earth Observation and Geoinformation, 2016, 50: 80. [本文引用:1]
[4]	Halperin J, Lemay V, Coops N, et al. Remote Sensing of Environment, 2016, 179: 170. [本文引用:1]
[5]	Cartus O, Kellndorfer J, Rombach M, et al. Remote Sensing, 2012, 44(11): 3320. [本文引用:1]
[6]	Castillo J A A, Apan A A, Maraseni T N, et al. ISPRS Journal of Photogrammetry and Remote Sensing, 2017, 134: 70. [本文引用:1]
[7]	Ploton P, Barbier N, Couteron P, et al. Remote Sensing of Environment, 2018, 200: 140. [本文引用:1]
[8]	Cao L, Coops N C, Innes J L, et al. Remote Sensing of Environment, 2016, 178: 158. [本文引用:1]
[9]	Yun T, An F, Li W, et al. Remote Sensing, 2016, 8(11): 942. [本文引用:2]
[10]	Pierzchała M, Giguère P, Astrup R. Computers and Electronics in Agriculture, 2018, 145: 217. [本文引用:1]
[11]	Jing L, Hu B, Noland T, et al. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 70(3): 88. [本文引用:1]
[12]	Hu B, Li J, Jing L, et al. International Journal of Applied Earth Observation and Geoinformation, 2014, 26(26): 145. [本文引用:1]
[13]	Cheein F A A, Guivant J. Computers and Electronics in Agriculture, 2014, 102(1): 19. [本文引用:1]

2014

0.0

... 因此, 如何对森林进行实时监测, 准确而又高效地获取森林相关资源信息对研究森林的生长状况及其动态变化具有非常重要的意义^[1] ...

2012

0.0

... (1) 在宏观尺度下, 运用星载影像数据获取大尺度林地影像并结合图像处理算法实现单株定位与树冠参数提取, 如采用QuickBird^[2]、 WorldView-2^[3]、 Ikonos卫星数据、 Landsat 8多波段成像数据^[4]、 ALOS PALSAR微波波段成像数据^[5]、 Sentinel-2多光谱遥感数据^[6]等 ...

2016

0.0

2016

0.0

2012

0.0

2017

0.0

2018

0.0

... 虽然宏观尺度下可获取大范围林木的波段影像信息, 但大尺度下精度低、高光谱敏感波段难找寻, 以及树种、树龄、环境等因素的影响, 都会造成该测量方法对不同种类的树冠识别与冠中心定位不准确^[7] ...

2016

0.0

... (2) 在中观尺度下, 如机载(light detection and ranging, LiDAR)^[8]、航空影像、移动式激光扫描、地基激光雷达^[9]、高光谱小型机载光谱成像仪等, 都应用于林业样地调查与参数反演中 ...

2016

0.0

... 因此, 采用文献[9]的方法对林木点云进行计算, 实现枝叶分离 ...

2018

0.0

... 对于整个系统数据的拼接采用的是simultaneous localization and mapping(SLAM)算法^[10] ...

2012

0.0

... 3节计算得到的单株橡胶树主干信息和树冠的中心点信息, 结合分水岭与Meanshift算法^[11,12], 对橡胶林段进行单株提取, 实现橡胶林的株株分离 ...

2014

0.0

... 3节计算得到的单株橡胶树主干信息和树冠的中心点信息, 结合分水岭与Meanshift算法^[11,12], 对橡胶林段进行单株提取, 实现橡胶林的株株分离 ...

2014

0.0

... 1 林木参量计算值与实测值对比从林段1和林段2中分别随机均匀选取十五棵树, 基于Convex Hull^[13]算法求得最小树冠凸包体积, 将凸包体积的平均值, 作为两个林段每株树的平均冠体积 ...