Forest tree diameter at breast height and tree height measurement method based on double-layer unmanned aerial vehicle laser radar point cloud fusion in complex under-forest environment

Abstract

The invention relates to a forest tree diameter at breast height and tree height measurement method based on double-layer unmanned aerial vehicle laser radar point cloud fusion in a complex under-forest environment, and belongs to the technical field of forestry measurement. According to the method, multi-source point cloud registration fusion is carried out based on feature point matching and a weighted fusion algorithm; generating a ground DEM and carrying out point cloud height normalization; carrying out individual tree segmentation by adopting an improved mark control watershed algorithm; constructing a multi-dimensional feature vector for each tree, and calculating each feature; extracting a tree height based on the fused point cloud; an effective point cloud coverage index is introduced to judge the effectiveness of a 1.3 m height point cloud, weighted least square circle fitting is adopted to directly extract the diameter at breast height, and a pre-trained multi-feature fusion gradient lifting model is adopted to estimate the diameter at breast height for a tree with point cloud missing. According to the method, complete stand three-dimensional data is obtained through multi-source point cloud fusion, a multi-dimensional feature system is constructed, a complete mathematical calculation framework is established, and the technical problem that the diameter at breast height cannot be accurately measured due to complex under-forest shrub and grass shielding is effectively solved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of forestry survey, and particularly relates to a method for measuring the diameter at breast height and tree height of a forest tree in a complex understory environment based on fusion of point clouds of dual-layer unmanned aerial vehicles (UAVs) and laser radars. BACKGROUND

[0002] The diameter at breast height (the diameter of a tree at 1.3 m above the ground, denoted as DBH) and the tree height (denoted as H) are the most important tree measurement factors in forest resource investigation and are basic data for calculating the volume, biomass and carbon storage of forest trees. Accurate acquisition of the diameter at breast height and the tree height of a stand is of great significance for forest resource management, ecological assessment and forestry production.

[0003] Traditional measurement of the diameter at breast height and the tree height mainly relies on manual field measurement using a girth tape, an altimeter and other tools. This method is not only time-consuming and labor-intensive, but also has low efficiency and high safety risks in complex terrain conditions. With the development of laser radar technology and unmanned aerial vehicle technology, the use of unmanned aerial vehicles carrying laser radar scanners for forest parameter extraction has become a research hotspot.

[0004] Currently, the following technical problems exist in the measurement of unmanned aerial vehicle laser radars:

[0005] (1) Data missing problem of single flight mode:

[0006] Flying above the canopy can obtain complete crown information and tree height data, but it is difficult to obtain point cloud data of the tree trunk under the canopy due to the canopy shelter. Low-altitude flight in the forest can obtain trunk point cloud data, but due to the flight height limitation, it is difficult to obtain complete crown and tree height information. Both modes have the problem of incomplete data.

[0007] (2) Complex understory shrub and grass sheltering problem:

[0008] In actual applications, there are dense shrub layers and herb layers under many plantations or natural forests, and the height is often between 0.5-2.0 m, which seriously blocks the position at the height of 1.3 m of the tree trunk. Even if low-altitude flight in the forest is used, it is difficult for the laser radar to directly obtain effective trunk point cloud data at this height.

[0009] (3) Limited accuracy of diameter at breast height estimation method:

[0010] Existing diameter at breast height estimation methods are mainly based on simple crown diameter-diameter at breast height or tree height-diameter at breast height binary regression relationships, such as:

[0011] (1)

[0012] (2)

[0013] This method utilizes only one or two tree features, failing to fully extract the rich three-dimensional structural information contained in point cloud data. The coefficient of determination R of the regression model is... 2 The value is usually between 0.6 and 0.75, which is insufficient to meet the needs of precision forestry.

[0014] (4) Lack of quantitative criteria for judging the validity of point clouds:

[0015] Existing methods lack clear quantitative criteria for determining whether point clouds at a height of 1.3m are sufficient for direct fitting of diameter at breast height. They often rely on human experience to make judgments, which affects the automation level of the measurement process and the consistency of the results.

[0016] Therefore, there is an urgent need for a comprehensive measurement method that can integrate multi-source point cloud data, establish a complete mathematical calculation framework, quantitatively determine the validity of point clouds, and fully mine multi-dimensional feature information. Summary of the Invention

[0017] The purpose of this invention is to overcome the shortcomings of the prior art and provide a method for measuring the diameter at breast height (DBH) and tree height of trees in complex forest environments based on the fusion of point clouds from dual-layer UAV lidar. It establishes a complete mathematical calculation framework from point cloud acquisition, registration and fusion, feature extraction to DBH measurement, effectively solving the technical problem of inaccurate DBH measurement caused by the obstruction of shrubs and grasses in complex forest environments.

[0018] To achieve the above objectives, the technical solution of the present invention is: a method for measuring tree diameter at breast height (DBH) and tree height in complex forest environments based on dual-layer UAV lidar point cloud fusion, comprising:

[0019] UAVs were used to collect point cloud data in and above the forest, and multi-source point cloud registration and fusion were performed based on feature point matching and weighted fusion algorithms.

[0020] Generate a ground DEM and perform point cloud height normalization;

[0021] An improved label-controlled watershed algorithm is used for single-tree segmentation;

[0022] For each tree, a multi-dimensional feature vector is constructed to calculate the tree crown geometric features, tree crown structural features, visible trunk features, point cloud distribution features, and neighborhood competition features.

[0023] Tree height extraction based on fused point cloud;

[0024] An effective point cloud coverage index is introduced to determine the effectiveness of point clouds at a height of 1.3m. For trees with effective point clouds, weighted least squares circle fitting is used to directly extract the diameter at breast height (DBH). For trees with missing point clouds, a pre-trained multi-feature fusion machine learning model is used to estimate the DBH.

[0025] Furthermore, drones are used to collect point cloud data while flying through the forest. Specifically, the method is as follows:

[0026] Drones equipped with LiDAR scanners were used to conduct low-altitude flights within the target forest stand, employing an AI visual obstacle avoidance system to achieve safe and autonomous flight in complex forest environments. The LiDAR performed 360-degree scans, using a multi-angle, multi-flight strategy to collect point cloud data from tree trunks in various directions, thus acquiring a forest point cloud dataset. , where the i-th point Includes three-dimensional coordinates (x) i y i , z i and reflection intensity information I i .

[0027] Furthermore, drones are used to collect point cloud data while flying over the forest. The specific method is as follows:

[0028] A drone equipped with a lidar scanner was used to fly over the canopy of the target forest stand, performing round-trip scans along a regular flight path. The spacing between the flight paths was calculated based on the lidar's field of view and flight altitude to ensure complete coverage of the canopy point cloud and acquire a forest point cloud dataset. ,in, This represents the j-th point data in the forest point cloud dataset.

[0029] Furthermore, multi-source point cloud registration and fusion are carried out in the following ways:

[0030] (1) Preprocess the forest point cloud dataset Pu and the forest point cloud dataset Po respectively, including denoising, filtering and coordinate system 1;

[0031] (2) Extract the same feature points from the two sets of point cloud datasets as registration control points. The feature points include tree top feature points, tree trunk branching points, or pre-deployed artificial target points.

[0032] (3) The optimal rigid transformation is solved using the ICP iterative nearest point algorithm to minimize the root mean square error after registration. Let the forest point cloud dataset be transformed as follows: The registration error is defined as:

[0033]

[0034] in, and R represents the corresponding points in the forest point cloud dataset and the transformed forest point cloud dataset, respectively. R is the rotation matrix, T is the translation vector, and n is the number of corresponding points. Iterative optimization is performed until the RMSE converges to below the set threshold.

[0035] (4) Adaptive weighted fusion is performed on the registered point clouds. In the overlapping region, the fusion weight is calculated based on the point cloud quality factor. The formula for calculating the point cloud quality factor Q is:

[0036]

[0037] Where ρ is the local point cloud density, I is the average reflection intensity, and σ z Let be the standard deviation in the height direction, and the corresponding parameter with the subscript max represents the maximum value of the corresponding parameter. α, β, and γ are weighting coefficients and satisfy α + β + γ = 1.

[0038] Calculate fusion weights based on point cloud quality factors:

[0039]

[0040] in, The fusion weight for cloud point analysis in the forest. The fusion weight for "Clouds Above the Forest" The quality factor for forest point clouds. The quality factor for clouds in the forest.

[0041] Fusion point cloud density The calculation formula is:

[0042] .

[0043] in, The density of the merged point cloud. , These represent the original densities of point clouds in the forest and above the forest at this location, respectively.

[0044] Furthermore, a ground DEM is generated and point cloud height is normalized, specifically as follows:

[0045] (1) The cloth simulation filtering algorithm CSF is used to extract the ground point cloud after the multi-source point cloud registration and fusion. The virtual cloth is dropped from above to simulate the process of the cloth sticking to the terrain, so as to realize the separation of ground point cloud and non-ground point cloud.

[0046] (2) Using ground point clouds, a continuous ground digital elevation model (DEM) is generated by kriging interpolation.

[0047] (3) Perform height normalization on all point clouds, with the normalized height h i The calculation is as follows:

[0048]

[0049] Among them, z iis the original elevation of the point, and DEM(xi, yi) is the DEM elevation value corresponding to the point location.

[0050] Furthermore, an improved label-controlled watershed algorithm is used for single-tree segmentation, specifically as follows:

[0051] (1) Generate a canopy height model CHM based on normalized point cloud, where CHM is the maximum normalized height of the point cloud within the corresponding grid;

[0052] (2) Gaussian smoothing is applied to the CHM, and the local maximum filtering method is used to detect the canopy apex; let the window size be W, if a pixel value is greater than all pixel values ​​in its W × W neighborhood, it is marked as a candidate point for the canopy apex; the window size is adaptively determined according to the average canopy width of the stand:

[0053]

[0054] in, is an estimate of the average crown width of the stand, and k is a proportionality coefficient;

[0055] (3) Using the detected tree crown vertices as seed points, the label control watershed algorithm is used to segment individual trees and obtain an independent point cloud subset for each tree.

[0056] Furthermore, a multi-dimensional feature vector is constructed for each tree to calculate the tree crown geometric features, tree crown structural features, visible trunk features, point cloud distribution features, and neighborhood competition features. The specific method is as follows:

[0057] For each segmented tree, 23 feature parameters across five categories are extracted from its point cloud data to construct a feature vector. ;

[0058] (1) Geometric characteristics of the tree crown:

[0059] Crown width CW:

[0060]

[0061] Among them, CW EW CW is the maximum width in the east-west direction. NS This represents the maximum width in the north-south direction.

[0062] The canopy projection area CA is calculated using a concave hull algorithm, with the concave hull parameter α adaptively determined based on the point cloud density.

[0063]

[0064] in, is the average spacing of the point cloud, and k is an empirical coefficient;

[0065] Crown volume CV:

[0066]

[0067] Where, N v s represents the number of voxels containing the point cloud, and s represents the voxel side length.

[0068] Canopy Surface Area (CSA): The three-dimensional convex hull surface area is calculated using the Alpha-Shape algorithm.

[0069] Crown length CL:

[0070]

[0071] Among them, H max H is the height of the highest point of the tree crown. cb The height below the branch is defined as the first time the accumulated point cloud density from the base of the tree upwards exceeds the density threshold ρ. th Height;

[0072] Coronal Index (CSI):

[0073]

[0074] (2) Tree crown structure characteristics:

[0075] Canopy density CD:

[0076]

[0077] Where, N c CV represents the total number of points in the canopy point cloud;

[0078] Vertical Distribution Ratio (VDR) of Tree Canopy Point Cloud: Dividing the tree canopy into m equal layers, the VDR of the point cloud in the j-th layer. j for:

[0079]

[0080] Where, N j Let j be the number of point clouds in the j-th layer;

[0081] Leaf Area Index (LAI): Based on point cloud penetration, a modified Beer-Lambert model is used.

[0082]

[0083] Where θ is the incident zenith angle of the laser, and k ext P is the extinction coefficient. gap The gap ratio is calculated using the following formula:

[0084]

[0085] Where, N ground N represents the number of laser pulses that penetrate the canopy and reach the ground. total This represents the total number of laser pulses entering the canopy.

[0086] Canopy porosity (GP):

[0087]

[0088] Where, N gap N represents the number of grid cells without point cloud coverage. grid This represents the total number of grid cells.

[0089] Crown Complexity Index (CCI): Defined as the ratio of the surface area of ​​a tree crown to the surface area of ​​a sphere of equal volume.

[0090]

[0091] (3) Visible features of the tree trunk:

[0092] The cross-sectional point cloud at height h is extracted from the tree trunk point cloud, and a circle is fitted using the weighted least squares method with diameter D. h The calculation is obtained by minimizing the objective function:

[0093] Among them, (x i y i (x) represents the point cloud coordinates of the cross section. c y c () represents the coordinates of the center of the circle, r represents the radius, and w represents the coordinates of the center of the circle. i Let D be the point weight. h = 2r;

[0094] Point weight w i Adaptive calculation based on the distance from the point to the fitted circle:

[0095]

[0096] Where, d i Let σ be the distance from point i to the circle in the current iteration. d This represents the distance from the standard deviation.

[0097] The diameter D at the lowest point is visible. L Searching for points in the tree trunk cloud at heights above the shrub layer (h). s minimum effective height h L Diameter at:

[0098]

[0099] Where, N h N represents the number of point clouds at height h.th The effective point threshold;

[0100] The diameter D at the highest point is visible. H Extract the diameter at the boundary between the trunk and the crown;

[0101] Visible length of tree trunk VL:

[0102]

[0103] Among them, h H h is the height of the highest visible point. L This is the lowest visible height.

[0104] Trunk taper T:

[0105]

[0106] (4) Point cloud distribution characteristics:

[0107] Point cloud height quantile H q Calculate the qth percentile of the point cloud height for the entire tree, where q takes values ​​of 25, 50, 75, 90, and 95.

[0108]

[0109] Where H(k) represents the k-th value in the height sequence, and n is the total number of point clouds;

[0110] Point cloud height variation coefficient HCV:

[0111]

[0112] Where, σ H The standard deviation of point cloud height. The mean height of the point cloud;

[0113] Point cloud height skewness H skew :

[0114]

[0115] Point cloud height kurtosis H kurt :

[0116]

[0117] Point cloud intensity statistical characteristics:

[0118] Mean Intensity :

[0119]

[0120] Strength standard deviation σI :

[0121]

[0122] Intensity coefficient of variation (ICV):

[0123]

[0124] (5) Neighborhood competition characteristics:

[0125] Number of neighboring trees NN: Statistical based on the target tree's location with a radius R. n Number of adjacent trees within the range:

[0126]

[0127] Where N is the total number of trees in the stand excluding the target tree, j is the tree traversal index, and dist j Let be the horizontal distance between the target tree and the j-th tree, 1(·) be the indicator function, and R be the horizontal distance between the target tree and the j-th tree. n = 5m;

[0128] Average Neighbor Distance AND: Calculates the average distance between the target tree and its k nearest neighboring trees.

[0129]

[0130] Among them, dist (j) This represents the distance to the j-th nearest neighboring tree after sorting by distance.

[0131] Hegyi Competition Index (CI):

[0132]

[0133] Among them, D i For the target tree's diameter at breast height (DBH) or crown width, D j For the diameter at breast height (DBH) or crown width of the j-th neighboring tree, dist ij The horizontal distance between the two trees;

[0134] Stand Density Index (SDI): The Reineke stand density index formula is used.

[0135]

[0136] Where N is the number of plants per unit area. D represents the average diameter at breast height (DBH) of the forest stand. ref For reference chest diameter.

[0137] Furthermore, tree height is extracted based on the fused point cloud, specifically in the following manner:

[0138] Based on the fused point cloud, the highest point of the point cloud for each tree is extracted. To eliminate the influence of outliers, the 99th percentile height is used as a robust estimate of the tree height.

[0139]

[0140] in, H represents the ground elevation at the base of the tree, (x0, y0) represents the coordinates of the tree base, and H represents the ground elevation at the base of the tree. (0 .99n) This is the 99th percentile value of the point cloud height.

[0141] Furthermore, an effective point cloud coverage index is introduced to determine the effectiveness of point clouds at a height of 1.3m. For trees with effective point clouds, weighted least squares circle fitting is used to directly extract the diameter at breast height (DBH). For trees with missing point clouds, a pre-trained multi-feature fusion machine learning model is used to estimate the DBH. The specific method is as follows:

[0142] (1) Calculate the Effective Point Cloud Cover Index (EPCI):

[0143] A sampling layer with a thickness of ∆h = 0.1m was set at a height of 1.3m, and the point cloud within the corresponding layer was extracted; the 360 ​​degrees were divided into 12 sectors, each sector being 30 degrees, and the angular coverage rate (ACR) and spatial distribution uniformity (SDU) were calculated.

[0144]

[0145] Where, N total sector = 12, N sector The number of sectors containing point clouds;

[0146]

[0147] Where, σ sector The standard deviation of the point cloud count for each sector. Let ϵ be the average number of point clouds in each sector, and let ϵ be a small constant to prevent division by zero.

[0148] The formula for calculating the Effective Point Cloud Cover Index (EPCI) is as follows:

[0149]

[0150] Where, N 1.3 N represents the number of point clouds within a sampling layer at a height of 1.3m. th ω1, ω2, and ω3 are the threshold values ​​for the number of points in the cloud, and ω1 + ω2 + ω3 = 1.

[0151] (2) Determination and execution of chest diameter measurement:

[0152] When EPCI ≥ EPCI th(EPCI stands for Effective Point Cloud Coverage Index) th (Assuming an effective point cloud coverage index threshold), the diameter at breast height (DBH) is directly extracted using a robust circle fitting method based on RANSAC.

[0153] Three points were randomly selected from the point cloud at a height of 1.3m to fit an initial circle;

[0154] Calculate the distance from all points to the initial circle, count the number of interior points, and determine the interior point as follows:

[0155]

[0156] Where, d i Let be the distance from point i to the center of the circle, r be the radius of the circle, and ϵ be the distance threshold;

[0157] After 100 iterations, the circle with the most interior points is selected as the optimal fitting result.

[0158] The final chest diameter was obtained by fine-fitting the inlier points using the weighted least squares method.

[0159]

[0160] Calculate the fit confidence score FC:

[0161]

[0162] Where, N inlier RMSE is the number of interior points. fit The root mean square error is the fitted value.

[0163] When EPCI < EPCI th At that time, a pre-trained multi-feature fusion machine learning model was used to estimate the diameter at breast height.

[0164] Furthermore, the multi-feature fusion machine learning model is as follows:

[0165] (1) Feature selection: A feature selection method based on mutual information is adopted to calculate the F of each feature. i Mutual information MI between diameter at breast height (DBH) and sternum (MI):

[0166]

[0167] in, Let f be the i-th feature variable, and f be the feature value. The value of d is the value of the diameter at breast height (DBH). Joint probability distribution Marginal probabilities of features Marginal probability of diaphragm diameter, screening mutual information value greater than threshold MI th Features;

[0168] (2) Collinearity diagnosis: Calculate the variance inflation factor (VIF) between features:

[0169]

[0170] in, The regression determination coefficient is the i-th feature as the dependent variable and the remaining features as independent variables, after removing redundant features with VIF>10.

[0171] (3) Model training: Gradient boosting decision tree algorithm (GBDT) is used, with Huber loss as the loss function.

[0172]

[0173] Where y is the measured chest diameter. To predict chest diameter, δ is a threshold parameter; the parameter is optimized through grid search and 5-fold cross-validation.

[0174] (4) Model evaluation indicators:

[0175] Coefficient of determination R 2 :

[0176]

[0177] Root Mean Square Error (RMSE):

[0178]

[0179] Mean Absolute Percentage Error (MAPE):

[0180]

[0181] The final model must satisfy R² > 0.85 and MAPE < 10%.

[0182] Compared with the prior art, the present invention has the following beneficial effects:

[0183] (1) Establishing a complete mathematical calculation framework: This invention establishes clear mathematical formulas and calculation methods for each processing step, from error assessment of point cloud registration and calculation of fusion weights to the definition of each parameter in feature extraction, and then to the validity judgment and fitting algorithm of diaphragm measurement, forming a complete and reproducible technical solution.

[0184] (2) Propose an effective point cloud coverage index: Innovatively propose the EPCI index, which quantitatively integrates three dimensions: angle coverage, spatial distribution uniformity and point cloud quantity. This provides an objective basis for automatically judging whether a 1.3m height point cloud is suitable for direct fitting and improves the automation level of the measurement process.

[0185] (3) Adaptive weighted fusion of multi-source point clouds: An adaptive weight calculation method based on point cloud quality factor Q is established so that different height layers can automatically select higher quality data sources and give full play to the complementary advantages of forest point clouds and forest top point clouds.

[0186] (4) Multidimensional feature system and quantitative calculation: Construct a feature system containing 23 parameters in 5 categories. Each feature is given a clear mathematical definition and calculation formula to fully explore the geometric, structural and competitive information contained in the point cloud.

[0187] (5) Robust method for measuring diameter at breast height: The circular fitting method combining RANSAC and weighted least squares is adopted, which has strong robustness against outliers and noise. At the same time, the fitting confidence index FC is introduced to provide a quantitative assessment of the reliability of the measurement results. Attached Figure Description

[0188] Figure 1 This is a flowchart illustrating the overall technical process of the method of the present invention.

[0189] Figure 2 This is a schematic diagram of a two-layer drone flight for data collection, demonstrating the data collection strategy of low-altitude flight in the forest and high-altitude flight over the forest.

[0190] Figure 3 This is a schematic diagram illustrating the principle of multi-source point cloud registration and fusion, showing the feature point matching and adaptive weighted fusion process.

[0191] Figure 4 This diagram illustrates the calculation of the Effective Point Cloud Coverage Index (EPCI), showcasing sector division and angle coverage calculation methods.

[0192] Figure 5 This is a schematic diagram of multidimensional feature extraction, illustrating the spatial meaning of the five major feature categories.

[0193] Figure 6 This is a flowchart of the thorax measurement process, illustrating the judgment and execution procedures based on EPCI. Detailed Implementation

[0194] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings.

[0195] This invention provides a method for measuring tree diameter at breast height (DBH) and tree height in complex forest environments based on dual-layer UAV lidar point cloud fusion, comprising:

[0196] UAVs were used to collect point cloud data in and above the forest, and multi-source point cloud registration and fusion were performed based on feature point matching and weighted fusion algorithms.

[0197] Generate a ground DEM and perform point cloud height normalization;

[0198] An improved label-controlled watershed algorithm is used for single-tree segmentation;

[0199] For each tree, a multi-dimensional feature vector is constructed to calculate the tree crown geometric features, tree crown structural features, visible trunk features, point cloud distribution features, and neighborhood competition features.

[0200] Tree height extraction based on fused point cloud;

[0201] An effective point cloud coverage index is introduced to determine the effectiveness of point clouds at a height of 1.3m. For trees with effective point clouds, weighted least squares circle fitting is used to directly extract the diameter at breast height (DBH). For trees with missing point clouds, a pre-trained multi-feature fusion machine learning model is used to estimate the DBH.

[0202] The following is a detailed implementation process of the present invention.

[0203] like Figures 1-6 As shown, the present invention provides a method for measuring tree diameter at breast height (DBH) and tree height in complex forest environments based on dual-layer UAV lidar point cloud fusion, comprising the following steps:

[0204] Step S1: Forest point cloud data collection:

[0205] Small drones equipped with LiDAR scanners were used to conduct low-altitude flights within the target forest stand, employing an AI visual obstacle avoidance system to achieve safe and autonomous flight in complex forest environments. The LiDAR performed 360-degree scans, using a multi-angle, multi-flight strategy to collect point cloud data from tree trunks in various directions, thus acquiring a forest point cloud dataset. , where the i-th point Includes three-dimensional coordinates (x) i y i , z i and reflection intensity information I i .

[0206] Step S2, Forest Point Cloud Data Acquisition:

[0207] A drone equipped with a lidar scanner was used to fly over the canopy of the target forest stand, performing round-trip scans along a regular flight path. The spacing between the flight paths was calculated based on the lidar's field of view and flight altitude to ensure complete coverage and appropriate overlap of the canopy point cloud, thereby acquiring a forest point cloud dataset. ,in, This represents the j-th point data in the forest point cloud dataset.

[0208] Step S3: Multi-source point cloud registration and fusion:

[0209] (1) Preprocess the forest point cloud dataset Pu and the forest point cloud dataset Po respectively, including denoising, filtering and coordinate system 1;

[0210] (2) Extract the same feature points from the two sets of point cloud datasets as registration control points. The feature points include tree top feature points, tree trunk branching points, or pre-deployed artificial target points.

[0211] (3) The optimal rigid transformation is solved using the ICP iterative nearest point algorithm to minimize the root mean square error after registration. Let the forest point cloud dataset be transformed as follows: The registration error is defined as:

[0212]

[0213] in, and R represents the corresponding points in the forest point cloud dataset and the transformed forest point cloud dataset, respectively. R is the rotation matrix, T is the translation vector, and n is the number of corresponding points. Iterative optimization is performed until the RMSE converges to below the set threshold.

[0214] (4) Adaptive weighted fusion is performed on the registered point clouds. In the overlapping region, the fusion weight is calculated based on the point cloud quality factor. The formula for calculating the point cloud quality factor Q is:

[0215]

[0216] Where ρ is the local point cloud density, I is the average reflection intensity, and σ z The standard deviation in the height direction (characterizing the dispersion of the point cloud) is denoted by 'max'. The corresponding parameter with the subscript 'max' represents the maximum value of the corresponding parameter. α, β, and γ are weighting coefficients and satisfy α + β + γ = 1. Recommended values ​​are α = 0.4, β = 0.3, and γ = 0.3.

[0217] Calculate fusion weights based on point cloud quality factors:

[0218]

[0219] in, The fusion weight for cloud point analysis in the forest. The fusion weight for "Clouds Above the Forest" The quality factor for forest point clouds. The quality factor for clouds in the forest.

[0220] In the understory (e.g., at a height of 0-3m), the quality of point clouds in the forest is generally better than that in the forest above, therefore w u > w o In the canopy area, the quality of overstory point clouds is better, therefore w o > w u .

[0221] Fusion point cloud density The calculation formula is:

[0222]

[0223] in, The density of the merged point cloud. , These represent the original densities of point clouds in the forest and above the forest at this location, respectively.

[0224] Step S4: Ground DEM generation and point cloud normalization:

[0225] (1) The cloth simulation filtering algorithm CSF is used to extract the ground point cloud after the multi-source point cloud registration and fusion. The virtual cloth is dropped from above to simulate the process of the cloth sticking to the terrain, so as to realize the separation of ground point cloud and non-ground point cloud.

[0226] (2) Using ground point cloud, a continuous ground digital elevation model (DEM) is generated by Kriging interpolation method, with a resolution of 0.5m;

[0227] (3) Perform height normalization on all point clouds, with the normalized height h i The calculation is as follows:

[0228]

[0229] Among them, z i is the original elevation of the point, and DEM(xi, yi) is the DEM elevation value corresponding to the point location.

[0230] Step S5, Single-log splitting:

[0231] (1) Generate a canopy height model CHM based on normalized point cloud, where CHM is the maximum normalized height of the point cloud within the corresponding grid;

[0232] (2) Gaussian smoothing is applied to the CHM, and the local maximum filtering method is used to detect the canopy apex; let the window size be W, if a pixel value is greater than all pixel values ​​in its W × W neighborhood, it is marked as a candidate point for the canopy apex; the window size is adaptively determined according to the average canopy width of the stand:

[0233]

[0234] in, This is an estimate of the average crown width of the forest stand, where k is a proportionality coefficient ranging from 0.6 to 0.8.

[0235] (3) Using the detected tree crown vertices as seed points, the label control watershed algorithm is used to segment individual trees and obtain an independent point cloud subset for each tree.

[0236] Step S6, Multidimensional Feature Extraction:

[0237] For each segmented tree, 23 feature parameters across five categories are extracted from its point cloud data to construct a feature vector. ;

[0238] (1) Geometric characteristics of the tree crown:

[0239] Crown width CW:

[0240]

[0241] Among them, CW EW CW is the maximum width in the east-west direction. NS This represents the maximum width in the north-south direction.

[0242] The canopy projection area CA is calculated using a concave hull algorithm, with the concave hull parameter α adaptively determined based on the point cloud density.

[0243]

[0244] in, is the average spacing of the point cloud, and k is an empirical coefficient;

[0245] Crown volume CV:

[0246]

[0247] Where, N v The number of voxels containing the point cloud is s, where s is the side length of the voxel and is taken as 0.1m.

[0248] Canopy Surface Area (CSA): The three-dimensional convex hull surface area is calculated using the Alpha-Shape algorithm.

[0249] Crown length CL:

[0250]

[0251] Among them, H max H is the height of the highest point of the tree crown. cb The height below the branch is defined as the first time the accumulated point cloud density from the base of the tree upwards exceeds the density threshold ρ. th Height;

[0252] Coronal Index (CSI):

[0253]

[0254] (2) Tree crown structure characteristics:

[0255] Canopy density CD:

[0256]

[0257] Where, N c CV represents the total number of points in the canopy point cloud;

[0258] Vertical Distribution Ratio (VDR) of Tree Canopy Point Cloud: Dividing the tree canopy into m equal layers, the VDR of the point cloud in the j-th layer. j for:

[0259]

[0260] Where, N j Let j be the number of point clouds in the j-th layer;

[0261] Leaf Area Index (LAI): Based on point cloud penetration, a modified Beer-Lambert model is used.

[0262]

[0263] Where θ is the incident zenith angle of the laser, and k ext P is the extinction coefficient, with a value of 0.5. gap The gap ratio is calculated using the following formula:

[0264]

[0265] Where, N ground N represents the number of laser pulses that penetrate the canopy and reach the ground. total This represents the total number of laser pulses entering the canopy.

[0266] Canopy porosity (GP):

[0267]

[0268] Where, N gap N represents the number of grid cells without point cloud coverage. grid This represents the total number of grid cells.

[0269] Crown Complexity Index (CCI): Defined as the ratio of the surface area of ​​a tree crown to the surface area of ​​a sphere of equal volume.

[0270]

[0271] (3) Visible features of the tree trunk:

[0272] The cross-sectional point cloud at height h is extracted from the tree trunk point cloud, and a circle is fitted using the weighted least squares method with diameter D. h The calculation is obtained by minimizing the objective function:

[0273]

[0274] Where (xi, yi) are the point cloud coordinates of the cross section, (xc, yc) are the coordinates of the center of the circle, r is the radius, and w i Let D be the point weight. h = 2r;

[0275] Point weight w i Adaptive calculation based on the distance from the point to the fitted circle:

[0276]

[0277] Where, d i Let σ be the distance from point i to the circle in the current iteration. d This represents the distance from the standard deviation.

[0278] The diameter D at the lowest point is visible. L Searching for points in the tree trunk cloud at heights above the shrub layer (h). s minimum effective height h L Diameter at:

[0279]

[0280] Where, N h N represents the number of point clouds at height h. th The effective point threshold;

[0281] The diameter D at the highest point is visible. H Extract the diameter at the boundary between the trunk and the crown;

[0282] Visible length of tree trunk VL:

[0283]

[0284] Among them, h H h is the height of the highest visible point. L This is the lowest visible height.

[0285] Trunk taper T:

[0286]

[0287] (4) Point cloud distribution characteristics:

[0288] Point cloud height quantile H qCalculate the qth percentile of the point cloud height for the entire tree, where q takes values ​​of 25, 50, 75, 90, and 95.

[0289]

[0290] Among them, H (k) This represents the k-th value in the height sequence, where n is the total number of point clouds;

[0291] Point cloud height variation coefficient HCV:

[0292]

[0293] Where, σ H The standard deviation of point cloud height. The mean height of the point cloud;

[0294] Point cloud height skewness H skew :

[0295]

[0296] Point cloud height kurtosis H kurt :

[0297]

[0298] Point cloud intensity statistical characteristics:

[0299] Mean Intensity :

[0300]

[0301] Strength standard deviation σ I :

[0302]

[0303] Intensity coefficient of variation (ICV):

[0304]

[0305] (5) Neighborhood competition characteristics:

[0306] NN (Nearest Trees): Counts the number of neighboring trees within a radius Rn centered on the target tree's location.

[0307]

[0308] Where N is the total number of trees in the stand excluding the target tree, j is the tree traversal index, and dist j Let be the horizontal distance between the target tree and the j-th tree, 1(·) be the indicator function, and R be the horizontal distance between the target tree and the j-th tree. n= 5m;

[0309] Average Neighbor Distance AND: Calculates the average distance between the target tree and its k nearest neighboring trees.

[0310]

[0311] Among them, dist (j) This represents the distance to the j-th nearest neighboring tree after sorting by distance.

[0312] Hegyi Competition Index (CI):

[0313]

[0314] Among them, D i For the target tree's diameter at breast height (DBH) or crown width, D j For the diameter at breast height (DBH) or crown width of the j-th neighboring tree, dist ij The horizontal distance between the two trees;

[0315] Stand Density Index (SDI): The Reineke stand density index formula is used.

[0316]

[0317] Where N is the number of plants per unit area. D represents the average diameter at breast height (DBH) of the forest stand. ref For reference, the diameter at breast height is usually taken as 25cm.

[0318] Furthermore, tree height is extracted based on the fused point cloud, specifically in the following manner:

[0319] Based on the fused point cloud, the highest point of the point cloud for each tree is extracted. To eliminate the influence of outliers, the 99th percentile height is used as a robust estimate of the tree height.

[0320]

[0321] in, H represents the ground elevation at the base of the tree, (x0, y0) represents the coordinates of the tree base, and H represents the ground elevation at the base of the tree. (0 .99n) This is the 99th percentile value of the point cloud height.

[0322] Furthermore, an effective point cloud coverage index is introduced to determine the effectiveness of point clouds at a height of 1.3m. For trees with effective point clouds, weighted least squares circle fitting is used to directly extract the diameter at breast height (DBH). For trees with missing point clouds, a pre-trained multi-feature fusion machine learning model is used to estimate the DBH. The specific method is as follows:

[0323] (1) Calculate the Effective Point Cloud Cover Index (EPCI):

[0324] A sampling layer with a thickness of ∆h = 0.1m was set at a height of 1.3m, and the point cloud within the corresponding layer was extracted; the 360 ​​degrees were divided into 12 sectors, each sector being 30 degrees, and the angular coverage rate (ACR) and spatial distribution uniformity (SDU) were calculated.

[0325]

[0326] Where, N total sector = 12, N sector The number of sectors containing point clouds;

[0327]

[0328] Where, σ sector The standard deviation of the point cloud count for each sector. Let ϵ be the average number of point clouds in each sector, and let ϵ be a small constant to prevent division by zero.

[0329] The formula for calculating the Effective Point Cloud Cover Index (EPCI) is as follows:

[0330]

[0331] Where, N 1.3 N represents the number of point clouds within a sampling layer at a height of 1.3m. th The threshold for the number of points in the cloud is ω1=0.4, ω2=0.3, and ω3=0.3, which are weighting coefficients and ω1 + ω2 + ω3 = 1.

[0332] (2) Determination and execution of chest diameter measurement:

[0333] When EPCI ≥ EPCI th (When the threshold is 0.6) (EPCI is the Effective Point Cloud Coverage Index, EPCI) th (Assuming an effective point cloud coverage index threshold), the diameter at breast height (DBH) is directly extracted using a robust circle fitting method based on RANSAC.

[0334] Three points were randomly selected from the point cloud at a height of 1.3m to fit an initial circle;

[0335] Calculate the distance from all points to the initial circle, count the number of interior points, and determine the interior point as follows:

[0336]

[0337] Where, d i Let be the distance from point i to the center of the circle, r be the radius of the circle, and ϵ = 0.02m be the distance threshold.

[0338] After 100 iterations, the circle with the most interior points is selected as the optimal fitting result.

[0339] The final chest diameter was obtained by fine-fitting the inlier points using the weighted least squares method.

[0340]

[0341] Calculate the fit confidence score FC:

[0342]

[0343] Where, N inlier RMSE is the number of interior points. fit The root mean square error is the fitted value.

[0344] When EPCI < EPCI th At that time, a pre-trained multi-feature fusion machine learning model was used to estimate the diameter at breast height.

[0345] Furthermore, the multi-feature fusion machine learning model is as follows:

[0346] (1) Feature selection: A feature selection method based on mutual information is adopted to calculate the F of each feature. i Mutual information MI between diameter at breast height (DBH) and sternum (MI):

[0347]

[0348] in, Let f be the i-th feature variable, and f be the feature value. The value of d is the value of the diameter at breast height (DBH). Joint probability distribution Marginal probabilities of features Marginal probability of diaphragm diameter, screening mutual information value greater than threshold MI th Features;

[0349] (2) Collinearity diagnosis: Calculate the variance inflation factor (VIF) between features:

[0350]

[0351] in, The regression determination coefficient is the i-th feature as the dependent variable and the remaining features as independent variables, after removing redundant features with VIF>10.

[0352] (3) Model training: Gradient boosting decision tree algorithm (GBDT) is used, with Huber loss as the loss function.

[0353]

[0354] Where y is the measured chest diameter. To predict chest diameter, δ is a threshold parameter; the parameter is optimized through grid search and 5-fold cross-validation.

[0355] (4) Model evaluation indicators:

[0356] Coefficient of determination R 2 :

[0357]

[0358] Root Mean Square Error (RMSE):

[0359]

[0360] Mean Absolute Percentage Error (MAPE):

[0361]

[0362] The final model must satisfy R 2 > 0.85 and MAPE < 10%.

[0363] The above are preferred embodiments of the present invention. Any changes made to the technical solution of the present invention that do not exceed the scope of the technical solution of the present invention shall fall within the protection scope of the present invention.

Claims

1. A method for measuring the diameter at breast height and tree height of a forest tree in a complex understory environment based on the fusion of dual-layer unmanned aerial vehicle LiDAR point clouds, characterized by, The application comprises the following steps: Respectively using unmanned aerial vehicle to collect point cloud data in forest and on forest, and using feature point matching and weighted fusion algorithm to perform multi-source point cloud registration and fusion; Generating ground DEM and performing point cloud height normalization; Using improved marker-controlled watershed algorithm to perform single tree segmentation; Constructing multi-dimensional feature vector for each tree, calculating crown geometric feature, crown structure feature, visible part of stem feature, point cloud distribution feature and neighborhood competition feature; Extracting tree height based on fused point cloud; Introducing effective point cloud coverage index to judge 1.3m height point cloud effectiveness, and using weighted least square circle fitting to directly extract diameter at breast height for trees with effective point cloud, and using pre-trained multi-feature fusion machine learning model to estimate diameter at breast height for trees with missing point cloud.

2. The method according to claim 1, wherein, Using unmanned aerial vehicle to collect point cloud data in forest, and the specific mode is as follows: An unmanned aerial vehicle equipped with a laser radar scanner is used to fly at low altitude in the forest of the target stand, and an AI vision obstacle avoidance system is used to achieve safe autonomous flight in a complex forest environment; the laser radar performs 360-degree scanning, and a multi-angle, multi-route flight strategy is used to collect point cloud data in all directions of the tree trunk and obtain a forest point cloud data set , wherein the ith point contains three-dimensional coordinates (x i , y i , z i ) and reflection intensity information I i .

3. The method according to claim 1, wherein, Using unmanned aerial vehicle to collect point cloud data on forest, and the specific mode is as follows: An unmanned aerial vehicle equipped with a laser radar scanner flies above the target forest canopy, adopts regular routes for round-trip flight scanning, the route spacing is calculated according to the field angle of the laser radar and the flight height, the complete coverage of the forest canopy point cloud is ensured, and the forest point cloud data set is obtained wherein, represents the jth point data of the forest point cloud data set.

4. The method according to claim 1, wherein, Multi-source point cloud registration and fusion, and the specific mode is as follows: (1) respectively pre-processing forest point cloud data set Pu and forest point cloud data set Po, including denoising, filtering and coordinate system unification; (2) extracting same name feature points in two groups of point cloud data sets as registration control points, and the feature points include tree top feature points, stem bifurcation points or pre-laid artificial target points; (3) The ICP (Iterative Closest Point) algorithm is used to solve the optimal rigid transformation, which minimizes the root mean square error of the registration. The point cloud dataset on the forest is transformed into , and the registration error is defined as: wherein, and are corresponding points in the in-forest point cloud dataset and the transformed on-forest point cloud dataset, respectively, R is a rotation matrix, T is a translation vector, and n is the number of corresponding points; the iterative optimization is performed until the RMSE converges to below a set threshold; (4) adaptively weighting and fusing the registered point cloud, in the overlapping area, calculating the fusion weight according to the point cloud quality factor, and the calculation formula of the point cloud quality factor Q is as follows: where p is the local point cloud density, I is the average reflection intensity, s z is the standard deviation in the height direction, the corresponding parameters with the subscript max represent the maximum value of the corresponding parameters, and a, b, and g are weight coefficients and satisfy a + b + g = 1. Calculating the fusion weight based on the point cloud quality factor: wherein, is a fusion weight for the inter-canopy point cloud, is a fusion weight for the on-canopy point cloud, is a quality factor for the inter-canopy point cloud, is a quality factor for the on-canopy point cloud; Fusion point cloud density The calculation formula is: wherein, is the density of the fused point cloud, , are the original densities of the in-forest and over-forest point clouds at the respective positions.

5. The method according to claim 1, wherein, Generating ground DEM and performing point cloud height normalization, and the specific mode is as follows: (1) using cloth simulation filtering algorithm CSF to extract ground point cloud from the point cloud after multi-source point cloud registration and fusion, simulating the fitting process of cloth and terrain from the top, and realizing the separation of ground point cloud and non-ground point cloud; (2) using ground point cloud, generating continuous ground digital elevation model DEM by using Kriging interpolation method; (3) Height normalization is performed on all point clouds, and the normalized height h i The calculation is: where z i is the original elevation of the point, DEM(xi, yi) is the DEM elevation value corresponding to the respective point position.

6. The method according to claim 1, wherein, Using improved marker-controlled watershed algorithm to perform single tree segmentation, and the specific mode is as follows: (1) generating canopy height model CHM based on normalized point cloud, and CHM is the maximum normalized height of point cloud in the corresponding grid; (2) performing Gaussian smoothing on CHM, and using local maximum filtering method to detect tree crown vertex; setting window size as W, if the pixel value is greater than all pixel values in the W*W neighborhood, it is marked as a tree crown vertex candidate point; the window size is adaptively determined according to the average crown width of the stand: wherein is an estimate of the stand average crown width, k is a proportionality factor; (3) taking the detected tree crown vertex as seed point, using marker-controlled watershed algorithm to perform single tree segmentation, and obtaining independent point cloud subset of each tree.

7. The method according to claim 1, wherein, Constructing multi-dimensional feature vector for each tree, calculating crown geometric feature, crown structure feature, visible part of stem feature, point cloud distribution feature and neighborhood competition feature, and the specific mode is as follows: For each segmented tree, 23 feature parameters in five categories are extracted from its point cloud data to construct a feature vector ; (1) crown geometric feature: Crown width CW: CW EW is the maximum width in the east-west direction, CW NS is the maximum width in the north-south direction; Crown projection area CA: calculated by using concave hull algorithm, and the concave hull parameter a is adaptively determined according to point cloud density: wherein, is the average distance of the point cloud, k is an empirical coefficient; Crown volume CV: where N v is the number of voxels comprising the point cloud, s is the voxel edge length; Crown surface area CSA: calculating three-dimensional convex hull surface area by using Alpha-Shape algorithm; Crown length CL: Among them, H max H is the height of the highest point of the tree crown. cb The height below the branch is defined as the first time the accumulated point cloud density from the base of the tree upwards exceeds the density threshold ρ. th Height; Crown shape index CSI: (2) crown structure feature: Crown density CD: where N c is the total number of crown point clouds, and CV is the crown volume. Crown point cloud vertical distribution ratio VDR: divide the crown into m layers, and the point cloud ratio VDR of the jth layer is j is: wherein N j is the number of point clouds for the jth layer; Leaf area index LAI: Based on point cloud transmittance, using modified Beer-Lambert model: where θ is the laser incident zenith angle, k ext is the extinction coefficient, P gap is the gap ratio, and the calculation formula of the gap ratio is: where N ground is the number of laser pulses that penetrated the canopy to reach the ground, N total is the total number of laser pulses that entered the canopy; Crown porosity GP: where N gap is the number of mesh points without point cloud coverage, N grid is the total number of mesh points; Crown complexity index CCI: Defined as the ratio of crown surface area to the surface area of an equivalent volume sphere: (3) Visible part of stem features: Extract the cross-section point cloud at height h in the trunk point cloud, fit a circle with weighted least squares method, diameter D h The calculation is obtained by minimizing the objective function: wherein (x i , y i ) are the cross-section point cloud coordinates, (x c , y c ) are the circle center coordinates, r is the radius, w i is the point weight, and D h = 2r. Point weight w i Adaptive calculation from point to fitted circle distance: where d i is the distance of point i to the current iteration circle, σ d is the distance standard deviation; Visible lowest diameter D L Search in the trunk point cloud for the lowest effective height h s above the height of the shrub layer h L Diameter at the lowest effective height h where N h is the number of point clouds at height h, N th is the effective point number threshold; The diameter D at the highest point H : Extract the diameter at the trunk-crown interface; Visible length of stem VL: wherein h H is the height of the highest visible point, h L is the height of the lowest visible point; Taper of stem T: (4) Point cloud distribution features: Point cloud height quantile H q : Calculate the qth quantile of the point cloud height of the whole tree, where q takes values 25, 50, 75, 90, 95: where H (k) represents the kth value in the height sequence, and n is the total number of point clouds; Height coefficient of variation of point cloud HCV: wherein σ H is the point cloud height standard deviation, is the point cloud height mean; Point cloud height skewness H skew : Point cloud height kurtosis H kurt : Point cloud intensity statistical features: intensity mean : Intensity standard deviation σ I : Intensity coefficient of variation ICV: (5) Neighborhood competition features: Number of neighboring trees NN: count number of trees within a radius R centered at the target tree location n Number of neighboring trees NN: count number of trees within a radius R centered at the target tree location where N is the total number of trees in the stand except the target tree, j is the tree traversal index, dist j is the horizontal distance between the target tree and the jth tree, 1(·) is the indicator function, R n = 5m; Average neighboring distance AND: Calculate the average distance between the target tree and the nearest k neighboring trees: where dist (j) is the jth closest distance to a neighboring tree after sorting by distance; Hegyi competition index CI: where D i is the target tree diameter or crown width, D j is the diameter or crown width of the jth neighboring tree, dist ij is the horizontal distance between the two trees; Stand density index SDI: Using Reineke stand density index formula: wherein N is the number of plants per unit area, D is the average diameter of the stand, ref D is the reference diameter.

8. The method according to claim 1, wherein, Based on the fusion point cloud, the height of each tree is extracted, and the 99th percentile height is used as the robust estimate of tree height to exclude the influence of outliers: Introduce effective point cloud coverage index to judge the effectiveness of point cloud at 1.3m height, for point cloud effective trees, use weighted least squares circle fitting to directly extract the diameter at breast height, for point cloud missing trees, use pre-trained multi-feature fusion machine learning model to estimate the diameter at breast height, the specific way is: wherein (x0, y0) is the tree base position coordinate, H is the ground elevation at the tree base position, (0 .99n) is the 99th percentile value of the point cloud height.

9. The method according to claim 1, wherein, (1) Calculate the effective point cloud coverage index EPCI: Set a sampling layer with thickness ∆h = 0.1m at 1.3m height, extract the point cloud in the corresponding layer; Divide 360 degrees into 12 sectors, each sector 30 degrees, calculate the angle coverage rate ACR and the spatial distribution uniformity SDU: The calculation formula of effective point cloud coverage index EPCI is: where N total sector = 12, N sector is the number of sectors containing the point cloud; where σ sector is the standard deviation of the number of point clouds in each sector, is the mean of the number of point clouds in each sector, and ε is a small constant to prevent division by zero. (2) Diameter at breast height measurement judgment and execution: wherein N 1.3 is the number of point clouds in the 1.3m height sampling layer, N th is the number of point clouds threshold, and ω1, ω2, ω3 are weight coefficients and ω1 + ω2 + ω3 = 1; Randomly select 3 points from the 1.3m height sampling layer point cloud to fit the initial circle; When EPCI ≥ EPCI th EPCI is an effective point cloud coverage index, EPCI th is an effective point cloud coverage index threshold, a robust circle fitting method based on RANSAC is used to directly extract the breast diameter: Calculate the distance of all points to the initial circle, count the number of inliers, the inlier determination condition is: Iterate 100 times, select the circle with the most inliers as the optimal fitting result; where d i is the distance from point i to the center of the circle, r is the radius of the circle, and e is the distance threshold. Use weighted least squares method to fine fit the inliers to get the final diameter at breast height: Calculate the fitting confidence FC: Multi-feature fusion machine learning model, as follows: where N inlier is the number of inliers, RMSE fit is the root mean square error of the fit. When EPCI < EPCI th a pre-trained multi-feature fusion machine learning model is used to estimate the chest diameter.

10. The method according to claim 9, wherein, (2) Collinearity diagnosis: Calculate the variance inflation factor VIF between features: (1) Feature screening: the mutual information-based feature selection method is adopted to calculate the mutual information MI between each feature F i and the breast diameter DBH: wherein is the ith feature variable, f is the value of the feature is the value of the diameter at breast height DBH, d is the value of the diameter at breast height DBH, the joint probability distribution, the marginal probability of the feature, the marginal probability of the diameter at breast height, the features with a screening mutual information value greater than a threshold value MI th . (3) Model training: Use gradient boosting decision tree algorithm GBDT, loss function is Huber loss: wherein, R2is the regression determination coefficient with the ith feature as the dependent variable and the remaining features as the independent variables, and redundant features with VIF > 10 are removed. (4) Model evaluation index: where y is the measured diameter at breast height, for predicting the diameter at breast height, δ is a threshold parameter; the parameters are optimized by grid search and 5-fold cross-validation; Root mean square error RMSE: Determination coefficient R 2 : Mean absolute percentage error MAPE: ​ The final model needs to satisfy R 2 > 0.85 and MAPE < 10%.

Images