CN109446691B - 基于激光点云与空气动力学的活立木抗风性能分析方法 - Google Patents

基于激光点云与空气动力学的活立木抗风性能分析方法 Download PDF

Info

Publication number
CN109446691B
CN109446691B CN201811322277.0A CN201811322277A CN109446691B CN 109446691 B CN109446691 B CN 109446691B CN 201811322277 A CN201811322277 A CN 201811322277A CN 109446691 B CN109446691 B CN 109446691B
Authority
CN
China
Prior art keywords
branch
point cloud
point
tree
branches
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201811322277.0A
Other languages
English (en)
Other versions
CN109446691A (zh
Inventor
云挺
黄笑
黄智先
曹林
薛联凤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing Maoting Information Technology Co ltd
Original Assignee
Nanjing Forestry University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing Forestry University filed Critical Nanjing Forestry University
Priority to CN201811322277.0A priority Critical patent/CN109446691B/zh
Publication of CN109446691A publication Critical patent/CN109446691A/zh
Application granted granted Critical
Publication of CN109446691B publication Critical patent/CN109446691B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/06—Power analysis or power optimisation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Abstract

本发明公开了一种基于激光点云与空气动力学的活立木抗风性能分析方法，包括获取林木点云数据并枝叶分离；采用拉普拉斯算法对枝干点云进行收缩；将枝干点云数据自下而上切分为不同的层；求取每个高度分层的聚类中心点，根据聚类中心点拟合每一个高度分层的枝干；将活立木不同枝干骨架分类为主枝干和次级枝干；完成对活立木叶子点云数据的归属；建立林分模型，将林分模型加载风力，根据湍流模型和流固耦合模型分析林分内部动态压力、风速度以及湍流动能强度分布，本发明计算复杂度适中，能更好地描述活立木的空间结构特征与生长参数变化，实现活立木在台风下抗风性能的定性定量化评估，准确度高，为树木的栽培种植与防风营造提供准确的理论依据。

Description

基于激光点云与空气动力学的活立木抗风性能分析方法
技术领域
本发明涉及活立木抗风性能分析技术领域，具体涉及一种基于激光点云与空气动力学的活立木抗风性能分析方法。
背景技术
中国是世界上受到台风影响最大的国家之一,强台风不但会造成人员的伤亡、摧毁建筑物和桥梁外,还会导致树木的断裂倒伏。在海南省,每次台风登陆都会导致橡胶树林出现大量的根倒、弯曲和干折现象,可能造成数十亿元的经济损失,因此开展橡胶树抗风的研究具有重要的意义。
现今，对树木抗风性研究的主要方法分为：
1.从林分角度讲，如：通过将林分看成的多孔介质模型，结合大涡模拟方法，研究大气稳定度和树冠阻力对风力的综合作用。该方法虽然在宏观尺度下可快速模拟计算大范围内不同林分在台风下的空气流场参数，但准确度低，很难描述林分内精细参数的变化以及风力场气动耦合之间的变迁。
2.在微观尺度下，进行单株活立木的精细建模，例如一些学者采样一定的树叶和树枝并结合空气动力学，研究其在风力作用下的摆动。该方法希望将树木仿真重建与抗风性研究相结合。然而，由于树木自身枝干树叶的数量繁杂及形态多样，很难精细重建单棵树中每片树叶与其对应的枝干。因此，加载风力后，不能准确分析各分支和各树叶受风力后形变的运动状态，同时难以研究树叶受不同攻角风力对树枝的作用力，导致此方法效率低，复杂度高，难以实现。
发明内容
本发明所要解决的技术问题是针对上述现有技术的不足提供一种基于激光点云与空气动力学的活立木抗风性能分析方法，本基于激光点云与空气动力学的活立木抗风性能分析方法计算复杂度适中，能更好地描述活立木的空间结构特征与生长参数变化，通过计算机仿真技术，实现活立木在台风下抗风性能的定性定量化评估，准确度高，为树木的栽培种植与防风营造提供准确的理论依据。
为实现上述技术目的，本发明采取的技术方案为：
一种基于激光点云与空气动力学的活立木抗风性能分析方法，包括：
步骤1、枝叶分离：获取林木点云数据，对林木点云数据进行计算，实现枝叶分离；
步骤2、去噪操作：对获得的单株活立木的n个枝干点云进行去除噪声点操作；
步骤3、枝干点云骨架收缩：采用拉普拉斯算法对枝干点云进行收缩；
步骤4、枝干高度分层与聚类中心点求取：根据单株活立木的垂直高度将收缩后的枝干点云数据自下而上切分为不同的层；根据欧式距离进行聚类，设定距离阈值，每一个高度分层中的枝干点云根据距离阈值进行聚类归属，即将每一个高度分层的枝干点云之间距离少于距离阈值的点云集归属为一类；求取每类的中心点，即每一个高度分层的局部枝干的聚类中心点；
步骤5、圆柱体拟合与骨架重建：遍历每一个高度分层的聚类中心点，根据聚类中心点并采用带有空间方向性的圆柱体进行拟合每一个高度分层的枝干骨架，从而实现单株活立木的完整骨架的重建；
步骤6、枝干骨架分类：通过遍历聚类中心点中的根节点、父节点、分支结点、边缘节点和儿子节点将整株活立木不同枝干骨架分类为主枝干和次级枝干；
步骤7、叶子点云归属：基于分类好的枝干骨架并结合空间分水岭算法完成对活立木叶子点云数据的归属；
步骤8、孔隙率求取：计算活立木的孔隙率；
步骤9、叶片归属度分析：根据主枝干和次级枝干的数量将活立木的枝干分为多类，统计各类分支的点云数据和叶片总数，对各类分支的叶片归属度进行分析从而计算各类分支的叶团簇体积、各类分支的叶子点云密度和各类分支的叶密度，测量树高和各类分支的胸径；
步骤10：模型构建与风力加载下的参数分析：根据林分参数建立树体三维模型，所述林分参数包括树高、主枝干与次级枝干的夹角、各类分支的点云数据、各类分支的叶片总数、各类分支的胸径、各类分支的叶团簇体积、各类分支的叶密度和各类分支的叶团簇孔隙度，由单株活立木的树体三维模型组合成多株活立木的林分模型，将林分模型加载风力，根据湍流模型和流固耦合模型分析冠内动力学参量变化以及流固耦合过程，从而得到林分内部动态压力、风速度以及湍流动能强度分布。
作为本发明进一步改进的技术方案，所述的步骤1包括：
采用激光雷达扫描技术获取林木点云数据，对林木点云数据进行计算，实现枝叶分离；
所述的步骤2包括：
设枝干点云中的任意一点为
Figure GDA0003924421620000021
求出其近邻点
Figure GDA0003924421620000022
记为：
Figure GDA0003924421620000023
其中i代表第i个枝干点云，j代表第i个枝干点云的第j个近邻点；
预先设定点云数目阈值num1和半径阈值dist1，计算每个近邻点
Figure GDA0003924421620000031
与点
Figure GDA0003924421620000032
的距离，若距离小于半径阈值dist1的近邻点数目小于等于num1，则该点
Figure GDA0003924421620000033
为噪声点，并将其去除。
作为本发明进一步改进的技术方案，所述的步骤3包括：
(1)对去噪操作后的原始枝干点云Pb中的任意一点
Figure GDA0003924421620000034
用最小平方面来拟合其近邻点
Figure GDA0003924421620000035
将点
Figure GDA0003924421620000036
和近邻点
Figure GDA0003924421620000037
在平面做投影从而得到点
Figure GDA0003924421620000038
的投影点
Figure GDA0003924421620000039
和其近邻点
Figure GDA00039244216200000310
的投影点
Figure GDA00039244216200000311
将这些投影点在投影平面上做Delaunay三角剖分；
(2)计算n×n点云拉普拉斯加权矩阵L，加权矩阵L的第i行y列的元素值Li，y为：
Figure GDA00039244216200000312
其中E是所有近邻点构成的边的集合，αi，y、βi，y是相应于点
Figure GDA00039244216200000313
和近邻点
Figure GDA00039244216200000314
在其Delaunay三角剖分平面投影点
Figure GDA00039244216200000315
和
Figure GDA00039244216200000316
构成的边所属两个三角形的对角；
(3)寻找一组新的顶点Pb′使得LPb′＝0，为了防止零解，将方程LPb′＝0转换为下式(2)：
Figure GDA00039244216200000317
其中Pb为原始枝干点云，WL与WH分别为控制收缩与控制保持原有位置的力度的n×n矩阵，且均为对角矩阵，WL的第i个对角线元素记为WL，i，WH的第i个对角线元素记为WH，i，因此将公式(2)等价于极小化下面函数：
Figure GDA00039244216200000318
其中||WLLPb′||2为点云收缩力项，
Figure GDA00039244216200000319
为形状保持约束项，通过公式(3)的反复迭代，每步调整力度控制矩阵WL与WH，进而将枝干点云进行收缩。
作为本发明进一步改进的技术方案，所述的步骤4包括：
(1)根据单株活立木的垂直高度将收缩后的枝干点云数据自下而上切分为不同的层layer，分层的数目LevelNum由树高Ht和高度间隔Δh决定，即为：
Figure GDA00039244216200000320
(2)根据欧式距离进行聚类，设定距离阈值dist2，每一个高度分层中的枝干点云根据距离阈值dist2进行聚类归属，即将每一个高度分层的枝干点云之间距离少于dist2的点云集归属为一类；
(3)求取每类的中心点
Figure GDA0003924421620000041
其中j代表第layer层的第j个中心点，即为每层局部枝干的聚类中心点，进而获取整株活立木的枝干骨架的聚类中心点。
作为本发明进一步改进的技术方案，所述的步骤5包括：
(1)遍历每一高度分层的聚类中心点
Figure GDA0003924421620000042
按照空间近邻的方法，找寻下一高度layer-1层的结点中与当前layer层的结点距离最近的结点，找寻的节点即为当前结点的父结点
Figure GDA0003924421620000043
(2)所有聚类中心点
Figure GDA0003924421620000044
连接到下一高度layer-1层中只有唯一父结点，根节点
Figure GDA0003924421620000045
没有父结点，根据上述聚类中心点与其父节点的连接性，构建整个枝干的Connection连通性矩阵，整个枝干骨架分成若干段，每段用一个带有空间方向性的圆柱体
Figure GDA0003924421620000046
进行拟合，r代表圆柱体的半径，从而实现单株活立木完整骨架的重建。
作为本发明进一步改进的技术方案，所述的圆柱体的半径的设置方法为：
从所有聚类中心点中确定边缘节点
Figure GDA0003924421620000047
即没有任何
Figure GDA0003924421620000048
与
Figure GDA0003924421620000049
相连接的聚类中心点，利用Dijkstra最短路径算法确定
Figure GDA00039244216200000410
到最远的边缘节点
Figure GDA00039244216200000411
的距离
Figure GDA00039244216200000412
进而决定对应的枝干粗度，即树干模型的圆柱体Cl的半径，树干模型的每个圆柱体半径估算方程为：
Figure GDA00039244216200000413
其中
Figure GDA00039244216200000414
为树干模型的每个圆柱体半径，λ为系数；
同时，具有多个接通的儿子结点
Figure GDA00039244216200000415
的聚类中心点
Figure GDA00039244216200000416
记为分支交叉点
Figure GDA00039244216200000417
即分支结点
Figure GDA00039244216200000418
作为本发明进一步改进的技术方案，所述步骤6包括：
通过遍历聚类中心点中的根节点、父节点、分支结点、边缘节点和儿子节点分别计算分支结点下一高度layer-1层中的聚类中心点
Figure GDA0003924421620000051
指向分支结点
Figure GDA0003924421620000052
的向量与分支结点
Figure GDA0003924421620000053
指向上一高度layer+1层中的聚类中心点
Figure GDA0003924421620000054
的向量之间的角度，角度小所对应的枝干为主枝干，角度大的枝干为次级枝干，即：
Figure GDA0003924421620000055
标记主枝干和次级枝干。
作为本发明进一步改进的技术方案，所述步骤7包括：
(1)计算叶子点云数据中每一点
Figure GDA0003924421620000056
的灰度值；
(2)在空间位置上相近并且灰度值相近的像素点互相连接起来构成封闭的轮廓；
(3)根据步骤6已分类好的枝干骨架进行叶子点云的归属。
作为本发明进一步改进的技术方案，所述步骤8中的计算活立木的孔隙率包括：
由布尔原理任一高度z处沿天顶角θ的孔隙率为：
P(z,θ)＝exp(-NSz,θ(1-az,θ)) (7)；
其中Sz,θ为树在高度z以上部分沿方向θ在水平面内的投影面积，az,θ为单棵活立木在高度z以上部分沿θ方向的孔隙率，N为单位面积内树木棵树的期望值，其中az,θ为：
Figure GDA0003924421620000057
其中Gl(θ)为对于入射角为θ的光束，第l层切片平面的平均消光系数，M为切片的总层数，LD为冠层内叶面积密度，V(z)为单个树冠在高度z以上的体积。
作为本发明进一步改进的技术方案，所述步骤10中的流固耦合模型的动量方程为：
Figure GDA0003924421620000058
公式(9)中ρ为流体密度，p为压强；μt为湍流粘性系数；
Figure GDA0003924421620000059
为雷诺应力项；xi和xj均代表直角坐标系中三个坐标分量且xi和xj为不同方向上的坐标分量，i＝1,2,3，j＝1,2,3，i≠j；ui和uj均代表直角坐标系中瞬时速度沿x，y，z方向的三个分量且ui和uj为不同方向上的速度分量，i＝1,2,3，j＝1,2,3，i≠j；Si表示动量源项沿x，y，z方向的分量；
其中Si为：
Figure GDA0003924421620000061
其中Cd为树冠阻力系数；其中Cd为：
Figure GDA0003924421620000062
其中A为叶面积密度；
所述步骤10中的湍流模型采用k-ε模型。
本发明的有益效果为：本发明采用在中观尺度下的树木仿真建模方法，运用叶团簇理论的相关算法，把主枝干和每个次级枝干上所归属的叶片整体看成一个叶团簇，通过计算该叶团簇的体积并结合实际手动测量的叶子数据，进而推算其孔隙率，从而构建基于真实叶团簇参数的树体空间模型。同时，运用计算机仿真技术加载流场，刻画风力流场在树内的分布特征，进而找寻树冠内的静压差、湍动能与流速变化，探寻面向台力侵害下不同品种活立木的抗风性能与安全稳定性。本发明一方面计算复杂度适中，另一方面又能更好地描述不同树种的空间结构特征与生长参数变化，通过计算机仿真技术，实现活立木在台风下抗风性能的定性定量化评估，准确率较高，为树木的栽培种植与防风营造提供准确的理论依据。
附图说明
图1为本实施例的工作流程图。
图2为本实施例的两种橡胶树点云数据。
图3为本实施例的骨架收缩算法所用到的权值角度。
图4为本实施例用聚类算法树干高度分层及聚类中心点结果图。
图5为本实施例枝干重建两种橡胶树的结果图。
图6为本实施例提取两种树干的边缘结点与分支交叉点的结果图。
图7为本实施例中采用沿着根节点的路径比较枝干间夹角的方法来判断主枝干与次级枝干的示意图。
图8为本实施例主枝干与次级枝干分类后结果图。
图9为本实施例叶子按最短距离归属到各枝干上的结果图。
图10为不同入射角下的将树逆时针旋转示意图。
图11为本实施例基于叶团簇的两个橡胶树品种的三维重建与空间网格剖分示意图。
图12为应用数值仿真计算显示橡胶树(PR107)组成的林分与橡胶树(CATAS7-20-59)组成的林分流场计算结果的横截面图。
具体实施方式
下面根据图1至图12对本发明的具体实施方式作出进一步说明：
本实施例提供一种基于激光点云与空气动力学的活立木抗风性能分析方法，具体方法阐述如下。
(一)不同橡胶树种的骨架提取。
本实施例将基于地面激光雷达数据，结合机器视觉与图像图形学算法，以两株橡胶树为研究对象，结合圆柱拟合和Dijkstra最短路径算法，实现了活立木的精准重建与叶面归属度计算，并对其归属度进行分析。基于激光点云与空气动力学的活立木抗风性能分析方法的流程如图1所示。
(1.1)枝叶分离：
激光扫描后，获取到的仅是稠密的林木点云数据，并非是橡胶林木的具体器官信息，如单株枝干的空间信息、树冠的叶子点云等。因此，本实施例对林木点云进行计算，实现枝叶分离。结果图如图2所示。图2为枝叶分离扫描获得的两种橡胶树点云数据，图2的(a)图为橡胶树(PR107)枝叶分离结果，图2的(b)图为橡胶树(CATAS7-50-29)枝叶分离结果。
(1.2)激光扫描数据噪声点去噪操作：
由于实验所获取的原始点云数据中含有大量对实验结果造成干扰的点(称之为噪声点)，因此，在提取树干骨架前，要对扫描得到的单株活立木的n个树干点云Pb进行去除噪声点操作。
首先对于点云数据中的任意一点
Figure GDA0003924421620000071
求出其近邻点
Figure GDA0003924421620000072
记为：
Figure GDA0003924421620000073
其中i代表第i个枝干点云，j代表第i个枝干点云的第j个近邻点。
其次计算每个近邻点
Figure GDA0003924421620000074
与点
Figure GDA0003924421620000075
的距离，如若邻域点距离小于半径阈值dist1的点云数目≤num1(num1可根据扫描精度与扫描范围设定，为点云数目阈值)，则判定此点
Figure GDA0003924421620000076
为噪声点，并将其去除。
(1.3)树木枝干点云初步骨架收缩：
对于扫描的橡胶树木骨架点云而言，在每个骨架的连接处可能存在一对多的关系(多个分支)，因此，本节运用拉普拉斯算法对枝干点云进行收缩。主要思路是：首先对原始点云构建邻域关系，并计算点云的拉普拉斯加权矩阵。在此基础上，通过迭代求解拉普拉斯方程，并不断更新矩阵的收缩权和顶点约束权，实现对点云的收缩。
具体步骤为：
(1)首先，用最小平方面来拟合近邻点
Figure GDA0003924421620000081
将点
Figure GDA0003924421620000082
和近邻点
Figure GDA0003924421620000083
在平面做投影分别记为中心点
Figure GDA0003924421620000084
和其近邻点的投影点
Figure GDA0003924421620000085
将这些投影点在投影平面上做Delaunay三角剖分。
(2)计算n×n点云拉普拉斯加权矩阵，加权矩阵L的第i行y列的元素值Li,y如下式：
Figure GDA0003924421620000086
其中E是所有近邻点构成的边的集合，αi,y、βi,y是相应于
Figure GDA0003924421620000087
和
Figure GDA0003924421620000088
在其Delaunay三角剖分平面投影点
Figure GDA0003924421620000089
和
Figure GDA00039244216200000810
构成的边所属两个三角形的对角,如图3所示。
(3)接着，需找出一组新的顶点Pb′使得LPb′＝0，然而L是奇异的，为了防止零解，将方程LPb′＝0转换为下式(2)：
Figure GDA00039244216200000811
其中Pb为原始枝干点云，WL与WH分别为控制收缩与控制保持原有位置的力度的n×n矩阵，且均为对角矩阵,WL(WH)的第i个对角线元素记为WL,i(WH,i)。因此，方程(2)式即可等价于极小化下面函数：
Figure GDA00039244216200000812
方程(3)式第一项为点云收缩力项，第二项为形状保持约束项。通过方程(3)式的反复迭代，每步调整力度控制矩阵WL与WH，进而将枝干进行合理的收缩。
(1.4)枝干高度分层及聚类中心点求取：
根据单株活立木的垂直高度将收缩后的枝干点云数据自下而上切分为不同的层；根据欧式距离进行聚类，设定距离阈值，每一个高度分层中的枝干点云根据距离阈值进行聚类归属，即将每一个高度分层的枝干点云之间距离少于距离阈值的点云集归属为一类；求取每类的中心点，即每一个高度分层的局部枝干的聚类中心点。
具体步骤为：
根据单株木的垂直高度，将点云数据自下而上切分为不同的层(layer)。其中分层的数目LevelNum是由树高Ht和高度间隔Δh决定，即为：
Figure GDA0003924421620000091
由于树木枝干存在分支情况，会造成同一层上存在多个分支，因此需要对存在分支的层级进行处理。本算法根据欧式距离进行聚类，设定一定的距离阈值dist2，每一个高度分层中的枝干点云根据dist2进行聚类归属，即每一个高度层的点云之间距离少于dist2的点云集即归属为一类。
然后求取每类的中心点
Figure GDA0003924421620000092
其中j代表第layer层的第j个中心点，即为每层局部枝干的聚类中心点。根据这些中心点，以此标记不同层次的枝干分布区域，即获取整株橡胶树干骨架的聚类中心点。橡胶树种林段(PR107)中树以第19、20和21层，橡胶树种林段(CATAS7-20-59)中树以第23、24和25层为例，分层结果和每类的聚类中心点结果如图4所示。图4为用聚类算法树干高度分层及聚类中心点结果图，图4的(a)图为橡胶树种林段(PR107)中单株树计算得到第19、20和21层的聚类中心点，分别用不同深度颜色标记，图4的(b)图为橡胶树种林段(CATAS7-20-59)中单株树计算得到23、24和25层的聚类中心点，分别用不同深度颜色标记。
(1.5)基于中心点的圆柱体拟合与骨架重建：
遍历每一个高度分层的聚类中心点，根据聚类中心点并采用带有空间方向性的圆柱体进行拟合每一个高度分层的枝干骨架，从而实现单株活立木的完整骨架的重建。
具体步骤为：
遍历每一层的聚类中心点
Figure GDA0003924421620000093
按照空间近邻的方法，找寻当前layer层的结点与下一高度layer-1层的结点中距离最近的结点，即当前结点的父结点
Figure GDA0003924421620000094
枝干重建依次根据
Figure GDA0003924421620000095
来重建树木骨架，具体如下：由于植物生长规律，所有聚类中心点
Figure GDA0003924421620000096
连接到
Figure GDA0003924421620000097
中只有唯一一个父结点(根节点
Figure GDA0003924421620000098
没有父结点)，因此，根据上述结点的连接性，构建整个树干的Connection连通性矩阵，整个树木枝干骨架分成若干段，每段用一个带有空间方向性的圆柱体
Figure GDA0003924421620000099
进行拟合，r代表圆柱体的半径，从而实现单株木完整骨架的重建。具体如图5所示。图5为本文方法枝干重建两种橡胶树的结果，图5的(a)图为橡胶树种林段(PR107)中单株树每层的聚类中心点；图5的(b)图为根据5的(a)图中的聚类中心生成不同半径的圆柱体拟合枝干，并顺枝干生长方向连接重建树干；图5的(c)图为橡胶树种林段(CATAS7-20-59)中单株树每层的聚类中心点；图5的(d)图为根据5的(c)图中的聚类中心生成不同半径的圆柱体拟合枝干，并顺枝干生长方向连接重建树干。
(1.6)圆柱体半径设置：
从树木生理的角度上讲，树木的上一级枝干维持着下一级枝干的生长，越靠近根部的枝干承担着更多的营养供给和重力支持，所相应的枝干半径也会越大。因此，本方法首先从所有聚类中心点中确定边缘点
Figure GDA0003924421620000101
即没有任何
Figure GDA0003924421620000102
与
Figure GDA0003924421620000103
相连接的聚类中心点，接着利用Dijkstra最短路径算法确定
Figure GDA0003924421620000104
到最远的边缘点
Figure GDA0003924421620000105
的距离
Figure GDA0003924421620000106
进而根据整体长度决定对应的枝干粗度，即树干模型的圆柱体Cl的半径r，树干模型的每个圆柱体半径估算方程为：
Figure GDA0003924421620000107
式中，λ为系数(与点云密度有关)。
同时，对于聚类中心点
Figure GDA0003924421620000108
其具有多个接通的儿子结点
Figure GDA0003924421620000109
的，则被认为是分支交叉点
Figure GDA00039244216200001010
将计算一并求取得到。具体找寻的枝干骨架边缘点与分支交叉点(分支结点)如图6所示。图6为本文方法提取两种树干的边缘结点与分支交叉点的结果图；图6的(a)图和6的(b)图分别为橡胶树(PR107)与橡胶树(CATAS7-20-59)边缘结点与分支交叉点(分支结点)提取的结果。
(二)树干分类及叶面归属。
(2.1)橡胶枝干骨架分类：
分叉处的骨架点通常是不同分枝的分离点。因而需要重点关注骨架点当中的特殊点，根据树型数据结构的基本特征，总结有如下三类特殊点：
1)根节点
Figure GDA00039244216200001011
没有父节点
Figure GDA00039244216200001012
2)分支结点
Figure GDA00039244216200001013
存在两个或以上相连通的儿子节点
Figure GDA00039244216200001014
3)边缘节点
Figure GDA00039244216200001015
没有相连通的儿子节点
Figure GDA00039244216200001016
依据以上三种特殊点性质能够很方便的通过遍历骨架点来完成对不同树枝骨架点的分离，从而进行整株树不同枝干的分类。橡胶树的结构一般是树干下部分分成多个次级分支，因此按Dijkstra搜索，从根结点
Figure GDA0003924421620000111
开始，将与根结点
Figure GDA0003924421620000112
相连通的结点入栈，直至到分支结点
Figure GDA0003924421620000113
结束，其次分别计算分支结点下一高度layer-1层中聚类中心点
Figure GDA0003924421620000114
指向分支结点
Figure GDA0003924421620000115
的向量，与分支结点
Figure GDA0003924421620000116
指向分支结点上一高度layer+1层中的聚类中心点
Figure GDA0003924421620000117
的向量之间的角度θ(如图7所示)，然后比较θ1和θ2大小，角度小的枝干被认为是主枝干，大的被认为是次级枝干，即：
Figure GDA0003924421620000118
图7为本文中采用沿着根节点的路径比较枝干间夹角的方法来判断主枝干与次级枝干。不同层聚类中心与分支点所成角度最小变化的路径为主枝干所在的路径。
最后将栈中的点云数据归为一类，用相同的颜色表示(如图8所示)。本实施例研究主枝与侧枝夹角对于树木抗风性的影响，因此在对拟合后枝干分类时只对主枝干和次级枝干分类标记。图8为主枝干与次级枝干分类后结果，图8的(a)图为橡胶树(PR107)主枝干与各次级枝干分别用不同颜色标记，图8的(b)图为橡胶树(CATAS7-20-59)主枝干与各次级枝干分别用不同颜色标记。
(2.2)叶子点云枝干归属度计算：
基于分类好的枝干骨架，结合空间分水岭算法，完成对叶子点云数据的归属。空间分水岭算法是一种用于图像分割的经典算法，是基于拓扑理论的数学形态学的分割方法。其基本思想是将图像看成是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。因此，根据算法，首先，计算叶子点云数据中每一点
Figure GDA0003924421620000119
的灰度值；其次将在空间位置上相近并且灰度值相近的像素点互相连接起来构成封闭的轮廓；最后根据已分类好的枝干骨架，进行叶子点云的归属。具体叶子点云数据的归属结果如图9所示。图9为叶子按最短距离归属到各枝干上，并与该枝干同色标记，图9的(a)图和9的(b)图分别为橡胶树(PR107)与橡胶树(CATAS7-20-59)叶子点云分类结果。
(2.3)不同橡胶树的孔隙率求取：
假设树冠在水平面内随机分布且相互独立，并忽略树干、树枝影响，由布尔原理任一高度z处沿天顶角θ的孔隙率为：
P(z，θ)＝exp(-NSz，θ(1-az，θ)) (7)；
(7)式中，Sz，θ为树在高度z以上部分沿方向θ在水平面内的投影面积。az，θ为单棵树在高度z以上部分沿θ方向的孔隙率。Sz，θ(1-az，θ)的几何意义为单个带孔树冠在在高度z以上部分沿θ方向在水平面内的投影面积。N为单位面积内树木棵树的期望值，与树木的冠幅有关。
在树冠内部为浑浊介质的假设下，az，θ的表示为下式：
Figure GDA0003924421620000121
(8)式中，(x，y，z)为树冠内的任意一点，s(x，y，z，θ)为光子沿θ方向到达该点经过的光学路径，G(θ)为投影函数，LD为冠层内叶面积密度。s是光子到达每一点(x，y，z)前的光学路径，其值很难确定，造成很难计算az，θ。本实施例计算其近似解，即：
Figure GDA0003924421620000122
(9)式中，
Figure GDA0003924421620000123
为光子沿θ方向到达z平面经过的近似光学路径，且：
Figure GDA0003924421620000124
(10)式中，V(z)为单个树冠在高度z以上的体积。
消光系数G(θ)是由入射角度方向和叶片分布方向所决定的，定义为垂直于光束方向的平面上的单位叶面积的平均投影。为了计算消光系数G(θ)：
首先，为了求解在天顶角θ方向上的G(θ)，把图10的(a)图中的三维树体逆时针旋转θ(如图10的(b)图所示)，再进行体素化(如图10的(c)图所示)。通过定义每个体素的宽度w，长度l和高度h，点云数据P将被分成m×n×p个体素，其中m＝(xmax-xmin)/w，n＝(ymax-ymin)/l，p＝(zmax-zmin)/h；w、l、h均为实数。将点云体素化处理后，定义每一个行或列层为提取了部分点云数据的平面，即切片平面。步骤如图10所示，以树(PR107)为例。图10为了计算不同入射角下的树冠消光系数，将树逆时针旋转，便于参数反演计算和体素化。
其次，计算平均消光系数：对于方向为θ的入射光束，第l层切片平面的平均投影系数可以通过将方位角
Figure GDA0003924421620000131
在[0,2π]范围内积分为式(11),其中
Figure GDA0003924421620000132
为天顶角θ与方位角
Figure GDA0003924421620000133
所构成的光线入射角：
Figure GDA0003924421620000134
再次，求取给定方向的消光系数：首先考虑方向为
Figure GDA0003924421620000135
的光束，取垂直于该光束的切片平面，设该切片平面里每个非空体素内的点数是n。分为以下几种情况：
(a)如果n＝1，激光扫描的采样分辨率为s，消光系数为
Figure GDA0003924421620000136
(b)如果n＝2，构造一条线来表示树叶面积，消光系数可以通过计算这条线的投影长度和原线段的长度之比来获得。
(c)如果n≥3，首先计算所有点是共面还是在同一线上。如果它们在同一条线上，它将与n＝2的情况归为一类。如果在同一平面上，将构造一个三角形，三角形的面积为A，其次根据三个点的投影坐标计算投影面积AP。消光系数为AP/A。
(d)如果体素内的所有原始点既不是共面也不是同一直线，则将构造三维凸包，并且该三维凸包的总表面积一半将用于表示叶面积A。然后，将所有原始点投影到与线性样方方向一致的平面上，由投影点构成的二维凸包的面积将代表投影的叶面积AP，消光系数为AP/A。
基于前面讨论的过程，计算了所有非空体素的消光系数。通过总结方向为
Figure GDA0003924421620000137
单个切片平面的消光系数，得到：
Figure GDA0003924421620000138
其中(12)式中，ξ代表第l层切片平面的第ξ个非空体素，S是第l层切片平面中非空体素的总数。将式(12)带入式(11)中，即可求得对于入射角为θ的光束，第l层切片平面的平均消光系数Gl(θ)。对于整个树冠而言，消光系数可由下式(13)得到：
Figure GDA0003924421620000139
(13)式中，M是指切片的总层数。将(13)式带入到(9)式中：
Figure GDA00039244216200001310
其中，冠层内叶面积密度LD由表4可知，单个树冠在高度z以上的体积V(z)由表3可知，Sz,θ是由图10的(c)图沿入射角为θ的投影计算得到。即可以求得az,θ，将az,θ的值带入式(7)：
Figure GDA00039244216200001311
N与树木的冠幅有关，可由表3得到。便可以分别求得两棵橡胶树的孔隙率。
(三)流体模型计算。
自然界中流体流动遵循质量守恒定律、动量守恒定律、能量守恒定律。流体运动形式分为两大类，层流和湍流。因此在构建计算模型时必须满足上述条件。
(3.1)模型构建：
质量守恒定律是单位时间内流体微元体中质量的增加，等于同一时间间隔内流入该微元体的净质量。质量守恒方程又称为连续性方程。
本实施例研究问题中空气流场流速均小于1000m/s，视为不可压缩流体，且大气环境为中性，忽略coriolis力(其中Coriolis力是对旋转体系中进行直线运动的质点由于惯性相对于旋转体系产生的直线运动的偏移，来自于流体运动的惯力)，因此选用稳态不可压缩流动模型模拟流畅内部流动。假设流体密度ρ不随时间变化，根据质量守恒定律，在流场中取任一流体微元则其连续性方程为：
Figure GDA0003924421620000141
式中xi代表直角坐标系中三个坐标分量，ui(i＝1,2,3)代表直角坐标系中瞬时速度沿x，y，z方向的三个分量。
动量守恒定律表述为微元体中流体的动量对时间的变化率等于外界作用在该微元体上的各种力之和。其动量方程为：
Figure GDA0003924421620000142
式中ρ为流体密度，t为时间。p为压强，v为空气动力学强度，
Figure GDA0003924421620000143
为雷诺应力项。xi和xj均代表直角坐标系中三个坐标分量且xi和xj为不同方向上的坐标分量，i＝1,2,3，j＝1,2,3，i≠j；ui和uj均代表直角坐标系中瞬时速度沿x，y，z方向的三个分量且ui和μj为不同方向上的速度分量，i＝1,2,3，j＝1,2,3，i≠j；
雷诺应力项根据网格变数解析定义为：
Figure GDA0003924421620000144
式中μt为湍流粘性系数，δij为Kroneckerδ函数
Figure GDA0003924421620000145
k为湍流动能。
Figure GDA0003924421620000146
式中ε为湍流耗散率，Cμ为经验常数，取值详见表1。
Figure GDA0003924421620000151
Figure GDA0003924421620000152
流体物体在不可压缩流体运动的形式分为层流与湍流。区分层流与湍流的参数为雷诺数Re。雷诺数是流体力学中表征粘性影响的相似准则数，用以判别粘性流体流动状态的一个无因次数群。
Re＝ρul/μt (21)；
式中ρ为流体密度，u为流体速度，l为流场长度，μt为湍流粘性系数。雷诺数较小时(<300)，粘滞力对流场的影响大于惯性，流场中流速的扰动会因粘滞力而衰减，流体流动稳定，为层流；反之，若雷诺数较大时，惯性对流场的影响大于粘滞力，流体流动较不稳定，流速的微小变化容易发展、增强，形成紊乱、不规则的湍流流场。因为强风压下风速过大，因此忽略层流，引入湍流模型。
(3.2)湍流方程与湍流耗散方程：
湍流模型有多种模拟方法，本算法基于对实验要求，选择了计算量与精度适中的标准k-ε模型。其输运方程为：
Figure GDA0003924421620000153
Figure GDA0003924421620000154
上式中Cε1，Cε2，σk，σε为模型经验取值，详见表1。表1为湍流模型中相关参数值：
C<sub>μ</sub> C<sub>ε1</sub> C<sub>ε2</sub> σ<sub>k</sub> σ<sub>ε</sub>
0.09 1.42 1.92 1.0 1.3
(3.3)流固耦合与多孔介质：
建立模型时，叶团簇采用多孔介质模拟，为了模拟强风压下林分内部流固耦合过程需在动量方程中添加动量源项Si，分析树冠内部的空隙部分阻力对流体运动的影响。则流固耦合模型的动量方程为：
Figure GDA0003924421620000155
式中p为压强；μt为湍流粘性系数；
Figure GDA0003924421620000156
为雷诺应力项；xi和xj均代表直角坐标系中三个坐标分量且xi和xj为不同方向上的坐标分量，i＝1,2,3，j＝1,2,3，i≠j，i和j均代表方向；ui和uj均代表直角坐标系中瞬时速度沿x，y，z方向的三个分量且ui和uj为不同方向上的速度分量，i＝1,2,3，j＝1,2,3，i≠j；Si表示动量源项沿x，y，z方向的分量。
Figure GDA0003924421620000161
uj为与Si不同方向上的速度分量,Cd为树冠阻力系数。
Figure GDA0003924421620000162
式中A为叶面积密度，其值由表4可得。
本实施例可将式(18)、式(19)、式(20)带入式(22)、式(23)中可得到空间中任一网格的速度u，将风速u带入式(18)、式(19)、式(20)中，可得到当前点湍流动能k、耗散率ε和湍流粘度μt。联立式(22)、式(23)、式(24)可根据风速计算出动态压力p。
(四)实验结果分析。
(4.1)树木重建精度验证：
上文主枝干重建与叶片归属是否符合客观事实，很大程度上取决于参数的选取的合理性。因此，要求所选用的参数必须能够正确地反映真实树木的情况，从而满足主枝干重建和叶片归属的精准要求。根据我们多次人工测算，本文算法涉及到的扫描参数与算法参数详细介绍如表2所示：表2为本文方法扫描与算法所用的相关参数：
Figure GDA0003924421620000163
另外，本文对两棵橡胶树的树冠体积、树高和冠幅等进行了参数分析，各参数的人工实测数据罗列在表3中。表3为人工测量橡胶树生长属性值：
Figure GDA0003924421620000164
Figure GDA0003924421620000171
为更好地研究树木叶片的生长情况，本实施例对橡胶树进行叶片枝干归属度的分析。将橡胶树的枝干分为五类(主枝干和次级枝干的数量决定了枝干的类别数，本实施例中的橡胶树的次级枝干共四个，主枝干1个，因此分为五类，即五类分支)，其分支A、B、C、D、E结果如图9的(a)图所示。另一棵橡胶树也分为了五类，如图9的(b)图所示。统计两棵橡胶树和各类分支的点云数据、叶片总数，对各类分支的叶片枝干归属度进行分析，以及对各类分支的胸径、叶团簇体积和叶密度等进行分析。两棵橡胶树的各参数的定量估计如表4所示，与人工实测数据相比较，其相关度也罗列在表4中。
以保证人工测量数据的准确度，我们通过反复测量与校正。由表4计算表明：主枝干上的叶子生长相对次级枝干较茂密；次级枝干上的叶团簇体积与胸径成正比。
表4为各枝干点云分类结果与人工测量值比较：
Figure GDA0003924421620000172
对各类分支的叶片归属度进行分析从而计算各类分支的叶团簇体积，叶子点云密度为叶子点云数除以叶团簇体积，叶密度为叶子面积除以叶团簇体积。树高和各分支的胸径通过人工测量获取。
(4.2)基于真实叶团簇参数的橡胶树空间模型构建：
根据前文求取的林分参数(具体包括树高、主枝干与次级枝干的夹角、各类分支的点云数据、叶片总数、胸径、叶团簇体积、叶密度和叶团簇孔隙度(即叶团簇的孔隙率)等)建立真实树体三维模型，由单棵树(PR107)的真实树体三维模型组合成多株(PR107)的林分模型；由单棵树(CATAS7-20-59)的真实树体三维模型组合成多株(CATAS7-20-59)的林分模型，将两个林分模型分别加载相同风力，并结合湍流模型和流固耦合模型，分析冠内动力学参量变化，以及流固耦合过程。树种模型参数如表5所示，枝干C为主枝干。
表5为两种橡胶树(PR107和CATAS7-20-59)枝干重建模型具体信息：
Figure GDA0003924421620000181
具体单株模型如图11的(a)图和11的(b)图所示，通过实际样例与真实种植间距，建立林分模型如图11的(c)图和11的(d)图所示。图11为基于叶团簇的两个橡胶树品种的三维重建与空间网格剖分，以便于流体计算仿真；图11的(a)图和11的(b)图分别表示两个橡胶树品种基于叶团簇参数重建的橡胶树空间模型。图11的(c)图和11的(d)图表示由重建模型构建的不同林段(两个橡胶树品种)的对应模型，并进行三维网格剖分。
(4.3)树木在风力加载下的参数分析：
将两个林分模型置于入强风力场，根据求解模型可得到林分内部动态压力，速度，湍流动能强度分布如图12所示，具体地，将式(18)、式(19)、式(20)带入式(22)、式(23)中可得到空间中任一网格的速度u，将风速u带入式(18)、式(19)、式(20)中，可得到当前点湍流动能k、耗散率ε和湍流粘度μt。联立式(22)、式(23)、式(24)可根据风速计算出动态压力p。图12为应用数值仿真计算显示橡胶树(PR107)组成的林分与橡胶树(CATAS7-20-59)流场计算结果的横截面。林分(PR107)中动态压力与湍流强度更加强烈会增加损伤树体的可能性。
由图12可知，将林段(PR107)与林段(CATAS7-20-59)数值计算结果置于同一图例下，可见其分布的差异性。由动态压力分布图可见林段(PR107)内的动压值明显高于林段(CATAS7-20-59)的动压值，林段(PR107)内部压差较大，且压力分布不均易造成林分迎风面树木倒折；由速度图可知在林段(PR107)内的速度波动较大而林段(CATAS7-20-59)速度均一性明显，针对林内风速而言，林段(PR107)内部由于树冠的孔隙度较小并且其主枝角小于林段(CATAS7-20-59)，易造成风阻，强风力与林冠易产生相互作用力，致使林冠内枝干易折断且受损严重；从湍动能强度分析可知，林段(PR107)湍动能强度明显强于林段(CATAS7-20-59)，林段(PR107)林分外部被不稳定的湍流包裹着，在林分背风面形成最大湍流团，此流极不稳定，易对整个林分造成损伤，林段(CATAS7-20-59)中背风面未形成较大湍流，证明其林分内部承受的强风力被减弱，对树体造成的伤害远小于林段(PR107)。在林分中，局部压力，速度及湍流的差异性会导致冠内枝干的折断及树体损伤，而林段(CATAS7-20-59)内部压力，速度与湍流分布较为均匀，因此林段(PR107)在强风压下更容易断折和损害。
本实施例针对地面激光雷达数据和移动激光雷达数据，分别对其进行了林分参数的提取和反演。对于地面激光雷达数据，本次采用单株完整的橡胶树活立木为研究对象，结合机器视觉与图像图形学算法，重建了单株树木骨架的三维模型，及叶子归属计算，精细而准确的再现了植株在空间中的生理特征。对于移动激光雷达数据，可以采用自下而上的移动扫描方式获取两个橡胶林段点云数据，并设计图形图像学算法实现对橡胶林段的枝叶分离和单株提取。并且通过对林段(PR107)和林段(CATAS7-20-59)中橡胶树倾斜角度的对比分析发现，其林段(PR107)受风力侵害程度要大于林段(CATAS7-20-59)，也更易倒伏。
针对单株活立木三维骨架的重建，本实施例提出了一种快速、方便的算法，更加直观的展现了树木骨架的拓扑结构。通过对叶子点云数据归属度分析，研究结果与原始扫描数据在形状和分布上有较好的相似度，保持了原始数据的大部分点云特征，证明了该算法的可行性。并且根据叶面归属计算后的结果，可以快速有效的了解到该树木每一个树枝上的叶子附属情况，便于了解树木的长势和叶子的生长情况。此外，由于树木枝繁叶茂，树冠郁闭度大，通常的扫描方式，会因为树木存在遮挡的问题，导致树冠的小分枝常常在扫描中丢失，使得在获取的林木参数常常存在误差。而本实施例在获取橡胶林数据时采用的是自下而上的方式，能够完全避免此问题，获取到全面并无遮挡的点云数据。本实施例实现了树木三维骨架的重建、叶子点云数据的归属计算、枝叶分离和单株树木的株株分离，并未对树叶进行建模，通过计算林分内部动态压力，速度，湍流动能强度分布从而实现对活立木抗风性能进行有效分析。
本实施例针对橡胶树的两个品种(PR107与CATAS7-20-59)，采用在中观尺度下的树木仿真建模方法，运用叶团簇理论的相关算法，把每个次级分枝上所归属的叶片整体看成一个叶团簇，通过计算该叶团簇的体积并结合实际手动测量的叶子数据，进而推算其孔隙率，从而构建基于真实叶团簇参数的两个品种橡胶树的树体空间模型。同时，运用计算机仿真技术加载流场，刻画风力流场在树内的分布特征，进而找寻树冠内的静压差、湍动能与流速变化，探寻面向台力侵害下不同品种的橡胶树的抗风性能与安全稳定性。本文的方法一方面计算复杂度适中，另一方面又能更好地描述不同树种的空间结构特征与生长参数变化，通过计算机仿真技术，实现橡胶树在台风下抗风性能的定性定量化评估，为橡胶树的栽培种植与防风营造提供准确的理论依据。
本发明的保护范围包括但不限于以上实施方式，本发明的保护范围以权利要求书为准，任何对本技术做出的本领域的技术人员容易想到的替换、变形、改进均落入本发明的保护范围。

Claims (10)

1.一种基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，包括：
步骤1、枝叶分离：获取林木点云数据，对林木点云数据进行计算，实现枝叶分离；
步骤2、去噪操作：对获得的单株活立木的n个枝干点云进行去除噪声点操作；
步骤3、枝干点云骨架收缩：采用拉普拉斯算法对枝干点云进行收缩；
步骤4、枝干高度分层与聚类中心点求取：根据单株活立木的垂直高度将收缩后的枝干点云数据自下而上切分为不同的层；根据欧式距离进行聚类，设定距离阈值，每一个高度分层中的枝干点云根据距离阈值进行聚类归属，即将每一个高度分层的枝干点云之间距离少于距离阈值的点云集归属为一类；求取每类的中心点，即每一个高度分层的局部枝干的聚类中心点；
步骤5、圆柱体拟合与骨架重建：遍历每一个高度分层的聚类中心点，根据聚类中心点并采用带有空间方向性的圆柱体进行拟合每一个高度分层的枝干骨架，从而实现单株活立木的完整骨架的重建；
步骤6、枝干骨架分类：通过遍历聚类中心点中的根节点、父节点、分支结点、边缘节点和儿子节点将整株活立木不同枝干骨架分类为主枝干和次级枝干；
步骤7、叶子点云归属：基于分类好的枝干骨架并结合空间分水岭算法完成对活立木叶子点云数据的归属；
步骤8、孔隙率求取：计算活立木的孔隙率；
步骤9、叶片归属度分析：根据主枝干和次级枝干的数量将活立木的枝干分为多类，统计各类分支的点云数据和叶片总数，对各类分支的叶片归属度进行分析从而计算各类分支的叶团簇体积、各类分支的叶子点云密度和各类分支的叶密度，测量树高和各类分支的胸径；
步骤10：模型构建与风力加载下的参数分析：根据林分参数建立树体三维模型，所述林分参数包括树高、主枝干与次级枝干的夹角、各类分支的点云数据、各类分支的叶片总数、各类分支的胸径、各类分支的叶团簇体积、各类分支的叶密度和各类分支的叶团簇孔隙度，由单株活立木的树体三维模型组合成多株活立木的林分模型，将林分模型加载风力，根据湍流模型和流固耦合模型分析冠内动力学参量变化以及流固耦合过程，从而得到林分内部动态压力、风速度以及湍流动能强度分布。
2.根据权利要求1所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，
所述的步骤1包括：
采用激光雷达扫描技术获取林木点云数据，对林木点云数据进行计算，实现枝叶分离；
所述的步骤2包括：
设枝干点云中的任意一点为
Figure FDA0003924421610000021
求出其近邻点
Figure FDA0003924421610000022
记为：
Figure FDA0003924421610000023
其中i代表第i个枝干点云，j代表第i个枝干点云的第j个近邻点；
预先设定点云数目阈值num1和半径阈值dist1，计算每个近邻点
Figure FDA0003924421610000024
与点
Figure FDA0003924421610000025
的距离，若距离小于半径阈值dist1的近邻点数目小于等于num1，则该点
Figure FDA0003924421610000026
为噪声点，并将其去除。
3.根据权利要求1所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，所述的步骤3包括：
(1)对去噪操作后的原始枝干点云Pb中的任意一点
Figure FDA0003924421610000027
用最小平方面来拟合其近邻点
Figure FDA0003924421610000028
将点
Figure FDA0003924421610000029
和近邻点
Figure FDA00039244216100000210
在平面做投影从而得到点
Figure FDA00039244216100000211
的投影点
Figure FDA00039244216100000212
和其近邻点
Figure FDA00039244216100000213
的投影点
Figure FDA00039244216100000214
将这些投影点在投影平面上做Delaunay三角剖分；
(2)计算n×n点云拉普拉斯加权矩阵L，加权矩阵L的第i行y列的元素值Li,y为：
Figure FDA00039244216100000215
其中E是所有近邻点构成的边的集合，αi,y、βi,y是相应于点
Figure FDA00039244216100000216
和近邻点
Figure FDA00039244216100000217
在其Delaunay三角剖分平面投影点
Figure FDA00039244216100000218
和
Figure FDA00039244216100000219
构成的边所属两个三角形的对角；
(3)寻找一组新的顶点Pb′使得LPb′＝0，为了防止零解，将方程LPb′＝0转换为下式(2)：
Figure FDA00039244216100000220
其中Pb为原始枝干点云，WL与WH分别为控制收缩与控制保持原有位置的力度的n×n矩阵，且均为对角矩阵,WL的第i个对角线元素记为WL,i，WH的第i个对角线元素记为WH,i，因此将公式(2)等价于极小化下面函数：
Figure FDA00039244216100000221
其中||WLLPb′||2为点云收缩力项，
Figure FDA00039244216100000222
为形状保持约束项，通过公式(3)的反复迭代，每步调整力度控制矩阵WL与WH，进而将枝干点云进行收缩。
4.根据权利要求1所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，所述的步骤4包括：
(1)根据单株活立木的垂直高度将收缩后的枝干点云数据自下而上切分为不同的层layer，分层的数目LevelNum由树高Ht和高度间隔Δh决定，即为：
Figure FDA0003924421610000031
(2)根据欧式距离进行聚类，设定距离阈值dist2，每一个高度分层中的枝干点云根据距离阈值dist2进行聚类归属，即将每一个高度分层的枝干点云之间距离少于dist2的点云集归属为一类；
(3)求取每类的中心点
Figure FDA0003924421610000032
其中j代表第layer层的第j个中心点，即为每层局部枝干的聚类中心点，进而获取整株活立木的枝干骨架的聚类中心点。
5.根据权利要求4所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，所述的步骤5包括：
(1)遍历每一高度分层的聚类中心点
Figure FDA0003924421610000033
按照空间近邻的方法，找寻下一高度layer-1层的结点中与当前layer层的结点距离最近的结点，找寻的节点即为当前结点的父结点
Figure FDA0003924421610000034
(2)所有聚类中心点
Figure FDA0003924421610000035
连接到下一高度layer-1层中只有唯一父结点，根节点
Figure FDA0003924421610000036
没有父结点，根据上述聚类中心点与其父节点的连接性，构建整个枝干的Connection连通性矩阵，整个枝干骨架分成若干段，每段用一个带有空间方向性的圆柱体
Figure FDA0003924421610000037
layer＝2,3,4，……,LevelNum进行拟合，r代表圆柱体的半径，从而实现单株活立木完整骨架的重建。
6.根据权利要求5所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，所述的圆柱体的半径的设置方法为：
从所有聚类中心点中确定边缘节点
Figure FDA00039244216100000314
Figure FDA0003924421610000038
即没有任何
Figure FDA0003924421610000039
与
Figure FDA00039244216100000310
相连接的聚类中心点，利用Dijkstra最短路径算法确定
Figure FDA00039244216100000311
到最远的边缘节点
Figure FDA00039244216100000312
的距离
Figure FDA00039244216100000313
进而决定对应的枝干粗度，即树干模型的圆柱体Cl的半径，树干模型的每个圆柱体半径估算方程为：
Figure FDA0003924421610000041
其中
Figure FDA0003924421610000042
为树干模型的每个圆柱体半径，λ为系数；
同时，具有多个接通的儿子结点
Figure FDA0003924421610000043
的聚类中心点
Figure FDA0003924421610000044
记为分支交叉点
Figure FDA0003924421610000045
即分支结点
Figure FDA00039244216100000413
Figure FDA0003924421610000046
7.根据权利要求1所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，所述步骤6包括：
通过遍历聚类中心点中的根节点、父节点、分支结点、边缘节点和儿子节点分别计算分支结点下一高度layer-1层中的聚类中心点
Figure FDA0003924421610000047
指向分支结点
Figure FDA0003924421610000048
的向量与分支结点
Figure FDA0003924421610000049
指向上一高度layer+1层中的聚类中心点
Figure FDA00039244216100000410
的向量之间的角度，角度小所对应的枝干为主枝干，角度大的枝干为次级枝干，即：
Figure FDA00039244216100000411
标记主枝干和次级枝干。
8.根据权利要求1所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，所述步骤7包括：
(1)计算叶子点云数据中每一点
Figure FDA00039244216100000412
的灰度值；
(2)在空间位置上相近并且灰度值相近的像素点互相连接起来构成封闭的轮廓；
(3)根据步骤6已分类好的枝干骨架进行叶子点云的归属。
9.根据权利要求1所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，所述步骤8中的计算活立木的孔隙率包括：
由布尔原理任一高度z处沿天顶角θ的孔隙率为：
P(z,θ)＝exp(-NSz,θ(1-az,θ)) (7)；
其中Sz,θ为树在高度z以上部分沿方向θ在水平面内的投影面积，az,θ为单棵活立木在高度z以上部分沿θ方向的孔隙率，N为单位面积内树木棵树的期望值，其中az,θ为：
Figure FDA0003924421610000051
其中Gl(θ)为对于入射角为θ的光束，第l层切片平面的平均消光系数，M为切片的总层数，LD为冠层内叶面积密度，V(z)为单个树冠在高度z以上的体积。
10.根据权利要求1所述的基于激光点云与空气动力学的活立木抗风性能分析方法，其特征在于，
所述步骤10中的流固耦合模型的动量方程为：
Figure FDA0003924421610000052
公式(9)中ρ为流体密度，p为压强；μt为湍流粘性系数；
Figure FDA0003924421610000053
为雷诺应力项；xi和xj均代表直角坐标系中三个坐标分量且xi和xj为不同方向上的坐标分量，i＝1,2,3，j＝1,2,3，i≠j；ui和uj均代表直角坐标系中瞬时速度沿x，y，z方向的三个分量且ui和uj为不同方向上的速度分量，i＝1,2,3，j＝1,2,3，i≠j；Si表示动量源项沿x，y，z方向的分量；
其中Si为：
Figure FDA0003924421610000054
其中Cd为树冠阻力系数；其中Cd为：
Figure FDA0003924421610000055
其中A为叶面积密度；
所述步骤10中的湍流模型采用k-ε模型。
CN201811322277.0A 2018-11-08 2018-11-08 基于激光点云与空气动力学的活立木抗风性能分析方法 Active CN109446691B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811322277.0A CN109446691B (zh) 2018-11-08 2018-11-08 基于激光点云与空气动力学的活立木抗风性能分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811322277.0A CN109446691B (zh) 2018-11-08 2018-11-08 基于激光点云与空气动力学的活立木抗风性能分析方法

Publications (2)

Publication Number Publication Date
CN109446691A CN109446691A (zh) 2019-03-08
CN109446691B true CN109446691B (zh) 2022-12-20

Family

ID=65551344

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811322277.0A Active CN109446691B (zh) 2018-11-08 2018-11-08 基于激光点云与空气动力学的活立木抗风性能分析方法

Country Status (1)

Country Link
CN (1) CN109446691B (zh)

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110211133B (zh) * 2019-05-27 2021-01-15 中国农业大学 带叶树木的安全防护策略获取方法、装置与电子设备
CN112052512B (zh) * 2020-07-23 2023-01-10 中国空气动力研究与发展中心计算空气动力研究所 湍流边界层分层判据的方法
CN113378427B (zh) * 2021-05-11 2023-03-10 三峡大学 一种评测乔木枝干抗风载折断性的计算方法
CN114494586B (zh) * 2022-01-10 2024-03-19 南京林业大学 晶格投影的深度学习网络阔叶树枝叶分离与骨架重建方法
CN114467533B (zh) * 2022-01-20 2023-04-07 深圳坤元生态科技有限公司 一种修剪树冠量对树木承受风荷载大小影响的实验方法
CN114600654A (zh) * 2022-01-24 2022-06-10 浙江机电职业技术学院 用于苗木盆景的自动修剪方法、修剪设备和自动修剪系统
CN114648566A (zh) * 2022-03-31 2022-06-21 西安石油大学 一种真实树木对街道污染物运移影响的模拟方法
CN115512121B (zh) * 2022-08-01 2023-05-16 南京林业大学 非完全模拟树木水分养分传输的枝干点云骨架提取方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FI117490B (fi) * 2004-03-15 2006-10-31 Geodeettinen Laitos Menetelmä puustotunnusten määrittämiseksi laserkeilaimen, kuvainformaation ja yksittäisten puiden tulkinnan avulla
CN101887596B (zh) * 2010-06-01 2013-02-13 中国科学院自动化研究所 树木点云数据基于分割和自动生长的三维模型重建方法
CN107705309B (zh) * 2017-10-15 2020-12-04 南京林业大学 激光点云中林木参数评估方法

Also Published As

Publication number Publication date
CN109446691A (zh) 2019-03-08

Similar Documents

Publication Publication Date Title
CN109446691B (zh) 基于激光点云与空气动力学的活立木抗风性能分析方法
CN110570428B (zh) 一种从大规模影像密集匹配点云分割建筑物屋顶面片的方法及系统
Zhou et al. 2.5 d dual contouring: A robust approach to creating building models from aerial lidar point clouds
CN110189412B (zh) 基于激光点云的多楼层室内结构化三维建模方法及系统
CN109509256B (zh) 基于激光雷达的建筑结构自动测量及3d模型生成方法
CN107146280B (zh) 一种基于切分的点云建筑物重建方法
CN106023312B (zh) 基于航空LiDAR数据的三维建筑物模型自动重建方法
CN111986322B (zh) 一种基于结构分析的点云室内场景布局重建方法
CN106651900B (zh) 一种基于轮廓分割的高架原位草莓三维建模方法
CN105787933B (zh) 基于多视角点云配准的岸线三维重建装置及方法
CN110992473B (zh) 一种基于车载激光扫描点云的树木枝干建模方法及系统
CN105844602B (zh) 一种基于体元的机载lidar点云三维滤波方法
CN109685914A (zh) 基于三角网格模型的剖切轮廓自动补面算法
CN111612896B (zh) 一种基于机载激光雷达树点云重建三维树模型的方法
CN106570468A (zh) 一种重建LiDAR原始点云建筑物轮廓线的方法
CN107657659A (zh) 基于长方体拟合扫描三维点云的曼哈顿结构建筑物自动建模方法
CN110794413B (zh) 线性体素分割的激光雷达点云数据电力线检测方法和系统
CN111508015B (zh) 一种基于三维实景数据的建筑物高度提取方法及其装置
CN111667574B (zh) 从倾斜摄影模型自动重建建筑物规则立面三维模型的方法
CN109165476A (zh) 一种模块化风场模型的建模方法及风场模拟方法
CN107545602B (zh) 基于LiDAR点云的空间拓扑关系约束下的建筑物建模方法
CN109493344A (zh) 一种大规模城市三维场景的语义分割方法
CN113066162A (zh) 一种用于电磁计算的城市环境快速建模方法
CN116883754A (zh) 一种面向地面LiDAR点云的建筑物信息提取方法
CN109215112A (zh) 一种单侧点云模型的标注方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230921

Address after: 210000 No.8 Huayuan Road, Xuanwu District, Nanjing, Jiangsu Province

Patentee after: Nanjing Maoting Information Technology Co.,Ltd.

Address before: Nanjing City, Jiangsu province 210037 Longpan Road No. 159

Patentee before: NANJING FORESTRY University