A multi-spectral spatial analysis algorithm for forest vegetation cover

Abstract

The present application relates to the technical field of remote sensing and artificial intelligence, and discloses a multispectral spatial analysis algorithm for forest vegetation coverage. The algorithm synchronously acquires LiDAR point cloud data and multispectral image data of a forest area, constructs a three-dimensional point cloud-spectrum fusion voxel grid after radiation calibration, geometric registration and coordinate unification; extracts multispectral reflectance features and three-dimensional geometric structure features at the voxel layer, and realizes joint modeling through cross-modal attention and graph convolution-Transformer hybrid coding; then, the aboveground biomass density field is predicted in combination with the crown effective thickness correction factor, and the aboveground biomass result of the pixel scale forest is obtained by integrating along the vertical direction. The method can improve the estimation stability and accuracy under complex canopy and cross-region conditions.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of remote sensing and artificial intelligence, specifically a multispectral spatial analysis algorithm for forest vegetation cover. Background Technology

[0002] With the improvement of multi-source remote sensing data acquisition capabilities, multispectral imagery can provide spectral response information such as vegetation chlorophyll, water content, and structure, while LiDAR point clouds can provide three-dimensional structural information such as canopy height, vertical stratification, and topography. Existing forest aboveground biomass inversion methods primarily rely on pixel-level two-dimensional features, typically concatenating LiDAR-derived structural indices with multispectral indices and then fitting the data using models such as random forests and support vector regression. While these methods can achieve a certain level of accuracy under single forest type and single sensor conditions, they are prone to problems such as insufficient feature representation, loss of spatial detail, and decreased generalization ability under conditions of significant spatial heterogeneity, severe canopy shading, multi-layered forest stands, or cross-regional observations.

[0003] For example, existing techniques propose first screening LiDAR 3D structural variables and multispectral features through correlation analysis, and then using regression models to estimate aboveground biomass. This process relies on human experience to determine features, scale, and regression form, making it difficult to characterize the strong nonlinear coupling relationship between the spectrum and 3D structure. Furthermore, canopy shading can lead to insufficient LiDAR acquisition of understory information, causing systematic bias. In addition, domain shifts caused by remote sensing imaging conditions, topography, and differences in forest type often require recalibration of the model using a large number of ground samples from the target area. Therefore, there is an urgent need for an intelligent analysis algorithm that can deeply integrate multispectral and point cloud information in voxel 3D space and possess cross-domain adaptive capabilities.

[0004] The above figures are merely illustrative representations of the process of this invention and do not limit the data source or sensor type. The voxel grid resolution, point cloud density, and number of bands in the figures can be adjusted according to actual task requirements. Data flow between modules is transmitted in vector / tensor mode; in actual implementation, sparse voxel representation can be used to reduce storage and computational overhead. Summary of the Invention

[0005] The purpose of this invention is to provide a multispectral spatial analysis algorithm for forest vegetation cover, addressing the problems of existing methods mentioned in the background, such as their inability to effectively integrate deep semantics of high-dimensional heterogeneous remote sensing data, insufficient ability to model complex canopy structures, and reliance on large amounts of ground samples. To achieve the above objective, this invention provides the following technical solution: A multispectral spatial analysis algorithm for forest vegetation cover includes the following steps: Step 1: Simultaneously acquire LiDAR point cloud data and multispectral image data of the forest area. Perform radiometric calibration and atmospheric correction on the multispectral image, coordinate unification and outlier removal on the LiDAR point cloud, and geometric registration and spatiotemporal alignment on the two types of data. Construct a 3D point cloud-spectral fusion voxel mesh under a unified geographic coordinate system. ; Step 2: In the fused voxel mesh In the process, the multispectral reflectance features of each voxel unit are extracted. With three-dimensional geometric features The three-dimensional geometric features include local normal vector entropy, rate of curvature change, vertical layer density, and echo intensity statistics; and a cross-modal attention mechanism is used to... and Interactive weight assignment is performed to obtain a high-dimensional heterogeneous feature tensor. ; The following is a relatively complete embodiment that enables the implementation of the technical solution of the present invention. It should be noted that this embodiment is used to explain the present invention and not to limit the scope of protection of the present invention; without departing from the spirit of the present invention, the order of related steps, parameter settings, and network structure can be equivalently replaced. To avoid ambiguity, several key symbols are defined as follows in this embodiment: the fused voxel mesh is denoted as... Voxel index is denoted as The voxel center coordinates are as follows: The center of the voxel height layer is denoted as Voxel multispectral characteristics are denoted as The geometric structural features of voxels are denoted as ; The fusion feature is denoted as The encoder output embedding is denoted as Voxel biomass density is denoted as ; Aboveground biomass at the pixel scale is denoted as .

[0006] Step 3: Convert the high-dimensional heterogeneous feature tensor The input is a graph convolutional-Transformer hybrid encoder, where the graph convolutional module captures the topological connectivity between neighboring voxels, and the Transformer module models long-range spatial dependencies through a self-attention mechanism. The output is a context-aware embedding vector. ; Step 4: Embedding Vectors Based on Context Awareness The continuously distributed aboveground biomass density field is predicted using voxel regression heads, and the estimation bias caused by canopy shading is suppressed by combining an effective canopy thickness correction factor. The calculation formula is as follows: in, Represents the coordinates of the voxel center Biomass density at the location Embed vectors for context awareness. and For learnable parameters, This represents the effective thickness of the canopy at the corresponding horizontal position. The attenuation coefficient is... As a stable term, This is used to ensure that the output is non-negative. Among them, A smaller value usually indicates that the low-quantile echo rise or insufficient penetration leads to stronger canopy shading. The correction factor attenuates the density prediction for this area more strongly to reduce the systematic overestimation bias caused by shading.

[0007] Step 5: Integrate the biomass density field along the vertical direction and aggregate it to the pixel scale to obtain the final forest aboveground biomass inversion result. Furthermore, unsupervised contrastive learning is used to perform domain adaptive optimization in unlabeled regions to improve the model's generalization ability under different forest types and observation conditions.

[0008] Compared with the prior art, the present invention has the following beneficial effects: (1) This invention uses a three-dimensional point cloud-spectral fusion voxel grid as a carrier to simultaneously model multispectral reflectance information and three-dimensional geometric structure information in the voxel layer. Combined with cross-modal attention and graph convolution-Transformer hybrid coding, it can effectively improve the feature expression ability in complex canopy, multi-layer forest stand and spatially heterogeneous scenarios, thereby improving the accuracy and stability of forest aboveground biomass inversion.

[0009] (2) The present invention introduces a canopy effective thickness correction factor to suppress systematic deviations caused by factors such as canopy shading and insufficient echo penetration, which can reduce the risk of overestimation under complex shading conditions and enhance the robustness of results under different point cloud densities and different stand structures.

[0010] (3) This invention combines unsupervised contrastive learning to carry out domain adaptive optimization, which can improve the generalization ability of the model under cross-forest type, cross-region and cross-sensor conditions while reducing the dependence on labeled samples in the target region. It has good engineering application value and promotion prospects. Attached Figure Description

[0011] Figure 1 This is a flowchart illustrating the overall process of the multispectral three-dimensional spatial intelligent analysis algorithm for forest vegetation cover scenarios described in this invention.

[0012] Figure 2A schematic diagram of a 3D point cloud-spectral fusion voxel mesh generated by spatiotemporal alignment and physical quantity normalization of synchronously acquired LiDAR point cloud data and multispectral image data.

[0013] Figure 3 This diagram illustrates the process of extracting multispectral reflectance features and 3D geometric features from individual voxel cells in a voxel mesh, and performing cross-modal attention-weighted fusion.

[0014] Figure 4 This is a schematic diagram of a lightweight graph convolution-Transformer hybrid encoder, illustrating the information interaction between the graph convolution module and the Transformer module on a voxel grid.

[0015] Figure 5 A schematic diagram of the computational process for generating an aboveground biomass density field for a voxel regression head, including context-aware embedding vectors, canopy effective thickness correction factors, and their fusion results.

[0016] Figure 6 The image shows the forest aboveground biomass inversion result obtained by integrating the biomass density field along the vertical direction and aggregating it to the pixel scale, as well as the schematic flowchart of domain adaptive optimization using unsupervised contrastive learning. Detailed Implementation

[0017] This embodiment uses a forest sample area as the object, and the overall process is as follows: Figure 1 As shown. By simultaneously acquiring airborne or ground-based LiDAR point cloud data and multispectral image data, and following steps 1 to 5, the inversion of ground biomass from three-dimensional voxel features to pixel-scale is completed.

[0018] To facilitate reproduction and verification, Table 1 provides the data scale and key parameters used in this embodiment (which can be replaced with the same type of measured data and parameter range in actual applications).

[0019] In step 1, the multispectral image is first radiometrically calibrated and atmospherically corrected, and the LiDAR point cloud is then subjected to coordinate unification and outlier removal. Next, the two types of data are precisely registered to achieve sub-pixel alignment within a unified projection coordinate system. Finally, the space is voxelized according to a given voxel resolution to obtain a 3D point cloud-spectral fusion voxel mesh. ( Figure 2 ).

[0020] Subsequently, the key processing steps in step 1 are further explained. After radiometric calibration and atmospheric correction, the multispectral imagery is converted to surface reflectance. The LiDAR point cloud undergoes coordinate unification, noise removal, and conversion to three-dimensional coordinates according to a unified elevation datum. During registration, coarse registration can be completed first using sensor exterior orientation elements or control points, followed by mutual information maximization or iterative nearest point methods. Optimize and complete fine registration to achieve sub-pixel-level alignment of the two types of data in a unified projection coordinate system. Then, adjust according to voxel resolution. The space is divided into voxel units, and the set of points falling into the same voxel is denoted as . Then, the imaging model is used to project the points onto the image coordinates to obtain the reflectivity of the corresponding band, and the reflectivity is then... The reflectance in the data is weighted by the number of points or the echo intensity to form the multispectral channel characteristics of the voxel.

[0021] First, multispectral radiometric calibration and atmospheric correction are performed to convert the digital quantization values. Converted to radiance For example, using Further calculation of surface reflectance Its formula is: in This is the Earth-Sun distance correction term. The solar irradiance constant for the specified band. This is the solar zenith angle.

[0022] For the point cloud portion, coordinate projection transformation and elevation benchmark unification can be performed first, and outliers can be removed; statistical filtering or radius filtering can be used to remove noise points; after classifying ground points, a digital elevation model can be obtained, and the point cloud height can be normalized to the relative ground height. .

[0023] In the initial registration stage, if the sensor's exterior orientation elements are known, the point cloud coordinates can be mapped to the image coordinate system through rigid body transformation, for example... Initial alignment can also be achieved using control points or features with the same name.

[0024] Based on coarse registration, methods such as mutual information maximization, feature matching, or iterative nearest point can be used. Methods such as [list of methods] can be used to optimize registration parameters; when using height / intensity gratings derived from imagery and point clouds as registration objects, a pyramid strategy can be adopted to improve robustness from coarse to fine.

[0025] Given spatial range With resolution Calculate the voxel index for each point: And denote the set of points falling into the same voxel as Project each point onto an image pixel to obtain reflectivity. and to The internal reflectivity and intensity are weighted statistically to form the multi-channel properties of the fused voxel.

[0026] Table 1. Data Scale and Key Parameters of the Example After proceeding to step 2, the fused voxel mesh is... For each voxel unit, multispectral reflectance features are extracted. With three-dimensional geometric features . The former consists of reflectance in five bands: blue, green, red, near-infrared, and short-wave infrared. The latter is composed of statistics such as local normal vector entropy, rate of curvature change, vertical layer density, and mean and standard deviation of echo intensity. Subsequently, a cross-modal attention mechanism is used to... and Interactive weighted fusion is used to enable the model to adaptively emphasize more discriminative modal information under different forest types and shading conditions, with an output dimension of High-dimensional heterogeneous feature tensors ( Figure 3 ).

[0027] Step 2, feature construction, can utilize the original reflectance vector or, without increasing sensor requirements, introduce commonly used spectral indices as a supplement, thereby enhancing sensitivity to chlorophyll, water, and structural information. For example, it can calculate... , etc., among which It is a stable term.

[0028] The three-dimensional geometric structure features can be obtained as follows: For the set of interior points of voxels Principal component analysis of its neighborhood points yields the covariance matrix. eigenvalues Then the normal vector can be taken as the corresponding The eigenvectors, whose curvature can be defined as The entropy of the normal vector can be estimated from the Shannon entropy of the distribution of the normal vector direction.

[0029] The relative height within a voxel or pixel Divide the data into segments based on several height intervals, and statistically analyze the percentage of points in each segment to form a vertical hierarchical density vector; the number of segments can be chosen as appropriate. to To balance expressive power and data sparsity.

[0030] Point cloud echo intensity Calculate the mean with standard deviation To characterize differences in scattering features, indicators such as echo count and first-to-last echo ratio can be used to enhance the characterization of canopy structure.

[0031] Cross-modal attention fusion can take the following general form: first, ... and Linearly mapped to the same dimension respectively ,get and Then calculate the attention weights. The fusion representation is obtained by weighted summation. .in can be and Combining and generating allows for adaptive interaction between "using spectra to guide structures" or "using structures to guide spectra".

[0032] The normal vector can be obtained by using local principal component analysis to calculate the three-dimensional geometric structure features. With curvature Vertical stratification density can be based on height ranges. The proportion of statistical points; the mean value can be taken as the statistical measure of echo intensity. with standard deviation Multispectral characteristics With geometric features After linear mapping, the signal enters the cross-modal attention module, where attention weights are applied via... Normalization enables adaptive weighted fusion.

[0033] In step 3, the high-dimensional heterogeneous feature tensor is... Input lightweight graph convolution-Transformer hybrid encoder ( Figure 4 First, an adjacency graph is constructed with voxel centers as nodes, and adjacency edges are determined by a spatial distance threshold. The graph convolution module uses the GraphSAGE aggregation strategy to aggregate features in the local neighborhood to represent the local continuity and boundaries of the canopy structure. Then, the multi-head self-attention mechanism of the Transformer encoding layer is used to model long-range dependencies across voxels to characterize the consistency and variation of the stand structure on a larger spatial scale. The encoder outputs a context-aware embedding vector for each voxel. Its dimensions are .

[0034] The following explains the adjacency graph construction method and encoder implementation details in step 3: Can be adopted Neighborhood or distance threshold Define the adjacency set of voxels in a certain way ,when When establishing edges, among which The coordinates are the center coordinates of the voxels; if the voxels are sparse, they can also be based on... Neighbors construct graphs in embedded space or geometric space.

[0035] For example, layer The update can be written as in Either mean aggregation or max pooling aggregation can be used. It is a nonlinear function.

[0036] Multi-head self-attention is introduced based on the graph convolution output, while single-head attention is... To preserve the three-dimensional spatial position, the voxel center coordinates can be... The input features are then added or concatenated after being encoded by sine and cosine positions.

[0037] To reduce computational overhead, the maximum number of adjacencies for each voxel can be limited and sparse attention can be used. When the number of voxels is large, they can also be grouped and encoded by block or by height layer, and then global information interaction can be achieved through cross-block attention.

[0038] The construction of an adjacency graph can be achieved by using... Neighborhood or distance threshold :when At the voxel node and Establish edges between them; where The coordinates are the voxel center coordinates. GraphSAGE aggregation can use mean aggregation or max pooling aggregation to improve robustness; the position encoding of the Transformer can be obtained from... The sine and cosine codes are used to preserve three-dimensional spatial location information.

[0039] In step 4, based on the context-aware embedding vector Generate a continuously distributed aboveground biomass density field. Specifically, this involves... Input a voxel regression head (composed of fully connected layers and nonlinear activations) to obtain the voxel center coordinates. Predict the base density at the location; simultaneously calculate the horizontal position. Effective thickness of the canopy And a correction factor is introduced to suppress systematic bias caused by canopy shading. Figure 5 The calculation formula is as follows: Effective canopy thickness in step 4 It can be determined by the height set within the pixel. Statistics show that: ,in This represents the quantile operator. For regions with low point density or voids, neighborhood interpolation or morphological imputation can be used to obtain robustness. .

[0040] A voxel regression head can consist of two or more fully connected networks, with the last layer using... To ensure non-negativity. For further stabilization of training, [the following can be done]: Introducing upper bound constraints or log-space regression, for example, for prediction. Then the index is restored.

[0041] Exponential form of the correction factor This reflects the priori understanding that stronger occlusion results in greater attenuation; in another optional implementation, the correction factor can also be written as... The form is used to learn more flexible monotonic mappings, in which for function, It can be learned or set as a parameter.

[0042] in, voxel center coordinates Aboveground biomass density at the location and These are the learnable parameters of the voxel regression head. Take the LiDAR echo height within this pixel. Quantiles and The difference in quantiles, The attenuation coefficient is... It is a stable term; This is used to ensure that the density prediction is non-negative. Through end-to-end training, the model can automatically learn the coupling relationship between spectral-structural information and the correction factor under different canopy shading intensities. A smaller value usually indicates that the low-quantile echo rise or insufficient penetration leads to stronger canopy shading. The correction factor attenuates the density prediction for this area more strongly to reduce the systematic overestimation bias caused by shading.

[0043] Subsequently, in step 5, the biomass density field is discretely integrated along the vertical direction and aggregated to the pixel scale to obtain the forest aboveground biomass inversion results. ( Figure 6 ).

[0044] During aggregation, the voxel-to-pixel mapping method can be selected based on the task scale: if voxels and pixels are aligned at horizontal resolution, they can be directly mapped to the same... The voxel layers at each location are summed; if the voxels are finer, the set of voxels covered by the pixel is summed first and then normalized; if the voxels are coarser, the voxel results can be distributed to multiple pixels using area weighting.

[0045] For unsupervised contrastive learning, data augmentation can be constructed by combining the characteristics of remote sensing data: for multispectral channels, random noise perturbation, channel dropping, illumination scaling, or slight geometric perturbation can be applied; for point clouds, random downsampling, random rotation, random jittering, or intensity perturbation can be applied; and for voxel features, random occlusion can be applied to simulate missing data. Positive sample pairs are constructed using these augmentations, and negative sample sets are constructed by combining them with embeddings from different locations or at different times. This allows for the alignment of feature distributions without relying on target domain annotations.

[0046] In terms of training strategy, a two-stage process of "supervised pre-training followed by unsupervised adaptive fine-tuning" can be adopted: the first stage involves training in the labeled source domain... The primary optimization stage; the second stage involves joint optimization on mixed batches of the source and target domains. and This allows for stronger cross-domain generalization capabilities.

[0047] To facilitate understanding of the inversion results, Table 2 provides example sample data for scatter regression analysis (which can be generated from measured plot data or public datasets); Table 3 provides example ablation comparison results under different module combinations (the experimental setup can be reproduced on target / source domain data). It is evident that cross-modal attention, canopy correction factors, and domain adaptation strategies all play a positive role in improving accuracy.

[0048] During model training, supervised regression loss is used for labeled samples, let the first... The reference aboveground biomass per pixel is Predicted as For unlabeled samples, a contrastive loss is used, and the embeddings obtained from the two enhancements of the same pixel are assumed to be... and And construct a negative sample set for each anchor sample. Cosine similarity can be used. Temperature coefficient With weighting coefficients The loss function example is as follows: Table 2. Sample data for scatter regression analysis (partial). Table 3 Comparison of Accuracy in Example Ablation Experiments (Cross-Domain Testing) In summary, this invention can achieve deep fusion of multispectral information and point cloud information in three-dimensional voxel space, and realize the fine expression of complex canopy structures through joint modeling of local topological relationships and global dependencies. At the same time, by using the effective canopy thickness correction and unsupervised domain adaptive strategy, the stability and transferability of the method are improved under conditions of strong shading, few samples and cross-regional conditions.

Claims

1. A multispectral spatial analysis algorithm for forest vegetation cover, characterized in that, Includes the following steps: Step 1: Simultaneously acquire LiDAR point cloud data and multispectral image data of the forest area, perform radiometric calibration on the multispectral image data, and perform geometric registration, coordinate unification, and physical quantity normalization on the two types of data. Construct a 3D point cloud-spectral fusion voxel mesh under a unified geographic coordinate system. ; Step 2: In the fused voxel mesh In the process, multispectral reflectance features are extracted for each voxel unit. With three-dimensional geometric features Furthermore, a high-dimensional heterogeneous feature tensor is obtained by weighted fusion through a cross-modal attention mechanism. ; Step 3: Convert the high-dimensional heterogeneous feature tensor The input is a graph convolution-Transformer hybrid encoder. The graph convolution module captures the topological connectivity between neighboring voxels, and the Transformer module models long-range spatial dependencies. The output is a context-aware embedding vector. ; Step 4: Embedding Vectors Based on Context Awareness Voxel center coordinates predicted by voxel regression head The aboveground biomass density at a given location is calculated using the formula: [Formula omitted for brevity]. A continuous aboveground biomass density field is generated by introducing a canopy effective thickness correction factor. in, voxel center coordinates Aboveground biomass density at the location Embed vectors for context awareness. and For learnable parameters, This represents the effective thickness of the canopy at the corresponding horizontal position. The attenuation coefficient is... It is a stable term; .in, A smaller value usually indicates that the low-quantile echo rise or insufficient penetration leads to stronger canopy shading. The correction factor attenuates the density prediction for this area more strongly to reduce the systematic overestimation bias caused by shading. Step 5: Integrate the aboveground biomass density field along the vertical direction and aggregate it to the pixel scale to obtain the forest aboveground biomass inversion result. Furthermore, unsupervised contrastive learning is employed to perform domain-adaptive optimization in unlabeled regions, thereby enhancing generalization capabilities across forest types and sensor conditions.

2. The multispectral spatial analysis algorithm for forest vegetation cover according to claim 1, characterized in that, The multispectral image data includes blue, green, red, near-infrared, and short-wave infrared bands; the voxel resolution of the 3D point cloud-spectral fusion voxel grid meets the following requirements. and Each voxel unit contains a multispectral reflectance channel and a LiDAR point cloud distribution.

3. The algorithm according to claim 1, characterized in that, The multispectral reflectance characteristics It is a 5-dimensional vector; the three-dimensional geometric structure features It is a 7-dimensional vector containing local normal vector entropy, rate of curvature change, vertical layer density, and mean and standard deviation of echo intensity; the high-dimensional heterogeneous feature tensor output by the cross-modal attention mechanism. The dimension is .

4. The multispectral spatial analysis algorithm for forest vegetation cover according to claim 1, characterized in that, The graph convolution-Transformer hybrid encoder comprises two graph convolutional modules and two Transformer encoding layers; the graph convolutional modules employ the GraphSAGE aggregation strategy to obtain neighborhood features. The Transformer encoding layer includes a multi-head self-attention mechanism and a feedforward neural network; the context-aware embedding vector The dimension is .

5. The multispectral spatial analysis algorithm for forest vegetation cover according to claim 1, characterized in that, The voxel regression head consists of at least two fully connected layers, and outputs via... The mapping is non-negative; the effective thickness of the canopy Defined as the LiDAR echo height within this pixel Quantiles and Difference of quantiles; the attenuation coefficient The stabilizing term .

6. The multispectral spatial analysis algorithm for forest vegetation cover according to claim 1, characterized in that, aboveground biomass at the pixel scale The formula is obtained by discretizing and integrating the aboveground biomass density field along the vertical direction: in, For the first The central height of the individual's pithelial layer, For voxel height, This represents the number of voxel layers from the ground to the top of the canopy.

7. A multispectral spatial analysis algorithm for forest vegetation cover according to claim 1, characterized in that, The unsupervised contrastive learning uses augmented data as positive sample pairs and embedding vectors from different spatial locations or at different times as negative samples. It optimizes the parameters of the graph convolution-Transformer hybrid encoder by maximizing the similarity of positive samples and minimizing the similarity of negative samples.

8. The multispectral spatial analysis algorithm for forest vegetation cover according to claim 1, characterized in that, The multispectral image data and LiDAR point cloud data use a unified projection coordinate system and elevation datum, and the voxel mesh is fused through geometric registration. The multispectral reflectance channels in the image are strictly aligned with the LiDAR point cloud distribution in the geographic coordinate system.

Images