A forest carbon sink change trend prediction method based on multi-source heterogeneous data analysis

Abstract

This invention relates to the field of ecological and environmental protection technology, and in particular to a method for predicting forest carbon sink change trends based on multi-source heterogeneous data analysis. The method includes: acquiring multi-source remote sensing images and completing geometric correction and registration to generate a spatially aligned raster set; extracting canopy boundaries to form a forest mask; extracting three-dimensional structural parameters of the tree canopy using LiDAR point cloud data; identifying tree species and ages by combining color histograms, texture features, and multispectral features; calculating static carbon storage and generating a carbon density map based on canopy area and carbon density coefficient; identifying transpiration anomalies and calculating leaf water potential by combining thermal infrared images to predict carbon absorption rates; and superimposing the carbon increment onto the carbon density map to form a time series, identifying carbon sink hotspots, transition zones, and carbon-depleted zones through difference and clustering. This invention, by integrating multi-source remote sensing data and efficient algorithms, achieves real-time monitoring, dynamic prediction, and refined management of forest carbon sinks, providing accurate data support for carbon markets, ecological compensation, and carbon neutrality decision-making.

Description

Technical Field

[0001] This invention relates to the field of ecological and environmental protection, and in particular to a method for predicting the trend of forest carbon sink changes based on multi-source heterogeneous data analysis. Background Technology

[0002] Traditional forest carbon sequestration monitoring relies on ground-based sample plot surveys and fixed flux observation stations. While these methods offer high accuracy, their spatial coverage is limited, making it difficult to achieve large-scale, real-time dynamic monitoring. Furthermore, carbon estimation methods that depend on a single remote sensing data source are affected by factors such as meteorological conditions, topographical obstruction, and mixed tree species, leading to significant fluctuations in monitoring accuracy.

[0003] Furthermore, multi-source remote sensing data often suffers from spatial registration difficulties, asynchronous acquisition times, and imperfect feature fusion, making carbon sink estimation results prone to bias and instability in practical applications, and failing to meet the needs of carbon markets, ecological compensation, and carbon neutrality pathway decisions. Most existing carbon monitoring methods remain at the regional average estimation level, lacking comprehensive analysis and dynamic updates of individual tree-scale structural characteristics, tree species attributes, and physiological states. Especially in thermal infrared remote sensing and transpiration monitoring, traditional methods cannot effectively identify tree transpiration anomalies and changes in carbon absorption status, resulting in insufficient timeliness of dynamic carbon sink assessments. Existing methods also lack the ability to perform long-term serialization and trend identification of carbon storage and carbon absorption rates, failing to achieve spatial cluster analysis of carbon sink dynamic evolution, making it difficult to locate carbon sink hotspots, transition zones, and carbon-depleted areas, and lacking technical support for early warning and regulation of carbon sink changes. Therefore, although remote sensing technology has provided new ideas for monitoring forest carbon sinks, existing methods are still unable to effectively solve the problem of monitoring large-scale, dynamically changing forest carbon sinks. There is an urgent need for an innovative technical solution that can combine multi-source remote sensing data to achieve accurate dynamic monitoring and prediction of forest carbon sinks. Summary of the Invention

[0004] This invention addresses the problems existing in the prior art by providing a method for predicting forest carbon sink change trends based on multi-source heterogeneous data analysis, mainly comprising: Raw multi-source remote sensing images are acquired using drones, and spatial alignment is generated by synchronizing with sensor calibration data using timestamps and performing geometric correction and pixel-level registration. Based on the NDVI index layers in the registered RGB and multispectral images, generate a binary forest canopy mask image of the canopy boundary; Based on the registered LiDAR point cloud data, extract the point cloud subset of the corresponding region of the canopy mask, generate the main direction projection image, obtain the segmentation contour and center point coordinate label of the individual tree, and construct the three-dimensional structural parameter set of the tree canopy. Based on RGB color histograms, LBP and GLCM texture features, multispectral features, and combined with crown width and height ratio data, tree species and age are identified. Based on the three-dimensional structural parameters of trees, tree species and age, the carbon density per unit crown width of trees is predicted, and the static carbon storage of each tree is calculated. Based on the acquisition of canopy surface temperature and identification of transpiration anomalies by thermal infrared images, combined with leaf water potential status indicators, tree species and age, three-dimensional structural parameters of the canopy and daily carbon absorption data, the carbon absorption rate of the canopy is predicted. Based on the static carbon storage and carbon absorption rate of each tree, a dynamic carbon storage time series layer is constructed. Combined with carbon storage difference monitoring and carbon sink status zoning layer, carbon sink hotspots and carbon deficit areas are identified.

[0005] Furthermore, the step of acquiring raw multi-source remote sensing images using a drone and generating spatial alignment by synchronizing them with sensor calibration data using timestamps and performing geometric correction and pixel-level registration includes: Based on the coordinate range information of the forest area to be monitored, a drone equipped with a high-resolution RGB camera, a multispectral camera, a hyperspectral imager, and a thermal infrared sensor was used to acquire original multi-source remote sensing image data covering the entire target forest area, including RGB images, multispectral images, and thermal infrared images. Time-stamp synchronization and sensor calibration data were used to perform time alignment and geometric correction on images from different sensors to obtain a spatially consistent remote sensing image dataset. An image registration algorithm was used to perform pixel-level joint registration of RGB, multispectral, and thermal infrared images, and nearest neighbor interpolation was used for resampling to obtain a spatially aligned multi-source sensing image raster set.

[0006] Further, the step of generating a binary forest canopy mask image of the canopy boundary based on the NDVI index layer in the registered RGB image and multispectral image includes: Based on the NDVI index layers in the registered RGB and multispectral images, the Sobel operator is used to perform edge gradient enhancement processing on the RGB images to extract high-frequency edge regions. Through a region growing algorithm, pixels with NDVI greater than a preset NDVI threshold are used as seed regions, and a growth threshold is set to expand adjacent pixels into 8-neighborhoods. By performing logical operations on the binary image of the Sobel edge detection result and the mask image generated by NDVI region growing at the pixel level, the output is set to 1 when the values ​​of the two layers are both 1 at the same pixel position, and 0 otherwise, a binary forest canopy mask image containing only the complete canopy boundary is obtained, and the minimum bounding rectangle boundary information of each connected component is preserved.

[0007] Further, based on the registered LiDAR point cloud data, a subset of point clouds corresponding to the canopy mask region is extracted to generate a main direction projection image, obtain the segmentation contour and center point coordinate labels of individual trees, and construct a three-dimensional structural parameter set for the tree canopy, including: Based on the registered LiDAR point cloud data, a subset of point clouds corresponding to the canopy mask region is extracted. Principal component analysis (PCA) is used to estimate the principal orientation of the point clouds within each connected canopy region, obtaining the principal axis direction vector of the point clouds in that region. A two-dimensional local orthographic projection coordinate system is constructed using the principal axis direction, and the point clouds are projected onto the principal orientation plane along the principal axis direction. Based on a preset rasterization resolution, the principal orientation projection image of the region is obtained. The Euclidean distance field image of each connected region is calculated through distance transformation, and local maxima are extracted as initial seed points for the watershed algorithm. The watershed segmentation algorithm is used to process the Euclidean distance field image, based on a preset damping threshold. The correspondence between the label layer and the original image coordinates is obtained, and the segmentation outline, polygon boundary coordinates, and cropped image blocks of each tree are obtained. The corresponding coordinate labels are recorded to obtain the spatial coordinates of the tree. The corresponding coordinate labels are the center point coordinates of the cropped image blocks. The crown width and crown area of ​​each tree are calculated using a boundary tracking algorithm. The crown height is extracted by combining the difference between the original digital surface model and the digital terrain model. A three-dimensional structural parameter set for each crown is constructed and saved to the forest area monitoring database. The three-dimensional structural parameters include crown area, crown height, crown width, crown index, and height ratio. The height ratio is the ratio of crown width to tree height.

[0008] Furthermore, the process of identifying tree species and age based on RGB color histograms, LBP and GLCM texture features, multispectral features, and combined with crown width and height ratio data includes: Color histograms are extracted from the RGB images of each tree to obtain color distribution feature vectors. Texture features are calculated based on LBP and GLCM, and energy, contrast, and entropy are extracted. Multispectral features, including NDVI, NDRE, and EVI, are extracted from multispectral images and saved to the forest area monitoring database. Historical data on color distribution feature vectors, texture features, multispectral features, crown width, and height ratio are obtained from the forest area monitoring database. Tree species and ages are labeled, and a tree species identification model is constructed using the random forest algorithm. Based on the color distribution feature vectors, texture features, multispectral features, crown width, and height ratio of the real-time acquired tree data, the tree species and age are identified using the tree species identification model.

[0009] Furthermore, the prediction of carbon density per unit canopy width of trees based on their three-dimensional structural parameters, tree species, and age, and the calculation of the static carbon storage of each tree, include: By using the forest area monitoring database, the three-dimensional structural parameters, tree species, and tree age of trees are obtained, and the corresponding unit crown carbon density of trees is labeled. A recurrent neural network is used to train the model and construct a unit crown carbon density prediction model. Based on the real-time acquired three-dimensional structural parameters, tree species, and tree age of trees, the unit crown carbon density prediction model is used to predict the current unit crown carbon density of trees. The static carbon storage of each tree is calculated by multiplying the crown area by the unit crown carbon density and combining the area data extracted from the segmentation boundary of a single tree. By traversing the spatial location and carbon storage of each tree, a carbon storage raster layer is constructed and reprojected onto the spatial coordinates of the trees to obtain a spatial distribution map of carbon storage in the forest area.

[0010] Furthermore, the step of obtaining canopy surface temperature and determining transpiration anomalies based on thermal infrared images, and then predicting the canopy carbon absorption rate by combining leaf water potential state indicators, tree species and age, three-dimensional structural parameters of the canopy, and daily carbon absorption data, includes: Based on the pixel region corresponding to the spatial location of each tree in the thermal infrared image, the canopy surface temperature value is obtained; by comparing the temperature value of each tree with the average temperature of the same type of trees within a preset distance, it is determined whether there is an abnormality in transpiration, and the location of the abnormal trees is marked; the average pixel temperature of the canopy of each tree with abnormal transpiration is calculated by the regional mean method, and the leaf surface water potential index of the canopy of each tree with abnormal transpiration is calculated; combined with the tree species, age, three-dimensional structural parameters of the canopy, and daily carbon absorption data measured by flux observation tower and ground, a recurrent neural network is used to train the model and construct a tree carbon absorption rate prediction model to predict the carbon absorption rate of each tree canopy; It also includes calculating the average pixel temperature of the canopy of each tree with abnormal transpiration using the regional mean method, and calculating the leaf surface water potential state index of the canopy of each tree with abnormal transpiration, specifically including: Based on the average pixel temperature of the canopy of each tree with abnormal transpiration, combined with the canopy temperature and transpiration model... Calculate latent heat flux ,in, Net radiation from the leaf surface is calculated using measured meteorological parameters or obtained from a drone's light sensor. The sensible heat flux was calculated using the temperature difference between the canopy and the air, as well as aerodynamic drag parameters. Soil heat flux; net canopy radiation flux was calculated using light intensity measurements and reflectance, and then expressed as a leaf water potential index. Calculate the leaf surface water potential state index of the canopy of each tree with abnormal transpiration. ,in, Leaf conductance was determined using tree species characteristics and temperature and humidity response curves. Δ is the pressure constant, determined by air pressure, and Δ is the slope of the saturated water vapor pressure curve, calculated based on air temperature. The average temperature of the canopy. The ambient air temperature is obtained from ground weather stations or drone sensors.

[0011] Furthermore, based on the static carbon storage and carbon absorption rate of each tree, a dynamic carbon storage time series layer is constructed, and combined with carbon storage difference monitoring and carbon sink status zoning layer, carbon sink hotspots and carbon deficit areas are identified, including: Based on the static carbon storage and carbon absorption rate of each tree, the carbon increments in preset time units are sequentially superimposed onto the static carbon density map to construct a dynamic carbon storage time series layer group for preset time units. A pixel-based back-difference calculation method is used to monitor the carbon storage difference at different time periods. By setting change threshold intervals, the carbon sink growth, carbon sink stability, or carbon sink decline status within a region is determined, resulting in a carbon sink status zoning layer. Based on the carbon sink status zoning layer, a density peak clustering algorithm is applied to spatially cluster the carbon sink growth, carbon sink stability, and carbon sink decline regions respectively. A sample feature space is constructed using the Euclidean distance matrix and carbon change density values. Local density and minimum distance indices are calculated, and cluster centers are determined using the truncation distance method. Clustering operations are then performed to obtain distribution maps of carbon sink hotspots, transition zones, and carbon deficit zones. The total area and monthly average carbon increment of each region are also statistically analyzed.

[0012] The technical solutions provided by the embodiments of the present invention may include the following beneficial effects: This invention provides a method for dynamic monitoring and prediction of forest carbon sinks based on multi-source remote sensing data. By generating a binary forest canopy mask image of the canopy boundary using an NDVI index layer, and extracting the three-dimensional structural parameters of individual trees using LiDAR point cloud data, this method can accurately identify the spatial location and canopy morphology of each tree. Furthermore, by combining RGB color histograms, LBP and GLCM texture features, and multispectral features, this method can efficiently identify tree species and ages, predict the carbon density per unit canopy width based on these features, and calculate the static carbon storage of each tree. By calculating the canopy temperature and latent heat flux of trees with abnormal transpiration, and combining thermal infrared imagery and canopy net radiation flux, this method monitors the leaf surface water potential status indicators of trees in real time, enabling timely detection of transpiration anomalies and prediction of tree carbon absorption rates, significantly improving the timeliness and accuracy of carbon sink assessment. This invention constructs a dynamic carbon storage time series layer based on the static carbon storage and carbon absorption rate of each tree. Combined with carbon storage difference monitoring and a carbon sink status zoning layer, it accurately identifies carbon sink hotspots and carbon-deficient areas. This method enables refined dynamic management of forest carbon sinks, optimizing carbon sink protection and carbon emission reduction strategies. The technical solution of this invention, by integrating multi-source remote sensing data and efficient algorithms, overcomes the limitations of traditional carbon sink monitoring, achieving real-time monitoring, dynamic prediction, and refined management of forest carbon sinks, providing precise data support for carbon markets, ecological compensation, and carbon neutrality decision-making. Attached Figure Description

[0013] Figure 1 This is a flowchart of a method for predicting forest carbon sink change trends based on multi-source heterogeneous data analysis according to the present invention. Detailed Implementation

[0014] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be described in detail below with reference to the accompanying drawings and specific embodiments.

[0015] like Figure 1 This embodiment of a method for predicting forest carbon sink change trends based on multi-source heterogeneous data analysis may specifically include: Step S101: Use a drone to acquire raw multi-source remote sensing images, and use timestamps to synchronize with sensor calibration data to perform geometric correction and pixel-level registration to generate spatial alignment.

[0016] Based on the coordinate range information of the forest area to be monitored, a drone equipped with a high-resolution RGB camera, a multispectral camera, a hyperspectral imager, and a thermal infrared sensor was used to acquire original multi-source remote sensing image data covering the entire target forest area, including RGB images, multispectral images, and thermal infrared images. Time-stamp synchronization and sensor calibration data were used to perform time alignment and geometric correction on images from different sensors to obtain a spatially consistent remote sensing image dataset. An image registration algorithm was used to perform pixel-level joint registration of RGB, multispectral, and thermal infrared images, and nearest neighbor interpolation was used for resampling to obtain a spatially aligned multi-source sensing image raster set.

[0017] For example, taking the mangrove forest area of ​​Futian Mangrove Ecological Park in Shenzhen as an example, carbon sink monitoring of this forest area was carried out using a UAV equipped with a high-resolution RGB camera (0.03 m resolution), a five-band multispectral camera (0.20 m resolution), a hyperspectral imager (400–1000 nm band resolution, 0.50 m resolution), and a thermal infrared sensor (8–14 μm band, temperature sensitivity 0.05 ℃). Full-coverage aerial survey was performed at a flight altitude of 120 m, with a flight path overlap rate of 80% and a lateral overlap rate of 70%. A total of 896 RGB images, 320 multispectral images, 205 hyperspectral stripes, and 240 thermal infrared images were collected, with a total raw image data volume of approximately 46 GB. The UAV is equipped with a unified timing module, synchronizing the acquisition time of different sensors to ±2 ms. Geometric correction of the acquired data is performed using 12 ground-based control points and pre-calibrated intrinsic and extrinsic parameters, ensuring that the spatial deviation of all sensor data in the CGCS2000 / UTM 49N coordinate system does not exceed 0.5 pixels. Subsequently, using RGB imagery as a reference, pixel-level registration of multispectral and thermal infrared imagery is achieved based on feature point matching, with the error controlled within 1 pixel. Nearest neighbor interpolation is then used to resample the multispectral and thermal infrared data to 0.03 m, ultimately yielding a coverage area of ​​approximately 6.5 km². 2 Multi-source remote sensing image raster set.

[0018] Step S102: Generate a binary forest canopy mask image of the canopy boundary based on the NDVI index layer in the registered RGB image and multispectral image.

[0019] Based on the NDVI index layers in the registered RGB and multispectral images, the Sobel operator is used to perform edge gradient enhancement processing on the RGB images to extract high-frequency edge regions. Through a region growing algorithm, pixels with NDVI greater than a preset NDVI threshold are used as seed regions, and a growth threshold is set to expand adjacent pixels into 8-neighborhoods. By performing logical operations on the binary image of the Sobel edge detection result and the mask image generated by NDVI region growing at the pixel level, the output is set to 1 when the values ​​of the two layers are both 1 at the same pixel position, and 0 otherwise, a binary forest canopy mask image containing only the complete canopy boundary is obtained, and the minimum bounding rectangle boundary information of each connected component is preserved.

[0020] For example, using registered RGB images with a resolution of 0.03m and multispectral images with a resolution of 0.20m, an NDVI index layer is calculated, and the NDVI threshold is set to 0.25 to distinguish the mangrove canopy area from the bare land and water background. The Sobel operator is applied to the RGB image for gradient enhancement, with the horizontal and vertical gradient thresholds set to 20 and 25 respectively when extracting edges to highlight high-frequency information at the canopy boundaries. Pixels with an NDVI greater than 0.25 are used as seed points, and a growth threshold of 0.05 is set. Region expansion is performed in an 8-neighborhood to obtain a continuous initial canopy mask area. The binary edge map obtained from Sobel edge detection is logically ANDed with the region growth result; only pixels with the same pixel position in both layers are retained, resulting in a binary forest mask map with closed canopy edges. To address reflective interference, the Otsu adaptive thresholding method is used to segment high-reflectivity areas, with a grayscale threshold of 180 for detecting reflective pixels, extracting solar reflective areas from the mangrove canopy. Subsequently, a morphological dilation operation was performed on the reflective areas. The structuring element was 3×3, and the dilation was 1 pixel, which moderately expanded the reflective area mask to ensure that the reflective edge areas were also occluded. This dilated mask was then subtracted from the forest binary mask to obtain a canopy mask that eliminated reflective interference and understory gaps. This mask contained approximately 18,500 connected regions. The minimum bounding rectangle was extracted for each region to obtain the boundary coordinate set.

[0021] Step S103: Based on the registered LiDAR point cloud data, extract the point cloud subset of the corresponding region of the canopy mask, generate the main direction projection image, obtain the segmentation contour and center point coordinate label of the single tree, and construct the three-dimensional structural parameter set of the tree canopy.

[0022] Based on the registered LiDAR point cloud data, a subset of point clouds corresponding to the canopy mask region is extracted. Principal component analysis (PCA) is used to estimate the principal orientation of the point clouds within each connected canopy region, obtaining the principal axis direction vector of the point clouds in that region. A two-dimensional local orthographic projection coordinate system is constructed using the principal axis direction, and the point clouds are projected onto the principal orientation plane along the principal axis direction. Based on a preset rasterization resolution, the principal orientation projection image of the region is obtained. The Euclidean distance field image of each connected region is calculated through distance transformation, and local maxima are extracted as initial seed points for the watershed algorithm. The watershed segmentation algorithm is used to process the Euclidean distance field image, based on a preset damping threshold. The correspondence between the label layer and the original image coordinates is obtained, and the segmentation outline, polygon boundary coordinates, and cropped image blocks of each tree are obtained. The corresponding coordinate labels are recorded to obtain the spatial coordinates of the tree. The corresponding coordinate labels are the center point coordinates of the cropped image blocks. The crown width and crown area of ​​each tree are calculated using a boundary tracking algorithm. The crown height is extracted by combining the difference between the original digital surface model and the digital terrain model. A three-dimensional structural parameter set for each crown is constructed and saved to the forest area monitoring database. The three-dimensional structural parameters include crown area, crown height, crown width, crown index, and height ratio. The height ratio is the ratio of crown width to tree height.

[0023] For example, based on UAV LiDAR point cloud data of a monitoring area in a mangrove ecological park, a subset of point clouds for the corresponding area was selected according to the previously extracted canopy mask, containing approximately 128 million 3D points. For each connected canopy region, PCA principal component analysis was used to calculate the principal direction, where the principal direction vector of a certain tree is (0.82, 0.57, 0.06), which represents the main extension direction of the canopy in the horizontal plane. A local 2D orthographic projection coordinate system was established based on this principal direction, and the corresponding point cloud was projected onto the principal direction plane. The rasterization resolution was set to 0.05 m to obtain the principal direction projection image of the canopy region. By performing a distance transformation on the projection image, a Euclidean distance field image was obtained, where the pixel coordinates corresponding to the local maximum point are (145, 212), which was used as the seed point for watershed segmentation. A damping threshold of 0.15 was used for watershed segmentation to accurately separate adjacent canopies and obtain the boundary contour of each individual tree. For example, for the single tree numbered T063, the polygon formed by its boundary points has a perimeter of 8.4 m and a corresponding area of ​​4.7 m². 2 The maximum crown width is 2.6 m. Calculations using the difference between DSM and DTM show a tree height of 3.9 m, a height ratio (crown width / tree height) of 0.67, and a crown index (area / height). 2The value is 0.31. For all segmented individual trees, the center point coordinates are extracted as corresponding labels. For example, the center point coordinates of T063 are (114.03652°, 22.52498°). Finally, the canopy area, canopy height, canopy width, canopy index, height ratio, and corresponding spatial coordinate information of all tree crowns are written into the forest area monitoring database, forming a three-dimensional structural parameter set for approximately 22,000 mangrove individuals.

[0024] Step S104: Based on the RGB color histogram, LBP and GLCM texture features, multispectral features, and combined with crown width and height ratio data, identify the tree species and age.

[0025] Color histograms are extracted from the RGB images of each tree to obtain color distribution feature vectors. Texture features are calculated based on LBP and GLCM, and energy, contrast, and entropy are extracted. Multispectral features, including NDVI, NDRE, and EVI, are extracted from multispectral images and saved to the forest area monitoring database. Historical data on color distribution feature vectors, texture features, multispectral features, crown width, and height ratio are obtained from the forest area monitoring database. Tree species and ages are labeled, and a tree species identification model is constructed using the random forest algorithm. Based on the color distribution feature vectors, texture features, multispectral features, crown width, and height ratio of the real-time acquired tree data, the tree species and age are identified using the tree species identification model.

[0026] For example, using RGB images of a mangrove monitoring area in a mangrove ecological park, a color histogram is extracted from the corresponding RGB image of a single tree, numbered T063. The histogram is divided into 256 brightness levels, with peak values ​​for the red, green, and blue channels appearing at 122, 136, and 115, respectively, corresponding to the typical green reflectance characteristics of the canopy leaves. These values ​​are combined into a 768-dimensional color distribution feature vector. The local texture structure of the image is calculated based on LBP local binary mode, and texture statistical features are calculated using the GLCM gray-level co-occurrence matrix. Texture indices of 0.413 energy, 6.27 contrast, and 1.98 are obtained at a distance of 2 pixels and an angle of 0°. These values ​​reflect the uniformity and directionality of the canopy leaf texture. Based on the corresponding multispectral image, the NDVI, NDRE, and EVI feature values ​​of T063 are calculated to be 0.76, 0.42, and 0.63, respectively, indicating high photosynthetic activity of its leaves. These feature data, along with their crown width of 2.6 m and height ratio of 0.67, were stored in the forest area monitoring database and compared with historical samples. The database contains feature samples of approximately 15,000 trees, covering major mangrove species such as Kandelia candel, Avicennia marina, and Excoecaria agallocha. Historical data on color distribution feature vectors, texture features, multispectral features, crown width, and height ratio were obtained from the forest area monitoring database. Tree species and ages were labeled, and a random forest algorithm was used to train the model, constructing a tree species identification model. During real-time identification, the model classified the input features of T063, ultimately outputting that the tree species was Kandelia candel, with a predicted age of 9 years, consistent with the ground survey record. This method enables automated identification and age estimation of individual tree species across an entire forest area.

[0027] Step S105: Based on the three-dimensional structural parameters of the tree, tree species and age, predict the carbon density per unit crown width of the tree, and calculate the static carbon storage of each tree.

[0028] By using the forest area monitoring database, the three-dimensional structural parameters, tree species, and tree age of trees are obtained, and the corresponding unit crown carbon density of trees is labeled. A recurrent neural network is used to train the model and construct a unit crown carbon density prediction model. Based on the real-time acquired three-dimensional structural parameters, tree species, and tree age of trees, the unit crown carbon density prediction model is used to predict the current unit crown carbon density of trees. The static carbon storage of each tree is calculated by multiplying the crown area by the unit crown carbon density and combining the area data extracted from the segmentation boundary of a single tree. By traversing the spatial location and carbon storage of each tree, a carbon storage raster layer is constructed and reprojected onto the spatial coordinates of the trees to obtain a spatial distribution map of carbon storage in the forest area.

[0029] For example, by retrieving the three-dimensional structural parameters of individual trees, tree species, and tree age from the forest area monitoring database, and labeling the corresponding estimated values ​​of carbon density per unit crown width, a recurrent neural network was used to train the model, constructing a carbon density prediction model per unit crown width. The mean absolute error was 0.9 kgC / m² on the validation set. 2 Subsequently, the parameters of individual trees entered into the database in real time were predicted. For example, tree number T063, species *Kandelia candel*, age 9 years, canopy area 4.7 m², was analyzed. 2 The predicted carbon density per unit crown width is 11.8 kgC / m³ for a crown width of 2.6 m and a tree height of 3.9 m. 2 Based on this, the static carbon storage was calculated to be 4.7 × 11.8 = 55.46 kgC. Using the same method, the static carbon storage and spatial coordinates of each of the 22,000 individual trees were obtained. The discrete individual tree carbon content was rasterized at a resolution of 0.5 m × 0.5 m in the CGCS2000 / UTM 49N coordinate system. Using area weighting, the carbon content of individual trees across pixels was allocated to corresponding pixels according to the canopy projection ratio, forming a carbon storage raster layer. A uniform reprojection and time labeling were performed on all layers, ultimately generating a coverage of approximately 6.5 km. 2 The spatial distribution map of carbon reserves in the region shows a total carbon sink of approximately 1,250 tC. Each pixel records the carbon density value, carbon sink amount, source tree ID, and calculation batch number at the corresponding time.

[0030] Step S106: Obtain the canopy surface temperature based on thermal infrared images and determine transpiration anomalies. Combine this with leaf water potential status indicators, tree species and age, three-dimensional structural parameters of the canopy, and daily carbon absorption data to predict the carbon absorption rate of the canopy.

[0031] Based on the pixel regions corresponding to the spatial location of each tree in the thermal infrared image, the canopy surface temperature value is obtained. By comparing the temperature value of each tree with the average temperature of the same species within a preset distance range, it is determined whether there is an abnormality in transpiration, and the location of the abnormal tree is marked. The average pixel temperature of the canopy of each tree with abnormal transpiration is calculated using the regional mean method, and the leaf surface water potential index of the canopy of each tree with abnormal transpiration is calculated. Combining the tree species, age, three-dimensional structural parameters of the canopy, and daily carbon absorption data measured by flux observation tower and ground, a recurrent neural network is used to train the model and construct a tree carbon absorption rate prediction model to predict the carbon absorption rate of each tree canopy.

[0032] For example, for a single Kandelia candel tree (T063), the pixel temperature of the corresponding canopy area was extracted using thermal infrared imagery. Each pixel was 0.03 m × 0.03 m, containing 245 pixels, with a temperature range between 29.6 ℃ and 31.2 ℃ and an average temperature of 30.3 ℃. Within a 10 m buffer zone, the average temperature of 38 reference trees of the same species was retrieved, with an average temperature of 28.9 ℃ and a temperature difference of 1.4 ℃, exceeding the set transpiration anomaly threshold of 1.0 ℃. Therefore, T063 was identified as a tree with abnormal transpiration and an anomaly label was added to the database. Subsequently, leaf water potential indicators were calculated based on the canopy pixel average of 30.3 ℃, combined with the measured air pressure of 101.3 kPa, ambient air temperature of 28.5 ℃, and stomatal conductance of 0.35 mol·m³ in the area. -2 ·s -2 With a saturated vapor pressure curve slope of 0.19 kPa / ℃, the leaf surface water potential state index was calculated to be -1.24 MPa. Based on the leaf surface water potential state index, combined with the three-dimensional structural parameters of tree species, age, and crown, as well as daily carbon uptake data measured by flux monitoring towers and the ground, a recurrent neural network was used to train the model and construct a tree carbon uptake rate prediction model. Real-time data obtained included tree species Kandelia candel, tree age of 9 years, crown width of 2.6 m, tree height of 3.9 m, and crown area of ​​4.7 m². 2 Structural parameters, and the daily carbon uptake flux data of 7.2 μmol CO2·m2 measured by the same flux observation tower. -2 ·s -1 These features were input into a pre-trained tree carbon uptake rate prediction model, and the model output that the carbon uptake rate corresponding to the T063 canopy was 6.83 μmol CO2·m⁻¹. -2 ·s -1 The corresponding carbon absorption per unit time is approximately 0.29 gC·m. -2 ·h -1 Using this method, carbon absorption rate distribution maps with anomaly markers were generated for each of the 22,000 trees in the area.

[0033] The average pixel temperature of the canopy of each tree with abnormal transpiration was calculated using the regional mean method, and the leaf surface water potential index of the canopy of each tree with abnormal transpiration was calculated.

[0034] Based on the average pixel temperature of the canopy of each tree with abnormal transpiration, combined with the canopy temperature and transpiration model... Calculate latent heat flux ,in, Net radiation from the leaf surface is calculated using measured meteorological parameters or obtained from a drone's light sensor. The sensible heat flux was calculated using the temperature difference between the canopy and the air, as well as aerodynamic drag parameters. Soil heat flux; net canopy radiation flux was calculated using light intensity measurements and reflectance, and then expressed as a leaf water potential index. Calculate the leaf surface water potential state index of the canopy of each tree with abnormal transpiration. ,in, Leaf conductance was determined using tree species characteristics and temperature and humidity response curves. Δ is the pressure constant, determined by air pressure, and Δ is the slope of the saturated water vapor pressure curve, calculated based on air temperature. The average temperature of the canopy. The ambient air temperature is obtained from ground weather stations or drone sensors.

[0035] For example, when analyzing a tree with abnormal transpiration, numbered T092, the average temperature of pixels within the canopy area was extracted using thermal infrared imagery to obtain the average canopy temperature. The temperature was 33.2℃. According to ground weather station records, the ambient air temperature during this period was... The temperature is 28.7℃, the actual air pressure is 98.2 kPa, and the corresponding pressure constant is... The saturated vapor pressure curve slope Δ is calculated to be 0.21 kPa / °C based on the air temperature, with a value of 0.066 kPa / °C. The net radiative flux of the tree canopy is calculated using solar irradiance recorded by a light sensor and visible light reflectance from a UAV. 410W / m 2 Sensible heat flux is estimated using temperature difference and measured wind speed. 126W / m 2 Soil heat flux The estimated value is 25 W / m 2 Substituting these data into the canopy temperature and transpiration model The latent heat flux was found to be 259 W / m. 2 Based on the characteristics of Kandelia candel tree species labeled in the database and matching their temperature and humidity response curves, leaf conductance was determined. 0.25 mol·m -2 ·s -1 According to the formula of leaf surface water potential state index The leaf surface water potential index of each tree with abnormal transpiration was calculated. The leaf surface water potential index of the abnormal tree was 101.08. This value deviates significantly from the upper limit of the normal range for this tree species. The normal range is usually 85~90. This indicates that its transpiration efficiency is abnormally high and the canopy water potential is relatively unbalanced. It can be used as an important reference indicator for abnormal carbon absorption rate.

[0036] Step S107: Based on the static carbon storage and carbon absorption rate of each tree, construct a dynamic carbon storage time series layer, and combine it with carbon storage difference monitoring and carbon sink status zoning layer to identify carbon sink hotspots and carbon deficit areas.

[0037] Based on the static carbon storage and carbon absorption rate of each tree, the carbon increments in preset time units are sequentially superimposed onto the static carbon density map to construct a dynamic carbon storage time series layer group for preset time units. A pixel-based back-difference calculation method is used to monitor the carbon storage difference at different time periods. By setting change threshold intervals, the carbon sink growth, carbon sink stability, or carbon sink decline status within a region is determined, resulting in a carbon sink status zoning layer. Based on the carbon sink status zoning layer, a density peak clustering algorithm is applied to spatially cluster the carbon sink growth, carbon sink stability, and carbon sink decline regions respectively. A sample feature space is constructed using the Euclidean distance matrix and carbon change density values. Local density and minimum distance indices are calculated, and cluster centers are determined using the truncation distance method. Clustering operations are then performed to obtain distribution maps of carbon sink hotspots, transition zones, and carbon deficit zones. The total area and monthly average carbon increment of each region are also statistically analyzed.

[0038] For example, based on the static carbon storage and carbon uptake rate of 22,000 mangrove individuals in the monitoring database, the monthly carbon increment was superimposed on the initial static carbon density map to construct a six-period carbon storage time series raster layer. For instance, the initial static carbon storage of a single Kandelia candel (T063) was 55.5 kgC, with a monthly carbon uptake increment of 4.3 kgC, increasing its cumulative carbon storage to 81.3 kgC by the end of June. Based on the pixel-based back-difference calculation method, the carbon storage difference between the monthly raster layers was calculated. A growth threshold range of [-0.5, 0.5] kgC / m² was set as the stable range; values ​​exceeding 0.5 were defined as growth zones, and values ​​below -0.5 were defined as decay zones. This resulted in a carbon sink status zoning layer, where the total area of ​​the growth zone was 2.18 km². 2 The stable zone is 3.74 km. 2 The attenuation zone is 0.58 km. 2 Density peak clustering algorithm was used to perform spatial clustering in the carbon sink growth, carbon sink stability, and carbon sink decline regions. For example, in the growth region, density peak clustering analysis was applied to construct the sample feature space by combining the pixel carbon increment value and Euclidean spatial distance. Local density and minimum distance indices were calculated, and eight cluster centers were determined by the truncation distance method. After clustering, three carbon sink hotspots, four transition zones, and one deficit edge zone were identified. The average monthly carbon increment in the hotspots was 0.87 kgC / m³. 2 It is significantly higher than the regional average of 0.34 kgC / m 2 These results were used to identify core areas for carbon sequestration growth in forest regions, providing a spatiotemporal basis for precise carbon sequestration enhancement and carbon compensation strategies.

[0039] The above description is merely a preferred embodiment of this application and an explanation of the technical principles employed. Those skilled in the art should understand that the scope of the invention involved in this application is not limited to technical solutions formed by specific combinations of the above-described technical features, but should also cover other technical solutions formed by arbitrary combinations of the above-described technical features or their equivalents without departing from the concept of this application. For example, technical solutions formed by substituting the above features with (but not limited to) technical features with similar functions disclosed in this application.

Claims

1. A method for predicting forest carbon sink change trends based on multi-source heterogeneous data analysis, characterized in that, include: Raw multi-source remote sensing images are acquired using drones, and spatial alignment is generated by synchronizing with sensor calibration data using timestamps and performing geometric correction and pixel-level registration. Based on the NDVI index layers in the registered RGB and multispectral images, generate a binary forest canopy mask image of the canopy boundary; Based on the registered LiDAR point cloud data, extract the point cloud subset of the corresponding region of the canopy mask, generate the main direction projection image, obtain the tree segmentation outline and center point coordinate label, and construct the three-dimensional structural parameter set of the tree canopy. Based on RGB color histograms, LBP and GLCM texture features, multispectral features, and combined with crown width and height ratio data, tree species and age are identified. Based on the three-dimensional structural parameters of trees, tree species and age, the carbon density per unit crown width of trees is predicted, and the static carbon storage of each tree is calculated. Based on the acquisition of canopy surface temperature and identification of transpiration anomalies by thermal infrared images, combined with leaf water potential status indicators, tree species and age, three-dimensional structural parameters of the canopy and daily carbon absorption data, the carbon absorption rate of the canopy is predicted. Based on the static carbon storage and carbon absorption rate of each tree, a dynamic carbon storage time series layer is constructed. Combined with carbon storage difference monitoring and carbon sink status zoning layer, carbon sink hotspots and carbon deficit areas are identified.

2. The method according to claim 1, characterized in that, The process of acquiring raw multi-source remote sensing images using drones and then performing geometric correction and pixel-level registration to generate spatial alignment by synchronizing with sensor calibration data using timestamps includes: Based on the coordinate range information of the forest area to be monitored, a drone equipped with a high-resolution RGB camera, a multispectral camera, a hyperspectral imager, and a thermal infrared sensor was used to acquire original multi-source remote sensing image data covering the entire target forest area, including RGB images, multispectral images, and thermal infrared images. Time-stamp synchronization and sensor calibration data were used to perform time alignment and geometric correction on images from different sensors to obtain a spatially consistent remote sensing image dataset. An image registration algorithm was used to perform pixel-level joint registration of RGB, multispectral, and thermal infrared images, and nearest neighbor interpolation was used for resampling to obtain a spatially aligned multi-source sensing image raster set.

3. The method according to claim 1, characterized in that, The step of generating a binary forest canopy mask image of the canopy boundary based on the NDVI index layer in the registered RGB and multispectral images includes: Based on the NDVI index layers in the registered RGB and multispectral images, the Sobel operator is used to perform edge gradient enhancement processing on the RGB images to extract high-frequency edge regions. Through a region growing algorithm, pixels with NDVI greater than a preset NDVI threshold are used as seed regions, and a growth threshold is set to expand adjacent pixels into 8-neighborhoods. By performing logical operations on the binary image of the Sobel edge detection result and the mask image generated by NDVI region growing at the pixel level, the output is set to 1 when the values ​​of the two layers are both 1 at the same pixel position, and 0 otherwise, a binary forest canopy mask image containing only the complete canopy boundary is obtained, and the minimum bounding rectangle boundary information of each connected component is preserved.

4. The method according to claim 1, characterized in that, The process involves extracting a subset of point clouds corresponding to the canopy mask region from the registered LiDAR point cloud data, generating a principal direction projection image, obtaining tree segmentation contours and center point coordinate labels, and constructing a three-dimensional structural parameter set for the tree canopy, including: Based on the registered LiDAR point cloud data, a subset of point clouds corresponding to the canopy mask region is extracted. Principal component analysis (PCA) is used to estimate the principal orientation of the point clouds within each connected canopy region, obtaining the principal axis direction vector of the point clouds in that region. A two-dimensional local orthographic projection coordinate system is constructed using the principal axis direction, and the point clouds are projected onto the principal orientation plane along the principal axis direction. Based on a preset rasterization resolution, the principal orientation projection image of the region is obtained. The Euclidean distance field image of each connected region is calculated through distance transformation, and local maxima are extracted as initial seed points for the watershed algorithm. The watershed segmentation algorithm is used to process the Euclidean distance field image, based on a preset damping threshold. The correspondence between the label layer and the original image coordinates is obtained, and the segmentation outline, polygon boundary coordinates, and cropped image blocks of each tree are obtained. The corresponding coordinate labels are recorded to obtain the spatial coordinates of the trees. The corresponding coordinate labels are the center point coordinates of the cropped image blocks. The crown width and crown area of ​​each tree are calculated using a boundary tracking algorithm. The crown height is extracted by combining the difference between the original digital surface model and the digital terrain model. A three-dimensional structural parameter set for each crown is constructed and saved to the forest area monitoring database. The three-dimensional structural parameters include crown area, crown height, crown width, crown index, and height ratio. The height ratio is the ratio of crown width to tree height.

5. The method according to claim 1, characterized in that, The method of identifying tree species and age based on RGB color histograms, LBP and GLCM texture features, multispectral features, and combined with crown width and height ratio data includes: Color histograms are extracted from the RGB images of each tree to obtain color distribution feature vectors. Texture features are calculated based on LBP and GLCM, and energy, contrast, and entropy are extracted. Multispectral features, including NDVI, NDRE, and EVI, are extracted from multispectral images and saved to the forest area monitoring database. Historical data on color distribution feature vectors, texture features, multispectral features, crown width, and height ratio are obtained from the forest area monitoring database. Tree species and ages are labeled, and a tree species identification model is constructed using the random forest algorithm. Based on the color distribution feature vectors, texture features, multispectral features, crown width, and height ratio of the real-time acquired tree data, the tree species and age are identified using the tree species identification model.

6. The method according to claim 1, characterized in that, The method predicts the carbon density per unit crown width of trees based on their three-dimensional structural parameters, tree species, and age, and calculates the static carbon storage of each tree, including: By using the forest area monitoring database, the three-dimensional structural parameters, tree species, and tree age of trees are obtained, and the corresponding unit crown carbon density of trees is labeled. A recurrent neural network is used to train the model and construct a unit crown carbon density prediction model. Based on the real-time acquired three-dimensional structural parameters, tree species, and tree age of trees, the unit crown carbon density prediction model is used to predict the current unit crown carbon density of trees. The static carbon storage of each tree is calculated by multiplying the crown area by the unit crown carbon density and combining the area data extracted from the segmentation boundary of a single tree. By traversing the spatial location and carbon storage of each tree, a carbon storage raster layer is constructed and reprojected onto the spatial coordinates of the trees to obtain a spatial distribution map of carbon storage in the forest area.

7. The method according to claim 1, characterized in that, The method of obtaining canopy surface temperature and identifying transpiration anomalies based on thermal infrared imaging, combined with leaf water potential state indicators, tree species and age, three-dimensional structural parameters of the canopy, and daily carbon uptake data, to predict the carbon uptake rate of the canopy includes: The canopy surface temperature value is obtained based on the pixel region corresponding to the spatial location of each tree in the thermal infrared image; By comparing the temperature value of each tree with the average temperature of the same type of trees within a preset distance, it is determined whether there is an abnormality in transpiration in the current tree, and the location of the abnormal tree is marked. The average pixel temperature of the canopy of each tree with abnormal transpiration is calculated by the regional mean method, and the leaf surface water potential index of the canopy of each tree with abnormal transpiration is calculated. Combining the tree species, tree age, three-dimensional structural parameters of the canopy, as well as the daily carbon absorption data measured by the flux observation tower and the ground, a recurrent neural network is used to train the model and construct a tree carbon absorption rate prediction model to predict the carbon absorption rate of the canopy of each tree.

8. The method according to claim 7, characterized in that, The calculation of the average pixel temperature of the canopy of each tree with abnormal transpiration using the regional mean method, and the calculation of the leaf surface water potential index of the canopy of each tree with abnormal transpiration, include: Based on the average pixel temperature of the canopy of each tree with abnormal transpiration, combined with the canopy temperature and transpiration model... Calculate latent heat flux ,in, Net radiation from the leaf surface is calculated using measured meteorological parameters or obtained from a drone's light sensor. The sensible heat flux was calculated using the temperature difference between the canopy and the air, as well as aerodynamic drag parameters. Soil heat flux; net canopy radiation flux was calculated using light intensity measurements and reflectance, and then expressed as a leaf water potential index. Calculate the leaf surface water potential state index of the canopy of each tree with abnormal transpiration. ,in, Leaf conductance was determined using tree species characteristics and temperature and humidity response curves. Δ is the pressure constant, determined by air pressure, and Δ is the slope of the saturated water vapor pressure curve, calculated based on air temperature. The average temperature of the canopy. The ambient air temperature is obtained from ground weather stations or drone sensors.

9. The method according to claim 1, characterized in that, The method constructs a dynamic carbon storage time series layer based on the static carbon storage and carbon absorption rate of each tree, and combines carbon storage difference monitoring and carbon sink status zoning layer to identify carbon sink hotspots and carbon deficit areas, including: Based on the static carbon storage and carbon absorption rate of each tree, the carbon increments in preset time units are sequentially superimposed onto the static carbon density map to construct a dynamic carbon storage time series layer group for preset time units. A pixel-based back-difference calculation method is used to monitor the carbon storage difference at different time periods. By setting change threshold intervals, the carbon sink growth, carbon sink stability, or carbon sink decline status within a region is determined, resulting in a carbon sink status zoning layer. Based on the carbon sink status zoning layer, a density peak clustering algorithm is applied to spatially cluster the carbon sink growth, carbon sink stability, and carbon sink decline regions respectively. A sample feature space is constructed using the Euclidean distance matrix and carbon change density values. Local density and minimum distance indices are calculated, and cluster centers are determined using the truncation distance method. Clustering operations are performed to obtain distribution maps of carbon sink hotspots, transition zones, and carbon deficit zones. The total area and monthly average carbon increment of each region are statistically analyzed.

Images