Satellite remote sensing-based multi-source data processing and crop dynamic yield estimation method

Abstract

The application provides a satellite remote sensing-based multi-source data processing and crop dynamic yield estimation method, and relates to the field of crop yield estimation, and comprises the following steps: S1, multi-source data acquisition: collecting remote sensing data composed of optical images and radar images of crops in a region, meteorological data and corresponding auxiliary data, and preprocessing the data; 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 crop growth status and environmental stress from the preprocessed data; S3, a dynamic graph neural network coupled with a physical mechanism is constructed to realize dynamic yield estimation of crops. The method can realize dynamic monitoring and prediction of crop yield in the whole growth period through multi-source data of satellite images, and 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 for guiding field management and decision-making of agricultural producers, optimizing field management, reducing "supply-demand mismatch" of the industry chain, improving industry efficiency, promoting variety selection and iterating 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, i.e. the method of manually sampling and drawing statistical report forms in the field 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 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, which 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 achieved by the following technical solutions:

[0005] A multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing, comprising:

[0006] 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;

[0007] Step S2, key feature extraction and quantification: biological physical feature inversion, structure 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;

[0008] Step S3, constructing a dynamic graph neural network coupled with physical mechanisms: modeling spatio-temporal correlations through dynamic graph structures, embedding the growth physical mechanisms of crops, capturing the lag effects of environmental factors, and realizing dynamic yield estimation of crops.

[0009] Based on 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 digital elevation model (DEM), land use map, soil texture data, crop variety parameters and the like.

[0010] Based on further optimization of the above scheme, the data preprocessing in step S1 includes optical and radar data registration and meteorological data scale registration.

[0011] Optical and radar data registration: mutual information maximization method is used to realize data spatial alignment and eliminate spatio-temporal deviation:

[0012]

[0013] In the formula: l o represents the optical image; l r represents the radar image; represents the registration offset;

[0014] H represents the information entropy, which is used to measure the uncertainty of the gray scale distribution in the image:

[0015]

[0016] In the formula: l represents a single image; G represents the total number of gray values; p(g) represents the probability of occurrence of the gray value g; represents that when the radar image l r offset After that, the joint entropy of the optical-radar pixel pair is calculated.

[0017] Meteorological data scale registration: based on the physical constraint graph convolution, the low-resolution meteorological data is down-scaled to the field scale:

[0018]

[0019] In the formula: X i,t represents the field i at time tDownscaled weather data X(t 0 ) denotes low resolution weather data at initial time t 0. A denotes the adjacency matrix constructed based on DEM and land use map. I denotes the identity matrix. denotes the degree matrix of . W q denotes the weight matrix, b q denotes the bias vector.

[0020] Based on further optimization of the above scheme, the biophysical feature inversion in step S2 adopts the coupling inversion of leaf optical model PROSPECT and canopy radiation transfer model SAIL, specifically as follows:

[0021] Taking the reflectance spectrum of satellite remote sensing image , solar zenith angle , observation zenith angle , relative azimuth angle Δφ as input data, and taking leaf structure parameters N , soil background reflectance as prior parameters;

[0022] The PROSPECT model is used to simulate the reflectance and transmittance of single leaf:

[0023]

[0024] In the formula: C ab denotes the chlorophyll content; C w denotes the leaf water content; C m denotes the dry matter content;

[0025] The SAIL model is used to simulate the reflectance of canopy direction:

[0026]

[0027] In the formula: denotes the leaf area index;

[0028] Then, the parameters are inverted by minimizing the residual error between the observed spectrum and the simulated spectrum, and the objective function is:

[0029]

[0030] In the formula: denotes the observed spectrum; represents the simulated spectrum, which is obtained by coupling the PROSPECT model with the PROSPECT model:

[0031]

[0032] The sampling Levenberg-Marquardt algorithm is combined with the range constraint of the prior parameters to iteratively solve the minimum residual chlorophyll content C ab , leaf area index , dry matter content C m , etc.

[0033] and obtain the biomass SW :

[0034]

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

[0036] Based on the further optimization of the above scheme, the structure characteristic inversion in the step S2 is inverted by the RVoG model to obtain the height and structure parameters, specifically:

[0037] The collected radar echo is converted into the backscattering coefficient , and is preprocessed by denoising and geometric correction;

[0038] The RVoG model is constructed:

[0039]

[0040] In the formula: represents the surface backscattering coefficient; 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;

[0041] The crop height and the vertical structure parameters are obtained by inversion:

[0042]

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

[0044] Based on further optimization of the above scheme, the environmental stress characteristic quantification in step S2 includes a water stress index SWD and a cumulative vapor pressure deficit VPD:

[0045] Obtaining effective soil moisture SM av

[0046]

[0047] In the formula: SM cu represents the current soil moisture (volume water content); SM wi represents the soil moisture at the wilting point;

[0048] Obtaining crop evapotranspiration demand using the Penman-Monteith model ET de

[0049]

[0050] 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 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 soil); represents the conversion coefficient of the difference between dry and wet bulb temperatures 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;

[0051]

[0052] In the formula: RH represents the relative humidity; P represents the atmospheric pressure;

[0053] Then the water stress index SWD is:

[0054]

[0055] In the formula: RD x represents the crop root zone depth; ​​

[0056] Accumulated vapor pressure deficit (VPD):

[0057]

[0058] In the formula: RJ(t) express t Rainfall during a given period; VPD(t) express t Vapor pressure deficit at any given moment;

[0059] The vapor pressure deficit during the crop growth period is calculated for each time period to obtain the cumulative vapor pressure deficit.

[0060] Based on further optimization of the above scheme, step S3 specifically includes:

[0061] First, construct a field-level dynamic graph neural network structure that evolves with the reproductive period:

[0062]

[0063] In the formula: V t Represents a set of nodes, which has a certain 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 characteristics (chlorophyll content C). ab Leaf area index Biomass (SW, etc.), structural characteristics (crop height, etc.) The volume scattering coefficient of crops (etc.) and environmental stress characteristics (water stress index SWD, cumulative vapor pressure deficit VPD, etc.); 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:

[0064]

[0065] In the formula: Indicates field i , j Soil texture gradient (reflecting soil heterogeneity); Q i,t , Q j,t Indicates fieldi 、 j weather vector at time t; represents a hyper-parameter, used to regulate the influence strength of soil texture gradient difference on edge weight;

[0066] Node features are aggregated by graph convolution to achieve spatio-temporal information aggregation, and the propagation equation is:

[0067]

[0068] In the formula: represents the hidden features of the l layer, and the initial layer ; represents the weight matrix of the l layer graph convolution;

[0069]

[0070] wherein, A t represents the original adjacency matrix of the field block at time t ; ;

[0071] Z t represents the historical feature weighting value output by the lag effect module, and a learnable time-varying convolution kernel is designed to capture the lag effect for the cumulative effect of environmental factors:

[0072]

[0073] In the formula: S(t-k) represents the field block features in the lag k period; W k represents the convolution kernel of the k period, used to extract the local correlation of historical features; represents the lag weight at time t;

[0074]

[0075] In the formula: U represents a row vector (dimension consistent with the tanh output), used to map the nonlinear transformed features to scalar values; Y represents a transformation matrix, mapping the input to the input dimension of tanh; E t represents the weather stress factor at time t ; h t-1 represents the GRU hidden state at time t-1 ; represents the weather stress factor Et GRU hidden state h t-1 concatenated into a vector;

[0076] hysteresis effect output Z t As a supplementary input of the dynamic graph propagation equation, the node feature aggregation contains both real-time information and historical cumulative impact, solving the limitation of traditional graph models relying only on current features;

[0077] Physical mechanism constraints:

[0078]

[0079] In the formula: The environmental stress and meteorological driving feature vector at time t is represented; NPP represents net primary productivity; HI represents the harvest index; f i,0 The initial value of the distribution coefficient (i.e., the carbon distribution proportion of crop organs in the early growth period) is represented by i f i,max The maximum value of the distribution coefficient (i.e., the upper limit of the carbon distribution proportion of crop organs in the late growth period) is represented by i The rate parameter of the distribution stage is represented by; GDD represents the growing degree day, i.e., the degree day number of cumulative temperature exceeding the base temperature, GDD s The starting GDD of crop organs i starting to dominate carbon distribution is represented by; GPP represents canopy primary productivity, The corrected maintenance respiration is represented by R g The base maintenance respiration is represented by; LAI represents the leaf area index, which is obtained by coupling the PROSPECT model and the SAIL model, V cmax The maximum carboxylation rate (which is positively correlated with chlorophyll content C ab The concentration of intercellular space of leaf cells is represented by The 2 compensation point (increasing with temperature rise) is represented by CO The Michaelis constant of The 2 is represented by CO The oxygen concentration in the atmosphere is represented by The Michaelis constant of CO The 2 is represented by O The concentration of oxygen in the atmosphere is represented by The Michaelis constant of O The 2 is represented by J() The electron transport efficiency is represented by; PAR represents photosynthetically active radiation, g s ​​() Indicates porosity, This indicates the concentration of carbon dioxide in the atmosphere. R d This indicates the crop's dark respiration rate (i.e., respiration and digestion in the absence of light). 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;

[0080] Closed-loop training of physical constraints and dynamic graphs:

[0081] Input: Dynamic graph node features: s i,t Integrating remote sensing inversion features, environmental stress features, and meteorological downscaling data;

[0082] State initialization: Biomass initialization was obtained by inverting the RVOG model. SW i,0 , as the initial condition for the differential equation;

[0083] Biomass at time t is predicted by aggregating spatiotemporal features using dynamic graphs. ;

[0084] Physical deviation loss L phys :

[0085]

[0086] In the formula: This represents the differential equation obtained through physical constraints, namely:

[0087]

[0088] The final total loss function is:

[0089]

[0090] In the formula: Indicates the model prediction t Output per moment; This indicates the actual measured yield.

[0091] The following are the technical effects of the present invention:

[0092] The present application comprehensively depicts the growth state of crops in the field by inverting biophysical characteristics through the PROSPECT model and the SAIL model, inverting structural characteristics by using the RVoG model, and quantifying environmental stress characteristics, thereby avoiding the limitations of single data; then, the growth process of crops is embedded in a dynamic graph neural network, the dependence on redundant empirical data is effectively reduced by 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 by using a time-varying convolution kernel and a GRU module, effectively solves the defects of traditional models relying on "current time data", fully considers the hysteresis of yield estimation, and thus conforms 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 images, radar images, 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

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

[0094] The technical solutions in the embodiments of the present application will be clearly and completely described below. In the following description, specific details such as specific system structures, technologies, etc. are proposed for the purpose of explanation but not for the purpose of limitation, so as to thoroughly understand the embodiments of the present application.

[0095] Embodiment 1:

[0096] A multi-source data processing and crop dynamic yield estimation method based on satellite remote sensing, comprising:

[0097] 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, the optical images are acquired by Landsat or Sentinel-2, the radar images are acquired by Sentinel-1, the optical images and the radar images cover the whole growth period of crops; the meteorological data includes time series data of temperature, precipitation, humidity, etc. in the regional scale (which can be acquired by meteorological monitoring stations and corresponding sensors); the auxiliary data includes digital elevation model (DEM), land use map, soil texture data, crop variety parameters, etc.

[0098] and pre-processing the data, including optical and radar data registration and meteorological data scale registration;

[0099] Optical and radar data registration: A mutual information maximization method is used to achieve spatial alignment of data and eliminate spatiotemporal bias.

[0100]

[0101] In the formula: l o Represents optical images; l r Represents radar imagery; Indicates the registration offset;

[0102] H Information entropy is used to measure the uncertainty of gray-level distribution in an image.

[0103]

[0104] 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;

[0105] Meteorological data scaling: Physically constrained graph convolution to downscale low-resolution meteorological data to the field scale.

[0106]

[0107] 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.

[0108] Step S2, key feature extraction and quantification: biological and physical features inversion, structure features inversion and environmental stress features quantification are used to extract biological and physical features reflecting crop growth status and environmental stress from pre-processed data;

[0109] Biophysical features inversion uses the coupling inversion of leaf optical model PROSPECT and canopy radiation transfer model SAIL, specifically as follows:

[0110] The reflectance spectrum of satellite remote sensing image (multi-spectral / high-spectral, covering 400-2500nm), solar zenith angle , observation zenith angle , relative azimuth angle Δφ (computed from satellite orbit and observation time) as input data, and 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;

[0111] The reflectivity and transmissivity of single leaf are simulated by PROSPECT model:

[0112]

[0113] In the formula: C ab Chlorophyll content is represented by Chl; C w Leaf water content is represented by Cw; C m Dry matter content is represented by Cd;

[0114] The reflectivity of canopy in different directions is simulated by SAIL model:

[0115]

[0116] In the formula: Leaf area index is represented by L;

[0117] Then, the parameters are inverted by minimizing the residual error between observed spectrum and simulated spectrum, and the objective function is as follows:

[0118]

[0119] In the formula: Observed spectrum is represented by Y; Simulated spectrum is represented by Y, which is obtained by coupling PROSPECT model and PROSPECT model:

[0120]

[0121] The Levenberg-Marquardt algorithm is used, combined with pre-defined range constraints on prior parameters, to iteratively solve for the chlorophyll content that minimizes the residual. C ab Leaf area index Dry matter content C m wait;

[0122] And obtain biomass SW :

[0123]

[0124] 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.

[0125] Structural feature inversion involves retrieving height and structural parameters using the RVOG model, specifically as follows:

[0126] The collected radar echoes are converted into backscattering coefficients. The data underwent preprocessing, including denoising (using Lee filtering to suppress speckle noise) and geometric correction (converting to a ground coordinate system and extracting pixels from the crop planting area).

[0127] Building the RVOG model:

[0128]

[0129] In the formula: This represents the surface backscattering coefficient (obtained through bare ground pixel inversion). The volume scattering coefficient of a crop (reflects the volume scattering ability of the crop's vertical structure and is related to stem density and distribution). Indicates the radar wavelength; Indicates the radar incident angle; Indicates crop height;

[0130] The crop height and vertical structure parameters were obtained through inversion:

[0131]

[0132] Crop height is solved by nonlinear optimization. Volume scattering coefficient of crops .

[0133] The quantification of environmental stress characteristics includes the water stress index (SWD) and the cumulative vapor pressure deficit (VPD):

[0134] Obtain effective soil moisture SMav :

[0135]

[0136] wherein: SM cu represents the current soil moisture (volume water content); SM wi represents the soil moisture at wilting point (obtained from empirical values and soil texture standards, e.g. clay 0.1, sandy soil 0.03);

[0137] Penman-Monteith model is used to obtain the crop evapotranspiration requirement ET de :

[0138]

[0139] wherein: K c represents the crop coefficient (adjusted according to the crop growth stage); represents the rate of change of saturated vapor pressure with temperature (reflecting the effect of temperature on saturated vapor pressure); R n represents the net radiation (representing the energy driving evapotranspiration); RG represents the soil heat flux (i.e. the amount of heat exchanged between the surface and the deeper layers of the soil); represents the conversion coefficient of the difference between dry and wet bulb temperature and vapor pressure (related to atmospheric pressure); v 2 represents the wind speed at two meters height; T represents the air temperature; e s represents the saturated vapor pressure (i.e. the maximum vapor pressure that the air can contain at the current temperature); e a represents the actual vapor pressure (obtained from the relative humidity);

[0140]

[0141] wherein: RH represents the relative humidity; P represents the atmospheric pressure;

[0142] The water stress index SWD is then:

[0143]

[0144] wherein: RD x represents the crop root zone depth;

[0145] Cumulative vapor pressure deficit VPD:

[0146]

[0147] In the formula: RJ(t) represents t the rainfall of the period; VPD(t) represents t the vapor pressure deficit at the moment;

[0148] The vapor pressure deficit in the growth period of the crop is calculated by time period to obtain the cumulative vapor pressure deficit.

[0149] 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 mechanism of the crop, capturing the lag effect of environmental factors, and realizing the dynamic yield estimation of the crop; specifically:

[0150] First, a field-level dynamic graph neural network structure evolving with the growth period is constructed:

[0151]

[0152] In the formula: V t represents a node set, which has several nodes v i,t ; each node v i,t corresponds to a field i ; at the moment t , the feature vector s i,t includes biophysical features (chlorophyll content C ab , leaf area index , biomass SW, etc.), structural features (crop height , crop volume scattering coefficient , etc.), and environmental stress features (water stress index SWD, cumulative vapor pressure deficit VPD, etc.); E t represents 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:

[0153]

[0154] In the formula: represents the soil texture gradient of the fields i , j (representing soil heterogeneity); Q i,t , Qj,t representing a field block i , j weather vector at time t; representing a hyperparameter for regulating the influence strength of soil texture gradient difference on edge weight;

[0155] Node features aggregate spatio-temporal information through graph convolution, and the propagation equation is:

[0156]

[0157] In the formula: representing the hidden features of the first l layer, the initial layer ; representing the weight matrix of the first l layer of graph convolution;

[0158]

[0159] wherein, A t representing the original adjacency matrix of the field block at time t , i.e. ;

[0160] Z t representing the historical feature weighting value output by the lag effect module, designed to capture the cumulative effect of environmental factors by learning time-varying convolution kernels:

[0161]

[0162] In the formula: S(t-k) representing the field features at the lag k period (for example: t-k SWD at time t); W k representing the convolution kernel of the first k period, used to extract local correlations of historical features; representing the lag weight at time t;

[0163]

[0164] In the formula: U represents a row vector (dimension consistent with the tanh output), used to map the nonlinearly transformed features to scalar values; Y represents a transformation matrix, mapping the input to the input dimension of tanh; E t representing the weather stress factor at time t ; h t-1 representing the weather stress factor at time t-1The 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;

[0165] 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.

[0166] Physical mechanism constraints:

[0167]

[0168] In the formula: The environmental stress and meteorological driving feature vectors at time t are represented; NPP represents net primary productivity. Indicates the allocation coefficient; HI represents the harvest index (dynamically adjusted based on the number of varieties adopted and the growth period). f i,0 This represents the initial value of the allocation coefficient (i.e., early growth period, crop organ). i (carbon allocation percentage) f i,max This represents the maximum value of the allocation coefficient (i.e., late growth period, crop organ). i (the upper limit of carbon allocation ratio). Rate parameters representing the distribution stage (set according to crop type and actual conditions, such as wheat grain distribution). GDD represents growth degree-days, which is 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 (which is related to chlorophyll content). C ab Positive correlation), that is V cmax =k C ab +b (k and b represent crop-specific parameters), Indicating intercellular spaces in leaf cells CO2 Concentration (affected by stomatal conductance and photosynthetic consumption, which can be obtained by hyperspectral red edge range characteristics of remote sensing inversion), denotes CO 2 Compensation point (increases with temperature rise), denotes CO 2 Michaelis constant, O 2 Denotes the oxygen concentration in the atmosphere, denotes O 2 Michaelis constant, J() denotes electron transport efficiency, and PAR represents photosynthetically active radiation, g s () denotes stomatal conductance, denotes the concentration of carbon dioxide in the atmosphere, R d denotes the crop dark respiration rate (i.e. respiration digestion under no light); k m denotes the crop base coefficient, Q 10 denotes the temperature coefficient, T denotes the real-time temperature, T ref denotes the reference temperature, denotes the water sensitive coefficient; denotes assimilate conversion efficiency (determined according to different varieties of crops);

[0169] LAI:

[0170] Using the PROSPECT model, input the 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:

[0171]

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

[0173]

[0174] According to the real crown layer reflectivity obtained from the remote sensing image, calculate the observation error between it and the simulated crown layer reflectivity :

[0175]

[0176] By iteratively adjusting LAI, minimize Error until the simulated value and the observed value optimal match, the LAI at this time is the inversion result;

[0177]

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

[0179]

[0180] In the formula: A represents net photosynthesis;

[0181] Physical constraints and dynamic graph closed-loop training:

[0182] Input: dynamic graph node features: s i,t Fusion remote sensing inversion characteristics, environmental stress characteristics and meteorological downscaling data;

[0183] State initialization: through the RVoG model inversion biomass initialization, get SW i,0 , as the initial condition of the differential equation;

[0184] Through the dynamic graph aggregation of spatio-temporal characteristics to predict the biomass at time t ;

[0185] Physical deviation loss L phys :

[0186]

[0187] In the formula: The differential equation obtained by physical mechanism constraint, that is:

[0188]

[0189] The final total loss function is:

[0190]

[0191] In the formula: Y (t) represents the yield predicted by the model at time t (obtained by predicted biomass t And the harvest index); Y (t) represents the measured yield (obtained by field sampling or statistics).

[0192] Example 2: ​

[0193] 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:

[0194] Taking t=0 as the starting point of the growth period, input the initialization features, including: soil texture at the sowing stage, initial biomass, initial meteorological values, etc.; and construct the initial graph: initialize the edge weight based on soil texture similarity w i,j,0 ;

[0195] Time series iteration:

[0196] 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, and update the hidden layer through the propagation equation ; then, calculate the biomass:

[0197]

[0198] In the formula: represents the output layer mapping function;

[0199] and correct the predicted value through the differential equation constraint;

[0200] Finally, based on the biomass and the crop harvest index, the yield prediction at time t is obtained:

[0201] .

Claims

1. A satellite remote sensing-based multi-source data processing and crop dynamic yield estimation method, characterized in that: Comprise: Step S1, multi-source data acquisition: collect remote sensing data composed of optical images and radar images of crops in the region, meteorological data and corresponding auxiliary data, and pre-process 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 pre-processed data; wherein, biophysical feature inversion uses the coupling inversion of leaf optical model PROSPECT and canopy radiation transfer model SAIL, structural feature inversion uses RVoG model to invert height and structural parameters, and environmental stress feature quantification includes water stress index SWD and cumulative vapor pressure deficit VPD; Step S3, constructing a dynamic graph neural network coupled with physical mechanism: through dynamic graph structure modeling, the growth physical mechanism of crops is embedded, the lagging effect of environmental factors is captured, and the dynamic yield estimation of crops is realized, specifically: Firstly, a field-level dynamic graph neural network structure evolving with growth period is constructed: In the formula: V t denotes a set of nodes, which is composed of several nodes v i,t ; each node v i,t corresponds to a field block i At time t , the feature vector s i,t , the feature vector includes biophysical characteristics, structural characteristics and environmental stress characteristics; E t denotes a set of edges, and the edge weight w i,j,t quantifies the spatio-temporal correlation between the field blocks i and j , fuses soil texture and real-time meteorological similarity, specifically: wherein: represents a field block i , j of a soil texture gradient; Q i,t , Q j,t represents a field block i , j a weather vector at time t; represents a hyperparameter; Node features are aggregated through graph convolution to realize temporal and spatial information aggregation, and the propagation equation is: In the formula: represents the first l layer hidden feature, the initial layer ; represents the first l weight matrix of layer graph convolution; wherein A t representing the time instant t Next, the original adjacency matrix of the field plot, i.e. ; Z t The historical feature weighting value representing the hysteresis effect module output is designed to capture the cumulative effect of the environmental factors by learning a time-varying convolution kernel: In the formula: S(t-k) represents a lag k period field characteristics of the field; W k represents the first k period convolution kernel for extracting local correlation of historical features; represents the lag weight at time t; Physical mechanism constraints: 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, f i,max represents the maximum value of the distribution coefficient, represents the rate parameter of the distribution stage; GDD represents the growing degree day, that is, the degree day number of cumulative temperature exceeding the base temperature, GDD s represents the crop organ i starting GDD at which the carbon distribution begins to be dominated; GPP represents canopy primary productivity, represents the corrected maintenance respiration, R g represents the base maintenance respiration; LAI represents the leaf area index, which is obtained by coupling the PROSPECT model and the SAIL model, V cmax represents the maximum carboxylation rate, C ab represents the chlorophyll content, represents the CO 2 concentration of leaf intercellular space, represents CO 2 compensation point, represents CO 2 Michaelis constant, O 2 represents the oxygen concentration in the atmosphere, represents O 2 Michaelis constant, J() represents the electron transport efficiency, and PAR represents the photosynthetic active radiation. g s () represents the stomatal conductance, represents the carbon dioxide concentration in the atmosphere, R d represents the crop dark respiration rate; 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 graph: Input: dynamic graph node features: s i,t Fusing remote sensing inversion features, environmental stress features, and meteorological downscaling data; State initialization: initialization of biomass by inversion of the RVoG model, yielding SW i,0 as initial conditions for the differential equations; Predicting biomass at time t by aggregating spatio-temporal features from dynamic graphs ; physical bias loss L phys : wherein: represents a differential equation obtained by constraining through a physical mechanism, i.e.: The final total loss function is: In the formulae: represents the model predicted output at time t; t represents the measured output.​ 2. The method according to claim 1, characterized in that: The optical image is obtained by Landsat or Sentinel-2, the radar image is obtained by Sentinel-1, and the optical image and the radar image cover the whole growth period of the crop; the meteorological data includes temperature, precipitation, and humidity time series data in the regional scale; the auxiliary data includes digital elevation model, land use map, soil texture data, and crop variety parameters.

3. The method according to claim 1 or 2, characterized in that: The data pre-processing in step S1 includes optical and radar data registration and meteorological data scale registration; Optical and radar data registration: mutual information maximization method is used to realize data spatial alignment and eliminate temporal and spatial deviation: In the formulae: l o represents an optical image; l r represents a radar image; represents a registration offset; H denotes the information entropy, which measures the uncertainty of the gray level distribution in the image: wherein: l denotes a single image; G denotes the total number of grey values; p(g) denotes the probability of occurrence of a grey value g; denotes when the radar image l r offset after which the joint entropy of the optical-radar pixel pairs is calculated; Meteorological data scale registration: based on physical constraint graph convolution, low-resolution meteorological data is scaled down to field scale: wherein: X i,t representing a field block i at time t downscaled weather data; X(t 0 ) representing an initial time t 0; A representing a matrix of DEM and land use map; I denotes the identity matrix; denotes the degree matrix of W q denotes a weight matrix, b q denotes a bias vector.

4. The method according to claim 3, characterized in that: In step S2, biophysical feature inversion uses the coupling inversion of leaf optical model PROSPECT and canopy radiation transfer model SAIL, specifically: with satellite remote sensing images of the reflected spectrum , solar zenith angle , observation zenith angle , relative azimuth angle Δφ as input data, with leaf structure parameters N , soil background reflectance as prior parameters; PROSPECT model is used to simulate the reflectivity and transmissivity of single leaf: wherein: C ab represents the chlorophyll content; C w represents the leaf water content; C m represents the dry matter content; SAIL model is used to simulate the reflectivity of crown direction: In the formulae: represents the 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: In the formulae: represents the observed spectrum; represents the simulated spectrum, obtained from the PROSPECT model coupled with the PROSPECT model: The sampling Levenberg-Marquardt algorithm, combined with the range constraint of prior parameters, iteratively solves the minimum residual chlorophyll content C ab , leaf area index , dry matter content C m ; and obtain biomass SW : where: LMA represents the leaf weight, which is derived from the dry matter content obtained by inversion of the PROSPECT model C m is derived from the leaf thickness.

5. The method according to claim 3, characterized in that: In step S2, structural feature inversion uses RVoG model to invert height and structural parameters, specifically: Converting the collected radar echoes into backscatter coefficients and pre-processed by de-noising, geometric correction; RVoG model is constructed: wherein: represents the surface backscatter coefficient; represents the volume scattering coefficient of the crop; represents the radar wavelength; represents the radar incidence angle; represents the crop height; Crop height and vertical structure parameters are obtained by inversion: Solving crop height by non-linear optimization with the volume scattering coefficient of the crop .

6. The method according to claim 1, characterized in that: In step S2, environmental stress feature quantification includes water stress index SWD and cumulative vapor pressure deficit VPD: Acquiring available soil moisture SM av : wherein: SM cu represents the current soil moisture; SM wi represents the soil moisture at the wilting point; Penman-Monteith model is used to obtain crop evapotranspiration requirement ET de : In the formula: K c represents a crop coefficient; represents a rate of change of saturated water vapor pressure with temperature; R n represents net radiation; RG represents soil heat flux; represents a conversion coefficient of the difference between dry-bulb temperature and wet-bulb temperature and water vapor pressure; v 2 represents a wind speed at a height of two meters; T represents air temperature; e s represents saturated water vapor pressure; e a represents actual water vapor pressure; Then the water stress index SWD is: In the formulae: RD x denotes the depth of the crop root zone; Cumulative vapor pressure deficit VPD: In the formulae: RJ(t) represents t the rainfall amount of the period; VPD(t) represents t the vapor pressure deficit at the time The vapor pressure deficit during the crop growth period is calculated by time period to obtain the cumulative vapor pressure deficit.

Images