2012, Vol. 55
2. 中国科学院研究生院地球科学学院, 北京 100049
2. College of Earth Science, Graduate University of Chinese Academy of Sciences, Beijing 100049, China
GOCE(Gravity field and steady-state Ocean Circulation Explorer)卫星是一颗用来反演静态高阶引力场的重力卫星[1-2]，其最基本的数据是梯度仪所观测的引力梯度数据.相对于其它观测数据而言，梯度数据对引力场高阶部分更加敏感，因此GOCE卫星的引力梯度观测数据对于解算引力场模型的高阶次项有着重要的作用.但由于仪器本身的缺陷，梯度仪无法给出全频段的高精度数据，仅在观测频段0.005~0.1Hz能够达到所需的精度，而在观测频段外，特别是低频部分含有非常大的误差，所以如何对引力梯度数据进行滤波在处理GOCE 数据时就显得十分重要，其原因是仅用有限的观测频段信号去反演全频段的引力场模型是不太可能的.如何评估和处理低频误差的影响是GOCE 数据处理最为关键的工作之一.
针对GOCE 数据如何滤波，已经有多篇文献进行阐述.Schuh等[3]讨论了在最小二乘解算模式下的去相关处理，其核心思想是根据误差功率谱，设计合理的离散滤波器系数，将有色噪声白化；同时为了解决滤波中的相位漂移问题，处理中不仅对观测值滤波，也同时对系数矩阵进行滤波.Haagmans等讨论分析了沿轨的滤波和对位系数的滤波两种方式[4].Reguzzoni 等[5] 主要介绍了威纳滤波器(Wienerfilter)在最小二乘配置方法恢复引力场位系数中的应用.Brockmann和Schuh等[6]基于已知的误差功率谱，先后对观测值进行高通滤波处理(Highpassfilter)､陷波滤波器(Notchfilter)滤波､自回归滑动平均滤波(ARMA)处理，从而达到了削弱有色噪声的目的.但实际的梯度仪误差功率谱是很难获知的，因此需要已有模型提供先验值然后作迭代处理.此外不少文章还讨论了梯度仪的校准问题[7-11]，但这些对处理引力梯度观测值中的低频误差的作用不大.由于原始数据的获取较为困难，国内在此领域的研究比较滞后.罗志才等[12]对GOCE 梯度数据的预处理工作进行了讨论；徐新禹等[13]分析了GOCE 梯度的测量误差，并基于AR 模型对其进行了模拟；吴云龙等[14]讨论了基于小波收缩阈值降噪的卫星重力梯度数据粗差探测方法；钟波[15]讨论了GOCE 沿轨重力梯度有色噪声的Wiener滤波预处理方法，并基于模拟数据予以了验证.综合已有的结果来看，Schuh 和Reguzzoni等所介绍的方法是当前最为流行的，前者主要应用于最小二乘求解，而后者则主要应用于最小二乘配置解算.事实上，官方所发布的Time-wise[16-17]和Space-wise[18-19]解算结果分别采用了这两种方法.总体来讲，上述滤波方法主要适应于最小二乘或配置算法解算重力场.
于锦海等[20]引入了利用引力梯度数据来恢复地球重力场的另一种途径，即:构建引力梯度不变量，并以此为基础建立了关于扰动位的边界条件，从而可以利用调和分析方法来恢复重力场，避免了大型方程组的解算.若将引力梯度不变量方法引入到GOCE 实际观测数据的处理中就必须讨论和建立与其相适应的滤波方法.本文正是基于这样的目标来研究针对GOCE 引力梯度数据的滤波方法.
在处理低频误差之前，首先要讨论低频部分对哪些阶次的引力场反演有影响.观测所说的低频误差主要是指时域上的误差.然而重力场模型的高低阶对应的是空间分辨率的高低频.越是低阶，空间上越是对应大尺度；越是高阶，空间上越是对应小尺度.这种空间上的尺度和卫星运行时的数据采集频率虽然有一定的关系，但并不是严格对应的.本文将试着从GOCE 卫星的运行规律入手，推导出引力梯度时间域上的表示公式，并粗略表示出频率与阶次的关系.然后通过模拟计算来分析相关特征.
在频率分析的基础上，如何进行滤波，显得极为重要和关键.本文将主要研究移去恢复法和零相位的数字滤波方法在GOCE 数据滤波中的应用.前者用以解决滤波中的低频信号损失问题，后者用以解决滤波过程中的相位漂移问题.最后本文将基于此种滤波方法对实际GOCE 引力梯度观测值进行滤波处理，得到相关的重力场模型.并与欧空局公布的结果进行对比，进一步地验证本文所提供的方法及相关结论.
2 频谱分析 2.1 基本分析引力位可用球谐展开表示如下[21]:
|
(1) |
其中R为地球的平均半径，(r，θ，λ)是计算点球坐标.考虑到GOCE 的实际情况，其位置可看作时间t的函数，从而θ，λ 可分别用ωθt，ωλt代替.由于GOCE 轨道接近球面，因此r可看作常数，不随时间变化.更进一步地考虑到GOCE 卫星的轨道倾角大致为96.7°，为讨论的方便现假定其轨道为严格的极地轨道，因此引力的径向梯度可表示为时间t的函数如下:
|
(2) |
其中ωλ 为卫星在东西方向的角速度ωλ ≈7.27220521×10-5rad/s，由于是极地轨道，其值为地球自转速度，ωθ 为卫星在南北方向运行的角速度ωθ=v/r≈1.212×10-3rad/s，v为卫星的运行速度.由于ωθ/ωλ=16.68，则从上式可以看出，在阶次n和m一定的情况下，径向梯度的时间频率主要受卫星运行速度的影响，而且各阶次的梯度含有各个频段的信号，换言之高阶中仍然含有低频信号.
为表明阶次与频率的对应关系，现做如下实验.利用EGM08[22]的第150 阶位系数，在GOCE 卫星轨道的实际位置处(来源于ESA 提供的2009年11月2日到5日之间的PKI数据)计算了该阶对应项的引力梯度在梯度仪坐标系下的六个分量，各分量的功率谱图形如图 1.
|
图 1 第150阶计算的引力梯度分量功率谱 Fig. 1 PSD of gravity gradient scalculated by 150th degree of EGM08 |
从图 1可以看出，第150阶引力梯度含有不同频率成分的梯度分量.由于150 阶可看作高阶，易见高阶引力场仍含有低频成分；同时也可以看出，与x轴相关的分量，高频更占优，这主要是由于x轴是卫星的运动方向，变化较快；从曲线的整体来看，频率与阶次虽不是一一对应的，但第150 阶球谐项的引力梯度分量在高频有明显的截止频率.下节将推导出特定阶次的球谐项与截止频率之间的对应关系.
2.2 阶次与截止频率的关系
|
(3) |
其中
|
|
而Flmk(I)为倾斜函数，ωs和ωe分别为卫星和地球的平均角速度.由此可知V(t)的第n阶径向梯度分量可表示如下:
|
(4) |
其中
|
很显然由于GOCE 卫星的轨道接近于一个球面，因此r(t)的时间变化不明显，故Vrr，nth(t)的时间频率主要由kωst+mωet来决定.根据前面的结论ωs$ \gg $ωe，则阶次与频率有如下近似关系:
|
(5) |
其中k为一中间量，取值范围为:[-n，n]，f为数据采集频率.进而得:
|
(6) |
易知，当阶次n固定时，最大频率为n/Ts，这就是截止频率.这表明各阶引力梯度分量含有不同频率成分的信号，其最大频率可用(6)式近似表示(k=n)，而最小频率很难给出明确表达，理论上可取0.因此第n阶引力梯度分量的频率分布范围为[0，n/Ts].具体到GOCE 卫星，其轨道周期Ts大约为5400s，因此可得表 1中结果.
|
表 1 关于阶数与频率之间的对应关系 Table 1 Example on the relationship between frequency and degree |
依据于锦海等[20, 25-26]提出的引力梯度不变量方法，可以利用下列的边值问题从GOCE 的引力梯度数据来恢复地球引力场
|
(7) |
其中T是扰动位，S是平均轨道球面，ΔB=B-B0是不变量扰动，B和B0 分别由观测值和正常位计算求得.由于问题(7)建立了关于扰动位的边界条件，因此根据球谐函数的正交性可以通过积分直接反演得到位系数.然而GOCE 引力梯度数据中存在着较大的低频误差，直接计算不可能得到很好的解，因而有必要对原始数据进行滤波处理.虽然Schuh 和Reguzzoni等所介绍的滤波技巧在相应的反演方法中十分有效，但却并不适用于此边值问题，例如Schuh方法中对方程两边的同步滤波，会导致球谐函数的正交性遭到破坏，调和分析方法将失效.因此本文研究了不同的滤波方法.
3.1 滤波方法介绍本文中，滤波主要采用移去恢复法和向前向后滤波.其中移去恢复法即首先用EGM08模型的前300阶在GOCE 的实际轨道处模拟出引力梯度观测值，然后将其与观测值作差，得到扰动值，并对此值进行滤波处理，滤波器为1000阶的带通有限脉冲响应滤波器[27](FIR)，然后用滤波后的扰动值进行引力场恢复.这样做的目的是用EGM08模型来提供滤波通带之外的信号.
由于通常的滤波方法会存在相位漂移，这对反演重力场十分不利.为了解决这个问题，本文采用了向前向后滤波处理，基本原理如下[28-30]:
|
(8) |
其中x(n)表示输入序列.滤波的原理即首先对x(n)顺序滤波得到y1(n)，然后将其反序滤波得到y3(n)，最终滤波结果即为y3(n)的逆序列.在频率域来看，可得:
|
(9) |
从此式可以看出，这种滤波方式不存在相位漂移[28-30]，同时滤波后的幅度变为之前的1/ H(ejω)2.
3.2 滤波效果现对ESA 提供的GOCE 卫星在2009年11月1日至11 月10 日时段的引力梯度观测值(EGG_NOM2)进行上述滤波处理，滤波器的带通频率为0.005~0.1Hz.滤波前后的功率谱如图 2所示.
|
图 2 滤波前和滤波后的功率谱对比 Fig. 2 PSD comparison.The filteris 1000 orders band pass FIR filter and time period of raw data is 2009-11-01-2009-11-102913 |
为了验证频率与阶次的关系式，本文也尝试通过改变滤波通带来检核公式(6)的合理性.实验中采用的通带频率为0.015~0.1 Hz，0.038~0.1 Hz.滤波后的功率谱如图 3所示.
|
图 3 通带频率为:0.015~0.1 Hz和0.038~0.1 Hz滤波后的功率谱 Fig. 3 PSD of residual data filtered by 1000 orders band pass FIR filter.Time period of raw data is 2009-11-01-2009-11-10 |
从图 2和图 3易知，滤波通带内的信号基本不变，而通带外的信号被滤除到仅占非常小的比例，基本可忽略.本文引力场反演所用的信号即为滤波后的扰动值，很显然此时的信号主要由滤波通带内的信号所控制，因此从滤波前后的功率谱可知，此种滤波方法达到了滤低频误差同时保留通带内信号的目的.
3.3 反演结果本文利用上述滤波方法对2009年11月1日到2010年1月1日的GOCE 引力梯度数据进行了滤波处理，并利用公式(7)所建立的观测方程进行调和分析得到位系数，现给出不同滤波下的反演结果.
图 4中的GOCE Time Wise Solution[31]是欧空局公布的结果之一，从此图可以明显看出，本文反演的模型与欧空局官方提供的模型在高阶是比较接近的，可以得出本文的滤波方法是有效的.特别需要说明的是图 4中显示的利用本文提供的滤波方法解算的重力场位系数是滤波后的直接结果，并没有做其它任何处理.对于图中的低阶部分与EGM08 差异较小，这恰好表现了移去恢复法的特点:低频信号主要由先验模型来控制；同时不同滤波通带下，低阶部分也明显不同，原因在于模型控制的低频部分随滤波通带的不同而不同.至于滤波通带下GOCE 引力梯度能够有效地反映出的阶次正好可从公式(6)计算得到，表 1列出的频率与阶次之间的关系也可从图 4中明确地表现出来，这也间接表明了公式(6)的正确性.
|
图 4 不同滤波通带下的反演结果 Fig. 4 Different results from different filtering pass band |
本文首先分析了GOCE 轨道处梯度数据的频谱特性，然后推导出球谐模型的阶次与最大截止频率之间的关系，并通过GOCE 实际引力梯度数据的处理验证了该组公式的正确性.
本文也介绍了一组针对GOCE 引力梯度观测数据的新的滤波方法，即低频部分可用某个先验模型替代，而在观测频段内保留GOCE 的观测信息.从数值实现来看，这样的滤波方法简单有效，并能有效保留观测频段内的信号.基于文中建立的滤波方法并结合引力梯度不变量使用调和分析对引力场位系数进行了恢复，最终结果与欧空局公布的GOCE引力场模型对比表明，本文建立的滤波方法能够有效地恢复引力场的位系数.
此外需要说明的是，本文建立的滤波方式同样可用于最小二乘法来解算引力场位系数.事实上，方程(7)本质上也是引力场位系数的法方程组，因此可以用矩阵计算来代替积分恢复引力场位系数.
与Schuh等提供的方法相比，本文中先验模型的使用主要是提供低频信号，Schuh 等提供的方法中则是用来提供先验误差功率谱；向前向后滤波充当的作用则相当于对方程两边的同步滤波，其作用均是消除相位漂移.但是Schuh等提供的方法很难用于调和分析，而本文的滤波方法同时适用于最小二乘和调和分析，因此通用性更好.
致谢感谢欧空局提供GOCE 数据!
| [1] | Balmino G, Rummel R, Visser P, et al. Gravity Field and Steady-State Ocean Circulation Mission. The Four Candidate Earth Explorer Core Missions. ESA Publications Division, 1999. |
| [2] | Albertella A, Migliaccio F, Sanso F. GOCE: The Earth gravity field by space gradiometry. Celestial Mechanics & Dynamical Astronomy , 2002, 83(1-4): 1-15. |
| [3] | Schuh W D. The processing of band-limited measurements: Filtering techniques in the least squares context and in the presence of data gaps. Space Science Reviews , 2003, 108(1-2): 67-78. |
| [4] | Haagmans R. GOCE Data and Gravity Field Model Filter Comparison. ESA: Stuttgart, 2007. |
| [5] | Reguzzoni M, Tselfes N. Optimal multi-step collocation: application to the space-wise approach for GOCE data analysis. Journal of Geodesy , 2009, 83(1): 13-29. DOI:10.1007/s00190-008-0225-x |
| [6] | Brockmann J M, Kargoll B, Krasbutter I, et al. GOCE data analysis: from calibrated measurements to the global earth gravity field. Advanced Technologies in Earth Sciences , 2010, 3: 213-229. |
| [7] | Bouman J, Koop R, Tscherning C C, et al. Calibration of GOCE SGG data using high low SST, terrestrial gravity data and global gravity field models. Journal of Geodesy , 2004, 78(1-2): 124-137. |
| [8] | Kern M, Preimesberger T, Allesch M, et al. Outlier detection algorithms and their performance in GOCE gravity field processing. Journal of Geodesy , 2005, 78(9): 509-519. DOI:10.1007/s00190-004-0419-9 |
| [9] | Visser P N A M. GOCE gradiometer validation by GPS. Advances in Space Research , 2007, 39(10): 1630-1637. DOI:10.1016/j.asr.2006.09.014 |
| [10] | Rispens S, Bouman J. Calibrating the GOCE accelerations with star sensor data and a global gravity field model. Journal of Geodesy , 2009, 83(8): 737-749. DOI:10.1007/s00190-008-0290-1 |
| [11] | Visser P N A M. GOCE gradiometer: estimation of biases and scale factors of all six individual accelerometers by precise orbit determination. Journal of Geodesy , 2009, 83(1): 69-85. DOI:10.1007/s00190-008-0235-8 |
| [12] | 罗志才, 吴云龙, 钟波, 等. GOCE卫星重力梯度测量数据的预处理. 武汉大学学报(信息科学版) , 2009(10): 1163–1167. Luo Z C, Wu Y L, Zhong B, et al. Pre-processing of the GOCE satellite gravity gradiometry data. Geomatics and Information Science of Wuhan University (in Chinese) , 2009(10): 1163-1167. |
| [13] | 徐新禹, 王正涛, 邹贤才, 等. GOCE卫星重力梯度测量误差分析及其模拟研究. 大地测量与地球动力学 , 2010, 30(2): 71–75. Xu X Y, Wang Z T, Zou X C, et al. Research on analysis and simulation of gravity gradiometry error of GOCE satellite. Journal of Geodesy and Geodynamics (in Chinese) , 2010, 30(2): 71-75. |
| [14] | 吴云龙, 罗志才, 李辉, 等. 基于小波收缩阈值降噪的卫星重力梯度数据粗差探测方法. 大地测量与地球动力学 , 2010, 30(4): 55–58. Wu Y L, Luo Z C, Li H, et al. Outlier detection algorithm for satellite gravity gradient data by using wavelet shrinkage denoising. Journal of Geodesy and Geodynamics (in Chinese) , 2010, 30(4): 55-58. |
| [15] | 钟波. 基于GOCE卫星重力测量技术确定地球重力场的研究. 武汉: 武汉大学, 2010. Zhong B. Study on the determination of the Earth's gravity field from satellite gravimetry mission GOCE (in Chinese). Wuhan: Wuhan University, 2010. http://www.oalib.com/references/18987354 |
| [16] | Pail R, Goiginger H, Schun W D, et al. Combined satellite gravity field model GOCO01S derived from GOCE and GRACE. Geophys. Res. Lett. , 2010, 37. DOI:10.1029/2010GL044906 |
| [17] | Milani A, Rossi A, Villani D. A timewise kinematic method for satellite gradiometry: GOCE simulations. Earth Moon and Planets , 2005, 97(1-2): 37-68. |
| [18] | Reguzzoni M. From the time-wise to space-wise GOCE observables. Advances in Geosciences , 2003, 1: 137-142. DOI:10.5194/adgeo-1-137-2003 |
| [19] | Sansò F, Migliaccio F, Reguzzoni M. Space-wise approach to satellite gravity field determination in the presence of coloured noise. Journal of Geodesy , 2004, 78(4-5): 304-313. DOI:10.1007/s00190-004-0396-z |
| [20] | Yu J H, Zhao D M. The gravitational gradient tensor's invariants and the related boundary conditions. Science China Earth Sciences , 2010, 53(5): 781-790. DOI:10.1007/s11430-010-0014-2 |
| [21] | Hofmann W B, Moritz H. Physical Geodesy. Wien: Springer Wien New York, 2005. http://www.oalib.com/references/19235317 |
| [22] | Pavlis N K, Holmes S A, Kenyon S C, et al. An Earth Gravitational Model to Degree 2160: EGM2008, presented at the 2008 General Assembly of the European Geosciences Union. Vienna, Austria, April 13-18, 2008. |
| [23] | Kaula W M. Theory of Satellite Geodesy. London: Blaisdell Publishing Company , 1966: 3-59. |
| [24] | Sneeuw N. A Semi-Analytical Approach to Gravity Field Analysis from Satellite Observations. Munich: Technique University Munich, 2000. |
| [25] | 于锦海, 万晓云. 引力梯度归算的模拟计算. 地球物理学报 , 2011, 54(5): 1182–1186. Yu J H, Wan X Y. Reduction for gradiometry and corresponding imitation. Chinese J. Geophys. (in Chinese) , 2011, 54(5): 1182-1186. |
| [26] | Yu J H, Wan X Y. The frequency analysis of gravity gradients and the methods of filtering processing. Proceedings of the 4th International GOCE User Workshop. Munich, Germany, 2011. |
| [27] | 维纳 K 英格尔, 约翰 G 普罗克斯著. 数字信号处理(MATLAB版). 刘树棠译. 西安: 西安交通大学出版社, 2009: 172-181. Vinay K I, John G P. Digital Signal Processing using MATLAB (in Chinese). Translated by Liu S T. Xi'an: Xi'an Jiaotong University Press, 2009: 172-181. |
| [28] | 邢园丁, 邹虹, 杨德猛, 等. 零相位滤波器在基于激光多普勒效应的动态位移信号处理中的应用. 电子测量技术 , 2009, 32(6): 48–50. Xing Y D, Zou H, Yang D M, et al. Application of zero-phase filter in dynamic displacement signal processing based on laser doppler effect. Electronic Measurement Technology (in Chinese) , 2009, 32(6): 48-50. |
| [29] | 吴国乔, 王兆华. 基于全相位的零相位数字滤波器的设计方法. 电子与信息学报 , 2007, 29(3): 574–577. Wu G Q, Wang Z H. Design method of digital filter with zero-phase based on all phase. Journal of Electronics & Information Technology (in Chinese) , 2007, 29(3): 574-577. |
| [30] | 纪跃波, 秦树人, 汤宝平. 零相位数字滤波器. 重庆大学学报(自然科学版) , 2000, 23(6): 4–7. Ji Y B, Qin S R, Tang B P. Digital filtering with zero phase error. Journal of Chongqing University (Natural Science Edition) (in Chinese) , 2000, 23(6): 4-7. |
| [31] | Pail R, Bruinsma S, Migliaccio F, et al. First GOCE gravity field models derived by three different approaches. Journal of Geodesy , 2011, 85(11): 819-843. DOI:10.1007/s00190-011-0467-x |