衛(wèi)星精密定軌與重力場(chǎng)建模的同解法

鄒賢才， 衷路萍， 李建成

武漢大學(xué)測(cè)繪學(xué)院，地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室， 武漢　430079

?

衛(wèi)星精密定軌與重力場(chǎng)建模的同解法

鄒賢才， 衷路萍， 李建成

武漢大學(xué)測(cè)繪學(xué)院，地球空間環(huán)境與大地測(cè)量教育部重點(diǎn)實(shí)驗(yàn)室， 武漢430079

摘要現(xiàn)代衛(wèi)星跟蹤衛(wèi)星重力測(cè)量技術(shù)顯著改善了地球重力場(chǎng)模型的中長波段信號(hào)，并拓展了重力場(chǎng)模型在相關(guān)科學(xué)研究中的應(yīng)用.同解法作為衛(wèi)星重力觀測(cè)數(shù)據(jù)的主要處理手段之一，國內(nèi)一直沒有實(shí)質(zhì)突破.本文從基本模型和關(guān)鍵技術(shù)的分析出發(fā)，剖析了同解法的特點(diǎn)，特別是在建立同解法與幾何法(運(yùn)動(dòng)學(xué))定軌、一般動(dòng)力學(xué)方法關(guān)系的基礎(chǔ)上，給出了一種同解法的實(shí)現(xiàn)路線.在已有的精細(xì)數(shù)據(jù)預(yù)處理和并行計(jì)算研究基礎(chǔ)上，結(jié)合GRACE衛(wèi)星的實(shí)測(cè)飛行數(shù)據(jù)，在國內(nèi)首次獲得了真實(shí)衛(wèi)星任務(wù)數(shù)據(jù)條件下的同解法結(jié)果，并進(jìn)行了動(dòng)力法軌道的外部質(zhì)量檢核、衛(wèi)星非保守力分析、重力場(chǎng)模型的GPS水準(zhǔn)檢驗(yàn)等.利用衛(wèi)星激光測(cè)距數(shù)據(jù)檢驗(yàn)，衛(wèi)星精密軌道的徑向精度優(yōu)于2 cm，同時(shí)建立了質(zhì)量可靠的衛(wèi)星重力場(chǎng)模型，充分展示了同解法的優(yōu)點(diǎn).數(shù)值結(jié)果及其分析表明，本文所提的同解法實(shí)施方案合理可行，已經(jīng)掌握了實(shí)現(xiàn)同解法的關(guān)鍵技術(shù)，獲得了從仿真研究到實(shí)際飛行數(shù)據(jù)處理的新進(jìn)展.最后，本文對(duì)同解法今后的發(fā)展思路，以及如何進(jìn)一步挖掘同解法的潛力，提出了見解和今后的工作方向.

關(guān)鍵詞衛(wèi)星重力； 同解法； 加速度計(jì)校準(zhǔn)； GRACE； 動(dòng)力法

1引言

進(jìn)入21世紀(jì)以來，歐美發(fā)達(dá)國家先后實(shí)施了CHAMP(CHAllenging Minisatellite Payload)、GRACE(Gravity Recovery and Climate Experiment)和GOCE(Gravity field and steady-state Ocean Circulation Explorer)衛(wèi)星重力測(cè)量任務(wù).在三個(gè)新一代衛(wèi)星重力計(jì)劃中，衛(wèi)星跟蹤衛(wèi)星技術(shù)得到廣泛使用，衛(wèi)星跟蹤衛(wèi)星觀測(cè)數(shù)據(jù)成為解算重力場(chǎng)模型的主要數(shù)據(jù)源之一，包含高-低和低-低模式的衛(wèi)星跟蹤衛(wèi)星數(shù)據(jù).應(yīng)用的基本步驟則是首先利用衛(wèi)星重力觀測(cè)數(shù)據(jù)確定全球重力場(chǎng)模型，然后以該模型為基礎(chǔ)，開展相關(guān)的科學(xué)研究和工程應(yīng)用.因此，如何利用衛(wèi)星跟蹤衛(wèi)星觀測(cè)數(shù)據(jù)解算高精度高分辨率的衛(wèi)星重力場(chǎng)模型成為衛(wèi)星重力測(cè)量任務(wù)中的關(guān)鍵環(huán)節(jié).

在衛(wèi)星重力場(chǎng)建模領(lǐng)域，有能量法(Wolff,1969; Jekeli, 1999)、動(dòng)力法(Reigber, 1989)(指一般意義上的兩步法)、加速度法(Reubelt et al., 2003)和短弧法(Mayer et al., 2005)等.這些方法的共同點(diǎn)是需要已知重力衛(wèi)星的精密軌道，然后由衛(wèi)星精密軌道，通過軌道攝動(dòng)原理解算重力場(chǎng)模型.與此相對(duì)，同解法是“一步法”，直接利用原始跟蹤數(shù)據(jù)同時(shí)確定重力衛(wèi)星的精密軌道以及衛(wèi)星重力場(chǎng)模型.同解法使觀測(cè)數(shù)據(jù)和觀測(cè)模型之間達(dá)到了更高程度的“同化”，理論上更加嚴(yán)密，可實(shí)現(xiàn)更優(yōu)的參數(shù)估計(jì).國際上的著名地學(xué)研究機(jī)構(gòu)GFZ(GeoForschungsZentrum)，就是通過直接處理原始的跟蹤數(shù)據(jù)，同時(shí)解算軌道和重力場(chǎng)模型，該方法較一般的“兩步法”，觀測(cè)模型更加嚴(yán)密，解算的精度也更高(Zhu et al., 2004).

同解法實(shí)現(xiàn)復(fù)雜，其具體技術(shù)路線并未見諸公開刊物，在這個(gè)特定領(lǐng)域，國內(nèi)一直處于跟蹤狀態(tài).目前，在衛(wèi)星精密定軌技術(shù)方面，國內(nèi)已經(jīng)開展了卓有成效的工作，在運(yùn)動(dòng)學(xué)定軌、簡化的動(dòng)力學(xué)定軌、自適應(yīng)定軌等方面，取得了較為全面的成果.國內(nèi)展開的定軌方法研究，覆蓋面廣、技術(shù)路線全面，包括了非差、單差等多種技術(shù)模式，徑向精度總體上已經(jīng)接近3～5 cm的水平.這些工作使得我國在低軌衛(wèi)星精密定軌的數(shù)據(jù)處理領(lǐng)域，達(dá)到了國際先進(jìn)水平(鄭作亞，2004；彭冬菊和吳斌，2007；趙齊樂等，2008；李建成等，2009；秦顯平，2009).在衛(wèi)星重力場(chǎng)建模方面，國內(nèi)多個(gè)單位采用“兩步法”研究了基于衛(wèi)星跟蹤衛(wèi)星技術(shù)的重力場(chǎng)建模(周旭華，2005；周旭華等，2006；肖云，2006；鄒賢才，2007；張興福，2007).在后續(xù)的模型研制上也取得了長足的進(jìn)步，并且在部分方法上做了一定的改進(jìn)(陳秋杰等，2013)，建立了靜態(tài)和時(shí)變重力場(chǎng)模型(沈云中等，2015；王長青等，2015).但是，在同時(shí)完成衛(wèi)星精密定軌與重力場(chǎng)解算的同解法上，目前國內(nèi)還沒有完全突破.而國際上，GRACE數(shù)據(jù)產(chǎn)品的官方機(jī)構(gòu)CSR與GFZ，采用的都是同解法提供的產(chǎn)品.可以說，同解法是衛(wèi)星重力場(chǎng)模型研制中的一項(xiàng)關(guān)鍵技術(shù)，國內(nèi)目前尚未完全掌握，當(dāng)前的工作主要是在仿真研究的基礎(chǔ)上(鄒賢才等, 2010a)，對(duì)關(guān)鍵問題分別進(jìn)行技術(shù)攻關(guān).作者通過多年工作，從仿真研究開始，到動(dòng)力法中的關(guān)鍵技術(shù)，比如并行計(jì)算、星載數(shù)據(jù)的精細(xì)處理等等，解決了一系列技術(shù)問題，不斷完善方法的同時(shí)，從無到有建立了軟件平臺(tái)，完成了同解法平臺(tái)的聯(lián)調(diào)聯(lián)試.

本文主要從觀測(cè)模型的特征分析、數(shù)據(jù)處理的主要方案等方面展開，通過實(shí)測(cè)數(shù)據(jù)處理剖析同解法的特點(diǎn).目的是希望能推動(dòng)同解法在海量實(shí)測(cè)衛(wèi)星觀測(cè)數(shù)據(jù)處理中的應(yīng)用，提高我國在衛(wèi)星重力場(chǎng)模型研制領(lǐng)域的技術(shù)水平，進(jìn)一步縮小與國際一流水平的差距，并促進(jìn)我國對(duì)國際上新一代重力衛(wèi)星數(shù)據(jù)資源的有效利用.上述研究對(duì)我國發(fā)展自主衛(wèi)星重力測(cè)量系統(tǒng)，以及衛(wèi)星重力模型的研制理論和實(shí)際應(yīng)用都有積極的促進(jìn)意義.

2理論模型與方法

在分別論述各關(guān)鍵技術(shù)之前，先對(duì)本文的術(shù)語做簡單約定，以免混淆.后文所述的動(dòng)力法，相對(duì)同解法(SSM，Simultaneous Solution Method)而言，可稱“兩步法”.它以衛(wèi)星精密軌道為觀測(cè)值，利用動(dòng)力學(xué)方法確定重力場(chǎng)模型和其他有關(guān)參數(shù)(周旭華，2005；周旭華等，2006；肖云,2006;鄒賢才等, 2010a).幾何法，則是利用星載GPS(Global Positioning System)觀測(cè)數(shù)據(jù)，按照動(dòng)態(tài)定位的原理所做的精密定軌，也稱為運(yùn)動(dòng)學(xué)定軌.另外，本文依然采用動(dòng)力法中的分弧段思想，即將觀測(cè)數(shù)據(jù)分解為多個(gè)弧段，每個(gè)弧段將軌道初值以及相關(guān)模型參數(shù)設(shè)置為未知參數(shù).通過多弧段法方程綜合獲得最終的參數(shù)解算結(jié)果，其基本數(shù)學(xué)模型參見鄒賢才等(2010b).

利用同解法處理星星跟蹤數(shù)據(jù)的實(shí)質(zhì)是如何利用星載GPS觀測(cè)數(shù)據(jù)和低低星間觀測(cè)值，同時(shí)確定衛(wèi)星的軌道初值、地球重力場(chǎng)模型位系數(shù)、加速度計(jì)校準(zhǔn)參數(shù)、模糊度、衛(wèi)星接收機(jī)鐘差、以及其他相關(guān)參數(shù).在分別組成了相位觀測(cè)數(shù)據(jù)以及星間距離變率的法方程后，可通過方差分量估計(jì)合理定權(quán)(崔希璋等，2009)，獲得最終的參數(shù)解算結(jié)果.下面分別闡述星載 GPS 載波相位以及星間距離變率的觀測(cè)方程.

2.1星載GPS載波相位觀測(cè)

基本的載波相位觀測(cè)方程可簡寫為(Xu,2003)

λ φ(tr)=ρ(ts,tr)-c δts+c δtr+δiono+δother-N λ+εL(tr),

(1)式中，tr和ts分別為信號(hào)的接收和發(fā)射時(shí)刻，δtr和δts分別為接收機(jī)鐘和星鐘的鐘差.φ(tr)為載波相位觀測(cè)值，λ為波長，N為對(duì)應(yīng)的整周模糊度.ρ(ts,tr)為發(fā)射天線相位中心到接收機(jī)天線相位中心間的幾何距離；c為真空中光速.δiono為電離層延遲改正；δother為其他各項(xiàng)改正，包括地球自轉(zhuǎn)、相位解纏(Wu et al.,1993)等，相關(guān)改正公式可參考文獻(xiàn).εL(tr)為載波相位觀測(cè)的測(cè)量噪聲.注意，本文處理的是星載GPS數(shù)據(jù)，因此無需考慮對(duì)流層改正.自然地，在實(shí)際的數(shù)據(jù)處理分析中，依然采用消電離層組合觀測(cè)值.

根據(jù)星載GPS載波相位的觀測(cè)方程及其非差組合觀測(cè)值，以及選定的同解法參數(shù)，可以獲得星載GPS載波相位觀測(cè)值對(duì)應(yīng)的設(shè)計(jì)矩陣.同時(shí)注意到，求偏導(dǎo)的基本原則是觀測(cè)泛函對(duì)所有選定的參數(shù)P進(jìn)行，因此設(shè)計(jì)矩陣有下列一般形式：

(2)

式中，F(xiàn)為觀測(cè)泛函，與具體的觀測(cè)值有關(guān).但選定某個(gè)類型的觀測(cè)值后，它的表達(dá)式一定是明確的，因此本文以F統(tǒng)一表達(dá).ΔT、X0、N、Acc、EGM分別為重力衛(wèi)星的鐘差、軌道初值(包括位置和速度)、模糊度、加速度計(jì)校準(zhǔn)參數(shù)以及重力場(chǎng)模型位系數(shù)參數(shù)組成的列向量.為了建立任意歷元觀測(cè)值與軌道初值的關(guān)系，借助了觀測(cè)歷元的軌道位置和速度XLEO.下標(biāo)1，2，…,n表示歷元編號(hào).

對(duì)于上述載波觀測(cè)數(shù)據(jù)的設(shè)計(jì)矩陣，它是由兩部分組成的，一部分是ΔT和N.這兩類參數(shù)加上衛(wèi)星的位置參數(shù)，構(gòu)成了幾何法精密定軌的參數(shù)(李建成等, 2009).另一部分是X0、Acc和EGM.它們屬于動(dòng)力法的參數(shù).因此，實(shí)現(xiàn)同解法的一個(gè)可能方案就是，通過幾何法精密定軌與動(dòng)力學(xué)方法兩項(xiàng)技術(shù)，分別攻關(guān)后集成.這也正是本文研究采取的技術(shù)路線.

這其中有幾個(gè)關(guān)鍵特征需要強(qiáng)調(diào).在同解法的待估參數(shù)中，接收機(jī)鐘差(后文簡稱為鐘差)為每個(gè)歷元設(shè)置一個(gè).因此，在一個(gè)弧段中，鐘差參數(shù)的個(gè)數(shù)等于有效觀測(cè)歷元數(shù).對(duì)某個(gè)歷元而言，對(duì)應(yīng)某個(gè)載波觀測(cè)值有且只有一個(gè)鐘差參數(shù)，因此式(2)中標(biāo)示為“1”的偏導(dǎo)數(shù)，除對(duì)應(yīng)歷元的值為常數(shù)外，其余值都為零.這種特殊結(jié)構(gòu)導(dǎo)致相位觀測(cè)值對(duì)應(yīng)的法方程系數(shù)陣具有一種稀疏結(jié)構(gòu)，參見圖 1.其中，①的對(duì)角部分表示與鐘差有關(guān)的對(duì)應(yīng)法方程，它是一個(gè)對(duì)角陣.按照本文的設(shè)置，弧段長度最大值設(shè)置為24h，根據(jù)GRACE星載載波數(shù)據(jù)采樣間隔10s計(jì)算，該對(duì)角陣的維數(shù)可以達(dá)到8640.②是與模糊度參數(shù)有關(guān)的部分.對(duì)模糊度參數(shù)而言，它具有和鐘差參數(shù)類似的特點(diǎn).對(duì)某個(gè)特定歷元的相位觀測(cè)值，一定依附于某顆衛(wèi)星，自然除了這顆衛(wèi)星的模糊度參數(shù)外，觀測(cè)值對(duì)其他模糊度參數(shù)的偏導(dǎo)數(shù)一定為零.因此也形成了對(duì)角陣.③是與軌道初值、重力場(chǎng)模型位系數(shù)、加速度計(jì)校準(zhǔn)參數(shù)等力學(xué)參數(shù)有關(guān)的部分，一般情況下是滿陣.對(duì)于②，一天的弧長包含的模糊度參數(shù)總數(shù)一般不超過500.因此，本文設(shè)計(jì)算法時(shí)利用①的稀疏性，而將②和③統(tǒng)一當(dāng)作一般的矩陣進(jìn)行處理，這樣在消參數(shù)與回代求解(鄒賢才,2007;鄒賢才等, 2010b)時(shí)，可以兼顧運(yùn)行效率與實(shí)現(xiàn)的難易度.上述幾類參數(shù)的非對(duì)角部分當(dāng)作一般的矩陣處理.法方程的系數(shù)陣是對(duì)稱陣，圖1展示了其下三角部分.

圖1　同解法GPS 觀測(cè)數(shù)據(jù)弧段法方程的結(jié)構(gòu)示意圖Fig.1　Sketch on the structure of the normal equation for the GPS data in SSM

2.2星間距離變率觀測(cè)值的觀測(cè)模型

(3)

和

(4)

式(3)和式(4)即所謂的觀測(cè)模型.根據(jù)觀測(cè)模型，可以獲得星間距離變率對(duì)衛(wèi)星狀態(tài)的偏導(dǎo)數(shù).由式(4)可得星間距離變率對(duì)兩顆衛(wèi)星狀態(tài)的偏導(dǎo)數(shù)為(周旭華,2005;鄒賢才, 2007)：

(5)

(6)

因此，可寫出星間距離變率對(duì)應(yīng)的設(shè)計(jì)矩陣Bkrr為：

(7)

式中雙星狀態(tài)對(duì)力學(xué)參數(shù)的偏導(dǎo)數(shù)，與公式(2)中的對(duì)應(yīng)部分計(jì)算方法相同.因此，星間距離變率觀測(cè)數(shù)據(jù)的法方程可以在動(dòng)力法工作的基礎(chǔ)上直接完成.由此，聯(lián)合解算這兩類觀測(cè)數(shù)據(jù)的法方程，可以獲得同解法中的各項(xiàng)參數(shù).

2.3模型設(shè)置以及參數(shù)配置

本節(jié)簡要介紹本文采用的基本參數(shù)、基本的力模型配置和參數(shù)估計(jì)方案.參見表1.

對(duì)于星載加速度計(jì)數(shù)據(jù)的處理，研究團(tuán)隊(duì)已經(jīng)在相關(guān)文章中探討過該問題，因而表1僅列舉了參數(shù)配置方案.有興趣的讀者可以參考鄒賢才等(2015)，后文將做進(jìn)一步的數(shù)值檢驗(yàn).

最后簡要討論一下并行計(jì)算.已有研究表明，衛(wèi)星重力數(shù)據(jù)處理中，分弧段處理以及法方程綜合技術(shù)非常適合在共享內(nèi)存的平臺(tái)上實(shí)現(xiàn)，并且無須設(shè)計(jì)并行算法，采用批處理即可實(shí)現(xiàn).在法方程綜合階段，計(jì)算熱點(diǎn)之一是大型線性方程組求解.注意到，同解法的參數(shù)解算中，法方程均為對(duì)稱陣，且正常條件下都是正定陣，自然可以直接采用已有成果(鄒賢才等,2010b).本文的所有計(jì)算工作都在武漢大學(xué)測(cè)繪學(xué)院的計(jì)算集群上實(shí)施(含7個(gè)Intel節(jié)點(diǎn)，共192核)，搭載CentOS LINUX操作系統(tǒng)，任務(wù)管理通過SLURM實(shí)現(xiàn).

表1　同解法采用的模型以及參數(shù)配置方案

3數(shù)值分析與討論

前已述及，同解法的最大特點(diǎn)是能以更統(tǒng)一的模型對(duì)參數(shù)進(jìn)行整體求解，理論模型更加嚴(yán)密.因而同解法的解算結(jié)果中包含了多種信息，包括精化的重力場(chǎng)模型、加速度計(jì)校準(zhǔn)參數(shù)、精化的軌道初值以及與GPS觀測(cè)有關(guān)的接收機(jī)鐘差、模糊度等.其中有些參數(shù)的解算結(jié)果目前沒有條件進(jìn)行檢核，比如模糊度和鐘差的解算結(jié)果.有些參數(shù)則可以設(shè)計(jì)合理方案進(jìn)行檢核.這正是本節(jié)的主要工作，通過科學(xué)合理的檢核方案，對(duì)同解法相關(guān)產(chǎn)品進(jìn)行評(píng)價(jià).

3.1軌道質(zhì)量評(píng)估

與星載GPS幾何法定軌不同(李建成等, 2009)，同解法并不是直接給出衛(wèi)星的精密軌道.它先是做力模型精化，然后通過軌道積分給出軌道.這一點(diǎn)與簡化的動(dòng)力學(xué)方法(Yunck et al.,1990)類似.同時(shí)考慮到JPL(Jet Prolusion Laboratory)的簡化動(dòng)力學(xué)軌道(GRACE官方發(fā)布的精密軌道)質(zhì)量很高，為國內(nèi)外許多研究群體采用，并且可以提取到時(shí)變信號(hào)(冉將軍等, 2014)，而且它與同解法一樣，能提供動(dòng)力學(xué)意義上的衛(wèi)星速度，所以本文將同解法給出的軌道(一個(gè)月的數(shù)據(jù))與JPL的軌道進(jìn)行了比較.給出軌道徑向(R: Radial)、切向(T: Tangent)和法向(N: Normal)差異的統(tǒng)計(jì)結(jié)果.參見表 2.從該表可以看出，兩套軌道在 RTN 三方向的均方根不超過3 cm.其中，軌道切向的精度在2.7 cm左右，另外兩個(gè)方向優(yōu)于2 cm.上述結(jié)論對(duì)兩顆衛(wèi)星均成立.相比已有的幾何法定軌結(jié)果(李建成等, 2009)，同解法的軌道精度更高，不僅如此，在顧及兩套軌道誤差的情況下，軌道差異的極值也只有11 cm.進(jìn)一步，兩套軌道的平均偏差不超過1 cm，軌道的一致性好.這些結(jié)果從側(cè)面表明本文在力模型配置上的合理性，說明動(dòng)力學(xué)模型在定軌中起到了良好的約束作用.

表2　同解法軌道與JPL 精密軌道的比較結(jié)果(單位： cm)

如果選擇一個(gè)單獨(dú)的弧段(弧長23小時(shí)59分)比較，不僅可以很清晰地看出與 JPL 軌道的差異，還可以發(fā)現(xiàn)差值的分布特點(diǎn).如圖 2所示.從絕對(duì)值上講，三個(gè)方向(RTN分別用紅、綠、藍(lán)表示)都在幾厘米左右，可以認(rèn)為兩者的軌道精度相當(dāng)，但是軌道差異中存在明顯的周期性現(xiàn)象，且與GRACE衛(wèi)星的軌道周期一致，這表明軌道殘差依然受到了類似重力場(chǎng)因素的影響.合理的解釋是簡化動(dòng)力學(xué)方法中采用的重力場(chǎng)沒有精化所致——JPL的精密軌道是采用簡化動(dòng)力學(xué)方法建立的，不是純粹意義上的動(dòng)力學(xué)定軌，沒有精化重力場(chǎng)模型這一過程.因此，在實(shí)現(xiàn)“最優(yōu)”軌道的同時(shí)，它在力模型上的缺點(diǎn)通過軌道積分，特別是較長時(shí)段的軌道積分放大，比較明顯地顯示出來了.

圖2　動(dòng)力學(xué)軌道與JPL精密軌道差異在RTN中的分布Fig.2　Difference between the SSM and JPL orbit products and its distribution in RTN

為了進(jìn)一步驗(yàn)證該思路，本文采用高精度的星間距離變率觀測(cè)值對(duì)軌道進(jìn)行檢核.具體過程如下.簡化動(dòng)力學(xué)軌道與同解法的動(dòng)力學(xué)軌道都能給出衛(wèi)星的位置和速度，因此可以利用軌道導(dǎo)出星間距離變率觀測(cè)值.通過導(dǎo)出量與實(shí)測(cè)值比較，顧及實(shí)測(cè)值達(dá)到10-6m·s-1左右的精度，可以對(duì)軌道和星間距離變率之間的一致性進(jìn)行檢驗(yàn).結(jié)果參見圖 3.該圖中，振幅較大(紅色)曲線是由JPL精密軌道導(dǎo)出的星間距離變率與觀測(cè)值之間的差異(JPLD)，振幅較小(綠色)曲線是由同解法精密軌道導(dǎo)出的星間距離變率與觀測(cè)值之間的差異(SSMD).很明顯，同解法給出的精密軌道與星間距離變率觀測(cè)值之間的一致性要好得多.數(shù)值上，JPLD分布在-7.3×10-6/6.5×10-6之間，RMS(root mean squares)為2.5×10-6，SSMD分布在-1.3×10-6/1.4×10-6之間，RMS為2.7×10-7.相較而言，SSMD不僅更加穩(wěn)定，而且與星間距離變率更吻合，RMS要高一個(gè)量級(jí).對(duì)于其他弧段，有類似的結(jié)果.因此，本文認(rèn)為同解法給出的軌道、速度與精密星間距離變率觀測(cè)值之間的一致性更好，這也從側(cè)面驗(yàn)證了同解法在精密定軌的同時(shí)精化重力場(chǎng)模型的重要優(yōu)勢(shì).

圖3　JPL與同解法精密軌道導(dǎo)出的星間距離變率與實(shí)測(cè)值之間的差異Fig.3　Consistency between the derived range rate from the JPL and SSM orbit products

上述軌道質(zhì)量檢核是利用星間距離變率對(duì)同解法和簡化動(dòng)力學(xué)軌道位置和速度的綜合檢驗(yàn).衛(wèi)星激光測(cè)距(SLR，Satellite Laser Ranging)是一種軌道質(zhì)量檢核的保留手段，因此本文也實(shí)施了SLR檢驗(yàn).做地面跟蹤站速度場(chǎng)、潮汐、偏心改正，衛(wèi)星質(zhì)心、相對(duì)論、大氣折射以及時(shí)距偏差等改正(Hulley and Pavlis, 2007)，得到星地距離檢核衛(wèi)星軌道，情況如表 3和表 4.以觀測(cè)數(shù)據(jù)最為豐富的跟蹤站7090(澳大利亞堪培拉站)為例，雙星的徑向精度在2 cm以內(nèi)，進(jìn)一步驗(yàn)證了表 2中結(jié)果的可靠性與軌道質(zhì)量.

表3　SLR 觀測(cè)檢核 GRACE-A衛(wèi)星的軌道

表4　SLR 觀測(cè)檢核GRACE-B衛(wèi)星的軌道

3.2加速度計(jì)校準(zhǔn)結(jié)果檢核

加速度計(jì)校準(zhǔn)是衛(wèi)星重力數(shù)據(jù)處理中的一項(xiàng)基本工作，但目前沒有準(zhǔn)確可靠的校準(zhǔn)參數(shù)作為外部檢核標(biāo)準(zhǔn).即使GRACE團(tuán)隊(duì)發(fā)布了技術(shù)文檔(Bettadpur,2009)，也只能作參考.根據(jù)GRACE官方手冊(cè)，衛(wèi)星星固坐標(biāo)系x軸實(shí)際上是兩顆衛(wèi)星之間的連線.因此這兩顆衛(wèi)星各自的星固系x軸重合，但是方向相反.沿軌的非保守力源主要是大氣阻力.這兩顆衛(wèi)星的空間距離在200 km左右，因此大氣阻力之間的相關(guān)性很強(qiáng).在準(zhǔn)備本文工作時(shí)，國外也正好有相關(guān)的成果發(fā)表(Kim and Tapley, 2015).為此，下面展示這種相關(guān)性，從側(cè)面驗(yàn)證同解法的非保守力模型精化效果.

利用同解法給出的加速度計(jì)校準(zhǔn)參數(shù)對(duì)原始的加速度計(jì)數(shù)據(jù)進(jìn)行校準(zhǔn)，并表示在衛(wèi)星的星固系中，這樣就得到了x方向的非保守力.兩顆衛(wèi)星的空間相關(guān)性使得它們各自受到的非保守力也存在很強(qiáng)的相關(guān)性，同時(shí)注意到兩顆衛(wèi)星星固系x軸方向相反，因此可以預(yù)見它們的x方向的非保守力大小接近，符號(hào)相反.這一特性，為實(shí)測(cè)數(shù)據(jù)的處理結(jié)果所證實(shí).參見圖 4.其中，紅、綠分別表示兩顆衛(wèi)星的非保守力，藍(lán)色為其差別.

圖4　GRACE雙星星固系沿跡方向的非保守力及其差異Fig.4　Non-gravitational forces along track for the GRACE twins and their differences

從該圖可以看出，兩顆衛(wèi)星沿軌方向的非保守力振幅接近，但是相位相反.通過求和之后的分布(藍(lán)色曲線)分析，可以更加清楚地看出這個(gè)特點(diǎn).結(jié)合它的基礎(chǔ)工作(鄒賢才等, 2015)，可以認(rèn)為同解法的加速度計(jì)校準(zhǔn)內(nèi)符合也非常好.這里也同時(shí)給出GRACE 雙星加速度計(jì)的比例系數(shù)供參考，其估計(jì)策略請(qǐng)參見表 1.GRACE-A 衛(wèi)星的三軸比例系數(shù)及其精度分別為(0.9701,0.9255,0.9374)、(7×10-5, 3×10-3, 4×10-4)，GRACE-B衛(wèi)星的比例系數(shù)為(0.9549,0.9418,0.9309)，精度指標(biāo)相同.

3.3動(dòng)力學(xué)模型配置檢驗(yàn)

與3.2節(jié)不同，本節(jié)在更宏觀的角度研究力模型配置的合理性，并主要通過重力場(chǎng)模型精度檢核展開.目前國外發(fā)布的月重力場(chǎng)模型并不是都有精度信息，很難簡單的直接應(yīng)用誤差階方差以及GPS水準(zhǔn)的檢核方法，因此也做了同解法重力場(chǎng)模型(SSM_EGM)與兩步法結(jié)果的對(duì)比.

基本考慮如下.同解法的觀測(cè)模型涵蓋了傳統(tǒng)的動(dòng)力學(xué)方法(“兩步法”)的所有參數(shù)，作者設(shè)計(jì)的軟件也支持運(yùn)行模式上的切換.如果將它的輸入數(shù)據(jù)修改為軌道，就可以用兩步法反演重力場(chǎng)模型及相關(guān)參數(shù).首先，在分析中測(cè)試?yán)猛夥ńo出的純動(dòng)力法軌道，按照兩步法反演重力場(chǎng).結(jié)果是反演的重力場(chǎng)模型與同解法確定的重力場(chǎng)模型一致，其精度指標(biāo)接近計(jì)算機(jī)的截?cái)嗾`差(在10-14～10-15之間).對(duì)于同一套軟件，力模型配置可以完全保持不變，參見表 1(注意，相關(guān)動(dòng)力學(xué)參數(shù)仍設(shè)置為待解參數(shù)進(jìn)行估計(jì)，并不是簡單的固定為同解法的結(jié)果).因此它只能看作是軟件測(cè)試中的正反算.但是這個(gè)現(xiàn)象說明，如果基于JPL精密軌道的兩步法反演結(jié)果不好，則說明雙方在力模型配置上的差異比較明顯.

檢驗(yàn)思路如下.由于無法確切知道JPL在定軌中采用的力學(xué)模型配置，特別是簡化動(dòng)力學(xué)定軌(Yunck et al., 1990)中有可能視需要加入一些所謂的脈沖加速度信號(hào)，所以力模型上的差異使得本文無法完全恢復(fù)它在定軌中采用的重力場(chǎng)模型.但是，根據(jù)簡化動(dòng)力學(xué)方法定軌的基本原理，它需要選定某個(gè)重力場(chǎng)模型(顯然它也是主要攝動(dòng)力)，輔以其他保守力模型(包括各種潮汐以及非潮汐時(shí)變模型)和非保守力模型(可由加速度計(jì)提供，或者結(jié)合經(jīng)驗(yàn)?zāi)Ｐ偷鹊?，完成精密定軌.因此，通過簡化動(dòng)力學(xué)軌道反演重力場(chǎng)(記為 JPL_EGM)時(shí)，一方面，不會(huì)獲得前述正反算中的內(nèi)符合精度.另一方面，如果在力模型的配置上，本文的配置方案與JPL的方案不出現(xiàn)大的偏差，按照前述同解法軌道反演重力場(chǎng)的分析，也應(yīng)該能很好地反演重力場(chǎng)模型，并能通過 GPS 水準(zhǔn)進(jìn)行獨(dú)立的外部檢核.本文依然采用美國的6169個(gè)GPS水準(zhǔn)點(diǎn).結(jié)果表明，SSM_EGM的精度為0.757 m，而JPL_EGM的精度為0.762 m(參與比較的截?cái)嚯A次為96階，參見表1). CSR最新的GGM05S模型對(duì)應(yīng)指標(biāo)為0.755 m，顧及到 GGM05S 使用了多年的數(shù)據(jù)，上述精度表明本文解算的重力場(chǎng)模型質(zhì)量也是可靠的，與國際機(jī)構(gòu)其他發(fā)布模型的分階次比較結(jié)果詳見表5，總體上看，隨著重力場(chǎng)模型階次的提高，模型大地水準(zhǔn)面的精度不斷提高；對(duì)不同的模型而言，在本文選擇的截?cái)嚯A上，精度相當(dāng).

最后，本文給出同解法重力場(chǎng)模型與國際上幾種主要的重力場(chǎng)模型之間差異的對(duì)比，以及與 CSR 月重力場(chǎng)模型(同階次模型UTCSR)之間的對(duì)比結(jié)果，參見圖5.基本的做法是假定EIGEN6C4的精度“最高”，然后各個(gè)模型分別與 EIGEN6C4作差，把差值看作位系數(shù)的誤差，計(jì)算所有階次上大地水準(zhǔn)面的誤差，從而表示重力場(chǎng)模型的精度.從作差后大地水準(zhǔn)面每階的差異來看，EGM2008與EIGEN6C4之間吻合較好，GGM05S 和 EIGEN6C4符合得最好，并且靜態(tài)模型之間的符合程度一致優(yōu)于月解模型SSMD和UTCSR.對(duì)于兩個(gè)月解模型，SSMD在13階以內(nèi)，特別是J2項(xiàng)上表現(xiàn)要好.在60階以內(nèi)，兩者與靜態(tài)模型的符合程度，各有優(yōu)勢(shì)，但是在20～35階左右，UTCSR的差異更加穩(wěn)定.60階之后，UTCSR 模型和 SSMD的表現(xiàn)比較接近.總體上講，結(jié)合GPS水準(zhǔn)檢驗(yàn)與大地水準(zhǔn)面誤差分析，再綜合前述的積分軌道的質(zhì)量檢驗(yàn)(參見圖3、表3和表4)，本文認(rèn)為SSMD與國際上的同類型月解模型是相當(dāng)?shù)?

表5　利用GPS水準(zhǔn)點(diǎn)的重力場(chǎng)模型精度檢核(單位：m)

圖5　同解法重力場(chǎng)模型(SSM_EGM)與其他重力場(chǎng)模型的差異對(duì)比Fig.5　Comparison on the differences between the SSM_EGM and other models

4結(jié)論與展望

本文以方法研究為主，提出了同解法的一種實(shí)現(xiàn)方案，并結(jié)合GRACE一個(gè)月的實(shí)測(cè)數(shù)據(jù)，在高性能計(jì)算平臺(tái)上實(shí)現(xiàn)了高質(zhì)量衛(wèi)星精密軌道與衛(wèi)星重力場(chǎng)模型的同時(shí)解算.對(duì)動(dòng)力法精密軌道、重力場(chǎng)模型以及相關(guān)附屬產(chǎn)品做了較為全面的評(píng)價(jià)與分析.結(jié)果表明，衛(wèi)星精密軌道與JPL軌道產(chǎn)品互檢，RTN三方向的精度均優(yōu)于3 cm；與SLR比較，徑向精度優(yōu)于2 cm.并且，本文通過星間距離變率檢核了同解法的純動(dòng)力學(xué)軌道與JPL簡化動(dòng)力學(xué)軌道，結(jié)果表明同解法的動(dòng)力學(xué)軌道與高精度距離變率之間符合更好.同時(shí)，本文還通過GRACE任務(wù)雙星之間非保守力的空間相關(guān)性，以及GNSS水準(zhǔn)對(duì)本文的力模型配置進(jìn)行了檢驗(yàn)，進(jìn)一步證明了力模型設(shè)置以及相關(guān)產(chǎn)品的可靠性.

綜上所述，本文的研究已經(jīng)涉及了同解法流程中的各個(gè)環(huán)節(jié).產(chǎn)品質(zhì)量分析表明，本文已經(jīng)突破了同解法的各項(xiàng)關(guān)鍵技術(shù)，在國內(nèi)率先實(shí)現(xiàn)了從仿真研究到實(shí)測(cè)飛行任務(wù)數(shù)據(jù)分析的重要進(jìn)展.

在同解法領(lǐng)域，未來仍然有許多問題可以做深入研究.比如，可以引入模糊度固定技術(shù)，探討如何進(jìn)一步提高軌道質(zhì)量以及重力場(chǎng)模型的精度，甚至從方法上對(duì)星星跟蹤觀測(cè)模型嚴(yán)密化，在更高的層次上建立統(tǒng)一的觀測(cè)模型，深化同解法及其產(chǎn)品在科學(xué)研究中的作用.上述工作涉及到理論、方法研究以及計(jì)算方案的大調(diào)整，值得做進(jìn)一步的深入研究.

致謝感謝JPL提供的GRACE任務(wù)觀測(cè)數(shù)據(jù)以及武漢大學(xué)測(cè)繪學(xué)院計(jì)算中心的大力支持.

References

Bettadpur S. 2009. Recommendation for A-Priori Bias & Scale Parameters for Level-1B ACC Data (Version 2). GRACE TN-02.

Chen Q J, Shen Y Z, Zhang X F. 2013. Linearization method of recovering earth′s gravity field with respect to gravity satellite′s kinematic orbits. Chinese Journal of Geophysics (in Chinese), 56(7): 2238-2244, doi: 10.6038/cjg20130711.

Cui X Z, Yu Z C, Tao B Z, et al. 2009. Generalized Surveying Adjustment. 2nd ed. (in Chinese). Wuhan: Wuhan University Press.

Hulley G C, Pavlis E C. 2007. A ray-tracing technique for improving satellite laser ranging atmospheric delay corrections, including the effects of horizontal refractivity gradients. Journal of Geophysical Research: Solid Earth, 112(B6): B06417. Jekeli C. 1999. The determination of gravitational potential differences from satellite-to-satellite tracking. Celestial Mechanics and Dynamical Astronomy, 75(2): 85-101. Kim J, Tapley B D. 2015. Estimation of non-gravitational acceleration difference between two co-orbiting satellites using single accelerometer data. Journal of Geodesy, 89(6): 537-550. Li J C, Zhang S J, Zou X C, et al. 2010. Precise orbit determination for GRACE with zero-difference kinematic method. Chinese Science Bulletin, 54(16): 2355-2362, doi: 10.1007/s11434-009-0286-0.Mayer-Gürr T, Ilk K H, Eicker A, et al. 2005. ITG-CHAMP01: a CHAMP gravity field model from short kinematic arcs over a one-year observation period. Journal of Geodesy, 78(7): 462-480.

Peng D J, Wu B. 2007. Discussion about precise orbit determination of LEO satellite using single-difference and indifference method. Chinese Science Bulletin (in Chinese), 52(6): 715-719.

Qin X P. 2009. Research on precision orbit determination theory and method of low earth orbiter based on GPS technique [Ph. D. thesis] (in Chinese). Zhengzhou: The PLA Information Engineering University. Ran J J, Xu H Z, Zhong M, et al. 2014. Global temporal gravity filed recovery using GRACE data. Chinese Journal of Geophysics (in Chinese), 57(4): 1032-1040, doi: 10.6038/cjg20140402. Reigber C. 1989. Gravity field recovery from satellite tracking data.∥ Sansò F, Rummel R eds. Theory of Satellite Geodesy and Gravity Field Determination. Berlin: Springer, 197-234.

Reubelt T, Austen G, Grafarend E W. 2003. Harmonic analysis of the earth′s gravitational field by means of semi-continuous ephemerides of a low earth orbiting GPS-tracked satellite. Case study: CHAMP. Journal of Geodesy, 77(5-6): 257-278.

Wang C Q, Xu H Z, Zhong M, et al. 2015. An investigation on GRACE temporal gravity field recovery using the dynamic approach. Chinese Journal of Geophysics (in Chinese), 58(3): 756-766, doi: 10.6038/cjg20150306.

Wolff M. 1969. Direct measurements of the earth′s gravitational potential using a satellite pair. Journal of Geophysical Research, 74(22): 5295-5300. Wu J T, Wu S C, Hajj G A, et al. 1993. Effects of antenna orientation on GPS carrier phase. Manuscripta Geodaetica, 18(2): 91-98.

Xiao Y. 2006. Analysis of earth gravity field recovery by satellite-to-satellite tracking data [Ph. D. thesis] (in Chinese). Zhengzhou: The PLA Information Engineering University. Xu G C. 2003. GPS: Theory, Algorithms and Applications. Berlin: Springer-Verlag.

Yunck T P, Wu S C, Wu J T, et al. 1990. Precise tracking of remote sensing satellites with the global positioning system. IEEE Transactions on Geoscience and Remote Sensing, 28(1): 108-116.

Zhang X F. 2007. The earth's field model recovery on the basis of satellite-to-satellite tracking mission [Ph. D. thesis] (in Chinese). Shanghai: Tongji University. Zhao Q L, Shi C, Liu X L, et al. 2008. Determination of precise orbit using onboard GPS data for gravity modeling oriented satellites. Geomatics and Information Science of Wuhan University (in Chinese), 33(8): 810-814. Zheng Z Y. 2004. Study and software implementation of GPS data pre-processing and onboard GPS kinematic orbit determination [Ph. D. thesis] (in Chinese). Shanghai: Shanghai Astronomical Observatory, Chinese Academy of Sciences. Zhou X H. 2005. Study on satellite gravity and its application [Ph. D. thesis] (in Chinese). Wuhan: Institute of Geodesy and Geophysics, Chinese Academy of Sciences.

Zhou X H, Hsu H, Wu B, et a1. 2006. Earth′s gravity field derived from GRACE satellite tracking data. Chinese Journal of Geophysics (in Chinese), 49(3): 718-723.

Zhu S, Reigber C, K?nig R. 2004. Integrated adjustment of CHAMP, GRACE, and GPS data. Journal of Geodesy, 78(1-2): 103-108.

Zou X C. 2007. Theory of satellite orbit and earth gravity field determination [Ph. D. thesis] (in Chinese). Wuhan: Wuhan University.

Zou X C, Li J C, Jiang W P, et al. 2010a. Research on the simultaneous solution method for satellite gravity data analysis and its simulation. Acta Geodaetica et Cartographica Sinica (in Chinese), 39(4): 344-348.

Zou X C, Li J C, Wang H H, et al. 2010b. Application of parallel computing with openMP in data processing for satellite gravity. Acta Geodaetica et Cartographica Sinica (in Chinese), 39(6): 636-641.

Zou X C, Li J C, Zhong L P, et al. 2015. Calibration of the accelerometers onboard GRACE with the dynamic method. Geomatics and Information Science of Wuhan University (in Chinese), 40(3): 357-360.

附中文參考文獻(xiàn)

陳秋杰, 沈云中, 張興福. 2013. 基于重力衛(wèi)星幾何軌道線性化的地球重力場(chǎng)反演方法. 地球物理學(xué)報(bào), 56(7): 2238-2244, doi: 10.6038/cjg20130711.

崔希璋, 於宗儔, 陶本藻等. 2009. 廣義測(cè)量平差. 2版. 武漢: 武漢大學(xué)出版社.

李建成, 張守建, 鄒賢才等. 2009. GRACE衛(wèi)星非差運(yùn)動(dòng)學(xué)厘米級(jí)定軌. 科學(xué)通報(bào), 54(16): 2355-2362.

彭冬菊, 吳斌. 2007. 非差和單差LEO星載GPS精密定軌探討. 科學(xué)通報(bào), 52(6): 715-719.

秦顯平. 2009. 星載GPS低軌衛(wèi)星定軌理論及方法研究[博士論文]. 鄭州: 中國人民解放軍信息工程大學(xué).

冉將軍, 許厚澤, 鐘敏等. 2014. 利用GRACE重力衛(wèi)星觀測(cè)數(shù)據(jù)反演全球時(shí)變地球重力場(chǎng)模型. 地球物理學(xué)報(bào), 57(4): 1032-1040, doi: 10.6038/cjg20140402.

沈云中, 陳秋杰, 許厚澤等. 2015. 基于GRACE觀測(cè)軌道線性化的衛(wèi)星重力反演模型與結(jié)果. 大地測(cè)量與導(dǎo)航專業(yè)委員會(huì)2015學(xué)術(shù)年會(huì)暨研究生學(xué)術(shù)論壇.

王長青, 許厚澤, 鐘敏等. 2015. 利用動(dòng)力學(xué)方法解算GRACE時(shí)變重力場(chǎng)研究. 地球物理學(xué)報(bào), 58(3): 756-766, doi: 10.6038/cjg20150306.

肖云. 2006. 基于衛(wèi)星跟蹤衛(wèi)星數(shù)據(jù)恢復(fù)地球重力場(chǎng)的研究[博士論文]. 鄭州: 中國人民解放軍信息工程大學(xué).

張興福. 2007. 應(yīng)用低軌衛(wèi)星跟蹤數(shù)據(jù)反演地球重力場(chǎng)模型[博士論文]. 上海: 同濟(jì)大學(xué).

趙齊樂, 施闖, 柳響林等. 2008. 重力衛(wèi)星的星載GPS精密定軌. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 33(8): 810-814.

鄭作亞. 2004. GPS數(shù)據(jù)預(yù)處理和星載GPS運(yùn)動(dòng)學(xué)定軌研究及其軟件實(shí)現(xiàn)[博士論文]. 上海: 中國科學(xué)院研究生院(上海天文臺(tái)).

周旭華. 2005. 衛(wèi)星重力及其應(yīng)用研究[博士論文]. 武漢: 中國科學(xué)院測(cè)量與地球物理研究所.

周旭華, 許厚澤, 吳斌等. 2006. 用GRACE衛(wèi)星跟蹤數(shù)據(jù)反演地球重力場(chǎng). 地球物理學(xué)報(bào), 49(3): 718-723.

鄒賢才. 2007. 衛(wèi)星軌道理論與地球重力場(chǎng)模型的確定[博士論文]. 武漢: 武漢大學(xué).

鄒賢才, 李建成, 姜衛(wèi)平等. 2010a. 衛(wèi)星重力資料分析的同解法研究及其仿真. 測(cè)繪學(xué)報(bào), 39(4): 344-348.

鄒賢才, 李建成, 汪海洪等. 2010b. OpenMP并行計(jì)算在衛(wèi)星重力數(shù)據(jù)處理中的應(yīng)用. 測(cè)繪學(xué)報(bào), 39(6): 636-641.

鄒賢才, 李建成, 衷路萍等. 2015. 動(dòng)力法校準(zhǔn)GRACE星載加速度計(jì). 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 40(3): 357-360.

(本文編輯何燕)

基金項(xiàng)目國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃(2013CB73301)，國家自然科學(xué)基金(41004007、41274033、41210006)，測(cè)繪地理信息公益性行業(yè)科研專項(xiàng)(201512001)資助.

作者簡介鄒賢才，男，1978年生，2007年畢業(yè)于武漢大學(xué)，博士、副教授，主要從事大地測(cè)量學(xué)方面的研究.E-mail:xczou@whu.edu.cn

doi:10.6038/cjg20160708 中圖分類號(hào)P223

收稿日期2015-08-27，2016-04-02收修定稿

Simultaneous solution for precise satellite orbit and earth gravity model

ZOU Xian-Cai, ZHONG Lu-Ping, LI Jian-Cheng

SchoolofGeodesyandGeomatics,WuhanUniversity,KeyLaboratoryofGeospaceEnvironmentandGeodesy,MinistryofEducation,Wuhan430079,China

AbstractThe modern satellite gravimetry missions have improved the mid-term and long wavelength parts of the earth gravity model significantly, and extended its usage in the scientific research in the relating subjects. The simultaneous solution method, which is one of the major methods used in the data processing for satellite gravimetry, is still under developing in China. Based on the basic model and some key technologies, especially on its relation with the kinematic orbit determination and the general dynamic method, the strategy on how to implement the simultaneous solution method is given in detail. And the real data from GRACE mission is analyzed with this method with the help of some work already done by our group. Some evaluations, including the external checks on the accuracy of the orbit, the analysis on the non-conservative forces and the earth gravity model validation by the GPS/Leveling benchmarks are conducted. The numerical results show that the strategy on the realization of the simultaneous solution method is valid, and reveal its advantages distinctly. It is assumed that the key technologies of the simultaneous solution method have been solved, and an advance from the simulation stage to the real data analyzing is achieved. Finally, some prospects and future work are discussed to exert the potential ability of this method.

KeywordsSatellite gravimetry； Simultaneous solution method； Accelerometer calibration; GRACE; Dynamic method

鄒賢才， 衷路萍， 李建成. 2016. 衛(wèi)星精密定軌與重力場(chǎng)建模的同解法. 地球物理學(xué)報(bào)，59(7):2413-2423,doi:10.6038/cjg20160708.

Zou X C, Zhong L P, Li J C. 2016. Simultaneous solution for precise satellite orbit and earth gravity model. Chinese J. Geophys. (in Chinese),59(7):2413-2423,doi:10.6038/cjg20160708.