The Potential of the Least-Squares Spectral and Cross-Wavelet Analyses for Near-Real-Time Disturbance Detection within Unequally Spaced Satellite Image Time Series

Abstract

Near-real-time disturbance detection within the remotely sensed time series has become a crucial task in many environmental applications that can help policymakers and responsible authorities to make rapid decisions and proper actions. Although there are several techniques for the near-real-time monitoring of time series, their reliability in regions with higher latitudes are not yet assessed, particularly in regions with consistent data gaps in certain time periods and with large observational uncertainties. A new method is proposed that determines a stable history period from which the least-squares spectral analysis can detect and classify the changes in newly acquired data. To validate the effectiveness of the method, both simulated and real-world vegetation time series obtained for a region in northern Alberta, Canada, are used, where there are consistent data gaps from November to April each year due to the availability of valid Landsat satellite imagery and climate conditions. Furthermore, the least-squares cross-wavelet analysis is applied to demonstrate how the temperature and precipitation time series can be used for assessment of the results. The proposed method is fast, does not rely on any interpolation methods, leaves the data gap as is, considers the observational uncertainties, and does not depend on thresholds.

1. Introduction

Real-time monitoring of the Earth’s future changes is crucial for raising awareness and decreasing negative effects on natural resources and humans [1,2,3]. For example, wildfire, a natural component of forest ecosystems, can create significant threats to public safety, property, and forest resources [4]. New advances in remote sensing technology using satellite imagery have made it possible to sense and track the changes spatially and temporally [5,6]. Depending on the purpose of monitoring, an appropriate spectral band combination from satellite imagery may be used to generate a time series for a region of interest from which an effective change detection algorithm can be used to detect and classify disturbances in newly acquired data [7,8].

Normalized Difference Vegetation Index (NDVI) [9] and Enhanced Vegetation Index (EVI) [10] time series are widely used for monitoring vegetation dynamics in order to detect disturbances caused by several factors, such as fires, deforestation, drought, flood, and insect attack [11,12,13,14]. EVI can be conceived as an improved NDVI which can attenuate canopy background and atmospheric noises and does not saturate in high biomass regions [15]. The EVI range is typically from $-1$ to 1 and varies between $0.2$ and $0.8$ for healthy vegetation [16]. Changes within vegetation time series can be categorized into three groups: (1) seasonal or periodic changes caused for instance by annual temperature and precipitation influences, (2) gradual trend changes caused for instance by inter-annual climate change or gradual recovery after fires, etc., and (3) abrupt (sudden) changes caused for instance by deforestation, wildfire, insect attack, etc. [11,12].

There are three major challenges when it comes to real-time monitoring of vegetation dynamics using satellite image time series. First, the cloud cover is the main factor that must be considered when acquiring optical images. Furthermore, the availability of satellite imagery can vary spatially and temporally due to various factors, such as the poleward convergence of satellite orbits, sensor defects, data reception, and acquisition [17,18,19]. Therefore, the vegetation time series are often unequally spaced with possibly long data gaps. In particular, in high latitude regions, vegetation time series usually have consistent long gaps, i.e., no valid observations are available in a certain season (e.g., November to April) every year due to the availability of satellite imagery, cloud cover, snow, etc. Interpolation and gap-filling methods based on the normal data variation assumption usually perform poorly and can prevent the change detection [20,21,22].

Second, the time series values have uncertainties that must be properly modeled to reduce the effect of noisy measurements. Observational uncertainties in the presence of undesired cloud/shadow, smoke, haze, snow, and ice are usually provided by an effective cloud-masking procedure [23,24,25]. This means that the vegetation time series are unequally weighted, i.e., the time series is associated with a covariance matrix. Considering the covariance matrix of observations can greatly improve the signal estimation [26,27].

Third, an effective change detection method needs to be threshold independent and be robust against random noise and seasonality to avoid the complexity of the method when applying it to various regions of interest [8,28,29]. An effective approach is to find a Stable History Period (SHP), i.e., a time period when the components of the vegetation time series remain stable (without any disturbance). The data in SHP can then be used to estimate the normal expected behaviour of the data in a Monitoring Period (MP) for disturbance detection [7,30,31].

To overcome these three challenges, a robust monitoring method is developed herein with the main focus on analyzing unequally spaced time series with consistent long gaps to detect and classify disturbances (abnormal changes) at the end of the time series. Toward this development, the following four important questions will also be answered:

• Can SHP representing both inter-annual (gradual) and intra-annual (seasonal) variations be determined within the time series?
• How reliable is SHP for detecting a disturbance in MP?
• How much the weights obtained from the observational uncertainties can increase the probability of a disturbance detection in MP?
• How additional data sets, such as temperature and precipitation, can be used for assessment?

The proposed monitoring method will be tested using simulated EVI time series to provide general solutions to the first three questions raised above. Then, using real-world time series these four questions will be discussed further. The simulated EVI time series are generated with consistent data gaps and many missing observations and by varying the amount of noise, seasonal variations, and adding abrupt changes with different magnitudes to represent a wide range of vegetation dynamics, where magnitude is defined as the vertical difference of the step change [7]. The real-world EVI time series are obtained from Landsat 8 composites since 2013 covering a region in northern Alberta, Canada (about $\mathrm{160,000}$ ha), where there are consistent data gaps from November to April each year with many other missing observations due to the climate and sensor conditions. Most parts of the region (toward the west) have been impacted by massive wildfires in May 2019, known as the Chuckegg Creek Fire [32]. Therefore, MP intentionally starts from April to assess the sensitivity and specificity of the proposed method.

2. Materials and Methods

2.1. Study Region

The study region is part of the Canadian boreal forest located in the southwest of High Level in northern Alberta, which has an area of approximately $\mathrm{160,000}$ ha, see Figure 1. The Canadian boreal forest stores approximately $11\%$ of the world’s total carbon, purifies the air and water and regulates the climate [33]. The predominant species in the region include lodgepole pine, black spruce, white spruce, balsam fir, trembling aspen, balsam poplar, and white birch. The prominent coniferous evergreen trees (trees with needles), such as lodgepole pine, black spruce, white spruce, and balsam fir typically burn five to ten times faster than other trees due to their tendency to grow much closer together [34]. Many trees toward the west and center of the study region burned in May 2019 according to the reports available at the Alberta’s wildfire website [32]. However, the northeast part of the study region, i.e., the northeast of AB-58 and east of AB-35 in the region, was not impacted by wildfire during 2019.

2.2. Data Sets and Pre-Processing

2.2.1. Satellite Imagery

Landsat 8 surface reflectance products (30 m spatial resolution) are downloaded from the U.S. geological survey (USGS) website [35]. They are atmospherically corrected using auxiliary data regained from the Moderate Resolution Imaging Spectroradiometer (MODIS), such as water vapor, ozone, and aerosol optical thickness, as well as a digital elevation derived model. Then the blue, red, and near-infrared bands are used to obtain the EVI images [10,36]. The start date of Landsat 8 images is April 2013. It is worth mentioning here that the time distribution of available satellite images covering this study region is every seven and nine days consecutively with many missing images. This naturally makes the vegetation time series unequally spaced, where almost all current change detection techniques can fail to correctly detect disturbances.

To increase the temporal resolution, other sensor data may be included, such as Landsat 7 Enhanced Thematic Mapper (ETM+). However, since Landsat 8 Operational Land Imager (OLI) has improved calibration and signal-to-noise characteristics compared to Landsat 7 (ETM+) and with different sensitivity to surface reflectance and atmospheric state, it is difficult to properly link their EVI values [37]. This may cause undesired biases in the EVI time series which prevents rigorous testing of the proposed monitoring method. Therefore, the use of other sensor data is not in the scope of this research and is subject to future studies.

The observational uncertainties herein are calculated from the pixel quality assessment (PQA) band obtained from the Function of Mask (FMask) algorithm for each satellite image [23,38]. According to the product guide for Landsat 8, if a PQA pixel value is not equal to 322, 386, 834, 898, or 1346, then there is a chance for existing cloud/shadow, snow, ice, smoke, or haze. Herein, the EVI images are downsampled to a resolution of 300 m for the sake of visualization purposes and computational efficiency. The downsampling of images is an optional process and should not be considered to be a limitation of the proposed monitoring method.

To downsample an EVI image, the average EVI for each non-overlapping spatial square window of size 100 pixels is computed to produce a lower resolution image with the same area coverage. A statistical weight, defined by $1/\left(1+\mathsf{\gamma }{\mathsf{\lambda }}^2\right)$, is also associated with each spatial window, where $\mathsf{\lambda }$ is defined as the number of contaminated pixels (from the PQA band) within each spatial window divided by the size of the window and $\mathsf{\gamma }$ is an adjustment factor. Therefore, a stack of weighted EVI images with a spatial resolution of 300 m will be obtained for the study region. By recording the pixel values along a vertical traverse of the stack, a weighted time series (a per-pixel time series) will be obtained to monitor the changes within a square region of 9 ha. The adjustment factor $\mathsf{\gamma }$ can be estimated from the response of a vegetation index to atmospheric and background noises. For example, since NDVI values are usually higher than EVI values for healthy vegetation, cloud contamination can reduce the values of NDVI significantly more than EVI [16,39]. Thus, the statistical weights must be adjusted accordingly for better performance when using these indices.

2.2.2. Validation Data Sets

The historical weather data, used for the coherency analyses, are provided by the Alberta Agriculture and Forestry [40]. The monthly temperature and precipitation data in Alberta are reported in township level (about 10 km resolution) and computed by a mathematical data interpolation procedure that weights up to the eight nearest station observations [41,42]. The weighted average of all weighted EVI images from 2013 to 2019 is computed and shown in Figure 2, representing the vegetation coverage. The areas with EVI less than or equal to $0.2$ are mainly in the town of High Level. Three pixels A, B, and C (300 m resolution) with known characteristics [32] are selected to assess the results of the monitoring method, displayed in Figure 2. Pixels A and B are in forested areas impacted by the Chuckegg Creek Fire and a wildfire started in August 2019, respectively. Pixel C is also in a forested area where no change is reported during 2019.

The geographic coordinates of the ignition areas of wildfires along with the start dates of wildfires, the burned areas in hectare, and the causes are listed in

Table 1. Validation data within the study region provided by the Government of Alberta. Each location (before 2019), is shown by an asterisk (labeled) in Figure 2, impacted by wildfire.

Label	Date of Wildfire	Latitude	Longitude	Area Burned (ha)	Cause
1	2014-05-10	${58}^{\circ }{26}^{\prime }{30}^{″}$	$-{117}^{\circ }{13}^{\prime }{14}^{″}$	0.7	Railroad
2	2015-04-28	${58}^{\circ }{31}^{\prime }{20}^{″}$	$-{117}^{\circ }{04}^{\prime }{59}^{″}$	0.5	Resident
3	2015-05-30	${58}^{\circ }{15}^{\prime }{01}^{″}$	$-{117}^{\circ }{20}^{\prime }{36}^{″}$	7.0	Recreation
4	2015-06-27	${58}^{\circ }{24}^{\prime }{53}^{″}$	$-{117}^{\circ }{17}^{\prime }{31}^{″}$	246.9	Lightning
5	2016-04-23	${58}^{\circ }{32}^{\prime }{03}^{″}$	$-{117}^{\circ }{19}^{\prime }{12}^{″}$	0.5	Recreation
6	2017-09-28	${58}^{\circ }{18}^{\prime }{41}^{″}$	$-{117}^{\circ }{09}^{\prime }{39}^{″}$	0.5	Recreation
7	2018-05-08	${58}^{\circ }{26}^{\prime }{23}^{″}$	$-{117}^{\circ }{09}^{\prime }{49}^{″}$	0.4	Agriculture
	2019-05-12	${58}^{\circ }{17}^{\prime }{16}^{″}$	$-{117}^{\circ }{16}^{\prime }{47}^{″}$	350,135.0	Lightning

. The locations listed in Table 1 are also highlighted by asterisks and labeled in Figure 2. The top-right panels in Figure 2 show parts of a satellite image clipped to the polygon shapefiles for locations 3 and 4, where the polygon shapefiles are provided by the Natural Resources Canada [43]. Due to a small size of burned areas in locations $1,2,5,6$ and 7, no polygon shapefiles are provided. The last row in Table 1 shows some information for the Chuckegg Creek Fire. An aerial view of this wildfire is also illustrated in Figure 2. No polygon shapefiles for the burned areas in 2019 are currently available for further comparisons and validations.

2.3. Monitoring Method

Many change detection methods, such as Breaks For Additive Seasonal and Trend (BFAST) [11,12], Detecting Breakpoints and Estimating Segments in Trend (DBEST) [44], and Seasonal-Trend decomposition procedure based on Loess (STL) [45,46], scan a time series $f\left(t_1\right),f\left(t_2\right),\dots ,f\left(t_n\right)$ to estimate a time $t_j$ when a change occurs, where $f\left(t_j\right)$ denotes the observation (measurement) at $t_j$. For the near-real-time monitoring; however, the goal here is to determine whether the newly acquired observations $f\left(t_{n+1}\right),f\left(t_{n+2}\right),\dots$ in MP are still in agreement with the expected variations of data (forecast) estimated from $f\left(t_1\right),f\left(t_2\right),\dots ,f\left(t_n\right)$ in SHP. Therefore, the monitoring method has two main steps. The first step is the selection of a reliable SHP, and the second step is the detection of a disturbance in the newly acquired observations using the observations in SHP. These two steps are described below in detail.

2.3.1. Stable History Period Selection

To detect possible disturbances at the end of the time series, the first step is the selection of SHP. This means that an integer $\ell \ge 1$ must exist such that the observations $f\left(t_{\ell }\right),f\left(t_{\ell +1}\right),\dots ,f\left(t_n\right)$ have stable components. Furthermore, SHP needs to be long enough (> two years) to consider both seasonal and gradual variations [7,31]. Several techniques may be applied for selecting SHP (choosing $t_{\ell }$), such as the reversed-ordered-cumulative sum of residuals [47], and BFAST which has a higher computational cost [7,11], or the expert knowledge. Watts and Laffan [13] shows that BFAST poorly performs for detecting disturbances caused by fires in semi-arid regions, likely due to a small difference in EVI values before and after the fire, and the disturbance detection is sensitive to the selection of the parameters in the BFAST decomposition which makes this method less reliable for selecting SHP.

Herein, the Anti-Leakage Least-Squares Spectral Analysis (ALLSSA) [48] is implemented to select SHP. ALLSSA is based on the Least-Squares Spectral Analysis (LSSA) [49,50] with an improved model that simultaneously determines an optimal set of functions, such as step, ramp, and harmonic, that fits best the time series. Thus, SHP is the most recent time period when the ALLSSA season-trend model does not break down statistically ([51] Chapter 3). ALLSSA is a robust minimizer technique that performs like the simulated annealing [52] and the Least Absolute Shrinkage and Selection Operator (LASSO) [53,54]. ALLSSA analyses unequally spaced time series without any need for interpolation and gap-filling, and it can also consider the observational uncertainties. ALLSSA accurately estimates the trend and seasonal components whose frequencies can be any real numbers without the problem of under/over-fitting unlike the Ordinary Least-Squares (OLS) regression method [51,55].

2.3.2. Near-Real-Time Disturbance Detection

Assume that SHP is determined, and the trend and seasonal components of the series within that period are successfully estimated. Now, the question here is: Do the new incoming observations in MP comply with the expected behaviour of the historical series? To answer this question, one may first obtain the residual series of the stable historical series and newly acquired observations, and then determine whether the residual series fluctuates randomly or not by a rigorous statistical inference as described below.

Suppose that ${\mathbf{f}}_n={\left[f\left(t_1\right),f\left(t_2\right),\dots ,f\left(t_n\right)\right]}^{\mathrm{T}}$ is the stable time series in the vector form, where $\mathrm{T}$ is transpose, and ${\mathbf{C}}_{\mathbf{f}n}$ is the covariance matrix of dimension n associated with ${\mathbf{f}}_n$. Consider ${\mathbf{D}}_n$ is a $n\times m$ matrix whose columns are the trend, cosine and sine functions with corresponding estimated coefficients $\widehat{\mathbf{c}}$ (an $m\times 1$ vector). Now let ${\mathbf{f}}_{n+d}$ ($d\ge 1$) be defined as ${\mathbf{f}}_n$ with additional newly acquired observations $f\left(t_{n+1}\right),f\left(t_{n+2}\right),\dots ,f\left(t_{n+d}\right)$, ${\mathbf{C}}_{\mathbf{f}n+d}$ be the covariance matrix associated with ${\mathbf{f}}_{n+d}$, and ${\mathbf{D}}_{n+d}$ be the same as ${\mathbf{D}}_n$ with additional d rows. Using the same estimated coefficients $\widehat{\mathbf{c}}$ for the time series in SHP, one may compute the residual series as ${\mathbf{r}}_{n+d}={\mathbf{f}}_{n+d}-{\mathbf{D}}_{n+d}\widehat{\mathbf{c}}$, for each $d\ge 1$.

Now if the residual series is random, then it means the new d observations are following the expectation, and there is no disturbance. To evaluate the randomness, one may compute the weighted least-squares spectrum of the residual series [26,49] to determine if there is any statistically significant spectral peak at a certain confidence level within a given frequency range, chosen from $0.1$ to $3.5$ (cycles/year) to account for disturbances in trend and seasonal components. The assumption here is that the residual series follows the multidimensional normal distribution with zero mean and covariance matrix ${\mathbf{C}}_{\mathbf{f}n+d}$ ([51] Chapter 2). Thus, the weighted least-squares spectrum of the residual series follows the beta distribution, and the critical value at $1-\mathsf{\alpha }$ confidence level will be $1-{\mathsf{\alpha }}^{2/(n+d-2)}$, where $\mathsf{\alpha }$ is the significance level, usually $0.01$ or $0.05$ [26,51].

Herein, the statistical weights are considered as a scaled inverse-square of the observational uncertainties (error bars), i.e., the covariance matrix of the observations is diagonal that can be treated as a vector for the sake of computational efficiency. Please note that the relative weights play a significant role in the proposed monitoring method, but the method is not sensitive to any covariance factors [26]. Please also note that defining the weights is not necessarily based on the downsampling, and an effective cloud masking technique can be used to define appropriate statistical weights associated with the original pixels. The main role of these weights is to control the effect of noisy observations during monitoring.

Figure 3 shows the flowchart of the near-real-time monitoring algorithm when applied to satellite imagery. The proposed method can also classify the detected disturbances into abrupt, gradual, and seasonal from the location of the statistically significant spectral peak in the frequency domain when more new observations become available. A statistically significant spectral peak at one cycle per year (annual peak) or after represents a seasonal change while before one cycle per year represents a gradual or abrupt change. If several peaks simultaneously become statistically significant in the spectrum, then the dominant peak may represent the type of a disturbance [48]. The magnitude and direction of the abrupt change can also be more accurately estimated when more new observations are acquired.

2.4. Least-Squares Cross-Wavelet Analysis as an Assessment Method

In many practical applications, one needs to investigate whether there is any relationship between different phenomena and how they interact with each other. For instance, questions such as how the temperature and precipitation can interact with the plant phenology in the forest and semi-arid ecosystems [56,57]. The coherency analysis of the time series can greatly help to answer such questions. The coherency is a measure in the frequency domain which shows how much the seasonal components are related which is different from the correlation that simply provides a possible relationship in the time/space domain [57].

The cross-wavelet transform (XWT) [58,59,60], and the least-squares cross-wavelet analysis (LSCWA) [61] are powerful tools that can be applied for coherency analysis. XWT and LSCWA are robust methods of analyzing two time series together which transform/decompose the time series into the time-frequency domain to search for possible coherency between the components under investigation in two series. Unlike XWT, LSCWA does not necessarily require equally spaced time series, and it can analyze time series with different sampling rates without any need for interpolations or gap-filling. LSCWA can also consider the covariance matrices associated with the time series. The main three properties of LSCWA are listed below.

First, the least-squares cross-wavelet spectrogram (LSCWS), computed by multiplying the spectrograms of the two time series, e.g., EVI and temperature. Each spectrogram is obtained by fitting the sinusoids to each time series segment while simultaneously accounting for other constituents, such as trends [26]. The segmentation depends on the frequency where the segment size decreases as the frequency increases and vice versa, enabling the detection of inter-annual and intra-annual components which may have variability of frequency and amplitude over time [27]. The least-squares cross-spectrum (LSCS) is a special case of LSCWS obtained by analyzing the entire series without any segmentation [61].

Second, the stochastic confidence level surface which can be derived when the time series are obtained from populations of normally distributed random variables. This surface shows whether the coherency between certain components of interest in LSCWS, expressed as percentage variance, is statistically significant at a certain confidence level. In other words, when a spectral peak within a local time-frequency neighborhood in the cross-spectrogram stands above the stochastic surface, the seasonal components of both time series are statistically coherent within that neighborhood.

Third, the phase differences, indicating how much the components in the first time series lag/lead the ones in the second time series. The phase differences are illustrated by arrows on LSCWS which follow the principle of the trigonometric circle [57]. In other words, if an arrow, displayed on LSCWS within a local time-frequency neighborhood, points to the east, then there is no phase difference between the seasonal components of both series in that neighborhood. In addition, if it points to the southeast (fourth quadrant) by an angle, then the seasonal component of the second time series leads the one in the first series by that angle. The angle can be converted to time depending on which cyclic frequency is used to derive it. For instance, an angle 30 degrees at 1 cycle per year is one month lag/lead, and at 2 cycles per year is two weeks lag/lead (approximately), etc.

3. Results

This section demonstrates the results of applying the proposed monitoring method to vegetation time series including assessments. First, the effectiveness of the method is evaluated by simulating EVI time series that have similar structures as the real-world EVI time series for the study region. Then, the method is tested by applying it to the collection of real-world EVI time series with a spatial resolution of 300 m as described in the previous section, and a comparison is made between NDVI and EVI time series. The SHP selection results are also validated by the per-pixel EVI time series of resolution 30 m, affected by the historical wildfires (see the details in Table 1). Finally, the monitoring results of three EVI time series, each corresponding to a square area of 9 ha, are assessed by LSCWA, also showing the possible coherency between the EVI time series and temperature/precipitation time series.

3.1. Simulation Experiment for Validation of the Monitoring Method

The sensitivity (true-positive rate) and specificity (true-negative rate) of the proposed monitoring method are examined by simulating EVI time series. First, the seasonal component of each time series is generated using an asymmetric Gaussian function for each season [11] whose amplitude is listed in

Table 2. Parameters for simulation of vegetation time series: amplitude of the asymmetric Gaussian function, noise level, magnitude of the abrupt change, and number of new observations to quantify detection delay (d).

Parameters		Values
amplitude		$0.1,0.2,0.3,0.4,0.5$
noise level		$0.005,0.01,0.015,\dots ,0.15$
magnitude		$0.0,-0.1,-0.2,-0.3$
d		$1,2,3,4,5,6$

. Second, a linear trend with an abrupt change is added such that the change occurs right after the start date of MP. The magnitudes are listed in Table 2, where zero magnitude means no abrupt change herein. Third, two independent noises are generated using random number generators following the normal distributions with zero means and standard deviations ${\mathsf{\sigma }}_1$ and ${\mathsf{\sigma }}_2$.

More precisely, let $f\left(t_{\ell }\right)=S\left(t_{\ell }\right)+\mathsf{\Gamma }\left(t_{\ell }\right)+{\mathsf{\eta }}_1\left(t_{\ell }\right)-\left|{\mathsf{\eta }}_2\left(t_{\ell }\right)\right|$, $1\le \ell \le n$, where S is an asymmetric Gaussian function, $\mathsf{\Gamma }$ is the trend function, ${\mathsf{\eta }}_1$ and ${\mathsf{\eta }}_2$ are two independent random noises of the same type as described above, and $|·|$ is the absolute value. The noise level less than $0.025$ means that ${\mathsf{\eta }}_2$ is omitted and the time series is treated as equally weighted. This is the approximate noise level that often remains after cloud-masking and atmospheric correction processes, making it very challenging to define appropriate statistical weights [12]. However, the noise level greater than $0.025$ means that ${\mathsf{\eta }}_1$ with ${\mathsf{\sigma }}_1=0.025$ is first added to the time series, and then $|{\mathsf{\eta }}_2|$ with ${\mathsf{\sigma }}_2>0$ is subtracted from the new time series, representing the uncertainties that one may obtain from a cloud-masking algorithm. For example, noise level $0.1$ means adding ${\mathsf{\eta }}_1$ with ${\mathsf{\sigma }}_1=0.025$ and then subtracting $|{\mathsf{\eta }}_2|$ with ${\mathsf{\sigma }}_2=0.075$ from the series. Since in forested and vegetated lands, the vegetation index values are positive, the presence of clouds typically decreases the index values, the reason $|{\mathsf{\eta }}_2|$ is subtracted from the simulated time series [16,39]. In this simulation experiment, the statistical weights are also defined by $1/\left(1+100{\mathsf{\eta }}_2^2\right)$, shown to be effective.

For each noise level listed in Table 2, 1000 unequally spaced time series are simulated with various SHPs and consistent data gaps from November to April to calculate the probability of change detection. Figure 4 shows one of such simulated time series when the amplitude is $0.3$, magnitude is $-0.2$, and the noise level is $0.05$. In this figure, SHP is accurately determined using ALLSSA, and then both OLS with three cyclic frequencies $1,2,3$ cycles per year (in red) and ALLSSA (in blue) are applied to estimate the season-trend components, clearly showing the over-fitting issue of OLS which causes significant biases in predicting values within the gaps (yellow backgrounds).

In general, it is not recommended to start MP within a long time period without any historical data available in that period. In this simulation experiment, since there are consistent gaps from November to April in SHP (yellow backgrounds), MP should not be started during November–April because it can easily increase the false-positive rate no matter which regression technique is used. The reason is due to the fact that the true missing values within the gap periods may be any numbers, e.g., because of harvesting activities, urbanization, etc.

Figure 4 also demonstrates how the predicted time series in MP is used to detect disturbances, compare the values of the solid red and blue lines (with the gray background, representing MP) with the values of black diamonds. Here an abrupt change with magnitude $-0.2$ exists right after the start of MP (June 2019) which after three observations ($d=3$) the change is detected (the green vertical line) from the process demonstrated in Figure 5. The residual of the time series shown in Figure 4 is illustrated in Figure 5a. All the spectral peaks, shown in Figure 5b,c, stand below the critical value at $99\%$ confidence level (the red lines), indicating the random behavior of the residual series. However, when three observations are available (i.e., $d=3$), the spectral peaks at the lower frequencies become statistically significant, see Figure 5d. Figure 5e is obtained by finding the maximum values of the spectral peaks for each residual series corresponding to $d=1,\dots ,6$, and it shows three observations (circled) are required in this case to detect the abrupt change.

Figure 6 shows the probability of change detection when the seasonal amplitude is $0.3$. The results for other amplitudes are very close and not shown here for brevity. One can observe that if the noise level is $0.1$, then an abrupt change with magnitude $-0.2$ can be detected by the weighted method with a probability of $0.7$ when three new observations are available, i.e., $d=3$ (see the solid green lines in the bottom-left panel). The top-left panel shows that the weighted method has a high true-negative rate (above $90\%$) for all d values. The results also signify that the sensitivity and specificity of the weighted method are higher than the unweighted method (without considering the statistical weights). Furthermore, the difference between the probabilities of change detection using the weighted and unweighted methods for each d generally increases as the noise level increases.

Various SHPs are randomly selected in these simulated time series to get a whole picture of a wide range of vegetation dynamics. However, the longer SHP is, the higher the method sensitivity will be, with an insignificant change in specificity. As the noise level increases, having a six-year SHP can increase the probability of disturbance detection in MP up to $30\%$ compared to having only a two-year SHP (the results are not shown here for brevity). This improvement is due to the larger degrees of freedom and more accurate prediction of gradual and seasonal components in MP [62]. Though SHP longer than six years may further improve the probability of disturbance detection, a six-year SHP is still pretty efficient.

3.2. Disturbance Detection in the Study Region

Using Landsat satellite imagery, it is shown how the weighted method detects the disturbances mainly caused by fire during MP (April–November, 2019) in the study region. The monitoring method is applied to all per-pixel weighted EVI time series with a spatial resolution of 300 m. Only the areas with the weighted average EVI greater than $0.2$ are selected for the near-real-time monitoring in 2019, covering approximately $99\%$ of the region, see Figure 2.

Figure 7a illustrates the dates of the most recent disturbances before 2019, estimated by the ALLSSA season-trend model. Twenty-seven hundred per-pixel EVI time series (30 m resolution) are obtained within the known burned areas for model validation. From the PQA bands, all clear observations are selected to obtain these per-pixel EVI time series within all the seven numerically labeled locations shown in Figure 2 to test the SHP selection method without weights. Since the burned areas in locations $1,2,5,6$, and 7 are small, only four pixels (30 m resolution) around the geographical coordinates of these locations are selected for the validation. The method could estimate the dates of disturbances caused by fire for $87\%$ of these EVI time series within two-month accuracy. Furthermore, the weighted EVI time series (300 m resolution) within the polygons are also analyzed separately using the weighted method, and the temporal accuracy of the detected changes was improved due to the inclusion of observations with higher uncertainties.

The pixels with SHPs shorter than two years are excluded from the near-real-time monitoring, see the blank areas in Figure 7b–d. No significant historical changes are detected for approximately $10\%$ of the study region, see the black pixels in Figure 7a. Though SHPs for these areas can be longer than six years, the trend and seasonal components can still be approximated very well for the six-year period, effective for disturbance detection during MP.

Figure 7b shows the spatial variation of the noise level for SHP, i.e., the standard deviation of the stable residual series. The stable historical series are used to detect possible disturbances in MP (April–November, 2019), and the dates of detected disturbances are shown in Figure 7c. Figure 7d shows the change magnitudes estimated by the ALLSSA season-trend model, i.e., the vertical difference of the step change. Approximately $62\%$ of the study region with the weighted average EVI greater than $0.2$ has SHPs greater than two years. No disturbance is detected for $13\%$ of the region with SHPs longer than two years, green areas shown in Figure 7c.

The more detailed analyses of the time series corresponding to the pixels A, B, and C (300 m resolution) are illustrated in Figure 8. The disturbance detected after one observation with magnitude $-0.21$ in Figure 8a (indicated by green line) is caused by the wildfire in May 2019. The disturbance detected after four observations with magnitude $-0.26$ in Figure 8b (green line) is also caused by the wildfire in August 2019. Indeed, the weighted method detected the disturbances in pixels A and B after one new observation is acquired, although the simulation experiment in Figure 6 shows approximately $25\%$ chance for disturbance detection from one new observation ($d=1$) when magnitude is $-0.3$. Low noise level (<0.05), long SHP (>four years), and low uncertainty in the new observation are main factors for quickly detecting the disturbances here. Most parts in the north and west sides of the study region, shown in orange and yellow in Figure 7c, are also impacted by the wildfire in May 2019, and the detection delays are due to the method sensitivity and missing values (poor temporal resolution). In general, the areas with lower noise levels and longer SHPs provide more reliable change detection results.

No disturbance is detected at pixel C by the monitoring method. The ALLSSA season-trend model provides a decent fit within the time periods when almost clear observations are available (April–November); however, the prediction within the gaps (yellow backgrounds) in Figure 8c does not agree with reality. Other season-trend models, such as OLS and LASSO also provide poor predictions in this case, not shown here. When all Landsat EVI images are used, four more observations fall within the gaps with relatively much larger uncertainties, shown in red in Figure 8c; however, the method cannot still provide a good prediction within the gaps (green dashed lines). Therefore, MP starts from April here (gray background).

To compare the change detection results using different vegetation indices, the weighted per-pixel NDVI time series are also calculated for pixels A, B, and C (300 m resolution), and the monitoring results are illustrated in Figure 9. SHPs and the detected disturbances shown in Figure 9a,c agree with the ones shown in Figure 8a,c. However, SHP is not correctly determined for pixel B. The ALLSSA season-trend model breaks down at the end of the historical time series, i.e., the model detects the observation pointed by the red arrow in Figure 9b as an outlier. This is mainly due to the presence of canopy background noise (e.g., snow and ice) that significantly reduces the NDVI value, and also the FMask algorithm has not properly identified the presence of clouds, snow, and ice on the 6 of November 2018, see a view of the clipped satellite image for pixel B in Figure 9b. Therefore, the statistical weight defined from the FMask algorithm could not properly reduce the effect of this observation. EVI, however, did not have this issue mainly because EVI can attenuate canopy background and atmospheric noises and remain sensitive to canopy variations [16,39]. After ignoring this observation, ALLSSA provided a new season-trend fit, shown by the red dashed lines in Figure 9b. Then, the disturbance in MP is currently detected.

The computational complexity of the proposed method depends on several factors including the time-frequency resolution, covariance matrix (scale, diagonal, fully populated), hardware, and implementation. The entire process of the monitoring method including the downsampling of the clipped EVI images to obtain Figure 7 took approximately 15 min in MATLAB on an average normal computer. The SHP selection process and monitoring a weighted time series of size 100 may take approximately $0.05$ s in MATLAB and $0.3$ s in Python. Also, only the monitoring (SHP is given) of a weighted time series of size 100 may take $0.001$ s in MATLAB and $0.006$ s in Python [63].

3.3. Assessment of the Results Using the Coherency Analyses

The results of the coherency analyses between the EVI time series and each of temperature and precipitation time series corresponding to pixels A, B, and C are illustrated in Figure 10. LSCS and LSCWS of the EVI and temperature time series (labeled) show significant annual coherency, see the peaks at one cycle per year. However, the percentage variances of the annual spectral peaks in LSCWS in Figure 10a,b significantly decrease in 2019, indicating that a disturbance has occurred in either the EVI or temperature time series. From the spectrogram of the temperature time series (not shown here), since the annual components of the temperature time series remain stable over time, the significant reduction of the percentage variance in 2019 is then due to a disturbance occurred in the EVI time series. Furthermore, the phase differences, represented by arrows in LSCWS, show that the annual components of the temperature series lead the ones in the EVI time series by approximately 15 degrees (two weeks).

The right panels in Figure 10 show that the annual and semi-annual components in the EVI time series are statistically coherent with the ones in the precipitation time series at $99\%$ confidence level. Similarly, the significant reduction of the percentage variances of the spectral peaks in LSCWS, corresponding to both annual and semi-annual components for pixels A and B, is due to the disturbances that occurred in 2019. The direction of the arrows indicates that the annual components of the precipitation time series lead the one in the temperature time series by approximately 15 degrees or two weeks. Also, the semi-annual components of the precipitation time series lead the ones in the EVI time series by approximately 60 degrees (one month) during 2017–2019. This loosely indicates that the peak of EVI comes a few weeks after the peak of precipitation possibly due to the plant root mechanism and water absorption [34,64].

Table 3. The ALLSSA simultaneous estimated components of the weather time series for pixels A, B, and C displayed in Figure 2 and Figure 7. The units of temperature, precipitation, and frequency are Celsius, mm, and cycles per year, respectively.

Pixel	Time Series	Frequency	Amplitude	Phase	Intercept	Slope
A	Temperature	$1.00$	$17.57\pm 0.42$	$-1.55\pm 0.00$	$-0.08\pm 0.59$	$-0.03\pm 0.15$
A	Precipitation	$1.00$	$28.03\pm 3.21$	$-1.59\pm 0.00$	$36.56\pm 4.52$	$-1.25\pm 1.13$
		$2.01$	$13.49\pm 3.20$	$1.64\pm 0.24$
B	Temperature	$1.00$	$17.99\pm 0.42$	$-1.55\pm 0.00$	$0.24\pm 0.59$	$-0.02\pm 0.15$
B	Precipitation	$1.00$	$18.30\pm 2.43$	$-1.68\pm 0.03$	$27.86\pm 3.41$	$-0.78\pm 0.85$
		$2.01$	$9.98\pm 2.41$	$1.70\pm 0.24$
C	Temperature	$1.00$	$18.05\pm 0.41$	$-1.54\pm 0.00$	$-0.34\pm 0.58$	$-0.02\pm 0.15$
C	Precipitation	$1.00$	$21.03\pm 2.84$	$-1.62\pm 0.01$	$31.79\pm 3.99$	$-1.28\pm 0.99$
		$2.01$	$11.95\pm 2.82$	$1.62\pm 0.24$

shows the results of the ALLSSA season-trend model applied to the temperature and precipitation time series for pixels A, B, and C displayed in Figure 2 and Figure 7. The errors are estimated by the covariance law of propagation [62,65]. The linear trends in the temperature time series for all the locations have almost zero slopes. However, the negative slope of the linear trend in each of the precipitation time series (gradual drought stress) is one of the reasons why the wildfire could progress so fast in the region.

The coherency analysis may also provide an idea of what EVI range one may expect within the gaps (November–April). In other words, since the annual components of the temperature and precipitation time series are coherent with the ones in the EVI time series, a relatively lower EVI range is expected from November to April.

4. Discussion

Based on the results in the previous section, the four questions raised in the introduction are discussed here. First, the simulated and real-world examples showed that SHP can represent the gradual and seasonal variations to predict the expected values from which a disturbance may be detected in the near future. Longer SHP may provide more reliable change detection results. SHP may not always be sufficiently long enough (e.g., <two years). For instance, the areas impacted by wildfire in 2019 cannot provide a long enough SHP for disturbance detection in 2020. The expected values, however, may be obtained by expert knowledge. For example, one may use the recovery time series from similar disturbance sites in previous years to use them for prediction.

Second, the reliability of SHP depends on the derivation methodology. Estimating the statistically significant season-trend components of time series by robust methods, such as ALLSSA and LASSO, improves the under/over-fitting issues which can increase the false-positive/negative rates. These methods can be applied to many applications since they are data-driven. It is shown in Figure 4 that when the historical data have consistent long gaps within a certain period, then the prediction of the values within the gaps can be improved by ALLSSA. However, it is still not recommended to start MP within that particular time period. The real-world EVI time series within the known burned regions also showed about $87\%$ accuracy for selecting SHPs by the ALLSSA season-trend model. The capacity of the real-time monitoring method in detecting a disturbance depends on the signal-to-noise ratio of the historical series. From Figure 7b–d, one can observe that the method could detect disturbances caused by wildfire, started in May 2019, within three new observations for most areas with lower estimated noise level (see toward the center of the region). Furthermore, the range of noise levels in Figure 7b, estimated by the ALLSSA season-trend model for SHP, confirms the simulated one in Figure 6 for appropriate testing of the proposed method. The simulated results in Figure 6 showed that at least three new observations are required to have a probability of more than $0.7$ to detect an abrupt change with magnitude $-0.2$ when the noise level is less than $0.1$.

Third, the simulated results in Figure 6 showed that the weighted method can increase the sensitivity and specificity of the method to disturbance detection up to $20\%$. The weights can be obtained by inverting the covariance matrix associated with the time series. Generally, if an observation has a higher uncertainty (lower weight), then its effect will be reduced in the estimation of components. The atmospheric processes in a region within a short time period tend to look alike (e.g., due to presence of clouds and smoke), and so the covariances between observations may also be estimated to further improve the results. As shown in Figure 8 and Figure 9, EVI time series appears to perform better than NDVI time series for near-real-time monitoring, yet both indices complement each other in global vegetation studies and improve upon the detection of vegetation changes.

Fourth, it is shown in Figure 10 that the temperature and precipitation time series can be used for assessment. When the coherency significantly deviates from its expectation, it means that a disturbance is present in either of the time series. The spectrograms can then show in which time series the disturbance occurs. LSCWA not only shows the coherency between the components of two time series but also shows how much the components of the two time series lead/lag from each other. The additional data sets, such as temperature and precipitation, can help to improve the prediction. However, note that the relationship between climate and vegetation may not always be straightforward especially when it comes to more complex land covers. LSCWA may help to investigate possible coherency between the components of two series, but it does not necessarily mean that such relationships must link different phenomena together.

To improve the temporal resolution of image time series using additional data sets, several methods have been proposed. For example, Zhu et al. [8,55] show how to combine all available Landsat data for change detection and classification and for generating synthetic Landsat images. Also, Reiche et al. [30,31] combine Sentinel-1 time series with Landsat and ALOS (Advanced Land Observing Satellite) PALSAR (Phased Array type L-band Synthetic Aperture Radar) to improve near-real-time deforestation monitoring and emphasize that change detection is always linked with an exchange between the temporal, producer’s, and user’s accuracy. However, the computational cost of these methods can be very expensive, and their performance is yet to be assessed in high latitude regions with consistent clouds, ice, and snow.

5. Conclusions

In this contribution, a near-real-time monitoring method is proposed which requires a stable history period to forecast the time series in the near future (expectation) when a new disturbance can be detected. It is a multi-purpose globally applicable and fast alerting method, does not depend on thresholds, and is capable of analyzing unequally spaced and weighted time series. The effectiveness of the method for detecting disturbances in near-real-time was assessed by both simulated and real-world vegetation time series, showing promising results. The MATLAB codes for LSSA, ALLSSA, and LSCWA are freely available [66]. The MATLAB and Python codes for the proposed monitoring method will also be freely available soon [67].

The proposed monitoring method may be applied to other applications. For example, in farming applications with crop rotation [68,69], agronomists may estimate what is expected in the upcoming season from auxiliary data sets, such as an elevation model, amount of applied fertilizers, and precipitation (historical and predicted). Thus, the proposed method can be applied to verify whether the incoming observations meet the expectation. In other applications, such as Interferometric Synthetic Aperture Radar (InSAR), the method may be applied to monitor interseismic, aseismic, and slow slip activities [70,71,72]. Considering the temporal covariances between observations may also improve retrieving small magnitude surface displacement signals from the InSAR time series and improve the disturbance detection. The coherency analyses presented herein may also be used for creating synthetic satellite images or other applications.