A Relative Radiometric Normalization Method for Enhancing Radiometric Consistency of Landsat Time-Series Imageries

Radiometric consistency of multitemporal satellite observations is affected by sensor stability and scene related issues. Relative radiometric normalization (RRN) is a widely used method to reduce these radiometric differences, its performance depends on the accurate identification of representative pseudoinvariant features (PIFs). However, existing RRN methods are mainly developed for bitemporal images and are limited to time-series imageries due to the complexity of identifying effective PIFs. In this study, we proposed a novel RRN method to enhance the radiometric consistency of Landsat time-series imageries. This method includes the following: first, a trend-based PIFs identification considering land cover changes and phenological trends from the entire time series; second, a PIFs optimization involving an automatic reference selection and a PIFs refining for each reference–target image pair; and third, a combined RRN modeling using the M-estimator sample consensus algorithm and robust linear regression. The Landsat surface reflectance products were used to validate the proposed method. The experimental results showed that the trend-based PIFs identification provided the consistent PIFs for all reference–target image pairs; aided by an automatic reference allocation, PIFs optimization filtered the proper PIFs with high spectral and spatial similarity for each image pair in monthly image stack; the proposed RRN method achieved good performance in model precision and radiance consistency improvement; the proposed RRN method outperformed seven commonly used RRN methods on majority images in image stack of December. The normalized images can help generate more comparable time-series analysis results by reducing the uncertainties from radiometric calibration, atmospheric correction, and sensor differences.

Radiometric normalization can be categorized as absolute (ARN) and relative (RRN) approaches. ARN uses physical parameters to convert digital values to surface reflectance in conjunction with sensor calibration and atmospheric correction [11]. Compared to the difficulties in collecting synchronous satellite data and building the radiative transfer model for ARN [12], RRN is easier to be implemented, which performs an image-toimage transformation by adjusting the radiometric properties of a target image to match those of a reference image [11], [13]. Because RRN cannot remove the differences caused by atmospheric conditions [14], combining RRN with ARN is used to generate normalized results with consistent physical meaning [10], [15].
RRN can be broadly divided into two categories: global statistics-based methods and radiometric control set samplebased methods [16], [17], [18], also called dense RRN (DRRN) and sparse RRN (SRRN) [19]. DRRN uses all pixels of reference and target images to determine the adjustment of radiometric properties, such as histogram matching [20], haze correction (HC) [21], minimum-maximum (MM), mean-standard deviation (MS) [22], and simple regression (SR) [23]. DRRN methods are of low complexity but sensitive to land cover changes and clouds of images [2], [24]. Comparatively, SRRN methods identify invariant pixels from a reference-target image pair and use them to establish a mathematical relationship for adjusting the target image, which can effectively reduce the negative effects of changed pixels on normalization. Some conventional methods include the dark and bright sets method [11], pseudoinvariant features (PIFs) regression [25], and the no-change sets method [13]. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ The basic hypothesis of SRRN methods is that the invariant pixels at time t1 are linearly related to the corresponding pixels at t2 [26], [27], although a nonlinear relationship is also applied for multisensor images in cases [14], [24], [28]. Therefore, the key to SRRN is to identify the representative and accurate invariant pixels to illumination and land cover changes, known as PIFs, from the images acquired on different dates. Considering the important impact of PIFs on the estimation of the mathematical relationship between a reference-target image pair, some semiautomatic methods have been proposed to ensure the quality of PIFs, such as the ridge method [5], principal component analysis based method [26], and temporal invariant cluster method [7], but these mentioned methods still face the challenges in determining threshold or cluster centers. Alternatively, represented by multivariate alteration detection (MAD) [29], [30] and iterative slow feature analysis [31], mathematical transformation-based RRN methods have been proposed to extract PIFs using change probability. MAD is invariant to linear and affine scaling but is unsuitable for dealing with images containing substantial land cover changes [32], [33], [34]. As a result, many MAD variants have been proposed, such as iteratively re-weighted MAD (IR-MAD) [33], [35], [36], [37], [38], multitemporal MAD [32], and KCCA-based MAD [39], [40]. However, their performance is highly dependent on the threshold segmentation result of the change probability map.
Enhancing radiance consistency is also necessary for applications using Landsat time-series images. Even though Landsat surface reflectance code (LaSRC) [47] and Landsat ecosystem disturbance adaptive processing system (LEDAPS) [48] are namely used to generate Level-2 surface reflectance products for Landsat 8 OLI and Landsat 5 TM/7 ETM+ data, the surface reflectance products' uncertainty is still increased by the uncertainties of radiometric calibration [47] and different atmospheric correction algorithm. Moreover, adjusting the spectral reflectance of similar bands is unneglectable to the application using multisensor observations of Landsat [49]. RRN can be an effective means to address uncertainties from the above factors. However, most RRN methods are developed for bitemporal images. For time-series images, the separated image-to-image PIFs identification hardly keeps consistent PIFs and leads to less comparable normalization models among all image pairs [50]. Coupling this with the reference selection under unstandardized criteria increases the risk to yield unstable normalization results [10]. To address mentioned limitations, we proposed a multi-rule-based RRN method to enhance the radiometric consistency of Landsat time-series imageries. The contributions of this study include the following.  1) We developed a novel trend based PIFs identification method, which provides consistent PIFs for all referencetarget image pairs by identifying land cover changes and phenological trends from entire time series. 2) We proposed a novel local PIFs optimization, including an automatic reference image allocation and a PIFs refining using spatial-spectral metrics for each reference-target image pair. 3) We achieved a good normalization performance aided by a combined RRN modeling using inlier data pair identification and robust linear regression (RLR). The rest of this article is organized as follows. Section II introduces the experimental data. Section III describes the proposed method. Section IV reports and analyzes the results. Section V discusses the performance differences compared to seven commonly used RRN methods and the pros/cons of the proposed method. Section VI draws the concluding remarks.

II. DATA
Landsat Collection 2 Level-2 surface reflectance product serves as ideal data for our study. Compared to Level-1 digital number data, Level-2 surface reflectance improves comparison among multiple images over the same region in detecting the Earth's surface changes by accounting for atmospheric effects. We collected Landsat Collection 2 Level-2 product from the United States Geological Survey (https://www.usgs.gov/corescience-systems/nli/landsat) via Google Earth Engine (GEE) platform [51], including the atmospherically corrected surface reflectance data from Landsat 5 TM, 7 ETM+, and 8 OLI. The main properties of the sensors are listed in Table I. We used blue, green, red, near-infrared (NIR), and short-wave infrared (SWIR) 1/2 bands of each image after removing the cloud and cloud shadows pixels using the C Function of Mask (CFMask) algorithm [52]. We selected a rectangle experimental area with 1100 × 1100 pixels in Nanjing, China. This area has various land covers and has experienced drastic urbanization over the past three decades [53]. Based on the experimental area (WRS-2 Path/Row: 120/038), we filtered a total of 580 images acquired from 2001-01-01 to 2020-12-31 from the GEE image collections of "LANDSAT/LT05/C02/T1_L2," "LAND-SAT/LE07/C02/T1_L2," and "LANDSAT/LC08/C02/T1_L2."

III. METHODOLOGY
The proposed RRN method includes three key steps: 1) trend-based PIFs identification from entire Landsat time-series images; 2) PIFs optimization for reference-target image pair; and 3) RRN modeling (see Fig. 1). First, we identified trend-based PIFs from annual harmonic-fitted spectral indices belonging to the unchanged area using the Mann-Kendall (MK) trend test method (see Section III-A). This aims to identify consistent PIFs for all reference-target image pairs. Then, after the automatic reference allocation for a monthly image stack, we refined the trend-based PIFs for each reference-target image pair using a spatial-spectral metric (see Section III-B).
Finally, we normalized target images in each monthly image stack using a combination of the M-estimator sample consensus (MSAC) algorithm and RLR (see Section III-C). To take full advantage of different processing platforms, Section III-A and the reference selection in Section III-B were performed on GEE. The remainder was implemented using MATLAB 2021a.

1) Unchanged Area Detection Using GEE-CCDC:
To identify potential PIFs from a steady area without land disturbances, we detected the unchanged area using the change detection module of the continuous change detection and classification (CCDC) [54] algorithm. CCDC is a statistical boundary method to detect abrupt and gradual changes in different land covers using all available Landsat images [55]. In data preparation, the cloud and cloud shadow pixels in each image are first masked using CFMask. Then, CCDC uses iterative Tmask (multiTemporal mask) cloud detection [56] to further mask other missed outliers. In detecting changes, CCDC first estimates a time series model using a given number of clear observations (i.e., 12) for each pixel position [54]. Then, CCDC flags a possible change by comparing model predictions with clear observations. When the required consecutive times of the possible changes are reached, this pixel is assigned as a break. After that, CCDC continues to estimate a new model until the next break has been identified or all observations have been exhausted. Finally, the time series at each pixel position are split into multiple temporal segments to indicate the changes deviate from the previous pattern.
The first occurred changes among multiple changes detected by CCDC in each pixel were used to obtain the unchanged area, aided by the repository of GEE-CCDC-Tools (a GEE version of CCDC) [57]. The main parameters were set as follows: 1) the surface reflectance bands of blue, green, red, NIR, SWIR1/2, normalized difference vegetation index (NDVI) [58], normalized burn ratio (NBR) [59], and normalized difference fraction index [60] were used as "breakpointBands" for change detection; 2) the "chiSquareProbability" (chi-square probability threshold for change detection) was set as 0.90; 3) the "minObservations" (number of observations required to flag a change) was set as 6. Other parameters were set as defaults [61].

2) Harmonic Fitting of Time-Series Spectral Index:
We temporally interpolated the masked pixels caused by cloud and shadow to obtain gapless time-series indices of the unchanged area. First, we filtered all observations of the unchanged area and removed the incomplete time series with less than 20% of the total observations. Then, we calculated three time-series spectral indices of NDVI, NBR, and modified normalized difference water index (MNDWI) [62] to highlight the pixel variations with different land covers. Finally, we applied the Harmonic Analysis of Time Series (HANTS) [63] algorithm to fit each time-series index because HANTS can decompose a time-dependent periodic phenomenon into a series of sinusoidal functions defined by unique amplitude and phase values [64].

3) Trend-Based PIFs Identification Using the Mann-Kendall Test:
We identified trend-based PIFs by distinguishing the pixels with nonsignificant trends from the harmonic-fitted time-series indices. First, we used the MK trend test [65], [66], [67], [68] to detect the upward or downward monotonic trend with a statistical significance in each fitted time-series index. Then, we selected the pixels where there are no significant monotonic trends over the years as the candidate PIFs of each index. Finally, we integrated the spatial union of the candidate PIFs of NDVI, NBR, and MNDWI as the trend-based PIFs, denoted as PIF tre . PIF tre are considered to be consistent radiometric-invariant in the entire time series. In detail, we selected the harmonic-fitted index in February, March, and April as the inputs because the MK test is not suitable for the data with seasonality or periodicities [69]. At a significance level α of the test, if |Z| ≥ Z 1−α/2 (Z is the MK test statistic), the null hypothesis that there is no monotonic trend is denied. Namely, there is either an upward (positive) or a downward (negative) monotonic trend when |Z| ≥1.65, 1.96, and 2.576 (at the α of 10%, 5%, and 1%). The trend categories were listed in Table II.

B. PIFs Optimization for the Reference-Target Image Pair
We first allocated all images with cloud coverage of less than 30% to 12 monthly image stacks. Then, the PIFs optimization was used to refine the trend-based PIFs for each reference-target image pair in the monthly image stack. It includes an automatic reference image selection and a PIFs refinement using spatial and spectral metrics.
1) Automatic Reference Image Selection for Monthly Image Stack: We used a pixel-based image quality evaluation to automatically select the reference image in each image stack. The sensor score (S S ), opacity score (S o ), and distance to cloud/cloud shadow score (S d ) were calculated for each pixel [70]. The sum of three scores was used to represent the pixel quality. We defined the better pixel P b as the pixel with a score greater than mean score of the current image. The image with the most P b in each monthly image stack was selected as reference image.
The S S is an image-level score to avoid selecting references from Landsat 7 ETM+ SLC-off data. The pixels of normal ETM+ and other sensors' data were allocated a score of 1. The pixels of ETM+ SLC-off data were assigned a score of 0.5.
We calculated the S o using the "SR_ATMOS_OPACITY" band. It is an estimator of aerosol optical thickness from the blue bands of Landsat 5 TM or 7 ETM+ images using the dark dense vegetation method [48], [71]. In general, an opacity value (O i ) less than or equal to 0.1 refers to a clear pixel, an O i ∈ (0.10.3] is an average pixel, and an O i > 0.3 is a hazy pixel. We assigned S o of a clear pixel (O i < 0.2) and a haze pixel as 1 and null, respectively. The S o of a pixel with O i ∈ [0.20.3] was calculated as follows [70]: (1) where O i is the opacity value of the ith pixel, O max and O min are namely the maximum opacity value (i.e., 0.3) and minimum opacity value (i.e., 0.2). A dummy band with an opacity score of 0.25 was used to replace the null "SR_ATMOS_OPACITY" band of Landsat 5 TM/7 ETM+ data and represent the opacity value for Landsat 8 OLI data.
We calculated the S d using the "QA_PIXEL" band. We assigned a score of 1 to the pixels that are greater than 30 pixels away from the cloud/cloud shadow. The rest pixels having a distance within 30 pixels were assigned a score as follows [72]: where D i is the distance to cloud/cloud shadow of the ith pixel, D r is the required minimum distance (i.e., 30 pixels), and D min is the minimum distance of the given pixels (i.e., 0 pixels).

2) PIFs Optimization Using Spatial and Spectral Metrics:
We refined the trend-based PIFs for each reference-target image pair based on their spatial and spectral differences. In the spatial domain, we first calculated the structural similarity index measure (SSIM) [73] for each band between reference and target images. Then, we sorted the SSIM values of each band in descending order and selected the top 80% of highly similar pixels as the candidate spatial similarity masks of each band, denoted as M spa_n . Finally, we obtained the spatial similarity mask (M spa ) using the spatial intersection part of all M spa_n In the spectral domain, we first calculated the spectral change magnitude (m) between reference and target images using the change vector analysis method where R i and T i are the ith band of reference and target images, m is the change magnitude, and n is the number of bands. Then, we converted m into a 1-D vector and sorted it in ascending order, denoted as vm. Next, we calculated the cumulative moving average of vm under K intervals (CMA K ) . . . , where N is the number of pixels within the interval, K is the number of intervals (i.e., 500), and vm i is the ith sorted vm value. Finally, we used the forward-difference method for CMA K to identify the inflection with minimum differential values [74]. All pixels with differential values (DV i ) less than the inflection value (DV inflection ) were determined as the spectral similarity mask, denoted as M spe (see Fig. 11) After converting PIF tre , M spa , and M spe to binary images, we used their spatial intersection as the optimized PIFs (PIF opt ) for current reference-target image pair

C. Relative Radiometric Normalization Modeling
We estimated RRN models for each band using the reference and target pixels corresponding to PIF opt . To achieve better fitting performance for each band, we first identified inlier data pairs of each band using MSAC [75] algorithm. Then, we used RLR to model the relationship between the identified inlier data pairs of each band. The above-mentioned processes aim to reduce the negative influence of outlier data pairs on the model accuracy. Finally, we normalized each band of the target images using the estimated models.

1) Identification of Inlier Data Pairs Using MSAC:
MSAC is an improved random sample consensus (RANSAC) [76] algorithm. The cost function of RANSAC is shown as where d is the selected subset, D is the dataset, Loss is the loss function, Err is the error function of geometric distance, and M is the estimated model parameter. The loss function is shown as where e is the error. T is the error threshold to determine inliers and outliers. However, higher T can lead to poor estimations (all solutions have the same cost as all the matches would be inliers) [75]. Therefore, MSAC uses the redescending M-estimator as a loss function to reduce the influence of the threshold on the model. The new robust loss function is given as In this way, outliers still can have a fixed penalty but inliers are scored to the extent they fit the data. For the linear regression in this article, the minimum number of the subset (n) was set to 2, the inlier proportion d p was set to 90%, and the error threshold t i for the ith band was set as follows: where θ is a proportion constant that was set to 0.3, R i and T i are the vectors of the ith band's pixel value for reference and target image, respectively, and n is the pixel number of each band.
2) Relative Radiometric Normalization Modeling Using RLR: After identifying inlier data pairs, we estimated RLR models for each band using the iteratively reweighted least square (IRLS) algorithm. The IRLS iteratively calculates weights to determine the influences of each response value on the final parameter estimates [77]. A lower weight is assigned to the point farther from model predictions in the previous iteration. Then, the IRLS solves model coefficients using weighted least squares. We used Tukey's bisquare function [78] as the weight function (w), which is given by where r is "residual." c is a positive parameter set to 4.685.

D. Evaluation of the Proposed RRN Method
We evaluated the performance of the proposed RRN method by comparing differences between the original and normalized reference-target image pairs. First, the improvement of image brightness intensity differences was visually checked using checkerboard visualization. Second, three qualitative evaluations were performed: 1) the precision of the estimated linear regression models was evaluated using a tenfold cross validation with root-mean-squared error (RMSE) and coefficient of determination (R 2 ); 2) the method's capability of minimizing the radiance differences in each band was evaluated using both the RMSE variations; and 3) using the difference histogram adjustment of the unchanged area derived by GEE-CCDC.

A. Trend-Based PIFs
A total of 106 632 pixels were identified as trend-based PIFs after integrating candidate PIFs of NDVI (26 080 pixels), NBR (62 300 pixels), and MNDWI (23 660 pixels) (see Fig. 2). Within trend-based PIFs identification, a total of 655 814 changed pixels were detected from 2001 to 2020 using GEE-CCDC [see Fig. 2(a)]. Most changed pixels were converted from vegetation to impervious. Different amounts and spatial distributions of the candidate PIFs showed the index capability biases in providing the trend-based PIFs from different land covers [see Fig. 2(h) and (i)]. NDVI and NBR provided most candidate PIFs corresponding to vegetation and the impervious, but NBR provided more candidate PIFs from forests. Integrating these candidate PIFs effectively improved the representativeness of trend-based PIFs.

B. Performance of PIFs Optimization
Within the PIFs optimization, the reference image of each monthly image stack was identified by automatically selecting the image with the greatest number of P b (see Fig. 3). Taking the image stack of December as an example, the image "LC08_20161209" was selected as a reference. Figs. 3 and 12 showed that the proposed reference image selection is effective to distinguish the reference with the best pixel quality from the images covered by haze/cloud or acquired by faulty Landsat 7 ETM+.
After identifying the trend-based PIFs that intersect the spatial and spectral similarity masks (M spa and M spe ) (see Fig. 4), an average of 44 330 pixels were obtained as optimized PIFs for each reference-target image pair in the image stack of December [see Fig. 4(g)]. The PIFs optimization further reduced the potential errors in trend-based PIFs to the linear regression models, according to RMSE and R 2 of the linear regression models using the PIFs obtained at different stages (see Fig. 5). The trend-based PIFs achieved a general fitting performance: the models of visible bands have a lower RMSE but the median R 2 lower than 0.65; the models of NIR and SWIR 1/2 bands have higher R 2 from 0.76 to 0.82. The implementation of optimized PIFs significantly improved the models' fitting performance for each band: the median of RMSE was decreased to a lower range of 0.0083-0.0130, and the median of R 2 was increased to a higher range over 0.85 except the blue band has a large variation around 0.80. We found that these models of the blue band with unsatisfying fitting performance are from hazy/cloudy images acquired by Landsat 5 TM/7 ETM+, such as "LT05_20081219" ( R 2 = 0.60), "LE07_20121206" ( R 2 = 0.61), and "LE07_20151231" ( R 2 = 0.63) (see Fig. 12). These images contain the pixels without homogeneous transparency that negatively impacted the fitting precision. Thus, it is necessary to eliminate these negative impacts on each band for a more precise RRN model estimation.

C. Performance of the Proposed RRN
The evaluation results showed that the proposed RRN method achieved a good performance in model precision and radiance consistency improvement. In terms of visual evaluation results, the checkerboard visualization showed that our RRN method has effectively reduced the brightness intensity differences between raw reference and target images, resulting in more seamless checkerboard cells (see Fig. 6). The normalized target image sequence has more consistent color tones with the reference image instead of the distinct fluctuations in the raw target image sequence.
The quantitative evaluation proved the effectiveness of our RRN method from the following aspects.
1) The RRN modeling combining MSAC with RLR better estimated the linear regression models for each band (see Fig. 5). The model precision was significantly improved: the median RMSE of each band was decreased to 0.0032, 0.0029, 0.0035, 0.0033, 0.0039, and 0.0037; the median R 2 of NIR and SWIR1/2 bands' models was increased to a steady range of around 0.99, and the R 2 of visible bands' models achieved a satisfying performance with a median R 2 of 0.97, 0.99, and 0.98. The model precision of hazy/cloudy images was improved using MSAC and RLR, for instance, the R 2 of each band for "LE07_20121206" was increased to 0.78, 0.94, 0.94, 0.98, 0.99, and 0.99. 2) The proposed RRN method significantly reduced the radiance differences of the pixels in the unchanged area. We used RMSE to represent the radiance differences between the pixels of reference and target images; the mean RMSE of each band was decreased by 29.37%, 20.77%, 15.37%, 7.87%, 6.15%, and 4.82% on average (see Fig. 7). Moreover, the radiance differences in visible bands of TM/ ETM+ data were effectively improved, especially for the hazy images. For example, the mean pixel differences in visible bands for "LT05_20081219," "LE07_20091230," "LE07_20121206," and "LE07_20151231" were decreased by 75 Fig. 8).
3) The proposed RRN method increased the number of pixels with values around zero. Based on the difference histograms of the unchanged area between the reference and the target images before (D R−T ) and after (D R−N ) performing RRN (see Fig. 9), the histograms of D R−N symmetrically follow a more concentrated normal Gaussian distribution centered at zero than that of D R−T . Our RRN method also successfully adjusted some bimodal distributions in D R−T caused by inconsistent radiance [see Fig. 9(b)] and maintained the bimodal parts caused by underdetected changes [see Fig. 9(c)]. It proved the effectiveness and importance of RRN in improving change detection accuracy by reducing the probability of pseudochanges [5], [79].

A. Comparison of RRN Methods
We compared the normalized image stacks of December between our method and seven commonly used RRN methods. The comparative methods include the DRRN methods of HC, MM, MS, and SR, the SRRN methods of conventional PIFs normalization, IR-MAD, and spectral angle mapping-based PIFs (SAM-PIF) [74] normalization. We visually interpreted 1670 unchanged pixels for the impervious (270), vegetation (900), and water (500) as ground truth to compare the mean RMSE of that between the reference and normalized target images.
Our RRN method outperformed the comparative RRN methods on 18 of 27 images, obtaining the most consistent normalization performance with a decreased mean RMSE of 19.84% and the lowest mean RMSE of 0.0170 (see Fig. 10). The methods of MS, SAM-PIFs, IR-MAD, and SR also performed well with Fig. 6. Checkerboard visualization of raw (a, c, e) and normalized (b, d, f) reference-target image pair in the image stack of December (SWIR 1, NIR, and green bands are displayed in red, green, and blue channels). More seamless cells indicate fewer differences between reference and target images. the decreased mean RMSE of 18.09% (4/27), 16.63% (2/27), 15.83% (2/27), and 7.56% (1/27). The pixel differences in Landsat Collection 2 Level-2 surface reflectance products are relatively small since the standardized ARN has reduced parts of the differences introduced by atmospheric conditions. Thus, the mean RMSE derived using each method is closer compared to the research using raw digital number (DN) values data. Inevitably, our method did not achieve the lowest mean RMSE on all images but reached relatively closer values to some lowest mean RMSE obtained using other methods, i.e., 0.0231 for "LE07_20021227," 0.0093 for "LE07_20161217," and 0.0111 for "LC08_20171212." In detail, compared to the results derived using comparative methods, our RRN method is more advantageous in identifying accurate PIFs from hazy/cloudy images and leading to better local normalization results. For example, we used the OTSU method to automatically extract the unchanged pixels from the IR-MAD chi-square change probability map. As a result, a greater number of pixels were used to build a regression model. Like the global method of MS, it resulted in some normalized target images with the lower statistical RMSE but probably unfits the actual scenario: for image "LC08_20191202," IR-MAD hardly excluded the impact of clouds even though it achieved the lowest mean RMSE [see Fig. 13(c)]. On the contrary, our method filtered the proper PIFs by identifying the cloud coverage and land cover changes [see Fig. 13(e)-(g)]. Two checkerboards of the normalized image showed that our method generated a more seamless image with less local inconsistency [see Fig. 13(d) and (h)]. Moreover, the MS and SAM-PIF achieved the lowest mean RMSE in cloudy images "LE07_20121206"/"LE07_20151231" and "LC08_20141220," respectively, but our proposed method outperformed MS and SAM-PIF in the improvement of local brightness differences (see Fig. 14). Therefore, we still need to focus on the specific purposes to normalize the targets instead of using a single evaluation metric to determine the best normalization method.

B. Discussion on the Proposed RRN Method
The proposed RRN method is a multi-rule-based method including three steps in order: trend-based PIFs identification, PIFs optimization, and RRN modeling. This method has been successfully applied to enhance the radiometric consistency of the Landsat Collection 2 Level-2 surface reflectance products, while there are some advantages and disadvantages in each step that need to be discussed.

1) Advantages:
In trend-based PIFs identification, we defined that PIFs are the pixels with steady physical attributes without experiencing land disturbances and significant phenological trends over the entire time series. Our trend PIFs identification addressed the influence of fluctuating numbers and distribution of PIFs on estimating RRN models. It provides consistent and representative trend-based PIFs for all referencetarget image pairs using a time series analysis strategy rather than separately comparing the instantaneous changes in radiometric features over two records. Two parameters need to be paid attention to when generating the unchanged area using the GEE-CCDC: "chiSquareProbability" and "minObservations." Both two are sensitive to the determination of changes. A lower chi-square probability threshold can better tolerate commission errors [61], [80], [81], and a lower "minObservations" can increase the number of detected changes by using fewer consecutive observations to flag a change [54]. To better balance the commission errors and integrity of the unchanged area, we namely assigned 0.90 and 6 to "chiSquareProbability" and "mi-nObservations," respectively. The involved time series analysis in this section is implemented on the GEE platform, alleviating the massive pressure of data processing on local devices, for example, a 500M CCDC-fitted result and a 1100 × 1100-pixel MK test result require 15 and 60 min to be outputted to user's Assets.
The PIFs optimization was designed on the assumption that the optimal PIFs for the current reference-target image pair are a subset of trend-based PIFs with similar spatial structures and spectral properties. The experimental results showed that the PIFs optimization is robust to hazy or cloudy images since both the spectral differences sensitive to radiance distortion and spatial differences insensitive to radiance distortion are exploited to yield M spa and M spe (see Fig. 13). Therefore, these optimized  PIFs are all from clear observations and are invariant to the gradual changes in the entire time series and the abrupt changes on two records. The percentage to determine the pixels with higher spatial similarity is the only required parameter during this step. It was allocated as 80% due to the higher local SSIM values caused by the medium-resolution image.
The quality of the PIFs to estimate RRN models is improved via the above processes. Therefore, the MSAC and RLR can effectively identify proper inliers and estimate the regression parameters of a precise RRN model without deviating from the actual relationship.
2) Limitations: During generating the trend-based PIFs, outputting a larger-scale MK-test result using a longer time series needs to overcome the limitation of running capacity. It can be addressed by splitting the image into M × N subsets to output but with the cost of more time. Besides, in reference identification, most references of the monthly image stacks are from Landsat 8 OLI images after 2013 and earlier Landsat 5 TM images due to high-penalty S S for Landsat 7 ETM+ SCL-off data in identifying reference (see Fig. 3). It brings possible unbalanced normalization performance caused by the uneven distribution of reference in the time series. It is noted that the hazy or cloudy images are still challenging to deal with because the images have heterogeneous transparency and hardly follow the assumption of a global affine modification of radiometric values [82].
3) Applications and Future Works: The proposed RRN method can eliminate the uncertainties caused by radiometric calibration, atmospheric correction, and sensor differences. The normalized images can be used as inputs for time-series change detection, vegetation monitoring, and long-term land cover classification to generate more comparable results. The steps of trend-based PIFs identification and automatic reference image selection require Landsat Collection 2 Level-2 products as input due to the usage of surface reflectance and quality assessment bands. Once the trend-based PIFs are obtained and the reference image is determined, the PIF optimization and RRN can also be performed to normalize the Collection 2 level-1 data belonging to the same spatial domain and acquisition time series, such as the calibrated top-of-atmosphere reflectance and DN values data.
This study also opens future research avenues to explore an effective RRN method between the Landsat series and Sentinel 2 MSI (MultiSpectral Instrument) data, which can contribute to improving the consistency of the harmonized Landsat 8/9 OLI and Sentinel-2 MSI surface reflectance products.

VI. CONCLUSION
We proposed a novel RRN method, including a global trendbased PIFs identification, a local PIFs optimization, and an RRN modeling, for enhancing the radiometric consistency of Landsat time-series imageries. The Landsat Collection 2 Level-2 surface reflectance data from 2001 to 2020 of a 1100 × 1100 pixels area in Nanjing, China, were used to validate our method. The results showed the following.
1) The proposed trend-based PIFs identification effectively identified the PIFs pixels without land cover changes and phenological trends from the entire time series. 2) The PIFs optimization is capable of automatically allocating the reference image to each monthly image stack and optimizing the PIFs based on specific spatial and spectral differences of the current reference-target image pair.
3) The RRN modeling using MSAC and RLR achieved a good performance in model precision (median R 2 > 0.96 and median RMSE < 0.0039) and radiance consistency improvement, resulting in a more consistent normalized series over time. 4) Our RRN method outperformed seven commonly used RRN methods on 18 images in the stack of December. As an essential preprocessing step, the proposed RRN method can effectively eliminate the uncertainties caused by radiometric calibration, atmospheric correction, and sensor property differences. The normalized images can be used as more reliable inputs for time-series change detection, vegetation monitoring, and long-term land cover classification.

APPENDIX
A. Tables   TABLE II  TREND CATEGORY OF THE MK TEST RESULTS B. Figures   Fig. 11. (a) Cumulative moving average and (b) forward-differential values for the change magnitude of the image pair of reference "LC08_20161209" and target "LC08_20201220."