Crop yield prediction method based on time series remote sensing images

Abstract

The present invention relates to the technical field of crop yield analysis, and specifically to a crop yield prediction method based on time-series remote sensing images. In the present invention, optical remote sensing images and synthetic aperture radar images during the crop growth period are acquired, and time-series images are constructed through complementary information fusion and cloud-underground information recovery. Cloud interference and data missing problems are addressed, and a multi-scale time-regulated convolution module is used to process the time-series images, extract multi-scale time-series features, and perform feature fusion in combination with attention differential jump connections to obtain a rice distribution map. Rice-planting areas are determined based on the rice distribution map, and an optimal assimilation period is selected within the rice-planting areas. A leaf area index obtained by remote sensing inversion is assimilated with crop growth model simulation results to generate an assimilated simulated yield. Based on the multi-scale time-series features, a parameter dynamic response matrix related to crop growth model parameters is established, and residual correction is performed on the assimilated simulated yield, thereby improving the accuracy of rice yield prediction.

Description

Technical Field

[0001] The present invention relates to the technical field of crop yield analysis, and in particular to a crop yield prediction method based on time-series remote sensing images. Background Art

[0002] Traditional crop yield prediction methods rely primarily on single remote sensing data sources (such as optical or radar imagery), which struggle to overcome cloud cover and missing spectral information, resulting in poor continuity in time series data. Existing technologies are susceptible to interference from clouds and rain, radar imagery has insufficient resolution, and multi-source data fusion methods fail to effectively extract multi-scale crop growth characteristics.

[0003] At the same time, crop model assimilation often uses fixed phenological parameters, ignoring the dynamic changes in the rice growth period and the temporal and spatial responses of remote sensing characteristics. This leads to biased assimilation period selection and accumulated yield residuals. Furthermore, existing models lack effective reconstruction of subcloud surface information and dynamic correction mechanisms for model parameters, resulting in errors in crop area identification and reduced yield prediction accuracy. Summary of the Invention

[0004] The purpose of the present invention is to provide a crop yield prediction method based on time series remote sensing images to solve the problems raised in the above-mentioned background technology. The core problems to be solved include how to construct high-precision time series images through complementary information fusion and under-cloud information recovery to solve the problems of cloud interference and data missing; how to achieve coordinated optimization of crop growth models and remote sensing inversion data based on multi-scale time series feature extraction and dynamic parameter correction, thereby improving the accuracy of rice yield prediction.

[0005] To achieve the above-mentioned object, the present invention provides the following technical solution: a method for predicting crop yield based on time-series remote sensing images, the method comprising the following steps:

[0006] S1. A dual-branch feature extraction module extracts the spectral features of optical images and the backscatter features of radar images. A cross-attention mechanism is used to calculate a complementary weight matrix. This supplements the optical spectral features in areas missing from radar images and reconstructs surface information in areas covered by optical clouds. A multi-scale feature fusion module is then used to integrate shallow details with deep semantic features to generate a fused image. This solution directly addresses the data loss issues caused by cloud occlusion in optical images and insufficient resolution in radar images, improving the integrity and spectral-texture consistency of time-series images.

[0007] By generating a cloud binary mask, we distinguish between areas completely obscured by thick clouds, areas covered by thin clouds, and areas of cloud shadows. In areas covered by thin clouds, spectral unmixing technology is used to separate the cloud layer and surface reflection signals. In areas completely obscured by thick clouds, the surface reflectivity beneath the clouds is reconstructed using a spatiotemporal series prediction model, combining cloud-free images from the previous phase with radar images from the same period. This scheme ensures the recovery of surface information in cloud-covered areas and guarantees the continuity of time series images.

[0008] S2: A multi-scale time-controlled convolution module is used to extract short-term phenological characteristics, medium-term trend characteristics, and long-term cycle characteristics using convolution kernels with different time windows. The gated recurrent unit is then used to model the periodicity of the rice growth period. This scheme enhances the ability to characterize crop growth dynamics.

[0009] The introduction of attention differential skip connections captures growth rate differences (such as the difference in growth rate between the tillering and heading periods) through element-by-element differential operations. This generates channel attention weights to enhance the feature response of key phenological periods (such as the heading period), fusing deep semantic features with shallow detail features across layers. This solution addresses the problems of traditional convolutional networks that ignore crop growth rate differences and blur field boundaries, improving the spatial accuracy of rice distribution maps.

[0010] S3. Remote sensing observation windows are divided based on rice phenological periods. Parameter sensitivity analysis is used to calculate the coefficient of determination and root mean square error (RMSE) between the assimilated yield and the measured yield for each period. The phenological period with the highest coefficient of determination, the lowest error, and dynamic synchronization with the critical yield period is selected as the optimal assimilation period. This solution addresses the period selection bias caused by traditional assimilation methods' reliance on fixed phenological periods.

[0011] During the heading period, the leaf area index is inverted based on the radiation transfer model, which drives the crop growth model to simulate the spatiotemporal distribution. The model state variables are updated through pixel-level assimilation iterations. This solution achieves dynamic coupling between remote sensing data and crop models, improving assimilation accuracy.

[0012] S4. Use the random forest algorithm to quantify the correlation strength between multi-scale time series features and crop model parameters, screen key time series features, and construct a dynamic response matrix with model parameters as row vectors and time series features as column vectors. This solution solves the problem of static parameters in traditional assimilation and realizes the dynamic correlation between parameters and remote sensing features.

[0013] By comparing the spatial residual distribution of assimilated and measured yields, and combining the dynamic response matrix to reversely deduce parameter deviations, the random forest feature contribution is used to identify high-contribution parameters. These are then input into the gated recurrent unit as correction factors to predict residual adjustments. The correction contributions of different growth stages are dynamically weighted through the attention mechanism. This scheme achieves spatially differentiated correction of residuals.

[0014] A correction coefficient for the maximum growth rate parameter of the leaf area index is generated at the pixel size. Multiple rounds of co-iterative optimization of parameters and residuals are performed using the alternating direction multiplier method until the coefficient of determination converges. This approach significantly reduces the systematic residuals of yield estimates through dynamic parameter correction and model rerunning.

[0015] Compared with the prior art, the present invention has the following beneficial effects:

[0016] Through the complementary fusion of optical and radar data and reconstruction of subcloud information, the integrity and availability of time series images are improved; multi-scale time series feature extraction and attention differential fusion mechanism enhance the ability to discriminate crop growth characteristics; dynamic parameter response matrix and residual correction optimization achieve deep collaboration between crop models and remote sensing inversion data, overcoming the assimilation period selection bias and parameter fixation problems in traditional methods; this method is used to improve the spatial accuracy of rice planting area identification and reduce yield prediction residuals, providing reliable technical support for agricultural precision management and food security decision-making. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 Schematic diagram of the method steps of the present invention. DETAILED DESCRIPTION

[0018] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making any creative efforts shall fall within the scope of protection of the present invention.

[0019] See also Figure 1 The present invention provides a technical solution: a method for predicting crop yield based on time series remote sensing images, comprising the following steps:

[0020] S1. Acquire optical remote sensing images and synthetic aperture radar images during the crop growth period, and construct time series images through complementary information fusion and cloud information recovery. Specifically, the following are involved:

[0021] Acquire time-series optical remote sensing images (such as Sentinel-2 multispectral data) and synthetic aperture radar images (such as Sentinel-1 data) throughout the crop growing season through satellite platforms; perform radiometric calibration and atmospheric correction on the optical images to eliminate the effects of atmospheric scattering; perform radiometric calibration, multi-look processing, and terrain correction on the radar images to eliminate geometric distortion caused by topographical fluctuations; and screen images covering the sowing to maturity period based on the crop phenological calendar, eliminating invalid images with cloud cover exceeding 50%.

[0022] A cross-modal fusion network based on the Swin-Transformer architecture is constructed to achieve feature-level complementary information fusion of optical and radar images. The complementary information fusion includes dual-branch feature extraction, complementary information enhancement module, and multi-scale feature fusion, among which:

[0023] Dual-branch feature extraction: Convolutional layers are used to extract the spectral features of optical images and the backscattering features of radar images, respectively. The spatial resolution differences between the two modalities are aligned through a spatial attention mechanism.

[0024] Complementary Information Enhancement Module: A cross-attention mechanism is designed in the Transformer encoder to calculate the complementary weight matrix of optical and radar features. For areas where radar features are missing (such as water bodies), optical spectral features are used to supplement them. For areas covered by optical clouds, radar penetration features are used to reconstruct surface information.

[0025] Multi-scale feature fusion: In the decoder stage, shallow detail features and deep semantic features are integrated to generate a fused image that preserves both spectral characteristics and texture details.

[0026] To address the cloud occlusion problem in optical images, a cloud-aware spatiotemporal information reconstruction mechanism is designed to achieve information recovery under the clouds. Specifically, the mechanism includes:

[0027] The shortwave infrared band of the optical image and the cloud detection algorithm are used to generate a cloud binary mask to distinguish between areas completely blocked by thick clouds, areas covered by thin clouds, and areas of cloud shadows.

[0028] For areas covered by thin clouds, spectral unmixing technology is used to separate the cloud layer and surface reflection signals. For areas completely obscured by thick clouds, the cloud-free image of the previous phase and the radar image of the same period are combined to reconstruct the surface reflectivity under the cloud through a spatiotemporal series prediction model (based on gated recurrent units). For cloud shadow areas, a spatiotemporal adaptive filtering algorithm is used to eliminate the reflectivity attenuation caused by shadows based on the backscattering coefficient of the radar image and the shortwave infrared band characteristics of the historical cloud-free image. Specifically, the cloud shadow boundary is detected by the ratio of the shortwave infrared band to the near-infrared band, and a shadow mask is generated. A linear regression model is established using the VV polarization backscattering coefficient of the radar image of the same period and the reflectivity of the corresponding phase of the historical cloud-free image to correct the reflectivity attenuation of the visible light band in the shadow area. The median filter of the time series sliding window is used to eliminate residual noise.

[0029] A cloud mask weighting term is introduced into the loss function, and the calculation formula is: total loss = image reconstruction loss + λ × cloud area reconstruction loss, where λ is dynamically adjusted according to the cloud thickness (for example, λ = 2.0 in thick cloud areas and λ = 0.5 in thin cloud areas) to ensure spectral continuity in the cloud edge area.

[0030] For time nodes with missing data, an adaptive weighted time series filtering algorithm was employed. Based on a cubic spline function, interpolation weights were dynamically adjusted based on the phenological characteristics of the crop growth curve (e.g., rapid growth during the tillering phase and steady changes during the heading phase) to ensure smoothness. Pixels with abnormal fluctuations were detected and corrected by calculating the difference in NDVI change rates (normalized difference in vegetation index change rates) between adjacent pixels over time. For areas with abnormal fluctuations over three consecutive time phases, a radar image-assisted verification mechanism was activated, using the temporal correlation of the backscatter coefficient to correct the data, ultimately constructing a time series image.

[0031] Step S1 synchronously acquires Sentinel-2 multispectral images and Sentinel-1 radar images covering the entire growth period of crops. First, atmospheric correction and geometric registration are performed on the optical image, and a cloud mask is generated using the shortwave infrared band to mark the cloud pollution area. At the same time, radiation calibration and terrain correction are performed on the radar image. Based on the dual-branch deep learning network, the spectral characteristics of the optical image and the scattering characteristics of the radar image are extracted respectively, and cross-modal associations are established through the cross-attention mechanism. In the vegetation coverage area, the spectral characteristics of the optical image are dominant and the backscattering characteristics of the radar image are supplementary. In the water area, the backscattering characteristics of the radar image are dominant and the backscattering characteristics of the optical image are supplementary. The short-wave infrared band features are used as a supplement to achieve the complementary advantages of cross-modal features; for cloud-blocked areas, thin cloud spectral unmixing and thick cloud spatiotemporal prediction (combining historical cloud-free images and radar data of the same period) are used to hierarchically repair the surface reflectivity, and a cloud mask weighted loss function is designed to optimize the reconstruction accuracy; finally, the missing phases are filled by adaptive cubic spline interpolation, combined with NDVI time series consistency test and radar backscatter verification, to generate spatially complete and temporally continuous time series images with a temporal resolution of 10 days. The reconstruction accuracy of the area under the cloud exceeds 90%, completely covering the key growth stages of crops from sowing to maturity, providing data support for precision agriculture monitoring.

[0032] S2. Use the multi-scale temporal control convolution module to process time series images, extract multi-scale temporal features, and combine them with attention differential jump connections for feature fusion to obtain a rice distribution map. Specifically,

[0033] Based on the time series images generated by S1, a three-level temporal feature extraction network is constructed using a multi-scale time-controlled convolutional module, including short-term phenological feature extraction, medium-term trend capture, and long-term cycle modeling.

[0034] Short-term phenological feature extraction: A 1×1 temporal convolution kernel is used to process adjacent frames of imagery (10-day window) to capture short-term variations such as the rapid growth of rice during the tillering period, and a 32-channel feature map is output.

[0035] Mid-term trend capture: 3×3 temporal convolution kernels are used to process three adjacent frames (30-day window) to extract progressive growth features from the jointing stage to the heading stage. Dilated convolution (dilation rate = 2) is used to expand the receptive field and output a 64-channel feature map.

[0036] Long-term cycle modeling: A 5×5 temporal convolution kernel is used to process five adjacent image frames (50-day window), combined with a gated recurrent unit (GRU) to memorize the periodicity of the entire rice growth period (from sowing to maturity), and output a 128-channel feature map.

[0037] Element-by-element differential operations are performed on feature maps of different scales (multi-scale time series features) to capture differences in crop growth rates. The sigmoid function is used to generate channel attention weights to dynamically enhance the feature responses of key phenological periods. The weighted deep semantic features and shallow detail features are fused across levels using jump connections. This enhances the spectral discrimination of highly sensitive periods such as the heading period while retaining the spatial details of the field boundaries. Finally, the rice distribution map is output through feature splicing and upsampling, achieving dual optimization of classification accuracy and phenological logic.

[0038] S3. Determine the rice planting area based on the rice distribution map, select the optimal assimilation period within the rice planting area, assimilate the leaf area index inverted by remote sensing with the crop growth model simulation results, and generate the assimilated simulated yield. Specifically, it includes:

[0039] Based on the rice distribution map output by S2, a rice planting binary mask with a spatial resolution of 10 meters was extracted, and non-planted areas (such as woodlands, bare land, and water bodies) were eliminated to generate a spatially continuous rice planting vector boundary. Combined with the administrative division data of the Chengdu Plain Economic Zone, the rice planting area was counted by county unit, which served as the basic spatial unit for regional yield estimation.

[0040] Based on the rice phenological periods (transplanting, tillering, jointing, heading, and milky stages), the time series images constructed by S1 were used to extract the remote sensing observation windows corresponding to each phenological period (e.g., heading corresponds to early to mid-August).

[0041] The leaf area index retrieved from remote sensing at different phenological periods was assimilated with the results of crop growth model simulation. The coefficient of determination between the assimilated yield and the measured yield at each period was calculated through parameter sensitivity analysis, and the model goodness of fit was evaluated using the root mean square error.

[0042] Based on the quantitative indicators of the coefficient of determination and the root mean square error, the phenological period with the highest coefficient of determination, the smallest model error, and the strongest synchronization with the dynamic changes of the leaf area index during the key period of yield formation is selected as the optimal assimilation period; for example, the heading period is determined as the optimal assimilation period because its leaf area index is highly correlated with the key period of yield formation;

[0043] Based on the radiation transfer model, the Sentinel-2 multispectral image and sun-sensor geometric parameters were input to invert the rice leaf area index at the heading stage using a lookup table method with a spatial resolution of 10 meters.

[0044] Drive a locally calibrated crop growth model, inputting meteorological data (sunshine, temperature and humidity), soil parameters (water retention and water conductivity), and management parameters (irrigation and fertilization) to simulate the spatiotemporal distribution of leaf area index at the heading stage;

[0045] The leaf area index retrieved from remote sensing is used as the observation value and assimilated with the crop growth model simulation results at the pixel level. The model state variables (leaf area index, biomass) are updated iteratively to generate the assimilated rice growth parameter set.

[0046] The assimilated leaf area index is input into the crop growth model to drive the model to complete the growth simulation from the heading period to the maturity period and output the yield per unit area. Based on the spatial weight of the rice distribution map, the total yield at the county scale is aggregated pixel by pixel and used as the simulated yield after assimilation.

[0047] S4. Based on the multi-scale time series features in step S2, a parameter dynamic response matrix related to the crop growth model parameters is established, residual correction is performed on the simulated yield after assimilation in step S3, and the optimized yield estimate of the rice planting area is output, specifically including:

[0048] Based on the multi-scale time series features extracted in step S2, a dynamic correlation between sensitive parameters of the crop growth model and remote sensing time series features is established by integrating short-term (tillering period), medium-term (jointing period), and long-term (heading period) time series features. Specifically, a random forest algorithm is used to quantify the nonlinear response intensity between time series features of different scales (such as the spectral variation coefficient of the tillering period and the texture entropy value of the jointing period) and crop growth model parameters (such as the maximum growth rate of the leaf area index and the accumulated temperature threshold). Through feature importance analysis and statistical significance testing, a set of key time series features sensitive to parameter changes is screened, and a multidimensional parameter-feature response matrix is ​​constructed; this matrix uses crop growth model parameters as row vectors and multi-scale time series features as column vectors, and the matrix element values ​​represent the dynamic response intensity of specific parameters to the corresponding features. During the matrix construction process, a gated recurrent unit network is introduced to encode the dynamic evolution pattern of the time series features, capturing the temporal dependency between parameters and features, thereby forming a parameter control knowledge base that supports dynamic optimization and provides a theoretical basis for parameter adjustment for subsequent residual correction.

[0049] Based on the spatial residual distribution of the simulated yield and the measured yield after assimilation in step S3, the key correction factors are screened in combination with the parameter dynamic response matrix. First, the yield residual value of each pixel is calculated by spatial statistical methods to generate a residual spatial distribution map. Then, the parameter-feature association relationship recorded in the parameter dynamic response matrix is ​​used to reversely deduce the potential model parameter deviation that causes the residual. Specifically, the corresponding multi-scale time series features of the high-value residual areas (such as spectral variation anomalies at the heading stage and sudden changes in texture entropy values ​​at the jointing stage) are extracted. Through random forest feature contribution analysis, the crop growth model parameters that significantly match the residual spatial distribution pattern (such as the parameter error of the maximum growth rate of leaf area index) are identified. The screened key correction factors (i.e., high-contribution parameters and their associated temporal features) are input into the gated recurrent unit network, which predicts spatially differentiated residual adjustments by learning the spatiotemporal coupling relationship between historical residual adjustments and multi-scale temporal features. In this process, the network dynamically weights the characteristic contributions of different growth periods through the attention mechanism to ensure the differentiated correction logic of the rapid growth characteristics of the tillering period and the stable period characteristics of the heading period, and finally outputs a residual adjustment field that matches the spatial heterogeneity of crop growth.

[0050] The predicted residual adjustment values ​​are applied to sensitive parameters of the crop growth model, enabling iterative optimization of model parameters and dynamic correction of yield estimates. For the maximum leaf area index growth rate parameter, a parameter correction coefficient is generated at the pixel level based on the spatial distribution of the residual adjustment field. The growth rate coefficient is reduced in areas with overestimated residuals to suppress excessive growth simulation, while it is increased in areas with underestimated residuals to enhance biomass accumulation. The corrected parameter set is then re-input into the crop growth model, driving the model's rerun from heading to maturity to generate growth simulation results with optimized parameters. This process employs the alternating direction multiplication method for parameter-residual collaborative optimization. Through multiple iterations, the true parameter values ​​are gradually approached until the coefficient of determination between the assimilated simulated yield and the measured yield reaches a convergence threshold. Finally, the optimized pixel-level yield estimates are spatially aggregated by county and, combined with administrative division weights, generate regional yield forecasts that meet agricultural management requirements. Furthermore, a leave-one-out cross-validation approach is used to evaluate the generalizability of the dynamic parameter response matrix, ensuring the stability of the optimization mechanism across different climate years and cropping systems.

[0051] The above shows and describes the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are merely preferred examples of the present invention and are not intended to limit the present invention. Various changes and improvements may be made to the present invention without departing from the spirit and scope of the present invention. Such changes and improvements fall within the scope of the present invention. The scope of protection claimed in the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for predicting crop yield based on time series remote sensing images, characterized in that: The method steps are as follows: S1. Acquire optical remote sensing images and synthetic aperture radar images during the crop growth period, and construct time series images through complementary information fusion and cloud information recovery; S2. Use the multi-scale temporal control convolution module to process time series images, extract multi-scale temporal features, and combine them with attention differential skip connections for feature fusion to obtain the rice distribution map; S3. Determine a rice planting area based on the rice distribution map, select an optimal assimilation period within the rice planting area, assimilate the leaf area index inverted by remote sensing with the crop growth model simulation results, and generate an assimilated simulated yield; S4. Based on the multi-scale time series characteristics in step S2, a parameter dynamic response matrix related to the crop growth model parameters is established, residual correction is performed on the simulated yield after assimilation in step S3, and the optimized rice planting area yield estimate is output; wherein: The construction process of the parameter dynamic response matrix specifically includes: A random forest algorithm was used to quantify the nonlinear correlation strength between each temporal feature and sensitive parameters of the crop growth model. A set of key temporal features sensitive to parameter changes was screened through feature importance analysis, and the dynamic evolution pattern of temporal features was encoded using a gated recurrent unit network. A parameter dynamic response matrix was constructed, with crop growth model parameters as row vectors and multi-scale temporal features as column vectors. The matrix element values ​​represent the dynamic response strength of the parameters to the features. The residual correction process specifically includes: By comparing the spatial residual distribution of simulated yield and measured yield after assimilation, the model parameter deviation leading to the residual is reversely deduced based on the parameter dynamic response matrix; the multi-scale time series features corresponding to the high-value residual areas are extracted, and the sensitive parameters matching the residuals are identified using random forest feature contribution analysis; the screened high-contribution parameters and their associated features are used as key correction factors and input into the gated recurrent unit network to predict the spatially differentiated residual adjustment amount, and the correction contribution of different growth period characteristics is dynamically weighted through the attention mechanism.

2. The method for predicting crop yield based on time series remote sensing images according to claim 1, characterized in that: The complementary information fusion process specifically includes: The spectral features of the optical image and the backscattering features of the radar image are extracted respectively through a dual-branch feature extraction module; The cross-attention mechanism is used to calculate the complementary weight matrix of optical and radar features, supplement optical spectral features in areas where radar features are missing, and reconstruct surface information in areas covered by optical clouds. A multi-scale feature fusion module is used to fuse shallow detail features and deep semantic features to generate a fused image that retains spectral characteristics and texture details.

3. The method for predicting crop yield based on time series remote sensing images according to claim 1, characterized in that: The process of recovering information off the cloud specifically includes: Generate cloud layer binary mask to distinguish the area completely blocked by thick clouds, the area covered by thin clouds and the area of ​​cloud shadow; In areas with thin cloud coverage, spectral unmixing technology is used to separate the cloud layer and surface reflection signals; For areas completely obscured by thick clouds, the cloud-free image of the previous phase and the radar image of the same period are combined to reconstruct the surface reflectivity under the clouds through a spatiotemporal series prediction model.

4. The method for predicting crop yield based on time series remote sensing images according to claim 1, characterized in that: The multi-scale time-controlled convolution module includes: Convolution kernels with different time windows are used to extract short-term phenological characteristics, medium-term trend characteristics and long-term periodic characteristics respectively; the periodic laws of the complete growth period of rice are modeled by gated recurrent units.

5. The method for predicting crop yield based on time series remote sensing images according to claim 1, characterized in that: The process of the attention differential jump connection specifically includes: Perform element-by-element difference operations on multi-scale time series features to capture differences in crop growth rates; Generate channel attention weights to enhance the characteristic responses of key phenological periods; Cross-level fusion of deep semantic features and shallow detail features preserves the spatial details of field boundaries and enhances spectral discrimination.

6. The method for predicting crop yield based on time series remote sensing images according to claim 1, characterized in that: The process of selecting the optimal assimilation period includes: Extract remote sensing observation windows corresponding to each growth stage based on rice phenological period; The leaf area index retrieved from remote sensing at different phenological periods was assimilated with the results of crop growth model simulation. The coefficient of determination between the assimilated yield and the measured yield at each period was calculated through parameter sensitivity analysis, and the model goodness of fit was evaluated using the root mean square error. According to the quantitative indicators of determination coefficient and root mean square error, the phenological period with the highest determination coefficient, the smallest model error and the strongest synchronization with the dynamic changes of leaf area index during the critical period of yield formation was selected as the optimal assimilation period.

7. The method for predicting crop yield based on time series remote sensing images according to claim 1, characterized in that: The assimilation process in step S3 includes: Inversion of rice leaf area index at heading stage based on radiation transfer model; Drive crop growth models to simulate the spatiotemporal distribution of leaf area index at heading stage; The remote sensing inversion results are assimilated with the crop growth model simulation results at the pixel level, and the crop growth model state variables are iteratively updated.

8. The method for predicting crop yield based on time series remote sensing images according to claim 1, characterized in that: The calculation of the optimized rice planting area yield estimate includes: Based on the spatial distribution of the residual adjustment amount, the correction coefficient of the maximum growth rate parameter of the leaf area index is generated according to the pixel granularity. The parameter value is lowered in the area where the residual is overestimated, and the parameter value is increased in the area where the residual is underestimated. The corrected parameter set is input into the crop growth model to drive the re-run, and the alternating direction multiplier method is used to perform multiple rounds of collaborative iterative optimization of parameters and residuals until the determination coefficient of the assimilated yield and the measured yield reaches the convergence threshold. Finally, the optimized yield estimate of the rice planting area is generated.

Images