ESTIMATING 3D LAND SUBSIDENCE FROM MULTI-TEMPORAL SAR IMAGES AND GNSS DATA USING WEIGHTED LEAST SQUARES

Analysis of multi-temporal synthetic aperture radar (SAR) satellite images using persistent scatterer interferometry is an effective approach for monitoring land subsidence, which is a serious issue in some urban areas. However, a drawback to this approach is that it is limited to displacement along the radar line-of-sight direction. An accurate understanding of land subsidence requires estimation of 3D displacement. One solution is to combine observations from multiple sources and directions, such as multi-temporal SAR images acquired on ascending and descending orbits, with global navigation satellite system (GNSS) data. While this approach estimates 3D displacement, other methods do not account for differences in data accuracy. Therefore, in this paper, we propose a method for estimating 3D land subsidence from multi-temporal SAR images and GNSS data by using the weighted least squares method. The weights for data sources are calculated from the PSI results and GNSS data. We apply the method to Kansai International Airport, using 13 ALOS-2/PALSAR-2 ascending images from 2014 to 2018 and 17 ALOS-2/PALSAR-2 descending images from 2015 to 2018. Root mean squared errors in the east–west, north–south and vertical directions are 6, 13, and 10 mm/year, respectively. These results demonstrate that combining PSI and geodetic results is effective for monitoring land deformation accurately with high spatial resolution.


INTRODUCTION
Monitoring environmental changes in urban areas is essential for local governments to maintain quality of life. Land subsidence can hinder urban management and may occur as a result of excessive extraction of groundwater or natural gas. Land subsidence can then lead to flooding and damage to infrastructure. Thus, it is important to detect subsidence early and take appropriate countermeasures. Leveling surveys using a global navigation satellite system (GNSS) such as the Global Positioning System (GPS) are a conventional approach for such monitoring. These point-based measurement approaches can measure subsidence with high accuracy, but not practical for wide-area monitoring. In contrast, satellite synthetic aperture radar (SAR) is an effective tool for a wide-area monitoring. Differential interferometric SAR (DInSAR) is a technique for observing displacement with high resolution; in particular, permanent scatterer interferometry (PSI) (Ferretti et al., 2000) estimates displacement with high accuracy using dozens of SAR images. However, with this technique, displacement can be estimated in only the direction along the radar's line of sight (LOS) and it is impossible to distinguish between the horizontal and vertical displacement. Therefore, it is necessary to separate the displacement in radar LOS direction into 3D components. Ito et al. (2019) proposed a method for estimating 3D displacement from SAR images. First, the method estimates the radar LOS direction velocity using PSI. Next, it estimates the 3D direction velocity from GPS and leveling survey data by interpolation (ordinary kriging). Finally, assuming equal weights for each observation, the ordinary least squares method (OLS) is used to combine the radar LOS direction velocity and the 3D direction velocity to estimate the 3D direction velocity with high accuracy and high resolution. However, using OLS to combine different data sources with equal weights may lower the accuracy of the displacement results because the final estimate will be strongly affected by any lower-accuracy data that are present. The weighted least squares method (WLS) can resolve this issue, but it is necessary to determine the optimal weights to apply WLS. In this paper, we propose a method for determining the weights to be used in WLS and estimate 3D displacement by combining PSI velocity and interpolated velocity by WLS.

Concept of the Proposed Method
A flowchart of the proposed method is shown in Figure 1. The proposed method is assumed to use multi-temporal SAR images and GNSS data. After estimating the displacement velocities by PSI, we apply OLS to both velocities from the PSI and GNSS data. The estimated 3D displacement field is regarded as an initial estimate. The weights required for WLS are derived from this initial estimate and the GNSS data. Once we derive the weights, we apply WLS to the velocities with the weights. Finally, we estimate the final 3D displacement field. We can calculate the displacement value by integrating the velocities from a designated reference to the permanent scatterer (PS) of interest.
. Flowchart of the proposed method for estimating 3D displacement.

Persistent Scatterer Interferometry
PSI involves selecting PSs as phase-stable pixels and using them alone 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 at 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 displacement is estimated relative to a certain point known as the zerodisplacement point. Thus, we convert the relative values to absolute values by adjusting the average PSI to those of the reference points to connect PSI to a reference frame. Specifically, we calculate the velocity difference between the leveling points and the nearest PS point and add the average of the difference to the PSI results.

Interpolation of Geodetic Data
We assume that temporal GNSS data are available. We generate the velocities along three directions, namely, the east-west, north-south, and vertical axes, in accordance with the velocities estimated from SAR images by PSI. Because the point data are spatially sparse, a majority of PSs have no neighboring GNSS data around them. Therefore, the GNSS velocities must be interpolated to match the spatial resolution of the PSI results. Several approaches can be used for interpolation. For example, inverse distance weighting (IDW) can interpolate by calculating the weight that is inversely proportional to the distance between available point data. In this research, we use ordinary kriging (Brooker, 1991;Deutsch et al., 1998), an interpolation method that allows for easy processing and intuitive understanding.

Estimation of Three-dimensional Displacement
The initial 3D displacement field is estimated from the PSI and interpolated results by applying OLS. The observation equations of the OLS model are shown as follows: Here, VASC and VDES are the estimated velocities of a PS along the radar LOS from images acquired on ascending and descending orbits, respectively.
is the heading angle of the satellite, and θ is the incident angle of the radar. ue, un, and uu are the unit vectors pointing from the PS toward the satellite. These three are the unknown components of the velocity vector.
Equation (1) can be rewritten as follows: Now we define the residual r for Equation (2): Applying OLS to the residual minimizes the sum of the squared residuals. The optimal estimate of vdis is given by Whereas Equation (4) assumes equal weights, WLS assumes non-uniform weights, expressed as P in the following equation: P is the following diagonal matrix: ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-3-2020, 2020 XXIV ISPRS Congress (2020 edition) Here, σ 2 represents the variance of the velocities estimated from PSI or obtained from GNSS data.

Determination of Weights
We now explain how to determine the weights in Equation (6). We assume that errors are included in the displacement velocity estimated by PSI and from the GNSS data, and that the former errors are larger than the latter errors. In addition, we also assume that no data other than temporal SAR images and GNSS data are available. Under these assumptions, it is impossible to determine the weights for PS and GNSS displacement rates independently. Instead, our approach is to utilize GNSS data to determine the weights of PSs. For these weights, we start matching the GNSS point with the nearest PS, as shown in Figure 2. The distance between the two points is measured as the Euclidean distance.
Then, we calculate the variance of the PS velocities as follows: Here, n denotes the number of GNSS data points.
Next, we explain the weights for GNSS data. Figure 3 represents an example of taking out GNSS data. Figure 3. Taking out GNSS data to calculate the weight of the GNSS data.
The velocity of the GNSS data of interest is compared with the interpolated velocity estimated from the surrounding GNSS data by interpolation. The difference is defined as = − As with the variance of PS velocities, we calculate the variance of the GPS velocities :

STUDY AREA AND DATA USED
We selected Kansai International Airport (KIX; Figure 4) in western Japan as the study area for this research. It has been experiencing serious land subsidence (Kansai Airports., n.d.). KIX is located 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.
We used level 1.1 (L1.1) satellite SAR images from the Phased Array type L-band SAR (PALSAR-2) instrument onboard the Advanced Land Observing Satellite 2 (ALOS-2). Table 1 lists the images used in this research: 13 acquired on the ascending orbit and 17 acquired on the descending orbit.
For the PSI analysis, we downloaded the 10-m-grid digital elevation model published by the Geospatial Information Authority of Japan (GSI, n.d.) to remove the effect of topography. We selected Island II because it has exhibited larger 3D displacement than Island I. Figure 4b shows the 54 leveling and 26 GPS locations measured at Kansai Airport in an annual survey, for a total of 80 GPS data points. The leveling data measured using GPS has only vertical displacement whereas the GPS data has 3D displacement. Because the leveling was conducted using GPS, we assumed that the accuracy of the leveling data was almost the same as that of the GPS data available in this research. Below, we incorporate the leveling data into the GPS data.

RESULTS
We divided the 80 GPS data points into two groups, one for interpolation (50 data points) and the other for validation (30 data points). Figure 5 shows the locations of geodetic data used for interpolation and validation. We used 32 leveling and 18 GPS locations (Figure 5a) for interpolation and the other 22 leveling locations and 8 GPS locations (Figure 5b) for validation of the land-subsidence results. We used ordinary kriging to interpolate the data. Figure 6 shows the semivariogram curves of the interpolation from GPS data. Curves were obtained for the east-west, north-south, and vertical displacement components, and the interpolated results from GPS data were generated.  Following the method described in Section 2, the variance of the velocities derived from PSI and obtained from the GPS data was calculated (Figure 8). Figure 9a-c shows the velocity field in all three directions in units of mm/year, estimated by combining the results from PSI with the terrestrial measurements at KIX. We compared the final vertical deformation velocity of the nearest PS point and the GPS data classified as validation data to assess the validity of the estimated results. Figure 9d-f shows the combined estimation versus the validation data. The validation was conducted by matching PS and GPS locations within a radius of 100 m (Euclidean distance). Whereas eight GPS locations were available for validating the east-west and north-south displacement, only seven data points were used for the actual validation of the east-west displacement (Figure 8d) because one location had no PS point within a radius of 100 m. The root mean squared errors (RMSEs) of the combined results were 6, 13, and 10 mm/year for the east-west, north-south, and vertical components, respectively.

DISCUSSION AND CONCLUSION
In this paper, we estimated 3D land-subsidence components at KIX by integrating the results estimated from SAR images observed on ascending and descending orbits with GPS data. We used WLS to examine the difference in accuracy between the SAR-derived estimates and the GPS data. Figure 10 shows the validation results of the vertical components obtained by OLS. The RMSEs of the results using OLS were 14, 16, and 14 mm/year for the east-west, north-south, and vertical components, respectively. The RMSEs of the results using WLS were 6, 13, and 10 mm/year for the east-west, northsouth, and vertical components, respectively. The accuracies of the vertical results are summarized in Table 2.
For reference, the results from the PSI ascending and descending images are relatively worse than those usually obtained using PSI. This was due mainly to the mean displacement velocity in the study area being relatively large at 34 cm/year, which accordingly resulted in a larger obtained RMSE. In addition, the result of 2D fusion of PSI ascending and descending results was 20 mm/year. This estimate was performed by assuming no north-south displacement and applying OLS. It was confirmed that combining the results from PSI and geodetic deformation measurements using WLS is more effective for land-subsidence monitoring with high spatial resolution compared with using PSI, interpolation only, or combining the results using OLS.  It has been reported that Islands I and II 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 not the seaward side (Furudoi et al., 2009 Next, we discuss the determination of the weights used in Equation (6). We resolved the observation equations using WLS and obtained the 3D displacement velocities. Several methods have been proposed for obtaining 3D displacements, such as combining DInSAR and GPS using WLS (Hu et al., 2012;Samsonov et al., 2007;Shi et al., 2015). The method proposed in this paper determines the weights by comparing the velocities of PSI with GNSS data. As shown in Figure 8, the variance of the vertical component was much larger than that of the others. This is a reasonable finding because the vertical displacement is much larger than the horizontal displacement. Compared with the results obtained using OLS, the proposed method using WLS helped to improve the accuracy of the 3D displacement.

Method
Finally, we discuss the transferability of the proposed method in terms of the number of GPS data points. In the study area, we had 54 leveling and 26 GPS locations. As previously mentioned, because the leveling data were measured using GPS, we had a total of 80 GPS data points. In actual application, it is not feasible to have this many GPS locations in such a relatively small area. As a result, it may be difficult to calculate the weight matrix P given by Equation (6). Therefore, a possible approach is to apply the weights obtained in this research to the area of interest. In the future, we aim to apply this method to areas with sparser GPS points and examine whether such an approach is effective.