A carbon fixation and yield increase synergistic optimization method based on crop root zone carbon and nitrogen regulation mechanism

Abstract

This invention relates to the field of computer-aided agriculture, specifically disclosing a synergistic optimization method for carbon sequestration and yield enhancement based on the carbon and nitrogen regulation mechanism of crop root zones. The method includes: utilizing a root zone carbon and nitrogen multi-stable-state system model that integrates the mapping relationship between functional gene expression and biochemical reaction rates to diagnose the current state and plan candidate state migration paths; performing inverse sensitivity analysis through a root zone dynamic response model, identifying key regulatory windows by combining microbial functional temporal characteristics, and querying a measure-process response knowledge graph to match precise leverage measures; generating a set of synergistic regulation schemes through multi-objective optimization including functional gene expression synergy; constructing a high-fidelity digital twin for deduction and early warning, and using ensemble Kalman filtering to assimilate measured data to drive model self-evolution. This invention achieves dynamic inverse design, intelligent path planning, and continuous adaptive optimization of the complex carbon and nitrogen system in the root zone, effectively synergistically improving farmland carbon sequestration capacity and crop yield.

Description

Technical Field

[0001] This invention relates to the field of computer-aided smart agriculture technology, specifically to a synergistic optimization method for carbon sequestration and yield increase based on the carbon and nitrogen regulation mechanism of crop root zones. Background Technology

[0002] Currently, agronomic management programs aimed at achieving carbon sequestration and emission reduction in farmland without reducing yields mainly rely on three technical approaches, but all of them have fundamental limitations.

[0003] The most common static fertilization schemes based on field trials pre-determine several ratios of exogenous carbon and nitrogen reduction, and then compare the static results at the end of the season to select the optimal combination. This type of method completely fails to make dynamic and reverse time-series designs based on the stage-specific needs of the crop throughout its entire growth period and the critical windows of root zone processes. This results in a mismatch between carbon and nitrogen supply in time and crop needs and microbial activity, making it difficult to synergistically optimize carbon sequestration and yield increase goals.

[0004] Some studies employ mechanistic or statistical models for simulation analysis. However, most of these models treat farmland ecosystems as homogeneous or monostable systems, failing to characterize their "multistable" characteristics as complex systems (such as chemical nitrogen fertilizer dependence and organic matter enrichment). Therefore, existing methods cannot diagnose the current state of the system, let alone plan a transition path from the current state to the ideal state that is least energy-intensive and most robust, resulting in management measures often being ineffective or costly.

[0005] Although some preliminary digital agricultural systems have attempted to integrate monitoring and models, their model parameters remain fixed once calibrated, lacking the ability to continuously evolve. Farmland environments and soil biological communities are constantly changing, causing static models and optimization schemes built based on historical data to rapidly decline in predictive and control effectiveness when applied in new environments or over a long period. This prevents the formation of an intelligent closed loop of "practice-learning-evolution," resulting in insufficient technological universality and long-term effectiveness.

[0006] Therefore, there is an urgent need for an innovative method that can dynamically reverse-design the multi-steady-state carbon and nitrogen processes in the root region, plan the system transformation path, and possess continuous self-evolution capabilities. Summary of the Invention

[0007] The purpose of this invention is to provide a synergistic optimization method for carbon sequestration and yield increase based on the carbon and nitrogen regulation mechanism of crop root zone, so as to solve the fundamental problems in the existing technology that static fertilization schemes cannot perform dynamic reverse design and transformation path planning for the multi-steady-state carbon and nitrogen system in the root zone, and that models and strategies cannot continuously self-evolve with changes in the farmland system.

[0008] To solve the above-mentioned technical problems, the present invention specifically provides the following technical solution: A method for synergistic optimization of carbon sequestration and yield enhancement based on the carbon and nitrogen regulation mechanism of crop rhizosphere includes the following steps: S1. Based on multi-source data, the current state is diagnosed by integrating the functional gene expression-biochemical reaction rate mapping relationship into a root region carbon and nitrogen multi-stable system model, and the optimal candidate state migration path is planned. S2. Using the root region dynamic response model, perform inverse sensitivity analysis on the candidate state migration path, and combine the temporal characteristics of functional gene expression to identify key regulatory windows and match combinations of leverage measures. S3. Integrate the candidate state transition paths, key regulatory windows and leveraging measures to optimize for multiple objectives, including functional gene expression synergy, and generate a set of synergistic regulatory schemes. S4. Construct a high-fidelity digital twin based on the aforementioned coordinated control scheme set for simulation, evaluate the scheme by monitoring critical transition early warning signals, and use the results to evolve the core model and feed it back to the previous steps.

[0009] As a preferred embodiment of the present invention, S1 specifically includes: S11. The functional gene expression-biochemical reaction rate mapping relationship established based on the current root zone carbon and nitrogen system status monitoring data, historical management data, and rhizosphere microbial metatranscriptome sequencing data of the target farmland is input into the root zone carbon and nitrogen multi-stable system model. S12. The root zone carbon and nitrogen multi-steady-state system model diagnoses the current system state and calculates its distance from the target steady-state state basin boundary by coupling the nonlinear dynamic equations and mapping relationships between the soil organic carbon pool, microbial functional groups and inorganic nitrogen pool. S13. Starting from the current state and ending at the preset carbon sequestration and production enhancement synergistic target state, in the model phase space, the state-action-state chain search algorithm is used to plan at least one candidate state transition path, with energy consumption quantified by agronomic energy equivalent and robustness measured by the average distance from the state basin boundary as the comprehensive cost function.

[0010] As a preferred embodiment of the present invention, S12 specifically includes: S121. The root zone carbon and nitrogen multi-steady-state system model couples the nonlinear dynamic relationship between the soil organic carbon pool, microbial functional groups and inorganic nitrogen pool through the first set of differential equations, and embeds the functional gene expression-biochemical reaction rate mapping relationship through the second set of functions, so as to dynamically transform the expression level of key functional genes into the constraint condition of corresponding biochemical reaction flux. S122. Based on the root region carbon and nitrogen multi-stable system model, input the current monitoring data, solve the coupling equations and functions by numerical integration, and determine the coordinate position of the current system state in the model phase space; S123. In the phase space, based on potential energy landscape analysis, the boundary of the state basin where the preset target steady state is located is identified, and the shortest Euclidean distance or potential energy difference from the current state coordinate point to the boundary of the state basin is calculated as a quantitative value of the distance of the state basin boundary, which is used to characterize the ease and stability of the system's migration from the current state to the target steady state.

[0011] As a preferred embodiment of the present invention, S2 specifically includes: S21. Input the candidate state transition path into the root region dynamic response model, which embeds agronomic process constraints and couples functional gene expression time series data; S22. Perform inverse sensitivity analysis on the root region dynamic response model, traverse potential control time points on the simulation timeline and inject virtual control pulses of unit intensity to generate a sensitivity matrix characterizing the influence intensity of control time point-state transition node. S23. Based on the sensitivity matrix and the preset control cost function, perform a benefit-cost ratio analysis to obtain a list of key control windows; S24. Based on the list of key control windows, query the pre-constructed measure-process response knowledge graph, match and output one or more specific ratios and application methods of exogenous carbon input and chemical nitrogen fertilizer regulation that can most efficiently drive the system state to evolve to the next key state transition node, thus forming a combination of leverage measure types.

[0012] As a preferred embodiment of the present invention, S22 specifically includes: S221. On the timeline of the dynamic response model of the root region, traverse all potential control time points that meet the feasibility of agricultural operations, as well as several key state transition nodes defined on the path. S222. In the root zone dynamic response model, at each of the regulation time points, a virtual regulation pulse of unit intensity is injected, defined by the combination of standardized exogenous carbon input and chemical nitrogen fertilizer regulation. S223. After injecting each of the virtual control pulses, run the root region dynamic response model forward to the end of the path; accurately record the actual simulated arrival state of each key state transition node, and compare it with the baseline arrival state when no pulse is injected, and calculate the change in state deviation caused by each pulse to each node. S224. Using the potential regulation time as the row index of the matrix, and each of the key state transition nodes as the column index of the matrix, the calculated deviation change is filled into the corresponding position to construct a sensitivity matrix, wherein the matrix element values ​​characterize the influence strength of applying standard intervention at a specific time point on achieving state transition at a specific node.

[0013] As a preferred embodiment of the present invention, S24 specifically includes: S241. For each window in the list of key control windows, extract its time location, the associated next key state transition node, and the direction and magnitude of the corresponding target state variable change; S242. The extracted information is used as a query request and input into a pre-constructed measure-process response knowledge graph; the knowledge graph stores different types of exogenous carbon input and chemical nitrogen fertilizer regulation methods as nodes, and stores their directional regulation efficiency parameters on soil humification, nitrification-denitrification, and microbial carbon pump activation of specific biogeochemical processes as edges; S243. Execute graph traversal and matching algorithms in the knowledge graph to find a combination of measures that can drive the system state from the current node to the next key state transition node with the highest prediction efficiency; S244. The matched one or more exogenous carbon inputs are parameterized and encapsulated with the specific ratio and application method of chemical nitrogen fertilizer regulation to form a combination of leverage measures that serve the window and output it.

[0014] As a preferred embodiment of the present invention, S3 specifically includes: S31. The candidate state transition paths, the list of key control windows corresponding to each path, and the associated leverage measures are paired and integrated to form several path-control decision chains; S32. Establish a multi-objective optimization model. Decision variables include the selection of candidate state transition paths and the determination of the specific types and amounts of measures in the associated leverage combination for each key control window. Optimization objectives include at least minimizing total economic cost, minimizing path tracking error, maximizing meteorological risk robustness, and maximizing the synergistic expression of carbon and nitrogen turnover functional genes. S33. Apply a multi-objective optimization algorithm to solve the multi-objective optimization model and output the Pareto optimal solution set; S34. Decode the Pareto optimal solution set to generate a set of coordinated control schemes that include specific migration paths, control timing, types of measures, and application amounts.

[0015] As a preferred embodiment of the present invention, S32 specifically includes: S321. The set of decision variables includes: a first variable, used to select a specific candidate state transition path from the path-regulation decision chain; and a second set of variables, used to determine the specific type, amount, method of application, and ratio and amount of chemical nitrogen fertilizer in the combination of leverage measures associated with each key regulation window of the selected path. S322. The optimization objective functions include: a first objective function, used to calculate and minimize the total economic cost of implementing all selected measures, including material costs, operating costs, and external costs of environmental risks; a second objective function, used to calculate and minimize the path tracking error, i.e., the weighted Euclidean distance between the simulated state of the system at each key state transition node and the preset target state of the selected path after simulating the decision variables through the root zone dynamic response model; a third objective function, used to evaluate and maximize the robustness to meteorological risks, i.e., the probability or average distance that the simulated system state remains within the target steady-state basin under preset extreme meteorological scenarios; and a fourth objective function, used to evaluate and maximize the synergy of carbon and nitrogen turnover functional gene expression, the value of which is calculated based on the rhizosphere microbial functional gene expression profile obtained from the simulation, reflecting the degree of synergy between enhanced expression of carbon fixation and nitrogen fixation genes and suppressed expression of carbon mineralization and denitrification genes. S323. The constraints of the multi-objective optimization model include at least the agronomic operation period, the upper limit of a single fertilization amount, and the total external carbon input resource limit.

[0016] As a preferred embodiment of the present invention, S4 specifically includes: S41. For at least one preferred scheme selected from the set of coordinated control schemes, and in combination with real-time meteorological and soil data, construct a high-fidelity digital twin that is highly mirrored with the target farmland; S42. Drive the high-fidelity digital twin to perform dynamic simulations under baseline and extreme weather scenarios, and monitor the critical transition warning signals of the system in real time to evaluate the robustness of the scheme; S43. Compare the extrapolated and predicted data with the measured root zone state and microbial macrotranscriptome data of the target farmland, and use the ensemble Kalman filter assimilation algorithm to automatically calibrate and evolve the core parameters of the root zone carbon and nitrogen multi-stable system model and the root zone dynamic response model. S44. Feedback the evolved model to conduct a new round of system diagnosis and path planning.

[0017] As a preferred embodiment of the present invention, S43 specifically includes: S431. Extract the observation data corresponding to the time and variables from the measured root zone state and microbial macrotranscriptome data of the target farmland to form an observation vector; S432. Extract the predicted data corresponding to the pilot monitoring time from the high-fidelity digital twin simulation results to form the model prediction vector; calculate the difference between the observation vector and the model prediction vector; S433. With minimizing the model prediction bias as the optimization objective, the core parameters to be optimized in the root zone carbon and nitrogen multi-steady-state system model and the root zone dynamic response model are defined as vectors to be estimated. The core parameters include at least the organic matter decomposition rate constant, microbial carbon use efficiency, nitrification activation energy and the regulatory coefficient of functional gene expression on the reaction rate. S434. The ensemble Kalman filter assimilation algorithm is used to iteratively perform the following operations: predict the model state based on the current set of parameter vectors; calculate the covariance matrix between the prediction set and the observation vectors; update the parameter set using the Kalman gain formula so that the updated parameter set minimizes the deviation between the prediction data and the observation data in a statistical sense. S435. The mean of the updated parameter set obtained after assimilation convergence is used as the new parameter values ​​after automatic calibration and evolution, and updated into the two core models; based on the residual analysis of the assimilated models, it is determined whether there is a systematic bias pattern, and if so, an instruction to optimize the model structure is triggered.

[0018] Compared with the prior art, the present invention has the following advantages: 1. By introducing state basin theory and inverse sensitivity analysis, a leap from static management to dynamic system path planning has been achieved. This method can accurately diagnose the steady state of the system and plan the migration path with the lowest energy consumption and the strongest robustness, solving the inertia problem of complex ecosystem transformation and realizing the active guidance of carbon and nitrogen processes in the root zone.

[0019] 2. By constructing a knowledge graph of the mapping relationship between functional gene expression and biochemical reaction rate, and between measures and process responses, a shift from extensive experience-based to precise reverse regulation was achieved. This method intelligently links microbial molecular functions with agronomic measures, enabling the matching of the most efficient leverage measures within the optimal time window, thus achieving the optimal allocation of management resources in time and space.

[0020] 3. By integrating critical transition early warning monitoring with an ensemble Kalman filter assimilation and evolution mechanism, an intelligent closed-loop system with forward-looking and self-evolving capabilities was constructed. This system can provide early warnings of system risks and continuously drive model self-optimization using multi-source measured data, ensuring the long-term adaptability and effectiveness of the technical solution. (See attached figures for details.) To more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings in the following description are merely exemplary, and those skilled in the art can derive other embodiments based on the provided drawings without creative effort.

[0021] Figure 1 This is a schematic diagram summarizing the process of the method described in Embodiment 1 of the present invention.

[0022] Figure 2 This is a detailed flowchart of the method described in Embodiment 1 of the present invention. Detailed Implementation

[0023] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0024] The concepts involved in this application will first be described with reference to the accompanying drawings. It should be noted that the following descriptions of various concepts are only for the purpose of making the content of this application easier to understand and do not constitute a limitation on the scope of protection of this application; furthermore, the embodiments and features in the embodiments of this application can be combined with each other unless otherwise specified. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0025] Example 1 like Figure 1 - Figure 2 As shown, this invention provides a synergistic optimization method for carbon sequestration and yield enhancement based on the crop root zone carbon and nitrogen regulation mechanism, comprising the following steps: S1. Based on multi-source data, a root region carbon and nitrogen multi-stable state system model integrating functional gene expression-biochemical reaction rate mapping relationships is used to diagnose the current state and plan the optimal candidate state transition path; specifically including: S11. Construct the comprehensive data input foundation required for the operation of the root region carbon-nitrogen multi-stable-state system model. The specific implementation process is as follows: S111. Real-time monitoring data on the current state of the carbon and nitrogen system in the root zone of the target farmland is collected through a multi-parameter sensor network deployed in the root zone. This monitoring data includes indicators characterizing the dynamics of the soil organic carbon pool, including soil organic carbon content, soluble organic carbon concentration, and soil respiration rate; indicators characterizing the activity of microbial functional groups, including microbial biomass carbon and microbial biomass nitrogen; and indicators characterizing the state of the inorganic nitrogen pool, including ammonium nitrogen concentration, nitrate nitrogen concentration, and the nitrogen fixation and release flux of soil microorganisms.

[0026] S112. Retrieve historical management data of the target farmland as supplementary input to the initial conditions of the model. Historical management data includes, but is not limited to, records of exogenous carbon inputs over the years, specifically the types, application amounts, and application time series of different carbon sources such as biochar, straw return to the field, and organic fertilizer; the history of chemical nitrogen fertilizer application, including the application amount, application period, and application method of urea, ammonium nitrogen fertilizer, and nitrate nitrogen fertilizer; and historical records of crop planting systems and irrigation systems.

[0027] S113. Obtain rhizosphere microbial metagenomic sequencing data and establish a mapping relationship between functional gene expression and biochemical reaction rates. Specifically, representative samples were collected from the rhizosphere soil region, total RNA was extracted, and a metagenomic sequencing library was constructed. High-throughput sequencing was used to obtain expression data of key functional genes involved in carbon and nitrogen cycling. This data specifically includes: expression levels of carbon fixation genes such as cbbL and cbbM, which regulate the conversion of inorganic carbon to organic carbon; expression levels of carbon mineralization genes such as amyA and lip, which regulate the decomposition and mineralization of organic carbon; expression levels of nitrogen fixation genes such as nifH, which regulate the conversion of atmospheric nitrogen to available nitrogen in the soil; expression levels of nitrification genes such as amoA and nxrA, which regulate the oxidation of ammonium nitrogen to nitrate nitrogen; and expression levels of denitrification genes such as nirK and nosZ, which regulate the reduction and loss of nitrate nitrogen to gaseous nitrogen.

[0028] S114. Based on the above metatranscriptome sequencing data, establish a mapping relationship between functional gene expression and biochemical reaction rates. This mapping relationship is achieved through the following mechanism: Based on the gene expression level TPM value obtained from metatranscriptome sequencing, the theoretical maximum reaction rate of the corresponding enzyme reaction is calculated by combining the Michaelis equation or flux balance analysis (FBA) method. The theoretical rate was then modified based on actual soil environmental factors to establish a dynamic correlation function between gene expression, environment, and response rate. Finally, this correlation function was used as a constraint on biochemical reaction flux and embedded into the differential equation set of the root zone carbon and nitrogen multi-steady-state system model to realize the quantitative transformation of microbial functional gene expression level into the rate of ecosystem material cycling process.

[0029] S115. Input the current carbon and nitrogen system status monitoring data collected above, the historical management data retrieved, and the functional gene expression and biochemical reaction rate mapping relationship established based on rhizosphere microbial metatranscriptome sequencing data into the root zone carbon and nitrogen multi-stable system model.

[0030] During input, monitoring data serves as the initial value of the state variable describing the current system state, historical management data serves as the prior parameter of the model boundary condition, and the mapping relationship between functional gene expression and biochemical reaction rate serves as the interface function that couples the nonlinear dynamic equations and microbial molecular biological processes within the model, thereby constructing a cross-scale data-driven foundation from gene expression to ecosystem function.

[0031] S12. The diagnostic mechanism of the root region carbon and nitrogen multi-stable system model for the current system state and the quantification of its topological relationship with the target steady state are as follows: S121. The core architecture for constructing a multi-steady-state carbon and nitrogen system model in the root zone is as follows: This model comprises three core state variable modules: soil organic carbon pool, microbial functional group, and inorganic nitrogen pool. The model couples the nonlinear dynamic relationships between these three modules through a first set of differential equations: The soil organic carbon pool module distinguishes between active organic carbon and inert organic carbon, and exchanges materials with the microbial functional group module through carbon mineralization and stabilization processes. The microbial functional group module distinguishes between bacterial biomass and fungal biomass, mediates the conversion of organic carbon, and regulates the fixation and release of inorganic nitrogen. The inorganic nitrogen pool module distinguishes between ammonium nitrogen and nitrate nitrogen, forming a feedback loop with the microbial functional group module through the nitrification-denitrification process. The aforementioned set of differential equations describes the nonlinear feedback dynamics between the three modules through carbon and nitrogen flux terms, microbial growth and death rate terms, and environmental response function terms.

[0032] The functional gene expression and biochemical reaction rate mapping relationship established in step S11 is embedded into the model as the second set of function constraints. This embedding mechanism dynamically calculates the kinetic parameters of the corresponding biochemical reaction process, such as the organic matter decomposition rate constant and the maximum nitrification rate, by using the expression levels of carbon fixation genes, carbon mineralization genes, nitrogen fixation genes, nitrification genes, and denitrification genes obtained from metatranscriptome sequencing, through a preset gene expression-enzyme activity-reaction rate conversion function.

[0033] These parameters, which are regulated in real time by gene expression levels, serve as dynamic constraints and are input into the first set of differential equations. This enables cross-scale coupling of microbial functional gene expression data with ecosystem material cycling processes, allowing the model to reflect the real-time impact of rhizosphere microbial community functional status on carbon and nitrogen transformation processes.

[0034] S122. Based on the input current monitoring data, the coupled differential equations and functional constraints mentioned above are solved by numerical integration to determine the coordinate position of the current system state in the high-dimensional phase space constructed by the model. This phase space uses the key state variables of the soil organic carbon pool, microbial functional groups, and inorganic nitrogen pool as orthogonal coordinate axes to form a multi-dimensional topological space characterizing the dynamic characteristics of the carbon and nitrogen system in the root zone.

[0035] After determining the current state coordinates, the model performs dynamic state diagnosis by calculating the eigenvalue distribution of the system's Jacobian matrix: if all the real parts of the eigenvalues ​​are negative, the current system is determined to be in a steady state, i.e., a stable state with self-recovery capability; if there are eigenvalues ​​with positive real parts, the current system is determined to be in a transient state, i.e., an unstable state facing a critical point of state transition.

[0036] Meanwhile, the state basin to which the current system state belongs is identified through potential energy landscape analysis. This potential energy landscape is constructed from the potential energy function of the model equation and characterizes the topological structure of the attraction domain.

[0037] S123. Calculate the state basin boundary distance between the current system state and the target steady state. In phase space, the preset carbon sequestration and production enhancement synergistic target state corresponds to an ideal state basin. The model determines the state basin boundary surface where the target steady state is located through potential energy landscape analysis.

[0038] The shortest Euclidean distance between the current state coordinates and the boundary of the basin in that state, or the potential energy difference between the two points, is calculated as the quantified value of the basin boundary distance. This distance value characterizes the ease and stability of the system's migration from the current state to the target steady state: a larger distance value indicates that the system is deep within the current state basin, requiring a higher energy barrier to migrate to the target state; a smaller distance value indicates that the system is close to the basin boundary or in a critical region, requiring less control energy for migration but facing the risk of transitioning to a non-target state. This quantified distance serves as a core indicator for evaluating the robustness of the path planning algorithm in subsequent steps, providing a topological basis for optimizing candidate state migration paths.

[0039] S13. The planning process for candidate state transition paths is implemented as follows: S131. Determine the spatial dimensions and start / end coordinates of the path search: The current system state determined in step S12 is used as the starting point of the path. This starting point is defined in phase space by the specific numerical coordinates of the current soil organic carbon pool (distinguishing between active organic carbon and inert organic carbon), the microbial functional group activity index (distinguishing between bacterial biomass and fungal biomass), and the inorganic nitrogen pool concentration (distinguishing between ammonium nitrogen and nitrate nitrogen).

[0040] The predetermined synergistic target state for carbon sequestration and yield enhancement is used as the endpoint of the path. This endpoint is defined by the central coordinates of the ideal range of the target soil organic carbon storage capacity interval, the functional microbial community structure index, and the mineral nitrogen supply intensity. In the multidimensional phase space composed of the above three modules, a high-dimensional topological space characterizing the dynamic features of the carbon and nitrogen system in the root zone is constructed with each state variable as an orthogonal coordinate axis.

[0041] S132. Apply the state-action-state chain search algorithm in phase space for path exploration: The algorithm abstracts agronomic measures into a discrete action space, which includes the application amount and timing of specific types of exogenous carbon inputs (biochar, straw return to the field, organic fertilizer), and the adjustment amount and application method of chemical nitrogen fertilizers (urea, ammonium nitrogen fertilizer, nitrate nitrogen fertilizer).

[0042] The algorithm starts from the path's origin and iteratively selects actions using a Monte Carlo tree search strategy. Each action acts on the current state node, driving quantitative transformations of the soil organic carbon pool, microbial functional groups, and inorganic nitrogen pool within the phase space described by the root zone carbon and nitrogen multi-stable system model. Based on the model's coupled differential equations, the algorithm calculates the evolution trajectory of state variables, generating new state nodes and forming a chain structure of state-action-new state. Through a search strategy combining depth-first and breadth-first search, multiple potential state transition chains are constructed from the path's origin to its destination. Each chain consists of a series of ordered state nodes and action edges connecting these nodes.

[0043] S133. Establish a dual-indicator quantitative system of energy consumption and robustness to evaluate the path: The quantification of energy consumption is achieved by summing the agronomic energy equivalents of the agronomic measures implemented at each state transition node along the calculation path. Specifically, it calculates the fossil energy consumption equivalent of exogenous carbon input materials from raw material production, processing and transportation to soil application, as well as the energy consumption equivalent of chemical nitrogen fertilizers from synthesis, manufacturing, storage and distribution to field application. Different forms of carbon and nitrogen inputs are uniformly converted into the standard energy unit megajoules per acre to form the cumulative energy cost of the path.

[0044] Robustness is quantified by calculating the average distance from each discrete state point on the migration path to the boundary of the nearest state basin. First, the potential landscape analysis of the root zone carbon-nitrogen multi-stable system model is used to determine the boundary surface of each state basin in phase space. Then, the Euclidean distance from each state point on the path to the boundary of its nearest state basin is calculated, and the arithmetic mean is taken as the robustness measure of the path. The larger the average distance, the stronger the system's resistance to random environmental disturbances (such as extreme precipitation, temperature fluctuations, and pest infestations) when migrating along the path, and the higher the overall stability of the path.

[0045] S134. Employ a weighted cost function to evaluate and optimize the path in real time: The calculated energy consumption and robustness metrics are input into the cost function of the path planning algorithm. This cost function is defined as a weighted sum of energy consumption and robustness metrics, specifically expressed as: the cost equals α multiplied by the normalized energy consumption plus β multiplied by the reciprocal of the normalized robustness metrics, where α and β are risk preference weighting coefficients. Weights are set according to the risk preferences of farmland operators: conservative strategies increase the β weight to enhance risk resistance, while economical strategies increase the α weight to reduce input costs.

[0046] The algorithm iteratively searches in the phase space, evaluates each exploration path in real time through the cost function, and prioritizes the state transition sequence with the lowest cost value. It adopts a pruning strategy to remove path branches whose energy consumption exceeds the economic threshold or whose robustness measure is lower than the safety threshold, and optimizes the selection of state-action-state chain through dynamic programming.

[0047] S135. Output candidate state transition paths and perform time series quantization: The algorithm eventually converges and outputs one or more continuous paths with the lowest overall cost as candidate state transition paths. Each candidate state transition path is mathematically represented as a continuous curve of the evolution of key state variables over time. Based on the crop planting system of the target farmland (such as one crop per year, two crops per year, or multi-year growth cycle) and the crop growth and development patterns (seedling stage, jointing stage, flowering stage, grain filling stage, maturity stage), the state transition process in phase space is allocated to each growth stage, forming a time series of key state variables spanning one or more complete crop growth stages.

[0048] This time series analysis clearly defines the target ranges and rates of change for key state variables such as active organic carbon, inert organic carbon, bacterial biomass, fungal biomass, ammonium nitrogen, and nitrate nitrogen in each growth stage. Simultaneously, key state transition nodes are automatically identified on the continuous trajectory based on the system dynamics rate of change and stability abrupt change points: the system dynamics rate of change is measured by calculating the absolute value of the derivative of the key state variables with respect to time; when the rate of change exceeds a preset threshold, a rapid transition phase is identified. Stability abrupt change points are identified by monitoring the distance change between the state trajectory and the boundary of the state basin; when the trajectory approaches the critical distance to the boundary, a stability abrupt change is identified. The extreme points of the system dynamics rate of change and stability abrupt change points are used as staged control targets, defined as key state transition nodes, providing objects for subsequent inverse sensitivity analysis.

[0049] S2. Utilize a root region dynamic response model to perform inverse sensitivity analysis on candidate state transition paths, combining functional gene expression temporal characteristics to identify key regulatory windows and match combinations of leverage measures; specifically including: S21. Data loading and constraint embedding of the root region dynamic response model, the specific implementation process is as follows: S211. Receive the candidate state transition path and its key state transition node sequence output in step S13. This path defines the complete evolution trajectory from the current system state to the desired endpoint state, and the key state transition node sequence identifies the time-state coordinate points on the trajectory that serve as staged control targets. Use the key state transition node sequence defined by the above candidate state transition path and the continuous change trajectory of key state variables as input data and load them into the root region dynamic response model.

[0050] S212. Based on the root zone carbon and nitrogen multi-steady-state system model, this model further couples agronomic process constraints such as soil moisture dynamics, temperature effect, and crop root absorption rate. Specifically, it includes soil moisture threshold constraints (field holding capacity and wilting coefficient), temperature effect constraints (microbial activity regulation function based on Q10 temperature coefficient), and crop nutrient absorption kinetic constraints.

[0051] S213. The key feature of the model is that it couples functional gene expression time-series data, which comes from the time-series sampling of rhizosphere microbial metatranscriptome sequencing in step S11. Specifically, it includes the dynamic changes in the expression levels of carbon fixation genes, carbon mineralization genes, nitrogen fixation genes, nitrification genes, and denitrification genes at different growth stages of crops.

[0052] S214. The model incorporates quantitative influence functions on soil organic carbon mineralization, microbial community succession, and nitrogen cycling processes by different exogenous carbon inputs (biochar, straw return to the field, organic fertilizer) and chemical nitrogen fertilizer management measures (application amount and method of urea, ammonium nitrogen fertilizer, and nitrate nitrogen fertilizer). It describes the dynamic response of carbon and nitrogen in the root zone under the regulation of exogenous input and gene expression through a set of differential equations, ensuring that the simulation results conform to the physicochemical and biological laws of actual farmland ecosystems.

[0053] S22. Perform inverse sensitivity analysis on the root region dynamic response model and construct the sensitivity matrix. The specific implementation process includes: S221. On the simulation timeline of the root zone dynamic response model, using the crop growth period as the basic time frame, systematically traverse all potential agronomical measures application points that meet the feasibility of agricultural operations. This traversal process specifically covers the pre-sowing basal application period, topdressing periods at each growth stage, and critical water management periods, ensuring that potential regulation time points cover the key agronomic operation windows throughout the entire crop growth period.

[0054] At the same time, several key state transition nodes defined on the candidate state transition path loaded in step S21 are extracted and identified. These nodes serve as response observation points for reverse analysis, representing the phased control objectives that must be accurately achieved during the system state evolution process.

[0055] S222. In the root zone dynamic response model, at each ergonomically determined potential intervention time point, a unit-intensity virtual intervention pulse defined by a standardized combination of exogenous carbon input and chemical nitrogen fertilizer regulation is injected. This virtual intervention pulse has strict definition specifications: the exogenous carbon input portion is limited to a standardized application unit with a fixed carbon equivalent of biochar, straw return to the field, or organic fertilizer; the chemical nitrogen fertilizer regulation portion is limited to a standardized application unit with a fixed nitrogen equivalent of urea, ammonium nitrogen fertilizer, or nitrate nitrogen fertilizer, and the carbon-nitrogen ratio is maintained at a fixed ratio (e.g., equivalent to applying 100 kg of exogenous carbon equivalent per acre combined with 30 kg of nitrogen fertilizer per acre). This standardized combination serves as a standard perturbation input to the root zone carbon-nitrogen system, used to quantify the potential effect of implementing standard intervention at a specific time point.

[0056] S223. After each virtual control pulse is injected, the root zone dynamic response model is run forward with the current state of the soil organic carbon pool (including active organic carbon and inert organic carbon), microbial functional groups (including bacterial biomass and fungal biomass), and inorganic nitrogen pool (including ammonium nitrogen and nitrate nitrogen) corresponding to the pulse injection time point as the initial conditions.

[0057] The model is based on embedded agronomic process constraints (including soil moisture dynamics, temperature effects, and crop root uptake rates) and quantitative influence functions of exogenous carbon input and chemical nitrogen fertilizer regulation on soil processes. It combines this with the dynamic regulation of biochemical reaction rates using time-series data on functional gene expression to track the dynamic response trajectories of each state variable over time, simulating the complete soil-crop system evolution from the current point in time until the end of the target growth period. During the simulation, the actual simulated state reached at each key state transition node is precisely recorded, specifically including the values ​​of active organic carbon concentration, inert organic carbon concentration, bacterial biomass, fungal biomass, ammonium nitrogen concentration, and nitrate nitrogen concentration at the corresponding moment.

[0058] For each of the aforementioned critical state transition nodes, a baseline simulation is performed to calculate the change in deviation. A baseline simulation refers to running a simulation from the starting point of time to the corresponding time of the critical state transition node in the root region dynamic response model under the same initial conditions without applying any virtual control impulses (i.e., natural succession conditions), to obtain the baseline simulated arrival state of that node.

[0059] The actual simulated state reached after the injection of the virtual control pulse is compared with the baseline simulated state reached without the injection pulse. The deviations of the two states relative to the target state value of the critical state transition node (defined in step S13 as the center or boundary of the ideal critical state variable value domain specified for the desired endpoint state) are calculated. The change in deviation is quantified by calculating the difference between the Euclidean distance between the baseline state and the target state and the difference between the state after the pulse and the target state. The larger the difference, the higher the contribution of injecting the virtual control pulse at the control time point to achieving the critical state transition node, that is, the stronger the ability to drive the system state closer to the target state.

[0060] S224. Construct a sensitivity matrix characterizing the relationship between the influence intensity of the control time point and the state transition node. Establish a two-dimensional matrix data structure, using the potential control time point of each injected virtual control pulse as the row index of the matrix, and each key state transition node defined on the candidate state transition path as the column index of the matrix. Fill the corresponding row and column intersection positions of the matrix with the previously calculated deviation changes between a specific control time point and a specific key state transition node, thereby generating a complete sensitivity matrix.

[0061] The matrix's dimension is the number of potential regulatory time points multiplied by the number of critical state transition nodes. Each element in the matrix quantitatively represents the contribution of applying a standard virtual regulatory pulse (a specific combination of exogenous carbon input and chemical nitrogen fertilizer regulation) at a specific row index time point to the realization of a critical state transition node at a specific column index. Each row vector represents the global regulatory effectiveness distribution of that time point across all critical state transition nodes along the entire candidate state migration path, and each column vector represents the temporal sensitivity distribution of that critical state transition node to different regulatory time points. This sensitivity matrix serves as the core input data for screening critical regulatory windows in subsequent steps, providing a quantitative basis for calculating the aggregation impact index and benefit-cost ratio.

[0062] S23. The specific implementation process of conducting benefit-cost ratio analysis and screening key control windows based on the sensitivity matrix and control cost function includes: S231. For the sensitivity matrix generated in step S22, perform the aggregated influence index calculation: Extract the potential control time point (i.e., the virtual control pulse injection time point) corresponding to each row in the matrix, and obtain the reduction in deviation of that row for all key state transition nodes (i.e., the element values ​​of each column in that row of the sensitivity matrix, representing the contribution of applying standard intervention at that time point to each node).

[0063] The weighted sum is calculated based on the importance weight of each key state transition node in the candidate state transition path. The importance weight is determined based on the topological distance of each node from the desired endpoint state (the closer the node is to the desired endpoint state, the higher the weight is given, because it plays a decisive role in achieving the final goal) and the magnitude of change of key state variables between nodes (the larger the magnitude of change of state variables, the higher the weight is given, because it is more difficult to control and has a wider range of influence).

[0064] The aggregated impact index for each potential regulatory time point is obtained by weighted summation. This index quantitatively characterizes the comprehensive driving effect of implementing a standardized virtual regulatory pulse (a specific combination of exogenous carbon input and chemical nitrogen fertilizer regulation) on the entire candidate state migration path at that time point. The higher the index value, the greater the global regulatory value at that time point.

[0065] S232. For each potential control time point, a preset control cost function is invoked to calculate the expected comprehensive cost corresponding to applying a virtual control pulse of unit intensity at that time point. This cost function is constructed as a multi-dimensional cost accounting system, specifically including: The cost of exogenous carbon materials is calculated based on the market purchase price and carbon equivalent content of different carbon sources such as biochar, straw return to the field, and organic fertilizer. The cost adjustment of chemical nitrogen fertilizer is calculated based on the purchase price and nitrogen equivalent content of urea, ammonium nitrogen fertilizer, and nitrate nitrogen fertilizer. Application operation costs are calculated based on the mechanical operation costs and labor costs required for different application methods such as deep tillage, surface application, strip application, hole application, or fertigation. The potential environmental risk costs are assessed based on soil moisture conditions, temperature conditions, and climate forecasts at that point in time, including the risk costs of nitrogen fertilizer runoff into water bodies and the externality costs of greenhouse gas (nitrous oxide) emissions (e.g., calculated using carbon trading prices).

[0066] By combining the above cost dimensions, the expected comprehensive cost for each potential intervention point is calculated. This cost value reflects the economic and resource costs of implementing standard interventions at that specific point in time.

[0067] S234. For each critical state transition node, iterate through all potential control time points and calculate the benefit-cost ratio for each time point: The benefit value is taken from the element value at the intersection of the corresponding row and column in the sensitivity matrix (i.e., the reduction in deviation of the key state transition node at that time point, which characterizes the control benefit), and the cost value is taken from the expected comprehensive cost at that time point calculated above. The benefit-cost ratio is obtained by dividing the two.

[0068] For each critical state transition node, among all its corresponding potential control time points, one or more time points with the highest benefit-cost ratio are selected (if multiple time points have similar benefit-cost ratios and are all significantly higher than other time points, they are selected simultaneously to provide control flexibility). These are identified as one or more critical control windows serving the realization of that critical state transition node. The critical control window is defined as the optimal timing for applying agronomic measures within a specific crop growth period to drive the system state to the target value of that critical state transition node with the highest efficiency and cost-effectiveness.

[0069] S235. Structure and integrate all the selected key control windows: The key control windows are sorted according to the time sequence of the key state transition nodes associated with each key control window on the candidate state migration path (i.e., the crop growth period process), forming a time-series control window sequence. At the same time, priority is marked according to the benefit-cost ratio value calculated for each key control window (e.g., high benefit-cost ratio windows are marked as preferred recommendations, and medium-value windows are as alternatives).

[0070] The final output is a structured list of key regulatory windows. This list clearly defines the unique identifier, time frame (specific crop growth period and date range), associated key state transition node identifiers, benefit-cost ratio, and recommended standardized virtual regulatory pulse type (a specific combination of exogenous carbon input and chemical nitrogen fertilizer regulation) for each key regulatory window. This list serves as the temporal coordinate input for subsequent steps in the correlation analysis of leverage action combinations, ensuring that the matching of subsequent measures has a clear temporal target.

[0071] S24. Based on the key control window list, query the measure-process response knowledge graph and construct a combination of leverage measure types. The specific implementation process includes: S241. Receive the list of key control windows output in step S23, and extract the temporal positioning and target features of each key control window: The time positioning specifically refers to the crop growth stage in which the window is located, such as the basal application period before sowing, the topdressing period during the jointing stage, the critical nutrient period during flowering, or the nutrient maintenance period during the grain filling stage.

[0072] The target features refer to the direction and magnitude of the key state variables required for the next key state transition node associated with the window. Specifically, these include the percentage or absolute amount that the active organic carbon concentration needs to be increased, the target ratio that needs to be optimized for the ratio of bacterial biomass to fungal biomass, the dynamic equilibrium range that needs to be maintained for the concentrations of ammonium nitrogen and nitrate nitrogen, or the target threshold that needs to be suppressed for excessive accumulation of nitrate nitrogen.

[0073] The extracted time features and target features are converted into query conditions for a knowledge graph, which are then used as input parameters for subsequent matching algorithms.

[0074] S242. Constructing the data structure for a pre-built measure-process response knowledge graph: This knowledge graph is constructed based on long-term field test data, research results on soil biogeochemical mechanisms, and process model calibration data, and contains two core elements: nodes and edges.

[0075] The nodes are divided into two categories: The first category is exogenous carbon input nodes, which store the directional regulatory efficacy parameters of biochar (including specific pyrolysis temperature range and particle size distribution specifications), straw return to the field (including specific crushing particle size and degree of decomposition), organic fertilizer (including specific carbon-nitrogen ratio and degree of humification) and their specific application methods (including deep application to the root dense layer, surface application and mulching, specific application ratio with chemical nitrogen fertilizer and stratified application strategy) on specific biogeochemical processes in the root zone; The second category is chemical nitrogen fertilizer regulation nodes, which store the regulatory efficacy parameters of specific ratios of ammonium nitrogen and nitrate nitrogen (such as 7:3 or 5:5), the ratio of slow-release nitrogen fertilizer to fast-acting nitrogen fertilizer, the timing interval of multiple applications, etc., on the dynamics of inorganic nitrogen pool and microbial nitrogen competition.

[0076] The system connects exogenous carbon input nodes and chemical nitrogen fertilizer regulation nodes, and stores parameters of the synergistic regulatory effectiveness of specific combinations on target processes, including quantitative indicators such as the humification process rate promotion coefficient, nitrification inhibition percentage, microbial carbon pump activation coefficient, organic matter-mineral complex formation efficiency, soil aggregate stabilization degree, the distribution ratio of microbial nitrogen fixation and crop root nitrogen absorption, and nitrate nitrogen leaching risk coefficient.

[0077] S243. Perform graph traversal and matching algorithms in a knowledge graph: For each key control window in the list of key control windows, based on the specific change targets of the key state variables (active organic carbon, inert organic carbon, bacterial biomass, fungal biomass, ammonium nitrogen, nitrate nitrogen) required for the key state transition nodes associated with that window, the knowledge graph is searched for node combinations with the following characteristics: capable of driving the target biogeochemical process with the highest directional control efficiency, capable of achieving effective evolution of the system state from the current key state transition node to the next key state transition node within a specific time window, and capable of minimizing the deviation from the target state value.

[0078] The contribution scores of different combinations of measures to the target biogeochemical processes are calculated using a quantitative matching algorithm. These scores are based on the fit between regulatory efficiency parameters stored in the knowledge graph and the target features within the current window. Simultaneously, the predictive efficiency of this combination in driving state transitions within a specific key regulatory window is calculated. This predictive efficiency is obtained through rapid simulation and verification using a root zone dynamic response model with knowledge graph parameters input. The specific ratio and application method of one or more exogenous carbon inputs and chemical nitrogen fertilizer regulation that best combine the contribution scores and predictive efficiency are selected as the matching results for that window.

[0079] S244. Parametrically encapsulate the matching results to form a combination of leverage measures and output it. The encapsulation includes: the type of exogenous carbon input substance, specifying the preparation temperature range of biochar (e.g., 350-500 degrees Celsius) and particle size distribution (e.g., 0.5-2 mm), the carbon-to-nitrogen ratio of organic fertilizer (e.g., 15-25) and the degree of humification; carbon equivalent input, measured per acre or hectare (e.g., kilograms of carbon per hectare or tons of carbon per hectare), and calculating its proportion of the total carbon input over the entire growth period; application method, specifying the detailed operating procedures and supporting machinery types for strip application, hole application, full-layer mixed application, or surface covering; and application depth, precisely measured in centimeters from the ground surface. (e.g., 10 to 15 cm or 20 to 30 cm), corresponding to the root-dense layer or specific soil layer; the ratio of nitrogen forms regulated by chemical nitrogen fertilizer, specifying the precise ratio of ammonium nitrogen to nitrate nitrogen (e.g., 7:3 or 5:5); the percentage of slow-release nitrogen in total nitrogen (e.g., 40% or 60%) and the expected release period (e.g., 60 days or 90 days); the amount of chemical nitrogen fertilizer applied, the precise amount measured in kilograms of nitrogen per hectare; and the timing of application, the number of days before or after the center of the key control window (e.g., three days before or two days after), to match the peak of crop nutrient demand.

[0080] The above parameters are encapsulated in a structured manner to form a combination of leverage measures that serve the key control window. This combination serves as the core measure parameter input for the subsequent multi-objective optimization model, ensuring that the control measures are accurately matched with the time targeting and state transition objectives of the key control window, thus achieving a complete process from time window identification to specific agronomic measure configuration.

[0081] S3. Integrate candidate state transition pathways, key regulatory windows, and combinations of leverage measures to optimize for multiple objectives, including functional gene expression synergy, and generate a set of synergistic regulatory schemes; specifically including: S31. The pairing and integration process of candidate state transition paths, key control windows, and leverage measure types is implemented as follows: S311. Receive the candidate state transition paths generated in steps S13 to S14. The path includes the sequence of key state transition nodes and the continuous change trajectory of key state variables with the crop growth period. At the same time, receive the list of key control windows corresponding to each candidate state transition path and the combination of leverage measures types associated with them, which are output in steps S23 to S24.

[0082] S312. Implementation of Path and Regulation Pairing and Integration: For each candidate state transition path, extract the key state transition nodes distributed over its time series, match them with the key regulation windows determined in step S23 that serve each node, and associate them with the combination of leverage measures implemented within each window as determined in step S24, including the specific ratio and application parameters of exogenous carbon input type and chemical nitrogen fertilizer regulation. Through this pairing operation, several complete path-regulation decision chains are formed. Each decision chain fully describes the entire process from the current system state to the desired endpoint state via a specific migration path, and the implementation of specific leverage measures through specific regulation windows at key nodes of the path.

[0083] S313. The resulting decision chain is loaded as structured input data into the subsequent multi-objective optimization model, providing a basic framework for defining decision variables and constructing constraints, and ensuring that the optimization process can simultaneously consider the synergistic effects of path selection, timing, and measure configuration.

[0084] S32. The construction process of the multi-objective optimization model is as follows: S321. Construct a decision variable system for a multi-objective optimization model, and adopt a hierarchical coding structure to fully represent all elements of the control scheme.

[0085] The first level is the candidate state transition path selection variable, which uses discrete integer encoding to select a specific candidate state transition path from several path-regulation decision chains formed in step S31. This selection determines the trajectory of key state variables (including the continuous change sequence of active organic carbon, inert organic carbon, bacterial biomass, fungal biomass, ammonium nitrogen, and nitrate nitrogen) and the spatiotemporal distribution of key state transition nodes as the system evolves from the current state to the target steady state.

[0086] The second level consists of a set of variables for leveraging measures. For each key control window on the selected path, a set of specific decision sub-variables are defined: an exogenous carbon input type selection variable, using integer codes to identify biochar, straw return to the field, organic fertilizer, or their specific ratio combinations; a continuous variable for exogenous carbon equivalent input, using real numbers and taking values ​​within preset resource constraints, with the unit of measurement being kilograms of carbon per hectare; a discrete variable for application method, using integer codes to identify specific operational modes such as deep application, surface application, strip application, hole application, or fertigation; a continuous variable for application depth, using real numbers to represent the number of centimeters from the ground surface, precisely corresponding to the root-dense layer or specific soil layer; a continuous variable for chemical nitrogen fertilizer adjustment ratio, using real numbers to define the ratio coefficient of ammonium nitrogen to nitrate nitrogen (e.g., 0.7:0.3 or 0.5:0.5) and the percentage of slow-release nitrogen in total nitrogen; and a continuous variable for chemical nitrogen fertilizer application rate, using real numbers and the unit of measurement being kilograms of nitrogen per hectare. These decision variables together constitute a high-dimensional optimization search space, whose feasible region is strictly limited by subsequent constraints.

[0087] S322. Construct four mutually constrained objective functions that need to be optimized simultaneously.

[0088] The first objective function is to minimize the total economic cost, which is obtained by summing the direct and indirect costs of implementing leverage measures at each key regulatory window throughout the entire growth period. Direct costs include the cost of purchasing exogenous carbon materials (calculated based on the market prices and carbon equivalent content of biochar, straw, and organic fertilizer), the cost of purchasing chemical nitrogen fertilizer (calculated based on the nitrogen equivalent price and application rate of urea, ammonium nitrogen fertilizer, and nitrate nitrogen fertilizer), and the cost of agronomic machinery and labor (calculated based on the application method and operating area). Indirect costs include the energy consumption cost calculated based on the agronomic energy equivalent calculation method in step S13 (the economic value after converting the fossil energy consumption equivalent of exogenous carbon and chemical fertilizer throughout the entire life cycle into megajoules per acre) and the potential nitrogen leaching and greenhouse gas emission externality costs calculated based on the environmental risk assessment method in step S23 (assessed based on soil moisture conditions, climate predictions, and carbon trading prices).

[0089] The second objective function is to minimize the path tracking error, defined as the degree of deviation between the trajectory of key state variables and the ideal trajectory defined by the selected candidate state transition path in phase space during the actual control process. Specifically, the decision variables are input into the root zone dynamic response model in step S21 for forward simulation. The weighted sum of squared Euclidean distances between the simulated state (concentration of active organic carbon, inert organic carbon, bacterial biomass, fungal biomass, ammonium nitrogen, and nitrate nitrogen) at each key state transition node and the preset target state of the selected path is calculated. The weights are set according to the relative importance of each state variable to the synergistic goal of carbon sequestration and production enhancement (e.g., inert organic carbon, which contributes more to the carbon sequestration goal, is given a higher weight, and mineral nitrogen, which contributes more to the production enhancement goal, is given a higher weight).

[0090] The third objective function is to maximize the robustness to meteorological risks. It quantifies the system stability of the scheme under preset extreme meteorological scenarios (such as soil moisture falling below the wilting coefficient due to continuous drought, or rapid nitrogen leaching due to rainstorm events). The probability of key state variables remaining within the boundary of the desired state basin (from step S12) under perturbation conditions is evaluated by Monte Carlo simulation or scenario analysis, or the average distance between the system state and the boundary of the nearest state basin is calculated (refer to the robustness measurement method in step S13). The greater the distance, the higher the value of the robustness objective function.

[0091] The fourth objective function is to maximize the synergistic expression of carbon and nitrogen turnover functional genes. Its value is calculated based on the rhizosphere microbial functional gene expression profile obtained by the rhizosphere dynamic response model simulation. Specifically, it assesses the degree of enhancement of carbon fixation genes (such as cbbL, cbbM) and nitrogen fixation genes (such as nifH), and the degree of inhibition of carbon mineralization genes (such as amyA, lip) and denitrification genes (such as nirK, nosZ). By calculating the synergistic coefficients of carbon and nitrogen metabolism functional gene expression (such as the positive correlation between carbon fixation gene expression and nitrogen fixation gene expression, and the negative correlation between carbon mineralization gene expression and denitrification gene expression), it comprehensively reflects the degree of synergistic optimization of microbial carbon and nitrogen metabolism processes, ensuring that regulatory measures not only achieve carbon sequestration and increased production at the ecosystem level, but also promote the synergistic gain of carbon accumulation and efficient nitrogen utilization at the molecular biology level.

[0092] S323. Set constraints for the multi-objective optimization model to ensure the feasibility and safety of the solution.

[0093] The first constraint is the feasibility constraint of the timing of agronomic operations, which requires that the implementation time of all leverage measures must strictly fall within the time interval defined by the corresponding key control window (such as the basal application period before sowing, the topdressing period during the jointing stage, and the key nutrient period during the flowering stage), and must conform to the agronomic operation window of the crop growth period (such as avoiding mechanical operations during the sensitive period of the crop).

[0094] The second constraint is the upper limit of a single application of fertilizer, which is set based on the crop's fertilizer tolerance threshold and environmental safety standards: the upper limit of a single application of exogenous carbon prevents excessive soil solution osmotic pressure from causing physiological drought in crops, and the upper limit of a single application of chemical nitrogen fertilizer prevents the risk of seedling burn and acute leaching.

[0095] The third constraint is the resource limit for total external carbon input. Based on the total amount of carbon sources available from farmland, the upper limit of the economic budget, or the energy consumption threshold in the cost function of step S13, an upper limit for the cumulative external carbon equivalent input throughout the entire growth period is set.

[0096] The fourth type of constraint is the implicit constraint of physical feasibility, including that the nitrogen form ratio of chemical nitrogen fertilizer adjustment must be within the range of zero to one, the application amount of each type must be a non-negative real number, and the application depth must be within the physical range of the soil profile.

[0097] By rigorously formalizing the decision variables, objective functions, and constraints mentioned above, a complete multi-objective optimization model is constructed that can simultaneously optimize economic cost, path accuracy, environmental robustness, and molecular synergy, providing a standardized mathematical framework for solving the algorithm in step S33.

[0098] S33. The specific implementation process of applying a multi-objective optimization algorithm to solve the model and output the Pareto optimal solution set includes: S331. The Non-Dominated Sorting Genetic Algorithm (NSGA-II) is used as the multi-objective optimization algorithm, and algorithm initialization and chromosome encoding are performed. The path-regulation decision chain integrated in step S31 and the decision variables defined in step S32 are encoded into chromosomes to construct a multi-layered gene structure: The first layer is the gene position for selecting candidate state transition paths, which uses integer encoding to identify the specific path number selected from multiple candidate state transition paths; The second layer configures gene segments for key regulatory windows. For each key regulatory window on the selected path, the following gene sites are encoded sequentially: exogenous carbon input type gene site (using integer encoding to identify biochar, straw return to the field, organic fertilizer or their specific ratio combination), exogenous carbon equivalent input gene site (using real number encoding, the value range is based on the resource constraints defined in step S32, the unit is kilograms of carbon per hectare), application method gene site (using integer encoding to identify deep application, surface application, strip application, hole application or fertigation), application depth gene site (real number encoding, the unit is centimeters from the ground surface), chemical nitrogen fertilizer ratio regulation gene site (real number encoding, defining the ratio coefficient of ammonium nitrogen to nitrate nitrogen and the percentage of slow-release nitrogen in total nitrogen), and chemical nitrogen fertilizer application amount gene site (real number encoding, the unit is kilograms of nitrogen per hectare).

[0099] Based on the above coding rules, an initial population is randomly generated. The population size is set to be ten to twenty times the dimension of the decision variables according to the complexity of the problem. A feasibility repair mechanism (such as correcting gene values ​​that exceed the constraint range and removing individuals that do not meet the feasibility constraints of the time period) is used to ensure that each individual meets the constraint conditions defined in step S32, thus forming an initial feasible solution set.

[0100] S332. Perform fitness evaluation and evolutionary operations in each generation of evolutionary iteration. Perform multi-objective fitness evaluation on each individual in the population (i.e., a complete decision-making scheme, including the selection of specific candidate state migration paths and the configuration of leverage measures for each key regulatory window): The cost accounting module is invoked to calculate the total economic cost based on the first objective function in step S32; Using the root zone dynamic response model of step S21 with the lever measures configuration encoded by the individual as input, perform forward simulation and calculate the path tracking error (weighted sum of squared Euclidean distance) of the second objective function in step S32. The system simulates extreme weather risk scenarios (extreme drought, rainstorm events) using Monte Carlo simulation, and evaluates the robustness of the weather risk (average distance of the system state from the state basin boundary or the probability of maintenance) based on the third objective function in step S32. Based on the rhizosphere microbial functional gene expression profile output by the root zone dynamic response model, the synergy of carbon and nitrogen turnover functional gene expression (synergy coefficient of carbon and nitrogen metabolism functional gene expression) is calculated according to the fourth objective function in step S32.

[0101] Based on the evaluation values ​​of the above four objectives, the population is non-dominated and sorted, and the solution set is divided into different non-dominated levels (Rank 1 is the optimal frontier, Rank 2 is next, and so on). The crowding distance of each solution within the same level is calculated to maintain the diversity of solutions.

[0102] Based on this, the selection operation is performed, using a tournament selection mechanism to prioritize individuals with low non-dominated levels and high crowding; the crossover operation is performed, using simulated binary crossover (SBX) to perform arithmetic recombination on real-number encoded gene positions; and the mutation operation is performed, using polynomial mutation to apply random perturbations to gene values, generate offspring populations, and iteratively evolve to approach the Pareto optimal frontier.

[0103] S333. Perform convergence judgment and Pareto optimal solution set extraction: The convergence condition is set as follows: the change in the Pareto front hypervolume index is less than a preset threshold (e.g., 0.1 percent) for several consecutive generations (e.g., ten to twenty generations), or the preset maximum number of generations (e.g., two hundred to five hundred generations) is reached.

[0104] When any convergence condition is met, the evolutionary process terminates, and all individuals with a non-dominated level of 1 in the current population are extracted to form a Pareto optimal solution set. This solution set contains multiple non-dominated solutions that achieve optimal trade-offs across four dimensions: total economic cost, path tracking accuracy, weather risk robustness, and functional gene expression synergy.

[0105] S34. The process of decoding the Pareto optimal solution set and generating a set of coordinated control schemes is implemented as follows: S341. Perform chromosome decoding on each solution in the Pareto optimal solution set extracted in step S33: The candidate state migration path selection gene loci are restored to specific candidate path numbers and descriptions, including the number of crop growth stages traversed and the spatiotemporal distribution characteristics of key state transition nodes. The exogenous carbon input type gene loci corresponding to each key regulatory window are restored to specific material types (e.g., biochar prepared at a specific pyrolysis temperature, crop straw with a specific particle size, organic fertilizer formulation with a specific degree of decomposition). The exogenous carbon equivalent input gene loci are restored to precise values. The application method gene loci are restored to specific operations (surface application, furrow application, hole application, or fertigation) and application depth. The chemical nitrogen fertilizer ratio adjustment gene loci are restored to the precise ratio of ammonium nitrogen to nitrate nitrogen (e.g., 7:3 or 5:5) and the percentage of slow-release nitrogen in total nitrogen (e.g., 40% or 60%). The chemical nitrogen fertilizer application rate gene loci are restored to the precise application rate and the fine-tuning amount of application timing relative to the central time point of the key regulatory window. The decoded parameter set forms specific decision parameters that can directly guide field operations.

[0106] S342. For each Pareto optimal solution after decoding, generate a structured, executable coordinated control scheme. Each scheme is explicitly listed in the form of a document or data interface: The selected candidate state transition path is numbered and described in detail, including the continuous change trajectory of the key state variables defined by the path, the reproductive period location of the key state transition nodes, and the target value range. For each critical control window of the selected path, clearly mark the crop growth period and specific date range, and the specific leverage measures to be taken, including the precise type of exogenous carbon (specify the preparation temperature range and particle size distribution of biochar, the carbon-nitrogen ratio and humification level of organic fertilizer), the application rate per unit area (accurate to kilograms per hectare or tons per hectare, and calculate the proportion of total carbon input during the entire growth period), the application method (specify the specific operating procedures for strip application, hole application, full-layer mixed application or surface covering), the application depth (precise centimeters from the ground surface, corresponding to the root dense layer or specific soil layer), and the precise chemical nitrogen fertilizer regulation strategy to be coordinated with it (specify the reduction ratio of urea or other nitrogen fertilizer varieties, the precise ratio of ammonium nitrogen to nitrate nitrogen, the ratio of slow-release nitrogen fertilizer to fast-acting nitrogen fertilizer, and the number of days before or after the topdressing time relative to the center time point of the critical control window).

[0107] The final set of collaborative regulation schemes is presented in the form of Pareto fronts, providing farm managers with multiple options to weigh optimization across four optimization objective dimensions. Each scheme meets the constraints of agronomic feasibility, resource limitations, and environmental safety, and can be directly used for digital twin simulation and field pilot implementation in step S4.

[0108] S4. Construct a digital twin based on the set of coordinated control schemes for simulation, and use the simulation and measured data to evolve the core model and feed it back to the previous steps; specifically including: S41. Constructing a high-fidelity digital twin for the preferred solution, the specific implementation process includes: S411. Select at least one preferred scheme from the set of executable synergistic regulation schemes generated in step S34 as the construction benchmark. The selection is based on the distribution characteristics of the Pareto front solutions or the decision-maker's preference, such as prioritizing the scheme with the highest robustness or the best synergistic expression of functional genes under the premise of controllable cost. The preferred scheme clearly includes the selected specific candidate state migration path (defining the continuous change trajectory of key state variables and the sequence of key state transition nodes), the time location of each key regulation window (clarifying the timing of regulation in each growth stage), and detailed leverage measures configuration (including the type of exogenous carbon input, carbon equivalent input, application method and depth, and the ratio and amount of nitrogen forms regulated by chemical nitrogen fertilizer).

[0109] Simultaneously, the system acquires field meteorological sequences (including driving variables such as daily or hourly precipitation, air temperature, solar radiation, wind speed, and relative humidity) provided by real-time monitoring or numerical weather forecasting, as well as initial soil conditions (including the initial spatial distribution of initial active organic carbon, inert organic carbon, bacterial biomass, fungal biomass, ammonium nitrogen, and nitrate nitrogen in the root zone profile, and physical properties such as initial soil moisture content, soil temperature field, and soil texture profile) extracted through field sampling or historical databases.

[0110] S412. Based on the above inputs, the parameters are refined and modules are enhanced using the root zone dynamic response model from step S21 as the basic framework. Regarding parameter refinement, based on the soil physical properties (texture, bulk density, porosity) and historical management records (cumulative effects of exogenous carbon input and chemical nitrogen fertilizer application over the years) of the target farmland, key parameters in the model are localized and calibrated. This includes adjusting the temperature sensitivity coefficient (Q10 value) of organic matter decomposition to match local climate conditions, correcting microbial carbon use efficiency to reflect the biological activity of specific soil types, calibrating nitrogen transformation kinetic parameters (such as nitrification rate constant and denitrification loss coefficient) to adapt to soil pH and texture characteristics, and calibrating the regulatory coefficients in the functional gene expression-biochemical reaction rate mapping relationship.

[0111] S413. Improve the accuracy of model physical process representation by coupling three high-resolution sub-modules: The water transport submodule uses Richards equations or simplified numerical schemes based on soil moisture characteristic curves to replace the lumped water balance calculation in the original model, realizing the dynamic simulation of soil moisture in the root zone profile. The spatial resolution is refined to the centimeter level of soil layers (such as one to five centimeter intervals), and the temporal resolution is refined to the hour level, accurately tracking the redistribution of soil moisture driven by irrigation, precipitation and crop evapotranspiration. The heat conduction submodule is coupled with the soil heat conduction equation to simulate the diurnal and seasonal dynamics of soil temperature profiles, including the effects of freeze-thaw processes on microbial community activity and organic carbon stability, and is bidirectionally coupled with the water transport submodule. Based on the original total parameters of microbial functional clusters, the microbial community dynamics submodule explicitly distinguishes the independent evolutionary trajectories of bacterial and fungal biomass, simulates their differentiated responses to exogenous carbon input and chemical nitrogen fertilizer regulation, and introduces the community structure index (fungus / bacteria ratio) as a state variable to dynamically adjust the microbial driving coefficient of carbon and nitrogen transformation process.

[0112] S414. Through the localization of parameters and the resolution of modules described above, a high-fidelity digital twin that highly mirrors the carbon and nitrogen system in the root zone of the target farmland is constructed. This twin can reproduce the evolution of active organic carbon, inert organic carbon, bacterial biomass, fungal biomass, ammonium nitrogen, nitrate nitrogen, and the dynamics of functional gene expression in the spatiotemporal domain of the root zone under specific meteorological driving and agronomic measures. It also has the ability to calculate the critical slowing down index in real time, providing a high-fidelity computing platform for the dynamic extrapolation simulation of the entire growth period in step S42.

[0113] S42. Drive a high-fidelity digital twin to perform dynamic simulations throughout the entire reproductive period and monitor critical transition signals in real time. The specific implementation process includes: S421. Initialize the computational environment of the digital twin. Load the specific candidate state migration paths (including key state variable target trajectories and key state transition node sequences) explicitly selected in the preferred scheme as the reference baseline trajectory of the digital twin; load the configuration of each key control window defined in the preferred scheme and its associated detailed leverage measures (including the specific types of exogenous carbon input, carbon equivalent input, application method and depth, and nitrogen form ratio, application amount and application timing of chemical nitrogen fertilizer regulation) as the external perturbation event sequence of the digital twin; load the real-time meteorological sequence and soil initial conditions as driving variables and initial states.

[0114] S422. Construct three types of simulation scenarios to evaluate the performance of the solution: The baseline climate scenario uses real-time acquired or historically averaged field weather sequences, representing typical climate year types. The first preset extreme weather scenario is a continuous drought scenario (such as 30 consecutive days of precipitation less than 50% of the evapotranspiration during the growing season and soil moisture continuously less than 40% of the field capacity) to test the resistance of the scheme to water stress. The second preset extreme weather scenario is set as an extreme precipitation scenario (such as a single day's precipitation exceeding 100 mm or a cumulative precipitation exceeding 200 mm over three consecutive days) to test the stability of the scheme under wet conditions.

[0115] S423. During the dynamic simulation execution phase of the entire growth period, the digital twin iteratively calculates from the initial conditions of the crop sowing period, with a time step of hours or days. During normal periods, the natural evolution of the soil organic carbon pool, microbial functional groups, and inorganic nitrogen pool is driven by the baseline climate scenario. When the simulation time enters the key control window defined by the preferred scheme, the system is injected with the leverage measures configured for that window at the corresponding time node. These include the instantaneous increase of active organic carbon and inert organic carbon and the shift of microbial community structure caused by specific types of exogenous carbon input, as well as the changes in ammonium nitrogen and nitrate nitrogen concentrations and the activation of subsequent nitrification-denitrification processes caused by chemical nitrogen fertilizer regulation. Through coupled sub-modules of water transport, heat conduction, and microbial community dynamics, the physical, chemical, and biological chain reactions triggered by each measure in the root zone profile are calculated in real time.

[0116] S424. During the simulation, monitor the system's critical slowing-down index in real time to assess the robustness of the scheme. Based on the time series of state variables output by the digital twin (including key state variables of soil organic carbon pool, microbial biomass, and inorganic nitrogen pool), calculate the dynamic changes of the system's state statistical characteristics: calculate the time trend of the variance of key state variables (such as the ratio of active organic carbon to total carbon, fungal / bacterial ratio), an increase in variance indicates a decrease in system stability; calculate the autocorrelation coefficient (especially the first-order lag autocorrelation), an increase in autocorrelation indicates a slowdown in the system's recovery rate to disturbances; calculate the recovery rate after the system state deviates from the equilibrium position, a significant decrease in the recovery rate indicates that the system is approaching a critical transition.

[0117] When an increase in variance exceeding a preset threshold, an increase in autocorrelation exceeding a preset threshold, or a decrease in recovery rate exceeding a preset threshold are detected, the system is determined to be approaching a critical state transition point, triggering an early warning signal. The digital twin dynamically adjusts the simulation parameters based on the early warning signal, simulating the system's critical response behavior under different disturbance intensities, and evaluating the robustness of the control scheme in the critical region.

[0118] S425. The simulation continues until crop harvest, outputting a complete prediction dataset after the optimal solution is implemented: the root zone state evolution trajectory (including the dynamic profile distribution of each state variable at each time point throughout the growth period), the deviation of the predicted values ​​of state variables at key state transition nodes from the target trajectory, the final carbon sequestration (the difference in soil organic carbon capacity between harvest and sowing periods), and crop yield (grain yield based on the coupled simulation of root zone nitrogen supply dynamics and crop growth model). Simultaneously, it outputs the critical transition warning signals detected during the simulation, along with their timing and severity assessments, providing a predictive benchmark for model validation and calibration in step S43.

[0119] S43. The specific implementation process of calculating model prediction bias and evolving the core model using deduction and measured data includes: S431. Establish limited pilot monitoring plots in the target farmland where the optimized scheme is implemented. Based on soil spatial variability characteristics (including soil texture profile, topographic relief, historical fertility levels, and microclimate differences), stratified random sampling is used to establish three to five replicate plots, each with an area of ​​no less than 30 square meters. Systematic root zone status monitoring data are collected before and after the start of each key regulatory window defined by the optimized scheme, as well as during key crop growth stages (such as jointing, flowering, and grain-filling stages). The collected actual root zone status monitoring data includes: soil soluble organic carbon (dynamic indicator of active organic carbon pool), microbial biomass carbon (microbial functional group biomass indicator), and profile data of ammonium nitrogen and nitrate nitrogen content (usually sampled in stratified layers of 0-10 cm, 10-20 cm, and 20-40 cm). Simultaneously, rhizosphere soil samples are collected for metagenomic sequencing to obtain expression levels of carbon fixation genes, carbon mineralization genes, nitrogen fixation genes, nitrification genes, and denitrification genes, forming a high-temporal-resolution measured dataset that strictly corresponds to the output variables of the digital twin.

[0120] S432. The predicted root zone state evolution trajectory and functional gene expression dynamics, obtained from high-fidelity digital twin simulations corresponding to the spatial location and soil depth of the pilot plot, are precisely aligned with actual monitoring data in the time dimension, matching them to the same monitoring date, the same soil depth, and the same growth stage. For each matched time point, each key state variable (soluble organic carbon, microbial biomass carbon, ammonium nitrogen, nitrate nitrogen), and the expression level of each functional gene, the absolute and relative differences between the predicted and measured values ​​are calculated.

[0121] We use the weighted root mean square error as a comprehensive quantitative index of model prediction bias. We assign weights to each key state variable (based on the relative importance of each variable to the synergistic goal of carbon sequestration and production enhancement and the uncertainty of the observation data, such as assigning a lower weight to microbial biomass carbon with large observational variability and a higher weight to the key nitrate nitrogen leaching risk index). We calculate the square mean of the difference between the weighted predicted value and the measured value, and then take the square root to obtain the overall bias index characterizing the prediction accuracy of the digital twin.

[0122] S433. Automated parameter calibration of the core model is performed using an ensemble Kalman filter assimilation algorithm. With minimizing the prediction bias of the above models as the optimization objective, the core parameters to be optimized in the root zone carbon and nitrogen multi-stable system model of step S12 and the root zone dynamic response model of step S21 are defined as vectors to be estimated. These include the organic matter decomposition rate constant, microbial carbon use efficiency, nitrification activation energy, denitrification loss coefficient, stability allocation coefficient of exogenous carbon input, and the regulatory coefficient of functional gene expression on the biochemical reaction rate.

[0123] S434. Data assimilation is performed using an ensemble Kalman filter algorithm: First, an initial set of parameters to be estimated is generated, and the possible distribution range of each parameter is set based on prior knowledge; the state of each set member is predicted using a root region dynamic response model; the collected actual root region state monitoring data and macrotranscriptome data are used as observation constraints to calculate the covariance matrix between the prediction set and the observations; the weights of set members and parameter estimates are updated using the Kalman gain matrix; the prediction-update loop is iteratively executed to gradually reduce the range of parameter uncertainty until the model prediction bias converges to below the preset threshold or the parameter update amount is less than the preset tolerance, thus completing the automated calibration of the core parameters.

[0124] S435. Model structure optimization based on residual analysis after parameter calibration. Analyze the spatiotemporal distribution patterns of systematic prediction biases that still exist after calibration, and identify uncharacterized key processes or mechanism defects: If residual analysis shows systematic biases under specific moisture conditions (such as extreme humidity or drought) (e.g., persistent underestimation of nitrate nitrogen concentration under extreme humidity), indicating insufficient characterization of the denitrification process, optimize the coupling relationship between the water transport submodule and the nitrogen transformation module in the rhizosphere dynamic response model, and add, delete, or modify the functional relationship regarding the denitrification process under anaerobic conditions (e.g., introducing a dynamic denitrification rate equation based on soil redox potential); if persistent biases are found in microbial biomass carbon at specific growth stages, adjust the carbon allocation function for crop root exudate input in the microbial community dynamic submodule; if a systematic shift is found in the mapping between functional gene expression and biochemical reaction rate, correct the parameters or function form in the gene expression-reaction rate conversion function; if a systematic shift is found in the conversion between active organic carbon and inert organic carbon, adjust the flux connection mode between different carbon and nitrogen pools. Through the above structural optimization, the model's ability to interpret and predict actual root region state monitoring data and metatranscriptome data is improved, thus achieving the evolution of the model structure.

[0125] S44. The specific implementation process of feeding back the evolved core model to form a self-evolutionary closed loop includes: S441. The evolutionary model system, after automated parameter calibration and structural optimization in step S43, is systematically packaged and updated. The evolved model system includes two core components: First, the updated root zone carbon and nitrogen multi-steady-state system model (derived from step S12) has been corrected based on actual farmland observation data. Its core parameters (such as the threshold for determining the state basin boundary and the quantitative index of steady-state stability) and structure (such as the nonlinear feedback equation between carbon and nitrogen pools) have been modified. Second, the updated root zone dynamic response model (derived from step S21) has adapted the parameters and structure of its response functions to exogenous carbon input and chemical nitrogen fertilizer regulation, agronomic process constraints, and high-resolution sub-modules (water transport, heat conduction, and microbial community dynamics) to the specific biogeochemical characteristics of the target farmland, and the functional gene expression-biochemical reaction rate mapping relationship has been calibrated based on metatranscriptome data.

[0126] S442. Feedback the evolved root zone carbon and nitrogen multi-steady-state system model to step S1 (specifically step S11) for a new round of system state diagnosis and candidate state migration path planning. At the start of a new application cycle, the current monitoring data of the target farmland (real-time status of soil organic carbon pool, microbial functional groups, and inorganic nitrogen pool) and the updated macrotranscriptome data are input into the evolved model. Based on the calibrated parameters and optimized structure, the model more accurately diagnoses the steady-state or transient characteristics of the current system (including the precise attribution of state basins, stability level, and state basin boundary distances), and based on the updated phase space topology, plans new candidate state migration paths from the current system state to the desired endpoint state with lower energy consumption and higher robustness (step S13), while more accurately identifying key state transition nodes (step S14).

[0127] S443. The evolved root region dynamic response model is fed back to step S2 (specifically step S21) for the next round of inverse identification of key regulatory windows and matching of leverage measures. In the new cycle, the evolved model, based on the calibrated response function and structural optimization, more accurately simulates the perturbation effects of different exogenous carbon inputs and chemical nitrogen fertilizer regulation at specific key regulatory windows (including the response patterns of microbial functional gene expression). This generates a more accurate sensitivity matrix in the inverse sensitivity analysis in step S22, screens out key regulatory windows with better benefit-cost ratios in step S23, and in step S24, matches more efficient combinations of leverage measures through the updated measure-process response knowledge graph (which incorporates new mechanisms and functional gene expression data identified during model structural optimization), including more precise selection of exogenous carbon input types, determination of carbon equivalent input, and chemical nitrogen fertilizer regulation ratios.

[0128] S444. Through the aforementioned two-way feedback mechanism, a self-evolutionary closed loop is achieved, from scheme deduction – pilot verification – model evolution – a new round of path planning and regulation identification. This closed loop enables the methodology to continuously improve its cognitive accuracy (by reducing uncertainty through parameter calibration) and mechanistic characterization ability (by incorporating newly discovered processes and metatranscriptome-level evidence through structural optimization) in continuous field application and data accumulation, thereby achieving adaptive evolution and long-term performance improvement of the carbon sequestration and yield-increasing synergistic optimization method.

[0129] As can be seen from the above description, the embodiments of the present invention achieve the following technical effects: This invention represents a leap from static agronomic management to dynamic system path planning, fundamentally solving the inertia problem of complex ecosystem transformation. Traditional methods rely on fixed fertilization schemes, failing to guide the active and controllable state transition of the multi-steady-state complex system of root zone carbon and nitrogen. This invention introduces the concept of state basins and calculates the boundary distances of these basins, providing a theoretical benchmark for quantifying the current stability of farmland ecosystems and the ease with which they can transition to a target steady state. Based on this, combined with inverse sensitivity analysis, candidate state transition paths with minimal energy consumption and maximum robustness can be deduced. This is equivalent to creating a dynamic navigation map for farmland management, leading from the current state to high-yield, low-carbon goals. It decomposes macro-level objectives into a series of controllable, phased objectives (key state transition nodes), making proactive guidance and precise transformation of complex ecological processes possible. This overcomes the problems of transformation failure or high costs caused by traditional methods neglecting system inertia.

[0130] This invention has developed a precise inverse regulation mechanism based on microbial functional responses, moving away from extensive empirical timing. This mechanism enables the optimal spatiotemporal allocation of management resources. Traditional fertilization relies on fixed phenological periods and cannot respond to the real-time dynamics of soil microbial function. This invention constructs a functional gene expression-biochemical reaction rate mapping relationship, directly linking macroscopic system behavior with microscopic microbial functional gene expression, giving the model cross-scale mechanistic characterization capabilities. Furthermore, the sensitivity matrix generated through inverse sensitivity analysis can precisely locate the key regulatory windows that have the greatest impact on system state transitions, much like medical diagnosis. More importantly, by querying a pre-constructed measure-process response knowledge graph, it can match each specific window with the most efficient combination of leverage measures that drive the expression of specific microbial functional genes, thereby influencing the system state to evolve in a predetermined direction. This mechanism ensures that the most appropriate and minimal intervention is applied at the right time, achieving precise spatiotemporal synchronization and efficient coupling between management measures and soil biological processes.

[0131] This invention constructs an intelligent closed-loop system, evolving from fixed models and schemes to systems with forward-looking early warning and continuous self-evolution capabilities, ensuring the long-term effectiveness and adaptability of the technology. Existing agricultural models, once established, often become rigid and difficult to adapt to changing environments. This invention, by monitoring critical transition early warning signals in a high-fidelity digital twin, can detect early signs of system near-collapse or sudden state changes, thus enabling forward-looking stress testing and optimization of the robustness of control schemes, transforming passive response into proactive defense. More importantly, by employing a Kalman filter assimilation algorithm, it can automatically and efficiently feed back and calibrate field-measured root zone state data and high-dimensional microbial macrotranscriptome data into the core model, driving the continuous evolution of model parameters and structure. This practice-learning-evolution intelligent closed loop allows the entire system to continuously self-optimize with seasonal changes and environmental shifts, ensuring that the generated optimized schemes maintain high scientific validity and effectiveness not only in the present but also in long-term applications, fundamentally solving the industry pain point of models and schemes becoming ineffective over time.

[0132] The embodiments and / or implementation methods described above are merely preferred embodiments and / or implementation methods for implementing the technology of the present invention, and are not intended to limit the implementation methods of the technology of the present invention in any way. Any person skilled in the art can make some modifications or alterations to other equivalent embodiments without departing from the scope of the technical means disclosed in the present invention, but these should still be regarded as the technology or embodiments that are substantially the same as the present invention.

[0133] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. The above descriptions are only preferred embodiments of this application. It should be noted that due to the limitations of written expression, while there are objectively infinite specific structures, those skilled in the art can make several improvements, modifications, or changes without departing from the principles of this application, and can also combine the above technical features in an appropriate manner. These improvements, modifications, changes, or combinations, or the direct application of the inventive concept and technical solution to other situations without modification, should all be considered within the scope of protection of this application.

Claims

1. A method for synergistic optimization of carbon sequestration and yield enhancement based on the carbon and nitrogen regulation mechanism of crop root zones, characterized in that, include: Based on multi-source data, the current state is diagnosed by integrating the functional gene expression-biochemical reaction rate mapping relationship into a root region carbon and nitrogen multi-stable system model, and the comprehensive optimal candidate state migration path is planned. The candidate state transition path was analyzed by using the root region dynamic response model, and key regulatory windows were identified and combinations of leverage measures were matched by combining the temporal characteristics of functional gene expression. By integrating the candidate state transition paths, key regulatory windows, and combinations of leverage measures, optimization is performed with multiple objectives, including functional gene expression synergy, to generate a set of synergistic regulatory schemes. A high-fidelity digital twin is constructed based on the aforementioned set of collaborative control schemes for simulation. The scheme is evaluated by monitoring critical transition early warning signals, and the results are used to evolve the core model and feed back to the preceding steps.

2. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 1, characterized in that, The process of diagnosing the current state based on monitoring data using a root-region carbon and nitrogen multi-stable-state system model and planning candidate state transition paths to the target state specifically includes: The current root zone carbon and nitrogen system status monitoring data, historical management data, and functional gene expression-biochemical reaction rate mapping relationship established based on rhizosphere microbial metatranscriptome sequencing data of the target farmland are input into the root zone carbon and nitrogen multistable system model. The root zone carbon and nitrogen multi-steady-state system model diagnoses the current system state and calculates its distance from the target steady-state state basin boundary by coupling the nonlinear dynamic equations and mapping relationships between the soil organic carbon pool, microbial functional groups and inorganic nitrogen pool. Starting from the current state and ending at the preset carbon sequestration and production enhancement synergistic target state, in the model phase space, a state-action-state chain search algorithm is used. The energy consumption quantified by agronomic energy equivalent and the robustness measured by the average distance from the state basin boundary are used as the comprehensive cost function to plan at least one candidate state transition path.

3. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 2, characterized in that, The root zone carbon and nitrogen multi-steady-state system model diagnoses the current system state and calculates its distance from the target steady-state basin boundary by coupling the nonlinear dynamic equations and mapping relationships between the soil organic carbon pool, microbial functional groups, and inorganic nitrogen pool. Specifically, this includes: The root zone carbon and nitrogen multi-steady-state system model couples the nonlinear dynamic relationship between the soil organic carbon pool, microbial functional groups and inorganic nitrogen pool through the first set of differential equations, and embeds the functional gene expression-biochemical reaction rate mapping relationship through the second set of functions, so as to dynamically transform the expression level of key functional genes into the constraint condition of corresponding biochemical reaction flux. Based on the root region carbon and nitrogen multi-stable system model, the current monitoring data is input, and the coupling equations and functions are solved by numerical integration to determine the coordinate position of the current system state in the model phase space. In the phase space, the state basin boundary where the preset target steady state is located is identified based on potential energy landscape analysis, and the shortest Euclidean distance or potential energy difference from the current state coordinate point to the state basin boundary is calculated as the quantified value of the state basin boundary distance.

4. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 3, characterized in that, The process of using a root region dynamic response model to perform inverse sensitivity analysis on the candidate state transition paths, combined with the temporal characteristics of functional gene expression, to identify key regulatory windows and match combinations of leverage measures, specifically includes: The candidate state transition path is input into the root region dynamic response model, which is embedded with agronomic process constraints and coupled with functional gene expression time series data; Inverse sensitivity analysis is performed on the root region dynamic response model. Potential control time points are traversed on the simulation timeline and virtual control pulses of unit intensity are injected to generate a sensitivity matrix characterizing the influence intensity of control time points-state transition nodes. Based on the sensitivity matrix and the preset control cost function, a benefit-cost ratio analysis is performed to obtain a list of key control windows. Based on the list of key control windows, the pre-constructed measure-process response knowledge graph is queried, and one or more specific ratios and application methods of exogenous carbon input and chemical nitrogen fertilizer regulation that can most efficiently drive the system state to the next key state transition node are matched and output, forming a combination of leverage measure types.

5. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 4, characterized in that, Inverse sensitivity analysis is performed on the root region dynamic response model. Potential control time points are traversed along the simulation timeline, and virtual control pulses of unit intensity are injected to generate a sensitivity matrix characterizing the influence intensity of the control time point-state transition node. Specifically, this includes: On the timeline of the root region dynamic response model, all potential control time points that meet the feasibility of agricultural operations are traversed, as well as several key state transition nodes defined on the path. In the root zone dynamic response model, at each of the regulation time points, a virtual regulation pulse of unit intensity is injected, defined by a combination of standardized exogenous carbon input and chemical nitrogen fertilizer regulation. After each virtual control pulse is injected, the root region dynamic response model is run forward to the end of the path; the actual simulated arrival state of each key state transition node is accurately recorded and compared with the baseline arrival state when no pulse is injected, and the change in state deviation caused by each pulse to each node is calculated. Using the potential regulation time as the row index of the matrix and each key state transition node as the column index of the matrix, the calculated deviation change is filled into the corresponding position to construct a sensitivity matrix, where the matrix element values ​​characterize the influence strength of applying standard intervention at a specific time point on achieving state transition at a specific node.

6. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 5, characterized in that, The combination of types of leverage measures specifically includes: For each window in the list of key control windows, extract its time location, the associated next key state transition node, and the direction and magnitude of the corresponding target state variable change; The extracted information is used as a query request and input into a pre-constructed measure-process response knowledge graph. The knowledge graph stores different types of exogenous carbon input and chemical nitrogen fertilizer regulation methods as nodes, and stores their directional regulation efficiency parameters on soil humification, nitrification-denitrification, and microbial carbon pump activation of specific biogeochemical processes as edges. Graph traversal and matching algorithms are performed in the knowledge graph to find a combination of measures that can drive the system state from the current node to the next critical state transition node with the highest prediction efficiency. The matched one or more exogenous carbon inputs are parameterized and encapsulated with specific ratios and application methods of chemical nitrogen fertilizer regulation to form a combination of leverage measures that serve the window and are then output.

7. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 6, characterized in that, By integrating the candidate state transition paths, key regulatory windows, and combinations of leverage measures, and optimizing for multiple objectives including functional gene expression synergy, a set of synergistic regulatory schemes is generated, specifically including: The candidate state transition paths, the list of key control windows corresponding to each path, and the associated leverage measures are paired and integrated to form several path-control decision chains. A multi-objective optimization model was established, with decision variables including the selection of candidate state transition paths and the determination of the specific types and amounts of measures in the associated leverage measure combinations for each key control window; the optimization objectives included at least minimizing total economic cost, minimizing path tracking error, maximizing meteorological risk robustness, and maximizing the synergistic expression of carbon and nitrogen turnover functional genes. The multi-objective optimization model is solved using a multi-objective optimization algorithm, and the Pareto optimal solution set is output. The Pareto optimal solution set is decoded to generate a set of coordinated control schemes that include specific migration paths, control timing, types of measures, and application amounts.

8. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 7, characterized in that, The establishment of the multi-objective optimization model specifically includes: The set of decision variables includes: a first variable, used to select a specific candidate state transition path from the path-regulation decision chain; and a second set of variables, used to determine the specific exogenous carbon type, application amount, application method, and chemical nitrogen fertilizer ratio and application amount in the associated combination of leverage measures at each key regulation window of the selected path. The optimization objective functions include: a first objective function, used to calculate and minimize the total economic cost of implementing all selected measures; a second objective function, used to calculate and minimize the path tracking error; a third objective function, used to evaluate and maximize the robustness to meteorological risks; and a fourth objective function, used to evaluate and maximize the synergistic expression of carbon and nitrogen turnover functional genes. The constraints of the multi-objective optimization model include at least the agronomic operation period, the upper limit of a single fertilization amount, and the total external carbon input resource limit.

9. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 8, characterized in that, A high-fidelity digital twin is constructed based on the aforementioned coordinated regulation scheme set for simulation. The scheme is evaluated by monitoring critical transition early warning signals, and the results are used to evolve the core model and feed back to the preceding steps. Specifically, this includes: For at least one preferred scheme selected from the set of coordinated control schemes, a high-fidelity digital twin that is highly mirrored with the target farmland is constructed by combining real-time meteorological and soil data. The high-fidelity digital twin is driven to perform dynamic simulations under baseline and extreme weather scenarios, and the critical transition warning signals of the system are monitored in real time to evaluate the robustness of the scheme. The predicted data are compared with the measured root zone state and microbial macrotranscriptome data of the target farmland. The core parameters of the root zone carbon and nitrogen multi-stable system model and the root zone dynamic response model are automatically calibrated and evolved using the ensemble Kalman filter assimilation algorithm. The evolved model is then fed back to facilitate a new round of system diagnosis and path planning.

10. The method for synergistic optimization of carbon sequestration and yield increase based on the crop root zone carbon and nitrogen regulation mechanism according to claim 9, characterized in that, The core parameters of the root-zone carbon-nitrogen multi-stable-state system model and the root-zone dynamic response model are automatically calibrated and evolved, specifically including: Observational data corresponding to time and variables are extracted from the measured root zone state and microbial macrotranscriptome data of the target farmland to form an observation vector; Predictive data corresponding to the pilot monitoring time are extracted from the high-fidelity digital twin simulation results to form a model prediction vector; the difference between the observation vector and the model prediction vector is calculated. With minimizing model prediction bias as the optimization objective, the core parameters to be optimized in the root region carbon-nitrogen multi-stable system model and the root region dynamic response model are defined as vectors to be estimated. An ensemble Kalman filter assimilation algorithm is employed to predict the model state based on the current set of parameter vectors; the covariance matrix between the prediction set and the observation vectors is calculated; and the parameter set is updated using the Kalman gain formula so that the updated parameter set statistically minimizes the deviation between the prediction data and the observation data. The mean of the updated parameter set obtained after assimilation convergence is used as the new parameter values ​​after automatic calibration and evolution, and then updated into the two core models.

Images