王 燚，姜效典

中國海洋大學(xué)海洋地球科學(xué)學(xué)院，山東 青島 266100

擾動重力梯度的球冠諧分析建模

王 燚，姜效典

中國海洋大學(xué)海洋地球科學(xué)學(xué)院，山東 青島 266100

從球冠諧理論出發(fā)，詳細(xì)推導(dǎo)了球冠坐標(biāo)系下擾動重力梯度的無奇異性計算公式?；赥ikhonov正則化方法，利用GOCE衛(wèi)星實際觀測數(shù)據(jù)解算局部重力場球冠諧模型。數(shù)值計算表明，基于擾動重力梯度的球冠諧分析建模方法能夠有效地恢復(fù)局部重力場中的短波信號，與GO_CONS_GCF_2_DIR_R5模型的差異在±0.3×10-5m/s2水平。

球冠諧和分析; 擾動重力; 重力梯度; Tikhonov正則化方法

隨著CHAMP、GRACE和GOCE等3顆重力衛(wèi)星的發(fā)射升空，利用衛(wèi)星重力觀測值恢復(fù)地球重力場的研究，已經(jīng)成為當(dāng)前地球物理學(xué)研究的熱點問題之一。文獻(xiàn)[1—2]的研究表明，重力梯度張量反映了重力位水準(zhǔn)面的曲率和力線彎曲，對重力場中的短波信號變化敏感，所以重力梯度張量更能夠反映重力場的精細(xì)結(jié)構(gòu)，包含有大量的局部重力場的有效信息。

國內(nèi)外眾多學(xué)者就利用衛(wèi)星重力觀測值恢復(fù)地球重力場的方法做了大量的研究工作。文獻(xiàn)[3]提出了廣義輪胎調(diào)和分析方法，利用衛(wèi)星重力梯度張量分量組合以較高的精度還原地球重力場。文獻(xiàn)[4]提出了使用最小二乘法恢復(fù)地球重力場的方法，通過數(shù)值模擬，證明重力異常的梯度對軌道高度的變化不敏感。在文獻(xiàn)[5]中研究了全張量重力梯度數(shù)據(jù)的全局和局部分量的廣義球諧譜表示和軌道根數(shù)表示，給出了廣義球諧函數(shù)與球諧函數(shù)之間的關(guān)系，從理論上得到了全張量重力梯度數(shù)據(jù)的描述方法和由全張量重力梯度網(wǎng)格數(shù)據(jù)恢復(fù)全球重力位譜系數(shù)的基本公式。然而這些研究工作在理論方面都是利用傳統(tǒng)的球諧分析法。在表達(dá)全球重力場時，球諧分析法有明顯的優(yōu)勢，但是對于局部區(qū)域重力場的精細(xì)結(jié)構(gòu)和重力場高頻信息的表述上，球諧分析法存在著不足和缺陷。例如球諧函數(shù)不能滿足局部重力場的邊界條件，在表述高分辨率重力場時，需要求解巨量的球諧系數(shù)等[6]。因此需要針對衛(wèi)星重力梯度恢復(fù)局部重力場的理論方法進(jìn)行系統(tǒng)的研究。

文獻(xiàn)[7]在前人的研究基礎(chǔ)上總結(jié)提出了球冠諧分析理論，對地球北極區(qū)域的磁異常做了數(shù)值分析。文獻(xiàn)[8]利用航磁異常觀測值，解算了中國大陸航磁異常的不同階次球冠諧模型，結(jié)合大地構(gòu)造理論對相關(guān)圖件進(jìn)行了分析解釋。文獻(xiàn)[9]研究了球冠諧函數(shù)在重力場中的表達(dá)式，并利用球冠諧理論建立了我國的局部重力場模型。文獻(xiàn)[10]提出了一種基于點質(zhì)量模型的多層球冠諧系數(shù)迭代構(gòu)造算法，并對求解矩陣A的病態(tài)性和球冠半角與點質(zhì)量埋藏深度和密度進(jìn)行了詳細(xì)的分析研究，通過數(shù)值模擬表明當(dāng)?shù)惴ㄊ諗繒r，球冠諧模型能夠以高精度擬合重力觀測值。文獻(xiàn)[11]提出了使用球冠諧模型的重力異常徑向?qū)?shù)和曲面導(dǎo)數(shù)建立曲面位場的方法。這些研究工作都取得了令人滿意的結(jié)果。然而利用擾動重力梯度恢復(fù)局部球冠諧重力場的研究，到目前為止還不多見。由此，本文從球冠諧理論出發(fā)，推導(dǎo)了無δ奇異性球冠諧梯度計算公式，使用GOCE衛(wèi)星數(shù)據(jù)產(chǎn)品中地球固定參考框架(terrestrial reference frame，TRF)下的重力梯度數(shù)據(jù)恢復(fù)了4個球冠諧重力場模型，最后在計算點上，與文獻(xiàn)[12]中使用GOCE衛(wèi)星梯度數(shù)據(jù)建立的GO_CONS_GCF_2_DIR_R5地球重力場模型的擾動重力數(shù)值進(jìn)行了比對分析，檢驗該方法的正確性。

1 球冠諧和分析擾動位重力梯度模型

球冠諧和分析是一種在球冠坐標(biāo)下的滿足位理論邊值條件的譜函數(shù)法，最早由加拿大人Haines在研究與地殼相關(guān)的局部區(qū)域磁異常時提出。通過利用非整階勒讓德函數(shù)代替整階勒讓德函數(shù)，余緯θ的定義域由[0,π]變換到[0，α](α是球冠半角)，同時將地球Z軸由北極點沿子午線旋轉(zhuǎn)到球冠中心，以保證球冠上點的經(jīng)度依然滿足周期性的邊界條件。采用分離變量法可以得到重力場擾動位在球冠坐標(biāo)下的解[13]

(1)

(2)

(3)

系數(shù)Aj(m,l)的遞推關(guān)系如下

(4)

定義

(5)

則由文獻(xiàn)[7]推導(dǎo)的近似式可得

(6)

由上述推導(dǎo)過程可知

(7)

如果令

為球冠諧系數(shù)，則式(1)可以改寫為

(8)

式中，G是萬有引力常數(shù)；M是地球質(zhì)量。因為δ與λc是球冠坐標(biāo)下的余緯和經(jīng)度(參見圖2)，所以將球冠中心Q看作是新的北極點，定義球冠坐標(biāo)系下的球冠局部指北坐標(biāo)系(spherical cap local north-oriented reference frame，SCLNOF)的x軸指向球冠極點，z軸平行地心向徑，方向向外，y軸與x和z軸構(gòu)成右手坐標(biāo)系，則與球諧函數(shù)的局部指北坐標(biāo)系(local north-oriented reference frame，LNOF )下重力梯度張量計算一樣，將式(8)代入局部指北坐標(biāo)系梯度計算公式[16]，可以得到球冠局部指北坐標(biāo)下的擾動重力和擾動重力梯度張量模型公式

(9)

(10)

式中

(11)

從式(9)、(10)中可以看出，在球冠極點附近，當(dāng)sinδ趨于零時，存在δ奇異性，因此根據(jù)非整階勒讓德函數(shù)的計算公式，本文推導(dǎo)了如下的去δ奇異性球冠諧重力梯度張量計算公式。記

(12)

(13)

(14)

(15)

將式(15)代入式(10)，就可以得到非奇異的球冠局部指北坐標(biāo)下擾動重力梯度張量的計算公式。根據(jù)球冠諧分析理論，整階勒讓德函數(shù)是非整階勒讓德函數(shù)的特例[17-18]，因此當(dāng)球冠半角α=π/2時，式(15)中擾動位梯度張量的各個階數(shù)l為整數(shù)，這時非整階勒讓德函數(shù)的計算結(jié)果與整階勒讓德函數(shù)的計算結(jié)果應(yīng)當(dāng)一致。下面取α=π/2，δ∈[0,π/2]，通過式(15)計算各個梯度中l(wèi)=2,3,m=0,1,2,3對應(yīng)的非整階勒讓德函數(shù)項，并與文獻(xiàn)[19]中的擾動重力梯度張量公式中n=2,3,m=0,1,2,3整階勒讓德函數(shù)項計算結(jié)果進(jìn)行比較，統(tǒng)計結(jié)果如表1所示，二者差值的均方根小于10-3數(shù)量級。從表1可見，當(dāng)m≠0且m≠n時，二者誤差較大，這是由于非整階勒讓德函數(shù)值使用近似式(6)計算，而當(dāng)m=0或m=n時，式(6)不可用，必須用式(5)計算，因此誤差范圍不一致。圖1是l=n=3,m=2時非整階勒讓德函數(shù)和整階勒讓德函數(shù)擾動重力梯度張量誤差。

表1 整階勒讓德函數(shù)和非整階勒讓德函數(shù)的誤差均方根Tab.1 RMS errors between Legendre function and Legendre function with non-integral degree

圖1 擾動重力梯度誤差(l=n=3、m=2)Fig.1 Error of disturbing gravity gradients(l=n=3、m=2)

2 不同坐標(biāo)系下重力梯度張量之間的轉(zhuǎn)換

如圖2所示，GOCE衛(wèi)星任務(wù)Level 2產(chǎn)品的EGG_TRF_2中重力梯度數(shù)據(jù)的軌道位置使用地心地固坐標(biāo)系(θ,λ,r)表示，其中θ表示地心余緯，λ表示地心經(jīng)度，r是衛(wèi)星軌道距地心距離；重力梯度張量使用的是局部指北坐標(biāo)系(LNOF)，因此使用GOCE衛(wèi)星重力梯度恢復(fù)球冠諧重力場模型時，需要進(jìn)行坐標(biāo)變換。首先是將軌道位置坐標(biāo)(θ,λ)轉(zhuǎn)換到球冠坐標(biāo)系[20]。定義計算點P在球坐標(biāo)下的位置用(θ,λ)表示，球冠坐標(biāo)系下的余緯和經(jīng)度用δ和λc表示，設(shè)球冠中心點Q在球坐標(biāo)系中為(θ0,λ0)，則根據(jù)球面三角公式可知，計算點P在球冠坐標(biāo)系中的余緯δ和λc經(jīng)度計算公式為

(16)

然后將重力梯度張量從局部指北坐標(biāo)系(LNOF)轉(zhuǎn)換到球冠指北坐標(biāo)系(SCLNOF)中。定義LNOF坐標(biāo)系到SCLNOF坐標(biāo)系的轉(zhuǎn)換矩陣。設(shè)β為球冠上點P到球冠中心點Q的方位角,如圖2。顯然SCLNOF坐標(biāo)系就是將點P處的LNOF坐標(biāo)系繞Z軸旋轉(zhuǎn)一個β角度得到，所以LNOF坐標(biāo)系到SCLNOF坐標(biāo)系的梯度張量轉(zhuǎn)換矩陣為

(17)

式中，方位角β的cosβ與sinβ的計算公式為

(18)

式中，δ為點P在球冠坐標(biāo)下的余緯。

圖2 球冠坐標(biāo)系Fig.2 Spherical cap coordinate system

3 球冠諧梯度公式恢復(fù)局部重力場的數(shù)值模擬計算和分析

(19)

(20)

表2 擾動重力梯度Txz分量統(tǒng)計結(jié)果Tab.2 Disturbing gravity gradient Txz survey results E

表3 擾動重力統(tǒng)計結(jié)果Tab.3 Disturbing gravity survey results 10-5 m/s2

圖3 球冠諧模型對應(yīng)的L曲線Fig.3 L-curve according to spherical cap harmonic models

從表6可以看出，4個球冠諧模型的Laplace方程的精度都在±10-8E，由于GOCE衛(wèi)星任務(wù)的重力梯度張量對角線分量的誤差PSD在測量帶寬范圍內(nèi)不大于±3.2×10-3E/Hz1/2(見文獻(xiàn)[16])，因此文中計算的球冠諧模型完全滿足衛(wèi)星重力梯度測量的精度要求，可以用于恢復(fù)地球局部重力場模型。

表4 球冠諧模型擾動重力梯度的誤差Tab.4 Disturbing gravity gradient error of the spherical cap harmonic models E

表5 球冠諧模型徑向擾動重力的誤差Tab.5 Disturbing gravity error of the spherical cap harmonic models 10-5 m/s2

表6 球冠諧模型計算的擾動重力梯度張量對角線分量所滿足Laplace方程的精度Tab.6 The diagonal components of the disturbing gravity gradients computed by the spherical cap harmonic models that fulfill the precision of Laplace equation E

圖4 球冠諧模型徑向擾動重力誤差Fig.4 Error of the spherical cap harmonic models

4 結(jié) 論

本文從球冠諧分析理論出發(fā)，詳細(xì)推導(dǎo)了球冠坐標(biāo)系下的擾動位重力梯度計算公式，將GOCE衛(wèi)星任務(wù)重力梯度觀測值延拓到平均軌道面上。網(wǎng)格化后，使用Tikhonov正則化算法反演恢復(fù)了4個球冠諧重力場模型，最后與使用GO_CONS_GCF_2_DIR_R5重力場模型計算的衛(wèi)星平均軌道面擾動重力進(jìn)行比較分析。計算結(jié)果表明，由于擾動重力梯度中各個分量對重力場頻譜的敏感性是不一樣的，因此單一徑向擾動重力梯度分量不能恢復(fù)全部擾動重力場的信號，所以聯(lián)合不同的重力梯度分量，利用其頻譜的互補性恢復(fù)的局部重力場模型，能夠包含更多的重力場信號。如果從擾動重力梯度數(shù)據(jù)中移除重力場的長波信號，這時無論是使用單一徑向梯度分量，還是使用多個梯度分量聯(lián)合的方式恢復(fù)的球冠諧模型，都能夠很好地擬合擾動重力場中剩余的重力信號，與GO_CONS_GCF_2_DIR_R5模型的差異不大于±0.3x10-5m/s2。

[1] 吳曉平.利用衛(wèi)星重力梯度與地面數(shù)據(jù)確定地球重力場的分析[J].測繪學(xué)院學(xué)報，1991(1)：19-32.

WU Xiaoping.Some Analyses of the Determination of the Earth’s Gravity Field from the Satellite Gravity Gradient and Terrestrial Data[J].Journal of the Institute of Surveying and Mapping，1991(1)：19-32.

[2] 孟嘉春，蔡喜楣.衛(wèi)星重力梯度測量及其應(yīng)用前景探討[J].地球物理學(xué)報，1991，34(3)：369-376.

MENG Jiachun，CAI Ximei.Approach on Satellite Gravity Gradiometry and Its Vistas of Applications[J].Acta Geophysica Sinica，1991，34(3)：369-376.

[3] 吳星，張傳定，趙東明.衛(wèi)星重力梯度分量的廣義輪胎調(diào)和分析[J].測繪學(xué)報，2009，38(2)：101-107.

WU Xing，ZHANG Chuanding，ZHAO Dongming.Generalized Torus Harmonic Analysis of Satellite Gravity Gradients Component[J].Acta Geodaetica et Cartographica Sinica，2009，38(2)：101-107.

[4] REED G B.Application of Kinematical Geodesy for Determining the Short Wave Length Components of the Gravity Field by Satellite Gradiometry[R].Columbus：Ohio State University，1973.

[5] 張傳定，吳曉平，陸仲連.全張量重力梯度數(shù)據(jù)的譜表示方法[J].測繪學(xué)報，2000，29(4)：297-304.

ZHANG Chuanding，WU Xiaoping，LU Zhonglian.Spectral Representation of the Full Gravity Tensor[J].Acta Geodaetica et Cartographica Sinica，2000，29(4)：297-304.

[6] 王燚，姜效典，李德勇.多源重力數(shù)據(jù)球冠諧模型抗差融合法[J].測繪學(xué)報，2015，44(9)：952-957.DOI：10.11947/j.AGCS.2015.20140345.

WANG Yi，JIANG Xiaodian，LI Deyong.The Robust Fusion of Multi-source Gravity Data Based on the Spherical Cap Harmonic Model[J].Acta Geodaetica et Cartographica Sinica，2015，44(9)：952-957.DOI：10.11947/j.AGCS.2015.20140345.

[7] HAINES G V.Spherical Cap Harmonic Analysis[J].Journal of Geophysical Research：Solid Earth，1985，90(B3)：2583-2591.

[8] 申寧華，王光杰，王喜臣.球冠諧和分析算法及其在中國大陸航磁異常分析中的應(yīng)用[J].物探化探計算技術(shù)，1998，20(1)：9-18.

SHEN Ninghua，WANG Guangjie，WANG Xichen.Spherical Cap Harmonic Analysis and Its Application to the Analysis of the Aeromagnetic Anomalies in China[J].Computing Techniques for Geophysical and Geochemical Exploration，1998，20(1)：9-18.

[9] LI Jiancheng，CHAO Dingbo,NING Jinsheng.Spherical Cap Harmonic Expansion for Local Gravity Field Representation[J].Manuscripta Geodaetica，1995，20(2)：265-277.

[10] 王燚，姜效典.地球重力場球冠諧模型的分層構(gòu)造和分析[J].大地測量與地球動力學(xué)，2014，34(5)：30-34，39.

WANG Yi，JIANG Xiaodian.Structure and Analysis of Multilayered Spherical Cap Harmonic Model[J].Journal of Geodesy and Geodynamics，2014，34(5)：30-34，39.

[11] 吳招才，劉天佑，高金耀.局部重力場球冠諧分析中的導(dǎo)數(shù)計算及應(yīng)用[J].海洋與湖沼，2006，37(6)：488-492.

WU Zhaocai，LIU Tianyou，GAO Jinyao.Derivative Calculation and Application in Spherical Cap Harmonic Analysis of Local Gravity Field[J].Oceanologia et Limnologia Sinica，2006，37(6)：488-492.

[12] BRUINSMA S L，F(xiàn)?RSTE C，ABRIKOSOV O，et al.The New ESA Satellite-only Gravity Field Model via the Direct Approach[J].Geophysical Research Letters，2013，40(14)：3607-3612.

[13] 曹月玲，王解先.利用球冠諧分析擬合GPS水準(zhǔn)高程[J].武漢大學(xué)學(xué)報(信息科學(xué)版)，2008，33(7)：740-743.

CAO Yueling，WANG Jiexian.Application of Spherical Cap Harmonic Analysis to Fit GPS Level Height[J].Geomatics and Information Science of Wuhan University，2008，33(7)：740-743.

[14] 彭富清，于錦海.球冠諧分析中非整階Legendre函數(shù)的性質(zhì)及其計算[J].測繪學(xué)報，2000，29(3)：204-208.

PENG Fuqing，YU Jinhai.The Characters and Computation of Legendre Function with Non-integral Degree in SCHA[J].Acta Geodaetica et Cartographica Sinica，2000，29(3)：204-208.

[15] HWANG C，CHEN S K.Fully Normalized Spherical Cap Harmonics：Application to the Analysis of Sea-level Data from TOPEX/Poseidon and ERS-1[J].Geophysical Journal International，1997，129(2)：450-460.

[16] 劉曉剛.GOCE衛(wèi)星測量恢復(fù)地球重力場模型的理論與方法[D].鄭州：信息工程大學(xué)，2011.

LIU Xiaogang.Theory and Methods of the Earth’s Gravity Field Model Recovery from GOCE Data[D].Zhengzhou：Information Engineering University，2011.

[17] HOBSON E W.The Theory of Spherical and Ellipsoidal Harmonics[M].Cambridge：Cambridge University Press，1931.

[18] 王竹溪，郭敦仁.特殊函數(shù)概論[M].北京：北京大學(xué)出版社，1989.

WANG Zhuxi，GUO Dunren.Introduction to Special Function[M].Beijing：Peking University Press，1989.

[19] 劉曉剛，吳曉平，趙東明，等.擾動重力梯度的非奇異表示[J].測繪學(xué)報，2010，39(5)：450-457.

LIU Xiaogang，WU Xiaoping，ZHAO Dongming，et al.Non-singular Expression of the Disturbing Gravity Gradients[J].Acta Geodaetica et Cartographica Sinica，2010，39(5)：450-457.

[20] 趙建虎，王勝平，劉輝，等.海洋局域地磁場球冠諧分析建模方法研究[J].測繪科學(xué)，2010，35(1)：50-52，9.

ZHAO Jianhu，WANG Shengping，LIU Hui，et al.Study on Establishing Local Geomagnetic Model Using Spherical Cap Harmonic Analysis[J].Science of Surveying and Mapping，2010，35(1)：50-52，9.

[21] 羅志才，鐘波，寧津生，等.衛(wèi)星重力梯度測量確定地球重力場的理論與方法[M].武漢：武漢大學(xué)出版社，2015.

LUO Zhicai，ZHONG Bo，NING Jinsheng，et al.Theory and Method for Determining the Earth’s Gravity Field from Satellite Gravity Gradiometry[M].Wuhan：Wuhan University Press，2015.

[22] 韋博成，魯國斌，史建清.統(tǒng)計診斷引論[M].南京：東南大學(xué)出版社，1991.

WEI Bocheng，LU Guobin，SHI Jianqing.Introduction to Statistical Diagnostics[M].Nanjing：Publishing House of Dongnan University，1991.

[23] 黃謨濤，歐陽永忠，劉敏，等.融合海域多源重力數(shù)據(jù)的正則化點質(zhì)量方法[J].武漢大學(xué)學(xué)報(信息科學(xué)版)，2015，40(2)：170-175.

HUANG Motao，OUYANG Yongzhong，LIU Min，et al.Regularization of Point-mass Model for Multi-source Gravity Data Fusion Processing[J].Geomatics and Information Science of Wuhan University，2015，40(2)：170-175.

[24] HANSEN P C.Analysis of Discrete Ill-posed Problems by Means of the L-curve[J].SIAM Review，1992，34(4)：561-580.

The Spherical Cap Harmonic Analysis Modeling Method Based on Disturbing Gravity Gradients

WANG Yi，JIANG Xiaodian

College of Marine Geo-science, Ocean University of China, Qingdao 266100，China

It is deduced that the non-singular computational formulae of the disturbing gravity gradients based on spherical cap harmonic analysis theory.On the basis of the Tikhonov regularization method，the spherical cap harmonic model of local gravity field is calculated using GOCE satellite real surveying data.It is shown that the short wavelength information of local gravity field can be recovered from the disturbing gravity gradient data using spherical cap harmonic analysis method，and the difference with the GO_CONS_GCF_2_DIR_R5 model is at ±0.3×10-5m/s2.

spherical cap harmonic analysis; disturbing gravity; gravity gradients;Tikhonov regulariza-tion method

The Oceanic Public Nonprofit Research Project(No.201305029-02)

WANG Yi(1974—)，male,PhD，lecturer，majors in geophysical exploration.

王燚，姜效典.擾動重力梯度的球冠諧分析建模[J].測繪學(xué)報，2017,46(11)：1802-1811.

10.11947/j.AGCS.2017.20160412.

WANG Yi，JIANG Xiaodian.The Spherical Cap Harmonic Analysis Modeling Method Based on Disturbing Gravity Gradients[J].Acta Geodaetica et Cartographica Sinica,2017,46(11):1802-1811.DOI:10.11947/j.AGCS.2017.20160412.

P223

A

1001-1595(2017)11-1802-10

海洋公益性行業(yè)科研專項經(jīng)費(201305029-02)

(責(zé)任編輯：叢樹平)

2016-08-18

修回日期：2017-07-26

王燚(1974—)，男，博士，講師，研究方向為地球物理勘探方法與信息技術(shù)。

E-mail：wyfox009@163.com