Assessment of Long-Term Sensor Radiometric ... - IEEE Xplore

2960

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 52, NO. 5, MAY 2014

Assessment of Long-Term Sensor Radiometric Degradation Using Time Series Analysis Wonkook Kim, Member, IEEE, Tao He, Dongdong Wang, Changyong Cao, and Shunlin Liang, Fellow, IEEE

Abstract—The monitoring of top-of-atmosphere (TOA) reflectance time series provides useful information regarding the long-term degradation of satellite sensors. For a precise assessment of sensor degradation, the TOA reflectance time series is usually corrected for surface and atmospheric anisotropy by using bidirectional reflectance models so that the angular effects do not compromise the trend estimates. However, the models sometimes fail to correct the angular effects, particularly for spectral bands that exhibit a large seasonal oscillation due to atmospheric variability. This paper investigates the use of time series algorithms to identify both the angular effects and the atmospheric variability simultaneously in the time domain using their periodical patterns within the time series. Two nonstationary time series algorithms were tested with the Landsat 5 Thematic Mapper time series data acquired over two pseudoinvariant desert sites, the Sonoran and Libyan Deserts, to compute a precise long-term trend of the time series by removing the seasonal variability. The trending results of the time series algorithms were compared to those of the original TOA reflectance time series and those normalized by a widely used bidirectional-reflectance-distribution-function model. The time series results showed an effective removal of seasonal oscillation, caused by angular and atmospheric effects, producing trending results that have a higher statistical significance than other approaches. Index Terms—Landsat 5, Libyan Desert, Sonoran Desert, time series analysis, vicarious calibration.

I. I NTRODUCTION

S

ATELLITE sensors are subject to the postlaunch degradation and aging of calibration devices [1]–[3] owing to their exposure to the harsh environment of space. Vicarious calibration, which refers to calibration techniques that do not depend on onboard calibration devices, is critical for quantitative remote sensing for which high measurement accuracy and long-term stability are necessary. Recently, many studies have investigated the long-term stability of satellite sensors in Manuscript received November 20, 2012; revised March 29, 2013; accepted May 30, 2013. Date of publication July 18, 2013; date of current version February 27, 2014. This work was supported in part by Geostationary Operational Environmental Satellite-R Series calibration funds from the National Oceanic and Atmospheric Administration and in part by the National Aeronautics and Space Administration under Award NNX09AN36G. W. Kim is with the Department of Geographical Sciences, University of Maryland, College Park, MD 20740 USA and also with the Korea Ocean Satellite Center, Korea Institute of Ocean Science and Technology, 425-600 Ansan, Korea (e-mail: [email protected]). T. He, D. Wang, and S. Liang are with the Department of Geographical Sciences, University of Maryland, College Park, MD 20740 USA (e-mail: the@ umd.edu; [email protected]; [email protected]). C. Cao is with the Center for Satellite Applications and Research, National Oceanic and Atmospheric Administration, College Park, MD 20740 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TGRS.2013.2268161

the context of relative calibration by using top of atmosphere (TOA) reflectance time series data acquired for pseudoinvariant calibration sites (PICSs). The satellite sensors validated with such methods include the Advanced Very High Resolution Radiometer (AVHRR) Metop-A [4], Moderate Resolution Imaging Spectroradiometer (MODIS) Terra and Aqua [5]–[7], Landsat 5 Thematic Mapper (TM) [6], Landsat 7 Enhanced TM [6], [7], and Along-Track Scanning Radiometer-2 [8]. The assessment of sensor degradation is usually performed after correcting the TOA reflectance for the angular variation caused by changes in sun–sensor–target geometry since the angular effects are not directly related to the sensor performance. In many studies [4], [5], [7], [9], the angular effects have been treated with semiempirical bidirectional reflectance distribution function (BRDF) models. The normalization of TOA reflectance by such directional models is, however, not always successful and fails to remove the angular effects, particularly for the spectral bands that are sensitive to water vapor absorption. For example, the time series of near-infrared bands (Bands 4, 5, and 7) in Landsat 5 TM shows that, in addition to the angular effects, it has a clear annual seasonal oscillation that occurs repeatedly over the entire time period of the data sets [9], [10]. Since regular BRDF models only take viewing geometry variables as inputs, angular effects are not corrected when there is a large seasonal oscillation in the time series. In this paper, a new data-driven normalization process is proposed based on time series algorithms that identify both types of variability (angular and atmospheric) using their temporal patterns in the time domain. The key idea is to identify shortterm variations in the time series based on their frequencies and to remove periodic variations that have frequencies relatively higher than the frequency of the long-term trend. The angular variation and the atmospheric variation have clearly shorter cycle periods (i.e., high frequencies) than the period of the longterm trend caused by sensor degradation. In the case of Landsat 5 TM, viewing geometric variables vary annually because the satellite has an almost nadir-view look angle, and the major changing factor for the reflectance is the solar zenith angle, which varies with seasons. For the atmospheric variability, a focus is placed on short-term variability such as the seasonal variation induced periodically by a local climate pattern at the target sites. A major difficulty in decomposing the time series of TOA reflectance into signals of different frequencies is the fact that the time series can be nonstationary in time because of the following: 1) There is a trend in the time series, and 2) the atmospheric variability itself is nonstationary, creating irregular periodic patterns every year.

0196-2892 © 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

KIM et al.: ASSESSMENT OF LONG-TERM SENSOR RADIOMETRIC DEGRADATION

Fig. 1.

2961

RGB composite images of (a) Sonoran Desert and (b) Libyan Desert.

Nonstationary time series algorithms are capable of decomposing such nonstationary time series into seasonal and trend components and have been used in cases where one needs to identify a long-term trend that is obscured in a time series with naturally occurring variability. Unlike algorithms with a stationary assumption such as the Fourier transform, nonstationary algorithms assume local stationarity and allow the secondorder structure of time series to vary over time. This property is suitable for analyzing potentially nonstationary time series data, the behavior of which is governed by a complex climate system. A primary example of time series algorithms related to sensor calibration is the detection of orbital drift artifacts in the normalized difference vegetation index (NDVI) product of AVHRR [11]. Artifacts associated with orbital drift are successfully extracted from the NDVI time series that exhibit high seasonal variation caused by vegetation phenology, using a time series algorithm without involving any physical models that simulate the orbit dynamics. Among many nonstationary time series algorithms, two are used in this paper to decompose the time series data for sensor degradation assessment: 1) seasonal trend decomposition procedures based on loess (STL) [12] and 2) discrete wavelet transform (DWT) [13]. The algorithms have shown the effective extraction of nonstationary seasonal oscillation from the time series with an increasing or decreasing trend, in many environmental studies such as the investigation of El Niño effects on vegetation variability [14], correlation analysis between NDVI and meteorological variables [15], and detection of trend and seasonal changes in MODIS time series data [16]. In this paper, several issues are addressed that occurred during the application of these algorithms to the degradation assessment of the Landsat 5 TM sensor, particularly for the data sets that were acquired over the two popular PICSs in desert areas. Although the proposed approach was tested with the Landsat TM data, the principle of the approach can also be used for the time series data of other polar-orbiting and geostationary sensors. This paper is outlined as follows. In Section II, the characteristics of the two desert calibration sites are briefly introduced. In Section III, the specification of the Landsat 5 TM sensor is presented with an emphasis on its calibration status, and the cloud screening procedure applied to the data sets is explained. In Section IV, the interpolation procedures applied prior to the

time series algorithms are described. Section V contains the formulations of the two time series algorithms, and Section VI contains the decomposition results of the two algorithms. In Section VII, the trend results of the time series algorithms are presented in comparison with those of other time series such as the original and BRDF-normalized time series. Section VIII concludes this paper with a proposal of future direction. II. P SEUDOINVARIANT S ITES The stability of the target sites is important for the long-term monitoring of sensor performance. The characteristics of the two popular PICSs used in this study are briefly described in the following. A. Sonoran Desert The Sonoran Desert is a large flat area of sand located near Yuma, Arizona, and the U.S.–Mexico border (31◦ 15 N− 32◦ 40 N, 114◦ 55 W−113◦ 15 W, WRS-2 path/row: 38/38). The desert area has a mean elevation of 0.2 km above sea level and is covered with sand dunes. The area has a biseasonal rainfall pattern. Frontal storms originating from the North Pacific Ocean sporadically bring gentle rain over the broad desert area from December to March, while the summer monsoon brings localized severe thunderstorms from July to mid-September. Temperatures may exceed 48 ◦ C in the hottest season, with surface temperatures reaching around 80 ◦ C. Humidity is less than 10% most of the time. Drought-tolerant desert vegetation such as Cresoto Bush (Larrea divaricata tridentate) and White Bursage (Ambrosia dumosa) are sparsely distributed over the entire region. To select an optimal location for longterm monitoring, a moving window experiment has been performed, where approximately 70 grid points were generated over the whole desert area, in which the radiometric uncertainty was computed for each grid point. The selected location (31◦ 55 N, 114◦ 32 W) [see Fig. 1(a)] is close to the area that had been identified in previous studies [7], [17]. B. Libyan Desert The location of the calibration site in the Libyan Desert, which is often referred to as Libya 4 (28◦ 33 N, 23◦ 23 W, WRS-2 path/row: 181/40), is a large area of sand dunes located

2962


in North Africa [see Fig. 1(b)]. The site has good long-term radiometric stability according to the stability analysis performed in [18]. The area is at an elevation of 118 m above sea level and lacks vegetation. The Libyan Desert site is widely used for the vicarious calibration of various satellite sensors due to its uniform spatial texture, excellent temporal stability, low aerosol loading, and minimal cloud cover with low annual precipitation [18]. III. S ENSOR AND C LOUD S CREENING A. Landsat 5 TM Landsat (LS) 5 TM data products have undergone continuous updates of the calibration process validated by independent cross-calibration and vicarious calibration activities [1], [19]. Initially, the calibration of the TM reflective bands was performed using three lamps loaded in the internal calibrator (IC), in which the radiometric gain and offsets could be computed on a scene-by-scene basis. However, related studies discovered that, due to their instability, the regression calibration based on IC lamps could not produce the actual detector gains for most of the mission. A new calibration process based on lifetime radiometric calibration curves was proposed in 2003 [1], [20] and was applied to the entire collection of the LS 5 TM data to remove the effects of the unstable calibration lamps. Another update was practiced in a similar manner in 2007 [19]. The major improvements by the new calibration procedure are primarily for data sets acquired during the first eight years of the mission (1984–1991), where the sensor gain changed by around 15%. The radiometric uncertainty of the sensor made less of an impact in the later years than in the earlier years. Since the lifetime adjustment is implemented in the product generation stage, end users can perform radiometric calibration by using the bias and the gain provided in individual Landsat metadata files (.MTL). In this paper, 13-year (1999–2011) and 7-year (2005–2011) level 1 products of LS 5 TM data were collected for the Sonoran Desert and Libyan Desert, respectively. The digital-number counts were converted to TOA reflectance using the radiometric calibration equation [21]. Conversion to TOA reflectance serves as normalization of data because it accomplishes the following: 1) removes the cosine effect of different solar zenith angles; 2) divides the radiance by in-band exoatmospheric solar irradiance computed from a spectral response function of each band (Thuillier’s model [22] is used for the solar irradiance spectrum); and 3) adjusts for the variation in the Earth–Sun distance between different data acquisition dates. Note that the surface and atmosphere anisotropy remain the factors causing seasonal variation in the TOA reflectance time series. With nearly nadir-view look angles, the LS 5 TM data have an advantage in monitoring the long-term variability of the site, minimizing the effect of the sensor’s azimuthal variation. In this paper, the long-term degradation of the LS 5 TM sensor was assessed for the six solar reflective bands—three visible bands (Bands 1, 2, and 3) and three near-infrared bands (Bands 4, 5, and 7). The relative spectral response functions of the reflective bands are plotted in Fig. 2 with a typical TOA reflectance signature of the desert area, which is retrieved from a Hyperion data set.

Fig. 2. Relative spectral response function of LS 5 TM (scaled by 0.5) and TOA reflectance signature of typical sand sea area in the Sonoran Desert.

B. Cloud Screening Cloud screening is a critical processing step in the evaluation of the stability of a calibration site because cloud pixels are a major source of abnormal observations or outliers in time series. Cloud coverage information is available in individual Landsat TM scenes, but per-pixel information (i.e., cloud mask) is not available in the current products. The automatic cloud cover assessment (ACCA) algorithm [23] is used to compute cloud mask automatically in every scene. The ACCA algorithm employs a two-pass algorithm; the first pass involves a series of filters that are designed to collect highly probable cloud pixels, and the second pass refines the results, minimizing the confusion with sand, vegetation, and snow by using the statistics of the cloud pixels collected in the first pass. The Landsat 7 Enhanced Thematic Mapper Plus ACCA algorithm proposed after the ACCA algorithm was not used since it is not capable of producing per-pixel cloud mask. We also compared the ACCA results with those obtained by the modified Luo–Trishchenko–Khlopenko algorithm [24] and found that the ACCA results are comparable or better. IV. I NTERPOLATION OF M ISSING DATA Missing data can occur commonly in the time series of Landsat data because of the following: 1) Scenes are not always routinely acquired for areas outside the U.S. continent; 2) acquired data are sometimes not archived due to unsuccessful downlink and processing; and 3) scenes with high cloud cover cannot be used for the analysis. Seasonal oscillation is a major difficulty in the interpolation of TOA reflectance time series over desert areas. Since the analytic form of the time series curves is not known a priori, it is difficult to estimate the nonlinear seasonality, particularly when data samples are missing consecutively for several time periods. For example, the biseasonal precipitation pattern in the Sonoran Desert sites sometimes makes data unavailable for 4–5 consecutive acquisition dates, making the seasonal oscillation pattern occluded


Fig. 3.

2963

Seasonal interpolation results for a missing entry in the time series of Band 4 for (a) the Sonoran Desert and (b) the Libyan Desert.

for that time period. To overcome the problem, a seasonal interpolation scheme was implemented which can better reflect the seasonality in the TOA reflectance time series. The basic idea is to interpolate the values of missing entries using a set of samples from a similar time of each year. Such collections of data are referred to as “cycle subseries.” For example, for daily data for 10 years (cycle period 365), there are 365 cycle subseries, and the first cycle subseries will be the collection from January 1 data each year, comprising 10 data samples. It is reasonable to assume that the data samples in the cycle subseries are not subject to seasonality, although there may be an interannual trend. First, one cycle subseries is constructed for each missing entry. If the date of the missing entry is dmiss , the cycle subseries Cdmiss consist of data from dates dmiss , dmiss + p, and dmiss + 2p, where p is the cycle period. For Landsat 5 TM data, p = 365.2425. A linear regression analysis was performed for the samples in the cycle subseries to estimate the values of the missing entries. A practical issue in the application of the seasonal interpolation method to our Landsat time series data is that the revisit period of Landsat data is 16 days and the image acquisition date for the following years may not fall on exactly the same day of the year as it did in the first year. To overcome this problem, the constraints of the cycle subseries were relaxed by allowing data samples close to the exact time period of each year to be eligible for the subseries. Such cycle subseries with the relaxed condition are referred to as “pseudo cycle subseries.” For a missing entry with a date dmiss , a pseudo cycle subseries for ∗ the date, Cdmiss , is defined as h ∗ Cdmiss = xd ||di − d| < 2 di = dmiss + i · p,

i = 0, 1, 2, . . .

(1)

where xd is the data sample of the date d, h is a window size for the relaxation, and di are the dates on the exact cycle period from the data xd . The deviation of the acquisition date of d from the date on the exact cycle period is v ≡ |di − d|.

(2)

As the samples in the pseudo cycle subseries are now not on the same dates each year as they are in the exact cycle period, the samples can no longer be considered as annual replications of the missing data. That is, the larger the deviation of the date from the exact cycle period, the less likely it is that the samples belong to the same season. A weighted least square regression is used to consider the effect of the differences. When the deviation from the exact cycle period is v, weights are defined as w = 1/(1 + v)

(3)

to make the influence of data on the exact cycle period largest (v = 0, w = 1) and that of the completely off-period samples insignificant (v = ∞, w = 0). Prior to the application of the weighted least square regression, outliers are removed to obtain more robust regression results. The standard deviation of the pseudo cycle subseries, ρc∗ , is first computed, and the samples that have deviation (from the mean value) of more than m · ρc∗ (m is data dependent) are considered as outliers and are then excluded from the regression analysis. Interpolation results are illustrated in Fig. 3 for both the Sonoran and Libyan Desert data. Because the Libyan Desert data have more missing data during the winter season, a higher m for outlier detection and a larger window size w are used (m = 2, w = 48). The parameters used for the Sonoran Desert are m = 1 and w = 16. Fig. 3(a) and (b) shows the results of the seasonal interpolation for one set of missing data in Band 4 for the Sonoran Desert and Libyan Desert, respectively. The samples in the pseudo cycle subseries are plotted as blue circles, and those that are actually used for the weighted regression are marked with an additional cross. Outliers are marked as empty circles. The interpolated missing data are colored in red with a cross enclosed in a circle. The regression lines are plotted in green. Both figures show that the samples in the pseudo cycle subseries have a decreasing trend over the entire time period, but no apparent seasonality is found. The time series curves in the original time domain are presented in Fig. 4 for the Sonoran Desert (a) and the Libyan Desert (b), respectively, where the missing entries are filled in with the interpolated values.

2964


Fig. 4. TOA reflectance time series whose missing entries are interpolated by the proposed approach are presented for Band 4 of (a) the Sonoran Desert and (b) the Libyan Desert. Interpolated values are marked as red crosses.

After this seasonal interpolation was made, in order to achieve a better interpretation of the time series analysis results, another run of interpolation was performed to make the current time series have a period of natural number. For LS 5 TM data with a revisit period of 16 days, the number of samples in a year is 22.828. A simple linear interpolation was performed to make the number of samples per year in a newly interpolated time series the greatest natural number less than the current one. The new time series now has the number of samples per year equal to 22, with a revisit period of 16.602 days.

Step 3)

V. T IME S ERIES A LGORITHMS

Step 4)

The STL method is a nonlinear time series analysis technique that uses a series of filtering procedures to decompose the series into three additive components: trend, seasonal, and remainder components [12]. A given time series (Yt ) is represented as three additive components: trend (Tt ), seasonal (St ), and remainder (Rt )

Step 5)

A. STL

Y t = Tt + St + R t .

Step 6)

(4)

The components in the STL model are computed through two recursive procedures: an inner loop and an outer loop. The inner loop is mainly dedicated to separating the seasonal and trend components from the original time series, whereas the outer loop that nests n(i) passes of the inner loop is incorporated to identify and remove abnormal observations that may possibly hamper the inner loop results. Each pass of the inner loop consists of six steps as follows. Step 1) Detrending step. The time series Y is detrended by subtracting the trend component, Y − T(k) , where k denotes the number of the pass. The initial trend component is set to be T(0) = 0. Step 2) Cycle-subseries smoothing step. Cycle subseries refer to a subseries of Y , which is a collection of values from the same timing in each period. In this

step, each cycle subseries of the detrended series is smoothed by using a loess filter with q = n(s) and d = 1, yielding the collection of smoothed values in all the cycle subseries, C(k) . Low-pass filter of smoothed cycle subseries. A low-pass filter is applied to C(k) to compute the trend component. The low-pass filter consists of two consecutive runs of a moving average of length n(p) , followed by another run of a moving run with an average of length 3, and a LOESS smoothing with d = 1 and q = n(l) , producing the low-frequency output series L(k) . Detrending of smoothed cycle subseries. The lowfrequency component L(k) is subtracted from the cycle subseries so that the seasonal component may not be affected by the low-frequency factor, S(k) = C(k) − L(k) . Deseasonalizing. The seasonal component S(k) is subtracted from Y to obtain a deseasonalized series Y − S(k) . Trend smoothing The trend component T(k) is computed by smoothing the deseasonalized series using a loess with q = n(t) and d = 1.

The outer loop is designed to increase the robustness of the STL algorithm by assigning smaller weights to outliers. After each run of the inner loop ends, samples with large remainder values are identified and assigned a small or zero weight. The weights are then applied to the smoothing procedure in Steps 2) and 6). The overall structure of the procedures is presented in a flowchart in Fig. 5. The parameter selection procedure is suggested in [12], and the selection results are presented in Section VI with the decomposition results.

B. DWT DWT decomposes a vector of observations, Yt , using a set of basis vectors that are localized both in location and in scale.


Fig. 5.

Flowchart of STL algorithm.

The decomposition is essentially expressed as a weighted sum of basis vectors as Yt − Y¯ 1 =

nj L

Dj,k φj,k

(5)

j=1 k=1

where φj,k denotes the wavelet basis vectors; Dj,k denotes the wavelet coefficients for scale level j and location index k, Y¯ is the mean of the input series Yt , and 1 is a vector of ones. The wavelet basis vectors φj,k have N elements as the observation vector Y and are uniquely defined for a level and a location. With Y being the samples in the first level (j = 1) at N locations, the number of samples is halved with each step up in the scale, yielding a total of L = log2N levels and nj = 2−j N locations in each level. An important characteristic of the wavelet basis vectors is that they are orthonormal across both locations and scales, satisfying φj,k , φJ,K ≡

N

(i)

(i)

φj,k φJ,K

i=1

= 1 if j = J

and

2965

data have been resampled to a new grid system that has 22 samples per year in the preprocessing step. There are suggested values for the number of passes of the inner loop, n(i) , and the number of outer loop iterations, n(o) , which provide a near certainty of convergence: n(i) = 1 and n(o) = 10. The smoothing parameters, n(l) , n(s) , and n(t) , need to be determined carefully to properly identify each component. Although the selection of the smoothing parameters poses intrinsic ambiguity, reasonable selection criteria are proposed in [12], based on the eigenvalue analysis and the frequency response functions. The low-pass smoothing parameter n(l) can be set as the least odd integer greater than or equal to n(p) . In our case of Landsat data, where n(p) = 22, the optimal number of n(l) is 23. The choice of n(s) significantly depends on the characteristics of the series. The intrinsic ambiguity in the definition of seasonal variation makes this selection procedure rather subjective, requiring the use of a diagnostic method that is based on a seasonal-diagnostic plot. For a cycle subseries, s(k) represents mean values of seasonal components of the kth subseries. Fig. 6 shows the seasonal-diagnostic plots for four different seasonal smoothing parameters, with n(s) for the first cycle subseries in the first band of the Sonoran Desert data. The seasonal-diagnostic plot consists of two plots: 1) the seasonal component of the subseries subtracted by s(k) (continuous line in Fig. 6) and 2) the seasonal component plus the remainder component, also subtracted by s(k) (circles in Fig. 6). The figure shows that n(s) = 3 is a rather small value that causes overfitting and that n(s) = 15 does not pick up the overall trend pattern well. Parameter values between 5 and 9 appear to be proper for the seasonal smoothing. The smoothing parameter n(s) is set to 7 for both the Sonoran and Libyan Desert data. The selection of the trend smoothing parameter n(t) is also important because exceedingly large n(t) values make major low-/middle-frequency effects enter into the remainder component, whereas too small n(t) values make the short-term trend change and the seasonal variation indistinguishable. A proper trend smoothing parameter can be decided based on the following formula:

k=K

= 0 otherwise

(6)

where φ(i) indicates the ith elements of the vector. The construction of wavelet basis vectors and computation of the corresponding wavelet coefficients are performed based on the multiresolution analysis [25]. In this paper, Symlet-8 wavelet is used, which has yielded successful results for atmospheric time series data [26]. VI. D ECOMPOSITION R ESULTS A. Seasonal Trend Decomposition Based on Loess (STL) Parameter Selection: The six free parameters in STL are determined based on the approach presented in [12]. The selection of the number of samples in each period, n(p) , is quite straightforward. Since the period of primary focus in this study is one year, n(p) is set to 22. Note that the 16-day-interval Landsat

n(t) ≥

1.5n(p) . 1 − 1.5n−1 (s)

(7)

The trend smoothing parameter n(t) is set as 41 both for the Sonoran Desert and Libyan Desert. Decomposition Results: The decomposition results are presented in Figs. 7–10 for the six spectral bands of the Sonoran Desert and the Libyan Desert data. In Fig. 7, the seasonal and the trend components for the Sonoran Desert are plotted over the TOA reflectance time series of each band. Despite the interannual variation in the oscillation pattern and the nonzero trend components, the seasonality in the time series is extracted in a robust way, producing fairly constant oscillation curves each year. For each band, the seasonal curves of all the years are overlaid in a single plot (see Fig. 8), which illustrates that the shapes and magnitude scales remain more or less constant over the 13 years, although the patterns for each band differ.

2966


Fig. 6. Seasonal-diagnostic plots for the first cycle subseries of the Sonoran Desert data (Band 1). Plots for four different seasonal smoothing parameters are presented: (a) n(s) = 3, (b) n(s) = 5, (c) n(s) = 9, and (d) n(s) = 15.

The seasonal patterns of the six spectral bands can be grouped into three types based on the fluctuation pattern of the curves. The seasonal curves of Band 1, which consist of the first group, have quite a different pattern than that of the other spectral bands, with the reflectance peaking at two seasons, summer and winter, unlike that of the other spectral bands. The other visible bands (Bands 2 and 3) belong to the second type, which has a peak season in early June and a minor reflectance rebound in the beginning of October. The third type, a group of near-infrared bands (Bands 4, 5, and 7), has a peak timing earlier than that of the visible bands (the peak is at May), which is followed by a dramatic decrease in reflectance in August and September. The sudden decrease in reflectance seems to be due to the influence of water vapor in the atmosphere, possibly associated with the summer monsoon occurring between July and mid-September in the Sonoran Desert region. The reflectance decrease in the late summer season can also be observed in Bands 2 and 3, although the effect is much smaller than that of the near-infrared bands. Clear trends exist in the visible bands while the near-infrared bands have a fairly flat trend component. Complete trending results for all the other time series with different processing methods are presented in the next section in more detail. STL decomposition results for the Libyan Desert are presented in Figs. 9 and 10. Data availability of the Libyan Desert

is much poorer than that of the Sonoran Desert, with a shorter data time period (7 years), rendering greater difficulty in the long-term trending analysis. Band 5 data are not available for the Libyan Desert site, as the radiance counts are saturated for the desert. Figs. 9 and 10 show that fairly regular seasonal patterns are obtained with robust trend components for all the bands, despite the extremely adverse data availability. As in the results of a previous study [10], the magnitude of the seasonal variability is smaller than that in the Sonoran Desert for most of the bands. The reflectance decrease in August and September is not as dramatic as in the Sonoran Desert data, implying that the effect of water vapor on the near-infrared bands is smaller in the Libyan Desert. Note that the Sonoran Desert site is geographically more sensitive to atmospheric water vapor, being adjacent to the Pacific Ocean and the Gulf of California. B. DWT Selection of Hierarchical Levels for Seasonal Component: DWT decomposes a time series into signals of different hierarchical levels. The Sonoran Desert time series are decomposed into 9 levels, and the Libyan Desert data are decomposed into 8 levels. Fig. 11 shows the decomposition result of Band 1 of the Sonoran Desert. Signals of shorter time periods are contained in the higher levels while the longer time period signals are in the


Fig. 7.

STL decomposition results for the Sonoran Desert.

Fig. 8.

Seasonal components extracted by STL: Sonoran Desert.

lower levels (the signal of the annual period is in Level 5). Since DWT does not explicitly provide a seasonal component as defined in the STL results, it has to be constructed by selecting levels that are relevant to the component. As the purpose of using time series algorithms is to remove periodic variability

2967

that might compromise the trending analysis for the composition of the seasonal component, the signals of Levels 1 and 2 were excluded from the construction of the seasonal component because the monotonic (Level 1) and unimodal (Level 2) signals can be directly related to the sensor degradation. The signals of

2968


Fig. 9. STL decomposition results for the Libyan Desert.

Fig. 10. Seasonal components extracted by STL: Libyan Desert.

Levels 3–7 were included since they have cycle periods larger than 6 months. Visual inspection of the Libyan Desert data leads to the same selection of levels. Decomposition Results: Decomposition results of DWT are presented in Figs. 12 and 13 for the Sonoran Desert and in

Figs. 14 and 15 for the Libyan Desert, respectively. The figures show that the decomposition results of DWT have trend components with a generally small fluctuation over the whole time period. Less fluctuation in the trend components is attributed to the fact that the seasonal component of DWT takes signals


2969

Fig. 11. DWT decomposition results of Band 1: Sonoran Desert.

that have a cycle period of longer than 1 year, preventing longer term annual variation (e.g., signals with 2-year/4-year period) from entering into the trend components. For example, the Band 1 result of the STL algorithm in Fig. 7 shows that there is local reflectance decrease at around day 4000, which makes the trend component of STL fluctuate around this time. However, the local fluctuation is not observed in the trend component of Band 1 of the DWT results (see Fig. 12) because it has been already captured in the signals of Levels 3 and 4 which have been removed from the time series as a part of the seasonal component. The seasonal curves of DWT are smoother, and the interannual variation of the components is significantly larger than that of STL both in the Sonoran Desert (see Fig. 13) and in the Libyan Desert (see Fig. 15). There are two reasons for the larger interannual variation in the DWT results. First, as explained earlier for the trend component, the inclusion of the long period (> 1 year) signals in the seasonal component makes the annual plot contain signals that are not necessarily of a period of 1 year or less. Second, DWT is more adaptive to the local nonstationarity of time series data, as it allows the annual periodic pattern to have varying magnitude and shape. This can be observed in the DWT decomposition results in Fig. 11, where the periodic patterns, particularly in Levels 3–5, have different oscillation magnitudes for each period. This adaptivity of DWT is useful for explaining the local reflectance variation that can be induced by irregular events such as El Niño, but this may not be desirable for data sets that have a number of missing

entries. The DWT decomposition results of the Libyan Desert (see Fig. 15) show that the interannual variation of the patterns is so high for the visible bands that the overall fluctuation patterns are significantly different from each other. For the Sonoran Desert data, the overall annual fluctuation patterns extracted by DWT are relatively consistent over the whole time period, although their intensities vary each year. VII. L ONG -T ERM T REND A NALYSIS The long-term trend of TOA time series is estimated by performing a linear regression analysis on the time series. A deseasonalized time series was constructed from the time series algorithm by subtracting the seasonal component from the initial time series data. The trending results of the deseasonalized time series from STL and DWT were compared to those of (a) the initial TOA reflectance time series and (b) the BRDF-normalized TOA reflectance time series. For the BRDF normalization, only the results of the Staylor model are presented here, although the experiment was performed with the three models, Staylor [27], Roujean [28], and RossLi [29]. This is because the normalization results were not significantly different from each other for the data sets used in this study [10]. A time series of those four processing types (original, BRDFnormalized, STL, and DWT) are illustrated for a near-infrared band (Band 4) of the Sonoran Desert in Fig. 16, with an estimated linear trend. The deseasonalized time series was plotted

2970


Fig. 12. DWT decomposition results for the Sonoran Desert.

Fig. 13. Seasonal components extracted by DWT: Sonoran Desert.

with interpolated entries. The plots show that the seasonal oscillation in the initial TOA reflectance time series was drastically decreased in the deseasonalized time series, which was not possible with the BRDF models. Root-mean-square errors (RMSEs), a measure of the variance caused by the seasonal

oscillation, are significantly reduced in the deseasonalized time series compared to the original TOA and the BRDF-normalized TOA time series (original: 9.74 × 10−3 , BRDF-normalized: 9.37 × 10−3 , and STL/DWT: 4.07 × 10−3 ). Complete trending results for all the spectral bands and the processing types for


2971

Fig. 14. DWT decomposition results for the Libyan Desert.

Fig. 15. Seasonal components extracted by DWT: Libyan Desert.

the Sonoran Desert and the Libyan Desert are presented in Tables I–III, respectively, together with the results of the statistical test for the slope. Linear regression is based on y = α · z + y0

(8)

where y is the TOA reflectance, α is the slope, z is the number of days since the first date of the beginning year, and y0 is the y-intercept of the regression line. A statistical test was performed for the slope α with 95% significance level, and the corresponding t-statistics and p-values are presented in the

2972


Fig. 16. Trending results for Band 4 of the Sonoran Desert: (a) Original TOA reflectance, (b) BRDF-normalized, (c) deseasonalized by STL, and (d) deseasonalized by DWT. TABLE I WAVELENGTH R ANGE AND THE C ENTER WAVELENGTH OF THE LS 5 TM BANDS

tables. Relative radiometric change in TOA reflectance was measured as a ratio of the decreased/increased amount during the whole time period to the initial TOA value y0 . Annual changes are calculated by dividing the total radiometric changes by 13 (for the Sonoran Desert) and 7 (for the Libyan Desert). In general, all the time series presented exhibit very small relative changes. The estimated annual trend (%/year) is at most 0.41% (Band 1 in the Libyan Desert with original TOA reflectance) and is generally smaller than 0.2% for the other bands of both data sets, reflecting the excellent calibration status of the TM sensor. The changes in reflectance for the deseasonalized time series are generally comparable to or slightly smaller than those of the original and the BRDF-normalized

TOA reflectance, which indicates that the removal of seasonal effects was conducted in a robust way without causing artifacts. A major advantage of the deseasonalization is that the slope estimates have a higher statistical significance than what they had before the deseasonalization. Trending results from the two time series algorithms, STL and DWT, tend to have smaller p-values than the other results, along with smaller RMSE values. For example, before the deseasonalization, the slope estimates of Bands 4 and 7 of the Libyan Desert data could not be said to be significant (p-values: 0.084 and 0.141 for original TOA time series), but they turned to “significant” (p-values < 0.05) after the removal of the seasonal oscillation. Trending results of the two test sites are generally consistent, producing annual changes less than 0.1% for Bands 2, 4, and 7 and larger degradation (> 0.1%) for Bands 1 and 3. To assess the impact of the water vapor absorption on the atmospheric variability, the transmittance coefficients of each scene were calculated by using National Centers for Environmental Prediction (NCEP) Reanalysis 2 data provided by the Physical Science Division in the National Oceanic and Atmospheric Administration [30]. The global 4-times daily precipitable water data were retrieved for 1999–2011, and the data corresponding to the calibration sites were extracted from the entire grid points which are interleaved by 2.5◦ in latitude


2973

TABLE II T RENDING R ESULTS FOR S ONORAN D ESERT

TABLE III T RENDING R ESULTS FOR L IBYAN D ESERT

and longitude. The temporally and spatially closest data were used to estimate the water vapor amount in the atmosphere at the time of each Landsat 5 data acquisition. The water vapor transmittance coefficient for a given air mass and water vapor amount can be computed based on the formula in [31]. Since the formula in [31] oversimplifies the relationship between the transmittance and water vapor amount for the data sets in this study, a second-order term has been added to better characterize the relationship τγ2 = a0 + a1 · X + a2 · X 2 X = log(γ · η)

(9) (10)

where τγ is the transmittance coefficient, γ is the water vapor amount in g/cm2 , η is the optical air mass at the surface, and a0 , a1 , and a2 are band-dependent coefficients. The atmospheric mass at 101.3 kPa is obtained by η=

35 1224 · cos(θ0 )2 + 1

(11)

where θ0 is the solar zenith angle. The Second Simulation of a Satellite Signal in the Solar Spectrum vector code (6S) [32] is repeatedly run to generate pairs of the variables (τγ , η, γ, and θs ) for each spectral band of Landsat 5 TM. The coefficients, a0 , a1 , and a2 , are then regressed by using the variable pairs generated by 6S. Now that (9) has been established, the

2974


Fig. 17. TOA reflectance curves for Band 4 of the Sonoran Desert before and after the water vapor correction, overlaid with the NCEP precipitable water data.

and the right column shows the curves after the correction. The comparison shows that the water vapor correction significantly reduced the abrupt decreases in the late summer season in all near-infrared bands, suggesting that the decrease patterns is clearly associated with the water vapor variation. However, even with the smaller magnitude, the drops in the late summer still remain in the seasonality curves, which may be attributed to the following: 1) temporal and spatial mismatch between the reanalysis data and actual water vapor amount that the incident and the reflected rays in each Landsat scene would have gone through and 2) effects of other precipitation-related factors such as phenology of residual desert plants. The RMSEs of the BRDF-normalized and the water vapor-corrected TOA reflectance for the Sonoran Desert are 2.27, 2.80, 3.56, 6.11, 9.92, and 12.2(×103 ) for Bands 1, 2, 3, 4, 5, and 7, respectively. The RSMEs of the near-infrared bands are still larger than those of STL and DWT but are significantly smaller than those before the water vapor correction, demonstrating that the water vapor absorption is a major factor in the failure of the BRDF normalization for the near-infrared bands.

VIII. C ONCLUSION Fig. 18. Seasonality curves derived from the STL results of the Sonoran Desert (left column) before the water vapor correction (right column) after the water vapor correction.

transmittance at each Landsat acquisition time can be calculated by using the NCEP precipitable water data. Once the transmittance is calculated, the TOA reflectance of each scene is corrected for water vapor absorption by dividing the uncorrected TOA reflectance by the transmittance. Fig. 17 shows TOA reflectance curves of Band 4 (Sonoran Desert) before and after the water vapor correction, overlaid with the time series of precipitable water from NCEP reanalysis data. The precipitable water curve clearly exhibits a seasonal oscillation, where the peaks are in the late summers and the drops are in the winter season. The time of abrupt reflectance drops generally matches the timing of the surge in precipitable water, and it is observed that the degree of drops in the late summer is reduced in the water vapor-corrected time series. Fig. 18 demonstrates the seasonality in the near-infrared bands extracted by the STL method for the Sonoran Desert. The left column shows the curves before the water vapor correction,

This paper has investigated the use of time series algorithms for assessing the long-term degradation of a satellite sensor. Two time series algorithms, STL and DWT, were used to remove the seasonal oscillation in the TOA reflectance time series that can lead to large uncertainty in the assessment of sensor degradation. It has been shown that the two time series algorithms characterized the seasonal oscillation of the time series with different emphases (more regularity in STL and more adaptivity in DWT), but both algorithms succeeded in extracting long-term trends of the TOA reflectance time series without involving physical models such as a bidirectional reflectance model and a radiative transfer model. The trend estimates from the deseasonalized time series have a higher statistical significance than those of the original TOA reflectance and the BRDF-normalized time series while maintaining similar trend estimates. Although the effect of seasonality is relatively small when long-term time series data are available, the effect becomes larger as the whole time period becomes shorter and the variability pattern becomes more irregular. The proposed approach yielded stable trending results for data sets of long


(13 years) and relatively short (7 years) time periods as well as for the data sets that had different levels of data availability. Future work includes the application of this data-driven approach to other airborne/spaceborne sensors that have different sun–target–sensor geometry variations and atmospheric effects, which stem from differences in orbit dynamics and relative spectral function, respectively. An interpolation technique also needs to be modified for different sensors since scenes from, for example, MODIS-type satellites are not necessarily nadir looking with a regular 16-day revisit period, as in Landsat. A correlation analysis between the TOA reflectance variation and other meteorological variables such as the amount of water vapor and aerosol loading is also warranted to further investigate a precise relationship between the time series. R EFERENCES [1] G. Chander and B. Markham, “Revised Landsat-5 TM radiometric calibration procedures and postcalibration dynamic ranges,” IEEE Trans. Geosci. Remote Sens., vol. 41, no. 11, pp. 2674–2677, Nov. 2003. [2] D. Helder and E. Micijevic, “Landsat-5 Thematic Mapper outgassing effects,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 12, pp. 2717– 2729, Dec. 2004. [3] D. Wang, D. Morton, J. Masek, A. Wu, J. Nagol, X. Xiong, R. Levy, E. Vermote, and R. Wolfe, “Impact of sensor degradation on the MODIS NDVI time series,” Remote Sens. Environ., vol. 119, pp. 55–61, Apr. 2012. [4] S. Uprety and C. Cao, “Radiometric and spectral characterization and comparison of the Antarctic Dome C and Sonoran Desert sites for the calibration and validation of visible and near-infrared radiometers,” J. Appl. Remote Sens., vol. 6, no. 1, p. 063541, Jul. 2012. [5] X. Xiong, A. Wu, and B. N. Wenny, “Using Dome C for moderate resolution imaging spectroradiometer calibration stability and consistency,” J. Appl. Remote Sens., vol. 3, no. 1, p. 033520, Mar. 2009. [6] A. Angal, G. Chander, T. Choi, A. Wu, and X. Xiong, “The use of the Sonoran Desert as a pseudo-invariant site for optical sensor crosscalibration and long-term stability monitoring,” in Proc. IEEE IGARSS, Jul. 2010, pp. 1656–1659. [7] A. Angal, X. Xiong, T. Choi, G. Chander, and A. Wu, “Using the Sonoran and Libyan Desert test sites to monitor the temporal stability of reflective solar bands for Landsat 7 Enhanced Thematic Mapper Plus and Terra moderate resolution imaging spectroradiometer sensors,” J. Appl. Remote Sens., vol. 4, no. 1, p. 043525, Jan. 2010. [8] D. Smith, C. Mutlow, and C. Nagaraja Rao, “Calibration monitoring of the visible and near-infrared channels of the Along-Track Scanning Radiometer-2 by use of stable terrestrial sites,” Appl. Opt., vol. 41, no. 3, pp. 515–523, Jan. 2002. [9] A. Angal, G. Chander, X. Xiong, T. Choi, and A. Wu, “Characterization of the Sonoran Desert as a radiometric calibration target for Earth observing sensors,” J. Appl. Remote Sens., vol. 5, no. 1, p. 059502, Jul. 2011. [10] W. Kim, S. Liang, and C. Cao, “On using BRDF models for assessment of radiometric stability of Sonoran Desert,” in Proc. IEEE IGARSS, 2012, pp. 4730–4733. [11] J. E. Pinzbn, M. E. Brown, and C. J. Tucker, “EMD correction of orbital drift artifacts,” in Hilbert-Huang Transformation and its Applications. Singapore: World Scientific Publishing, 2005. [12] R. Cleveland, “STL: A seasonal-trend decomposition procedure based on loess,” J. Official Stat., vol. 6, no. 1, pp. 3–73, Jan. 1990. [13] S. Mallat, “Multifrequency channel decompositions of images and wavelet models,” IEEE Trans. Acoust., Speech Signal Process., vol. 37, no. 12, pp. 2091–2110, Dec. 1989. [14] Z. Li and M. Kafatos, “Interannual variability of vegetation in the United States and its relation to El Nino/Southern Oscillation,” Remote Sens. Environ., vol. 71, no. 3, pp. 239–247, Mar. 2000. [15] S. Sarkar and M. Kafatos, “Interannual variability of vegetation over the Indian sub-continent and its relation to the different meteorological parameters,” Remote Sens. Environ., vol. 90, no. 2, pp. 268–280, Mar. 2004. [16] J. VerBesselt, R. Hyndman, G. Newnham, and D. Culvenor, “Detecting trend and seasonal changes in satellite image time series,” Remote Sens. Environ., vol. 114, no. 1, pp. 106–115, Jan. 2010. [17] D. L. Morstad and D. L. Helder, “Use of pseudo-invariant sites for longterm sensor calibration,” in Proc. IEEE IGARSS, 2008, pp. I-253–I-256.

2975

[18] H. Cosnefroy, M. Leroy, and X. Briottet, “Selection and characterization of Saharan and Arabian desert sites for the calibration of optical satellite sensors,” Remote Sens. Environ., vol. 58, no. 1, pp. 101–114, Oct. 1996. [19] G. Chander, B. L. Markham, and J. A. Barsi, “Revised Landsat-5 Thematic Mapper radiometric calibration,” IEEE Geosci. Remote Sens. Lett., vol. 4, no. 3, pp. 490–494, Jul. 2007. [20] P. M. Teillet, S. Falls, and F. D. Palluconi, “A lifetime radiometric calibration record for the Landsat Thematic Mapper,” in Proc. Can. Symp. Remote Sens., 2001, pp. 17–25. [21] G. Chander, B. L. Markham, and D. L. Helder, “Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors,” Remote Sens. Environ., vol. 113, no. 5, pp. 893–903, May 2009. [22] G. Thuillier, M. Hersé, D. Labs, and T. Foujols, “The solar spectral irradiance from 200 to 2400 nm as measured by the SOLSPEC spectrometer from the ATLAS and EURECA missions,” Solar Phys., vol. 214, no. 1, pp. 1–22, May 2003. [23] R. Irish, “Landsat 7 automatic cloud cover assessment,” in Proc. SPIE, Algorithms Multispectr., Hyperspectr., Ultraspectr. Imagery VI, 2000, pp. 348–355, International Society for Optical Engineering. [24] L. Oreopoulos, M. Wilson, and T. Várnai, “Implementation on Landsat data of a simple cloud-mask algorithm developed for MODIS land bands,” IEEE Geosci. Remote Sens. Lett., vol. 8, no. 4, pp. 597–601, Jul. 2011. [25] C. S. Burrus, R. A. Gopinath, and H. Guo, “Introduction to wavelets and wavelet transforms : A primer,” in Recherche. Englewood Cliffs, NJ, USA: Prentice Hall, 1998. [26] B. Whitcher, P. Guttorp, and D. Percival, “Wavelet analysis of covariance with application to atmospheric time series,” J. Geophys. Res., vol. 105, no. D11, pp. 14 941–14 962, Jun. 2000. [27] W. Staylor and J. Suttles, “Reflection and emission models for deserts derived from Nimbus-7 ERB scanner measurements,” J. Appl. Meteorol., vol. 25, no. 2, pp. 196–202, Feb. 1986. [28] J. Roujean, “A bidirectional reflectance model of the Earth’s surface for the correction of remote sensing data,” J. Geophys. Res., vol. 97, no. D18, pp. 20 455–20 468, Dec. 1992. [29] W. Wanner, X. Li, and A. Strahler, “On the derivation of kernels for kernel-driven models of bidirectional reflectance,” J. Geophys. Res, vol. 100, no. D10, pp. 21 077–21 089, Oct. 1995. [30] M. Kanamitsu, W. Ebisuzaki, J. Woollen, S.-K. Yang, J. J. Hnilo, M. Fiorino, and G. L. Potter, “NCEP-DOE AMIP-II Reanalysis (R-2),” Bull. Amer. Meteorol. Soc., vol. 83, no. 11, pp. 1631–1643, Nov. 2002. [31] H.-Y. Kim and S. Liang, “Development of a hybrid method for estimating land surface shortwave net radiation from MODIS data,” Remote Sens. Environ., vol. 114, no. 11, pp. 2393–2402, Nov. 2010. [32] S. Y. Kotchenova, E. F. Vermote, R. Matarrese, and F. J. Klemm, “Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: Path radiance,” Appl. Opt., vol. 45, no. 26, pp. 6762–6774, Sep. 2006.

Wonkook Kim (S’07–M’11) received the B.S. degree in civil engineering from Seoul National University, Seoul, Korea, in 2004 and the M.S. and Ph.D. degrees also in civil engineering from Purdue University, West Lafayette, IN, USA, in 2005 and 2011, respectively. He served the Korea Army as infantry from 1999 to 2001. He has been with the Department of Geographical Sciences, University of Maryland, College Park, MD, USA, since 2011 and has also been a Visiting Scientist in the National Oceanic and Atmospheric Administration, College Park, MD, USA. His research subjects include the following: 1) vicarious calibration of satellite sensors; 2) characterization of bidirectional reflectance distribution function of calibration sites; 3) robust land cover classification with machine learning algorithms; and 4) hyperspectral data analysis. He is currently working on vicarious calibration of the Advanced Baseline Imager in GOES-R.

2976


Tao He received the B.E. degree in photogrammetry and remote sensing from Wuhan University, Wuhan, China, in 2006 and the Ph.D. degree in geography from the University of Maryland, College Park, MD, USA, in 2012. He has been a Research Associate with the Department of Geographical Sciences, University of Maryland. His areas of interest include surface anisotropy and albedo modeling, data fusion of satellite products, and long-term regional and global surface radiation budget analysis.

Dongdong Wang received the Ph.D. degree in geography from the University of Maryland, College Park, MD, USA, in 2009. He is currently a Research Assistant Professor with the Department of Geographical Sciences, University of Maryland. His primary research interest includes the retrieval of biophysical and flux parameters and spatiotemporal analysis of these variables.

Changyong Cao received the B.S. degree in geography from Peking University, Beijing, China, in 1982 and the Ph.D. degree specializing in remote sensing from the Louisiana State University, Baton Rouge, LA, USA, in 1992. He is a Research Physical Scientist at National Oceanic and Atmospheric Administration (NOAA)/The National Environmental Satellite, Data, and Information Service (NESDIS)/Center for Satellite Applications and Research. He specializes in the calibration of radiometers on board NOAA’s Operational Environmental Satellites and is currently leading the Suomi National Polar-orbiting Partnership (NPP) Visible/Infrared Imager/Radiometer Suite (VIIRS) Sensor Data Records (SDR) team and the GOES-R calibration working group. He is responsible for developing and refining the methodology for intersatellite calibration using the simultaneous-nadir-overpass method, which has been used for the long-term performance monitoring of all radiometers on NOAA’s polar orbiting satellites and is being used by scientists for quantifying intersatellite calibration biases in developing long-term time series for climate change detection studies. Before joined NOAA in 1999, he had five years of aerospace industry experience supporting National Aeronautics and Space Administration (NASA) projects and three years of university teaching experiences.

Shunlin Liang (M’94–F’13) received the Ph.D. degree from Boston University, Boston, MA, USA. He is currently a Professor with the Department of Geographical Sciences, University of Maryland, College Park, MD, USA, and the College of Global Change and Earth System Sciences, Beijing Normal University, Beijing, China. His main research interests focus on the estimation of land surface variables from satellite data, earth energy balance, and assessment of environmental impacts of vegetation changes. He published about 180 peer-reviewed journal papers, authored the book Quantitative Remote Sensing of Land Surfaces (Wiley, 2004), edited the book Advances in Land Remote Sensing: System, Modeling, Inversion and Application (Springer, 2008), and coedited the books Advanced Remote Sensing: Terrestrial Information Extraction and Applications (Academic Press, 2012) and Land Surface Observation, Modeling, Data Assimilation (World Scientific, 2013). Dr. Liang is an Associate Editor of the IEEE T RANSACTIONS ON G EO SCIENCE AND R EMOTE S ENSING and also a guest editor of several remote sensing-related journals.