徐新禹，李建成，姜衛(wèi)平，鄒賢才

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

基于空域最小二乘法求解GOCE衛(wèi)星重力場的模擬研究

徐新禹1，2，李建成1，2，姜衛(wèi)平1，2，鄒賢才1，2

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

論述最小二乘過程中有色噪聲的處理方法，提出使用自回歸模型對GOCE梯度觀測值中的有色噪聲進(jìn)行時(shí)域?yàn)V波，數(shù)值模擬結(jié)果驗(yàn)證該方法的有效性。利用數(shù)值模擬驗(yàn)證直接求逆方法和PCCG法求解大型法方程的有效性，后者的效率遠(yuǎn)遠(yuǎn)高于前者。聯(lián)合加入噪聲（有色噪聲和白噪聲）的衛(wèi)星重力梯度張量徑向分量觀測值Vzz和SST觀測值，分別使用空域最小二乘法和半解析法恢復(fù)180階全球重力場模型，前者求解重力場模型的大地水準(zhǔn)面和重力異常在180階次的精度分別為3.01 cm和0.75 mGal（1 mGal＝10－5m／s2），優(yōu)于半解析法求解模型的精度。

GOCE衛(wèi)星；空域最小二乘法；預(yù)處理共軛梯度法；半解析法；衛(wèi)星重力梯度；自回歸模型

1 引 言

21世紀(jì)初成功實(shí)施的CHAMP（challenging minisatellite payload）和GRACE（gravity recovery and climate experiment）衛(wèi)星重力探測計(jì)劃以前所未有的精度獲得了重力場的中、長波信息，為實(shí)現(xiàn)本世紀(jì)大地測量的宏偉目標(biāo)——厘米級大地水準(zhǔn)面奠定了堅(jiān)實(shí)的基礎(chǔ)［1］。但由于已實(shí)施的衛(wèi)星重力探測模式以及重力場信號隨高度上升快速衰減特性的限制，恢復(fù)重力場模型的精度和分辨率仍然很難滿足許多地球相關(guān)學(xué)科的科學(xué)目標(biāo)需求［2］。ESA負(fù)責(zé)實(shí)施的GOCE（gravity field and steady－state ocean circulation explorer）衛(wèi)星于2009年3月17日成功發(fā)射（http：∥www.esa. int／esaLP／ESAYEK1VMOC＿LPgoce＿0.html），其采用SGG（satellite gravity gradiometry）技術(shù)和SST－h(huán)l（satellite－to－satellite tracking in the high－low mode）技術(shù)相結(jié)合的測量模式，SGG技術(shù)的觀測量為引力位的二階導(dǎo)數(shù)，能在一定程度上補(bǔ)償重力場信號隨高度產(chǎn)生的衰減，因此期望獲得更高精度和分辨率的重力場模型［3］。

自20世紀(jì)70年代以來，國內(nèi)外學(xué)者和機(jī)構(gòu)對利用衛(wèi)星重力梯度觀測值恢復(fù)地球重力場的相關(guān)理論與方法進(jìn)行了大量的研究，并進(jìn)行了大量的數(shù)值分析與驗(yàn)證［4－15］，主要是從空域法和時(shí)域法出發(fā)討論?？沼蚍ò延^測值看做位置的函數(shù)，其解算方法主要有積分法、最小二乘配置法、最小二乘譜組合法、空域最小二乘法［4－5］。時(shí)域法是將觀測值看做沿軌道的時(shí)間序列，使用傾角函數(shù)表示球諧函數(shù)，解算方法一般有兩種：時(shí)間域時(shí)域法和半解析法（semi－analytical，SA）［4－5，9，11］?？沼蜃钚《朔ㄖ苯永醚匦l(wèi)星軌道的觀測值建立觀測方程，由經(jīng)典的最小二乘方法求解。該方法不對數(shù)據(jù)進(jìn)行歸算、格網(wǎng)化等處理，同時(shí)不需考慮數(shù)據(jù)中斷、跳躍等問題，因而可獲得較精確的結(jié)果，并提供求解參數(shù)完全的方差－協(xié)方差矩陣，是一種較為嚴(yán)密的方法。但該方法在處理海量觀測值，且觀測值中包含有色噪聲時(shí)，恢復(fù)高階次重力場模型對計(jì)算機(jī)性能要求較高。空域最小二乘法是ESA最終利用GOCE衛(wèi)星任務(wù)實(shí)測數(shù)據(jù)確定高精度高分辨率全球重力場的方法之一［13－14］。本文主要研究其由SGG觀測數(shù)據(jù)求解高階衛(wèi)星重力場模型的相關(guān)理論、方法及算法。

2 數(shù)學(xué)模型

2.1 引力位的二階導(dǎo)數(shù)

地心球坐標(biāo)系下引力位的二階導(dǎo)數(shù)f（r，θ，λ）可表示成傅里葉級數(shù)形式［5］

式中，

r、θ、λ分別表示地心球坐標(biāo)的半徑、地心余緯（90°－地心緯度φ）、地心經(jīng)度；n、m是球諧展開的階和次；ˉCnm、ˉSnm是正規(guī)化的n階m次引力位球諧系數(shù)；ˉPnm（cosθ）為完全規(guī)格化締合勒讓德函數(shù)；Afm（r，θ）、Bfm（r，θ）為傅里葉系數(shù)；Hfnm（r，θ）為其與引力位系數(shù)之間的轉(zhuǎn)換系數(shù)，不同引力位函數(shù)對應(yīng)不同形式的傅里葉系數(shù)和轉(zhuǎn)換系數(shù)。根據(jù)局部指北坐標(biāo)系（x，y，z）的定義（x軸指向北，y軸指向西，z軸徑向朝外）及其與地心球坐標(biāo)（r，θ，λ）之間的關(guān)系，可推導(dǎo)出其中引力梯度張量球諧級數(shù)展開式中傅里葉系數(shù)和轉(zhuǎn)換系數(shù)的形式，如表1所示，其中為地球萬有引力常數(shù)與地球質(zhì)量的乘積，R是地球平均半徑。本文的模擬試算均在局部指北坐標(biāo)系下進(jìn)行，該坐標(biāo)系是GOCE任務(wù)實(shí)測梯度數(shù)據(jù)處理中采用的坐標(biāo)系之一。

表1 局部指北坐標(biāo)系中引力梯度張量的6個(gè)分量所對應(yīng)傅里葉系數(shù)和轉(zhuǎn)換系數(shù)的形式Tab.1 Fourier coefficients and transformation coefficients of the six components of gravitational gradient tensor in the local north－oriented coordinate system

2.2 空域最小二乘法的數(shù)學(xué)模型及有色噪聲處理方法

從式（1）可看出衛(wèi)星引力梯度觀測值是衛(wèi)星位置和重力場模型參數(shù)的函數(shù)。從嚴(yán)格意義上說，衛(wèi)星位置和重力場參數(shù)都是未知的，因此觀測方程是非線性的，在實(shí)際使用中需要線性化以便同時(shí)求解位置和重力場參數(shù)。考慮到GOCE任務(wù)的精密定軌技術(shù)已經(jīng)可達(dá)厘米水平，經(jīng)驗(yàn)證5cm的軌道誤差對引力梯度張量觀測值影響的最大量級為0.04mE，這相對于GOCE衛(wèi)星重力梯度儀在測量帶寬內(nèi)3mE的精度可忽略不計(jì)［3］，因此可認(rèn)為衛(wèi)星位置已知。在引入正常引力場的情況下，可建立如下觀測方程

式中，ΔΓij（P）表示引力梯度異常，是引力位梯度與正常位梯度之差；Tij（P）為擾動引力梯度，εij表示觀測誤差以及模型誤差；P表示精密定軌技術(shù)確定的衛(wèi)星位置。式（2）中僅包含重力場位系數(shù)參數(shù)，可由衛(wèi)星軌道上每個(gè)點(diǎn)的觀測值形成方程組，將Tij（P）展開如式（1）的形式，則式（2）對應(yīng)的線性Gauss－Markov觀測模型和統(tǒng)計(jì)模型為

式中，y表示梯度觀測值向量，既可是引力梯度異常的單個(gè)分量，也可是多個(gè)分量的組合；A為設(shè)計(jì)矩陣；x表示位系數(shù)參數(shù)向量；ε為觀測誤差向量；為單位權(quán)方差；Q為觀測值的誤差協(xié)方差矩陣；P為觀測值的權(quán)。根據(jù)經(jīng)典線性最小二乘原理可求得式（3）的最優(yōu)無偏估計(jì)解

GOCE任務(wù)梯度觀測值中包含有色噪聲，將導(dǎo)致形成的誤差協(xié)方差矩陣Q不再是對角陣，且隨著觀測量的增多，該矩陣維數(shù)不斷增大，給直接處理協(xié)方差矩陣（求逆、乘積等）帶來困難，因而需要找到合適的方法處理有色噪聲帶來的計(jì)算問題。根據(jù)GOCE任務(wù)重力梯度測量誤差的特性，誤差時(shí)間序列為廣義平穩(wěn)隨機(jī)信號，若假設(shè)觀測值為周期的等間隔（空間）采樣或者為格網(wǎng)點(diǎn)值，則Q是一個(gè)循環(huán)Toeplitz矩陣，但是Q的逆將破壞其原有的循環(huán)Toeplitz結(jié)構(gòu)。若將Q分解為兩個(gè)三角矩陣的乘積Q＝RTR，根據(jù)矩陣運(yùn)算規(guī)律（RTR）－1＝R－1（R－1）T、式（3）和式（4）可推導(dǎo)出如下法方程

該卷積過程可看做一個(gè)離散線性移不變系統(tǒng)，向量f對應(yīng)濾波器的脈沖響應(yīng)。該過程可利用FFT算法快速計(jì)算（時(shí)域卷積在頻域內(nèi)為簡單的乘積），首先將時(shí)域內(nèi)的時(shí)間序列轉(zhuǎn)換到頻域內(nèi)，濾波后再轉(zhuǎn)換到時(shí)域內(nèi)，該方法可為頻域?yàn)V波；根據(jù)信號處理理論中的相關(guān)定理，信號的自相關(guān)函數(shù)與PSD（power spectral density）成傅里葉變換對，因此還可通過噪聲特性構(gòu)造時(shí)域?yàn)V波器，例如基于AR（autoregressive）模型的高通濾波器，這樣空域最小二乘中有色噪聲的處理過程就變成了一個(gè)時(shí)域?yàn)V波過程。

需要指出，求解高階重力場模型的最小二乘方法非常耗時(shí)，尤其是法方程的解算，本文將對直接求逆法和PCCG（preconditioned conjugate gradient）進(jìn)行模擬比較，前者是對法方程矩陣嚴(yán)格求逆，優(yōu)點(diǎn)是可以獲得參數(shù)的方差－協(xié)方差矩陣（可表征求解參數(shù)的精度），后者通過迭代快速解算法方程，但無法嚴(yán)格給出求解參數(shù)的精度信息。

3 數(shù)值分析

3.1 濾波方法的分析與驗(yàn)證

為驗(yàn)證基于AR模型的時(shí)域?yàn)V波方法的有效性，本文基于先驗(yàn)誤差PSD模擬了30d5s采樣的有色噪聲時(shí)間序列，其統(tǒng)計(jì)特性如表2所示，統(tǒng)計(jì)量的單位為E（E?tv?s），先驗(yàn)誤差PSD的解析形式為［16］

表2 模擬有色噪聲誤差時(shí)間序列的統(tǒng)計(jì)結(jié)果Tab.2 Statistical results of the simulated colored noise time series

本文取先驗(yàn)誤差PSD的倒數(shù)來構(gòu)造AR高通濾波器，AR模型的最高階取8 640，首先用先驗(yàn)解析形式的誤差PSD的倒數(shù)計(jì)算自相關(guān)函數(shù)，再利用得到的自相關(guān)函數(shù)計(jì)算AR模型的參數(shù)以形成時(shí)域?yàn)V波器，利用其對模擬的噪聲時(shí)間序列進(jìn)行濾波處理，為分析濾波的效果，驗(yàn)證該方法的有效性，圖1、圖2和圖3分別為濾波前后噪聲時(shí)間序列、分布圖和協(xié)方差函數(shù)圖。從圖中可明顯看出，濾波前噪聲時(shí)間序列具有一定的相關(guān)性，且不具備正態(tài)分布特性，時(shí)間序列曲線呈現(xiàn)一定的規(guī)律性，具有明顯的有色噪聲特性；濾波后噪聲時(shí)間序列呈現(xiàn)正態(tài)分布，且不存在相關(guān)性，其時(shí)間序列分布也顯示出隨機(jī)噪聲的特性。結(jié)果充分說明采用構(gòu)造的AR濾波器對GOCE梯度觀測值有色噪聲的濾波非常有效，濾波過程是白噪聲化的過程，即消除噪聲時(shí)間序列的相關(guān)性，這與2.2節(jié)中的理論分析一致。

圖1 濾波前（左）、后（右）的噪聲時(shí)間序列Fig.1 Noise time series before（left）and after（right）filtering

圖2 濾波前（左）、后（右）噪聲時(shí)間序列的分布Fig.2 Distribution of noise time series before（left）and after（right）filtering

圖3 濾波前（左）、后（右）噪聲時(shí)間序列的協(xié)方差函數(shù)Fig.3 Covariance function of noise time series before（left）and after（right）filtering

3.2 恢復(fù)180階重力場模型的數(shù)值分析

為了分析空域最小二乘法求解GOCE衛(wèi)星重力場的精度及運(yùn)行效率，并將其與SA方法進(jìn)行比較，本文模擬了GOCE任務(wù)29d的觀測數(shù)據(jù)，采樣間隔為5s，EGM96為軌道模擬的參考模型，軌道傾角約為96.7°，偏心率為0.001，軌道的平均高度為245km，是一個(gè)近似重復(fù)軌道，選擇GFZ發(fā)布的EIGEN＿CG03C重力場模型模擬引力梯度張量觀測值，模型最高階取至180。SST觀測值的模擬采用EIGEN＿CG03C模型計(jì)算沿模擬衛(wèi)星軌道的引力位來實(shí)現(xiàn)，模型的最高階為80階，加入方差為0.5m2s－2的白噪聲，其相應(yīng)于方差為5cm的位置誤差和10－4m／s的速度誤差。下面的數(shù)值分析僅以重力梯度的徑向分量Vzz為例，其在引力梯度張量中量級最大，加入有色噪聲的統(tǒng)計(jì)特性如表2所示。

在使用空域最小二乘法求解地球重力場模型之前，首先對直接求逆法和PCCG方法求解法方程的有效性及軟件模塊性能進(jìn)行測試，測試平臺為IBM P5－575，該小型機(jī)共有24個(gè)64位1.9GHz的POWER5CPU，測試數(shù)據(jù)選用不包含誤差的引力梯度張量觀測值的徑向分量Vzz。需要指出本文所采用的直接求逆法和PCCG方法均是在已形成法方程的基礎(chǔ)上進(jìn)行的，對于由大量的觀測數(shù)據(jù)求解高階重力場模型來說，形成法方程也是一個(gè)非常耗時(shí)的過程，根據(jù)法方程的計(jì)算特性，可采用多線程提高計(jì)算效率?；谥苯忧竽娣椒ê蚉CCG方法由無誤差的觀測值求得模型的階誤差均方差如圖4所示，圖中實(shí)線表示EIGEN＿CG03C模型的信號階次方差均方差（下同），點(diǎn)虛線表示直接求逆得到模型的階誤差均方差，短劃線虛線表示PCCG方法求得模型的階誤差均方差，前者計(jì)算需要約13d，后者僅需30min經(jīng)過25次迭代即可收斂。需要指出本文求解結(jié)果的比較均是相對于EIGEN＿CG03C模型，對于模擬試算來說，其相當(dāng)于“真”實(shí)場，階誤差均方差的計(jì)算公式為

式中，n，m為模型的階次；上標(biāo)‘EST’表示求解的重力場模型，‘REF’表示“真”實(shí)重力場模型。

從圖4中可看出，兩種方法均能以較高的精度（10－15～10－13）恢復(fù)位系數(shù)，求解模型的大地水準(zhǔn)面和重力異常誤差均方差在緯度±80°范圍內(nèi)量級分別為10－3cm和10－4mGal，在緯度±90°范圍內(nèi)量級分別為10－2cm和10－2mGal，這說明求解結(jié)果受到了極區(qū)無觀測數(shù)據(jù)的影響，導(dǎo)致求解過程不穩(wěn)定。從求解結(jié)果還可看出直接求逆方法比PCCG方法更精確，但PCCG方法的速度要快很多，雖然精度略有降低，但足以滿足求解的精度要求。需要指出，運(yùn)算中并未使用任何正則化技術(shù)，圖中的階誤差均方差在接近180階處有跳躍，應(yīng)該是高階重力場數(shù)值不穩(wěn)定性引起的。

圖4 利用直接求逆方法和PCCG方法求得模型的階誤差均方差Fig.4 Degree error RMS of the gravity model estimated from the direct inverse and the PCCG method

基于包含誤差（有色噪聲和白噪聲）的引力梯度張量徑向分量Vzz和SST觀測值，分別利用SA方法和空域最小二乘法恢復(fù)重力場模型的階誤差均方差，如圖5所示，其中點(diǎn)虛線表示SA方法求得模型的階誤差均方差，短劃線虛線表示采用了AR濾波的空域最小二乘法求得模型的階均方差，短劃線－點(diǎn)虛線表示未采用濾波的空域最小二乘法求得模型的階均方差，對于GOCE衛(wèi)星高階地球重力場的確定，由于極空白、重力場信號向下延拓的誤差放大、有色噪聲等原因都將使最小二乘的法方程具有較強(qiáng)的病態(tài)性，需要采用正則化技術(shù)來穩(wěn)定求解過程，本文采用Kaula正則化方法。形成法方程系數(shù)陣時(shí)采用時(shí)域?yàn)V波策略，分別采用PCCG方法和直接求逆法求解，經(jīng)比較求解結(jié)果，發(fā)現(xiàn)在觀測值包含誤差的情況下，兩種方法解算的結(jié)果幾乎一致，這進(jìn)一步表明PCCG方法求解法方程過程的快速有效性。限于篇幅，本文未給出SA方法的數(shù)學(xué)模型及解算方法，其具體實(shí)現(xiàn)及模擬計(jì)算結(jié)果可參閱文獻(xiàn)［17］。本文所用觀測值聯(lián)合求解的方法是通過法方程矩陣加權(quán)疊加實(shí)現(xiàn)的，即每一類觀測值單獨(dú)形成法方程矩陣，最后進(jìn)行法方程疊加。

從圖中可看出，采用加權(quán)策略（使用AR濾波實(shí)現(xiàn)），可大大改善求解結(jié)果的精度，在整個(gè)頻譜范圍內(nèi)均有很大提高，這也說明衛(wèi)星重力梯度觀測值中大的低頻誤差不僅影響重力場求解的低頻部分。采用空域最小二乘方法求解的結(jié)果優(yōu)于SA方法的結(jié)果，位系數(shù)精度在20到120階的范圍有明顯提高，顯示出空域最小二乘法求解過程的嚴(yán)密性。為了進(jìn)一步比較求解模型的精度，表3給出了兩種方法求解模型的大地水準(zhǔn)面和重力異常誤差統(tǒng)計(jì)結(jié)果。從表中可看出，不論是全球還是局部范圍內(nèi)，空域最小二乘法求解模型的精度均優(yōu)于SA方法，尤其是全球的統(tǒng)計(jì)結(jié)果，這說明了空域最小二乘法較SA方法求解更穩(wěn)定，受極空白區(qū)域的影響相對較小。雖然空域最小二乘法求解模型的精度較SA方法有提高，但并沒有量級上的差異，大地水準(zhǔn)面在±80°范圍內(nèi)均為幾個(gè)厘米，重力異常的精度在1mGal以內(nèi)，基本達(dá)到GOCE任務(wù)的目標(biāo)要求。

圖5 由空域最小二乘法和SA方法求解模型的階誤差均方差Fig.5 Degree error RMS of the gravity model estimated from space－wise LS method and SA approach

表3 由空域最小二乘法和SA方法求解模型的大地水準(zhǔn)面和重力異常誤差統(tǒng)計(jì)Tab.3 Statistics of geoid error and gravity anomaly error of the gravity model estimated from space－wise LS method and SA approach

此外，在本文模擬計(jì)算的情況下，空域最小二乘法以11個(gè)CPU同時(shí)計(jì)算大約需兩周時(shí)間，同樣的求解條件下，SA方法使用1個(gè)CPU僅需3h，所需內(nèi)存也遠(yuǎn)小于空域最小二乘法，可見SA方法是一種快速有效的解算方法，但前面的分析比較是在觀測值具有連續(xù)采樣、近似重復(fù)等特性的基礎(chǔ)上討論的，若考慮數(shù)據(jù)有間斷、非嚴(yán)格周期等情況，SA方法受到模型誤差的影響會增大，此時(shí)可以預(yù)見空域最小二乘法相比SA方法應(yīng)該有更好的穩(wěn)定性。

4 結(jié) 論

本文研究了空域最小二乘法求解高階GOCE衛(wèi)星重力場的相關(guān)理論、方法與算法，提出基于AR模型對最小二乘過程中的有色噪聲進(jìn)行濾波處理，數(shù)值模擬結(jié)果說明了該方法的有效性。直接求逆法和PCCG法用于法方程求解的數(shù)值比較結(jié)果說明兩者均能獲得較高精度的解，PCCG效率更高，但無法提供完全的方差－協(xié)方差矩陣?；诤姓`差（有色噪聲和白噪聲）的GOCE衛(wèi)星任務(wù)模擬數(shù)據(jù)Vzz和SST，分別利用空域最小二乘法和SA方法分別恢復(fù)了180階全球重力場模型，前者求解模型的精度優(yōu)于后者，但SA方法所需的計(jì)算時(shí)間和資源都遠(yuǎn)小于空域最小二乘法。在未知實(shí)測數(shù)據(jù)真實(shí)特性的情況下，很難預(yù)見SA方法處理實(shí)測數(shù)據(jù)的性能，空域最小二乘法是一種可用于GOCE實(shí)測數(shù)據(jù)處理的穩(wěn)定和較為嚴(yán)密的方法。需要指出，空域最小二乘法的快速解算問題需要作進(jìn)一步研究，在當(dāng)前擁有的計(jì)算機(jī)資源的條件下，若要對將來實(shí)測數(shù)據(jù)進(jìn)行處理（更長周期、多種類型的引力梯度觀測值，包括Vxx、Vyy等），研究并行算法在空域最小二乘解算（法方程形成、濾波、求逆以及正則化等運(yùn)算）中的應(yīng)用非常必要。

［1］ ZOU Xiancai.Theory of Satellite Orbit and Earth Gravity Field Determination［D］.Wuhan：Wuhan University，2007.（鄒賢才.衛(wèi)星軌道理論與地球重力場模型的確定［D］.武漢：武漢大學(xué)，2007.）

［2］ XU Xinyu，LI Jiancheng，ZOU Xiancai，et al.GOCE Gravity Explorer Mission［J］.Journal of Geodesy and Geodynamics，2006，26（4）：49－55.（徐新禹，李建成，鄒賢才，等.GOCE衛(wèi)星重力探測任務(wù)［J］.大地測量與地球動力學(xué)，2006，26（4）：49－55.）

［3］ ESA.Gravity Field and Steady－state Ocean Circulation Mission：Report for Mission Selection，the Four Candidate Earth Explorer Missions［R］.Noordwijk：ESA Publications Division，1999.

［4］ RUMMEL R，VAN GELDEREN M，KOOP R，et al.Spherical Harmonic Analysis of Satellite Gradiometry［R］.Delft：Netherlands Geodetic Commission，1993.

［5］ KOOP R.Global Gravity Field Modeling Using Satellite Gravity Gradiometry［R］.Delft：Netherlands Geodetic Commission，1993.

［6］ LUO Zhicai.The Theory and Methodology of Determining the Earth’s Gravity Field Using Satellite Gravity Gradiometry Data［D］.Wuhan：Wuhan Technical University of Surveying and Mapping，1996.（羅志才.利用衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場的理論和方法［D］.武漢：武漢測繪科技大學(xué)，1996.）

［7］ TSCHERNING C C.Computation of Spherical Harmonic Coefficients and Their Error Estimates Using Least Squares Collocation［J］.Journal of Geodesy，2001，75（1）：12－18.

［8］ WENZEL H G.Geoid Computation by Least Squares Spectral Combination Using Integral Kernels［C］∥Proceedings of the International Association of Geodesy General Meeting.Tokyo：IAG，1982：438－453.

［9］ KLEES R，KOOP R，VISSER P，et al.Efficient Gravity Field Recovery from GOCE Gravity Gradient Observations［J］.Journal of Geodesy，2000，74（7－8）：561－571.

［10］ SCHUH W D，PAIL R，PLANK G.Assessment of Different Numerical Solution Strategies for Gravity Field Recovery［C］∥Proceedings of the“International GOCE User Workshop”.Noordwijk：ESA Publications Division，2001：87－95.

［11］ SNEEUW N J.A Semi－analytical Approach to Gravity Field Analysis from Satellite Observations［D］.Munich： Institute of Astronomy and Physical Geodesy，Technical University of Munich，2000.

［12］ 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.

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

［14］ RUMMEL R，GRUBER T，KOOP R.High Level Processing Facility for GOCE：Products and Processing Strategy［C］∥Proceedings of the Second International GOCE User Workshop，"GOCE，The Geoid and Oceanography".Frascati：ESA－ESRIN，2004.

［15］ 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（4）：465－470.（徐新禹，李建成，王正濤，等.Tikhonov正則化方法在GOCE重力場求解中的模擬研究［J］.測繪學(xué)報(bào)，2010，39（4）：465－470.）

［16］ XU Xinyu，WANG Zhengtao，ZOU Xiancai，et al.The Research on Analysis and Simulation of Satellite Gravity Gradiometry Error of GOCE［J］.Journal of Geodesy and Geodynamics，2010，30（2）：71－75.（徐新禹，王正濤，鄒賢才，等.GOCE衛(wèi)星重力梯度測量誤差分析及其模擬研究［J］.大地測量與地球動力學(xué)，2010，30（2）：71－75.）

［17］ XU Xinyu，LI Jiancheng，JIANG Weiping，et al.The Fast Analysis of GOCE Gravity Field［C］∥International Association of Geodesy Symposia：Observing Our Changing Earth.Berlin：Springer，2009，133：379－385.

Simulation Study for Recovering GOCE Satellite Gravity Model Based on Space－wise LS Method

XU Xinyu1，2，LI Jiancheng1，2，JIANG Weiping1，2，ZOU Xiancai1，2
1.School of Geodesy and Geomatics，Wuhan University，Wuhan 430079，China；2.Key Laboratory of Geospace Environment and Geodesy Ministry of Education，Wuhan University，Wuhan 430079，China

The approach to deal with the colored noise in space－wise LS is discussed.The autoregressive（AR）filter is proposed to process the colored noise in the GOCE satellite gravity gradiometry，which is validated by the numerical results.The preconditioned conjugate gradient（PCCG）method and the direct inverse for solving the large normal equation are also verified by the numerical analysis.And the former approach is much faster than the latter one.Combing the radial component Vzzof the simulated gravity gradient tensor with colored noise and SST observation with white noise，the global gravity field model up to degree and order 180 is estimated by the spacewise LS method and the semi－analytical（SA）approach respectively.The accuracies of geoid and gravity anomaly of the global gravity field model estimated from the space－wise LS method are 3.01 cm and 0.75 mGal（1 Gal＝10－2m／s2）respectively，which are better than them from the SA approach.

GOCE satellite；space－wise LS method；preconditioned conjugate gradient（PCCG）；semi－analytical（SA）；satellite gravity gradiometry；autoregressive（AR）model

XU Xinyu（1978－），male，PhD，lecturer，majors in satellite geodesy.

1001－1595（2011）06－0697－06

P223

A

國家自然科學(xué)基金（40904003；40804004；40637034）；武漢大學(xué)自主科研項(xiàng)目（4082011）

雷秀麗）

2010－12－06

2011－02－23

徐新禹（1978－），男，博士，講師，主要研究方向?yàn)樾l(wèi)星大地測量。

E－mail：xyxu＠sgg.whu.edu.cn