朱廣彬，李建成，文漢江，常曉濤，王正濤，鄒賢才

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

衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場的Slepian局部譜分析方法

朱廣彬1，李建成2，文漢江3，常曉濤1，王正濤2，鄒賢才2

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

在引入Slepian局部譜分析方法的基礎(chǔ)上，詳細分析Slepian函數(shù)的數(shù)學(xué)特性，采用Grünbaum算子提高Slepian方法求解的穩(wěn)定性和效率，推導(dǎo)衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場的Slepian方法表達式。通過仿真分析，就Slepian方法在衛(wèi)星重力梯度數(shù)據(jù)確定地球重力位模型中的應(yīng)用和前景進行分析和討論。研究表明，Slepian函數(shù)在整個球面和球帶上具有雙正交性，其頻譜能量分布特性與衛(wèi)星軌道的測量特點具有很好的一致性。Slepian低次項系數(shù)精度受到極空白影響很小，較之球諧系數(shù)低次項明顯改善。Slepian方法對大地水準(zhǔn)面空間分布恢復(fù)精度的直接貢獻不明顯。

衛(wèi)星重力梯度；重力場；Slepian局部譜分析；極空白；面球諧函數(shù)；Slepian函數(shù)

1 引 言

隨著空間技術(shù)的發(fā)展，利用衛(wèi)星重力技術(shù)獲取地球重力場精細結(jié)構(gòu)的探測方法得到了不斷發(fā)展。CHAMP（chanllenging mini－satellite payload）、GRACE（gravity recovery and climate experiment）以及GOCE（gravity field and steadystate ocean circulation explorer）［1－3］衛(wèi)星任務(wù)就是此技術(shù)背景下的產(chǎn)物。全球重力位模型可以用球諧系數(shù)進行表達，可以看做是觀測數(shù)據(jù)以面球諧函數(shù)為基函數(shù)在球面上的正交展開。這就需要覆蓋全球的重力場數(shù)據(jù)，即衛(wèi)星以極軌飛行，否則面球諧函數(shù)將不滿足其球面正交的特性。但是由于各種工程技術(shù)問題，目前的衛(wèi)星任務(wù)均是以近極軌飛行，這就帶來了地球兩極各有一個小的數(shù)據(jù)空白區(qū)域，也就是極空白問題。

極空白問題導(dǎo)致了在地球重力位模型的解算中，低次球諧系數(shù)的估計精度較低［4］。目前解決極空白問題的方案，一是利用已有的重力位模型、航空和地面等重力數(shù)據(jù)對極空白區(qū)域進行數(shù)據(jù)填充［5］，使其滿足球諧系數(shù)的正交性，但這會帶來先驗信息的引入；二是利用數(shù)理統(tǒng)計方法降低或消除極空白問題的影響，包括正則化方法等［6－8］；三是對原有基函數(shù)的非正交特性進行數(shù)學(xué)處理或者引入新的正交函數(shù)，從本質(zhì)上解決極空白問題所帶來的影響，包括Cholesky分解、特征值分析、Gram矩陣法等正交處理策略以及Slepian［9］局部譜分析方法等［10－14］。

本文推導(dǎo)利用衛(wèi)星重力梯度數(shù)據(jù)反演地球位模型的Slepian局部譜分析方法表達式，研究Slepian函數(shù)的數(shù)學(xué)特性。在此基礎(chǔ)上，利用衛(wèi)星重力梯度模擬數(shù)據(jù)分析Slepian方法在求解地球重力位模型中的優(yōu)點和不足。

2 Slepian方法的引入與求解

地球重力場量y（r）（r表示坐標(biāo)向量）在頻域上是具有無限帶寬的，但是在進行參數(shù)估計時，僅僅能獲得一定帶寬內(nèi)的估計值，即估計量是具有頻域帶限特征的。構(gòu)建一種新的帶限基函數(shù)g（r），表示為面球諧函數(shù)的線性組合形式

式中，gnm為待求系數(shù)。

地球重力場量的帶限估計值可以寫為

式中，Ω代表整個球面；δ為Kronecker符號。

在利用衛(wèi)星觀測數(shù)據(jù)求解地球重力位模型的過程中，兩極數(shù)據(jù)空白帶來了面球諧函數(shù)的不正交問題。如果新的帶限函數(shù)能最大限度地獲取球面中帶為極空白球冠的余緯半徑，對應(yīng)于GOCE任務(wù)，約為6.7°）上的觀測信息，即在球面上具有局部最優(yōu)特性，則該問題可以得到較為圓滿的解決。為此，令新的帶限函數(shù)的頻譜能量在球帶B上最大，即

此即所謂的Slepian問題［9，12－14］。

在Slepian問題中，比值0＜λ＜1作為頻域帶限函數(shù)g（r）空域集中度的一種度量。聯(lián)合式（1）、式（3）、式（4），可將Slepian問題轉(zhuǎn)化為N＋1個（N－m＋1）×（N－m＋1）維實特征值問題［13－14］

式中，矩陣元素利用式（7）計算［7，15］

式中，ˉFnmk為傾角函數(shù)；E［·］代表取整算子。

對每一個特征值問題，將特征值λ1、λ2、…、λ（N－m＋1）和相應(yīng)的特征向量g1、g2、…、g（N－m＋1）按照1＞λ1≥λ2≥…≥λ（N－m＋1）＞0進行排序，由于不存在能量完全集中于球帶B上的帶限基函數(shù)，故最大特征值λ1非常接近于1但不等于1。同樣，矩陣Dm的正定性保證了最小特征值λ（N－m＋1）非常接近于零但不等于零。相應(yīng)的帶限函數(shù)g1（r）、g2（r）、…、g（N－m＋1）（r）滿足雙正交性

寫成矩陣形式為

Slepian函數(shù)在整個球面和球帶上均是正交的，其本質(zhì)可看做是同次面球諧函數(shù)的線性組合形式，次數(shù)m對于Slepian函數(shù)和面球諧函數(shù)具有相同意義。

3 Slepian函數(shù)的性質(zhì)

Slepian函數(shù)中，在球帶B上能量較大的特征函數(shù)gα（r）所對應(yīng)的特征值λα接近于1，稱之為有效特征值；而在球帶B上能量較小的特征函數(shù)的特征值λα則接近于零，稱之為無效特征值。這些特征值的大小反映了其所對應(yīng)的特征向量的能量大小。與最大特征值相對應(yīng)的基函數(shù)g1（r）是在球帶B上能量最大的帶限函數(shù)，g2（r）則是僅次于g1（r）且與其正交的在球帶B上能量最大的帶限函數(shù)，依此類推。當(dāng)極空白區(qū)域較小時，有效特征值較多，無效特征值相對很少，反之亦然。圖1給出了N＝100，θ0＝6.7°時的特征函數(shù)gα（r）的特征值分布。圖中顯示，由于極空白區(qū)域很小，絕大多數(shù)的特征值均接近于1，僅有小部分的特征值接近于零，而位于這兩種情況中間的過渡帶相對較短，僅有極少數(shù)特征值。

圖1 Slepian函數(shù)的特征值分布ig.1 The distribution of Slepian function eigenvalues

圖2給出了N＝18、極空白半徑θ0＝40°、次數(shù)0≤m≤3時的能量占前3位和后3位的Slepian函數(shù)緯向分布圖（不考慮經(jīng)向變化）。特征值截至小數(shù)點后6位。具有最大權(quán)重（λ≈1）的前3位特征函數(shù)中，能量非常好地集中于40°≤θ≤140°的范圍上，在其他區(qū)域上則基本為零。相反的，具有最小權(quán)重（λ≈0）的前3位特征函數(shù)中，能量非常好地集中于0°≤θ≤40°以及140°≤θ≤180°的范圍上。通過這種分布，Slepian函數(shù)實現(xiàn)了將衛(wèi)星觀測信息最大化地聚集在非極空白區(qū)域，即保證了衛(wèi)星軌道的測量特點與數(shù)據(jù)處理方法的一致性。

圖2 前3位和后3位特征函數(shù)緯向變化圖Fig.2 Co－latitudinal dependence of the first three and last three eigenfunctions

此外，從圖2可以觀察到，當(dāng)α為奇數(shù)時，特征函數(shù)為一以赤道為對稱軸的偶函數(shù)；當(dāng)α為偶數(shù)時，特征函數(shù)為一奇函數(shù)。當(dāng)m為偶數(shù)時，具有最小權(quán)重的特征函數(shù)的奇偶性與N的奇偶性相同；當(dāng)m為奇數(shù)時，具有最小權(quán)重的特征函數(shù)的奇偶性與N的奇偶性相反。也就是說，特征函數(shù)的奇偶性與α－m的奇偶性相反。當(dāng)α＝1時，特征函數(shù)不具有零點；當(dāng)α＝2時，只有一個零點；當(dāng)α＝3時，特征函數(shù)有兩個零點，依此類推。也就是說隨著α的增加，特征函數(shù)的零點個數(shù)逐漸增加。

4 數(shù)值計算效率和穩(wěn)定性的改進

當(dāng)極空白的球冠半徑較小時，式（7）的對角化運算會變得很不穩(wěn)定，特征向量的求解存在不唯一性。另一方面，在求解特征值時，需要通過數(shù)值積分或者其他方法計算每一個矩陣元素，總的運算次數(shù)為o（N2），這大大增加了運算時間。這里引入一個二階差分算子Grünbaum算子［13－14，16－17］

引入一（N－m＋1）×（N－m＋1）矩陣Tm，矩陣元素為

矩陣Dm與Tm均為對稱陣，且滿足DmTm＝TmDm，即二者為可交換矩陣。由于可交換矩陣具有相同的特征向量，故式（5）等價于特征值問題

式中，I為Grünbaum算子；Λ≠λ為相應(yīng)的Grünbaum特征值。

Grünbaum矩陣Tm可以利用下式進行計算得到［13－14］

Tm＝為一三對角陣，這在很大程度上降低了計算的復(fù)雜度，另一方面，Grünbaum算子的運用提高了數(shù)值計算的穩(wěn)定性。由此實現(xiàn)了Slepian問題的穩(wěn)定快速求解。

圖3顯示了采用Grünbaum矩陣進行求解后的Slepian函數(shù)的二維空間分布圖（紅藍分別代表正負值）。從中可以發(fā)現(xiàn)，由于經(jīng)向的正負交叉點與三角函數(shù)sin mλ或者cos mλ的性質(zhì)有關(guān)，因此當(dāng)m＝0時，經(jīng)向上不存在零點；m＝±1時，經(jīng)線上存在兩個零點；m＝±2時，經(jīng)線上存在4個零點，依此類推?？臻g分布圖的緯向分布特征與圖2的分析一致。

圖3 前4位和后4位特征函數(shù)的空間分布圖（90°＜θ＜180°）Fig.3 Spatial dependence of the first four and last four eigenfunctions（90°＜θ＜180°）

圖4為勒讓德函數(shù)和Slepian函數(shù)緯向空間分布圖（不考慮經(jīng)向變化），最高階數(shù)為100，次數(shù)固定為m＝0。圖中可以看出，與勒讓德函數(shù)的頻譜分布較為均勻不同，Slepian函數(shù)的頻譜能量隨著緯度增加逐漸減小，在兩極地區(qū)達到最小，這種數(shù)學(xué)特性與衛(wèi)星大地測量中的極空白問題符合很好。這也是引入Slepian函數(shù)的主要初衷。而且，頻譜能量與階數(shù)密切相關(guān)，階數(shù)越高，能量逐漸變小甚至消失。

圖4 勒讓德函數(shù)與Slepian函數(shù)的緯向頻譜分布Fig.4 Co－latitudinal dependence of Legendre function and Slepian function

5 Slepian方法在衛(wèi)星重力梯度數(shù)據(jù)恢復(fù)地球重力場中的應(yīng)用研究

5.1 引力梯度分量的Slepian函數(shù)表達

由式（11）可知，Gm為解特征值問題（式（5））所得的特征向量矩陣，為一正交矩陣，故

則對于平方可積空間L2（σ）內(nèi)相應(yīng)的線性延拓算子有

可見，Slepian函數(shù)的延拓矩陣為非對角陣形式。

將引力位在地心球坐標(biāo)系（r，θ，λ）下利用Slepian函數(shù)進行展開得到（截斷至N階）

根據(jù)引力梯度分量在地心球坐標(biāo)（r，θ，λ）與局部指北坐標(biāo)（x，y，z）之間的轉(zhuǎn)換關(guān)系［18］，可得到引力梯度分量的Slepian函數(shù)表達式為（截斷至N階）

根據(jù)式（23）在最小二乘框架下求解即可得到Slepian系數(shù)

式中，P為權(quán)陣；A為相應(yīng)的設(shè)計矩陣；l為梯度觀測值向量。

表1 局部指北坐標(biāo)系下的引力梯度分量Slepian函數(shù)表達式Tab.1 Slepian expression of gravity gradient exponents in local north－oriented coordinates

5.2 衛(wèi)星重力梯度數(shù)據(jù)恢復(fù)地球重力場的Slepian方法仿真分析

基于德國波恩大學(xué)提供的GOCE衛(wèi)星軌道仿真數(shù)據(jù)，利用EIGEN－GL04C位模型的前300階模擬了30d、5s采樣的衛(wèi)星重力梯度觀測值。利用梯度徑向分量Vzz，基于空域最小二乘方法［8，19］進行地球重力位模型的解算。為使得模擬方案更符合GOCE重力位模型實際解算情況，階數(shù)僅取至200階。圖5給出了以面球諧函數(shù)為基函數(shù)的球諧系數(shù)解的相對誤差以及以Slepian函數(shù)為基函數(shù)的Slepian解的相對誤差。

圖5 球諧解的相對誤差與Slepian解的相對誤差比較Fig.5 The relative errors of spherical harmonic solution and Slepian solution

從圖5可以看出，由于混頻效應(yīng)的影響，球諧系數(shù)解的低次項的估計精度較低，此外高階項的求解精度亦受到很大影響。將求解過程轉(zhuǎn)換至Slepian域后，求解結(jié)果明顯變優(yōu)。低次項的估計精度有了大幅度的提高，高階項的求解精度也有所提高，但是扇諧項及其周圍系數(shù)的精度要低一些。對于球諧系數(shù)解來講，混頻效應(yīng)主要包括兩部分：第1部分是200階以上的重力場信號混疊進200階以下的重力場信號引起；第2部分是地球重力場非極空白區(qū)域的觀測信息混疊進極空白區(qū)域引起的重力場頻譜的變化，由于低次勒讓德函數(shù)在高緯度地區(qū)包含更多的能量信號，故而低次項球諧系數(shù)的估計精度較低［4］。就Slepian系數(shù)解而言，混頻效應(yīng)主要是高階項（大于200階）的重力場信號混疊進低階項所致，且主要影響Slepian函數(shù)的扇諧項及其周圍系數(shù)。

圖6顯示了球諧系數(shù)解和Slepian解的空間域結(jié)果。為了便于量級上的比較，這里對低次項進行了截斷。圖中可以發(fā)現(xiàn)，Slepian解較之球諧系數(shù)解的大地水準(zhǔn)面誤差均方差并沒有明顯的改善。這說明，Slepian解對于大地水準(zhǔn)面的求解并沒有直接的幫助，大地水準(zhǔn)面精度的提高還是需要依賴于高精度的海量觀測值的獲取。當(dāng)存在噪聲時，需要借助數(shù)理統(tǒng)計方法，特別是正則化方法進行處理。Slepian方法的主要貢獻在于其開辟了Slepian頻域的新范疇，且模型在該頻域內(nèi)的求解精度較之傳統(tǒng)的球諧域要高。如何結(jié)合數(shù)理統(tǒng)計方法在Slepian域內(nèi)對求解結(jié)果進行有效處理，以提高大地水準(zhǔn)面的求解精度，有待進一步的研究和分析。

圖6 球諧解與Slepian解的大地水準(zhǔn)面誤差緯向分布圖Fig.6 Co－latitudinal dependence of geoid error RMS of spherical harmonic solution and Slepian solution

6 結(jié) 論

利用衛(wèi)星重力梯度數(shù)據(jù)恢復(fù)地球重力場的過程中，極空白問題使得球諧系數(shù)的低次項估計精度較低。Slepian函數(shù)能夠在很大程度上解決面球諧函數(shù)的不正交問題，實現(xiàn)對地球重力場的有效恢復(fù)。本文對Slepian函數(shù)的數(shù)學(xué)性質(zhì)進行詳細分析，改進了求解的穩(wěn)定性和求解效率，針對衛(wèi)星重力梯度數(shù)據(jù)反演地球重力場，推導(dǎo)了引力梯度分量的Slepian函數(shù)實用展開式，最后通過仿真分析，就Slepian局部譜分析方法在衛(wèi)星重力梯度數(shù)據(jù)確定地球重力位模型中的應(yīng)用和前景進行了分析和討論。研究表明：

Slepian函數(shù)在整個球面和球帶上具有雙正交性，可看做是同次面球諧函數(shù)的線性組合形式，其奇偶性與α－m的奇偶性相反。Slepian頻譜能量隨著緯度的增加逐漸減小，在兩極地區(qū)達到最小，實現(xiàn)了將衛(wèi)星觀測信息最大化地聚集在非極空白區(qū)域，即保證了衛(wèi)星軌道的測量特點與數(shù)據(jù)處理方法的一致性。Grünbaum算子的引入有效提高了Slepian方法求解的穩(wěn)定性和效率。

Slepian方法的優(yōu)良特性體現(xiàn)在Slepian頻域內(nèi)的地球重力場求解精度較之球諧域內(nèi)的結(jié)果要高，但其對大地水準(zhǔn)面空間分布的恢復(fù)精度貢獻不明顯。地球重力場和大地水準(zhǔn)面精度的改進依賴于觀測數(shù)據(jù)精度的提高以及觀測數(shù)據(jù)頻譜結(jié)構(gòu)的改善。

［1］ REIGBER C，LTIHR H，SCHWINTZER P.CHAMP Mission Status［J］.Advances in Space Research，2002，30（2）：129－134.

［2］ TAPLEY B D，BETTADPUR S.The Gravity Recovery and Climate Experiment：Mission Overview and Early Results［J］.Geophysical Research Letters，2004，31（9）：1－6.

［3］ ESA.Gravity Field and Steady－state Ocean Circulation Mission［M］.Paris：ESA Publications Division，1999.

［4］ SNEEUW N J，GELDEREN M V.The Polar Gap［J］.Lecture Notes in Earth Sciences，1997，65：559－568.

［5］ TSCHERNING C C，F(xiàn)ORSBERG R，ALBERTELLA A，et al.The Polar Gap Problem：Space－wise Approaches to Gravity Field Determination in Polar Areas，from E?tv?s to mGal［R］.Copenhagen：University of Copenhagen，2000，331－336.

［6］ DITMAR P，KUSCHE J，KLEES R.Computation of Spherical Harmonic Coefficients from Gravity Gradiometry Data to Be Acquired by the GOCE Satellite：Regularization Issues［J］.Journal of Geodesy，2003，77（7－8）：465－477.

［7］ METZLER B.PAIL R.GOCE Data Processing：The Spherical Cap Regularization Approach［J］.Studia Geophysica et Geodaetica，2005，49（4）：441－462.

［8］ 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.）

［9］ SLEPIAN D.Some Comments on Fourier－analysis，Uncertainty and Modeling［J］.SIAM Rev，1983，25（3）：379－393.

［10］ ALBERTELLA A，SANSòF，SNEEUW N.Band－limited Functions on a Bounded Spherical Domain：The Slepian Problem on the Sphere［J］.Journal of Geodesy，1999，73（9）：436－447.

［11］ PAIL R，PLANKG，SCHUH W D.Spatially Restricted Data Distribution on the Sphere：The Method of Orthonormalized Functions and Applications［J］.Journal of Geodesy，2001，75（1）：44－56.

［12］ MARK A W，SIMONS F J.Localized Spectral Analysis on the Sphere［J］.Geophysical Journal International，2005，162（3）：655－675.

［13］ SIMONS F J，DAHLEN F A.Spherical Slepian Functions and the Polar Gap in Geodesy［J］.Geophysical Journal International，2006，166（3）：1039－1061.

［14］ SIMONS F J，DAHLEN F A，WIECZOREK M A.Spatiospectral Concentration on a Sphere［J］.SIAM Review，2006，48（3）：504－536.

［15］ SNEEUW N.Inclination Functions：Group Theoretical Background and a Recursive Algorithm［D］.Delft：Delft University of Technology，1991.

［16］ GILBER E N，SLEPIAN D.Doubly Orthogonal Concentrated Polynomials［J］.SIAM Journal on Mathematical Analysis，1977，8（2）：290－319.

［17］ GRüNBAUM F A，LONGHI L，PERLSTADT M.Differential Operators Commuting with Finite Convolution Integral Operators：Some Non－Abelian Examples［J］.SIAM Journal on Applied Mathematics，1982，42（5）：941－955.

［18］ KOOP R.Global Gravity Field Modeling Using Satellite Gravity Gradiometry［D］.Delft：Neth Geod Comm，1993.

［19］ XU Xinyu.Study of Determining the Earth’s Gravity Field from Satellite Gravity Gradient and Satellite－to－Satellite Tracking Data［D］.Wuhan：Wuhan University，2008.（徐新禹.衛(wèi)星重力梯度及衛(wèi)星跟蹤衛(wèi)星數(shù)據(jù)確定地球重力場的研究［D］.武漢：武漢大學(xué)，2008.）

Slepian Localized Spectral Analysis of the Determination of the Earth’s Gravity Field Using Satellite Gravity Gradiometry Data

ZHU Guangbin1，LI Jiancheng2，WEN Hanjiang3，CHANG Xiaotao1，WANG Zhengtao2，ZOU Xiancai2
1.Satellite Surveying and Mapping Application Center，National Administration of Surveying，Mapping and Geoinformation，Beijing 100830，China；2.School of Geodesy and Geomatics，Wuhan University，Wuhan430079，China；3.Chinese Academy of Surveying and Mapping，Beijing100830，China

The Slepian localized spectral analysis method is introduced and its mathematical characteristics are analyzed.Grünbaum operator is introduced to improve the computational efficiency and stability.Then，the formulas of Slepian method for recovering the gravity field with satellite gravity gradiometry data are educed and numerical analysis are done.The results show that，Slepian function is orthogonal in both the sphere and the belt，and its spectral energy distribution characteristic is in accord with the orbit characteristic.Also，the accuracy of low order coefficients in Slepian domain is less influenced by polar gaps than that in spherical harmonic domain.However，the direct contribution of Slepian method to the improvement of geoid is limited.

satellite gravity gradiometry；gravity field；Slepian localized spectral analysis；polar gaps；surface spherical harmonics；Slepian function

ZHU Guangbin（1981—），male，PhD，majors in satellite geodesy.

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.）

P223

A

1001－1595（2012）01－0001－07

國家863計劃（2007AA12Z346）；國家自然科學(xué)基金（40874012；40974016；41004007）

叢樹平）

2011－01－13

2011－05－31

朱廣彬（1981—），男，博士，研究方向為衛(wèi)星大地測量。

E－mail：whu＿gbzhu＠hotmail.com