CN104134095A - Crop yield estimation method based on scale transformation and data assimilation - Google Patents

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

A kind of crop yield estimating and measuring method based on yardstick conversion and data assimilation

Technical field

The invention belongs to agricultural remote sensing field, be specifically related to a kind of crop yield estimating and measuring method based on yardstick conversion and data assimilation.

Background technology

Utilizing data assimilation that remote-sensing inversion parameter is fused to crop mechanism process model, is the important channel of current improvement area crops growth simulation precision and raising agricultural output assessment precision.On the yardstick of field, crop growth model based on mechanism processes such as crop photosynthesis, breathing, transpiration, nutrition relies on its inherent physical process and kinetic mechanism, can accurately simulate crop object take " my god " growth and development state and the output of crop on the continuous evolution of time step and single-point yardstick.And when crop modeling is extended to regional scale, due to earth's surface, near surface environment heterogeneity, the uncertainty of space distribution of the starting condition in crop modeling, soil parameters, crop parameter, meteorological forcing factors and the difficulty that data is obtained have been caused.Satellite remote sensing has space continuously and the advantage of time dynamic, can effectively solve this bottleneck of area crops parameter acquiring difficulty.Yet remote sensing earth observation is owing to being subject to the restriction of the factors such as satellite spatial and temporal resolution, can't really disclose intrinsic procedure mechanism that crop growth and output forms, individual growth development condition and with the relation of Meteorological Elements in China, and this advantage place of crop modeling just.Data assimilation, by coupling remote sensing observations and crop modeling, can be realized both mutual supplement with each other's advantages, to improving the estimation of area crops output, has huge potentiality.

Because high time resolution remotely-sensed data can be flutterred crop growth variation and the phenology information of catching, therefore crop growth monitoring and the yield by estimation are had very important significance, yet, often spatial resolution is all very low for current high time resolution remote sensing satellite data in orbit, AVHRR for example, MODIS, MERIS and SPOT vegetation spatial resolution are 250m-1km yardstick.This low spatial resolution remotely-sensed data, has further increased pixel heterogeneous, has caused remote sensing observations yardstick to simulate not mating between yardstick with crop modeling, has greatly limited the precision of data assimilation model.In addition, the MODIS LAI product of standard is towards all vegetation pattern exploitations of Global Scale, is not to design and develop for crops.Research shows, MODIS LAI exists the serious phenomenon of underestimating in winter wheat region.Therefore, directly assimilating standard MODIS LAI can cause LAI and output after assimilation unpractical on the low side.The method that yardstick falls in the method for correcting based on geo-statistic (variation function) yardstick and pixel decomposition is to process two kinds of means of pixel heterogeneity and scale effect.Yet, a common weak point of these two kinds of methods is resolution datas between the requirement high and medium of a series of crop growth phase, and is a very complicated problem to heterogeneous pixel modeling, requires strict approximate and some prioris, this is difficult for obtaining in actual applications.Potential method is based on ground surface sample LAI, middle high-resolution remotely-sensed data and a low resolution remotely-sensed data, builds secondary scale-transformation method, generates the revised LAI of yardstick.Assimilation yardstick revised LAI, estimation that will larger raising Regional Fall Wheat output.To overcome the error that direct assimilation MODIS LAI product brings, it is key issue urgently to be resolved hurrily in application data assimilation model the yield by estimation practice.

Summary of the invention

For solving the following problem existing in prior art: " how building the scale-transformation method in data assimilation model; the revised LAI of assimilation time sequence scale is to crop growth model; do not mate to reduce yardstick between remote sensing observations and crop modeling the data assimilation model error causing ", the invention provides and a kind ofly change and the crop yield estimating and measuring method of data assimilation based on yardstick.

The invention provides a kind of crop yield estimating and measuring method based on yardstick conversion and data assimilation, concrete steps are as follows:

Meteorologic parameter in S1 collection research district, crop parameter, soil parameters and management parameters the single-point yardstick that carries out WOFOST crop modeling with it are demarcated; Based on meteorologic parameter, carry out compartmentalization afterwards, complete the spatialization of WOFOST crop modeling;

The remote sensing image of S2 based on multidate, the agrotype sample point that combined ground investigation obtains, in conjunction with the phenology feature of crops to be measured, builds classifying rules, obtains agrotype distribution plan and crops purity chart of percentage comparison to be measured;

S3 carries out filtering to crops to be measured time series breeding time MODIS LAI curve, the impact of polluting to eliminate shortage of data and cloud; Actual measurement LAI based on ground actual measurement sample point and multidate vegetation index (TM VI) build and return statistical model, obtain the TM LAI in 30 meters, region;

S4 merges ground actual measurement sample point LAI, multidate TM LAI and filtered MODIS LAI information, builds secondary yardstick transformation model, and rise time sequence scale is adjusted LAI;

S5, based on the ground sample point LAI in crucial phenological period, calculates the standard deviation of remote sensing observations and crop modeling simulation, obtains crops to be measured remote sensing observations error and the crop modeling error in crucial phenological period;

The rescaling LAI that S6 obtains with S4 and crop modeling simulation LAI, and introduce remote sensing observations error and crop modeling error, build four-dimensional variation cost function, uses optimization algorithm to minimize cost function, by iteration repeatedly, obtains and optimize crop modeling parameter;

The optimization crop modeling parameter substitution crop modeling that S7 obtains S6, the pixel of selecting crops purity to be measured to be greater than 60%, by pixel unit operation crop modeling simulation output, is then aggregated into administrative unit, territory, county, output territory, county crop yield prediction to be measured, instructs grain-production.

Wherein, described in S1, meteorologic parameter is the parameters such as the highest temperature, the lowest temperature, total radiation, vapour pressure, wind speed, precipitation.

Wherein, described in S1, crop parameter is phenology feature of crops distribution parameter, crops etc.

Wherein, described in S1, based on meteorologic parameter, carry out compartmentalization, the spatialization that completes crop modeling refers to carries out the coupling of locus to the parameter of collecting in conjunction with remote sensing image, to meteorologic parameter and the required accumulated temperature parameter of crop modeling, adopt Kriging interpolation algorithm, complete parameter regionization and demarcate.

Wherein, the Landsat TM that described in S2, the remote sensing image of multidate is multidate.

Wherein, described in S2, obtain agrotype distribution plan and crops purity chart of percentage comparison to be measured for adopting the Decision Tree Algorithm of C5.0, phenology feature and spectral signature in conjunction with crops, build classifying rules, obtain agrotype distribution plan, adopt 1km graticule mesh fit, generate the crops purity number percent spatial distribution map to be measured of 1km, in order to reduce the impact of non-crops to be measured region on data assimilation result, the pixel of selecting crops purity to be measured to be greater than 60%, as the pixel of follow-up County Scale output estimation.

Wherein, carry out filtering described in S3, adopt coenvelope line (Savitzky-Golay, S-G) filtering algorithm, its expression formula is shown in formula (1):

$y_J^{*}=\frac{\underset{i=-m}{\overset{i=m}{\mathrm{\Σ}}}{C}_i{Y}_{j+i}}{N}---\left(1\right)$

In formula: Y _j+irepresent the value in a window on original LAI curve, the radius that m is window, N is convolution number, the width of window is 2m+1; Y _j ^*represent filtered LAI value; C represents the filter factor of i LAI value.

Wherein, described in S3, carry out filtering, on to the improved basis of S-G filtering, form coenvelope line filtering method, by LAI timing variations trend, LAI value is divided into true value and falsity, by the mode of local circulation iteration, S-G filter value is replaced to falsity, the LAI smooth curve synthetic new with true value, progressively matching forms the coenvelope line of LAI time series data.

Wherein, described in S4, building secondary yardstick transformation model consists of 2 committed steps:

1) from 30 meters of crops TM LAI to be measured of S3 formation zone, select the TM LAI of 3 key developmental stages;

2) utilize the ratio coefficient between the S-G MODIS LAI in the key developmental stages TM LAI selected and corresponding period, adjust other key developmental stages S-G MODIS LAI, in crops to be measured ascent stage and the decline stage of breeding time, utilize respectively Logistic curve to carry out matching

Logistic curvilinear equation is shown in formula (2):

$y\left(t\right)=\frac{c}{1+e^{a+\mathrm{bt}}}+d---\left(2\right)$

Wherein, t is MODIS LAI seasonal effect in time series index, and y (t) is LAI value corresponding to t time, and a and b are fitting parameters, and c+d is maximum LAI value, and d is LAI initial value, i.e. first value in LAI time series.

Wherein, described in S5, the remote sensing observations error in crucial phenological period of crops to be measured is standard variance calculating based between actual measurement LAI and yardstick correction LAI; Described in S5, crop modeling source of error is in 2 model parameters optimizing, TDWI and SPAN.

Wherein, use described in S6 optimization algorithm to minimize cost function and adopt SCE_UA algorithm to minimize four-dimensional variation cost function, four-dimensional variation cost function press formula (3) calculating:

$J\left(x\right)={\mathrm{\Σ}}_{k=1}^2{(x_k-x_{k0})}^T{B}^{-1}(x_k-x_{k0})+{\mathrm{\Σ}}_{i=1}^N{(y_i-H_i\left(x\right))}^T{Q}_0^{-1}(y_i-H_i\left(x\right))---\left(3\right)$

Wherein, k represents the number of Optimized model parameter in cost function, x _kthe WOFOST parameter of representing optimized is at the numerical value of interval range, x _k0the empirical value of the WOFOST number of parameters of representing optimized, B represents model error, with the error of 2 Optimized model parameters, represents, the transposition of T representing matrix, N represents the number of times of time series remote sensing observations data, y _irepresent remote sensing observations LAI, H _i(x) represent WOFOST modeling LAI, Q ₀represent the error of remote sensing LAI.

Wherein, S7 concrete steps are: select crops purity to be measured to be greater than 60% pixel by continuous iteration, reinitialize 2 crop modeling parameters, LAI and the output of crop modeling output are constantly changed, then in cost function, contrast under the minimum difference of MODIS LAI and crop modeling output LAI, when following three conditions of convergence meet one, can finish assimilation.Obtain the numerical value of Optimal Parameters:

1. the rear parameter value to be optimized of continuous 5 circulations has been retracted to the codomain scope of appointment;

2. target function value cannot improve 0.0001% after 5 circulations;

3. the number of times of calculation cost function is over 10000 times;

Model parameter substitution crop modeling operation by after optimizing, obtains output, and the crops distribution number percent to be measured in conjunction with pixel, gathers the per unit area yield result on output area by territory, county Administrative boundaries.

Wherein, described crops to be measured are preferably winter wheat.

The present invention also provides the application of the described crop yield estimating and measuring method based on yardstick conversion and data assimilation in instructing crops production.

Compared with prior art, beneficial effect is in the present invention:

The present invention has overcome the unmatched problem of yardstick between remote sensing observations pixel and crop modeling analogue unit, has improved the precision of data assimilation model, is suitable for the winter wheat yields estimation of County Scale.

Embodiment

Below in conjunction with embodiment, the specific embodiment of the present invention is described in further detail.Following examples are used for illustrating the present invention, but are not used for limiting the scope of the invention.

Embodiment 1

The winter wheat yields estimation of Baoding and Hengshui Prefecture of take is that example is further set forth technical scheme of the present invention.The flow process of the Regional Fall Wheat output estimating and measuring method based on yardstick conversion and data assimilation of the present embodiment comprises:

Meteorologic parameter in S1 collection research district, crop parameter, soil parameters and management parameters the single-point yardstick that carries out WOFOST crop modeling with it are demarcated; Based on meteorologic parameter, carry out compartmentalization afterwards, complete the spatialization of crop modeling.Specific as follows:

Selecting Baoding and Hengshui Prefecture is survey region, and this district is important winter wheat producing region, Hebei province.Obtain following data: choose 6 meteorological site the highest/lowest temperature, built-up radiation, vapour pressure, wind speed, 6 required meteorological elements of dewatering model; The soil parameters that in study area, agricultural weather testing station gathers and crop parameter etc.; The Fen Xian of Hebei province winter wheat yields statistics in 2009; MODIS15A2 data product, the coordinate projection of unified all data is Albers Conical Equal Area, georeferencing is WGS1984.

For insensitive crop parameter in WOFOST model, but or the less parameter of the higher variation range of sensitivity, by By consulting literatures or directly adopt the default value of WOFOST model to determine (as grown the long DLO of the suitableeest light, the minimum temperature of emerging lower limit TBASE, leaf maintain breathing with RML etc.); The crop parameter high for sensitivity, variation range is large, can be by consulting the pertinent literature of this study area or utilizing measured data to calculate (as the accumulated temperature TSUM1, TSUM2, leaf assimilation quotient conversion ratio CVL, specific leaf area SLATB etc.) obtaining.Also the span that definite parameters such as priori that can study area are possible and numerical value (as the lifetime SPAN of blade 35 ℃ time, maximum CO2 Assimilation rate AMAXTB, the efficiency of light energy utilization EFFTB of single blade assimilation CO2 etc.).

Utilize the observation data of 6 meteorological site in study area, in conjunction with gram in golden spatial interpolation methods to 2008～2009 years, the highest temperature, the lowest temperature, total radiation, vapour pressure, wind speed, six key elements of quantity of precipitation were carried out interpolation day by day, and be converted into the standard format of WOFOST model, obtain the corresponding weather data of each lattice point of study area.For the parameter of easily obtaining in spatial dimension, such as the accumulated temperature in sowing time to seeding stage, the accumulated temperature in seeding stage to florescence, the accumulated temperature in florescence to maturity stage etc., utilize the daily mean temperature of meteorological observation station acquisition and the biological zero point of winter wheat to subtract each other the accumulated temperature that obtains every day, accumulated temperature every day by the cumulative different phenology stage obtains above three accumulated temperature variablees, same adopt golden method of interpolation in identical with weather data interpolation gram, obtain the accumulated temperature parameter of spatial variations in whole study area.

S2 utilizes U.S.'s thematic mapper (Thematic Mapper of study area multidate, TM) image, the agrotype sample point that combined ground investigation obtains, adopt the Decision Tree Algorithm of C5.0, phenology feature in conjunction with crops, build classifying rules, obtain agrotype distribution plan, with 1km graticule mesh fit, generate the winter wheat purity number percent spatial distribution map of 1km, in order to reduce the impact of non-winter wheat region on data assimilation result, the pixel of selecting winter wheat purity to be greater than 60%, as the basis of follow-up County Scale output estimation.

S3 to test block the MODIS LAI data in whole breeding time synthetic by time series, to each grid cell rise time sequence curve.The MRT projection crossover tool specifically providing by NASA, by the data of Hebei Province inlay, projection conversion and format conversion, choose 1 to 177 day (Julian date), phase image superposes totally 23 time, generates the time-serial position based on grid cell.Time series MODIS LAI curve carries out filtering, to eliminate cloud, pollutes the shortage of data causing.The concrete discontinuous problem of coenvelope line (Savitzky-Golay, S-G) filtering solution leaf area index change procedure of using, its Principle representation is as follows:

$y_J^{*}=\frac{\underset{i=-m}{\overset{i=m}{\mathrm{\Σ}}}{C}_i{Y}_{j+i}}{N}---\left(1\right)$

In formula: Y _j+irepresent the value in a window on original LAI curve, the radius that m is window, N is convolution number, the width of window is 2m+1; Y _j ^*represent filtered LAI value; C represents the filter factor of i LAI value.

On to the improved basis of S-G filtering, form coenvelope line filtering method, by LAI timing variations trend, LAI value is divided into true value and falsity, by the mode of local circulation iteration, S-G filter value is replaced to falsity, the LAI smooth curve synthetic new with true value, progressively matching forms the coenvelope line of LAI time series data.Concrete steps are as follows:

(1) the LAI value point that has cloud to pollute is carried out to linear interpolation.Because cloud covers or the impact of atmospheric effect, can make some pixel values with larger noise, to these pixels, adopt the credible pixel value of its periphery to carry out space linear interpolation, generate initial LAI data.

(2) with S-G filtering, carry out the extraction of LAI timing variations trend.According to hypothesis, LAI time series data should be followed the progressive feature of vegetation phenology dynamic change, and in LAI sequential, sudden change is considered to cause due to noise that cloud or atmosphere produce.Therefore, each winter wheat pixel application S-G filtering in study area is carried out smoothly, obtained LAI seasonal effect in time series variation tendency data LAItr.

(3) determine the weights W i of each value point in LAI time series, for follow-up reconstruction result provides basis for estimation.The computing formula of Wi is:

$W_i=\left\{\begin{array}{cc}1& {\mathrm{LAI}}_i\&GreaterEqual;{\mathrm{LAI}}_i^{\mathrm{tr}}\\ 1-d_i/d_{\max }& {\mathrm{LAI}}_i<{\mathrm{LAI}}_i^{\mathrm{tr}}\end{array}\right.$

Wherein, $d_i=|{\mathrm{LAI}}_i-{\mathrm{LAI}}_i^{\mathrm{tr}|},$ D _maxmaximal value for di.

(4) generate new LAI time series.From above theory, less than the LAItr value on time trend curve containing noisy LAI value, these points are carried out to simple replacement and generate new LAI1 timing curve.

${\mathrm{LAI}}_i^1=\left\{\begin{array}{cc}{\mathrm{LAI}}_i& {\mathrm{LAI}}_i\&GreaterEqual;{\mathrm{LAI}}_i^{\mathrm{tr}}\\ {\mathrm{LAI}}_i^{\mathrm{tr}}& {\mathrm{LAI}}_i<{\mathrm{LAI}}_i^{\mathrm{tr}}\end{array}\right.$

(5) LAI1 is carried out to S-G filtering, then the LAI value after filtered data and space interpolation is compared.Equally, if the LAI value after space interpolation is less than newly-generated LAI j, with LAI j, replace, generate new sequence LAI j+1, iteration is carried out filtering reconstruct, and filter effect is by F _ksize evaluate:

$F_k=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\left(\right|{\mathrm{LAI}}_i^{k+1}-{\mathrm{LAI}}_i|\&times;W_i)$

LAIik is the LAI value of i pixel after the k time filtering reconstruct, and Wi is the weight that step (3) calculates, and Fk is that the k time filtered iteration stops judgement parameter.The condition that exits iteration is F _k-1>=F _k≤ F _k+1.

Obtained the cloudless TM data of test block 3 scapes, in time, is respectively mutually: 2009-3-14,2009-5-17,2009-6-10, build respectively recurrence statistical model with corresponding actual measurement LAI.

On March 14th, 2009: due to the middle ten days in March winter wheat in the stage of turning green, exposed soil is larger on the impact of vegetation index, the statistical model of therefore selecting SAVI vegetation index to set up between LAI and TM VI is:

$\mathrm{LAI}=\frac{\ln [(1-\mathrm{SAVI}/1.2581)/0.9130]}{-0.8377}$

Coefficient of determination R ²=0.849.

On May 17th, 2009, the statistical model of setting up between actual measurement LAI and TM NDVI is:

$\mathrm{LAI}=\frac{\ln [(1-\mathrm{NDVI}/1.0866)/3.3790]}{-0.3994},$

Coefficient of determination R ²=0.742.

On June 10th, 2009, the statistical model of setting up between actual measurement LAI and TM NDVI is:

$\mathrm{LAI}=\frac{\ln [(1-\mathrm{NDVI}/9.7639)/1.0081]}{-0.0155}$

Coefficient of determination R ²=0.874.

3 recurrence statistical models based on above-mentioned, the spatial distribution map of the 30 meters of winter wheat in acquisition study area.

S4 merges ground actual measurement sample point LAI, multidate TM LAI and filtered MODISLAI information, builds secondary yardstick transformation model, and rise time sequence scale is adjusted LAI.

Secondary yardstick transformation model comprises 2 steps, first, by the actual measurement LAI in 3 phenological periods that generate in S3 step and TM VI, returns statistical model, 30 meters of formation zone winter wheat TMLAI; Second step, utilize the ratio coefficient between the S-G MODIS LAI in existing key developmental stages TM LAI and corresponding period, adjust other key developmental stages S-G MODIS LAI, ascent stage (turning green to the jointing stage) at During Growing Period of Winter Wheat, obtain the adjustment LAI(boot stage in 4 phenological periods, jointing stage, stand up the phase, period of seedling establishment); In the decline stage (jointing is to the maturity stage) of winter wheat, obtain the adjustment LAI(boot stage in 4 phenological periods, heading stage, florescence, maturity stage), utilize respectively Logistic curve to carry out matching.

Logistic curvilinear equation is as follows:

$y\left(t\right)=\frac{c}{1+e^{a+\mathrm{bt}}}+d$

Wherein, t is MODIS LAI seasonal effect in time series index, and y (t) is LAI value corresponding to t time, and a and b are fitting parameters, and c+d is maximum LAI value, and d is LAI initial value, i.e. first value in LAI time series.

S5, based on the ground sample point LAI in crucial phenological period, calculates the standard deviation of remote sensing observations and crop modeling simulation, obtains winter wheat remote sensing observations error and the crop modeling error in crucial phenological period.

Model error thinks to derive from 2 parameters (TDWI and SPAN) of optimization, and calculating by analysis and obtaining the observational error of yardstick correction LAI 7 phenological periods is 0.34; Model error is set to 0.52.

The rescaling LAI that S6 obtains with S4 and crop modeling simulation LAI, and introduce remote sensing observations error and crop modeling error, build four-dimensional variation cost function, uses optimization algorithm to minimize cost function, by iteration repeatedly, obtains and optimize crop modeling parameter.

At S1, carry out on the basis of crop modeling demarcation, in Experimental Area, operation WOFOST crop modeling simulation LAI, in conjunction with the yardstick correction LAI generating in S4, introduces remote sensing observations and modeling error in S5, builds four-dimensional variation cost function, as follows:

$J\left(x\right)={\mathrm{\&Sigma;}}_{k=1}^2{(x_k-x_{k0})}^T{B}^{-1}(x_k-x_{k0})+{\mathrm{\&Sigma;}}_{i=1}^N{(y_i-H_i\left(x\right))}^T{Q}_0^{-1}(y_i-H_i\left(x\right))$

Wherein, k represents the number of Optimized model parameter in cost function, x _kthe WOFOST parameter of representing optimized is at the numerical value of interval range, x _k0the empirical value of the WOFOST number of parameters of representing optimized, B represents model error, with the error of 2 Optimized model parameters, represents, the transposition of T representing matrix, N represents the number of times of time series remote sensing observations data, y _irepresent remote sensing observations LAI, H _i(x) represent WOFOST modeling LAI, Q ₀represent the error of remote sensing LAI.

The pixel of selecting winter wheat purity to be greater than 60%, use optimized algorithm, constantly iteration reinitializes 2 WOFOST model parameters, LAI and the output of model output are changed, then in cost function, contrast under the minimum difference of MODIS LAI and WOFOST model output LAI, obtain the most optimized parameter.SCE_UA algorithm is used for searching for globally optimal solution in initial parameter space, makes cost function Fast Convergent.When following three conditions of convergence meet one, can finish assimilation, obtain the numerical value of Optimal Parameters.

(1) the rear parameter value to be optimized of continuous 5 circulations has been retracted to the codomain scope of appointment;

(2) target function value cannot improve 0.0001% after 5 circulations;

(3) number of times of calculation cost function is over 10000 times.

The optimization crop modeling parameter substitution crop modeling that S7 obtains S6, the pixel of selecting winter wheat purity to be greater than 60%, by pixel unit operation crop modeling simulation output, is then aggregated into administrative unit, territory, county, output territory, county winter wheat yield.TDWI and SPAN parameter substitution WOFOST model running by after optimizing, obtain output.In conjunction with the winter wheat number percent of pixel, by territory, county Administrative boundaries, gather the winter wheat yield result on output area.

The present invention is by utilizing S-G filtering algorithm to carry out filtering to During Growing Period of Winter Wheat time series MODISLAI curve, eliminated the impact that shortage of data and cloud pollute, generated level and smooth LAI time-serial position, for secondary yardstick transformation model provides high-precision data source.Then by secondary yardstick transformation model, generated time series rescaling LAI, effective integration MODIS LAI seasonal effect in time series tendency information and TM LAI information, improved the precision of LAI, for assimilation process provides high-precision data source.Finally, the LAI of the LAI after adjusting by assimilation and crop modeling simulation, has improved the precision of winter wheat yields estimation on region.

The producer of grain and consumer need to understand timely and accurately grain yield information, according to this method, can obtain yield data in large area in the winter wheat maturity stage, for national departments concerned, carrying out science decisions such as the judgement of grain feelings, grain regulation and control etc. provides important scientific basis, and can be used as the important evidence of grain trade.

Method of the present invention also can be for the estimation of other crop yields.

Although above the present invention is described in detail with a general description of the specific embodiments, on basis of the present invention, can make some modifications or improvements it, this will be apparent to those skilled in the art.Therefore, these modifications or improvements, all belong to the scope of protection of present invention without departing from theon the basis of the spirit of the present invention.

Claims (10)

1. the crop yield estimating and measuring method based on yardstick conversion and data assimilation, is characterized in that, concrete steps are as follows:

Meteorologic parameter in S1 collection research district, crop parameter, soil parameters and management parameters the single-point yardstick that carries out WOFOST crop modeling with it are demarcated; Based on meteorologic parameter, carry out compartmentalization afterwards, complete the spatialization of WOFOST crop modeling;

The remote sensing image of S2 based on multidate, the agrotype sample point that combined ground investigation obtains, in conjunction with the phenology feature of crops to be measured, builds classifying rules, obtains agrotype distribution plan and crops purity chart of percentage comparison to be measured;

S3 carries out filtering to crops to be measured time series breeding time MODIS LAI curve, the impact of polluting to eliminate shortage of data and cloud; Actual measurement LAI based on ground actual measurement sample point and multidate vegetation index TM VI build and return statistical model, obtain the TM LAI in 30 meters, region;

S4 merges ground actual measurement sample point LAI, multidate TM LAI and filtered MODISLAI information, builds secondary yardstick transformation model, and rise time sequence scale is adjusted LAI;

S5, based on the ground sample point LAI in crucial phenological period, calculates the standard deviation of remote sensing observations and crop modeling simulation, obtains crops to be measured remote sensing observations error and the crop modeling error in crucial phenological period;

The rescaling LAI that S6 obtains with S4 and crop modeling simulation LAI, and introduce remote sensing observations error and crop modeling error, build four-dimensional variation cost function, uses optimization algorithm to minimize cost function, by iteration repeatedly, obtains and optimize crop modeling parameter;

The optimization crop modeling parameter substitution crop modeling that S7 obtains S6, the pixel of selecting crops purity to be measured to be greater than 60%, by pixel unit operation crop modeling simulation output, is then aggregated into administrative unit, territory, county, output territory, county crop yield prediction to be measured, instructs grain-production.

2. the crop yield estimating and measuring method based on yardstick conversion and data assimilation described in claim 1, it is characterized in that, described in S1, based on meteorologic parameter, carry out compartmentalization, the spatialization that completes crop modeling refers to carries out the coupling of locus to the parameter of collecting in conjunction with remote sensing image, to meteorologic parameter and the required accumulated temperature parameter of crop modeling, adopt Kriging interpolation algorithm, complete parameter regionization and demarcate.

3. the crop yield estimating and measuring method based on yardstick conversion and data assimilation described in claim 1, is characterized in that, the Landsat TM that the remote sensing image of multidate is multidate described in S2.

4. the crop yield estimating and measuring method based on yardstick conversion and data assimilation described in claim 1, it is characterized in that, described in S2, obtain agrotype distribution plan and crops purity chart of percentage comparison to be measured for adopting the Decision Tree Algorithm of C5.0, phenology feature and spectral signature in conjunction with crops, build classifying rules, obtain agrotype distribution plan, adopt 1km graticule mesh fit, generate the crops purity number percent spatial distribution map to be measured of 1km, in order to reduce the impact of non-crops to be measured region on data assimilation result, the pixel of selecting crops purity to be measured to be greater than 60%, pixel as follow-up County Scale output estimation.

5. the crop yield estimating and measuring method based on yardstick conversion and data assimilation described in claim 1, is characterized in that, carries out filtering described in S3, employing coenvelope line (Savitzky-Golay, S-G) filtering algorithm, and its expression formula is shown in formula (1):

$y_J^{*}=\frac{\underset{i=-m}{\overset{i=m}{\mathrm{\&Sigma;}}}{C}_i{Y}_{j+i}}{N}---\left(1\right)$

In formula: Y _j+irepresent the value in a window on original LAI curve, the radius that m is window, N is convolution number, the width of window is 2m+1; Y _j ^*represent filtered LAI value; C represents the filter factor of i LAI value.

6. crop yield estimation side 0 method based on yardstick conversion and data assimilation described in claim 1, is characterized in that, builds secondary yardstick transformation model described in S4 to consist of 2 committed steps:

1) from 30 meters of crops TM LAI to be measured of S3 formation zone, select the TM LAI of 3 key developmental stages;

2) utilize the ratio coefficient between the S-G MODIS LAI in the key developmental stages TM LAI selected and corresponding period, adjust other key developmental stages S-G MODIS LAI, in crops to be measured ascent stage and the decline stage of breeding time, utilize respectively Logistic curve to carry out matching, Logistic curvilinear equation is shown in formula (2):

$y\left(t\right)=\frac{c}{1+e^{a+\mathrm{bt}}}+d---\left(2\right)$

Wherein, t is MODIS LAI seasonal effect in time series index, and y (t) is LAI value corresponding to t time, and a and b are fitting parameters, and c+d is maximum LAI value, and d is LAI initial value, i.e. first value in LAI time series.

7. the crop yield estimating and measuring method based on yardstick conversion and data assimilation described in claim 1, is characterized in that, the remote sensing observations error in crucial phenological period of crops to be measured is standard variance calculating based between actual measurement LAI and yardstick correction LAI described in S5; Described in S5, crop modeling source of error is in 2 model parameters optimizing, TDWI and SPAN.

8. the crop yield estimating and measuring method based on yardstick conversion and data assimilation described in claim 1, it is characterized in that, described in S6, use optimization algorithm to minimize cost function and adopt SCE_UA algorithm to minimize four-dimensional variation cost function, four-dimensional variation cost function is pressed formula (3) and is calculated:

$J\left(x\right)={\mathrm{\&Sigma;}}_{k=1}^2{(x_k-x_{k0})}^T{B}^{-1}(x_k-x_{k0})+{\mathrm{\&Sigma;}}_{i=1}^N{(y_i-H_i\left(x\right))}^T{Q}_0^{-1}(y_i-H_i\left(x\right))---\left(3\right)$

Wherein, k represents the number of Optimized model parameter in cost function, x _kthe WOFOST parameter of representing optimized is at the numerical value of interval range, x _k0the empirical value of the WOFOST number of parameters of representing optimized, B represents model error, with the error of 2 Optimized model parameters, represents, the transposition of T representing matrix, N represents the number of times of time series remote sensing observations data, y _irepresent remote sensing observations LAI, H _i(x) represent WOFOST modeling LAI, Q ₀represent the error of remote sensing LAI.

9. the crop yield estimating and measuring method based on yardstick conversion and data assimilation described in claim 1, it is characterized in that, S7 concrete steps are: select crops purity to be measured to be greater than 60% pixel by continuous iteration, reinitialize 2 crop modeling parameters, LAI and the output of crop modeling output are constantly changed, then in cost function, contrast under the minimum difference of MODIS LAI and crop modeling output LAI, when following three conditions of convergence meet one, can finish assimilation.Obtain the numerical value of Optimal Parameters:

1. the rear parameter value to be optimized of continuous 5 circulations has been retracted to the codomain scope of appointment;

2. target function value cannot improve 0.0001% after 5 circulations;

3. the number of times of calculation cost function is over 10000 times;

Model parameter substitution crop modeling operation by after optimizing, obtains output, and the crops distribution number percent to be measured in conjunction with pixel, gathers the per unit area yield result on output area by territory, county Administrative boundaries.

10. the application of the crop yield estimating and measuring method based on yardstick conversion and data assimilation in instructing crops production described in claim 1-9 any one.