CN105699624B - A kind of Soil Carbon Stock evaluation method based on soil genetic horizon thickness prediction - Google Patents

Abstract

The invention relates to a method for estimating soil organic carbon storage based on the prediction of the thickness of the soil occurrence layer, which provides a better technology for the prediction of soil continuity considering the horizontal dimension spatial distribution of soil attributes by encapsulating the unstructured information of the occurrence layer Among them, the technology of "merging occurrence layers, predicting and recalculating" is used to ensure that the characteristic information of soil occurrence layers is not missing, and at the same time, the limitations of traditional prediction methods in the description of soil continuity are corrected, and the "description specification, prediction Accurate" general-purpose soil organic carbon storage estimation technology has broad industrial application prospects in engineering surveys in agricultural applications, environmental protection, land resources and other related departments.

Description

A Soil Organic Carbon Storage Estimation Method Based on Prediction of Soil Genetic Layer Thickness

technical field

The invention relates to a method for estimating soil organic carbon storage based on the prediction of the thickness of soil origin layer, and belongs to the technical field of soil measurement.

Background technique

Soil carbon mainly includes soil organic carbon and soil inorganic carbon, and the soil inorganic carbon pool (carbonate carbon) is relatively stable. Soil organic carbon is mainly distributed in the soil at a depth of 1 meter, which has an important impact on the physical and chemical properties of the soil, and directly affects the soil quality. As one of the important indicators of soil fertility evaluation, the continuous reduction of soil organic carbon will directly lead to the impoverishment of agricultural soil. Soil organic carbon storage is not only one of the social issues that the ecological public interest survey focuses on, but also a global basic topic. Many application production and environmental protection departments at home and abroad have begun to pay extensive attention to the global warming problem faced by human survival and development. Relevant international regulations have also put forward different requirements for national-level carbon budget calculation technologies. The estimation of soil organic carbon stocks has become a major ecological and environmental production technology issue that affects the national economy and diplomacy.

Soil organic carbon density refers to the carbon storage at a fixed soil depth in a certain area, and the unit is kg/m2. The soil organic carbon storage in a certain area is the product of organic carbon density and area, and the unit is kg. Conventional carbon storage estimation techniques mainly include soil type method, vegetation type method, life zone type method, model method, correlation statistical method, and GIS (Geographic Information System) spatial prediction method. The spatial prediction method based on GIS is a spatial prediction mode widely used in modern digital soil mapping. Different from traditional soil survey and mapping techniques, this method adopts the quantitative soil-landscape model widely recognized by soil scientists at home and abroad, and can effectively integrate remote sensing image processing technology, digital terrain analysis technology, GIS spatial analysis technology and soil survey technology. Analysis of landscape information to predict the spatial distribution of soil organic carbon density. The main operation process of the GIS spatial prediction method is to estimate the organic carbon storage in the target area through field sampling in different geographical locations, laboratory analysis of soil physical and chemical properties, and establishment of a carbon storage estimation model based on environmental variables and known soil property data of sample points.

At present, commonly used techniques for estimating organic carbon stocks are mainly based on the following two production calculation models: calculate first, then predict (CTM), and predict first, then calculate (MTC). The "calculation" here is to calculate the soil organic carbon density ( Based on soil organic carbon content, soil bulk density, soil gravel content, soil thickness), "prediction" is spatial interpolation to predict the spatial distribution of the target soil properties. The application of the prior art mainly adopts one calculation mode, and there are few reference materials for comparing the application mode.

In recent years, with the development of computer, remote sensing, soil survey, and surveying and mapping technologies, especially the rapid development of high-resolution remote sensing and digital surveying and mapping technologies, the acquisition of regional-level environmental variable information including more heterogeneous details of the measured data has become a Possibly, the quantitative, objective, real-time, and accurate simulation of the spatial distribution of soil attributes in complex landscape areas from spatially discrete soil sample point data puts forward higher requirements for the specific production of soil information related departments, which has also developed into an international high-precision soil An important development direction of organic carbon storage estimation.

The occurrence layer is an important basis for identifying soil types, and has a series of quantitative descriptions on its properties. As the prediction basis for carbon storage estimation, the soil-landscape model assumes that the actual distribution of soil properties (including horizontal and vertical dimensions) is closely related to landscape properties. Especially in terms of the vertical dimension distribution of soil properties, the occurrence layers often have homogeneous soil property characteristics. However, the existing soil organic carbon storage estimation models are all oriented to the prediction of soil properties at a fixed depth. The joint agreement of the Global Soil Digital Mapping Project (GSM) stipulates that the soil layer thickness of the production of thematic soil maps adopts a fixed method of 0-5, 5-15, 15-30, 30-60, 60-100, and 100-200cm.

At present, there are certain limitations in the estimation of soil organic carbon storage in complex landscape areas and production techniques. Specifically, the following points can be summarized:

(1) The calculation model of fixed layer thickness ignores the theoretical model of soil genesis to some extent, and has certain limitations in production and processing. Due to topography, climate, and human disturbance, the complexity of soil vertical spatial variation far exceeds the simulation capabilities of existing computing technologies. If the attribute information of soil sampling points is missing, especially when the total number of sampling points is small and the spatial representation of the sampling points is poor, it is difficult to accurately estimate the regional SOC storage in practical applications. In fact, this is also one of the main reasons for the serious regional uncertainty in the existing digital soil mapping representation. The Nanjing Institute of Soil Science, Chinese Academy of Sciences and ISRIC - World Soil Information Reference Center (Netherlands) have pointed out that the application framework of modern digital soil mapping technology should not ignore the soil information characteristics of the occurrence layer.

(2) In the process of actual soil information production and application, it is difficult for the existing production technology to comprehensively consider the thickness of the soil genesis layer. Existing soil survey projects still adopt the traditional field survey mode, which not only cannot obtain real-time soil physical and chemical property data and descriptive information in large-scale areas, but is also often limited by the support of survey funds and the technical capabilities of personnel for survey operations. Soil system classification involves many soil diagnostic characteristics, and different diagnostic characteristics have quantitative and qualitative soil genetic model theory, which brings great difficulties to the actual production of technicians. Therefore, the existing basic soil organic carbon storage estimation projects often lack merging techniques that comprehensively consider the characteristics of soil genetic layers.

(3) Most of the existing technologies focus on a fixed computing mode, and lack of intensive comparative research. A large number of application cases at home and abroad have shown that different production models (CTM, MTC) often show very different regional accuracy problems in different model assumptions and landscape abrupt changes. Therefore, strengthening the comparison of different calculation modes has important guiding significance for specific production links. At the same time, the unstructured data structure of soil data and the quantitative soil characteristics of occurrence layers are important reasons for the imperfection of existing production technologies.

The insufficiency of the existing soil organic carbon storage estimation technology mentioned above has brought great difficulties in the production and processing of soil information products and engineering applications in different application departments, such as remote mountainous areas, areas with complex human influence, and highly variable forest areas. The engineering application of carbon stock estimation will bring potential wrong decision-making support, which will directly cause losses to the national economic planning.

Contents of the invention

The technical problem to be solved by the present invention is to provide a method for estimating soil organic carbon storage based on the prediction of soil genesis layer thickness, which adopts a new architecture design and can effectively improve the estimation accuracy and estimation work efficiency of soil organic carbon storage.

The present invention adopts the following technical solutions in order to solve the above-mentioned technical problems: the present invention designs a method for estimating soil organic carbon storage based on the prediction of the thickness of the soil formation layer, comprising the following steps:

Step 001. For the target soil area, each sample point location is set, and the number of statistical sample point locations is 1, and then enters step 002;

Step 002. In the target soil area, obtain the environmental information of the sample point for each sample point, and if the depth of the sample point vertically downward soil area is greater than or equal to the preset depth L, then obtain the sample point The profile area of the vertically downward preset depth L at the position; if the depth of the vertically downward soil area at the sample point is less than the preset depth L, then the profile area of the vertically downward soil at the sample point is obtained; and according to China According to the operation requirements of soil system classification, the soil occurrence layer is divided according to the section area, and each soil occurrence layer corresponding to the sample point position is obtained; then the thickness information and soil form information of each soil occurrence layer corresponding to the sample point position are respectively obtained; Then obtain the environmental information of each sample point position in the target soil area, and the thickness information and soil form information of each soil occurrence layer corresponding to each sample point position respectively, and then enter step 003;

Step 003. For each sample point position in the target soil area, respectively collect preset quality soil samples in each soil occurrence layer corresponding to the sample point position, and respectively measure and obtain the soil property information of each soil generation layer, Then obtain the soil attribute information of each soil occurrence layer corresponding to each sample point position in the target soil area, and then enter step 004; wherein, the soil attribute information includes soil organic carbon content, gravel content greater than a preset diameter, and soil bulk density;

Step 004. For each sample point position in the target soil area, according to the soil attribute information and thickness information of each soil occurrence layer corresponding to the sample point position, obtain the soil organic carbon density measured value SOCD _i of each sample point position, i= {1,...,I}, then go to step 005;

Step 005. For each sample point position in the target soil area, perform merge processing for each soil occurrence layer corresponding to the sample point position, obtain each merged layer corresponding to the sample point position, and then merge according to the sample point position For each soil occurrence layer corresponding to the corresponding layer, the weighted calculation is performed on the soil attribute information of each soil occurrence layer at the sampling point position to obtain the soil attribute information of each merged layer corresponding to the sampling point position, and the soil attribute information for each soil occurrence layer at the sampling point position The thickness information of the occurrence layer is summed and calculated to obtain the thickness of each merged layer corresponding to the sample point position; and then each merged layer corresponding to each sample point position in the target soil area, as well as the soil attribute information of each merged layer and Thickness information, then enter step 006;

Step 006. For each sample point position in the target soil area, determine whether the sum of the thicknesses of the merged layers corresponding to the sample point position is less than the preset depth L, if so, under each merged layer corresponding to the sample point position, Set the merged layer with the bedrock layer, so that the sum of the thicknesses of the merged layers corresponding to the sample point position is equal to the preset depth L, and enter step 007; otherwise, directly enter step 007;

Step 007. Obtain the types of merged layers corresponding to all sample point positions in the target soil area to form a set of merged layer types in the target soil area; according to the set of merged layer types in the target soil area, merge Layers are uniformly operated, so that the types of merged layers corresponding to each sample point position in the target soil area are the same as each other, and enter step 008;

Step 008. Using the linear congruence algorithm, divide the sample point locations in the target soil area into a set of predicted sample point locations and a set of verified sample point locations, and according to the measured value SOCD of soil organic carbon density at each sample point location in the target soil area _i , to obtain the measured value of soil organic carbon density at each verification sample point location in the verification sample point location set Constitute the set V of measured values of soil organic carbon density at each verification sample point location in the verification sample point location set; and enter step 009, where i ₂ ={1,...,I ₂ }, where I ₂ is the verification sample point location The number of verification sample point positions in the set, the number of predicted sample point positions in the set of predicted sample point positions is I ₁ , I ₁ >I ₂ ;

Step 009. According to the environmental information of each predicted sample point location in the predicted sample point location set, and the thickness information of each merged layer corresponding to each predicted sample point location, the random forest method is used to train and obtain the thickness information of each merged layer as Each prediction model of the target constitutes a first prediction model set;

At the same time, according to the environmental information of each predicted sample point location in the predicted sample point location set, and the soil attribute information of each merged layer corresponding to each predicted sample point location, the random forest method is used to train and obtain the soil attribute information of each merged layer respectively. Each prediction model for the target constitutes a second prediction model set; then enter step 010;

Step 010. For each verification sample point location in the verification sample point location set, according to the environmental information of the verification sample point location, through each prediction model in the first prediction model set, respectively obtain the corresponding each of the verification sample point locations. The first predicted thickness information of the merging layer; at the same time, for each verification sample point position in the verification sample point position set, according to the environmental information of the verification sample point position, through each prediction model in the second prediction model set, obtain the respective Verify the first predicted soil attribute information of each merged layer corresponding to the position of the sample point; and then obtain the first predicted thickness information and the first predicted soil attribute information of each merged layer corresponding to each verified sample point position in the verified sample point position set ; Then enter step 011; wherein, the first predicted soil attribute information includes the first predicted organic carbon content, the first predicted soil bulk density and the first predicted gravel content greater than the preset diameter;

Step 011. For each verification sample point position in the verification sample point position set, obtain the stretch coefficient Rc of the sum of the first predicted thickness information of each merging layer corresponding to the L relative to the verification sample point position; then according to the verification sample point position set The first predicted soil attribute information and the first predicted thickness information of each merged layer corresponding to each verification sample point position in , and obtain the first predicted value of soil organic carbon density of each verification sample point position in the verification sample point position set Constitute the first soil organic carbon density prediction value set PMTC-D of each verification sample point position set in the verification sample point position set; then enter step 012;

Step 012. According to the environmental information of each predicted sample point location in the predicted sample point location set, and the measured value of soil organic carbon density at each predicted sample point location, use the random forest method to train and obtain the target soil organic carbon density measured value The third prediction model; then, for each verification sample point location in the verification sample point location set, according to the environmental information of the verification sample point location, through the third prediction model, obtain the second soil organic carbon density prediction of the verification sample point location value, and then obtain the second predicted value of soil organic carbon density of each verification sample point location in the verification sample point location set, and form the second soil organic carbon density prediction value set PCTM of each verification sample point location in the verification sample point location set; then Go to step 013;

Step 013. For all sample point positions in the target soil area, according to the preset division rules, each merged layer corresponding to the sample point position is divided into each fitting layer based on L, and the fitting layer corresponding to each sample point position, and the number of fitting layers is the same, and the thickness information of each fitting layer is obtained; then, for each predicted sample point position in the predicted sample point position set, according to each merging layer corresponding to each fitting layer of the predicted sample point position , fitting the soil attribute information of each merged layer sampling occurrence layer at the predicted sample point position to obtain the soil attribute information of each fitted layer at the predicted sample point position; Corresponding to the soil attribute information of each fitting layer, and enter step 014;

Step 014. According to the environmental information of each predicted sample point location in the predicted sample point location set, and the soil attribute information of each fitting layer corresponding to each predicted sample point location, use the random forest method to train and obtain the soil of each fitting layer Each prediction model whose attribute information is the target constitutes a fourth prediction model set; then, for each verification sample point position in the verification sample point position set, according to the environmental information of the verification sample point position, through each of the fourth prediction model set The prediction model is used to obtain the third soil attribute prediction information of each fitting layer corresponding to the location of the verification sample point, and then obtain the third soil attribute prediction information of each fitting layer corresponding to each verification sample point location in the verification sample point location set information, and then enter step 015; wherein, the third predicted soil property information includes the third predicted information of organic carbon content, the third predicted information of soil bulk density and the third predicted content of gravel greater than the preset diameter;

Step 015. According to the third soil attribute prediction information and thickness information of each fitting layer corresponding to each verification sample point position in the verification sample point position set, obtain the third soil organic value of each verification sample point position in the verification sample point position set. Predicted value of carbon density Constitute the third soil organic carbon density prediction value set PMTC-F of each verification sample point position in the verification sample point position set; then enter step 016;

Step 016. Perform accuracy test according to the set of measured values of soil organic carbon density V, and obtain the first set of predicted values of soil organic carbon density PMTC-D, the second set of predicted values of soil organic carbon density PCTM, and the third set of predicted values of soil organic carbon density The optimal soil organic carbon density prediction value set in PMTC-F, and enter step 017;

Step 017. Discretize the spatial raster data of the target soil area, and use the sampling data of the occurrence layer corresponding to all sample points in the target soil area as the prediction data set. If the optimal soil organic carbon density prediction value set is the first soil For the set of predicted values of organic carbon density PMTC-D, the method from step 009 to step 011 is used to obtain the grid data of spatial distribution of soil organic carbon density in the target soil area; if the optimal set of predicted values of soil organic carbon density is the second soil organic carbon Density prediction value set PCTM, the method of step 012 is used to obtain the spatial distribution raster data of soil organic carbon density in the target soil area; if the optimal soil organic carbon density prediction value set is the third soil organic carbon density prediction value set PMTC-F , then use the method from step 013 to step 015 to obtain the grid data of the spatial distribution of soil organic carbon density in the target soil area; then go to step 018;

Step 018. According to the grid data of spatial distribution of soil organic carbon density in the target soil area, obtain the soil organic carbon storage of the target soil area.

As a preferred technical solution of the present invention, the step 001 specifically includes the following steps:

Step 00101. Obtain the soil type distribution map, soil use distribution map and soil geological distribution map of the target soil area, and spatially superimpose the soil type distribution map, soil use distribution map and soil geological distribution map to obtain the multi-layer superposition of the target soil area Figure, then enter step 00102;

Step 00102. Obtain the area of each patch area in the multi-layer overlay map of the target soil area, and obtain each patch area whose area ratio exceeds the threshold of the patch area; then, for each patch area, according to the field accessibility Detect requirements, set sample point position, the quantity of statistical sample point position is 1, then enter step 002.

As a preferred technical solution of the present invention: the environmental information includes vegetation coverage, rock outcrop area ratio, topography, important markers, size of coarse fragments on the surface, cracks on the surface, and salt spots on the surface.

As a preferred technical solution of the present invention: the soil morphology information represents soil dry and wet conditions, soil color, root system information, pore information, sample structure, speckle constituents, nodular nodules, cementation degree, and lime reaction information.

As a preferred technical solution of the present invention, the step 004 specifically includes the following operations:

For each sample point position in the target soil area, according to the soil attribute information and thickness information of each soil occurrence layer corresponding to the sample point position, according to the following formula:

Obtain the measured value of soil organic carbon density SOCD _i at each sampling point, where i={1,...,I}, n _i ={1,...N _i }, and n _i represents the value corresponding to the i-th sampling point The nth soil occurrence layer, N _i represents the total number of soil occurrence layers corresponding to the i-th sampling point position; Indicates the organic carbon content of the n-th soil layer corresponding to the i-th sample point position, Indicates the soil bulk density of the nth soil occurrence layer corresponding to the i-th sampling point position, Indicates the gravel content greater than the preset diameter in the n-th soil occurrence layer corresponding to the i-th sample point position, Indicates the thickness information of the nth soil occurrence layer corresponding to the ith sample point position, and then enters step 005.

As a preferred technical solution of the present invention, in the step 005, for each sample point position in the target soil area, each soil occurrence layer corresponding to the sample point position is merged and processed, and the sample point corresponding to the sample point position is obtained. The process of each merge layer, specifically includes the following steps:

Step 00501. For each sampling point position in the target soil area, according to the soil form information of each soil occurrence layer corresponding to the sampling point position, and the Chinese soil system classification, obtain the respective soil occurrence layers corresponding to the sampling point position Character expression symbols; and then obtain the characteristic expression symbols corresponding to each soil occurrence layer at each sample point position in the target soil area, and then enter step 00502;

Step 00502. According to the Chinese Soil System Classification, establish a merge table of occurrence layer characteristics, as shown in Table 1 below:

Table 1

Among them, the Arabic numerals are the same soil property expression symbols, and there are further divisions of each soil occurrence layer arranged in order from top to bottom from the soil profile;

Then, for each sample point position in the target soil area, and for each soil occurrence layer corresponding to the sampling point location, divide according to humus surface layer, sedimentary layer, and parent material layer, and merge the layers corresponding to each merged layer in the table based on the characteristics of the occurrence layer The characteristic expression symbol of the soil occurrence layer, from the section area of the sample point position from top to bottom, divide each soil generation layer corresponding to the sample point position in turn, and obtain each merged layer corresponding to the sample point position;

At the same time, in the process of performing the above merging operation for each sample point position in the target soil area, if the same soil occurrence layer characteristic expression symbol appears in the two merging levels corresponding to the sample point position, enter step 00503;

Step 00503. Determine whether there is a soil occurrence layer characteristic expression symbol located in the first merged level from top to bottom in the occurrence layer characteristic merger table among the respective soil occurrence layer characteristic expression symbols corresponding to the sample point position, and if so, put The soil occurrence layers respectively corresponding to the same soil occurrence layer characteristic expression symbols are merged, and then step 00502 is performed for merging processing, and each merged layer corresponding to the sample point position is obtained; otherwise, enter step 00504;

Step 00504. Among the same soil generatrix layer characteristic expression symbols, the soil generative layer characteristic expression symbols located in the lower merged layer are divided into the merged layer where the same soil generative layer characteristic expression symbols are located, and the sample point position The corresponding expression symbols for the characteristics of each soil occurrence layer move up one merged level according to the order of the merged levels in the merged table of the occurrence layer characteristics, and obtain the respective merged layers corresponding to the sample point positions.

As a preferred technical solution of the present invention, in the step 007, according to the target soil area merged layer type set, a unified operation is performed on the merged layer corresponding to each sample point position in the target soil area, so that each in the target soil area The types of merging layers corresponding to the sample point positions are the same, including the following operations:

For each sampling point position in the target soil area, judge whether the merged layer corresponding to the sample point position is equal to the type set of the merged layer in the target soil area, and if yes, do nothing; otherwise, set its relative target soil area merged The merged layer that is missing in the layer type set, and the thickness information of the merged layer is set to 0.00001cm, and the soil attribute information of the merged layer is set to 0.00001; and then the merged layers corresponding to each sample point position in the target soil area are equal to The target soil area is merged with a collection of layer types.

As a preferred technical solution of the present invention, the step 011 specifically includes the following operations:

For each verification sample point position in the verification sample point position set, obtain the stretch coefficient Rc of the sum of the first predicted thickness information of each merged layer corresponding to the L relative to the verification sample point position; then according to each verification sample point position set in the verification sample point position The first predicted soil attribute information and the first predicted thickness information of each merged layer corresponding to the sample point positions are according to the following formula:

Obtain the first predicted value of soil organic carbon density for each verification sample point location in the verification sample point location set Constitute the first set of predicted values of soil organic carbon density PMTC-D for each verification sample point location in the verification sample point location set; where, Indicates the total number of merging layers corresponding to the ith ₂ verification sample point positions in the verification sample point position set, Indicates the n-th merge layer corresponding to the i _- th verification sample point position in the verification sample point location set, Indicates the prediction information of the first organic carbon content of the n-th merged layer corresponding to the i _- th verification sample point position in the verification sample point position set, Indicates the first soil bulk density prediction information of the n-th merged layer corresponding to the i _- th verification sample point position in the verification sample point location set, Indicates the predicted content of the first gravel larger than the preset diameter in the n-th merged layer corresponding to the i- _2th verification sample point position in the verification sample point location set, Indicates the first predicted thickness information of the nth merged layer corresponding to the _i2th verification sample point position in the verification sample point position set; then enter step 012.

As a preferred technical solution of the present invention, in the step 013, the soil attribute information of each fitting layer corresponding to each predicted sample point position in the predicted sample point position set is obtained, specifically including the following operations:

For each predicted sample point position in the predicted sample point position set, according to the equal-area Spline function, the following formula is adopted:

Obtain the soil attribute information of each fitting layer at the predicted sample point location; and then obtain the soil attribute information of each fitting layer corresponding to each predicted sample point location in the predicted sample point location set, wherein, i ₁ ={1,..., I ₁ }, I ₁ represents the number of predicted sample point locations in the predicted sample point location set, Indicates the number of layers of the fitting layer corresponding to the i ₁ predicted sample point position in the predicted sample point position set, Indicates the soil attribute information of the k-th layer fitting layer corresponding to the i- _1th predicted sample point position in the set of predicted sample point positions, express function in layer with mean of layer fitting results, Indicates the laboratory measurement and analysis error of the k-th layer fitting layer corresponding to the i- _1th predicted sample point position in the set of predicted sample point positions, The function is an equal-area Spline function.

As a preferred technical solution of the present invention, the step 015 specifically includes the following operations:

According to the third soil attribute prediction information and thickness information of each fitting layer corresponding to each verification sample point position in the verification sample point position set, according to the following formula:

Obtain the third predicted value of soil organic carbon density for each verification sample point location in the verification sample point location set Constitute the third predicted value set PMTC-F of soil organic carbon density in each verification sample point location set in the verification sample point location set; where, Indicates the total number of fitting layers corresponding to the i _2th verification sample point position in the verification sample point position set, Indicates the kth fitting layer corresponding to the i ₂ verification sample point position in the verification sample point position set, Indicates the prediction information of the third organic carbon content of the k-th fitting layer corresponding to the i- _2th verification sample point position in the set of verification sample point positions, Indicates the third soil bulk density prediction information of the k-th fitting layer corresponding to the i- _2th verification sample point position in the verification sample point location set, Indicates the predicted content of the third gravel greater than the preset diameter in the k-th fitting layer corresponding to the i- ₂ verification sample point position in the verification sample point location set, Indicates the thickness information of the k-th fitting layer corresponding to the i- _2th verification sample point position in the verification sample point position set; then enter step 016.

Compared with the prior art, a method for estimating soil organic carbon storage based on the prediction of soil genesis layer thickness according to the present invention has the following technical effects:

(1) A method for estimating soil organic carbon storage based on the prediction of the thickness of the soil genesis layer is designed in the present invention. By encapsulating the unstructured information of the genesis layer, it provides a basis for the prediction of soil continuity considering the horizontal dimension spatial distribution of soil attributes. A good technical idea; Among them, the technology of "merging occurrence layers, forecasting and recalculation" is adopted to ensure that the characteristic information of soil occurrence layers is not missing, and at the same time, the limitations of traditional prediction methods in the description of soil continuity are corrected, and " The general-purpose soil organic carbon storage estimation technology with standardized description and accurate prediction has broad industrial application prospects in engineering surveys in agricultural applications, environmental protection, land resources and other related departments;

(2) In the method for estimating soil organic carbon storage based on the prediction of the thickness of the soil genetic layer designed by the present invention, the proposed thickness prediction of the genetic layer and the calculation of the stretching coefficient have certain universality, and its technical scheme is not only oriented to the soil The prediction of organic carbon density can also be combined with the prediction of occurrence layer type to constitute the technical process of soil type prediction. Therefore, the prediction technology proposed in the present invention also has good stability, and the dynamically constructed generation layer symbol merging mechanism is expected to ensure While improving the prediction accuracy, it also reduces the error in the prediction process of the soil occurrence layer thickness;

(3) In the method for estimating soil organic carbon storage based on the prediction of soil genetic layer thickness, the proposed genetic layer merging technology can simplify the structural application of soil information to the greatest extent, and is expected to be used in other engineering applications. For example, it provides technical guidance for three-dimensional soil mapping, global digital soil mapping, and regional carbon storage spatio-temporal evolution analysis, and the enhanced comparison model used takes into account the advantages of different calculation models to ensure that the optimal calculation model can be used to analyze soil organic carbon density. Carry out spatial prediction, and then be able to more quantitatively and objectively evaluate the organic carbon storage of the target area.

Description of drawings

Fig. 1 is the schematic flow chart of the method for estimating soil organic carbon storage based on the thickness prediction of the soil genesis layer designed by the present invention;

Fig. 2a is a schematic diagram of a soil profile area collected at a depth of 1m when the soil depth is greater than or equal to a preset depth of 1m in an embodiment of the present invention;

Fig. 2b is a schematic diagram of the vertically downward soil profile area collected when the soil depth is less than the preset depth 1m in the embodiment of the present invention;

Fig. 3 is a schematic diagram of the spatial variation of different soil layer thicknesses in the embodiment of the present invention;

Fig. 4 is a schematic diagram of obtaining a fixed layer thickness merged layer for 1m depth profile area division in an embodiment of the present invention;

Figures 5a to 5i are spatial distribution diagrams of the predicted thickness of each merged layer in the embodiment of the present invention;

Fig. 6 is a schematic diagram of the distribution of the tensile coefficient of the 1m depth profile area relative to the predicted thickness in the embodiment of the present invention;

7a to 7c are schematic diagrams of comparative analysis of prediction results of different methods in the embodiment of the present invention.

detailed description

The specific implementation manners of the present invention will be further described in detail below in conjunction with the accompanying drawings.

Aiming at the above-mentioned technical defects, the present invention designs a new set of engineering technical schemes for estimating soil organic carbon storage based on the prediction of soil genesis layer thickness, guided by the expert knowledge of soil evolution. This technology systematically gives the technical process of the estimation and prediction of organic carbon storage. At the same time, considering the high spatial variation characteristics of soil organic carbon density, the present invention innovatively proposes to compare and use different carbon storage estimation calculation models. Enhanced comparison strategies are expected to provide optimal computational models that are neglected by traditional GIS-based estimation techniques.

The basic idea of the invention is to systematically record the occurrence layer characteristics and description information of different soil bodies during the soil survey and sampling process. After laboratory soil property assay analysis, use the structured soil property data structure to input unstructured soil information. Taking the expression symbol of occurrence layer characteristics as the operation unit, the standardized soil profile data with similar occurrence layers and equal number of occurrence layers were constructed. The machine learning algorithm is used to predict the occurrence layer thickness and soil properties of the target point, and the calculation of soil organic carbon density is updated based on the stretch coefficient of the predicted occurrence layer thickness. By comparing and analyzing other calculation models, the optimal calculation model is selected, and the organic carbon storage in the target area is finally estimated.

As shown in Fig. 1, a kind of method for estimating soil organic carbon storage based on the prediction of soil genesis layer thickness designed by the present invention, in the actual application process, specifically includes the following steps:

Step 001. For the target soil area, set each sample point location, and count the number of sample point locations as I, and then enter step 002.

Wherein, step 001 specifically includes the following steps:

Step 00101. Obtain the soil type distribution map, soil use distribution map and soil geological distribution map of the target soil area, and spatially superimpose the soil type distribution map, soil use distribution map and soil geological distribution map to obtain the multi-layer superposition of the target soil area Figure, and then go to step 00102.

Step 00102. Obtain the area of each patch area in the multi-layer overlay map of the target soil area, and obtain each patch area whose area ratio exceeds the threshold of the patch area; then, for each patch area, according to the field accessibility Detect requirements, set sample point position, count the quantity of sample point position as 1, then enter step 002.

Step 002. In the target soil area, obtain the environmental information of each sample point respectively for each sample point location. The environmental information includes vegetation coverage, rock outcrop area ratio, topography, important markers, size of coarse fragments on the surface, and cracks on the surface situation, surface salt patch information; at the same time, as shown in Figure 2a, if the depth of the vertically downward soil area of the sampling point is greater than or equal to the preset depth L, then the vertically downward preset depth L of the sampling point is obtained Profile area; as shown in Figure 2b, if the depth of the vertically downward soil area at the sample point is less than the preset depth L, then the profile area of the vertically downward soil at the sample point is obtained; and according to the Chinese soil system classification According to the operation requirements, the soil occurrence layer is divided according to the section area, and each soil occurrence layer corresponding to the sampling point position is obtained; then the thickness information and soil form information of each soil occurrence layer corresponding to the sampling point position are obtained respectively, and the soil form information represents Soil wet and dry conditions, soil color, root system information, pore information, sample structure, speckle constituents, nodular nodules, cementation degree, and lime reaction information; and then obtain the environmental information of each sample point in the target soil area, and each The thickness information and soil form information of each soil occurrence layer corresponding to the sample point positions respectively, and then enter step 003 .

Among them, the Chinese Soil System Classification Table is shown in Table 2 below:

Occurrence layer symbol Occurrence layer designation information description Occurrence layer symbol Occurrence layer designation information description A humus surface l Reticulate B deposition layer m strong cement C parent layer no Sodium accumulation R bedrock o root knot a highly decomposed organic matter p farming impact b buried layer q secondary silicon accumulation c crust r Redox d freeze-thaw sheet the s Iron manganese accumulation e semi-decomposed organic matter t Clay accumulation f Permafrost u man-made accumulation g Lading spots v Transgender features h humus accumulation w Colored, structured layer formed by in situ weathering i Low and undecomposed organic matter x Solid and hard cement, not forming a pan j jarosite the y Gypsum accumulation k Carbonate accumulation z Soluble salt accumulation

Table 2

Step 003. For each sample point position in the target soil area, respectively collect preset quality soil samples in each soil occurrence layer corresponding to the sample point position, and respectively measure and obtain the soil property information of each soil generation layer, Then obtain the soil attribute information of each soil occurrence layer corresponding to each sample point position in the target soil area, and then enter step 004; wherein, the soil attribute information includes soil organic carbon content (unit: g/kg), greater than the preset diameter of 0.2 Gravel content in mm (volume percentage, %) and soil bulk density (unit: g/cm ³ ).

Step 004. For each sample point position in the target soil area, according to the soil attribute information and thickness information of each soil occurrence layer corresponding to the sample point position, obtain the soil organic carbon density measured value SOCD _i of each sample point position, i= {1,...,I}, then go to step 005.

The step 004 specifically includes the following operations:

For each sample point position in the target soil area, according to the soil attribute information and thickness information of each soil occurrence layer corresponding to the sample point position, according to the following formula:

Obtain the measured value of soil organic carbon density SOCD _i at each sampling point, where i={1,...,I}, n _i ={1,...N _i }, and n _i represents the value corresponding to the i-th sampling point The nth soil occurrence layer, N _i represents the total number of soil occurrence layers corresponding to the i-th sampling point position; Indicates the organic carbon content of the n-th soil layer corresponding to the i-th sample point position, Indicates the soil bulk density of the nth soil occurrence layer corresponding to the i-th sampling point position, Indicates the gravel content greater than the preset diameter in the n-th soil occurrence layer corresponding to the i-th sample point position, Indicates the thickness information of the nth soil occurrence layer corresponding to the ith sample point position, and then enters step 005.

Step 005. For each sample point position in the target soil area, perform merge processing for each soil occurrence layer corresponding to the sample point position, and obtain each merged layer corresponding to the sample point position; and then merge according to the sample point position For each soil occurrence layer corresponding to the corresponding layer, the weighted calculation is performed on the soil attribute information of each soil occurrence layer at the sampling point position to obtain the soil attribute information of each merged layer corresponding to the sampling point position, and the soil attribute information for each soil occurrence layer at the sampling point position The thickness information of the occurrence layer is summed and calculated to obtain the thickness of each merged layer corresponding to the sample point position; and then each merged layer corresponding to each sample point position in the target soil area, as well as the soil attribute information of each merged layer and Thickness information, then enter step 006.

In step 005, for each sample point position in the target soil area, perform merge processing for each soil occurrence layer corresponding to the sample point position, and obtain each merged layer corresponding to the sample point position, specifically including the following steps:

Step 00501. For each sampling point position in the target soil area, according to the soil form information of each soil occurrence layer corresponding to the sampling point position, and the Chinese soil system classification, obtain the respective soil occurrence layers corresponding to the sampling point position Character expression symbols; and then obtain the characteristic expression symbols corresponding to each soil occurrence layer at each sample point position in the target soil area, and then go to step 00502.

Step 00502. According to the Chinese Soil System Classification, establish a merge table of occurrence layer characteristics, as shown in Table 1 below:

Table 1

Among them, the Arabic numerals represent the same soil characteristic expression symbols, and there are further divisions of soil occurrence layers arranged in order from top to bottom of the soil profile; Note: (1) The main occurrence layers or characteristic occurrence layers can be classified according to their occurrence degree The differences are further subdivided into several sublayers. They are all represented by Arabic numerals and capital letters in parallel, such as C1, C2, Bt1, Bt2, Bt3; (2) Representation of heterogeneous parent material soil layer: represented by Arabic numerals placed in front of the occurrence layer symbol, such as 2C.

Then, for each sample point position in the target soil area, and for each soil occurrence layer corresponding to the sampling point location, divide according to humus surface layer, sedimentary layer, and parent material layer, and merge the layers corresponding to each merged layer in the table based on the characteristics of the occurrence layer The characteristic expression symbol of the soil occurrence layer, from the section area of the sample point position from top to bottom, divide each soil generation layer corresponding to the sample point position in turn, and obtain the merged layers corresponding to the sample point position.

At the same time, in the process of performing the above merging operation on each sample point position in the target soil area, if the same soil occurrence layer characteristic expression symbol appears in the two merging levels corresponding to the sample point position, go to step 00503.

Step 00503. Determine whether there is a soil occurrence layer characteristic expression symbol located in the first merged level from top to bottom in the occurrence layer characteristic merger table among the respective soil occurrence layer characteristic expression symbols corresponding to the sample point position, and if so, put The soil occurrence layers corresponding to the same soil occurrence layer characteristic expression symbols are merged, and then step 00502 is performed for merging processing to obtain each merged layer corresponding to the sampling point position; otherwise, go to step 00504.

Step 00504. Among the same soil generatrix layer characteristic expression symbols, the soil generative layer characteristic expression symbols located in the lower merged layer are divided into the merged layer where the same soil generative layer characteristic expression symbols are located, and the sample point position The corresponding expression symbols for the characteristics of each soil occurrence layer move up one merged level according to the order of the merged levels in the merged table of the occurrence layer characteristics, and obtain the respective merged layers corresponding to the sample point positions.

Step 006. As shown in Figure 3, for each sample point position in the target soil area, determine whether the sum of the thicknesses of each merged layer corresponding to the sample point position is less than the preset depth L, if so, at the corresponding sample point position Under each merged layer, the bedrock layer is used to set the merged layer so that the sum of the thicknesses of each merged layer corresponding to the sample point position is equal to the preset depth L, and proceed to step 007; otherwise, directly proceed to step 007.

Step 007. Obtain the types of merged layers corresponding to all sample point positions in the target soil area to form a set of merged layer types in the target soil area; according to the set of merged layer types in the target soil area, merge Perform unified operations on the layers so that the types of merged layers corresponding to each sample point position in the target soil area are the same as each other, and enter step 008.

In step 007, according to the type set of the merged layer in the target soil area, a unified operation is performed on the merged layer corresponding to each sample point position in the target soil area, so that the types of the merged layer corresponding to each sample point position in the target soil area are the same, Specifically include the following operations:

For each sampling point position in the target soil area, judge whether the merged layer corresponding to the sample point position is equal to the type set of the merged layer in the target soil area, and if yes, do nothing; otherwise, set its relative target soil area merged The merged layer that is missing in the layer type set, and the thickness information of the merged layer is set to 0.00001cm, and the soil attribute information of the merged layer is set to 0.00001; and then the merged layers corresponding to each sample point position in the target soil area are equal to The target soil area is merged with a collection of layer types.

Step 008. Using the linear congruence algorithm, divide the sample point locations in the target soil area into a set of predicted sample point locations and a set of verified sample point locations, and according to the measured value SOCD of soil organic carbon density at each sample point location in the target soil area _i , to obtain the measured value of soil organic carbon density at each verification sample point location in the verification sample point location set Constitute the set V of measured values of soil organic carbon density at each verification sample point location in the verification sample point location set; and enter step 009, where i ₂ ={1,...,I ₂ }, where I ₂ is the verification sample point location The number of verification sample point positions in the set is equal to 25%I, and the number of predicted sample point positions in the set of predicted sample point positions is I ₁ , which is equal to 75%I.

Step 009. According to the environmental information of each predicted sample point location in the predicted sample point location set, and the thickness information of each merged layer corresponding to each predicted sample point location, the random forest method is used to train and obtain the thickness information of each merged layer as Each prediction model of the target constitutes a first prediction model set.

At the same time, according to the environmental information of each predicted sample point location in the predicted sample point location set, and the soil attribute information of each merged layer corresponding to each predicted sample point location, the random forest method is used to train and obtain the soil attribute information of each merged layer respectively. For each prediction model of the target, form a second prediction model set; then go to step 010.

Step 010. Using the MTC-D method, for each verification sample point location in the verification sample point location set, according to the environmental information of the verification sample point location, through each prediction model in the first prediction model set, obtain the verification sample respectively. The first predicted thickness information of each merged layer corresponding to the point positions, as shown in Figure 5a to Figure 5i; at the same time, for each verification sample point position in the verification sample point position set, according to the environmental information of the verification sample point position, Through each prediction model in the second prediction model set, the first predicted soil attribute information of each merged layer corresponding to the verification sample point position is obtained respectively; The first predicted thickness information and the first predicted soil attribute information of the merged layer; then enter step 011; wherein, the first predicted soil attribute information includes the first predicted organic carbon content information, the first soil bulk density predicted information and the Predicted content of the first gravel.

Step 011. For each verification sample point position in the verification sample point position set, obtain the stretch coefficient Rc of the sum of the first predicted thickness information of each merging layer corresponding to the L relative to the verification sample point position, as shown in Figure 6; then According to the first predicted soil attribute information and the first predicted thickness information of each merged layer corresponding to each verification sample point position in the verification sample point position set, obtain the first soil organic carbon of each verification sample point position in the verification sample point position set Density Prediction Constitute the first set of predicted values of soil organic carbon density PMTC-D for each verification sample point location set in the verification sample point location set; then go to step 012.

Wherein, step 011 specifically includes the following operations:

For each verification sample point position in the verification sample point position set, obtain the stretch coefficient Rc of the sum of the first predicted thickness information of each merged layer corresponding to the L relative to the verification sample point position; then according to each verification sample point position set in the verification sample point position The first predicted soil attribute information and the first predicted thickness information of each merged layer corresponding to the sample point positions are according to the following formula:

Obtain the first predicted value of soil organic carbon density for each verification sample point location in the verification sample point location set Constitute the first set of predicted values of soil organic carbon density PMTC-D for each verification sample point location in the verification sample point location set; where, Indicates the total number of merging layers corresponding to the ith ₂ verification sample point positions in the verification sample point position set, Indicates the n-th merge layer corresponding to the i _- th verification sample point position in the verification sample point location set, Indicates the prediction information of the first organic carbon content of the n-th merged layer corresponding to the i _- th verification sample point position in the verification sample point position set, Indicates the first soil bulk density prediction information of the n-th merged layer corresponding to the i _- th verification sample point position in the verification sample point location set, Indicates the predicted content of the first gravel larger than the preset diameter in the n-th merged layer corresponding to the i- _2th verification sample point position in the verification sample point location set, Indicates the first predicted thickness information of the nth merged layer corresponding to the _i2th verification sample point position in the verification sample point position set; then enter step 012.

Step 012. Using the CTM method, according to the environmental information of each predicted sample point location in the predicted sample point location set, and the measured value of soil organic carbon density at each predicted sample point location, the random forest method is used to train and obtain the measured soil organic carbon density The value is the third prediction model of the target; then, for each verification sample point location in the verification sample point location set, according to the environmental information of the verification sample point location, through the third prediction model, the second soil of the verification sample point location is obtained The predicted value of organic carbon density, and then obtain the second predicted value of soil organic carbon density of each verification sample point in the verification sample point set, which constitutes the second predicted value of soil organic carbon density of each verification sample point in the verification sample point set Assemble PCTM; then go to step 013.

Step 013. As shown in Figure 4, for all sample point positions in the target soil area, according to the preset division rules, each merged layer corresponding to the sample point position is divided into each fitting layer based on L, and each sample point position The corresponding fitting layer and the number of fitting layers are the same, and the thickness information of each fitting layer is obtained. In this embodiment, for the section area with a depth of 1m, it is specifically divided into 0-5cm, 5-15cm, 15cm Each fitting layer with a fixed depth of -30cm, 30-60cm and 60-100cm is divided into 5 fitting layers; For each merged layer corresponding to the fitting layer, the soil attribute information of each merged layer sampling occurrence layer at the predicted sample point position is fitted to obtain the soil attribute information of each fitted layer at the predicted sample point position; and then the predicted sample point is obtained The soil attribute information of each fitting layer corresponding to each predicted sample point position in the position set, and enter step 014.

Wherein, in the step 013, the soil attribute information of each fitting layer corresponding to each predicted sample point position in the predicted sample point position set is obtained, specifically including the following operations:

For each predicted sample point position in the predicted sample point position set, according to the equal-area Spline function, the following formula is adopted:

Obtain the soil attribute information of each fitting layer at the predicted sample point location; and then obtain the soil attribute information of each fitting layer corresponding to each predicted sample point location in the predicted sample point location set, wherein, i ₁ ={1,..., I ₁ }, I ₁ represents the number of predicted sample point locations in the predicted sample point location set, Indicates the number of layers of the fitting layer corresponding to the i ₁ predicted sample point position in the predicted sample point position set, Indicates the soil attribute information of the k-th layer fitting layer corresponding to the i- _1th predicted sample point position in the set of predicted sample point positions, express function in layer with mean of layer fitting results, Indicates the laboratory measurement and analysis error of the k-th layer fitting layer corresponding to the i- _1th predicted sample point position in the set of predicted sample point positions, The function is an equal-area Spline function.

Step 014. Using the MTC-F method, according to the environmental information of each predicted sample point location in the predicted sample point location set, and the soil attribute information of each fitting layer corresponding to each predicted sample point location, the random forest method is used to train and obtain Each prediction model with the soil property information of each fitting layer as the target constitutes the fourth prediction model set; For each prediction model in the prediction model set, the third soil attribute prediction information corresponding to each fitting layer corresponding to the verification sample point position is obtained, and then each corresponding fitting layer corresponding to each verification sample point position in the verification sample point position set is obtained The third predicted soil attribute information, and then enter step 015; wherein, the third predicted soil attribute information includes the third predicted information of organic carbon content, the third predicted information of soil bulk density and the third predicted content of gravel larger than the preset diameter.

Step 015. According to the third soil attribute prediction information and thickness information of each fitting layer corresponding to each verification sample point position in the verification sample point position set, obtain the third soil organic value of each verification sample point position in the verification sample point position set. Predicted value of carbon density Constitute the third set of predicted values of soil organic carbon density PMTC-F for each verification sample point location in the verification sample point location set; then go to step 016.

The step 015 specifically includes the following operations:

According to the third soil attribute prediction information and thickness information of each fitting layer corresponding to each verification sample point position in the verification sample point position set, according to the following formula:

Obtain the third predicted value of soil organic carbon density for each verification sample point location in the verification sample point location set Constitute the third predicted value set PMTC-F of soil organic carbon density in each verification sample point location set in the verification sample point location set; where, Indicates the total number of fitting layers corresponding to the i _2th verification sample point position in the verification sample point position set, Indicates the kth fitting layer corresponding to the i ₂ verification sample point position in the verification sample point position set, Indicates the prediction information of the third organic carbon content of the k-th fitting layer corresponding to the i- _2th verification sample point position in the set of verification sample point positions, Indicates the third soil bulk density prediction information of the k-th fitting layer corresponding to the i- _2th verification sample point position in the verification sample point location set, Indicates the predicted content of the third gravel greater than the preset diameter in the k-th fitting layer corresponding to the i- ₂ verification sample point position in the verification sample point location set, Indicates the thickness information of the k-th fitting layer corresponding to the i- _2th verification sample point position in the verification sample point position set; then enter step 016.

Step 016. Use Lins's concordance correlation coefficient (Lins's concordance correlation coefficient) as the accuracy verification index to evaluate the accuracy of the prediction results:

in, Indicates the measured value of soil organic carbon density at the i _2th verification sample point location in the verification sample point location set, that is, the measured value of soil organic carbon density Indicates the predicted value of soil organic carbon density at the i _2th verification sample point location in the set of verification sample point locations, that is, corresponding to the first predicted value of soil organic carbon density, the second predicted value of soil organic carbon density, and the third predicted value of soil organic carbon density Predictive value; Indicates the measured values of soil organic carbon density at all verification sample locations in the verification sample location set average value; Indicates the predicted value of soil organic carbon density at all verification sample locations in the verification sample location set The average value of ; Lins's consistency correlation coefficient has nothing to do with the target data magnitude of the precision test, the value range is [-1,1], and the maximum value of 1 indicates the best prediction result, which is completely consistent with the observed data. In the process of practical application, the accuracy evaluation above is used respectively for the first set of predicted values of soil organic carbon density PMTC-D, the second set of predicted values of soil organic carbon density PCTM and the third set of predicted values of soil organic carbon density PMTC-F. method, wherein, when using the first soil organic carbon density predicted value set PMTC-D, each value in the first soil organic carbon density predicted value set PMTC-D corresponds to When the second soil organic carbon density predicted value set PCTM is used, each value in the second soil organic carbon density predicted value set PCTM corresponds to When using the third soil organic carbon density predicted value set PMTC-F, each value in the third soil organic carbon density predicted value set PMTC-F corresponds to

According to the accuracy test of the measured value set V of soil organic carbon density, the first set of predicted values of soil organic carbon density PMTC-D, the second set of predicted values of soil organic carbon density PCTM and the third set of predicted values of soil organic carbon density PMTC-F are obtained The optimal soil organic carbon density prediction value set in , and enter step 017.

Step 017. Discretize the spatial raster data of the target soil area, and use the sampling data of the occurrence layer corresponding to all sample points in the target soil area as the prediction data set. If the optimal soil organic carbon density prediction value set is the first soil For the set of predicted values of organic carbon density PMTC-D, the method from step 009 to step 011 is used to obtain the grid data of spatial distribution of soil organic carbon density in the target soil area; if the optimal set of predicted values of soil organic carbon density is the second soil organic carbon Density prediction value set PCTM, the method of step 012 is used to obtain the spatial distribution raster data of soil organic carbon density in the target soil area; if the optimal soil organic carbon density prediction value set is the third soil organic carbon density prediction value set PMTC-F , then use the method from step 013 to step 015 to obtain the grid data of the spatial distribution of soil organic carbon density in the target soil area; then go to step 018.

Step 018. According to the grid data of spatial distribution of soil organic carbon density in the target soil area, obtain the soil organic carbon storage of the target soil area.

The method for estimating soil organic carbon storage based on the above-mentioned design of the present invention based on the prediction of the thickness of the soil origin layer is applied in practice, such as taking the estimation of soil organic carbon storage in Liaoning Province as an example.

Liaoning Province is located in the south of Northeast China, bordering Jilin Province in the northeast and Inner Mongolia Autonomous Region in the northwest. Liaoning Province has insufficient land resources, fewer cultivated land resources, and many types of land use. Taking the soil survey data in 2009 as input, on the basis of the present invention, using different prediction methods and merging strategies of occurrence strata, the estimated value of organic carbon storage based on the optimal method can be obtained. In the actual application in Liaoning Province, as shown in Figure 7a to Figure 7c, for the prediction results of different prediction methods PMTC-D, PCTM, PMTC-F, Lins's concordance correlation coefficient (Lins's concordance correlation coefficient) is used as the accuracy verification indicators to evaluate the accuracy of prediction results. The accuracy evaluation results showed that CTM was the highest (ρ _c =0.21), followed by MTC-D (ρ _c =0.18), and MTC-F was the worst (ρ _c =0.13). Therefore, for the estimation of organic carbon storage in Liaoning Province, it is recommended to use the CTM method to predict the organic carbon storage in this area; then the target soil area is discretized spatial raster data, and all sample points in the target soil area correspond to the occurrence The sampling data of the layer is used as the prediction data set, and the method of step 012 is used to obtain the grid data of the spatial distribution of soil organic carbon density in the target soil area; finally, according to the grid data of the spatial distribution of soil organic carbon density in the target soil area, the target soil area is obtained soil organic carbon storage.

Different from conventional GIS-based methods for estimating soil organic carbon storage, the present invention fully considers the homogeneity of soil properties on the basis of occurrence layers and the spatial continuity of horizontal dimensions. The use of cross-validation provides the prediction accuracy and optimal calculation mechanism guarantee for the carbon storage estimation project in the target research area, and effectively overcomes the technical bottleneck of the limitations of a spatial prediction technology. This technology has very important guiding significance for scientifically formulating an effective management mechanism, realizing the sustainable development of resources and giving full play to the ecological benefits of the ecosystem.

The embodiments of the present invention have been described in detail above in conjunction with the accompanying drawings, but the present invention is not limited to the above embodiments, and can also be made without departing from the gist of the present invention within the scope of knowledge possessed by those of ordinary skill in the art. Variations.

Claims (10)

1. A soil organic carbon reserve estimation method based on soil occurrence layer thickness prediction is characterized by comprising the following steps:

001, setting each sampling point position aiming at the target soil area, counting the number of the sampling point positions to be I, and then entering the step 002;

step 002, in the target soil area, respectively aiming at each sampling point position, obtaining the environmental information of the sampling point position, and meanwhile, if the depth of the sampling point position in the soil area vertically downward is greater than or equal to the preset depth L, obtaining a profile area of the sampling point position vertically downward with the preset depth L; if the depth of the soil area with the sampling point vertically downward is smaller than the preset depth L, acquiring a profile area of the soil with the sampling point vertically downward; dividing soil generation layers aiming at the profile area according to the operation requirements of Chinese soil system classification to obtain each soil generation layer corresponding to the sampling point position; respectively obtaining thickness information and soil morphology information of each soil occurrence layer corresponding to the sampling point position; further, respectively obtaining environmental information of each sampling point position in the target soil area, and thickness information and soil morphology information of each soil occurrence layer corresponding to each sampling point position, and then entering step 003;

step 003, respectively collecting soil samples with preset quality in each soil occurrence layer corresponding to the sampling point positions aiming at the sampling point positions in the target soil area, respectively measuring and obtaining soil attribute information of each soil occurrence layer, further obtaining the soil attribute information of each soil occurrence layer corresponding to each sampling point position in the target soil area, and then entering step 004; the soil attribute information comprises soil organic carbon content, gravel content larger than a preset diameter and soil volume weight;

step 004, aiming at each sampling point position in the target soil area, respectively, obtaining soil organic carbon density measured value SOCD (soil organic carbon density) of each sampling point position according to soil attribute information and thickness information of each soil occurrence layer corresponding to the sampling point position_iI ═ 1, …, I }, and then proceed to step 005;

005, merging each soil occurrence layer corresponding to the sampling point position respectively aiming at each sampling point position in the target soil area to obtain each merging layer corresponding to the sampling point position, then carrying out weighted calculation aiming at the soil attribute information of each soil occurrence layer of the sampling point position according to each soil occurrence layer corresponding to each merging layer of the sampling point position respectively to obtain the soil attribute information of each merging layer corresponding to the sampling point position, and carrying out summation calculation aiming at the thickness information of each soil occurrence layer of the sampling point position to obtain the thickness of each merging layer corresponding to the sampling point position; further, the merging layers corresponding to the positions of the sampling points in the target soil area, and the soil property information and the thickness information of the merging layers are obtained, and the process then proceeds to step 006;

step 006, respectively aiming at each sampling point position in the target soil area, judging whether the sum of the thicknesses of the merging layers corresponding to the sampling point position is smaller than the preset depth L, if so, arranging the merging layers below the merging layers corresponding to the sampling point position by using a base rock layer to ensure that the sum of the thicknesses of the merging layers corresponding to the sampling point position is equal to the preset depth L, and entering step 007; otherwise, directly entering step 007;

step 007, acquiring the types of merging layers corresponding to all the positions of the sampling points in the target soil area to form a merging layer type set of the target soil area; according to the type set of the merging layer of the target soil area, carrying out unified operation on the merging layers corresponding to the positions of the sampling points in the target soil area, enabling the types of the merging layers corresponding to the positions of the sampling points in the target soil area to be the same, and entering the step 008;

step 008, dividing the positions of the sampling points in the target soil area into a prediction sampling point position set and a verification sampling point position set by adopting a linear congruence algorithm, and according to the soil organic carbon density measured value SOCD of each sampling point position in the target soil area_iObtaining the measured value of the organic carbon density of the soil at each verification sampling point position in the verification sampling point position setForming a soil organic carbon density measured value set V of each verification sampling point position in the verification sampling point position set; and proceeds to step 009, where i₂＝{1,…,I₂In which I₂In order to verify the number of the verification sampling point positions in the sampling point position set, the number of the prediction sampling point positions in the prediction sampling point position set is I₁，I₁>I₂；

Step 009, training to obtain each prediction model respectively taking the thickness information of each merging layer as a target by adopting a random forest method according to the environment information of each predicted sampling point position in the predicted sampling point position set and the thickness information of each merging layer corresponding to each predicted sampling point position respectively, and forming a first prediction model set;

meanwhile, according to the environmental information of each predicted sampling point position in the predicted sampling point position set and the soil attribute information of each merging layer corresponding to each predicted sampling point position, a random forest method is adopted to train and obtain each prediction model which takes each merging layer soil attribute information as a target, and a second prediction model set is formed; then, the step 010 is executed;

step 010, aiming at each verification sampling point position in the verification sampling point position set, respectively obtaining first prediction thickness information of each merging layer corresponding to the verification sampling point position respectively through each prediction model in the first prediction model set according to the environment information of the verification sampling point position; meanwhile, respectively aiming at each verification sampling point position in the verification sampling point position set, respectively obtaining first prediction soil attribute information of each merging layer corresponding to the verification sampling point position through each prediction model in the second prediction model set according to the environment information of the verification sampling point position; further acquiring first predicted thickness information and first predicted soil attribute information of each merging layer corresponding to each verification sampling point position in the verification sampling point position set; then entering step 011; the first predicted soil attribute information comprises first organic carbon content prediction information, first soil volume weight prediction information and first gravel content which is larger than a preset diameter;

step 011, respectively aiming at each verification sampling point position in the verification sampling point position set, obtaining a stretching coefficient Rc of the sum of the first predicted thickness information of each merging layer corresponding to the L relative verification sampling point position; then, according to first predicted soil attribute information and first predicted thickness information of each merging layer corresponding to each verification sampling point position in the verification sampling point position set, obtaining a first soil organic carbon density predicted value of each verification sampling point position in the verification sampling point position setComposing individual verifications in a set of verification sample point locationsA first soil organic carbon density prediction value set PMTC-D of a sampling point position; then, go to step 012;

step 012, training to obtain a third prediction model with the soil organic carbon density measured value as a target by adopting a random forest method according to the environmental information of each predicted sampling point position in the predicted sampling point position set and the soil organic carbon density measured value of each predicted sampling point position; then, respectively aiming at each verification sampling point position in the verification sampling point position set, obtaining a second soil organic carbon density predicted value of the verification sampling point position through a third prediction model according to the environmental information of the verification sampling point position, further obtaining the second soil organic carbon density predicted value of each verification sampling point position in the verification sampling point position set, and forming a second soil organic carbon density predicted value set PCTM of each verification sampling point position in the verification sampling point position set; then, go to step 013;

step 013, dividing each merging layer corresponding to the sampling point positions into each matching layer based on L according to a preset division rule aiming at all the sampling point positions in the target soil area, wherein the matching layers corresponding to the sampling point positions and the matching layers are the same in number, and obtaining thickness information of the matching layers; then, fitting the soil attribute information of each merging layer sampling occurrence layer of the predicted sampling point position according to each merging layer corresponding to each simulated layer of the predicted sampling point position respectively aiming at each predicted sampling point position in the predicted sampling point position set, and obtaining the soil attribute information of each simulated layer of the predicted sampling point position; further obtaining soil attribute information of each simulated layer corresponding to each predicted sampling point position in the set of predicted sampling point positions, and entering step 014;

step 014, training to obtain each prediction model taking soil attribute information of each fitting layer as a target by adopting a random forest method according to the environment information of each prediction sampling point position in the prediction sampling point position set and the soil attribute information of each fitting layer corresponding to each prediction sampling point position, and forming a fourth prediction model set; then, respectively aiming at each verification sampling point position in the verification sampling point position set, according to the environmental information of the verification sampling point position, obtaining third soil attribute prediction information of each fitting layer corresponding to the verification sampling point position respectively through each prediction model in the fourth prediction model set, further obtaining the third soil attribute prediction information of each fitting layer corresponding to each verification sampling point position in the verification sampling point position set respectively, and then entering step 015; the third predicted soil attribute information comprises third organic carbon content prediction information, third soil volume weight prediction information and third gravel content predicted to be larger than a preset diameter;

step 015, according to the third soil attribute prediction information of each fitting layer corresponding to each verification sampling point position in the verification sampling point position set and the thickness information, obtaining a third soil organic carbon density prediction value of each verification sampling point position in the verification sampling point position setForming a third soil organic carbon density predicted value set PMTC-F of each verification sampling point position in the verification sampling point position set; then step 016 is entered;

step 016, carrying out precision test according to the soil organic carbon density measured value set V to obtain an optimal soil organic carbon density predicted value set in a first soil organic carbon density predicted value set PMTC-D, a second soil organic carbon density predicted value set PCTM and a third soil organic carbon density predicted value set PMTC-F, and entering step 017;

step 017, discretizing spatial grid data of the target soil area, taking sampling data of occurrence layers corresponding to all sampling point positions in the target soil area as a prediction data set, and if the optimal soil organic carbon density prediction value set is the first soil organic carbon density prediction value set PMTC-D, obtaining soil organic carbon density spatial distribution grid data of the target soil area by adopting the method from step 009 to step 011; if the optimal soil organic carbon density prediction value set is the second soil organic carbon density prediction value set PCTM, acquiring soil organic carbon density spatial distribution grid data of the target soil area by adopting the method of the step 012; if the optimal soil organic carbon density prediction value set is a third soil organic carbon density prediction value set PMTC-F, obtaining soil organic carbon density spatial distribution grid data of the target soil area by adopting the method from the step 013 to the step 015; then proceed to step 018;

and 018, obtaining the soil organic carbon reserve of the target soil area according to the soil organic carbon density space distribution grid data of the target soil area.

2. The soil organic carbon storage estimation method based on soil occurrence layer thickness prediction as claimed in claim 1, wherein the step 001 specifically comprises the following steps:

00101, obtaining a soil type distribution map, a soil utilization distribution map and a soil geological distribution map of the target soil area, spatially stacking the soil type distribution map, the soil utilization distribution map and the soil geological distribution map to obtain a multilayer stack map of the target soil area, and then entering 00102;

00102, respectively obtaining the area of each pattern spot area in the multilayer superposition map of the target soil area, and counting to obtain each pattern spot area with the area ratio exceeding a pattern spot area threshold; then, for each of the pattern spot areas, sampling point positions are set in accordance with the real-area accessibility detection request, the number of sampling point positions is counted as I, and the process proceeds to step 002.

3. The soil organic carbon storage estimation method based on soil occurrence layer thickness prediction as claimed in claim 1, wherein: the environmental information comprises vegetation coverage, rock outcrop area ratio, terrain, important markers, ground surface coarse fragment size, ground surface fracture conditions and ground surface salt spot information.

4. The soil organic carbon storage estimation method based on soil occurrence layer thickness prediction as claimed in claim 1, wherein: the soil morphology information represents soil dry and wet conditions, soil color, root system information, pore information, sample structure, stripe composition substances, nodular tuberculosis substances, cementation degree and lime reaction information.

5. The soil organic carbon storage estimation method based on soil occurrence layer thickness prediction according to claim 1, wherein the step 004 specifically comprises the following operations:

respectively aiming at each sampling point position in the target soil area, and according to soil attribute information and thickness information of each soil occurrence layer corresponding to the sampling point position, the following formula is adopted:

${\mathrm{SOCD}}_i=\frac{\underset{n_i=1}{\overset{N_i}{\Σ}}({\mathrm{SOC}}_{n_i}\×{\mathrm{BD}}_{n_i}\×(1-{\mathrm{Gr}}_{n_i})\×T_{n_i})}{L}$

obtaining soil organic carbon density measured value SOCD of each sampling point position_iWherein i ═ {1, …,I}，n_i＝{1,…N_i}，n_iRepresents the nth soil occurrence layer corresponding to the ith sampling point position, N_iRepresenting the total number of soil occurrence layers corresponding to the ith sampling point position;showing the organic carbon content of the n-th soil occurrence layer corresponding to the ith sampling point position,showing the soil volume weight of the nth soil occurrence layer corresponding to the ith sampling point position,indicating the content of gravel with diameter larger than the preset diameter in the nth soil occurrence layer corresponding to the ith sampling point position,and (4) thickness information of the nth soil occurrence layer corresponding to the ith sampling point position is shown, and then the process goes to step 005.

6. The method for estimating the organic carbon content in soil based on the soil occurrence layer thickness prediction according to claim 1, wherein in the step 005, the merging process is performed for each soil occurrence layer corresponding to the sampling point position respectively for each sampling point position in the target soil area, and the process of obtaining each merging layer corresponding to the sampling point position specifically includes the following steps:

00501, aiming at each sampling point position in the target soil area, respectively, obtaining characteristic expression symbols corresponding to each soil occurrence layer of the sampling point position according to the soil morphological information of each soil occurrence layer corresponding to the sampling point position and Chinese soil system classification; further acquiring characteristic expression symbols corresponding to each soil occurrence layer of each sampling point position in the target soil area, and then entering step 00502;

00502, establishing a occurrence layer characteristic merging table according to Chinese soil system classification;

then, aiming at each sampling point position in the target soil area, aiming at each soil occurrence layer corresponding to the sampling point position, dividing according to a humus surface layer, a deposition layer and a matrix layer, and according to a soil occurrence layer characteristic expression symbol corresponding to each merging layer in a generation layer characteristic merging table, sequentially dividing each soil occurrence layer corresponding to the sampling point position from top to bottom in a section area of the sampling point position to obtain each merging layer corresponding to the sampling point position;

meanwhile, in the process of performing the merging operation on each sampling point position in the target soil area, if the same soil occurrence layer characteristic expression symbol appears in two merging layers corresponding to the sampling point position, the step 00503 is performed; 00503, judging whether the soil occurrence layer characteristic expression symbols in the first merging level from top to bottom in the occurrence layer characteristic merging table exist in the soil occurrence layer characteristic expression symbols corresponding to the sampling point positions, if so, merging the soil occurrence layers corresponding to the same soil occurrence layer characteristic expression symbols respectively, and then performing 00502 to merge to obtain the merging layers corresponding to the sampling point positions respectively; otherwise, entering step 00504;

00504 dividing the soil occurrence layer characteristic expression symbols in the lower merging layer in the same soil occurrence layer characteristic expression symbols into the merging layer where the same soil occurrence layer characteristic expression symbols are located, and moving each soil occurrence layer characteristic expression symbol corresponding to the sampling point position upwards by a merging layer according to the merging layer sequence in the occurrence layer characteristic merging table to obtain each merging layer corresponding to the sampling point position.

7. The method for estimating organic carbon storage in soil according to claim 1, wherein in the step 007, according to the set of types of merging layers of the target soil region, a unified operation is performed on the merging layers corresponding to the positions of the sampling points in the target soil region, so that the types of the merging layers corresponding to the positions of the sampling points in the target soil region are the same, and specifically, the method includes the following operations:

respectively judging whether the merging layer corresponding to the sampling point position is equal to the merging layer type set of the target soil area or not aiming at each sampling point position in the target soil area, if so, not doing any operation; otherwise, setting a merging layer which is lack relative to the merging layer type set of the target soil area to the sampling point position, setting the thickness information of the merging layer to be 0.00001cm, and setting the soil attribute information of the merging layer to be 0.00001; and then the merging layers respectively corresponding to the positions of the sampling points in the target soil area are equal to the type set of the merging layers in the target soil area.

8. The soil organic carbon storage estimation method based on soil occurrence layer thickness prediction according to claim 1, wherein the step 011 specifically comprises the following operations:

respectively aiming at each verification sampling point position in the verification sampling point position set, obtaining a stretching coefficient Rc of the sum of the first predicted thickness information of each merging layer corresponding to the L relative verification sampling point position; then according to the first predicted soil attribute information and the first predicted thickness information of each merging layer corresponding to each verification sampling point position in the verification sampling point position set, the following formula is adopted:

${\mathrm{SOCD}}_{i_2}^{\′}=\frac{\underset{n_{i_2}=1}{\overset{N_{i_2}}{\mathrm{\Σ}}}\left({\mathrm{SOC}}_{n_{i_2}}^{\′}\×{\mathrm{BD}}_{n_{i_2}}^{\′}\×\left(1-{\mathrm{Gr}}_{n_{i_2}}^{\′}\right)\×Rc\×T_{n_{i_2}}^{\′}\right)}{L}$

obtaining a first soil organic carbon density predicted value of each verification sampling point position in the verification sampling point position setForming a first of each verified sample position in a set of verified sample positionsCollecting a soil organic carbon density predicted value PMTC-D; wherein, representing the ith in a verified sample location set₂The total number of merging layers corresponding to each verified sampling point position,representing the ith in a verified sample location set₂The nth merging layer corresponding to the verification sampling point position,representing the ith in a verified sample location set₂The first organic carbon content prediction information of the nth merging layer corresponds to the position of each verification sampling point,representing the ith in a verified sample location set₂The first soil volume weight prediction information of the nth merging layer corresponding to the position of each verification sampling point,representing the ith in a verified sample location set₂The position of each verification sampling point corresponds to the predicted content of the first gravel with the diameter larger than the preset diameter in the nth merging layer,representing the ith in a verified sample location set₂The first predicted thickness information of the nth merging layer corresponds to the verification sampling point position; then, the process proceeds to step 012.

9. The method for estimating soil organic carbon reserves based on soil occurrence layer thickness prediction according to claim 1, wherein in step 013, soil attribute information of each simulated layer corresponding to each predicted sampling point position in the set of predicted sampling point positions is obtained, specifically comprising the following operations:

aiming at each predicted sampling point position in the predicted sampling point position set, the following formula is adopted according to an equal-area Spline function:

$S\left(k_{i_1}\right)=\stackrel{\‾}{f}\left(k_{i_1}\right)+e\left(k_{i_1}\right)$

obtaining soil attribute information of each fitting layer of the predicted sampling point position; and further acquiring soil attribute information of each fitting layer corresponding to each predicted sampling point position in the predicted sampling point position set, wherein i₁＝{1,…,I₁}，I₁Representing the number of predicted sample positions in the set of predicted sample positions, representing the ith in a set of predicted sample positions₁The number of layers of the fitting layer corresponding to the position of each predicted sampling point,representing the ith in a set of predicted sample positions₁Soil of k-th layer fitting layer corresponding to each predicted sampling point positionThe information of the soil property is obtained,to representFunction is asLayer andthe average of the results of the layer fit,representing the ith in a set of predicted sample positions₁The laboratory measurement analysis error of the k-th fitting layer corresponding to the position of each predicted sampling point,the function is an equal area Spline function.

10. The soil organic carbon storage estimation method based on soil occurrence layer thickness prediction as claimed in claim 1, wherein the step 015 specifically comprises the following operations:

according to the third soil attribute prediction information of each fitting layer corresponding to each verification sampling point position in the verification sampling point position set and the thickness information, the following formula is adopted:

${\mathrm{SOCD}}_{i_2}^{\′\′}=\frac{\underset{k_{i_2}=1}{\overset{K_{i_2}}{\mathrm{\Σ}}}\left({\mathrm{SOC}}_{k_{i_2}}^{\′\′}\×{\mathrm{BD}}_{k_{i_2}}^{\′\′}\×\left(1-{\mathrm{Gr}}_{k_{i_2}}^{\′\′}\right)\×Rc\×T_{k_{i_2}}^{\′\′}\right)}{L}$

obtaining a third soil organic carbon density predicted value of each verification sampling point position in the verification sampling point position setForming a third soil organic carbon density predicted value set PMTC-F of each verification sampling point position in the verification sampling point position set; wherein, representing the ith in a verified sample location set₂The total number of the fitting layers corresponding to the positions of the verification sampling points,representing the ith in a verified sample location set₂The k-th fitting layer corresponding to the verification sampling point position,representing the ith in a verified sample location set₂The third organic carbon content prediction information of the kth fitting layer corresponding to the verification sampling point position,representing the ith in a verified sample location set₂The third soil volume weight prediction information of the kth fitting layer corresponding to the verification sampling point position,representing a set of verified sample locationsMiddle (i)₂The position of each verification sampling point corresponds to the predicted content of the third gravel with the preset diameter in the kth fitting layer,representing the ith in a verified sample location set₂Thickness information of a kth fitting layer corresponding to the verification sampling point position; then proceed to step 016.