Karst forest soil carbon cycle feature extraction method based on multi-modal data fusion

Abstract

The invention discloses a multi-modal data fused karst forest soil carbon cycle feature extraction method. The method comprises the following steps: acquiring soil physical and chemical protection component carbon data, bioactive component carbon data, isotopic tracing data and the like of different succession stages of a karst forest to establish multi-source heterogeneous data; integrating the obtained multi-source heterogeneous data into a'physical-chemical-biological-isotope 'four-dimensional feature space through feature tensor decomposition; constructing a carbon saturation capacity prediction hybrid optimization model, and outputting lt; the maximum carbon sequestration capacity and component carbon proportion of 53 [mu] m microaggregates; the method comprises the following steps: performing mineralization culture on soil in different succession stages of the karst forest, analyzing a dynamic proportion of exogenous carbon to original soil carbon and an input / output path of soil carbon in combination with a 13C / 15N double-standard-sample time-varying calibration curve, and calculating a carbon stability disturbance threshold. According to the method, collaborative analysis of multi-source data is realized, and the temporal-spatial resolution and the model precision of carbon cycle feature extraction are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the fields of environmental science and agricultural ecological technology, and specifically relates to a method for extracting the characteristics of the karst forest soil carbon cycle based on multi-modal data fusion, which is applicable to the assessment of the carbon sink potential of the karst ecosystem, the analysis of the carbon cycle path, and the formulation of sustainable carbon management strategies. Background Art

[0002] Soil carbon is a core component of the carbon cycle in terrestrial ecosystems, and its dynamic changes directly affect global climate change and soil quality. As an important ecosystem in southern China, the soil of karst forests is characterized by high heterogeneity, low carbon storage, and easy loss. Traditional research mainly relies on single data sources (such as physical component analysis or static models), making it difficult to comprehensively reveal the complex mechanisms of the carbon cycle.

[0003] Existing methods are mostly based on single-dimensional data (such as soil physical classification or chemical protection component analysis), ignoring the synergistic effects of multi-modal data (such as hyperspectral remote sensing, microbial communities, bioactive component carbon, isotope tracing). For example, estimating carbon storage only through soil bulk density and total organic carbon content cannot analyze the carbon sequestration differences of different aggregates (such as >250μm, <53μm) and the microbial regulation effect, resulting in insufficient accuracy in quantifying the carbon cycle path.

[0004] Traditional carbon cycle models (such as the Century model) are mostly static empirical formulas, lacking the ability to respond in real time to dynamic environmental factors (such as temperature, humidity, pH value) and spatial heterogeneity. For example, existing models are difficult to capture the boundary effect of carbon saturation capacity caused by terrain undulation in karst areas, and do not combine the boundary effect method with the least square method for hybrid optimization, restricting the accuracy of carbon sequestration potential prediction.

[0005] Although 13 C isotope labeling technology has been used to analyze carbon sources, existing methods mostly rely on a single isotope (such as 13 C) and lack dynamic calibration. For example, traditional mineralization experiments only analyze the ratio of exogenous carbon to native soil carbon, without combining 15 N double standard samples to construct a time-varying calibration curve, resulting in large quantization calculation errors in the priming effect (such as positive priming, negative priming), and unable to accurately quantify the carbon stability perturbation threshold.

[0006] Existing research has limited understanding of the regulatory role of microbial communities in the carbon cycle. Most are based on simple statistical analysis of community abundance and do not construct a dynamic regulatory network model. For example, high-throughput sequencing data combined with graph convolutional algorithms are not used to capture microbial synergistic effects, resulting in the inability to reveal the spatio-temporal dynamic laws of microbial-carbon interactions, restricting the development of intelligent monitoring and prediction of the carbon cycle.

[0007] In summary, due to the single data dimension, insufficient model dynamics, low precision of isotope analysis, and weak analysis of microbial mechanisms in the existing technology, it is difficult to meet the needs of refined research on the carbon cycle in karst forests. The present invention breaks through the above technical bottlenecks through multi-modal data fusion and dynamic modeling. Summary of the Invention

[0008] The purpose of the present invention is to solve the problems in the research on soil carbon cycle during the restoration of karst forests. The present invention integrates multi-modal data (soil physically protected fraction carbon data, chemically protected fraction carbon data, biologically active fraction carbon data, and isotope tracer data), constructs a four-dimensional feature space and a dynamic model, and realizes the accurate quantification of soil carbon cycle paths, saturation capacities, and priming effects. Combining carbon-nitrogen coupling kinetics analysis and digital twin simulation provides a scientific basis for the assessment of carbon sink potential, early warning of carbon instability risks, and formulation of sustainable carbon management strategies in karst ecosystems.

[0009] The present invention is implemented as follows: A method for extracting karst forest soil carbon cycle characteristics by multi-modal data fusion includes the following steps: S1. Obtain soil physically protected fraction carbon data, chemically protected fraction carbon data, biologically active fraction carbon data, and isotope tracer data at different succession stages of karst forests, and establish multi-source heterogeneous data; S2. Integrate the multi-source heterogeneous data obtained in step S1 into a "physical-chemical-biological-isotope" four-dimensional feature space through feature tensor decomposition, and dynamically adjust the contribution degrees of each modal data using an adaptive weight allocation algorithm; S3. Based on the spatial constraint of the boundary effect method and the statistical constraint of the least squares method, using the four-dimensional feature tensor and the weight allocation result output in step S2 as inputs, construct a carbon saturation capacity prediction hybrid optimization model, and output the maximum carbon sequestration capacity of <53μm micro-aggregates and the proportion of component carbon; S4. Design a mineralization incubation experiment based on the maximum carbon capacity predicted in step S3, conduct mineralization incubation on soils at different succession stages of karst forests using isotope tracer, and 13 C / 15 Analyze the dynamic ratio of exogenous carbon to native carbon and the input-output path of soil carbon by the time-varying calibration curve of the C / N double standard sample, and calculate the carbon stability perturbation threshold.

[0010] In the above solution, in step S2, the multi-source heterogeneous data (physical, chemical, biological, isotope) obtained in the previous step is mapped to a unified structured space to solve data heterogeneity (differences in units, dimensions, and formats), facilitating subsequent model processing. The contribution degrees of each modality (physical / chemical / biological / isotope) are dynamically adjusted through an adaptive algorithm. For example, the weight of the chemically protected fraction is higher at the climax community stage, and the weight of the physically protected fraction is higher at the shrub stage, ensuring that the model focuses on key influencing factors.

[0011] Among them, the four-dimensional feature tensor provides high-information-density input after dimensionality reduction, avoiding redundancy of the original data; and it can guide the model optimization in subsequent step S3 through weights, because the weight assignment result directly affects the penalty coefficients of each modality in the optimized model. For example, a high-weight modality occupies a larger proportion in the objective function, and the constrained model pays more attention to its contribution.

[0012] In addition, in the experimental design stage of subsequent step S4, the weight assignment result guides the addition ratio of isotope standard samples (for example, a high-weight isotope modality requires a higher sampling frequency). And it is convenient for dynamic threshold adjustment. For example, when the weight of the physical modality is high, the stability threshold is more stringent (because the physical protection effect is strong).

[0013] The reasons why step S4 needs to be based on the maximum carbon capacity of S3 mainly consider the following aspects: Firstly, based on ecological rationality, carbon capacity overload can be avoided, because exceeding the maximum capacity will lead to carbon saturation, causing unnatural loss and distorted experimental results; considering that the experiment simulates the real scenario, 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 that of the bare rock period.

[0014] Secondly, it is beneficial to parameter standardization. Unifying the benchmark with the maximum capacity ensures the comparability of experimental results at different succession stages (for example, adding 5% in the shrub stage vs adding 5% in the climax community represents different absolute amounts); and it is also beneficial to the association of dynamic thresholds. The excitation effect intensity needs to be adjusted according to the capacity ratio. For example, when the addition amount is 0.3x3, the threshold is reduced by 20%.

[0015] Finally, it is based on technical necessity. When calibrating the curve accuracy, the analysis of isotope signals (such as δ¹³C) needs to have a linear response within the carbon capacity range. Beyond the range may lead to non-linear errors; and it also takes into account 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.

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

[0017] The reasons why this step is a non-essential step in the present invention are as follows: If the research only needs to evaluate the carbon capacity at a certain succession stage (such as the shrub stage) and does not require cross-stage dynamic simulation, the complex modeling of the digital twin may be redundant. Moreover, the construction of the digital twin is also limited by data and resources because building a high-precision digital twin requires multi-dimensional real-time data (such as high-frequency isotope monitoring, UAV remote sensing). If the data acquisition cost is too high or the technical conditions are insufficient, it can be simplified to a static model.

[0018] If only the carbon saturation deficit value needs to be predicted, traditional empirical models (such as the Century model) or regression analysis can already meet the requirements, and there is no need to introduce digital twin technology. If the management goal is limited to a specific area (such as a single peak cluster depression), strategies can be directly formulated through the outputs of steps S3 - S4 (maximum capacity, threshold), without global simulation.

[0019] On the other hand, the construction of the digital twin requires a large amount of computing resources (such as GPU clusters, cloud platforms), and continuous calibration is needed (such as updating microbial community data annually), which may not be economical for small and medium-sized research teams. When the carbon management strategy approaches the optimal state, the improvement brought by the twin simulation may be limited and difficult to cover its development cost.

[0020] However, the addition of this step can also bring very positive effects. First, the karst forest soil carbon cycle involves multi-modal interactions of physics, chemistry, biology, and isotopes. The digital twin can integrate all the data and models of steps S1 - S4 to achieve cross-scale (microscopic - macroscopic) dynamic simulation of carbon flow, revealing non-linear and hysteresis effects. The micro-topographic differences in karst landforms (such as peak cluster depressions, karst caves) significantly affect carbon distribution. The digital twin quantifies the spatial differentiation law of carbon saturation deficit values through three-dimensional space mapping (such as Kriging interpolation, point cloud rendering) to guide precise management.

[0021] Second, due to the significant differences in the carbon cycle mechanisms among the grassland → shrub → tree stages, the digital twin can dynamically update model parameters (such as microbial diversity, aggregate structure) to adapt to the changes in carbon sink functions during the succession process. It can provide a simulation platform for policy designs such as carbon trading and ecological compensation. For example, it can evaluate the economic cost and ecological benefits of sequestering 1 ton of carbon.

[0022] In step S5, by embedding the carbon stability perturbation threshold of step S4, it can real-time simulate the carbon instability risk under different fertilization or land use scenarios and generate a risk heat map (such as the red warning area concentrated in the foot of the slope). Based on the input-output path (carbon loss / fixation dominant type), it can simulate the carbon sequestration effect of adding biochar or adjusting vegetation cover, providing a quantitative decision-making basis for ecological restoration.

[0023] In most cases, in step S1, the carbon data of the soil physical protection component is the carbon content of 53 - 250 μm, including the physical protection component carbon iPOM; The carbon data of the chemical protection component is the carbon content of <53μm, including the carbon of non-acid-dissolved free silt component NH-dSilt, the carbon of non-acid-dissolved free clay component NH-dClay, the carbon of non-acid-dissolved occluded silt component NH-μSilt, the carbon of non-acid-dissolved occluded clay component NH-μClay, the carbon of acid-dissolved free silt component H-dSilt, the carbon of acid-dissolved free clay component H-dClay, the carbon of acid-dissolved occluded silt component H-μSilt and the carbon of acid-dissolved occluded clay component H-μClay; The carbon data of the soil biological activity component is the carbon content of >250μm, including the carbon of coarse particle component cPOM and the carbon of free particulate component fPOM; The isotope tracer data includes 13 the abundance of 15 C, the abundance of

[0024] N, the content and proportion of the carbon of the physical protection component, the carbon of the chemical protection component, and the carbon of the biological activity component derived from the original soil carbon and exogenous carbon, the direction and intensity of the priming effect, the mineralization amount and the mineralization rate.

[0025] As a further preferred solution, in step S2, the characteristic tensor decomposition adopts Tucker decomposition or PARAFAC decomposition to decompose the multi-source heterogeneous data into a core tensor and factor matrices, and maps the data to a four-dimensional characteristic space of "physical-chemical-biological-isotope" through the transformation of the factor matrices. 1) Receive the four-dimensional characteristic tensor and the adaptive weight allocation result output in step S2, and extract the physical modal feature vector and the chemical modal feature vector; 2) Generate a particle size constraint matrix according to the physical modal feature vector, and calculate the theoretical saturation threshold of the carbon content of each particle size level based on the chemical modal feature vector and the weight; 3) Construct the following objective function: ; where A is the particle size constraint matrix, = 1, 2, 3] T is the prediction variable of the carbon capacity of the particle size level, b is the vector of the measured total organic carbon, is the penalty coefficient, is the theoretical saturation threshold of the carbon content of each particle size level, is the index parameter of the particle size classification, =1, 2, 3, corresponding to the index parameters of the particle size levels of >250μm, 53-250μm, and <53μm respectively, is the prediction variable of the carbon capacity of the th level of aggregates; 4) Output the maximum carbon capacity of <53μm micro-aggregates and the proportion of component carbon using a hybrid strategy of particle swarm optimization and genetic algorithm.

[0026] As a further preferred solution, in step S4, the 13 C / 15 Construction method of time-varying calibration curve for C / N double standard samples includes: Inject 13 C-glucose and 15 N-ammonium sulfate mixed standard samples into the mineralization culture system every 12 hours, with each injection time being 5 minutes and their concentrations being 15% of the original soil carbon content and 3% of the nitrogen content respectively; Collect gas samples at 1, 3, 6, and 12 hours after each injection, and measure 13 C and 15 N abundances by isotope ratio mass spectrometer; Establish a carbon-nitrogen coupling kinetic model; among them, 13 C abundance uses an exponential decay equation to characterize carbon turnover, 15 N abundance uses a logarithmic growth equation to characterize nitrogen fixation; Synchronously fit the carbon turnover rate constant and the nitrogen fixation characteristic time through a simulated annealing-least squares hybrid algorithm to complete the construction of 13 C / 15 Time-varying calibration curve for C / N double standard samples.

[0027] As a further preferred solution, the carbon-nitrogen coupling kinetic model includes the following coupling equations: 13 C decay equation: ; Among them, is the carbon turnover rate constant, which is inversely solved by the simulated annealing algorithm; 15 N adsorption equation: ; Among them, is the nitrogen fixation characteristic time, , are soil type parameters.

[0028] As a further preferred solution, in step S4, based on the maximum carbon capacity of <53μm micro-aggregates output in step S3, add exogenous 13 C-labeled glucose and 15 N-labeled ammonium sulfate to soil samples at different succession stages for mineralization culture; During the mineralization culture process, collect gas and soil samples every 12 hours, and synchronously measure 13 C and 15N dynamic value; establishment 13 C attenuation and 15 N adsorption-coupled double standard sample time-varying calibration curve to separate the contribution degree of exogenous carbon and nitrogen transformation interference; Calculate the exogenous carbon input rate and the native carbon output rate based on the calibration curve, define the carbon flow balance ratio, and determine the carbon net loss path and carbon net fixation path through the carbon flow balance ratio; Based on the dynamic monitoring data of the mineralization rate, perform time integration on the difference between the mineralization rates of exogenous carbon and native carbon to generate the cumulative priming effect, and normalize the cumulative priming effect in combination with the microbial diversity index; Determine the carbon stability perturbation threshold based on the carbon proportion of the <53μm fraction, and determine whether the carbon pool enters an unstable state through the carbon stability perturbation threshold and the normalized cumulative priming effect; Among them, the microbial diversity index includes the Shannon index based on 16S rRNA sequencing.

[0029] As a further preferred solution, the digital twin construction in step S5 includes: Establish a three-dimensional carbon pool interaction map, and map the carbon in the physical protection component, chemical protection component, and biological activity component in the three-dimensional space according to the soil in different succession stages of herbs, shrubs, and trees; Embed the LSTM-Transformer hybrid neural network. The input layer receives the cumulative carbon input, the carbon mass of the <53μm fraction, and the carbon content data of the <53μm fraction. The output layer predicts the carbon saturation deficit value, and combines the existing carbon content and carbon sequestration rate to calculate the carbon saturation capacity time limit.

[0030] As a further preferred solution, the hybrid strategy of particle swarm optimization and genetic algorithm is used to output the maximum carbon capacity of the <53μm microaggregate and the carbon proportion of the components, including the following steps: Generate multiple groups of initial solutions, limit the search range in the computer system, and define the following weighted fitness function: ; Among them, is the fitness, is the physical modal weight, is the chemical modal weight; Through gradient search of the particle swarm along the physical modal eigenvector direction in the physical modal space, the genetic algorithm performs arithmetic crossover and Gaussian mutation on the chemical modal parameters to complete coevolution; when the change rate of the fitness of the optimal solution is less than 0.05% for 5 consecutive generations, terminate and output the result.

[0031] Compared with the current situation in the prior art where the analysis mainly relies on a single data dimension (such as only physical particle size or static total carbon content) and it is difficult to analyze the multi-factor coupling mechanism of the carbon cycle, the present invention realizes the collaborative analysis of multi-source data, significantly improving the spatio-temporal resolution and model accuracy of carbon cycle feature extraction.

[0032] First of all, the present invention pioneered the "physical-chemical-biological-isotope" four-dimensional feature space fusion technology. By using Tucker / PARAFAC tensor decomposition to integrate particle size components (such as iPOM, H-μClay), chemically protected forms (NH-dSilt, etc.), bioactive carbon (cPOM / fPOM) and dual-standard isotope data, and developing an adaptive weight algorithm to dynamically allocate the contribution degree of each modality (such as the chemical weight is increased to 0.6 at the climax community stage). This collaborative modeling of multi-source heterogeneous data has greatly improved the carbon sequestration prediction accuracy compared with the traditional single-modal method, solving the problem of model distortion caused by the high heterogeneity of karst soil.

[0033] Secondly, most of the existing carbon capacity models adopt fixed empirical formulas (such as the Century model), ignoring the coordination of boundary effects and statistical constraints. The hybrid optimization model constructed by the present invention innovatively integrates the three-dimensional constraint surface of the boundary effect method and the least squares statistical fitting, and realizes the adaptive adjustment of the constraint strength driven by the physical weight ( ( =0.2 + 0.6*(1 - ))). For example, when the physical weight at the shrub stage is 0.3, =0.62 to strengthen the boundary constraint and prevent overloading prediction; while at the climax community =0.5, =0.5 focuses on statistical fitting. =0.5 focuses on statistical fitting.

[0034] On the other hand, the traditional priming effect analysis is limited to qualitative description. The present invention realizes the precise analysis of the dynamic ratio of exogenous carbon to native carbon during the mineralization process through the time-varying calibration curve of 13C / 15N dual-standard samples, combined with the weighted integration of the microbial diversity index (Shannon index). Through the carbon-nitrogen coupling kinetic model and the quantization algorithm of the priming effect, the traditionally qualitatively described priming effect is transformed into a computable integral quantity (such as the cumulative priming effect), and the cross-interference is eliminated by combining the isotope offset. This method breaks through the limitations of single-isotope analysis, reducing the calculation error of the carbon stability perturbation threshold to within ±5%, significantly improving the accuracy of carbon cycle mechanism research. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 The working flow chart in an embodiment of the present invention is shown. DETAILED DESCRIPTION OF THE INVENTION

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

[0037] In the following description, certain specific details are set forth for the purpose of explaining various disclosed embodiments 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 instances, well-known technologies associated with the present application may not be shown or described in detail so as to avoid unnecessarily obscuring the description of the embodiments.

[0038] References to "an embodiment" or "one embodiment" throughout the specification mean 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 an embodiment" or "in one embodiment" throughout the specification need not all refer to the same embodiment. Additionally, the particular features, structures, or characteristics may be combined in any manner in one or more embodiments.

[0039] As Figure 1 shown, a method for extracting the characteristics of the karst forest soil carbon cycle with multi-modal data fusion includes the following steps: S1. Obtain the carbon data of the physical protection component, chemical protection component, biological activity component and isotope tracer data of the soil in different succession stages of the karst forest, and establish multi-source heterogeneous data.

[0040] During specific implementation, sampling areas are respectively set in different succession stages of the karst forest (such as the herb stage, shrub stage, tree stage, etc.). A plurality of sampling points are set in each area according to a certain grid layout to ensure that the samples can represent the soil characteristics of this succession stage. For each sampling point, soil samples with a depth of 0-20 cm are collected, and the samples of multiple sampling points in the same area are mixed into a comprehensive sample to reduce sampling errors. 3-5 comprehensive samples are collected for each succession stage, impurities such as visible plant residues, roots and stones are removed, and they are taken back to the laboratory to be air-dried, and then gently ground with a grinding rod and passed through a 2 mm sieve for standby.

[0041] (1) Acquisition of carbon data of physical protection components The sieved soil is subjected to particle grading by the wet sieving method to obtain soil particles with three particle sizes of >250 μm, 53-250 μm, and <53 μm respectively. The specific operation is to put a certain amount of soil into sieves with different pore sizes and oscillate in water for a certain time to separate particles with different particle sizes.

[0042] For each particle size fraction of soil particles, the organic carbon content was determined by the potassium dichromate volumetric method - external heating method. The specific steps were as follows: Weigh a certain amount of soil sample and put it into a test tube, add a certain amount of potassium dichromate solution and concentrated sulfuric acid, heat in an oil bath for a certain period of time, cool, and then titrate the remaining potassium dichromate with a ferrous sulfate standard solution. Calculate the soil organic carbon content based on the titration results.

[0043] The iPOM organic carbon, namely physically protected organic matter, is the component carbon existing inside micro-aggregates (particle size ranging from 53 - 250 μm). Since it is wrapped inside the micro-aggregates and physically protected, its decomposition rate is relatively slow, and it has relatively high stability in the soil carbon pool, playing an important role in the long-term storage of soil carbon.

[0044] In the actual implementation process, the density flotation method was used to separate the iPOM component carbon. Put the soil sample into a solution with a certain density (such as NaI solution), after centrifugation, the iPOM component carbon will float on the surface of the solution. After collection and washing, use the above potassium dichromate volumetric method - external heating method to determine its carbon content.

[0045] (2) Obtaining data of chemically protected component carbon Chemically protected component carbon includes chemically protected carbon and biochemically protected carbon. Among them, chemically protected carbon was separated into H-silt carbon and H-clay carbon by chemical extraction method. Mix the soil sample with a sodium hydroxide solution of a certain concentration, oscillate for a certain period of time, and then centrifuge to separate the supernatant and the precipitate. The supernatant mainly contains H-silt carbon and H-clay carbon. Use the potassium dichromate volumetric method - external heating method to determine the carbon content in the supernatant, which is the total amount of chemically protected carbon. In order to distinguish H-silt carbon and H-clay carbon, further fractionation methods can be used, such as filtering through filter membranes with different pore sizes and measuring the carbon content of different particle size parts respectively.

[0046] H-silt carbon is the carbon combined with silt particles and protected by chemical action. The particle size of silt is between sand particles and clay particles, and this combination makes the organic carbon have a certain stability in the soil and is not easily decomposed by microorganisms. H-clay carbon is the chemically protected carbon combined with clay particles. Clay has a large specific surface area and strong adsorption ability, and can combine tightly with organic carbon through chemical bonds and other interactions, further enhancing the stability of organic carbon.

[0047] Biochemical protected carbon is separated into NH-silt carbon and NH-clay carbon by enzymatic hydrolysis combined with chemical extraction. First, the soil sample is mixed with a specific enzyme solution (such as cellulase, protease, etc.), reacted for a certain time under certain temperature and pH conditions, then extracted with sodium hydroxide solution, and the supernatant is separated by centrifugation. The potassium dichromate volumetric method - external heating method is used to determine the carbon content in the supernatant, which is the total amount of biochemical protected carbon. Similarly, NH-silt carbon and NH-clay carbon can be distinguished by further fractionation methods.

[0048] NH-silt carbon and NH-clay carbon are respectively carbon combined with silt and clay particles and protected by biochemical processes. Biochemical protection may involve the interaction between microbial metabolites, soil particles, and organic carbon, forming more complex structures that make organic carbon more difficult to decompose, which is of great significance for the long-term storage of soil carbon and the maintenance of soil fertility.

[0049] (3)Obtaining data on bioactive component carbon Bioactive component carbon includes cPOM and fPOM. Among them, cPOM component carbon usually refers to coarse particulate organic carbon, which has relatively large particles, mainly derived from fresh plant residues, etc., has high biological activity, is easily decomposed and utilized by microorganisms, and has a relatively fast turnover in the soil carbon cycle. fPOM component carbon is free organic matter in aggregates, that is, carbon in aggregates larger than 250 μm, and also has high activity, which can quickly provide energy and nutrients for soil microorganisms.

[0050] Physical separation combined with chemical analysis is used to separate and determine bioactive component carbon. That is, the soil sample is sieved through a 53 μm sieve, and the part on the sieve is cPOM component carbon and fPOM component carbon. The potassium dichromate volumetric method - external heating method is used to determine the carbon content of the separated cPOM component carbon and fPOM component carbon.

[0051] To distinguish cPOM component carbon and fPOM component carbon, they can be further separated according to their physical properties (such as color, texture, etc.), and then their carbon contents are measured separately.

[0052] (4)Obtaining isotope tracer data Respectively place 13 C / 15 N double-labeled herb, shrub, and tree soil 50 g soil samples in 50 mL beakers, and inject 13 C-glucose and 15 N-ammonium sulfate mixed standard samples into the mineralization culture system at a cycle of 12 hours. The injection time for each time is 5 minutes, and their concentrations are 15% of the original soil organic carbon content and 3% of the total nitrogen content respectively; gas samples are collected at 1, 3, 6, and 12 hours after each injection, and measured by isotope ratio mass spectrometry 13C 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 of which was dried and ground for soil carbon grouping.

[0053] The mineralized soil was separated into physical protection components, chemical protection components and biologically active components by wet sieving, and the content of each component in the soil 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 the element analyzer to convert the carbon and nitrogen in it into carbon dioxide and nitrogen, and then introduce the generated gas into the isotope ratio mass spectrometer through the gas transmission system to measure its isotope abundance.

[0054] 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. 15 The ratio of N to the total amount of nitrogen. It can be used to study the soil nitrogen cycle, such as nitrogen fixation, mineralization, nitrification and denitrification.

[0055] (5) Establish multi-source heterogeneous data The obtained soil physical protection component carbon (iPOM), chemical protection component carbon data (non-acid-degradable free silt component carbon (NH-dSilt), non-acid-degradable free clay component carbon (NH-dClay), non-acid-degradable closed silt component carbon (NH-μSilt), non-acid-degradable closed clay component carbon (NH-μClay), acid-degradable free silt component carbon (H-dSilt), acid-degradable free clay component carbon (H-dClay), acid-degradable closed silt component (H-μSilt) carbon and acid-degradable closed 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 The carbon content and proportion of carbon in the original soil and exogenous carbon, the direction and intensity of the excitation effect, the amount of mineralization and the rate of mineralization) are sorted into a table to record the collection location, succession stage, measurement time and other information of each sample. These data are then stored in a database to establish multi-source heterogeneous data, providing a basis for subsequent analysis and modeling.

[0056] S2. Integrate the multi-source heterogeneous data obtained in step S1 into the four-dimensional feature space of "physical-chemical-biological- 13 C abundance" through feature tensor decomposition, and dynamically adjust the contribution degrees of each modal data by using an adaptive weight allocation algorithm.

[0057] Preprocess the multi-source heterogeneous data obtained in the previous step and construct the four-dimensional feature space of "physical-chemical-biological- 13 C abundance", and the specific modal division is as follows: Physical modality: carbon content of iPOM components, 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 components, carbon content of fPOM components, isotope modality: 13 C abundance, 15 N abundance, the carbon of physical protection components, chemical protection components, and biological activity components respectively come from the content and proportion of native soil carbon and exogenous carbon, the direction and intensity of the priming effect, the mineralization amount and mineralization rate (time series data).

[0058] Data classification and standardization are as follows: The physical protection components include the iPOM organic carbon content of the 53-250 μm particle size fraction. After separation by dry sieving-sedimentation method, it is measured by an elemental analyzer, and the data is standardized to g C / kg; the chemical protection components are obtained by extracting 8 types of chemically bound carbon such as NH-dSilt, NH-dClay, and H-μSilt of the <53 μm particle size fraction, and then measured by hydrofluoric acid treatment and oxidation titration method, and the unit is unified as mgC / g; the biological activity components include the carbon content of cPOM and fPOM of the >250 μm particle size fraction, which are separated by density flotation method and measured by a TOC analyzer, and the data is normalized to a percentage; the isotope data include δ 13 C, δ 15 N abundance and the proportion of exogenous carbon, etc., are measured by an isotope ratio mass spectrometer and baseline correction is performed (such as based on the VPDB standard).

[0059] The process of constructing the four-dimensional tensor is as follows: First, the following dimension definitions are made: Physical modality (I1 = 3): >250 μm, 53-250 μm, <53 μm particle size fractions; Chemical modality (I2 = 8): NH-dSilt, NH-dClay, NH-μSilt, NH-μClay, H-dSilt, H-dClay, H-μSilt, H-μClay; Biological modality (I3 = 2): cPOM, fPOM; Isotope modality (I4 = 4): δ¹³C, δ 15 N, proportion of exogenous carbon, intensity of priming effect.

[0060] Then, the data at each sampling point is filled in as a tensor according to the four-dimensional coordinates. For example: X(i1, i2, i3, i4) = the value of the i1-th grain size - the i2-th chemical component - the i3-th biological component - the i4-th isotope parameter.

[0061] Then, the feature tensor is decomposed. In this embodiment, two feature tensor decomposition methods are provided, namely the Tucker decomposition method or the PARAFAC decomposition method, which decompose the multi-source heterogeneous data into a core tensor and factor matrices, and map the data to the four-dimensional feature space of "physics-chemistry-biology-isotope" through the transformation of the factor matrices.

[0062] Taking the Tucker decomposition method as an example: First, the tensor is decomposed into a core tensor G and factor matrices U (i) : ; Core tensor , representing the interaction intensity between modalities and the feature space after dimensionality compression. Factor matrix , indicating that the column vectors are orthogonal.

[0063] The specific steps are as follows: First, perform initialization and randomly generate ; Then, perform iterative optimization: fix other modalities and update : ; Similarly, update ; Set the convergence condition. For example, in this example, set the change rate of the objective function < 0.1% or reach the maximum number of iterations (100 times).

[0064] Taking PARAFAC (CP) decomposition as an example: PARAFAC decomposes the tensor into the sum of rank-1 components. For a four-dimensional tensor: ; Where: is the weight coefficient; are the rank-1 vectors of each modality; represents the outer product. represents the set of real numbers, i.e., the set that contains all rational and irrational numbers. R represents a tensor of a specific structure.

[0065] The specific decomposition steps are as follows: (1) Construction of the objective function Minimize the sum of squared residuals: ; (2) Optimization by Alternating Least Squares (ALS) Fix all variables except and update : ; Update , , , （1） represents the first-mode unfolding matrix of the tensor X.

[0066] (3) Normalization and weight calculation Normalize each vector: ; The weight .

[0067] Among them, each rank-1 component represents an independent influencing factor. For example, the coupling effect between the aggregate wrapping mode (physical mode) and a specific chemical protection mechanism (chemical mode). The weight reflects the contribution intensity of this factor to the overall data. represents the tensor product.

[0068] Example: For the four-dimensional tensor , with R = 3, the decomposition result: (Physical mode: dominated by iPOM); (Chemical mode: dominated by H-clay); = 12.5, indicating the proportion of the contribution of this factor to the data variation.

[0069] When implementing the dynamic adjustment of the contribution degrees of each modal data using the adaptive weight allocation algorithm, first dynamically adjust the weights based on the information entropy of the factor matrix: Calculate the variance contribution rate of each mode: ; Among them, is the variance of the k-th column of the factor matrix U (i) ; Modal information entropy: ; Adjust the weights according to the succession stage: ; Among them, S i is the stage adjustment coefficient, which can be set in the following way, for example: 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; When the succession stage is the arbor 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; The fused feature tensor after weighting

[0070] Among them, X i is the sub-tensor after standardization of each mode.

[0071] Project the original data into a four-dimensional space: F = ; The output tensor F ∈ , retaining more than 95% of the original information.

[0072] For the physical mode principal component, the first principal component usually corresponds to the synergistic effect of <53μm particle size and iPOM; for the chemical mode principal component, the first principal component can be interpreted as the adsorption competition relationship between H-μClay and NH-dSilt; The above method accurately depicts the chemical protection mechanism by subdividing subclasses such as NH-silt / NH-clay, and combines information entropy and succession stage coefficients to achieve intelligent adaptation of herbaceous stage heavy biological activity and arbor stage heavy chemistry.

[0073] S3, based on the spatial constraint of the boundary effect method and the statistical constraint of the least squares method, constructs a carbon saturation capacity prediction hybrid optimization model with the four-dimensional feature tensor and weight allocation results output in step S2 as the input, and outputs the maximum carbon sequestration capacity of <53μm microaggregates and the component carbon proportion; The specific implementation process is as follows: Receive the four-dimensional feature tensor F ∈ output in step S2 and the adaptive weight allocation result and extract the physical mode feature vector and the chemical mode feature vector ; ii Generate a particle size constraint matrix A according to the physical mode feature vector, and its diagonal element a ; Among them, is the physical modal weight; Based on the chemical modal eigenvector and weight , calculate the theoretical saturation threshold of the carbon content of each particle size :

[0074] Among them, β k is the regression coefficient, and C base is the laboratory calibration reference value; The sigmoid function is a non-linear activation function, and its mathematical expression is usually: ; Its role 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 convert the linear combination of chemical features into a proportionality coefficient, thereby dynamically adjusting the laboratory calibration reference value C base , and finally obtaining the theoretical saturation threshold .

[0075] When the chemical feature combination approaches positive infinity, sigmoid(z) → 1, and the threshold ≈ C base ; When z approaches negative infinity, sigmoid(z) → 0, and the threshold → 0 (in practice, negative values are avoided through experimental design); Through this function, the threshold is adaptively adjusted within the range of (0, C base ), realizing the inhibitory or enhancing effect of the chemical protection mechanism on the carbon capacity.

[0076] is the chemical modal eigenvector extracted by tensor decomposition (Tucker / PARAFAC) in step S2, representing the dimensionality-reduced features of the chemical protection components (such as the synergistic effects of 8 types of components such as NH-dSilt and H-μClay).

[0077] is the transpose of this eigenvector, used for dot product operation with the regression coefficient β k : ; Among them, 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 (such as <53μm) through linear combination.

[0078] C base To calibrate the reference value for the laboratory, the experiment is designed to determine through the following steps: Select typical karst soil samples (such as climax communities, shrub stages), and through constant temperature and humidity cultivation experiments, measure the carbon adsorption saturation values of each particle size (>250μm, 53 - 250μm, <53μm); control the environmental conditions (temperature 25°C, humidity 60%, no external carbon input) to ensure that the measurement results reflect the intrinsic carbon-holding capacity of the soil.

[0079] Take the average value of multiple experiments for the same particle size. For example, the C of the <53μm micro-aggregate base = 18.2 g·C / kg; through independent sample verification (such as cross-validation), ensure the generalization of the reference value (error <5%).

[0080] β k is the regression coefficient vector, determined by fitting historical data, and the specific process is as follows: Collect chemical characteristic data of multiple groups of karst soil samples and the measured saturation thresholds of the corresponding particle sizes ; the samples need to cover different succession stages (grassland, shrub, tree) and geomorphic types (peak cluster depressions, canyons, etc.).

[0081] Use Logistic regression or maximum likelihood estimation to optimize the following objective function: ; Introduce regularization (such as L2 regularization) to prevent overfitting. Select the optimal β through cross-validation k , requiring the coefficient of determination R2≥0.85; verify the prediction error on the independent test set (such as <10%).

[0082] The above steps are carried out through dynamic adjustment , so that the threshold changes adaptively with the intensity of the chemical protection mechanism; based on experimental calibration and regression fitting β k , ensure the ecological rationality of the model; The chemical characteristic vector covers data of multiple succession stages, supporting the migration application of the model from the grassland to the tree stage.

[0083] For example, in a certain karst peak cluster depression, through the above method, the of <53μm is calibrated to be 18.2 g·C / kg, and after fitting β k , the predicted threshold of the model = 16.5 g·C / kg (measured value 16.8 g·C / kg), with an error of only 1.8%.

[0084] Construct the following hybrid objective function: ; where A is the particle size constraint matrix, = 1, 2, 3] T is the predicted variable of carbon capacity of particle size fraction, b is the measured total organic carbon vector, is the penalty coefficient, is the theoretical saturation threshold of carbon content in each particle size fraction, is the particle size classification index parameter, = 1, 2, 3, corresponding to the particle size fraction index parameters of > 250μm, 53 - 250μm, < 53μm respectively, is the predicted variable of carbon capacity of the th aggregate; The dynamic penalty coefficient , is adaptively adjusted with the physical weight ; Adopt the hybrid strategy of Particle Swarm Optimization (PSO) and Genetic Algorithm (GA): First, initialize the data, generate 10 × R1 groups of initial solutions, and limit the search range to [0.5 , 1.2 in the computer system; Calculate the fitness, and define the weighted fitness function: ; where, is the fitness, is the physical mode weight, is the chemical mode weight; The particle swarm searches along the direction gradient in the physical mode space f phy ; the genetic algorithm performs arithmetic crossover and Gaussian mutation on the chemical mode parameter β k ; terminate when the change rate of the optimal solution fitness is < 0.05% for 5 consecutive generations; Output the maximum carbon capacity of < 53μm micro - aggregates and the proportion of component carbon: ; In the above method, PSO searches along the principal component direction in the physical mode space, improving the prediction efficiency of particle size capacity; GA mutates the chemical mode parameters, enhancing the self - adaptability of boundary constraints; multi - objective optimization is achieved through the weighted fitness function, and the weight allocation result { } directly regulates the balance between the statistical term and the boundary term.

[0085] Among them, the proportion output of component carbon provides a basis for adjusting the stability threshold in step S4 (for example, when the proportion > 50%, the critical value is reduced).

[0086] S4. Design a mineralization culture experiment based on the maximum carbon capacity predicted in step S3. Use isotope tracing to conduct mineralization culture on soils at different succession stages of karst forests. Through 13 C / 15 The dynamic ratio of exogenous carbon to native carbon and the input-output path of soil carbon are analyzed by the time-varying calibration curve of the C / N double standard sample, and the carbon stability disturbance threshold is calculated.

[0087] Among them, 13 C / 15 The construction method of the time-varying calibration curve of the C / N double standard sample includes: Based on the maximum carbon capacity of the <53μm micro-aggregates output in step S3, inject a mixed standard sample of 13C-glucose and 15N-ammonium sulfate into the mineralization culture system at a cycle of 12 hours. The injection time for each time is 5 minutes, and their concentrations are 15% of the original soil carbon content and 3% of the nitrogen content respectively; Set the exogenous carbon addition amount to 0.2 -0.5 . For example, = 20 g·C / kg, then add 4 - 10 g·C / kg of 13 C-glucose. Under normal circumstances, the nitrogen addition amount is 3% - 5% of the total nitrogen (TN) of the original soil. Use 15 N-ammonium sulfate to match the carbon-nitrogen stoichiometric ratio (such as C:N≈10:1).

[0088] Collect gas samples at 1, 3, 6, and 12 hours after each injection, and measure the 13 C and 15 N abundances by an isotope ratio mass spectrometer; collect soils on the 0th, 7th, 14th, 21st, and 28th days of cultivation, separate different particle sizes (>250μm, 53 - 250μm, <53μm), and measure the δ 13 C and δ 15 N values of each component.

[0089] Establish a carbon-nitrogen coupled kinetic model; among them, 13 The C abundance uses an exponential decay equation to characterize carbon turnover, 15 The N abundance uses a logarithmic growth equation to characterize nitrogen fixation; Synchronously fit the carbon turnover rate constant and the characteristic time of nitrogen fixation through a simulated annealing-least squares hybrid algorithm to complete the construction of the 13 C / 15 time-varying calibration curve of the C / N double standard sample.

[0090] The carbon-nitrogen coupled kinetic model includes the following coupled equations: 13 C decay equation:

[0091] where, is the carbon turnover rate constant, which is inversely solved by the simulated annealing algorithm; 15N adsorption equation:

[0092] where, is the characteristic time of nitrogen retention, the smaller it is, the faster the retention rate (e.g., clay may have <sand). , are soil type parameters (such as organic matter content, mineral composition). α is positively correlated with the nitrogen adsorption capacity, and β represents the initial adsorption background value.

[0093] The process of analyzing the carbon flow path and quantifying the priming effect is carried out according to the following steps: First, it is necessary to determine the input and output paths.

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

[0095] The native carbon output rate (R out ) represents the total mineralization amount minus the exogenous part.

[0096] The carbon flow balance ratio is calculated according to the following formula: ; If η > 1.2, it is marked as a carbon net loss path; if η < 0.8, it is marked as a carbon net fixation path.

[0097] The priming effect quantization is carried out according to the following steps: Calculate the cumulative priming effect (CPE): ; is the microbial diversity index (Shannon index based on 16S rRNA sequencing). The lower the diversity, the more significant the priming effect. The integration interval is from the initial time t0 to the critical time t c (usually at the end of the incubation period, such as 28 days).

[0098] Finally, the carbon stability perturbation threshold is determined. First, a dynamic threshold needs to be set, and the carbon stability perturbation threshold T is adjusted according to the <53μm component carbon ratio (P<53)crit : ; Among them, P<53 represents the carbon proportion of the <53 μm microaggregate component in the soil (unit: %), which characterizes the physical protection intensity of carbon. The stronger the physical protection ability of the microaggregates (<53 μm) to carbon (the higher P<53), the more stable the carbon pool, and the lower the carbon stability perturbation threshold T crit The lower it is. 2.5 is a coefficient determined by historical data regression, reflecting the empirical relationship between the carbon loss rate and the microaggregate protection ability.

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

[0100] In this embodiment, when the following conditions are simultaneously met, it is determined that the carbon pool is unstable: The priming effect exceeds the limit: CPE > Tcrit; Isotope offset verification: The δ13C offset exceeds the background value by ±2‰ (excluding instrument error interference).

[0101] To ensure the reliability of the threshold, the following verification is required: For example, sensitivity analysis, by perturbing the physical protection weight and the chemical protection weight by ±10%, it is required that the coefficient of variation of the threshold CV < 8%. Field evidence can also be combined for synchronous verification in peak cluster depressions, canyons, and basins, and it is required that the spatial matching degree between the predicted and measured carbon loss hotspots ≥ 85%≥85%.

[0102] S5. Construct a digital twin of the carbon cycle to simulate the carbon saturation deficit value and spatial heterogeneity of soil at different succession stages of karst forests.

[0103] The construction of the digital twin includes: Establish a three-dimensional carbon pool interaction map, and map the carbon in the physical protection component, the carbon in the chemical protection component, and the carbon in the bioactive component in a vector field in three-dimensional space according to the soil at different succession stages of herbs, shrubs, and trees; embed the LSTM-Transformer hybrid neural network, the input layer receives data on cumulative carbon input, the carbon mass of the <53μm component, and the carbon content of the <53μm component, and the output layer predicts the carbon saturation deficit value. Combining the existing carbon content and carbon sequestration rate, the carbon saturation capacity time limit is calculated.

[0104] The specific implementation process is as follows: First, a three-dimensional carbon pool interactive map was constructed. The main goal was to dynamically map the physical, chemical, and biologically active components of the soil carbon pool in three-dimensional space and quantify the correlation between carbon distribution and succession stage. Three-dimensional topographic data of karst landforms (resolution ≤ 1m) were obtained using UAV LiDAR or high-precision terrain scanners. Sample areas were divided according to succession stages (herbs, shrubs, and trees), and at least 50 sampling points were set up in each sample area.

[0105] 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.

[0106] Then, a three-dimensional space grid is constructed based on Python's PyVista or Paraview, with a grid unit size of ≤10m³, to achieve three-dimensional vector field mapping. The mapping method is 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 the physical protective carbon of a point = 15g / kg, the chemical protective carbon = 8g / kg, and the biologically active carbon = 2g / kg, then its vector coordinates are (15, 8, 2).

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

[0108] 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 ).

[0109] 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.

[0110] For the input layer, its characteristic dimensions include cumulative carbon input (total historical carbon input, unit: t·C / ha), carbon mass of components <53μm (unit: g·C / kg) and carbon proportion of components <53μm (unit: %). For time series processing, the input data are arranged by time step (such as year), with the longest retrospective period of 20 years.

[0111] For the LSTM layer (Long Short-Term Memory Network) in the network structure, set the number of neurons to 64 to capture the temporal dependence of carbon input (such as the carbon sequestration lag effect). For the Transformer layer in the network structure, set the multi-head attention mechanism (4 heads) to focus on spatial heterogeneity (such as the impact of micro-aggregate distribution on carbon saturation). For the fusion module in the network structure, fuse the temporal features output by LSTM and the spatial features of Transformer through Cross-Attention.

[0112] For the output layer, define the carbon saturation deficit value (ΔC) as the difference between the predicted current carbon content and the saturation capacity (unit: g·C / kg). The calculation of the carbon saturation capacity time limit (T sat ) is combined with the carbon sequestration rate (k sat , unit: g·C / kg / yr) to calculate the time required to reach saturation: ; Model training and validation are carried out according to the following steps: First, perform dataset division. The training set (70%) covers typical sample plots at each succession stage of herbs, shrubs, and trees. The validation set (15%) is used for hyperparameter tuning. The test set (15%) evaluates the model's generalization ability.

[0113] Adopt the Huber loss function to balance the mean squared error (MSE) and the mean absolute error (MAE): ; In this embodiment, 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.

[0114] Set the verification accuracy requirements as follows: the prediction error of the carbon saturation deficit value ≤ 10% (R² ≥ 0.85); the error of the carbon saturation time limit ≤ 15% (RMSE ≤ 2 years). For the sensitivity test, set the perturbed input parameters (such as the proportion of <53μm carbon ± 5%), and the model output fluctuation should ≤ 8%.

[0115] When deploying and applying the digital twin, first perform real-time data fusion. Real-time collect soil humidity, temperature, and microbial activity data through Internet of Things (IoT) sensors, and dynamically update the bioactive carbon component in the three-dimensional map. Synchronize the model prediction results every 6 hours to generate a carbon saturation heat map (red: high deficit; green: low deficit).

[0116] Set the warning threshold. For example, when ΔC > the critical value (such as ΔC > 5g·C / kg) and T satWhen it is < 10 years, trigger the alert for the priority area of carbon sequestration. If ΔC is high and the proportion of microaggregates is low, the system automatically recommends the application of biochar (to enhance physical protection); if ΔC is low but the proportion of chemical protection is high, the system automatically suggests planting deep-rooted plants (to enhance mineral binding).

[0117] The above shows and describes the basic principles, main features and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited by the above embodiments. What is described in the above embodiments and the specification only illustrates the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention will have various changes and improvements, and these changes and improvements all fall within the scope of the present invention claimed. The scope of protection claimed by the present invention is defined by the appended claims and their equivalents.

Claims

1. A multimodal data fusion method for extracting soil carbon cycle features in karst forests, characterized in that: The following steps are involved: S1, obtain the carbon data of soil physical protection components, chemical protection components, biologically active components and isotope tracer data at different succession stages of karst forests, and establish multi-source heterogeneous data; S2, integrating the multi-source heterogeneous data obtained in step S1 into the "physical-chemical-biological-isotopic" four-dimensional feature space through feature tensor decomposition, and dynamically adjusting the contribution of each modal data using an adaptive weight allocation algorithm; S3, based on the spatial constraints of the boundary effect method and the statistical constraints of the least squares method, the four-dimensional feature tensor and weight distribution results output by step S2 are used as input to construct a hybrid optimization model for predicting carbon saturation capacity, and output the maximum carbon sequestration capacity of microaggregates <53 μm and the carbon proportion of the components; S4, designing a mineralization culture experiment based on the maximum carbon capacity predicted in step S3, and using isotope tracing to mineralize soils at different succession stages of karst forests. 13 C / 15 The time-varying calibration curve of the N double standard samples 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.

2. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 1 is characterized in that: The following steps are also included: 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 in different succession stages of karst forest.

3. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 1 or 2, characterized in that: In step S1, the soil physical protection component carbon data is 53-250 μm carbon content, including physical protection component carbon iPOM; The chemical protection component carbon data is the carbon content of <53μm, including non-acid decomposition free powder component carbon NH-dSilt, non-acid decomposition free clay component carbon NH-dClay, non-acid decomposition closed accumulation powder component carbon NH-μSilt, non-acid decomposition closed accumulation clay component carbon NH-μClay, acid decomposition free powder component carbon H-dSilt, acid decomposition free clay component carbon H-dClay, acid decomposition closed accumulation powder component carbon H-μSilt and acid decomposition closed accumulation clay component carbon H-μClay; 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; The isotope tracing data includes 13 C abundance, 15 N abundance, physical protection component carbon, chemical protection component carbon, and biologically active component carbon respectively come from the content and proportion of original soil carbon and exogenous carbon, the direction and intensity of the excitation effect, and the mineralization amount and mineralization rate.

4. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 1 or 2, 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, and the data is mapped to the "physical-chemical-biological-isotope" four-dimensional feature space through the transformation of the factor matrix.

5. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 1 or 2, characterized in that: The carbon saturation capacity prediction hybrid optimization model described in step S3 is constructed according to the following steps: 1) receiving the four-dimensional feature tensor and the adaptive weight allocation 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 eigenvector, and calculate the theoretical saturation threshold of carbon content in each particle size based on the chemical modal eigenvector and weight; 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 index parameters of particle size >250μm, 53-250μm, <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.

6. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 1 or 2, characterized in that: In step S4 13 C / 15 The method for constructing a time-varying calibration curve of N double 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 simulated annealing-least squares hybrid algorithm was used to simultaneously fit the carbon turnover rate constant and the nitrogen fixation characteristic time. 13 C / 15 N time-varying calibration curves of double standards.

7. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 6 is 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.

8. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 7, 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 the microbial diversity index. The carbon stability disturbance threshold is determined based on the carbon proportion of components <53 μm, and whether the carbon pool has entered an unstable state is determined by the carbon stability disturbance threshold and the normalized cumulative excitation effect; Wherein, the microbial diversity index includes the Shannon index based on 16S rRNA sequencing.

9. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 2, 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 according to the soil at different succession stages of herbs, shrubs and trees; 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.

10. The method for extracting soil carbon cycle characteristics of karst forests by multimodal data fusion according to claim 5, characterized in that: The hybrid strategy of particle swarm optimization and genetic algorithm to output the maximum carbon capacity of microaggregates <53 μm and the carbon proportion of components includes the following steps: Generate multiple sets of initial solutions, limit the search range in the computer system, and define the weighted fitness function as follows: ； 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