引言

高精度、高分辨率地球重力场为全球或区域高程基准统一、大地水准面确定、海洋洋流确定、强地震监测等提供更加精细的地球重力场信息(许厚泽等,2012;Becker et al.,2014; Gerlach and Rummel,2013;Fuchs et al.,2013).GOCE(Gravity field and steady state Ocean Circulation Explorer)卫星于2009年3月发射,2013年11月结束任务,共获取了约42个月的对地观测数据,其主要科学目标是确定高精度的中高阶静态地球重力场,其对应的空间分辨率约为100 km,精度为1～2 cm的大地水准面和1～2 mGal的重力异常.GOCE重力卫星计划的成功实施,对提高全球重力场模型的精度和分辨率取得了显著作用(Pail et al.,2011).国内外科研机构采用GOCE卫星观测数据或联合其他卫星(如GRACE卫星等)观测数据解算了多个版本的全球重力场模型(Brockmann et al.,2014,2019;Bruinsma et al.,2014;Förste et al.,2019;Kvas et al.,2019;Wu et al.,2016;Xu et al.,2017;Lu et al.,2018;苏勇等, 2018;Zhou et al.,2019),近来国际上融合新版GOCE梯度数据和GRACE卫星数据解算了精度最高的DIR-R6等系列模型,而在我国,独立研制与国际精度相当的重力场模型对于自主重力卫星计划的实施具有重要的科学意义和现实意义.

利用GOCE卫星数据进行重力解算主要有3种方法,即直接法、时域法和空域法(Pail et al.,2011;郑伟等,2014;于锦海和万晓云,2012),欧空局(ESA)官方发布的重力场模型均由这3种方法解算,相应成果已在ICGEM (International Centre for Global Earth Models)网站(http://icgem.gfz-potsdam.de/home)发布.直接法是选定某一参考模型,直接在卫星观测轨道处建立扰动重力梯度与重力位系数观测方程,可通过联合GRACE卫星数据及LAGEOS 1/2等SLR数据反演重力场(Pail et al.,2011;Bruinsma et al.,2013,2014;Förste et al.,2019).时域法是将观测值看作时间序列,根据卫星轨道摄动理论构建观测方程,由于时域法未联合GRACE等其他卫星数据,其解算的模型也称为纯GOCE卫星重力场模型(Pail et al.,2011;Brockmann et al.,2014,2019),通过进一步融合两极地面重力数据反演了TIM-R6e模型(Zingerle et al.,2019).空域法是基于大地测量边值理论构建梯度观测值与位系数关系式,并通过多步最小二乘配置法解算位系数,解算的模型也为纯GOCE卫星重力场模型(Reguzzoni and Tselfes,2009;Pail et al.,2011;Gatti and Reguzzoni,2017).为了克服GOCE卫星极区空白及法方程病态等问题,通常采用球冠正则化、Kaula正则化等正则化方法(Förste et al.,2019).综合看,直接法解算重力位系数的难度和计算量最大,而时域法和空域法通过一些假设或近似,减小了计算量,但直接法解算的DIR_Rx系列模型精度要明显优于时域法的TIM_Rx系列和空域法的SPW_Rx系列模型.目前,由GOCE完整观测周期的数据求得的第5、6代(R5、R6)重力场模型已基本实现GOCE卫星的设计目标(Brockmann et al.,2014;Bruinsma et al.,2014;Gruber,2019),而第6代模型由于采用了新版GOCE卫星重力梯度数据(Gruber,2019;Siemes et al.,2019),模型精度优于第5代.此外,GOCO团队(Gravity Observation Combination)联合GOCE、GRACE等卫星数据也获得了与同期官方模型精度相当的GOCO01S、GOCO03s、GOCO05s以及GOCO06s等模型(Pail et al.,2010; Kvas et al.,2019),其中,GOCO06s模型还采用了Swarm等低轨卫星数据以及Starlette等SLR数据,是目前采用卫星数据最多的重力场模型;Schall等(2014)和Yi等(Yi, 2012a,b; Yi et al.,2013)采用GOCE卫星数据分别反演了240阶次的ITG-Goce02模型和230阶次的JYY_GOCE04S模型,JYY_GOCE04S模型在极空白区利用ITG-Grace2010s模型信息进行了约束,而ITG-Goce02模型自5阶起采用kaula正则化进行约束.在国内,Zhou等(2019)、Lu等(2018)、苏勇等(2018)、梁建青(2016)和Xu等(2017)等采用GOCE卫星数据(或者联合GRACE卫星数据)分别反演了210阶次的HUST-GOGRA2018s模型、240阶次的IGGT_R1模型、210阶次的SWJTU-GOGR01S模型、250阶次的Tongji-GRACE-GOCE01模型以及220阶次的GOSG01S模型,但这些模型并未使用新版的GOCE卫星重力梯度数据,模型精度也稍逊于GOCE卫星第6代模型.罗志才等(2015)、于锦海和万晓云(2012)、徐天河和贺凯飞(2009)、钟波(2010)、郑伟等(2011)、刘晓刚(2011)和黄强(2014)等对GOCE卫星数据的重力解算所涉及的预处理方法、解算理论与方法等进行了相关研究和探讨,所有这些都促进了我国GOCE卫星重力反演理论的发展.

在国内,融合GRACE和GOCE卫星数据自主研制可达国际主流模型精度的重力场模型迫在眉睫,鉴于直接法具有解算精度高,且便于联合多种类型的观测数据进行重力场模型解算等优点,并考虑到GOCE卫星任务目标(确定分辨率约为100 km,大地水准面精度为1～2 cm的重力场),本文选择直接法作为重力反演方法,反演220阶次的重力场.基于直接法和IIR带通滤波技术,利用约24个月的GOCE卫星重力梯度数据构建了反演220阶次重力场模型的法方程,并与由13.5年GRACE卫星数据构建的180阶次的Tongji-Grace02s模型的法方程融合,解算了220阶次的Tongji-GOGR2019S重力场模型,最后采用多源数据对Tongji-GOGR2019S模型精度进行全面检核.

1、GOCE和GRACE卫星观测模型与联合反演

1.1基于GOCE卫星重力梯度观测数据的重力反演

在球坐标系中,地球外部空间任意一点的引力位可表示为(Pail et al.,2011)

V(θ,λ,r)=GMR∑n=0Nmax(Rr)n+1∑m=0n(C¯¯¯n mcosm λ   +S¯¯n msinm λ)P¯¯¯n m(cosθ),         (1)

其中,GM为引力常数与地球质量的乘积,R为地球平均半径,θ为地心余纬,λ为地心经度,r为地心向径,(C¯¯¯n m,S¯¯n m)为完全规格化的重力位系数,P¯¯¯n m(cosθ)为完全规格化的n阶m次缔合勒让德函数.

在梯度仪坐标系下(GRF),用梯度观测值进行重力反演的公式可写为

其中,RGRFLNOF是局部指北坐标系到重力梯度仪坐标系的转换矩阵,δu为重力位系数的改正数向量,Vij(i,j=x,y,z)为梯度仪坐标系下的经过时变重力场变化改正后的新版梯度观测值,V0ij为由参考重力位系数近似值u0计算得到的梯度值的近似值(LNOF坐标系),∂Vij∂u为重力梯度张量对重力位系数的偏导数(LNOF坐标系),E为梯度残差值矩阵.对(2)式中满足精度要求的4个重力梯度观测量(Vxx、Vyy、Vzz和Vxz)表示成向量,并将各历元的误差方程合并在一起,表示成如下向量形式:

vSGG=ASGGδu−ySGG,         (3)

其中,ASGG为各梯度观测量对重力位系数的偏导数组成的设计矩阵,ySGG为梯度观测值与重力位系数近似值u0计算的梯度值之差向量,vSGG为残差向量.

由于GOCE重力梯度观测值信号主要集中在5～100 mHz测量带宽内,为了保留频带内信号,并尽量减小频带外噪声对重力反演的影响,本文采用IIR数字带通滤波器对观测方程两边同步滤波,其中滤波的带宽范围可根据实际反演阶次来确定,则利用重力梯度数据解算重力场位系数的误差方程可写为

F{vSGG}=F{ASGG}δu−F{ySGG},         (4)

其中,F{}为IIR滤波算子.

1.2基于GOCE和GRACE观测数据的联合重力反演

由于GOCE观测数据在南北两极存在约±6.7°的极区空白、梯度观测值的低频噪声大以及带通滤波特性等原因,导致仅利用GOCE重力梯度直接反演重力场过程严重病态,严重影响了重力位系数的解算精度,而GRACE数据反演中长波重力位系数的精度较高,联合GOCE和GRACE可有效克服仅利用GOCE重力梯度数据进行重力反演的弊端,两类卫星观测数据联合重力反演的法方程可表示为

(σ20σ2SGGNSGG+σ20σ2KBRNKBR+σ20σ2OrbNOrb)δu=  (σ20σ2SGGwSGG+σ20σ2KBRwKBR+σ20σ2OrbwOrb),         (5)

其中,Ni、wi、σi(i=SGG,KBR,Orb)分别为GOCE梯度观测值、GRACE星间距离变率观测值和轨道观测值所对应的系数矩阵、常数项和精度,σ0为单位权中误差,本文取单位权中误差为GRACE卫星轨道精度,即σ0=σOrb为0.02 m,与Tongji-Grace02s模型单位权一致,梯度分量精度采用相对于参考模型的统计结果,与GRACE卫星有关的法方程等信息由Tongji-Grace02s重力场模型对应的法方程提供(Chen et al.,2018).

2、Tongji-GOGR2019S重力场模型解算

2.1GOCE梯度观测值滤波策略

本文采用GOCE卫星的Vxx、Vyy、Vzz和Vxz 4个重力梯度分量观测值反演重力场,考虑到GOCE重力梯度观测值的噪声特性和实际反演的模型阶数,本文采用带宽为5～41 mHz的IIR带通滤波器对扰动重力梯度观测值进行滤波,截断频率41 mHz大致对应能反演重力场模型的最高阶次为220(万晓云等,2012;Xu et al.,2017).由于不同阶的IIR滤波器,其滤波效果会有一定的差异.图1给出了利用IIR带通滤波器对2013年9月4个扰动梯度分量观测值进行滤波的结果,滤波的阶数分别取8、10和12阶,图1结果显示,IIR带通滤波可以很好地保留了各扰动重力梯度观测值在给定频带内的信号,抑制了各重力梯度观测值在频带外的噪声.由于计算效率随滤波阶数增加而降低,综合考虑滤波效果和计算效率,本文选择8阶IIR带通滤波器作为Tongji-GOGR2019S模型反演滤波器.

2.2Tongji-GOGR2019S重力场模型解算

图1不同阶IIR带通滤波器滤波结果比较

由于GOCE梯度观测值中不可避免的含有粗差,为了减小粗差对重力反演的影响,本文在所选滤波器下采用阀值法和基于移动窗口的阀值法探测各重力梯度分量观测值粗差(徐天河和贺凯飞,2009;Yi,2012a;梁建青,2016).在剔除GOCE梯度观测值粗差后,本文基于直接法和8阶IIR带通滤波技术,利用GOCE卫星2011年2月至2013年10月共约24个月经过时变重力场变化改正后的Level 1b重力梯度观测值,并联合13.5年的GRACE卫星数据(180阶次的Tongji-Grace02s模型的法方程)反演了220阶次的Tongji-GOGR2019S重力场模型,其中,Tongji-Grace02s模型是基于改进的短弧积分法,采用多项式对GRACE卫星加速度计的尺度和偏差进行校正,并顾及加速度计和姿态测量误差后解算的,在同类模型中,精度较高(Chen et al.,2016,2018).GOCE重力场模型解算所采取的主要措施如下:(1)按照重力梯度观测值连续性将整个观测数据序列划分为45个弧段,每个弧段连续滤波,各扰动重力梯度分量法方程加权累加;(2)4个重力梯度分量观测值的精度是采用GIF48为参考模型截断至300阶次统计获得.Tongji-GOGR2019S模型计算信息见表1.

Tongji-GOGR2019S等模型与EIGEN-6C4模型比较,大地水准面阶误差见图2,其中,Tongji-GOCE2019S是仅利用GOCE卫星重力梯度数据反演的无约束解模型.从图2可以看出,无约束解模型Tongji-GOCE2019S的精度很差,这主要由于极区空白及带通滤波等问题引起的法方程严重病态,导致位系数的解算误差很大.将Tongji-GOCE2019S法方程与Tongji-Grace02s法方程融合后,解决了GOCE卫星的极区空白等问题,从72阶起有效提高了Tongji-Grace02s模型位系数精度;Tongji-GOGR2019S模型前188阶比参考模型GIF48更接近EIGEN-6C4模型(该模型融合了DIR-R5模型信息)(Förste et al.,2014),说明本文GOCE卫星重力反演方法是有效的.

表1 Tongji-GOGR2019S模型信息

图2 Tongji-GOGR2019S等模型与EIGEN-6C4模型比较大地水准面阶误差

3、Tongji-GOGR2019S重力场模型的精度评定

目前,评价重力场模型精度的主要方法有:(1)与现有的高精度重力场模型在阶方差及全球大地水准面或重力异常等方面进行全面比较;(2)选用独立的重力观测数据进行比较,如卫星测高海洋重力异常数据、区域实测重力异常数据、区域(似)大地水准面模型、区域GNSS水准数据等,下面采用上述两种方法来全面评估Tongji-GOGR2019S模型精度.

3.1模型阶方差分析

为了全面分析Tongji-GOGR2019S模型的精度,我们选择官方发布的第5/6代直接法DIR-R5/R6模型、时域法TIM-R5/R6/R6e模型以及武汉大学在ICGEM发布的GOSG01S模型进行比较分析,各模型与EIGEN-6C4模型比较大地水准面阶误差见图3.图3结果表明:(1)Tongji-GOGR2019S模型比GOSG01S模型更接近EIGEN-6C4模型,主要原因是GOSG01S模型只用GOCE卫星数据解算,没有结合GRACE卫星数据或其他卫星数据(Xu et al.,2017);(2)在133阶前,DIR-R5模型比其他模型更接近EIGEN-6C4模型,主要原因是EIGEN-6C4模型融合了DIR-R5模型信息(Förste et al.,2014);Tongji-GOGR2019S模型和DIR-R6模型差异不大,但要优于TIM-R5/R6/R6e等模型;(3)在133～180阶范围内Tongji-GOGR2019S模型精度优于DIR-R5、TIM-R5/R6等模型,与DIR-R6和TIM-R6e模型相当;(4)超过180阶后,除GOSG01S模型外,其他模型间差异不大,Tongji-GOGR2019S模型精度整体接近DIR-R6模型.

3.2模型全球大地水准面精度分析

为进一步分析Tongji-GOGR2019S模型的精度,我们选择DIR-R5模型、TIM-R5/R6/R6e模型以及GOSG01S模型计算全球1°×1°格网高程异常值,并与同阶次的DIR-R6模型进行比较,截断至220阶次的计算结果见图4,其中截断至100～220阶次差值的具体统计结果见图5.比较分两种方案,方案1:经度-180°～180°,纬度-90°～90°(如图5a);方案2:经度-180°～180°,纬度-83°～83°(如图5b).

图4以及图5结果表明:(1)全球范围内截断至210阶,Tongji-GOGR2019S模型精度优于DIR-R5、TIM-R5/R6/R6e及GOSG01S等模型,与DIR-R6模型最接近;截断至220阶,Tongji-GOGR2019S模型稍逊于DIR-R5和TIM-R6e模型,原因是Tongji-GOGR2019S模型没进行任何约束;(2)扣除极空白区后,截断至160阶,Tongji-GOGR2019S模型与GOCE第6代模型精度一致;161～180阶,Tongji-GOGR2019S模型精度稍优于TIM-R5及DIR-R5模型,与TIM-R6及TIM-R6e模型精度差异很小;超过180阶后,Tongji-GOGR2019S模型与DIR-R5及TIM-R5模型差异不大,稍逊于TIM-R6/R6e模型;(3)Tongji-GOGR2019S模型在南北两极区域精度优于TIM-R5/R6.

3.3利用区域GNSS水准数据评定模型精度

选用8006个美国GNSS水准数据、536个墨西哥GNSS水准数据和212个加拿大GNSS水准数据对Tongji-GOGR2019S等模型的精度进行评价,模型最高阶次取2160,截断阶次以后的位系数用EIGEN-6C4模型补充(如截断阶次为220,则221～2160阶次的位系数采用EIGEN-6C4模型,模型阶次处理方法同样用于3.4节和3.5节).各模型截断至180～220阶次的精度统计结果见表2、表3和表4,表2和表3结果表明,截断至220阶次,Tongji-GOGR2019S模型精度与DIR-R6等GOCE第6代模型精度基本一致,优于GOSG01S模型;表4结果表明,截断至200阶,Tongji-GOGR2019S模型精度与TIM-R5/R6/R6e及DIR-R5/R6模型精度相当,截断至220阶,Tongji-GOGR2019S模型与TIM-R5/R6/R6e及DIR-R5/R6模型之间的差异不大,基本与TIM-R6/R6e相当.从三个区域的GNSS水准数据检核结果来看,Tongji-GOGR2019S模型精度是可靠的.

图3 Tongji-GOGR2019S及其他GOCE模型与EIGEN-6C4模型比较大地水准面阶误差

图4截断至220阶各重力场模型与DIR-R6模型全球大地水准面差值

图5各模型与DIR-R6模型高程异常差值统计结果

表2用美国GNSS水准数据检核各模型精度统计结果(单位:m)

表3用墨西哥GNSS水准数据检核各模型精度统计结果(单位:m)

表4用加拿大GNSS水准数据检核各模型精度统计结果(单位:m)

3.4利用DTU15海洋重力异常数据评定模型精度

选用七个区域的DTU15海洋重力异常数据(分辨率为5′×5′)(Andersen et al.,2019)对Tongji-GOGR2019S等模型的精度进行评价,其区域划分见表5.各模型精度统计结果见图6,整体上,Tongji-GOGR2019S模型精度与同阶次的GOCE卫星第5、6代模型精度相当,说明Tongji-GOGR2019S模型在这些海洋区域精度是可靠的.进一步从图6G可得,Tongji-GOGR2019S模型精度在北极区域要优于DIR-R5/R6和TIM-R5/R6模型精度,与TIM-R6e模型精度相当.

3.5利用EGG2015似大地水准面数据评定模型精度

选用欧洲大陆EGG2015似大地水准面数据(Denker,2016)(范围为40°N—60°N,0°E—30°E,分辨率为10′×15′,约14641个格网值)对Tongji-GOGR2019S等模型的精度进行评价,统计结果见图7.从图7可以看出,在200阶前,Tongji-GOGR2019S与TIM-R5/R6/R6e和DIR-R6模型精度是一致的;截断至220阶,Tongji-GOGR2019S模型精度稍好于GOCE卫星第5代模型精度,比较接近GOCE卫星第6代模型精度.

表5 DTU15海洋重力异常数据全球区域划分一览表

图6利用DTU15海洋重力异常数据检核各模型精度统计结果

图7采用EGG2015似大地水准面数据检核各模型精度统计结果

4、结语

本文基于直接法和IIR带通(5～41 mHz)滤波技术,采用约24个月的GOCE重力梯度数据和13.5 a的GRACE数据(Tongji-Grace02s模型的法方程)联合反演了220阶次的无约束重力场模型Tongji-GOGR2019S,并采用多源数据对Tongji-GOGR2019S模型精度进行全面检核,结论如下:

(1)综合滤波效果和计算效率,IIR带通滤波器阶数选为8阶较为适合;融合GRACE卫星观测数据可有效克服GOCE卫星独立反演重力场的极空白和梯度观测值的低频噪声大等问题,并通过合理确定GOCE卫星各梯度分量的权比,确保了Tongji-GOGR2019S模型的解算精度.

(2)各类检核结果表明,全球范围内,Tongji-GOGR2019S模型精度优于DIR-R5等GOCE卫星第5代模型,比较接近DIR-R6等GOCE卫星第6代模型,主要原因是Tongji-GOGR2019S模型仅使用了约24个月的GOCE梯度数据,且没做任何约束;从南北两极区域来看,Tongji-GOGR2019S模型精度与DIR-R6和TIM-R6e模型相当.综上所述,220阶次的Tongji-GOGR2019S模型在全球范围内精度是可靠的.

陈鑑华,张兴福,陈秋杰,梁建青,沈云中.融合GOCE和GRACE卫星数据的无约束重力场模型Tongji-GOGR2019S[J].地球物理学报,2020,63(09):3251-3262.

基金:国家自然科学基金(41674006,41731069);高分遥感测绘应用示范系统(一期);NSFC-广东联合基金(第二期)(U1501501)资助