## 3 Method
We then set up the cost function related to the observed differences and modeled the differences [26; 27]:  
\[R^{2}=\sum_{i=1}^{N}[\Delta d_{i}-(m_{A}+m_{B}x_{i}+m_{C}y_{i})]^{2} \tag{2}\]  
where \(\Delta d_{i}\) is the observed difference between GPS and InSAR at the GPS station, \(i,x_{i}\), and \(y_{i}\) are the East and North coordinates of the GPS station, \(i\) and \(N\) are the number of GPS stations, and \(m_{A},m_{B},\) and \(m_{C}\) are the coefficients for the two-dimensional trend surface that absorbs the systematic errors in the initial InSAR result. Using the method of least-squares, we estimated \(m_{A},m_{B},\text{and }m_{C}\) by minimizing \(R^{2}\). Finally, we removed the linear trend from the initial InSAR displacement rates. After removing the trend, we used the GPS station LNJS as a reference station, and shifted all the InSAR displacement rates by the offset value between the InSAR velocity and the GPS velocity at LNJS to obtain the calibrated InSAR result (displacement rates) (see Section 4.1). As shown in Table 1, only a few images from the descending track were available and our tests showed that these images contained relatively large noises. In addition, our focus was to see whether a one-year Sentinel-1A result was sufficient to see the annual subsidence in the contemporary leveling results [6] and to experiment with detecting subsidence in wet and dry seasons. As such, we excluded the images from the descending track, and the InSAR result presented in this paper only used images from the ascending track (Table 1). The InSAR result only contained displacement rates along the one-dimensional, LOS direction, which was not able to fully resolve the displacement rates in the three-dimension. While it is possible to resolve the two-horizontal components using SAR images from the ascending track and descending track, this was not carried out in this study due to the low number of images from the descending track and their noises (see the previous paragraph and Table 1), that could degrade our final InSAR result. Such LOS displacement rates have been converted into vertical displacement rates by neglecting horizontal velocities [12]. [PERSON] and [PERSON] [28] have shown that the error in the vertical velocity, when converting LOS to vertical, depends on the incidence angle and the amount of horizontal motion. From a real-world local displacement phenomenon measured with Envisat data from multiple viewing geometries (four ascending and three descending passes), they found that an error of up to 67% of the maximum vertical motion was introduced into the projected vertical component when neglecting the horizontal component. However, the horizontal velocity field was almost of the same magnitude as the vertical in their case. In the Yunlin area, the vertical velocities and the horizontal velocities from April 2016 to April 2017 from GPS observations are shown in Figure 5. The vertical velocities were mostly over 20 mm/a and the horizontal velocities were mostly less than 5 mm/a. The maximum error resulting from the projection of LOS into the vertical velocity induced by the horizontal velocity component was [28]:  
\[\Delta V_{e}=\tan\theta(V_{E}\cos\alpha-V_{N}\sin\alpha) \tag{3}\]  
which could be obtained by dividing the contributions of the \(V_{N}\) and \(V_{E}\) terms in Equation (1) by \(\cos\theta\). Using an average incidence angle of 33.9\({}^{\circ}\) and a heading angle of -12.4\({}^{\circ}\), we found that the average \(\Delta V_{e}\) values at the GPS stations (Figure 5) were -2.4 mm/a, which was the average error due to neglecting the horizontal velocity components. As this error was relatively small, we neglected the horizontal velocities, i.e., \(V_{N}=V_{E}=0\), and converted the calibrated InSAR velocities to vertical velocities as (see Equation (1)):  
\[V_{h}=V_{LOS}/\cos\theta \tag{4}\]  
Figure 5: (**a**) The horizontal velocities and (**b**) the vertical velocities from April 2016 to April 2017 in Yunlin at the 14 GPS stations.