Land Subsidence Susceptibility Mapping in Jakarta Using Functional and Meta-Ensemble Machine Learning Algorithm Based on Time-Series InSAR Data

Abstract

Areas at risk of land subsidence in Jakarta can be identified using a land subsidence susceptibility map. This study evaluates the quality of a susceptibility map made using functional (logistic regression and multilayer perceptron) and meta-ensemble (AdaBoost and LogitBoost) machine learning algorithms based on a land subsidence inventory map generated using the Sentinel-1 synthetic aperture radar (SAR) dataset from 2017 to 2020. The land subsidence locations were assessed using the time-series interferometry synthetic aperture radar (InSAR) method based on the Stanford Method for Persistent Scatterers (StaMPS) algorithm. The mean vertical deformation maps from ascending and descending tracks were compared and showed a good correlation between displacement patterns. Persistent scatterer points with mean vertical deformation value were randomly divided into two datasets: 50% for training the susceptibility model and 50% for validating the model in terms of accuracy and reliability. Additionally, 14 land subsidence conditioning factors correlated with subsidence occurrence were used to generate land subsidence susceptibility maps from the four algorithms. The receiver operating characteristic (ROC) curve analysis showed that the AdaBoost algorithm has higher subsidence susceptibility prediction accuracy (81.1%) than the multilayer perceptron (80%), logistic regression (79.4%), and LogitBoost (79.1%) algorithms. The land subsidence susceptibility map can be used to mitigate disasters caused by land subsidence in Jakarta, and our method can be applied to other study areas.

1. Introduction

Several cities in Indonesia suffer from degradation at the ground level of buildings, known as land subsidence [1,2]. In Jakarta, this process has had a severe impact on urban infrastructure, leading to cracks in buildings, roads and damage to drainage systems [3]. These conditions are problematic because land subsidence may expand coastal flood areas due to sea-level rise [4]. Heavy monsoon rainfall [5] has caused frequent river flooding; if this occurs again, Jakarta could be submerged entirely underwater [1,4].

Recent studies of land subsidence in Jakarta have used various geodetic measurement methods, such as leveling surveys [3] and a global positioning system (GPS) surveys [6,7]. These studies indicated that excessive groundwater extraction is the leading cause of land subsidence and compaction to the vulnerable aquifer system [8]. This compaction may be exacerbated by natural consolidation since Jakarta’s landform mostly comprises young alluvial soils that cannot support the weight of human-made structures [9,10]. Therefore, it is essential to monitor land subsidence in Jakarta to predict further possible occurrences and mitigate damage [11,12].

Over the last decade, land subsidence susceptibility maps have been generated using geological, geomorphological, topographical, and hydrological data; these are considered the main factors influencing land subsidence [11,13,14]. Various methods are used to generate land subsidence susceptibility maps, including frequency ratios (FR) [12,15], weight of evidence (WOE) [16,17], logistic regression (LR) [18], evidential belief functions [11], analytical hierarchy processes (AHP) [19], support vector machines (SVM) [14,20], decision trees [21,22], fuzzy logic [12,23], adaptive neuro-fuzzy inference systems (ANFIS) [24,25], and artificial neural networks (ANN) [26,27]. In general, using a single modeling method leads to lower predictive accuracy than an ensemble method that uses a combination of models and machine learning algorithms [13,28]. Machine learning algorithms have the advantage of finding unpredictable relationships in datasets at multiple scales and have, thus, been recommended to obtain accurate land subsidence susceptibility maps [22].

Another challenge in generating land subsidence susceptibility maps is the low availability of land subsidence inventory maps. In this study, a land subsidence inventory map was generated via time-series interferometry synthetic aperture radar (InSAR) analysis. This technique can be applied to measure the displacement of the earth’s surface with an accuracy of up to millimeters per year by improving the selection of coherent pixels and reducing atmospheric noise [29,30,31,32,33]. It has been widely used to monitor land subsidence and generate land subsidence maps, for example, in large cities in Mexico [33,34], Kurdistan, and Iran [35], in Yuncheng Basin, China [36], and in coastal cities and areas such as Venice, Italy [37], New Orleans, United States [38], and Shanghai, China [39]. The InSAR algorithm has been successfully applied to earthquakes [40], volcanic activities [41,42], crustal deformation [31], landslides [43], manmade deformations [44], excessive groundwater extraction [3], mining activities [45], and natural consolidation of young alluvial soil [9].

The recent studies of monitoring land subsidence in Jakarta using interferometry synthetic aperture radar (InSAR) techniques was conducted using the Small Baseline Subset (SBAS) algorithm from 2007 to 2009 [1] and using Geodesy and Earth Observing Systems-Persistent Scatterer Interferometry (GEOS-PSI) algorithm from 2007 to 2010 [46]. Both studies utilized Advanced Land Observing Satellite (ALOS) phased array type-L synthetic aperture radar (PALSAR) data to produce a land subsidence map, and both compared their results with GPS survey data. However, Jakarta’s land subsidence map requires updating as it has remained unchanged for the past 10 years, and research on land subsidence susceptibility maps was not found in Jakarta.

Therefore, this study’s objective was to generate the updated land subsidence map in Jakarta using the Stanford Method for Persistent Scatterers (StaMPS) with the Sentinel-1 synthetic aperture radar (SAR) dataset from March 2017 to May 2020 in ascending track and March 2017 to April 2020 in a descending track. The mean vertical deformation maps in ascending and descending tracks were compared to validate Jakarta’s land subsidence location. After that, the land subsidence map obtained with this method was used as an inventory map to generate a land subsidence susceptibility map in Jakarta that predicted the areas at risk of land subsidence in the future. Two meta-ensemble machine learning algorithms (adaptive boosting (AdaBoost) and Logit Boost), and two functional machine learning algorithms (logistic regression and multilayer perceptron) were used. The result of the land subsidence susceptibility map produced by these algorithms was evaluated to compare all models’ accuracy and reliability using receiver operating characteristic (ROC) curve analysis.

2. Study Area

Jakarta is Indonesia’s capital city (Figure 1a), located at 6°17’ south (S) and 106°82’ east (E), on the northern coast of western Java. It is considered a lowland area, with an average altitude of ±7 m above sea level [3,47]. In 2019, its population was 10.5 million, with a population growth rate of about 1.19% per year. The population density was 15,900 people per km², within a total land area of 662.33 km² [47]. Historically, the population of Jakarta gradually migrated to the other districts and municipalities that are part of the Jakarta Metropolitan Region (JMR), which covers a total area of 5897 km² and includes Bogor, Depok, Tangerang, and Bekasi [48] (Figure 1b). This migration to Jakarta’s outskirts led to increased urban development that could cause land subsidence in these areas [3]. Our study area covers the JMR; however, for simplicity, we refer to it as Jakarta. The three administrative areas were chosen as study area due to land subsidence reported on those areas from the recent study using GPS surveys, leveling surveys, and InSAR techniques. Training and testing data in Figure 1b used in this study to generate a land subsidence susceptibility map were assessed from the land subsidence inventory map.

The geological and geomorphological area in Jakarta lithologically was dominated by alluvium landform (50.20%) and alluvium fans (19.66%), followed by Serpong form (7.03%) that dominated by fragmented pumice sandstones, some limestone, and andesite. A volcanic pyroclastic flow formed tuff Banten (6.77%) and upper Banten tuff (4.88%), Sandstone unit (4.05%), swamp deposit (1.75%), Cihoe form (1.61%), beach ridge deposit (1.53%), and old alluvium (0.74%). Subang form (0.64%) was dominated by layered claystone lithology with limestone and marl found locally, marine deposit (0.52%), and young volcanic rocks (0.22%); Bojongmanik form (0.22%) was dominated by alternating sandstones and claystone inserted by limestones and coastal deposit (0.13%); Parigi form (0.05%) was dominated by medium limestone, lake (0.01%), and sandstone tuff (0.01%). The domination of alluvium landform has a risk of land subsidence due to the compaction of natural consolidation worsened by a human-made structure [3,9] covering the Jakarta area’s land use. Land use in Jakarta was dominated by settlement area with 45.72% and followed by 39.79% of rice field, 5.94% of dryland agriculture, 4.75% of fish pond, 1.98% of shrub-mixed dryland farms, 0.84% of airport, 0.51% of swamp, 0.35% of estate crop plantation, 0.09% of barren land, 0.03% of secondary mangrove forest, and 0.004% of swamp shrub.

Research on land subsidence in Jakarta has been conducted since 1982 using leveling surveys. There were two main monitoring periods: (1) 1982–1991 and (2) 1991–1997. Land subsidence was found in the period 1982–1991 in three regions with the highest accumulated subsidence value compared to other regions, namely, two regions in the northwestern part (Cengkareng and Kalideres districts) and the third in the northeastern part of Jakarta (Kemayoran-Sunter district); accumulated subsidence was found in Kalideres district with up to 68.5 cm between 1982 and 1991 with an annual rate of subsidence around 8 cm/year. Accumulated land subsidence was found in the Cengkareng district, with up to 60 cm between 1982 and 1991, with an annual subsidence rate of around 7 cm/year. Accumulated land subsidence was found in the Kemayoran-Sunter district, with up to 70 cm with an annual subsidence rate of around 6 cm/year. The second period of monitoring land subsidence using a leveling survey between 1991 and 1997 highlighted one region with the highest accumulated subsidence value than other regions in the Kalideres district, with up to 154.1 cm with an annual subsidence rate of around 23 cm/year [7].

Following this research, land subsidence in Jakarta was monitored using GPS surveys from 1997 to 2005, and two regions affected by land subsidence were reported in two stations: (1) Kwitang district and (2) Pantai Mutiara district. The accumulated land subsidence was found in the Kwitang district, with up to 48 cm with an annual rate of subsidence rate of 5 cm/year; accumulated land subsidence in Pantai Mutiara district was around 50 cm with an annual rate of subsidence around 4.6 cm/year. [6]. Land subsidence in Jakarta measured using leveling and GPS surveys positively correlated with excessive groundwater extraction and sea-level rise [6,7,49].

Land subsidence was also reported using InSAR techniques based on ALOS PALSAR satellite data in 2007–2009 using the SBAS method, with land subsidence found in Pluit district with an annual subsidence rate of 21.6 cm/year. Land subsidence was found in Cengkareng district with an annual subsidence rate around 21.8 cm/year, in Bekasi district with an annual rate of subsidence of around 10.6 cm/year, and in Karawang district with an annual rate of around 16.4 cm/year [1]. Research was also conducted with ALOS PALSAR satellite data within the period 2007–2010 using the GEOS-PSI method, finding land subsidence in the coastal area and lowland area in northwestern Jakarta within Penjaringan and Cengkareng districts with an annual subsidence rate of up to 26 cm/year, with accumulated subsidence rates of up to 86.5 cm between 2007 and 2010. Land subsidence was observed in the Bekasi district with an annual subsidence rate of up to 11.5 cm/year [46].

3. Material and Methods

3.1. SAR Datasets

A land subsidence inventory map for generating the land subsidence susceptibility map in Jakarta was generated using Sentinel-1 SAR C-band data (5.5 cm wavelength) provided by the European Space Agency (ESA). The SAR data were acquired from March 2017 to May 2020 (91 datasets in the ascending track with path number 98 and frame number 1160, with vertical–vertical (VV) polarization) and in the period of March 2017 to April 2020 (89 datasets in the descending track with path number 47 and frame number 614, with vertical–vertical (VV) polarization). The ascending and descending datasets are listed in Table 1 and Table 2, and the reference dates with zero delta day and zero perpendicular baselines from the ascending track are shown on 15 October 2018, whereas those from the descending track are shown on 16 November 2018; both reference dates are shown in bold text in table number 45. A perpendicular baseline graph (generated from Table 1 and Table 2 and shown in Figure 2) was used to visualize the temporal baseline from the reference date.

3.2. Land Subsidence Conditioning Factors

In total, 14 land subsidence conditioning factors consisting of geological, geomorphological, topographical, and hydrological data were chosen as the conditioning factors (Table 3) in this study as they are considered the main factors influencing land subsidence [11,13,14]. Each factor was classified using the quantile method, and the factors with different cell size resolutions were resampled into raster datasets with 30 m cell size from each conditioning factor to standardize each factor’s resolution.

A recent study in Jakarta found that the leading causes of land subsidence in Jakarta are groundwater extraction, the load of buildings and construction, natural consolidation of alluvium soil, and tectonic activities [3]. In this study, groundwater drawdown data were collected from the Groundwater Conservation Center, The Ministry of Energy and Mineral Resources of Indonesia. The observation of groundwater drawdown data is carried out periodically through an automatic water level record (AWLR) system. The annual change in groundwater level from 2019 to 2020 was calculated from 15 monitoring wells, and the map was constructed (Figure 3a) using the inverse distance weighting (IDW) interpolation method from the annual change in groundwater level data. The obtained groundwater drawdown data in this study were insufficient compared to the study area due to the limited monitoring wells over the study area. Although the available data are few, the use of groundwater drawdown data is essential to determine the relationship between land subsidence and groundwater extraction.

The groundwater level can increase due to conditional factors indirectly associated with land subsidence, such as rainfall, distance from a river, and river density, which can recharge the groundwater level [13,50]. Four years of daily rainfall data (from 2017 to 2020) from seven weather stations were acquired from the Indonesian Meteorology, Climatology, and Geophysical Agency of Indonesia, and the annual rainfall intensity was calculated and interpolated using the IDW tool in geographic information system (GIS) software (ArcGIS; ESRI, Redlands, CA, USA) (Figure 3b). However, the availability of rainfall data was limited compared to the study area due to the cost and constraints of acquiring data. Nevertheless, the utilization of the rainfall data in this study is important to correlate the topographical factor related to the infiltration flow affecting the soil’s strength.

Road, river, and fault networks were acquired from the Geospatial Information Agency of Indonesia from the main road map, main river map, and fault map at a scale of 1:250,000 of polygonal-shaped data using the Atlas map of Indonesia. The information on the location of roads, rivers, and faults was used to construct a distance map (Figure 3c, Figure 3d, and Figure 3e, respectively) using the Euclidean distance tool, and the maps were classified using the quantile method to provide suitable classes within a 30 m cell size. The drainage density or river map density was estimated using the kernel density tool (Figure 3f). Distance to the road and land use are related to urban development in Jakarta, which can affect land subsidence [3]. The land-use map (Figure 3g) was acquired from the Ministry of Environment and Forestry of Indonesia that used Landsat data to generate a land cover map. The geological parameters describe the spatial correlation of lithological landform and land subsidence in Jakarta as being caused by the compaction of the alluvium soil landform [3]. The lithology map (Figure 3h) was acquired from the Geological Atlas map of Indonesia from the Ministry of Energy and Mineral Resources, and the polygonal-shaped map was converted into a raster map with 30 m cell size using GIS tools.

The topographical map, which includes elevation, slope, aspect, profile curvature, plan curvature, and topographic wetness index (TWI) (Figure 3i, Figure 3j, Figure 3k, Figure 3l, Figure 3m, and Figure 3n, respectively) data extracted from the digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) [51], was constructed using the basic terrain analysis tools. Elevation influences the hydrological properties and soil moisture, whereby a lower-elevation area possibly gains more precipitation than a higher-elevation area [13,14]. The slope is associated with land subsidence because it affects the infiltration of rainfall (a steeper surface slope decreases infiltration) [13,21], and the aspect is the second derivative of the slope that has a relationship with land subsidence because the slope aspect affects the strength of the soil due to the moisture preservation and the amount of vegetation [11]. Profile curvature is associated with flow speed, sediment, and erosion quantity, while plan curvature is perpendicular to the slope and indirectly affects land subsidence by influencing convergence and divergence of flow across the surface [13]. The TWI defines the degree of deposition of water at a specific site [52]. The topographical factors extracted from the SRTM DEM are widely used as conditioning factors in land subsidence susceptibility mapping [11,12,13,28].

The relationship of land subsidence occurrence with the conditioning factors in each class is described in Table 4. A ratio greater than 1 denotes that the class in the conditioning factor is more correlated with the land subsidence occurrence [11]. The calculation of frequency ratio shows that the land subsidence occurred in areas with groundwater level data between 20.27 and 21.00 and between 23.31 and 34.55 m below ground level, while land subsidence also occurred in areas with more annual rainfall intensity due to the recharge of groundwater level and the usage of groundwater. Areas between 0 and 126 m with roads were more correlated with land subsidence occurrence. The land subsidence areas correlated with fault distance were between 15,944 and 68,975 m from the fault location. There were three drainage density classes correlated with land subsidence occurrence and settlement areas correlated with land subsidence. In terms of lithological factors, alluvium, alluvium fans, beach ridge deposits, and sandstone landforms were considered more correlated with the land subsidence occurrences. There were three classes in the elevation map, four classes in the aspect map, two classes in the slope map, one class in the plan curvature, one class in the profile curvature, and three classes in the topographic wetness index map correlated with land subsidence occurrences.

3.3. Illustration of Methodology

The methods to generate land subsidence susceptibility maps using functional and meta-ensemble algorithms are described below and illustrated in Figure 4.

• Land subsidence occurrences were identified by exploiting Sentinel-1 SAR datasets from 2017 to 2020 from both ascending and descending tracks using time-series InSAR techniques based on the StaMPS algorithm. The persistent scatterer points from co-registered single master images showing a deformation value were used as the land subsidence inventory map.

• Preparation of training and testing datasets was conducted by randomly dividing the persistent scatterer (PS) points of time-series InSAR showing a vertical deformation into 50% training data to generate land subsidence susceptibility models and 50% testing data to validate the land subsidence susceptibility map, as done in other studies finding optimal results [28,53]. The distribution of training and test data is shown in Figure 1b.

• Preparation of land subsidence conditioning factors for spatial correlation analysis was done using the frequency ratio method to find the correlation between each factor and land subsidence occurrence [53]. We used each model’s ratio value and then used as the conditioning factors related to land subsidence occurrences. First, the conditioning factors were classified using quantile methods in GIS tools with a similar environment of 30 m cell size for each factor; then, the number of subsidence occurrences in each class was calculated using the cross-tabulation tool in GIS. Next, we calculated the ratio between the percentage of pixels of each conditioning factor class and the percentage of subsidence occurrence pixels to obtain the FR value as follows:

(1) $\mathrm{Frequency} \mathrm{Ratio}=\frac{\% \mathrm{of} \mathrm{class} \mathrm{conditioning} \mathrm{factor}}{\% \mathrm{of} \mathrm{land} \mathrm{subsidence} \mathrm{occurrence}}$

• 4.. The conditioning factors consisting of frequency ratio values were used to generate land subsidence susceptibility models using two functional algorithms (logistic regression and multilayer perceptron) and two meta-ensemble algorithms (AdaBoost and LogitBoost).

• 5.. After all land subsidence susceptibility maps were generated, all maps were validated using the test data prepared before and analyzed using ROC curve analysis.

3.4. StaMPS Processing

StaMPS (Stanford Method for Persistent Scatterers) is an analysis method used to facilitate the generation of time-series deformation images of all terrains, including nonurban areas. The StaMPS algorithm uses the spatial correlation of phase measurements rather than a functional temporal model to identify PS pixels. StaMPS processing does not use a model to describe how the displacement rate varies with time. To identify PS pixels in a single master stack of interferograms, StaMPS uses the phase characteristics from the dominant point scatterer in each area and creates interferograms from SAR images. It also reduces decorrelation [54]. Thus, the StaMPS algorithm can identify and extract the deformation signal from stable pixels in all terrains.

The slave images in the acquired SAR datasets were resampled to perfectly match the master image through a co-registration process before generating an interferogram. The co-registration process was applied to two different images to produce refined SAR images, which were then cropped to focus only on the area of interest. Next, differential InSAR (DInSAR) images were generated by subtracting the topographic InSAR images generated using the interferograms from the topographic phases of the SRTM DEM [55,56]. We used the PS method to measure the displacement of the earth’s surface [57] to generate the time-series deformation map. The main processes generated a single master stack of interferograms and topographic phase removal [33].

The StaMPS algorithm is shown in Figure 5. Multiple images were co-registered to generate a single master image for the ascending track on 15 October 2018 and the descending track on 16 November 2018; co-registered images with the topographic phases removed were subjected to amplitude and phase noise analysis to derive a subset comprising all PS pixels, with weeding performed using a threshold value of 0.4. The wrapped phase of the selected pixels was corrected for the spatially uncorrelated look-angle error. The phase was then unwrapped, and PS outputs were generated. The deformation map from the line of sight (LOS) displacement could be converted into vertical deformation data [58,59] by assuming the horizontal deformation as very small compared to the vertical deformation caused by land subsidence [60,61,62]. Recent studies using GPS and leveling surveys reported that the land subsidence in Jakarta shows a vertical deformation [7], and a vertical deformation pattern was also found in research using InSAR to monitor land subsidence in Jakarta [1,46]; hence, the deformation from the line of sight (LOS) in this study could be assumed as negligible and could be converted directly into vertical deformation value using Equation (2) by dividing the displacement or deformation from the line of sight (${\mathrm{d}}_{\mathrm{LOS}}$) by the cosine of incident angle (θ) from the radar signal. The results of vertical deformation were assigned a negative value from the initial ground-level observation point, indicating that the land subsidence occurred vertically at that point [9].

(2)$\mathrm{V}= \frac{{\mathrm{d}}_{\mathrm{LOS}}}{\cos \mathsf{\theta }}$

3.5. AdaBoost

AdaBoost is a machine learning algorithm introduced by Freund and Schapire (1997) [63]. AdaBoost’s classifier uses an adaptive resampling technique that produces a series of individual classifiers to classify training samples accurately. The frequency of variables selected by a weak learner was examined, and the relative importance of the variables could be determined [64]. AdaBoost combines multiple weak learners to derive a single strong learner by repeatedly calling a weak classifier and adjusting the attributed weight to the sample. The new classifier focuses more on the misclassified sample as the misclassified sample is increased compared to the correct sample. The final weight is obtained by adding or subtracting the updated weight in each iteration. The final model is obtained by dividing each final weight by the total adjusted weight [65]. The general form of the AdaBoost algorithm is as follows [66]:

• Start with weights ${\mathrm{w}}_{\mathrm{i}}=\frac{1}{\mathrm{N}}$ for $\mathrm{i} = 1, \dots , \mathrm{N};$

• Repeat this step for  $\mathrm{m} = 1, \dots , \mathrm{M}:$

  • Fit the classifier  ${\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right) \in \left\{-1,1\right\}$ using weights ${\mathrm{w}}_{\mathrm{i}}$ with the training data;

  • Compute  ${\mathrm{err}}_{\mathrm{m}}{= \mathrm{E}}_{\mathrm{w}}\left[1_{(\mathrm{y}\ne {\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right))}\right]{, \mathrm{c}}_{\mathrm{m}}=\log \left(\frac{{1 - \mathrm{err}}_{\mathrm{m}}}{{\mathrm{err}}_{\mathrm{m}}}\right);$

  • Set  ${\mathrm{w}}_{\mathrm{i}}\leftarrow {\mathrm{w}}_{\mathrm{i}} \exp \left[{\mathrm{c}}_{\mathrm{m}}{1}_{(\mathrm{y}\ne {\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right))}\right], \mathrm{i} = 1, 2, \dots , \mathrm{N}$, and renormalize so that ${{\displaystyle \sum }}_{\mathrm{i}}{\mathrm{w}}_{\mathrm{i}}= 1;$

• Output the classifier:  $\mathrm{sign}\left[{{\displaystyle \sum }}_{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{c}}_{\mathrm{m}}{\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right)\right]$.

3.6. LogitBoost

LogitBoost is a modified version of the AdaBoost algorithm, introduced by Friedman, Hastie, and Tibshirani [66], where the exponential loss function is replaced with the log-likelihood loss function. This method reduces classification error and is less sensitive to noise [28,67]. LogitBoost can handle multiple class problems and uses a regression scheme as the base learner for classification [64]. The general form of the LogitBoost algorithm is as follows [66]:

• Start with weights ${ \mathrm{w}}_{\mathrm{i}}=\frac{1}{\mathrm{N}}$ for $\mathrm{i} = 1, 2, \dots , \mathrm{N}, \mathrm{F}\left(\mathrm{x}\right) = 0$, and probability estimates $\mathrm{p}\left({\mathrm{x}}_{\mathrm{i}}\right) = \frac{1}{2}$;

• Repeat this step for $\mathrm{m} = 1, \dots , \mathrm{M}:$

  • Compute the working response and weights:

    (3)${\mathrm{l}}_{\mathrm{i}} =\frac{{\mathrm{y}}_{\mathrm{i}}^{\ast }-\mathrm{p}\left({\mathrm{x}}_{\mathrm{i}}\right)}{\mathrm{p}\left({\mathrm{x}}_{\mathrm{i}}\right) \left(1 - \mathrm{p}\left({\mathrm{x}}_{\mathrm{i}}\right)\right)},$

    (4)${\mathrm{w}}_{\mathrm{i}} = \mathrm{p}\left({\mathrm{x}}_{\mathrm{i}}\right) \left(1 - \mathrm{p}\left({\mathrm{x}}_{\mathrm{i}}\right)\right);$

  • Fit the function by weighted least-squares regression of ${\mathrm{l}}_{\mathrm{i}}$ to ${\mathrm{x}}_{\mathrm{i}}$ using weight ${\mathrm{w}}_{\mathrm{i}};$

  • Update the function as follows:

    (5)$\mathrm{f}\left(\mathrm{x}\right) \leftarrow \mathrm{f}\left(\mathrm{x}\right) + \frac{1}{2}{\mathrm{f}}_{\mathrm{m}} \left(\mathrm{x}\right),$

    (6)$\mathrm{p}\left(\mathrm{x}\right) \leftarrow \frac{{\mathrm{e}}^{\mathrm{f}\left(\mathrm{x}\right)}}{{\mathrm{e}}^{\mathrm{f}\left(\mathrm{x}\right)}{+\mathrm{e}}^{-\mathrm{f}\left(\mathrm{x}\right)}};$

• Output the classifier: $\mathrm{sign}\left[\mathrm{F}\left(\mathrm{x}\right)\right] = \mathrm{sign}\left[{{\displaystyle \sum }}_{\mathrm{m}=1}^{\mathrm{M}}{\mathrm{f}}_{\mathrm{m}}\left(\mathrm{x}\right)\right]$.

3.7. Logistic Regression

Logistic regression is a statistical method used to find the best model to describe the correlation between a dependent variable and several independent variables. This method’s advantage is that the variables do not need to be normally distributed [68]. It also offers several ways of selecting the best predictor for use in the P probability model [69,70]. The equations describing the logistic regression are as follows [28,69,70,71]:

(7)$\mathrm{f}\left(\mathrm{x}\right)=\mathrm{logit}\left(\mathrm{P}\right)=\ln \left[\frac{\mathrm{P}}{1-\mathrm{P}}\right]{=\mathrm{c}}_0{+\mathrm{c}}_1{\mathrm{x}}_1+\dots +{\mathrm{c}}_{\mathrm{n}}{\mathrm{x}}_{\mathrm{n}},$

(8)$\mathrm{P} = \frac{1}{{1 + \mathrm{e}}^{-\mathrm{f}\left(\mathrm{x}\right)}} = \frac{1}{{1 + \mathrm{e}}^{{-(\mathrm{c}}_0{+\mathrm{c}}_1{\mathrm{x}}_1{+\dots +\mathrm{c}}_{\mathrm{n}}{\mathrm{x}}_{\mathrm{n}})}},$

where $\mathrm{f}\left(\mathrm{x}\right)$ is a linear combination function called $\mathrm{logit}\left(\mathrm{P}\right)$, P is the probability of subsidence occurrence, $1-\mathrm{P}$ is the probability that subsidence will not occur, ${\mathrm{x}}_1{, \mathrm{x}}_2{,\dots , \mathrm{x}}_{\mathrm{n}}$ are input variables, ${\mathrm{c}}_0$ is the model intercept, and ${\mathrm{c}}_1{,\dots , \mathrm{c}}_{\mathrm{n}}$ are the approximated coefficients of regression.

3.8. Multilayer Perceptron

The multilayer perceptron is a machine learning algorithm based on the ANN technique that consists of input layers, hidden layers, and output layers [72]. The advantages of the multilayer perceptron algorithm are as follows: the distribution of the training data does not require any assumptions, most inputs are selected during the training process based on the weight adjustment, and the relative importance of the various input measures does not need to be determined [73]. The equation representing the multilayer perceptron for land subsidence classification is as follows [28,72,73]:

(9)$\mathrm{m} = \mathrm{f}\left(\mathrm{c}\right),$

where f(c) is a hidden function that is optimized during the training process for a given network architecture via the adjustable weights, and ${\mathrm{c} = \mathrm{c}}_{\mathrm{i}}$ for $\mathrm{i}=1,\dots ,14$, which is a vector containing 14 land subsidence conditioning factors.

4. Results

4.1. Land Subsidence Inventory Map

The mean vertical deformation maps for Jakarta shown in Figure 6a,b were derived by time-series InSAR analysis based on the StaMPS algorithm. We overlaid a true color red/green/blue (RGB) composite image from Sentinel-2 taken on 28 August 2020. Six areas from both ascending and descending tracks showed high deformation; thus, we plotted a deformation trend for these areas.

The vertical deformation rate in the ascending track at Point P1 (Figure 6c) represents the Kosambi area, which subsided 189.48 mm from 2017 to 2020. Point P2 (Figure 6c) represents the Cengkareng area, which subsided 184.02 mm from 2017 to 2020. Point P3 (Figure 6e) represents the Ciledug area, which subsided up to 155.22 mm from 2017 to 2020. Point P4 (Figure 6e) represents the Penjaringan area, which subsided 148.36 mm from 2017 to 2020. Point P5 (Figure 6g) represents the Bekasi area, which subsided 128.17 mm from 2017 to 2020. Point P6 (Figure 6g) represents the Cikarang area, which subsided 271.84 mm from 2017 to 2020.

The vertical deformation rate in the descending track at Point P1 (Figure 6d) represents the Kosambi area, which subsided 210.07 mm from 2017 to 2020. Point P2 (Figure 6d) represents the Cengkareng area, which subsided 216.19 mm from 2017 to 2020. Point P3 (Figure 6f) represents the Ciledug area, which subsided up to 155.92 mm from 2017 to 2020. Point P4 (Figure 6f) represents the Penjaringan area, which subsided 148.31 mm from 2017 to 2020. Point P5 (Figure 6h) represents the Bekasi area, which subsided 107.19 mm from 2017 to 2020. Point P6 (Figure 6h) represents the Cikarang area, which subsided 257.94 mm from 2017 to 2020.

Figure 6g,h show that the deformation pattern was mostly linear. However, Figure 6c–f, representing the Kosambi, Cengkareng, Ciledug, and Penjaringan areas, show quite periodic subsidence with the standard deviation of the vertical deformation rate being higher than the vertical deformation rate in Figure 6g,h. These results occurred due to the seasonal effect of groundwater extraction, and those areas were surrounded by a residential area that mostly used groundwater as the water source. Meanwhile, P6 is one of the most significant industrial areas in Indonesia. The mean vertical deformation maps from ascending and descending tracks were compared to validate the accuracy of the land subsidence inventory map using the StaMPS algorithm, and the result showed a good correlation in Figure 6i with a coefficient of correlation (R²) up to 0.9584 between ascending and descending tracks.

The persistent scatterer density of both tracks in the east and west areas was relatively low due to them being wetland areas and more vegetated than other areas. The SAR dataset used in this study from Sentinel-1 SAR C-band data with 5.5 cm wavelength could not deeply penetrate beneath the trees. Thus, to overcome that limitation, the interpolation of the persistent scatterer points from the descending track was constructed using kriging interpolation in GIS tools to provide land subsidence information over the study area as shown in Figure 6j [74].

4.2. Land Subsidence Susceptibility Map

The land subsidence susceptibility model’s performance depended on the calculated parameters (Table 5) for optimization used in this study.

A land subsidence susceptibility map was generated using 14 land subsidence conditioning factors, training data from our land subsidence inventory map, and four different algorithms: LogitBoost (Figure 7a), AdaBoost (Figure 7b), logistic regression (Figure 7c), and multilayer perceptron (Figure 7d). Land subsidence susceptibility indices were generated for all unique pixels in the study area. The susceptibility indices were reclassified using the quantile method to identify feature pairs in five susceptibility classes: very low, low, moderate, high, and very high [28,64].

The proportion of very high susceptibility land was quite similar for all four algorithms. However, the map generated by AdaBoost differed from those of the other three methods because the AdaBoost model could only classify susceptibility into four classes due to the limits of the probability range. Therefore, the very low susceptibility class was excluded. The susceptibility class proportions (pixel distributions) are shown in Figure 8. For the Adaboost algorithm, the proportions were 0%, 57.88%, 2.54%, 22.02%, and 17.56% for the very low, low, moderate, high, and very high classes, respectively; the respective values for the LogitBoost algorithm were 33.33%, 27.09%, 7.03%, 16.77%, and 15.78%, and those for the logistic regression algorithm were 32.62%, 13.86%, 16.44%, 18.51%, and 18.57%. Finally, the multilayer perceptron algorithm’s respective values were 40.64%, 18.70%, 14.02%, 13.42%, and 13.23%.

The distribution of pixels in the very high class and high class in Figure 8 showed similar results for each algorithm, resulting in the maps of these classes being quite similar. The very high class was considered to be the land subsidence areas shown in the mean vertical deformation map in Figure 6a,b, with similarity seen because of the training data used in this study being acquired from the land subsidence inventory map with a large spatial resolution, which allowed more effectively defining the land subsidence area in the susceptibility map. The moderate and high classes were considered land subsidence areas in the future, and the very low and low classes were areas with the lowest probability of land subsidence in the future. The difference between the moderate class from AdaBoost and that from other algorithms is that AdaBoost did not consider other areas far from the land subsidence occurrences. Other algorithms showed that residential areas located in alluvium landforms had a reasonable possibility of land subsidence.

4.3. Model Validation

The accuracy of all four algorithms in this study was evaluated by ROC curve analysis [12,17]. ROC curve analysis is a standard way of validating the probability models used to generate land subsidence susceptibility maps, according to the area under the curve (AUC) [22,28]. Higher values indicate more accurate and reliable models. If the AUC, which ranges from 0 to 1, is lower than 0.5, the model is considered unacceptably inaccurate [75].

Land subsidence susceptibility maps produced using functional (AdaBoost, LogitBoost) and meta-ensemble (logistic regression and multilayer perceptron) algorithms were compared. The ROC curves for the four algorithms are shown in Figure 9. The largest AUC of 0.811 was from the AdaBoost algorithm (blue line in Figure 9). The multilayer perceptron algorithm had the next largest AUC (0.800; purple line in Figure 9), followed by the logistic regression (0.794; green line in Figure 9) and LogitBoost algorithms (0.791; red line in Figure 9).

Variables having at least one tie between the positive actual state group and the negative actual state group can bias the results. Since all of the AUC values were higher than 0.5, the land subsidence susceptibility maps produced by all algorithms used in this study are acceptable for predicting high-risk subsidence areas in Jakarta [75]. However, the AdaBoost algorithm map had the best performance.

5. Discussion

5.1. Land Subsidence Inventory Map

A land subsidence inventory map was successfully created through time-series InSAR analysis of Sentinel-1 datasets, using the StaMPS algorithm. The map is displayed as a vertical deformation map in Figure 6a,b for descending and ascending tracks. The comparison was made to validate the mean vertical deformation’s accuracy used as a land subsidence inventory map. A coefficient of correlation of 0.9584 was found. As stated above, there were six high-deformation areas (Figure 6c–h) within ascending and descending tracks, which were equally affected (in terms of subsidence) by increased groundwater usage and existing buildings [3].

Figure 6g,h show areas with largely linear deformation rates, the standard deviations of which were smaller than those of the area shown in Figure 6c–f. The periodic subsidence shown was affected due to the variable climate in Jakarta, which is generally characterized by high rainfall but a dry summer. In the rainy season, the groundwater level of all aquifers beneath Jakarta increases, which manifests as deformation, i.e., an uplift in the ground level due to a change in the underlying material’s thickness [34,76,77,78,79].

As shown by the land subsidence occurrence (Figure 6a,b) and groundwater drawdown (Figure 3a) data, land subsidence for P1–P5, which are between 20.58 and 34.55 m below ground level, correlated with the groundwater level. The correlation between land subsidence and groundwater level might have been stronger if groundwater level data were available for the whole study area. The groundwater level is deepest for P1 and P2, ranging from 22.47 to 34.55 m. P6 has more significant deformation because it contains the largest industrial estate in Indonesia.

Jakarta is mostly situated on young alluvial soil (Figure 3h), which cannot tolerate the maximum compression force of many buildings [9]. Thus, compaction of the unconsolidated alluvial soil occurs and is exacerbated by groundwater extraction (because pore pressure is reduced, leading to further clay compaction) [80,81].

5.2. Land Subsidence Susceptibility Map

The method used to create land subsidence susceptibility maps strongly affects the quality of the mapping. Machine learning techniques are effective [12,13,28,64]. In particular, the method used to generate training and testing data is important. Accurate land subsidence inventory maps can be obtained using InSAR; we combined InSAR and GIS spatial data to produce an accurate land subsidence susceptibility map [12]. Nevertheless, the distribution data of groundwater drawdown and rainfall intensity data from the interpolation process might not represent the whole study area. They could have affected the spatial distribution of the raster map. Access to the monitoring wells and rainfall station outside the Jakarta area could provide more accurate analysis to determine the relationship between those factors and land subsidence.

The land subsidence susceptibility maps of all four algorithms used in this study could be ordered according to the accuracy and time taken to build the model due to the similarities in the conditioning factors, training and testing data, and study area [28]. We used ROC curve analysis to assess the accuracy of the maps. The AUC data showed that the AdaBoost algorithm had the highest susceptibility class predictive accuracy of 81.1%, which was 1.1% higher than the multilayer perceptron algorithm, 1.7% higher than the logistic regression algorithm, and 2% higher than the LogitBoost algorithm. Since the accuracy of the algorithms used was so closely related, we analyzed the time consumption of data preparation needed to build the model. The AdaBoost algorithm needed 218.63 s to build the model, the LogitBoost algorithm needed 174.94 s to build the model, logistic regression needed 16667.96 s or 4.63 h to build the model, and multilayer perceptron needed 47,538.65 s or 13.20 h to build the model, thus identifying the LogitBoost algorithm as the fastest.

However, the LogitBoost algorithm had the lowest prediction accuracy (79.10%). The Adaboost algorithm needed 218.63 s or 3.64 min to build the model, along with having the highest prediction accuracy of 81.10% in the ROC curve analysis. Logistic regression needed 4.63 h to build the model, along with a prediction accuracy of about 79.40%, and the multilayer perceptron needed 13.20 h, along with a prediction accuracy of 80%. The maps generated using LogitBoost, logistic regression, and multilayer perceptron (Figure 7b–d, respectively) showed similar weights and indicated that the settlement area was at a high risk of land subsidence. The land-use map generated in this study indicated a relatively strong correlation between land subsidence and urban development, likely due to excessive groundwater extraction in urban areas [3,6,7]. Our map could be used by local environmental authorities and those in charge of urban development to identify areas with high subsidence risk. The method of combining SAR and GIS spatial data to generate land subsidence susceptibility maps employed in this study could be applied to other regions.

6. Conclusions

We generated a land subsidence inventory map using a time-series InSAR method based on the StaMPS algorithm from Sentinel-1 SAR datasets in ascending and descending tracks. The comparison of both tracks was conducted, finding a coefficient of correlation between the two tracks of 0.9548. The inventory map could be used as a training and testing dataset to create a land subsidence susceptibility map for Jakarta using meta-ensemble (AdaBoost and LogitBoost) and functional (logistic regression and multilayer perceptron) machine learning algorithms. We created a land subsidence inventory map through time-series InSAR analysis of Sentinel-1 SAR datasets using the StaMPS algorithm. The land subsidence susceptibility map produced by the AdaBoost machine learning algorithm had higher accuracy (AUC = 0.811) and only need 3.64 min to build the model compared with the maps created using the other algorithms (multilayer perceptron, AUC = 0.800; logistic regression, AUC = 0.794; LogitBoost, AUC = 0.791). LogitBoost was the fastest algorithm in building the model but had the lowest predictive accuracy. Logistic regression and multilayer perceptron needed 4.63 and 13.20 h, respectively, to build the model. All of the maps showed acceptable accuracy, as the AUC values were all higher than 0.5; thus, they can all be used for analyzing land subsidence susceptibility in Jakarta. Our approach based on time-series InSAR analysis, machine learning, and GIS spatial data yielded reasonable predictions of areas with high risk of land subsidence. Further research could use alternative algorithms and conditioning factors to generate land subsidence susceptibility maps in other regions.

Author Contributions

Conceptualization, C.-W.L., and W.L.H.; methodology, C.-W.L., A.R.A., and W.L.H.; software, A.R.A. and W.L.H.; validation, C.-W.L. and W.L.H.; formal analysis, W.L.H.; investigation, C.-W.L.; resources, C.-W.L.; data curation, C.-W.L. and W.L.H.; writing—original draft preparation, W.L.H.; writing—review and editing, C.-W.L., A.R.A., and W.L.H.; visualization, C.-W.L., A.R.A., and W.L.H.; supervision, C.-W.L.; project administration, C.-W.L.; funding acquisition, C.-W.L. All authors read and agreed to the published version of the manuscript.

Funding

This research was supported by a grant from the National Research Foundation of Korea provided by the government of Korea (No. 2019R1A2C1085686).

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. (a) Study area location in Indonesia and (b) training and testing data sets from land subsidence location within Jakarta Metropolitan Region (JMR) depicted from Sentinel-2 satellite imagery taken on 28 August 2020, divided into three administrative areas: Jakarta, Tangerang, and Bekasi.

Figure 2. Perpendicular baseline graph from (a) ascending track and (b) descending track.

Figure 3. Land subsidence conditioning factors: (a) groundwater drawdown; (b) rainfall intensity; (c) distance to roads; (d) distance to rivers; (e) distance to faults; (f) drainage density; (g) land use; (h) lithology; (i) elevation; (j) slope; (k) aspect; (l) profile curvature; (m) plan curvature; (n) topographic wetness index.

Figure 3. Land subsidence conditioning factors: (a) groundwater drawdown; (b) rainfall intensity; (c) distance to roads; (d) distance to rivers; (e) distance to faults; (f) drainage density; (g) land use; (h) lithology; (i) elevation; (j) slope; (k) aspect; (l) profile curvature; (m) plan curvature; (n) topographic wetness index.

Figure 4. The study workflow.

Figure 5. Flowchart of Stanford Method for Persistent Scatterers (StaMPS) Processing.

Figure 6. Mean vertical deformation map for Jakarta depicted using Sentinel-2 image in (a) ascending tracks and (b) descending tracks; the vertical deformation rate at P1 (Kosambi) and P2 (Cengkareng) in (c) ascending and (d) descending tracks; the vertical deformation rate at P3 (Ciledug) and P4 (Penjaringan) in ascending (e) and descending (f) tracks; the vertical deformation at P5 (Bekasi) and P4 (Cikarang) in ascending (g) and descending tracks (h). (i) The comparison of mean vertical deformation between ascending and descending tracks. (j) Kriging interpolation of time-series deformation map from mean vertical deformation map of descending track, resulting in the extension of the land subsidence inventory map; the blue area is considered as the nonoccurrence area of land subsidence.

Figure 6. Mean vertical deformation map for Jakarta depicted using Sentinel-2 image in (a) ascending tracks and (b) descending tracks; the vertical deformation rate at P1 (Kosambi) and P2 (Cengkareng) in (c) ascending and (d) descending tracks; the vertical deformation rate at P3 (Ciledug) and P4 (Penjaringan) in ascending (e) and descending (f) tracks; the vertical deformation at P5 (Bekasi) and P4 (Cikarang) in ascending (g) and descending tracks (h). (i) The comparison of mean vertical deformation between ascending and descending tracks. (j) Kriging interpolation of time-series deformation map from mean vertical deformation map of descending track, resulting in the extension of the land subsidence inventory map; the blue area is considered as the nonoccurrence area of land subsidence.

Figure 7. Land subsidence susceptibility map for Jakarta, generated using four different algorithms: (a) AdaBoost, (b) LogitBoost, (c) logistic regression, and (d) multilayer perceptron.

Figure 8. Proportions of susceptibility classes for land subsidence susceptibility maps generated using four different machine learning algorithms.

Figure 9. Receiver operating characteristic (ROC) curves for the land subsidence susceptibility maps produced by functional (logistic regression and multilayer perceptron) and meta-ensemble (AdaBoost and LogitBoost) algorithms.

Table 1

Acquisition dates (format ddmmyyyy) of the Sentinel-1 synthetic aperture radar (SAR) datasets in ascending track. Delta days = number of days between each acquisition date. B⊥ = perpendicular baseline. The reference dates in ascending tracks are shown in bold text.

No.	Acquisition Date(ddmmyyyy)	DeltaDays	B⊥(m)	No.	AcquisitionDate(ddmmyyyy)	DeltaDays	B⊥(m)	No.	Acquisition Date(ddmmyyyy)	DeltaDays	B⊥(m)
1	18032017	−576	77	32	30042018	−168	64	62	07052019	204	93
2	30032017	−564	65	33	12052018	−156	−2	63	19052019	216	16
3	11042017	−552	3	34	24052018	−144	28	64	31052019	228	20
4	23042017	−540	20	35	05062018	−132	127	65	12062019	240	166
5	05052017	−528	−25	36	17062018	−120	103	66	06072019	264	104
6	17052017	−516	17	37	11072018	−96	95	67	18072019	276	44
7	29052017	−504	110	38	23072018	−84	62	68	30072019	288	90
8	10062017	−492	21	39	04082018	−72	98	69	11082019	300	−9
9	22062017	−480	21	40	16082018	−60	75	70	23082019	312	1
10	04072017	−468	112	41	28082018	−48	61	71	04092019	324	46
11	09082017	−432	54	42	09092018	−36	60	72	16092019	336	106
12	21082017	−420	91	43	21092018	−24	55	73	28092019	348	−14
13	02092017	−408	50	44	03102018	−12	115	74	10102019	360	−110
14	14092017	−396	22	45	15102018	0	0	75	22102019	372	−125
15	26092017	−384	43	46	27102018	12	53	76	03112019	384	19
16	08102017	−372	48	47	08112018	24	85	77	15112019	396	−2
17	20102017	−360	72	48	20112018	36	85	78	27112019	408	38
18	01112017	−348	43	49	02122018	48	85	79	09122019	420	98
19	13112017	−336	91	50	14122018	60	142	80	21122019	432	87
20	25112017	−324	30	51	26122018	72	6	81	02012020	444	92
21	07122017	−312	140	52	07012019	84	74	82	14012020	456	40
22	19122017	−300	60	53	19012019	96	46	83	26012020	468	22
23	31122017	−288	135	54	31012019	108	40	84	07022020	480	36
24	12012018	−276	81	55	12022019	120	103	85	19022020	492	74
25	24012018	−264	78	56	24022019	132	29	86	02032020	504	93
26	05022018	−252	37	57	08032019	144	−26	87	14032020	516	31
27	17022018	−240	15	58	20032019	156	26	88	26032020	528	41
28	01032018	−228	30	59	01042019	168	8	89	07042020	540	40
29	13032018	−216	101	60	13042019	180	75	90	19042020	552	14
30	06042018	−192	139	61	25042019	192	−31	91	01052020	564	67
31	18042018	−180	126

Table 2

Acquisition dates (format ddmmyyy) of the Sentinel-1 Sentinel-1 synthetic aperture radar SAR datasets in descending track. Delta days = number of days between each acquisition date. B⊥ = perpendicular baseline. The reference dates in descending tracks are shown in bold text.

No.	Acquisition Date(ddmmyyyy)	DeltaDays	B⊥(m)	No.	AcquisitionDate(ddmmyyyy)	DeltaDays	B⊥(m)	No.	Acquisition Date(ddmmyyyy)	DeltaDays	B⊥(m)
1	26032017	−600	22	31	02042018	−228	−20	61	27052019	192	24
2	07042017	−588	76	32	14042018	−216	40	62	08062019	204	0
3	19042017	−576	83	33	26042018	−204	0	63	20062019	216	−7
4	01052017	−564	60	34	08052018	−192	26	64	02072019	228	16
5	13052017	−552	0	35	20052018	−180	0	65	14072019	240	29
6	06062017	−528	0	36	13062018	−156	26	66	07082019	264	87
7	18062017	−516	48	37	25062018	−144	26	67	19082019	276	84
8	30062017	−504	0	38	31072018	−108	−42	68	31082019	288	−5
9	12072017	−492	0	39	24082018	−84	−1	69	12092019	300	0
10	24072017	−480	−10	40	17092018	−60	0	70	24092019	312	117
11	05082017	−468	41	41	29092018	−48	30	71	06102019	324	94
12	17082017	−456	0	42	11102018	−36	51	72	18102019	336	16
13	29082017	−444	57	43	23102018	−24	0	73	30102019	348	19
14	10092017	−432	0	44	04112018	−12	0	74	11112019	360	45
15	22092017	−420	−4	45	16112018	0	0	75	23112019	372	0
16	04102017	−408	0	46	28112018	12	−18	76	05122019	384	5
17	16102017	−396	16	47	10122018	24	3	77	17122019	396	54
18	28102017	−384	96	48	22122018	36	11	78	29122019	408	63
19	09112017	−372	0	49	03012019	48	0	79	10012020	420	35
20	21112017	−360	62	50	15012019	60	24	80	22012020	432	25
21	03122017	−348	−2	51	27012019	72	0	81	03022020	444	11
22	15122017	−336	−23	52	08022019	84	1	82	15022020	456	24
23	27122017	−324	−50	53	20022019	96	34	83	27022020	468	−8
24	08012018	−312	0	54	04032019	108	101	84	10032020	480	91
25	20012018	−300	0	55	16032019	120	23	85	22032020	492	72
26	01022018	−288	14	56	28032019	132	0	86	03042020	504	33
27	13022018	−276	39	57	09042019	144	−12	87	15042020	516	0
28	25022018	−264	0	58	21042019	156	133	88	27042020	528	0
29	09032018	−252	−36	59	03052019	168	82	89	09052020	540	0
30	21032018	−240	−68	60	15052019	180	42

Table 3

Land subsidence conditioning factors in the study area. DEM, digital elevation model; SRTM, Shuttle Radar Topography Mission.

Category	Factor	Source
Hydrological factors	Groundwater drawdown	Groundwater Conservation Center of Indonesia, The Ministry of Energy and Mineral Resources
Hydrological factors	Rainfall intensity	Meteorology, Climatology, and Geophysical Agency of Indonesia
Land cover factors	Road network	Geospatial Information Agency of Indonesia
Hydrological factors	River network	Geospatial Information Agency of Indonesia
Geological factors	Faults	Geospatial Information Agency of Indonesia
Land cover factors	Land use	The Ministry of Environment and Forestry of Indonesia
Geological factors	Lithology	The Ministry of Energy and Mineral Resources
Topographical factors	Elevation	DEM SRTM 1 Arc-Second Global
Topographical factors	Slope	DEM SRTM 1 Arc-Second Global
Topographical factors	Aspect	DEM SRTM 1 Arc-Second Global
Geomorphological factors	Profile curvature	DEM SRTM 1 Arc-Second Global
Geomorphological factors	Plan curvature	DEM SRTM 1 Arc-Second Global
Hydrological factors	Topographic wetness index	DEM SRTM 1 Arc-Second Global

Table 4

Relationship between land subsidence occurrence and conditioning factors using frequency ratio (FR) model.

No.	Conditioning Factor	Class/Category	Ratio each Class	Ratio of Occurrence	FR
1	Groundwater drawdown(m below ground level)	7.77–20.27	0.1831	0.1201	0.6561
		20.27–21.00	0.2021	0.2953	1.4616
		21.00–21.84	0.2361	0.2127	0.9009
		21.84–23.31	0.1914	0.1282	0.6697
		23.31–34.55	0.1874	0.2437	1.3003
2	Rainfall intensity map (mm/year)	1,549–1,781	0.1999	0.1009	0.5046
		1,781–1,874	0.1945	0.1635	0.8404
		1,874–1,908	0.1977	0.2241	1.1338
		1,908–1,975	0.2090	0.3069	1.4680
		1,975–2,124	0.1989	0.2047	1.0290
3	Distance to road map (m)	0–126	0.2114	0.2201	1.0412
		126–328	0.1978	0.1964	0.9931
		328–632	0.1972	0.1943	0.9853
		632–1,163	0.1968	0.1964	0.9979
		1,163–6,451	0.1968	0.1928	0.9795
4	Distance to river map (m)	0–340	0.2009	0.2262	1.1262
		340–1,020	0.1999	0.2369	1.1848
		1,020–2,254	0.1998	0.2030	1.0159
		2,254–4,465	0.1997	0.1422	0.7120
		4,465–10,845	0.1997	0.1917	0.9601
5	Distance to fault map (m)	0–15,944	0.2000	0.0386	0.1929
		15,944–29,115	0.2000	0.2029	1.0145
		29,155–53,378	0.2000	0.3375	1.6874
		53,378–68,975	0.2000	0.2845	1.4223
		68,975–88,386	0.2000	0.1366	0.6829
6	Drainage density (km/km2)	0	0.2331	0.2408	1.0331
		0–6	0.1917	0.1655	0.8632
		6–16	0.1917	0.2071	1.0801
		16–48	0.1917	0.2101	1.0960
		48–157	0.1917	0.1765	0.9206
7	Land-use map	Airport	0.0084	0.0022	0.2644
		Barren land	0.0009	0.0004	0.4666
		Dryland agriculture	0.0594	0.0076	0.1275
		Estate crop plantation	0.0035	0.0003	0.0739
		Fish pond	0.0475	0.0012	0.0250
		Rice field	0.3979	0.0650	0.1634
		Secondary mangrove forest	0.0003	0.0001	0.3295
		Settlement area	0.4572	0.9194	2.0110
		Shrub-mixed dryland farms	0.0198	0.0028	0.1440
		Swamp	0.0051	0.0009	0.1831
		Swamp shrub	0.0000	0.0000	0.1441
8	Lithology map	Alluvium	0.5020	0.5577	1.1111
		Alluvium fans	0.1966	0.3200	1.6274
		Beach ridge deposits	0.0153	0.0341	2.2385
		Bojongmanik form	0.0022	0.0000	0.0000
		Cihoe form	0.0161	0.0002	0.0149
		Coastal deposit	0.0013	0.0000	0.0000
		Lake	0.0001	0.0000	0.1178
		Marine deposits	0.0052	0.0000	0.0024
		Old alluvium	0.0074	0.0002	0.0239
		Parigi form	0.0005	0.0000	0.0000
		Sandstone tuff	0.0001	0.0000	0.0000
		Sandstone unit	0.0405	0.0440	1.0873
		Serpong form	0.0703	0.0014	0.0198
		Subang form	0.0064	0.0000	0.0030
		Swamp deposits	0.0175	0.0005	0.0300
		Tuff banten	0.0677	0.0391	0.5779
		Upper banten tuff	0.0488	0.0026	0.0543
		Young volcanic rocks	0.0022	0.0000	0.0000
9	Elevation map (m)	0–3	0.2227	0.1088	0.4883
		3–9	0.2029	0.2170	1.0699
		9–20	0.1976	0.3767	1.9061
		20–37	0.1968	0.2637	1.3397
		37–156	0.1800	0.0338	0.1878
10	Slope (degree)	0	0.2011	0.1727	0.8587
		0–1.36	0.1997	0.1657	0.8294
		1.36–3.19	0.1997	0.2399	1.2009
		3.19–5.47	0.1997	0.2470	1.2367
		>5.47	0.1997	0.1748	0.8753
	Aspect	Flat	0.1214	0.1339	1.1030
		North	0.1490	0.1778	1.1937
		Northeast	0.1091	0.1323	1.2128
		East	0.1177	0.1232	1.0468
11		Southeast	0.1006	0.0888	0.8829
		South	0.1212	0.0924	0.7626
		Southwest	0.0937	0.0711	0.7592
		West	0.0937	0.0880	0.9386
		Northwest	0.0937	0.0925	0.9869
12	Profile curvature	Concave	0.3332	0.3253	0.9764
		Flat	0.3336	0.3017	0.9045
		Convex	0.3332	0.3729	1.1192
13	Plan curvature	Concave	0.3332	0.3253	0.9764
		Flat	0.3336	0.3017	0.9045
		Convex	0.3332	0.3729	1.1192
	Topographic wetness index	2.52–6.81	0.1430	0.1533	1.0722
		6.81–8.00	0.1939	0.2252	1.1614
14		8.00–10.14	0.2140	0.2279	1.0647
		10.14–11.96	0.2258	0.2098	0.9293
		11.96–22.90	0.2233	0.1839	0.8232

Table 5

Calculated parameters for the algorithms used in this study.

Algorithm	Parameters
AdaBoost	The number of iterations: 10; seed: 1; weight threshold: 100.
LogitBoost	Number of iterations: 10; Seed: 1; weight threshold: 100; likelihood threshold: -1.7976E308; shrinkage: 1.0; max threshold: 3; thread pool: 1; thread to batch prediction: 1.
Logistic Regression	Ridge: 1.0E-8; max iterations: -1; number of decimal places: 4.
Multilayer Perceptron	Hidden layers: a; learning rate: 0.3; momentum: 0.2; number of decimal places: 2; seed: 0; training time: 500; validation set size: 0; validation threshold: 20.