GRACE和GOCE卫星数据融合下无约束重力场模型Tongji-GOGR2019S研究.docx

GRACE和GOCE卫星数据融合下无约束重力场模型Tongji-GOGR2019S研究引言高精度、高分辨率地球重力场为全球或区域高程基准统一、大地水准面确定、海洋洋流确定、强地震监测等提供更加精细的地球重力场信息(许厚泽等,2012;Beckeretal.,2014;GerlachandRummel,2013;Fuchsetal.,2013).GOCE(GravityfieldandsteadystateOceanCirculationExplorer)卫星于2009年3月发射,2013年11月结束任务,共获取了约42个月的对地观测数据,其主要科学目标是确定高精度的中高阶静态地球重力场,其对应的空间分辨率约为100km,精度为1～2cm的大地水准面和1～2mGal的重力异常.GOCE重力卫星计划的成功实施,对提高全球重力场模型的精度和分辨率取得了显著作用(Pailetal.,2011).国内外科研机构采用GOCE卫星观测数据或联合其他卫星(如GRACE卫星等)观测数据解算了多个版本的全球重力场模型(Brockmannetal.,2014,2019;Bruinsmaetal.,2014;Försteetal.,2019;Kvasetal.,2019;Wuetal.,2016;Xuetal.,2017;Luetal.,2018;苏勇等,2018;Zhouetal.,2019),近来国际上融合新版GOCE梯度数据和GRACE卫星数据解算了精度最高的DIR-R6等系列模型,而在我国,独立研制与国际精度相当的重力场模型对于自主重力卫星计划的实施具有重要的科学意义和现实意义.利用GOCE卫星数据进行重力解算主要有3种方法,即直接法、时域法和空域法(Pailetal.,2011;郑伟等,2014;于锦海和万晓云,2012),欧空局(ESA)官方发布的重力场模型均由这3种方法解算,相应成果已在ICGEM(InternationalCentreforGlobalEarthModels)网站(http://icgem.gfz-potsdam.de/home)发布.直接法是选定某一参考模型,直接在卫星观测轨道处建立扰动重力梯度与重力位系数观测方程,可通过联合GRACE卫星数据及LAGEOS1/2等SLR数据反演重力场(Pailetal.,2011;Bruinsmaetal.,2013,2014;Försteetal.,2019).时域法是将观测值看作时间序列,根据卫星轨道摄动理论构建观测方程,由于时域法未联合GRACE等其他卫星数据,其解算的模型也称为纯GOCE卫星重力场模型(Pailetal.,2011;Brockmannetal.,2014,2019),通过进一步融合两极地面重力数据反演了TIM-R6e模型(Zingerleetal.,2019).空域法是基于大地测量边值理论构建梯度观测值与位系数关系式,并通过多步最小二乘配置法解算位系数,解算的模型也为纯GOCE卫星重力场模型(ReguzzoniandTselfes,2009;Pailetal.,2011;GattiandReguzzoni,2017).为了克服GOCE卫星极区空白及法方程病态等问题,通常采用球冠正则化、Kaula正则化等正则化方法(Försteetal.,2019).综合看,直接法解算重力位系数的难度和计算量最大,而时域法和空域法通过一些假设或近似,减小了计算量,但直接法解算的DIR_Rx系列模型精度要明显优于时域法的TIM_Rx系列和空域法的SPW_Rx系列模型.目前,由GOCE完整观测周期的数据求得的第5、6代(R5、R6)重力场模型已基本实现GOCE卫星的设计目标(Brockmannetal.,2014;Bruinsmaetal.,2014;Gruber,2019),而第6代模型由于采用了新版GOCE卫星重力梯度数据(Gruber,2019;Siemesetal.,2019),模型精度优于第5代.此外,GOCO团队(GravityObservationCombination)联合GOCE、GRACE等卫星数据也获得了与同期官方模型精度相当的GOCO01S、GOCO03s、GOCO05s以及GOCO06s等模型(Pailetal.,2010;Kvasetal.,2019),其中,GOCO06s模型还采用了Swarm等低轨卫星数据以及Starlette等SLR数据,是目前采用卫星数据最多的重力场模型;Schall等(2014)和Yi等(Yi,2012a,b;Yietal.,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卫星任务目标(确定分辨率约为100km,大地水准面精度为1～2cm的重力场),本文选择直接法作为重力反演方法,反演220阶次的重力场.基于直接法和IIR带通滤波技术,利用约24个月的GOCE卫星重力梯度数据构建了反演220阶次重力场模型的法方程,并与由13.5年GRACE卫星数据构建的180阶次的Tongji-Grace02s模型的法方程融合,解算了220阶次的Tongji-GOGR2019S重力场模型,最后采用多源数据对Tongji-GOGR2019S模型精度进行全面检核.1、GOCE和GRACE卫星观测模型与联合反演1.1基于GOCE卫星重力梯度观测数据的重力反演在球坐标系中,地球外部空间任意一点的引力位可表示为(Pailetal.,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¯¯¯nm,S¯¯nm)为完全规格化的重力位系数,P¯¯¯nm(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～100mHz测量带宽内,为了保留频带内信号,并尽量减小频带外噪声对重力反演的影响,本文采用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.02m,与Tongji-Grace02s模型单位权一致,梯度分量精度采用相对于参考模型的统计结果,与GRACE卫星有关的法方程等信息由Tongji-Grace02s重力场模型对应的法方程提供(Chenetal.,2018).2、Tongji-GOGR2019S重力场模型解算2.1GOCE梯度观测值滤波策略本文采用GOCE卫星的Vxx、Vyy、Vzz和Vxz4个重力梯度分量观测值反演重力场,考虑到GOCE重力梯度观测值的噪声特性和实际反演的模型阶数,本文采用带宽为5～41mHz的IIR带通滤波器对扰动重力梯度观测值进行滤波,截断频率41mHz大致对应能反演重力场模型的最高阶次为220(万晓云等,2012;Xuetal.,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个月经过时变重力场变化改正后的Level1b重力梯度观测值,并联合13.5年的GRACE卫星数据(180阶次的Tongji-Grace02s模型的法方程)反演了220阶次的Tongji-GOGR2019S重力场模型,其中,Tongji-Grace02s模型是基于改进的短弧积分法,采用多项式对GRAC。