Change Detection within Remotely Sensed Satellite Image Time Series via Spectral Analysis

: Jump or break detection within a non-stationary time series is a crucial and challenging problem in a broad range of applications including environmental monitoring. Remotely sensed time series are not only non-stationary and unequally spaced (irregularly sampled) but also noisy due to atmospheric effects, such as clouds, haze, and smoke. To address this challenge, a robust method of jump detection is proposed based on the Anti-Leakage Least-Squares Spectral Analysis (ALLSSA) along with an appropriate temporal segmentation. This method, namely, Jumps Upon Spectrum and Trend (JUST), can simultaneously search for trends and statistically signiﬁcant spectral components of each time series segment to identify the potential jumps by considering appropriate weights associated with the time series. JUST is successfully applied to simulated vegetation time series with varying jump location and magnitude, the number of observations, seasonal component, and noises. Using a collection of simulated and real-world vegetation time series in southeastern Australia, it is shown that JUST performs better than Breaks For Additive Seasonal and Trend (BFAST) in identifying jumps within the trend component of time series with various types. Furthermore, JUST is applied to Landsat 8 composites for a forested region in California, U.S., to show its potential in characterizing spatial and temporal changes in a forested landscape. Therefore, JUST is recommended as a robust and alternative change detection method which can consider the observational uncertainties and does not require any interpolations and/or gap ﬁllings.


1.
The detection of jumps (sudden changes) in the trend component of an unequally spaced time series;

2.
Accounting for uncertainties in the time series values (observational uncertainties) to improve the estimation of trend and seasonal components, and jump locations; and 3.
The characterization of the gradual and sudden changes in the ecosystem by estimating the jump direction and magnitude.
The ALLSSA season-trend model, the core of JUST proposed herein, was used for selecting a stable history period for near-real-time monitoring [21]. In this article, the focus is detecting disturbances within the historical time series like BFAST [7,8]. JUST segments the entire historical time series and applies ALLSSA to each segment to find potential jumps. JUST then finds an optimal set of jump locations for the entire time series and finally divides the time series into trend, seasonal, and remainder components like BFAST. JUST is assessed by simulating EVI time series and changing the number of observations, seasonality, noise level, and by adding sudden changes with different magnitudes. JUST is also compared with BFAST using simulated time series with various types and two NDVI time series for a forested region in southeastern Australia. Another reason for choosing BFAST for the comparison herein is that BFAST is a popular change detection method whose advantages and weaknesses are shown by many researchers [35][36][37][38][39], and herein it is also aimed to show the source of some of its weaknesses. The potential of JUST for jump detection and classification is also shown on a collection of Landsat 8 images covering a forested region in the northeast part of Santa Rosa in California, U.S., where BFAST cannot be directly applied as the time series are (inherently) unequally spaced, and BFAST does not consider the observational uncertainties.

Study Regions
The top right panel of Figure 1 illustrates two different study regions chosen herein to test and validate the performance of JUST in detecting and classifying changes within real-world vegetation time series. The first region is a forested area (Pinus radiata plantation) in southeastern Australia, selected herein to compare BFAST with JUST [8]. The second study region is mainly a forested region (≈1600 km 2 ) located in the northeast part of Santa Rosa, the main city in California's famous Wine Country. The second region is covered by a mixture of shrubs and trees, such as oak, madrone, alder, and redwood [40,41]. The ground reference areas, displayed as polygons with semi-transparent fill color toward the north and east of Santa Rosa in the magnified section of Figure 1, are mainly impacted by the massive wildfires in October 2017, known as Tubbs Fire, Nuns Fire, and Atlas Fire [42]. The polygon shapefiles for the burned areas within the second study region, used for validation herein, are obtained from the Monitoring Trends in Burn Severity, conducted through a partnership between the U.S. Geological Survey National Center for Earth Resources Observation and Science and Forest Service Remote Sensing Applications Center, the U.S. Department of Agriculture [43]. Another reason for selecting the second study region is the availability of Landsat 8 satellite imagery every seven and nine days (because of the spatial overlapping of satellite images in higher latitudes), making the image time series (inherently) unequally spaced, see below and [21,44] for more details.

Data Sets and Pre-Processing
Two NDVI time series from 2000 to 2009 are used for the first study region to compare BFAST with JUST. The 16-day Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI composites (250 m spatial resolution) are used to obtain these two time series which are the same as the ones shown in the middle and bottom panels of Figure 8 in [8], downloaded from [45]. The binary MODIS Quality Assurance flags are used to choose the cloud-free data with optimal quality. The NDVI time series corresponding to the bottom panel of Figure 8 in [8] has several missing values, replaced by linear interpolation between neighboring values within the time series before applying BFAST [8,45]. However, no interpolation technique is required when applying JUST. Study regions: the first region is a forested region in New South Wales, Australia, the same region studied in [8]. The second region is a forested region in California, U.S. The magnified section shows a view of the second study region in January 2020 by Google Earth. The transparent areas in the magnified section are generated by adding a layer containing the polygon shapefiles for the burned areas in October 2017 obtained from [43]. The topographic maps are generated by the SimpleMappr online tool [46].
For the second study region, Landsat 8 surface reflectance data products from the U.S. Geological Survey (USGS) with 30 m spatial resolution are used [47]. The Landsat 8 satellite is in a near-polar, sun-synchronous, circular orbit with 705 km altitude and has a 16-day repeat cycle [48]. The Landsat orbits are poleward convergent, and thus Landsat 8 sensor can acquire more frequent scene coverage at higher latitudes [49]. The amount of Landsat data in the U.S. Landsat archive is not constant temporally and spatially among Landsat sensors due to various factors, such as sensor health, and the capability of data reception and acquisition [48][49][50]. According to the product guide for the Landsat 8 surface reflectance code, the calibration parameters from the metadata files associated with the Landsat 8 data are used to generate Top Of Atmosphere (TOA) reflectance and TOA brightness temperature. Then an atmospheric correction procedure is applied to Landsat 8 TOA reflectance data, using additional input data, such as water vapor, ozone, and aerosol optical thickness, as well as a digital elevation derived model to generate the surface reflectance data products.
The surface reflectance bands 2 (Blue), 4 (Red), and 5 (Near Infrared) are used to calculate EVI defined as follows where ρ NIR , ρ RED , and ρ BLUE signify the spectral reflectance measurements acquired in the Near Infrared, Red, and Blue regions, respectively. Coefficients G, L, C 1 , and C 2 are the gain factor, canopy background adjustment, and aerosol resistance terms, respectively [10]. Herein, the adopted coefficients G = 2.5, L = 1, C 1 = 6, and C 2 = 7.5 for the Landsat 8 imagery are used [51]. The scale factor 0.0001 is multiplied by each band value to convert their values to floating point values from 0.0 to 1.0 before calculating EVI. This scale multiplication is not necessary for NDVI because NDVI defined by with the value range from −1 to 1, does not have the adjustment factor L [11,34,52]. The value range for EVI is usually from −1 to 1 (may exceed) and varies between 0.2 and 0.8 for healthy vegetation [52].

Weighted Vegetation Time Series
A per-pixel vegetation time series is obtained from a set of spatially aligned pixels in the vegetation images, also referred to as vegetation time series [21]. Figure 2 demonstrates the process of obtaining an EVI time series from a set of spatially aligned pixels (shown in solid colors). The Pixel Quality Assessment (PQA) band obtained for each Landsat 8 satellite image is used to define the weights. According to the product guide for Landsat 8, if a pixel value in the PQA band is 322, 386, 834, 898, or 1346, then the pixel is clear, otherwise, there is a high chance that the pixel is contaminated with clouds, cloud shadows, water, snow, ice, smoke, or haze. The following technique shows how one may obtain a weighted vegetation time series using the PQA band. An illustration of obtaining per-pixel weighted EVI time series. EVI images are acquired over unequally spaced time intervals due to availability of valid satellite imagery that provide unequally spaced EVI image time series. In this illustration, each EVI image consists of 100 pixels, and so in total there are 100 unequally spaced time series whose sizes may not be the same due to several reasons, such as cloud contamination and high aerosol content. The cloud score for each spatial window of size nine is used to define an appropriate weight for the pixel located at the window center (see the shaded gray pixels around the pixel in a solid color for each image).
The key point here is to assign lower weights to the observations with higher uncertainties, and vice versa. For a given Landsat 8 EVI image, first count the number of pixels for each spatial window of size nine whose PQA band values are not equal to 322, 386, 834, 898, or 1346 and denote it by Λ. To account more for the observational uncertainty corresponding to the pixel located at the window center, add one unit to Λ if that pixel is contaminated, so 0 ≤ Λ ≤ 10. If Λ = 10, then there is a high chance that the pixel in the center of the spatial window is contaminated, and so one may mask out that pixel (ignore the observation). Otherwise, assign an appropriate weight to each pixel located at the window center, herein defined by 1/ 1 + Λ 2 and multiplied by 100 to express it as percentage. Therefore, for all spatially aligned windows with the same region coverage, a weighted vegetation time series will be obtained for the spatially aligned pixels located at centers of the spatial windows, see the shaded gray and solid pixels in Figure 2 used to define weights for the solid pixels. Continuing this process for all the spatially aligned pixels located at the centers of the spatially aligned windows, a set of weighted vegetation time series will be obtained covering a region of interest (the second study region herein). For each pixel, considering the uncertainties in its neighboring pixels can reduce the noise effect of the pixel during the season-trend analysis because it is possible that the PQA band incorrectly states that the pixel is clear (cloud-free, etc.) while the PQA band shows some of the neighboring pixels are in fact contaminated. Choosing larger sizes for the spatial window (e.g., 25 pixels) may also provide a more reliable weight for the pixel in the center of the spatial window. The weighting procedure presented here could be applied to image time series prior to any change detection method. Therefore, this is an option for a generally improved pre-processing of optical satellite data.

Jumps Upon Spectrum and Trend (JUST)
This section focuses on analyzing a time series that can be a weighted vegetation time series. JUST is established by simultaneous estimation of spectrum and trend for a given time series with an associated covariance matrix whose inverse, denoted by P (weight matrix), can be used to control the effect of noisy measurements [30]. Please note that the relative weights play a significant role for improving the estimation of trend and seasonal components, but the proposed method is not sensitive to any covariance factor [53]. In the next section, the weight matrix is a diagonal matrix whose diagonal entries are obtained from the observational uncertainties and is treated as a vector for the sake of computational efficiency. The seasonal component of the time series may not be stable over time, i.e., its frequencies and amplitudes may change over time, and so performing the spectral analysis on the entire time series can potentially reduce the accuracy of estimating the trend and seasonal components [30]. To detect jumps, a temporal segmentation technique is proposed to account for changes within the trend and seasonal components. A sequential approach is then performed that scans each time series segment using ALLSSA to determine statistically at which location there exists a significant jump in the segment.
For a better understanding of the segmentation of the time series, the term window (temporal window) is used. In this terminology, a window contains a time series segment with its associated statistical weights as well as the constituents of known forms (e.g., linear trend) and sinusoidal functions being fitted to the segment [30]. The window size R is defined as the number of observations within the window or simply the segment size. However, the window length is the time difference between the last and the first observations within the window [30]. For an unequally spaced time series, a translating window with a fixed size rather than a fixed length may contain additional observations from the neighboring seasons which may allow better estimation of trend and seasonal components, a useful approach when there are few observations available within a particular season [54]. Please note that the window size selected herein is fixed during translation by δ observations over time, and it is independent of frequency, a special case of LSWA when the window size parameter L 1 is zero [30].
For each segment of size R within the translating window, ALLSSA can be set up to find an optimal piece-wise linear function with two pieces along with a set of estimated sinusoidal functions (the sequential approach). In other words, the constituents of known forms in ALLSSA may be selected as the linear functions with two pieces whose intercepts and slopes will be simultaneously estimated with the sinusoidal functions. Therefore, sequentially two linear pieces will be selected that along with the estimated sinusoids minimize the weighted L2 norm of the residual segment, where the residual segment is the original segment minus the estimated linear trend with two pieces and the estimated sinusoids. In ALLSSA, for a given linear trend with two pieces as the constituents of known forms, the frequencies and amplitudes of the sinusoidal functions will be estimated simultaneously with the intercepts and slopes of the two linear pieces until there is no more statistically significant peak in the spectrum at a certain confidence level (usually selected at 99% which corresponds to the significance level α = 0.01). Now, letâ 1 andâ 2 be the estimated slopes andb 1 andb 2 be the estimated intercepts of the two linear pieces for the segment, respectively. These estimated coefficients can be used to derive the direction and magnitude of the potential jump respectively as where the hat symbol stands for estimation, and t is the time of the potential jump determined by the sequential approach. Please note that JUST can also detect t as the time when a jump in trend occurs even if MAG = 0 because the intercepts and slopes of the two linear trends can be different (DIR = 0) but still making the magnitude equal to zero (the word "jump" is used loosely in this case). Figure 3a demonstrates three time series such that each has a jump at t 9 (the ninth observation). The time series whose observations are shown by the blue circles has a weak jump at t 9 with a small negative direction and a small positive magnitude. For example, a very small and negligible part of an arid region, corresponding to a per-pixel NDVI time series, starts growing plants around the jump location. One may disregard this jump by applying certain thresholds (e.g., |DIR| < 0.01 and |MAG| < 0.05, where | · | is the absolute value). The time series, shown by red diamonds, has a jump with MAG = −0.2. For example, in real-world applications, a jump with a negative magnitude may have been created by harvesting, wildfire, insect attack, etc. The time series, shown by green squares, has a jump with positive direction and zero magnitude, e.g., it may have been created indirectly by irrigation after a gradual drought. DIR and MAG denote the true direction and magnitude of the jump located at t 9 , respectively, and (b) the translating window of size R = 21 observations, translating (sliding) by δ = 6 observations. JUST has detected a jump within the first window (see the red arrow) and another jump within the second window (see the green arrow). The solid red and green lines are the estimated linear trends (each trend has two pieces) corresponding to the segments within the red and green windows, respectively.
JUST evaluates a determined jump location t (breakpoint) by its occurrence that is the number of times that t is determined as a jump location within all the windows containing t . When there are several identical or close enough jump locations, each detected within a window, JUST selects the one that has the maximum occurrence, then the one whose time index is closer to the time index of its corresponding window location [30]. The latter is because the selected t has data support from both sides of the window, making the magnitude and direction of the jump more accurate [30,31]. Certain thresholds for the direction and magnitude of a jump may be applied at this point to ignore weak changes and to reduce false-positive rates in real-world applications.
Given the jump locations, ALLSSA may optionally be applied to a new set of segments to estimate the trend and seasonal components within each segment. A weighted moving average method, whose weights are obtained from a symmetric Gaussian function [30], may then be applied to estimate the trend and seasonal components of the entire time series. Alternatively, given the jump locations, the LSWAVE software may be implemented to study seasonal changes in the time-frequency domain [32].
To understand the JUST algorithm better, an NDVI time series along with two windows are illustrated in Figure 3b. The red window translates to the green window by δ = 6 observations (the windowing or segmentation approach). Within each window, the sequential approach is performed which finds a potential jump in the segment, see the red and green arrows pointing to the jumps detected within the red and green windows, respectively. Since the two detected jump locations are close to each other, JUST selects the point in the segment within the green window (see the green arrow) because it has more observations on its right hand side within the green window to support a more reliable simultaneous estimation of trend and seasonal components. Having the list of all detected jump locations within translating windows is useful as it may also reveal the potential outliers. For example, the observation shown by the red arrow in Figure 3b, may be an outlier due to the presence of atmospheric or background noise that could not be detected from a cloud masking method, neither an appropriate weight could be assigned to this observation to mitigate its effect [21,55]. Please note that within each window, ALLSSA estimates a separate set of sinusoidal functions that can change as the window translates to account for the phenological changes [7,8,56]. Figure 4 shows the flowchart of JUST for the sake of completeness.
The Ordinary Least-Squares (OLS) estimation simply fits the constituents of known forms, such as trends and sinusoids of known frequencies to the time series to estimate their coefficients [22]. For each time series segment, the OLS-based method proposed herein (the simple case of JUST) finds an optimal piece-wise linear function with two pieces along with sinusoidal functions of fixed cyclic frequencies (1, 2, 3, 4 cycles per year) like in BFAST. The flowchart of the OLS-based method for a given time series is the same as the one shown in Figure 4 but instead of applying ALLSSA, OLS is applied to estimate the trend and sinusoids of fixed cyclic frequencies (1, 2, 3, 4 cycles per year) for each segment. The OLS estimation not only may increase the false-positive rate for change detection but it may also cause the over/under-fitting problem [21,57]. Unlike the OLS-based method, JUST searches each segment to find statistically significant spectral peaks without any assumptions on the existence of seasonal component and its frequencies. JUST simultaneously estimates the trend and statistically significant seasonal components within each segment, introducing a robust change detection method.

JUST
Time series: , · · · , w(t n ) Window size: R; Translation step: δ Initial set of cyclic frequencies: Ω Significance level: α Apply ALLSSA to simultaneously estimate the trend Γ and seasonal S components and obtain the residual segment r : The time index of the potential jump in the segment:  Figure 4. Flowchart of JUST. The weights are optional inputs and in general can also be entered in the form of a matrix P, the inverse of the covariance matrix associated with the time series. The constituents of known forms are the linear trends by default (Γ ). Symbol | · | denotes the absolute value. The flowchart of ALLSSA used in JUST is in Figure 3.4 of [22], or see [31].
To validate the sequential approach and also to seek an effective window size for the segmentation approach, in Sections 3.1.1 and 3.1.2, the results of jump detection methods using two window sizes R = 2M and R = 3M are compared, where M is an average sampling rate. For an unequally spaced time series, the sampling rate cannot be a fixed number [58], and so an average sampling rate may be computed as M = n/(t n − t 1 ), where n is the number of observations, t 1 is the start time, and t n is the end time, e.g., M = 23 observations per year, where t 23 = 1 and t 1 = 0. Both window sizes can account for inter-annual variations caused by several factors, such as climate change, land management or land degradation [7,59]. It is shown that selecting window size R = 3M performs better in identifying true jumps within time series with poorer signal-to-noise ratio than R = 2M, however, R = 2M can be a better choice for less noisy time series whose components have more variability of amplitude and frequency over time. Please note that noise in this context refers to components other than trend and seasonal, such as random noise from cloud contamination, and other high-frequency signatures. For the automation purposes in environmental applications, R = 3M and δ = M are recommended. This selection is capable of detecting jumps occurred approximately one year apart and is robust in capturing inter-annual and seasonal variations within satellite image time series. The main inputs and outputs of the JUST code are listed in Table 1. This section briefly describes BFAST due to its similarities to JUST and the numerical comparisons made in the next section. BFAST is an iterative algorithm that decomposes a time series into trend, seasonal, and remainder components [7]. The trend component between any two consecutive jump locations is assumed to be linear and the seasonal component is estimated by summing the cosine and sine (sinusoidal) functions of fixed frequencies (e.g., 1, 2, 3 cycles per year). The OLS residual segment is obtained by subtracting the estimated trend and/or sinusoidal functions of fixed frequencies via the Least-Squares fit from the segment [22,60].

Validation through Descriptive Statistics
Suppose that there are K time series where each time series has a known jump at a specific location, and a change detection method is set to detect a maximum of one jump in each time series. The following statistical measures are used for validation in the next section: 1.
Jump error for the K time series: the number of incorrectly detected jump locations divided by K, i.e., normalized.

2.
Root Mean Square Error (RMSE) for the jump magnitude when the jump location is correctly detected: where m is the number of correctly detected jump locations (m ≤ K), MAG k and MAG k are the true and estimated magnitudes of a correctly detected jump, respectively, see Equation (4) where f is a time series with its associated weights w (if available), n is the number of observations, and r is the residual series obtained after subtracting the estimated trend and seasonal components from f .

Results and Discussion
JUST can be applied to any non-stationary time series, and it is not limited to satellite image time series. In JUST, a time series can be equally or unequally spaced and does not require any interpolations or gap fillings, neither it manipulates the observations (i.e., JUST does not repair the data and analyzes them as they are provided). JUST determines whether there are any statistically significant components in the time series, and it accurately estimates the components by considering the observational uncertainties. In this section, JUST and its simple case (the OLS-based method) are validated using simulated and real-world vegetation time series with varying trend and seasonal components, and noises. Furthermore, JUST is compared with BFAST using both simulated and real-world vegetation time series.

Simulation Experiment
Simulation of satellite image time series is very challenging due to its immense dependence on various factors, such as vegetation phenology, inter-annual climate variations, sensor conditions (e.g., viewing angle), and signal contamination (e.g., clouds) [20]. Therefore, if there are any periodic signals in the time series, their cyclic frequencies can be any real numbers, and so it is necessary to estimate them accurately using statistical inference in order to detect jumps rigorously. Below, the advantages of JUST over the OLS-based method and BFAST are shown by simulating noisy time series with unknown seasonality, i.e., when the seasonal component of the time series cannot be well estimated by the sum of sinusoidal functions whose cyclic frequencies are integers.

Simulation of Time Series with Unknown Seasonality
To show the advantages of JUST over the OLS-based method, an EVI time series is simulated by summing a seasonal component S: and a piece-wise liner trend function Γ with two pieces: To represent cloud contamination that may remain after atmospheric and cloud masking procedures, a noise component η 1 is added using a random number generator following a normal distribution with mean zero and standard deviation σ [7]. Thus, f (t j ) = S(t j ) + Γ (t j ) + η 1 (t j ), 1 ≤ j ≤ n. Hereafter, the noise level is expressed as 4σ, which is approximately 99% of the noise range to allow one to compare the noise level with magnitude.
Since in the Landsat 8 data example, the images are acquired every nine and seven days consecutively, in the simulated examples, the distribution of the t j are set up accordingly unless stated otherwise. In other words, t j follow the sequence of (0, 9, 16, 25, 32, . . . ) divided by 365.25. Note that JUST does not depend on the distribution of t j , and it can be directly applied to any equally or unequally spaced time series without the need for any interpolations.
For some given sets of amplitudes A i from 0 to 0.5, phases θ i from 0 to 2π, and cyclic frequencies ω i from 0.5 to 4 (i = 1, 2, 3), 1000 time series are generated for each noise level in {0, 0.008, 0.016, . . . , 0.24}, while randomly removing up to 50% of observations from each time series and changing the jump location. The locations of true jumps in all the simulated time series are recorded to check them against the estimated ones (model sensitivity test). For each noise level, the OLS-based method and JUST are applied to the same 1000 time series to estimate the jump locations, magnitudes, and residual series. Therefore, for each noise level and for each change detection method, the jump error, RMSE, and MNRN for the 1000 time series are also computed, see Section 2.6 where K = 1000 and the weights w are not considered. Figure 5 shows the results of analyzing the 1000 time series for each noise level when the parameters in Equations (7) and (8) are given by A 1 = 0.1, A 2 = A 3 = 0.05, θ 1 = π/4, θ 2 = π/3, θ 3 = 0, ω 1 = 1.1, ω 2 = 2.2, and ω 3 = 3.3. Each panel in Figure 5 shows four different plots for each method, obtained for four different magnitudes listed in Table 2 Figure 5c shows how MNRN in JUST is less than MNRN in the OLS-based method, signifying that the accuracy of estimated trend and seasonal components using OLS is generally poorer than ALLSSA.
The same analyses with the same values for the parameters mentioned above are performed on a longer time series (a three-year period), and the results are illustrated in Figure 6. The results show that when the noise level is less than the magnitude, the jump error in JUST is less than 0.1. Comparing Figure 6 with Figure 5, one can clearly observe that the jump error and RMSE of magnitudes decrease significantly when three years of data are used. The increase in MNRN values using OLS for the three years case is mainly due to the additional errors coming from the additional observations. Therefore, when a time series has stable seasonal component over time and three years of data are used to estimate the jump locations and magnitudes, the results are more resistant to poorer signal-to-noise ratio. However, this conclusion is more valid when the components are almost stable over time, and the time series does not have multiple jumps within the specified time period. JUST is also tested for other values of amplitudes, phases, frequencies, and trend coefficients, and similar results are obtained, but not shown here for brevity. When the seasonal component can be well estimated by sinusoidal functions whose cyclic frequencies are integers, the OLS-based method and JUST produced very similar results (e.g., when setting ω 1 = 1, ω 2 = 2, and ω 3 = 3). However, in practice, the cyclic frequencies can be any real numbers.
To compare JUST with BFAST, since BFAST R-code [45] does not support unequally spaced time series, for each noise level, 1000 equally spaced time series are simulated for a three-year period (16-day, about 23 observations per year) with the same set of coefficients for the sinusoidal functions (except A 3 = 0 in Equation (7)) and with the trend coefficients listed in Table 3. Table 3. Parameter values for simulation of the trend in Equation (8) to compare JUST with BFAST.  Table 4. As the results in this table show, JUST performed much better than BFAST in detecting the true jump and estimating its magnitude particularly when MAG = −0.1. The higher the magnitude gets, the closer the results of BFAST become to the ones in JUST because the seasonal component is less dominant at higher magnitudes. In fact, since BFAST is designed according to OLS-MOSUM with fixed frequencies [7], it cannot rigorously estimate the true seasonal component, causing significant biases into its model, and thus performs poorly in cases when the seasonal component is unknown (e.g., waves whose frequencies are any real numbers and can change over time). The jump error in BFAST appears to slightly decrease with noise for MAG = 0.0 (with DIR = a 2 − a 1 = 0.19 = 0) and MAG = −0.1. This is likely due to the accidental detection of jumps at higher noise levels as the unknown seasonal component possibly becomes less dominant. The jump error in BFAST increases with noise when MAG = −0.2 as the effect of seasonal component reduces at higher magnitudes and the increased noise level decreases the accurate detection of jumps. However, the jump error in BFAST is significantly higher than JUST in this example. Table 4. Comparison between JUST and BFAST for the time series simulated by the seasonal component given in Equation (7) with trend coefficients given in Table 3 and added noise. N/A means no answer. Please note that in this example the jump (break) with zero magnitude has positive direction.

Simulation of Time Series with Two Noises of the Same Type
To demonstrate another advantage of JUST, the performance of JUST is examined when considering the weights. To do so, in addition to η 1 , another noise component η 2 is added to the time series as described below. To show the possibility of the existence of other types of seasonal component in certain land covers, the following asymmetric Gaussian function G is used: where A is the amplitude, τ is the position of the maximum with respect to the independent time variable t, and σ 1 and σ 2 are positive numbers that determine the width of the left and right hand sides, respectively [7]. Herein, only a case is shown when A = 0.2, and S is defined as where | · | is the absolute value and Γ is given by Equation (8). The noise level for η 1 is set to 0.1 to account for noise that may remain after the cloud masking and atmospheric correction procedures, see Figure 7. Since in forested and vegetated lands, the EVI values are positive, cloud contamination will typically decrease the EVI values, the reason why the absolute value of η 2 is subtracted from the simulated EVI time series [21]. For each η 2 noise level in {0, 0.02, 0.04, . . . , 0.4}, 1000 time series are generated while randomly removing up to 50% of observations from each time series and changing the jump locations. The true jump locations are recorded to test the model sensitivity. Noise η 2 is defined after the cloud masking procedure that can be partially controlled by introducing appropriate weights, e.g., 1/ 0.01 + η 2 2 (t j ) , 1 ≤ j ≤ n, so that the noisier measurements get lower weights during the analysis [21,30].  Table 2, (c) the noise η 1 with 0.1 level and noise η 2 with 0.4 level, (d) the weights defined by η 2 , and (e) the simulated EVI time series that is the sum of the seasonal, trend, and η 1 components minus the absolute value of η 2 . To help comparisons, the solid bars with the same data range are displayed on the right-hand side of the plot.
JUST is applied with and without considering the weights to estimate the jump locations, magnitudes of change, and residual series. Figure 8 shows the results of JUST with (solid lines) and without (dashed lines) considering the weights for the four different magnitudes listed in Table 2. From Figure 8, one can clearly observe that considering the weights in JUST (weighted JUST) reduces the jump error, RMSE of magnitudes and MNRN, and the differences between the results of weighted and unweighted JUST increase gradually as the η 2 noise level increases.
while other components have the same structures as before. Figure 9 shows the new results, signifying that using three years of data in this example can also decrease the jump error and RMSE of magnitudes significantly for noisy observations. JUST has been successfully tested for some other values for A, τ, σ 1 , σ 2 , and trend coefficients, and similar results are obtained, not illustrated here. Generally, for time series with almost stable seasonality over time, using three years of data improves the accuracy of jump locations and magnitudes, and it is more resistant to the poorer signal-to-noise ratio (having higher true positive rate). This accuracy will improve when the weights associated with the time series values are also considered in JUST.
To compare BFAST with the unweighted JUST when the seasonal component is defined by Equation (11), for each noise level η 1 , 1000 equally spaced time series are simulated (16-day, for a three-year period) with the trend coefficients listed in Table 3. The results of both methods on the same set of time series are shown in Table 5. Comparing Tables 4 and 5, the BFAST results are improved in this case because the asymmetric Gaussian components can be estimated better using the harmonic functions with fixed frequencies. However, JUST still performed better than BFAST in most cases. Table 5. Comparison between JUST and BFAST for the time series simulated by the asymmetric Gaussian function in Equation (11) with trend coefficients given in Table 3, and added noise. Please note that in this example the jump (break) with zero magnitude has positive direction.

A Simulated EVI Time Series with Multiple Jumps
To show how JUST can detect multiple jumps in the trend component of an inherently unequally spaced time series, an EVI time series is simulated that has one observation every nine and seven days since 2013 with 60% missing values. The seasonal component of this time series is defined using Equation (7), where A 1 = 0.1, A 2 = 0.05, A 3 = 0, ω 1 = 1, ω 2 = 2.3, θ 1 = π/2, and θ 2 = π/4. Two jumps are introduced in years 2016 and 2018 with magnitudes 0.23 and −0.1 (EVI), and directions 0.15 and 0.025 (EVI per year), respectively. The seasonal, trend, and noise η 1 (0.1 noise level) components are added to generate the time series. Then the absolute value of noise η 2 (0.4 noise level) is subtracted from half of randomly selected time series values, see Figure 10a. The JUST decomposition results are shown in Figure 10, signifying that the proposed approach can rigorously considers the weights to accurately detect the jump locations and estimate the trend and seasonal components without any need for interpolation.
As described earlier, noise η 2 is used to define the weights. As it can be seen in Figure 10b, half of the observations are associated with weights equal to 100 because they are assumed clean observations (cloud-free, etc.), although they are still contaminated by noise η 1 (the noise that often remains after the atmospheric correction and cloud detection procedures). The other half of observations are associated with weights less than 100 according to noise η 2 . In real-world applications, these weights may be defined from a cloud detection algorithm (see Section 2.3). Many researchers mask the noisy observations using a cloud score threshold (e.g., remove the observations whose corresponding cloud scores are greater than 30%). However, removing such observations may result in getting a very coarsely sampled time series. JUST has this feature to keep these observations, but instead it considers their uncertainties because these observations may still be clean and valid.

Detecting and Characterizing Jumps in the NDVI Time Series for the First Study Region
For the first NDVI time series shown in Figure 11, one can observe that JUST and BFAST with their default parameters have detected four jumps. JUST has accurately detected a jump corresponding to the harvesting operation activity in August 2004, however, BFAST detected a jump in September 2004. When the number of jump locations is set to two, BFAST detects a jump in October 2004 that is more than one month ahead of the harvesting period. Furthermore, in Figure 11c, BFAST shows an unusual seasonal anomaly in the early 2005 which disappears when the number of jump locations is set to two, signifying the model sensitivity to the input parameters. BFAST has detected another jump in December 2005 with a small magnitude but large direction while JUST has detected a jump in May 2006. Setting two jumps as the input, BFAST detects a jump in February 2006, see the middle panel in Figure 8 in [8] for the demonstration.
For the second NDVI time series shown in Figure 12, JUST with its default parameters has detected four jumps while BFAST with h = 0.13 has detected three jumps, where h is the minimal segment size between potentially detected jumps in the trend model given as fraction relative to the sample size [45]. BFAST with its default parameter (h = 0.15) detects three jumps that are in late 2002, early 2004, and late 2006 as shown in the bottom panel in Figure 8 in [8]. In 2002, the region experienced a severe drought which stressed on the pine plantations, and so it resulted in decreasing the NDVI values. The jump detected in late 2006 is due to severe tree mortality as trees were drought-stressed and could not protect themselves against insect attacks [8]. The linear interpolation used before applying BFAST introduced some biases in the BFAST model which may have altered the jump detection results [39,62]. The slight curvature in the trend component in JUST from 2003 to 2007 as seen in Figure 12b is due to the smoothing operation, i.e., when the Gaussian moving average is used for all the estimated linear trends within the translating windows. This example shows how JUST and BFAST can detect and characterize jumps related to forest health.

Detecting and Characterizing Jumps in the Landsat 8 Image Time Series for the Second Study Region
JUST is only applied to all the weighted EVI time series (30 m spatial resolution) whose weighted averages are above 0.2 EVI, see Figure 13a. Figure 13b demonstrates the year of the major jump, i.e., when JUST has detected a jump whose absolute value of magnitude is the largest among the magnitudes of other detected jumps within the EVI time series. The magnitude and direction maps of the major jumps are also illustrated in Figure 13c,d, respectively. In environmental applications, it is interesting to see whether a jump corresponds to a sudden increase or decrease in the signal (the sign of the magnitude). Furthermore, the direction of jump, as expressed here, shows the change in the slope of two linear trends. The per-pixel EVI time series inside the polygon shapefiles shown in Figure 1 are analyzed by JUST and the OLS-based method. JUST and the OLS-based method could successfully detect the locations of the jumps with negative magnitudes, caused by wildfires in October 2017, for approximately 94% and 86% of these time series within one-month accuracy, respectively. To see more details of the analyses and compare JUST with the OLS-based method, three locations (pixels) on the map are selected, labeled by A, B, and C in Figure 13c. Locations A and B are known to be impacted by the Nuns Fire and Tubbs Fire, respectively. The decomposition results for time series corresponding to pixels A, B, and C are illustrated in Figures 14-16, respectively. The time series have gaps mainly during the winter, but neither of JUST and the OLS-based method requires gap-filling. From Figure 14c, one can observe that JUST and the OLS-based method have truly detected the jump in October 2017 with MAG = −0.06, explaining why the magnitude of the sudden change is negative and the slope of the gradual change is positive afterward (recovery). Figure 15 shows an example when the OLS-based method could not detect the jump with a small magnitude caused by the wildfire in October 2017. Furthermore, this example shows why the absolute value of η 2 noise was subtracted from the simulated time series, see the downward spikes in the time series whose weights are close to zero, panels (a) and (b) in Figure 15. Please note that the cloud contamination may not always reduce the EVI values [55]. Finally, Figure 16 shows how both methods have detected identical jumps, except for the jump in 2016, in the time series corresponding to pixel C in the northeast part of the second study region. The overall estimated trend and seasonal components using both JUST and the OLS-based method in these three time series are approximately the same, although JUST has detected the jump locations more accurately.
As mentioned above, Figure 13 shows the JUST results for only the jumps with the largest absolute values of magnitudes, e.g., out of the few jumps detected in the time series corresponding to pixels A, B, and C (see blue lines in Figures 14c, 15c and 16c), only the details of the jumps in 2017 for A and B and the jump in 2014 for C are illustrated in Figure 13b-d. Overall, the results illustrated in Figure 13 not only show the robustness of JUST in detecting and characterizing changes but also emphasize the significance of simulating time series in a controlled environment because it is not easy to find validation data for all land cover changes.
The computational complexity of JUST depends on several factors, including hardware, implementation, the covariance matrix (fully populated or diagonal), and the number of statistically significant spectral peaks that are generally unknown. A fully populated covariance matrix contains the variances of observations in its diagonal and the covariances between the observations in its off-diagonal. When the observations are statistically independent of each other (zero covariances), then the diagonal entries of the covariance matrix can be treated as an array for the sake of computational efficiency. In the simulated results shown in Tables 4 and 5, the processing of 1000 time series of size 69 (three-year period) using JUST on a normal desktop computer took approximately 75 s in MATLAB and 460 seconds in Python, while the processing of the same time series using BFAST took 225 s in R [63]. The OLS-based method proposed herein generally performs faster than JUST (about ten times or more). The MATLAB and Python codes for JUST and the OLS-based method are available and will appear freely online soon [64].

Conclusions and Future Work
In this paper, a robust jump detection method is proposed for non-stationary satellite image time series, namely, JUST (with its simple case, the OLS-based method). JUST is successfully applied to estimate the jump locations and magnitudes in many unequally spaced and weighted time series. The results showed that the jump detection depends on the signal-to-noise ratio, but JUST is more resistant to poorer signal-to-noise ratio than the change detection using the OLS-based method and BFAST, especially when the seasonal component cannot be well estimated by sinusoidal functions of integer cyclic frequencies (see Figures 5 and 6 and Tables 4 and 5). It was shown how one may use the uncertainties associated with likely contaminated pixels (e.g., cloudy, smoky, hazy) to define appropriate weights. It was demonstrated in Figures 8-10, 14 and 15 how JUST can take these weights into account to mitigate the effect of noisy measurements and can consequently improve the estimation of the jump locations and components in the time series.
JUST can rigorously analyze time series corresponding to any geographical location without any assumption on the existence and type of seasonal component. JUST is not limited to any particular time series, and it can be applied to any type of time series without the need for interpolation or editing the original measurements. JUST uses the original data to determine whether there are any statistically significant spectral peaks in the least-squares spectrum. JUST is able to accurately estimate the frequencies and amplitudes of sinusoidal functions which results in a more accurate estimation of jump locations and magnitudes. Section 3.1.1 discussed cases when the seasonal component cannot be estimated well using the sinusoidal functions of fixed cyclic frequencies (e.g., 1, 2, 3 cycles per year) assumed in the OLS-based method or BFAST. The poor estimation of the seasonal component results in propagating error in the BFAST model when the time series is subtracted from its estimated seasonal to search for jumps in the trend component, failing to accurately estimate the jump locations. Postulating a finite set of possible frequencies in a remotely sensed time series, as in the case of the OLS-based method and BFAST, is however a fair assumption in majority of cases. It is worth mentioning that the statistical testing in JUST is not limited to the harmonic functions and linear trend, and other types of constituents, such as quadratic trend, cubic trend, and wavelets, can be integrated into the model. Another advantage of JUST over BFAST is that JUST simultaneously estimates the trend and statistically significant seasonal components for each time series segment, and it does not simply remove the seasonal or trend components from the time series without considering the effect of removing one or another in the residual time series as done in BFAST, see Appendices A and B in [31]. This means that the correlation between the trend and seasonal components is not carried forward in BFAST, making it less accurate and less compatible with unequally spaced time series or even equally spaced time series with unknown signal frequencies. Therefore, we recommend JUST as an alternative method for detecting and classifying jumps within non-stationary satellite image time series, knowing the fact that BFAST with its recent modifications (BFAST monitor) is considered as a powerful change detection method in a variety of ecosystems and applications [65]. JUST promises nice features and great opportunities, but like any other method, its usability in real-world applications still has to be further investigated. Such investigation is strongly recommended as it will provide the opportunity to improve the methods and learn more about the ecosystems.
In JUST, the temporal window size is fixed during translation, and two effective and reasonable sizes (double and triple the sampling rate) are investigated herein. A common shortcoming of all the three change detection methods (JUST, the OLS-based method, and BFAST) is their sensitivity to the segment size. In certain applications, when there is significant variability of frequency and amplitude within the seasonal component over time, the window size may be defined to have variable sizes for different frequencies to account for all the irregularities like in LSWA [22,30,32]. Further research and development need to be done to adapt JUST to the spectrogram rather than the spectrum. The proposed windowing approach may also be used to find the most recent significant jump to determine a stable historical period in a time series. Then, the stable historical series can be used to forecast the series from which a disturbance can be detected in near-real-time [21,62,65].

Acknowledgments:
The authors acknowledge the Landsat mission scientists, and associated USGS and NASA personnel for the production of the data used in this research. The authors would like to thank Spiros Pagiatakis, Tony Ware, Babak Farjad, David Bekaert, Zhe Zhu, and the anonymous reviewers for their time and comments that significantly improved the presentation of the paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this article: