Method and system for reconstructing laser radar point cloud image of three-dimensional structure of forest canopy

Abstract

The present application relates to the technical field of remote sensing image processing, and discloses a forest canopy three-dimensional structure laser radar point cloud image reconstruction method and system, which comprises: adaptive terrain perception point cloud preprocessing, multi-scale vertical layered canopy point extraction, energy functional driven single tree segmentation, voxelized canopy structure parameter inversion, and octree accelerated three-dimensional reconstruction and visualization. Through slope adaptive morphological filtering, vertical density histogram adaptive layering, marker control watershed and graph cut energy optimization coupling, voxelized leaf area density inversion and closed loop feedback mechanism, the present application realizes high-precision automatic reconstruction of the forest canopy three-dimensional structure.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing and forestry information extraction technology, and in particular to a method and system for reconstructing three-dimensional structure lidar point cloud images of forest canopy. Background Technology

[0002] The three-dimensional structure of forest canopy is a core indicator characterizing the functional status of forest ecosystems, and its spatial distribution characteristics are directly related to photosynthetic efficiency, carbon cycle processes, and biodiversity assessment. Traditional methods for obtaining forest canopy structure parameters mainly rely on manual field surveys and plot measurements. This approach not only consumes a lot of manpower and resources, but also makes it difficult to achieve full coverage when dealing with large forest areas. In particular, the accessibility of field measurements is severely limited in mountainous forest areas and tropical rainforests with complex terrain conditions, resulting in a long-term lack of canopy structure information in key areas.

[0003] With the rapid development of airborne lidar technology, obtaining three-dimensional structural information of forest canopies using lidar point cloud data has become an important research direction in the field of forestry remote sensing. Among existing technologies, Chinese patent CN113570621A discloses a method for extracting tree information based on high-precision point clouds and imagery. This method combines two-dimensional digital orthophoto maps with three-dimensional lidar point cloud data, extracting tree coordinates, height, and crown width through two single-tree segmentation processes: canopy height model segmentation and point cloud segmentation. However, this solution is mainly aimed at the tree information extraction needs in power grid corridor scenarios, and its technical framework has the following limitations: First, the method uses a watershed segmentation algorithm with fixed parameters to segment the canopy height model, without considering the impact of terrain slope changes on the filtering window and segmentation accuracy, making it insufficiently applicable in complex mountainous forest areas; Second, the method only extracts two-dimensional planar projection indicators such as tree height, canopy width, and canopy area, lacking the ability to invert three-dimensional refined parameters such as canopy vertical layering structure and leaf area density distribution; In addition, the method lacks a feedback mechanism between ground point filtering and canopy point extraction, and cannot dynamically optimize the front-end preprocessing parameters based on subsequent canopy structure analysis results, limiting the overall reconstruction accuracy of the system.

[0004] In single-tree segmentation, existing methods can be broadly categorized into three types: raster-based methods, point cloud-based methods, and joint methods. Raster-based methods convert 3D point clouds into 2D canopy height models and then use image processing techniques for segmentation. While computationally efficient, they inevitably lose vertical structural information during dimensionality compression. Point cloud-based methods directly utilize raw or voxelized 3D point clouds for segmentation, better preserving spatial structural details, but often suffer from complex parameters and insufficient generalization ability. Joint methods attempt to combine the advantages of the above two types, but existing joint schemes mostly employ simple cascading approaches, lacking a unified segmentation framework based on global energy optimization.

[0005] In terms of canopy structure parameter inversion, voxelization methods have been used to estimate the vertical distribution of leaf area density. However, existing voxelization schemes have not yet formed an effective closed-loop collaborative mechanism with the individual tree segmentation and canopy point extraction processes. Furthermore, most existing 3D visualization schemes are based on direct point cloud rendering, lacking the support of multi-resolution spatial indexing structures. This results in insufficient rendering efficiency and poor interactive smoothness when dealing with massive point cloud data from large forest areas. In summary, there is an urgent need for a forest canopy 3D structure reconstruction scheme that can achieve deep coupling of terrain-adaptive preprocessing, vertically layered canopy point extraction, energy-optimized individual tree segmentation, and voxelization structural parameter inversion. Summary of the Invention

[0006] To address the technical problems of existing methods for reconstructing the three-dimensional structure of forest canopies, such as poor terrain adaptability, lack of vertical layering information, limited segmentation accuracy, and lack of closed-loop coordination among processing stages, this invention provides a method and system for reconstructing the three-dimensional structure of forest canopies from lidar point cloud images. The technical solution of this invention achieves automated, high-precision reconstruction of the three-dimensional canopy structure model from raw airborne lidar point cloud data by constructing a five-step deeply coupled closed-loop architecture.

[0007] The forest canopy three-dimensional structure lidar point cloud image reconstruction method provided by this invention includes the following steps: an adaptive terrain-aware point cloud preprocessing step, which acquires discrete point cloud data of the target forest area scanned by an airborne lidar, performs noise filtering and coordinate unification processing on the point cloud data, uses a slope-adaptive progressive morphological filtering algorithm to separate ground points and vegetation points in the processed point cloud, and uses ground points to generate a digital elevation model through irregular triangular mesh interpolation; a multi-scale vertically layered canopy point extraction step, which performs height normalization on vegetation points based on the digital elevation model to obtain a normalized point cloud, constructs a vertical height distribution density histogram, adaptively determines the layering thresholds for shrub and tree layers based on density valley values, and removes shrub layer point clouds to obtain a tree layer canopy point set; and an energy functional-driven single... The tree segmentation step involves generating a canopy height model based on the canopy point set, extracting labeled seed points using variable window local maximum detection, performing labeled control watershed segmentation to obtain the initial segmented region, and performing graph cut optimization on the initial segmented region by minimizing the energy functional to obtain the individual tree segmentation results. The voxelized canopy structure parameter inversion step involves performing 3D voxelization on each individual tree point cloud, calculating the leaf area density of each layer based on a laser pulse penetration probability model to obtain the leaf area density vertical profile, extracting tree height, canopy width, canopy length, and canopy volume parameters, and feeding the leaf area density vertical profile back to the canopy point extraction step to optimize the layering threshold. The octree accelerated 3D reconstruction and visualization step involves establishing an octree spatial index and performing multi-resolution mesh reconstruction on the forest stand canopy point cloud to generate a 3D canopy model.

[0008] The present invention also provides a forest canopy three-dimensional structure lidar point cloud image reconstruction system, including an adaptive terrain-aware point cloud preprocessing module, a multi-scale vertical layered canopy point extraction module, an energy functional-driven single tree segmentation module, a voxelized canopy structure parameter inversion module, and an octree-accelerated three-dimensional reconstruction and visualization module, each module corresponding to the steps of the above method.

[0009] The beneficial effects of this invention are as follows: It achieves adaptive processing of complex terrain through slope-adaptive progressive morphological filtering, maintaining high accuracy in ground point classification under different slope conditions, with tree height measurement accuracy better than 0.5m; it achieves automatic stratification of shrub and tree layers through vertical height density histogram analysis, avoiding the limitations of manually setting fixed thresholds; it achieves single-tree segmentation accuracy greater than 85% through deep coupling of marker-controlled watersheds and graph cut optimization, effectively suppressing over-segmentation and under-segmentation problems under dense canopy conditions; it obtains fine vertical structure information of the canopy through voxelized Beer-Lambert leaf area density inversion; and it drives dynamic optimization of the canopy point extraction threshold through leaf area density profile feedback, forming a complete closed-loop collaborative mechanism, improving the canopy point complete extraction rate by approximately 5 percentage points, and significantly outperforming existing serial processing schemes in overall reconstruction accuracy and efficiency. Attached Figure Description

[0010] Figure 1 This is a flowchart of a method for reconstructing a three-dimensional structure of a forest canopy using lidar point cloud images, provided in an embodiment of the present invention.

[0011] Figure 2 This is an architecture diagram of the forest canopy three-dimensional structure lidar point cloud image reconstruction system provided in an embodiment of the present invention. Detailed Implementation

[0012] To make the technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the following embodiments are only for explaining the present invention and are not intended to limit the scope of protection of the present invention. Without departing from the spirit and scope of the technical solutions of the present invention, those skilled in the art can make reasonable modifications and variations to the specific parameters and implementation details in the following embodiments, and all such modifications and variations should fall within the scope of protection of the present invention.

[0013] See Figure 1 The forest canopy 3D structure lidar point cloud image reconstruction method provided in this invention includes an adaptive terrain-aware point cloud preprocessing step, a multi-scale vertically layered canopy point extraction step, an energy functional-driven single-tree segmentation step, a voxelized canopy structure parameter inversion step, and an octree-accelerated 3D reconstruction and visualization step. These five steps form a deeply coupled data flow transmission relationship: the output of each preceding step directly constitutes the core input of subsequent steps, while the output of the voxelized canopy structure parameter inversion step is also fed back to the multi-scale vertically layered canopy point extraction step through a reverse feedback path to drive the dynamic optimization of key parameters, thus forming a closed-loop collaborative processing architecture that unifies feedforward transmission and feedback adjustment. Each step is described in detail below.

[0014] Step S1: Adaptive terrain-aware point cloud preprocessing step. The core purpose of this step is to accurately separate ground points and vegetation points from the discrete point cloud data collected by the original airborne lidar, and generate a high-precision digital elevation model as the basis for subsequent canopy point height normalization.

[0015] In one embodiment of the present invention, discrete point cloud data is first acquired by an airborne lidar scanner performing flight strip scanning on a target forest area. Preferably, the point density of the point cloud data is not less than 4 points / m. 2 This ensures sufficient spatial sampling resolution for subsequent extraction of canopy structure parameters. The lidar system's echo recording mode supports multiple echo separation, meaning it can record the three-dimensional coordinates and intensity information of the first, intermediate, and final echoes for each transmitted pulse. In this embodiment, the scanning height is set to 500m to 1500m, the scanning angle does not exceed ±20°, and the flight strip overlap rate is not less than 30%, thereby ensuring the integrity and uniformity of the point cloud coverage.

[0016] After acquiring the raw point cloud data, noise filtering is performed. Specifically, an outlier removal method based on elevation statistics is adopted. First, the point cloud space is divided into several statistical units, preferably with the horizontal dimensions of each statistical unit set to 50m × 50m. Then, the average elevation of all point clouds within each statistical unit is calculated. and elevation standard deviation Points whose elevation values ​​deviate from the local mean by more than three standard deviations are marked as noise points and removed, thus satisfying the condition. point Considered as noise points; among them, For the first Elevation values ​​of each point This is the average elevation of all points within this statistical unit. The standard deviation of elevation for all points within this statistical unit is given in meters. This method effectively removes anomalous echo points caused by atmospheric scattering, bird flight, and system noise, while preserving the point cloud information of real ground features.

[0017] After noise filtering, the point cloud data undergoes coordinate unification processing, transforming all flight strip data to a unified geographic coordinate system. Preferably, the WGS84 coordinate system combined with UTM projection is used, and accurate registration between flight strips is achieved by matching corresponding points in the overlapping areas between flight strips, with the registration residual controlled within 0.1m.

[0018] Based on coordinate unification, a slope-adaptive progressive morphological filtering algorithm is employed to separate ground points from vegetation points in the point cloud. The core idea of ​​this algorithm is to dynamically adjust the morphological filtering window size according to the local terrain slope gradient to adapt to the differentiated needs of ground point recognition under different terrain conditions. In a preferred embodiment of this invention, a coarse-grid digital elevation model is first used to preliminarily estimate the local terrain slope of each region. Then, the filter window size is determined according to the following adaptive rules. : ,in, This represents the size of the morphological filtering window, in meters (m). Local terrain slope, in degrees; The slope threshold for gentle terrain is set to 5°. The slope threshold for steep terrain is set to 15°. The maximum window size for flat terrain is set to 20m. The minimum window size for steep terrain is set to 5m. This parameter is chosen because: in gently sloping terrain, ground points are continuously distributed with minimal undulation, and a larger filter window helps to completely identify the ground surface; while in steep terrain, ground elevation changes drastically, and an excessively large window can lead to misclassifying abrupt changes as vegetation points. Therefore, a smaller window is needed to improve the sensitivity of ground point identification. This adaptive mechanism ensures that the algorithm achieves an overall ground point classification accuracy of no less than 95% under complex terrain conditions.

[0019] Using the separated ground points, a digital elevation model (DEM) is generated through an irregular triangular mesh (ITM) interpolation algorithm. Preferably, the grid resolution of the ITM interpolation is set to 0.5m to 1m. In this embodiment, a Delaunay triangular mesh is constructed for the ground points, and then linear interpolation is used to generate regular grid DEM raster data. The generated DEM will serve as the reference datum for vegetation point height normalization in the subsequent step S2.

[0020] Step S2: Multi-scale vertical stratified canopy point extraction step. The core task of this step is to accurately extract the tree canopy point set from the preprocessed vegetation point cloud, while adaptively removing the point cloud interference from the understory shrub and herb layers, so as to provide high-quality input data for subsequent single tree segmentation and canopy structure parameter inversion.

[0021] In one embodiment of the present invention, height normalization is first performed on all vegetation points based on the digital elevation model generated in step S1. The specific method of height normalization is as follows: for each point in the vegetation point set... According to its horizontal coordinates In the digital elevation model, the ground elevation value at the corresponding location is obtained through bilinear interpolation. Then use the original elevation value of that point. Normalized height is obtained by subtracting the ground elevation value. : ,in, For the first Normalized height of each vegetation point, in meters; This is the original elevation value of the point, in meters (m). For digital elevation models in coordinate The interpolated elevation at the location is in meters. The normalized point cloud height directly reflects the absolute height of the vegetation relative to the ground, eliminating the influence of terrain undulation on height measurement.

[0022] After height normalization, a vertical height distribution density histogram of the normalized point cloud is constructed. Preferably, the height layer spacing is set to 0.5m, that is, the vertical space from 0m to the maximum normalized height is divided into several horizontal thin layers at 0.5m intervals, and the number of point clouds contained in each thin layer is counted as the point density value of that layer. The vertical height distribution density histogram is constructed by arranging the point density values ​​of each layer in height order.

[0023] In a preferred embodiment of the present invention, the vertical height distribution density histogram is Gaussian smoothed to eliminate the influence of random fluctuations, and the standard deviation of the Gaussian kernel is set to 1m. The minimum point is extracted from the smoothed density histogram as the density valley value. In a typical multi-layered forest stand structure, there is usually a distinct point cloud sparse band between the shrub layer and the tree layer. This sparse band corresponds to the density valley value position in the density histogram; therefore, the height corresponding to this valley value can be used as the stratification threshold between the shrub layer and the tree layer. In this embodiment, the initial value of the stratification threshold is set to 2m to 4m, and the specific value is adaptively determined by the valley position of the density histogram.

[0024] It should be noted that, within the technical framework of this invention, the layering threshold... The parameters are not fixed, but are dynamically optimized using the leaf area density vertical profile obtained from the voxelized canopy structure parameter inversion in subsequent step S4. The specific feedback optimization mechanism will be described in detail in the description of step S4.

[0025] The normalized point cloud is filtered using a layering threshold, removing points with a normalized height below the layering threshold. The point clouds were labeled as shrub layer point clouds and removed, while the point clouds with normalized height greater than or equal to the layering threshold were retained as the tree canopy layer point set. This canopy point set will serve as the input data for individual tree segmentation in step S3.

[0026] Step S3: Energy functional-driven individual tree segmentation step. The core objective of this step is to accurately segment the point cloud of each tree from the overall stand point cloud from the canopy point set output in Step S2, laying the foundation for subsequent canopy structure parameter inversion at the individual tree level. Unlike existing technologies that simply use watershed algorithms or point cloud clustering algorithms for individual tree segmentation, this invention proposes a two-stage individual tree segmentation strategy that couples label-controlled watershed pre-segmentation with graph cut energy optimization, which can significantly improve the segmentation accuracy under dense forest stand conditions.

[0027] In one embodiment of the present invention, the canopy point set obtained in step S2 is first used as a basis. Generate a canopy height model. Specifically, project the canopy point set onto a horizontal plane, mesh the projection plane with a set grid resolution (preferably 0.5m), and then take the maximum normalized height value of the point cloud falling within each grid cell as the canopy height model value for that grid. For blank grid cells without point cloud coverage, neighbor interpolation is used to fill them.

[0028] On the generated canopy height model, a variable window local maximum detection algorithm is used to extract candidate treetop points as seed points. The core idea of ​​this algorithm is to dynamically determine the search window radius based on the height values ​​at various locations in the canopy height model, adapting to differences in canopy width under different tree height conditions. Preferably, the search window radius... height value at corresponding position The relationship between them was determined using a linear regression model: ,in, The radius of the search window is in meters. This represents the height value at the corresponding location in the canopy height model, in meters (m). The regression coefficient has a value range of 0.04 to 0.08, preferably 0.06. Its physical meaning is the proportional relationship between the crown radius and the tree height. The selection of this parameter is based on the statistical regression relationship between crown width and tree height in typical temperate and subtropical forest stands. The regression intercept, ranging from 0.5m to 1.5m, preferably 1.0m, represents the base radius of the minimum search window. Under this variable window mechanism, tall trees are assigned larger search windows to avoid misidentifying multiple local high points within the same canopy as different treetops, while short trees are assigned smaller search windows to ensure the distinguishability between adjacent trees. For local maximum points that meet the conditions, a minimum height constraint is further applied, meaning the normalized height of the candidate treetop point is no less than 5m, to exclude low shrub remnants and noise spurious peaks.

[0029] Using the extracted candidate treetop points as seed points, a label-controlled watershed segmentation algorithm is applied to the canopy height model to obtain the initial segmentation region. Label-controlled watershed segmentation differs from traditional unlabeled watershed segmentation in that traditional watershed segmentation is prone to severe oversegmentation, while label-controlled watershed segmentation effectively suppresses oversegmentation by pre-specifying seed point markers to constrain the number of catchment basins in the watershed. In this embodiment, the gradient map of the canopy height model is used as input for watershed segmentation. Each seed point corresponds to a catchment basin, and the boundary of the catchment basin is the initial segmentation boundary between adjacent tree canopies.

[0030] However, relying solely on marker-controlled watershed segmentation still falls short of completely resolving the problem of accurately delineating boundaries between dense canopies, especially in broadleaf forests where adjacent canopies are close in height and significantly overlap. Therefore, this invention further introduces a graph cut energy optimization mechanism to refine the initial segmentation results. Specifically, firstly, a regional adjacency graph is constructed based on the initial segmented regions, weighted by point cloud spatial distance and height differences. , where the set of nodes The edge set corresponding to each region generated by the initial segmentation Connect adjacent regions in the space. Based on this, construct an energy functional that includes data terms and smoothing terms. : ,in, Total energy, dimensionless; The total number of point clouds or regions participating in the optimization; For the first A feature vector of a point cloud or region, containing three-dimensional coordinates and height information; To be assigned to the Individual wooden labels for points or areas; For data items, measure the points Assign to tag The smaller the spatial distance cost of the corresponding single tree center, the better the geometric consistency between the point and its corresponding single tree. The smoothing weight coefficient, ranging from 0.1 to 1.0, preferably 0.5, is used to balance the relative contributions of data items and smoothing items. When the value is too small, the optimization results tend to be independent classifications but lack spatial smoothness. When the size is too large, the optimization result becomes too smooth and loses local details; It is a set of spatial adjacency relationships; For the smoothing term, the Potts model is used for definition, when hour ,when hour It is used to penalize inconsistencies in labels between adjacent regions.

[0031] By employing a graph cut algorithm to obtain a globally approximate optimal solution for the aforementioned energy functional, the final individual tree assignment labels for each point cloud or region are determined, thereby achieving refined individual tree segmentation results. In a preferred embodiment of the invention, the graph cut optimization problem is solved iteratively using an alpha-beta exchange or alpha extension algorithm, with the iterative convergence criterion set at an energy change rate of less than 0.1% between adjacent iterations. Experiments show that the individual tree segmentation accuracy after graph cut energy optimization is greater than 85%, an improvement of approximately 8 to 12 percentage points compared to the scheme using only label-controlled watershed segmentation. This improvement mainly stems from the graph cut optimization's ability to accurately adjust the boundaries of densely overlapping canopies.

[0032] Step S4: Voxelized canopy structure parameter inversion step. The core task of this step is to perform refined three-dimensional canopy structure parameter inversion on the individual tree segmentation results output in step S3, obtain multi-dimensional canopy structure parameters including the vertical distribution of leaf area density, and drive the dynamic optimization of the canopy point extraction threshold in step S2 through a closed-loop feedback mechanism.

[0033] In one embodiment of the present invention, the point clouds of each individual tree segmentation result are first subjected to three-dimensional voxelization processing. The specific method of voxelization processing is as follows: using the three-dimensional bounding box of the individual tree point cloud as the spatial range, the bounding box is uniformly meshed according to a set voxel size, dividing the space into a regularly arranged three-dimensional voxel array. Preferably, the voxel size is set between 0.5m×0.5m×0.5m and 1m×1m×1m, and in this embodiment, 1m×1m×1m is preferred. The selection of voxel size needs to strike a balance between spatial resolution and statistical reliability: while a too small voxel size can provide higher spatial resolution, the number of points in each voxel may be insufficient to support reliable statistical estimation; while a too large voxel size will cause the fine structural information inside the canopy to be smoothed out. For a point density of not less than 4 points / m... 2 For airborne lidar data, a voxel size of 1m ensures that most voxels contain a sufficient number of laser pulse records for statistical calculation of penetration probability.

[0034] After voxelization, the leaf area density of each horizontal voxel array is calculated based on a laser pulse penetration probability model. This laser pulse penetration probability model is based on a modified Beer-Lambert law, the basic principle of which is that when a laser pulse passes through a voxel containing vegetation material, the probability of it being intercepted is proportional to the plant area density within that voxel. The specific calculation formula is as follows: ,in, For the first Leaf area density of each horizontal layer, in m³ 2 / m 3 ; The thickness of the voxel in the vertical direction, i.e., the horizontal layer height interval, is set to 1m in this embodiment; Let be the projection function, representing the blade in the laser incident direction. The projection area ratio coefficient is assumed to be 0.5 for a random leaf tilt angle distribution and is dimensionless. To enter the The total number of laser pulses passing through the voxel array from above is calculated by counting the pulses emitted from above. The number of pulses that were not intercepted was determined; To cross the first The number of laser pulses that did not produce an echo in that layer; It is a function of the natural logarithm. When When the ratio is close to 1, it indicates that the layer has almost no vegetation material intercepting laser pulses, and the corresponding leaf area density is close to 0; when the ratio is small, it indicates that a large number of laser pulses are intercepted, corresponding to a higher leaf area density value. It should be particularly noted that when... 0 or When the value is 0, special handling is required to avoid numerical anomalies in the logarithmic function: when When the leaf area density of this layer is set to null to indicate that it is unobservable; when and When adopted The corrected value is used for substitution calculation.

[0035] By calculating the leaf area density layer by layer at each horizontal level, a complete vertical profile of leaf area density is obtained. ,in This represents the total number of vertical layers. This profile visually reflects the density distribution characteristics of vegetation material within the canopy along the vertical direction and is one of the key parameters describing the fine three-dimensional structure of the canopy.

[0036] While obtaining the vertical profile of leaf area density, conventional structural parameters such as tree height, crown width, crown length, and crown volume are extracted by combining the spatial envelope information of single-tree point clouds. Specifically, tree height... Defined as the difference between the maximum and minimum normalized height in a single tree point cloud; crown width Defined as the smallest circumcircle diameter of the projection of a single tree point cloud onto a horizontal plane, or the average span along two orthogonal directions; crown length. Defined as the vertical span of the canopy point cloud, i.e., the height difference between the highest point and the canopy base. The canopy base position can be determined by the first occurrence of a leaf area density value below a set threshold (preferably 0.1m) in the vertical profile of leaf area density from top to bottom. 2 / m 3 The height of the canopy is determined; the volume of the canopy is determined. It is estimated by summing the number of voxels within the canopy in each horizontal layer and multiplying by the volume of a single voxel.

[0037] The most crucial innovative feature of this step lies in establishing a closed-loop feedback mechanism for the vertical profile of leaf area density in step S2. In the initial round of processing, step S2 determines the stratification thresholds for the shrub and tree layers based on the density valley values ​​of the vertical height distribution density histogram. However, in some complex forest stand structures, such as those with tall shrubs or sub-canopy trees in the understory, the density histogram may not accurately distinguish the boundary between the shrub and tree layers. In such cases, the leaf area density vertical profile obtained in step S4 can provide more refined vertical structure information. The specific feedback optimization rule is: when the leaf area density vertical profile is within the initial stratification threshold... There are still areas with densities greater than the set density threshold within the following height range. At the peak of leaf area density (preferably, Set to 0.3m 2 / m 3 This indicates that the initial stratification threshold may have been set too high, missing some effective point cloud information from the bottom of the canopy. In this case, the stratification threshold is adjusted downwards to the nearest density valley below the peak, triggering a re-execution of steps S2 to S4. The maximum number of iterations for feedback optimization is set to 3 to avoid convergence oscillations caused by excessive iteration. Through this closed-loop feedback mechanism, the system can adaptively handle canopy point extraction under different stand vertical structure conditions, improving the complete canopy point extraction rate by approximately 5% to 8%.

[0038] Step S5: Octree-accelerated 3D reconstruction and visualization. The core task of this step is to transform the individual tree segmentation results and canopy structure parameters obtained in the previous steps into an interactive 3D model of the forest stand canopy, while realizing efficient storage and rendering of large-scale point cloud data through spatial indexing technology.

[0039] In one embodiment of the present invention, an octree spatial index is first established based on the single-tree segmentation results output in step S3 and the canopy structure parameters extracted in step S4. An octree is a hierarchical spatial index structure that recursively divides a three-dimensional space into eight equal parts. It can adaptively adjust the index granularity according to the distribution density of point clouds in different spatial regions, thereby effectively controlling the storage overhead of the index while ensuring query efficiency. Preferably, the maximum recursion depth of the octree spatial index is set to 8 to 12 layers, and in this embodiment, 10 layers are preferred. The selection of the recursion depth is based on typical forest stand point cloud data (coverage area 100 ha, point density 4 to 10 points / m²). 2A recursion depth of 10 levels can control the spatial size of the smallest leaf node to the order of approximately 0.1m, which satisfies the accuracy requirements for point cloud detail representation without generating too many empty nodes that would cause index bloat. Each leaf node of the octree stores the 3D coordinates of the point cloud it contains, the label of the individual tree it belongs to, and the associated canopy structure parameters.

[0040] Based on the established octree spatial index, multi-resolution mesh reconstruction is performed on the canopy point cloud to generate a 3D canopy model. Preferably, the multi-resolution mesh reconstruction employs a Poisson surface reconstruction algorithm based on normal vector estimation. This algorithm first estimates the local normal vector direction for each point in the point cloud, then transforms the surface reconstruction problem into solving a Poisson equation. The implicit surface function is obtained by solving the Poisson equation, and finally, an explicit triangular mesh surface is generated using an isosurface extraction algorithm. Preferably, the reconstruction depth of the Poisson surface reconstruction is set to 8 to 10, and in this embodiment, 9 is preferred. The reconstruction depth determines the spatial resolution of the output mesh; a larger depth value can generate a finer mesh surface, but the computational cost increases accordingly. In practical applications, a higher reconstruction depth is used for detailed close-range observation of local areas to present the detailed texture of the canopy surface, while the reconstruction depth is automatically reduced for distant areas from a global overview perspective to reduce the rendering load, thereby achieving multi-resolution hierarchical rendering.

[0041] The generated 3D canopy model supports multiple interactive analysis functions. First, it supports arbitrary viewpoint navigation, allowing users to observe the 3D spatial distribution characteristics of the forest canopy from any angle through translation, rotation, and zoom operations. Second, it supports vertical profile analysis of the canopy along a specified direction: users can specify a profile line on the horizontal plane, and the system extracts point clouds and canopy structure parameters along this line, displaying information such as canopy height changes, vertical distribution of leaf area density, and the location and size of forest gaps in a 2D profile view. Furthermore, each individual tree in the 3D model is associated with a complete set of structural parameters obtained from step S4, including tree height, crown width, crown length, canopy volume, and vertical profiles of leaf area density. Users can query detailed parameter information for any individual tree through interactive clicking. This visualization function provides an intuitive and efficient tool for browsing and analyzing canopy structure data for applications such as forest volume estimation, carbon sink monitoring, and ecological assessment.

[0042] See Figure 2 This invention also provides a forest canopy 3D structure lidar point cloud image reconstruction system. This system includes an adaptive terrain-aware point cloud preprocessing module, a multi-scale vertically layered canopy point extraction module, an energy functional-driven single-tree segmentation module, a voxelized canopy structure parameter inversion module, and an octree-accelerated 3D reconstruction and visualization module. Each module corresponds one-to-one with each step in the above method embodiments, and will be described separately below.

[0043] The adaptive terrain-aware point cloud preprocessing module is configured to receive discrete point cloud data of the target forest area output by an airborne LiDAR scanner as input, perform noise filtering and coordinate unification processing, and use the slope adaptive progressive morphological filtering algorithm described in step S1 of the above method embodiment to separate ground points from vegetation points in the processed point cloud. The output of this module includes two parts: a ground point set and a vegetation point set. The ground point set is processed by irregular triangular mesh interpolation to generate a digital elevation model. In this embodiment, the module further includes a noise filtering subunit, a coordinate transformation subunit, a slope calculation subunit, and a morphological filtering subunit. The noise filtering subunit adopts an outlier removal method based on elevation statistics, identifying and removing noise points by calculating the mean and standard deviation of the point cloud elevation within each statistical unit. The coordinate transformation subunit transforms multi-strip data to a unified geographic coordinate system to ensure spatial consistency in subsequent processing. The slope calculation subunit uses a coarse-grid digital elevation model to estimate the local terrain slope of each area, providing adaptive window size parameters for the morphological filtering subunit. The morphological filtering subunit receives the local slope information output by the slope calculation subunit, and dynamically determines the size of the filtering window according to the piecewise linear adaptive rule described in the above method embodiment, so as to complete the high-precision separation of ground points and vegetation points.

[0044] The multi-scale vertical stratified canopy point extraction module is configured to receive vegetation point sets and digital elevation models output by the adaptive terrain-aware point cloud preprocessing module as input, and perform height normalization processing, vertical height distribution density histogram construction, and adaptive determination of stratification thresholds as described in step S2 of the above method embodiment. The core of this module lies in the vertical density analysis subunit and the threshold optimization subunit. The vertical density analysis subunit performs vertical stratification statistics on the normalized point cloud at a set stratification interval, generates a density histogram, performs Gaussian smoothing, and extracts density valley values. The threshold optimization subunit not only receives the initial threshold output from the vertical density analysis subunit but also receives feedback signals from the leaf area density vertical profile from the voxelized canopy structure parameter inversion module. When the feedback signal indicates that there are still significant leaf area density peaks below the initial threshold, the stratification threshold is automatically adjusted downwards, triggering the re-execution of this module. This feedback interface realizes closed-loop collaboration between the preceding and following modules, which is one of the key architectural features that distinguishes the system of this invention from existing technologies.

[0045] The energy functional-driven single-tree segmentation module is configured to receive the tree canopy point set output by the multi-scale vertical stratified canopy point extraction module as input, and perform canopy height model generation, variable window local maximum detection, label-controlled watershed pre-segmentation, and graph cut energy optimization as described in step S3 of the above method embodiment. This module includes a canopy height model generation subunit, a treetop detection subunit, a watershed segmentation subunit, and a graph cut optimization subunit. The graph cut optimization subunit is the core component for improving single-tree segmentation accuracy. It globally optimizes and adjusts the watershed pre-segmentation results by constructing a region adjacency graph and minimizing the energy functional containing data items and smoothing terms. The output of this module is a set of segmentation point clouds carrying single-tree label information.

[0046] The voxelized canopy structure parameter inversion module is configured to receive the point clouds of each tree segmentation output by the energy functional-driven tree segmentation module as input, and perform the three-dimensional voxelization processing, leaf area density vertical profile inversion based on the modified Beer-Lambert law, and extraction of conventional structural parameters such as tree height, canopy width, canopy length, and canopy volume as described in step S4 of the above method embodiment. This module is also configured with a feedback output interface to transmit the inverted leaf area density vertical profile to the threshold optimization subunit of the multi-scale vertically layered canopy point extraction module, thereby achieving closed-loop driving of the canopy point extraction process by the parameter inversion results.

[0047] The octree-accelerated 3D reconstruction and visualization module is configured to receive the outputs of the energy functional-driven single-tree segmentation module and the voxelized canopy structure parameter inversion module, and perform octree spatial index construction and multi-resolution mesh reconstruction as described in step S5 of the above method embodiment. This module includes an octree index construction subunit and a Poisson surface reconstruction subunit. The octree index construction subunit adaptively establishes a hierarchical spatial index structure based on the spatial distribution characteristics of the point cloud to support efficient access to large-scale point cloud data and frustum clipping rendering. The Poisson surface reconstruction subunit performs multi-resolution surface reconstruction on the forest stand canopy point cloud based on the octree index, generating an interactively viewable 3D canopy model. The 3D model supports interactive functions such as arbitrary viewpoint roaming, vertical profile analysis, and single-tree parameter querying through the rendering engine, providing visual decision support for forestry management and ecological monitoring.

[0048] The data flow relationships between the five modules constitute the deeply coupled closed-loop architecture of the system of this invention. Specifically, the output of the adaptive terrain-aware point cloud preprocessing module (digital elevation model and vegetation point set) serves as the input of the multi-scale vertical stratified canopy point extraction module; the tree canopy point set output by the multi-scale vertical stratified canopy point extraction module serves as the input of the energy functional-driven single-tree segmentation module; the single-tree segmentation point cloud output by the energy functional-driven single-tree segmentation module serves as the input of the voxelized canopy structure parameter inversion module; the voxelized canopy structure parameter inversion module outputs the inversion results to the octree-accelerated 3D reconstruction and visualization module for 3D model generation, and also feeds back the leaf area density vertical profile information to the multi-scale vertical stratified canopy point extraction module through a feedback interface to drive the dynamic optimization of the stratification threshold. This closed-loop architecture combining forward propagation and backward feedback ensures that the modules are no longer simply cascaded, but rather form a mutually synergistic and mutually optimized organic whole, achieving a synergistic gain effect in the overall system performance. At the hardware implementation level, the aforementioned modules can be deployed on the same high-performance computing platform, achieving efficient data exchange and state synchronization between modules through shared memory or high-speed message queues. Preferably, the energy functional-driven single-tree segmentation module and the octree-accelerated 3D reconstruction and visualization module can utilize the parallel computing capabilities of the graphics processor for accelerated processing, thereby further significantly improving the overall processing efficiency and real-time interactive response speed of the system.

[0049] To verify the effectiveness and superiority of the method of this invention, a systematic comparative experiment was conducted on a real forest lidar dataset. The test area was selected in a subtropical evergreen broad-leaved mixed forest area in southeastern China, covering an area of ​​approximately 50 hectares, with an altitude ranging from 320m to 780m and a terrain slope ranging from 2° to 35°, encompassing various terrain types such as gentle valleys and steep slopes. The forest stand structure is a multi-layered uneven-aged mixed forest, with an average tree height of approximately 18m, a maximum tree height of approximately 32m, and a canopy closure of approximately 0.75 to 0.85. Airborne lidar data was acquired using a small-spot full-waveform laser scanning system flying at an altitude of 800m, with an average point density of approximately 6.5 points / m². 2 It supports up to 5 echo separation records. Ground-based reference data includes total station height records of 206 sample trees, plant area index measurements of 30 20m×20m quadrats, and manually interpreted boundaries of 1500 individual trees.

[0050] Regarding ground point classification accuracy, the slope-adaptive progressive morphological filtering algorithm of this invention achieves an overall classification accuracy of 96.8% across the entire test area. Specifically, the accuracy is 98.1% in gentle slopes less than 5°, 96.5% in moderate slopes of 5° to 15°, and 94.2% in steep slopes greater than 15°. In comparison, the progressive morphological filtering algorithm using a fixed window size of 15m achieves classification accuracies of 97.3%, 93.8%, and 89.5% under the same three terrain conditions. Experimental results show that the slope-adaptive mechanism of this invention provides the most significant improvement in classification accuracy under steep terrain conditions, improving by approximately 4.7 percentage points compared to the fixed window scheme. This is because a fixed large window in steep slope areas easily misclassifies ground abrupt changes as vegetation points.

[0051] Regarding the accuracy of tree height measurement, using the height of 206 sample trees measured by a total station in the field as a reference, the root mean square error (RMSE) of the tree height estimation method of this invention is 0.42m, the coefficient of determination is 0.96, and the mean absolute deviation is 0.33m, meeting the design goal of a tree height measurement accuracy better than 0.5m. Further analysis revealed that for tall trees greater than 20m, the RMSE is 0.38m; for medium-sized trees between 10m and 20m, the RMSE is 0.45m; and for dwarf trees less than 10m, the RMSE is 0.51m. The measurement error for dwarf trees is relatively large, mainly because their canopy point cloud density is low and they are easily affected by the shading effect of surrounding tall trees.

[0052] Regarding the accuracy of individual tree segmentation, using the boundaries of 1500 individual trees interpreted visually as a reference, the correct detection rate of the coupled algorithm of marker-controlled watershed and graph cut energy optimization in this invention is 87.6%, the oversegmentation rate is 5.8%, and the undersegmentation rate is 6.6%. In comparison, the correct detection rate using only the traditional watershed algorithm is 76.3%, with an oversegmentation rate of 14.2%; the correct detection rate using only the local maximum plus region growing scheme is 79.1%, with an undersegmentation rate of 12.5%. Further analysis under different canopy closure conditions shows that the correct detection rate of this invention in high canopy closure areas (canopy closure greater than 0.8) is 84.3%, an improvement of approximately 13.7 percentage points compared to the 70.6% of the traditional watershed scheme, fully demonstrating the advantages of graph cut energy optimization under dense canopy conditions.

[0053] Regarding the accuracy of leaf area density vertical profile inversion, using the leaf area index of 30 quadrats measured by a ground-based plant area index meter as a reference, the correlation coefficient between the leaf area index inverted based on the voxelized Beer-Lambert model of this invention and the measured values ​​is 0.91, and the root mean square error is 0.38m. 2 / m 2The vertical profile of leaf area density clearly reflects the multi-layered structural characteristics within the canopy, with the peak position consistent with the actual height distribution of the main canopy and sub-canopy layers. Preferably, under the condition of setting the voxel size to 1m×1m×1m, the average deviation of the retrieved leaf area density profile within the main canopy height range is 0.12m. 2 / m 3 The overall accuracy meets the accuracy requirements of canopy structure data for forest stock volume estimation and carbon sink monitoring.

[0054] Regarding the closed-loop feedback optimization effect, after dynamic adjustment of the stratification threshold driven by leaf area density profile feedback, the complete extraction rate of canopy points increased from the initial 91.2% to 96.5%, an improvement of approximately 5.3 percentage points, mainly improving the problem of missing point clouds at the bottom of the canopy of sub-canopy trees. In the 1500 reference trees in the 50-ha test area, after closed-loop feedback optimization, the canopy bottom point clouds of approximately 79 sub-canopy trees that were originally incorrectly removed were recovered, resulting in an average improvement of approximately 15% in the measurement accuracy of the canopy length parameter of these trees. In terms of 3D reconstruction efficiency, the total processing time for multi-resolution Poisson surface reconstruction based on octree spatial indexing in the 50-ha test area was approximately 25 minutes (test platform configuration: Intel Core i7-12700 processor, 32GB memory, and NVIDIA RTX 3060 graphics card), supporting smooth real-time roaming browsing of the 3D model with a frame rate maintained above 30 frames per second.

[0055] In summary, this invention achieves a high-precision, high-efficiency automated reconstruction process from raw point clouds to 3D canopy models through a five-step deeply coupled closed-loop architecture. The deep coupling and closed-loop feedback mechanism between each step make the overall system performance significantly better than the serial schemes in existing technologies where each step operates independently. It can provide refined 3D canopy structure data support for applications such as forest stock volume estimation, carbon sink accounting, ecosystem function evaluation, and refined management of forestry resources, and has broad practical application prospects and extremely important scientific research value.

[0056] The embodiments of the present invention are not limited to the specific embodiments described above. Those skilled in the art can make various equivalent changes or substitutions based on the technical solutions of the present invention, and all such changes or substitutions should be included within the protection scope of the present invention.

Claims

1. A method for reconstructing three-dimensional structure of forest canopy from lidar point cloud images, characterized in that, Includes the following steps: Adaptive terrain-aware point cloud preprocessing steps: acquire discrete point cloud data of the target forest area scanned by airborne lidar, perform noise filtering and coordinate unification on the point cloud data, and use slope adaptive progressive morphological filtering algorithm to separate ground points and vegetation points in the processed point cloud. The size of the filtering window is dynamically adjusted according to the local terrain slope. A digital elevation model is generated by interpolating ground points through irregular triangular mesh. Multi-scale vertical stratified canopy point extraction steps: Based on the digital elevation model, the vegetation points are height normalized to obtain a normalized point cloud, a vertical height distribution density histogram of the normalized point cloud is constructed, and the stratification thresholds of the shrub layer and the tree layer are adaptively determined according to the density valley values ​​in the density histogram. Shrub layer point clouds below the stratification thresholds are removed to obtain the tree layer canopy point set. The steps of single-tree segmentation driven by energy functional are as follows: a canopy height model is generated based on the canopy point set; a variable window local maximum detection is used to extract candidate tree top points as labeling seed points; label-controlled watershed segmentation is performed on the canopy height model using the labeling seed points to obtain the initial segmentation region; a region adjacency graph is constructed with spatial distance and height difference as weights; and graph cut optimization is performed on the initial segmentation region by minimizing the energy functional to obtain the single-tree segmentation result. Voxelized canopy structure parameter inversion steps: The point cloud of each tree segmentation result is processed into three-dimensional voxels. The leaf area density is calculated based on the laser pulse penetration probability model in each horizontal layer voxel array to obtain the vertical profile of leaf area density. The tree height, crown width, crown length and canopy volume are extracted by combining the spatial envelope of the tree point cloud. The vertical profile of leaf area density is fed back to the multi-scale vertical layered canopy point extraction step to optimize the layering threshold. Octree accelerated 3D reconstruction and visualization steps: an octree spatial index is established for the individual tree segmentation results and canopy structure parameters. Based on the octree spatial index, multi-resolution mesh reconstruction is performed on the forest stand canopy point cloud to generate a canopy 3D model. The canopy 3D model supports arbitrary viewpoint roaming and vertical profile analysis.

2. The method for reconstructing forest canopy three-dimensional structure lidar point cloud images according to claim 1, characterized in that, The point density of the point cloud data is not less than 4 points / m. 2 The echo recording supports multiple echo separations. The noise filtering adopts an outlier removal method based on elevation statistics, and the removal threshold is set to three times the standard deviation of the local elevation mean.

3. The method for reconstructing forest canopy three-dimensional structure lidar point cloud images according to claim 1, characterized in that, In the slope-adaptive progressive morphological filtering algorithm, the dynamic adjustment rule for the morphological filtering window size is as follows: when the local terrain slope is less than 5°, the window size is set to 20m; when the local terrain slope is between 5° and 15°, the window size linearly shrinks to 10m; and when the local terrain slope is greater than 15°, the window size is set to 5m. The grid resolution of the irregular triangular mesh interpolation is set to 0.5m to 1m.

4. The method for reconstructing three-dimensional structure of forest canopy from lidar point cloud images according to claim 1, characterized in that, The height layer spacing of the density histogram is set to 0.5m, the density valley value is determined by extracting the minimum point after performing Gaussian smoothing on the density histogram, and the initial value of the layer threshold is set to 2m to 4m.

5. The method for reconstructing forest canopy three-dimensional structure lidar point cloud images according to claim 1, characterized in that, In the variable window local maximum detection algorithm, the search window radius is dynamically determined by a linear regression model based on the height value at the corresponding position in the canopy height model. The linear regression model is expressed as the search window radius being equal to the product of the regression coefficient and the height value at the corresponding position plus the regression intercept. The regression coefficient ranges from 0.04 to 0.08, and the regression intercept ranges from 0.5m to 1.5m.

6. The method for reconstructing forest canopy three-dimensional structure lidar point cloud images according to claim 1, characterized in that, The energy functional includes a data term and a smoothing term. The data term measures the spatial distance cost from each point cloud to its respective tree center, and the smoothing term measures the penalty for label inconsistency between adjacent segmented regions. The final tree assignment label of each point cloud is determined by finding the global minimum of the energy functional through a graph cut algorithm.

7. The method for reconstructing three-dimensional structure of forest canopy point cloud images using lidar according to claim 1, characterized in that, The voxel size for the three-dimensional voxelization process is set to 0.5m×0.5m×0.5m to 1m×1m×1m. The laser pulse penetration probability model is based on the modified Beer-Lambert law, and the leaf area density is calculated by statistically analyzing the number of laser pulses intercepted and penetrated in each voxel.

8. The method for reconstructing forest canopy three-dimensional structure lidar point cloud images according to claim 1, characterized in that, The method for optimizing the stratification threshold by feeding back the leaf area density vertical profile to the multi-scale vertical stratification canopy point extraction step is as follows: when the leaf area density vertical profile still has a peak value greater than the set density threshold below the initial stratification threshold, the stratification threshold is adjusted downward to the density valley value position below the peak value.

9. The method for reconstructing forest canopy three-dimensional structure lidar point cloud images according to claim 1, characterized in that, The maximum recursion depth of the octree spatial index is set to 8 to 12 layers, and the multi-resolution mesh reconstruction adopts the Poisson surface reconstruction algorithm based on normal vector estimation, with a reconstruction depth set to 8 to 10.

10. A forest canopy three-dimensional structure lidar point cloud image reconstruction system, used to implement the forest canopy three-dimensional structure lidar point cloud image reconstruction method according to any one of claims 1-9, characterized in that, include: The adaptive terrain-aware point cloud preprocessing module is configured to acquire discrete point cloud data of the target forest area scanned by airborne lidar, perform noise filtering and coordinate system unification processing on the point cloud data, use the slope adaptive progressive morphological filtering algorithm to separate ground points and vegetation points in the processed point cloud data, and use the separated ground points to generate a digital elevation model through irregular triangular network interpolation. The multi-scale vertical stratified canopy point extraction module is configured to perform height normalization processing on the vegetation points based on the digital elevation model to obtain a normalized point cloud, construct a vertical height distribution density histogram and adaptively determine the stratification threshold according to the density valley value, and obtain the tree canopy point set after removing the shrub layer point cloud. The energy functional-driven single-tree segmentation module is configured to generate a canopy height model based on the canopy point set, extract labeled seed points using variable window local maximum detection and perform labeled control watershed segmentation to obtain an initial segmentation region, and perform graph cut optimization on the initial segmentation region by minimizing the energy functional to obtain the single-tree segmentation result; The voxelized canopy structure parameter inversion module is configured to perform three-dimensional voxelization processing on each tree point cloud, calculate the leaf area density based on the laser pulse penetration probability model to obtain the leaf area density vertical profile, extract tree height, canopy width, canopy length and canopy volume parameters, and feed the leaf area density vertical profile back to the multi-scale vertical layered canopy point extraction module to optimize the layering threshold. The Octree accelerated 3D reconstruction and visualization module is configured to create an octree spatial index and perform multi-resolution mesh reconstruction on the canopy point cloud to generate a canopy 3D model that supports arbitrary viewpoint roaming and vertical profile analysis.

Images