CN101718866A - Improved physical method for topographic correction of remote sensing images - Google Patents

Abstract

The invention relates to an improved physical method for topographic correction of remote sensing images, which considers the influence of the topography on incident irradiance and earth surface reflectivity received by the sloping surface, more accurately computes sky radiation received by the sloping surface pixel, and meanwhile considers the influence of atmosphere during the topographic correction. The method of the invention can obtain the better result of the topographic correction; and the method belongs to the physical method, has universality and can be applied to various optical remote sensing images.

Description

A kind of improved physical method for topographic correction of remote sensing images

Technical field

The invention belongs to the remote sensing technology field, relate to the processing and the application of remotely-sensed data.

Background technology

In rugged mountain region, landform shows two aspects to the influence of satellite sensor received signal: the one, and the incident irradiance that influence of topography terrain object is received; The 2nd, topographic change the geometry that constituted of the sun, the face of land and satellite sensor three, and then have influence on the reflectivity of terrain object on the satellite sensor direction.Corresponding topographic correction also should comprise two parts, at first the built-up radiation that the surface level pixel is received is converted to the built-up radiation that domatic pixel receives, and (built-up radiation that domatic pixel received comprises three parts: the reflected radiation of beam radia, sky radiation and surrounding terrain to be called the correction of incident irradiance; Accordingly, the incident irradiance is proofreaied and correct and should be comprised that also beam radia is proofreaied and correct, sky radiation is proofreaied and correct and the reflected radiation of surrounding terrain is proofreaied and correct).Be the surface level reflectivity with domatic reflectivity conversion then, be called reflectivity correction.The influence of elimination or minimizing mountain area remote sensing image mesorelief factor is the inevitable requirement that improves the spectral reflectivity inversion accuracy.

Since the seventies in 20th century, Chinese scholars is just launched research at the orographic effect problem of remote sensing image, has proposed many topographic correction of remote sensing images methods, proofreaies and correct as cosine, and C proofreaies and correct, SCS bearing calibration etc.Wherein the most frequently used is the C landform correcting method of lane.But the obvious deficiency of C method is it to be supposed based on lambert's body, only considered that topographic relief receives the influence of incident irradiance to pixel, and it is a kind of empirical model, have subjectivity, the calibration model that obtains for different remote sensing images has nothing in common with each other, do not have versatility, this has just caused difficulty to its widespread use.Development along with the quantitative remote sensing technology, requirement to topographic correction is also more and more higher, this paper has proposed a kind of improved topographic correction of remote sensing images method on the basis of traditional landform correcting method of lane, the advantage of this method is that it has considered the influence of landform to incident irradiance and earth surface reflection rate simultaneously, calculated to degree of precision the sky radiation that domatic pixel received, in topographic correction, considered atmospheric effect, and it is a kind of physical model, have universality, be applicable to all optical remote sensing images.

Summary of the invention

The object of the present invention is to provide a kind of improved physical method for topographic correction of remote sensing images, to overcome the deficiency of existing landform correcting method of lane.Utilize this method can obtain better topographic correction effect, thereby promote the development of association area remote sensing application.

For achieving the above object, the improved topographic correction of remote sensing images method of the present invention's proposition is:

The first step, the sky radiation that utilizes beam radia that 6S radiation delivery Model Calculation surface level pixel received and surface level pixel to be received.

Second step, the beam radia that utilizes the domatic pixel of cosine correction calculation to receive.

The 3rd step, the sky radiation that receives by the domatic pixel of Perez (1990) Model Calculation.

The 4th step, the reflected radiation of calculating surrounding terrain.

The 5th the step, utilize 6S radiation delivery model that remote sensing image is carried out atmospheric correction.

Formula is converted to radiance on the star radiance of atmosphere bottom below the 6th step, the utilization.

$L=\frac{L_{\mathrm{sat}}-L_p-L_{\mathrm{env}}}{\mathrm{\τ}}$

The 7th step, to utilize Dymond (1999) formula be the surface level reflectivity with domatic reflectivity conversion, finishes topographic correction.

Embodiment

1. beam radia calculates

The beam radia that domatic pixel received can be obtained through the cosine calibration shift by the direct solar radiation that its corresponding surface level pixel is received, and concrete conversion formula is as follows:

$E_d=\mathrm{\Θ}{E}_d^h\frac{\cos i}{\cos z}$ (formula 1)

In the formula: E _dBe the beam radia that domatic pixel received, E _d ^hBeam radia for the surface level pixel is received is obtained by the 6S model.I is domatic solar incident angle (angle of direct sunlight and domatic normal), and z is a solar zenith angle.Θ is the topographic shadowing coefficient, if domatic be shadow region (cosi＜0), Θ is 0, otherwise is 1; Domatic solar incident angle i is obtained by formula (2):

(formula 2)

S is the angle of gradient of face of slope,

Be solar azimuth, A is the aspect of face of slope.Solar zenith angle and solar azimuth can be obtained by the image header file, and the domatic angle of gradient and aspect can be obtained by dem data.

2. sky radiation calculates

The sky radiation that the quantitative Analysis face of slope is received and the reflected radiation of surrounding terrain are often relatively more difficult, in a lot of topographic correction models, all do not consider the influence of these two factors.Document [high Yongnian, Zhang Wanchang. topographic correction of remote sensing images progress and comparative experiments thereof. Geographical Study, 2008,27 (2): 467-477.] calibration result to the different terrain bearing calibration compares analysis, found that the reflected radiation of considering sky radiation and surrounding terrain, the topographic correction effect is just better; Consider that not effect is then relatively poor.Therefore in order to improve the correction accuracy of topographic correction model, need the influence of the reflected radiation of consideration sky radiation and surrounding terrain.

Noorian [Noorian A M such as (2008), Moradi I, Kamali G A.Evaluation of 12 models to estimatehourly diffuse irradiation on inclined surfaces.Renewable Energy, 2008,33:1406-1412.] 12 kinds of sky radiation computation models have been carried out comparative study, and utilize the field survey data to verify.Found that Perez model [Perez R, Ineichen P, Seals R, et al.Modeling daylight availability and irradiance componentsfrom direct and global irradiance.Solar Energy, 1990,44 (5): 271-289.] conversion effect optimum, the present invention uses for reference the newest research results of Noorian (2008) etc., utilizes the Perez model to calculate the sky radiation that domatic pixel receives.

$E_f=E_f^h\×\{0.5\×(1-F_1)[1+\cos \left(S\right)]+F_1\×(a/b)+F_2\×\sin \left(S\right)\}$ (formula 3)

A=max (0, cosi) (formula 4)

B=max (0.087, cos z) (formula 5)

E _fBe the sky radiation that domatic pixel receives, E _f ^hBe the sky radiation that the surface level pixel is received, obtain by the 6S modeling.A and b are intermediate variables.S is the angle of gradient of face of slope, and z and i are respectively solar zenith angle and domatic solar incident angle, F ₁And F ₂Be the coefficient that characterizes sky anisotropy degree, they are functions of weather condition.

F ₁=F ₁₁+ F ₁₂* Δ+F ₁₃* z (formula 6)

F ₂=F ₂₁+ F ₂₂* Δ+F ₂₃* z (formula 7)

Δ is a sky brightness, calculates by formula (8),

$\mathrm{\Δ}=E_f^h\×m/E_0$ (formula 8)

E ₀It is the exosphere solar irradiance, m is relative atm number [Kasten F.A new table and approximate formula forrelative optical air mass.Archives of Meteorology Geophysics and Bioclimatology SeriesB, 1966,14 (2): 206-223.]

M=[cos z+0.15 * (93.885-z) ^-1.253] ^-1(formula 9)

F ₁₁, F ₁₂, F ₁₃And F ₂₁, F ₂₂, F ₂₃Table look-up according to sky sharpness (μ) and to obtain, sky sharpness (μ) is pressed following formula and is determined [Perez R, Ineichen P, Seals R, et al.Modeling daylight availability and irradiance components from direct andglobal irradiance.Solar Energy, 1990,44 (5): 271-289.]:

$\mathrm{\μ}=[(E_f^h+E_d^h/\cos z)/E_f^h+1.041\×z^3]/(1+1.041\×z^3)$ (formula 10)

E _d ^hBe the beam radia that is received on the surface level, z is a solar zenith angle.

3. calculate the reflected radiation of surrounding terrain

Reflected radiation (the E of surrounding terrain _a) calculate [Sandmeier S according to following formula, Itten K I.A physically-based modelto correct atmospheric and illumination effects in optical satellite data of rugged terrain.IEEETransactions on Geoscience and Remote Sensing, 1997,35 (3): 708-717.]:

E _a=E ^h* ρ _Mean* V _t(formula 11)

E wherein ^hBe the solar global irradiance (E that surface level received _d ^hAnd E _f ^hSum), ρ _MeanThe average reflectance of expression surrounding terrain, the average reflectance of the image behind the 6S atmospheric correction of learning from else's experience, V _tBe the visible factor of landform, obtain by following formula:

$V_t=\frac{1-\cos \left(S\right)}{2}$ (formula 12)

S is the angle of gradient of face of slope.

4 topographic corrections

Based on earth surface reflection is anisotropic supposition, and topographic relief can cause the variation of the sun-face of land-sensor geometric relationship, thereby has changed the incident angle and the emergence angle of sunray, makes the reflectivity generation marked change of pixel.And for surface level, the incident angle of different pixels is identical with emergence angle, thereby the surface level reflectivity of similar atural object is identical.Therefore, in order to eliminate the influence of topographic relief, need be the surface level reflectivity with domatic reflectivity conversion to the pixel reflectivity.

Surface radiation brightness calculation formula through topographic correction is as follows:

$L=\frac{{\mathrm{\ρ}}_d{E}_d+{\mathrm{\ρ}}_f(E_f+E_d)}{\mathrm{\π}}$ (formula 13)

Wherein, E _d, E _fAnd E _aRepresent the reflected radiation of beam radia, sky radiation and the surrounding terrain of domatic pixel reception respectively.ρ _dBe domatic reflectivity, ρ to direct solar radiation _fBe domatic reflectivity to scattered radiation, L is the radiance of atmosphere bottom, is obtained by following formula:

$L=\frac{L_{\mathrm{sat}}-L_p-L_{\mathrm{env}}}{\mathrm{\τ}}$ (formula 14)

L _SatBe radiance on the star, can obtain L by the sensor radiation calibration _pBe journey radiation, L _EnvBe environmental radiation, τ is an atmosphere optical thickness, and back three can obtain from the 6S model.

[Dymond J R such as Dymond, Shepherd J D.Correction of the topographic effect in remote sensing.IEEE Transactions on Geoscience and Remote Sensing, 1999,37 (5): 2618-2620.] proposed following formula and considered the influence of landform the beam radia reflectivity:

$\frac{{\mathrm{\ρ}}_d^h}{{\mathrm{\ρ}}_d}=\frac{\cos \left(i\right)+\cos \left(e\right)}{\cos \left(i_h\right)+\cos \left(e_h\right)}$ (formula 15)

ρ _dBe domatic reflectivity, ρ to direct solar radiation _d ^hBe the reflectivity of surface level to direct solar radiation, i and e are respectively incident angle and the emergence angle of direct solar radiation on domatic, i _hAnd e _hBe respectively incident angle and the emergence angle of direct solar radiation on surface level.

Because the information major part that satellite sensor received derives from the reflection of the face of land to beam radia, the reflected radiation of sky radiation and surrounding terrain is relative less, and landform is also less to the influence of their reflectivity.In order to simplify calculating, suppose

So:

$L=\frac{{\mathrm{\ρ}}_d^h{E}_d/\mathrm{\δ}+{\mathrm{\ρ}}_d^h(E_f+E_a)}{\mathrm{\π}}$ (formula 16)

$\mathrm{\δ}=\frac{\cos \left(i\right)+\cos \left(e\right)}{\cos \left(i_h\right)+\cos \left(e_h\right)}$ (formula 17)

Final computing formula through the earth surface reflection rate behind the topographic correction is:

${\mathrm{\ρ}}_d^h=\frac{\mathrm{\πL}}{E_d/\mathrm{\δ}+E_f+E_a}$ (formula 18)

The test of 5 topographic corrections

In the topographic correction test below, the inventive method and the C landform correcting method of lane of using always have been carried out comparative analysis.

Test figure is utilized the Landsat 5TM image of Beijing Mountainous Area of North on July 12nd, 2006, and the resolution of testing used DEM is 30 meters.

To the checking of landform calibration result, the present invention uses for reference two kinds of method for quantitatively evaluating commonly used in the practice, i.e. " correlativity of earth surface reflection rate and corresponding solar incident angle cosine " and " standard deviation ".Before the topographic correction, the earth surface reflection rate of pixel and the correlativity of solar incident angle cosine are obvious, and pixel earth surface reflection rate changes along with the variation of solar incident angle significantly.Through behind the topographic correction, if the influence that the solar incident angle that the earth surface reflection rate of pixel no longer is subjected to cause because of topographic relief changes, correlativity between them should weaken so, and correlativity is weak more, shows that the effect of topographic correction is good more.After passing through topographic correction in addition, the earth surface reflection rate value of similar atural object is more approaching, and atural object becomes more " homogeneous ", so standard deviation should reduce, and standard deviation is more little, shows that the effect of topographic correction is good more

In order to calculate the correlativity of earth surface reflection rate and corresponding solar incident angle cosine, with earth surface reflection rate image and the cosi image overlay before and after the topographic correction, picked at random sampled point then, then the slope of these sampled point equations of linear regression and related coefficient are represented the correlativity of earth surface reflection rate and cosi.Slope is big more, correlativity strong more (referring to that the earth surface reflection rate is subjected to influence that cosi changes greatly); Related coefficient is big more, and correlativity is strong more.

The correlativity (is example with TM the 4th wave band) as shown in table 1 of each face of land reflectivity and solar incident angle cosine cosi before and after the topographic correction:

The correlativity of each face of land reflectivity and cosi before and after table 1 topographic correction

	The slope of linear equation	Related coefficient
			Before the topographic correction	??0.1415	??0.61
After the inventive method is proofreaied and correct	??-0.0008	??0.004
			After the C bearing calibration is proofreaied and correct	??-0.0511	??0.232

As can be seen from Table 1, utilize the inventive method to carry out topographic correction after, the correlativity between earth surface reflection rate and the cosi (slope and related coefficient) obviously weakens, and has shown the validity of the inventive method.After utilizing the C bearing calibration to proofread and correct, the correlativity between earth surface reflection rate and the cosi also obviously weakens, but slope and related coefficient show that greater than the inventive method the calibration result of the inventive method is better than the C bearing calibration.

In order further to verify the effect of topographic correction, on the TM image, choose forest sample district, add up the standard deviation (table 2) of topographic correction front and back TM the 4th wave band reflectivity then respectively.

Standard deviation relatively before and after table 2 topographic correction

	Standard deviation
		Before the topographic correction	??6.99
After the inventive method is proofreaied and correct	??4.38
		After the C bearing calibration is proofreaied and correct	??5.01

A good landform correcting method of lane makes that the standard deviation of similar atural object reduces behind the topographic correction, the inventive method gained result has just in time verified this point, through behind the topographic correction, standard deviation obviously reduces, the inside that shows similar atural object diminishes, the actual state that more can reflect the face of land, and the standard deviation of the inventive method correspondence shows that less than the C bearing calibration calibration result of the inventive method is better than the C bearing calibration.

Claims (5)

1. an improved physical method for topographic correction of remote sensing images the steps include:

The first step, the sky radiation that utilizes beam radia that 6S radiation delivery Model Calculation surface level pixel received and surface level pixel to be received;

Second step, the beam radia that utilizes the domatic pixel of cosine correction calculation to receive;

The 3rd step, the sky radiation that receives by the domatic pixel of Perez (1990) Model Calculation;

The 4th step, the reflected radiation of calculating surrounding terrain;

The 5th the step, utilize 6S radiation delivery model that remote sensing image is carried out atmospheric correction;

Formula is converted to radiance on the star radiance of atmosphere bottom below the 6th step, the utilization;

$L=\frac{L_{\mathrm{sat}}-L_p-L_{\mathrm{env}}}{\mathrm{\τ}}$

The 7th step, to utilize Dymond (1999) formula be the surface level reflectivity with domatic reflectivity conversion, finishes topographic correction.

2. physical method for topographic correction as claimed in claim 1, wherein the concrete conversion formula in second step is as follows:

$E_d=\mathrm{\Θ}{E}_d^h\frac{\cos i}{\cos z}$

In the formula: E _dBe the beam radia that domatic pixel received, E _d ^hBeam radia for the surface level pixel is received is obtained by the 6S model; I is domatic solar incident angle, and z is a solar zenith angle, and Θ is the topographic shadowing coefficient, if domatic be the shadow region, Θ is 0, otherwise is 1; Domatic solar incident angle i is obtained by following formula:

S is the angle of gradient of face of slope,

Be solar azimuth, A is the aspect of face of slope, and solar zenith angle and solar azimuth can be obtained by the image header file, and the domatic angle of gradient and aspect can be obtained by dem data.

3. physical method for topographic correction as claimed in claim 2, wherein the detailed calculated process in the 3rd step is:

$E_f=E_f^h\×\{0.5\×(1-F_1)[1+\cos \left(S\right)]+F_1(a/b)+F_2\×\sin \left(S\right)\}$

a＝max(0，cosi)

b＝max(0.087，cosz)

E _fBe the sky radiation that domatic pixel receives, E _f ^hBe the sky radiation that the surface level pixel is received, obtained by the 6S modeling that a and b are intermediate variables, S is the angle of gradient of face of slope, and z and i are respectively solar zenith angle and domatic solar incident angle, F ₁And F ₂Be the coefficient that characterizes sky anisotropy degree, they are functions of weather condition;

F ₁＝F ₁₁+F ₁₂×Δ+F ₁₃×z

F ₂＝F ₂₁+F ₂₂×Δ+F ₂₃×z

Δ is a sky brightness, calculates by following formula,

$\mathrm{\Δ}=E_f^h\×m/E_0$

E ₀Be the exosphere solar irradiance, m is relative atm number,

m＝[cos?z+0.15×(93.885-z) ^-1.253] ^-1

F ₁₁, F ₁₂, F ₁₃And F ₂₁, F ₂₂, F ₂₃Tabling look-up according to sky sharpness (μ) obtains, and sky sharpness (μ) is pressed following formula and determined:

$\mathrm{\μ}=[(E_f^h+E_d^h/\cos z)/E_f^h+1.041\×z^3]/(1+1.041\×z^3)$

E _d ^hBe the beam radia that is received on the surface level, z is a solar zenith angle.

4. physical method for topographic correction as claimed in claim 3, wherein the concrete computing formula in the 4th step is as follows:

E _a＝E ^h×ρ _mean×V _t

E wherein _aThe reflected radiation of the surrounding terrain that receives for domatic pixel, E ^hSolar global irradiance (the E that is received for the surface level pixel _d ^hAnd E _f ^hSum), ρ _MeanThe average reflectance of expression surrounding terrain, the average reflectance of the image behind the 6S atmospheric correction of learning from else's experience, V _tBe the visible factor of landform, obtain by following formula:

$V_t=\frac{1-\cos \left(S\right)}{2}$

S is the angle of gradient of face of slope.

5. physical method for topographic correction as claimed in claim 4, wherein the detailed calculated process in the 7th step is:

$L=\frac{{\mathrm{\ρ}}_d{E}_d+{\mathrm{\ρ}}_f(E_f+E_a)}{\mathrm{\π}}$

Wherein, E _d, E _fAnd E _aRepresent the reflected radiation of beam radia, sky radiation and the surrounding terrain of domatic pixel reception respectively, ρ _dBe domatic reflectivity, ρ to direct solar radiation _fBe domatic reflectivity to scattered radiation, L is the radiance of atmosphere bottom, is obtained by following formula:

$L=\frac{L_{\mathrm{sat}}-L_p-L_{\mathrm{env}}}{\mathrm{\τ}}$

L _SatBe radiance on the star, can obtain L by the sensor radiation calibration _pBe journey radiation, L _EnvBe environmental radiation, τ is an atmosphere optical thickness, and back three can obtain from the 6S model;

Consider the influence of landform according to following formula to the beam radia reflectivity:

$\frac{{\mathrm{\ρ}}_d^h}{{\mathrm{\ρ}}_d}=\frac{\cos \left(i\right)+\cos \left(e\right)}{\cos \left(i_h\right)+\cos \left(e_h\right)}$

ρ _dBe domatic reflectivity, ρ to direct solar radiation _d ^hBe the reflectivity of surface level to direct solar radiation, i and e are respectively incident angle and the emergence angle of direct solar radiation on domatic, i _hAnd e _hBe respectively incident angle and the emergence angle of direct solar radiation on surface level.

Suppose So:

$L=\frac{{\mathrm{\ρ}}_d^h{E}_d/\mathrm{\δ}+{\mathrm{\ρ}}_d^h(E_f+E_a)}{\mathrm{\π}}$

$\mathrm{\δ}=\frac{\cos \left(i\right)+\cos \left(e\right)}{\cos \left(i_h\right)+\cos \left(e_h\right)}$

Final computing formula through the earth surface reflection rate behind the topographic correction is:

${\mathrm{\ρ}}_d^h=\frac{\mathrm{\πL}}{E_d/\mathrm{\δ}+E_f+E_a}$