CN106169058A - Pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information - Google Patents

Abstract

The invention discloses pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information, it is characterized in that, this method considers concordance and the seasonal effect in time series periodicity thereof of surface temperature spatial distribution, and with time series filtering, statistical model has been carried out effective combination.Carry out just estimating to pixel LST value under cloud initially with multichannel statistical model based on passive microwave remote sensing, in this course, NDVI is carried out ground mulching classification as classification foundation, abandon using existing ground surface type Remote Sensing Products to improve nicety of grading；Then estimated value based on statistical model is filled to LST time series as background value, it is considered to front and back the period is on the impact without the LST period, select mobile weighting filter to correct estimated value and obtain pixel LST reconstructed results under cloud.

Description

Pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information

Technical field

The present invention relates to pixel LST evaluation method technical field under MODIS cloud, particularly relate to a kind of based on microwave remote sensing with Pixel LST evaluation method under the cloud of space time information.

Background technology

Surface temperature (Land Surface Temperature, LST) is to weigh the crucial ginseng of earth surface Heat And Water Balance Number, the correlational study for scientific domains such as weather, the hydrology, geophysicses is significant.Since 20 century 70s, state Inside and outside scholar has carried out numerous studies with regard to how utilizing thermal infrared remote sensing to obtain surface temperature, it is proposed that including Split window algorithms Multiple inversion algorithm.Although thermal infrared remote sensing wide coverage, real-time, can quick obtaining Land surface emissivity and temperature Degree information, but it is easily affected by atmosphere vapour, it is impossible to through cloud layer, the round-the-clock monitoring of surface temperature therefore can not be carried out.When Before the MODIS surface temperature product that is widely used, partial image missing value rate is up to 60%, has had a strong impact on product and has used effect Rate, it is achieved the reconstruction to pixel surface temperatures such as cloud coverings is significant.And passive microwave remote sensing is little by atmospheric interference, can Penetrate cloud layer and obtain surface radiation information, and there is round-the-clock, multipolarization feature, Surface Temperature Retrieval has the excellent of uniqueness More property, provides possibility for surface temperature products such as reconstruct MODIS.Inquire into the dependency of the bright temperature in earth's surface and surface temperature, set up Therebetween relational model inverting surface temperature becomes the focus and emphasis of current research, and the microwave being widely used at present passes Sensor mainly has three kinds: SMMR (the Scanning Multichannel Microwave being mounted on Nimbus-7 satellite Radiometer) sensor, the SSM/I that is mounted in the serial polar-orbiting satellite of U.S. national defense meteorological satellite plan (DMSP) Terra and Aqua of (Special Sensor Microwave Imager) sensor and earth observing system (EOS) defends AMSR_E (the Advanced Microwave Scanning Radiometer-EOS) sensor that star is carried.Wherein, SMMR The sensor time in orbit is 1978～1987, and the SSM/I time in orbit is 1987～2005, AMSR_E sensor Serviced so far from 2002, and the spatial resolution that its spatial resolution is higher than SMMR and SSM/I, the bright temperature product of this sensor It it is the preference data source being currently used in inverting surface temperature；

The algorithm being currently based on passive microwave remote sensing inverting surface temperature mainly includes three classes, be respectively statistical model method, Physical model method and neural network algorithm.Wherein, physical model method, with radiation transfer equation as theoretical basis, has clear and definite thing Reason meaning, is not limited by space-time, but model parameter is many and complicated, mostly uses the method for condition hypothesis to simplify, and impact is anti- Drill precision；And neural network algorithm can be with Parallel implementation nonlinear problem, it is not necessary to accurately portray inversion equation, it is to avoid right The simulation of uncertain factor in physical parameter, but this algorithm is to weaken surface temperature ill-posed inversion from mathematics aspect to ask Topic, lacks actual physics meaning；Statistical model rule is to utilize internal relation between bright temperature and surface temperature, by a large amount of numbers According to returning Rule Summary, algorithm mechanism is simple, it is simple to realize, it has also become current application is most commonly used based on passive microwave remote sensing The method of inverting surface temperature.According to the bright temperature number of channels of participation matching, statistical model method can be divided into again one channel model And multi-channel model, wherein the efficiency of inverse process of multi-channel model is far better than single pass efficiency of inverse process, and its algorithm idea is: base Realizing pixel classification in ground surface type, category sets up the regression model of LST and the bright temperature in earth's surface, according to the ground without LST value pixel Table Properties selects suitable regression model estimation surface temperature value.But, this algorithm is in actual applications with existing ground surface type On the basis of product (such as MOD12), inversion accuracy is low, and mean error is at about 3K；Inversion result is substantially by environmental condition and meteorology Condition limits, and missing value rate is the highest, and inversion error is the biggest；Inverting region is in units of pixel, it is impossible to realize the earth's surface temperature in big region Degree is rebuild, and the practicality of algorithm is low.

Summary of the invention

The technical problem existed based on background technology, the present invention proposes a kind of cloud based on microwave remote sensing Yu space time information Lower pixel LST evaluation method.

Pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information that the present invention proposes, including as follows Step:

S1, ground mulching classification: the pre-places such as Landsat/TM remotely-sensed data collection corrects via radiation, atmospheric correction, resampling After reason, utilizing ENVI its NDVI of band math function calculating, formula is as follows:

NDVI=(b4-b3)/(b4+b3) 1

In formula, b3, b4 are respectively the reflectance of third channel and fourth lane.

According to NDVI value, ground mulching is divided into 6 big classes, respectively: NDVI ＜ 0；0≤NDVI ＜ 0.1；0.1≤NDVI ＜ 0.3；0.3≤NDVI ＜ 0.5；0.5≤NDVI ＜ 0.7；NDVI≥0.7；

S2, pixel is split: select 18.7,23.8,36.5 to participate in Surface Temperature Retrieval, profit with 8 passages of 89GHz (H/V) With ENVI software, 18.7,23.8 and 36.5GHz (H/V) data are carried out pixel segmentation so that it is consistent with 89GHz Pixel size, In LST process of reconstruction, coordinates restriction is utilized to realize the spatial match of the bright temperature in earth's surface and surface temperature；

S3, MODIS LST Yu AMSR_E brightness temperature regression analysis: the spoke that microwave radiometer receives under given frequency range Penetrate bright temperature to be expressed as:

T_f=τ_fε_fT_s+Δ 2

In formula, T_fFor bright temperature, T_sFor surface temperature, ε_fFor surface radiation rate, τ_fFor atmospheric transmittance, Δ is that other improve , including air uplink radiation and downlink radiation；Above formula is reduced to:

T_f=ε_fT_s 3

S4, sets up the multivariate regression models of surface temperature and the bright temperature in earth's surface: MODIS surface temperature and the bright temperature in AMSR_E earth's surface Multivariate regression models as shown in Equation 3:

${\mathrm{LST}}_{n\×1}^{MODIS}=T_{n\×8}^{AMSR\_E}\×A_{8\×1}+B_{n\×1}---4$

$\begin{array}{ccc}LST=\left[\begin{array}{c}{\mathrm{lst}}_1\\ {\mathrm{lst}}_2\\ \mathrm{...}\\ {\mathrm{lst}}_{n-1}\\ {\mathrm{lst}}_n\end{array}\right]& A=\left[\begin{array}{c}{a}_1\\ a_2\\ \mathrm{...}\\ a_7\\ a_8\end{array}\right]& B=\left[\begin{array}{c}{b}_1\\ b_2\\ \mathrm{...}\\ b_{n-1}\\ b_n\end{array}\right]\end{array}$

$T=\left[\begin{array}{cccccccc}{t}_1^{18.7H}& t_1^{18.7V}& t_1^{23.8H}& t_1^{23.8V}& t_1^{36.5H}& t_1^{36.5V}& t_1^{89H}& t_1^{89V}\\ t_2^{18.7H}& t_2^{18.7V}& t_2^{23.8H}& t_2^{23.8V}& t_2^{36.5H}& t_2^{36.5V}& t_2^{89H}& t_2^{89V}\\ \mathrm{...}& \mathrm{...}& \mathrm{....}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ t_{n-1}^{18.7H}& t_{n-1}^{18.7V}& t_{n-1}^{23.8H}& t_{n-1}^{23.8V}& t_{n-1}^{36.5H}& t_{n-1}^{36.5V}& t_{n-1}^{89H}& t_{n-1}^{89V}\\ t_n^{18.7H}& t_n^{18.7V}& t_n^{23.8H}& t_n^{23.8V}& t_n^{36.5H}& t_n^{36.5V}& t_n^{89H}& t_n^{89V}\end{array}\right]$

In formula, LST is MODIS LST column vector, and T is the AMSR_E same period each passage earth's surface bright temperature Mean Matrix, and A, B are for treating Estimating parameter, n is the pixel quantity that MODIS participates in inverting；

S5, estimated value based on weighted moving average filtering corrects: setting independent variable x and do equidistant observation with step-length h, correspondence is seen Survey result y, it may be assumed that

x_i=x₀+ih 5

Table 1 observation sequence explanation

Formula 5 is converted, obtains:

$t=\frac{x-x_i}{h}---6$

Then:

Table 2 observation sequence deforms

Structure smoothing formula is as follows:

y_i'_+t=A_o+A₁t+A₂t²+...+A_mt^m 7

The coefficient of smooth polynomial is determined by the principle of least square, it may be assumed that

$min\left(\mathrm{\δ}\right)=min\[\underset{t}{\Σ}{(y_{i+t}^{\′}-y_{i+t})}^2\]---8$

That is:

$min\left(\mathrm{\δ}\right)=min\underset{t}{\Σ}{\left(\left(A_0+A_{1}t+A_2{t}^2+\mathrm{...}+A_m{t}^m\right)-y_{i+t}\right)}^2---9$

In formula, t takes 2n+1 integer value near i (i.e. smooth count), i.e. t=-n ,-n+1 ... 0 ... n-1, n, And should ensure that m ＜ 2n+1 ＜ N,

According toSolve each coefficient A_i, according to the principle of least square:

$\left\{\begin{array}{c}\frac{\∂\mathrm{\δ}}{\∂A_0}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)=0\\ \frac{\∂\mathrm{\δ}}{\∂A_1}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)t=0\\ \frac{\∂\mathrm{\δ}}{\∂A_2}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)t^2=0\end{array}\right.---10;$

By analyzing, when using weight secondary smothing filtering methods to MODIS surface temperature containing initial estimate at 5 Between sequence be reconstructed, solve weights coefficient according to formula (10), be applied in MODIS LST time series, obtain public affairs Formula (11):

$\left\{\begin{array}{c}\begin{array}{cc}y\left(i\right)=(-3\×y_0(i-2)+12\×y_0(i-1)+17\×y_0(i)+12\×y_0(i+1)-3\×y_0(i+2\left)\right)/35& (3\≤i\≤44)\end{array}\\ y\left(1\right)=\left(31\×y_0\left(1\right)+9\×y_0\left(2\right)-3\×y_0\left(3\right)-5\×y_0\left(4\right)+3\×y_0\left(5\right)\right)/35\\ y\left(2\right)=\left(9\×y_0\left(1\right)+13\×y_0\left(2\right)+12\×y_0\left(3\right)+6\×y_0\left(4\right)-5\×y_0\left(5\right)\right)/35\\ y\left(45\right)=\left(-5\×y_0\left(42\right)+6\×y_0\left(43\right)+12\×y_0\left(44\right)+13\×y_0\left(45\right)+9\×y_0\left(46\right)\right)/35\\ y\left(46\right)=\left(3\×y_0\left(42\right)-5\×y_0\left(43\right)-3\×y_0\left(44\right)+9\×y_0\left(45\right)+31\×y_0\left(46\right)\right)/35\end{array}\right.---11$

In formula, y (i) is the i-th phase LST reconstruction value, y₀I () is for containing the i-th phase temperature value being worth according to a preliminary estimate；

S6, LST reconstruction precision is evaluated: calculating each region LST and rebuild equal error, formula is as follows:

$ME=\underset{i=1}{\overset{n}{\Σ}}|P_i-O_i|/n---12$

In formula, ME is equal error, P_iIt is_iIndividual pixel LST reconstructed value, O_iFor original value, n adds up pixel amount.

Preferably, in 1 formula in described S1, b3, b4 are respectively the reflectance of third channel and fourth lane.

Preferably, in 10 formulas in described S5, the numerical value of employing is, m=2,2n+1=5, t=-2 ,-1,0,1,2；Formula 11 is the computing formula for MODIS LST 8Day product one annual data.

The present invention stating, the algorithm of passive microwave remote sensing inverting surface temperature mainly includes three classes, is statistical model respectively Method, physical model method and neural network algorithm.Wherein, physical model method, with radiation transfer equation as theoretical basis, has clearly Physical significance, do not limited by space-time, but model parameter be many and complicated, mostly use the method for condition hypothesis to simplify, shadow Ring inversion accuracy；And neural network algorithm can be with Parallel implementation nonlinear problem, it is not necessary to accurately portray inversion equation, it is to avoid To the simulation of uncertain factor in physical parameter, but this algorithm simply weakens surface temperature ill-posed inversion from mathematics aspect Problem, lacks actual physics meaning；Statistical model rule is to utilize internal relation between bright temperature and surface temperature, by a large amount of Data regression Rule Summary, algorithm mechanism is simple, it is simple to realize, it has also become current application is most commonly used distant based on passive microwave The method of sense inverting surface temperature.Statistical model method is divided into one channel model and multi-channel model, and its algorithm idea is based on ground Table type realizes pixel classification, and category sets up the surface temperature value of LST and the regression model estimation disappearance of the bright temperature in earth's surface.So And, this algorithm is in actual applications on the basis of existing ground surface type product (such as MOD12), and inversion accuracy is low, and mean error exists About 3K；Inversion result is substantially limited by environmental condition and meteorological condition, and missing value rate is the highest, and inversion error is the biggest；Inverting region In units of pixel, area the biggest error in study area is the biggest, it is impossible to the surface temperature realizing big region is rebuild, the practicality of algorithm Low.Present invention introduces time series filtering, before and after utilization, the impact without the LST period is corrected statistical model estimated value by the period, logical Cross prolongation information and reduce the inversion result dependence to spatial information, improve the space relevance grade of algorithm itself.Current open It is relatively big that the global seismic issued covers Satellite Product error, when carrying out ground mulching classification easily by different attribute or attribute phase Difference thing mistake significantly is divided into a class, affects the matching of model, reduces inversion accuracy.The present invention proposes to utilize NDVI, and (normalization is planted By index) as the benchmark of terrain classification.The means of vegetation index quantification describe vegetation coverage, different types of ground objects NDVI have higher can discrimination, be widely used in Land_use change as a kind of remote sensing and cover detection, vegetative coverage The aspects such as density evaluation, crop identification and crop forecast.Multiple, such as simple vegetation for explaining that the vegetation index of remotely-sensed data has Index (LCI), ratio vegetation index (RVI), normalized differential vegetation index (NDVI), green degree vegetation index (GVI) and soil regulation Vegetation indexs (SAVI) etc., most common of which is normalized differential vegetation index (NDVI), and the present invention is using NDVI as terrain classification Reference standard carry out ground mulching classification, set up based on underlying surface type, portray surface temperature temperature bright to earth's surface relevant close The multivariate regression models of system, it is achieved the first estimation to the surface temperature value without LST pixel；Consider that the periodicity of surface temperature is special Levying, introduce time series filtering, before and after utilization, the estimated value of statistical model is corrected in the impact without the LST period by the period, improves anti- Drill precision.

Accompanying drawing explanation

Fig. 1 is the LST of pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information that the present invention proposes Rebuild flow chart；

Fig. 2 is the LST of pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information that the present invention proposes Time series fitting of a polynomial exemplary plot；

Fig. 3 is the base of pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information that the present invention proposes LST reconstructed results comparison diagram in the secondary weighted filtering algorithm of different time window.

Detailed description of the invention

Below in conjunction with specific embodiment, the present invention is explained orally further.

With reference to Fig. 1-3, pixel LST estimation side under a kind of cloud based on microwave remote sensing with space time information that the present invention proposes Method, comprises the steps:

After the pretreatment such as S1, Landsat/TM remotely-sensed data collection corrects via radiation, atmospheric correction, resampling, utilize ENVI Its NDVI of band math function calculating, formula is as follows:

NDVI=(b4-b3)/(b4+b3) 1

In formula, b3, b4 are respectively the reflectance of third channel and fourth lane.

According to NDVI value, ground mulching is divided into 6 big classes by the present invention, respectively: NDVI ＜ 0；0≤NDVI ＜ 0.1；0.1 ≤ NDVI ＜ 0.3；0.3≤NDVI ＜ 0.5；0.5≤NDVI ＜ 0.7；NDVI≥0.7；

S2, pixel is split: studies and shows that 37GHz is less by the scatterings such as steam, cloud, rain and inhalation effects, suitable for anti- Drill surface temperature；Additionally, 37GHz Yu 22GHz combination of channels is conducive to eliminating atmosphere vapour impact, 37GHz Yu 19GHz passage group Close and be conducive to eliminating soil moisture impact, and 89GHz with MODIS surface temperature dependency is the highest.Therefore, the present invention selects 18.7,23.8,36.5 Surface Temperature Retrieval is participated in 8 passages of 89GHz (H/V).But, the AMSR_E earth's surface of different channel is bright Temperature data spatial resolution differs greatly, and for convenience of calculating, utilizes ENVI software to 18.7,23.8 and 36.5GHz (H/V) data Carry out pixel segmentation so that it is consistent with 89GHz Pixel size, in LST process of reconstruction, utilize coordinates restriction to realize the bright temperature in earth's surface Spatial match with surface temperature；

S3, MODIS LST Yu AMSR_E brightness temperature regression analysis: according to radiation transfer equation, microwave radiometer is being given The radiation brightness determining to receive under frequency range can be expressed as:

T_f=τ_fε_fT_s+Δ 2

In formula, T_fFor bright temperature, T_sFor surface temperature, ε_fFor surface radiation rate, τ_fFor atmospheric transmittance, Δ is that other improve , including air uplink radiation and downlink radiation etc..Ignoring atmospheric effect, above formula can be reduced to:

T_f=ε_fT_s 3

It is known that how many dielectric constants of major part atural object is mainly determined by its contained humidity.Utilize passive microwave bright Temperature data acquisition land table temperature is exactly in the case of surface radiation rate the unknown, by the most like-polarized bright temperature group of different frequency Close, set up the empirical relation between each Channels Brightness Temperature and surface temperature；

S4, surface temperature and the multivariate regression models of the bright temperature in earth's surface: MODIS surface temperature is many with the bright temperature in AMSR_E earth's surface Unit regression model as shown in Equation 3:

${\mathrm{LST}}_{n\×1}^{MODIS}=T_{n\×8}^{AMSR\_E}\×A_{8\×1}+B_{n\×1}---4$

$\begin{array}{ccc}LST=\left[\begin{array}{c}{\mathrm{lst}}_1\\ {\mathrm{lst}}_2\\ \mathrm{...}\\ {\mathrm{lst}}_{n-1}\\ {\mathrm{lst}}_n\end{array}\right]& A=\left[\begin{array}{c}{a}_1\\ a_2\\ \mathrm{...}\\ a_7\\ a_8\end{array}\right]& B=\left[\begin{array}{c}{b}_1\\ b_2\\ \mathrm{...}\\ b_{n-1}\\ b_n\end{array}\right]\end{array}$

$T=\left[\begin{array}{cccccccc}{t}_1^{18.7H}& t_1^{18.7V}& t_1^{23.8H}& t_1^{23.8V}& t_1^{36.5H}& t_1^{36.5V}& t_1^{89H}& t_1^{89V}\\ t_2^{18.7H}& t_2^{18.7V}& t_2^{23.8H}& t_2^{23.8V}& t_2^{36.5H}& t_2^{36.5V}& t_2^{89H}& t_2^{89V}\\ \mathrm{...}& \mathrm{...}& \mathrm{....}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ t_{n-1}^{18.7H}& t_{n-1}^{18.7V}& t_{n-1}^{23.8H}& t_{n-1}^{23.8V}& t_{n-1}^{36.5H}& t_{n-1}^{36.5V}& t_{n-1}^{89H}& t_{n-1}^{89V}\\ t_n^{18.7H}& t_n^{18.7V}& t_n^{23.8H}& t_n^{23.8V}& t_n^{36.5H}& t_n^{36.5V}& t_n^{89H}& t_n^{89V}\end{array}\right]$

In formula, LST is MODIS LST column vector, and T is the AMSR_E same period each passage earth's surface bright temperature Mean Matrix, and A, B are for treating Estimating parameter, n is the pixel quantity that MODIS participates in inverting；

S5, estimated value based on weighted moving average filtering corrects: it is known that surface temperature has cyclically-varying spy Levying, introduce time series filtering, in available short time window, the estimated value of previous step is entered by the concordance of surface temperature change Row corrects, i.e. prolongation information reduces estimation difference.In general, areal, the time closer to, surface temperature value gets over phase Seemingly.Based on this premise, the present invention uses weighted moving average filter method to containing the LST time series that is worth according to a preliminary estimate in addition Process；The basic thought of weighting is the maximum weight of centre data in average area, and the more data weights at off-center are the least, Reduce the smoothing effect to actual signal itself the most to a certain extent.If independent variable x does equidistant observation with step-length h, Corresponding observed result y, it may be assumed that

x_i=x₀+ih 5

Table 1 observation sequence explanation

Formula 5 is converted, obtains:

$t=\frac{x-x_i}{h}---6$

Then:

Table 2 observation sequence deforms

Structure smoothing formula is as follows:

y_i'_+t=A_o+A₁t+A₂t²+...+A_mt^m 7

The coefficient of smooth polynomial is determined by the principle of least square, it may be assumed that

$min\left(\mathrm{\δ}\right)=min\[\underset{t}{\Σ}{(y_{i+t}^{\′}-y_{i+t})}^2\]---8$

That is:

$min\left(\mathrm{\δ}\right)=min\underset{t}{\Σ}{\left(\left(A_0+A_{1}t+A_2{t}^2+\mathrm{...}+A_m{t}^m\right)-y_{i+t}\right)}^2---9$

In formula, t takes 2n+1 integer value near i (i.e. smooth count), i.e. t=-n ,-n+1 ... 0 ... n-1, n, And should ensure that m ＜ 2n+1 ＜ N,

According toSolve each coefficient A_i, according to the principle of least square:

$\left\{\begin{array}{c}\frac{\∂\mathrm{\δ}}{\∂A_0}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)=0\\ \frac{\∂\mathrm{\δ}}{\∂A_1}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)t=0\\ \frac{\∂\mathrm{\δ}}{\∂A_2}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)t^2=0\end{array}\right.---10$

Analyze earth's surface temperature-time sequence character, find temperature value and to have higher second nonlinear between the time relevant Feature, as shown in figure a (LST time series example).Additionally, smooth window is the biggest, smooth effect is the best, but distorted signals is also got over Greatly.According to analysis, the MODIS surface temperature time series containing initial estimate is entered by 5 smooth window of the most final employing Row filtering reconstruct.Solve weights coefficient according to formula 10, be applied in MODIS LST time series, with an annual data be Example, obtains formula (11).

$\left\{\begin{array}{c}\begin{array}{cc}y\left(i\right)=(-3\×y_0(i-2)+12\×y_0(i-1)+17\×y_0(i)+12\×y_0(i+1)-3\×y_0(i+2\left)\right)/35& (3\≤i\≤44)\end{array}\\ y\left(1\right)=\left(31\×y_0\left(1\right)+9\×y_0\left(2\right)-3\×y_0\left(3\right)-5\×y_0\left(4\right)+3\×y_0\left(5\right)\right)/35\\ y\left(2\right)=\left(9\×y_0\left(1\right)+13\×y_0\left(2\right)+12\×y_0\left(3\right)+6\×y_0\left(4\right)-5\×y_0\left(5\right)\right)/35\\ y\left(45\right)=\left(-5\×y_0\left(42\right)+6\×y_0\left(43\right)+12\×y_0\left(44\right)+13\×y_0\left(45\right)+9\×y_0\left(46\right)\right)/35\\ y\left(46\right)=\left(3\×y_0\left(42\right)-5\×y_0\left(43\right)-3\×y_0\left(44\right)+9\×y_0\left(45\right)+31\×y_0\left(46\right)\right)/35\end{array}\right.---11$

In formula, y (i) is the i-th phase LST reconstruction value, y₀I () is for containing the i-th phase temperature value being worth according to a preliminary estimate.

S6, LST reconstruction precision is evaluated: successively the value pixel that has of every scape image is regarded as pixel under cloud, uses this patent to carry Its surface temperature value is rebuild by the algorithm for reconstructing gone out, and contrasts with original value, and statistics absolute error is at 0～1K, 0～2K respectively Pixel scale with 0～3K, calculates each region LST and rebuilds equal error, and formula is as follows.

$ME=\underset{i=1}{\overset{n}{\Σ}}|P_i-O_i|/n---12$

In formula, ME is equal error, P_iIt is_iIndividual pixel LST reconstructed value, O_iFor original value, n adds up pixel amount.

Owing to this algorithm adds time dimension information, reduce the inversion result degree of dependence to spatial information, rebuild Process is significantly reduced by the constraint of environment and meteorological condition, rebuilds efficiency height, and the space that reconstructed results meets surface temperature is consistent Property.

With MOD11A2 surface temperature product as object of study, with Xinjiang, Qinghai, Sichuan, Yunnan, Henan, Anhui, Hubei, Pixel under the MODIS cloud in these regions, as study area, is carried out rebuilding research, and contrasts base by totally 9 provinces such as Hunan, Jiangxi In the inversion result of the spatial domain multichannel statistical model of MODIS ground mulching product MOD12Q1, reconstruction precision significantly improves: The average inversion accuracy of statistical model based on MOD12Q1 classification is 5.14K, and averagely rebuilds essence based on time-space domain unified algorithm Degree is 1.24K, improves 76 percentage points；Therefore, compared to conventional research, this algorithm has upper zone universality and weight Build precision, pixel LST under the cloud in big region can be realized well and rebuild, practical.

The above, the only present invention preferably detailed description of the invention, but protection scope of the present invention is not limited thereto, Any those familiar with the art in the technical scope that the invention discloses, according to technical scheme and Inventive concept equivalent or change in addition, all should contain within protection scope of the present invention.

Claims (3)

1. a pixel LST evaluation method under cloud based on microwave remote sensing and space time information, comprises the steps:

S1, ground mulching classification: after the pretreatment such as Landsat/TM remotely-sensed data collection corrects via radiation, atmospheric correction, resampling, Utilizing ENVI its NDVI of band math function calculating, formula is as follows:

NDVI=(b4-b3)/(b4+b3) 1

In formula, b3, b4 are respectively the reflectance of third channel and fourth lane.

According to NDVI value, ground mulching is divided into 6 big classes, respectively: NDVI ＜ 0；0≤NDVI ＜ 0.1；0.1≤NDVI ＜ 0.3；0.3≤NDVI ＜ 0.5；0.5≤NDVI ＜ 0.7；NDVI≥0.7；

S2, pixel is split: select in AMSR E data set 18.7,23.8,36.5 to participate in earth's surface with 8 passages of 89GHz (H/V) Temperature retrieval, utilizes ENVI software that 18.7,23.8 and 36.5GHz (H/V) data are carried out pixel segmentation so that it is with 89GHz picture Unit is in the same size, in LST process of reconstruction, utilizes coordinates restriction to realize the spatial match of the bright temperature in earth's surface and surface temperature；

S3, MODIS LST Yu AMSR_E brightness temperature regression analysis: the radiation that microwave radiometer receives under given frequency range is bright Temperature is expressed as:

T_f=τ_fε_fT_s+Δ 2

In formula, T_fFor bright temperature, T_sFor surface temperature, ε_fFor surface radiation rate, τ_fFor atmospheric transmittance, Δ is that other improve item, bag Include air uplink radiation and downlink radiation；Above formula is reduced to:

T_f=ε_fT_s 3

S4, sets up the multivariate regression models of surface temperature and the bright temperature in earth's surface: MODIS surface temperature is many with the bright temperature in AMSR_E earth's surface Unit regression model as shown in Equation 3:

${\mathrm{LST}}_{n\×1}^{MODIS}=T_n^{AMSR\_E}\×A_{8\×1}+B_{n\×1}---4$

$\begin{array}{ccc}LST=\left[\begin{array}{c}ls{t}_1\\ {\mathrm{lst}}_2\\ ...\\ ls{t}_{n-1}\\ {\mathrm{lst}}_n\end{array}\right]& A=\left[\begin{array}{c}{a}_1\\ a_2\\ ...\\ a_7\\ a_8\end{array}\right]& B=\left[\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_{n-1}\\ b_n\end{array}\right]\end{array}$

$T=\left[\begin{array}{cccccccc}{t}_1^{18.7H}& t_1^{18.7V}& t_1^{23.8H}& t_1^{23.8V}& t_1^{36.5H}& t_1^{36.5V}& t_1^{89H}& t_1^{89V}\\ t_2^{18.7H}& t_2^{18.7V}& t_2^{23.8H}& t_2^{23.8V}& t_2^{36.5H}& t_2^{36.5V}& t_2^{89H}& t_2^{89V}\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}& \mathrm{...}\\ t_{n-1}^{18.7H}& t_{n-1}^{18.7V}& t_{n-1}^{23.8H}& t_{n-1}^{23.8V}& t_{n-1}^{36.5H}& t_{n-1}^{36.5V}& t_{n-1}^{89H}& t_{n-1}^{89V}\\ t_n^{18.7H}& t_n^{18.7V}& t_n^{23.8H}& t_n^{23.8V}& t_n^{36.5H}& t_n^{36.5V}& t_n^{89H}& t_n^{89V}\end{array}\right]$

In formula, LST is MODIS LST column vector, and T is the AMSR_E same period each passage earth's surface bright temperature Mean Matrix, and A, B are ginseng to be estimated Number, n is the pixel quantity that MODIS participates in inverting；

S5, estimated value based on weighted moving average filtering corrects: sets independent variable x and does equidistant observation, corresponding observation knot with step-length h Really y, it may be assumed that

x_i=x₀+ih 5

Table 1 observation sequence explanation

Formula 5 is converted, obtains:

$t=\frac{x-x_i}{h}---6$

Then:

Table 2 observation sequence deforms

Structure smoothing formula is as follows:

y′_i+t=A_o+A₁t+A₂t²+...+A_mt^m 7

The coefficient of smooth polynomial is determined by the principle of least square, it may be assumed that

$min\left(\mathrm{\δ}\right)=\min \[\underset{t}{\Σ}{(y_{i+t}^{\′}-y_{i+t})}^2\]---8$

That is:

$min\left(\mathrm{\δ}\right)=\min \underset{t}{\Σ}{\left(\right(A_0+A_{1}t+A_2{t}^2+...+A_m{t}^m)-y_{i+t})}^2---9$

In formula, t takes 2n+1 integer value near i (i.e. smooth count), i.e. t=-n ,-n+1 ... 0 ... n-1, n, and should Ensure m ＜ 2n+1 ＜ N,

According toSolve each coefficient A_i, according to the principle of least square:

$\left\{\begin{array}{c}\frac{\∂\mathrm{\δ}}{\∂A_0}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)=0\\ \frac{\∂\mathrm{\δ}}{\∂A_1}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)t=0\\ \frac{\∂\mathrm{\δ}}{\∂A_2}=\underset{t=-2}{\overset{2}{\Σ}}(y_{i+t}-A_0-A_{1}t-A_2{t}^2)t^2=0\end{array}\right.---10;$

By analyzing, use 5 weighting secondary smothing filtering methods to the MODIS surface temperature time sequence containing initial estimate Row are reconstructed, and solve weights coefficient according to formula 10, and are applied to MODIS LST time series (with 8Day data product Product, as a example by the term) in, obtain formula 11:

$\left\{\begin{array}{c}\begin{array}{cc}y\left(i\right)=(-3\×y_0(i-2)+12\×y_0(i-1)+17\×y_0(i)+12\×y_0(i+1)-3\×y_0(i+2\left)\right)/35& (3\≤i\≤44)\end{array}\\ y\left(1\right)=(31\×y_0(1)+9\×y_0(2)-3\×y_0(3)-5\×y_0(4)+3\×y_0(5\left)\right)/35\\ y\left(2\right)=(9\×y_0(1)+13\×y_0(2)+12\×y_0(3)+6\×y_0(4)-5\×y_0(5\left)\right)/35\\ y\left(45\right)=(-5\×y_0(42)+6\×y_0(43)+12\×y_0(44)+13\×y_0(45)+9\×y_0(46\left)\right)/35\\ y\left(46\right)=(3\×y_0(42)-5\×y_0(43)-3\×y_0(44)+9\×y_0(45)+31\×y_0(46\left)\right)/35\end{array}\right.---11$

In formula, y (i) is_iPhase LST reconstruction value, y₀I () is for containing the be worth according to a preliminary estimate_iPhase temperature value；

S6, LST reconstruction precision is evaluated: calculating each region LST and rebuild equal error, formula is as follows:

$ME=\underset{i=1}{\overset{n}{\Σ}}\left|P_i-O_i\right|/n---12$

In formula, ME is equal error, P_iIt is_iIndividual pixel LST reconstructed value, O_iFor original value, n adds up pixel amount.

Pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information the most according to claim 1, it is special Levying and be, in 1 formula in described S1, b3, b4 are respectively the reflectance of third channel and fourth lane.

Pixel LST evaluation method under a kind of cloud based on microwave remote sensing with space time information the most according to claim 1, it is special Levying and be, in 10 formulas in described S5, the numerical value of employing is, m=2,2n+1=5, t=-2 ,-1,0,1,2；Formula 11 be for The computing formula of MODIS LST8Day product one term data.