Potato fertilizer automatic application control program based on digital twinning

Abstract

The present application relates to the field of digital twinning, in particular to a potato fertilizer automatic application control program based on digital twinning. The physiological characteristic parameters of potato growth period, soil moisture parameters and historical fertilization data are obtained and input into the pre-constructed coupling digital twin model of potato root growth morphology and soil three-dimensional nutrient transport; the growth trajectory of the root system in the three-dimensional soil space and the nutrient absorption interception radius of the next growth cycle are deduced forward; the spatial coincidence degree of the root interception domain and the current soil nutrient concentration distribution is compared to generate a nutrient deficiency spatial distribution map; and the same is mapped to the fertilizer execution mechanism to calculate the fertilizer discharge time delay and the fertilizer discharge pulse width with the root interception domain center as the target coordinate to generate a targeted fertilization control instruction. The present application matches the fertilizer application position with the actual spatial development morphology of the potato root system, overcomes the defects of fertilizer deviating from the effective absorption area of the root system, and improves the interception rate of fertilizer in the rhizosphere micro domain.

Description

Technical Field

[0001] This invention relates to the field of digital twins, specifically to an automatic fertilization control program for potatoes based on digital twins. Background Technology

[0002] Existing potato fertilization control programs mostly rely on variable fertilization based on surface soil monitoring data and the timeline of potato growth stages. During program execution, the fertilization actuator is controlled to apply fertilizer at a fixed opening or based on uniform nutrient concentration values ​​collected by surface soil sensors, according to a preset growth schedule and corresponding basal and topdressing ratio parameters. Fertilizer is directly applied to the surface or into single-depth furrows created by furrow openers. The location and timing of fertilization are limited by the monitoring results of the surface time dimension, failing to incorporate the dynamic development process of underground roots and the physical transport of soil nutrients into the closed-loop control of the program.

[0003] The aforementioned existing technologies suffer from a core technical problem: the fertilization control program cannot dynamically match the fertilizer application location with the actual spatial development morphology of the potato's underground root system, causing the fertilizer to deviate from the effective absorption zone of the roots. The growth trajectory and nutrient absorption interception range of potato roots in the three-dimensional soil space dynamically change with the soil physical environment and growth stage. Existing technologies, which implement fertilization control based on time nodes and surface concentrations, do not target the application based on the specific location and interception capacity of the roots in the underground space. As a result, the fertilizer remains in non-rhizosphere areas after application and cannot be intercepted and absorbed by the roots. Summary of the Invention

[0004] The purpose of this invention is to provide an automatic fertilization control program for potatoes based on digital twins, which can effectively solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0006] The digital twin-based automated fertilization control program for potatoes includes:

[0007] Obtain physiological parameters of potato growth period, soil moisture parameters, and historical fertilization data;

[0008] The potato growth period physiological characteristic parameters, the soil moisture parameters, and the historical fertilization data are input into a pre-constructed coupled digital twin model of potato root growth morphology and soil three-dimensional nutrient transport.

[0009] The coupled digital twin model is driven to forward extrapolate the growth trajectory and nutrient absorption interception radius of the potato root system in three-dimensional soil space in the next growth cycle based on the potato growth period physiological characteristic parameters, the soil moisture parameters and the historical fertilization data.

[0010] By comparing the spatial overlap between the root interception domain corresponding to the nutrient absorption interception radius in the coupled digital twin model and the current soil nutrient concentration distribution, a nutrient deficit spatial distribution map is generated.

[0011] The nutrient deficiency spatial distribution map is mapped to the walking path and fertilization opening of the fertilization actuator. Using the center of the root interception domain as the target coordinate, the fertilization time delay and fertilization pulse width of the fertilization actuator are calculated to generate targeted fertilization control commands.

[0012] Preferably, the acquisition of potato growth period physiological characteristic parameters, soil moisture parameters and historical fertilization data also includes acquiring soil bulk density parameters and soil porosity parameters in three-dimensional soil space;

[0013] The driving force of the coupled digital twin model, based on the potato growth period physiological characteristic parameters, soil moisture parameters, and historical fertilization data, forward extrapolates the growth trajectory and nutrient absorption interception radius of the potato root system in three-dimensional soil space for the next growth cycle, including:

[0014] The soil bulk density parameter and the soil porosity parameter are introduced into the coupled digital twin model to construct a root growth spatial resistance field;

[0015] Based on the root tropism growth trend determined by the physiological characteristic parameters of potato growth period, the optimal descent path of the root tip growth potential energy is solved in the root growth space resistance field, and the optimal descent path is determined as the growth trajectory.

[0016] Based on the soil porosity parameters and soil moisture parameters of the region where the growth trajectory is located, the root hair branch density is dynamically adjusted, and the nutrient absorption interception radius is calculated based on the root hair branch density.

[0017] Preferably, the coupled digital twin model embeds soil solute transport convection-dispersion equations;

[0018] By comparing the spatial overlap between the root interception domain corresponding to the nutrient absorption interception radius in the coupled digital twin model and the current soil nutrient concentration distribution, a nutrient deficit spatial distribution map is generated, including:

[0019] The historical fertilization data and the soil moisture parameters are input into the soil solute transport convection-diffusion equation to calculate the time-varying nutrient concentration gradient of each spatial voxel in the three-dimensional soil space at the current moment.

[0020] Extract the target voxel set covered by the root interception domain, and calculate the spatial distribution of the difference between the integral value of the time-varying nutrient concentration gradient within the target voxel set and the potato target fertilizer concentration threshold.

[0021] The voxel regions with differences greater than zero in the difference space distribution are extracted as nutrient deficiency spaces. The nutrient deficiency spaces are then subjected to three-dimensional gridded interpolation fitting to generate the nutrient deficiency space distribution map.

[0022] Preferably, the nutrient deficit spatial distribution map is mapped to the walking path and fertilization opening of the fertilization execution mechanism, and the fertilization time delay of the fertilization execution mechanism is calculated using the center of the root interception domain as the target coordinate, including:

[0023] The current walking speed of the fertilization actuator and the soil penetration depth of the furrow opener are obtained, and the relative offset of the target coordinates in the horizontal direction and the target fertilization depth in the vertical direction are determined by combining the nutrient deficiency spatial distribution map.

[0024] The gravity settling time required for fertilizer particles to fall is calculated based on the soil penetration depth of the trencher and the target fertilization depth.

[0025] Calculate the mechanism displacement time required for the fertilizer application actuator to reach the horizontal projection position of the target coordinates based on the current walking speed and the relative offset;

[0026] The sum of the gravity settling time and the mechanism displacement time is set as the fertilizer discharge time delay, so that the fertilizer application mechanism performs the fertilizer discharge action according to the fertilizer discharge time delay when it reaches the target coordinate.

[0027] Preferably, calculating the fertilizer discharge pulse width of the fertilizer application actuator includes:

[0028] Extract the difference between the deficit volume and the deficit concentration of the deficit grid to which the target coordinate belongs in the nutrient deficit spatial distribution map, and calculate the initial number of fertilizer discharge pulses by combining the preset single-pulse fertilizer discharge calibration value.

[0029] Obtain the variance of the discharge velocity fluctuation at the discharge port of the fertilization actuator, input the variance of the discharge velocity fluctuation into the pre-constructed pulse width modulation compensation model, and output the pulse width correction coefficient for the initial number of discharge pulses;

[0030] Multiply the initial number of fertilizer discharge pulses, the reference pulse width corresponding to the variance of fertilizer discharge velocity fluctuation, and the pulse width correction coefficient to obtain the target fertilizer discharge pulse width for the target coordinates;

[0031] The fertilization actuator controls the opening and closing duration of the fertilization valve according to the target fertilization pulse width, and delivers the corresponding amount of fertilizer to the target coordinate position.

[0032] Preferably, after generating the targeted fertilization control instruction, the method further includes:

[0033] The real-time feedback data of the actual fertilizer discharge from the fertilization execution mechanism and the real-time multispectral nutrient detection data of the soil profile after fertilization are obtained.

[0034] The real-time fertilizer discharge feedback data and the real-time multispectral nutrient detection data are input into the coupled digital twin model to drive the coupled digital twin model to update the current soil nutrient concentration distribution and obtain the actual soil nutrient distribution state after fertilization.

[0035] Compare the actual distribution of soil nutrients after fertilization with the target nutrient distribution corresponding to the spatial distribution map of nutrient deficit, and calculate the spatial concentration residual.

[0036] If the spatial concentration residual exceeds a preset residual threshold, the root growth rate parameter and soil nutrient diffusion coefficient in the coupled digital twin model are corrected based on the spatial concentration residual, and the targeted fertilization control command is regenerated based on the corrected coupled digital twin model.

[0037] Preferably, the soil bulk density parameter and the soil porosity parameter are introduced into the coupled digital twin model to construct a root growth spatial resistance field, including:

[0038] The distribution value of soil penetration resistance is calculated based on the soil bulk density parameter and the soil porosity parameter;

[0039] The soil penetration resistance distribution value is compared with the root critical penetration stress to delineate the root impenetrable boundary;

[0040] The method of dynamically adjusting the root hair branch density based on the soil porosity parameter and soil moisture parameter of the region where the growth trajectory is located, and calculating the nutrient absorption interception radius based on the root hair branch density, includes:

[0041] Extract the path length of the growth trajectory through the pores between the impenetrable boundaries of adjacent roots, and calculate the seepage pressure on the surface of the main root based on the path length and the soil moisture parameters.

[0042] When the seepage pressure on the main root surface exceeds the osmotic pressure threshold of the root cortex cells, the density of the hairy root branches is triggered to increase exponentially along the seepage pressure gradient direction on the main root surface. The radial radius of the increased hairy root branch density is set as the nutrient absorption interception radius.

[0043] Preferably, the historical fertilization data and the soil moisture parameters are input into the soil solute transport convection-diffusion equation to calculate the time-varying nutrient concentration gradient of each spatial voxel in the three-dimensional soil space at the current moment, including:

[0044] The isothermal adsorption and desorption parameters of soil particles for nutrient molecules are obtained, and the isothermal adsorption and desorption parameters are introduced as source and sink terms into the soil solute transport convection-dispersion equation.

[0045] By combining the water infiltration flux in the soil moisture parameters, the soil solute transport convection-dispersion equation after introducing the source-sink term is solved to obtain the time-varying nutrient concentration gradient of free nutrient molecules in each spatial voxel;

[0046] The calculation of the spatial distribution of the difference between the integral value of the time-varying nutrient concentration gradient within the target voxel set and the target fertilizer requirement threshold for potatoes includes:

[0047] Within the target voxel set, the time-varying concentration gradient of the free nutrient ions is integrally divided along the normal direction of the nutrient absorption interception radius to obtain the effective nutrient integral that the root system can actually absorb. The spatial difference between the effective nutrient integral and the potato target fertilizer concentration threshold is calculated, and the spatial difference is mapped back to a three-dimensional mesh to generate the spatial distribution of the difference.

[0048] Preferably, the calculation of the gravity settling time required for fertilizer particles to fall, based on the trencher's penetration depth and the target fertilization depth, includes:

[0049] Obtain the particle size distribution parameters and particle density parameters of fertilizer particles, and construct a particle settling resistance model by combining the soil moisture parameters between the trencher's depth of penetration and the target fertilization depth.

[0050] The particle size distribution parameter, the particle density parameter, and the soil moisture parameter are input into the particle settling resistance model to calculate the terminal settling velocity of fertilizer particles in the soil gaps. The gravity settling time is calculated based on the difference between the terminal settling velocity and the soil penetration depth of the furrow opener and the target fertilization depth.

[0051] The calculation of the mechanism displacement time required for the fertilization actuator to reach the horizontal projection position of the target coordinates based on the current walking speed and the relative offset includes:

[0052] The vibration frequency of the fertilizing machinery chassis and the elevation difference of the ground surface undulation are obtained. The current walking speed is corrected based on the vibration frequency of the chassis and the elevation difference of the ground surface undulation to obtain the actual effective walking speed. The displacement time of the mechanism is calculated based on the actual effective walking speed and the relative offset.

[0053] Preferably, the variance of the discharge velocity fluctuation at the discharge port of the fertilization actuator is obtained, and the variance of the discharge velocity fluctuation is input into a pre-constructed pulse width modulation compensation model to output a pulse width correction coefficient for the initial number of discharge pulses, including:

[0054] The fluctuation values ​​of the discharge shaft speed and the change value of the discharge port opening gap at the discharge port are collected in real time, and the variance of the discharge flow velocity fluctuation is calculated based on the fluctuation values ​​of the discharge shaft speed and the change value of the discharge port opening gap.

[0055] The variance of the fertilizer discharge velocity fluctuation is matched and queried with a preset variance-coefficient mapping table. The variance-coefficient mapping table is constructed based on the nonlinear velocity compensation coefficients corresponding to different variance intervals. The pulse width correction coefficient is obtained by querying the table.

[0056] The fertilization actuator controls the opening and closing duration of the fertilization valve according to the target fertilization pulse width, delivering the corresponding amount of fertilizer to the target coordinate position, including:

[0057] The target fertilizer discharge pulse width is converted into a control level signal for the fertilizer discharge valve. When the fertilizer discharge time delay arrives, the control level signal is output to drive the fertilizer discharge valve to open for the specified opening and closing duration, so that the fertilizer particles are transported to the target coordinate position under the combined action of gravity and forced airflow.

[0058] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0059] 1. This invention constructs a coupled digital twin model of potato root growth morphology and three-dimensional nutrient transport in soil. Potato physiological parameters during the growth period, soil moisture parameters, and historical fertilization data are input into the model. The model forward extrapolates the growth trajectory and nutrient absorption interception radius of the potato root system in three-dimensional soil space for the next growth cycle. The spatial overlap between the root interception domain and the current soil nutrient concentration distribution is compared to generate a nutrient deficit spatial distribution map. This map is mapped to the fertilization execution mechanism, using the center of the root interception domain as the target coordinate to calculate the fertilizer dispensing time delay and pulse width. This ensures that the fertilizer application location matches the actual spatial development morphology of the potato underground root system, and fertilizer particles are transported to the target coordinates for future root development. This overcomes the defect of fertilizer deviating from the effective absorption zone of the root system and improves the fertilizer interception rate in the rhizosphere microdomain.

[0060] 2. A root growth spatial resistance field is constructed by introducing soil bulk density and soil porosity parameters. The optimal descent path of the root tip growth potential energy is calculated within the resistance field as the growth trajectory. Based on porosity and moisture parameters, the root branch density and nutrient absorption interception radius are dynamically corrected to ensure that the predicted root spatial distribution conforms to the soil physical resistance environment. Historical fertilization data and soil moisture parameters are input into the soil solute transport convection-dispersion equation to calculate the time-varying nutrient concentration gradient. The difference is extracted to generate a deficit map, ensuring that the fertilization amount matches the actual free nutrient deficit state in the rhizosphere. The gravity settling time of fertilizer particles is calculated based on the trencher's penetration depth and the target fertilization depth. Combined with the current walking speed, the displacement time of the mechanism is calculated to generate the fertilizer discharge time delay. The fertilizer discharge pulse width is corrected by combining the variance of the fertilizer discharge flow rate fluctuation to eliminate the spatial positioning deviation caused by the mechanical movement and fertilizer settling physical process, ensuring that fertilizer particles fall to the target coordinates. Attached Figure Description

[0061] Figure 1 This is the overall execution flowchart of the automatic fertilization control program for potatoes of the present invention;

[0062] Figure 2 This is a flowchart illustrating the derivation of root growth trajectory and nutrient absorption interception radius in the coupled digital twin model of the present invention.

[0063] Figure 3 This is a flowchart illustrating the process of generating a spatial distribution map of soil nutrient deficiency according to the present invention.

[0064] Figure 4 This is a flowchart of the fertilizer application execution mechanism for calculating the fertilizer discharge time delay according to the present invention;

[0065] Figure 5 This is a flowchart illustrating the calculation of the fertilizer discharge pulse width of the fertilizer application actuator of the present invention;

[0066] Figure 6 This is a flowchart of the targeted fertilization closed-loop feedback correction process of the present invention. Detailed Implementation

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

[0068] Please refer to Figure 1This embodiment provides an automatic fertilization control program for potatoes based on digital twins, acquiring physiological characteristic parameters of potato growth stages, soil moisture parameters, and historical fertilization data. The physiological characteristic parameters of potato growth stages include plant height, stem diameter, leaf area index, taproot length, number of lateral roots, and growth stage identifiers. Soil moisture parameters are collected using soil moisture sensors buried at different depths in the three-dimensional soil space, ranging from 0-60 cm, with a collection interval of 24 hours. Historical fertilization data includes the fertilization time, location, fertilizer type, amount, and depth for the past three growth cycles. The physiological characteristic parameters of potato growth stages, soil moisture parameters, and historical fertilization data are input into a pre-constructed coupled digital twin model of potato root growth morphology and three-dimensional nutrient transport in the soil. The coupled digital twin model is formed by coupling the root growth morphology sub-model and the soil three-dimensional nutrient transport sub-model through a data interface. The root spatial distribution parameters output by the root growth morphology sub-model serve as the boundary conditions of the soil three-dimensional nutrient transport sub-model, and the nutrient concentration distribution parameters output by the soil three-dimensional nutrient transport sub-model serve as the growth driving force parameters of the root growth morphology sub-model.

[0069] The driven coupled digital twin model, based on potato growth period physiological parameters, soil moisture parameters, and historical fertilization data, forward extrapolates the growth trajectory and nutrient absorption interception radius of potato roots in three-dimensional soil space for the next growth cycle. The next growth cycle spans 7 days, consistent with the typical renewal cycle of potato root growth. The forward extrapolation process employs a time-stepping method with a time step of 1 hour. Within each time step, the state parameters of the root growth morphology sub-model are updated first, followed by the state parameters of the three-dimensional soil nutrient transport sub-model. Finally, data interaction between the two sub-models is completed through a coupling interface. The root growth trajectory is represented by a three-dimensional spatial curve, with each node containing spatial coordinates, growth direction, growth rate, and growth time information. The nutrient absorption interception radius is represented by the radial distance centered on the main root axis, and varies at different depths.

[0070] By comparing the spatial overlap between the root interception domain corresponding to the nutrient absorption interception radius in the coupled digital twin model and the current soil nutrient concentration distribution, a nutrient deficit spatial distribution map is generated. The root interception domain is a cylindrical spatial region with the root growth trajectory as the axis and the corresponding nutrient absorption interception radius as the radius. The current soil nutrient concentration distribution is calculated using a three-dimensional soil nutrient transport sub-model and stored in the form of a three-dimensional voxel grid, with each voxel measuring 1cm × 1cm × 1cm. The spatial overlap is calculated using the voxel intersection method, which counts the ratio of the number of voxels covered by the root interception domain to the number of voxels whose nutrient concentration reaches the potato target fertilizer requirement threshold. The nutrient deficit spatial distribution map is represented in the form of a three-dimensional grid, with each grid node storing the nutrient deficit value at that location.

[0071] The spatial distribution map of nutrient deficit is mapped to the walking path and fertilizer discharge opening of the fertilization actuator. Using the center of the root interception domain as the target coordinate, the fertilizer discharge time delay and fertilizer discharge pulse width of the actuator are calculated to generate targeted fertilization control commands. The walking path of the fertilization actuator is pre-planned as a straight line parallel to the potato planting rows, with the distance between two adjacent walking paths matching the potato planting row spacing. The fertilizer discharge opening is achieved by adjusting the opening angle of the fertilizer discharge valve, and there is a linear correspondence between the fertilizer discharge opening and the fertilizer flow rate. The target coordinate is the projection point of the geometric center of the root interception domain onto the horizontal plane. The fertilizer discharge time delay is the time required from the current moment for the fertilization actuator to reach the horizontal projection position of the target coordinate and complete the fertilizer discharge action. The fertilizer discharge pulse width is the duration for which the fertilizer discharge valve remains open. The targeted fertilization control commands include commands for adjusting the walking speed of the fertilization actuator, the fertilizer discharge valve opening time, and the fertilizer discharge valve opening duration.

[0072] In this embodiment, the construction process of the coupled digital twin model is as follows. First, potato root growth sample data from different growth stages and soil environments are collected. The sample data includes the three-dimensional spatial coordinates of the roots, root diameter, root length, number of branches, and corresponding soil physical and nutrient parameters. Then, a root growth morphology sub-model is trained based on the collected sample data. The root growth morphology sub-model adopts a parameterized root growth model based on the L-system, and the model parameters include the taproot growth rate, lateral root germination angle, lateral root growth rate, branch interval, and geotropism coefficient. Next, a three-dimensional soil nutrient transport sub-model is constructed. The soil three-dimensional nutrient transport sub-model is established based on the soil solute transport theory, considering the convection, dispersion, adsorption and desorption of nutrients, and root absorption processes. Finally, the root growth morphology sub-model and the soil three-dimensional nutrient transport sub-model are coupled, and the coupling interface enables real-time bidirectional transmission of root spatial distribution parameters and soil nutrient concentration distribution parameters.

[0073] Specifically, the process of forward extrapolating the growth trajectory and nutrient absorption interception radius of potato roots in three-dimensional soil space for the next growth cycle is as follows: First, based on the current physiological characteristics of the potato growth stage, the initial state parameters of the root growth morphology sub-model are determined. These initial state parameters include the taproot length, the number of lateral roots, the growth rate of each root segment, and the growth direction. Then, based on the current soil moisture parameters, the influence coefficient of soil moisture on root growth is calculated using the following formula:

[0074]

[0075] in, This represents the coefficient of influence of soil moisture on root growth. This represents the current soil volumetric water content. Soil wilting coefficient, It refers to the soil field water holding capacity.

[0076] Next, the influence coefficient of soil moisture on root growth is multiplied by the basic root growth rate to obtain the corrected root growth rate. Then, based on the corrected root growth rate, the growth length and spatial coordinates of each root segment at the next time step are calculated. Finally, based on the spatial coordinates of each root segment and its corresponding nutrient absorption capacity, the nutrient absorption interception radius at each location is calculated. The nutrient absorption interception radius is calculated using the following formula:

[0077]

[0078] in, For nutrient absorption interception radius, This is the root nutrient absorption coefficient. This refers to the root length density per unit volume of soil.

[0079] The process of generating a nutrient deficit spatial distribution map by comparing the spatial overlap between the root interception domain corresponding to the nutrient absorption interception radius in the coupled digital twin model and the current soil nutrient concentration distribution is as follows: First, the three-dimensional soil space is divided into a uniform voxel grid, with each voxel measuring 1cm × 1cm × 1cm. Then, based on the root growth trajectory and nutrient absorption interception radius, it is determined whether each voxel belongs to the root interception domain. For voxels belonging to the root interception domain, the current soil nutrient concentration value within that voxel is extracted. Next, the extracted soil nutrient concentration value is compared with the potato target fertilizer requirement threshold to calculate the nutrient deficit for each voxel. The nutrient deficit is calculated using the following formula:

[0080]

[0081] in, coordinates Nutritional deficiency of the body's own body. The target fertilizer concentration threshold for potatoes. coordinates The current soil nutrient concentration value of the local element.

[0082] Finally, the nutrient deficit of all voxels is fitted using three-dimensional grid interpolation to generate a spatial distribution map of the nutrient deficit. The three-dimensional grid interpolation uses inverse distance weighted interpolation, and the interpolation formula is as follows:

[0083]

[0084] in, coordinates Nutrient deficiency at the site, For the first Nutrient deficit at a known point coordinates To the Distance between known points The power exponent, The number of known points involved in the interpolation.

[0085] The process of mapping the spatial distribution map of nutrient deficit to the walking path and fertilization opening of the fertilization actuator, using the root interception domain center as the target coordinate, and calculating the fertilization time delay and fertilization pulse width of the fertilization actuator to generate targeted fertilization control commands is as follows: First, the spatial distribution map of nutrient deficit is projected onto a horizontal plane to obtain the horizontal nutrient deficit distribution. Then, based on the pre-planned walking path of the fertilization actuator, the nutrient deficit amount corresponding to each position on the walking path is determined. Next, based on the nutrient deficit amount at each position on the walking path, the corresponding fertilization opening is determined. There is a linear correspondence between the fertilization opening and the nutrient deficit amount, which is obtained through pre-calibration. Then, the projection point of the center of each root interception domain on the horizontal plane is extracted as the target coordinate. Next, the time required for the fertilization actuator to reach the horizontal projection position of each target coordinate from the current position is calculated as the fertilization time delay. The fertilization time delay is calculated using the following formula:

[0086]

[0087] in, To delay the excretion of fertilizer, This represents the distance from the current location of the fertilization actuator to the horizontal projection of the target coordinates. The walking speed of the fertilization implementation agency.

[0088] Then, based on the nutrient deficit at each target coordinate location, the corresponding nutrient excretion pulse width is calculated. The nutrient excretion pulse width is calculated using the following formula:

[0089]

[0090] in, The width of the fertilizer excretion pulse. The amount of fertilizer required for the target coordinate location. This refers to the amount of fertilizer discharged per unit time by the fertilizer discharge valve.

[0091] Finally, the fertilization time delay, fertilization pulse width, and corresponding target coordinate information were integrated to generate targeted fertilization control instructions. Table 1 shows the comparison of potato root growth parameters and nutrient absorption interception radius at different growth stages.

[0092] Table 1 Comparison of potato root growth parameters and nutrient absorption interception radius at different growth stages

[0093] Seedling stage 1.2 45 0.8 0.3 1.5 tuber formation period 2.5 60 1.5 0.8 2.8 tuber enlargement period 1.8 75 1.2 1.2 3.5 Maturity 0.5 90 0.3 0.6 2.2

[0094] Table 1 shows the correlation between potato root growth parameters and nutrient absorption interception radius at different growth stages. In practical applications, based on the potato's growth stage, the corresponding root growth parameters and nutrient absorption interception radius are retrieved from the table and used as input parameters for the coupled digital twin model. Table 1 allows for the rapid determination of potato root growth characteristics and nutrient absorption capacity at different growth stages, improving the inference efficiency of the coupled digital twin model.

[0095] In this embodiment, the fertilization actuator includes a chassis, a furrow opener, a fertilizer discharge box, a fertilizer discharge valve, and a control system. The chassis is a tracked chassis, capable of adapting to different terrain conditions in the field. The furrow opener is installed at the front of the chassis and is used to create fertilization furrows in the soil. The fertilizer discharge box is installed in the middle of the chassis and is used to store fertilizer. The fertilizer discharge valve is installed at the bottom of the fertilizer discharge box and is used to control the amount of fertilizer discharged. The control system is installed in the driver's cab of the chassis and is used to receive targeted fertilization control commands and, based on these commands, control the chassis's travel speed, the furrow opener's penetration depth, and the opening and closing of the fertilizer discharge valve.

[0096] In this embodiment, the physiological characteristics of potatoes during their growth period, soil moisture parameters, and historical fertilization data are first collected using a field sensor network. Then, the collected data is transmitted to a remote server, which hosts a coupled digital twin model. The remote server drives the coupled digital twin model to perform forward inference, generating a spatial distribution map of nutrient deficit. Next, the remote server calculates the fertilization time delay and fertilization pulse width based on the nutrient deficit spatial distribution map, generating targeted fertilization control commands. Finally, the targeted fertilization control commands are transmitted to the control system of the fertilization actuator via a wireless communication network. The control system then controls the fertilization actuator to complete the targeted fertilization operation according to the commands.

[0097] In a preferred embodiment, reference Figure 2 This study acquired physiological parameters of potatoes during their growth period, soil moisture parameters, and historical fertilization data. It also included obtaining soil bulk density and porosity parameters in a three-dimensional soil space. Soil bulk density and porosity parameters were obtained through soil sampling analysis, with a sampling depth range of 0-60 cm and a sampling interval of 10 cm. Three replicate samples were collected from each sampling point, and the average value was used as the soil bulk density and porosity parameters for that sampling point. These parameters were then incorporated into a coupled digital twin model to construct a root growth spatial resistance field. This root growth spatial resistance field is a three-dimensional scalar field, where the value at each point represents the magnitude of soil resistance to root growth at that location.

[0098] Based on the root tropisms determined by physiological parameters during the potato growth period, the optimal descent path of the root apex growth potential energy is determined within the root growth space resistance field. This optimal descent path is then defined as the growth trajectory. Root tropisms include geotropism, hydrotropism, and nutrient tropism. Geotropism causes roots to grow in the direction of gravity, hydrotropism causes roots to grow towards areas with high soil moisture content, and nutrient tropism causes roots to grow towards areas with high soil nutrient concentration. The root apex growth potential energy is the potential energy possessed by the root apex in the growth space resistance field, and its magnitude is directly proportional to the magnitude of the soil resistance at that location. The optimal descent path is the path from the current position where the potential energy decreases most rapidly, while satisfying the constraints of the tropisms.

[0099] Based on soil porosity and soil moisture parameters in the region where the growth trajectory is located, the root hairy root branching density is dynamically adjusted, and the nutrient absorption interception radius is calculated based on the hairy root branching density. Hairy root branching density is the number of hairy roots growing per unit length of taproot. Soil porosity parameters affect the growth space of hairy roots, while soil moisture parameters affect the growth rate of hairy roots. When soil porosity is high and soil moisture content is suitable, the hairy root branching density is high; when soil porosity is low or soil moisture content is too low or too high, the hairy root branching density is low. The nutrient absorption interception radius is proportional to the square root of the hairy root branching density.

[0100] Specifically, the process of introducing soil bulk density and soil porosity parameters into a coupled digital twin model to construct the root growth space resistance field is as follows: First, the soil penetration resistance distribution value is calculated based on the soil bulk density and soil porosity parameters. The soil penetration resistance distribution value is calculated using the following formula:

[0101]

[0102] in, This represents the distribution value of soil penetration resistance. For soil bulk density, For soil porosity, , , This is an empirical coefficient, determined through soil penetration tests.

[0103] Then, the soil penetration resistance distribution value is compared with the root critical penetration stress to delineate the root-impenetrable boundary. The root critical penetration stress is the maximum resistance value that the root system can withstand when penetrating the soil, and the critical penetration stress varies at different growth stages of potato roots. When the soil penetration resistance distribution value at a certain location is greater than the root critical penetration stress, that location is designated as the root-impenetrable boundary; when the soil penetration resistance distribution value at a certain location is less than or equal to the root critical penetration stress, that location is designated as the root-penetrable region. Finally, the soil penetration resistance distribution value is used as the value of a three-dimensional scalar field to construct the root growth space resistance field.

[0104] Based on the root tropism growth trend determined by the physiological characteristics parameters of potato growth period, the optimal descent path of the root apical growth potential energy is solved in the root growth space resistance field, and the process of determining the optimal descent path as the growth trajectory is as follows: First, based on the physiological characteristics parameters of potato growth period, the root geotropism coefficient, hydrotropism coefficient, and fertility coefficient are determined. The values ​​of the geotropism coefficient, hydrotropism coefficient, and fertility coefficient range from 0 to 1, with larger coefficients indicating a stronger corresponding tropism growth trend. Then, the root apical growth potential energy function is constructed, which includes soil resistance potential energy terms, geotropism potential energy terms, hydrotropism potential energy terms, and fertility potential energy terms. The growth potential energy function is calculated using the following formula:

[0105]

[0106] in, coordinates The apical growth potential of the root system at that location. The soil resistance potential energy coefficient, coordinates Soil penetration resistance distribution at the location, The geotropic potential coefficient, Let be the geotropic potential function. The hydrotropic potential coefficient, Let be the hydrotropic potential function. The coefficient of fertilization potential energy. Let be the tropism potential function.

[0107] Geotropic potential function With depth The geotropic potential energy is directly proportional to the depth; the greater the depth, the smaller the geotropic potential energy. Hydrotropic potential energy function. The hydrotropism potential energy is inversely proportional to soil moisture content; the higher the soil moisture content, the lower the hydrotropism potential energy. It is inversely proportional to the soil nutrient concentration; the higher the soil nutrient concentration, the lower the nutrient tropism potential.

[0108] Next, the rapid descent method is used to solve for the optimal descent path of the root tip growth potential. The rapid descent method is a numerical method for solving the static Hamilton-Jacobi equation, which can efficiently calculate the shortest path from the starting point to the ending point. In this embodiment, the starting point is the current position of the root tip, and the ending point is the expected position of the root tip at the end of the next growth cycle. The optimal descent path is the path from the starting point along the direction of the fastest decrease in growth potential to the ending point. Finally, the optimal descent path obtained is determined as the growth trajectory of the potato root system in the next growth cycle.

[0109] Based on the soil porosity and soil moisture parameters of the growth trajectory area, the root hairline branch density is dynamically adjusted, and the nutrient absorption interception radius is calculated based on the hairline branch density as follows: First, the path length of the growth trajectory through the pores between the impenetrable boundaries of adjacent roots is extracted. Then, the seepage pressure on the main root surface is calculated based on the path length and soil moisture parameters. The seepage pressure on the main root surface is calculated using the following formula:

[0110]

[0111] in, The seepage pressure on the main root surface, The density of water, It is the acceleration due to gravity. For soil water potential, This is the path length of the growth trajectory through the pores between the impenetrable boundaries of adjacent roots. This is the reference path length.

[0112] When the seepage pressure on the main root surface exceeds the osmotic pressure threshold of the root cortex cells, it triggers an exponential increase in the density of hairy root branches along the seepage pressure gradient on the main root surface. The correction formula for hairy root branch density is as follows:

[0113]

[0114] in, The corrected hair root branch density, Based on the density of basal hair root branches, The coefficient representing the influence of seepage pressure is denoted as . This represents the osmotic pressure threshold of the root cortical cells.

[0115] Finally, the radial radius of the envelope surface of the increased root branch density is set as the nutrient absorption interception radius. The nutrient absorption interception radius is calculated using the following formula:

[0116]

[0117] in, For nutrient absorption interception radius, For the density of hair root branches, The average length of a single hairy root is given. Table 2 shows the distribution of soil penetration resistance and the comparison with the critical penetration stress of the root system under different soil bulk density and porosity conditions.

[0118] Table 2. Comparison of soil penetration resistance distribution and root critical penetration stress under different soil bulk density and porosity conditions.

[0119] 1.1 58 120 150 200 1.2 55 180 150 200 1.3 51 250 150 200 1.4 47 330 150 200 1.5 43 420 150 200

[0120] Table 2 shows the correlation between the distribution values ​​of soil penetration resistance under different soil bulk density and porosity conditions and the critical penetration stress of potato roots at different growth stages. In practical applications, based on soil bulk density and porosity parameters, the corresponding soil penetration resistance distribution value can be found in the table and compared with the critical penetration stress of the root system at the corresponding growth stage to delineate the non-penetrating boundary of the root system. Table 2 allows for the rapid determination of the root's scalable area under different soil physical conditions, improving the accuracy of root growth trajectory prediction.

[0121] In this embodiment, the resolution of the root growth spatial resistance field is 1cm×1cm×1cm, consistent with the voxel grid resolution of the soil three-dimensional nutrient transport sub-model. When constructing the root growth spatial resistance field, Kriging interpolation is used to interpolate discrete soil sampling point data into a continuous three-dimensional scalar field. Kriging interpolation is a spatial interpolation method based on variograms, which can fully consider the correlation of spatial data and improve the accuracy of the interpolation results.

[0122] In this embodiment, soil bulk density and porosity parameters within a three-dimensional soil space are first obtained through soil sampling analysis. Then, these parameters are input into a coupled digital twin model to construct a root growth spatial resistance field. Next, based on the root tropism determined by the physiological characteristics of potato growth stages, the optimal descent path of the root apical growth potential is calculated within the root growth spatial resistance field, yielding the root trajectory for the next growth cycle. Then, based on the soil porosity and moisture parameters of the region containing the growth trajectory, the root hair branch density is dynamically adjusted, and the nutrient absorption interception radius is calculated. Finally, based on the root growth trajectory and the nutrient absorption interception radius, the root interception domain is determined, providing a basis for subsequently generating a nutrient deficit spatial distribution map.

[0123] In a preferred embodiment, reference Figure 3A coupled digital twin model embeds a soil solute transport convection-dispersion equation. This equation describes the transport patterns of solutes in soil under the influence of convection, dispersion, adsorption-desorption, and root uptake. Historical fertilization data and soil moisture parameters are input into the equation to calculate the time-varying nutrient concentration gradient of each spatial voxel in the three-dimensional soil space at the current moment. The time-varying nutrient concentration gradient represents the rate of change of nutrient concentration over time, reflecting the dynamic trend of nutrient concentration in the soil.

[0124] Extract the target voxel set covered by the root interception domain, and calculate the spatial distribution of the difference between the integral value of the time-varying nutrient concentration gradient within the target voxel set and the potato target fertilizer requirement threshold. The target voxel set is the set of all voxels covered by the root interception domain. The integral value of the time-varying nutrient concentration gradient is the total change in nutrient concentration within the target voxel over a certain period of time. The spatial distribution of the difference is the spatial distribution of the difference between the potato target fertilizer requirement threshold and the integral value of the time-varying nutrient concentration gradient within each target voxel.

[0125] Voxel regions with differences greater than zero in the difference spatial distribution are extracted as nutrient-deficient spaces. Three-dimensional gridded interpolation is then performed on these nutrient-deficient spaces to generate a nutrient-deficient spatial distribution map. A difference greater than zero indicates that the nutrient concentration within that voxel will fall below the target fertilizer requirement threshold for potatoes in the near future, necessitating fertilization. The three-dimensional gridded interpolation using Kriging interpolation improves the resolution and accuracy of the nutrient-deficient spatial distribution map.

[0126] Specifically, the process of inputting historical fertilization data and soil moisture parameters into the soil solute transport convection-dispersion equation to calculate the time-varying nutrient concentration gradient of each spatial voxel in the three-dimensional soil space at the current moment is as follows: First, the isothermal adsorption and desorption parameters of soil particles for nutrient ions are obtained, and these parameters are introduced as source and sink terms into the soil solute transport convection-dispersion equation. The form of the soil solute transport convection-dispersion equation is as follows:

[0127]

[0128] in, This refers to the soil volumetric water content. This represents the concentration of the solute in the soil solution. For time, The hydrodynamic dispersion coefficient, Soil moisture flux. It is a source-sink item, including the adsorption and desorption processes of nutrients and the root absorption process.

[0129] The isothermal adsorption and desorption process is described by the Langmuir isothermal adsorption equation, which is as follows:

[0130]

[0131] in, This represents the amount of nutrients adsorbed per unit mass of soil particles. The Langmuir adsorption constant is given. This represents the maximum adsorption capacity of soil particles.

[0132] The root absorption process is described by the Michaelis-Menten equation, which has the following form:

[0133]

[0134] in, This refers to the amount of nutrients absorbed by the roots per unit volume of soil. This represents the maximum absorption rate of the root system. It is the Michaelis constant. This refers to the root length density per unit volume of soil.

[0135] Source and Exchange The combined effect of adsorption / desorption and root absorption processes is expressed as follows:

[0136]

[0137] in, It is the bulk density of the soil.

[0138] Then, combining the water infiltration flux in the soil moisture parameters, the soil solute transport convection-dispersion equation after introducing source and sink terms is solved to obtain the time-varying nutrient concentration gradient of free nutrient ions in each spatial voxel. The solution process uses the finite difference method, dividing the three-dimensional soil space into a uniform voxel grid and the time into uniform time steps, updating the solute concentration of each voxel within each time step. The discretization form of the finite difference method is as follows:

[0139]

[0140] Among them, superscript Indicates the first Each time step, index Coordinates representing voxels For time step, , , Voxels in , , Dimensions in the direction.

[0141] Finally, the time-varying nutrient concentration gradient was calculated based on the change in the concentration of free nutrient ions in each spatial voxel over time. The time-varying nutrient concentration gradient was calculated using the following formula:

[0142]

[0143] in, coordinates The time-varying concentration gradient of nutrients in the voxel.

[0144] The process of extracting the target voxel set covered by the root interception domain and calculating the spatial distribution of the difference between the integral value of the time-varying nutrient concentration gradient within the target voxel set and the target nutrient requirement threshold for potatoes is as follows: First, based on the root growth trajectory and nutrient absorption interception radius, it is determined whether each voxel belongs to the root interception domain, and voxels belonging to the root interception domain are extracted to form the target voxel set. Then, within the target voxel set, the volume integral of the time-varying nutrient concentration gradient of free nutrient ions is performed along the normal direction of the nutrient absorption interception radius to obtain the effective nutrient integral that the root system can actually absorb. The formula for calculating the volume integral is as follows:

[0145]

[0146] in, This represents the actual amount of effective nutrients that the root system can absorb. The volume of the target voxel.

[0147] Next, the spatial difference between the effective nutrient integral and the target fertilizer concentration threshold for potatoes is calculated. The spatial difference is calculated using the following formula:

[0148]

[0149] in, coordinates Spatial difference of voxels, To meet the target fertilizer requirements for potatoes coordinates The effective nutrient content of the body element.

[0150] Finally, the spatial difference is mapped back to the 3D mesh to generate the spatial distribution of the difference.

[0151] The process of extracting voxel regions with differences greater than zero from the difference spatial distribution as nutrient-deficient spaces, and then performing 3D gridded interpolation fitting on these nutrient-deficient spaces to generate a nutrient-deficient spatial distribution map is as follows: First, all voxels in the difference spatial distribution are traversed, and voxels with differences greater than zero are extracted to form nutrient-deficient spaces. Then, 3D gridded interpolation fitting is performed on these nutrient-deficient spaces to transform the discrete voxel data into a continuous 3D spatial distribution. The 3D gridded interpolation fitting uses the Kriging interpolation method, and the interpolation formula is as follows:

[0152]

[0153] in, coordinates Nutrient deficiency at the site, These are the Kriging weighting coefficients. For the first Nutrient deficit at a known point The number of known points involved in the interpolation.

[0154] The Kriging weighting coefficients are obtained by solving the following system of equations:

[0155]

[0156] in, For point With point The variogram values ​​between For point with interpolation point The variogram values ​​between It is a Lagrange multiplier.

[0157] Finally, the interpolated continuous three-dimensional spatial distribution is used as a spatial distribution map of nutrient deficiency.

[0158] In this embodiment, reference Figure 4 The nutrient deficit spatial distribution map is mapped onto the walking path and fertilizer dispensing opening of the fertilization actuator. Using the center of the root interception domain as the target coordinate, the fertilizer dispensing time delay of the actuator is calculated. This includes acquiring the current walking speed of the actuator and the trencher's penetration depth. Combined with the nutrient deficit spatial distribution map, the relative offset of the target coordinate in the horizontal direction and the target fertilization depth in the vertical direction are determined. The trencher's penetration depth is acquired in real time by a depth sensor installed on the trencher. The relative offset of the target coordinate in the horizontal direction is the vertical distance between the horizontal projection position of the target coordinate and the walking path of the fertilization actuator. The target fertilization depth in the vertical direction is the vertical coordinate of the center of the root interception domain.

[0159] The gravity settling time required for fertilizer granules to fall is calculated based on the depth of the furrow opener and the target fertilization depth. After being discharged from the discharge port, fertilizer granules first fall into the air to the bottom of the furrow opened by the furrow opener, and then continue to fall through the soil gaps to the target fertilization depth. The gravity settling time is the time required for fertilizer granules to reach the target fertilization depth from the discharge port.

[0160] The displacement time required for the fertilization actuator to reach the horizontal projection position of the target coordinates is calculated based on the current travel speed and relative offset. The displacement time is the time required for the fertilization actuator to travel from its current position to the horizontal projection position of the target coordinates.

[0161] The sum of the gravity settling time and the mechanism displacement time is set as the fertilizer discharge time delay, so that the fertilizer application actuator performs the fertilizer discharge action according to the fertilizer discharge time delay when it reaches the target coordinate. By setting the fertilizer discharge time delay, the time difference between the fertilizer particle falling process and the fertilizer application actuator traveling process can be compensated, ensuring that the fertilizer particles fall accurately at the target coordinate position.

[0162] Specifically, the process of calculating the gravity settling time required for fertilizer particles to fall, based on the trencher's penetration depth and the target fertilization depth, is as follows: First, the particle size distribution and particle density parameters of the fertilizer particles are obtained. Then, a particle settling resistance model is constructed by combining this with soil moisture parameters between the trencher's penetration depth and the target fertilization depth. The particle settling resistance model considers both the frictional resistance of soil particles to fertilizer particles and the viscous resistance of soil moisture to fertilizer particles. The particle settling resistance is calculated using the following formula:

[0163]

[0164] in, For particle settling resistance, The drag coefficient, The density of the soil medium, This represents the cross-sectional area of ​​the fertilizer granules. This refers to the settling velocity of the fertilizer particles.

[0165] Then, the particle size distribution parameters, particle density parameters, and soil moisture parameters are input into the particle settling resistance model to calculate the terminal settling velocity of fertilizer particles in the soil interstitial spaces. The terminal settling velocity is the falling velocity of fertilizer particles when gravity and settling resistance are in equilibrium. The terminal settling velocity is calculated using the following formula:

[0166]

[0167] in, For terminal settlement velocity, It is the acceleration due to gravity. The diameter of the fertilizer granules. This refers to the density of the fertilizer particles.

[0168] Finally, the gravity settlement time is calculated based on the difference between the terminal settlement velocity and the trencher's penetration depth and the target fertilization depth. The gravity settlement time is calculated using the following formula:

[0169]

[0170] in, For gravity settling time, This refers to the depth to which the trencher penetrates the soil. Target fertilization depth.

[0171] The process of calculating the displacement time required for the fertilization actuator to reach the horizontal projection position of the target coordinates based on the current walking speed and relative offset is as follows: First, obtain the chassis vibration frequency and the elevation difference of the ground surface undulations of the fertilization machinery. Then, correct the current walking speed based on the chassis vibration frequency and the elevation difference of the ground surface undulations to obtain the actual effective walking speed. The actual effective walking speed is calculated using the following formula:

[0172]

[0173] in, The actual effective walking speed, The current walking speed, This is the vibration influence coefficient. The chassis vibration frequency, This refers to the elevation difference caused by surface undulations.

[0174] Then, the mechanism displacement time is calculated based on the actual effective travel speed and relative offset. The mechanism displacement time is calculated using the following formula:

[0175]

[0176] in, For the displacement time of the mechanism, This represents the relative offset of the target coordinates in the horizontal direction.

[0177] Finally, the sum of the gravity settling time and the mechanism displacement time is set as the fertilizer discharge time delay. The fertilizer discharge time delay is calculated using the following formula:

[0178]

[0179] In this embodiment, reference Figure 5The calculation of the fertilizer discharge pulse width of the fertilization actuator involves extracting the deficit volume and deficit concentration difference of the deficit grid corresponding to the target coordinate in the nutrient deficit spatial distribution map, and combining this with a preset single-pulse fertilizer discharge calibration value to calculate the initial number of fertilizer discharge pulses. The deficit volume is the volume of the deficit grid corresponding to the target coordinate. The deficit concentration difference is the ratio of the nutrient deficit amount of the deficit grid corresponding to the target coordinate to the soil volume. The single-pulse fertilizer discharge calibration value is the amount of fertilizer discharged when the fertilizer discharge valve opens for one reference pulse width, which is obtained through pre-calibration.

[0180] The variance of the fertilizer discharge velocity fluctuation at the fertilizer discharge outlet of the fertilization actuator is obtained. This variance is then input into a pre-constructed pulse width modulation (PWM) compensation model, which outputs a pulse width correction coefficient for the initial number of fertilizer discharge pulses. The variance of the fertilizer discharge velocity fluctuation reflects the instability of the discharge velocity; the larger the variance, the more unstable the discharge velocity. The PWM compensation model is used to compensate for the error in fertilizer discharge volume caused by the fluctuation of the discharge velocity.

[0181] The target fertilization pulse width for the target coordinates is obtained by multiplying the initial number of fertilization pulses and the baseline pulse width corresponding to the variance of fertilization flow rate fluctuation by the pulse width correction coefficient. The baseline pulse width is the standard duration when the fertilization valve opens with one pulse.

[0182] The fertilization actuator controls the opening and closing duration of the fertilizer discharge valve based on the target fertilizer discharge pulse width, delivering the corresponding amount of fertilizer to the target coordinate position. By correcting the fertilizer discharge pulse width, errors in fertilizer discharge caused by fluctuations in fertilizer flow rate can be compensated, thus improving the accuracy of fertilizer application.

[0183] Specifically, the process of extracting the difference between the deficit volume and deficit concentration of the deficit grid corresponding to the target coordinate in the nutrient deficit spatial distribution map, and calculating the initial number of fertilizer discharge pulses based on the preset single-pulse fertilizer discharge calibration value, is as follows: First, determine the deficit grid corresponding to the target coordinate and extract the difference between the deficit volume and deficit concentration of that grid. Then, calculate the total fertilizer application required for that deficit grid. The total fertilizer application is calculated using the following formula:

[0184]

[0185] in, This is the total amount of fertilizer applied. For the missing volume, This represents the deficit concentration difference.

[0186] Next, the total fertilizer application rate is divided by the calibrated value of the single-pulse fertilizer discharge rate to obtain the initial number of fertilizer discharge pulses. The initial number of fertilizer discharge pulses is calculated using the following formula:

[0187]

[0188] in, This represents the initial number of fertilizer excretion pulses. This is the calibration value for single-pulse fertilizer excretion.

[0189] The process of obtaining the variance of fertilizer flow velocity fluctuation at the fertilizer discharge port of the fertilization actuator, inputting this variance into a pre-constructed pulse width modulation compensation model, and outputting the pulse width correction coefficient for the initial number of fertilizer discharge pulses is as follows: First, the fluctuation values ​​of the fertilizer discharge shaft speed and the change values ​​of the fertilizer discharge port opening gap are collected in real time. The variance of fertilizer flow velocity fluctuation is calculated based on these values. The variance of fertilizer flow velocity fluctuation is calculated using the following formula:

[0190]

[0191] in, For the variance of fertilizer discharge velocity fluctuation, The variance of the discharge shaft speed fluctuation. The variance of the variation in the opening gap of the fertilizer discharge outlet. , These are the weighting coefficients.

[0192] Then, the variance of the fertilizer discharge velocity fluctuation is matched with a preset variance-coefficient mapping table. This table is constructed based on the nonlinear velocity compensation coefficients corresponding to different variance intervals, and the pulse width correction coefficient is obtained from the query. The variance-coefficient mapping table is established through a large amount of fertilizer discharge experimental data, and it stores the pulse width correction coefficients corresponding to different fertilizer discharge velocity fluctuation variance intervals.

[0193] The target fertilization pulse width is obtained by multiplying the initial number of fertilization pulses and the variance of fertilization flow rate fluctuations by the baseline pulse width correction coefficient. The target fertilization pulse width is calculated using the following formula:

[0194]

[0195] in, The target nutrient excretion pulse width, The reference pulse width, This is the pulse width correction factor.

[0196] The fertilization actuator controls the opening and closing duration of the fertilizer discharge valve based on the target fertilizer discharge pulse width to deliver the corresponding amount of fertilizer to the target coordinate position as follows: First, the target fertilizer discharge pulse width is converted into a control level signal for the fertilizer discharge valve. When the control level signal is high, the fertilizer discharge valve opens; when the control level signal is low, the fertilizer discharge valve closes. Then, when the fertilizer discharge time delay arrives, the output control level signal drives the opening and closing duration of the fertilizer discharge valve, allowing the fertilizer particles to be delivered to the target coordinate position under the combined action of gravity and forced airflow. The forced airflow is provided by a fan installed at the bottom of the fertilizer discharge box to improve the conveying speed and uniformity of the fertilizer particles. The pulse width correction coefficients corresponding to different fertilizer discharge flow rate fluctuation variance ranges are shown in Table 3.

[0197] Table 3. Comparison of Pulse Width Correction Coefficients for Different Variance Ranges of Fertilizer Discharge Velocity Fluctuations

[0198] Variance range of fertilizer discharge velocity fluctuation ((m / s)²) Pulse width correction factor 0-0.01 1.00 0.01-0.03 1.05 0.03-0.05 1.12 0.05-0.07 1.20 0.07-0.10 1.30

[0199] Table 3 shows the correlation between pulse width correction coefficients corresponding to different variance intervals of fertilizer discharge flow rate fluctuations. In practical applications, based on the calculated variance of the fertilizer discharge flow rate fluctuations, the corresponding variance interval is determined, and then the corresponding pulse width correction coefficient is looked up from the table. Table 3 allows for the rapid determination of the correction amount for the fertilizer discharge pulse width, improving the accuracy of fertilizer application.

[0200] In this embodiment, reference Figure 6 After generating the targeted fertilization control command, the system also includes acquiring real-time feedback data on the actual fertilizer discharge from the fertilization execution mechanism and real-time multispectral nutrient detection data of the soil profile after fertilization. The real-time fertilizer discharge feedback data is collected in real-time by a flow sensor installed at the discharge port. The real-time multispectral nutrient detection data is collected in real-time by a multispectral sensor installed at the rear of the fertilization execution mechanism; the multispectral sensor can detect the nutrient content at different depths of the soil profile.

[0201] Real-time fertilizer application feedback data and real-time multispectral nutrient detection data are input into a coupled digital twin model, driving the model to update the current soil nutrient concentration distribution and obtain the actual soil nutrient distribution state after fertilization. The coupled digital twin model updates the fertilizer application input parameters based on the real-time fertilizer application feedback data and updates the soil nutrient concentration distribution parameters based on the real-time multispectral nutrient detection data.

[0202] The spatial concentration residual is calculated by comparing the actual distribution of soil nutrients after fertilization with the target nutrient distribution corresponding to the spatial distribution map of nutrient deficit. The spatial concentration residual is the difference between the actual soil nutrient concentration and the target nutrient concentration within each voxel after fertilization.

[0203] If the spatial concentration residual exceeds a preset residual threshold, the root growth rate parameters and soil nutrient diffusion coefficient in the coupled digital twin model are corrected based on the spatial concentration residual, and targeted fertilization control instructions are regenerated based on the corrected coupled digital twin model. The preset residual threshold is set according to the fertilization accuracy requirements. When the spatial concentration residual exceeds the preset residual threshold, it indicates that there is an error in the inference results of the coupled digital twin model, and the model parameters need to be corrected. Through closed-loop feedback correction, the inference accuracy of the coupled digital twin model and the accuracy of fertilization control can be continuously improved.

[0204] Specifically, the process of calculating the spatial concentration residual by comparing the actual distribution of soil nutrients after fertilization with the target nutrient distribution corresponding to the nutrient deficit spatial distribution map is as follows: First, calculate the target nutrient concentration for each voxel based on the nutrient deficit spatial distribution map. The target nutrient concentration is the threshold value for the potato's target nutrient requirement. Then, compare the actual concentration of each voxel in the actual distribution of soil nutrients after fertilization with the corresponding target nutrient concentration to calculate the spatial concentration residual. The spatial concentration residual is calculated using the following formula:

[0205]

[0206] in, coordinates Spatial concentration residuals of voxels coordinates The actual concentration of soil nutrients after fertilization with phytonutrients. coordinates The target concentration of nutrients in the body.

[0207] If the spatial concentration residual exceeds a preset residual threshold, the process of correcting the root growth rate parameter and soil nutrient diffusion coefficient in the coupled digital twin model based on the spatial concentration residual is as follows. First, the root mean square error of the spatial concentration residual is calculated. The root mean square error is calculated using the following formula:

[0208]

[0209] in, The root mean square error of the spatial concentration residuals. For the number of voxels, For the first Spatial concentration residuals of individual elements.

[0210] Then, if the root mean square error exceeds the preset residual threshold, the gradient descent method is used to correct the root growth rate parameter and soil nutrient diffusion coefficient in the coupled digital twin model. The correction formula for the gradient descent method is as follows:

[0211]

[0212] in, These are the corrected model parameters. These are the model parameters before correction. For learning rate, This represents the partial derivative of the root mean square error with respect to the model parameters.

[0213] Finally, based on the revised coupled digital twin model, a forward deduction is performed again to generate a new spatial distribution map of nutrient deficit, and the new fertilizer excretion time delay and fertilizer excretion pulse width are calculated to generate new targeted fertilizer application control instructions.

[0214] In this embodiment, historical fertilization data and soil moisture parameters are first input into the soil solute transport convection-dispersion equation to calculate the time-varying nutrient concentration gradient of each spatial voxel in the three-dimensional soil space at the current moment. Then, the target voxel set covered by the root interception domain is extracted, and the spatial distribution of the difference between the integral value of the time-varying nutrient concentration gradient within the target voxel set and the potato target fertilizer requirement threshold is calculated. Next, voxel regions with differences greater than zero in the spatial distribution of differences are extracted as nutrient-deficient spaces, generating a nutrient-deficient spatial distribution map. Then, the fertilizer discharge time delay and fertilizer discharge pulse width are calculated based on the nutrient-deficient spatial distribution map to generate targeted fertilization control instructions. Next, the fertilization execution mechanism completes the fertilization operation according to the targeted fertilization control instructions. Finally, real-time fertilizer discharge feedback data and real-time multispectral nutrient detection data are acquired to drive the coupled digital twin model to update the soil nutrient concentration distribution. Finally, the actual distribution of soil nutrients after fertilization is compared with the target distribution of nutrients, and the spatial concentration residual is calculated. If the residual exceeds the threshold, the model parameters are corrected and the targeted fertilization control instructions are regenerated.

Claims

1. An automatic fertilization control program for potatoes based on digital twins, characterized in that, include: Obtain physiological parameters of potato growth period, soil moisture parameters, and historical fertilization data; The potato growth period physiological characteristic parameters, the soil moisture parameters, and the historical fertilization data are input into a pre-constructed coupled digital twin model of potato root growth morphology and soil three-dimensional nutrient transport. The coupled digital twin model is driven to forward extrapolate the growth trajectory and nutrient absorption interception radius of the potato root system in three-dimensional soil space in the next growth cycle based on the potato growth period physiological characteristic parameters, the soil moisture parameters and the historical fertilization data. By comparing the spatial overlap between the root interception domain corresponding to the nutrient absorption interception radius in the coupled digital twin model and the current soil nutrient concentration distribution, a nutrient deficit spatial distribution map is generated. The nutrient deficiency spatial distribution map is mapped to the walking path and fertilization opening of the fertilization actuator. Using the center of the root interception domain as the target coordinate, the fertilization time delay and fertilization pulse width of the fertilization actuator are calculated to generate targeted fertilization control commands.

2. The automatic fertilization control program for potatoes based on digital twins according to claim 1, characterized in that, The acquisition of potato growth period physiological characteristic parameters, soil moisture parameters and historical fertilization data also includes the acquisition of soil bulk density parameters and soil porosity parameters in three-dimensional soil space. The driving force of the coupled digital twin model, based on the potato growth period physiological characteristic parameters, soil moisture parameters, and historical fertilization data, forward extrapolates the growth trajectory and nutrient absorption interception radius of the potato root system in three-dimensional soil space for the next growth cycle, including: The soil bulk density parameter and the soil porosity parameter are introduced into the coupled digital twin model to construct a root growth spatial resistance field; Based on the root tropism growth trend determined by the physiological characteristic parameters of potato growth period, the optimal descent path of the root tip growth potential energy is solved in the root growth space resistance field, and the optimal descent path is determined as the growth trajectory. Based on the soil porosity parameters and soil moisture parameters of the region where the growth trajectory is located, the root hair branch density is dynamically adjusted, and the nutrient absorption interception radius is calculated based on the root hair branch density.

3. The automatic fertilization control program for potatoes based on digital twins according to claim 1, characterized in that, The coupled digital twin model incorporates soil solute transport convection-dispersion equations; By comparing the spatial overlap between the root interception domain corresponding to the nutrient absorption interception radius in the coupled digital twin model and the current soil nutrient concentration distribution, a nutrient deficit spatial distribution map is generated, including: The historical fertilization data and the soil moisture parameters are input into the soil solute transport convection-diffusion equation to calculate the time-varying nutrient concentration gradient of each spatial voxel in the three-dimensional soil space at the current moment. Extract the target voxel set covered by the root interception domain, and calculate the spatial distribution of the difference between the integral value of the time-varying nutrient concentration gradient within the target voxel set and the potato target fertilizer concentration threshold. The voxel regions with differences greater than zero in the difference space distribution are extracted as nutrient deficiency spaces. The nutrient deficiency spaces are then subjected to three-dimensional gridded interpolation fitting to generate the nutrient deficiency space distribution map.

4. The automatic fertilization control program for potatoes based on digital twins according to claim 1, characterized in that, Mapping the nutrient deficit spatial distribution map to the walking path and fertilizer discharge opening of the fertilization execution mechanism, and using the center of the root interception domain as the target coordinate, the fertilizer discharge time delay of the fertilization execution mechanism is calculated, including: The current walking speed of the fertilization actuator and the soil penetration depth of the furrow opener are obtained, and the relative offset of the target coordinates in the horizontal direction and the target fertilization depth in the vertical direction are determined by combining the nutrient deficiency spatial distribution map. The gravity settling time required for fertilizer particles to fall is calculated based on the soil penetration depth of the trencher and the target fertilization depth. Calculate the mechanism displacement time required for the fertilizer application actuator to reach the horizontal projection position of the target coordinates based on the current walking speed and the relative offset; The sum of the gravity settling time and the mechanism displacement time is set as the fertilizer discharge time delay, so that the fertilizer application mechanism performs the fertilizer discharge action according to the fertilizer discharge time delay when it reaches the target coordinate.

5. The automatic fertilization control program for potatoes based on digital twins according to claim 1, characterized in that, Calculating the fertilizer discharge pulse width of the fertilizer application actuator includes: Extract the difference between the deficit volume and the deficit concentration of the deficit grid to which the target coordinate belongs in the nutrient deficit spatial distribution map, and calculate the initial number of fertilizer discharge pulses by combining the preset single-pulse fertilizer discharge calibration value. Obtain the variance of the discharge velocity fluctuation at the discharge port of the fertilization actuator, input the variance of the discharge velocity fluctuation into the pre-constructed pulse width modulation compensation model, and output the pulse width correction coefficient for the initial number of discharge pulses; Multiply the initial number of fertilizer discharge pulses, the reference pulse width corresponding to the variance of fertilizer discharge velocity fluctuation, and the pulse width correction coefficient to obtain the target fertilizer discharge pulse width for the target coordinates; The fertilization actuator controls the opening and closing duration of the fertilization valve according to the target fertilization pulse width, and delivers the corresponding amount of fertilizer to the target coordinate position.

6. The automatic fertilization control program for potatoes based on digital twins according to claim 1, characterized in that, Following the generation of the targeted fertilization control command, the following is also included: The real-time feedback data of the actual fertilizer discharge from the fertilization execution mechanism and the real-time multispectral nutrient detection data of the soil profile after fertilization are obtained. The real-time fertilizer discharge feedback data and the real-time multispectral nutrient detection data are input into the coupled digital twin model to drive the coupled digital twin model to update the current soil nutrient concentration distribution and obtain the actual soil nutrient distribution state after fertilization. Compare the actual distribution of soil nutrients after fertilization with the target nutrient distribution corresponding to the spatial distribution map of nutrient deficit, and calculate the spatial concentration residual. If the spatial concentration residual exceeds a preset residual threshold, the root growth rate parameter and soil nutrient diffusion coefficient in the coupled digital twin model are corrected based on the spatial concentration residual, and the targeted fertilization control command is regenerated based on the corrected coupled digital twin model.

7. The automatic fertilization control program for potatoes based on digital twins according to claim 2, characterized in that, The soil bulk density parameter and the soil porosity parameter are introduced into the coupled digital twin model to construct a root growth spatial resistance field, including: The distribution value of soil penetration resistance is calculated based on the soil bulk density parameter and the soil porosity parameter; The soil penetration resistance distribution value is compared with the root critical penetration stress to delineate the root impenetrable boundary; The method of dynamically adjusting the root hair branch density based on the soil porosity parameter and soil moisture parameter of the region where the growth trajectory is located, and calculating the nutrient absorption interception radius based on the root hair branch density, includes: Extract the path length of the growth trajectory through the pores between the impenetrable boundaries of adjacent roots, and calculate the seepage pressure on the surface of the main root based on the path length and the soil moisture parameters. When the seepage pressure on the main root surface exceeds the osmotic pressure threshold of the root cortex cells, the density of the hairy root branches is triggered to increase exponentially along the seepage pressure gradient direction on the main root surface. The radial radius of the increased hairy root branch density is set as the nutrient absorption interception radius.

8. The automatic fertilization control program for potatoes based on digital twins according to claim 3, characterized in that, The historical fertilization data and soil moisture parameters are input into the soil solute transport convection-diffusion equation to calculate the time-varying nutrient concentration gradient of each spatial voxel in the three-dimensional soil space at the current moment, including: The isothermal adsorption and desorption parameters of soil particles for nutrient molecules are obtained, and the isothermal adsorption and desorption parameters are introduced as source and sink terms into the soil solute transport convection-dispersion equation. By combining the water infiltration flux in the soil moisture parameters, the soil solute transport convection-dispersion equation after introducing the source-sink term is solved to obtain the time-varying nutrient concentration gradient of free nutrient molecules in each spatial voxel; The calculation of the spatial distribution of the difference between the integral value of the time-varying nutrient concentration gradient within the target voxel set and the target fertilizer requirement threshold for potatoes includes: Within the target voxel set, the time-varying concentration gradient of the free nutrient ions is integrally divided along the normal direction of the nutrient absorption interception radius to obtain the effective nutrient integral that the root system can actually absorb. The spatial difference between the effective nutrient integral and the potato target fertilizer concentration threshold is calculated, and the spatial difference is mapped back to a three-dimensional mesh to generate the spatial distribution of the difference.

9. The automatic fertilization control program for potatoes based on digital twins according to claim 4, characterized in that, The calculation of the gravity settling time required for fertilizer particles to fall, based on the trencher's penetration depth and the target fertilization depth, includes: Obtain the particle size distribution parameters and particle density parameters of fertilizer particles, and construct a particle settling resistance model by combining the soil moisture parameters between the trencher's depth of penetration and the target fertilization depth. The particle size distribution parameter, the particle density parameter, and the soil moisture parameter are input into the particle settling resistance model to calculate the terminal settling velocity of fertilizer particles in the soil gaps. The gravity settling time is calculated based on the difference between the terminal settling velocity and the soil penetration depth of the furrow opener and the target fertilization depth. The calculation of the mechanism displacement time required for the fertilization actuator to reach the horizontal projection position of the target coordinates based on the current walking speed and the relative offset includes: The vibration frequency of the fertilizing machinery chassis and the elevation difference of the ground surface undulation are obtained. The current walking speed is corrected based on the vibration frequency of the chassis and the elevation difference of the ground surface undulation to obtain the actual effective walking speed. The displacement time of the mechanism is calculated based on the actual effective walking speed and the relative offset.

10. The automatic fertilization control program for potatoes based on digital twins according to claim 5, characterized in that, The variance of the discharge velocity fluctuation at the discharge port of the fertilization actuator is obtained, and the variance of the discharge velocity fluctuation is input into a pre-constructed pulse width modulation compensation model. The output is a pulse width correction coefficient for the initial number of discharge pulses, including: The fluctuation values ​​of the discharge shaft speed and the change value of the discharge port opening gap at the discharge port are collected in real time, and the variance of the discharge flow velocity fluctuation is calculated based on the fluctuation values ​​of the discharge shaft speed and the change value of the discharge port opening gap. The variance of the fertilizer discharge velocity fluctuation is matched and queried with a preset variance-coefficient mapping table. The variance-coefficient mapping table is constructed based on the nonlinear velocity compensation coefficients corresponding to different variance intervals. The pulse width correction coefficient is obtained by querying the table. The fertilization actuator controls the opening and closing duration of the fertilization valve according to the target fertilization pulse width, delivering the corresponding amount of fertilizer to the target coordinate position, including: The target fertilizer discharge pulse width is converted into a control level signal for the fertilizer discharge valve. When the fertilizer discharge time delay arrives, the control level signal is output to drive the fertilizer discharge valve to open for the specified opening and closing duration, so that the fertilizer particles are transported to the target coordinate position under the combined action of gravity and forced airflow.

Images