地球重力場球冠諧模型的分層構造和分析*

王 燚 姜效典

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

地球重力場球冠諧模型的分層構造和分析*

王 燚 姜效典

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

根據球冠諧系數與點質量模型的關系，提出一種基于多層點質量模型分層構造球冠諧系數的方法。以32°N～34°N，102°E～104°E為計算區域，利用EGM2008模型和實測觀測值構造三層球冠諧模型系數，逼近該區域的重力異常場。結果表明，使用本方法構造的球冠諧模型和實測觀測值的誤差平均值小于0.5×10－5ms－2，當擬合區域的球冠半角為 0.71°時，模型的內符合精度為 ±4.65 ×10－5ms－2。

球冠諧和分析;EGM2008模型;點質量模型;迭代算法;病態性

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

1 球冠諧系數的點質量模型

點質量模型是地球重力場的一種數值逼近方法。Paul［6］和 Sunkel等［7］對點質量模型的構建作了深入研究。根據 Keldysh-Lavrentiev 定理［8－9］，擾動位可以用球外的一個調和函數逼近。設虛擬場的點質量mi(i=1，2，…，N)分布在一個半徑為RB的Bjerhanmmar球層上，由位疊加原理，點質量的擾動場源在邊界上產生的擾動位可表示為:

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

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

2 使用迭代法構造球冠諧系數

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

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

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

3 球冠半角與點質量密度和埋藏深度的關系

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

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

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

圖3 矩陣A的條件數Fig.3 Condition number of matrix A

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

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

4 應用點質量組模型構造三層球冠諧系數

使用點質量組模型分層構造球冠諧系數可以有多種組合，本文將虛擬球層設置在上地幔、殼幔交界處和地殼中，相應地將點質量分為A、B、C 3組模型，其埋藏深度分別設置為93、55和8 km。根據分析，選取31.5°N ～34°N，102°E ～104.5°E 為覆蓋區域計算球冠諧模型。這里位于四川和甘肅交界，地勢復雜，海拔落差大，西部地形以褶皺系為主，地殼厚度50～72 km，東部以臺干為主，地殼厚度38～50 km。因為殼幔密度的差異，低頻重力場變化顯著［18］。實測重力異常觀測值分辨率為3'×3'網格，插值出網格分辨率為1°×1°、20'×20'和3'×3'的重力異常。設置球冠諧半徑為1.76°，點質量網格分辨率與重力異常觀測值一致，為1°×1°、20'×20'和3'×3'網格分辨率。通過點質量組解算得到三層球冠諧系數，具體計算過程如下。

1)在覆蓋范圍內使用EGM2008模型生成180階次位模型重力異常值，然后使用迭代算法計算180階次位模型對應的球冠諧系數，記為。

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

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

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

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

表4 三層球冠諧模型內符合精度(單位: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可以看出，本模型擬合得到的重力異常場精度隨著球冠諧階數的提高而提高;在相同的球冠諧階數下，擬合的重力異常場范圍越小，計算得到的重力異常場精度越高。當球冠諧階數設置為20階、重力數據擬合范圍的球冠半角為0.71°時，本模型擬合的重力異常場內符合精度最高，均方根誤差為 ±4.65 ×10－5ms－2。實測重力異常圖和本模型計算的重力異常如圖5所示。

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

從實驗結果看，通過分層方式構造球冠諧模型，根據不同層點質量分辨率選取合適的埋藏深度，其系數矩陣A穩定性好，保證迭代過程收斂，可以很好地逼近局部地球重力場。

1 彭富清，于錦海.球冠諧分析中非整階Legendre函數的性質及其計算［J］.測繪學報，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］.測繪學報，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 吳星.衛星重力梯度數據處理理論與方法［D］.鄭州:信息工程大學，2009.(Wu Xing.Research on determining the global earth gravity model from satellite gravity gradients［D］.Zhengzhou:Information Engineering University，2009)

9 吳星，張傳定，王凱.衛星重力梯度邊值問題的點質量調和分析［J］.測繪學報，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 王建強.重力三層點質量模型的構造與分析［J］.測繪學報，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 王竹溪，郭敦仁.特殊函數概論［M］.北京:北京大學出版，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水準問題［J］.測繪科學，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水準高程［J］.武漢大學學報:信息科學版，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 吳曉平.局部重力場的點質量模型［J］.測繪學報，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］.地球學報，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］.世界地質，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 翟振和，孫中苗.渤海灣多源重力數據的自適應融合處理［J］.測繪學報，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

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

王燚，男，1974年生，博士，講師，主要研究方向為地球物理勘探方法與信息技術。E－mail:wyfox009@163.com。