Land Subsidence Susceptibility Mapping Using Persistent Scatterer SAR Interferometry Technique and Optimized Hybrid Machine Learning Algorithms

: In this paper, land subsidence susceptibility was assessed for Shahryar County in Iran using the adaptive neuro-fuzzy inference system (ANFIS) machine learning algorithm. Another aim of the present paper was to assess if ensembles of ANFIS with two meta-heuristic algorithms (imperialist competitive algorithm (ICA) and gray wolf optimization (GWO)) would yield a better prediction performance. A remote sensing synthetic aperture radar (SAR) dataset from 2019 to 2020 and the persistent-scatterer SAR interferometry (PS-InSAR) technique were used to obtain a land subsidence inventory of the study area and use it for training and testing models. Resulting PS points were divided into two parts of 70% and 30% for training and testing the models, respectively. For susceptibility analysis, eleven conditioning factors were taken into account: the altitude, slope, aspect, plan curvature, proﬁle curvature, topographic wetness index (TWI), distance to stream, distance to road, stream density, groundwater drawdown, and land use/land cover (LULC). A frequency ratio (FR) was applied to assess the correlation of factors to subsidence occurrence. The prediction power of the models and their generated land subsidence susceptibility maps (LSSMs) were validated using the root mean square error (RMSE) value and area under curve of receiver operating characteristic (AUC-ROC) analysis. The ROC results showed that ANFIS-ICA had the best accuracy (0.932) among the models (ANFIS-GWO (0.926), ANFIS (0.908)). The results of this work showed that optimizing ANFIS with meta-heuristics considerably improves LSSM accuracy although ANFIS alone had an acceptable result.


Introduction
Land subsidence (LS) is one of the most challenging catastrophic geohazards due to its potential consequences, including damage to infrastructures, power lines, and buildings, causing sinkholes, floods in coastal areas, and soil degradation [1][2][3].Land subsidence is a gradual and slow deformation or sudden collapse of the Earth's surface, which is caused by numerous natural and human-induced factors [4][5][6][7].The ground subsiding movement can be the result of natural causes such as floods, ground lithology, dissolution of carbonated rocks (e.g., limestone), sediment compaction, and tectonic motions of faults [8][9][10][11].Further, anthropogenic activities that alleviate these geological factors, including underground excavations (e.g., mining and tunneling), underground resource withdrawal (gas or oil), overloading of the land surface through road construction and extending the built environment [12][13][14], and most importantly, over-exploitation of underground aquifers [15][16][17].
In the last few decades, the land subsidence phenomenon has widely increased in Iran [18][19][20], and therefore, growing research interest has focused on studying this geological problem [2,13,17,[21][22][23].One of the most important causes of the LS in Iran with an arid and semi-arid environment is the excessive groundwater extraction for agricultural usage [13,24].Therefore, modeling factors affecting land subsidence and land subsidence susceptibility mapping (LSSM) is vital for the environment, safety, economy, and human well-being.
However, these methods are mainly based on human assumptions and need expert knowledge.Recently, GIS-based machine learning algorithms (MLAs) have become a favorite in modeling and analyzing environmental hazards, especially LS.They can cope with data peculiarities, reveal complex relationships between data, and produce high accuracy and close-to-real world results [10,23].Lee and Park [12] conducted a comparative investigation between the decision tree (DT) algorithm and FR model in estimation of LS and its causing factors.Abdollahi et al. [34] applied a support vector machine (SVM) to predict LSS using water table drawdown and other influential factors.Taravatrooy et al. [35] used a hybrid clustering method based on k-means, genetic optimization, and several soft computing algorithms to examine subsided zones.Tien Bui et al. [10] compared four MLAs (Bayesian logistic regression (BLR), SVM, logistic model tree (LMT), and alternate decision tree (ADT)) in assessing LSS near abandoned mining areas in South Korea.In a study in Kerman, Iran, the random forest (RF) algorithm showed superior capability in LSS mapping [36].Ebrahimy et al. [23] performed a comparative study using three tree-based MLAs, a boosted regression tree (BRT), RF, and classification and regression tree (CART), for studying land susceptibility in Tasuj plane, Iran.Evaluation results revealed that BRT had the best performance.In another study by Rahmati et al. [37], four tree-based MLAs, a rule-based decision tree (RDT), RF, CART, and BRT, were compared for generating LS hazard maps in Hamedan Province, Iran.The results indicated that RF had the best accuracy amongst the employed methods.
Despite the better performance and accuracy of MLAs, all the above-mentioned approaches are dependent on the availability and precision of the subsidence inventory data, which is a serious challenge in developing countries [2,11,37].On the other hand, interferometric synthetic aperture radar (InSAR) has been utilized in land displacement measurement and demonstrated promising results with millimetric precision [38,39].SAR is satellite data so there is no need for time-consuming field survey data acquisition; therefore, it is superior to other approaches such as leveling data and is denser than ground positioning system (GPS) station data.Furthermore, radar data are functional in all-time all-weather conditions, making this a cost-efficient method to obtain land subsidence measurements.Recently, InSAR methods were employed in LSS studies as reliable input data along with other data to achieve finer accuracies [6,40].In this paper, we used the PS-InSAR method to obtain land subsidence inventory data and utilize them among other subsidence triggering factors for land subsidence susceptibility mapping in the study area.As a novel methodology in LSSM, we used ANFIS optimized with two meta-heuristic algorithms: (1) imperialist competitive algorithm (ICA) and (2) gray wolf optimization (GWO).ANFIS uses hybrid learning of ANN in adjusting its membership functions (MF) with output data [41,42].Further, by taking advantage of meta-heuristic algorithms, the weight parameters of MFs were optimized.The results of the method were compared using the statistical approach of the root mean square error (RMSE) value.Furthermore, the accuracy of LSSMs was evaluated by the area under the receiver operating characteristic (ROC) curve.

Materials and Methods
The methodology applied in this research (summarized in Figure 1) is to generate an updated land subsidence inventory through the PSInSAR technique using a Sentinel-1 SAR dataset spanning the period from 2019 to 2020.Moreover, the aim of this paper is to develop and use an ensemble of ANFIS and meta-heuristic algorithms in modeling land subsidence susceptibility.The main steps of the study are as follows.First, a spatial database was created using generated land subsidence inventory and the layers of the conditioning factors.In the second step, PS-InSAR-derived subsidence inventory data were divided into training (70%) and testing (30%) data.Next, MF parameters were optimized in ANFIS using ICA and GWO meta-heuristic algorithms, and then LSSMs were produced using ANFIS, ANFIS-ICA, and ANFIS-GWO individually.Finally, the produced land susceptibility maps were compared and evaluated using the area under the ROC curves.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 25 area.As a novel methodology in LSSM, we used ANFIS optimized with two meta-heuristic algorithms: (1) imperialist competitive algorithm (ICA) and (2) gray wolf optimization (GWO).ANFIS uses hybrid learning of ANN in adjusting its membership functions (MF) with output data [41,42].Further, by taking advantage of meta-heuristic algorithms, the weight parameters of MFs were optimized.The results of the method were compared using the statistical approach of the root mean square error (RMSE) value.Furthermore, the accuracy of LSSMs was evaluated by the area under the receiver operating characteristic (ROC) curve.

Materials and Methods
The methodology applied in this research (summarized in Figure 1) is to generate an updated land subsidence inventory through the PSInSAR technique using a Sentinel-1 SAR dataset spanning the period from 2019 to 2020.Moreover, the aim of this paper is to develop and use an ensemble of ANFIS and meta-heuristic algorithms in modeling land subsidence susceptibility.The main steps of the study are as follows.First, a spatial database was created using generated land subsidence inventory and the layers of the conditioning factors.In the second step, PS-InSAR-derived subsidence inventory data were divided into training (70%) and testing (30%) data.Next, MF parameters were optimized in ANFIS using ICA and GWO meta-heuristic algorithms, and then LSSMs were produced using ANFIS, ANFIS-ICA, and ANFIS-GWO individually.Finally, the produced land susceptibility maps were compared and evaluated using the area under the ROC curves.

Study Area
The region of interest in this paper is Shahryar, the central city of Shahryar County within 35  6 longitudes (Figure 2), with the elevation ranging from 1081 to 1222 m.This county is located west of Tehran, the capital of Iran.In recent years, there has been an increase in population migration to cities near Tehran, including Shahryar, for better jobs and income.According to Iran census data in 2016, the county is the 12th largest in the country with a population of more than 700,000 people.This has become a serious problem in urban environment management and food production.Shahryar County is known for its green and beautiful landscape and the major income of the people in the area originates from gardening and agriculture.Owing to population increase, the demand for food has grown dramatically.Therefore, more illegal wells are dug.As a result, the county has suffered from severe land subsidence (250 to 310 mm/year) due to exhaustion of underground water aquifers and water table dropdown.

Study Area
The region of interest in this paper is Shahryar, the central city of Shahryar County within 35° 35′ to 35° 42′ latitudes and 50° 59′ to 51° 6′ longitudes (Figure 2), with the elevation ranging from 1081 to 1222 m.This county is located west of Tehran, the capital of Iran.In recent years, there has been an increase in population migration to cities near Tehran, including Shahryar, for better jobs and income.According to Iran census data in 2016, the county is the 12th largest in the country with a population of more than 700,000 people.This has become a serious problem in urban environment management and food production.Shahryar County is known for its green and beautiful landscape and the major income of the people in the area originates from gardening and agriculture.Owing to population increase, the demand for food has grown dramatically.Therefore, more illegal wells are dug.As a result, the county has suffered from severe land subsidence (250 to 310 mm/year) due to exhaustion of underground water aquifers and water table dropdown.

Date Used
As outlined earlier, the methodology has two main parts.One is the procedure of generating a land subsidence inventory using the PS-InSAR technique, which needs a satellite SAR dataset.The other is the conditioning factors that are taken into consideration for hazard modeling.In the following section, the process of acquiring and preparing these datasets is discussed.

SAR Data
To generate a land subsidence inventory in order to acquire the training and testing data necessary for LSSM, sentinel-1A single look complex (SLC) SAR data provided by European Space Agency (ESA) were used.In total, 31 SAR images were obtained from January 2019 to January 2020 in ascending track with dual polarization (vertical-vertical and vertical-horizontal) (Table 1).

Date Used
As outlined earlier, the methodology has two main parts.One is the procedure of generating a land subsidence inventory using the PS-InSAR technique, which needs a satellite SAR dataset.The other is the conditioning factors that are taken into consideration for hazard modeling.In the following section, the process of acquiring and preparing these datasets is discussed.

SAR Data
To generate a land subsidence inventory in order to acquire the training and testing data necessary for LSSM, sentinel-1A single look complex (SLC) SAR data provided by European Space Agency (ESA) were used.In total, 31 SAR images were obtained from January 2019 to January 2020 in ascending track with dual polarization (vertical-vertical and vertical-horizontal) (Table 1).To date, there is no single guideline unanimously applied for selecting subsidenceaffecting factors.However, based on the literature of LSSM studies carried out in Iran [13,23,33,37], data that were directly and indirectly related to the subsidence-inducing factors are gathered and used for mapping LSS in the study area.These input data include the altitude, slope, aspect, plan curvature, profile curvature, topographic wetness index (TWI), distance to stream, distance to road, stream density, groundwater drawdown, and land use/land cover (LULC) (Table 2).Altitude (Figure 3a) and its derivative factors such as slope (Figure 3b), aspect (Figure 3c), TWI (Figure 3d), plan (Figure 3e), and profile (Figure 3f) are among crucial topo-hydrological criteria of ground subsidence [2].An advanced space borne thermal emission and reflection radiometer (ASTER) digital elevation model (DEM) was obtained (1 arc second or approximately 30 m resolution) and processed in the GIS environment to produce topography-related factor layers.Land subsidence can cause deformations in the Earth's surface slope and topography [43].Therefore, altitude and its derivatives were considered because they can directly affect LS.The slope can have a potential impact on runoff infiltration since steep slopes bring about less recharge due to limited infiltration of rainfall [2,37].As mentioned above, the main cause of land subsidence in Iran is excessive underground water extraction.Undue groundwater extraction results in pore water pressure (PWP) decrease and aquifer compaction increase [44].Therefore, the well inventory of the study area was acquired from the Iranian department of water resource management (IDWRM) to generate groundwater drawdown (Figure 3k).Further, distance to stream (Figure 3i) and stream density (Figure 3j) were used, as these factors can impact the groundwater level by recharging the groundwater tables [45].The network of streams was extracted from the ASTER DEM.Another geomorphological influential factor is the ground lithology of the area.However, we could not take the lithological layer into account since this factor did not have substantial variation in the region of interest.Road data were obtained from the open street map (OSM) at the scale of 1:100,000, and distance to a road was calculated in the GIS environment (Figure 3h).Google Earth Engine (GEE) cloud computing, gathering massive volumes of various satellite imagery alongside popular machine learning algorithms, is a suitable platform for analyzing geo big data and monitoring the environment [46,47].The available 10-m Iran-wide LULC map was used.The map was generated in the GEE platform using Sentinel-1 and Sentinel-2 images and object-based random forest classifier with 95% overall accuracy [48].All the data were generated or resampled to a 30 m pixel size.Groundwater drawdown Well inventory of the study area

Methodology
In this section, the methods and models used in different parts of the approach are presented in detail.First, the PS-InSAR technique used for generating land subsidence map of the study area and the LS inventory are discussed.Then, the ANFIS model and evolutionary algorithms, GWO and ICA, are stated.Moreover, FR analysis and the approach to optimize ANFIS using meta-heuristics are presented.Finally, the validation methodology of this paper is outlined.

PS-InSAR technique
SAR Interferometry (InSAR) is a well-established technique for measuring and monitoring ground deformation with millimetric accuracy.This is mainly based on the phase difference between SAR images acquired at different times and slightly different sensor

Methodology
In this section, the methods and models used in different parts of the approach are presented in detail.First, the PS-InSAR technique used for generating land subsidence map of the study area and the LS inventory are discussed.Then, the ANFIS model and evolutionary algorithms, GWO and ICA, are stated.Moreover, FR analysis and the approach to optimize ANFIS using meta-heuristics are presented.Finally, the validation methodology of this paper is outlined.

PS-InSAR Technique
SAR Interferometry (InSAR) is a well-established technique for measuring and monitoring ground deformation with millimetric accuracy.This is mainly based on the phase difference between SAR images acquired at different times and slightly different sensor positions.Time-series InSAR analysis aims to identify coherent image pixels (persistent scatterers (PSs)), which have high phase stability and reflect strong backscatter to the satellite over a long time period.A baseline configuration can determine the set of interferometric image pairs, which is used in the time series analysis.The baseline shows the distance between two images, in terms of antenna position (spatial baseline), acquisition time (temporal baseline), or Doppler centroid (Doppler baseline).The single-master configuration, where each image is co-registered to a unique master image, is the most common one for PSI analysis [49,50].The master image is chosen in the middle of the 2D spatio-temporal space so that the high coherence of all formed interferometric pairs (interferograms) is guaranteed.The interferogram contains the ground deformation phase component as well as some other distinct contributions, such as atmospheric disturbance, topographic, and flat-Earth terms.These components are removed in the next step from the interferometric phase using an external DEM.
A spatial network is formed using a primary set of points as PS candidates (PSCs) to estimate the Atmospheric Phase Screen (APS) and further densification of PS points.Since the original interferometric phase is wrapped (i.e., phase observations in the [−π, +π)), and it is composed of a large number of phase contributions, the PSCs cannot be selected based on the phase.Thus, the amplitude dispersion index (ADI) based on Equation ( 1) is used as an approximation of phase stability.A point can be selected as a PSC if it always has a higher amplitude value than a suitable threshold.It was proved that assuming sufficient data images, the phase behavior with the standard deviation (σ v ) lower than a threshold of 0.25 is similar to the trend of ADI.Hence, this index, which represents the phase stability of points, can be used for PSCs selection.
where σ A is the standard deviation and m A is the average amplitude value of each pixel over time.Next, the spatial network is used to estimate the unknown parameters, DEM error (residual topographic phase component), and the deformation rate, along with each connection between two adjacent PSCs through the maximization of a periodogram (Equation ( 2)) [51].All PSI methods are based on assumptions regarding the spatial and temporal smoothness of the deformation signal, expressed by a model.Here, the model is considered as a linear deformation trend in time.
where p ij demonstrates the connection between adjacent PSCs p i and p j .N is the number of interferograms.The term ∆ϕ s.k is the double difference interferometric phase in image pairs s and k, while B t.s and B n.s are the temporal and interferometric normal baselines, respectively.θ refers to the incidence angle of the SAR signal.For each PSC, the average residual phase after correction for the modeled parameters is taken to obtain an estimate for the atmospheric signal in the master image.Then, the atmospheric signal of the slave images as well as phase noise is separated from un-modeled deformation based on high-pass filtering.After APS removal, to increase point density, the second set of PSs is selected using a higher threshold for the ADI criterion.The unknown parameters are re-estimated for all pixels based on another maximization of the periodogram [49].Eventually, temporal coherence, which is a function of residual phase noise, is used to determine the final PSs which build the land deformation map.

Adaptive Neuro-Fuzzy Inference System (ANFIS)
Fuzzy inference systems (FIS) are capable of depicting multifaceted processes using the concepts and if-then rules.However, they are incapable of learning [52].Furthermore, if the number of input variables is too large, then selecting the membership functions and setting the fuzzy rules will become challenging [53].On the other hand, learning algorithms automatically choose the suitable set of parameters for fuzzy membership functions despite their inability to explain the system under study.Thus, the ANFIS model [54] is a mixture of ANN and FL benefiting from both ANN's computation capability and FL's decision making.The structure of the ANFIS model contains five layers, called adaptive and fixed [55].The ANFIS model employs the Takagi-Sugeno-Kang fuzzy algorithm in two rules of 'if-then' with two inputs, x and y, and one output f for both as follows [56]: Each node contains adaptive nodes, and input variables are fuzzified in first layer (Equations ( 5) and ( 6)): where, x and y are the input nodes, A and B are the linguistic variables, and µA i(x) and µB i(y) are membership functions for that node.The second layer contains fixed nodes denoted as Fuzzy inference systems (FIS) are capable of depicting multifaceted processes using the concepts and if-then rules.However, they are incapable of learning [52].Furthermore, if the number of input variables is too large, then selecting the membership functions and setting the fuzzy rules will become challenging [53].On the other hand, learning algorithms automatically choose the suitable set of parameters for fuzzy membership functions despite their inability to explain the system under study.Thus, the ANFIS model [54] is a mixture of ANN and FL benefiting from both ANN's computation capability and FL's decision making.The structure of the ANFIS model contains five layers, called adaptive and fixed [55].The ANFIS model employs the Takagi-Sugeno-Kang fuzzy algorithm in two rules of 'if-then' with two inputs, and , and one output for both as follows [56]: 1: , , Each node contains adaptive nodes, and input variables are fuzzified in first layer (Equations ( 5) and ( 6)): where, x and y are the input nodes, A and B are the linguistic variables, and μ ( ) and μ ( ) are membership functions for that node.The second layer contains fixed nodes denoted as ᴫ to compute the strength of the rules.The output of each node is the product of all input signals to that node (Equation ( 7)): where is the output for each node.The third layer encompasses fixed nodes denoted as N.The nodes in this layer are the normalized outputs of the second layer, which are referred to as the normal firepower (Equation ( 8)): All nodes in the fourth layer are adaptive and associated with a node function described by the following equation: where is the normalized firepower of third layer and , , and are node parameters.The parameters of this layer are can be interpreted as the result parameters.
The final layer has only one node denoted as Ʃ, which represents the summation of all the input signals:

Imperialist Competitive Algorithm
The imperialist competitive algorithm (ICA) is a novel evolutionary algorithm based on human social evolution, developed by Atashpaz-Gargari and Lucas [57].The ICA belongs to the group of swarm intelligence, which provides a powerful algorithm for solving to compute the strength of the rules.The output of each node is the product of all input signals to that node (Equation ( 7)): where W i is the output for each node.
The third layer encompasses fixed nodes denoted as N.The nodes in this layer are the normalized outputs of the second layer, which are referred to as the normal firepower (Equation ( 8)): All nodes in the fourth layer are adaptive and associated with a node function described by the following equation: where w i is the normalized firepower of third layer and p i , q i , and r i are node parameters.The parameters of this layer are can be interpreted as the result parameters.
The final layer has only one node denoted as Σ, which represents the summation of all the input signals:

Imperialist Competitive Algorithm
The imperialist competitive algorithm (ICA) is a novel evolutionary algorithm based on human social evolution, developed by Atashpaz-Gargari and Lucas [57].The ICA belongs to the group of swarm intelligence, which provides a powerful algorithm for solving NP-hard problems through its capability of dealing with continuous optimization.In this algorithm, the primary population is composed of several countries, and they interact with each other to form empires. Assuming the value of the objective function, colonialist and colonial groups are formed based on the existing countries.After ascertaining a colonialist, other countries are randomly allocated to one of the colonizers [57,58].Every colonialist and its associated colony is called an empire.The algorithm then simulates the competition among imperialists in order to acquire more colonies.The best colonialist typically has more chance to occupy more colonies.
Another way of allocating the colonies to each colonialist is based on their normalized cost, which is calculated via Equation (11): where T.c n is the empire's total cost n, and N.T.C n is the total normalized cost value of that empire.Possession eventuality of the colonization competition by each empire is calculated by Equation ( 12): The next phase is to attempt to approach a colonial country to analyze the colonies' cultural and social structures in different political and social layers.The colonies then move to the colonialist country.The colonialist and colony will change their positions; the new colonialist position will continue with the algorithm.The new colonialist country will start applying adjustment to its colonies this time.To calculate the cost function, the total empire cost is given by Equation ( 13) as follows: T.c n = f (imp) + algorithm, the primary population is composed of several countries, and they interact with each other to form empires. Assuming the value of the objective function, colonialist and colonial groups are formed based on the existing countries.After ascertaining a colonialist, other countries are randomly allocated to one of the colonizers [57,58].Every colonialist and its associated colony is called an empire.The algorithm then simulates the competition among imperialists in order to acquire more colonies.The best colonialist typically has more chance to occupy more colonies.
Another way of allocating the colonies to each colonialist is based on their normalized cost, which is calculated via Equation (11): . .= .− ., (11) where .
is the empire's total cost n, and . . is the total normalized cost value of that empire.Possession eventuality of the colonization competition by each empire is calculated by Equation ( 12): The next phase is to attempt to approach a colonial country to analyze the colonies' cultural and social structures in different political and social layers.The colonies then move to the colonialist country.The colonialist and colony will change their positions; the new colonialist position will continue with the algorithm.The new colonialist country will start applying adjustment to its colonies this time.To calculate the cost function, the total empire cost is given by Equation ( 13) as follows: . = ( where f(imp) is the value of the cost function for the colonialist, f(colony) represents the mean values of the cost function for the colonies, and the constant ᵹ is considered a value between 0 and 1.
Finally, the cost of each empire is calculated, and the colonies of weak empires are abolished and join to stronger empires.This process of recruitment or competition between colonialists is continued.In the next stage, empires that have lost all their colonies will be eliminated and will join other colonies.The process is repeated until a single universal empire in the globe is built that is very close to the empire with colonial nations [59].

Grey Wolf Optimization
The novel grey wolf optimization (GWO) algorithm, presented by Mirjalili et al. [60], is an inspiration from the hunting behavior of grey wolves and the social hierarchy in nature.Wolves are social animals that live in packs, and they have a hierarchy consisting of four groups.The leader of each group, the alpha wolf (α), makes decisions about hunting, sleeping, and walking time, and all the other group members must follow its directives.In terms of hierarchy, the other wolves fall into three levels, called beta (β), delta (δ), and omega (ω).The beta wolves at the second level assist the alpha in making decisions and devise them.They are the best candidates for alpha replacement.Another notable characteristic is their group hunting, which can be summarized in four stages: (1) encircling prey, (2) hunting, (3) attacking prey (exploitation), and (4) searching for prey (exploring) [61].The hunting process (optimization) is led by α, β, and δ wolves, and ω wolves have to abide by these three groups.
1. Encircling prey: where f (imp) is the value of the cost function for the colonialist, f (colony) represents the mean values of the cost function for the colonies, and the constant NP-hard problems through its capability of dealing with continuous op algorithm, the primary population is composed of several countries, with each other to form empires. Assuming the value of the objective fu and colonial groups are formed based on the existing countries.After as nialist, other countries are randomly allocated to one of the colonizers [5 nialist and its associated colony is called an empire.The algorithm th competition among imperialists in order to acquire more colonies.Th typically has more chance to occupy more colonies.Another way of allocating the colonies to each colonialist is based ized cost, which is calculated via Equation (11): . .= .− ., where .
is the empire's total cost n, and . . is the total normal that empire.Possession eventuality of the colonization competition by e culated by Equation ( 12): The next phase is to attempt to approach a colonial country to ana cultural and social structures in different political and social layers.T move to the colonialist country.The colonialist and colony will change th new colonialist position will continue with the algorithm.The new colon start applying adjustment to its colonies this time.To calculate the cost f empire cost is given by Equation ( 13) as follows: . = ( ) + ᵹ ( ) where f(imp) is the value of the cost function for the colonialist, f(colon mean values of the cost function for the colonies, and the constant ᵹ is co between 0 and 1.
Finally, the cost of each empire is calculated, and the colonies of w abolished and join to stronger empires.This process of recruitment o tween colonialists is continued.In the next stage, empires that have lost will be eliminated and will join other colonies.The process is repeated versal empire in the globe is built that is very close to the empire with [59].

Grey Wolf Optimization
The novel grey wolf optimization (GWO) algorithm, presented by M is an inspiration from the hunting behavior of grey wolves and the so nature.Wolves are social animals that live in packs, and they have a hie of four groups.The leader of each group, the alpha wolf (α), makes deci ing, sleeping, and walking time, and all the other group members mus tives.In terms of hierarchy, the other wolves fall into three levels, called and omega (ω).The beta wolves at the second level assist the alpha in and devise them.They are the best candidates for alpha replacement.characteristic is their group hunting, which can be summarized in four cling prey, (2) hunting, (3) attacking prey (exploitation), and (4) search ploring) [61].The hunting process (optimization) is led by α, β, and wolves have to abide by these three groups.

Encircling prey:
is considered a value between 0 and 1.
Finally, the cost of each empire is calculated, and the colonies of weak empires are abolished and join to stronger empires.This process of recruitment or competition between colonialists is continued.In the next stage, empires that have lost all their colonies will be eliminated and will join other colonies.The process is repeated until a single universal empire in the globe is built that is very close to the empire with colonial nations [59].

Grey Wolf Optimization
The novel grey wolf optimization (GWO) algorithm, presented by Mirjalili et al. [60], is an inspiration from the hunting behavior of grey wolves and the social hierarchy in nature.Wolves are social animals that live in packs, and they have a hierarchy consisting of four groups.The leader of each group, the alpha wolf (α), makes decisions about hunting, sleeping, and walking time, and all the other group members must follow its directives.In terms of hierarchy, the other wolves fall into three levels, called beta (β), delta (δ), and omega (ω).The beta wolves at the second level assist the alpha in making decisions and devise them.They are the best candidates for alpha replacement.Another notable characteristic is their group hunting, which can be summarized in four stages: (1) encircling prey, (2) hunting, (3) attacking prey (exploitation), and (4) searching for prey (exploring) [61].The hunting process (optimization) is led by α, β, and δ wolves, and ω wolves have to abide by these three groups.

1.
Encircling prey In the first stage, the grey wolves encircle and surround the prey during hunting.To define this phase mathematically, the following Equations ( 14) and ( 15) are proposed.The parameter D measures the distance between the grey wolf and the prey, and → X represents the location of the prey: where t denotes the current iteration, and → X p and → X denote the position vectors of the prey and the grey wolves, respectively.The → A and → C coefficient vectors are defined as follows: where components of a are linearly decreased from 2 to 0 over the course of iterations, and r 1 and r 2 are random vectors between [0,1].

Hunting
After the encircling of the prey, the hunting phase is guided by α, β, and δ since they are supposed to have compressive knowledge about the prey's position.This can be computed using following formulas: where → X 3 denote the position of α, β and δ, wolves respectively.
→ C 3 are the respective coefficient vectors.The position of a grey wolf in the search space can be updated as follows: The other wolves update their positions randomly according to the position of the prey.

Attacking prey
The process of hunting ends when the prey stops moving and gray wolves attack the prey.It is important to note that the fluctuation range of where a is linearly decreased from 2 to 0. The exploration trend happens when → A < 1 and → C < 1.At this moment, the wolves attack the prey.

Search for prey
The grey wolves track and chase the prey.The pursuing of the prey is known as the exploration phase in the GWO algorithm [62].The parameters α, β, and δ have guidance responsibility in this process.If

Frequency Ratio
One of the statistical bivariate models is FR, which is widely used in modeling environmental hazards as a geospatial assessment tool for quantifying the potential relationship between dependent and independent variables [63].The FR value for a certain class from a given factor can be calculated using: where N pix (X i ) is the number of pixels in each class of each factor with land subsidence locations.X.N pix X j is the number of pixels of X j factor, m is the number of classes in the X i factor, and n is the number of factors in the study area [64].

ANFIS with Meta-Heuristic Algorithms
In ANFIS, parameter adjustment and the creation of a basic fuzzy system are done by combining traditional methods and then back error propagation.In this research, ICA and GWO were used as meta-heuristic algorithms to enhance the results of the ANFIS system and also to tweak the parameters of membership functions [52,65].First, using input and target data, the FIS is created by the ANFIS model.Next, the membership functions are optimized and adjusted by meta-heuristic algorithms, and the output for the ANFIS (y) model is computed by [66]: where t is the target data, y is a function of input data and optimized FIS, and e is the error function that should be minimized.When the final conditions are met with the best output, the optimization process stops; otherwise, the membership function optimization is repeated.

Validation
In this research, the ROC curves were used for the accuracy assessment of the LSS models employed [6,10,39].The ROC curve analysis is a common method to evaluate the goodness-of-fit and prediction power of models regarding the area under the curve (AUC) [2,67].Ranging from 0 to 1, higher AUC values represent more reliable and accurate model performance.According to Yesilnacar [53], the qualitative relationship between AUC and the prediction accuracy of a model can be classified into the following categories: 0.5-0.6 (poor), 0.6-0.7 (average), 0.7-0.8(good), 0.8-0.9(very good), and 0.9-1 (excellent).

Results
By utilizing the spatial data and subsidence inventory generated and the methods discussed above, the mapping and assessment of land subsidence susceptibility for Shahryar County were conducted.In the following sub-sections, the results of the various parts of the methodology are thoroughly discussed.

Result of PS-InSAR
The ground deformation rate along the line of sight (LOS) direction during the acquisition time interval was obtained based on the PS-InSAR technique (see Section 3.1).SARPROZ (SAR PROcessing tool by periZ) [51] software was used to implement the PS-InSAR technique in the current research paper.The star configuration of the S-1A SLC time-series data stack is shown in Figure 4.The master image for the dataset was selected based on maximizing the stack coherence [68].All slave images were co-registered with respect to the master image (2019/09/11).The shuttle radar topography mission (SRTM) DEM was used to remove topographic-related phase components from the interferometric phases.After the selection of 7881 PSCs, a spatial network was created by Delaunay Triangulation, connecting each point to the other.The unknown LOS velocity and DEM error were calculated along with the connections by maximizing a periodogram.All the obtained parameters were then integrated into the absolute values with respect to a reference point, which had no subsidence rate.The atmospheric phase for PSCs can be resampled on the uniform image grid as the APS.With having APS compensated for all slave images, the unknowns were estimated again for more PS points, selected by a lower threshold on ADI to obtain a more dense subsidence map. Figure 5a shows the LOS deformation map.According to Figure 5a, the maximum velocity was about −175 mm/year, which occurred in the southern agricultural part of the region of interest, where more ground-water extraction was observed.
Further, the LOS deformation rates should be decomposed into horizontal and vertical deformation components as it is the inherently vertical movement of the Earth's surface with a slight horizontal displacement.It has been proved that the horizontal deformation is a very small portion of motion compared to the vertical deformation [67,69].Hence, the LOS deformation could be assumed as negligible and converted simply into the vertical deformation rates using the cosine of the incidence angle of the radar signal.The interpolated land subsidence inventory map was designed based on the vertical velocity deformation map for Shahryar County, depicted in Figure 5b.ADI was more successful in identifying PS points in man-made areas with stable targets than agricultural areas [70].Since vegetated regions are the main land cover in our case study, an interpolation was applied to the obtained vertical map to extend the deformation information for the whole study area.The inverse distance weighted (IDW) interpolation was used with a weighting power of 2 and neighboring radius of 12 for calculating the vertical velocity interpolation.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 24 obtained parameters were then integrated into the absolute values with respect to a reference point, which had no subsidence rate.The atmospheric phase for PSCs can be resampled on the uniform image grid as the APS.With having APS compensated for all slave images, the unknowns were estimated again for more PS points, selected by a lower threshold on ADI to obtain a more dense subsidence map.Further, the LOS deformation rates should be decomposed into horizontal and vertical deformation components as it is the inherently vertical movement of the Earth's surface with a slight horizontal displacement.It has been proved that the horizontal deformation is a very small portion of motion compared to the vertical deformation [67,69].Hence, the LOS deformation could be assumed as negligible and converted simply into the vertical deformation rates using the cosine of the incidence angle of the radar signal.The interpolated land subsidence inventory map was designed based on the vertical velocity deformation map for Shahryar County, depicted in Figure 5(b).ADI was more suc-

Result of FR
The results of the FR analysis in identifying the relationship of land subsidence occurrence with the conditioning factors are summarized in Table 3.Two out of five altitude classes had the highest probability (FR > 1.0), with 1119 to 1137 meters being the most correlated class with land subsidence, followed by the altitude class lower than 1119 m.The results of the slope angle analysis showed that slopes ranging between 4.5 and 6.8 degrees had the highest FR (1.11).Further, a TWI class lower than 4.84, profile curvature higher than 0.0029, convex plan curvature class, and flat (F) slope aspect had the most influence on LSS for each corresponding factor.Land cover analysis results indicated that forest and urban classes had the highest probability of land subsidence occurrence, with FR values of 1.17 and 1.04, respectively.Distance to a stream of between 50 to 100 meters had the highest FR, and the class of lower than 50 meters had a considerable correlation; a stream density higher than 2.68 had the highest correlation with land subsidence and the 1.23 to 1.92 class, which also had a considerable FR value.For distance to road, the 0 to 100 meter class had the highest FR followed by the 100 to 200 m classes.Finally, groundwater drawdown ranging from 55 to 83 m and from 28 to 55 m had higher impacts on land subsidence occurrence.Table 3. Relationship between land subsidence occurrence and conditioning factors using the FR model.

Result of FR
The results of the FR analysis in identifying the relationship of land subsidence occurrence with the conditioning factors are summarized in Table 3.Two out of five altitude classes had the highest probability (FR > 1.0), with 1119 to 1137 m being the most correlated class with land subsidence, followed by the altitude class lower than 1119 m.The results of the slope angle analysis showed that slopes ranging between 4.5 and 6.8 degrees had the highest FR (1.11).Further, a TWI class lower than 4.84, profile curvature higher than 0.0029, convex plan curvature class, and flat (F) slope aspect had the most influence on LSS for each corresponding factor.Land cover analysis results indicated that forest and urban classes had the highest probability of land subsidence occurrence, with FR values of 1.17 and 1.04, respectively.Distance to a stream of between 50 to 100 m had the highest FR, and the class of lower than 50 m had a considerable correlation; a stream density higher than 2.68 had the highest correlation with land subsidence and the 1.23 to 1.92 class, which also had a considerable FR value.For distance to road, the 0 to 100 m class had the highest FR followed by the 100 to 200 m classes.Finally, groundwater drawdown ranging from 55 to 83 m and from 28 to 55 m had higher impacts on land subsidence occurrence.

Result of Hybrid Models
In the course of implementation of the hybrid models, 70% of the land subsidence points (210 locations) were used for training with values 1, and the same number of randomly selected non-subsidence points were taken into account with 0 values in the training phase.For the test dataset, 30% of the subsidence inventory (90 locations) with a value of 1 was used, with 90 randomly assigned points with a value of 0. The training datasets were used to calibrate the weights of the membership functions.The testing dataset was used to evaluate the performance of the trained ANFIS ensemble models.Hybrid models were implemented in MATLAB 2017b software.The parameters used in meta-heuristic algorithms are presented in Table 4.The prediction power of ANFIS and the two hybrid models with the training dataset (target) along with the comparison of the output and target testing dataset is shown in Figure 6.The convergence results of the two ANFIS-ICA and ANFIS-GWO ensemble models up to 1000 iterations are shown in Figure 7. ANFIS-ICA had a better convergence value (0.276) than ANFIS-GWO (0.313).The lowest amount of the cost function (RMSE) indicates the best cost and thus the best performance in predicting the results.

LSSM using ANFIS and its optimized models
The original ANFIS and its optimized ensembles in this research were trained and used to estimate land subsidence susceptibility in the study area.Susceptibility modelling and estimation were all carried out in MATLAB 2017b and were then exported to ArcGIS 10.3 software to classify and generate the susceptibility maps.Land susceptibility index was classified into five classes, very high, high, moderate, low, and very low, based on a natural break classification scheme [22,41].Figure 8 presents the generated classified subsidence susceptibility maps obtained from ANFIS, ANFIS-GWO, and ANFIS-ICA.As can be seen, all the output subsidence susceptibility maps are similar and consistent with each other, particularly the ones for ANFIS-ICA and ANFIS-GWO.Moreover, the map based on ANFIS-ICA is much smoother than the others.

LSSM Using ANFIS and Its Optimized Models
The original ANFIS and its optimized ensembles in this research were trained and used to estimate land subsidence susceptibility in the study area.Susceptibility modelling and estimation were all carried out in MATLAB 2017b and were then exported to ArcGIS 10.3 software to classify and generate the susceptibility maps.Land susceptibility index was classified into five classes, very high, high, moderate, low, and very low, based on a natural break classification scheme [22,41].Figure 8 presents the generated classified subsidence susceptibility maps obtained from ANFIS, ANFIS-GWO, and ANFIS-ICA.As can be seen, all the output subsidence susceptibility maps are similar and consistent with each other, particularly the ones for ANFIS-ICA and ANFIS-GWO.Moreover, the map based on ANFIS-ICA is much smoother than the others.

Validation
The ROC curves were calculated for all LSS maps using the test data.Figure 9 demonstrates the comparison of AUC for all the models used.The results showed that the ANFIS-ICA had the highest prediction accuracy (0.932), followed by the ANFIS-GWO (0.926) and ANFIS (0.908).This proves that the combination of the ANFIS model with meta-heuristic algorithms such as GWO and ICA can significantly improve the output land subsidence susceptibility maps in comparison to ANFIS alone.The ROC curves were calculated for all LSS maps using the test data.Figure 9 demonstrates the comparison of AUC for all the models used.The results showed that the AN-FIS-ICA had the highest prediction accuracy (0.932), followed by the ANFIS-GWO (0.926) and ANFIS (0.908).This proves that the combination of the ANFIS model with meta-heuristic algorithms such as GWO and ICA can significantly improve the output land subsidence susceptibility maps in comparison to ANFIS alone.

Discussions
Land subsidence is the slow vertical lowering of the Earth's surface, posing a serious threat to both the environment and human life.Recently, there has been an increasing interest in land subsidence analysis and monitoring in Iran as it is one of the highest subsidence-prone countries [17,18,20,24].Natural hazard vulnerability analysis using machine learning algorithms (i.e., ANFIS) has shown promising results.Therefore, in this research, the focus was to employ the ANFIS model in combination with meta-heuristics in land subsidence susceptibility mapping.
Land subsidence inventories are necessary for accurate subsidence susceptibility analysis.The use of remote sensing SAR data is suitable for providing subsidence inventories due to their wide availability, independence from fieldwork, time and cost efficiency, frequent repeatability over time, and, especially, high precision [6].In this work, the PS-InSAR technique with its millimetric precision was employed to determine the subsided areas in the region of interest to form the inventory data used for training and testing the LSS models.

Discussion
Land subsidence is the slow vertical lowering of the Earth's surface, posing a serious threat to both the environment and human life.Recently, there has been an increasing interest in land subsidence analysis and monitoring in Iran as it is one of the highest subsidence-prone countries [17,18,20,24].Natural hazard vulnerability analysis using machine learning algorithms (i.e., ANFIS) has shown promising results.Therefore, in this research, the focus was to employ the ANFIS model in combination with meta-heuristics in land subsidence susceptibility mapping.
Land subsidence inventories are necessary for accurate subsidence susceptibility analysis.The use of remote sensing SAR data is suitable for providing subsidence inventories due to their wide availability, independence from fieldwork, time and cost efficiency, frequent repeatability over time, and, especially, high precision [6].In this work, the PS-InSAR technique with its millimetric precision was employed to determine the subsided areas in the region of interest to form the inventory data used for training and testing the LSS models.
Important conditioning factors for determining land subsidence prone areas were identified and collected based on either the literature or availability of data.The FR model was used to evaluate the correlation and influence of the factors.The results showed that all the factors employed in this paper had a considerable effect on LSS in Shahryar County.Among all the factors, the flat (slope aspect) area had the highest FR value (3.7), indicating high subsidence susceptibility in flat areas.The slope angle is related to the hydro physiographic characteristics that can influence the water infiltration rate and the volume and velocity of the Earth's surface flow [13].Altitude and groundwater drawdown were the best predictors of land subsidence in this study, followed by stream density and distance to stream.Rahmati et al. and Arabameri et al. [2,37] also found that groundwater drawdown had a greater impact on land subsidence.Other topo-hydrographic factors, such as stream density and distance to stream are indirectly related to LS as they impact groundwater recharge and infiltration [2,40] and, as can be seen in the results, the land areas closer to streams and with a certain stream density were more susceptible to LS.In terms of altitude, the lower lands were more prone to subsidence as the class 1137 to 1119 m and those lower than 1119 m had the highest FR.TWI, plan curvature, and profile curvature are among secondary topographic derivatives indirectly influencing LS [2,40,71].These factors were not among best predictors of subsidence in the study area, which may be due to smooth and low altitude changes in the study area.The FR analysis showed that a TWI lower than 4.84 was strongly correlated with LS.In a similar study [40], lower TWI values (i.e., 2.54 to 8) have been reported to be more prone to subsidence.Positive and convex plan and profile curvatures had the highest FR value, as reported in [40].Cropland and urban land cover types exist in lower altitude and flat areas.The main water source of the area is groundwater; therefore, more extraction of water in recent years as a result of population increase has caused the subsidence rate to increase.Previous studies have stressed the impact of groundwater extraction on subsidence occurrence [72,73].Regarding the distance to road factor, the closer to a road, the greater the land subsidence risk, which can be due to closeness to urban land cover and thus indirectly related to the subsidence phenomenon.
Two novel meta-heuristic algorithms, GWO and ICA, were used to optimize the rules and parameters of the ANFIS model.Both of these evolutionary algorithms belong to swarm intelligence.The results showed that ICA had a slower convergence rate than GWO; however, it had better performance.In order to evaluate the prediction power and accuracy of the models, RMSE and ROC criteria were used.The RMSE is simply based on error assessment, whereas ROC is based on true positive (TP), false positive (FP), true negative (TN), and false negative (FN), which is more appropriate for comparison [42].ANFIS-ICA had the lowest RMSE in both the training (0.276) and testing (0.3199) phases, followed by ANFIS-GWO and ANFIS alone.According to the AUC-ROC results, the ANFIS-ICA model was more accurate (0.932), followed by the ANFIS-GWO model (0.926) and the ANFIS model (0.908).It can be seen that the use of machine learning algorithms resulted in higher prediction accuracy since ANFIS alone yielded a suitable performance compared to other statistical methods in other studies.It can also be concluded that optimization of the ANFIS algorithm by meta-heuristics improves its results considerably.This was also reported in cases of other applications [64,74].The results showed that the ICA algorithm was more accurate than the GWO algorithm in optimization of the ANFIS model.The advantages of the ICA algorithm are high convergence speed and the ability to optimize functions with a large number of variables [75].The GWO algorithm has a small number of disadvantages, including a low solving accuracy, poor local searching ability, and slow convergence rate [60].
The output land subsidence susceptibility maps of the three models were similar and in line with each other.However, the map produced by the ANFIS-ICA was smoother than that of the other two.As could be observed, the high-risk areas were predicted where the groundwater extraction was higher, elevation was lower, and agricultural land use was higher.This is because the main source of income in the study area is agriculture.Further, the population has increased; therefore, food production has stressed the groundwater, the main water source of the area, and thus the land subsidence risk has become higher in those areas.Further, it is evident that the subsidence trend is gradually reaching towards the urban part of the Shahryar County, posing a serious threat to settlements and human life.The generated LSSMs in this paper can benefit authorities and decision-makers to identify subsidence-prone areas regarding environmental and urban management.

Conclusions
Land subsidence is an important issue in Iran due to the semi-arid and arid climate and excessive groundwater extraction.Therefore, modeling, simulation, and risk mapping offer valuable knowledge of environmental geohazards.GIS-based predictions have proved to be essential for authorities in terms of planning and decision-making.In this work, we used remote sensing SAR data and the PSInSAR technique to create a land subsidence inventory of the study area as a high-precision tool with a low cost and frequent reproducibility.Since machine learning tools have shown appropriate performance in modeling and mapping hazard susceptibility, the ANFIS model was used in this research to map the land subsidence risk in Shahryar County, Tehran province, Iran.Another objective of this paper was to investigate the effect of optimization of the ANFIS model through meta-heuristics.Two novel evolutionary algorithms, namely, GWO and ICA, were used to create ensemble models.The results of the three models in both training and testing phases were assessed by RMSE.In both phases, ANFIS-ICA had the lowest RMSE, followed by ANFIS-GWO and ANFIS alone.AUC-ROC analysis was also used for model evaluation, and its results indicated that ANFIS-ICA had the best prediction performance (0.932), followed by ANFIS-GWO (0.926) and ANFIS (0.908).To conclude, the results overall showed the applicability of the ANFIS machine learning algorithm in land subsidence susceptibility mapping and the effectiveness of its ensembles with metaheuristic algorithms.The methodology used is reproducible and can be applied to other regions with different environmental parameters to test the modelling performance.Further studies should be applied using other machine learning and deep learning algorithms to compare their prediction accuracy.In addition, future research can focus on developing risk monitoring and early-warning frameworks.

Figure 1 .
Figure 1.Flowchart of the overall methodology.Figure 1. Flowchart of the overall methodology.

Figure 1 .
Figure 1.Flowchart of the overall methodology.Figure 1. Flowchart of the overall methodology.

Figure 2 .
Figure 2. The location of the study area along with the extracted land subsidence inventory.

Figure 2 .
Figure 2. The location of the study area along with the extracted land subsidence inventory.

→A > 1 ,
it means the grey wolves are split and distributed in diverse ways for searching of the prey.After finding it, they congregate to attack.The coefficient → C provides a random weight for the prey while → C > 1 and promotes the exploration phase.In addition, → C models the natural hindrances in hunting for the grey wolves.
Figure 5(a) shows the LOS deformation map.According to Figure 5(a), the maximum velocity was about -175 mm/year, which occurred in the southern agricultural part of the region of interest, where more ground-water extraction was observed.

Figure 4 .
Figure 4.The perpendicular baseline graph for the time-series data stack.The dots represents the images, and the edges denote the interferograms.

Figure 4 .Figure 5 .
Figure 4.The perpendicular baseline graph for the time-series data stack.The dots represents the images, and the edges denote the interferograms.

Figure 5 .
Figure 5. (a) The annual velocity map based on the PSInSAR technique.The map is superimposed on the Google Earth (GE) imagery; (b) IDW interpolation raster of the velocity map.

Table 4 .The
Parameters used in hybrid algorithms.RMSEs of the training and testing phases were calculated and are shown in Table5.The two ensemble models enhanced the ANFIS model, and the ANFIS-ICA outperformed the ANFIS-GWO with an RMSE of 0.276 in the training phase and 0.3199 in the validation and testing phase.The ANFIS-GWO yielded an RMSE of 0.313 and 0.3217 in training and validation phases, respectively.Finally, the ANFIS model resulted in 0.323 in training and 0.34 in the validation phase.The convergence results of the two ANFIS-ICA and ANFIS-GWO ensemble models up to 1000 iterations are shown in Figure7.ANFIS-ICA had a better convergence value (0.276) than ANFIS-GWO (0.313).The lowest amount of the cost function (RMSE) indicates the best cost and thus the best performance in predicting the results.

Figure 6 .
Figure 6.The target and output values for training and validation datasets of (a) ANFIS, (b) ANFIS-GWO, and (c) ANFIS-ICA.The RMSEs of the training and testing phases were calculated and are shown in Table 5.The two ensemble models enhanced the ANFIS model, and the ANFIS-ICA outperformed the ANFIS-GWO with an RMSE of 0.276 in the training phase and 0.3199 in the validation and testing phase.The ANFIS-GWO yielded an RMSE of 0.313 and 0.3217 in training and validation phases, respectively.Finally, the ANFIS model resulted in 0.323 in training and 0.34 in the validation phase.

Figure 7 .
Figure 7.The convergence graph of the objective functions.

Figure 7 .
Figure 7.The convergence graph of the objective functions.

Figure 9 .
Figure 9.The ROC curves for the LSSMs of the models and their AUC.

Figure 9 .
Figure 9.The ROC curves for the LSSMs of the models and their AUC.

Table 1 .
The details of the SAR data used.

Table 2 .
The details of the input conditioning factors.

Table 3 .
Relationship between land subsidence occurrence and conditioning factors using the FR model.

Table 5 .
The comparison of model performance

Table 5 .
The comparison of model performance.