A new GIS-based technique using an adaptive neuro-fuzzy inference system for land subsidence susceptibility mapping

ABSTRACT

In this study, we evaluated the predictive performance of an adaptive neuro-fuzzy inference system (ANFIS) with six different membership functions (MFs). Using a geographic information system (GIS), we applied ANFIS to land subsidence susceptibility mapping (LSSM) in the study area of Amol County, northern Iran. As a novelty, we derived a land subsidence inventory from the differential synthetic aperture radar interferometry (DInSAR) of two Sentinel-1 images. We used 70% of surface subsidence deformation areas for training, while 30% were reserved for testing and validation. We then investigated regions that are susceptible to subsidence via the ANFIS method and evaluated the resulting prediction maps using receiver operating characteristics (ROC) curves. Out of the six different versions, the most accurate map was generated with a Gaussian membership function, yielding an accuracy of 84%.

KEYWORDS:

• Sentinel-1
• land subsidence
• adaptive neuro-fuzzy inference system
• Amol County
• Iran

1. Introduction

Surface deformations like land subsidence are natural hazard phenomena, which can cause vast damage to the land (Feizizadeh et al. 2014a). Land subsidence moves in a vertical direction within a small angle (Xue et al. 2005, Shi et al. 2016). It can occur suddenly or gradually. The most prominent cause of this type of deformation is a structural weakness of the underground layers (Ashraf and Cawood 2015). This weakness is usually caused by human activities, such as mineral extraction from underground sources (Dehghani et al. 2014). In our study area, a decline in the aquifer level caused by several wells is one of the main reasons for land subsidence. This surface deformation is harmful to human lives, properties and the environment.

In this study, we used the DInSAR data for measuring surface subsidence in Amol County. This type of data has been used to measure/monitor the land subsidence in several previous studies; for instance, Ge et al. (2007) used DInSAR data to measure the land subsidence caused by mining activities. Amelung et al. (1999) used InSAR data to measure the subsidence rate before and after artificial recharge in the Las Vegas Valley. Bell et al. (2008) used the same data to detect ground surface variations caused by groundwater pumping and artificial recharge, both seasonally and in the long term. Kim et al. (2010) measured and predicted land subsidence areas using InSAR data in conjunction with hyperbolic models. Qu et al. (2015) used multi-temporal InSAR data to map the spatial and temporal differences of the land surface deformation in the Houston–Galveston region in Texas. In the case of Iran, Anderssohn et al. (2008) and Dehghan-Soraki et al. (2015) applied the same data to analyse and measure the surface deformation in the northeastern and northwestern parts of the country, respectively. Nowadays, DInSAR is considered as a sophisticated method of space geodesy to monitor different aspects of land deformations. DInSAR’s wide spatial coverage and its sensitivity to vertical surface deformations make it a popular method for measuring land subsidence (Karimzadeh 2016). One of the most recent works about land subsidence was carried out by Shrestha et al. (2017), who used the Shetran hydrological modelling system (http://research.ncl.ac.uk/shetran/), which is a system of integrated surface and subsurface modelling for land subsidence estimation.

The simultaneous use of remote sensing (RS), GIS and surface deformation observations, along with an appropriate model, has been shown to produce satisfactory results for monitoring natural hazards (Contreras et al. 2016). Generating models to identify areas of possible future land subsidence through susceptibility mapping is a widely accepted approach to hazard mitigation (Lee et al. 2012). Such models are useful for environmental planning managers and stakeholders to control the adverse impacts of this phenomenon (VaeziNejad et al. 2011). Land subsidence susceptibility mapping is a practical approach to identify high-risk regions. The spatial models with a variety of related factors and parameters have produced fairly good results for this aim (Feizizadeh and Blaschke 2013). GIS as a spatial model and well-organized datasets such as satellite images are necessary to monitor any natural hazard (Janssen 2007). Several recent studies have been conducted to map regions susceptible to natural hazards based on surface deformation. These studies used a range of methods, such as global positioning system survey methods (GPS, Abidin et al. 2007, Ge et al. 2014, 2017), artificial neural networks (ANN, Lee et al. 2012, Chen et al. 2017a), multi-criteria decision analysis (Feizizadeh and Blaschke 2014, Ashraf and Cawood 2015, Ghorbanzadeh et al. 2017), logistic regression analysis (Oh and Lee 2011, LR, Pham et al. 2016), decision trees (Lee and Park 2013), k-nearest-neighbour analysis (k-NN, Pradhan and Jebur 2017), weights of evidence and support vector machines (WOE and SVM, Goetz et al. 2015, Okwuashi and Ndehedehe 2017), neuro-fuzzy systems (Vahidnia et al. 2010) and adaptive neuro-fuzzy inference systems (ANFIS, Sezer et al. 2011, Bui et al. 2012, Basser et al. 2014, Polykretis et al. 2017, Chen et al. 2017b). A literature review indicates that a variety of methods have been used for natural hazard mapping and susceptibility prediction.

Combining new methods is a common approach to obtaining more accurate results of susceptibility mapping (Feizizadeh and Kienberger 2017). To this end, our study combines the ANFIS method and DInSAR data to achieve more accurate results. Once provided with input and training data, ANFIS can run a fuzzy inference structure (Cakıt and Karwowski 2017). In this structure, the membership functions(MFs) take advantage of the hybrid learning algorithm of the ANN to adjust themselves to output function parameters (Bui et al. 2017). The ANFIS method with different MFs has previously been applied to model natural hazards such as landslides. We used this approach along with the DInSAR data for land subsidence prediction. The proper land subsidence investigation technique depends on the subsidence’s spatial extent and its magnitude (Koros and Agustin 2017). Our study utilized GIS capabilities along with the ANFIS method for LSSM in a three-stage research design. First, we prepared input layers from different sources and generated a training data layer from DInSAR data. Second, we produced LSSMs using input and training data with different MFs. Third, we evaluated the resulting displacement map and the final results generated from each MF.

2. Study area

The study area was Amol County, which is situated in the Mazandaran province in northern Iran (see Figure 1). This county is located on the southern coastline of the Caspian Sea, on the Mazandaran plain of the Haraz river delta. This area is bound to the south by the Alborz mountain range, which is characterized by steep slopes, with angles reaching over 47%. However, the vast majority of slopes in the centre and the north of the study area have gradients of less than 2%. The vegetation cover is almost homogeneous in the whole study area, except for the mountainous border areas in the south, which are covered by forest.

Figure 1. Study area: (a) north of Iran; (b) Amol County; (c) surface deformation map.

In general, SAR images are sensitive to the change of vegetation cover. Using a pair of SAR images acquired at times of varying vegetation cover can lead to a low coherence, even within a period of one month (Ge et al. 2007). To avoid vegetation cover changes, we selected two SAR images that were captured on the same day of August in two separate years for the present study.

Our study area is suitable for a wide range of agricultural activities such as rice production and orchards, which are the main sources of income for the local population. Although agriculture in this area is primarily dependent on rainfall, there are several wells that have been used for groundwater extraction. The widespread exploitation of groundwater resources makes Amol County subject to land subsidence. According to the Iranian Bureau of Statistics 2017 census dataset, the population rate in Amol County has increased dramatically over the past three decades. This growth can also be seen by the footprint of Amol city. As a result of the population growth in our study area, agricultural activities have increased and, consequently, water demand and groundwater extraction have greatly increased.

3. Input data

3.1. Related input factors

There are nine factors related to land subsidence in our study area: groundwater level, slope, elevation, distance to streams, well density, rainfall, land use, aspect and lithology (see Figure 2). The lithological elements of the lithology layer were aggregated to five groups that are summarized in Table 1. The Landsat satellite data that served as the basis of the land cover map have a spatial resolution of 30 m. Similarly, the digital elevation model (DEM) of the SRTM, which was used for the generation of derivatives like the slope and aspect layers, also has a 30 m spatial resolution. However, in this study, we aggregated pixels to a pixel size of 90 m to accelerate computations and analysis. The aggregation is suitable for the subsidence study, considering the influence of each layer on this phenomenon. The layers were also prepared in raster format for further analysis.

Table 1. Lithology of Amol County.

Groups	Geo-units	Descriptions
Group 1	$Q_2$, $Q_{Haraz}$, $Q_1$	Old/young fans and terraces.
Group 2	${\mathrm{Q}}^{\mathrm{m}\mathrm{b}}$, ${\mathrm{Q}}^{\mathrm{c}\mathrm{p}}$	Coastal deposit with gravel and shell fragments.
Group 3	${\mathrm{Q}}^{\mathrm{a}\mathrm{l}}$	River beds and flood plain deposits.
Group 4	$Jd$, ${\mathrm{M}}^{\mathrm{m}\mathrm{s}}$, ${\mathrm{P}}^{\mathrm{l}\mathrm{q}}$	Marl, limestone and sandstone.
Group 5	$Tr{e}_2$, $K_2^{l,m}$, ${\mathrm{P}}^{\mathrm{e}\mathrm{m}}$	Thick bedded dolomitic, limestone/gypsum with marl and sandstone.

Figure 2. Nine input land subsidence factors used in the ANFIS model: (a) land use/cover; (b) ground water level; (c) aspect; (d) DEM; (e) distance to streams; (f) rainfall; (g) density of wells; (h) slope; (i) lithology.

3.2. Data generation for training and testing

In this subsection, we summarize the generation process of our training and test data. There are some common surface deformation measurement systems that are used to monitor and measure land subsidence, such as precise theodolite vertical angles and electronic levels (Savvaidis 2003). These common systems are really useful, and the approximated accuracies are close to 10 mm. They are straightforward and suitable for small regions (Koros and Agustin 2017). Furthermore, since these systems are ground point-based, they are usually affected by some problems in the process of field surveying and information gathering (Nagaomo et al. 2007). Nowadays, monitoring of ground deformation hazards such as land subsidence is possible using remote sensing satellites. The land subsidence inventory database of the study area was created using DInSAR data. In this study, the interferometric wide swath (IW) modes of Sentinel-1 images were used for further processing. This satellite (launched by the European environmental monitoring program Copernicus in 2014) presents synthetic aperture radar (SAR) data for ground deformation applications (Yu et al. 2017). Two sets of SAR images acquired at different times are the main input for the approach (Barra et al. 2017). To avoid the impacts of vegetation cover changes our first and second selected images were taken on 22 August 2015 and 22 August 2017, respectively. The Sentinel-1 IW interferometric pair was coregistered by the precise orbit state vectors and DEM (Jiang et al. 2017). The next step was to generate an interferogram. Because Sentinel-1 applies the TOPS (Terrain Observation by Progressive Scan) imaging mode (De Zan and Guarnieri 2006), it is able to cover a wide area (Yagüe-Martínez et al. 2016). This leads to images being captured within a series of overlapping regions (Czikhardt et al. 2017). Therefore, the small difference within overlap bursts is advantageous for retrieving the horizontal motion of the ground parallel to the satellite path (Grandin et al. 2016). The noise of these bursts was removed in a de-bursting process. De-bursted stacks were used to subtract the topographic phase. A Goldstein filter (Goldstein and Werner 1998) was applied to the results to improve the phase measurement. This filter has also been effective in reducing phase noise and before the implementation of unwrapping process (Notti et al. 2015). The resulting deformation map is presented in Figure 1, and all these steps are shown in Figure 4. The resulting deformation map provides detailed information on surface deformation for landslides and land subsidence. However, in this study, we focused on just land subsidence. For more certainty, the subsidences that occurred with a magnitude of over 0.02 metres were selected for both training and test data. Seventy percent (6654 pixels) of these areas were randomly selected as training data.

4. Workflow

4.1. Land subsidence susceptibility mapping using ANFIS

The ANFIS method contains a hybrid model with different nodes and layers in a neural network, which can estimate the fuzzy parameters of the model (Polykretis et al. 2017). This method is a combination of both fuzzy logic and artificial neural networks (Dehnavi et al. 2015). Figure 3 presents two input values of x and y, and f as an output (Jang 1993, Armaghani et al. 2015).

Figure 3. ANFIS architecture (Asadpour et al. 2013).

Rule 1: If x is A₁ and y is B₁ then

f ₁ = p₁x+ q₁y+ r₁

Rule 2: If x is A₂ and y is B₂ then

f ₂ = p₂x+ q₂y+ r₂

The membership functions of x and y are$A_1$, $A_2$, $B_1$, $B_2$, and output function parameters are $p_1$, $q_1$, $r_1$, $p_2$, $q_2$,$r_2$.

ANFIS consists of five different layers and two types of nodes, namely adaptive and fixed. The differences between the layers are illustrated as follows.

The input variables were fuzzified in the first layer

(1) $O_{1,i}=\mu A_{i\left(x\right)}$(1)

(2) $O_{1,i}=\mu B_{i\left(y\right)}$(2)

All nodes in the second layer are used to compute the strength of the rules. All nodes are fixed (shown by ∏ in Figure 3), and each of them determines the input and output signals (Zhang et al. 2017).

(3) $O_{2,i}={\omega }_i=\mu A_i\left(y\right)B_i\left(y\right)i=1,2$(3)

In the third layer, all nodes that are presented with an N are fixed (Asadpour et al. 2013). Further, the outputs are normalized for this layer.

(4) $O_{3,i}={\bar{\omega }}_i=\frac{w_i}{\sum _{j=1}^2{w}_j}i=1,2$(4)

All nodes in layer four are adaptive and described by the following equation:

(5) $O_{4,i}={\bar{\omega }}_i{f}_i={\bar{\omega }}_i(p_{i}x+q_{i}y+r_i)$(5)

This layer is the final layer. It is shown as Ʃ, which represents the summation. (Basser et al. 2014).

(6) $O_{5,i}=\sum _i{\bar{\omega }}_i{f}_i=\frac{\sum _i{\omega }_i{f}_i}{\sum _i{\omega }_i}i=1,2$(6)

4.2. The hybrid learning algorithm

The ANFIS method utilizes an ANN learning algorithm to set up the fuzzy inference system (Tahmasebi and Hezarkhani 2012). This algorithm was used for the training part of the process. The main aim of training is to discover the minimum value of the error for the fuzzy inference system. The error is the difference between the target and the results of the method. The ANN learning algorithm is fully described by Bui et al. (2012). The whole workflow is illustrated in Figure 4.

Figure 4. The flowchart of the integrated methodology for land subsidence susceptibility mapping and accuracy assessment.

5. Results

5.1. The interferometric performances

To validate the displacement map resulting from the interferometry process, we compared it against the ground truth data of the land subsidence inventory database, obtained from the Ministry of Water Resources of the Mazandaran Province (MWRMP, 2017) and the Ministry of Agricultural Resources for Mazandaran Province (MARMP, 2017). The land subsidence inventory database includes 19 land subsidence occurrences, or initial signs of land subsidence, within the study area; it consists of point data collected during field surveys using GPS locations (see Figure 5). The values of the resulting surface deformation map were extracted to the GPS points and are shown in Figure 6.

Figure 5. The extracted values of the surface deformation map to the GPS points. The first number represents the point number and the second number indicates the surface deformation map’s value at its corresponding point.

The validation metric of precision (Guirado et al. 2017) was used to assess the reliability of the resulting map. In this regard, the ground truth GPS points located in the subsidence areas depicted by the surface deformation map are considered as true positives (TP), and the rest are false positives (FP). The precision can be calculated from the Equation (7)(7) $\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}=\mathrm{T}\mathrm{P}\mathrm{s}/\left(\mathrm{T}\mathrm{P}\mathrm{s}+\mathrm{F}\mathrm{P}\mathrm{s}\right)$(7) :

(7) $\mathrm{P}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}=\mathrm{T}\mathrm{P}\mathrm{s}/\left(\mathrm{T}\mathrm{P}\mathrm{s}+\mathrm{F}\mathrm{P}\mathrm{s}\right)$(7)

It can be seen from Figure 5 that there is just one GPS point (point 9) from the ground truth inventory database that does not match our resulting deformation map. According to the analysis of these GPS points, the precision for the displacement map is over 94%.

5.2. The predictive performances

In this study, LSSM was carried out for the Amol County using DInSAR and ANFIS techniques. To identify a suitable mapping technique for land subsidence, we chose the nine most related factors as input data, and the surface displacement layer was defined as the training and test dataset. All input factors were included in the GIS environment as raster layers, and each layer is composed of 97,645 pixels in total. The surface displacement layer was prepared using two different SAR images of the study area from the years 2015 and 2017. The DInSAR approach was used to generate the surface displacement layer. The whole process of the model, including training and testing, was implemented in MATLAB. Seventy percent of land subsidence areas served as the training data, while the rest were divided randomly into two separate parts consisting of 2852 pixels for testing and the same number of pixels for evaluation (see Figure 6). Our ANFIS method was implemented with six different types of MFs. Fifty epochs were used for training.

5.3. Comparison of land subsidence susceptibility maps

In this subsection, we summarize the performance assessments of the method. The performance assessments were tested with two validating statistical parameters, namely, the root mean square error (RMSE; Equation (8)(8) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{N}{\sum }_{i=1}^n{\left(fou{t}_i-t_i\right)}^2}$(8) )(8) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{N}{\sum }_{i=1}^n{\left(fou{t}_i-t_i\right)}^2}$(8) and the coefficient of determination (R²; Equation (9)(9) ${\mathrm{R}}^2=1-\left[\frac{\sum _{i=1}^n{\left(fou{t}_i-t_i\right)}^2}{\sum _{i=1}^{n}fou{t_i}^2}\right]$(9) ), as recommended by other researchers, such as Bui et al. (2012), Folorunsho et al. (2012), Chai and Draxler (2014) and Singh et al. (2017).

(8) $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{N}{\sum }_{i=1}^n{\left(fou{t}_i-t_i\right)}^2}$(8)

(9) ${\mathrm{R}}^2=1-\left[\frac{\sum _{i=1}^n{\left(fou{t}_i-t_i\right)}^2}{\sum _{i=1}^{n}fou{t_i}^2}\right]$(9)

where $fou{t}_i$ and $t_i$ are the ith observed and result of the model respectively. Low RMSE values indicate better results (Popoola et al. 2017), while R² values closer to 1 show a better fit between the results and training data (Bui et al. 2012). According to Table 2, the approach revealed the best performance when using a Gaussian MF, which yielded a RMSE value of more than 0.15 and an R² value of 0.9253. The values of the resulting maps vary from 0 to 1, where a value of 0 indicates a very low susceptibility to land subsidence and 1 a very high susceptibility (see Figure 7)

Table 2. Description of MFs and statistical criteria for RMSE, R².

No	Type of MF	Descriptions	R^2	RMSE
1	Gauss MF	Gaussian curve membership function.	0.9253	0.1531
2	Pi MF	Pi-shaped membership function.	0.9037	0.1569
3	Tri MF	Triangular-shaped membership function.	0.9195	0.1673
4	Trap MF	Trapezoidal-shaped membership function.	0.9048	0.1507
5	Dsig MF	Difference of two sigmoid membership functions.	0.9032	0.1702
6	Gbell MF	Generalized bell curve membership function.	0.913	0.1701

6. Validation and discussion

Validation is considered to be an important phase in preparing susceptibility maps to predict possible future natural hazards (Pourghasemi et al. 2012, Feizizadeh et al. 2014b). In our study, the layer of surface deformation generated from Sentinel-1 data was used for this aim (see Figure 1). As previously mentioned, 2852 pixels of land subsidence areas were reserved for the validation section. The receiver operating characteristics (ROC) curves (Shadman Roodposhti et al. 2016) were used to validate the resulting maps within validation data. This validation approach shows the full scene of trade-off between the rate of TPs and that of FPs at all possible levels of positivity for the studied case (Liu et al. 2005, Baird et al. 2013). These curves were calculated for all six resulting LSSMs. The vertical axis presents the rate of true positives, while the horizontal axis shows the rate of false positives (see Figure 8) (Sezer et al. 2011, Feizizadeh and Ghorbanzadeh 2017). TPs are the pixels that correctly referred to land subsidence and FPs are the pixels incorrectly labelled as land subsidence. Generating ROC curves involves first obtaining the correctly and incorrectly labelled pixels in the validation data and then plotting TPs vs. FPs across the whole range of values (Linden 2006). Based on the theory behind the ROC, the area under the curve (AUC) is a single measure that represents the accuracy of the results of the prediction methods (Oh and Pradhan 2011). The AUC indicates the probability that the value of a correctly labelled pixel will exceed the value of an incorrectly labelled pixel. One of the benefits of this measure is that its output values are easy to interpret. A greater AUC value shows a higher accuracy of the model. A value close to 1 represents a perfect model with the highest accuracy; a value of 0.5 indicates a worthless model (the results of the model are the same as chance). Although there is no defined standard for comparing the results of the AUC, there are some studies that have employed a number of intervals with linguistic measurements or meaning (Baird et al. 2013). For instance:

• 0.90–1 = excellent

• 0.80–0.90 = good

• 0.70–0.80 = fair

• 0.60–0.70 = poor

• 0.50–0.60 = fail

According to these intervals, almost all of our results are categorized as good performance.

Although all the results of the ANFIS method with different membership functions are relatively close to each other and categorized on the same level, Gbell MF (0.882) yielded the best result (see Figure 8). This MF approach also achieved high R² values and low values of RMSE in the comparison of susceptibility map appearances. This matter indicates a meaningful relationship between the accuracy of the model within the training and testing process and also the overall accuracy of the model’s outputs. In this study, using the ANFIS method with the Gbell MF approach resulted in a fairly good model for land subsidence prediction. The profiles of cumulative vertical surface deformation predicted a probability represented in Figure 7, which also shows that the pixel values of LSSM derived from the Gbell MF vary at a moderate interval. Such moderate variation is more logical for our case study, since, based on the subsidence information from DInSAR (see Figure 1), there is not a severe deformation in this plain. The produced LSSMs show that the northeastern parts of the study area are more susceptible to land subsidence. This shows the impact of groundwater extraction on this phenomenon. According to the layer of well-density, there are more wells in the northeastern areas, which have led to a reduction in the groundwater level. Therefore, it is fair to say that land subsidence in Amol County is caused by groundwater extraction.

Figure 6. Training and validation data, with the location of GPS points.

Figure 7. Land subsidence susceptibility maps using the ANFIS model with six different MFs, namely: (a) Dsig MF; (b) Gauss MF; (c) Gbell MF; (d) Pi MF; (e) Trap MF; (f) Tri MF.

Figure 8. ROC curves for the resulting maps of the six MFs.

7. Conclusion

Land subsidence is one of the most severe natural hazards in Iran. Amol County is particularly susceptible to surface deformation with potentially severe consequences. GIS techniques along with ANFIS were used for LSSM in Amol County using six different MFs. The ANFIS method was structured for land subsidence susceptibility mapping. Sentinel-1 data were used for generating both training and test data, and a validation step was done for the final results. As the quality of the results of the ANFIS method highly depends on the quality of the training data, the surface deformation layer generated within the DInSAR approach provided fairly appropriate data for LSSM. Moreover, the accuracy of this method depends on the selection of the MFs. In summary, the ANFIS method is a technical method which is sensitive to the quality and quantity of training data and parameter selection during the whole process. However, the method is less dependent on user decisions and expert knowledge, which are usually prone to uncertainties. The drawback of the method is that the absolute importance of each factor is relatively hazy. The results of this work can be useful for prioritization of efforts to mitigate adverse consequences of surface deformations, specially land subsidence.

Acknowledgements

We would like to thank the Ministries of Water Resource and Agriculture of the Mazandaran Province, Iran, for their support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research is partly funded by the Austrian Science Fund (FWF) through the GIScience Doctoral College (DK W 1237-N23).