CN105760978A - Agricultural drought grade monitoring method based on temperature vegetation drought index (TVDI) - Google Patents

Abstract

The invention discloses an agricultural drought grade monitoring method based on a temperature vegetation drought index (TVDI). The method comprises the following steps of: 1, data preparation; 2, land surface temperature (LST) data reconstruction; 3, construction of crop plantation area normalized difference vegetation index-land surface temperature (NDVI-LST) feature space; 4, TVDI calculation; and 5, drought grade monitoring based on the TVDI. The invention brings forward an LST data reconstruction method based on multi-year background values and area fluctuation values, a cultivated land area multi-year NDVI-LST feature space is constructed for crops, the TVDI is calculated, a crop drought grade monitoring model based on a supervision classification idea is designed for drought grade remote sensing monitoring, the model can quite accurately reflect drought threat degrees of the crops under difference conditions in real time, and the method has great significance in monitoring, early warning and prevention of agricultural disasters.

Description

A kind of agricultural drought disaster grade monitoring method based on temperature vegetation drought index (TVDI)

Technical field

The present invention relates to a kind of with MODIS remotely-sensed data for key data, by reconstructing land surface temperature data and setting up crops temperature arid vegetation index in conjunction with normalized differential vegetation index, agricultural drought disaster grade monitoring model based on crops temperature vegetation drought index carries out agricultural arid grade monitoring method, similarity evaluating model realizes the method for Natural Disaster rapid evaluation, is specially a kind of agricultural drought disaster grade monitoring method based on temperature vegetation drought index (TVDI).

Background technology

Drought is one of topmost natural disaster in human being's production life, and compared with other natural disasters, drought has occurrence frequency height, persistent period length, involves the feature that scope is wide.Arid is caused by water deficit, and under hazard-affected body existent condition, large range of long-time arid just can cause drought.Drought is the process of a complexity, accumulation, its performance is mainly that surface water is dry, level of ground water declines, surface vegetation growth is affected, its impact mainly has crop production reduction or total crop failure, human livestock drinking water difficulty, and the commercial production that continuous high temperature hydropenia causes is obstructed.China is population and large agricultural country, and Food Security is the major issue concerning social stability and economic development.Harm agricultural production brought at China's drought is the most serious.Since 1949, arid the grain disaster area caused, every year on average up to 3.2 hundred million mu, rate of wherein causing disaster nearly 30%, bring loss not only to China's agricultural, also seriously constrain social and economic development.

The features such as drought rank is many, complicated mechanism, has and involves scope greatly, and the persistent period is long, and space and time above variability is strong, do not have ripe effective method at present on agricultural drought disaster is monitored；Traditional meteorological drought monitoring: meteorological site is sparse, and skewness, and data space precision is not high, it is impossible to effectively damage caused by a drought and agricultural disaster are set up corresponding relation；Agricultural drought disaster the statistics of geological disaster situation: although the Disaster degree of agricultural arid can be obtained more accurately, but the statistical data acquisition cycle is longer, expends a large amount of manpower and materials, and data are many in units of administrative division, spatial accuracy is limited, additionally, there is also the problems such as Statistical Criteria disunity.For the problems referred to above, remote sensing technology has the features such as multi-platform, multisensor, high-resolution, it is possible to obtains damage caused by a drought information quickly, accurately, can being repeated property of devastated be observed at short notice；Spectral characteristic of ground according to different-waveband, remote sensing technology can obtain the information such as surface water change, vegetation growth, soil moisture content directly or indirectly.Compare meteorological drought monitoring and agricultural drought disaster the statistics of geological disaster situation, agricultural drought disaster remote sensing monitoring can obtain soil moisture content transformation and crop growthing state in high precision, on a large scale, rapidly, the drought stress degree that under different condition, crop is subject to can be reflected in real time, exactly, significant in the monitoring of agricultural drought disaster, early warning and strick precaution.

Summary of the invention

Remote sensing technology can obtain earth's surface information quickly, repeatedly, on a large scale, significant in Monitoring of drought.Wherein utilize temperature vegetation drought index (TVDI) of optics and IRMSS thermal band, it is effectively combined the indicator normalized differential vegetation index (NDVI) of vegetation state, important parameter land surface temperature (LST) thermally equilibrated with surface water, is widely used in remote sensing drought study on monitoring.But TVDI is in actual Agriculture Drought is monitored, and there are the following problems: (1) missing values is more makes LST data discontinuous on space-time, it is impossible to extract TVDI further；(2) in NDVI-LST feature space builds, there is the impact of bare place pixel, and traditional TVDI lacks comparability between many annual datas.For as above two problems, the present invention proposes the LST data reconstruction method based on background value for many years Yu region undulating value, compare existing reconstructing method, the space missing value of large area in restructural list scape image, with the shortage of data of single pixel long-term sequence, and the change details of LST can be retained to a certain extent；With crops for object of study, build region, arable land NDVI-LST feature space for many years, and crops temperature vegetation drought index (C-TVDI) are proposed, on this basis, calculate the C-TVDI of monitoring section crops Critical growing period, and design and carry out Henan Province's drought loss remote sensing monitoring based on the crops drought loss monitoring model of supervised classification thought, in order to fight calamities and provide relief, decision-making provides reference.

The present invention adopts the following technical scheme that a kind of agricultural drought disaster grade monitoring method based on temperature vegetation drought index (TVDI) for achieving the above object, comprises the steps:

(1) preparation of data:

The present invention utilizes the data that remote sensing technology obtains to carry out agricultural drought disaster study on monitoring, wherein remotely-sensed data used is NASA (http://ladsweb.nascom.nasa.gov/) MODISLST product, vegetation index product, Land_use change covering product and the dem data provided, and arable land multiple crop index data；Statistical data is the Model on Sown Areas of Farm data that provide of the cultivated area data in China Statistical Yearbook, disaster area data and plant husbandry management department of the Ministry of Agriculture (http://www.zzys.moa.gov.cn/) and crops phenological calendar.

(2) based on the land surface temperature data reconstruction of background value Yu undulating value:

Land surface temperature data (being called for short LST below) are reconstructed according to corresponding method based on the MODISLST data product obtained.

(3) structure of proportion of crop planting district normalized differential vegetation index-land surface temperature feature space:

Extract draught monitor region to plough, contemporaneous data for many years is utilized jointly to build crops each trophophase proportion of crop planting district normalized differential vegetation index-land surface temperature (being called for short NDVI-LST below) scatterplot, and matching each phase is in wet limit equation, build NDVI-LST feature space.

(4) calculating of crops temperature vegetation drought index:

Based on the result of step 3, utilize Price to propose temperature vegetation drought index in nineteen ninety and calculate proportion of crop planting district, monitored area crops temperature vegetation drought index (being called for short C-TVDI below).

(5) monitor based on the drought loss of crops temperature vegetation drought index:

Based on the design philosophy of Supervised classification device, determine the monitored area parameter based on the drought loss monitoring model of crops temperature vegetation drought index by historical data, carry out the monitoring of drought.

As the further scheme of the present invention, the land surface temperature data reconstruction based on background value Yu undulating value of described step (2) includes reconstructing method and reconstructed operation step, concrete grammar and operating process are as follows:

(1) LST reconstructing method

The present invention is based on the thought utilizing initial fields value in Cressman objective analysis method with correct value and jointly approach observation, and by many yearly mean levels of pixel LST, namely background value regards as the initial fields value of pixel；By the undulating value of pixel LST, namely in the periphery radius of influence, interpolation pixel LST, as correcting value, is once corrected interpolation with this by the observation of pixel and the difference of background value.Pixel a (i, j) specific algorithm of interpolation can be expressed as follows:

LST_insert=LST_background+LST_varianceFormula 1-1

${\mathrm{LST}}_{\mathrm{background}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{LST}}_{\mathrm{Gaussian}}$ Formula 1-2

${\mathrm{LST}}_{\mathrm{variance}}=\frac{\mathrm{\Σ}{W}_K\mathrm{\Δ}{\mathrm{LST}}_K}{\mathrm{\Σ}{W}_K}$ Formula 1-3

In formula, LST_insertFor a point LST interpolation result, LST_backgroundFor a point background value for many years.

Owing to LST data existing missing values and low quality data, these points should be removed when calculating background value, it is therefore desirable to LST time series data is reconstructed.The present invention selects required setup parameter less asymmetric Gaussian function fitting method that LST time series data is reconstructed, and to the time series data LST after matching_GaussianAsking for background value for many years of each phase, n is that in data set, a certain issue is according to contained year number, as shown in above formula 1-2；If some pixel missing values or low quality data in time series data is too much, cannot be carried out matching, this pixel value is then substituted by the meansigma methods of ground mulching type similar in the periphery radius of influence (in historical years identical at least over half ground mulching type) pixel, the problem wherein with periphery pixel, elevation difference being existed for disappearance pixel, the present invention utilizes height above sea level often raise 1000m temperature to decline the relation of about 6K, remove the elevation participating in the calculating pixel impact on LST.

LST_varianceFor the undulating value at a point place, by ground mulching type and radius of influence decision, Δ LST is the difference of the observation of high-quality pixel and background value in the radius of influence, and K is the high-quality pixel number of identical earth's surface cover type in the radius of influence, W_KFor its respective weights, following formula calculate and obtain:

$W_{\mathrm{ijK}}=\left\{\begin{array}{cc}\frac{R^2-d_{\mathrm{ijK}}^2}{R^2+d_{\mathrm{ijK}}^2}& (d_{\mathrm{ijK}}<R)\\ 0& (d_{\mathrm{ijK}}\&GreaterEqual;R)\end{array}\right.$ Formula 1-4

In formula, d_ijkDistance for the high-quality pixel of interpolation pixel to identical ground surface type；R is the radius of influence, owing to temperature change within the scope of horizontal range 6000m is generally less than 0.6K^[57], therefore, herein radius of influence R is set to 3, namely the high-quality pixel in 7x7 window participates in interpolation calculation.

(2) LST reconstructed operation step

1) data prediction.The present invention chooses MODIS product for draught monitor district LST data reconstruction.MOD11A2, MOD12Q1 product used in reconstruct first passes through MODIS re-projection instrument (MODISReprojectionTool, MRT) re-projection is WGS84 coordinate system, and extracts the data of LST round the clock in product, round the clock LST quality control file and LAI/fPAR system Land Use/Cover Classification result；Research is also used the digital elevation data (DigitalElevationModel, DEM) provided by NASA, first dem data is carried out registration with MODIS data, then be 1000m by the spatial resolution of dem data from 30m resampling.

Owing to kind and the quantity of institute data are more, before treatment the studies above data are cut out, consider restructing algorithm relates to spatial window convolution algorithm, therefore select rectangular mask file that data are cut, spatial resolution is 1000m, guarantees that monitored area is in rectangular area central authorities as far as possible.After above-mentioned pretreatment, data used is: LST data, round the clock LST quality control file, dem data and MOD12Q1 data round the clock, data above is carried out Band fusion, constitutes day night LST data set, day night LST quality control file data collection, ground mulching categorical data collection and dem data.

2) many years background values of LST are calculated.Asymmetric Gaussian function fitting will be carried out respectively by LST data set round the clock.Afterwards, the time series data collection after inspection matching: due to, in data set, the pixel value of LST is open type temperature, so the time series all values by successful for non-matching pixel place is set to 0；Thereafter, all pixels are calculated this pixel long-time average annual value same period, as this phase background value, it is thus achieved that background value data set, each 46 wave bands round the clock.It is the pixel of 0 to pixel value in background value data set, it is interpolated in 7x7 window around this pixel: first, with reference to ground mulching categorical data collection, choose pixel identical with interpolation pixel ground mulching type in 7x7 window (exceed with interpolation pixel type same number in historical years half namely think identical earth's surface cover type), if not, choose whole pixel；Subsequently, the weighing computation method in Cressman objective analysis method is utilized to calculate each pixel weight according to choosing the pixel distance with interpolation pixel；Afterwards, utilize dem data, calculate the depth displacement choosing pixel with interpolation pixel, utilize height above sea level often to raise 1000m, the relation of temperature decline 6k, interpolation pixel place elevation temperature is arrived in all LST unifications choosing pixel；Finally, the LST after choosing pixel " levelling " in window is weighted summation with its weight, obtains the many years background value data sets of LST of space and time continuous.

3) LST quality control file data collection round the clock is utilized, respectively LST data round the clock are carried out unit's screening item by item: according to LST quality control file round the clock, only retain high-quality pixel value in each wave band, all the other pixel values are set to 0, it is thus achieved that interpolation LST data set round the clock.

4) calculating gained LST background value data set, interpolation LST data set and ground mulching categorical data collection is utilized to be interpolated, concretely comprise the following steps a. and utilize ground mulching categorical data collection, choose in the pixel periphery 7x7 window that pixel value is 0, it is worth the identical pixel of pixel ground mulching type with 0, if not, chooses whole pixel；B. the weighing computation method in Cressman objective analysis method is utilized to calculate the weight choosing pixel；C. calculate and choose the actual value (i.e. high-quality pixel observation) of pixel and the difference of background value, as undulating value；D. the corresponding weight of the undulating value choosing pixel in window is weighted summation, obtains final LST interpolation result.

As the further scheme of the present invention, the structure of described step (3) proportion of crop planting district normalized differential vegetation index-land surface temperature feature space, arable land, draught monitor region is extracted including utilizing, contemporaneous data for many years is utilized jointly to build crops each trophophase proportion of crop planting district normalized differential vegetation index-land surface temperature (being called for short NDVI-LST below) scatterplot, and matching each phase dry and wet limit equation, build NDVI-LST feature space.Specific as follows:

(1) invention utilizes 16 days sintetics MOD13A2 of MODIS vegetation index herein, and the LST8 day data and monitored area crops multiple crop index data on daytime after the inventive method reconstructs, each data spatial resolution is 1000m.

1) MODIS vegetation index product MOD13A2 MODIS re-projection instrument is carried out projection transform, and extract NDVI data therein and quality assessment data (QualityAssurance, QA)；

2) monitoring section data year by year are carried out arable land image element extraction, wherein due in the multiple crop index data of arable land, select the pixel value that pixel area coincide the most with cultivated area data in statistical yearbook, each year, Area distortion controlled below 10%, plough using the pixel of this multiple crop index pixel value as monitored area pixel, is ploughed in monitored area and extract；

3) with reference to 16 grades of criteria for classifications of NDVI mass in QA data, as shown in table 1, the screening NDVI quality of data is medium above pixel, namely at the value of QA data 2～5 pixel less than 4, all the other pixels in NDVI data is assigned to 0；

4) utilizing Timesat3.1.1 to the NDVI data after QA data screening, carry out time series reconstruct, reconstructing method selects Savitzky-Golay filter method, and filter window is set to 5；

5) every for LST8 day data on daytime 2 scapes are merged into 16 days LST data of 1 scape, in merging, the high-quality pixel of original LST data substitutes interpolation pixel, when being all high-quality pixel or interpolation pixel when two pixels, merge algorithm with reference to MOD11A2 user's manual, to two pixel value averaged, as merging data pixel value.

(2) first, many annual datas are grouped on schedule, NDVI, LST data are divided into some groups, each group comprises the equal number same period NDVI, LST data, all crops pixels of each group of data are built NDVI-LST feature space；Then, one of percentage of difference of, minima maximum with the NDVI value in each group of data is for step-length, and progressively in the long group obtained corresponding to different NDVI value, LST is maximum, minima；Finally utilize the dry and wet limit equation that NDVI and LST is maximum, minima matching is respectively organized in data.

As the further scheme of the present invention, the calculating of described step (4) crops temperature vegetation drought index, calculate monitored area C-TVDI value including proposing temperature vegetation drought index based on Price in nineteen ninety.The computational methods of concrete temperature vegetation drought index are expressed as follows:

$\mathrm{TVDI}=\frac{T_s-T_{s\_\min }}{T_{s\_\max }-T_{s\_\min }}$ Formula 1-5

Wherein,

T_{s_max}=a₁+b₁NDVI formula 1-6

T_{s_min}=a₂+b₂NDVI formula 1-7

In formula, T_sFor pixel (i, LST value j), T_{s_max}With T_{s_min}Respectively pixel (i, LST value on the dry limit of NDVI value correspondence j) and wet limit, dry limit, by the scatterplot that the maximum that all kinds of NDVI values in feature space are corresponding is constituted, carries out linear regression and tries to achieve, and wet limit is then tried to achieve by the linear regression of corresponding minima scatterplot；a₁, a₂, b₁, b₂The respectively undetermined coefficient in the equation of dry and wet limit.

As the further scheme of the present invention, described step (5) is monitored based on the drought loss of crops temperature vegetation drought index, including the design philosophy based on Supervised classification device, determine the monitored area parameter based on the drought loss monitoring model of crops temperature vegetation drought index by historical data, carry out the monitoring of drought.Specific as follows:

Based on the design philosophy of Supervised classification device, using former years crop growth period remotely-sensed data with actual measurement of each year drought loss data as training sample, by the study to training sample, it is determined that the parameter in drought loss monitoring model；Recycle the crop growth period remotely-sensed data of up-to-date a year afterwards, this year year-end drought loss is monitored.This model is mainly monitored three parts by disaster area evaluation function, parameter optimization and drought loss and is formed:

(1) disaster area evaluation function: disaster area is defined as the Model on Sown Areas of Farm suffering various natural disaster, and suffer several or natural disaster several times in the same year, wherein to endanger maximum once calculating disaster area, do not repeat meter calamity, the maximum Model on Sown Areas of Farm that affected area by drought was then affected by drought for this year one or many.Consider that disaster area value is for drought maximum effect area to crops in year, if considering the cropping pattern that monitoring section is repeatedly multiple cropping, and crops impact of a drought after harvesting terminates, therefore, disaster area evaluation function is divided into different planting season, is set as follows:

S_i=Max [Mean (s_ik..., s_il) ..., Mean (s_im..., s_in)] formula 1-8

In formula, S_iIt is 1 year disaster area estimated value；K～1 is harvesting of annual first phase crop Critical growing period, and m～n is annual Final Issue harvesting crop Critical growing period, s_ij(j=k～1) is the area suffered from drought of 1 year first phase harvesting each Critical growing period of crop, s_ij(j=m～n) gathers in the area suffered from drought of each trophophase of crop for Final Issue；Owing to needing the accumulation of certain time from drought of suffering from drought, the damage caused by a drought of single issue evidence may not cause drought, and the Spatial Variability of arid is little in a short time, therefore using the meansigma methods of each phase area suffered from drought as disaster area；Annual disaster area value is then the maximum in each phase disaster area.

Area suffered from drought s_ij(j=k～1 or m～n) is by 1 year crops Critical growing period agricultural remote sensing exponent data and this trophophase threshold parameter t_j(j=k～1 or m～n) determines.For crops temperature vegetation drought index, exponential quantity is more arid closer to 1 expression, then the area suffered from drought S of 1 year jth Critical growing period_ijFor this year this issue according in C-TVDI value more than this phase threshold value t_jPixel shared by area, threshold value t_jFor model parameter, parameter and trophophase one_to_one corresponding used.

(2) parameter optimization, namely utilizes optimizing algorithm that the object function set is asked for the process being worth most.Drought loss monitoring model carries out the parameter optimization of model using remotely-sensed data for many years and agricultural disaster area data as training sample.The optimal solution of model parameter is regarded as under this group parameter, calculates and obtains disaster area estimated value and the minimum average B configuration deviation in actual disaster area in years of training sample for many years, is therefore set as by object function:

$\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}|S_i-S_{i0}|$ Formula 1-9

In formula, S_iIt is 1 year disaster area estimated value, S_i0Being 1 year actual disaster area value, n is training sample number, and object function is tried to achieve model parameter t corresponding during minima_j(j=k～1 or m～n) is optimized parameter.

This research selects grid-search algorithms to carry out parameter optimization, and the speed of searching optimization of grid-search algorithms is very fast, it is possible to obtain globally optimal solution, will not be absorbed in locally optimal solution.But wanting setup parameter scope before search, relatively big at parameter area, in the search that step-size in search is less, time consumption for training is longer.

Parameter optimization process is expressed as follows: the remotely-sensed data in training sample is carried out crops Indices extraction by (i), and selects to participate in the crop growth period that model calculates；(ii) minimum target functional value, each trophophase parameter t are set_i0Hunting zone, step-size in search and parameter t_i0Initial value；(iii) the agricultural remote sensing exponent data of extraction is input in the evaluation function of disaster area, it is thus achieved that disaster area estimated value；(iv) disaster area estimated value is brought in object function with disaster area data in training sample, calculating target function value, if this value is less than current minimum target functional value, then retains parameter current, and minima is updated；V () obtains new parameter value by grid-search algorithms, and repeat above-mentioned (iii), (iv) step, if all parameter all travels through, completes optimizing, records optimized parameter.

(3) drought loss monitoring, namely by crop growth year, remotely-sensed data phase to be monitored, is monitored this year drought loss.Sample number owing to participating in training in this research is less, and training sample is limited to the representativeness of the condition of a disaster in various degree, hence with the linear relationship between training sample disaster area estimated value and actual disaster area, the disaster area estimated value in year to be monitored is adjusted.What drought loss was monitored specifically comprises the following steps that (i) calculates actual disaster area and the equation of linear regression of disaster area estimated value under optimized parameter in training sample；(ii) year to be monitored crops Critical growing period Indices is extracted；(iii) disaster area evaluation function is utilized to ask for the disaster area estimated value under optimized parameter；(iv) using estimated value in previous step as linear equation independent variable, it is thus achieved that disaster area estimated value after adjustment；V () utilizes disaster area estimated value after adjusting to calculate disaster-stricken rate (I), and with reference to drought loss criteria table, it is thus achieved that drought loss monitoring result.Drought loss criteria table is as follows:

Table 2 drought loss criteria table

Compared with prior art, the present invention has carried out following improvement: (1) proposes the LST data reconstruction method based on background value for many years Yu region undulating value, compare existing reconstructing method, the space missing value of large area in restructural list scape image, with the shortage of data of single pixel long-term sequence, and the change details of LST can be retained to a certain extent；(2) with crops for object of study, build region, arable land NDVI-LST feature space for many years, and crops temperature vegetation drought index (C-TVDI) are proposed, on this basis, calculate the C-TVDI of monitoring section crops Critical growing period, and design and carry out Henan Province's drought loss remote sensing monitoring based on the crops drought loss monitoring model of supervised classification thought, in order to fight calamities and provide relief, decision-making provides reference.

Accompanying drawing explanation

Fig. 1 is the land surface temperature data reconstruction Technology Roadmap based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI) of the one according to description of the present invention and specific embodiment；

Fig. 2 is the Henan land surface temperature data reconstruction Comparative result figure on daytime 7 day January in 2005 based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI) of the one according to the specific embodiment of the invention；

Fig. 3 is the Henan land surface temperature data reconstruction Comparative result figure at night 1 day January in 2009 based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI) of the one according to the specific embodiment of the invention；

Fig. 4 is the land surface temperature data reconstruction original value reconstruction result scatterplot based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI) of the one according to the specific embodiment of the invention；

Fig. 5 is the NDVI-LST part scatterplot of crops for many years based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI) of the one according to the specific embodiment of the invention；

Fig. 6 is the C-TVDI and agricultural disaster area administrative division profiles versus figure based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI) of the one according to the specific embodiment of the invention；

Fig. 7 is the drought loss monitoring flow chart based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI) of the one according to the specific embodiment of the invention；

Detailed description of the invention

Below in conjunction with the drawings and specific embodiments, the present invention is further elaborated.

Present case, with Henan for study area, mainly uses vegetation index data product, LST data product, land cover pattern/land cover pattern delta data product.All kinds of MODIS product such as following table used herein:

Table 1MODIS land product

Arable land multiple crop index data are mainly used in obtaining the spatial distribution of arable land, Henan pixel, these data are provided by doctor Liu Jianhong, data spatial resolution is 500m, present case is used is Henan area 2001-2011 totally 11 scape data, pixel value is the proportion of crop planting number of times of this pixel annual, bare place pixel is labeled as 0, and the arable land concrete extraction algorithm of multiple crop index data is as follows:

(1) according to average meteorological data for many years, the whole nation is divided a ripe district, Liang Shu district and three ripe districts；

(2) in each processed district, select suitable training sample, determine the shortest Length of growing season of crops, the longest Length of growing season, minimum growth amplitude in each processed district according to training sample；

(3) utilize MODIS data acquisition each pixel enhancement mode meta file (EnhancedVegetationIndex, EVI) time-serial position, and extract thing marquis's parameters such as vegetation growing season number, Length of growing season, growth amplitude；

(4) thing marquis's parameter of image element extraction is compared with the shortest Length of growing season of crops in processed district, pixel place, the longest Length of growing season, minimum growth amplitude, differentiate whether a vegetation growing season belongs to crop growth season, if it is retain, if not then deleting；

(5) the crop growth season number finally given is the multiple crop index of this pixel.

In present case, statistical yearbook data are also used, farming counts evidence one by one and all kinds of Model on Sown Areas of Farm data 3 class statistical data；Wherein farming counts evidence and all kinds of Model on Sown Areas of Farm data one by one, for determining the Critical growing period of staple crops, and screens the remotely-sensed data consistent with crops Critical growing period for agricultural drought disaster study on monitoring；Statistical yearbook data construct remote sensing drought loss monitoring model with model monitoring result is tested.Institute statistical yearbook data include Henan Province's affected area by drought data of each year and cultivated area data, and these data all are from 2001-2011 China Statistical Yearbook；Farming counts evidence one by one and all kinds of Model on Sown Areas of Farm data come from plant husbandry management department of the Ministry of Agriculture (http://www.zzys.moa.gov.cn/), sown area according to all kinds of crops, Crops in Henan Province is winter wheat, summer corn, semilate rice, Semen sojae atricolor, Cotton Gossypii, Brassica campestris L and Semen arachidis hypogaeae, and all kinds of staple crops farming in the research time limit is gone through such as following table:

Present case main flow includes: 1) utilize the LST data reconstruction method of background value and region undulating value for many years, Henan 2000-2011 LST data are reconstructed, 5000 high-quality pixels of random selection, and reconstruct data are carried out precision test；2) utilize the LST data after reconstruct and MODISNDVI data, build Henan area crops NDVI-LST feature space；3) using 2001-2007 crops temperature vegetation drought index data and actual disaster area data as training sample, using 2008-2011 crops temperature vegetation drought index data, actual disaster area data and drought loss data as test samples, crops drought monitoring model is evaluated.

Present case is chosen MODIS product and is studied for Henan area LST data reconstruction.MOD11A2, MOD12Q1 product used in reconstruct first passes through MODIS re-projection instrument (MODISReprojectionTool, MRT) re-projection is WGS84 coordinate system, and extracts the data of LST round the clock in product, round the clock LST quality control file and LAI/fPAR system Land Use/Cover Classification result；Research is also used the digital elevation data (DigitalElevationModel, DEM) provided by NASA, first dem data is carried out registration with MODIS data, then be 1000m by the spatial resolution of dem data from 30m resampling.

Owing to kind and the quantity of institute data are more, before treatment the studies above data are cut out, consider restructing algorithm relates to spatial window convolution algorithm, therefore select rectangular mask file that data are cut, it is sized to 641x622, spatial resolution is 1000m, and Henan, study area is in rectangular area central authorities.After above-mentioned pretreatment, data used is: each 506 scapes of LST data round the clock, round the clock each 506 scapes of LST quality control file, wherein annual 46 scapes, totally 11 years；Dem data 1 scape；MOD12Q1 is totally 11 scape, annual 1 scape, and data above carries out Band fusion, constitutes day night LST data set, day night LST quality control file data collection, ground mulching categorical data collection and dem data, it is simple to the process of restructing algorithm further.

Present case land surface temperature data reconstruction (LST) includes three steps (flow chart is shown in Figure of description 1):

(1) many years background values of LST are calculated.To carry out asymmetric Gaussian function fitting respectively by LST data set round the clock, this step is realized by Timesat3.1.1.Afterwards, the time series data collection after inspection matching: due to, in data set, the pixel value of LST is open type temperature, so the time series all values by successful for non-matching pixel place is set to 0；Thereafter, all pixels are calculated this pixel long-time average annual value same period, as this phase background value, it is thus achieved that background value data set, each 46 wave bands round the clock.It is the pixel of 0 to pixel value in background value data set, it is interpolated in 7x7 window around this pixel: first, with reference to ground mulching categorical data collection, choose pixel identical with interpolation pixel ground mulching type in 7x7 window (namely thought more than 6 times with interpolation pixel type same number in 11 years identical earth's surface cover type), if not, choose whole pixel；Subsequently, the weighing computation method in Cressman objective analysis method is utilized to calculate each pixel weight according to choosing the pixel distance with interpolation pixel；Afterwards, utilize dem data, calculate the depth displacement choosing pixel with interpolation pixel, utilize height above sea level often to raise 1000m, the relation of temperature decline 6k, interpolation pixel place elevation temperature is arrived in all LST unifications choosing pixel；Finally, the LST after choosing pixel " levelling " in window is weighted summation with its weight, obtains the many years background value data sets of LST of space and time continuous.

(2) LST quality control file data collection round the clock is utilized, respectively LST data round the clock are carried out unit's screening item by item: the division of pixel credit rating in saving according to 3.1, only retain high-quality pixel value in each wave band, all the other pixel values are set to 0, obtain interpolation LST data set round the clock, each 506 scapes.

(3) calculating gained LST background value data set, interpolation LST data set and ground mulching categorical data collection is utilized to be interpolated, concretely comprise the following steps 1. and utilize ground mulching categorical data collection, choose in the pixel periphery 7x7 window that pixel value is 0, it is worth the identical pixel of pixel ground mulching type with 0, if not, chooses whole pixel；2. utilize the weighing computation method in Cressman objective analysis method to calculate the weight choosing pixel；3. calculate and choose the actual value (i.e. high-quality pixel observation) of pixel and the difference of background value, as undulating value；4. the corresponding weight of undulating value choosing pixel in pair window is weighted summation, obtains final LST interpolation result (Figure of description 2).

Present case is verify the interpolation precision based on background value Yu undulating value LST data reconstruction algorithm, design accuracy confirmatory experiment, the original data set of LST round the clock has randomly selected 250 scape LST data, the 250 scape data chosen randomly select 5000 high-quality pixels again, pixel value is assigned to 0, generating LST round the clock and verify data set, except checking data set is modified except pixel value with original LST data set, all the other are all identical.Checking data set is carried out based on background value and undulating value LST data reconstruction, and the interpolation result in reconstruction result and original pixel value are contrasted.In comparing result, the maximum deviation of LST interpolation and original pixel value is 15.44K, and average deviation is 0.81K, and the deviation pixel more than 2K accounts for the 5.2% of pixel sum；Visible, the interpolation precision based on background value Yu undulating value LST data reconstruction algorithm is higher, can MODISLST data be repaired preferably.Figure of description 4 is the scatterplot of interpolation result and initial data, and original pixel value and interpolation result have stronger linear dependence, and wherein, linear gradient is 0.986, and intercept is 1.21, R2 is 0.960.(Figure of description 2, Fig. 3, Fig. 4)

Present case utilizes 16 days sintetics MOD13A2 of MODIS vegetation index, and the LST8 day data and Henan area 2001-2011 crops multiple crop index data on daytime after above-mentioned reconstruct, each data spatial resolution is 1000m.Data handling procedure is as follows:

(1) MODIS vegetation index product MOD13A2 MODIS re-projection instrument is carried out projection transform, and extract NDVI data therein and quality assessment data (QualityAssurance, QA)；

(2) study area data year by year are carried out arable land image element extraction, wherein due in the multiple crop index data of arable land, pixel value be 2 pixel area and Henan statistical yearbook in cultivated area data the most identical, each year Area distortion 5%～10% not etc., therefore the pixel being 2 using multiple crop index value is ploughed as study area pixel, is ploughed in study area and extracts；

(3) with reference to 16 grades of criteria for classifications of NDVI mass in QA data, as shown in table 6, the screening NDVI quality of data is medium above pixel, namely at the value of QA data 2～5 pixel less than 4, all the other pixels in NDVI data is assigned to 0；

(4) utilizing Timesat3.1.1 to the NDVI data after QA data screening, carry out time series reconstruct, reconstructing method selects Savitzky-Golay filter method, and filter window is set to 5；

(5) every for LST8 day data on daytime 2 scapes are merged into 16 days LST data of 1 scape, in merging, the high-quality pixel of original LST data substitutes interpolation pixel, when being all high-quality pixel or interpolation pixel when two pixels, merge algorithm with reference to MOD11A2 user's manual, to two pixel value averaged, as merging data pixel value.

Present case utilizes pretreated NDVI, LST data, and proportion of crop planting district, Henan is built NDVI-LST feature space.First, many annual datas are grouped on schedule, NDVI, LST data are divided into 23 groups, each group comprises the 11 scape same period NDVI, LST data, all crops pixels of each group of data are built NDVI-LST feature space；Then, one of percentage of difference of, minima maximum with the NDVI value in each group of data is for step-length, and progressively in the long group obtained corresponding to different NDVI value, LST is maximum, minima；Finally utilizing the dry and wet limit equation that NDVI and LST is maximum, minima matching is respectively organized in data, partial results is such as shown in Figure of description 5.

In result of many phases, crops NDVI-LST scatterplot distribution for many years is broadly divided into three classes: wherein, angular distribution relation on the NDVI-LST scatterplot of the crops for many years distribution coincidence theory in most of period, such as Figure of description 5 (a) o. 11th (late June), 5 (b) the 18th phase (mid-October), 5 (c) the 21st phase (early November), wherein dry limit slope is respectively less than 0, and wet limit be not all the straight line being parallel to NDVI axle, as at Figure of description 5 o. 11th, 18 is interim, wet limit slope value is all higher than 0, at the interim slope of Figure of description 5 the 21st less than 0, from overall distribution, wet limit slope value is all higher than the slope value on dry limit, dry and wet limit increases with NDVI to be drawn close gradually；The NDVI-LST scatter plot distributions of crops for many years of some period is contrary with theoretical Triangle-Profile, as shown in Figure of description 5 (d) the 8th phase (the first tenday period of a month in May), 5 (e) the 13rd phase (late July), wherein dry limit slope is all higher than 0, wet limit slope is respectively less than 0, and dry and wet limit deviates from gradually with the growth of NDVI；The NDVI-LST scatter plot distributions of crops for many years separately having indivedual period is a rectangle, and dry and wet limit is parallel to each other, and slope is all close to 0, such as Figure of description 5 (f) the 3rd phase (mid-February).

In the above results, the different decisions of the difference main crop growth phase corresponding to scatterplot each phase of crops NDVI-LST scatterplot distribution for many years: in the scatterplot of theoretical Triangle-Profile, such as 11 phases, 18 phases and 21 phases, respectively late June time, mid-October and early November, correspond to autumn grain crops sowing/seeding stage, autumn grain crops harvest time and summer grain crops sowing time respectively, in the above three period, proportion of crop planting region NDVI Distribution value is comparatively uniform, bare area, low nurse crop and high nurse crop coexist, and therefore NDVI-LST scatterplot is close with theoretical triangle；And the scatterplot in anti-triangle, such as 8 phases, 13 phases, time is the first tenday period of a month in May and late July respectively, correspond to summer grain crops trophophase and autumn grain crops trophophase respectively, the two Grain Growth Situation in period reaches maximum, most pixels in proportion of crop planting district are high nurse crop, only having minority pixel is bare area or the mixed pixel for bare area with crop, in NDVI-LST scatterplot, point focuses mostly in the middle high level region of NDVI, NDVI low value region only has base point, corresponding LST value does not possess representativeness and change limited difference, therefore dry and wet limit in NDVI low value region closer to；The scatterplot of rectangular distribution, such as 3 phases, time is mid-February, this period is the Wintering Period of summer grain crops winter wheat, now proportion of crop planting region pixel has certain vegetation coverage, scatterplot midpoint to concentrate on the middle low value district of NDVI, but it is very weak to be in Wintering Period Crop transpirstion effect, not possessing the ability regulating canopy surface temperature, therefore dry and wet limit is all in the straight line being parallel to NDVI axle.

Present case utilizes the dry and wet limit equation that said process is tried to achieve to calculate the C-TVDI value in proportion of crop planting region for many years, Henan further, and computational methods are as follows:

$\mathrm{TVDI}=\frac{T_s-T_{s\_\min }}{T_{s\_\max }-T_{s\_\min }}$ Formula 1

Wherein,

T_{s_max}=a₁+b₁NDVI formula 2

T_{s_min}=a₂+b₂NDVI formula 3

In formula, T_sFor pixel (i, LST value j), T_{s_max}With T_{s_min}Respectively pixel (i, LST value on the dry limit of NDVI value correspondence j) and wet limit, dry limit, by the scatterplot that the maximum that all kinds of NDVI values in feature space are corresponding is constituted, carries out linear regression and tries to achieve, and wet limit is then tried to achieve by the linear regression of corresponding minima scatterplot；a₁, a₂, b₁, b₂The respectively undetermined coefficient in the equation of dry and wet limit.

Present case calculates the C-TVDI value part result in proportion of crop planting region for many years, Henan and sees Figure of description 6. with agricultural disaster area administrative division in the same year scattergram

Each city of present case Henan Province disaster area spatial distribution is visible with C-TVDI index spatial distribution comparing result (Figure of description 6), C-TVDI can correctly reflect the spatial distribution of agricultural drought disaster on the whole, but the accuracy on prefecture-level region is not enough, and its reason is likely drought the condition of a disaster and is together decided on by many factors；For further C-TVDI being used for agricultural drought disaster remote sensing monitoring, and by remotely-sensed data and agricultural drought disaster the condition of a disaster opening relationships, present case proposes a kind of drought loss monitoring model based on remotely-sensed data of many phases (C-TVDI), and in conjunction with Henan Province's affected area by drought data, model monitoring result is verified.

Present case drought loss monitoring model is based on the design philosophy of Supervised classification device, using former years crop growth period remotely-sensed data and each year actual measurement drought loss data as training sample, by the study to training sample, it is determined that the parameter in drought loss monitoring model；Recycle the crop growth period remotely-sensed data of up-to-date a year afterwards, this year year-end drought loss is monitored.This model is mainly monitored three parts by disaster area evaluation function, parameter optimization and drought loss and is formed:

(1) disaster area evaluation function: disaster area is defined as the Model on Sown Areas of Farm suffering various natural disaster, and suffer several or natural disaster several times in the same year, wherein to endanger maximum once calculating disaster area, do not repeat meter calamity^[65], maximum Model on Sown Areas of Farm that affected area by drought was then affected by drought for this year one or many.Consider that disaster area value is for drought maximum effect area to crops in year, Henan, study area is mainly the cropping pattern of twice multiple cropping, and crops impact of a drought after harvesting terminates, therefore, disaster area evaluation function is divided into summer grain crops disaster area and autumn grain crops disaster area two parts, is set as follows:

S_i=Max [Mean (s_ik..., s_il), Mean (s_im..., s_in)] formula 4

In formula, S_iIt is 1 year disaster area estimated value；K～1 is summer grain crops crop Critical growing periods, and m～n is autumn grain crop Critical growing period, s_ti(j=k～1) is the area suffered from drought of 1 year each Critical growing period of summer grain crops crop, s_ij(j=m～n) is the area suffered from drought of each trophophase of autumn grain crop；Owing to needing the accumulation of certain time from drought of suffering from drought, the damage caused by a drought of single issue evidence may not cause drought, and the Spatial Variability of arid is little in a short time, therefore using summer/each phase area suffered from drought of autumn grain crops meansigma methods as summer/disaster area of autumn grain crops；Annual disaster area value is then the maximum of summer grain crops and autumn grain crops disaster area.

Area suffered from drought s_ij(j=k～1 or m～n) is by 1 year crops Critical growing period agricultural remote sensing exponent data and this trophophase threshold parameter t_j(j=k～1 or m～n) determines.For crops temperature vegetation drought index, exponential quantity is more arid closer to 1 expression, then the area suffered from drought S of 1 year jth Critical growing period_ijFor this year this issue according in C-TVDI value more than this phase threshold value t_jPixel shared by area, threshold value t_jFor model parameter, parameter and trophophase one_to_one corresponding used.

(2) parameter optimization, namely utilizes optimizing algorithm that the object function set is asked for the process being worth most.Drought loss monitoring model carries out the parameter optimization of model using remotely-sensed data for many years and agricultural disaster area data as training sample.The optimal solution of model parameter is regarded as under this group parameter, calculates and obtains disaster area estimated value and the minimum average B configuration deviation in actual disaster area in years of training sample for many years, is therefore set as by object function:

$\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}|S_i-S_{i0}|$ Formula 5

In formula, S_iIt is 1 year disaster area estimated value, S_i0Being 1 year actual disaster area value, n is training sample number, and object function is tried to achieve model parameter t corresponding during minima_j(j=k～1 or m～n) is optimized parameter.

Conventional optimizing algorithm has genetic algorithm, simulated annealing, particle cluster algorithm and grid-search algorithms.Due to this paper model, to relate to parameter less, and training sample amount is less, therefore selects grid-search algorithms to carry out parameter optimization, and the speed of searching optimization of grid-search algorithms is very fast, it is possible to obtain globally optimal solution, will not be absorbed in locally optimal solution.But wanting setup parameter scope before search, relatively big at parameter area, in the search that step-size in search is less, time consumption for training is longer.

Parameter optimization process is expressed as follows: the remotely-sensed data in training sample is carried out crops Indices extraction by (i), and selects to participate in the crop growth period that model calculates；(ii) minimum target functional value, each trophophase parameter t are set_i0Hunting zone, step-size in search and parameter t_i0Initial value；(iii) the agricultural remote sensing exponent data of extraction is input in the evaluation function of disaster area, it is thus achieved that disaster area estimated value；(iv) disaster area estimated value is brought in object function with disaster area data in training sample, calculating target function value, if this value is less than current minimum target functional value, then retains parameter current, and minima is updated；V () obtains new parameter value by grid-search algorithms, and repeat above-mentioned (iii), (iv) step, if all parameter all travels through, completes optimizing, records optimized parameter.

(3) drought loss monitoring, namely by crop growth year, remotely-sensed data phase to be monitored, is monitored this year drought loss.Sample number owing to participating in training in this research is less, and training sample is limited to the representativeness of the condition of a disaster in various degree, hence with the linear relationship between training sample disaster area estimated value and actual disaster area, the disaster area estimated value in year to be monitored is adjusted.What drought loss was monitored specifically comprises the following steps that (i) calculates actual disaster area and the equation of linear regression of disaster area estimated value under optimized parameter in training sample；(ii) year to be monitored crops Critical growing period Indices is extracted；(iii) disaster area evaluation function is utilized to ask for the disaster area estimated value under optimized parameter；(iv) using estimated value in previous step as linear equation independent variable, it is thus achieved that disaster area estimated value after adjustment；V () utilizes disaster area estimated value after adjusting to calculate disaster-stricken rate (I), and with reference to drought loss criteria table, it is thus achieved that drought loss monitoring result^[64].Drought loss criteria table is as follows:

Table 4 drought loss criteria table

The registration monitoring of present case drought utilizes 2001-2011 Henan Province crops temperature vegetation drought index data and disaster area, China Statistical Yearbook Henan Province, cultivated area data, and drought loss monitoring model is evaluated；According to Henan summer autumn grain crops chief crop, winter wheat and summer corn are in the different demands to moisture of day part, by the jointing of winter wheat, heading and pustulation period, and summer corn take out male and milk stage as the crops Critical growing period in drought loss monitoring model, 6 phase remotely-sensed datas during this trophophase correspondence is annual, each issue corresponding threshold parameter tj (j=1～6) of data.For determining 6 threshold parameters of model, using 2001-2007 annual data as model training sample, carry out parameter optimization and determine the regression equation of the initial estimated value in disaster area and actual disaster area；For the monitoring result of evaluation model, using 2008-2011 annual data as test samples, utilize the optimized parameter of model and equation of linear regression that disaster area of each year is estimated, and drought loss is monitored, and contrastive detection result and actual drought loss.Idiographic flow is shown in Figure of description 7.

Present case utilizes training sample 2001-2007 Henan Province crops temperature vegetation drought index data and disaster area, Henan data, and drought loss monitoring model each main growing period threshold parameter optimizing result is as follows:

The each phase threshold value optimizing result of table 5 model

Utilize optimized parameter to try to achieve training sample disaster area estimated value, and to carry out linear regression, regression equation and the coefficient of determination as follows with actual disaster area value:

S_{Actual measurement}=0.8448 × S_Estimation+ 2693.4 formulas 6

R²=0.9370

According to above-mentioned parameter and regression equation, test samples 2008-2011 crops temperature vegetation drought index data are inputted in drought loss monitoring model respectively, try to achieve test samples each year disaster area estimated value, and obtain each year drought loss monitoring result according to drought loss criteria table, as shown in the table:

Table 6 drought loss monitoring result

In table 6, disaster-stricken rate is the ratio in disaster area and cultivated area, and disaster area deviation ratio is the ratio that disaster area estimated value and actual value absolute deviation account for actual disaster area.Visible by table 6, as the 2008-2011 of test samples, cover from without drought to weight four grades of drought, wherein disaster area deviation ratio increases with the decline of actual drought loss, and as there is weight non-irrigated 2009, actual disaster area is more than 1500khm², and disaster area deviation ratio is only 6%, and at 2010 occurred without drought, actual disaster area was less than 100khm², disaster area deviation is for actual disaster area more than 4 times, from disaster area deviation value it can be seen that the estimated value that the above results is because disaster area of each year all exists 200khm with actual value²The deviation of left and right, deviation ratio then raises along with the minimizing in actual disaster area；From the monitoring result of drought loss, the drought monitoring grade of test samples is all consistent with drought actual grade, has certain feasibility in the qualitative monitoring of disaster loss grade.

Present case describes the agricultural drought disaster grade monitoring technology flow process based on crops temperature vegetation drought index in detail, and for Henan Province, have collected this province's agricultural drought disaster related data of calendar year 2001 to 2011, using 2008-2011 annual data as test samples, disaster area of each year is estimated by the optimized parameter utilizing model with equation of linear regression, and drought loss is monitored, and contrastive detection result and actual drought loss.From the monitoring result of drought loss, the drought monitoring grade of test samples is all consistent with drought actual grade, has certain feasibility in the qualitative monitoring of disaster loss grade, it is possible to meet the basic demand of Droughts grade monitoring.Based on crops temperature vegetation drought index agricultural drought disaster grade monitoring method data can availability strong, method is simple, and processing ease has certain operational use prospect.

Above content is in conjunction with concrete preferred implementation further description made for the present invention, it is impossible to assert that specific embodiment of the invention is confined to these explanations.For general technical staff of the technical field of the invention, without departing from the inventive concept of the premise, it is also possible to make some simple deduction or replace, protection scope of the present invention all should be considered as belonging to.

Claims (5)

1. the agricultural drought disaster grade monitoring method based on temperature vegetation drought index (TVDI), it is characterised in that comprise the steps:

(1) preparation of data:

The present invention utilizes the data that remote sensing technology obtains to carry out agricultural drought disaster study on monitoring, wherein remotely-sensed data used is NASA (http://ladsweb.nascom.nasa.gov/) MODISLST product, vegetation index product, Land_use change covering product and the dem data provided, and arable land multiple crop index data；Statistical data is the Model on Sown Areas of Farm data that provide of the cultivated area data in China Statistical Yearbook, disaster area data and plant husbandry management department of the Ministry of Agriculture (http://www.zzys.moa.gov.cn/) and crops phenological calendar.

(2) based on the land surface temperature data reconstruction of background value Yu undulating value:

Land surface temperature data (being called for short LST below) are reconstructed according to corresponding method based on the MODISLST data product obtained.

(3) structure of proportion of crop planting district normalized differential vegetation index-land surface temperature feature space:

Extract draught monitor region to plough, contemporaneous data for many years is utilized jointly to build crops each trophophase proportion of crop planting district normalized differential vegetation index-land surface temperature (being called for short NDVI-LST below) scatterplot, and matching each phase dry and wet limit equation, build NDVI-LST feature space.

(4) calculating of crops temperature vegetation drought index:

Based on the result of step 3, utilize Price to propose temperature vegetation drought index in nineteen ninety and calculate proportion of crop planting district, monitored area crops temperature vegetation drought index (being called for short C-TVDI below).

(5) monitor based on the drought loss of crops temperature vegetation drought index:

Based on the design philosophy of Supervised classification device, determine the monitored area parameter based on the drought loss monitoring model of crops temperature vegetation drought index by historical data, carry out the monitoring of drought.

2. one according to claim 1 is based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI), it is characterized in that, the land surface temperature data reconstruction based on background value Yu undulating value of described step (2) includes reconstructing method and reconstructed operation step, concrete grammar and operating process are as follows:

(1) LST reconstructing method

The present invention is based on the thought utilizing initial fields value in Cressman objective analysis method with correct value and jointly approach observation, and by many yearly mean levels of pixel LST, namely background value regards as the initial fields value of pixel；By the undulating value of pixel LST, namely in the periphery radius of influence, interpolation pixel LST, as correcting value, is once corrected interpolation with this by the observation of pixel and the difference of background value.Pixel a (i, j) specific algorithm of interpolation can be expressed as follows:

LST_insert=LST_background+LST_varianceFormula 1-1

${\mathrm{LST}}_{\mathrm{background}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{LST}}_{\mathrm{Gaussian}}$ Formula 1-2

${\mathrm{LST}}_{\mathrm{variance}}=\frac{\mathrm{\&Sigma;}{W}_K\mathrm{\&Delta;}{\mathrm{LST}}_K}{\mathrm{\&Sigma;}{W}_K}$ Formula 1-3

In formula, LST_insertFor a point LST interpolation result, LST_backgroundFor a point background value for many years.

Owing to LST data existing missing values and low quality data, these points should be removed when calculating background value, it is therefore desirable to LST time series data is reconstructed.The present invention selects required setup parameter less asymmetric Gaussian function fitting method that LST time series data is reconstructed, and to the time series data LST after matching_GaussianAsking for background value for many years of each phase, n is that in data set, a certain issue is according to contained year number, as shown in above formula 1-2；If some pixel missing values or low quality data in time series data is too much, cannot be carried out matching, this pixel value is then substituted by the meansigma methods of ground mulching type similar in the periphery radius of influence (in historical years identical at least over half ground mulching type) pixel, the problem wherein with periphery pixel, elevation difference being existed for disappearance pixel, the present invention utilizes height above sea level often raise 1000m temperature to decline the relation of about 6K, remove the elevation participating in the calculating pixel impact on LST.

LST_varianceFor the undulating value at a point place, by ground mulching type and radius of influence decision, Δ LST is the difference of the observation of high-quality pixel and background value in the radius of influence, and K is the high-quality pixel number of identical earth's surface cover type in the radius of influence, W_KFor its respective weights, following formula calculate and obtain:

$W_{\mathrm{ijK}}=\left\{\begin{array}{cc}\frac{R^2-d_{\mathrm{ijK}}^2}{R^2+d_{\mathrm{ijK}}^2}& (d_{\mathrm{ijK}}<R)\\ 0& (d_{\mathrm{ijK}}\&GreaterEqual;R)\end{array}\right.$ Formula 1-4

In formula, d_ijkDistance for the high-quality pixel of interpolation pixel to identical ground surface type；R is the radius of influence, and owing to temperature change within the scope of horizontal range 6000m is generally less than 0.6K, therefore, radius of influence R is set to 3 by the present invention, and namely the high-quality pixel in 7x7 window participates in interpolation calculation.

(2) LST reconstructed operation step

1) data prediction.The present invention chooses MODIS product for draught monitor district LST data reconstruction.MOD11A2, MOD12Q1 product used in reconstruct first passes through MODIS re-projection instrument (MODISReprojectionTool, MRT) re-projection is WGS84 coordinate system, and extracts the data of LST round the clock in product, round the clock LST quality control file and LAI/fPAR system Land Use/Cover Classification result；Research is also used the digital elevation data (DigitalElevationModel, DEM) provided by NASA, first dem data is carried out registration with MODIS data, then be 1000m by the spatial resolution of dem data from 30m resampling.

Used by the present invention, the kind of data and quantity are more, before treatment the studies above data are cut out, consider restructing algorithm relates to spatial window convolution algorithm, therefore select rectangular mask file that data are cut, spatial resolution is 1000m, guarantees that monitored area is in rectangular area central authorities as far as possible.After above-mentioned pretreatment, data used is: LST data, round the clock LST quality control file, dem data and MOD12Q1 data round the clock, data above is carried out Band fusion, constitutes day night LST data set, day night LST quality control file data collection, ground mulching categorical data collection and dem data.

2) many years background values of LST are calculated.Asymmetric Gaussian function fitting will be carried out respectively by LST data set round the clock.Afterwards, the time series data collection after inspection matching: due to, in data set, the pixel value of LST is open type temperature, so the time series all values by successful for non-matching pixel place is set to 0；Thereafter, all pixels are calculated this pixel long-time average annual value same period, as this phase background value, it is thus achieved that background value data set, each 46 wave bands round the clock.It is the pixel of 0 to pixel value in background value data set, it is interpolated in 7x7 window around this pixel: first, with reference to ground mulching categorical data collection, choose pixel identical with interpolation pixel ground mulching type in 7x7 window (exceed with interpolation pixel type same number in historical years half namely think identical earth's surface cover type), if not, choose whole pixel；Subsequently, the weighing computation method in Cressman objective analysis method is utilized to calculate each pixel weight according to choosing the pixel distance with interpolation pixel；Afterwards, utilize dem data, calculate the depth displacement choosing pixel with interpolation pixel, utilize height above sea level often to raise 1000m, the relation of temperature decline 6k, interpolation pixel place elevation temperature is arrived in all LST unifications choosing pixel；Finally, the LST after choosing pixel " levelling " in window is weighted summation with its weight, obtains the many years background value data sets of LST of space and time continuous.

3) LST quality control file data collection round the clock is utilized, respectively LST data round the clock are carried out unit's screening item by item: according to LST quality control file round the clock, only retain high-quality pixel value in each wave band, all the other pixel values are set to 0, it is thus achieved that interpolation LST data set round the clock.

4) calculating gained LST background value data set, interpolation LST data set and ground mulching categorical data collection is utilized to be interpolated, concretely comprise the following steps a. and utilize ground mulching categorical data collection, choose in the pixel periphery 7x7 window that pixel value is 0, it is worth the identical pixel of pixel ground mulching type with 0, if not, chooses whole pixel；B. the weighing computation method in Cressman objective analysis method is utilized to calculate the weight choosing pixel；C. calculate and choose the actual value (i.e. high-quality pixel observation) of pixel and the difference of background value, as undulating value；D. the corresponding weight of the undulating value choosing pixel in window is weighted summation, obtains final LST interpolation result.

3. one according to claim 1 is based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI), it is characterized in that, the structure of described step (3) proportion of crop planting district normalized differential vegetation index-land surface temperature feature space, arable land, draught monitor region is extracted including utilizing, contemporaneous data for many years is utilized jointly to build crops each trophophase proportion of crop planting district normalized differential vegetation index-land surface temperature (being called for short NDVI-LST below) scatterplot, and matching each phase dry and wet limit equation, build NDVI-LST feature space.Specific as follows:

(1) invention utilizes 16 days sintetics MOD13A2 of MODIS vegetation index herein, and the LST8 day data and monitored area crops multiple crop index data on daytime after the inventive method reconstructs, each data spatial resolution is 1000m.

1) MODIS vegetation index product MOD13A2 MODIS re-projection instrument is carried out projection transform, and extract NDVI data therein and quality assessment data (QualityAssurance, QA)；

2) monitoring section data year by year are carried out arable land image element extraction, wherein due in the multiple crop index data of arable land, select the pixel value that pixel area coincide the most with cultivated area data in statistical yearbook, each year, Area distortion controlled below 10%, plough using the pixel of this multiple crop index pixel value as monitored area pixel, is ploughed in monitored area and extract；

3) with reference to 16 grades of criteria for classifications of NDVI mass in QA data, as shown in table 1, the screening NDVI quality of data is medium above pixel, namely at the value of QA data 2～5 pixel less than 4, all the other pixels in NDVI data is assigned to 0；

4) to the NDVI data after QA data screening, carrying out time series reconstruct, reconstructing method selects Savitzky-Golay filter method, and filter window is set to 5；

Every for LST8 day data on daytime 2 scapes are merged into 16 days LST data of 1 scape, in merging, the high-quality pixel of original LST data substitutes interpolation pixel, when being all high-quality pixel or interpolation pixel when two pixels, merge algorithm with reference to MOD11A2 user's manual, to two pixel value averaged, as merging data pixel value.

(2) first, many annual datas are grouped on schedule, NDVI, LST data are divided into some groups, each group comprises the equal number same period NDVI, LST data, all crops pixels of each group of data are built NDVI-LST feature space；Then, one of percentage of difference of, minima maximum with the NDVI value in each group of data is for step-length, and progressively in the long group obtained corresponding to different NDVI value, LST is maximum, minima；Finally utilize the dry and wet limit equation that NDVI and LST is maximum, minima matching is respectively organized in data.

4. one according to claim 1 is based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI), it is characterized in that, the calculating of described step (4) crops temperature vegetation drought index, calculates monitored area C-TVDI value including proposing temperature vegetation drought index based on Price in nineteen ninety.The computational methods of concrete temperature vegetation drought index are expressed as follows:

$\mathrm{TVDI}=\frac{T_s-T_{s\_\min }}{T_{s\_\max }-T_{s\_\min }}$ Formula 1-5

Wherein,

T_{s_max}=a₁+b₁NDVI formula 1-6

T_{s_min}=a₂+b₂NDVI formula 1-7

In formula, T_sFor pixel (i, LST value j), T_{s_max}With T_{s_min}Respectively pixel (i, LST value on the dry limit of NDVI value correspondence j) and wet limit, dry limit, by the scatterplot that the maximum that all kinds of NDVI values in feature space are corresponding is constituted, carries out linear regression and tries to achieve, and wet limit is then tried to achieve by the linear regression of corresponding minima scatterplot；a₁, a₂, b₁, b₂The respectively undetermined coefficient in the equation of dry and wet limit.

5. one according to claim 1 is based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI), it is characterized in that, described step (5) is monitored based on the drought loss of crops temperature vegetation drought index, including the design philosophy based on Supervised classification device, determine the monitored area parameter based on the drought loss monitoring model of crops temperature vegetation drought index by historical data, carry out the monitoring of drought.Specific as follows:

Based on the design philosophy of Supervised classification device, using former years crop growth period remotely-sensed data with actual measurement of each year drought loss data as training sample, by the study to training sample, it is determined that the parameter in drought loss monitoring model；Recycle the crop growth period remotely-sensed data of up-to-date a year afterwards, this year year-end drought loss is monitored.This model is mainly monitored three parts by disaster area evaluation function, parameter optimization and drought loss and is formed:

(1) disaster area evaluation function: disaster area is defined as the Model on Sown Areas of Farm suffering various natural disaster, and suffer several or natural disaster several times in the same year, wherein to endanger maximum once calculating disaster area, do not repeat meter calamity, the maximum Model on Sown Areas of Farm that affected area by drought was then affected by drought for this year one or many.Consider that disaster area value is for drought maximum effect area to crops in year, if considering the cropping pattern that monitoring section is repeatedly multiple cropping, and crops impact of a drought after harvesting terminates, therefore, disaster area evaluation function is divided into different planting season, is set as follows:

S_i=Max [Mean (s_ik..., s_il) ..., Mean (s_im..., s_in)] formula 1-8

In formula, S_iIt is 1 year disaster area estimated value；K～l is harvesting of annual first phase crop Critical growing period, and m～n is annual Final Issue harvesting crop Critical growing period, s_ij(j=k～l) is the area suffered from drought of 1 year first phase harvesting each Critical growing period of crop, s_ij(j=m～n) gathers in the area suffered from drought of each trophophase of crop for Final Issue；Owing to needing the accumulation of certain time from drought of suffering from drought, the damage caused by a drought of single issue evidence may not cause drought, and the Spatial Variability of arid is little in a short time, therefore using the meansigma methods of each phase area suffered from drought as disaster area；Annual disaster area value is then the maximum in each phase disaster area.

Area suffered from drought s_ij(j=k～l or m～n) is by 1 year crops Critical growing period agricultural remote sensing exponent data and this trophophase threshold parameter t_j(j=k～l or m～n) determines.For crops temperature vegetation drought index, exponential quantity is more arid closer to 1 expression, then the area suffered from drought S of 1 year jth Critical growing period_ijFor this year this issue according in C-TVDI value more than this phase threshold value t_jPixel shared by area, threshold value t_jFor model parameter, parameter and trophophase one_to_one corresponding used.

(2) parameter optimization, namely utilizes optimizing algorithm that the object function set is asked for the process being worth most.Drought loss monitoring model carries out the parameter optimization of model using remotely-sensed data for many years and agricultural disaster area data as training sample.The optimal solution of model parameter is regarded as under this group parameter, calculates and obtains disaster area estimated value and the minimum average B configuration deviation in actual disaster area in years of training sample for many years, is therefore set as by object function:

$\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}|S_i-S_{i0}|$ Formula 1-9

In formula, S_iIt is 1 year disaster area estimated value, S_i0Being 1 year actual disaster area value, n is training sample number, and object function is tried to achieve model parameter t corresponding during minima_j(j=k～l or m～n) is optimized parameter.

This research selects grid-search algorithms to carry out parameter optimization, and the speed of searching optimization of grid-search algorithms is very fast, it is possible to obtain globally optimal solution, will not be absorbed in locally optimal solution.But wanting setup parameter scope before search, relatively big at parameter area, in the search that step-size in search is less, time consumption for training is longer.

Parameter optimization process is expressed as follows: the remotely-sensed data in training sample is carried out crops Indices extraction by (i), and selects to participate in the crop growth period that model calculates；(ii) minimum target functional value, each trophophase parameter t are set_i0Hunting zone, step-size in search and parameter t_i0Initial value；(iii) the agricultural remote sensing exponent data of extraction is input in the evaluation function of disaster area, it is thus achieved that disaster area estimated value；(iv) disaster area estimated value is brought in object function with disaster area data in training sample, calculating target function value, if this value is less than current minimum target functional value, then retains parameter current, and minima is updated；V () obtains new parameter value by grid-search algorithms, and repeat above-mentioned (iii), (iv) step, if all parameter all travels through, completes optimizing, records optimized parameter.

(3) drought loss monitoring, namely by crop growth year, remotely-sensed data phase to be monitored, is monitored this year drought loss.Sample number owing to participating in training in this research is less, and training sample is limited to the representativeness of the condition of a disaster in various degree, hence with the linear relationship between training sample disaster area estimated value and actual disaster area, the disaster area estimated value in year to be monitored is adjusted.What drought loss was monitored specifically comprises the following steps that (i) calculates actual disaster area and the equation of linear regression of disaster area estimated value under optimized parameter in training sample；(ii) year to be monitored crops Critical growing period Indices is extracted；(iii) disaster area evaluation function is utilized to ask for the disaster area estimated value under optimized parameter；(iv) using estimated value in previous step as linear equation independent variable, it is thus achieved that disaster area estimated value after adjustment；V () utilizes disaster area estimated value after adjusting to calculate disaster-stricken rate (I), and with reference to drought loss criteria table, it is thus achieved that drought loss monitoring result.Drought loss criteria table is as follows:

Table 2 drought loss criteria table