Multi-source remote sensing data forestry carbon sink metering and monitoring method and device and storage medium

Abstract

The invention provides a multi-source remote sensing data forestry carbon sink metering and monitoring method and device and a storage medium, and the method comprises the steps: 1, designing a monitoring scheme; step 2, processing and calculating data of manual scale measurement; step 3, unmanned aerial vehicle data processing and calculation; step 4, satellite data processing and calculation; and 5, based on the steps 2, 3 and 4, calculating to obtain the final target carbon reserve and carbon sink quantity, thereby completing the monitoring work of carbon sink metering. According to the method, the reliability of forest land carbon reserve monitoring is improved, the monitoring cost is reduced, the data processing efficiency is improved, the monitoring period is shortened, various technical means are combined, the comprehensive utilization rate of data is improved, an innovative stratified sampling method is provided, and the method is suitable for popularization and application. The carbon sink metering and monitoring precision is improved, efficient carbon layer division is achieved, and the overall operation efficiency is improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of environmental monitoring and assessment, and in particular relates to a multi-source remote sensing data forestry carbon sink measurement monitoring method, equipment and storage medium. Background Art

[0002] Today, global climate change and carbon emissions are becoming increasingly prominent. As an important carbon sink, the accurate calculation of forest carbon storage is of great significance for global response to climate change. Traditional forest carbon storage calculation methods rely on sample site surveys, single data sources, and systematic sampling methods, which are not sufficient to provide high-precision and efficient monitoring results.

[0003] Moreover, the current monitoring of forest carbon storage mostly uses methods such as manual measurement and remote sensing technology, but these methods have problems such as low accuracy, high cost and low efficiency. Summary of the invention

[0004] In view of the above problems, the present invention provides a forestry carbon sink measurement and monitoring method based on multi-source remote sensing data.

[0005] The multi-source remote sensing data forestry carbon sink measurement and monitoring method provided by the present invention comprises:

[0006] Step 1: Design of monitoring plan;

[0007] Step 2: Manual measurement data processing and calculation;

[0008] Step 3: UAV data processing and calculation;

[0009] Step 4: Satellite data processing and calculation;

[0010] Step 5: Based on the steps 2, 3, and 4, the final target carbon storage and carbon sink are calculated, thereby completing the monitoring of carbon sink measurement.

[0011] in,

[0012] The step 1 comprises:

[0013] Step 1.1, small class vector data management, including: obtaining vector data, i.e. small class vector data, entering, storing, managing, visualizing and applying shape data files composed of the vector data, including obtaining a vector data small class attribute table;

[0014] Step 1.2, carbon layer factor screening and carbon layer division;

[0015] Step 1.3: Calculate the number of carbon layer sample plots and systematically randomly distribute the sample plots, determine the sample plots, and generate the sample site locations.

[0016] further,

[0017] The step 1.2 comprises:

[0018] The factor information in the vector data small class attribute table: tree species, age group, slope and the statistical information of the factors: data distribution map, factor maximum value, factor average value are used for small class classification. Specifically, the factors also include: county, township, village, forest farm, forest team, forest class, land type, forest species, terrain, slope direction, slope position, dominant tree species, age group, and canopy density. Different factors correspond to their own enumeration value groups one by one.

[0019] The carbon layer division includes: selecting enumeration values ​​from different enumeration value groups in turn for combination, selecting one enumeration value from each enumeration value group for each combination, until all the enumeration values ​​are involved in the combination, and an enumeration value combination obtained by combination corresponds to a carbon layer. When combining, if N factors are involved, the N factors correspond to enumeration value group 1, ..., enumeration value group N one by one, and the enumeration value group 1, ..., enumeration value group N respectively include n1, ..., nn enumeration values, then this step divides n1×...×nn carbon layers, and outputs the carbon layer division result, i.e., the carbon layer result, which includes the carbon layer code, the carbon layer name, the carbon layer characteristics, i.e., the factor information and the statistical information of the factor, the number of small classes contained in the carbon layer, the carbon layer area and other carbon layer data, wherein N, n1, ..., nn are integers greater than 0.

[0020] In addition, the plots within the target boundary need to be divided into different levels according to different stratification factors.

[0021] further,

[0022] The step 1.3 comprises:

[0023] According to the carbon layer division results, sampling formulas and statistical data in vector data, the number of sample plots and the area of ​​small classes in each carbon layer are calculated, and the sample plot layout is completed according to the method of random starting point and systematic layout. The sample plots are determined, and the sample plot distribution map, UAV flight area map and sample plot information table are output, including the sample plot number and the longitude and latitude coordinates of the sample plot location, that is, the sample plot location. According to the selected target tree species, for each target tree species, at least 60 individual trees are randomly selected as modeling samples at the target location. The sample tree diameter at breast height is not less than 2 cm, and they are divided according to the diameter class interval of 2 cm to 3 cm. There are no less than 2 to 3 sample trees in each diameter class to cover all diameter classes.

[0024] further,

[0025] The step 2 comprises:

[0026] 2.2.1. Calculation of biomass and carbon storage of sample plots:

[0027] Based on the data obtained from the sample plot monitoring: tree species, breast diameter, tree height, crown width, coordinates, combined with the biomass equation method, the biomass of individual trees in the sample plot is calculated, and the biomass at the sample plot level is accumulated. Then, combined with the biomass carbon content, the biomass carbon storage at the sample plot level in each year and each target carbon layer is calculated;

[0028] 2.2.2. Calculation of estimated value of average carbon storage per unit area of ​​carbon layer and statistical estimation of target carbon storage, including the following steps:

[0029] Step 1: Calculate the estimated value of the average biomass carbon storage per unit area for each carbon layer:

[0030]

[0031] Where:

[0032] is the average biomass carbon stock per unit area of ​​the i-th target carbon layer in the T-th year, in tons of carbon per hectare, T is an integer, and i is an integer greater than 0;

[0033] is the biomass carbon storage per unit area of ​​the i-th target carbon layer plot p in the T-th year, in tons of carbon per hectare, The biomass carbon stock per unit area of ​​the plot is obtained by dividing the biomass carbon stock at the plot level of the i-th target carbon layer in the T-th year obtained in step 2.2.1 by the area of ​​the plot p, in tons of carbon per hectare;

[0034] In formula (1), n i is the number of monitoring plots of the i-th target carbon layer, which is a positive integer and dimensionless; m represents the data directly obtained through field measurements of the plots or the value calculated based on the data and dimensionless;

[0035] Step 2: Calculate the estimated value of the average biomass carbon stock per unit area of ​​the target population:

[0036]

[0037] Where:

[0038] In year T, the target is the average biomass carbon stock per unit area within the target area boundary, in tons of carbon per hectare;

[0039] w i is the area weight of the i-th target carbon layer within the target boundary, dimensionless, w i =A i / A, where A is the total area of ​​each carbon layer within the target boundary, in hectares, A i is the area of ​​the i-th target carbon layer, in hectares;

[0040] Step 3: Calculate the total biomass carbon stock within the target boundary in year T:

[0041]

[0042] Where:

[0043] It is the biomass carbon stock within the target boundary in year T, in tons of carbon.

[0044] further,

[0045] The step 3 comprises:

[0046] Data Registration:

[0047] The single tree samples of the sample plot are matched according to the single tree position distance error and the single tree height accuracy. That is, for the sample plot, the single tree monitored by the drone is matched with the actual single tree. At the same time, the biomass of the single tree, that is, the reference aboveground biomass of the single tree, is obtained from the step 2.2.1. Then, the single tree laser radar characteristic parameters are extracted from the laser radar data obtained by scanning the target area by the three-dimensional laser radar equipment on the drone: tree height, percentile tree height, crown width, laser radar biomass index and percentile tree height data.

[0048] Specifically, the matching is to match the data obtained by manual measurement with the data obtained by scanning or monitoring by the three-dimensional laser radar equipment on the drone, find the same tree in the same position, and use the biomass of the single tree as the prior value of the biomass monitored by the drone to form a first data set, which includes at least 20 data.

[0049] further,

[0050] The step 3 includes: 3.3, constructing a laser radar biomass inversion model:

[0051] 3.3.1.1. Construct a linear regression model:

[0052] Randomly select 3 / 4 samples in the prior biomass data set as the training set. In this step, the prior biomass data set is the first data set, and the biomass model is constructed using the linear regression method:

[0053]

[0054] Where:

[0055] is the measured whole plant or aboveground biomass of the kth tree, in kilograms. In this step, the measured whole plant or aboveground biomass is the biomass measured by the manual measuring method;

[0056] is the monitored whole plant or aboveground biomass of the kth tree, in kilograms. In this step, the monitored whole plant or aboveground biomass is the laser radar whole plant or aboveground biomass;

[0057] is the monitoring characteristic parameter of the kth tree, in meters or dimensionless. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter;

[0058] SSE is the residual sum of squares, SSE = ∑ k (b1 k -b2 k ) 2 , where b1 k is the first biomass of the kth tree, b2 k is the second biomass of the kth tree. In this step, b1 k for b2 k for

[0059] is the assumed linear function between the measured whole plant or aboveground biomass and the monitored characteristic parameters of the sample plot; the parameters are solved by minimizing the SSE in the survey plot data, and the linear function after solving the parameters is recorded as The fitted function is used to calculate the biomass of the whole plant (or aboveground).

[0060] l represents data directly obtained through monitoring technology or a value calculated based on the data. In this step, the monitoring technology is laser radar monitoring;

[0061] m represents the data directly obtained through field measurements or the value calculated based on the data;

[0062] k represents a single sample tree and is an integer greater than 0.

[0063] Based on formula (4), a model equation for the relationship between monitoring characteristic parameters and monitoring whole plant or aboveground biomass was established, where the univariate linear regression model is:

[0064]

[0065] Where:

[0066] a1 and b1 are model parameters, dimensionless;

[0067] The binary linear regression model is:

[0068]

[0069] Where:

[0070] a2, b2 and c2 are model parameters, dimensionless;

[0071] is the i1th monitoring characteristic parameter of the kth tree. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter, i1 is an integer greater than 0, and the unit is meter or dimensionless;

[0072] is the j1th monitoring characteristic parameter of the kth tree. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter, j1 is an integer greater than 0, and j1≠i1, and the unit is meter or dimensionless.

[0073] Model verification:

[0074]

[0075] Where:

[0076] R 2 is the coefficient of determination of the model, dimensionless;

[0077] b1 k The average value, in kilograms,

[0078] In this step, b1 k for b2 k for Calculate the coefficient of determination of equations (5) and (6) and compare them, and select R 2 The equation with the best value, i.e. the largest value, is taken as the final biomass prediction equation.

[0079] 3.3.1.2. Construct a machine learning algorithm based on a random forest regression model: The machine learning algorithm in this step uses the Bootstrap resampling technology to extract n Bootstrap samples from the first data set to form a data set, where n is an integer greater than 1, and performs decision tree modeling on each Bootstrap sample to obtain n decision trees: decision tree-1, ..., decision tree-n, and then the prediction results of the n decision trees, i.e., result-1, ..., result-n, are used to obtain the final prediction result through majority voting or an average mechanism. The prediction result is the biomass. When constructing the model, 2 / 3 of the samples in the first data set are randomly selected as a training set, and the randomForest function under the randomForest data package of the R language is used to implement the construction of the laser radar biomass model, including the steps of:

[0080] A) Selection of laser radar characteristic parameters:

[0081] When using the random forest method to build a biomass model,

[0082] Aa), the random forest algorithm automatically evaluates the influence of each lidar feature parameter in the regression process based on the model mean square error increment IncMSE and the influence on the purity of the model tree node IncNodePurity;

[0083] Ab), select parameters TreeHeight, H99, H98, H90, H97 to construct the lidar biomass model, where the parameters mean: TreeHeight—tree height, H99—99th percentile tree height, H98—98th percentile tree height, H90—90th percentile tree height, H97—97th percentile tree height,

[0084] B) Selection of the number of regression trees:

[0085] Ba), use the parameter ntree to define the number of decision trees, and the parameter mtry to define the number of randomly selected features;

[0086] Bb), ntree = 200 as the number of regression trees;

[0087] Bc) mtry is not assigned a value, and the system default value of 1 / 3 of the number of input variables is used. If the number of characteristic parameters is less than 3, the value of mtry is 1.

[0088] Model verification:

[0089] According to formula (7), if the model meets the minimum requirement, the model accuracy test index is the determination coefficient R 2 If the minimum value of is 0.7, the model passes the test and is used for biomass calculation.

[0090] further,

[0091] The step 3 also includes: 3.4, calculating the target biomass and carbon storage: the estimated value of the average biomass carbon storage per unit area of ​​each carbon layer obtained by airborne laser radar monitoring is combined with the area of ​​the target area to calculate the biomass carbon storage within the target boundary, including the steps of:

[0092] 3.4.1. Calculate the biomass and carbon storage of individual trees in the target area covered by the drone:

[0093]

[0094] In the above formula:

[0095] is the LiDAR biomass carbon stock of the s2th tree of the j2th tree species in the i2th target carbon layer in the T2th year, in tons of carbon, T2 is an integer, and i2, j2 and s2 are integers greater than 0;

[0096] is the total plant or aboveground biomass of the s2th tree species in the i2th target carbon layer in the t2th year, in kilograms, which is obtained by bringing the characteristic parameter data monitored by the airborne lidar into the lidar biomass inversion model. The lidar biomass inversion model uses the machine learning algorithm based on the linear regression model and the random forest regression model. 2 the larger one;

[0097] CF j2 is the carbon content of the biomass of the j2-type tree species, dimensionless,

[0098] 3.4.2. Calculate the estimated value of the average carbon storage per unit area of ​​the carbon layer and the target carbon storage:

[0099] The first step is to calculate the biomass carbon storage per unit area of ​​the j2 type of tree species:

[0100]

[0101] Where:

[0102] is the LiDAR biomass carbon stock per unit area of ​​the j2th tree species in the i2th target carbon layer in year T2, in tons of carbon;

[0103] A i2,j2,t2 is the forest area of ​​the j2-type tree species in the i2-target carbon layer in year T2, in hectares;

[0104] Step 2: Calculate the estimated value of the average forest biomass carbon storage per unit area in each carbon layer:

[0105]

[0106] Where:

[0107] is the average per unit area lidar biomass carbon stock of the i2th target carbon layer in year T2, in tons of carbon;

[0108] nn i2 is the number of tree species in the i2th target carbon layer, a positive integer, dimensionless;

[0109] Step 3: Combined with the area weight of the i2th target carbon layer within the target boundary, calculate the estimated value of the target overall average biomass carbon storage per unit area:

[0110]

[0111] Where:

[0112] is the average biomass carbon stock per unit area within the target boundary in year T2, in tons of carbon per hectare;

[0113] w i2 =A i2 / A, where, A i2 is the area of ​​the i2th target carbon layer, in hectares;

[0114] Step 4: Calculate the total biomass carbon stock within the target boundary based on the total target area:

[0115]

[0116] Where:

[0117] It is the biomass carbon stock obtained through drone monitoring within the target boundary in year T2, in tons of carbon.

[0118] further,

[0119] The step 4 comprises:

[0120] Data Registration:

[0121] According to the position distance error of the sample plot, that is, the distance error between the position of the UAV laser radar sample plot and the position of the satellite laser radar sample plot, the first sample plot measured by the satellite laser radar, that is, the satellite remote sensing data segmentation, and the second sample plot sample obtained by the UAV laser radar measurement data segmentation are matched. At the same time, the biomass of the individual trees obtained in step 3.4.1 is accumulated as the sample plot biomass obtained by laser radar segmentation, that is, the reference aboveground biomass of the sample plot.

[0122] Among them, the matching is to match the data obtained by the drone lidar monitoring with the satellite remote sensing data, find the same sample at the same position, and use the biomass of the same sample calculated based on the data obtained by the drone lidar monitoring, recorded as the lidar biomass, as the biomass monitored by the satellite lidar, that is, the prior value of the satellite remote sensing biomass, to form a second data set. The second data set includes at least 20 data. The data obtained based on the drone lidar monitoring is processed into the lidar characteristic parameters, and the satellite remote sensing data is processed into the satellite remote sensing characteristic parameters.

[0123] further,

[0124] The step 4 also includes: 4.3.1.1, constructing a linear regression model,

[0125] The difference between step 4.3.1.1 and step 3.3.1.1 is that the second data set is used as a priori biomass data set;

[0126] In formula (5) and (6), Replace with is the satellite remote sensing biomass of the kkth plot, in kilograms, where kk represents the plot and is an integer greater than 0;

[0127] In formula (5) Replace with is the satellite remote sensing data of the kkth sample plot, in kilograms;

[0128] In formula (6) and Replace with and is the i3th satellite remote sensing characteristic parameter of the kkth sample plot, i3 = 1, 2, 3, ..., in meters or dimensionless units, is the j3th satellite remote sensing characteristic parameter of the kkth sample plot, j3=1, 2, 3, ... (j3≠i3), in meters or dimensionless;

[0129] s represents the data directly obtained through satellite remote sensing or the value calculated based on the data, dimensionless,

[0130] The step 4 also includes: 4.3.1.2, a machine learning algorithm based on a random forest regression model,

[0131] The difference between step 4.3.1.2 and step 3.3.1.2 is:

[0132] Taking the second data set as the prior biomass data set, the randomForest function under the randomForest data package of R language is used to realize the construction of the satellite remote sensing biomass model;

[0133] Replace the lidar characteristic parameters of the drone coverage area in step 3.3.1.2 with satellite remote sensing characteristic parameters.

[0134] further,

[0135] The step 4 further includes: 4.4, calculating the target biomass and carbon stock: based on the characteristic parameter data of satellite remote sensing, combined with the satellite remote sensing biomass inversion model constructed in step 4.3, calculating the biomass of the trees in the target area, and then calculating the biomass carbon stock within the target boundary in combination with the biomass carbon content, including:

[0136] 4.4.1. Calculate the tree biomass and carbon storage in the satellite coverage area, i.e. the target area, including:

[0137] First, based on the characteristic parameter data of satellite remote sensing, combined with the satellite remote sensing biomass inversion model constructed in step 4.3, the biomass of each sample plot of the i4th target carbon layer in the target area is calculated to form a satellite remote sensing sample plot biomass data set of the i4th target carbon layer in the T4th year, and then the sample plot biomass is converted into sample plot biomass carbon storage in combination with the biomass carbon content, T4 is an integer, i4 and q4 are both integers greater than 0, wherein the satellite remote sensing biomass inversion model uses the R in the machine learning algorithm based on the linear regression model and the random forest regression model 2 The larger one,

[0138]

[0139] In formula (13):

[0140] is the satellite remote sensing biomass carbon stock of the q4th plot of the i4th target carbon layer in the T4th year, in tons of carbon;

[0141] is the satellite remote sensing biomass of the q4th plot of the i4th target carbon layer in the T4th year, in kilograms, which is obtained by bringing the characteristic parameter data of satellite remote sensing monitoring into the satellite remote sensing biomass inversion model;

[0142] CF is the carbon content of biomass, dimensionless,

[0143] 4.4.2. Calculation of carbon reserves in carbon layer and target carbon reserves:

[0144] Step 1: Calculate the biomass carbon storage in each carbon layer:

[0145]

[0146] Where:

[0147] is the satellite remote sensing biomass carbon stock of the i4th target carbon layer in the T4th year, in tons of carbon,

[0148] Step 2: Calculate the biomass carbon stock within the target boundary:

[0149]

[0150] It is the satellite remote sensing biomass carbon stock within the target boundary in year T4, in tons of carbon.

[0151] further,

[0152] The step 5 comprises:

[0153] The final target carbon storage and annual change, i.e. carbon sink, are calculated by weighted average:

[0154]

[0155] In the above formula:

[0156] C Biomass,T5 is the final biomass carbon stock within the target boundary in year T5, in tons of carbon;

[0157]

[0158] In the above formula:

[0159] ΔC Biomass,T5 is the annual change in the target biomass carbon stock in year T5 within the target boundary, in tons of carbon dioxide equivalent per year;

[0160] t1, t2 – the t1th year and t2th year after the start of monitoring, expressed in years. t1, t2 and T5 are all integers, and t1≤T5≤t2.

[0161] The present invention also provides a computer device, which includes a memory, a first processor, and a first computer program stored in the memory and executable on the first processor. When the first computer program is executed by the first processor, the multi-source remote sensing data forestry carbon sink measurement and monitoring method is implemented.

[0162] The present invention also provides a computer-readable storage medium, which is used to store a second computer program. The second computer program can be executed by at least one second processor so that the at least one second processor executes the multi-source remote sensing data forestry carbon sink measurement and monitoring method.

[0163] Compared with the prior art, the multi-source remote sensing data forestry carbon sink measurement and monitoring method provided by the present invention has the following advantages and beneficial effects:

[0164] 1. Improved the reliability of forest carbon storage monitoring.

[0165] 2. Reduce monitoring costs, especially in forest areas that are difficult to reach manually.

[0166] 3. Improved data processing efficiency and shortened monitoring cycle.

[0167] 4. Combining multiple technical means to improve the comprehensive utilization rate of data.

[0168] 5. Through the collection, storage, analysis, integration and optimization of vector data, an innovative stratified sampling method is proposed to improve the accuracy of carbon sink measurement and monitoring. This method analyzes multiple factors in the forest sub-class attribute table to accurately divide the carbon layer, ensuring the scientific layout and effective monitoring of the sampling points. Since the factors contained in the sub-class attribute table can reflect the ecological characteristics of the forest area in detail, this method uses these key information to divide the target covered plots into multiple levels, thereby ensuring that the selection of samples can more accurately represent the overall forest area. This hierarchical sampling strategy not only enhances the representativeness of the sample, but also significantly improves the characterization ability and generalization of the monitoring results, allowing monitoring program compilers to better understand the sub-class attribute characteristics of their target forest area, reduce sampling errors, ensure that the accuracy of carbon sink measurement is significantly improved, and formulate more effective carbon layer division strategies based on this information, achieving efficient carbon layer division and improving overall work efficiency.

[0169] Other features and advantages of the present invention will be described in the following description or understood by practicing the present invention. The purpose and other advantages of the present invention can be realized and obtained by the structures pointed out in the description, claims and drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0170] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative work.

[0171] Figure 1 A basic flow chart of a method for measuring and monitoring forestry carbon sinks using multi-source remote sensing data according to an embodiment of the present invention is shown;

[0172] Figure 2 A schematic diagram of a decision tree of a random forest regression model in a forestry carbon sink measurement and monitoring method using multi-source remote sensing data according to an embodiment of the present invention is shown;

[0173] Figure 3 The IncMSE value in the random forest regression model in the multi-source remote sensing data forestry carbon sink measurement and monitoring method according to an embodiment of the present invention is shown;

[0174] Figure 4 The IncNodePurity value in the random forest regression model in the forestry carbon sink measurement and monitoring method using multi-source remote sensing data according to an embodiment of the present invention is shown;

[0175] Figure 5A schematic diagram of the relationship between ntree and error in a random forest regression model in a forestry carbon sink measurement and monitoring method using multi-source remote sensing data according to an embodiment of the present invention is shown. DETAILED DESCRIPTION

[0176] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0177] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by technicians in the technical field to which this application belongs; the terms used in the specification of the application herein are only for the purpose of describing specific embodiments and are not intended to limit this application; the terms "including" and "having" in the specification and claims of this application and the above-mentioned figure descriptions and any variations thereof are intended to cover non-exclusive inclusions. The terms "first", "second", "third", etc. in the specification and claims of this application or the above-mentioned figures are used to distinguish different objects, not to describe a specific order or a primary and secondary relationship. "Multiple" appearing in this application refers to more than two (including two).

[0178] Reference to "embodiments" herein means that a particular feature, structure, or characteristic described in conjunction with the embodiments may be included in at least one embodiment of the present application. The appearance of the phrase in various locations in the specification does not necessarily refer to the same embodiment, nor is it an independent or alternative embodiment that is mutually exclusive with other embodiments. It is explicitly and implicitly understood by those skilled in the art that the embodiments described herein may be combined with other embodiments.

[0179] Figure 1 The basic flow chart of the multi-source remote sensing data forestry carbon sink measurement and monitoring method provided by the present invention includes the following steps.

[0180] 1. Monitoring plan design. This step, based on vector data and on-site surveys, deeply analyzes the vegetation type, distribution and quantity of the forest land for the target area, and then completes the carbon layer division. On this basis, the idea of ​​sampling survey is used to complete the sampling and point distribution as well as the selection of single tree samples. This step finally outputs the monitoring plots and drone flight areas, providing a solid foundation for subsequent monitoring work. This step further includes the following steps:

[0181] 1.1. Vector data management for small classes.

[0182] Obtain vector data (i.e. small class vector data), enter, store, manage, visualize and apply shape data files composed of vector data, including obtaining small class attribute tables of vector data. The vector data is provided by the forest management unit.

[0183] 1.2. Carbon layer factor screening and carbon layer division (including small class classification).

[0184] The small class classification is performed using the tree species, age group, slope and other factor information (i.e., carbon layer factor screening) and the factor statistical information (such as data distribution map, factor maximum value, factor average value, etc.) in the vector data small class attribute table. Specifically, the factors may also include: county, township, village, forest farm, forest team, forest class, land type, forest species, terrain, slope direction, slope position, dominant tree species, age group, and canopy density. Among them, different factors correspond to their own enumeration value groups one by one. For example, the three factors of tree species, age group, and slope correspond to tree species enumeration value groups, age group enumeration value groups, and slope enumeration value groups, respectively. The tree species enumeration value group may include enumeration values ​​such as Masson pine and eucalyptus, the age group enumeration value group may include enumeration values ​​such as young forest, middle-aged forest, near-mature forest, and over-mature forest, and the slope enumeration value group may include enumeration values ​​such as 0-4, 16-25, and 26-35 (in degrees). The factors and enumeration values ​​are determined by technical personnel based on specific measurement and monitoring. Then the carbon layer division is performed: enumeration values ​​are selected from different enumeration value groups in turn for combination, and one enumeration value is selected from each enumeration value group each time the combination is performed, until all enumeration values ​​are involved in the combination, and one enumeration value combination obtained by the combination corresponds to a carbon layer. When combining, if N factors are involved, the N factors correspond to enumeration value groups 1, ..., and enumeration value groups N one by one, and enumeration value groups 1, ..., and enumeration value groups N include n1, ..., nn enumeration values ​​respectively, then this step divides n1×...×nn carbon layers, where N, n1, ..., nn are integers greater than 0. If the tree species enumeration value group, age group enumeration value group, and slope enumeration value group contain n1, n2, and n3 enumeration values ​​respectively, then a total of n1×n2×n3 carbon layers are obtained by combination, and the carbon layer division result, i.e., the carbon layer result, is output. The carbon layer results include carbon layer data such as carbon layer code, carbon layer name, carbon layer characteristics (i.e. the aforementioned factor information and factor statistical information), the number of small classes contained in the carbon layer, and the carbon layer area.

[0185] In addition, in order to improve the accuracy of calculating carbon stock changes and reduce the number of monitoring plots under certain accuracy requirements, the plots within the target boundary need to be divided into different layers according to different stratification factors. The divided carbon layers are used to calculate carbon stocks. The site conditions of the land within the target boundary before afforestation (such as soil type, slope aspect, altitude, etc.) and the proposed afforestation time, afforestation tree species, afforestation density and other factors are considered to divide the carbon layers. Afforestation plots with no significant differences are divided into the same carbon layer. Adjustments and determinations are made based on the actual situation of afforestation activities.

[0186] 1.3. Calculate the number of carbon layer sample plots and randomly distribute them systematically, determine the sample plots, and generate the sample site locations.

[0187] According to the results of carbon layer division, sampling formulas such as the calculation formula for the number of monitoring sample plots and statistical data in vector data, the number of sample plots for each carbon layer, small class area and other variables are calculated, and the sample plot layout is completed according to the method of random starting point and system layout (i.e., system random layout), the sample plot is determined, and the sample plot distribution map, drone flight area map and sample plot information table (including sample plot number, sample plot location latitude and longitude coordinates, i.e., sample site location) are output. According to the selected target tree species, no less than 60 single trees are randomly selected for each target tree species at the target location as modeling samples. The sample single tree has a breast diameter of no less than 2 cm, and is divided into diameter classes of 2 cm to 3 cm. There are no less than 2 to 3 sample trees in each diameter class, and all diameter classes are covered as much as possible. All selected sample single trees should be easier to segment from the lidar point cloud data.

[0188] 2. Manual measurement data processing and calculation. Including steps:

[0189] 2.1、Manual measurement data acquisition and processing, this step includes:

[0190] 2.1.1. Data acquisition.

[0191] The traditional manual measurement method is used to measure and record the individual trees in the sample plot, obtain the basic information of the trees in the sample plot, such as tree species, breast diameter, height and other data, and collect the aforementioned factors of the sample plot.

[0192] 2.1.2. Data processing.

[0193] Enter the acquired data, bind the data to the small class, complete the data integrity and compliance check and adjustment, and use it for subsequent data calculations. Data integrity check and adjustment refers to determining whether the acquired data, such as plot number, tree number, plot location, tree height, and DBH are missing. If missing, continue to acquire data to supplement until the required data is complete or sufficient. Compliance check and adjustment is to check and adjust whether the data format meets the system requirements, such as whether the coordinate system of the forest location uses the same coordinate system, and if not, adjust it to the same coordinate system.

[0194] 2.2. Calculation of target carbon storage.

[0195] The calculation of carbon storage based on the manual measurement method is usually to obtain the carbon storage of the sample plots measured in the target area, and then calculate the estimated value of the average annual biomass carbon storage per unit area of ​​the carbon layer, and finally calculate the biomass carbon storage within the target boundary. The specific steps include:

[0196] 2.2.1. Calculation of biomass and carbon storage of sample plots.

[0197] Based on the data of tree species, breast diameter, tree height, crown width, coordinates, etc. obtained from sample plot monitoring, combined with the biomass equation method, the biomass of individual trees in the sample plot is calculated, and the biomass at the sample plot level is accumulated. Then, combined with the biomass carbon content, the biomass carbon storage at the sample plot level for each year and each target carbon layer is calculated.

[0198] 2.2.2. Calculation of estimated value of average carbon storage per unit area of ​​carbon layer and statistical estimation of target carbon storage. Including steps:

[0199] Step 1: Calculate the estimated value of the average biomass carbon storage per unit area for each carbon layer:

[0200]

[0201] Where:

[0202] is the average biomass carbon stock per unit area of ​​the i-th target carbon layer in the T-th year, in tons of carbon per hectare, T is an integer, and i is an integer greater than 0;

[0203] is the biomass carbon storage per unit area of ​​the i-th target carbon layer plot p in the T-th year, in tons of carbon per hectare, The biomass carbon stock per unit area of ​​the plot is obtained by dividing the biomass carbon stock at the plot level of the i-th target carbon layer in the T-th year obtained in the above step 2.2.1 by the area of ​​the plot p, in tons of carbon per hectare;

[0204] In formula (1), n i is the number of monitoring plots for the i-th target carbon layer, which is a positive integer and dimensionless; m represents the data directly obtained through field measurements of the plots or the value calculated based on the data and is dimensionless.

[0205] Step 2: Calculate the estimated value of the average biomass carbon stock per unit area of ​​the target population:

[0206]

[0207] Where:

[0208] is the average biomass carbon stock per unit area within the target (i.e. target area) boundary in year T, in tons of carbon per hectare;

[0209] w i is the area weight of the i-th target carbon layer within the target boundary, dimensionless, w i =A i / A, where A is the total area of ​​each carbon layer within the target boundary, in hectares (hm2 ), A i is the area of ​​the i-th target carbon layer, in hectares (hm 2 );

[0210] Step 3: Calculate the total biomass carbon stock within the target boundary in year T:

[0211]

[0212] Where:

[0213] It is the biomass carbon stock within the target boundary in year T, in tons of carbon.

[0214] 3. UAV data processing and calculation. Including:

[0215] 3.1. UAV data acquisition and processing. This step includes:

[0216] 3.1.1. Data acquisition.

[0217] The data set of spatial points obtained by scanning the target area (i.e., the area covered by the drone) with a three-dimensional laser radar device on the drone is a laser point cloud, in which the data, i.e., the data points, are recorded as laser radar data. In the present invention, the laser point cloud data is mainly used to extract forest characteristic parameters such as single tree scale and sample plot scale, and the forest characteristic parameters are used as independent variables of subsequent models.

[0218] 3.1.2. Data processing.

[0219] LiDAR data processing includes preprocessing and postprocessing. LiDAR data preprocessing includes track solution and data solution. After preprocessing, a LiDAR point cloud file is generated. The data in the file is the original point cloud data. Postprocessing includes point cloud denoising, point cloud data filtering, digital elevation model generation, digital surface model generation, canopy height model generation, feature parameter extraction (if single tree scale modeling is required, single tree segmentation is required), etc.

[0220] (1) Point cloud denoising refers to the removal of outliers in the original point cloud data. Usually, noise points appear as isolated points or isolated point clusters, with a small number of laser points and data values ​​that are significantly higher or lower than the laser points of the ground objects.

[0221] (2) Point cloud data filtering refers to separating ground points and non-ground points from the original point cloud data. The accurate separation of ground points and non-ground points is the basis for generating digital elevation models (DEM), digital surface models (DSM), and canopy height models (CHM).

[0222] in,

[0223] DEM is a digital model used to express surface undulations. In the XYZ three-dimensional coordinate system, the elevation Z value is expressed as a continuous function of the plane coordinate variables X and Y. There are many methods for generating DEM, such as spline function interpolation, nearest point interpolation, Kriging interpolation, Inverse Distance Weighted (IDW), irregular triangulation and other methods.

[0224] DSM is a collective model that expresses the surface morphology of objects above the ground, and can also reflect the canopy conditions of trees at different heights.

[0225] CHM is a model that reflects the vertical height of the forest canopy surface from the ground, and can reflect the horizontal and vertical structure of the forest canopy in the spatial structure.

[0226] (3) Single tree segmentation uses the relative spacing between tree canopies and uses the maximum point as the seed point to perform single tree point cloud segmentation from top to bottom.

[0227] (4) Sample segmentation is to use the relative spacing between tree canopies, take the maximum point as the seed point, and perform sample point cloud segmentation from top to bottom.

[0228] (5) Feature parameter extraction, that is, extracting forest structure parameters from the preprocessed LiDAR data to obtain LiDAR feature parameters, such as tree height, crown width, LiDAR Biomass Index (LBI), and percentile tree height.

[0229] After the acquired data is processed, data entry is performed, the data is bound to the small class, and data quality inspection and adjustment are completed. The standardized data is output for subsequent data calculation.

[0230] 3.2 Data registration.

[0231] According to the distance error of the single tree position (the distance error between the measured single tree position obtained by the manual measurement method and the single tree segmentation position obtained by single tree segmentation), and the single tree height accuracy (the ratio of the error between the measured single tree height obtained by the manual measurement method and the single tree segmentation height obtained by single tree segmentation to the single tree segmentation height), the single tree sample matching of the sample plot is performed, that is, for the sample plot, the single tree monitored by the drone is matched with the actual single tree. At the same time, the biomass of the single tree, that is, the reference aboveground biomass of the single tree, and the single tree laser radar characteristic parameters extracted in 3.1.2, such as tree height, percentile tree height, crown width, laser radar biomass index and percentile tree height, are obtained from step 2.2.1. Specifically, the matching is to match the data obtained by manual measurement with the data obtained by scanning, that is, monitoring, of the three-dimensional laser radar device (i.e., laser radar) on the drone, find the same tree at the same position, and use the biomass of the single tree as the biomass prior value monitored by the drone to form the first data set for the subsequent laser radar biomass inversion model training. The first data set includes at least 20 data.

[0232] 3.3. Construct a LiDAR biomass inversion model. This includes the following steps:

[0233] 3.3.1. Model construction, including the following steps:

[0234] 3.3.1.1. Build a linear regression model.

[0235] The present invention adopts a linear regression model to explain the relationship between the laser radar characteristic parameters of the UAV and the biomass, that is, the linear regression model fits a straight line by establishing a linear equation to minimize the residual square sum between the predicted biomass and the actually observed biomass.

[0236] (1) Model construction:

[0237] Randomly select 3 / 4 samples in the prior biomass data set as the training set. In this step, the prior biomass data set is the first data set. Use linear regression method to construct the biomass model:

[0238]

[0239] Where:

[0240] is the measured whole plant (or aboveground) biomass of the kth tree, in kilograms (kg dm), where dm represents dry matter. In this step, the measured whole plant (or aboveground) biomass is the biomass measured by the manual measuring method;

[0241] is the monitored whole-plant (or above-ground) biomass of the k-th tree, in kilograms (kg dm). In this step, the monitored whole-plant (or above-ground) biomass is the whole-plant (or above-ground) biomass measured by the laser radar;

[0242] is the monitoring characteristic parameter of the kth tree, in meters or dimensionless. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter;

[0243] SSE is the residual sum of squares, SSE = ∑ k (b1 k -b2 k ) 2 , where b1 k is the first biomass of the kth tree, b2 k is the second biomass of the kth tree. In this step, b1 k for b2 k for

[0244] is the assumed linear function between the measured whole plant (or aboveground) biomass of the sample plot and the monitored characteristic parameters; by minimizing the SSE in the survey plot data, the parameters can be solved, and the linear function after solving the parameters is recorded as The fitted function is used to calculate the biomass of the whole plant (or aboveground).

[0245] l represents data directly obtained through monitoring technology or a value calculated based on the data. In this step, the monitoring technology is laser radar monitoring;

[0246] m represents the data directly obtained through field measurements or the value calculated based on the data;

[0247] k represents a single sample tree and is an integer greater than 0.

[0248] Based on formula (4), a model equation for the relationship between monitoring characteristic parameters and monitoring whole plant (or aboveground) biomass is established, where the univariate linear regression model is:

[0249]

[0250] Where:

[0251] a1 and b1 are model parameters, dimensionless;

[0252] The binary linear regression model is:

[0253]

[0254] Where:

[0255] a2, b2 and c2 are model parameters, dimensionless;

[0256] is the i1th monitoring characteristic parameter of the kth tree. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter, i1 is an integer greater than 0, and the unit is meter or dimensionless;

[0257] is the j1th monitoring characteristic parameter of the kth tree. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter, j1 is an integer greater than 0, and j1≠i1, and the unit is meter or dimensionless;

[0258] (2) Model verification:

[0259]

[0260] Where:

[0261] R 2 is the coefficient of determination of the model, dimensionless;

[0262] b1 k The average value of the unit is kilograms (kg dm),

[0263] In this step, b1 k for b2 k for Calculate the coefficient of determination of equations (5) and (6) and compare them, and select R 2 The equation with the best (i.e. maximum) value is taken as the final biomass prediction equation.

[0264] 3.3.1.2. Construct a machine learning algorithm based on the random forest regression model. Random forest can explain the effect of LiDAR feature parameters on biomass. Figure 2 The decision tree of the random forest regression model used in the present invention is shown in Figure 2 The machine learning algorithm in this step extracts n Bootstrap samples from the first data set to form a data set through Bootstrap resampling technology, where n is an integer greater than 1, and performs decision tree modeling on each Bootstrap sample to obtain n decision trees (i.e., decision tree-1, ..., decision tree-n), and then obtains the final prediction result through majority voting or averaging mechanism using the prediction results of the n decision trees (i.e., result-1, ..., result-n). The prediction result is the biomass.

[0265] (1) Model construction:

[0266] Randomly select 2 / 3 of the samples in the first data set as the training set, and use the randomForest function under the randomForest data package of R language to implement the lidar biomass model construction. The steps include:

[0267] A) Selection of laser radar characteristic parameters.

[0268] When using the random forest method to build a biomass model,

[0269] Aa), the random forest algorithm automatically evaluates the influence of each lidar feature parameter in the regression process based on two indicators: the model mean square error increment (IncMSE) and the influence on the purity of the model tree node (IncNodePurity).

[0270] Ab) Considering the values ​​of IncMSE and IncNodePurity comprehensively, the lidar characteristic parameters with higher values ​​are selected for model construction. Figure 3 and Figure 4 , parameters TreeHeight, H99, H98, H90, and H97 can be selected to construct the lidar biomass model, where the parameters mean: TreeHeight—tree height (maximum height); H99—99th percentile tree height; H98—98th percentile tree height; H90—90th percentile tree height; H97—97th percentile tree height.

[0271] B) Selection of the number of regression trees.

[0272] Ba), use the parameter ntree to define the number of decision trees, and the parameter mtry to define the number of randomly selected features;

[0273] Bb), determine the value of ntree by obtaining the relationship diagram of ntree and error of the model through different numbers of regression trees. According to the relationship diagram, at the point where the error tends to be stable (such as Figure 5 As shown in , ntree=200) is the number of regression trees;

[0274] Bc) mtry is not assigned a value, and the system default value of 1 / 3 of the number of input variables is used. If the number of characteristic parameters is less than 3, the value of mtry is 1.

[0275] (2) Model verification:

[0276] According to formula (7), if the model meets the minimum requirement, the model accuracy test index R 2 When the minimum value of is 0.7, the model passes the test and can be used for biomass calculation.

[0277] 3.4. Calculate target biomass and carbon storage.

[0278] The estimated value of the average biomass carbon storage per unit area of ​​each carbon layer obtained by airborne lidar monitoring is combined with the area of ​​the target area to calculate the biomass carbon storage within the target boundary. The steps include:

[0279] 3.4.1. Calculate the biomass and carbon storage of individual trees in the area covered by the drone (i.e. the target area).

[0280] First, based on the characteristic parameter data monitored by the airborne LiDAR and combined with the LiDAR biomass inversion model constructed in step 3.3, the whole plant (or aboveground) biomass of each single tree in the i2 target carbon layer is calculated to form a LiDAR single tree whole plant (or aboveground) biomass dataset of the j2 type of tree species in the i2 target carbon layer in the T2 year. Then, the single tree biomass is converted into single tree biomass carbon storage in combination with the carbon content of the biomass of each type of tree species. The LiDAR biomass inversion model uses the R algorithm in the machine learning algorithm based on the linear regression model and the random forest regression model. 2 The larger one.

[0281]

[0282] In the above formula:

[0283] is the LiDAR biomass carbon stock of the s2th tree of the j2th tree species in the i2th target carbon layer in the T2th year, in tons of carbon, T2 is an integer, and i2, j2 and s2 are integers greater than 0;

[0284] is the total plant (or aboveground) biomass of the s2th tree species in the i2th target carbon layer in the T2th year, in kilograms (kg dm), which is obtained by bringing the characteristic parameter data monitored by the airborne lidar into the model;

[0285] CF j2 is the carbon content of biomass of the j2th tree species, dimensionless.

[0286] 3.4.2. Calculate the estimated value of the average carbon storage per unit area of ​​the carbon layer and the target carbon storage.

[0287] The first step is to calculate the biomass carbon storage per unit area of ​​the j2 type of tree species:

[0288]

[0289] Where:

[0290] is the LiDAR biomass carbon stock per unit area of ​​the j2th tree species in the i2th target carbon layer in year T2, in tons of carbon;

[0291] Ai2,j2,t2 =The forest area of ​​the j2-type tree species in the i2-target carbon layer in year T2, in hectares (hm 2 ).

[0292] Step 2: Calculate the estimated value of the average forest biomass carbon storage per unit area in each carbon layer:

[0293]

[0294] Where:

[0295] is the average per unit area lidar biomass carbon stock of the i2th target carbon layer in year T2, in tons of carbon;

[0296] nn i2 is the number of tree species in the i2th target carbon layer, a positive integer, dimensionless.

[0297] Step 3: Combined with the area weight of the i2th target carbon layer within the target boundary, calculate the estimated value of the target overall average biomass carbon storage per unit area:

[0298]

[0299] Where:

[0300] is the average biomass carbon stock per unit area within the target boundary in year T2, in tons of carbon per hectare;

[0301] w i2 =A i2 / A, where, A i2 is the area of ​​the i2th target carbon layer, in hectares (hm 2 ).

[0302] Step 4: Calculate the total biomass carbon stock within the target boundary based on the total target area:

[0303]

[0304] Where:

[0305] It is the biomass carbon stock obtained through drone monitoring within the target boundary in year T2, in tons of carbon.

[0306] 4. Satellite data processing and calculation. This includes the following steps:

[0307] 4.1 Satellite data acquisition and processing (i.e. satellite characteristic parameter extraction), including the following steps:

[0308] 4.1.1. Data acquisition.

[0309] In order to obtain carbon storage information over a large area (all target areas), it is necessary to introduce ICESAT-2 satellite data, Global Ecosystem Dynamics Investigation LiDAR (GEDI) data, Landsat data, Sentinel-2 satellite data, etc. to build an integrated satellite-aircraft-ground carbon storage estimation framework.

[0310] After the acquired data is processed, data entry is performed, the data is bound to the small class, and data quality inspection and adjustment are completed. The standardized data is output for subsequent data calculation.

[0311] 4.1.2 Data processing. Including the following steps:

[0312] (1) ICESat-2 satellite data preprocessing operations:

[0313] ICESat-2 satellite provides 21 standard data products ATL00~ATL21 (no ATL05), which are divided into four levels: Level 0, Level 1, Level 2 and Level 3.

[0314] The present invention mainly uses the ATL08 land vegetation height product to extract data from the target area and obtain vegetation height information (maximum forest height, average forest height, median forest height, minimum forest height and forest height percentile). To ensure the quality of data products, data screening operations are performed to remove daytime data, data affected by clouds, and data with large uncertainty in forest height.

[0315] (2) GEDI data preprocessing operations:

[0316] GEDI is carried on the International Space Station, and its data products include four levels: L1, L2, L3, and L4.

[0317] The present invention mainly uses GEDI L2 products to extract data products of the target area and obtains the longitude, latitude, ground elevation and relative forest height, canopy coverage, leaf area index and other structural parameters of the light spot. In order to ensure data quality, invalid waveform data is removed through data screening operations.

[0318] (3) Landsat 8 satellite data preprocessing operations:

[0319] The present invention is based on the Google Earth Engine (GEE) platform to use Landsat 8 satellite data to synthesize cloud-free or low-cloud images of the target area, and extract various relevant characteristic parameters such as spectral bands and vegetation indexes for subsequent model construction and training.

[0320] (4) Sentinel 2 satellite data preprocessing operations:

[0321] Based on the GEE platform, the present invention uses Sentinel 2 satellite data to synthesize cloud-free or low-cloud images of the target area, and extracts various relevant characteristic parameters such as spectral bands and vegetation index for subsequent model construction and training.

[0322] After preprocessing the acquired data, the vegetation structure characteristic parameters at the sample scale, i.e., satellite remote sensing characteristic parameters, such as tree height, crown width, percentile tree height and other data (including maximum height, average height, maximum intensity, average intensity, etc. in point cloud data) are extracted (here, the satellite remote sensing data of the target area is divided into sample scale data, hereinafter referred to as sample data), and data entry is performed, and the data is bound to the small class to complete data quality inspection and adjustment. The standardized data is output for subsequent data calculation.

[0323] 4.2. Data registration.

[0324] According to the position distance error of the sample plot, that is, the distance error between the position of the drone laser radar sample plot and the position of the satellite laser radar sample plot, the first sample plot obtained by segmenting the satellite laser radar (i.e., satellite remote sensing) data and the second sample plot obtained by segmenting the drone laser radar data are matched. At the same time, the biomass of the individual trees obtained from step 3.4.1, the biomass of all the individual trees covered by the sample plot is accumulated as the sample plot biomass obtained by laser radar segmentation, that is, the reference aboveground biomass of the sample plot. The matching is to match the data obtained by drone laser radar monitoring with the satellite remote sensing data (which is the satellite remote sensing characteristic parameters after processing), find the same sample plot at the same position, and use the biomass of the same sample plot calculated based on the data obtained by drone laser radar monitoring (which is the laser radar characteristic parameters after processing) as the prior value of the biomass monitored by the satellite laser radar (i.e., satellite remote sensing biomass), forming a second data set for subsequent model training. The second data set includes at least 20 data.

[0325] 4.3. Constructing a satellite remote sensing biomass inversion model, including the following steps:

[0326] 4.3.1. Model construction, including the following steps:

[0327] 4.3.1.1. Construct a linear regression model.

[0328] Referring to step 3.3.1.1, the difference between this step and step 3.3.1.1 is that: the second data set is used as the prior biomass data set;

[0329] In formula (5) and (6), Replace with is the satellite remote sensing biomass of the kkth plot, in kilograms (kgd.m.), where kk represents the plot and is an integer greater than 0;

[0330] In formula (5) Replace with is the satellite remote sensing data of the kkth sample plot, in kilograms (kg dm);

[0331] In formula (6) and Replace with and is the i3th satellite remote sensing characteristic parameter of the kkth sample plot, i3 = 1, 2, 3, ..., in meters or dimensionless units, is the j3th satellite remote sensing characteristic parameter of the kkth sample plot, j3=1, 2, 3, ... (j3≠i3), in meters or dimensionless;

[0332] s represents data directly obtained through satellite remote sensing or a value calculated based on the data and is dimensionless.

[0333] 4.3.1.2. Machine learning algorithm based on random forest regression model.

[0334] See step 3.3.1.2. The difference between this step and step 3.3.1.2 is:

[0335] Taking the second data set as the prior biomass data set, the randomForest function under the randomForest data package of R language is used to realize the construction of the satellite remote sensing biomass model;

[0336] Replace the lidar feature parameters in step 3.3.1.2 with the satellite remote sensing feature parameters (UAV coverage area).

[0337] 4.4. Calculate target biomass and carbon storage.

[0338] Based on the characteristic parameter data of satellite remote sensing, combined with the satellite remote sensing biomass inversion model constructed in step 4.3, the biomass of the trees in the target area is calculated, and then the biomass carbon storage within the target boundary is calculated in combination with the biomass carbon content. Including:

[0339] 4.4.1. Calculate the tree biomass and carbon storage in the sample plots in the satellite coverage area (i.e. the target area).

[0340] First, based on the characteristic parameter data of satellite remote sensing, combined with the satellite remote sensing biomass inversion model constructed in step 4.3, the biomass of each sample plot in the i4th target carbon layer in the target area is calculated to form the satellite remote sensing sample plot biomass data set of the i4th target carbon layer in the T4th year, and then the sample plot biomass is converted into sample plot biomass carbon storage in combination with the biomass carbon content, T4 is an integer, and i4 and q4 are both integers greater than 0. Among them, the satellite remote sensing biomass inversion model uses the R in the machine learning algorithm based on the linear regression model and the random forest regression model 2 The larger one.

[0341]

[0342] In formula (13):

[0343] It is the satellite remote sensing biomass carbon stock of the q4th plot of the i4th target carbon layer in the T4th year, in tons of carbon.

[0344] The satellite remote sensing biomass of the q4th plot of the i4th target carbon layer in the T4th year, in kilograms (kgd.m.), is obtained by bringing the characteristic parameter data of satellite remote sensing monitoring into the satellite biomass model;

[0345] CF is the carbon content of biomass, dimensionless.

[0346] 4.4.2. Calculation of carbon reserves in carbon layer and target carbon reserves.

[0347] Step 1: Calculate the biomass carbon storage in each carbon layer:

[0348]

[0349] Where:

[0350] It is the satellite remote sensing biomass carbon stock of the i4th target carbon layer in the T4th year, in tons of carbon.

[0351] Step 2: Calculate the biomass carbon stock within the target boundary:

[0352]

[0353] It is the satellite remote sensing biomass carbon stock within the target boundary in year T4, in tons of carbon.

[0354] 5. Through comprehensive calculation, the final target carbon storage and carbon sink are obtained, thus completing the monitoring of carbon sink measurement. This step ensures the comprehensiveness and accuracy of the monitoring results and provides strong support for the decision-making of carbon sink measurement. This step includes:

[0355] The final target carbon storage and annual change (i.e. carbon sink) were obtained by weighted average comprehensive calculation:

[0356]

[0357] In the above formula:

[0358] C Biomass,T5 It is the final biomass carbon stock within the target boundary in year T5, in tons of carbon.

[0359]

[0360] In the above formula:

[0361] ΔC Biomass,T5 is the annual change in the target biomass carbon stock (carbon sink) in year T5 within the target boundary, in tons of carbon dioxide equivalent per year;

[0362] t1, t2 – the t1th year and t2th year after the start of monitoring, in years (a). t1, t2 and T5 are all integers, and t1≤T5≤t2.

[0363] The multi-source remote sensing data forestry carbon sink measurement and monitoring method provided by the present invention specifically includes the steps of technical specifications for data collection, geographic information service (GIS) and remote sensing data processing flow, mathematical model for biomass and carbon storage calculation, etc.

[0364] The core of the entire multi-source remote sensing data forestry carbon sink measurement and monitoring method is to effectively divide the carbon layer and design the sampling method, then collect data through different channels, and use a series of calculation methods to estimate biomass and carbon storage, and finally obtain the carbon storage of the entire target forest.

[0365] The present invention also provides a computer device, which includes a memory, a first processor, and a first computer program stored in the memory and executable on the first processor. When the first computer program is executed by the first processor, the multi-source remote sensing data forestry carbon sink measurement and monitoring method is implemented.

[0366] The present invention also provides a computer-readable storage medium, which is used to store a second computer program. The second computer program can be executed by at least one second processor so that the at least one second processor executes the multi-source remote sensing data forestry carbon sink measurement and monitoring method.

[0367] It is worth noting that, in terms of the calculation of biomass and carbon storage, the present invention first uses artificial gauge data to invert forest biomass and carbon storage through the biomass equation method. Secondly, drones are used to obtain sample data from a larger range for forest biomass and carbon storage inversion. During the inversion process, the gauge data is used to perform modeling and accuracy verification. Subsequently, the forest biomass and carbon storage are inverted by obtaining full coverage of sample data in the target area through satellites. During the inversion process, the gauge data and drone lidar data are used to perform modeling and accuracy verification. Finally, the final target carbon storage is obtained after cross-validation of multi-source remote sensing data.

[0368] The present invention effectively solves the problems of high measurement cost and low data quality in traditional manual measurement calculation methods through GIS technology, biomass calculation model, application of carbon coefficient and fusion processing of multi-source data, and can be widely used in carbon sink resource monitoring services and carbon sink target development monitoring technology services.

[0369] Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent substitutions for some of the technical features therein; and these modifications or substitutions do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for measuring and monitoring forestry carbon sinks using multi-source remote sensing data, characterized in that: include: Step 1: Design of monitoring plan; Step 2: Manual measurement data processing and calculation; Step 3: UAV data processing and calculation; Step 4: Satellite data processing and calculation; Step 5: Based on the steps 2, 3, and 4, the final target carbon storage and carbon sink are calculated, thereby completing the monitoring of carbon sink measurement. in, The step 1 comprises: Step 1.1, small class vector data management, including: obtaining vector data, i.e. small class vector data, entering, storing, managing, visualizing and applying shape data files composed of the vector data, including obtaining a vector data small class attribute table; Step 1.2, carbon layer factor screening and carbon layer division; Step 1.3: Calculate the number of carbon layer sample plots and systematically randomly distribute the sample plots, determine the sample plots, and generate the sample site locations.

2. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 1 is characterized in that: The step 1.2 comprises: The factor information in the vector data small class attribute table: tree species, age group, slope and the statistical information of the factors: data distribution map, factor maximum value, factor average value are used for small class classification. Specifically, the factors also include: county, township, village, forest farm, forest team, forest class, land type, forest species, terrain, slope direction, slope position, dominant tree species, age group, and canopy density. Different factors correspond to their own enumeration value groups one by one. The carbon layer division includes: selecting enumeration values ​​from different enumeration value groups in turn for combination, selecting one enumeration value from each enumeration value group for each combination, until all the enumeration values ​​are involved in the combination, and an enumeration value combination obtained by combination corresponds to a carbon layer. When combining, if N factors are involved, the N factors correspond to enumeration value group 1, ..., enumeration value group N one by one, and the enumeration value group 1, ..., enumeration value group N respectively include n1, ..., nn enumeration values, then this step divides n1×...×nn carbon layers, and outputs the carbon layer division result, i.e., the carbon layer result, which includes the carbon layer code, the carbon layer name, the carbon layer characteristics, i.e., the factor information and the statistical information of the factor, the number of small classes contained in the carbon layer, and the carbon layer area, wherein N, n1, ..., nn are integers greater than 0. In addition, the plots within the target boundary need to be divided into different levels according to different stratification factors.

3. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 2 is characterized in that: The step 1.3 comprises: According to the carbon layer division results, sampling formulas and statistical data in vector data, the number of sample plots and the area of ​​small classes in each carbon layer are calculated, and the sample plot layout is completed according to the method of random starting point and systematic layout. The sample plots are determined, and the sample plot distribution map, UAV flight area map and sample plot information table are output, including the sample plot number and the longitude and latitude coordinates of the sample plot location, that is, the sample plot location. According to the selected target tree species, for each target tree species, at least 60 individual trees are randomly selected as modeling samples at the target location. The sample tree diameter at breast height is not less than 2 cm, and they are divided according to the diameter class interval of 2 cm to 3 cm. There are no less than 2 to 3 sample trees in each diameter class to cover all diameter classes.

4. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 3 is characterized in that: The step 2 comprises: 2.2.

1. Calculation of biomass and carbon storage of sample plots: Based on the data obtained from the sample plot monitoring: tree species, breast diameter, tree height, crown width, coordinates, combined with the biomass equation method, the biomass of individual trees in the sample plot is calculated, and the biomass at the sample plot level is accumulated. Then, combined with the biomass carbon content, the biomass carbon storage at the sample plot level in each year and each target carbon layer is calculated; 2.2.

2. Calculation of estimated value of average carbon storage per unit area of ​​carbon layer and statistical estimation of target carbon storage, including the following steps: Step 1: Calculate the estimated value of the average biomass carbon storage per unit area for each carbon layer: Where: is the average biomass carbon stock per unit area of ​​the i-th target carbon layer in the T-th year, in tons of carbon per hectare, T is an integer, and i is an integer greater than 0; is the biomass carbon storage per unit area of ​​the i-th target carbon layer plot p in the T-th year, in tons of carbon per hectare, The biomass carbon stock per unit area of ​​the plot is obtained by dividing the biomass carbon stock at the plot level of the i-th target carbon layer in the T-th year obtained in step 2.2.1 by the area of ​​the plot p, in tons of carbon per hectare; In formula (1), n i is the number of monitoring plots of the i-th target carbon layer, which is a positive integer and dimensionless; m represents the data directly obtained through field measurements of the plots or the value calculated based on the data and dimensionless; Step 2: Calculate the estimated value of the average biomass carbon stock per unit area of ​​the target population: Where: In year T, the target is the average biomass carbon stock per unit area within the target area boundary, in tons of carbon per hectare; w i is the area weight of the i-th target carbon layer within the target boundary, dimensionless, w i =A i / A, where A is the total area of ​​each carbon layer within the target boundary, in hectares, A i is the area of ​​the i-th target carbon layer, in hectares; Step 3: Calculate the total biomass carbon stock within the target boundary in year T: Where: It is the biomass carbon stock within the target boundary in year T, in tons of carbon.

5. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 4 is characterized in that: The step 3 comprises: Data Registration: The single tree samples of the sample plot are matched according to the single tree position distance error and the single tree height accuracy. That is, for the sample plot, the single tree monitored by the drone is matched with the actual single tree. At the same time, the biomass of the single tree, that is, the reference aboveground biomass of the single tree, is obtained from the step 2.2.

1. Then, the single tree laser radar characteristic parameters are extracted from the laser radar data obtained by scanning the target area by the three-dimensional laser radar equipment on the drone: tree height, percentile tree height, crown width, laser radar biomass index and percentile tree height. Specifically, the matching is to match the data obtained by manual measurement with the data obtained by scanning or monitoring by the three-dimensional laser radar equipment on the drone, find the same tree in the same position, and use the biomass of the single tree as the prior value of the biomass monitored by the drone to form a first data set, which includes at least 20 data.

6. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 5 is characterized in that: The step 3 includes: 3.3, constructing a laser radar biomass inversion model: 3.3.1.

1. Construct a linear regression model: Randomly select 3 / 4 samples in the prior biomass data set as the training set. In this step, the prior biomass data set is the first data set, and the biomass model is constructed using the linear regression method: Where: is the measured whole plant or aboveground biomass of the kth tree, in kilograms. In this step, the measured whole plant or aboveground biomass is the biomass measured by the manual measuring method; is the monitored whole plant or aboveground biomass of the kth tree, in kilograms. In this step, the monitored whole plant or aboveground biomass is the laser radar whole plant or aboveground biomass; is the monitoring characteristic parameter of the kth tree, in meters or dimensionless. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter; SSE is the residual sum of squares, SSE = ∑ k (b1 k -b2 k ) 2 , where b1 k is the first biomass of the kth tree, b2 k is the second biomass of the kth tree. In this step, b1 k for b2 k for is the assumed linear function between the measured whole plant or aboveground biomass and the monitored characteristic parameters of the sample plot; the parameters are solved by minimizing the SSE in the survey plot data, and the linear function after solving the parameters is recorded as The fitted function is used to calculate the biomass of the whole plant (or aboveground). l represents data directly obtained through monitoring technology or a value calculated based on the data. In this step, the monitoring technology is laser radar monitoring; m represents the data directly obtained through field measurements or the value calculated based on the data; k represents a single sample tree and is an integer greater than 0. Based on formula (4), a model equation for the relationship between monitoring characteristic parameters and monitoring whole plant or aboveground biomass was established, where the univariate linear regression model is: Where: a1 and b1 are model parameters, dimensionless; The binary linear regression model is: Where: a2, b2 and c2 are model parameters, dimensionless; is the i1th monitoring characteristic parameter of the kth tree. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter, i1 is an integer greater than 0, and the unit is meter or dimensionless; is the j1th monitoring characteristic parameter of the kth tree. In this step, the monitoring characteristic parameter is the laser radar characteristic parameter, j1 is an integer greater than 0, and j1≠i1, and the unit is meter or dimensionless. Model verification: Where: R 2 is the coefficient of determination of the model, dimensionless; b1 k The average value, in kilograms, In this step, b1 k for b2 k for Calculate the coefficient of determination of equations (5) and (6) and compare them, and select R 2 The equation with the best value, i.e. the largest value, is taken as the final biomass prediction equation. 3.3.1.

2. Construct a machine learning algorithm based on a random forest regression model: The machine learning algorithm in this step uses the Bootstrap resampling technology to extract n Bootstrap samples from the first data set to form a data set, where n is an integer greater than 1, and performs decision tree modeling on each Bootstrap sample to obtain n decision trees: decision tree-1, ..., decision tree-n, and then the prediction results of the n decision trees, i.e., result-1, ..., result-n, are used to obtain the final prediction result through majority voting or an average mechanism. The prediction result is the biomass. When constructing the model, 2 / 3 of the samples in the first data set are randomly selected as a training set, and the randomForest function under the randomForest data package of the R language is used to implement the construction of the laser radar biomass model, including the steps of: A) Selection of laser radar characteristic parameters: When using the random forest method to build a biomass model, Aa), the random forest algorithm automatically evaluates the influence of each lidar feature parameter in the regression process based on the model mean square error increment IncMSE and the influence on the purity of the model tree node IncNodePurity; Ab), select parameters TreeHeight, H99, H98, H90, H97 to construct the lidar biomass model, where the parameters mean: TreeHeight—tree height, H99—99th percentile tree height, H98—98th percentile tree height, H90—90th percentile tree height, H97—97th percentile tree height, B) Selection of the number of regression trees: Ba), use the parameter ntree to define the number of decision trees, and the parameter mtry to define the number of randomly selected features; Bb), ntree = 200 as the number of regression trees; Bc) mtry is not assigned a value, and the system default value of 1 / 3 of the number of input variables is used. If the number of characteristic parameters is less than 3, the value of mtry is 1. Model verification: According to formula (7), if the model meets the minimum requirement, the model accuracy test index is the determination coefficient R 2 If the minimum value of is 0.7, the model passes the test and is used for biomass calculation.

7. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 6 is characterized in that: The step 3 also includes: 3.4, calculating the target biomass and carbon storage: the estimated value of the average biomass carbon storage per unit area of ​​each carbon layer obtained by airborne laser radar monitoring is combined with the area of ​​the target area to calculate the biomass carbon storage within the target boundary, including the steps of: 3.4.

1. Calculate the biomass and carbon storage of individual trees in the target area covered by the drone: In the above formula: is the LiDAR biomass carbon stock of the s2th tree of the j2th tree species in the i2th target carbon layer in the T2th year, in tons of carbon, T2 is an integer, and i2, j2 and s2 are integers greater than 0; is the total plant or aboveground biomass of the s2th tree species in the i2th target carbon layer in the t2th year, in kilograms, which is obtained by bringing the characteristic parameter data monitored by the airborne lidar into the lidar biomass inversion model. The lidar biomass inversion model uses the machine learning algorithm based on the linear regression model and the random forest regression model. 2 the larger one; CF j2 is the carbon content of the biomass of the j2-type tree species, dimensionless, 3.4.

2. Calculate the estimated value of the average carbon storage per unit area of ​​the carbon layer and the target carbon storage: The first step is to calculate the biomass carbon storage per unit area of ​​the j2 type of tree species: Where: is the LiDAR biomass carbon stock per unit area of ​​the j2th tree species in the i2th target carbon layer in year T2, in tons of carbon; A i2,j2,t2 is the forest area of ​​the j2-type tree species in the i2-target carbon layer in year T2, in hectares; Step 2: Calculate the estimated value of the average forest biomass carbon storage per unit area in each carbon layer: Where: is the average per unit area lidar biomass carbon stock of the i2th target carbon layer in year T2, in tons of carbon; nn i2 is the number of tree species in the i2th target carbon layer, a positive integer, dimensionless; Step 3: Combined with the area weight of the i2th target carbon layer within the target boundary, calculate the estimated value of the target overall average biomass carbon storage per unit area: Where: is the average biomass carbon stock per unit area within the target boundary in year T2, in tons of carbon per hectare; w i2 =A i2 / A, where, A i2 is the area of ​​the i2th target carbon layer, in hectares; Step 4: Calculate the total biomass carbon stock within the target boundary based on the total target area: Where: It is the biomass carbon stock obtained through drone monitoring within the target boundary in year T2, in tons of carbon.

8. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 7 is characterized in that: The step 4 comprises: Data Registration: According to the position distance error of the sample plot, that is, the distance error between the position of the UAV laser radar sample plot and the position of the satellite laser radar sample plot, the first sample plot measured by the satellite laser radar, that is, the satellite remote sensing data segmentation, and the second sample plot sample obtained by the UAV laser radar measurement data segmentation are matched. At the same time, the biomass of the individual trees obtained in step 3.4.1 is accumulated as the sample plot biomass obtained by laser radar segmentation, that is, the reference aboveground biomass of the sample plot. Among them, the matching is to match the data obtained by the drone lidar monitoring with the satellite remote sensing data, find the same sample at the same position, and use the biomass of the same sample calculated based on the data obtained by the drone lidar monitoring, recorded as the lidar biomass, as the biomass monitored by the satellite lidar, that is, the prior value of the satellite remote sensing biomass, to form a second data set. The second data set includes at least 20 data. The data obtained based on the drone lidar monitoring is processed into the lidar characteristic parameters, and the satellite remote sensing data is processed into the satellite remote sensing characteristic parameters.

9. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 8 is characterized in that: The step 4 also includes: 4.3.1.1, constructing a linear regression model, The difference between step 4.3.1.1 and step 3.3.1.1 is that: the second data set is used as the prior biomass data set; In formula (5) and (6), Replace with is the satellite remote sensing biomass of the kkth plot, in kilograms, where kk represents the plot and is an integer greater than 0; In formula (5) Replace with is the satellite remote sensing data of the kkth sample plot, in kilograms; In formula (6) and Replace with and is the i3th satellite remote sensing characteristic parameter of the kkth sample plot, i3 = 1, 2, 3, ..., in meters or dimensionless units, is the j3th satellite remote sensing characteristic parameter of the kkth sample plot, j3=1, 2, 3, ... (j3≠i3), in meters or dimensionless; s represents the data directly obtained through satellite remote sensing or the value calculated based on the data, dimensionless, The step 4 also includes: 4.3.1.2, a machine learning algorithm based on a random forest regression model, The difference between step 4.3.1.2 and step 3.3.1.2 is: The second data set was used as the prior biomass data set and the randomForest function was used in R language. The randomForest function under the package is used to build the satellite remote sensing biomass model; Replace the lidar characteristic parameters of the drone coverage area in step 3.3.1.2 with satellite remote sensing characteristic parameters.

10. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 9 is characterized in that: The step 4 further includes: 4.4, calculating the target biomass and carbon stock: based on the characteristic parameter data of satellite remote sensing, combined with the satellite remote sensing biomass inversion model constructed in step 4.3, calculating the biomass of the trees in the target area, and then calculating the biomass carbon stock within the target boundary in combination with the biomass carbon content, including: 4.4.

1. Calculate the tree biomass and carbon storage in the satellite coverage area, i.e. the target area, including: First, based on the characteristic parameter data of satellite remote sensing, combined with the satellite remote sensing biomass inversion model constructed in step 4.3, the biomass of each sample plot of the i4th target carbon layer in the target area is calculated to form a satellite remote sensing sample plot biomass data set of the i4th target carbon layer in the T4th year, and then the sample plot biomass is converted into sample plot biomass carbon storage in combination with the biomass carbon content, T4 is an integer, i4 and q4 are both integers greater than 0, wherein the satellite remote sensing biomass inversion model uses the R in the machine learning algorithm based on the linear regression model and the random forest regression model 2 The larger one, In formula (13): is the satellite remote sensing biomass carbon stock of the q4th plot of the i4th target carbon layer in the T4th year, in tons of carbon; is the satellite remote sensing biomass of the q4th plot of the i4th target carbon layer in the T4th year, in kilograms, which is obtained by bringing the characteristic parameter data of satellite remote sensing monitoring into the satellite remote sensing biomass inversion model; CF is the carbon content of biomass, dimensionless, 4.4.

2. Calculation of carbon reserves in carbon layer and target carbon reserves: Step 1: Calculate the biomass carbon storage in each carbon layer: Where: is the satellite remote sensing biomass carbon stock of the i4th target carbon layer in the T4th year, in tons of carbon, Step 2: Calculate the biomass carbon stock within the target boundary: It is the satellite remote sensing biomass carbon stock within the target boundary in year T4, in tons of carbon.

11. The multi-source remote sensing data forestry carbon sink measurement and monitoring method according to claim 10 is characterized in that: The step 5 comprises: The final target carbon storage and annual change, i.e. carbon sink, are calculated by weighted average: In the above formula: C Biomass,T5 is the final biomass carbon stock within the target boundary in year T5, in tons of carbon; In the above formula: ΔC Biomass,T5 is the annual change in the target biomass carbon stock in year T5 within the target boundary, in tons of carbon dioxide equivalent per year; t1, t2 – the t1th year and t2th year after the start of monitoring, expressed in years. t1, t2 and T5 are all integers, and t1≤T5≤t2.

12. Computer device, characterized in that It includes a memory, a first processor, and a first computer program stored in the memory and executable on the first processor, wherein the first computer program, when executed by the first processor, implements the multi-source remote sensing data forestry carbon sink measurement and monitoring method described in any one of claims 1 to 11.

13. A computer-readable storage medium, characterized in that The computer-readable storage medium is used to store a second computer program, and the second computer program can be executed by at least one second processor so that the at least one second processor executes the multi-source remote sensing data forestry carbon sink measurement and monitoring method described in any one of claims 1-11.

Images