CN104049256A - Physical method for computing atmospheric path radiance of satellite remote sensing images through picture elements one by one - Google Patents

Abstract

A dark object method is mainly utilized when path radiance of satellite remote sensing images is computed traditionally. However, the dark object method is an empirical method; by means of the dark object method, the whole scene of images has one path radiance value only, and the error is big. Therefore, the invention provides a physical method for computing the atmospheric path radiance of the satellite remote sensing images through picture elements one by one. The physical method is based on the remote sensing images completely, and no auxiliary data is needed. Accordingly, the physical method is strong in universality, and a new method is provided for computation of a large area of atmospheric path radiance.

Description

A kind of physical method that calculates satellite remote-sensing image atmospheric path radiation by pixel

Technical field

The present invention relates to a kind of method that satellite remote-sensing image atmospheric path radiation obtains, can be applied in the industry departments such as agricultural, forestry, meteorology, ecologic environment.

Background technology

Solar radiation in atmospheric propagation process, run into atmosphere in the particle such as gas molecule, gasoloid, ice crystal, make part light change direction, and scatter to all directions, be called scattering.Atmospheric path radiation is that solar radiation directly arrives the radiation of sensor in propagation in atmosphere process after each component scattering in atmosphere.In remote sensing image, atmospheric path radiation does not comprise spectrum and the spatial information of any atural object, and atmospheric path radiation value has reduced the contrast of atural object in remote sensing image, the fuzzy detailed information of atural object, affects the quality of remote sensing image, is one of important content of remote sensing images atmospheric correction.

Traditional remote sensing images atmospheric path radiation calculates and mainly utilizes dark goal method, on the method supposition remote sensing images, there is dark target (water body, shade etc.), the actual reflectance of supposing dark target is 0, impact due to atmospheric path radiation, causing the reflectivity of dark target on image is not 0, and this part of increase is atmospheric path radiation, and this method is experimental, and whole scape image only has a journey radiation value, there is larger error.The theoretical foundation that the present invention calculates from journey radiation, a kind of physical method that calculates satellite remote-sensing image atmospheric path radiation by pixel has been proposed, the method is completely based on remote sensing image itself, do not need auxiliary data, there is stronger universality, for a kind of new method that provides is provided in large regional atmospheric journey radiation.

Summary of the invention

The object of the present invention is to provide a kind of physical method that calculates satellite remote-sensing image atmospheric path radiation by pixel, to overcome now methodical deficiency, thus the development of promotion association area remote sensing application.

For achieving the above object, the method that the present invention proposes comprises the following steps:

The first step, satellite image is carried out to geometry correction;

Second step, according to the longitude and latitude of Landsat remote sensing image and video imaging time, the MODIS aerosol optical depth that obtain same geographic range, synchronously passes by (0.47 μ m and 0.66 μ m), total atmospheric water steam content, earth's surface air pressure and ozone concentration image, and each image is done to geometry correction, be converted to the projection pattern consistent with Landsat remote sensing image and resolution, and image ranks number is consistent;

The 3rd step, by total atmospheric water steam content W, calculate atmospheric water steam optical thickness;

${\mathrm{\τ}}_w=\frac{0.2385a_{\mathrm{w\λ}}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\λ}}\mathrm{wM})}^{0.45}}$

Wherein, τ _watmospheric water steam optical thickness, a _{w λ}be water vapor absorption coefficient, w is precipitable water vapor, precipitable water vapor w and total atmospheric water steam content W numerically equal, and M is relative atm number;

The 4th step, calculating Ozone Absorption optical thickness

τ _o＝C _ozone*A _ozone(λ)

τ _oozone Absorption optical thickness, C _ozonefor MODIS ozone concentration (unit is Dobson), A _ozone(λ) be the Ozone Absorption coefficient of wave band λ;

The 5th step, by the aerosol optical depth of MODIS0.47 μ m and 0.66 μ m, utilize following formula to calculate Wavelength Indices (α) and atmospheric turbidity coefficient (β)

$\mathrm{\α}=-\frac{\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_1\right)-\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_2\right)}{\ln {\mathrm{\λ}}_1-\ln {\mathrm{\λ}}_2}$

β＝τ _a(λ ₁)λ ₁ ^α＝τ _a(λ ₂)λ ₂ ^α

λ ₁and λ ₂be respectively 0.47 μ m and 0.66 μ m, τ _a(λ ₁) and τ _a(λ ₂) be respectively the aerosol optical depth of 0.47 μ m and 0.66 μ m;

After obtaining α and β, can utilize following formula to obtain the aerosol optical depth τ of any wavelength X _a(λ):

τ _a(λ)＝βλ ^-α

Wherein, λ is wavelength (μ m), τ _a(λ) be the aerosol optical depth of this wavelength;

The 6th step, calculating Rayleigh scattering optics thickness;

${\mathrm{\τ}}_r=\frac{p}{1013.25}0.008569{\mathrm{\λ}}^{-4}(1+0.0133{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})$

Wherein, τ _rfor Rayleigh scattering optics thickness, λ is the centre wavelength (μ m) of each wave band of image; P is the earth's surface air pressure product of MODIS;

The 7th step, calculating shine upon the atmospheric transmittance of direction and sensor observed ray;

T _z＝exp(-τ/cosθ _z)＝exp{(-τ _r-τ _a-τ _o-τ _w)/cosθ _z}

T _v＝exp(-τ/cosθ _v)＝exp{(-τ _r-τ _a-τ _o-τ _w)/cosθ _v}

Wherein, T _zand T _vbe respectively the atmospheric transmittance that shines upon direction and sensor observed ray, τ is atmosphere optical thickness, θ _zsolar zenith angle, θ _vsensor observation zenith angle, τ _r, τ _a, τ _oand τ _wrespectively Rayleigh scattering optics thickness, aerosol optical depth, Ozone Absorption optical thickness and atmospheric water steam optical thickness;

The 8th step, calculating solar distance

Solar distance d ²(astronomical unit) is calculated as follows:

$\mathrm{\Γ}=\frac{360}{365}(\mathrm{dn}-1)$

$d^2=\frac{1}{1.00011+0.034221\cos \mathrm{\Γ}+0.00128\sin \mathrm{\Γ}+0.000719\cos \mathrm{\Γ}+0.00077\sin 2\mathrm{\Γ}}$

Wherein, dn is Julian date, and the image capturing date is apart from the number of days on January 1;

The 9th step, calculating scattering angle

θ _p＝180°-θ _z-θ _v

Wherein, θ _pfor scattering angle;

The tenth step, suppose that atmosphere is that homogeneous atmosphere, scattering angle are greater than 90 while spending, the computing formula of atmospheric path radiation is:

$L_p=\frac{E_0}{d^2}{\mathrm{\τ}}_r\frac{3}{16\mathrm{\π}}(1+{\cos }^2{\mathrm{\θ}}_p)(T_z{T}_v-1)/\ln \left(T_z{T}_v\right)$

Wherein, L _pfor journey radiation, E ₀it is the solar spectrum irradiancy of exoatmosphere respective wavelength.

Embodiment

Atmospheric path radiation is the result of Rayleigh scattering and gasoloid Mie scattering combined action, and the total atmospheric path radiation that arrives sensor can be regarded the integration (Fig. 1) of the up radiation of atmosphere in earth's surface-sensor path as, and the theoretical calculation formula of journey radiation is as follows:

$L_p=\underset{0}{\overset{Z_p}{\∫}}{E}_0/d^2{T}_{z1}{\mathrm{\β}}_{\mathrm{sca}}({\mathrm{\θ}}_p,\mathrm{\λ})T_{v1}\mathrm{dr}---\left(1\right)$

Wherein, L _pfor journey radiation, Z _pfor total air path of earth's surface-sensor orientation, E ₀be the solar spectrum irradiancy of exoatmosphere respective wavelength, can be calculated by detector response function d ²for solar distance, T _z1representative shines upon the path transmitance in direction, T _v1represent the path transmitance in earth's surface-sensor orientation, the air path that r is earth's surface-sensor orientation, θ _pfor scattering angle, shine upon the angle of direction and earth's surface-sensor orientation, λ is wavelength, β _sca(θ _p, λ) being scattered through angles coefficient, i.e. the solar radiation of wavelength X is dispersed into the ratio of earth's surface-sensor orientation.

Scattered through angles factor beta _sca(θ _p, λ) be Rayleigh scattering scattered through angles factor beta _r(θ _p, λ) with Mie scattering scattered through angles factor beta _m(θ _p, λ) sum.Scatteringangleθ _p=180 °-θ _z-θ _v, θ _zsolar zenith angle, θ _vit is sensor observation zenith angle.Rayleigh scattering is relevant with scattering angle, and the relation between them can be used phase function (p (θ _p)) expression (Schott J R, 1997, Remote Sensing-The Image Chain Approach (NewYork:Oxford University Press)):

$p\left({\mathrm{\θ}}_p\right)=\frac{{\mathrm{\β}}_r({\mathrm{\θ}}_p,\mathrm{\λ})}{({\mathrm{\β}}_r\left(\mathrm{\λ}\right)/4\mathrm{\π})}=\frac{3}{4}(1+{\cos }^2{\mathrm{\θ}}_p)$

Wherein, β _r(λ) be β _r(θ _pλ) the integration on all angle directions, Rayleigh scattering can be to all directions scattering, and the actual phase function of the Mie scattering that gasoloid etc. cause mainly concentrates on forward scattering, and back scattering is considerably less, therefore, when scattering angle is greater than 90 while spending, Mie scattering can be ignored, and journey radiation is only caused by Rayleigh scattering, like this, scattered through angles coefficient can be expressed as:

${\mathrm{\β}}_{\mathrm{sca}}({\mathrm{\θ}}_p,\mathrm{\λ})={\mathrm{\β}}_r({\mathrm{\θ}}_p,\mathrm{\λ})=\frac{3}{4}(1+{\cos }^2{\mathrm{\θ}}_p)\frac{{\mathrm{\β}}_r\left(\mathrm{\λ}\right)}{4\mathrm{\π}}---\left(2\right)$

Fig. 1 atmospheric path radiation diagram

Suppose that atmosphere is homogeneous atmosphere, attenuation coefficient is k _λ, air path is S, atmospheric transmittance T (λ) can be expressed as:

T(λ)＝exp[-k _λS]

In Fig. 1, S _z1for shining upon the air path in direction, T _z1for air path is S _z1time atmospheric transmittance; S ₁air path while arriving earth's surface for sunray, T ₂for air path is S ₁time atmospheric transmittance,

T _z1＝exp[-k _λS _z1]

T _z＝exp[-k _λS ₁]

: $T_{z1}={T_z}^{\frac{S_{z1}}{S_1}}$

In Fig. 1, r is the air path in earth's surface-sensor orientation, T _v1atmospheric transmittance while being r for air path; Z _pfor the total air path in earth's surface-sensor orientation, T _vfor air path is Z _ptime atmospheric transmittance,

T _v1＝exp[-k _λr]

T _v＝exp[-k _λZ _p]

: $T_{v1}={T_v}^{\frac{r}{Z_p}}---\left(3\right)$

Because earth radius is much larger than atmospheric envelope thickness, therefore AB section camber line can be similar to and sees and be in line in Fig. 1, by triangle similarity relation, can be obtained:

$\frac{S_{z1}}{S_1}=\frac{r}{Z_p}$

Therefore: $T_{z1}={T_z}^{\frac{r}{Z_p}}---\left(4\right)$

(3) (4) are brought into formula (1), obtain:

$L_p=\underset{0}{\overset{Z_p}{\∫}}{E}_0/d^2{T_z}^{\frac{r}{Z_p}}{\mathrm{\β}}_{\mathrm{sca}}({\mathrm{\θ}}_p,\mathrm{\λ}){T_v}^{\frac{r}{Z_p}}\mathrm{dr}$

Under homogeneous atmosphere condition, β _sca(θ _p, λ) can regard constant as:

$L_p=E_0/d^2{\mathrm{\β}}_{\mathrm{sca}}({\mathrm{\θ}}_p,\mathrm{\λ})\underset{0}{\overset{Z_p}{\∫}}{\left(T_z{T}_v\right)}^{\frac{r}{Z_p}}\mathrm{dr}$

L _p＝E ₀/d ²Iβ _sca(θ _p，λ)Z _p(T _zT _v-1)/ln(T _zT _v)

By formula (2) substitution above formula, obtain

$L_p=E_0/d^2{\mathrm{\β}}_r\left(\mathrm{\λ}\right)Z_p\frac{3}{16\mathrm{\π}}(1+{\cos }^2{\mathrm{\θ}}_p)(T_z{T}_v-1)/\ln \left(T_z{T}_v\right)$

β _r(λ) Z _pthe product that represents rayleigh scattering coefficient and air path, is Rayleigh scattering optics thickness:

$L_p=E_0/d^2{\mathrm{\τ}}_r\left(\mathrm{\λ}\right)\frac{3}{16\mathrm{\π}}(1+{\cos }^2{\mathrm{\θ}}_p)(T_z{T}_v-1)/\ln \left(T_z{T}_v\right)$

T _z＝exp(-τ/cosθ _z)＝exp{(-τ _r-τ _a-τ _o-τ _w)/cosθ _z}

T _v＝exp(-τ/cosθ _v)＝exp{(-τ _r-τ _a-τ _o-τ _w)/cosθ _v}

Wherein, τ is atmosphere optical thickness, θ _zsolar zenith angle, θ _vsensor observation zenith angle, τ _r, τ _a, τ _oand τ _wrespectively Rayleigh scattering optics thickness, aerosol optical depth, Ozone Absorption optical thickness and atmospheric water steam optical thickness;

${\mathrm{\τ}}_r=\frac{p}{1013.25}0.008569{\mathrm{\λ}}^{-4}(1+0.0133{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})$

Wherein, τ _rfor Rayleigh scattering optics thickness, λ is the centre wavelength of each wave band of image; P is the earth's surface air pressure product of MODIS;

By the aerosol optical depth of MODIS0.47 μ m and 0.66 μ m, utilize following formula to calculate Wavelength Indices (α) and atmospheric turbidity coefficient (β)

$\mathrm{\α}=-\frac{\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_1\right)-\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_2\right)}{\ln {\mathrm{\λ}}_1-\ln {\mathrm{\λ}}_2}$

β＝τ _a(λ ₁)λ ₁ ^α＝τ _a(λ ₂)λ ₂ ^α

λ ₁and λ ₂be respectively 0.47 μ m and 0.66 μ m, τ _a(λ ₁) and τ _a(λ ₂) be respectively the aerosol optical depth of 0.47 μ m and 0.66 μ m;

After obtaining α and β, can utilize following formula to obtain the aerosol optical depth τ of any wavelength X _a(λ):

τ _a(λ)＝βλ ^-α

Wherein, λ is wavelength (μ m), τ _a(λ) be the aerosol optical depth of this wavelength;

Ozone Absorption opticalthicknessτ _opressing following formula obtains:

τ _o＝C _ozone*A _ozone(λ)

τ _oozone Absorption optical thickness, C _ozonefor MODIS ozone concentration (unit is Dobson), Aozone (λ) is the Ozone Absorption coefficient of wave band λ, can obtain according to document Aerosol Optical Depth Value-Added Product (Koontz A et al.2013).

https：//www.arm.gov/publications/tech_reports/doe-sc-arm-tr-129.pdf

Atmospheric water steam opticalthicknessτ _wby total atmospheric water steam content W, calculated:

${\mathrm{\τ}}_w=\frac{0.2385a_{\mathrm{w\λ}}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\λ}}\mathrm{wM})}^{0.45}}$

Wherein, τ _watmospheric water steam optical thickness, a _{w λ}water vapor absorption coefficient, a of each wavelength _{w λ}value can read up the literature and obtain [Bird R E, Riordan C.Simple solar spectral model for direct and diffuse irradiance on horizontal and tilted planes at the earth ' s surface for cloudless atmospheres.Journal of climate and applied meteorology, 1986,25: 87-97.].W is precipitable water vapor (cm), precipitable water vapor w and total atmospheric water steam content W numerically equal, and M is relative atm number, by following formula, is obtained:

M＝[cosθ _z+0.15(93.885-θ _z) ^-1.253] ^-1

θ _zbe solar zenith angle, by camera file, obtained.

Solar distance d ²(astronomical unit) is calculated as follows:

$\mathrm{\Γ}=\frac{360}{365}(\mathrm{dn}-1)$

$d^2=\frac{1}{1.00011+0.034221\cos \mathrm{\Γ}+0.00128\sin \mathrm{\Γ}+0.000719\cos \mathrm{\Γ}+0.00077\sin 2\mathrm{\Γ}}$

Wherein, dn is Julian date, and the image capturing date is apart from the number of days on January 1.

The IRS P6LISS3 image in Yongan City, Fujian Province area on March 24th, 2008 of take is data source, utilizes method of the present invention to calculate atmospheric path radiation (take triband as example), and (unit is wm to result as shown in Figure 2 ^-2sr ^-1μ m ^-1).

Fig. 2 atmospheric path radiation figure.

Claims (1)

1. by pixel, calculate a physical method for satellite remote-sensing image atmospheric path radiation, the steps include:

The first step, satellite image is carried out to geometry correction;

Second step, according to the longitude and latitude of Landsat remote sensing image and video imaging time, the MODIS aerosol optical depth that obtain same geographic range, synchronously passes by (0.47 μ m and 0.66 μ m), total atmospheric water steam content, earth's surface air pressure and ozone concentration image, and each image is done to geometry correction, be converted to the projection pattern consistent with Landsat remote sensing image and resolution, and image ranks number is consistent;

The 3rd step, by total atmospheric water steam content W, calculate atmospheric water steam optical thickness;

${\mathrm{\τ}}_w=\frac{0.2385a_{\mathrm{w\λ}}\mathrm{wM}}{{(1+20.07a_{\mathrm{w\λ}}\mathrm{wM})}^{0.45}}$

Wherein, τ _watmospheric water steam optical thickness, α _{w λ}be water vapor absorption coefficient, w is precipitable water vapor, precipitable water vapor w and total atmospheric water steam content W numerically equal, and M is relative atm number;

The 4th step, calculating Ozone Absorption optical thickness

τ _o＝C _ozone*A _ozone(λ)

τ _oozone Absorption optical thickness, C _ozonefor MODIS ozone concentration (unit is Dobson), A _ozone(λ) be the Ozone Absorption coefficient of wave band λ;

The 5th step, by the aerosol optical depth of MODIS0.47 μ m and 0.66 μ m, utilize following formula to calculate Wavelength Indices (α) and atmospheric turbidity coefficient (β)

$\mathrm{\α}=-\frac{\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_1\right)-\ln {\mathrm{\τ}}_a\left({\mathrm{\λ}}_2\right)}{\ln {\mathrm{\λ}}_1-\ln {\mathrm{\λ}}_2}$

β＝τ _a(λ ₁)λ ₁ ^α＝τ _a(λ ₂)λ ₂ ^α

λ ₁and λ ₂be respectively 0.47 μ m and 0.66 μ m, τ _a(λ ₁) and τ _a(λ ₂) be respectively the aerosol optical depth of 0.47 μ m and 0.66 μ m;

After obtaining α and β, can utilize following formula to obtain the aerosol optical depth τ of any wavelength X _a(λ):

τ _a(λ)＝βλ ^-α

Wherein, λ is wavelength (μ m), τ _a(λ) be the aerosol optical depth of this wavelength;

The 6th step, calculating Rayleigh scattering optics thickness;

${\mathrm{\τ}}_r=\frac{p}{1013.25}0.008569{\mathrm{\λ}}^{-4}(1+0.0133{\mathrm{\λ}}^{-2}+0.00013{\mathrm{\λ}}^{-4})$

Wherein, τ _rfor Rayleigh scattering optics thickness, λ is the centre wavelength (μ m) of each wave band of image; P is the earth's surface air pressure product of MODIS;

The 7th step, calculating shine upon the atmospheric transmittance of direction and sensor observed ray;

T _z＝exp(-τ/cosθ _z)＝exp{(-τ _r-τ _a-τ _o-τ _w)/cosθ _z}

T _v＝exp(-τ/cosθ _v)＝exp{(-τ _r-τ _a-τ _o-τ _w)/cosθ _v}

Wherein, T _zand T _vbe respectively the atmospheric transmittance that shines upon direction and sensor observed ray, τ is atmosphere optical thickness, θ _zsolar zenith angle, θ _vsensor observation zenith angle, τ _r, τ _a, τ _oand τ _wrespectively Rayleigh scattering optics thickness, aerosol optical depth, Ozone Absorption optical thickness and atmospheric water steam optical thickness;

The 8th step, calculating solar distance

Solar distance d ²(astronomical unit) is calculated as follows:

$\mathrm{\Γ}=\frac{360}{365}(\mathrm{dn}-1)$

$d^2=\frac{1}{1.00011+0.034221\cos \mathrm{\Γ}+0.00128\sin \mathrm{\Γ}+0.000719\cos \mathrm{\Γ}+0.00077\sin 2\mathrm{\Γ}}$

Wherein, dn is Julian date, and the image capturing date is apart from the number of days on January 1;

The 9th step, calculating scattering angle

θ _p＝180°-θ _z-θ _v

Wherein, θ _pfor scattering angle;

The tenth step, suppose that atmosphere is that homogeneous atmosphere, scattering angle are greater than 90 while spending, the computing formula of atmospheric path radiation is:

$L_p=\frac{E_0}{d^2}{\mathrm{\τ}}_r\frac{3}{16\mathrm{\π}}(1+{\cos }^2{\mathrm{\θ}}_p)(T_z{T}_v-1)/\ln \left(T_z{T}_v\right)$

Wherein, L _pfor journey radiation, E ₀it is the solar spectrum irradiancy of exoatmosphere respective wavelength.