Monitoring Ground Subsidence in Hong Kong via Spaceborne Radar: Experiments and Validation

Abstract

The persistent scatterers interferometry (PSI) technique is gradually becoming known for its capability of providing up to millimeter accuracy of measurement on ground displacement. Nevertheless, there is still quite a good amount of doubt regarding its correctness or accuracy. In this paper, we carried out an experiment corroborating the capability of the PSI technique with the help of a traditional survey method in the urban area of Hong Kong, China. Seventy three TerraSAR-X (TSX) and TanDEM-X (TDX) images spanning over four years are used for the data process. There are three aims of this study. The first is to generate a displacement map of urban Hong Kong and to check for spots with possible ground movements. This information will be provided to the local surveyors so that they can check these specific locations. The second is to validate if the accuracy of the PSI technique can indeed reach the millimeter level in this real application scenario. For validating the accuracy of PSI, four corner reflectors (CR) were installed at a construction site on reclaimed land in Hong Kong. They were manually moved up or down by a few to tens of millimeters, and the value derived from the PSI analysis was compared to the true value. The experiment, carried out in unideal conditions, nevertheless proved undoubtedly that millimeter accuracy can be achieved by the PSI technique. The last is to evaluate the advantages and limitations of the PSI technique. Overall, the PSI technique can be extremely useful if used in collaboration with other techniques, so that the advantages can be highlighted and the drawbacks avoided.

1. Introduction

Ground subsidence has gradually become one of the most common threats to lives and properties in the fast-booming mega-cities in China. Along with the rapid development of cities often comes the excessive withdrawal of underground resources and ceaseless new construction above and under the ground surface, all of which are potential risks that cause ground subsidence. Monitoring the stability over urban areas thus becomes very important for preventing ground subsidence and to ensure surface stability [1,2,3,4,5].

Up to now, the most common approach for observing and measuring ground subsidence remains traditional ground leveling and GPS measurements. On the one hand, when combining with the advanced GPS technique and the precise coordinates of control point grids, ground leveling makes it possible to directly measure ground motion with quite high accuracy. On the other hand, however, as cities expand rapidly over time and space, the traditional methods are now facing many challenges, among which are the shortage in manpower, the cost required for performing more measurements, as well as the inability to get full coverage of time and space over large areas. To overcome these shortcomings, a new technique that calls for all-time, wide-area and high-accuracy ground monitoring is needed, and that is where synthetic aperture radar interferometry (InSAR) comes in.

Since its first illumination of the Earth from the 1980s, the InSAR technique has been successful drawing enormous attention from scholars and researchers in the field of geoscience. Nowadays, this technology has already been widely acknowledged for its ability to continuously monitor ground motions and to recover the digital elevation model (DEM) [6,7]. Furthermore, for the purpose of overcoming the two major limiting factors of InSAR, namely geometric and temporal decorrelation along with the atmospheric phase screen (APS) disturbance, the persistent scatterers interferometry (PSI) technique, also known as the multi-temporal InSAR (MTI) technique, was proposed and has since become the core concept in precisely measuring ground displacement from a stack of multi-SAR-images [8,9,10]. Among many of the technique’s strengths, the most important is detecting displacement within millimeter to sub-millimeter accuracy. In the past, a number of projects were carried out to demonstrate the sensitivity of PSI technique under ideal circumstances. For example, the experiment carried out by Ferretti [11], using two dihedral corner reflectors (CR) with a specifically-designed size and placed on the top of a building roof (so that a good signal-to-clutter ratio could be guaranteed), mounted on two sliding platforms controlled by a robust clamping system (so that the CR can be manually and precisely adjusted by the millimeter) and close to each other (so that the atmospheric effect can be neglected), demonstrated the sub-millimeter accuracy of the PSI technique with the C-band synthetic aperture radar (SAR) dataset. Other similar experiments also showed a similar conclusion [12,13].

However, in a real scenario, the accuracy of the PSI technique can vary with a blending of different factors, including the accuracy of the satellite’s spatial baseline, atmospheric conditions [14], the physical nature of persistent scatterer (PS) points [15], the coherence between adjacent PS points [16], etc. For example, our area of interest (AOI), Hong Kong, consists mostly of ocean, vegetation and dense skyscrapers, together with high humidity all year round. From the radar point of view, the complicated ground features make it very difficult to generate high quality interferograms with high spatial resolution and wide spatial coverage if using the conventional InSAR technique. Even applying the PSI technique, a large portion of the ocean, mountainous area and isolated islands can mean a low temporal/spatial coherence, bringing down the accuracy of PSI technique. In addition, the fact that skyscrapers usually look like a complex mix of point-wise targets, but vegetated areas absolutely lack these brings more challenges to PSI processing. From the atmospheric point of view, the humid subtropical climate of the city and extremely high water vapor density in the air result in a greater extent of water vapor delay than in general, causing significant distortion in detected ground displacement. To successfully apply the PSI technique, the atmospheric water vapor delay must be carefully estimated and removed [1].

In this project, we try to carry out the experiment of generating a ground displacement map in the urban area of Hong Kong, China, with the PSI technique and compare the outcome with traditional survey methods so as to demonstrate the millimeter accuracy of PSI in practical applications. To prove the accuracy of the PSI technique, four corner reflectors (CR) were installed inside the AOI. In order to compare the deferential displacements of each CR with regard to the others, one CR was chosen as the reference and kept stable all of the time, while the others were either lifted or lowered manually to an extent of some millimeters, and their displacements were recorded by optical leveling. The accuracy of PSI would then be evaluated by comparing the displacement computed by PSI and recorded by optical leveling. The experiment, carried out in unideal conditions (since the CR validation test was located at a construction site on reclaimed land), nevertheless proved that millimeter accuracy can definitely be achieved. At last, the advantages and limitations of the PSI technique applied in this project were fully evaluated, aiming at providing some beneficial references for similar applications in the future.

This paper is organized as follows. Firstly, PSI processing over the Hong Kong urban area was carried out, bringing three outcomes: the reflectivity map, the 3D location of PS points with 1-m precision and the estimated millimeter displacement along the satellite line of sight (LOS) of PS points. For the second part, a set of CR was deployed for a time period of a half year for the purpose of validation testing, and the displacement values of CR from both ground leveling and PSI processing were compared. A discussion of the validation outcome and the assessment of the accuracy will follow at the end of the paper.

2. Study Area and Dataset

The study area of interest belongs to Kowloon, which is the majority region of urban Hong Kong, China. The area is marked with the red rectangle in Figure 1. In this study, a total of 73 TerraSAR-X (TSX) and TanDEM-X (TDX) images acquired between October 2008 and September 2012 was used. The SAR images are in the X-band with a wavelength of 31 mm and a frequency of 9.6 GHz, with a resolution of up to 2 m. The minimum interval between the two successive acquisitions is 11 days. The images were acquired at about 6:25 pm each day, along an ascending orbit with an incidence angle of approximately 37°.

Figure 1. The area of interest (AOI) for the study, which covers the majority of urban Hong Kong.

Figure 1. The area of interest (AOI) for the study, which covers the majority of urban Hong Kong.

Figure 2. This is the figure showing TSX and TDX images used in our analysis, where the x-axis indicates the temporal baseline and the y-axis indicates the perpendicular baseline. It is worth noting that nearly all of the baselines for TSX images are precisely controlled to within a distance of about 200 m, with only a few exceptions that are mostly TDX images, which requires a larger baseline for its DEM mapping mission. The color bar indicates the average spatial coherence of the interferogram given one master and one slave image. It gives a general idea that with a smaller spatial/temporal baseline usually comes a higher spatial coherence.

Figure 2. This is the figure showing TSX and TDX images used in our analysis, where the x-axis indicates the temporal baseline and the y-axis indicates the perpendicular baseline. It is worth noting that nearly all of the baselines for TSX images are precisely controlled to within a distance of about 200 m, with only a few exceptions that are mostly TDX images, which requires a larger baseline for its DEM mapping mission. The color bar indicates the average spatial coherence of the interferogram given one master and one slave image. It gives a general idea that with a smaller spatial/temporal baseline usually comes a higher spatial coherence.

In Figure 2, the interferometric configuration is shown. Each dot in the image represents a SAR acquisition in the normal baseline vs. the temporal baseline axis, and each line depicts one interferogram. The star-like graph is also usually the standard and optimized configuration for PSI processing. In this star-like connection, all images are referred to one master image taken on 11 June 2010. In general, the image acquired midway through the monitoring period will be selected as the master image. In general, the selected master will not impact the final output.

3. Methodology

The processing methodology can be divided into two major steps: PSI processing and CR validation testing.

3.1. PSI Processing

For a detailed description of the PSI approach, please refer to [9] and [8]. In short, the applied processing steps are as follows: master image selection, SAR data coregistration, generating reflectivity map and amplitude stability index map, persistent scatterers candidate (PSC) selection, multi-image sparse grid phase unwrapping, atmospheric phase screen (APS) estimation and removal, PS phase reading and displacement estimation. The movements of the PS targets are retrieved as a function of time with respect to the master image selected. By processing the phase of SAR images with the PSI technique, two results can be retrieved: the three-dimensional localization of PS points with metric precision and the estimated millimeter displacement of PS along the satellite line of sight. The process was conducted in SARPROZ [17] software (developed by Dr. Perissin, Purdue University, West Lafayette, IN, USA) and the final outcome was geocoded to and displayed in Google Earth. It is worth mentioning that selected PSCs in AOI are mostly man-made objects, including buildings (edges and corners on buildings), lamp posts, road signs and other concrete structures, the majority of which are considered typical examples that reflect the radar signal well [18].

3.2. Corner Reflector Validation Test

In order to evaluate the performance of PSI techniques, a CR validation test was carried out in our AOI. For the validation test, a set of reflectors were deployed, with one of them selected as the reference CR and kept stable, while the other reflectors were either lifted or lowered manually from a few to tens of millimeters, and their displacements were recorded by a survey team. We would then compare the displacement value calculated by PSI and measured by optical leveling. For this study, the validation period lasted five months from January 2012 to June 2012, and 11 images were used for validation testing. The details of the data can be found in Table 1. The CR prototype installed for the validation test is shown in Figure 3. Four CRs were installed for the test and were sequentially named CR4, CR3, CR2 and CR1, where CR4 was the reference.

Table 1. List of TSX and TDX images used for the CR validation test.

Scenes	Date	Perpendicular Baseline (m)	Type	Incidence Angle	Orbit	Image Mode	Local Visiting Time (Local Time)
1	6 February 2012	–238	TDX	37.32	ASCENDING	Strip map	06:25 pm
2	17 February 2012	–216	TDX	37.32
3	21 March 2012	–222	TDX	37.32
4	1 April 2012	–18	TSX	37.36
5	12 April 2012	49	TSX	37.36
6	23 April 2012	–114	TSX	37.35
7	4 May 2012	–19	TSX	37.35
8	15 May 2012	–131	TSX	37.35
9	26 May 2012	–11	TSX	37.35
10	6 June 2012	–179	TSX	37.35
11	17 June 2012	66	TSX	37.36

Table 1. List of TSX and TDX images used for the CR validation test.

Scenes	Date	Perpendicular Baseline (m)	Type	Incidence Angle	Orbit	Image Mode	Local Visiting Time (Local Time)
1	6 February 2012	–238	TDX	37.32	ASCENDING	Strip map	06:25 pm
2	17 February 2012	–216	TDX	37.32
3	21 March 2012	–222	TDX	37.32
4	1 April 2012	–18	TSX	37.36
5	12 April 2012	49	TSX	37.36
6	23 April 2012	–114	TSX	37.35
7	4 May 2012	–19	TSX	37.35
8	15 May 2012	–131	TSX	37.35
9	26 May 2012	–11	TSX	37.35
10	6 June 2012	–179	TSX	37.35
11	17 June 2012	66	TSX	37.36

Figure 3. Corner reflectors installed on site for the validation test.

Figure 3. Corner reflectors installed on site for the validation test.

The characteristics of different types and sizes of CR will generally have very different impacts on CR validation tests; thus, it is important to carefully design and deploy CR prototypes. Prior to this study, the authors have conducted comprehensive experiments for the design and deployment of an optimized prototype specifically for this case, and a detailed discussion was carried out in [19]. In summary, the CRs were designed as rectangular trihedrals that were 0.5 m by 0.5 m wide and 0.75 m high, installed on an adjustable base. The size of the CRs were calculated to balance between a reasonable radar signal reflectance and the weight. To further reduce the weight of reflectors and the impact of wind resistance, the metallic plates were perforated regularly with a circular hole array. For deploying, the CRs were installed in a region that presents a low backscatter to guarantee a good signal-to-noise ratio (SNR), thus guaranteeing a higher accuracy. The relationship between SNR, phase noise and the dispersion of the displacement measurements along the satellite line of sight (LOS) can be approximated as follows: [20]:

$${\sigma }_{LOS}=\frac{\lambda }{4\pi }·{\sigma }_{\varphi }\approx \frac{\lambda }{4\pi }·\sqrt{\frac{1}{2·SNR}}$$

(1)

For example, a CR that represents an SNR of 100 (20 dB) corresponds to a dispersion of the displacement in the LOS direction of about 0.17 mm at the X-band (λ = 3.1 cm). In addition, in order to avoid possible unexpected relative motions and to keep the atmospheric noise as limited as possible, the CRs were kept within a close distance (50 m) from each other [11].

The image on the right of Figure 4 shows the reflectivity map of the area where the corner reflectors were deployed, which is generated by taking the average of 73 SAR images. It is worth mentioning that, although the CRs only appear on 11 out of 73 images, the SNR values are so strong that a short period (five months) of signal availability can be clearly seen on the average over a four-year period.

Figure 4. Left: the location of CRs geocoded to the optical image. The location is near the sea that is identical for its low intensity in radar images. Right: the location and intensity of the four CRs on TSX images. From the image, we can see that the backscatter in that area is very low, and the intensity of the CRs is considerably high.

Figure 4. Left: the location of CRs geocoded to the optical image. The location is near the sea that is identical for its low intensity in radar images. Right: the location and intensity of the four CRs on TSX images. From the image, we can see that the backscatter in that area is very low, and the intensity of the CRs is considerably high.

3.3. Some Details about Ground Survey Data

The locations of four CRs were measured before and after the validation test by local surveyors. Three CRs were moved perpendicularly up or down from a few to ten millimeters, and their new heights were recorded. For measuring the coordinates of four CRs, the real-time kinematic (RTK) technique was used, the typical nominal accuracy for which is around 1 centimeter horizontally [21]. For lifting/lowering CRs, the surveyors carefully screwed the nuts that were installed on the base of CR to adjust the height and then measured the new height with the optical leveling method with the nearest control point available (in this case, the nearest control point is approximately 200 m away). The accuracy of height measured by optical leveling for this experiments is controlled to 1 mm [22].

3.4. InSAR Principles for Estimating CR Displacements

Before we estimate the displacement of the CRs, the other factors that might bias the displacement should be considered and removed. First of all, the interferometric phase can be decomposed into the following terms:

$$\Delta \phi =\Delta {\phi }_{displacement}+\Delta {\phi }_{height}+\Delta {\phi }_{atmosphere}+\Delta {\phi }_{noise}$$

(2)

In sequence, the four parameters are possible displacement, height of targets, atmosphere and noise. The atmospheric delay and the noise (assuming a low noise level here within a small area) can be neglected in this case, since the four CRs were placed within 50 m of each other. Therefore, the phase of displacement can be computed as:

$$\Delta {\phi }_{displacement}\cong \Delta \phi -\Delta {\phi }_{height}$$

(3)

and the relationship between displacement in length and phase in radians is:

$$\Delta {\phi }_{displacement}=\frac{4\pi }{\lambda }d$$

(4)

where d is the displacement along the LOS. Hence:

$$d=\frac{\Delta {\phi }_{displacement}}{4\pi }·\lambda \cong \frac{\Delta \phi -\Delta {\phi }_{height}}{4\pi }·\lambda$$

(5)

In this case, the height of the target was provided by the ground survey. For TerraSAR-X, given the wavelength λ equal to 3.1 cm, when we read the phase difference $\Delta \phi$, we can thus calculate the displacement value d by subtracting $\Delta {\phi }_{height}$ from $\Delta \phi$. Meanwhile, the displacement retrieved from the interferometric phase is always ambiguous. In particular, the following equations hold:

$$\phi =\phi \pm 2k\pi (k=1,2,3,...)$$

(6)

$$d=d\pm k\frac{\lambda }{2}(k=1,2,3,...)$$

(7)

For TerraSAR-X, the equation equals:

$$d=d\pm k·15.5mm(k=1,2,3,...)$$

(8)

It is worth noticing that, from Equation (8), the SAR system can only detect a change within a phase cycle between the intervals, 15.5 mm in this case. However, in this project, all of the movement values for the CR test will be smaller than the value of the phase cycle, so we do not have to worry about the phase ambiguity term.

3.5. Direction of Displacement

It has to be emphasized here that, while the lifting/lowering experiments carried out by ground survey are only for the vertical direction, the PSI estimated displacement is sensitive along its LOS direction; thus, it can be affected by both vertical and horizontal components. As shown in Figure 5, the movement of one target is always projected onto the LOS direction in InSAR process, and in turn, the detected velocity along the LOS can be decomposed into an orthogonal coordinated system, where the three directions represent the movement in the east-west, north-south and vertical directions. With only a single geometry, we cannot distinguish between vertical and horizontal movements. This will be an important factor that will later be discussed in the CR validation test.

Figure 5. SAR side-looking geometry and angular parameter.

Figure 5. SAR side-looking geometry and angular parameter.

4. Results

4.1. The Reflectivity Map and Displacement Velocity Map

The reflectivity map is the first processing result that shows how the SAR image looks, and it gives a general idea of the SAR signal backscatter characteristics inside one’s AOI. The reflectivity map is generated by averaging the amplitude of all images within the stack. Figure 6 shows the reflectivity map of one whole frame of the reflectivity map geocoded to Google Earth.

Figure 7 shows the annual velocity of displacement for PS points that have a height ranging from −5 m to 5 m compared to the reference point located at zero elevation. PS points closer to the ground are chosen here because they are more relevant for demonstrating the ground displacement. The overall annual velocity around our AOI is mostly close to zero with partial regions varying from −10 mm/year to 5 mm/year, for example a number of spots inside the west corridor shown in Figure 7 has an approximately −5 mm/year displacement trend.

Figure 6. The reflectivity map for one whole frame of the TSX/TDX dataset.

Figure 6. The reflectivity map for one whole frame of the TSX/TDX dataset.

Figure 7. The annual LOS velocity of the linear displacement map for PS points that have a height within the region of [−5,5] meters located in the Hong Kong urban area.

Figure 7. The annual LOS velocity of the linear displacement map for PS points that have a height within the region of [−5,5] meters located in the Hong Kong urban area.

4.2. Analysis of Ground Displacement

In this section, we try to analyze the cause of displacement in the AOI from the points of geological formations and human activities, which are usually the two dominant factors of ground displacement, and it is common that the two factors mutually promote each other to cause displacement.

To evaluate the geological factors, a geological map geocoded to Google Earth is shown in Figure 8. Most of our AOI is covered in granites, alluvium/colluvium and reclamation land. Granites are known as an intrusive igneous rock, which are generally considered relatively stable as the layer. On the other side, the perimeter of granite land is a large portion of reclamation land, which constitutes the majority of the Kowloon urban area that is flat and concentrated with skyscrapers. Reclamation land is generally not as stable as igneous rocks, and previous studies have shown ground displacement on reclamation land in other parts of Hong Kong [1].

Figure 8. Geological map of the AOI. Granites, alluvium and reclamation land are the three main types of geological formations inside the AOI. Geological data are provided by Hong Kong Geological Survey.

Figure 8. Geological map of the AOI. Granites, alluvium and reclamation land are the three main types of geological formations inside the AOI. Geological data are provided by Hong Kong Geological Survey.

On the other side, human activities might also play an important role. It is well known that the west corridor depicted in Figure 7, standing on reclamation land, has a massive underground railway system, some of it still under construction in our study period. The corner reflector testing site, which shows a displacement from 0 to −12 mm, is also standing at an ongoing construction site where the new West Kowloon Terminus was being constructed. As a matter of fact, there have been reports, including [23], that talked about ground subsidence and building cracks near this construction site and at other places along the west corridor. The ground displacement concentrated within this area is likely the interaction between the geological formation and underground construction projects. For the next step, it is important to conduct more investigation, specifically at those displacements locations that we can see in Figure 7.

4.3. Two Examples of PS Points

As stated above, PSI processing gives two outcomes, the three-dimensional localization of PS points with meter precision and the estimated millimeter displacement of PS along the LOS direction. Here, we present two typical PS points in the area. Figure 9 and Figure 10 show the 3D representations of the PS points’ locations geocoded to Google Earth, where the color of the points indicates the heights of the points, and estimated displacements are shown in time series.

The first PS, displayed in Figure 9, is located at the high elevation of the façade of a high building. The detected height is 116 m; the linear displacement is close to zero (1 mm/year); it is worth noticing that a strong seasonal movement is evident. For estimating the movement of a point on a high building, besides the linear movement model, we use a simple linear thermal expansion model dependent on the height of the point and local temperature. However, the true movement on a high building is generally much more complicated, and even from Figure 9, we can see that the displacement time series is rather noisy apart from the periodic seasonal trend, which means more complicated movements are not fitting well into the linear model. This also explains why we choose points with a height close to the ground for drawing the deformation map as in Figure 7.

The second target, shown in Figure 10, is located at ground level along the coast, close to the actual position where the corner reflectors were installed. The target shows a moving trend at a speed of −5.1 mm/year along the LOS direction. Please notice this displacement, as it will come in handy in the following discussion of the CR validation test.

Figure 9. Displacement time series for a PS point detected on the façade of the building. The point’s elevation is about 116 m. The displacement time series show an evident seasonal moving trend, but the overall average velocity is only 1.0 mm/year, which is close to zero.

Figure 9. Displacement time series for a PS point detected on the façade of the building. The point’s elevation is about 116 m. The displacement time series show an evident seasonal moving trend, but the overall average velocity is only 1.0 mm/year, which is close to zero.

Figure 10. A PS point detected along the coastal area and near the CR validation test test. The target is located at ground level with an estimated elevation of 0.3 m. The displacement time series show a very smooth moving trend with an estimated velocity of −5.1 mm/year along the LOS direction.

Figure 10. A PS point detected along the coastal area and near the CR validation test test. The target is located at ground level with an estimated elevation of 0.3 m. The displacement time series show a very smooth moving trend with an estimated velocity of −5.1 mm/year along the LOS direction.

5. Validation

5.1. Result of the Detected Displacement

The final results of displacement for CRs are reported in Table 2 and Figure 11. For each CR, the time series of displacement measured by the ground survey and calculated by the PSI technique are being compared. The reference CR (denoted as CR4) was kept still all of the time; thus, the displacement time series shows zero. The displacement detected by SAR images was provided by the processing software SARPROZ. By comparing the deviation between the two measurements, we can have the first insight that deviations are mostly within one millimeter. It is worth mentioning that the detected movement may have been caused by either a lifting or a lowering, and the SAR system can only detect a change within a phase cycle between the intervals, 15.5 mm in this case. For this case, we only report the results close to the ground survey. It is true that without the help of the ground truth, it is impossible to draw the lines in a specific way; however, for natural ground targets, the movements always follow a trend, mostly a linear trend or a seasonal trend. If one ground target follows a certain pattern of moving, then it is possible to derive the movements. However, for random moving ground objects, such as vegetation or a water body, it is then impossible for the SAR system to detect the change once the phase change is out of the phase cycle.

Table 2. The comparison of vertical displacements between SAR images and those surveyed.

Date of Satellite Images	Vertical Displacement Detected by SAR Images (mm)				True Vertical Displacement Detected by Survey Team (mm)				Deviation between Survey and SAR Images (mm)
	CR1	CR2	CR3	CR4	CR1	CR2	CR3	CR4	CR1	CR2	CR3	CR4
6 February 2012	0	0	0	0	0	0	0	0	0	0	0	0
17 February 2012	8.4	8.9	4.4	0	10	10	5	0	1.6	1.1	0.6	0
21 March 2012	−5.2	19.6	9.6	0	−5	19	8	0	0.2	−0.6	−1.6	0
1 April 2012	−5.7	18.8	5.6	0	−5	18	5	0	0.7	−0.8	−0.6	0
12 April 2012	−0.3	26.9	−0.2	0	−1	26	−1	0	−0.7	−0.9	−0.8	0
23 April 2012	2.1	19.3	8.5	0	3	20	8	0	0.9	0.7	−0.5	0
4 May 2012	2.0	26.5	17.9	0	1	25	15	0	−1	−1.5	−2.9	0
15 May 2012	2.0	9.3	37.7	0	1	9	35	0	−1	−0.3	−2.7	0
26 May 2012	13.5	−0.6	37.7	0	12	−1	35	0	−1.5	−0.4	−2.7	0
6 June 2012	18.0	−0.8	33.3	0	16	−3	30	0	−2	−2.2	−3.3	0
17 June 2012	21.6	−4.0	30.4	0	18	−5	27	0	−3.6	−1.0	−3.4	0

Table 2. The comparison of vertical displacements between SAR images and those surveyed.

Date of Satellite Images	Vertical Displacement Detected by SAR Images (mm)				True Vertical Displacement Detected by Survey Team (mm)				Deviation between Survey and SAR Images (mm)
	CR1	CR2	CR3	CR4	CR1	CR2	CR3	CR4	CR1	CR2	CR3	CR4
6 February 2012	0	0	0	0	0	0	0	0	0	0	0	0
17 February 2012	8.4	8.9	4.4	0	10	10	5	0	1.6	1.1	0.6	0
21 March 2012	−5.2	19.6	9.6	0	−5	19	8	0	0.2	−0.6	−1.6	0
1 April 2012	−5.7	18.8	5.6	0	−5	18	5	0	0.7	−0.8	−0.6	0
12 April 2012	−0.3	26.9	−0.2	0	−1	26	−1	0	−0.7	−0.9	−0.8	0
23 April 2012	2.1	19.3	8.5	0	3	20	8	0	0.9	0.7	−0.5	0
4 May 2012	2.0	26.5	17.9	0	1	25	15	0	−1	−1.5	−2.9	0
15 May 2012	2.0	9.3	37.7	0	1	9	35	0	−1	−0.3	−2.7	0
26 May 2012	13.5	−0.6	37.7	0	12	−1	35	0	−1.5	−0.4	−2.7	0
6 June 2012	18.0	−0.8	33.3	0	16	−3	30	0	−2	−2.2	−3.3	0
17 June 2012	21.6	−4.0	30.4	0	18	−5	27	0	−3.6	−1.0	−3.4	0

Figure 11. The vertical displacements time series retrieved by InSAR in comparison to the result retrieved by the leveling survey. Note that there are biases between the two values that increase with time. This point will be discussed later.

Figure 11. The vertical displacements time series retrieved by InSAR in comparison to the result retrieved by the leveling survey. Note that there are biases between the two values that increase with time. This point will be discussed later.

5.2. Deviation of the Detected Displacement

Figure 12 shows the linear regression of the displacement between the surveyed value and the PSI-derived value. For each target in the figure, the x-axis value indicates the vertical displacement measured by the ground survey, and the y-axis value indicates the vertical displacement detected by PSI.

Figure 12. The statistical analyses of the deviation between the displacements computed from PSI and by the leveling survey measurement.

Figure 12. The statistical analyses of the deviation between the displacements computed from PSI and by the leveling survey measurement.

The equation of the linear regression is:

$$y=1.06x+0.30$$

(9)

It shows a bias of about 0.3 millimeters and a deviation from the linear relation of about 6%. The linear correlation coefficient, $R^2$, is equal to 0.992, very close to one, which means that the linear correlation between the displacements detected by PSI and by the leveling survey is very high. The root mean square error (RMSE) has a value of 1.20 mm, indicating that the accuracy of PSI is close to one millimeter, confirming the fact that the accuracy of the PSI technique can indeed reach the millimeter level in practical applications.

6. Discussion

6.1. Analysis of the Bias

From Figure 11, we can see that the results retrieved from PSI and the ground survey in May and June were not matching as well as the previous months; the biases have a trend of growing bigger. However, the discrepancy between InSAR and leveling in May and June is not irregular, if we notice that PSI-derived data are always higher than ground survey data. In other words, the result appears shifted. This suggests that some factors occurred somehow and biased the results.

It is important to investigate the reasons for such biases. First of all, we can rule out the possibilities of precipitation and thermal expansion after some calculations. Hence, the most possible explanation is the movement of the ground where the CRs were located. As a matter of fact, we know that the area where the CRs were installed (which were located in the area of Figure 10) was moving.

Here, we try to carry out a small area PSI process to investigate the differential velocity between the four CRs. During the process, we will choose a new reference point close to the test site with a high temporal coherence. Here, the temporal coherence of one PS point is defined as [24]:

$$\xi =max\left|\frac{1}{N}\sum _{m;i=1}^N{e}^{j(\Delta {\varphi }^{i,m}-\Delta {\varphi }_{height}^{i,m}-\Delta {\varphi }_{deform}^{i,m})}\right|$$

(10)

where $\Delta {\varphi }^{i,m}$ is the acquired interferometric phase, $\Delta {\varphi }_{height}^{i,m}$ is the elevation-dependent phase term, ${\varphi }_{deform}^{i,m}$ is the deformation-dependent phase term, N is the total number of scenes for PSI processing and $<i,m>$ means the phase difference between the master image and the i-th image for the selected point. Points with higher temporal coherence can be seen as potential PS candidates. In this case, when choosing a new reference point with a high temporal coherence, we can also observe a moving trend of the points around the CRs, including the reference reflector in May and June. In particular, by taking the four PS points closest to each of the four CRs with a high temporal coherence value, we calculate the annual velocity during the period of the validation test, and we can see from Figure 13 that the reference CR4 has a very different velocity, −7.9 mm/year, compared to the rest, which all have a velocity close to 1 mm/year. The differential movement between the reference CR and the rest explains the biases we observed in Figure 11.

Figure 13. Here, we selected four PS points that were closest to each of the four CRs. The time series analysis of the four PS points shows a moving trend. In this case, since the vertical displacements were biased, this moving trend is very likely to be an eastward horizontal movement.

Figure 13. Here, we selected four PS points that were closest to each of the four CRs. The time series analysis of the four PS points shows a moving trend. In this case, since the vertical displacements were biased, this moving trend is very likely to be an eastward horizontal movement.

Hence, a possible explanation here of the bias could be a horizontal west-east direction movement between the reference CR and the others, recorded by InSAR, but not by ground leveling. We will verify the assumption by using ground GPS records.

6.2. Verification of CR Movement

The coordinates of the four CRs were measured by GPS RTK with the help of a ground survey team at the beginning and end of the experiment. As shown in Figure 5, the SAR system is more sensitive in the west-east direction movement and less sensitive in the north-south direction; hence, here, we only analyze the west-east movement. The west-east coordinates measured by GPS RTK are listed in Table 3. According to Table 3, the reference CR, CR4, was moving towards the other three, while the others remain relatively stable. This moving east is equivalent to a bias that will result in a larger displacement value than it should be. This explains why the SAR results are always higher than those of the ground survey, because the SAR results include this horizontal differential motion, while the ground survey’s measurement does not.

Table 3. CRs’ reference mark coordinates in the east-west direction at the beginning and ending of the validation test.

	CR1 Easting (m)	CR2 Easting (m)	CR3 Easting (m)	CR4 Easting (m)
2 June 2012	834,831.51	834,810.39	834,791.74	834,773.77
11 September 2012	834,831.50	834,810.39	834,791.74	834,773.81
Δ(m)	−0.01	0.00	0.00	0.04

Table 3. CRs’ reference mark coordinates in the east-west direction at the beginning and ending of the validation test.

	CR1 Easting (m)	CR2 Easting (m)	CR3 Easting (m)	CR4 Easting (m)
2 June 2012	834,831.51	834,810.39	834,791.74	834,773.77
11 September 2012	834,831.50	834,810.39	834,791.74	834,773.81
Δ(m)	−0.01	0.00	0.00	0.04

Since the accuracy of GPS RTK is far lower than PSI and only two sets of data were recorded, no quantitative analysis can be carried out. However, we can carry out the quantitative analysis in another way. If we neglect the reference CR in this case and take any other CR as the new reference CR, we can calculate the relative movement between the new reference CR and the rest using exactly the same method shown in Figure 12. Since the other three CRs are relatively stable, the results between SAR measurement and ground survey measurement should be well matched.

For example, if we take CR3 as the new reference, we can get the relative movement of the other two compared to the new reference. The result is illustrated in Figure 14 and Figure 15, from which we can see that the PSI result and optical leveling result are better matching each other. A statistical analysis is also carried out. The outcome of the linear regression is:

$$y=1.03x-0.82$$

(11)

It shows a bias of less than one millimeter and a deviation from the linear relation of about 3%. Comparing with Figure 12, the new one shows a $R^2$ of 0.999 and an RMSE of 0.67 mm, which is a great improvement. This new RMSE is significantly lower than 1 mm, meaning that the accuracy can even be better than 1 mm. This verifies the previous assumption that the default reference CR is not stable during the period of experiment with an eastward horizontal movement. Furthermore, by ruling out this bias, the PSI technique can indeed reach a millimeter to sub-millimeter level accuracy.

Figure 14. The relative vertical displacements time series compared to CR3 retrieved by InSAR in comparison to the result retrieved by the leveling survey.

Figure 14. The relative vertical displacements time series compared to CR3 retrieved by InSAR in comparison to the result retrieved by the leveling survey.

Figure 15. The statistical analyses on the deviation between the displacement detected by InSAR and by the leveling survey measurement when CR3 is taken as the reference CR.

Figure 15. The statistical analyses on the deviation between the displacement detected by InSAR and by the leveling survey measurement when CR3 is taken as the reference CR.

6.3. Strengths and Weaknesses of the PSI Technique

Here, we try to discuss some of the strengths and weaknesses of the PSI technique for monitoring urban areas that came out of this study. As demonstrated in the project, PSI is specialized in wide area monitoring. For example, for TerraSAR-X strip map mode, the image can cover a wide area of approximately 1500 km². Besides, the PSI technique can derive the extremely high density of PSCs from images. Most importantly, PSI is capable of achieving the same level of accuracy of the ground survey, while being much more convenient for monitoring large areas.

However, there are still some limitations. In the first place, the method does not work with vegetated area. As shown in Figure 7, there is no information in the vegetated area.

Secondly, even in non-vegetated areas, PS points are not everywhere. For example, no PS points can be found in the shadowed-by-skyscrapers area, which is often the case in Hong Kong. In order to overcome this problem, it is suggested to use both ascending and descending SAR images.

Thirdly, the continuity of the data relies on the stability of satellite acquisition. If the satellite fails to acquire data for a long period, then no information can be retrieved.

Moreover, if the monitoring area undergoes a subsidence greater than the length of one phase cycle, then PSI is unable to read this fast movement. The solution for this kind of problem in the future is to maintain stable and regular acquisition and to shorten the revisiting time to get more available data.

Last, but not the least, as already shown in Figure 5, PSI is only sensitive along the LOS direction. With one single geometry, we cannot decompose vertical and horizontal movements. One way to solve the problem is again to use both ascending and descending track images.

In conclusion, PSI is a very powerful technique for wide areas and high accuracy monitoring tasks, where it is possible to reach millimeter accuracy in monitoring the displacement of millions of targets all at once. PSI can be extremely useful if used in collaboration with other techniques, so that the advantages can be highlighted and the drawbacks avoided. Future efforts should be spent to find the best way to combine PSI with other techniques for the purpose of monitoring ground displacement more accurately.

7. Conclusions

In this study, we analyzed the ground displacement over the Hong Kong urban area using the PSI technique. Seventy three TSX and TDX images spanning over a four-year time period are used for the data process. The PSI process generated a map of displacement annual velocity along the LOS direction that shows an overall velocity close to zero in most of the regions with partial regions varying from −10 mm/year to 5 mm/year. In particular, a number of spots inside the west corridor shown in Figure 7 and the corner reflector validation test site (Figure 10) are showing some displacement trends that exceed −5 mm/year; both of the sites stand on reclamation land and near ongoing underground/ground construction sites. It is likely that both the geological features of surface land and ongoing human activities resulted in ground displacement in these areas.

To validate if the accuracy of the PSI technique can indeed reach the millimeter level in this study, a corner reflector validation test was carried out. The reflectors were manually moved up or down, and their true values were compared to the value derived from the PSI process. A linear regression of the displacement between the surveyed value and the PSI-derived value was conducted, and the result shows a very high correlation between PSI processing and the ground leveling recorded value, with an RMSE of 1.20 mm, very close to 1 mm. The validation test presents the fact that the PSI technique can reach a millimeter-level accuracy.

As a further effort to the analysis and discarding a trend of bias between PSI and ground leveling that happened in the latter phase of the validation test, we carried out another small area PSI process to investigate the ground movement where the four corner reflectors were located. The analysis shows a horizontal moving trend between the reference reflector and the other three. The fact that the reference reflector was moving horizontally towards the other three explains the bias well, because ground leveling can only record vertical movement, while PSI records both horizontal and vertical movement. After this factor was ruled out in calculating the deviation of displacement between PSI and ground leveling, the new RMSE significantly dropped to 0.67 mm, far lower than 1 mm, showing that the PSI technique can indeed reach up to sub-millimeter level accuracy when no bias is introduced.