地球重力場球冠諧模型的分層構(gòu)造和分析*

王 燚 姜效典

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

地球重力場球冠諧模型的分層構(gòu)造和分析*

王 燚 姜效典

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

根據(jù)球冠諧系數(shù)與點(diǎn)質(zhì)量模型的關(guān)系，提出一種基于多層點(diǎn)質(zhì)量模型分層構(gòu)造球冠諧系數(shù)的方法。以32°N～34°N，102°E～104°E為計(jì)算區(qū)域，利用EGM2008模型和實(shí)測觀測值構(gòu)造三層球冠諧模型系數(shù)，逼近該區(qū)域的重力異常場。結(jié)果表明，使用本方法構(gòu)造的球冠諧模型和實(shí)測觀測值的誤差平均值小于0.5×10－5ms－2，當(dāng)擬合區(qū)域的球冠半角為 0.71°時，模型的內(nèi)符合精度為 ±4.65 ×10－5ms－2。

球冠諧和分析;EGM2008模型;點(diǎn)質(zhì)量模型;迭代算法;病態(tài)性

EGM2008模型通過球諧函數(shù)可以較好地描述全球重力場，但在局部重力場描述上，球諧函數(shù)不能滿足局部重力場的邊界條件，對局部重力場的描述只是理論上的一種近似［1］。Haines等［2］指出了運(yùn)用球冠諧展開逼近局部重力場的可能性。Hwang［3］使用球冠諧理論針對海面地形進(jìn)行研究;張勤等［4］推導(dǎo)了球冠諧分析法逼近區(qū)域地殼垂直形變場的模型;李建成［5］運(yùn)用球冠諧理論建立我國的局部重力場模型，都取得比較滿意的結(jié)果。然而在球冠諧系數(shù)計(jì)算過程中，當(dāng)觀測數(shù)小于待求參數(shù)時，基函數(shù)組成的求解球冠諧系數(shù)的矩陣是秩虧的，當(dāng)觀測數(shù)遠(yuǎn)遠(yuǎn)多于待求的球冠諧系數(shù)時，矩陣的條件數(shù)比較大，求逆運(yùn)算的結(jié)果不穩(wěn)定。為了克服該問題，本文引入點(diǎn)質(zhì)量模型，利用逐級余差思想，分層計(jì)算殘差重力異常對應(yīng)的球冠諧系數(shù)，并通過系數(shù)疊加獲得重力場的球冠諧模型。

1 球冠諧系數(shù)的點(diǎn)質(zhì)量模型

點(diǎn)質(zhì)量模型是地球重力場的一種數(shù)值逼近方法。Paul［6］和 Sunkel等［7］對點(diǎn)質(zhì)量模型的構(gòu)建作了深入研究。根據(jù) Keldysh-Lavrentiev 定理［8－9］，擾動位可以用球外的一個調(diào)和函數(shù)逼近。設(shè)虛擬場的點(diǎn)質(zhì)量mi(i=1，2，…，N)分布在一個半徑為RB的Bjerhanmmar球?qū)由希晌化B加原理，點(diǎn)質(zhì)量的擾動場源在邊界上產(chǎn)生的擾動位可表示為:

式中，G是萬有引力常數(shù)，lpi是計(jì)算點(diǎn)P到點(diǎn)質(zhì)量mi的距離。根據(jù)第三類邊值問題的邊界方程［10］，地面上觀測點(diǎn)的重力異常值為:

根據(jù)非整階勒讓德函數(shù)的加法公式［11－12］，球冠外點(diǎn)P到球冠點(diǎn)Q距離的倒數(shù)的展開式為:

2 使用迭代法構(gòu)造球冠諧系數(shù)

設(shè)地面上有一系列重力異常觀測值ΔG={…，Δgj，…}T(j=1，…，M)，根據(jù)第三類邊值問題，球冠諧重力異常計(jì)算公式為:

式中(θj，λj)是球冠上第j個點(diǎn)的余緯和經(jīng)度，是第j個點(diǎn)的重力異常觀測值。應(yīng)用式(9)、(10)，使用迭代方法計(jì)算球冠諧系數(shù)CS的計(jì)算步驟如圖1所示。

圖1 迭代法計(jì)算球冠諧系數(shù)Fig.1 Spherical cap harmonic coefficient calculated by using the iterative method

3 球冠半角與點(diǎn)質(zhì)量密度和埋藏深度的關(guān)系

Bowin［16］利用等效質(zhì)點(diǎn)模型，給出球諧階數(shù) n與點(diǎn)質(zhì)量埋藏深度的關(guān)系式為H=R/(n－1)。由球冠諧階數(shù)lk和球冠半角α的關(guān)系式lk=90°(k+0.5)/a－0.5，得球冠半角與點(diǎn)質(zhì)量最大埋藏深度的關(guān)系為 H=R/(90°(k+0.5)/a －0.5 －1)。與球諧階數(shù)類似，球冠諧階數(shù)能夠擬合的網(wǎng)格間距最小為Δθ=180°/lk。選擇北緯31.5°～34°，東經(jīng)102°～104.5°區(qū)域，球冠半角約 1.76°的覆蓋范圍作為數(shù)值計(jì)算區(qū)域，選取球冠諧階數(shù)k=20，計(jì)算得lk=1 047.8，因此點(diǎn)質(zhì)量模型的最小網(wǎng)格分辨率為Δθ=0.17°。使用EMG2008模型計(jì)算網(wǎng)格分辨率分別為0.1°×0.1°、0.2°× 0.2°、0.5°×0.5°和 0.8°×0.8°的180階重力異常觀測值，選取點(diǎn)質(zhì)量的網(wǎng)格分辨率與重力異常網(wǎng)格分辨率一致，埋藏深度從1 km到80 km，使用迭代法計(jì)算球冠諧系數(shù)，得到不同網(wǎng)格分辨率下的重力異常誤差的均方根如圖2和表1所示;點(diǎn)質(zhì)量求解系數(shù)矩陣A的條件數(shù)如圖3和表2所示。

可以看到，當(dāng)點(diǎn)質(zhì)量網(wǎng)格分辨率小于0.17°時，隨著點(diǎn)質(zhì)量埋藏深度的增加，其求解系數(shù)矩陣A的條件數(shù)迅速上升。雖然重力異常誤差均方根在埋深25 km 時達(dá)到最小0.08 ×10－5ms－2，但這時求解系數(shù)矩陣A的條件數(shù)是13 894，矩陣A是病態(tài)的，其求逆運(yùn)算結(jié)果不穩(wěn)定。當(dāng)埋藏深度小于15 km時，求解系數(shù)矩陣A的條件數(shù)小于200，其求逆運(yùn)算結(jié)果穩(wěn)定，與文獻(xiàn)［17］結(jié)論一致。當(dāng)點(diǎn)質(zhì)量網(wǎng)格分辨率大于0.17°時，求解系數(shù)矩陣A的條件數(shù)增長緩慢，而且不同埋藏深度下對應(yīng)的重力異常誤差均方根 都小于0.01 ×10－5ms－2。因此在給定球冠半角和階數(shù)的條件下，當(dāng)點(diǎn)質(zhì)量網(wǎng)格分辨率小于最小間距時，點(diǎn)質(zhì)量的埋深應(yīng)當(dāng)選取與重力異常觀測網(wǎng)格分辨率相當(dāng)?shù)某叨取．?dāng)設(shè)置的點(diǎn)質(zhì)量網(wǎng)格分辨率大于最小間距時，在點(diǎn)質(zhì)量一定埋深范圍內(nèi)，其求解系數(shù)矩陣A的條件數(shù)增長緩慢，所以能根據(jù)數(shù)據(jù)擬合需要設(shè)置合理的埋藏度，并保證迭代算法收斂。

圖2 重力異常誤差均方根(單位:10－5ms－2)Fig.2 RMS errors of gravity anomaly(unit:10 －5ms－2)

圖3 矩陣A的條件數(shù)Fig.3 Condition number of matrix A

表1 重力異常誤差均方根(單位:10－5ms－2)Tab.1 RMS errors of gravity anomaly(unit:10 －5ms－2)

表2 不同埋藏深度下矩陣A的條件數(shù)Tab.2 Condition numbers of matrix Ain different buried depth

4 應(yīng)用點(diǎn)質(zhì)量組模型構(gòu)造三層球冠諧系數(shù)

使用點(diǎn)質(zhì)量組模型分層構(gòu)造球冠諧系數(shù)可以有多種組合，本文將虛擬球?qū)釉O(shè)置在上地幔、殼幔交界處和地殼中，相應(yīng)地將點(diǎn)質(zhì)量分為A、B、C 3組模型，其埋藏深度分別設(shè)置為93、55和8 km。根據(jù)分析，選取31.5°N ～34°N，102°E ～104.5°E 為覆蓋區(qū)域計(jì)算球冠諧模型。這里位于四川和甘肅交界，地勢復(fù)雜，海拔落差大，西部地形以褶皺系為主，地殼厚度50～72 km，東部以臺干為主，地殼厚度38～50 km。因?yàn)闅めＣ芏鹊牟町悾皖l重力場變化顯著［18］。實(shí)測重力異常觀測值分辨率為3'×3'網(wǎng)格，插值出網(wǎng)格分辨率為1°×1°、20'×20'和3'×3'的重力異常。設(shè)置球冠諧半徑為1.76°，點(diǎn)質(zhì)量網(wǎng)格分辨率與重力異常觀測值一致，為1°×1°、20'×20'和3'×3'網(wǎng)格分辨率。通過點(diǎn)質(zhì)量組解算得到三層球冠諧系數(shù)，具體計(jì)算過程如下。

1)在覆蓋范圍內(nèi)使用EGM2008模型生成180階次位模型重力異常值，然后使用迭代算法計(jì)算180階次位模型對應(yīng)的球冠諧系數(shù)，記為。

5)將4組球冠諧系數(shù)疊加得到重力異常觀測值的球冠諧模型。用上述方法解算出的球冠諧系數(shù)CAS，CBS和CCS，分別計(jì)算與之對應(yīng)的3'×3'網(wǎng)格分辨率殘差重力異常值如圖4。

圖4(a)是球冠諧系數(shù)CAS對應(yīng)的殘差重力異常，(b)是球冠諧系數(shù)CBS對應(yīng)的殘差重力異常，(c)是球冠諧系數(shù)CCS對應(yīng)的殘差重力異常，可以看到3層球冠諧模型能夠從低頻到高頻逐步逼近地球重力場。

在給定球冠半角的條件下，球冠諧階數(shù)決定了能夠描述的重力場最小波長［19－20］。三層球冠諧模型的內(nèi)符合精度是利用實(shí)測重力異常觀測值與球冠諧模型計(jì)算值之差來衡量的。下面給定不同的球冠半角，對應(yīng)的重力異常覆蓋范圍如表3。選取球冠諧階數(shù) k=10，12，14，16，20，計(jì)算得本模型內(nèi)符合精度如表4。

表3 三層球冠諧模型范圍Tab.3 The coverage of three layers spherical cap harmonic model

表4 三層球冠諧模型內(nèi)符合精度(單位:10－5ms－2)Tab.4 Internal consistency accuracy of three layers spherical cap harmonic model(unit:10 －5ms－2)

圖4 殘差重力異常等值線圖Fig.4 Contour map of residual gravity anomaly

由表4可以看出，本模型擬合得到的重力異常場精度隨著球冠諧階數(shù)的提高而提高;在相同的球冠諧階數(shù)下，擬合的重力異常場范圍越小，計(jì)算得到的重力異常場精度越高。當(dāng)球冠諧階數(shù)設(shè)置為20階、重力數(shù)據(jù)擬合范圍的球冠半角為0.71°時，本模型擬合的重力異常場內(nèi)符合精度最高，均方根誤差為 ±4.65 ×10－5ms－2。實(shí)測重力異常圖和本模型計(jì)算的重力異常如圖5所示。

圖5 實(shí)測與模型計(jì)算的重力異常Fig.5 Comparison of the gravity anomaly from actual survey and the gravity anomaly computed with spherical cap harmonic model

從實(shí)驗(yàn)結(jié)果看，通過分層方式構(gòu)造球冠諧模型，根據(jù)不同層點(diǎn)質(zhì)量分辨率選取合適的埋藏深度，其系數(shù)矩陣A穩(wěn)定性好，保證迭代過程收斂，可以很好地逼近局部地球重力場。

1 彭富清，于錦海.球冠諧分析中非整階Legendre函數(shù)的性質(zhì)及其計(jì)算［J］.測繪學(xué)報(bào)，2000，29(3):204 －208.(Peng Fuqing，Yu Jinhai.The characters and computation of Legendre functionwith non-integral degree in SCHA［J］.Acta Geodaetica et Cartographica Sinica，2000，29(3):204 －208)

2 Haines G V.Spherical cap harmonic analysis［J］.J geophys Res，1985，90:2 583 – 2 591.

3 Hwang C，Chen S K.Fully normalized spherical cap harmonics:App lication to the analysis of sea 2 level data from TO PEX/POSEIDON and ERS-1［J］.Geophys J Int，1997，129:450－460.

4 張勤，趙超英.地殼垂直形變場逼近的球冠諧分析法［J］.測繪學(xué)報(bào)，2004，33(1):39 －42.(Zhang Qin，Zhao Chaoying.The spherical cap harmonic analysis methodfor crust vertical deformation field fitting［J］.Acta Geodaetica et Cartographica Sinica 2004，33(1):39 －42)

5 Li Jiancheng.Spherical cap harmonic expansion for local gravity represention［J］.Manuscripta Geodaetica，1995，(20):265－277.

6 Needham P E.The formation and evaluation of detailed geopotential models based on point masses［R］.Ohio:Ohio State University，1970.

7 Sunkel H.The generation of a mass point model from surface gravity data［R］.Ohio:Ohio State University，1983.

8 吳星.衛(wèi)星重力梯度數(shù)據(jù)處理理論與方法［D］.鄭州:信息工程大學(xué)，2009.(Wu Xing.Research on determining the global earth gravity model from satellite gravity gradients［D］.Zhengzhou:Information Engineering University，2009)

9 吳星，張傳定，王凱.衛(wèi)星重力梯度邊值問題的點(diǎn)質(zhì)量調(diào)和分析［J］.測繪學(xué)報(bào)，2011，40(2):213 －219.(Wu Xing，Zhang Chuanding，Wang Kai.Point-mass harmonic analysis of satellite gradiometry boundary value problem［J］.Acta Geodaetica et Cartographica Sinica，2011，40(2):213－219)

10 王建強(qiáng).重力三層點(diǎn)質(zhì)量模型的構(gòu)造與分析［J］.測繪學(xué)報(bào)，2010，39(5):503 － 507.(Wang Jianqiang.The construction and analysis for three-tier point mass model of gravity［J］.Acta Geodaetica et Cartographica Sinica，2010，39(5):503－507)

11 Hobson E W.The theory of spherical and ellipsoidal harmonics［M］.Cambridge University Press，1931.

12 王竹溪，郭敦仁.特殊函數(shù)概論［M］.北京:北京大學(xué)出版，1989.(Wang Zhuxi，Guo Dunren.Special functions［M］.Beijing:Peking University Press，1989)

13 Ghadi K A.Transformation of global spherical harmonic models of the gravity field to a local adjusted spherical cap harmonic model［J］.Arab J Geosci，2013(6):375 － 381)

14 王三軍.利用球冠諧模型研究GPS水準(zhǔn)問題［J］.測繪科學(xué)，2008，33(2):13 －14.(WangSanjun.Thespherical cap harmonic analysis of GPS leveling height［J］.Science of Surveying and Mapping，2008，33(2):13 －14)

15 曹月玲，王解先.利用球冠諧分析擬合GPS水準(zhǔn)高程［J］.武漢大學(xué)學(xué)報(bào):信息科學(xué)版，2008，33(5):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(5):740－743)

16 Bowin C.Depth of principal mass anomalies contributing to the earth’s geoidal undulation and gravity anomalies［J］.Marine Geodesy，1983，21(7):61 －100.

17 吳曉平.局部重力場的點(diǎn)質(zhì)量模型［J］.測繪學(xué)報(bào)，1984，13(4):249 － 258.(Wu Xiaoping.Point-mass model of local gravity field［J］.Acta Geodaetica et Cartographica Sinica，1984，13(4):249 －258)

18 高銳.由地震探測揭示的青藏高原莫霍面深度［J］.地球?qū)W報(bào)，2009，30(6):761 － 773.(Gao Rui.The Moho depth of Qinghai-Tibet plateau revealed by seismic detection［J］.Acta Geoscientica Sinica，2009，30(6):761 －773)

19 王喜臣.利用球冠諧分析方法提取不同波長重力場異常［J］.世界地質(zhì)，1996，15(3):80 － 83.(Wang Xichen.Extracting different wavelength gravity anomaly by the spherical cap harmonic analysis method［J］.Global Geology，1996，15(3):80 －83)

20 翟振和，孫中苗.渤海灣多源重力數(shù)據(jù)的自適應(yīng)融合處理［J］.測繪學(xué)報(bào)，2010，39(5):444 －449.(Zhai Zhenhe，Sun Zhongmiao.The adaptive fusion of multi-source gravity data in Bohai Gulf［J］.Acta Geodaetica et Cartographica Sinica，2010，39(5):444 －449)

STRUCTURE OF MULTILAYERED SPHERICAL CAP HARMONIC MODEL

Wang Yi and Jiang Xiaodian
(College of Marine Geo-science，Ocean University of China，Qingdao 266100)

According to the relationship between the spherical cap function and the point mass model，a model was proposed based on multilayer point mass model structure coefficient of spherical cap harmonic method.Taking a region(32°N-34°N，102°E-104°E)as the calculating area，the Gravity anomaly was calculated using the three layers spherical cap harmonic coefficient which constructed by EGM2008 model and the Gravity measured observations.Analysis of fitting precision of gravitational field under different spherical cap harmonic order indicates that the average error between the spherical cap harmonic model and the measured observations is less than 0.5 ×10－5ms－2，and internal consistency accuracy of the model is 4.56 ×10－5ms－2when Spherical cap half angle is 0.71°.

spherical cap harmonic analysis;EGM2008 model;point mass model;Iterative algorithm;ill-condition

P312

A

1671-5942(2014)05-0030-05

2014-01-03

海洋公益性項(xiàng)目(201305029-02)。

王燚，男，1974年生，博士，講師，主要研究方向?yàn)榈厍蛭锢砜碧椒椒ㄅc信息技術(shù)。E－mail:wyfox009@163.com。