劉煥玲，文漢江，徐新禹，朱廣彬

1.武漢大學(xué)測繪學(xué)院，湖北 武漢430079；2.中國測繪科學(xué)研究院地理空間信息工程國家測繪地理信息局重點實驗室，北京100830；3.國家測繪地理信息局衛(wèi)星測繪應(yīng)用中心，北京101300

1 引 言

2009年3月17日，歐空局成功發(fā)射了GOCE（Gravity field and steady－state Ocean Circulation Explorer）衛(wèi)星，旨在探測高精度高分辨率的地球重力場和穩(wěn)態(tài)海洋環(huán)流，具體目標(biāo)為：利用衛(wèi)星引力梯度數(shù)據(jù)和衛(wèi)星跟蹤衛(wèi)星（satellite－satellite tracking，SST）數(shù)據(jù)分別以1mGal和1～2cm的精度恢復(fù)地球重力異常和大地水準面，同時要求空間分辨率達到100km，即重力場模型至少為200階次［1－2］。2013年10月，GOCE 衛(wèi)星結(jié)束使命，共測得約42個月的衛(wèi)星引力梯度數(shù)據(jù)。利用該類數(shù)據(jù)反演地球重力場模型的方法大致可分為3類：直接法、空域法和時域法［3－5］。直接法利用衛(wèi)星位置處觀測值與球諧系數(shù)的嚴密數(shù)學(xué)關(guān)系式進行解算?？沼蚍▽⒂^測值看作是衛(wèi)星位置的函數(shù)，可直接利用最小二乘法求解球諧系數(shù)，也可將觀測值歸算并內(nèi)插成平均軌道球面上的均勻格網(wǎng)值，進而利用解析法等方法進行解算。時域法將觀測值看作是沿軌時間序列，按解算方法的不同，又可分為時間域時域法和頻率域時域法，前者直接通過沿軌觀測值與球諧系數(shù)的線性關(guān)系利用最小二乘法進行解算，而后者又稱半解析法，即SA（semianalytical）方法［6］，根據(jù)觀測值與集總系數(shù)的關(guān)系式，通過1D－FFT或2D－FFT技術(shù)得到頻域內(nèi)的集總系數(shù)，再根據(jù)集總系數(shù)與球諧系數(shù)的線性關(guān)系用塊對角最小二乘法解算球諧系數(shù)，2D－FFT對應(yīng)的SA 方法稱為 Torus方法［6－7］。

國內(nèi)外學(xué)者對如何利用SA方法解算地球重力場模型進行了大量研究，快速傅里葉技術(shù)和塊對角最小二乘法的使用極大地提高了計算效率。對于利用1D－FFT的SA方法解算地球重力場模型已有大量相關(guān)研究，所用數(shù)據(jù)為GOCE引力梯度數(shù)據(jù)和SST 數(shù)據(jù)［8－11］。而對于利用 Torus方法解算模型的研究中，所用數(shù)據(jù)為模擬的GOCE軌道數(shù)據(jù)、衛(wèi)星引力梯度數(shù)據(jù)和SST數(shù)據(jù)，利用模擬引力梯度數(shù)據(jù)恢復(fù)模型的最高階次為120，受格網(wǎng)化誤差的影響，利用模擬GOCE衛(wèi)星引力梯度分量Vzz解算模型的精度優(yōu)于Vxx、Vyy分量的解算結(jié)果。Torus方法的解算精度與直接法相當(dāng)，但Torus方法所需時間僅為直接法的1%［12］。

本文介紹了利用衛(wèi)星引力梯度數(shù)據(jù)反演地球重力場模型的Torus法的理論和算法，利用GOCE實測軌道數(shù)據(jù)模擬無噪聲和含白噪聲情況下的衛(wèi)星引力梯度數(shù)據(jù)，分析Torus方法解算地球重力場模型的精度，極空白、歸算誤差、格網(wǎng)化誤差以及白噪聲等對結(jié)果的影響。

2 求解地球重力場模型的Torus方法

根據(jù)衛(wèi)星軌道攝動理論，局部軌道坐標(biāo)系（local orbit reference frame，LORF）下引力位與衛(wèi)星6個軌道根數(shù)之間的關(guān)系可表示為［13］

式中，a為衛(wèi)星軌道長半軸；I為軌道傾角；u為升交角距，u＝ω＋M，ω和M分別為近地點角距和平近點角；Λ為升交點經(jīng)度，Λ＝Ω－θG，其中Ω為升交點赤經(jīng)，θG為格林尼治恒星時角；μ為地心引力常數(shù)，即萬有引力常數(shù)與地球質(zhì)量的乘積；R為地球平均半徑；n、m為球諧函數(shù)的階和次，第3個加和下標(biāo)中的［2］表示k的變化步長為2。（I）和Gnkq（e）分別為正規(guī)化的傾角函數(shù)和偏心率函數(shù)，e為偏心率，ψmkq為u、Λ和M的函數(shù)，Ψmkq＝ku＋mΛ＋qM。球諧系數(shù)αnm和βnm隨n－m的奇偶性取不同的值，具體關(guān)系如下

對于e≈0的近圓軌道，式（1）可以簡化為［4，13－15］

式中，r為地心向徑；ψmk＝ku＋mΛ，此時，u＝ω＋f，f為真近點角。

式（4）中級數(shù)求和的n應(yīng)取至無窮大，但實際計算中常截止到某最高階N。經(jīng)過符號代換，式（4）可寫為［6］

式中，Amk、Bmk為集總系數(shù)，是轉(zhuǎn)換系數(shù)和球諧系數(shù)（αnm和βnm）的線性組合

轉(zhuǎn)換系數(shù)建立了球諧系數(shù)與集總系數(shù)的線性關(guān)系。6個引力梯度分量對應(yīng)的轉(zhuǎn)換系數(shù)如下［6］

式（5）是SA方法的基本表達式，是傅里葉級數(shù)的形式，集總系數(shù)Amk和Bmk為傅里葉系數(shù)。SA方法要求使用名義軌道上的觀測值，名義軌道具有3個特點：軌道為圓形、軌道傾角為常量、軌道受J2項影響而長周期進動［6－7］，同時，為便于FFT快速計算，要求等間隔采樣的觀測值。式（5）中變量u和Λ分別為定義在兩個不同方向的圓周，其變化范圍均為［0，2π），這兩個圓周可形成一個封閉的圓環(huán)，若觀測值均勻分布在該圓環(huán)面上，則可以利用2D－FFT方法解算集總系數(shù)，因此D－FFT對應(yīng)的SA方法稱為Torus方法。

采用Torus方法由GOCE沿軌衛(wèi)星引力梯度觀測值確定地球重力場模型的流程如圖1所示，其中ai表示第i次迭代的集總系數(shù)，Ki＋1為第i次迭代后的球諧系數(shù)。

圖1 利用Torus方法解算GOCE衛(wèi)星重力場模型的流程Fig.1 Flow chart for GOCE satellite gravity field determination using Torus approach

3 試驗與分析

選擇2009年11月01日至2009年12月31日，共61d（GOCE衛(wèi)星軌道重復(fù)周期）的SST＿PRD＿2數(shù)據(jù)產(chǎn)品，該產(chǎn)品提供了地固系（earth－fixed reference frame，EFRF）下采樣間隔為10s的簡化動力學(xué)軌道，計算6個軌道根數(shù)［13，16－17］，確定名義軌道的高度r0＝6 637.655km和傾角I0＝96.628°。由于在地球重力場模型解算中，引力梯度分量Vzz的影響較大，所以僅以該分量為例。

3.1 名義軌道格網(wǎng)引力梯度模擬數(shù)據(jù)

為驗證無誤差情況下Torus方法反演地球重力場模型的精度，分別利用EGM2008模型［18］的前60、80、90、110、150和200階次系數(shù)，采用2D－FFT方法由式（5）計算名義軌道格網(wǎng)點的引力梯度分量Vzz。

利用2D－FFT方法和按次排列的塊對角最小二乘法求解地球重力場模型系數(shù)，并計算模型的階誤差［19］，如圖2所示，其中60、110、200階模型的球諧系數(shù)誤差譜如圖3所示。僅從階誤差曲線來看，60階模型的階誤差小于10－21，而200階模型的階誤差達10－16，且隨著模型最高階次的增大，階誤差有增大的趨勢，且階誤差呈現(xiàn)奇偶性的差異。圖3中橫軸中心對應(yīng)的次為0，左右兩邊分別表示ˉSnm和ˉCnm的誤差譜。圖左右兩邊基本對稱，說明同階次的ˉCnm和ˉSnm的誤差基本一致。另外，較大的誤差集中在低次項系數(shù)，這主要是受GOCE任務(wù)極空白的影響［20］，且極空白的影響隨著最高階次的增大而增大，這是導(dǎo)致模型階誤差隨著最高階增加而增大的主要原因。從200階重力場模型的誤差譜可以看出，盡管受到極空白的影響，球諧系數(shù)誤差均小于10－16，說明當(dāng)觀測值滿足Torus方法的理論要求時可以完全恢復(fù)重力場模型。

圖2 不同階次模型的階誤差比較Fig.2 Comparisons of degree error among different models

圖3 最高階為60、110和200的模型誤差譜Fig.3 Error spectrum of models with maximum degree 60，110and 200

3.2 無誤差的沿軌GOCE衛(wèi)星引力梯度模擬觀測值

Torus方法要求使用圓環(huán)面均勻格網(wǎng)點的觀測值，而GOCE衛(wèi)星的軌道并非圓形，且傾角也不是常數(shù)，為了分析在GOCE衛(wèi)星實測軌道情況下Torus方法解算地球重力場的能力，利用EGM2008模型前200階次的系數(shù)模擬了局部軌道坐標(biāo)系下的沿軌衛(wèi)星引力梯度分量Vzz，將其作為實際觀測值，以驗證方法的正確性和嚴密性。

為滿足Torus方法的要求，GOCE衛(wèi)星的沿軌引力梯度觀測值應(yīng)歸算并內(nèi)插至名義軌道的格網(wǎng)點。設(shè)點P為GOCE衛(wèi)星軌道上的任意一點，其軌道高度和傾角分別為h和I，Pt為名義軌道上的點，其軌道高度和傾角分別為常數(shù)h0、I0。利用泰勒級數(shù)展開可以得到Pt點處的引力梯度分量Vij（h0，I0）為［12］

式中，Vij（h，I）為P點的引力梯度分量，下標(biāo)ij表示局部軌道坐標(biāo)系下引力梯度張量的6個分量（xx，xy，xz，yy，yz，zz），式（8）中引力梯度分量在軌道高度和傾角方向的改正量可以利用式（5）求偏導(dǎo)得到。

表1給出僅考慮引力梯度分量Vzz的情況下61d軌道對應(yīng)的式（8）中6項偏導(dǎo)分量的STD以及它們與Vzz分量的比值。其中，關(guān)于h的偏導(dǎo)分量較大，一階偏導(dǎo)分量達9.746E，占Vzz分量的3.580×10－3，二階偏導(dǎo)分量達26.84mE，與Vzz分量的比值為9.681×10－6，三階偏導(dǎo)分量小于1mE；關(guān)于I的偏導(dǎo)分量較小，其中一階偏導(dǎo)分量最大，為1.352mE，僅占Vzz分量的4.969×10－7，二階偏導(dǎo)分量和交叉偏導(dǎo)分量更小。相對于GOCE衛(wèi)星引力梯度分量Vzz中的噪聲而言，軌道高度方向三階以上的偏導(dǎo)分量，以及軌道傾角三階及三階以上的分量均可忽略。

表1 6項關(guān)于h和I的偏導(dǎo)改正量的大小以及與Vzz分量的比值Tab.1 Six derivatives and their percentage of Vzzwith respect to hand I

利用Kriging方法將歸算后的離散Vzz分量內(nèi)插至格網(wǎng)點，其中，半變異函數(shù)模型選用球狀模型，搜索半徑為1°。利用模型 EGM2008前200階次的系數(shù)采用2D－FFT方法模擬格網(wǎng)點的引力梯度值，并計算格網(wǎng)點引力梯度的差值，統(tǒng)計結(jié)果見表2，格網(wǎng)點引力梯度差值反映了歸算誤差和格網(wǎng)化誤差的影響，由表2可知Vzz的誤差約為0.110mE，與GOCE衛(wèi)星重力梯度儀精度指標(biāo)（5mE）相比，歸算和格網(wǎng)化的精度滿足要求。

表2 格網(wǎng)點的衛(wèi)星引力梯度觀測值與模擬值的差異Tab.2 Difference between the gridded GOCE satellite gradiometry observations and the simulations E

為了消除歸算誤差并減小格網(wǎng)化誤差，本文選擇EGM96為參考模型，并利用其前200階次的球諧系數(shù)模擬了局部軌道坐標(biāo)系下的沿軌引力梯度參考值（僅指Vzz分量），將殘差格網(wǎng)化，并利用2D－FFT和按次排列的塊對角最小二乘法求解球諧系數(shù)改正量，對參考模型進行改正，并將改正后的重力場模型作為新的參考模型進行迭代。

3次迭代后，除小于10次的球諧系數(shù)受極空白的影響改正量較大外，其他系數(shù)的改正量小于10－12。圖4和圖5分別顯示了模型的誤差階中值和誤差譜，受低次系數(shù)誤差較大的影響，低階部分的誤差階中值也較大。

圖4 無噪聲情況下模型階中值Fig.4 Degree median of models from simulated noise－free data

圖5 無噪聲情況下模型誤差譜Fig.5 Error spectrum of models from simulated noise－free data

為了分析Torus方法解算模型在不同階的大地水準面精度，利用式（9）和（10）計算模型的大地水準面階誤差和累積誤差

圖6給出了EGM2008模型的大地水準面信號、無噪聲時Torus方法解算的重力場模型相對于EGM2008的大地水準面階誤差和累積誤差。由于次小于10的系數(shù)受極空白影響較大，所以在計算大地水準面階誤差和累積誤差時未考慮低次部分。在格網(wǎng)化誤差的影響下，隨著階的增大，大地水準面階誤差有增大的趨勢，200階時階誤差達到最大值，為0.022mm，累積誤差為0.099mm，遠小于EGM2008模型的大地水準面信號。另外，相對于GOCE任務(wù)對大地水準面的精度要求（1～2cm），格網(wǎng)化誤差基本可以忽略。

圖6 無噪聲情況下模型的大地水準面階誤差和累積誤差Fig.6 Degree and cumulative geoid error of models from simulated noise－free data

3.3 有誤差的沿軌GOCE衛(wèi)星引力梯度模擬觀測值

為探究含白噪聲時Torus方法解算地球重力場模型的精度，按照GOCE衛(wèi)星任務(wù)的設(shè)計要求［1］，在模擬數(shù)據(jù)中加入了功率譜密度為5mE/Hz1/2的白噪聲，然后解算位系數(shù)，需要注意的是，利用無誤差的模擬觀測值解算地球重力場模型時，迭代可以減小格網(wǎng)化的影響，但當(dāng)觀測值中含5mE/Hz1/2白噪聲時，噪聲遠大于格網(wǎng)化誤差，所以不需要迭代。為了與Torus方法的解算結(jié)果作比較，利用相同觀測數(shù)據(jù)，采用空域最小二乘法解算地球重力場模型［21］，圖7顯示了兩種方法解算模型的誤差譜，相比空域最小二乘法，Torus方法受極空白影響相對更大，模型解算的精度略低，特別是140階以內(nèi)，空域最小二乘法的精度明顯更高。與無噪聲的結(jié)果相比，含白噪聲時誤差大了3～4個數(shù)量級，同時，極空白影響的范圍變大，主要是次小于13的系數(shù)誤差較大，誤差隨著階的增大而變大。

圖7 觀測值中含白噪聲時模型誤差譜Fig.7 Error spectrum of models from observations with white noise

為了進一步比較兩種方法的解算精度，分別利用式（9）和式（10）計算大地水準面階誤差和累積誤差，其中mmin＝13，計算結(jié)果見圖8，隨著階的增大，兩種方法解算模型的階誤差和累積誤差不斷增大，Torus方法和空域最小二乘法解算模型的累積誤差曲線分別在約193階和196階處與EGM2008模型曲線相交。受格網(wǎng)化誤差和極空白的影響，Torus方法的大地水準面階誤差和累積誤差略大于空域最小二乘法，140階以內(nèi)尤為明顯，這與圖7有良好的一致性。Torus模型的階誤差和累積誤差的最大值分別為1.58cm和6.37cm，而空域最小二乘法對應(yīng)的最大值分別為1.45cm和5.55cm。但就計算效率而言，利用空域最小二乘法解算200階重力場模型，使用106個CPU同時計算，所需時間約為9.4h，而由于采用了2D－FFT技術(shù)和按次排列的塊對角最小二乘法，使用1個CPU的情況下，Torus方法僅需要51min，計算效率大大提高。

目前發(fā)布的GOCE衛(wèi)星重力場模型中，go＿cons＿gcf＿2＿dir＿r1（直接法）和go＿cons＿gcf＿2＿tim＿r1（時域法）均聯(lián)合采用71d（2009－11－01—2010－01－11）的實測GOCE衛(wèi)星引力梯度數(shù)據(jù)和衛(wèi)星跟蹤衛(wèi)星數(shù)據(jù)［5］，不考慮次小于13的系數(shù)，計算兩個第1代模型的大地水準面階誤差和累積誤差，結(jié)果見圖8，200階時，直接解和時域解的大地水準面階誤差分別為1.06cm、2.63cm，累積誤差分別為7.00cm、12.27cm。

圖8 觀測值中含白噪聲時模型的大地水準面階誤差和累積誤差Fig.8 Degree and cumulative geoid error of models from observations with white noise

4 結(jié) 論

本文在介紹Torus方法解算地球重力場模型的基本原理和方法步驟的基礎(chǔ)上，利用GOCE衛(wèi)星實測軌道和EGM2008模型，分別模擬了名義軌道上的格網(wǎng)引力梯度值、實測軌道上無誤差和含5mE/Hz1/2白噪聲的GOCE衛(wèi)星引力梯度觀測值，利用這3類數(shù)據(jù)分別解算了200階次地球重力場模型，分析了無誤差、受格網(wǎng)化誤差影響以及觀測值中含白噪聲3種情況下模型的精度，結(jié)果表明在無誤差的情況下，Torus方法解算模型的精度較高，受極空白的影響，低次項系數(shù)誤差較大，且隨著模型最高階次的增大，階誤差有增大的趨勢。200階模型低階項的階誤差優(yōu)于10－16，隨著階的增大，階誤差減小，200階時，階誤差僅10－21；觀測值中不含白噪聲的情況下，迭代3次后，受極空白和格網(wǎng)化誤差的影響，次小于10的系數(shù)誤差較大，階中值達10－10。在不考慮低次系數(shù)的情況下，大地水準面階誤差隨著階的增大而增大，200階時，階誤差為0.022mm，累積誤差為0.099mm。觀測值中含有5mE/Hz1/2白噪聲時，解算模型比無噪聲模型的精度差了3～4個數(shù)量級。模型的大地水準面階誤差和累積誤差（未考慮次小于13的系數(shù)）隨著階的增大而增大，200階時，兩種誤差分別為1.58cm和6.37cm，而空域最小二乘解的誤差分別為1.45cm和5.55cm，受格網(wǎng)化誤差和極空白的影響，Torus方法的誤差略大于空域最小二乘法的誤差，但由于在Torus方法的解算過程中采用了2D－FFT技術(shù)和按次排列的塊對角最小二乘法，在計算效率方面，Torus方法明顯占優(yōu)。

模擬表明利用Torus方法可以高精度地恢復(fù)地球重力場，為利用GOCE引力梯度觀測值反演地球重力場模型提供了一種獨立有效的方法，由于采用了二維快速傅里葉技術(shù)和塊對角最小二乘法，極大地提高了解算效率。但受極空白的影響，重力場模型低次項的精度較差，可采用Kaula正則化等方法［22－24］解決。另外，GOCE衛(wèi)星實測引力梯度數(shù)據(jù)中存在有色噪聲，可采用移動恢復(fù)法或維納濾波方法［5，25］進行處理。為滿足Torus方法的要求，需將梯度儀坐標(biāo)系下的衛(wèi)星引力梯度觀測值轉(zhuǎn)換至局部軌道坐標(biāo)系，轉(zhuǎn)換過程中低精度分量會影響高精度分量的精度，可利用先驗?zāi)Ｐ突虻庥嬎愕挠^測值代替低精度分量，以減小坐標(biāo)轉(zhuǎn)換產(chǎn)生的影響。

［1］ESA.Gravity Field and Steady－state Ocean Circulation Mission.Reports for Mission Selection of the Four Candidate Earth Explorer Core Missions［R］.ESA Publications Division，ES SP－1233（1），1999.

［2］DRINKWATER M R，HAAGMANS R，MUZI D，et al.The GOCE Gravity Mission：ESA’s First Core Earth Explorer［R］.Proceedings of the 3rd International GOCE User Workshop，F(xiàn)rascati，Italy，ESA Special Publication，SP－627，ISBN 92－9092－938－3，2006：1－8.

［3］RUMMEL R，VAN GELDEREN M，KOOP R，et al.Spherical Harmonic Analysis of Satellite Gradiometry［M］.Netherlands Geodetic Commission：Publications on Geodesy，New Series 39，1993.

［4］KOOP R.Global Gravity Field Modeling Using Satellite Gravity Gradiometry［M］.Netherlands Geodetic Commission：Publications on Geodesy，New Series 38，1993.

［5］PAIL R，BRUINSMA S，MIGLIACCIO F，et al.First GOCE Gravity Field Models Derived by Three Different Approaches［J］.Journal of Geodesy，2011，85（11）：819－843.

［6］SNEEUW N J.A Semi－analytical Approach to Gravity Field Analysis from Satellite Observations［D］.Munich，Germany：Institut für Astronomische und Physikalische Geod?sie，Technische Universit?t München，2000.

［7］SNEEUW N J.Space－wise，Time－wise，Torus and Rosborough Representations in Gravity Field Modeling［J］.Space Science Reviews，2003，108（1－2）：37－46.

［8］KLEES R，DITMAR P.The Performance of the Time－wise Semi－analytical Inversion of Satellite Gravity Gradients［M］∥áDáM J，SCHWARZ K P.Vistas for Geodesy in the New Millennium，International Association of Geodesy Symposia.Berlin Heidelberg：Springer，2001，125：253－258.

［9］SCHUH W D，PAIL R，PLANK G.Assessment of Different Numerical Solution Strategies for Gravity Field Recovery［C］∥Proceedings of the 1st International GOCE User Workshop，ESA WPP－188，87－95，ESA/ESTEC，2001.

［10］PAIL R，WERMUTH M.GOCE SGG and SST Quicklook Gravity Field Analysis［J］.Advances in Geosciences，2003，1：5－9.

［11］PAIL R，PLANK G.GOCE Gravity Field Processing Strategy［J］.Studia Geophysica et Geodaetica，2004，48（2）：289－309.

［12］XU Chen.The Torus－based Semi－analytical Approach in Spaceborne Gravimetry［D］.Calgary：Department of Geomatics Engineering，University of Calgary，2008.

［13］KAULA W M.Theory of Satellite Geodesy：Applications of Satellites to Geodesy［M］.Waltham Massachusetts：Blaisdell Publishing Company，1966.

［14］SCHRAMA E J O.The Role of Orbit Errors in Processing of Satellite Altimeter Data［M］.Netherlands Geodetic Commission：Publications on Geodesy，New Series 33，1989.

［15］PAIL R，PLANK G.Assessment of Three Numerical Solution Strategies for Gravity Field Recovery from GOCE Satellite Gravity Gradiometry Implemented on a Parallel Platform［J］.Journal of Geodesy，2002，76（8）：462－474.

［16］ESA.GOCE HPF：GOCE Level 2Product Data Handbook［R］.Technical Note，GO－MA－HPF－GS－0110，2010.

［17］CAPITAINE N，WALLACE P T，MCCARTHY D D.Expressions to Implement the IAU 2000Definition of UT1［J］.Astronomy ＆ Astrophysics，2003，406（3）：1135－1149.

［18］PAVLIS N K，HOLMES S A，KENYON S C，et al.An Earth Gravitational Model to Degree 2160：EGM2008［J］.EGU General Assembly，2008，10：13－18.

［19］ESA.GOCE L1bProducts User Handbook［R］.Technical Note，GOCE－GSEG－EOPGTN－06－0137，2006.

［20］SNEEUW N J，VAN GELDEREN M.The Polar Gap［M］∥SANSóF，RUMMEL R.Geodetic Boundary Value Problems in View of the One Centimeter Geoid.Lecture Notes in Earth Sciences.Berlin Heidelberg：Springer，1997，65：559－568.

［21］XU Xinyu，LI Jiancheng，JIANG Weiping，et al.Simulation Study for Recovering GOCE Satellite Gravity Model Based on Space－wise LS Method［J］.Acta Geodaetica et Cartographica Sinica，2011，40（6）：697－702.（徐新禹，李建成，姜衛(wèi)平，等.基于空域最小二乘法求解GOCE衛(wèi)星重力場的模擬研究［J］.測繪學(xué)報，2011，40（6）：697－702.）

［22］RUDOLPH S，KUSCHE J，IIK K H.Investigations on the Polar Gap Problem in ESA’s Gravity Field and Steady－state Ocean Circulation Explorer Mission（GOCE）［J］.Journal of Geodynamics，2002，33（1－2）：65－74.

［23］XU Xinyu，LI Jiancheng，WANG Zhengtao，et al.The Simulation Research on the Tikhonov Regularization Applied in Gravity Field Determination of GOCE Satellite Mission［J］.Acta Geodaetica et Cartographica Sinica，2010，39（5）：465－470.（徐新禹，李建成，王正濤，等.Tikhonov正則化方法在GOCE重力場求解中的模擬研究［J］.測繪學(xué)報，2010，39（5）：465－470.）

［24］ZHU Guangbin，LI Jiancheng，WEN Hanjiang，et al.Slepian Localized Spectral Analysis of the Determination of the Earth’s Gravity Field Using Satellite Gravity Gradiometry Data［J］.Acta Geodaetica et Cartographica Sinica，2012，41（1）：1－7.（朱廣彬，李建成，文漢江，等.衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場的Slepian局部譜分析方法［J］.測繪學(xué)報，2012，41（1）：1－7.）

［25］BAUR O，SNEEUW N，GRAFAREND E W.Methodology and Use of Tensor Invariants for Satellite Gravity Gradiometry［J］.Journal of Geodesy，2008，82（4－5）：279－293.