A method and system for hyperspectral fertilization control of greenhouse crops

Abstract

The present application relates to the field of facility agriculture and intelligent control technology, and particularly relates to a greenhouse crop hyperspectral fertilization control method and system, the method comprising: collecting hyperspectral images and illumination coding information, constructing a differential reflection cube, a canopy height map and a porosity map through radiation correction and registration; obtaining leaf surface nitrogen, leaf area index, leaf inclination angle distribution and cluster index through polarized spectral unmixing and radiation transmission inversion and giving parameter posterior covariance; applying a pseudo-random binary sequence perturbation within a control period, combining a root zone nitrogen budget model and a Gaussian process regression to obtain a growth gain function and a marginal response and generate a scene distribution; accordingly, a distribution robust probability constraint model is established to predict control output zoning fertilization settings, and a closed-loop evaluation data is used to update an active spectral coding matrix. The present application realizes stable and robust quality improvement and yield increase, reduces liquid discharge and nitrogen loss, and improves model distinguishability and long-term maintainability.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of facility agriculture and intelligent control technology, and particularly relates to a greenhouse crop hyperspectral fertilization control method and system. BACKGROUND

[0002] Facility planting is under triple pressure of yield stability, quality consistency and input cost. Excessive nitrogen input leads to increased conductivity and increased liquid discharge, resulting in resource waste and environmental burden; insufficient input causes growth stagnation and commodity decline. The production line generally relies on experience and fixed period irrigation formula, lacking simultaneous perception and verifiable causal response of canopy structure and root zone nitrogen dynamics, so it is difficult to maintain stable income and risk boundary under changing light and heat conditions and variety differences. The industry urgently needs to integrate high-information phenotyping and interpretable control into a deployable closed-loop system to achieve strong constraint fertilization decision-making and continuous self-calibration by partition and period. SUMMARY

[0003] In view of the problems existing in the prior art, the present application provides a greenhouse crop hyperspectral fertilization control method and system. The present application constructs a differential reflection cube based on active illumination and hyperspectral imaging, obtains leaf surface nitrogen, leaf area index, leaf inclination angle distribution and cluster index by combining polarization spectrum unmixing and radiation transfer inversion, and forms parameter posterior covariance; injects a pseudo-random binary sequence perturbation in the control period, determines the growth gain function and marginal response by the root zone nitrogen budget model and Gaussian process regression; generates partition fertilization settings in the distribution robust probabilistic constraint model predictive control and updates the active spectral coding matrix in a closed loop, achieving stable yield increase and salt control.

[0004] One or more embodiments of the present specification provide a greenhouse crop hyperspectral fertilization control method, comprising the following steps:

[0005] Collect hyperspectral images and illumination coding information, perform radiation correction and registration, construct a differential reflection cube, a canopy height map and a porosity map;

[0006] Based on polarization spectrum unmixing and radiation transfer inversion, determine leaf surface nitrogen, leaf area index, leaf inclination angle distribution and cluster index, and generate parameter posterior covariance;

[0007] Apply a pseudo-random binary sequence perturbation in the control period, combine the root zone nitrogen budget model and Gaussian process regression, determine the growth gain function and marginal response and generate the scene distribution;

[0008] Establish a distribution robust probabilistic constraint model predictive control according to the growth gain function, the scene distribution and the parameter posterior covariance, output partition fertilization settings, and update the active spectral coding matrix based on closed-loop evaluation data.

[0009] According to the method described in one or more embodiments of the present specification, the illumination coding information is composed of the combination of narrowband, polarization and illumination angle, the difference cube is constructed by the difference normalization of the reference illumination intensity and the no-illumination reference frame, and the canopy height map is reconstructed by three-dimensional imaging, and the porosity map is calculated by the area ratio of the local window through the vegetation mask.

[0010] According to the method described in one or more embodiments of the present specification, the polarization spectral unmixing adopts non-negative matrix decomposition and introduces polarization difference constraints to separate leaf end members, and the radiation transfer inversion is solved based on the joint two-layer model of leaf radiation transfer model and canopy radiation transfer model.

[0011] According to the method described in one or more embodiments of the present specification, the parameter posterior covariance is estimated based on the linearized sensitivity of the radiation transfer inversion and the observation noise to obtain the uncertainty measure of each inversion parameter.

[0012] According to the method described in one or more embodiments of the present specification, the timing of the pseudo-random binary sequence perturbation is designed under the constraints of amplitude and daily total input, and optimized to improve parameter distinguishability.

[0013] According to the method described in one or more embodiments of the present specification, the Gaussian process regression introduces the pseudo-random binary sequence as instrumental variables, and eliminates the correlation between the input and the external environment by constraining the correlation between the instrumental variables and the regression residuals to zero, thereby determining the growth gain function and the marginal response.

[0014] According to the method described in one or more embodiments of the present specification, the root zone nitrogen budget model combines nonlinear Kalman filtering to jointly estimate the state and update the covariance of the root zone inorganic nitrogen inventory and nitrogen uptake rate.

[0015] According to the method described in one or more embodiments of the present specification, the distribution robust probabilistic constraint model predictive control adopts uncertainty set formation probability constraints based on distribution distance metrics to cope with the uncertainty of scene distribution.

[0016] According to the method described in one or more embodiments of the present specification, the update of the active spectral coding matrix is based on the index reflecting the constraint tension degree in the closed-loop evaluation data and the difference between the growth gain prediction and the observed growth gain to construct the update target, and adjusts the parameters of the active spectral coding matrix according to the gradient rule.

[0017] Some embodiments of the present specification also provide a greenhouse crop hyperspectral fertilization control system for implementing the greenhouse crop hyperspectral fertilization control method, which comprises:

[0018] The acquisition processing module is configured to acquire hyperspectral images and illumination encoding information, perform radiometric correction and registration, construct a differential reflectance cube, a canopy height map, and a porosity map, and provide the differential reflectance cube, the canopy height map, and the porosity map to the inversion estimation module;

[0019] The inversion estimation module is configured to determine leaf area nitrogen, leaf area index, leaf inclination distribution, and clumping index based on polarized spectral unmixing and radiative transfer inversion, generate parameter posterior covariance, and provide the leaf area nitrogen, the leaf area index, the leaf inclination distribution, the clumping index, and the parameter posterior covariance to the identification modeling module and the control optimization module.

[0020] The identification modeling module is configured to apply a pseudo-random binary sequence perturbation within a control period, combine a root zone nitrogen budget model and Gaussian process regression, determine growth gain functions and marginal responses based on the leaf area nitrogen, the leaf area index, the leaf inclination distribution, the clumping index, and the parameter posterior covariance, generate a scene distribution, and provide the growth gain functions, the marginal responses, and the scene distribution to the control optimization module.

[0021] The control optimization module is configured to establish a distribution robust probabilistic model predictive control based on the growth gain functions, the marginal responses, and the scene distribution, as well as the parameter posterior covariance, generate a partitioned fertilization setting, and issue an execution, receive closed-loop evaluation data, and update an active spectral encoding matrix based on the closed-loop evaluation data.

[0022] Compared with the prior art, the application has the following advantages and beneficial effects:

[0023] By means of the differential reflectance cube and the polarized spectral unmixing, effective suppression of mixed pixels and specular components is achieved, and a leaf net reflectance spectrum that can be used for inversion is obtained, thereby solving the problem of sensitivity of traditional vegetation indices to structure and illumination.

[0024] By means of the radiative transfer inversion and the parameter posterior covariance, physical consistent estimation and uncertainty quantification of leaf area nitrogen, leaf area index, leaf inclination distribution, and clumping index are achieved, thereby solving the problems of uninterpretable black-box estimation and lack of credibility.

[0025] By means of the pseudo-random binary sequence perturbation, the root zone nitrogen budget model, and the Gaussian process regression, causal identification of growth gain functions and marginal responses is achieved, thereby solving the problem of bias estimation caused by input and environmental correlation.

[0026] By means of the distribution robust probabilistic model predictive control, robust fertilization settings under the boundaries of electrical conductivity and drainage are achieved, thereby solving the problems of over-limit and failure of traditional control under scene distribution drift.

[0027] By means of the active spectral encoding matrix update driven by closed-loop evaluation, synergistic evolution of a perception link and a control link is achieved, thereby solving the problems of model aging and information efficiency decline in long-term operation. Attached Figure Description

[0028] Figure 1 This is a schematic flowchart of the method of the present invention;

[0029] Figure 2 This is a structural block diagram of the system of the present invention. Detailed Implementation

[0030] The embodiments of the present disclosure will now be described with reference to the accompanying drawings. However, it should be understood that these descriptions are exemplary only and are not intended to limit the scope of the disclosure. In the following detailed description, numerous specific details are set forth to provide a thorough understanding of the embodiments of the present disclosure for ease of explanation. However, it will be apparent that one or more embodiments may be practiced without these specific details. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concepts of the present disclosure.

[0031] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit this disclosure. The terms “comprising,” “including,” etc., as used herein indicate the presence of the stated features, steps, operations, and / or components, but do not exclude the presence or addition of one or more other features, steps, operations, or components.

[0032] All terms used herein (including technical and scientific terms) have the meanings commonly understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein are to be interpreted in a manner consistent with the context of this specification, and not in an idealized or overly rigid way.

[0033] Greenhouse crops are grown in a controllable microclimate and substrate, whose production characteristics are not a passive process of "natural response", but an engineered process of "setting-feedback-correction". Crop growth is influenced by spectral absorption, scattering and canopy structure, and physiological parameters such as leaf nitrogen content, water content and leaf area index are observable in the spectral domain; at the same time, root zone inorganic nitrogen is actively regulated by fertilizer solution ratio and irrigation rhythm, and directly affects photosynthesis and assimilate accumulation. Unlike open fields, the covering material, lighting geometry and air flow organization of greenhouses make the reflectance spectrum and process variables such as conductivity and drainage rate strongly coupled and strongly time-varying, and it is difficult to obtain a stable input-output relationship by relying on empirical thresholds or a single index. Therefore, a closed-loop path is needed to convert "observable" into "controllable": high-information phenotyping of canopy and leaves using hyperspectral, obtaining interpretable nitrogen-related parameters using physically consistent inversion, characterizing the marginal response of input to gain using small controllable perturbations and statistical learning, and then rolling optimization under safety and resource constraints. The present invention is around the engineering properties of greenhouse crops, and proposes a hyperspectral-driven fertilization control system, so that spectral information is not only for monitoring, but also for zoning, time-based and verifiable decision-making basis.

[0034] As shown in Figure 1 A greenhouse crop hyperspectral fertilization control method, comprising the following steps:

[0035] Collecting hyperspectral images and lighting encoding information, performing radiometric correction and registration, constructing a differential reflectance cube, a canopy height map and a porosity map;

[0036] The lighting encoding information is composed of a combination of narrowband, polarization and illumination angle, the differential reflectance cube is constructed by differential normalization of reference illumination intensity and no-illumination reference frame, the canopy height map is reconstructed by three-dimensional imaging, and the porosity map is calculated by the area ratio of the local window under the vegetation mask.

[0037] The present invention acquires basic data for subsequent inversion and control in a greenhouse in an active lighting and synchronous imaging manner. The lighting encoding information is composed of three elements of narrowband, polarization and illumination angle. The narrowband is used to select the reflection response of a specific wavelength, the polarization is used to suppress mirror interference and distinguish the structure, and the illumination angle is used to change the incident geometry to enhance the observability of the canopy structure. The control end triggers the light source and polarizer one by one according to the preset encoding sequence, synchronously triggers the camera to collect the hyperspectral image, records the reference illumination intensity and the no-illumination reference frame, and ensures that each exposure has complete radiation and geometric metadata.

[0038] The radiation correction link in turn carries out dark field deduction, response linearization, reference illumination intensity normalization and whiteboard consistency calibration. The dark field deduction is used to eliminate the background signal of the photosensitive device, the response linearization is used to eliminate the proportional deviation introduced by gain and exposure, the reference illumination intensity normalization is used to offset the short-term fluctuations of the light source, and the un-illumination reference frame is used to eliminate the background and stray light. The whiteboard consistency calibration obtains the gain coefficient across the narrow band by comparing the whiteboard reflection of the current narrow band with the reference whiteboard reflection, so as to keep the consistency of the data between different days.

[0039] The cross-band and cross-polarization registration takes a stable near-infrared narrow band as a reference. First, the initial affine relationship is estimated by using the feature matching based on the corner points and textures, and then the sub-pixel refinement is performed by using the phase correlation, and the image sequence in the unified coordinates is output. Subsequently, the differential reflection at the same pixel position under different narrow bands, different polarizations and different illumination angles is stacked in a fixed order to form a differential reflection cube. The differential reflection cube and the illumination encoding information are stored one by one, including the timestamp and the illumination intensity record, which is convenient for tracing and recalculation.

[0040] The canopy height map is derived from a three-dimensional imaging device. After the time-of-flight camera or the structured light camera outputs the depth map, the depth map is projected to the imaging coordinates according to the external parameter relationship, and is aligned with the registered image to obtain the canopy height map with the same grid as the differential reflection cube. For the depth hole region, the neighborhood interpolation is used to generate a quality mask to guide the subsequent weighted processing. The porosity map is used to describe the proportion of non-vegetation voids in the canopy. First, the vegetation mask is generated by threshold segmentation and morphological operation on the near-infrared narrow band, and then the non-vegetation area ratio in a fixed window centered on each pixel is calculated, and the result is taken as the porosity map. The porosity map and the canopy height map together constitute the prior constraint on the canopy structure, and directly serve the subsequent radiation transfer inversion and control modeling.

[0041] The data processing module encapsulates the differential reflection cube, the canopy height map, the porosity map, the quality mask and the timestamp into a single data object, and the coordinate system, the resolution and the time axis in the object are completely consistent, so that the subsequent algorithms do not need to be aligned again. The processing flow is realized in a pipeline manner, supports batch processing and online incremental processing, triggers a degradation strategy when the illumination is abnormal or the image is saturated, automatically removes the abnormal narrow band or replaces it with a backup narrow band, and records the alarm information in the metadata.

[0042] Embodiment: Nine narrowband light sources, one motorized polarizer and one snapshot hyperspectral camera are arranged on a greenhouse trolley, the three-dimensional imaging device is fixed on the same bracket. The inspection timing is that the nine narrowband top light triggers in turn, the two groups of angle side light triggers in turn, and the reference frame is collected after each trigger ends and the reference illumination intensity is recorded. After the collection is completed, the data processing module performs dark field subtraction, response linearization, reference illumination intensity normalization and whiteboard consistency calibration, completes cross-band and cross-polarization registration, and outputs the differential reflection cube. The depth map generated by the three-dimensional imaging device is projected to the imaging coordinate through the calibration external parameter, and the crown layer height map is output. The near-infrared narrowband generates a vegetation mask, and the non-vegetation area ratio is calculated in a fixed-size sliding window to obtain the porosity map. Finally, the differential reflection cube, the crown layer height map and the porosity map are stored together with the quality mask and the time stamp as the standard input for subsequent inversion and control.

[0043] Based on the polarization spectrum unmixing and radiation transmission inversion, the leaf surface nitrogen, leaf area index, leaf inclination distribution and cluster index are determined, and the parameter posterior covariance is generated;

[0044] After the differential reflection cube, the crown layer height map and the porosity map are completed, the parameter estimation stage based on the polarization spectrum unmixing and the radiation transmission inversion is entered. The target is to obtain the leaf surface nitrogen, the leaf area index, the leaf inclination distribution and the cluster index, and to give the parameter posterior covariance for the risk measurement of subsequent identification and control.

[0045] The polarization spectrum unmixing is used to separate the leaf net reflectance spectrum from the mixed pixels. The specific method is as follows: the differential reflection data of two polarization states are read at the same spatial position, a non-negative decomposition model with leaf, soil and non-leaf vegetation as end members is constructed, the end member coefficient is iteratively optimized, and the specular component is suppressed through polarization difference constraint to obtain the leaf net reflectance spectrum and the corresponding net spectrum weight. The net spectrum weight is used as the pixel reliability weight for subsequent inversion, which is used to reduce the influence of mixed and shadow areas on the results.

[0046] The radiation transmission inversion adopts a two-layer combination of leaf optical model and crown scattering model. In order to improve the efficiency, the spectral library covering chlorophyll content, dry matter content, leaf area index, leaf inclination distribution and cluster index is generated offline, the coarse initial value is obtained by nearest neighbor search in the online stage, and the weighted least squares search is carried out near the initial value, and the constraint range is given by the crown height map and the porosity map. The parameters output by the inversion include the leaf area index, the leaf inclination distribution and the cluster index, and the leaf optical parameters are also output.

[0047] The leaf surface nitrogen is obtained through the leaf optical parameters and experimental calibration. The calibration establishes a mapping by periodically collecting leaf samples, and the core expression is:

[0048] ;

[0049] for leaf area nitrogen, for chlorophyll content, for dry matter content. The functional form is determined by sampling regression of the same greenhouse variety and management condition, and is reviewed by batch in subsequent operation. To quantify the reliability of the parameters, the present application calculates the posterior covariance of the parameters, and the core expression is:

[0050]

[0051] for the posterior covariance of the parameters, for the sensitivity matrix of the parameters, for the observation noise covariance. The sensitivity matrix is obtained by numerical differentiation of the forward model, and the observation noise covariance is obtained by converting the imaging noise and the net spectrum weight. The posterior covariance of the parameters is summarized at the partition scale, which is used for risk constraints and weight setting in subsequent control.

[0052] Embodiment: Take the near-infrared main band as the registration reference, select four red edge and near-infrared narrow bands as the inversion main band, and initialize the polarization unmixing end member by small block statistics of net leaf area, bare soil area and flower organ area. The offline spectrum library covers the common cultivation range, and the Euclidean distance is used to retrieve a number of candidates in the online stage, and then the local least squares search is performed with the net spectrum weight as the weighted factor and the structure prior is applied to obtain the leaf area index, leaf inclination distribution and cluster index. Then read the leaf sample regression model of the same batch to calculate the leaf area nitrogen and output the posterior covariance of the parameters. Finally, a rasterized product containing leaf area nitrogen, leaf area index, leaf inclination distribution, cluster index and posterior covariance of the parameters is formed, which is used as the input of causal identification and model predictive control.

[0053] Polarized spectral unmixing uses non-negative matrix decomposition and introduces polarization difference constraints to separate leaf end members, and radiation transfer inversion is based on the joint two-layer model of leaf radiation transfer model and canopy radiation transfer model for solving.

[0054] After obtaining the differential reflection cube, canopy height map and porosity map, the present application enters the joint stage of polarized spectral unmixing and radiation transfer inversion, realizes the separation of leaf end members and estimates the canopy parameters. The overall process is: first, complete the polarized spectral unmixing at the pixel or small block scale to obtain the leaf net reflectance and net spectrum weight; then, perform weighted inversion by using the two-layer radiation transfer model to obtain the leaf area index, leaf inclination distribution and cluster index, and derive the leaf area nitrogen accordingly.

[0055] Polarized spectral unmixing uses constrained non-negative decomposition. For two polarization state data at the same spatial position, pure leaf and bare soil small blocks are used to initialize end members, and pixels with saturation and shadow are removed according to the quality mask. The core optimization problem is: ​

[0056] ;

[0057] is the stacked spectrum observation matrix; is the vertical polarization observation matrix; is the parallel polarization observation matrix; is the endmember matrix; is the endmember response in vertical polarization; is the endmember response in parallel polarization; is the endmember coefficient matrix; is the polarization difference constraint weight. The solution adopts alternating minimization and non-negative projection, and the iteration threshold and maximum round are uniformly managed by the quality monitoring module. The output leaflet net reflectance spectrum is recorded by waveband, and the net spectrum weight is used for weighted and elimination in subsequent inversion.

[0058] The radiation transfer inversion adopts a joint two-layer model of leaflet radiation transfer model and canopy radiation transfer model. In the offline stage, a spectrum library is generated based on the cultivation range, and an index is established to support fast retrieval. In the online stage, the following steps are implemented: first, a number of candidates are retrieved in the spectrum library based on the leaflet net reflectance spectrum to obtain initial values; second, the search range of leaf area index, leaf inclination angle distribution and cluster index is limited according to the canopy height map and porosity map; third, local weighted minimization is performed with the net spectrum weight as the observation weight. The necessary objective function is:

[0059] ;

[0060] is the set of estimated parameters; is the leaflet net reflectance value of the i-th waveband; is the simulation value of the two-layer radiation transfer model in the i-th waveband; is the waveband weight; is the prior constraint weight; is the range penalty derived from the canopy height map and the porosity map. The optimization adopts local search and takes the residual threshold and the upper limit of iteration as the stopping condition. The output parameter is rasterized and stored, and is accompanied by a quality identifier for use by the control stage. The mapping of leaflet surface nitrogen through field calibration is obtained from the leaflet optical parameters. Samples are periodically collected under the same greenhouse conditions, and a regression mapping is established through laboratory determination, which is reviewed by batch in operation. To ensure implementability, the mapping model, spectrum library version and inversion hyperparameters are uniformly recorded in the metadata, supporting traceability and version rollback.

[0061] The mapping of leaflet surface nitrogen through field calibration is obtained from the leaflet optical parameters. Samples are periodically collected under the same greenhouse conditions, and a regression mapping is established through laboratory determination, which is reviewed by batch in operation. To ensure implementability, the mapping model, spectrum library version and inversion hyperparameters are uniformly recorded in the metadata, supporting traceability and version rollback.

[0062] ​The system architecture is composed of an unmixing sub-module and an inversion sub-module. The unmixing sub-module completes data loading, quality mask application, endmember initialization and constraint decomposition; the inversion sub-module completes candidate retrieval, prior range setting and local weighted optimization. The two sub-modules share a unified coordinate system, grid and timestamp, avoiding repeated registration and interpolation. In terms of exception handling, a degradation strategy is adopted for low net spectral weight or high residual pixels, and only partition statistics are output, avoiding the introduction of unstable pixels.

[0063] Embodiment: In a multi-span greenhouse, four main bands of red edge and near infrared are selected as the inversion main bands. The pure leaf and bare soil endmember initial values are obtained from a small block in the near-infrared stable region, and the leaf area index, leaf inclination distribution and cluster index are generated by unmixing and generating leaf net reflectance and net spectral weight in 32x32 pixel blocks. Then ten nearest candidates are retrieved from the spectral library as initial values, combined with the canopy height map and porosity map to limit the parameter range, and the weighted minimization is performed to obtain the leaf area index, leaf inclination distribution and cluster index. The leaf surface nitrogen is calculated by the regression model of the sample at that time. Finally, the data object containing the parameter grid and quality identification is output, which is directly called for subsequent identification and model prediction control.

[0064] The parameter posterior covariance is estimated based on the linearized sensitivity of the radiation transfer inversion and the observation noise to obtain the uncertainty measure of each inversion parameter.

[0065] After obtaining the leaf net reflectance, canopy height map and porosity map and completing the two-layer radiation transfer inversion, the uncertainty measure of the inversion parameter for each grid position is needed. For this purpose, the parameter posterior covariance estimation method based on linearized sensitivity and observation noise is used, so that the subsequent identification and control can reduce the weight or tighten the constraint in the high uncertainty area.

[0066] First, the input data and parameter set for uncertainty evaluation are determined. The input data includes the reflectance value of each band of the leaf net reflectance, the band signal-to-noise ratio, the pixel quality mask and the inversion residual. The parameter set is the chlorophyll content, the dry matter content, the leaf area index, the leaf inclination distribution parameter and the cluster index. For each inversion pixel or superpixel, an observation noise covariance matrix is established. The method is as follows: the variance of each band is estimated by the camera read noise and illumination intensity fluctuation, and the variance is converted by combining the weight of the leaf net reflectance and the inversion residual, to obtain the diagonal matrix form of the observation noise covariance; when the pixel quality mask shows saturation, shadow or endmember impurity, the variance of the corresponding band is enlarged according to the rules, and its influence is weakened.

[0067] The sensitivity matrix of the reflectance to the parameters is then calculated. The method is as follows: in the vicinity of the obtained parameter point in the inversion, each parameter is perturbed one by one according to a physically feasible small step, the two-layer radiation transfer forward model is called, and the derivative of each band reflectance to each parameter is obtained by using the difference method to form the sensitivity matrix. The step size follows a fixed proportion of the parameter value range, and the boundary is clipped to ensure that the perturbation does not go out of bounds. The calculation of the sensitivity matrix uses the same band combination and geometric settings as the inversion to ensure consistency.

[0068] After obtaining the observation noise covariance and the sensitivity matrix, the parameter posterior covariance is calculated. The core calculation expression is:

[0069]

[0070] is the parameter posterior covariance, is the sensitivity matrix of the reflectance to the parameters, is the observation noise covariance. To improve numerical stability, the sensitivity matrix is first normalized by the band scale, and then the inverse operation is performed; when the condition number exceeds the threshold, a small diagonal regularization term is introduced and the quality identifier is recorded for subsequent control stage degradation processing.

[0071] The output results include: the variance of each parameter and the covariance information of each pair of parameters, the confidence interval based on the covariance, and the partition-level summary value for control. The pixel-level variance is robustly smoothed by a spatial sliding window to form a partition uncertainty layer; the partition-level summary adopts area-weighted average and robust quantile parallel storage to meet the needs of different control strategies. All results and parameter inversion layers are aligned on the same coordinate and time axis, and the metadata is written, including the sensitivity step size, normalization scale and regularization term record, to ensure traceability.

[0072] Embodiment: In a management partition of a multi-span greenhouse, four red edge and near-infrared bands are selected as evaluation bands. First, the observation variance of each band is calculated according to the camera calibration data and illumination monitoring value of the day, and the corresponding variance is amplified for end-member impure pixels. Then the chlorophyll content, dry matter content, leaf area index, leaf inclination distribution parameter and cluster index are perturbed one by one, the forward model is called to obtain the reflectance change, and the sensitivity matrix is constructed. Finally, the parameter posterior covariance is calculated according to the above expression, the variance map of each parameter and the partition-level variance statistics are output, and the condition number and regularization term usage status are attached. The results are used as sample weights by the subsequent identification module, and as probability constraints and risk weights by the control module, so as to maintain the yield target while suppressing excessive fertilization in high-uncertainty areas.

[0073] ​A pseudo-random binary sequence perturbation is applied to the fertilizer input within a control period, combined with a root zone nitrogen budget model and Gaussian process regression to determine the growth gain function and marginal response and generate scenario distributions;

[0074] The present application applies a light perturbation to the fertilizer input within a control period, combined with a root zone nitrogen budget and regression learning to obtain the growth gain function and marginal response, and simultaneously constructs scenario distributions for prediction. The perturbation uses a pseudo-random binary sequence, which switches within the constraints of crop safety and daily input to increase or decrease by a small amount, and is injected by the controller according to the plan, with timestamps, execution amounts, and sensor data being recorded.

[0075] The fertilizer input is represented by the product of flow rate and concentration:

[0076] ;

[0077] The nitrogen application flux at time ; The irrigation flow rate; The inorganic nitrogen concentration in the nutrient solution. The control module triggers the perturbation when the light is sufficient, the substrate water content is stable, and the electrical conductivity is within the safe range, and the amplitude and frequency are determined by the safety constraints and information benefits. All execution deviations and alarms are written into the metadata.

[0078] The root zone nitrogen budget is used to estimate the root zone inorganic nitrogen inventory and nitrogen uptake rate, and the discrete expression based on mass conservation is:

[0079] ;

[0080] The root zone inorganic nitrogen inventory; The return flow rate; The inorganic nitrogen concentration in the return liquid; The nitrogen uptake rate. The return liquid conductivity and sampling test are used to correct the return liquid concentration sequence, and recursive filtering is used to estimate the inventory and nitrogen uptake rate, and consistency check is performed with the partitioned statistics of leaf surface nitrogen; when abnormal, the weight of this partition is reduced and the perturbation is delayed.

[0081] The growth response uses the increment of short-term growth indicators as the output, and the input is the fertilizer input and environmental state, and Gaussian process regression is used to learn the function relationship:

[0082] ;

[0083] ;

[0084] The growth indicator increment within the next observation window from time ; The growth gain function; marginal response; are environmental and state vectors, including photosynthetically active radiation, air temperature, air humidity, root zone inorganic nitrogen stock, and leaf area nitrogen. The training data is weighted by pixel quality and parameter posterior covariance, and pseudo-random binary sequences are used to provide independent excitation to reduce the correlation between inputs and the environment. The marginal response is directly calculated by the derivative interface of the model and is output together with the stock estimate as a time series.

[0085] Scenario distribution is used for control optimization in the prediction period. The system extracts external condition trajectories from recent residuals and environmental history, and performs perturbation sampling on root zone inorganic nitrogen stock and model parameters to generate a limited number of representative scenarios; each scenario contains an environmental sequence, a root zone inorganic nitrogen stock trajectory, and a corresponding uncertainty weight. The scenario and the growth gain function are packaged together as a standardized data object for subsequent model predictive control.

[0086] Embodiment: The control cycle is set to 15 minutes, and the fertilization input is switched at a small amplitude of no more than the baseline percentage during the sunshine period, each perturbation lasts for 30 to 45 minutes and is separated from adjacent perturbations; the flow and conductivity of irrigation and return liquid are recorded by online sensors, and the return liquid concentration is sampled daily for bias correction; recursive filtering outputs the root zone inorganic nitrogen stock and nitrogen uptake rate, and abnormal partitions automatically reduce the disturbance frequency; the regression model is updated with data from the past two weeks, with the input being the fertilization input and environmental state, and the output being the growth gain and marginal response; scenario distribution is generated daily for the next day's control solution. The above processes are aligned on the same coordinate and time axis, forming a traceable and reproducible data link and model link.

[0087] The timing of the pseudo-random binary sequence perturbation is designed under the constraints of amplitude and daily total input, and is optimized to improve parameter identifiability.

[0088] The present application superimposes a pseudo-random binary sequence perturbation on the fertilization input within the control cycle, aiming to improve parameter identifiability under safety and input constraints. The perturbation is generated by the planner, issued by the executor, and closed-loop verified by the monitor, forming a complete link of "design - execution - evaluation". The design principle is: amplitude limited, daily total input limited, executor change rate limited, only excitation within the safety window, and the safety window is determined by conductivity, substrate moisture content, and photosynthetically active radiation.

[0089] The fertilization input is in the form of baseline plus perturbation:

[0090] ;

[0091] is the nitrogen application flux at time is the nitrogen application flux at time ​​the baseline nitrogen flux, is the perturbation amplitude, is a binary sequence taking values -1 and +1. The perturbation amplitude does not exceed the allowed variation range of the device and the crop, and the flipping interval of the binary sequence satisfies the minimum residence time. The daily total input satisfies the constraint, and if necessary, it is checked using the following formula: where, is the control period, is the upper limit of the daily total nitrogen.

[0092] To maximize the recognition benefit of the perturbation, the planner calculates the information gain score for each candidate period and optimizes the binary sequence under the constraints. The information gain score is based on the sensitivity and noise level of the current learning model:

[0093] ;

[0094] is the information gain score, is the time window covered by this design, is the safety weight, is the first-order sensitivity of the growth gain function to the nitrogen flux, is the corresponding observation noise standard deviation. The safety weight is calculated from the conductivity margin and the substrate water content margin, and the sensitivity and noise are given by the trained regression model and the sensor historical residual. The planner selects the flipping time and duration in a greedy or rolling search manner, maximizing the score without violating the amplitude, total amount and rate constraints. When the score is below the threshold or the safety window is insufficient, the planner skips the perturbation of this partition.

[0095] In the execution phase, the controller switches the flow and concentration settings according to the binary sequence, and all actual execution amounts, alarms and loop pressures are returned in real time. The monitor performs out-of-limit judgment on the conductivity and effluent nitrogen concentration, and if the trend approaches the threshold, it will pause the subsequent flip and revert to the baseline. After each perturbation, the system records the effective excitation duration and actual amplitude, and packages the data for subsequent modeling.

[0096] To adapt to multiple partitions, the system allocates daily budgets and available windows by partition, and prioritizes scheduling perturbations in partitions with high information gain scores, scarce past data or large model uncertainty. Cross-partition perturbations are avoided from occurring simultaneously to prevent coupling effects on the liquid supply circuit. The planner maintains a blacklist period, including greenhouse opening and closing shed, transplanting and pesticide spraying periods, to ensure that the perturbation does not interfere with production.

[0097] Embodiment: In a management partition, the system first calculates the safety window and noise level from the data of the past week to form the available time slices within the day; then based on the sensitivity curve output by the current learning model, it calculates the information gain score by time slice and generates a binary sequence, sets the minimum residence time and the maximum number of flips; during the execution, the conductivity and liquid return flow trigger threshold value pauses the perturbation and records the reason; at the end of the day, the actual executed binary sequence, the effective incentive duration and the safety events are summarized for the next day planner to update. This process provides samples with discrimination for modeling without changing the production rhythm, making subsequent parameter estimation and marginal response fitting more stable.

[0098] Gaussian process regression introduces pseudo-random binary sequence as instrumental variable, eliminates the correlation between input and external environment by constraining the correlation between instrumental variable and regression residual to zero, and determines the growth gain function and marginal response.

[0099] The present application uses pseudo-random binary sequence as instrumental variable in the control cycle to inject small perturbations to the fertilizer input, and trains Gaussian process regression in a two-stage idea to obtain growth gain function and marginal response, while avoiding the bias caused by the correlation between fertilizer input and external environment. The data is first aligned on the same time axis, including nitrogen flux, environmental state vector and short-term growth index increment. Abnormal samples are removed or down-weighted according to quality mask, and all quantities are standardized by partition.

[0100] The two-stage training steps are as follows. In the first stage, the instrumental variable and the environmental state are used as independent variables to fit the nitrogen flux by linear regression or regular regression to obtain the external component estimate, denoted as predicted nitrogen flux. This estimate is only used to drive the second stage, and the original nitrogen flux does not directly enter the regression to cut off the inherent correlation between input and environment. In the second stage, the predicted nitrogen flux and the environmental state are used as independent variables, and the growth index increment is used as the dependent variable to train the Gaussian process regression. The kernel function uses the three-second order form of automatic correlation length radial basis or Mahalanobis kernel, and the noise term is independent and identically distributed Gaussian. The hyperparameters are obtained by maximizing the marginal likelihood, and are jointly optimized with the following orthogonal constraint penalty:

[0101] ;

[0102] is the growth index increment, is the predicted nitrogen flux, is the environmental state vector, is the Gaussian process hyperparameter, is the penalty weight, is the number of samples, is the instrumental variable sequence item, The regression residual is returned. The penalty will press the correlation between instrumental variables and regression residuals to zero, thus eliminating the influence of input and external environment.

[0103] The marginal response is obtained by the input derivative of the Gaussian process:

[0104]

[0105] The marginal response at time t is, The growth gain function is, The predicted nitrogen flux at time t is, The environmental state vector at time t is. The derivative is calculated by the analytical gradient of the kernel function, and is used together with the parameter posterior covariance to form the uncertainty index for subsequent control weighting.

[0106] The system architecture includes data management, two-stage modeling and model evaluation. The data management is responsible for alignment and standardization, and outputs the training window. The two-stage modeling module is updated in batch at each day end, and automatically reduces the model weight and prompts to increase the disturbance period when the instrumental variable intensity is insufficient. The model evaluation includes prediction residual monitoring and drift detection, triggering the rollback to the previous version model or shortening the training window.

[0107] In the embodiment, the control period is set to fifteen minutes, and the training window is nearly two weeks. The instrumental variable sequence comes from the predetermined disturbance plan, the nitrogen flux is fitted by the sequence and the environmental state to obtain the predicted nitrogen flux, and the growth gain function and the marginal response are obtained by training the Gaussian process. The model and uncertainty summary generated every day are written into the metadata for model predictive control and disturbance plan update the next day.

[0108] The root zone nitrogen budget model combines nonlinear Kalman filtering to jointly estimate and update the covariance of the root zone inorganic nitrogen inventory and nitrogen uptake rate.

[0109] After obtaining the leaf net reflectance spectrum, canopy height map and porosity map and completing parameter inversion, the application establishes a joint estimation process of “nitrogen budget + nonlinear Kalman filtering” for root zone nutrient dynamics, outputs the root zone inorganic nitrogen inventory and nitrogen uptake rate, and gives the covariance for subsequent control weighting and risk constraint. The data channels include: flow and concentration of fertilizer input, return liquid flow and conductivity, experimental correction results of return liquid nitrogen concentration, time series of partition leaf surface nitrogen, and environmental quantities (photosynthetically active radiation, temperature and humidity, and substrate water content). The online link first converts the return liquid conductivity into the return liquid nitrogen concentration by the calibration curve of the day, and aligns it with other channels according to the time axis, and removes saturated and sensor packet loss samples.

[0110] The process model uses discrete nitrogen budget, and the core expression is:

[0111] ​ ;

[0112] where, is the root-zone inorganic nitrogen stock at time , is the nitrogen application flux, is the return flow rate, is the return flow nitrogen concentration, is the nitrogen uptake rate. The nitrogen uptake rate is treated as a slowly varying process, evolving smoothly with environmental and physiological states, and is adaptively tracked by a filter with process noise.

[0113] The observation model consists of three types of available information: first, the return flow concentration and rate, directly constraining the budget terms; second, the substrate conductivity and water content, as indirect observations of the stock, mapped to the stock via daily calibration curves as soft constraints; third, the leaf area nitrogen partition increments, integrated over a sliding window with the nitrogen uptake rate to form consistency constraints, used to suppress long-term drift. Each observation channel is weighted by a quality weight, which is automatically updated by the sampling confidence and residual history.

[0114] The filter implementation uses a nonlinear Kalman filter. The state quantities are the root-zone inorganic nitrogen stock and the nitrogen uptake rate; the initialization is based on the previous day's end value and the opening shed baseline; the time update is advanced according to the nitrogen budget, with the covariance growing with the process noise; the measurement update fuses the three types of observations and performs threshold discrimination, with the out-of-limit samples being downgraded with robust weights. To ensure physical feasibility, the stock and the nitrogen uptake rate are non-negatively constrained by a log or soft threshold function, and are projected back to the allowed range after the update. The innovation statistics are used to adaptively tune the process noise and the observation noise; when the innovation variance is consistently large, the process noise is automatically relaxed to enhance the tracking, and vice versa to reduce the chattering.

[0115] In system implementation, the control period is taken as ten to fifteen minutes. Each period completes the data alignment, curve conversion, one-time time update and measurement update, and writes the state, covariance, innovation and quality identifier into the metadata. At the end of the day, the return flow nitrogen concentration curve is corrected for bias using laboratory tests, and the data of the day is played back to fine-tune the noise level. The abnormal processing includes: when the sensor is interrupted, only the available channel is retained and the covariance is increased; when the continuous out-of-limit occurs, the change of the nitrogen uptake rate is frozen and the manual review is prompted.

[0116] Example: In a zone of a tomato greenhouse, a baseline irrigation strategy is established, and small disturbances are inserted according to the planner. The liquid nitrogen concentration is converted back online using the daily electrical conductivity-concentration curve. A nonlinear Kalman filter is run according to the control cycle to obtain the root zone inorganic nitrogen inventory, nitrogen uptake rate, and their covariance. When the afternoon light intensity rapidly increases, leading to a larger innovation variance, the system adaptively increases the process noise, and the inventory estimation converges within two to three cycles. At the end of the day, the conversion curve is corrected using laboratory data, and the data is replayed. The next day, estimation is initiated using updated noise parameters. This process ensures that the inventory and nitrogen uptake rate are traceable, quantifiable, and meet physical constraints.

[0117] Based on the growth gain function, scene distribution, and posterior covariance of parameters, a bibliometric probability constraint model is established for predictive control, outputting zonal fertilization settings, and updating the active spectral coding matrix based on closed-loop evaluation data.

[0118] This invention, after obtaining the growth gain function, scene distribution, and posterior covariance of parameters, employs rolling time-domain model predictive control to output zonal fertilization settings, and drives the update of the active spectral coding matrix with closed-loop evaluation results. All data are aligned over the control cycle, with nitrogen flux as the control variable, and the execution quantity determined by pump speed and nutrient solution ratio.

[0119] The scene generator extracts representative scenes from the scene distribution, including photosynthetically active radiation, air temperature, air humidity, root zone inorganic nitrogen inventory trajectory, and model noise in the future time domain. For each scene, the controller calls the growth gain function to obtain gain prediction and process variable prediction, and uses the posterior covariance of the parameters to weight the uncertainty of the gain prediction. To ensure online solvability, the process variables are linearized with respect to the nitrogen flux near the solution of the previous time step to form an approximate simulation model.

[0120] The core objective is to achieve a risk-compensated stage return, and the necessary expression is:

[0121] ;

[0122] For the rolling time domain objective function, To predict the number of steps, For a moment Gain prediction, For a moment Nitrogen flux, For the benefit weight, To assign weights, For risk weights. The standard deviation is a composite of the regression model's predicted variance and the posterior covariance of the parameters.

[0123] The scenario approximation of safety and resource constraints using probabilistic form has the necessary expression as follows: ,in, For a moment conductivity prediction, is an upper limit of the conductivity, is an allowed violation probability. In implementation, the constraint is approximated by scenario quantile, while device boundary, rate of change and daily total constraints are imposed. Control variables are parameterized by piecewise constant to reduce the scale of solving, and a solver with warm start quadratic programming or second order cone programming is used with maximum iteration and failure fallback strategy.

[0124] The execution stage maps the nitrogen application flux of the last beat to pump speed and ratio instructions, and reads the flow and conductivity for deviation correction. All predictions and measurements are timestamped to form closed-loop evaluation data. Evaluation data is used for both controller self-checking and active spectral encoding matrix updating.

[0125] The update of the active spectral encoding matrix takes prediction consistency and constraint tension as indicators, and the necessary expression is:

[0126] ;

[0127] is an evaluation loss, is a gain prediction, is an observation gain, is a constraint relaxation, is a weight. The encoding matrix is fine-tuned along the negative gradient direction according to the step strategy, only adjusting the narrowband and polarization combination involved in imaging, and taking effect in the next day's inspection.

[0128] Embodiment: In a tomato greenhouse, the control period is fifteen minutes, and the prediction step is sixteen steps. A batch of scenes is generated every morning, and the controller uses the solution at the last time as the linearization point to obtain the zoned fertilization settings and issue them. When the boundary is triggered, the feasible region is automatically contracted and the conservative solution is degraded. At the end of the day, the prediction error and the constraint relaxation are summarized, and the narrowband and polarization weights of the active spectral encoding matrix are optimized, and the next day the imaging and inversion are carried out according to the new encoding. The above process stably improves the gain under the premise of ensuring the conductivity and the discharge boundary, and gradually reduces the model uncertainty.

[0129] The distributed robust probabilistic constraint model predictive control uses an uncertainty set based on distance measurement to form a probability constraint to cope with the uncertainty of the scene distribution.

[0130] The present application introduces the growth gain function, scene distribution and parameter posterior covariance into model predictive control, and constructs a probability constraint based on the uncertainty set of distance measurement, so as to maintain the robustness of the fertilization decision under external conditions and model residual fluctuations. The implementation process includes five links of scene generation, distance radius calibration, constraint transcription, online solving and closed-loop evaluation.

[0131] First, the scenario distribution fitted by recent data is taken as the benchmark, and a sample set containing environmental trajectories and root zone states is generated, with the source time, residual, weight, and parameter posterior covariance information attached to each sample. Then, according to the prediction residual statistics and parameter posterior covariance within the window, the distribution distance radius is calibrated, which is used to describe the upper limit of the "true distribution relative to the benchmark distribution". The uncertainty set is defined under this radius and used to limit the worst-case scenario of the probability constraint.

[0132] The core constraint takes "still meet the safety probability under the worst distribution" as the criterion, and the necessary expression is: where, is the distribution to be evaluated, is the uncertainty set defined by the distribution distance metric and radius, is the constraint expression of the safety function (such as conductivity, upper limit of drainage, etc.) composed of nitrogen flux and scenario variables, is the allowed violation probability. To facilitate online calculation, this constraint is approximated by the "worst quantile": re-weighting on samples and adding a conservative inflation amount to ensure that the probability threshold is met. The re-weighting is allocated according to the sample residual and the freshness of the source, and the inflation amount is given by the parameter posterior covariance and historical fluctuations.

[0133] After the target and constraint are converted, a standard rolling optimization problem is formed: the target uses the weighted sum of revenue, input, and risk compromise; the constraints include device boundaries, change rates, daily total, and the distribution robust approximation of the above probability constraints. To adapt to real-time, the process variables are linearized around the last round solution, the control variables are parameterized as piecewise constants, and the solver uses a hot-start quadratic programming or second-order cone programming; if the feasible region shrinks leads to failure, then automatically switch to a conservative alternative strategy and record the reason.

[0134] In terms of system architecture, the controller contains four sub-modules: the sample management module is responsible for sample storage and weight update; the radius calibration module calculates the distribution distance radius according to the residual and parameter posterior covariance and outputs the conservative inflation amount; the constraint generation module converts the safety function and sample weight into a solvable form; the solving and issuing module completes the optimization, maps the first-in-command to pump speed and ratio, and reads the execution data. The closed-loop evaluation updates the sample weight and radius according to the prediction and observation difference, constraint tightness, and alarm record, achieving consistent improvement of "data-model-control".

[0135] Embodiment: In a tomato greenhouse, the prediction time domain is four hours, and the control cycle is fifteen minutes. Five hundred samples are extracted from the latest two weeks of data every morning, resampled by residual stratification, and the weight is calculated; the distance radius is calibrated according to the prediction-observation difference of yesterday and the parameter posterior covariance, and the conservative inflation quantity is generated; the distribution robust probability constraint is constructed for the conductivity and liquid discharge two safety functions respectively, and the process relationship is linearized, and the nitrogen flux sequence of each partition is obtained after solving and issued. If the conductivity approaches the upper limit during operation, the system automatically shrinks the radius and increases the inflation quantity, and switches to the conservative solution; the prediction error and constraint relaxation are summarized at the end of the day, and the sample weight and radius are updated for the next day control. The above process ensures that when the scene distribution drifts and the model uncertainty exists, the fertilization setting still meets the preset safety probability and maintains stable income.

[0136] The update of the active spectral coding matrix is based on the index reflecting the constraint tension degree in the closed-loop evaluation data and the difference between the growth gain prediction and the observed growth gain to construct an update target, and the parameters of the active spectral coding matrix are adjusted according to the gradient rule.

[0137] After completing a prediction and execution, the active spectral coding matrix is updated in batches by using the closed-loop evaluation data. The active spectral coding matrix defines the combination of narrowband, polarization and illumination angle and its weight, and the update target is to improve the recognition of the observed control quantity and key physiological parameters without changing the safety and device boundaries, and to reduce the deviation between prediction and observation.

[0138] The data summary link collects prediction gain and observation gain, constraint relaxation, conductivity and liquid return record, execution deviation, parameter posterior covariance and quality mask according to the partition. Abnormal samples are excluded or weighted according to the quality mask; the parameter posterior covariance is converted into observation reliability weight for subsequent weighting.

[0139] The target construction adopts the aforementioned evaluation loss, which is summarized as a global target according to the area and reliability weight of each partition, and a change amplitude regular is added to limit the coding change of the day to avoid exposure and energy consumption exceeding the limit. The constraint set is given by the device and physiological boundaries: the upper limit of total illumination for single imaging, the range of duty cycle for each narrowband, the upper limit of switching times for polarization and illumination angle, and the upper limit of difference from historical coding.

[0140] The parameter solving adopts the aforementioned gradient update rule, and the formula is not repeated. The implementation is: fine-tuning the weight along the negative gradient direction of the evaluation loss with respect to the active spectral coding matrix in the feasible region; the gradient is obtained through the link of "prediction error-key parameter sensitivity-channel weight", and the sensitivity is estimated by the parameter posterior covariance and model residual estimated by the parameter posterior covariance and model residual. If any constraint is touched after one update, project back along the normal direction of the feasible region and shrink the step size; if it still does not meet the condition for several times, roll back to the previous version and record the reason.

[0141] To enhance the feasibility, the optimizer adopts a group update and sequential freezing strategy: first, fine-tune the weights of the red edge and near-infrared related narrow bands without changing the illumination and switching times, and then decide whether to make small step adjustments to the polarization and illumination angle according to the evaluation loss reduction amplitude; when the constraints are tight at any time of the day, the combined weight that causes the tightness is reduced first. All updates have a time stamp, only affect the encoding of the next day's inspection, and the acquisition and registration processes remain unchanged.

[0142] The system architecture includes an evaluator, a constraint manager, and an encoding optimizer. The evaluator generates a prediction error and a constraint tightness time series, and accordingly gives a partition weight; the constraint manager generates a feasible region based on device logs and daily operation boundaries, and monitors out-of-bound events; the encoding optimizer minimizes the global objective within the feasible region, and outputs a new matrix and a change summary. To ensure traceability, the system retains at least two versions of historical encoding and corresponding evaluation indicators, which can be rolled back daily.

[0143] Embodiment: In tomato greenhouse applications, the prediction and observation gain, constraint relaxation, execution deviation, and parameter posterior covariance within a four-hour rolling control window are summarized at the end of the day; first, fine-tune the weights of the red edge and near-infrared narrow bands to enhance the sensitivity to leaf surface nitrogen and leaf area index; if the conductivity tightness events are concentrated in a specific illumination angle and polarization combination, reduce its weight and increase the proportion of the top illumination combination; the whole process performs at most three rounds of step backtracking search, and if the evaluation index does not decrease enough or the constraint is approaching, the original matrix is maintained. The next day, active imaging is performed according to the updated encoding, and it is observed that the prediction error and tightness events are reduced, and the contribution of observation to control is gradually improved with iteration.

[0144] As shown in Figure 2 A greenhouse crop hyperspectral fertilization control system for implementing the greenhouse crop hyperspectral fertilization control method, the system comprises:

[0145] An acquisition and processing module is configured to acquire hyperspectral images and illumination encoding information, perform radiation correction and registration, construct a differential reflectance cube, a canopy height map, and a porosity map, and provide the differential reflectance cube, the canopy height map, and the porosity map to an inversion estimation module. The acquisition and processing module is composed of a snapshot hyperspectral camera, an electrically controllable polarizer, a narrow-band programmable light source array, a variable illumination angle mechanism, an irradiance monitoring probe, a white board reference target, and a three-dimensional imaging device. The camera, the light source, and the three-dimensional imaging device are linked through a hardware trigger and a time synchronization signal, and are connected to an edge computing unit through Ethernet or optical fiber. The edge computing unit is configured with a GPU and a high-speed NVMe storage, performs dark field subtraction, response linearization, intensity normalization, and cross-band registration, and outputs the differential reflectance cube, the canopy height map, and the porosity map; a temperature control and stray light shielding structure is used to stabilize the radiation response. The data is published to the inversion estimation module through gigabit Ethernet.

[0146] The inversion estimation module, based on polarization spectral unmixing and radiative transfer inversion, determines leaf surface nitrogen, leaf area index, leaf tilt angle distribution, and cluster index, generates the posterior covariance of the parameters, and provides these parameters to the identification modeling module and the control optimization module. The inversion estimation module consists of an edge server or rack-mounted workstation equipped with a GPU and large-capacity memory, and has a built-in high-reliability solid-state array for storing the spectral library and calibration files. This module receives data from the acquisition and processing module via industrial Ethernet, utilizes the GPU to accelerate polarization spectral unmixing and two-layer radiative transfer forward calculations, and runs continuously under local redundant power supply and UPS protection. The results are output in shared memory or message queue format, including leaf surface nitrogen, leaf area index, leaf tilt angle distribution, cluster index, and posterior covariance of the parameters, for downstream applications to subscribe to in real time.

[0147] The identification and modeling module applies pseudo-random binary perturbations within the control cycle. Combining the root zone nitrogen budget model with Gaussian process regression, it determines the growth gain function and marginal response based on leaf surface nitrogen, leaf area index, leaf tilt angle distribution, cluster index, and posterior covariance of parameters. It then generates a scene distribution and provides the growth gain function, marginal response, and scene distribution to the control optimization module. The identification and modeling module consists of a greenhouse fertilization and irrigation control cabinet, a variable frequency metering pump, a proportional valve, a flow meter, a pressure sensor, a return liquid metering and sampling unit, online EC and pH sensors, a substrate moisture content sensor, a photosynthetically active radiation sensor, and an industrial controller. The industrial controller injects pseudo-random binary perturbations as planned, records actual flow rates and concentrations, aligns them with environmental sensor times, and sends the data via Ethernet to an embedded computing unit. The embedded unit runs root zone nitrogen budget estimation and Gaussian process regression, generating the growth gain function, marginal response, and scene distribution in real time. All execution and alarm events are written to both local logs and a central database.

[0148] The control optimization module is used to establish a bibliometric probabilistic constraint model for predictive control based on the growth gain function, marginal response, scene distribution, and posterior covariance of parameters. It generates and executes zonal fertilization settings, receives closed-loop evaluation data, and updates the active spectral coding matrix accordingly. The control optimization module consists of a real-time industrial computer or PLC with an HMI panel, communicating with the fertilization and irrigation control cabinet via Modbus or OPCUA to issue zonal fertilization settings to the variable frequency pumps and valves. This module loads the growth gain function, marginal response, scene distribution, and posterior covariance of parameters, solves the model predictive control online, and performs closed-loop correction based on the read-back flow rate and EC. The encoding / decoding unit, directly connected to the lighting controller, is responsible for updating the active spectral coding matrix and generating a narrowband, polarization, and lighting angle command table for the next day's inspection. The safety link includes an emergency stop circuit, solenoid valve cutoff, and over-limit interlocking to ensure that the strategy degrades to a conservative approach when constraints are approaching and that closed-loop evaluation data is fully recorded.

[0149] Those skilled in the art will appreciate that embodiments of the present application can be readily used as a method, a system or a computer program product. Accordingly, the present application can take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects.

[0150] The embodiments of the present application are only intended to illustrate the present application, and not intended to limit the present application. Various modifications and changes can be made by those skilled in the art without departing from the spirit and scope of the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should be included in the scope of the claims of the present application.

Claims

1. A method for controlling hyperspectral fertilization of greenhouse crops, characterized in that, Includes the following steps: Acquire hyperspectral images and illumination coding information, perform radiometric correction and registration, and construct differential reflectance cubes, canopy height maps, and porosity maps; Based on polarization spectral unmixing and radiative transfer inversion, the leaf surface nitrogen, leaf area index, leaf tilt angle distribution and cluster index are determined, and the posterior covariance of the parameters is generated. The polarization spectrum unmixing is performed by nonnegative matrix decomposition and polarization difference constraint to separate the blade endmembers. The radiative transfer inversion is solved by a joint two-layer model of the blade radiative transfer model and the canopy radiative transfer model. The posterior covariance of the parameters is estimated based on the linearization sensitivity of the radiative transfer inversion and the observation noise to obtain a measure of the uncertainty of each inversion parameter. A pseudo-random binary sequence perturbation is applied within the control period. The growth gain function and marginal response are determined by combining the root region nitrogen budget model with Gaussian process regression, and the scene distribution is generated. The temporal model of pseudo-random binary sequence perturbation is designed under amplitude constraints and daily total input constraints, and optimized with the goal of improving parameter identifiability. Gaussian process regression introduces pseudo-random binary sequences as instrumental variables. By constraining the correlation between the instrumental variables and the regression residuals to zero, the correlation between the input and the external environment is eliminated, thereby determining the growth gain function and the marginal response. Based on the growth gain function, scene distribution, and posterior covariance of parameters, a bibliometric probability constraint model is established for predictive control, outputting zonal fertilization settings, and updating the active spectral coding matrix based on closed-loop evaluation data.

2. The method according to claim 1, characterized in that, The illumination coding information consists of a combination of narrowband, polarization, and illumination angle. The differential reflection cube is constructed by differential normalization of the reference illumination intensity and the unilluminated reference frame. The canopy height map is reconstructed from 3D imaging. The porosity map is calculated by the area ratio of the vegetation mask in the local window.

3. The method according to claim 1, characterized in that, The root zone nitrogen budget model, combined with nonlinear Kalman filtering, performs joint state estimation and covariance update of the root zone inorganic nitrogen inventory and nitrogen uptake rate.

4. The method according to claim 1, characterized in that, The predictive control model of the distributed bar probabilistic constraint uses a set of uncertainties built based on the distribution distance metric to form probabilistic constraints in order to cope with the uncertainty of the scene distribution.

5. The method according to claim 1, characterized in that, The update of the active spectral coding matrix is ​​based on the indicators reflecting the degree of constraint tension in the closed-loop evaluation data and the difference between the predicted growth gain and the observed growth gain. The update target is constructed and the parameters of the active spectral coding matrix are adjusted according to the gradient rule.

6. A greenhouse crop hyperspectral fertilization control system, used to implement the greenhouse crop hyperspectral fertilization control method according to any one of claims 1 to 5, characterized in that, The system includes: The acquisition and processing module is used to acquire hyperspectral images and illumination coding information, perform radiometric correction and registration, construct differential reflectance cubes, canopy height maps and porosity maps, and provide the differential reflectance cubes, canopy height maps and porosity maps to the inversion estimation module; The inversion estimation module is used to determine leaf surface nitrogen, leaf area index, leaf tilt angle distribution, and cluster index based on polarization spectrum unmixing and radiative transfer inversion, generate parameter posterior covariance, and provide leaf surface nitrogen, leaf area index, leaf tilt angle distribution, cluster index, and parameter posterior covariance to the identification modeling module and control optimization module; the time series of pseudo-random binary sequence perturbation is designed under amplitude constraints and daily total input constraints, and optimized with the goal of improving parameter identifiability; Gaussian process regression introduces pseudo-random binary sequence as instrumental variable, and eliminates the correlation between input and external environment by constraining the correlation between instrumental variable and regression residual to zero, thereby determining the growth gain function and marginal response; The identification and modeling module applies pseudo-random binary sequence perturbations within the control cycle. Combining the root zone nitrogen budget model with Gaussian process regression, it determines the growth gain function and marginal response based on leaf surface nitrogen, leaf area index, leaf tilt angle distribution, cluster index, and parameter posterior covariance, generating scene distributions. The growth gain function, marginal response, and scene distribution are then provided to the control optimization module. The temporal series of the pseudo-random binary sequence perturbations is designed under amplitude and daily total input constraints and optimized to improve parameter identifiability. Gaussian process regression introduces the pseudo-random binary sequence as an instrumental variable. By constraining the correlation between the instrumental variable and the regression residual to zero, the correlation between the input and the external environment is eliminated, thereby determining the growth gain function and marginal response. The control optimization module is used to establish a predictive control model based on the growth gain function, marginal response, scene distribution, and parameter posterior covariance, generate zonal fertilization settings and issue execution orders, receive closed-loop evaluation data and update the active spectral coding matrix accordingly.

Images