An agricultural nitrogen and phosphorus full flux assimilation and scenario difference potential accounting method

Abstract

The application discloses a kind of agricultural nitrogen phosphorus full flux assimilation and scenario difference potential accounting method, belong to agricultural ecosystem nutrient cycle accounting and data assimilation field, including: obtaining standardized interface data package;Build assimilable agricultural nitrogen phosphorus state space model and forward prediction interface;Establish observation operator library and unified error statistics framework;Carry out information amount diagnosis and generate assimilation configuration;Build and implement material conservation constraints and physical feasible region constraint set;Constrained joint data assimilation is solved, obtains unified posteriori datum package;Based on the benchmark, build counterfactual baseline and measure scenario, carry out scenario difference potential accounting, output includes the difference potential result of uncertainty.This application solves the problem that multi-source data is difficult to fuse, full path flux is not closed, uncertainty lacks quantification and counterfactual baseline comparability is insufficient in the prior art, realizes the traceability of agricultural nitrogen phosphorus full process flux, verifiable accounting and precise attribution of emission reduction potential.

Description

Technical Field

[0001] This invention belongs to the field of agricultural ecosystem nutrient cycling accounting and data assimilation technology, and particularly relates to a method for agricultural nitrogen and phosphorus total flux assimilation and scenario differential potential accounting. Background Technology

[0002] Excessive input of nitrogen (N) and phosphorus (P) into agricultural ecosystems is a major contributing factor to complex environmental risks, including nutrient imbalances in farmland, eutrophication of water bodies, nitrate pollution of groundwater, and greenhouse gas emissions. Therefore, establishing traceable and verifiable methods for calculating nitrogen and phosphorus emission fluxes and assessing emission reduction potential is crucial for policy formulation and optimization of governance measures. Currently, the mainstream methods in this field can be summarized into six categories: ① Nutrient budget accounting and material flow analysis, which establishes a budget based on system inputs and outputs, has low data requirements, and is suitable for rapid regional accounting, but is sensitive to boundary delineation and struggles to handle "implicit fluxes" such as deposition and nitrogen fixation; ② Emission factor or empirical coefficient method, which quickly estimates loss fluxes through "activity data × factor parameters," is computationally simple, but its parameters are regionally dependent and difficult to characterize nonlinear responses; ③ Field process mechanism model, which outputs fluxes with organic constraints by coupling soil hydrothermal and biogeochemical processes. ① Quantitative analysis can be used for process attribution, but the parameters have high dimensionality and insufficient parameter transferability when applied across regions; ② Watershed hydrological and water quality models connect the entire process of "source-transport-sink" and can quantify the load of receiving water bodies, but they require high data accuracy and errors are easily amplified by accumulation; ③ Life cycle assessment method expands the system boundary and is suitable for cross-product comparison, but it is difficult to provide direct attribution for optimizing measures at the field level; ④ Data-driven remote sensing fusion and statistical inference can identify risk hotspots with the help of spatial statistics, but the consistency and verifiability of flux estimation are difficult to guarantee when there is a lack of material conservation constraints.

[0003] Although the above methods provide technical support at different levels, there are still common bottlenecks in the accounting of agricultural nitrogen and phosphorus emission reduction potential, mainly manifested as: (1) difficulty in integrating multi-source data: data such as fertilizer statistics, soil monitoring, remote sensing inversion, water quality and greenhouse gas observation have significant differences in spatiotemporal scale and error structure, lacking integration and mutual verification under a unified framework, resulting in poor consistency of accounting conclusions; (2) the flux chain of the whole path is not closed: existing studies often focus on a single loss flux or surplus indicator, which is difficult to simultaneously depict the whole process of "input-conversion-transfer-output-sink change". (3) Lack of systematic quantification of uncertainty: The superposition of errors in input data, parameters and model structure makes the robustness of cross-regional promotion insufficient. Most studies only provide point estimates, lack confidence intervals and risk probabilities, which are difficult to meet the reliability requirements of policy evaluation. (4) Insufficient comparability between counterfactual baseline and measure scenario: If the baseline and measure scenario are inconsistent in initial state, parameter set or external driving force, the difference results will be mixed with the influence of non-measure factors, which weakens the attribution ability and verifiability of emission reduction potential assessment. Summary of the Invention

[0004] The main objective of this invention is to propose a method for agricultural nitrogen and phosphorus total flux assimilation and scenario-based differential potential calculation, which aims to overcome the shortcomings of existing agricultural nitrogen and phosphorus emission accounting and emission reduction potential assessment technologies, such as the difficulty in integrating multi-source observation data under a unified framework, the lack of a closed-loop flux accounting chain and incomplete flux, the complexity of uncertainty sources and the lack of systematic quantification, and the inconsistency between the counterfactual baseline and the action scenario setting, leading to a mixture of action effects and background differences. By constructing a unified error statistics and assimilation estimation framework that integrates multi-source data such as fertilization statistics, soil monitoring, remote sensing inversion, water quality, groundwater, and greenhouse gases, and by introducing material conservation constraints and physical feasible domain constraints to achieve flux closure and result consistency, the method conducts counterfactual baseline-action scenario differential calculation under a unified posterior benchmark, thereby achieving traceability, verifiability, and attributability of emission reduction potential results.

[0005] To achieve the above objectives, this invention provides a method for agricultural nitrogen and phosphorus total flux assimilation and scenario-difference potential calculation, comprising: Based on multi-source data, obtain standardized interface data packets; Based on the standardized interface data package, an assimilable agricultural nitrogen and phosphorus state-space model and forward prediction interface are constructed. Based on the standardized interface data package and the forward prediction interface, an observation operator library and a unified error statistics framework are established. Based on the aforementioned observation operator library and unified error statistics framework, information content diagnosis is performed and assimilation configuration is generated; Based on the aforementioned forward prediction interface, observation operator library, unified error statistics framework, and assimilation configuration, a constraint set of matter conservation constraints and physical feasible region constraints is constructed and implemented. Based on the forward prediction interface, observation operator library, unified error statistics framework, assimilation configuration and constraint set, perform constraint joint data assimilation solution to obtain unified posterior baseline data package; Based on the unified posterior baseline data package, a counterfactual baseline and action scenarios are constructed, scenario differential potential is calculated, and differential potential results containing uncertainty are output.

[0006] Preferably, the process of obtaining the standardized interface data packet includes: Determine the system boundaries, accounting units, and time base of the accounting objects; Collect and organize the driving input library, observation library, and prior library, and perform element coding and unit unification; Spatiotemporal scale matching and quality control are performed on the data in the driving input library and the observation library; The data that has undergone spatiotemporal scale matching and quality control is assembled into a standardized interface data package containing observations, driving inputs, quality labels, and metadata.

[0007] Preferably, the process of constructing an assimilable agricultural nitrogen and phosphorus state-space model and forward prediction interface includes: Construct an extended state vector to describe key nitrogen and phosphorus pools, crop nutrient status, and optional deviation states; Define a full flux vector that explicitly decomposes the input, transformation, transfer, loss, and output processes into traceable components; Define the driving inputs as the entry point for management measures, the static covariates characterizing spatial heterogeneity, and the set of parameters to be estimated; Based on the extended state vector, full flux vector, driving input, static covariates, and parameter set, a forward prediction interface is established to advance the state and synchronously output the full flux at time steps.

[0008] Preferably, the process of establishing the observation operator library and the unified error statistics framework includes: Define corresponding observation mapping rules for the classification of multi-source observation data, and construct an observation operator library; For each type of observation data, an observation error model is established, and the observation error covariance matrix is ​​constructed. Establish a process error or background error model to describe the uncertainty of model structure, driving forces, and parameters; By introducing a bias state or reducing the weight of representative errors, the inconsistency between system bias and scale support is explicitly addressed.

[0009] Preferably, the process of performing information content diagnosis and generating assimilation configuration includes: The contribution of different observation sources to the constraints of state, flux and parameters is evaluated, and observability analysis is performed to obtain the observability analysis results. Based on the observability analysis results, assimilation variables were screened and stratified. Develop assimilation and update strategies in groups or phases; Regularization and flux decomposition rules are introduced for weakly identifiable flux components; The weights of soft and hard constraints are configured collaboratively and interfaced with the constraint set and error statistics framework. Perform a consistency check before assimilation and output an assimilation configuration file containing a list of variables, update order, regularization rules, and weight configuration.

[0010] Preferably, the process of constructing and implementing the constraint set includes: Construct equation constraints for the conservation of nitrogen and phosphorus and the total flux closure; Construct inequality constraints to ensure that state variables, fluxes, and key derived variables are non-negative and within the physical boundary; Construct a parameter feasible region and scale constraints to ensure that process parameters are physically feasible and spatially reasonable; Construct process coupling constraints to ensure the consistency of substrate constraints, carrier constraints, and morphological transformation in order to guarantee the mechanistic self-consistency of the process; Define a constraint implementation operator for performing feasible region correction on candidate solutions and outputting a closure verification diagnostic quantity.

[0011] Preferably, the process of performing constraint joint data assimilation solution includes: Generate an initialization set of joint assimilation variables containing states and parameters based on the prior distribution; The forward prediction interface is used to advance the set members to generate prior states and prior total flux. Based on the aforementioned observation operator library and unified error statistics framework, multi-source observations are fused to update the prior set and obtain candidate posterior solutions. The constraint set and its constraint implementation operators are invoked to correct the candidate posterior solution, thereby obtaining a feasible posterior solution that satisfies conservation closure and physical feasibility. The feasible posterior solutions are subjected to set statistics to quantify the posterior uncertainty and perform attribution decomposition. Perform consistency diagnosis and quality classification on the posterior results, and output a unified posterior benchmark data package with admission criteria.

[0012] Preferably, the process of performing scenario differential potential calculation includes: Using the unified posterior baseline data package as a common starting point, baseline scenarios and action scenarios are constructed, wherein the action scenarios are implemented by modifying the driving input; Under the same model, constraint set, and initial state, forward extrapolation is performed on the baseline scenario and the action scenario respectively, and their respective state trajectories and full flux ledgers are output; By pairing and differencing members of the same posterior set, the flux difference and state difference between the action scenario and the baseline scenario are calculated, and a potential accounting index system is constructed. The posterior uncertainty is propagated in the scenario difference, the difference results are aggregated and statistically analyzed, the potential confidence interval and robustness probability are output, and the risk screening and feasibility assessment are combined with engineering constraints. The differential potential results of the accounting units are spatially aggregated according to the specified partition boundaries, and the statistical and priority calculation results at the partition scale are output.

[0013] Preferably, the construction of the measure scenario includes: by adjusting the fertilizer application amount, fertilizer application time, fertilizer application method, fertilizer form, irrigation system, and organic material management in the driving input, nitrogen control measures, phosphorus control measures, water and transport control measures, and organic fertilizer and straw management measures with different intensity levels.

[0014] Preferably, the output of the difference potential result containing uncertainty includes: Output the posterior state and the full-process accounting ledger of total flux; Output the scenario potential list and spatial distribution product of the differential potential results; Output probabilistic expressions and risk characterization indicators for posterior uncertainty and difference potential; Output the verification indicators and log records for conservation closure and consistency; Based on feasibility assessment and robustness probability, the results of optimal measures and comprehensive evaluation are output.

[0015] Compared with the prior art, the present invention has the following advantages and technical effects: This invention addresses the aforementioned problems of difficulty in fusing multi-source data, incomplete flux chains, lack of quantification of uncertainty, and insufficient comparability of counterfactual baselines, achieving the following beneficial effects: This invention achieves the fusion of multi-source data under a unified conservation framework, improving the internal consistency of accounting conclusions: By constructing a unified error statistics framework and observation operator library, this invention jointly assimilates and estimates multi-source observation data such as fertilization statistics, soil monitoring, remote sensing inversion, water quality and greenhouse gas, and introduces material conservation constraints, enabling the flux of nitrogen and phosphorus throughout the entire process of "input-conversion-transport-output-sink change" to achieve closure and consistency verification within the same framework, significantly reducing the "gap" between fluxes and improving the physical rationality and internal consistency of accounting results.

[0016] This invention achieves full flux chain closure and flux completeness, enhancing the traceability of accounting: By constructing a state-space model containing key quantities and the full flux vector, and implementing strict conservation closure constraints, this invention ensures that each flux component within the system is traceable, aggregated, and closed. Simultaneously, by recording constraint residuals and error convergence information, the entire accounting chain possesses traceability, reproducibility, and verifiability, facilitating third-party auditing and result review.

[0017] The system quantifies uncertainty, provides probabilistic accounting results, and improves decision reliability: This invention uses an ensemble method to systematically quantify the combined impact of input errors, parameter uncertainties, and structural errors, and ensures their consistent propagation during assimilation and scenario extrapolation. The final output of key state quantities, fluxes, and potential results is accompanied by statistical information such as confidence intervals and risk probabilities, avoiding decision biases caused by providing only point estimates and better meeting the requirements of policy evaluation for reliability and risk representation.

[0018] This invention ensures the comparability of the counterfactual baseline and the action scenario, enhancing the attribution capability of emission reduction potential. It constructs both the counterfactual baseline and the action scenario from a unified posterior benchmark (i.e., a consistent set of assimilated states and parameters) and conducts scenario differential accounting under consistent initial states, parameter sets, and external driving settings. This effectively separates the confounding effects of the action effect from differences in climate, soil, and management backgrounds, giving the calculated emission reduction potential a clear causal explanation and attribution capability, resulting in more explanatory and credible results. Attached Figure Description

[0019] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings: Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention. Detailed Implementation

[0020] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0021] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0022] like Figure 1 As shown, this embodiment provides a method for agricultural nitrogen and phosphorus total flux assimilation and scenario-differential potential calculation, including: Based on standardized interface data packages, an assimilable agricultural nitrogen and phosphorus state-space model and forward prediction interface are constructed. Based on standardized interface data packages and forward prediction interfaces, an observation operator library and a unified error statistics framework are established. Based on the observation operator library and unified error statistics framework, information content diagnosis is performed and assimilation configuration is generated; Based on the forward prediction interface, observation operator library, unified error statistics framework and assimilation configuration, a constraint set of matter conservation constraints and physical feasible region constraints is constructed and implemented. Based on the forward prediction interface, observation operator library, unified error statistics framework, assimilation configuration and constraint set, perform constraint joint data assimilation solution to obtain unified posterior baseline data package; Based on the unified posterior baseline data package, a counterfactual baseline and action scenarios are constructed, scenario differential potential is calculated, and differential potential results containing uncertainties are output.

[0023] Furthermore, the process of obtaining standardized interface data packets includes: Determine the system boundaries, accounting units, and time base of the accounting objects; Collect and organize the driving input library, observation library, and prior library, and perform element coding and unit unification; Spatiotemporal scale matching and quality control are performed on the data in the driving input library and the observation library; The data that has undergone spatiotemporal scale matching and quality control is assembled into a standardized interface data package containing observations, driving inputs, quality labels, and metadata.

[0024] Furthermore, regarding data acquisition and standardized preprocessing in this embodiment, under a unified spatiotemporal benchmark and accounting caliber, an observation library, a driving input library, and a priori library for multi-source data are constructed, forming a standard interface data package that can be directly called by the assimilation solver. The standard interface data package is defined as follows: These represent: observation, driving force (management input), static covariate, quality identifier, and metadata, respectively; where i is the accounting unit and t is the time step or window index, thus providing a consistent data foundation, traceable data links, and reproducible data support for subsequent conservation constraint assimilation updates and scenario difference accounting. Specifically, this includes the following steps: S1.1 Define the system boundaries and accounting units; The system boundaries, accounting units, and time bases of the accounting objects are determined to unify subsequent data alignment, model operation, and flux closure. Specifically, spatial boundaries and accounting unit types (e.g., fields, grids, watersheds) are defined, and each accounting unit is assigned a unique identifier; vertical boundaries (soil profile depth and stratification) are determined, and receiving ports (surface water, groundwater, and atmospheric exchange boundaries) are clarified to support subsequent cross-port flux accounting; a list of accounting elements (key nitrogen and phosphorus reservoirs and major input, transformation, transport, loss, and output fluxes) is formed to establish a variable dictionary; at the same time, the time step or statistical window (in days, ten-day periods, months, or years) is determined, and alignment and aggregation rules (e.g., cumulative, average, flow-weighted average, etc.) and the recording method of the window start and end times are pre-defined when the observation time resolution and model time step are inconsistent, to ensure that the observation and model outputs are comparable and consistent.

[0025] S1.2 Multi-source data acquisition and element coding; Multi-source data were collected and organized, and divided into three categories according to their purpose: driving input database, observation database, and prior database. At the same time, the element coding, unit dimension, and morphological caliber were unified.

[0026] Driver input library This includes, but is not limited to, the amount and method of input and application of chemical or organic fertilizers, the amount and method of straw returning to the field, the amount of irrigation and the nitrogen and phosphorus concentration of irrigation water, the input of livestock and poultry manure, dry and wet deposition, as well as crop type and key management events, and can be expanded to include external driving variables such as meteorology and hydrology. Observation library At least including: soil ammonium nitrogen ( ), nitrate nitrogen ( Total nitrogen (TN), total phosphorus (TP) (or available phosphorus Olsen-P); water environment flow rate (Q) and concentrations of total nitrogen (TN), total phosphorus (TP), and dissolved active phosphorus (DRP); nitrous oxide in the gas ( ),ammonia( The flux of soil organic carbon (SOC) and soil organic nitrogen (SON) is considered. Other variables (such as soil organic carbon (SOC) and soil organic nitrogen (SON)) are used as optional weak constraints or covariates. The prior library includes basic soil properties and static covariates, initial library priors, empirical intervals and probability distributions of model parameters, and prior characterization information of input errors and process uncertainties, which are used to initialize and constrain the assimilation update range. At the same time, an element code and unit or morphology dictionary are established to uniformly label nitrogen and phosphorus forms, unit conversions, statistical calibers and data versions to ensure semantic consistency, traceability and reusability of cross-source data.

[0027] S1.3 Spatiotemporal consistency and scale matching; This step is used to uniformly map observation and driving data from different sources and supporting scales to the target accounting unit, time step, and soil stratification structure to ensure scale consistency and comparability in subsequent assimilation updates. At the spatial scale, a unified coordinate reference system and spatial representation are used, and point (or site) observations are mapped and merged into the target accounting unit. Mapping methods include, but are not limited to, nearest neighbor, area weighting, and spatial interpolation. Supporting scale information such as representative range, coverage, and buffer radius are recorded. At the temporal scale, observations are aligned and aggregated according to the model time step or window rules. Flux-related observations are preferably accumulated or integrated by window, while concentration-related observations are preferably averaged by window or flow-weighted average. When estimating load or flux from concentration, flow-weighted or concentration-flow integration methods are preferred, and the source of the hydrological data and window rules are recorded. At the vertical scale, soil observations are matched to the model stratification structure according to sampling depth. Interpolation or stratification weighting methods are used to convert the observations into equivalent state quantities corresponding to the model layers, and layer thickness weights, conversion coefficients, and methodological basis are recorded, thus forming a data representation that can be directly compared with the model output.

[0028] S1.4 Quality Control and Handling of Missing Test Items; This step is used to implement hierarchical quality control on the data entering the assimilation system and create traceable records to reduce the interference of erroneous data on assimilation and inversion and to provide a basis for subsequent error modeling. The quality control process includes at least four stages: rule verification, anomaly identification, anomaly handling, and missing data processing. Rule verification is used to unify unit dimensions, verify field meanings and sampling and inversion calibers, and set physical reasonable ranges and engineering thresholds, while also performing logical consistency checks (e.g., matching fertilization events with crop type or growth period, and the rationality of the time sequence of management events). Anomaly identification preferably uses a combination of robust statistical tests and physical tests. Robust statistical methods include, but are not limited to, IQR, robust z-score, or grouping. The process is hierarchical, categorized by "variable—season (or month)—spatial unit (or station)" to avoid misjudging seasonal peaks as anomalies. Physical inspection focuses on identifying negative values, out-of-range values, and orders-of-magnitude errors. For anomaly handling, data with clearly identifiable errors are removed or corrected and replaced. Observations that may represent genuine extreme events are retained but subject to reduced weighting or event window aggregation. The handling type, threshold basis, and result label are recorded. For missing data handling, continuous time-series variables are interpolated using time interpolation (e.g., linear interpolation, spline interpolation, or smoothing algorithms). Spatial variables are handled using borrowing from similar units or spatial interpolation. When key variables are difficult to interpolate effectively, prior model predictions are used for placeholder filling, but this must be clearly labeled as "prior filling" and given a large uncertainty indicator. The entire process must record the processing method, key parameters, and version information to support subsequent audit verification and uncertainty quantification.

[0029] S1.5 standard interface data packet assembly and output; This step is used to structure the multi-source data after it has been standardized and quality controlled into a standard interface data package that can be directly called by the assimilation solver, so as to ensure that the data is traceable, weightable, reproducible and supports subsequent scenario simulation. Specifically, a unified indexing system is established for each observation and driving data, including at least spatial indexes (e.g., accounting unit ID), temporal indexes (e.g., time step t or window start and end times), and vertical indexes (e.g., layer number or sampling depth). Data source, acquisition method, and product version information are also recorded. The observation data is standardized and organized, including at least variable type, observation value, unit, timestamp or window, spatial unit, layer depth information, sample size or coverage, and data aggregation rules (e.g., cumulative, average, flow-weighted average, etc.) to ensure comparability with model output in terms of time scale and spatial support. Quality and error attribute labels are added to each observation data, including error source type (e.g., instrument error, inversion error, representativeness error, interpolation error, etc.) and quality level (e.g., raw pass, converted, reduced-weight retention, interpolation, removal and archiving, etc.). Prior descriptions of observation uncertainty (e.g., standard deviation, interval or level mapping parameters) can be selectively recorded for subsequent construction of the error statistics framework R. Finally, observation vectors are assembled according to time step or window. (Grouping can be done by variable type or data source to support batch assimilation), and synchronous output drives the input sequence. Static covariates Quality labeling With metadata This forms a standard interface data package for assimilation solvers.

[0030] Furthermore, the process of constructing an assimilable state-space model for agricultural nitrogen and phosphorus and an interface for forward prediction of total flux includes: Construct an extended state vector to describe key nitrogen and phosphorus pools, crop nutrient status, and optional deviation states; Define a full flux vector that explicitly decomposes the input, transformation, transfer, loss, and output processes into traceable components; Define the driving inputs as the entry point for management measures, the static covariates characterizing spatial heterogeneity, and the set of parameters to be estimated; Based on the extended state vector, full flux vector, driving input, static covariates, and parameter set, a forward prediction interface is established to advance the state and synchronously output the full flux at time steps.

[0031] Furthermore, this embodiment constructs an assimilable agricultural nitrogen and phosphorus state-space model and a total flux forward prediction interface. This is used to represent the agricultural nitrogen and phosphorus cycle process as an assimilable state-space model under a defined accounting unit, vertical stratification, and time step, and to clarify the state variables that need to be jointly estimated. Total flux vector With parameter set And establish a mechanism that can output the state of the next time step at each time step. and the total flux at the current step Forward prediction interface The model is organized around the core principle of "storage-flux" consistency: on the one hand, changes in storage are determined by traceable inputs, transformations, transports, losses, and output fluxes; on the other hand, each flux has a clear port affiliation (e.g., surface water, groundwater, atmosphere, harvest removal, etc.) and statistical caliber (step rate or window accumulation), thus providing a unified model foundation for subsequent conservation closure constraints, observation operator mapping, and scenario difference accounting.

[0032] Specifically, the following steps are included: S2.1 Constructing the Extended State Vector ; This sub-step is used to define the extended state vector. This enables it to carry the key nitrogen and phosphorus stockpile and the necessary process control states in the assimilation and update process, and to correspond one-to-one with the accounting unit or hierarchical structure. It must contain at least one of the key states: nitrogen pool or phosphorus pool; preferably, the nitrogen pool includes the soil inorganic nitrogen pool. , The phosphorus library (expressed by layer or by cross-section) and optional organic nitrogen libraries (e.g., active, inert organic nitrogen, mineralizable nitrogen, etc.); the phosphorus library preferably includes available phosphorus, dissolved phosphorus, adsorbed or precipitated phosphorus, and optionally includes organic phosphorus.

[0033] To ensure that agricultural output can be mapped to remote sensing or statistical observation data, Preferably, crop nutrient status (e.g., crop nitrogen and phosphorus accumulation, crop nitrogen and phosphorus uptake pool, and biological nitrogen and phosphorus content) is included. To improve the interpretability and observability of the process calculations, A limited number of auxiliary parameters can be added (such as soil moisture content, water availability index, temperature index, permeable water volume, runoff indicator, etc.).

[0034] Furthermore, to address the systematic biases of multi-source observations and support a subsequent unified error framework, Optional extensions include biased or stochastic effects states, which introduce slowly varying bias terms over time for specific observation source categories (e.g., remote sensing inversion, statistical inputs, cross-sectional water quality and gas inventories), allowing them to be estimated along with the main state during assimilation, thereby avoiding unreasonable mutual interference between observations from different sources during updates.

[0035] S2.2 Defines the total flux vector ; The input, transformation, transport, loss and output processes of nitrogen and phosphorus in the current step are explicitly decomposed into a set of traceable flux components and form a closed interface with state updates.

[0036] It should include at least external input fluxes (nitrogen and phosphorus from chemical fertilizers, organic fertilizers or manure, irrigation, dry and wet deposition, nitrogen fixation, or external recharge), agricultural output fluxes (e.g., nitrogen and phosphorus removed during harvest and straw removal), and major loss fluxes. Nitrogen loss fluxes should include at least leaching, seepage, and runoff into the aquatic environment, as well as volatilization and greenhouse gas emissions into the atmosphere (e.g., nitrogen loss). , (etc.); phosphorus loss flux includes at least dissolved and particulate nitrogen and phosphorus carried out by runoff and erosion.

[0037] To ensure the feasibility of subsequent conservation closure, each flux component should be clearly defined: (1) The port of origin (e.g., surface water, groundwater, atmosphere, or harvest output, etc.); (2) Statistical caliber (e.g., instantaneous rate, window cumulative rate, etc.); (3) Units and conversion rules; (4) Expressions that can be mapped to observations (e.g., concentration-flow integral to obtain load or flux).

[0038] For internal transformation processes (such as mineralization, nitrification and denitrification, adsorption-desorption, precipitation-dissolution, etc.), you can choose to "explicitly list and output" or "output as a net effect term" depending on the model complexity. However, you must ensure that it does not disrupt the closure relationship between the stock volume and flux and can be connected with subsequent constraints and diagnostic indicators (such as closure residuals).

[0039] S2.3 Setting Drive Input Static covariate s and parameter set Distinguish between "measure entry points" and "estimateable parameters"; This sub-step is used to inject management measures and external drivers into the model in a unified manner, and to clarify which parameters need to be assimilated and estimated, and which should be kept constant, thus providing a clear and actionable entry point for subsequent scenario construction. Driver Input It should include at least the sequence of management events and external input information, such as fertilization events (type, amount, frequency, timing, application method and depth), organic fertilizer or manure input, irrigation volume and nitrogen and phosphorus concentration in irrigation water, straw return amount and method, tillage and cover crop management, etc., and can be extended to meteorological, hydrological and other driving sequences; static covariates s should include at least basic soil properties and spatial background information (such as texture, bulk density, pH, SOC, topography, land use type, field, watershed characteristics, etc.), used to characterize spatial heterogeneity and participate in process calculation and parameter constraints. Parameter set This is used to characterize key process parameters (such as reaction rate coefficient, partition coefficient, migration coefficient, adsorption isotherm-related parameters, crop uptake efficiency parameters, etc.), some of which can be set as constants, while others can be set to change slowly over time. To express seasonal differences or management changes; for each parameter, at least the prior range and probability distribution type should be given so that subsequent joint updates can be performed under a unified error statistics framework and constrained by the physical feasible domain.

[0040] S2.4 Establish the forward prediction interface Output at each step and And retain closed-loop verification quantities; This sub-step is used to advance the system state to the desired state at time step t. Forward prediction interface And calculate the total flux at each step simultaneously during the advancement process. This creates an assimilable computing engine that integrates state and flux. The interface can be described as follows: ; in, This represents process noise or structural error terms, used to characterize the impact of undescribed processes, driving errors, and model structural uncertainties on state progression, and serves as the entry point for subsequent process error modeling and uncertainty propagation. Indicates by The next predicted state vector obtained from the propagation, Represents the total flux vector. This represents the system state vector at time step t. The external driving and management input vector represents time step t, and s represents the static covariate. Represents a set of parameters. This indicates process noise or structural error.

[0041] To ensure the feasibility of the project and the traceability of flux, Modular organization is preferred: The input allocation module is used to allocate inputs such as fertilization, irrigation, and sedimentation to the soil layer and morphology library; The crop uptake module is used to calculate crop nitrogen and phosphorus uptake and harvest removal and to maintain mapping compatibility with yield and biomass-related observations; The soil transformation module is used to calculate the net effects of processes such as mineralization, nitrification and denitrification, phosphorus adsorption-desorption and precipitation-dissolution; The transport and loss module is used to calculate loss fluxes such as leaching, seepage, runoff erosion, and gas emissions under given hydrological conditions and to determine port attribution. The Status Update and Throughput Summary module is used to summarize all throughput components according to a unified standard and update the status of each inventory.

[0042] To enhance the executability of subsequent constraints and diagnostics, the interface can optionally output closure check quantities (such as the nitrogen and phosphorus mass balance residuals or branch-end income and expenditure items in the current step), which can be directly called in subsequent conservation constraint execution and consistency checks; at the same time, It may provide only the forward simulation interface, or provide linearization information (such as incremental operators or Jacobi information) when needed to adapt to different assimilation solver implementations, but in any form, it should meet the closed interface requirements of "traceable full flux components, closed-loop inventory updates, and accountable port affiliation".

[0043] Furthermore, the process of establishing the observation operator library and the unified error statistics framework includes: Define corresponding observation mapping rules for the classification of multi-source observation data, and construct an observation operator library; For each type of observation data, an observation error model is established, and the observation error covariance matrix is ​​constructed. Establish a process error or background error model to describe the uncertainty of model structure, driving forces, and parameters; By introducing a bias state or reducing the weight of representative errors, the inconsistency between system bias and scale support is explicitly addressed.

[0044] Furthermore, this embodiment establishes an observation operator library H(·) and a unified error statistics framework, which are used to process the output standard observation package. With the constructed state-space model (state) Flux ,drive ,parameter With forward prediction interface Align and connect them within the same error statistics framework.

[0045] Specifically, this step involves establishing a reusable library of observation operators for observations from different sources, with different support scales, and different statistical calibers. This ensures that each type of observation can be mapped to the model state or flux space with explicit rules and is comparable to the model output. On the other hand, a unified error statistics framework is constructed, including the observation error covariance R and its structure, process error, background error covariance Q and its prior structure. Systematic bias and representative error are incorporated into the assimilation system in a "modelable, weightable, and traceable" manner, thereby providing the necessary key inputs for subsequent constraint joint assimilation solutions. Consistent weighting criteria.

[0046] S3.1 Classification of Observation Types and Definition of Observation Mapping Rules: Family of Observation Operators This is used to classify multi-source observations according to their "physical meaning and statistical caliber" and to define a corresponding observation operator for each type of observation. This ensures that observations and model outputs are strictly comparable in terms of spatial support, time window, and vertical level.

[0047] The observation operator must cover at least the following general mapping patterns: One type is the state-observable operator, used to analyze soil conditions. , Monitoring values ​​such as Olsen-P are mapped to the corresponding library quantities or their equivalent state representations in the model state vector; The second is a vertical layer weighting or interpolation operator, which is used to convert soil observations at different sampling depths into equivalent state quantities corresponding to model layers; The third is the time window aggregation operator, which is used to convert the model output from the time step scale to the cumulative or average quantity corresponding to the observation window; The fourth type is the concentration-flow or confluence coupling operator, which is used to combine the area source output generated by the model with the hydrological process to calculate the cross-sectional concentration, load or flux; The fifth is the remote sensing inversion or proxy variable operator, which is used to map crop biomass, nitrogen and phosphorus absorption status or related ecological process quantities in the model into remote sensing inversion quantities (such as leaf area index LAI, photosynthetically active radiation absorption ratio FAPAR, yield proxy variables, etc.), and to clarify the version of the inversion product and the statistical scope. The sixth type is the weakly constrained operator for inventory or inversion, which is used to weakly constrain greenhouse gas inventory or inversion results to the corresponding gas loss flux components in the model.

[0048] Through the above classification and definition, each observation entry can be associated with a clear operator type identifier and scale processing rules, thereby forming a reusable, scalable, and auditable observation operator library.

[0049] S3.2 Observation Operator Implementation and Interface Specification; This sub-step is used to engineer the family of observation operators defined in S3.1 into an interface that can be directly called by the assimilation solver, and to ensure that it is consistent with the observation package index system (spatial cell ID, time index or window, vertical layer number).

[0050] The interface is preferably satisfied as follows: given the model state or flux at any time step t. and necessary driving and static covariates The simulated values ​​corresponding to the observations can be calculated for all of them. It returns alignment results corresponding one-to-one with the observation entries (including statistical calibers for window cumulative or average, layer weighting coefficients, and key intermediate quantities such as confluence or flow weights for traceability). To support "batch assimilation" or "group update," the operator interface can be called in groups according to the observation source category or variable category, and the output of each group is mapped to the comparable simulated observation vector and its index. For implementation routes that adopt variational or optimization-based assimilation or require sensitivity analysis, the operator interface can optionally provide linearization information (such as local incremental mapping or Jacobian information). However, if a set-based assimilation method is used, only the forward operator output can be provided without requiring explicit linearization, thereby reducing the engineering implementation threshold while ensuring versatility.

[0051] S3.3 Construction of the observation error model and the observation error covariance R; This sub-step is used to establish an observation error model for each type of observation within a unified error statistics framework, and to construct the observation error covariance matrix R or its equivalent representation accordingly, so as to achieve comparable weighting of observations across data sources, scales, and calibers. The observation error model considers at least the sources of measurement error, sampling error, inversion error, preprocessing-introduced errors (e.g., interpolation, aggregation, layer weighting conversion errors), and representativeness errors, and allows the error to exhibit heteroscedasticity characteristics with the size of the observation, season, spatial unit, or quality level. Specifically, the error of remote sensing inversion quantities can be derived from product quality labels and inversion uncertainty estimation; the water quality concentration-load conversion error can be characterized by a combination of flow uncertainty and window integration error; and the error of gas flux observations or inventories can be given a larger variance according to the weak constraint level to reflect its representative uncertainty. The construction of R and the output quality label are optimized. Direct correlation involves applying different error amplification rules to different levels of "original pass, conversion, weight reduction retention, interpolation, and prior imputation," thereby achieving weight adjustment without altering the observations themselves. Simultaneously, R can optionally support temporal or spatial correlation structures (e.g., temporal correlation of continuous observations at the same station, spatial correlation of adjacent sections within the same watershed), recorded with structure labels or sparse representations for efficient solver retrieval. Through these rules, a correlation with the observation vector can be obtained. A consistent, traceable, and scalable statistical description of observation errors provides a unified weighting basis for subsequent joint assimilation updates.

[0052] S3.4 Process error or background error model and Q setting; Specifically, statistical representations of process error or background error are constructed to describe the impact of model structural uncertainty, driving uncertainty, and parameter uncertainty on state progression, and to provide a prior covariance structure for subsequent assimilation updates. For filtering or ensemble assimilation frameworks, process error can be represented as... Q is used to control the variability of the state and the magnitude of error spread during time progression, and different variance scales and correlation structures can be set according to variable type and hierarchical structure. For implementation using a static background field, the uncertainty of the prior state and parameters can be represented as background covariance Q, and prior samples can be generated by set perturbation, parameter perturbation or random walk.

[0053] The error model preferably distinguishes between "state error" and "parameter error": state error reflects the impact of unmodeled processes and driving errors on inventory and throughput, while parameter error allows key process parameters to be updated within a reasonable range or drift slowly with the seasons; simultaneously, error inflation or variance lower bound mechanisms can be introduced to avoid set collapse or overconfidence. Through these settings, a model can be formed that... The error propagation entry point of the docking allows the assimilation update to both utilize observations for correction and maintain the consistency of system dynamics within a reasonable range of uncertainty.

[0054] S3.5 Explicit handling of systematic bias and representativeness error; This sub-step is used to explicitly address the systematic bias and scale inconsistency issues in multi-source observations within a unified framework, avoiding conflicts between different data sources during assimilation updates and preventing flux allocation distortion. For observation sources with significant systematic bias risks (e.g., remote sensing inversion, statistical inputs, differences in cross-sectional water quality caliber, systemic differences in gas inventories, etc.), a "biased state incorporated into the state vector" approach is preferred for modeling. This involves the bias term being slowly changed over time as part of the extended state and updated along with the main state, thus separating the "bias" from the "real process changes." For representative errors caused by scale mismatch, it is preferable to include them as a component of R or as random effects in the observation model, and associate them with representative metadata (coverage, buffer radius, supporting scale level, etc.) to form interpretable error amplification or weighting rules. If necessary, a "soft constraint" mechanism (increasing variance, grouping updates, or event window aggregation) can be applied to specific observations to ensure that the constraint strength matches their representativeness, thereby improving the stability of cross-source fusion and the physical consistency of posterior results.

[0055] Furthermore, the process of performing information content diagnosis and generating assimilation configurations includes: The contribution of different observation sources to the constraints of state, flux and parameters is evaluated, and observability analysis is performed to obtain the observability analysis results. Based on the observability analysis results, assimilation variables were screened and stratified. Develop assimilation and update strategies in groups or phases; Regularization and flux decomposition rules are introduced for weakly identifiable flux components; Coordinate the configuration of the weights of soft and hard constraints, and interface with the constraint set and error statistics framework; Perform a consistency check before assimilation and output an assimilation configuration file containing a list of variables, update order, regularization rules, and weight configuration.

[0056] Furthermore, this embodiment also involves information content and identifiability diagnosis and assimilation configuration, used to systematically evaluate the effect of existing observation combinations on state variables before performing joint assimilation with conservation constraints. flux vector With parameter set The constraint capability identifies key variables that are "observable and identifiable" and variables that are "weakly constrained and unidentifiable," and generates feasible assimilation configuration schemes accordingly. The core purpose of this step is to avoid non-unique solutions, flux allocation drift, or overfitting caused by "high dimensionality of total flux variables but insufficient observational information." Through strategies such as group updates, regularization, and synergy between soft and hard constraints, subsequent assimilation solutions obtain stable, interpretable, and engineering-usable posterior results while satisfying conservation closure and the physical feasible region.

[0057] Specifically, it includes: S4.1 Observability and sensitivity assessment to determine the contribution of observations to state-flux-parameter constraints; This sub-step is used to quantitatively evaluate different observation sources and their observation operators. The sensitivity and information contribution of each state variable, flux component, and parameter are assessed to determine which variables can be effectively constrained. Evaluation methods can include local perturbation sensitivity, ensemble perturbation response, incremental mapping analysis, or empirical information content indicators based on innovation statistics. That is, given a prior sample or typical state, the sensitivity and information contribution of each state variable are assessed. , and Small perturbations are applied to each component and the observed simulated values ​​are calculated. The response magnitude and direction are analyzed to characterize the "sensitivity or identifiability of observations to variables." Simultaneously, a stratified assessment is conducted according to "variable—season (or growing season)—spatial unit" to reflect the impact of management events and seasonal processes on observability. The assessment results are used to form an observation contribution matrix or contribution score, clarifying the strength of constraints imposed by different observation sources (remote sensing, soil, water quality, groundwater, gas, etc.) on different process chains (crop uptake, soil transformation, hydrological loss, gas loss, etc.), providing a basis for subsequent variable selection and grouping updates.

[0058] S4.2 Assimilation variable selection and hierarchical expression: determining the main assimilation variable, auxiliary variable and derived output; This sub-step is used to screen and stratify the set of states, fluxes, and parameters to be assimilated based on the observability assessment results of S4.1, forming a structured list of "assimilation main variables - weakly constrained variables - derived output variables". For variables with sufficient observational support and critical to the accounting objective (e.g., crop uptake and harvest removal related states, soil inorganic nitrogen, key pools such as Olsen-P, total water environmental output load, etc.), they are preferably included in the set of assimilation main variables and participate in direct updates. For variables with weak observational constraints but necessary components of the flux chain (e.g., subdivided gas loss components, some internal transformation components, etc.), weakly constrained updates can be used, or they can be derived and allocated only under the condition of conservation closure by the main variables and total constraints. For variables that are difficult to identify and have limited contribution to the objective, they are preferably set as derived outputs or fixed as prior empirical values ​​to reduce dimensionality and improve posterior stability. Through the above stratification, without sacrificing flux closure and accounting integrity, too many weakly constrained variables can be avoided from directly entering the assimilation, which could lead to instability and decreased interpretability of the solution.

[0059] S4.3 Grouping or phased assimilation strategy and update order; This sub-step is used to formulate grouped or phased assimilation and update strategies to leverage the differences in information structure and process time scales of different observation sources, achieving a robust update path of "stabilizing key reservoirs and total amounts first, then refining allocation." Preferred strategies include: first, using observations with the strongest system constraints and relatively stable caliber (e.g., yield, biomass, fertilization statistics, key soil reservoir amounts) to update crop uptake and the main soil reservoir, ensuring the system is within the physically feasible domain and under reasonable water and fertilizer conditions; second, introducing water environment observation indicators (e.g., cross-sectional concentrations and loads, groundwater nitrates, etc.) to constrain hydrological transport and total water environment output; and finally, introducing gas flux observations or inventory-type weak constraints to refine gas loss allocation or provide boundary checks. For processes with strong seasonality and event-response characteristics (e.g., short-term gas peaks after fertilization, peak runoff loads from heavy rainfall), aggregation can be performed by event window and processed separately using an "event update" method to avoid high-frequency spikes unreasonably dominating long-term state updates. The aforementioned phased strategy should clearly define the observation group, update frequency, window rules, and quality threshold for each phase, and output the update order configuration that can be directly called by the solver.

[0060] S4.4 regularization and flux decomposability rules solve the "closed but not unique" flux allocation problem; This sub-step introduces implementable regularization and decomposition rules for weakly identifiable flux components to obtain stable and interpretable flux allocation results while satisfying conservation closure and observation fitting. Regularization methods may include, but are not limited to: applying hierarchical prior constraints with shared hyperparameters to similar flux components; using contraction priors or penalty terms to suppress unfounded fluctuations in expected sparse or low-amplitude fluxes; and introducing smoothing or random walk constraints to parameters to limit unreasonable jumps. The flux decomposition rules preferably adopt a "total first, then allocation" structure, i.e., first estimating a certain type of total flux (e.g., total nitrogen loss or total water environment output) using strongly constrained observations, and then combining auxiliary observations such as water quality and gas observations to decompose the total flux into sub-fluxes such as "hydrological loss or gas loss" according to interpretable allocation coefficients, and applying physical feasible region constraints (0–1, non-negative, sum to 1, etc.) to the allocation coefficients. The introduction of regularization and decomposition rules effectively avoids multiple solutions or flux allocation drift when observational information is insufficient, thereby improving the robustness and engineering usability of the posterior results.

[0061] S4.5 Soft and Hard Constraint Coordination and Weight Configuration; This sub-step connects the error statistics framework with the constraint set at the solver configuration level, clarifying which constraints are executed as hard constraints and which are incorporated into the objective function or update rules as soft constraints, and providing a traceable weight configuration scheme. For constraints that are clearly expressible and crucial to physical consistency, such as mass conservation closure, non-negativity, and key boundary conditions, hard constraints are preferred, or they are enforced through projection and variable transformation. For constraints significantly affected by observation caliber and representativeness errors (e.g., the correspondence between certain cross-sectional loads and area source outputs, inventory-type gas emissions, etc.), soft constraints can be set, and "weak constraints" can be achieved by increasing uncertainty or reducing weights. Simultaneously, the weights for different observation sources can be stratified based on quality level, representative metadata, and innovation statistics, and the threshold basis and version information are recorded in the configuration file. This configuration output will serve as the parameter input for the subsequent assimilation solver, ensuring a controllable balance between "observation fit—prior rationality—constraint feasibility" in the update process.

[0062] S4.6 Pre-assimilation consistency check and entry threshold; This sub-step performs a pre-assimilation consistency check and admission determination before entering the joint assimilation solution, to avoid solution failure or posterior anomalies caused by obviously inconsistent data combinations or configuration errors. The check includes at least: whether the observation and model scales are aligned; whether the observation operators can generate simulated values ​​for all observations; whether the error variance meets the non-negativity and reasonable range; whether the prior range of key variables is consistent with the physical feasible region; and whether the grouping update order covers the key variables required for the computational objective. For cases that fail the admission threshold, data correction or error parameter adjustment can be performed, or the assimilation dimension can be reduced, regularization enhanced, or the soft and hard constraint configuration adjusted. This pre-assimilation threshold mechanism significantly improves the stability and reproducibility of subsequent solutions and provides a clear quality control point for engineering applications.

[0063] S4.7 outputs the assimilation configuration file for the solver to use directly; This sub-step is used to summarize and output the assimilation configuration file, which includes at least: a list of assimilation main variables, weak constraint variables, and derived output variables; observation grouping and update order, update frequency, and window rules; regularization and flux decomposition rules and their parameters; soft and hard constraint co-setting and weight configuration; and pre-assimilation entry thresholds and checking rules. The configuration file, along with the observation package, model interface, and error statistics description, constitutes a reproducible input for subsequent conserved constraint joint assimilation solutions, enabling the overall scheme to be quickly migrated and stably run under different regions and data conditions by adjusting the configuration.

[0064] Furthermore, the process of constructing and enforcing the constraint set includes: Construct equation constraints for the conservation of nitrogen and phosphorus and the total flux closure; Construct inequality constraints to ensure that state variables, fluxes, and key derived variables are non-negative and within the physical boundary; Construct a parameter feasible region and scale constraints to ensure that process parameters are physically feasible and spatially reasonable; Construct process coupling constraints to ensure the consistency of substrate constraints, carrier constraints, and morphological transformation in order to guarantee the mechanistic self-consistency of the process; Define a constraint implementation operator for performing feasible region correction on candidate solutions and outputting a closure verification diagnostic quantity.

[0065] Furthermore, this embodiment constructs and operatorizes a set C of matter conservation constraints and physically feasible region constraints; specifically, after completing the observation operator library... After constructing the error statistics framework, to ensure the assimilation and update of the state variable x, total flux f, and parameters... To simultaneously satisfy the principles of matter conservation, physical feasibility, and consistency across processes and scales, this embodiment constructs a set of constraints. This involves transforming the operators into callable constraint enforcement interfaces. The constraint set is preferably represented as follows: ,in Characterizing the conservation of nitrogen and phosphorus and the closed-loop equation constraint of total flux. It characterizes nonnegativity and boundary constraints, parameter feasible region constraints, and necessary process coupling or scale consistency constraints. The constraint enforcement operator is used to perform feasible region correction and closure verification on the candidate solutions obtained for each assimilation window (or after each update), and outputs diagnostic quantities such as closure residuals and default degree. In addition to observational statistical constraints (characterized by R) and model uncertainties (characterized by Q or set perturbations), it introduces an auditable conservation and feasible region constraint mechanism to improve the reliability, interpretability, and verifiability of the verification results.

[0066] S5.1 defines the mass conservation and flux closure equation constraints for nitrogen and phosphorus: ; For each accounting unit (e.g., field, grid, watershed, and administrative region) and each assimilation window We construct closed-loop material conservation relationships for nitrogen and phosphorus respectively to constrain the updated state quantities, flux quantities, and parameter quantities to satisfy the closed-loop rule of "inventory change = total input within the window - total output within the window".

[0067] Taking nitrogen as an example, define the total nitrogen inventory of the system within the window. The sum of the stockpiles of each nitrogen form in the state vector and the measurable stockpiles such as crop nitrogen stock; the total input within the window is defined. With total output For the accumulation of input or output components within a window in the total flux set (using discrete summation or continuous integration), the following condition is satisfied: ; Similarly, define the total phosphorus inventory in the system. Total Input With total output ,satisfy: ; in, This is a closed-loop residual term used to accommodate discretization errors, small throughput not explicitly characterized, and statistical caliber differences, and serves as the output for closure verification and quality diagnosis; preferably... Set a threshold or variance to ensure closure while avoiding infeasible solutions due to representativeness errors.

[0068] and The flux component list and port attribution should be consistent with the total flux vector. One-to-one mapping and embedding in the constraint configuration. For nitrogen, It includes at least nitrogen from chemical fertilizers, available nitrogen from organic fertilizers or manure, nitrogen from straw returned to the field, nitrogen from dry and wet sedimentation, nitrogen from biological nitrogen fixation, and nitrogen input from irrigation. At least including nitrogen carried away during harvest, volatiles and gaseous emissions ( , Optional ), leaching and leakage (with (primarily) and runoff erosion output. Regarding phosphorus, It includes at least phosphorus from chemical fertilizers, available phosphorus from organic fertilizers or manure, phosphorus from straw returned to the field, and phosphorus input through irrigation. It includes at least harvest-carried phosphorus, runoff dissolved phosphorus (DRP), and erosion particulate phosphorus (PP) outputs, and optionally includes seepage-output phosphorus.

[0069] Preferably, when data permits, decomposition closure constraints can be added for sub-ports or hierarchical structures, but should remain consistent with the overall system closure to enhance identifiability and suppress non-physical flux compensation.

[0070] S5.2 Defines the physically feasible region and boundary inequality constraints: ; To ensure the state variables are updated after assimilation , full flux And given that the key derived quantities satisfy physical feasibility, this invention constructs an inequality-constrained subset. Nonnegativity and boundary constraints are applied to each accounting unit and each assimilation window (or time step). The nonnegativity constraints include at least: ; Wherein, the state vector Preferably, it should contain at least a stock of nitrogen and phosphorus forms from each layer; for example etc.; among them, This represents the stock of ammonium nitrogen, nitrate nitrogen, and organic nitrogen in the z-th soil layer. This represents the active or available phosphorus stock in the z-th soil layer. This represents the inorganic phosphorus buffer stock on the surface of soil minerals in the z-th layer (such as iron / aluminum oxides, clay minerals, etc.) that can be reversibly adsorbed and desorbed. This indicates the organic phosphorus stock, etc., in the z-th soil layer that participates in mineralization or fixation processes. Preferably, it should include at least the main loss or output flux component, for example This avoids non-physical outcomes such as negative inventory, negative losses, or negative output. This represents the nitrogen loss flux that escapes into the atmosphere due to ammonia volatilization. This represents the nitrogen loss flux emitted into the atmosphere via N2O. This represents the nitrogen output flux that enters deeper layers or groundwater receptors through leaching or deep percolation. This represents the nitrogen output flux entering surface water receptors via surface runoff. This represents the dissolved reactive phosphorus flux exported with surface runoff. This represents the particulate phosphorus flux transported by erosion.

[0071] Unless negative, boundary constraints should preferably cover key moisture states and hydrological transport conditions to ensure that process drivers are within a reasonable range, such as setting limits on soil moisture content or equivalent moisture states. ; It also constrains the transport of water quantities such as runoff, seepage, and erosion to be non-negative to avoid abnormal flux amplification or sign errors caused by water state exceeding limits. This refers to soil moisture content or equivalent water state. Saturated soil moisture content, z is the soil depth.

[0072] For key process parameter vectors For parameters such as rate coefficient, migration coefficient, adsorption parameters, and efficiency parameters, it is preferable to apply prior feasible domain boundary (upper and lower bound) constraints; for flux distribution coefficients and efficiency coefficients (e.g., volatiles ratio, etc.), (Allocation ratio, crop absorption efficiency, etc.) Preferably, interval constraints are applied: ; If necessary, an allocation constraint of "summing to 1" can be added to ensure the interpretability and stability of the total allocation.

[0073] Optionally, to avoid the generation of obvious non-physiological crop states by assimilation updates, reasonable upper and lower limits or growth boundaries can be set for crop physiological derivatives (such as LAI, biomass, and yield proxy variables), and these can be correlated with the crop type's growth period or historical statistical range to suppress abnormal growth states caused by local observation errors.

[0074] To facilitate the implementation and verification of constraints, this step defines a uniform degree of breach for inequality constraints. : The location and magnitude of the breach are recorded; the diagnostic parameters are used by the subsequent constraint implementation operator to perform feasible region correction and quality classification on the candidate solutions.

[0075] S5.3 Feasible region and scale constraints for parameters; To ensure the physical feasibility and transferability of the updated process parameters and to avoid achieving surface fitting through unreasonable parameter drift, this invention optimizes the process parameter vector. Construct the feasible region and scale constraints, and incorporate them into the constraint set. Constraints can be used as hard constraints to directly limit the range of parameter values, or as soft constraints to participate in updates in the form of penalty terms or prior constraints. In addition to satisfying observational and statistical constraints and conservation closure constraints, these constraints further ensure the interpretability and engineering usability of parameters.

[0076] First, apply basic feasible region boundary constraints to the parameter vector, satisfying the following for each parameter component: ; in This is the minimum value of the parameter vector, i.e., the lower bound; The upper bound is the maximum value of the parameter vector. The upper and lower bounds can be determined based on literature review, experimental calibration results, soil property estimation, or regional empirical statistics, and can be stratified according to soil type, crop system, or management zone. For parameters such as rate constants, migration coefficients, and distribution coefficients that span orders of magnitude and are strictly positive, it is preferable to impose constraints in the logarithmic space to improve numerical stability, i.e., let... and in Applying boundaries or regularization to the space achieves a unified expression of "positive constraints + orders of magnitude constraints".

[0077] Secondly, to ensure the consistency and spatial rationality of parameters across different accounting scales, this invention constructs hierarchical consistency and spatial continuity constraints. For parameters with spatial continuity or zone-shared attributes (such as runoff control parameters, adsorption buffer parameters, and leakage control parameters), a "zone consistency or neighborhood smoothing" constraint can be optionally introduced. This allows parameters to be shared within the same soil type zone or the same management zone, or allows parameters to change smoothly between adjacent units, thereby reducing pixel-level noise and enhancing cross-regional generalization ability. This constraint can be expressed in the form of "bounded adjacent differences" or "smoothing regularization," for example, for adjacent units... satisfy: ; Or use a hierarchical representation And for the disturbance term Apply small-scale constraints to achieve scale consistency of "regional dominance + local fine-tuning".

[0078] Furthermore, to prevent assimilation from replacing mechanistic explanations by infinitely amplifying management input deviation terms, this invention sets a truncated feasible region centered at 1 for management input deviation terms or correction factors. For management-related multiplicative factors such as fertilizer application deviation coefficient, manure effectiveness coefficient, and irrigation input deviation coefficient, the following is preferably satisfied: ; in, Uncertainty can be set according to data quality level and statistical scope (the range can be appropriately widened for lower quality) to avoid masking process structural defects or causing non-transferable posterior solutions through "unbounded amplification of deviation terms".

[0079] By introducing the aforementioned feasible region of parameters and scale constraints, in addition to the constraints of conservation closure and physical feasible region, the degrees of freedom and spatial form of parameter updates can be further restricted. This ensures that the assimilation results maintain mechanistic consistency and scenario extrapolation usability while fitting observations, and provides auditable constraint basis for subsequent output parameter fields and uncertainty assessments.

[0080] S5.4 Process Consistency and Coupling Constraints; To prevent assimilation updates from obtaining surface fit through unreasonable flux compensation, and to ensure the updated state variables... Flux With parameters Mechanistically consistent, this invention introduces process consistency and coupling constraints, consistent with the agricultural nitrogen and phosphorus cycle mechanism, as a constraint set, in addition to the matter conservation closure (S5.1), physical feasible region (S5.2), and parameter feasible region (S5.3). An important component. Constraint optimization is expressed in the form of inequalities. The expression (which can be embedded in assimilation as a hard constraint or as a soft constraint in the form of a penalty term when necessary) can be organized into three categories: "substrate or inventory constraints, carrier constraints and morphological transformation consistency" to enhance the identifiability of assimilation inversion and suppress non-physical compensation between different fluxes.

[0081] (1) Substrate or stockpiling constraints. Any output or loss flux within the assimilation window should not exceed the available stockpiles or migratable share to avoid the non-physical outcome of "small stockpiles but large losses." Preferably, substrate constraints are imposed on leaching and seepage outputs, such that nitrate nitrogen leaching loss fluxes within the window are jointly limited by deep nitrate nitrogen stockpiles and leaching carriers; active phosphorus stock constraints are imposed on runoff dissolved phosphorus outputs; and soil particulate phosphorus content and erosion amount constraints are imposed on erosion particulate phosphorus outputs. These constraints can be expressed in the form of "cumulative fluxes within the window do not exceed the available supply," for example, for any flux component. satisfy: ; in, This represents the upper limit of available supply or the share of available migration determined by inventory and parameters, and can be implemented in layers z of soil to improve the interpretability and stability of the constraints.

[0082] (2) Carrier Constraints. The transport flux coupled with water or sediment carriers should be constrained by carrier strength to avoid flux compensation through abnormally high concentrations when the carrier is small. Preferably, the nitrogen and phosphorus output carried by runoff or seepage satisfies the bounded relationship of "flux = carrier quantity × effective concentration (or unit output coefficient)", and the phosphorus output from erosion particulate matter satisfies the bounded relationship of "flux = erosion quantity × particulate matter content × enrichment coefficient". This type of constraint can be achieved by limiting the "unit carrier output strength", for example, by satisfying the requirement that seepage flux meets the following condition. ; For dissolved phosphorus in runoff, satisfy ; in This refers to the leakage and runoff volume (or its cumulative amount) within the window. This represents the upper limit of concentration or unit output within the a priori feasible region, and can be stratified according to soil type, season, or management zone. For erosion-induced particulate phosphorus output, the following can be selected: ; in, This refers to erosion volume or sediment flux. This refers to the particulate phosphorus content in the soil. This is the upper bound of the enrichment ratio or its equivalent parameter bound.

[0083] (3) Consistency of Form Transformation and Directional Constraint. The internal transformation process of nitrogen and phosphorus should be consistent with the direction of change in the corresponding stock volume and meet the necessary allocation and proportional relationships to avoid sign mismatch or "cyclic compensation". For nitrogen, nitrification flux consumption and generate Denitrification flux consumption And gas loss occurs; for the distribution of gases produced by denitrification, it is preferable to impose distribution consistency constraints, for example: ; in, This represents the total denitrification flux. for Allocation ratio (subject to boundary constraints S5.2 and S5.3). Further constraints may be imposed if necessary. The allocation should satisfy a complementary relationship to ensure that the gas allocation does not exceed the total denitrification. For phosphorus, the adsorption-desorption exchange between the active phosphorus pool and the adsorption buffer pool should be constrained by capacity and parameter ranges. Preferably, a feasible region restriction is imposed on the exchange process to ensure... The direction of the exchange flux is consistent with the relationship between the coulomb difference or isotherm, thus avoiding non-physical behaviors such as negative adsorption or instantaneous infinite desorption.

[0084] The aforementioned process consistency and coupling constraints can be selectively enabled based on the model structure and available observations. They can be enforced as hard constraints when data support is sufficient, and embedded as soft constraints in the assimilation solution in the form of penalty terms when representativeness errors are large or observations are sparse. By introducing substrate constraints, carrier constraints, and morphological transformation consistency constraints, non-physical compensation between different fluxes in assimilation and inversion can be significantly suppressed, improving the mechanistic consistency and scenario extrapolation reliability of the full flux closure results.

[0085] S5.5 Constraint Enforcement Operator and Closure Verification Output; To ensure the enforceability of the conservation constraints and feasible regions or coupled constraints defined in S5.1–S5.4, this invention further defines constraint enforcement operators and closure verification output interfaces to perform constraint correction on candidate solutions obtained in each assimilation window (or each update), and to form auditable closure and consistency diagnostic quantities. The constraint enforcement operator can be formally represented as a mapping from candidate solutions to feasible solutions: ; in, Candidate solutions generated by assimilation update or forward push. The set of constraints that are satisfied after the constraints are implemented. (or satisfy within the tolerance range) feasible solutions; and simultaneously output the conserved closure residuals, the degree of failure of the feasible region and the classification index for subsequent quality control, verifiable record keeping and scenario simulation access.

[0086] (1) Projection-type constraint implementation operator (hard constraint preferred). For nonnegativity and boundary constraints (S5.2) and parameter feasible region constraints (S5.3), the feasible region projection method is preferred to achieve fast correction, that is, to map the candidate solution components one by one into their allowable intervals. For example, for any component v, the following is satisfied: ; Where v can be inventory, throughput, parameters, or proportionality coefficients, etc.; constraints on proportionality or allocation coefficients (e.g.) and The preferred approach is to use simplex projection or normalized projection, ensuring that both interval and sum constraints are satisfied. For terms in the process coupling constraints (S5.4) that have a clear upper bound relationship (…),… Sequential projection or group projection strategies can be used, prioritizing the correction of the decisive upper bound ( ). , ), and then correct the controlled quantity ( This is to avoid projection conflicts and maintain mechanistic consistency.

[0087] (2) Penalty function or augmented Lagrangian constraint implementation operator (soft or hard switchable). For the matter conservation closed equation constraint (S5.1) and the partially coupled consistency constraint (S5.4), this invention can use the constraint residual as an additional term of the objective function, and use a penalty function or augmented Lagrangian form to achieve constraint satisfaction. Preferably, the equation constraint residual is written as: ; and incorporate it as a penalty term into the assimilation objective (or into the update correction process), for example, using a quadratic penalty function. Apply a penalty for breach of inequality constraints. ,in When a "hard constraint" effect is required, the residuals can be converged to within the preset tolerance by increasing the penalty coefficient or using augmented Lagrange multipliers for updating. When the representativeness error of the observations is large or the model structure is simplified, a smaller penalty coefficient can be used as a soft constraint to improve robustness and avoid infeasible solutions.

[0088] (3) Variable transformation constraint implementation and automatic satisfaction methods. For strictly non-negative variables (e.g., inventory, throughput, positive parameters), logarithmic or softplus transformations can be used; for proportional variables, softmax transformations can be used, so that the updated variable is unconstrained in the transformation space and automatically satisfies the feasible region requirements in the original space. For example, let... To ensure ,make To ensure and This method can be used in conjunction with projection or penalty functions to reduce projection conflicts and improve numerical stability.

[0089] (4) Closure Verification Output and Quality Grading. This invention outputs a unified closure verification index system after constraint implementation to determine whether the closure and consistency of each accounting unit and each window meet the admission requirements. Preferred outputs include: quality conservation closure residuals. (Including absolute and relative residuals, such as normalized by total input within the window or inventory scale), inequality constraint default rate (Including maximum default rate and weighted default rate), and dedicated diagnostic quantities for process consistency constraints (such as unit carrier output intensity, whether the allocation coefficient exceeds the limit, whether the substrate upper limit is touched, etc.). Furthermore, the results can be graded according to preset thresholds (e.g., "pass, warning, or fail" or "A, B, C grades"), and the triggering reasons, involved variables, and correction ranges can be recorded to form a "verifiable" chain of evidence and engineering traceability.

[0090] By implementing the above-mentioned constraint operators and closure verification output, this invention realizes a closed loop from "constraint definition" to "constraint execution, diagnosis, and traceability," enabling the assimilation results to maintain conservation closure and consistency with the mechanism while satisfying observational statistical fitting. This provides a stable and auditable constraint basis for subsequent joint assimilation updates, scenario difference accounting, and uncertainty assessment.

[0091] Furthermore, the process of performing constraint-based joint data assimilation and solution includes: Generate an initialization set of joint assimilation variables containing states and parameters based on the prior distribution; The forward prediction interface is used to advance the set members to generate prior states and prior total flux. Based on the observation operator library and unified error statistics framework, multi-source observations are fused to update the prior set and obtain candidate posterior solutions; By invoking the constraint set and its constraint implementation operators, the candidate posterior solutions are corrected to obtain feasible posterior solutions that satisfy conservation closure and physical feasibility. Perform set statistics on feasible posterior solutions to quantify posterior uncertainty and perform attribution decomposition. Perform consistency diagnosis and quality classification on the posterior results, and output a unified posterior benchmark data package with admission criteria.

[0092] Furthermore, this embodiment constrains the joint data assimilation solution and posterior full throughput output (performs assimilation and includes a built-in quality control admission module); This step involves multi-source observation. Observation operator library With error statistics description ( Based on (Q or set perturbation), according to the assimilation window The process of "prediction-update-constraint correction" is repeated to achieve nitrogen and phosphorus status. ,parameter With total flux Joint inversion. First, using a mechanistic model. With drive input Forward propagation yields prior predictions ( Subsequently, multi-source observations were integrated to form an innovation quantity and updated to obtain candidate posteriors. .

[0093] To ensure that the results satisfy the laws of matter conservation, flux closure, and physical feasible region, this step calls the constraint set after each update. and constraint enforcement operators Correct the candidate solutions and output the feasible posterior. Simultaneously generate closed residuals With default As a verifiable diagnostic metric. The system has a built-in quality control access module: comprehensive. The window and spatial units are graded for quality using key fit or stability indices, providing an admission judgment of "enterable scenario simulation, cautious, or not enterable". Finally, a unified posterior baseline data package (state, parameters, full throughput time series + uncertainty + quality control level) is output as a reference for subsequent scenario difference potential calculation.

[0094] S6.1 initialization; This step generates the prior starting point for the assimilation solution, completing the initial binding of observations, error models, and constraint systems. Preferably, the joint assimilation variables are defined as augmented vectors: ; in, This is the augmented state vector (S3.2). This is the process parameter vector (S3.3). For managing input or observation system bias correction terms (optional). The driving input is denoted as... Observations are recorded as The observation operators and error models are indexed by data source as follows: .

[0095] Based on the prior distribution and perturbation rules of S2.3, from The sampled set prior is denoted as . The set members are: ; For parameters that are strictly positive and span orders of magnitude, sampling in logarithmic space is preferred; for proportional or efficiency parameters, constraints in [0,1] are preferred; for management bias multiplicative factors, a cutoff interval centered at 1 is preferred, consistent with the parameter feasible region in S5.3. Optionally, hierarchical perturbations or shared bias term perturbations are applied to key drivers such as fertilizer application rate, manure effectiveness, and irrigation input, so that the ensemble prior reflects an uncertainty representation equivalent to Q.

[0096] After initialization, it is preferable to perform feasible region preprocessing once for each set member and call the constraint set. Projection or transformation operators (S5.5), for example: ; This ensures that basic constraints such as non-negativity, boundary conditions, and upper and lower bounds of parameters are satisfied at the initial stage.

[0097] Simultaneously load the observation package Observation operators corresponding to each observation source With error covariance And set the window length Together with the observation assembly method (window batch or sequential), it forms an initialization configuration that can directly enter the S6.2 forward prediction.

[0098] S6.2 Forward prediction (by Advancement is based on prior knowledge and )); This step is used in the assimilation window. Internally, based on drive input Mechanistic advancement of the ensemble priors generates prior states and prior total fluxes that are aligned with observations, providing a comparable prior background field for subsequent multi-source observation updates. Specifically, for each ensemble member... (Optional) ), using mechanistic models Prior predictions are obtained by moving forward within the window: ; Among them, the superscript " "Indicates prior knowledge, This represents process noise or structural error disturbance (set according to the rules in S2.2; if Q is represented using a set approach, then...). (Given by the corresponding perturbation term). If If set as a static parameter, it remains unchanged within the window; if gradual change is allowed, its prior evolution can be updated according to preset random walk or small perturbation rules during advancement.

[0099] The model is advanced by synchronously outputting or diagnosing the total flux set within the window. The fluxes are then converted into observation-consistent statistics (e.g., windowed cumulative or windowed average) according to the time alignment rules in S1.3. The total flux includes at least the input flux, the main internal transformation fluxes, and the port output fluxes (surface water, groundwater, atmosphere, harvest removal, etc.), and can be decomposed by soil layer or subsystem to support subsequent conservation closure checks (S5.1) and process consistency constraints (S5.4).

[0100] To ensure that the prior predictions meet basic physical feasibility, this step may optionally perform a rapid feasibility check on the results: when obvious non-physical values ​​appear (e.g., negative inventory, key ratio exceeding limits, exceeding boundary ranges, etc.), it is preferable to first perform feasible region projection or variable transformation (S5.5) to... Pull back into the feasible region and record the correction markers; this process is only used to ensure that the forward prediction can continue into the observation update, and does not replace the subsequent constraint execution in S6.4. After completing the window advance, ensemble priors are used. This serves as the background information input for the S6.3 multi-source observation fusion update.

[0101] S6.3 Multi-source observation fusion update (by...) Update to obtain candidate posterior ); This step is used to move the window The available multi-source observations are introduced into the assimilation system to update the prior set under unified error statistical weights, yielding candidate posteriors (solutions before applying conservation or feasible region corrections). Specifically, based on the observation index and metadata in S1.5, observation vectors are assembled within a window. (Can be used for concatenating multiple variables or grouping by data source or variable), and binds corresponding observation operators to each type of observation. Covariance of observation error .in Map prior states or parameters to comparable quantities consistent with the observation scope (such as mapping from points to cells, window accumulation or averaging, flow weighted averaging, layer weighted conversion, etc.). The error model and quality level rules of S2.1 determine the weighting (the variance of interpolated, underrepresented or outlier observations is automatically increased to reduce their weighting).

[0102] For each set member, first calculate the prior simulation value of the observation space: ; And generate innovation (residual): ; Subsequently, an update is performed in the set space to obtain candidate posterior augmented variables. (Optional) The update method can be implemented using ensemble Kalman filtering or equivalent incremental maximum a posteriori estimation. Without explicitly expanding the algorithm details, this step emphasizes that: through the state-parameter cross term in the ensemble covariance, observation information can be propagated from the observable state to the relevant parameters, realizing the joint update of state and parameters, thereby improving the consistency and transferability of cross-flux and cross-scale inversion.

[0103] To improve the stability and identifiability of multi-source fusion, this step preferably supports two update organization methods: one is "window batch update," which updates observations within a window all at once; the other is "group or sequential update," which updates observations sequentially by observation type or port (e.g., remote sensing—soil—water environment—gas), and performs a lightweight consistency check after each group update to avoid non-physical drift caused by single-type observation dominance. For observations with obvious anomalies or significant differences from prior knowledge, this step may introduce an innovation gating strategy: when the standardization innovation exceeds a threshold, the gating value of that observation is increased. Alternatively, it could be temporarily suspended from updating, thereby enhancing robustness to spike events and representativeness errors.

[0104] After completing this step, the candidate posterior is obtained. and their corresponding candidate flux (Based on the throughput diagnostic module) (Calculated or output synchronously by the model), the candidate posterior will undergo unified correction for conserved closure and feasible region or coupling constraints in S6.4 to form the final feasible posterior.

[0105] S6.4 Constraint Enforcement and Closure Correction (Call) Obtaining a feasible posterior ); This step transforms the candidate posterior solutions obtained in S6.3 into consistent solutions that satisfy conservation closure and the physical feasible region, and is a key step in the "conservation constraint total flux assimilation" of this invention. Specifically, for each set member's candidate posterior; (Optional) ) and its candidate flux Call the constraint set constructed by S5 With constraint enforcement operator Perform constraint corrections and output feasible posteriors: ; And based on this, (Optional) ).

[0106] Constraint execution is preferably performed in the order of "feasible region first, then closure, then coupling" to reduce constraint conflicts and maintain mechanism consistency: First, the nonnegativity and boundary constraints of S5.2 and the parameter feasible region constraints of S5.3 are projected or transformed to bring inventory, flux, and parameters back to the allowable range; then, the nitrogen and phosphorus mass conservation and total flux closure equation constraints of S5.1 are closed to ensure that the relationship "inventory change = total input - total output" within the window is valid within the preset tolerance; finally, the process consistency and coupling constraints of S5.4 (substrate constraints, carrier constraints, morphological transformation consistency, etc.) are consistent to perform consistency correction or penalty convergence to suppress flux compensation and sign mismatch.

[0107] To achieve "verifiable and traceable" results, this step simultaneously calculates and outputs the closure residual and the default rate index. For equality constraints, the conserved residual is calculated: ; It can also output normalized residuals (e.g., normalized by total input within the window or inventory scale); and calculate the default rate for inequality constraints. ; It also records the maximum default rate and the weighted default rate. and It serves as the basis for quality grading and admission judgment in subsequent S6.6, and also provides diagnostic clues for identifying problems such as "inconsistent observation caliber, missing model structure, or excessive input deviation".

[0108] When constraint correction results in excessive adjustment (e.g., the correction magnitude exceeds a preset threshold or the residual fails to converge), this step may optionally trigger a rollback or weight reduction mechanism: including increasing the weight of relevant observation sources. Increase the representation of process noise or shrink the parameter update amplitude to avoid the diffusion of non-physical updates driven by a single observation or scale mismatch. The result after constraint execution... This is the posterior result that satisfies conservation closure and physical consistency, and serves as the input for uncertainty quantification and attribution decomposition in S6.5.

[0109] S6.5 Quantification and Attribution of Posterior Uncertainty; This step is used to obtain posterior results that satisfy the constraints. Subsequently, a unified uncertainty characterization is given for key states, parameters, and fluxes, and the sources of uncertainty are decomposed into three categories: input uncertainty, parameter uncertainty, and structural and process uncertainty. This improves the credibility and interpretability of the calculation conclusions and provides comparable confidence intervals for subsequent scenario-differential potential calculations. Uncertainty quantification is mainly based on set statistics: for any output quantity q (which can be a component of x, ... Components or f-components and their summary indices), by set It calculates the posterior mean, variance, and quantile intervals (e.g., 5%–95%), and can further output uncertainty distribution maps on spatial units and time windows to reflect the spatial and temporal heterogeneity of uncertainty.

[0110] To facilitate engineering accounting, this step prioritizes the unified output and dissemination of "accounting-sensitive" summary indicators, including but not limited to: total inputs and outputs within the window, port outputs (surface water, groundwater, atmosphere, harvest removal), and key loss fluxes (e.g., This includes crop absorption and yield-related derivatives, as well as nitrogen and phosphorus budget deficits and closure diagnostics (e.g., normalized residuals). For indicators that need to be aggregated across scales (from grids to watersheds or administrative regions), it is preferable to use area-weighted or carrier-weighted rules consistent with S1.3 for aggregation, and simultaneously propagate uncertainties (e.g., aggregate set members first and then perform statistics) to ensure that scale transformation does not introduce additional caliber errors.

[0111] Regarding uncertainty attribution, this step preferably uses the concepts of "factor set set" or "variance decomposition" to break down the sources of uncertainty: the set disturbances are divided into input or driving disturbance sets according to their sources (corresponding to...). With optional deviation items ), parameter disturbance group (corresponding) Prior and update uncertainties, as well as structural and process disturbances (corresponding to process noise or equivalent representations of model structural simplification), are propagated and their output variances calculated while keeping other disturbance groups constant, thus obtaining the contribution ratio of each source to the total uncertainty. For key outputs, sensitivity measures (such as ensemble regression or local linear approximation) can be further combined to give a ranking of "major contributing factors," which explains why uncertainty is higher in certain regions or windows and provides a priority basis for subsequent data supplementation and model improvement.

[0112] The output of this step is a posterior uncertainty data package, which includes confidence intervals for key outputs, spatial-temporal uncertainty distributions, and source contribution decomposition results. Together with the quality control grading in S6.6, it constitutes a posterior benchmark quality description that can be used for scenario extrapolation, thereby ensuring that subsequent scenario differencing not only provides point estimates but also provides potential intervals and explanations of their uncertainty sources.

[0113] S6.6 Consistency Diagnosis and Quality Grading; This step is used to perform consistency diagnosis and quality classification on the posterior results of each assimilation window and accounting unit, and to provide an admission judgment of "entering scenario simulation, requiring caution, or not entering," to ensure that subsequent scenario difference potential accounting is based on a closed, verifiable, mechanistically consistent, and numerically stable posterior benchmark. The quality control admission threshold is expressed as the conserved residuals output by S6.4. With default Based on the core criteria, and combined with the observation fit consistency and parameter stability index, a comprehensive judgment is formed.

[0114] First, a consistency index between conservation closure and feasible region is constructed. For equality constraints, the normalized closure residual is used as the main index (e.g., using...). Normalize by total input or inventory scale within the window, and set a threshold. Determine whether the closure is successful; for inequality constraints, use the maximum default rate and weighted default rate as the main indicators, and set a threshold. Determine whether the feasible region and coupling constraints pass. If persistent exceedances occur, prioritize outputting port decomposition or process decomposition markers of residuals or default values ​​(e.g., surface water port, groundwater port, or atmospheric port contribute more) to pinpoint the source of the problem.

[0115] Secondly, an observation fit consistency index is constructed. Based on the innovation statistics of each observation source within the window (such as standardized innovation, weighted residuals, or equivalent goodness-of-fit index), it is checked whether there is systematic bias or a single observation source dominating the update; when a certain type of observation shows a long-term same-sign bias or significantly large standardized innovation, it is preferable to output a mark indicating that "observation caliber, observation operator, and error variance may be mismatched", providing a basis for subsequent rollback or parameterization adjustment.

[0116] Next, construct indices for parameter portability and numerical stability. (Check...) Does the boundary frequently touch (S5.3)? Are there unreasonable window jumps or spatial jaggedness? And what are the key deviation items? Whether there is an unbounded drift tendency; if the relevant threshold is triggered, it is preferable to mark the window or cell as "caution required" and indicate that it may be necessary to reduce the degree of freedom of parameter update, increase process noise, or strengthen partition sharing or smoothing constraints.

[0117] Based on the above three types of indicators, this step establishes quality grading and admission rules. A preferred definition is a three-level judgment: when... When both the fit and parameter stability indices meet the thresholds, the result is judged as "Pass (can proceed to scenario simulation)". When only a slight exceedance of the threshold occurs or a single warning exists, the result is judged as "Warning (use with caution, can proceed but with uncertainty and marking)". When the closure residuals or default rate significantly exceed the threshold, or when there is obvious systematic fitting distortion or parameter drift, the result is judged as "Fail (cannot proceed to scenario simulation)". All judgment results should output the corresponding triggering reasons, involved variables, ports, observation sources and their numerical evidence (e.g., residuals, default rates, standardization innovations), and form an auditable and traceable record.

[0118] When a result is deemed "failed" or a series of "warnings" appear, this step may trigger a rollback or deweighting strategy: for example, increasing the weight of relevant observation sources. (Reducing its weight), gating abnormal observations, increasing process noise characterization to absorb structural errors, or shrinking parameter update amplitude and strengthening the parameter feasible region or smoothing constraints; after rollback, S6.3–S6.4 can be re-executed to obtain a posterior benchmark that meets the admission requirements. By making closure and consistency determination the quality control admission threshold in advance, this invention realizes closed-loop control of "assimilation-constraint-verification", ensuring that the benchmark for subsequent scenario difference potential calculation has reliability, interpretability and engineering verifiability.

[0119] Furthermore, the process of performing scenario differential potential assessment includes: Starting from a unified posterior baseline data package, baseline scenarios and action scenarios are constructed. The action scenarios are implemented by modifying the driving inputs. Under the same model, constraint set and initial state, forward extrapolation is performed on the baseline scenario and the action scenario respectively, and the respective state trajectory and full flux ledger are output; By pairing and differencing members of the same posterior set, the flux difference and state difference between the action scenario and the baseline scenario are calculated, and a potential accounting index system is constructed. The posterior uncertainty is propagated in the scenario difference, the difference results are aggregated and statistically analyzed, the potential confidence interval and robustness probability are output, and the risk screening and feasibility assessment are combined with engineering constraints. The differential potential results of the accounting units are spatially aggregated according to the specified partition boundaries, and the statistical and priority calculation results at the partition scale are output.

[0120] Furthermore, the construction of the action scenarios includes: by adjusting the amount of fertilizer applied, the timing of fertilizer application, the method of fertilizer application, the form of fertilizer, the irrigation system, and the management of organic materials in the driving input, to obtain nitrogen control measures, phosphorus control measures, water and transport control measures, and organic fertilizer and straw management measures with different intensity levels.

[0121] Furthermore, this embodiment unifies the construction of counterfactual baselines and action scenarios under a unified posterior benchmark, as well as the differential potential calculation and uncertainty output. Specifically, using the "unified posterior baseline" output by S6 as a common starting point, counterfactual baseline scenarios and action scenarios are constructed under a consistent model structure, parameter system, and constraint system. Potential results such as emission reduction, loss reduction, and efficiency improvement are obtained through scenario differencing. Simultaneously, posterior uncertainties are propagated and quantified during scenario extrapolation to ensure that the results of each scenario are comparable, interpretable, and verifiable. Inputs include: the unified posterior baseline data package from S6 (…). and Uncertainty and quality control level) and the set of action plans to be evaluated (by adjusting the driving inputs) Some parameters (Or implemented through process settings). Preferably, scenario differential calculations are only performed on spatial units or time windows that pass quality control gating. The "caution required and fail" sections are marked and can be optionally excluded from potential summaries or reported separately.

[0122] More specifically, it includes the following steps: S7.1 Building a Baseline Scenario Set of scenarios for measures ; This step uses the unified posterior benchmark output by S6 as a common starting point, while maintaining the model structure. Constraint Sets (S5) Under the premise that the posterior state-parameter information is consistent, construct the baseline scenario. Set of scenarios for measures This ensures that the scenario difference has a "counterfactual" meaning and strict comparability. For consistency with the notation above, the joint assimilation variable is denoted as... (Optional) The driving input sequence is denoted as Multi-source observation is The size of the set is .

[0123] Starting with a unified posterior support, for each set member All scenarios share the same posterior initial values ​​and parameter set, that is: ; in, This is the starting point for scenario simulation (e.g., the start of a seasonal or annual window). This setting ensures that the "difference" is caused only by changes in the scenario inputs, and not by incomparable errors due to different initial values ​​or parameters.

[0124] Baseline scenario : Represents the current management and background conditions during the study period, using a baseline input sequence consistent with assimilation as the driving force, denoted as: ; in It can be directly taken from the management input library built by S1 (and can use the input version estimated or corrected during the assimilation process or with...) (A consistent revised definition) should include at least fertilization, organic fertilizer input, straw return to the field, irrigation volume and concentration, sedimentation, and other management and external inputs; external drivers such as meteorology and hydrology should also be considered. Components or parallel driving inputs. Model structure and constraint set in the baseline scenario. Keep it unchanged to ensure output flux ledger closure and physical feasibility.

[0125] Measure Scenario Without altering the external natural driving parameters (preferably meteorological, hydrological, etc., to maintain consistency with the baseline) and sharing the same posterior initial values ​​and parameter set, only the sub-vectors related to the measures in the input vector are replaced or adjusted to obtain the measure input sequence: ; in, This is a measure mapping operator used to transform "measure schemes" into pairs of measures. Structural modifications (such as adjusting fertilizer application rate, timing of application, application method, changing fertilizer form ratio, changing manure effectiveness coefficient, changing irrigation system or drainage management, introducing cover crops and field engineering measures, etc.).

[0126] The key constraint is: for the same set member i, only between scenarios... The difference is: ; This ensures the comparability of scenario differences and the semantics of "counterfactual".

[0127] Classification and Intensity Levels of Implementation Scenarios: Implementation Scenarios Scenario families can be organized according to the controlled object and the action path, and their functions can be explicitly described. Implementation method in: (1) Nitrogen control measures: total nitrogen application reduction, multiple nitrogen applications, deep application, side-deep application, controlled-release fertilizer, slow-release fertilizer, nitrification inhibitors, optimization of fertilization timing, etc. (corresponding to) (Adjustment of nitrogen input intensity, time allocation and speciation). (2) Phosphorus control measures: phosphorus reduction and substitution, precise phosphorus application, formulation optimization, regulation of soil phosphorus activity and buffering (corresponding to Adjustments to phosphorus input intensity, application method, and related management factors (for medium phosphorus input). (3) Water and transport control: controlled irrigation, intermittent irrigation, drainage management, field engineering measures (buffer zones, intercepting ditches, ecological ditches, etc.), and cover crops to reduce runoff erosion (corresponding to...) Adjustments to irrigation systems, drainage parameter inputs, and management events. (4) Organic fertilizer and straw management: manure substitution for chemical fertilizer, composting and optimization of application timing, optimization of straw return to the field, etc. (corresponding to) Adjustments to the amount, effectiveness, and timing of organic inputs.

[0128] Ideally, each type of measure should have intensity levels (e.g., 10%, 20%, 30% reduction, or different coverage ratios or engineering interception efficiency levels) to form a comparable family of scenarios. And record the modifications for each scenario. The components, time windows, and applicable conditions serve as metadata inputs for subsequent differential accounting and verifiable record keeping.

[0129] S7.2 Scenario Simulation Execution; This step is used in the same mechanistic model. Same set of constraints Under condition (S5), for the baseline scenario Scenarios for each measure A forward inference is performed, outputting the state trajectories and full flux ledger for each scenario, providing rigorously comparable scenario results for subsequent differential potential calculations (S7.3 and S7.4). The inference uses the unified posterior set of S6 as a common starting point: for each set member... Using the same initial posterior (Optional) Only change the context input sequence and This ensures that the difference has counterfactual semantics.

[0130] In the window Within, for any scenario Perform set deduction according to the model advancement operator: ; Among them, the superscript " "This indicates the predicted state obtained from scenario deduction." For process noise or structural error perturbations (if ensemble perturbations are used to characterize process uncertainty, the same set of perturbations or the same random seed strategy should be used in all scenarios to reduce differential noise and maintain comparability). If parameters If set as static posterior parameters, then each scenario remains unchanged. If gradual change is allowed, all scenarios will proceed using the same gradual change rules, and their evolutionary mechanisms will not change due to changes in scenarios.

[0131] The full throughput ledger within the window is obtained synchronously during model advancement or through diagnostics. The data is then converted into statistics (window cumulative or average, flow weighted, etc.) consistent with the accounting caliber according to the alignment rules in S1.3. The ledger should at least cover: input fluxes, key internal transformation fluxes, and port output fluxes (surface water, groundwater, atmosphere, harvest removal, etc.), and can be decomposed by soil layer or subsystem to support subsequent conservation closure verification and differential interpretation.

[0132] To ensure that the scenario simulation meets the requirements of material conservation and physical feasibility throughout, this step invokes constraint verification and necessary correction mechanisms after each window or key event: for situations such as nonnegativity, boundary out-of-bounds, parameter out-of-bounds, and conservation closure residual exceeding the threshold, feasible domain projection, variable transformation, or closure correction is preferably performed, and the correction magnitude and marking are recorded; if a certain measure scenario continuously triggers non-convergence default or causes a significant increase in closure residual in a certain accounting unit, then the "scenario-unit-window" combination is marked as infeasible (for use in S7.5 risk screening), avoiding the inclusion of infeasible simulation results in the potential summary.

[0133] After completing this step, output each scenario. and Collective results during the study period: and It also carries basic information on quality control indicators (closure residuals, default rate, feasibility markers) and uncertainty statistics; these results will serve as direct inputs for subsequent scenario difference potential calculations and uncertainty propagation.

[0134] S7.3 Differential Potential Calculation ( ) and the construction of the indicator system; This step, after completing each scenario simulation (S7.2), calculates the scenario difference between the "baseline" and the "measures" under the condition of one-to-one correspondence among members of the same posterior set. This yields potential results such as nitrogen and phosphorus emission reduction, loss reduction, and efficiency improvement, and organizes the difference results into a deliverable indicator system. To ensure counterfactual comparability, this step strictly employs "paired differencing within the same member": for each set member i, it is paired with each accounting unit and window. Baseline scenario Context of Measures Using the same initial values Same parameters (Same option) And consistent process perturbation is achieved, only This allows the difference components to have a clear causal explanation.

[0135] (1) Definition and symbol conventions of difference components; For any output variable q (which can be an element of the state component x, a flux of the flux ledger f, or an index derived therefrom), define the action scenario. The difference potential is: ; in, The measures indicate that the measures result in a "reduction or decrease" relative to the baseline (e.g., loss of flux, port output, emissions, etc.). If the indicator is "benefit-oriented" (e.g., crop absorption, yield, utilization efficiency), the same difference sign can be explicitly used in the indicator definition (e.g., "increase" is positive), or a unified sign explanation can be given in the output to avoid confusion between the positive and negative meanings of different indicators.

[0136] (2) Total flux differential ( The accounting objects of ) Ideally, the differential accounting object should cover "port output + key loss + key transformation + input change" to form a closed interpretable chain. It should include at least: Differential port outputs: surface water output, groundwater output, atmospheric exchange (e.g.) , The difference between port fluxes such as harvest and removal; Key loss differentials: nitrogen leaching, leakage, runoff, erosion, volatilization and gas emissions, etc.; dissolved phosphorus (DRP) and particulate phosphorus (PP) emissions, etc. Key process differentials (optional): Flux differentials for processes such as nitrification, denitrification, mineralization-retention, adsorption-desorption, and sedimentation-resuspension, used to explain why port outputs change; Input-side differential: Changes in management inputs such as fertilizer or organic fertilizer input, irrigation input, and sedimentation input between the baseline and the intervention (directly determined by...). (Obtained), used to establish a correspondence between "measure intensity" and "response magnitude".

[0137] The flux of the cumulative caliber of the window, Take the cumulative value within the window; for concentration or proportion outputs, use window averaging or carrier weighted averaging to form comparable quantities according to the rules in S1.3.

[0138] (3) State difference ( ) and inventory response; In addition to flux differential, the preferred output is the differential response of key inventory or status to identify whether "emission reduction or loss reduction comes at the cost of inventory accumulation" and "soil nitrogen and phosphorus speciation migration." For example: soil inorganic nitrogen pool ( The differences in the stockpiles of organic nitrogen, active phosphorus, adsorbed phosphorus, and organic phosphorus, as well as the differences in the state of nitrogen and phosphorus absorption by crops, should be included. The accounting scope for the state differences should be consistent with the conservation closed window in S5.1 to form a closed-loop interpretation with the total flux differences.

[0139] (4) Construction of the indicator system; To facilitate engineering accounting and policy evaluation, the differential results are preferably organized into a "three-tiered indicator system": Target-oriented indicators: Targeting specific ports or environmental objects (e.g., reduction in surface water TN and TP output, reduction in groundwater nitrate load, atmospheric...). (Emission reduction).

[0140] Path-based indicators: These are process-path oriented indicators (e.g., leaching loss, runoff reduction, erosion output reduction, volatilization reduction, etc.) used to explain the sources of change in object-based indicators.

[0141] Comprehensive indicators: used for cross-regional comparisons and rankings, such as potential per unit area (kg N ha). -1 ,kg P ha -1 ), unit yield potential (kg N t) -1 ,kg P t -1 Changes in intensity indicators (pollution intensity, emission intensity), improvements in utilization efficiency (NUE, PUE), and the "marginal emission reduction effect of changes in unit input" ).

[0142] Meanwhile, to meet the "feasibility" screening requirements, the indicator system should include at least one type of constraint-related indicator (such as whether production reduction is not required, whether the closure residual is passed, and whether the default rate is zero or below the threshold) for subsequent risk screening and access judgment in S7.4 and S7.5.

[0143] (5) Quality labeling and usability output of the difference results; For each "measure k - accounting unit - window", this step should inherit the feasibility marking and quality control indicators from S7.2: if the combination is determined to be infeasible (e.g., frequent defaults, closed residual exceeding the threshold, leading to obvious non-physical states), the difference results should be marked as "unavailable" or "requires caution" and removed or listed separately during spatial aggregation. For combinations that pass quality control, the difference results proceed to the subsequent uncertainty propagation and risk screening steps.

[0144] S7.4 Uncertainty Transmission, Risk Screening and Feasibility Assessment; This step is used to consistently propagate the posterior uncertainty of S6 through scenario simulation and differential accounting, and to perform risk screening and feasibility assessment for each "measure k – accounting unit – window" combination, outputting a comparable and implementable potential result (including confidence intervals and pass / fail indicators). The input to this step is the differential set from S7.3. (where q can be) , (or derived indicators), as well as quality control indicators (closure residuals, default rate, observation consistency, parameter stability, etc.) and engineering constraints (no production reduction, threshold, boundary conditions, etc.) generated by S7.2 and S6.6. The output is: the expected value and quantile interval of the potential, the robustness probability, and the judgment result of "can be included in the summary, require caution, or is not available".

[0145] (1) Difference uncertainty statistics; For any indicator q and any action scenario Based on set pairing difference To calculate the statistic for the posterior difference distribution, at least the following should be included: ; in For quantile functions (e.g.) (Corresponding to the 5%–95% range). When outputting the aggregation potential across windows and spatial units, it is preferable to adopt the principle of "aggregating each set member first, then performing statistics" to avoid underestimating uncertainty (i.e., Then perform statistics on i.

[0146] (2) Robust potential and the "efficiency probability" indicator; To avoid judging potential based solely on point estimates, robustness indicators are preferred. For "emission reduction or loss reduction type" indicators, exceeding a threshold is defined. The probability is: ; in, For indicator functions, It can be set to 0 (representing "net benefit probability") or an engineering threshold (representing "probability of achieving the target"). For "benefit-oriented" metrics (such as NUE improvement or output change), a consistent threshold probability is also output (e.g., ...). (This indicates the probability of no production cut).

[0147] (3) Feasibility constraints and risk screening rules; This invention preferably employs a "two-layer screening" mechanism: first, hard constraints are used to determine feasibility; then, within the feasible set, the results are ranked according to robustness and potential. Hard constraints include at least: Conservative closure and physical feasibility: Inherit and reuse the scenario deduction quality control of S7.2, requiring that the closure residuals and default degree do not exceed the threshold or that there is no persistent default; No production cuts or no harm constraints (optional but recommended): e.g., requirements (e.g., 0.8 or other thresholds), or require that the lower quantile of output not be lower than the baseline; Engineering thresholds and boundary conditions: for example, fertilizer reduction does not exceed the feasible upper limit, coverage ratio or engineering efficiency is within a reasonable range, parameters do not touch the boundary and have no drift (consistent with S5.3). Side effect constraints (optional): for example, nitrogen control leading to When rising or controlled water levels lead to increased groundwater risk, a "cooperative non-deterioration" threshold should be set or multi-objective screening should be used.

[0148] If a combination does not meet the hard constraints, it is marked as "unavailable or infeasible"; if it meets the hard constraints but lacks robustness (e.g., the interval crosses zero), it is marked as "unavailable or infeasible". Those with lower values ​​are marked as "caution required"; those that meet hard constraints and have high robustness are marked as "pass or can enter the summary".

[0149] (4) Overall score and recommendation ranking; To generate an actionable list of potential projects, a comprehensive score can be built within the set. For example, a weighted score considering "potential size + robustness + feasibility" can be used simultaneously: ; in, The difference index of measure k The posterior expectation or set mean; Norm(⋅) is a normalization function that normalizes values ​​of different dimensions or magnitudes. Convert to comparable dimensionless values ​​(usually 0–1); The probability that measure k will achieve "benefit or meet the target" (robustness); , , The importance weights of the three types of objectives; It can be constructed from factors such as interval width, number of defaults, and number of warnings triggered. This score is used to prioritize measures within the same area and outputs a combined list of "recommended measures - intensity level - applicable conditions".

[0150] (5) Output structure; This step outputs the following for each "Measure k – Accounting Unit – Window": Potential point estimates (mean or median) and uncertainty intervals (quantiles or confidence intervals); Robustness probability indicators (such as) , , ); Feasibility assessment labels (e.g., approved, require caution, unusable) and triggering reasons (e.g., closed residual, default degree, production constraints, threshold exceedance, etc.). Overall score and ranking results.

[0151] S7.5 Spatial aggregation, partition statistics and priority calculation; This step aggregates the differential potential results obtained in S7.3–S7.4 from the accounting unit scale to the watershed or administrative region scale, forming statistical and priority calculation results at the zoning scale. Before aggregation, the zoning boundary G and weights are determined. (e.g., area weight, crop area weight, yield weight, runoff contribution weight, etc.), and consistent with the aforementioned accounting methods.

[0152] During aggregation, feasibility filtering is preferred: only "measure-unit-window" combinations that meet the S7.4 "pass or can enter the summary" criteria are included in the main summary; combinations that "require caution or are not feasible" are marked and removed from the summary or counted separately, while outputting the pass rate and coverage rate (e.g., percentage of passed area, percentage of passed unit, etc.) to characterize the feasibility and stability of the measures at the zonal scale.

[0153] To maintain consistency in uncertainty, this step adopts the rule of "aggregating each set member first, then statistically analyzing the aggregation result interval": For any partition G and each set member i within measure k, calculate the partition aggregation difference: ; For aggregate indicators, weighted summation can be used; for intensity indicators, weighted average can be used, selected according to the indicator's scope. Then... Calculate the point estimate and quantile interval of the partition potential, and output the robustness probability (e.g.) Or, according to the threshold probability agreed upon by the indicator "direction of improvement", it is used to characterize the statistical significance and robustness of potential.

[0154] The final output is a statistical result of "measures – intensity level – scope of application" at the zoning scale (including point estimates, intervals, robustness probabilities, and summaries of pass rates or coverage rates), which can be used to form a priority ranking.

[0155] Furthermore, the output of the difference potential results, which include uncertainties, includes: Output the posterior state and the full-process accounting ledger of total flux; Output a list of scenario potentials and spatial distribution products for the differential potential results; Output probabilistic expressions and risk characterization indicators for posterior uncertainty and difference potential; Output the verification indicators and log records for conservation closure and consistency; Based on feasibility assessment and robustness probability, the results of optimal measures and comprehensive evaluation are output.

[0156] Compared with the prior art, the present invention has the following advantages and technical effects: Furthermore, this embodiment outputs results and a verifiable indicator system. After completing scenario differential deduction and potential calculation under a unified posterior benchmark (S7), this step standardizes, summarizes, structurally encapsulates, spatially aggregates, and verifies the posterior results and scenario differential results, forming a results output package and indicator system oriented towards management decision-making and engineering applications. The output follows the principle of "full throughput closure—traceability—verifiability," and includes at least: a full-process throughput accounting ledger, a scenario potential list and spatial products, uncertainty intervals and probabilistic expressions, conservation closure and consistency verification indicators, as well as a unified data structure and output format specification.

[0157] S8.1 Full-Process Throughput Accounting Ledger Output; This step is based on the posterior "state-parameter-flux" set obtained from S6, and is performed according to the accounting unit. With time window Output nitrogen and phosphorus full-process flux accounting ledger It is used to support "traceable accounting" and cross-regional comparison. The ledger records posterior statistics (such as set mean or quantiles) and necessary set member numbers i in a standardized field format, and decomposes and summarizes the throughput according to the path of "input-transformation-loss-output".

[0158] The ledger contains at least two types of entries: one is inventory or status-based post-hoc output, corresponding to... Key reservoirs (e.g., soil nitrogen multiforms) Soil phosphorus multiform pool Crop storage and key hydrological conditions, soil moisture content The first type of output includes surface runoff (R), infiltration (D), evapotranspiration (ET), etc.; the second type is flux-related a posteriori output, corresponding to the cumulative flux in the window. It should at least cover nitrogen input, mineralization-nitrification-denitrification, Volatilization The data includes emissions, leaching, seepage, runoff erosion, crop absorption and harvesting, as well as phosphorus input, adsorption-desorption exchange, dissolved phosphorus (DRP) and particulate matter (PP) runoff erosion and harvesting. Preferably, the ledger simultaneously outputs receptor decomposition results (including water receptors, atmospheric receptors, and agricultural product carryover) and intensity indicators (per unit area, per unit yield) for use in engineering evaluation.

[0159] S8.2 Scenario Potential List and Spatial Distribution Products; This step outputs the scenario difference results calculated in S7 as a unified expression: "Potential List + Spatial Distribution Product". For each action scenario... , with baseline For reference, output differential flux and differential state. It also includes path decomposition items, which are bound to scenario numbers, intensity levels, and applicable condition tags, thus forming a checklist-based delivery that can be directly used for decision-making and project implementation.

[0160] The list preferably includes both the accounting unit scale (corresponding to S7.3–S7.4) and the zonal scale summary results (corresponding to S7.5): at the zonal scale, it references the aggregation point estimate, quantile interval, and robustness probability (e.g., from the output of S7.5). (This includes) pass rate or coverage summaries to characterize the implementability and effectiveness stability across regions. Differential metrics can be derived by path decomposition, providing typical key terms and corresponding intensity expressions (e.g., per unit area, per unit output, or synergistic benefits). For example, output can be categorized by process or medium. , , , wait, This represents the difference in nitrogen leaching or leakage flux. express Emission flux differential This represents the difference in dissolved reactive phosphorus output flux in runoff. This represents the difference in phosphorus output flux from eroded particulate matter.

[0161] Spatial products should include at least: a potential spatial distribution layer, the results of hotspot or priority area division, and a graded expression layer under different measure intensity levels; and can simultaneously provide spatial distribution information such as "approved, requires caution, and not feasible" to support implementation screening and risk warnings.

[0162] S8.3 Quantitative Output and Probabilistic Expression of Uncertainty; This step solidifies the output of the posterior uncertainty and scenario difference uncertainty generated in S6–S7, avoiding the unverifiable and indecisive nature of providing only single-point conclusions. The uncertainty output is based on ensemble statistics and includes at least: the mean, standard deviation, and quantile range (e.g., 5%–95%) of key state quantities, key fluxes, and potential differences; and probability indicators for meeting engineering thresholds or policy objectives, such as: ; in Any differential flux, differential load, or composite index can be selected. Let Y be the target threshold, and Y be the output or a production proxy indicator. Preferably, uncertainty source hints (e.g., observations, inputs, parameters, structures) are also output as explanatory fields.

[0163] S8.4 Conservation Closure and Consistency Verification Indicators; To demonstrate the core innovation of this invention—"conservation constraint—total flux closure"—this step structures the output of the quality control threshold results and constraint execution information of S6, forming a verification indicator system and record-keeping mechanism that can be verified by a third party. The verification indicators include at least: nitrogen and phosphorus closure residuals (e.g., The log includes time-series and statistical summaries (mean, maximum, number of times exceeding the threshold, and relative closure error); constraint triggering and correction logs (including the number of non-negative or boundary projections, the number of parameter projections, the number of closure corrections, and the correction magnitude, recording the adjusted variable and the adjustment magnitude); observation fit consistency indicators (e.g., the mean or variance of innovation statistics, the proportion of exceeding the threshold) to identify overfitting or underfitting; and optional physical consistency check results (e.g., "output flux does not exceed available supply inventory" and "water quantity is consistent"). This log information is output as part of the result delivery to support audit traceability and recalculation verification.

[0164] S8.5 Optimization and Comprehensive Evaluation Output of Measures; This step, without repeating the scenario deduction and difference calculations of S7, summarizes and expresses the "passable, require caution, infeasible" judgments, robustness probabilities, and potential results already formed in S7, outputting a feasible solution set and recommended solutions. The output should include at least: a set of measures that meet the constraints (e.g., no production reduction, closure residuals meeting the target, risk thresholds met, etc.); ranking and grouping results according to dimensions such as total potential, potential per unit area, synergistic benefits, and risk reduction; and an optional summary of implementation recommendations (recommendation strength level, expected effect range, and key sources of uncertainty).

[0165] S8.6 Standardized Data Structures and Output Formats; To improve project usability, this step adopts a unified data structure and format specification for all outputs, including at least: spatial index (accounting unit ID, administrative region code, watershed code), time index (year, quarter, window number, start and end time), scenario number k, ensemble statistical fields (mean, quantiles, variance, probability), unit and caliber identifiers (cumulative, average, flow-weighted, etc.), and flux type identifiers (input, transformation, loss, harvest carryover, receptor decomposition). Preferably, the output medium includes a tabular list, raster or vector layers, log files, and metadata descriptions (model version, data version, constraint thresholds, and assimilation configuration) to ensure consistency and recalculation across batches.

[0166] In summary, this embodiment performs joint assimilation estimation of multi-source observations such as fertilization statistics, soil monitoring, remote sensing inversion, water quality, groundwater and greenhouse gases under a unified error statistics framework. It also introduces material conservation constraints and physical feasible region constraints, so that the fluxes of nitrogen and phosphorus throughout the entire process of "input-conversion-transport-output-sink change" can be closed and consistent within the same framework. This significantly reduces inconsistencies between fluxes and difficult-to-explain "gaps", and improves the physical rationality and internal consistency of the calculation results.

[0167] This embodiment fully utilizes the complementary information from multiple data sources, reducing the impact of single data source bias and model structure errors on the results, and improving the estimation accuracy and robustness of key state quantities and major losses or emission fluxes. Simultaneously, the assimilation process forms a unified set of posterior states, parameters, and fluxes, and can record constraint residuals and error convergence information, making the accounting chain "input-estimation-constraint-output" traceable, reproducible, and verifiable, facilitating third-party auditing and result review.

[0168] In terms of emission reduction potential assessment, this embodiment constructs counterfactual baselines and action scenarios using a unified posterior benchmark, and conducts scenario differential accounting under consistent initial states, parameter sets, and external driving settings. This effectively reduces the confounding effects caused by differences in climate, soil, and management backgrounds, thereby more reliably separating and attributing the effects of the measures, and enhancing the comparability, explanatory power, and attribution ability of the potential accounting results.

[0169] This embodiment can systematically quantify the combined impact of input errors, parameter uncertainties, and structural errors, outputting information such as confidence intervals, risk probabilities, or sensitivity contributions. This avoids decision-making biases caused by providing only point estimates and better aligns with the reliability and risk representation requirements of policy evaluation. Simultaneously, the method is compatible with multi-scale applications such as field plots, grids, watersheds, and administrative regions, facilitating cross-regional promotion and engineering deployment, and providing stable quantitative support for priority zoning, precise governance, and measure optimization. Standardized interface data packages are obtained based on multi-source data. The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for agricultural nitrogen and phosphorus total flux assimilation and scenario-difference potential calculation, characterized in that, include: Based on multi-source data, obtain standardized interface data packets; Based on the standardized interface data package, an assimilable agricultural nitrogen and phosphorus state-space model and forward prediction interface are constructed. Based on the standardized interface data package and the forward prediction interface, an observation operator library and a unified error statistics framework are established. Based on the aforementioned observation operator library and unified error statistics framework, information content diagnosis is performed and assimilation configuration is generated; Based on the aforementioned forward prediction interface, observation operator library, unified error statistics framework, and assimilation configuration, a constraint set of matter conservation constraints and physical feasible region constraints is constructed and implemented. Based on the forward prediction interface, observation operator library, unified error statistics framework, assimilation configuration and constraint set, perform constraint joint data assimilation solution to obtain unified posterior baseline data package; Based on the unified posterior baseline data package, a counterfactual baseline and action scenarios are constructed, scenario differential potential is calculated, and differential potential results containing uncertainty are output.

2. The method according to claim 1, characterized in that, The process of obtaining the standardized interface data packet includes: Determine the system boundaries, accounting units, and time base of the accounting objects; Collect and organize the driving input library, observation library, and prior library, and perform element coding and unit unification; Spatiotemporal scale matching and quality control are performed on the data in the driving input library and the observation library; The data that has undergone spatiotemporal scale matching and quality control is assembled into a standardized interface data package containing observations, driving inputs, quality labels, and metadata.

3. The method according to claim 1, characterized in that, The process of constructing an assimilable state-space model of agricultural nitrogen and phosphorus and a forward prediction interface includes: Construct an extended state vector to describe key nitrogen and phosphorus pools, crop nutrient status, and optional deviation states; Define a full flux vector that explicitly decomposes the input, transformation, transfer, loss, and output processes into traceable components; Define the driving inputs as the entry point for management measures, the static covariates characterizing spatial heterogeneity, and the set of parameters to be estimated; Based on the extended state vector, full flux vector, driving input, static covariates, and parameter set, a forward prediction interface is established to advance the state and synchronously output the full flux at time steps.

4. The method according to claim 1, characterized in that, The process of establishing an observation operator library and a unified error statistics framework includes: Define corresponding observation mapping rules for the classification of multi-source observation data, and construct an observation operator library; For each type of observation data, an observation error model is established, and the observation error covariance matrix is ​​constructed. Establish a process error or background error model to describe the uncertainty of model structure, driving forces, and parameters; By introducing a bias state or reducing the weight of representative errors, the inconsistency between system bias and scale support is explicitly addressed.

5. The method according to claim 1, characterized in that, The process of performing information content diagnosis and generating assimilation configuration includes: The contribution of different observation sources to the constraints of state, flux and parameters is evaluated, and observability analysis is performed to obtain the observability analysis results. Based on the observability analysis results, assimilation variables were screened and stratified. Develop assimilation and update strategies in groups or phases; Regularization and flux decomposition rules are introduced for weakly identifiable flux components; The weights of soft and hard constraints are configured collaboratively and interfaced with the constraint set and error statistics framework. Perform a consistency check before assimilation and output an assimilation configuration file containing a list of variables, update order, regularization rules, and weight configuration.

6. The method according to claim 1, characterized in that, The process of constructing and enforcing a set of constraints includes: Construct equation constraints for the conservation of nitrogen and phosphorus and the total flux closure; Construct inequality constraints to ensure that state variables, fluxes, and key derived variables are non-negative and within the physical boundary; Construct a parameter feasible region and scale constraints to ensure that process parameters are physically feasible and spatially reasonable; Construct process coupling constraints to ensure the consistency of substrate constraints, carrier constraints, and morphological transformation in order to guarantee the mechanistic self-consistency of the process; Define a constraint implementation operator for performing feasible region correction on candidate solutions and outputting a closure verification diagnostic quantity.

7. The method according to claim 1, characterized in that, The process of performing constraint-based joint data assimilation and solving includes: Generate an initialization set of joint assimilation variables containing states and parameters based on the prior distribution; The forward prediction interface is used to advance the set members to generate prior states and prior total flux. Based on the aforementioned observation operator library and unified error statistics framework, multi-source observations are fused to update the prior set and obtain candidate posterior solutions. The constraint set and its constraint implementation operators are invoked to correct the candidate posterior solution, thereby obtaining a feasible posterior solution that satisfies conservation closure and physical feasibility. The feasible posterior solutions are subjected to set statistics to quantify the posterior uncertainty and perform attribution decomposition. Perform consistency diagnosis and quality classification on the posterior results, and output a unified posterior benchmark data package with admission criteria.

8. The method according to claim 1, characterized in that, The process of performing scenario differential potential assessment includes: Using the unified posterior baseline data package as a common starting point, baseline scenarios and action scenarios are constructed, wherein the action scenarios are implemented by modifying the driving input; Under the same model, constraint set, and initial state, forward extrapolation is performed on the baseline scenario and the action scenario respectively, and their respective state trajectories and full flux ledgers are output; By pairing and differencing members of the same posterior set, the flux difference and state difference between the action scenario and the baseline scenario are calculated, and a potential accounting index system is constructed. The posterior uncertainty is propagated in the scenario difference, the difference results are aggregated and statistically analyzed, the potential confidence interval and robustness probability are output, and the risk screening and feasibility assessment are combined with engineering constraints. The differential potential results of the accounting units are spatially aggregated according to the specified partition boundaries, and the statistical and priority calculation results at the partition scale are output.

9. The method according to claim 8, characterized in that, The construction of the proposed measures includes: by adjusting the amount of fertilizer applied, the timing of fertilizer application, the method of fertilizer application, the form of fertilizer, the irrigation system, and the management of organic materials in the driving input, to obtain nitrogen control measures, phosphorus control measures, water and transport control measures, and organic fertilizer and straw management measures with different intensity levels.

10. The method according to claim 1, characterized in that, The outputs of the difference potential results that include uncertainty include: Output the posterior state and the full-process accounting ledger of total flux; Output the scenario potential list and spatial distribution product of the differential potential results; Output probabilistic expressions and risk characterization indicators for posterior uncertainty and difference potential; Output the verification indicators and log records for conservation closure and consistency; Based on feasibility assessment and robustness probability, the results of optimal measures and comprehensive evaluation are output.

Images