CN102176002A - Surface water heat flux remote sensing inversion-based drought monitoring method and system - Google Patents

Abstract

The invention discloses drought monitoring methods and system based on earth's surface water and heat remote-sensing inversion. The method, including the following steps: calculate surface net radiation flux

Based on Remote sensing parameters inverting, soil heat flux is calculated

Using iterative calculation method, sensible heat flux is calculated

According to the calculating of the surface flux, zoning water shortage index

According to region water shortage index RWSI, the monitoring of regional drought degree is carried out.

Description

Drought monitoring method and system based on surface water thermoflux remote-sensing inversion

Technical field

The present invention relates to the monitoring technical field of the ecosystem, particularly relate to drought monitoring method and system based on surface water thermoflux remote-sensing inversion.

Background technology

China's northern natural ecosystem Agro-ecological System of unifying is coerced by moisture condition generally, by means of quick, macroscopical observation advantage of remote sensing means, the ecosystem surface water thermoflux of simulation, monitored area yardstick is the problem in science with theory significance and actual value.Remote sensing observations is unique estimation means that regional surface water thermoflux may be provided.

Evapotranspire and comprised the transpiration of moisture in evaporation from topsoil and the plant, it belongs to the part of surface water thermoflux, being the important component part of keeping land face water balance, also is the major part of keeping face of land energy equilibrium, is the important step in the surface water cyclic process.Measure exactly and estimate that evapotranspiration not only has significance to research global climate differentiation, ecological environment problem and water resources assessment etc., and of far-reaching significance to draining and the irrigation of instructing agricultural, the utilization factor of monitoring agriculture damage caused by a drought, raising agricultural water resources etc.

Numerous to the achievement in research of evapotranspiring in the face of land, many scholars have extracted different evaluation methods and model from different perspectives.The difference that trial such as Seguin and Itier and Hatfield utilizes the satellite thermal infrared data to calculate terrestrial radiation temperature and temperature is estimated the evapotranspiration in large scale zone; Meneti etc. also propose to calculate face of land evapotranspiration based on face of land energy equilibrium index SEBI; Roerink etc. have proposed the face of land energy equilibrium index of simplifying, and (Simplified Surface Energy Balance Index S-SEBI), calculates evaporite ratio earlier, and further estimation is evapotranspired then; Proposition face of land energy equilibrium systems (SEBS) such as Su evapotranspire in order to the estimation face of land, this method is by determining a series of faces of land physical parameter and setting up heat conduction roughness model, utilize similarity theory to determine face of land dynamic characteristic, calculate evaporite ratio based on face of land energy equilibrium index (SEBI) then.Systems such as Bastiaanssen have delivered and have utilized the serial paper of Landsat TM data based on the face of land, the SEBAL model assessment zone evapotranspiration of face of land energy-balance equation, and are widely used in the countries and regions under the different weather conditions such as the U.S., Spain, Pakistan, Egypt, Sri Lanka, Turkey and China; Kustas et al etc. has also carried out utilizing the remote sensing Application and Development research of model of evapotranspiring.

At present, soil moisture and the damage caused by a drought of utilizing the thermal infrared technology to carry out the zone are monitored mainly according to soil water balance equation and principle of energy balance, relation estimation soil moisture by soil surface emissivity and surface temperature mainly comprises thermal inertia method, moisture index method and lack of water index method etc.

The thermal inertia method is that Wang et al has developed index of water deficit (WDI) index draught monitor and grade classification have been carried out in the zone, loess plateau with the method for remote sensing techniques monitoring soil moisture.The moisture index method is to utilize the strong absorption characteristic of moisture at short-wave infrared, and the formation subindex carries out the Regional Drought monitoring.Mcffters et al. is combined into normalization difference moisture index (NDWI, Normalized Difference Water Index) by TM green glow and near-infrared band the earliest; Kogan proposition vegetation state indices (Vegetation Condition Index, VCI).The lack of water index method is people the measure index into crop water shortage of likening to of vegetation actual evapotranspiration and potential evapotranspiration amount.Jackson et al proposed the crop water shortage index (Crop Water StressIndex, notion CWSI), Moran et al. proposed the water deficit index (Water Deficit Index, WDI).

According to the classification of Courault et al. (2005), the model that evapotranspires can divide four classes: the empirical model method, energy residue modelling is determined model method, the vegetation index modelling.Wherein energy residue class model is to utilize remotely-sensed data to carry out regional surface flux and evapotranspire estimating the superior in operability model I at present.But there is following deficiency in above-mentioned energy residue class model:

1. only can and evapotranspire and estimate, not be suitable for complicated topographical features even, smooth regional surface flux;

2. the artificial affirmation of doing wet pixel need be carried out, the automatic examination of doing wet pixel of simulated domain can't be realized;

3. only can carry out the face of land energy of short time and evapotranspire the estimation simulation the zonule.

Summary of the invention

The object of the present invention is to provide drought monitoring method and system based on surface water thermoflux remote-sensing inversion.It can estimate that the face of land energy in the zone under different weathers, the topographic condition and ET distribute in conjunction with conventional ground meteorological data, for the monitoring of regional agriculture damage caused by a drought provides the practical technique support.

Be the drought monitoring method based on surface water thermoflux remote-sensing inversion of realizing that purpose of the present invention provides, described method comprises the following steps:

Step 100. is calculated surface net radiation flux

Wherein, R _{S ↓}Be the sun shortwave radiation of incident, α is a surface albedo, ε _sBe face of land emissivity; σ is a Stefan-Boltzmann constant (5.6696 * 10 ^-8Wm ^-2K ^-4), T _sIt is surface temperature; T _aBe the air themperature at reference altitude place, ε _aIt is the air ratio radiance;

Step 200. is calculated soil heat flux based on the remote sensing parametric inversion:

$G=\left[\frac{(T_s-273.15)}{\mathrm{\α}}(0.0038\mathrm{\α}+0.0074{\mathrm{\α}}^2)(1-0.98\mathrm{NDV}{I}^4)\right]R_n$

Wherein, α is a surface albedo, T _sBe surface temperature, NDVI is the water body of vegetation normalization index for study area, because of its between limpid deep water body and wetland, get G _Water=0.5R _n

Step 300. is utilized iterative calculation method, calculates sensible heat flux:

$H={\mathrm{\ρ}}_a{c}_p\frac{(T_s-T_a)}{r_{\mathrm{ah}}}={\mathrm{\ρ}}_a{c}_p\frac{\mathrm{dT}}{r_{\mathrm{ah}}}$

Wherein, ρ _aBe atmospheric density (kgm ^-3); c _pBe the specific heat capacity at constant pressure (1004.07Jkg of air ^-1K ^-1); T _sBe the face of land of remote sensing or the surface temperature of canopy; T _aAir themperature (K) for the reference altitude place; DT is two height Z ₂With Z ₁Between the temperature difference; r _AhBe two aerodynamic resistance (ms between height ^-1);

Step 400. is according to the calculating of described surface flux, zoning lack of water index

According to regional lack of water index RWSI, carry out the monitoring of regional degree of drought.

Described step 100 comprises the following steps: step 110. calculating surface albedo α;

MODIS The data following formula calculates surface albedo:

α ^MODIS＝0.160α ₁+0.291α ₂+0.243α ₃+0.116α ₄

+0.112α ₅+0.081α ₇-0.0015

In the formula, α _i, i=1 ... each narrow wave band surface albedo on the corresponding satellite sensor of n difference;

Step 120. is calculated face of land emissivity ε _s

Extract water body earlier, compose typical emissivity value 0.995 with water body; The emissivity estimation equation of remaining mixed pixel is as follows:

${\mathrm{\&epsiv;}}_s=\left\{\begin{array}{c}0.9624744+0.0652704P_v-0.0461286P_v^2,P_v<0.5\\ 0.9643744+0.0614704P_v-0.0461286P_v^2,P_v=0.5\\ 0.9662744+0.0576704P_v-{0.0461286P}_v^2,P_v>0.5\end{array}\right.$

Vegetation coverage is provided by following formula:

$P_v=\frac{(\mathrm{NDVI}-{\mathrm{NDVI}}_{\min })}{({\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min })}$

In the formula, NDVI _Min, NDVI _Max, NDVI corresponds to pure exposed soil pixel respectively, the NDVI value of pure vegetation pixel and mixed pixel correspondence is got NDVI _Min=0.099, NDVI _Max=0.77; When NDVI smaller or equal to 0.099 the time, vegetation coverage P _v=0; When NDVI more than or equal to 0.77 the time, vegetation coverage P _v=1;

Step 130. is calculated the sun shortwave radiation R of incident _{S ↓}

The sun shortwave radiation in a certain moment can be regarded direct solar radiation, scattered radiation and reflected radiation three parts as and form:

R _s↓＝R _s+R _d+R _r

In the formula, R _sBe direct radiation; R _dBe scattered radiation; R _rBe reflected radiation.Wherein:

$R_s={\left(\frac{1}{\mathrm{\&rho;}}\right)}^2{\mathrm{<figure-callout\; id="0"\; label="I"\; filenames="HSA00000414253100042.png"\; state="\{\{state\}\}">I</figure-callout>}}_0\cos {\mathrm{\&theta;\&tau;}}_b={\mathrm{<figure-callout\; id="0"\; label="R"\; filenames="HSA00000414253100042.png"\; state="\{\{state\}\}">R</figure-callout>}}_0{\mathrm{\&tau;}}_b\cos \mathrm{\&theta;}$

Wherein, I ₀Be solar constant, get 1367 (Wm ^-2); R ₀Solar radiation for the straight incident of atmospheric envelope roof pendant; θ is the solar zenith angle on domatic; τ _bBe atmosphere direct projection transmissivity; (1/ ρ) ²Be called solar distance correction coefficient or Earth's orbit and correct the factor, about 10 ^-4During precision, calculate by fourier progression expanding method:

R _d＝R ₀τ _d?cos ²β/2sina

In the formula, β is the gradient, and a is a sun altitude.

R _r＝rR ₀τ _rsin ²β/2sina

In the formula, r is the first reflectivity in slope, the face of land, can simplify being taken as 0.2; τ _rTransmissivity for reflected radiation.

Described step 300 comprises the following steps:

Step 310. provides degree of stability correction function Ψ based on similarity theory _m(x), Ψ _h(x), respectively to friction velocity u _*, aerodynamic resistance r _AhRevise:

$u_{*}=\frac{\mathrm{ku}\left(z\right)}{\ln \left(\frac{z-d}{z_{\mathrm{om}}}\right)-{\mathrm{\&Psi;}}_m\left(\frac{z-d}{L}\right)+{\mathrm{\&Psi;}}_m\left(\frac{z_{\mathrm{om}}}{L}\right)}$

$r_{\mathrm{ah}}=\frac{1}{{\mathrm{ku}}_{*}}\left[\ln \left(\frac{Z_2}{Z_1}\right)-{\mathrm{\&Psi;}}_h\left(\frac{Z_2}{L}\right)+{\mathrm{\&Psi;}}_h\left(\frac{Z_1}{L}\right)\right]$

Wherein, k is von Karman constant (getting 0.41, dimensionless), and z is the height of mean wind speed correspondence, z _OmBe face of land power roughness, Z ₂It is the reference altitude of measuring temperature; Z ₁=z _Oh, the heat delivered roughness; u _*Be friction velocity (ms ^-1), L is a Monin-Obukhov length,

In the formula, ρ _aBe atmospheric density; c _pIt is the specific heat capacity at constant pressure of air; T _sBe the face of land of remote sensing or the surface temperature of canopy; K is a von Karman constant; u _*It is friction velocity; G is an acceleration of gravity; H is the sensible heat flux;

Step 320. is calculated based on the surface temperature under the DEM, T _{S_dem}=T _s+ 0.0065 (h-h _Mean), wherein, h is a sea level elevation, h _MeanSea level on the average for survey region;

Step 330. is chosen extremely dried pixel and extremely wet pixel, with described T _{S_dem}Substitution formula dT=aT _{S_dem}+ b, wherein,

T _HotThe temperature of the extremely dried pixel of expression, T _ColdThe temperature of the extremely wet pixel of expression.

In the step 330, the described criterion of choosing extremely dried pixel and extremely wet pixel is: the MSAVI about 0.1 reaches wherein combination correspondence " Hot " point of maximum LST; 0.8 about MSAVI and combination correspondence " Cold " point of wherein minimum LST.

Described step 400 comprises the following steps:

Based on the Regional Drought index of RWSI, utilize the national standard of the damage caused by a drought grade classification of soil moisture, the arid ranking score grade standard of the RWSI of each time period of converting carries out regional arid remote sensing monitoring.

For realizing that purpose of the present invention also provides the draught monitor system based on surface water thermoflux remote-sensing inversion, described system comprises:

The surface net radiation flux computing module is used to calculate surface net radiation flux $R_n=(1-\mathrm{\&alpha;})R_{s\&DownArrow;}+{\mathrm{\&epsiv;}}_s\mathrm{\&sigma;}({\mathrm{\&epsiv;}}_a{T}_a^4-T_s^4);$

Wherein, R _{S ↓}Be the sun shortwave radiation of incident, α is a surface albedo, ε _sBe face of land emissivity; σ is a Stefan-Boltzmann constant (5.6696 * 10 ^-8Wm ^-2K ^-4), T _sIt is surface temperature; T _aBe the air themperature at reference altitude place, ε _aIt is the air ratio radiance;

The soil heat flux computing module is used for based on the remote sensing parametric inversion, calculates soil heat flux:

$G=\left[\frac{(T_s-273.15)}{\mathrm{\&alpha;}}(0.0038\mathrm{\&alpha;}+0.0074{\mathrm{\&alpha;}}^2)(1-0.98{\mathrm{NDVI}}^4)\right]R_n$

Wherein, α is a surface albedo, T _sBe surface temperature, NDVI is the water body of vegetation normalization index for study area, because of its between limpid deep water body and wetland, get G _Water=0.5R _n

The sensible heat flux computing module; Be used to utilize iterative calculation method, calculate sensible heat flux:

$H={\mathrm{\&rho;}}_a{c}_p\frac{(T_s-T_a)}{r_{\mathrm{ah}}}={\mathrm{\&rho;}}_a{c}_p\frac{\mathrm{dT}}{r_{\mathrm{ah}}}$

Wherein, ρ _aBe atmospheric density (kgm ^-3); c _pBe the specific heat capacity at constant pressure (1004.07Jkg of air ^-1K ^-1); T _sBe the face of land of remote sensing or the surface temperature of canopy; T _aAir themperature (K) for the reference altitude place; DT is two height Z ₂With Z ₁Between the temperature difference; r _AhBe two aerodynamic resistance (ms between height ^-1);

Zone lack of water Index for Calculation module is used for the calculating according to described surface flux, zoning lack of water index According to regional lack of water index RWSI, carry out the monitoring of regional degree of drought.

Described surface net radiation flux computing module comprises:

The surface albedo computing module is used to calculate surface albedo α;

MODIS The data following formula calculates surface albedo:

α ^MODIS＝0.160α ₁+0.291α ₂+0.243α ₃+0.116α ₄

+0.112α ₅+0.081α ₇-0.0015

In the formula, α _i, i=1 ... each narrow wave band surface albedo on the corresponding satellite sensor of n difference;

Face of land emissivity computing module is used to calculate face of land emissivity ε _s

Extract water body earlier, compose typical emissivity value 0.995 with water body; The emissivity estimation equation of remaining mixed pixel is as follows:

${\mathrm{\&epsiv;}}_s=\left\{\begin{array}{c}0.9624744+0.0652704P_v-0.0461286{\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v<0.5\\ 0.9643744+0.0614704P_v-0.0461286{\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v=0.5\\ 0.9662744+0.0576704P_v-{0.0461286\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v>0.5\end{array}\right.$

Vegetation coverage is provided by following formula:

$P_v=\frac{(\mathrm{NDVI}-{\mathrm{NDVI}}_{\min })}{({\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min })}$

In the formula, NDVI _Min, NDVI _Max, NDVI corresponds to pure exposed soil pixel respectively, the NDVI value of pure vegetation pixel and mixed pixel correspondence is got NDVI _Min=0.099, NDVI _Max=0.77; When NDVI smaller or equal to 0.099 the time, vegetation coverage P _v=0; When NDVI more than or equal to 0.77 the time, vegetation coverage P _v=1;

Sun shortwave radiation computing module is used to calculate the sun shortwave radiation R of incident _{S ↓}

The sun shortwave radiation in a certain moment can be regarded direct solar radiation, scattered radiation and reflected radiation three parts as and form:

R _s↓＝R _s+R _d+R _r

In the formula, R _sBe direct radiation; R _dBe scattered radiation; R _rBe reflected radiation.Wherein:

$R_s={\left(\frac{1}{\mathrm{\&rho;}}\right)}^2{I}_0\cos {\mathrm{\&theta;\&tau;}}_b={\mathrm{<figure-callout\; id="0"\; label="R"\; filenames="HSA00000414253100042.png"\; state="\{\{state\}\}">R</figure-callout>}}_0{\mathrm{\&tau;}}_b\cos \mathrm{\&theta;}$

Wherein, I ₀Be solar constant, get 1367 (Wm ^-2); R ₀Solar radiation for the straight incident of atmospheric envelope roof pendant; θ is the solar zenith angle on domatic; τ _bBe atmosphere direct projection transmissivity; (1/ ρ) ²Be called solar distance correction coefficient or Earth's orbit and correct the factor, about 10 ^-4During precision, calculate by fourier progression expanding method:

R _d＝R ₀τ _dcos ²β/2sina

In the formula, β is the gradient, and a is a sun altitude.

R _r＝rR ₀τ _rsin ²β/2sina

In the formula, r is the first reflectivity in slope, the face of land, can simplify being taken as 0.2; τ _rTransmissivity for reflected radiation.

Described sensible heat flux computing module comprises:

Correcting module is used for providing degree of stability correction function Ψ based on similarity theory _m(x), Ψ _h(x), respectively to friction velocity u _*, aerodynamic resistance r _AhRevise:

$u_{*}=\frac{\mathrm{ku}\left(z\right)}{\ln \left(\frac{z-d}{z_{\mathrm{om}}}\right)-{\mathrm{\&Psi;}}_m\left(\frac{z-d}{L}\right)+{\mathrm{\&Psi;}}_m\left(\frac{z_{\mathrm{om}}}{L}\right)}$

$r_{\mathrm{ah}}=\frac{1}{{\mathrm{ku}}_{*}}\left[\ln \left(\frac{Z_2}{Z_1}\right)-{\mathrm{\&Psi;}}_h\left(\frac{Z_2}{L}\right)+{\mathrm{\&Psi;}}_h\left(\frac{Z_1}{L}\right)\right]$

Wherein, k is von Karman constant (getting 0.41, dimensionless), and z is the height of mean wind speed correspondence, z _OmBe face of land power roughness, Z ₂It is the reference altitude of measuring temperature; Z ₁=z _Oh, the heat delivered roughness; u _*Be friction velocity (ms ^-1), L is a Monin-Obukhov length,

In the formula, ρ _aBe atmospheric density; c _pIt is the specific heat capacity at constant pressure of air; T _sBe the face of land of remote sensing or the surface temperature of canopy; K is a von Karman constant; u _*It is friction velocity; G is an acceleration of gravity; H is the sensible heat flux;

The surface temperature computing module is used to calculate the surface temperature based under the DEM,

T _{S_dem}=T _s+ 0.0065 (h-h _Mean), wherein, h is a sea level elevation, h _MeanSea level on the average for survey region;

Pixel is chosen module, is used to choose extremely dried pixel and extremely wet pixel, with described T _{S_dem}Substitution formula dT=aT _{S_dem}+ b, wherein,

T _HotThe temperature of the extremely dried pixel of expression, T _ColdThe temperature of the extremely wet pixel of expression.

Pixel is chosen in the module, and the described criterion of choosing extremely dried pixel and extremely wet pixel is: the MSAVI about 0.1 reaches wherein combination correspondence " Hot " point of maximum LST; 0.8 about MSAVI and combination correspondence " Cold " point of wherein minimum LST.

Described regional lack of water Index for Calculation module based on the Regional Drought index of RWSI, is utilized the national standard of the damage caused by a drought grade classification of soil moisture, and the arid ranking score grade standard of the RWSI of each time period of converting carries out regional arid remote sensing monitoring.

The invention has the beneficial effects as follows:

Drought monitoring method and system based on surface water thermoflux remote-sensing inversion of the present invention, be based on principle of energy balance and boundary layer similarity theory, set up regional surface water thermoflux remote sensing algorithm, utilize the MODIS satellite remote sensing date, introduce the influence that the landform and the face of land cover, the regional face of land parameter (surface temperature that inverting is relevant, surface albedo etc.) and vegetation structure parameter (vegetation coverage, normalized differential vegetation index NDVI etc.), in conjunction with conventional ground measured data and weather data, the face of land, the zoning energy and the parameter of evapotranspiring are to carry out draught monitor, for the monitoring of regional agriculture damage caused by a drought provides the practical technique support.

Description of drawings

Fig. 1 is the flow chart of steps of the drought monitoring method based on surface water thermoflux remote-sensing inversion of the present invention;

Fig. 2 is a flow chart of steps of calculating the sensible heat flux among the present invention;

Fig. 3 is meant by the VITT space synoptic diagram of index and surface temperature formation;

Fig. 4 is that the draught monitor system that the present invention is based on surface water thermoflux remote-sensing inversion gets structural representation;

Fig. 5 is the graph of a relation of considering (a) and not considering elevation and day evapotranspiration under (b) action of topography;

Fig. 6 is the influence curve figure of aspect to face of land water and heat.

Embodiment

In order to make purpose of the present invention, technical scheme and advantage clearer,, drought monitoring method and the system based on surface water thermoflux remote-sensing inversion of the present invention is further elaborated below in conjunction with drawings and Examples.Should be appreciated that specific embodiment described herein only in order to explanation the present invention, and be not used in qualification the present invention.

Drought monitoring method and system based on surface water thermoflux remote-sensing inversion of the present invention, be based on energy residual error principle, with reference to the related algorithm in SEBAL, SEBS and the S-SEBI model, utilize the spatial analysis functions of GIS, the remote sensing algorithm evapotranspires in the zone based on digital elevation model (DEM) parametric calibration of foundation.Utilize DEM (gradient, aspect and sea level elevation) that the model correlation parameter is carried out topographic correction, determine to do, wet pixel automatically, and pass through the dynamic surface effective height of estimation and then calculate roughness of ground surface.Can estimate that the face of land energy in the zone under different weathers, the topographic condition and ET distribute in conjunction with conventional ground meteorological data, expanded the regional scope of application, not only can carry out the face of land energy in zone, Plain (evenly topographical features) and the calculating of evapotranspiring, can also carry out the energy of knob, mountain area (complicated earth surface feature) and the calculating of evapotranspiring, for the monitoring of regional agriculture damage caused by a drought provides the practical technique support.

Fig. 1 is the flow chart of steps of the drought monitoring method based on surface water thermoflux remote-sensing inversion of the present invention, and as shown in Figure 1, the foundation basis of the drought monitoring method based on surface water thermoflux remote-sensing inversion of the present invention is the surface radiation energy-balance equation:

R _n＝H+G+LE+PH (1)

In the formula, R _nBe net radiation flux, G is a soil heat flux, and H is the sensible heat flux, and LE is a latent heat flux, and PH represents the energy that crop photosynthesis effect and biomass increase, because its value is too little, calculating is ignored often usually.As long as therefore know R _n, G, H value, just can obtain the LE value.

1. surface net radiation flux R _n, the clean revenue and expenditure of expression face of land shortwave and long-wave radiation is the main energy sources that energy between the face of land, atmosphere, momentum, moisture and other material molecules exchange.The net radiation that the unit surface area receives can be expressed as:

$R_n=(1-\mathrm{\&alpha;})R_{s\&DownArrow;}+{\mathrm{\&epsiv;}}_s\mathrm{\&sigma;}({\mathrm{\&epsiv;}}_a{T}_a^4-T_s^4)---\left(2\right)$

In the formula, R _{S ↓}Sun shortwave radiation for incident is also named total solar radiation; α is a surface albedo.ε _sBe face of land emissivity; σ is a Stefan-Boltzmann constant (5.6696 * 10 ^-8Wm ^-2K ^-4); T _sBe surface temperature (K); T _aBe reference altitude (Z ₂) air themperature located; ε _aIt is the air ratio radiance.

Wherein each Parameters Calculation is as follows in (2):

A) surface albedo α calculates

Liang etc. have drawn the experimental formula of the broadband ground albedo of multiple sensors based on the modeling effort of radiation delivery, and obtain effect (Liang, 2001 preferably in checking; Liang, et al., 2003).

The formula that this paper adopts Liang etc. to provide calculates surface albedo, and MODIS The data following formula calculates surface albedo:

α ^MODIS＝0.160α ₁+0.291α ₂+0.243α ₃+0.116α ₄ (2a)

+0.112α ₅+0.081α ₇-0.0015

In the formula, α _i, i=1 ... each narrow wave band surface albedo on the corresponding satellite sensor of n difference.

B) face of land emissivity ε _sCalculate

With reference to research (Qin Zhihao, et al., 2004; Qin Zhihao, et al., 2005), by determining the emissivity of typical feature, determine that with NDVI the composition of atural object (is vegetation coverage P then earlier _v), thereby determine the face of land emissivity that study area is interior.Specific practice is to extract water body earlier, composes the typical emissivity value 0.995 with water body; The emissivity estimation equation of remaining mixed pixel is as follows:

${\mathrm{\&epsiv;}}_s=\left\{\begin{array}{c}0.9624744+0.0652704P_v-0.0461286{\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v<0.5\\ 0.9643744+0.0614704P_v-0.0461286{\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v=0.5\\ 0.9662744+0.0576704P_v-{0.0461286\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v>0.5\end{array}\right.---\left(2b\right)$

Vegetation coverage provides (Zhang Renhua, 1996) by following formula:

$P_v=\frac{(\mathrm{NDVI}-{\mathrm{NDVI}}_{\min })}{({\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min })}---\left(2c\right)$

In the formula, NDVI _Min, NDVI _Max, NDVI corresponds to pure exposed soil pixel respectively, the NDVI value of pure vegetation pixel and mixed pixel correspondence.Get NDVI with reference to (Zhang Renhua, et al., 2004) _Min=0.099, NDVI _Max=0.77.When NDVI smaller or equal to 0.099 the time, vegetation coverage P _v=0; When NDVI more than or equal to 0.77 the time, vegetation coverage P _v=1.

C) sun shortwave radiation R of incident _{S ↓}Calculate

The sun shortwave radiation in a certain moment can be regarded direct solar radiation, scattered radiation and reflected radiation three parts as and form:

R _s↓＝R _s+R _d+R _r (2d)

In the formula, R _sBe direct radiation; R _dBe scattered radiation; R _rBe reflected radiation.

Consider under the atmospheric attenuation situation that arriving the arbitrary instantaneous direct radiation on domatic of earth surface can be expressed as:

$R_s={\left(\frac{1}{\mathrm{\&rho;}}\right)}^2{I}_0\cos {\mathrm{\&theta;\&tau;}}_b={\mathrm{<figure-callout\; id="0"\; label="R"\; filenames="HSA00000414253100042.png"\; state="\{\{state\}\}">R</figure-callout>}}_0{\mathrm{\&tau;}}_b\cos \mathrm{\&theta;}---\left(2e\right)$

Wherein, I ₀Be solar constant, get 1367 (Wm ^-2); R ₀Solar radiation for the straight incident of atmospheric envelope roof pendant; θ is the solar zenith angle on domatic; τ _bBe atmosphere direct projection transmissivity; (1/ ρ) ²Be called solar distance correction coefficient or Earth's orbit and correct the factor, about 10 ^-4During precision, can calculate (Liou, 2004) by fourier progression expanding method:

The scattering of sky is a homogeneous scattering under the ceiling unlimited condition, has linear relationship between direct radiation and scattered radiation, the computing formula of scattered radiation:

R _d＝R ₀τ _dcos ²β/2sina (2f)

In the formula, β is the gradient, and a is a sun altitude.

The calculating of reflected radiation is relevant with the reflectivity on the aspect on the face of land and the gradient and the face of land.Suppose that the reflection face of land is lambert's body, the computing formula that draws reflected radiation is:

R _r＝rR ₀τ _rsin ²β/2sina (2g)

In the formula, r is the first reflectivity in slope, the face of land, can simplify being taken as 0.2; τ _rTransmissivity for reflected radiation.

2. soil heat flux G is meant the heat-exchange power that enters soil inside, and is relevant with water, the heat interchange of soil inside.Bastiaanssen etc. think that light is unfavorable with the ratio relation estimation area thermoflux of vegetation characteristics, therefore propose a parametrization formula of taking all factors into consideration surface albedo, vegetation index and surface temperature:

$G=\left[\frac{(T_s-273.15)}{\mathrm{\&alpha;}}(0.0038\mathrm{\&alpha;}+0.0074{\mathrm{\&alpha;}}^2)(1-0.98{\mathrm{NDVI}}^4)\right]R_n---\left(3\right)$

In the formula, α is a surface albedo.T _sBe surface temperature (K); NDVI is the water body of vegetation normalization index for study area, because of its between limpid deep water body and wetland, get G _Water=0.5R _n

3. sensible heat flux (Sensible Heat Flux) H is the heat interchange that characterizes turbulence form between underlying surface and atmosphere, and it is the effect owing to the temperature difference, and energy is passed to atmosphere in the mode of turbulent flow, promptly is used to add that part of energy of hot-air.Fig. 2 is a flow chart of steps of calculating the sensible heat flux among the present invention, and as shown in Figure 2, sensible heat flux (H) with the expression formula of turbulence form is:

$H={\mathrm{\&rho;}}_a{c}_p\frac{(T_s-T_a)}{r_{\mathrm{ah}}}={\mathrm{\&rho;}}_a{c}_p\frac{\mathrm{dT}}{r_{\mathrm{ah}}}---\left(4\right)$

In the formula, ρ _aBe atmospheric density (kgm ^-3); c _pBe the specific heat capacity at constant pressure (1004.07Jkg of air ^-1K ^-1); T _sBe the face of land of remote sensing or the surface temperature of canopy (K); T _aBe (the Z of reference altitude place ₂) air themperature (K); DT is two height Z ₂With Z ₁Between the temperature difference (K); r _AhBe two aerodynamic resistance (ms between height ^-1).

A) calculate each pixel aerodynamic force impedance r _Ah

The neutral atmosphere layer is forged, and aerodynamic resistance can be expressed as:

$r_{\mathrm{ah}}=\frac{1}{{\mathrm{ku}}_{*}}\ln \left(\frac{Z_2}{Z_1}\right)---\left(5\right)$

In the formula, k is von Karman constant (getting 0.41, dimensionless); Z ₂It is the reference altitude (generally getting 2m) of measuring temperature; Z ₁=z _Oh, i.e. heat transmission roughness (m); u _*Be friction velocity (ms ^-1), neutral line is forged and can be provided by the logarithmic velocity profile equation:

$\frac{1}{u_{*}}=\frac{1}{k\stackrel{\&OverBar;}{u}\left(z\right)}\ln \left(\frac{z}{z_{\mathrm{om}}}\right)---\left(6\right)$

In the formula,

Be mean wind speed; Z is the height of mean wind speed correspondence; z _OmIt is face of land power roughness, it is the null height of mean wind speed in the logarithm wind profile formula, theoretical, with the contacted fluid motion speed of solid boundaries speed only being arranged just on solid boundaries is zero, but actual natural terrain fluctuations, making that the boundary layer wind speed just is zero at face of land certain altitude place, also is face of land power roughness.

The real atmosphere situation all is non-neutral line knot by day usually, but not wind number that neutral line is forged and warm and humid profile all are different from the profile that neutral line is forged, because degree of stability is influential to profile.Mo Ning-Ao Bu Hough (Monin-Obukhov) has proposed now the still Monin-Obukhov similarity theory of widespread use, method with similarity theory and dimensional analysis, discuss the influence that shearing stress and buoyancy are carried the ground layer turbulent flow, set up the general expression formula of ground layer meteorological element profile rule.

Dyer (1974) should consider the influence of atmospheric buoyancy force effect (Buoyancy Effect) when pointing out to calculate the sensible heat flux.The atmospheric buoyancy force effect is meant when the air of ground layer is heated, volume produces expansion and rises to higher position, when the temperature lapse rate when if this air mass rises and the temperature lapse rate of external environment produce difference, will produce so-called buoyancy effect: when the temperature lapse rate of ascending air mass and external environment is identical, be called neutral atmosphere layer knot (neutral state); And when the temperature lapse rate of ascending air mass during greater than the temperature lapse rate of external environment, the air mass that is heated will continue to rise and form non-steady state (unstable state), and this state can quicken the speed of vertical convection down; Otherwise if the temperature lapse rate of ascending air mass during less than the temperature lapse rate of external environment, then forms steady state (SS) (stable state), ascending air mass is forced to horizontal motion, and vertical convection speed is suppressed and weakens.

B) correcting friction wind speed and aerodynamic resistance

The present invention adopts Dyer (1974) to provide degree of stability correction function Ψ based on similarity theory _m(x), Ψ _h(x) respectively to friction velocity u _*, aerodynamic resistance r _AhCorrect:

$u_{*}=\frac{\mathrm{ku}\left(z\right)}{\ln \left(\frac{z-d}{z_{\mathrm{om}}}\right)-{\mathrm{\&Psi;}}_m\left(\frac{z-d}{L}\right)+{\mathrm{\&Psi;}}_m\left(\frac{z_{\mathrm{om}}}{L}\right)}---\left(7\right)$

$r_{\mathrm{ah}}=\frac{1}{{\mathrm{ku}}_{*}}\left[\ln \left(\frac{Z_2}{Z_1}\right)-{\mathrm{\&Psi;}}_h\left(\frac{Z_2}{L}\right)+{\mathrm{\&Psi;}}_h\left(\frac{Z_1}{L}\right)\right]---\left(8\right)$

In the formula, L is a Monin-Obukhov length, and it has reflected the situation of ground layer atmospheric turbulence, is used for judging the stability state of atmospheric stratification, and d is the zero plane displacement height, and u (z) is the wind speed at computed altitude place.All the other each parameters together above.Can calculating of Monin-Obukhov length by following formula:

$L=-\frac{{\mathrm{\&rho;}}_a{c}_p{u}_{*}^3{T}_s}{\mathrm{kgH}}---\left(9\right)$

In the formula, ρ _aBe atmospheric density (kgm ^-3); c _pBe the specific heat capacity at constant pressure (1004.07Jkg of air ^-1K ^-1); T _sBe the face of land of remote sensing or the surface temperature of canopy (K); K is von Karman constant (getting 0.41, dimensionless); u _*Be friction velocity (ms ^-1); G is acceleration of gravity (9.807ms ^-2); H is sensible heat flux (Wm ^-2).

Power roughness correction function Ψ _m(x) and heat delivered degree of stability correction function Ψ _h(x) as follows respectively when (L＜0) under unstable atmospheric condition:

${\mathrm{\&Psi;}}_m\left(x\right)=2\ln \left[\frac{1+{(1-16x)}^{1/4}}{2}\right]+\ln \left[\frac{1+{(1-16x)}^{1/2}}{2}\right]---\left(7a\right)$

$-2\arctan \left[{(1-16x)}^{1/4}\right]+\frac{\mathrm{\&pi;}}{2}$

${\mathrm{\&Psi;}}_h\left(x\right)=2\ln \left[\frac{1+{(1-16x)}^{1/2}}{2}\right]---\left(7b\right)$

Stablizing (L＞0) under the atmospheric conditions, correction function is expressed as:

Ψ _m(x)＝-5x (7c)

Ψ _h(x)＝-5x (7d)

Under the neutral atmosphere condition (| L| → ∞, x → 0), then:

ψ _m(0)＝Ψ _h(0)＝0 (7e)

More than various in, x=z/L.

C) sensible heat flux under the impedance of iterative computation stable air dynamics

Can be known by (4) formula, ask sensible heat flux (H) such situation to occur: the molecule and the denominator on equation the right respectively have a unknown parameter dT and r _Ah, in order to find the solution H, there is extreme (extremely dried, extremely wet) pixel on the model hypothesis satellite image, be called " Hot " point and " Cold " point; And have a kind of linear relationship between the temperature at supposition surface temperature and reference altitude place, as boundary condition, introducing Monin-Obukhov similarity theory is tried to achieve H under the impedance of stable air dynamics by iteration with this.Below introduce in detail the step of finding the solution:

Based on the preamble analysis, the wind speed that at first selected 200m highly locates is the wind speed in the narrow sense mixolimnion, its value is relatively stable, therefore supposes that its size is not subjected to the underlay The surface roughness affected, and the view picture image can be tried to achieve same value u for the meteorological website actual measurement of people wind speed according to logarithmic velocity profile ₂₀₀(as have ready conditions also can Mesoscale Meteorology output mixolimnion height place wind field replace):

$u_{200}=\frac{{\mathrm{<figure-callout\; id="10"\; label="u"\; filenames="HSA00000414253100021.png"\; state="\{\{state\}\}">u</figure-callout>}}_{10}\ln (200/z_{\mathrm{om}\_\mathrm{station}})}{\ln (10/z_{\mathrm{om}\_\mathrm{station}})}---\left(10\right)$

In the formula, u ₁₀The wind speed of measuring for meteorological station (is generally 10m and highly locates, ms ^-1); z _{Om_station}For the power roughness around the meteorological site, get 0.06m.

Calculate the power roughness z that image obtains the face of land pixel of correspondence in period _Om(parametric method is face as follows) supposes that the neutral atmosphere layer forges, can be in the hope of the friction wind speed u of all pixels according to (7) formula _*, and then (8) formula of utilization can get the aerodynamic resistance r under the neutrallty condition _Ah

Then, suppose in the master mould:

dT＝aT _s+b (11)

By artificially selecting water body pixel or vegetation coverage height and the very low pixel of temperature is " Cold " point, and supposes that all effective energies of this point all are used to (the LE that evapotranspires _Cold=R _Ncold-G _Cold), i.e. H _Cold=0, dT _Cold=0; The selected pixel that vegetation coverage is low and temperature is higher is " Hot " point, supposes this some generation (LE that do not evapotranspire _Hot=0), H then _Hot=R _Nhot-G _Hot, dT _Hot=(R _Nhot-G _Hot) r _Ahhot/ ρ _ac _pCan obtain then:

$a=\frac{{\mathrm{dT}}_{\mathrm{hot}}}{T_{\mathrm{hot}}-T_{\mathrm{cold}}},b=\frac{-{\mathrm{dT}}_{\mathrm{hot}}{T}_{\mathrm{cold}}}{T_{\mathrm{hot}}-T_{\mathrm{cold}}}---\left(12\right)$

But in fact, there are some problems in above-mentioned hypothesis, " Cold " some maximum LE that evapotranspires _ColdShould equal potential evapotranspiration λ ET _r, research experiment is also pointed out under vegetation growth busy season and water supply sufficiency, LE _ColdUsually also can be than potential evapotranspiration λ ET _rBig approximately by about 5%, so supposition LE _Cold=1.05 λ ET _rConsidering that image obtains constantly may vegetation be in the non-season of growth or image resolution too big " Cold " point and " Hot " point has mixed mixed pixel, maximum is evapotranspired and minimum is evapotranspired can not reach above-mentioned supposition, can suitably adjust by survey region NDVI value:

$\left\{\begin{array}{c}{H}_{\mathrm{cold}}=R_{\mathrm{ncold}}-G_{\mathrm{cold}}-k_{\mathrm{cold}}{\mathrm{\&lambda;ET}}_r\\ H_{\mathrm{hot}}=R_{\mathrm{nhot}}-G_{\mathrm{hot}}-k_{\mathrm{hot}}{\mathrm{\&lambda;ET}}_r\end{array}\right.---\left(13\right)$

Wherein:

$k_{\mathrm{cold}}=\left\{\begin{array}{cc}1.05-\frac{0.85-\mathrm{<figure-callout\; id="2"\; label="NDVI"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">NDVI</figure-callout>}}{2},& \mathrm{NDVI}<0.85\\ 1.05,& \mathrm{NDVI}\&GreaterEqual;0.85\end{array}\right.$ $k_{\mathrm{hot}}=\left\{\begin{array}{cc}\mathrm{NDVI}-0.15,& \mathrm{NDVI}>0.15\\ 0,& \mathrm{NDVI}\&le;0.15\end{array}\right.$

Potential evapotranspiration ET in the formula _rThe meteorological data that can utilize satellite to pass by constantly obtains by the Penman-Monteith Equation for Calculating:

${\mathrm{ET}}_r=\frac{0.408\mathrm{\&Delta;}(R_n-G)+\mathrm{\&gamma;}\frac{900}{T_a+273}{u}_2(e_s-e_a)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}(1+0.34u_2)}---\left(13a\right)$

In the formula, Δ is saturation vapour pressure-temperature curve curvature; γ is the wet and dry bulb constant; u ₂Be that 2m highly locates wind speed; e _sBe saturation vapour pressure; e _aFor actual vapor is pressed; The same preamble of all the other parameters; It should be noted that need be potential evapotranspiration λ ET _rUnit be converted to (Wm ^-2), λ=[2.501-0.00236 (T _a-273.15)] * 10 ⁶(Jkg ^-1).As obtaining corresponding meteorological data constantly, then calculate (Priestley and Taylor, 1972) with the Priestley-Taylor formula:

${\mathrm{\&lambda;ET}}_r=\mathrm{\&alpha;}\frac{\mathrm{\&Delta;}(R_n-G)}{\mathrm{\&Delta;}+\mathrm{\&gamma;}}---\left(13b\right)$

In the formula, α is the Priestley-Taylor coefficient, gets 1.26～1.29; Δ is saturation vapour pressure-temperature curve curvature (PaK ^-1):

$\mathrm{\&Delta;}=\frac{4098\left[610.8\exp \left(\frac{17.27T_a}{T_a+237.3}\right)\right]}{{(T_a+237.3)}^2}---\left(13c\right)$

γ is wet and dry bulb constant (PaK ^-1):

$\mathrm{\&gamma;}=\frac{1004.07\&times;P}{0.622\&times;[2.501-0.00236(T_a-273.15)]\&times;{10}^6}---\left(13d\right)$

P is atmospheric pressure (Pa):

$P=1013.25\&times;100/{10}^{\frac{h}{18410(T_a/273.15)}---\left(13e\right)}$

So dT is revised as:

$\left\{\begin{array}{c}{\mathrm{dT}}_{\mathrm{cold}}=H_{\mathrm{cold}}{r}_{\mathrm{ahcold}}/{\mathrm{\&rho;}}_{\mathrm{acold}}{c}_p\\ {\mathrm{dT}}_{\mathrm{hot}}=H_{\mathrm{hot}}{r}_{\mathrm{ahhot}}/{\mathrm{\&rho;}}_{\mathrm{ahot}}{c}_p\end{array}\right.---\left(14\right)$

Wherein, atmospheric density can be obtained by the equation of gas state:

${\mathrm{\&rho;}}_a=\frac{P}{287.05\&times;(1+0.608q)\&times;T_a}\&ap;\frac{P}{287.05\&times;1.01\&times;T_a}---\left(14a\right)$

A and b also are revised as thereupon:

$a=\frac{{\mathrm{dT}}_{\mathrm{hot}}-{\mathrm{dT}}_{\mathrm{cold}}}{T_{\mathrm{hot}}-T_{\mathrm{cold}}},b=\frac{{\mathrm{dT}}_{\mathrm{cold}}{T}_{\mathrm{hot}}-{\mathrm{dT}}_{\mathrm{hot}}{T}_{\mathrm{cold}}}{T_{\mathrm{hot}}-T_{\mathrm{cold}}}---\left(12a\right)$

D) based on the calculating of the surface temperature under the DEM

Simultaneously, consider surface temperature that remote sensing obtains with height change trend substantially with the lapse rate of atmospheric temperature unanimity, increase with height reduces, and the surface temperature of whole survey region is all adjusted to same reference altitude (regional average height) change the influence of calculating evapotranspiring at last to eliminate the surface temperature that causes owing to elevation:

T _{s_dem}＝T _s+0.0065(h-h _mean) (15)

In the formula, h is a sea level elevation, h _MeanBe the sea level on the average of survey region, T _sAir themperature for the reference altitude place.

Can determine two height Z by above calculating ₂With Z ₁Between temperature difference dT:

dT＝aT _{s_dem}+b (16)

The air themperature that satellite passes by is constantly also determined thereupon:

T _a＝T _s-dT (17)

So far, dT and r _AhAll definite, just can obtain the sensible heat flux according to (4) formula.But just as mentioned before, atmosphere is the unstable stratification condition usually, therefore adopts the Monin-Obukhov similarity theory, introduces the Monin-Obukhov length L, with r _Ah, H, dT and atmospheric density ρ _aConnect, by iteration repeatedly, until r _AhReach stable, obtain final sensible heat flux H.

4. calculate instantaneous latent heat (LE) based on energy-balance equation

Above calculation procedure can reduce following iterative computation process flow diagram.

At last, can obtain instantaneous latent heat flux based on energy-balance equation:

LE＝R _n-G-H (18)

5. based on the focus of feature space and determining automatically of cold spot

In the said process, need artificial selected " Hot " point and " Cold " point, have artificial subjectivity, and increased workload, be not easy to be used for handling daily Monitoring Data (as MODIS).(Vegetation Index/Temperature Trapezoid VITT) helps solve the selected automatically of " Hot " point and " Cold " point can to utilize the feature space of the vegetation index-surface temperature of the popular research of Chinese scholars in recent years.

Fig. 3 is meant the VITT space synoptic diagram that is constituted by index and surface temperature, vegetation index in the corresponding diagram 3, different scholar's research the space characteristics that constitutes of many different vegetation indexs (NDVI, LAI, SAVI etc.) and surface temperature (LST).Therefore can select " Hot " some pixel and " Cold " some pixel in conjunction with vegetation index and surface temperature.NDVI and LST combination are adopted in suggestion such as Verstraeten, and Allen etc. then utilize the combination of LAI and LST to help reconnaissance.Work (the corresponding exposed soil fully in suggestion SAVI～0.1 with reference to Moran etc., corresponding vegetation fully covers in SAVI～0.8), the present invention adopts and suppresses the more satisfactory adjustable vegetation index of modified soil (MSAVI) of Soil Background, in conjunction with LST, provide the criterion of the automatic selection of " Hot " point and " Cold " point: the MSAVI about 0.1 reaches wherein combination correspondence " Hot " point of maximum LST; 0.8 about MSAVI and combination correspondence " Cold " point of wherein minimum LST.Extract the correlation parameter of corresponding point then, realize the automatic operation under unmanned the participation.

6. dynamically estimate face of land kinetic parameter

In the surface water thermal diffusion process, face of land kinetic parameter mainly is meant face of land power roughness z _Om, zero plane displacement d (Hussein et al.2001).The land Data Assimilation system of U.S. NASA Goddard space research center (calls in the following text: assimilation system) at the seasonal variations of surface vegetation, adopt the parameter look-up table of every month feature of a kind of vegetation.Usually, available experimental formula (19), (20) and (21) estimate zero plane displacement and roughness height by the surface effective height:

d＝0.667H _eff (19)

z _om＝0.136H _eff (20)

z _oh＝0.1z _om (21)

In the formula, H _EffBe the pixel significant height, promptly considered the interior atural object significant height of pixel of factors such as component height, proportion of composing and spatial distribution characteristic such as vegetation, atural object in the pixel.

Can obtain heat delivered roughness z by (21) formula then _OhThe present invention introduces remote sensing vegetation index, a kind of empirical method of dynamically estimating pixel effective coverage height of structure on this look-up table method basis.

$H_{\mathrm{eff}}=H_{\min }+\frac{\mathrm{VI}-{\mathrm{VI}}_{\min }}{{\mathrm{VI}}_{\max }-{\mathrm{VI}}_{\min }}\&times;\mathrm{\&Delta;H}---\left(22\right)$

In the formula, H _EffBe pixel effective coverage height, VI is the pixel vegetation index, H _Max, H _MinBe respectively maximum, minimum effective height, Δ H is the poor of minimax significant height, VI _Max, VI _MinThen be respectively maximum, minimum vegetation index.Each vegetation pattern significant height and vegetation index are worth with reference to U.S. NASA land face Data Assimilation system face of land parameter look-up table most to be determined, to then not getting constant surface effective height with other atural objects of seasonal variations, this research is adopted is 1,/10 ten thousand land use data.

The table 1 dynamically vegetation index and the significant height of estimation face of land kinetic parameter is worth look-up table most

By the empirical parameter of formula (21) and table 1, utilize remote sensing MSAVI can estimate dynamic surface effective height, and then can calculate zero plane displacement and roughness of ground surface by experimental formula (19), (20).Can just make full use of remote sensing observations broad perspectives and ageing advantage like this improves the estimation of face of land kinetic parameter and the realization of SEBTA algorithm.

7. based on surface flux zoning lack of water index

The present invention is on crop water shortage index (CWSI) mechanism based, the surface flux that utilization calculates, designed regional lack of water index (the Regional Water Stress Index in the suitable China north, RWSI), as the index and the reference frame of regional water surplus situation, carry out the monitoring of regional soil moisture and damage caused by a drought.Its form is:

$\mathrm{RSWI}=1-\frac{\mathrm{ET}}{{\mathrm{ET}}_{\mathrm{wet}}}---\left(23a\right)$

In the formula, ET is regional actual evapotranspiration, ET _WetBe evapotranspiring under the regional water supply sufficiency, promptly regional potential evapotranspiration.Potential evapotranspiration with reference to evapotranspiring, is that the maximum on the face of land under the desirable water supply conditions is evapotranspired promptly, supposes the sensible heat flux minimum (H ≈ 0) under this state, and all effective energies that the face of land receives all are used to evapotranspire, i.e. λ ET _Wet=R _n-G, the simultaneous energy-balance equation obtains:

$\mathrm{RWSI}=1-\frac{\mathrm{\&lambda;ET}}{{\mathrm{\&lambda;ET}}_{\mathrm{wet}}}=\frac{H}{R_n-G}---\left(23b\right).$

In the formula, H is the sensible heat flux, R _nBe net radiation flux, G is a soil heat flux, all can be by obtaining in the region in front surface water thermoflux Remote Sensing Model.Therefore can monitor in real time the regional moisture status of profit and loss from the angle of remote sensing, and then can monitor developing of arid.

Utilize the national standard of the damage caused by a drought grade classification of soil moisture, the arid ranking score grade standard of the RWSI of each time period of converting is used for regional arid remote sensing monitoring.

With reference to agriculture damage caused by a drought grading standard: relative moisture of the soil＜40% drought of attaching most importance to, 40%～50% is middle drought, 50%～60% light drought, 60%～80% is normal, 80%～100% is moistening, is 5 grade with the RWSI value with the degree of water shortage grade classification of survey region according to this standard, and the criteria for classifying sees Table in detail.

The regional arid grading standard of table 2

Corresponding to the drought monitoring method based on surface water thermoflux remote-sensing inversion of the present invention, draught monitor system based on surface water thermoflux remote-sensing inversion also is provided, Fig. 4 is that the draught monitor system that the present invention is based on surface water thermoflux remote-sensing inversion gets structural representation, as shown in Figure 4, described system comprises:

Surface net radiation flux computing module 1 is used to calculate surface net radiation flux $R_n=(1-\mathrm{\&alpha;})R_{s\&DownArrow;}+{\mathrm{\&epsiv;}}_s\mathrm{\&sigma;}({\mathrm{\&epsiv;}}_a{T}_a^4-T_s^4);$

Wherein, R _{S ↓}Be the sun shortwave radiation of incident, α is a surface albedo, ε _sBe face of land emissivity; σ is a Stefan-Boltzmann constant (5.6696 * 10 ^-8Wm ^-2K ^-4), T _sIt is surface temperature; T _aBe the air themperature at reference altitude place, ε _aIt is the air ratio radiance;

Soil heat flux computing module 2 is used for based on the remote sensing parametric inversion, calculates soil heat flux:

$G=\left[\frac{(T_s-273.15)}{\mathrm{\&alpha;}}(0.0038\mathrm{\&alpha;}+0.0074{\mathrm{\&alpha;}}^2)(1-0.98{\mathrm{NDVI}}^4)\right]R_n$

Wherein, α is a surface albedo, T _sBe surface temperature, NDVI is the water body of vegetation normalization index for study area, because of its between limpid deep water body and wetland, get G _Water=0.5R _n

Sensible heat flux computing module 3. is used to utilize iterative calculation method, calculates sensible heat flux:

$H={\mathrm{\&rho;}}_a{c}_p\frac{(T_s-T_a)}{r_{\mathrm{ah}}}={\mathrm{\&rho;}}_a{c}_p\frac{\mathrm{dT}}{r_{\mathrm{ah}}}$

Wherein, ρ _aBe atmospheric density (kgm ^-3); c _pBe the specific heat capacity at constant pressure (1004.07Jkg of air ^-1K ^-1); T _sBe the face of land of remote sensing or the surface temperature of canopy; T _aAir themperature (K) for the reference altitude place; DT is two height Z ₂With Z ₁Between the temperature difference; r _AhBe two aerodynamic resistance (ms between height ^-1);

Zone lack of water Index for Calculation module 4. is used for the calculating according to described surface flux, zoning lack of water index

According to regional lack of water index RWSI, carry out the monitoring of regional degree of drought.

Wherein, described surface net radiation flux computing module 1 comprises:

Surface albedo computing module 11 is used to calculate surface albedo α;

MODIS The data following formula calculates surface albedo:

α ^MODIS＝0.160α ₁+0.291α ₂+0.243α ₃+0.116α ₄

+0.112α ₅+0.081α ₇-0.0015

In the formula, α _i, i=1 ... each narrow wave band surface albedo on the corresponding satellite sensor of n difference;

Face of land emissivity computing module 12 is used to calculate face of land emissivity ε _s

Extract water body earlier, compose typical emissivity value 0.995 with water body; The emissivity estimation equation of remaining mixed pixel is as follows:

${\mathrm{\&epsiv;}}_s=\left\{\begin{array}{c}0.9624744+0.0652704P_v-0.0461286{\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v<0.5\\ 0.9643744+0.0614704P_v-0.0461286{\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v=0.5\\ 0.9662744+0.0576704P_v-{0.0461286\mathrm{<figure-callout\; id="2"\; label="P\; v"\; filenames="HSA00000414253100011.png,HSA00000414253100021.png"\; state="\{\{state\}\}">P</figure-callout>}}_v^{},P_v>0.5\end{array}\right.$

Vegetation coverage is provided by following formula:

$P_v=\frac{(\mathrm{NDVI}-{\mathrm{NDVI}}_{\min })}{({\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min })}$

In the formula, NDVI _Min, NDVI _Max, NDVI corresponds to pure exposed soil pixel respectively, the NDVI value of pure vegetation pixel and mixed pixel correspondence is got NDVI _Min=0.099, NDVI _Max=0.77; When NDVI smaller or equal to 0.099 the time, vegetation coverage P _v=0; When NDVI more than or equal to 0.77 the time, vegetation coverage P _v=1;

Sun shortwave radiation computing module 13 is used to calculate the sun shortwave radiation R of incident _{S ↓}

The sun shortwave radiation in a certain moment can be regarded direct solar radiation, scattered radiation and reflected radiation three parts as and form:

R _s↓＝R _s+R _d+R _r

In the formula, R _sBe direct radiation; R _dBe scattered radiation; R _rBe reflected radiation.Wherein:

$R_s={\left(\frac{1}{\mathrm{\&rho;}}\right)}^2{I}_0\cos {\mathrm{\&theta;\&tau;}}_b={\mathrm{<figure-callout\; id="0"\; label="R"\; filenames="HSA00000414253100042.png"\; state="\{\{state\}\}">R</figure-callout>}}_0{\mathrm{\&tau;}}_b\cos \mathrm{\&theta;}$

Wherein, I ₀Be solar constant, get 1367 (Wm ^-2); R ₀Solar radiation for the straight incident of atmospheric envelope roof pendant; θ is the solar zenith angle on domatic; τ _bBe atmosphere direct projection transmissivity; (1/ ρ) ²Be called solar distance correction coefficient or Earth's orbit and correct the factor, about 10 ^-4During precision, calculate by fourier progression expanding method:

R _d＝R ₀τ _dcos ²β/2sina

In the formula, β is the gradient, and a is a sun altitude.

R _r＝rR ₀τ _rsin ²β/2sina

In the formula, r is the first reflectivity in slope, the face of land, can simplify being taken as 0.2; τ _rTransmissivity for reflected radiation.

Wherein, described sensible heat flux computing module 3 comprises:

Correcting module 31 is used for providing degree of stability correction function Ψ based on similarity theory _m(x), Ψ _h(x), respectively to friction velocity u _*, aerodynamic resistance r _AhRevise:

$u_{*}=\frac{\mathrm{ku}\left(z\right)}{\ln \left(\frac{z-d}{z_{\mathrm{om}}}\right)-{\mathrm{\&Psi;}}_m\left(\frac{z-d}{L}\right)+{\mathrm{\&Psi;}}_m\left(\frac{z_{\mathrm{om}}}{L}\right)}$

$r_{\mathrm{ah}}=\frac{1}{{\mathrm{ku}}_{*}}\left[\ln \left(\frac{Z_2}{Z_1}\right)-{\mathrm{\&Psi;}}_h\left(\frac{Z_2}{L}\right)+{\mathrm{\&Psi;}}_h\left(\frac{Z_1}{L}\right)\right]$

Wherein, k is von Karman constant (getting 0.41, dimensionless), and z is the height of mean wind speed correspondence, z _OmBe face of land power roughness, Z ₂It is the reference altitude of measuring temperature; Z ₁=z _Oh, the heat delivered roughness; u _*Be friction velocity (ms ^-1), L is a Monin-Obukhov length, In the formula, ρ _aBe atmospheric density; c _pIt is the specific heat capacity at constant pressure of air; T _sBe the face of land of remote sensing or the surface temperature of canopy; K is a von Karman constant; u _*It is friction velocity; G is an acceleration of gravity; H is the sensible heat flux;

Surface temperature computing module 32 is used to calculate the surface temperature based under the DEM,

T _{S_dem}=T _s+ 0.0065 (h-h _Mean), wherein, h is a sea level elevation, h _MeanSea level on the average for survey region;

Pixel is chosen module 33, is used to choose extremely dried pixel and extremely wet pixel, with described T _{S_dem}Substitution formula dT=aT _{S_dem}+ b, wherein,

T _HotThe temperature of the extremely dried pixel of expression, T _ColdThe temperature of the extremely wet pixel of expression.

Described pixel is chosen in the module 33, and the described criterion of choosing extremely dried pixel and extremely wet pixel is: the MSAVI about 0.1 reaches wherein combination correspondence " Hot " point of maximum LST; 0.8 about MSAVI and combination correspondence " Cold " point of wherein minimum LST.

Described regional lack of water Index for Calculation module 4 based on the Regional Drought index of RWSI, is utilized the national standard of the damage caused by a drought grade classification of soil moisture, and the arid ranking score grade standard of the RWSI of each time period of converting carries out regional arid remote sensing monitoring.

DEM is to the sensitivity analysis of surface flux influence

In SEBTA model of the present invention, DEM is introduced the computation process of correlation parameter, taken into full account the influence of the topographical features exterior heat amount over the ground and the ET that evapotranspires.Sea level elevation, the gradient, the aspect influence degree to SEBTA modeling result of the present invention will be described respectively below.Concrete way is respectively with considering the influence of topography and do not consider that the Model Calculation of action of topography goes out surface water thermoflux (LE), then in conjunction with the two difference of landform factorial analysis.The height of test site is between 8m～1532m, with every 100m classification, calculate the consideration action of topography of the correspondence in each elevation segment limit respectively and do not consider the average of the following day evapotranspiration of action of topography, Fig. 5 is the graph of a relation of considering (a) and not considering elevation and day evapotranspiration under (b) action of topography.Can find that from Fig. 5 the two is more similar with the curve distribution of height, be one " spoon " type and distribute that this distribution trend is relevant with the soil utilization/soil cover type of test site.Owing to the plant growth in season, water consumption is big in the plains region; And mostly the sea level elevation scope about 200m is that vegetation covers the exposed soil ground of non-constant, thereby the minimum of evapotranspiring; Between 200m～1000m, the vegetation coverage height, vegetation distributes than comparatively dense, and surface temperature reduces gradually, and solar radiation adds wind action on the highland with increasing highly gradually, evapotranspires substantially with highly increasing.This shows that the model of not considering action of topography obviously over-evaluated evapotranspiring of mountain area.

The gradient and aspect are the important indicators of topographic relief feature, aspect mainly influences the solar radiation of ground acceptance and the angle of cut of ground and prevailing wind direction, this makes and has significant hydro-thermal difference between the different aspects, the increase of the gradient can aggravate this influence, thereby causes the difference of surface water thermoflux space distribution.Mainly studied the test site gradient below greater than the influence of aspect domatic more than 5 degree to face of land water and heat space distribution.

Fig. 6 is the influence curve figure of aspect to face of land water and heat, net radiation flux is influenced by aspect obviously as can be seen from Figure 6, mxm. appears at 120 ° aspect approximately, it is very to coincide that this and image obtain solar azimuth (118.7 °) constantly, because sun direct projection this moment is on aspect is 120 ° domatic, radiant quantity maximum, the domatic radiation that receives of the 300 ° aspects relative with it are then minimum.The net radiation flux of not considering influence of topography estimation is big or small homogeneous on all aspects almost, and this obviously is incorrect in the mountain area.The difference of received radiation energy also directly is reflected on the evapotranspiring of different aspects, and the beam radia that tailo (Nan Po is a north slope in the Southern Hemisphere) is accepted is many, causes domatic temperature height, water evaporates strong.On the contrary, evapotranspiring of Schattenseite (north slope, the Southern Hemisphere are Nan Po) is slightly littler than Nan Po, and this considers also obviously to embody in the day distribution of evapotranspiration with aspect of action of topography in Fig. 4.This shows that the flat earth model of not considering aspect, gradient effect over-evaluated evapotranspiring of Schattenseite.

This shows that terrain factor (sea level elevation, the gradient, aspect) is remarkable to SEBTA modeling result's influence, will not over-evaluate the face of land heat and the evapotranspiration of high height above sea level place and Schattenseite, cause very big estimation error if consider the influence of topography.The SEBTA model is that SEBTA not only can use in the plains region because DEM is introduced the CALCULATION OF PARAMETERS process, also can use in the mountain area, has enlarged the usable range of SEBTA model.

In sum, because drought monitoring method and the system based on surface water thermoflux remote-sensing inversion of the present invention introduced the factor of influence of DEM and face of land covering, make the present invention can carry out heat under the complicated earth surface of big zone (the validation test zone is 660,000 Km2) and the simulation estimate precision of ET can reach the practicability precision, expanded the usable range and the degree of being practical of SEBTA model greatly.

Beneficial effect of the present invention is:

1. DEM and face of land covering (mainly are meant face of land power roughness z in conjunction with carrying out face of land kinetic parameter _Om, zero plane displacement d) dynamic calculation;

2. the space characteristics that utilizes adjustable vegetation index of modified soil (MSAVI) and LST to constitute is realized the automatic examination of doing wet pixel of simulated domain;

3. in the estimation of DEM and the face of land Cover Characteristics factor being introduced face of land energy and evapotranspiring, make SEBTA (Surface Energy Balance with Topography Algorithm) can carry out the face of land energy of the long-time sequence big regional complicated earth surface under and the estimation of evapotranspiring is simulated;

Drought monitoring method and system employs principle of energy balance based on surface water thermoflux remote-sensing inversion of the present invention, taking into full account the complex-terrain and the face of land covers face of land energy and the influence of evapotranspiring, the DEM and the face of land are covered in the introducing algorithm, make this algorithm not be only applicable to the calculating of subdued topography, also applicable to hills, mountain area face of land water and heat and draught monitor CALCULATION OF PARAMETERS, and by the dynamic surface effective high computational roughness of ground surface of estimation, the influence that complicated earth surface covers is introduced in the model, make algorithm can carry out the different Parameters Calculation that cover under (comprising water body), determine automatically to do, wet pixel.So not only expand the usable range of algorithm, and improved the practicality and the operability of algorithm.

In conjunction with the drawings to the description of the specific embodiment of the invention, others of the present invention and feature are conspicuous to those skilled in the art.

More than specific embodiments of the invention are described and illustrate it is exemplary that these embodiment should be considered to it, and be not used in and limit the invention, the present invention should make an explanation according to appended claim.

Claims (10)

1. based on the drought monitoring method of surface water thermoflux remote-sensing inversion, it is characterized in that described method comprises the following steps:

Step 100. is calculated surface net radiation flux

Wherein, R _{S ↓}Be the sun shortwave radiation of incident, α is a surface albedo, ε _sBe face of land emissivity; σ is a Stefan-Boltzmann constant (5.6696 * 10 ^-8Wm ^-2K ^-4), T _sIt is surface temperature; T _aBe the air themperature at reference altitude place, ε _aIt is the air ratio radiance;

Step 200. is calculated soil heat flux based on the remote sensing parametric inversion:

$G=\left[\frac{(T_s-273.15)}{\mathrm{\&alpha;}}(0.0038\mathrm{\&alpha;}+0.0074{\mathrm{\&alpha;}}^2)(1-0.98\mathrm{NDV}{I}^4)\right]R_n$

Wherein, α is a surface albedo, T _sBe surface temperature, NDVI is the water body of vegetation normalization index for study area, because of its between limpid deep water body and wetland, get G _Water=0.5R _n

Step 300. is utilized iterative calculation method, calculates sensible heat flux:

$H={\mathrm{\&rho;}}_a{c}_p\frac{(T_s-T_a)}{r_{\mathrm{ah}}}={\mathrm{\&rho;}}_a{c}_p\frac{\mathrm{dT}}{r_{\mathrm{ah}}}$

Wherein, ρ _aBe atmospheric density (kgm ^-3); c _pBe the specific heat capacity at constant pressure (1004.07Jkg of air ^-1K ^-1); T _sBe the face of land of remote sensing or the surface temperature of canopy; T _aAir themperature (K) for the reference altitude place; DT is two height Z ₂With Z ₁Between the temperature difference; r _AhBe two aerodynamic resistance (ms between height ^-1);

Step 400. is according to the calculating of described surface flux, zoning lack of water index

According to regional lack of water index RWSI, carry out the monitoring of regional degree of drought.

2. the drought monitoring method based on surface water thermoflux remote-sensing inversion according to claim 1 is characterized in that, described step 100 comprises the following steps: step 110. calculating surface albedo α;

MODIS The data following formula calculates surface albedo:

α ^MODIS＝0.160α ₁+0.291α ₂+0.243α ₃+0.116α ₄

+0.112α ₅+0.081α ₇-0.0015

In the formula, α _i, i=1 ... each narrow wave band surface albedo on the corresponding satellite sensor of n difference;

Step 120. is calculated face of land emissivity ε _s

Extract water body earlier, compose typical emissivity value 0.995 with water body; The emissivity estimation equation of remaining mixed pixel is as follows:

${\mathrm{\&epsiv;}}_s=\left\{\begin{array}{c}0.9624744+0.0652704P_v-0.0461286P_v^2,P_v<0.5\\ 0.9643744+0.0614704P_v-0.0461286P_v^2,P_v=0.5\\ 0.9662744+0.0576704P_v-{0.0461286P}_v^2,P_v>0.5\end{array}\right.$

Vegetation coverage is provided by following formula:

$P_v=\frac{(\mathrm{NDVI}-{\mathrm{NDVI}}_{\min })}{({\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min })}$

In the formula, NDVI _Min, NDVI _Max, NDVI corresponds to pure exposed soil pixel respectively, the NDVI value of pure vegetation pixel and mixed pixel correspondence is got NDVI _Min=0.099, NDVI _Max=0.77; When NDVI smaller or equal to 0.099 the time, vegetation coverage P _v=0; When NDVI more than or equal to 0.77 the time, vegetation coverage P _v=1;

Step 130. is calculated the sun shortwave radiation R of incident _{S ↓}

The sun shortwave radiation in a certain moment can be regarded direct solar radiation, scattered radiation and reflected radiation three parts as and form:

R _s↓＝R _s+R _d+R _r

In the formula, R _sBe direct radiation; R _dBe scattered radiation; R _rBe reflected radiation.Wherein:

$R_s={\left(\frac{1}{\mathrm{\&rho;}}\right)}^2{I}_0\cos {\mathrm{\&theta;\&tau;}}_b=R_0{\mathrm{\&tau;}}_b\cos \mathrm{\&theta;}$

Wherein, I ₀Be solar constant, get 1367 (Wm ^-2); R ₀Solar radiation for the straight incident of atmospheric envelope roof pendant; θ is the solar zenith angle on domatic; τ _bBe atmosphere direct projection transmissivity; (1/ ρ) ²Be called solar distance correction coefficient or Earth's orbit and correct the factor, about 10 ^-4During precision, calculate by fourier progression expanding method:

R _d＝R ₀τ _d?cos ²β/2sina

In the formula, β is the gradient, and a is a sun altitude.

R _r＝rR ₀τ _rsin ²β/2sina

In the formula, r is the first reflectivity in slope, the face of land, can simplify being taken as 0.2; τ _rTransmissivity for reflected radiation.

3. the drought monitoring method based on surface water thermoflux remote-sensing inversion according to claim 1 is characterized in that described step 300 comprises the following steps:

Step 310. provides degree of stability correction function Ψ based on similarity theory _m(x), Ψ _h(x), respectively to friction velocity u _*, aerodynamic resistance r _AhRevise:

$u_{*}=\frac{\mathrm{ku}\left(z\right)}{\ln \left(\frac{z-d}{z_{\mathrm{om}}}\right)-{\mathrm{\&Psi;}}_m\left(\frac{z-d}{L}\right)+{\mathrm{\&Psi;}}_m\left(\frac{z_{\mathrm{om}}}{L}\right)}$

$r_{\mathrm{ah}}=\frac{1}{{\mathrm{ku}}_{*}}\left[\ln \left(\frac{Z_2}{Z_1}\right)-{\mathrm{\&Psi;}}_h\left(\frac{Z_2}{L}\right)+{\mathrm{\&Psi;}}_h\left(\frac{Z_1}{L}\right)\right]$

Wherein, k is von Karman constant (getting 0.41, dimensionless), and z is the height of mean wind speed correspondence, z _OmBe face of land power roughness, Z ₂It is the reference altitude of measuring temperature; Z ₁=z _Oh, the heat delivered roughness; u _*Be friction velocity (ms ^-1), L is a Monin-Obukhov length,

In the formula, ρ _aBe atmospheric density; c _pIt is the specific heat capacity at constant pressure of air; T _sBe the face of land of remote sensing or the surface temperature of canopy; K is a von Karman constant; u _*It is friction velocity; G is an acceleration of gravity; H is the sensible heat flux;

Step 320. is calculated based on the surface temperature under the DEM, T _{S_dem}=T _s+ 0.0065 (h-h _Mean), wherein, h is a sea level elevation, h _MeanSea level on the average for survey region;

Step 330. is chosen extremely dried pixel and extremely wet pixel, with described T _{S_dem}Substitution formula dT=aT _{S_dem}+ b, wherein,

T _HotThe temperature of the extremely dried pixel of expression, T _ColdThe temperature of the extremely wet pixel of expression.

4. the drought monitoring method based on surface water thermoflux remote-sensing inversion according to claim 3, it is characterized in that, in the step 330, the described criterion of choosing extremely dried pixel and extremely wet pixel is: the MSAVI about 0.1 reaches wherein combination correspondence " Hot " point of maximum LST; 0.8 about MSAVI and combination correspondence " Cold " point of wherein minimum LST.

5. the drought monitoring method based on surface water thermoflux remote-sensing inversion according to claim 1 is characterized in that described step 400 comprises the following steps:

Based on the Regional Drought index of RWSI, utilize the national standard of the damage caused by a drought grade classification of soil moisture, the arid ranking score grade standard of the RWSI of each time period of converting carries out regional arid remote sensing monitoring.

6. based on the draught monitor system of surface water thermoflux remote-sensing inversion, it is characterized in that described system comprises:

The surface net radiation flux computing module is used to calculate surface net radiation flux $R_n=(1-\mathrm{\&alpha;})R_{s\&DownArrow;}+{\mathrm{\&epsiv;}}_s\mathrm{\&sigma;}({\mathrm{\&epsiv;}}_a{T}_a^4-T_s^4);$

Wherein, R _{S ↓}Be the sun shortwave radiation of incident, α is a surface albedo, ε _sBe face of land emissivity; σ is a Stefan-Boltzmann constant (5.6696 * 10 ^-8Wm ^-2K ^-4), T _sIt is surface temperature; T _aBe the air themperature at reference altitude place, ε _aIt is the air ratio radiance;

The soil heat flux computing module is used for based on the remote sensing parametric inversion, calculates soil heat flux:

$G=\left[\frac{(T_s-273.15)}{\mathrm{\&alpha;}}(0.0038\mathrm{\&alpha;}+0.0074{\mathrm{\&alpha;}}^2)(1-0.98\mathrm{NDV}{I}^4)\right]R_n$

Wherein, α is a surface albedo, T _sBe surface temperature, NDVI is the water body of vegetation normalization index for study area, because of its between limpid deep water body and wetland, get G _Water=0.5R _n

The sensible heat flux computing module; Be used to utilize iterative calculation method, calculate sensible heat flux:

$H={\mathrm{\&rho;}}_a{c}_p\frac{(T_s-T_a)}{r_{\mathrm{ah}}}={\mathrm{\&rho;}}_a{c}_p\frac{\mathrm{dT}}{r_{\mathrm{ah}}}$

Wherein, ρ _aBe atmospheric density (kgm ^-3); c _pBe the specific heat capacity at constant pressure (1004.07Jkg of air ^-1K ^-1); T _sBe the face of land of remote sensing or the surface temperature of canopy; T _aAir themperature (K) for the reference altitude place; DT is two height Z ₂With Z ₁Between the temperature difference; r _AhBe two aerodynamic resistance (ms between height ^-1);

Zone lack of water Index for Calculation module is used for the calculating according to described surface flux, zoning lack of water index

According to regional lack of water index RWSI, carry out the monitoring of regional degree of drought.

7. the draught monitor system based on surface water thermoflux remote-sensing inversion according to claim 6 is characterized in that, described surface net radiation flux computing module comprises:

The surface albedo computing module is used to calculate surface albedo α;

MODIS The data following formula calculates surface albedo:

α ^MODIS＝0.160α ₁+0.291α ₂+0.243α ₃+0.116α ₄

+0.112α ₅+0.081α ₇-0.0015

In the formula, α _i, i=1 ... each narrow wave band surface albedo on the corresponding satellite sensor of n difference;

Face of land emissivity computing module is used to calculate face of land emissivity ε _s

Extract water body earlier, compose typical emissivity value 0.995 with water body; The emissivity estimation equation of remaining mixed pixel is as follows:

${\mathrm{\&epsiv;}}_s=\left\{\begin{array}{c}0.9624744+0.0652704P_v-0.0461286P_v^2,P_v<0.5\\ 0.9643744+0.0614704P_v-0.0461286P_v^2,P_v=0.5\\ 0.9662744+0.0576704P_v-{0.0461286P}_v^2,P_v>0.5\end{array}\right.$

Vegetation coverage is provided by following formula:

$P_v=\frac{(\mathrm{NDVI}-{\mathrm{NDVI}}_{\min })}{({\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min })}$

In the formula, NDVI _Min, NDVI _Max, NDVI corresponds to pure exposed soil pixel respectively, the NDVI value of pure vegetation pixel and mixed pixel correspondence is got NDVI _Min=0.099, NDVI _Max=0.77; When NDVI smaller or equal to 0.099 the time, vegetation coverage P _v=0; When NDVI more than or equal to 0.77 the time, vegetation coverage P _v=1;

Sun shortwave radiation computing module is used to calculate the sun shortwave radiation R of incident _{S ↓}

The sun shortwave radiation in a certain moment can be regarded direct solar radiation, scattered radiation and reflected radiation three parts as and form:

P _s↓＝R _s+R _d+R _r

In the formula, R _sBe direct radiation; R _dBe scattered radiation; R _rBe reflected radiation.Wherein:

$R_s={\left(\frac{1}{\mathrm{\&rho;}}\right)}^2{I}_0\cos {\mathrm{\&theta;\&tau;}}_b=R_0{\mathrm{\&tau;}}_b\cos \mathrm{\&theta;}$

Wherein, I ₀Be solar constant, get 1367 (Wm ^-2); R ₀Solar radiation for the straight incident of atmospheric envelope roof pendant; θ is the solar zenith angle on domatic; τ _bBe atmosphere direct projection transmissivity; (1/ ρ) ²Be called solar distance correction coefficient or Earth's orbit and correct the factor, about 10 ^-4During precision, calculate by fourier progression expanding method:

R _d＝R ₀τ _dcos ²β/2sina

In the formula, β is the gradient, and a is a sun altitude.

R _r＝rR ₀τ _rsin ²β/2sina

In the formula, r is the first reflectivity in slope, the face of land, can simplify being taken as 0.2; τ _rTransmissivity for reflected radiation.

8. the draught monitor system based on surface water thermoflux remote-sensing inversion according to claim 6 is characterized in that, described sensible heat flux computing module comprises:

Correcting module is used for providing degree of stability correction function Ψ based on similarity theory _m(x), Ψ _h(x), respectively to friction velocity u _*, aerodynamic resistance r _AhRevise:

$u_{*}=\frac{\mathrm{ku}\left(z\right)}{\ln \left(\frac{z-d}{z_{\mathrm{om}}}\right)-{\mathrm{\&Psi;}}_m\left(\frac{z-d}{L}\right)+{\mathrm{\&Psi;}}_m\left(\frac{z_{\mathrm{om}}}{L}\right)}$

$r_{\mathrm{ah}}=\frac{1}{{\mathrm{ku}}_{*}}\left[\ln \left(\frac{Z_2}{Z_1}\right)-{\mathrm{\&Psi;}}_h\left(\frac{Z_2}{L}\right)+{\mathrm{\&Psi;}}_h\left(\frac{Z_1}{L}\right)\right]$

Wherein, k is von Karman constant (getting 0.41, dimensionless), and z is the height of mean wind speed correspondence, z _OmBe face of land power roughness, Z ₂It is the reference altitude of measuring temperature; Z ₁=z _Oh, the heat delivered roughness; u _*Be friction velocity (ms ^-1), L is a Monin-Obukhov length,

In the formula, ρ _aBe atmospheric density; c _pIt is the specific heat capacity at constant pressure of air; T _sBe the face of land of remote sensing or the surface temperature of canopy; K is a von Karman constant; u _*It is friction velocity; G is an acceleration of gravity; H is the sensible heat flux;

The surface temperature computing module is used to calculate the surface temperature based under the DEM,

T _{S_dem}=T _s+ 0.0065 (h-h _Mean), wherein, h is a sea level elevation, h _MeanSea level on the average for survey region;

Pixel is chosen module, is used to choose extremely dried pixel and extremely wet pixel, with described T _{S_dem}Substitution formula dT=aT _{S_dem}+ b, wherein,

T _HotThe temperature of the extremely dried pixel of expression, T _ColdThe temperature of the extremely wet pixel of expression.

9. the draught monitor system based on surface water thermoflux remote-sensing inversion according to claim 8, it is characterized in that, pixel is chosen in the module, and the described criterion of choosing extremely dried pixel and extremely wet pixel is: the MSAVI about 0.1 reaches wherein combination correspondence " Hot " point of maximum LST; 0.8 about MSAVI and combination correspondence " Cold " point of wherein minimum LST.

10. the draught monitor system based on surface water thermoflux remote-sensing inversion according to claim 6, it is characterized in that, described regional lack of water Index for Calculation module, Regional Drought index based on RWSI, utilize the national standard of the damage caused by a drought grade classification of soil moisture, the convert arid ranking score grade standard of RWSI of each time period carries out the arid remote sensing monitoring in zone.

Images