2016, Vol. 59
2. 武汉大学测绘遥感信息工程国家重点实验室, 武汉 430079
2. State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, Wuhan 430079, China
全球连续GPS网络在不同的空间和时间尺度上为形变模式的监测提供了统一框架.连续GPS测量能够在全球板块运动、板块边界区域地壳形变以及局部同震变形等一系列空间尺度上揭示形变模式.形变模式的量化有助于从冰后回弹、季节性形变、瞬变运动等不同时间尺度来研究地球物理过程.在区域连续GPS网络中，存在着一种与时空相关的误差，称为共性误差(也称共模误差，Common Mode Error).共性误差的有效提取与剔除对于提高区域连续GPS坐标时间序列精度、改进形变模式分析具有重要的作用.
Dong等(2006)使用PCA和KLE联合滤波方法对南加利福尼亚综合GPS网络(SCIGN)5年单天坐标时间序列进行滤波，比stacking方法更有效地提取了共性误差.伍吉仓等(2008)和胡守超等(2009)也对该区域的共性误差进行了提取，前者使用了PCA和KLE法，后者还对stacking、PCA和KLE三种方法进行了比较.蒋志浩等(2010)采用PCA方法提取了CGCS2000框架下国家连续进行参考站系统(CORS)网络坐标时间序列中的共性误差，其结果表明使用空间滤波能够有效降低N、E、U三个方向的平均坐标重复性和U方向周年项振幅.Serpelloni等(2013)利用800多个连续GPS站点的坐标时间序列来研究欧洲地中海地区的垂向变形速率，利用PCA方法估计并剔除了坐标时间序列中的共性误差，使位移数据的信噪比有了明显提高.Shen等(2014)针对不连续的GPS时间序列提出了改进的PCA方法，从中国地壳运动观测网络(CMONOC)不完整的时间序列中提取了共性误差.虽然目前有关共性误差的研究已日趋成熟，但仍然存在一些问题值得探讨，例如当测站较少时，PCA或KLE方法得到的高阶主分量模式是否可以作为共有模式来提取共性误差.
一直以来，南极区域连续GPS网络由于受各种不利因素影响，测站稀少、数据质量差且空间分布十分不均匀.2007/08年国际极地年间，由多国合作开始在南极实施POLENET计划(The Polar Earth Observing Network)，该计划的主要任务是在极地布设GPS与地震仪等应用于地球科学研究的常年连续观测仪器(鄂栋臣和张胜凯，2006).近年来，随着POLENET计划的实施和南极IGS测站的不断扩充，南极半岛地区的连续GPS站在数量和空间分布上有了很大大改善，相对南极其他区域更为密集，研究该地区GPS网络的共性误差成为可能.
2 共性误差提取方法目前提取共性误差一般使用空间滤波方法，主要包括stacking、PCA和KLE方法.提取共性误差需要残差时间序列，各测站的残差时间序列由坐标时间序列剔除线性项和周期项得到，公式如下：
|
(1) |
其中，vi为第i个历元的观测噪声，也就是残差；y(ti)为站坐标时间序列；ti代表第i个历元时间，以小数年的形式表示；a代表初始历元测站坐标，b代表线性速度，c、d代表全年周期运动系数，e、f代表半年周期运动系数，g表示历元Tg处由各种原因引起的阶跃式偏移量(如天线位置变化，同震位移等)；H表示Heaviside阶梯函数.
2.1 堆栈法(stacking)假设一个区域网络有n个测站，站坐标时间序列长度为m天，建立一个m×n的矩阵 X，其元素x(ti，sj)代表第i天第k个测站的坐标残差值，那么该网络的共性误差可以表示为
|
(2) |
其中，δi，k是第i天第k个测站的中误差(Nikolaidis，2002). 由(2)式可知，该方法提取的共性误差即为各测站残差的加权平均值，故也称其为加权平均法.
堆栈法假设共性误差在空间均匀分布，能够对数百公里范围GPS网络有很好的近似.但随着空间范围扩大，共性误差会逐渐减弱至消失.
2.2 PCA方法(PCA)PCA方法是一种广义的空间滤波方法，在同震变形的研究中该方法已成功用于大地测量数据的分解(Scherneck et al.，2000)和区域滤波(Johansson et al.，2002).某一区域残差时间序列 X 为m×n矩阵，可定义其协方差阵 B 为
|
(3) |
B 是一个n×n的对称矩阵，对 B 进行正交分解可得
|
(4) |
其中， V T是一个行向量正交的n×n特征向量矩阵，Λ 是具有r个非零特征值的对角矩阵(n≥r).对于大多数大地测量数据，矩阵 B 是满秩的，即n=r.由线性代数可知，任何秩为n的矩阵都由n维正交基扩充得到.因此，可以使用正交基 V 来扩充矩阵 X ：
|
(5) |
即
|
(6) |
其中，A 是m×n矩阵，vk，j是第j个特征向量第k个元素.由于正交矩阵的性质，(5)式可以变为
|
(7) |
即
|
(8) |
(6)式中ak(t)称为矩阵 X 的第j个主分量，代表时间变化；vk(s)是其对应特征向量，代表对主分量ak(t)的空间响应.上述分解方法称为经验正交函数(EOF)分析，也称为主分量分析(PCA)(Preisendorfer，1988). 将特征向量按照特征值大小进行降序排列，排在前面的几个主分量对残差时间序列 X 的贡献最大，通常与共源时间函数有关；而排在后面的主分量通常与当地或个别测站的影响有关.那么，PCA方法的共性误差可以定义为
|
(9) |
其中，p为定义共性误差的主分量个数.由(9)式可知，PCA方法定义的共性误差在不同的测站上是不同的，它舍弃了堆栈法空间均匀分布的假设，而让数据本身去揭示共性误差的空间分布.
2.3 Karhunen-Loeve展开法(KLE)该方法与PCA方法类似，只是将协方差矩阵 B 标准化得到相关系数矩阵 C，其元素：
|
(10) |
其中，δ是残差时间序列矩阵 X 的标准差向量，即
|
(11) |
相关系数矩阵 C 也可以正交分解为
|
(12) |
Λ C、 W 分别是 C 的特征值组成的对角矩阵和特征向量组成的正交矩阵.
进而得到主分量矩阵 A ：
|
(13) |
同理，Karhunen-Loeve展开法的共性误差可以定义为
|
(14) |
为了提取和研究南极半岛区域的共性误差，本 文选取了该区域11个GPS连续跟踪站，包括5个IGS测站和6个POLENET测站，时间跨度为2010—2014年(5年)，测站分布如图 1所示.采用GAMIT/ GLOBK10.5对GPS数据进行解算和平差.GAMIT 基线解算策略如下：解算类型为松弛解(RELAX)，同时解算卫星轨道、测站坐标；观测值使用无电离层组合、自动修复周跳模式(LC_AUTCLN)； 测 站N、E、U方向先验坐标约束：IGS站约束为0.05 m、0.05 m、0.10 m，其他站约束为100 m、100 m、100 m；截止高度角为15°；历元间隔为30 s；对流层延迟改正使用 Saastamoinen模型；潮汐改正使用FES2004模型；太阳光压摄动改正使用BERNE模型；参考框架为ITRF2008；其他均采用GAMIT10.5默认设置.
|
图 1 南极半岛区域GPS测站分布图，圆点表示IGS站，三角形表示POLENET站 Fig. 1 The distribution of GPS stations in Antarctic Peninsula，the dots represent the IGS stations， and the triangles represent the POLENET stations |
基线处理完成后，再使用GLOBK在ITRF2008 框架下进行网平差.由于本文解算的是区域GPS网络，为了维持框架稳定，平差时加入斯克里普斯轨道和永久阵列中心(SOPAC)提供的全球子网的H文件；选用IGS重处理结果(IGb08)中91个核心站作为框架稳定站；对稳定站坐标N、E、U方向分别约束为0.01、0.01、0.05 m；对极点水平位置和误差约束为0.25、0.25、0.1、0.1 m，UT1及误差约束为0.25、0.1 s.部分GPS测站因为地震或天线变动等原因发生了阶跃，使用GAMIT10.5最新的地震重命名文件itrf08_comb.eq进行改正，从而得到各跟踪站单天站坐标时间序列.
3.2 残差时间序列得到单天坐标时间序列之后，通过公式(1)给出的模型来剔除趋势项和周期项等，从而得到单天残差时间序列.具体方法是：
1) 将GLOBK平差结果(org文件)利用sh_plotcrd 的-mb命令生成各测站的mb文件，再利用sh_mb2cats命令将mb文件转换成cats格式文件.
2) 将cats格式文件输入到CATS软件，噪声模型选择白噪声+闪烁噪声+随机漫步噪声(李伟伟等，2014)，估计各测站的6个系数(公式(1)中a、b、c、d、e、f).
3) 将各测站的系数带入公式(1)，由于在平差中使用地震重命名文件改正了阶跃项，所以阶跃项不再考虑，从而得到各测站的单天残差时间序列.
3.3 异常值剔除和时间序列插值时间序列中的异常值会对时间序列分析造成偏差，所以有必要对异常值进行剔除.首先对坐标时间序列使用三倍四分位数差法则进行孤立值剔除(Nikolaidis，2002)，然后剔除线性项和周期项，得到残差时间序列.由于残差时间序列是在0附近波动，因此采用阈值法，超过某一阈值的的残差被当作粗差予以剔除.剔除标准是：N、E、U三个分量的中误差阈值分别设为50、50、100 mm，残差阈值设为100、100、300 mm.残差时间序列中还存在间隙，无法直接提取共性误差，需要进行插值.本文参考Dong等(2006)的插值方法，结合实际情况(测站较少)对其插值方法进行了相应修改，以减少迭代次数.具体方法是：首先使用三点拉格朗日法(间隙小 于等于3天)和空间平均法(间隙大于三天)对站坐标时间序列进行初步插值；然后使用PCA方法对初步插值的时间序列进行主分量分析，取前三个主分量与缺失测站空间响应乘积的和作为该测站新的插值；迭代这一过程直至前、后结果之差的平均值小于0.01 mm.
3.4 共性误差提取堆栈法一般不考虑测站的特异性，所有测站的残差按中误差加权平均即可得到共性误差，接下来主要介绍PCA和KLE方法的提取过程.为了使不同主分量之间便于比较，我们把每个特征向量除以其绝对值最大元素，得到标准化空间特征向量，也称标准化空间响应，正值代表对 主分量的积极响应，负值代表对主分量的消极响应.将主分量乘以对应标准化因子，可得到标准化的主分量.
使用PCA方法对11个测站的残差时间序列进行主分量分析，分别得到了N、E、U方向第一、二、三的标准化主分量时间序列及其对应的标准化空间向量，如图 2中a、b、c所示.通过对比三个方向的标 准化空间向量，发现N方向第三主分量(图 2c)OHI2、OHI3、ROTH三个测站具有异常大的空间特征向量，其绝对值比其他测站平均大了约80%，且OHI2、OHI3与ROTH空间特征向量的方向相反，在E方向第三主分量模式以及U方向第二、三主分量模式也发现了类似的现象.考虑到OHI2、OHI3两个测站位于南极半岛东北端，而ROTH位于南极半岛西南部的阿德莱德岛，在空间上与其他测站的距离都相距较远，因此这三个测站有可能存在与其他测站不一致的强烈局部效应，这里先将它们从GPS网络中剔除，以便进行下一步工作.利用 KLE方法进行主分量分析，得到的N、E、U三个方向的第三主分量模式(图 2d)与PCA第三主分量模式十分相似，从某种程度上也验证了PCA方法的适用性.
|
图 2 11个测站N、E、U方向的标准化主分量序列(上)及其对应的标准化空间响应(下)，其中a、b、c列为PCA分析得到的第一、二、三主分量模式，d列为KLE方法分析的第三主分量模式. Fig. 2 The scaled principal components and normalized spatial eigenvectors of 11 stations in north，east and up component，which a，b，c are the first，second，third mode of PCA solution respectively，and d is the third mode of KLE solution |
利用PCA和KLE方法对剩余8个测站重新进行主分量分析，得到N、E、U方向的第一、二标准化主分量序列及其标准化空间响应，如图 3所示.在剔 除OHI2、OHI3、ROTH之后，两种方法的结果中N、E、U三个方向的第一标准化主分量序列振幅都有所减小，但标准化空间响应更加一致.在PCA方法得到的第二主分量模式中，南极半岛东部和西部测站呈现出方向相反的空间响应，KLE方法也发现类似的现象，说明南极半岛可能存在着对东部和西部影响相反的某种误差源，再加上该模式中位于内陆的FONP测站的标准化空间响应远小于其他沿海测站，本文认为这种现象可能是该区域海洋潮汐未完全改正或改正过度引起的.
|
图 3 8个测站N方向的标准化主分量(上)及其对应的标准化空间向量(下)， 其中a、c列为PCA方法得到的第一、二主分量模式，b、d列为KLE方法得到的第一、二主分量模式 Fig. 3 The scaled principal components and its normalized spatial eigenvectors of 8 stations in north component，which a，c are the first，second mode of PCA solution respectively，and b，d are the first，second mode of KLE solution respectively |
对共性误差的定义，尚无明确的标准和共识.Dong等(2006)利用共有模式来定义共性误差，即某一主分量模式中大多数测站(50%)具有明显(>25%)的标准化空间响应，且该模式的特征值超过了所有特征值总和的1%，可认为是共有模式.根据这个定义，前几阶主分量模式都有可能是共有模式，而高阶模式通常只与部分测站有关并且可能反映局部效应，因此不会成为共有模式.为了确定PCA和KLE方法提取共性误差要使用的共有模式主分量个数p，本文参考Dong等(2006)共有模式标准，对各主分量模式从标准化空间响应和特征值两方面分别进行显著性检验.在标准化空间响应方面，由图 3也可以看出，利用PCA和KLE方法得到的第一、二主分量均能够满足50%以上的测站标准化空间响应大于25%的条件，而更高阶主分量并不满足该条件，这里不再列出.在特征值方面，表 1给出了第一、二主分量特征值占总特征值的百分比，可以看出PCA和KLE分析得到的第一、二主分量模式均满足特征值百分比大于1%的条件.为了更加客观、准确地提取共性误差，本文还分析了各主分量对残差的贡献，主分量模式特征值的百分比反映了该 主分量模式对整个残差时间序列的贡献，图 4给出了两种方法得到的各主分量特征值的累计百分比，两种方法特征值累计百分比都在80%左右，且前两个主分量的特征值远大于其他主分量的特征值，由此可见前两阶主分量对整个残差时间序列的贡献十分可观.总体来说，在PCA和KLE方法中第一、二主分量模式都可以作为共有模式.将PCA和KLE方法获取的第一、二主分量序列及其对应空间响应分别代入公式(9)和(14)，其中p=2，就得到了相应的共性误差.
|
表 1 第一、二主分量特征值占总特征值的百分比 Table 1 the percentage of the first and the second principal component eigenvalues |
|
图 4 各主分量特征值占总特征值的累计百分比，左为PCA分析的结果，右为KLE分析的结果 Fig. 4 Cumulative percentage of PC eigenvalues：(left)PCA result；(right)KLE result |
由于堆栈法提取的共性误差是空间均匀分布，而PCA方法和KLE方法得到的共性误差在空间上并非均匀分布，为了便于与堆栈法进行比较，对后两者得到的共性误差在空间上求平均，从而得到了与堆栈法一致的共模时间函数，即一维共性误差.图 5对三种方法提取的一维共性误差进行了比较，可以看出三种方法提取的一维共性误差随时间变化趋势十分相近，特别是PCA和KLE两种方法的结果几乎一致，水平(N、E)方向的一维共性误差最大量级约为7 mm，振幅绝对值的平均量级约为1.5 mm；而垂直方向一维共性误差的最大量级约为24 mm，振幅绝对值的平均量级约为5.2 mm.三种空间滤 波方法振幅绝对值的平均量级差别很小(＜0.1 mm)，PCA结果比KLE和stacking的结果要稍微大一些.
|
图 5 分别使用stacking、PCA和KLE方法提取的共模时间函数(一维共性误差).其中，灰色实线代表stacking提取结果，黑色圆点代表PCA提取结果，黑色”+”号代表KLE提取结果 Fig. 5 Comparison of the common mode time function extracted by stacking、PCA and KLE approach. The gray solid line represents stacking extraction results，black dots represent PCA extraction results，black “+” represents the KLE extraction result |
从残差时间序列 X 中减去共性误差ε，对于堆栈法，由于其共性误差在空间均匀分布，公式如下：
|
(15) |
对于PCA方法和KLE方法，由于其共性误差并非均匀分布，公式如下：
|
(16) |
为了统计和对比各方法共性误差的剔除效果并减弱个别测站异常残差值的影响，对剔除前后残差的绝对值求和再除以测站数，得到单天L1范数离散度序列；对剔除前后的残差平方和除以测站数再开平方，得到单天L2范数离散度序列，残差的振幅定义为L1范数的离散度，残差的功率定义为L2范数离散度的平方.图 6给出了不同方法剔除共性误差前后残差L1范数离散度(即残差振幅)的比较，由图 6可知PCA和KLE方法滤波后的残差振幅明显小于Stacking方法的结果，说明PCA和KLE方法对南极半岛GPS数据的适应性更好.N、E、U三 个方向滤波后残差L1范数离散度序列中冬季(南极)的振幅明显大于夏季的，这可能与受温度影响的误差源有关，例如低层大气环流、小尺度温度扰动、天线热噪声以及热弹性张力(Prawirodirdjo et al.，2006)等.
|
图 6 stacking、PCA和KLE三种方法剔除共性误差后的L1范数离散度序列比较 Fig. 6 Comparison of the daily L1 scatters of the regionally filtered time series by the stacking，PCA and KLE approach |
表 2比较了不同方法剔除共性误差前后残差的平均振幅、平均功率(对历元求平均)和平均RMS(对测站求平均)，表 3统计了残差平均振幅、平均功率和平均RMS的下降幅度.由表 2、3可知，三种空间滤波方法在滤波后残差的平均振幅、平均功率以及平均RMS都有了明显的下降，PCA方法的效果最好，其次是KLE方法.
|
|
表 2 不同方法滤波后残差平均振幅、平均功率和平均RMS的比较 Table 2 Comparison of the daily L1 and L2 scatters and average RMS of the regionally filtered time series by different approaches |
|
|
表 3 不同方法滤波后残差的平均振幅、平均功率和平均RMS的降幅比较 Table 3 Decline comparison of amplitude and power of the daily scatters and average RMS of the regionally filtered time series by different approaches |
从站坐标时间序列中剔除共性误差的公式与从残差时间序列中剔除共性误差类似，只需将公式(15)、(16)中残差时间序列替换为坐标时间序列，共性误差采用PCA方法提取的结果.图 7给出了剔除共性误差前、后各测站坐标时间序列，其中各测站纵轴相对坐标是指相对于各测站的2010年1月1日初始坐标的坐标.由图 7可以看出，滤波后的坐标时间序列(黑点)相比滤波前的序列更加收敛，其中U方向最为明显，说明使用PCA方法对坐标时间序列进行空间滤波能够有效降低各站点坐标时间序列的不确定性.
|
图 7 PCA滤波前(灰点)、后(黑色点)的各测站相对坐标时间序列 Fig. 7 The unfiltered and filtered position time series by PCA |
利用CATS软件重新估计空间滤波前、后的坐标时间序列的线性项和周期项参数，噪声模型选择白噪声+闪烁噪声+随机漫步噪声，各测站的参数 估计结果在表 4中给出.表 4中各测站时间序列参 数的误差均有明显的下降，表明坐标时间序列剔除共性误差后能够大大减小估计线性项、周期项的误差，从而提高坐标时间序列的精度，使原来淹没在误差里的信号更加明显.
|
|
表 4 测站垂向坐标时间序列使用PCA方法滤波前后斜率、振幅和初相对比 Table 4 The unfiltered and filtered slope，amplitude and initial phase by PCA in vertical position time series |
共性误差作为一种区域时空相关的误差，并不是特指某一种误差，它是某一区域大部分站点都存在的一种或多种误差的集合，通过空间滤波，我们仅能得到这些误差的总和，却无法确定它由哪些误差构成以及哪些原因造成的.对时间序列进行谱分析，可以对不同周期的信号进行分离，有助于分析和解释共性误差的构成.由于前面三种滤波方法中PCA滤波后的残差离散度功率降幅最大，这里选用PCA方法得到的共性误差进行分析.
为了便于在时间域进行分析，对该方法提取的二维共性误差在空间域取平均(即在测站间取平均)，得到一维共性误差，然后再进行傅里叶变换，得到N、E、U三个分量的振幅谱.为了便于统计周期，以周期取代频率作为X轴坐标，单位为天，如图 8所示.我们发现N、E、U三个方向的最大振幅的周期分别出现在89.1、170.8和227.7天，对应的最大 振幅分别为0.56、0.46和1.71 mm·a-1，而且89.1 天的周期同时存在于N、E、U三个方向，170.8天的周期同时存在于E、U两个方向，227.7天的周期存在于N、U两个方向，另外仅U方向存在周期50天振幅1.5 mm·a-1的信号，这些周期由于振幅足够大，可以认为是在共性误差中的主要成分.其中部分周期可能是由于参考框架、大气质量负荷、未建模的卫星轨道、EOP和交点周期(Wdowinski et al.，1997)等不准确而引入的系统误差，这些周期的具体 成因目前尚无法解释，还需要进一步的研究和探 讨.对于U方向，还存在周期为9.4天(振幅0.69 mm·a-1)和13.66天(振幅0.90 mm·a-1)的短周期信号，考虑到海洋潮汐中存在9.1、9.6和13.7天(例如Mf分潮、半月太阴潮周期为13.7天)等周期，推测共性误差组成成分中可能含与海潮相关的误差源，该误差源可能是由南极半岛区域海洋潮汐未建模或模型不准确引起的.
|
图 8 PCA方法提取共性误差的幅值谱分析 Fig. 8 The amplitude spectrum analysis of the common mode error extracted by PCA method |
本文使用了stacking、PCA和KLE三种空间滤波方法对南极半岛区域GPS网络坐标时间序列进行区域滤波，其中PCA和KLE方法使用了前两阶主分量.通过比较发现三种方法在南极半岛区域均能有效提取共性误差，但PCA方法的滤波效果明显优于stacking，且略微优于KLE方法.使用PCA方法进行空间滤波后，残差时间序列的振幅、功率和RMS均有大幅下降，坐标时间序列的线性项和周期项误差也明显减小，说明进行空间滤波能够有助于提高坐标时间序列参数估计的精度.
通过对共性误差进行分析，本文总结了共性误差的一些性质：共性误差是一种与区域空间尺度相关的误差，随着区域范围的扩大，原本的部分共性误差因不适用整个区域会变为局部效应引起的特性误差；共性误差并非恒定不变，而是随时间不断发生变化；共性误差并非均匀分布，因为同一因素引起的误差对不同站点的影响大小并不相同.共性误差的谱分析结果表明其U方向存在明显的9.4、13.7天短周期信号，可能与南极半岛区域海潮未建模或模型不准确有关.