GISdevelopment.net ---> AARS ---> ACRS 1990 ---> Poster Session Q

Atmospheric and topographic correction of multitemporal Landsat TM imagery

Zhang Yugui
The Research Institute of Forest resource Information
Chinese Academy of Forestry


Abstract
The paper briefly introduced a developing system for comprehensive correction satellite imagery. Standarized reflectance images were produced by atmospheric correction, topographic correction, sensor calibration etc. The atmospheric parameters parameters were interpolated by using local meteorological records. The path radiance derived from visual range is adjusted with dark pixels in the image data. A digital topographic model was used both to correct for the georadiometric effects and to monitor the optical thickness. The normalized reflectance image were enhanced with linear or non-linear transformation. A set of multitemporal subscenes of four dates in Trolleba, Sweden were used as examples to illustrate the good correspondence of correction results in quite different atmospheric conditions.

Introduction
The main factors that seriously affect the radiometric quality of digital satellite data are atmospheric attenuation and topographic distortion, in addition to sensor variation. All these factors are related to and restrictive with each other. Therefore, a correction algorithm that considers only a single digital topographic model. The paper described a developing correction system which made comprehensive correction of Landsat TM and other satellite data, including sensor calibration, path radiance evaluation, the calculation of attenuation coefficients of atmosphere, the variations of direct solar illumination and diffuse illumination caused by topography, as well as minor modifications of the components of optical thickness due to calculate the needed atmospheric parameters by using metrological records, and presented result ofthe calculation after these parameters were input to the radiative transfer model of Turner.

The multitemporal TM images used in this study consist of 4 subscenes of 400 X 400 pixels of the same Trollebo scene, located in south-central part of sweden. The elevation of the area is from 164 to 285 meters ASL. Clear water bodies are scattered over the whole area. The trollebo area contains mainly forest lands with a small portion of agricultural fields. The four dates constitute the whole growing season of forest vegetation , Pine, spruce , birch and broadleaves make up the major forest type categories. Table 1 listed the data of the 4 dates.

Table 1. Parameters for trollebo subscence. The scene center coordinates,
sun elevation azimuth and view angles are in decimal degrees.
Scene Over pass time scene center coordinates sun azim sun elev view angle Sun-earth dist. (au)
1 194/20:2 1984.4.27 09 35 39 (57.317,14.756) 133.42 39.17 2.058 0.9932
2 194/20:2 1984.5.13 09 36 12 (57.819,14.718) 131.47 43.47 2.733 0.9794
3 193/20 : 1984.7.09 09 27 05 (57.317,15.713) 123.88 44.98 1.324 0.9670
4 193/20:2 1984.9.0209 39 38 (57.317,14.493) 136.92 33.66 2.987 0.9827

The correction procedures
Conceptually, the reflectance images can be obtained through three sequential linear transformations.
  1. A Transformation From Digital Numbers (DN) To observed Spectral Radiance Values (Lobs)
    This is accomplished by using published radiance calibration constants (Barker 1987), which is usually termed absolute calibration. Spectral radiance and irradiance were used in the study instead of in-band units throughout the computation procedures to facilitate comparison between sensor (Price 1987).

    The Lobs can be obtained by

    Lobs = DN* (Lmax - Lmin) / 255 +Lmin ----------------(1)

    Lmax - Lmin are constants representing the radiance at DN= 0 and DN=255, respectively. These constants may vary with date, ground receiving and processing station. Before this tranformation, a relative compensation calobration between etectors were performed for most scenes.

  2. A Transformation From Radiance Lobs to Spectral Reflectance R
    It is complicated and critical. The assumption of a Lambertian scattering behaviour for most targets would simplify the view angle effect as only differential atmospheric effect. Therefore, the correction of main interest is atmospheric and topographic corrections, illumination and view angle corrections, in addition to sensor calibration. The radiance observed by a sensor in a given band can be expressed as

    Lobs = 1 / p (Eo cosi Exp(-t/cosqo) +gEd)R Exp(-t/cosqv ) +Lp ------(2)

    or in a simplified form

    Lobs = 1/ p R Eg Tv +Lp

    and the transformed reflectance R for a target can be expressed as

    R= p (Lobs- Lp ) / Eg Tv) ..............(3)

    The meanings of various factors are defined as follows:

    Lobs = radiance observed at the sensor (mWcm-2 Sr-1)

    1 / p = The bidirectional reflectance distribution function (BRDF) for a Lambertian surface

    Eo= solar irradiance at the top of the atmosphere (m Wcm-2),

    Eo was calibrated with published solar exoatmospheric spectral irradiances in TM bands

    Eo was calibrated with published solar exoatmospheric spectral irradiances in TM bands

    i= solar incidence angle with respect to surface normal

    t= atmospheric extinction coefficient (optical thickness)

    qo, qv = solar zenith angle and observation zenith angle

    g = a geometric factor taking into account of the part of sky radiance not contributing in the illumination of inclined terrain , g= 0.5+0.5coss, s=terrain slope angle

    Ed = diffuse sky irradiance on a horizontal surface (m Wcm-2) Ed is modified by atmosphere and topography.

    Eg= (Eo cosi Tq+g Ed) = global irradiance on the surface (m Wcm-2)

    Tq = exp(-t/cosqo) = atmospheric downward transmittance

    Tv = Exp (-t/ cosqv) = atmospheric up eard transmittance

    Lp = path radiance from multiple scattering (m Wcm-2 sr-1 )

    This reflectance tranformation involves the approximation of some varibles by using different methods according to parameters available . Combining equation (1) and (3) we obtain

    R= p / (Eg Tv) *(DN(Lmax - Lmin) / 255 +Lmin -Lp) -----(3a)

  3. A Transformation of R to an Eight-bit Digital Number (DN')
    It may be a linear or non-linear transformation, incorporating appropriate contrast stretch between predefined minimum and maximum reflectance limits Rmax - and Rmin

    DN' = 255 F[(R-Rmin) / (Rmax - Rmin)] ................(4)

    The linear transformation is convenient if the enhanced reflectance images are to undergo further processing, but for the hard copy products, a non-linear transformation can be applied to TM bands. The square root function increases the contrast of low reflectance portion, and compress the contrast of high reflectance portion (Ahern 1989).

    Although each step involves much consideration and techniques, the second step requires more investigation and exploration. The actual scattering properties of natural surfaces differ from the Lambertian assumption to varying extent (Teillet et al. 1982) The diffuse approximation is more likely to be valid for terrain slopes of less than 25 degrees and effective illumination angles of less than 45 degrees (Smith et al. 1980) or further less (Cavayas 1987). In our study, the solar zenith angle exceeded these limits, which may cause a system error of the reflectance values. Since we have no sufficient knowledge about the BRDF for natural targets, such assumption is the only convenient solution at present, espescially for forest canopy complexities which are considered to be Lambertian surface in small solar zenith angle.
Evaluation and computation of atmospheric parameters
The global irradiance Eg, diffuse irradiance Ed and path radiance Lp for the center wavelength of TM bands were calculated with the Turner and Spencer (1972) reductive transfer model, and with Forster's modification (1984) for the computation of the single scattering phase function, after appropriate TM band adjustment. This requires estimated or calculated optical thickness t,solar zenith qo solar azimuth, nadir view angle qv , view azimuth as input. The scene center coordinates and scene center overpass time were used to derive the components of optical thickness t. The radiative transfer model provided convenience for establishing relationship between atmospheric parameters. In the process of computing the parameters, the scale heights for Ed, Eg and Lp were determined, which facilitate the computation procedures of corresponding values for each pixel.
  1. The Estimation of Path Radiance Lp
    Lp is an important factor. In our study, it was derived from radiative transfer model after computation of the input meteorological data. However, several approaches can be applied to evaluate Lp for the purpose of reference. Since clear water bodies exist in the studied area, the dark pixel method (Ahern 1977) was used due to its simplicity . The relative scattering function method (chavez 1987 , 1988) is also convenient . The covariance matrix method (Switzer 1981) and RIM ( Crippen 1987) require that spectrally homogeneous targets being distributed over slopes of various aspects. In this study, the Lp values were used to compare with those calculated by using meteorological data, and to monitor and adjust the latter if they were unreliable.

  2. Use of Local metrological records to Estimate optical Thickness
    The optical thickness t in equation (1) is a critical coefficient and is affected by many factors. since the measurement of at overpass time is seldom done, it is usually estimate some atmospheric various methods meteorological records could be utilized to estimate some atmospheric parameters. The daily records at 1,7,13 and 19 hours include temperature, air pressures, water vapour pressure, relative humidity, visibility, wind speed etc. The overpass time data were derived from the interpolated data. The physical atmospheric thickness was accounted for as function of topographic elevation. There fore the station elevation ASL (above sea level) and the parameter values at the station elevation were used to derive the scale heights of parameters, which were applied to each pixel of the image when a topographic model was available. The normal optical thickness t equals the sum of the attenading coefficients, and they are independent. Therefore t is made up of Rayleigh component tr due to molecular scattering, aerosol component tp due to particle scattering, and components due to selective absortption of gases, such as water wapour component tw, ozone component toz, and carbon dioxide component tco2. All these are wavelength dependent.

    t = tr + tp + tw + toz + tco2 +.......(5)

    we deal with the t components separately. The Rayeigh molecular scattering is regarded as of dry air, and the water absorption as of partial pressure of water vapour, and the aersol scattering as of suspending particles.

    The visibility may be a measure of tp for a particular vertical distribution . Turner (1972) graphed tp as a function of visibility for continental atmosphere. We have interpolated and extrapolated the graphs for TM bands with the central wavelength as reference, and then an adjustment or alignment can be applied to the interpolated tp values with shaded clear water pixel values in TM band 4 as reference in the process of calculation.
The double usage of topographic model
The inclusion of even a rough elevation data could decisively help to stress the difference of the physical thickness of the atmosphere and consequently help to adjust the optical thickness components tr tp and tw for each pixel as well as to monitor the changes of the scaleheights of Eg Ed and Lp with the elevation of each pixel. A course topographic model as we have used would help to modify and improve the enhanced reflectance images to 3-4 counts.

The coarseness of spatial and quantization resolution of the DTM data may not help to raise the correction accuracy for rugged mountain area, since the slope the slope-aspect geometry is not a Gaussian distribution (Teillet 1986). The digital elevation data used in this experiment were created by digitizing contour lines of the UTM map sheet at 1:50 000 scale and linearly interpolated between contour lines. The terrain slope s and aspect a can be derived with reference to the neibouring pixels. The solar incidence angle to the surface normal is expressed by:

cosi =cosqo coss +sin qosincos(fo-a)

s= terrain slope , a = terrain aspect (azimuth)

There are three optional applications for the elevation data in the algorithms:
  1. In absence of elevation data the average elevation of the subscene as well the elevation of the reference water bodies are used.

  2. when a coarse elevation model is available and the topography is rather complicated used for the adjustment of the optical thickness and illumination conditions for each pixel.

  3. when accurate elevation data could be obtained or the areas of internet are not so complicated as to demand an accurate topographic model it is worth to transfer the data in to a slope aspect format and both the optical thickness and illumination geometry expressed by cosi and g could be considered in our study after such correction the mean square of spectral features of the some targets were reduced by about 30%.
Result and discussion
The aerosol content, its thickness and distribution atmosphere are the main problem atmospheric correction .The visibility codes 1,2,4 indicates the stability of atmosphere at the time of the observation and the tp and Lp values hence obtained were quite reliable but the adjustment tp are over estimated the zero reflectance pixels were not exactly found because of negligence of wind speed and diffuse irradiance on the water interpreted due to partial shading of the pixels date 3 had pixel values approached and the interpreted visibility is only 7.5 Km the adjustment of dark pixel values approached the real value and improved the result for a clear stable atmosphere the visual range codes seemed to be more reliable but if the recorded visual ranges changed abruptly the dark pixel helped to approach the real value than interpolated value.

Table 2. The extroatmospheric solar irradiance Eo are sin earth distance adjusted. The global spectral irradiance
Eg diffuse spectral irradiance Ed and spectral path radiance Lp* are the values at sea level the scale
heights of Eg* Ed Lp are also justed with the pixel respectively
Parameter TM1 TM2 TM3 TM4 TM5 TM7
date 1 (4/27)
Eo(mWcm-2um-1) 193.21 180.57 153.72 103.36 21.65 7.35
Eg(mWcm-2um-1) 106.03 102.49 91.11 62.61 12.81 4.45
Ed(mWcm-2um-1) 27.03 21.49 16.22 9.24 1.30 0.38
Lp*(**um-1sr-1) 1.1814 0.4247 0.1811 0.1213 0.0099 0.0019
Lp*(**um-1sr-1) 1.9415 0.7823 0.3619 0.2623 0.0261 0.0052
Egscale(m) -38569 -50997 -65349 -89163 -101761 -127130
Edscale(m) 1754 1542 1439 1366 1283 1267
Lp Scale(m) 1482 1384 1341 1395 1287 1265
date 2 (5/13)
Eo(mWcm-2um-1) 191.67 179.14 152.49 102.54 21.47 7.49
Eg (mWcm-2um-1) 115.77 111.43 98.80 67.38 13.61 4.76
Ed (mWcm-2um-1) 31.55 24.93 18.68 10.53 1.46 0.43
Lp(**um-1sr-1) 1.7235 0.6813 0.3016 0.2514 0.0208 0.0039
Lp*(**um-1sr-1) 2.1155 0.8612 0.3887 0.3241 0.0277 0.0053
Egs scale (m) -40923 -53979 -68640 -92234 -91231 -117572
Edscale (m) 1806 1576 1461 1381 1289 1271
Lp Scale(m) 1459 1372 1406 1406 1294 1270
date 3 (7/9)
Eo(mWcm-2um-1) 189.23 176.86 150.55 101.24 21.20 7.20
EgmWcm-2um-1) 117.20 112.71 99.85 65.76 12.37 4.50
EdmWcm-2um-1) 39.74 31.85 23.89 13.36 1.72 0.53
Lp(**um-1sr-1) 5.4500 2.1598 0.8805 0.9670 0.0480 0.0071
Lp*(**um-1sr-1) 2.6805 1.0984 0.4522 0.5544 0.0269 0.0039
Egscale(m) -36276 -46297 -56874 -69934 -49391 -66263
Edscale(m) 1837 1609 1492 1424 1328 13.01
Lp Scale(m) 1302 1267 1228 1424 1319 1278
date 4 (9/2)
Eo (mWcm-2um-1) 192.31 179.73 153.00 102.89 21.55 7.32
Eg (mWcm-2um-1) 91.18 88.86 79.27 52.92 10.20 3.67
Ed (mWcm-2um-1 17.01 13.43 10.12 5.81 0.74 0.22
Lp (**um-1sr-1 1.1681 0.4255 0.1659 0.1785 0.0077 0.0011
Lp* (**um-1sr-1) 1.4018 0.5345 0.2173 0.2373 0.0110 0.0015
Eg scale(m) -39656 -54371 -70899 -83554 -63013 -86801
Ed scale(m) 1661 1465 1375 1348 1272 1255
Lp scale(m) 1670 1509 1392 1457 1286 1254

Table 2 listed the sea level values of Ed and Eg and Lp and their scale heights diffuse spectral irradiance Ed and spectral path radiance Lp will decrease with elevation while the global spectral irradiance Eg increases with elevation. The introduction of scale height for these three parameters simplified the computation of their modification due to elevation variation it has been examined that error committed by using scale height method is only one per thousand compared with applying more complicated algorithms.

Clear lakes of relative stable reflectance were selected for comparison of correction result the obvious atmospheric between original counts were greatly compressed(Table 3).

Table 3. Comparison of original values and corrected reflectance for relative stable object (clear lakes)
  TM1 TM2 TM3 TM4 TM5 TM7
average original counts
date 1 53.8 16.6 11.5 9.2 4.7 3.6
date 2 57.7 18.2 13.7 10.6 6.1 3.8
date3 63.5 19.7 15.9 11.4 7.2 3.9
date 4 41.1 13.5 9.4 6.3 3 3.0
Corrected reflectance
date 1 0.0539 0.425 0.0266 0.0000 0.0001 0.0000
date 2 0.533 0.414 0.0266 0.0000 0. 0000 0.0000
date 3 0.549 0.0403 0.0259 0.0000 0.0000 0.0000
date 4 0.0556 0.0401 0.273 0.0000 0.0000 0.0000

The following conclusion can be derived from this study:
  1. The comprehensive correction to remotely sensed data for atmospheric and georadiometric distortion is feasible the converted reflectance stable reflectance objects showed an excellent conformity how ever the degree of correction accuracy on the precision of the input data if a even radioactive transfer model and a precise correction equation are used.

  2. The meteorological record could be used to estimate the components of optical thickness and to make clear distribution of both direct solar and diffuse irradiance the path radiance Lp should be estimated by using multiple methods for comparison since the visibility observation in comparison among transfer models is also beneficial.

  3. A georadiometric correction may be an option if a reliable digital terrain model is not available such a model could be used for both adjusting the topography dependent atmospheric coefficients and for correcting the surface illumination conditions.

  4. For the enhancements of reflectance images the forest reflectance limits fora biophysical areas need to be recertified from sufficient TM scenes and added to the procedure the correct decision of the reflectance limits would guarantee color balance color products.

  5. Only for a solar zenith angle of less than 33 degree or so the topographic effect can be removed by including a slope model relative standardization of reflectance images can be expected at present due to the estimation and approximation nature of the parameters and due to the gap between the prerequisites for an accurate and the data available.
References
  • Abern , F.J et al.. 1989 Reflectance enhancements for the TM images Ph. Eng. Rem. Sen. V-55 No.1 pp.- 61-67.
  • Chavez P.S. 1989 Radiometric Calibration of land sat TM images Ph. Eng. Rem Sen V.55 No 9 pp. 1285 -1294.
  • Civco D.L 1989 Topographic normalization of land sat TM imagery Ph Eng Rem. Sen. V.55 No.9 pp 1393 1309.
  • Crippen R..E.. 1987 The RIM of adjusting images data for band rationing int.. J.R.S. V.8 No.2 pp. 137-155.
  • Forester B.C. 1984 Derivation of atmosphere correction procedures Int. J.r.S. V.5 no.5 pp. 799 817
  • Kawata Y. et.. al. 1987 Radiometric correction for atm. and topogr. effects on MSS images int. J.R.S. pp 730 748.
  • Markham, B.L. and Barker, J. L., 1987, TM bandpass solar irradiances, Int. J. R. S., V.8, No. 3, pp.517-523
  • Markham, B. L. et al., 1985, Spectral characterization of TM sensors, Int. J. R. S., V.6, No. 5, pp. 697--716
  • Woodham, R. J. et al., 1987, An analytic method for radiometric correction, IEEE transactions on geoscience and R. s., V.GE-25, No.3.
  • Zhang Y. G., 1989, Normalization of Landsat TM imagery by atm. Correction, Report 1, R. S. lab., Swedish Univ. Agr. Science.