Statistical Analysis of Changes in Sentinel-1 Time Series on the Google Earth Engine

Abstract

Time series analysis of Sentinel-1 SAR imagery made available by the Google Earth Engine (GEE) is described. Advantage is taken of a recent modification of a sequential complex Wishart-based algorithm which is applicable to the dual polarization intensity data archived on the GEE. Both the algorithm and a software interface to the GEE Python API for convenient data exploration and analysis are presented; the latter can be run from a platform independent Docker container and the source code is available on GitHub. Application examples are given involving the monitoring of anthropogenic activity (shipping, uranium mining, deforestation) and disaster assessment (flash floods). These highlight the advantages of the good temporal resolution resulting from cloud cover independence, short revisit times and near real time data availability.

1. Introduction

The Sentinel-1a and -1b Synthetic Aperture Radar (SAR) satellite platforms, with spatial resolutions as low as 20 meters depending on acquisition mode, combined revisit times of the order of six days and complete independence from solar illumination and cloud cover, provide an attractive source of Earth observation data for change detection tasks. The Google Earth Engine (GEE) [1,2] includes, in its extensive and up-to-date archive, Sentinel-1 data in dual polarization (vertical transmission and vertical and horizontal reception over land surfaces) multi-look intensity format. The GEE makes available not only near real time data access but also a very powerful application programming interface (API) for processing and visualizing the data. The GEE API is presently written in JavaScript for direct interaction with the web-based GEE Code Editor and in Python for data analysis outside the GEE web environment. Relatively little work [3] has been published till now involving the use of the GEE’s Sentinel-1 archive for change detection. Since single look complex (SLC) data are not included in the archive for technical reasons, commonly used iterferometric coherence methods are unavailable. Change detection with GEE Sentinel-1 imagery in forest land cover based on Otsu thresholding [4] and K -means clustering is investigated in [5]. The use of multitemporal SAR for automated monitoring of coastal structures with the COSMOS-SkyMed platform is described in [6], and in [7] long-term flood monitoring in Greece is investigated with time series involving a combination of COSMO-SkyMed, Sentinel-1 and Landsat 8 imagery.

In their original publication on polarimetric SAR data analysis, Conradsen et al. [8] introduced a change detection procedure for multi-look SAR data involving a test statistic for the equality of polarimetric covariance matrices assumed to follow a complex Wishart distribution. The procedure is capable of determining, on a per-pixel basis, if a change at any prescribed significance level has occurred in two SAR images. Single polarization (power data, dimensionality p = 1), dual polarization (for example vertically polarized transmission, vertical and horizontal reception, p = 2) and full or quad polarimetric data (all four combinations of vertical and horizontal transmission/reception, p = 3) can be analyzed with this method. More recently [9,10] the procedure was extended significantly to deal with multi-temporal polarimetric SAR data in order to determine the points in time at which change occurs by applying omnibus tests and their factorizations. These latter publications describe test statistics and their distributions for time series of full (2 × 2 or 3 × 3) covariance matrices. For alternative approaches to polarimetric SAR change detection, see [11,12].

The GEE Sentinel-1 intensity data can only be cast into 2 × 2 covariance matrix format without the complex off diagonal elements, as these are not available in the archive. Since the diagonal form is not Wishart distributed, the originally developed algorithms are not immediately applicable. This restriction was removed in [13] where the theory was extended to include covariance matrices in block diagonal form, the GEE data being a special case. In addition, an interpretation of direction of change based on the so-called Loewner order was recently added to the method [14,15].

In this contribution we first describe briefly, in Section 2 below, the relevant aspects of the multi-temporal Wishart-based change detection algorithm, with emphasis on the dual polarization, diagonal only covariance matrices. In Section 3, details of the analysis software (available on GitHub) are given, including a platform independent graphical interface for convenient data exploration and generation of change maps. Then in Section 4 we present some examples intended to demonstrate the advantages of the algorithm’s application to long time series of SAR imagery extracted from the GEE archive. Some conclusions are drawn in Section 5 and a critical discussion is given in Section 6.

2. Theory

The observed fully polarimetric SAR signals, when expressed in covariance matrix form and multiplied by the (equivalent) number of looks, are complex Wishart distributed for fully developed speckle. This distribution is the multivariate complex analogue of the chi squared distribution and is completely determined by the dimensionality p, the equivalent number of looks n, and $\mathsf{\Sigma }$, the complex variance-covariance matrix. In a time series of k observations ${\widehat{\mathsf{\Sigma }}}_i$ of the covariance matrix, the likelihood ratio test statistic Q for rejecting the null hypothesis that all observations were sampled from the same distribution (no-change hypothesis) in favour of the alternative that at least one change occurred is given by [9]

$$lnQ=n\left(pklnk+\sum _{i=1}^{k}ln\left|{\mathit{X}}_i\right|-kln\left|\sum _{i=1}^k{\mathit{X}}_i\right|\right),$$

(1)

where ${\mathit{X}}_i=n{\widehat{\mathsf{\Sigma }}}_i$. This is referred to as an omnibus test. The expression $|·|$ denotes determinant. The probability of finding a smaller value of $-2\rho lnQ$ than an observed value z can then be shown to be approximately

$$\mathrm{Pr}(-2\rho lnQ\le z)\simeq \mathrm{Pr}({\chi }^2\left(f\right)\le z)+{\omega }_2\left(\mathrm{Pr}({\chi }^2(f+4)\le z)-\mathrm{Pr}({\chi }^2\left(f\right)\le z)\right),$$

(2)

where ${\chi }^2\left(f\right)$ is the chi square probability distribution with f degrees of freedom and where the parameters $f,\rho$ and ${\omega }_2$ are determined by $p,n$ and k. For the GEE Sentinel-1 imagery, the diagonal-only covariance matrix is not Wishart distributed, however the test statistic Q is still given by Equation (1) but with special values for the parameters $f,\rho$ and ${\omega }_2$, see [13] for the details.

The omnibus test statistic Q can be factored into a sequence of test statistics $R_j$ for null hypotheses of the form $H_{0,j}:{\mathsf{\Sigma }}_1={\mathsf{\Sigma }}_2=\cdots ={\mathsf{\Sigma }}_{j-1}={\mathsf{\Sigma }}_j$ against the alternative $H_{1,j}:{\mathsf{\Sigma }}_1={\mathsf{\Sigma }}_2=\cdots ={\mathsf{\Sigma }}_{j-1}\ne {\mathsf{\Sigma }}_j$, thus testing whether the first change occurred between the $(j-1)$th and jth observations, $j=2\dots k$. The likelihood ratio test statistic $R_j$ is given by

$$ln{R}_j=n\left(p(jlnj-(j-1)ln(j-1)\right)+(j-1)ln|\sum _{i=1}^{j-1}{\mathit{X}}_i|+ln|{\mathit{X}}_j|-jln|\sum _{i=1}^j{\mathit{X}}_i\left|\right)$$

(3)

and $Q={\prod }_{j=2}^k{R}_j$. Furthermore the statistics $R_j$ are stochastically independent under the null hypothesis. Similarly to the omnibus test statistic, the approximate probability distribution of $-2{\rho }_{j}ln{R}_j$ is given by

$$\mathrm{Pr}(-2{\rho }_{j}ln{R}_j\le z)\simeq \mathrm{Pr}({\chi }^2\left(f\right)\le z)+{\omega }_{2j}\left(\mathrm{Pr}({\chi }^2(f+4)\le z)-\mathrm{Pr}({\chi }^2\left(f\right)\le z)\right).$$

The quantities $f,{\rho }_j$ and ${\omega }_{2j}$ are functions of $p,n$ and j again with particular values for the GEE diagonal-only covariance matrices [13].

Rejection of the null hypothesis for a particular time interval signals a change, but says nothing about the nature of that change. In [14] a characterization of significant change is suggested based on the Loewner order [16], which calculates the definiteness of the difference between the covariance matrices at successive times. Changes are thus classified as positive definite, negative definite or indefinite.

3. Materials and Methods

Software for the sequential SAR change detection algorithm described in the preceding section is provided in Matlab (https://people.compute.dtu.dk/~alan/software.html) and in Python (https://mortcanty.github.io/CRC4Docker/), see also [17]. An additional Python implementation written specifically against the GEE Python API is available on GitHub (https://github.com/mortcanty/EESARDocker) and runs in a Docker container serving Jupyter notebooks to the user’s local browser. The pre-built Docker image can be pulled and run directly from Dockerhub and is platform independent. Since the algorithm runs on the GEE servers, the client-side requirements are light and a Raspberry Pi Docker image is also available. Installation instructions are given in the Github ReadMe file.

Let the index ℓ denote subsequences of the images ${\mathrm{Im}}_j$ beginning with the ℓth image in the time series and ending with the last, that is,

$${\displaystyle \phantom{(}\begin{array}{ll}& \ell =1:\{{\mathrm{Im}}_1,{\mathrm{Im}}_2,{\mathrm{Im}}_3,\dots {\mathrm{Im}}_k\}\\ & \ell =2:\{{\mathrm{Im}}_2,{\mathrm{Im}}_3,\dots {\mathrm{Im}}_k\}\\ \ell =& k-1:\{{\mathrm{Im}}_{k-1},{\mathrm{Im}}_k\}.\end{array}}$$

Denote with $Q^{\ell }$ and $R_j^{\ell }$ the test statistics associated with the ℓth subsequence, see Equations (1) and (3). The evaluation strategy for the GEE API is to pre-calculate, for each realization (observation) $z_j^{\ell }$ of $-2{\rho }_{j}ln{R}_j^{\ell }$, an array of associated p-values

$$p_{\ell j}=1-\mathrm{Pr}(-2{\rho }_{j}ln{R}_j^{\ell }\le z_j^{\ell }),j=2\dots k-\ell +1,\ell =1\dots k-1,$$

and, for realizations $y^{\ell }$ of $-2\rho ln{Q}^{\ell }$, the p-value array

$$p_{\ell }=1-\mathrm{Pr}(-2\rho ln{Q}^{\ell }\le y^{\ell })\ell =1\dots k-1.$$

The arrays are calculated via nested iterations over the subsequences.

In a second iteration step a change map is generated by testing the p-values against the desired significance level. The $p_{\ell j}$-values are only tested within a subsequence if the p-value $p_{\ell }$ rejects the null-hypothesis for the entire subsequence (the omnibus test). This guarantees a fixed false positive probability for all registered changes. If, for a given pixel location, $p_{\ell j}$ rejects the null hypothesis, the change map is updated correspondingly and the iteration continues on the subsequence ${\mathrm{Im}}_j,\dots {\mathrm{Im}}_k$. This allows registration of successive changes at a single pixel location, at a spatial resolution of 20 m for the Sentinel-1 data. In a third and final iteration the average covariance matrix over the no-change observations up to but not including the time of a change is accumulated with the provisional means method, see e.g., [17]. This average is subtracted from the covariance matrix immediately following the signaled change and the positive/negative definiteness (or indefiniteness) of the difference is ascertained to establish the Loewner order. (In [14,15] the opposite convention is adopted: the matrix following change is subtracted from the average prior to change). A particularly efficient method involving matrix pivots [15] is applied which avoids determining, in the case of full $(2\times 2)$ and $(3\times 3)$ covariance matrices, the eigenvalues of the complex matrix differences.

A Jupyter notebook (https://jupyter.org/) included with the software displays a simple and convenient graphical interface (Figure 1) allowing easy data exploration of the entire GEE Sentinel-1 archive together with the previewing of change results. Provision is also made for processing the user’s own image collections that have been uploaded to the GEE, e.g., Radarsat-2 quad pol or Sentinel-1 full dual pol data. Change maps and selected spatial subsets of the images can be exported to the user’s assets on the GEE Code Editor or to his or her Google Drive. Detailed instructions are included in a text cell in the notebook.

4. Results

One of the characteristics of the sequential omnibus algorithm when applied to long time series is its ability to pinpoint regions of high anthropogenic activity. Figure 2 is a color-coded change map of the total number of significant changes in a sequence of 74 Sentinel-1 images over the NATO Airbase near Geilenkirchen, Germany, acquired from January 2018 through June 2019. The aircraft movements to and from their parking positions on the tarmac are emphasized. Similar observations of frequent, localized changes are the subject of the first two examples discussed in the following. The other examples involve flood monitoring and deforestation. For investigations of wetlands monitoring with the complex Wishart omnibus method both in full dual polarimteric form as well as with GEE imagery, see [3,18].

4.1. Libyan Maritime Port Activity

In part due to the refugee crisis in the Mediterranean and also because of political instability, ports along the coast of Libya have been the subject of recent attention, many of them having been closed periodically to international traffic. Arrivals and departures of large sea vessels are easily recorded with the algorithm, as can be seen for example from Figure 3 for the port of Tripolis. Positive definite significant differences (arrivals) are in red, negative definite (departures) in green. The number of detected events represents a lower bound on the vessel movements; a single vessel arrival/departure within a six-day interval between consecutive acquisitions will go undetected.

According to marine insurer North P&I Club, (http://www.nepia.com/insights/industry-news/status-of-libyan-ports-starupdatestar/) the Libyan Ministry of Transport closed the port of Tobruk from 27 October 2017 until 22 February 2018, although the nearby oil terminal at Marsa Elhariga remained open. By setting polygons around the two jetty berths at the oil terminal (Figure 4) and measuring the fraction of changed pixels within them, tanker activity was monitored from January 2017 through June 2019, see Figure 5. Two arrivals were recorded during the closure period, namely in the intervals ending 12 December 2017 and 29 January 2018.

One further example: The port of Benghazi was closed in the first part of 2017, reopening on 5 October 2017. This is confirmed in the chart of Figure 6.

4.2. Arms Control and Verification of Non-Proliferation

In the context of international treaties on nuclear non-prolifieration, satellite sensors can be used to monitor peaceful nuclear fuel cycle activities, especially to detect undeclared activity at sites which, due to their remoteness, are not easily subject to regular on-site inspection. We first consider an example of a declared, operational uranium site, the McArthur River Uranium Mine in northern Saskatchewan, Canada, see Figure 7. A change profile for the core site can be seen in Figure 8. It is evident that changes are occurring throughout the observation time interval. This is contrasted with a similar profile for the decommissioned Cluff Lake mining site (Figure 7), also located in northern Saskatchewan. This mine ceased uranium production at the end of 2002 when the ore reserves were depleted. Apart from the two spikes in June and early August, much less activity is evident in its change profile. By observing the spatial distribution of change pixels, the spikes can be attributed to environmental changes, most probably in ground moisture following precipitation.

4.3. Flood Monitoring

SAR imagery, because of its insensitivity to cloud cover, can often complement rapid response disaster management procedures, typically in situations involving sudden, large scale flooding. However many of the water detection methods proposed in the literature using reflected SAR radar signals require a considerable amount of user intervention. For example [19] gives an overview of the variety of radar signal/water surface interactions in flooded areas that must be accounted for. These depend on wavelength, polarization, incidence angle, wind/wave conditions, presence of partially flooded vegetation or man-made structures, etc. A recent example of the use of Sentinel-1 SAR imagery for flood monitoring and integration of results into rapid response concepts is given in [20].

In some cases, short-period change detection with polarimetric SAR may offer a simpler alternative to water detection/classification methods. Sudden wide-spread inundation will generally lead to a significant change in the polarimteric matrices characterizing the reflected signal. If the observation interval is short enough, most of the observed changes will be attributable to flooding, and afflicted areas quickly identified.

Cyclone Idai, recorded as the worst weather-based event to ever occur in the southern hemisphere, made landfall near the city of Beira, Mozambique on 15 March 2019 causing widespread inundation in Mozambique and Tanzania and a death toll of more than 1300. Figure 9 shows a sequence of six change maps generated with a GEE time series of Sentinel-1a and Sentinel-1b images. The reduction of reflectance from the water surfaces in both the VV and VH bands corresponds to a negative definite covariance matrix difference (green pixels, center left), the rapid receding of flood waters to a positive definite change (red pixels, center right).

A similar series of change maps is shown in Figure 10 for Golestan province, Iran, where from mid-March to April, 2019 widespread flash flooding affected large parts of the province. Here again one sees the initial flooding (green) followed by receding water (red), clearly delineating the areas most seriously affected.

4.4. Clear Cut Logging

Another common application for satellite Earth observation is deforestation assessment/control, not only in the Amazon rain forest, but elsewhere. For a multi-sensor and multi-temporal approach to deforestation monitoring see for example [21].

Logging, and especially clear cutting, in Canadian boreal and temperate rain forests [22] continues to be controversial, with the last few percent of remaining old growth under threat. For example in the Nahmint Valley near Port Alberni, Vancouver Island, centuries-old red cedar and Douglas fir trees are still being cut down. Figure 11 shows two Sentinel-2 acquisitions from September, 2017 and June, 2018, where three new clear cuts are visible. A corresponding change map in Figure 12 obtained from a series of 33 Sentinal-1 images from September, 2017 through September, 2018 confirms that the clear cuts can be easily located, both geographically and temporally, by the sequential omnibus method. In [23] an investigation of ALOS-PALSAR polarimetric radar reflectance from reference forest stands and clear cuts in Swedish forests indicates a significant reduction in intensity in the HH and HV bands after cutting. Although this is for L-band, this might be expected to correspond to similar reductions in the VV and VH bands for the GEE Sentinel-1 C-band data, or equivalently to a negative definite change direction. This is confirmed in Figure 13, which also indicates that the bulk of harvesting occurred within just two consecutive 12-day intervals in December, 2017.

5. Conclusions

We have described the fusion of a statistically sound sequential change detection method with the power, convenience and vast data resources of the Google Earth Engine platform. Long time series of pre-processed dual polarization SAR imagery from the GEE archive from all over the globe can be easily accessed and quickly analyzed in a convenient graphical environment. Several illustrative examples were discussed. The software, including source code, is documented and freely available to anyone registered as a GEE user.

6. Discussion

Although we feel that they are outweighed by the advantages mentioned in the preceding Section, there are some drawbacks in the methodology that should be pointed out, and some possible suggestions for future work:

• First of all, and most significantly, the sequential omnibus tests on the GEE are carried out at the nominal scale of the archived Sentinel-1 data (10 m). This is because of the dependence of the Wishart distribution on the equivalent number of looks (ENL). Confining analysis to a single scale precludes leveraging one of the great advantages of the Earth Engine, namely up-scaling to very large geographical regions. One way to mitigate this in future might be to download representative images with well-developed speckle statistics at different scales and then estimate the ENL values off-line, e.g., with the methods given in [24]. Then those values could be hard wired into the GEE code to allow running the algorithm at coarser scales and on larger scenes.
• The change detection algorithm is purely data driven and unsupervised: The physical cause of detected changes must be inferred from the context. Here, the Loewner order discussed in the text can offer additional information.
• It is our experience that very long time series, typically 75 images or more, can lead to stack overflow on the GEE servers. With typically a 6-day temporal resolution this still allows well over a year of continuous observation at any given location.
• The diagonal-only dual polarization matrix format necessitates resorting to the block diagonal version of the algorithm as discussed in the theory section and in [13]. It would be desirable to have access to the full $2\times 2$ dual polarization matrix. We understand that the GEE developers are considering ways to ingest single look complex (SLC) Sentinel-1 imagery, which would solve this problem: The multi-look dual polarization matrix format could then be constructed from the SLC data.
• The GEE archive is updated very quickly, the Sentinel-1 images are available within a few days of acquisition. But for timely disaster assessment this may not be good enough. Thus the tools described here will be useful only in situations which are not extremely time critical.

Comparison of the performance of our methodology with other approaches cited here, for example using confusion matrices or ROC curves, presupposes the availability of reliable ground reference data. Such information is often very hard or even impossible to come by, especially in disaster assessment situations where remote sensing data are perhaps the only data source available to assess changes/damages. We have however made considerable effort to provide the computer software tools in a form (Docker container) which guarantees reproducibility, verifiability and ease of use in any field of application for which multi-temporal polarimetric SAR change detection is considered to be useful.