## 2 Data Processing and Analysis Method
### OTL Model and Tidal Signal Analysis
The phasor vector relationship between the GPS estimate and OTL model is shown in Figure 2.  
The phasor difference between the GPS estimate and OTL model for the \(j\)-th tidal constituent at station \(k\) is as follows:  
\[Z_{k,j}=Z_{GPS,k,j}-Z_{OTL,k,j} \tag{3}\]  
where \(Z_{GPS,k,j}\) is the tidal displacement estimated by GPS, and \(Z_{OTL,k,j}\) denotes the tidal deformation predicted by the OTL model. The complex form of \(Z_{k,j}\) can be written as:  
\[Z_{k,j}=A_{k,j}(\cos\Phi_{k,j}+i\sin\Phi_{k,j}) \tag{4}\]  
where \(i\) is the imaginary unit; and \(A_{k,j}\) and \(\Phi_{k,j}\) represent the residual amplitude and phase of the \(j\)-th tidal constituent, respectively. The RMS of the phasor difference for the \(j\)-th tidal constituent between the GPS estimate and OTL model for \(N\) stations is calculated as follows:  
\[\text{RMS}_{j}=\left(\frac{1}{N}{\sum_{k=1}^{N}|Z_{k,j}|^{2}}\right)^{\frac{1} {2}} \tag{5}\]  
In a small region, the phasor difference for the \(j\)-th tidal constituent between the GPS estimate and OTL model at station \(k\) can be divided into two kinds of components. The first kind of component is related to the specific station, i.e., the multipath effect on GPS observation. The second kind of component is applicable to all stations, i.e., the unmodeled body tide, the uncertainties of the OTL model caused by ocean tide model errors and mantle anelasticity, the GPS orbit errors, and the signal propagation delay. Therefore, the phasor difference \(Z_{k,j}\) of the \(j\)-th tidal constituent at station \(k\) can be decomposed into \(Z_{residual,k,j}\), which is related to the specific station, and \(Z_{system,j}\), which is applicable  
Figure 2: Phasor vector diagram of the relationship between the GPS estimate and OTL model. to all stations. For the \(N\) stations, the phasor difference of the \(j\)-th tidal constituent can be expressed as follows:  
\[Z_{k,j}=Z_{system,j}+Z_{residual,k,j}\ \ \ (k=1,2,\cdots,N) \tag{6}\]  
Next, we explain how to obtain \(Z_{residual,k,j}\) and \(Z_{system,j}\) on the complex plane by the least-squares method. Equation (6) is written in matrix form:  
\[\left[ \begin{array}{c}Z_{1,j}\\ \vdots\\ Z_{k,j}\\ \vdots\\ Z_{N,j}\end{array} \right]=\left[ \begin{array}{c}1\\ \vdots\\ 1\\ \vdots\\ 1\end{array} \right]_{Xystem,j}+\left[ \begin{array}{c}Z_{residual,1,j}\\ \vdots\\ Z_{residual,k,j}\\ \vdots\\ Z_{residual,N,j}\end{array} \right] \tag{7}\]  
Equation (7) can be abbreviated as:  
\[ \begin{array}{c}\mathbf{Z}_{j}=\underset{N\times 1}{\mathbf{I}}Z_{system,j}+\underset{N \times 1}{\mathbf{Z}_{residual,j}}\\ \end{array}  \tag{8}\]  
where \(\mathbf{I}\) denotes the unit matrix, and \(Z_{system,j}\) is the parameter to be estimated. It is assumed that \(Z_{residual,k,j}(k=1,\ldots,N)\) obey a Gaussian distribution and are uncorrelated with each other. The expected value and variance of \(Z_{residual,k,j}(k=1,\ldots,N)\) are \(E\left([Z_{residual,k,j}]\right)=0\) and \(D\left([Z_{residual,k,j}]\right)=\sigma^{2}\), respectively. Therefore, the variance-covariance matrix of \(Z_{j}\) can be written as:  
\[D\left(\underset{N\times 1}{\mathbf{Z}_{j}}\right)=D\left(\underset{N\times 1}{ \mathbf{Z}_{residual,j}}\right)=\sigma^{2}\underset{N\times N}{\mathbf{I}} \tag{9}\]  
where \(D()\) is the variance-covariance matrix. If we let \(\mathbf{V}\) donate the correction value of \(\mathbf{Z}_{residual,j}\) and \(\hat{Z}_{system,j}\) represent the estimated value of \(Z_{system,j}\), then formula (8) can be transformed as:  
\[\underset{N\times 1}{\mathbf{V}}=\underset{N\times 1}{\mathbf{I}}\hat{Z}_{system,j}- \underset{N\times 1}{\mathbf{Z}_{j}} \tag{10}\]  
From formula (9), we know that \(Z_{k,j}(k=1,\ldots,N)\) are equal-weight observations. Therefore, to obtain \(\hat{Z}_{system,j}\) based on the least-squares principle, we only need to satisfy \(\mathbf{V}^{T}\mathbf{V}=min\), in which \(\mathbf{V}^{T}\) is the transpose matrix of \(\mathbf{V}\).