A Nonlinear Radiometric Normalization Model for Satellite Imgaes Time Series Based on Artiﬁcial Neural Networks and Greedy Algroithm

: Satellite Image Time Series (SITS) is a data set that includes satellite images across several years with a high acquisition rate. Radiometric normalization is a fundamental and important preprocessing method for remote sensing applications using SITS due to the radiometric distortion caused by noise between images. Normalizing the subject image based on the reference image is a general strategy when using traditional radiometric normalization methods to normalize multitemporal imagery (usually two or three scenes in different time phases). However, these methods are unsuitable for calibrating SITS because they cannot minimize the radiometric distortion between any pair of images in SITS. The existing relative radiometric normalization methods for SITS are based on linear assumptions, which cannot effectively reduce nonlinear radiometric distortion caused by continuously changing noise in SITS. To overcome this problem and obtain a more accurate SITS, we propose a nonlinear radiometric normalization model (NMAG) for SITS based on Artiﬁcial Neural Networks (ANN) and Greedy Algorithm (GA). In this method, GA is used to determine the correction order of SITS and calculate the error between the image to be corrected and normalized images, which avoids the selection of a single reference image. ANN is used to obtain the optimal solution of error function, which minimizes the radiometric distortion between different images in SITS. The SITS composed of 21 Landsat-8 images in Tianjin, China, from October 2017 to January 2019 was selected to test the method. We compared NMAG with other two contrast methods (Contrast Method 1 (CM1) and Contrast Method 2 (CM2)), and found that the average root mean square error ( µ RMSE ) of NMAG ( 497.22 ) is signiﬁcantly smaller than those of CM1 (641.39) and CM2 (543.47), and the accuracy of normalized SITS obtained using NMAG increases by 22.4% and 8.5% compared with CM1 and CM2, respectively. These experimental results conﬁrm the effectiveness of NMAG in reducing radiometric distortion caused by continuously changing noise between images in SITS.


Introduction
Satellite Image Time Series (SITS) is a data set that includes satellite images across several years with a high acquisition rate. SITS can provide abundant information to describe temporal changes in the generation and development of ground features in an area [1]. Thus, it has been used as an important data source in many fields such as environmental monitoring, land cover change monitoring, crop growth monitoring, and so on [2][3][4][5]. However, the temporal information extracted from SITS is inevitably disturbed by noise unrelated to ground features, such as atmospheric absorption and scattering, sensortarget illumination geometry, sensor calibration, etc., which lead to inaccurate results in remote sensing applications [6]. Thus, radiometric normalization is required prior to any applications using SITS.
Radiometric calibration can be classified into absolute radiometric calibration and relative radiometric calibration (also named relative normalization). Because it is usually difficult to obtain atmospheric properties, performing relative radiometric calibration based on the inherent radiometric information of images provides an alternative when absolute surface radiances are not needed [7]. In early research methods, global image statistics were used to directly establish the gray value mapping relationship between the subject image and the reference image, such as Histogram Matching (HM), Mean-Standard deviation (MS), and so on [8,9]. However, these methods reduce the radiometric difference resulting from the variation in ground features while reducing the radiometric distortion caused by noise. Therefore, these methods are only suitable for image mosaicking due to the ambiguous result of the physical meaning of normalized images [10]. To effectively address this issue, methods based on the regression model have been developed, which is established using a set of invariant pixels [11] (named Pseudo-Invariant Features (PIFs)) from the subject and reference images.
Linear regression models are simple and effective, so they have been widely used to minimize radiometric distortion. Jenson et al. proposed a method based on a simple linear regression (SR) model, which obtains the linear relation between two scenes using the least square method [12]. Du et al. proposed an objective normalization procedure to find the most likely linear relation between the subject image and reference image through a principal component analysis (PCA) [13]. Olthof et al. developed a Theil-Sen regression model to avoid the error propagation caused by outliers [14]. Ghanbari et al. obtained the linear function between subject and reference images through the Error Ellipse (EE) process with a determined confidence level, which can removes the outliers from the linear model fitting [15]. However, relative normalization methods based on linear assumptions are inappropriate for dealing with complex nonlinear radiometric distortions between the subject image and the reference image. Radiometric normalization models based on artificial intelligence method, such as genetic algorithms [16], artificial neural network [17], random forest [18], kernel canonical correlation analysis [19], etc., have been used as an effective tool to deal with this problem.
The relative radiometric methods mentioned above are unsuitable for calibrating SITS, because they are based on the reference image to implement the radiometric normalization of the subject image, which cannot minimize the radiometric distortion between any pair of images in SITS. Wu et al. developed a new radiometric normalization procedure for SITS to effectively solve this issues, obtaining more objective and accurate correction results than previous methods [20]. However, continuously changing noise in SITS usually results in nonlinear radiometric distortion between images in reality, and the method proposed by Wu et al. cannot effectively reduce such nonlinear radiometric distortion.
An optimum method for minimizing nonlinear radiometric distortion between images in SITS is needed to facilitate remote sensing application. Thus, in this study we constructed a nonlinear radiometric normalization model for SITS, in which an artificial neural network (ANN) and greedy algorithm (GA) are combined for the nonlinear radiometric normalization of SITS. Here, GA is used for determining the correction order of SITS and calculating the error function of the image to be corrected, and ANN is used for obtaining the optimal solution of the error function. Hereafter, this nonlinear radiometric normalization procedure is called NMAG model.

Description of Methodology
The existing radiometric normalization approaches are based on the reference image X r to correct the subject image X c . The key of these methods is finding a mapping function f () to minimize the radiometric distortion between the subject image f (X c ) and reference image X r . The objective function Q r can be expressed as: where S represents Pseudo-Invariant Features (PIFs), s ∈ S represents the Pseudo-Invariant Feature (PIF) in the PIFs, f (X s c ) represents the Digital Number (DN) value of PIF in the normalized subject image f (X c ) , and X s r represents the DN value of the PIF in the reference image. When Equation (1) is used for the correction of SITS X, an image X r ∈ X is selected as the reference image, and then other images X i in X are normalized to the reference image X r one by one. Thus, the objective function Q r here can be written as: where n represents the number of images in X, and f (X s i ) represents the DN value of the PIF in the image to be corrected f (X i ). It can be easily found from Equation (2) that Q r cannot ensure that the radiometric distortion between each pair of images in normalized SITS is minimum. Therefore, Wu et al. proposed an improved objective function Q G for the relative radiometric normalization of X, and Q G is expressed as [20]: Our NMAG method is based on this objective function to achieve nonlinear radiometric normalization for SITS data by using artificial neural networks and greedy algorithms. The NMAG approach can be implemented using the following steps as shown in Figure 1, and the code of entire experiments are available from the public repositories (named NormSITS) (https://github.com/zhaohyin/NormSITS (accessed on 26 February 2021)) in GitHub.
Step 1. Create a greedy algorithm [21]. A greedy algorithm is an intuitive, well-tested algorithm used in optimization problems [22,23]. The algorithm makes the optimal choice at each step as it attempts to find the overall optimal method to solve the entire problem. To ensure that Q G can obtain the optimal solution, the greedy algorithm should be created to adopt the most greedy solution when implementing the rediometric normalization of each image in SITS. SITS X is divided into two groups, one as reference image dataset R and the other as image dataset U to be corrected: When NMAG has not yet started normalization, R and U can be expressed as: A clear and cloudless image X r is selected from SITS as the reference image, and then X r is added to the reference image dataset R. Update U according to Equation (6), then an image with the smallest rediometric distortion from R is selected from U as the next image to be corrected U c . The selection of U c can be given by Equation (7).
where k represents the number of normalized images in R, m represents the number of uncorrected images in U, and s ∈ S represents the PIF in PIFs. R s i represents the DN value of PIF in the i-th image of R, and U s x represents the DN value of the PIF in the x-th image in U.
Next, implement the radiometric normalization of U c , and the local optimal solution Q c can be calculated by: where f (U c ) represents the radiometric normalization result for U c , and f (U s c ) represents the DN value of the PIF in f (U c ) .
Step 2. Generate an ANN regression model. The linearity assumption of radiometric distortion in SITS is imprecise, and the nonlinear regression model may achieve more accurate normalization results of SITS. As a representative method of machine learning, ANN regression model has been widely used to solve nonlinear regression problems in remote sensing applications and has achieved significant results [24,25]. The schematic diagram of ANN [26] is shown in Figure 2. In this study, PIFs were randomly divided into a training sample set S l and a test sample set S t , accounting for 70% and 30% of total PIFs respectively. S l was used for the training of ANN regression model, and S t was used for the accuracy assessment of ANN regression model after training. As shown in Figure 2, ANN includes an input layer, an output layer, and several hidden layers. The input layer is the DN value of the training samples in the image U c to be corrected, which is expressed as U s c , s ∈ S l . The output layer is the radiometric normalization results f (U s c ), s ∈ S l . The experimental details of training and testing ANN, including the number of training samples and testing samples, the parameter of ANN, the ANN error convergence graphs for training etc., are recorded in detail in the Technical Document. This document is also available from the the public repositories (named NormSITS) (https://github.com/zhaohyin/NormSITS/blob/master/ README.md (accessed on 26 February 2021)) in GitHub.
The ANN regression model, which consisted of an input-output pair, can be adjusted by the connection weights between the nodes by learning to memorize every one of the network learning and training samples. The magnitude of weights represents the importance of the link between the neurons to f (U s c ). The value of weights is initialized randomly, thus, the initial output is also random.The difference (i.e., error) between the generated output and a training set output is calculated and is fed back to the network, where it is used for connection-weight readjustment by Gradient Descent Algorithm (GDA) to minimize the error to within a predefined tolerance. In this way, the iteration calculation can reduce the error to a predetermined allowable range. The error estimation function can be expressed as: where k represents the number of images in R, and R s i represents the DN value of PIFs in the i-th image in R.
Step 3. Loop execution. Add the corrected image f (U c ) to reference image dataset R and remove U c from uncorrected image dataset U. The update of R and U can be expressed as: where f (U c ) is the radiometric normalization results of image U c . If images still exist in the updated U, we execute the relative radiometric correction of the next image in the same way.

Study Area and Satellite Data
Here, to illustrate the application of NMAG modle, SITS with 21 high-quality Landsat-8 satellite images covering an area of Tianjin, China, were used in this study. The acquisition time of these images ranges from October 2017 to January 2019 and the average cloud cover of images is 3.58%. The cost time of performing relative radiometric correction on entire images in SITS is too long; thus, only 1000*1000 pixels were selected as the study area. The false color composite image of the study area is shown in Figure 3. No further geometric correction is needed during image preprocessing because the positioning accuracy of these images obtained from the United States Geological Survey (USGS) (https://glovis.usgs.gov (accessed on 26 February 2021)) is less than one pixel. The false-color composite images of each temporal datum in the SITS is given in Figure 4, and we can find from this figure that this area is rich in ground features, of which artificial buildings are mainly concentrated in the urban area of Tianjin, bare land is mainly distributed in the suburbs of Tianjin, and road pixels are evenly distributed in the study area. The PIFs were mainly selected from pixels corresponding to these three types of ground features. Vegetation with intra-year seasonal changes are mainly distributed in suburbs of Tianjin, which can be selected to assess the accuracy of the radiometric normalization results of SITS.  Table 1 shows that the range of Pearson Correlation Coefficient between images both acquired from October 2017 to December 2017 ranges from 0.71~1.00, and the average value is 0.88, indicating that a strong positive linear correlation exists between images with close acquisition times. Table 2 shows that the range of Pearson Correlation Coefficient between images acquired from October 2017 to December 2017 and images acquired from October 2018 to December 2018 is 0.61~0.75, and the average value is 0.67, indicating that the linear correlation between images which acquisition time is far apart is weak. Due to both linear and nonlinear radiation distortions existing between images in SITS in this area, the use of SITS in this area is more beneficial for testing the effectiveness of NMAG model.

The Preparation of PIFs
PIFs refer to pixels with constant radiometric value in SITS, which were mainly selected from artificial buildings, roads or bare ground pixels. Using high-quality Landsat-8 images of the study area can avoided the interference of clouds and shadows while selecting PIFs, which help us to select high-quality PIFs. The SITS used here can be expressed as: The least square method was used to estimate the best-fitting parameters (k p and b p ) of the time-series' DN value of pixel p sorted in ascending order DN X p . Figure 5a shows the distribution of k p values in the study area. The slope of fitted line of PIFs (k s ) is higher than that of water bodies and lower than that of vegetation [20]. Thus, the time-series' DN values of pixel p sorted in ascending order can be given by Equation (12) and the range of k s is expressed as: where P represents all pixels in images in this area, S represents selected PIFs, and DN X p represents the time-series' DN value of pixel p sorted in ascending order. As seen in Figure 5a, water bodies are displayed in light blue, artificial buildings and bare lands are shown in gray, and vegetation is expressed in orange or dark brown. The results indicated that the k p value of the water bodies is the lowest (k p < 250), the k p value of the vegetation is the highest (k p > 500), and the k p values of the artifical buildings and bare lands are between those of the water bodies and the vegetation (250 < k p < 500). PIFs represent pixels with unchanged radiometric values over time, they were mainly selected from artificial buildings and bare lands. Thus, we set the thresholds of k s used in this paper to extract PIFs were 275 for k Low and 300 for k High . The extraction results of PIFs are shown in Figure 5b, and white pixels represent PIFs, with a total of 12,515 pixels. Of these, 8760 pixels are training samples and 3755 pixels are testing samples.

Experimental Results
As an illustrative example, we applied NMAG model to the relative radiometric normalization of Landast-8 SITS. Figure 6 shows the normalized results for SITS.    These results indicate that the radiometric distortion between different images caused by noise in SITS can be effectively suppressed by NMAG method. Figure 7d shows the time-series' curve of DN values (DN-TSC) from the near-infrared Band (Band 5) of road pixels before and after radiometric normalization. The untreated DN-TSC (marked by a black line) fluctuates greatly, which indicates that the noise can conceal the actual variation in the spectral characteristics from ground features. After radiometric normalization, the DN-TSC fluctuates slightly and nearly lies to a straight line. This indicates that NMAG effectively suppresses the radiometric distortion between different images in SITS and results in DN values between images in SITS become more comparable. Figure 8c shows the DN-TSC from the near-infrared Band (Band 5) of pixels in cropland images before (black) and after(red) radiometric normalization. In general, a single peak of DN-TSC exists during the entire corn growth period. According to information obtained from the website of Tianjin Agriculture and Rural Committee (http://nync.tj.gov.cn (accessed on 26 February 2021)) and crop images during the corn growth period (see Figure 8b), we know that corn in this area are generally sown in early June and harvested at the end of October. This means that DN values of cropland pixel is similar to that of bare land pixel from November of the previous year to June of the current year, which should theoretically maintain a fixed value. After radiometric normalization, the DN-TSC (marked in red line) fluctuated slightly from November 2017 to June 2018, demonstrating that NMAG method effectively minimizes the disturbance of noise. A sigle peak (8 September 2018) existed in the normalized DN-TSC during the entire growth period, which indicated that NMAG can enhance the real time-series' characteristics for cropland pixels.

Comparison with Other Methods
We compared NMAG module with other two contrast methods using the same data. Contrast Method 1 (CM1) proposed by Sadeghi et al. is a nonlinear radiometric normalization method based on the reference image to normalize the subject image [17]. When CM1 is used for the relative normalization of SITS, one image in SITS should be selected as a reference image first and then the other images are normalized to the reference image one by one. The Contrast Method 2 (CM2) proposed by Wu et al. is a radiometric normalization method for SITS based on linearity assumption [20]. This method takes all the corrected images as reference images, and then normalize the subject image by using the Least Squares Method (LSM). To evaluate the accuracy of radiometric normalization results obtained using different methods for SITS, we computed the root mean square error RMSE( f (X i ), f (X j )) to measure the radiometric distortion between the normalized image f (X i ) and the normalized image f (X j ), which is can be calculated from Equation (14). The smaller the RMSE value, the more accurate the radiometric normalization results of images.
where X i represents i-th image in SITS, X j represents j-th image in SITS. S t represents the test samples of PIFs, N represents the number of PIF in S t , s represents the PIF in S t , and f (X s i ) represents the DN value of PIF in the radiometric normalization result for image X i . Figure 9a-c shows error matrices of RMSE for the normalized results obtained using NMAG model and two contrasting methods with the same SITS data. In this figure, as the color transitions from blue to red, the RMSE increases. Compared with the error matrices of RMSE obtained using NMAG (Figure 9c), more red pixels are observed in the error matrices of RMSE obtained using CM1 (Figure 9a) and CM2 (Figure 9b). This qualitative analysis result shows that the error matrices of RMSE obtained using NMAG is generally distributed at low values. Figure 9. The error matrices and the frequency distribution curve of the root mean square error (RMSE) for the normalized results obtained using NMAG and two contrasting methods (refered as Contrast Method 1 (CM1) and Contrast Method 2 (CM2)). As the color transitions from blue to red, the radiometric distortion increases. The error matrices of RMSE obtained using (a) CM1, (b) CM2, and (c) NMAG. (d) The frequency distribution curve of RMSE obtained using NMAG and two constrasting methods. Figure 9d shows the RMSE frequency distribution curve of NMAG (marked by a red line) and other two methods (marked by black line and blue line, respectively). The RMSE between various images of CM1, CM2 and NMAG is concentrated around 700, 550, and 400, respectively. Thus, the error of our method between images in SITS is significantly lower than that of CM1 and CM2. The comparison result demonstrates that NMAG method can reduce the radiometric distortion caused by noise between images in SITS more effectively than two contrasting methods, and obtains more accurate normalized results.
We further compared NMAG method with CM1 and CM2 by calculating the average RMSE( f (X i ), f (X j )). This parameter can be calculated by: where n represents the number of images in SITS. f (X i ) and f (X j ) are respectively i-th and j-th image in normalized SITS. Table 3 shows that the µ RMSE of NMAG (497.22) is significantly smaller than those of CM1 (641.39) and CM2 (543.47), and the accuracy of SITS obtained by using NMAG has increased by 22.4% and 8.5% toward CM1 and CM2, respectively. These results indicate that NMAG model can obtain more accurate SITS.

Application of NMAG to Vegetation Index
Vegetaion index (VI) obtained from multi-band image data can better reflect the green vegetation status than DN value from single band data. Therefore, the time-series' vegetation index is often used in the field of land cover change monitoring, environmental monitoring and so on [27,28]. As the normalized difference vegetation index (NDVI) is the most frequently used VI in remote sensing applications [29], we further examined the application of NMAG to NDVI. The NDVI can be expressed as: where ρ N IR is the top-of-atmosphere (TOA) reflectance from the near-infrared Band and ρ R is the top-of-atmosphere (TOA) reflectance from the red Band. With the help of metadata file (_MTL.txt), we can easily transform DN value into TOA reflectance. We used NMAG model to normalize the near-infrared band and the red band separately, and normalized NDVI could be calculated by using normalized result of TOA reflectance of near-infrared Band and red Band. Figure 10 shows the comparison of time-series' curve of NDVI values(NDVI-TSC) calculated from pixels in corpland images before (black) and after (red) the radiometric normalization. These two NDVI-TSC exist significant differences in details though they have the same trends. The fluctuation amplitude of NDVI-TSC obtained using NMAG is lower than that of untreated curve from November 2017 to May 2018. As mentioned in Section 4.1, the DN values of cropland pixel are similar to those of bare land pixels from November of the previous year to June of the current year, which should maintain a fixed value, theoretically. This means the NDVI value of pixels should be also theoretically maintain a fixed value from November 2017 to May 2018. The reduction in the fluctuation amplitude of NDVI-TSC indicates that NMAG can effectively decrease noises contained in different bands of SITS.
The times corresponding to the turning point of NDVI-TSC before and after the radiometric normalization are 3 and 19 May 2018. This means that the corn sowing time corresponding to the untreated NDVI-TSC is between 3 and 19 May 2018, and the corn sowing time corresponding to the normalized NDVI-TSC using NMAG is between 19 May and 4 June 2018. In reality, the corn sowing occurred in early June, 2018 according to the information obtained from the website of Tianjin Agriculture and Rural Committee. The normalized result by using NMAG fully matchs the actual situation, showing that the timer-series' NDVI values obtained by using NMAG to normalize SITS data more accurately reflect the vegetation coverage, which can improve the accuracy of remote sensing applications.

Discussion
The objective of this study was to supply a model to reduce the nonlinear radiometric distortion caused by continuously changing noise between images in SITS, thereby obtaining a more accurate SITS. So we proposed a nonlinear radiometric normalization method (named NMAG) for SITS based on GA and ANN. A SITS composed of 21 Landsat-8 images in Tianjin, China from October 2017 to January 2019 was used to test NMAG. The experimental results confirmed the effectiveness of NMAG in the reduction of radiometric distortion caused by continuously changing noise between images in SITS. In addition, the performance of NMAG surpasses that of the existing methods (CM1 and CM2) upon comparison.
In theory, NMAG is also suitable for SITS acquired by other satellites, provided that enough relatively invariant pixels can be selected as training samples in the study area. If the number of training samples is insufficient, the model training of NMAG may be insufficient, and the performance of NMAG would be negatively impacted. We recommend that the total number of selected PIFs exceeds 10,000, and the proportion of training samples should not be less than 70%. In the same research area, the higher the spatial resolution of satellite images is, the more convenient the selection of PIFs would be. This means that NMAG may be more convenient for the relative radiometric normalization of satellite image data with higher spatial resolution than Landsat-8. Notably, with increasing number number of PIFs, NMAG model training will take longer, and the performance of NMAG will be increase. Please choose an appropriate number of PIFs according to the hardware performance and the actual needs of the research.
However, NMAG has a certain limitation that may hinder its application. As mentioned above, NMAG needs enough PIFs in the study area to train the model. Thus, NMAG has weak application potential in study areas with high vegetation coverage or water coverage. Additional studies are needed to find effective models to solve the problem of the scarcity of PIFs that may occur in practical applications.

Conclusions
In this paper, we propose a nonlinear radiometric normalization method NMAG for SITS based on GA and ANN. This method can effectively suppress the noise contained in SITS, resulting in the gray values of images acquired at different times being more comparable. In the method, GA was used to determine the correction order of SITS and calculate the error between the image to be corrected and corrected images, which avoided the selection of a single reference image. ANN was used to obtain the optimal solution of error function, which minimized the radiometric distortion between images in SITS.
SITS composed of 21 Landsat-8 images in Tianjin, China from October 2017 to January 2019 were used to test our method. The resultant images have similar color and color contrast to each other, and the normalized DN-TSC of the near-infrared Band (B5) obtain by using NMAG can better reflect the actual change of different features than the original DN-TSC. In addition, we compared NMAG with other two existing methods (CM1 and CM2) using the same data. The result shows that the µ RMSE of NMAG (497.22) is significantly smaller than those of CM1 (641.39) and CM2 (543.47), and the accuracy of SITS obtained by using NMAG has increased by 22.4% and 8.5% toward CM1 and CM2, respectively. This indicates that NMAG can obtain more accurate SITS than other two contrasting methods.
Because the NDVI obtained from multi-band image data can more accurately reflect the green vegetation status than DN value from single band data, we further analyzed the application of NMAG to Vegetation Index. The NDVI-TSC obtained using NMAG to normalize SITS data is fully matched the actual situation, indicating that NMAG can effectively reduce the rediometric distortion caused by noise in SITS so that we can obtain more precise time-series' NDVI values.