Radiometric Correction of Multitemporal Satellite Imagery 1

Problem statement: Repeated observation of a given area over time yiel ds potential for many forms of change detection analysis. These repe ated observations are confounded in terms of radiometric consistency due to changes in sensor ca libration over time, differences in illumination, observation angles and variation in atmospheric eff ects. Also major problem with satellite images is that regions below clouds are not covered by sensor . Cloud detection, removal and data prediction in cloudy region is essential for image interpretation . Approach: This study demonstrated applicability of empirical relative radiometric normalization met hods to a set of multitemporal cloudy images acquired by Resourcesat-1 LISS III sensor. Objectiv e of this study was to detect and remove cloud cover and normalize an image radiometrically. Cloud detection was achieved by using Average Brightness Threshold (ABT) algorithm. The detected cloud removed and replaced with data from another images of the same area. We proposed a new method in which cloudy pixels are replaced with predicted pixel values obtained by regression. Afte r cloud removal, the proposed normalization method was applied to reduce the radiometric influe nce caused by non surface factors. This process identified landscape elements whose reflectance val ues are nearly constant over time, i.e., the subset of non-changing pixels are identified using frequency based correlation technique. Further, we proposed another method of radiometric correction in frequen cy domain, Pseudo-Invariant Feature regression and this process removed landscape elements such as vegetation whose reflectance values are not constant over time. It takes advantage of vegetatio n being typically high frequency area, can be removed by low pass filter. Results: The quality of radiometric normalization is statis tically assessed by R 2 value and Root Mean Square Error (RMSE) between each pair of analogous band. Further we verified that difference in mean and sta ndard deviation is reduced after normalization of subject image with respect to reference image. Resu lts are compared with commonly used No Change regression method in spatial domain. Conclusion: Cloud removal depends on spatial registration between the two images (reference imag e and cloudy subject image). Visual inspection shows proposed cloudy pixel prediction method performs better than replacing cloudy pixels with another image of the same area. Statistical analysi s also shows that average RMSE of all bands is more in No Change method. Correlation in Fourier domain does not require water body in the scene, while No Change method does require.


INTRODUCTION
Satellite images are useful for monitoring changes in land use and land cover. But major problem with these images is that regions below clouds are not covered by sensor. The image distortion due to cloud cover is a classical problem of visible band of remote sensing imagery. Especially, for non-stationary satellite, it is commonly found in the earth resource observation application. Removing cloud cover from satellite imagery is very useful for assisting image interpretation. Hence cloud detection and removal is very vital in processing of satellite imagery. Further it is more difficult to quantify and interpret changes on multitemporal images under different illumination, atmospheric or sensor conditions without radiometric calibration. The relative approach to radiometric correction, known as relative radiometric normalization, is preferred because no in situ atmospheric data at the time of satellite overpasses are required (Yang and Lo, 2000). Relative radiometric normalization is a method of correction that applies one image as a reference and adjusts the radiometric properties of subject images to match the reference (Hall et al., 1991;Yuan and Elvidge, 1996). Thus, normalized images appear to have been acquired with the reference image sensor under atmospheric and illumination conditions equal to those in the reference scene (Hall et al., 1991). From the above descriptions, a key difference between the methods is that absolute normalization removes atmospheric errors that exist in both images, while relative normalization does not actually remove atmospheric errors. The subject image keeps the same atmospheric errors as the reference image. A variety of papers are available on image normalization methods (Schott et al., 1988;Eckhardt et al., 1990;Hall et al., 1991;Elvidge et al., 1995;Yuan and Elvidge, 1996;Heo and FitzHugh, 2000;Yang and Lo, 2000) reviewed existing image normalization methods.

Radiometric normalization methods:
Image normalization reported in literature can be classified into three general categories: statistical methods (i.e., standard deviation method); the histogram matching method; linear regression methods (i.e., PIF, DB and NC). These normalization methods are based on simple linear regression.
Normalization target selection: Several methods have been introduced by different authors showing how to select ideal targets for estimating the normalization transformation coefficients (Schott et al., 1988;Eckhardt et al., 1990;Hall et al., 1991;Elvidge et al., 1995). In general, for all normalization methods, the targets that meet the following criteria are selected as ideal targets for normalization (Eckhardt et al., 1990): • The targets should be approximately at the same elevation so that the thickness of the atmosphere over each target is approximately the same • The targets should contain only minimal amounts of vegetation because vegetation spectral reflectance is subject to change over time • The targets must be in relatively flat areas • When viewed on the image display screen, the patterns on the normalization targets should not change over time • A set of targets must have a wide range of grey values for the regression model to be reliable The mathematical equation for image normalization: Caselles and Garcia (1989) verified the linear relationship between reference image y and subject image x. The objective of linear spectral normalization is to rectify subject image x to reference Image y through a linear transformation. The common form for linear radiometric image normalization is: Where: x k = The DN of band k in image is x on date1 k y′ = The normalized DN of band k on date1 a k , b k = Are normalization constants for band k Image normalization can be divided into two steps: (1) selecting normalization targets and (2) determining normalization coefficients.
The methods reported in literature can be summarized as.
Haze Correction (HC) normalization: One approach to Relative Radiometric Normalization (RRN) is to remove the atmospheric Haze difference between two images. A Haze Correction (HC) is typically accomplished by subtracting the digital count associated with the darkest material present in scene. The concept is that the Haze value will be equal to the DN count observed for a ground surface having zero reflectance. Haze correction assumes that pixels having zero reflectance (darkest 0.1% of the image pixels) should have the same minimum DN values on both subject and reference images. The HC normalization coefficients are:  Simple Regression (SR) normalization: In this method, the subject image is regressed against the reference image in each band. Simple regression normalization (Jensen, 1983) uses least-squares to derive the normalization coefficients. The SR normalization coefficients are solved from the leastsquares regression equation: where the summation runs the whole scene. To solve this equation, one obtains the normalization coefficients as: Where: and: k k where, N is number of pixels in the scene.
Dark set-Bright set (DB) normalization: Hall et al. (1991) used the average of a set of dark and bright pixels (Dark-Bright set-simply called DB), which are extracted from the subject and reference image through Kauth-Thomas greenness-brightness transformation, to derive the normalization coefficients. It is assumed that an image always contains at least some pixels that have the same average surface reflectance among images acquired at different dates. The DB normalization coefficients are: x are the means of the bright set (b) dark set (d) and of band k in the reference and subject images respectively.
Pseudo-Invariant (PI) normalization: Schott et al. (1988) presented pseudoinvariant feature normalization, which analyzes the elements whose reflection distribution has statistical invariance, such as concrete, asphalt and rooftops. The pseudo-invariant features are extracted by analyzing the infrared to red ratio of the subject and reference images to identify pixels having low green vegetation cover and a NIR threshold to eliminate water pixels. Differences in the gray-level distributions of these invariant objects are supposed to be linear. Let the means and standard deviations of the selected pseudo-invariant sets for the two dates to be No-Change regression (NC) normalization: Elvidge et al. (1995) developed a radiometric normalization method, which locates the statistical centers for stable land and stable water data clusters using the near infrared date 1 versus date 2 scattergram to establish an initial regression line. The near-infrared data were used because at these wavelengths the spectral clusters for water and land are clearly separated and a distinct axis of "no change" can be observed.
A no-change set determined from the scattergram between near-infrared bands of the subject image and the reference image has used to compute normalization coefficients for all bands. Thus if the no-change subset NC is identified, the normalization coefficients are derived from following equations: where, (nc) k x and (nc) k y are the means. Sample variance and covariance for subset NC on two dates can be determined using Eq. 12 and 13: Here two satellite images are used, which are acquired at different time and different cloud distribution. The cloud is detected by Average Brightness Thresholding (ABT) algorithm. This is followed by relative radiometric normalization of the image. Here correlation procedure implemented in frequency domain is used for identification of nonchanging pixels. This approach is chosen because of its intrinsic simplicity than other methods. Also we suggested second method in frequency domain.
Cloud removal: The practical approach for identifying clouds in remotely sensed images of earth from satellite is studied. In tropical areas where cloud cover is common throughout the year, it is nearly impossible to obtain a cloud free image of even a moderately-sized land area. With these imperfect cloud contaminated image sources, we may still interpret images by switching to data from a second source (of the same geographic area) where image data of the main source is corrupted by cloud cover. The only requirement is that images should be acquired at different times; so that they have different cloud cover patterns.
Generally, bright areas are assumed to be clouds. Detection of cloud is achieved by Average Brightness Thresholding (ABT) algorithm (Leung and Jordan, 1995).
Average brightness thresholding algorithm is based on three observations: • Amplitude thresholding is simple to implement and provides quick processing • When correct threshold level (s) are chosen, amplitude thresholding is highly effective • Clouds are usually the brightest objects in satellite image The algorithm is quite simple to implement. First, the average brightness of the grayscale image is calculated. Next, threshold brightness is chosen based on the average brightness. Finally this threshold is applied to the image to divide it into cloudy and cloud free regions. The detected cloud is removed and replaced with data from another image of same area. The procedure is repeated for all bands.
ABT uses an average-cutoff function to determine appropriate brightness threshold level. This function has the characteristic that at low average brightness levels, cutoff is very much above the average, whereas at high average brightness levels, the cutoff is marginally above average brightness. We define cutoff function as:  Figure 1a shows original Resources at LISS III reference image and Fig. 1b shows cloudy subject image which is to be radiometrically corrected with respect to the reference image. Figure 2a shows detected cloud region in image acquired on 26 May 2006 using ABT algorithm. In Fig. 2b the detected cloud is removed and replaced with data from image acquired on 11 May 2004 of the same area. We suggested new method in which cloudy pixels are replaced with predicted pixel values obtained by regression. Result image is presented in Fig. 2c. Visual inspection shows that later performs better. Figure 3 shows histograms of image (band 4) (26 May 2006) before and after cloud removal, showing bright cloudy pixels is removed.
Radiometric correction: Radiometric correction of remotely sensed data normally involves the processing of digital images to improve the fidelity of the brightness value magnitudes (as opposed to geometric correction which involves improving the fidelity of relative spatial or absolute locational aspects of image brightness values). The main purpose for applying radiometric corrections is to reduce the influence of errors or inconsistencies in image brightness values that may limit one's ability to interpret or quantitatively process and analyze digital remotely sensed images. When a sensor records the solar energy on Earth's surface, the atmosphere affects both target radiance and irradiance. As sunlight pierces through atmosphere, it is both attenuated and scattered, reducing target illumination and making it diffuse. The atmosphere also acts as a scattering reflector, adding extra radiance directed back to sensor. When expressing these two atmospheric effects mathematically, total radiance recorded for the sensor can be related to object's reflectance at the surface and to irradiance: Where: L T = Total radiance is measured by sensor T R = Atmospheric transmittance L P = The radiance of atmosphere the target to sensor trajectory (and not of the object) from the scattering effect R = The object reflectance E T = Total irradiance reaching the earth There are number of important reasons to calibrate remote sensing data. The raw sensor DNs are simply numbers, without physical units. Each sensor has its own gains and offsets applied to the recorded signals to create DNs, they must be converted to at-sensor radiances. This step is sensor calibration. If we desire to compare surface features over time, or to laboratory or field reflectance data, corrections must be made for atmospheric conditions, solar angles and topography. This is atmospheric, solar and topographic correction. There are several levels of radiometric correction. The first converts the DNs to at sensor radiance and requires sensor calibration information.
The second is transformation of the at-sensor radiances to radiances at the earth's surface. This level is much more difficult to achieve, since it requires information about the view path atmospheric conditions at the time and locations of the image and sensor. That information can be in different forms, ranging from a simple categorization of the atmospheric conditions as one of several "standard atmospheres", to estimate of certain parameters such as path radiance from the image data itself, to coincident ground measurements. The final level of correction to surface reflectance is achieved by correction for topographic slope and aspect, atmospheric slope path length variation due to topographic relief, solar spectral irradiance, solar path atmospheric transmittance and down scattered skylight radiance.
Because of the complexity of full remote sensing image correction, there has been considerable interest in image based techniques that provide relative normalization in certain applications, such as multitemporal comparisons among images from the same sensor, or classification of multitemporal images before change detection.

Data processing:
The image data set used in this study is obtained from Resourcesat-1 LISS III sensor. The Resourcesat-1 (IRS-P6) satellite launched on October 17, 2003 is designed to provide multispectral, monoscopic and stereoscopic images of the earth's surface (Horne, 2003). It has LISS-III multi-spectral camera operating in four spectral bands, two in the visible and near infrared and one in SWIR region, as in the case of IRS-1C/1D. The new feature in LISS-III camera is the SWIR band (1.55-1.7 microns), which provides data with a spatial resolution of 23.5 m unlike IRS-1C/1D in which the spatial resolution is 70 m (courtesy: www.nrsa.gov.in).
Technique 1: Correlation procedure using Fourier transform: Relative radiometric correction is a method of correction that applies one image as a reference and adjusts the radiometric properties of subject image to match the reference (Chavez and Mackinnon, 1994;Yuan and Elvidge, 1996).
After removal of cloud from subject image and replacing detected cloud with data from other images of the same area, radiometric normalization is applied on cloud free subject image. Following assumptions are made.
• Camera parameters are not known • Imaging conditions are not known Hall et al. (1991) concept of radiometric rectification is based on stable reflectance control sets derived from the extremes of the image scattergram (Heo and FitzHugh, 2000). Other methods use stable reflectance targets or so-called pseudo-invariant features (Hall et al., 1991). Some approaches interpret changed pixels as outliers and adopt a strategy to eliminate, or attenuate them.
To identify no-change pixels, correlation in frequency domain is used. The reference image and subject image are divided into rectangular blocks of size 16×16 pixels. Blocks or extracts of image were used in a procedure that processes images from different dates band to band. A block of reference image is placed over block centered on the same coordinates in the other image. Then, normalized correlation between two corresponding blocks is calculated. This operation is repeated for all blocks. This procedure is repeated for all bands. After this, we applied a threshold criterion, in order to select nochange pixels, used to find normalization coefficients. A block is assumed to belong to the no-change set if it has normalized correlation in all bands greater than 0.9. Correlation can be calculated as: Where: f(m,n) = Block of 16×16 pixels of reference image w(m,n) = Block of 16×16 pixels of subject image m,n = Are spatial co-ordinates F(u,v) and W(u, v) = Fourier transform of f(m, n) and w(m, n)respectively F -1 = Inverse Fourier transform O = Correlation * = The complex conjugate Normalized Correlation is defined as: The DN values of no-change set are then used in a linear model as: To obtain the coefficients: Where: and: are the sample variance and covariance for the subset NC on two dates. |NC| is the number of pixels in the set NC.
Technique 2: Pseudo-invariant feature regression using Fourier transforms: Pseudo-invariant objects are the targets that have not experienced any significant change from date 1 to date 2 in terms of reflectivity (Jensen, 1983). Some of the Pseudo-invariant objects may be taken as roads, lakes and beaches, urban area, low vegetation area, industrial centers. This technique focuses on spatial information content of high resolution satellite images. It takes advantage of vegetation being typically high frequency area, can be removed by low pass filter.
First, Fourier transform is applied to nir band (band4) of multispectral images of two different dates. Fourier Coefficients of both images are then passed through low pass filter to remove vegetation area. Cut off is calculated from spatial frequency of vegetation region.
The spatial frequency for a given image is defined as follows: Consider an M × N image, where M = Number of rows and N = Number of columns. The row and column frequencies are given by: The total frequency is then: Here we selected cut off frequency of 26 to remove vegetation.
After filtering inverse Fourier transform is applied to obtain PIF set, which is then used for obtaining normalization coefficients using linear least squares regression. Normalization coefficients obtained are used for normalization.

Algorithm:
1: Decompose nir band of reference image using Fourier transform.

RESULTS
The quality of radiometric normalization can be statistically assessed with the use of R 2 value and Root Mean Square Error (RMSE) between each pair of analogous bands.
The RMSE is defined as: The R 2 value is defined as: where, n is total number of pixels in the scene.    Closer the R 2 value to one, better the radiometric process. RMSEs and R 2 values are listed in Table 1. Results show that RMSEs of all bands of normalized image are less than uncorrected image and improvement in R 2 value. This implies that normalized image is more similar to the reference image. Table 2 gives statistical results for correction of another set of multitemporal images which are not contaminated by cloud, in which Subject image (8 Jan. 2007) is corrected by reference image (2 March 2006).
Further we verified that difference in mean and standard deviation is reduced after normalization of subject image with respect to reference image, listed in Table 2 and 3 for two sets of multitemporal images. Figure 4a shows reference image acquired on 11 May 2004. Figure 4b is subject image (26 May 2006) after cloud removal and cloudy pixels are replaced by data from another image of same area. Here we used three different techniques for normalization of image.
Raw means between subject and reference image before normalization.
Result of applying radiometric normalization using correlation in Fourier domain (technique 1) is presented in Fig. 4c. Figure 4d is Corrected image using commonly used NC method in spatial domain. Figure 4e corrected image using proposed (technique 2) PIF in Frequency domain method.
In PIF method we used low pass filter to remove vegetation, as it is high frequency region. Figure 5 a shows Low pass filter with cut off 26. Figure 5b shows Filtered reference image (2 March Fig. 6a and b. Figure 6c shows result of applying PIF using Fourier transform radiometric correction on subject image. Figure 6d is normalized image (8 Jan. 2007) using Correlation in Frequency domain method. Figure 6e shows corrected image using NC Method in spatial domain.

DISCUSSION
The detected cloud if removed and replaced with data from image acquired on 11 May 2004 of the same area, patches are seen in the result image. We suggested new method in which cloudy pixels are replaced with predicted pixel values obtained by regression. Visual inspection shows that later performs better.
Visually much difference is not seen in the results of proposed two radiometric normalization using Fourier transform methods. However, result of spatial domain method is different. Statistical analysis also shows that average RMSE of all bands is more in NC spatial domain method.
We tested the proposed methods on another set of multitemporal images. These images are not corrupted by cloud cover so we applied radiometric normalization techniques only. Average RMSE after normalization is less in correlation in frequency domain than other two techniques. Improvement in R 2 value is also better than other methods. Proposed PIF method work very well visually as well as statistically. But this requires water bodies in the scene. Visual and statistical analysis shows that proposed methods work very well.

CONCLUSION
In this study, two techniques for relative radiometric correction of cloudy multitemporal satellite images are proposed. The proposed techniques successfully normalize cloudy subject image. However, Cloud removal depends on spatial registration between the two images (reference image and cloudy subject image). We suggested two radiometric normalization methods in frequency domain. Visually much difference is not seen in the results of proposed normalization methods. However, result of spatial domain method is different. Statistical analysis also shows that average RMSE of all bands is more in NC method. Correlation in Fourier domain does not require water body in the scene, while NC method does require. Subject image looks close to reference image after normalization.
We tested the proposed radiometric correction techniques for normalization of another set of images which do not contain cloud, visual inspection and statistical analysis shows that proposed methods work very well.
Proposed methods of radiometric normalization have following advantages: • We suggested PIF method based on spatial frequency. We observed that proposed PIF method work very well visually as well as statistically. But these require water bodies in the scene Computing correlation in the frequency domain using the fast Fourier transform is more efficient than in spatial domain • Distribute normalization error among different landcover types • Eliminate the necessity of identifying bright and dark radiometric control pixels • Accelerate the speed of the normalization procedure.
Correlation technique does not require the presence of both land and water areas in the satellite images • The results are not sensitive to outliers in the data