MONITORING OF LAND SURFACE DISPLACEMENT BASED ON SBAS-INSAR TIME-SERIES AND GIS TECHNIQUES: A CASE STUDY OVER THE SHIRAZ METROPOLIS, IRAN

more, the amount of coherence was 0.8 in estimating displacements. Based on the results, the highest displacements occurred in districts 2, 4, 7, and 11. The maximum value of the displacement (30 mm per year) was seen in the southeast of the Shiraz metropolis, in district 7. Moreover, the surface displacement was mainly observed in the lands of the Shiraz international airport and some of the agricultural lands located in district 2. The reasons for the displacements in districts 4 and 11 are related to the high density of the building coverages in these regions. Also, the increased population density was another reason for the subsidence in district 4. It can be concluded that the suggested framework is efficient in determining time-series surface movements in extensive urban and non-urban areas.


INTRODUCTION
Subsidence is one of the types of natural and geological dangers that can happen naturally or human activities such as long-term extraction of subsurface water and traditional agriculture Irrigation patterns (Lall et al., 2020;Salehi et al., 2022).In the last two decades, due to the occurrence of climate changes and consecutive droughts, corrupt management of water resources, indiscriminate extraction of underground water, and increasing population growth, the land subsidence phenomenon has occurred in Fars province, especially in the southeast region (Farhadi and Najafzadeh, 2021;Najafzadeh et al., 2021).As an effective tool for studying all phenomena that cause changes to the land surface, interferometric synthetic aperture radar (InSAR) was proposed (Xu et al., 2019).InSAR-based radar remote sensing (RS) has been proposed as one of the geodetic techniques for identifying the altitude changes on the earth's surface, which is very interesting due to its advantages over other methods (Berardino et al., 2002;Lohman and Simons, 2005).This method is considered the most efficient technique among space methods for measuring the changes in the land surface with high spatial resolution and accuracy (Foroughnia et al., 2019;Khorrami et al., 2020;Khoshlahjeh Azar et al., 2021).Among the advantages of this method, we can mention very high accuracy, wide coverage, high spatial resolution, cost-effectiveness, and the possibility of accessing information in any weather condition (Seymour and Cumming, 1994;Tough et al., 1995;Biggs et al., 2007).
To describe the changes in land surface subsidence in the Yazd-Ardakan Plain between 2003 and 2020, Mirzadeh et al., (2021), applied the InSAR time-series technique and statistical models (Mirzadeh et al., 2021).Based on their study, it was discovered that a region of 234.45 square kilometers has surface subsidence at rates up to 15 cm per year in the northwest-southeast direction.Also, data have confirmed that groundwater levels have dropped by 18 m between the years 1974 and 2018 that driving the compression of sediments within the underlying aquifer system.Sorkhabi et al. (2022) have measured Isfahan land subsidence by 30 Sentinel-1 SAR data by using the InSAR technique from 2020 to 2021 (Sorkhabi et al., 2022).To assessment the suggested techniques, three levelling points were measured in two months of 2020 by Specific root-mean-square error (RMSE).The maximum amount of surface subsidence gets greater than 170 mm per year.The total average subsidence of the region of interest in Isfahan was 110.9 mm per year.Based on 89 Sentinel-1 images prepared between December 2017 and December 2020, Liao et al. (2021) proposed an improved time-series InSAR technique using Mintpy open-source software for experimental purposes.An inversion of the interferometric network was performed by using the weighted least squares (WLS) solutions in this research using the GAMMA software.They found that the subsidence funnels widened and land subsidence velocity increased between 2017 and 2020 (Liao et al., 2021).A significant amount of land subsidence was observed along Wuyue Plaza, Yanzhou Avenue, and in the Eastern Development Zone in YAND.Improved time-series InSAR (TS-InSAR) technique, detected maximum subsidence rates of 21.1mm per year, 37.5mm per year, and 50.2mm per year, respectively, from 2017 to 2018, 2018 to 2019, and 2019 to 2020.The final monitoring outcomes demonstrated that the YAND area, with deformation amount mostly in the range of −10 to 10 mm per year, is generally relatively stable (Liao et al., 2021).Nevertheless, the maximum subsidence rate of 100 mm per year can be observed in three meaningful subsidence funnels in the fill region.In another study, Zhang et al (2019) examined long-term monitoring of subsidence in Wuhan city using the SBAS-InSAR technique with high-resolution Radarsat-2 images and the reasons for subsidence (Zhang et al., 2019).For their work, they investigated the potentiality of 20 Radarsat-2 images for descending passes acquired between 17 October 2015 to 3 June 2018 to derive subsidence in Wuhan city.The InSAR results were validated by 56 levelling points.At the end of the study, they studied the interconnection between natural conditions and human activities in influencing land subsidence.Furthermore, Wu et al (2019) applied the SBAS-InSAR technique to delineate the subsidence process in Yan'an city.In total, 55 Sentinel-1A images were processed to obtain the land subsidence in Yan'an city and Yan'an new airport between 2015 and 2019, related before and after the earthwork (Wu et al., 2019).A ground-based measurement of the new Yan'an airport validated the results obtained.A subsidence rate of -70 to 30 mm per year was recorded in the YND, and -50 to 25 mm per year in the YNA.In the other research, Li et al. (2020) investigated the relationship between surface subsidence and groundwater changes.Between June 2003 and November 2015, 47 ENVISAT ASAR images and 48 RADARSAT-2 images (descending track) were used to acquire the land deformation of the Beijing Plain and then validated with levelling benchmark measurements (Li et al., 2020).Also, they innovatively estimated auxiliary stresses to obtain information on static and dynamic loads.Using extremely randomized trees (ERT) machine learning approach and spearman's rank correlation coefficient (SRCC) technique, they evaluated the contribution rate of the influencing factors to land subsidence.As a result, between 2003 and 2010, the maximum land subsidence rate reached 110.7 mm per year, and from 2010 to 2015, it reached 144.4 mm per year.Yao et al (2019) employed ENVISAT ASAR data from 2004 to 2010 and Sentinel-1A data for 2015 to 2017 in ascending orbit pass.They used the SBAS InSAR technique in Shanghai, China, to monitor land deformation variations on a spatial and temporal scale (Yao et al., 2019).Deformation rates in the urban area gradually decreased from 14 to 8 mm per year due to the gradual rebound.An analysis of the Spatio-temporal pattern of ground subsidence in Tehran was conducted by Alipour et al. (2008) (2009) used ENVISAT ASAR data by small spatial and temporal baselines to monitor subsidence in Neyshabour and LOS deformation computed from the time-series analysis shows a remarkable subsidence rate of up to 190 mm per year.they confirmed their results with GPS measurements.
Previous studies show that the main focus was on just one track of sensors.In the current research, we used both the ascending and descending images to evaluate our estimations.Also, the subsidence phenomenon was investigated for 4 years which was a suitable period.Furthermore, estimation of the average annual velocity of subsidence has not been done yet for every district in the Shiraz metropolis.In this research, we employed the SBAS method to assess the subsidence in 10 districts of the Shiraz metropolis.Since this method mainly has some information gaps for subsidence in the orchards areas.Therefore, the novelty of our study was using the derived points from SBAS in the kriging method for the land surface subsidence estimation in these areas.The calculation of kernel estimation density (KED) for the production of the population map and built-up density map is one of the other innovations of this research.In addition, another novelty of this research was the production of fuzzy maps to decrease the uncertainty in built-up and density maps.
The main aims of this research are to (i) an estimation of the average annual velocity of subsidence between 2017 to 2021 using the SBAS technique, (ii) identification of the areas of exposure to the subsidence disaster in the different districts of the Shiraz, and (iii) investigate the reasons for subsidence in the different districts of the Shiraz.
In general, the current research organization begins with an overview of the study area in section 2. A detailed description of the methodology is provided in the section 3 of the paper.The results and discussion are then qualitative and quantitatively analysed in section 4. Ultimately, the major results are provided in the conclusions section.

Case Study:
Shiraz is one of Iran's five largest cities in Fars Province, the southwest of Iran.The city is geographically situated between the coordinates 29°53'59" N, 52°16'51" E, and 29°32'7" N, 52°39'54" E, and near Maharloo Lake.The height above sea level is 1488 meters in the east of the city and 1700 meters in the west.The city has been a regional trade center for over a thousand years and has a moderate climate.As the capital of Fars province, Shiraz is known for its literary history and many gardens.Nowadays, there is a subsidence phenomenon in some cities of Iran, especially in the Shiraz metropolis.The high population of the metropolis and the consequence of high-water consumption can be one of the reasons for the subsidence.Recently, many wells were excavated in and around the urban environments to harvest the water, which caused the plains and cities built on plains to encounter serious threats.The existence of this high population, water harvesting of wells, and the density of buildings in Shiraz make prone this metropolitan to subsidence that can threaten the citizens of this city (Sadeghi et al., 2015).The location of Shiraz city is presented in Figure 1.

Data:
Data from Sentinel-1 satellites can be acquired in four different acquisition modes: interferometric wide swath (IW), extra-wide swath (EW), strip map (SM), and wave mode (WV).The single look complex (SLC) products are composed of three sub-swaths on each Sentinel-1 satellite image.In this research, two sets of Sentinel-1 C-band data were obtained between February 2017 and August 2021, in the ascending orbits number 28 and, descending orbits number 137, (see Table 1 for detail), with the imaging mode IW swath mode, polarization mode of VV + VH, and spatial resolution of 5 m by 20 m (azimuth and range).Therefore, the total dataset consisted of 33 scenes captured by Sentinel-1A satellite, in ascending and descending orbits.For ascending orbit, the investigated region is located in the sub-swaths of IW2 and for descending orbit in the sub-swaths of IW1.In addition, the DEM of Shuttle Radar Topographic Mapping Mission (SRTM) with 30 m spatial resolution was used to remove the topographic undulation error in differential interferometry of two dense networks of SBAS interferometric pairs (Figure 2) with both perpendicular and temporal baselines in the data processing.To generate satellite image products, the processor requires additional data not included in the satellite images.Therefore, it was essential to retrieve auxiliary information from the https://scihub.copernicus.eu/gnss/#/home.The orbit files should correspond to the formerly downloaded SLC images for the time information.Multi-orbit datasets are available freely online, allowing inter-comparison of the monitoring outcomes to validate self-consistency.

Method:
The current research investigated the SBAS interferometry method using InSAR Scientific Computing Environment (ISCE) and Miami InSAR Time-series (MintPy) software.First, ISCE software was used for data pre-processing, and in the next step, MintPy software was used to produce the displacement map.The procedures for carrying out the proposed framework are depicted in Figure 3.According to the flowchart in Figure 3, the current research methodology is included in three parts.In the first stage, the data preprocessing process is accomplished by ISCE software.In this stage, Sentinel-1 satellite images, corresponding orbital files, SRTM DEM, number of connections (three for both passes), and the window dimensions of the multi-looks along distance and azimuth direction have entered this software as input.Then, image co-registration, interferogram generation, and phase unwrapping will be done.The interferograms were unwrapped, geocoded, and prepared for the next step.In the second stage, and after the preprocessing stages, the MintPy software we used.In this step, some reference points were identified for the estimation of displacements and corrected the grids of interferograms.Then, the network inversion stage was done by using the inversion of the interferograms network into time-series with weighted least square (WLS).After removing the unwrapping, atmospheric and topographic errors, velocity and time series maps were produced for two passes.The final step was the evaluation phase of the displacements for the estimations in both passes.Then the reasons for the subsidence were investigated in different districts using the maps of well, land cover, population density, and builtup.

SBAS-InSAR Data Processing Technique:
The SBAS approach is one of the most recent techniques used to monitor the expansion of land surface displacements over time.This technique was suggested by (Berardino et al., 2002).The SBAS InSAR technique minimizes topography errors and spatial decorrelation effects with interferograms with a small baseline, which is the primary benefit of the SBAS technique.In the other words, this method increases the temporal sampling rate and effectively overcomes the spatial uncorrelated problem caused by long spatial baselines and extracts the atmospheric delay effects as well as noise and topography components while removing the time-series deformation information of the ground surface and atmospheric delay effects as well as topography and noise (Berardino et al., 2002).
The following sequential acquisition process was followed for each SAR image.All interferometric pairs had a perpendicular baseline of lower than 100 meters.Afterwards, a multi-looking technique was performed with 40 pixels in the range direction and 10 pixels in the azimuth direction to form a square pixel (200 m) and reduce the speckle effect.Based on 33 radar images in ascending and descending orbit, a total of 95 interferograms have been computed.At first, the SRTM DEM was used to simulate the topographic phase and then subtracted from the interferogram.An unwrapping method based on the minimum cost flow (Snaphu) was used in the current research to unwrap the differential phase.Implementation of these processes was conducted by ISCE software (https://github.com/isceframework/isce2,accessed on 5 May 2022).Therefore, the interferometric processing conducted by ISCE software is orbit correction, de-burst, co-registration, generation of interferogram and applying the adaptive filtering, topographic phase removal using the DEM, and 2D unwrapping technique.
In time series land deformation estimation, the line-of-sight (LOS) displacement are calculated by inverting the unwrapped interferogram network.In this research, we use the WLS estimator to convert the network of interferograms into time series.The WLS estimator was used to minimize the phase residual for the time series raw phase.For ascending datasets, the image of 2019/04/08, and for descending datasets, the image of 2019/4/18 was processed as reference data.The spatial differences in land surface deformation were calculated based on pixels with high temporal coherence more than 0.8.Then, to reduce residual errors in the phase of deformation, tropospheric delay correction, phase deramping, and topographic residual correction were used.By using reliable pixels at each acquisition, linear phase ramps, which might be caused by the residual troposphere and ionosphere delays, were calculated and eliminated from the displacement time series using Python-based atmospheric phase screen estimation (PyAPS) software.Based on the proportionality of the perpendicular baseline time series, we calculated the systematic topographic phase residual caused by the DEM error.This processing was carried out using MintPy software.(https://github.com/insarlab/MintPy,accessed on 5 May 2022).

Data derived from GIS:
In this study, we prepared the wells, population density, and built-up density layers into a point layer in ArcMap software.The wells layer was used for visual interpretation and overlay on Sentinel-2 false color composite imagery with the subsidence map.The population and building density layers were converted to raster layers using KED (Węglarczyk, 2018).The KDE estimates the probability density function (PDF) of a random variable in a non-parametric manner.Raster-prepared layers (population and buildings) were converted into fuzzy maps using the fuzzy method of MS Large.In this step, values from 0 to 1 are assigned, where 0 is unlikely or not the target and 1 is very likely or the target.Therefore, a higher fuzzy membership value indicates a better target (the area affected by the subsidence disaster in this study).Understanding the membership types is important when assigning fuzzy membership values (Ahamed et al., 2000).As follows are the types of memberships: • Fuzzy Linear Membership (FLM)-High fuzzy membership (HFM) is assigned to large or small values and the fuzzy membership decreases at a constant rate.• Fuzzy Small Membership (FSM)-The HFM is assigned to small values.• Fuzzy Large Membership (FLM)-The HFM is assigned to large values.• Fuzzy MS Small Membership (FMSSM)-The HFM is assigned to values smaller than the mean.• Fuzzy MS Large Membership (FMSLM)-The HFM is assigned to values larger than the mean.• Fuzzy Near Membership (FNM)-Values in the middle range are assigned to the HFM.
In the current research, the Fuzzy MS Large (FMSL) membership function is used for the population and building density layers.
The FMSL transformation function is similar to the FMSL function, except that the definition of the function is based on a specific mean and standard deviation (Su et al., 2012;Farhadi et al., 2022).In general, the difference between the two functions is that the FMSL function is more applicable when a large value is more likely to be a set member.The result can be similar to the large function depending on the mean and standard deviation.By using a function of mean and standard deviation, a fuzzy membership is defined, where larger values show a closer membership to 1. Depending on the product of a×m, there are two Equation (1) and Equation (2) for FMSL.
Where m, s, a, b, x, and u(X) indicated the mean, standard deviation, coefficient of the mean, coefficient of the standard deviation, fuzzy membership, and fuzzy membership function, respectively.The FMSL is useful when there is a higher membership ratio for large input values.Input values may be integers or floating-point positive values.
After fuzzification of the population and building layers, we used the ordinary kriging method for estimating small sites (mainly orchards in the Ghasreh Dasht) that has not been calculated in the SBAS method to produce the final map of land subsidence in the two-orbit path.Due to the reduction of calculation and the stability of the candidate Permanent Scatterers (PSs) in the urban region as well as the reduction of the phase bias effect, we used long-term interferograms (Maghsoudi et al., 2022).For the points calculated by the SBAS method, we fitted the appropriate model based on the default settings of ArcMap.The root means square error (RMSE) of the Kriging method was very small, about 1.Then, the gap points of the SBAS technique were filled and a continuous map of the ground subsidence was created.Finally, the maps of wells, population, and development density were overlaid with the map of ground subsidence in ArcMap and Google Earth software, and then evaluated by visual interpretation.

Result and Discussion
The current study uses the SBAS technique to perform a continuous time-series examination, which helps to reduce undesired phase elements from differential interferograms.For the estimation of ground displacements in the Shiraz metropolis, Sentinel-1A data were collected along ascending and descending orbits between 2017 and 2021 (4 years).The location and extent of the observed subsidence between the acquisition times of the radar images are shown in Figure 4.The interferograms from both ascending and descending orbits show a clear deformation pattern.Furthermore, the generated interferograms thus visually confirm that the ground subsidence occurred between 2017 and 2021.Based on Figure 4, for both ascending and descending orbits, the maximum land displacement occurred in the southeast of the Shiraz metropolis, District 7. The displacements were mainly in districts 2, 4, 7, and 11.The pattern of subsidence was seen along the west to east (Figures 4 and 5) and increased in the east (District 7) and south (District 2), where there were irrigated agricultural lands and many wells (Figure 5).As seen in Figure 5, in districts 2 and 7, there are agricultural lands that are irrigated with water from the wells.Based on Figure 7, District 4 has the maximum population density in the shiraz metropolis.Therefore, the existence of high population density is another challenge in identifying areas exposed to subsidence in the Shiraz metropolis.Moreover, Figure 7 shows the fuzzy kernel memberships of population density.The highest values of fuzzy memberships (0.74) are related to the center and around it (districts 3, 8, 4, 5, 11).The values of population density were 123, 100, 95, 94, 91, 84, 81, 71, 66, and 57 for the districts of 4, 5, 2, 8, 11, 7, 3, 1, 6, 9, respectively.

2.2.1
Reasons Assessment for the Land Subsidence: Wu et al. (2019) observed that many residential and commercial buildings are built-up along the subsidence zone.Their results showed that human activities, including land creation and urban construction, might be the primary human factors contributing to severe local subsidence (Wu et al., 2019).Therefore, it seems the reason for the subsidence in district 4 is related to the high density of built-up (fuzzy membership>0.5) in these areas (Figure 6) and also the harvesting of groundwater (see wells in Figure 5) and high population density.Anjasmara et al. (2020) found that most of the land subsidence in Surabaya city has occurred in the industrial and residential areas (Anjasmara et al., 2020).A majority of Surabaya's residences are situated near mangroves and the coast.Zhou et al., 2020).They found many water-consuming plants and crops planted in these areas.Also, they showed a large amount of groundwater overexploitation resulted in severe surface subsidence in the study area.These researchers pointed out the dense population as well as agricultural and industrial lands as the factors that had many effects on increasing demand for water and further aggravation of land subsidence.Meanwhile, Shiraz international airport in district 2 and the active mines in district 7 seem to be other reasons for the subsidence in these areas.The density of the built-up and population in the districts also was investigated in this study.Our results about the fuzzy built-up density map showed that fuzzy memberships had values between 0 to 0.85, and were high in the center, east, and west of the shiraz metropolis.

Evaluation of the Result:
Since in the studied period (2017 to 2021) there is no levelling or GPS data present in the studied area, it is impossible to carry out an absolute assessment using ground reality data.Therefore, to properly evaluate the processing results in this region, the relative validation method was applied.So that the radar interferometric processing was done using a path with various passes.It should be noted that, if the results of the two paths (descending and ascending) confirm each other, the process will be relatively valid (Tough et al., 1995;Biggs et al., 2007).Therefore, the results showed that the maximum displacement was approximately 30 mm in descending and 25 mm in ascending orbits (Figure 4).According to the obtained results in the two paths, the results confirm each other relatively.

CONCLUSION
In the current research, we investigated the phenomenon of ground subsidence in the Shiraz metropolis between 2017 and 2021 using the SBAS InSAR method.The aim of this research was to an estimation of the average annual velocity of subsidence and determine the vulnerability to land subsidence for the phases of mitigation/prevention in crisis management between 2017 and 2021 using the SBAS method in the different districts of Shiraz.For this reason, to detect the deformation rate due to rapid urban growth in urban areas, the SBAS technique was used in the current investigation.However, because of the decorrelation effect, the SBAS technique is not effective in densely vegetated areas.Since the SBAS method can estimate millions of points for ground subsidence, these data can be used as input to the kriging method to estimate the unknown points, especially in the orchards where the SBAS method could not estimate.Also, in this research, to reduce the uncertainty of some layers such as population density and building density applied the fuzzy layers to identify the reasons for the surface subsidence.Based on the result, the most severe surface subsidence occurred in the west, east, and center of Shiraz, i.e., districts 7, 2, 4, and 11.Moreover, land subsidence was maximum in district 7. Our maps show that land subsidence occurred in agricultural areas (District 7), residential areas (districts 4 and 11), and near Shiraz International Airport in district 2. Visual interpretation of the pattern of ground subsidence in the metropolis of Shiraz showed a pattern from west to east.The high population density and a large number of built-up areas in districts 4 and 11, as well as some wells for water consumption in these regions, were the reasons for the ground subsidence in these districts.The extraction of water for agricultural lands and mining of some mines in district 7 can be other reasons for the land subsidence in these areas.Shiraz international airport was also a cause of the land subsidence in some areas of district 2. We concluded that the LOS displacement in Shiraz from 2017 to 2021 was nearly -30 mm per year.The results of this study have many benefits for the administrators and urban planners.Furthermore, the proposed framework in this research can also be used for estimating ground subsidence and its causes in other metropolitan cities.According to the main findings of the study, the integration of remote sensing techniques and GIS is very suitable for investigating land subsidence and its causes.Ultimately, it is suggested that, the use of the different theatrical semi-variogram models as well as different Kriging methods to interpolate the lost SBAS points in future research.

Figure 1 .
Figure 1.The geographical location of Shiraz in Iran and Fars province, and districts of 1 to 9, and 11 in this city.

Figure 2 .
Figure 2. The SBAS networks of the used datasets.(a) Path 28, (b) Path137.Each black point represents one image of data.The connected graph between every two points depicts one interferometric pair, whose color demonstrates the moderate coherence of the corresponding interferogram.

Figure
Figure 3. SBAS processing framework

Figure 4 .
Figure 4.The average annual velocity of the subsidence in 10 districts of Shiraz, is derived from the descending (a) and ascending (b) orbits; The reference point coordinate is 29.5516, 52.166.

Figure 5 .Figure 6 .
Figure 5. Agricultural lands and wells (green points) in district 7 and other districts

Figure 7 .
Figure 7.The fuzzy density of population in districts In addition, Zhou et al. (2020) investigated the subsidence associated with different land-use types in Beijing-Tianjin-Hebei and the use of water resources(Zhou et al., 2020).They found many water-consuming plants and crops planted in these areas.Also, they showed a large amount of groundwater overexploitation resulted in severe surface subsidence in the study area.These researchers pointed out the dense population as well as agricultural and industrial lands as the factors that had many effects on increasing demand for water and further aggravation of land subsidence.Meanwhile, Shiraz international airport in district 2 and the active mines in district 7 seem to be other reasons for the subsidence in these areas.The density of the built-up and population in the districts also was investigated in this study.Our results about the fuzzy built-up density map showed that fuzzy memberships had values between 0 to 0.85, and were high in the center, east, and west of the shiraz metropolis.

Table 1 .
Sentinel-1A data characteristic on descending and ascending orbits.