INTEGRATING MULTI-TEMPORAL SAR IMAGES AND GPS DATA TO MONITOR THREE-DIMENSIONAL LAND SUBSIDENCE

Synthetic aperture radar (SAR) is an effective means of monitoring land subsidence, and differential interferometric SAR (DInSAR) is commonly used to acquire the necessary data. In particular, persistent scatterer interferometry (PSI) can be used to measure land subsidence accurately over a wide area from multi-temporal SAR images. However, the estimated displacement is obtained only in the radar line-of-sight (LOS) direction, making it necessary to develop a method for measuring three-dimensional displacements by combining multidirectional observations. Therefore, we propose herein a method for estimating three-dimensional displacement velocities by combining the results from PSI and geodetic deformation measurements, namely, Global Positioning System and leveling data. We apply the least-squares method to Kansai International Airport in Japan by using 13 ALOS-2/PALSAR-2 ascending images from 2014 to 2018 and 17 ALOS-2/PALSAR-2 descending images from 2015 to 2018. In validation, the rootmean-square errors are 14, 16, and 14 mm/year for the east–west, north–south, and vertical components, respectively, showing that combining PSI results and geodetic deformation measurements is effective for monitoring land subsidence.


INTRODUCTION
Land subsidence is now a serious problem globally because of climate change and ground-water extraction, making inexpensive techniques for monitoring land subsidence highly necessary.In particular, we are interested in monitoring land subsidence in and around civil infrastructure such as airports built on landfill.Synthetic aperture radar (SAR) is an effective means of doing so, and differential interferometric SAR (DInSAR) can be used to measure land subsidence over a wide area.In particular, persistent scatterer interferometry (PSI) (Ferretti et al., 2000) extracts and exploits stable scatterers, known as persistent or permanent scatterers (PSs), to generate accurate land-subsidence maps by excluding erroneous components and extracting the actual land-subsidence components from multi-temporal SAR images.However, the estimated displacement is obtained only in the radar line-ofsight (LOS) direction.We can convert that one-dimensional displacement into vertical displacement by assuming that no horizontal displacement occurs and multiplying the factor calculated using the incidence angle.However, that assumption may not be applicable to civil infrastructure that is prone to three-dimensional (3D) displacement.In addition, earthquake engineers and researchers require 3D displacement maps to understand better the actual motion caused by earthquakes.Therefore, a method is needed for measuring 3D displacements by combining multidirectional observations.A simple way to realize 3D displacement monitoring is to synthesize the one-dimensional displacement vectors along the LOS direction.Fujiwara et al. (2000) reported a two-and-a-halfdimensional (2.5D) surface-deformation analysis method that merges the deformation vectors obtained from ascending and descending SAR images and decomposes them into horizontal and vertical deformation components.Similar approaches have been used to analyze the vertical and horizontal displacements induced by earthquakes.However, SAR observations are insensitive to horizontal displacements in the north-south direction because the radar signals are emitted and received by on-board satellite sensors in the east-west direction, meaning that the horizontal displacement estimated by 2.5D analysis is equivalent to that in the east-west direction.Manzo et al. (2006) reported a method for estimating such east-west and vertical displacements with the additional constraint that no north-south displacement occurs.Consequently, combining satellite images has the limitation that it is difficult to estimate the actual 3D displacement.To overcome this north-south limitation, multiple-aperture interferometry (MAI) (Bechor and Zebker, 2006) has been developed.This method allows us to measure the along-track displacement by generating forwardand backward-looking interferograms and creating a multipleaperture interferogram from them.By using MAI, we can obtain an azimuth displacement that is almost in the north-south direction.Jung et al. (2011) reported a method for estimating the 3D displacement due to the 2007 eruption of the Kilauea volcano in Hawaii by combining DInSAR and MAI measurements from ascending and descending SAR images.However, MAI is less precise than standard InSAR because of phase decorrelation.Consequently, no studies have used MAI to monitor civil infrastructure because MAI is suitable for monitoring ground displacements that are larger than those normally associated with such infrastructure.
Meanwhile, because global satellite navigation systems, including the Global Positioning System (GPS), provide 3D displacement with sub-centimeter accuracy, point-based and spatially sparse observations can be merged with estimates from satellite images.For example, Shi et al. (2015) integrated GPS and DInSAR data to estimate the 3D deformation of the Kilauea volcano.However, such previous studies dealt with displacement on a large scale (e.g., dozens of square kilometers) and did not address whether such approaches are effective for monitoring land subsidence in and around civil infrastructure.Therefore, we propose herein a method for estimating 3D displacement velocities by combining PSI data with external geodetic displacement measurements.We report results of doing so, and we discuss the utility and difficulty of the method when applied to the monitoring of land subsidence in and around civil infrastructure.

STUDY AREA AND DATA USED
For the study area, we selected Kansai International Airport (KIX) in western Japan, which is reported to have been suffering from land subsidence (Kansai Airports., n.d.).KIX sits on two islands, Islands I and II, and in 2015 the subsidence velocities of those two islands were 6 cm/year and 34 cm/year, respectively.In this research, we used level 1.1 (L1.1) images from the Phased Array type L-band SAR (PALSAR-2) instrument onboard the Advanced Land Observing Satellite 2 (ALOS-2).As listed in Table 1, we had access to 13 images acquired on the ascending orbit from September 2014 to March 2018 and 17 images acquired on the descending orbit from March 2015 to May 2018 were available.For the PSI analysis, we downloaded the 10-m-grid digital elevation model published by the Geospatial Information Authority of Japan (GSI) (GSI, n.d.) to remove the effect of topography.As the area of interest within the study area, we selected only Island II because there are no GPS data for Island I. We used 54 leveling and 26 GPS stations measured by Kansai Airports in an annual survey (shown in Fig. 2), and in Fig. 3 we present a box plot showing the difference in vertical velocity between the GPS and leveling data.
Table 1.Observation dates of PALSAR-2 images used for the experiment (yyyy/mm/dd)

Persistent scatterer Interferometry
PSI involves selecting PSs as phase-stable pixels and using only them to reduce the main errors (e.g., temporal and geometrical decorrelation, atmospheric artifacts) in conventional processing methods.In the present study, we used PSI software developed by the infrastructure monitoring group in the Earth Observation Research Center (EORC) of the Japan Aerospace Exploration Agency (JAXA) (JAXA/EORC., n.d.).Once the user has selected the region of interest, the software processes the PSI data from the raw data.The estimated displacement is a displacement relative to a certain point known as a zerodisplacement point.Thus, we convert the relative values to absolute values by adjusting the average PSI to that of the reference points to connect PSI to a reference frame.To be more specific, we calculate the velocity difference between the leveling points and the nearest PS point and add the average of the difference to the PSI result.
The velocity estimated by PSI, VPS, is that in the radar LOS direction and is related to the radar LOS unit vector as follows: where VPS is the estimated velocity of a PS along the radar LOS, are the unknown components of the velocity vector, and ue, un, and uu are the unit vectors pointing from the PS toward the satellite.

Interpolation of Geodetic Data
We assume that spatially sparse geodetic data are available.In this research, GPS data representing 3D displacement and leveling data representing vertical displacement were available.The leveling and GPS velocities must be interpolated to match the spatial resolution of the PSI results.Herein, we use the ordinary kriging (Brooker, 1991;Deutsch et al., 1998) interpolation method.The interpolated velocity V is related to the velocity components by  = (2) where are the east, north, and vertical interpolated velocity components, respectively.Note that the vertical velocity component was generated from both the GPS and leveling data, and therefore the suffix of the variable is expressed as "GPS + leveling."

Estimation of Three-dimensional Displacement
The final 3D displacement field is estimated from the PSI and interpolated results by applying the least-squares (LS) method.The observation equations of the LS model are (3) where is the vector of displacement velocities, are the PSI velocities estimated along the ascen d i n g an d d escending orbits, respectively, and and are the unit vectors along the ascending and descending orbits, respectively.We resolve the observation equations using the LS approach and obtain the 3D displacement velocities.

RESULTS
Figure 3 shows that the deformation characteristics differ between the GPS and leveling data, which is due to locational differences.The GPS positions are mostly at the embankment, whereas the leveling positions are on the reclaimed land.Thus, the GPS velocities corresponding to the embankment are lower than the leveling velocities corresponding to the reclaimed land.As such, it is necessary to use both the GPS data and the leveling data for vertical interpolation.In the present study, we used 32 Figure 6(a)-(c) show the results of the ordinary kriging in units of mm/year.The horizontal results show that Island II is contracting toward its center, excluding the northwest part.We assessed the validity of the vertical interpolated results using the same method, and the RMSE was 32 mm/year.Figure 7(a)-(c) show the velocity field in all three directions in units of mm/year, estimated by combining the results from PSI and terrestrial measurements at KIX.The north-south results are the same as the interpolated results because the SAR measurements were less sensitive to the north-south displacement because of the orientation of the SAR track.The vertical results show that the estimated settlement at the embankment is slower than that at the island interior.We also assessed the validity of the estimated result using the same method as above.

DISCUSSION AND CONCLUSION
In this paper, we estimated the 3D land-subsidence components at KIX by combining GPS data, leveling data, and PSIestimated results from SAR images observed on ascending and descending orbits, and we compared the combined results with validation data.The RMSEs were 14, 16, and 14 mm/year for the east-west, north-south, and vertical components, respectively, showing the reliability of the results for the landsurface estimation.The accuracies of all the vertical results are summarized in Table 2. From these, it seems that combining the results from PSI and geodetic deformation measurements is more effective for land-subsidence monitoring with high spatial resolution than is using PSI or interpolation only.Regarding the horizontal displacements, it has been reported that Islands I and II of KIX are contracting toward their respective centers because of unequal settlement of the seawall, which means that there is a load on the reclaimed side but none on the seaward side (Furudoi et al., 2009).The horizontal displacement measured herein, excluding the north-south displacement in the northwest part, is consistent with the report by Furodoi et al. because it shows that Island II is contracting toward its center.As mentioned above, the SAR measurements did not pertain to the north-south displacement because of the orientation of the SAR track.Thus, the north-south component of the combined results in the northwest part is inconsistent with the report by Furodoi et al. because of incorrect interpolation there due mainly to there being no reference GPS points for interpolation.
We examined the effect of the number of points used for interpolation on land displacement accuracy.We changed the numbers of both GPS and leveling points, and implemented the interpolation.The final results obtained were validated and the RMSEs were calculated.Figure 8 shows that even when 4 GPS points and 4 leveling points were used, the RMSE of vertical displacement was almost the same as that in the case that 18 GPS points and 4 leveling points were used.Thus, the number of GPS points used for interpolation is not critical to the accuracy of the vertical displacement.On the other hand, the RMSE in the north-south direction is sensitive to the number of GPS points.This may be because the SAR observations are insensitive to the north-south displacement due to the observation geometry, and the estimation of the north-south displacement strongly depends on the interpolation results of GPS data.
We now discuss the interpolation technique.The reason why the RMSE of the interpolated results is larger than those of the PSI results is probably due mainly to the difference in the deformation characteristics between GPS and leveling.As described above, most of the GPS points were located at the embankment whereas the leveling measured the surface deformation.Thus, where the GPS measurements dominated the interpolation, the interpolated velocity was different from the leveling one.
Finally, we turn our attention to the optimization technique.We resolved the observation equations using the LS approach and obtained the 3D displacement velocities.However, with this approach, the variances of the GPS, leveling, and PSI data are the same.Recently, several methods have been proposed to obtain 3D displacements by combining DInSAR and GPS using the weighted LS approach (Hu et al., 2012;Samsonov et al., 2007).Thus, in future work we will use the weighted LS approach to obtain 3D displacement velocities.

Figure 2 .
Figure 2. Study area: Kansai International Airport (KIX), Osaka, Japan.(a) The white rectangle represents the location of KIX in Japan.(b) Islands I (lower) and II (upper) of KIX.In (b), the red and green points represent leveling and GPS locations, respectively.The images are from Google Earth.Map data ©2018 Google, ZENRIN

Figure 5 .Figure 6 .
Figure 5. PSI results: (a) mean velocity field from ascending data in mm/year, (b) validation results from ascending data, (c) mean velocity field from descending data, and (d) validation results from descending data.Map data ©2018 Google, ZENRIN Figure6(a)-(c) show the results of the ordinary kriging in units of mm/year.The horizontal results show that Island II is contracting toward its center, excluding the northwest part.We assessed the validity of the vertical interpolated results using the same method, and the RMSE was 32 mm/year.Figure7(a)-(c) show the velocity field in all three directions in units of mm/year, estimated by combining the results from PSI and terrestrial measurements at KIX.The north-south results are the same as the interpolated results because the SAR measurements were less sensitive to the north-south displacement because of the orientation of the SAR track.The vertical results show that the estimated settlement at the embankment is slower than that at the island interior.We also assessed the validity of the estimated result using the same method as above.Figure7(d)-(f) show the combined estimation versus the validation data.The RMSEs of the combined results are 14, 16, and 14 mm/year for the east-west, north-south, and vertical components, respectively.Method RMSE PSI ascending 22 mm/year PSI descending 18 mm/year Interpolated GPS+Leveling 32 mm/year Fusion of GPS, leveling, and PSI 14 mm/year

Figure 7 .
Figure 7. Effect of of the number of GPS and leveling points used for interpolation on land displacement accuracy.

Table 2 .
Comparison between estimated vertical results and validation data