Moho面擾動(dòng)重力梯度信息的提取

葉周潤(rùn)，柳林濤，梁星輝，3

1.中國(guó)科學(xué)院大學(xué)，北京100049；2.中國(guó)科學(xué)院測(cè)量與地球物理研究所大地測(cè)量與地球動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢430077；3.地理信息工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安710054

1 引 言

GOCE（Gravity Field and Steady－state Ocean Circulation Explorer）是一顆用于恢復(fù)地球重力場(chǎng)和探測(cè)靜態(tài)洋流的重力梯度衛(wèi)星，其優(yōu)于100km的空間分辨率觀測(cè)值能改善地球重力場(chǎng)的中高頻信息，因此在精化全球大地水準(zhǔn)面［1］的同時(shí)也可以為地球物理學(xué)研究帶來更多信息。GOCE被認(rèn)為可以捕捉到殼幔邊界信號(hào)［2］，對(duì)巖石圈的密度敏感［3］，因此在盆地密度、地質(zhì)構(gòu)造、Moho面［4－6］等具體研究中具有巨大潛力。

殼幔起伏和巖石圈密度是重要的地球內(nèi)部結(jié)構(gòu)研究，而Moho面的確定可以為地殼均衡和巖石圈有效彈性厚度研究提供相應(yīng)依據(jù)［7］，并且重點(diǎn)地區(qū)（如青藏高原等）的Moho面結(jié)構(gòu)可以為此區(qū)域的地球動(dòng)力學(xué)研究提供參考［8－9］。相對(duì)于地震方法，利用衛(wèi)星重力確定地殼厚度可以彌補(bǔ)地面測(cè)量重力和地震資料的缺陷［10－12］，尤其是重力梯度可以改善特定頻段信號(hào)在地球物理應(yīng)用中的精度［13］，因此在GOCE衛(wèi)星發(fā)射后非洲等地區(qū)成為研究熱點(diǎn)［3］。另外重力測(cè)量值可以和地震資料進(jìn)行聯(lián)合反演，使解算結(jié)果相互約束從而提高反演精度。同時(shí)重力資料的分析方法可以用于其他星球，如月球等無(wú)地震數(shù)據(jù)的星球地殼厚度和地質(zhì)構(gòu)造研究［14］。本文的研究重點(diǎn)主要在于如何從GOCE觀測(cè)值中提取出用于Moho面反演的初值。在理論上借鑒和改進(jìn)了前人的研究思路：①伴隨GPS定位技術(shù)發(fā)展，利用重力擾動(dòng)（純重力異常）分析擾動(dòng)源更加方便［15］，相對(duì)于重力異常，重力擾動(dòng)用于Moho面反演也顯現(xiàn)出更高精度［15－18］，因此本文研究中所提取的信號(hào)都是擾動(dòng)重力梯度值；②比點(diǎn)質(zhì)量模型正演精度更高的空間域Tesseroid單元體和頻譜域球諧分析與綜合方法被應(yīng)用于重力梯度正演［19－21］；③重力觀測(cè)值結(jié)合地殼先驗(yàn)?zāi)Ｐ偷闹鸩絼冸x方案被應(yīng)用于重力梯度信號(hào)的提取中［22－23］，該處理策略減少了 Moho面反演中由于地殼密度不均帶來的誤差影響。在具體計(jì)算中采用了最新CRUST1.0先驗(yàn)?zāi)Ｐ停?4］。本文最后運(yùn)用上述方法提取的擾動(dòng)垂直重力梯度值，進(jìn)行了全球Moho面的恢復(fù)工作，并與現(xiàn)有模型進(jìn)行了對(duì)比驗(yàn)證。

2 方法介紹

2.1 Moho面信息提取策略

如圖1所示，Moho面以上的地球結(jié)構(gòu)大體分層主要有陸地巖石（1）、海水（2）、冰層（3）、沉積層（4）和除上述之外的密度不均勻地殼（5）。

如果要用Vening Meinesz－Moritz［15－18］方法解算出地殼厚度，首先需要得到如圖2所示的一個(gè)單密度界面。

為此，本文采取結(jié)合地殼先驗(yàn)?zāi)Ｐ偷闹鸩絼冸x方案，這種處理策略得到的最終結(jié)果被認(rèn)為和由地震數(shù)據(jù)解算的Moho面具有很強(qiáng)的相關(guān)系數(shù)［22－23］。對(duì)于 GOCE 觀測(cè)數(shù)據(jù)而言，則可采用式（1）進(jìn)行有效信息逐步提取

圖1 Moho面上部地殼結(jié)構(gòu)Fig.1 Crustal structure above Moho

圖2 Moho面反演初始信息Fig.2 Initial information for Moho inversion

在殼幔信號(hào)分離中，依托于GOCE重力梯度模型，根據(jù)文獻(xiàn)［25］結(jié)論在此扣除17階以內(nèi)主要存在于上地幔以下的長(zhǎng)波信息。對(duì)于重力梯度正演部分，采用球坐標(biāo)系下的空間域Tesseroid和頻域球諧方法［19－21］，以克服意大利 GEMMA Moho團(tuán)隊(duì)的點(diǎn)質(zhì)量模型計(jì)算誤差［26］。在地殼信息校正部分，采用2013年7月最新發(fā)布的CRUST 1.0模型。該模型的陸地地形，海洋和冰層資料來自ETOPO1［27］，同時(shí)沉積層和地殼密度的空間域分辨率也從2°提高到1°［24］。

2.2 GOCE模型計(jì)算擾動(dòng)重力梯度值

從GOCE球諧位系數(shù)模型恢復(fù)擾動(dòng)重力位信息可以采用式（2）的形式

式中，r、θ和λ分別表示計(jì)算點(diǎn)的球心距離、地心緯度和經(jīng)度；GM代表萬(wàn)有引力常數(shù)和地球總質(zhì)量的乘積（3.986 005e14m3/s2）；R是地球平均半徑（6371km）；Cnm和Snm是 GOCE模型位系數(shù)和正常位系數(shù)相減后得到的擾動(dòng)位系數(shù)；是完全正則化的連帶勒讓德函數(shù)；Nmin和Nmax分別為起始和最大階數(shù)。

考慮到便于后續(xù)空間域和頻譜域方法比較，這里選用當(dāng)?shù)赜沂肿鴺?biāo)系LNOF（local northoriented frame）。如果要轉(zhuǎn)換到此坐標(biāo)系，則可運(yùn)用式（3）表述［28］

式（3）中對(duì)應(yīng)的關(guān)于位的一階和二階導(dǎo)數(shù)詳見文獻(xiàn)［28］。

2.3 空間域和頻譜域正演擾動(dòng)重力梯度

2.3.1 空間域Tesseroid單元體方法

在空間域計(jì)算中，筆者選取與地球表面地形非常吻合的Tesseroid單元體作為積分模型。根據(jù)文獻(xiàn)［20］模擬試驗(yàn)結(jié)論，相對(duì)于點(diǎn)質(zhì)量模型，Tesseroid單元體正演精度可以提高3個(gè)量級(jí)；另外文獻(xiàn)［14］中試驗(yàn)表明0.25°的 Tesseroid單元體對(duì)理論月球模型的擬合誤差為0.001%，并且文獻(xiàn)［19］和［20］中已經(jīng)驗(yàn)證了該單元體用于扣除GOCE衛(wèi)星擾動(dòng)物質(zhì)影響的可行性。

文中計(jì)算公式來自泰勒級(jí)數(shù)展開方法［19］，若取Tesseroid單元體的中心點(diǎn)并展開到4階無(wú)窮小，則經(jīng)過相應(yīng)數(shù)學(xué)計(jì)算，可以得到

式中，計(jì)算基準(zhǔn)面為半徑R球面；G表示萬(wàn)有引力常數(shù)；ρ表示擾動(dòng)物質(zhì)密度；（r，φ，λ）、表示流動(dòng)計(jì)算點(diǎn)、積分點(diǎn)、Tesseroid中心點(diǎn)、上限和下限點(diǎn)的球心距離、余緯和經(jīng)度；Δr、Δφ、Δλ分別表示Tesseroid在球心、緯度和經(jīng)度方向的間距分別表示在球心、緯度和經(jīng)度方向的Γij函數(shù)二階導(dǎo)數(shù)位于Tesseroid中心點(diǎn)時(shí)的數(shù)值。o（Δ4）表示泰勒展開的截?cái)嗾`差。

由于從引力位表達(dá)式泰勒展開到二階導(dǎo)數(shù)（引力梯度）比較繁瑣，文獻(xiàn)［19］將核函數(shù)分開求導(dǎo)，則有

式中，φ表示緯度，且φ＝90°－θ；i，j，k∈｛1，2，3｝；Δxi、Δxj、Δxk分別代表x、y、z方向的間距；并且當(dāng)i＝j(luò)時(shí)，δij＝1；i≠j時(shí)，δij＝0。

式（6）中其他表達(dá)式有

式（7）—式（10）中，l和ψ分別表示積分點(diǎn)和計(jì)算點(diǎn)之間的球面距離和夾角；其他變量具體表達(dá)詳見文獻(xiàn)［19］。

2.3.2 頻譜域球諧分析與綜合方法

設(shè)一物體的內(nèi)部質(zhì)點(diǎn)分布連續(xù)，則該物體外一點(diǎn)的引力位為

可將距離l函數(shù)進(jìn)行球諧展開，然后對(duì)進(jìn)行積分，進(jìn)行級(jí)數(shù)3次項(xiàng)展開，并假定用球面代替大地水準(zhǔn)面，則經(jīng)過運(yùn)算后有如文獻(xiàn)［29］形式的表達(dá)

式中

式中，i∈｛1，2，3｝；ρ表示地球平均密度（5500kg/m3）；hup和hlw分別表示積分體高程上下限，具體數(shù)值以球面（大地水準(zhǔn)面）為基準(zhǔn)，即球面以上為正，反之為負(fù)；表示球諧譜系數(shù)，可由式（15）計(jì)算得到，將代入式（14）即可得到最終計(jì)算需要譜系數(shù)。則上述就是由擾動(dòng)物質(zhì)通過球諧譜域計(jì)算擾動(dòng)位的一般表達(dá)形式，之后的各重力梯度分量則可由式（3）得到。

2.4 Moho面恢復(fù)譜域計(jì)算公式

為了驗(yàn)證前述Moho面擾動(dòng)重力梯度信息的提取結(jié)果，筆者利用Γzz分量進(jìn)行了試驗(yàn)計(jì)算。從Γzz分量進(jìn)行Moho面恢復(fù)的頻譜域公式為

式中，dσ′＝sinθ′dθ′dλ′；D（θ，λ）和D（θ′，λ′）分別為計(jì)算點(diǎn)和流動(dòng)點(diǎn)的Moho面深度為頻譜域的Moho面解算初值。

3 試驗(yàn)計(jì)算及其結(jié)果分析

本文試驗(yàn)原始重力梯度數(shù)據(jù)來源于GOCO03S模型，該模型在低階以GRACE為主，高階采用GOCE從2009年11月開始?xì)v時(shí)18個(gè)月的中高頻數(shù)據(jù)［30］。在球諧合成時(shí)，正常重力場(chǎng)參考GRS80模型。選定GOCE軌道平均值255km為計(jì)算高度，并忽略空氣微小值影響［21］。巖石、海水和冰層密度分別取值2670kg/m3、1025kg/m3和920kg/m3。沉積層和地殼密度數(shù)據(jù)全部來自CRUST1.0模型。綜合GOCE觀測(cè)數(shù)據(jù)的空間分辨率和地殼先驗(yàn)?zāi)Ｐ途?，文中?jì)算結(jié)果分辨率設(shè)為1°×1°?？紤]到空間域和球諧譜對(duì)應(yīng)關(guān)系，式（2）和式（12）中n最大值取180。

圖3表示以Γzz為例的空間域與頻譜域方法在隨網(wǎng)格分辨率改變時(shí)的計(jì)算所需時(shí)間和兩者結(jié)果差異變化圖（648個(gè)10°×10°等間距采樣點(diǎn)）。左側(cè)縱軸表示計(jì)算時(shí)間，右側(cè)縱軸表示兩種方法絕對(duì)值之差的平均值。如圖3所示，隨著網(wǎng)格分辨率的加密，空間域方法時(shí)間需求遠(yuǎn)高于頻譜域方法，并且兩者的計(jì)算時(shí)間差異呈非線性增長(zhǎng)的方式。同時(shí)，計(jì)算方法的差異表現(xiàn)出隨格網(wǎng)分辨率加密而減少的趨勢(shì)。當(dāng)格網(wǎng)點(diǎn)加密到1°×1°時(shí)，采樣點(diǎn)絕對(duì)值之差的平均值僅為0.003 8E（1E＝10－9/s2）。表1所展示的是空間域 Tesseroid和頻譜域球諧分析與綜合方法的結(jié)果比較。表2所示為Moho面擾動(dòng)重力梯度信息逐步提取的結(jié)果統(tǒng)計(jì)。表中T、B、I、SED（S）和CRD（C）分別代表陸地地形、海水、冰層、沉積層和不均勻地殼。GTBISC表示扣除長(zhǎng)波影響后的擾動(dòng)重力梯度（G）通過逐步剝離步驟得到的含Moho面信息的擾動(dòng)重力梯度值。如表1所示，兩種方法的平均值差異在10－3量級(jí)左右。此外，筆者統(tǒng)計(jì)了表1中的差異值對(duì)于最終擾動(dòng)值（表2中GTBISC項(xiàng)）的相對(duì)平均幅度，其結(jié)果是：2.2%（Γxx）、3.4%（Γyy）和1.3%（Γzz）。

圖3 空間域與頻譜域方法計(jì)算時(shí)間和差異變化圖（Γzz）Fig.3 Comparison of difference on computation time and accuracy between spatial and spectral methodology（Γzz）

圖4表示在GOCE軌道平均高度扣除長(zhǎng)波影響后的Γzz分量結(jié)果（1E＝10－9/s2）。如圖4所示，雖然在球諧合成時(shí)進(jìn)行了長(zhǎng)波信息濾除，但圖4中的擾動(dòng)重力梯度圖形輪廓變化依舊分明，尤其在青藏高原和安第斯山脈地區(qū)表現(xiàn)明顯。這說明擾動(dòng)信息的來源主要是陸地地形起伏，海水補(bǔ)償密度和地幔及其上部的密度異常等。如表2所示，上面提及的擾動(dòng)源造成的Γzz影響接近1.5E。

表1 空間域與頻譜域方法結(jié)果比較Tab.1 Comparison between spatial and spectral methodology 10－9/s2

表2 Moho面擾動(dòng)重力梯度信息逐步提取結(jié)果統(tǒng)計(jì)Tab.2 Statistics of gravity gradient from Moho using step by step correction 10－9/s2

圖5反映的是經(jīng)過陸地地形、海水和冰層校正后的結(jié)果。圖形輪廓相對(duì)比圖4已然發(fā)生較大改變，這些顯著變化表明陸地和海洋地形是大成分?jǐn)_動(dòng)源。其實(shí)在地殼厚度粗略估算中，一般也只扣除陸地和海洋地形影響。這些在表2中也有所體現(xiàn)，ΓzzGTBI項(xiàng)極值較于G項(xiàng)數(shù)值已擴(kuò)大3倍以上。

圖6和圖7表現(xiàn)的是沉積層和地殼補(bǔ)償密度的正演結(jié)果。表3是相應(yīng)的結(jié)果統(tǒng)計(jì)。如圖6所示，沉積層不僅在全球范圍內(nèi)分布廣泛，而且極值區(qū)間也可達(dá)2E左右。這是由于雖然沉積層密度相對(duì)于巖石密度波動(dòng)不大，但厚度可達(dá)數(shù)千米以上。所以該項(xiàng)改正不僅在地殼厚度研究中要加以重視，而且在地幔密度反演中也是一項(xiàng)重要研究。考慮到Moho面實(shí)質(zhì)是同一密度單層厚度反演，地殼中殘余密度差影響要盡量消除，而圖7計(jì)算就是為解決此問題。如圖7所示，以青藏高原為例，此區(qū)平均密度值明顯小于正常巖石密度。如表3所示，由于地殼密度殘差造成全球重力梯度影響最大值達(dá)到近4E。

表3 沉積層和地殼補(bǔ)償密度正演結(jié)果統(tǒng)計(jì)Tab.3 Statistics of forward modelling from sediment and uneven density of crust 10－9/s2

圖8所示是經(jīng)過全部校正后的結(jié)果，也即筆者采用重力梯度測(cè)量值進(jìn)行Moho面反演的初始值。如表2所示，在GOCE平均軌道高度最終Moho面信息Γzz分量位于區(qū)間［－8.7，4.9］E。圖9表現(xiàn)的是采用文中計(jì)算的初始值解算的全球Moho面。結(jié)果顯示陸地地區(qū)Moho面Γzz擾動(dòng)多為負(fù)值，解算的Moho面多大于正常地殼厚度，海洋地區(qū)則反之。其中青藏地區(qū)由于印度板塊和歐亞板塊的擠壓，呈現(xiàn)明顯的邊緣淺，中部深的特點(diǎn)。這些結(jié)果符合均衡理論，使得Moho面的深度和地形負(fù)載呈現(xiàn)強(qiáng)相關(guān)性。

4 結(jié) 論

依據(jù)球坐標(biāo)系的空間域和頻譜域正演方法，結(jié)合地殼先驗(yàn)?zāi)Ｐ停疚挠懻摿藦腉OCE重力梯度模型提取Moho面擾動(dòng)信息的步驟和方法，并提供了適合該層面反演的最終結(jié)果。從文中的討論和數(shù)據(jù)分析可以得到如下結(jié)論：

（1）在文中1°×1°空間分辨率下，空間域和頻譜域方法結(jié)果接近，同時(shí)也驗(yàn)證了程序和計(jì)算的正確性，這兩種方法除了適合本文信息提取工作，同時(shí)也適用于其他擾動(dòng)源改正。

（2）在Moho面擾動(dòng)信息剝離過程中，陸地地形和海水的影響最大，并且沉積層和地殼密度差的校正不可忽略。

（3）在通過頻域公式恢復(fù)的Moho面結(jié)果與地震模型和基于GOCE數(shù)據(jù)已解算的Moho面結(jié)果的對(duì)比后驗(yàn)證了文中的擾動(dòng)信息提取計(jì)算的可靠性。

致謝：感謝斯圖加特大學(xué)F.Wild博士、中國(guó)地質(zhì)大學(xué)杜勁松博士、米蘭理工大學(xué)D.Sampietro博士、武漢大學(xué)Robert Tenzer教授和中科院測(cè)地所方劍教授等在該論文上的幫助。

圖4 經(jīng)過長(zhǎng)波扣除的Γzz結(jié)果Fig.4 Result after the reduction of long－length waves（Γzz）

圖5 經(jīng)過陸地地形，海水，冰層校正后的Γzz結(jié)果Fig.5 Result after the correction of topography，bathymetry and ice（Γzz）

圖6 沉積層正演結(jié)果（Γzz）Fig.6 Forward modelling from sediment（Γzz）

圖7 地殼密度差正演結(jié)果（Γzz）Fig.7 Forward modelling from inhomogeneous density of crust（Γzz）

圖8 經(jīng)過陸地地形、海水、冰層、沉積層和地殼密度差校正后的結(jié)果Fig.8 Result after the correction from topography，bathymetry，ice，sediment，inhomogeneous density of crust

圖9 Moho面解算結(jié)果Fig.9 Result of global Moho

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

［2］HIRT C，KUHN M，F(xiàn)EATHERSTONE W E，et al.Topographic/Isostatic Evaluation of New－generation GOCE Gravity Field Models［J］.Journal of Geophysical Research：Solid Earth（1978－2012），2012，117（B5），DOI：10.1029/2011JB008878.

［3］EBBING J，BOUMAN J，F(xiàn)UCHS M，et al.Sensitivity of GOCE Gravity Gradients to Crustal Thickness and Density Variations：Case Study for the Northeast Atlantic Region［M］∥ MARTI U.Gravity，Geoid and Height Systems.Berlin：Springer，2014：291－298.

［4］EVERTON P B.Ouso dos Dados da Miss?o GOCE Para a Caracteriza??o ea Investiga??o das Implica??es na Estrutura de Densidade das Bacias Sedimentares do Amazonas e Solim?es，Brasil［D］.S?o Paulo：Universidade de S?o Paulo，2012.

［5］áLVAREZ O，GIMENEZ M，BRAITENBERG C，et al.GOCE Satellite Derived Gravity and Gravity Gradient Corrected for Topographic Effect in the South Central Andes Region［J］.Geophysical Journal International，2012，190（2）：941－959.

［6］REGUZZONI M，SAMPIETRO D.Moho Estimation Using GOCE Data：A Numerical Simulation［M］∥KENYON S，PACINO M C，MARTI U.Geodesy for Planet Earth.Berlin：Springer，2012：205－214.

［7］YANG Ting，F(xiàn)U Rongshan，HUANG Jinshui.On the Inversion of Effective Elastic Thickness of the Lithosphere with Moho Relief and Topography Data［J］.Chinese Journal of Geophysics，2012，55（11）：3671－3680.（楊亭，傅容珊，黃金水.利用 Moho面起伏及地表地形數(shù)據(jù)反演巖石圈有效彈性厚度的莫霍地形導(dǎo)納法（MDDF）［J］.地球物理學(xué)報(bào)，2012，55（11）：3671－3680.）

［8］KE Xiaoping，WANG Yong，XU Houze.Moho Depths Inversion of Qinghai－Tibet Plateau with Variable Density Model［J］.Geomatics and Information Science of Wuhan University，2006，31（4）：289－292.（柯小平，王勇，許厚澤.用變密度模型反演青藏高原的莫霍面深度［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2006，31（4）：289－292.）

［9］XING Lelin，SUN Wenke，LI Hui，et al.Present－day Crust Thickness Increasing Beneath the Qinghai－Tibetan Plateau by Using Geodetic Data at Lhasa Station［J］.Acta Geodaetica et Cartographica Sinica，2011，40（1）：41－44.（邢樂林，孫文科，李輝，等.用拉薩點(diǎn)大地測(cè)量資料檢測(cè)青藏高原地殼的增厚［J］.測(cè)繪學(xué)報(bào)，2011，40（1）：41－44.）

［10］ZHANG Chijun，REN Kang.The Undulation between Core and Mantle of Earth from Disturbing Potential［J］.Chinese Journal of Geophysics，1994，37（1）：115－119.（張赤軍，任康.由擾動(dòng)位確定核幔起伏［J］.地球物理學(xué)報(bào)，1994，37（1）：115－119.）

［11］FANG Jian.Global Crustal and Lithospheric Thickness Inversed by Using Satellite Gravity Data［J］.Crustal Deformation and Earthquake，1999，19（1）：26－31.（方劍.利用衛(wèi)星重力資料反演地殼及巖石圈厚度［J］.地殼形變與地震，1999，19（1）：26－31.）

［12］SHIN Y H，XU H Z，BRAITENBERG C，et al.Moho Undulations beneath Tibet from GRACE－integrated Gravity Data［J］.Geophysical Journal International，2007，170（3）：971－985.

［13］HU Minzhang，LI Jiancheng，XING Lelin.Global Bathymetry Model Predicted from Vertical Gravity Gradient Anomalies［J］.Acta Geodaetica et Cartographica Sinica，2014，43（6）：558－574.（胡敏章，李建成，邢樂林.由垂直重力梯度異常反演全球海底地形模型［J］.測(cè)繪學(xué)報(bào)，2014，43（6）：558－574.）

［14］LIANG Q，CHEN C，LI Y G.3DInversion of Gravity Data in Spherical Coordinates with Application to the GRAIL Data［J］.Journal of Geophysical Research：Planets，2014，119（6）：1359－1373.

［15］TENZER R，BAGHERBANDI M.Reformulation of the Vening Meinesz－Moritz Inverse Problem of Isostasy for Isostatic Gravity Disturbances［J］.International Journal of Geosciences，2012，3（5A）：918－929.

［16］WANG Hansheng，CHEN Xue，YANG Hongzhi.An Iterative Method for Inversion of Deep Large－scale Single Density Interface by Using Gravity Anomaly Data［J］.Chinese Journal of Geophysics，1993，36（5）：643－650.（汪漢勝，陳雪，楊洪之.深部大尺度單一密度界面重力異常迭 代 反 演［J］.地 球 物 理 學(xué) 報(bào)，1993，36（5）：643－650.）

［17］SJ?BERG L E.Solving Vening Meinesz－Moritz Inverse Problem in Isostasy［J］.Geophysical Journal International，2009，179（3）：1527－1536.

［18］SJ?BERG L E.On the Isostatic Gravity Anomaly and Disturbance and Their Applications to Vening Meinesz－Moritz Gravimetric Inverse Problem［J］.Geophysical Journal International，2013，193（3）：1277－1282.

［19］GROMBEIN T，SEITZ K，HECK B.Untersuchungen zur Effizienten Berechnung Topographischer Effekte auf den Gradiententensor am Fallbeispiel der Satellitengradiometriemission GOCE［M］.［S.l.］：KIT Scientific Publishing，2010.

［20］WILD F，HECK B.Topographic and Isostatic Reductions for Use in Satellite Gravity Gradiometry［C］∥ XU P L，LIU J N，DERMANIS A.VI Hotine－Marussi Symposium on Theoretical and Computational Geodesy.Berlin：Springer，2008：49－55.

［21］NOVáK P，GRAFAREND E W.The Effect of Topographical and Atmospheric Masses on Spaceborne Gravimetric and Radiometric Data［J］.Studia Geophysica et Geodaetica，2006，50（4）：549－582.

［22］TENZER R，HAMAYUN K，VAJDA P.Global Maps of the CRUST 2.0Crustal Components Stripped Gravity Disturbances［J］.Journal of Geophysical Research：Solid Earth（1978－2012），2009，114（B5），DOI：10.1029/2008 JB006016.

［23］TENZER R，HAMAYUN，NOVáK P，et al.Global Crust－Mantle Density Contrast Estimated from EGM2008，DTM2008，CRUST2.0，and ICE－5G［J］.Pure and Applied Geophysics，2012，169（9）：1663－1678.

［24］LASKE G，MASTERS G，MA Z T，et al.Update on CRUST1.0－A 1－degree Global Model of Earth’s Crust［C］∥Proceedings of EGU General Assembly Conference Abstracts.Vienna，Austria：［s.n.］，2013：15－26.

［25］BAGHERBANDI M，SJ?BERG L E.Non－isostatic Effects on Crustal Thickness：A Study Using CRUST2.0in Fennoscandia［J］.Physics of the Earth and Planetary Interiors，2012，200－201（1）：37－44.

［26］REGUZZONI M，SAMPIETRO D.Moho Estimation Using GOCE Data：A Numerical Simulation［M］∥KENYON S，PACINO M C，MARTI U.Geodesy for Planet Earth.Berlin：Springer，2012：205－214.

［27］AMANTE C，EAKINS B W.ETOPO1 1Arc－minute Global Relief Model：Procedures，Data Sources and Analysis［R］.National Geophysical Data Center.Washington D C：NOAA，2009.

［28］ZHU L Z.Gradient Modelling with Gravity and DEM［D］.Columbus，Ohio： The Ohio State University，2007.

［29］TSOULIS D.A Comparison between the Airy/Heiskanen and the Pratt/Hayford Isostatic Models for the Computation of Potential Harmonic Coefficients［J］.Journal of Geodesy，2001，74（9）：637－643.

［30］MAYER－GüRR T，The GOCO Consortium：The New Combined Satellite only Model GOCO03s［J］.Journal of Geodesy，2012，87（9）：843－867.