A forest gap fraction prediction method based on multi-source data fusion

Abstract

This invention relates to the field of forestry remote sensing and forest structure parameter inversion technology, specifically a forest porosity prediction method based on multi-source data fusion, comprising the following steps: S1: Collecting UAV-borne lidar point cloud data and concurrent multispectral image data of the target forest area, and constructing a digital elevation model (DEM) based on ground points in the UAV-borne lidar point cloud data; simultaneously, setting up typical sample plots within the target forest area, and collecting ground-based lidar TLS point cloud data within the typical sample plot area; S2: Performing spatial registration, denoising, ground point filtering, and elevation normalization processing on the UAV-borne lidar point cloud data, multispectral image data, and DEM to construct a standardized forest structure dataset. This invention uses UAV-borne lidar point clouds, multispectral image data, and DEM as the main data sources, and introduces measured TLS porosity of typical sample plots as supervised learning labels, constructing a complete technical closed loop from sample plot determination to regional prediction.

Description

Technical Field

[0001] This invention belongs to the field of forestry remote sensing and forest structure parameter inversion technology, specifically involving a method for predicting forest porosity based on multi-source data fusion. Background Technology

[0002] Forest porosity is an important parameter characterizing the openness, spatial heterogeneity, and structural state of a forest canopy, and it plays a significant role in forest structure assessment, regeneration monitoring, and management regulation. Existing methods for obtaining forest porosity mainly include ground quadrat surveys, canopy porosity extraction based on two-dimensional remote sensing images, and canopy structure inversion based on lidar data.

[0003] Existing methods have the following shortcomings: First, although ground survey methods can obtain local sample point information, their coverage is limited and their efficiency is low, making it difficult to meet the needs of continuous estimation of porosity at the regional scale. Second, methods based on two-dimensional remote sensing images mainly rely on canopy projection information, which is difficult to effectively characterize the hierarchical structure of forest gaps in the vertical direction. Third, although single lidar data can reflect the three-dimensional structure of forest stands, if synchronous spectral information and topographic factors are lacking, the porosity estimation results are prone to deviation under complex forest stand types and undulating terrain conditions.

[0004] Furthermore, existing technologies often use lidar data primarily to generate digital elevation models or canopy height models, while the extraction of vertical structural features, horizontal connectivity features, void volume distribution features, and terrain correction features required for porosity prediction is not systematic enough. At the same time, most methods lack a model training step that uses high-precision 3D ground observation results as supervisory labels, resulting in insufficient verifiability and stability of regional-scale prediction results.

[0005] Therefore, it is necessary to propose a forest porosity prediction method based on multi-source data fusion, which uses UAV-borne lidar point cloud, multispectral image data, and digital elevation model (DEM) as the main data sources, and combines the measured porosity of ground-based lidar (TLS) in typical sample plots as training labels, so as to achieve regionalized, quantitative, and verifiable prediction of forest porosity. Summary of the Invention

[0006] The purpose of this invention is to provide a forest porosity prediction method based on multi-source data fusion to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A forest porosity prediction method based on multi-source data fusion includes the following steps: S1: Collect UAV-borne lidar point cloud data and concurrent multispectral image data of the target forest area, and construct a digital elevation model (DEM) based on the ground points in the UAV-borne lidar point cloud data; at the same time, set up typical sample plots in the target forest area and collect ground-based lidar TLS point cloud data within the typical sample plot area. S2: Spatial registration, denoising, ground point filtering, and elevation normalization are performed on UAV-borne lidar point cloud data, multispectral image data, and DEM to construct a standardized forest structure dataset; S3: Perform registration, denoising, ground normalization, and voxelization on the TLS point cloud data of typical sample plots, and calculate the reference porosity of each typical sample plot as training labels for supervised learning. S4: At the scale of plot units or grid units corresponding to typical plots, extract forest porosity prediction features from the standardized forest structure dataset. The prediction features include at least lidar structure features, multispectral features, topographic features, and porosity structure features. S5: Based on the DEM, calculate the slope, aspect and related terrain parameters, perform terrain correction on the lidar structural features, multispectral features and void structure features, and construct the corrected input feature vector; S6: Using the corrected input feature vector as input and the reference porosity of typical sample plots as supervision label, establish a forest porosity random forest regression prediction model. S7: Apply the trained random forest regression model to all grid cells to be predicted in the target forest area to obtain the forest porosity prediction value; S8: Identify areas of abnormal porosity based on the degree of deviation between the predicted forest porosity value and the preset threshold range, and output reference information for forest stand structure regulation.

[0008] Preferably, step S2 includes: Unify UAV-borne lidar point cloud data, multispectral image data, and DEM into the same coordinate reference system; Outlier noise points can be removed using statistical filtering or radius filtering methods. A ground reference surface is constructed based on the ground point classification results, and the canopy height is normalized for the point cloud. Spatial matching is used to assign band reflectance information or vegetation index information from multispectral images to corresponding plot units, raster units or point cloud units to form spatially aligned multi-source fusion feature data.

[0009] Preferably, the calculation method for the reference porosity of each typical sample plot in step S3 is as follows: Using the boundary of a typical sample plot as the horizontal range of the TLS analysis space, the TLS point cloud after ground normalization is divided into three-dimensional voxels according to the preset voxel sizes Δx, Δy, and Δz. The lower boundary of the analysis space is the normalized ground elevation 0, and the upper boundary is the upper limit of the canopy determined by the percentile value of the height of the normalized point cloud in the sample plot. For the t-th individual element in the i-th typical sample plot, its occupancy state is defined as: When the t-th genus contains at least one vegetation point; When the t-th individual element does not contain any vegetation points; Will If a voxel is defined as a void voxel, then the number of void voxels in the i-th typical plot is... for:

[0010] in, To analyze the total number of prime numbers in the spatial analysis of the i-th typical sample plot; Reference porosity of the i-th typical sample plot for:

[0011] in, Indicates the TLS reference porosity of the i-th typical sample site, and the above... This serves as a supervised learning label corresponding to the typical sample plot.

[0012] Preferably, the lidar structural features in step S4 include at least the mean, maximum, standard deviation, coefficient of variation, quantiles at different heights, canopy coverage, point cloud penetration features, and statistical features of laser reflection intensity.

[0013] Preferably, the forest porosity prediction features in step S4 further include: Vertical void profile, which is a sequence of the proportion of void voxels to the total voxels of the layer at different heights; Horizontal void connectivity refers to the number, area, or connectivity ratio of connected patches formed by adjacent void units within the same height layer or the same grid layer. The void volume distribution density is the proportion of void voxels per unit volume. Multispectral features, including at least the reflectance of each multispectral band and the vegetation index constructed therefrom, are used as input features of the forest porosity prediction model to supplement the characterization of canopy cover, vegetation growth status and canopy heterogeneity. Topographic features include at least slope, aspect, relative elevation difference, and topographic relief.

[0014] Preferably, the forest gap prediction model in step S4 is a random forest regression model; cross-validation is used to optimize parameters during model training, and one or more of the following indicators are used to evaluate model performance: coefficient of determination R², root mean square error RMSE, and mean absolute error MAE.

[0015] Preferably, the terrain correction in step S5 includes: Based on DEM, calculate topographic parameters such as slope, aspect, relative elevation difference, and topographic relief of each sample plot or grid unit; Terrain normalization processing is performed on the point cloud density characteristics, penetration characteristics, and intensity characteristics in the structural features of lidar; Topographic illumination correction is performed on band reflectance and vegetation index in multispectral features; Based on the corrected point cloud occupancy state, the vertical void profile, horizontal void connectivity, and void volume distribution density are recalculated to form the corrected input feature vector. The corrected input feature vector is used as input to the random forest regression model, rather than as a post-hoc correction to the forest gap prediction value output by the model.

[0016] Preferably, the method for identifying abnormal porosity regions in step S8 is as follows: Compare the predicted forest porosity values ​​with a preset threshold range; When the porosity of a certain plot unit or grid unit is lower than the target lower limit, it is judged as an area of ​​excessive forest density; When the porosity of a certain plot unit or grid unit is higher than the target upper limit, it is determined to be an area of ​​excessive forest stand; Areas with excessively dense or sparse forest stands were identified as areas with abnormal porosity.

[0017] Preferably, the forest stand structure regulation reference information includes: The location, area, and degree of low porosity of areas with excessively dense forest stands; The location, area, and degree of high porosity of areas with excessively sparse forest stands; Suggestions for thinning, selective cutting, or replanting based on the results of identifying areas with abnormal porosity.

[0018] Preferably, the method further includes: During subsequent monitoring periods, the target forest area will be re-measured to obtain new UAV-borne lidar point cloud data and / or typical sample plot TLS point cloud data. The actual porosity obtained from the remeasurement is compared with the predicted porosity to analyze the prediction deviation. The forest porosity prediction model is updated based on the prediction bias to improve the prediction accuracy of the model under different forest stand types or different terrain conditions.

[0019] Compared with the prior art, the beneficial effects of the present invention are: (1) This invention uses UAV-borne lidar point cloud, multispectral image data and DEM as the main data sources, and introduces the measured porosity of typical sample plots as a supervised learning label to construct a complete technical closed loop from sample plot determination to regional prediction, thereby improving the verifiability and reliability of regional estimation of forest porosity.

[0020] (2) The present invention comprehensively extracts lidar structural features, multispectral features and topographic features at the scale of plot unit or grid unit, which can simultaneously characterize the vertical structure, horizontal heterogeneity and topographic differences of forest stand, which is beneficial to improving the stability and applicability of porosity prediction under complex forest stand conditions.

[0021] (3) The present invention uses a random forest regression model to establish a forest gap prediction relationship, and combines cross-validation and accuracy evaluation index to constrain the model performance, which is conducive to improving the repeatability of the model training process and the robustness of the prediction results.

[0022] (4) Based on the digital elevation model (DEM), the present invention performs terrain correction on the structural features, multispectral features and void structure features of lidar, and performs forest void prediction based on the corrected input feature vector, thereby reducing the systematic deviation caused by slope, aspect, local incident angle and light difference under complex terrain conditions, and improving the stability, adaptability and consistency of void prediction results.

[0023] (5) This invention can identify areas with abnormal porosity and output reference information for forest stand structure regulation, providing a quantitative basis for forest management and structural optimization. Attached Figure Description

[0024] Figure 1 This is a flowchart illustrating the steps of a forest porosity prediction method based on multi-source data fusion according to the present invention. Detailed Implementation

[0025] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0026] See Figure 1This embodiment provides a forest porosity prediction method based on multi-source data fusion. The method uses UAV-borne lidar point cloud data, multispectral image data, and a digital elevation model (DEM) as the main data sources, and uses reference porosity obtained from ground-based lidar TLS point cloud computing of typical sample plots as a supervision label. The method includes the following steps: S1: Collect UAV-borne lidar point cloud data and concurrent multispectral image data of the target forest area, and construct a digital elevation model (DEM) based on the ground points in the UAV-borne lidar point cloud data; at the same time, set up typical sample plots in the target forest area and collect ground-based lidar TLS point cloud data within the typical sample plot area. UAV-borne lidar point cloud data is acquired through a lidar system mounted on the UAV platform. The lidar system emits laser pulses towards the target forest area and receives the echo signals, obtaining three-dimensional point cloud data of targets such as the canopy, branches, understory vegetation, and ground surface, and recording the spatial coordinates and laser reflection intensity information of each point. UAV-borne lidar point cloud data is used to characterize the three-dimensional structural features of forest vegetation and the ground surface.

[0027] Multispectral imagery data was acquired by a multispectral sensor simultaneously mounted on an unmanned aerial vehicle (UAV) platform. The multispectral sensor can acquire imagery data in red, green, blue, and near-infrared bands. The multispectral imagery data is primarily used to extract spectral features characterizing vegetation cover and growth status, including reflectance in each band and vegetation indices constructed from them. It should be noted that in this invention, multispectral imagery data is primarily used as a source of spectral features, not as the main data source for the three-dimensional structure.

[0028] The DEM is preferably generated by interpolation of ground points in the UAV-borne lidar point cloud. It is used to characterize the surface elevation information and serves as the basis for subsequent point cloud elevation normalization and terrain factor extraction.

[0029] To obtain highly reliable reference values ​​required for supervised learning, several typical sample plots were established within the target forest area, and ground-based lidar TLS point cloud data were collected within these typical sample plots. Typical sample plots were selected based on differences in forest stand type, density level, slope position, and canopy structure. TLS point clouds were used to accurately characterize the internal canopy structure at the sample plot scale, and reference porosity was further calculated for each typical sample plot. This reference porosity served as the training label for the forest porosity prediction model. It should be noted that TLS in this invention is used as a means of sample plot identification and label acquisition, rather than as the primary data source for regional-scale porosity prediction.

[0030] The target forest area can be a mixed coniferous and broad-leaved forest in a mountainous and hilly region. A UAV-borne lidar system will be used to acquire lidar point cloud data within the forest area, and a synchronous multispectral sensor will be used to acquire multispectral image data of the corresponding area. A DEM will be generated based on the lidar ground points. Simultaneously, several representative typical plots within the forest area will be selected for TLS scanning. The plot area can be set according to the research object and flight resolution, such as 20 m × 20 m, 25 m × 25 m, or other suitable scales.

[0031] In the data preprocessing stage, the UAV-borne LiDAR point cloud data, multispectral image data, and DEM are first uniformly registered in coordinates. Registration can be based on the flight platform's positioning and attitude information, ground control points, or corresponding feature points to ensure that all types of data correspond under the same spatial reference system. Subsequently, the UAV-borne LiDAR point cloud data undergoes denoising processing to remove isolated noise points and outliers. Statistical filtering or radius filtering methods are preferred for denoising.

[0032] S2: Spatial registration, denoising, ground point filtering, and elevation normalization are performed on UAV-borne lidar point cloud data, multispectral image data, and DEM to construct a standardized forest structure dataset; Unify UAV-borne lidar point cloud data, multispectral image data, and DEM into the same coordinate reference system; Outlier noise points can be removed using statistical filtering or radius filtering methods. A ground reference surface is constructed based on the ground point classification results, and the canopy height is normalized for the point cloud. Spatial matching is used to assign band reflectance information or vegetation index information from multispectral images to corresponding plot units, raster units or point cloud units to form spatially aligned multi-source fusion feature data.

[0033] The statistical filtering method is implemented as follows: for each point in the point cloud, search for several nearest neighbor points within its neighborhood, calculate the average distance from the point to its neighbors, and statistically analyze the distribution of the average neighborhood distances of all points. When the average neighborhood distance of a point is greater than the mean of the entire distribution plus a preset multiple of the standard deviation, it is identified as an outlier and deleted. The number of nearest neighbor points, the search radius, and the standard deviation multiple can be adjusted according to the point cloud density and the complexity of the forest stand.

[0034] After denoising, ground point filtering is performed on the UAV-borne LiDAR point cloud, and elevation normalization is applied to the point cloud using a DEM to obtain the normalized altitude of each point relative to the ground surface. The formula for calculating the normalized altitude is:

[0035] in, The normalized height of the point, This represents the original elevation value of the target point in the point cloud. This refers to the surface elevation value provided by the DEM at the corresponding location.

[0036] After the above processing, the original point cloud elevation can be converted into structural parameters that reflect the relative height of trees to the ground surface, thereby eliminating the direct impact of topographic relief on canopy structure analysis.

[0037] For multispectral image data, band reflectance and vegetation index features can be further extracted. In one embodiment, the Normalized Difference Vegetation Index (NDVI) can be calculated using the following formula:

[0038] in, For near-infrared reflectivity, This represents the reflectance in the red light band. NDVI and other multispectral features can be spatially matched with lidar structural features at the scale of plot units, grid units, or point cloud units for feature construction of subsequent forest porosity prediction models. The vegetation index threshold can be adjusted according to forest stand type, sensor parameters, and the collection season, rather than being limited to a single fixed value.

[0039] During the data fusion phase, the registered and preprocessed UAV-borne lidar point cloud data, multispectral image features, and DEM-derived topographic factors are spatially mapped to form a standardized forest structure dataset. This standardized forest structure dataset can use plot units, raster units, or point cloud units as basic organizational units and includes at least the following information: spatial location, normalized altitude, laser reflectance intensity, multispectral reflectance or vegetation index, and topographic factors such as slope, aspect, and relative elevation difference extracted from the DEM.

[0040] The construction of the standardized forest structure dataset does not require all data sources to correspond one-to-one. Instead, it requires effective matching of structural features, spectral features, and topographic features under a unified spatial reference, thereby providing a consistent data foundation for subsequent typical sample plot reference porosity calculation, feature extraction, and training of forest porosity prediction models.

[0041] S3: Perform registration, denoising, ground normalization, and voxelization on the TLS point cloud data of typical sample plots, and calculate the reference porosity of each typical sample plot as training labels for supervised learning. Typical sample plot TLS point clouds are used only for calculating reference porosity labels required for supervised learning, and are not used as the main data source for regional prediction of the study area. To ensure a one-to-one correspondence between TLS labels and UAV-borne LiDAR point clouds, multispectral image data, and DEM features, the boundaries of typical sample plots adopt a planar range consistent with the plot layout as the analysis boundary, and are unified with the UAV-borne data under the same coordinate reference system.

[0042] For the i-th typical sample plot, the TLS point cloud is first subjected to multi-station stitching and registration, outlier removal, ground point identification, and elevation normalization. Then, the sample plot's planar boundary is used as the horizontal range, the normalized ground elevation 0 is used as the lower boundary of the analysis space, and the canopy's upper limit height, determined by the percentile value of the normalized point cloud height, is used as the upper boundary of the analysis space. Next, the sample plot's analysis space is divided into three dimensions using preset voxel sizes Δx, Δy, and Δz. In this embodiment, the voxel side length can be set to Δx = Δy = Δz = 0.5m.

[0043] For the t-th individual element in the i-th typical sample plot, its occupancy state is defined as: When the t-th genus contains at least one vegetation point; When the t-th individual element does not contain any vegetation points; Will If a voxel is defined as a void voxel, then the number of void voxels in the i-th typical plot is... for:

[0044] in, To analyze the total number of prime numbers in the spatial analysis of the i-th typical sample plot; Reference porosity of the i-th typical sample plot for:

[0045] in, Indicates the TLS reference porosity of the i-th typical sample site, and... This serves as a supervised learning label corresponding to the typical sample plot.

[0046] S4: At the scale of plot units or grid units corresponding to typical plots, extract forest porosity prediction features from the standardized forest structure dataset. The prediction features include at least lidar structure features, multispectral features, topographic features, and porosity structure features. The structural characteristics of lidar include at least the mean, maximum, standard deviation, coefficient of variation, quantiles at different heights, canopy coverage, point cloud penetration characteristics, and statistical characteristics of laser reflection intensity.

[0047] The forest porosity prediction features in step S4 also include: Vertical void profile, which is a sequence of the proportion of void voxels to the total voxels of the layer at different heights; Horizontal void connectivity refers to the number, area, or connectivity ratio of connected patches formed by adjacent void units within the same height layer or the same grid layer. The void volume distribution density is the proportion of void voxels per unit volume. Multispectral features, which include at least the reflectance of each multispectral band and the vegetation index constructed from them, serve as input features for the forest porosity prediction model and are used to supplement the characterization of canopy cover, vegetation growth status and canopy heterogeneity. Topographic features include at least slope, aspect, relative elevation difference, and topographic relief.

[0048] The calculation steps for horizontal porosity and vertical porosity profile are as follows: Feature extraction was performed on the standardized forest structure dataset after coordinate unification, elevation normalization, and vegetation point selection. The extracted features included at least lidar structure features, multispectral features, topographic features, and void structure features, among which void structure features included vertical void profiles, horizontal void connectivity, and void volume distribution density.

[0049] A standardized forest structure dataset must include at least the three-dimensional coordinates of points, normalized height, echo intensity, multispectral reflectance or vegetation index, and DEM-derived topographic factors. The calculation object is the sample plot area or the area of ​​the grid cell to be predicted, and calculations are performed at fixed height intervals along the vertical direction. Δh The point cloud is sliced ​​into layers to obtain multiple horizontal slice layers.

[0050] In this embodiment, a mixed coniferous and broad-leaved forest in a mountainous and hilly area is used as an example. The vertical stratification interval is set to Δh = 0.5 m, starting from the normalized ground elevation z = 0 m to the highest point of the forest canopy H. max Up to 25.2 m, a total of K=51 horizontal slice layers were obtained. The k-th horizontal slice layer is denoted as:

[0051] Each horizontal slice layer Sk is projected onto the XY plane and divided into regular grid cells with side length r. In this embodiment, r = 0.1m. For the (u,v)th grid cell in the kth layer, the vegetation occupancy indicator variable is defined as:

[0052] Let Nk be the total number of grid cells in the k-th layer, then the horizontal coverage Ck of this layer is:

[0053] Accordingly, the horizontal porosity Gk of the k-th layer is defined as:

[0054] Where Gk represents the proportion of voids in the projection plane of that height layer. Arranging the horizontal porosity of each slice layer in order of height yields the basic sequence of vertical void profiles:

[0055] Vertical void profiles are used to characterize the void variation patterns from the understory to the canopy in a forest stand.

[0056] To characterize the vertical connectivity of voids between adjacent height layers, adjacent slice layers S... k and S k+1 The corresponding grids are compared. A vertical gap continuity index C is defined between layer k and layer (k+1). v,k for:

[0057] Among them, (1−B k (u,v))=1 indicates that the grid cell corresponding to the k-th layer is a gap grid cell; the numerator represents the number of gap grid cells that maintain vertical continuity between adjacent layers; the denominator represents the total number of gap grid cells in the k-th layer. C v,k The value range is [0,1]. The larger the value, the higher the degree of vertical void penetration.

[0058] Arrange the vertical void continuity indices of each layer to obtain:

[0059] The steps for calculating the connectivity of horizontal gaps are as follows: Perform connectivity analysis on the void grid in each horizontal slice layer. This will satisfy B... k A grid with (u,v)=0 is defined as a gapped grid, and a gapped connected domain is constructed based on the eight-neighbor rule. Let M be the total number of grids identified at the k-th layer. k There are j-th interstitial connected components with gaps, and the area of ​​the j-th interstitial connected component is a. k,j The perimeter is p k,j Then the total void area of ​​the kth layer is:

[0060] To characterize the connectivity of horizontal voids, the maximum connected patch proportion H in the k-th layer is defined. k for:

[0061] Among them, H k The larger the value, the more concentrated the voids in the layer are, forming a larger connected opening area.

[0062] Furthermore, to describe the complexity of the gap morphology, the shape index of the j-th gap connected region is defined as:

[0063] Then the area-weighted average shape index S of the k-th layer k for:

[0064] Among them, a k,j It is obtained by multiplying the number of connected graticles by the area of ​​a single graticle. Simultaneously, the average connected region area can be calculated:

[0065] Therefore, the horizontal void connectivity feature of the k-th layer can be expressed as:

[0066] The horizontal void connectivity features of each height layer are combined in height order to form a horizontal void connectivity feature sequence.

[0067] S5: Based on the DEM, calculate the slope, aspect and related terrain parameters, perform terrain correction on the lidar structural features, multispectral features and void structure features, and construct the corrected input feature vector; The forest porosity prediction model in step S4 is a random forest regression model. During the model training process, cross-validation is used to optimize parameters, and one or more of the following indicators are used to evaluate the model performance: coefficient of determination R², root mean square error RMSE, and mean absolute error MAE.

[0068] The terrain correction in step S5 includes: Based on DEM, calculate topographic parameters such as slope, aspect, relative elevation difference, and topographic relief of each sample plot or grid unit; Terrain normalization processing is performed on the point cloud density characteristics, penetration characteristics, and intensity characteristics in the structural features of lidar; Topographic illumination correction is performed on band reflectance and vegetation index in multispectral features; Based on the corrected point cloud occupancy state, the vertical void profile, horizontal void connectivity, and void volume distribution density are recalculated to form the corrected input feature vector. The corrected input feature vector is used as the input to the random forest regression model, rather than as a post-hoc correction to the forest gap prediction value output by the model.

[0069] The steps for calculating the void volume distribution density are as follows: In three-dimensional space, the sample plot or the grid cell to be predicted is divided into regular voxel cells. Let the side lengths of the voxels be Δx, Δy, and Δz, respectively; in this embodiment, Δx = Δy = Δz = 0.5 m. Since all voxels have the same size, the volume fraction is characterized by the proportion of voxels. For any voxel V...t Its occupied state is defined as:

[0070] Within a larger-scale three-dimensional statistical unit, such as a volume block Qm of 2 m × 2 m × 2 m, the void volume distribution density Dm is defined as:

[0071] Where, N m For statistical unit Q m The total number of voxels contained in D. m The value range is [0,1], representing the proportion of void voxels within this statistical unit. The D of all statistical units... m Together they constitute the feature set of void volume distribution density.

[0072] The steps for constructing the feature vector of the void structure are as follows: In one embodiment of the present invention, the input variables for a forest porosity prediction model are constructed based on the lidar structural features, multispectral features, terrain features, and porosity structural features extracted in the aforementioned steps. The porosity structural features include vertical porosity profile features, horizontal porosity connectivity features, and porosity volume distribution density features. To meet the requirements of the random forest regression model for fixed-length input variables, statistics are further extracted from various sequence or stratified features to form feature vectors at the plot or grid unit scale.

[0073] Therefore, the input feature vector of the i-th sample plot or grid cell can be expressed as:

[0074] Among them, f LiDAR,i At least including: point cloud density features ρ i =N i / A i , where N i Let A be the number of laser points within the i-th sample plot or grid cell. i Its horizontal projected area; point cloud penetration feature Pi=N below,i / N i , where N below,i The number of echo points below a preset height threshold or entering the canopy; the laser intensity characteristic includes at least the average echo intensity I. mean,i =(1 / Ni)∑I p Sum of standard deviations I std,i f MS,i At least include reflectivity R for each band b,i and the vegetation indices constructed from them, among which the characteristic vegetation index can be the Normalized Difference Vegetation Index (NDVI). i =(R NIR,i-R Red,i ) / (R NIR,i +R Red,i f DEM,i It should include at least elevation, slope, aspect, relative elevation difference, topographic relief, and curvature. Gap,i It should include at least the mean, standard deviation, extreme values ​​and corresponding height positions of the vertical void profile, the statistics of horizontal void connectivity, and the statistics of void volume distribution density.

[0075] The lidar terrain correction steps are as follows: To reduce the impact of complex terrain on the sampling density of UAV-borne lidar point clouds, the spectral response of multispectral images, and the characterization results of void structures, terrain correction is performed on the input features.

[0076] Taking any pixel or grid cell j within the study area as the object, its local topographic parameters are calculated based on the DEM. Let the DEM elevation surface be Z(x,y), then the first-order topographic gradient at that location is denoted as:

[0077] Then the local slope α j and slope β j Calculated separately as follows:

[0078] Wherein, the slope angle α is defined as the angle between the surface normal vector and the vertical direction, α j Indicates the local slope angle, β j This represents the local slope aspect angle. Slope and aspect can be calculated using the 3×3 neighborhood finite difference method.

[0079] For the UAV-borne lidar point cloud, the laser zenith angle θ at position j for each laser pulse is obtained based on the flight path and scanning parameters. l,j and scanning azimuth angle ϕ l,j Further construct the local surface normal vector n j and the incident direction vector l of the laser j :

[0080] Therefore, the local effective incident angle γ is calculated. j :

[0081] Where, γ j This represents the angle between the laser beam and the local surface normal. When the terrain undulation increases or the scanning direction is inconsistent with the slope aspect, γ... j As the density increases, the local sampling density and penetration performance will change.

[0082] Terrain normalization processing is performed on the point cloud density features, penetration features, and intensity features in the lidar structural features.

[0083] To normalize point cloud densities under different terrain conditions to a horizontal reference condition, a laser point cloud density terrain correction factor K is defined. L,j for:

[0084] Where ε is a very small positive number to prevent the denominator from being zero. If the vegetation point density observed within a statistical unit at position j is... Then the corrected vegetation point density is:

[0085] Accordingly, for the (u,v)th grid cell of the kth horizontal slice layer in the calculation of void structure characteristics, let its number of observed vegetation points be... The corrected number of vegetation points is: .

[0086] The steps for multispectral image topographic correction are as follows: For multispectral image data, the solar incidence angle ij at position j is calculated based on the solar zenith angle θs and solar azimuth angle ϕs at the time of acquisition, combined with the local slope αj and aspect βj:

[0087] For the b-th spectral band, the C-correction model is used to correct for topographic irradiance differences. If the original reflectance is... Then the corrected reflectivity for:

[0088] Where Cb is an empirical parameter for the b-th band, which can be obtained from the reflectivity and cosine of that band. j The linear regression relationship was determined.

[0089] Based on the corrected band reflectance, the vegetation index is recalculated. Taking the Normalized Difference Vegetation Index (NDVI) as an example, its corrected expression is:

[0090] in, To correct reflectivity in the near-infrared band, To correct reflectivity in the red light band. Based on Alternatively, other corrected vegetation indices can be used to re-distinguish between vegetation and non-vegetation, thereby reducing the impact of light differences between shady and sunny slopes on the accuracy of vegetation identification.

[0091] Topographic illumination correction is performed on the band reflectance and vegetation index in the multispectral features.

[0092] The steps for feature reconstruction after terrain correction are as follows: Based on the above correction results, the lidar structural features, multispectral features, and void structure features are reconstructed respectively. For the j-th grid cell to be predicted, its terrain-corrected input feature vector is denoted as:

[0093] Among them, f' LiDAR,j f′ represents the structural features of the lidar after terrain correction. MS,j This represents the multispectral characteristics after illumination-topography correction, f DEM,j f′ represents the terrain features extracted from the DEM. Gap,j This represents the void structure features, such as the vertical void profile, horizontal void connectivity, and void volume distribution density, obtained by recalculating based on the corrected point cloud occupancy state.

[0094] Based on the corrected number of vegetation points, re-determine whether the grid cell is a vegetation-occupied grid cell:

[0095] Where, τ n The threshold for determining vegetation occupation is set. Based on this, the horizontal cover, horizontal porosity, vertical porosity continuity, and horizontal porosity connectivity characteristics of each height layer are recalculated.

[0096] Corrected horizontal coverage of the kth layer and horizontal porosity They are respectively:

[0097] Where Nk is the total number of grid cells in the k-th layer.

[0098] S6: Using the corrected input feature vector as input and the reference porosity of typical sample plots as supervision label, establish a forest porosity random forest regression prediction model. The training steps for a random forest regression model are as follows: In this invention, typical sample plot TLS point clouds are used to calculate reference porosity labels. Let yi be the reference porosity calculated by TLS voxelization for the i-th typical sample plot, and xi be the corresponding input feature vector. The random forest regression model can be expressed as:

[0099] in, i represents the model prediction porosity of the i-th typical sample plot, T represents the number of regression trees, and ht(xi) represents the output of the t-th regression tree to the input feature vector xi.

[0100] In practice, all typical sample plots can be divided into training and validation sets, or k-fold cross-validation can be used to train and validate the model. After model training is complete, the model accuracy is evaluated using validation set samples. The model evaluation metric should at least include the coefficient of determination R0. 2 The root mean square error (RMSE) and the mean absolute error (MAE) are calculated using the following formulas:

[0101] Where n is the number of sample plots, and yi is the TLS reference porosity of the i-th typical sample plot. i represents the model-predicted porosity of the i-th typical sample plot. This represents the average TLS reference porosity across all typical sample plots.

[0102] The supervision label yi of the i-th typical sample plot is the reference porosity obtained based on TLS point cloud computing. The training sample set consists of all typical sample plots.

[0103] Where n is the number of typical sample plots.

[0104] S7: Apply the trained random forest regression model to all grid cells to be predicted in the target forest area to obtain the forest porosity prediction value; The steps for predicting porosity after terrain correction are as follows: The terrain-corrected input feature vector x′j is input into the trained random forest regression model to obtain the predicted forest gap density at location j:

[0105] Where x′j represents the terrain-corrected input feature vector of the j-th grid cell to be predicted, and yij represents the predicted porosity value of the corresponding grid cell. This represents the random forest regression function after training is complete.

[0106] After the model training is completed, the trained model is applied to all the grid cells to be predicted in the study area to obtain regionalized forest porosity prediction results. Finally, the porosity prediction values ​​of each grid cell obtained after terrain correction are spatially stitched and mapped to generate a forest porosity distribution map of the study area.

[0107] S8: Identify areas of abnormal porosity based on the degree of deviation between the predicted forest porosity value and the preset threshold range, and output reference information for forest stand structure regulation; The method for identifying porosity anomaly regions in step S8 is as follows: Compare the predicted forest porosity values ​​with a preset threshold range; When the porosity of a certain plot unit or grid unit is lower than the target lower limit, it is judged as an area of ​​excessive forest density; When the porosity of a certain plot unit or grid unit is higher than the target upper limit, it is determined to be an area of ​​excessive forest stand; Areas with excessively dense or sparse forest stands were identified as areas with abnormal porosity.

[0108] Preferably, the reference information for forest stand structure regulation includes: The location, area, and degree of low porosity of areas with excessively dense forest stands; The location, area, and degree of high porosity of areas with excessively sparse forest stands; Suggestions for thinning, selective cutting, or replanting based on the results of identifying areas with abnormal porosity.

[0109] The steps for identifying abnormal regions are as follows: Based on the porosity prediction results at the region scale, regions with abnormal porosity can be further identified. In one embodiment, let the predicted porosity of the j-th grid cell be... j, with an upper threshold of T high The lower threshold is T low .when j>T high When, it is determined to be a region with high porosity; when j <T low When the porosity is low, the area is identified as a region with low porosity; the remaining areas are identified as regions with normal porosity.

[0110] To improve the interpretability of the results of anomaly region identification, the structural type of the anomaly region can be analyzed by combining the vertical void profile features, horizontal void connectivity features, and void volume distribution density features of the corresponding grid cells.

[0111] Based on the predicted porosity and the results of anomaly region identification, key porosity regions can be identified in three-dimensional space. A key porosity region is a set of grid cells that meet preset anomaly judgment conditions and have a continuous spatial distribution. Optionally, key porosity regions can be stored as vector surfaces, raster layers, or three-dimensional bounding boxes, and their spatial location, area, corresponding height layer range, predicted porosity, and anomaly type can be recorded for reference in subsequent forest stand structure regulation.

[0112] Based on the porosity prediction results and anomaly area identification results, the study area can be divided into zones and mapped. For each anomaly area, its spatial location, area, predicted porosity value, anomaly type, main structural features, and forest stand structure regulation reference information can be recorded, and the relevant information can be stored as a raster layer or vector layer. Optionally, the corresponding height layer range can be overlaid in three-dimensional space to form key porosity area identification results for subsequent forest stand structure analysis and management decision-making reference.

[0113] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A forest porosity prediction method based on multi-source data fusion, characterized in that, Includes the following steps: S1: Collect UAV-borne lidar point cloud data and concurrent multispectral image data of the target forest area, and construct a digital elevation model (DEM) based on the ground points in the UAV-borne lidar point cloud data; at the same time, set up typical sample plots in the target forest area and collect ground-based lidar TLS point cloud data within the typical sample plot area. S2: Spatial registration, denoising, ground point filtering, and elevation normalization are performed on UAV-borne lidar point cloud data, multispectral image data, and DEM to construct a standardized forest structure dataset; S3: Perform registration, denoising, ground normalization, and voxelization on the TLS point cloud data of typical sample plots, and calculate the reference porosity of each typical sample plot as training labels for supervised learning. S4: At the scale of plot units or grid units corresponding to typical plots, extract forest porosity prediction features from the standardized forest structure dataset. The prediction features include at least lidar structure features, multispectral features, topographic features, and porosity structure features. S5: Based on the DEM, calculate the slope, aspect and related terrain parameters, perform terrain correction on the lidar structural features, multispectral features and void structure features, and construct the corrected input feature vector; S6: Using the corrected input feature vector as input and the reference porosity of typical sample plots as supervision label, establish a forest porosity random forest regression prediction model. S7: Apply the trained random forest regression model to all grid cells to be predicted in the target forest area to obtain the forest porosity prediction value; S8: Identify areas of abnormal porosity based on the degree of deviation between the predicted forest porosity value and the preset threshold range, and output reference information for forest stand structure regulation.

2. The forest porosity prediction method based on multi-source data fusion according to claim 1, characterized in that, Step S2 includes: Unify UAV-borne lidar point cloud data, multispectral image data, and DEM into the same coordinate reference system; Outlier noise points can be removed using statistical filtering or radius filtering methods. A ground reference surface is constructed based on the ground point classification results, and the canopy height is normalized for the point cloud. Spatial matching is used to assign band reflectance information or vegetation index information from multispectral images to corresponding plot units, raster units or point cloud units to form spatially aligned multi-source fusion feature data.

3. The forest porosity prediction method based on multi-source data fusion according to claim 1, characterized in that, The calculation method for the reference porosity of each typical sample plot in step S3 is as follows: Using the boundary of a typical sample plot as the horizontal range of the TLS analysis space, the TLS point cloud after ground normalization is divided into three-dimensional voxels according to the preset voxel sizes Δx, Δy, and Δz. The lower boundary of the analysis space is the normalized ground elevation 0, and the upper boundary is the upper limit of the canopy determined by the percentile value of the height of the normalized point cloud in the sample plot. For the t-th individual element in the i-th typical sample plot, its occupancy state is defined as: When the t-th genus contains at least one vegetation point; When the t-th individual element does not contain a vegetation point; Will If a voxel is defined as a void voxel, then the number of void voxels in the i-th typical plot is... for: in, To analyze the total number of prime numbers in the spatial analysis of the i-th typical sample plot; Reference porosity of the i-th typical sample plot for: in, Indicates the TLS reference porosity of the i-th typical sample site, and the above... This serves as a supervised learning label corresponding to the typical sample plot.

4. The forest porosity prediction method based on multi-source data fusion according to claim 1, characterized in that, The lidar structural features in step S4 include at least the mean, maximum, standard deviation, coefficient of variation, quantiles at different heights, canopy coverage, point cloud penetration features, and statistical features of laser reflection intensity.

5. The forest porosity prediction method based on multi-source data fusion according to claim 4, characterized in that, The forest porosity prediction features in step S4 also include: Vertical void profile, which is a sequence of the proportion of void voxels to the total voxels of the layer at different heights; Horizontal void connectivity refers to the number, area, or connectivity ratio of connected patches formed by adjacent void units within the same height layer or the same grid layer. The void volume distribution density is the proportion of void voxels per unit volume. Multispectral features, including at least the reflectance of each multispectral band and the vegetation index constructed therefrom, are used as input features of the forest porosity prediction model to supplement the characterization of canopy cover, vegetation growth status and canopy heterogeneity. Topographic features include at least slope, aspect, relative elevation difference, and topographic relief.

6. The forest porosity prediction method based on multi-source data fusion according to claim 5, characterized in that, The forest porosity prediction model in step S4 is a random forest regression model. During the model training process, cross-validation is used to optimize parameters, and one or more of the following indicators are used to evaluate the model performance: coefficient of determination R², root mean square error RMSE, and mean absolute error MAE.

7. The forest porosity prediction method based on multi-source data fusion according to claim 1, characterized in that, The terrain correction in step S5 includes: Based on DEM, calculate topographic parameters such as slope, aspect, relative elevation difference, and topographic relief of each sample plot or grid unit; Terrain normalization processing is performed on the point cloud density characteristics, penetration characteristics, and intensity characteristics in the structural features of lidar; Topographic illumination correction is performed on band reflectance and vegetation index in multispectral features; Based on the corrected point cloud occupancy state, the vertical void profile, horizontal void connectivity, and void volume distribution density are recalculated to form the corrected input feature vector. The corrected input feature vector is used as input to the random forest regression model, rather than as a post-hoc correction to the forest gap prediction value output by the model.

8. The forest porosity prediction method based on multi-source data fusion according to claim 1, characterized in that, The method for identifying abnormal porosity regions in step S8 is as follows: Compare the predicted forest porosity values ​​with a preset threshold range; When the porosity of a certain plot unit or grid unit is lower than the target lower limit, it is judged as an area of ​​excessive forest density; When the porosity of a certain plot unit or grid unit is higher than the target upper limit, it is determined to be an area of ​​excessive forest stand; Areas with excessively dense or sparse forest stands were identified as areas with abnormal porosity.

9. A forest porosity prediction method based on multi-source data fusion according to claim 8, characterized in that, The forest stand structure regulation reference information includes: The location, area, and degree of low porosity of areas with excessively dense forest stands; The location, area, and degree of high porosity of areas with excessively sparse forest stands; Suggestions for thinning, selective cutting, or replanting based on the results of identifying areas with abnormal porosity.

10. The forest porosity prediction method based on multi-source data fusion according to claim 1, characterized in that, The method also includes: During subsequent monitoring periods, the target forest area will be re-measured to obtain new UAV-borne lidar point cloud data and / or typical sample plot TLS point cloud data. The actual porosity obtained from the remeasurement is compared with the predicted porosity to analyze the prediction deviation. The forest porosity prediction model is updated based on the prediction bias to improve the prediction accuracy of the model under different forest stand types or different terrain conditions.

Images