2. 同济大学测绘与地理信息学院, 上海 200092;
3. 波恩大学大地测量与地理信息学院, 波恩 53121
2. College of Surveying and Geo-Informatics, Tongji University, Shanghai 200092, China;
3. Institute of Geodesy and Geo-information, University of Bonn, Bonn 53121, Germany
2009年，欧空局(European Space Agency，ESA)发射了GOCE(gravity field and steady-state ocean circulation explorer)卫星，其主要科学目标是确定约100 km分辨率、1~2 cm精度的大地水准面，可为固体地球物理学、海洋学、大地测量学和冰川学等科学研究提供高精度的中高阶静态地球重力场模型[1]。GOCE卫星在任务期间共获取了约42个月的有效重力梯度观测值，由于重力梯度观测值包含了时变重力场变化信号，在高阶静态重力场反演中需考虑时变重力场变化改正的影响。ESA在Level 2级别的GOCE卫星重力梯度数据中给出了ESA计算的时变重力场变化改正值[2]，目前国内外绝大多数科研机构均直接利用GOCE卫星Level 2重力梯度数据进行重力场反演[3-5]，但时变重力场变化改正是实现由GOCE卫星Level 1b级别数据到Level 2级别数据转换的重要过程，也是完整GOCE卫星数据处理过程中的一个重要环节，而ESA有关文档中并未给出详细的数据处理过程，自主实现由GOCE卫星Level 1b重力梯度数据直接进行重力场反演具有重要的科学意义，可为我国独立开展相关GOCE类卫星任务的数据处理研究提供重要参考。
GOCE卫星重力梯度观测值的时变重力场变化改正主要包括直接潮汐改正、固体潮汐改正、海洋潮汐改正、极潮改正和非潮汐质量变化改正等[6-7]。鉴于GOCE卫星重力梯度观测值的时变重力场变化改正的重要性，国内外已有学者对重力梯度的时变重力场变化改正展开研究[7-16]，分别从功率谱密度特性[13-15]以及时变重力场变化的高频信号源[9-10]等方面分析了时变重力场变化改正对GOCE卫星重力梯度观测值的影响。但目前暂未见从重力场反演结果层面分析时变重力场变化改正及标准和背景模型更新对高阶静态重力场反演的影响等相关研究成果。基于上述背景，本文采用直接法进行GOCE卫星重力反演，并进行有关分析，首先选用和ESA一致或同类的标准和背景模型计算出时变重力场变化改正值并与ESA结果进行比较，验证本文计算的正确性；然后采用新标准和背景模型对GOCE卫星重力梯度观测值进行时变重力场变化改正，自主实现了由Level 1b重力梯度数据直接进行重力反演；最后通过3种时变重力场变化改正方案分别反演220阶次的重力场模型，从大地水准面阶误差和不同方案模型间全球大地水准面差异等方面来分析时变重力场变化改正对高阶静态重力场反演结果的影响，并利用实测大地水准面模型对结果进行检核。
1 原理与方法GOCE卫星的主要任务是反演高精度的静态地球重力场，但重力梯度观测值包含了时变重力场变化信号，需要考虑对其进行时变重力场变化改正，主要包括：直接潮汐改正、固体潮汐改正、海洋潮汐改正、极潮改正和非潮汐质量变化改正等，其中，直接潮汐改正先计算引潮位，再将其转换为重力梯度改正值，而固体潮汐改正、海洋潮汐改正、极潮改正、非潮汐质量变化改正是先确定对应的位系数改正数，然后将其转换为重力梯度改正值。
1.1 直接潮汐改正直接潮汐是由日月及其他天体对卫星的引潮力引起的，对应的引潮力位Vj可表示为[17]
(1)
式中，G为万有引力常数；M⊕为地球质量；R为地球平均半径；(θ, λ, r)分别为计算点的地心余纬、地心经度和地心向径；Pnm为完全规格化的n阶m次缔合勒让德函数；Mj为天体j的质量；(θj, λj, rj)分别为天体j的地心余纬、地心经度和地心向径；Nmax为计算的最大阶数，在此项改正中本文取至10阶，可采用DE430等JPL星历文件计算天体位置。将引潮力位Vj对(θ, λ, r)分别求一阶偏导数和二阶偏导数后，再利用旋转矩阵进行转换即可得到梯度仪坐标系下的重力梯度改正值。
1.2 固体潮汐改正固体潮汐引起的重力位系数改正可分为两步计算[18]。
第1步：计算固体潮汐引起的2阶和3阶重力位系数改正，以及2阶对4阶重力位系数的改正，公式为
(2)
(3)
第2步：计算2阶项的日潮、半日潮及长周期潮改正，其中计算日潮改正的公式为
(4)
2阶项的潮汐校正包含48个日潮、两个半日潮以及21个长周期潮(其他具体计算公式可参考文献[18])，f为各潮波分量, Ampf(ip/op)和其他相关参数含义可通过文献[18]获得。
1.3 海洋潮汐改正由海洋潮汐引起的重力位系数改正数可用式(5)计算[18]
(5)
式中，Cf, nm±、Sf, nm±为分潮波f的海潮系数，θf为潮波分量f的相位角。常用的海潮模型有FES2004、EOT11a和FES2014b等。
1.4 极潮改正极潮改正中主要考虑固体极潮，其引起的重力位系数改正数中，主要考虑ΔC21和ΔS21的影响[6, 17]，其计算公式为[18]
(6)
式中，m1、m2参数含义及其具体计算过程可参考文献[18]。
1.5 非潮汐质量变化改正非潮汐质量变化改正包括大气与海洋非潮汐质量变化改正、陆地水质量变化改正和冰雪质量变化改正等地表质量季节性变化的改正。
1.5.1 大气与海洋非潮汐质量变化改正ESA提供的大气质量变化改正采用ECMWF(European Centre for Medium-range Weather Forecast)模型计算，其海洋非潮汐质量变化改正采用OMCT(ocean model for circulation and tides)模型数据计算。由于GFZ研发的AOD1B(Atmosphere and Ocean De-aliasing level-1B)模型给出了大气和海洋非潮汐变化模型，因此，本文直接采用AOD1B模型中提供的有关重力位系数的改正数进行计算。
1.5.2 地表质量季节性变化改正地表质量季节性变化改正主要为陆地水质量变化改正和冰雪质量变化改正等，而GRACE月时变重力场模型已扣除了各种潮汐及大气和海洋非潮汐质量变化影响，剩余的时变重力场变化信号主要来自地表质量季节性变化信号，因此，可以采用GRACE月时变重力场模型计算地表质量季节性变化影响的改正值[6, 10]。
2 计算结果验证与更新比较为验证本文时变重力场变化改正的正确性，并分析标准和背景模型更新前后时变重力场变化改正结果的差异，本文选择了2010年6月1日的GOCE卫星实测数据，分别计算出各项潮汐改正值以及非潮汐质量变化改正值，并进行有关比较与分析，所采用的标准和背景模型信息见表 1，其中"ESA"表示ESA标准，"Tongji"表示本文采用的新标准。后文计算结果中，"ESA"表示ESA发布结果；"Tongji"表示本文按照新标准计算结果；"Tongji(ESA)"表示本文按照ESA标准计算结果，其中计算大气与海洋非潮汐质量变化改正时，采用了AOD1B RL05模型来代替ECMWF和OMCT模型，其他与ESA标准一致。
2.1 时变重力场变化改正结果验证
ESA和Tongji(ESA)比较结果见图 1—图 5，两者差值标准差的统计结果见表 2。综合图 1—图 4以及表 2可得，当采用的标准和背景模型与ESA一致时，本文的计算结果与ESA的结果是吻合的；本文计算值与ESA结果差值的标准差比对应潮汐改正值的最大量级低1~3个量级，说明本文的潮汐改正结果是正确的。图 5结果显示，本文非潮汐质量变化计算值与ESA发布结果在数值上稍有差异，产生差异的主要原因是所使用的背景模型不完全一致(本文计算大气与海洋非潮汐质量变化采用的背景模型为AOD1B RL05模型)，但两者在趋势和量级上具有较好的一致性，说明本文非潮汐质量变化改正值计算结果应该是可靠的。综合看，本文时变重力场变化改正的结果是正确的，在所有改正项中，海潮改正的量级最大，其次为非潮汐质量变化改正，极潮改正的量级最小。
|
| 图 1 直接潮汐改正值与ESA结果比较 Fig. 1 Comparisons between direct tide corrections and ESA results |
|
| 图 2 固体潮汐改正值与ESA结果比较 Fig. 2 Comparisons between solid tide corrections and ESA results |
|
| 图 3 海洋潮汐改正值与ESA结果比较 Fig. 3 Comparisons between ocean tide corrections and ESA results |
|
| 图 4 极潮改正值与ESA结果比较 Fig. 4 Comparisons between pole tide corrections and ESA results |
|
| 图 5 非潮汐改正值与ESA值的比较 Fig. 5 Comparisons between non-tidal corrections and ESA results |
| 改正项 | Vxx | Vyy | Vzz | Vxy | Vxz | Vyz |
| 直接潮汐 | 1.168E-04 | 1.869E-05 | 1.124E-04 | 8.645E-05 | 1.140E-04 | 7.712E-05 |
| 固体潮汐 | 1.912E-03 | 1.446E-03 | 3.341E-03 | 2.335E-03 | 2.227E-03 | 9.592E-03 |
| 海洋潮汐 | 7.180E-03 | 6.119E-03 | 1.181E-02 | 3.002E-03 | 8.129E-03 | 6.549E-03 |
| 极潮 | 1.141E-05 | 8.699E-06 | 1.986E-05 | 3.580E-06 | 1.305E-05 | 1.467E-05 |
| 非潮汐质量变化 | 1.113E-02 | 1.106E-02 | 1.908E-02 | 5.433E-03 | 1.293E-02 | 1.260E-02 |
2.2 ESA标准和背景模型更新比较
由于ESA计算重力梯度观测值的时变重力场变化改正所采用的标准和背景模型较老，理论上，新标准和背景模型能更好地模型化时变重力场变化信号，因此本节探讨标准和背景模型更新前后时变重力场变化改正值的差异(新标准和背景模型采用表 1中"Tongji"项)，计算结果比较见图 6。由图 6可知，更新的标准和背景模型会影响时变重力场变化改正的结果，但更新前后两者的时变重力场变化改正值趋势一致，大小稍有差异但差异不大，总量在同一量级，下文将详细分析标准和背景模型更新前后对高阶静态重力场反演的影响。
|
| 图 6 Tongji和ESA的时变重力场变化改正值比较 Fig. 6 Comparisons of temporal gravity field variations corrections between Tongji and ESA |
3 时变重力场变化改正对高阶静态重力场反演的影响分析
本节设计了3种方案来分析时变重力场变化改正对利用GOCE卫星重力梯度数据反演高阶静态重力场的影响，分别为：①无改正方案(Case 1)，不进行任何时变重力场变化改正；②ESA方案(Case 2)，采用ESA发布的时变重力场变化改正值；③Tongji方案(Case 3)，采用新标准和背景模型计算时变重力场变化改正值，3种方案差异仅在于时变重力场变化改正值，其他数据处理方法完全一致。
考虑到低轨道GOCE卫星对重力场信号更敏感，本文采用GOCE卫星2012年10月—2013年4月期间新版的Level 1b重力梯度观测数据[23]，按照上述3种方案进行高阶静态重力场反演。首先选择GIF48模型截断至300阶次作为本次反演的参考模型，利用频带范围为5~100 mHz的IIR带通滤波器对扰动重力梯度值进行带通滤波，采用直接法构建反演220阶次重力场模型的法方程；然后与180阶次的Tongji-Grace02s模型[24]法方程进行融合，反演220阶次的无约束重力场模型。在数据融合过程中单位权中误差取GRACE卫星的轨道精度，GRACE KBR观测值精度取验后统计精度，而重力梯度分量的精度根据滤波之后的扰动重力梯度值进行统计确定。3种方案解算模型和Tongji-Grace02s模型的大地水准面阶误差如图 7所示，比对模型为EIGEN-6C4模型。
|
| 图 7 各方案模型与EIGEN-6C4模型比较大地水准面阶误差 Fig. 7 Geoid degree error of different models with respect to EIGEN-6C4 |
图 7结果显示，3种方案解算模型的大地水准面阶误差几乎一致，不能明显地分析时变重力场变化改正对静态重力场反演的影响。融合了GOCE卫星重力梯度数据后，模型大地水准面阶误差相对于Tongji-Grace02s模型在72阶次后逐渐有了改善，在108阶次后其改善效果更为明显，说明相比GRACE卫星，GOCE卫星重力梯度数据对中高阶重力位系数的贡献更大，因此在中高阶的位系数精度改善效果会比较明显。为了更清晰地分析时变重力场变化对静态重力场反演的影响，本文计算了全球1°×1°的大地水准面(220阶次)，其比较结果见图 8(a)、图 8(b)和图 8(c)，同时图 8(d)给出了100~220阶次不同方案反演模型所对应大地水准面差异的统计结果。
|
| 图 8 不同模型间的全球大地水准面比较 Fig. 8 Comparison of global geoid height between different models |
由图 8(a)和图 8(b)可知，时变重力场变化改正对重力场反演是有影响的，其在局部区域的大地水准面差异的绝对值最大达到1.4 cm，差异主要体现在南北两极、阿根廷东部沿海和欧洲等地区，这些地区也是该时间段时变重力场变化信号较为明显的区域，如果采用的数据观测时间段不同，时变重力场变化影响会有所差异。但由图 8(c)可知，采用不同的标准和背景模型同样也会影响高阶静态重力场反演结果，两者在南北两极、亚马孙河口附近和欧洲等区域中稍有差异。由图 8(d)可知，从全球大地水准面差异看，在220阶次，时变重力场变化改正对静态重力场反演的整体影响超过1 mm，影响主要体现在中高阶次位系数，低于120阶次的影响很小。综合图 8可知：GOCE卫星重力梯度观测值的时变重力场变化改正是需要的，且不同标准和背景模型对静态重力场反演的影响稍有差异。
为进一步检核时变重力场变化改正对利用GOCE卫星重力梯度数据反演高阶静态重力场的影响，同时考虑到南北两极可能会受到极区空白的影响，而极区空白对远离两极的区域大地水准面的精度影响较小[25]，因此本文利用收集到的欧洲区域(图 8中时变重力场影响较大的区域)约1 cm精度的France水准面模型RAF18(图 8(b)中的A区域)和Great Britain水准面模型OSGM15(图 8(b)中的B区域)对本文3种反演方案获得的重力场模型进行精度评定，精度统计结果见表 3(其中截断阶次后采用EIGEN-6C4模型填充至2160阶次)。
| 似大地水准面模型 | 阶数 | Case 1 | Case 2 | Case 3 |
| France(RAF18)(42°N—51.5°N，5.5°W—8.5°E)1.5'×2' | 180 | 0.139 55 | 0.139 56 | 0.139 54 |
| 190 | 0.141 19 | 0.141 22 | 0.141 17 | |
| 200 | 0.144 33 | 0.144 36 | 0.144 28 | |
| 210 | 0.150 61 | 0.150 45 | 0.150 20 | |
| 220 | 0.154 13 | 0.154 01 | 0.153 75 | |
| Great Britain(OSGM15)(51°N—56°N，11.5°W—5°W)0.8'×1.2' | 180 | 0.131 34 | 0.131 33 | 0.131 31 |
| 190 | 0.130 30 | 0.130 25 | 0.130 22 | |
| 200 | 0.138 30 | 0.138 13 | 0.138 00 | |
| 210 | 0.144 18 | 0.144 09 | 0.143 48 | |
| 220 | 0.162 39 | 0.162 40 | 0.161 51 |
表 3结果表明，虽然不同方案模型间差异不大，但综合看，进行时变重力场变化改正(Case 2和Case 3)方案比无改正方案(Case 1)解算的模型精度稍好，说明重力梯度观测值的时变重力场变化改正是有效的，同时Tongji方案(Case 3)精度稍优于ESA方案(Case 2)，说明背景模型和标准的更新更有利于扣除时变重力场变化信号，从而提高高阶静态重力场反演的精度。
4 结束语本文研究了GOCE卫星重力梯度观测值的时变重力场变化改正方法，更新了ESA标准和背景模型，并从重力场反演角度分析了时变重力场变化改正对高阶静态重力场反演的影响，实现了从GOCE卫星Level 1b重力梯度数据直接进行重力场反演，主要结论如下。
(1) 在时变重力场变化改正中，海洋潮汐改正值的量级最大，其次为非潮汐质量变化改正，极潮改正值的量级最小。利用3种时变重力场变化改正方案进行高阶静态重力场模型反演，并进行分析与比较，结果表明：从全球大地水准面比较可看出时变重力场变化改正对重力场反演结果是有影响的，其影响在局部区域的大地水准面差异绝对值最大达到1.4 cm；在全球范围内对重力场反演的整体影响超过1 mm，因此利用GOCE卫星重力梯度观测值反演高阶静态重力场时需要扣除时变重力场变化的影响。
(2) 利用France水准面模型RAF18和Great Britain水准面模型OSGM15检核3种方案模型的精度，结果表明进行时变重力场变化改正方案比无改正方案解算的模型精度稍好，且更新后的标准和背景模型更有利于扣除时变重力场变化的影响。
| [1] |
DRINKWATER M R, FLOBERGHAGEN R, HAAGMANS R, et al. GOCE: ESA's first earth explorer core mission[J]. Space Science Reviews, 2003, 108(1): 419-432. DOI:10.1023/A:1026104216284 |
| [2] |
ESA. GOCE HPF: GOCE level 2 product data handbook[R]. Technical Note, GO-MA-HPF-GS-0110. [S. l. ]: ESA, 2014.
|
| [3] |
于锦海, 万晓云. 利用引力梯度不变量解算的GOCE引力场模型[J]. 中国科学(地球科学), 2012, 42(9): 1450-1458. YU Jinhai, WAN Xiaoyun. Recovery of the gravity field from GOCE data by using the invariants of gradient tensor[J]. Science China Earth Sciences, 2012, 42(9): 1450-1458. |
| [4] |
XU Xinyu, ZHAO Yongqi, REUBELT T, et al. A GOCE only gravity model GOSG01S and the validation of GOCE related satellite gravity models[J]. Geodesy and Geodynamics, 2017, 8(4): 260-272. DOI:10.1016/j.geog.2017.03.013 |
| [5] |
MARCHENKO A N, MARCHENKO D A, LOPUSHANSKY A N. Gravity field models derived from the second degree radial derivatives of the GOCE mission: a case study[J]. Annals of Geophysics, 2017, 59(6): s0649-s0659. |
| [6] |
GRUBER T, ABRIKOSOV O, HUGENTOBLER U. Doc No: GO-TN-HPF-GS-0111 GOCE standards[S]. 2014.
|
| [7] |
BOUMAN J, RISPENS S, GRUBER T, et al. Preprocessing of gravity gradients at the GOCE high-level processing facility[J]. Journal of Geodesy, 2009, 83(7): 659-678. DOI:10.1007/s00190-008-0279-9 |
| [8] |
BOUMAN J, FIOROT S, FUCHS M, et al. GOCE gravitational gradients along the orbit[J]. Journal of Geodesy, 2011, 85(11): 791-805. DOI:10.1007/s00190-011-0464-0 |
| [9] |
HAN S C, SHUM C K, DITMAR P, et al. Aliasing effect of high-frequency mass variations on GOCE recovery of the earth's gravity field[J]. Journal of Geodynamics, 2006, 41(1-3): 69-76. DOI:10.1016/j.jog.2005.08.005 |
| [10] |
JARECKI F, MÜLLER J, PETROVIC S, et al. Temporal gravity variations in GOCE gradiometric data[C]//Proceedings of 2005 IAG International Symposium Porto. Berlin, Heidelberg: Springer, 2005, 129: 333-338.
|
| [11] |
邓文彬, 许闯, 阿力甫·努尔买买提. 卫星重力梯度观测数据的时变信号影响分析[J]. 武汉大学学报(信息科学版), 2016, 41(1): 72-78. DENG Wenbin, XU Chuang, NURMMAMAT A. Effect of time-varying gravity signal on satellite gravity gradient observations[J]. Geomatics and Information Science of Wuhan University, 2016, 41(1): 72-78. |
| [12] |
吴云龙. GOCE卫星重力梯度测量数据的预处理研究[D]. 武汉: 武汉大学, 2010. WU Yunlong. Study on the pre-processing of GOCE satellite gravity gradiometry data[D]. Wuhan: Wuhan University, 2010. |
| [13] |
周睿. GOCE卫星SGG数据处理及空域法恢复地球重力场模型研究[D]. 郑州: 信息工程大学, 2013. ZHOU Rui. Research on GOCE SGG data processing and the space-wise approach to recover gravity field model[D]. Zhengzhou: Information Engineering University, 2013. |
| [14] |
YI Weiyong. The earth's gravity field from GOCE[D]. Germany: Technische Universität Mänchen, 2012.
|
| [15] |
QU Qingliang, CHANG Xiaotao, YU Shengwen, et al. Calibration for GOCE gradiometer data based on the prior gravity models[J]. Journal of Geodesy and Geoinformation Science, 2019, 2(4): 21-30. |
| [16] |
WU Hu. Gravity field recovery from GOCE observations[D]. Germany: Fachrichtung Geodäsie und Geoinformatik der Leibniz Universität Hannover, 2016.
|
| [17] |
SÜNKEL H. From Eötvös to Mgal, final report[R]. Noordwijk: European Space Agency, 2000.
|
| [18] |
PETIT G, LUZUM B. IERS Conventions. IERS 2010[R]. Frankfurt: Verlag des Bundesamts für Kartographie und Geodäsie, 2010.
|
| [19] |
MCCARTHY D D, PETIT G. IERS conventions (2003)[M]. Frankfurt: IERS Central Bureau, 2004.
|
| [20] |
CARRERE L, LYARD F, CANCET M, et al. FES 2014, A new tidal model on the global ocean with enhanced accuracy in shallow seas and in the Arctic region[C]//EGU general Assembly Conference Abstracts. Vienna, Austria: EGU, 2015.
|
| [21] |
DOBSLAW H, BERGMANN-WOLF I, DILL R, et al. A new high-resolution model of non-tidal atmosphere and ocean mass variability for de-aliasing of satellite gravity observations: AOD1B RL06[J]. Geophysical Journal International, 2017, 211(1): 263-269. DOI:10.1093/gji/ggx302 |
| [22] |
CHEN Qiujie, SHEN Yunzhong, CHEN Wu, et al. An optimized short-arc approach: methodology and application to develop refined time series of Tongji-Grace2018 GRACE monthly solutions[J]. Journal of Geophysical Research: Solid Earth, 2019, 124(6): 6010-6038. DOI:10.1029/2018JB016596 |
| [23] |
SIEMES C, REXER M, SCHLICHT A, et al. GOCE gradiometer data calibration[J]. Journal of Geodesy, 2019, 93(9): 1603-1630. DOI:10.1007/s00190-019-01271-9 |
| [24] |
CHEN Qiujie, SHEN Yunzhong, FRANCIS O, et al. Tongji-Grace02s and Tongji-Grace02k: high-precision static GRACE-only global earth's gravity field models derived by refined data processing strategies[J]. Journal of Geophysical Research: Solid Earth, 2018, 123(7): 6111-6137. DOI:10.1029/2018JB015641 |
| [25] |
万晓云, 于锦海. 极地空白对GOCE引力场恢复的影响[J]. 测绘学报, 2013, 42(3): 317-322. WAN Xiaoyun, YU Jinhai. Influence of polar gaps on gravity field recovery using GOCE data[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(3): 317-322. |