A method for extracting soil carbon cycle characteristics in karst forests based on multimodal data fusion

Abstract

The present invention discloses a method for extracting carbon cycle characteristics of karst forest soil by multimodal data fusion, including obtaining soil physical and chemical protection component carbon data, biological active component carbon data and isotope tracer data at different succession stages of karst forest to establish multi-source heterogeneous data; integrating the obtained multi-source heterogeneous data into the "physical-chemical-biological-isotope" four-dimensional feature space through feature tensor decomposition; constructing a carbon saturation capacity prediction hybrid optimization model, outputting the maximum carbon sequestration capacity and component carbon proportion of microaggregates <53μm; mineralizing and cultivating soil at different succession stages of karst forest, combining with the above methods, to obtain the carbon cycle characteristics of karst forest soil by multimodal data fusion. 13 C / 15 The method uses a time-varying calibration curve for a double-standard sample of N, analyzes the dynamic ratio of exogenous carbon to native soil carbon, analyzes the input and output pathways of soil carbon, and calculates the carbon stability perturbation threshold. This method enables the collaborative analysis of multi-source data, significantly improving the temporal and spatial resolution of carbon cycle feature extraction and model accuracy.

Description

Technical Field

[0001] The present invention relates to the fields of environmental science and agricultural ecological technology, and specifically to a karst forest soil carbon cycle feature extraction method based on multimodal data fusion, which is suitable for karst ecosystem carbon sequestration potential assessment, carbon cycle path analysis and sustainable carbon management strategy formulation. Background Art

[0002] Soil carbon is a core component of the terrestrial ecosystem carbon cycle, and its dynamics directly impact global climate change and soil quality. Karst forests, a crucial ecosystem in southern my country, possess soils characterized by high heterogeneity, low carbon storage, and susceptibility to erosion. Traditional studies have relied on single data sources (such as physical component analysis or static models), making it difficult to fully understand the complex mechanisms of the carbon cycle.

[0003] Existing methods often rely on single-dimensional data (such as soil physical fractionation or analysis of chemically protected components), ignoring the synergistic effects of multimodal data (such as hyperspectral remote sensing, microbial communities, bioactive carbon components, and isotope tracing). For example, estimating carbon storage solely based on soil bulk density and total organic carbon content fails to account for differences in carbon sequestration among different aggregates (e.g., >250μm and <53μm) and the effects of microbial regulation, resulting in insufficient quantification accuracy of carbon cycle pathways.

[0004] Traditional carbon cycle models, such as the Century model, are mostly static empirical formulas that lack the ability to respond in real time to dynamic environmental factors (such as temperature, humidity, and pH) and spatial heterogeneity. For example, existing models struggle to capture the carbon saturation capacity boundary effect caused by topographical fluctuations in karst regions, and they fail to combine boundary effect methods with least squares optimization for hybrid optimization, limiting the accuracy of carbon sequestration potential predictions.

[0005] although 13 C isotope labeling technology has been used to analyze carbon sources, but existing methods mostly rely on single isotopes (such as 13 C) and lack of dynamic calibration. For example, traditional mineralization experiments only analyze the ratio of exogenous carbon to native soil carbon, without considering 15 The construction of a time-varying calibration curve using N double standards results in large errors in the quantization calculation of excitation effects (such as positive excitation and negative excitation), making it impossible to accurately quantify the carbon stability perturbation threshold.

[0006] Existing research has limited understanding of the role of microbial communities in regulating carbon cycling, often relying on simple statistical analyses of community abundance without constructing dynamic regulatory network models. For example, high-throughput sequencing data combined with graph convolution algorithms have not been used to capture microbial synergies. Consequently, the spatiotemporal dynamics of microbial-carbon interactions cannot be revealed, hindering the development of intelligent carbon cycle monitoring and prediction.

[0007] In summary, existing technologies struggle to meet the needs of in-depth research on the carbon cycle in karst forests due to their limited data dimensionality, insufficient model dynamics, low isotope analysis precision, and limited understanding of microbial mechanisms. This present invention overcomes these technical bottlenecks through multimodal data fusion and dynamic modeling. Summary of the Invention

[0008] This invention aims to address the challenges of studying soil carbon cycling during karst forest restoration. By integrating multimodal data (carbon data from soil physical protection components, chemical protection components, biologically active components, and isotope tracer data), this method constructs a four-dimensional feature space and dynamic model, enabling precise quantification of soil carbon cycle pathways, saturation capacity, and priming effects. Combined with carbon-nitrogen coupling kinetic analysis and digital twin simulation, this approach provides a scientific basis for assessing the carbon sink potential of karst ecosystems, early warning of carbon instability risks, and the development of sustainable carbon management strategies.

[0009] The present invention is implemented as follows: a multimodal data fusion method for extracting carbon cycle characteristics of karst forest soil, comprising the following steps:

[0010] S1: Obtain carbon data of soil physical protection components, chemical protection components, biologically active components, and isotope tracer data at different succession stages of karst forests to establish multi-source heterogeneous data;

[0011] S2, integrating the multi-source heterogeneous data obtained in step S1 into the four-dimensional feature space of "physical-chemical-biological-isotopic" through feature tensor decomposition, and dynamically adjusting the contribution of each modal data using an adaptive weight allocation algorithm;

[0012] S3, based on the spatial constraints of the boundary effect method and the statistical constraints of the least squares method, takes the four-dimensional feature tensor and weight distribution results output in step S2 as input, constructs a hybrid optimization model for predicting carbon saturation capacity, and outputs the maximum carbon sequestration capacity of microaggregates <53 μm and the carbon proportion of the components;

[0013] S4, based on the maximum carbon capacity predicted in step S3, a mineralization culture experiment was designed, and the soils of different succession stages of the karst forest were mineralized using isotope tracing. 13 C / 15 The time-varying calibration curve of the N double standard sample analyzes the dynamic ratio of exogenous carbon to native carbon and the input and output pathways of soil carbon, and calculates the carbon stability disturbance threshold.

[0014] In the above scheme, step S2 maps the multi-source heterogeneous data (physical, chemical, biological, and isotopic) obtained in the previous step into a unified structured space, resolving data heterogeneity (unit, dimension, and format differences) and facilitating subsequent model processing. An adaptive algorithm dynamically adjusts the contribution of each modality (physical, chemical, biological, and isotopic). For example, chemical protection components are weighted more highly in the climax community stage, while physical protection components are weighted even more highly in the shrub stage, ensuring that the model focuses on key influencing factors.

[0015] The four-dimensional feature tensor provides high-information-density input after dimensionality reduction, avoiding redundancy in the original data. Furthermore, it can guide model optimization in the subsequent step S3 through weights. This is because the weight distribution directly affects the penalty coefficients of each mode in the optimization model. For example, high-weight modes occupy a larger proportion in the objective function, and the constraint model pays more attention to their contribution.

[0016] Furthermore, during the experimental design phase in step S4, the weight assignment results guide the addition ratio of isotope standards (e.g., high-weighted isotope modes require higher sampling frequencies). They also facilitate dynamic threshold adjustment. For example, when the physical mode weight is high, the stability threshold is more stringent (due to the strong physical protection effect).

[0017] The reasons why step S4 needs to be based on the maximum carbon capacity of S3 are mainly due to the following aspects:

[0018] First, based on ecological rationality, carbon capacity overload can be avoided, because exceeding the maximum capacity will lead to carbon saturation, causing unnatural loss and distorting the experimental results; considering that the experiment simulates real scenarios, the experimental load needs to reflect the carbon sequestration potential of the actual ecosystem, for example, the soil of the climax community can carry more exogenous carbon than the bare rock period.

[0019] Secondly, it is conducive to parameter standardization, unifying the maximum capacity as the benchmark to ensure that the experimental results of different succession stages are comparable (for example, adding 5% in the shrub stage vs. adding 5% in the climax community represent different absolute amounts); and it is conducive to dynamic threshold association. The intensity of the excitation effect needs to be adjusted according to the capacity ratio. For example, when the addition amount is 0.3x3, the threshold is reduced by 20%.

[0020] Finally, based on technical necessity, when it comes to calibration curve accuracy, the interpretation of isotope signals (such as δ¹³C) needs to respond linearly within the carbon capacity range. Exceeding the range may lead to nonlinear errors. In addition, considering the authenticity of microbial responses, the microbial metabolic rate is affected by the carbon saturation state, and the experimental design needs to simulate the community behavior under real carbon loads.

[0021] As a further preferred solution, it also includes: S5, constructing a digital twin of the carbon cycle, embedding the carbon stability disturbance threshold and input and output paths of step S4, and simulating the soil carbon saturation deficit value and spatial heterogeneity in different succession stages of karst forests (grassland stage, shrub stage, and tree stage).

[0022] The reasons why this step is not necessary in the present invention are as follows:

[0023] If a study only needs to assess the carbon capacity of a particular successional stage (such as the shrub stage) and does not require dynamic simulation across these stages, the complex modeling of a digital twin may be redundant. Furthermore, building a digital twin is also limited by data and resources. High-precision digital twins require multidimensional, real-time data (such as high-frequency isotope monitoring and drone remote sensing). If data acquisition costs are prohibitive or technical requirements are insufficient, a static model can be used.

[0024] If only carbon saturation deficit prediction is needed, traditional empirical models (such as the Century model) or regression analysis can meet this requirement, eliminating the need for digital twin technology. If the management objective is limited to a specific region (such as a single peak-cluster depression), the outputs of steps S3-S4 (maximum capacity, threshold) can be used to directly formulate strategies, eliminating the need for global simulation.

[0025] On the other hand, building digital twins requires significant computing resources (e.g., GPU clusters and cloud platforms) and ongoing calibration (e.g., annual updates of microbial community data), which may not be economical for small and medium-sized research teams. As carbon management strategies approach optimality, the improvements brought by digital twin simulations may be limited and difficult to justify their development costs.

[0026] However, the inclusion of this step can also have significant positive effects. First, the soil carbon cycle in karst forests involves multimodal interactions involving physics, chemistry, biology, and isotopes. Digital twins can integrate all the data and models from steps S1-S4, enabling dynamic simulation of carbon flows across scales (micro-macro), revealing nonlinear and hysteresis effects. Furthermore, microtopographic variations in karst landforms (such as peak clusters and karst caves) significantly influence carbon distribution. Digital twins use three-dimensional spatial mapping (such as kriging interpolation and point cloud rendering) to quantify the spatial variation of carbon saturation deficit values, guiding precise management.

[0027] Secondly, because carbon cycle mechanisms differ significantly between grassland, shrub, and tree stages, digital twins can dynamically update model parameters (such as microbial diversity and aggregate structure) to adapt to changes in carbon sequestration during succession. This can provide a simulation platform for policy design such as carbon trading and ecological compensation, for example, by assessing the economic costs and ecological benefits of sequestering one ton of carbon.

[0028] In step S5, the carbon stability perturbation threshold from step S4 is embedded, enabling real-time simulation of carbon instability risks under different fertilization or land use scenarios, generating a risk heat map (e.g., red alert zones are concentrated in foothills). Based on the input-output pathway (carbon loss / fixation-dominated), the carbon sequestration effects of increasing biochar application or adjusting vegetation cover can be simulated, providing a quantitative basis for decision-making in ecological restoration.

[0029] In most cases, in step S1, the soil physical protection component carbon data is the 53-250 μm carbon content, including the physical protection component carbon iPOM;

[0030] The chemical protection component carbon data is the carbon content of <53μm, including non-acid-decomposable free powder component carbon NH-dSilt, non-acid-decomposable free clay component carbon NH-dClay, non-acid-decomposable closed accumulation powder component carbon NH-μSilt, non-acid-decomposable closed accumulation clay component carbon NH-μClay, acid-decomposable free powder component carbon H-dSilt, acid-decomposable free clay component carbon H-dClay, acid-decomposable closed accumulation powder component carbon H-μSilt and acid-decomposable closed accumulation clay component carbon H-μClay;

[0031] The soil bioactive component carbon data refers to the carbon content >250 μm, including coarse particle component carbon cPOM and free particle component carbon fPOM;

[0032] The isotope tracing data includes 13 C abundance, 15 N abundance, the content and proportion of physical protection component carbon, chemical protection component carbon, and biologically active component carbon derived from original soil carbon and exogenous carbon, the direction and intensity of the excitation effect, and the mineralization amount and mineralization rate.

[0033] As a further preferred solution, in step S2, the feature tensor decomposition adopts Tucker decomposition or PARAFAC decomposition to decompose the multi-source heterogeneous data into core tensors and factor matrices, and the data is mapped to the "physical-chemical-biological-isotopic" four-dimensional feature space through the transformation of the factor matrix.

[0034] As a further preferred solution, the carbon saturation capacity prediction hybrid optimization model in step S3 is constructed according to the following steps:

[0035] 1) receiving the four-dimensional feature tensor and the adaptive weight distribution result outputted from step S2, and extracting the physical modal feature vector and the chemical modal feature vector;

[0036] 2) Generate a particle size constraint matrix based on the physical modal eigenvectors and calculate the theoretical saturation threshold of carbon content in each particle size based on the chemical modal eigenvectors and weights;

[0037] 3) Construct the following objective function:

[0038] ;

[0039] Where A is the particle size constraint matrix, =[ 1, 2, 3] T is the particle size carbon capacity prediction variable, b is the measured total organic carbon vector, is the penalty coefficient, is the theoretical saturation threshold of carbon content in each particle size, is the particle size classification index parameter, =1, 2, 3, corresponding to the particle size index parameters of >250μm, 53-250μm, and <53μm, respectively. For the Predictors of carbon capacity of aggregates;

[0040] 4) A hybrid strategy of particle swarm optimization and genetic algorithm was used to output the maximum carbon capacity and component carbon proportion of microaggregates <53 μm.

[0041] As a further preferred solution, the step S4 13 C / 15 The method for constructing a time-varying calibration curve of N dual standard samples includes:

[0042] Inject into the mineralization culture system in a 12-hour cycle 13 C-glucose and 15 N-ammonium sulfate mixed standard sample, each injection time is 5 minutes, its concentration is 15% of the original soil carbon content and 3% of the nitrogen content;

[0043] Gas samples were collected 1, 3, 6, and 12 hours after each injection and measured by isotope ratio mass spectrometry. 13 C and 15 N abundance;

[0044] A carbon-nitrogen coupling kinetic model was established; 13 C abundance uses an exponential decay equation to characterize carbon turnover. 15 The logarithmic growth equation was used to characterize nitrogen fixation for N abundance;

[0045] The carbon turnover rate constant and nitrogen retention characteristic time were fitted simultaneously by the simulated annealing-least squares hybrid algorithm to complete the construction. 13 C / 15 Time-varying calibration curve of N dual standards.

[0046] As a further preferred solution, the carbon-nitrogen coupling kinetic model includes the following coupling equations:

[0047] 13 C decay equation: ;

[0048] in, is the carbon turnover rate constant, which is solved by inversion using the simulated annealing algorithm;

[0049] 15 N adsorption equation: ;

[0050] in, is the characteristic time of nitrogen fixation, 、 is the soil type parameter.

[0051] As a further preferred solution, in step S4, based on the maximum carbon capacity of microaggregates <53 μm output in step S3, exogenous 13 C-labeled glucose and 15 N-labeled ammonium sulfate for mineralization culture;

[0052] During the mineralization incubation process, gas and soil samples were collected every 12 hours and measured simultaneously by isotope ratio mass spectrometry. 13 C and 15 N dynamic value; establish 13 C decay and 15 The double standard time-varying calibration curve coupled with N adsorption separates the contribution of exogenous carbon and the interference of nitrogen transformation;

[0053] Based on the calibration curve, the exogenous carbon input rate and the native carbon output rate are calculated, and the carbon flow balance ratio is defined. The carbon net loss path and the carbon net fixation path are determined by the carbon flow balance ratio.

[0054] Based on the dynamic monitoring data of mineralization rate, the difference between the mineralization rate of exogenous carbon and native carbon was integrated over time to generate the cumulative stimulation effect, which was then normalized by combining it with the microbial diversity index.

[0055] The carbon stability perturbation threshold is determined based on the proportion of carbon in the component <53 μm, and whether the carbon pool has entered an unstable state is determined by the carbon stability perturbation threshold and the normalized cumulative excitation effect;

[0056] The microbial diversity index includes the Shannon index based on 16S rRNA sequencing.

[0057] As a further preferred solution, the digital twin construction in step S5 includes:

[0058] Establish a three-dimensional carbon pool interactive map, and map the physical protection component carbon, chemical protection component carbon, and biological active component carbon in three-dimensional space using vector fields according to the different succession stages of herbaceous, shrubby, and tree soils.

[0059] Embedded in the LSTM-Transformer hybrid neural network, the input layer receives the cumulative carbon input, carbon mass of components <53μm and carbon content of components <53μm. The output layer predicts the carbon saturation deficit value and calculates the carbon saturation capacity time limit based on the existing carbon content and carbon retention rate.

[0060] As a further preferred solution, the hybrid strategy of particle swarm optimization and genetic algorithm to output the maximum carbon capacity and component carbon proportion of microaggregates <53 μm includes the following steps:

[0061] Generate multiple sets of initial solutions, limit the search range in the computer system, and define the weighted fitness function as follows:

[0062] ;

[0063] in, For fitness, is the physical mode weight, is the chemical modal weight;

[0064] By making the particle swarm perform gradient search in the physical modal space along the direction of the physical modal eigenvector, the genetic algorithm performs arithmetic crossover and Gaussian mutation on the chemical modal parameters to complete collaborative evolution; when the fitness change rate of the optimal solution is <0.05% for five consecutive generations, the algorithm is terminated and the result is output.

[0065] Compared with the current state of the art, which relies heavily on a single data dimension (such as only physical particle size or static total carbon content) and has difficulty analyzing the multi-factor coupling mechanism of the carbon cycle, the present invention realizes the collaborative analysis of multi-source data, significantly improving the spatiotemporal resolution and model accuracy of carbon cycle feature extraction.

[0066] First, this invention pioneers a four-dimensional feature space fusion technique combining "physical-chemical-biological-isotopic" features. This technique integrates particle size fractions (e.g., iPOM, H-μClay), chemically protected forms (e.g., NH-dSilt), bioreactive carbon (cPOM / fPOM), and dual-standard isotopic data through Tucker / PARAFAC tensor decomposition. Furthermore, an adaptive weighting algorithm is developed to dynamically allocate the contribution of each modality (e.g., increasing the chemical weight of the climax community stage to 0.6). This collaborative modeling of multi-source heterogeneous data significantly improves the accuracy of carbon sequestration predictions compared to traditional single-modal methods, addressing the model distortion caused by the high heterogeneity of karst soils.

[0067] Secondly, most existing carbon capacity models use fixed empirical formulas (such as the Century model), ignoring the synergy between boundary effects and statistical constraints. The hybrid optimization model constructed by this invention innovatively integrates the three-dimensional constraint surface of the boundary effect method with the least squares statistical fitting, and uses the dynamic penalty coefficient to ( =0.2+0.6*(1- )) Implement physical weight ( ) driven by adaptive adjustment of constraint strength. For example, when the physical weight of the shrub stage =0.3, =0.62 Strengthen boundary constraints to prevent overload prediction; while the top community =0.5, =0.5 focuses on statistical fitting.

[0068] On the other hand, traditional excitation effect analysis is limited to qualitative descriptions. This method, through a time-varying calibration curve of a 13C / 15N dual standard sample combined with a weighted integral of the microbial diversity index (Shannon index), accurately analyzes the dynamic ratio of exogenous carbon to native carbon during mineralization. Through a carbon-nitrogen coupling kinetic model and an excitation effect quantization algorithm, the traditional qualitative description of the excitation effect is converted into a calculable integral (such as the cumulative excitation effect), and isotope offsets are combined to eliminate cross-interference. This method overcomes the limitations of monoisotopic analysis, reducing the calculated error of the carbon stability perturbation threshold to within ±5%, significantly improving the accuracy of carbon cycle mechanism research. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 A flowchart of a working process in one embodiment of the present invention is shown. DETAILED DESCRIPTION

[0070] The preferred embodiments of the present invention will be described in detail below so that the purpose, features and advantages of the present invention can be more clearly understood. It should be understood that the following embodiments are not intended to limit the scope of the present invention, but are only intended to illustrate the essential spirit of the technical solution of the present invention.

[0071] In the following description, for the purpose of illustrating the various disclosed embodiments, certain specific details are set forth in order to provide a thorough understanding of the various disclosed embodiments. However, those skilled in the relevant art will recognize that the embodiments may be practiced without one or more of these specific details. In other cases, well-known techniques associated with this application may not be shown or described in detail to avoid unnecessarily obscuring the description of the embodiments.

[0072] Reference throughout this specification to "one embodiment" or "an embodiment" means that a particular feature, structure, or characteristic described in connection with the embodiment is included in at least one embodiment. Thus, the appearances of "in one embodiment" or "in an embodiment" in various places throughout this specification are not necessarily all referring to the same embodiment. Furthermore, the particular features, structures, or characteristics may be combined in any manner in one or more embodiments.

[0073] like Figure 1 As shown in FIG, a multimodal data fusion method for extracting soil carbon cycle characteristics in karst forests includes the following steps:

[0074] S1. Obtain carbon data of soil physical protection components, chemical protection components, biologically active components and isotope tracer data at different succession stages of karst forests to establish multi-source heterogeneous data.

[0075] In practice, sampling areas were established at different succession stages of the karst forest (e.g., herbaceous, shrubby, and arborescent). Multiple sampling points were set up in each area according to a grid layout to ensure that the samples were representative of the soil characteristics of that succession stage. Soil samples were collected at a depth of 0-20 cm at each sampling point. Samples from multiple points within the same area were combined into a single composite sample to reduce sampling error. Three to five composite samples were collected from each succession stage. Visible impurities such as plant debris, roots, and rocks were removed, and the samples were brought back to the laboratory to air-dry. Then, they were gently ground with a grinding rod and passed through a 2 mm sieve for later use.

[0076] (1) Acquisition of carbon data of physical protection components

[0077] The sieved soil is then subjected to wet sieving to separate particles into three sizes: >250μm, 53-250μm, and <53μm. This involves placing a certain amount of soil through sieves of varying pore sizes and shaking them in water for a specified period of time to separate the particles.

[0078] The organic carbon content of soil particles of each size fraction was determined using the potassium dichromate volumetric method coupled with external heating. The specific steps are: a certain amount of soil sample was weighed and placed into a test tube. A certain amount of potassium dichromate solution and concentrated sulfuric acid were added. The sample was heated in an oil bath for a predetermined period of time. After cooling, the remaining potassium dichromate was titrated with a standard ferrous sulfate solution. The soil organic carbon content was calculated based on the titration results.

[0079] iPOM, or physically protected organic matter, is the carbon component found within microaggregates (particle size 53-250 μm). Because it is physically protected within microaggregates, its decomposition rate is relatively slow, resulting in a relatively high stability within the soil carbon pool and a significant role in the long-term storage of soil carbon.

[0080] In practice, density flotation is used to separate the iPOM carbon fraction. A soil sample is placed in a solution of a certain density (such as NaI solution) and centrifuged. The iPOM carbon fraction floats to the surface of the solution. After collection and washing, its carbon content is determined using the potassium dichromate volumetric method coupled with external heating, as described above.

[0081] (2) Acquisition of carbon data of chemical protection components

[0082] Chemically protected component carbon includes chemically protected carbon and biochemically protected carbon. Chemically protected carbon is separated into H-silt carbon and H-clay carbon by chemical extraction. The soil sample is mixed with a certain concentration of sodium hydroxide solution, shaken for a certain period of time, and then the supernatant and precipitate are centrifuged. The supernatant mainly contains H-silt carbon and H-clay carbon. The carbon content in the supernatant is determined by the potassium dichromate volumetric method-external heating method, which is the total amount of chemically protected carbon. In order to distinguish H-silt carbon from H-clay carbon, further classification methods can be used, such as filtering through filter membranes with different pore sizes to separately determine the carbon content of different particle size fractions.

[0083] H-silt carbon is carbon that is chemically bound to silt particles and protected by chemical reactions. Silt particles have a particle size between sand and clay. This binding gives organic carbon a certain stability in the soil, making it less susceptible to microbial decomposition. H-clay carbon is chemically bound to clay particles. Clay particles have a large surface area and strong adsorption capacity, allowing them to bind tightly to organic carbon through chemical bonds, further enhancing its stability.

[0084] Biochemically protected carbon is separated from NH-silt carbon and NH-clay carbon using enzymatic hydrolysis combined with chemical extraction. A soil sample is first mixed with a specific enzyme solution (such as cellulase or protease) and reacted for a specified time at a specific temperature and pH. The sample is then extracted with sodium hydroxide solution, and the supernatant is separated by centrifugation. The carbon content in the supernatant is determined using the potassium dichromate volumetric method followed by external heating, representing the total amount of biochemically protected carbon. Similarly, further classification methods can be used to distinguish between NH-silt carbon and NH-clay carbon.

[0085] NH-silt carbon and NH-clay carbon are carbon bound to silt and clay, respectively, and protected by biochemical processes. This biochemical protection likely involves interactions between microbial metabolites and soil particles and organic carbon, forming a more complex structure that makes organic carbon more difficult to decompose. This is crucial for the long-term storage of soil carbon and the maintenance of soil fertility.

[0086] (3) Acquisition of carbon data of bioactive components

[0087] Bioactive carbon includes cPOM and fPOM. cPOM carbon generally refers to coarse organic carbon, which is relatively large and primarily derived from fresh plant debris. It has high bioactivity, is easily decomposed and utilized by microorganisms, and has a rapid turnover in the soil carbon cycle. fPOM carbon, which is free organic matter within aggregates—that is, carbon in particles larger than 250 μm—is also highly active and can quickly provide energy and nutrients to soil microorganisms.

[0088] Physical separation method combined with chemical analysis method was used to separate and determine the bioactive component carbon. That is, the soil sample was passed through a 53μm sieve, and the part above the sieve was cPOM component carbon and fPOM component carbon. The separated cPOM component carbon and fPOM component carbon were determined by potassium dichromate volumetric method-external heating method to determine their carbon content.

[0089] In order to distinguish the cPOM component carbon from the fPOM component carbon, they can be further separated according to their physical properties (such as color, texture, etc.), and then their carbon contents can be measured separately.

[0090] (4) Isotope tracer data acquisition

[0091] Respectively 13 C / 15 50 g of N-doubly labeled soil samples of herbs, shrubs, and trees were placed in a 50 mL beaker and injected into the mineralization culture system at a cycle of 12 h. 13 C-glucose and 15 The N-ammonium sulfate mixed standard sample was injected for 5 minutes each time, and its concentration was 15% of the original soil organic carbon content and 3% of the total nitrogen content; gas samples were collected 1, 3, 6, and 12 hours after each injection and determined by isotope ratio mass spectrometry. 13 C and 15 N abundance; soil samples were collected at the same time, part of which was placed in a -70℃ refrigerator for the determination of soil microorganisms, and part was dried and ground for soil carbon grouping.

[0092] The mineralized soil was separated into physical protection components, chemical protection components, and biologically active components by wet sieving, and the content of each soil component was determined by elemental analyzer-isotope ratio mass spectrometry (EA-IRMS). 13 C abundance and 15 N abundance. The specific operation is to burn the soil sample at high temperature in an elemental analyzer to convert the carbon and nitrogen in it into carbon dioxide and nitrogen. The generated gas is then introduced into the isotope ratio mass spectrometer through a gas transmission system to measure its isotope abundance.

[0093] 13 C abundance refers to the stable isotope of carbon in the soil 13 The proportion of C in the total amount of carbon. In the study of soil carbon cycle, carbon from different sources (such as plant residues, microbial metabolites, etc.) has different 13 C abundance characteristics were determined by measuring the 13 Changes in C abundance can track the source, transformation, and migration of carbon. 15 N abundance refers to the stable isotopes of nitrogen in the soil. 15The ratio of nitrogen to the total nitrogen content can be used to study soil nitrogen cycle processes, such as nitrogen fixation, mineralization, nitrification and denitrification.

[0094] (5) Establish multi-source heterogeneous data

[0095] The obtained soil physical protection component carbon (iPOM), chemical protection component carbon data (non-acid decomposable free silt component carbon (NH-dSilt), non-acid decomposable free clay component carbon (NH-dClay), non-acid decomposable silt component carbon (NH-μSilt), non-acid decomposable clay component carbon (NH-μClay), acid decomposable free silt component carbon (H-dSilt), acid decomposable free clay component carbon (H-dClay), acid decomposable silt component carbon (H-μSilt) carbon and acid decomposable clay component carbon (H-μClay)), biologically active component carbon data (coarse particle component carbon (cPOM) and free particle component carbon (fPOM)), isotope tracer data ( 13 C abundance, carbon from physical protected components, carbon from chemical protected components, and carbon from biologically active components (including the content and ratio of native soil carbon and exogenous carbon, the direction and intensity of the stimulating effect, and the amount and rate of mineralization) are organized into tables, recording information such as the collection location, succession stage, and measurement time of each sample. This data is then stored in a database, creating multi-source heterogeneous data sets that provide a foundation for subsequent analysis and modeling.

[0096] S2, integrate the multi-source heterogeneous data obtained in step S1 into the “physical-chemical-biological- 13 C abundance” four-dimensional feature space, and an adaptive weight distribution algorithm is used to dynamically adjust the contribution of each modal data.

[0097] The multi-source heterogeneous data obtained in the previous step are preprocessed to construct the "physical-chemical-biological- 13 The four-dimensional feature space of "C abundance" is divided into the following modes:

[0098] Physical modality: carbon content of iPOM component, chemical modality: NH-dSilt, NH-dClay, NH-μSilt, NH-μClay, H-dSilt, H-dClay, H-μSilt, H-μClay, biological activity modality: carbon content of cPOM component, carbon content of fPOM component, isotope modality: 13 C abundance, 15 N abundance, the content and proportion of physical protection component carbon, chemical protection component carbon, and biologically active component carbon derived from native soil carbon and exogenous carbon, respectively, the direction and intensity of the stimulation effect, and the mineralization amount and mineralization rate (time series data).

[0099] Data classification and standardization are as follows:

[0100] The physical protection component includes the organic carbon content of iPOM in the 53-250μm particle size, which is separated by dry sieving and sedimentation, measured by an elemental analyzer, and the data are normalized to g C / kg. The chemical protection component is obtained by extracting eight types of chemically bound carbon, such as NH-dSilt, NH-dClay, and H-μSilt, from the particle size <53μm, and then measuring them by hydrofluoric acid treatment and oxidation titration. The unit is unified as mgC / g. The biologically active component includes the carbon content of cPOM and fPOM in the particle size >250μm, which is separated by density flotation and measured by a TOC analyzer. The data are normalized to percentages. The isotope data include δ 13 C, δ 15 N abundance and the proportion of exogenous carbon are measured by isotope ratio mass spectrometry and baseline correction is performed (e.g., using the VPDB standard as a benchmark).

[0101] The four-dimensional tensor construction process is as follows:

[0102] First, define the following dimensions:

[0103] Physical mode (I1=3): >250μm, 53-250μm, <53μm particle size;

[0104] Chemical modality (I2=8): NH-dSilt, NH-dClay, NH-μSilt, NH-μClay, H-dSilt, H-dClay, H-μSilt, H-μClay;

[0105] Biological modality (I3=2): cPOM, fPOM;

[0106] Isotope mode (I4=4): δ¹³C, δ 15 N, proportion of exogenous carbon, and intensity of excitation effect.

[0107] Then, the data of each sampling point are tensor-filled according to the four-dimensional coordinates, for example: X(i1, i2, i3, i4) = i1th particle size - i2th chemical component - i3th biological component - i4th isotope parameter value.

[0108] Then, the feature tensor is decomposed. In this embodiment, two feature tensor decomposition methods are provided: the Tucker decomposition method or the PARAFAC decomposition method. These methods decompose multi-source heterogeneous data into a core tensor and a factor matrix. The data is then mapped into a four-dimensional feature space of "physical-chemical-biological-isotopic" through transformation of the factor matrix.

[0109] Take Tucker decomposition method as an example:

[0110] First, decompose the tensor into the core tensor G and the factor matrix U(i) :

[0111] ;

[0112] Core Tensor , represents the interaction strength between each mode, and the feature space after dimension compression. Factor matrix , indicating that the column vectors are orthogonal.

[0113] The specific steps are as follows:

[0114] First initialize and randomly generate ;

[0115] Then perform iterative optimization: fix other modes, update :

[0116] ;

[0117] Similarly, update in sequence ;

[0118] Set the convergence conditions, for example, in this example, set the objective function change rate to <0.1% or reach the maximum iteration of this tree (100 times).

[0119] Take PARAFAC (CP) decomposition as an example:

[0120] PARAFAC decomposes a tensor into a sum of rank 1 components. For a 4D tensor:

[0121] ;

[0122] in: is the weight coefficient; is the rank 1 vector of each mode; represents the outer product. Represents the set of real numbers, which is the set of all rational and irrational numbers. R represents a tensor of a specific structure.

[0123] The specific decomposition steps are as follows:

[0124] (1) Objective function construction

[0125] Minimize the residual sum of squares:

[0126] ;

[0127] (2) Alternating Least Squares (ALS) Optimization

[0128] Fixed Update all variables except :

[0129] ;

[0130] Update in sequence , , , the method is similar. Among them, X （1） Represents the first modal expansion matrix of the tensor X.

[0131] (3) Normalization and weight calculation

[0132] Normalize each vector: ;

[0133] Weight .

[0134] Each rank 1 component represents an independent influencing factor, such as the coupling of the aggregate wrapping mode (physical mode) and a specific chemical protection mechanism (chemical mode). Reflects the contribution strength of the factor to the overall data. Represents a tensor product.

[0135] Example: For a 4D tensor , take R=3, decomposition result:

[0136] (Physical mode: iPOM dominant);

[0137] (Chemical mode: H-clay dominant);

[0138] =12.5, indicating the contribution of this factor to the data variation.

[0139] When the adaptive weight allocation algorithm is used to dynamically adjust the contribution of each modal data, the weight is first adjusted dynamically based on the information entropy of the factor matrix:

[0140] Calculate the variance contribution of each mode:

[0141] ;

[0142] in, is the factor matrix U (i) The variance of the kth column;

[0143] Modal information entropy:

[0144] ;

[0145] Adjust weights based on succession stage:

[0146] ;

[0147] Among them, S i is the stage adjustment coefficient, which can be set as follows:

[0148] When the succession stage is the herbaceous stage, the biological activity mode S3 = 0.5, the physical mode S1 = 0.8, the chemical mode S2 = 0.6, and the isotope mode S4 = 0.7;

[0149] When the succession stage is the tree stage, the biological activity mode S3 = 0.8, the physical mode S1 = 0.5, the chemical mode S2 = 1.2, and the isotope mode S4 = 0.9;

[0150] Weighted fusion feature tensor

[0151] where X i is the normalized sub-tensor of each mode.

[0152] Project the original data into four-dimensional space:

[0153] F= ;

[0154] Output tensor F∈ , retaining more than 95% of the original information.

[0155] For the physical modal principal components, the first principal component usually corresponds to the synergistic effect between the <53 μm particle size and iPOM. For the chemical modal principal components, the first principal component can be interpreted as the adsorption competition between H-μClay and NH-dSilt.

[0156] The above method accurately characterizes the chemical protection mechanism by subdividing NH-silt / NH-clay and other subcategories, and combines information entropy with succession stage coefficients to achieve intelligent adaptation that focuses on biological activity in the herbaceous stage and on chemistry in the arbor stage.

[0157] S3, based on the spatial constraints of the boundary effect method and the statistical constraints of the least squares method, takes the four-dimensional feature tensor and weight distribution results output in step S2 as input, constructs a hybrid optimization model for predicting carbon saturation capacity, and outputs the maximum carbon sequestration capacity of microaggregates <53 μm and the carbon proportion of the components;

[0158] The specific implementation process is as follows:

[0159] Receive the four-dimensional feature tensor F∈ output from step S2 And adaptive weight distribution results , extract the physical modal eigenvector and chemical mode eigenvectors ;

[0160] Generate the particle size constraint matrix A according to the physical mode eigenvector, and its diagonal element a ii Calculated by the following formula:

[0161] ;

[0162] in, is the physical modal weight;

[0163] Based on chemical modal eigenvector and weight , calculate the theoretical saturation threshold of carbon content in each particle size :

[0164]

[0165] Among them, β k is the regression coefficient, C base To calibrate the reference value for the laboratory;

[0166] The sigmoid function is a nonlinear activation function, and its mathematical expression is usually:

[0167] ;

[0168] Its function is to map any real number z to the interval (0,1). In the model of the present invention, the sigmoid function is used to linearly combine the chemical features. Converted into a proportional coefficient to dynamically adjust the laboratory calibration reference value C base , and finally get the theoretical saturation threshold .

[0169] When the chemical features are combined As it approaches positive infinity, sigmoid(z)→1, the threshold ≈C base ;

[0170] When z approaches negative infinity, sigmoid(z)→0, the threshold →0 (in practice, negative values ​​are avoided through experimental design);

[0171] Through this function, the threshold In (0,C base ) range to achieve the inhibitory or enhancement effect of the chemical protection mechanism on carbon capacity.

[0172] is the chemical modal feature vector extracted by tensor decomposition (Tucker / PARAFAC) in step S2, representing the dimensionality reduction characteristics of the chemical protective components (such as the synergistic effects of eight types of components such as NH-dSilt and H-μClay).

[0173] is the transpose of the eigenvector and is used to compare with the regression coefficient β k Perform a dot product operation:

[0174] ;

[0175] Where R2 is the number of principal components of the chemical mode. The above steps quantify the contribution of chemical features to the carbon saturation threshold of a specific particle size (e.g., <53 μm) through linear combination.

[0176] C base To calibrate the reference value for the laboratory, design an experiment to determine it through the following steps:

[0177] Typical karst soil samples (e.g., climax community and shrub stage) were selected, and the carbon adsorption saturation values ​​of each particle size (>250μm, 53-250μm, and <53μm) were measured through constant temperature and humidity incubation experiments. Environmental conditions were controlled (temperature 25°C, humidity 60%, and no exogenous carbon input) to ensure that the measurement results reflected the intrinsic carbon holding capacity of the soil.

[0178] Take the average value of multiple experiments for the same particle size, for example, the C value of microaggregates <53 μm is base =18.2 g·C / kg; ensure the generalizability of the benchmark value (error <5%) through independent sample verification (e.g., cross-validation).

[0179] β k is the regression coefficient vector, which is determined by fitting historical data. The specific process is as follows:

[0180] Collect chemical characteristics data of multiple sets of karst soil samples and the measured saturation threshold of the corresponding particle size ; Samples must cover different succession stages (grassland, shrubs, trees) and landform types (peak cluster depression, canyon, etc.).

[0181] Use logistic regression or maximum likelihood estimation to optimize the following objective function:

[0182] ;

[0183] Introduce regularization (such as L2 regularization) to prevent overfitting. Select the optimal β through cross-validation k , the determination coefficient R2 is required to be ≥ 0.85; the prediction error is verified on an independent test set (such as < 10%).

[0184] The above steps are passed Dynamic Adjustment , so that the threshold value changes adaptively with the strength of the chemical protection mechanism; based on experimental calibration and regression fitting β k , ensuring the ecological rationality of the model;

[0185] Chemical feature vectors Covers data from multiple succession stages and supports the application of the model in the grassland → tree stage migration.

[0186] For example, in a karst peak cluster depression, the above method was used to calibrate the <53μm =18.2 g·C / kg, and fitting β k After that, the model predicts the threshold =16.5 g·C / kg (measured value 16.8 g·C / kg), with an error of only 1.8%.

[0187] Construct the following hybrid objective function:

[0188] ;

[0189] Where A is the particle size constraint matrix, =[ 1, 2, 3] T is the particle size carbon capacity prediction variable, b is the measured total organic carbon vector, is the penalty coefficient, is the theoretical saturation threshold of carbon content in each particle size, is the particle size classification index parameter, =1, 2, 3, corresponding to the particle size index parameters of >250μm, 53-250μm, and <53μm, respectively. For the Predictors of carbon capacity of aggregates;

[0190] Dynamic penalty coefficient , with physical weight Adaptive adjustment;

[0191] A hybrid strategy of particle swarm optimization (PSO) and genetic algorithm (GA) is used:

[0192] First, data initialization is performed to generate a 10×R1 group of initial solutions. In the computer system, the search range is limited to [0.5 ,1.2 ];

[0193] Perform fitness calculation and define the weighted fitness function:

[0194] ;

[0195] in, For fitness, is the physical mode weight, is the chemical modal weight;

[0196] Particle swarm in physical modal space f phy Gradient search along the direction; genetic algorithm for chemical modal parameter β k Perform arithmetic crossover and Gaussian mutation; terminate when the fitness change rate of the optimal solution is <0.05% for five consecutive generations;

[0197] Output <53μm microaggregates maximum carbon capacity And the proportion of carbon components:

[0198] ;

[0199] In the above methods, PSO searches along the principal component direction in the physical modal space, which improves the efficiency of particle size capacity prediction; GA enhances the adaptability of boundary constraints for chemical modal parameter variation; multi-objective optimization is achieved through weighted fitness function, and the weight distribution result { }Directly regulate the balance between statistical items and boundary items.

[0200] The output of the component carbon percentage provides a basis for adjusting the stability threshold in step S4 (e.g., lowering the critical value when the percentage is >50%).

[0201] S4, based on the maximum carbon capacity predicted in step S3, a mineralization culture experiment was designed, and the soils of different succession stages of the karst forest were mineralized using isotope tracing. 13 C / 15 The time-varying calibration curve of the N double standard sample analyzes the dynamic ratio of exogenous carbon to native carbon and the input and output pathways of soil carbon, and calculates the carbon stability disturbance threshold.

[0202] in, 13 C / 15 The method for constructing a time-varying calibration curve of N dual standard samples includes:

[0203] Based on the maximum carbon capacity of microaggregates <53 μm output in step S3, a mixed standard of 13C-glucose and 15N-ammonium sulfate was injected into the mineralization culture system every 12 hours, with each injection lasting 5 minutes. The concentrations of the mixed standard samples were 15% of the carbon content and 3% of the nitrogen content of the original soil, respectively.

[0204] Set the amount of exogenous carbon added to 0.2 -0.5 .For example, =20 g·C / kg, then add 4-10 g·C / kg 13 C-glucose. Under normal circumstances, the amount of nitrogen added is 3%-5% of the total nitrogen (TN) of the original soil. 15 N-ammonium sulfate to match the carbon-nitrogen stoichiometric ratio (e.g., C:N ≈ 10:1).

[0205] Gas samples were collected 1, 3, 6, and 12 hours after each injection and measured by isotope ratio mass spectrometry. 13 C and 15 N abundance; soil was collected on days 0, 7, 14, 21, and 28 of incubation, and different particle sizes (>250 μm, 53-250 μm, and <53 μm) were separated and the δ 13 C and δ 15 N value.

[0206] A carbon-nitrogen coupling kinetic model was established; 13 C abundance uses an exponential decay equation to characterize carbon turnover. 15 The logarithmic growth equation was used to characterize nitrogen fixation for N abundance;

[0207] The carbon turnover rate constant and nitrogen retention characteristic time were fitted simultaneously by the simulated annealing-least squares hybrid algorithm to complete the construction. 13 C / 15 Time-varying calibration curve of N dual standards.

[0208] The carbon-nitrogen coupling kinetic model includes the following coupling equations:

[0209] 13 C decay equation:

[0210] in, is the carbon turnover rate constant, which is solved by inversion using the simulated annealing algorithm;

[0211] 15N adsorption equation:

[0212] in, is the characteristic time of nitrogen fixation, The smaller the size, the faster the fixation rate (e.g. clay may have a high specific surface area). <sand soil). 、 is a soil type parameter (such as organic matter content and mineral composition). α is positively correlated with nitrogen adsorption capacity, and β represents the initial adsorption background value.

[0213] The carbon flow path analysis and excitation effect quantification process is carried out in the following steps:

[0214] First, the input and output paths need to be determined.

[0215] The rate of exogenous carbon input (R in ) indicates that based on f exo (t) Calculate the amount of exogenous carbon mineralization per unit time. exo (t) represents the cumulative contribution ratio of exogenous carbon to the total soil carbon mineralization at time t.

[0216] Native carbon output rate (R out ) represents the total mineralization minus the exogenous part.

[0217] The carbon flow balance ratio is calculated according to the following formula:

[0218] ;

[0219] If η>1.2, it is marked as a net carbon loss path; if η<0.8, it is marked as a net carbon fixation path.

[0220] The excitation effect is quantized as follows:

[0221] Calculate the cumulative priming effect (CPE):

[0222] ;

[0223] is the microbial diversity index (Shannon index based on 16S rRNA sequencing). The lower the diversity, the more significant the stimulation effect. The integral interval is from the initial time t0 to the critical time t c (Usually the end of the culture period, such as 28 days).

[0224] Finally, the carbon stability disturbance threshold is determined. First, a dynamic threshold needs to be set, and the carbon stability disturbance threshold T is adjusted according to the carbon proportion of the <53μm component (P<53). crit :

[0225] ;

[0226] Among them, P<53 represents the carbon proportion of microaggregates <53 μm in the soil (unit: %), which represents the physical protection strength of carbon. The stronger the physical protection ability of microaggregates (<53 μm) for carbon (the higher the P<53), the more stable the carbon pool is, and the carbon stability disturbance threshold T crit The lower the carbon loss rate, the lower the carbon loss rate. 2.5 is a coefficient determined by regression of historical data, reflecting the empirical relationship between carbon loss rate and microaggregate protection capacity.

[0227] For example, if P < 53 = 60%, Tcrit = 2.5 × (1 − 0.6) = 1.0 .

[0228] In this embodiment, the carbon pool is determined to be unstable when the following conditions are met simultaneously:

[0229] Excitation effect exceeds limit: CPE>Tcrit;

[0230] Isotope shift verification: δ13C shift exceeds the background value by ±2‰ (excluding instrument error interference).

[0231] To ensure the reliability of the threshold, the following verification is required:

[0232] For example, sensitivity analysis can be performed by perturbing the physical protection weights and chemical protection weights ±10%, with a threshold coefficient of variation (CV) of <8%. Field evidence can also be combined with simultaneous verification in peak-cluster depressions, canyons, and basins, requiring a spatial match of ≥85% between the predicted and measured carbon loss hotspots.

[0233] S5, construct a carbon cycle digital twin to simulate the carbon saturation deficit and spatial heterogeneity of soils at different succession stages of karst forests.

[0234] Digital twin construction includes:

[0235] A three-dimensional carbon pool interactive map was established, and the physical protection component carbon, chemical protection component carbon, and biologically active component carbon were mapped in three-dimensional space as vector fields according to the different succession stages of herbs, shrubs, and trees. An LSTM-Transformer hybrid neural network was embedded, and the input layer received the cumulative carbon input, carbon mass of components <53μm, and carbon content of components <53μm. The output layer predicted the carbon saturation deficit value, and combined with the existing carbon content and carbon retention rate, the carbon saturation capacity time limit was estimated.

[0236] The specific implementation process is as follows:

[0237] First, a three-dimensional carbon pool interactive map was constructed. The primary goal was to dynamically map the physical, chemical, and biologically active components of the soil carbon pool in three dimensions and quantify the correlation between carbon distribution and successional stage. Three-dimensional topographic data of the karst landscape (resolution ≤ 1m) was acquired using drone LiDAR or a high-precision terrain scanner. Sample plots were then divided according to successional stage (herbs, shrubs, and trees), with at least 50 sampling points deployed in each plot.

[0238] For physically protected carbon, microaggregates <53 μm were separated by density flotation, and their carbon content (unit: g·C / kg) was determined. For chemically protected carbon, mineral-bound carbon was extracted by acid hydrolysis (6 M HCl). For biologically activated carbon, microbial biomass carbon (MBC) was extracted by chloroform fumigation.

[0239] Then, using Python's PyVista or ParaView, a 3D spatial grid with a cell size of ≤10 m³ was constructed to implement 3D vector field mapping. The mapping method was to assign the physical, chemical, and biologically active carbon content of each sampling point to the X, Y, and Z components of the vector field, respectively. For example, if a point has a physical protective carbon content of 15 g / kg, a chemical protective carbon content of 8 g / kg, and a biologically active carbon content of 2 g / kg, its vector coordinates would be (15, 8, 2).

[0240] Time series data (such as annual updates) can be used to reflect the spatial migration patterns of carbon components during the succession process.

[0241] Then the LSTM-Transformer hybrid neural network design was completed, the main goal of which was to predict the carbon saturation deficit value (ΔC) and estimate the carbon saturation capacity time limit (T sat ).

[0242] The model architecture includes at least an input layer, a network structure, and an output layer, where the network structure includes an LSTM layer (long short-term memory network), a Transformer layer, and a fusion module.

[0243] For the input layer, feature dimensions include cumulative carbon input (total historical carbon input, unit: t·C / ha), carbon mass of fractions <53 μm (unit: g·C / kg), and carbon percentage of fractions <53 μm (unit: %). Time series processing is performed, and input data is arranged by time step (e.g., year), with a maximum lookback period of 20 years.

[0244] The LSTM layer (Long Short-Term Memory) in the network structure has a neuron count of 64 to capture temporal dependencies of carbon inputs (such as the lag effect of carbon sequestration). A multi-head attention mechanism (four heads) is used in the Transformer layer to focus on spatial heterogeneity (such as the impact of microaggregate distribution on carbon saturation). The fusion module in the network structure uses cross-attention to fuse the temporal features of the LSTM output with the spatial features of the Transformer.

[0245] For the output layer, the carbon saturation deficit value (ΔC) is defined as the difference between the predicted current carbon content and the saturation capacity (unit: g·C / kg). The carbon saturation capacity time limit (T sat ) is calculated by combining the carbon retention rate (k sat , unit: g·C / kg / yr), calculate the time required to reach saturation:

[0246] ;

[0247] Model training and validation are carried out in the following steps:

[0248] First, the dataset was divided into a training set (70%) covering representative plots of herbaceous, shrub, and tree succession stages. The validation set (15%) was used for hyperparameter tuning, and the test set (15%) was used to evaluate the model's generalization ability.

[0249] Using the Huber loss function, we can balance the mean square error (MSE) and the mean absolute error (MAE):

[0250] ;

[0251] In this embodiment, it is set = 1.0, to adapt to the noise interference of carbon content prediction. Finally, the AdamW optimizer (learning rate = 1e-4, weight decay = 1e-5) can be used.

[0252] Verification accuracy requirements are set as follows: carbon saturation deficit prediction error ≤ 10% (R² ≥ 0.85); carbon saturation time limit error ≤ 15% (RMSE ≤ 2 years). Sensitivity testing sets perturbation input parameters (e.g., ±5% of carbon <53μm), and model output fluctuation should be ≤ 8%.

[0253] The deployment and application of the digital twin begins with real-time data fusion. IoT sensors collect soil moisture, temperature, and microbial activity data in real time, dynamically updating the bioactive carbon content in the three-dimensional map. Model predictions are synchronized every six hours to generate a carbon saturation heat map (red: high deficit; green: low deficit).

[0254] Set an early warning threshold, for example, when ΔC > critical value (such as ΔC > 5g·C / kg) and T sat If the ΔC is less than 10 years, an alert is triggered for the carbon sequestration priority zone. If the ΔC is high and the proportion of microaggregates is low, the system automatically recommends applying biochar (to enhance physical protection). If the ΔC is low but the proportion of chemical protection is high, the system automatically recommends planting deep-rooted plants (to enhance mineral binding).

[0255] The basic principles, main features, and advantages of the present invention are shown and described above. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are merely illustrative of the principles of the present invention. Various changes and modifications may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and modifications are intended to fall within the scope of the present invention. The scope of protection claimed in the present invention is defined by the appended claims and their equivalents.

Claims

1. A multimodal data fusion method for extracting carbon cycle characteristics of karst forest soil, characterized by: The following steps are involved: S1: Obtain carbon data of soil physical protection components, chemical protection components, biologically active components, and isotope tracer data at different succession stages of karst forests to establish multi-source heterogeneous data; S2, integrates the multi-source heterogeneous data obtained in step S1 into the "physical-chemical-biological-isotopic" four-dimensional feature space through feature tensor decomposition, and uses an adaptive weight allocation algorithm to dynamically adjust the contribution of each modality data; S3, based on the spatial constraints of the boundary effect method and the statistical constraints of the least squares method, takes the four-dimensional feature tensor and weight distribution results output in step S2 as input, constructs a hybrid optimization model for predicting carbon saturation capacity, and outputs the maximum carbon capacity of microaggregates <53 μm and the carbon proportion of the components; S4, based on the maximum carbon capacity predicted in step S3, a mineralization culture experiment was designed, and the soils of different succession stages of the karst forest were mineralized using isotope tracing. 13 C / 15 The time-varying calibration curve of the N double standard sample analyzes the dynamic ratio of exogenous carbon to native carbon and the input and output pathways of soil carbon, and calculates the carbon stability disturbance threshold; S5, construct a carbon cycle digital twin, embed the carbon stability disturbance threshold and input and output paths of step S4, and simulate the soil carbon saturation deficit value and spatial heterogeneity at different succession stages of karst forest.

2. The multimodal data fusion karst forest soil carbon cycle feature extraction method according to claim 1 is characterized in that: In step S1, the soil physical protection component carbon data is the 53-250 μm carbon content, including the physical protection component carbon iPOM; The chemical protection component carbon data is the carbon content of <53μm, including non-acid-decomposable free powder component carbon NH-dSilt, non-acid-decomposable free clay component carbon NH-dClay, non-acid-decomposable closed accumulation powder component carbon NH-μSilt, non-acid-decomposable closed accumulation clay component carbon NH-μClay, acid-decomposable free powder component carbon H-dSilt, acid-decomposable free clay component carbon H-dClay, acid-decomposable closed accumulation powder component carbon H-μSilt and acid-decomposable closed accumulation clay component carbon H-μClay; The bioactive component carbon data refers to the carbon content >250 μm, including coarse particle component carbon cPOM and free particle component carbon fPOM; The isotope tracing data includes 13 C abundance, 15 N abundance, the content and proportion of physical protection component carbon, chemical protection component carbon, and biologically active component carbon derived from original soil carbon and exogenous carbon, the direction and intensity of the excitation effect, and the mineralization amount and mineralization rate.

3. The multimodal data fusion karst forest soil carbon cycle feature extraction method according to claim 1 is characterized in that: In step S2, the feature tensor decomposition uses Tucker decomposition or PARAFAC decomposition to decompose the multi-source heterogeneous data into core tensors and factor matrices. The data is mapped to the "physical-chemical-biological-isotopic" four-dimensional feature space through transformation of the factor matrix.

4. The multimodal data fusion karst forest soil carbon cycle feature extraction method according to claim 1 is characterized in that: The carbon saturation capacity prediction hybrid optimization model in step S3 is constructed according to the following steps: 1) receiving the four-dimensional feature tensor and the adaptive weight distribution result outputted from step S2, and extracting the physical modal feature vector and the chemical modal feature vector; 2) Generate a particle size constraint matrix based on the physical modal eigenvectors and calculate the theoretical saturation threshold of carbon content in each particle size based on the chemical modal eigenvectors and weights; 3) Construct the following objective function: ； Where A is the particle size constraint matrix, =[ 1, 2, 3] T is the particle size carbon capacity prediction variable, b is the measured total organic carbon vector, is the penalty coefficient, is the theoretical saturation threshold of carbon content in each particle size, is the particle size classification index parameter, =1, 2, 3, corresponding to the particle size index parameters of >250μm, 53-250μm, and <53μm, respectively. For the Predictors of carbon capacity of aggregates; 4) A hybrid strategy of particle swarm optimization and genetic algorithm was used to output the maximum carbon capacity and component carbon proportion of microaggregates <53 μm.

5. The method for extracting carbon cycle characteristics of karst forest soil using multimodal data fusion according to claim 1 is characterized in that: As described in step S4 13 C / 15 The method for constructing a time-varying calibration curve of N dual-standard samples includes: Inject into the mineralization culture system in a 12-hour cycle 13 C-glucose and 15 N-ammonium sulfate mixed standard sample, each injection time is 5 minutes, its concentration is 15% of the original soil carbon content and 3% of the nitrogen content; Gas samples were collected 1, 3, 6, and 12 hours after each injection and measured by isotope ratio mass spectrometry. 13 C and 15 N abundance; A carbon-nitrogen coupling kinetic model was established; 13 C abundance uses an exponential decay equation to characterize carbon turnover. 15 The logarithmic growth equation was used to characterize nitrogen fixation for N abundance; The carbon turnover rate constant and nitrogen retention characteristic time were simultaneously fitted by the simulated annealing-least squares hybrid algorithm to complete the construction. 13 C / 15 Time-varying calibration curve of N dual standards.

6. The multimodal data fusion karst forest soil carbon cycle feature extraction method according to claim 5, characterized in that: The carbon-nitrogen coupling kinetic model includes the following coupling equations: 13 C decay equation: ; in, is the carbon turnover rate constant, which is solved by inversion using the simulated annealing algorithm; 15 N adsorption equation: ; in, is the characteristic time of nitrogen fixation, 、 is the soil type parameter.

7. The multimodal data fusion karst forest soil carbon cycle feature extraction method according to claim 6 is characterized in that: In step S4, based on the maximum carbon capacity of microaggregates <53 μm outputted in step S3, exogenous 13 C-labeled glucose and 15 N-labeled ammonium sulfate for mineralization culture; During the mineralization incubation process, gas and soil samples were collected every 12 hours and measured simultaneously by isotope ratio mass spectrometry. 13 C and 15 N dynamic value; establish 13 C decay and 15 The double standard time-varying calibration curve coupled with N adsorption separates the contribution of exogenous carbon and the interference of nitrogen transformation; Based on the calibration curve, the exogenous carbon input rate and the native carbon output rate are calculated, and the carbon flow balance ratio is defined. The carbon net loss path and the carbon net fixation path are determined by the carbon flow balance ratio. Based on the dynamic monitoring data of mineralization rate, the difference between the mineralization rate of exogenous carbon and native carbon was integrated over time to generate the cumulative stimulation effect, which was then normalized by combining it with the microbial diversity index. The carbon stability perturbation threshold is determined based on the proportion of carbon in the component <53 μm, and whether the carbon pool has entered an unstable state is determined by the carbon stability perturbation threshold and the normalized cumulative excitation effect; The microbial diversity index includes the Shannon index based on 16S rRNA sequencing.

8. The multimodal data fusion karst forest soil carbon cycle feature extraction method according to claim 1, characterized in that: The digital twin construction in step S5 includes: Establish a three-dimensional carbon pool interactive map, and map the physical protection component carbon, chemical protection component carbon, and biological active component carbon in three-dimensional space using vector fields according to the different succession stages of herbaceous, shrubby, and tree soils. Embedded in the LSTM-Transformer hybrid neural network, the input layer receives the cumulative carbon input, carbon mass of components <53μm and carbon content of components <53μm. The output layer predicts the carbon saturation deficit value and calculates the carbon saturation capacity time limit based on the existing carbon content and carbon retention rate.

9. The method for extracting carbon cycle characteristics of karst forest soil using multimodal data fusion according to claim 4, characterized in that: The hybrid strategy of particle swarm optimization and genetic algorithm to output the maximum carbon capacity and component carbon proportion of microaggregates <53 μm includes the following steps: Generate multiple sets of initial solutions, limit the search range in the computer system, and define the following weighted fitness function: ； in, For fitness, is the physical mode weight, is the chemical modal weight; By making the particle swarm perform gradient search in the physical modal space along the direction of the physical modal eigenvector, the genetic algorithm performs arithmetic crossover and Gaussian mutation on the chemical modal parameters to complete collaborative evolution; when the fitness change rate of the optimal solution is <0.05% for five consecutive generations, the algorithm is terminated and the result is output.

Images