Time series analysis of L-band PALSAR-2 images in Istanbul and Kocaeli, Turkey

ABSTRACT

Conducting long measurements of infrastructure deformation is a critical engineering task. Conventional methods are both time-consuming and expensive, limiting their use for large-scale applications. The synergy of synthetic aperture radar (SAR) and geographic information systems (GIS) offers a complementary approach. This study focuses on the feasibility of using time series analysis of L-band PALSAR-2 images to discover land displacements in Istanbul and Kocaeli, significant industrial and residential areas in Turkey. PALSAR-2 phase and intensity information were analyzed. For phase analysis, 14 L-band images from 2014 to 2021 were taken into account. Small baseline subset (SBAS) analysis was performed using 44 pairs, and results of the velocity, coherence and backscattering values are presented. Coherence of all pairs and their correlations were calculated. Principal Component Analysis (PCA) reduced the dimension of coherence pairs, enhancing feature extraction and the final geocoded velocity map revealed a fastest subsidence rate of −58 mm/yr and a mean subsidence of −20 mm/yr. These findings were confirmed through mean vertical velocity from Sentinel-1 datasets and field observations. The results showed that immature land subsidence in the mentioned areas are growing slowly, which can be taken as a serious risk in future.

KEYWORDS:

• InSAR time series
• PALSAR-2
• SAR intensity
• infrastructure deformation

1. Introduction

Conducting long and effective measurements of the displacement of infrastructures are among required critical engineering tasks. The expansion of industrial and tourist activities has resulted in higher numbers of factories, roads and railway networks (Chang et al., 2017; Hu et al., 2019; Karimzadeh & Matsuoka, 2021; Qin et al., 2017; Wang et al., 2021). The conventional methods for discovering and monitoring structural health basically use in situ measurements by experts. Conventional methods, such as GPS (Lovse et al., 1994; Lepadatu & Tiberius, 2014) and precise leveling (Chen et al., 2012), are also time-consuming, and costly for the employers who cannot apply them at any scale. These limitations make the conventional methods difficult to utilize for real applications in a rather large scale. The synergy of synthetic aperture radar (SAR) remote sensing and geographic information systems (GIS) is not only a complementary method for conventional techniques, but it also can be used as a useful technique for discovering hidden crises in earlier stages, so that the necessary steps can be taken to control or manage such crises (Garg et al., 2021; Karimzadeh & Matsuoka, 2018).

Land subsidence is a gentle downward movement of the Earth’s surface, which is identified as one of the problematic environmental issues in recent decades (Bull & Polland, 1975). This is not only causing serious damage to agricultural lands, but it can also lead to irreparable results in settlement and industrial areas. Typically, the early stages of land subsidence cannot be scaled by the human eyes. Thus, special attention to such long-term phenomena requires special tools other than normal human visual capabilities. PALSAR-2 SAR images have a potential to penetrate tall vegetation and collect dataset from the bare land. Along with high penetration capability of the L-band sensors, high-resolution imagery mode of ALOS-2 satellite allows us to study urban areas together with other C-band and X-band sensors. PALSAR-2 SAR remote-sensing datasets also have shown an acceptable ground coverage (swath) in the recent years, ranging from between 15 km × 15 km (range × azimuth) to 490 km × 490 km (range × azimuth). However, the revisit interval of L-band satellites such as ALOS-2 is higher than the recent C-band and X-band missions such as Sentinel-1 and TerraSAR-X. The main question here is whether we can develop a semi-automated approach to detect land subsidence or changes in urban areas at early stages to mitigate their consequences. As the revisit period of the SAR satellite has been decreasing in recent years (Cigna et al., 2021; Karimzadeh & Matsuoka, 2018, 2020b), the expectation for more effective, semi-kinematic displacement analysis on a centimeter or millimeter scale has increased (Ezquerro, Del Soldato, et al., 2020; Ezquerro, Tomás, et al., 2020), not only in regional scales but also in nationwide scales (Lazecký et al., 2020; Morishita et al., 2020). For example, LiCSAR (Looking Into Continents from Space with Synthetic Aperture Radar) is an operational system which is designed for SAR interferometric (InSAR) analysis of Sentinel-1 data in a large scale (Morishita et al., 2020). Such continuous spatial and temporal coverage of displacement analysis is an invaluable source of information for monitoring of various short-term (e.g. earthquake) and long-term (e.g. land subsidence) phenomena on the Earth. Aiming for effective infrastructure health monitoring systems, automated and semi-automated data analysis techniques have also been presented recently with very-high-resolution images such as WorldView-2 and object-based image analysis (OBIA) (Miura et al., 2020; Omarzadeh et al., 2021; Xiu et al., 2022). Chen et al. (2019) employed a combination of Persistent Scatterer Interferometric Synthetic Aperture Radar (PS-InSAR) and tomographic SAR techniques. Their research utilized TerraSAR-X and CosmoSkyMed imagery to delve into the displacement analysis of structures including skyscrapers, bridges, and historical landmarks situated in Nanjing City, China. Mele et al. (2022) introduced a novel approach involving joint Differential Interferometric Synthetic Aperture Radar (DInSAR) and Density-Based Spatial Clustering of Applications with Noise (DBSCAN) techniques. This innovative methodology was designed to identify clusters of pixel values associated with buildings. In a subsequent study, Mele et al. (2023) studied health of structures in Barcelona, Spain, using a basic concept for active displacement areas (ADA) using a semi-automated ADAfinder tool for the identified PSs in the study area.

The necessity of this study is to establish a new geodatabase incorporating diverse spatial layers, consolidating them to provide comprehensive insights into infrastructure changes. Here, we introduce a semi-automated GIS-based analysis of spatiotemporal ground characteristics using InSAR time series analysis, focusing on areas within Istanbul and Kocaeli Provinces. The method utilizes coherence and intensity information from PALSAR-2 (StripMap 3) images (with a wavelength of approximately 24 cm) from JAXA (the Japan Aerospace Exploration Agency). This approach is suitable for establishing temporal or spatial relationships among processed intensity, principal components of coherence values, and InSAR displacement maps, providing valuable insights into changes in the Earth’s surface.

2. Study area and datasets

2.1. Study area description

The present study mainly focused in one frame of the ALOS-2 PALSAR-2 images, which cover major parts of Istanbul and Kocaeli Provinces. Both provinces have experienced rapid development in terms of urbanization and industrial activities. About 16 million people are settled in Istanbul, and this number is expected to increase in the coming years. Furthermore, about 2 million people inhabit Kocaeli (Güven, 2022; Halicioglu et al., 2021). The heavy industrial projects west of Kocaeli and all over Istanbul meet the demands of not only domestic people, but also of international markets. According to Halicioglu et al. (2021), such parallel growth of population and industrial projects manifests itself in the future with serious difficulties and consequences, as both provinces are located in an earthquake-prone area (Aslan et al., 2018; Güven, 2022; Özarpacı et al., 2022). Thus, intensive monitoring of the infrastructures for tectonic and non-tectonic hazards must gain more value for sustainable urban developments in the future. According to Akarvardar et al. (2009), Avcilar district in Istanbul experienced ground subsidence and the same area damaged by the M 7.4 Izmit (Kocaeli) earthquake of 1999 which was mainly associated with compressible material in shallow subsurface layers. Here, we only selected two candidate areas in Istanbul and Kocaeli Provinces, which shows land subsidence displacement as follows: 1—Kadikoy, which is a filled land near a water canal in Istanbul; and 2—Balcik, or Balçık in Kocaeli. Moreover, spatiotemporal changes in the newly inaugurated third bridge of Istanbul (i.e. Yavuz Sultan Selim Bridge), together with the connection roads in the Asian and European sides, are taken into account to investigate the length of the roads constructed within the period of about 7 years, from November 2014 to March 2021.

2.2. Datasets description

Using ERS-1/2, Envisat and Sentinel-1, all in C-band domain (λ ~ 6 cm), and a persistent scatterer interferometric (PS-InSAR) analysis, the ground displacement rate in Istanbul is well studied, especially from the point of view of geological descriptions in previous studies from 1992 to 2017 (Aslan et al., 2018; Stramondo et al., 2010, 2011). From 2015 to 2018, using the same method and Sentinel-1 images, other research was carried out for Istanbul. Alternatively, we used PALSAR-2 images of the ALOS-2 satellite with its better spatial accuracy (~10 m) to provide more information about the (in)stability of the study area. Since X-band SAR datasets are not available for the study area, the finer spatial resolution of the L-band images gain importance for monitoring of urban areas. The PALSAR-2 aboard the ALOS-2 satellite is a SAR sensor, which transmits and receives the backscattered microwave energy from the ground features to provide information about the stability, soil moisture, surface roughness, surface changes including structural damage and fire. This sensor is independent of the source of light such as the sun, and accordingly, PALSAR-2 benefits from the advantage of providing satellite images, both day and night. The wavelength of L-band (~24 cm) used in this sensor is less affected by cloudy, rainy or snowy weather, which provides a suitable capability for rapidly monitoring various objects, structures and phenomena (Aimaiti et al., 2018; Natsuaki, 2012; Watanabe et al., 2016; ALOS-2 Product FormatDescription, 2014). Additionally, L-band wavelength penetration is fair enough to reach the ground surface through vegetation, in order to obtain information about the ground surface. PALSAR-2 has SpotLight imagery mode, with a resolution of 1 m × 3 m in the azimuth and range directions, respectively. It allows for the collection of detailed information for small-scale analysis, such as urban analysis. Additionally, StripMap and ScanSAR modes are available for large-scale projects. The spatial resolution for StripMap1 is 3 m (Ultrafine); for StripMap 2, it is 6 m (High sensitive); for StripMap 3, it is 10 m (Fine) (Karimzadeh & Matsuoka, 2020a; Natsuaki, 2012; ALOS-2 Product Format Description, 2014). In this study, only StripMap 3 was available, and accordingly, for Istanbul and Kocaeli, 14 L-band PALSAR-2 images from ascending orbits between 2014 and 2021 were gathered. Table 1 shows the detailed information of the used images, and Figure 1 also shows the location of the study area and the used SAR track. The total area analyzed by this study was roughly 7300 km². The path and row numbers of the acquired PALSAR-2 scenes were 188 and 810, respectively. The central latitude and longitude of the scenes were 29.382° and 41.083°, respectively.

Figure 1. The study area with the locations of the investigated sites in Istanbul and Kocaeli provinces. The green box shows the track 188 PALSAR-2 (ascending orbit). Red circles show the investigated sites kadikoy and Balcik. The location of the Yavuz Sultan Selim Bridge is also shown by the red line. The inset is the map of Turkey.

Table 1. SAR images and detailed characteristics of the L-band PALSAR-2 dataset used in this study. “*” indicates super master image, “A” indicates ascending orbit, respectively.

Label (#)	Date (YYYY/MM/DD)	Polarization	Incidence Angle (°)	Band	Orbit	Path Number
0	2014/11/16	HH/HV	40.57	L-band	A	188
1	2015/01/25	HH/HV	40.57	L-band	A	188
2	2015//10/08	HH/HV	40.57	L-band	A	188
3	2015/11/15	HH/HV	40.57	L-band	A	188
4	2016/06/26*	HH/HV	40.57	L-band	A	188
5	2016/10/30	HH/HV	40.57	L-band	A	188
6	2017/09/03	HH/HV	40.57	L-band	A	188
7	2017/10/29	HH/HV	40.57	L-band	A	188
8	2018/01/07	HH/HV	40.57	L-band	A	188
9	2018/04/29	HH/HV	40.57	L-band	A	188
10	2018/12/23	HH/HV	40.57	L-band	A	188
11	2019/03/17	HH/HV	40.57	L-band	A	188
12	2020/03/15	HH/HV	40.57	L-band	A	188
13	2021/03/14	HH/HV	40.57	L-band	A	188

Although the main objective of this study was related with time-series interferometric SBAS analysis, intensity and GIS-based analyses were included as well. Since the main targets in the study area were infrastructures (roads, buildings and railways), only HH polarization of the dual polarized PALSAR-2 data (HH+HV) were used. The choice between HH and HV polarizations in PALSAR-2 data for interferometric analysis depends on the specific objectives of the study. In this study, HH polarization is preferred over HV for interferometric analysis due to higher coherence results (Vaglio Laurin et al., 2022). The incidence angle for all images was 40.57 degrees, and the difference of incidence angle during the SBAS analysis (Δθ) was not more than 0.02°, which implied that the acquisition positions or normal baselines were quite acceptable.

2.2.1.

3. Methodology

Although image classification techniques for optical images can provide useful information for urban areas (Ghasemi et al., 2021), monitoring changes and displacement is more sensitive in SAR images since they deal with phase and intensity information (Karimzadeh & Matsuoka, 2020b; Karimzadeh, Ghasemi, et al., 2022). However, the interpretation of SAR analysis is sometimes complex and difficult, due to the complex behavior of backscattering coefficients and the unique characteristics of side-looking SAR systems. Thus, this study used both phase and intensity information analyzed by Sarscape 5.5 for all possibilities in a QGIS database, in order to achieve a reliable analysis of urban elements. The methods and techniques used in this study can be divided into three main branches, as follows: 1—intensity analysis; 2—interferometric analysis; and 3—GIS analysis (Figure 2).

Figure 2. Flowchart of the proposed research.

3.1. Intensity analysis

In this study, the PALSAR-2 images were gathered in SLC (single look complex) format in level 1.1. In order to visualize the high-resolution L1.1 images, a pre-processing step was necessary to bring both phase and intensity information to an actionable manner. The L1.1 is a kind of raw data, right after the imagery, and it is the fundamental product for the main processes (ALOS-2 Product Format Description, 2014). The geometry of SAR imaging, rather than a geographic concept, is based on slant range and azimuth coordinates. Thus, in a side-looking SAR sensor, the distortions mainly based on topographic changes are considerable, and an image calibration process and coordinate conversion (from SAR coordinate to geographic coordinate) should be conducted. This conversion is possible using an inverse technique and a digital elevation model (DEM). The Doppler frequencies, the time for pulse transmission time of the SAR image, and the position and velocity vectors of both the sensor and the backscatter feature are used for calibration, as follows (Meier et al., 1993):

(1) $\mathit{R}=\mathit{S}-\mathit{P}$(1)

(2) $f_D=\mathit{\hspace{0.17em}}\frac{2f_0\left(v_p-v_s\right)R_s}{\mathit{c}\left|R_s\right|}$(2)

where $\mathit{R}$ is the slant range, $\mathit{S}$ is the satellite position, $\mathit{P}$ is the backscatterer position, $v_s$ is the velocity of the satellite, $v_p$ is the velocity of the backscatterer feature, c is the speed of light, $f_0$ is the frequency of the carrier, and $f_D$ is the Doppler frequency. During geocoding processing and after processing the Doppler frequencies, the output can be transformed to any reference system. To convert digital numbers (e.g. imaginary and real values) to sigma naught (${\sigma }^0$) or backscattering coefficients, JAXA provides the following equation for PALSAR-2 ALOS-2 products (ALOS-2 Product Format Description, 2014; Karimzadeh et al., 2022):

(3) ${\sigma }^0=10.\mathit{lo}{g}_{10}⟨\mathit{D}{N}^2⟩+\mathit{C}{F}_1$(3)

where ${\sigma }^0$ is the backscattering coefficient in the dB unit, $\mathit{DN}$ is a digital number (or raw pixel value), and $\mathit{C}{F}_1$ is the calibration constant adjusted for the PALSAR-2 ALOS-2 (−83 dB).

To suppress speckle noise on the geocoded intensity images, we used the Lee filter with a window size of 3 × 3 to the geocoded image, in order to reduce the level of noise, especially in the urban areas (Lee, 1980; Miura et al., 2016). We also applied a multilooking operation to reduce the speckle noise as much as possible by setting a 1 × 3 look size in the range and azimuth directions, respectively. It must be noted that applying large look sizes may reduce the spatial resolution dramatically, which is not ideal for the purpose of this study. Different look sizes can also change the dimensions of the pixels negatively if the pixels are “rectangle” instead of “square”. Scatter plots of the intensity pixel values, with respect to the super master image (2016.06.26) for Kadikoy and Balcik regions, are shown in Figures 3 and 4, respectively. A comparison of the intensity values shows that the correlation coefficient (CC) is high when the time interval of each image with the super master image is short for both Kadikoy and Balcik.

Figure 3. Scatter plots and correlation coefficients (CCs) of the intensity pixel values of the PALSAR-2 images with respect to the super master image (2016.06.26) in Kadikoy.

Figure 4. Scatter plots and correlation coefficients (CCs) of the intensity pixel values of the PALSAR-2 images with respect to the super master image (2016.06.26) in Balcik.

3.2. SBAS interferometric analysis

Several InSAR studies using ERS-1/2, Envisat, Sentinel-1 datasets and persistent scatterer interferometry (PS-InSAR) were carried out in the study area (Aslan et al., 2018; Halicioglu et al., 2021; Stramondo et al., 2010, 2011). The performance of permanent scatterers (PS) interferometry technique is superior to SBAS-InSAR when the number of SAR images is relatively high, but the SBAS-InSAR is less sensitive to temporal decorrelations compared to PS-InSAR (Aimaiti et al., 2018; Karimzadeh, Matsuoka, & Ogushi, 2018).

Range and azimuth looks are considered as the same of the intensity analysis (i.e. 1 range look and 3 azimuth look). The selection of candidate interferometric pairs was based on their short temporal and spatial baselines. In order to form an SBAS network, we assumed the PALSAR-2 images to be an array of SLC images ($\mathit{N}$) for the same area, in which all of the images were taken by the same geometry, orbit and path at the ordered acquisition dates (t1, t2, t3, … , tn). The outcome of the interferometric correlation analysis between the candidate SBAS pairs was directly related with the spatial (normal) baseline of the two images ($\mathit{B}$), Doppler centroid of the two images ($\mathit{D}$), temporal gap of the two images ($\mathit{T}$), and miscellaneous sources of the noise generated by the sensor. Assuming the noise is negligible, the general SAR interferometic correlation is as follows (Berardino et al., 2002; Karimzadeh & Matsuoka, 2020a):

(4) $\mathit{\rho }={\rho }_{\mathit{temporal}}\cdot {\rho }_{\mathit{spatial}}\cdot {\rho }_{\mathit{Doppler}}\approx (1-\mathit{f}\left(T_{\mathrm{\perp }}-T_{\perp }^c\right)).\left(1-\mathit{f}\left(B_{\mathrm{\perp }}-B_{\perp }^c\right)\right)\cdot \left(1-\mathit{f}\left(D_{\mathrm{\perp }}-D_{\perp }^c\right)\right)$(4)

where $\mathit{\rho }$ is the interferometric joint correlation; $T_{\mathrm{\perp }}$, $B_{\mathrm{\perp }}$ and $D_{\mathrm{\perp }}$ are time interval, spatial interferometric baseline and Doppler difference, respectively; and $T_{\perp }^c$, $B_{\perp }^c$ and $D_{\perp }^c$ are critical time interval, critical normal baseline and Doppler centroid difference, respectively. The interferometric correlation ($\mathit{\rho }$) is approximated by considering the product of individual correlations, where each component is influenced by differences from reference values through the function ($\mathit{f}$). The criteria defined by expert or user were taken into account, based on the Equation (4)(4) $\mathit{\rho }={\rho }_{\mathit{temporal}}\cdot {\rho }_{\mathit{spatial}}\cdot {\rho }_{\mathit{Doppler}}\approx (1-\mathit{f}\left(T_{\mathrm{\perp }}-T_{\perp }^c\right)).\left(1-\mathit{f}\left(B_{\mathrm{\perp }}-B_{\perp }^c\right)\right)\cdot \left(1-\mathit{f}\left(D_{\mathrm{\perp }}-D_{\perp }^c\right)\right)$(4) . If the criteria matched with the SBAS requirements, the candidate SAR pair could be included into the SBAS analysis. Here, the maximum time interval was 730 days, and the maximum spatial baseline was 45% of the critical baseline. One image was selected as a “super master” image. The selection of the super master image was based on its fair spatial and temporal baselines with all the other images. The image with label 4 (2016.06.26) in Figure 5 was ideal to be a super master image. L-band PALSAR-2 images are in the long wavelength range of images that pose a higher phase correlation compared to C-band and X-band products. In this study, in spite of large temporal gaps between the images in the L-band PALSAR-2 images, the phase correlation was still in good condition, and we assumed that the acquired images could be contributed into the SBAS analysis. As shown in Figure 5(a,b), forty-four pairs formed a solid SBAS network to reduce the possible errors. The number of PALSAR-2 images was not too high in this study, but was still enough to form the SBAS analysis, as discussed in previous studies (Karimzadeh & Matsuoka, 2020a, 2020b). Due to the relatively long time interval of the interferometric pairs, finding the optimal coherence threshold required multiple attempts. After careful consideration, we selected the value 0.7 as the optimal threshold to achieve the desired balance and consequently all the low values were removed from the analysis. The minimum and maximum normal baselines were 25 m and 720 m, respectively.

Figure 5. (a) SBAS SAR network of the analysis for the PALSAR-2 images (black dots). The super master image is highlighted by a yellow dot; (b) arrangement of InSAR pairs and the corresponding pair number over each link.

In this study, SBAS analysis proposed by Berardino et al. (2002) was carried out for entire processes. The minimum number of pairs was one pair with two images, and the phase correlation for N as the odd number of images and M as the number of interferograms could be calculated with Equation (5)(5) $\frac{\mathit{N}+1}{2}\le \mathit{M}\le \mathit{N}\left(\frac{\mathit{N}+1}{2}\right)$(5) :

(5) $\frac{\mathit{N}+1}{2}\le \mathit{M}\le \mathit{N}\left(\frac{\mathit{N}+1}{2}\right)$(5)

Pixel values for the unwrapped interferograms’ ($\mathit{j}$) radar coordinate system (in azimuth and range), between the two acquisitions ($t_B\mathit{\hspace{0.17em}}$ and $t_A$), could be determined as follows (Karimzadeh & Matsuoka, 2020a):

(6) ${\phi }_j\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)=\mathit{\phi }\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)-\mathit{\phi }\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$(6)

Equation (6)(6) ${\phi }_j\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)=\mathit{\phi }\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)-\mathit{\phi }\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$(6) can be rewritten as follows:

(7) ${\phi }_j\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)=\frac{4\mathit{\pi }}{\lambda }\left[d^{\mathit{LP}}\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)-d^{\mathit{LP}}\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)\right]+\mathrm{\Delta }{\phi }_j^{atm}\left(t_B,\mathit{\hspace{0.17em}}{t}_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)+\mathrm{\Delta }{\phi }_j^{\mathit{topo}}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$(7)

where $\mathit{j}$ is an integer from 1 to $\mathit{M}$; $\mathit{\phi }\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$ and $\mathit{\phi }\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$ are phase values for the multilooked SAR images in two acquisitions ($t_B\mathit{\hspace{0.17em}}$ and $t_A$); $d^{\mathit{LP}}\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$ and $d^{\mathit{LP}}\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$ are displacement values of low-pass components between $t_A\mathit{\hspace{0.17em}}$ to $t_B$ in the line-of-sight (LOS) of the satellite; $\mathrm{\Delta }{\phi }_j^{atm}\left(t_B,\mathit{\hspace{0.17em}}{t}_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$ is the related atmospheric between $t_A\mathit{\hspace{0.17em}}$ and $t_B$; $\mathit{\lambda }$ is the used wavelength associated with satellites (e.g. for ALOS-2, the wavelength was 23.6 cm); $\mathrm{\Delta }{\phi }_j^{\mathit{topo}}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$ is the phase of topography for the Earth, which is not flat in reality. The phase of topography was defined as follows:

(8) $\mathrm{\Delta }{\phi }_j^{\mathit{topo}}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)\approx \frac{4\mathit{\pi }}{\lambda }\frac{{\mathit{B}}_{\mathrm{\perp j}}\mathrm{\Delta z}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)}{\mathit{rsin\theta }}$(8)

Here, $B_{\mathrm{\perp j}}$ is the normal baseline of the two SAR acquisitions, as shown in Figure 5, and $\mathit{\theta }$ is the incidence angle of the SAR images. As mentioned in Table 1, the incidence angle for the used images was 40.57°, and $\mathrm{\Delta z}$ was the exposed topographic effect.

In Equation (7)(7) ${\phi }_j\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)=\frac{4\mathit{\pi }}{\lambda }\left[d^{\mathit{LP}}\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)-d^{\mathit{LP}}\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)\right]+\mathrm{\Delta }{\phi }_j^{atm}\left(t_B,\mathit{\hspace{0.17em}}{t}_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)+\mathrm{\Delta }{\phi }_j^{\mathit{topo}}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$(7) , the noise effects (which are primarily associated with instrumental performance) are assumed to be negligible. Refinement and re-flattening of the interferograms to remove the turbulent phases were performed by selecting approximately 100 ground control points in motionless and highly coherent areas and using a 3 arc-second Shuttle Radar Topography Mission (SRTM) DEM and a third-order polynomial refinement equation. Equations (7)(7) ${\phi }_j\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)=\frac{4\mathit{\pi }}{\lambda }\left[d^{\mathit{LP}}\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)-d^{\mathit{LP}}\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)\right]+\mathrm{\Delta }{\phi }_j^{atm}\left(t_B,\mathit{\hspace{0.17em}}{t}_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)+\mathrm{\Delta }{\phi }_j^{\mathit{topo}}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$(7) and (8(8) $\mathrm{\Delta }{\phi }_j^{\mathit{topo}}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)\approx \frac{4\mathit{\pi }}{\lambda }\frac{{\mathit{B}}_{\mathrm{\perp j}}\mathrm{\Delta z}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)}{\mathit{rsin\theta }}$(8) ) describe a system with infinite solutions because of the N + 1 unknowns. The singular value decomposition (SVD) method, which does not require prior information about the temporal behaviour of the deformed area, provides a simple solution to invert the system. In the final step of refining the surface displacement data ($d^{\mathit{LP}}\left(t_i,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$, i = 0, … , N), the phase related to the atmospheric effects ($\mathrm{\Delta }{\phi }_j^{atm}\left(t_B,\mathit{\hspace{0.17em}}{t}_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$) must be detected and mitigated by considering that the possible atmospheric phase components are clustered primarily in space and time (Berardino et al., 2002). Pixels with coherence values smaller than our defined threshold (0.7) were set to dummy variables (NaN).

Once the displacement map in the LOS was estimated with Equation (7)(7) ${\phi }_j\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)=\frac{4\mathit{\pi }}{\lambda }\left[d^{\mathit{LP}}\left(t_B,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)-d^{\mathit{LP}}\left(t_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)\right]+\mathrm{\Delta }{\phi }_j^{atm}\left(t_B,\mathit{\hspace{0.17em}}{t}_A,\mathit{\hspace{0.17em}}\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)+\mathrm{\Delta }{\phi }_j^{\mathit{topo}}\left(\mathit{x},\mathit{\hspace{0.17em}}\mathit{r}\right)$(7) , the LOS displacement could be converted to vertical displacement. If the descending and ascending datasets were available simultaneously, 3 displacement components could be retrieved (E – W, N – S, and Vertical). Due to a lack of the descending datasets, only vertical values could be estimated by considering some assumptions: 1—the LOS results contained a few N – S movements, since the SAR system was not sensitive for the N – S direction; thus, we assumed N – S movement did not exist; and 2—the E – W movement was trivial, since field observation proved vertical displacement displacement due to long subsidence; thus, only the vertical displacement could be inferred as follows (Karimzadeh & Matsuoka, 2018; Su et al., 2017):

(9) $V_v=\frac{{\mathit{V}}_{\mathit{LOS}}+\mathit{\hspace{0.17em}}{V}_{\mathit{EW}}\mathit{sin\theta cos\alpha }}{\mathit{cos\theta }}$(9)

Here, $V_v$ is the vertical velocity, $V_{\mathit{LOS}}$ is the LOS velocity, $V_{\mathit{EW}}$ is the east‒west velocity, $\mathit{\alpha }$ is the azimuth angle, and $\mathit{\theta }$ is the incidence angle of the satellite.

Coherence of all pairs, their correlation and principal component analysis (PCA) were calculated after the velocity conversion. Coherence, as an index for the quality of interferometric pairs and a measure for explaining the Earth’s changes, is defined as follows (Karimzadeh & Matsuoka, 2018):

(10) $\mathit{\gamma }=\frac{\mathit{E}⟨\mathit{ab}\ast ⟩}{\sqrt{\mathit{E}⟨\mathit{aa}\ast ⟩E⟨\mathit{bb}\ast ⟩}}$(10)

where $\mathit{E}$ is the expected value, $\mathit{a}$ and$\mathit{\hspace{0.17em}}\mathit{b}$ are master and slave images, $\ast$ is the complex conjugate of the used images. As shown in Figure 6, the histogram of coherence values in different epochs shows different peaks for Kadikoy and Balcik, ranging between 0 and 1. In both regions, a large temporal gap manifested itself in a normal distribution histogram, in which the peak formed in lower coherence values. In contrast, the short temporal gap manifested itself in a normal distribution histogram, in which the peak formed in higher coherence values. In Figure 6(a,b), the highest peak belonged to the SAR pair with 2015.10.18 and 2015.11.15, for both the Kadikoy and Balcik regions. The peak coherences for Kadikoy and Balcik were 0.83 and 0.82, respectively. Furthermore, the average coherence values in the shortest pairs for Kadikoy and Balcik were 0.74 and 0.70, respectively.

Figure 6. (a) Histogram of coherence values for different pairs in Kadikoy; (b) histogram of coherence values for different pairs in Balcik.

3.3. GIS and PCA analyses

In order to turn the results of the SAR analysis into an actionable status, geodatabase creation and GIS analysis were necessary. Geodatabase is a collection of layers with a geographical concept in a common management system such as Microsoft SQL (Mohamadzadeh et al., 2020). Geodatabase contains not only the results of intensity and interferometric analyses, but also geospatial layers and information, including building inventories, roads and railways. Geodatabase is a comprehensive information model that keeps the attribute tables, and a series of raster values (e.g. velocity values) in the integrated folder accessible for eligible users. Thus, generation of such a database increases the value of the remote-sensing data by providing easy access for the required information and statistics of GIS analyses. We used the road, railway and building information of the OpenStreetMap (OSM) in the GIS Geodatabase. Since the OSM layers sometimes have geometric inconsistencies, the calculation of area (for parcels or buildings) and length (for railways and lines) may produce negative values (Haklay & Weber, 2008; Ram, 2022; Sehra et al., 2016) because outer lines in OSM shapefile (*.shp) that are snapped in a clockwise direction and inner rings that are snapped in a counter-clockwise direction generate a polygon with some errors. This is also called “incorrect ring order”, which creates negative values in the shapefile. In order to fix the errors, the geometric correction was carried out in the GIS environment. Once the area and length for the corresponding layers were calculated, the displacement thresholding for Kadikoy and Balcik was applied in a 5 km × 5 km area. The error sources in OSMs can also be reduced if a unique georeferenced system such as WGS-84 (World Geodetic System) is exploited together with the positioning of the candidate SBAS pixels (Miano et al., 2022). Using the spatial queries, the mean velocity values of the SBAS analysis were extracted to the centroid of the buildings, while for roads and railways, the 100 m buffer was used to extract the displacement values.

To further investigate, we applied PCA for the generated coherence pairs. PCA is a quick and a reliable visualization technique for image interpretation, especially when the dimensions of the dataset are rather high (Karimzadeh & Matsuoka, 2018; Kim et al., 2021). It offers several contributions for:

1. Dimensionality reduction: coherence pairs involve multi-temporal pixel values, leading to high-dimensional datasets. PCA helps in reducing the dimensionality by transforming the original data into a new set of uncorrelated variables called principal components. This simplifies the dataset while retaining the most critical information, making it easier to visualize and analyze. Here, the forty-four coherence maps are a large source of information, and sometimes redundant information makes interpretation difficult. The coherence maps (from SBAS analysis) were assumed into a matrix of scores of different locations, and the variables were different time spans.

2. Visualization: PCA facilitates the visualization of SAR datasets. The reduced-dimensional principal components can be plotted, providing a visual representation of the primary variation in the data. This visualization aids in better understanding the relationships and structures within the dataset.

3. Feature Extraction: coherence results might contain subtle variations (e.g. non-urban changes) that are challenging to detect directly. PCA can uncover hidden features by highlighting variations that might otherwise go unnoticed. This can be particularly useful for uncovering localized displacements or changes in different parts of the dataset.

The PCA calculation was defined as follows (Karimzadeh & Matsuoka, 2018):

(11) $\mathit{X}=\mathit{TP}{ }^{\mathrm{\prime }}+\mathit{r}$(11)

where $\mathit{X}$ is the PCA matrix of different time spans; $\mathit{T}$ is a matrix of the scores; $\mathit{P}$ is a matrix containing the loading factors (transposed matrix); and $\mathit{r}$ is the residuals of the results. The matrix of principal component scores ($\mathit{T}$) represents the transformed data in the new coordinate system defined by the principal components. Each row in $\mathit{T}$ corresponds to an observation (i.e. coherence map), and each column corresponds to a principal component. It must be noted that PCA originates from a projection method. The size of the dataset was ideal for PCA analysis, and we only created six principal components (PCs), since the first six components can generally explain the whole dataset.

4. Results

Using PALSAR-2 SAR images and the 30 m SRTM (shuttle radar topography mission) DEM of the study area, we performed SAR interferometry and PCA analyses.

4.1. SBAS displacement and field observations

Figure 7 shows the results of the SBAS InSAR for Istanbul and Kocaeli, from November 2014 to March 2021, using HH PALSAR-2 data. Although positive values may be indicative of natural processes such as tectonic uplifts, the predominant cause of positive values in this study is related with tall vegetation and different noises. Thus, the noisy and positive values were removed from the final velocity map. Warm colors are negative values, which indicate downward displacement (mostly land subsidence or landslide movement), while the green color represents motionless regions. The fastest subsidence rate reached −58 mm/yr in landslide and hill slope areas (dashed polygons in Figure 7). These areas are located in a military training area and forests. However, gentle subsidence in Kadikoy and Balcik was confirmed by a field visit in November 2022. Figure 8 shows buildings, roads and railways in Kadikoy overlaid on the InSAR vertical velocity map. The maximum vertical velocity rate in Kadikoy was −30 mm/yr in the central parts of the oval-shaped region, while the rates of subsidence in the marginal parts were less than −30 mm/yr. The number of buildings in the Kadikoy (5 km × 5 km) region was 9691 buildings, in which 1754 buildings (parcel) were located in a gap area without InSAR displacement values. The threshold for the critical subsidence rate was calculated from the sample displacement data of the geodatabase, the number of whole buildings and GIS query analysis. The number of buildings in the critical area (red area in Figure 8(a)) with a subsidence rate more than −5 mm/yr was found to be 594. The total area of the critical buildings was 121,875 m² in Kadikoy. More details about the length of the exposed roads and railways in both Kadikoy and Balcik are given in Table 2.

Figure 7. Mean vertical velocity map of Istanbul and Kocaeli. Black dashed polygons are mainly landslide areas or forest. The black rectangles are Kadikoy and Balcik.

Figure 8. (a) Kadikoy buildings overlaid on the mean vertical velocity (MVV) map; (b) roads overlaid on the MVV map; (c) railways overlaid on the MVV map. The locations of field observations and the water canal (blue polygon) are shown by the numbers in (c).

Table 2. Detailed information of affected length of railways and roads in PALSAR-2 and Sentinel-1 displacement maps.

PALSAR-2 displacement map	Total Road Length (km)	Total Railway Length (km)	Mean Road Subsidence Rate (mm/year)	Maximum Road Subsidence Rate (mm/year)	Mean Railway Subsidence Rate (mm/year)	Maximum Railway Subsidence Rate (mm/year)
Kadikoy	268	144	−8	−28.8	−5	−28
Balcik	330	-	−5	−14.7	-	-

Figure 9 shows buildings and roads in Balcik overlaid on the InSAR vertical velocity map. The fastest vertical velocity rate in Bacik was found to be −15 mm/yr, while the rates of subsidence in the other areas were trivial. The total number of buildings found in this area (7 km × 7 km rectangle) was 1258, in which 607 buildings were located in a gap area without InSAR velocity values. The 651 remaining buildings were mainly located in stable areas. Only 47 buildings were located in the critical red zone and the total area of these buildings was 12,000 m². If progressive land subsidence in the Balcik continues in the future, towards the industrial areas (south and west of the core of the subsidence), the risk of damage to these buildings may increase. Although the number of deep and semi-deep water wells illegally excavated in Balcik is unknown, the presence of illegal wells in the industrial and agricultural area was confirmed by regional administration. Table 2 also presents the total road length and the maximum and mean subsidence rates for Balcik.

Figure 9. (a) Balcik buildings overlaid on the mean vertical velocity (MVV) map; (b) roads overlaid on the MVV map. The locations of field observations are shown by the numbers in (b).

The field observations of Kadikoy and Balcik were carried out in November 2022, and signs of progressive land subsidence were confirmed in the visited areas. Figure 10 shows the locations of subsiding areas in the Kadikoy region. The majority of buildings had serious cracks in outside and inside walls. However, most of these cracks were repaired before our field visit (Figures 10(1, 5)). As a result of progressive land subsidence, most outer doors cannot be locked (Figures 10(3)). Short interviews with local people and the administrative office in Kadikoy revealed that the oval-shaped area near the water canal is filled land. Thus, the land subsidence due to soil consolidation and traffic effects near the football stadium (Kadikoy stadium) could be the main reasons for subsidence in the urban area. However, the subsidence in Balcik region could be different. Figure 11 shows the locations of subsiding areas in Balcik. Most of the affected buildings are in rural areas (Balcik village). This could be partially because of illegal water extraction in this area, since most of the lands in this area are muddy, and their fertility is good enough for agricultural activities. Our field observations showed that because of local subsidence, cracks are continuously appearing on the inner and outer walls of buildings in this area, despite the local people fixing these cracks from time to time. This subsidence area is also located near an industrial area. Thus, this is a call for action to make necessary decisions, in order to avoid possible damage to large infrastructures west and south of the subsidence zone. Table 3 shows the geographic location of observed points in the Kadikoy and Balcik regions. Additionally, we also analyzed 7 co-polarized (VV) ascending (orbit 58) images of Sentinel-1 from 1 May 2015 to 18 March 2021 using stacked InSAR technique for the study area. Detailed information of the Sentinel-1 images are presented in Table A1. The mean vertical velocity maps deduced from the ALOS-2 and Sentinel-1 datasets present similar patterns for Kadikoy region (Figure 12(a)). The correlation between vertical velocities of the two maps for the location of 595 buildings is 0.84 (Figure 12(b)).

Figure 10. Subsidence evidence in different locations in the Kadikoy region. The locations of the images are shown in figure 8(c).

Figure 11. Subsidence evidence in different locations in the Balcik region. The locations of the images are shown in figure 9(a).

Figure 12. (a) Sentinel-1 mean vertical velocity in the Kadikoy region; (b) scatter plot of mean vertical velocities.

Table 3. Geographic coordinates of the observed areas in Kadikoy and Balbik (see locations in Figure 8(c) and Figure 10).

Site number (Kadikoy)	Latitude (°)	Longitude (°)	Site number (Balcik)	Latitude (°)	Longitude (°)
1	40.994	29.040	1	40.872	29.429
2	40.993	29.042	2	40.873	29.429
3	40.993	29.041	3	40.875	29.431
4	40.992	29.039	4	40.875	29.430
5	40.986	29.037	-	-	-
6	40.989	29.040	-	-	-
7	40.993	29.041	-	-	-
8	40.986	29.035	-	-	-

4.2. Coherence changes

Figure 13 shows the coherence histograms over Kadikoy and Balcik for the near-to-one-year temporal pairs: time span one from 16 November 2014 to 15 November 2015 (364 days); time span two from 15 November 2015 to 30 October 2016 (350 days); time span three from 30 October 2016 to 29 October 2017 (364 days); time span four from 20 October 2017 to 23 December 2018 (420 days); time span five from 17 March 2019 to 15 March 2020 (364 days); and time span six from 15 March 2020 to 14 March 2021. The mean coherence values in Kadikoy for time spans one, two, three, four, five and six are 0.605, 0.610, 0.620, 0.580, 0.620 and 0.631, respectively. For Balcik, the mean coherence values for time spans one, two, three, four, five and six are 0.547, 0.558, 0.555, 0.543, 0.545 and 0.549, respectively. Although the general behavior of the histograms is similar, the peak position of time span four (black curve in Figure 13) is shorter than the other time spans. This shows that time span four in both regions has a lower value, due to its longer time interval.

Figure 13. (a) Statistical results of coherence changes for the SAR pairs with a time interval near to one year for Kadikoy and (b) for Balcik.

In Kadikoy, for time spans one and two, the correlation coefficient is 0.738; for time spans two and three, the correlation coefficient is 0.739; for time spans three and four, the correlation coefficient is 0.706; for time spans four and five, the correlation coefficient is 0.629; and for time spans five and six, the correlation coefficient is 0.726. For the Balcik area, for the same time spans, the correlation coefficients are 0.734, 0.770, 0.765, 0.697 and 0.715, respectively. Low correlation coefficients between time span three with time span four, and also time span four with time span five, are not surprising; since the temporal gap increases, the correlation decreases because of the urban/land use changes and vegetation growth.

4.3. RGB and PCA visualization

Using RGB coherence combination (Karimzadeh & Matsuoka, 2018), we tracked surficial changes in Istanbul. As shown in Figure 14, RGB combinations in different time spans were applied. Figure 14(a) shows the changes between 2014.11.16 and 2016.10.30 for almost two years, when the Yazuz Sultan Selim Bridge and the North Marmaray highway (Otoyol 7) were under construction. As a result of this, the bridge and highway are shown in black color, which means that all the channels were at the minimum values, or the values continuously were changing in the corresponding time span. In Figure 14(b), the RGB color composite covers between 2014.11.16 and 2021.03.14. The highway (except the southeast segment) is bluish, which means that the blue channel (coherence values from 2020.03.15 to 2021.03.14) had a large contribution, or the major parts of the bridge and the highway were made in this period, with a minimum contribution from the green and red channels. However, in Figure 14(c), the contribution of the green channel (coherence values from 2016.10.30 to 2017.10.29) manifested itself by changing the color of the bridge and highway from blue to cyan. In Figure 14(d), the color of the bridge and highway are “as white as” the built-up areas, because all of the channels cover from 2017.10.29 to 2021.03.14. This means that construction of the bridge and the highway was almost completed in this period, and the area was not changing anymore.

Figure 14. (a-d) RGB color composite for SBAS coherence pairs with one-year time spans.

As shown in Figure 15(a), PC1 (mainly urban regions) almost explains 92% of the data. PC2, PC3, PC4, PC5 and PC6 explain 5%, 2%, 0.5%, 0.5% and 0.0% of the data, respectively. Figure 15(b), shows that PC2 belongs mainly to the bridge and highway. The object-based image analysis shows that in the covered footprint, the total length of the constructed road and the bridge together is 55 km.

Figure 15. (a) Principal components (PCs) for the coherence pairs. Red box shows the location of Yavuz Sultan Selim Bridge; (b) RGB color composite of the first three PCs.

5. Discussion

The SBAS PALASR-2 analysis over the study area revealed new displacement areas. The reliability of the subsidence areas was confirmed by field investigations, and a comparative analysis using Sentinel-1 C-band images was carried out. In this study, the changes such as cracks in the walls are different from the earthquake damage. For example, in EMS-98 (European Macroseismic Scale) the damaged buildings with cracks on the walls are mainly categorized in Grade1. This means that all the damaged buildings due to the ground subsidence fall at the same grade of damage. In addition, the damaged buildings (i.e. buildings with cracks on the wall) are quickly repaired by the owners and the large temporal gap of the SAR images makes the establishment of an empirical relationship difficult. ALOS-2 satellite has a limited revisit frequency, meaning that it does not acquire data for a specific location as frequently as other active SAR satellites such as Sentinel-1. This limitation can affect our ability to monitor dynamic processes over short acquisition dates and further analysis such as X-band InSAR analysis may reveal other aspects of land subsidence in the urban areas. Although there are PALSAR-2 datasets available for Istanbul and Kocaeli, the accessibility of PALSAR-2 data for different regions or time periods may be limited due to factors like satellite orbit, acquisition scheduling, or operational considerations.

The PALSAR-2 images provided higher spatial resolution compared to the others, which was helpful to create a database for Istanbul and Kocaeli, in terms of the past, current and ongoing changes in the earth, especially in the urban-scale projects. For example, its efficiency for monitoring of Yavuz Sultan Selim Bridge before, during and after construction were evaluated with high-resolution coherence images. Coherence PCA provided invaluable information of L-bands regarding the amount of work that had been done by construction companies. Thus, monitoring the performance of different construction companies in terms of volume of work is possible. Dense vegetation around the highway was not problematic for producing high-quality coherence maps, since the penetration of PALSAR-2 images was at an acceptable level. The average coherence values in whole analyzed pairs were more than 0.5, despite the large time intervals. Full coverage of the study area with high spatial resolution remains a question for further investigations. Since there is a trade-off between spatial resolution and spatial coverage, here we used StripMap mode and high resolution, but the ALOS-2 footprint (i.e. track 188) was not able to cover the entire Marmaray highway. On the other hand, the resolution of ALOS-2 ScanSAR mode or Sentinel-1 may be unsatisfactory for calculations of changes over the highway. Further synergetic studies in the future may satisfy the purpose of this research for large areas.

6. Conclusion

In this study, we demonstrated the performance of ALOS-2 PALSAR-2 images for monitoring infrastructures such as roads and railways in Istanbul and Kocaeli. We applied SAR interferometry and SAR intensity analysis for 14 acquisitions that provided surface velocity and change indices for a large area. This approach together with an updated geodatabase allows us to classify the displacement behavior (e.g. mean, maximum and minimum displacements), and total length of displacement for the entire study area. In the central parts of the oval-shaped region in Kadikoy, the maximum vertical velocity rate was −30 mm/yr, whereas the rates of subsidence in the marginal parts were less than −30 mm/yr. The Kadikoy region (Figure 7) comprised a total of 9691 buildings. Among these, 1754 buildings (parcels) were situated in a gap area without InSAR displacement values, and there were 594 buildings located in areas experiencing high subsidence rates, exceeding −5 mm/yr. The total area of the critical buildings was 121,875 m² in Kadikoy. In Bacik, the maximum vertical velocity rate was −15 mm/yr. The total number of buildings in this area (Figure 7) was 1258, with 607 buildings in a gap area lacking InSAR velocity values. The remaining 651 buildings were mainly situated in stable zones. Only 47 buildings were in the critical red zone. Continued land subsidence towards industrial areas may pose a risk to these structures. Critical construction such as archeological sites and large structures/infrastructures can be selected automatically in the GIS environment based on attributes and descriptions. For example, if the coherence threshold or deformation rate of a certain region is available, the total number of constructions at risk can be initially identified using the floor area, proximity to the possible landslides and multicriteria decision systems.

Thus, the stakeholders in an urban area can derive some benefits from the results of the proposed method. Here are some ways stakeholders can benefit from the results:

1. Enhanced risk assessement and updates: the results facilitate the definition and creation of the risk maps for the experts who are interested in enriching their damage/risk maps through remote sensing techniques.

2. Urban planning and infrastructure development: InSAR results identify areas prone to subsidence, landslides, or other ground movements as shown in Figure 7. Stakeholders, such as local governments and disaster management agencies, can use these insights to plan evacuation routes, and develop strategies to mitigate potential hazards.

3. Infrastructure maintenance and management: InSAR results highlight infrastructural changes in roads and bridges (as shown in Figure 14). Municipalities and transportation authorities can use these insights to prioritize maintenance efforts, identify areas in need of repairs, and ensure the safety of residents and travelers.

Overall, the results offer an enhanced understanding of ground movement and urban changes, enabling stakeholders to make informed decisions and contribute to achieving sustainable development goals.

Acknowledgements

We thank JAXA for providing the PALSAR-2 images under the agreements via the 3rd Research Announcement on the Earth Observations (EO-RA3) for principal investigator (PI) Sadra Karimzadeh (ER3A2N105).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available on request from the corresponding author.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/20964471.2024.2320466

Additional information

Funding

This research was funded by the TUBITAK #2221 project and the University of Tabriz, and the Japanese Society for the Promotion of Science (JSPS) Grants-in-Aid for Scientific Research (KAKENHI), grant number #23H01654.

Notes on contributors

Sadra Karimzadeh

Sadra Karimzadeh received the Ph.D. degree in engineering with a certificate of commendation from Kanazawa University, Kanazawa, Japan, in 2015. He was a Postdoctoral Researcher with the Department of Geographical Information System (GIS) and Remote Sensing, University of Tabriz, Tabriz, Iran, in 2016. From 2016 to 2018, he was a JSPS Researcher with the Department of Architecture and Building Engineering, Tokyo Institute of Technology, Tokyo, Japan. From January 2019 and July 2019, he was a researcher with Geoinformatics Unit, RIKEN, in Japan. He was also an invited research scientists in Gebze Technical University, Turkey from October 2022 to September 2023. He is currently a research fellow in Tokyo Institute of Technology and an Assistant Professor with the Department of GIS and Remote Sensing, University of Tabriz, Tabriz, Iran. His research interests include SAR interferometry, radar imagery for earthquake damage assessment, infrastructure monitoring using SAR, GIS applications for disaster management, and seismic site characterizations based on GIS and remote sensing technologies. Dr. Karimzadeh is a Reviewer for most of the international journals in remote sensing and earthquake engineering.

Abdullah Can Zulfikar

Abdullah Can Zulfikar received the M.Sc. and Ph.D. degrees in civil or earthquake engineering. He is now with the Disaster Management Institute, Istanbul Technical University, Turkey. His research interests include seismic hazard, risk assessment, earthquake risk reduction technologies, such as near-real-time damage/risk assessment, and early warning systems.

Masashi Matsuoka

Masashi Matsuoka received the B.S. degree in architectural engineering from the Muroran Institute of Technology, Muroran, Japan, in 1990, and the M.S. and Ph.D. degrees from the Tokyo Institute of Technology (Tokyo Tech), Tokyo, Japan, in 1992 and 1996, respectively. He is currently a Professor with the Department of Architecture and Building Engineering, Tokyo Tech. He was with Tokyo Tech, Remote Sensing Technology Center of Japan, Earthquake Disaster Mitigation Research Center, Kobe, and the National Institute of Advanced Industrial Science and Technology, Tsukuba, prior to returning to Tokyo Tech. He is a Member of the Earthquake Engineering Research Institute and the Architectural Institute of Japan. His current research interests include earthquake engineering, geomorphology, GIS, and application of remote sensing technology in disaster management. Dr. Matsuoka is a Member of the Earthquake Engineering Research Institute, Tokyo, Japan and the Architectural Institute of Japan, Tokyo, Japan.