Analysis of Land Deformation using Small Baseline Subset (SBAS) INSAR Technique in Pokhara, Nepal

---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- Abstract Pokhara Valley lies in the geographic region rich in karst and the increase in urbanization with natural processes like sink-holes has led to the deformation of land causing socio-economic and physical effects on the people and the geography. Likewise, minimal research using the InSAR approach has been done in the context of the Pokhara Valley. So, this study focuses on the determination of land deformation in Pokhara Metropolitan from 2017 to 2022 using the SBAS technique. Sentinel 1A SLC images were adopted using ascending data acquisition mode. The result shows a cumulative displacement of 1578.97 mm to -1052.65 mm in ascending mode. For validation, SBAS-generated data points of the WRC regions were validated with GNSS value. The Pearson R value was found significant at a 5% significance level and the RMSD value was between 0.8 to 1.73. Thus, the results and the validation process indicate the suitability of SBAS for the determination of land displacement using SBAS.


Introduction
The geological study of the Pokhara Valley highlights the fact that most of the areas of the valley lie above the Ghachwok formation which is rich in calcareous sediments making it prone to karstification and eventually leading to land subsidence and sinkhole formation [1].Furthermore, the valley possesses a low bearing capacity and rapid urbanization growth, changes in land uses, and haphazard building construction further increase the risks associated with land subsidence [1,2].The fragile geology of Pokhara Valley makes this region susceptible to land deformation processes, including subsidence, uplift, and lateral movements.These hazards greatly impact the local environment, infrastructure, and community well-being, causing damage to structures, increasing flood and landslide risks, disrupting drainage systems, and affecting public health and sanitation [3].There is limited research exploring land deformation in this region.Traditional techniques such as spirit-leveling, extensometers, and continuous GPS measurements offer high precision but have limited spatial coverage [4].In contrast, space-borne instruments such as Interferometric Synthetic Aperture Radar (InSAR) offer comprehensive coverage and continuous measurements, serving as a cost-effective and efficient method for monitoring and analyzing surface subsidence.These advantages make InSAR widely preferred for regional subsidence monitoring purposes.InSAR utilizes satellite data to accurately track and analyze surface deformation over large areas, offering a valuable tool for assessing and managing subsidence [3,5].Differential Synthetic Aperture Radar Interferometry (DInSAR), Ground-Based Synthetic Aperture Radar (GBInSAR), Persistent Scatterer Interferometric Synthetic Aperture Radar (PSInSAR), Multi-Temporal Interferometric Synthetic Aperture Radar (MTInSAR), QuasiPersistent Scatterer Interferometric Synthetic Aperture Radar (QPSInSAR) and Small BAseline Subset (SBAS) are the deformation monitoring tool [6].They each have both advantages and disadvantages, however, PSInSAR, GBInSAR, MTInSAR, and DInSAR techniques were quite commonly used for deformation studies [6].The SBAS technique links independent unwrapped interferograms in time using distributed scatterers and singular value decomposition [7].The SBAS approach has been extensively employed in many geodetic domains, including ground subsidence, landslides, and seismic activity, due to its capacity to monitor large-scale deformation with millimeter accuracy [8].In China, SBAS was used to measure the land subsidence in Beijing, and further validation with leveling data the mean difference was found to be -1.6 mm/year with a standard deviation of 2.5 mm/year [9].Similarly, SBAS-based approach data was quantitatively assessed for the measurement of land subsidence in Napoli Bay (ITALY) and Los Angeles (USA), results show the accuracy of sub-centimeter and standard deviation of about 5 mm [10].SBAS was also used for processing the Sentinel-1A images for the identification of causes of the land subsidence [11].This study was conducted in Konya, Turkey which has a geological structure prone to the occurrence of land subsidence and sinkhole formation.Notably [12] used the SBAS InSAR technique for monitoring deformation velocity in both ascending and descending modes.Not only this, [13] also reported that the SBAS-based measurement generally coincides with leveling data though some difference between these methods still exists.This finding was further supported by the study conducted in Wuhan City to measure the land subsidence which resulted in lower variability of the deformation velocity when compared with ground leveling points [14].The validation of the time series SAR data with the available global navigation satellite system (GNSS) was used by [11].Further, [15] has observed the correlation coefficient value of 0.96 and RMSE value of 5.9 mm/year under validation of SBAS-based land subsidence with ground levelling points.This indicates the reliability of SBAS in measuring land deformation.Hence, the SBAS approach is applicable for the longtime monitoring of land subsidence in larger areas with lower variability.With ongoing urbanization in Pokhara, it is crucial to investigate the relationship between urban development and land deformation.This research aims to provide valuable insights into the dynamics of land deformation contributing to a better understanding of the associated risks and supporting sustainable development practices in the region.

Research Location
The research site for this study is Pokhara, where a comprehensive analysis of land deformation was conducted.Pokhara Valley is located at 28º 12'30''N and 83 º 59'20''E with the Seti River flowing through the valley.The geological characteristics reveal that Pokhara Valley is abundant in calcareous sediments [1].This composition renders the area susceptible to karstification, eventually resulting in land subsidence and the creation of sinkholes [1,2].The occurrence of sinkholes in 2013 serves as a significant event to examine the subsequent impacts on land deformation.

Methodological Framework
The fundamental methodological framework for deformation analysis is based on SARPROZ which is presented in the provided Figure 2. According to the figure, the methods adopted are categorized into distinct phases, such as:

Phase I: Data Preparation
In this study, we used Ascending track SAR acquisitions of Sentinel 1-A from 2/14/2017 TO 12/21/2022.
Sub-swath 3 with VV polarization was selected.The Master image of the date 18/01/2020 was selected for ascending.the rest were considered slave images.The selected Master area has a latitude of 28.2245 and a longitude of 83.9550 with a radius of 15km.The master area was selected such that it represents the area of Pokhara Valley.ii.Selection of Ground Control Points (GCPs): The task of selecting ground control points was performed to geocode the dataset using the Shuttle Radar Topography Mission (SRTM).

Phase III: InSAR Processing
Full graph coherence estimation and interferogram processing were done.For this process, a modified Goldstein filter with 15x15 coherence window size, Flattening, DEM Removal and unwrapping was done.
In this study, Small Temporal Baselines were used for image graph combination in the processing step for SBAS as shown in Figure 3.The parameter for selecting sparse points was "Amplitude Stability Index 1-Sigma/Mu.".The threshold is adjusted to optimize results with diverse PS points where the threshold value for the parameter was 0.7.Additionally, the downsample (DS) value was set at 10, and the downline (DL) value was set to 0. The Delaunay graph was generated to link different scattered points.
In the processing steps, parameters for the height and linear trend need to be set for the estimation of deformation, For SBAS linear trend from -150 to 150 was selected for the estimated deformation rate in the study area while the height parameter from -200 to 200 was considered for the average height of the study area to approximate calculation.Coherence as weight has been considered as the weight value for the SBAS.The second parameter, integrated residual height, arises from the disparity between integrated height and external DEM.The histogram of integrated residual height predominantly centers around zero.This distribution pattern signifies that a significant majority of points possess nearly identical relative heights in relation to the reference point.This outcome strongly indicates that the reference point is situated at ground level, as the prevalent scenario anticipates most points to be grounded.After the selection of PS, the Atmospheric disturbance and temporal coherence were evaluated, and to look at the fitness of the model.

PHASE IV: Validation with Continuous GNSS Observation
The GNSS data provided by the Septentrional Po-laRx5e receiver was compared with SAR-derived deformation measurements.The GNSS data acquisition was done from April 2022 to December 2022.Furthermore, Pearson's correlation coefficient (R) and Root mean square deviation (RMSD) of GNSS value and SBAS generated point were calculated for validation.An assessment was conducted to determine the reference point for linking scatterer points.The reference point, characterized as the most stable within the entire process, was identified.Integrated cumulative displacement and integrated residual height represent two crucial metrics that offer insights into the stability of our chosen reference PS point.An assessment was conducted to determine the reference point for linking scatterer points.The reference point, characterized as the most stable within the entire process, was identified.Integrated cumulative displacement and integrated residual height represent two crucial metrics that offer insights into the stability of our chosen reference PS point.Figure 6 The temporal coherence after removing APS and Temporal Coherence after Graph inversion and APS removal.

Reflectance Map
The presented reflectivity map, generated through SBAS in InSAR, shows the characteristics of the study area's radar backscatter.The intensity of radar backscatter is influenced by the surface properties of the observed area including the terrain type, vegetation cover and surface roughness.The brightness of the feature is dependent on the portion of transmitted energy that is the second parameter, integrated residual height, returned from the targets on the surface and the intensity of this backscatter energy.In the figure, brighter areas in the yellow circle indicate higher reflectivity while the dark portion in the red and blue lines represents lower reflectivity.The yellow region represents the urban areas while the red region represents lakes.Due to the surface roughness of the urban region and lake, it can be said that there is volumetric scattering in the yellow region and specular reflection in the dark region.The blue region represents the airport runway.As the runway's surface is smooth, there is specular reflection causing the area to appear dark.Similarly, a brighter portion is seen in the figure.This is because in SBAS more scatter points are formed in non-urban and vegetative areas.

SBAS-Based Land Deformation in Ascending Mode
The SBAS processing resulted in a cumulative displacement ranging from 1578.97 mm to -1052.65 mm.
The figure shows the cumulative displacement of the selected samples over the six years i.e., 2017 to 2022.On the x-axis, points are spaced at intervals of thirty, representing the displacement, while the corresponding displacement of each point is shown on the y-axis.Figure 10 The rate of displacement in ascending mode using SBAS in 2019 and 2020 Figure 11 The rate of displacement in ascending mode using SBAS in 2021 and 2022.

Deformation in ascending mode in a specific site
Furthermore, the data obtained through SBAS from 2017 to 2023 indicates distinct displacement patterns across several regions.The Armala area showed a gradual decrease in land displacement, evident in the declining red line on the graph.Lakeside experienced some displacement during the study period, while the Airport, WRC, and Basundhara regions showed relatively minimal displacement in comparison.
Figure 12 Land deformation in specific regions   The cumulative displacement of the different points taken from the period of 2017 to 2021 shows higher displacement, especially in the hilly areas with green vegetation and agricultural land.Furthermore, time series analysis of the cumulative displacement from the year 2017 to 2021 shows a similar displacement curve but varies in the range of displacement.Similarly, specific locations were analyzed using SBAS.The sample points are selected randomly in four main areas which are the airport, lakeside, WRC region, and Armala.The line of the airport (green line) showed minimal displacement very close to the reference line.This point was formed on the airport's runway.In the case of the WRC region, the upward displacement was observed because of the earth-filling process carried out in this region.However, the Armala region showed maximum downward displacement compared to the other three.Armala region has experienced sinkhole formation in the past few years which is one of the major reasons for the occurrence of the maximum downward displacement [16].In addition to this, the points of the Armala region consist of mostly agricultural fields and green vegetation while the settlement area is also gradually increasing which has resulted in downward displacement.

Validation with GNSS
For the validation, the SBAS displacement value is compared with the coarse station of GNSS.Similarly, the SBAS-generated data points also have a stronger relationship with the GNSS points with Pearson r value ranging between 0.61 to 0.69.[17]

Conclusion
In this study, multi-temporal InSAR observations to investigate the land surface deformation in Pokhara city were adopted.InSAR observations acquired by Sentinel-1 from 2017 to 2022 were processed and the time series analysis of land deformation was accomplished by adopting SBAS techniques.The cumulative displacement of the Armala, Lakeside, WRC and Airport was studied, and subsequently, higher downward deformation was observed in the Armala Region while minimal deformation was observed in the Airport.Thus, with the support of the validation results we can conclude that SBAS is suitable for land deformation analysis in areas covered with vegetation.However, the limitation of availability of an updated level network, the high level of specialized knowledge regarding the processing chain that is needed, selection and processing of the imagery are also a challenging step to get effective results.

Recommendations
This study provides land deformation trends over the years, which further supports planning in infrastructural activities, agriculture, and natural disaster occurrence.This study provides baseline information on the deformation for urban planners and agriculturalists.Furthermore, more academic studies focusing on the INSAR techniques in land deformation of other regions of Pokhara Valley are recommended.Similarly, the findings of this study give us an idea of the displacement of Pokhara Valley over the years and along with that it further supports the use of INSAR techniques in studies and urban planning.Future research on the establishment of GNSS reference points, and timely updating of the level network especially focusing on the Armala region is required.

Figure 1
Figure 1 Research Location 2.2 Data and Software 174 sentinel 1 SLC SAR images were acquired in ascending passes from 2/14/2017 to 12/21/2022.The images were downloaded from Sentinel-1 Scientific Data Hub.All the images used are in VV polarization mode.

Figure 2
Figure 2 Methodological Framework Based on SAPROZPhase II: Preliminary Analysis:A preliminary analysis was carried out, including the inclusion of the Reflectivity map, the Amplitude Stability Index, and the Mask of Sparse point selection.To obtain the temporal average of all the images in the collection, we built the Reflectivity Map, for which the temporal Standard Deviation has been analyzed.Mask generation involves creating a mask to identify sparse points.The chosen parameter for creating this mask was the extraction of local maxima from the reflectivity map.The following steps were performed for Preliminary Analysis i. Selection of External DEM: Shuttle Radar Topography was employed to select the DEM for the study area.

Figure 3
Figure 3 Different sets of graphs for Image Combination showing normal baseline vs Coherence in SBASAfter the connection processing with the following values, temporal coherence was estimated in the form of a histogram as shown in Figure4.The histogram ranges from 0 to 1.This indirectly reveals spatial coherence crucial for deformation estimation.This leads to the ultimate step: APS calculation.An assessment was conducted to determine the reference point for linking scatterer points.The reference point, characterized as the most stable within the entire process, was identified.Integrated cumulative displacement and integrated residual height represent two crucial metrics that offer insights into the stability of our chosen reference PS point.The zenith of the histogram for integrated cumulative displacement resides at zero, implying that a substantial portion of PS points-relative to the reference point-exhibit zero velocity.This observation aligns with the expectation that a considerable number of PS points should indeed possess a velocity of zero.

Figure 4
Figure 4 Temporal Coherence vs Connection graph

Figure 5
Figure 5 The integrated cumulative displacement graph and Integrated Residual Height graph of the selected Reference PS point in SBAS mode.

Figure 8 :
Figure 8: Cumulative Displacement using SBAS An analysis of the year-wise displacement of the land

Figure 9
Figure 9 The rate of displacement in ascending mode using SBAS.2017 and 2018

Figure 13
Figure 13 Graphical relationship of different SBAS generated data with GNSS value.The graphical relationship with the GNSS value and SBAS-generated data points was performed as shown in Figure.The scattered diagram shows a linear but inverse relation of SBAS generated value and GNSS value.The cumulative displacement of the different points taken from the period of 2017 to 2021 shows higher displacement, especially in the hilly areas with green vegetation and agricultural land.Furthermore, time series analysis of the cumulative displacement from the year 2017 to 2021 shows a similar displacement curve but varies in the range of displacement.Similarly, specific locations were analyzed using SBAS.The sample points are selected randomly in four main areas which are the airport, lakeside, WRC region, and Armala.The line of the airport (green line) showed minimal displacement very close to the reference line.This point was formed on the airport's runway.In the case of the WRC region, the upward displacement was observed because of the earth-filling process carried out in this region.However, the Armala region showed maximum downward displacement compared to the other three.Armala region has experienced sinkhole formation in the past few years which is one of the major reasons for the occurrence of the maximum downward displacement[16].In addition to this, the points of the Armala region consist of mostly agricultural

Table 1
Summary of Sentinel Data for Processing GNSS Data was acquired from the station located at Pashchimanchal Campus Pokhara.SARPROZ and QGIS software packages were used for processing and analyzing InSAR data.Excel and Google Colab were used for statistical analysis of data.
TableIIrepresents the results obtained from Bivariate Correlation and RMSE of the elevation of the different points obtained SBAS based technique with GNSS value.Similarly, for the SBAS-based measurement except for point 12913, all other points have a significant positive and higher correlation with GNSS value at a 5 % level of significance.

Table 2 .
Pearson's R and RMSE value for SBAS generated points with GNSS value