## 2 Methodology
### Least Squares Model with Outlier Detector  
From wrapped phase differences at arcs, the parameters of interest, i.e., the deformation and topographic residuals, can be estimated, where the phase ambiguity is usually resolved by solution space searching [8], 3D phase unwrapping [10], or integer least squares [9]. However, the success rate can never be guaranteed. In TCPInSAR, instead of resolving the ambiguity, we detect and remove it as an outlier. The strategy is mainly based on the fact that when we conduct least squares directly on the observations at arcs, arcs having phase ambiguity in certain observations will have abnormally large residuals. Therefore, we can identify and delete arcs with non-zero phase integers by diagnosing the least squares residuals. Moreover, since we just need to know whether the arc has phase ambiguity or not, rather than the exact location of the ambiguity, the outlier detection procedure can be simplified and implemented very efficiently. The details of the model used in TCPInSAR can be found in [31]. Here we just recall the model briefly. Supposing that \(N\) differential interferograms with temporal baselines of \(t_{i}(i=1,\ 2,\cdots,\ N)\) have been generated from a set of \(M\) SAR images over the study area, for an arc formed by connecting TCP \(x\) and TCP \(y\), observation equations can be written with \(N\) interferometric pairs as follows:  
\[\left\{ \begin{array}{l}\Delta\Phi(x,y,t_{1})=\alpha(x,t_{1})\cdot\Delta h_{x,y}+\beta(t_{1})\cdot\Delta v_{x,y}^{1}+\Delta\omega(x,y,t_{1})\\ \Delta\Phi(x,y,t_{2})=\alpha(x,t_{2})\cdot\Delta h_{x,y}+\beta(t_{2})\cdot \Delta v_{x,y}^{2}+\Delta\omega(x,y,t_{2})\\ \vdots\\ \Delta\Phi(x,y,t_{N})=\alpha(x,t_{N})\cdot\Delta h_{x,y}+\beta(t_{N})\cdot \Delta v_{x,y}^{N}+\Delta\omega(x,y,t_{N})\end{array} \right., \tag{4}\]  
where \(\alpha(x,t_{i})\) is the conversion factor for the differential elevation error (\(\Delta h_{x,y}\)), \(\beta(t_{i})\) is the conversion factor for the relative deformation rate (\(\Delta v_{x,y}^{i}\)), and \(\omega(x,y,t_{i})\) is the differential part of the phase components related to atmospheric delay and other noises. The differential deformation rates of the arc for all the temporal baselines can be treated as constant, i.e., \(\Delta\mathbf{V}_{x,y}\) represents the deformation rate during the whole temporal interval. The Equation (4) can be rewritten in the matrix form as follows:  
\[\Delta\mathbf{\Phi}=\mathbf{A}\left[ \begin{array}{c}\Delta h_{x,y}\\ \Delta\mathbf{V}_{x,y}\end{array} \right]+\mathbf{W}, \tag{5}\]  
Figure 6: TCP network constructed by (**a**) global Delaunay triangulation and (**b**) local Delaunay triangulation. Black dots and gray lines indicate the temporarily coherent points (TCPs) and the arcs, respectively. The figure is plotted in the World Geodetic System (WGS) 84 coordinate system. where:  
\[\mathbf{\Delta\Phi}=\Big{[} \begin{array}{cccc}\Delta\Phi(x,y,t_{1})&\Delta\Phi(x,y,t _{2})&\cdots&\Delta\Phi(x,y,t_{N})\end{array} \Big{]}^{T}, \tag{6}\]  
\[\mathbf{W}=\Big{[} \begin{array}{cccc}\Delta\omega(x,y,t_{1})&\Delta\omega(x,y,t _{2})&\cdots&\Delta\omega(x,y,t_{N})\end{array} \Big{]}^{T}, \tag{7}\]  
\[\mathbf{A}=\Bigg{[} \begin{array}{cccc}\alpha(x,t_{1})&\alpha(x,t_{2})&\cdots& \alpha(x,t_{N})\\ \beta(t_{1})&\beta(t_{2})&\cdots&\beta(t_{N})\end{array} \Bigg{]}^{T}, \tag{8}\]  
The unknown parameters \(\mathbf{\hat{X}}\) can be solved by the LS approach, which can be expressed by:  
\[\mathbf{\hat{X}}=\Bigg{[} \begin{array}{c}\Delta t_{x,y}\\ \Delta\mathbf{V}_{x,y}\end{array} \Bigg{]}=\Big{(}\mathbf{A}^{T}\mathbf{P}^{dd}\mathbf{A} \Big{)}^{-1}\mathbf{A}^{T}\mathbf{P}^{dd}\mathbf{\Delta\Phi}, \tag{9}\]  
where \(\hat{}\) denotes the estimated quantities, and \(\mathbf{P}^{dd}\) is the prior weight matrix, which can be obtained by taking the inverse of the variance matrix \(\mathbf{Q}^{dd}\) of the double-difference phases.