Performance Study of Landslide Detection Using Multi-Temporal SAR Images

Abstract

This study addresses one of the most commonly-asked questions in synthetic aperture radar (SAR)-based landslide detection: How the choice of datatypes affects the detection performance. In two examples, the 2018 Hokkaido landslides in Japan and the 2017 Putanpunas landslide in Taiwan, we utilize the Growing Split-Based Approach to obtain Bayesian probability maps for such a performance evaluation. Our result shows that the high-resolution, full-polarimetric data offers superior detection capability for landslides in forest areas, followed by single-polarimetric datasets of high spatial resolutions at various radar wavelengths. The medium-resolution single-polarimetric data have comparable performance if the landslide occupies a large area and occurs on bare surfaces, but the detection capability decays significantly for small landslides in forest areas. Our result also indicates that large local incidence angles may not necessarily hinder landslide detection, while areas of small local incidence angles may coincide with layover zones, making the data unusable for detection. The best area under curve value among all datatypes is 0.77, suggesting that the performance of SAR-based landslide detection is limited. The limitation may result from radar wave’s sensitivity to multiple physical factors, including changes in land cover types, local topography, surface roughness and soil moistures.

1. Introduction

According to the Global Fatal Landslide Database, Asia suffers the greatest impact of fatal landslides among all continents [1,2]. This fact has to do with the physiographical environment of this region, including active tectonics, frequent typhoons/tropical cyclones, as well as socioeconomic factors, such as rapid economic growth, human population increase, habitat expansion and even loose regulations [1,3]. Among all questions associated with landslides, where and when they occur and how big they are remain the first information people demand to know. From a rapid response perspective, landslide sizes and locations are the key information needed by the ground crews to ensure the safety of human lives and the transportation of aids and supplies. From a policy-making perspective, long-term spatiotemporal evolution of landslide hotspots impacts the formulation of management strategies and mitigation plans [4]. From a scientific perspective, landslide volumes and frequencies may provide information about rock strength and denudation rates [5,6], the influence of climate changes [2], and the interaction among the lithosphere, hydrosphere and biosphere [7,8].

Given the large-area imaging capability from the sky, remote sensing has been the most widely-used tool in landslide detection during the past decades. For example, based on Formosat-2 satellite images, Lin et al. (2017) found that about 70% of the mountainous area in Taiwan has experienced at least 1 landslide within the decade between 2003 and 2012 [4]. In Japan, the National Research Institute for Earth Science and Disaster Prevention (NIED) also conducts regular landslide mappings based on aerial photos (https://www.bosai.go.jp/e/research/database/earth-and-sand.html, accessed on 10 March 2022; a study using the national landslide database can be found in [8]). The optical image-based landslide mapping has high quality and accuracy, although the availability usually depends on the weather condition and the source of sun lights after a landslide occurs. In some cases, the latency between the event and the first usable image can be weeks. This is where the synthetic aperture radar (SAR)-based mapping may serve as an alternative solution particularly during the early phase of the hazards. Some recent examples include the earthquake-triggered landslides after the 1999 Chi-Chi earthquake in Taiwan [9], the 2015 Gorkha earthquake in Nepal [10,11], the 2015 Mt. Kinabalu earthquake in Malaysia [12], the 2018 Lombok earthquakes in Indonesia [10], and the 2018 Hokkaido Eastern Iburi earthquake in Japan [10,13,14,15,16,17]. Examples of rainfall-triggered landslides include the 2009 Typhoon Morakot in Taiwan [18], the 2011 Typhoon Talas in Honshu, Japan [19], the 2015 heavy rain in Chin State, Myanmar [20], and the 2017 heavy rain in Kyushu, Japan [16].

In general, SAR-based change detection can be classified into two categories–the coherent change detection (CCD) and incoherent change detection (ICD)–depending on whether interferometric phase is used [15]. Speaking of change detection for landslides, CCD includes the pre-failure landslide monitoring by using interferometric phase time-series, as well as the post-failure landslide detection by comparing the change of the interferometric phase qualities (interferometric coherence). In the post-failure landslide detection, a landslide patch is usually depicted by low coherence due to the changes in surface geometry, roughness and dielectric properties. It works particularly well in regions with intermediate-to-high pre-event coherence [10,12,15]. However, in places where phase decorrelation occurs constantly, CCD may fail to provide accurate landslide information. One such a place is the forest, where volumetric decorrelation and temporal decorrelation prevail due to frequent changes in vegetation structures and dielectric properties, as well as the disturbance by winds, rains, water vapors and other atmospheric conditions [21]. It has been shown that coherence-based landslide detection methods may yield unreliable results in the forest areas where the pre-event coherence is too low [10,15].

On the other hand, ICD compares SAR backscattering amplitudes or intensities before and after the landslide event. Similar to interferometric coherence, changes in intensity are also associated with variations in surface geometry, roughness and dielectric properties. Backscattering intensities are, however, less sensitive to atmospheric conditions, and they also appear in a wider range of values that allow a better separation of major changes (such as from vegetation to bare surface) from minor variations (such as tree growth). In addition, most ICD methods are relatively easy to implement as they only involve image-wise calculations (except for the intensity correlation method in [15,19,20]). ICD has therefore been adopted more widely for landslide detection than CCD, especially for forest areas [13,15,17,19,22].

Single-polarimetric (single-pol) SAR data is by far the most commonly-used datatype in ICD. Having said that, multi-polarimetric (multi-pol) datasets and their decomposition parameters can also be used in ICD. SAR polarimetry, or PolSAR, takes into account of both the intensities and phases from the same image epoch acquired at different transmitting-receiving polarizations (HH, HV, VV and VH). Polarimetric decomposition then recombines the complex scattering coefficients to extract parameters that can directly infer physical properties of the scatterers. These decomposition parameters, even from a single image epoch, have been proved to be efficient in differentiating land cover types including landslides [9,16,23,24]. Some studies also use PolSAR datasets and decomposition parameters in the dual-temporal (1 pre-event and 1 post-event image) change detection to improve the result accuracy and to gain physical insights [16,18,25]. In comparison, PolSAR’s detection capability in a multi-temporal context (≥2 pre-event images and 1 post-event image) has not received much attention yet and deserves more exploration.

This study attempts to incorporate the multi-temporal PolSAR decomposition parameters into SAR-based landslide detection in forest areas, and compares the performance with those from single-pol datasets. This comparison is out of a practical consideration: Given more and more public and private players in SAR space missions, we expect the list of SAR sensors to increase and the overpass latency to shorten substantially in the near future. Soon SAR-based landslide detection may be not so much limited by data availability, and hence our knowledge about how different data properties affect the detection performance will be pivotal. Such knowledge will help researchers and government agencies make decisions before choosing the best dataset for landslide mapping. It will also provide information about the uncertainties and limitations especially for a rapid-response product.

The goal of this study is therefore to evaluate the influence of different SAR data properties on landslide detection, including radar wavelengths, spatial resolutions, polarizations and viewing geometry. Different datatypes, or combinations of the properties listed above, are produced from some of the most commonly-used sensors (the L-band ALOS-2, C-band Sentinel-1 and X-band COSMO-SkyMed) to facilitate this performance study. To unify the comparison basis, we carry out change detection using a newly-designed algorithm Growing Split-Based Approach (GSBA). GSBA takes in a SAR-derived value, either a backscattering intensity or decomposition parameter, normalized over its time-series variations and computes the Bayesian probability of landslides given that value. We compare the detection performance in two landslide cases, the earthquake-triggered landslides due to the 2018 Hokkaido Eastern Iburi Earthquake in Japan, and the rainfall-triggered landslide caused by the 2017 heavy rain in the catchment of Putanpunas River, southern Taiwan. With results from these two cases, we discuss how different data properties affect the landslide detection, and what may limit the detection efficacy.

2. SAR Data Processing

To achieve our objectives, we need to produce multiple datatypes from the same SAR data (see

Table 1. List of SAR data used in this study.

Sensor & Track *1	Pre-Event Epochs	Post-Event Epoch	Average Look Angle ($\mathsf{\theta }$)	Mode and Resolution *2	Wavelength	Polarization *3
Hokkaido Landslides (Japan), 2018-09-06, earthquake-triggered
ALOS-2 A122	2018-08-25 2017-08-26 2016-08-27	2018-09-08	30°	High-Sensitive 6 m (HR)	L-band 22.9 cm	Full-pol HH, HV, VV, VH
S-1 A68	2018-09-01 2018-08-20 2018-08-08	2018-09-13	39°	Interferometric Wide 15 m (MR)	C-band 5.6 cm	Dual-pol VV, VH
CSK A	2018-06-04 2017-07-16	2018-09-08	37°	StripMap 3 m (UHR)	X-band 3.1 cm	Single-pol HH
Putanpunas Landslide (southern Taiwan), 2017-06-07, rainfall-triggered
ALOS-2 A137	2016-12-22 2016-08-18 2016-06-09 2016-04-14 2016-03-03	2017-08-03	33°	ScanSAR 60 m (LR)	L-band 22.9 cm	Dual-pol HH, HV (HV-mode is missing on the post-event epoch)
ALOS-2 D27	2017-05-21 2017-04-23 2017-01-01 2016-12-04 2016-10-09	2017-07-02	44°	ScanSAR 60 m (LR)	L-band 22.9 cm	Dual-pol HH, HV
S-1 A69	2017-05-27 2017-05-15 2017-05-03 2017-04-21 2017-04-09	2017-06-08	35°	Interferometric Wide 15 m (MR)	C-band 5.6 cm	Dual-pol VV, VH
S-1 D105	2017-05-29 2017-05-17 2017-05-05 2017-04-23 2017-04-11	2017-06-10	38°	Interferometric Wide 15 m (MR)	C-band 5.6 cm	Dual-pol VV,VH
CSK D	2017-06-01 2017-05-24 2017-05-08 2017-04-22 2017-04-14	2017-06-09	27°	StripMap 3 m (UHR)	X-band 3.1 cm	Single-pol HH

*¹ S-1 = Sentinel-1; CSK = COSMO-SkyMed; A = ascending; D = descending. *² UHR = ultra-high resolution; HR = high resolution; MR = medium resolution; LR = low resolution. *³ Full-pol = full-polarimetric; Dual-pol = dual-polarimetric; Single-pol = single-polarimetric. H = horizontally polarized; V = vertically polarized. The first letter in the combination stands for the polarization of the transmitted wave, and the second is for the received wave.

for the list of SAR data used in this study). Next we describe the processing flows for two major datatype categories and the generation of Z-score maps, which are the input to the GSBA algorithm.

2.1. Single-Polarization: Backscattering Coefficient (${\sigma }^0$)

Backscattering coefficient ${\sigma }^0$ is the normalized measure of the radar signal’s strength reflected by a distributed target. From the single-look complex (SLC) stack, the general processing steps to obtain ${\sigma }^0$ include radiometric calibration [26], speckle noise attenuation [27,28], multilooking and geocoding (Figure 1). For Sentinel-1 (S-1) data, thermal noise removal is carried out concurrently with radiometric calibration before converting the digital numbers to ${\sigma }^0$ [29]. We only process the co-polarized ${\sigma }^0$ stacks (HH or VV) for their higher sensitivity to surface-related scattering [30]. All data processing is carried out using the graph processing tool (gpt) in SeNtinel Application Platform (SNAP) and the InSAR Scientific Computing Environment (ISCE, for ALOS-2 ScanSAR data only) built in a high-performance computing cluster. Two auxiliary datasets, local incidence angle (LIA) and layover-shadow mask, are also generated during data processing [31]. LIA is the angle between the radar incidence direction and the slope normal vector. Small LIAs indicate slopes facing the satellite, while large LIAs indicate either slopes facing the satellite but significantly deviating from the line-of-sight (LOS) direction, or slopes facing away from the satellite.

2.2. Multi-Polarization: Degree of Polarization ($m_{DP}$) and Scattering Powers

For a dual-polarimetric SAR data, we follow the same steps as those in the ${\sigma }^0$ processing flow to generate geocoded complex scattering coefficients $S_{ij} \left(i, j=H\mathrm{or}V\right)$(Figure 1). We then construct the $2\times 2$ covariance matrix $C_{\mathbf{2}}$ [32],

$$C_{\mathbf{2}}=⟨k·k^{*T}⟩=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]=\left[\begin{array}{cc}⟨{\left|S_{ii}\right|}^2⟩& ⟨S_{ii}{S}_{ij}^{*}⟩\\ ⟨S_{ij}{S}_{ii}^{*}⟩& ⟨{\left|S_{ij}\right|}^2⟩\end{array}\right]$$

(1)

where $⟨·⟩$ indicates spatial ensemble averaging (using a 5 × 5 window in our case) and $k$ is the target vector,

$$k=\left[\begin{array}{c}{S}_{ii}\\ S_{ij}\end{array}\right]$$

(2)

From $C_{\mathbf{2}}$ we can derive the 2D Barakat degree of polarization $m_{\mathrm{DP}}$, which is defined as the ratio between the intensity of the polarized portion to that of the total intensity [33],

$$m_{\mathrm{DP}}=\sqrt{1-\frac{4\times \det \left(C_{\mathbf{2}}\right)}{{\left(\mathrm{Tr}\left(C_{\mathbf{2}}\right)\right)}^2}}$$

(3)

$m_{\mathrm{DP}}$ represents the anisotropy from polarization structures, whereas the scattering randomness, $1-m_{\mathrm{DP}}\beta$ ($\beta$ is a measure of the relative dominance of polarized scattering), is considered as the dual-pol radar vegetation index [34].

For the full-pol SAR data (ALOS-2 A122 in Table 1), we form the $3\times 3$ coherency matrix $T_{\mathbf{3}}$ [30]:

$$\begin{gathered} T_{\mathbf{3}}=⟨\underset{\_}{k}·{\underset{\_}{k}}^{*T}⟩=\left[\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right] \\ =\frac{1}{2}\left[\begin{array}{ccc}⟨{\left|S_{HH}+S_{VV}\right|}^2⟩& ⟨\left(S_{HH}+S_{VV}\right){\left(S_{HH}-S_{VV}\right)}^{*}⟩& 2⟨\left(S_{HH}+S_{VV}\right)S_{HV}^{*}⟩\\ ⟨\left(S_{HH}-S_{VV}\right){\left(S_{HH}+S_{VV}\right)}^{*}⟩& ⟨{\left|S_{HH}-S_{VV}\right|}^2⟩& 2⟨\left(S_{HH}-S_{VV}\right)S_{HV}^{*}⟩\\ 2⟨S_{HV}{\left(S_{HH}+S_{VV}\right)}^{*}⟩& 2⟨S_{HV}{\left(S_{HH}-S_{VV}\right)}^{*}⟩& 4⟨{\left|S_{HV}\right|}^2⟩\end{array}\right] \end{gathered}$$

(4)

where the target vector $k$ is defined as [30]

$$k=\frac{1}{\sqrt{2}}\left[\begin{array}{c}{S}_{HH}+S_{VV}\\ S_{HH}-S_{VV}\\ 2S_{HV}\end{array}\right]$$

(5)

From $T_{\mathbf{3}}$, we carry out a Model-Free 3-Component decomposition for Full-pol data (MF3CF) which jointly considers the Barakat degree of polarization and the received wave information to allow the estimation of scattering powers without any assumption of scattering models [35]. The odd-bounce surface scattering power $P_s$, even-bounce scattering power $P_d$ and the diffused (volumetric) scattering power $P_v$ can be estimated as [35]:

$$\left\{\begin{array}{c}{P}_s=\frac{m_{\mathrm{FP}}\mathrm{Span}}{2}\left(1+\sin 2{\theta }_{\mathrm{FP}}\right)\hfill \\ P_d=\frac{m_{\mathrm{FP}}\mathrm{Span}}{2}\left(1-\sin 2{\theta }_{\mathrm{FP}}\right)\hfill \\ P_v=\mathrm{Span}\left(1-m_{\mathrm{FP}}\right)\hfill \end{array}\right.$$

(6)

where $m_{\mathrm{FP}}$ is the 3D Barakat degree of polarization,

$$m_{\mathrm{FP}}=\sqrt{1-\frac{27\times \det \left(T_{\mathbf{3}}\right)}{{\left(\mathrm{Tr}\left(T_{\mathbf{3}}\right)\right)}^3}}$$

(7)

$\mathrm{Span}=T_{11}+T_{22}+T_{33}$ , and the scattering type parameter ${\theta }_{\mathrm{FP}}$ is defined as,

$$\tan {\theta }_{\mathrm{FP}}=\frac{m_{\mathrm{FP}}\mathrm{Span}\left(T_{11}-T_{22}-T_{33}\right)}{T_{11}\left(T_{22}+T_{33}\right)+m_{\mathrm{FP}}^2{\mathrm{Span}}^2}$$

(8)

The calculation of 2D degree of polarization and scattering powers is carried out in PolSAR-tools available at https://github.com/Narayana-Rao/PolSAR-tools (accessed on 10 March 2022).

2.3. Generating Z-Score Maps

To standardize the change detection flow for various data values (${\sigma }^0$, $m_{\mathrm{FP}}$ and scattering powers), we adopt a dimensionless Z-score (Z) of the post-event value normalized by the statistical information obtained from its pre-event time-series. It is calculated as [36],

$$Z=\frac{y_{post}-{\overline{y}}_{pre}}{{\sigma }_{pre}}$$

(9)

where $y_{post}$ is the data value on the post-event epoch, and [${\overline{y}}_{pre}$, ${\sigma }_{pre}$] are the mean and standard deviation calculated from the pre-event time-series. The Z-score value represents the difference between a pixel’s change value and its background mean value in the unit of background standard deviations. The background mean and standard deviation are derived from multiple pre-event epochs, and therefore may avoid the issues associated with a single biased pre-event image, a consideration commonly seen in a dual-temporal change detection scheme [37]. It also allows a more robust detection of minor changes if the pre-event time-series contains relatively stable values [36].

Among different datatypes, however, we caution that the Z-score map for CSK data in the Hokkaido case (Table 1) may not be as accurate due to the insufficient number of pre-event images–no acquisition exists before 16 July 2017. The temporal standard deviations thus derived tend to be too large and yield Z-score values smaller than those generated from other datasets. To work around this problem, we estimate the spatial standard variation within a 21 × 21 window centered at each pixel, and use it as ${\sigma }_{pre}$ in Equation (9) when the value is smaller than the temporal standard deviation. The Z-score values thus generated are visually more comparable to other Z-score maps. The detection results, however, may still be less accurate because of incomplete pre-event information.

In addition to Z-score maps generated from the aforementioned datasets, we further generate a Z-score map that combines the positive Z-score values in $P_s$ ($Z_{{\mathrm{P}}_s}$) and negative values in $P_v$ ($Z_{{\mathrm{P}}_v}$),

$$Z_{{\mathrm{P}}_c}=\left\{\begin{array}{c}{Z}_{{\mathrm{P}}_v} \mathrm{if} Z_{{\mathrm{P}}_v}<0 \mathrm{and} \left|Z_{{\mathrm{P}}_v}\right|>\left|Z_{{\mathrm{P}}_s}\right|\hfill \\ Z_{{\mathrm{P}}_s} \mathrm{otherwise}\hfill \end{array}\right.$$

(10)

where $Z_{{\mathrm{P}}_c}$ stands for the combined Z-score map. The necessity and performance of such a combined Z-score map will be further demonstrated in Section 4. Next, we describe how these Z-score maps are used in the GSBA change-detection algorithm.

3. Change Detection Method

The method proposed in this study is a variation of the split-based approach (SBA), which was first proposed for SAR-based flood mapping [38]. SBA is designed for the generalization of a mapping algorithm regardless of the image’s spatial resolution, swath size and the target’s geospatial distribution. The basic idea is to separate the image into multiple, non-overlapping tiles (also called splits). The tiles are then checked one by one to identify the existence of a certain proxies that signify the changes. Some suggested proxies include standard deviation [38], coefficient of variation [39,40], and the ratio between the tile mean and the global mean [40]. Tiles with proxies above a given threshold are selected to jointly determine the global threshold either via a non-parametric approach such as the Otsu method [41] or the KI algorithm [42], or through a parameterized fitting to the data histogram in order to determine the best cutoff point [43,44].

One variation of SBA is the hierarchical split-based approach (HSBA) [45]. Different from the conventional SBA which adopts a uniform tile size [38,40,46,47], HSBA adaptively and hierarchically splits the image to variable sizes until a bimodal histogram (representing the change and non-change classes) can be identified in the tile or when a minimal tile size is reached. This way it avoids the need of a pre-defined tile size, which in some cases may compromise the detection if the change area within a tile is extremely large or small. The issue with HSBA is that the splitting process is nonlinear-every split depends on the result of the previous split, and hence the computation time can be long when the image is large.

Instead of the top-down splitting strategy of HSBA, here we propose a bottom-up approach called Growing Split-Based Approach (GSBA) (Figure 2). We initialize the image splitting as in the conventional SBA. Once the tiles with changes are identified, we “grow” a patch within each tile cluster until the maximum patch area with a consistent bimodal histogram is reached. This growing step produces patches of different areas that mimic the variable tile size obtained by HSBA. The second variation in GSBA is that rather than obtaining a global threshold from the patches to generate a binary map, we calculate the Bayesian probability of changes instead [48]. Given the probability map, we can obtain a binary change map at a given cutoff probability (by default 0.5).

The last variation in GSBA is that instead of using a single tile size, the flow described above is repeated at different tile sizes. This practice is to acknowledge the observation that depending on the size and spatial distribution of changes, a single arbitrary tile size may in some cases produce results that fall into a local minimum or maximum of change area. With multiple binary maps generated at different tile sizes, we can select the one with intermediate spatial clustering (using Ripley’s K, see step (g) in Appendix A) as the final output. In other cases where different tile sizes produce similar binary maps, this practice also offers reassurance regarding the robustness of the detection output. The processing flow in GSBA is linear and can be fully parallelized, so the increase in computation time can be minor in a multi-processor computing system.

To avoid distraction from the main focus of this paper, details of the GSBA algorithm are given in the Appendix. The output Bayesian probability map is used in the following Receiver Operating Characteristics (ROC) curve analysis. To generate the ROC curves, we calculate the true positive rate (TPR, or recall) and false positive rate (FPR) at different cutoff probabilities. They are defined as follows:

$$\left\{\begin{array}{c}\mathrm{TPR}=\mathrm{TP}/\left(\mathrm{TP}+\mathrm{FN}\right)\hfill \\ \mathrm{FPR}=\mathrm{FP}/\left(\mathrm{FP}+\mathrm{TN}\right)\hfill \end{array}\right.$$

(11)

where [T, F] stand for true and false, [P, N] stand for positive and negative, and any two-letter combination means the number of such pixels identified through validation. We further analyze the area under curve (AUC) to evaluate the overall detection performance. This continuous tracking of trade-off effects between TPR and FPR allows users to have a complete and visual overview of the detection performance [15]. To compare with other studies using the same landslide cases, we also report the overall accuracy (OA) by using the final binary maps from GSBA. The overall accuracy is defined as:

OA = (TP + TN)/(TP + TN + FP + FN)

(12)

Note that these metrics are calculated at each datatype’s spatial resolutions, which means, the validation dataset is resampled to match the resolution of the SAR images.

4. Results

4.1. Case Study 1: Earthquake-Triggered Hokkaido Landslides in Japan

Widely-distributed landslides occurred due to seismic shaking during the 6 September 2018 Mw 6.6 Hokkaido Eastern Iburi Earthquake (Figure 3). After the earthquake, the Geospatial Information Authority of Japan (GSI) acquired aerial photos on 6 and 11 September over the landslide areas and identified more than 3300 landslide patches manually (https://www.gsi.go.jp/BOUSAI/H30-hokkaidoiburi-east-earthquake-index.html#1, accessed on 10 March 2022). Visual comparison between the aerial photos taken on 6 and 11 September shows that only minor changes occurred between these two dates, and hence the majority of the landslides have existed since the earthquake.

We define a large area of interest (AOI) of 480 km² that covers most of the landslide patches (Figure 3). Within the large AOI we further define four test areas, three of which are the same as those chosen by Jung and Yun (2020) [15]. This choice allows us to compare the effect of different spatial resolutions-they adopt the ALOS-2 PALSAR-2 Ultra-Fine HH-pol SAR data of 3-m resolution, while the data used in this study is the 6-m High-Sensitive full-pol SAR. In addition, we define a fourth test area around the Atsuma Reservoir in order to examine the effect of water bodies during landslide detection.

4.1.1. Qualitative Comparison

From the three sensors used in this case study (Table 1), we produce the following seven SAR datatypes and one combined Z-score map to carry out the performance comparison:

(1)

High-res L-band HH-pol ${\sigma }^0$

(2)

High-res L-band VV-pol ${\sigma }^0$

(3)

Medium-res C-band VV-pol ${\sigma }^0$

(4)

Ultra-high-res X-band HH-pol ${\sigma }^0$

(5)

High-res L-band dual-pol (VV + VH) $m_{\mathrm{DP}}$

(6)

Medium-res C-band dual-pol (VV + VH) $m_{\mathrm{DP}}$

(7)

High-res L-band full-pol $P_s$

(8)

High-res L-band full-pol $Z_{{\mathrm{P}}_c}$ (denoted as $P_c$ datatype hereafter)

In Figure 4 we choose Test Area 2 to demonstrate the differences among some of the selected datatypes. Both the high-res L-band HH-pol ${\sigma }^0$ and the ultra-high-res X-band HH-pol ${\sigma }^0$ provide sharp outlines for landslides compared to the multi-pol datatypes ($m_{\mathrm{DP}}$ and $P_s$). This difference is caused by the spatial ensemble averaging carried out on the polarimetric datasets. The C-band VV-pol ${\sigma }^0$, acquired at a lower spatial resolution, does not capture the landslide boundaries as clearly. In addition, many pixels in the landslide patches contain low Z-score values, suggesting that these pixels are indistinguishable from non-landslides. The same phenomenon is also observed more in the X-band ${\sigma }^0$ than in the L-band ${\sigma }^0$, implying that shorter radar wavelengths may not perform as efficient in landslide detection within forest areas, possibly due to their higher sensitivity to small-scale changes in vegetation during the pre-event periods.

Next, we examine the effect of viewing geometry on different polarimetric combinations in the L-band datatypes. The $P_s$ datatype shows strong positive Z-scores for the landslides on the slopes with small LIAs. On the slopes with large LIAs, the landslide patches contain low to nearly zero $P_s$ Z-scores (pointed by yellow arrows in Figure 4). The same is also seen on the $m_{\mathrm{DP}}$ Z-score map. In comparison, the L-band ${\sigma }^0$ may still show clear and predominantly negative Z-score values on the slopes with large LIAs. We investigate the Z-score maps from other scattering powers and find that, instead of a positive increase in $P_s$ Z-score, landslide patches with larger LIAs tend to show a stronger decrease in $P_v$ Z-score (yellow arrows in Figure 5). This different dependency on LIA between $P_s$ and $P_v$ is also observed in other landslide cases [23]. Numerical simulation in [49] offers some possible explanation to this phenomenon, in which backscattering energy due to direct-ground reflection and crown-ground interaction decreases with LIA, while the energy due to direct-crown backscattering increases. In other words, landslides on the slopes with larger LIAs are sensed as “loses in tree crowns” instead of “increases in bare surfaces”. That means $P_s$ or $P_v$ each only carries half of the information about landslides. It is therefore necessary to consider both values when carrying out landslide detection, such as by combining them into a joint Z-score map $Z_{{\mathrm{P}}_c}$ (Equation (10) and Figure 5).

The aforementioned complementary effect between $P_s$ or $P_v$ is not seen between the dual-pol $m_{\mathrm{DP}}$ and the dual-pol radar vegetation index [34]. In fact, when comparing the Z-score maps generated from these dual-pol parameters, one appears as a sign-flipped image of the other–pixels in the landslide patches display similar amplitudes but opposite signs on the two Z-score maps. This is probably due to the incomplete information of scattering properties in dual polarization. We therefore do not attempt to combine these two parameters and remain with the $m_{\mathrm{DP}}$-only Z-score map. Next we will look into the quantitative comparison of detection performance among the eight listed datatypes.

4.1.2. Quantitative Comparison

Figure 6 shows the Bayesian probability maps calculated from different datatypes. Overall, the full-pol $P_c$ depicts the most complete shape of the landslide bodies, followed by the single-pol ${\sigma }^0$ at different radar wavelengths. The dual-pol datatypes detect only few landslides. The X-band ${\sigma }^0$ captures an additional high-probability patch in Test Area 1 with sharp outlines. This patch is possibly related to crops, and the false detection is likely associated with the lower number of pre-event images. In areas where water bodies exist (Test Area 4), the single-pol ${\sigma }^0$ may detect changes associated with both the landslides and the water bodies, while the multi-pol datatypes are relatively free of such confusions.

The performance of each datatype is shown in the ROC curves (Figure 7). In Test Area 1, the ROC curves for the L-band ${\sigma }^0$ and $P_c$ show a steep increase of TPR at low FPR, yielding AUC values as high as 0.83 to 0.86 (Figure 7f). These values are comparable to the detection results based on the ALOS-2 Ultra-Fine HH-pol multi-temporal intensities of 3-m spatial resolution (AUC = 0.79) [15]. The AUC for X-band ${\sigma }^0$ and L-band $P_s$ rank the second and third at 0.77 and 0.73. The L-band $m_{\mathrm{DP}}$, C-band $m_{\mathrm{DP}}$ and C-band ${\sigma }^0$ perform poorly, yielding AUC values lower than 0.7.

The landslide patches in Test Area 2, 3 and 4 result in slow-rising ROC curves and low AUC values in almost all datatypes (Figure 7b–d,f). Among them, the full-pol $P_c$ still gives the highest and more stable AUC values, above 0.7 for all three test areas. All the other datatypes yield AUC values lower than 0.7. The C-band ${\sigma }^0$ yields the lowest AUC in Test Area 2 and 3–lower than 0.6. In Test Area 4, the L-band ${\sigma }^0$ gives the lowest AUC values as a result of false detection over the reservoir lake (Figure 3 and Figure 6).

For the full area, the L-band full-pol $P_c$ has superior performance in landslide detection with AUC = 0.77 (Figure 7f). The comparison between $P_c$ and $P_s$-only AUC values again confirms that the combination of $P_s$ and $P_v$ is necessary for landslide detection. The full-area AUC value for the high-res L-band ${\sigma }^0$ and ultra-high-res X-band ${\sigma }^0$ ranks the second and third, higher than the medium-res C-band ${\sigma }^0$. We can therefore confirm that spatial resolutions plays a role as important as radar wavelengths in SAR-based landslide detection. The performance of the dual-pol $m_{\mathrm{DP}}$ is least favored because of its lower full-area AUC and its tendency to detect fewer landslides.

In Figure 7a–e, we also mark the positions of the binary change maps determined by GSBA on each ROC curve. These solutions are located close to the turning points of the curves, but slightly leaning towards the lower-FPR end. These positions suggest that the default cutoff probability of 0.5 tends to create conservative change maps with a higher positive likelihood ratio (TPR-to-FPR ratio).

4.2. Case Study 2: Rainfall-Triggered Putanpunas Landslide in Southern Taiwan

In the second case, we look into a different scenario–a rainfall-triggered landslide. During the first 4 days in June 2017, a total of nearly 900 mm of precipitation occurred in the upstream of the Laonong River (according to the records from weather station C0V210, available at https://e-service.cwb.gov.tw/HistoryDataQuery/index.jsp, accessed on 10 March 2022) (Figure 8). Three days later on 7 June 2017, a landslide-associated earthquake (landquake) was detected through the Real-time Landquake Monitoring System (RLAMS, http://collab.cv.nctu.edu.tw/older/catalog_20171231.html, accessed on 10 March 2022) [50]. Given the order of occurrence between the torrential rain and the landquake, we attribute this event to a rainfall-triggered landslide.

To determine the landslide area associated with the 7 June 2017 landquake, we use SPOT-6/7 images acquired on 18 April 2017 and 5 July 2017 to carry out manual detection. Within a 15-km radius (approximately the estimation error of landquake relocation) of the epicenter, we locate a new sliding patch in the existing landslide area of the Putanpunas River catchment. Actually, there have been repeating landslides in this catchment since typhoon Morakot in the year of 2009 [18,51], making this catchment one of the most actively evolving landslide area in southern Taiwan. Different from the multiple small-scale landslides triggered by the Hokkaido Eastern Iburi Earthquake, this rainfall-triggered Putanpunas landslide contains a single patch of a relatively large area, up to 400,000 m² (Figure 8). We caution that the actual landslide patch may differ from what we map here due to the longer latency between the event date and the post-event optical image.

4.2.1. Qualitative Comparison

In this case study, we produce the following five datatypes for comparison (note that different viewing geometry is involved):

(1)

Low-res L-band HH-pol ${\sigma }^0$, ascending

(2)

Low-res L-band HH-pol ${\sigma }^0$, descending

(3)

Medium-res C-band VV-pol ${\sigma }^0$, ascending

(4

Medium-res C-band VV-pol ${\sigma }^0$, descending

(5)

Ultra-high-res X-band, HH-pol ${\sigma }^0$, descending

Some datasets allow dual-pol polarimetric combinations, such as HH-HV for ALOS-2 track D27 and VV-VH for both Sentinel-1 tracks (Table 1). However, in the Hokkaido case we have demonstrated that the use of dual-pol $m_{\mathrm{DP}}$ datatype tends to capture fewer landslides, so we decide not to adopt the dual-pol datatype in this case study. Here we wish to focus on the detection performance among the ${\sigma }^0$ datatypes of different wavelengths, spatial resolutions and viewing geometry.

Figure 9 shows how different ${\sigma }^0$ datatypes compare visually. The ascending tracks show more prominent landslide signatures on the Z-score maps compared to those from the descending tracks, regardless of the wavelengths and spatial resolutions. The intriguing point is that the viewing geometry from ascending tracks yield larger LIAs compared to those from descending tracks. In addition, we find that the viewing geometry from descending tracks produces large layover zones that tend to overlap with the slopes of small LIAs. Among all sensors, the X-band CSK descending track produces the largest layover area, which also corresponds to the smallest look angle (27°) among all three sensors. Within the layover zones, landslide-related changes may still be recorded but are mixed with energies from multiple ground targets of the same range distance, leading to brighter pixels and stretched patterns after geocoding (Figure 10). We also notice that the current layover-shadow mask does not necessarily mask out all layover areas (Figure 10), which is possibly caused by errors in DEM or limitations in the SAR geometric distortion simulation [52].

We have to emphasize that the same layover-shadow masks are also calculated and applied on the SAR images in the previous Hokkaido landslide case. However, given the smaller slope angles (Figure 4), the layover and shadow effects are nearly negligible. In comparison, the Putanpunas River catchment has a large topographic relief and hence greater slope angles (Figure 9), resulting in a larger portion of unusable data within the layover zone.

4.2.2. Quantitative Comparison

Figure 11 shows the Bayesian probability, ROC curves and AUC values for different datatypes. The ascending tracks show better performance than descending tracks, with the highest AUC of 0.78 from the medium-res C-band VV ${\sigma }^0$. The ascending L-band HH ${\sigma }^0$, albeit its low spatial resolution (60 m), still yields an AUC value of 0.71. On the other hand, the descending tracks do not detect as many changes, with an AUC value of 0.65 for the medium-res C-band VV ${\sigma }^0$ and 0.64 for the ultra-high-res X-band HH ${\sigma }^0$. No change is detected from the descending L-band HH ${\sigma }^0$. We should point out that despite the similar AUC values from the descending C-band VV ${\sigma }^0$ and the X-band HH ${\sigma }^0$, their effective area (unmasked area) is different–more than 50% of the catchment is under the layover-shadow mask of the X-band data.

5. Discussion

5.1. How Data Properties Affect Detection Performance

In

Table 2. Summary of full-area AUC, OA, TPR at PRF = 0.1 and the ratio of effective area (A_e) for the two case studies.

Event and AOI Area	Sensor-Track	Datatype *1		AUC	OA	TPRFPR = 0.1 *2		Ae *3 Ratio
Hokkaido 480 km2	ALOS-2 A122	Single-pol HR L-band	HH ${\sigma }^0$	0.67	0.89	0.37	(0.66)	0.99
		Single-pol HR L-band	VV ${\sigma }^0$	0.68	0.89	0.39	(0.70)	0.99
		Dual-pol HR L-band	VV + VH $m_{\mathrm{DP}}$	0.65	0.89	0.42	(0.75)	0.99
		Full-pol HR L-band	$P_s$	0.66	0.90	0.40	(0.71)	0.99
		Full-pol HR L-band	$P_c$ Z-score	0.77	0.90	0.56	(1.00)	0.99
	S-1 A68	Single-pol MR C-band	VV ${\sigma }^0$	0.63	0.87	0.27	(0.48)	0.99
		Dual-pol MR C-band	VV + VH $m_{\mathrm{DP}}$	0.62	0.87	0.27	(0.48)	0.99
	CSK-A	Single-pol UHR X-band	HH ${\sigma }^0$	0.67	0.88	0.34	(0.61)	0.99
Putanpunas 20 km2	ALOS-2 A137	Single-pol LR L-band	HH ${\sigma }^0$	0.71	0.98	0.50		0.97
	ALOS-2 D27	Single-pol LR L-band	HH ${\sigma }^0$	-	-	-		0.97
	S-1 A69	Single-pol MR C-band	VV ${\sigma }^0$	0.78	0.92	0.61		0.87
	S-1 D105	Single-pol MR C-band	VV ${\sigma }^0$	0.65	0.90	0.41		0.81
	CSK-D	Single-pol UHR X-band	HH ${\sigma }^0$	0.64	0.98	0.37		0.54

*¹ UHR: ultra-high resolution; HR: high resolution; MR: medium resolution; LR: low resolution. *² The TPR value at FPR = 0.1. Values in parentheses represent the normalized value by the best performance. *³ The ratio between the effective area (area outside the layover-shadow mask) and the full AOI.

we summarize the performance of different datatypes, including the full-area AUC and OA. Note that AUC does not represent the performance of any particular detection outcome but a full spectrum of outcomes, whereas OA is calculated at a specific selection (e.g., at the cutoff probability of 0.5) to mimic the operator’s choice. The OA values among different datatypes in the Hokkaido case, however, are all very similar and cannot reflect their actual relative performance. This similarity results from the large number of non-landslide pixels (TN in Equation (12)) when computing the OA values. To allow a better judgement of the relative performance at a specific landslide mapping outcome, we compute the full-area TPR at a fixed FPR of 0.1 (TPR_{FPR = 0.1}). This value indicates the ratio of detectable landslides at the cost of 10% false positive detection. In the Hokkaido case, we further normalize the TPR_{FPR = 0.1} values by using the best TPR_{FPR = 0.1} from $P_c$ (Table 2). With the AUC and TPR_{FPR = 0.1} values, we discuss how the following data properties affect the performance of SAR-based landslide detection.

Radar wavelengths and spatial resolutions. In the Hokkaido case, the L-band datatypes have better AUC and TPR_{FPR = 0.1} values than the other two radar wavelengths. This result seems to suggest that longer wavelengths work better in landslide detection. However, spatial resolutions can be an equally important factor. This inference is made from the poor performance of the medium-res C-band single-pol datatype–its result is worse than that from the X-band single-pol data at a higher spatial resolution. At the same time, the medium-res C-band single-pol data seems to perform relatively well in the Putanpunas case. This better performance is probably associated with the geometry of the landslide patches (small and distributed landslides in Hokkaido vs. one single large patch in Putanpunas), and the fact that the Putanpunas landslide is a repeated landslide on a bare surface instead of on a forested land (Figure 8).

Polarizations. In the Hokkaido case, the high-res L-band full-pol data can offer the best landslide detection capability in forest areas. It can even avoid false detection over water bodies. Even though the landslide boundaries become slightly blurry compared to those detected by using single-pol datatypes, the amount of information contained in a full-pol dataset and the detection performance thereof is indeed unparalleled. The L-band dual-pol datatype, despite a slightly better TPR_{FPR = 0.1}, gives an AUC value lower than those from many single-pol datatypes. It also tends to detect fewer landslide patches. As dual polarization will be a major observation mode for some upcoming missions such as NISAR, more efforts are needed to explore a better utilization strategy for dual-pol data in landslide mapping.

Viewing Geometry. In the Hokkaido case, we do not see a notable difference in LIA between the detected and the missed pixels (Figure 12a). This similarity in LIA distributions suggests that large LIAs do not necessarily hinder landslide detection. In the Putanpunas case, landslides are even better detected on datatypes with larger LIAs (Figure 9). On the other hand, slopes with small LIAs may coincide with overlay zones, which are also seen in the Putanpunas case. The smaller (steeper) the satellite look angle (θ), the larger the overlay zone and the smaller the effective area. The percentage of effective area drops from 97% at θ = 44° (ALOS-2 D27) to 54% at θ = 27° (CSK-D) (Figure 9 and Table 2). As the area of layover and shadow depends on the specific topography of a region, a DEM-based SAR geometric distortion simulation should be executed before planning for a mission or submitting a tasking request in order to determine the best viewing geometry [52].

5.2. Limitations in SAR-Based Landslide Detection

In the Hokkaido case, the best full-area AUC is 0.77 and the best OA is 0.90 (Table 2). These values are close to the results produced by using the ultra-high-res L-band HH-pol multi-temporal intensity with different detection algorithms [15]. The similarity in limited performance from different studies suggests that there may exist some factors that prevent the SAR-based landslide mapping from achieving the high accuracy attained in optical image-based mapping [53]. We plot the histogram of the post-to-pre-event difference in scattering powers ($\mathrm{\Delta }\mathrm{Power}$) for the detected and the missed pixels, respectively, within the Hokkaido AOI (Figure 12b,c). While the histograms for the detected pixels are clearly skewed, the histograms for the missed pixels are centered and symmetric around zero, indicating no significant difference in any of the scattering powers before and after the landslide. This phenomenon is bewildering and needs some explanation.

The first possibility is the change of local topography. From a broader scale, there does not seem to be a systematic difference in LIA between the detected and the missed pixels (Figure 12a). However, the LIA is calculated using the pre-event global 30-m DEM [54]. The topography must have changed locally after the landslide. In fact, according to the post-event DEM obtained by airborne laser survey, the surface morphology within the landslide patches has changed substantially [55]. Features such as scarps and crown cracks can reach meters tall with very steep (nearly vertical) facets [55]. Field photos also reveal that large boulders, huge piles of dead trees and meter-scale surface undulations exist on the ground [55,56]. The chaotic distribution of these features may result in scattering properties considerably deviating from what is expected from a pure land cover-type change (forest to bare surface).

Another equally important factor to consider is the spatial variation of water contents. Several studies show that the water contents can vary remarkably in the Hokkaido area, from 30% to 280% among different geological materials sampled at the same site [57,58]. As these materials are spread out during the landslide, the randomness in surface soil moisture within the landslide patches increases. To sum up, soil moistures, local topography, surface roughness and land cover changes together form a complicated backscattering field in the landslide patches, leading to clear changes of scattering powers in some places and no clear changes in other places. We may even hypothesize that it is radar wave’s sensitivity to multiple physical properties that limits its performance in landslide detection. This hypothesis needs further validations though, such as through numerical simulation of 3D backscattering fields.

The increased randomness in backscattering fields within a landslide patch may explain why intensity correlation method can yield better detection results than those from pixel-by-pixel detection methods [15,19]. Intensity correlation is calculated based on a number of pixels within a moving window, which offers contextual information about the objects and their changes. Most important of all, it can potentially average out the randomness in the complicated scattering field within a landslide patch, leading to a smoother and less noisy detection result.

6. Conclusions

By applying a newly-designed change detection algorithm Growing Split-Based Approach (GSBA) on two different landslide cases, this study examines how different SAR data properties affect the performance of landslide detection. Our result shows that the high-resolution, full-polarimetric SAR datatype has unparalleled performance in landslide detection over forest areas. Single polarimetric datasets of high or ultra-high spatial resolution rank the next, regardless of their radar wavelengths. This result suggests that high spatial resolution is critical especially for detecting small and distributed landslides in forest areas. Datatypes of medium or low spatial resolution work better in detecting large landslide patches over bare surfaces; their detection performance decays significantly over small landslides in forest areas. Dual polarimetric datatypes have the worst performance among all; a better utilization strategy may be needed for their use in landslide detection. Different viewing geometry mainly impacts the effective detection area by creating layover and shadow zones. This problem is more severe in areas of large slopes (≥40–50°), in which a steep viewing angle (<30°) may render half of the image unusable for landslide detection. SAR geometric distortion simulation is recommended before planning for a mission or submitting a tasking request for landslide mapping purposes. Given the limited performance of SAR-based landslide detection (both in this study and in previous studies) as compared to that from optical image-based landslide detection, we propose that other confounding factors, including but not limited to local topography, surface roughness and soil moistures, are all contributing to the randomness in the backscattering field and hinder the detection of land cover changes. Such limitations need to be properly acknowledged when adopting SAR-based landslide mapping for emergency responses or post-hazard assessments.