衛星重力梯度觀測數據L1級構建方法

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

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

0 引言

全球重力場和海洋環流探測器(Gravity field and steady-state Ocean Circulation Explorer, GOCE)由歐洲空間局(European Space Agency，ESA)于2009年3月17日成功發射，在軌運行4年后于2013年11月11日結束衛星任務.GOCE衛星主要任務目標是在100 km的空間分辨率下確定精度為1～2 cm的全球大地水準面和精度為1 mGal的全球重力異常(ESA, 1999).為實現這一極具難度的科學目標，GOCE衛星搭載了高精度靜電重力梯度儀(Electrostatic Gravity Gradiometry，EGG)、高精度全球導航衛星系統(Global Navigation Satellite System, GNSS)星載接收機和星敏感器(Star Sensor，STR)等多個關鍵載荷(Cesare, 2008; Rummel et al., 2011).

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

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

本文面向我國發展的梯度衛星的任務需要，系統研究并初步實現了自主的衛星重力梯度觀測數據L1級構建方法，主要包括加速度計電壓數據轉換、多星敏感器聯合姿態數據的角速度重建、衛星重力梯度分量構建等數據處理技術環節.研究工作可為下一步我國推進實施民用重力梯度測量衛星任務提供自主的原始數據處理技術支撐與儲備.

1 衛星重力梯度觀測數據構建原理

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

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

(1)

圖1 梯度儀坐標系內梯度儀的6個加速度計位置Fig.1 Position of six accelerometers of gradiometer in gradiometer reference frame

GOCE衛星所搭載的6個加速度計都是通過剛性連接固定在梯度儀內，因此它們在衛星質心處的非保守力加速度相同，加速度計對距離矢量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級觀測數據構建方法

重力梯度測量衛星觀測數據L1級構建主要包括電壓數據轉換、姿態數據重建、角速度重建和重力梯度分量構建等四個步驟，具體計算方法如下：

圖2 單個加速度計電極結構Fig.2 Electrode structure of single accelerometer

2.1 加速度計電壓數據轉換

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

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

(11)

(12)

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

2.2 多星敏感器聯合重建姿態數據

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

(13)

2.3 角速度重建

衛星角速度與姿態四元數的精度直接影響地球重力場模型反演精度.星敏感器、重力梯度儀觀測數據確定角速度時，會出現兩者角速度噪聲分別在高頻、低頻增大的影響.為了精確地確定衛星在空間中的慣性角速度，其重建過程中應充分考慮到EGG與STR的誤差特性.基于Stummer等(2011)提出的維納濾波重建角速度方法，可將EGG與STR的角速度根據其精度在頻域內聯合.表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

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

(14)

(15)

所有頻率的權重之和等于1，由此可得到：

(16)

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

(17)

(18)

圖4 EGG和STR維納濾波權重Fig.4 EGG and STR Wiener filter weight

2.4 重力梯度數據構建

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

(19)

(20)

衛星星體在3個方向上所受到的非保守力可通過加速度計測量的共模加速度表示，其共模加速度存在以下關系：

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

(21)

即共模加速度間存在以下線性組合關系：

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)中得到重力梯度，應從差模加速度中分離出重力梯度.即V的主對角線重力梯度分量計算式為：

(23)

(24)

(25)

非對角線重力梯度分量計算式為：

(26)

(27)

(28)

式中，ωx、ωy、ωz為衛星的角速度.高精度的衛星角速度ωx、ωy和ωz為解算梯度V過程中的關鍵參數.

3 計算結果與分析

3.1 數據

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

3.2 加速度計電壓數據轉換

本文針對GOCE衛星6個加速度計(A1—A6)展開分析計算，通過對電壓數據轉換得到加速度數據，并與ESA發布的L1b數據進行對比分析.表3為解算的加速度計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.以加速度計A1計算結果為例，圖5為加速度計A1各個軸線性加速度與ESA發布的L1b數據中NLA_A1(Acceleration Nominal Linear, NLA)差值，圖6為加速度計A1的計算值Cal_A1與ESA中NLA_A1各軸線性加速度功率譜密度(Power Spectral Density，PSD).表4給出了A1計算值與NLA_A1各個軸加速度的統計特性.

從圖5可以看出，計算所得的A1線性加速度值與官方發布加速度數據，在x、y、z軸上分量ax,ay,az的差值均在非常小的范圍內.圖5中黑線為各軸差值的均值線，x、y、z軸對應均值分別為-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，其差值量級體現出很強的相似性.如圖6所示為A1各軸加速度功率譜密度，表明計算結果與NLA_A1中對應數據在功率譜密度上均呈現一致性(Welch, 1967)，其中ax、az分別在低頻段2～7 mHz、高頻段300～500 mHz內精度稍低，其余頻段各軸精度幾乎相同.表4中各軸角速度的統計特性也體現了計算結果與GOCE官方發布數據在各個軸上高度相似性，其標準差均非常接近，差異極小，達到了1.1×10-14～2.1×10-12m·s-2量級，精度結果十分相近.

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

表4 加速度計A1與NLA_A1各個軸加速度的統計Table 4 Statistics of triaxial acceleration of accelerometer A1 and NLA_A1

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

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

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

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

(29)

(30)

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

σa_comb≈σa_US.

(31)

圖7 x、y、z軸上加速度計對差分(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軸加速度的真實噪聲，在測量帶寬5～100 mHz內各US軸噪聲達到了10-10～10-11m·s-2·Hz-1/2范圍內的良好指標，符合GOCE衛星重力梯度儀的設計要求.

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

3.3 角速度重建

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

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

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

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

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

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

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

表5 梯度各分量計算值與L1b_GGT的統計Table 5 Calculation values of gradient components and L1b_GGT statistics

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

4 結論

本文針對重力梯度測量衛星L1級數據處理，以歐洲空間局發布的GOCE衛星數據資料為參考，初步實現了重力梯度測量衛星L1級數據全流程構建計算方法.建立了重力梯度儀觀測電壓數據到加速度數據轉換，提出了多星敏感器聯合法，有效地抑制了坐標系變換導致的精度較低角速度分量噪聲傳播，重建了高精度衛星姿態角速度，并最終構建了符合精度要求的重力梯度分量.

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

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

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

(4)本文初步系統實現了重力梯度測量衛星L1級數據構建計算方法，能夠為我國未來開展重力衛星任務提供自主的數據處理技術積累和科學參考.