CN105760978B

CN105760978B - One kind being based on the agricultural drought disaster grade monitoring method of temperature vegetation drought index (TVDI)

Info

Publication number: CN105760978B
Application number: CN201510430499.4A
Authority: CN
Inventors: 李天祺; 朱秀芳; 潘耀忠; 范一大; 李素菊; 王志强; 和海霞
Original assignee: MINISTRY OF CIVIL AFFAIRS NATIONAL DISASTER REDUCTION CENTER; Beijing Normal University
Current assignee: MINISTRY OF CIVIL AFFAIRS NATIONAL DISASTER REDUCTION CENTER; Beijing Normal University
Priority date: 2015-07-22
Filing date: 2015-07-22
Publication date: 2019-11-19
Anticipated expiration: 2035-07-22
Also published as: CN105760978A

Abstract

The invention discloses a kind of agricultural drought disaster grade monitoring methods for being based on temperature vegetation drought index (TVDI), include the following steps: (1) data preparation；(2) land surface temperature (LST) data reconstruction；(3) crop planting area normalized differential vegetation index-land surface temperature (NDVI-LST) feature space building；(4) TVDI is calculated；(5) the drought loss monitoring based on TVDI.The invention proposes the LST data reconstruction methods based on many years background value Yu region undulating value, and for crops building arable land region many years NDVI-LST feature space, calculate temperature vegetation drought index, it devises the crops drought loss monitoring model based on supervised classification thought and carries out drought loss remote sensing monitoring, the model can in real time, relatively accurately reflect the drought stress degree that crop is subject under different condition, be of great significance in the monitoring of agricultural drought disaster, early warning and prevention.

Description

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

Technical field

The present invention relates to one kind using MODIS remotely-sensed data as key data, by reconstruct land surface temperature data and ties It closes normalized differential vegetation index and establishes crops temperature arid vegetation index, the agricultural drought based on crops temperature vegetation drought index Calamity grade monitoring model carries out agricultural arid grade monitoring method, and similarity evaluating model realizes Natural Disaster rapid evaluation Method, specially a kind of agricultural drought disaster grade monitoring method for being based on temperature vegetation drought index (TVDI).

Background technique

Drought is one of most important natural calamity in human production life, and compared with other natural calamities, drought has Occurrence frequency is high, the duration is long, involves the wide feature of range.Arid is as caused by water deficit, the item existing for hazard-affected body Under part, large range of long-time arid can cause drought.Drought is a complexity, the process of accumulation, and performance is mainly Surface water is dry, level of ground water declines, surface vegetation growth is affected, and influence mainly has crop production reduction or total crop failure, Industrial production caused by human livestock drinking water difficulty and continuous high temperature water shortage is obstructed.China is population and large agricultural country, grain peace Full problem is the major issue concerning social stability and economic development.In China, drought endangers agricultural production bring the tightest Weight.Since 1949, the grain disaster area as caused by arid, every year on average up to 3.2 hundred million mu, wherein rate of causing disaster is close 30%, it is not only lost to China's agricultural belt, also seriously constrains social and economic development.

Drought rank is more, complicated mechanism, has and involves range is big, and the duration is long, and space and time upper variability are strong etc. Feature, at present without mature effective method in agricultural drought disaster monitoring；Traditional meteorological drought monitoring: meteorological site is sparse, And be unevenly distributed, data space precision is not high, damage caused by a drought and agricultural disaster effectively cannot be established corresponding relationship；Agricultural drought disaster calamity Feelings statistics: although the Disaster degree of agricultural arid can be obtained more accurately, the statistical data acquisition period is longer, expends a large amount of Manpower and material resources, and data, mostly as unit of administrative division, spatial accuracy is limited, in addition, there is also Statistical Criteria disunities etc. to ask Topic.In view of the above-mentioned problems, remote sensing technology has the characteristics that multi-platform, multisensor, high-resolution, it can quickly, accurately Damage caused by a drought information is obtained, repeated observation can be carried out to devastated in a short time；According to the spectral characteristic of ground of different-waveband, Remote sensing technology can directly or indirectly obtain the information such as surface water variation, vegetation growth, soil moisture content.Compared to meteorological drought Monitoring and agricultural drought disaster the statistics of geological disaster situation, agricultural drought disaster remote sensing monitoring with high precision, on a large scale, rapidly can obtain soil moisture and become Change and crop growthing state, can in real time, accurately reflect the drought stress degree that crop is subject under different condition, in agriculture drought It is of great significance in monitoring, early warning and the prevention of calamity.

Summary of the invention

Remote sensing technology can quickly, repeatedly, on a large scale obtain earth's surface information, be of great significance in Monitoring of drought. The temperature vegetation drought index (TVDI) for wherein utilizing optics and IRMSS thermal band, has been effectively combined the finger of vegetation state Show factor normalized differential vegetation index (NDVI), it is dry in remote sensing with the thermally equilibrated important parameter land surface temperature (LST) of surface water It is widely used in non-irrigated study on monitoring.But there are the following problems in the monitoring of practical Agriculture Drought by TVDI: (1) missing values are more makes LST data are discontinuous on space-time, can not further extract TVDI；(2) it is ploughed in the building of NDVI-LST feature space there are non- The influence of ground pixel, and traditional TVDI lacks comparativity between many years data.For two problems as above, the invention proposes LST data reconstruction method based on many years background value Yu region undulating value, compared to existing reconstructing method, in restructural list scape image The space missing value of large area, and the shortage of data of single pixel long-term sequence, and the variation of LST can be retained to a certain extent Details；Using crops as research object, building arable land region many years NDVI-LST feature space, and propose crops temperature vegetation Drought index (C-TVDI) herein on basis, calculates the C-TVDI of monitoring section crops Critical growing period, and design based on prison The crops drought loss monitoring model for superintending and directing classificating thought carries out Henan Province's drought loss remote sensing monitoring, and to fight calamities and provide relief, decision is mentioned For reference.

To achieve the above object the present invention adopts the following technical scheme: one kind is based on temperature vegetation drought index (TVDI) Agricultural drought disaster grade monitoring method, includes the following steps:

(1) preparation of data:

The present invention carries out agricultural drought disaster study on monitoring using the data that remote sensing technology obtains, wherein remotely-sensed data used is The MODIS LST product of NASA (http://ladsweb.nascom.nasa.gov/) offer, vegetation index product, soil benefit With covering product and dem data, and arable land multiple crop index data；Statistical data is the cultivated area number in China Statistical Yearbook The agriculture provided according to, disaster area data and planting industry management department, the Ministry of Agriculture (http://www.zzys.moa.gov.cn/) Agricultural crop sown area data and crops phenological calendar.

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

MODIS LST data product based on acquisition according to corresponding method to land surface temperature data (hereinafter referred to as LST it) is reconstructed.

(3) crop planting area normalized differential vegetation index-land surface temperature feature space building:

Draught monitor region arable land is extracted, constructs each growth period crop planting of crops jointly using many years contemporaneous data Area's normalized differential vegetation index-land surface temperature (hereinafter referred to as NDVI-LST) scatter plot, and each phase is fitted in wet side equation, structure NDVI-LST feature space is built in construction.

(4) calculating of crops temperature vegetation drought index:

It is based on step 3 as a result, proposing in nineteen ninety temperature vegetation drought index using Price calculates monitoring region Crop planting area crops temperature vegetation drought index (hereinafter referred to as C-TVDI).

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

Based on the design philosophy of Supervised classification device, determine that monitoring region is based on crops temperature vegetation by historical data The parameter of the drought loss monitoring model of drought index, carries out the monitoring of drought.

As a further solution of the present invention, the land surface temperature based on background value and undulating value of the step (2) Data reconstruction includes reconstructing method and reconstructed operation step, and specific method and operating process are as follows:

(1) LST reconstructing method

The present invention is based in Cressman objective analysis method using initial fields value with correct the think of that value approaches observation jointly Think, by many years average level of pixel LST, i.e. the background value initial fields value that is regarded as pixel；It is by the undulating value of pixel LST, i.e., all The observation of pixel and the difference of background value once correct interpolation pixel LST with this as value is corrected in the radius of influence of side Interpolation.The specific algorithm of pixel a (i, j) interpolation can be expressed as follows:

LST_insert=LST_background+LST_varianceFormula 1-1

Formula 1-2

Formula 1-3

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

Since there are missing values and low quality data in LST data, these points should be removed when calculating background value, therefore It needs that LST time series data is reconstructed.The less asymmetric Gaussian function fitting of setup parameter needed for the present invention selects LST time series data is reconstructed in method, and to the time series data LST after fitting_GaussianSeek each many years phase background Value, n is that a certain issue is according to contained year in data set, as shown in above formula 1-2；If certain pixels lack in time series data Value or low quality data are excessive, can not be fitted, the pixel value is then by ground mulching type (In similar in the radius of influence of periphery It is identical at least over half ground mulching type in historical years) average value of pixel substitution, wherein for missing pixel and week There is elevation difference in side pixel, utilize the every relationship for increasing 1000m temperature decline about 6K of height above sea level, removal ginseng in the present invention With influence of the elevation to LST for calculating pixel.

LST_varianceIt for the undulating value at a point, is determined by ground mulching type and the radius of influence, Δ LST is the radius of influence The observation of interior high quality pixel and the difference of background value, K are the high quality pixel of identical earth's surface cover type in the radius of influence Number, W_KFor its respective weights, it is calculated by following formula:

Formula 1-4

In formula, d_ijkFor the distance of the high quality pixel of interpolation pixel to identical ground surface type；R is the radius of influence, due to temperature The variation spent within the scope of horizontal distance 6000m is generally less than 0.6K^[57], therefore, will affect radius R herein and be set as 3, that is, exist High quality pixel in 7x7 window participates in interpolation calculation.

(2) LST reconstructed operation step

1) data prediction.The present invention chooses MODIS product and is used for draught monitor area LST data reconstruction.It is used in reconstruct To MOD11A2, MOD12Q1 product pass through first MODIS re-projection tool (MODIS Reprojection Tool, MRT) weight It is projected as WGS84 coordinate system, and extracts the data of LST round the clock in product, round the clock LST mass control file and LAI/fPAR body It is Land Use/Cover Classification result；Digital elevation data (the Digital provided by NASA is also provided in research Elevation Model, DEM), first dem data is registrated with MODIS data, then by the spatial resolution of dem data from 30m resampling is 1000m.

Since the type and quantity of studying data used are more, the studies above data are cut out before treatment, are considered Spatial window convolution algorithm involved in restructing algorithm, therefore rectangular mask file is selected to cut data, spatial resolution For 1000m, ensures to monitor region as far as possible and be in rectangular area center.The data used after above-mentioned pretreatment are as follows: round the clock LST data, the file of LST mass control round the clock, dem data and MOD12Q1 data, carry out Band fusion to above data, constitute Day night LST data set, day night LST mass control file data collection, ground mulching categorical data collection and dem data.

2) more years background values of LST are calculated.Asymmetric Gaussian function fitting will be carried out respectively by LST data set round the clock.Later, it examines Time series data collection after testing fitting: since the pixel value of LST is open type temperature in data set, so success will be fitted Pixel where time series all values be set as 0；Thereafter, pixel same period long-time average annual value is calculated to all pixels, as The phase background value obtains background value data set, round the clock each 46 wave bands.The pixel for being 0 to pixel value in background value data set, In Interpolation is carried out around the pixel in 7x7 window: firstly, referring to ground mulching categorical data collection, being chosen in 7x7 window and interpolation The identical pixel of pixel ground mulching type (is thought with interpolation pixel type same number more than half in historical years Identical earth's surface cover type), if choosing whole pixels without if；Then, the weight meter in Cressman objective analysis method is utilized Calculation method calculates each pixel weight according to pixel is chosen at a distance from interpolation pixel；Later, it using dem data, calculates and chooses The depth displacement of pixel and interpolation pixel, using the every raising 1000m of height above sea level, temperature declines the relationship of 6k, by all selection pixels The unified pixel to interpolation of LST where elevation temperature；Finally, to LST and its weight after selection pixel " levelling " in window It is weighted summation, obtains the more years background value data sets of LST of space and time continuous.

3) file data collection is controlled using LST mass round the clock, member item by item is carried out to LST data round the clock respectively and is screened: foundation LST mass controls file round the clock, only retains high quality pixel value in each wave band, remaining pixel value is set as 0, is obtained to be inserted round the clock Value LST data set.

4) using calculate gained LST background value data set, interpolation LST data set and ground mulching categorical data collection into Row interpolation, the specific steps are a. to utilize ground mulching categorical data collection, chooses in the pixel periphery 7x7 window that pixel value is 0, with The identical pixel of 0 value pixel ground mulching type, if choosing whole pixels without if；B. using in Cressman objective analysis method Weighing computation method calculates the weight for choosing pixel；C. the actual value (i.e. high quality pixel observation) and back for choosing pixel are calculated The difference of scape value, as undulating value；D. summation is weighted to the corresponding weight of undulating value for choosing pixel in window, obtained most Whole LST interpolation result.

As a further solution of the present invention, step (3) crop planting area normalized differential vegetation index-land surface The building in temperature profile space constructs farming using many years contemporaneous data including the use of draught monitor region arable land is extracted jointly Each growth period crop planting of object area normalized differential vegetation index-land surface temperature (hereinafter referred to as NDVI-LST) scatter plot, and It is fitted each phase dry and wet side equation, building building NDVI-LST feature space.It is specific as follows:

(1) invention utilizes 16 days sintetics MOD13A2 of MODIS vegetation index herein, after the method for the present invention reconstructs Daytime LST 8 day data and monitoring region crops multiple crop index data, each data spatial resolution is 1000m.

1) MODIS vegetation index product MOD13A2 is subjected to projection transform with MODIS re-projection tool, and extracted wherein NDVI data and quality assessment data (QualityAssurance, QA)；

2) to monitoring section, data carry out arable land image element extraction year by year, wherein due in the multiple crop index data of arable land, selection Cultivated area data pixel value the most identical in pixel area and statistical yearbook, each year Area distortion control 10% hereinafter, Using the pixel of the multiple crop index pixel value as monitoring region arable land pixel, monitoring region arable land is extracted；

3) referring to 16 grades of classification standards of NDVI mass in QA data, as shown in table 1, the screening NDVI quality of data is medium Remaining pixel in NDVI data is assigned to 0 that is, in QA data 2~5 pixels of the value less than 4 by above pixel；

4) time series reconstruct, reconstruct side are carried out to the NDVI data after QA data screening using Timesat3.1.1 Method selects Savitzky-Golay filter method, and filter window is set as 5；

5) every 2 scape of 8 day data of LST on daytime is merged into 1 scape, 16 days LST data, original LST data is high-quality in merging It measures pixel and substitutes interpolation pixel, when two pixels are all high quality pixel or interpolation pixel, merge algorithm and refer to MOD11A2 user Handbook, to two pixel value averageds, as merging data pixel value.

(2) firstly, be grouped to many years data on schedule, NDVI, LST data is divided into several groups, include identical number in each group The same period NDVI, LST data are measured, NDVI-LST feature space is constructed to all crops pixels of each group of data；Then, with each group 1 percent of maximum, minimum value the difference of NDVI value in data are step-length, gradually long to obtain group corresponding to different NDVI values Interior LST maximum, minimum value；Dry and wet side equation finally maximum using NDVI and LST, in minimum value fitting each group of data.

As a further solution of the present invention, the calculating of step (4) the crops temperature vegetation drought index, including base Temperature vegetation drought index, which is proposed, in nineteen ninety in Price calculates monitoring region C-TVDI value.Specific temperature vegetation arid The calculation method of index is expressed as follows:

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 the LST value of pixel (i, j), T_{s_max}With T_{s_min}The corresponding dry side of the NDVI value of respectively pixel (i, j) With the LST value on wet side, dry side is carried out linear by the scatter plot constituted to the corresponding maximum value of NDVI values all kinds of in feature space Recurrence acquires, and wet side is then acquired by the linear regression of corresponding minimum value scatter plot；a₁, a₂, b₁, b₂Respectively dry and wet side equation In undetermined coefficient.

As a further solution of the present invention, the drought etc. of the step (5) based on crops temperature vegetation drought index Grade monitoring determines that monitoring region is based on crops temperature by historical data including the design philosophy based on Supervised classification device The parameter of the drought loss monitoring model of vegetation drought index, carries out the monitoring of drought.It is specific as follows:

Based on the design philosophy of Supervised classification device, with the remotely-sensed data of former years crop growth period and actual measurement of each year drought etc. Grade data determine the parameter in drought loss monitoring model by the study to training sample as training sample；It is sharp again later With newest 1 year crop growth period remotely-sensed data, the year-end drought loss in this year is monitored.The model is mainly by disaster-stricken face Product evaluation function, parameter optimization and drought loss monitoring three parts composition:

(1) disaster area evaluation function: disaster area is defined as the Model on Sown Areas of Farm by various natural calamities, and In the same year by several or several times natural calamity does not repeat meter calamity, drought wherein to endanger maximum primary calculating disaster area Calamity disaster area is then this year one or many maximum Model on Sown Areas of Farm influenced by drought.Consider that disaster area value is Drought is to the biggest impact areas of crops in year, considers monitoring section if it is the cropping pattern of multiple multiple cropping, and crops pass through The influence of a drought terminates after harvesting, therefore, disaster area evaluation function is divided into different planting seasons, is set as follows:

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

In formula, S_iFor 1 year disaster area estimated value；K~1 is the annual first phase to gather in crop Critical growing period, m~n Crop Critical growing period, s are gathered in for annual Final Issue_ij(j=k~1) is 1 year first phase to gather in each Critical growing period of crop Area suffered from drought, s_ij(j=m~n) is the area suffered from drought that Final Issue gathers in crop each growth period；Due to from suffering from drought drought The accumulation of certain time is needed, the damage caused by a drought of single issue evidence may not cause drought, and Spatial Variability arid in a short time is not Greatly, therefore using the average value of each phase area suffered from drought as disaster area；Annual disaster area value is then each phase disaster area Maximum value.

Area suffered from drought s_ij(j=k~1 or m~n) was by 1 year crops Critical growing period agricultural remote sensing exponent data and was somebody's turn to do Growth period threshold parameter t_j(j=k~1 or m~n) is determined.By taking crops temperature vegetation drought index as an example, index value is closer 1 expression is more arid, then the area suffered from drought S of 1 year j-th Critical growing period_ijFor the year issue, C-TVDI value is greater than in Phase threshold value t_jPixel occupied area, threshold value t_jFor model parameter, parameter and growth period used are corresponded.

(2) parameter optimization seeks the process being most worth using objective function of the optimizing algorithm to setting.Drought loss monitoring Model carries out the parameter optimization of model using many years remotely-sensed data and agricultural disaster area data as training sample.Model parameter Optimal solution is regarded as under this group of parameter, and many years disaster area estimated value is calculated and reality in years of training sample is disaster-stricken The minimum average B configuration deviation of area, therefore objective function is set are as follows:

Formula 1-9

In formula, S_iFor 1 year disaster area estimated value, S_i0For 1 year practical disaster area value, n was training sample number, Corresponding model parameter t when acquiring minimum value to objective function_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, can obtain Globally optimal solution is obtained, locally optimal solution will not be fallen into.But setup parameter range is wanted before search, search larger in parameter area In the lesser search of step-length, time consumption for training is longer.

Parameter optimization process is expressed as follows: the remotely-sensed data in training sample is carried out the extraction of crops Indices by (i), And it selects to participate in the crop growth period that model calculates；(ii) minimum target functional value, each growth period parameters t are set_i0Search model It encloses, step-size in search and parameter t_i0Initial value；(iii) the agricultural remote sensing exponent data of extraction is input to disaster area estimation letter In number, disaster area estimated value is obtained；(iv) disaster area data in disaster area estimated value and training sample are brought into target In function, calculating target function value retains parameter current if the value is less than current minimum target functional value, and by minimum value It updates；(v) new parameter value is obtained by grid-search algorithms, and repeats above-mentioned (iii), (iv) step, if whole parameters are equal It has been traversed that, complete optimizing, recorded optimized parameter.

(3) drought loss monitors, i.e., by year crop growth phase remotely-sensed data to be monitored, carries out to this year drought loss Monitoring.Since the sample number for participating in trained in this research is less, training sample is limited to the representativeness of different degrees of the condition of a disaster, therefore Using the linear relationship between training sample disaster area estimated value and practical disaster area, the disaster area in year to be monitored is estimated Value is adjusted.Specific step is as follows for drought loss monitoring: (i) calculates practical disaster area and optimized parameter in training sample The equation of linear regression of lower disaster area estimated value；(ii) year crops Critical growing period Indices to be monitored are extracted；(iii) The disaster area estimated value under optimized parameter is sought using disaster area evaluation function；(iv) using estimated value in previous step as line Property equation independent variable, disaster area estimated value after being adjusted；(v) disaster-stricken rate is calculated using disaster area estimated value after adjustment (I), table and referring to drought loss standard is divided, obtains drought loss monitoring result.It is as follows that drought loss standard divides table:

2 drought loss standard of table divides table

Compared with prior art, the present invention has carried out following improvement: (1) proposing and fluctuated based on many years background value and region The LST data reconstruction method of value compares existing reconstructing method, the space missing value of large area in restructural list scape image, and single The shortage of data of pixel long-term sequence, and the variation details of LST can be retained to a certain extent；It (2) is research with crops Object, building arable land region many years NDVI-LST feature space, and propose crops temperature vegetation drought index (C-TVDI), In On this basis, the C-TVDI of monitoring section crops Critical growing period is calculated, and design the crops based on supervised classification thought Drought loss monitoring model carries out Henan Province's drought loss remote sensing monitoring, and to fight calamities and provide relief, decision provides reference.

Detailed description of the invention

Fig. 1 is the agriculture that one kind of specification and specific embodiment is based on temperature vegetation drought index (TVDI) according to the present invention The land surface temperature data reconstruction Technology Roadmap of industry drought loss monitoring method；

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

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

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

Fig. 5 is that one kind of specific embodiment according to the present invention is based on the agricultural drought disaster etc. of temperature vegetation drought index (TVDI) The many years crops NDVI-LST part scatter plot of grade monitoring method；

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

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

Specific embodiment

The present invention is further elaborated in the following with reference to the drawings and specific embodiments.

Present case is research area with Henan, mainly uses vegetation index data product, LST data product, land cover pattern/soil Ground covers delta data product.All kinds of MODIS products used herein are as follows:

1 land MODIS product of table

Arable land multiple crop index data are mainly used for obtaining the spatial distribution of Henan arable land pixel, and the data are by doctor Liu Jianhong It provides, data spatial resolution 500m, present case is used, and for Henan area 2001-2011, totally 11 scape data, pixel value are The crop planting number of the annual pixel, bare place pixel are labeled as 0, and the arable land specific extraction algorithm of multiple crop index data is such as Under:

(1) whole nation is divided by a ripe area, the area Liang Shu and three ripe areas according to many years averagely meteorological data；

(2) training sample appropriate is selected in each shortening area, farming in each shortening area is determined according to training sample The most short Length of growing season of object, longest Length of growing season, minimum growth amplitude；

(3) using MODIS data obtain each pixel enhancement mode meta file (Enhanced Vegetation Index, EVI) time-serial position, and objects marquis's parameter such as extract vegetation growing season number, Length of growing season, growth amplitude；

(4) by the most short Length of growing season of the crops where the object marquis parameter of image element extraction and pixel in shortening area, most Long Length of growing season, minimum growth amplitude are compared, and differentiate whether a vegetation growing season belongs to crop growth season, if It is to retain, if not then deleting；

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

In present case, statistical yearbook data have also been used, farming counts evidence one by one and all kinds of 3 classes of Model on Sown Areas of Farm data Statistical data；Wherein farming counts evidence and all kinds of Model on Sown Areas of Farm data one by one, for determining the critical developmental of staple crops Phase, and screen and be used for agricultural drought disaster study on monitoring with the consistent remotely-sensed data of crops Critical growing period；Statistical yearbook data structure Remote sensing drought loss monitoring model is built to test with to model monitoring result.Studying statistical yearbook data used includes Henan Province Each year affected area by drought data and cultivated area data, the data all are from 2001-2011 China Statistical Yearbook；Farming Count one by one evidence and all kinds of Model on Sown Areas of Farm data from planting industry management department, the Ministry of Agriculture (http: // Www.zzys.moa.gov.cn/), according to the sown area of all kinds of crops, Crops in Henan Province is winter wheat, Xia Yu Rice, semilate rice, soybean, cotton, rape and peanut, farming of all kinds of staple crops in the research time limit, which is gone through, to be as follows:

Present case main flow includes: the LST data reconstruction method for 1) utilizing many years background value and region undulating value, to river Southern 2000-2011 LST data are reconstructed, 5000 high quality pixels of random selection, carry out precision to reconstruct data and test Card；2) using the LST data and MODISNDVI data after reconstruct, Henan area crops NDVI-LST feature space is constructed；3) Using 2001-2007 crops temperature vegetation drought index data and practical disaster area data as training sample, with 2008- Crops temperature vegetation drought index data in 2011, practical disaster area data and drought loss data are as inspection sample This, evaluates crops drought monitoring model.

Present case is chosen MODIS product and is studied for Henan area LST data reconstruction.Used MOD11A2 in reconstruct, It is WGS84 that MOD12Q1 product passes through MODIS re-projection tool (MODIS Reprojection Tool, MRT) re-projection first Coordinate system, and extract the data of LST round the clock in product, the file and LAI/fPAR system land use of LST mass control round the clock/ Cover classification result；Also used in research by NASA provide digital elevation data (Digital Elevation Model, DEM), first dem data is registrated with MODIS data, then is from 30m resampling by the spatial resolution of dem data 1000m。

Since the type and quantity of studying data used are more, the studies above data are cut out before treatment, are considered Spatial window convolution algorithm involved in restructing algorithm, therefore rectangular mask file is selected to cut data, it is sized to 641x622, spatial resolution 1000m, research area Henan are in rectangular area center.The research number used after above-mentioned pretreatment According to are as follows: each 506 scape of LST data round the clock, LST mass controls each 506 scape of file round the clock, wherein annual 46 scape, totally 11 years；DEM number According to 1 scape；Totally 11 scape, annual 1 scape carry out Band fusion to above data, constitute day night LST data set, day night MOD12Q1 LST mass controls file data collection, ground mulching categorical data collection and dem data, convenient for the processing of further restructing algorithm.

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

(1) more years background values of LST are calculated.Asymmetric Gaussian function fitting will be carried out respectively by LST data set round the clock, the step It is realized by Timesat3.1.1.Later, the time series data collection after fitting is examined: since the pixel value of LST is in data set Open type temperature, so the time series all values where not being fitted successful pixel are set as 0；Thereafter, all pixels are calculated The pixel same period long-time average annual value obtains background value data set, round the clock each 46 wave bands as the phase background value.To background value The pixel that pixel value is 0 in data set carries out interpolation in 7x7 window around the pixel: firstly, referring to ground mulching number of types According to collection, the interior pixel identical with interpolation pixel ground mulching type of 7x7 window is chosen (in 11 Nian Zhongyu interpolation pixel types Same number is more than to think identical earth's surface cover type 6 times), if choosing whole pixels without if；Then, Cressman is utilized Weighing computation method in objective analysis method calculates each pixel weight according to pixel is chosen at a distance from interpolation pixel；Later, Using dem data, the depth displacement for choosing pixel and interpolation pixel is calculated, using the every raising 1000m of height above sea level, temperature declines 6k's Relationship, by elevation temperature where the unified pixel to interpolation of all LST for choosing pixel；Finally, " looking for pixel is chosen in window It is flat " after LST and its weight be weighted summation, obtain the more years background value data sets of LST of space and time continuous.

(2) file data collection is controlled using LST mass round the clock, member item by item is carried out to LST data round the clock respectively and is screened: according to The division of pixel credit rating in 3.1 sections only retains high quality pixel value in each wave band, remaining pixel value is set as 0, obtains daytime Night interpolation LST data set, each 506 scape.

(3) calculating gained LST background value data set, interpolation LST data set and ground mulching categorical data collection are utilized Interpolation is carried out, the specific steps are 1. to utilize ground mulching categorical data collection, it chooses in the pixel periphery 7x7 window that pixel value is 0, Pixel identical as 0 value pixel ground mulching type, if choosing whole pixels without if；2. using in Cressman objective analysis method Weighing computation method calculate choose pixel weight；3. calculate choose pixel actual value (i.e. high quality pixel observation) with The difference of background value, as undulating value；4. the corresponding weight of undulating value for choosing pixel in pair window is weighted summation, obtain Final LST interpolation result (Figure of description 2).

Present case is interpolation precision of the verifying based on background value Yu undulating value LST data reconstruction algorithm, design accuracy verifying Experiment, has randomly selected 250 scape LST data in the original data set of LST round the clock, the random choosing again in 250 scape data of selection 5000 high quality pixels are taken, pixel value is assigned to 0, generates LST validation data set round the clock, validation data set and original LST number According to collection in addition to pixel value is by modification, remaining is all the same.Validation data set is carried out based on background value and undulating value LST data weight Structure, and the interpolation result in reconstruction result is compared with original pixel value.In comparing result, LST interpolation and original pixel The maximum deviation of value is 15.44K, and average deviation 0.81K, deviation is more than that the pixel of 2K accounts for the 5.2% of pixel sum；As it can be seen that Interpolation precision based on background value and undulating value LST data reconstruction algorithm is higher, can preferably be repaired to MODIS LST data It is multiple.Figure of description 4 is the scatter plot of interpolation result and initial data, and original pixel value has stronger linear phase with interpolation result Guan Xing, wherein linear gradient 0.986, intercept 1.21, R2 0.960.(Figure of description 2, Fig. 3, Fig. 4)

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

(1) MODIS vegetation index product MOD13A2 is subjected to projection transform with MODIS re-projection tool, and extracted wherein NDVI data and quality assessment data (Quality Assurance, QA)；

(2) to research area, data carry out arable land image element extraction year by year, wherein due in the multiple crop index data of arable land, pixel Value for 2 pixel area and Henan statistical yearbook in cultivated area data the most coincide, each year Area distortion 5%~10% not Deng, therefore pixel of being ploughed using the pixel that multiple crop index value is 2 as research area, research area arable land is extracted；

(3) referring to 16 grades of classification standards of NDVI mass in QA data, as shown in table 6, during the screening NDVI quality of data is Remaining pixel in NDVI data is assigned to 0 that is, in QA data 2~5 pixels of the value less than 4 by the pixel Deng more than；

(4) time series reconstruct, reconstruct side are carried out to the NDVI data after QA data screening using Timesat3.1.1 Method selects Savitzky-Golay filter method, and filter window is set as 5；

(5) every 2 scape of 8 day data of LST on daytime is merged into 1 scape, 16 days LST data, original LST data is high-quality in merging It measures pixel and substitutes interpolation pixel, when two pixels are all high quality pixel or interpolation pixel, merge algorithm and refer to MOD11A2 user Handbook, to two pixel value averageds, as merging data pixel value.

Present case utilizes pretreated NDVI, LST data, empty to Henan crop planting area building NDVI-LST feature Between.Firstly, be grouped to many years data on schedule, NDVI, LST data are divided into 23 groups, include 11 the scape same period NDVI, LST in each group Data construct NDVI-LST feature space to all crops pixels of each group of data；Then, with the NDVI value in each group of data 1 the percent of maximum, minimum value difference are step-length, and gradually LST is maximum, minimum in group corresponding to the different NDVI values of long acquisition Value；Dry and wet side equation finally maximum using NDVI and LST, in minimum value fitting each group of data, partial results such as specification are attached Shown in Fig. 5.

In more phase results, the distribution of many years crops NDVI-LST scatterplot is broadly divided into three classes: where most of period Many years crops NDVI-LST scatterplot distribution meet theoretic angular distribution relationship, such as Figure of description 5 (a) o. 11th (6 Month last ten-days period), 5 (b) the 18th phases (mid-October), 5 (c) the 21st phases (early November), wherein dry side slope is respectively less than 0, and wet side And it is not all the straight line for being parallel to NDVI axis, such as in 5 o. 11th of Figure of description, 18 interim, wet side slope value is all larger than 0, In The interim slope of Figure of description 5 the 21st is less than 0, in terms of overall distribution, the wet slope value when slope value is all larger than dry, and dry and wet side It is gradually drawn close with NDVI growth；The many years crops NDVI-LST scatter plot distributions of some period and theoretical Triangle-Profile phase Instead, as shown in the 8th phase of Figure of description 5 (d) (the first tenday period of a month in May), 5 (e) the 13rd phases (late July), wherein dry side slope is all larger than 0, wet side slope is respectively less than 0, and dry and wet side gradually deviates from the growth of NDVI；Separately there is many years crops NDVI-LST in individual periods Scatter plot distributions are a rectangle, and dry and wet side is parallel to each other, and slope is all close to 0, if the 3rd phase of Figure of description 5 (f) is (in 2 months Ten days).

In the above results, the difference of many years crops NDVI-LST scatterplot distribution is mainly as corresponding to scatter plot each phase The different of crop growth phase determine: in the scatter plot of theoretical Triangle-Profile, such as 11 phases, 18 phases and 21 phases, the time is respectively 6 It the last ten-days period moon, mid-October and early November, respectively corresponds as autumn grain crops sowing/seeding stage, autumn grain crops harvest time and summer grain crops sowing time, In The above three period, crop planting region NDVI Distribution value is more uniform, and bare area, low nurse crop and high nurse crop are total It deposits, therefore NDVI-LST scatter plot is close with theoretical triangle；And it is in the scatter plot of anti-triangle, and such as 8 phases, 13 phases, the time point Not Wei the first tenday period of a month in May and late July, respectively correspond as summer grain crops growth period and autumn grain crops growth period, the two period Grain Growth Situations reach To maximum, most pixels in crop planting area are high nurse crop, only a small number of pixels be bare area or be bare area with The mixed pixel of crop, in NDVI-LST scatter plot, point focuses mostly in the middle high level region of NDVI, and NDVI low value region is only few Several, it is limited that corresponding LST value does not have representative and variation, therefore dry and wet side is more close in NDVI low value region； In the scatter plot of distributed rectangular, such as 3 phases, the time is mid-February, which is the Wintering Period of summer grain crops winter wheat, at this time crops Planting area pixel has certain vegetation coverage, and scatter plot midpoint concentrates on the middle low value area of NDVI, but makees in Wintering Period Object transpiration is very weak, does not have the ability for adjusting canopy surface temperature, therefore dry and wet Bian Juncheng is parallel to the straight line of NDVI axis.

The dry and wet side equation that present case is acquired using the above process further calculates Henan many years crop planting region C-TVDI value, calculation method are as follows:

Formula 1

Wherein,

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

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

Present case calculates the C-TVDI value part result and the same year agricultural disaster area in Henan many years crop planting region Administrative division distribution map is shown in Figure of description 6.

(specification is attached with C-TVDI index spatial distribution comparing result for each city disaster area spatial distribution in present case Henan Province Fig. 6) as it can be seen that C-TVDI can correctly reflect the spatial distribution of agricultural drought disaster on the whole, but on prefecture-level region Accuracy is insufficient, and reason may be codetermined for drought the condition of a disaster by many factors；For C-TVDI is further used for agriculture drought Calamity remote sensing monitoring, and by remotely-sensed data and the opening relationships of agricultural drought disaster the condition of a disaster, present case proposes a kind of based on more phase remotely-sensed datas (C-TVDI) drought loss monitoring model, and model monitoring result is tested in conjunction with Henan Province's affected area by drought data Card.

Present case drought loss monitoring model is the design philosophy based on Supervised classification device, with former years crop growth period Remotely-sensed data and each year survey drought loss data as training sample and determine drought loss by the study to training sample Parameter in monitoring model；The crop growth period remotely-sensed data for recycling newest 1 year later, to the year-end drought loss in this year into Row monitoring.The model is mainly made of disaster area evaluation function, parameter optimization and drought loss monitoring three parts:

(1) disaster area evaluation function: disaster area is defined as the Model on Sown Areas of Farm by various natural calamities, and In the same year by several or several times natural calamity does not repeat to count wherein to endanger maximum primary calculating disaster area Calamity^[65], affected area by drought is then this year one or many maximum Model on Sown Areas of Farm influenced by drought.Consider disaster-stricken Area value is drought in year to the biggest impact area of crops, and research area Henan is mainly the cropping pattern of multiple cropping twice, and Crops influence of a drought after gathering in terminates, therefore, by disaster area evaluation function be divided into summer grain crops disaster area and Autumn grain crops disaster area two parts are set as follows:

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

In formula, S_iFor 1 year disaster area estimated value；K~1 is summer grain crops crop Critical growing period, 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 autumn The area suffered from drought in grain crop each growth period；Due to needing the accumulation of certain time from suffering from drought drought, single issue evidence Damage caused by a drought may not cause drought, and Spatial Variability arid in a short time is little, therefore being averaged each phase area suffered from drought of summer/autumn grain crops Value is as summer/autumn grain crops disaster area；Annual disaster area value is then the maximum value of summer grain crops Yu autumn grain crops disaster area.

Formula 5

Common optimizing algorithm has genetic algorithm, simulated annealing, particle swarm algorithm and grid-search algorithms.Due to It is less that this paper model is related to parameter, and training sample amount is smaller, therefore grid-search algorithms is selected to carry out parameter optimization, and grid is searched The speed of searching optimization of rope algorithm is very fast, can obtain globally optimal solution, will not fall into locally optimal solution.But ginseng is set before search Number range, larger in parameter area, in the lesser search of step-size in search, time consumption for training is longer.

(3) drought loss monitors, i.e., by year crop growth phase remotely-sensed data to be monitored, carries out to this year drought loss Monitoring.Since the sample number for participating in trained in this research is less, training sample is limited to the representativeness of different degrees of the condition of a disaster, therefore Using the linear relationship between training sample disaster area estimated value and practical disaster area, the disaster area in year to be monitored is estimated Value is adjusted.Specific step is as follows for drought loss monitoring: (i) calculates practical disaster area and optimized parameter in training sample The equation of linear regression of lower disaster area estimated value；(ii) year crops Critical growing period Indices to be monitored are extracted；(iii) The disaster area estimated value under optimized parameter is sought using disaster area evaluation function；(iv) using estimated value in previous step as line Property equation independent variable, disaster area estimated value after being adjusted；(v) disaster-stricken rate is calculated using disaster area estimated value after adjustment (I), table and referring to drought loss standard is divided, obtains drought loss monitoring result^[64].It is as follows that drought loss standard divides table:

4 drought loss standard of table divides table

The registration monitoring of present case drought is using the Henan Province 2001-2011 crops temperature vegetation drought index data in Statistical yearbook Henan Province, state disaster area, cultivated area data, evaluate drought loss monitoring model；According to Henan summer and autumn Grain chief crop, winter wheat and summer corn day part the different demands to moisture, by the jointing of winter wheat, heading and grouting The pumping of phase and summer corn is male and milk stage is as the crops Critical growing period in drought loss monitoring model, the growth period Correspond to 6 phase remotely-sensed datas in annual, the corresponding threshold parameter tj (j=1~6) of each issue of data.For 6 thresholds for determining model Value parameter carries out parameter optimization and determines the initial estimated value in disaster area using 2001-2007 annual data as model training sample With the regression equation in practical disaster area；For the monitoring result of evaluation model, using 2008-2011 annual data as test samples, Each year disaster area is estimated using the optimized parameter and equation of linear regression of model, and drought loss is monitored, And contrastive detection result and practical drought loss.Detailed process is shown in Figure of description 7.

Present case using the Henan Province training sample 2001-2007 crops temperature vegetation drought index data and Henan by Calamity area data is as follows to each main growing period threshold parameter optimizing result of drought loss monitoring model:

Each phase threshold value optimizing result of 5 model of table

Training sample disaster area estimated value is acquired using optimized parameter, and is linearly returned with practical disaster area value Return, regression equation and the coefficient of determination are as follows:

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

R²=0.9370

According to above-mentioned parameter and regression equation, by test samples 2008-2011 crops temperature vegetation drought index number According to inputting in drought loss monitoring model respectively, each year disaster area estimated value of test samples is acquired, and according to drought loss mark Standard divides table and obtains each year drought loss monitoring result, as shown in the table:

6 drought loss monitoring result of table

In table 6, disaster-stricken rate be the ratio between disaster area and cultivated area, disaster area deviation ratio be disaster area estimated value and Actual value absolute deviation accounts for the ratio in practical disaster area.By table 6 as it can be seen that 2008-2011 as test samples, covers Non-irrigated four grades from no drought to weight, wherein disaster area deviation ratio increases with the decline of practical drought loss, and weight such as occurs 2009 of drought, practical disaster area is more than 1500khm², and disaster area deviation ratio is only 6%, and occur in no drought 2010, practical disaster area was less than 100khm², disaster area deviation is practical disaster area more than 4 times, from disaster area deviation Value is as can be seen that the above results are because the estimated value and actual value in each year disaster area have 200khm²The deviation of left and right, Deviation ratio is then increased with the reduction in practical disaster area；From the monitoring result of drought loss, the drought of test samples Monitoring grade and drought actual grade are consistent, have certain feasibility in the qualitative monitoring of disaster loss grade.

Present case describes the agricultural drought disaster grade monitoring technology process based on crops temperature vegetation drought index in detail, And by taking Henan Province as an example, 2001 to 2011 province's agricultural drought disaster related data is had collected, 2008-2011 annual data is made For test samples, each year disaster area is estimated using the optimized parameter and equation of linear regression of model, and to drought etc. Grade is monitored, and contrastive detection result and practical drought loss.From the monitoring result of drought loss, the drought of test samples Calamity monitoring grade and drought actual grade are consistent, have certain feasibility, Neng Gouman in the qualitative monitoring of disaster loss grade The basic demand of sufficient Droughts grade monitoring.Agricultural drought disaster grade monitoring method based on crops temperature vegetation drought index Data accessibility is strong, and method is simple, and operation is easy, and has certain operational use prospect.

The above content is a further detailed description of the present invention in conjunction with specific preferred embodiments, and it cannot be said that Specific implementation of the invention is only limited to these instructions.For those of ordinary skill in the art to which the present invention belongs, In Under the premise of not departing from present inventive concept, a number of simple deductions or replacements can also be made, all shall be regarded as belonging to of the invention Protection scope.

Claims

1. a kind of agricultural drought disaster grade monitoring method based on temperature vegetation drought index TVDI, which is characterized in that including as follows Step:

One, the preparation of data:

Agricultural drought disaster study on monitoring is carried out using the data that remote sensing technology obtains, wherein remotely-sensed data used is NASA, http: // Ladsweb.nascom.nasa.gov/ provide MODIS LST product, vegetation index product, land use covering product and Dem data, and arable land multiple crop index data；Statistical data is the cultivated area data in China Statistical Yearbook, disaster area Data and planting industry management department, People's Republic of China's agricultural rural area portion, what http://www.zzys.moa.gov.cn/ was provided Model on Sown Areas of Farm data and crops phenological calendar；

Two, the land surface temperature data reconstruction based on background value and undulating value:

MODIS LST product based on acquisition is reconstructed land surface temperature data LST according to corresponding method；

The land surface temperature data reconstruction based on background value and undulating value of the step 2 includes reconstructing method and reconstruct behaviour Make step, specific method and operating process are as follows:

(1) LST reconstructing method

Based in Cressman objective analysis method using initial fields value with correct the thought that value approaches observation jointly, by pixel The many years average level of LST, i.e. background value are regarded as the initial fields value of pixel；By the undulating value of pixel LST, i.e. the periphery radius of influence The observation of interior pixel and the difference of background value carry out interpolation pixel LST with this once to correct interpolation as value is corrected；Pixel a The specific algorithm of (i, j) interpolation can be expressed as follows:

LST_insert=LST_background+LST_varianceFormula 1-1

In formula, LST_insertFor pixel a LST interpolation result, LST_backgroundFor more years background values of pixel a；

Since there are missing values and low quality data in LST data, these points should be removed when calculating background value, it is therefore desirable to LST time series data is reconstructed；The less asymmetric Gaussian function fitting method of setup parameter needed for selecting is to the LST time Sequence data is reconstructed, and to the time series data LST after fitting_GaussianEach many years phase background value is sought, n is data set In a certain issue according to contained year, as shown in above formula 1-2；If certain pixels missing values or low quality number in time series data According to excessive, can not be fitted, which is then replaced by the average value of ground mulching type pixel similar in the radius of influence of periphery It is generation, identical at least over half ground mulching type in historical years, wherein existing for missing pixel and periphery pixel high The different problem of path difference, using the every relationship for increasing 1000m temperature decline about 6K of height above sea level, removal participates in calculating the elevation pair of pixel The influence of LST；

LST_varianceIt for the undulating value at pixel a, is determined by ground mulching type and the radius of influence, Δ LST is in the radius of influence The observation of high quality pixel and the difference of background value, K are the high quality pixel number of identical earth's surface cover type in the radius of influence, W_KFor its respective weights, it is calculated by following formula:

In formula, d_ijkFor the distance of the high quality pixel of interpolation pixel to identical ground surface type；R is the radius of influence, since temperature exists Variation within the scope of horizontal distance 6000m is generally less than 0.6K, therefore, will affect radius R and is set as 3, the i.e. height in 7x7 window Quality pixel participates in interpolation calculation；

(2) LST reconstructed operation step

1) data prediction；It chooses MODIS product and is used for draught monitor area LST data reconstruction；Used in reconstruct MOD11A2, MOD12Q1 product pass through MODIS re-projection tool MODISReprojectionTool, MRT re-projection first WGS84 coordinate system, and extract the data of LST round the clock in product, round the clock LST mass control file and LAI/fPAR system soil Utilize/cover classification result；The digital elevation data Digital ElevationModel provided by NASA is also provided in research, Dem data is first registrated by DEM with MODIS data, then by the spatial resolution of dem data from 30m resampling be 1000m；

Since the type and quantity of data used are more, the studies above data are cut out before treatment, consider restructing algorithm Involved in spatial window convolution algorithm, therefore select rectangular mask file data are cut, spatial resolution 1000m, Ensure to monitor region as far as possible and is in rectangular area center；The data used after above-mentioned pretreatment are as follows: round the clock LST data, The file of LST mass control round the clock, dem data and MOD12Q1 data, carry out Band fusion to above data, constitute LST number round the clock According to collection, the file data collection of LST mass control round the clock, ground mulching categorical data collection and dem data；

2) more years background values of LST are calculated；Asymmetric Gaussian function fitting will be carried out respectively by LST data set round the clock；Later, it examines quasi- Time series data collection after conjunction: since the pixel value of LST is open type temperature in data set, so successful picture will be fitted Time series all values where first are set as 0；Thereafter, pixel same period long-time average annual value is calculated to all pixels, as the phase Background value obtains LST background value data set, round the clock each 46 wave bands；The pixel for being 0 to pixel value in LST background value data set, Carry out interpolation in 7x7 window around the pixel: firstly, referring to ground mulching categorical data collection, choose in 7x7 window with it is to be inserted It is worth the identical pixel of pixel ground mulching type, recognizes with interpolation pixel type same number more than half in historical years For identical earth's surface cover type, if choosing whole pixels without if；Then, the weight meter in Cressman objective analysis method is utilized Calculation method calculates each pixel weight according to pixel is chosen at a distance from interpolation pixel；Later, it using dem data, calculates and chooses The depth displacement of pixel and interpolation pixel, using the every raising 1000m of height above sea level, temperature declines the relationship of 6K, by all selection pixels The unified pixel to interpolation of LST where elevation temperature；Finally, to LST and its weight after selection pixel " levelling " in window It is weighted summation, obtains the LST background value data set of space and time continuous；

3) file data collection is controlled using LST mass round the clock, member item by item is carried out to LST data round the clock respectively and is screened: according to round the clock LST mass controls file, only retains high quality pixel value in each wave band, remaining pixel value is set as 0, obtains interpolation round the clock LST data set；

4) using calculate gained LST background value data set, round the clock interpolation LST data set and ground mulching categorical data collection into Row interpolation, the specific steps are a. to utilize ground mulching categorical data collection, chooses in the pixel periphery 7x7 window that pixel value is 0, with The identical pixel of 0 value pixel ground mulching type, if choosing whole pixels without if；B. using in Cressman objective analysis method Weighing computation method calculates the weight for choosing pixel；C. the actual value for choosing pixel, i.e. high quality pixel observation and back are calculated The difference of scape value, as undulating value；D. summation is weighted to the corresponding weight of undulating value for choosing pixel in window, obtained most Whole LST interpolation result；

Three, crop planting area normalized differential vegetation index-land surface temperature feature space building:

Draught monitor region arable land is extracted, each growth period crop planting of crops area is constructed jointly using many years contemporaneous data and returns One changes vegetation index-land surface temperature NDVI-LST scatter plot, and is fitted each phase dry and wet side equation, constructs NDVI-LST feature Space；

Four, the calculating of crops temperature vegetation drought index:

It is based on step 3 as a result, proposing in nineteen ninety temperature vegetation drought index using Price calculates monitoring region farming Species growing area crops temperature vegetation drought index C-TVDI；

Five, the drought loss monitoring based on crops temperature vegetation drought index:

Based on the design of Supervised classification device, determine that monitoring region is based on crops temperature vegetation drought index by historical data Drought loss monitoring model parameter, carry out the monitoring of drought；

The step 5 is monitored based on the drought loss of crops temperature vegetation drought index, including based on Supervised classification device Design determines the ginseng of monitoring drought loss monitoring model of the region based on crops temperature vegetation drought index by historical data Number, carries out the monitoring of drought；It is specific as follows:

Based on the design of Supervised classification device, made with the remotely-sensed data of former years crop growth period and actual measurement of each year drought loss data The parameter in drought loss monitoring model is determined by the study to training sample for training sample；Newest one is recycled later The crop growth period remotely-sensed data in year is monitored the year-end drought loss in this year；The model mainly estimates letter by disaster area Number, parameter optimization and drought loss monitoring three parts composition:

(1) disaster area evaluation function: disaster area is defined as the Model on Sown Areas of Farm by various natural calamities, and same Year by several or natural calamity several times, wherein to endanger maximum primary calculating disaster area, do not repeat meter calamity, drought by Calamity area is then this year one or many maximum Model on Sown Areas of Farm influenced by drought；Consider that disaster area value is in year Drought considers monitoring section if it is the cropping pattern of multiple multiple cropping, and crops are through gathering in the biggest impact areas of crops The influence of a drought terminates later, therefore, disaster area evaluation function is divided into different planting seasons, is set as follows:

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

In formula, S_iFor 1 year disaster area estimated value；K~l is the annual first phase to gather in crop Critical growing period, and m~n is annual Final Issue gathers in crop Critical growing period, s_ij(j=k~l) is 1 year first phase to gather in suffering from drought for each Critical growing period of crop Area, s_ij(j=m~n) is the area suffered from drought that Final Issue gathers in crop each growth period；Since the generation from suffering from drought drought needs The accumulation of certain time is wanted, the damage caused by a drought of single issue evidence may not cause drought, and Spatial Variability arid in a short time is little, therefore Using the average value of each phase area suffered from drought as disaster area；Annual disaster area value is then the maximum value 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 the growth Phase threshold parameter t_j(j=k~l or m~n) is determined；By taking crops temperature vegetation drought index as an example, index value is closer to 1 Indicate more arid, then the area suffered from drought S of 1 year j-th Critical growing period_ijFor the year issue, C-TVDI value is greater than the phase in Threshold value t_jPixel occupied area, threshold value t_jFor model parameter, parameter and growth period used are corresponded；

(2) parameter optimization seeks the process being most worth using objective function of the optimizing algorithm to setting；Drought loss monitoring model The parameter optimization of model is carried out using many years remotely-sensed data and agricultural disaster area data as training sample；Model parameter it is optimal Solution is regarded as under this group of parameter, and practical disaster area in many years disaster area estimated value and years of training sample is calculated Minimum average B configuration deviation, therefore objective function is set are as follows:

In formula, S_iFor 1 year disaster area estimated value, S_i0For 1 year practical disaster area value, n was training sample number, to mesh Scalar functions acquire corresponding model parameter t when minimum value_j(j=k~l or m~n) is optimized parameter；

Grid-search algorithms are selected to carry out parameter optimization, the speed of searching optimization of grid-search algorithms is very fast, obtains globally optimal solution, no Locally optimal solution can be fallen into；But setup parameter range is wanted before search, step-size in search lesser search larger in parameter area In, time consumption for training is longer；

Parameter optimization process is expressed as follows: the remotely-sensed data in training sample is carried out the extraction of crops Indices by (i), and is selected It selects and participates in the crop growth period that model calculates；(ii) minimum target functional value, each growth period parameters t are set_i0Search range, search Suo Buchang and parameter t_i0Initial value；(iii) the agricultural remote sensing exponent data of extraction is input in the evaluation function of disaster area, Obtain disaster area estimated value；(iv) disaster area data in disaster area estimated value and training sample are brought into objective function In, calculating target function value retains parameter current, and minimum value is updated if the value is less than current minimum target functional value； (v) new parameter value is obtained by grid-search algorithms, and repeats above-mentioned (iii), (iv) step, if whole parameters have traversed Optimizing is then completed, optimized parameter is recorded；

(3) drought loss monitors, i.e., by year crop growth phase remotely-sensed data to be monitored, supervises to this year drought loss It surveys；Since the sample number for participating in training is less, training sample is limited to the representativeness of different degrees of the condition of a disaster, therefore utilizes training sample Linear relationship between this disaster area estimated value and practical disaster area adjusts the disaster area estimated value in year to be monitored It is whole；Specific step is as follows for drought loss monitoring: (i) calculates practical disaster area and disaster-stricken face under optimized parameter in training sample The equation of linear regression of product estimated value；(ii) year crops Critical growing period Indices to be monitored are extracted；(iii) using disaster-stricken Area reckoning function seeks the disaster area estimated value under optimized parameter；(iv) certainly using estimated value in previous step as linear equation Variable, disaster area estimated value after being adjusted；(v) disaster-stricken rate I is calculated using disaster area estimated value after adjustment, and referring to drought Calamity classification standard divides table, obtains drought loss monitoring result；Drought loss standard divides as follows:

Without drought: I≤5%；Mild drought: 5% < I≤10%；Mild drought: 10% < I≤20%；Severe drought: 20% < I ≤ 30%；Serious Drought Event: I >=30%.

2. a kind of agricultural drought disaster grade monitoring method based on temperature vegetation drought index TVDI according to claim 1, It is characterized in that, the building of step 3 crop planting area normalized differential vegetation index-land surface temperature feature space, packet It includes using draught monitor region arable land is extracted, constructs each growth period crop planting of crops area jointly using many years contemporaneous data Normalized differential vegetation index-land surface temperature NDVI-LST scatter plot, and it is fitted each phase dry and wet side equation, building NDVI-LST is special Levy space；It is specific as follows:

(1) using the sintetics of MODIS vegetation index 16 days MOD13A2, it is reconstructed after daytime LST 8 day data and monitoring Region crops multiple crop index data, each data spatial resolution is 1000m；

1) MODIS vegetation index product MOD13A2 is subjected to projection transform with MODIS re-projection tool, and extracted therein NDVI data and quality assessment data QualityAssurance, QA；

2) to monitoring section, data carry out arable land image element extraction year by year, wherein due in monitoring region crops multiple crop index data In, the pixel value that cultivated area data are the most identical in pixel area and statistical yearbook is selected, each year Area distortion control exists 10% hereinafter, pixel of being ploughed using the pixel of multiple crop index pixel value as monitoring region, extracts monitoring region arable land；

3) referring to 16 grades of classification standards of NDVI mass in QA data, the screening NDVI quality of data is the medium above pixel, i.e., In QA data 2~5 pixels of the value less than 4, remaining pixel in NDVI data is assigned to 0；

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

8 day data of LST on daytime is merged into 16 days LST data, the high quality pixel of original LST data substitutes interpolation in merging Pixel merges algorithm and refers to MOD11A2 user's manual, to two pixels when two pixels are all high quality pixel or interpolation pixel It is worth averaged, as merging data pixel value；

(2) firstly, be grouped to many years data on schedule, NDVI, LST data is divided into several groups, include that identical quantity is same in each group Phase NDVI, LST data construct NDVI-LST feature space to all crops pixels of each group of data；Then, with each group of data In maximum, minimum value the difference of NDVI value 1 percent be step-length, it is gradually long to obtain LST in group corresponding to different NDVI values Maximum, minimum value；Dry and wet side equation finally maximum using NDVI and LST, in minimum value fitting each group of data.

3. a kind of agricultural drought disaster grade monitoring method based on temperature vegetation drought index TVDI according to claim 1, It is characterized in that, the calculating of the step 4 crops temperature vegetation drought index, including proposed based on Price in nineteen ninety Temperature vegetation drought index calculates monitoring region C-TVDI value；The calculation method of specific temperature vegetation drought index indicates such as Under:

Wherein,

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

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

In formula, T_sFor the LST value of pixel (i, j), T_{s_max}With T_{s_min}The corresponding dry side of the NDVI value of respectively pixel (i, j) and wet LST value on side, dry side carry out linear regression by the scatter plot constituted to the corresponding maximum value of NDVI values all kinds of in feature space It acquires, wet side is then acquired by the linear regression of corresponding minimum value scatter plot；a₁, a₂, b₁, b₂Respectively in the equation of dry and wet side Undetermined coefficient.