A method for inverting leaf area index of mountain forests based on ground-based laser point cloud data

Abstract

The present invention discloses a mountain forest leaf area index inversion method based on ground-based laser point cloud data, comprising the following steps: Step 1, using a ground-based laser radar to set up multiple stations to collect dense laser point cloud data in a mountain forest area; Step 2, point cloud data registration, cropping, ground point filtering, and height normalization preprocessing; Step 3, point cloud voxelization, and calculation of the effective leaf area index (eLAI) by statistical collision probability; Step 4, correction for clustering effects, estimation of the canopy clustering index by the Pielou separation index method (PCS), and correction of the effective leaf area index (eLAI) to obtain the leaf area index (LAI); Step 5, verification in different forest scenarios. The present invention simultaneously considers the impact of terrain undulation and clustering effects on mountain forest leaf area index estimation, and can maintain high accuracy in scenarios with different vegetation structures and terrain conditions. Compared with traditional inversion methods, it has better versatility and higher inversion accuracy.

Description

Technical Field

[0001] The present invention relates to the field of forest leaf area index inversion, and in particular to a mountain forest leaf area index inversion method based on ground-based laser point cloud data. Background Art

[0002] The Leaf Area Index (LAI) is a key parameter describing forest canopy structure and function, directly reflecting forest productivity and health. Accurately measuring and retrieving LAI is crucial for understanding and predicting changes in forest ecosystems and formulating sound forest management policies (Tang et al., 2016). Traditional LAI measurement methods rely primarily on field measurements. While these methods offer high accuracy, they face numerous challenges in practical application. In recent years, with the rapid development of remote sensing technology, LAI retrieval methods based on remotely sensed data have become a hot topic of research. Remote sensing technology can rapidly acquire surface information over a wide area, unrestricted by topography, effectively addressing the shortcomings of traditional measurement methods (Zheng & Moskal, 2009). Terrestrial Laser Scanning (TLS), with its ground-based data acquisition capabilities and ability to obtain detailed three-dimensional structural information, offers unique advantages in forest ecosystem monitoring and provides new possibilities for LAI retrieval.

[0003] Methods used to invert LAI using ground-based lidar mainly include the gap ratio-based inversion method and the collision probability-based voxel method (Wang & Fang, 2020). Many researchers at home and abroad have used ground-based lidar data to invert LAI. For example, Takeda et al. divided the forest horizontally and vertically into voxels, calculated the gap ratio and leaf area for each voxel, and determined the three-dimensional distribution of canopy leaf area. The resulting vertical distribution of leaf area correlated well with the true value, but slightly overestimated the leaf area density at the top of the canopy (Takeda et al., 2008). Chen et al. used the true leaf area index obtained from stratified sampling to verify the leaf area index inverted at different layer thicknesses and voxel sizes. They determined the optimal combination of layer thickness and voxel size for accurate LAI estimation, but the general applicability is limited (Chen et al., 2022).

[0004] Some current studies assume a simple random leaf distribution, and the resulting LAI values ​​are actually effective LAI. However, vegetation clumping is common in forests (Fang, 2021; García et al., 2015), including clumping within leaf clusters, clumping of leaves within tree crowns, and clumping between tree crowns. When two vegetation canopies have the same LAI, the clumping canopy will contain more gaps than the randomly distributed canopy. Therefore, accounting for clumping effects is crucial when accurately estimating LAI using indirect optical methods. Failure to consider the clumping index (CI) can result in LAI errors of up to 70% (Fang, 2021). Previous methods for calculating the clumping index from ground-based lidar data include the finite length averaging method (LX), the gap size distribution method (CC), a correction method combining the finite length averaging method and the gap size distribution method (CLX), and the path length distribution method (Kuusk et al., 2018). For example, Zhu et al. estimated the leaf angle distribution of deciduous and coniferous forests to obtain the extinction coefficient, calculated the aggregation index using the gap size distribution method, and further improved the classification method of wood components to obtain the wood proportion. Using these parameters, they corrected the LAI and VFP (Vertical Foliage Profile) to correct for the underestimation caused by the aggregation effect and the height caused by the wood component, effectively improving the estimation accuracy of LAI and VFP (Zhu et al., 2018).

[0005] Therefore, existing research methods have made significant progress in LAI inversion. However, due to the undulating terrain and complex and changeable canopy structure of mountain forests, the accuracy and versatility of LAI inversion are still lacking. Summary of the Invention

[0006] Purpose of the invention: The present invention provides a mountain forest leaf area index inversion method based on ground-based laser point cloud data, which solves the problem that traditional LAI inversion methods have low accuracy and versatility and do not meet actual needs.

[0007] Technical solution: The present invention provides a method for inverting the leaf area index of mountain forests based on ground-based laser point cloud data, comprising the following steps:

[0008] Step 1: Use ground-based lidar to set up multiple stations to collect dense laser point cloud data in mountainous forest areas;

[0009] Step 2: Preprocess the laser point cloud data;

[0010] Step 3: voxelize the point cloud and calculate the effective leaf area index (eLAI) by statistically analyzing the collision probability;

[0011] Step 4: Correct the aggregation effect, estimate the canopy aggregation index using the Pielou separation index method (PCS), and correct the eLAI to obtain the LAI;

[0012] Step 5: Verify in different forest scenarios.

[0013] Furthermore, in step 1, the specific steps of setting up multiple stations to collect dense laser point cloud data in mountain forest areas using ground-based lidar are as follows: determining the sample area that needs TLS scanning, setting up 9 TLS instruments around the area to minimize occlusion and ensure the integrity of the point cloud data.

[0014] Furthermore, in step 2, the laser point cloud data is preprocessed by registration, cropping, ground point filtering, and height normalization.

[0015] Furthermore, in step 2, point cloud data registration, cropping, ground point filtering, and height normalization preprocessing specifically include the following steps:

[0016] Step 21: Accurately register the TLS point cloud data of the nine stations, align them to the same coordinate system, and crop the point cloud of the sample area;

[0017] Step 22: The registered and cropped point cloud data is filtered using the CSF cloth simulation filtering method to perform ground point filtering, and the point cloud is classified into vegetation points and ground points.

[0018] Step 23: Perform height normalization based on the filtered ground points and extract the vegetation point cloud of the sample area.

[0019] Furthermore, in step 3, the point cloud is voxelized and the eLAI is calculated by statistically analyzing the collision probability, which specifically includes the following steps:

[0020] Step 31: voxelize vegetation point cloud;

[0021] Step 32: Calculate the probability of laser interception based on the voxelized point cloud model;

[0022] Step 33: Based on the probability of laser beam being intercepted, the eLAI calculation formula of the k-th horizontal layer is:

[0023] l(k)=α(θ)N(k)

[0024] Where l(k) represents the leaf area density of the kth horizontal layer, α(θ) represents the correction factor of the laser beam direction and the leaf inclination angle. Assuming that the canopy leaves in the plot are randomly distributed, the correction factor of the leaf inclination angle is taken as 1.1;

[0025] Step 34: Accumulate the eLAI of each horizontal layer to obtain the total eLAI of the plot. The calculation expression is as follows:

[0026]

[0027] Furthermore, in step 31, the vegetation point cloud voxelization is specifically as follows: determining the range of the sample point cloud, taking the minimum value of the three-dimensional coordinates XYZ of the bounding box of the sample point cloud (X min , Y min , Z min ) is used as the starting point of voxel segmentation, and the point cloud data is divided into i×j×k grids with a voxel size of Δi×Δj×Δk. The voxelization calculation formula is:

[0028]

[0029] Among them, (i, j, j) is the number of voxels in the X, Y, and Z axis directions of the point cloud data, int is the rounding function, (x max ,y max ,z max ) is the maximum value of the sample point cloud bounding box, (x min ,y min ,z min ) is the minimum value of the bounding box of the sample point cloud, and (Δi, Δj, Δk) is the voxel size in the X, Y, and Z axis directions.

[0030] Furthermore, in step 32, the probability of laser interception is calculated based on the voxelized point cloud model as follows: if there is a point cloud distributed within the voxel, it means that the laser beam is intercepted and marked as a point cloud voxel; otherwise, it is a canopy gap voxel, indicating that the laser beam is not intercepted. After the sample point cloud is voxelized, it is divided into k layers in the Z-axis direction. By statistically analyzing the laser collision probability of each horizontal layer after voxel stratification, the ratio of voxels with point clouds to those without point clouds is calculated to obtain the gap rate of the layer:

[0031]

[0032] Among them, N k represents the collision probability of the laser beam at the kth level, n I (k) represents the number of laser interceptions at the kth level, i.e., the number of voxels with point clouds, n T (k) represents the total number of voxels in the k-th level layer.

[0033] Furthermore, in step 4, the aggregation effect is corrected, and the canopy aggregation index is estimated by the Pielou separation index method PCS. The LAI is obtained by correcting the eLAI. Specifically, the aggregation index is calculated by the Pielou separation index method (PCS). This method is a coefficient used to estimate the degree of separation between two co-growing plants. This method provides a randomness measure for each species relative to other species, where a and b are the probabilities of encountering species A and species B in the two species populations, respectively. Assuming a + b = 1, for two randomly distributed populations, there is a 95% probability:

[0034]

[0035] in, and is the maximum likelihood estimate of A and B, μ a and μ b is the average length of A and B, and is the variance of a and b; if PCS < 1, the species are clustered, if PCS = 1, the species are randomly distributed, and if PCS > 1, the species are uniformly distributed;

[0036] Correcting eLAI to obtain plot LAI, the correction expression is:

[0037] LAI=eLAI / λ PCS .

[0038] Furthermore, in step 5, verification was carried out in different forest scenarios. A three-dimensional radiation transfer model and computer simulation technology were used to construct forest scenarios with different vegetation types (including broad-leaved forests and coniferous forests) and terrain conditions (including flat land, sloping land, and irregular undulating terrain). TLS point cloud data were simulated and the LAI of the sample plots was estimated according to the above method. The data was then compared with the LAI provided by the scenario to verify its estimation accuracy.

[0039] Beneficial effects: Compared with the existing technology, the present invention has the following significant advantages: In the estimation of mountain forest leaf area index (LAI), the present invention fully considers the dual effects of terrain undulation and vegetation aggregation effect, thereby significantly improving the accuracy and stability of LAI inversion under complex vegetation structure and changeable terrain conditions; compared with traditional LAI inversion methods, the present invention not only has stronger versatility and adaptability, but also can significantly improve the accuracy of estimation, reduce deviation, and provide more reliable data support for ecological monitoring and forest resource management. BRIEF DESCRIPTION OF THE DRAWINGS

[0040] Figure 1 Schematic diagram of the method of the present invention.

[0041] Figure 2 This is a scatter plot of the LAI accuracy results of the LAI estimation method of the present invention in all scenarios. DETAILED DESCRIPTION

[0042] like Figure 1 As shown in FIG, a mountain forest leaf area index inversion method based on ground-based laser point cloud data includes the following steps:

[0043] Step 1: Use ground-based lidar to set up multiple stations to collect dense laser point cloud data in mountainous forest areas;

[0044] Step 2: Preprocessing of point cloud data including registration, cropping, ground point filtering, height normalization, etc.

[0045] Step 3: voxelize the point cloud and calculate the eLAI by statistically analyzing the collision probability;

[0046] Step 4: Correct the aggregation effect, estimate the canopy aggregation index using the Pielou separation index method (PCS), and correct the eLAI to obtain the LAI;

[0047] Step 5: Verify this method in different forest scenarios.

[0048] Preferably, the specific steps of collecting point cloud data of mountain forest areas by ground-based lidar in step 1 are: determining the sample area that needs TLS scanning, setting up 9 TLS instruments around the area to minimize occlusion and ensure the integrity of the point cloud data.

[0049] To eliminate interference from factors other than terrain and vegetation structure, this paper validates the aforementioned calculation method using a three-dimensional radiative transfer model combined with computer simulation technology. First, a virtual forest scene library containing 480 different forest scenes was constructed, with the scene information shown in Table 1. TLS point cloud data was generated through simulation, and the accuracy of the aforementioned TLS inversion method was verified using scene ground truth. The results demonstrate the accuracy and versatility of this method for LAI inversion in different scenarios.

[0050] Table 1 Virtual forest scene information

[0051]

[0052] Preferably, the specific steps of pre-processing such as point cloud data registration, cropping, ground point filtering, and height normalization in step 2 are as follows:

[0053] Step 2.1: Accurately register the TLS point cloud data of the nine stations, align them to the same coordinate system, and crop the point cloud of the sample area;

[0054] Step 2.2: Perform CSF ground point filtering on the registered and cropped point cloud data, and classify the point cloud into vegetation points and ground points;

[0055] Step 2.3: Perform height normalization based on the filtered ground points and extract the vegetation point cloud of the sample area.

[0056] Preferably, the specific steps of voxelizing the point cloud in step 3 and calculating the eLAI by statistically analyzing the collision probability are as follows:

[0057] Step 3.1: Voxelize vegetation point cloud. Determine the size of the sample point cloud, and use the minimum value of the three-dimensional coordinates XYZ of the bounding box of the sample point cloud (X min , Y min , Z min ) is used as the starting point of voxel segmentation, and the point cloud data is divided into i×j×k grids with a voxel size of Δi×Δj×Δk. The voxelization calculation formula is:

[0058]

[0059] Among them, (i, j, k) is the number of voxels in the X, Y, and Z axis directions of the point cloud data, int is the rounding function, (x max ,y max ,z max ) is the maximum value of the sample point cloud bounding box, (x min ,y min ,z min ) is the minimum value of the bounding box of the sample point cloud, and (Δi, Δj, Δk) is the voxel size in the X, Y, and Z axis directions;

[0060] Step 3.2: Calculate the probability of laser interception based on the voxelized point cloud model. If there is a point cloud distributed within the voxel, it means that the laser beam is intercepted and is marked as a point cloud voxel. Otherwise, it is a canopy gap voxel, indicating that the laser beam is not intercepted. After the sample point cloud is voxelized, it is divided into k layers in the Z-axis direction. By statistically analyzing the laser collision probability of each horizontal layer after voxel stratification, the ratio of voxels with point clouds to those without point clouds is calculated to obtain the gap rate of this layer:

[0061]

[0062] Among them, N k represents the collision probability of the laser beam at the kth level, n I (k) represents the number of laser interceptions at the kth level, i.e., the number of voxels with point clouds, n T (k) represents the total number of voxels in the k-th level;

[0063] Step 3.3: Based on the probability of laser beam being intercepted, the eLAI calculation formula for the k-th horizontal layer is:

[0064] l(k)=α(θ)N(k)

[0065] Where l(k) represents the leaf area density of the kth horizontal layer, α(θ) represents the correction factor of the laser beam direction and the leaf inclination angle. Assuming that the canopy leaves in the plot are randomly distributed, the correction factor of the leaf inclination angle is taken as 1.1;

[0066] Step 3.4: Accumulate the eLAI of each horizontal layer to obtain the total eLAI of the plot. The calculation expression is as follows:

[0067]

[0068] Preferably, in step 4, the aggregation effect is corrected, the canopy aggregation index is estimated by PCS, and the eLAI is corrected to obtain LAI. The specific steps are as follows:

[0069] The clustering index was calculated using the PCS method, which is a coefficient used to estimate the degree of separation between two co-growing plants. This method provides a measure of randomness for each species relative to the other species, where a and b are the probabilities of encountering species A and species B in the two populations, respectively, assuming a + b = 1. For two randomly distributed populations, there is a 95% probability that:

[0070]

[0071] in, and is the maximum likelihood estimate of A and B, μ a and μ B is the average length of A and B, and is the variance of a and b. If PCS < 1, the species are clustered, if PCS = 1, the species are randomly distributed, and if PCS > 1, the species are uniformly distributed;

[0072] Correcting eLAI to obtain plot LAI, the correction expression is:

[0073] LAI=eLAI / λ PCS

[0074] Preferably, the specific steps for verifying this method in different forest scenarios in step 5 are as follows: using a three-dimensional radiation transfer model and computer simulation technology to construct forest scenarios with different vegetation types (including broad-leaved forests and coniferous forests) and terrain conditions (including flat land, sloping land, and irregular undulating terrain), simulating to obtain TLS point cloud data, estimating the LAI of the sample plot according to the above method, and comparing it with the LAI provided by the scenario to verify its estimation accuracy.

[0075] To evaluate the above calculation method, we chose R 2(coefficient of determination), RMSE (root mean square error) and nRMSE (normalized root mean square error) were used as evaluation indicators. 2 The mathematical formulas for , RMSE, and nRMSE are as follows:

[0076]

[0077] Among them, i represents the sequence number of the i-th sample point, Represents the predicted value of the i-th sample point, y i represents the true value of the i-th sample point, Represents the average value of the true value of all sample points, m represents the number of sample points, and mean(y) represents the average value of the true value of the sample points.

[0078] The PCS clustering index estimation method is used to estimate the clustering index of the sample plot and perform eLAI correction to obtain LAI. The accuracy of the estimated LAI result is evaluated using the reference true value of the virtual sample plot. Considering the impact of terrain undulation on the accuracy of TLS inversion LAI, the present invention compares the estimation accuracy of the PCS method under different vegetation types and different slope conditions. The results are shown in Table 2. It can be seen that the PCS index correction for clustering effect can obtain better estimation accuracy (R 2 =0.9061(broad-leaved forest), R 2 =0.8907 (coniferous forest). Although the LAI estimation method of the present invention decreases with the increase of terrain undulation, it can still maintain a good accuracy (R 2 >0.7350). The LAI estimation results under all scenarios are summarized as follows: Figure 2 As shown (R 2 =0.9295), which can meet the actual needs.

[0079] Table 2 Summary of LAI estimation accuracy under different vegetation types and terrain conditions

[0080]

[0081] In the table, Slope0-Slope50 indicates regular terrain with a slope of 0°-50°, Slopeterrain1 and Slopeterrain2 indicate irregular terrain, and All indicates the summary of all terrain scenarios.

Claims

1. A mountain forest leaf area index inversion method based on ground-based laser point cloud data, characterized in that: The steps include: Step 1: Use ground-based lidar to set up multiple stations to collect dense laser point cloud data in mountainous forest areas; Step 2: Preprocess the laser point cloud data; Step 3: Voxelize the point cloud and calculate the eLAI by statistically analyzing the collision probability. This includes the following steps: Step 31: voxelize vegetation point cloud; specifically, determine the size of the sample point cloud, and use the minimum value of the three-dimensional coordinates XYZ of the bounding box of the sample point cloud (X min , Y min , Z min ) is used as the starting point of voxel segmentation, and the point cloud data is divided into i×j×k grids with a voxel size of Δi×Δj×Δk. The voxelization calculation formula is: Among them, (i, j, k) is the number of voxels in the X, Y, and Z axis directions of the point cloud data, int is the rounding function, (x max ,y max ,z max ) is the maximum value of the sample point cloud bounding box, (x min ,y min ,z min ) is the minimum value of the bounding box of the sample point cloud, and (Δi, Δj, Δk) is the voxel size in the X, Y, and Z axis directions; Step 32: Calculate the probability of laser interception based on the voxelized point cloud model. Specifically, if there is a point cloud distributed within the voxel, it means that the laser beam is intercepted and is marked as a point cloud voxel. Otherwise, it is a canopy gap voxel, indicating that the laser beam is not intercepted. After the sample point cloud is voxelized, it is divided into k layers in the Z-axis direction. By statistically analyzing the laser collision probability of each horizontal layer after voxel stratification, the ratio of voxels with point clouds to those without point clouds is calculated to obtain the gap rate of the layer: Among them, N k represents the collision probability of the laser beam at the kth level, n I (k) represents the number of laser interceptions at the kth level, i.e., the number of voxels with point clouds, n T (k) represents the total number of voxels in the k-th level; Step 33: Based on the probability of laser beam being intercepted, the eLAI calculation formula of the k-th horizontal layer is: l(k)=α(θ)N(k) Where l(k) represents the leaf area density of the kth horizontal layer, α(θ) represents the correction factor of the laser beam direction and the leaf inclination angle. Assuming that the canopy leaves in the plot are randomly distributed, the correction factor of the leaf inclination angle is taken as 1.1; Step 34: Accumulate the eLAI of each horizontal layer to obtain the total eLAI of the plot. The calculation expression is as follows: Step 4: Correct the aggregation effect, estimate the canopy aggregation index using the Pielou separation index method (PCS), and correct the eLAI to obtain the LAI; Step 5: Verify in different forest scenarios.

2. The mountain forest leaf area index inversion method based on ground-based laser point cloud data according to claim 1, characterized in that: In step 1, the ground-based lidar is used to set up multiple stations to collect dense laser point cloud data in mountain forest areas. The specific steps are as follows: determine the sample area that needs TLS scanning and set up 9 TLS instruments around the area.

3. The mountain forest leaf area index inversion method based on ground-based laser point cloud data according to claim 1, characterized in that: In step 2, the laser point cloud data is preprocessed by registration, cropping, ground point filtering, and height normalization.

4. The mountain forest leaf area index inversion method based on ground-based laser point cloud data according to claim 1, characterized in that: In step 2, point cloud data registration, cropping, ground point filtering, and height normalization preprocessing specifically include the following steps: Step 21: Accurately register the TLS point cloud data of the nine stations, align them to the same coordinate system, and crop the point cloud of the sample area; Step 22: The registered and cropped point cloud data is filtered using the CSF cloth simulation filtering method to perform ground point filtering, and the point cloud is classified into vegetation points and ground points. Step 23: Perform height normalization based on the filtered ground points and extract the vegetation point cloud of the sample area.

5. The mountain forest leaf area index inversion method based on ground-based laser point cloud data according to claim 1, characterized in that: In step 4, the aggregation effect is corrected and the canopy aggregation index is estimated by the Pielou separation index method PCS. The LAI is obtained by correcting the eLAI. Specifically, the aggregation index is calculated by the Pielou separation index method PCS, which is a coefficient used to estimate the degree of separation between two co-growing plants. It provides a randomness measure for each species relative to other species, where a and b are the probabilities of encountering species A and species B in the two populations, respectively. Assuming a + b = 1, for two randomly distributed populations, there is a 95% probability: in, and is the maximum likelihood estimate of A and B, μ a and μ b is the average length of A and B, and is the variance of a and b; if PCS < 1, the species are clustered, if PCS = 1, the species are randomly distributed, and if PCS > 1, the species are uniformly distributed; Correct eLAI to get plot LAI. The corrected expression is: LAI=eLAI / λ PCS 。 6. The mountain forest leaf area index inversion method based on ground-based laser point cloud data according to claim 1, characterized in that: In step 5, verification was performed in different forest scenarios. A three-dimensional radiation transfer model and computer simulation technology were used to construct forest scenarios with different vegetation types and terrain conditions. TLS point cloud data were simulated and the LAI of the sample plots was estimated according to the above method. The LAI was compared with the LAI provided by the scenario to verify its estimation accuracy.

7. The mountain forest leaf area index inversion method based on ground-based laser point cloud data according to claim 6, characterized in that: Different vegetation types include broad-leaved forests and coniferous forests, and terrain conditions include flat land, sloping land, and irregular undulating terrain.

Images