A Statistical Approach for the Integration of Multi-Temporal InSAR and GNSS-PPP Ground Deformation Measurements

Determining and monitoring ground deformations is critical for hazard management studies, especially in megacities, and these studies might help prevent future disaster conditions and save many lives. In recent years, the Golden Horn, located in the southeast of the European part of Istanbul within a UNESCO-protected region, has experienced significant changes and regional deformations linked to rapid population growth, infrastructure work, and tramway construction. In this study, we used Interferometric Synthetic Aperture Radar (InSAR) and Global Navigation Satellite System (GNSS) techniques to investigate the ground deformations along the Golden Horn coastlines. The investigated periods are between 2015 and 2020 and 2017 and 2020 for InSAR and GNSS, respectively. For the InSAR analyses, we used sequences of multi-temporal synthetic aperture radar (SAR) images collected by the Sentinel-1 and ALOS-2 satellites. The ground displacement products (i.e., time series and velocity maps) were then cross-compared with those achievable using the Precise Point Positioning (PPP) technique for the GNSS solutions, which can provide precise positions with a single receiver. In the proposed analysis, we compared the ground displacement velocities obtained by both methods by computing the standard deviations of the difference between the relevant observations considering a weighted least square estimation procedure. Additionally, we identified five circle buffers with different radii ranging between 50 m and 250 m for selecting the most appropriate coherent points to conduct the cross-comparison analysis. Moreover, a vertical displacement rate map was produced. The comparison of the vertical ground velocities derived from PPP and InSAR demonstrates that the PPP technique is valuable. For the coherent stations, the vertical displacement rates vary between −4.86 mm/yr and −23.58 mm/yr and −9.50 and −27.77 mm/yr for InSAR and GNSS, respectively.


Introduction
Istanbul is the province with the lowest per capita area in Turkey.With a population of more than 15 million, it is called Turkey's megacity.This city, which has a large population in a limited area, is vulnerable to natural and artificial threats.The most important of these are due to earthquakes, which are still in memory and cause severe socio-economic damage when they occur.In parallel with the city's development, local deformations that may increase the severity of possible disasters also have appeared.In the literature, Sensors 2024, 24, 43 3 of 16 The most suitable radius was determined by looking at the correlation values of the GNSS station with different radii.Additionally, one of the outlier detection methods was performed on the PS points included in the analysis with the aim to remove the inconsistent points.The weight of the points was calculated with the standard deviations of the obtained PS points.The weighted velocity estimation of the PS points was performed.Moreover, we calculated the vertical displacement rate of the stations based on the three InSAR datasets to determine the subsidence and uplift zones.The vertical displacement rates obtained from InSAR and GNSS methods were compared at the final stage.

Study Area
In this study, we focused on determining the ground deformations in the Golden Horn, one of the coastal regions of Istanbul.The city, with a population of over 15 million, plays a crucial role in culture, politics, and the economy.Specifically, the study area is the Golden Horn coastline, a natural harbor with a key position for transportation.Moreover, Golden Horn, known as the trading center of old Istanbul, is one of the most popular tourist areas of the city.The study area is also intertwined with the region under UNESCO protection.It is crucial to note that the coastlines of the Golden Horn do not have a solid geological structure due to the alluvial structure and artificial ground fillings.Hence, the region is vulnerable to slowly developing landslides due to its topography and geological structure and is in the southeast of the European part of Istanbul (Figure 1).Within the transportation network of Istanbul, Golden Horn has critical importance since the region serves as a hub for both sea and railway transportation, highlighting its strategic significance in facilitating urban connectivity.

SAR Data
Sentinel-1A and ALOS-2 SAR data were used in this study (Table 1).ALOS-2, operated by the Japan Space Agency (JAXA), is the latest developed and operating satellite among Japanese earth observation satellites.An active microwave radar (PALSAR-2) using the 1.2 GHz frequency range is mounted aboard this satellite [17].Sentinel-1 is one of Sentinel's mission satellites, which provides all-weather, day, and night radar imaging for land and ocean services at the C-band under the European Copernicus program.SAR images of Sentinel 1A were collected through both descending and ascending orbits.We used 28 ALOS-2 PALSAR-2 L band images obtained from ascending trajectories.In addition, we utilized 62 Sentinel-1 SAR data acquired with the TOPSAR Interferometric Wide (IW) Swath Mode in Single Look Complex (SLC) format from ascending orbits and 61 Sentinel-1 IW SLC data from descending orbits.For both ascending and descending datasets, VV polarization bands were processed.

SAR Data
Sentinel-1A and ALOS-2 SAR data were used in this study (Table 1).ALOS-2, operated by the Japan Space Agency (JAXA), is the latest developed and operating satellite among Japanese earth observation satellites.An active microwave radar (PALSAR-2) using the 1.2 GHz frequency range is mounted aboard this satellite [17].Sentinel-1 is one of Sentinel's mission satellites, which provides all-weather, day, and night radar imaging for land and ocean services at the C-band under the European Copernicus program.SAR images of Sentinel 1A were collected through both descending and ascending orbits.We used 28 ALOS-2 PALSAR-2 L band images obtained from ascending trajectories.In addition, we utilized 62 Sentinel-1 SAR data acquired with the TOPSAR Interferometric Wide (IW) Swath Mode in Single Look Complex (SLC) format from ascending orbits and 61 Sentinel-1 IW SLC data from descending orbits.For both ascending and descending datasets, VV polarization bands were processed.

GNSS Observations
The Golden Horn GNSS (GHGNSS) network was designed to determine the displacements along Golden Horn.GNSS measurements were started in August 2017 and repeated until March 2020.In August 2017, 14 stations were established according to the deformation areas given by Calò et al. [4] (Figure 2).However, due to the construction activities in the Golden Horn region, the GH12 station was re-established and named GX12.Moreover, two new stations named GH14 and GH17 were established after the fifth campaign.The session duration of the observations is about 6 h, with an interval of 10 s for all stations.The details of the campaigns are shown in Table 2.

Analysis of the GNSS Data
GNSS observations were processed based on the PPP technique, as stated before.This technique, which can provide position information with a single receiver, was first introduced at the end of the 1990s as an alternative Global Navigation Satellite Systems (GNSS)-based positioning technique [18].Several error sources and correction models should be considered to obtain precise positions when using PPP.Accurate position information and displacements can be obtained using precise orbit and clock products.Since this technique does not need a network structure and provides global reference frame solutions, local deformations can be determined with a single receiver.
PPP solutions were carried out using Trimble CenterPoint RTX online post-processing service.The service supports the processing of dual-frequency code and carrier phase observations of multi-GNSS systems based on the following equations [19].

Analysis of the GNSS Data
GNSS observations were processed based on the PPP technique, as stated before.This technique, which can provide position information with a single receiver, was first introduced at the end of the 1990s as an alternative Global Navigation Satellite Systems (GNSS)-based positioning technique [18].Several error sources and correction models should be considered to obtain precise positions when using PPP.Accurate position infor-Sensors 2024, 24, 43 6 of 16 mation and displacements can be obtained using precise orbit and clock products.Since this technique does not need a network structure and provides global reference frame solutions, local deformations can be determined with a single receiver.
PPP solutions were carried out using Trimble CenterPoint RTX online post-processing service.The service supports the processing of dual-frequency code and carrier phase observations of multi-GNSS systems based on the following equations [19].In the processing, we used both GPS and GLONASS observations.This service uses precise orbits and clocks produced by Trimble RTX.The update rate of the products is 1 Hz [20].For static positioning, the convergence time is 15 min [19].Moreover, linear combinations are used to eliminate ionospheric delays.Tropospheric delays and receiver clock errors are estimated in the processing step.In the service, a coordinate system can be selected by the user.We selected the ITRF2014 datum.The coordinates and their covariance matrix were obtained in a 3D Cartesian coordinate system.We then transformed each campaign's coordinates and covariance matrices to a topocentric coordinate system with respect to each station's first epoch.Moreover, a 6 h observation session duration, as indicated by Sezer et al. [20], employing the accuracy function provided by Saracoglu and Sanli [21], can yield millimeter-level accuracy.

Analysis of the SAR Data
Ground deformations in the Golden Horn were found using a multi-temporal InSAR analysis.SAR data processing was performed using SarProz software [22].Before beginning the image processing, a subset selection of the study area was carried out.Considering the definition of the orbit and atmospheric information, the image dated July 2017 was selected as the primary image.Then, the co-registration process was applied to the primary and secondary images.
At the stage of processing the ALOS-2 data, although the temporal and spatial bases of the images dated 5 October 2014, 14 December 2014, and 22 February 2015 were suitable, they were not included in the PSI analysis process due to low coherence and high level of errors in their trajectories.Provided analyses were carried out with 24 images.
Sentinel-1 images acquired at C-band on different orbits (descending and ascending) were obtained from the multi-temporal InSAR analysis conducted within the scope of the study.A set of coherent PS points were independently obtained after the PSI analysis was carried out using the collected SAR datasets.The ground deformations of these points were determined along the Line of Sight (LOS) direction.From the ALOS-2 dataset, 29,401 PS points were obtained.From the Sentinel-1 dataset, 33,501 and 33,552 PS points were obtained for ascending and descending trajectories, respectively.Figure 3 shows the spatial distribution of the PS points obtained for ALOS-2, Sentinel-1 descending, and ascending orbit.For all three datasets, PS points were obtained along with GH.The distribution of the PS points for the region was predominantly clustered in the urban area (road, buildings, port, bridge, etc.).For the whole coastal region of the GH, the mean displacement rates are −1.10 mm/yr, −3.82 mm/yr, and −5.08 mm/yr for ALOS-2, Sentinel-1 ascending, and

Statistical Evaluation of InSAR and GNSS-PPP Results
For the InSAR analysis, the ground displacements were obtained along the LOS direction.Calò et al. [4] employed one and a half years of TerraSAR-X data from November 2010 to June 2012, utilizing the Small Baseline Subset (SBAS) method.Their analysis revealed deformation in the coastal areas of the GH, resulting in a displacement of approximately 4-5 cm (an average of ~3 cm/year).Moreover, Imamoglu et al. [2] reported a displacement rate of 8 mm/year based on Sentinel-1 data collected between 2014 and 2017.Additionally, Aslan et al. [1] determined a maximum displacement of 10 mm/year using Sentinel-1 data for the same period.The present study corroborates these findings, indicating displacement values in the same region, particularly within the latter section of the GH.Furthermore, it is observed that the observed movement has persisted in this area after 2017.
From the GNSS measurements, the ground displacements are obtained in the topocentric coordinate system, namely, the east, north, and up directions.To compare the results obtained from InSAR and GNSS, we calculated the LOS displacements of the GNSS analysis by using the equation given by Hanssen [23] in Equation ( 3): where θ and α denote the angle of incidence of the radar wave and the direction angle of the satellite, U V , U N , and U E are the displacements obtained from the PPP solutions using the observations from the GNSS campaigns.To calculate the velocities of the stations for the InSAR and GNSS techniques, Equation (4) was used.
Here, LOS t and LOS 0 denote the displacements in the LOS direction in epoch t and the reference epoch, respectively, and v is the velocity.In the velocity estimation step, we also applied the Pope test to exclude outliers [24].
Before comparing the velocities, we checked the stations' velocities obtained by the GNSS analysis to see whether the movements were significant.In this step, considerable displacement was seen in nine of the seventeen stations.The subsequent analysis was performed for the nine stations.To compare the InSAR velocities and velocities obtained from the GNSS observations, we calculated the GHGNSS network stations' InSAR-based velocities using the PS points' velocities.Here, we created five different circle areas for each station.The centers of the circles were chosen as the GNSS stations, and the radii were determined as 50, 100, 150, 200, and 250 m.Then, we used the PS points within the determined circles.In Tables 3 and 4, the numbers of the PS points are shown for each station and each circle.As expected, and in the case of the ALOS-2 data (Table 3), the number of PS points increases when the radius increases for the Sentinel-1 data (Table 4).As shown in Tables 3 and 4, there are no PS points in some areas within the first and second-smallest circles, i.e., in the 50 m and 100 m radii.As a result, statistical analyses and estimates were carried out for the other circles (150, 200, and 250 m).Next, we calculated the correlations of the displacements of the PS points within the circles and the GNSS stations' displacements.The correlations between the ALOS-2 and GNSS data are up to values of 0.72, 0.84, and 0.78, respectively.The correlations between Sentinel-1 descending and GNSS are 0.65, 0.75, and 0.61; between Sentinel-1 ascending and GNSS, they are 0.65, 0.82, and 0.75.As can be seen from the correlation values, the best values were reached at a radius of 200 m for all three datasets.Hence, the statistical evaluations for this study will be on the results obtained in the 200 m radius since the highest correlation values were obtained in the 200 m radius circle.
The correlation values for the GH02, GH13, and IGNA stations were close to zero between the GNSS stations and PS points.Thus, these three stations were not included in the evaluation scope due to the errors caused by the observation periods and the possible problems that may have occurred in the campaign measurements.The sources of these errors will be examined in future research.
Under the assumption that the PS points within the circles reflect the GNSS stations' movements, we took the average of the PS points' displacements.We calculated the standard deviations for each station.However, some of the PS points within the circles could be outliers.Hence, we applied the median approach to the PS points' displacements to exclude outliers before taking averages in each epoch for each GNSS station.With a 50% breakdown point, the median approach is one of the most reliable outlier detection methods [20,25,26].Additionally, the median and median absolute deviation (MAD) are more efficient in outlier detection than the mean and standard deviation.The approach was applied using Equation (5).
where X denotes the displacements of the PS points.And n is the number of the PS points.In this approach, the |X − median(X)| values are compared with 3 × mad.
In the next step, the InSAR velocities of the GNSS stations for ALOS-2 and Sentinel-1 were estimated by least square estimation (LSE) based on Equation (4) [27].In the velocity estimation, we also considered the standard deviations of the displacements for each technique.Thus, we took into account the weights of each PS point.
As stated, six stations and three circles were analyzed using statistical methods.However, we simply reported the data of the 200 m radius circle for explanation purposes (Tables 5 and 6).In Table 5, the velocities obtained from the InSAR methods are shown.The GNSS velocities, transformed into LOS directions, have been presented in Table 6.As seen from Table 5, the velocities of the Sentinel-1 ascending and descending stations fit very well for most of the stations, and their formal errors are below 1 mm/yr.The correlation of the velocities is "0.96".Moreover, almost all the formal errors of the velocities from ALOS-2 are sub-mm/yr.The correlations of the velocities from ALOS-2 and Sentinel-1 ascending and from ALOS-2 and Sentinel-1 descending are 0.98 and 0.93, respectively.However, the formal errors of the GNSS velocities are higher than the InSAR's (see Table 6).This might be related to the number of observations.Meanwhile, in the InSAR techniques, there are numerous observations; seven campaigns were carried out for GNSS observations.Nevertheless, the GNSS velocities can be compared with the InSAR technique since, in almost all stations, the directions of the movements are the same for the GNSS and InSAR methods.Moreover, we calculated the correlation values of the displacements between the InSAR and GNSS techniques for all three datasets for 150 m, 200 m, and 250 m.During the correlation calculation, the displacements of the PS points within the distinct buffer zone radii and the displacements derived from the GNSS analyses, which were transformed into the LOS direction, were utilized.These analyses were performed across various buffer zones (Figure 4).For the 200 m radius circle, the highest correlation between the GNSS and InSAR techniques was found at GX12 for the ALOS-2 and Sentinel-1 datasets.One should note that since the incidence and heading angles are different for the InSAR datasets, the LOS values should be compared separately.Hence, we plo ed the time series for each dataset, as in Figure 5.It can be seen from the figure that since the LOS values derived from the GNSS perfectly fit its trend, the correlation values are also high.One should note that it is possible to compare the combination/integration of the InSAR One should note that since the incidence and heading angles are different for the InSAR datasets, the LOS values should be compared separately.Hence, we plotted the time series for each dataset, as in Figure 5.It can be seen from the figure that since the LOS values derived from the GNSS perfectly fit its trend, the correlation values are also high.One should note that it is possible to compare the combination/integration of the InSAR analysis with the GNSS.We also integrated the three InSAR datasets in this study.The results of the integration are presented in the next section.One should note that since the incidence and heading angles are di InSAR datasets, the LOS values should be compared separately.Hence, w time series for each dataset, as in Figure 5.It can be seen from the figure tha values derived from the GNSS perfectly fit its trend, the correlation values One should note that it is possible to compare the combination/integration analysis with the GNSS.We also integrated the three InSAR datasets in t results of the integration are presented in the next section.the next step, we tested the significance of the velocity difference GNSS and InSAR techniques to determine whether the differences are stat ent from zero.To carry this out, we calculated the ratios of the velocity dividing their standard deviations (Equation ( 6)).
where  is the test value,  and  are the velocities obtained from GN and  and  are the standard deviations of the velocities, respectively.were compared with the critical value with a 95% confidence level.
According to Table 7, the velocity differences for station GH19 are s other stations, all differences are insignificant.Although it can be seen from 6 that the velocity differences are up to 1-1.5 cm for GH06, the differences since the formal error of the GNSS is high compared to the InSAR.The read that the test values are the ratios of the velocity differences obtained by standard deviations.The test values are significantly lower because the G In the next step, we tested the significance of the velocity differences between the GNSS and InSAR techniques to determine whether the differences are statistically different from zero.To carry this out, we calculated the ratios of the velocity differences by dividing their standard deviations (Equation ( 6)).
where T is the test value, v G and v I are the velocities obtained from GNSS and InSAR, and s v G and s v I are the standard deviations of the velocities, respectively.The test values were compared with the critical value with a 95% confidence level.
According to Table 7, the velocity differences for station GH19 are significant.For other stations, all differences are insignificant.Although it can be seen from Tables 5 and 6 that the velocity differences are up to 1-1.5 cm for GH06, the differences are negligible since the formal error of the GNSS is high compared to the InSAR.The reader should note that the test values are the ratios of the velocity differences obtained by dividing their standard deviations.The test values are significantly lower because the GNSS velocities have more considerable standard deviations than InSARs.Nevertheless, in the stations (except GH06 and GH19), the trends of the GNSS and InSAR techniques show the same direction and fit.The maximum velocity difference is 3.68 mm/yr.In this subsection, we finally address the problem of decomposing the LOS-projected ground displacement rates into the region's 3D (up-down, east-west, and north-south) displacement rates.Generally speaking, at least three LOS-projected displacement rates recovered from the complementary (ascending/descending) orbits are needed to solve this problem.Our work combined the ALOS-2, Sentinel-1 ascending, and Sentinel-1 descending SAR datasets.Nonetheless, since the satellites fly near-polar orbits, the sensitivity of the LOS displacements to the north-south components of the deformation is significantly limited, and only the vertical and east-west displacement rates can accurately be retrieved.Indeed, it is known that the sensitivities of the SAR system to the 3D ground motions are as follows: where [U v , U E , U N ] are the 3D components of the ground deformation, θ is the local incidence angle, and α is the azimuth angle.Accordingly, the observation equation can be written as follows: where i, j, and k represent ALOS-2, Sentinel-1 ascending, and Sentinel-1 descending, respectively, [LOS i , LOS j , LOS k ] T is the vector of the relevant LOS-projected ground-displacement rates, and R is the transformation matrix: Instead of applying more sophisticated, time-consuming approaches for the calculation of the 3D ground displacement rates, relying on the combination of the LOS-projected ground displacement time series (e.g., see [28][29][30]), in this work, we followed the straightforward strategy outlined in [31], which assumes the deformations are almost linear over time (and this assumption is reliable in our case) and only requires LOS-projected ground displacement rates to be available for the three datasets over a common spatial grid.To this aim, the obtained ground displacement rates of the detected PS points were calculated over 50 m × 50 m grid points after masking the incoherent areas.Then, Equation (8) was applied to recover the 3D ground deformation rates.
Figure 6 shows the region's vertical (up-down) displacement rate map.As can be seen from the figure, the subsidence was obtained mainly in the coastal areas and the northern side of the GH.Moreover, Figure 7 shows a map of the region's horizontal (east-west) displacement rates.west) displacement rates.
In the final step, we calculated the vertical displacement rates of the GNSS stations based on the PS points within the 200 m radius circle (Table 8).According to Table 8, it can be concluded that the vertical movements of the stations show the same direction for almost all stations.The velocities obtained from the PPP solutions are consistent with those based on the PS points.The correlation of the velocities is 0.96 (except GH06 and GH19).For GH06 and GH19, we found relatively higher displacement rates with PPP.This might be related to the PS point number within the circle and the displacements of the PS points.Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Figure 1 .
Figure 1.Topography and geography of the studied area.Red, blue, and purple boxes highlight the swaths of the used Sentinel-1 (descending), Sentinel-1 (ascending), and ALOS-2 SAR datasets.A zoomed in view of the height profile of the zone identified by the black rectangle in (a) is shown in (b).

Figure 1 .
Figure 1.Topography and geography of the studied area.Red, blue, and purple boxes highlight the swaths of the used Sentinel-1 (descending), Sentinel-1 (ascending), and ALOS-2 SAR datasets.A zoomed in view of the height profile of the zone identified by the black rectangle in (a) is shown in (b).

Figure 2 .
Figure 2. GHGNSS network.Fourteen established stations in August 2017 are colored red, and the others are colored blue.The red star represents the study area.

Figure 2 .
Figure 2. GHGNSS network.Fourteen established stations in August 2017 are colored red, and the others are colored blue.The red star represents the study area.
are code and carrier phase observations, ρ j i is the geometric range, c is speed of light, ∆t i and ∆t j are the receiver and satellite clock errors, respectively.T j represents tropospheric delay.I j P,i,k and I j ϕ,i,k (= −I j P,i,k ) are the ionospheric delay, b P,i,k and b ϕ,i,k code and carrier phase biases for the receiver, b j P,k are b j ϕ,k code and carrier phase biases for satellite, λ k is carrier wavelength of frequency k, N j i,k is integer ambiguity, m j P,i,k and m j ϕ,i,k are code and carrier multipaths, ϵ j P,i,k and ϵ j ϕ,i,k are the measurement errors for code and carrier phase measurements, respectively.

Sentinel- 1 Figure 3 .
Figure 3. (a-c) Maps of the ground deformation velocities as obtained for ALOS-2, Sentinel-1 (ascending), and Sentinel-1 (descending) SAR datasets, respectively.The insets in (d-f) show the distribution of the detected PS points highlighted by the black rectangle in (a-c).

Figure 3 .
Figure 3. (a-c) Maps of the ground deformation velocities as obtained for ALOS-2, Sentinel-1 (ascending), and Sentinel-1 (descending) SAR datasets, respectively.The insets in (d-f) show the distribution of the detected PS points highlighted by the black rectangle in (a-c).

Figure 4 .
Figure 4.The graph of correlation values.Blue, green and yellow histograms refer to spatial boxes on the ground with a radius of 150 m, 200 m, and 250 m, respectively.

Figure 4 .
Figure 4.The graph of correlation values.Blue, green and yellow histograms refer to spatial boxes on the ground with a radius of 150 m, 200 m, and 250 m, respectively.

Figure 4 .
Figure 4.The graph of correlation values.Blue, green and yellow histograms refer on the ground with a radius of 150 m, 200 m, and 250 m, respectively.

Figure 5 .
Figure 5.The time series of GX12.

Figure 5 .
Figure 5.The time series of GX12.

Figure 6 .
Figure 6.Map of the vertical displacement rate computed over a 50 m × 50 m grid after masking the incoherent areas, interpolating PSI-derived LOS ground deformation velocity values and applying the adopted, simple combination method in[31].

Figure 6 . 17 Figure 7 .
Figure 6.Map of the vertical displacement rate computed over a 50 m × 50 m grid after masking the incoherent areas, interpolating PSI-derived LOS ground deformation velocity values and applying the adopted, simple combination method in [31].Sensors 2024, 24, x FOR PEER REVIEW 15 of 17

Figure 7 .Funding:
Figure 7. Same as Figure 6 for the east-west ground displacement rates of the region.

Table 1 .
The SAR data used in the study.

Table 2 .
The details of the campaigns.

Table 3 .
Number of PS points for the ALOS-2 dataset.

Table 4 .
Number of PS points for Sentinel-1 datasets.

Table 5 .
Estimated velocities and formal errors of the PS points within the 200 m buffer zone of the GNSS stations using different InSAR datasets obtained from InSAR analysis (mm/yr).
i : velocity, s i : standard deviation of velocity, i: dataset.

Table 6 .
Estimated velocities and formal errors of the GNSS stations, transformed into LOS direction for different orbit trajectories, obtained from GNSS analysis (mm/yr).
i : velocity, s i : standard deviation of velocity, i: dataset.

Table 7 .
The test and critical values of the significance test.The bold values indicate significant velocities.

Table 8 .
Vertical displacement rates of the stations (mm/yr).