Method and apparatus for simulating agricultural biogeochemical processes coupling a microenable model with machine learning
Patent Information
- Application Number
- CN202611076674.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-07-20
- Publication Date
- 2026-09-22
AI Technical Summary
这些硬阈值(Hard Thresholds)和开关逻辑(Switch Logic)导致模型输出对输入变量的梯度在临界点不连续或为零,从而阻断机器学习算法的梯度优化路径,而无法利用农田观测数据联合优化机器学习算法,从而无法自动校准用于驱动农田生物地球化学模型的物理参数,进而影响农田生物地球化学模型在农田区域尺度的模拟精度
[0019]根据本公开的各方面,通过对生物地球化学模型中反应过程中的离散逻辑进行连续化处理所得到的可微化模型可以与参数生成网络形成端到端的耦合优化,且能够使目标损失梯度穿过可微化模型至参数生成网络,指导参数生成网络和可微化模型的参数更新,有效得到适配目标农田的优化后的参数生成网络和优化后的可微化模型,这样利用优化后的参数生成网络可以实现自动校准用于控制可微化模型的物理参数,且优化后的可微化模型中部分模型参数更加精确,相当于实现针对可微化模型中各类参数的参数校准,从而可以利用优化后的可微化模型以及优化后的参数生成网络针对目标农田进行更精确地在线演化模拟,提高模拟精度。
Smart Images

Figure CN122797152A_ABST
Abstract
Description
Technical Field
[0001] This disclosure relates to the intersection of smart agriculture and scientific computing, and in particular to a method and apparatus for simulating farmland biogeochemical processes by coupling a differential model with machine learning. Background Technology
[0002] Existing biogeochemical models, such as the DNDC (Denitrification-Decomposition) model, the CENTURY model, and the DAYCENT model, are core tools for understanding the exchange of matter and energy in agricultural ecosystems and assessing greenhouse gas emissions and crop yields. These models are based on thermodynamic laws and reaction kinetic equations (such as the Arrhenius equation and the Michaelis-Menten equation) to quantitatively simulate carbon and nitrogen transformation processes in the soil-atmosphere-vegetation continuum.
[0003] With the development of the AI for Science paradigm, researchers are committed to using massive amounts of multi-source data such as satellite remote sensing and flux towers, combined with the powerful feature extraction capabilities of machine learning, to drive the above-mentioned biogeochemical cycle model in order to achieve high-precision simulation at the farmland scale.
[0004] However, existing biogeochemical models generally employ discrete empirical rules to describe complex biochemical reactions. For example, the control of soil moisture on nitrification / denitrification is often triggered by specific water-filled porosity (WFPS) thresholds (e.g., denitrification is initiated when WFPS > 60%), and the effect of temperature on enzyme activity is often expressed as a piecewise function. These hard thresholds and switch logic cause the gradient of the model output with respect to the input variables to be discontinuous or zero at critical points, thus blocking the gradient optimization path of the machine learning algorithm. Consequently, it is impossible to jointly optimize the machine learning algorithm using farmland observation data, and therefore, it is impossible to automatically calibrate the physical parameters used to drive the farmland biogeochemical model, thereby affecting the simulation accuracy of the farmland biogeochemical model at the farmland regional scale. Summary of the Invention
[0005] In view of this, this disclosure proposes a method and apparatus for simulating farmland biogeochemical processes by coupling a differentiable model with machine learning, which can realize end-to-end optimization of the coupling between a parameter generation network based on machine learning and a differentiable model, thereby improving the simulation accuracy of the differentiable model.
[0006] According to one aspect of this disclosure, a method for simulating farmland biogeochemical processes is provided, comprising: acquiring spatial data, meteorological data, agricultural management data, and crop data of a target farmland, wherein the spatial data includes at least soil attribute data, topographic data, and land use data of the target farmland; generating physical parameters using a parameter generation network based on machine learning based on the spatial data, the meteorological data, and the crop data; the physical parameters being used to control reaction processes in a differential model, wherein the differential model is a model obtained by continuously processing the discrete logic in the reaction process of a general biogeochemical model, and the biogeochemical model being used to simulate carbon and nitrogen transformation, greenhouse gas emissions, and crop growth in an agricultural ecosystem; and using the differential model to generate physical parameters based on the meteorological data, the soil attribute data, and the crop data. The agricultural management data, crop data, and physical parameters are used to perform time-step evolution simulations to obtain simulation results. The simulation results include at least one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions of the target farmland. Based on the difference between the simulation results and the observation results for the target farmland, a target loss is determined. The observation results include one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions obtained from actual observations of the target farmland. The model parameters of the differentiable model and the network parameters of the parameter generation network are optimized using the target loss, so that online evolution simulations can be performed using the optimized differentiable model based on the physical parameters regenerated by the optimized parameter generation network.
[0007] In one possible implementation, the reaction process includes at least: soil organic matter decomposition, nitrification and denitrification, crop growth, and soil hydrothermal processes; the discrete logic includes: environmental response logic and inventory limit constraint logic; the environmental response logic represents the threshold control logic of environmental control variables on the process rate of the reaction process, wherein the threshold control logic is manifested as a sudden change or truncation of the process rate when the environmental control variable reaches the threshold; the process rate includes at least: nitrification rate, denitrification rate, and organic matter decomposition rate; the environmental control variables include at least: soil water porosity, temperature, and soil pH or pH state; the inventory limit constraint logic includes logic that limits the actual reaction flux during the reaction process to not exceed the corresponding inventory limit; the actual reaction flux includes at least: nitrification flux, denitrification flux, and organic matter decomposition flux; wherein, for the reaction process in a general biogeochemical model, the discrete logic includes: environmental response logic and inventory limit constraint logic; the discrete logic includes: environmental response logic and inventory limit constraint logic; the discrete logic includes: environmental response logic and environmental response ... The continuous processing of discrete logic includes: for the environmental response logic, replacing it with smoothed continuous logic constructed using a smoothing function, wherein the smoothed continuous logic is used to make the process rate continuously change around a threshold; for the inventory upper limit constraint logic, replacing it with flux scaling logic with non-negative constraints or smoothed minimum logic; wherein the flux scaling logic with non-negative constraints is used to scale the theoretical reaction flux using a scaling factor constructed based on the inventory upper limit, so that the scaled actual reaction flux is gradually limited by the inventory upper limit and does not exceed the inventory upper limit; the smoothed minimum logic is used to make the actual reaction flux approximately the minimum of the inventory upper limit and the theoretical reaction flux, wherein the theoretical reaction flux is the theoretical value generated by the reaction process in the differentiable model, and the actual reaction flux is the actual value after the inventory upper limit constraint logic imposes an upper limit constraint on the theoretical reaction flux.
[0008] In one possible implementation, the differentiable model is still obtained by reconstructing the dispersed state variables in the biogeochemical model into a specified array format; wherein, reconstructing the dispersed state variables in the biogeochemical model into a specified array format includes: reconstructing the state variables used to characterize soil attribute data into a two-dimensional array of soil layer dimension × grid dimension; reconstructing the state variables used to characterize water state into a three-dimensional array of water state type dimension × soil layer dimension × grid dimension; reconstructing the state variables used to characterize carbon pool state and nitrogen pool state into a two-dimensional array of soil layer dimension × grid dimension; reconstructing the state variables used to characterize crop state into a two-dimensional array of crop type dimension × grid dimension; and reconstructing the state variables used to characterize agricultural management measures into a two-dimensional array of date dimension × grid dimension; wherein, the soil layer dimension represents different soil layers divided downwards from the surface in the vertical dimension, the grid dimension represents different geographical locations, plots, subdivided sampling points or raster cells in the horizontal spatial dimension; the water state type dimension represents different types of water state; and the crop type dimension represents different types of crops.
[0009] In one possible implementation, during the time-step evolutionary simulation of the differentiable model, the differentiable model undergoes parallel time-step evolution based on a state array, parameter array, driving array, and process mask, along with a grid dimension. The state array, parameter array, driving array, and process mask maintain a uniform grid dimension. The state array characterizes the state data generated by the differentiable model during the evolutionary simulation, which changes with time steps. This state data includes water state, carbon pool state, nitrogen pool state, and crop state. The parameter array indicates the parameter data characterizing various reaction processes in the differentiable model. This parameter data includes soil property data, crop data, and physical parameters output by the parameter generation network. The driving array characterizes meteorological data and agricultural management data. The process mask controls the start and stop of any reaction process under any grid, any date, any soil layer, any crop type, and any water state type.
[0010] In one possible implementation, the soil property data includes at least: soil volume, soil clay content, field water holding capacity, wilting point, soil pH, and soil structure; the meteorological data includes at least: temperature, precipitation, humidity, wind speed, and radiation varying over time; the agricultural management data includes management measures taken for the target farmland from sowing, fertilization, irrigation to harvest; the crop data includes: the type of crop to be sown in the target farmland and crop growth information; wherein, the step of using the differentiable model to perform time-step evolutionary simulation based on the meteorological data, soil property data, agricultural management data, crop data, and physical parameters to obtain simulation results includes: using a general-purpose parallel computing platform to perform time-step parallel computation on various reaction processes in the differentiable model in the grid dimension according to the state array, the parameter array, the driving array, and the process mask until the evolution termination condition is reached to obtain simulation results, wherein the state array in each time step is updated in parallel by the image processor thread in the general-purpose parallel computing platform.
[0011] In one possible implementation, optimizing the model parameters of the differentiable model and the network parameters of the parameter generation network using the target loss includes: backpropagating the loss gradient of the target loss relative to the simulation results through the differentiable model to the parameter generation network to update the network parameters of the parameter generation network and the model parameters of the differentiable model; generating new physical parameters using the updated parameter generation network, and regenerating the simulation results using the updated differentiable model based on the new physical parameters, until the difference between the simulation results regenerated by the updated differentiable model and the observation results is minimized, thereby obtaining the optimized parameter generation network and the optimized differentiable model.
[0012] In one possible implementation, the method further includes: during the time-step evolutionary simulation of the differentiable model, saving checkpoints at preset time intervals to obtain a checkpoint sequence, and using the checkpoint sequence to perform gradient backpropagation; wherein each checkpoint in the checkpoint sequence includes at least: a time point, the corresponding soil hydrothermal state, carbon pool state, nitrogen pool state, crop state, physical parameters, driving data location, and process mask; the driving data location represents the read position of the driving array; wherein, backpropagating the target loss relative to the simulation result as a gradient to the parameter generation network through the differentiable model includes: extracting gradients from the checkpoint sequence and the evolutionary model... During the simulation process, starting from the nearest checkpoint at the last time step, the following calculations are performed sequentially until all checkpoints in the checkpoint sequence are traversed: the differentiable model is restored to the current checkpoint, and the evolution simulation within a preset time interval is re-executed based on the current checkpoint to obtain a forward computation path, which represents the complete computation process in the evolution simulation within the preset time interval; gradient calculation is performed based on the forward computation path to obtain the local gradient; the target gradient is determined based on the local gradient and the loss gradient, and the target gradient is passed to the previous checkpoint, the differentiable model, and the parameter generation network; the previous checkpoint is the checkpoint saved before the current checkpoint in the checkpoint sequence.
[0013] In one possible implementation, the step of calculating gradients based on the forward computation path to obtain local gradients includes: for the locally defined response functions in the differentiable model, calculating gradients using the derivative formulas of the pre-derived locally defined response functions based on the forward computation path to obtain a first type of local gradient; wherein the locally defined response functions include at least: soil water-filling porosity response function, temperature response function, pH response function, smoothing continuous logic, and smoothing minimum logic; the derivative formula is implemented as a custom gradient kernel on a general-purpose parallel computing platform; for the complex reaction processes in the differentiable model, establishing local computation graphs corresponding to these complex reaction processes based on the forward computation path, and calculating gradients based on the local computation graphs to obtain a second type of local gradient; wherein the complex reaction processes include at least: organic matter decomposition process, nitrogen transformation process, and crop growth process; wherein the step of determining the target gradient based on the local gradients and the loss gradient includes: combining and accumulating the first type of local gradient, the second type of local gradient, and the loss gradient based on the chain rule to obtain the target gradient.
[0014] In one possible implementation, the physical parameters include at least: organic matter decomposition rate parameters, crop parameters, nitrification rate parameters, denitrification rate parameters, water response parameters, temperature response parameters, acid-base response parameters, and soil water-filled porosity response parameters; the crop parameters include at least crop physiological parameters, phenological parameters, and biomass allocation parameters; wherein, after generating the physical parameters, the method further includes: mapping the physical parameters to corresponding physical constraint boundaries to obtain mapped physical parameters, so that the differentiable model can perform time-step evolution simulation based on the mapped physical parameters; wherein, mapping the physical parameters to corresponding physical constraint boundaries includes: mapping the water response parameters, temperature response parameters, acid-base response parameters, soil water-filled porosity response parameters, and crop parameters to their respective corresponding preset reasonable intervals; and mapping the organic matter decomposition rate parameters, nitrification rate parameters, and denitrification rate parameters to non-negative values.
[0015] According to another aspect of this disclosure, a farmland biogeochemical process simulation device is provided, comprising: a data acquisition module for acquiring spatial data, meteorological data, agricultural management data, and crop data of a target farmland, wherein the spatial data includes at least soil attribute data, topographic data, and land use data of the target farmland; a parameter generation module for generating physical parameters based on the spatial data, the meteorological data, and the crop data using a machine learning-based parameter generation network; wherein the physical parameters are used to control the reaction process in a differential model, the differential model being a model obtained by continuously processing the discrete logic in the reaction process of a general biogeochemical model, the biogeochemical model being used to simulate carbon and nitrogen transformation, greenhouse gas emissions, and crop growth in an agricultural ecosystem; and an evolution simulation module for using the differential model to acquire spatial data, meteorological data, agricultural management data, and crop data of a target farmland, wherein the spatial data includes at least soil attribute data, topographic data, and land use data of the target farmland; a parameter generation module for generating physical parameters based on the spatial data, the meteorological data, and the crop data using the machine learning-based parameter generation network; wherein the physical parameters are used to control the reaction process in a differential model, wherein the differential model is obtained by continuously processing the discrete logic in the reaction process of a general biogeochemical model, and the biogeochemical model is used to simulate carbon and nitrogen transformation, greenhouse gas emissions, and crop growth in an agricultural ecosystem; and an evolution simulation module for using the differential model to acquire physical parameters based on the meteorological data, the soil attribute data, and the crop data. The system performs time-step evolutionary simulations using the target farmland's data, agricultural management data, crop data, and physical parameters to obtain simulation results. These simulation results include at least one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions. A loss determination module is used to determine a target loss based on the difference between the simulation results and observational results for the target farmland. The observational results include one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions obtained from actual observations of the target farmland. A model optimization module is used to optimize the model parameters of the differentiable model and the network parameters of the parameter generation network using the target loss, so that the optimized differentiable model can be used to perform online evolutionary simulations based on the physical parameters regenerated by the optimized parameter generation network.
[0016] According to another aspect of this disclosure, an electronic device is provided, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above-described method.
[0017] According to another aspect of this disclosure, a non-volatile computer-readable storage medium is provided, on which a computer program is stored, which, when executed by a processor, implements the steps of the above-described method.
[0018] According to another aspect of this disclosure, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps of the above-described method.
[0019] According to various aspects of this disclosure, the differentiable model obtained by continuously processing the discrete logic in the reaction process of the biogeochemical model can form an end-to-end coupled optimization with the parameter generation network. Furthermore, the target loss gradient can pass through the differentiable model to the parameter generation network, guiding the parameter updates of both the parameter generation network and the differentiable model. This effectively yields an optimized parameter generation network and an optimized differentiable model adapted to the target farmland. The optimized parameter generation network can then be used to automatically calibrate the physical parameters controlling the differentiable model. Moreover, some model parameters in the optimized differentiable model are more accurate, effectively achieving parameter calibration for various parameters within the differentiable model. Therefore, the optimized differentiable model and the optimized parameter generation network can be used to perform more accurate online evolution simulations of the target farmland, improving simulation accuracy.
[0020] Other features and aspects of this disclosure will become clear from the following detailed description of exemplary embodiments with reference to the accompanying drawings. Attached Figure Description
[0021] The accompanying drawings, which are included in and form part of this specification, illustrate exemplary embodiments, features, and aspects of this disclosure together with the specification and serve to explain the principles of this disclosure.
[0022] Figure 1 A flowchart is shown for a method for simulating farmland biogeochemical processes according to an embodiment of the present disclosure.
[0023] Figure 2 The diagram illustrates a discrete logic according to an embodiment of the present disclosure and a schematic diagram of the input-output curves obtained by performing continuous processing on the discrete logic.
[0024] Figure 3 A schematic diagram of the calculation flow of a nitration process according to an embodiment of the present disclosure is shown.
[0025] Figure 4The diagram shows an evolution simulation process of forward computation and hybrid backpropagation at the end of a miniaturized DNDC device according to an embodiment of the present disclosure.
[0026] Figure 5 A schematic diagram of the architecture of an evolution simulation system that couples a differentiable model and a parameter generation network according to an embodiment of the present disclosure is shown.
[0027] Figure 6 A block diagram of a farmland biogeochemical process simulation apparatus according to an embodiment of the present disclosure is shown.
[0028] Figure 7 A block diagram of an electronic device according to an embodiment of the present disclosure is shown. Detailed Implementation
[0029] Various exemplary embodiments, features, and aspects of this disclosure will now be described in detail with reference to the accompanying drawings. The same reference numerals in the drawings denote elements that have the same or similar functions. Although various aspects of the embodiments are shown in the drawings, they are not necessarily drawn to scale unless specifically indicated otherwise.
[0030] As used herein, the terms “comprising,” “including,” “having,” or variations thereof are open-ended and include one or more of the stated features, integrals, elements, steps, components, or functions, but do not exclude the presence or addition of one or more other features, integrals, elements, steps, components, functions, or groups thereof.
[0031] When an element is referred to as “connected,” “coupled,” “responding,” or a variation thereof relative to another element, it may be directly connected, coupled, or responding to another element, or there may be an intermediate element present.
[0032] Although the terms first, second, third, etc., may be used herein to describe various elements / operations, these elements / operations should not be limited by these terms. These terms are only used to distinguish one element / operation from another. Therefore, without departing from the teachings of the inventive concept, a first element / operation in some embodiments may be referred to as a second element / operation in other embodiments.
[0033] The term “exemplary” as used herein means “serving as an example, embodiment, or illustration.” Any embodiment illustrated herein as “exemplary” is not necessarily to be construed as superior to or better than other embodiments.
[0034] Furthermore, to better illustrate this disclosure, numerous specific details are set forth in the following detailed description. Those skilled in the art will understand that this disclosure can be practiced without certain specific details. In some instances, methods, means, components, and circuits well known to those skilled in the art have not been described in detail in order to highlight the main points of this disclosure.
[0035] It should be noted that the information (including but not limited to user device information, user personal information, etc.), data (including but not limited to data used for analysis, data stored, data displayed, etc.) and signals involved in this application are all authorized by the user or fully authorized by all parties, and the collection, use and processing of related data must comply with the relevant laws, regulations and standards of the relevant regions.
[0036] The simulation method of this disclosure can be deployed on various terminal devices through software or hardware modifications. The terminal devices involved in this disclosure can refer to devices with wireless and / or wired connection functions. Wireless connection means that they can connect to other devices via wireless connection methods such as Wi-Fi and Bluetooth. The terminal devices involved in this disclosure can also communicate with other devices via wired connection functions. The terminal devices involved in this disclosure can be touchscreen, non-touchscreen, or screenless. Touchscreen devices can be controlled by clicking or swiping on the display screen using fingers or styluses. Non-touchscreen devices can connect to input devices such as mice, keyboards, and touch panels to control the terminal device. Screenless devices can be, for example, screenless Bluetooth speakers. For example, the terminal devices in this application can include, but are not limited to, user equipment (UE), mobile devices, user terminals, terminals, handheld devices, tablet computers, laptops, PDAs, computing devices, etc.
[0037] The simulation method of this disclosure can also be deployed on a server, which can be located in the cloud or locally, and can be a physical device or a virtual device, such as a virtual machine or container, with wireless communication capabilities. These wireless communication capabilities can be configured in the server's chip (system) or other components. This can refer to a device with wireless connectivity, meaning it can connect to other servers or terminal devices via wireless connections such as Wi-Fi or Bluetooth. The server involved in this disclosure can also have wired communication capabilities. For example, the server in this disclosure can be located in the cloud, receiving spatial data, meteorological data, agricultural management data, and crop data of the target farmland sent by the terminal device. It then uses the simulation method deployed on the server to optimize the parameter generation network or physical parameters, performs online evolution simulation based on the regenerated physical parameters using a differentiable model, or directly performs online evolution simulation using the optimized physical parameters. The online simulation results are then returned to the terminal device for display to the user.
[0038] Figure 1 A flowchart illustrating a method for simulating farmland biogeochemical processes according to an embodiment of this disclosure is shown. Figure 1 As shown, the method includes steps S11 to S15.
[0039] In step S11, spatial data, meteorological data, agricultural management data and crop data of the target farmland are acquired. The spatial data includes at least soil attribute data, topographic data and land use data of the target farmland.
[0040] It should be understood that those skilled in the art can use any data acquisition method and data acquisition channel known in the art to obtain spatial data, meteorological data, agricultural management data and crop data for the target farmland. The target farmland can be any farmland to be simulated for biogeochemical processes, and this disclosure does not limit this.
[0041] The soil property data may include at least: soil volume, soil clay content, field water holding capacity, wilting point, soil pH, and soil structure. Topographic data may include at least: elevation and slope. Land use data describes the spatial layout of the target farmland, including the location of major crop areas, buffer zones, ditches, and field ridges. Optionally, the spatial data may also include remote sensing data of the target farmland, which may be data collected by remote sensing satellites; this embodiment of the disclosure does not limit this.
[0042] The meteorological data includes at least the following: temperature, precipitation, humidity, wind speed, and radiation (i.e., solar radiation) that vary over time; agricultural management data includes the management measures taken for the target farmland from sowing, fertilization, irrigation to harvest, such as: sowing date, fertilization date and type and amount of fertilizer, irrigation date and amount, amount of crop residue returned to the field, and harvest date; crop data includes: the type of crop to be sown in the target farmland and crop growth information. The crop growth information is used to describe the crop growth characteristics of different crops during the growth cycle from sowing to harvest (such as germination and seedling stage, vegetative growth stage, reproductive growth stage, and maturity stage).
[0043] In step S12, a parameter generation network based on machine learning is used to generate physical parameters based on spatial data, meteorological data, and crop data. The physical parameters are used to control the reaction process in the minimizable model. The minimizable model is a model obtained by processing the discrete logic in the reaction process of the general biogeochemical model into a continuous model. The biogeochemical model is used to simulate carbon and nitrogen conversion, greenhouse gas emissions, and crop growth in agricultural ecosystems.
[0044] In practical applications, the parameter generation network based on machine learning can be constructed using one or more of the following: Convolutional Neural Network (CNN), Artificial Neural Network (ANN), and Long Short-Term Memory (LSTM). This disclosure does not limit the type or structure of the parameter generation network. For example, the parameter generation network may include an encoding module and a decoding module. The encoding module extracts data features from the input spatial data, meteorological data, and crop data. The decoding module generates physical parameters based on the data features extracted by the encoding model. The encoding module may include, for example, a Convolutional Neural Network (CNN) and a Long Short-Term Memory (LSTM) network. The CNN can be used to process spatial data, and the LSTM network can be used to process meteorological and crop data. This disclosure does not limit the specific implementation of the LSTM network. In some embodiments, the parameter generation network may be a pre-trained network model that initially generates physical parameters based on spatial data, meteorological data, and crop data. Subsequently, the parameter generation network can be optimized using observation results from the target farmland to obtain a more accurate parameter generation network adapted to the target farmland for online evolution simulation of the target farmland.
[0045] Physical parameters can be understood as parameters required in a differential model (or biogeochemical model) to control the reaction process. These physical parameters typically have clear physical or biogeochemical meanings in the biogeochemical model and are difficult to measure directly. In some embodiments, physical parameters include at least: organic matter decomposition rate parameters, crop parameters, nitrification rate parameters and denitrification rate parameters, water response parameters, temperature response parameters, acid-base response parameters, soil water-filled porosity response parameters, etc. Among them, the organic matter decomposition rate parameter is used to control the decomposition rate during the organic matter decomposition process; crop parameters may include crop physiological parameters (such as crop water requirement, carbon-nitrogen ratio, and nitrogen fixation constant), phenological parameters (controlling the crop growth stage response to temperature, time, or management measures), and crop physiological parameters. The parameters for controlling crop growth rate include: growth rate (controlling the rate of crop growth), biomass allocation (controlling the allocation of biomass among grains, stems, and roots), and root parameters (controlling the absorption of water and nutrients by crop roots from the soil, as well as the carbon and nitrogen cycling process within the root system); nitrification rate parameters are used to control the nitrification rate during nitrification; denitrification rate parameters are used to control the denitrification rate during denitrification; water response parameters are used to control the effect of soil moisture content on the rates of nitrification, denitrification, and organic matter decomposition; temperature response parameters are used to control the effect of soil temperature on the rates of decomposition, mineralization, and nitrification; pH response parameters are used to control the effect of soil pH on nitrogen transformation; and soil water-filled porosity response parameters are used to control the effect of soil water-filled porosity (WFPS) on the rates of nitrification, denitrification, and organic matter decomposition.
[0046] The parameter generation network can obtain parameters such as nitrification rate, denitrification rate, water response, acid-base response, and soil water-filled porosity response by processing spatial data. It can also obtain parameters such as organic matter decomposition rate, water response, temperature response, and crop parameters by processing meteorological and crop data. Different physical parameters can be obtained by processing different input data. Those skilled in the art can design a parameter generation network to output corresponding types of physical parameters based on the input data type; this disclosure does not limit such design.
[0047] The general biogeochemical models used here refer to existing traditional biogeochemical models, such as the aforementioned DNDC model, CENTURY model, and DAYCENT model. Any model capable of simulating carbon and nitrogen conversion, greenhouse gas emissions, and crop growth in an agricultural ecosystem is acceptable, and this disclosure does not impose any limitations on these aspects. The reaction processes in the biogeochemical model may include, for example, at least: soil organic matter decomposition, nitrification and denitrification, crop growth, and soil hydrothermal processes. These reaction processes involve discrete logic, typically manifested as hard thresholds, hard stages, or on / off logic, making the biogeochemical model non-differentiable.
[0048] The discrete logic in the reaction process within the biogeochemical model can include: environmental response logic and inventory ceiling constraint logic. Environmental response logic represents the threshold control logic of environmental control variables on the process rate of the reaction (e.g., the control of WFPS, temperature, pH, or acid-base state on the rates of nitrification, denitrification, and organic matter decomposition). The threshold control logic manifests as a sudden change or truncation of the process rate when the environmental control variable reaches a threshold (e.g., when WFPS > 0.05). Denitrification is initiated at 60%, which is a hard threshold or hard cutoff. The process rates include at least: nitrification rate, denitrification rate, and organic matter decomposition rate. Environmental control variables include at least: soil water-filled porosity, temperature, and soil pH or pH state. The stock cap constraint logic includes logic that limits the actual reaction flux during the reaction process to not exceeding the corresponding stock cap (e.g., nitrification and denitrification fluxes cannot exceed the stock caps of available ammonium nitrogen pool NH4, nitrate nitrogen pool NO3, easily decomposable residual carbon pool RCVL, or urea nitrogen pool urea), which usually exhibits a slight hard cutoff or hard minimum. The actual reaction flux includes at least: nitrification flux, denitrification flux, and organic matter decomposition flux. Therefore, different treatments can be adopted for different types of discrete logic in the reaction process in the biogeochemical model. Specifically, the continuous processing of discrete logic in the reaction process of a general biogeochemical model can include:
[0049] For environmental response logic, smooth continuous logic constructed using a smoothing function is used to replace the environmental response logic. The smooth continuous logic is used to make the process rate change continuously around a threshold.
[0050] For the inventory ceiling constraint logic, use flux scaling logic with non-negative constraints or smooth minimum value logic to replace the inventory ceiling constraint logic.
[0051] Among them, the flux scaling logic with non-negative constraints is used to scale the theoretical reaction flux using a scaling factor constructed based on the inventory upper limit, so that the scaled actual reaction flux is gradually limited by the inventory upper limit and does not exceed the inventory upper limit; the smoothing minimum value logic is used to make the actual reaction flux approximate the minimum value between the inventory upper limit and the theoretical reaction flux, wherein the theoretical reaction flux is the theoretical value generated by the reaction process in the differentiable model, and the actual reaction flux is the actual value after the inventory upper limit constraint logic applies an upper limit constraint to the theoretical reaction flux.
[0052] In practical applications, smoothing functions may include, but are not limited to, the Sigmoid function, the Softplus function, the smoothing truncation function smooth_clip(x, min, max, sharpness), or various continuous weighting functions. This disclosure does not limit the scope of these functions, where sharpness is the sharpness parameter in the smoothing truncation function. It should be understood that those skilled in the art can construct corresponding smooth continuous logic based on environmental response logic using smoothing functions, ensuring that the process rate changes continuously around a threshold. This disclosure does not limit the scope of these logics. For example, Figure 2 The diagram illustrates discrete logic and the input-output curves obtained by performing continuous processing on discrete logic, as shown below. Figure 2 (a) shows the threshold control logic (i.e., step function) in traditional DNDC, which is characterized by a sudden change (i.e., a step) in the output (e.g., rate factor, i.e., process rate) when the input (e.g., soil water-filled porosity) reaches the threshold T. This results in a hard threshold (non-differentiable). This type of step function can be approximated using a combination of smooth_clip and Sigmoid / Softplus, thereby replacing discrete logic (e.g., if wfps > threshold) with continuous logic. For example, it can generate a continuous and non-zero gradient of nitrification / denitrification rates with respect to variables such as soil moisture and temperature when crossing the threshold. Figure 2 (b) shows that the process rate of the input and output after using smooth continuous logic changes continuously at the threshold T, that is, it has a smooth gradient (i.e., differentiable) at the threshold.
[0053] As mentioned above, the original form of inventory cap constraint logic is usually hard truncation or hard minimum. For example, in the nitrification process, the model calculates the theoretical nitrification flux F_theory, but the actual nitrification flux cannot exceed the inventory cap S_NH4 of the available ammonium nitrogen pool NH4 in the current soil layer and current grid. The inventory cap can be understood as the available inventory. The traditional inventory cap constraint logic can be written as: F_actual = min(F_theory, S_NH4), or, if F_theory > S_NH4, then F_actual = S_NH4, otherwise F_actual = F_theory. Although this type of hard minimum or "if-condition logic" can ensure that the inventory is not deducted to a negative value, there are non-differentiability or gradient discontinuity problems when F_theory is close to S_NH4.
[0054] Therefore, the first approach is to replace the inventory cap constraint logic with flux scaling logic with non-negative constraints. Specifically, a continuous scaling factor `scale` can be constructed based on the inventory cap `S_NH4`, with `scale` ranging from 0 to 1, and the actual nitrification flux being the theoretical nitrification flux multiplied by this scaling factor, i.e., `F_actual = F_theory × scale`. Thus, when the theoretical nitrification flux is much smaller than the inventory cap, `scale` approaches 1, and the actual nitrification flux approaches the theoretical nitrification flux. When the theoretical nitrification flux approaches or exceeds the inventory cap, `scale` continuously decreases, gradually limiting the actual nitrification flux to the inventory cap, while ensuring that `F_actual` does not exceed the inventory cap of available NH4. The inventory status can be updated, such as `NH4_new = NH4_old - F_actual`, and the non-negative constraint ensures that `NH4_new ≥ 0`.
[0055] The second approach is to replace the inventory ceiling constraint logic with a smoothed minimum logic. Specifically, a continuously differentiable smoothed minimum function (smooth_min) replaces the hard minimum function (min). For example, F_actual = min(F_theory, S_NH4) is replaced with F_actual = smooth_min(F_theory, S_NH4). Here, the smooth_min function makes F_actual numerically approximate the minimum of F_theory and S_NH4, while maintaining continuous change and a computable gradient even when F_theory and S_NH4 are close. This preserves the physical constraint that the actual nitrification flux cannot exceed the inventory ceiling while avoiding gradient discontinuities caused by hard truncation.
[0056] For example, taking the nitrification process as an example, the original DNDC model may contain the following discontinuous discrete logic during digestion: 1. When WFPS, temperature, soil pH, or acid-base state reach a certain condition, the nitrification rate is started, stopped, or cut off; 2. The theoretical nitrification flux cannot exceed the upper limit of the current available ammonium nitrogen pool NH4, so a min or if judgment is used for hard truncation. Therefore, continuous processing mainly occurs in the following two locations: First, in the environmental response calculation stage, the hard threshold control of environmental control variables such as WFPS, temperature, and pH on the nitrification rate is replaced with continuous smooth logic (such as a continuous response function or a smooth gating function). For example, the hard judgment that changes the nitrification rate when WFPS exceeds a certain threshold is replaced with a response factor that changes continuously with WFPS. In this way, the nitrification rate will not suddenly jump around the threshold, but will change continuously, and the gradient with respect to variables such as WFPS, temperature, and pH will be preserved. Second, in the inventory constraint stage, the nitrification process first calculates the theoretical nitrification flux F_theory based on the environmental response and the upper limit of the NH4 inventory. Then, it is constrained by the available NH4 inventory, i.e., as mentioned above, F_actual = min(F_theory, NH4_S) or the flux is truncated to within the NH4 inventory using an if condition. For continuous processing, smooth_min or continuous flux scaling can be used to smoothly limit the actual nitrification flux as it approaches the upper limit of the NH4 inventory, rather than hard truncating it. This ensures that the actual flux does not exceed the NH4 inventory and that the flux maintains a calculable gradient with respect to the theoretical rate parameters, environmental response parameters, and NH4 state.
[0057] Furthermore, after performing the aforementioned continuum transformation process on the discrete logic of the nitration process, the following can be obtained: Figure 3 The calculation flow of the nitration process is shown below, such as Figure 3As shown, the nitrification process takes ammonium nitrogen pool NH4, nitrate nitrogen pool NO3, soil water-filled porosity WFPS, temperature, clay, soil pH or acid-base state, and flooding / tillage indicators as inputs to calculate environmental control variables such as WFPS response, temperature response, and pH response. Then, a smooth continuous logic that replaces the original hard threshold response logic outputs the nitrification rate based on these environmental control variables, and the theoretical nitrification flux is calculated based on the nitrification rate. Subsequently, a smoothing inventory constraint is applied to the flux based on the upper limit of the available ammonium nitrogen pool NH4. That is, a smoothing minimum logic or a flux scaling logic with non-negative constraints (F_actual = smooth_min(F_theory, NH4_S) or F_actual = F_theory× scale, and NH4_new ≥ 0) is used to ensure that the actual nitrification flux does not exceed the substrate inventory. Then, based on the actual nitrification flux, the daily soil nitrification amount (day_soil_nitrify), the daily nitrification nitric oxide NO flux (day_nitrify_NO), and the daily nitrification nitrous oxide N2O flux (day_nitrify_N2O) are calculated and output, and the ammonium nitrogen pool NH4 status, NO flux status, N2O flux status, and nitrate nitrogen pool NO3 status are updated simultaneously.
[0058] It should be understood that the above-described implementation method of continuous processing of discrete logic in the reaction process, taking the nitrification process as an example, is a feasible embodiment provided by this disclosure. In fact, for discrete logic in other reaction processes in the biogeochemical model, continuous processing can be carried out by referring to the implementation method for the nitrification process. For example, discrete logics such as the decomposition of easily decomposable residual carbon pool RCVL, the turnover of microbial pool CRB1 / CRB2, the mineralization or fixation of ammonia nitrogen pool NH4, and inventory scaling involved in the organic matter decomposition process can all be classified and continuously processed. This disclosure does not limit this.
[0059] It should be noted that the various smoothing functions mentioned above are not the core innovation of this disclosure. The key to this disclosure is to embed them into the specific carbon pool, nitrogen pool, environmental response, and material conservation relationship of the biogeochemical model in order to achieve continuous processing of discrete logic.
[0060] In practical applications, biogeochemical models (such as the DNDC model) also contain a third type of discrete logic for managing phenological events, such as flooding, irrigation, fertilization, harvesting, and crop development stage switching logic. These are triggered by integer or discrete events such as date, crop type, sowing date, harvest date, fertilization date, and irrigation date. For this type of discrete logic, its discrete scheduling method is retained. The variables in this type of logic have clear discrete meanings, such as the day of fertilization, the day of harvest, and which crop to plant. Directly converting date numbers or crop type numbers to continuous values could easily disrupt the original DNDC model. While managing event scheduling logic, it can perform candidate gradient calculations on parameters that participate in continuous process calculations after an event occurs, such as event quantity, response parameters, or physical parameters. In other words, the continuous physical quantities corresponding to the event can be used as candidate optimization objects. For example, the fertilization date itself remains discrete, but the fertilization amount, fertilizer conversion rate parameters, and nitrogen utilization-related parameters can be used as continuous variables; the irrigation date remains discrete, but the irrigation amount and water response parameters can be used as continuous variables; the crop type remains discrete, but the crop growth parameters, phenological response parameters, and biomass allocation parameters can be used as continuous variables.
[0061] Considering that the physical parameters generated by the parameter generation network may not satisfy the physical constraints of the biogeochemical model for these parameters, a mapping function can be used to limit them to a reasonable range before they are written into the parameter array or state array of the differentiable model (i.e., before they are input into the differentiable model for evolutionary simulation). For example, bounded mapping is used for parameters with physical upper and lower bounds; non-negative mapping is used for parameters that require non-negativity, so that the output of the parameter generation network satisfies the corresponding physical constraints. Parameters with physical upper and lower bounds refer to parameters that not only require being greater than or equal to zero, but also have a clear upper limit or a reasonable range of values. Parameters that require non-negativity refer to parameters that cannot be less than zero in a physical sense, but do not necessarily have a fixed upper limit. Therefore, after generating the physical parameters, the method further includes:
[0062] The physical parameters are mapped to their corresponding physical constraint boundaries to obtain the mapped physical parameters, so that the differentiable model can perform time-step evolution simulation based on the mapped physical parameters. The mapping of the physical parameters to their corresponding physical constraint boundaries includes mapping the water response parameters, temperature response parameters, acid-base response parameters, soil water-filled porosity response parameters, and crop parameters to their respective preset reasonable ranges; and mapping the organic matter decomposition rate parameters, nitrification rate parameters, and denitrification rate parameters to non-negative values.
[0063] It should be understood that the preset reasonable ranges corresponding to the above-mentioned different parameters can be customized based on practical experience and different parameter types, and this disclosure does not impose any restrictions on this. Bounded mapping and non-negative mapping can ensure that the physical parameters input to the differentiable model satisfy the physical meaning and model computational constraints.
[0064] As is known, existing farmland biogeochemical models typically run as independent compilers or standalone software. The model's input parameters, intermediate states, and simulation results are mainly exchanged through configuration files and result files, lacking fine-grained state access and dynamic control interfaces for external programs. Therefore, external machine learning models struggle to dynamically inject parameters, read intermediate states, adjust computational processes, or receive gradient information during model operation. This makes it difficult to construct an end-to-end parameter inversion and data assimilation system that allows the process model and machine learning model to run together. Therefore, this embodiment of the present disclosure also reconstructs the state variables stored in the original biogeochemical model code, categorized by plot, soil layer, crop type, and daily events, into a device-side array (i.e., an array format specified by the device, which is the device executing the evolutionary simulation calculation of the differentiable model). In other words, a differentiable model can also be obtained by reconstructing the state variables stored in the biogeochemical model into a specified array format.
[0065] This involves reconstructing the dispersed state variables in the biogeochemical model into a specified array format. This includes reconstructing the state variables characterizing soil properties into a two-dimensional array of soil layer dimensions × grid dimensions. Specifically, this means reconstructing state variables characterizing soil volume, soil clay content, field capacity, wilting point, soil pH, and soil structure into a two-dimensional array of soil layer dimensions × grid dimensions. For example, soil clay content can be reconstructed as clay[layer, ... [grid] represents the clay content of the th layer in the th spatial grid; where the soil layer dimension represents the different soil layers divided from the surface downwards in the vertical dimension, and the grid dimension represents the different geographical locations, plots, sample points, or grid cells in the horizontal spatial dimension; and, the state variables used to characterize the water state are reconstructed into a three-dimensional array of water state type dimension × soil layer dimension × grid dimension. The water state type dimension represents different types of water states, or the numbering dimension used to distinguish different water-related state quantities or process quantities. Soil moisture state refers to the dynamic changes in water content within the soil and its existing forms. Soil moisture is classified according to stress conditions into hygroscopic water, film water, capillary water, and gravity water, etc. For example, water[type, layer, [grid] represents the soil moisture content of the type in the layer of the grid spatial grid; and the state variables used to characterize the carbon pool and nitrogen pool are reconstructed into a two-dimensional array of soil layer dimension × grid dimension. That is, the state variables such as easily decomposable residual carbon pool, microbial carbon pool, soil organic carbon pool, nitrate nitrogen pool, ammonium nitrogen pool, urea nitrogen pool, nitric oxide and nitrous oxide can be reconstructed into a two-dimensional array of soil layer dimension × grid dimension. For example, the state of the easily decomposable residual carbon pool can be reconstructed as rcvl[layer, grid] represents the state of the easily decomposable residual carbon pool in the layer of the grid spatial grid.
[0066] Furthermore, the state variables used to characterize crop states are reconstructed into a two-dimensional array of crop type dimension × grid dimension. Crop states can include crop biomass and crop physiological states, i.e., crop biomass such as leaf biomass, stem biomass, root biomass, and grain biomass, and crop physiological states such as crop growth index or crop physiological growth index, crop development stage or crop phenological stage, etc., are reconstructed into a two-dimensional array of crop type dimension × grid dimension. Different crop types, rotation crops, multiple cropping crops, or multiple crop instances may need to be sown in the same grid. Therefore, corresponding crop numbers can be assigned to different crops to characterize crop types. The crop type dimension represents different types of crops. For example, Leaf_Wt[crop, grid] represents the leaf biomass corresponding to the crop number in the 'grid' spatial grid. Also, the state variables used to characterize agricultural management measures are reconstructed into a two-dimensional array of date dimension × grid dimension. The date dimension can specifically correspond to a date number, day step number, or day of the year. Therefore, the date dimension × grid dimension represents the management measures for different spatial grids on different simulated dates. For example, Matrix_Planting[day, [grid] represents the sowing measures for the grid-th spatial grid on the day-th simulation day; Matrix_Harvest[day, grid], Matrix_FertDay[day, grid], and Matrix_IrrDay[day, grid] represent the harvesting measures, fertilization measures, and irrigation measures for the grid-th spatial grid on the day-th simulation day, respectively.
[0067] In this context, reconstructing the dispersed state variables in the biogeochemical model into a specified array format can be understood as reconstructing the scattered scalar variables in the original C++ code into cupy.ndarray vectors. This allows for parallel processing of large-scale grids using CuPy's broadcasting mechanism, eliminating the need for serial loops. Furthermore, it enables the construction of native Python interfaces, directly exposing all physical states in memory. This transforms the model into a "library" that can be arbitrarily called and rolled back by scripts, providing a direct memory access (DMA) channel for machine learning. In other words, machine learning tensors and device-side arrays can share a view using device memory sharing or zero-copy tensor conversion.
[0068] Back Figure 1 In step S13, a differentiable model is used to perform time-step evolution simulation based on meteorological data, soil property data, agricultural management data, crop data, and physical parameters to obtain simulation results. The simulation results include at least one or more of the following: crop yield of the target farmland, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emission results.
[0069] As described above, a differential model is a model obtained by continuously processing the discrete logic in the reaction process of a general biogeochemical model. Therefore, a differential model also has the ability of a biogeochemical model to simulate carbon and nitrogen transformation, greenhouse gas emissions, and crop growth in an agricultural ecosystem. Thus, the differential model receives physical parameters output by the parameter generation network and evolves step by step according to meteorological data, soil property data, crop data, and agricultural management data, and outputs one or more simulation results among crop yield, soil hydrothermal state, carbon and nitrogen pool state, and greenhouse gas emission results. It can also output various gas or nitrogen transformation-related results, which are not limited in this embodiment.
[0070] Understandably, a differentiable model, during runtime, not only receives the physical parameters output by the parameter generation network but also the fundamental input data necessary for the model's normal time-step progression. Meteorological data, soil properties, crop data, and agricultural management data constitute the fundamental inputs or driving conditions for the model's forward evolution simulation. For example, meteorological data such as temperature, precipitation, and radiation drive daily hydrothermal and biogeochemical processes; soil properties determine soil moisture, temperature, and reaction environment; crop data is used for crop growth and carbon and nitrogen uptake calculations; and management data such as sowing, fertilization, irrigation, and harvesting trigger corresponding management events. Therefore, the raw data provides the environmental and management conditions required for the differentiable model to run, the parameter generation network provides the physical parameters to be inverted, and the differentiable model uses both as input, evolving step-by-step and outputting results such as crop yield, soil hydrothermal state, carbon and nitrogen pool status, and greenhouse gas emissions.
[0071] As mentioned above, a differentiable model can also be obtained by reconstructing the dispersed state variables in the biogeochemical model into a specified array form. Thus, the physical parameters output by the parameter generation network can be mapped to the device-side array dimension inside the differentiable model (that is, the physical parameters output by the parameter generation network can be mapped into the corresponding array form) and then input into the differentiable model for calculation. For grid-scale or global-scale parameters, they can be extended to the corresponding dimension through broadcasting. For parameters with physical constraints, bounded mapping or non-negative mapping is performed before writing them into the differentiable model.
[0072] Furthermore, in some embodiments, during time-step evolutionary simulations of the differentiable model, the differentiable model undergoes parallel time-step evolution in the grid dimension based on a state array, parameter array, driving array, and process mask; wherein the state array, parameter array, driving array, and process mask maintain a uniform grid dimension; wherein the state array is used to characterize the state data generated by the differentiable model in the evolutionary simulation that changes with time step progression, and the state data includes water state, carbon pool state, nitrogen pool state, and crop state; the parameter array is used to indicate the parameter data characterizing various reaction processes in the differentiable model, and the parameter data includes soil property data, crop data, and physical parameters output by the parameter generation network; the driving array is used to characterize meteorological data and agricultural management data; and the process mask is used to control the start and stop of any reaction process under any grid, any date, any soil layer, any crop type, and any water state type.
[0073] Parallel evolution at the grid dimension can be understood as performing time-step evolution simulations in parallel across multiple grids based on state data, parameter data, driving data, and process masks within multiple grids. This is equivalent to dividing the evolution simulation of the entire target farmland into evolution simulations of multiple spatial grid units, which helps improve the efficiency of evolution simulation. In practical applications, the batch_size can be pre-set to indicate the number of grids for parallel computation at one time, and the parallel evolution of multiple networks can be performed according to the number of grids. This disclosure does not limit this aspect.
[0074] It should be understood that the state array primarily represents various device-side arrays formed after state vectorization. These arrays represent the changing states of the differentiable model as the simulation progresses over time steps. The state array serves as both the input for a specific time step and the updated output after that time step, continuously updating in the time-step evolutionary simulation. The parameter array, on the other hand, represents the parameters used to control the physical, biochemical, or crop growth processes in the differentiable model, formed at the device level. These parameters can originate from soil parameter libraries, crop parameter libraries, model default parameters, or physical parameters obtained by converting the output of a parameter generation network through bounded and nonnegative mapping. Parameter arrays typically do not represent daily dynamic states but are used to control process rates, response functions, or physical constraints. For example, static soil properties such as soil clay content, field capacity, wilting point, and soil pH can be used as parameter arrays in calculations; decomposition rate parameters, nitrification / denitrification response parameters, crop physiological parameters, and water response parameters can also be used as parameter arrays, providing process control parameters for the model's forward calculations. A driving array refers to an array formed on the device from externally input meteorological sequence data or management event data. Its source is typically meteorological data, agricultural management records, etc. Driving arrays are not generated internally by the model but rather drive the model's progression over time. For example, daily temperature, precipitation, and radiation data, as well as two-dimensional arrays mapped from sowing, harvesting, fertilization, and irrigation measures, can all be categorized as driving arrays. Their function is to provide external boundary conditions or management events for the corresponding grid at each daily step. A process mask is a Boolean or 0 / 1 array used to mark whether a reaction process participates in the calculation at a specific time, grid, soil layer, or crop number. It can also be a continuous weighted array. Process masks are typically generated by conditions such as crop presence, management measures occurrence, soil layer effectiveness, grid effectiveness, flooding or cultivation status, and crop phenological stage during grid evolution. Their function is to control the start and stop of reaction processes under different grids, dates, soil layers, crop types, or moisture states, reducing grid-by-grid "if" statements and ensuring that all arrays maintain a consistent grid dimension during parallel computation on the device.
[0075] Therefore, maintaining a unified grid dimension for the state array, parameter array, driver array, and process mask can be understood as follows: regardless of dynamic states, static or physical parameters, external driver inputs, or flags indicating whether a reaction process is enabled, they are all organized according to the same set of spatial grid indices. Thus, on the th grid, the differentiable model can simultaneously read the corresponding state, parameters, driver data, and process mask for that grid and complete the parallel computation of that grid at the current time step. Furthermore, both single-time-step advancement and continuous multi-time-step advancement of the differentiable model can be accomplished through array slicing, broadcast computation, and parallel computing kernels. External programs can control the model evolution process through interfaces such as parameter writing, state reading, output snapshots, checkpoint saving and restoration. A shared view can be formed between the parameter generation network and the device-side array using device memory sharing or zero-copy tensor conversion.
[0076] In some embodiments, based on the reconstructing of the state vector into a device-side array, a general-purpose parallel computing platform (such as Compute Unified Device Architecture (CUDA)) can be used to perform device-side parallel computing of the differentiable model, thereby accelerating the forward evolution calculation of the differentiable model. Thus, the simulation results obtained by performing time-step evolutionary simulations using the differentiable model based on meteorological data, soil property data, agricultural management data, crop data, and physical parameters can include:
[0077] Using a general-purpose parallel computing platform, various reaction processes in a differentiable model are computed step-by-step in parallel at the grid dimension based on the state array, parameter array, driving array, and process mask, until the evolution termination condition is reached, yielding simulation results. In each time step, the state array is updated in parallel by the GPU thread within the general-purpose parallel computing platform. In other words, the parallel computing capabilities of a general-purpose parallel computing platform can be used to achieve the time-step parallel evolution of the aforementioned differentiable model based on the state array, parameter array, driving array, and process mask at the grid dimension.
[0078] Specifically, based on device-side arrays such as state arrays, parameter arrays, driver arrays, and process masks, as well as CUDA forward compute kernels, time-step parallel computation can be performed on various reaction process models in the differentiable model at the grid dimension, including soil hydrothermal processes, organic matter decomposition, nitrification, denitrification, crop growth, and gas emissions. In other words, data corresponding to multiple grids can be obtained at once from the state array, parameter array, driver array, and process mask, and time-step parallel computation of various reaction processes can be performed on the data corresponding to these multiple networks using CUDA forward compute kernels. In each time step, the state arrays corresponding to different grids, soil layers, or crop types can be updated in parallel by GPU threads, thereby reducing the computational overhead caused by the serial looping of the original biogeochemical model on a plot-by-plot and soil-by-soil-layer basis.
[0079] The evolution termination condition can be a condition customized by technicians according to actual needs, such as reaching a specified evolution duration or completing crop harvesting, etc., and this embodiment of the disclosure does not impose any limitations on this. As mentioned above, the state array represents the state data generated in the evolution simulation that changes with the progress of time steps. Therefore, the state array corresponding to each time step evolution can be updated in parallel by the image processor thread in the general parallel computing platform, which is beneficial to improving the efficiency of subsequent evolution.
[0080] It should be understood that the differentiable model can provide initial running conditions (or initial state) based on existing input files, default configurations or measured data during runtime, so as to start evolution simulation from the initial state and update the state array. This disclosure does not limit this.
[0081] Back Figure 1 In step S14, the target loss is determined based on the difference between the simulation results and the observation results for the target farmland. The observation results include one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions obtained from actual observations of the target farmland.
[0082] It should be understood that those skilled in the art can use observation methods known in the art or any source of observation data, such as ground observation, flux observation, remote sensing inversion or statistical yield data, to obtain the above-mentioned observation results for the target farmland, and this disclosure does not limit such observations.
[0083] In some embodiments, loss functions known in the art, such as mean squared error, can be used to determine the target loss based on the difference between the simulation results and the observation results for the target farmland (i.e., the observation error). For example, the target loss determined based on the difference between the simulation results and the observation results can be expressed as: L_obs = w_yield ||Y_yield - Y_yield_obs||²+w_N2O ||Y_N2O - Y_N2O_obs||²+w_SOC ||Y_SOC - Y_SOC_obs||², where Y_yield represents the simulated crop yield, Y_yield_obs represents the observed crop yield, w_yield represents the weight corresponding to the crop yield, Y_N2O represents the simulated nitrous oxide gas flux (i.e., one of the greenhouse gas emission results), Y_N2O_obs represents the observed nitrous oxide gas flux, w_N2O represents the weight corresponding to the nitrous oxide N2O gas flux, Y_SOC represents the simulated organic carbon pool SOC state, and Y_SOC_obs represents the observed organic carbon pool state. It should be understood that the target loss can be constructed by combining the simulation results with one or more of the actual observed crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions, and this disclosure does not limit this.
[0084] In some embodiments, the loss function can also be composed of an observation error term and a physical constraint term. That is, the target loss can be determined based on the difference between the simulation results and the observation results, as well as at least one physical constraint on the input and output. Thus, the loss function can be expressed as: L = L_obs + λ_prior L_prior + λ_nonneg L_nonneg + λ_cons L_cons, where L_obs represents the observation error term, used to measure the difference between the simulation results and the observation results; L_prior represents the parameter prior constraint term, used to limit the physical parameters output by the parameter generation network from deviating from empirical values or reasonable ranges; L_nonneg represents the non-negative constraint term, used to ensure that physical quantities such as carbon pool, nitrogen pool, crop biomass, and reaction flux are not less than zero; and L_cons represents the conservation constraint term, used to constrain the carbon, nitrogen, or water budget. The prior constraint term can be expressed as: L_prior = ||θ - θ_prior||², where θ represents the physical parameter, and θ_prior represents the empirical parameter, default parameter, or prior estimate of the physical parameter. This prior constraint term ensures that the physical parameters generated by the parameter generation network closely match actual experience. The nonnegation constraint term can be expressed as: L_nonneg = Σ ReLU(-S)², where S represents the state variable or flux that should remain nonnegative, such as the carbon and nitrogen pool state, crop biomass, and actual reaction flux during the simulation. The conservation constraint term can be expressed as: L_cons = ||Input - Output - State Change||², where input refers to the amount input during the reaction, output refers to the amount output during the reaction, and state change refers to the change in the state of the carbon and nitrogen pools. For example, for nitrogen processes, the balance between fertilizer input, mineralization input, nitrification / denitrification output, gas emissions, leaching losses, and nitrogen pool changes can be constrained; for carbon processes, the balance between residue input, decomposition output, respiration losses, and carbon pool changes can be constrained.
[0085] Back Figure 1 In step S15, the model parameters of the differentiable model and the network parameters of the optimized parameter generation network are optimized using the target loss, so as to perform online evolution simulation using the optimized differentiable model based on the physical parameters regenerated by the optimized parameter generation network.
[0086] As can be seen, biogeochemical models require input of the aforementioned physical parameters (which can be understood as dynamic parameters, depending on spatial, meteorological, and crop data), as well as a portion of static model parameters. Therefore, the model parameters of a differentiable model can be understood as the trainable static model parameters other than the aforementioned physical parameters. These static model parameters may be empirical or prior values, and cannot be generated by a parameter generation network combined with actual data (spatial, meteorological, and crop data). They may not have a clear physical or biogeochemical meaning, but they affect the computational process. For example, the sharpness parameter in the smoothing truncation function smooth_clip(x, min, max, sharpness) used in the aforementioned continuous processing, or the function parameters (or coefficients, factors, etc.) of various response functions, or the initialization parameters or default configurations required for model operation.
[0087] Understandingly, differentiability does not alter the physical equation structure of the original biogeochemical model. Instead, it enables the model's output simulation results to have a computable gradient relative to the input physical parameters. This allows the target loss to be passed to the parameter generation network through the differentiable model, thereby updating the network parameters of both the parameter generation network and the differentiable model. In other words, based on differentiable gradient optimization, the model parameters in the differentiable model and the network parameters of the parameter generation network can be optimized as a whole. By optimizing the network parameters of the parameter generation network, the physical parameters can be regenerated using the optimized parameter generation network, which is equivalent to optimizing some of the physical parameters in the differentiable model. At the same time, the target loss is used to optimize another part of the model parameters in the differentiable model, thus achieving optimization of all parameters in the entire differentiable model. This allows the differentiable model to perform accurate online evolution simulations based on the optimized physical and model parameters.
[0088] In practical applications, backpropagation and gradient descent algorithms known in the art can be used to optimize the model parameters of a differentiable model using the target loss and the network parameters of the optimization parameter generation network. Specifically, the above-mentioned optimization of the model parameters of a differentiable model using the target loss and the network parameters of the optimization parameter generation network may include:
[0089] The target loss gradient relative to the simulation results is backpropagated through the differentiable model to the parameter generation network to update the network parameters of the parameter generation network and the model parameters of the differentiable model.
[0090] New physical parameters are generated using the updated parameter generation network, and simulation results are regenerated based on the new physical parameters using the updated differentiable model until the difference between the simulation results regenerated by the updated differentiable model and the observation results is minimized, thus obtaining the optimized parameter generation network and the optimized differentiable model.
[0091] As mentioned above, differentiability enables the model output simulation results to have a computable gradient with respect to the input physical parameters. In other words, the simulation results output by a differentiable model (e.g., crop yield, N2O flux, SOC state, etc.) can be differentiated with respect to the input physical parameters, meaning it's possible to calculate: change in simulation result / change in parameter. For example, the gradient of simulated N2O flux with respect to nitrification response parameters or water response parameters can be calculated. This allows us to determine how an increase or decrease in a certain parameter affects the error when there is an error between the simulation and observation results. The gradient of the target loss relative to the simulation result is backpropagated through a differentiable model to the parameter generation network to update the gradient chain of the parameter generation network. This can be understood as: target loss L → simulation result Y output by the differentiable model → physical parameter θ input by the differentiable model → network parameter W of the parameter generation network. That is: dL / dW = dL / dY × dY / dθ × dθ / dW. Here, dY / dθ must be obtained through a differentiable model. If the model is not differentiable, it is impossible to know what impact the output θ of the parameter generation network has on the final simulation error, and therefore it is impossible to correctly update the network parameters of the parameter generation network used to generate θ.
[0092] It should be understood that during the backpropagation of the target loss gradient relative to the simulation results to the parameter generation network through the differentiable model, the target loss gradient relative to the simulation results is also transmitted to the differentiable model, thereby updating the model parameters of the differentiable model. In other words, this embodiment of the disclosure outputs the physical parameters of the differentiable model through the parameter generation network; the differentiable model receives these physical parameters and performs forward simulation; the loss function calculates the target loss based on the simulation results and observation results; the target loss gradient relative to the simulation results is transmitted back to the parameter generation network through the differentiable model to update the network parameters of the parameter generation network and the model parameters in the differentiable model other than the physical parameters; since the differentiable model is differentiable, the loss gradient is transmitted to the parameter generation network through the differentiable model, achieving end-to-end optimization of the differentiable model and the parameter generation network.
[0093] It should be understood that the above-described update process for the differentiable model and parameter generation network can be iteratively executed in multiple rounds. In each round of the update process, new physical parameters can be generated using the parameter generation network updated in the previous round, and simulation results can be regenerated using the differentiable model updated in the previous round based on the new physical parameters. Then, the target loss between the simulation results generated in the current round and the observation results can be calculated, and the network parameters of the parameter generation network updated in the previous round and the model parameters of the differentiable model can be updated using the target loss calculated in the current round, until the difference between the simulation results regenerated by the differentiable model updated in a certain round and the observation results is minimized (that is, the target loss converges or is set to zero), thus obtaining the optimized parameter generation network and the optimized differentiable model.
[0094] Then, based on the optimized parameter generation network (i.e., the parameter generation network with optimized network parameters), more accurate and suitable physical parameters for the target farmland can be regenerated according to the spatial data, meteorological data, agricultural management data, and crop data of the target farmland obtained online (i.e., the physical parameters are inverted or corrected using the optimized parameter generation network). Then, the optimized differentiable model (i.e., the differentiable model with optimized model parameters) is used to perform online evolution simulation based on the meteorological data, soil property data, agricultural management data, crop data, and the regenerated physical parameters obtained online, so as to obtain accurate online simulation results for the target farmland.
[0095] Considering that directly rewriting complex farmland biogeochemical models process-by-process using a general automatic differential framework would require expanding numerous time-step advancements, grid state updates, and biochemical reaction processes into fine-grained computational graphs, and given the long model evolution chain, numerous state variables, and complex inter-process dependencies, a large number of intermediate states typically need to be saved during forward computation for gradient calculation during backpropagation. This can easily lead to increased operator scheduling overhead, excessive GPU memory usage, and decreased computational efficiency, making it difficult to meet the needs of regional-scale, multi-grid, long-term series simulations and iterative training. Therefore, this embodiment of the present disclosure, for masked simulations of long-term series in differentiable models, does not permanently save the complete computational graph for all time steps. Instead, it saves state checkpoints at set intervals. Each checkpoint in the saved checkpoint sequence includes at least a time point, soil hydrothermal state, carbon pool state, nitrogen pool state, crop state, physical parameters, driving data location, and process mask; where the driving data location represents the read position of the driving array. In this way, during backpropagation, the state can be recovered from adjacent checkpoints, and the forward process within that interval can be re-executed, reducing the GPU memory overhead caused by long-term caching of intermediate states.
[0096] Therefore, the method may further include: in the time-step evolution simulation of the differentiable model, saving checkpoints at preset time intervals to obtain a checkpoint sequence, and using the checkpoint sequence to perform gradient backpropagation. For example, checkpoint saving can be performed every 100 time steps. The time point saved in the checkpoint refers to a specific simulation time point, usually corresponding to the current simulation day, day step number, or a day of the year. For example, if a checkpoint is saved on the 120th simulation day, then the time point saved in that checkpoint is 120. The purpose of saving the time point is to ensure that the corresponding meteorological data and management measures can be read from that day during reverse recalculation. Saving the physical parameters ensures that the parameters used during reverse recalculation are consistent with those used in the forward calculation, and calculates the gradient of the loss function on these parameters. The driving data location can be understood as the external input data index corresponding to the checkpoint, such as the location of the meteorological data array or management measures array corresponding to the current simulation day. Saving the driving data location ensures that the meteorological and management inputs corresponding to the same time period and the same grid can be reread during reverse recalculation, ensuring that the recalculation process is consistent with the original forward process. The purpose of preserving soil hydrothermal, carbon pool, nitrogen pool, and crop status data at checkpoints is to ensure that the differentiable model can be restored to the state corresponding to the checkpoint during reverse recalculation, thus guaranteeing consistency between the recalculation process and the original forward process. Similarly, preserving the process mask also ensures that the same reaction processes in the differentiable model as in the forward calculation can be restored during reverse recalculation, thereby guaranteeing consistency between the recalculation process and the original forward calculation process.
[0097] Therefore, based on the checkpoint sequence stored during the forward evolution simulation, the aforementioned backpropagation of the target loss relative to the simulation result's loss gradient to the parameter generation network via a differentiable model includes:
[0098] Starting from the checkpoint in the checkpoint sequence that is closest to the last time step in the evolution simulation, perform the following calculations sequentially until all checkpoints in the checkpoint sequence have been traversed:
[0099] The differentiable model is restored to the current checkpoint, and the evolution simulation within the preset time interval is re-executed based on the current checkpoint to obtain the forward computation path, which represents the complete computation process in the evolution simulation within the preset time interval.
[0100] Gradient calculation is performed based on the forward computation path to obtain the local gradient. The target gradient is determined based on the local gradient and the loss gradient, and the target gradient is passed to the previous checkpoint, the differentiable model, and the parameter generation network. The previous checkpoint is the checkpoint saved before the current checkpoint in the checkpoint sequence.
[0101] In this context, the checkpoint closest to the last time step in the evolutionary simulation can be understood as the last saved checkpoint in the checkpoint sequence, and the current checkpoint is the checkpoint at which recalculation is currently performed. Restoring the differentiable model to the current checkpoint means rolling back the differentiable model to the state at that current checkpoint, so as to re-execute the forward evolutionary simulation within the interval indicated by the preset time interval from the state at the current checkpoint. In other words, recalculation is performed to obtain the complete computation process in the evolutionary simulation within the preset time interval, that is, the complete and detailed forward computation process within the preset time interval (such as the inputs, outputs, and operators of each step of the calculation).
[0102] It should be understood that those skilled in the art can use gradient calculation methods known in the art to perform gradient calculation based on the forward calculation path, that is, calculate the derivative of the output with respect to the input to obtain the local gradient within the preset time interval; then, based on the chain rule, the local gradient dY / dθ and the loss gradient dL / dY of the current checkpoint can be accumulated (that is, the local gradient of the current checkpoint and the loss gradient are multiplied and then added to the target gradient passed from the upstream checkpoint) to obtain the target gradient corresponding to the current checkpoint and pass it to the previous checkpoint (to continue to propagate the gradient to an earlier time), the differentiable model (to calculate the gradient of the target loss value with respect to each model parameter of the differentiable model and update the trainable model parameters in the differentiable model) and the parameter generation network (so that the target gradient continues to be backpropagated to the parameter generation network, and based on the internal gradient dθ / dW of the parameter generation network, the gradient dL / dW of the target loss with respect to the internal network parameters of the parameter generation network is finally calculated to update the network parameters of the parameter generation network).
[0103] In some embodiments, considering that differentiable models contain computations of varying complexity—for example, some local response functions with stable formulas, well-defined derivatives, and frequent calls, or complex computational processes involving state coupling such as organic matter decomposition, nitrogen transformation, and crop growth—the aforementioned gradient calculation based on the forward computation path to obtain local gradients, in order to improve gradient calculation efficiency, includes:
[0104] For the locally defined response functions in the differentiable model, gradient calculation is performed using the pre-derived derivative formula of the locally defined response function based on the forward computation path to obtain the first type of local gradient; wherein, the locally defined response function includes at least: soil water-filling porosity response function, temperature response function, pH response function, smoothing continuous logic, and smoothing minimum logic; the derivative formula is implemented as a custom gradient kernel on a general parallel computing platform;
[0105] For complex reaction processes in differentiable models, local computational graphs corresponding to these complex reaction processes are established based on forward computational paths, and gradient calculations are performed based on the local computational graphs to obtain the second type of local gradients; wherein, the complex reaction processes include at least: organic matter decomposition process, nitrogen transformation process, and crop growth process.
[0106] Among them, determining the target gradient based on local gradients and loss gradients includes: combining and accumulating the first type of local gradient, the second type of local gradient, and the loss gradient based on the chain rule to obtain the target gradient.
[0107] Local response functions typically describe the continuous response relationship of a specific local process in a differentiable model. Examples include the WFPS response function (used to regulate the nitrification rate based on soil moisture conditions), the temperature response function (used to determine the effect of different soil temperatures on nitrifying bacteria activity, simulating how temperature changes accelerate or decelerate the nitrification reaction rate), the pH response function (used to simulate the promoting or inhibiting effect of soil pH on microorganisms, especially nitrifying bacteria activity), and the smooth_clip and smooth_min smoothing minimum functions used in the aforementioned continuous processing. These functions usually have few inputs, clear formulas, and act on a specific local process, such as controlling the rates of nitrification, denitrification, or organic matter decomposition. This disclosure does not limit the scope of these functions. Because these local response functions have well-defined formulas and require relatively few inputs, their derivative formulas can be derived in advance. For example, if the nitration reaction includes a WFPS smoothing response function, the derivative formula of that response function with respect to WFPS or parameters can be directly derived. These derivative formulas of the local response functions can then be implemented as custom gradient kernels (such as custom CUDA gradient kernels) on a general-purpose parallel computing platform. These gradient kernels are used to perform the calculation of the derivative formulas, allowing for parallel computation of the derivative formulas within the image processor. Therefore, based on the inputs and outputs of the local response functions used in the forward computation path, gradient calculation can be performed directly using the derivative formulas of the local response functions (i.e., by calling custom gradient kernels on a general-purpose parallel computing platform) to obtain the first type of local gradient.
[0108] For complex state-coupled reaction processes in differentiable models, such as organic matter decomposition, nitrogen transformation, and crop growth, the derivative formulas are difficult to derive directly due to the multiple variables involved and their mutual coupling. Therefore, a local computational graph can be established within the recalculation interval of the checkpoint. The automatic differentiation capability of this local computational graph can then be used to calculate the gradients of these complex reaction processes (i.e., automatically differentiating along the local computational graph), yielding the second type of local gradient. It should be understood that knowing the complete computational process represented by the forward computational path allows us to know the computational path of the complex reaction processes within that path, thus enabling the construction of the corresponding local computational graph. These local computational graphs can be temporarily stored and automatically deleted after calculating the second type of local gradient.
[0109] In some embodiments, for complex state-coupled reaction processes in differentiable models, a second type of local gradient can be obtained by employing a modular analytical approximation within the recalculation interval of checkpoints. This means that the complex reaction process can be approximated as a mathematical function throughout the entire recalculation interval. During reaction propagation, the second type of local gradient of the complex reaction process within the recalculation interval is obtained by differentiating this approximated mathematical function. This modular approximation sacrifices the fine-grained gradient changes for each day within the interval, resulting in lower gradient calculation accuracy but higher speed.
[0110] It should be understood that gradients generated by different backpropagation paths (i.e., the first and second type of local gradients mentioned above) can converge and accumulate at the interval boundary. Specifically, based on the chain rule in the backpropagation algorithm, the first and second type of local gradients can be combined and accumulated with the loss gradient to obtain the target gradient. In other words, within the recalculation interval of the current checkpoint, the first and second type of local gradients generated by different processes (modules) can be combined, multiplied with the loss gradient, and then added with the target gradient passed from the upstream checkpoint to obtain the target gradient accumulated at the current checkpoint. The target gradient accumulated at the current checkpoint can continue to be passed to the previous checkpoint, the differentiable model, and the parameter generation network.
[0111] In practical applications, after all gradients in the entire checkpoint sequence have been backpropagated (i.e., all intervals have been rerun and accumulated), an optimizer known in the art (such as the Adam optimizer) can be called to update the model parameters and network parameters based on the accumulated total gradients and the learning rate. This disclosure does not limit this process.
[0112] This disclosure discloses embodiments that divide the long-term sequence reaction process in the minimizable model into modules such as soil hydrothermal processes, organic matter decomposition, nitrification, denitrification, and crop growth, and select different reverse paths based on the formula complexity of each module. For example, such as... Figure 4This illustrates an evolutionary simulation process for forward computation and hybrid backpropagation at a miniaturizable DNDC device end, as shown below. Figure 4 As shown, the host side is responsible for providing the daily sequence index, meteorological data-driven operation, managing the event matrix, physical parameters to be written, and process masks; the device side maintains various state arrays such as water state, easily decomposable residual carbon library rcvl, organic carbon library soc, ammonia nitrogen library nh4, nitrate nitrogen library no3, urea nitrogen library urea, and crop state crop satates; and uses DLPack technology for device memory sharing / tensor conversion. The forward simulation phase calls the N2O emission module (N2O_balloon), the organic matter decomposition module (dndc_decomposition), and forward kernels for hydrothermal and crop simulations to complete parallel daily progress. It outputs snapshots / observational comparisons (i.e., nitrate nitrogen (NO3) state, ammonia nitrogen (NH4) state, nitrous oxide (N2O) gas flux, ammonia (NH3) gas flux, crop yield, hydrothermal state, etc.) and saves output snapshots or state checkpoints. Each checkpoint saves the current daily sequence (i.e., time point), soil hydrothermal state, carbon pool, nitrogen pool, crop state, physical parameters, driving data index, and process mask. The forward phase can be implemented using CUDA C++ Raw. Kernels reconstructs core modules (such as heat conduction and biochemical reactions) for parallel computation on GPU threads, eliminating loop overhead generated in the original model. During the backpropagation phase, the recovery state is restored from the nearest checkpoint before the target output, and the corresponding time interval is recalculated. This means restoring the soil hydrothermal state, carbon and nitrogen pools, crop state, physical parameters, and driving data index of the checkpoint, and rereading the meteorological and management events for the corresponding time period based on the driving data index. Then, the forward process between the checkpoint and the target time point (i.e., within the preset time interval) is re-executed (i.e., the forward interval is recalculated). For local functions with well-defined formulas, analytical gradients (i.e., pre-derived derivative formulas) or device-side gradient kernels are directly used. For multi-state coupled processes such as organic matter decomposition, nitrogen conversion, and crop growth, a local computation graph is temporarily established within the recalculation interval to calculate the gradient of the output result with respect to the input physical parameters within that interval, or a modular analytical approximation is used. This avoids writing complex inverse kernel functions by hand and significantly reduces GPU memory usage. Finally, after the gradients of each module converge at the interval boundary, the gradients generated by different modules are accumulated and continued to be passed to the previous checkpoint, the differentiable model, the physical parameters, and the parameter generation network.
[0113] According to the method of this disclosure, the differentiable model obtained by continuously processing the discrete logic in the reaction process of the biogeochemical model can form an end-to-end coupled optimization with the parameter generation network. It can also allow the target loss gradient to pass through the differentiable model, guiding the parameter generation network and the parameter update of the differentiable model. This allows for the rapid and accurate acquisition of an optimized parameter generation network and an optimized differentiable model adapted to the target farmland. In this way, the optimized parameter generation network can be used to automatically calibrate the physical parameters used to control the differentiable model. Furthermore, some model parameters in the optimized differentiable model are more accurate, which is equivalent to achieving parameter calibration for various parameters in the differentiable model. As a result, the optimized differentiable model and the optimized parameter generation network can be used to perform more accurate online evolution simulation of the target farmland, thereby improving the simulation accuracy.
[0114] Based on the methods proposed in the embodiments of this disclosure above, Figure 5 This diagram illustrates the architecture of an evolutionary simulation system that couples a differentiable model and a parameter generation network according to an embodiment of the present disclosure. Figure 5 As shown, the input data sources include meteorological data, soil property data, crop data, and agricultural management data. Optionally, remote sensing data may also be included. The input data undergoes data preprocessing (such as data cleaning) and feature encoding (i.e., the parameter generation network encodes the preprocessed input data) to obtain spatial features, time series features, crop features, and other data features. These data features are then input into the parameter generator (i.e., the parameter generation network decodes the extracted data features) to obtain physical parameters, such as WFPS response parameters, temperature / pH response parameters, decomposition rate parameters, crop parameters, water response parameters, and nitrification / denitrification related parameters. After bounded or non-negative mapping, these parameters are written into the device-side parameter array. The device-side forward engine of the differentiable model advances daily based on the state array and parameter array, outputting simulation results such as crop yield, gas flux, soil hydrothermal state, and carbon and nitrogen pool state. The simulation results and observation data constitute a loss function. The loss gradient is fed back through the simulation results, physical parameters, and parameter generation network output by the differentiable model to update the network parameters of the parameter generation network and / or the aforementioned physical parameters, and can also update the model parameters of the differentiable model.
[0115] According to the method of this disclosure, a farmland biogeochemical model, taking DNDC as an example, is endowed with "learning ability" for the first time. The non-differentiable problem is solved through mathematical reconstruction, achieving gradient-based data-driven parameter inversion. The heterogeneous fusion operator architecture achieves a 10-50 times speed improvement and reduces memory usage by more than 60% compared to pure PyTorch, making long-duration training at the regional level possible. The fused model (i.e., the differentiable model and the parameter generation network) combines the interpretability of physical models (following mass conservation) with the fitting ability of machine learning (handling nonlinearity), significantly improving simulation accuracy.
[0116] According to the method of this disclosure, a continuous simulation process of farmland biogeochemical processes based on smoothing operators is achieved, that is, a method that uses a custom smoothing function (such as smooth_clip) to replace the discrete decision logic in soil carbon and nitrogen cycling and water transport processes (such as WFPS threshold and temperature response) in agricultural models. Furthermore, a system coupling machine learning with differentiable farmland biogeochemical data assimilation and parameter inversion is implemented, that is, a system architecture that uses neural networks to process agricultural spatiotemporal data to predict physical parameters and performs end-to-end joint optimization by calculating the target loss through the physical layer of a differentiable agricultural biogeochemical model. Also, a heterogeneous differentiable operator construction architecture for large-scale farmland biogeochemical grid simulation is implemented, that is, a hybrid computing architecture that uses DLPack for zero-copy and calls CUDA computing kernels for parallel computation during forward propagation, and uses automatic differentiation or recomputation strategies for gradient backpropagation during backpropagation. Furthermore, a vectorized reconstruction method for farmland biogeochemical models oriented towards vector computation was implemented, which is a method for reconstructing traditional agricultural models into script-controlled parallel vector architectures. This method includes: mapping discrete objects to memory vectors, vectorizing physical process logic, and a dynamic binding mechanism between the interpreter and the underlying CUDA kernel functions.
[0117] Figure 6 A block diagram of a farmland biogeochemical process simulation apparatus according to an embodiment of the present disclosure is shown, such as Figure 6 As shown, the device includes:
[0118] The data acquisition module 601 is used to acquire spatial data, meteorological data, agricultural management data and crop data of the target farmland. The spatial data includes at least soil attribute data, topographic data and land use data of the target farmland.
[0119] The parameter generation module 602 is used to generate physical parameters based on the spatial data, the meteorological data, and the crop data using a machine learning-based parameter generation network. The physical parameters are used to control the reaction process in the differential model. The differential model is a model obtained by continuously processing the discrete logic in the reaction process of a general biogeochemical model. The biogeochemical model is used to simulate carbon and nitrogen conversion, greenhouse gas emissions, and crop growth in an agricultural ecosystem.
[0120] The evolution simulation module 603 is used to perform time-step evolution simulation using the differentiable model based on the meteorological data, the soil property data, the agricultural management data, the crop data and the physical parameters to obtain simulation results. The simulation results include at least one or more of the following: crop yield of the target farmland, soil carbon and nitrogen pool status, soil hydrothermal status and greenhouse gas emission results.
[0121] The loss determination module 604 is used to determine the target loss based on the difference between the simulation results and the observation results for the target farmland. The observation results include one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions obtained from actual observations of the target farmland.
[0122] The model optimization module 605 is used to optimize the model parameters of the differentiable model and the network parameters of the parameter generation network using the target loss, so as to perform online evolution simulation using the optimized differentiable model based on the physical parameters regenerated by the optimized parameter generation network.
[0123] In some embodiments, the reaction process includes at least: soil organic matter decomposition, nitrification and denitrification, crop growth, and soil hydrothermal processes; the discrete logic includes: environmental response logic and inventory limit constraint logic; the environmental response logic characterizes the threshold control logic of environmental control variables on the process rate of the reaction process, the threshold control logic being manifested as a sudden change or truncation of the process rate when the environmental control variable reaches the threshold; the process rate includes at least: nitrification rate, denitrification rate, and organic matter decomposition rate; the environmental control variables include at least: soil water porosity, temperature, and soil pH or pH state; the inventory limit constraint logic includes logic that limits the actual reaction flux during the reaction process to not exceed the corresponding inventory limit; the actual reaction flux includes at least: nitrification flux, denitrification flux, and organic matter decomposition flux; wherein, for the discrete logic of the reaction process in a general biogeochemical model... The continuous processing of the logic includes: replacing the environmental response logic with a smoothed continuous logic constructed using a smoothing function, wherein the smoothed continuous logic is used to make the process rate change continuously around a threshold; replacing the inventory upper limit constraint logic with flux scaling logic with non-negative constraints or smoothed minimum logic; wherein the flux scaling logic with non-negative constraints is used to scale the theoretical reaction flux using a scaling factor constructed based on the inventory upper limit, so that the scaled actual reaction flux is gradually limited by the inventory upper limit and does not exceed the inventory upper limit; the smoothed minimum logic is used to make the actual reaction flux approximately the minimum of the inventory upper limit and the theoretical reaction flux, wherein the theoretical reaction flux is the theoretical value generated by the reaction process in the differentiable model, and the actual reaction flux is the actual value after the inventory upper limit constraint logic imposes an upper limit constraint on the theoretical reaction flux.
[0124] In some embodiments, the differentiable model is obtained by reconstructing the dispersed state variables in the biogeochemical model into a specified array format; wherein, reconstructing the dispersed state variables in the biogeochemical model into a specified array format includes: reconstructing the state variables used to characterize soil attribute data into a two-dimensional array of soil layer dimension × grid dimension; reconstructing the state variables used to characterize water state into a three-dimensional array of water state type dimension × soil layer dimension × grid dimension; reconstructing the state variables used to characterize carbon pool state and nitrogen pool state into a two-dimensional array of soil layer dimension × grid dimension; reconstructing the state variables used to characterize crop state into a two-dimensional array of crop type dimension × grid dimension; reconstructing the state variables used to characterize agricultural management measures into a two-dimensional array of date dimension × grid dimension; wherein, the soil layer dimension represents different soil layers divided downwards from the surface in the vertical dimension, the grid dimension represents different geographical locations, plots, subdivided sampling points or raster cells in the horizontal spatial dimension; the water state type dimension represents different types of water state; and the crop type dimension characterizes different types of crops.
[0125] In some embodiments, during time-step evolutionary simulations of the differentiable model, the differentiable model undergoes parallel time-step evolution in a grid dimension based on a state array, a parameter array, a driving array, and a process mask; the state array, the parameter array, the driving array, and the process mask maintain a uniform grid dimension; wherein, the state array is used to characterize the state data generated by the differentiable model in the evolutionary simulation that changes with time steps, the state data including water state, carbon pool state, nitrogen pool state, and crop state; the parameter array is used to indicate the parameter data characterizing various reaction processes in the differentiable model, the parameter data including soil property data, crop data, and physical parameters output by the parameter generation network; the driving array is used to characterize the meteorological data and the agricultural management data; the process mask is used to control the start and stop of any reaction process under any grid, any date, any soil layer, any crop type, and any water state type.
[0126] In some embodiments, the soil property data includes at least: soil volume, soil clay content, field water holding capacity, wilting point, soil pH, and soil structure; the meteorological data includes at least: temperature, precipitation, humidity, wind speed, and radiation varying over time; the agricultural management data includes management measures taken for the target farmland from sowing, fertilization, irrigation to harvest; the crop data includes: the type of crop to be sown in the target farmland and crop growth information; wherein, the step of using the differentiable model to perform time-step evolutionary simulation based on the meteorological data, soil property data, agricultural management data, crop data, and physical parameters to obtain simulation results includes: using a general-purpose parallel computing platform to perform time-step parallel computation on various reaction processes in the differentiable model in the grid dimension according to the state array, the parameter array, the driving array, and the process mask until the evolution termination condition is reached to obtain simulation results, wherein the state array in each time step is updated in parallel by the image processor thread in the general-purpose parallel computing platform.
[0127] In some embodiments, optimizing the model parameters of the differentiable model and the network parameters of the parameter generation network using the target loss includes: backpropagating the loss gradient of the target loss relative to the simulation results through the differentiable model to the parameter generation network to update the network parameters of the parameter generation network and the model parameters of the differentiable model; generating new physical parameters using the updated parameter generation network, and regenerating the simulation results using the updated differentiable model based on the new physical parameters, until the difference between the simulation results regenerated by the updated differentiable model and the observation results is minimized, thereby obtaining the optimized parameter generation network and the optimized differentiable model.
[0128] In some embodiments, the apparatus further includes: a checkpoint saving module, configured to save checkpoints at preset time intervals during time-step evolutionary simulation of the differentiable model to obtain a checkpoint sequence, and to use the checkpoint sequence to perform gradient backpropagation; wherein each checkpoint in the checkpoint sequence includes at least: a time point, the corresponding soil hydrothermal state, carbon pool state, nitrogen pool state, crop state, physical parameters, driving data location, and process mask; the driving data location represents the read position of the driving array; wherein, backpropagating the target loss relative to the simulation result as a gradient to the parameter generation network through the differentiable model includes: extracting the gradient from the checkpoint sequence and... In the evolutionary simulation process, starting from the nearest checkpoint at the last time step, the following calculations are performed sequentially until all checkpoints in the checkpoint sequence are traversed: the differentiable model is restored to the current checkpoint, and the evolutionary simulation within a preset time interval is re-executed based on the current checkpoint to obtain a forward computation path, which represents the complete computation process in the evolutionary simulation within the preset time interval; gradient calculation is performed based on the forward computation path to obtain local gradients, and a target gradient is determined based on the local gradients and the loss gradient, and the target gradient is passed to the previous checkpoint, the differentiable model, and the parameter generation network; the previous checkpoint is the checkpoint saved before the current checkpoint in the checkpoint sequence.
[0129] In some embodiments, the step of calculating gradients based on the forward computation path to obtain local gradients includes: for the locally defined response functions in the differentiable model, calculating gradients based on the forward computation path using the derivative formula of the locally defined response functions to obtain a first type of local gradient; wherein the locally defined response functions include at least: soil water-filling porosity response function, temperature response function, pH response function, smoothing continuous logic, and smoothing minimum logic; the derivative formula is implemented as a custom gradient kernel on a general-purpose parallel computing platform; for complex reaction processes in the differentiable model, establishing local computation graphs corresponding to these complex reaction processes based on the forward computation path, and calculating gradients based on the local computation graphs to obtain a second type of local gradient; wherein the complex reaction processes include at least: organic matter decomposition process, nitrogen transformation process, and crop growth process; wherein the step of determining the target gradient based on the local gradients and the loss gradient includes: combining and accumulating the first type of local gradient, the second type of local gradient, and the loss gradient based on the chain rule to obtain the target gradient.
[0130] In some embodiments, the physical parameters include at least: organic matter decomposition rate parameters, crop parameters, nitrification rate parameters and denitrification rate parameters, water response parameters, temperature response parameters, acid-base response parameters, and soil water-filled porosity response parameters; the crop parameters include at least crop physiological parameters, phenological parameters, and biomass allocation parameters; wherein, after generating the physical parameters, the device further includes: a mapping module, used to map the physical parameters to corresponding physical constraint boundaries to obtain mapped physical parameters, so that the differentiable model can perform time-step evolution simulation based on the mapped physical parameters; wherein, mapping the physical parameters to corresponding physical constraint boundaries includes: mapping the water response parameters, temperature response parameters, acid-base response parameters, soil water-filled porosity response parameters, and crop parameters to their respective corresponding preset reasonable ranges; and mapping the organic matter decomposition rate parameters, nitrification rate parameters, and denitrification rate parameters to non-negative values.
[0131] In some embodiments, the functions or modules of the apparatus provided in this disclosure can be used to perform the methods described in the above method embodiments. The specific implementation can be referred to the description of the above method embodiments, and for the sake of brevity, it will not be repeated here.
[0132] This disclosure also provides an electronic device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the above method.
[0133] This disclosure also provides a non-volatile computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the steps of the above-described method.
[0134] This disclosure also provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the above-described method.
[0135] Figure 7 A block diagram of an electronic device 1900 according to an embodiment of the present disclosure is shown. For example, the electronic device 1900 may be provided as a server or a terminal device. (Refer to...) Figure 7 The electronic device 1900 includes a processing component 1922, which further includes one or more processors, and memory resources represented by memory 1932 for storing instructions, such as application programs, that can be executed by the processing component 1922. The application programs stored in memory 1932 may include one or more modules, each corresponding to a set of instructions. Furthermore, the processing component 1922 is configured to execute instructions to perform the methods described above.
[0136] Electronic device 1900 may also include a power supply component 1926 configured to perform power management of electronic device 1900, a wired or wireless network interface 1950 configured to connect electronic device 1900 to a network, and an input / output interface 1958 (I / O interface). Electronic device 1900 can operate on an operating system, such as Windows Server, stored in memory 1932. TM Mac OS X TM Unix TM Linux TM FreeBSD TM Or similar.
[0137] In an exemplary embodiment, a non-volatile computer-readable storage medium is also provided, such as a memory 1932 including computer program instructions that can be executed by a processing component 1922 of an electronic device 1900 to perform the above-described method.
[0138] Computer-readable storage media can be tangible devices capable of holding and storing programs / instructions used by instruction execution devices. Computer-readable storage media can be, for example—but not limited to—electrical storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of computer-readable storage media include: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), static random access memory (SRAM), portable compact disc read-only memory (CD-ROM), digital multifunction disc (DVD), memory sticks, floppy disks, mechanical encoding devices, such as punch cards or recessed protrusions storing instructions thereon, and any suitable combination of the foregoing. The computer-readable storage media used herein are not to be construed as transient signals themselves, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through waveguides or other transmission media (e.g., light pulses through fiber optic cables), or electrical signals transmitted through wires.
[0139] The computer program (or computer-readable program instructions) described herein can be downloaded from a computer-readable storage medium to various computing / processing devices, or downloaded via a network, such as the Internet, local area network, wide area network, and / or wireless network, to an external computer or external storage device. The network may include copper transmission cables, fiber optic transmission, wireless transmission, routers, firewalls, switches, gateway computers, and / or edge servers. A network adapter card or network interface in each computing / processing device receives the computer-readable program instructions from the network and forwards them to the computer-readable storage medium in the respective computing / processing device.
[0140] The computer program (or computer program instructions) used to perform the operations of this disclosure may be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, state setting data, or source code or object code written in any combination of one or more programming languages, including object-oriented programming languages such as Smalltalk, C++, etc., and conventional procedural programming languages such as the "C" language or similar programming languages. The computer-readable program instructions may execute entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving a remote computer, the remote computer may be connected to the user's computer via any type of network—including a local area network (LAN) or a wide area network (WAN)—or may be connected to an external computer (e.g., via the Internet using an Internet service provider). In some embodiments, electronic circuitry, such as programmable logic circuitry, field-programmable gate arrays (FPGAs), or programmable logic arrays (PLAs), is personalized by utilizing state information from the computer-readable program instructions to implement various aspects of this disclosure.
[0141] Various aspects of this disclosure are described herein with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this disclosure. It should be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer-readable program instructions.
[0142] These computer-readable program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing apparatus to produce a machine such that, when executed by the processor of the computer or other programmable data processing apparatus, they create means for implementing the functions / actions specified in one or more blocks of the flowchart and / or block diagram. These computer-readable program instructions can also be stored in a computer-readable storage medium that causes a computer, programmable data processing apparatus, and / or other device to operate in a particular manner; thus, the computer-readable medium storing the instructions comprises an article of manufacture that includes instructions for implementing aspects of the functions / actions specified in one or more blocks of the flowchart and / or block diagram.
[0143] Computer-readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable data processing apparatus, or other device to produce a computer-implemented process, thereby causing the instructions executed on the computer, other programmable data processing apparatus, or other device to perform the functions / actions specified in one or more boxes of a flowchart and / or block diagram.
[0144] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present disclosure. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of an instruction containing one or more executable instructions for implementing a specified logical function. In some alternative implementations, the functions marked in the blocks may occur in a different order than those shown in the drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in the block diagrams and / or flowcharts, and combinations of blocks in the block diagrams and / or flowcharts, may be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions.
[0145] The various embodiments of this disclosure have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical application, or technical improvements to the embodiments in the market, or to enable others skilled in the art to understand the embodiments disclosed herein.
Claims
1. A method for simulating biogeochemical processes in farmland, characterized in that, include: Acquire spatial data, meteorological data, agricultural management data, and crop data of the target farmland. The spatial data shall include at least soil attribute data, topographic data, and land use data of the target farmland. A parameter generation network based on machine learning is used to generate physical parameters based on the spatial data, the meteorological data, and the crop data. The physical parameters are used to control the reaction process in a differential model. The differential model is a model obtained by continuous processing of the discrete logic in the reaction process of a general biogeochemical model. The biogeochemical model is used to simulate carbon and nitrogen conversion, greenhouse gas emissions, and crop growth in an agricultural ecosystem. The differentiable model is used to perform time-step evolution simulation based on the meteorological data, soil property data, agricultural management data, crop data, and physical parameters to obtain simulation results. The simulation results include at least one or more of the following: crop yield of the target farmland, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emission results. The target loss is determined based on the difference between the simulation results and the observation results for the target farmland. The observation results include one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions obtained from actual observations of the target farmland. The model parameters of the differentiable model and the network parameters of the parameter generation network are optimized using the target loss, so as to perform online evolution simulation using the optimized differentiable model based on the physical parameters regenerated by the optimized parameter generation network.
2. The method according to claim 1, characterized in that, The reaction process includes at least: soil organic matter decomposition process, nitrification and denitrification process, crop growth process, and soil hydrothermal process; The discrete logic includes: environmental response logic and inventory upper limit constraint logic; the environmental response logic represents the threshold control logic of environmental control variables on the process rate of the reaction process, and the threshold control logic is manifested as the process rate undergoing a sudden change or truncation when the environmental control variable reaches the threshold; the process rate includes at least: nitrification rate, denitrification rate, and organic matter decomposition rate; the environmental control variables include at least: soil water porosity, temperature, and soil pH or pH state. The inventory limit constraint logic includes logic that restricts the actual reaction flux during the reaction process from not exceeding the corresponding inventory limit; the actual reaction flux includes at least: nitrification flux, denitrification flux, and organic matter decomposition flux; The process of making the discrete logic in the reaction process of a general biogeochemical model continuous includes: For the aforementioned environmental response logic, a smooth continuous logic constructed using a smoothing function is used to replace the environmental response logic. The smooth continuous logic is used to make the process rate change continuously around a threshold. For the aforementioned inventory upper limit constraint logic, use flux scaling logic with non-negative constraints or smoothing minimum value logic to replace the inventory upper limit constraint logic; The flux scaling logic with non-negative constraints is used to scale the theoretical reaction flux using a scaling factor constructed based on the inventory cap, so that the scaled actual reaction flux is gradually limited by the inventory cap and does not exceed the inventory cap; the smoothing minimum logic is used to make the actual reaction flux approximate the minimum of the inventory cap and the theoretical reaction flux, where the theoretical reaction flux is the theoretical value generated by the reaction process in the differentiable model, and the actual reaction flux is the actual value after the inventory cap constraint logic imposes an upper limit constraint on the theoretical reaction flux.
3. The method according to claim 1, characterized in that, The differentiable model is obtained by reconstructing the dispersed state variables in the biogeochemical model into a specified array form; The process of reconstructing the dispersed state variables in the biogeochemical model into a specified array format includes: reconstructing the state variables used to characterize soil attribute data into a two-dimensional array of soil layer dimension × grid dimension; reconstructing the state variables used to characterize water state into a three-dimensional array of water state type dimension × soil layer dimension × grid dimension; reconstructing the state variables used to characterize carbon pool state and nitrogen pool state into a two-dimensional array of soil layer dimension × grid dimension; reconstructing the state variables used to characterize crop state into a two-dimensional array of crop type dimension × grid dimension; and reconstructing the state variables used to characterize agricultural management measures into a two-dimensional array of date dimension × grid dimension. Wherein, the soil layer dimension represents different soil layers divided downwards from the ground surface in the vertical dimension; the grid dimension represents different geographical locations, plots, sample points, or grid cells in the horizontal spatial dimension; the moisture state type dimension represents different types of moisture states; and the crop type dimension represents different types of crops.
4. The method according to claim 3, characterized in that, In the time-step evolution simulation of the differentiable model, the differentiable model performs time-step parallel evolution based on the state array, parameter array, driving array, and process mask in the grid dimension; the state array, the parameter array, the driving array, and the process mask maintain a unified grid dimension; The state array is used to characterize the state data generated by the differentiable model in the evolutionary simulation that changes with the time step. The state data includes water state, carbon pool state, nitrogen pool state, and crop state. The parameter array is used to indicate parameter data characterizing various reaction processes in the differentiable model. The parameter data includes the soil property data, the crop data, and the physical parameters output by the parameter generation network. The driving array is used to represent the meteorological data and the agricultural management data; The process mask is used to control the start and stop of any reaction process under any grid, any date, any soil layer, any crop type, and any moisture state.
5. The method according to claim 4, characterized in that, The soil property data includes at least: soil volume, soil clay content, field capacity, wilting point, soil pH, and soil structure; the meteorological data includes at least: temperature, precipitation, humidity, wind speed, and radiation over time; the agricultural management data includes the management measures taken for the target farmland from sowing, fertilization, irrigation to harvest; and the crop data includes: the types of crops to be sown in the target farmland and crop growth information. The step of using the differentiable model to perform time-step evolutionary simulations based on the meteorological data, soil property data, agricultural management data, crop data, and physical parameters to obtain simulation results includes: Using a general-purpose parallel computing platform, various reaction processes in the differentiable model are computed in parallel step by step in the grid dimension according to the state array, the parameter array, the driving array and the process mask, until the evolution termination condition is reached, and the simulation results are obtained. In each time step, the state array is updated in parallel by the image processor thread in the general-purpose parallel computing platform.
6. The method according to claim 4, characterized in that, The process of optimizing the model parameters of the differentiable model and the network parameters of the parameter generation network using the target loss includes: The loss gradient of the target loss relative to the simulation result is backpropagated to the parameter generation network through the differentiable model to update the network parameters of the parameter generation network and the model parameters of the differentiable model. New physical parameters are generated using the updated parameter generation network, and simulation results are regenerated based on the new physical parameters using the updated differentiable model until the difference between the simulation results regenerated by the updated differentiable model and the observed results is minimized, thus obtaining the optimized parameter generation network and the optimized differentiable model.
7. The method according to claim 6, characterized in that, The method further includes: In the time-step evolution simulation of the differentiable model, checkpoints are saved at preset time intervals to obtain a checkpoint sequence, which is then used to perform gradient backpropagation during backpropagation. Each checkpoint in the checkpoint sequence includes at least: a time point, the corresponding soil hydrothermal state, carbon pool state, nitrogen pool state, crop state, physical parameters, driving data location, and process mask; the driving data location represents the read position of the driving array. Wherein, the step of backpropagating the loss gradient of the target loss relative to the simulation result to the parameter generation network through the differentiable model includes: Starting from the checkpoint in the checkpoint sequence that is closest to the last time step in the evolution simulation, perform the following calculations sequentially until all checkpoints in the checkpoint sequence have been traversed: The differentiable model is restored to the current checkpoint, and the evolution simulation within the preset time interval is re-executed based on the current checkpoint to obtain the forward computation path, which represents the complete computation process in the evolution simulation within the preset time interval. Gradient calculation is performed based on the forward computation path to obtain the local gradient. The target gradient is determined based on the local gradient and the loss gradient, and the target gradient is passed to the previous checkpoint, the differentiable model, and the parameter generation network. The previous checkpoint is the checkpoint saved before the current checkpoint in the checkpoint sequence.
8. The method according to claim 7, characterized in that, The step of calculating the gradient based on the forward computation path to obtain the local gradient includes: For the locally defined response functions in the differentiable model, gradient calculation is performed using the pre-derived derivative formula of the locally defined response function based on the forward computation path to obtain the first type of local gradient; wherein, the locally defined response function includes at least: soil water-filling porosity response function, temperature response function, pH response function, smoothing continuous logic, and smoothing minimum logic; the derivative formula is implemented as a custom gradient kernel on a general-purpose parallel computing platform; For the complex reaction process in the differentiable model, a local computation graph corresponding to the complex reaction process is established based on the forward computation path, and gradient calculation is performed based on the local computation graph to obtain a second type of local gradient; wherein, the complex reaction process includes at least: organic matter decomposition process, nitrogen transformation process, and crop growth process; The step of determining the target gradient based on the local gradient and the loss gradient includes: combining and accumulating the first type of local gradient, the second type of local gradient, and the loss gradient based on the chain rule to obtain the target gradient.
9. The method according to claim 1, characterized in that, The physical parameters include at least: organic matter decomposition rate parameters, crop parameters, nitrification rate parameters, denitrification rate parameters, water response parameters, temperature response parameters, acid-base response parameters, and soil water-filled porosity response parameters; the crop parameters include at least crop physiological parameters, phenological parameters, and biomass allocation parameters; wherein, after generating the physical parameters, the method further includes: The physical parameters are mapped to the corresponding physical constraint boundaries to obtain the mapped physical parameters, so that the differentiable model can perform time-step evolution simulation based on the mapped physical parameters. The step of mapping the physical parameters to the corresponding physical constraint boundaries includes: The moisture response parameters, temperature response parameters, acid-base response parameters, soil water-filling porosity response parameters, and crop parameters are each mapped to their respective preset reasonable ranges. The organic matter decomposition rate parameter, the nitrification rate parameter, and the denitrification rate parameter are each mapped to non-negative values.
10. A device for simulating farmland biogeochemical processes, characterized in that, include: The data acquisition module is used to acquire spatial data, meteorological data, agricultural management data and crop data of the target farmland. The spatial data includes at least soil attribute data, topographic data and land use data of the target farmland. The parameter generation module is used to generate physical parameters based on the spatial data, meteorological data, and crop data using a machine learning-based parameter generation network. The physical parameters are used to control the reaction process in the differential model, which is a model obtained by continuous processing of the discrete logic in the reaction process of a general biogeochemical model. The biogeochemical model is used to simulate carbon and nitrogen conversion, greenhouse gas emissions, and crop growth in an agricultural ecosystem. The evolution simulation module is used to perform time-step evolution simulations using the differentiable model based on the meteorological data, the soil property data, the agricultural management data, the crop data, and the physical parameters to obtain simulation results. The simulation results include at least one or more of the following: crop yield of the target farmland, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emission results. The loss determination module is used to determine the target loss based on the difference between the simulation results and the observation results for the target farmland. The observation results include one or more of the following: crop yield, soil carbon and nitrogen pool status, soil hydrothermal status, and greenhouse gas emissions obtained from actual observations of the target farmland. The model optimization module is used to optimize the model parameters of the differentiable model and the network parameters of the parameter generation network using the target loss, so as to perform online evolution simulation using the optimized differentiable model based on the physical parameters regenerated by the optimized parameter generation network.