Enhanced Feature Extraction From Assimilated VTCI and LAI With a Particle Filter for Wheat Yield Estimation Using Cross-Wavelet Transform

To further reveal the relationships between different variables and yield at each growth stage of winter wheat, an approach for estimating regional yields of winter wheat at multiple time scales was developed by assimilating the CERES-Wheat model simulations and remotely sensed observations. Specifically, the particle filter assimilation algorithm was chosen to assimilate the simulated soil moisture at the depth of 0–20 cm and leaf area index (LAI) and MODIS retrieved vegetation temperature condition index (VTCI) and LAI. The resonance periods of time series assimilated VTCIs and LAIs at different growth stages of winter wheat with crop yield were analyzed separately using the cross-wavelet transform to determine the variation relationships between the assimilated variables and yield at multiple time scales, and the calculated weights of assimilated VTCI and LAI at each growth stage of winter wheat were used to establish a yield estimation model. Both assimilated VTCI and LAI could comprehensively integrate the effects of the CERES-Wheat model simulations and remotely sensing observations, and cross-wavelet transformed time series VTCIs and LAIs at each growth stage had specific resonance periods with the time series yields, regardless of whether they were assimilated or not. Compared with the unassimilated variables, assimilated VTCI and LAI were given greater weights at the key growth stages, namely VTCI at the jointing and heading-filling stages and LAI at the heading-filling and milk maturity stages, enhancing feature extraction and the accuracy of yield estimation was improved.


I. INTRODUCTION
A S ONE of the major food crops, timely and accurate estimation of regional wheat yields can provide effective support for agricultural management and food policy formulation [1]. In recent years, remote sensing technology has been widely used for regional yield estimation because of its advantages of being able to acquire data over large areas and rapidly. At the same time, it also has the shortcoming of lacking mechanical descriptions of crop growth, crop development, and yield formation as well as the relationships between crop growth and meteorological, soil, and environmental factors. Significantly, crop growth models can continuously simulate the dynamic response of crop growth and development by inputting data, including weather, soil, crop cultivars, and agronomic management practices [2]. However, limitations in data availability and quality have restricted the ability of crop growth models to precisely quantify the spatial heterogeneity of environment and crop growth, and consequently affected the accuracy of yield estimation at the regional scale, particularly for the simulation of crop canopy development and soil moisture [3], [4], [5]. Based on previous work, the real-time, macroscopic characteristics of remotely sensed data and the continuity and mechanics of crop growth models are complementary. Data assimilation, as an effective method to combine remote sensing information and crop growth models [6], has been recently used in regional yield estimation studies with high accuracy [7], [8], [9], [10].
Data assimilation methods can mainly be divided into variational assimilation approaches and sequential assimilation approaches. In contrast to variational assimilation approaches, which solve for the global optimum, sequential assimilation approaches emphasize solving for the optimal estimate at a single moment, i.e., the forecast field of the model is continuously updated with new observations so that the initial field at the next moment is formed and updated sequentially to obtain the overall optimal solution. The ensemble Kalman filter (EnKF) algorithm and the particle filter (PF) algorithm are representative methods in sequential assimilation approaches. In particular, PF utilizes the entire probability density function of the model states given the observations in computing the posterior weights for resampling prior states and parameters [11], [12], [13], [14]. Compared with the EnKF, PF is not strictly limited by Gaussian distributions [15]. Currently, the PF algorithm has been widely This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. used in crop growth monitoring and yield estimation studies [16], [17].
In recent years, when using the assimilation method for regional yield estimation, variables that are closely related to crop growth are often selected. Among them, the leaf area index (LAI) is an important variable to characterize light interception and absorption, vegetation growth conditions, and productivity [18]. Many studies have demonstrated that the assimilated LAI could reduce model uncertainty and improve the accuracy of yield estimates [19]. In addition to LAI, soil moisture is an important state variable and a major controlling factor for crop water stress. In the assimilation system, observations of the volumetric moisture content in the top few centimeters of the soil layer could be taken from remote sensing sources (e.g., ESA CCI, SMOS, and SMAP) [20]. However, the accuracy of these remotely sensed soil moisture products decreases as the moisture content of vegetation increases as the measurement is unable to sufficiently penetrate through dense vegetation to detect surface soil moisture [20], [21], [22]. Based on the triangular distribution of land surface temperature (LST) and normalized difference vegetation index (NDVI) at the region level, Wang et al. [23] proposed the vegetation temperature condition index (VTCI) to monitor the crop water stress. Sun et al. [24] explored the relationship between the soil moisture and VTCI in the Guanzhong Plain, China, and obtained the conclusion that VTCI is a near real-time drought monitoring index, which has a significant linear correlation with soil moisture at the depth of 0-20 cm. Xie et al. [25] obtained remotely observed soil moisture in the assimilated systems based on the linear relationship between VTCI and soil moisture at 0-20 cm depth to achieve regional yield estimation, further demonstrating the feasibility of using the relationship between VTCI and soil moisture for research. Han et al. [26] further integrated the time series LAIs and VTCIs into the crop growth model to achieve high accuracy wheat yield estimation, demonstrating the advantage of using time series data of the two variables for yield estimation.
When using the assimilated variables for regional yield estimation, fewer studies have explored the relationships between the time series of the assimilation variables at each growth stage and the time series yields, while multiple time scales help to analyze the effects of characteristic variables on final yield at different time scales, thus revealing underlying rules hidden in the time series. As an effective time series analysis method, wavelet transform has the feature of time-frequency multiresolution, which allows a more detailed description of the internal structure of time series. Wavelet transform has been widely used in the field of remote sensing such as vegetation phenology analysis [27], image fusion and compression [28], and classification [29]. With the continuous development of the wavelet transform, the cross-wavelet transform has emerged, a technique for multisignal, multiscale analysis based on the wavelet transform [30]. The cross-wavelet transform combines wavelet transform and cross-spectrum analysis to characterize the correlation degree and phase relationship of two time series at different time scales, and the analysis method based on cross-wavelet transform has been widely used in drought [31], and hydrology [32]. Characteristic variables at different growth stages respond to different degrees of crop yield, and the cross-wavelet transform can analyze the relationship between time series at multiple time scales. Importantly, to the best of our knowledge, there are few reports of regional yield estimation using a combination of data assimilation and cross-wavelet transform.
This article proposed a regional yield estimation approach that combined data assimilation and cross-wavelet transform. In this study, the PF assimilation method was used to assimilate soil moisture at the depth of 0-20 cm and LAI from the CERES-Wheat model simulations with remotely sensed retrieved VTCI and LAI, allowing the assimilation results to fully integrate the effects of the crop growth model simulations and remote sensing observations. The cross-wavelet transform was used to analyze the correlation between time series assimilated variables at each growth stage of winter wheat and time series yields to achieve a high accuracy regional yield estimation. In summary, the objectives are to 1) achieve two-parameter assimilation of VTCI and LAI, and to further explore the effect of data assimilation on yield estimation results; 2) determine the period variation relationships that exist between time series assimilated variables and time series yields at the temporal scale; 3) explore the applicability of combining the data assimilation method with cross-wavelet transform for regional yield estimation.

A. Study Area
The study area is the Guanzhong Plain, which is located in the central part of Shaanxi province, China (see Fig. 1). The climate is a temperate continental monsoon climate, with an annual average temperature ranging from 6°to 13°C. It is warm and rainy in summer, and cold and dry in winter, with average annual precipitation ranging from 500 to 700 mm. The soil types of Guanzhong Plain include manural loessial soil, cinnamon soil, dark loessial soil, and yellow-cinnamon soil. Among them, manural loessial soil is the main cultivation soil. It is formed in natural cinnamon soil through the long-term application of manure and maturation of tillage and is a powdered loam with high organic matter and nitrogen and phosphorus nutrient contents in the surface layer, which is fertile and has a long cultivation period and good arable properties. Due to its suitable geographical conditions, it has become an important commercial grain production base. The main type of cultivation is the rotation of winter wheat and summer maize in irrigated areas and there is only one wheat cultivation in a year in rain-fed areas. Winter wheat is sown in early or mid-October and harvested in early June of the following year. Considering the growth of winter wheat, this article focused on the overwintering stages, i.e., March to May. It can be divided into four growth stages [33]: the green-up stage is from early March to mid-March, the jointing stage is from late March to mid-April, the heading-filling stage is from late April to early May, and the milk maturity stage is from mid-May to late May. B. Remote Sensing Data 1) MODIS-retrieved VTCI: Daily LST and NDVI from MODIS MYD11A1 and MYD09GA products, with a spatial resolution of 1 km from March to May in the years 2011 to 2020 were calculated and retrieved. Next, using the maximum value composite and minimum value composite methods, the VTCI for the 10-day interval was calculated using a 10-day time scale according to the following formulas [24], [34]: where LST max (NDVI i ) and LST min (NDVI i ) are the maximum and minimum LST values of all pixels in the study area when the NDVI value is equal to NDVI i , which are called the warm and cold edges, respectively. LST(NDVI i ) is the LST value when the NDVI value is equal to NDVI i . The coefficients a, b, a , and b are approximated by the scatter plots of NDVI and LST in the study area. The average of the VTCI for the 10-day intervals included at each growth stage was calculated as the VTCI value for that growth stage. Based on the distribution map of winter wheat planting areas in each county of the Guanzhong Plain, the average VTCI value of the pixels contained in the growing areas in each county was taken as the VTCI value for that growth stage in that region for that year.
2) MODIS-retrieved LAI: LAI was extracted using the MODIS MCD15A3H product with a spatial resolution of 500 m and a temporal resolution of 4 days. Due to the influence of clouds and shadows, the MCD15A3H product has data gaps. To address this issue, the original time series LAIs were smoothed by using the upper envelope Savitzky-Golay (S-G) filter [35]. The LAI smoothed by the upper envelope S-G filter is more consistent with the actual growth of winter wheat. To ensure that LAI and VTCI had the same value range and spatial resolution, the LAI after the S-G filter was normalized to 0-1, and the VTCI was resampled to 500 m by using the nearest neighbor method. The maximum value of LAI contained every 10 days was taken as the LAI at 10-day intervals, and the LAI of each growth stage was calculated from the LAI of multiple 10-day intervals contained in that stage. Similar to VTCI, based on the distribution map of winter wheat planting areas in the counties of the Guanzhong Plain, the average LAI of the pixels in the planting areas in each county was taken as the LAI value of the region for the winter wheat growth stages of the year.

C. Field Measured Data
In this article, 12 typical planting sites in the Guanzhong Plain, all located in the main winter wheat growing areas from 2011 to 2020, were selected as sampling sites, and their distribution was shown in Fig. 1. The measurements obtained from each sampling site were used to calibrate and run the CERES-Wheat model. Specifically, during the sowing period, soil physicochemical parameters such as pH value and soil moisture were measured in each of the seven layers of the soil profile (0-12 cm, 12-20 cm, 20-50 cm, 50-80 cm, 80-120 cm, 120-160 cm, 160-200 cm) at each sampling site. At the jointing and heading stages, record the cultivar, sowing method, height, planting density, coverage rate, fertilization date, fertilization amount, irrigation date, irrigation amount, and other information of winter wheat at each sampling site, and LAI and soil moisture at the depth of 0-20 cm were measured.

D. CERES-Wheat Simulated Data
The CERES-Wheat model, as an individual crop growth module of Decision Support System for Agrotechnology Transfer, can specifically simulate the development, growth and senescence, biomass accumulation and partitioning, and leaf and grain growth of wheat [36]. The input parameters of the CERES-Wheat model mainly included crop management information, soil characteristics, weather data, and crop genetic information [37]. Among them, field management information was obtained through field investigation and actual measurements, including crop cultivar, sowing date, planting density and depth, fertilizer application date and quantity, irrigation date and quantity, and previous crop type. Soil data include those measured in Section II-C above. Weather data were measured at 43 weather stations distributed in the Guanzhong Plain, and the weather data of each sampling site were measured by the nearest weather station, including daily maximum and minimum air temperatures, precipitation, and solar radiation. Genetic parameters control growth and development processes and are closely related to plant morphological development and grain yield [25]. Before applying the model to the Plain, it is necessary to calibrate crop genetic parameters based on the "trial and error" method to obtain more accurate simulation results [37]. The main aspects were 1) to verify whether the simulated LAI and soil moisture of winter wheat at the jointing and heading stages were close to the measured LAI and soil moisture; 2) to verify whether the simulated harvest date of winter wheat was close to the field survey. To verify the calibration accuracy of the CERES-Wheat model, soil moisture at the depth of 20 cm and the LAI from field measured sample sites in 2014 were compared with the model simulations. The results indicated that the root-mean-square error (RMSEs) between the simulated and measured LAIs and between the simulated and measured soil moisture contents were 1.15 m 2 ·m −2 and 0.02 mm 3 ·mm −3 , respectively, and the differences between the simulated and surveyed harvesting dates were all less than 7 days. Considering that the calibration was conducted over different years and that the wheat cultivars differed slightly among fields, these errors were considered acceptable.
Considering that the simulated VTCI of the sampling sites required for this article could not be obtained directly from the CERES-Wheat model, and based on the results of Sun et al. [24], a similar method was used to obtain the simulated VTCI. A linear regression relationship was established between soil moisture simulated with the CERES-Wheat model at 10-day intervals from 0-20 cm depth and the corresponding remotely sensed VTCI at 10-day intervals determined by site latitude and longitude. The linear regression relationship showed a Pearson correlation coefficient of 0.354, p < 0.01, and based on this linear relationship, the simulated VTCI of the sampling sites were obtained from the soil moisture simulated. The LAI simulated by the sampling sites with the step size of days can be obtained directly from the model output, and the maximum value of the simulated LAI data at a 10-day interval in the winter wheat growing season is taken as the simulated LAI at the 10-day intervals.

E. Winter Wheat Planting Area and Yield Data
Accurate winter wheat planting area and distribution are crucial for winter wheat yield estimation. However, reliable, published data sets are lacking. The winter wheat planting area records obtained from the National Bureau of Statistics (http: //www.stats.gov.cn/) indicate that the winter wheat planting area has been relatively stable in the last decade. In addition, the planting area of winter wheat accounts for more than 90% of the total planting area in the Guanzhong Plain, and some studies had also shown that although the use of a generalized cropland mask may affect the yield estimation accuracy, it was feasible to use a generalized cropland mask for regional crop yield estimation [38]. Therefore, the MODIS land cover type product MCD12Q1 was used for the extraction of planting areas. In this article, the International Geosphere Biosphere Programme (IGBP) classification scheme included in the MCD12Q1 product was specifically selected. We used the MODIS reprojection tool to process MCD12Q1, including conversion of remote sensing image format and map projection, and overlaying with the administrative boundary vector of the study area. Then, the results were extracted according to the IGBP classification scheme to obtain the winter wheat growing areas. Yield data were obtained from the official release of the Shaanxi municipal statistical yearbook of each county data.

F. Construction of Time Series Data
Considering that the time series constructed for each variable at each growth stage covered data of 24 counties from 2011 to 2020 and the time series was too long, this article divided the time series into three groups according to the distribution of counties in the Guanzhong Plain. Eight counties were selected as a group without replication, and each group had 10 a of data from 2011 to 2020, so the length of each time series was 80 a. And each group of data had four growth stages of VTCI, LAI, and corresponding yield, respectively, for a total of 9 time series, so there were 27 time series in 3 groups.

A. PF Assimilation Algorithm
The PF algorithm used in this article was based on the implementation of Nagarajan et al. [14] and Bi et al. [12]. The basic implementation process can be divided into two phases; forecast and update. In the forecast stage, the initial state variable x k of the model is perturbed with a given random distributed random noise v i k (i = 1, 2, …, N) to obtain the initial state x i k at time k. N particles propagate through the state equation, and the forecast state x i k+1 is obtained. The state equation for each particle in the filter is defined as follows: where f ( · ) is the state equation, which represents the CERES-Wheat model in this article; u k represents the model drive parameters. Subsequently, the prior state (x i k+1 ) is converted into observations y i k+1 through the measurement operators H( · ) where ε k+1 represents the observed noise sequence at time k+1.
In the update stage, the importance weight of each particle (wî k+1 ) is calculated by using the model predicted state, remote sensing observations, and the initial importance probability den- Next, the weights (wî k+1 ) are normalized to obtain the normalized particle weights w i k+1 . At the same time, considering the problem of sample weight degradation, this article used a less computationally intensive residual resampling method [39] to resample the posterior probability density of states, removing particles with smaller weights and retaining or replicating particles with larger weights. Setting the weight of each particle after normalized resampling is w i k+1,res , and the state estimate value x k+1 at the k + 1 is as follows: where N was set to 200 based on sensitivity analysis [16].

B. Cross-Wavelet Transform
The Morlet wavelet was chosen to meet the requirements of this article for time series analysis. Morlet wavelet is expressed as [40] The corresponding frequency domain wavelet function is where t denotes time, ω 0 and ω represent the dimensionless angular frequency and angular frequency, respectively, s stands for the scaling scale and H(ω) represents the Heavyside function that satisfies the following conditions when ω > 0, H(ω) = 1, The cross-wavelet power spectrum focuses on the interrelationship between the sequences in the high-energy region in the time-frequency domain [41], and reflects the resonance periods between the two time series. The wavelet cross-spectrum W x n y (s) between the time series VTCIs or LAIs at each growth stage of winter wheat and the time series yields is [40] W x n y (s) = W x n (s)W y (s) * (10) where W x n (s) represents the wavelet transform coefficients of the time series VTCIs or LAIs at four growth stages, and W y (s) * stands for the complex conjugate of the wavelet transform coefficient of the time series yields. The larger the crosswavelet power spectrum(|W x n y (s)|), the stronger the correlation between the two series at this time scale, which is used to determine the resonance periods between the time series VTCIs or LAIs at different growth stages of winter wheat and the time series yields.

C. Wavelet Cross-Correlation Analysis
This study utilized wavelet cross-correlation analysis to further quantify the correlation between time series. The wavelet cross-correlation coefficient C W R n (s) between the time series VTCIs or LAIs at each growth stage and the time series yields can be expressed as [42] where σ 2 (W x n (s)) and σ 2 (W y (s)) denote the variance of the time series VTCIs or LAIs and time series yields wavelet coefficients corresponding to scale s, respectively. Cov(W x n , W y n ) s is the covariance of the time series VTCIs or LAIs and time series yields wavelet coefficients at scale s.
Based on this C W R n (s), the wavelet cross-correlation degree of time series VTCIs or LAIs and time series yields was defined as where f (C W R n (s)) represents the weight coefficient of time series VTCIs or LAIs and time series yields on the time scale s.

D. Construction of Regional Yield Estimation Model
Since the selected sampling sites are mainly distributed in winter wheat planting areas of the Guanzhong Plain, it is considered to represent the growth situation of winter wheat in the whole Plain. Therefore, the regional assimilated VTCI and LAI were obtained by inputting regional MODIS-retrieved VTCI and LAI into the linear relationships that the sampling site scale assimilated VTCI and the MODIS-retrieved VTCI, and the sampling site scale assimilated LAI and the MODIS-retrieved LAI, respectively. The cross-wavelet transform was applied to the regional-scale time series assimilated VTCIs and LAIs with time series yields, respectively, and the wavelet cross-correlation degrees were calculated and normalized to obtain the weights of VTCI and LAI at each growth stage of winter wheat. As this article divided the time series into three groups, the weights of the same growth stage in the three groups were averaged to obtain the final weights of VTCI and LAI at each growth stage, thus constructing a bivariate yield estimation model based on assimilation weighting. In this article, the first group would be analyzed as an example.

A. Sampling Site Scale Assimilation 1) Assimilation of VTCI at Different Types of Sampling Sites:
To evaluate the accuracy of the assimilation, specific sampling sites were selected for analysis. Taking the rain-fed sampling site in Shiniu, Qian County, and the irrigated sampling site in Luqiao, Sanyuan County in 2015 as examples, the impact of assimilation on VTCI was analyzed by comparing the changes of observed (retrieved) VTCI, simulated VTCI, and assimilated VTCI, respectively (see Fig. 2). At the rain-fed sampling site [see Fig. 2(a)], from early to mid-March, both the observed VTCIs and the simulated VTCIs showed an upward trend, but the observed VTCI values were significantly higher than those of the simulated. The assimilated VTCIs can comprehensively consider the influence of observations and simulations, and the values were in between. In mid-April, the observed VTCI value was high, and the VTCI value after assimilation decreased due to the influence of the simulated value. At the irrigated sampling site [see Fig. 2(b)], from early to late March, the VTCIs after the assimilation can make corresponding adjustments according to the changes in simulated and observed values, so that they can objectively reflect the changes of VTCI. At the same time, in the process of analysis, it was found that some assimilation results deviated from the simulated and observed values. For example, in Shiniu, Qian County, the assimilated VTCI in early April was significantly high. The reason for this result was that the limited sample size during the use of the assimilation method tended to cause sample depletion in the PF [43], [44]. Overall, the assimilated VTCI can synthesize the influence of the observation values and the simulated values.
To further validate the assimilation of VTCI and considering that VTCI can reflect the effects of water stress during crop growth in real-time, the relationships between observed VTCI, assimilated VTCI and 10-day cumulative precipitation were further analyzed. Here, the analysis was still performed using the two sampling sites from 2015 mentioned previously as examples, as shown in Fig. 3. When the accumulated precipitation in Shiniu, Qian County in early March was less than 10 mm, the observed VTCI was significantly higher, and the assimilated VTCI was significantly decreased, which was more consistent with the changes in precipitation. The cumulative precipitation for early April reached a peak of over 80 mm for the ten days. At this time, the VTCI before and after assimilation was above 0.8, but compared with the observed VTCI, the assimilated VTCI was significantly improved and more realistically reflected the actual precipitation variability. In early March, at the irrigated sampling site in Luqiao, Sanyuan County, the cumulated precipitation was less than 5 mm. The assimilated VTCI was able to respond to changes in the environment on time, and its value was significantly lower. Cumulative precipitation peaked in early May, with VTCI above 0.7 before and after assimilation. Based on the abovementioned analysis, it was found that the assimilated VTCI values were able to respond better to changes in precipitation compared to the observed VTCI values.
2) Assimilation of LAI at Different Types of Sampling Sites: For the analysis of the results of the assimilation of LAI, we still used the sampling sites of two counties in 2015 as examples to compare to the observed (retrieved) and simulated LAIs and analyze the changes in LAI after assimilation, as shown in Fig. 4. At the Shiniu, the rain-fed sampling site [see Fig. 4(a)], the simulated LAI peaked in late April, and the assimilated LAI can comprehensively consider the effects of the observed LAI and the simulated LAI, which can reduce the impact of excessively high simulated value. At the Luqiao, the irrigated sampling site [see Fig. 4(b)], both the simulated LAI and the observed LAI reached their peaks in late April, but the observed value was higher than the simulated value. The assimilated LAI can effectively synthesize the effects of simulated and observed values. To summarize, the assimilated LAI can better take into account the changes in both the observed LAI and the simulated LAI.
The assimilation results were analyzed at individual sampling sites, in order to verify the overall assimilation accuracy, and the errors between assimilated and observed values were calculated for all sampling sites from 2011 to 2020. Specifically, for VTCI, the mean absolute error (MAE) between the assimilated VTCI and the observed VTCI was 0.07, and the RMSE was 0.09. For LAI, the MAE is 0.45 m 2 ·m −2 , and RMSE is 0.59 m 2 ·m −2 . To sum up, it can be concluded that the assimilated VTCIs (LAIs) were closer to the observed VTCIs (LAIs). Based on the abovementioned analysis, the assimilated sampling site scale VTCI and LAI were linearly regressed with the remote sensing observation sampling sites VTCI and LAI, and the coefficients of determination (R 2 ) were 0.78 and 0.77, respectively, to realize the extension from sampling site scale to regional scale, i.e., regional assimilation of VTCI and LAI.

B. Multiple Time Scale Resonance Periods Analysis 1) Analysis of Resonance Periods Between Time Series VTCIs and Time Series Yields:
Based on the Morlet wavelet function, the first group was taken as an example to analyze the resonance periods between the time series assimilated VTCIs at each growth stage and time series yields, respectively. For assimilated VTCIs (see Fig. 5), there was a resonance period of mainly 2-3 a at the green-up stage, which showed a positive correlation with VTCI lagging yield change, indicating that there was a lagging effect of water stress on yield at the green-up stage. At the jointing stage, there were resonance periods of 5-6 a, 10-13 a, and 22-26 a between the time series assimilated VTCIs and the time series yields. At the 10-13 a resonance period, a negative correlation was observed where assimilated VTCI overtook the yield change. It can be understood that at the jointing stage, the effect of water stress on yield was not fully evident in time during the shorter resonance period, and as the resonance period expanded, when the water stress situation eased, i.e., rehydration, it can have different degrees of compensatory effects on dry matter accumulation in different organs of winter wheat, but to a limited extent, which would eventually lead to a negative correlation between the two. In the resonance period of 22-26 a, the assimilated VTCI and yield were positively correlated in the same direction, indicating that assimilated VTCI and yield varied simultaneously in this resonance period. The analysis of the resonance periods between VTCI and yield at the jointing stage showed that the resonance periods between VTCI and yield were complex and had a relatively large impact on yield during this growth stage. At the heading-filling stage, there were mainly resonance periods of 2-3 a, 14-16 a, and 22-26 a between the time series assimilated VTCIs and the time series yields. In the 2-3 a period, assimilated VTCI was negatively correlated with yield change, with assimilated VTCI ahead of yield change. At the 14-16 a period, it showed a positive correlation lagging behind the change in yield. The resonance period was more pronounced at 22-26 a and passed the standard background spectrum test with a confidence level of 95%, showing a positive correlation of VTCI ahead of the change in yield. At the milk maturity stage, assimilated VTCI and yield had resonance periods of 5-6 a, 10-12 a, and 22-26 a. Among them, the resonance period of 22-26 a was more pronounced and passed the standard background spectrum test with a confidence level of 95%. It also showed a positive correlation between VTCI and yield changed simultaneously, i.e., VTCI can characterize yield change at the milk maturity stage in real time. And the resonance periods between VTCI retrieved by MODIS and yield at each growth stage were consistent with those between assimilated VTCI and yield.
The abovementioned analysis showed that the resonance characteristics of the time series assimilated VTCIs and time series yields were more complex at the jointing and heading-filling stages. In contrast, at the green-up and milk maturity stages, the resonance characteristics showed relatively stable periods of temporal resonance with VTCI lagging behind and changing in synchrony with yield, respectively. It showed that the resonance changes between VTCI and yield contained more attributes that influence final yield at the jointing and headingfilling stages, and therefore, VTCI was more important at the jointing and heading-filling stages relative to the green-up and milk maturity stages. The results obtained previously for the relative importance of the various growth stages of VTCI with winter wheat yield were consistent with those of previous research [45].
Note: The thick black solid line enclosed area in the figure indicates that the power spectrum value has passed the standard background spectrum test with a confidence level of 95%, and the thin black solid line is the influence cone of the wavelet boundary effect and its envelope area is the effective value. The arrows in the figure indicate the phase relationship, → indicates that the yield and the assimilated VTCI are in the same phase, representing that the two are positively correlated; ← indicates that the yield and assimilated VTCI are in antiphase, representing that the two are negatively correlated; and indicate the positive and negative correlations of assimilated VTCI lagging yield change, respectively; and and indicate the positive and negative correlations of the assimilated VTCI ahead of yield change, respectively.  I  WEIGHTS OF THE VARIABLES OBTAINED UNDER THE OBSERVED AND ASSIMILATED SOURCES AT EACH GROWTH STAGE BASED ON THE CHARACTERISTIC TIME  SCALES 2) Analysis of Resonance Periods Between Time Series LAIs and Time Series Yields: By analyzing the results of the crosswavelet power spectrum of the time series LAIs at different growth stages and the time series yields, the following results were found. For assimilated LAI, there were resonance periods of 14-18 a and 24-30 a with the yield at the green-up stage. The positive correlation of LAI ahead of yield change passed the standard background spectrum test with a confidence level of 95%, indicating that the change in crop growth at green-up had a significant effect on later growth. At the jointing stage, the assimilated LAI resonated with the yield at 14-18 a and 26-30 a, and both showed positive correlations with LAI ahead of yield change, with the resonance being more significant at 26-30 a. There was a resonance period of 26-30 a between assimilated LAI at the heading-filling and milk maturity stages and yield, and LAI showed a positive correlation with simultaneous change in yield. And the resonance periods between LAI retrieved by MODIS and yield at each growth stage were consistent with those between assimilated LAI and yield. The abovementioned analysis showed that the early growth of LAI had a strong influence on the later growth, with later growth stages of LAI having a greater impact on yield. Therefore, the effect of LAI on yield was greater in the last two growth stages compared to the first two, and this result was in line with the findings by Tian et al. [45].

C. Wavelet Cross-Correlation Analysis 1) Cross-Correlation Analysis of Time Series VTCIs and Time
Series Yields: Wavelet cross-correlation coefficient can reflect the degree of correlation between two time series at different scales in the whole time domain. The wavelet cross-correlation coefficients between the assimilated VTCI at each growth stage and winter wheat yield were obtained using the Morlet wavelet. The results showed that the wavelet cross-correlation between assimilated VTCI and yield alternated with peaks and troughs at each growth stage as the time scale changed and there were negative correlation coefficients between VTCI and yield, indicating that there was a clear periodic characteristic between VTCI and yield at each growth stage, and VTCI showed a negative correlation with yield. According to the definition of VTCI that the smaller the VTCI value, the more severe the drought and the lower the yield. The resonance periods corresponding to the positive correlation between assimilated VTCI at each growth stage and yield were selected as the characteristic time scales for analyzing the relative importance of VTCI at each growth stage to yield. The determined characteristic time scales were 2-3 a at the green-up stage, 5-6 a and 22-26 a at the jointing stage, 14-16 a and 22-26 a at the heading-filling stage, and 5-6 a and 22-26 a at the milk maturity stage.
Wavelet cross-correlation degrees at characteristic time scales were calculated and normalized to obtain the VTCI weights at each growth stage (see Table I) while comparing the weights obtained at each stage to explore the effect of assimilation or not on the analysis of the relationship between VTCI and yield, i.e., a comparative analysis of the weights obtained from MODIS observed (retrieved) VTCI and assimilated VTCI, which were both cross-wavelet transformed. The analysis showed that before and after assimilation, the relative weights of VTCI at each growth stage were consistent, that was, the weight at the jointing stage was the largest, followed by the heading-filling stage, and the weights at the green-up and milk maturity stages were relatively small. The jointing stage is the main stage of the growth of winter wheat roots, stems, and leaves, and the absorption and utilization of water in the soil are the most urgent. When soil water is deficient, it will significantly affect the accumulation rate of dry matter mass, thereby affecting the final growth and yield of winter wheat. At the heading-filling stage, the growth mode of winter wheat is mainly changed from nutritional to reproductive growth. The crop water stress due to inadequate water supply significantly affects the rate of photosynthesis, and reduces the synthesis of starch, protein, and organic matter, resulting in a significant decrease in the grain weight of winter wheat. At the stage of milk maturity, the grain structure of the wheat ear has been formed and shows strong tolerance to certain water deficits. At the green-up stage, the leaf, stem, and root organs of winter wheat grow more slowly and do not accumulate much dry matter mass, so its water deficit has relatively little effect on it. In conclusion, the weight of VTCI at each growth stage determined by the resonance periods was reasonable. On this basis, comparing the weights obtained from assimilated and MODIS observed VTCI after the cross-wavelet transform, it can be seen that the weights of VTCI at the jointing and headingfilling stages increased after assimilation, from 0.277 and 0.268 to 0.279 and 0.274, respectively. It was further demonstrated that assimilation could further enhance the ability to extract the importance of VTCI at the key stages of winter wheat.

2) Cross-correlation Analysis of Time Series LAIs and Time Series Yields:
For LAI, we found that wavelet cross-correlation coefficients between the assimilated LAI and the yield, in some stages, were negative. According to the prior knowledge that LAI of winter wheat is positively correlated with yield within a certain range, resonance periods corresponding to the positive correlation between LAI at each growth stage and yield were selected as the characteristic time scales to analyze the relative importance of LAI to yield at each growth stage. The characteristic time scales determined by the resonance periods were 14-18 a and 24-30 a at the green-up stage, 14-16 a and 26-30 a at the jointing stage, and 26-30 a at the heading-filling and milk maturity stages.
According to the determined time scales, the wavelet crosscorrelation degrees of the LAI of the four growth stages were obtained, which were normalized to obtain the LAI weights of each growth stage (see Table I). It can be seen that the weights of LAI based on the resonance periods were greater at the heading-filling and milk maturity stages than those at the greenup and jointing stages, indicating that LAI has a stronger effect on yield during the heading-filling and milk maturity stages of winter wheat growth, while LAI has a less effect on yield at the green-up and jointing stages. This is because, in contrast to other stages, winter wheat undergoes mainly reproductive growth at the heading-filling and milk maturity stages, which determine the grain weight of winter wheat. Comparing the weights of assimilated and MODIS observed LAI, it can be seen that the weights of assimilated LAI added up to higher weights at the latter two growth stages, which were more closely related to yield. It was further shown that the assimilation method allowed greater weights to be given to variables at growth stages that were more closely related to yield, and enabled the effective extraction of key features.

D. Construction and Application of Yield Estimation Models
Based on the obtained weights of VTCI and LAI at each growth stage, the weighted observed VTCI and LAI, and the weighted assimilated VTCI and LAI were calculated, respectively, and bivariate regressions based on weighted VTCI and LAI were performed on the yield to obtain the estimated yield models. The constructed yield estimation models using observed variables and assimilated variables had normalized root mean square errors (NRMSEs) of 13.29% and 13.22%, R 2 s of 0.50 and 0.50, and mean relative errors (MREs) of 10.68% and 10.58%, respectively. From the accuracy analysis, it was easy to see that the estimated yield model constructed after assimilation had slightly improved accuracy compared to the unassimilated (observed) model, indicating the model based on PF assimilation performed better. To further analyze the performance of the models, the probability distribution of the estimated yields of the above models was plotted and compared with the official yields (see Fig. 6). It can be seen that the estimated yield probability distribution curve of the yield estimation model constructed using assimilated data was a good fit close to the probability distribution curve of the official yields. Specifically, the distribution probability curve of the model constructed by assimilation had a slight shift to the left around 3000 kg/ha, a slight decrease around 4000 kg/ha, and a slight upward shift around 5000 kg/ha compared to the model constructed from MODIS observations. All of these changes resulted in assimilation results that were closer to the probability distribution curve of the official yields. Therefore, we concluded that the bivariate yield estimation model constructed with assimilated VTCI and LAI by PF was superior to the unassimilated (observed) model.

V. DISCUSSION
The growth of wheat is staged; this article focusing on exploring the multiple time scale relationships that exist between variables and yield at different growth stages using the crosswavelet transform. It was found that for VTCI, compared to the lagging behind and changing in synchrony with the yield at the green-up and milk maturity stages, respectively, at the jointing and heading-filling stages, VTCI and yield exhibited more complex and diverse relationships with ahead and lagging as the time scale changes. In terms of LAI, there was a positive correlation between LAI and yield at the heading-filling and milk maturity stages with simultaneous changes, further demonstrating that LAI in the latter two growth stages can characterize yield changes on time. When the weights of each growth stage were obtained by calculating the wavelet cross correlation, it was further found that VTCI at the jointing and filling-heading stages and LAI at the filling-heading and milk maturity stages were more strongly correlated with yield than the other stages. This finding was consistent with Tian's study [45], which focused on the relationships between VTCI and LAI and yield at different growth stages, except that the study used a typical machine learning approach of back propagation neural network. Furthermore, the construction process was more straightforward and simpler. In contrast, this article constructed a yield estimation model based on the investigation of the resonance periods between VTCI and LAI and yield, which made the construction process relatively complicated, but at the same time, it led to the model having better theoretical support and interpretability. In the process of comparison, it was found that machine learning methods, especially deep learning methods, can better extract the nonlinear relationships existing with variables and improve the estimation accuracy, but have the problem of poor model interpretability. Therefore, in the next step of research, we can try to combine wavelet analysis methods, which have the advantage of interpretability, and deep learning methods to improve the model in terms of both accuracy and interpretability.
In this article, a method combining data assimilation and cross-wavelet transform was proposed for the regional winter wheat yield estimation study. First, the data assimilation method is used to combine remote sensing data with the crop growth model, and second, the cross-wavelet transform is used to analyze the periodic variation relationships existing between time series, thus, yield estimation was achieved. It was found that the assimilated remote sensing variables were assigned higher weights during the critical growth stages, thus explicitly demonstrating the ability of assimilation to help achieve feature extraction. However, when comparing the estimation accuracy of the constructed models before and after assimilation, we found that this study suffered from limited improvement in yield estimation accuracy after assimilation. This may be due to the fact that the Guanzhong Plain includes rain-fed and irrigated croplands and when scaling conversions are made from the sampling site scale to the regional scale, the scaling conversion relationships are different for different cropland types. This article only used a unified and simple linear relationship for the scale transformation, which may result in the nonlinear relationship features present in the transformation relationship being ignored, thus weakening the prominence of assimilation for feature extraction. Therefore, in future research, different nonlinear relationships will be tried for different types of land cover to perform the scale conversion process, thus improving the accuracy of the assimilation-based yield estimation model. Also, it is worth noting that remote sensing observations can obtain the spatial distribution of environmental factors associated with crop growth, whereas crop growth models cannot simulate the spatial distribution of environmental factors in the field, thus leading to a spatial mismatch between remote sensing observations and crop growth simulations [7], [21]. Therefore, the existing research can be improved by increasing the spatial density of sampling data and improving the spatial analysis function of the crop model from the perspective of spatial matching.

VI. CONCLUSION
In this article, a multiple time scale approach for regional yield estimation of winter wheat has been developed by combining the crop model with remote sensing through assimilation and crosswavelet transform. Among them, the bivariate assimilation of VTCI and LAI was achieved using the PF algorithm and further enhanced the feature extraction of key growth stages, i.e., VTCI at jointing and heading-filling stages and LAI at heading-filling and milk maturity stages.
When the time series were analyzed using the cross wavelet transform, it was further understood that there were period specific variation relationships between the time series VTCIs, LAIs and yields for different growth stages. The temporal variations between VTCI and yield at the jointing and heading-filling stages were more complex, while the changes at the green-up and milk maturity stages were more stable, i.e., the variation of VTCI lagged and synchronized with the yield variation, respectively. For LAI, it took priority over yield changes in the first two growth stages of winter wheat and changed in parallel with yield in the last two growth stages.
Compared with the yield estimation model constructed from unassimilated variables, the accuracy of the model constructed by assimilated VTCI and LAI has been improved and the distribution of estimated yields was more consistent with the official statistical yields. Therefore, the yield estimation model constructed by combining data assimilation and cross-wavelet transform further clarified the multitime scale variation relationships between variables and yield based on combining remote sensing and crop growth models, which can provide a reference for regional yield estimation. Her research interests include quantitative remote sensing and its application in agriculture.