Methods

   Landsat 5 Thematic Mapper and Landsat 7 Enhanced Thematic Mapper imagery were collected for a small area (~121 km²) in northern British Columbia. The satellite scenes were captured on: 20 July 1985, 24 August 1992, 29 July 1994, 23 September 1997, 4 August 1999, 13 September 1999, 23 September 2000, 22 July 2003. The 4 August 1999, and 23 September 2000 scenes are the only Landsat 7 Enhanced Thematic Mapper scenes. The imagery was image to image georectified to a single reference image (the 1997 image).

   Two of the radiometric correction methods required a working definition for what would constitute a significant digital number for this study. A significant digital number is a digital number that occurs in such proportion within an image that it is highly unlikely to be an outlier, which is an extreme deviation from the mean. In previous studies, a significant digital number was arbitrarily defined as any digital number having a pixel count of at least 1000 pixels (Tiellet and Fedosejevs, 1995; McDonald et al., 1998; Song et al., 2001). However, this does not account for the digital number spread of the dataset. In this study a significant digital number was mathematically defined by equation 1, i.e. pixel count was defined as an inverse proportion of the spread of the data.

   The value for σij uses old-growth spruce stands because it has been found that spruce dominated stands in regions similar to the study area reach an old-growth stage where growth is equal to mortality, and therefore most stand structural attributes are in equilibrium provided the system is not changed by any kind of disturbance (Kneeshaw and Burton, 1998). British Columbia’s Ministry of Forests working definition of old-growth for any spruce stand in the Sub-Boreal Spruce Biogeoclimatic zone, is any spruce stand which has not experienced a stand-initiating event in 150 years or more. In addition, it is noted that the older the stand the more likely it is to be an old-growth stand and that coupling a minimum stand height of 30m with age may be the best possible working definition (Kneeshaw and Burton, 1998; Wells et al., 1998). Using a forest cover data set, the spruce stands in the study area that were greater than 210 years of age and greater than 30m in height for all imagery in the study were identified as old-growth spruce stands. The underlying assumption for the use of old-growth spruce stands is, in terms of Landsat sensors, that old-growth spruce stands do not change in reflectance over time.

   Pj was defined such that the resulting mean PCij for all imagery of a given band was equal to 0.1% of the entire image. Significant pixel counts were defined in this way to ensure that imagery with a much higher range of digital numbers had a PCij inversely proportional to that range. Therefore, in this study, PCij values were lower for imagery with a greater range of data, and thus, the range of significant digital numbers in any image was more highly correlated to the true range of digital numbers than was the case if PCij was defined as simply equal to Pj. In any case, Pj should be defined such that imagery outliers are removed but a healthy proportion of the true range of the imagery is maintained. For the remainder of this paper, a significant digital number is a digital number that meets a minimum pixel count of PCij.

ATCOR1
   The majority of absolute radiometric correction methods rely on data that was not available for this time series. The algorithms available from PCI Geomatics for absolute radiometric correction (ATCOR1) were used because the only ancillary data required is the solar zenith angle of each image and the location of old-growth spruce stands.

   The atmospheric models requires the selection of atmospheric properties. The atmospheric properties are predefined and are either tropical, mid-latitude, or the US standard atmosphere and are also either rural, urban, desert, or maritime. For this area the US standard rural atmosphere best described the study area. The atmospheric model uses atmospheric properties, sensor calibration defaults, and solar zenith angles to calculate reflectance values. A variable, optical visibility, has to be calculated for each image to perform the final algorithms in the atmospheric correction package. To derive the optical visibility, the algorithm requires the selection of a target with known reflectance values, which are compared to the reflectance values calculated within each image. For this study the known target was old-growth spruce. The optical visibility is adjusted to account for the variation between known reflectance and calculated reflectance. These algorithms are based on the Richter atmospheric model (Richter, 1990). In the remainder of this paper this radiometric correction method will be referred to as ATCOR1. The calculation to determine optical visibility requires the sensor gain and offsets for each TM band. In the majority of archived data, this is not known and therefore established averages for sensor gains and offsets were used for this study. The gains and offsets could have been adjusted to better match the known reflectance of old-growth spruce stands to the calculated reflectance. This was not done for two reasons. First, the adjustment of gains and offsets would introduce subjectivity by matching calculated reflectance to the desired known reflectance, and not necessarily to an accurate reflectance. Second, the adjustment of gains and offsets would have greatly increased user processing time and detract from one purpose of this study, which was to test or develop an automated or semi-automated radiometric correction method.

ATCOR2
   The atmospheric correction package also allows for the adjacency effect, which is the effect that radiance from neighbouring pixels may have on the measured radiance of a target pixel. To account for the adjacency effect, an average filter is passed over the results from ATCOR1. The filtered data is incorporated into the results from ATCOR1 to minimize the adjacency effect. This form of radiometric correction will be referred to as ATCOR2. A systematic or software error altered the results of ATCOR2 on TM band 7. This error resulted in only four of the eight images having meaningful results, the other four having a large proportion of the minimum and maximum digital numbers, 0 and 255. Results for ATCOR2 on TM band 7 incorporate only those four channels with meaningful output.

OFS
   The origin and range for each set of TM bands in the imagery were found to be quite variable. The significant origin for each band was simply the minimum significant digital number in the band. The range was defined as the maximum significant digital number minus the minimum significant digital number. The true origin and range of digital numbers were not used because extremely bright and dark outlier pixels were present in some images but not in others. Table 1 shows the significant origin and range and the true origin and range for each TM band 4 in the image series. The final column in the figure shows the true range divided by the significant range, indicating the variation in the imagery of outlier pixels. The image data was transformed using equation 2.

Table 1. True origin and range, and significant origin and range for Landsat Thematic Mapper band 4.

   Rmax j and Amax j are used for quality control to ensure that the resulting imagery does not get compressed and does not contain values less than zero. Before each band is transformed an inspection is made to ensure that equation 2 does not result in the original maximum pixel value for each band of imagery being scaled to greater than 255. In this study this was not the case for any band of imagery, but if it were, the output should be to 32-bit channels. For the remainder of this paper, this radiometric correction method will be referred to as origin fix with scaling (OFS).

PIFR
   A method for non-subjective PIF selection was developed and showed a significant decrease in radiometric inconsistency (Du et al., 2002; Du et al., 2001; Song et al., 2001). Consider the scatterplot of TM band 4 for two images (Fig. 1). The major axis is the solid black line and the two dotted parallel lines are thresholds defined by a deviation l from the major axis. For this method, the variations in pixel values during the period represented by the scatterplot are assumed to be linear, spatially homogenous, and normally distributed. All pixels that fall within the threshold are considered PIFs for that image pair. For each TM band a scatterplot was created using the 23 September 1997 image as the y axis and every other image as an x axis in a separate scatterplot. The time series has eight images, so seven comparisons with the 23 September 1997 image were possible. Any pixel that fell within the thresholds for every scatterplot is considered a PIF. Thresholds were determined by ensuring that the scatterplot of each image pair under the resultant PIFs had a correlation coefficient greater that 0.9. If the correlation coefficient was less than 0.9 the deviation l from the major axis was adjusted to determine new thresholds.

Figure 1. Scatterplot of Landsat Thematic Mapper band 4 for two images, showing the major axis and the deviation l from the major axis to determine thresholds.

   Each TM band has a single dataset describing the PIFs for every image in that band. For each TM band, the appropriate PIF dataset is used to calculate the standard deviation and mean for the pixel values for every image. The standard deviations and means for TM band 4, as well as the gains and offsets used in the image transformations are shown in table 2.

   In table 2, gain for any given image is defined as the maximum standard deviation for all imagery for a given TM band divided by the standard deviation of the TM band for that given image. This calculation for gain ensures that the gain is greater than or equal to one, which ensures the data is stretched and not compressed. Sref is defined as the gain multiplied by the mean. The offset for any given image is defined as the maximum Sref for all imagery of a given TM band minus the Sref of the TM band for that given image. This calculation for offset ensures that the offset is greater than or equal to zero which ensures that the output data cannot be negative. Due to the large variation in the standard deviation and mean of the imagery under the PIFs in this study, the resulting gain and offset also have large variation. When transformed using the gain and offset, the imagery falls outside of the eight bit range (0-255) and therefore 32 bit storage was used for the transformation output. The output was then scaled as 32 bit data to match the 0-255 range of the other radiometric correction methods. This method of radiometric correction will be referred to as pseudo invariant features regression (PIFR).

Table 2. Pseudo-invariant features statistics for Landsat Thematic Mapper band 4, used to determine gains and offsets for PIFR.

MBDS
   The large variation in gain and offset from the PIFR method resulted in the minimum significant digital number and maximum significant digital number being highly variable across the imagery for a given TM band using PIFR radiometric correction. To account for this, median based directional scaling was performed on the data set. The median was chosen over the mean as the median had slightly less variation for all imagery of a given band and is less sensitive to extreme outliers. For each channel, two range calculations were performed, the median minus the minimum significant digital number and the maximum significant digital number minus the median. Any pixel less than the median was adjusted using equation 3 and any pixel greater than or equal to the median was adjusted using equation 4. The data was already 32 bit, so negative numbers were allowable as the end result would be scaled to 0-255 as 32 bit data.

   This method of radiometric correction will be referred to as median based directional scaling of pseudo-invariant features regression (MBDS).

Results