梁磊,閆易浩,王長(zhǎng)青,朱紫彤,高銘,鐘敏,于錦海,徐煥
1 中國(guó)科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,武漢 430077 2 地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,西安 710054 3 中山大學(xué)物理與天文學(xué)院,廣東珠海 519082 4 中山大學(xué)測(cè)繪科學(xué)與技術(shù)學(xué)院,廣東珠海 519082 5 中國(guó)科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院,北京 100049
GRACE-FO(GRACE Follow-On)是GRACE(Gravity Recovery and Climate Experiment)的后續(xù)衛(wèi)星,于2018年5月22日成功發(fā)射(Kornfeld et al.,2019).二者是同類型衛(wèi)星,通過(guò)精確測(cè)量?jī)深w衛(wèi)星之間的距離反演地球時(shí)變重力場(chǎng)模型,研究地球質(zhì)量變化(Tapley et al.,2004).GRACE-FO相較于GRACE而言除了搭載K波段微波測(cè)距系統(tǒng)外,還首次實(shí)現(xiàn)了激光干涉測(cè)距,實(shí)現(xiàn)了以nm級(jí)的精度測(cè)量?jī)深w衛(wèi)星之間的距離(Goswami et al.,2021).GRACE與GRACE-FO的數(shù)據(jù)產(chǎn)品分為L(zhǎng)evel-0、Level-1A、Level-1B和Level-2,其中Level-0至Level-1B稱為載荷數(shù)據(jù)處理,Level-1B至Level-2稱為時(shí)變重力場(chǎng)反演(Wu et al.,2006;Case et al.,2010;Wen et al.,2019).
在載荷數(shù)據(jù)處理過(guò)程中,GNSS(Global Navigation Satellite System)接收機(jī)、加速度計(jì)和星間測(cè)距等數(shù)據(jù)處理與姿態(tài)數(shù)據(jù)有關(guān),主要為:第一個(gè)是計(jì)算KBR(K-Band Ranging)天線相位中心改正和GNSS天線相位中心改正,將測(cè)量的信息改正到衛(wèi)星質(zhì)心上;第二個(gè)是將加速度計(jì)(ACC,Accelerometer)測(cè)量的非保守力轉(zhuǎn)換到慣性坐標(biāo)系;第三個(gè)是標(biāo)定KBR天線相位中心位置和衛(wèi)星質(zhì)心位置;第四個(gè)是評(píng)估激光測(cè)距系統(tǒng)的耦合噪聲.姿態(tài)數(shù)據(jù)誤差會(huì)通過(guò)這四個(gè)方面?zhèn)鬟f到重力場(chǎng)模型中,例如Inácio等(2015)分析GRACE星間距變率殘差與KBR天線相位中心改正時(shí),指出姿態(tài)誤差占總誤差的18%左右;Horwath等(2011)通過(guò)KBR天線相位中心改正,指出姿態(tài)誤差導(dǎo)致了星間指向存在偏差,通過(guò)對(duì)姿態(tài)誤差進(jìn)行建模提高了重力場(chǎng)模型解算精度;Wegener等(2020)通過(guò)姿態(tài)數(shù)據(jù)評(píng)估了激光反射棱鏡交點(diǎn)與衛(wèi)星質(zhì)心不重合對(duì)激光測(cè)距的影響.因此確定高精度的重力衛(wèi)星姿態(tài)數(shù)據(jù)是原始載荷數(shù)據(jù)處理的重要內(nèi)容之一.
GRACE與GRACE-FO均搭載了星敏感器和慣性測(cè)量單元測(cè)量衛(wèi)星的姿態(tài),但二者又有較大區(qū)別.對(duì)于GRACE而言,每顆衛(wèi)星上均搭載了兩顆星敏感器和一個(gè)三軸的慣性測(cè)量單元.星敏感器固定安裝于衛(wèi)星左右兩個(gè)面上,三軸的慣性測(cè)量單元固定安裝于衛(wèi)星的內(nèi)部.然而,在發(fā)射之初GRACE-A星的慣性測(cè)量單元一個(gè)軸出現(xiàn)故障,未有冗余觀測(cè)量,GRACE-B星的慣性測(cè)量單元僅在備份模式下開(kāi)啟(Bandikova,2015),所以兩顆衛(wèi)星均未有IMU的測(cè)量數(shù)據(jù),因此GRACE衛(wèi)星姿態(tài)主要由星敏感器提供.受日月光照的影響GRACE只有一顆星敏感器或者無(wú)星敏感器測(cè)量衛(wèi)星姿態(tài),從而導(dǎo)致衛(wèi)星姿態(tài)精度降低或出現(xiàn)大量間斷.借鑒于GRACE上的經(jīng)驗(yàn),GRACE-FO在設(shè)計(jì)之初采用了三顆星敏感器和一個(gè)四軸光纖陀螺儀設(shè)計(jì)方案.相較于GRACE而言,三顆星敏感器設(shè)計(jì)方案減少了姿態(tài)數(shù)據(jù)間斷,并且較多的時(shí)間存在兩顆星敏感器測(cè)量衛(wèi)星的姿態(tài),部分時(shí)間有三顆星敏感器測(cè)量衛(wèi)星姿態(tài);四軸光纖陀螺儀具有更好的測(cè)量性能,固定安裝于衛(wèi)星的內(nèi)部,同時(shí)整個(gè)科學(xué)任務(wù)階段均可觀測(cè)(Kornfeld et al.,2019).
融合多個(gè)載荷數(shù)據(jù),提高最終數(shù)據(jù)產(chǎn)品處理精度,是目前數(shù)據(jù)處理常用的方法.對(duì)于GRACE而言,除了星敏感器測(cè)量衛(wèi)星姿態(tài)數(shù)據(jù)外,加速度計(jì)可以同步測(cè)量衛(wèi)星的角加速度數(shù)據(jù),融合星敏感器數(shù)據(jù)和角加速度數(shù)據(jù),可以提高姿態(tài)數(shù)據(jù)精度,彌補(bǔ)星敏感器數(shù)據(jù)的缺失,減少姿態(tài)數(shù)據(jù)間斷.例如Klinger(2018)在反演重力場(chǎng)模型時(shí),首先基于最小二乘配置法融合了這兩類數(shù)據(jù),構(gòu)建了高精度的姿態(tài)數(shù)據(jù)產(chǎn)品;Goswami等(2018)組合了星敏感器與角加速度計(jì)數(shù)據(jù),提高了姿態(tài)數(shù)據(jù)產(chǎn)品的精度;官方機(jī)構(gòu)JPL(Jet Propulsion Laboratory)基于Kalman濾波融合了這兩類數(shù)據(jù),發(fā)布了RL03版本的姿態(tài)數(shù)據(jù)(Harvey and Sakumura,2019).GOCE衛(wèi)星是另一顆用于地球重力場(chǎng)研究的衛(wèi)星,其上搭載了引力梯度儀,測(cè)量地球引力梯度數(shù)據(jù),可同步獲得測(cè)量頻段內(nèi)的角加速度.由于在處理梯度數(shù)據(jù)時(shí),需要改正的離心力項(xiàng)與衛(wèi)星的角速度是非線性關(guān)系(Stummer et al.,2011;Pail,2005;郭澤華等,2021;趙宇鵬等,2021),這就需要高精度的姿態(tài)信息,由此需要融合星敏感器和梯度儀測(cè)量的角加速度數(shù)據(jù).例如Cesare等(2008)在時(shí)域上通過(guò)Kalman濾波融合了這兩類數(shù)據(jù);Stummer等(2011)和郭澤華等(2021)在頻率域上,基于Wiener濾波融合了這兩類數(shù)據(jù).對(duì)于GRACE-FO而言,星敏感器和IMU是測(cè)量衛(wèi)星姿態(tài)的主要載荷,其中星敏感器對(duì)姿態(tài)的低頻部分敏感,IMU對(duì)姿態(tài)的高頻部分敏感,融合這兩類數(shù)據(jù)可以獲得高精度的姿態(tài)數(shù)據(jù).在時(shí)域上進(jìn)行數(shù)據(jù)融合具有一定的便利性,可以同步評(píng)估IMU的漂移參數(shù),彌補(bǔ)微小的間斷,因此本文選擇在時(shí)域下進(jìn)行星敏感器數(shù)據(jù)與IMU數(shù)據(jù)融合.
關(guān)于姿態(tài)Kalman濾波,Lefferts等(1982)基于狀態(tài)變量(四元數(shù)和陀螺儀的漂移參數(shù))的一階導(dǎo)形式構(gòu)建了姿態(tài)Kalman濾波算法,此時(shí)在計(jì)算狀態(tài)轉(zhuǎn)移矩陣和狀態(tài)轉(zhuǎn)移方程時(shí),其形式了變成了Riccati方程;Trawny和Roumeliotis(2005)與Lefferts方法類似,只是推導(dǎo)過(guò)程將Riccati方程變成級(jí)數(shù)形式.目前關(guān)于這種方法尚未有研究應(yīng)用于實(shí)際的GRACE-FO數(shù)據(jù)處理中.Harvey和Sakumura(2019)將角速度近似處理轉(zhuǎn)換成四元數(shù)形式微小變化量作為Kalman濾波的狀態(tài)轉(zhuǎn)移矩陣,并基于此方法解算了GRACE和GRACE-FO的姿態(tài)數(shù)據(jù),這種近似處理過(guò)程理論上具有一定的瑕疵.本文將以四元數(shù)和陀螺儀的漂移參數(shù)作為狀態(tài)變量,直接構(gòu)建另一種姿態(tài)Kalman濾波融合算法,處理GRACE-FO的姿態(tài)數(shù)據(jù).鑒于此,本文開(kāi)展了GRACE-FO高精度的姿態(tài)數(shù)據(jù)處理的研究.在本文的第1節(jié)是介紹星敏感器處理算法、IMU處理算法和姿態(tài)Kalman濾波算法;第2節(jié)是分析姿態(tài)數(shù)據(jù)處理結(jié)果;第3節(jié)是評(píng)估姿態(tài)數(shù)據(jù)產(chǎn)品對(duì)時(shí)變重力場(chǎng)的影響;第4節(jié)是結(jié)論和總結(jié).
GRACE-FO搭載的三顆星敏感器和IMU記錄的時(shí)標(biāo)是OBC(Onboard Computer)時(shí)標(biāo),首先需要利用TIM1B和CLK1B數(shù)據(jù)將OBC時(shí)標(biāo)改正到GPS(Global Position System)時(shí)標(biāo)上.其次對(duì)2 Hz采樣的多星敏感器數(shù)據(jù)進(jìn)行組合并降采樣至1 Hz,然后通過(guò)姿態(tài)Kalman濾波與8 Hz采樣的IMU數(shù)據(jù)進(jìn)行融合,輸出1 Hz姿態(tài)數(shù)據(jù)產(chǎn)品.下面分別介紹星敏感器數(shù)據(jù)處理、IMU數(shù)據(jù)處理和姿態(tài)Kalman濾波算法,其中星敏感器數(shù)據(jù)處理和IMU數(shù)據(jù)處理中的詳細(xì)細(xì)節(jié)可參考Yang等(2022).
星敏感器測(cè)量的噪聲具有各向異性的特點(diǎn),即繞垂直于星敏感器視軸方向旋轉(zhuǎn)的精度要比繞視軸方向旋轉(zhuǎn)的精度大約高8~10倍(Bandikova and Flury,2014;Stanton,2000).由于在重力場(chǎng)反演過(guò)程中,需要的是SRF (Science Reference Frame)坐標(biāo)系下的衛(wèi)星姿態(tài)數(shù)據(jù).為此,需要將星敏感器坐標(biāo)系下測(cè)量的四元數(shù)轉(zhuǎn)換到SRF坐標(biāo)系中,其中星敏感器坐標(biāo)系與SRF的之間的轉(zhuǎn)換矩陣(即星敏感器安裝矩陣)記錄于QSA1B或者SOE文件中.轉(zhuǎn)換過(guò)程可看成是一個(gè)線性組合的形式,繞星敏感器視軸旋轉(zhuǎn)的精度低會(huì)傳遞到繞科學(xué)坐標(biāo)系Y軸和Z軸旋轉(zhuǎn)中,導(dǎo)致繞這兩個(gè)軸旋轉(zhuǎn)的精度降低.融合多個(gè)星敏感器數(shù)據(jù)可以降低這一影響.目前,多星敏感器融合主要有兩種方法.第一種是根據(jù)星敏感器測(cè)量的噪聲構(gòu)建權(quán)陣,通過(guò)加權(quán)的形式融合多星敏感器數(shù)據(jù)(Romans,2003).第二種是利用星敏感器的視軸構(gòu)建一個(gè)共同參考框架融合多星敏感器數(shù)據(jù)(Mandea et al.,2010).這兩種方法均已用于多種衛(wèi)星,例如CHAMP(Challenging Minisatellite Payload)、GRACE、GRACE-FO和GOCE(Gravity field and steady-state Ocean Circulation Explorer)(Bandikova and Flury,2014;Siemes,2011;Stummer et al.,2011).對(duì)于這兩種方法而言,融合后的姿態(tài)數(shù)據(jù)精度基本一致(Bandikova,2015).本文選擇使用第一種方法,融合多個(gè)星敏感器的測(cè)量值,相關(guān)算法具體參考Romans(2003).
四元數(shù)運(yùn)算有別于實(shí)數(shù)運(yùn)算,需要采用合適的插值算法保證插值后的姿態(tài)具有很好的光滑性.本文采用單位四元素插值樣條曲線對(duì)四元數(shù)進(jìn)行插值,保證插值后的衛(wèi)星姿態(tài)具有二階導(dǎo)連續(xù)特點(diǎn)(邢燕等,2017),具體公式為
GRACE-FO的慣性測(cè)量單元是由4個(gè)四面體形式的光纖陀螺組成,每個(gè)光纖陀螺在其自身的陀螺儀框架下記錄濾波角(filtered angle),采樣率為8 Hz(Wen et al.,2019).由于IMU記錄的每個(gè)軸采樣時(shí)刻并不一致,需要對(duì)每個(gè)軸的時(shí)間進(jìn)行統(tǒng)一.IMU記錄的濾波角的范圍是從-5758°到5758°,當(dāng)記錄到5758°時(shí)自動(dòng)調(diào)整到-5758°開(kāi)始,或與之相反.然而這導(dǎo)致對(duì)濾波角進(jìn)行數(shù)值差分計(jì)算角速度時(shí),存在異常情況,需要對(duì)異常情況進(jìn)行處理.對(duì)于數(shù)值差分計(jì)算,為了保證數(shù)值差分計(jì)算的精度,本文選擇使用Diebel(2006)的方法,具體形式如下:
(1)
由于IMU測(cè)量存在冗余觀測(cè),需要對(duì)冗余觀測(cè)量進(jìn)行處理,根據(jù)Jafari(2015)指出同時(shí)融合冗余觀測(cè)量,可以提高角速度的精度,處理方法如下:
ω=(HTH)-1HTm,
(2)
其中,H4×3是一個(gè)轉(zhuǎn)換矩陣,表示的是4個(gè)光纖陀螺儀測(cè)量軸在IMU坐標(biāo)系下的矢量,記錄于GRACE-FO Level-1A數(shù)據(jù)使用手冊(cè)中(Wen et al.,2019),m是陀螺儀的角速度測(cè)量值,ω是IMU坐標(biāo)系下角速度.
四元數(shù)運(yùn)動(dòng)方程可寫(xiě)成如下形式(秦永元,2014):
(3)
其中q是四元數(shù)矢量,其運(yùn)算法則按照附錄A的形式計(jì)算,ω=[ωx,ωy,ωz],對(duì)四元數(shù)微分方程組(3)按照角增量的形式進(jìn)行求解,省略中間求解過(guò)程,可得:
q(t+1)=Φ(ωΔt)q(t),
(4)
其中
(5)
(6)
陀螺儀的噪聲模型可以寫(xiě)成如下的形式(Lefferts et al.,1982)
(7)
E[εω]=0,
(8)
E[εω(t1),εω(t2)]=Q1(t)δ(t1-t2),
(9)
(10)
(11)
δ是克萊羅符號(hào),Q1與Q2是相應(yīng)的協(xié)方差.將公式(7)代入到公式(4)中,可以寫(xiě)成
(12)
(13)
將公式(12)和(13)代入到公式(4)中,有
X(t+1)=f(t,X(t))+G(t)W(t),
(15)
那么f(t,X(t))可以寫(xiě)成如下形式:
(16)
根據(jù)附錄公式(A3),G(t)W(t)可以寫(xiě)成
(17)
對(duì)公式(17)進(jìn)行泰勒展開(kāi),保留線性形式有
(18)
(19)
四元數(shù)的觀測(cè)方程來(lái)源于星敏感器,有
Z(t)=HX(t)+vk,
(20)
其中H=[I,0],狀態(tài)轉(zhuǎn)移矩陣可以寫(xiě)成
(21)
基于擴(kuò)展Kalman濾波,姿態(tài)Kalman濾波算法過(guò)程按照以下五步進(jìn)行:
第一步:狀態(tài)預(yù)測(cè)方程
(22)
第二步:協(xié)方差矩陣預(yù)測(cè)
P(t+1|t)=Φ(t+1|t)P(t|t)ΦT(t+1|t)+Q(t+1);
(23)
第三步:計(jì)算姿態(tài)Kalman濾波增益
K(t+1)=P(t+1|t)HT(t+1)[H(t+1)P(t+1|t)×HT(t+1)+R(t+1)]-1
(24)
第四步:狀態(tài)更新
(25)
第五步:協(xié)方差更新
P(t+1)=[I-K(t+1)H(t+1)]P(t+1|t).
(26)
為了對(duì)處理后的姿態(tài)數(shù)據(jù)的精度進(jìn)行評(píng)估,本文將與JPL發(fā)布的GRACE-FO Level-1B姿態(tài)產(chǎn)品進(jìn)行對(duì)比分析,即SCA1B數(shù)據(jù).對(duì)比分析的方法主要為星間指向(inter-satellite pointing)和衛(wèi)星角速度矢量.星間指向?yàn)镵BR天線相位中心矢量與LOS(Light of Sight)矢量的偏差并用角來(lái)表示,其中LOS為兩顆衛(wèi)星質(zhì)心之間的連線方向上的矢量通過(guò)軌道數(shù)據(jù)計(jì)算(Bandikova et al.,2012;Bandikova,2015).星間指向主要是用于衛(wèi)星在軌運(yùn)行期間,監(jiān)測(cè)兩顆衛(wèi)星的對(duì)準(zhǔn)情況,是判斷是否需要進(jìn)行在軌控制的依據(jù).由于直接計(jì)算慣性坐標(biāo)系下衛(wèi)星的姿態(tài)角,會(huì)受到較大值的影響難以展示相關(guān)細(xì)節(jié),因此為了便于直觀展示解算的最終結(jié)果,本文采用星間指向進(jìn)行對(duì)比分析.衛(wèi)星角速度矢量用于評(píng)估融合后的姿態(tài)數(shù)據(jù)在中高頻部分是否完全融合了IMU信息.
星間指向計(jì)算過(guò)程如下:
(27)
(28)
求解的φ、θ和ψ即為星間指向上的偏差角.
由于官方機(jī)構(gòu)并未公布Level-0數(shù)據(jù),僅公布了Level-1A和Level-1B數(shù)據(jù),因此本文將從Level-1A數(shù)據(jù)出發(fā)構(gòu)建GRACE-FO的姿態(tài)數(shù)據(jù).本文選擇以2019年1月1日的數(shù)據(jù)為例,對(duì)星敏感器數(shù)據(jù)處理以及與IMU數(shù)據(jù)進(jìn)行姿態(tài)Kalman濾波融合,并對(duì)濾波后的結(jié)果進(jìn)行對(duì)比分析.
首先選擇2019年1月1日中三顆星敏感器同時(shí)存在觀測(cè)數(shù)據(jù)的時(shí)間段進(jìn)行分析.單個(gè)星敏感器解算的衛(wèi)星姿態(tài)結(jié)果如圖1所示,其中SCA1、SCA2和SCA3表示星敏感器的編號(hào)(Kornfeld et al.,2019).從星間指向的Yaw方向可見(jiàn),即圖1中第一行,三個(gè)星敏感器的時(shí)間序列變化并不一致但基本維持在±1×10-3rad之間,從功率譜上看三個(gè)星敏感器在高頻部分的精度基本一致.從星間指向的Pitch方向可見(jiàn),即圖1中第二行,三個(gè)星敏感器的時(shí)間序列變化基本維持在±1×10-3rad之間,并且SCA2和SCA3基本一致且有較多毛刺,而SCA1變化更加平滑.從功率譜上看SCA2和SCA3在高頻部分的精度基本一致,而SCA1具有更低的高頻噪聲.這主要是由于SCA1安裝在衛(wèi)星的天頂方向,星敏感器視軸方向剛好與科學(xué)坐標(biāo)系的Z軸平行,而SCA2和SCA3安裝在衛(wèi)星左右兩個(gè)面方向上,從星敏感器坐標(biāo)系轉(zhuǎn)換到SRF坐標(biāo)系過(guò)程中,受到繞視軸旋轉(zhuǎn)不精確的影響傳遞到Pitch方向上.從星間指向的Roll方向可見(jiàn),即圖中第三行,三個(gè)星敏感器的時(shí)間序列變化趨勢(shì)基本一致,從功率譜上看三個(gè)星敏感器在高頻部分精度處于同一水平.
融合多個(gè)星敏感器計(jì)算的星間指向的結(jié)果如圖2所示,其中SCA12、SCA13、SCA23和SCA123表示融合不同的星敏感器,例如SCA12表示融合SCA1和SCA2兩個(gè)星敏感器,SCA123表示融合SCA1、SCA2和SCA3三個(gè)星敏感器.從星間指向的Yaw、Pitch和Roll方向的時(shí)間序列可見(jiàn),三個(gè)時(shí)間序列的變化趨勢(shì)基本一致.從功率譜上看,與圖1中的功率譜對(duì)比可見(jiàn),融合多個(gè)星敏感器可以提高高頻部分的精度,并且融合三顆星敏感器解算的姿態(tài)的高頻部分的精度要略優(yōu)于融合兩顆星敏感器的結(jié)果.此外,從星間指向的Pitch方向可發(fā)現(xiàn)一個(gè)有趣的現(xiàn)象,融合不同星敏感器計(jì)算的Pitch角存在一個(gè)偏差,最大相差約3.0×10-4rad.我們檢查了相應(yīng)的計(jì)算流程,發(fā)現(xiàn)QSA1B中記錄的星敏感器坐標(biāo)系與SRF坐標(biāo)系的轉(zhuǎn)換四元數(shù)僅提供了地面標(biāo)定值,未提供在軌標(biāo)定值.對(duì)比GRACE的SOE文件可見(jiàn),在GRACE入軌后進(jìn)行了數(shù)次在軌標(biāo)定,提供了相應(yīng)的星敏感器坐標(biāo)系與SRF坐標(biāo)系轉(zhuǎn)換四元數(shù).目前,由于缺少相關(guān)的遙測(cè)數(shù)據(jù),無(wú)法驗(yàn)證這樣的偏差是否來(lái)源于QSA1B,因此在本文的后續(xù)處理過(guò)程中,暫不考慮這樣的微小偏差影響.
融合星敏感器數(shù)據(jù)與IMU數(shù)據(jù)解算衛(wèi)星姿態(tài)數(shù)據(jù),其中首先對(duì)星敏感器數(shù)據(jù)進(jìn)行融合并降采樣至1 Hz,然后再與8 Hz的IMU數(shù)據(jù)進(jìn)行姿態(tài)Kalman濾波融合,在姿態(tài)Kalman濾波過(guò)程中僅保留整數(shù)秒的姿態(tài)數(shù)據(jù).融合后的衛(wèi)星姿態(tài)數(shù)據(jù)結(jié)果如圖2所示,其中APM表示本文所處理的姿態(tài)數(shù)據(jù)結(jié)果,JPL表示官方機(jī)構(gòu)解算的結(jié)果.從星間指向的時(shí)間序列上看,Yaw、Pitch和Roll的時(shí)間序列趨勢(shì)變化基本一致,僅在Pitch方向上存在一個(gè)小的偏差,約為0.6 mrad.從功率譜看上基本一致,僅在截尾部分(即0.2~0.5 Hz處)JPL略低,這里將在后面進(jìn)一步分析.
圖1 單個(gè)星敏感器解算的衛(wèi)星姿態(tài)(a) 左列是星間指向的時(shí)間序列,(b) 右列是對(duì)應(yīng)的功率譜,其中第一至三行分別表示星間指向的Yaw、Pitch和Roll三個(gè)方向.Fig.1 Satellite attitude calculated by a single star camera(a) The left plots are the time series of inter-satellite pointing,and (b) the right plots are the corresponding amplitude spectral density,where the first to third rows represent the Yaw,Pitch,and Roll directions of the inter-satellite pointing,respectively.
圖2 融合多個(gè)星敏感器解算的衛(wèi)星姿態(tài)(a) 左列是星間指向的時(shí)間序列,(b) 右列是對(duì)應(yīng)的功率譜,其中第一至三行分別表示星間指向的Yaw、Pitch和Roll三個(gè)方向.Fig.2 Satellite attitude calculated by fusing multiple star cameras(a) The left plots are the time series of inter-satellite pointing,and (b) the right plots are the corresponding amplitude spectral density,where the first to third rows represent the Yaw,Pitch,and Roll directions of the inter-satellite pointing,respectively.
圖3 星間指向的分析(a) 左列表示相應(yīng)星間指向的時(shí)間序列,(b) 右列表示相應(yīng)的功率譜,其中第一至三行分別表示星間指向的Yaw、Pitch和Roll三個(gè)方向.Fig.3 Analysis of inter-satellite pointing(a) The left plots are the time series of inter-satellite pointing,and (b) the right plots are the corresponding amplitude spectral density,where the first to third rows represent the Yaw,Pitch,and Roll directions of the inter-satellite pointing,respectively.
圖4 衛(wèi)星角速度矢量的分析(a) 左列表示時(shí)間序列,(b) 右列表示相應(yīng)的功率譜,其中第一至三行分別表示ωx、ωy和ωz三個(gè)方向.Fig.4 Analysis of satellite angular velocity vector(a) The left plots are the time series of inter-satellite pointing,and (b) the right plots are the corresponding amplitude spectral density,where the first to third rows represent the ωx,ωy,and ωz.
圖5 慣性坐標(biāo)系下非保守力的比較(a) 左列表示APM與JPL姿態(tài)產(chǎn)品計(jì)算的非保守力差異的時(shí)間序列,(b) 右列是相應(yīng)的功率,其中第一行至第三行分別表示慣性坐標(biāo)系下X,Y和Z方向.Fig.5 Comparison of non-conservative forces in an inertial coordinate system(a) The left plots represent the time series of non-conservative force differences calculated by the APM and JPL attitude data,(b) the right plots represent the corresponding amplitude spectral density,where the first to third rows represent the X,Y,and Z directions of the inertial coordinate system.
圖6 KBR天線相位中心改正變率(a) 時(shí)間序列;(b) 相對(duì)應(yīng)的功率譜.Fig.6 KBR antenna phase correction rate(a) The time-series;(b) The corresponding amplitude spectral density.
圖7 時(shí)變重力場(chǎng)模型的階方差(a) 2018年6月;(b) 2019年3月.Fig.7 The degree variance of the time-variable gravity field model(a) The June 2018;(b) The March 2019.
圖8 質(zhì)量變化等效水高(a) 2018年6月;(b) 2019年3月.Fig.8 The equivalent water height for mass change(a) June 2018;(b) March 2019.
利用公式(3)計(jì)算衛(wèi)星的角速度驗(yàn)證進(jìn)行比較,結(jié)果如圖4所示.圖4中APM、JPL和IMU1B分別表示本文解算的結(jié)果、官方機(jī)構(gòu)結(jié)果和IMU測(cè)量的角速度.APM-IMU1B、JPL-IMU1B分別表示APM和 JPL計(jì)算的角速度與IMU1B測(cè)量值的差異.為了更好的展示細(xì)節(jié),在圖4的左列中僅展示了2019年1月1日前10000秒的時(shí)間序列,其中相對(duì)應(yīng)的功率譜則使用完整的一天數(shù)據(jù)進(jìn)行計(jì)算.從圖4的時(shí)間序列可見(jiàn),{ωx,ωy,ωz}變化趨勢(shì)基本一致,APM與IMU1B的差異在三個(gè)軸的標(biāo)準(zhǔn)差分別為1.89×10-7rad·s-1,1.56×10-7rad·s-1,1.63×10-7rad·s-1,JPL與IMU1B的差異在三個(gè)軸的標(biāo)準(zhǔn)差分別為6.09×10-7rad·s-1,3.62×10-7rad·s-1,4.01×10-7rad·s-1,可見(jiàn)APM相較于JPL而言精度至少提高了3倍,具有更高的精度.從功率譜上看,APM、JPL和IMU1B變化基本一致,但是從與IMU1B產(chǎn)品的差異上看,APM在0.004~0.2 Hz上具有更低的噪聲,特別在0.01~0.1 Hz相較于JPL而言噪聲水平低一個(gè)量級(jí)左右,即圖中青色虛線低于品紅色虛線.換句話說(shuō)本文所解算的姿態(tài)數(shù)據(jù)充分融合IMU1B數(shù)據(jù),具有更低的噪聲水平.在截尾部分(即0.2~0.5 Hz處)可見(jiàn),APM與JPL都有明顯的下降趨勢(shì),這與濾波所選取的截止頻率有關(guān).
正如引言所述,姿態(tài)數(shù)據(jù)主要通過(guò)非保守力轉(zhuǎn)換和KBR天線相位中心改正影響重力場(chǎng)模型解算精度.因此,本節(jié)將對(duì)比分析APM和JPL姿態(tài)數(shù)據(jù)對(duì)非保守力轉(zhuǎn)換以及KBR天線相位中心改正的影響,并進(jìn)一步反演時(shí)變重力場(chǎng)模型分析姿態(tài)數(shù)據(jù)對(duì)時(shí)變重力場(chǎng)模型精度的影響.時(shí)變重力場(chǎng)反演的過(guò)程及等效水高計(jì)算過(guò)程本文參考Liang等(2021)、梁磊等(2019)和Zhou等(2018).
計(jì)算慣性坐標(biāo)系下非保守力的結(jié)果如圖5所示.由于只需要評(píng)估APM與JPL姿態(tài)數(shù)據(jù)計(jì)算慣性坐標(biāo)系下非保守力的差異,因此在圖5左列中的時(shí)間序列只展示了轉(zhuǎn)換后的非保守力的差異,在圖5的右列中展示了JPL與APM轉(zhuǎn)換后的非保守力及其差異的功率譜.對(duì)圖5中的時(shí)間序列進(jìn)行統(tǒng)計(jì),結(jié)果如表1所示.從圖5左列的時(shí)間序列和表1可見(jiàn),通過(guò)APM和JPL姿態(tài)數(shù)據(jù)轉(zhuǎn)換后的非保守力差異的標(biāo)準(zhǔn)差在10-10m·s-2左右,非保守力差異的均值在X和Y軸的方向?yàn)?0-11m·s-2左右,在Z軸的方向均值差異為4.61×10-10m·s-2,也就是說(shuō)在6 h積分弧長(zhǎng)上,非保守力的差異累積不超過(guò)2 cm,低于定軌精度.從圖5右列的功率譜上看高頻部分的差異主要受姿態(tài)數(shù)據(jù)的影響,低頻部分主要是集中于1CPR和2CPR.對(duì)于1CPR和2CPR的低頻信號(hào),在重力場(chǎng)反演過(guò)程中可看成是低頻誤差(Liang et al.,2021;Kim,2000;Zhao et al.,2011;Zhou et al.,2018,2019),通過(guò)低頻誤差處理或者加速度標(biāo)定進(jìn)行吸收,并不影響重力場(chǎng)模型解算精度.
表1 轉(zhuǎn)后非保守力差異的標(biāo)準(zhǔn)差和均值Table 1 Standard deviation and mean of the non-conservative force difference
由于在反演重力場(chǎng)過(guò)程中,主要使用的是星間距變率作為觀測(cè)值,因此這里僅計(jì)算KBR天線相位中心改正變率進(jìn)行分析,計(jì)算方法參考Bandikova(2015),結(jié)果如圖6所示.從圖中可見(jiàn),KBR天線相位中心改正的時(shí)間變率量級(jí)在10-9m·s-1左右,遠(yuǎn)低于KBR測(cè)距精度1 μm·s-1的要求,因此KBR天線相位中心改正的時(shí)間變率對(duì)重力場(chǎng)反演影響很小.此外,可見(jiàn)GRACE-FO對(duì)兩顆衛(wèi)星的對(duì)準(zhǔn)情況要求要更為嚴(yán)格.
選取2018年6月和2019年3月數(shù)據(jù)反演時(shí)變重力場(chǎng)模型進(jìn)行比較,其中反演重力場(chǎng)過(guò)程中只考慮姿態(tài)數(shù)據(jù)的差異,其他數(shù)據(jù)則使用JPL Level-1B產(chǎn)品.反演時(shí)變重力場(chǎng)模型的階方差和等效水高結(jié)果如圖7—8所示.從圖7可見(jiàn),反演時(shí)變重力場(chǎng)模型的階方差基本一致.圖8表示的是兩個(gè)模型計(jì)算的等效水高差異,從圖中可見(jiàn)在全球范圍內(nèi)利用APM與JPL姿態(tài)數(shù)據(jù)計(jì)算的等效水高基本一致,僅在GRACE-FO任務(wù)開(kāi)始階段(2018年6月),在質(zhì)量變化較大的亞馬遜流域上有略微差異,約占最大信號(hào)的5%左右,進(jìn)入到任務(wù)穩(wěn)定時(shí)期(2019年3月),計(jì)算的全球等效水高基本一致.
本文研究了GRACE-FO姿態(tài)數(shù)據(jù)Level-1A至Level-1B的處理算法及其對(duì)時(shí)變重力場(chǎng)的影響.首先根據(jù)GRACE-FO載荷設(shè)計(jì)特點(diǎn),從四元數(shù)運(yùn)動(dòng)微分方程出發(fā),以四元數(shù)和陀螺儀漂移參數(shù)為狀態(tài)變量,構(gòu)建了一種新的姿態(tài)Kalman濾波融合算法.然后以GRACE-FO數(shù)據(jù)為基礎(chǔ),分析了單個(gè)星敏感及多星敏感器組合對(duì)衛(wèi)星姿態(tài)解算的影響,可見(jiàn)融合多星敏感器數(shù)據(jù)可以抑制高頻部分的噪聲.另外,星敏感器數(shù)據(jù)之間的偏差可能來(lái)源于官方機(jī)構(gòu)未提供在軌標(biāo)定的星敏感器安裝矩陣,需要進(jìn)一步利用遙測(cè)數(shù)據(jù)進(jìn)行分析驗(yàn)證.最后基于本文推導(dǎo)的姿態(tài)Kalman濾波算法,處理實(shí)際的GRACE-FO Level-1A數(shù)據(jù),從星間指向和角速度矢量可見(jiàn),本文解算的姿態(tài)數(shù)據(jù)充分發(fā)揮了IMU載荷的性能,具有更低的噪聲水平,優(yōu)于JPL解算的結(jié)果.
對(duì)時(shí)變重力場(chǎng)模型精度的影響,本文首先分析了姿態(tài)數(shù)據(jù)對(duì)轉(zhuǎn)換后的非保守力和KBR天線相位中心改正的影響.從轉(zhuǎn)換后的非保守力看,本文解算的姿態(tài)數(shù)據(jù)與JPL的姿態(tài)數(shù)據(jù)引起的差異在10-10m·s-2左右,從KBR天線相位中心改正上看,本文解算的姿態(tài)數(shù)據(jù)與JPL的姿態(tài)數(shù)據(jù)計(jì)算的KBR天線相位中心改正變率量級(jí)在10-9m·s-1,遠(yuǎn)低于KBR星間變率1 μm·s-1的精度要求.進(jìn)一步解算了2018年6月和2019年3月的時(shí)變重力場(chǎng)進(jìn)行比較,從階方差和等效水高差異上看,本文解算的姿態(tài)數(shù)據(jù)與JPL的姿態(tài)數(shù)據(jù)對(duì)時(shí)變重力場(chǎng)模型的精度影響很小.總的來(lái)說(shuō),本文構(gòu)建了一種的新的姿態(tài)Kalman濾波算法,解算的姿態(tài)數(shù)據(jù)充分發(fā)揮了IMU的性能精度優(yōu)于JPL,但對(duì)時(shí)變重力場(chǎng)模型解算精度的影響很小,目前限制時(shí)變重力場(chǎng)模型的精度主要來(lái)源于其他誤差.
致謝感謝GFZ網(wǎng)站提供的GRACE-FO Level-1A和Level-1B數(shù)據(jù)(ftp:∥isdcftp.gfz-potsdam.de/grace-fo/).感謝審稿專家對(duì)本文提出的寶貴修改意見(jiàn).
附錄A 四元數(shù)的定義和運(yùn)算
稱q=q0+q1i+q2j+q3k為一個(gè)四元數(shù),其中q0,q1,q2,q3∈R,i,j,k是不同的虛數(shù)單位,滿足i2=j2=k2=-1,ij=-ji=k,jk=-kj=i,ki=-ik=j.四元數(shù)也可以表示成(q0,q1,q2,q3),其中q0是標(biāo)量.
令任意三個(gè)四元數(shù)為q=(q0,q1,q2,q3),p=(p0,p1,p2,p3)和s=(s0,s1,s2,s3),那么四元數(shù)的運(yùn)算法則定義如下.
加法:
q±p=(q0±p0,q1±p2,q1±p2,q1±p2),
(A1)
點(diǎn)乘:
q·p=q0p0+q1p1+q2p2+q3p3,
(A2)
乘法:
s=qp,
(A3)
共軛:
(A4)
模:
(A5)
若‖q‖=1,則稱q是單位四元數(shù).
逆:
(A6)
單位四元數(shù)三角函數(shù)形式及指數(shù)形式:
q=cosθ+nsinθ,
(A7)
其中n∈S2為單位矢量.