1. Introduction
Remote sensing satellites observe the land surface of the Earth at regular time intervals with the same observation geometry and obtain time-series of images, which record the occurrence and development patterns of land surface phenomena and thus have been widely applied for land change detection [
1], crop yield estimation [
2], urban sprawl analyses [
3], land-cover transition evaluations [
4], and forest succession analyses [
5], and achieving great success. Despite the great success of applications based on time-series of images, the physical signal recorded by a remote sensor (such as gray-scale value or reflectance) at different dates is inevitably contaminated by noise unrelated to the land surface, including different atmospheric conditions at the time of image acquisition and sensor distortion, which can cause variations in radiometric features between images and decrease the comparability between different images over the same study area [
6]. Contaminated signals can lead to sharp increases or decreases in the profiles of time-series indicators such as vegetation indices, which conceals the actual changes of the land surface and hinders information extraction. Thus, removing radiometric distortion is urgently needed to facilitate remote sensing applications.
Radiometric calibration attempts to eliminate the radiometric distortion caused by non-surface factors and correct radiometric differences between different images. Based on the transformation of gray-scale values to physical signals, calibrations can be classified into absolute radiometric calibrations and relative radiometric calibrations [
1]. An absolute radiometric calibration establishes the relationship between the measurement values from a remote sensor and the reflectance of the land surface to eliminate the radiometric distortion between images. This method needs to establish an “atmosphere–land surface–sensor” interaction model, involving certain environmental parameters (such as the atmosphere) at acquisition time [
7]. Additionally, some pseudo invariant calibration sites are required to calibrate the on-orbit sensors [
8]. However, not all archived historical data were recorded with environmental information, which restricts the practicability of this method [
5,
9]. The relative radiometric calibration uses a certain image as a reference and corrects another image based on the reference; thus, the processed image will have a similar radiometric condition as the reference for the same land surface, namely, radiometric normalization. Radiometric normalization directly establishes the mapping relationship of radiometric features between different images and can obtain an application effect comparable to that of absolute radiometric calibration [
5,
10]. Currently, there are primarily two types of methods for radiometric normalization; mapping and regression methods.
The mapping method directly establishes a gray-scale mapping equation between images and uses the mapping value to replace the gray-scale value of the input image. For instance, using a linear equation, the mean and standard deviation of a reference image can be assigned to the image to be processed, and the processed image therefore has the same average and standard deviation as the reference image, which eliminates the radiometric distortion. Another widely used mapping method is histogram specification, which assigns the histogram of the reference image to the input image so that the processed image has the same gray-scale distribution as the reference image [
11]. Because different bands of a multispectral image are usually correlated, the differences in the radiometric features can be eliminated by defining a high-dimensional rotation matrix to match the density function in multidimensional space [
12]. Radiometric normalization based on histogram specifications has been used for land surface change detection [
13], gap filling [
14], and image mosaicking [
11].
The regression method establishes a regression model to describe the radiometric distortion relationship between the images through pseudo-invariant features (PIFs), which represent pixels with radiometric features that do not change with time, such as those of buildings and bare land [
15]. As the radiometric difference between PIFs are mainly attributed to noise factors, we can establish a regression model to quantify the radiometric difference between different images and eliminates the radiometric difference according to the obtained regression model. As a result, the selection of high quality PIFs is the key for radiometric normalization. The PIF selection method mainly includes principal component analysis [
16], weighted principal component analysis [
17], multivariate alteration detection (MAD) [
18], improved Iteratively Re-Weighted MAD (IRMAD) [
9], and iterative slow feature analysis [
19]. The PIF extraction methods based on categories, for example, the temporally invariant cluster (TIC) method [
1], are more robust than that of pixel-level methods. For radiometric normalization models, linear regression models are simple and effective, and later improved linear regression methods such as the ordinary least-squares regression, reduced major axis regression [
20] and Theil-Sen regression [
21] have been successively developed and widely used to correct radiometric distortions. In addition, artificial intelligence methods, for example, genetic algorithms, can be used to not only optimize the regression parameters but also eliminate nonlinear distortions [
22].
According to the comparison of the two methods above, the mapping method changes variations in the radiometric features caused by atmospheric and sensor factors as well as land surface information. The physical meaning of the corrected image is not definitive; therefore, the mapping method is mainly suitable for generating visually seamless image mosaics [
11] or change detection [
13]. However, the regression method selects PIFs that are affected only by atmospheric and sensor factors and establishes a model to describe land surface radiometric variation [
15]. The regression method can maintain land surface information and eliminate image-related radiometric distortion, and these merits make it an ideal method for radiometric normalization of long time-series of images for land surface applications.
However, traditional PIF extraction and regression equation establishment methods are developed for two images and operate in an image-to-image manner [
1,
9,
15]. Although virtual reference image determination methods [
23] and overall adjustments [
16] have been developed, when these methods are extended to radiometric normalization for long time-series of images, they also operate in an image-to-image manner [
10]. For a radiometric normalization problem involving time-series of images with
n scenes, we should subdivide it into a group of
n − 1 image-to-image subproblems. Obviously, these extensions have two deficiencies. (1) The PIFs are extracted for each group of images separately, which may not be consistent between different groups of images; therefore, the correction models among different groups may not be comparable. Additionally, the invariance of different surface features has different time scales. The traditional selection of PIFs considers the changes of radiometric features at only two times and causes a potential error in PIF selection. (2) When an image is processed based on a reference image, although the radiometric residual in between can be minimized, but the residual for the correction between any two scenes of the images to be processed cannot be minimized. The time-series analysis views all images in the time-series as a continuously varying entity, thereby minimizing the radiometric residual for the overall correction between any two image scenes represents an optimal solution.
In recent years, with the development of computing technology and the free distribution of medium resolution images through the internet, applications utilizing long time-series of images can provide more temporal details with high accuracy and satisfy various demands; thus, they have become a popular research topic in remote sensing fields [
24,
25,
26]. However, according to the above analysis, traditional radiometric normalization is unsuitable for long time-series of images, and, therefore, a radiometric normalization method suitable for long time-series of images must be developed. To solve this problem and meet application demands, we will try to identify a better method for radiometric normalization of long time-series of remote sensing images in this paper. The main contribution of our method is twofold:
we developed a PIF selection method, which can consider all images in time-series for PIF selection and automatically suppress the negative effective of outliers, for example, clouds and cloud shadows; and
a novel optimization strategy is proposed to minimize the residual between the image to be processed and the images that have been processed previously, which can avoid the problem of reference image selection and obtain a smoother time-series profile.
The remaining part of the paper is structured as follows:
Section 2 introduces the research area and the experimental data;
Section 3 describes the principle of our method;
Section 4 introduces the implementation approach in detail;
Section 5 introduces the experimental results and their performance comparison;
Section 6 discusses the applicability of our method and its uncertainty; and
Section 7 presents the conclusions.
2. Materials
Hangzhou City in Zhejiang Province, China, is selected as our study area because of its rapid land use and cover changes and intensive urbanization in the last thirty years [
27]. The research area is indicated in
Figure 1b. To demonstrate our radiometric normalization method, we selected Thematic Mapper (TM) images from Landsat 5 with Path = 119 and Row = 38 in the Worldwide Reference System (WRS) for the experiment. From the website of the United States Geological Survey (USGS), we obtained a total of 438 scenes of image covering Hangzhou City from 1984 to 2010, the entire operational period of Landsat 5 (the distribution of the imaging dates is shown in
Figure 1c). Then, we conducted preprocessing that included image decompression and band synthesis, and a subregion of 2000 × 2000 (pixels) was clipped for our experiment. No further geometric processing was conducted because the geolocation accuracy of these images from the USGS was better than one pixel. The goal of this paper is to find a better method of radiometric normalization; therefore, we directly adopted the gray-scale values of different bands to process without converting them to physical signals such as reflectance. The gray-scale depths of the points in
Figure 1c represent the coverage proportion of noise (such as clouds) in each image (the detailed mask method will be described in the
Section 4.2). We can see that the proportion of clouds is approximately 50% on average, and the identification of cloud noise in the image is thus the basis for subsequent analyses of the time-series of remote sensing images.
After looking through all the images, the typical false-color composite images over Hangzhou City during four different time periods (1990, 1998, 2003, and 2008) are shown in
Figure 2. The built-up area of Hangzhou City in 1990 was mainly concentrated in the periphery of West Lake, and some villages and towns were distributed discretely in the research area. However, the land cover type in the suburb of Hangzhou City was mainly agricultural land. In approximately 2000, the large amount of farmland around Hangzhou City was transformed into cultivated land, and there was only a small amount of farmland remaining around Hangzhou City. In 2010, the extent of the built-up district expanded. Except for the mountains in the west of the research area, where it was difficult to use the land for construction, almost all the land was converted to urban buildings [
27]. The land use/cover changes of Hangzhou City indicates that the selection of PIFs is difficult for time-series of remote sensing images because the number of PIFs continues to decrease during long time periods. A better method of finding PIFs over long time-series of remote sensing images will benefit radiometric normalization and land use/cover analysis.
4. Methods
According to the analysis above, after segmenting the Arc into segments, our method can be implemented using the following steps (as
Figure 5). First, we classified the time-series’ observation values for every pixel into inliers (namely, normal observation values with good acquisition condition) and outliers (namely, abnormal observation values caused by cloud and cloud shadow noise). Second, PIFs were extracted based on the smaller and larger variation magnitudes of PIFs relative to those of vegetation and water, respectively. Next, we established the radiometric calibration equation for images acquired at different time points to eliminate the radiometric distortion between images. Therefore, we will introduce our method with four main steps as follows.
4.1. Arc Segmentation
For the time-series of observations
Pl for pixel
l (
) (as shown in
Figure 3a), we first segmented the Arc into segments using the following method.
Step 1.1: Sorting the time-series’ gray-scale values: The time-series’ gray-scale values of
Pl for pixel
were sorted in ascending order to generate a series of sequence values, which is expressed as
:
The sorting results are shown in
Figure 3b with the black solid line. This arc is denoted as
, and the points in
can be expressed as
The coordinates of points on can be expressed as . Here, represents the image sequence after the observation values are sorted and is a natural number in the range of [1, n]; and represents the gray-scale value corresponding to , which ranges from 0 to 255 for the TM image.
Step 1.2: Arc Segmentation with the inflection points: Connect point
, which is marked by
A in
Figure 3b and has a minimum gray-scale value, and
, which is marked by
B in
Figure 3b and has a maximum gray-scale value, to determine the straight line
. The equation of the straight line can be expressed as
According to Equation (6), we calculate the distance
for point
in the arc
to the straight line
:
The distances of all the points in the
form a vector:
Obviously, the point of inflection on
corresponds to the maximum value in the distance vector. Accordingly, the inflection point
C can be obtained by extracting the maximum distance in
Connect points
A and
C, and then points
B and
C, and then repeat the step mentioned above. Then, we can obtain two other inflection points named
D and
E. According to the inflection points obtained, we can split
into four segments as shown in
Figure 3c, which can be expressed as,
AD,
DC,
CE, and
EB.
4.2. Outliers and Inlier Identification
Step 2.1: Marking outliers and inliers: We took the segment comprising the minimum values (namely, the blue segment in
Figure 3c) as the abnormal observation value affected by a cloud shadow. The segment in the middle part (namely, the black segment in
Figure 3c) represents the normal observation value, and we denoted it as a clear line segment. The two segments of and with maximum gray values (namely, the red and pink segment in
Figure 3c) are marked as the abnormal observation values affected by cloud.
4.3. Extraction of Pseudo-Invariant Features
Step 3.1: Slope calculation for the clear line segment: First, we conducted a least squares fitting on the clear line segment and then used the slope to express the inclination degree of various line segments. Then, the research area can be expressed with one band of BS indicating the slope of the clear line segment of the pixel, which describes the gray-scale value variation magnitude.
Step 3.2: Extraction of PIFs: We set two thresholds TLow and THig and took the values of BS greater than TLow and smaller than THig as the PIFs. Here, the selection of the thresholds was not obtained by a universal criterion or an automated method. Instead, we selected the threshold manually to ensure the representativeness and integrity of the PIFs by statistically analyzing and observing the characteristics of BS, adjusting TLow and THig, and using the trial-and-error method.
More specifically, by analyzing the histogram of the Band BS (transformed to bins of integers), we defined two thresholds to exclude the land surfaces with low slopes (such as water) and land surfaces with high slopes (such as vegetation and changed land cover) first. We then segmented the remaining pixels into certain predefined groups. Then, through investigating the pixels in the various groups, we selected minimum and maximum gray-scale values of the group of the most invariant pixels (PIFs candidates) as the initial TLow and THig , and then adjusted the TLow and THig to exclude or include some pixels as PIFs.
Step 3.3: Exclusion of cloudy images with small amounts of clear PIFs: It is worth noting that images containing high proportions of clouds and cloud shadows can bring in noisy pixels. When using these cloudy images, the numbers of PIFs from different images substantially differ, and it is difficult to obtain sufficient PIFs to estimate a reliable radiometric normalization equation. In this study, we excluded these cloudy images by using two conditions, minimum clear PIFs number () and the coefficient of determination () between the image to be processed and the first image (the method for determination of the first image will be introduced in Step 4.2). If the was smaller than the threshold, the corresponding images were excluded from further radiometric normalization.
According to this method, we could obtain PIFs with relatively small changes in land surface features during the entire image acquisition period, and the obtained m PIFs are expressed as a set S. Obviously, S is a subset of L, namely, .
4.4. Radiometric Normalization Optimization
When the least squares method is adopted to normalize the image
based on the reference image
, the objective is to obtain normalized images denoted by
, which minimizes the correction error between
and
. This process is equivalent to solving the optimal solution under the constraint of the objective function
:
where
represents the radiometric normalization result of the image
,
represents the 2-norm of the vector, and
represents the selected PIFs. Clearly, this model cannot ensure an optimized solution with the smallest residual for all the corrected images.
Hence, in this paper, we propose a time-series normalization strategy under the constraint of the objective function
:
This objective optimization model not only minimizes the error between the normalized image and the reference image but also minimizes the error between the image to be processed and the image these were processed previously. As a result, a radiometric consistent time-series of images will be obtained. To find the optimal solution for the objective function , we designed the following algorithm.
Step 4.1: Sorting the standard deviations for PIFs: We calculated the standard deviation of the gray-scale value for the sample set
S for the
n scenes of time-series images and sorted them in descending order. The standard deviation of the gray-scale value after sorting can be expressed as
Obviously, the standard deviation satisfies
. The sequence numbers of the images are sorted according to standard deviation in descending order, which can be expressed as
Step 4.2: Correction of the first image: We first corrected image
to image
. For implementation, we obtained the correction parameters of
and
under the restricted condition
.
The restricted condition minimizes the error between corrected image of image and the corrected image of image ; this problem can be solved using the least squares method. Note that we denoted the image as for descriptive simplicity, the correction parameters of which are and .
Because the standard deviation for the gray-scale values of and satisfies , the correction coefficient of is expected to be greater than or equal to 1 in most situations; namely, the gray-scale value of the correction results is relatively stretched. Therefore, the compression of the gray-scale value was avoided in the correction process.
Step 4.3: Correction of other images in the time series: For image
in the time-series, we set the reference image
for
as follows:
The implementation method takes all the corrected images as reference images. The correction parameters of
and
for image
can be obtained under the restricted condition:
Obviously, this implementation method ensures the minimum error between the image to be processed and all the corrected images and is thus a greedy algorithm. Accordingly, the correction coefficient of image can be obtained using the least squares method.
We repeated this step for all the images in the time-series with the order determined by (12), and we obtained the correction coefficients for various images. The correction vector can be expressed as follows:
Step 4.4: Adjustment of correction parameters: As mentioned above, in the linear model, when the slope is
, the radiation resolution will be compressed. When
, the gray value will be negative [
16]. Assume that
is the minimum of
From this, we can calculate
If
< 1, we adjust
K according to the following method and derive the new slope vector
:
If
< 0, we adjust
B and obtain the new intercept vector
:
It is worth noting that our method performs sorting according to the standard deviation to ensure that the slope value is greater than 1, and a slope less than 1 is therefore less likely to appear. Additionally, the parameter adjustment will have a negative effect on calculating the error (the details will be analyzed in
Section 5.2.2). Therefore, we could omit this step decided by the user according to that
was near 1 or.
According to the correction coefficients and after adjustment, the near-infrared band of the image was corrected, and the resultant image of the radiometric normalization could obtain a smaller residual. According to the PIFs marked by this process and the corrected order, we then normalized the other bands of the image to obtain the resultant image from radiometric normalization.
6. Discussion
Our method for the radiometric normalization of long time-series Landsat images has the following innovations.
The identification of inliers (namely, normal observation values with good acquisition conditions) and outliers (namely, abnormal observation values caused by noise such as clouds and cloud shadows) in time-series of observation values for various pixels forms the basis for the radiometric normalization and further applications. As the number of images used for long time-series is very large, cloud detection is a time- and labor-intensive process. We captured the distinctive feature that variations in the gray-scale values of a clear pixel are concentrated in a narrow range, whereas the variations in the gray-scale values of outliers are much larger than those of clear pixels. We introduced the IBCD method for outlier identification, which alleviates the time and labor costs and obtains an acceptable result. Similarly, we selected clear-sky PIFs according to a novel measurement, namely, the slope of the clear line segment in the sorted time-series profile, which can exclude the negative effect of noisy pixels. This method can consider all the observations in the time-series images instead of just a small number of pixels, which increases the robustness of our method for undetected noise.
One objective of the time-series analysis is to maintain consistency among the observation values obtained at different times while reducing rapid up-and-down fluctuations between adjacent observations. The traditional correction method reduces the residual between the reference image and the image to be corrected, without constraining the residual between other images to be corrected. Our method, which is a typical greedy algorithm, takes all the corrected images as a reference for the images to be corrected, thus minimizing the residual for correction after adding each new image. The results also indicate that the residual obtained by our method is smaller than that of the contrasting method, and the observation values of the time-series are also smoother, which indicates the effectiveness of our method.
Our method can automatically determine the sequence of radiometric correction, which is very important for radiometric correction. If the gray-scale value distribution of the reference image is concentrated, while the gray-scale value distribution of the image to be processed is scattered, the gray-scale resolution of the images to be processed will be compressed compared with that of the reference image, and radiation information may be lost. However, in this paper, we developed a method that sorts the standard deviations of the gray-scale values of the PIFs in descending order. This approach guarantees that images with a wide gray-scale distribution are always corrected earlier and images with a narrow gray-scale distribution are corrected later. The image corrected later is stretched relative to that corrected earlier. The results of this strategy also indicate that the correction slope obtained by our method is generally greater than 1 (i.e., is relatively stretched). Finally, for the possible situation with a slope smaller than 1, we adjusted the parameters to ensure that the radiation resolution maintained a relatively large gray scale without compressing any pixels.
Additionally, some uncertainty may induce a decrease in the time-series normalization performance. Clouds and cloud shadows are important sources of noise in passive remote sensing and identifying abnormal observation values from time-series observation data is the basis of various applications. In this paper, according to the large fluctuation range of the gray-scale value caused by clouds, we developed a method based on the IBCD method to automatically identify the inliers and outliers to suppress the negative effect of contaminated pixels. Water bodies exhibit low reflectance in the near-infrared band and separating them from cloud shadows may be difficult; similarly, dense vegetation exhibits a high reflectance in the near-infrared band, and separating such vegetation from clouds may be difficult. For the PIFs (such as bare land, buildings, and roads), the gray-scale values in the near-infrared band are between those of low-reflectance water bodies and high-reflectance vegetation. Our method can easily differentiate abnormal observation values when the land cover type is PIFs. Therefore, the influence of these two deficiencies on the selection of PIFs is small.
From the image of the slope of the clear line segment, we can extract the PIFs using the threshold segmentation. At present, identifying PIFs requires manual trial-and-error to adjust the threshold of the slope, and they are evaluated and adjusted according to the obtained result. The unappropriated threshold setting may lead to failure in PIF selection. The time interval of the time-series of remote sensing images is long, and urban sprawl is dramatic in the research area. In a strict sense, no pixel is free from any changes. A PIF is a variable that cannot be precisely defined but instead refers to a type of pixel with a small variation magnitude.
Additionally, due to clouds, cloud shadows, and other noise factors, images have different numbers of clear PIFs, which may induce uncertainty in radiometric normalization. In the worst case, an image may not have sufficient PIFs to estimate a normalization equation, which will cause the image to not be processed. However, images that cannot be processed have large proportions of clouds or fewer clear pixels, and abandonment of these images will not result in a loss of much information.
Additionally, our method has certain drawbacks that may impede its application as follows: (1) this method needs relatively invariant pixels in the research area and thus has poor application potential in research areas that are completely covered by forests and farmlands; (2) in a region of rapid urbanization, the PIFs also experience variations in various radiometric characteristics, and many uncertainties are observed when there are changes in the entire time-series; (3) a long time-series is required to present the gray-scale curve of the images, which is unsuitable for short time-series; and (4) when one or more images have been accurately calibrated to reflectance, one intuitive method is to normalize other images to the corrected images, which provides the explicit physical meaning to the radiometric normalized result. However, our method cannot directly normalize the other images to the corrected image.
7. Conclusions
In this paper, we have proposed a radiometric normalization method for long time-series of remote sensing images, which exhibits favorable merits such as automatic outlier exclusion, PIF selection, and a novel strategy to minimize the RMSE between the image to be processed and the previously corrected images. In addition, we tested the method using long time-series of remote sensing data acquired by Landsat 5 TM for Hangzhou City. For experiments 1 and 2, the mean RMSEs of the images in the time series dropped from 22.97 and 17.73 (by contrasting method 1) to 17.39 and 13.87 (by our method), respectively, the standard deviations dropped from 8.51 and 6.12 to 5.93 and 4.51, respectively, and the means of the correlation coefficients between the time series gray-scale values increased from 0.508 and 0.562 (by contrasting method 1) to 0.781 and 0.793 (by our method), respectively, reflecting a significant performance gain by our method.
Additionally, the result indicates that our method can effectively eliminate differences in radiometric features between images and improve comparability between images. Moreover, the biophysical information from image time-series is well preserved, showing a smooth gray-scale value curve after radiometric normalization. The comparison between our method and the radiometric calibrated image demonstrates that our method provides a promising alternative method for radiometric normalization, especially when the parameters needed for absolute radiometric corrections are absent.