陳秋杰，沈云中＊，張興福

1 同濟(jì)大學(xué)測繪與地理信息學(xué)院，上海 200092

2 同濟(jì)大學(xué)空間信息與可持續(xù)發(fā)展應(yīng)用中心，上海 200092

3 廣東工業(yè)大學(xué)測繪工程系，廣州 510006

1 引 言

隨著低軌重力衛(wèi)星CHAMP、GRACE和GOCE衛(wèi)星陸續(xù)發(fā)射及相關(guān)數(shù)據(jù)的釋放，國際相關(guān)研究機(jī)構(gòu)已使用該類數(shù)據(jù)反演了多個高精度的地球重力場模型.其主要代表是GFZ（Geo Forschungs Zentrum Potsdam）、CSR（Center for Space Research）、伯爾尼大學(xué)、波恩大學(xué)所分別反演的 EIGEN［1－2］、GGM［3－4］、AIUB［5］和 ITG－GRACE［6－8］系 列 的 重 力場模型以及多家機(jī)構(gòu)聯(lián)合（慕尼黑工業(yè)大學(xué)、波恩大學(xué)、格拉茨工業(yè)大學(xué)、奧地利科學(xué)院和伯爾尼大學(xué)）所反演的GOCO［9］系列重力場模型，其中，EIGEN和GGM模型均是采用動力學(xué)法，ITG－GRACE模型采用短弧長積分法，AIUB模型采用天體力學(xué)方法，GOCO模型是聯(lián)合GRACE和GOCE以及其他重力數(shù)據(jù)所反演的地球重力場模型.動力學(xué)法基于衛(wèi)星初始狀態(tài)向量、先驗地球重力場模型及其他攝動力模型計算出參考軌道后，以參考軌道為初值對非線性衛(wèi)星運動方程進(jìn)行線性化［10－11］，涉及變分方程積分計算，運算量較大.天體力學(xué)法同樣要用先驗地球重力場模型計算參考軌道，它與動力學(xué)法的區(qū)別在于采用了不同的參數(shù)化方法，其衛(wèi)星狀態(tài)采用開普勒軌道根數(shù)，并增加了經(jīng)驗加速度參數(shù)和其他動力學(xué)參數(shù)［5］.短弧長積分法將弧段的軌道向量表示成邊界軌道向量、重力場待估參數(shù)和其他力模型參數(shù)相關(guān)的積分方程［6－8，12］，其線性化的參考軌道由邊界軌道參數(shù)和已知力模型計算，為減小線性化誤差，還可直接用作了梯度改正的幾何觀測軌道［7］.Ditmar利用三點數(shù)值微分公式，根據(jù)3個相鄰點的衛(wèi)星軌道觀測值計算衛(wèi)星加速度，并將求得的加速度表示成3點軌道弧段上衛(wèi)星加速度的加權(quán)平均值，建立待估重力場參數(shù)與平均加速度之間的關(guān)系進(jìn)行求解［13］，但其解算模型并沒有以幾何軌道作為初值進(jìn)行線性化.Xu（2011）給出的模型直接以衛(wèi)星觀測軌道進(jìn)行線性化［14］，但該模型將線性化軌道參數(shù)的改正量重新表示成初始狀態(tài)參數(shù)和重力場參數(shù)的函數(shù)，這需要進(jìn)行4重積分.由于GNSS（Global Navigation Satellite System）確定的低軌重力衛(wèi)星幾何軌道精度可達(dá)1～2cm，因此直接以幾何軌道進(jìn)行線性化時，其二階以上的改正量可以忽略.考慮到線性化的一階改正量完全與定軌精度有關(guān)，本文將該改正量作為間接觀測值的誤差直接由最小二乘求解，除了求解的重力場模型參數(shù)外，只需要考慮非保守力加速度計觀測值的尺度和偏差參數(shù).與天體力學(xué)法和動力學(xué)法不同，本文方法不需要衛(wèi)星的初始狀態(tài)參數(shù)、初始重力場模型和解算變分方程，因此解算模型非常簡單、便于程序?qū)崿F(xiàn).Mayer－Gürr等（2005）采用短弧長積分法解算ITG－CHAMP01模型時，直接用衛(wèi)星軌道計算重力場參數(shù)的系數(shù)矩陣［6］，其本質(zhì)是忽略軌道誤差對系數(shù)矩陣影響的一種線性化方案，系數(shù)矩陣誤差將被觀測誤差吸收，影響模型的解算精度.因此，Mayer－Gürr（2006）對軌道進(jìn)行梯度改正［7］，減小其誤差對系數(shù)矩陣影響后，才用于GRACE高精度的低低衛(wèi)衛(wèi)跟蹤數(shù)據(jù)的解算.本文所提出方法本質(zhì)上是對短弧長積分法的改進(jìn)，不需要梯度改正就能用GRACE低低衛(wèi)衛(wèi)跟蹤數(shù)據(jù)解算地球重力場模型.利用2008年全年的2顆GRACE衛(wèi)星軌道數(shù)據(jù)解算了90階次的地球重力場模型TJGRACE01S，并與其他模型進(jìn)行比較，結(jié)果顯示，該模型精度優(yōu)于同階次的EIGENCHAMP03S和EIGEN－CHAMP05S模型，其前27階位系數(shù)整體精度甚至優(yōu)于EIGEN－GRACE01S，其前15階位系數(shù)整體精度與EIGEN－GRACE02S模型精度大致相當(dāng).

2 觀測與解算模型

對于由m＋1個離散點組成的衛(wèi)星軌道弧段，其相鄰點坐標(biāo)與加速度間的關(guān)系可表示為［15］

其中，r為衛(wèi)星坐標(biāo)，ti為第i歷元時刻，Δt為積分步長，等于相鄰離散點的時間間隔，a（tj）為衛(wèi)星加速度，等于單位質(zhì)量作用力，kj為Cowell積分器系數(shù)，n為積分器長度，J（i）定義如下：

衛(wèi)星的加速度a（tj）與衛(wèi)星位置r（tj）、重力位系數(shù)u和加速度計尺度和偏差參數(shù)w有關(guān).考慮到

其中，ag（tj）為保守力加速度，af（tj）為非保守力加速度.若直接以衛(wèi)星軌道觀測量rg（tj）以及重力場位系數(shù)和加速度計尺度和偏差參數(shù)的先驗值u0，w0為初值對ag（tj）和af（tj）進(jìn)行線性化，則

其中，δμ，δw分別為重力位系數(shù)參數(shù)和加速度計尺度與偏差參數(shù)的改正數(shù)，vr（tj）為軌道觀測值改正數(shù)，

順便指出，由于利用GNSS定軌結(jié)果rg（tj）作為初始軌道，引力向量ag（tj）是引力位系數(shù)的線性函數(shù)，因此系數(shù)矩陣H（tj）與先驗位系數(shù)u0無關(guān)，即使沒有先驗重力場模型（即u0＝0）時，展開式（3）依然成立.將（3）、（4）式代入（2）式，得

其中

將（5）代入（1），則有

對于有m＋1個觀測值的軌道弧段，可建立m－n＋1個觀測方程，將其表示為

相對于傳統(tǒng)動力學(xué)法，觀測模型（7）的主要特點如下：

（1）將位系數(shù)直接在重力衛(wèi)星觀測軌道處展開，而不需表示為衛(wèi)星初始狀態(tài)參數(shù)及變分方程解算的形式.

（2）由于H（tj）與先驗重力場模型的位系數(shù)參數(shù)無關(guān)，因此（7）式是位系數(shù)參數(shù)的線性觀測方程，引入先驗重力場模型僅改變常數(shù)項，其解算結(jié)果不受影響.

如果共有K個弧段，每個弧段k（k＝1，2，…，K）都可以組成如（7）式所示的觀測方程：

其中，下標(biāo)k表示與該弧段相關(guān)的量，其他含義與（7）相同.若Qk為觀測向量lk的協(xié)因數(shù)陣，wk為第k弧段加速度計尺度和偏差參數(shù)，則第k弧段法方程為

（8）式法方程可分塊表示為

位系數(shù)改正量δu為全局參數(shù)，加速度計尺度和偏差參數(shù)改正量δwk為局部參數(shù)，Nukuk為位系數(shù)參數(shù)在該弧段的法方程塊，Nwkwk為尺度和偏差參數(shù)的法方程塊，Nukwk和Nwkuk為全局參數(shù)與局部參數(shù)間聯(lián)系數(shù)的法方程塊，luk和lwk為相應(yīng)的常數(shù)項.考慮到法方程的可疊加性，消去每個弧段的局部參數(shù)后，得全局參數(shù)的法方程［16］

將每一弧段的全局參數(shù)法方程進(jìn)行疊加，得到求解位系數(shù)的總法方程式：

由（11）式可求得位系數(shù)參數(shù)的改正量.由于本文方法直接以幾何觀測軌道為初值進(jìn)行線性化，不需要初始狀態(tài)向量參數(shù)，因此避免了初始狀態(tài)向量參數(shù)與重力場模型參數(shù)相關(guān)性的影響.但本文的誤差方程（7）式的殘差向量具有系數(shù)矩陣B，導(dǎo)致形成（8）式的法矩陣和常數(shù)項時，均需的逆陣，該矩陣的階數(shù)約為弧段歷元數(shù)的3倍，因此考慮到高維逆矩陣求解的數(shù)值穩(wěn)定性與解算效率，弧段長度不宜設(shè)置得過長.同時，采用Cholesky矩陣分解方法進(jìn)行求逆，可以提高逆矩陣求解的穩(wěn)定性，相應(yīng)算法可參考相關(guān)文獻(xiàn)［17－19］.

3 反演弧段長度

為了選擇合適的弧段長度，采用官方JPL（Jet Propulsion Laboratory）所公布的2008年1月份5s采樣率的RL02約化動力學(xué)軌道數(shù)據(jù)（RL02版本軌道計算采用GIF48為先驗地球重力場模型）和1s采樣率的加速度計數(shù)據(jù)，軌道精度約為2cm，加速度計數(shù)據(jù)采樣率進(jìn)一步壓縮為5s.分別采用不同弧段長度，解算60階次的地球重力場模型.為減小重力場高階信號截斷對位系數(shù)解算精度的影響，第61階至150階采用ITG－GRACE2010模型的重力場位系數(shù)，要反演的前60階系數(shù)不用先驗重力場模型.ITG－GRACE2010模型為波恩大學(xué)采用7.5年GRACE低低衛(wèi)衛(wèi)跟蹤數(shù)據(jù)反演的靜態(tài)地球重力場模型.衛(wèi)星軌道、加速度及姿態(tài)觀測數(shù)據(jù)先要進(jìn)行預(yù)處理，包括粗差探測、數(shù)據(jù)內(nèi)插、數(shù)據(jù)間斷處理、坐標(biāo)系轉(zhuǎn)換等.考慮了日月引力、海潮、固體潮、極潮、大氣與海洋的非潮汐變化影響等保守力和加速度計測定的非保守力（具體信息見表1），非保守力的尺度和偏差參數(shù)先驗值采用GFZ公布值，每個弧段估計一組尺度和偏差參數(shù).形成總法方程后采用Cholesky矩陣分解法進(jìn)行重力位系數(shù)的解算，根據(jù)解算結(jié)果與EGM2008模型的差值，進(jìn)行相關(guān)分析.不同弧段長度求得的位系數(shù)與EGM2008模型位系數(shù)差值的階方差和相應(yīng)的大地水準(zhǔn)面累積誤差分別如圖1和圖2所示，其結(jié)果表明，30min弧段長度的位系數(shù)差值的階方差和大地水準(zhǔn)面累積誤差最小，因此實際數(shù)據(jù)處理時，采用30min弧段長度進(jìn)行地球重力場反演.

表1 攝動力模型信息Table 1 The information of perturbation force models

4 反演結(jié)果與分析

實際數(shù)據(jù)處理采用 JPL（Jet Propulsion Laboratory）官方公布的2008年全年的RL02版本的GRACE約化動力學(xué)軌道和加速度計數(shù)據(jù).根據(jù)上一節(jié)結(jié)果，我們采用30min的弧段長度反演90階次地球重力場模型.其數(shù)據(jù)預(yù)處理方法、保守力與非保守力計算模型如前所述，為減小重力場高階信號截斷對重力場解算精度的影響，第91階至150階的重力場位系數(shù)采用高精度的ITG－GRACE2010模型系數(shù)，前90階重力位系數(shù)的初值設(shè)置為零，迭代解算了90階次的地球重力場模型TJGRACE01S，其大地水準(zhǔn)面累積誤差為17.6cm.為了分析比較TJGRACE01S模型的精度，本文給出TJGRACE01S、EIGEN－CHAMP03S、EIGEN－CHAMP05S、EIGENGRACE01S和EIGEN－GRACE02S模型與EGM2008模型差值的階方差如圖3所示，相應(yīng)的大地水準(zhǔn)面累積誤差如圖4所示.

其中，EIGEN－CHAMP03S、EIGEN－CHAMP05S分別為GFZ采用3年、6年的CHAMP衛(wèi)星軌道解算的地球重力場模型，EIGEN－GRACE01S、EIGENGRACE02S分別為GFZ采用39天、110天的低低衛(wèi)星跟蹤數(shù)據(jù)解算的地球重力場模型.圖3和圖4表明，采用一年的GRACE軌道所解算出的重力場模型TJGRACE01S的精度比使用多年CHAMP軌道數(shù)據(jù)所反演的EIGEN－CHAMP03S和EIGENCHAMP05S模型都要高.從大地水準(zhǔn)面累積誤差看，前27階的TJGRACE01S模型整體精度比EIGEN－GRACE01S模型的高；從位系數(shù)差值階方差看，前15階的TJGRACE01S模型精度與EIGENGRACE02S模型大致相當(dāng).顯然TJGRACE01S模型的低階位系數(shù)精度比較高，其高階位系數(shù)精度不如 EIGEN－GRACE01S和 EIGEN－GRACE02S，主要原因為GRACE軌道數(shù)據(jù)適合反演重力場的長波長位系數(shù)，而中階位系數(shù)主要依靠GRACE KBR觀測數(shù)據(jù).

為了進(jìn)一步分析TJGRACE01S模型的精度與可靠程度，選擇TJGRACE01S、EIGEN－CHAMP03S、EIGEN－CHAMP05S、EIGEN－GRACE01S、EIGENGRACE02S和EGM2008模型，分別取不同階次計算區(qū)域：緯度－90°～90°，經(jīng)度－180°～180°，間隔為1°的網(wǎng)格點的大地水準(zhǔn)面高.以EGM2008模型計算得到的大地水準(zhǔn)面高為基準(zhǔn)值，其他模型求得的結(jié)果與基準(zhǔn)值之差的統(tǒng)計結(jié)果如表2所示.由表2可見，對于不同地球重力場模型，截斷至相同階次（10、30、50、70、90），從所計算的大地水準(zhǔn)面高差值的標(biāo)準(zhǔn)差來看，TJGRACE01S模型的精度都比相同階次的EIGEN－CHAMP03S和EIGEN－CHAMP05S模型的高.

為了檢驗?zāi)Ｐ偷耐夥暇?，選擇美國8221個A級GPS水準(zhǔn)點的實測高程異常作為基準(zhǔn)，將EIGEN－CHAMP03S、 EIGEN－CHAMP05S、 EIGENGRACE01S、EIGEN－GRACE02S 和 TJGRACE01S模型所解算的這些點的高程異常與基準(zhǔn)值進(jìn)行比較.為避免各模型的譜泄露，采用譜組合技術(shù)，即被分析模型所選擇階次以外的位系數(shù)均用EGM2008模型擴(kuò)充到2160階次，如被分析模型選擇階次為N，則N＋1階至2160階用EGM2008位系數(shù)補(bǔ)充，得到各模型的外符合精度如表3所示.表3結(jié)果表明，對于不同地球重力場模型，截斷至相同階次（10、30、50、70、90）并與 EGM2008模型進(jìn)行譜組合，從所計算的模型高程異常與實測高程異常差的標(biāo)準(zhǔn)差來看，除了階次截斷為70外，TJGRACE01S模型所計算的標(biāo)準(zhǔn)差均小于EIGEN－CHAMP03S和EIGEN－CHAMP05S模型，因此可以說TJGRACE01S模型的精度總體高于EIGEN－CHAMP03S和EIGEN－CHAMP05S模型，驗證了 TJGRACE01S模型的可靠性以及本文算法的有效性.

表2 不同模型大地水準(zhǔn)面高度差統(tǒng)計結(jié)果Table 2 The statistical results of geoid height difference between different models

5 結(jié) 論

本文提出的基于重力衛(wèi)星幾何軌道線性化的地球重力場反演方法不需要衛(wèi)星初始軌道參數(shù)及先驗地球重力場模型和變分方程解算，方法簡單，解算效率高.利用2008年全年的GRACE雙星衛(wèi)星軌道數(shù)據(jù)，采用本文方法解算了90階次的地球重力場模型TJGRACE01S，結(jié)果表明，TJGRACE01S模型在90階次的大地水準(zhǔn)面累積誤差為17.6cm，其精度優(yōu)于同階次的EIGEN－CHAMP03S和EIGEN－CHAMP05S模型；從大地水準(zhǔn)面累積誤差看，前27階位系數(shù)精度優(yōu)于EIGEN－GRACE01S；從位系數(shù)差值階方差看，前15階位系數(shù)精度與EIGEN－GRACE02S模型精度大致相當(dāng).利用EGM2008模型和美國8221個GPS水準(zhǔn)數(shù)據(jù)進(jìn)一步驗證了本文方法的有效性.

表3 各模型確定的高程異常與美國實測GPS／水準(zhǔn)結(jié)果的比較Table 3 The statistic results of the height anomaly differences between GPS／leveling of the USA and the height anomaly from all models

TJGRACE01S模型僅僅采用1年的GRACE雙星衛(wèi)星軌道數(shù)據(jù)，若采用更長時間的觀測數(shù)據(jù)，并與GRACE KBR數(shù)據(jù)進(jìn)行融合處理，TJGRACE01S模型的精度將會更進(jìn)一步得到提高.

（References）

［1］Reigber C，Schmidt R，F(xiàn)lechtner F，et al.First GFZ GRACE gravity field model IGEN－GRACE01S.http：／／op.gfz－potdam.de／grace／results.

［2］Reigber C，Schmidt R，F(xiàn)lechtner F，et al.An Earth gravity field model complete to degree and order 150from GRACE EIGEN－GRACE02S.Journal of Geodynamics，2005，39（1）：1－10.

［3］Tapley B，Ries J，Bettadpur S，et al.GGM02－An improved Earth gravity field model from GRACE.Journal of Geodesy，2005，79（8）：467－478.

［4］Chambers D P.Evaluation of new GRACE time－variable gravity data over the ocean.Geophys.Res.Lett.，33，L17603，doi：10.1029／2006GL027296.

［5］Beutler G，Jaggi A， Mervart L，et al.The celestial mechanics approach：theoretical foundations.Journal of Geodesy，2010，84（10）：605－624.

［6］Mayer－Gürr T，Ilk KH，Eicker A，et al.ITG－CHAMP01：a CHAMP gravity field model from short kinematic arcs over a one－year observation period.Journal of Geodesy，2005，78（7－8）：462－480.

［7］Mayer－Gürr T.Gravitationsfeldbestimmung aus der Analyse kurzer Bahnboegen am Beispiel der Satellitenmissionen CHAMP und GRACE.Dissertation at the Institute of Theoretical Geodesy，University Bonn，2006.

［8］Mayer－Gürr T，Eicker A，Kurtenbach E，et al.ITGGRACE：Global Static and Temporal Gravity Field Models from GRACE Data.Earth and Environmental Science，2010，（Part 2）：159－168.

［9］Pail R，Bruinsma S，Migliaccio F，et al.First GOCE gravity field models derived by three different approaches.Journal of Geodesy，2011，85（11）：819－843.

［10］沈云中.應(yīng)用CHAMP衛(wèi)星星歷精化地球重力場模型的研究.武漢：中國科學(xué)院測量與地球物理研究所，2000.Shen Y Z.Study of recovering gravitational potential model from the ephemerides of CHAMP（in Chinese）.Wuhan：The Institute of Geodesy and Geophysics，Chinese Academy of Sciences，2000.

［11］周旭華，吳斌，許厚澤等.數(shù)值模擬估算低低衛(wèi)衛(wèi)跟蹤觀測技術(shù)反演地球重力場的空間分辨率.地球物理學(xué)報，2005，48（2）：282－287.Zhou X H，Wu B，Xu H Z，et al.Resolution estimation of Earth gravity field recovery through the low－low satellite to satellite technology by numerical simulation.Chinese J.Geophys.（in Chinese），2005，48（2）：282－287.

［12］Ilk K H， Mayer－Gürr T，F(xiàn)euchtinger M.Gravity field recovery by analysis of short arcs of CHAMP.Earth Observation with CHAMP，2005：127－132.

［13］Ditmar P，Kuznetsov V，van Eck van der Sluijs A A，et al.＇DEOS＿CHAMP－01C＿70＇：a model of the Earth＇s gravity field computed from accelerations of the CHAMP satellite.Journal of Geodesy，2006，79（10－11）：586－601.

［14］Xu P L. Position and velocity perturbations for the determination of geopotential from space Geodetic measurements.Celestial Mechanics and Dynamical Astronomy，2008，100（3）：231－249.

［15］王正濤.衛(wèi)星跟蹤衛(wèi)星測量確定地球重力場的理論與方法.武漢：武漢大學(xué)，2005.Wang Z T.Theory and methodology of Earth gravity field recovery by satellite to satellite tracking data （in Chinese）.Wuhan：Wuhan University，2005.

［16］張興福.應(yīng)用低軌衛(wèi)星跟蹤數(shù)據(jù)反演地球重力場模型.上海：同濟(jì)大學(xué)，2007.Zhang X F.The earth′s field model recovery on the basis of satellite－to－satellite tracking missions （in Chinese）.Shanghai：Tongji University，2007.

［17］Davis T A，Hager W W.Dynamic supernodes in sparse Cholesky update／downdate and triangular solves.Math Software，2009，35（4）：1－17.

［18］Chen Y，Davis T A，Hager W W，et al.Algorithm 887：CHOLMOD，supernodal sparse Cholesky factorization and update／downdate.Math Software，2009，35（3）：2－14.

［19］Davis T A，Hager W W.Row modifications of a sparse Cholesky factorization.SIAM Journal on Matrix Analysis and Applications，2005，26（3）：621－639.