Quantifying horizontal and vertical movements in Ho Chi Minh city by Sentinel-1 radar interferometry

Ho Chi Minh City (HCMC), the most populous city and the economic center of Viet Nam, has faced ground subsidence in recent decades. This work aims at providing an unprecedented spatial extent coverage of the subsidence in HCMC in both horizontal and vertical components using Interferometric Synthetic Aperture Radar (InSAR) time series. For this purpose, an advanced InSAR technique PSDS (Permanent Scatterers and Distributed Scatterers) was applied to two big European Space Agency (ESA) Sentinel-1 datasets composed of 96 ascending and 202 descending images, acquired from 2014 to 2020 over HCMC area. A time series of 33 Cosmos SkyMED images was also used for comparison purpose. The combination of ascending and descending satellite passes allows the decomposition of the light of sight velocities into horizontal East-west and vertical components. By taking into account the presence of the horizontal East-west movement, our finding indicates that the precision of the decomposed vertical velocity can be improved up to 3 mm/year for Sentinel-1 data. The obtained results revealed that subsidence is most severe in areas along the Sai Gon river in the northwest-southeast axis and the southwest of the city with the maximum value up to 80 mm/year, consistent with findings in the literature. The magnitude of horizontal East-West velocities is relatively small and a large-scale westward motion can be observed in the northwest of the city at a rate of 2-5 mm/year. Together, these results reinforced the remarkable suitability of ESA’s Sentinel-1 SAR for subsidence applications even for non-Europe countries such as Vietnam and Southeast Asia.


Introduction
Over the last 20 years, the population of Ho Chi Minh City (HCMC) has doubled to 9 million [1], making the city one of the fastest-growing economic centers of Asia [2]. However, land subsidence, sealevel rise, and extensive delta-based water infrastructure have raised concern due to global climate change [3]. In this context, monitoring the spatial extent of surface subsidence is critical to mitigating related hazards, since it may affect population, infrastructures, and buildings.
Since the 2000's, time series SAR interferometry (InSAR) has been exploited to map large scale surface subsidence with high accuracy from space [4,5]. Particularly, the recently advanced technique PSDS (permanent scatterers and distributed scatterers) becomes the main InSAR tool for many deformation applications due to its superior performance [6,7]. The radar interferometric measurement can only detect millimetric motion in the direction linking between the SAR's sensor and the sensed area (so-called light of sight -LOS) [8]. This leads to the challenging for InSAR measurements to do interpretation and communication [9]. To meet this challenge, there is a common assumption on neglecting horizontal component motions when projecting the LOS measurement to vertical direction [9,10]. However, it is possible to decompose it to horizontal and vertical products in certain areas of interest where two LOS observations are available (e.g., from ascending and descending satellite passes) [9]. Unfortunately, this is not the case in Vietnam and many non-European countries. Up to now, no precise information is available concerning the decomposed horizontal and vertical motions from InSAR measurements in Vietnam.
In Vietnam, the InSAR technique is an interesting approach applied in many deformation studies using X-band [11][12][13], C-band [14,15], and L-band [3,16,17]. Particularly, in [3], the authors have reported a detailed subsidence map that covers urban regions of the Ho Chi Minh city area for the period 2006-2010 using L-band ALOS PalSAR data. It was reported that subsidence areas are severe along the Sai Gon River and in the southwest of the city where their subsurface soils are mostly Holocene silt loam. However, the study so far has estimated vertical subsidence using the classical assumption on the projection of the LOS measurement. Fortunately, since 2014, the European Space Agency (ESA) Sentinel-1 constellation provides open-access and global SAR data with over 10 TB and more than 4700 products per day [18,19]. Big Data of SAR images regularly acquired from such satellite mission with short repeat cycle (e.g., 6 days over Europe and 12 days for the rest of the world), and a large swath width of 250 km on the ground open new opportunities for the monitoring of the large scale Earth's surface. For example, to cover the HCMC from 2014 to 2020, there are over 200 images in descending and 100 images in ascending passes, respectively. Thus, Sentinel-1 Big Data offers the best chance to quantify the horizontal and vertical movements. In this paper, we aim to provide a better understanding of the capabilities of Sentinel-1 Cband in estimating decomposed horizontal and vertical components for subsidence phenomena in Ho Chi Minh city. This paper is organized as follows: the study site and SAR data are shown in Section 2; the PSDS InSAR algorithm and the combination of two satellite passes are presented in Section 3; the results are examined in Section 4; and finally the discussion is provided in Section 5.

Study site
The study area is the whole area of Ho Chi Minh City, Vietnam. The population of this megacity is over 9 million (in 2019), which greatly contributes to the economic, financial, cultural, and scientific developments of the country [1]. It is located at 106.65 E longitude and 10.85 N latitude, southeastern Vietnam. The city is surrounded by the Sai Gon River system. The whole study area is approximately 90 km x 70 km, including the HCMC area and its surrounding, see Figure 1. The city consists of three regions based on the urban growth: (1) the core region, where urbanization is completed; (2) the fringe region, where urban growth is undertaking place; and (3) the suburban region for the rest of the city [20].

SAR data
The SAR stacks are comprised of C-band Sentinel-1 (∼ 5.55 cm) and X-band Cosmos SkyMED (∼ 3.13 cm) data. Their selected footprints are shown in Figure 1. For Sentinel-1 data, both ascending pass track 128 (Sentinel-1 A128) and descending pass track 18 (Sentinel-1 D18) are considered. They are imaged in TOPS mode (i.e., Terrain Observation by Progressive Scans [21]), whereas Comos SkyMED stack is imaged in stripmap mode. In each stack, each Single Look Complex (SLC) image was resampled on a reference grid. The characteristic's three datasets is reported in Table 1 and Figure 2 for the detailed information. For Sentinel-1 TOPS mode, the coregistration is done by using 6 bursts of swath 1 for both A128 and D18 tracks. These bursts are then merged and cropped to generates stacks of coregistered SLC datasets as shown in Figure 1. The subset operator is not applied to Cosmos SkyMED scenes to maximize their coverages. Although we will exploit all possible information in the PSDS technique as shown in Section 3.1, we selected the Cosmos SkyMed reference date in the middle of the geometric baseline distribution to minimize mis-coregistration errors due to long baselines. For Sentinel-1, the baseline tube is relatively small, but the coregistation has been carried out carefully to avoid a phase jump in the azimuth direction [22].

Method
The rationale of InSAR time series is to take advantage of multiple flight tracks to minimize issues relative to atmospheric, spatial, and temporal decorrelations to estimate ground deformations robustly [4,23] . Let us suppose that N SAR images are available for a certain area of interest. The images are coregistered on a reference grid and phase contributions due to terrain topography and orbit have been compensated. The residual phase of the differential interferogram ϕ n is composed of the phase contributions related to residual topography, deformation, atmosphere, and noise [8]: where W is the wrapping operator; ϕ n = 4πbn ε z is the residual topographic phase; ϕ n is the deformation phase; n topo λsinθRn de f o ϕ aps is the atmospheric phase screen, which represents the signal delay due to weather conditions; and ϕ n is the phase noise due to temporal decorrelation, mis-coregistration, uncompensated spectral shift decorrelation, orbital errors, soil moisture, and thermal noise. ε z is the elevation error; bn is the normal baseline relative to the n-th image with respect to the reference image; θ is the local incidence angle; λ is the carrier wavelength; Rn is the (zero-Doppler) distance between the target and the n-th orbit acquisition.
The phase component due to the displacement can be usefully split into two components: where dn is the light of sight distance change of the target; v is the constant velocity; tn is the temporal baseline; and µNL is the phase term due to a possible non-linear motion. In Figure 3, an example of target motion between the n-th and reference acquisition time is shown. Our goal is to estimate the deformed phase from the interferometric phase. In details, given N − 1 differential interferograms, for each pixel, from Eq. 1, we can have a nonlinear system of N 1 equations and two parameters of interest ε z and v. The interferometric phase is naturally given as a wrapped number and needs to be unwrapped (i.e., by adding a certain integer ambiguity number 2π [24]) to obtain ε z and v. In the InSAR literature, there are two strategies to unwrap the Eq. 1. In the first approach, the interferograms are spatially multilooked, filtered, and unwrapped [25]. In this way, the model in Eq. 1 becomes linear and the inversion is straightforward by using a Singular Value Decomposition [25,26]. To account for signal decorrelations, only several subsets with short spatial and temporal baseline (i.e., Small Baseline Subset -SBAS) possible are considered. The main drawback is that the result can be biased because the unwrapping error is introduced (i.e., by adding a wrong integer ambiguity number 2π during the phase unwrapping). Furthermore, it may be not suitable for studying at local-scale deformations [25,27] and it also has a non-closure phase biased due to working on multilook imagery [28][29][30]. In contrary, in the second approach, the analysis can be carried out using the wrapped phase in model Eq. 1, avoiding the introduction of the unwrapping error. Indeed, the unwrapping and estimation of parameters ε z and v are possible only if the signal-to-noise ratio (SNR) is high enough. In detail, only stable point-like pixels, characterized by low noise, are considered in the analysis. Such stable targets are named Permanent/persistent Scatter (PS) which is a dominant scatterer within resolution elements [4,5,23]. In practice, the relative ε z and v for every PS pair are jointly estimated by maximizing a complex ensemble coherence in time [7]. We are then able to calculate ε z and v over the whole network by integrating their relative values from pairs of neighboring pixels for a reference pixel. The robustness of this method is highly dependent on PS availability. The PS technique works on a time series dataset, − resulting in an averaging effect. This can lead to a loss of information because it can not be able to consider slowly decorrelating targets (i.e., Distributed Scatterers (DS)) which are usually corresponding to non-cultivated and bare land areas. Recently, the motivation of using both PS and DS is emphasized, aiming to maximize the amount of information available from interferometric phases and to improve the density of selected measurement points [6,[31][32][33][34].

The PSDS technique
The main limitation of DS targets is that their signal-to-noise ratio is typically low because of geometrical and temporal decorrelations [35]. The principle idea of the combination of PS and DS targets is to enhance SNR of DS so that it can become PS. To do it, a Maximum Likelihood estimation framework is employed to extract the best estimates N − 1 DS interferometric phases from all the N(N − 1)/2 interferometric phases available of N images [6,[31][32][33]. For this purpose, to select a DS candidate, a brotherhood that contained similarly statistical homogeneous pixels (SHP) are identified. To test their similarity, a traditional two-sample Kolmogorov-Smirnov test can be applied on the reflectivity time series of the interested DS candidate and neighbor pixels inside a certain window [6]. Base on this SHP family, a DS sample coherence matrix (Γ ) (i.e, which contains N(N − 1)/2 phase values) is generated and is exploited by Phase Triangulation Algorithm [6] or Phase Linking [36] to estimate for N 1 linked phase (λ ) as : where Λ = exp(jλ); H is Hermitian conjugation; and • indicates the Hadamard entry-wise product. After this procedure, the DS phase values are filtered, putting the linked phase N − 1 phase values λ equaled to the quality of the PS. Hence, the enhanced DS targets can be combined with PS using the same PS algorithm. For more details on the technique, the reader can be found in [7]. The PSDS, nowadays, is a replacement of PS technique due to its stable performance particularly in nonurban areas [6]. Recently, this technique is applied as the main tool for many deformation applications [3,13,[37][38][39]. In this work, we use the TomoSAR platform to carry out PSDS processing as demontrated in [3,13,15] with an update on using Baumgartner-Weiss-Schindler test statistic for DS selections [40]. As supported in [41], two-sample test based on the Baumgartner-Weiss-Schindler statistic is more robust than traditional Kolmogorov-Smirnov [6] and Anderson-Darling tests [37].

The combination of ascending and descending satellite passes
The interferometric observation is naturally only sensitive to motion in the SAR LOS direction. However, LOS measurements can result in some of ambiguities. By combining both ascending and descending information (see Figure 4), we can decompose the LOS velocities into East-west and vertical components as follow [9]: where vlos is the LOS velocity, α is the azimuth of the satellite heading vector, and vE, vN, vU are the East, North, and vertical components of motion, respectively. The Eq. 4 has three unknowns but we have only two observational constraints from ascending and descending LOS measurements. It is then possible to calculate the three-dimensional velocity field only if there is an additional constraint. We observe that in the case of our Sentinel-1 dataset (i.e., swath 1), both ascending and descending directions are quite steep angles (about 34 degree from the normal). As a result, the north component is insensitive and we can only measure horizontal motion along the East-west axis more accurately. Consequently, we add the additional constraint to the Eq. 4 is that the North-south motion is assumed to be zero. In this way, we can calculate velocity in East-west and Up-down direction. For this purpose, in this work, both ascending and descending measurements were resampled to a common grid with a 50 m spatial resolution. We note that there is no interpolation performed at the grid where PSDS measure points are not available.

Metric assessment
For the comparison analysis, the velocities derived by projecting LOS to vertical are compared with the decomposed vertical velocity calculated by using both ascending and descending passes as in section 3.2. The assessments are carried out through standard deviation, bias, correlation coefficient (R 2 ) and Root Mean Square Difference (RMSD): where Yi and Ȳ i are the estimated variables and their mean values, Xi and X i are the reference measurement variables and their mean, and m is the number of data.

Surface subsidence results
The Sentinel-1 A128, Sentinel-1 D128, and Cosmos SkyMED data were treated by using PSDS processing. In this section, to have the same metric for comparison among velocity's results, LOS measurements are converted to vertical direction by dividing the cosine of the incidence angle (i.e., neglecting the presence of horizontal motions). The reference point was set at the same location in the geological sandstone area in District 1, at 10.782N, 106.698E, as suggested in [3]. The result was reported in Figure 5, reflecting similar spatial distributed subsidence patterns. This shows that subsidence is most severe in areas along the river and the southwest of the city, consistent with findings reported in [3].
To make it comparable among three velocities, a zoom focused on a red-bordered rectangle area (see Figure 5c) is shown in Figure 6a, b and c. More than 35000, 29000, and 87000 PSDS measuring points in Sentinel-1 A128, Sentinel-1 D18, and Cosmos SkyMED datasets, respectively, were identified within an area of about 42 km 2 . The lower measurement points in Sentinel-1 is expected due to the higher spatial resolution of Cosmos SkyMed (see Table 1).
A buffer of 100 m diameter centered on a certain line AB in Figure 6c was associated with a cluster of PSDS. All the PSDS measuring points inside the 300 m bin cluster were used to calculate the average velocity for each dataset. The profile was shown in Figure 6d. It is noted that the maximum value is up to -80 mm/year, representing a highly vulnerable area.
Finally, to appreciate our long term datasets, a time series plot is reported in Figure 6e. Two Sentinel-1 datasets are referred to the same starting temporal reference in Comos SkyMED (i.e., 30 March 2014). The displacement history at two points representing moderate (P1) and severe (P2) subsidence is highlighted. It is clear that subsidence is primarily linear, which agrees well with findings presented in [3].

East-west and vertical velocities
The ascending and descending velocities (see Figure 5) can provide a detailed picture of Ho Chi Minh city area subsidence. However, it maybe some of ambiguities associated with LOS measurements due to non-negligable influences of horizontal displacement. To remove it, we decompose the LOS velocities into East-west and vertical components using Eq. 4. The decomposed horizontal East-west and vertical velocities are reported in Figure 7. To facilitate visualization of the subsidence results, Figure 8 and 9 are highlighted in two severe areas. It is expected to observe that there is no significant difference in Figure  7a with respect to Figure 5 as tectonics are not found in southern Vietnam [42].
The resulting decomposed East-west velocity is more difficult to interpret than the vertical component. It is evident that the magnitude of east-west components is relatively small with respect to vertical contributions. The most prominent feature in East-west velocity is the large-scale westward motion in the northwest of the city at a rate of 2-5 mm/year.
To quantify the relationship between horizontal and vertical motions, we plotted their joint distribution and the associated standard deviation in Figure 10. The maximum magnitude of East-west velocity is about 10 mm/year. The joint distribution is normalized between 0 and 1 to better do the interpretation. It is evident that they are independent random variables (i.e., R 2 = 10 −3 ). For vertical velocities varying from -40 mm/year to zero, we can observe that the mean of East-west velocity is approximately close to zero, whereas its standard deviation is decreased from 4.5 mm/year to 0.9 mm/year.

Comparison vertical velocities
The impact of horizontal motion on the projection of LOS measurements into the vertical direction is investigated in this section. To do it, all measurements (see Figure 5) were resampled to a common grid with a 50 m spatial resolution as the decomposed vertical velocity product (see Figure 7a). A total of 101620 50-m pixels is found for the same area of interest between Sentinel-1 and Cosmos SkyMED data. The pixel-wise joint distributions between decomposed and projected vertical results were shown in Figure 11. The distribution is normalized between 0 and 1 to facilitate visualization. Strong correlation coefficients (i.e., R 2 ≥ 0.67 for 101620 samples) were found. The RMSD is from 2.9 mm/year to 3.1 mm/year for Sentinel-1 and 3.9 mm/year for Cosmos SkyMED.
To quantify the detail of the magnitude of error, the standard deviation and bias were calculated and reported in Figure 10. The results of projecting LOS to vertical would erroneously indicate that the bias value of ±2 mm/year and ±4 mm/year can be made for Sentinel-1 and Cosmos SkyMED datasets, respectively.

Discussion
In this work, we showed that horizontal East-west and vertical motion can be decomposed in Ho Chi Minh City area using the Sentinel-1 PSDS InSAR technique. The relative comparison with Cosmos SkyMED indicated good agreement with the strong correlation for vertical velocities. The unprecedented spatial extent coverage of the vertical velocity reflects similar spatial subsidence patterns with respect to previous studies. For example, along the river in the northwest-southeast axis and the southwest of the city, subsidence is most severe at a rate up to 80 mm/year. The study shows that the magnitude of horizontal East-west velocities is relatively small, showing the large-scale westward motion in the northwest of the city at a rate of 2-5 mm/year. Together, these results confirm the remarkable suitability of ESA's Sentinel-1 SAR for continuously monitoring ground subsidence applications with the frequent and reliable coverage even for non-Europe countries such as Vietnam and Southeast Asia.    First of all, this is the first work (to our knowledge) using data from ascending and descending measurements to provide decomposed horizontal East-west and vertical velocities in the Vietnam area. This is not possible for many non-Europe countries before ESA's Sentinel-1 era. In fact, in SAR archive data, there are very few ESA SAR data (e.g., C-band ERS-1/2 and ENVISAT ASAR) exist that are suitable for such applications. This is the reason why previous studies are based on L-band from Japan Aerospace Exploration Agency ALOS PalSAR [3,16] and X-band (e.g., Cosmos SkyMED from Italian Space Agency and TerraSAR-X from German Space Agency) [12,13]. They were all based on the assumption of neglecting the presence of horizontal motions, projecting LOS measurements into the vertical values.
It is important to point out that our subsidence results were analyzed only on pixels where their time series signals are stable or high signal-to-noise ratio (see Section 3). It does not mean the deformation does not happen at other pixels. For example, although we mapped subsidence for the whole city area, there is an area that strongly lacks PSDS points in the Southeast part (see Figure 5), which is mainly Can Gio mangrove forest's area. This is clearly due to temporal decorrelation caused by vegetation [35,43,44]. Future missions on longer wavelength such as NISAR L-band [45], Tandem-L [46], ROSE L-band [47], and BIOMASS P-band [48,49] are expected to minimize temporal decorrelation and hence yield better performances in vegetation areas.
To provide an independent measurement for comparison, we exploited Cosmos SkyMED X-band SAR data as shown in Figures 5c, 6, 11c, and 12. It is worth pointing out that this stack is characterized by higher spatial resolution and shorter wavelength (see Table 1). Consequently, denser pixels were identified (i.e., 87000 pixels vs 29000-35000 pixels from Sentinel-1 data, see Section 4.1). This is the reason why the RMSD, standard deviation, and bias of projected Comos SkyMED results are slightly bigger than the projected Sentinel-1 A128 and D18 measurements. In other words, it does not mean Comos SkyMED data is of minority quality to Sentinel-1. In fact, with X-band (e.g., Comos SkyMED and TerraSAR-X), it is possible to study at a target-scale level such as bridges and buildings [50]. Particularly, Ho Tong Minh et al., [13] have demonstrated the feasibility of subsidence estimates by X-band InSAR even in conditions impacted by a strong atmosphere in northern Vietnam, reinforcing the scientific potential for a new Vietnamese X-band SAR space-borne mission (JV-LOTUSat).
Finally, given ascending and descending SAR data systematically available from Sentinel-1, nowadays, neglecting the impact of horizontal motion in converting LOS measurements into the vertical direction is generally not recommended. We have shown that the improvement of 3 mm/year for Sentinel-1 data (see Figure 11a and b) is evident with the combination. However, since only the East-west and vertical components are estimated, a slight error can be introduced because of neglecting the North-south component. Fortunately, this error is very small (i.e., 0.1 mm/year [10]) in many isotropic deformation cases, whereas it becomes critical for non-isotropic deformation scenarios (e.g., landslide and earthquake deformation phenomena) [9]. Finally, the Sentinel-1 constellation is the only system that can provide the systematic survey, open-access, and long term SAR data. Therefore, it is recommended for using Sentinel-1 as a SAR data resource for operational subsidence monitoring tasks, particularly for large-scale, delta-wide, and nation-wide coverages.
To conclude, the PSDS InSAR technique is a powerful approach to provide unprecedented spatial extent ground subsidence in Ho Chi Minh City area using Sentinel-1 SAR data. Horizontal East-west and vertical motions can now be mapped on a large scale with the reliable coverage provided by the ESA's Sentinel-1 constellation. Thus, it is recommended to develop operational projects such as delta-wide and national wide initiatives for long term motion monitoring. Consequently, this will foster the developments of downstream solutions by the final beneficiaries: infrastructure managers (e.g. roads, railways, pipelines), civil engineering companies, regional and municipality administrators, and public policymakers.