CN102103076B - Surface albedo inversion method and system - Google Patents

Abstract

The invention discloses a surface albedo inversion method, which comprises the following steps of: 1, acquiring multi-angle surface two-direction albedo data by using each image element on remotely sensed data of an observation satellite; 2, selecting corresponding component wave spectrum data from a prior wave spectrum knowledge base according to the surface type of each image element in the remotely sensed data; 3, integrating the component wave spectrum data to the corresponding wave band according to wave band setting of different sensors; 4, reading observation geometrical data from the remotely sensed data; 5, constructing a linear equation set; 6, calculating a black hemisphere wave spectrum albedo, a white hemisphere wave spectrum albedo and an actual surface wave spectrum albedo; and 7, calculating a black hemisphere wide wave band albedo, a white hemisphere wide wave band albedo and an actual surface wide wave band albedo in any wave band. The method solves the problems that measured information quantity is insufficient, the error when a narrow wave band is converted into a wide wave band in the conventional algorithm is large, and the utilization rate of the remotely sensed data is low.

Description

Surface albedo inversion method and system

Technical field

The present invention relates to the remote sensing technology field, particularly a kind of surface albedo inversion method and system.

Background technology

In the parameter of numerous faces of land; Face of land bidirectional reflectance distribution function (BidirectionalReflectance Distribution Function, BRDF)/the accurate remote-sensing inversion of albedo usually is the condition precedent of estimating that accurately other top parameters such as vegetation parameter, soil utilization and soil cover in the narrow sense.Therefore, as the key variables that influence the terrestrial climate system, surface albedo is an important parameter in numerical value climate model and the face of land energy-balance equation.Accurately calculate albedo of underlying surface and can disclose the local inherent mechanism that changes with regional climate, can better carry out parametrization, in the raising, the level of long-term climatic prediction ecological models.Though face of land BRDF/ albedo remote-sensing inversion has very important significance, in BRDF/ albedo refutation process, it is present to face main issue table:

1) face of land two of the multi-angle that obtains of sensor is low to the reflectivity data utilization ratio, how to work in coordination with that to utilize RS data to have to be strengthened

The BRDF inverting needs the information of a plurality of angles, and albedo comprises the reflected energy of all directions of hemisphere space, the face of land as the ratio of all shortwave reflected energies with projectile energy.The reason that in most cases reflects owing to direction; Different directions earth surface reflection rate is different; Different surface is formed, and is also variant at the reflectivity of all directions, so accurately carry out the face of land two of multi-angle that face of land BRDF/ albedo remote-sensing inversion needs all hemisphere directions to reflectivity data.And this realizes than difficult for remote sensing technology.Present most of sensor is only observed in a direction.Even for big visual field sensor, under the little situation of supposition face of land short term variations, utilize the repeated measures of the multidate of same face of land target can obtain the information of a plurality of angles.But the required cycle is longer, and (Moderate-Resolution ImagingSpectroradiometer MODIS) needs time of two weeks to generate the product of BRDF/ albedo like Moderate Imaging Spectroradiomete.And, be difficult to guarantee the unchangeability on the face of land because time span is big.Even so, directly carry out the demand that albedo is extracted even all angles are all used also can't satisfy.The multi-source data fusion has to be strengthened, makes full use of multispectral information, and multi-angle sensor information is combined with multispectral remote sensing information.

2) spectral information utilizes form single, is combined with to be strengthened with data of multiple angles

Surface albedo still is an all band spectrum albedo sum, but the wave band data on the moment sensor observation face of land is limited, and sensor can't cover all band, and therefore going to calculate full wave albedo with limited wave band has bigger error.When all band albedo is obtained, scanning radiometer AVHRR the earliest.Two available band can only be provided, four wave bands of MSS, (Thematic Mapper TM) can bring up to 6 wave bands to TM, even MODIS also can only provide 20 wave bands at present.Present BRDF/ albedo (ALBEDO) model mainly is to sensor wave band to be set; Also be to implement respectively during inverting; And do not make full use of the advantage of multispectral information; And regional also the giving with the gap that responds of band setting brought difficulty between the different sensors when RS data is worked in coordination with inverting, and feasible data from the different platform sensor can not combine effectively, how effectively to fully utilize the multispectral remotely-sensed data of multi-angle and awaits to strengthen; Because both combinations can be us and provide more about the information of vegetation itself, also can in inverting, reduce the requirement of the face of land two of multi-angle to the reflectivity data amount.

3) two of the face of land obtain to reflection

Though have at present many models describe the face of land two to reflection characteristic.But because the finiteness of observed ray is difficult to accurately carry out the face of land two obtaining to reflection at present.Some is too complicated for present physical model, and parametric inversion is difficulty comparatively.The face of land two that also needs multi-angle as semiempirical nuclear driving model is to reflectivity data, and carries out the parameter that interative computation obtains model, the big and less stable of operand.

Generally speaking, face of land BRDF/ albedo remote-sensing inversion needs the support of multi-angle information, can only pass through realizations such as multi-angle/big visual field sensor such as MODIS, POLDER, MISR at present.The BRDF/ albedo remote-sensing inversion that single-sensor utilizes multi-angle information to carry out often has the quantity of information deficiency; Inversion accuracy can not meet the demands; Particularly on the face of land two of multi-angle under the unreasonable condition of reflectivity data angle distribution; There is great risk in inverting, generates extreme wrong result sometimes; Owing to the limitation of sensor band setting, single sensor often can not cover all band simultaneously, asks for full wave surface albedo with extremely limited wave band and also has bigger error.When obtaining the broadband albedo, need the narrow wave band albedo that inverting obtains be transformed on the broadband; At present to the most algorithms that adopted given experience weight of studying of this problem; This process can be brought extra error again undoubtedly, and weight coefficient also changes with different surface type and different satellite, and also there is bigger difference in the result that different researchers provides; Region is very strong; Be difficult to promote, this otherness is also restricting the collaborative inverting of comprehensive utilization RS data simultaneously, can not better synthesis utilize limited data resource.

Summary of the invention

The technical matters that (one) will solve

The technical matters that the present invention will solve is how to solve the observation information quantity not sufficient, and narrow wave band is very big to broadband conversion error in the traditional algorithm, and the low problem of remotely-sensed data utilization factor.

(2) technical scheme

For solving the problems of the technologies described above, the invention provides a kind of surface albedo inversion method, may further comprise the steps:

S1: the face of land two that obtains multi-angle through each pixel on the remotely-sensed data of observation satellite is to reflectivity data;

S2: based on the face of land type of each pixel in the said remotely-sensed data; Select corresponding component spectral data from priori wave spectrum knowledge base; Said priori wave spectrum knowledge base is the database of the storage face of land type and the component spectral data of correspondence, and said component spectral data is continuous;

S3: the band setting to different sensors is integrated to corresponding wave band with said component spectral data;

S4: from said remotely-sensed data, read the observation geometric data, said observation geometric data comprises: the relative bearing of solar zenith angle, observation zenith angle and the sun and satellite;

S5: component spectral data and the said observation geometric data substitution linear equation of the face of land two after reflectivity data and process integral processing according to said multi-angle, make up system of linear equations, comprise karyonide number and kernel function in the said linear equation;

S6: solve the karyonide number of said system of linear equations through least square method,, calculate black hemisphere wave spectrum albedo, white hemisphere wave spectrum albedo and real surface wave spectrum albedo according to the karyonide number and the kernel function of said linear equation;

S7: according to the karyonide number of said linear equation and kernel function, and the wave spectrum of total descending radiant quantity distribute, calculate the black hemisphere broadband albedo in the random wave segment limit, white hemisphere broadband albedo and true broadband albedo.

Wherein, face of land type comprises described in the step S2: vegetation, soil, water body and snow.

Wherein, in step S3, carry out integration through following formula:

$R_i={\∫}_{{\mathrm{\λ}}_{\mathrm{si}}}^{{\mathrm{\λ}}_{\mathrm{ei}}}f\left(\mathrm{\λ}\right)\·R\left(\mathrm{\λ}\right)\mathrm{d\λ}$

Wherein, i representes i wave band, λ _SiAnd λ _EiBe respectively the initial wavelength and termination wavelength of i wave band, f (λ) is the wave spectrum response function of respective sensor, and R (λ) is continuous component spectral data, R _iBe component spectral data through integral processing.

Wherein, the linear equation described in the step S5 is:

R(θ _i，θ _v，φ，λ)＝c ₁·k ₁(λ)+c ₂·k ₂(λ)+c ₃·k ₃(θ _i，θ _v，φ，λ)+c ₄·k ₄(θ _i，θ _v，φ，λ)+c ₅·k ₅(θ _i，θ _v，φ，λ)

Wherein,

$k_1\left(\mathrm{\λ}\right)=\frac{1}{\mathrm{\π}}{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_2\left(\mathrm{\λ}\right)=-\frac{1}{\mathrm{\π}}{A}_w\left(\mathrm{\λ}\right)\·{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_3({\mathrm{\θ}}_i,{\mathrm{\θ}}_v,\mathrm{\φ},\mathrm{\λ})=\frac{2}{3\mathrm{\π}}{\mathrm{\ρ}}_c\left(\mathrm{\λ}\right)\·(k_{\mathrm{geo}}^c+1)+\frac{1}{\mathrm{\π}}{k}_{\mathrm{geo}}^g\·{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_4({\mathrm{\θ}}_i,{\mathrm{\θ}}_v,\mathrm{\φ},\mathrm{\λ})=-\frac{1}{\mathrm{\π}}{k}_{\mathrm{geo}}^g\·A_w\left(\mathrm{\λ}\right)\·{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_5({\mathrm{\θ}}_i,{\mathrm{\θ}}_v,\mathrm{\φ},\mathrm{\λ})=\frac{2}{3{\mathrm{\π}}^2}{\mathrm{\ρ}}_c\left(\mathrm{\λ}\right)\·k_{\mathrm{vol}}^{\mathrm{\ρ}}+\frac{2}{{3\mathrm{\π}}^2}{\mathrm{\τ}}_c\left(\mathrm{\λ}\right)\·k_{\mathrm{vol}}^{\mathrm{\τ}}+\frac{1}{3\mathrm{\π}}{\mathrm{\ρ}}_c\left(\mathrm{\λ}\right)+\frac{1-\sqrt{1-\mathrm{\ω}\left(\mathrm{\λ}\right)}}{\mathrm{\π}(1+2\cos \left({\mathrm{\θ}}_i\right)\sqrt{1-\mathrm{\ω}\left(\mathrm{\λ}\right)})}$

$k_{\mathrm{geo}}^g=\frac{1}{\mathrm{\π}}(\sec {\mathrm{\θ}}_i+\sec {\mathrm{\θ}}_v)\·(t-\cos t\·\sin t)-\sec {\mathrm{\θ}}_i-\sec {\mathrm{\θ}}_v+1$

$k_{\mathrm{geo}}^c=\sec {\mathrm{\θ}}_i\sec {\mathrm{\θ}}_v{\cos }^2(\mathrm{\ξ}/2)-1$

$k_{\mathrm{vol}}^{\mathrm{\ρ}}=\frac{(\mathrm{\π}-\mathrm{\ξ})\cos \mathrm{\ξ}+\sin \mathrm{\ξ}}{\cos {\mathrm{\θ}}_i+\cos {\mathrm{\θ}}_v}-\frac{\mathrm{\π}}{2}$

$k_{\mathrm{vol}}^{\mathrm{\τ}}=\frac{\mathrm{\ξ}\·\cos \mathrm{\ξ}+\sin \mathrm{\ξ}}{\cos {\mathrm{\θ}}_i+\cos {\mathrm{\θ}}_v}$

ξ＝arccos(cosθ _i?cosθ _v+sinθ _isinθ _vcosφ)

$t=\arccos \frac{\sqrt{D^2+{\left(\tan {\mathrm{\θ}}_i\tan {\mathrm{\θ}}_v\sin \mathrm{\φ}\right)}^2}}{\sec {\mathrm{\θ}}_i+\sec {\mathrm{\θ}}_v}$

D ²＝tan ²θ _i+tan ²θ _v-2tanθ _itanθ _vcosφ

ω(λ)＝ρ _c(λ)+τ _c(λ)

k ₁(λ), k ₂(λ), k ₃(θ _i, θ _v, φ, λ), k ₄(θ _i, θ _v, φ, λ), k ₅(θ _i, θ _v, φ λ) is kernel function, c ₁, c ₂, c ₃, c ₄, c ₅Be the karyonide number, ρ _c(λ) be the reflectivity of scatterer, τ _c(λ) be the transmitance of scatterer, λ is a wavelength, A _w(λ) be the moisture content absorption factor, ρ _g(m ₀, be λ) with reference to the soil reflectivity, θ _iBe said solar zenith angle, θ _vBe said observation zenith angle, φ is the relative bearing of the said sun and satellite.

Wherein, comprise among the step S6:

S61: adopt least square method through cost function to said system of linear equations, calculate said karyonide and count c ₁, c ₂, c ₃, c ₄, c ₅

S62: through the black hemisphere wave spectrum albedo of computes, white hemisphere wave spectrum albedo and real surface wave spectrum albedo,

${\mathrm{\α}}_{\mathrm{a\λ}}=\underset{x}{\mathrm{\Σ}}{c}_{x}H{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})$

${\mathrm{\α}}_{\mathrm{b\λ}}\left({\mathrm{\θ}}_i\right)=\underset{x}{\mathrm{\Σ}}{c}_x{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})$

${\mathrm{\α}}_{\mathrm{w\λ}}=\underset{x}{\mathrm{\Σ}}{c}_x{H}_x\left(\mathrm{\λ}\right)$

Wherein,

$h_x({\mathrm{\θ}}_i,\mathrm{\λ})=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}{\∫}_0^{\frac{\mathrm{\π}}{2}}{k}_x\sin {\mathrm{\θ}}_v\cos {\mathrm{\θ}}_{v}d{\mathrm{\θ}}_v\mathrm{d\φ}$

$H_x\left(\mathrm{\λ}\right)=2{\∫}_0^{\frac{\mathrm{\π}}{2}}{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})\sin {\mathrm{\θ}}_i\cos {\mathrm{\θ}}_{i}d{\mathrm{\θ}}_i$

Hh _x(θ _i，λ)＝d ₀h _x(θ _i，λ)+(1-d ₀)H _x(λ)

α _{B λ}(θ _i) for deceiving hemisphere wave spectrum albedo, α _{W λ}Be white hemisphere wave spectrum albedo, α _{A λ}Be real surface wave spectrum albedo, k _xBe x kernel function, c _xBe x karyonide number, d ₀For beam radia accounts for the ratio of built-up radiation, 1-d ₀Account for the ratio of built-up radiation for the skylight radiation.

Wherein, pass through the black hemisphere broadband albedo in the computes random wave segment limit, white hemisphere broadband albedo and real surface broadband albedo among the step S7,

${\mathrm{\α}}_{\mathrm{bb}}\left({\mathrm{\θ}}_i\right)=\underset{x}{\mathrm{\Σ}}{c}_x\·h_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)$

${\mathrm{\α}}_{\mathrm{bw}}=\underset{x}{\mathrm{\Σ}}{c}_x\·H_{x\_\mathrm{new}}$

${\mathrm{\α}}_{\mathrm{ba}}\left({\mathrm{\θ}}_i\right)=\underset{x}{\mathrm{\Σ}}{c}_x\·{\mathrm{Hh}}_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)$

Wherein,

$h_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})\mathrm{d\λ}$

$H_{x\_\mathrm{new}}={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{H}_x\left(\mathrm{\λ}\right)\mathrm{d\λ}$

${\mathrm{Hh}}_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}H{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})\mathrm{d\λ}$

α _Bb(θ _i) for deceiving hemisphere broadband albedo, α _BwBe white hemisphere broadband albedo, α _Ba(θ _i) be real surface broadband albedo, w _λBe the wave spectrum distribution of total descending radiant quantity, λ ₁Be the lower limit of wave band, λ ₂Higher limit for wave band.

Wherein, said observation satellite is A/B star of environment, and said observation satellite remotely-sensed data is the data that the charge coupled cell CCD camera of an A/B star of said Chinese environmental obtains.

Wherein, the wavelength band of said sensor is 0.3～3 μ m.

The invention also discloses a kind of surface albedo inverting system, comprising:

Acquisition module, the face of land two that is used for obtaining multi-angle through each pixel on the remotely-sensed data of observation satellite is to reflectivity data;

Select module; Be used for face of land type based on said each pixel of remotely-sensed data; Utilize the corresponding component spectral data of priori wave spectrum knowledge base selection, said priori wave spectrum knowledge base is the database of the storage face of land type and the component spectral data of correspondence, and said component spectral data is continuous;

Integration module is used for the band setting to different sensors, and said component spectral data is integrated to corresponding wave band;

Read module is used for reading the observation geometric data from said remotely-sensed data, and said observation geometric data comprises: the relative bearing of solar zenith angle, observation zenith angle and the sun and satellite;

Make up module, be used for component spectral data and the said observation geometric data substitution linear equation of the face of land two after reflectivity data and process integral processing, make up system of linear equations, comprise karyonide number and kernel function in the said linear equation according to said multi-angle;

Computing module is used for solving through least square method the karyonide number of said system of linear equations, according to the karyonide number and the kernel function of said linear equation, calculates black hemisphere wave spectrum albedo, white hemisphere wave spectrum albedo and real surface wave spectrum albedo;

The broadband computing module; Be used for according to the karyonide number of said linear equation and kernel function, and the wave spectrum of total descending radiant quantity distribute, calculate the black hemisphere broadband albedo in the random wave segment limit, white hemisphere broadband albedo and true broadband albedo.

(3) beneficial effect

The invention solves the measurement information quantity not sufficient, narrow wave band is very big to broadband conversion error in the traditional algorithm, and the low problem of remotely-sensed data utilization factor.

Description of drawings

Fig. 1 is the process flow diagram according to the surface albedo inversion method of one embodiment of the present invention;

Fig. 2 is the structured flowchart according to the surface albedo inverting system of one embodiment of the present invention.

Embodiment

Below in conjunction with accompanying drawing and embodiment, specific embodiments of the invention describes in further detail.Following examples are used to explain the present invention, but are not used for limiting scope of the present invention.

The present invention counts the wave spectrum independence face of land two to reflection model through setting up a kind of new face of land two to reflection model-karyonide; Its center not only is the function of the sun and observation angle; Be the function of wavelength simultaneously; The karyonide number only is expressed as the amount with canopy or surface infrastructure parameter correlation, realizes multi-angle and multiband joint inversion thus, solve the problem of observation information quantity not sufficient; In addition; Can kernel function carried out integration on the angle and on the wavelength simultaneously when calculating all band albedo utilizing the face of land two to ask to reflectivity; Integral result is counted weighted sum with karyonide just obtain all band albedo; Thereby the difficult problem that narrow wave band transforms to broadband in the solution traditional algorithm improves the precision and the algorithm versatility of inverting; Support to reflection model inverting surface albedo algorithm that RS data is collaborative based on newly-built two in addition and utilize, can solve the low problem of mass remote sensing data utilization factor.

Fig. 1 is according to the process flow diagram of the surface albedo inversion method of one embodiment of the present invention, may further comprise the steps:

S1: the face of land two that obtains multi-angle through each pixel on the remotely-sensed data of observation satellite is to reflectivity data.

S2: according to the face of land type of each pixel in the said remotely-sensed data; Select corresponding component spectral data from priori wave spectrum knowledge base; Said priori wave spectrum knowledge base is the database of the storage face of land type and the component spectral data of correspondence; Said component spectral data is continuous, and said face of land type comprises vegetation, soil, water body, snow etc., and said component spectral data comprises reflectivity, transmitance and soil, the water body of vegetation blade, reflectivity of snow etc.

S3: the wave band (in this embodiment, the wavelength band that sensor is set is 0.3～3 μ m) to different sensors is provided with, and said component spectral data is integrated to corresponding wave band, is specially, and carries out integration through following formula:

$R_i={\∫}_{{\mathrm{\λ}}_{\mathrm{si}}}^{{\mathrm{\λ}}_{\mathrm{ei}}}f\left(\mathrm{\λ}\right)\·R\left(\mathrm{\λ}\right)\mathrm{d\λ}$

Wherein, i representes i wave band, λ _SiAnd λ _EiBe respectively the initial wavelength and termination wavelength of wave band i, f (λ) is the wave spectrum response function of respective sensor, and R (λ) is continuous component spectral data, R _iBe component spectral data through integral processing.

S4: from the header file of said remotely-sensed data, read the observation geometric data, said observation geometric data comprises: the relative bearing of solar zenith angle, observation zenith angle and the sun and satellite;

S5: component spectral data and the said observation geometric data substitution linear equation of the face of land two after reflectivity data and process integral processing according to said multi-angle, make up system of linear equations, said linear equation is:

R(θ _i，θ _v，φ，λ)＝c ₁·k ₁(λ)+c ₂·k ₂(λ)+c ₃·k ₃(θ _i，θ _v，φ，λ)+c ₄·k ₄(θ _i，θ _v，φ，λ)+c ₅·k ₅(θ _i，θ _v，φ，λ)

Wherein, $k_1\left(\mathrm{\λ}\right)=\frac{1}{\mathrm{\π}}{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_2\left(\mathrm{\λ}\right)=-\frac{1}{\mathrm{\π}}{A}_w\left(\mathrm{\λ}\right)\·{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_3({\mathrm{\θ}}_i,{\mathrm{\θ}}_v,\mathrm{\φ},\mathrm{\λ})=\frac{2}{3\mathrm{\π}}{\mathrm{\ρ}}_c\left(\mathrm{\λ}\right)\·(k_{\mathrm{geo}}^c+1)+\frac{1}{\mathrm{\π}}{k}_{\mathrm{geo}}^g\·{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_4({\mathrm{\θ}}_i,{\mathrm{\θ}}_v,\mathrm{\φ},\mathrm{\λ})=-\frac{1}{\mathrm{\π}}{k}_{\mathrm{geo}}^g\·A_w\left(\mathrm{\λ}\right)\·{\mathrm{\ρ}}_g(m_0,\mathrm{\λ})$

$k_5({\mathrm{\θ}}_i,{\mathrm{\θ}}_v,\mathrm{\φ},\mathrm{\λ})=\frac{2}{3{\mathrm{\π}}^2}{\mathrm{\ρ}}_c\left(\mathrm{\λ}\right)\·k_{\mathrm{vol}}^{\mathrm{\ρ}}+\frac{2}{{3\mathrm{\π}}^2}{\mathrm{\τ}}_c\left(\mathrm{\λ}\right)\·k_{\mathrm{vol}}^{\mathrm{\τ}}+\frac{1}{3\mathrm{\π}}{\mathrm{\ρ}}_c\left(\mathrm{\λ}\right)+\frac{1-\sqrt{1-\mathrm{\ω}\left(\mathrm{\λ}\right)}}{\mathrm{\π}(1+2\cos \left({\mathrm{\θ}}_i\right)\sqrt{1-\mathrm{\ω}\left(\mathrm{\λ}\right)})}$

ω(λ)＝ρ _c(λ)+τ _c(λ)

R (θ _i, θ _v, φ, λ) be the face of land two to reflectivity, k ₁(λ), k ₂(λ), k ₃(θ _i, θ _v, φ, λ), k ₄(θ _i, θ _v, φ, λ), k ₅(θ _i, θ _v, φ λ) is kernel function, c ₁, c ₂, c ₃, c ₄, c ₅Be the karyonide number, ρ _c(λ) be the reflectivity of scatterer, τ _c(λ) be the transmitance of scatterer, λ is a wavelength, A _w(λ) be the moisture content absorption factor, ρ _g(m ₀, be λ) with reference to the soil reflectivity, θ _iBe said solar zenith angle, θ _vBe said observation zenith angle, φ is the relative bearing of the said sun and satellite.

S6:, specifically comprise according to said linear equation set of calculated real surface albedo:

S61: adopt least square method through cost function to said linear equation, calculate said karyonide and count c ₁, c ₂, c ₃, c ₄, c ₅, particularly, refutation process employing method is a least square method, obtains optimum solution through the calculation cost function, through M observation geometric data and N wave band of participating in inverting, said system of linear equations is here:

R(θ _im，θ _vm，φ _m，λ _n)＝c ₁k ₁(λ _n)+c ₂k ₂(λ _n)+c ₃k ₃(θ _im，θ _vm，φ _m，λ _n)+c ₄k ₄(θ _im，θ _vm，φ _m，λ _n)+c ₅k ₅(θ _im，θ _vm，φ _m，λ _n)

Wherein, m=1～M, m are integer, θ _ImBe m the solar zenith angle in the observation geometric data, θ _VmBe m the observation zenith angle in the observation geometric data, φ _mBe the sun and the relative bearing of satellite of m observation in the geometric data, n=1～N, n are integer, λ _nBe n wave band.

Said cost equation is:

$\mathrm{COST}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\frac{{({R^{\mathrm{obs}}}_{m,n}-{R^{\mathrm{simu}}}_{m,n})}^2}{W_m\·W_n}$

In the formula, R ^Obs _{M, s}The face of land two that is m observation geometric data n wave band is to the reflectivity measured value, R ^Simu _{M, n}The face of land two of n wave band of m observation geometric data is to the reflectivity analogue value, and the said face of land two is obtained by the sensor on the observation satellite to the reflectivity measured value, and the said face of land two is to the value of the reflectivity analogue value for obtaining through said linear equation calculating, W _mBe m the weight coefficient on the observation geometric data, W _nBe the weight coefficient on n the wave band, W _mAnd W _nBe provided with and can confirm that for example, the refutation strategy final purpose that this paper proposes is the extraction of accomplishing surface albedo to the quality condition and the inverting purpose of reflectivity data according to the face of land two of multi-angle; Its final calculating effect quality depend on the spectral energy of skip band distribute relevant, under the atmospheric conditions of drying, visible light, near infrared and in infrared in total sun power shared proportion basic fixed; Wherein visible light partly accounts for 52.6, and near infrared partly accounts for 36.2, and middle infrared part accounts for 11.2; It is bigger that the face of land two of the wave band that energy distribution is concentrated is extracted influence to reflectivity inverting order of accuarcy to albedo, for this reason, and at our given bigger weight coefficient when the inverting karyonide is counted of such wave band; Make its match residual error COST that tries one's best keep little value, through utilizing the big institute of MODTRAN radiation delivery simulation softward, the descending built-up radiation that the ground surface that we to clear sky under the typical atmospheric state are receives is simulated; Account for the coefficient [21.8,24.7,35.3 of the given corresponding weights of ratio value of whole skip band according to the radiant quantity of each wave band of MODIS; 29.0,5.4,7.3; 1]; Obtain the weight coefficient matrix of the first seven wave band of MODIS, to the observed reading on the different observation angles, we think that its quality of data is suitable in observation process; Therefore given identical weighted value is not done special processing. and can be to other satellite sensor band setting according to the weight coefficient of the given different-waveband data of similar approach.

S62: through the black hemisphere wave spectrum albedo of computes, white hemisphere wave spectrum albedo and real surface wave spectrum albedo,

${\mathrm{\α}}_{\mathrm{a\λ}}=\underset{x}{\mathrm{\Σ}}{c}_{x}H{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})$

${\mathrm{\α}}_{\mathrm{b\λ}}\left({\mathrm{\θ}}_i\right)=\underset{x}{\mathrm{\Σ}}{c}_x{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})$

${\mathrm{\α}}_{\mathrm{w\λ}}=\underset{x}{\mathrm{\Σ}}{c}_x{H}_x\left(\mathrm{\λ}\right)$

Wherein,

$h_x({\mathrm{\θ}}_i,\mathrm{\λ})=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}{\∫}_0^{\frac{\mathrm{\π}}{2}}{k}_x\sin {\mathrm{\θ}}_v\cos {\mathrm{\θ}}_{v}d{\mathrm{\θ}}_v\mathrm{d\φ}$

$H_x\left(\mathrm{\λ}\right)=2{\∫}_0^{\frac{\mathrm{\π}}{2}}{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})\sin {\mathrm{\θ}}_i\cos {\mathrm{\θ}}_{i}d{\mathrm{\θ}}_i$

Hh _x(θ _i，λ)＝d ₀h _x(θ _i，λ)+(1-d ₀)H _x(λ)

α _{B λ}(θ _i) for deceiving hemisphere wave spectrum albedo, α _{W λ}Be white hemisphere wave spectrum albedo, α _{A λ}Be real surface wave spectrum albedo, k _xBe x kernel function, c _xBe x karyonide number, d ₀For beam radia accounts for the ratio of built-up radiation, 1-d ₀Account for the ratio of built-up radiation for the skylight radiation.

S7: calculate the black hemisphere broadband albedo in the random wave segment limit according to said real surface albedo; Bai Banqiu broadband albedo; And real surface broadband albedo; Be specially; The energy distribution (by different atmospheric condition decisions) of total descending wave spectrum radiant quantity on different wave spectrums that broadband albedo (broadband albedo) can be received by the face of land of wave spectrum albedo (spectral albedo) and atmosphere bottom according to its definition; Through computes broadband albedo

${\mathrm{\α}}_{\mathrm{broad}}=\frac{{\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}G\left(\mathrm{\λ}\right){\mathrm{\α}}_{\mathrm{\λ}}\mathrm{d\λ}}{{\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}G\left(\mathrm{\λ}\right)\mathrm{d\λ}}={\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}{w}_{\mathrm{\λ}}{\mathrm{\α}}_{\mathrm{\λ}}\mathrm{d\λ}$

Wherein,

$w_{\mathrm{\λ}}=\frac{G\left(\mathrm{\λ}\right)}{{\∫}_{{\mathrm{\λ}}_1}^{{\mathrm{\λ}}_2}G\left(\mathrm{\λ}\right)\mathrm{d\λ}}$

w _λBe the wave spectrum distribution of total descending radiant quantity, G (λ) is the total descending wave spectrum radiant quantity that the face of land of atmosphere bottom receives, λ ₁Be the lower limit of wave band, λ ₂Be the higher limit of wave band, α _BroadBe the broadband albedo.

In this embodiment, through the black hemisphere broadband albedo in the computes random wave segment limit, white hemisphere broadband albedo and real surface broadband albedo,

${\mathrm{\α}}_{\mathrm{bb}}\left({\mathrm{\θ}}_i\right)={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{\mathrm{\α}}_{\mathrm{b\λ}}\left({\mathrm{\θ}}_i\right)\mathrm{d\λ}=\underset{x}{\mathrm{\Σ}}{c}_x\·{\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})\mathrm{d\λ}=\underset{x}{\mathrm{\Σ}}{c}_x\·h_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)$

${\mathrm{\α}}_{\mathrm{bw}}={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{\mathrm{\α}}_{\mathrm{w\λ}}\mathrm{d\λ}=\underset{x}{\mathrm{\Σ}}{c}_x\·{\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{H}_x\left(\mathrm{\λ}\right)\mathrm{d\λ}=\underset{x}{\mathrm{\Σ}}{c}_x\·H_{x\_\mathrm{new}}$

${\mathrm{\α}}_{\mathrm{ba}}\left({\mathrm{\θ}}_i\right)={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{\mathrm{\α}}_{\mathrm{a\λ}}\mathrm{d\λ}=\underset{x}{\mathrm{\Σ}}{c}_x\·{\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}H{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})\mathrm{d\λ}=\underset{x}{\mathrm{\Σ}}{c}_x\·{\mathrm{Hh}}_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)$

Wherein,

$h_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{h}_x({\mathrm{\θ}}_i,\mathrm{\λ})\mathrm{d\λ}$

$H_{x\_\mathrm{new}}={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{H}_x\left(\mathrm{\λ}\right)\mathrm{d\λ}$

${\mathrm{Hh}}_{x\_\mathrm{new}}\left({\mathrm{\θ}}_i\right)={\∫}_{{\mathrm{\λ}}_2}^{{\mathrm{\λ}}_1}{w}_{\mathrm{\λ}}{\mathrm{Hh}}_x({\mathrm{\θ}}_i,\mathrm{\λ})\mathrm{d\λ}$

α _Bb(θ _i) for deceiving hemisphere broadband albedo, α _BwBe white hemisphere broadband albedo, α _Ba(θ _i) be real surface broadband albedo, w _λBe the wave spectrum distribution of total descending radiant quantity, λ ₁Be the lower limit of wave band, λ ₂Higher limit for wave band.

Because Hh _{X_new}(θ _i), h _{X_new}(θ _i) and H _{X_new}All be and the irrelevant amount of the angle distribution of the observation of wavelength and multi-angle albedo data to be stored as look-up table shown in therefore can calculating in advance and according to face of land type, atmospheric condition and observation zenith angle, like this α _Ba(θ _i) the karyonide number that can directly utilize inverting to obtain obtains the integral kernel numerical evaluation by look-up table and obtain so can significantly improving counting yield; And when nuclear is carried out the wavelength integration; Can realize obtaining the broadband albedo by continuous wave spectrum albedo integration, the wave spectrum integration interval during calculating has identical spectral resolution with the component wave spectrum, and whole process has tight Fundamentals of Mathematics; Go for the remotely-sensed data in various sources, broken away from the traditional algorithm restriction observation platform and sensing.

The invention also discloses a kind of surface albedo inverting system, as shown in Figure 2, comprising:

Acquisition module, the face of land two that is used for obtaining multi-angle through each pixel on the remotely-sensed data of observation satellite is to reflectivity data;

Select module; Be used for face of land type based on said each pixel of remotely-sensed data; Utilize the corresponding component spectral data of priori wave spectrum knowledge base selection, said priori wave spectrum knowledge base is the database of the storage face of land type and the component spectral data of correspondence, and said component spectral data is continuous;

Integration module is used for the band setting to different sensors, and said component spectral data is integrated to corresponding wave band;

Read module is used for reading the observation geometric data from said remotely-sensed data, and said observation geometric data comprises: the relative bearing of solar zenith angle, observation zenith angle and the sun and satellite;

Make up module, be used for component spectral data and the said observation geometric data substitution linear equation of the face of land two after reflectivity data and process integral processing, make up system of linear equations, comprise karyonide number and kernel function in the said linear equation according to said multi-angle;

Computing module is used for solving through least square method the karyonide number of said system of linear equations, according to the karyonide number and the kernel function of said linear equation, calculates black hemisphere wave spectrum albedo, white hemisphere wave spectrum albedo and real surface wave spectrum albedo;

The broadband computing module; Be used for according to the karyonide number of said linear equation and kernel function, and the wave spectrum of total descending radiant quantity distribute, calculate the black hemisphere broadband albedo in the random wave segment limit, white hemisphere broadband albedo and true broadband albedo.

Above embodiment only is used to explain the present invention; And be not limitation of the present invention; The those of ordinary skill in relevant technologies field under the situation that does not break away from the spirit and scope of the present invention, can also be made various variations and modification; Therefore all technical schemes that are equal to also belong to category of the present invention, and scope of patent protection of the present invention should be defined by the claims.

Claims (5)

1. a surface albedo inversion method is characterized in that, may further comprise the steps:

S1: the face of land two that obtains multi-angle through each pixel on the remotely-sensed data of observation satellite is to reflectivity data;

S2: based on the face of land type of each pixel in the said remotely-sensed data; Select corresponding component spectral data from priori wave spectrum knowledge base; Said priori wave spectrum knowledge base is the database of the storage face of land type and the component spectral data of correspondence, and said component spectral data is continuous;

S3: the band setting to different sensors is integrated to corresponding wave band with said component spectral data;

S4: from said remotely-sensed data, read the observation geometric data, said observation geometric data comprises: the relative bearing of solar zenith angle, observation zenith angle and the sun and satellite;

S5: component spectral data and the said observation geometric data substitution linear equation of the face of land two after reflectivity data and process integral processing according to said multi-angle, make up system of linear equations, comprise karyonide number and kernel function in the said linear equation;

S6: solve the karyonide number of said system of linear equations through least square method,, calculate black hemisphere wave spectrum albedo, white hemisphere wave spectrum albedo and real surface wave spectrum albedo according to the karyonide number and the kernel function of said linear equation;

S7: according to the karyonide number of said linear equation and kernel function, and the wave spectrum of total descending radiant quantity distribute, calculate the black hemisphere broadband albedo in the random wave segment limit, white hemisphere broadband albedo and true broadband albedo;

Wherein, in step S3, carry out integration through following formula:

Wherein, i representes i wave band, λ _SiAnd λ _EiBe respectively the initial wavelength and termination wavelength of i wave band, f (λ) is the wave spectrum response function of respective sensor, and R (λ) is continuous component spectral data, R _iBe component spectral data through integral processing;

Linear equation described in the step S5 is:

R(θ _i，θ _v，φ，λ)＝c ₁·k ₁(λ)+c ₂·k ₂(λ)+c ₃·k ₃(θ _i，θ _v，φ，λ)+c ₄·k ₄(θ _i，θ _v，φ，λ)+c ₅·k ₅(θ _i，θ _v，φ，λ)

Wherein,

ξ＝arccos(cosθ _icosθ _v+sinθ _isinθ _vcosφ)

D ²＝tan ²θ _i+tan ²θ _v-2tanθ _itanθ _vcosφ

ω(λ)＝ρ _c(λ)+τ _c(λ)

R (θ _i, θ _v, φ, λ) be the face of land two to reflectivity, k ₁(λ), k ₂(λ), k ₃(θ _i, θ _v, φ, λ), k ₄(θ _i, θ _v, φ, λ), k ₅(θ _i, θ _v, φ λ) is kernel function, c ₁, c ₂, c ₃, c ₄, c ₅Be the karyonide number, ρ _c(λ) be the reflectivity of scatterer, τ _c(λ) be the transmitance of scatterer, λ is a wavelength, A _w(λ) be the moisture content absorption factor, ρ _g(m ₀, be λ) with reference to the soil reflectivity, θ _iBe said solar zenith angle, θ _vBe said observation zenith angle, φ is the relative bearing of the said sun and satellite;

Comprise among the step S6:

S61: adopt least square method through cost function to said system of linear equations, calculate said karyonide and count c ₁, c ₂, c ₃, c ₄, c ₅

S62: through the black hemisphere wave spectrum albedo of computes, white hemisphere wave spectrum albedo and real surface wave spectrum albedo,

Wherein,

Hh _x(θ _i，λ)＝d ₀h _x(θ _i，λ)+(1-d ₀)H _x(λ)

α _{B λ}(θ _i) for deceiving hemisphere wave spectrum albedo, α _{W λ}Be white hemisphere wave spectrum albedo, α _{A λ}Be real surface wave spectrum albedo, k _xBe x kernel function, c _xBe x karyonide number, d ₀For beam radia accounts for the ratio of built-up radiation, 1-d ₀Account for the ratio of built-up radiation for the skylight radiation;

Pass through the black hemisphere broadband albedo in the computes random wave segment limit, white hemisphere broadband albedo and real surface broadband albedo among the step S7,

Wherein,

α _Bb(θ _i) for deceiving hemisphere broadband albedo, α _BwBe white hemisphere broadband albedo, α _Ba(θ _i) be real surface broadband albedo, w _λBe the wave spectrum distribution of total descending radiant quantity, λ ₁Be the lower limit of wave band, λ ₂Higher limit for wave band.

2. surface albedo inversion method as claimed in claim 1 is characterized in that, the type of the face of land described in the step S2 comprises: vegetation, soil, water body and snow.

3. like each described surface albedo inversion method of claim 1-2, it is characterized in that said observation satellite is A/B star of environment, said observation satellite remotely-sensed data is the data that the charge coupled cell CCD camera of an A/B star of said environment obtains.

4. like each described surface albedo inversion method of claim 1-2, it is characterized in that the wavelength band of said sensor is 0.3～3 μ m.

5. a surface albedo inverting system is characterized in that, comprising:

Acquisition module, the face of land two that is used for obtaining multi-angle through each pixel on the remotely-sensed data of observation satellite is to reflectivity data;

Select module; Be used for face of land type based on said each pixel of remotely-sensed data; Utilize the corresponding component spectral data of priori wave spectrum knowledge base selection, said priori wave spectrum knowledge base is the database of the storage face of land type and the component spectral data of correspondence, and said component spectral data is continuous;

Integration module is used for the band setting to different sensors, and said component spectral data is integrated to corresponding wave band;

Read module is used for reading the observation geometric data from said remotely-sensed data, and said observation geometric data comprises: the relative bearing of solar zenith angle, observation zenith angle and the sun and satellite;

Make up module, be used for component spectral data and the said observation geometric data substitution linear equation of the face of land two after reflectivity data and process integral processing, make up system of linear equations, comprise karyonide number and kernel function in the said linear equation according to said multi-angle;

Computing module is used for solving through least square method the karyonide number of said system of linear equations, according to the karyonide number and the kernel function of said linear equation, calculates black hemisphere wave spectrum albedo, white hemisphere wave spectrum albedo and real surface wave spectrum albedo;

The broadband computing module; Be used for based on the karyonide number of said linear equation and kernel function, and the wave spectrum of total descending amount of radiation distribute, calculate the black hemisphere broadband albedo in the random wave segment limit, white hemisphere broadband albedo and true broadband albedo;

Wherein, in integration module, carry out integration through following formula:

Wherein, i representes i wave band, λ _SiAnd λ _EiBe respectively the initial wavelength and termination wavelength of i wave band, f (λ) is the wave spectrum response function of respective sensor, and R (λ) is continuous component spectral data, R _iBe component spectral data through integral processing;

Making up the linear equation described in the module is:

R(θ _i，θ _v，φ，λ)＝c ₁·k ₁(λ)+c ₂·k ₂(λ)+c ₃·k ₃(θ _i，θ _v，φ，λ)+c ₄·k ₄(θ _i，θ _v，φ，λ)+c ₅·k ₅(θ _i，θ _v，φ，λ)

Wherein,

ξ＝arccos(cosθ _icosθ _v+sinθ _isinθ _vcosφ)

D ²＝tan ²θ _i+tan ²θ _v-2tanθ _itanθ _vcosφ

ω(λ)＝ρ _c(λ)+τ _c(λ)

R (θ _i, θ _v, φ, λ) be the face of land two to reflectivity, k ₁(λ), k ₂(λ), k ₃(θ _i, θ _v, φ, λ), k ₄(θ _i, θ _v, φ, λ), k ₅(θ _i, θ _v, φ λ) is kernel function, c ₁, c ₂, c ₃, c ₄, c ₅Be the karyonide number, ρ _c(λ) be the reflectivity of scatterer, τ _c(λ) be the transmitance of scatterer, λ is a wavelength, A _w(λ) be the moisture content absorption factor, ρ _g(m ₀, be λ) with reference to the soil reflectivity, θ _iBe said solar zenith angle, θ _vBe said observation zenith angle, φ is the relative bearing of the said sun and satellite;

Comprise in the computing module:

Karyonide is counted calculating sub module, is used for adopting least square method through cost function to said system of linear equations, calculates said karyonide and counts c ₁, c ₂, c ₃, c ₄, c ₅

The albedometer operator module is used for through the black hemisphere wave spectrum albedo of computes, white hemisphere wave spectrum albedo and real surface wave spectrum albedo,

Wherein,

Hh _x(θ _i，λ)＝d ₀h _x(θ _i，λ)+(1-d ₀)H _x(λ)

α _{B λ}(θ _i) for deceiving hemisphere wave spectrum albedo, α _{W λ}Be white hemisphere wave spectrum albedo, α _{A λ}Be real surface wave spectrum albedo, k _xBe x kernel function, c _xBe x karyonide number, d ₀For beam radia accounts for the ratio of built-up radiation, 1-d ₀Account for the ratio of built-up radiation for the skylight radiation;

Pass through the black hemisphere broadband albedo in the computes random wave segment limit, white hemisphere broadband albedo and real surface broadband albedo in the broadband computing module,

Wherein,

α _Bb(θ _i) for deceiving hemisphere broadband albedo, α _BwBe white hemisphere broadband albedo, α _Ba(θ _i) be real surface broadband albedo, w _λBe the wave spectrum distribution of total descending radiant quantity, λ ₁Be the lower limit of wave band, λ ₂Higher limit for wave band.

Images