GPS坐标时间序列具有时间相关性和空间相关性^[1-2]。时间相关性一般采用最大似然估计法进行计算，而空间域误差(共模误差，common mode error，CME)则采用数据后处理方法予以削弱，称之为区域时空滤波^[3-7]。文献[3]最早提出堆栈滤波法来减弱空间相关性。但堆栈滤波假设CME在区域内均匀分布，对于空间跨度较大的观测网该假设不成立^[6]。文献[4]针对堆栈滤波的不足提出利用主成分分析法(PCA)来去除CME。相对堆栈滤波，PCA在理论上更加严密，但仍存在不足。首先，PCA基于二阶统计量(方差和协方差)，其潜在假设为CME服从高斯分布，但CME包含有色噪声，并不满足正态分布^{[1, 8]}。因此，PCA没有充分利用CME的高阶统计信息。其次，PCA滤波时，由于是将各正交分量的方差最大化作为准则，可能会导致提取的某一“数学上”的主成分中包含不同的物理模式，进而导致虚假特征。该问题称为“混合”或“聚类”现象^[9]。最后，全球各区域GPS观测网的CME特征并不一致，采用前多少个主分量(principal component，PC)目前仍无客观标准，这就给计算CME带来较大的不确定性。

本文将GPS网的CME提取看作一盲信号分离问题，引入独立分量分析(ICA)方法提取CME。相对于PCA，由于ICA引入高阶统计量，能够分离出统计独立的非高斯信号。本文首先介绍PCA和ICA的理论背景，然后利用GPS数据分析结果进行模拟并加入已知的共模误差对ICA提取共模误差的有效性进行验证，最后与PCA结果进行对比。

1 分析方法 1.1 主成分分析(PCA)

设观测网某一区域m×n维矩阵为X，X的列向量为各测站坐标时间序列，其中m为观测历元数，n为测站数，m>n。设X去趋势、去均值后的矩阵为，将进行奇异值分解可得：

(1)

式中，S为奇异值s_i构成的对角阵，i=1, 2, …, k，k为奇异值个数；V的列向量称为主轴或主方向，中的列向量则称为主分量。定义协方差阵其特征值为λ_i=s_i²/(m－1)，对应主分量的方差。将U和V按C的特征值由大到小对应排列，则前几个方差较大的主分量代表了中的“信号”，也即CME；高阶的主分量则主要代表了的噪声。假设前l(l＜k)个主分量代表CME，则CME按下式计算^[4]：

(2)

1.2 独立分量分析(ICA)

假设有m个不同的时间序列x，为n个相互独立的未知信号s的线性叠加，用混合矩阵A表示。假设该混合模型为一个瞬时模型，也即忽略该模型中任何时间延迟的信号混合，在某一历元t(t=1, …, T)的瞬时数学形式为：

(3)

其中，x_t为t时刻m×1维观测向量；s_j^t为t时刻第j个信号源；a_ij为第j个信号源对第i个观测值的贡献系数，且满足1≤i≤m，1≤j≤n。考虑所有观测时刻T，写成矩阵形式为：

(4)

式中，X=[x₁ x₂ … x_T]，S=[s₁ s₂ … s_T]。

为了恢复原始信号，需要求出A的逆，也称解混矩阵W：

(5)

由于CME在空间域具有相关性，本文采用tICA来提取CME。tICA假设各信号源在时间域相互统计独立，但对应的空间模式并不施加统计约束。本文采用FastICA算法来计算tICA^[10]。该算法以各估计的信号源之间的互信息最小化为目标函数，使用基于牛顿迭代的固定点优化方案，其计算效率比传统的梯度法快约10~100倍，且稳健性好。

由(4)和(5)式估计出以及，则观测网CME依据下式计算：

(6)

其中，R是中具有“共同空间特征”的IC子集，为与对应的空间响应(SR)。

2 数值分析 2.1 模拟数据

为了验证ICA法提取CME的有效性，本文利用CMONOC I(crustal movement observation network of China) 26个基准站信息生成的模拟数据进行分析。模拟GPS时间序列采用Fakenet软件包生成，其中包括长期速度、周年振幅信号^[11]。模拟CME时，假设其为白噪声和闪烁噪声的组合，其噪声模型参数值来自文献[8]，各测站速度、周年振幅等参数来源于文献[12]。为了专注于对CME的分析，本实验没有考虑粗差、阶跃以及数据缺失等因素，也没有模拟系统误差以及非线性瞬时信号。将模拟的CME分别加入到26个基准站时间序列上，最终得到包含CME的坐标时间序列。模拟的CME和BJFS站的时间序列如图 1所示。

将模拟数据E、N、U分量的残差序列(线性速度+周年项模型)分别去均值后，对其协方差阵进行奇异值分解，其特征值谱和累积方差如图 2所示。可以看出，E、N、U分量前3个PC的方差之和分别占总方差的36.57%、36.51%和37.03%，其中PC₁的方差对总方差的贡献分别为19.90%、19.33%和22.30%，表明各分量的“信号”主要集中在PC₁上。

为进一步评估E、N、U分量前多少个PC具有统计显著性，本文采用PA(parallel analysis)进行检验。PA本质上是一蒙特卡洛方法，通过随机模拟大量服从正态分布的观测数据，求得各特征值在一定显著性水平α下的边界值。若观测值矩阵第i个特征值大于对应的边界值，就认为该特征值是显著的^[13]，则取前i个PC及对应的SR进行分析。经PA检验表明，在1%的显著性水平下E、N、U3个分量前3、5、5个PC是统计显著的，如图 3所示。

2.2 结果分析 2.2.1 PCA法的结果

基于上述检验结果，分别对E、N、U分量进行PCA滤波。本文仅给出U分量结果，如图 4所示(SR是归一化后的结果，即绝对值最大的响应为±100%, 下同)。可以看出，SR₁的符号均相同，且其量级几乎一致，而其他SR的符号和量级则表现出一定的随机性，因而高阶的PC主要包含了与测站相关的局部效应的影响。虽然在本实验中，E、N、U分量SR₁中大于25%的测站刚达到50%，但具有明显的“共同特征”，因此本文分别取E、N、U分量的PC₁计算CME。该结果与此模拟实验的设置值相吻合，同时也间接证实了文献[4]的结论，即当区域GPS网的PC₁占主导时，PCA法的结果与堆栈滤波相一致。

2.2.2 ICA法的结果

分别取E、N、U分量统计显著的PC作为输入数据，采用FastICA法进行ICA分解，其中U分量恢复的独立信号IC和其SR如图 5所示。与PCA不同的是，ICA本身无法确定各IC的符号以及各IC之间的次序，但本文的目的是利用ICA提取CME，ICA的不确定性并不影响最终的结果。从图 5可以看出，U分量的IC₁对应的测站具有几乎一致的SR，而其他IC对应的SR则无明显空间分布特征。据此分析，IC₁对应的SR分别代表了CME在U分量的“共同空间”特征。因此，本文采用IC₁来计算CME。

2.2.3 PCA与ICA结果对比

表 1给出了PCA、ICA提取的CME与真值的Pearson相关系数。可以看出，PCA和ICA提取的CME与真值差异不大。该结果表明，当CME在残差中占主导时，ICA能够有效地提取CME，得到和PCA几乎一致的结果。表 2列出了PCA、ICA时空滤波前后残差时间序列RMS的变化情况，括号中的数值为RMS变化百分比。可以看出，PCA滤波后的残差RMS均值在E、N、U分量均比ICA滤波后的RMS要小，其中U分量差异较大。由前文所述，由于PCA是取残差序列方差最大的PC计算CME，因而滤波后的残差序列二阶统计量(RMS)要小于ICA法。

3 结语

CME是区域GPS观测网中一项重要误差源，目前广泛应用的PCA仅基于观测数据二阶统计量，丢失了高阶统计信息。本文假设CME与观测网其他误差统计独立，将CME提取看作是一个盲信号分离问题，利用ICA将其与其他误差源进行分离，并用半模拟、半实际数据验证ICA提取CME的有效性和正确性。本文结果为区域GPS观测网时空滤波提供了一个新的思路。