Multi-level dense tree crown segmentation method based on unmanned aerial vehicle image and laser point cloud

Abstract

The application discloses a multi-level dense tree crown segmentation method based on unmanned aerial vehicle images and laser point clouds, and comprises the following steps: collecting forest plot point cloud data and forest plot orthographic image data; filtering the forest plot orthographic image and filtering the forest plot point cloud data; performing normalization processing on the filtered point cloud and removing ground point clouds by using a plane grid method; segmenting the filtered orthographic image by using an improved superpixel algorithm and an energy function; extracting mask coordinates after segmentation, aligning the mask coordinates of a two-dimensional image with normalized point cloud coordinates, and then completing three-dimensional segmentation; determining the number of tree crown vertexes in each initial segmentation region by using a kernel density estimation algorithm based on a long axis; and performing point cloud segmentation by using a Gaussian mixture model to obtain a tree crown separation result. The application can accurately detect a single tree crown in a dense forest plot by combining unmanned aerial vehicle point cloud and image data, is a scientific research method for evaluating a prerequisite of a forest ecosystem, and provides an important driving force for precision forestry.

Description

Technical Field

[0001] This invention relates to the field of forestry technology, specifically a multi-layered dense canopy segmentation method based on UAV imagery and laser point clouds. Background Technology

[0002] Trees are the basic units of forests, and their spatial structure, biological characteristics, and composition are key factors in forest ecological environment modeling, biomass estimation, carbon storage estimation, and resource integration. Individual tree detection is an ideal product for integrating remote sensing data into forest inventory, and individual tree canopy segmentation (ITCs) is an important component of precision forestry. Once a tree is accurately segmented, attributes such as height, species, canopy size, timber volume, and biomass can be obtained. With industrial development and deforestation, forest protection and sustainable development have become important issues in contemporary society. Compared with natural forests, plantations have advantages such as faster growth, higher growth rates, and easier development; therefore, replacing natural forest resources with plantations is an important way to protect forests. Although traditional manual surveying is widely used in forestry, it requires a large amount of time and human resources, has poor real-time performance, and cannot obtain forest land data for large areas.

[0003] The rapid development of remote sensing technology has created favorable conditions for monitoring forest resources. However, the low resolution of satellite images greatly limits the accuracy of canopy segmentation. Meanwhile, UAV remote sensing, due to its ability to fly at low altitudes and acquire high-resolution data, has become the mainstream tool for forestry data acquisition. High-resolution two-dimensional images obtained through RGB sensors possess rich shape and texture information. Segmenting high-canopy-density forest stands using fine spatial resolution images is challenging because a single pixel encompasses different components, especially given the complexity of segmenting individual canopies; edge segmentation is typically at the pixel level. Traditional single-canopy segmentation assumes uniformity in canopy size and shape, but blurs the boundaries between canopies. In multi-layered forests, traditional deep learning methods, such as image segmentation, segment canopies by adjusting parameters. Identifying different locations with different canopy sizes is an iterative and time-consuming process, and the results require visual evaluation and fine-tuning. Its generalization ability is also insufficient. In recent years, deep learning methods using CNNs have been used for vegetation remote sensing. This method does not require image segmentation and feature selection by machine learning classifiers, but it requires a large amount of labeled data and has poor generalization performance. The model performance degrades after changing the data for a different forest stand structure. In addition, its accuracy in forest canopy segmentation is low.

[0004] With the increasing demand for precise forestry, methods for acquiring canopy information from two-dimensional images are no longer sufficient. While 2.5D depth cameras (RGBD) are widely used for indoor modeling and acquiring information over small areas, forest resource integration requires large-area data. Therefore, the focus has shifted to higher dimensions, namely 3D UAV radar, for acquiring forest information. Although this point cloud data possesses good spatial information and three-dimensional structure, and data acquisition is unaffected by the environment, segmenting point clouds generated by airborne lidar in forest areas remains a significant challenge. Existing methods can be broadly categorized into two types:

[0005] The first approach is direct segmentation based on point cloud data. Numerous convolutional neural network segmentation algorithms, such as PointNet, have emerged, demonstrating excellent performance in detecting and segmenting individual tree canopies in sparse forests. In dense forests and complex stands, traditional deep learning algorithms are more effective, combining tree canopy phenotypes and structural geometry. For example, the Random Sample Consensus (RANSAC) algorithm first fits a model to a subset of the point cloud data, then uses this model to classify the remaining points. This process is repeated until a satisfactory segmentation result is achieved. There are also contour model-based segmentation algorithms, which achieve segmentation by minimizing an energy function on the point cloud data. The energy function is defined based on the similarity between neighboring points and can segment point cloud data into different categories. Watershed algorithms treat point cloud data as a terrain surface, simulate flooding, and then segment the point cloud data based on the boundaries of the flooded areas. Other methods include region growing algorithms, voxel-based point cloud analysis algorithms, K-means algorithms, and Support Vector Machine (SVM) algorithms. These methods identify points and classify them into single trees based on simple proximity rules and reliable tree shapes, focusing on the relationships between points rather than their shapes. They are inefficient when processing 3D point clouds and require more computational power compared to image processing-based methods.

[0006] Another approach is based on a two-step segmentation strategy. First, the point cloud is converted into a Digital Surface Model (DSM) or Canopy Cover Model (CHM) and mapped onto a planar raster to construct a two-dimensional digital morphology to describe the forest. Then, computer vision and image processing are combined to segment individual tree canopies. The segmentation part of this method mainly uses gradient-based methods. For example, Li et al. used a watershed algorithm to segment individual tree canopies in coniferous forests using DSMs obtained from LiDAR data. Some researchers have combined region growing and watershed algorithms, using a combination of watershed segmentation and shape-based modeling to segment individual tree canopies from high-density LiDAR data. When the data is more complex, such as when the boundaries of the segmented regions overlap, it is necessary to determine the location of the tree before segmenting the boundaries. A common method is the region growing algorithm with a local search strategy, which determines the location of the tree by finding seed points. Another method is the fixed kernel bandwidth average displacement algorithm, which sets the kernel bandwidth to aggregate the color information of pixels to segment individual tree canopies from complex data. The core concept of these algorithms is mostly based on computer vision, achieving better segmentation by effectively processing point cloud data.

[0007] Although many methods have been proposed for treetop detection and canopy segmentation of laser point cloud data in forests, these methods have certain shortcomings for different types of forests. Canopy segmentation in most plantations is difficult, mainly due to: (1) overlapping canopies, making it difficult to define boundary points, and irregular canopy shapes increase the complexity of segmentation; (2) UAV radar provides a top-down scanning perspective, and its data may be affected by occlusion from other objects in the scene, thus reducing the accuracy of segmentation; (3) variable tree density leads to increased segmentation complexity, and the density of trees may vary throughout the scene, thus affecting the performance of the segmentation algorithm. Therefore, the segmentation algorithm needs to be adaptive to handle trees of different densities. In order to meet the needs of precision forestry and solve the problems of difficult segmentation of ITCs and low overall segmentation accuracy under dense forest conditions, a more accurate method is needed to achieve accurate segmentation of individual canopies and to detect and measure individual canopies in complex forests. Summary of the Invention

[0008] The technical problem to be solved by this invention is to provide a multi-level dense canopy segmentation method based on UAV imagery and laser point cloud to address the difficulties in segmenting ITCs and the low overall segmentation accuracy under dense forest conditions. This multi-level dense canopy segmentation method based on UAV imagery and laser point cloud combines UAV point cloud and image data to accurately detect individual canopies in dense forest stands. It is a scientific research method that is a prerequisite for assessing forest ecosystems and provides an important impetus for precision forestry.

[0009] To achieve the above-mentioned technical objectives, the technical solution adopted by the present invention is as follows:

[0010] A multi-level dense canopy segmentation method based on UAV imagery and laser point clouds includes:

[0011] Step 1: Collect stand point cloud data and stand orthophoto data;

[0012] Step 2: Filter the forest stand orthophoto using Voronoi division, and filter the forest stand point cloud data using bilateral filtering;

[0013] Step 3: Normalize the filtered point cloud and remove the ground point cloud using the planar grid method;

[0014] Step 4: Segment the filtered orthophoto using an improved superpixel algorithm and energy function;

[0015] Step 5: Extract the mask coordinates after segmentation, align the mask coordinates of the 2D image with the normalized point cloud coordinates to complete the 3D segmentation, and obtain multiple initial segmentation segments;

[0016] Step 6: Use the kernel density estimation algorithm to determine the number of crown vertices in each initial segment;

[0017] Step 7: Based on the number of tree canopy vertices, perform point cloud segmentation using a Gaussian mixture model to obtain the tree canopy separation results.

[0018] As a further improvement to the present invention, step 4 specifically comprises:

[0019] Step 4.1: Extract the color distance component, spatial distance component, and texture similarity component from the filtered orthophoto image;

[0020] Step 4.2: Based on the color distance component, spatial distance component, and texture similarity component, and using an improved superpixel algorithm, the orthophoto is segmented into multiple region blocks;

[0021] Step 4.3: Merge region blocks based on energy function;

[0022] Step 4.3.1: Calculate the energy value of each region block, where the region representative point P of the j-th region block... j The local average energy value is simply referred to as the energy value of region j. The energy value of region j is:

[0023]

[0024] in, This represents the energy value of pixel r″ within region block j;

[0025] Step 4.3.2: Select cluster centers:

[0026]

[0027]

[0028] in, It is the energy value of region block j with region block j as the cluster center, where region block j is the representative point P. j Let P be the representative point of the cluster center region. j , It is the average energy value of all regions. The distance between the representative point of a region block in the same cluster and the representative point P of the cluster center region. j distance, It is all Expectations This represents the distance P between the representative point of a region block within the same cluster and the representative point of the cluster center region. j The farthest distance;

[0029] The region block j that satisfies formulas (2) and (3) is the cluster center;

[0030] Conditions for initially determining that they belong to the same cluster:

[0031] The region representative point P of the region block with the highest energy value j The search is performed in eight directions around the center. The region n located in the farthest direction from the center is used as the reference. i The local average energy value of the representative point in the region decreases to a certain value ε, or when the region block n is located at the farthest point... i When the local average energy value of the representative point in the region no longer continues to decrease, then the region block n is taken as... i The regional representative point and the regional representative point P of the region block with the highest energy value j The distance is the radius, and the region represents point P. j Draw a circle with the center as the center, and the regions within the circle belong to the same cluster;

[0032] Optimized clusters:

[0033] With region j as the cluster center, when region b i Energy similarity value with region block j Greater than or equal to the region block n located furthest away in the same cluster i Energy similarity value with region block j At that time, region block b i It belongs to the cluster centered on region block j;

[0034] With region j as the cluster center, when region b iEnergy similarity value with region block j Smaller than the farthest region n in the same cluster i Energy similarity value with region block j At that time, region block b i It does not belong to the cluster centered on region block j;

[0035] Where region block j and region block b i Energy similarity value for: Represents region block b i Energy value; Block n i Energy similarity value with region block j for: Represents region block n i The energy value, K1 is the set of all regions that belong to the same cluster as the cluster center region j;

[0036] Step 4.3.3: Calculate the energy value of the remaining regions that do not belong to a cluster. Then, select cluster centers in the remaining regions according to the method in Step 4.3.2, determine the regions that belong to the same cluster, and optimize the cluster.

[0037] Step 4.3.4: Select cluster centers for the remaining regions that do not belong to a cluster, determine the regions that belong to the same cluster, and optimize the cluster, until all regions belong to the corresponding cluster.

[0038] Step 4.3.5: If a region block belongs to two or more different clusters, calculate and compare the energy similarity values ​​of the region block with the cluster center region blocks in the different clusters to which it belongs. The region block belongs to the cluster to which the cluster center region block with the highest energy similarity value belongs.

[0039] Step 4.3.6: Merge regions belonging to the same cluster.

[0040] As a further improvement to the present invention, in step 6, the formula for the kernel density estimation algorithm is:

[0041]

[0042]

[0043] Where n′ is the number of points in each initial segment, h>0 is a smoothing parameter called bandwidth, K(x′) is the kernel function, and x′ represents the sample point. γ Indicates the observation point.

[0044] As a further improvement to the present invention, in step 7, a Gaussian mixture model is used for point cloud segmentation, and its formula is:

[0045]

[0046] Where N′ represents the number of canopy points estimated by kernel density, X represents the points projected onto the profile, and λ k′ denoted by μ, the proportionality coefficient, representing the prior probability of each mixture component. k′ δ k′ Let represent the initial parameters of the Gaussian distribution, and let represent the mean and variance of the Gaussian function, respectively.

[0047] The beneficial effects of this invention are as follows:

[0048] This invention addresses the challenges of segmenting ITCs (Indoor Tree Canopies) in dense forest conditions, including low overall segmentation accuracy. It develops a high-precision segmentation algorithm by combining the characteristics of UAV orthophotos and point cloud data. First, an improved superpixel algorithm and energy function are developed using UAV orthophoto data for pixel block segmentation and merging, achieving two-dimensional canopy segmentation. The proposed energy function-based image block merging algorithm solves the common oversegmentation problem and achieves better segmentation results. Second, the coordinates of the two-dimensional mask pixels are aligned with the coordinates of the UAV point cloud data to complete three-dimensional dimensionality reduction segmentation. Then, the number of cluster components within each initial segment is determined by kernel density function estimation, and a Gaussian mixture model is used to separate the canopy layer of each tree. Finally, experimental comparisons verify the accuracy and generalization of the proposed segmentation method. This algorithm achieves high segmentation accuracy not only in sparse forests but also in dense forests, demonstrating high segmentation precision. In summary, the method of this invention improves the segmentation accuracy of canopy edges by using an improved superpixel algorithm and energy function, refines point cloud segmentation to improve accuracy by using kernel density function estimation and Gaussian mixture model, and promotes canopy segmentation towards a better solution by detecting each mask. Attached Figure Description

[0049] Figure 1 A location map of the study site.

[0050] Figure 2 This is a flowchart of the initial segmentation process for orthophotos based on an improved superpixel algorithm and energy function.

[0051] Figure 3 This is a segmentation map based on an improved superpixel algorithm (with the addition of an edge probability function).

[0052] Figure 4 This is an existing superpixel segmentation map.

[0053] Figure 5 This is a graph of the normalized function.

[0054] Figure 6 This is a schematic diagram of image block (region block) merging.

[0055] Figure 6 (a) in the diagram is a schematic diagram of two different clusters that do not have a common region block.

[0056] Figure 6 (b) in the diagram is a schematic diagram of two different clusters sharing a common region block.

[0057] Figure 7 This is a fine segmentation diagram of the point cloud.

[0058] Figure 7 (a) in the diagram is a schematic diagram of the normal single-crown point cloud.

[0059] Figure 7 (b) in the diagram is a schematic diagram of an abnormal single crown point cloud.

[0060] Figure 7 (c) in the diagram is a schematic diagram of the treetop detection point cloud.

[0061] Figure 7 (d) in the figure is a schematic diagram of the treetop detection histogram.

[0062] Figure 7 (e) in the figure is the result of segmenting a single tree canopy point cloud.

[0063] Figure 8 The graph shows the histogram of point density distribution and the estimated curve of Gaussian kernel density.

[0064] Figure 8 (a) in the figure is the histogram of point density distribution corresponding to the statistical interval of 0.3.

[0065] Figure 8 (b) in the figure is the histogram of point density distribution corresponding to the statistical interval of 0.5.

[0066] Figure 8 (c) in the figure is the histogram of point density distribution corresponding to the statistical interval of 0.7.

[0067] Figure 8 (d) in the figure is the Gaussian kernel density estimation curve corresponding to a bandwidth of 0.3.

[0068] Figure 8 (e) in the figure is the Gaussian kernel density estimation curve corresponding to a bandwidth of 0.5.

[0069] Figure 8 (f) in the figure is the Gaussian kernel density estimation curve corresponding to a bandwidth of 0.7.

[0070] Figure 9 This is a schematic diagram of the experimental results.

[0071] Figure 9 (a) in the figure is a schematic diagram of the canopy segmentation results of the low-density sample plot using the normalized UAV point cloud method.

[0072] Figure 9 (b) in the diagram is a schematic diagram of the canopy segmentation results of the low-density sample plot using the local maximum algorithm.

[0073] Figure 9 (c) in the figure is a schematic diagram of the tree canopy segmentation result using the point cloud segmentation algorithm in low-density sample plots.

[0074] Figure 9 (d) in the diagram is a schematic diagram of the canopy segmentation results of the seed point stacking algorithm in the low-density sample plot.

[0075] Figure 9 (e) in the figure is a schematic diagram of the canopy segmentation results of the low-density sample plot using the algorithm studied in this paper. Figure 9 (f) in the figure is a schematic diagram of the canopy segmentation results of the medium-density plot using the normalized UAV point cloud method.

[0076] Figure 9 (g) in the diagram represents the canopy segmentation result of the local maximum algorithm in the medium-density plot.

[0077] Figure 9 (h) in the figure is a schematic diagram of the canopy segmentation result of the point cloud segmentation algorithm in the medium density plot.

[0078] Figure 9 (i) in the diagram is a schematic diagram of the canopy segmentation result of the seed point stacking algorithm in the medium-density sample plot.

[0079] Figure 9 In the diagram, (j) represents the canopy segmentation results of the medium-density sample plot using the algorithm studied in this paper. Figure 9 In the diagram, (k) represents the canopy segmentation results of high-density plots using the normalized UAV point cloud method.

[0080] Figure 9 (l) in the figure is a schematic diagram of the canopy segmentation result of the local maximum algorithm in the high-density sample plot.

[0081] Figure 9 The (m) in the diagram represents the canopy segmentation result of the point cloud segmentation algorithm used in the high-density sample plot.

[0082] Figure 9 In the diagram, (n) represents the canopy segmentation result of the seed point stacking algorithm in high-density sample plots.

[0083] Figure 9 (o) in the diagram represents the canopy segmentation results of high-density sample plots using the algorithm studied in this paper.

[0084] Figure 10 Scatter plots of measured crown width and extracted crown width obtained by four different methods.

[0085] Figure 10 (a) in the figure is a scatter plot of the measured crown width and the extracted crown width obtained by using the local maximum algorithm in the low-density plot.

[0086] Figure 10 (b) in the figure is a scatter plot of the measured crown width and the extracted crown width obtained by the point cloud segmentation algorithm for low-density plots.

[0087] Figure 10 (c) in the figure is a scatter plot of the measured crown width and the extracted crown width obtained by using the seed point stacking algorithm in low-density plots.

[0088] Figure 10 In the diagram, (d) is a scatter plot of the measured crown width and the extracted crown width obtained by the algorithm studied in this paper for low-density plots.

[0089] Figure 10 In the diagram, (e) is a scatter plot of the measured crown width and the extracted crown width obtained by using the local maximum algorithm for medium-density plots.

[0090] Figure 10 In the diagram, (f) represents the scatter plot of the measured crown width and the extracted crown width obtained by the point cloud segmentation algorithm for medium-density plots.

[0091] Figure 10 In the diagram, (g) represents a scatter plot of the measured crown width and the extracted crown width obtained using the seed point stacking algorithm for medium-density plots.

[0092] Figure 10 In the diagram, (h) represents a scatter plot of the measured crown width and the extracted crown width obtained using the algorithm studied in this paper for medium-density plots.

[0093] Figure 10 In the diagram, (i) is a scatter plot of the measured crown width and the extracted crown width obtained by using the local maximum algorithm in the high-density plot.

[0094] Figure 10 In the diagram, (j) represents the scatter plot of the measured crown width and the extracted crown width obtained by the point cloud segmentation algorithm for high-density sample plots.

[0095] Figure 10 In the diagram, (k) represents a scatter plot of the measured crown width and the extracted crown width obtained by using the seed point stacking algorithm in high-density plots.

[0096] Figure 10In the diagram (l), the measured crown width and the extracted crown width are scatter plots obtained by using the algorithm studied in this paper on high-density plots. Detailed Implementation

[0097] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0098] To address the challenges of segmenting ITCs in dense forest conditions and the overall low segmentation accuracy, this paper develops a high-precision segmentation algorithm that combines the characteristics of UAV orthophotos and point cloud data.

[0099] (1) First, an improved superpixel algorithm and energy function were developed using UAV orthophoto data to segment two-dimensional tree canopies by pixel block segmentation and merging.

[0100] (2) Next, the coordinates of the two-dimensional mask points are aligned with the coordinates of the UAV point cloud data to complete the three-dimensional dimensionality reduction segmentation.

[0101] (3) Then, the number of cluster components in each initial segment is determined by kernel density function estimation, and the canopy of each tree is separated by Gaussian mixture model separation method.

[0102] (4) Finally, the accuracy and generalization of the segmentation method in this study were verified through experimental comparison. The algorithm not only has a high segmentation accuracy in sparse forest land, but also has a high segmentation result in dense forest stands.

[0103] 1. Materials and Methods:

[0104] 1.1 Study Area:

[0105] The data used in this study came from Hongze Lake (31°32″N, 118°89″E), collected in August 2023. This region belongs to the subtropical monsoon climate zone, with abundant rainfall, ample sunshine, an average annual temperature of 15.5℃, and an average annual precipitation of 1019.5 mm, suitable for plant growth. The region is dominated by plantations with high canopy density and complex tree species composition, featuring two very distinct vegetation layers: a tree layer, mainly composed of willows (average height 6.72m), and a grass layer with a variety of annual organisms. After decades of forest management, logging, and reforestation, the original tree distribution, stand age, and vertical structure, which were relatively uniform, have evolved into a multi-layered and complex forest structure. This different stand structure presents a significant challenge to the segmentation of ITCs (Individual Tree Cities). To ensure the validity of the experiment, this study used data from areas with shorter planting times and simpler stand structures for comparison, selecting three densities of willows as experimental subjects in both areas. The selected areas are as follows: Figure 1 As shown.

[0106] 1.1.1 Remote sensing data and field data:

[0107] This study used the DJI M600 Pro to acquire imagery and point cloud data. Equipped with a triple-redundant A3 Pro flight controller, Lightbridge 2 high-definition digital image transmission, intelligent flight battery pack, and battery management system, it provides centimeter-level precise positioning. Data was acquired using the gAirHawk GS-100C acquisition system, which integrates a Livox laser scanner, an airborne positioning and orientation (POS) system, and a visible light camera, weighing less than 1050g. Remote sensing imagery data was acquired using the SONY A6000 camera built into the M3M drone system manufactured by DJI Innovations Technology Co., Ltd. in Shenzhen. The camera's main parameters are: 20 million effective pixels, 83° field of view; maximum resolution 5280*3956; image size: per... Figure 6-2 4MB; Focal length: 25mm; Acquired along an S-shaped trajectory for subsequent registration and correction. UAV point cloud data was simultaneously acquired using a GS-100C laser sensor system, which has a ranging accuracy of ±20mm and a laser pulse frequency of 240kHz. The wavelength and beam divergence are approximately 1550nm and 0.5mrad, respectively. A frequency of 720,000 points / second for three echoes ensures good vegetation penetration. The Positioning and Orientation (POS) system uses position data acquired from the Global Positioning and Navigation Satellite System (GNSS) as initial values, combined with an Inertial Measurement Unit (IMU) to acquire attitude change increments. Through Kalman filters and feedback error control, iterative calculations generate navigation data in real time, achieving a positioning accuracy of 0.02-0.03m. It can acquire the exterior orientation elements of the camera and the absolute position of the aircraft, enabling high-precision direct ground positioning without ground control, facilitating further data processing. The ground-based measurement data for the study area were completed collaboratively by our team in August 2023. The measured factors included tree species, diameter at breast height (DBH), tree height, and crown width in the east-west and north-south directions.

[0108] The UAV data acquisition process involves importing photos and POS data into Pix4Dmapper software, performing image registration, spatial calculation, adjustment, and geometric correction, and finally performing image mosaicking and orthorectification to complete the production of orthophotos. To ensure the validity of the high-resolution digital orthophoto results, the UAV flight trajectory designed in this study has a heading overlap of 60% and a lateral overlap of 20%.

[0109] After completing outdoor scanning operations, the airborne laser scanning system acquired GNSS data from the airborne unit and base station, as well as IMU inertial navigation data. It then obtained POS data corresponding to the spatial position and attitude at each moment. The POS data was then fused with the original waveform data to generate a point cloud LAS file, using WGS84 as the coordinate system. The processing was performed using the point cloud computing software gAirHawk.

[0110] 1.2 Methodological Framework:

[0111] The overall workflow of this study comprises three main steps: First, the acquisition and preprocessing of UAV imagery data, LiDAR data, and field data at three different densities. Second, an improved superpixel segmentation method is used to segment the images into small pixel blocks, and an energy function is used to merge these blocks. Then, the segmentation mask is aligned with the pixel coordinates of the simultaneously acquired 3D point cloud data to complete 3D segmentation. Next, a kernel density estimation algorithm is used to detect the tree canopy position, and a Gaussian mixture model is combined to complete the fine segmentation of the 3D tree canopy. Finally, the 3D segmentation results are evaluated, and the understory trees are segmented using the spatial advantages of the 3D point cloud. The effectiveness of the model is verified by comparing the results with LM, PCS, and LS algorithms, and the generalization ability of the model is verified by comparing it at different densities.

[0112] 1.2.1 Data Processing:

[0113] Image preprocessing:

[0114] The watershed algorithm, most commonly used in optical remote sensing imagery for segmenting individual trees in coniferous forests or sparse woodlands, is a terrain-based image segmentation algorithm capable of identifying regions and boundaries in the image. This algorithm avoids the influence of noise and typically first uses local maximum filtering to extract feature points in the canopy region, or uses Gaussian filtering to eliminate internal canopy texture and noise. Then, a multi-scale watershed algorithm is used to segment mixed deciduous forests. These methods have good segmentation results for forests with low canopy closure. However, local maximum filtering is often limited by the size of the filtering window; a window that is too large will cause small trees to be missed, while a window that is too small will cause a single canopy to be divided into multiple canopies. Applying this method often requires smoothing the original image to eliminate the influence of canopy internal structure and noise. However, this processing method, like direct Gaussian filtering, will blur canopy boundaries, easily causing oversegmentation and undersegmentation errors. Currently, there is little research on the segmentation of tree canopies in plantations with high tree outline clarity. Dense forest canopies are characterized by complex canopy textures, irregular arrangement, lack of prominent local high points, and similar canopy widths. In addition, the segmentation of high-density forests involves confusion between the ground and tree canopy outlines, severe canopy overlap, and difficulty in clearly defining individual trees.

[0115] To address the above issues, this paper references the principle of Voronoi diagram proximity property partitioning, dividing pixels into regions based on color threshold constraints for filtering. First, the image is converted to Lab color space, generating coordinates for each pixel. Then, a color range within each region is defined. Pixels are traversed, and their respective regions are determined based on the color range constraints. The specific method is as follows: a) Set upper and lower limits for the color threshold based on image characteristics; b) Check if the color threshold of each pixel in the image falls within the set threshold range. Pixels that simultaneously meet both conditions are grouped into the same region. This method can quickly separate grassland from forest areas and also makes boundary features more distinct.

[0116] ALS preprocessing:

[0117] Before fine segmenting the point cloud, it needs to be normalized. Therefore, a Canopy Elevation Model (CHM) needs to be generated first. The CHM is a two-dimensional surface model representing the height of the upper canopy surface above the ground. It eliminates the influence of topographic relief on canopy height, thus reflecting the fluctuations in forest canopy height, and plays an important role in the inversion of forest parameters or forest biomass. This paper uses ArcGIS software to perform bilateral filtering denoising on the imported point cloud and uses planar grid segmentation to segment the ground. Then, a Digital Elevation Model (DEM) and a Digital Surface Model (DSM) are generated. The CHM model is generated by subtracting the elevation value of each point in the DSM from the corresponding ground elevation values, and finally, the normalized point cloud is obtained.

[0118] 1.2.2 Extraction of region center points based on boundary energy function:

[0119] This paper proposes a novel UAV image segmentation algorithm. Within each sub-region of a UAV image, samples share certain common attributes, such as color, brightness, and texture features. Based on different methods, image segmentation algorithms can be broadly categorized into four types: (a) clustering-based algorithms; (b) region merging-based algorithms; (c) graph theory-based algorithms; and (d) classification-based algorithms. This paper employs a region-based method. The algorithm first performs superpixel pre-segmentation, then introduces an energy function for block merging, effectively addressing the problem of blurred edge feature extraction in traditional algorithms. Its edge energy function is expressed as:

[0120] EE=B Energy {j D(Edge)} (1);

[0121] Where, j D(Edge) This represents the improved superpixel algorithm segmenting region block j, and then merging multiple region blocks into the final result based on the energy function.

[0122] The initial segmentation flowchart of the unmanned orthophoto image is as follows: Figure 2 As shown.

[0123] Image patch pre-segmentation:

[0124] This paper extracts Lab color space and spatial location features to construct a 5-dimensional feature vector to measure similarity, thus solving the segmentation boundary blurring problem based on pixel local clustering. It also utilizes the image's color, coordinates, and texture information to construct a 7-dimensional vector.

[0125] First, assume the original input image I has N pixels, and it is expected to be segmented into K pixel blocks, with a block size of N / K. The distance between each cluster center is defined as... During the initialization of clustering, the gradient value of each pixel is calculated within a window of size 3×3 pixels. The pixel with the smallest gradient value is found within this window, and the cluster center is placed at that position.

[0126] Then, each superpixel is treated as a basic unit of image processing, and the Euclidean distance d between any two superpixels is calculated. c (x,k) represents the Euclidean distance between pixel k and cluster center x in the Lab color space; the closer the two are, the smaller the value of C(x,k); d l (x,k) represents the spatial distance from pixel k to cluster center x.

[0127] Texture similarity is determined by the x-axis of its texture histogram. 2 (c x ,c k Distance, defined as the similarity of textures, also known as the histogram distance of texture elements, is:

[0128]

[0129] Where c x (i)=∑ k∈w(x) I[T(k)=i],c x (i) represents the frequency of a certain type of fringes within a certain window w(x) (fringe histogram), where w(x) represents a window, T(k) represents a pointer function, and c x c k These represent two histograms of the two elements, where x 2 (c x ,c k The smaller the value, the stronger the c. x c k The closer.

[0130] Define d w (x,k) is a measure of texture similarity:

[0131] d w (x,k)=x 2 (c x ,c k(3);

[0132] c x c represents the histogram of the texels representing the cluster centers. k This represents the histogram of pixels k within a certain range centered on the cluster center.

[0133] This paper incorporates texture feature components into the overall metric calculation, using the boundary probability function g(x) to represent the probability that pixel x lies on the image boundary. The constructed overall objective function formula is as follows:

[0134]

[0135]

[0136] Where D is the similarity (i.e., distance) between a pixel and a cluster center, d w (x,k) represents the texture similarity (i.e., texture distance) between pixel k and cluster center x; d c (x,k) represents the Euclidean distance between the color vectors of pixel k and cluster center x in the Lab color space; d l (x,k) represents the Euclidean distance of the spatial vector from pixel k to cluster center x, where d c (x,k) and d l (x,k) is calculated using existing formulas, d w (x,k) is calculated using formulas (2) and (3); λ1 is a constant. g(x) is the boundary probability function, representing the probability that pixel x is on the image boundary. This term assigns different weights to the color distance component C(x,k), spatial distance component L(x,k), and texture similarity component W(x,k) for different pixels, and the weights are determined by g(x). When the probability of a pixel being on the boundary is high, the color distance between the pixel and the seed point should be considered more during the clustering process to ensure that the superpixels generated by the segmentation fit the image edge better. Therefore, the weight of the texture component W(x,k) should be increased, and the weights of the spatial distance component L(x,k) and the color similarity component C(x,k) should be decreased.

[0137] Here, the edge probability g(x) is obtained by mapping the image gradient magnitude f(x) through an edge probability function. The gradient magnitude of a pixel can reflect, to some extent, the probability that the pixel falls on the image boundary. Generally, the gradient is smaller in flat areas of the image and larger at the image edges. Although the gradient magnitude can characterize the probability that a pixel falls on the image boundary, the gradient information varies greatly between different images, and the range of gradient magnitude values ​​is large, making it difficult to directly use the gradient magnitude to characterize the probability that a pixel is on the edge. Therefore, this study needs a normalization function to reasonably map the gradient magnitude to obtain the edge probability of the pixel, using a function such as... Figure 3 The normalization function shown maps the gradient magnitude to the range (0, 1), i.e. The function is centrally symmetric about (m, 0.5). When the gradient magnitude is in the range (0, t), it indicates that the gradient magnitude of the pixel is too small, and the pixel can be considered not to be at the image edge. In this case, the boundary probability mapping function maps the gradient magnitude to a value close to 0. When the gradient magnitude of the pixel is in the range (T, +∞), it indicates that the gradient magnitude of the pixel is large enough, and the pixel can be considered to be at the image edge. The boundary probability function maps it to a value close to 1. Here, a and m are hyperparameters. In this study, the values ​​of the two thresholds t and T can be adjusted by adjusting the values ​​of a and b.

[0138] Using formula (4), calculate the similarity between the cluster center and the pixels within its 2S×2S window, and assign the label of the most similar cluster center to the pixel corresponding to the smallest D value (dividing each pixel into the region with the highest similarity). Update the cluster centers iteratively until convergence. To increase coherence, smaller superpixel blocks are incorporated into adjacent pixel blocks.

[0139] The main steps of the improved superpixel segmentation in this paper are as follows:

[0140] Step 1: Initialize cluster centers. Input image I, set the expected number of superpixels for segmentation N, the size of each superpixel block is N / K, and the cluster center distance is...

[0141] Step 2: Adjust cluster centers. To prevent the cluster centers to be located at the image edges, the gradient values ​​of the pixels are calculated within a 3×3 pixel window of the initial cluster center point. The pixel with the smallest gradient value within this window is then placed as the cluster center.

[0142] Step 3: Calculate similarity. Use formula (4) to calculate the similarity (i.e., distance) between the cluster center and every pixel within its 2S×2S window, and assign the label of the most similar cluster center to the pixel corresponding to the smallest D value.

[0143] Step 4: Update the cluster centers. Repeat step (3) until convergence.

[0144] Step 5: Increase connectivity. For smaller superpixel blocks, the principle is to include them in the nearest adjacent superpixel blocks.

[0145] The segmentation results are as follows Figure 4 As shown. Figure 4 The boundary point segmentation effect at point a is significantly better than that at point 'a'. Figure 5 The value of a′ reflects the effectiveness of incorporating the marginal probability function. Figure 4 The pixel segmentation effect at point b is significantly better than Figure 5 b′, when a pixel is located at a boundary point, separates the overlapping edges of the tree canopy based on the difference in texture components. This also ensures the fit between the superpixel boundary and the image boundary.

[0146] Image patch merging based on energy function:

[0147] This study found two characteristics in the tree canopy data: 1) Although the tree canopy itself is relatively complex, it contains rich local information, making the tree canopy layers relatively easy to identify; 2) The regions with rich local information are far apart; 3) The tree canopy layers are quasi-circular, and this study can define the merged shape. Therefore, to more clearly identify local values, this study enhanced the image and converted it into a normalized energy map.

[0148] In pre-segmentation, an image is divided into K regions of varying shapes. Assume that the cluster center of region j is P. j This image can be represented as: P = (p1, p2, ..., p j ,…,p M The set of points in region j is in Let n be the nth point in region j. Then the representative point of this region is P. j In this study, the energy value at each point is defined as... Then the cluster center of region block j is the region representative point P of region block j. j The energy value of region block j is defined as the region representative point P. j The local average energy value, the energy value of region block j is:

[0149]

[0150] This represents the energy value of pixel r″ within region block j. In other words, the energy value of region block j is the sum and average of the energy values ​​of all pixels in region block j.

[0151] The steps of energy function-based image patch merging include:

[0152] Step a: Calculate the energy value of each region block according to formula (5); where the energy value of a region block is the sum and average of the energy values ​​of all pixels in the region block;

[0153] Step b, Step b1, Selecting cluster centers:

[0154]

[0155]

[0156] in, It is the energy value of region block j with region block j as the cluster center, where region block j is the representative point P. j Let P be the representative point of the cluster center region. j , It is the average energy value of all regions. The distance between the representative point of a region block in the same cluster and the representative point P of the cluster center region. j distance, It is all Expectations This represents the distance P between the representative point of a region block within the same cluster and the representative point of the cluster center region. j The farthest distance;

[0157] The region block j that satisfies formulas (6) and (7) is the cluster center;

[0158] Conditions for initially determining that they belong to the same cluster:

[0159] From all regions, select the representative point P of the region with the highest energy value. j The search is performed in eight directions around the center. The region n located in the farthest direction from the center is used as the reference. i The local average energy value of the representative point in the region decreases to a certain value ε, or when the region block n is located at the farthest point... i When the local average energy value of the representative point in the region no longer continues to decrease, then the region block n is taken as... i The regional representative point and the regional representative point P of the region block with the highest energy value j The distance is the radius, and the region represents point P. j Draw a circle with the center as the center, and the regions within the circle belong to the same cluster;

[0160] Step b2, optimize the cluster:

[0161]

[0162] B jThis indicates that multiple regions centered at region j belong to the same cluster;

[0163] Using region block j as the cluster center, when step b1 determines that a certain region block b belongs to the same cluster as region block j... i Energy similarity value with region block j Greater than or equal to the region block n located furthest away in the same cluster i Energy similarity value with region block j hour, A value of 1 indicates that the region is block b. i It belongs to the cluster centered on region block j;

[0164] Using region block j as the cluster center, when step b1 determines that a certain region block b belongs to the same cluster as region block j... i Energy similarity value with region block j Smaller than the farthest region n in the same cluster i Energy similarity value with region block j hour, A value of 0 indicates that the region is block b. i If a region does not belong to a cluster centered on region block j, remove region block b from the cluster. i ,like Figure 6 In the diagram, the area marked with a cross is the area that has been deleted.

[0165] Where region block j and region block b i Energy similarity value for: Represents region block b i Energy value, This represents the energy value of region block j; region block n i Energy similarity value with region block j for: Represents region block n i Energy value, Middle molecule Represents region block n i The energy difference value between region block j and region block j, where K1 represents the set of all region blocks that belong to the same cluster as the cluster center region block j, as determined by step b1. This represents the energy value of region block f1, where region block f1 belongs to K1; The meaning of the denominator is the energy difference value between all regions belonging to K1 and region j. Perform summation.

[0166] Step c: Calculate the energy value of the remaining regions that do not belong to a cluster; select cluster centers in the remaining regions according to the method in step b, determine the regions that belong to the same cluster, and optimize the cluster.

[0167] Step d: Select cluster centers for the remaining regions that do not belong to a cluster, determine the regions that belong to the same cluster, and optimize the cluster until all regions belong to the corresponding cluster.

[0168] Step e: If a region block belongs to two or more different clusters, such as Figure 6 (b) shows that if two circles have overlapping regions, the energy similarity values ​​of the region block and the cluster center regions of different clusters to which they belong are calculated and compared. The region block belongs to the cluster to which the cluster center region with the highest energy similarity value belongs.

[0169] Step f: Merge multiple region blocks belonging to the same cluster. After merging multiple region blocks, a large region block is formed. Therefore, after merging in step f, multiple region blocks are obtained. However, the number of region blocks is greatly reduced compared with the number of region blocks obtained based on the improved superpixel segmentation, thus solving the problem of oversegmentation in the improved superpixel algorithm.

[0170] In the above process, the high-energy region is used as a representative point P. j Centered on the target area, searches are conducted in eight directions, and the furthest distance is selected for merging. Among the four directions, the furthest direction is used as the benchmark. When the energy value of a certain region drops to a certain value ε, or when the energy value of a certain region stops decreasing, the region's representative point and the high-energy center point P are used as the reference. j The distance is the radius. Within the pre-merging circle, there may be other areas where the energy values ​​do not meet the requirements. Therefore, further filtering is performed according to the above conditions. Figure 6 As shown in (a); there may also be a situation where a region is within the pre-merging range of two canopies simultaneously. In this case, the region is divided according to the energy similarity value between the region block and the cluster center region blocks (region blocks that are the cluster centers) in different clusters (different circles) to which it belongs. Figure 6 As shown in (b) of the diagram. Figure 6 (a) in the text indicates a merging method for objects that are far apart. Figure 6 (b) in the text represents the area within which two tree canopies merge simultaneously.

[0171] The proposed image patch merging algorithm based on energy function solves the common oversegmentation problem and achieves better segmentation results.

[0172] 1.2.3 Radar point cloud segmentation:

[0173] The two-dimensional tree canopy ring mask was extracted using ArcGIS software. Based on the characteristic that the point cloud and image coordinate points correspond to each other under the same spatial and temporal shooting, the raster data of the mask was matched with the CHM of the UAV radar point cloud for CHM segmentation. Points within the boundary range in the CHM were picked with the vector boundary as a reference and randomly assigned different colors. In this way, two-dimensional segmentation promotes three-dimensional segmentation, thereby achieving the effect of dimensionality reduction segmentation.

[0174] Kernel density estimation function for canopy number detection:

[0175] During the field survey, the location of each tree was accurately measured at the base of its trunk. Because the point density of airborne lidar point clouds typically decreases from top to bottom, this results in a relatively small number of points representing the structural details of the trunk and below compared to other parts of the tree—a lack of understory information. This is particularly problematic in overlapping forests, making it difficult for airborne lidar point clouds to accurately pinpoint the location of individual tree trunks. Similar to the trunk, the tree vertices are also used as markers for individual trees, representing a specific number of individual trees within the point cloud. Figure 7 As shown in (c), the density of points at the center of each tree is usually very high, and the point density decreases from the center to both sides. The projection contour of a single tree crown on the horizontal plane (i.e., the orthophoto segmentation contour) is generally approximately circular, while the projection contour and spatial distribution of multiple tree point clouds are approximately elliptical. Therefore, segmenting point clusters according to the major axis can better detect multiple trees. Thus, this paper proposes kernel density estimation along the major axis to locate the potential tree tops, i.e., finding the major axis of the ellipse projected along the tree crown, and using the major axis as the x-axis to perform kernel density estimation to find the tree crown location. To accurately detect local maxima of point density, kernel density estimation is used to calculate the probability density function distribution of each initial segment. Kernel density estimation is defined as:

[0176]

[0177]

[0178] Where n′ is the number of points in each initial segment, h>0 is a smoothing parameter called bandwidth, K(x′) is the kernel function, and x′ represents the sample point. γ To represent the observation point, this paper uses the Gaussian kernel function for density estimation.

[0179] Gaussian mixture model for tree canopy segmentation:

[0180] After coarse segmentation, fine segmentation is performed on the clusters under each mask. The K-means algorithm is commonly used for 3D point cloud segmentation, but it has significant limitations for segmenting non-convex data. Therefore, this paper uses a Generic Model (GMM) suitable for complex numbers and non-convex data to model the canopy point cloud. This method classifies data points into different categories based on their attribute columns to achieve fine segmentation of the 3D point cloud. When projecting onto the 3D point cloud, this study constructs a local coordinate reference frame for each cluster under the mask. After constructing the local coordinate reference frame for the current point cluster, the clusters are divided into three groups from 0° to 180° in a counter-clockwise direction at 60° intervals for further analysis. Starting from 0° (the x-axis direction), the line connecting the farthest points along the x-axis on the xoy projection plane is denoted as d1. The distance to the farthest point in the 60° direction is denoted as d2, and the distance to the farthest point in the 120° direction is denoted as d3. This means that this study will analyze the distance of each cluster in the 0°, 60°, and 120° directions during implementation. This study selects the horizontal direction of the profile projection with the maximum distance in three directions as the x-axis, i.e., max{d1,d2,d3}. Figure 7 (b) is a multi-directional three-dimensional spatial structure diagram when rotated 0° counterclockwise, with kernel density estimation performed along the x-axis.

[0181] Then, point cloud segmentation is performed based on the crown vertices estimated by kernel density (referred to as crown points) and combined with a Gaussian mixture model. The formula can be expressed as:

[0182]

[0183] Where N′ represents the number of canopy points estimated by the kernel density above, and X represents the point projected onto the profile (x p ,z p ), λ k′ denoted by μ, the proportionality coefficient, representing the prior probability of each mixture component. k′ δ k′ Let represent the initial parameters of the Gaussian distribution, and let represent the mean and variance of the Gaussian function, respectively. The algorithm calculates the component probabilities using the E-step and updates the Gaussian mixture parameters μ using the M-step. k′ δ k′ The mixture distribution parameters are calculated iteratively. When the parameter value is less than the threshold or the number of iterations reaches the maximum, the EM step converges, thereby classifying the point cloud. Figure 7 (e) represents the result of optimization to two tree canopies.

[0184] Figure 7 This is a fine-grained segmentation image of the point cloud, in which Figure 7 In the text, (a) represents a normal single crown. Figure 7 (b) in the text indicates an anomalous single crown, meaning it includes two crowns. Figure 7 In the diagram, (c) and (d) represent treetop detection, and (e) represents the result of individual timber segmentation.

[0185] 1.2.4 Verification Indicators:

[0186] To avoid errors in manual canopy measurement caused by complex terrain, high canopy density, and irregular canopy growth, this study used digital orthophoto data as the base map and manually delineated the canopy boundaries to assess accuracy. The position of individual trees obtained from ITCs segmentation was compared with the true canopy overlap rate to determine the number of true positives (TP), false positives (FP), and false negatives (FN). Segmentation was considered correct when the overlap rate between the reference canopy and the segmented canopy was ≥70%, and incorrect when the overlap rate was ≤70%. This study calculated the ITS results for three sample blocks, obtaining TP, FN, and FP. The F_score (f, overall accuracy) was also calculated to assess the overall accuracy considering errors and omissions.

[0187]

[0188]

[0189]

[0190] In the formula, r represents the detection rate of a single tree; p represents the accuracy of detection of a single tree; and f is the overall accuracy of ITCs segmentation.

[0191]

[0192]

[0193] In the formula, g is the number of correctly split individual trees; X τ This indicates the crown width of a correctly segmented individual tree; x represents the average width of a single tree. τ This represents the measured crown width of a single segmented tree. The measured average crown width, representing the crown width of a single segmented tree, is compared with the measured crown width obtained from the segmentation method used in this study. A relative coefficient R0 is used. 2 The root mean square error (RMSE) is used to evaluate the width of the extracted crown.

[0194] 2. Results:

[0195] 2.1 Selection of kernel density estimation parameters:

[0196] The kernel density distribution curves of different trees clustered together change depending on the size of the kernel bandwidth. Therefore, appropriate kernel bandwidth parameters can be selected through cross-validation or based on domain knowledge. This paper conducts experiments on 13 tree canopy groups and selects the optimal kernel bandwidth parameter based on classification accuracy.

[0197] Table 1:

[0198]

[0199] This article is based on Figure 7 Taking the two canopy divisions in the example, Figure 8 The density histograms for statistical intervals of 0.3, 0.5, and 0.7, and the kernel density estimation curves for this bandwidth, show that the canopy point can be found correctly when the bandwidth is 0.5, i.e., the optimal bandwidth is 0.5.

[0200] Figure 8 This represents the histogram of point density distribution and the Gaussian kernel density estimation curve. Figure 8 The statistical interval in (a) is 0.3. Figure 8 The statistical interval in (b) is 0.5. Figure 8 The statistical interval in (c) is 0.7. Figure 8 The bandwidth in (d) is 0.3. Figure 8 The bandwidth in (e) is 0.5. Figure 8 The bandwidth in (d) is 0.7.

[0201] 2.2 Comparison of ITCS segmentation under different point cloud densities:

[0202] The samples were divided into three different densities: low-density samples, medium-density samples, and high-density samples. Figure 9 In the previous experiment, four different segmentation methods were used to segment trees of three densities. The Local Maximum Algorithm (LM algorithm) and the Seed Stacking Algorithm (LS algorithm) can also accurately identify the position of the tree top. This experiment, however, accurately identifies the boundaries and segments according to their vertical characteristics, and uses different colors to visualize the results. Figure 9 The low-to-medium density plots were sparse and evenly planted willow groves, with some trees being saplings and having small individual canopies. The medium-density plots were dense but evenly planted willow groves, with more mature willows. The high-density plots were dense but unevenly planted willow groves, with some trees surrounded by newly sprouted saplings. Among the four individual canopy segmentation methods, Local Maximum (LM) and Seed Point Stacking (LS) methods achieved satisfactory tree location detection results. In groves with large canopies, Seed Point Stacking...

[0203] The LS algorithm detected multiple trees with small canopies, but its segmentation of some adjacent trees was insufficient. The point cloud segmentation method (PCS algorithm) can detect tree canopies and accurately represent the actual size of individual canopies, but this method segments the canopy size of individual trees relatively uniformly, meaning it cannot segment unevenly distributed forests well, i.e., its segmentation effect is poor in high-density plots. The method in this study, however, achieves both a better detection rate and excellent canopy boundary delineation.

[0204] Figure 9 This is a schematic diagram of the experimental results. Figure 9 In the middle, (a)-(e) represent low-density plots, (f)-(j) represent medium-density plots, and (k)-(o) represent high-density plots. (a), (f), and (k) represent normalized UAV point clouds. (b), (g), and (l) represent local maximum algorithms. (c), (h), and (m) represent point cloud segmentation algorithms. (d), (i), and (n) represent seed point stacking algorithms. (d), (j), and (n) represent the algorithms used in this study.

[0205] The accuracy of segmentation was evaluated by matching the segmentation results of individual trees obtained by four segmentation methods in three sample plots with the measured location data of each tree in the sample plots, as shown in Table 2. Significant differences in segmentation accuracy were observed between different methods for ITCs. The fluctuation ranges of r, p, and f for the three different tree densities were 0.81–0.98, 0.71–0.95, and 0.76–0.97, respectively. Experiments showed that the method presented in this study was the best. The algorithm segmented the tree canopies relatively completely, with clear boundaries and fewer broken canopies. The PCS algorithm produced many undersegments, and while the LM algorithm could detect small canopies, its segmentation accuracy was not high. In high-density data, compared with the LS algorithm, the r, p, and f values ​​of this study were improved by 0.02, 0.04, and 0.03, respectively. The segmentation accuracy of ITCs varied significantly across different plot densities. In low-density areas, the differences in r, p, and f among the four methods were minimal, but PCS and the algorithm proposed in this study showed slightly higher segmentation accuracy than LM and LS. In high-density areas, segmentation was more difficult due to the presence of some saplings closely adjacent to large trees, especially for the PCS method, which performed poorly in dense forests. Because the algorithm in this study refined the edges during two-dimensional segmentation, it achieved the best canopy segmentation accuracy for dense canopies.

[0206] Table 2:

[0207]

[0208] Since this algorithm involves the optimized segmentation of derived seedlings, precision reflects the segmentation accuracy of the trees. Therefore, it is clear that the precision improvement of this study is the greatest, followed by the f-score which considers both r and precision. The recall improvement is the smallest, as it is related to the number of oversegmented samples.

[0209] 2.3 Accuracy of tree parameters:

[0210] Based on the plot type, linear relationships were established between the extracted tree height and the measured tree height, and between the extracted tree crown and the measured tree crown. The R-squared values ​​for these relationships were calculated. 2 The calculation results for RMSE are shown in Table 3.

[0211] Table 3:

[0212]

[0213] LM: Local Maximum Algorithm; PCS: Point Cloud Segmentation Algorithm; LS: Seed Point Stacking Algorithm; OURS: Algorithm of this study.

[0214] Accuracy of tree height parameters:

[0215] Experimental results show that different segmentation methods exhibit significant differences in segmentation accuracy among plots of three densities, but the difference between the segmentation accuracy and the tree height extraction accuracy calculated after actual matching is small. 2 The difference between the maximum and minimum values ​​is within 0.06, and the difference in RMSE is less than 0.12m. These results indicate that the tree height extraction results of the four methods are stable, and the four methods can extract individual tree height parameters involving UAV-LiDAR point clouds.

[0216] Accuracy of crown width parameter:

[0217] Figure 10 This shows the linear regression of willow crown width estimates and measurements for three samples using different methods, based on R... 2 RMSE was used to verify the segmentation quality. Although actual location measurements were not possible due to practical limitations, the canopy of most experimental plots had a distinct top, yet the lower canopy stands were still clearly visible. Information on individual trees obtained through visualization and manual point cloud measurements, while having limited accuracy, can still serve as a reference value for ITCs results. Linear regression results for canopy width at three forest densities are shown below. Figure 10 As shown. From R 2 In terms of values, there is a strong correlation between the estimated values ​​calculated by the segmentation algorithm.

[0218] Figure 10 Scatter plots of measured crown width and extracted crown width obtained by four different methods. Figure 10In the diagram, (a), (b), (c), and (d) represent low-density plots, (e), (f), (g), and (h) represent medium-density plots, and (i), (j), (k), and (l) represent high-density plots. (a), (e), and (i) represent local maximum algorithms, (b), (f), and (j) represent point cloud segmentation algorithms, (c), (g), and (k) represent seed point stacking algorithms, and (d), (h), and (l) represent the algorithms used in this study.

[0219] In low-density plots, the tree canopy distribution is relatively uniform, and the trees are relatively young and have consistent growth, making them easy to segment and measure. The R-value of the tree canopy under four segmentation methods was determined. 2 Similar to RMSE, in medium-density plots, trees are relatively evenly distributed. However, due to the long-term and dense planting, some trees have their growth space occupied by surrounding trees, resulting in uneven canopy sizes. In high-density plots, under uneven and dense planting conditions, smaller trees emerge, making the canopies more compact and difficult to segment. The ITCs segmentation accuracy of LM, PCS, and LM algorithms is affected by the environment to varying degrees, with PCS performing the worst. The algorithm in this study performs the best, with R... 2 The value was 0.93 and the RMSE was 0.42m. This is mainly because the willow canopy properties were taken into account, which reduced the impact of the derivative tree on ITCs segmentation to a certain extent.

[0220] 3. Conclusion:

[0221] In recent years, numerous methods for single-tree segmentation using UAV laser point clouds have emerged, but segmentation of single trees with high canopy closure remains challenging. This study proposes a high-precision ITCs segmentation algorithm by combining the characteristics of UAV orthophotos and point cloud data. Based on canopy feature analysis, this method improves the segmentation accuracy of canopy edges by using an improved superpixel algorithm and energy function. It further refines the point cloud segmentation by using the local maximum method to find profile seed points, thus improving accuracy. The detection of each mask promotes the development of canopy segmentation towards a better solution. Validation was performed on forests of three density types using manual measurements. The canopy detection rate R ≥ 0.81, while the coefficient of determination R for canopy width estimation was also high. 2 ≥0.93, mean error RMSE≤0.42m, coefficient of determination R for height estimation 2 The mean square error (RMSE) was ≤0.36m. In low-density plots, the four methods showed little difference in performance. In high-density plots with high canopy closure, the LM and LS algorithms detected more tree canopies, but their canopy segmentation accuracy was poor. The PCS method is suitable for plots with uniform canopy size, therefore it performed the worst in high-canopy-closure plots. The method in this study has a high detection rate and superior segmentation of canopy boundaries. Furthermore, the projection of 2D segmentation onto 3D reduces the time and complexity of point cloud segmentation to some extent.

[0222] The scope of protection of this invention includes, but is not limited to, the above embodiments. The scope of protection of this invention is defined by the claims. Any substitutions, modifications, or improvements to this technology that are easily conceived by those skilled in the art fall within the scope of protection of this invention.

Claims

1. A multi-level dense canopy segmentation method based on UAV imagery and laser point clouds, characterized in that: include: Step 1: Collect stand point cloud data and stand orthophoto data; Step 2: Filter the forest stand orthophoto using Voronoi division, and filter the forest stand point cloud data using bilateral filtering; Step 3: Normalize the filtered point cloud and remove the ground point cloud using the planar grid method; Step 4: Segment the filtered orthophoto using an improved superpixel algorithm and energy function; Step 5: Extract the mask coordinates after segmentation, align the mask coordinates of the 2D image with the normalized point cloud coordinates to complete the 3D segmentation, and obtain multiple initial segmented regions; Step 6: Use a kernel density estimation algorithm based on the major axis to determine the number of crown vertices in each initial segmentation region; Step 7: Based on the number of tree canopy vertices and using a Gaussian mixture model, perform point cloud segmentation to obtain the tree canopy separation results; Step 4 specifically includes: Step 4.1: Extract the color distance component, spatial distance component, and texture similarity component from the filtered orthophoto image; Step 4.2: Based on the color distance component, spatial distance component, and texture similarity component, and using an improved superpixel algorithm, the orthophoto is segmented into multiple region blocks; Step 4.3: Merge region blocks based on energy function; Step 4.3.1: Calculate the energy value of each region block, where the region representative point of the j-th region block... The local average energy value is simply referred to as the energy value of region j. The energy value of region j is: （1）； in, Represents the number of pixels within region block j. Energy value; Step 4.3.2: Select cluster centers: （2）； （3）； in, It is the energy value of region block j with region block j as the cluster center, where the representative point of the region with region block j as the cluster center is... The representative point of the cluster center region is denoted as [the point of cluster center region]. , It is the average energy value of all regions. The distance between the representative point of a region block in the same cluster and the representative point of the cluster center region. distance, It is all Expectations This represents the distance between the representative point of a region block within the same cluster and the representative point of the cluster center. The farthest distance; The region block j that satisfies formulas (2) and (3) is the cluster center; Conditions for initially determining that they belong to the same cluster: The region representative point of the region block with the highest energy value The search is performed in eight directions around the center, with the direction furthest from the center as the reference. The region located at the furthest direction is then searched. The local average energy value of the representative points in the region decreases to a certain value. Or when the region block is located at the farthest position When the local average energy value of the representative point in the region no longer continues to decrease, then the region block is considered as such. Regional representative point and regional representative point of the highest energy value The distance is the radius, and the representative point of the region is... Draw a circle with the center as the center, and the regions within the circle belong to the same cluster; Optimized clusters: Using region block j as the cluster center, when region block Energy similarity value with region block j Greater than or equal to the farthest region block in the same cluster Energy similarity value with region block j At that time, the area block It belongs to the cluster centered on region block j; Using region block j as the cluster center, when region block Energy similarity value with region block j Smaller than the farthest region block in the same cluster Energy similarity value with region block j At that time, the area block It does not belong to the cluster centered on region block j; Where region block j and region block Energy similarity value for: ; Represents a region block Energy value; region block Energy similarity value with region block j for: ; Represents a region block Energy value, It is the set of all regions that belong to the same cluster as the cluster center region j; Step 4.3.3: Calculate the energy value of the remaining regions that do not belong to a cluster. Then, select cluster centers in the remaining regions according to the method in Step 4.3.2, determine the regions that belong to the same cluster, and optimize the cluster. Step 4.3.4: Select cluster centers for the remaining regions that do not belong to a cluster, determine the regions that belong to the same cluster, and optimize the cluster, until all regions belong to the corresponding cluster. Step 4.3.5: If a region block belongs to two or more different clusters, calculate and compare the energy similarity values ​​of the region block with the cluster center region blocks in the different clusters to which it belongs. The region block belongs to the cluster to which the cluster center region block with the highest energy similarity value belongs. Step 4.3.6: Merge regions belonging to the same cluster.

2. The multi-layered dense canopy segmentation method based on UAV imagery and laser point clouds according to claim 1, characterized in that: In step 6, the formula for the kernel density estimation algorithm is: （4）； （5）； In formula (4), The number of points in each initial segment. This is a smoothing parameter, called bandwidth. For kernel function, Represents sample points, Indicates the observation point.

3. The multi-level dense canopy segmentation method based on UAV imagery and laser point clouds according to claim 1, characterized in that: In step 7, a Gaussian mixture model is used for point cloud segmentation, and its formula is: （6）； In formula (6), This represents the number of canopy points estimated by kernel density, and X represents the points projected onto the profile. is the proportionality coefficient, representing the prior probability of each mixture component. Let represent the initial parameters of the Gaussian distribution, and let represent the mean and variance of the Gaussian function, respectively.

Images