CN105321120A - Assimilation evapotranspiration and LAI (leaf area index) region soil moisture monitoring method - Google Patents

Abstract

The invention discloses an assimilation evapotranspiration and LAI (leaf area index) region soil moisture monitoring method, and the method comprises the specific steps: S1, calibrating a crop growth model, and completing the determination of spatial parameters; S2, generating an MODIS ET and LAI time sequence curve through superposition; S3, constructing an LAI time sequence curve through filtering; S4, calculating the one-order differential monotonicity of the curve; S5, operating an SWAP model grid by grid, outputting the ET and LAI, and calculating the one-order differential monotonicity; S6, building a cost function according to the difference between the one-order differential monotonicity obtained at step S4 and step S5, taking an seedling emergence date and irrigation time as to-be-optimized parameters, and obtaining an optimal value through employing an optimization algorithm; S7, inputting the optimal parameter into a model, and simulating soil moisture grid by grid. The method estimates the soil moisture through a data assimilation method, introduces two types of remote sensing data, improves the precision of large-area monitoring of soil moisture, and is suitable for the estimation of the soil moisture of regional-scale of farmland.

Description

The regional soil moisture monitoring method of assimilation evapotranspiration and leaf area index

Technical field

The present invention relates to agricultural remote sensing technical field, particularly a kind of regional soil moisture monitoring method of assimilating evapotranspiration and leaf area index.

Background technology

Soil moisture is the chief component of surface water cycle, directly affects the energy cycle of matter process on earth's surface, is the material base of crop growth and growth, especially modern agricultural production and the important indicator parameter in studying.Therefore, the Soil Water evaluation method developing large regions has important practical significance for precision agriculture, yield forecast and agriculture feelings monitoring field.

Traditional soil moisture evaluation method mainly contains gravimetric method, tensiometer method, time domain reflectometry etc., but these methods can only gather the discrete data of sampling point, need to drop into a large amount of manpower and materials, the low and representative difference of efficiency, the soil moisture be difficult on regional scale obtains continuously.The remote sensing inversion method of soil moisture has the advantage of quick obtaining regional scale data, compensate for the deficiency of classic method, mainly contains thermal inertia method, temperature-vegetation index, active microwave remote sensing method and passive microwave remote sensing method etc.But the restriction that satellite remote sensing platform transit time is short, cannot observe, be subject to weather condition as required, thus can only obtain limited amount, time-discrete valid data.

MODIS data are owing to having higher temporal resolution and spectral resolution, in crop growth monitoring and water flux estimation etc., there is advantage, it is the most frequently used data source of Monitoring of Soil Moisture By Remote Sensing, but adaptive different model is needed for Cover treatment degree condition, computation process is very complicated, and the soil moisture information of shallow surface can only be obtained, cause its application to be also very restricted.

Mechanism model can provide simulated object continuous evolution process over time and space, SWAP soil-water-air-plant integration model can simulate the soil water movement in farmland under Different Irrigation level, and exportable relevant Farmland Water and ET component, but require to describe and irrigate data, this is difficult to accurately obtain in extensive area simulation.Data assimilation method can bonding mechanism model on point, remote sensing observations advantage on the whole, the estimation precision of crop soil moisture is improved.

Evapotranspiration comprises transpiration and soil evaporation, is the Important Parameters of surface water capacity balance and heat balance, is also the important indicator of soil moisture, directly reflects surface water thermal equilibrium state.The ET of assimilation regional scale, having very important value for estimation farmland water consumption with carrying out irrigating studying, is one of emphasis approach of the estimation improving crop soil moisture.Assimilation method more common at present minimizes difference between farmland ET that the high spatial resolution ET data of SEBAL model inversion and modeling obtain by genetic algorithm, finally obtain full seasonal effect in time series ET numerical value, but high spatial resolution data not easily obtain owing to being subject to the impact of cloud, and its temporal resolution can not meet assimilation requirement.

The ET product of MODIS is the large regions ET data product that uniquely can obtain at present.But in the winter wheat main producing region of North China of China, because MODIS data spatial resolution is lower, mixed pixel effect is serious, causes product to there is system error on the low side, directly assimilation MODISET and LAI product will cause worse assimilation effect.And a large amount of field inspection experiments finds, the variation tendency of During Growing Period of Winter Wheat time series MODISET and LAI product is very identical with actual measurement variation tendency.Therefore, assimilation MODISET and LAI time series trend change information, to mechanism model, to overcome the space scale mismatch problem between remote sensing observations parameter and modeling, improves regional scale soil moisture estimation precision.

Summary of the invention

(1) technical matters that will solve

The technical problem to be solved in the present invention is: how to do accurate description to the irrigation parameter in SWAP model, the problem that the numerical value simultaneously avoiding MODISET and LAI data to cause due to mixed pixel effect is on the low side.

(2) technical scheme

For solving the problems of the technologies described above, the invention provides a kind of regional soil moisture monitoring method of assimilating evapotranspiration and leaf area index, comprising the following steps:

S1: calibration experiment district crop growth model, geometric match is carried out to remotely-sensed data and ground parameter, use geo-statistic Spatial Interpolation Method to realize the solid formatted meteorology of SWAP model, soil, crop and controling parameters each grid cell in region, complete the spatialization of crop modeling;

S2: to MODISET and the LAI product temporally its superimposition in test block whole breeding time, to each grid cell rise time sequence curve;

S3: based on filtering algorithm, the time series MODISLAI curve obtained in reconstruct S2, to eliminate the impact of shortage of data and cloud pollution;

S4: MODISLAI timing curve filtered in the MODISET timing curve in S2 and S3 is carried out first order difference calculating and extracts monotonicity information;

S5: carry out on the basis of crop modeling demarcation at S1, runs SWAP crop modeling by pixel cell, carries out first order difference to simulating MODISET and the MODISLAI time series obtained and extracts monotonicity information;

S6: MODISET and the MODISLAI curve first order difference monotonicity similarities and differences obtained separately by S4 and S5 set up cost function, emerge date and irrigation time of selection is Optimal Parameters, SCE-UA optimized algorithm is adopted to make cost function Fast Convergent, final when the condition of convergence is satisfied, obtain optimized parameter;

S7: by emerge date and irrigation time substitution crop modeling after optimization, obtain soil moisture by pixel cell simulation.

In said method:

In described step S1:

Described geo-statistic Spatial Interpolation Method is anti-distance weighting (IDW) method of interpolation;

The crop modeling adopted is SWAP model, gather soil, meteorology and the crop parameter in study area, the parameter of remote sensing image and collection is carried out to the coupling of locus, for the direct use experience value of insensitive model parameter, only responsive to soil moisture, use each pixel of geo-statistic anti-distance weighting method of interpolation to remote sensing image to give parameter value to the insensitive model parameter of LAI and ET timing curve feature, forming region parameter, completes the spatialization of crop modeling;

Described meteorological element is daily maximum temperature, the lowest temperature, built-up radiation at sunshine, vapour pressure, wind speed and precipitation;

MODISET and LAI product in described step S2 is MOD16A2 and MCD15A3.

Filtering algorithm in described step S3 adopts coenvelope line filtering algorithm, and filtering is by following 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 one piece of window on original LAI curve, m is the radius of window, and N is convolution number, and the width of window is 2m+1, represent the LAI value at filtering rear hatch center, C _irepresent the filter factor of i-th LAI value.

The first order difference monotonicity adopted in described step S4 or S5 is calculated as follows formula (2):

In formula, Dif _irepresent the curve first order difference monotonicity at timing node i place, whether representative is increasing function herein, if 1 growth, if-1 decline, if 0 represents not change; Jnd represents that the just noticeable difference on curve is other; Subscript Obs is expressed as the MODIS product curve of observation, and subscript Sim then represents the numerical curve of SWAP modeling; ET is the value at timing node i place in ET time series, carries out after the same method for LAI.

The cost function will set up in described step S6 calculates by following formula (3):

In formula, J is the target function value of cost function, and i is the observation date of MODIS data.

In described step S6, variable to be optimized is the irrigation time of emerging date and 2 periods in SWAP model, assimilation ET and LAI seasonal effect in time series first order difference monotonicity obtains the optimum irrigation time of emerging date and 2 periods, cost function comprises the first order difference monotonicity similarities and differences of ET and LAI time series shape as formula (3) simultaneously, obtains optimized parameter during total cost function convergence.

The condition of convergence in described step S6 has one of following condition:

(1) after continuous 5 circulations, model parameter value to be optimized has been retracted to the codomain scope of specifying;

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

(3) number of times of calculation cost function is more than 1000 times,

The model parameter obtained according to assimilation is run by SWAP model, exports soil moisture result.

Method provided by the invention is for monitoring the change of soil moisture in process of crop growth.

(3) beneficial effect

The regional soil moisture monitoring method of assimilation evapotranspiration provided by the present invention and leaf area index, by setting up ET and the LAI time-serial position of winter wheat in breeding time, adopt the strategy different from assimilation method in the past, extract the first order difference monotonicity information of its timing curve to build cost function, thus efficiently avoid numerical value that MODIS product data cause due to scale effect problem on the low side, and using emerge date and irrigation time as optimized variable, solve in SWAP model the problem of irrigating data and being difficult to obtain, for the estimation of large area region yardstick soil moisture provides effective way.

Compared with CN201210133136.0 (publication number is CN102651096A), the present invention makes a big difference, and is specially:

1) test crop modeling, the present invention does not need to carry out parameter calibration;

2) data that the present invention uses also use ET product data except MODISLAI data;

3) the present invention does not need to carry out curve fitting;

4) model of the present invention is SWAP model, and the model of CN102651096A is WOFOST crop modeling;

5) what the present invention assimilated is not curvilinear characteristic point but curve monotonicity information;

6) what the present invention exported is the soil moisture content of designated depth (10cm);

7) obtain the Error weight of remote sensing observations error at three some places by formula (6), the present invention of this step does not need to carry out.

Sum up, the present invention is compared with the document, and the model difference adopted, assimilation strategy is different, and Optimal Parameters is different, and final object is also different, and the model in CN102651096A can only be used for assessing, and cannot export soil moisture.

Accompanying drawing explanation

Fig. 1: the regional soil moisture monitoring method flow diagram of assimilation evapotranspiration provided by the invention and leaf area index;

Fig. 2-1 and Fig. 2-2:ET and LAI timing curve first order difference monotonicity information extraction effect schematic diagram, wherein Fig. 2-1 is the schematic diagram of ET product, and Fig. 2-2 is the schematic diagram of LAI product;

Fig. 3-1,3-2: actual measurement soil moisture (weight method) and monitoring effect of the present invention are verified, wherein Fig. 3-1 is SWAP direct modeling (assimilation), and Fig. 3-2 is assimilation MODISET and LAI tendency information;

Fig. 4: Guanzhong，Shanxi plain soils moisture monitoring schematic diagram, upper figure are the effect of SWAP simulation before assimilation, and figure below is the effect after assimilation.

Embodiment

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

Embodiment 1:

To estimate that technical scheme of the present invention is set forth in Guanzhong，Shanxi Plain further.The flow process in the Guanzhong，Shanxi Plain of the present embodiment as shown in Figure 1, comprising:

Step S1, carries out the demarcation of test block crop modeling, carries out space geometry registration, and complete crop modeling parameter region to remotely-sensed data and ground parameter.

The central Shaanxi plain selecting Shaanxi Province is test block, this area comprises Fengxiang, Qishan, Fufeng, Mei County, Yang Ling, acrobatic skill Deng31Ge county, be under the jurisdiction of Xi'an, Xianyang, Baoji, Weinan and five, Tongchuan prefecture-level city respectively, coverage is 106 ° of 18 ' E--110 ° 30 ' E, 33 ° of 40 ' N--35 ° of 28 ' N, as shown in Fig. 2-1.This ground weather belongs to temperate zone monsoon sex climate, year samming 6 ~ 13 DEG C.Annual precipitation 500 ~ 800 millimeters, wherein accounts for 60%, mostly be short-time storm 6 ~ September, and Winter-Spring precipitation is less.Owing to having a moderate climate, make a clear distinction between the four seasons, hydrothermal condition is suitable for the plantation of crops.Obtain following data: according to study area external envelope scope, the day choosing 33 National Meteorological stations is the highest/lowest temperature, built-up radiation at sunshine, vapour pressure, wind speed, 6 meteorological elements needed for dewatering model; The soil parameters that in study area, agricultural weather testing station gathers and crop parameter; The controling parameters such as longitude and latitude, elevation; MOD16A2 and MCD15A3 data product, by the data uniform coordinate containing geographical location information, completes spatial match.

Research on utilization district is with the U.S. thematic mapper (ThematicMapper of phase, TM) image interpretation obtains winter wheat planting area, generate the data mask of 1 kilometer, winter wheat plantation number percent in computation grid unit, sets 20% threshold value and rejects the too low unit of winter wheat planting proportion.To output and the weather data unified metric of model, weather data and crop parameter interpolation are generated the raster data of 1 kilometer of every pixel.

According to " soil-water-air-plant integration model " (Soil-Water-Atmosphere-Plant of Wageningen University's exploitation, be called for short SWAP model), at the preliminary calibration result of Guanzhong，Shanxi plains region winter wheat, for insensitive crop parameter, or the higher but parameter that parameter variation range is less of susceptibility, the default value of direct employing SWAP model is determined, and parameter that parameter variation range larger higher for susceptibility, by By consulting literatures or utilize measured value to calculate, also can according to the scope of document or priori determination parameter, then determined by the optimizer (FSEOPT) that model is subsidiary, the localization completing model is demarcated.The parameter that soil moisture and ET/LAI timing curve feature are all responsive is carried out to the assimilation process of subsequent step.

Step S2, to the temporally its superimposition synthesis of MODISET and the LAI data in test block whole breeding time, to each grid cell rise time sequence curve.Especially by the MRT projection transform instrument that NASA provides, the data of Hebei Province are carried out inlay, projection transform and format conversion, in choose 1 to 177 days (Julian date), totally 23 scape images superpose, and generate the time-serial position of each grid cell.

Step S3, carries out filtering to the time series MODISLAI curve obtained in step S2, pollutes to eliminate cloud the shortage of data caused.Specifically utilize coenvelope line (Savitzky-Golay, SG) filtering to solve the discontinuous problem of leaf area index change procedure, 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 (1): Y _j+irepresent the value in one piece of window on original LAI curve, m is the radius of window, and N is convolution number, and the width of window is 2m+1; represent filtered LAI value; C represents the filter factor of i-th LAI value.This polynomial design reduces exceptional value to retain high numerical value.Concrete operations are as follows:

(1) SG filtering is carried out to initial LAI time series, obtain level and smooth after result, and preserve respectively level and smooth before and level and smooth after sequence.

(2) contrast two sequences of previous step, generate new sequence with following formula, as initiation sequence,

$N_i^t=\left\{\begin{array}{cc}{O}_i& \mathrm{if}{O}_i^{t-1}\&GreaterEqual;N_i^{t-1}\\ N_i^{t-1}& \mathrm{if}{O}_i^{t-1}<N_i^{t-1}\end{array}\right.---\left(2\right)$

In formula (2), O and N is respectively initial and filtered LAI value, and subscript t represents iterations, and subscript i is LAI time series index.

(3) repetition (1) (2) two step, until whole sequence Σ N _i ^tbe less than the threshold value 0.1 of specifying.

Step S4, the MODISET time-serial position in the MODISLAI time-serial position obtained by filtering algorithm in step S3 and step S2 carries out the calculating (see accompanying drawing 2) of first order difference monotonicity,

Its equation expression is as follows:

In formula, Dif _irepresent the curve first order difference monotonicity at timing node i place, whether representative is increasing function herein, if 1 growth, if-1 decline, if 0 represents not change; Jnd represents that the just noticeable difference on curve is other; Subscript Obs is expressed as the MODIS product curve of observation, and subscript Sim then represents the numerical curve of SWAP modeling; ET is the value at timing node i place in ET time series, carries out after the same method for LAI.

Step S5, carries out on the basis of crop modeling demarcation at S1, running SWAP crop modeling, and S5 is similar by grids of pixels unit, to the calculating (see accompanying drawing 2) of simulating ET and the LAI time series obtained and carry out equally first order difference monotonicity.

Step S6, and the total cost function of ET curve and LAI curve is minimized, the cost function that set up calculates by following formula (4):

In formula, J is the target function value of cost function, and i is the observation date of MODIS data.

SCE-UA algorithm is used for searching for globally optimal solution in initial parameter space simultaneously, make it constantly reinitialize to emerge date and irrigation time parameter, ET and the LAI sequential Long-term change trend characteristic sum soil moisture that SWAP model is exported constantly changes, finally make cost function Fast Convergent, assimilation can be terminated when following three conditions of convergence meet one:

(1) after continuous 5 circulations, parameter value to be optimized has been retracted to the codomain scope of specifying;

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

(3) number of times of calculation cost function is more than 1000 times.

Step S7, by emerge date and irrigation time parameter substitution SWAP model after optimization, obtains soil moisture by pixel cell simulation.

The regional soil moisture monitoring method of the assimilation evapotranspiration described in the embodiment of the present invention and leaf area index, merge the advantage of remotely-sensed data and crop modeling, the time series ET variation characteristic of remote sensing observations is extracted as assimilation variable, to be undertaken emerging by optimized algorithm the parameter estimation of date and irrigation time, achieve the assimilation of remote sensing and model, solve MODISET and the LAI problem that can not directly apply in data assimilation on the low side.By assimilating the soil moisture data that ET and LAI time series trend change information obtains simultaneously, according to the central Shaanxi plain area Heyang, Pucheng, tramplle on the measured data of shop, three former, six eyeballs in Fufeng and Lantian, can obtain surveying the relation between soil water content and modeling soil moisture content.Wherein, measured data comprises the dark soil moisture content of underground 10cm of six eyeballs, four periods (on November 11st, 2011, on April 3rd, 2012, on May 5th, 2012 and on June 6th, 2012).

Fig. 3 is the statistical relationship between soil moisture content and actual measurement soil moisture content that before and after assimilation, model exports, as can be seen from the figure after assimilation, the precision of the soil moisture content that model exports, higher than the precision of not assimilating the soil moisture content that situation drag exports, illustrates the validity of assimilation method in soil moisture content simulation.In non-assimilation situation, the coefficient of determination R2=0.661 between the soil moisture content that SWAP modeling obtains and actual measurement soil moisture content, root-mean-square error RMSE=0.041.After the assimilation of first order difference cost function method, R2 raising is reduced to 0.034 from 0.661 to 0.809, RMSE from 0.041.

Fig. 4 is the soil moisture spatial distribution map of regional scale, and result shows, the soil moisture Spatial Difference after assimilation increases, the more realistic rule of space distribution.

Above embodiment is only for illustration of the present invention; and be not limitation of the present invention; the those of ordinary skill of relevant technical field; without departing from the spirit and scope of the present invention; can also make a variety of changes and modification; therefore all equivalent technical schemes also belong to category of the present invention, and scope of patent protection of the present invention should be defined by the claims.

Claims (9)

1. assimilate a regional soil moisture monitoring method for evapotranspiration and leaf area index, it is characterized in that, comprise the following steps:

S1: calibration experiment district crop growth model, geometric match is carried out to remotely-sensed data and ground parameter, use geo-statistic Spatial Interpolation Method to realize the solid formatted meteorology of SWAP model, soil, crop and controling parameters each grid cell in region, complete the spatialization of crop modeling;

S2: to MODISET and the MODISLAI product temporally its superimposition in test block whole breeding time, to each grid cell rise time sequence curve;

S3: based on filtering algorithm, the time series MODISLAI curve obtained in reconstruct S2, to eliminate the impact of shortage of data and cloud pollution;

S4: MODISLAI timing curve filtered in the MODISET timing curve in S2 and S3 is carried out first order difference and calculates monotonicity;

S5: carry out on the basis of crop modeling demarcation at S1, runs SWAP crop modeling by pixel cell, carries out first order difference to simulating MODISET and the MODISLAI time series obtained and calculates monotonicity;

S6: MODISET and the MODISLAI curve first order difference monotonicity similarities and differences obtained separately by S4 and S5 set up cost function, emerge date and irrigation time of selection is parameter to be optimized, SCE-UA optimized algorithm is adopted to make cost function Fast Convergent, final when the condition of convergence is satisfied, obtain the optimal value of parameter;

S7: by emerge date and irrigation time substitution crop modeling after optimization, obtain soil moisture by pixel cell simulation.

2. monitoring method according to claim 1, it is characterized in that, the crop modeling adopted in described step S1 is SWAP model, gather the soil in study area, meteorology and crop parameter, by same place, space geometry registration is carried out to the model parameter of remote sensing image and collection, for the direct use experience value of insensitive model parameter, for only responsive to soil moisture, the each pixel of anti-distance weighting method of interpolation to remote sensing image is used to give parameter value to the insensitive model parameter of LAI and ET timing curve feature, complete parameter regionization to demarcate, the parameter of compartmentalization is substituted into the spatialization that model completes crop modeling,

Described meteorological element is daily maximum temperature, the lowest temperature, built-up radiation at sunshine, vapour pressure, wind speed and precipitation.

3. monitoring method according to claim 1, is characterized in that, MODISET and the LAI product data in described step S2 are respectively MOD16A2 and MCD15A3.

4. monitoring method according to claim 1, is characterized in that, adopts the filtering of coenvelope line in described step S3, and filtering is by following 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 one piece of window on original LAI curve, m is the radius of window, and N is convolution number, and the width of window is 2m+1, represent the LAI value at filtering rear hatch center, C _irepresent the filter factor of i-th LAI value.

5. monitoring method according to claim 1, is characterized in that, the first order difference monotonicity adopted in described step S4 or S5 is calculated as follows formula (2):

In formula, Dif _irepresent the curve first order difference monotonicity at timing node i place, whether representative is increasing function herein, if 1 growth, if-1 decline, if 0 represents not change; Jnd represents the minimum differentiation threshold value on curve; Subscript Obs is expressed as the MODIS product curve of observation, and subscript Sim then represents the numerical curve of SWAP modeling; ET is the value at timing node i place in ET time series, builds after the same method for LAI.

6. monitoring method according to claim 1, is characterized in that, the cost function will set up in described step S6 calculates by following formula (3):

In formula, J is the target function value of cost function, and i is the observation date of MODIS data.

7. monitoring method according to claim 1, is characterized in that, in described step S6, variable to be optimized is emerge date and the irrigation time in SWAP model; With ET and LAI for assimilation link variable, shape, as the first order difference monotonicity similarities and differences of formula (3), obtains optimized parameter during total cost function convergence.

8. monitoring method according to claim 1, is characterized in that, the condition of convergence in described step S6 has one of following condition:

(1) after continuous 5 circulations, model parameter value to be optimized has been retracted to the codomain scope of specifying;

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

(3) number of times of calculation cost function is more than 1000 times,

The model parameter obtained according to assimilation is run by SWAP model, exports soil moisture result.

9. the application in the change of the soil moisture in monitoring crop growth course of the monitoring method described in any one of claim 1-8.