A fine site index precision evaluation method for Chinese fir plantations based on UAV-LiDAR

Abstract

The application discloses a kind of based on UAV-LiDAR's Chinese fir plantation site index precision evaluation method, belongs to the evaluation method field of forest site condition.The present application is based on the UAV-LiDAR point cloud data of target stand, carries out digital elevation model, digital surface model and canopy height model construction, adopts marker watershed control method and local maximum value method to the single tree segmentation and height extraction of the constructed canopy height model.AnselinLocal Moran's I is used to analyze the spatial aggregation distribution degree of the whole forest single tree height, and the spatial heterogeneity grade of the stand is divided, and the grid size is set according to the spatial heterogeneity grade, and the microtopography is referred to for stand gridding, and the site index of grid scale is accurately evaluated by comparing the site index table of the region Chinese fir plantation, and the fine site index spatial distribution map is drawn.The present application uses airborne laser radar technology and computer-aided means, and is a more efficient, meticulous and accurate Chinese fir plantation site index division and evaluation method.

Description

Technical Field

[0001] This invention belongs to the field of forest site condition evaluation methods, specifically involving a precise evaluation method for the fine position index of Chinese fir plantations based on UAV-LiDAR. Background Technology

[0002] Forest site quality not only affects the functioning of forests but also the sustainable development of the ecological environment and forest resources. Accurate assessment of forest site quality is a crucial basis for the precise cultivation of plantations and the foundation for achieving the principle of planting the right trees in the right locations. It provides vital information support for effectively improving the quality of plantation cultivation and forest productivity. Among commonly used methods for site quality assessment, the status index is more intuitive than the status grade, and the average height of dominant trees is less affected by silvicultural practices and density, making it relatively stable and the mainstream method for site assessment. The status index of a target stand is usually obtained by experienced forestry technicians after a field survey of the entire stand, visually estimating the dominant trees, calculating the average height of dominant trees, and finding the corresponding status index in a compiled status index table, combined with the stand's age. Due to the interactions between various natural factors, complex and changing topographic factors not only directly affect site quality but also indirectly affect it by influencing other factors. Therefore, for plain forests with relatively flat terrain, low spatial heterogeneity, and small differences in site quality, the status index obtained using traditional methods is relatively representative of the overall forest site quality. However, in complex mountainous terrain, the accuracy and representativeness of the status index obtained through traditional methods are limited. Furthermore, traditional methods rely heavily on the experience of surveyors, whose subjectivity can significantly impact the accuracy of the standard plot representation. Continuing to use traditional methods to survey status indices is insufficient to meet the needs of precision forestry development.

[0003] The aforementioned limitations mainly stem from the lack of an efficient and systematic site quality assessment technique that can accurately obtain fine-grained position indices for target forest stands. With technological advancements, high-efficiency and high-precision UAV remote sensing technology has been widely applied in forestry management and production. UAV-LiDAR can eliminate various interference factors in complex mountainous terrain, rapidly extract stand growth parameters at the individual tree scale, and perform high-precision positioning.

[0004] To address the limited accuracy of status index surveys in current forest site condition assessments, this invention fully considers the impact of spatial heterogeneity within forest stands on site quality. It quantifies stand spatial heterogeneity using the local Moran index and selectively chooses different grid sizes to extract fine-grained status indices at the grid scale. Combining airborne LiDAR data acquisition and computer-aided decision-making, this invention provides a more efficient, detailed, and accurate method for evaluating the fine-grained status index of Chinese fir plantations. Summary of the Invention

[0005] The purpose of this invention is to improve the accuracy of status index surveys in current forest site condition assessment work, and to provide a precise evaluation method for the refined status index of Chinese fir plantations based on UAV-LiDAR. This invention mainly focuses on the spatial heterogeneity of forest site quality, quantifying the levels, selecting appropriate dimensions, referencing forest stand micro-topography, performing stand gridding, refining the stand status index, and optimizing the accuracy of stand status index extraction.

[0006] The specific technical solution adopted in this invention is as follows:

[0007] This invention provides a precise evaluation method for the fine-grained status index of Chinese fir plantations based on UAV-LiDAR, as detailed below:

[0008] S1. Use UAV-LiDAR technology to collect LiDAR point cloud data of the target Chinese fir plantation, and preprocess the collected LiDAR point cloud data.

[0009] S2. Separate ground points from non-ground points in the preprocessed LiDAR point cloud data from S1; construct a model for the separated ground points to obtain a digital elevation model of the target forest stand; construct a digital surface model of the target forest stand based on the separated non-ground points; obtain a canopy height model based on the elevation difference between the digital elevation model and the digital surface model at the same coordinate point; obtain normalized LiDAR point cloud data based on the preprocessed LiDAR point cloud data from S1 and the digital elevation model.

[0010] S3. Analyze and create positive, negative, and reverse terrain data for the digital elevation model described in S2, and extract the ridges and valleys of the target forest and vectorize them.

[0011] S4. Perform single-tree segmentation on the canopy height model described in S2 based on the marked watershed control method and the local maximum method, and extract the single-tree height to obtain the spatial latitude and longitude projection coordinates and tree height of the target forest stand.

[0012] S5. Vectorize the spatial latitude and longitude projection coordinates of individual trees described in S4 to establish a database of individual tree spatial locations and tree heights in the target forest stand; use the local Moran index method to analyze the spatial clustering degree of individual tree heights in the target forest stand and conduct a graded evaluation of the spatial heterogeneity of the target forest stand.

[0013] S6. Based on the spatial heterogeneity rating results obtained in S5, and taking into account the spatial requirements of the target forest stand and the impact of micro-topography on the trees, set different grid sizes for different levels of fineness for the target forest stand.

[0014] S7. Using the vectorized valleys and ridges in S3 as micro-topographic boundary areas, overlay the refined grids of the corresponding size in S6. Referring to the normalized LiDAR point cloud data obtained in S2, perform gridding of the valleys and ridges within the target forest stand range, and vectorize the slope areas based on the valleys and ridges. Obtain and label the gridded valley areas, ridge areas, and slope areas of the target forest stand, and obtain the vectorized grid of the target forest stand.

[0015] S8. Based on the individual tree spatial location and tree height database described in S5 and the target forest stand vectorized grid described in S7, establish a database of individual tree height and number within each grid, and further statistically select the average height SH of the dominant trees within each grid. t ;

[0016] S9, based on the average height of the dominant trees in each grid in S8, SH t By referring to the status index table of Chinese fir plantations in the relevant province and city, the grid-scale fine status index of the target forest stand is obtained; then, combined with the remote sensing image of the target forest stand, the grid-scale fine status index is visualized and mapped.

[0017] Preferably, in step S1, the acquisition of LiDAR point cloud data is achieved by using a DJI Matrice 300RTK drone equipped with a lidar sensor to perform flight operations on the target forest stand.

[0018] Preferably, in step S1, the preprocessing method for LiDAR point cloud data includes 3D reconstruction and denoising, as detailed below:

[0019] The advanced version of DJI Terra software was used for 3D reconstruction of LiDAR point cloud data; CloudCompare software was used for manual denoising of non-target ground features; and LiDAR 360 software was used for denoising of drift noise points.

[0020] Preferably, in step S2, the separation of ground points and non-ground points in the LiDAR point cloud data is achieved using a cloth-based analog filtering algorithm, as detailed below:

[0021] (1) Invert the preprocessed LiDAR point cloud data in S1 and simulate a sufficiently soft cloth to be laid on top of the inverted LiDAR point cloud data. The smallest fiber molecules that make up the simulated cloth are regarded as countless cloth particles. The cloth particles are limited to moving only in the vertical direction. The position and velocity of the cloth particles are determined by the internal and external forces between the cloth particles. The ground points in the LiDAR point cloud data are classified by comparing the height values ​​of the cloth particles and the terrain.

[0022] (2) Set the grid size and mesh the simulated LiDAR point cloud in step (1) to perform point cloud filtering with the grid as the basic unit.

[0023] (3) Based on the results of step (2), all LiDAR point clouds and cloth particles in the target area are projected onto the same two-dimensional horizontal plane. The particles on the plane projected by the LiDAR point cloud are called point cloud particles. According to the two-dimensional coordinate distance between the point cloud particles and the cloth particles, the LiDAR point corresponding to each cloth particle in the two-dimensional plane is found. The height of the LiDAR point is defined as IH, and IH is the minimum height that the cloth particle can fall.

[0024] (4) Based on the results of step (3), calculate the displacement of the fabric particles under the action of external force. Since the particles are not subject to collision force in the air, but only to gravity, calculate the distance the fabric particles move under the action of gravity according to the following formula:

[0025]

[0026] In the formula, t represents a certain time node; Δt represents the time step; P(t) represents the position of the cloth particle at time node t under the influence of gravity alone; P(t+Δt) represents the position of the cloth particle after time Δt under the influence of gravity alone; P(t-Δt) represents the position of the cloth particle before time Δt at that node under the influence of gravity alone; m represents the mass of the cloth particle, and g is a constant.

[0027] (5) Based on the results of step (3), calculate the displacement distance of the fabric particles under the influence of internal forces, as shown in the following formula:

[0028]

[0029] In the formula, This represents the displacement vector of the cloth particle; b is a constant, with a value of 1 when the cloth particle can move and a value of 0 when it cannot move. The current coordinate vector representing the cloth particle; Represents the coordinate vector of adjacent particles connected to the cloth particle; A unit vector representing the vertical direction;

[0030] (6) Based on the results of steps (4) and (5), compare the height CH of the cloth particles that have been displaced by external and internal forces with IH; if CH≤IH, place them at the height of IH and set them as immovable particles; if CH>IH, treat them as movable particles and repeat steps (4) to (5) until the maximum change height of all cloth particles is small enough or immovable or exceeds the set maximum number of iterations. At this time, the plane composed of all cloth particles is an approximate real terrain, and the cloth simulation iteration is completed.

[0031] (7) Based on the results of step (6), calculate the distance MH between the cloth particles and the corresponding LiDAR points after the cloth simulation is completed;

[0032] (8) Based on the results of step (7), set the distance threshold OH; compare the magnitudes of MH and OH for each LiDAR point. When MH≤OH, the LiDAR point is marked as a ground point, and when MH>OH, the LiDAR point is marked as a non-ground point.

[0033] Preferably, in step S2, the digital elevation model and digital surface model are constructed using the Kriging interpolation method, as shown in the following formula:

[0034]

[0035] In the formula, n is the number of points that can determine the attribute value of that point; Z′0 is the estimated value of the point Z0 to be estimated, Z i W represents the attribute value of the i-th known point; i This represents the weight of the i-th known point;

[0036]

[0037] In the formula, γ ij Z represents the semivariance between the i-th and j-th neighborhood points; n The attribute value representing the nth known point;

[0038] Choose an appropriate function form and fit the function γ = f(d); the location of the point to be found is known, and d can be obtained. 10 , ..., d n0 Substituting into the function yields γ 10 , ..., γ n0 Calculate matrix multiplication:

[0039]

[0040] In the formula, λ is a constant; the 0th neighborhood point represents the point to be solved;

[0041] Finally, the valuation of the point Z0 to be valued is completed.

[0042] Preferably, S3 is as follows:

[0043] (1) Using ArcGIS software, based on the digital elevation model described in S2, statistical raster data Focal_St is obtained through the focus statistics method in terrain domain analysis; the difference between the digital elevation model and Focal_St at the same spatial location is calculated and raster data R1_Ca is obtained; R1_Ca is reclassified with 0 as the critical value, and R1_Ca above 0 is regarded as positive terrain raster Reclassify1, and R1_Ca below 0 is regarded as negative terrain raster Reclassify2;

[0044] (2) Based on the digital elevation model described in S2, perform raster calculations and obtain reverse terrain raster data Reclassify3;

[0045] (31) Based on the digital elevation model described in S2, perform depression filling, flow direction, flow rate analysis and raster calculation of the topographic hydrology to obtain the raster area data R2_Co with a flow rate of 0; determine the ridge classification threshold of R2_Co by analyzing the surface contour lines and mountain shadows of the digital elevation model data; reclassify R2_Co according to the classification threshold and obtain Reclassify4;

[0046] (32) Multiply Reclassify1 and Reclassify4 to obtain the raster data Reclassify; the raster with a reclassification retention value of 1 is the ridge raster area, and the ridge can be extracted by vectorizing the result;

[0047] (41) Based on Reclassify3, perform depression filling, flow direction, flow rate analysis and raster calculation for the topographic hydrology to obtain the raster area data R3_Co with a flow rate of 0; determine the valley classification threshold of R3_Co by analyzing the surface contour lines and mountain shadows of the digital elevation model data; reclassify R3_Co according to the classification threshold and obtain Reclassify5;

[0048] (42) Multiply Reclassify2 and Reclassify5 to obtain raster data Reclassify6; the raster with a reclassification retention value of 1 is the valley raster area, and the valley can be extracted by vectorizing the result;

[0049] (5) Multiply Reclassify by Reclassify6 to obtain raster data Reclassify7; the raster with a reclassification retention value of 0 is the hillside raster area, and the hillside can be extracted by vectorizing the result.

[0050] Preferably, S4 is as follows:

[0051] (1) Based on the height variation characteristics of a single tree crown, the local maximum algorithm is used to identify the treetop point of a single tree; the identification result is regarded as the X and Y coordinates of the tree, and the foreground is marked with it to optimize the watershed segmentation effect;

[0052] (2) The local maximum detection method based on the canopy height model is used to identify individual trees. During the iteration process, the window size is dynamically adjusted to determine the optimal value. The sliding window with the optimal window size is used to search for local maximum values. Finally, the coordinates of the identified individual trees are output and marked on the DOM.

[0053] (3) Invert the pixel values ​​of the canopy height model, mark the treetop points identified by the local maximum algorithm as the foreground, mark the non-canopy areas determined by threshold segmentation as the background, and mark the remaining uncertain parts as 0; then use the watershed algorithm, the labels will be updated each time water is poured, and when two different colored labels meet, a watershed dam will be built, eventually forming a closed outline; each local minimum and its influence range form a water basin, and the boundary of the water basin is the watershed;

[0054] (4) Take the projection and affine transformation of the input canopy height model, add projection coordinates to the image, and save it as a raster image in tif format; convert the raster data of the raster image into vector surface data according to the tree ID. There will be some small patches in the converted data. Merge these small patches into the nearest large patch.

[0055] Preferably, S5 is as follows:

[0056] Using ArcGIS software, local Moran's index analysis was performed on the obtained tree height values ​​of the target forest stand to obtain the clustering patterns of high and low tree height values. Based on the proportion of non-clustered tree points in non-high-high and low-low clustering patterns to the total number of tree points, the spatial heterogeneity of the target forest stand site quality was rated according to five levels.

[0057] Furthermore, in S6, five different grid sizes with varying degrees of fineness are customized to correspond to the five rating results in S5.

[0058] Preferably, in step S8, the individual trees in each grid are arranged by height, and the top 20% are selected as the dominant trees of the target tree species. The average height SH of the dominant trees in that grid is then calculated. t .

[0059] Compared with the prior art, the present invention has the following advantages:

[0060] 1) Comprehensive and scientific. Based on UAV-LiDAR data, the location and height parameters of all individual trees in the target forest stand are presented in a comprehensive and scientific manner, which makes up for the shortcomings of traditional manual survey methods, which can only obtain empirical values ​​based on the growth of the forest stand within the field of view.

[0061] 2) Quantitative Standards. Based on the local Moran index of single tree height in the whole stand, the spatial heterogeneity of site quality in the target stand is quantified into five levels. The grid size is selected based on this, which solves the problems of difficulty in quantifying spatial heterogeneity of site quality and difficulty in determining the degree of refinement of the site index.

[0062] 3) High efficiency and cost-effectiveness. The use of airborne LiDAR for data acquisition and computer-aided decision-making effectively shortened the time for field surveys and saved manpower and material costs. Attached Figure Description

[0063] Figure 1 A schematic diagram of the process of this invention;

[0064] Figure 2 3D reconstruction of UAV-LiDAR point cloud of the target forest stand;

[0065] Figure 3 DEM map (A); DSM map (B); CHM map (C); TNPC map (D) of the target forest stand;

[0066] Figure 4 Schematic diagram of the valleys and ridges of the target forest stand;

[0067] Figure 5 Individual tree segmentation results and spatial distribution map of the target forest stand;

[0068] Figure 6 Local Moran index plot of individual tree height in the target forest stand;

[0069] Figure 7 A grid diagram of the target forest stand;

[0070] Figure 8 Fine-scale position index thematic map of target forest stands based on UAV-LiDAR point clouds at the grid scale. Detailed Implementation

[0071] The present invention will be further described and illustrated below with reference to the accompanying drawings and specific embodiments. The technical features of each embodiment of the present invention can be combined accordingly, provided that there is no mutual conflict.

[0072] like Figure 1 As shown, this invention provides a precise evaluation method for the fine-grained status index of Chinese fir plantations based on UAV-LiDAR. The method is as follows:

[0073] S1. Use UAV-LiDAR technology (i.e., UAV equipped with a lidar gimbal) to collect LiDAR point cloud data of the target Chinese fir plantation, and perform preprocessing such as three-dimensional reconstruction and noise reduction on the collected LiDAR point cloud data.

[0074] In practical use, LiDAR point cloud data is collected by using a DJI Matrice 300RTK drone equipped with a lidar sensor (DJI Zenmuse L1 gimbal) to fly over the target forest stand.

[0075] In practical applications, the preprocessing methods for LiDAR point cloud data mainly include 3D reconstruction and denoising, as detailed below:

[0076] (1) 3D reconstruction of LiDAR point cloud data: DJI Terra software (advanced version) was used;

[0077] (2) Denoising of LiDAR point cloud after 3D reconstruction: CloudCompare software was used to manually denoise non-target ground objects (high voltage lines, iron towers, utility poles, etc.); LiDAR 360 software was used to denoise drift noise points.

[0078] S2. Separate ground points from non-ground points in the preprocessed LiDAR point cloud data from S1; construct a model for the separated ground points to obtain the digital elevation model (DEM) of the target forest stand; construct the digital surface model (DSM) of the target forest stand based on the separated non-ground points; obtain the canopy height model (CHM) based on the elevation difference between the DEM and the DSM at the same coordinate point; obtain normalized LiDAR point cloud data (TNPC) based on the preprocessed LiDAR point cloud data from S1 and the DEM.

[0079] In practical applications, the cloth-based analog filtering (CSF) algorithm is used to separate ground points from non-ground points in LiDAR point cloud data, as detailed below:

[0080] (1) Invert the preprocessed LiDAR point cloud data in S1 and simulate a sufficiently soft cloth to be laid on top of the inverted LiDAR point cloud data. The smallest fiber molecules that make up the simulated cloth are considered as countless particles, called cloth particles. The cloth particles are restricted to move only in the vertical direction. The position and velocity of the cloth particles are determined by the internal force (spring force) and external force (gravity and collision force of the particles themselves) between the cloth particles. The ground points in the LiDAR point cloud data are classified by comparing the height values ​​of the cloth particles and the terrain.

[0081] (2) Set the grid size and mesh the simulated LiDAR point cloud in step (1) to perform point cloud filtering with the grid as the basic unit.

[0082] (3) Based on the results of step (2), all LiDAR point clouds and cloth particles in the target area are projected onto the same two-dimensional horizontal plane. The particles on the plane projected by the LiDAR point cloud are called point cloud particles. According to the two-dimensional coordinate distance between the point cloud particles and the cloth particles, the LiDAR point corresponding to each cloth particle in the two-dimensional plane is found. The height of the LiDAR point is defined as IH, and IH is the minimum height that the cloth particle can fall.

[0083] (4) Based on the results of step (3), calculate the displacement of the fabric particles under the action of external force. Since the particles are not subject to collision force in the air, but only to gravity, calculate the distance the fabric particles move under the action of gravity according to the following formula:

[0084]

[0085] In the formula, t represents a certain time node; Δt represents the time step; P(t) represents the position of the cloth particle at time node t under the influence of gravity alone; P(t+Δt) represents the position of the cloth particle after time Δt under the influence of gravity alone; P(t-Δt) represents the position of the cloth particle before time Δt at that node under the influence of gravity alone; m represents the mass of the cloth particle (all cloth particles can be defined as having a constant mass of 1); g is a constant.

[0086] (5) Based on the results of step (3), calculate the displacement distance of the fabric particles under the influence of internal force (spring force between two fabric particles), as follows:

[0087]

[0088] In the formula, This represents the displacement vector of the cloth particle; b is a constant, with a value of 1 when the cloth particle can move and a value of 0 when it cannot move. The current coordinate vector representing the cloth particle; Represents the coordinate vector of adjacent particles connected to the cloth particle; The unit vector (0,0,1) represents the vertical direction;

[0089] (6) Based on the results of steps (4) and (5), compare the height CH of the cloth particles that have been displaced by external and internal forces with IH; if CH≤IH, place them at the height of IH and set them as “immovable” particles; if CH>IH, treat them as movable particles and repeat steps (4)~(5) until the maximum change height of all cloth particles is small enough or immovable or exceeds the set maximum number of iterations. At this time, the plane composed of all cloth particles is an approximate real terrain, and the cloth simulation iteration is completed.

[0090] (7) Based on the results of step (6), calculate the distance MH between the cloth particles and the corresponding LiDAR points after the cloth simulation is completed;

[0091] (8) Based on the results of step (7), set the distance threshold OH; compare the magnitudes of MH and OH for each LiDAR point. When MH≤OH, the LiDAR point is marked as a ground point, and when MH>OH, the LiDAR point is marked as a non-ground point.

[0092] In practical use, the construction of DEM, DSM, and CHM is as follows:

[0093] (1) The obtained ground points are used as the basic LiDAR point cloud set for constructing the DEM, and the non-ground points are used as the basic LiDAR point cloud set for constructing the DSM.

[0094] (2) The digital elevation model and digital surface model are constructed using the Kriging interpolation method, as shown in the following formula:

[0095]

[0096] In the formula, n is the number of points that can determine the attribute value of the point (generally, a certain number of known points around the point are taken); Z′0 is the estimated value of the point Z0 to be estimated; Z i W represents the attribute value of the i-th known point; i This represents the weight of the i-th known point;

[0097]

[0098] In the formula, γ ij Z represents the semivariance between the i-th and j-th neighborhood points; n The attribute value representing the nth known point;

[0099] Choose an appropriate function form and fit the function γ = f(d); the location of the point to be found is known, and d can be obtained. 10 , ..., d n0 Substituting into the function yields γ 10 , ..., γ n0 Calculate matrix multiplication:

[0100]

[0101] In the formula, λ is a constant; the 0th neighborhood point represents the point to be solved;

[0102] Finally, the valuation of the point Z0 to be valued is completed.

[0103] (3) Constructing CHM: CHM is obtained by calculating the difference between DSM and DEM at the same spatial location.

[0104] (4) Obtain normalized LiDAR point cloud data (TNPC): This is obtained by subtracting the elevation of the same latitude and longitude spatial location from the LiDAR point cloud data obtained in S1.

[0105] S3. Analyze and create positive, negative, and reverse terrain data from the digital elevation model obtained in S2, extract the ridges and valleys of the target forest and vectorize them.

[0106] In practical use, the steps are as follows:

[0107] (1) Analysis and creation of positive and negative terrain. Using ArcGIS software, based on the digital elevation model obtained by S2, statistical raster data Focal_St was obtained through the focus statistics method in terrain domain analysis; the difference between the digital elevation model and Focal_St at the same spatial location was calculated and raster data R1_Ca was obtained; R1_Ca was reclassified with 0 as the threshold value, and R1_Ca above 0 was regarded as positive terrain raster Reclassify1, and R1_Ca below 0 was regarded as negative terrain raster Reclassify2.

[0108] (2) Analysis and creation of reverse terrain. Based on the digital elevation model obtained from S2, raster calculations are performed to obtain reverse terrain raster data Reclassify3.

[0109] (3) Extraction of ridges

[0110] (31) Ridge generation. Based on the digital elevation model obtained in S2, the topographic hydrology of this terrain is analyzed for depression filling, flow direction, flow rate and raster calculation to obtain the raster area data R2_Co with a flow rate of 0; by analyzing the surface contour lines and mountain shadows of the digital elevation model data, the ridge classification threshold of R2_Co is determined; R2_Co is reclassified according to the classification threshold and Reclassify4 is obtained.

[0111] (32) Ridge extraction. Multiply Reclassify1 and Reclassify4 to obtain the raster data Reclassify; the raster with a reclassification retention value of 1 is the ridge raster area, and the ridge can be extracted by vectorizing the result.

[0112] (4) Extraction of valleys

[0113] (41) Valley formation. Based on Reclassify3, the topographic and hydrological data of this terrain are analyzed for depression filling, flow direction, flow rate and raster calculation to obtain the raster area data R3_Co with a flow rate of 0; by analyzing the surface contour lines and mountain shadows of the digital elevation model data, the valley classification threshold of R3_Co is determined; R3_Co is reclassified according to the classification threshold and Reclassify5 is obtained.

[0114] (42) Valley extraction. Multiply Reclassify2 and Reclassify5 to obtain raster data Reclassify6; the raster with a reclassification retention value of 1 is the valley raster area, and the valley can be extracted by vectorizing the result.

[0115] S4. The canopy height model obtained in S2 is segmented into individual trees based on the marked watershed control method and the local maximum method, and the individual tree height is extracted to obtain the spatial latitude and longitude projection coordinates and tree height of the target forest stand.

[0116] In practical use, the steps are as follows:

[0117] (1) Based on the height variation characteristics of a single tree crown, the local maximum (LM) algorithm is used to identify the treetop point of a single tree; the identification result is regarded as the X and Y coordinates of the tree, and the foreground is marked with it to optimize the watershed segmentation effect.

[0118] (2) The local maximum detection method based on the canopy height model is used to identify individual trees. During the iteration process, the window size is dynamically adjusted to determine the optimal value. The sliding window with the optimal window size is used to search for local maximum values. Finally, the coordinates of the identified individual trees are output and marked on the DOM.

[0119] (3) Invert the pixel values ​​of the canopy height model, mark the treetop points identified by the local maximum algorithm as the foreground, mark the non-canopy areas determined by threshold segmentation as the background, and mark the remaining uncertain parts as 0; then use the watershed algorithm. Each time water is poured, the labels will be updated. When two different colored labels meet, a watershed dam will be built, eventually forming a closed outline; each local minimum and its influence range (single tree canopy area) form a water accumulation basin, and the boundary of the water accumulation basin is the watershed (single tree canopy boundary).

[0120] (4) Take the projection and affine transformation of the input canopy height model, add projection coordinates to the image, and save it as a raster image in tif format; convert the raster data of the raster image into vector surface data according to the tree ID. There will be some small patches in the converted data. Merge these small patches into the nearest large patch.

[0121] S5. Vectorize the spatial latitude and longitude projection coordinates of individual trees obtained in S4 (i.e., the spatial location of individual trees in the whole stand) to establish a database of the spatial location and tree height of individual trees in the target stand; use the Anselin Local Moran's I method to analyze the spatial clustering of individual tree height in the target stand and evaluate the spatial heterogeneity of the target stand in a graded manner (the evaluation level can be adjusted according to the actual situation, and this embodiment is customized to 5 levels).

[0122] In practical use, the steps are as follows:

[0123] Using ArcGIS software, local Moran's index analysis was performed on the obtained tree height values ​​of the target forest stand to obtain the clustering patterns of high and low tree height values. Based on the proportion of non-clustered tree points in non-high-high and low-low clustering patterns to the total number of tree points, the spatial heterogeneity of the target forest stand site quality was rated according to five levels.

[0124] S6. Based on the spatial heterogeneity rating results obtained in S5, and taking into account the spatial requirements such as the number of trees and the size of the tree crowns in the target forest stand, as well as the influence of micro-topography (i.e., topography at a relatively small scale) on the trees, different grid sizes with different levels of fineness are set for the target forest stand.

[0125] In practical use, the grid size division in this step needs to be the same as the number of levels in S5. In this embodiment, it is customized to five grid sizes corresponding to the S5 rating results. For example, the spatial heterogeneity of the target forest stand site quality is rated according to I, II, III, IV, and V using the proportions [0%, 20%), [20%, 40%), [40%, 60%), [60%, 80%), and [80%, 100%). The corresponding forest stand division uses five grid sizes: 15m×15m, 12m×12m, 9m×9m, 6m×6m, and 3m×3m.

[0126] S7. Using the vectorized valleys and ridges from S3 as micro-topographic boundary areas, overlay the corresponding refined grids from S6. Referring to the normalized LiDAR point cloud data obtained in S2, perform gridding of the valleys and ridges within the target forest stand range, and vectorize the slope areas based on the valleys and ridges. Obtain and label the gridded valley areas, ridge areas, and slope areas of the target forest stand, and obtain the vectorized grid of the target forest stand.

[0127] S8. Based on the individual tree spatial location and tree height database obtained in S5 and the target forest stand vectorized grid obtained in S7, establish a database of individual tree height and number within each grid, and further statistically select the average height SH of the dominant trees within each grid. t .

[0128] In practical application, the individual trees within each grid are ranked by height, and the top 20% are selected as the dominant trees of the target species. The average height SH of the dominant trees in that grid is then calculated. t .

[0129] S9, based on the average height of the dominant trees in each grid in S8, SH t By referring to the status index table of Chinese fir plantations in the relevant province and city, the grid-scale fine status index of the target forest stand was obtained; then, combined with the remote sensing image of the target forest stand, the grid-scale fine status index was visualized and mapped.

[0130] The methods and effects of the present invention will be specifically illustrated below through examples.

[0131] Example

[0132] This embodiment provides a precise evaluation method for the fine-grained status index of Chinese fir plantations based on UAV-LiDAR. The specific steps are as follows:

[0133] S1. Collection of LiDAR point cloud data for the target area: The study area is located in Quankeng Village, Kaihua Forest Farm, Kaihua County, Quzhou City, Zhejiang Province, a 25-year-old plantation of Chinese fir, with a study area of ​​9.10 hm². 2 On February 14, 2023, a DJI Matrice 300RTK equipped with a Zenmuse L1 lidar system was used to conduct field operations in the studied forest stand area to obtain lidar point cloud data. The flight parameters were: terrain-following flight mode, relative flight altitude of 150m, lidar forward overlap of 20%, visible light lateral overlap of 70%, FOV of 70°×74.5°, and point cloud density of 360 points / m². 2 DJI Terra software (Advanced version) was used to perform 3D reconstruction of the point cloud data of the target area collected in S1, obtaining LiDAR point cloud data (such as...). Figure 2 CloudCompare software was used for manual noise reduction of non-target ground features (high-voltage lines, towers, utility poles, etc.), and LiDAR 360 software was used for noise reduction of drift noise points.

[0134] S2 employs the Cloth Simulation Filter (CSF) algorithm to separate ground points from non-ground points in the 3D LiDAR point cloud data obtained in S1. The kriging interpolation method is then used to generate a digital elevation model (DEM) from the separated ground points in the point cloud. Figure 3 A) and Digital Surface Models (DSMs, such as Figure 3 B); Based on the generated DEM and DSM, the canopy height model (CHM) of the sample plot is obtained by performing interpolation calculation. Figure 3 C) Normalize the LiDAR point cloud to obtain normalized point cloud data (TNPC, e.g.) Figure 3 D).

[0135] S3. Analyze and create positive, negative, and inverted terrain features from the DEM obtained in S2, extracting and vectorizing the ridges and valleys of the target forest area (e.g., Figure 4 ).

[0136] S4. Using the marked watershed control method and the local maximum method, the canopy height model obtained in S2 is segmented into individual trees (until the segmentation accuracy reaches 90%) to obtain the spatial latitude and longitude projection coordinates and tree height of individual trees in the entire stand. The spatial location of individual trees in the entire stand is visualized using ArcGIS software, and the spatial locations of individual trees in the entire stand are marked as points (e.g., ...). Figure 5 ).

[0137] S5. Vectorize the spatial locations of individual trees in the entire stand (i.e., latitude and longitude projected coordinates) obtained in S4 to establish a database of spatial locations and tree heights of individual trees in the entire stand. Use the Anselin Local Moran's I method to analyze the spatial clustering of tree heights in the target Chinese fir forest, obtaining a local Moran's I map of tree heights in the target stand (e.g., ...). Figure 6 Then, the spatial heterogeneity of the target forest stand was evaluated at different levels (defined as 5 levels).

[0138] S6. Statistical analysis of the target forest stand revealed 6168 clusters of high and low tree height values, and 9619 no-cluster points, accounting for 60.98% of all individual tree points. This falls within the range of [60%, 80%). Therefore, the site quality spatial heterogeneity rating of the target forest stand is IV, and a 6m × 6m grid should be used for stand division.

[0139] S7. Using the valleys and ridges obtained in S3 as the micro-topographic boundary areas, and employing a 6m×6m grid, and referring to the normalized point cloud data obtained in S2, the valleys and ridges within the target forest stand are gridded to obtain the valley and ridge grids of the target forest stand; the target forest stand vectors are then cut according to the valley and ridge grids to obtain the slope vectors. A schematic diagram of the target forest stand grid is obtained (e.g., ...). Figure 7 The target forest stand comprises 2,800 grids.

[0140] S8. Based on the individual tree heights and stand grids obtained in S4 and S6, rank the individual trees in each grid by height, select the top 20% as the dominant trees in that grid, and calculate the average height (SH) of the dominant trees in each grid. t .

[0141] S9, based on the average height of the dominant trees in each grid obtained from S8, SH tThe target forest is located in Kaihua County, Zhejiang Province, with a forest age of 25 years. Referring to the "Status Index Table of Chinese Fir Seedling Forests in Zhejiang Province" on page 130 of the "Commonly Used Number Tables for Forestry Exploration and Design," the fine-grained status index of the target forest stand at the grid scale was obtained. For the 6m×6m grid scale, there are 11 instances with a status index of 6, 13 instances with 8, 28 instances with 10, 136 instances with 12, 370 instances with 14, 672 instances with 16, 726 instances with 18, 529 instances with 20, 179 instances with 22, 14 instances with 24, and 122 instances with "no standing trees."

[0142] S10. Based on the position index of each grid obtained in S9, and combined with the remote sensing image of the target forest stand, produce a fine-grained position index thematic map at the grid scale (e.g., Figure 8 ).

[0143] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. A precise evaluation method for the fine-grained status index of Chinese fir plantations based on UAV-LiDAR, characterized in that, Specifically as follows: S1. Use UAV-LiDAR technology to collect LiDAR point cloud data of the target Chinese fir plantation, and preprocess the collected LiDAR point cloud data. S2. Perform the operation of separating ground points from non-ground points on the preprocessed LiDAR point cloud data in S1. A digital elevation model (DEM) of the target forest stand is obtained by modeling the separated ground points; a digital surface model of the target forest stand is constructed based on the separated non-ground points; a canopy height model is obtained based on the elevation difference between the DEM and the digital surface model at the same coordinate point; and normalized LiDAR point cloud data is obtained based on the preprocessed LiDAR point cloud data in S1 and the DEM. S3. Analyze and create positive, negative, and reverse terrain data for the digital elevation model described in S2, and extract the ridges and valleys of the target forest and vectorize them. S4. Perform single-tree segmentation on the canopy height model described in S2 based on the marked watershed control method and the local maximum method, and extract the single-tree height to obtain the spatial latitude and longitude projection coordinates and tree height of the target forest stand. S5. Vectorize the spatial latitude and longitude projection coordinates of individual trees described in S4 to establish a database of individual tree spatial locations and tree heights in the target forest stand; use the local Moran index method to analyze the spatial clustering degree of individual tree heights in the target forest stand and conduct a graded evaluation of the spatial heterogeneity of the target forest stand. S6. Based on the spatial heterogeneity rating results obtained in S5, and taking into account the spatial requirements of the target forest stand and the impact of micro-topography on the trees, set different grid sizes for different levels of fineness for the target forest stand. S7. Using the vectorized valleys and ridges in S3 as micro-topographic boundary areas, overlay the refined grids of the corresponding size in S6. Referring to the normalized LiDAR point cloud data obtained in S2, perform gridding of the valleys and ridges within the target forest stand range, and vectorize the slope areas based on the valleys and ridges. Obtain and label the gridded valley areas, ridge areas, and slope areas of the target forest stand, and obtain the vectorized grid of the target forest stand. S8. Based on the individual tree spatial location and tree height database described in S5 and the target forest stand vectorized grid described in S7, establish a database of individual tree height and number within each grid, and further statistically select the average height SH of the dominant trees within each grid. t ; S9, based on the average height of the dominant trees in each grid in S8, SH t By referring to the status index table of Chinese fir plantations in the relevant province and city, the grid-scale fine status index of the target forest stand is obtained; then, combined with the remote sensing image of the target forest stand, the grid-scale fine status index is visualized and mapped. S5 is specifically as follows: Using ArcGIS software, local Moran's index analysis was performed on the obtained tree height values ​​of the target forest stand to obtain the clustering patterns of high and low tree height values. Based on the proportion of non-clustered tree points in non-high-high and low-low clustering patterns to the total number of tree points, the spatial heterogeneity of the target forest stand site quality was rated according to five levels.

2. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR as described in claim 1, characterized in that, In S1, the acquisition of LiDAR point cloud data is achieved by using a DJI Matrice 300 RTK drone equipped with a lidar sensor to perform flight operations on the target forest stand.

3. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR as described in claim 1, characterized in that, In step S1, the preprocessing method for LiDAR point cloud data includes 3D reconstruction and denoising, as detailed below: The advanced version of DJI Terra software was used for 3D reconstruction of LiDAR point cloud data; CloudCompare software was used for manual denoising of non-target ground features; and LiDAR 360 software was used for denoising of drift noise points.

4. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR as described in claim 1, characterized in that, In step S2, the separation of ground points and non-ground points in LiDAR point cloud data is achieved using a cloth-based analog filtering algorithm, as follows: (1) Invert the preprocessed LiDAR point cloud data in S1 and simulate a sufficiently soft cloth to be laid on top of the inverted LiDAR point cloud data. The smallest fiber molecules that make up the simulated cloth are regarded as countless cloth particles. The cloth particles are limited to moving only in the vertical direction. The position and velocity of the cloth particles are determined by the internal and external forces between the cloth particles. The ground points in the LiDAR point cloud data are classified by comparing the height values ​​of the cloth particles and the terrain. (2) Set the grid size and mesh the simulated LiDAR point cloud in step (1) to perform point cloud filtering with the grid as the basic unit; (3) Based on the results of step (2), all LiDAR point clouds and cloth particles in the target area are projected onto the same two-dimensional horizontal plane. The particles on the plane projected by the LiDAR point cloud are called point cloud particles. According to the two-dimensional coordinate distance between the point cloud particles and the cloth particles, the LiDAR point corresponding to each cloth particle in the two-dimensional plane is found. The height of the LiDAR point is defined as IH, and IH is the minimum height that the cloth particle can fall. (4) Based on the results of step (3), calculate the displacement of the fabric particles under the action of external force. Since the particles are not subject to collision force in the air, but only to gravity, calculate the distance the fabric particles move under the action of gravity according to the following formula: ； In the formula, Represents a specific point in time; Represents the time step; This represents the fabric particles under the influence of gravity alone. The location of the time node; This represents the fabric particles under the influence of gravity alone. The location after the time; This represents the value of a fabric particle at that node under the influence of gravity alone. The location at the time; This represents the mass of the fabric particles; g is a constant. (5) Based on the results of step (3), calculate the displacement distance of the fabric particles under the influence of internal forces, as shown in the following formula: ； In the formula, This represents the displacement vector of the cloth particle; b is a constant, with a value of 1 when the cloth particle can move and a value of 0 when it cannot move. The current coordinate vector representing the cloth particle; Represents the coordinate vector of adjacent particles connected to the cloth particle; A unit vector representing the vertical direction; (6) Based on the results of steps (4) and (5), compare the height CH of the cloth particles that have been displaced by external and internal forces with IH; if CH≤IH, place them at the height of IH and set them as immovable particles; if CH>IH, treat them as movable particles and repeat steps (4)~(5) until the maximum change height of all cloth particles is small enough or immovable or exceeds the set maximum number of iterations. At this time, the plane composed of all cloth particles is an approximate real terrain, and the cloth simulation iteration is completed. (7) Based on the results of step (6), calculate the distance MH between the cloth particles and the corresponding LiDAR points after the cloth simulation is completed; (8) Based on the results of step (7), set the distance threshold OH; compare the magnitudes of MH and OH for each LiDAR point. When MH≤OH, the LiDAR point is marked as a ground point; when MH>OH, the LiDAR point is marked as a non-ground point.

5. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR as described in claim 1, characterized in that, In S2, the digital elevation model and digital surface model are constructed using the Kriging interpolation method, as shown in the following formula: ； In the formula, The number of points that can determine the attribute value of a point; Points to be estimated The estimated value; Representing the The attribute values ​​of a known point; Representing the The weights of the known points; ； In the formula, Indicates the first , The semivariance between the domain points; Representing the The attribute values ​​of a known point; Choose an appropriate function form and fit the function. ; Given the location of the point to be found, we can obtain... Substituting into the function yields Calculate matrix multiplication: ； In the formula, It is a constant; the 0th neighborhood point represents the point to be found itself; Finally, the points to be estimated were completed. Valuation.

6. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR as described in claim 1, characterized in that, S3 is specifically as follows: (1) Using ArcGIS software, based on the digital elevation model described in S2, obtain statistical raster data Focal_St through the focus statistics method in terrain domain analysis; calculate the difference between the digital elevation model and Focal_St at the same spatial location and obtain raster data R1_Ca; reclassify R1_Ca with 0 as the critical value, and treat R1_Ca above 0 as positive terrain raster Reclassify1, and R1_Ca below 0 as negative terrain raster Reclassify2; (2) Based on the digital elevation model described in S2, perform raster calculations and obtain reverse terrain raster data Reclassify3; (31) Based on the digital elevation model described in S2, perform depression filling, flow direction, flow rate analysis and raster calculation of the topographic hydrology to obtain the raster area data R2_Co with a flow rate of 0; determine the ridge classification threshold of R2_Co by analyzing the surface contour lines and mountain shadows of the digital elevation model data; reclassify R2_Co according to the classification threshold and obtain Reclassify4; (32) Multiply Reclassify1 and Reclassify4 to obtain the raster data Reclassify; the raster with a reclassification retention value of 1 is the ridge raster area, and the ridge can be extracted by vectorizing the result; (41) Based on Reclassify3, perform depression filling, flow direction, flow rate analysis and raster calculation for the topographic hydrology to obtain the raster area data R3_Co with a flow rate of 0; determine the valley classification threshold of R3_Co by analyzing the surface contour lines and mountain shadows of the digital elevation model data; reclassify R3_Co according to the classification threshold and obtain Reclassify5; (42) Multiply Reclassify2 and Reclassify5 to obtain raster data Reclassify6; the raster with a reclassification retention value of 1 is the valley raster area, and the valley can be extracted by vectorizing the result; (5) Multiply Reclassify by Reclassify6 to obtain raster data Reclassify7; the raster with a reclassification retention value of 0 is the hillside raster area, and the hillside can be extracted by vectorizing the result.

7. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR as described in claim 1, characterized in that, S4 is specifically as follows: (1) Based on the height variation characteristics of a single tree crown, the local maximum algorithm is used to identify the treetop point of a single tree; the identification result is regarded as the X and Y coordinates of the tree, and the foreground is marked with it to optimize the watershed segmentation effect; (2) The local maximum detection method based on the canopy height model is used to identify individual trees. During the iteration process, the window size is dynamically adjusted to determine the optimal value. The sliding window with the optimal window size is used to search for local maximum values. Finally, the coordinates of the identified individual trees are output and marked on the DOM. (3) Invert the pixel values ​​of the canopy height model, mark the treetop points identified by the local maximum algorithm as the foreground, mark the non-canopy areas determined by threshold segmentation as the background, and mark the remaining uncertain parts as 0; then use the watershed algorithm. Each time water is irrigated, the labels will be updated. When two different colored labels meet, a watershed dam will be built, eventually forming a closed outline; each local minimum and its influence range form a water basin, and the boundary of the water basin is the watershed; (4) Take the projection and affine transformation of the input canopy height model, add projection coordinates to the image, and save it as a raster image in tif format; convert the raster data of the raster image into vector surface data according to the tree ID. There will be some small patches in the converted data. Merge these small patches into the nearest large patch.

8. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR as described in claim 1, characterized in that, In S6, five different grid sizes with varying degrees of fineness are defined to correspond to the five rating results in S5.

9. The method for precise evaluation of the fine position index of Chinese fir plantation land based on UAV-LiDAR according to claim 1, characterized in that, In step S8, the individual trees in each grid are arranged by height, and the top 20% are selected as the dominant trees of the target tree species. The average height SH of the dominant trees in that grid is then calculated. t .

Images