## 3 Methodology  
### SBAS-InSAR Method  
The SARscape Modules for ENVI software suite (HARRIS Geospatial Solutions, Broomfield, CO, USA) were used to perform the Small Baseline Subset InSAR (SBAS-InSAR) processing [32; 33]. \begin{table}
\begin{tabular}{c c} \hline \hline
**Parameters** & **Description** \\ \hline Product type & Sentinel-1A TOPS SLC \\ Band & C \\ Wavelength (cm) & 5.63 \\ Beam mode & IW \\ Orbit direction & Ascending \\ Track No. & 142 \\ Polarization & VV \\ Acquisition time & 20180103-20190416 \\ Azimuth/Range pixel spacing (m) & 2.33/13.99 \\ Incidence angle (\({}^{\circ}\)) & 33.99 \\ Number of data & 36 \\ Revisit frequency (day) & 12 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Main parameters of Sentinel-1 Data. Figure 4: The spatial and temporal baselines of SAR datasets. The main steps were as follows. There are \(N\)+1 SLC (Single Look Complex) images acquired according to time series \(t_{0}\dots\)\(t_{i}\dots\)\(t_{N}\). The first step of the SBAS method was \"connection graph\", which selected SAR image pairs with spatiotemporal baseline constraints. \(M\) interferograms were generated after the removal of the topographic phase based on orbital information and external DEM data in the step of \"differential interferogram generation\". Multi-looking in range and azimuth direction and the Goldstein filtering method was used to improve the signal/noise ratio of the interferograms [34]. Then followed the completion of phase unwrapping of the wrapped phase signal calculated from the selected interferograms by Minimum Cost Flow (MCF) method [35; 36; 37]. Several interferogram pairs with unwrapped phase error or low coherence were be discarded [38]. The next step was \"refinement and re-flattening\". By selecting and refining the stable Ground Control Point (GCP), the residual phase content and phase ramps were calculated to correct the unwrapped phase [39]. Then came the key step of \"first estimation of residual height and the displacement velocity\". In this step, assuming that the time \(t_{0}\) is used as the reference time, and the subsidence of the study area is regarded as zero, the phase of a pixel in the \(j\)-\(th\) interference image relative to the unwrapping start point can be expressed as follows:  
\[\delta\varphi_{j}=\varphi_{t_{1}}-\varphi_{t_{2}}\approx\delta\varphi_{j}^{def} +\delta\varphi_{j}^{topo}+\delta\varphi_{j}^{atm}+\delta\varphi_{j}^{noise}, \tag{1}\]  
where, \(1\leq j\leq M\), \(\delta\varphi_{j}^{def}\) is subsidence phase, \(\delta\varphi_{j}^{topo}\) represents topography phase, \(\delta\varphi_{j}^{atm}\) represents atmosphere phase, and \(\delta\varphi_{j}^{noise}\) represents noise phase. Those could be further expressed as follows:  
\[\left\{ \begin{array}{ll}\delta\varphi_{j}^{disp}(x,r)=\frac{4\pi}{\lambda} [d(t_{2})-d(t_{1})],&\forall j=1,\dots,M\\ \delta\varphi_{j}^{topo}(x,r)=\frac{4\pi}{\lambda}\cdot\frac{B_{\perp}\Delta h }{r\sin\theta}\\ \delta\varphi_{j}^{atm}(x,r)=\varphi_{atm}(t_{2})-\varphi_{atm}(t_{1})\end{array}  \right., \tag{2}\]  
where, \(\lambda\) is the radar wavelength, \(\theta\) represents the radar angle, \(\Delta h\) accounts for the error of DEM, \(r\) accounts for the slant range. \(d(t_{1})\) and \(d(t_{2})\) represent the cumulative amount of subsidence compared to the initial time (\(t_{0}\)) in the direction of the radar sight. We generally assume that \(d(t_{0})=0\). In the SBAS method, it is assumed that the surface subsidence in the interval between the acquisition of two adjacent SAR images conforms to the linear cumulative trend [39], and the solution of the phase time series is converted into the solution of the phase change rate.