Multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing

Abstract

The invention provides a multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing, which relates to the field of crop yield estimation, and comprises the following steps: S1, multi-source data acquisition: respectively acquiring remote sensing data consisting of optical images and radar images of crops in a region, meteorological data and corresponding auxiliary data, the data is preprocessed; s2, key feature extraction and quantification: biophysical feature inversion, structural feature inversion and environmental stress feature quantification are adopted to extract biophysical features reflecting crop growth states and environmental stress from the preprocessed data; and S3, constructing a dynamic graph neural network of a coupling physical mechanism to realize dynamic yield estimation of crops. According to the method, the yield of crops in the whole growth period can be dynamically monitored and predicted through the multi-source data of the satellite image, and the method has high precision and high timeliness.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of crop yield estimation, and particularly relates to a multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing. BACKGROUND

[0002] As a core link of agricultural production management, food security early warning and agricultural product market regulation, crop yield estimation is of great significance to guide field management and decision-making of agricultural production operators, optimize field management, reduce "supply-demand mismatch" of the industry chain, improve industry efficiency, promote variety selection and iteration of planting technology. At present, the yield estimation of crops mainly relies on long-term statistical investigation and experience model. However, the long-term statistical investigation, that is, the method of manually sampling and drawing statistical report forms to estimate the yield of crops in the region, has the disadvantages of poor timeliness, insufficient spatial coverage, high labor cost and the like, and it is difficult to meet the needs of modern agricultural efficient, low-cost and fine management. At the same time, the traditional experience model relies on a large number of and complex experience data, and the data is large and has many redundant data, and is greatly affected by external environmental factors (for example, if the crops in the field have great differences in growth due to diseases, pests, uneven water and fertilizer, etc., the data sample size needs to be increased), resulting in slow data updating and large estimation error. SUMMARY

[0003] In view of the problems existing in the prior art, the purpose of the present application is to provide a multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing. The method is based on multi-source data of satellite remote sensing, fuses dynamic graph neural network and growth law of crops, realizes dynamic monitoring and prediction of crop yield in the whole growth period, has high precision and high timeliness, and effectively avoids the problems of large estimation error and poor timeliness caused by sampling cycle, external environmental factors and human influence factors.

[0004] The purpose of the present application is realized by the following technical scheme: A multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing, comprising: Step S1, multi-source data acquisition: acquiring remote sensing data composed of optical images and radar images of crops in the region, meteorological data and corresponding auxiliary data, and pre-processing the data; Step S2, key feature extraction and quantification: biological physical feature inversion, structural feature inversion and environmental stress feature quantification are used to extract biological physical features reflecting the growth state and environmental stress of crops from the pre-processed data; Step S3, constructing a dynamic graph neural network coupled with physical mechanism: through dynamic graph structure modeling of space-time correlation, embedding the growth physical mechanism of crops, capturing the lagging effect of environmental factors, realizing dynamic yield estimation of crops.

[0005] Based on the further optimization of the above scheme, the optical image is obtained by Landsat or Sentinel-2, the radar image is obtained by Sentinel-1, the optical image and the radar image cover the whole growth period of the crop; the meteorological data includes time series data of temperature, precipitation, humidity and the like within the regional scale; the auxiliary data includes a digital elevation model (DEM), a land use map, soil texture data, crop variety parameters and the like.

[0006] Based on the further optimization of the above scheme, the data preprocessing in the step S1 includes optical and radar data registration and meteorological data scale registration. Optical and radar data registration: mutual information maximization method is adopted to realize data space alignment and eliminate space-time deviation:

[0007] In the formula: l o The optical image is represented by I; l r The radar image is represented by R; The registration offset is represented by D; H The information entropy is represented by H, which is used to measure the uncertainty of the gray scale distribution in the image:

[0008] In the formula: l The single image is represented by I; G The total number of gray values is represented by N; p(g) The probability of the occurrence of the gray value g is represented by P; The radar image is represented by R; l r Offset After the offset, the joint entropy of the optical-radar pixel pair is calculated; Meteorological data scale registration: based on the graph convolution of physical constraints, the low-resolution meteorological data is down-scaled to the field scale:

[0009] In the formula: X i,t The field is represented by F; i The down-scaled meteorological data at time t is represented by M; t The low-resolution meteorological data at the initial time t0 is represented by M0; X(t 0 ) The adjacency matrix constructed based on the DEM and the land use map is represented by A; t 0; A The unit matrix is represented by I; DEM The adjacency matrix constructed based on the DEM and the land use map is represented by A; I The unit matrix is represented by I; The a degree matrix; W q denotes a weight matrix, b q denotes a bias vector.

[0010] Based on further optimization of the above scheme, the biophysical characteristic inversion in the step S2 adopts a leaf optical model PROSPECT and a canopy radiation transfer model SAIL coupled inversion, specifically: taking the reflectance spectrum of a satellite remote sensing image , a solar zenith angle , an observation zenith angle , a relative azimuth angle Δφ as input data, and taking a leaf structure parameter N , a soil background reflectance as a prior parameter; adopting the PROSPECT model to simulate the reflectance and transmittance of a single leaf:

[0011] In the formula: C ab denotes a chlorophyll content; C w denotes a leaf water content; C m denotes a dry matter content; adopting the SAIL model to simulate the reflectance of a canopy in a direction:

[0012] In the formula: denotes a leaf area index; Then, the parameters are inverted by minimizing the residual error between the observed spectrum and the simulated spectrum, and the objective function is:

[0013] In the formula: denotes an observed spectrum; denotes a simulated spectrum, which is obtained by coupling the PROSPECT model and the PROSPECT model:

[0014] The Levenberg-Marquardt algorithm is adopted, the range constraint of the prior parameter is combined, and the chlorophyll content C ab , the leaf area index , the dry matter content C m and the like are iteratively solved to minimize the residual error; and the biomass is obtainedSW :

[0015] wherein: LMA represents leaf weight, which is derived from dry matter content obtained by PROSPECT model inversion C m and leaf thickness.

[0016] Further optimization based on the above scheme, the structural characteristics inversion in step S2 is inverted by RVoG model to obtain height and structural parameters, specifically: Convert the collected radar echo to backscattering coefficient , and pre-process by denoising and geometric correction; Construct the RVoG model:

[0017] wherein: represents the backscattering coefficient of the ground; represents the volume scattering coefficient of the crop (reflecting the volume scattering ability of the vertical structure of the crop, related to the stem density and distribution); represents the radar wavelength; represents the radar incidence angle; represents the crop height; Invert the crop height and vertical structure parameters:

[0018] Solve the crop height by nonlinear optimization and the volume scattering coefficient of the crop .

[0019] Further optimization based on the above scheme, the environmental stress characteristic quantification in step S2 includes water stress index SWD and cumulative vapor pressure deficit VPD: Obtain the effective soil moisture SM av :

[0020] wherein: SM cu represents the current soil moisture (volume water content); SM wi represents the soil moisture at the wilting point; Obtain the crop evapotranspiration requirement by Penman-Monteith model ET de :

[0021] In the formula: K c represents the crop coefficient; represents the rate of change of saturated water vapor pressure with temperature (reflecting the effect of temperature on saturated water vapor pressure); R n represents net radiation (representing the energy driving of evapotranspiration); RG represents soil heat flux (i.e., the amount of heat exchange between the soil surface and the deep layer); represents the conversion coefficient of dry-wet bulb temperature difference and water vapor pressure (related to atmospheric pressure); v 2 represents the wind speed at two meters high; T represents the air temperature; e s represents the saturated water vapor pressure (i.e., the maximum water vapor pressure that the air can contain at the current temperature); e a represents the actual water vapor pressure;

[0022] In the formula: RH represents the relative humidity; P represents the atmospheric pressure; Then the water stress index SWD is:

[0023] In the formula: RD x represents the crop root zone depth; Cumulative vapor pressure deficit VPD:

[0024] In the formula: RJ(t) represents t the rainfall of the period; VPD(t) represents t the vapor pressure deficit at the moment; The vapor pressure deficit in the crop growth period is calculated by time period to obtain the cumulative vapor pressure deficit.

[0025] Based on the further optimization of the above scheme, the step S3 is specifically: First, a field-level dynamic graph neural network structure evolving with the growth period is constructed:

[0026] In the formula: V t represents a node set, which has a plurality of nodes v i,t ; each node v i,t corresponds to a field i at the momentFeature vector of t s i,t The feature vector includes biophysical features (chlorophyll content C ab , leaf area index , biomass SW, etc.), structural features (crop height , volume scattering coefficient of crops , etc.), and environmental stress features (water stress index SWD, cumulative vapor pressure deficit VPD, etc.); E t Indicates the edge set, and the edge weight w i,j,t Quantify the spatio-temporal correlation between the field block and i and the real-time meteorological similarity, specifically: j

[0027] In the formula: Indicates the soil texture gradient of the field block i , j (representing soil heterogeneity); Q i,t , Q j,t Indicates the meteorological vector of the field block i , j at time t; Indicates a hyperparameter for regulating the influence strength of the soil texture gradient difference on the edge weight; The node feature aggregates spatio-temporal information through graph convolution, and the propagation equation is:

[0028] In the formula: Indicates the hidden feature of the l layer, and the initial layer ; Indicates the weight matrix of the l layer graph convolution;

[0029] Wherein, A t Indicates the original adjacency matrix of the field block at time t , that is, ; Z t Indicates the historical feature weighting value output by the lag effect module, which is designed to capture the cumulative effect of environmental factors by using a learnable time-varying convolution kernel to capture the lag effect:

[0030] ​In the formula: S(t-k) represents the lag k period field characteristics; W k represents the first k period convolution kernel for extracting the local correlation of historical features; represents the lag weight at time t;

[0031] In the formula: U represents a row vector (dimension consistent with the tanh output) for mapping the nonlinear transformed features to scalar values; Y represents a transformation matrix for mapping the input to the input dimension of tanh; E t represents the meteorological stress factor at time t ; h t-1 represents the GRU hidden state at time t-1 ; represents the meteorological stress factor E t and the GRU hidden state h t-1 spliced into a vector; lag effect output Z t as a supplementary input of the dynamic graph propagation equation, so that the node feature aggregation contains real-time information and historical cumulative influence at the same time, solving the limitation of traditional graph models relying only on current time features; Physical mechanism constraints:

[0032] In the formula: represents the environmental stress and meteorological driving feature vector at time t; NPP represents net primary productivity; represents the distribution coefficient; HI represents the harvest index; f i,0 represents the initial value of the distribution coefficient (i.e., the carbon distribution proportion of crop organs i in the early growth period), f i,max represents the maximum value of the distribution coefficient (i.e., the upper limit of the carbon distribution proportion of crop organs i in the late growth period), represents the rate parameter of the distribution phase; GDD represents the growing degree day, i.e., the degree day number of cumulative temperature exceeding the base temperature, GDD s represents the starting GDD at which crop organs i begin to dominate carbon distribution; GPP represents canopy primary productivity, represents the corrected maintenance respiration, Rg represents the base maintenance respiration; LAI represents the leaf area index, which is obtained by inversion of the PROSPECT model coupled with the SAIL model, V cmax represents the maximum carboxylation rate (which is proportional to the chlorophyll content C ab positively correlated), represents the CO 2concentration of the intercellular space of the leaf, represents CO 2the compensation point (which increases with the increase of temperature), represents CO 2the Michaelis constant of Rubisco, O 2represents the oxygen concentration in the atmosphere, represents O 2the Michaelis constant of Rubisco, J() represents the electron transport efficiency, and PAR represents the photosynthetically active radiation, g s () represents the stomatal conductance, represents the carbon dioxide concentration in the atmosphere, R d represents the crop dark respiration rate (i.e. the respiration digestion under no light); k m represents the crop base coefficient, Q 10 represents the temperature coefficient, and T represents the real-time temperature, T ref represents the reference temperature, represents the water sensitivity coefficient; represents the assimilate conversion efficiency; Physical constraints and closed-loop training of dynamic graphs: Input: dynamic graph node features: s i,t Fusion of remote sensing inversion features, environmental stress features and meteorological downscaling data; State initialization: initialization by the RVoG model inversion of the biomass, to obtain SW i,0 as the initial conditions of the differential equation; Prediction of the biomass at time t by dynamic graph aggregation of spatio-temporal features ; Physical bias loss L phys :

[0033] In the formula: represents the differential equation obtained by physical mechanism constraint, i.e.

[0034] The final total loss function is:

[0035] In the formula: Y represents the model prediction of the yield at the moment t; t Y represents the measured yield at the moment t. Y represents the measured yield at the moment t.

[0036] The technical effects of the technical scheme of the present application are as follows: The present application inverses biophysical characteristics by PROSPECT model and SAIL model, inverses structural characteristics by RVoG model, and quantifies environmental stress characteristics, so as to comprehensively depict the growth state of field crops, avoiding the limitation of single data; then, the growth process of crops is embedded into a dynamic graph neural network, the dependence on redundant empirical data is effectively reduced through the constraint of physical deviation loss, the estimation deviation caused by external environmental factors is reduced, and the accuracy of crop yield estimation is improved. At the same time, the present application extracts the historical cumulative influence of environmental factors through time-varying convolution kernel and GRU module, effectively solves the defects of traditional models relying on "current time data", fully considers the hysteresis of yield estimation, so as to adapt to the continuity law of crop growth. In addition, the present application relies on existing satellite remote sensing technology, through multi-source data complementation (i.e. optical image, radar image, etc.) and quantification of environmental stress characteristics, reduces the influence of single data on weather, land heterogeneity (such as soil texture difference), etc., improves the stability of yield estimation under complex environmental light, and efficiently and accurately realizes dynamic yield estimation of crops. BRIEF DESCRIPTION OF DRAWINGS

[0037] Figure 1 The structural block diagram of the dynamic yield estimation of crops in the embodiment of the present application is shown. DETAILED DESCRIPTION

[0038] The technical scheme in the embodiment of the present application will be described clearly and completely below, and in the following description, specific details such as specific system structure, technology, etc. are proposed for the purpose of explanation rather than limitation, so as to thoroughly understand the embodiment of the present application.

[0039] Embodiment 1: A multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing, comprising: Step S1: Multi-source data acquisition: Collect remote sensing data, meteorological data, and corresponding auxiliary data, consisting of optical and radar images of crops in the region. Optical images are acquired through Landsat or Sentinel-2, and radar images are acquired through Sentinel-1. The optical and radar images cover the entire growth period of crops. Meteorological data includes time-series data such as temperature, precipitation, and humidity at the regional scale (which can be acquired through meteorological monitoring stations and corresponding sensors). Auxiliary data includes digital elevation models (DEM), land use maps, soil texture data, crop variety parameters, etc.

[0040] The data is preprocessed, including optical and radar data registration and meteorological data scale registration. Optical and radar data registration: A mutual information maximization method is used to achieve spatial alignment of data and eliminate spatiotemporal bias.

[0041] In the formula: l o Represents optical images; l r Represents radar imagery; Indicates the registration offset; H Information entropy is used to measure the uncertainty of gray-level distribution in an image.

[0042] In the formula: l Represents a single image; G This represents the total number of grayscale values ​​(e.g., in an 8-bit image). G =256); p(g) This represents the probability of the grayscale value g appearing. Indicates when radar image l r Offset Then, the joint entropy of the optical-radar pixel pairs is calculated; Meteorological data scaling: Physically constrained graph convolution to downscale low-resolution meteorological data to the field scale.

[0043] In the formula: X i,t Indicates field i At any moment t Downscaled meteorological data; X(t 0 ) Indicates the initial time. t Low-resolution meteorological data at 0; A Indicates based on DEMand a land use map; I denotes an identity matrix; denotes a degree matrix of W q denotes a weight matrix, b q denotes a bias vector.

[0044] Step S2, key feature extraction and quantification: biological and physical feature inversion, structural feature inversion and environmental stress feature quantification are used to extract biological and physical features reflecting crop growth status and environmental stress from preprocessed data; Biophysical feature inversion uses the coupling of leaf optical model PROSPECT and canopy radiation transfer model SAIL, specifically: The reflectance spectrum of satellite remote sensing image (multi-spectral / hyperspectral, covering 400-2500nm), solar zenith angle , observation zenith angle , relative azimuth angle Δφ (computed from satellite orbit and observation time) as input data, leaf structure parameters N (adjusted according to crop type, generally 1.2-2.5), soil background reflectivity (obtained by bare ground pixel inversion or a large number of empirical values) as prior parameters; The reflectivity and transmittance of single leaf are simulated by PROSPECT model:

[0045] In the formula: C ab denotes chlorophyll content; C w denotes leaf water content; C m denotes dry matter content; The reflectivity of the canopy in the direction is simulated by SAIL model:

[0046] In the formula: denotes leaf area index; Then, the parameters are inverted by minimizing the residual error between the observed spectrum and the simulated spectrum, and the objective function is:

[0047] In the formula: denotes the observed spectrum; denotes the simulated spectrum, which is obtained by coupling the PROSPECT model and the PROSPECT model:

[0048] The sampling Levenberg-Marquardt algorithm, combined with the range constraint (pre-set) of prior parameters, iteratively solves the minimum residual chlorophyll content C ab , leaf area index , dry matter content C m , etc. and obtain biomass SW :

[0049] In the formula: LMA represents the leaf weight, which is the dry matter content obtained by the inversion of the PROSPECT model C m and the leaf thickness is derived.

[0050] The structural feature inversion is to invert the height and structural parameters through the RVoG model, specifically: The collected radar echo is converted into the backscattering coefficient , and is preprocessed by noise removal (Lee filter is used to suppress speckle noise) and geometric correction (converted into the ground object coordinate system, and the crop planting area pixel is extracted); The RVoG model is constructed:

[0051] In the formula: represents the surface backscattering coefficient (obtained by inversion of the bare soil pixel); represents the volume scattering coefficient of the crop (reflecting the volume scattering ability of the vertical structure of the crop, related to the stem density and distribution); represents the radar wavelength; represents the radar incidence angle; represents the crop height; The crop height and the vertical structure parameters are obtained by inversion:

[0052] The crop height and the volume scattering coefficient of the crop are obtained by nonlinear optimization .

[0053] The environmental stress quantification includes the water stress index SWD and the cumulative vapor pressure deficit VPD: The effective soil water is obtained SM av :

[0054] ​In the formula: SM cu represents the current soil moisture (volume water content); SM wi represents the soil moisture at the wilting point (obtained according to empirical values and soil texture standards, for example: clay 0.1, sandy soil 0.03); Penman-Monteith model is used to obtain crop evapotranspiration demand ET de :

[0055] In the formula: K c represents the crop coefficient (adjusted according to the crop growth period); represents the rate of change of saturated water vapor pressure with temperature (reflecting the effect of temperature on saturated water vapor pressure); R n represents the net radiation (representing the energy driving of evapotranspiration); RG represents the soil heat flux (i.e. the heat exchange between the soil surface and the deep layer); represents the conversion coefficient of dry-wet bulb temperature difference and water vapor pressure (related to atmospheric pressure); v 2 represents the wind speed at two meters high; T represents the air temperature; e s represents the saturated water vapor pressure (i.e. the maximum water vapor pressure that the air can contain at the current temperature); e a represents the actual water vapor pressure (obtained by relative humidity);

[0056] In the formula: RH represents the relative humidity; P represents the atmospheric pressure; Then the water stress index SWD is:

[0057] In the formula: RD x represents the crop root zone depth; Cumulative vapor pressure deficit VPD:

[0058] In the formula: RJ(t) represents t the rainfall of the period; VPD(t) represents t the vapor pressure deficit at the moment; The vapor pressure deficit during the crop growth period is calculated by time period to obtain the cumulative vapor pressure deficit.

[0059] Step S3, constructing a dynamic graph neural network coupled with physical mechanisms: modeling the spatio-temporal correlation through a dynamic graph structure, embedding the growth physical mechanisms of crops, capturing the lag effects of environmental factors, and realizing the dynamic yield estimation of crops; specifically: First, construct a field-level dynamic graph neural network structure evolving with growth period:

[0060] In the formula: V t denotes a node set, which has several nodes v i,t consisting of; each node v i,t corresponds to a field i At time t , the feature vector s i,t , the feature vector includes biophysical characteristics (chlorophyll content C ab , leaf area index , biomass SW, etc.), structural characteristics (crop height , crop volume scattering coefficient , etc.), and environmental stress characteristics (water stress index SWD, cumulative vapor pressure deficit VPD, etc.); E t denotes an edge set, and the edge weight w i,j,t quantifies the spatio-temporal correlation between the fields i and j , and fuses the soil texture and real-time meteorological similarity, specifically:

[0061] In the formula: denotes the soil texture gradient of the fields i , j (representing soil heterogeneity); Q i,t , Q j,t denotes the meteorological vector of the fields i , j at time t; denotes a hyperparameter for regulating the influence strength of the soil texture gradient difference on the edge weight; The node features are aggregated and propagated through graph convolution, and the propagation equation is:

[0062] In the formula: denotes the hidden features of the l layer, and the initial layer ; Indicates the first l The weight matrix of layer graph convolution;

[0063] in, A t Indicates time t The original adjacency matrix of the lower and upper plots, i.e. ; Z t This represents the weighted value of historical features output by the lag effect module. To capture the cumulative effects of environmental factors, a learnable time-varying convolution kernel is designed.

[0064] In the formula: S(t-k) Indicates lag k Field characteristics during the period (e.g.: t-k (SWD at any moment) W k Indicates the first k The convolutional kernels of the period are used to extract local correlations of historical features; The lag weight at time t;

[0065] In the formula: U represents a row vector (with the same dimension as the tanh output), used to map the nonlinearly transformed features to scalar values; Y represents the transformation matrix, which transforms the input... The input dimension mapped to tanh; E t Indicates time t Meteorological stress factors; h t-1 Indicates time t-1 The hidden state of GRU; This indicates the meteorological stress factors. E t With GRU hidden state h t-1 Concatenate them into a single vector; Lag effect output Z t As a supplementary input to the dynamic graph propagation equation, it enables the aggregation of node features to simultaneously include real-time information and historical cumulative influence, thus overcoming the limitation of traditional graph models that only rely on features at the current moment. Physical mechanism constraints:

[0066] In the formula: Environmental stress and weather driven feature vector at time t; NPP represents net primary productivity; Distribution coefficient; HI represents harvest index (adjusted according to variety adoption number and growth period dynamics); f i,0 Distribution coefficient initial value (i.e. carbon distribution proportion of crop organs in early growth period), i f i,max Distribution coefficient maximum value (i.e. upper limit of carbon distribution proportion of crop organs in late growth period), i Rate parameter of distribution stage (set according to crop variety and actual situation, such as wheat grain distribution rate), GDD represents growing degree day, i.e. degree day number of cumulative temperature exceeding base temperature, GDD s Starting GDD of crop organ i starting to dominate carbon distribution; GPP represents canopy primary productivity, Corrected maintenance respiration, R g Base maintenance respiration; LAI represents leaf area index, obtained by coupling PROSPECT model and SAIL model, V cmax Maximum carboxylation rate (positively correlated with chlorophyll content C ab ), V cmax =k C ab +b (k, b represent crop-specific parameters), Concentration of CO 2intercellular space of leaf (affected by stomatal conductance and photosynthetic consumption, can be obtained by remote sensing inversion of hyperspectral red edge range characteristics), Compensation point of CO 2, Michaelis constant of CO 2, O 2represents oxygen concentration in the atmosphere, Michaelis constant of O 2, J() Electron transport efficiency, PAR represents photosynthetic active radiation, g s () Stomatal conductance, Carbon dioxide concentration in the atmosphere, R d Crop dark respiration rate (i.e. respiration digestion under no light); k m ​​Crop base coefficient, Q 10 Temperature coefficient, T Real-time temperature, T ref Reference temperature, Moisture sensitivity coefficient; Assimilate conversion efficiency (determined according to different varieties of crops); LAI: Using the PROSPECT model, input biochemical parameters of leaves (such as chlorophyll content C ab , leaf structure parameters N , water content C w , etc.), simulate the reflectivity and transmittance of leaves at different wavelengths:

[0067] Take the reflectivity and transmittance output by the PROSPECT model as input, combine with canopy parameters (LAI, leaf angle distribution LAD, sun / observation angle, etc.), simulate canopy reflectivity:

[0068] According to the real canopy reflectivity obtained from remote sensing images, calculate the observation error between the simulated canopy reflectivity :

[0069] Adjust LAI by iteration to minimize Error until the simulated value and the observed value are optimally matched, and the LAI at this time is the inversion result;

[0070] In the formula: m, n respectively represent the corresponding crop coefficient;

[0071] In the formula: A represents net photosynthesis; Physical constraints and closed-loop training of dynamic graph: Input: dynamic graph node features: s i,t Fuse remote sensing inversion features, environmental stress features and meteorological downscaling data; State initialization: initialize biomass through RVoG model inversion, obtain SW i,0 as the initial condition of the differential equation; Predicting biomass at time t by dynamic graph aggregation of spatio-temporal features Physical bias loss L phys

[0072] wherein: is the differential equation obtained by constraining through physical mechanisms, i.e.:

[0073] The final total loss function is:

[0074] wherein: is the yield predicted by the model at time t (obtained through predicted biomass t and harvest index); is the measured yield (obtained through field sampling or statistics).

[0075] Embodiment 2: As another preferred embodiment of the present application, on the basis of the scheme of embodiment 1, the dynamic yield estimation process in step S3 is specifically: Take t=0 as the starting point of the growth period, input the initial features, including: soil texture at sowing period, initial biomass, initial values of meteorology, etc.; and construct the initial graph: initialize the edge weight based on soil texture similarity w i,j,0 Time sequence iteration: Take t =1, 2, …, t T as the total length of the growth period, process the historical features through the time-varying convolution kernel, and output the lagging features Z t update the edge weight based on real-time meteorology w i,j,t obtain the adjacency matrix A t complete the dynamic graph update; input Z t and the current node features s i,t into the dynamic graph, update the hidden layer through the propagation equation ; then, calculate the biomass:

[0076] wherein: is the output layer mapping function; ​​​​and the predicted value is corrected by the differential equation constraint; Finally, the biomass The yield prediction at time t is obtained with the crop harvest index: .

Claims

1. A method for multi-source data processing and dynamic crop yield estimation based on satellite remote sensing, characterized in that: include: Step S1: Multi-source data acquisition: Collect remote sensing data composed of optical images and radar images of crops in the area, meteorological data, and corresponding auxiliary data, and preprocess the data. Step S2, Key Feature Extraction and Quantification: Biophysical feature inversion, structural feature inversion, and environmental stress feature quantification are used to extract biophysical features reflecting crop growth status and environmental stress from the preprocessed data; Step S3: Construct a dynamic graph neural network with coupled physical mechanisms: Model spatiotemporal correlations through dynamic graph structures, embed the physical mechanisms of crop growth, capture the lag effects of environmental factors, and realize dynamic yield estimation of crops.

2. The method for multi-source data processing and dynamic crop yield estimation based on satellite remote sensing according to claim 1, characterized in that: The optical images are acquired using Landsat or Sentinel-2, and the radar images are acquired using Sentinel-1. The optical and radar images cover the entire growth period of the crop. The meteorological data includes time-series data of temperature, precipitation, and humidity at the regional scale. The auxiliary data include digital elevation models, land use maps, soil texture data, and crop variety parameters.

3. A method for multi-source data processing and dynamic crop yield estimation based on satellite remote sensing according to claim 1 or 2, characterized in that: The data preprocessing in step S1 includes optical and radar data registration and meteorological data scale registration. Optical and radar data registration: A mutual information maximization method is used to achieve spatial alignment of data and eliminate spatiotemporal bias. In the formula: l o Represents optical images; l r Represents radar imagery; Indicates the registration offset; H Information entropy is used to measure the uncertainty of gray-level distribution in an image. In the formula: l Represents a single image; G Indicates the total number of grayscale values; p(g) This represents the probability of the grayscale value g appearing. Indicates when radar image l r Offset Then, the joint entropy of the optical-radar pixel pairs is calculated; Meteorological data scaling: Physically constrained graph convolution to downscale low-resolution meteorological data to the field scale. In the formula: X i,t Indicates field i At any moment t Downscaled meteorological data; X(t 0 ) Indicates the initial time. t Low-resolution meteorological data at 0; A Indicates based on DEM Adjacency matrix constructed from land use map; I Represents the identity matrix; express The degree matrix; W q Represents the weight matrix. b q This represents the bias vector.

4. A method for multi-source data processing and dynamic crop yield estimation based on satellite remote sensing according to claim 1 or 3, characterized in that: In step S2, the biophysical feature inversion employs a coupled inversion method using the leaf optical model PROSPECT and the canopy radiative transfer model SAIL, specifically as follows: using the reflectance spectrum of satellite remote sensing images Solar zenith angle Observing the zenith angle Relative azimuth Δφ As input data, blade structural parameters N Soil background reflectance As a priori parameter; The reflectivity and transmittance of a single blade were simulated using the PROSPECT model: In the formula: C ab Indicates chlorophyll content; C w Indicates the water content of the leaves; C m Indicates dry matter content; The SAIL model was used to simulate the reflectivity in the canopy direction: In the formula: Indicates leaf area index; Then, by minimizing the residual inversion parameters between the observed and simulated spectra, the objective function is: In the formula: Indicates the observed spectrum; The simulated spectrum is obtained by coupling the PROSPECT model with another PROSPECT model: The Levenberg-Marquardt algorithm is sampled and, combined with the range constraints of prior parameters, iteratively solves for the chlorophyll content that minimizes the residual. C ab Leaf area index Dry matter content C m ; And obtain biomass SW : In the formula: LMA This represents the leaf weight, which is the dry matter content obtained by inversion using the PROSPECT model. C m The thickness was derived from the blade thickness.

5. The method for multi-source data processing and dynamic crop yield estimation based on satellite remote sensing according to claim 3, characterized in that: In step S2, the structural feature inversion involves retrieving height and structural parameters using the RVOG model, specifically as follows: The collected radar echoes are converted into backscattering coefficients. And it has undergone preprocessing including noise reduction and geometric correction; Building the RVOG model: In the formula: Indicates the backscattering coefficient of the Earth's surface; Represents the volume scattering coefficient of the crop; Indicates the radar wavelength; Indicates the radar incident angle; Indicates crop height; The crop height and vertical structure parameters were obtained through inversion: Crop height is solved by nonlinear optimization. Volume scattering coefficient of crops .

6. The method for multi-source data processing and dynamic crop yield estimation based on satellite remote sensing according to claim 1, characterized in that: The environmental stress characteristics quantification in step S2 includes the moisture stress index SWD and the cumulative vapor pressure deficit VPD: Obtain effective soil moisture SM av : In the formula: SM cu Indicates the current soil moisture; SM wi Indicates soil moisture at the point of wilting; Crop evapotranspiration demand was obtained using the Penman-Monteith model. ET de : In the formula: K c Represents the crop coefficient; This represents the rate of change of saturated vapor pressure with temperature. R n Indicates net radiation; RG Indicates soil heat flux; This represents the conversion factor between the wet-bulb and dry-bulb temperature difference and the water vapor pressure. v 2 indicates the wind speed at a height of two meters; T Indicates air temperature; e s Indicates saturated water vapor pressure; e a Indicates the actual water vapor pressure; The water stress index (SWD) is: In the formula: RD x Indicates the depth of the crop root zone; Accumulated vapor pressure deficit (VPD): In the formula: RJ(t) express t Rainfall during the period; VPD(t) express t Vapor pressure deficit at any given moment; The vapor pressure deficit during the crop growth period is calculated for each time period to obtain the cumulative vapor pressure deficit.

7. The method for multi-source data processing and dynamic crop yield estimation based on satellite remote sensing according to claim 1, characterized in that: Step S3 specifically involves: First, construct a field-level dynamic graph neural network structure that evolves with the reproductive period: In the formula: V t Represents a set of nodes, which has a number of nodes. v i,t Composition; each node v i,t Corresponding field i At any moment t eigenvectors s i,t The feature vector includes biophysical features, structural features, and environmental stress features; E t Represents an edge set and its weight. w i,j,t Dynamic quantification of fields i and j The spatiotemporal correlation between them, integrating soil texture and real-time meteorological similarities, is as follows: In the formula: Indicates field i , j Soil texture gradient; Q i,t , Q j,t Indicates field i , j The meteorological vector at time t; Indicates hyperparameters; Node features are aggregated into spatiotemporal information through graph convolution, and the propagation equation is: In the formula: Indicates the first l Hidden features in the first layer ; Indicates the first l The weight matrix of layer graph convolution; in, A t Indicates time t The original adjacency matrix of the lower and upper plots, i.e. ; Z t This represents the weighted value of historical features output by the lag effect module. To capture the cumulative effects of environmental factors, a learnable time-varying convolution kernel is designed. In the formula: S(tk) Indicates lag k Field characteristics during the period; W k Indicates the first k The convolutional kernels of the period are used to extract local correlations of historical features; The lag weight at time t; Lag effect output Z t As a supplementary input to the dynamic graph propagation equation, it enables the aggregation of node features to simultaneously include real-time information and historical cumulative influence, thus overcoming the limitation of traditional graph models that only rely on features at the current moment. Physical mechanism constraints: In the formula: The environmental stress and meteorological driving feature vectors at time t are represented; NPP represents net primary productivity. HI represents the distribution coefficient; HI represents the harvest index. f i,0 This represents the initial value of the allocation coefficient. f i,max This represents the maximum value of the allocation coefficient. The rate parameter represents the distribution phase; GDD represents growth degree-days, i.e., the number of degree-days in which the cumulative temperature exceeds the baseline temperature. s Indicates crop organs i The initial GDD that begins to dominate carbon allocation; GPP represents canopy primary productivity. This indicates the corrected method for maintaining respiration. R g The LAI represents basal respiration; LAI represents leaf area index, obtained through inversion by coupling the PROSPECT model and the SAIL model. V cmax Indicates the maximum carboxylation rate. Indicating intercellular spaces in leaf cells CO 2 concentrations, express CO 2 compensation points express CO Michaelis constant of 2 O 2 represents the oxygen concentration in the atmosphere. express O Michaelis constant of 2 J() PAR represents photosynthetically active radiation, indicating electron transport efficiency. g s () Indicates porosity, This indicates the concentration of carbon dioxide in the atmosphere. R d Indicates the crop's dark respiration rate; k m Represents the basic coefficient of crops. Q 10 This represents the temperature coefficient, where T represents the real-time temperature. T ref Indicates reference temperature. Indicates the moisture sensitivity coefficient; Indicates the assimilate conversion efficiency; Closed-loop training of physical constraints and dynamic graphs: Input: Dynamic graph node features: s i,t Integrating remote sensing inversion features, environmental stress features, and meteorological downscaling data; State initialization: Biomass initialization was obtained by inverting the RVOG model. SW i,0 , as the initial condition for the differential equation; Biomass at time t is predicted by aggregating spatiotemporal features using dynamic graphs. ; Physical deviation loss L phys : In the formula: This represents the differential equation obtained through physical constraints, namely: The final total loss function is: In the formula: Indicates the model prediction t Output at any given moment; This indicates the actual measured yield.

Images