## 2 Spatiotemporal Filtering Considering Formal Errors  
In a regional GNSS network of \(n\) stations, the residual position time series of the \(i\)th station at the \(k\)th epoch is expressed as:  
\[\mathbf{X}(t_{k},i)=\mathbf{S}(t_{k},i)+\mathbf{e}(t_{k},i)\qquad i=1,2,\cdots,n \tag{1}\]  
where \(\mathbf{S}\) represents the CME and \(\mathbf{e}\) denotes the remaining error. In order to extract the CME from the time series, the covariance of the residual time series is first generated with [24]:  
\[\left. \begin{array}{l}\mathbf{B}(i,i)=\frac{1}{m_{i}-1}\underset{k\in M_{i}}{ \sum}\mathbf{X}(t_{k},i)\mathbf{X}(t_{k},i)\\ \mathbf{B}(i,j)=\frac{1}{m_{ij}-1}\underset{k\in M_{i}\cap M_{j}}{\sum}\mathbf{X}(t_{k },i)\mathbf{X}(t_{k},j)\end{array} \right\} \tag{2}\]  
where \(m_{i}\) is the number of epochs of the \(i\)th station; and \(m_{ij}\) denotes the number of epochs of the intersect set \(M_{i}\cap M_{j}\). From this covariance construction, we can see that unit weight is adopted for all epochs of the observations of each station. Due to different surroundings, the formal errors derived from daily GNSS solutions are different epoch by epoch and station by station. Therefore, it is not reasonable to regard all epochs of observations with unit weight in forming the covariance with Equation (2). Since formal errors are related to observation error, but not necessarily related to CME [1], in order to directly use Equation (2) to form covariance, the time series is generated as follows:  
\[\mathbf{X}^{\prime}(t_{k},i)=\mathbf{S}(t_{k},i)+\frac{\sigma_{0i}\mathbf{e}(t_{k},i)}{ \mathbf{\sigma}(t_{k},i)} \tag{3}\]  
where \(\mathbf{S}\) and \(\mathbf{e}\) have the same meaning as Equation (1); \(\mathbf{\sigma}(t_{k},i)\) represents the formal error of the \(i\)th station at the \(k\)th epoch; and \(\sigma_{0i}\) is the standard deviation of unit weight of the \(i\)th station, which is determined based on the principal that the total power energy of \(\mathbf{e}\) is not changed. Thereby, the sum of the weight for the \(i\)th station must be equal to the number of epochs. For example:  
\[\sum\mathbf{P}(t_{k},i)=\sum\frac{\sigma_{0i}^{2}}{\sigma^{2}(t_{k},i)}=m_{i} \tag{4}\]  
where \(\mathbf{P}(t_{k},i)\) is the weight of the \(i\)th station at the \(k\)th epoch. From Equation (4), it is easy to derive the standard deviation of unit weight as follows:  
\[\sigma_{0i}=\sqrt{\frac{m_{i}}{\sum^{1}/\sigma^{2}(t_{k},i)}} \tag{5}\]  
Once \(\sigma_{0i}\) is determined, the weight factor \(\sigma_{0i}/\sigma_{i}(t_{k},i)\) in Equation (3) is obtained. The initial values of \(\mathbf{S}\) and \(\mathbf{e}\) in Equation (3) are determined using a traditional PCA approach. Then replacing \(\mathbf{X}(t_{k},i)\) with \(\mathbf{X}^{\prime}(t_{k},i)\), the covariance of the generated time series is formed with Equation (2). Through the eigen-decomposition of covariance, we obtain the eigenvectors and eigenvalues and rearrange them in descending order. With the sorted eigenvectors, the principal components (PCs) calculated in Equation (6) and the time series can be completely reconstructed from the PCs using Equation (7). \[\mathbf{a}_{k}^{\prime}(t_{i})=\sum_{j=1}^{n}\mathbf{X}^{\prime}(t_{i},j)\mathbf{v}_{k}^{ \prime}(j)\qquad k=1,2,\cdots,n \tag{6}\]\[\mathbf{X}^{\prime}(t_{i},j)=\sum_{k=1}^{n}\mathbf{a}^{\prime}_{k}(t_{i})\mathbf{v}^{\prime}_{ k}(j)\qquad j=1,2,\cdots,n \tag{7}\]  
where, \(\mathbf{a}^{\prime}_{k}(t_{i})\) is the \(k\)th PC at the \(i\)th epoch; and \(\mathbf{v}^{\prime}_{k}(j)\) denotes the \(k\)th element of the eigenvector of the \(j\)th station. However, Equation (6) is only appropriate for the non-missing data. When the missing data exists in the time series, according to [PERSON] et al.