CN110766308B - A regional crop yield estimation method based on set assimilation strategy - Google Patents

Abstract

The invention discloses a method for estimating regional crop yield based on a set assimilation strategy. The scale simulation yield, as well as the corresponding multiple sets of correlation coefficients r and root mean square error rmse; calculate the fitting index F according to r and rmse, and obtain multiple fitting indices; select the largest fitting index as the optimal result, and Record the corresponding assimilation weight as the optimal assimilation weight H _{op_y-1} in the y-1 year; Repeat steps S2 to S5 for y-2, y-3...y-M years to obtain the corresponding The optimal assimilation weights H _{op_y‑2} , H _{op_y‑3} ......H _{op_y‑M} in each year; expand the assimilation weight value interval H _low ‑H _up , and then use the assimilation weight value interval H _low ‑ H _up conducts simulated production estimates for the forecast year Y.

Description

A regional crop yield estimation method based on set assimilation strategy

technical field

The present invention relates to the technical field of agricultural remote sensing, in particular to a method for estimating regional crop yield based on a set assimilation strategy.

Background technique

Crop models often face the problem of insufficient input data when making regional-scale crop yield estimates. There are often obvious spatial differences in surface characteristics, near-surface environment and crop management measures, and some necessary model input data, such as crop phenology data, meteorological data, crop management information, etc., are recorded at the site scale, which leads to the application of crop models to At the regional scale, it is difficult to obtain enough data to represent the spatial heterogeneity of key factors such as initial conditions, crop parameters, and growth processes. Satellite remote sensing data provides continuous monitoring data of large-scale surface information, which can reflect the spatial continuity and time-series change characteristics of surface information. This advantage effectively complements the weakness of crop model application at the regional scale. Therefore, crop model and remote sensing data assimilation technology has become an important way to improve the accuracy of crop model regional simulation, and it is also a hot direction in the field of crop yield estimation in recent years. one.

In the process of assimilating crop models with remote sensing data, a key factor that has a key impact on the assimilation results is the uncertainty of remote sensing data. The main role of this uncertainty is to quantify the assimilation weights of model simulation values and remote sensing data. Whether the assimilation weight is reasonable or not directly affects the accuracy of the final assimilation result. Once the assimilation weight setting deviates, the assimilation result will also have a significant deviation, which will seriously affect the application effect of the assimilation technology. The best way to quantify the uncertainty of remote sensing data is to compare and analyze the remote sensing data and the surface observation results. However, this method requires surface observation data with temporal and spatial density, and the labor cost and time cost are very high, which is difficult to promote. Existing crop model and remote sensing data assimilation technologies often use simplified methods, such as directly setting the assimilation weight value, to reduce the need for surface observations. This method has achieved some positive research results, but is of limited help for practical assimilation yield estimation applications. The reason is that in practical applications, the assimilation weight needs to exist as a prior parameter to drive the assimilation process, and the existing research does not deal with how to perform an accurate prior on the assimilation weight. In fact, due to the many influencing factors of the assimilation weight and the complex interannual variation, the optimal assimilation weight cannot be directly deduced through historical experience, which greatly reduces the practicability of the assimilation technology. At present, the fact that the optimal assimilation weight cannot be known a priori is a major limitation of the practical application of assimilation production estimation technology. How to break through this limitation, carry out assimilation production estimation under the premise of unknown optimal assimilation weight and obtain production estimation results with sufficient accuracy is still a problem. The technological gap is in urgent need of breakthrough.

Accordingly, new technologies are needed to at least partially address the above-mentioned limitations of the prior art.

SUMMARY OF THE INVENTION

The present invention proposes an assimilation production estimation scheme based on a set assimilation strategy. The production estimation accuracy of the scheme is similar to that based on the optimal assimilation weight, that is, a sufficiently accurate production estimation result can be obtained even when the optimal assimilation weight is unknown, greatly reducing the Improves the practicability of assimilation yield estimation technology.

According to an aspect of the present invention, there is provided a regional crop yield estimation method based on a set assimilation strategy, comprising the following steps:

S1. Obtain remote sensing LAI data of forecast year Y and previous M years in the study area, where M is a natural number greater than 3;

S2. For the Y-1 year, use the crop model and the remote sensing LAI data to perform a spatial difference assimilation operation, and the spatial difference assimilation operation equation is:

是同化时间点k的归一化分析矩阵，

和

分别表示同化时间点k的归一化模拟LAI矩阵与归一化遥感反演LAI矩阵，H是同化权重，取值范围为0.01-1.00；in

is the normalized analysis matrix of assimilation time point k,

and

represent the normalized simulated LAI matrix and the normalized remote sensing inversion LAI matrix of the assimilation time point k, respectively, and H is the assimilation weight, which ranges from 0.01 to 1.00;

Among them, within the value range of 0.01-1.00, H takes the value with the step size S respectively, and in the case of the obtained 1/S different assimilation weight values, the assimilation estimation operation is performed respectively, thereby obtaining the corresponding 1/S grid-scale simulation yield;

S3, based on the corresponding yield record data in the study area, evaluate the accuracy of each of the grid-scale simulated yields, characterize them with the correlation coefficient r and the root mean square error rmse, and obtain a total of 1/S group r and rmse;

S4. Calculate the fitting index F according to r and rmse obtained in S3, and obtain a total of 1/S fitting indices;

The calculation formula of the fitting index F is:

Among them, Fi is the fitting index of the ith assimilation yield estimate, ri and rmse _i represent the correlation coefficient and root mean square error of the _ith assimilation yield estimate, r _max , r _min , rmse _max and rmse _min represent the 1 /The maximum value and minimum value of r in group r and rmse and the maximum value and minimum value of rmse;

S5, select the largest fitting index as the optimal result, and record the corresponding assimilation weight as the optimal assimilation weight H _{op_y-1} in the y-1 year;

S6. Repeat steps S2 to S5 for y-2, y-3...yM years to obtain the optimal assimilation weights H _{op_y-2} , H _{op_y-3} ...H _{op_y} for each year _-M ;

S7. Based on H _{op_y-1} , H _{op_y-2} , H _{op_y-3} ...... the maximum value H _opmax and the minimum value H _opmin in H _{op_y-M} delimit the value range of the expanded assimilation weight, and the expanded assimilation weight The lower limit of the value range is called H _low , the upper limit is called H _up , H _low is the result of expanding the minimum value by 10%, H _up is the result of expanding the maximum value by 10%, and the expansion assimilation weight takes the value The range is not greater than 0.99, not less than 0.01;

S8. Based on the value interval H _low -H _up of the expanded assimilation weight obtained in S7, perform simulated production estimation for the forecast year Y, including: using the value interval H _low - H _up of the expanded assimilation weight to replace the value range of 0.01-1.00 In this case, repeat step S2, expand the spatial difference assimilation operation, and perform (H _up -H _low )/S+1 assimilation operations in total, thereby obtaining the corresponding (H _up -H _low )/S+1 grid scales Simulate the output and take the average value as the simulated output of the forecast year Y.

According to an embodiment of the present invention, in step S1, M is 4, of course, other values can be appropriately selected.

According to an embodiment of the present invention, the crop model is the MCWLA series model, or other crop models can also be selected.

According to an embodiment of the present invention, the step size S is 0.01, or other values, such as 0.02, are selected according to specific conditions.

According to one embodiment of the present invention, the remote sensing LAI data may be selected from Copernicus LAI, GLASS LAI and GLOBMAP LAI, or other remote sensing LAI products.

According to one embodiment of the present invention, the crop is wheat, or other crops such as maize.

According to an embodiment of the present invention, in step S2, the normalized analysis matrix, the normalized simulated LAI matrix and the normalized remote sensing inversion LAI matrix are obtained by a normalization method, so that the original matrix is normalized to take A normalized matrix with values in the range (0,1), represented by the following equations (3) and (4):

表示同化时间点k上的模拟LAI矩阵，D_k代表同一时间点k上对应的遥感LAI矩阵；

和

分别表示由

和D_k转化而来的归一化模拟LAI矩阵和归一化遥感LAI矩阵；S^f和S^D是用来对原矩阵进行归一化的归一化函数。in,

represents the simulated LAI matrix at the assimilation time point k, and D _k represents the corresponding remote sensing LAI matrix at the same time point k;

and

respectively represented by

The normalized analog LAI matrix and the normalized remote sensing LAI matrix transformed from D _k ; S ^f and S ^D are the normalization functions used to normalize the original matrix.

通过归一化函数S^f的反函数

进行反归一化转化以得到分析LAI矩阵

并由

驱动模型向下一同化时间点k+1运行：According to an embodiment of the present invention, in step S2, the differential assimilation further includes:

By the inverse of the normalized function S ^f

Perform an inverse normalization transformation to get the analytical LAI matrix

and by

The drive model runs downward at time point k+1 of the assimilation:

Compared with the prior art, the present invention integrates the idea of collective operation on the basis of the assimilation technology of remote sensing data and crop model. It provides a solution to the practical problem of "how to carry out assimilation estimation and obtain yield estimation results with sufficient accuracy when the optimal assimilation weight cannot be known a priori". Using the assimilation yield estimation scheme based on the set assimilation strategy proposed by the present invention, the yield estimation accuracy is similar to the yield estimation accuracy based on the optimal assimilation weight, that is, a sufficiently accurate yield estimation result can be obtained even when the optimal assimilation weight is unknown. The invention makes up for the deficiencies of the assimilation production estimation technology when it is oriented to business application, and greatly improves the practical application capability of the assimilation production estimation technology.

Description of drawings

The same reference numbers in the figures designate the same or similar parts or parts. Objects and features of the present invention will become more apparent in view of the following description taken in conjunction with the accompanying drawings, in which:

FIG. 1 is a schematic flowchart of a method for estimating regional crop yield based on a set assimilation strategy according to an embodiment of the present invention.

Figure 2 is a schematic diagram of a study area for implementing the method of the present invention according to one embodiment of the present invention;

FIG. 3 is an implementation effect diagram of the method according to the present invention when different remote sensing LAI products are applied.

Detailed ways

In order to clearly illustrate the solutions in the present invention, preferred embodiments are given below and described in detail with reference to the accompanying drawings. The following description is merely exemplary in nature and is not intended to limit the application or uses of the present disclosure

It should be understood that the crop model and remote sensing data referenced in the present invention are known per se, such as various sub-modules of the model, various parameters, operating mechanisms, etc. Therefore, the present invention focuses on the collection between the crop model and the remote sensing data. assimilation process.

Figure 1 is a schematic diagram according to one embodiment of the present invention. Referring to Fig. 1 , the method for estimating yield of regional crops based on a set assimilation strategy according to the present invention includes the following steps:

S1. Obtain the remote sensing LAI data of the predicted year Y in the study area and the previous M years, where M is a natural number greater than 3; the figure shows that M is 4, that is, the y-1, y-2, y-3, The remote sensing LAI data of y-4, wherein the remote sensing LAI data can be selected from Copernicus LAI, GLASS LAI and GLOBMAP LAI, or other appropriate data. That is to say, this method can be applied to a variety of remote sensing LAI data, and the spatial resolution and temporal frequency of the assimilation production estimation operation should be consistent with the temporal and spatial resolution of the remote sensing LAI data.

S2. For the Y-1 year, use the crop model (the crop model is shown as the MCWLA-Wheat model) and the remote sensing LAI data to perform differential assimilation, including performing spatial differential assimilation operations, and the spatial differential assimilation operational equation for:

是同化时间点k的归一化分析矩阵，

和

分别表示同化时间点k的归一化模拟LAI矩阵与归一化遥感反演LAI矩阵，H是同化权重，取值范围为0.01-1.00；应该理解，式中的上标是为了区分矩阵，a是analysis缩写，表示分析矩阵；f是forecast缩写，表示预报矩阵即模型模拟矩阵；n是normalization缩写，表示归一化。in

is the normalized analysis matrix of assimilation time point k,

and

Represents the normalized simulated LAI matrix and the normalized remote sensing inversion LAI matrix of the assimilation time point k, respectively, and H is the assimilation weight, which ranges from 0.01 to 1.00; it should be understood that the superscript in the formula is to distinguish the matrix, a is the abbreviation of analysis, which means the analysis matrix; f is the abbreviation of forecast, which means that the forecast matrix is the model simulation matrix; n is the abbreviation of normalization, which means normalization.

Wherein, within the value range of 0.01-1.00, H takes a value with a step size S (for example, 0.01), and in the case of the obtained 1/S (for example, 100) different assimilation weight values, respectively perform assimilation. According to the calculation of yield estimation, the corresponding 1/S(100) grid-scale simulated yield is obtained.

More specifically, an assimilation process of differential assimilation includes:

1) Normalize the simulated LAI matrix and the corresponding remote sensing LAI matrix at the time point when remote sensing data exists (called assimilation time point) to normalize the original matrix to a normalized value range of (0,1]. The normalization matrix (the following equations (3) and (4)) has many normalization operation methods. For example, it can be obtained by dividing the LAI matrix by the maximum value of its matrix elements:

表示同化时间点k上的模拟LAI矩阵，D_k代表同一时间点k上对应的遥感LAI矩阵；

和

分别表示由

和D_k转化而来的归一化模拟LAI矩阵和归一化遥感LAI矩阵；S^f和S^D是用来对原矩阵进行归一化的归一化函数。in,

represents the simulated LAI matrix at the assimilation time point k, and D _k represents the corresponding remote sensing LAI matrix at the same time point k;

and

respectively represented by

The normalized analog LAI matrix and the normalized remote sensing LAI matrix transformed from D _k ; S ^f and S ^D are the normalization functions used to normalize the original matrix.

2) Use the assimilation equation (1) to assimilate the normalized analog LAI matrix and the normalized remote sensing LAI matrix:

是同化时间点k的归一化分析矩阵，H是同化权重，取值范围为0.01-1.00，在此H以步长0.01分别取值，得到100个不同的同化权重值。in

is the normalized analysis matrix of the assimilation time point k, and H is the assimilation weight, which ranges from 0.01 to 1.00. Here, H takes values with a step size of 0.01, and 100 different assimilation weight values are obtained.

通过归一化函数S^f的反函数

进行反归一化转化以得到分析LAI矩阵

并由

驱动模型向下一同化时间点k+1运行：3)

By the inverse of the normalized function S ^f

Perform an inverse normalization transformation to get the analytical LAI matrix

and by

The drive model runs downward at time point k+1 of the assimilation:

Among them, a complete assimilation and yield estimation calculation process in this step includes: the model starts daily simulation from the sowing period, starts to assimilate with the remote sensing data at the assimilation starting point, the assimilation process continues to the assimilation end point, and the model continues to run until the mature stage to output grid-scale simulated yield .

S3, based on the corresponding yield record data in the study area, the accuracy of each of the grid-scale simulated yields is evaluated, characterized by the correlation coefficient r and the root mean square error rmse, and a total of 1/S groups r and rmse are obtained; the The evaluation of yield estimation accuracy is based on the spatial scale of the actual statistical yield, and the yield estimation result obtained by the assimilation operation is at the grid scale, which needs to be aggregated to the scale of statistical yield, such as the county level.

S4. Calculate the fitting index F according to r and rmse obtained in S3, and obtain a total of 1/S fitting indices;

The calculation formula of the fitting index F is:

Among them, F _i is the fitting index of the ith assimilation yield estimate, ri and rmse _i represent the correlation coefficient and root mean square error of the _ith assimilation yield estimate, r _max , r _min , rmse _max and rmse _min represent the The maximum and minimum values of r in 1/S group r and rmse, as well as the maximum and minimum values of rmse; the larger the F value, the higher the accuracy of the yield estimation result.

S5, select the largest fitting index as the optimal result, and record the corresponding assimilation weight as the optimal assimilation weight H _{op_y-1} in the y-1 year;

S6. Repeat steps S2 to S5 for y-2, y-3...yM years to obtain the optimal assimilation weights H _{op_y-2} , H _{op_y-3} ...H _{op_y} for each year _-M , the figure shows that M is 4, that is, H _{op_y-4} ;

S7. Based on H _{op_y-1} , H _{op_y-2} , H _{op_y-3} ...... the maximum value H _opmax and the minimum value H _opmin in H _{op_y-M} delimit the value range of the expanded assimilation weight, and the expanded assimilation weight The lower limit of the value range is called H _low , the upper limit is called H _up , H _low is the result of expanding the minimum value by 10%, H _up is the result of expanding the maximum value by 10%, and the expansion assimilation weight takes the value The range is not greater than 0.99 and not less than 0.01; that is,

H _low =max(0.9*H _opmin ,0.01) (6)

H _up =min(1.1*H _opmax ,0.99) (7)

S8. Based on the value interval H _low -H _up of the expanded assimilation weight obtained in S7, perform simulated production estimation for the forecast year Y, including: using the value interval H _low - H _up of the expanded assimilation weight to replace the value range of 0.01-1.00 In this case, repeat step S2, expand the spatial difference assimilation operation, and perform a total of (H _up -H _low )/S+1(101) assimilation operations, thereby obtaining the corresponding (H _up -H _low )/S+1 Grid-scale simulated yield (that is, n sets of results in the attached figure), and the average is taken as the simulated yield of the forecast year Y.

Below in conjunction with specific embodiment, further set forth the inventive method:

The specific application of the method of the present invention is exemplified below by taking the estimation of winter wheat yield in the North China Plain as an example.

The central part of the North China Plain was selected as the study area, and the study period was from 2008 to 2015. Using the MCWLA-Wheat model, the remote sensing LAI data used include Copernicus LAI (spatial resolution 1km × 1km, time frequency is every 10 days), GLASS LAI (spatial resolution 1km × 1km, time frequency is every 8 days) ), and GLOBMAP LAI (spatial resolution 0.08°×0.08°, temporal frequency is every 8 days), each remote sensing LAI product was assimilated with the MCWLA-Wheat model to estimate winter wheat yield. The study area is shown in Figure 2.

Based on the set assimilation strategy proposed by the present invention, the implementation process of assimilation production estimation is as follows (taking Copernicus LAI as an example, the implementation process of GLASS LAI and GLOBMAP LAI is consistent with Copernicus LAI):

In this example, 2012-2015 is used as the estimated production period. For each year's estimated yield, the first 4 years are used as the historical period for assimilation yield estimation to provide prior knowledge about the optimal assimilation weight. That is, for the assimilation and yield estimation operation of year y, the first four years (ie y-1, y-2, y-3, y-4) need to be assimilated yield estimation calculation to provide the first year y with the optimal assimilation weight. test knowledge. For example, for the assimilation production estimation in 2012, it is necessary to use the assimilation production estimation operation for 2008 to 2011 to obtain the optimal assimilation weight year by year. The assimilation is performed using the spatial difference assimilation algorithm.

The operation steps of the present invention will be described in detail below with reference to FIG. 1 .

As shown in the figure, in year y-1, the MCWLA-Wheat model and remote sensing LAI data are used to perform spatial difference assimilation operation. It can drive 100 assimilation yield estimation operations and obtain 100 sets of yield estimation results. The complete process of each assimilation yield estimation operation is as follows:

Daily simulations of winter wheat growth using MCWLA-Wheat from the winter wheat sowing date. The starting point of assimilation is the greening period of winter wheat, and the end point of assimilation is the mature period of winter wheat. Between the rejuvenation period and the mature period, the model-simulated LAI and the remote sensing LAI data were assimilated using the spatial difference assimilation algorithm at each time point with remote sensing observations (called assimilation time points). At the mature stage, the assimilation and model simulation process were completed, and the simulated yield of winter wheat was output grid-by-grid.

After an assimilation yield estimation operation is completed, the yield estimation results at the grid scale are aggregated to the county level and compared with the statistical yield, and the yield estimation accuracy is characterized by the correlation coefficient r and the root mean square error rmse. Using different assimilation weights H to repeat 100 assimilation operations, 100 sets of r and rmse can be obtained, and each assimilation weight H corresponds to a set of r and rmse values. On this basis, the fitting index _F of each _{assimilation yield estimation} result is calculated. The maximum value, the minimum value, and the maximum value and minimum value of rmse.

According to formula (2), 100 fitting indices F corresponding to 100 assimilation weights H can be obtained. The largest fitting index is selected as the optimal fitting index, and its corresponding assimilation weight H is used as the optimal assimilation weight for year y-1, named H _{op_y-1} .

The process of calculating Hop_y-1 for y-1 year needs to be repeated in y-2, y-3, and y-4 three years respectively to obtain the optimal assimilation weights H _{op_y-2} , H _{op_y-3} , H _{op_y- 4} .

The maximum value H _opmax and the minimum value H _opmin in _{op_y-1} , H _{op_y-2} , H _{op_y-3} , H _{op_y-4} delimit the value range of the assimilation weight, and the value range is set as the maximum and minimum value ranges outward It is formed by expanding by 10%, and the value range is not greater than 0.99 and not less than 0.01. The lower limit of the interval is called H _low , and the upper limit is called H _up . For details, please refer to the above formulas (6) and (7).

After completing the above calculation, perform a set assimilation calculation of yield estimation for year y. The assimilation weight H used in the calculation is the value range [H _low , H _up ] based on the obtained H _low and H _up , and the value is taken in steps of 0.01. long increment. A total of (H _up -H _low )/0.01+1 assimilation operations were run. The mean of the production estimation results obtained by the collective operation is used as the final production estimation result, and the simulation results are shown in Figure 3.

For different remote sensing LAI data, the production estimation method is the same as the above process, the difference is that when using different remote sensing LAI data, the spatial resolution of the grid scale of the model operation needs to be consistent with the remote sensing LAI data, and the assimilation time point interval needs to be consistent with the remote sensing LAI data. The time frequency of the data remains the same. In addition, when using different remote sensing LAI data, the optimal assimilation weights obtained will change each year.

In order to verify the technical effect of the present invention, adopt the following method to evaluate the implementation effect of the present invention:

First, the optimal assimilation weight calculation method applied to years y-1 to y-4 is also applied to year y, and the following parameters of year y are obtained: including r _max , r _min , rmse _max , rmse _min and the maximum fitting index F, called F _op . Next, compare the results of the set assimilative yield estimate in year y with the actual statistical yield, and calculate its r and rmse, and combine the aforementioned r _max , r _min , rmse _max , and rmse _min to calculate the F value of the set assimilative yield estimate result, which is called F _s . According to F _op and F _s , calculate the accuracy index pc:

The accuracy index pc is used to characterize the percentage of the yield estimation accuracy of the set assimilation compared to the yield estimation accuracy of applying the optimal assimilation weight. The larger the pc, the higher the accuracy.

In addition, in order to compare the accuracy of production estimation, the accuracy index pc of some typical production estimation methods can be calculated. These methods include:

a. Use the optimal assimilation weight H of year y-1 to drive the assimilation estimation operation of year y;

b. Use the mean value of the optimal assimilation weight H of year y-1 to year y-2 to drive the assimilation estimation operation of year y;

c. Use the mean value of the optimal assimilation weight H of year y-1 to year y-3 to drive the assimilation estimation operation of year y;

d. Use the mean value of the optimal assimilation weight H of year y-1 to year y-4 to drive the assimilation estimation operation of year y;

e. Use the optimal assimilation weight H of year y-1 to year y-4 to drive the assimilation yield estimation operation of year y and take the mean of the results as the yield estimation result;

f. Use the maximum and minimum interval of the optimal assimilation weight H from year y-1 to year y-4, namely [H _opmin , H _opmax ], with a step size of 0.01, to drive the assimilation yield estimation operation of year y and take the mean of the results as the yield estimation result;

g. Use the method proposed in the present invention to estimate the production of year y.

Using the above methods a to g, the assimilation production estimation calculation was carried out respectively from 2012 to 2015, and the average value of the accuracy index pc of the four-year production estimation results was used to measure the comprehensive implementation effect of the method. The results are shown in FIG. 3 , which shows that the integrated integration strategy proposed by the present invention shows excellent implementation effects when applying different remote sensing LAI products. For example, for CopernicusLAI, the yield estimation accuracy can reach 98.1% of the yield estimation accuracy when using the optimal assimilation weight. When using GLASS LAI and GLOBMAP LAI, the yield estimation accuracy reaches 90.8% and 85.0% of the yield estimation accuracy using the optimal assimilation weight, respectively. . Relatively better accuracy can be obtained by using remote sensing products with higher resolution, which is closer to the yield estimation effect when using optimal assimilation weights.

In addition, the yield estimation accuracy of the collective assimilation strategy proposed in the present invention is stably superior to other yield estimation methods in assimilation yield estimation using different remote sensing LAI data, and has universality. The above results show that the set assimilation strategy proposed in the present invention can carry out assimilation estimation and obtain yield estimation results with sufficient accuracy without knowing the optimal assimilation weight. The invention makes up for the deficiencies of the assimilation production estimation technology when it is oriented to business application, and greatly improves the practical application capability of the assimilation production estimation technology.

In this paper, specific examples are used to illustrate the principles and implementations of the present invention. The descriptions of the above embodiments are only used to help understand the device and the core idea of the present invention; meanwhile, for those skilled in the art, according to the present invention There will be changes in the specific implementation and application scope. To sum up, the content of this specification should not be construed as a limitation to the present invention.

Claims (8)

1. a regional crop yield estimation method based on a set assimilation strategy, comprising the steps: S1. Obtain the remote sensing LAI data of the forecast year y in the study area and the previous M years, where M is a natural number greater than 3; S2, for the y-1st year, carry out differential assimilation, including using the crop model and the remote sensing LAI data to perform spatial differential assimilation operations, and the spatial differential assimilation operation equation is:

in

is the normalized analysis matrix of assimilation time point k,

and

represent the normalized simulated LAI matrix and the normalized remote sensing inversion LAI matrix of the assimilation time point k, respectively, and H is the assimilation weight, which ranges from 0.01 to 1.00; Among them, within the value range of 0.01-1.00, H takes the value with the step size S respectively, and in the case of the obtained 1/S different assimilation weight values, the assimilation estimation operation is performed respectively, thereby obtaining the corresponding 1/S grid-scale simulation yield; S3, based on the corresponding yield record data in the study area, evaluate the accuracy of each of the grid-scale simulated yields, characterize them with the correlation coefficient r and the root mean square error rmse, and obtain a total of 1/S group r and rmse; S4. Calculate the fitting index F according to r and rmse obtained in S3, and obtain a total of 1/S fitting indices; The calculation formula of the fitting index F is:

Among them, Fi is the fitting index of the ith assimilation yield estimate, ri and rmse _i represent the correlation coefficient and root mean square error of the _ith assimilation yield estimate, r _max , r _min , rmse _max and rmse _min represent the 1 /The maximum value and minimum value of r in group r and rmse and the maximum value and minimum value of rmse; S5, select the largest fitting index as the optimal result, and record the corresponding assimilation weight as the optimal assimilation weight H _{op_y-1} in the y-1 year; S6. Repeat steps S2 to S5 for y-2, y-3...yM years to obtain the optimal assimilation weights H _{op_y-2} , H _{op_y-3} ...H _{op_y} for each year _-M ; S7. Based on H _{op_y-1} , H _{op_y-2} , H _{op_y-3} ...... the maximum value H _opmax and the minimum value H _opmin in H _{op_y-M} delimit the value range of the expanded assimilation weight, and the expanded assimilation weight The lower limit of the value range is called H _low , the upper limit is called H _up , H _low is the result of expanding the minimum value by 10%, H _up is the result of expanding the maximum value by 10%, and the expansion assimilation weight takes the value The range is not greater than 0.99 and not less than 0.01, that is, H _low =max(0.9*H _opmin ,0.01) H _up =min(1.1*H _opmax ,0.99); S8. Based on the value interval H _low -H _up of the expanded assimilation weight obtained in S7, perform a simulated production estimation for the forecast year y, including: using the value interval H _low - H _up of the expanded assimilation weight to replace the value range of 0.01-1.00 In this case, repeat step S2, expand the spatial difference assimilation operation, and perform (H _up -H _low )/S+1 assimilation operations in total, thereby obtaining the corresponding (H _up -H _low )/S+1 grid scales Simulate the yield and take the average as the simulated yield for the forecast year y. 2 . The method for estimating yield of regional crops based on a set assimilation strategy according to claim 1 , wherein, in step S1 , M is 4. 3 . 3. The regional crop yield estimation method based on a set assimilation strategy according to claim 1, wherein the crop model is an MCWLA series model. 4 . The method for estimating regional crop yield based on a set assimilation strategy according to claim 1 , wherein the step size S is 0.01. 5 . 5 . The method for estimating regional crop yield based on a set assimilation strategy according to claim 4 , wherein the remote sensing LAI data is selected from Copernicus LAI, GLASS LAI and GLOBMAP LAI. 6 . 6 . The method for estimating yield of regional crops based on a set assimilation strategy according to claim 1 , wherein the crop is wheat. 7 . 7. the regional crop yield estimation method based on set assimilation strategy according to claim 1, is characterized in that, in step S2, described normalization analysis matrix, normalization simulation LAI matrix and normalization remote sensing inversion LAI matrix Obtained by the normalization method, the original matrix is normalized into a normalized matrix with a value range of (0, 1], which is represented by the following formulas (3) and (4):

in,

represents the simulated LAI matrix at the assimilation time point k, and D _k represents the corresponding remote sensing LAI matrix at the same time point k;

and

respectively represented by

The normalized analog LAI matrix and the normalized remote sensing LAI matrix transformed from D _k ; S ^f and S ^D are the normalization functions used to normalize the original matrix. 8. The regional crop yield estimation method based on set assimilation strategy according to claim 1, is characterized in that, described differential assimilation further comprises:

By the inverse of the normalized function S ^f

Perform an inverse normalization transformation to get the analytical LAI matrix

and by

The drive model runs downward at time point k+1 of the assimilation:

Images