Construction method of digital canopy height model of artificial forest based on high-precision point cloud data
By combining high-precision point cloud data and RGB imagery, and adaptively adjusting the grid size and stiffness, the CHM accuracy problem caused by the complexity of artificial forest terrain was solved, and a high-precision and robust canopy height model was constructed.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- GUIZHOU NORMAL UNIVERSITY
- Filing Date
- 2026-04-27
- Publication Date
- 2026-06-02
AI Technical Summary
Existing technologies, when constructing digital canopy height models for plantations, struggle to adapt to the spatial heterogeneity of plantation terrain complexity and canopy disturbance intensity using fixed-stiffness fabric simulation filtering algorithms. This results in large ground point classification errors, affecting the accuracy and reliability of DEM and CHM.
By acquiring high-precision point cloud data and RGB visible light images, a spatial mapping is established, a grid is divided, local RGB image features are analyzed, core candidate points are screened, the grid size is adjusted, an adaptive stiffness CSF algorithm is constructed, point cloud data is classified, and a high-precision CHM model is generated.
It improves the accuracy and robustness of the plantation CHM model, enabling it to accurately reflect canopy height distribution in complex environments, reduce the influence of noise points, and enhance the model's identification ability and accuracy.
Smart Images
Figure CN122134944A_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of image processing technology, specifically to a method for constructing a digital canopy height model of plantations based on high-precision point cloud data. Background Technology
[0002] Plantations are important ecosystems, widely used in timber production, carbon sequestration assessment, biomass estimation, and forest structure regulation. They are typically characterized by monoculture, regular rows and columns, and controllable density. Canopy height is a key parameter for measuring the growth and spatial structure of plantations. Digital canopy height models (CHMs) describe the horizontal and vertical distribution of tree canopies and are an important data source for forest structure analysis. With the development of UAV lidar and high-resolution photography technologies, high-precision 3D point cloud data has become the primary data source for constructing CHMs.
[0003] CHM (Content Surface Model) is typically obtained by subtracting the Digital Surface Model (DSM) from the Digital Elevation Model (DEM), where the DEM is generated by ground point cloud interpolation, and its accuracy has a decisive impact on the CHM results. The Fabric Simulation Filter (CSF) algorithm is widely used for ground point cloud extraction and DEM construction due to its intuitive principle, simple implementation, and relatively easy parameter setting. However, this method uses a globally uniform fabric stiffness parameter, essentially assuming that surface undulation features and non-ground target disturbances are spatially consistent. This assumption is often difficult to satisfy in plantation forest environments.
[0004] Plantations are often distributed in areas with complex terrain, such as slopes and valleys, where surface slopes and micro-topographic features vary greatly. Fixed-stiffness fabric models struggle to accurately identify ground points on steep slopes or in depressions, easily leading to missed detections; in gentler areas, they may misclassify low vegetation as ground points. Furthermore, the dense canopy of plantations results in fewer and unevenly distributed ground points in lidar echoes, further exacerbating ground point classification errors and causing DEM distortion, thus affecting the accuracy of CHM and the reliability of subsequent applications. Summary of the Invention
[0005] In view of the above, it is necessary to provide a method for constructing a digital canopy height model of plantations based on high-precision point cloud data to solve the above problems.
[0006] One embodiment of this application provides a method for constructing a digital canopy height model of plantations based on high-precision point cloud data. The method includes: Point cloud data and RGB visible light images of the target plantation area were collected; spatial mapping from point cloud data to RGB visible light images was established using collinearity equations, and the local RGB image of each point cloud data was determined. The point cloud data is initially divided into grids of a preset size to determine the target analysis points. The edge distribution features and regional color and texture features in the local RGB image of each target analysis point are analyzed to obtain the canopy region feature coefficients. The target analysis points are screened based on all the canopy region feature coefficients obtained within the grid to obtain core candidate points. The elevation data of the core candidate points corresponding to the grid are analyzed to determine the elevation undulation coefficient. The initial grid size is adaptively adjusted based on the canopy region feature coefficients and the elevation undulation coefficients to extract the seed terrain skeleton. By utilizing the shape features of each grid after size adjustment and its corresponding seed terrain skeleton, and combining the canopy region feature coefficients corresponding to all target analysis points in each grid, the adaptive stiffness of the CSF algorithm is constructed, and the point cloud data is classified by the CSF algorithm. A CHM model of planted forests was constructed based on the classified point cloud data.
[0007] Specifically, determining the local RGB image of each point cloud data point is as follows: Based on the spatial location of the point cloud data and the camera position, combined with the laser divergence angle, the physical surface diameter of each point cloud data is calculated; based on the point cloud data and camera parameters, the sampling distance of the ground corresponding to each point cloud data is calculated; the product of the negative correlation mapping result of the sampling distance of each point cloud data and the physical surface diameter is used as the diameter of the local RGB image of each point cloud data. Using the projection coordinates of each point cloud data in the RGB visible light image as the center, draw a circle based on the diameter, and use all the pixels inside the circle as the local RGB image corresponding to each point cloud data.
[0008] Specifically, determining the target analysis point involves: Three-dimensional coordinates In Point cloud data falling within each grid area are used as target analysis points; if in There are multiple point cloud data points at the location, The point cloud data with the smallest value is used as the target analysis point at that location.
[0009] Specifically, the canopy region characteristic coefficients are obtained as follows: The proportion of edge pixels in the local RGB image of each target analysis point is used as the edge density, denoted as . ; Extract the principal axis direction of the edge lines in the local RGB image of each target analysis point, and calculate the information entropy of all principal axis directions, denoted as . ; Through formula Calculate the vegetation texture features of the local RGB image corresponding to each target analysis point. Where n represents the number of species along the principal axis. Represents the logarithmic function with base 2; It is a preset constant; Calculate the mean difference between the green over-green index and the red over-red index of all target analysis points in the local RGB image of each target analysis point, and then normalize the result. ; Obtain the mean value of the S-channel pixel value of all pixels in the local RGB image corresponding to each target analysis point in HSV space, and denot it after normalization. ; Will , , The mean value is used as the canopy region feature coefficient of the local RGB image of each target analysis point.
[0010] The specific process for obtaining the core candidate points is as follows: Sort the canopy region feature coefficients of all target analysis points corresponding to each grid region in descending order, select the smaller of the two elements with the largest difference between adjacent elements, and select the point cloud data corresponding to all elements sorted after this element as the core candidate points for each grid.
[0011] Specifically, determining the elevation fluctuation coefficient involves: Based on all core candidate points corresponding to each grid, a fitting plane is obtained, and the angle between the fitting plane and the plane where the grid is located is obtained. Combining the dispersion of the elevation values of all core candidate points, the elevation fluctuation coefficient of each grid is obtained.
[0012] The specific process of adaptively adjusting the initial mesh size is as follows: Based on the elevation undulation coefficient of each grid and the average canopy region characteristic coefficient of all core candidate points corresponding to the grid, the decomposition decision coefficient of each grid is determined; wherein, the decomposition decision coefficient is positively correlated with the elevation undulation coefficient and negatively correlated with the average canopy region characteristic coefficient. Set a decision threshold. When the decomposition decision coefficient of a grid is greater than the decision threshold, perform quadtree splitting on the grid; otherwise, stop splitting.
[0013] Specifically, the extraction of the seed terrain skeleton includes: Among all core candidate points corresponding to the grid after the splitting stops, the core candidate point with the smallest sum of the normalized elevation value and the canopy region feature coefficient is selected as the seed point of the corresponding grid. The difference between the natural number 1 and the canopy region feature coefficient of each target analysis point in each grid after the splitting stops is calculated. The mean of all the differences is used as the constraint weight of the corresponding seed point. The seed terrain skeleton is obtained by interpolation.
[0014] Specifically, the process of constructing the adaptive stiffness of the CSF algorithm and classifying point cloud data using the CSF algorithm is as follows: Based on the interpolated continuous seed terrain skeleton surface, for each grid after the splitting stops, the standard deviation of all curvature values within that grid is calculated; the mean of the minimum angle difference between the average aspect of that grid and the average aspect of all adjacent grids is calculated; the average slope within that grid is calculated; and the standard deviation, the mean of the minimum angle difference, and the average slope are all globally normalized and then averaged to obtain the comprehensive terrain complexity index of the corresponding grid, denoted as . ; By interpolating the canopy region characteristic coefficients of all target analysis points within each grid where splitting stops, the canopy disturbance intensity field is obtained, denoted as... ; The specific formula for the adaptive stiffness K of the spatial point corresponding to the seed terrain skeleton surface within each grid of the cloth node in the CSF algorithm is as follows: Where norm() represents the normalization function, This is the difference between the preset maximum stiffness and the preset minimum stiffness. Indicates the preset minimum stiffness; Represents a preset constant; Using the fabric surface after CSF algorithm iteration as a reference surface, the vertical distance between the fabric nodes and the corresponding original point cloud data is calculated. Point cloud data with a vertical distance greater than a set threshold are classified as non-ground point cloud data; point cloud data with a vertical distance less than or equal to the set threshold are classified as ground point cloud data.
[0015] The construction of the plantation forest CHM model based on the classified point cloud data includes: Interpolation calculations are performed on ground point cloud data to generate a DEM; Non-ground point cloud data is divided based on DEM raster to obtain the top canopy point cloud, and the top canopy point cloud is interpolated to generate DSM; The DSM and DEM are subtracted pixel by pixel to generate the original CHM; the negative height values in the original CHM are uniformly set to 0; the average height and maximum height of the canopy are statistically analyzed, and abnormal high values that exceed the preset quantile height by a preset multiple are replaced with the height mean of their preset neighborhood units to obtain the corrected CHM. After filtering, the plantation forest CHM model is obtained.
[0016] This application has at least the following beneficial effects: This application collects point cloud data and RGB visible light images of the target plantation area and establishes a spatial mapping through collinearity equations. This allows for the assignment of corresponding RGB image information to each point cloud data point, aiding in the accurate identification of the spatial structure of the canopy area. Preliminary grid-based subdivision of the point cloud data refines the data processing area, making the point cloud data within each grid more concentrated and uniform, ensuring that subsequent analysis focuses on local features and reduces errors from global processing. By analyzing the edge distribution and color texture features of the local RGB images of each target analysis point, canopy area feature coefficients are obtained, helping to distinguish canopy area feature coefficients, enhancing the model's identification ability in complex forest environments, and improving the accuracy of canopy capture during CHM construction. By screening core candidate points, more representative points can be prioritized. By identifying key ground points and important areas, and reducing the impact of noise points, the accuracy of point cloud classification is improved. Analyzing elevation data within the grid and determining the elevation undulation coefficient helps adjust the grid size to adapt to different terrain features, ensuring high-precision point cloud classification under various terrain conditions and thus improving the accuracy of the DEM. Combining the canopy region feature coefficients within the grid with the adjusted seed terrain skeleton, the adaptive stiffness optimization CSF algorithm makes the data distribution model more flexible and accurate, adapting to complex terrain changes and further improving the accuracy of ground point classification. Finally, a CHM model is constructed based on the classified point cloud data, accurately reflecting the spatial distribution of plantation canopy height, effectively improving the accuracy of plantation CHM construction, and making the final canopy model more robust and accurate in complex environments. Attached Figure Description
[0017] Figure 1 A flowchart illustrating the method for constructing a digital canopy height model of planted forests based on high-precision point cloud data, as provided in this application. Detailed Implementation
[0018] In the description of the embodiments in this application, the words "exemplary," "or," and "for example" are used to indicate examples, illustrations, or descriptions. Any embodiment or design scheme described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or design schemes. Specifically, the use of the words "exemplary," "or," and "for example" is intended to present the relevant concepts in a specific manner.
[0019] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The terminology used in this application's specification is for the purpose of describing particular embodiments only and is not intended to be limiting of the application.
[0020] It should also be noted that the terms "first" and "second" in this application and its accompanying drawings are used to distinguish similar objects, rather than to describe a specific order or sequence. The methods disclosed in the embodiments of this application or the methods shown in the flowcharts include one or more steps for implementing the method. Without departing from the scope of protection of this application, the execution order of multiple steps can be interchanged, and some steps can also be deleted.
[0021] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application pertains.
[0022] This application proposes a method for constructing a digital canopy height model of plantations based on high-precision point cloud data, which is applied to the field of image processing technology. (See attached document.) Figure 1 The method includes the following steps: S1: Collect point cloud data and RGB visible light images of the target artificial forest area; establish a spatial mapping from point cloud data to RGB visible light images through collinearity equations, and determine the local RGB image of each point cloud data.
[0023] In plantation forest scenarios, the quality of high-precision point cloud data directly determines the reliability of all subsequent processing steps. Since plantations are typically located in areas with complex terrain, the drone's lidar flight is affected by factors such as terrain undulations and airflow disturbances, leading to systematic errors in the point cloud data. Simultaneously, lidar generates a large number of noise points when penetrating high-canopy layers. If this noise is not preprocessed, it will be misclassified in subsequent ground point classification, further amplifying the classification error of the CSF algorithm. Therefore, the first step is to collect plantation forest data.
[0024] A drone platform equipped with a lidar, high-resolution RGB camera, and GNSS / IMU integrated navigation module was used to collaboratively collect data on the target plantation area. By unifying the time reference and extrinsic parameter calibration, strict spatiotemporal alignment of LiDAR point clouds, RGB visible light imagery, and pose information was achieved, providing a multi-source data foundation for subsequent point cloud semantic enhancement and terrain analysis. Flight path planning employed a high forward and lateral overlap rate, and a multi-echo mode was enabled to record multiple reflections of the laser pulse during its penetration of the canopy. During the drone's flight scan, 3D point cloud data, high-resolution RGB visible light imagery, and GNSS / IMU attitude and position information were simultaneously acquired. A timestamp alignment mechanism ensured that each frame of RGB visible light imagery could be correlated with the corresponding point cloud data. The GNSS, IMU, and LiDAR point clouds were jointly calculated, and a seven-parameter coordinate transformation method was used to project all point cloud data onto the Universal Transverse Mercator (UTM) projection coordinate system under the WGS-84 datum, completing spatial registration between the sensors.
[0025] Point cloud data is preprocessed using point cloud quality control methods such as statistical filtering and radius filtering to obtain denoised, coordinate-unified, and structurally complete LAS format 3D point cloud data.
[0026] The detection beam of a lidar is not an ideal line, but a cone with a divergence angle. Therefore, each point in the point cloud actually represents a circular light spot in physical space; at the same location, the physical size represented by a pixel on the camera increases with the distance between the UAV and the target. Therefore, it is necessary to analyze the size of the lidar light spots in different point clouds and the corresponding ground resolution of the images to establish the correspondence between point cloud data and RGB images.
[0027] First, using the pose data and camera intrinsic parameters obtained during the data acquisition phase, a collinearity equation is established, and each point cloud data is projected onto the corresponding RGB image plane to obtain the corresponding image pixels and their coordinates. Based on the instantaneous distance between the UAV and the target (calculated by the spatial distance between the point cloud coordinates and the camera center) and the divergence angle of the laser pulse, the diameter of the point cloud data on the physical surface of the target object is calculated. .
[0028] By combining the instantaneous distance between the drone and the target with the RGB camera lens parameters (specifically including pixel size and lens focal length), the sampling distance of the ground corresponding to this point cloud data is calculated. This value represents the actual size of the object corresponding to each pixel in the image (e.g., ...). When the value is 5, the actual object size corresponding to that pixel is 5×5 (the unit of this value is determined by the instantaneous distance and the lens parameter unit). The calculation process of the above parameters is all existing technology and will not be elaborated on here.
[0029] For a given point cloud data, the diameter of the pixel range it represents in the corresponding image should be: ,in, As a preset constant, in this embodiment Its function is to prevent the denominator from being 0; This represents the negative correlation mapping result of the sampling distance in the point cloud data.
[0030] Using the image pixel coordinates corresponding to the current point cloud data as the center, a preset neighborhood expansion coefficient (10 in this embodiment) is introduced, and the... Multiply by this coefficient to obtain the final local RGB image diameter, ensuring that sufficient context pixels are included for edge detection. All pixels within the circle are used as the local RGB image corresponding to the current point cloud data.
[0031] S2: Initially divide the point cloud data into grids of preset size to determine target analysis points; analyze the edge distribution features and regional color texture features in the local RGB image of each target analysis point to obtain canopy region feature coefficients; filter the target analysis points based on all canopy region feature coefficients obtained within the grid to obtain core candidate points; analyze the elevation data of the core candidate points corresponding to the grid to determine the elevation undulation coefficient; adaptively adjust the initial grid size based on the canopy region feature coefficients and the elevation undulation coefficients to extract the seed terrain skeleton.
[0032] Fixed-stiffness fabric models lack adaptability under conditions of complex plantation terrain and spatial heterogeneity of canopy disturbance intensity, leading to systematic biases in ground point classification across different terrain units. Traditional CSF algorithms assume globally fixed fabric stiffness parameters, implicitly assuming uniform surface morphology and consistent canopy disturbance intensity. However, in plantation scenarios, this assumption can affect the accuracy of DEM models because the terrain and canopy features of plantations often exhibit significant spatial variations.
[0033] To improve the accuracy of the CHM layer in plantations, the traditional CSF algorithm requires parameter separation to distinguish the ground. However, the optimal parameters depend on the ground morphology. Therefore, before constructing the DEM, ground points in the point cloud data need to be roughly extracted. By using the lower envelope operator, while maintaining topological coherence, interference from high-level canopies is filtered out, and extremely sparse but relatively close to the ground surface seed points are extracted. Although these seed points are insufficient to construct the DEM, they can express the slope variation trend of the terrain.
[0034] In densely wooded areas, the ground points of the drone's lidar are extremely sparse. In the corresponding image area, the main features are dense canopy areas with high greenness and high texture. At this time, the laser penetration rate is extremely low. Therefore, it is necessary to increase the grid size to improve the probability of capturing ground echo points within the grid and avoid mistaking the canopy area for the ground. At the same time, the brown understory ground area often shows low greenness and low saturation.
[0035] The target artificial forest area is initially divided into grids of a preset size. In this embodiment, the initial grid size is 10m × 10m. The three-dimensional coordinates of each point cloud data are compared with the grid boundary to determine its grid affiliation, and the point cloud set within each grid is recorded. In this embodiment, the grid coordinate range is set to left-closed and right-open, and bottom-closed and top-open. Next, all point cloud data within a specific grid are acquired, and the point cloud data with the lowest height within each point on the grid plane is selected. It should be noted that in the point cloud space, the grid is defined as a subspace within the x-coordinate range and the y-coordinate range, while each point within each grid is a point with fixed x and y coordinates in that subspace, and the z-coordinate is variable. The three-dimensional coordinates are... In Point cloud data falling within each grid area are used as target analysis points; if in There are multiple point cloud data at the location. The point cloud data with the smallest z value is taken as the target analysis point at that location to obtain the target analysis point of the current grid. The two-dimensional local RGB image corresponding to the target analysis point is extracted. The purpose is to evaluate the degree of canopy closure in the vertical space above the ground point and the intensity of interference to laser penetration.
[0036] Calculate the mean difference between the green over-green index and the red over-red index of all target analysis points in the local RGB image corresponding to each target analysis point, and normalize it using the sigmoid function. This result is denoted as... The larger this value, the greater the probability that the current image locally represents a vegetation area. Canny edge detection is performed on the local RGB image of the target analysis point to obtain several edge lines in the region. The ratio of the number of pixels on all edge lines to the number of pixels in the local RGB image of the target analysis point is calculated and used as the edge density of the local RGB image of the target analysis point, denoted as . The PCA algorithm is used to extract the principal axis direction of each edge line, and the information entropy of the principal axis direction in the local RGB image of the target analysis point is calculated, denoted as . When the number of vegetation texture features extracted along the main axis is less than or equal to 1, the vegetation texture feature is directly set. If the value is 0, skip the subsequent formula calculations; otherwise, calculate the vegetation texture features of the local RGB image corresponding to each target analysis point. The specific formula is as follows: Where n represents the number of species along the principal axis. Represents the logarithmic function with base 2; As a preset constant, in this embodiment Its function is to prevent the denominator from being 0.
[0037] It should be noted that only through It cannot be directly determined whether it is a vegetated area. In plantation areas, some man-made structures may also have similar color characteristics. Further analysis is needed, combining edge information within the image: plantation vegetation often has more texture, and the direction of the texture is more chaotic compared to man-made structures. Therefore, through calculation... and The density and orientation disorder of the texture were evaluated separately, and then... Will The value of is set between (0, 1). The larger the value, the more vegetation texture features the corresponding region of the local RGB image of each target analysis point possesses.
[0038] Obtain the pixel value in HSV space for each pixel in the local RGB image corresponding to each target analysis point; calculate the mean value of the S-channel pixel value of the current local image region, and normalize it, denoted as . ; Calculate the canopy region feature coefficients of the current local image. The specific formula is as follows: The larger this value, the greater the probability that the current target analysis point represents the canopy region in the image.
[0039] For the current grid region, the canopy region characteristic coefficients of all target analysis points are calculated. Sort the elements in descending order, calculate the difference between adjacent elements, select the smaller of the two elements with the largest difference, and then select... Point cloud data corresponding to all elements less than or equal to the given element are used as core candidate points for the current grid. It should be noted that if the maximum difference between adjacent elements in the sequence is less than a preset minimum value, then all point cloud data of the top preset proportion (30% in this embodiment) of the canopy region feature coefficients are directly selected as core candidate points.
[0040] Meanwhile, in areas with fragmented terrain or steep slopes, such as the stepped land preparation commonly seen in plantations, the surface elevation changes drastically over short distances. In such cases, it is necessary to reduce the mesh size to prevent large-size meshes from smoothing the terrain by cutting off curves, which could lead to the loss of key terrain features. A fitting plane for the core candidate points within each mesh is constructed using the least squares method, and the angle between this fitting plane and the XOY plane is obtained. Through this included angle and The ratio is normalized and denoted as . Calculate the standard deviation of the elevation values of all core candidate points, denoted as . to After multiplying, the product is normalized and used as the elevation undulation coefficient of the current grid. The normalization method used is the maximum-minimum normalization method.
[0041] It should be understood that, This indicates the overall tilt angle of the core candidate points within the current grid. This indicates the degree of unevenness of the core candidate points within the current grid. The combination of the two indicates the overall and local undulation of the plane represented by the core candidate points within the grid. The larger this value is, the more drastic the change in the elevation of the plane within the current grid.
[0042] Based on the elevation undulation coefficient of each grid Mean of canopy region feature coefficients of all core candidate points corresponding to the grid Determine the decomposition decision coefficient for each grid. The specific formula is as follows: norm() represents the normalization function. As a preset constant, in this embodiment Its function is to prevent the denominator from being 0; a decision threshold N=0.6 is set. When the decomposition decision coefficient of the grid is greater than the decision threshold, the grid is split into a quadtree, otherwise the splitting is stopped.
[0043] After the grid stops splitting, select the normalized value of the elevation value within the grid and... The core candidate point with the smallest sum is used as the seed point for the current grid; the seed points for all grids are obtained, and all point cloud data within each grid are processed. The mean value is used as the constraint weight, and the seed terrain skeleton is obtained by TPS interpolation.
[0044] S3: Utilize the shape features of each grid after size adjustment and its corresponding seed terrain skeleton, and combine the canopy region feature coefficients corresponding to all target analysis points in each grid to construct the adaptive stiffness of the CSF algorithm, and classify the point cloud data through the CSF algorithm.
[0045] The core of the CSF algorithm is to simulate a cloth covering a point cloud from the air. The cloth adheres to the ground surface under gravity; the greater the stiffness, the less likely the cloth is to bend (suitable for flat terrain, avoiding embedding in low vegetation); the smaller the stiffness, the easier it is for the cloth to conform to terrain undulations (suitable for broken, steep slopes, capturing terrain details). Traditional CSF uses a globally fixed stiffness, which cannot adapt to heterogeneous terrains such as flat areas, steep, broken slopes, and densely canopied areas in plantations. By extracting grid-level terrain undulations, canopy interference, slope, and other features based on a seed terrain skeleton, a dynamic mapping relationship between terrain features and stiffness parameters can be constructed, enabling adaptive stiffness adjustment and precise separation of ground points and canopy points.
[0046] Based on the interpolated continuous seed terrain skeleton surface, for each grid after the splitting stops, the standard deviation of all curvature values within that grid is calculated; the mean of the minimum angle difference between the average aspect of that grid and the average aspect of all adjacent grids is calculated; the average slope within that grid is calculated; and the standard deviation, the mean, and the average slope are all globally normalized and then averaged to obtain the comprehensive terrain complexity index for the current grid. The higher the value, the more complex and drastic the terrain changes; the normalization method uses the maximum-minimum normalization method.
[0047] All target analysis points within the grid As a canopy disturbance attribute, a continuous canopy disturbance intensity field is generated through interpolation. Therefore, for a certain point on the current mesh, the specific formula for its adaptive stiffness K is: Where, norm() represents the normalization function, specifically a robust truncation normalization function based on a preset quantile interval (5% to 95% in this embodiment), used to eliminate the influence of extreme outliers generated when the denominator approaches the minimum value on the overall mapping interval. This is the difference between the preset maximum stiffness and the preset minimum stiffness. Indicates the preset minimum stiffness. The value range is 1~10 (low stiffness, suitable for steep slopes). This represents the maximum stiffness, with a value range of 50~200 (high stiffness, suitable for gentle regions). As a preset constant, in this embodiment Its function is to prevent the denominator from being zero; in this embodiment, both the maximum stiffness and the minimum stiffness are preset engineering standard values. The value is 5. The value is 100.
[0048] It should be understood that, The larger the value, the steeper the terrain, and the more priority should be given to fitting it. The larger the value, the denser the vegetation, and the more necessary it is to prevent fabric penetration. By normalizing this ratio and linearly integrating it with the stiffness setting range, K is used as the adaptive stiffness at a point in the current mesh.
[0049] The point cloud is inverted and flipped, and the cloth falls from above (below the original ground). In the mass-spring model iteration, the displacement of each cloth node is constrained by the cloth stiffness parameter corresponding to its position. If the cloth node falls below the ground point, its height is fixed at the ground point position. After several iterations, the final shape of the cloth is an accurate fit to the complex surface.
[0050] Using the fabric surface after CSF algorithm iteration as a reference surface, the point cloud data is divided into ground point cloud and non-ground point cloud based on the vertical distance between the fabric node and the corresponding original point cloud (the threshold is set to 0.1m in this embodiment). Point cloud data with a vertical distance greater than the set threshold is non-ground point cloud data, and point cloud data with a vertical distance less than or equal to the set threshold is ground point cloud data.
[0051] S4: Construct a CHM model of artificial forests based on the classified point cloud data.
[0052] Based on the acquired ground point cloud data, ordinary kriging interpolation is used. Assuming that the spatial variation of terrain elevation has randomness and structure, the spatial correlation of ground points is fitted by a variogram function. Interpolation calculation is performed on the ground point cloud to generate a raster DEM, which is output in TIFF format and contains the absolute ground elevation value of each raster.
[0053] Based on the acquired non-ground point cloud data, grid cells are divided to be fully aligned with the DEM grid. Within each cell, the point with the highest elevation is selected as the top point cloud of the canopy. For grid cells without a top point of the canopy, inverse distance weighted interpolation is used to supplement the top point elevation of neighboring cells. Interpolation calculations are performed on the top point cloud of the canopy to generate a grid DSM, which is output in TIFF format, containing the absolute elevation value of the top of the canopy for each grid.
[0054] According to the core definition of CHM, the original CHM is obtained by subtracting DSM from DEM pixel by pixel, and its physical meaning is the vertical height from the top of the canopy to the ground.
[0055] Using a raster calculation tool, the DSM and DEM are subtracted pixel by pixel to generate the original CHM. Negative height values in the original CHM are uniformly set to 0. The average and maximum canopy heights are statistically analyzed. Abnormal height values exceeding 1.5 times the 99th percentile height are replaced with the average of non-abnormal height values within their 3×3 neighborhood. If all values within that neighborhood are considered abnormally high, the neighborhood is progressively expanded until non-abnormal height values are included in the calculation. It should be noted that negative height values mainly originate from DEM elevations being too high or DSM elevations being too low, while abnormally high values mainly originate from noise points in the DSM. Median filtering is used to smooth the corrected CHM, eliminating jagged edges on the raster edges. A TIFF format plantation CHM model is output, containing the canopy height value for each raster.
[0056] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to embodiments of this application. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of code containing one or more executable instructions for implementing a specified logical function. In some alternative implementations, the functions marked in the blocks may occur in a different order than that shown in the drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. In the descriptions corresponding to the flowcharts and block diagrams in the accompanying drawings, the operations or steps corresponding to different blocks may also occur in a different order than disclosed in the description; sometimes there is no specific order between different operations or steps. For example, two consecutive operations or steps may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. Each block in a block diagram and / or flowchart, and combinations of blocks in a block diagram and / or flowchart, can be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions.
[0057] The above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.
Claims
1. A method for constructing a digital canopy height model of plantations based on high-precision point cloud data, characterized in that, The method includes the following steps: Point cloud data and RGB visible light images of the target plantation area were collected; spatial mapping from point cloud data to RGB visible light images was established using collinearity equations, and the local RGB image of each point cloud data was determined. The point cloud data is initially divided into grids of a preset size to determine the target analysis points. The edge distribution features and regional color and texture features in the local RGB image of each target analysis point are analyzed to obtain the canopy region feature coefficients. The target analysis points are screened based on all the canopy region feature coefficients obtained within the grid to obtain core candidate points. The elevation data of the core candidate points corresponding to the grid are analyzed to determine the elevation undulation coefficient. The initial grid size is adaptively adjusted based on the canopy region feature coefficients and the elevation undulation coefficients to extract the seed terrain skeleton. By utilizing the shape features of each grid after size adjustment and its corresponding seed terrain skeleton, and combining the canopy region feature coefficients corresponding to all target analysis points in each grid, the adaptive stiffness of the CSF algorithm is constructed, and the point cloud data is classified by the CSF algorithm. A CHM model of planted forests was constructed based on the classified point cloud data.
2. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The determination of the local RGB image for each point cloud data specifically involves: Based on the spatial location of the point cloud data and the camera position, combined with the laser divergence angle, the physical surface diameter of each point cloud data is calculated; based on the point cloud data and camera parameters, the sampling distance of the ground corresponding to each point cloud data is calculated. The product of the negative correlation mapping result of the sampling distance of each point cloud data and the diameter of the physical surface is used as the diameter of the local RGB image of each point cloud data. Using the projection coordinates of each point cloud data in the RGB visible light image as the center, draw a circle based on the diameter, and use all the pixels inside the circle as the local RGB image corresponding to each point cloud data.
3. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The determination of the target analysis point specifically involves: Three-dimensional coordinates In Point cloud data falling within each grid area are used as target analysis points; if in There are multiple point cloud data points at the location, The point cloud data with the smallest value is used as the target analysis point at that location.
4. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The obtained canopy region characteristic coefficients are specifically as follows: The proportion of edge pixels in the local RGB image of each target analysis point is used as the edge density, denoted as . ; Extract the principal axis direction of the edge lines in the local RGB image of each target analysis point, and calculate the information entropy of all principal axis directions, denoted as . ; Through formula Calculate the vegetation texture features of the local RGB image corresponding to each target analysis point. Where n represents the number of species along the principal axis. Represents the logarithmic function with base 2; It is a preset constant; Calculate the mean difference between the green over-green index and the red over-red index of all target analysis points in the local RGB image of each target analysis point, and then normalize the result. ; Obtain the mean value of the S-channel pixel value of all pixels in the local RGB image corresponding to each target analysis point in HSV space, and denot it after normalization. ; Will , , The mean value is used as the canopy region feature coefficient of the local RGB image of each target analysis point.
5. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The specific process for obtaining the core candidate points is as follows: Sort the canopy region feature coefficients of all target analysis points corresponding to each grid region in descending order, select the smaller of the two elements with the largest difference between adjacent elements, and select the point cloud data corresponding to all elements sorted after this element as the core candidate points for each grid.
6. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The determination of the elevation fluctuation coefficient specifically involves: Based on all core candidate points corresponding to each grid, a fitting plane is obtained, and the angle between the fitting plane and the plane where the grid is located is obtained. Combining the dispersion of the elevation values of all core candidate points, the elevation fluctuation coefficient of each grid is obtained.
7. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The process of adaptively adjusting the initial mesh size is as follows: Based on the elevation undulation coefficient of each grid and the average canopy region characteristic coefficient of all core candidate points corresponding to the grid, the decomposition decision coefficient of each grid is determined; wherein, the decomposition decision coefficient is positively correlated with the elevation undulation coefficient and negatively correlated with the average canopy region characteristic coefficient. Set a decision threshold. When the decomposition decision coefficient of a grid is greater than the decision threshold, perform quadtree splitting on the grid; otherwise, stop splitting.
8. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The extraction of the seed terrain skeleton specifically involves: Among all the core candidate points corresponding to the grid after the splitting stops, the core candidate point with the smallest sum of the normalized elevation value and the characteristic coefficient of the canopy region is selected as the seed point of the corresponding grid. Calculate the difference between the natural number 1 and the canopy region feature coefficient of each target analysis point in each grid after stopping the split. Use the mean of all the differences as the constraint weight of the corresponding seed point, and obtain the seed terrain skeleton by interpolation.
9. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 1, characterized in that, The process of constructing the adaptive stiffness of the CSF algorithm and classifying point cloud data using the CSF algorithm is as follows: Based on the interpolated continuous seed terrain skeleton surface, for each grid after the splitting stops, the standard deviation of all curvature values within that grid is calculated; the mean of the minimum angle difference between the average aspect of that grid and the average aspect of all adjacent grids is calculated; the average slope within that grid is calculated; and the standard deviation, the mean of the minimum angle difference, and the average slope are all globally normalized and then averaged to obtain the comprehensive terrain complexity index of the corresponding grid, denoted as . ; By interpolating the canopy region characteristic coefficients of all target analysis points within each grid where splitting stops, the canopy disturbance intensity field is obtained, denoted as... ; The specific formula for the adaptive stiffness K of the spatial point corresponding to the seed terrain skeleton surface within each grid of the cloth node in the CSF algorithm is as follows: Where norm() represents the normalization function, This is the difference between the preset maximum stiffness and the preset minimum stiffness. Indicates the preset minimum stiffness; Represents a preset constant; Using the fabric surface after CSF algorithm iteration as a reference surface, the vertical distance between the fabric nodes and the corresponding original point cloud data is calculated. Point cloud data with a vertical distance greater than a set threshold are classified as non-ground point cloud data; point cloud data with a vertical distance less than or equal to the set threshold are classified as ground point cloud data.
10. The method for constructing a digital canopy height model of plantations based on high-precision point cloud data as described in claim 9, characterized in that, The construction of the plantation forest CHM model based on the classified point cloud data includes: Interpolation calculations are performed on ground point cloud data to generate a DEM; Non-ground point cloud data is divided based on DEM raster to obtain the top canopy point cloud, and the top canopy point cloud is interpolated to generate DSM; The DSM and DEM are subtracted pixel by pixel to generate the original CHM; the negative height values in the original CHM are uniformly set to 0; the average height and maximum height of the canopy are statistically analyzed, and abnormal high values that exceed the preset quantile height by a preset multiple are replaced with the height mean of their preset neighborhood units to obtain the corrected CHM. After filtering, the plantation forest CHM model is obtained.