1 引言

CG-5相对重力仪是加拿大Scintrex公司制造生产的新型数字重力仪，采用微处理器装置，可以实现自动测量。传感器为无静电熔石英弹簧，设计精度为 5×10-8m/s2，读数分辨率为 1×10-8m/s2[1]。

相对重力测量计算中的改正包括固体潮、大气压、水文影响[2]、仪器高、零点漂移和格值改正[3]，潮汐改正为重力测量中最重要的改正，能达到±300×10-8m/s2[4][5]。海洋负荷潮对沿海地区、岛屿之间的高精度相对重力测量影响较大，能达到±10×10-8ms-2量级，并随着离海岸线距离的增加影响逐步减小。然而海洋负荷潮对内地影响较小，为1×10-8ms-2量级，可忽略不计[6][7]。地球固体潮是地球对月亮和太阳起潮力的响应。地球外部和内部形状时刻发生着周期性变化,对固体潮的研究可为检测地球内部动力学、反演地球内部构造提供有效途径，为研究地球内部振荡、地球液核共振参数和地球自转等提供有效约束，还可对环境效应和地表质量迁移做出解释，如陆地水、大气和海洋等[8]~[12]。计算重力潮汐参数一般采用Eterna标准方法[13][14]。Tsoft潮汐数据预处理软件,直接对各种干扰进行修正,并带有许多数据分析的工具,可进行滤波、谱分析、回归分析等,使观测数据预处理变得便捷直观[15]。另外，在重力测量中使用的弹性重力仪测量准确度主要受金属或石英弹簧的不稳定性限制,存在零点漂移问题。CG-5 熔凝石英弹性系统的漂移线性较好，可以精确确定其静态零漂率，从而可以有效改正零漂影响[16]。因此CG-5静态重力观测中的线性漂移及固体潮的研究和获取就显得尤为重要。

本文通过实验研究利用主成分分析（PCA）来求取CG-5静态观测数据的零点漂移和固体潮。将CG-5静态数据作为统计量，假设信号与噪声不相关，提取线性漂移常数和固体潮的影响。实验中使用了两组实验观测数据。首先用进行固体潮改正过的9天静态数据作为原统计量，使用主成分分析提取日漂移常数，并与最小二乘线性拟合求取的值做对比分析。然后用含有固体潮和零点漂移影响的观测数据，在进行线性改正的基础上使用主成分分析提取固体潮，并与CG-5重力仪自带的固体潮模型计算值对比，验证该方法的可靠性。

2 主成分分析

重力观测的主成分分析是基于协方差矩阵将重力数据作为原统计变量通过线性变换得到新变量，新变量之间统计不相关[17]~[21]。假设新变量是由CG-5重力观测数据的原统计变量经过线性变换得到，即：

   (1)

式中，E是由重力观测的统计变量G的协方差矩阵计算得到的线性变换矩阵。假设协方差矩阵为,则它是一个n×n阶对称矩阵。通过求解矩阵的特征根和特征向量，可以分解为

(2)

式中，D是对角阵，对角线上的元素λ1≥λ2…≥λn是矩阵的特征根，同时等于对应重力观测数据主成分的方差。E是由特征向量构成的正交矩阵[20]。因此，可由(1)得到

             (3)

前k个主成分共有多大的综合能力以及所包含的重力值信息量可用累积贡献率α表示

(4)

式中，λi表示特征值，p表示特征值总的个数，k表示选取的个数。

如果选取k个特征值较大的新变量重构原重力观测值统计变量，可以达到提取主要信息的目的[20][21]：

(5)

式中，表示由选择的特征向量构成的k阶矩阵。

3 静态观测数据

实验中静态重力数据是由CG-5重力仪(系列号140541221)于某一观测站的9天观测所得，时间从2014年8月15日至23日，以每天9点21分到次日9点20分为一个时间序列，即作为一个统计观测量。采用6HZ的采样率，加之仪器自身具有的地震滤波功能，对低频噪声能够进行过滤，并能舍弃高于6倍标准偏差的高频噪声，因此可对1分钟的360个样本数据进行平均得到一个读数，并进行了温度补偿和倾斜改正[1]。因此，每天可以采集同一个时间段（24小时）的1440个重力观测数据。

连续9天的重力数据由于来自每天的同一时间观测，必然含有类似的成分，包括信号和噪声。因为所含信号比如固体潮每天会有所不同，所以包含的信号尺度和噪声水平必然有所区别。假设信号与噪声统计不相关，尝试利用主成分分析，把15号到23号的时间序列作为9个统计量，来分离信号和噪声。

实验中使用了两组实验观测数据。第一组数据是进行了固体潮改正的静态观测数据，尝试通过主成分分析将信号和噪声分离来估算CG-5的静态零漂率；第二组数据是未进行固体潮改正的静态数据，尝试通过主成分分析将信号和噪声分离来提取固体潮等影响，用来对含有这些影响的静态数据进行改正。

由于观测站环境干扰较多，采集得到的重力观测数据一般含有少量高频干扰信号和毛刺，这些突变点造成数据的不稳定性，不利于观测数据的研究，故需要对数据进行预处理，即滑动平均。滑动平均相当于低通滤波器，就是让数据通过一个n点的滤波器。滑动平均法的参数选取直接影响着对数据的平滑效果。参数值过小则可能对低频随机起伏的误差未作平均而减小，不利于抑制随机误差，参数过大可能将高频变化的成分一起被平均而削弱[22]。通过实验并根据重力数据的实际变化情况，选取大小为5的移动窗口进行平滑处理能更好的达到平滑的目的，较好地解决了数据的突变点等不稳定问题。　

<!--[if !supportLists]-->4 <!--[endif]-->CG-5观测数据的主成分分析

4.1静态零漂率

首先用第一组数据进行实验研究， 9天数据作为原统计量，用主成分分析方法将噪声和信号进行分离。每天的数据作为一个统计量，每个统计量为1440×1个元素的列向量，这样9天的观测值就组成了1440×9的矩阵,即作为统计量。在得到统计量G之后，通过公式（2）计算其特征向量构成的线性变换矩阵E，再由公式（1），得到主成分的9个模式。研究可得，主成分分析实现了信号和噪声的分离，图1给出了所得到的9个特征根的值，它们等于对应重力数据主成分的方差,特征值很小甚至接近于零则认为不含有主要的信息，因此可以寻找到一个点来确定主成分的个数，使得剩下的特征根具有接近相同且相对很小的值[20]。因此从图1中可以看出相对于第一个特征值外其它的都非常小，均接近于零。并且通过计算得到第一个模式的贡献率为99.8%，包含了数据几乎所有的信息。再结合图1分析得出，第一模式可以认为是信号部分，而其他8个模式为噪声部分。

图1 九个模式的特征值

通过主成分分析第一模式PCA1为分离得到的信号，因此，我们选择该模式来重构统计量G，分别得到PCA估算9天的线性漂移率。由公式(5)，选择第一模式进行重构：

（6）

将其详细展开，则得到：

（7）

式中，分别为PCA估算15-23日9天的线性漂移，是矩阵E第一列向量中的元素,为主成分的第一模式，上标PCA表示主成分估算得到的零点漂移，可以看出，PCA对于15-23日的估算结果只有尺度大小的区别。为了说明PCA估算的线性漂移常数的可靠性，利用最小二乘拟合（Least Squares Linear Fitting,简LSLF）计算值作比较分析，即对原始观测数据使用最小二乘拟合原理求出每天的零点漂移值。表1给出了PCA与LSLF的对比。

表115-23日9天数据零点漂移计算结果表/（mgal/day）

注1 PCA、LSLF分别表示使用PCA分析和LSLF所得的日漂移常数；P-L为PCA与 LSLF作差。

从表1中可以看出，PCA与LSLF得出9天日漂移差值均在10-2μgal/day量级，可以说是相差甚微，具有非常好的一致性。但是PCA剔除了噪声，效果更好，因此可以说，PCA更适用于CG-5静态数据的日漂移常数的估算。从表1中还可以看出，从15日到23日的日漂移常数有明显增加，说明CG-5重力仪的零点漂移具有随时间在不断增加的特点。

4.2 固体潮

根据前面PCA所得的日漂移常数先对第二组数据进行线性改正，得到经线性改正但未进行固体潮改正的一组数据，作为原统计量，对其尝试利用主成分分析，分离信号和噪声，来实现估算固体潮。每个统计量为1440×1个元素的列向量，这样9天的观测值就组成了1440×9的矩阵，即作为统计量。

在得到统计量G之后，通过公式（2）计算其特征向量构成的线性变换矩阵E，再由公式（1），得到主成分的9个模式，研究可得，主成分分析实现了信号和噪声的分离，图2给出了所得到的9个特征根的值，它们等于对应重力数据主成分的方差[20]。因此从图2中根据特征值的相对大小可以确定前三个为主要成分，并通过进一步计算，第一模式的贡献率为65.2%，前两个模式的累积贡献率为94.8%，前三个模式的累积贡献率达到99.9%，几乎包含了重力数据的全部信息。

图2 九个模式的特征值

因此，主成分的前三个模式为分离得到的信号，我们选择该模式来重构统计量G，分别得到PCA对应的每一个模式的值。图3给出了线性改正后的9天数据时序变化图。从图3中可以看到，15和16号的波形变化相近，17和18号的波形变化相近，19到23号的波形变化相近.。又因为受固体潮影响,重力观测曲线起伏波动,每天可观测到两个波峰、两个波谷[23]，因此可以初步推测波形不同是这9天的固体潮、海潮等影响不同所导致。然而通过计算PCA的前三个模式计算的累积量（具体计算见公式9），与CG-5自带固体潮模型计算的固体潮值(Solid Tide Model,简STMOD)比较发现，前三个模式叠加并进行中心化处理后的结果（PCA），正好对应于固体潮影响量。

于是，由公式（5），选择前三个模式进行重构：

(8)

将其详细展开，则得到：

(9)

式中，分别表示PCA估算的15-23日9天的固体潮改正量，是矩阵E第i列向量中的元素，PCAi为主成分的第i模式，i=1,2,3。上标PCA表示主成分估算得到的固体潮值，可以看出，PCA对于15-23日的估算结果只有尺度大小的区别。表2给出了15-23号9天两种固体潮统计和比较结果。

图3 原始重力数据时序图（线性已做改正）

注2 D15-D23表示15-23日每天的数据

表2 固体潮于15-23日的PCA与STMOD统计结果/mgal

注3PCA为主成分分析计算结果；STMOD为CG-5内部软件固体潮模型计算结果；P-S为PAC与STMOD作差

从表2可以看出，PCA和STMOD得到的固体潮值的统计结果中，相差最大的仅为5μgal，具有较好的一致性，可以说PCA可用于静态观测中的固体潮的估算。15-23日固体潮的平均值均在零附近，STMOD计算的均值稍微偏离零轴，但在整体上正负值具有较好的对称性；15-18日的固体潮高潮最大值达到90μgal，低潮最大值达到80μgal；19-23日的固体潮高潮最大值达到124μgal，低潮最大值高于70μgal。另外，结合图3中可以看出，固体潮具有一定的规律性，每天都会有两个高潮和低潮，并且每天的高潮和低潮都会有往后一段时间的推迟。PCA和STMOD在均方根和标准差方面几乎接近，仅相差在10-1μgal量级。另外，对于PCA和STMOD结果的差值P-S，最大值均不大于8μgal，最小值都在5μgal内，均值在零附近波动，波动范围均在3μgal内，均方根(RMS)均小于5μgal，说明PCA与STMOD计算结果具有很好的一致性。

经过以上的分析，并采用以上的计算结果对观测数据进行线性改正和固体潮改正，得到改正后的时序图如图4。CG-5重力仪的重复性标准差为5μgal[1]，而从图4中可以看出，经过改正后的重力数据变化浮动很小，均在5μgal范围内变化。另外，导致这种变化范围的可能还有比如海潮，零点漂移的非线性等其他一些微小影响量。最后可对数据进行平均得到该测站点的重力读数值。

图4 固体潮改正后的重力观测值时序图

注4 D15-D23表示15-23日的每天数据

5 结论

本文利用主成分分析，把9天的静态观测数据作为9个统计量，并假定信号与噪声之间统计不相关，从而对信号和噪声进行分离。文中采用了两组实验数据，首先对数据进行了平滑处理移除了高频干扰信号和大的突变点。实验中使用主成分分析提取线性漂移常数，第一个模式的贡献率为99.8%，包含了数据几乎所有信息。通过分析并由计算结果得出，第一模式可以认为是信号部分，而其他8个模式为噪声部分。用第一主成分进一步估算了线性漂移，并与最小二乘线性拟合求取的数值对比，结果仅相差在10-2微伽量级。但主成分分析剔除了噪声，效果更好。然后，数据在做了线性漂移改正的基础上求取固体潮，使用主成分分析实现了信号和噪声的分离，其第一模式的贡献率为65.2%，前两个模式的累积贡献率为94.8%，前三个模式的累积贡献率达到99.9%，几乎包含了重力数据的全部信息。因此前三个模式为信号部分，其余六个模式为噪声部分。并通过实验分析和计算得出，前三个模式的叠加量正好对应固体潮，并与CG-5内部软件提供的固体潮计算值对比，结果较为一致。两者计算结果作差的统计值均不大于8μgal，其中RMS都小于5μgal。因此主成分分析用于估算日漂移常数和固体潮等影响有显著的效果。最后，根据以上计算结果对观测数据做了零漂和固体潮改正后，变化均在5μgal范围内，原因可能是仪器的测量精度、海潮和零漂的二次项等微小影响，还需要将来进行深入研究。

参考文献

[1] SCINTREX LIMITED. CG-5 ScintrexAutograv System OPERATION MANUALV5.0, 2009: 99-172.

[2] Lederer M. Accuracy of the Relative Gravity Measurement[J]. ActaGeodyn. Geomater, 2009, 6(3):383-390.

[3] 肖凡,何志堂，张宏伟, 等.CG-5型相对重力仪测量精度分析[J]. 测绘技术装备, 2011,13(3):6-8.

[4] 王勇, 张为民, 王虎彪, 等．绝对重力观测的潮汐改正[Ｊ]. 大地测量与地球动力学, 2003, 23(2):65-68.

[5] 马玄龙, 肖毅, 刘磊. CG-5重力仪静态观测残差分析[J]. 资源环境与工程, 2010, 24(6): 698-700.

[6] 周江存, 徐建桥, 孙和平. 中国大陆精密重力潮汐改正模[J]. 地球物理学报, 2009,52(6):1474-1482.

[7] 张宏伟, 张松堂, 肖凡, 等. 海岛礁相对重力测量的潮汐影响[J].测绘信息与工程, 2012, 37(1):43-45.

[8] Crossley D, Hinderer J, Casula G, et al. Network of superconducting gravimeters beneLSLFs a number of disciplines[J]. Eos, Transactions American Geophysical Union, 1999, 80(11): 121-126.

[9] Sun H P, Hsu H T, Jentzsch G, et al. Tidal gravity observations obtained with a superconducting gravimeter at Wuhan/China and its application to geodynamics[J]. Journal of Geodynamics, 2002, 33(1): 187-198.

[10] 孙和平, 徐建桥. 基于全球超导重力仪观测资料考虑液核近周口共振效应的固体潮实验模型[J]. 科学通报, 2004, 48(6): 610-614.

[11] 孙和平, 许厚泽, 陈武, 等.香港地区重力固体潮和海潮负荷特征研究[J].地球物理学报,2006,49(3): 724-734.

[12] Dehant V, Defraigne P, Wahr J M.Tides for a convective Earth[J] . Journal of Geophysical Research: Solid Earth (1978–2012), 1999, 104(B1): 1035-1058.

[13] Wenzel H G. The nanogal software: Earth tide data processing package ETERNA 3.30[J]. Bull. Inf. MaréesTerrestres, 1996, 124: 9425-9439.

[14] 孙和平, 许厚泽, 陈武,等. 香港地区重力固体潮和海潮负荷特征研究[J]. 地球物理学报, 2006, 49(3): 724-733.

[15] Van Camp M, Vauterin P. Tsoft: graphical and interactive software for the analysis of time series and Earth tides[J]. Computers & Geosciences, 2005, 31(5): 631-640.

[16] 邢乐林, 李辉, 夏正超, 等.CG-5重力仪零漂特性研究[J].地震学报, 2010,32(3): 369-373.

[17] Wold S, Esbensen K, Geladi P. Principal component analysis[J]. Chemometrics and intelligent laboratory systems, 1987, 2(1): 37-52.

[18] Abdi H, Williams L J. Principal component analysis[J]. Wiley Interdisciplinary Reviews: Computational Statistics, 2010, 2(4): 433-459.

[19] 穆大鹏,郭金运,孙中昶, 等.基于主成分分析的GRACE重力场模型等效水高[J].地球物理学进展, 2014, 29(4): 1512-1517.

[20] Johnson R A, Wichern D W. Applied Multivarariate Statistical Analysis[M].Beijing: Tsinghua university press, 2008. 430-459.

[21] Jolliffe I T. Principal component analysis[M]. Springer, New York. 2012.

[22] 史文海,李正农,吴建佳.近地面强风不同间隔滑动平均统计特性的对比分析[J]. 空气动力学学报, 2013, 31(5):611-615.

顾兆峰, 张志珣, 杨慧良,等. KSS 31M 海洋重力仪静态观测结果及分析[J]. 海洋测绘, 2005, 25(2): 66-68.