衛(wèi)星重力梯度觀測(cè)數(shù)據(jù)L1級(jí)構(gòu)建方法

吳云龍， 郭澤華， 肖云， 馬林

1 地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室, 中國(guó)地震局地震研究所, 武漢 430071 2 防災(zāi)科技學(xué)院, 河北三河 065201 3 引力與固體潮國(guó)家野外科學(xué)觀測(cè)研究站， 武漢 430071 4 西安測(cè)繪研究所， 西安 710054 5 航天東方紅衛(wèi)星有限公司， 北京 100094

0 引言

全球重力場(chǎng)和海洋環(huán)流探測(cè)器(Gravity field and steady-state Ocean Circulation Explorer, GOCE)由歐洲空間局(European Space Agency，ESA)于2009年3月17日成功發(fā)射，在軌運(yùn)行4年后于2013年11月11日結(jié)束衛(wèi)星任務(wù).GOCE衛(wèi)星主要任務(wù)目標(biāo)是在100 km的空間分辨率下確定精度為1～2 cm的全球大地水準(zhǔn)面和精度為1 mGal的全球重力異常(ESA, 1999).為實(shí)現(xiàn)這一極具難度的科學(xué)目標(biāo)，GOCE衛(wèi)星搭載了高精度靜電重力梯度儀(Electrostatic Gravity Gradiometry，EGG)、高精度全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System, GNSS)星載接收機(jī)和星敏感器(Star Sensor，STR)等多個(gè)關(guān)鍵載荷(Cesare, 2008; Rummel et al., 2011).

高精度GOCE重力梯度觀測(cè)數(shù)據(jù)在地球科學(xué)的多個(gè)領(lǐng)域應(yīng)用廣泛.國(guó)際地球重力場(chǎng)模型中心(International Centre for Global Earth Models, ICGEM)發(fā)布的全球最優(yōu)重力場(chǎng)模型產(chǎn)品系列中，相當(dāng)一部分產(chǎn)品來(lái)源于衛(wèi)星重力梯度觀測(cè)數(shù)據(jù)(Drewes et al.，2016).國(guó)內(nèi)外研究同行將GOCE重力梯度觀測(cè)數(shù)據(jù)及其模型產(chǎn)品在大地測(cè)量學(xué)、固體地球物理學(xué)、海洋學(xué)以及冰川學(xué)等領(lǐng)域開(kāi)展了持續(xù)深入的應(yīng)用研究(Sun and Okubo, 2004; Tapley et al., 2004; 孫文科, 2008; 鐘敏等, 2009; Knudsen et al., 2011; Gruber et al., 2012; Becker et al., 2014; Fuchs et al., 2014; Hirt, 2014; Bouman et al., 2014).但是，所有這些地學(xué)領(lǐng)域應(yīng)用都取決于能否獲取高精度重力梯度測(cè)量觀測(cè)數(shù)據(jù).此外，為更好發(fā)揮重力梯度儀的觀測(cè)性能，GOCE衛(wèi)星設(shè)計(jì)運(yùn)行在極低的軌道高度(～259 km)，以最大效應(yīng)感知重力梯度信息，從而對(duì)衛(wèi)星軌道高度維持和高精度姿態(tài)控制都提出了更高要求(Floberghagen et al., 2011; Pail et al., 2011).GOCE衛(wèi)星官方數(shù)據(jù)處理部門將構(gòu)建高精度重力梯度數(shù)據(jù)作為L(zhǎng)1級(jí)數(shù)據(jù)處理流程的核心環(huán)節(jié)，主要為將靜電重力梯度儀得到的加速度數(shù)據(jù)、星敏感器得到的姿態(tài)數(shù)據(jù)，聯(lián)合構(gòu)建高精度重力梯度分量觀測(cè)數(shù)據(jù)(ESA, 2009).因此，高精度的重力梯度數(shù)據(jù)L1級(jí)預(yù)處理可為恢復(fù)高精度靜態(tài)地球重力場(chǎng)模型提供數(shù)據(jù)質(zhì)量保障，是衛(wèi)星重力梯度測(cè)量數(shù)據(jù)處理及應(yīng)用中的重要環(huán)節(jié)，也是實(shí)現(xiàn)GOCE衛(wèi)星預(yù)期科學(xué)目標(biāo)的關(guān)鍵任務(wù)之一(吳云龍，2010).

國(guó)際上，歐洲空間局GOCE衛(wèi)星任務(wù)的高級(jí)數(shù)據(jù)處理部門(High Level Processing Facility，HPF)負(fù)責(zé)L1級(jí)和L2級(jí)官方數(shù)據(jù)產(chǎn)品的處理和發(fā)布(Frommknecht et al., 2011).其分布在歐洲多個(gè)研究機(jī)構(gòu)的10個(gè)研究小組圍繞整個(gè)數(shù)據(jù)處理流程進(jìn)行了系統(tǒng)深入的研究，主要包括角速度重建、梯度觀測(cè)數(shù)據(jù)構(gòu)建、科學(xué)數(shù)據(jù)預(yù)處理、衛(wèi)星軌道精密定軌、重力場(chǎng)模型恢復(fù)(快速產(chǎn)品和科學(xué)產(chǎn)品)以及數(shù)據(jù)極空白改善等(Siemes et al., 2012; Wan and Yu, 2013; Baur et al., 2014; Visser et al., 2014, 2016; Visser, 2017; Siemes, 2018a, b; Siemes et al., 2019a, b).近年來(lái)，我國(guó)圍繞發(fā)展民用重力梯度測(cè)量衛(wèi)星也開(kāi)展了相關(guān)預(yù)研究工作，仿真設(shè)計(jì)了新型重力測(cè)量衛(wèi)星模式，模擬分析了不同重力測(cè)量衛(wèi)星系統(tǒng)配置條件、關(guān)鍵載荷技術(shù)指標(biāo)以及噪聲水平下的重力場(chǎng)模型反演能力(冉將軍等，2015; 徐新禹等，2018)；基于地面實(shí)測(cè)重力數(shù)據(jù)給出了重力梯度測(cè)量衛(wèi)星極空白的彌補(bǔ)改善策略(Lu et al.，2020)；系統(tǒng)研究了重力梯度測(cè)量衛(wèi)星數(shù)據(jù)預(yù)處理方法，包括重力梯度觀測(cè)數(shù)據(jù)預(yù)處理中的潮汐/非潮汐效應(yīng)改正、粗差探測(cè)方法、外部校準(zhǔn)方法等(吳云龍, 2010)；分析了星敏感器噪聲特性，提出了聯(lián)合兩個(gè)或多個(gè)星敏感器姿態(tài)數(shù)據(jù)求解最佳姿態(tài)四元數(shù)方法(郭澤華等, 2021)；提出了一種基于地面重力的衛(wèi)星在軌檢校方法，從地面先驗(yàn)重力數(shù)據(jù)的時(shí)空精度、重力梯度儀觀測(cè)噪聲等預(yù)處理中的外部檢校環(huán)節(jié)開(kāi)展了分析研究(瞿慶亮等，2021).總體而言，高精度重力梯度觀測(cè)數(shù)據(jù)L1級(jí)構(gòu)建的系統(tǒng)方法與技術(shù)，受限于國(guó)外原始數(shù)據(jù)和具體技術(shù)文獻(xiàn)公開(kāi)，國(guó)內(nèi)少有學(xué)者在此研究領(lǐng)域開(kāi)展較為系統(tǒng)的研究工作.隨著我國(guó)推進(jìn)自主重力衛(wèi)星任務(wù)的發(fā)展需要，突破國(guó)外機(jī)構(gòu)對(duì)重力衛(wèi)星L0-L1級(jí)數(shù)據(jù)預(yù)處理的技術(shù)封鎖，實(shí)現(xiàn)自主衛(wèi)星重力梯度觀測(cè)數(shù)據(jù)L1級(jí)構(gòu)建能力，具有迫切技術(shù)需求(許厚澤, 2001; 寧津生, 2002).

本文面向我國(guó)發(fā)展的梯度衛(wèi)星的任務(wù)需要，系統(tǒng)研究并初步實(shí)現(xiàn)了自主的衛(wèi)星重力梯度觀測(cè)數(shù)據(jù)L1級(jí)構(gòu)建方法，主要包括加速度計(jì)電壓數(shù)據(jù)轉(zhuǎn)換、多星敏感器聯(lián)合姿態(tài)數(shù)據(jù)的角速度重建、衛(wèi)星重力梯度分量構(gòu)建等數(shù)據(jù)處理技術(shù)環(huán)節(jié).研究工作可為下一步我國(guó)推進(jìn)實(shí)施民用重力梯度測(cè)量衛(wèi)星任務(wù)提供自主的原始數(shù)據(jù)處理技術(shù)支撐與儲(chǔ)備.

1 衛(wèi)星重力梯度觀測(cè)數(shù)據(jù)構(gòu)建原理

從重力梯度測(cè)量衛(wèi)星原始觀測(cè)數(shù)據(jù)構(gòu)建重力梯度分量，主要以重力梯度儀和星敏感器兩種關(guān)鍵載荷得到的觀測(cè)數(shù)據(jù)為基礎(chǔ)(Rispens and Bouman, 2009).在衛(wèi)星重力梯度測(cè)量中，觀測(cè)量同時(shí)包括了離心加速度和角加速度.通過(guò)觀測(cè)量對(duì)稱性和非對(duì)稱性的表達(dá)式，將角加速度分離.通過(guò)衛(wèi)星搭載的星敏感器提供重力梯度儀的初始加速度矢量，對(duì)其積分可求得重力梯度儀的角速度分量，從而得到離心加速度，最終可將重力梯度從觀測(cè)值中分離出來(lái)(Rummel et al., 2011).

為構(gòu)建重力梯度分量，需要將其他信號(hào)由加速度計(jì)觀測(cè)值中分離.根據(jù)衛(wèi)星重力梯度測(cè)量原理(Stummer et al., 2012)，重力梯度儀中加速度計(jì)觀測(cè)的加速度值ai為：

(1)

圖1 梯度儀坐標(biāo)系內(nèi)梯度儀的6個(gè)加速度計(jì)位置Fig.1 Position of six accelerometers of gradiometer in gradiometer reference frame

GOCE衛(wèi)星所搭載的6個(gè)加速度計(jì)都是通過(guò)剛性連接固定在梯度儀內(nèi)，因此它們?cè)谛l(wèi)星質(zhì)心處的非保守力加速度相同，加速度計(jì)對(duì)距離矢量ri+rj≈0，ij∈{14，25，36}，可通過(guò)建立共模加速度確定.

定義共模加速度acij和差分加速度adij：

(2)

非保守力加速度d在作差分時(shí)相互抵消：

(3)

Ad=[ad,14ad,25ad,36],

(4)

(5)

利用V與Ω2的對(duì)稱性可得：

AdL-1+(AdL-1)T=-2 (V-Ω2),

(6)

(7)

由式(7)可提取角加速度:

(8)

(9)

(10)

2 L1級(jí)觀測(cè)數(shù)據(jù)構(gòu)建方法

重力梯度測(cè)量衛(wèi)星觀測(cè)數(shù)據(jù)L1級(jí)構(gòu)建主要包括電壓數(shù)據(jù)轉(zhuǎn)換、姿態(tài)數(shù)據(jù)重建、角速度重建和重力梯度分量構(gòu)建等四個(gè)步驟，具體計(jì)算方法如下：

圖2 單個(gè)加速度計(jì)電極結(jié)構(gòu)Fig.2 Electrode structure of single accelerometer

2.1 加速度計(jì)電壓數(shù)據(jù)轉(zhuǎn)換

重力梯度儀中單個(gè)加速度計(jì)參考系(Accelerometer Reference Frame, ARF)的電極結(jié)構(gòu)如圖2所示，在加速度計(jì)的較小壁上各有兩個(gè)電極(Y1,Y2,Z1,Z2)，在加速度計(jì)的較大壁上各有四個(gè)電極(X1,X2,X3,X4)，其中兩個(gè)相對(duì)的電極形成一對(duì).檢測(cè)質(zhì)量為扁平立方體(尺寸為4 cm×4 cm×1 cm)，在1 g環(huán)境下進(jìn)行地面試驗(yàn)標(biāo)定，利用較大壁上的四對(duì)電極完成對(duì)檢測(cè)質(zhì)量的懸浮.

Frommknecht等(2011)對(duì)控制電壓-加速度轉(zhuǎn)換的處理步驟做出了詳細(xì)闡述.其關(guān)鍵環(huán)節(jié)是對(duì)控制電壓觀測(cè)數(shù)據(jù)施加靜電增益(即增益因子)，再轉(zhuǎn)化為加速度觀測(cè)數(shù)據(jù).需要注意的是，每個(gè)電極對(duì)增益因子均略有不同.將加速度計(jì)中8個(gè)極板的電壓觀測(cè)值通過(guò)與適當(dāng)?shù)撵o電增益重組矩陣相乘，將控制電壓Vci轉(zhuǎn)換為線性加速度與角加速度.

(11)

(12)

其中，VP為極化電壓；ε0為真空介電常數(shù)；m為檢測(cè)質(zhì)塊的質(zhì)量；e為檢測(cè)質(zhì)塊與極板間的距離.

2.2 多星敏感器聯(lián)合重建姿態(tài)數(shù)據(jù)

根據(jù)郭澤華等(2021)提出的多星敏感器聯(lián)合算法，構(gòu)建噪聲分布加權(quán)矩陣，可獲得衛(wèi)星最佳慣性姿態(tài)四元數(shù).聯(lián)合解算得到的角速度不會(huì)受到單個(gè)星敏感器的視軸測(cè)量精度較低的影響，能夠有效抑制由于坐標(biāo)系變換(SSRF-GRF)導(dǎo)致的精度較低角速度誤差傳播到其他分量(Siemes, 2018b).將星敏感器測(cè)得的四元數(shù)建模，計(jì)算得到最佳四元數(shù)q*為：

(13)

2.3 角速度重建

衛(wèi)星角速度與姿態(tài)四元數(shù)的精度直接影響地球重力場(chǎng)模型反演精度.星敏感器、重力梯度儀觀測(cè)數(shù)據(jù)確定角速度時(shí)，會(huì)出現(xiàn)兩者角速度噪聲分別在高頻、低頻增大的影響.為了精確地確定衛(wèi)星在空間中的慣性角速度，其重建過(guò)程中應(yīng)充分考慮到EGG與STR的誤差特性.基于Stummer等(2011)提出的維納濾波重建角速度方法，可將EGG與STR的角速度根據(jù)其精度在頻域內(nèi)聯(lián)合.表1、表2分別為EGG和STR的噪聲模型，STR和EGG三個(gè)角速度分量具有不同的噪聲功率頻譜密度，如圖3所示.

表1 EGG角速度噪聲模型Table 1 EGG angular velocity noise model

表2 STR角速度噪聲模型Table 2 STR angular velocity noise model

圖3 STR和EGG角速度噪聲的PSD1/2Fig.3 PSD1/2 of STR and EGG angular velocity noise

由于重力梯度儀和星敏感器角速度的精度與頻率有關(guān)，而功率譜密度P(f)可表示在頻率f處的精度，可根據(jù)公式(14)、(15)計(jì)算維納濾波的權(quán)重(Stummer et al., 2012)：

(14)

(15)

所有頻率的權(quán)重之和等于1，由此可得到：

(16)

如圖4所示，實(shí)線、虛線分別表示STR、EGG的維納濾波權(quán)重.由于維納濾波是在時(shí)域中聯(lián)合角速度，濾波系數(shù)可通過(guò)離散傅里葉逆變換(Inverse Discrete Fourier Transform, IDFT)從權(quán)重中獲得：

(17)

(18)

圖4 EGG和STR維納濾波權(quán)重Fig.4 EGG and STR Wiener filter weight

2.4 重力梯度數(shù)據(jù)構(gòu)建

基于上述公式得到的加速度計(jì)觀測(cè)數(shù)據(jù)和角速度觀測(cè)數(shù)據(jù)，綜合式(6)和(7)：

(19)

(20)

衛(wèi)星星體在3個(gè)方向上所受到的非保守力可通過(guò)加速度計(jì)測(cè)量的共模加速度表示，其共模加速度存在以下關(guān)系：

E{ac,ij-ac,kl}=0;ij=14,25,36;ij≠kl

(21)

即共模加速度間存在以下線性組合關(guān)系：

ac,14,x-ac,25,x=0,

ac,14,x-ac,36,x=0,

ac,25,y-ac,36,y=0,

ac,25,y-ac,14,y=0,

ac,36,z-ac,25,z=0,

ac,36,z-ac,14,z=0.

(22)

通過(guò)式(3)得到差分加速度后，只有角速度保留在差分加速度中.為了從式(7)中得到重力梯度，應(yīng)從差模加速度中分離出重力梯度.即V的主對(duì)角線重力梯度分量計(jì)算式為：

(23)

(24)

(25)

非對(duì)角線重力梯度分量計(jì)算式為：

(26)

(27)

(28)

式中，ωx、ωy、ωz為衛(wèi)星的角速度.高精度的衛(wèi)星角速度ωx、ωy和ωz為解算梯度V過(guò)程中的關(guān)鍵參數(shù).

3 計(jì)算結(jié)果與分析

3.1 數(shù)據(jù)

考慮到GOCE衛(wèi)星已有官方發(fā)布的觀測(cè)數(shù)據(jù)，可作為本文構(gòu)建的算法結(jié)果進(jìn)行比對(duì).本文選擇GOCE衛(wèi)星關(guān)鍵載荷的原始數(shù)據(jù)進(jìn)行計(jì)算，數(shù)據(jù)長(zhǎng)度為2013年10月8日至10月14日(共7天)，ESA(2006)提供了GOCE觀測(cè)數(shù)據(jù)內(nèi)容和格式的詳細(xì)描述，文中所用的數(shù)據(jù)文件類型如下：EGG_NOM：加速度計(jì)控制電壓數(shù)據(jù)；STR_VC2, STR_VC3：星敏感器姿態(tài)數(shù)據(jù)；AUX_EGG_DB: 重力梯度儀臂長(zhǎng)及SSRF至GRF轉(zhuǎn)換矩陣.

3.2 加速度計(jì)電壓數(shù)據(jù)轉(zhuǎn)換

本文針對(duì)GOCE衛(wèi)星6個(gè)加速度計(jì)(A1—A6)展開(kāi)分析計(jì)算，通過(guò)對(duì)電壓數(shù)據(jù)轉(zhuǎn)換得到加速度數(shù)據(jù)，并與ESA發(fā)布的L1b數(shù)據(jù)進(jìn)行對(duì)比分析.表3為解算的加速度計(jì)A1—A6增益因子.

從表3可以看出，如果不考慮漂移，Gy、Gz的增益因子變化很小，Gx的增益因子變化則比較明顯.Gx的增益因子約為(3.64～4.51)×10-4m·s-2·V-1，Gy、Gz的增益因子約為(9.67～9.89)×10-7m·s-2·V-1.以加速度計(jì)A1計(jì)算結(jié)果為例，圖5為加速度計(jì)A1各個(gè)軸線性加速度與ESA發(fā)布的L1b數(shù)據(jù)中NLA_A1(Acceleration Nominal Linear, NLA)差值，圖6為加速度計(jì)A1的計(jì)算值Cal_A1與ESA中NLA_A1各軸線性加速度功率譜密度(Power Spectral Density，PSD).表4給出了A1計(jì)算值與NLA_A1各個(gè)軸加速度的統(tǒng)計(jì)特性.

從圖5可以看出，計(jì)算所得的A1線性加速度值與官方發(fā)布加速度數(shù)據(jù)，在x、y、z軸上分量ax,ay,az的差值均在非常小的范圍內(nèi).圖5中黑線為各軸差值的均值線，x、y、z軸對(duì)應(yīng)均值分別為-8.6×10-14m·s-2、-3.2×10-11m·s-2、1.9×10-13m·s-2，各軸最大差值分別為-12.7×10-14m·s-2、-4.2×10-11m·s-2、2.5×10-13m·s-2.由于NLA_A1加速度分量ax,ay,az的量級(jí)分別為10-8m·s-2，10-6m·s-2，10-7m·s-2，其差值量級(jí)體現(xiàn)出很強(qiáng)的相似性.如圖6所示為A1各軸加速度功率譜密度，表明計(jì)算結(jié)果與NLA_A1中對(duì)應(yīng)數(shù)據(jù)在功率譜密度上均呈現(xiàn)一致性(Welch, 1967)，其中ax、az分別在低頻段2～7 mHz、高頻段300～500 mHz內(nèi)精度稍低，其余頻段各軸精度幾乎相同.表4中各軸角速度的統(tǒng)計(jì)特性也體現(xiàn)了計(jì)算結(jié)果與GOCE官方發(fā)布數(shù)據(jù)在各個(gè)軸上高度相似性，其標(biāo)準(zhǔn)差均非常接近，差異極小，達(dá)到了1.1×10-14～2.1×10-12m·s-2量級(jí)，精度結(jié)果十分相近.

表3 加速度計(jì)A1—A6增益因子(單位：m·s-2·V-1)Table 3 Gain factors of accelerometers A1—A6 (unit：m·s-2·V-1)

表4 加速度計(jì)A1與NLA_A1各個(gè)軸加速度的統(tǒng)計(jì)Table 4 Statistics of triaxial acceleration of accelerometer A1 and NLA_A1

圖5 加速度計(jì)A1三軸線性加速度差值Fig.5 Triaxial line acceleration difference of accelerometer A1

圖6 A1三軸線性加速度計(jì)算結(jié)果與NLA_A1的PSD1/2對(duì)比Fig.6 PSD1/2 comparison between triaxial acceleration calculation results and NLA_A1

各個(gè)加速度計(jì)觀測(cè)的加速度值通過(guò)式(2)、(3)進(jìn)一步計(jì)算出共模與差分加速度，結(jié)合式(19)—(21)逆校準(zhǔn)矩陣ICM，得到校正共模加速度(Calibrated Common Mode, CCM)和校正差分加速度(Calibrated Differential Mode, CDM).圖7為重力梯度儀x、y、z軸校正后共模與差分加速度平方根功率譜密度.其中差分加速度在30～200 mHz內(nèi)呈現(xiàn)白噪聲特性，CDM-Y25與CDM-Z36在該頻段范圍內(nèi)精度最高，達(dá)到了2×10-11、5×10-11m·s-2·Hz-1/2量級(jí)，其余各差分加速度約在10-9m·s-2·Hz-1/2量級(jí).共模加速度呈現(xiàn)出隨著頻率升高其精度越高的特點(diǎn)，各軸共模加速度在測(cè)量帶寬(Measurement Band Width, MBW)范圍內(nèi)(5～100 mHz)趨勢(shì)一致且精度相似.

重力梯度儀由6個(gè)三軸加速度計(jì)組成，每個(gè)加速度計(jì)都有2個(gè)超靈敏軸(Ultra Sensitive, US)和1個(gè)低靈敏度軸(Less Sensitive, LS).由于所有6個(gè)加速度計(jì)超靈敏軸均固定安裝在GRF中的x方向上(如圖1)，對(duì)x軸方向共模加速度作差, 可檢驗(yàn)計(jì)算共模加速度是否達(dá)到重力梯度儀的設(shè)計(jì)精度.由誤差傳播公式：

(29)

(30)

式中，σa_comb為共模加速度差值的標(biāo)準(zhǔn)差；σa_c/d_US為共模/差分加速度的標(biāo)準(zhǔn)差；σa_US為單個(gè)US軸標(biāo)準(zhǔn)差.將式(30)代入式(29)可得：

σa_comb≈σa_US.

(31)

圖7 x、y、z軸上加速度計(jì)對(duì)差分(a)與共模(b)加速度PSD1/2Fig.7 x, y and z axis accelerometers couples differential mode (a) and common mode (b) acceleration PSD1/2

圖8為x軸共模加速度作差后的功率譜密度，三條曲線為US軸加速度的真實(shí)噪聲，在測(cè)量帶寬5～100 mHz內(nèi)各US軸噪聲達(dá)到了10-10～10-11m·s-2·Hz-1/2范圍內(nèi)的良好指標(biāo)，符合GOCE衛(wèi)星重力梯度儀的設(shè)計(jì)要求.

圖8 x軸共模加速度的差值PSD1/2Fig.8 The difference of x axis common mode acceleration PSD1/2

3.3 角速度重建

圖9為GRF下多星敏感器聯(lián)合前后各軸角速度平方根的功率譜密度.從圖9a與圖9b可以看出，與單個(gè)STR相比，多個(gè)STR聯(lián)合后顯示出y、z軸的角速度分量wy、wz的總體精度提升明顯；在10～100 mHz處逐步下降了約1個(gè)量級(jí)，精度達(dá)到10-6rad·s-1·Hz-1/2量級(jí)；在3～30 mHz內(nèi)wy的精度最高，而在30 mHz以上wx精度最高，wy與wz精度幾乎相同，且wx與wy、wz均呈現(xiàn)出精度變化趨勢(shì)相似.通過(guò)對(duì)圖9a與圖9b的分析可得，聯(lián)合多星敏感器能有效抑制由于SSRF-GRF坐標(biāo)系變換導(dǎo)致精度較低的角速度噪聲傳播到其他分量.

圖9 GRF下多星敏感器聯(lián)合前后各軸角速度的PSD1/2Fig.9 PSD1/2 of angular velocity of each axis before and after multiple star sensor combination under GRF

圖10為基于維納濾波的重建角速度與EGG各軸角速度平方根功率譜密度，結(jié)果表明在100 mHz頻率范圍內(nèi)，重建角速度相較于EGG角速度的精度有明顯提升.在測(cè)量帶寬5～100 mHz內(nèi)，角速度分量ωx、ωy、ωz平方根功率譜密度均約在32 mHz處精度改進(jìn)最大，其平方根功率譜密度最大改進(jìn)值范圍是(5.21～6.56)×10-11rad·s-1·Hz-1/2，其中ωy的精度提升改進(jìn)最大約為6.56×10-11rad·s-1·Hz-1/2.

圖10 基于維納濾波重建角速度與EGG各軸角速度的PSD1/2Fig.10 Angular velocity reconstructed based on Wiener filter and EGG angular velocity PSD1/2

3.4 重力梯度分量構(gòu)建結(jié)果

式(6)計(jì)算結(jié)果中仍含有離心角速度Ω，應(yīng)被減去才能得到重力梯度分量Vij；再將校準(zhǔn)差分加速度與恢復(fù)得到的角速度代入式(23)—(28)計(jì)算得到重力梯度分量；最后基于拉普拉斯跡準(zhǔn)則進(jìn)行檢驗(yàn).圖11為計(jì)算所得重力梯度分量與GOCE任務(wù)發(fā)布的重力梯度分量EGG_GGT (Gravity Gradients Tensor, GGT)的平方根功率譜密度分析，表5給出了重力梯度各分量計(jì)算值與GOCE任務(wù)重力梯度分量的統(tǒng)計(jì)特性.

從圖11中可以看出，6個(gè)梯度分量的平方根功率譜密度均約在1.85×10-4Hz處達(dá)到峰值，即在該頻率處的精度達(dá)到最低，Vxy約為1.7×105mE·Hz-1/2，其余分量達(dá)到了(1.5～5.0)×106mE·Hz-1/2，其原因是受到STR姿態(tài)信號(hào)噪聲的影響.STR姿態(tài)噪聲以每轉(zhuǎn)周期(cycle per revolution, cpr)的頻率(1 cpr≈1.85×10-4Hz)表現(xiàn)出明顯的重復(fù)性，即在kcpr頻率處(k為整數(shù))6個(gè)梯度分量精度均受到的影響最大，并且在相應(yīng)的頻率下向較高的頻率呈1/f幅度減小.圖11中6個(gè)重力梯度分量中對(duì)應(yīng)的平方根功率譜密度均顯示出整體較為吻合，趨勢(shì)較為一致.在5 mHz以上頻率范圍呈現(xiàn)出一致性，精度量級(jí)相當(dāng).在30 mHz以上頻率范圍，計(jì)算得到的Vxx、Vyy、Vzz、Vxz的噪聲約為10 mE·Hz-1/2,Vxy、Vxz的噪聲約為200 mE·Hz-1/2，明顯體現(xiàn)了加速度計(jì)US軸和LS軸不同敏感度對(duì)重力梯度分量的影響，即Vxx、Vyy、Vzz和Vxz的測(cè)量精度較高，而Vxy和Vyz的測(cè)量精度較低.Vxx、Vyy、Vzz、Vxy對(duì)應(yīng)的平方根功率譜在0.02～5 mHz低頻率段內(nèi)，其計(jì)算分量精度略低于EGG_GGT中對(duì)應(yīng)分量(如圖11a—d所示).

圖11 計(jì)算所得梯度分量與L1b_GGT的PSD1/2Fig.11 PSD1/2 of the calculated gradient component and L1b_GGT

由表5可得，計(jì)算得到的Vxx與官方發(fā)布的GGT_Vxx相比較，其標(biāo)準(zhǔn)差分別為6791.54 mE、6714.88 mE，兩者相差約為76.7 mE；而Vyy、Vzz、Vxy的標(biāo)準(zhǔn)差均在20.4～34.2 mE范圍內(nèi).如圖11e、f所示，計(jì)算得到的Vxz、Vyz與GGT對(duì)應(yīng)分量在全頻段內(nèi)的精度呈現(xiàn)出極強(qiáng)的相似性，Vxz、Vyz與GGT對(duì)應(yīng)分量的標(biāo)準(zhǔn)差分別為11.5 mE、17.5 mE.測(cè)量帶寬范圍內(nèi)，重力梯度對(duì)角線分量在5 mHz處噪聲約為1 E·Hz-1/2，0.1 Hz處約為10 mE·Hz-1/2，并在20～40 mHz范圍內(nèi)降低到8 mE·Hz-1/2.圖12中計(jì)算得到梯度分量的跡與GGT相比，在低于10 mHz范圍內(nèi)受到的噪聲影響更大，約低于1個(gè)量級(jí).GOCE任務(wù)設(shè)計(jì)重力梯度張量跡在20～100 mHz內(nèi)精度為11 mE·Hz-1/2，本文計(jì)算結(jié)果在該范圍內(nèi)約為10 mE·Hz-1/2，達(dá)到設(shè)計(jì)精度.

表5 梯度各分量計(jì)算值與L1b_GGT的統(tǒng)計(jì)Table 5 Calculation values of gradient components and L1b_GGT statistics

圖12 計(jì)算所得梯度跡與L1b_GGT跡的PSD1/2Fig.12 PSD1/2 of calculated gradient trace and L1b_GGT trace

4 結(jié)論

本文針對(duì)重力梯度測(cè)量衛(wèi)星L1級(jí)數(shù)據(jù)處理，以歐洲空間局發(fā)布的GOCE衛(wèi)星數(shù)據(jù)資料為參考，初步實(shí)現(xiàn)了重力梯度測(cè)量衛(wèi)星L1級(jí)數(shù)據(jù)全流程構(gòu)建計(jì)算方法.建立了重力梯度儀觀測(cè)電壓數(shù)據(jù)到加速度數(shù)據(jù)轉(zhuǎn)換，提出了多星敏感器聯(lián)合法，有效地抑制了坐標(biāo)系變換導(dǎo)致的精度較低角速度分量噪聲傳播，重建了高精度衛(wèi)星姿態(tài)角速度，并最終構(gòu)建了符合精度要求的重力梯度分量.

(1)通過(guò)解算電壓增益因子得到的加速度觀測(cè)數(shù)據(jù)與歐空局官方數(shù)據(jù)在各軸標(biāo)準(zhǔn)差達(dá)到1.1×10-14～2.1×10-12m·s-2量級(jí)，在測(cè)量帶寬5～100 mHz內(nèi)，共模加速度計(jì)算精度達(dá)到了梯度儀超靈敏軸10-10～10-11m·s-2·Hz-1/2的精度要求.

(2)聯(lián)合多星敏感器姿態(tài)數(shù)據(jù)能有效抑制SSRF-GRF轉(zhuǎn)換過(guò)程中導(dǎo)致的噪聲傳播，角速度重建精度達(dá)到10-6rad/s量級(jí)，其中分量wy、wz相比聯(lián)合前約提升1個(gè)量級(jí)；基于維納濾波重建角速度方法，將EGG與STR的角速度根據(jù)其精度在頻域內(nèi)進(jìn)行聯(lián)合，在測(cè)量帶寬5～100 mHz內(nèi)角速度分量wy精度最大提升了6.56×10-11rad·s-1·Hz-1/2.

(3)構(gòu)建了重力梯度各分量，其中Vxx、Vyy、Vzz在5 mHz處約為1 E·Hz-1/2，0.1 Hz處約為10 mE·Hz-1/2，與歐空局發(fā)布結(jié)果的標(biāo)準(zhǔn)差相差為20.4～76.7 mE；在20～100 mHz頻率范圍內(nèi)，重力梯度張量的跡約為10 mE·Hz-1/2，符合官方精度要求，驗(yàn)證了本文構(gòu)建方法的有效性.

(4)本文初步系統(tǒng)實(shí)現(xiàn)了重力梯度測(cè)量衛(wèi)星L1級(jí)數(shù)據(jù)構(gòu)建計(jì)算方法，能夠?yàn)槲覈?guó)未來(lái)開(kāi)展重力衛(wèi)星任務(wù)提供自主的數(shù)據(jù)處理技術(shù)積累和科學(xué)參考.