## 3 Methodology
### Time Series Deformation Processing  
After the phase unwrapping, we selected highly coherent targets (coherence \(>\) 0.7) to extract the time series deformation. Let us assume that there are N + 1 scene SAR images coregistered to the master image, and the time order is \((t_{0}\dots t_{N})\) which can generate M interferograms. The original phase of each SAR image is \(\phi=\left[0,\phi^{1}\dots\phi^{N}\right]\), where we assume \(\phi^{0}=0\). \(\phi\) contains the ground deformation phase, atmospheric phase, topographic phase and noise phase. For each pixel, the model is described as:  
\[\Delta\phi=\Delta\phi \tag{4}\]  
where \(\Delta\phi\) represents the interferometric phase and A is the M \(\times\) N \(-\) 1 coefficient matrix. In matrix A, 1 represents the master image, and \(-\)1 represents the slave image. Then, the coherence weighted least square method is used to obtain raw phase \(\phi\)[58; 59]:  
\[\hat{\phi}=argmin\ \Big{|}\Big{|}W^{\frac{1}{2}}(\Delta\phi-A\phi)\Big{|} \Big{|}_{2}=\Big{(}A^{T}WA\Big{)}^{-1}A^{T}WA\phi \tag{5}\]  
The estimated original phase \(\hat{\phi}\) includes the following items:  
\[\hat{\phi}=\hat{\phi}_{dis}+\hat{\phi}_{trope}+\hat{\phi}_{geom}+\hat{\phi}_ {resid} \tag{6}\]where \(\hat{\phi}_{dis}\) represents the ground deformation phase, \(\hat{\phi}_{tropo}\) represents the atmospheric phase, \(\hat{\phi}_{geom}\) represents the topographic phase, and \(\hat{\phi}_{resid}\) represents the residual phase. W is a M \(\times\) M diagonal weight matrix:  
\[\text{W}=\text{diag}\Bigg{\{}\frac{2L\gamma^{1^{2}}}{1-\gamma^{1^{2}}},\dots, \frac{2L\gamma^{M^{2}}}{1-\gamma^{M^{2}}}\Bigg{\}} \tag{7}\]  
\[\gamma^{\prime^{2}}=\frac{1}{\sigma_{\Delta\Phi^{j}}^{2}} \tag{8}\]  
where \(\sigma_{\Delta\Phi^{j}}^{2}\) is the phase variance of the \(j_{\text{th}}\) interferogram calculated through the integration of the phase probability distribution function and L is the number of looks. After the raw phase is calculated by the least squares method, the displacement time series can be obtained:  
\[\hat{\phi}_{dis}=\hat{\phi}-\hat{\phi}_{tropo}-\hat{\phi}_{geom}-\hat{\phi}_{ resid} \tag{9}\]  
For atmospheric phase removal, global atmospheric models (GAMs), particularly ERA5, have great potential in atmospheric phase screen (APS) correction for InSAR applications on the Tibetan plateau [60]. The ERA5 model has been successfully applied in the APS correction in the QTP [41; 47; 54]. Therefore, in this study, we used ERA-5 reanalysis data to estimate the tropospheric delay [61]. For the topographic residual correction, the topographic phase residual caused by a DEM error was estimated based on the proportionality with the perpendicular baseline time series [62]:  
\[\hat{\phi}^{i}-\hat{\phi}^{i}_{tropo}=\frac{-4\pi}{\lambda}\Bigg{(}\frac{B^{i }_{\perp}}{rsin(\theta)}z_{\epsilon}+\sum_{k=0}^{2}\frac{c_{k}(t_{i}-t_{1})^{k }}{k!}+\sum_{l\in I_{s}}s_{l}H(t_{i}-t_{1})\Bigg{)}+\hat{\phi}^{i}_{resid} \tag{10}\]  
where \(B^{i}_{\perp}\) is the perpendicular baseline, \(\theta\) is the incidence angle, \(H(t_{i}-t_{1})\) is a Heaviside step function and \(z_{\epsilon}\), \(c_{k}\) and \(s_{l}\) are the unknown parameters, which can be estimated by minimizing the \(L^{2}\)-norm of the residual phase. The time series deformation processing was implemented by MintPy ([[https://github.com/insarlab/MintPy](https://github.com/insarlab/MintPy)]([https://github.com/insarlab/MintPy](https://github.com/insarlab/MintPy)), accessed on 20 March 2022) [63].