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

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

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

0 引言

全球重力場和海洋環(huán)流探測器(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)等多個關(guān)鍵載荷(Cesare, 2008; Rummel et al., 2011).

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

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

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

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

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

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

(1)

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

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

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

(2)

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

(3)

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

(4)

(5)

利用V與Ω2的對稱性可得：

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

(6)

(7)

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

(8)

(9)

(10)

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

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

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

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

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

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

(11)

(12)

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

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

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

(13)

2.3 角速度重建

衛(wèi)星角速度與姿態(tài)四元數(shù)的精度直接影響地球重力場模型反演精度.星敏感器、重力梯度儀觀測數(shù)據(jù)確定角速度時，會出現(xiàn)兩者角速度噪聲分別在高頻、低頻增大的影響.為了精確地確定衛(wèi)星在空間中的慣性角速度，其重建過程中應(yīng)充分考慮到EGG與STR的誤差特性.基于Stummer等(2011)提出的維納濾波重建角速度方法，可將EGG與STR的角速度根據(jù)其精度在頻域內(nèi)聯(lián)合.表1、表2分別為EGG和STR的噪聲模型，STR和EGG三個角速度分量具有不同的噪聲功率頻譜密度，如圖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)重.由于維納濾波是在時域中聯(lián)合角速度，濾波系數(shù)可通過離散傅里葉逆變換(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ì)觀測數(shù)據(jù)和角速度觀測數(shù)據(jù)，綜合式(6)和(7)：

(19)

(20)

衛(wèi)星星體在3個方向上所受到的非保守力可通過加速度計(jì)測量的共模加速度表示，其共模加速度存在以下關(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)

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

(23)

(24)

(25)

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

(26)

(27)

(28)

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

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

3.1 數(shù)據(jù)

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

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

本文針對GOCE衛(wèi)星6個加速度計(jì)(A1—A6)展開分析計(jì)算，通過對電壓數(shù)據(jù)轉(zhuǎn)換得到加速度數(shù)據(jù)，并與ESA發(fā)布的L1b數(shù)據(jù)進(jìn)行對比分析.表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各個軸線性加速度與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各個軸加速度的統(tǒng)計(jì)特性.

從圖5可以看出，計(jì)算所得的A1線性加速度值與官方發(fā)布加速度數(shù)據(jù)，在x、y、z軸上分量ax,ay,az的差值均在非常小的范圍內(nèi).圖5中黑線為各軸差值的均值線，x、y、z軸對應(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的量級分別為10-8m·s-2，10-6m·s-2，10-7m·s-2，其差值量級體現(xiàn)出很強(qiáng)的相似性.如圖6所示為A1各軸加速度功率譜密度，表明計(jì)算結(jié)果與NLA_A1中對應(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ù)在各個軸上高度相似性，其標(biāo)準(zhǔn)差均非常接近，差異極小，達(dá)到了1.1×10-14～2.1×10-12m·s-2量級，精度結(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各個軸加速度的統(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對比Fig.6 PSD1/2 comparison between triaxial acceleration calculation results and NLA_A1

各個加速度計(jì)觀測的加速度值通過式(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量級，其余各差分加速度約在10-9m·s-2·Hz-1/2量級.共模加速度呈現(xiàn)出隨著頻率升高其精度越高的特點(diǎn)，各軸共模加速度在測量帶寬(Measurement Band Width, MBW)范圍內(nèi)(5～100 mHz)趨勢一致且精度相似.

重力梯度儀由6個三軸加速度計(jì)組成，每個加速度計(jì)都有2個超靈敏軸(Ultra Sensitive, US)和1個低靈敏度軸(Less Sensitive, LS).由于所有6個加速度計(jì)超靈敏軸均固定安裝在GRF中的x方向上(如圖1)，對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為單個US軸標(biāo)準(zhǔn)差.將式(30)代入式(29)可得：

σa_comb≈σa_US.

(31)

圖7 x、y、z軸上加速度計(jì)對差分(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í)噪聲，在測量帶寬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可以看出，與單個STR相比，多個STR聯(lián)合后顯示出y、z軸的角速度分量wy、wz的總體精度提升明顯；在10～100 mHz處逐步下降了約1個量級，精度達(dá)到10-6rad·s-1·Hz-1/2量級；在3～30 mHz內(nèi)wy的精度最高，而在30 mHz以上wx精度最高，wy與wz精度幾乎相同，且wx與wy、wz均呈現(xiàn)出精度變化趨勢相似.通過對圖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角速度的精度有明顯提升.在測量帶寬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個梯度分量的平方根功率譜密度均約在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)信號噪聲的影響.STR姿態(tài)噪聲以每轉(zhuǎn)周期(cycle per revolution, cpr)的頻率(1 cpr≈1.85×10-4Hz)表現(xiàn)出明顯的重復(fù)性，即在kcpr頻率處(k為整數(shù))6個梯度分量精度均受到的影響最大，并且在相應(yīng)的頻率下向較高的頻率呈1/f幅度減小.圖11中6個重力梯度分量中對應(yīng)的平方根功率譜密度均顯示出整體較為吻合，趨勢較為一致.在5 mHz以上頻率范圍呈現(xiàn)出一致性，精度量級相當(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軸不同敏感度對重力梯度分量的影響，即Vxx、Vyy、Vzz和Vxz的測量精度較高，而Vxy和Vyz的測量精度較低.Vxx、Vyy、Vzz、Vxy對應(yīng)的平方根功率譜在0.02～5 mHz低頻率段內(nèi)，其計(jì)算分量精度略低于EGG_GGT中對應(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對應(yīng)分量在全頻段內(nèi)的精度呈現(xiàn)出極強(qiáng)的相似性，Vxz、Vyz與GGT對應(yīng)分量的標(biāo)準(zhǔn)差分別為11.5 mE、17.5 mE.測量帶寬范圍內(nèi)，重力梯度對角線分量在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個量級.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é)論

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

(1)通過解算電壓增益因子得到的加速度觀測數(shù)據(jù)與歐空局官方數(shù)據(jù)在各軸標(biāo)準(zhǔn)差達(dá)到1.1×10-14～2.1×10-12m·s-2量級，在測量帶寬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)換過程中導(dǎo)致的噪聲傳播，角速度重建精度達(dá)到10-6rad/s量級，其中分量wy、wz相比聯(lián)合前約提升1個量級；基于維納濾波重建角速度方法，將EGG與STR的角速度根據(jù)其精度在頻域內(nèi)進(jìn)行聯(lián)合，在測量帶寬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)了重力梯度測量衛(wèi)星L1級數(shù)據(jù)構(gòu)建計(jì)算方法，能夠?yàn)槲覈磥黹_展重力衛(wèi)星任務(wù)提供自主的數(shù)據(jù)處理技術(shù)積累和科學(xué)參考.