CN104134095B - Crop yield estimation method based on scale transformation and data assimilation - Google Patents
Crop yield estimation method based on scale transformation and data assimilation Download PDFInfo
- Publication number
- CN104134095B CN104134095B CN201410155913.0A CN201410155913A CN104134095B CN 104134095 B CN104134095 B CN 104134095B CN 201410155913 A CN201410155913 A CN 201410155913A CN 104134095 B CN104134095 B CN 104134095B
- Authority
- CN
- China
- Prior art keywords
- crop
- lai
- model
- parameters
- scale
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 230000009466 transformation Effects 0.000 title abstract 2
- 238000006243 chemical reaction Methods 0.000 claims abstract description 29
- 238000009826 distribution Methods 0.000 claims abstract description 17
- 238000004088 simulation Methods 0.000 claims abstract description 15
- 208000036855 Left sided atrial isomerism Diseases 0.000 claims description 135
- 238000001914 filtration Methods 0.000 claims description 24
- 238000004422 calculation algorithm Methods 0.000 claims description 14
- 238000011160 research Methods 0.000 claims description 13
- 238000005457 optimization Methods 0.000 claims description 9
- 239000002689 soil Substances 0.000 claims description 9
- 238000013179 statistical model Methods 0.000 claims description 8
- 238000005259 measurement Methods 0.000 claims description 4
- 230000007261 regionalization Effects 0.000 claims description 4
- 238000007635 classification algorithm Methods 0.000 claims description 3
- 238000003066 decision tree Methods 0.000 claims description 3
- 238000011835 investigation Methods 0.000 claims description 3
- 238000004519 manufacturing process Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 3
- 230000001174 ascending effect Effects 0.000 claims description 2
- 230000003595 spectral effect Effects 0.000 claims description 2
- 238000010276 construction Methods 0.000 claims 1
- 241000209140 Triticum Species 0.000 abstract description 28
- 235000021307 Triticum Nutrition 0.000 abstract description 28
- 230000008901 benefit Effects 0.000 abstract description 4
- 238000012272 crop production Methods 0.000 abstract description 2
- 230000001105 regulatory effect Effects 0.000 abstract 1
- 238000011161 development Methods 0.000 description 6
- 230000018109 developmental process Effects 0.000 description 6
- 230000008569 process Effects 0.000 description 6
- 238000012937 correction Methods 0.000 description 5
- 235000013339 cereals Nutrition 0.000 description 4
- 230000008859 change Effects 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 230000002354 daily effect Effects 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 230000007246 mechanism Effects 0.000 description 3
- 238000001556 precipitation Methods 0.000 description 3
- 230000017260 vegetative to reproductive phase transition of meristem Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000029058 respiratory gaseous exchange Effects 0.000 description 2
- 230000000630 rising effect Effects 0.000 description 2
- 241001420622 Meris Species 0.000 description 1
- 208000027066 STING-associated vasculopathy with onset in infancy Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000008278 dynamic mechanism Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 230000003203 everyday effect Effects 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 235000016709 nutrition Nutrition 0.000 description 1
- 230000035764 nutrition Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 230000029553 photosynthesis Effects 0.000 description 1
- 238000010672 photosynthesis Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000009331 sowing Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002194 synthesizing effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000005068 transpiration Effects 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention belongs to the field of agricultural remote sensing, and relates to a crop yield estimation method based on scale transformation and data assimilation and application of the crop yield estimation method to the guiding of crop production. The method comprises the following concrete steps that: 1, parameters are collected for completing spatialization of a WOFOST crop model; 2, a crop type distribution map and a purity percentage map of the crops to be tested are obtained; 3, TMLAI in a range of 30m is obtained; 4, a secondary scale conversion model is built, and a time sequence scale regulating LAI is generated; 5, the crop model error and the remote sensing observation error in a key phonological period of the crops to be tested are obtained; 6, a four-dimensional variation cost function is built, and an optimized crop model parameter is obtained; and 7, the per unit yield of the crops to be tested in a county region is output. The method provided by the invention has the advantages that the problem of scale mismatching between a remote sensing observation picture element and a crop model simulation unit is solved; the precision of a data assimilation model is improved; and the method is suitable for the crop yield estimation in the county region scale, and is particularly suitable for the winter wheat yield estimation in the county region scale.
Description
Technical Field
The invention belongs to the field of agricultural remote sensing, and particularly relates to a crop yield estimation method based on scale conversion and data assimilation.
Background
The remote sensing inversion parameters are fused to a crop mechanism process model by utilizing a data assimilation technology, and the method is an important way for improving the regional crop growth simulation precision and improving the crop yield estimation precision at present. On the field scale, a crop growth model based on the mechanism processes of crop photosynthesis, respiration, transpiration, nutrition and the like can accurately simulate the continuous evolution of crop objects with the time step length of 'day' and the growth and development conditions and the yield of crops on the single-point scale by means of the internal physical process and the dynamic mechanism. When the object model is applied to regional scale in an expansion mode, uncertainty of the initial conditions, soil parameters, crop parameters and spatial distribution of meteorological compulsive factors in the crop model and difficulty of data acquisition are caused due to the heterogeneity of the earth surface and the near-earth surface environment. The satellite remote sensing has the advantages of space continuity and time dynamic change, and can effectively solve the bottleneck that regional crop parameters are difficult to obtain. However, the remote sensing earth observation is restricted by factors such as satellite space-time resolution, and the like, and the internal process mechanism of crop growth and development and yield formation, individual growth and development conditions and the relation between the individual growth and development conditions and environmental meteorological conditions cannot be really disclosed, which is the advantage of a crop model. The data assimilation technology can realize the advantage complementation of the remote sensing observation and the crop model through coupling, and has great potential for improving the regional crop yield estimation.
The high-time-resolution remote sensing data can capture the growth and development changes and the phenological information of crops, so that the high-time-resolution remote sensing data has very important significance for monitoring and estimating the growth of the crops, however, the high-time-resolution remote sensing satellite data running on the orbit at present are often low in spatial resolution, for example, the spatial resolution of AVHRR, MODIS, MERIS and SPOT prediction is 250m-1km in scale. The low spatial resolution remote sensing data further increases the heterogeneity of pixels, causes mismatching between a remote sensing observation scale and a crop model simulation scale, and greatly limits the precision of a data assimilation model. In addition, the standard MODIS LAI product is developed for all vegetation types on a global scale and is not developed for crop design. Research shows that MODIS LAI has serious underestimation phenomenon in winter wheat area. Therefore, direct assimilation of the standard MODIS LAI resulted in an unrealistic low level of post-assimilation LAI and yield. The geostatistical (variogram) scale correction based method and the pel decomposition downscaling method are two means of dealing with the heterogeneity and scale effect of pels. However, a common disadvantage of both methods is that a series of medium and high spatial resolution data of the crop growth period is required, and modeling heterogeneous pixels is a very complex problem, requiring strict approximation and some a priori knowledge, which is not easy to obtain in practical applications. A potential method is to construct a two-level scale conversion method based on ground sampling LAI, medium-high resolution remote sensing data and low resolution remote sensing data to generate the LAI after scale correction. The assimilation scale modified LAI will greatly improve the estimation of regional winter wheat yield. The method aims to overcome the problem that errors caused by direct assimilation of MODIS LAI products are a key problem to be solved urgently in the practice of estimating production by using a data assimilation model.
Disclosure of Invention
In order to solve the following problems in the prior art: the invention provides a crop yield estimation method based on scale conversion and data assimilation, and aims to construct a scale conversion method in a data assimilation model, assimilate LAI after time sequence scale correction to a crop growth model so as to reduce data assimilation model errors caused by scale mismatching between remote sensing observation and the crop model.
The invention provides a crop yield estimation method based on scale conversion and data assimilation, which comprises the following specific steps of:
s1, collecting meteorological parameters, crop parameters, soil parameters and management parameters in a research area and carrying out single-point scale calibration of a WOFOST crop model according to the meteorological parameters, the crop parameters, the soil parameters and the management parameters; then regionalization is carried out based on meteorological parameters, and the spatialization of the WOFOST crop model is completed;
s2, based on the remote sensing images of multiple time phases, combining crop type sample points obtained by ground investigation and combining the phenological characteristics of crops to be detected, constructing a classification rule, and obtaining a crop type distribution map and a purity percentage map of the crops to be detected;
s3, filtering a time sequence MODIS LAI curve of the growth period of the crop to be detected so as to eliminate the influence of data loss and cloud pollution; constructing a regression statistical model based on the actually-measured LAI of the ground actually-measured sample points and the multi-temporal vegetation index (TM VI) to obtain the TM LAI of the area of 30 meters;
s4, fusing the LAI, multi-temporal TM LAI and MODIS LAI information after filtering of the ground actual measurement sample points, constructing a two-stage scale conversion model, and generating a time series scale adjustment LAI;
s5, calculating standard deviation of remote sensing observation and crop model simulation based on LAI of key phenological period of ground sample point, and obtaining remote sensing observation error and crop model error of key phenological period of crop to be detected;
s6, adjusting LAI and crop model simulation LAI according to the scale obtained in S4, introducing remote sensing observation error and crop model error, constructing a four-dimensional variational cost function, minimizing the cost function by using an optimization algorithm, and obtaining optimized crop model parameters through multiple iterations;
and S7, substituting the optimized crop model parameters obtained in S6 into a crop model, selecting pixels of which the purity of the crop to be detected is more than 60%, operating the crop model by pixel units to simulate the yield, summarizing the crop model to county administrative units, outputting the unit yield of the crop to be detected in counties, and guiding grain production.
Wherein, the meteorological parameters of S1 are parameters such as the highest air temperature, the lowest air temperature, the total radiant quantity, the water vapor pressure, the wind speed and the precipitation.
Wherein, the crop parameters of S1 are crop distribution parameters, climatic characteristics of crops and the like.
The step of performing regionalization based on meteorological parameters and completing spatialization of the crop model in step S1 is to perform spatial position matching on collected parameters by combining remote sensing images, and perform regional parameter calibration on the meteorological parameters and accumulated temperature parameters required by the crop model by using a Kriging interpolation algorithm.
And S2, wherein the multi-temporal remote sensing image is a multi-temporal Landsat TM.
In S2, the obtaining of the crop type distribution map and the percentage of purity of the crop to be measured is to adopt a C5.0 decision tree classification algorithm, combine the phenological characteristics and the spectral characteristics of the crop to construct a classification rule, obtain the crop type distribution map, use a 1km grid to fit to generate a 1km spatial distribution map of percentage of purity of the crop to be measured, and select a pixel with a purity of the crop to be measured greater than 60% as a pixel for estimating the scale yield of the subsequent county area in order to reduce the influence of the non-crop area to be measured on the data assimilation result.
Wherein, the filtering is performed as described in S3, and an upper envelope (Savitzky-Golay, S-G) filtering algorithm is adopted, and the expression is shown in formula (1):
in the formula: y isj+iRepresenting the value in a window on the original LAI curve, wherein m is the radius of the window, N is the convolution number, and the width of the window is 2m + 1; y isj *Represents a filtered LAI value; c denotes a filter coefficient of the ith LAI value.
And S3, performing filtering, forming an upper envelope filtering method on the basis of improvement of S-G filtering, dividing the LAI value into a true value and a false value through the LAI time sequence variation trend, replacing the false value with the S-G filtering value in a local loop iteration mode, synthesizing a new LAI smooth curve with the true value, and gradually fitting to form an upper envelope of the LAI time sequence data.
Wherein, the step of constructing the two-level scale conversion model in S4 comprises 2 key steps:
1) selecting TM LAIs in 3 key growth periods from TM LAIs of crops to be detected in the 30 m generation area of S3;
2) adjusting other key growth periods S-G MODIS LAI by utilizing the ratio coefficient between the selected key growth period TM LAI and the S-G MODIS LAI of the corresponding period, respectively fitting by utilizing Logistic curves at the ascending stage and the descending stage of the growth period of the crop to be detected,
the Logistic curve equation is shown in formula (2):
where t is the index of the MODIS LAI time series, y (t) is the LAI value corresponding to t time, a and b are fitting parameters, c + d is the maximum LAI value, and d is the LAI initial value, i.e. the first value in the LAI time series.
S5, calculating the remote sensing observation error of the key phenological period of the crop to be detected based on the standard variance between the actual measurement LAI and the scale correction LAI; the crop model error of S5 is derived from the optimized 2 model parameters, TDWI and SPAN.
Wherein, in S6, the minimization of the cost function by using the optimization algorithm minimizes a four-dimensional variational cost function by using SCE _ UA algorithm, and the four-dimensional variational cost function is calculated according to formula (3):
wherein k represents the number of optimization model parameters in the cost function, xkValues, x, representing the optimized WOFOST parameter over a range of intervalsk0An empirical value representing the number of optimized WOFOST parameters, B represents a model error represented by an error of 2 optimized model parameters, T represents a transposition of a matrix, N represents the number of times of time series remote sensing observation data, yiRepresenting telemetric observations LAI, Hi(x) Representing WOFOST model simulation LAI, Q0Representing the error of remotely sensed LAI.
The S7 specific steps are as follows: and selecting pixels with the crop purity of more than 60% to be detected, reinitializing 2 crop model parameters through continuous iteration to enable LAI and yield output by the crop model to be continuously changed, and then, under the condition that the minimum difference between MODIS LAI and crop model output LAI is compared in a cost function, finishing assimilation when the following three convergence conditions meet one of the three convergence conditions. Obtaining the value of the optimized parameter:
firstly, after 5 continuous cycles, the parameter value to be optimized is contracted to a specified value range;
secondly, the objective function value cannot be improved by 0.0001 percent after 5 cycles;
thirdly, calculating the times of the cost function to exceed 10000 times;
and substituting the optimized model parameters into a crop model to operate to obtain the yield, summarizing according to administrative boundaries of county areas by combining the distribution percentage of the crops to be detected of the pixels, and outputting a unit yield result on the area.
Wherein, the crop to be detected is preferably winter wheat.
The invention also provides application of the crop yield estimation method based on scale conversion and data assimilation in guiding crop production.
Compared with the prior art, the invention has the beneficial effects that:
the method solves the problem of scale mismatching between the remote sensing observation pixel and the crop model simulation unit, improves the precision of the data assimilation model, and is suitable for winter wheat yield estimation of county scale.
Detailed Description
The following examples are given to further illustrate the embodiments of the present invention. The following examples are intended to illustrate the invention but are not intended to limit the scope of the invention.
Example 1
The technical scheme of the invention is further illustrated by taking the estimation of the yield of the winter wheat in Hebei Baoding and constant water areas as an example. The flow of the regional winter wheat yield estimation method based on scale conversion and data assimilation of the embodiment comprises the following steps:
s1, collecting meteorological parameters, crop parameters, soil parameters and management parameters in a research area and carrying out single-point scale calibration of a WOFOST crop model according to the meteorological parameters, the crop parameters, the soil parameters and the management parameters; and then regionalizing based on the meteorological parameters to complete the spatialization of the crop model. The method comprises the following specific steps:
the area for fixing and balancing water in Hebei province is selected as a research area, and the area is an important winter wheat producing area in Hebei province. The following data were obtained: selecting 6 meteorological elements required by the highest/lowest air temperature, total radiation, water vapor pressure, wind speed and precipitation model of 6 meteorological sites; soil parameters and crop parameters collected by an agricultural meteorological test station in a research area; statistics data of winter wheat yield in Hebei province and county in 2009; MODIS15A2 data product, unifying the coordinate projection of all data as Albers ConicalEqualarea, space referenced as WGS 1984.
For the crop parameters which are not sensitive in the WOFOST model, or the parameters which are sensitive to a higher degree but have a smaller variation range, the parameters are determined by looking up the literature or directly adopting the default values of the WOFOST model (such as optimal development light length DLO, lowest emergence temperature lower limit TBASE, most used RML for maintaining respiration of leaves, and the like); the crop parameters with high sensitivity and wide variation range can be obtained by consulting relevant documents in the research area or calculating by using measured data (such as accumulated temperature TSUM1, TSUM2, leaf assimilation conversion rate CVL, specific leaf area SLATB and the like). The possible value range and value of the parameters (such as SPAN in the life SPAN of the blade at 35 ℃, AMAXTB in the maximum CO2 assimilation rate, EFFTB in the light energy utilization rate of single-blade assimilation CO2 and the like) can also be determined by the priori knowledge of the research area.
And (3) interpolating six elements of the highest air temperature, the lowest air temperature, the total radiant quantity, the vapor pressure, the wind speed and the precipitation quantity every day in 2008-2009 by utilizing the observation data of 6 meteorological stations in the research area and combining a Krigin space interpolation method, and converting the six elements into a standard format of a WOFOST model to obtain the meteorological data corresponding to each grid point in the research area. For parameters which are easy to obtain in the spatial range, such as accumulated temperature from a sowing period to a seedling emergence period, accumulated temperature from the seedling emergence period to a flowering period, accumulated temperature from the flowering period to a mature period and the like, daily accumulated temperature is obtained by subtracting biological zero temperature of winter wheat from daily average temperature obtained by a meteorological observation station, the three accumulated temperature variables are obtained by accumulating the daily accumulated temperature at different climatic stages, and the accumulated temperature parameter of spatial variation in the whole research area is obtained by adopting a kriging interpolation method which is the same as meteorological data interpolation.
S2, using multi-temporal U.S. special cartographer (TM) images of a research area, combining crop type sample points obtained by ground investigation, adopting a C5.0 decision tree classification algorithm, combining the phenological characteristics of crops, constructing a classification rule, obtaining a crop type distribution diagram, matching the crop type distribution diagram with a 1km grid to generate a 1km winter wheat purity percentage spatial distribution diagram, and selecting pixels with winter wheat purity of more than 60% as a basis for subsequent county scale yield estimation in order to reduce the influence of non-winter wheat regions on data assimilation results.
S3 time-series data of MODIS LAI over the entire growth period of the experimental area were synthesized, and a time-series curve was generated for each grid element. Specifically, by using an MRT projection conversion tool provided by NASA, the data in the Hebei area are subjected to mosaic, projection conversion and format conversion, 1 to 177 days (julian days) are selected, and 23 time phase images are superposed to generate a time sequence curve based on a grid unit. And filtering the time sequence MODIS LAI curve to eliminate data loss caused by cloud pollution. The problem of discontinuous change process of the leaf area index is solved by using upper envelope (Savitzky-Golay, S-G) filtering, and the principle is expressed as follows:
in the formula: y isj+iRepresenting the value in a window on the original LAI curve, wherein m is the radius of the window, N is the convolution number, and the width of the window is 2m + 1; y isj *Represents a filtered LAI value; c denotes a filter coefficient of the ith LAI value.
On the basis of improving S-G filtering, an upper envelope filtering method is formed, an LAI value is divided into a true value and a false value through an LAI time sequence variation trend, the S-G filtering value replaces the false value in a local loop iteration mode, a new LAI smooth curve is synthesized with the true value, and an upper envelope of LAI time sequence data is formed through gradual fitting. The method comprises the following specific steps:
(1) and carrying out linear interpolation on the LAI value points with cloud pollution. Due to the influence of cloud coverage or atmospheric effect, some pixel values have larger noise, and the pixel values adopt the peripheral credible pixel values to carry out spatial linear interpolation to generate initial LAI data.
(2) And extracting the LAI time sequence variation trend by S-G filtering. According to the hypothesis, LAI time series data should follow the progressive nature of vegetation phenological dynamics, where abrupt changes are thought to be due to cloud or atmospheric noise. Therefore, S-G filtering is applied to each winter wheat pixel in the research area for smoothing, and the change trend data LAItr of the LAI time sequence is obtained.
(3) And determining the weight Wi of each value point in the LAI time sequence, and providing a judgment basis for a subsequent reconstruction result. The formula for Wi is:
wherein,dmaxis the maximum value of di.
(4) A new LAI time series is generated. From the above theory, the LAI value with noise is smaller than the LAItr value on the time variation trend curve, and a simple substitution of these points generates a new LAI1 time sequence curve.
(5) The LAI1 is S-G filtered and the filtered data is compared to the spatially interpolated LAI value. Similarly, if the LAI value after the spatial interpolation is smaller than the newly generated LAI j, the LAI j is used for replacing the LAI j to generate a new sequence LAI j +1, the filtering reconstruction is carried out iteratively, and the filtering effect is FkTo assess the following:
LAIik is an LAI value of the ith pixel after the kth filtering reconstruction, Wi is the weight calculated in the step (3), and Fk is an iteration termination judgment parameter after the kth filtering. The condition for exiting iteration is Fk-1≥Fk≤Fk+1。
The method comprises the steps of obtaining cloud-free TM data of 3 scenes in an experimental area, wherein the time phases are as follows: 2009-3-14, 2009-5-17 and 2009-6-10, and respectively constructing a regression statistical model with the corresponding actually measured LAI.
3, month and 14 days 2009: because winter wheat is in the striking-back stage in the middle of 3 months and the influence of bare soil on the vegetation index is large, the SAVI vegetation index is selected to establish a statistical model between LAI and TM VI:
determining the coefficient R2=0.849。
In 2009, 5, 17 days, the statistical model between the actually measured LAI and TM NDVI was established as:
determining the coefficient R2=0.742。
In 2009, 6, month 10, the statistical model between the actually measured LAI and TM NDVI was established as:
determining the coefficient R2=0.874。
Based on the above 3 regression statistical models, a spatial distribution map of 30 m winter wheat in the study area was obtained.
S4, fusing the LAI, multi-temporal TM LAI and the filtered MODISLAI information of the ground actual measurement sample point, constructing a two-stage scale conversion model, and generating a time sequence scale adjustment LAI.
The second-order scale conversion model comprises 2 steps, firstly, regression statistical models of actually measured LAI and TM VI of 3 phenological periods generated in the step S3 are used for generating a regional 30-meter winter wheat TMLAI; secondly, adjusting other key growth periods S-G MODIS LAI by utilizing a ratio coefficient between the existing key growth period TM LAI and the S-G MODIS LAI of the corresponding period, and obtaining adjusted LAI (a booting period, an elongation period, a rising period and a green returning period) of 4 phenological periods in the rising stage (green returning to the elongation period) of the growth period of the winter wheat; in the descent stage (jointing to maturity) of winter wheat, adjusted LAIs (booting stage, heading stage, flowering stage, maturity stage) of 4 phenological stages were obtained, and fitting was performed using Logistic curves, respectively.
The Logistic curve equation is as follows:
where t is the index of the MODIS LAI time series, y (t) is the LAI value corresponding to t time, a and b are fitting parameters, c + d is the maximum LAI value, and d is the LAI initial value, i.e. the first value in the LAI time series.
S5, calculating standard deviation of remote sensing observation and crop model simulation based on LAI of key phenological period of ground sample point, and obtaining remote sensing observation error and crop model error of key phenological period of winter wheat.
The model error is considered to be derived from 2 optimized parameters (TDWI and SPAN), and the observation error of the scale-corrected LAI in 7 phenological periods is 0.34 through analysis and calculation; the model error is set to 0.52.
S6, adjusting LAI and crop model simulation LAI according to the scale obtained in S4, introducing remote sensing observation error and crop model error, constructing a four-dimensional variational cost function, minimizing the cost function by using an optimization algorithm, and obtaining optimized crop model parameters through multiple iterations.
On the basis of crop model calibration in S1, in an experimental area, operating WOFOST crop model simulation LAI, combining with scale correction LAI generated in S4, introducing remote sensing observation and model simulation error in S5, and constructing a four-dimensional variational cost function as follows:
wherein k represents the number of optimization model parameters in the cost function, xkValues, x, representing the optimized WOFOST parameter over a range of intervalsk0Empirical value representing the number of optimized WOFOST parameters, B representsModel error is expressed by 2 errors of optimized model parameters, T represents transposition of matrix, N represents the number of times of time series remote sensing observation data, yiRepresenting telemetric observations LAI, Hi(x) Representing WOFOST model simulation LAI, Q0Representing the error of remotely sensed LAI.
Selecting pixels with the purity of winter wheat more than 60%, continuously iterating and reinitializing 2 WOFOST model parameters by using an optimization algorithm to change LAI and output of the model, and then obtaining optimized parameters under the condition of comparing the minimum difference between MODIS LAI and WOFOST model output LAI in a cost function. The SCE _ UA algorithm is used to search the initial parameter space for a globally optimal solution, so that the cost function converges quickly. And ending assimilation when the following three convergence conditions are met, and obtaining the numerical value of the optimized parameter.
(1) After 5 continuous cycles, the parameter value to be optimized is contracted to the designated value range;
(2) the objective function value can not be improved by 0.0001% after 5 cycles;
(3) the number of times the cost function is calculated exceeds 10000 times.
And S7, substituting the optimized crop model parameters obtained in S6 into the crop model, selecting pixels with the winter wheat purity of more than 60%, operating the crop model by pixel units to simulate the yield, summarizing to county administrative units, and outputting the county winter wheat unit yield. And substituting the optimized TDWI and SPAN parameters into a WOFOST model to operate, and obtaining the yield. And combining the percentage of the winter wheat of the pixels, summarizing according to administrative boundaries of counties, and outputting a single-yield result of the winter wheat on the region.
According to the method, the winter wheat growth period time sequence MODISLAI curve is filtered by using an S-G filtering algorithm, the influence of data loss and cloud pollution is eliminated, a smooth LAI time sequence curve is generated, and a high-precision data source is provided for a two-level scale conversion model. And then, a time series scale adjustment LAI is generated through a two-stage scale conversion model, the trend information of the MODIS LAI time series and the TM LAI information are effectively fused, the LAI precision is improved, and a high-precision data source is provided for the assimilation process. And finally, the adjusted LAI and the LAI simulated by the crop model are assimilated, so that the precision of the yield estimation of the winter wheat in the region is improved.
According to the method, yield data can be obtained in a large area in the winter wheat mature period, important scientific bases are provided for scientific decisions such as grain condition judgment, grain regulation and control and the like of relevant departments of China, and the method can be used as an important basis for grain trade.
The method of the invention can also be used for the estimation of the yield of other crops.
Although the invention has been described in detail hereinabove with respect to a general description and specific embodiments thereof, it will be apparent to those skilled in the art that modifications or improvements may be made thereto based on the invention. Accordingly, such modifications and improvements are intended to be within the scope of the invention as claimed.
Claims (8)
1. A crop yield estimation method based on scale conversion and data assimilation is characterized by comprising the following specific steps:
s1, collecting meteorological parameters, crop parameters, soil parameters and management parameters in a research area and carrying out single-point scale calibration of a WOFOST crop model according to the meteorological parameters, the crop parameters, the soil parameters and the management parameters; then regionalization is carried out based on meteorological parameters, and the spatialization of the WOFOST crop model is completed;
s2, based on the remote sensing images of multiple time phases, combining crop type sample points obtained by ground investigation and combining the phenological characteristics of crops to be detected, constructing a classification rule, and obtaining a crop type distribution map and a purity percentage map of the crops to be detected;
s3, filtering a time sequence MODIS LAI curve of the growth period of the crop to be detected so as to eliminate the influence of data loss and cloud pollution; constructing a regression statistical model based on the actually-measured LAI of the ground actually-measured sample points and the multi-temporal vegetation index TM VI to obtain the TM LAI of the area of 30 meters;
s4, fusing the LAI and TM LAI of the ground actual measurement sample points and the MODIS LAI information after filtering, constructing a two-level scale conversion model, and generating a time sequence scale adjustment LAI; the construction of the two-level scale conversion model comprises 2 key steps:
1) selecting TM LAIs in 3 key growth periods from TM LAIs of crops to be detected in the 30 m generation area of S3;
2) adjusting other key growth periods S-G MODIS LAI by utilizing a ratio coefficient between the selected key growth period TM LAI and the S-G MODIS LAI of the corresponding period, and respectively utilizing a Logistic curve to carry out fitting in an ascending stage and a descending stage of the growth period of the crop to be detected, wherein the Logistic curve equation is shown in a formula (2):
wherein t is an index of the MODIS LAI time sequence, y (t) is an LAI value corresponding to t time, a and b are fitting parameters, c + d is a maximum LAI value, d is an LAI initial value, namely a first value in the LAI time sequence;
s5, calculating standard deviation of remote sensing observation and crop model simulation based on LAI of key phenological period of ground sample point, and obtaining remote sensing observation error and crop model error of key phenological period of crop to be detected;
s6, adjusting LAI according to the time sequence scale obtained by crop model simulation LAI and S4, introducing remote sensing observation error and crop model error, constructing a four-dimensional variational cost function, minimizing the cost function by using an optimization algorithm, and obtaining optimized crop model parameters through multiple iterations;
and S7, substituting the optimized crop model parameters obtained in S6 into a crop model, selecting pixels of which the purity of the crop to be detected is more than 60%, operating the crop model by pixel units to simulate the yield, summarizing the crop model to county administrative units, outputting the unit yield of the crop to be detected in counties, and guiding grain production.
2. The method for estimating crop yield based on scale conversion and data assimilation of claim 1, wherein the step of performing regionalization based on meteorological parameters and the spatialization of the crop model in step S1 is performed by matching spatial positions of collected parameters with remote sensing images and performing regional calibration of the meteorological parameters and accumulated temperature parameters required by the crop model by using a Kriging interpolation algorithm.
3. The method for crop yield estimation based on scale conversion and data assimilation of claim 1, wherein the multi-temporal remote sensing image of S2 is a multi-temporal Landsat TM.
4. The method for estimating crop yield based on scale conversion and data assimilation of claim 1, wherein the step of obtaining the crop type distribution map and the percentage of crop purity to be tested is performed at S2, a C5.0 decision tree classification algorithm is adopted, a classification rule is constructed by combining the phenological characteristics and the spectral characteristics of the crops, the crop type distribution map is obtained, a 1km spatial distribution map of the percentage of crop purity to be tested is generated by using a 1km grid, and in order to reduce the influence of non-to-be-tested crop areas on the data assimilation result, the pixel with the purity of the to-be-tested crop higher than 60% is selected as the pixel for estimating the yield of the scale in the subsequent county area.
5. The method for estimating crop yield based on scale conversion and data assimilation of claim 1, characterized in that the filtering of S3 is performed by using the upper envelope (Savitzky-Golay, S-G) filtering algorithm, which is expressed by formula (1):
in the formula: y isj+iRepresenting the value in a window on the original LAI curve, wherein m is the radius of the window, N is the convolution number, and the width of the window is 2m + 1; y isj *Represents a filtered LAI value; c denotes a filter coefficient of the ith LAI value.
6. The method for estimating crop yield based on scale conversion and data assimilation of claim 1, characterized in that the remote sensing observation error of the key phenological period of the crop to be tested is calculated based on the standard deviation between the measured LAI and the scale-corrected LAI at S5; the crop model error of S5 is derived from the optimized 2 model parameters, TDWI and SPAN.
7. The method for crop yield estimation based on scale conversion and data assimilation of claim 1, characterized in that the minimization of the cost function using the optimization algorithm is performed at S6
The SCE _ UA algorithm minimizes a four-dimensional variational cost function, which is calculated according to equation (3):
wherein k represents the number of optimization model parameters in the cost function, xkValues, x, representing the optimized WOFOST parameter over a range of intervalsk0An empirical value representing the number of optimized WOFOST parameters, B represents a model error represented by an error of 2 optimized model parameters, T represents a transposition of a matrix, N represents the number of times of time series remote sensing observation data, yiRepresenting telemetric observations LAI, Hi(x) Representing WOFOST model simulation LAI, Q0Representing the error of remotely sensed LAI.
8. The method for estimating crop yield based on scale conversion and data assimilation of claim 1, characterized by the specific steps of S7: selecting pixels with crop purity higher than 60% to be detected, reinitializing 2 crop model parameters through continuous iteration to enable LAI and yield output by a crop model to be changed continuously, and then ending assimilation when the following three convergence conditions are met under the condition that the minimum difference between MODIS LAI and crop model output LAI is compared in a cost function, so as to obtain a numerical value of an optimized parameter:
firstly, after 5 continuous cycles, the parameter value to be optimized is contracted to a specified value range;
secondly, the objective function value cannot be improved by 0.0001 percent after 5 cycles;
thirdly, calculating the times of the cost function to exceed 10000 times;
and substituting the optimized model parameters into a crop model to operate to obtain the yield, summarizing according to administrative boundaries of county areas by combining the distribution percentage of the crops to be detected of the pixels, and outputting a unit yield result on the area.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410155913.0A CN104134095B (en) | 2014-04-17 | 2014-04-17 | Crop yield estimation method based on scale transformation and data assimilation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410155913.0A CN104134095B (en) | 2014-04-17 | 2014-04-17 | Crop yield estimation method based on scale transformation and data assimilation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104134095A CN104134095A (en) | 2014-11-05 |
| CN104134095B true CN104134095B (en) | 2017-02-15 |
Family
ID=51806768
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410155913.0A Active CN104134095B (en) | 2014-04-17 | 2014-04-17 | Crop yield estimation method based on scale transformation and data assimilation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104134095B (en) |
Families Citing this family (28)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106780323B (en) * | 2016-11-23 | 2020-03-17 | 湖北大学 | Agricultural condition acquisition and real-time updating method and system based on smart phone |
| CN107122739B (en) * | 2017-01-23 | 2020-07-17 | 东北农业大学 | Crop yield estimation model for reconstructing VI time series curve based on Extreme mathematical model |
| CN106951979A (en) * | 2017-02-20 | 2017-07-14 | 中国农业大学 | The crop maturity phase Forecasting Methodology that remote sensing, crop modeling are merged with weather forecast |
| CN107016466B (en) * | 2017-04-11 | 2020-12-22 | 北京林业大学 | Landscape pattern optimization construction method and system |
| CN107274036B (en) * | 2017-07-25 | 2020-04-03 | 中国农业科学院农业信息研究所 | Crop yield prediction method and system |
| CN108509836B (en) * | 2018-01-29 | 2021-10-08 | 中国农业大学 | Crop yield estimation method based on double-polarized synthetic aperture radar and crop model data assimilation |
| CN108304973A (en) * | 2018-02-11 | 2018-07-20 | 中国农业大学 | Area crops maturity period prediction technique based on accumulated temperature, radiation and soil moisture content |
| CN108133298B (en) * | 2018-03-08 | 2022-04-19 | 河南工业大学 | National grain consumption prediction method based on multiple regression model |
| CN108982369B (en) * | 2018-04-28 | 2020-09-01 | 中国农业大学 | Plot scale crop growth monitoring method integrating GF-1WFV and MODIS data |
| CN109033539A (en) * | 2018-07-02 | 2018-12-18 | 河北省科学院地理科学研究所 | The calculation method influenced based on plant growth mechanism model separation key factor pair crop phenology |
| CN108921351A (en) * | 2018-07-06 | 2018-11-30 | 北京兴农丰华科技有限公司 | Crop production forecast method based on trend yield and Meteorological Output |
| CN108984803B (en) * | 2018-10-22 | 2020-09-22 | 北京师范大学 | Method and system for spatializing crop yield |
| CN109800921B (en) * | 2019-01-30 | 2019-12-20 | 北京师范大学 | Regional winter wheat yield estimation method based on remote sensing phenological assimilation and particle swarm optimization |
| CN110008621B (en) * | 2019-04-15 | 2023-01-10 | 中国农业科学院农业资源与农业区划研究所 | Crop model remote sensing assimilation estimation method based on dual stream dependence set square root filtering assimilation algorithm |
| CN110378521B (en) * | 2019-06-28 | 2022-04-22 | 河南农业大学 | Construction and application of yield prediction model of winter wheat in northeast China of Henan province |
| CN110633841B (en) * | 2019-08-13 | 2022-04-01 | 中国农业大学 | Provincial range plot scale data assimilation yield prediction method based on set sampling |
| CN110766308B (en) * | 2019-10-17 | 2020-07-10 | 中国科学院地理科学与资源研究所 | Regional crop yield estimation method based on set assimilation strategy |
| CN110751094B (en) * | 2019-10-21 | 2020-08-28 | 北京师范大学 | Crop yield estimation method based on GEE comprehensive remote sensing image and deep learning method |
| CN113033262B (en) * | 2019-12-25 | 2022-08-05 | 中移(成都)信息通信科技有限公司 | Model training method and crop yield estimation method |
| US11935242B2 (en) | 2020-03-09 | 2024-03-19 | International Business Machines Corporation | Crop yield estimation |
| CN111915096B (en) * | 2020-08-14 | 2021-03-09 | 中国科学院地理科学与资源研究所 | Crop yield early-stage forecasting technology based on crop model, remote sensing data and climate forecasting information |
| CN112052988B (en) * | 2020-08-18 | 2024-03-22 | 中国农业大学 | Crop yield estimation method coupling multi-objective optimization and collection assimilation and application |
| CN115018105A (en) * | 2021-03-03 | 2022-09-06 | 中国农业大学 | Winter wheat meteorological yield prediction method and system |
| CN114528282B (en) * | 2022-01-20 | 2024-07-19 | 重庆邮电大学 | Time-space Kriging improvement-based method for reconstructing MODIS NDVI time sequence in cloudy region |
| CN114529097B (en) * | 2022-02-26 | 2023-01-24 | 黑龙江八一农垦大学 | Multi-scale crop phenological period remote sensing dimensionality reduction prediction method |
| CN115753625B (en) * | 2022-11-02 | 2024-04-19 | 中国农业大学 | Yield estimation method and device for regional crops, electronic equipment and storage medium |
| CN116432859B (en) * | 2023-05-08 | 2023-10-27 | 中山大学 | Crop yield statistical data downscaling method |
| CN117830860A (en) * | 2024-03-06 | 2024-04-05 | 江苏省基础地理信息中心 | Remote sensing automatic extraction method of winter wheat planting structure |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7702597B2 (en) * | 2004-04-20 | 2010-04-20 | George Mason Intellectual Properties, Inc. | Crop yield prediction using piecewise linear regression with a break point and weather and agricultural parameters |
| CN102323987A (en) * | 2011-09-09 | 2012-01-18 | 北京农业信息技术研究中心 | Crop leaf area index assimilation method |
| CN102651096A (en) * | 2012-04-28 | 2012-08-29 | 中国农业大学 | Method for estimating yield of winter wheat by assimilating characteristics of leaf area index time-sequence curve |
| JP2012203875A (en) * | 2011-03-28 | 2012-10-22 | Tokyo Electric Power Co Inc:The | Yield estimation device and computer program |
| CN103345707A (en) * | 2013-06-04 | 2013-10-09 | 中国科学院遥感与数字地球研究所 | Crop maturation stage remote sensing prediction method based on multi-source remote sensing data |
-
2014
- 2014-04-17 CN CN201410155913.0A patent/CN104134095B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7702597B2 (en) * | 2004-04-20 | 2010-04-20 | George Mason Intellectual Properties, Inc. | Crop yield prediction using piecewise linear regression with a break point and weather and agricultural parameters |
| JP2012203875A (en) * | 2011-03-28 | 2012-10-22 | Tokyo Electric Power Co Inc:The | Yield estimation device and computer program |
| CN102323987A (en) * | 2011-09-09 | 2012-01-18 | 北京农业信息技术研究中心 | Crop leaf area index assimilation method |
| CN102651096A (en) * | 2012-04-28 | 2012-08-29 | 中国农业大学 | Method for estimating yield of winter wheat by assimilating characteristics of leaf area index time-sequence curve |
| CN103345707A (en) * | 2013-06-04 | 2013-10-09 | 中国科学院遥感与数字地球研究所 | Crop maturation stage remote sensing prediction method based on multi-source remote sensing data |
Non-Patent Citations (3)
| Title |
|---|
| Estimating regional winter wheat yield by assimilation of time series of HJ-1 CCD NDVI into WOFOST-ACRM model with Ensemble Kalman Filter;Hongyuan Ma etal.;《Mathematical and Computer Modeling》;20130831;第58卷;第759-770页 * |
| 基于遥感信息与作物模型集合卡尔曼滤波同化地区域冬小麦产量预测;黄健熙 等;《农业工程学报》;20120229;第28卷(第4期);第142-148页 * |
| 条件植被温度指数的四维变分与集合卡尔曼同化方法;王维 等;《农业工程学报》;20111231;第27卷(第12期);第184-190页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104134095A (en) | 2014-11-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104134095B (en) | Crop yield estimation method based on scale transformation and data assimilation | |
| CN110751094B (en) | Crop yield estimation method based on GEE comprehensive remote sensing image and deep learning method | |
| Huang et al. | Assimilating a synthetic Kalman filter leaf area index series into the WOFOST model to improve regional winter wheat yield estimation | |
| CN102651096B (en) | Method for estimating yield of winter wheat by assimilating characteristics of leaf area index time-sequence curve | |
| Huang et al. | The improved winter wheat yield estimation by assimilating GLASS LAI into a crop growth model with the proposed Bayesian posterior-based ensemble Kalman filter | |
| Ma et al. | Estimating regional winter wheat yield by assimilation of time series of HJ-1 CCD NDVI into WOFOST–ACRM model with Ensemble Kalman Filter | |
| Stathopoulos et al. | Wind power prediction based on numerical and statistical models | |
| CN109344865B (en) | Data fusion method for multiple data sources | |
| Hazarika et al. | Estimation of net primary productivity by integrating remote sensing data with an ecosystem model | |
| CN106951979A (en) | The crop maturity phase Forecasting Methodology that remote sensing, crop modeling are merged with weather forecast | |
| CN107423850B (en) | Regional corn maturity prediction method based on time series LAI curve integral area | |
| CN110633841B (en) | Provincial range plot scale data assimilation yield prediction method based on set sampling | |
| CN103345707A (en) | Crop maturation stage remote sensing prediction method based on multi-source remote sensing data | |
| CN101480143A (en) | Method for predicating single yield of crops in irrigated area | |
| CN106483147B (en) | Long-time sequence passive microwave soil moisture precision improvement research method based on multi-source data | |
| CN102346808B (en) | Method for inverting LAI (leaf area index) from HJ-1 satellite data | |
| CN108537679B (en) | Remote sensing and crop model fused region scale crop emergence date estimation method | |
| Zhang et al. | Developing a process-based and remote sensing driven crop yield model for maize (PRYM–Maize) and its validation over the Northeast China Plain | |
| Pereira et al. | Solar irradiance modelling using an offline coupling procedure for the Weather Research and Forecasting (WRF) model | |
| CN115375036A (en) | Crop maturity prediction method based on fusion of remote sensing and light energy utilization rate model and weather | |
| CN116522145B (en) | Drought prediction method considering space-time constraint and vegetation condition | |
| CN117611997A (en) | Interpretable forest overground biomass remote sensing estimation method and device | |
| Zhang et al. | Enhanced Feature Extraction From Assimilated VTCI and LAI With a Particle Filter for Wheat Yield Estimation Using Cross-Wavelet Transform | |
| Xia et al. | A spatial frequency/spectral indicator-driven model for estimating cultivated land quality using the gradient boosting decision tree and genetic algorithm-back propagation neural network | |
| AU2021105536A4 (en) | A High Spatial-Temporal Resolution Method for Near-Surface Air Temperature Reconstruction |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant |