葉周潤, 柳林濤, 梁星輝, 郎駿健

1 合肥工業(yè)大學土木與水利工程學院，測繪工程系， 合肥　230009 2 中國科學院測量與地球物理研究所大地測量與地球動力學國家重點實驗室， 武漢　430077

?

垂直重力梯度反演Moho面的頻譜域公式及其應用

葉周潤1,2, 柳林濤2, 梁星輝2, 郎駿健2

1 合肥工業(yè)大學土木與水利工程學院，測繪工程系， 合肥230009 2 中國科學院測量與地球物理研究所大地測量與地球動力學國家重點實驗室， 武漢430077

摘要通過求解引力相等原則下的Fredholm積分方程，可以得到不規(guī)則單一密度界面(Moho面)的起伏.本文充分參考了前人的理論研究，推導出擾動垂直重力梯度確定Moho面深度的頻譜域表達式，該式具有二次項迭代精度.運用此公式進行了全球Moho面的恢復計算，并將該結果與CRUST1.0模型和GEMMA Moho模型進行了對比和驗證.

關鍵詞垂直重力梯度； Moho面反演； 頻譜域公式

In this work, the Fredholm integral equation of Moho depth is converted into the one in the spectrum domain. By considering the second-item influence, we provide a simple and general expression that is applicable for the vertical gravity gradients. In the investigation, we utilize the new methodology to calculate the global Moho depth. The method of the spherical harmonic synthesis is applied to raw gravity gradient data and step-by-step stripping of Moho signal from the GOCO03s model and the prior crust information, respectively. The global crustal thickness is from 6 km to 79 km. The land crust is usually thicker than the normal value. The minimum thickness of the crust is located in the ocean. Especially, because of the collision between in the India and the Eurasian plates, the Tibetan plateau has a very thick crust. With respect to the African region, the crust in the north continent is thinner than the one in the south. Finally, by comparison with CRUST1.0 and GEMMA model, it is demonstrated that our Moho solutions have a good consistency on the global scale with standard deviations of 4.6 km and 3.4 km, respectively.KeywordsVertical gravity gradient; Moho inversion; Spectral domain formula

1引言

重力和重力梯度擾動反映了物質(zhì)的密度變化和界面的起伏.在經(jīng)典重力學中，Moho面(地殼厚度)的確定主要依據(jù)Vening Meinesz-Moritz(VMM)方法(Moritz，1990).在該方法中，地形或海水等質(zhì)量被認為是一種加在不斷裂而有彈性的地殼上的負載，通過求解引力相等原則下的Fredholm積分方程即可反推出地殼的厚度.Moritz(1990)首先進行了理論推導，但在該表達式中重力異常求導的方向并沒有嚴格指向球心，從而不太符合大尺度的界面反演(Moritz，1990；汪漢勝等，1993).Sj?berg(2009)和汪漢勝等(1993)修正了導數(shù)指向的不足，并分別推導出不同形式的顧及二次項影響的迭代算法.另外，在FFT正演位場界面的基礎上(Parker，1973)， Oldenburg(1974)提出確定起伏密度界面深度的迭代法，該方法被廣泛應用于局部層面反演中，如海洋深度計算和青藏地區(qū)的地殼厚度計算(Oldenburg，1974；Braitenberg et al.，2000；柯小平，2006；Shin et al.，2007).此外，重力反演地殼界面的主要進展還包括：(1) 伴隨GPS 定位技術的發(fā)展，重力擾動逐漸取代重力異常并被應用于地球物理學研究中(Sj?berg，2013).Tenzer和Bagherbandi(2012)使用重力擾動數(shù)據(jù)得到了結果更優(yōu)的全球Moho面.Sj?berg(2013)推導出聯(lián)合重力異常和重力擾動反演Moho面的公式.(2) 隨著GOCE衛(wèi)星的成功發(fā)射，空間分辨率更高的衛(wèi)星重力梯度觀測值被用于全球地殼厚度的恢復研究(Ebbing et al.，2011；Hirt et al.，2012).Reguzzoni等(2013)首先給出一階近似下的球面迭代算法.Bagherbandi和Eshagh(2012)采用離散最小二乘方法模擬反演了伊朗地區(qū)的地殼厚度.(3) 地殼先驗模型的聯(lián)合應用.Tenzer等(2009)首先提出基于地殼先驗信息的逐步剝離策略，該方法經(jīng)過地形、海水、冰層、沉積層和地殼殘余密度差改正后，可以有效地減小密度差對地殼界面恢復的影響.它被認為是重震數(shù)據(jù)的一種融合解算方法，用來改善地震先驗模型在資料欠缺地區(qū)的不足.

Reguzzoni等(2013)雖然給出了基于垂直重力梯度反演全球地殼厚度的迭代法，但在其公式推導中的截斷誤差只考慮了一階近似.同時，關于Moho面信號的提取采用了精度較差的點質(zhì)量正演模型.為改善上述提及工作的不足，本文在Moritz，Wang，Sj?berg和Reguzzoni 等研究思路的基礎上，推導出顧及二次項影響的頻譜域迭代公式，并在正演計算過程中采用了精度更高的球諧譜方法.

2計算方法

2.1理論推導

如圖1所示，在大地水準面以下存在一個深度為D的不規(guī)則界面.假設該層面的正常厚度為D0，同時設定該界面與外界的密度差為Δρc/m，則其產(chǎn)生的重力位擾動可以表達為(通常忽略離心力的影響)(Sj?berg，2009)：

(1)

根據(jù)Heiskanen和Moritz(1967)可知，當r≥R時，式(1)中關于l-1(r,θ,λ)的表達式可以用球諧函數(shù)展開為：

(2)

式中，Pn(cosψ)表示面諧函數(shù).將式(2)代入式(1)，則有：

×Pn(cosψ)dr′dσ′.

(3)

根據(jù)定義，垂直重力梯度為重力位在球心方向的二階導數(shù)，則可以得到：

×(n+1)(n+2)Pn(cosψ)dr′dσ′.

(4)

考慮到定積分的可加性質(zhì)，(4)式中的積分項可以改寫為：

(5)

(6)考慮到：(1)τ的最大值是100 km/6371 km≈0.0157(100 km表示最大Moho面深度，6371 km表示地球平均半徑)，故進行有限級數(shù)展開即可滿足計算精度.(2)據(jù)Sj?berg (2009)指出，τ級數(shù)展開的三次項最大近似誤差是1003/63712km≈25 m，略去此項可滿足大多數(shù)情況的實際應用.因此對1-(1-τ)n+3展開到二次項即可，具體展開如下：

根據(jù)球諧譜分析理論，式(8)可以等價轉換為：

(9)

(10)考慮到有如下數(shù)學關系式成立(Heiskanen and Moritz，1967)：

(11)

式中，Yn(θ,λ)是球諧分析與綜合理論中的球面函數(shù)；V(θ,λ)和V(θ′,λ′)是對應積分點和流動點在球面上的普通變量函數(shù).綜合式(10)和(11)，進行球諧系數(shù)恢復計算，則可以得到擾動垂直重力梯度恢復Moho面的頻譜域公式：

(12)

式中，D(θ,λ)和D(θ′,λ′)分別為計算點和流動點的界面深度.雖然(12)式是個迭代公式，但由于它有類似于Molodensky公式的快速收斂性質(zhì)，Sj?berg(2009)認為該類型公式在具體實用計算時可以把它當作直接表達式.

2.2Moho面計算流程

如圖2所示，本文進行Moho面反演的基本流程如下：首先從原始觀測數(shù)據(jù)中扣除球諧函數(shù)低階項的影響，由此得到的結果主要包含Moho面以上的擾動信息；再結合地殼先驗密度模型和重力梯度正演方法，通過扣除(補償)計算得到單一密度界面的擾動值；之后加上正常地殼厚度的重力梯度值即反演初值；最后利用前面介紹的頻譜域反演公式便可解算出Moho面深度.

圖2　Moho面反演流程Fig.2　Process of Moho inversion

圖2中所提及的Moho面信號提取策略介紹如下：如圖3a所示，Moho面以上的地球結構大體分層(依次)有陸地巖石(1)、海水(2)、冰層(3)、沉積層(4)和密度不均勻的地殼(5).如果我們要用Vening Meinesz-Moritz方法解算出地殼厚度，首先需要得到如圖3b所示的一個單密度界面.為此，這里采取結合地殼先驗模型的逐步剝離方案，這種處理策略得到的最終結果被認為和地震數(shù)據(jù)解算的Moho面具有很強的相關性(Tenzer et al.，2009).對于GOCE觀測數(shù)據(jù)而言，則可采用(13)式進行有效信息逐步提?。?/p>

(13)

2.3重力梯度球諧合成和正演公式

(14)

(15)

(16)

(17)

圖3　(a) Moho面上部地殼結構示意圖； (b) Moho面示意圖Fig.3　(a) Crustal structure above Moho； (b) Schematic of Moho

3實驗計算與分析

本文原始重力梯度數(shù)據(jù)來源于GOCO03s模型，該模型具有低階GRACE和中高頻GOCE數(shù)據(jù)互補優(yōu)勢(Mayer-Gürr et al.，2012). 在球諧合成時，選取GRS80模型為正常重力場參數(shù)(Mortiz，1980)，選定GOCE軌道平均值255 km為計算高度，并考慮到低階信號主要包含地核和下地幔物質(zhì)的密度不均勻信息，依據(jù)Bagherbandi和Sj?berg(2012)、方劍(1999)結論，在此扣除16階以內(nèi)的球諧合成值.在地殼信息校正部分，采用2013年7月最新發(fā)布的CRUST1.0模型(Laske et al.，2013).巖石、海水、冰層、沉積層和地殼密度數(shù)據(jù)全部來自該模型，并取殼幔密度差為485 kg·m-3.綜合GOCE觀測數(shù)據(jù)的空間分辨率和地殼先驗模型精度，文中計算結果分辨率設為1°×1°.考慮到空間域和球諧譜對應關系，式(12)和(14)中球諧函數(shù)階次的最大值取180.

表1　擾動垂直重力梯度逐步剝離結果統(tǒng)計

圖5展示的是Moho面反演公式二次項的影響.總體而言，二次項的貢獻最大可達14 km，尤其在重力梯度變化比較劇烈的青藏高原等地區(qū).圖6展示的是全球Moho面分布圖.本文解算得到的全球地殼厚度分布區(qū)間為[6～79] km，陸地地區(qū)的地殼厚度多大于正常地殼，最小的地殼厚度區(qū)域位于海洋，該區(qū)域結果多小于正常值.其中，青藏地區(qū)由于印度板塊和歐亞板塊的擠壓，呈現(xiàn)明顯的邊緣淺和中部深特點；非洲地區(qū)則呈現(xiàn)出南部地殼厚于北部的特征.表2展示的是本文解算結果分別與CRUST1.0和GEMMA Moho模型的對比統(tǒng)計，相應的模型對比偏差平均值是-3.9 km和-2.7 km，模型之間差異的標準差是4.7 km和3.5 km.與CRUST1.0模型的差異主要集中在非洲、南美和兩極地區(qū)，這是由于該模型在這些區(qū)域缺少準確的地震數(shù)據(jù).與GEMMA Moho模型的差異主要集中在青藏高原地區(qū)，該模型在此區(qū)域解算的地殼厚度值較大，最大值可達100 km左右，該區(qū)域值均大于本文和CRUST1.0模型結果.雖然存在局部差異，但本文解算結果與CRUST1.0模型對比的平均值和標準偏差均優(yōu)于Tenzer和Bagherbandi(2012)、Reguzzoni等(2013).

表2　模型結果對比統(tǒng)計

圖5　Moho面反演公式二次項的影響Fig.5　Effect of second item in the formula of Moho inversion

圖6　全球Moho面解算結果Fig.6　Result of global Moho

4結論

本文的主要工作是推導了新的擾動垂直重力梯度確定界面深度的頻譜域算法，并用該算法進行了全球地殼厚度的反演工作.在理論創(chuàng)新部分，將求解界面深度的Fredholm積分方程轉為易于求解的頻譜域等式進行解算，所推導的界面深度反演表達形式簡潔而通用.由于考慮了級數(shù)展開二次項的影響，所得公式保證了理論上的精度.在實驗部分，結合最新的地殼先驗模型CRUST1.0，采用球諧譜正演重力梯度的方法對Moho面的信號進行了逐步剝離；并將反演得到的全球Moho面深度分別與CRUST1.0模型和GEMMA Moho模型進行了對比分析，從文中相關統(tǒng)計顯示，本文解算結果和兩種模型在全球范圍內(nèi)具有較良好的一致性.

Moho面的恢復計算在地球物理學研究中是比較重要和復雜的工作.本文的研究只是基于傳統(tǒng)重力學的研究思路，在理論算法和正演模型方面提出了合理的改進措施，在全球反演過程的其他細節(jié)方面借鑒了前人的較成熟結論，比如低階長波項的扣除，Moho面信號的提取和殼幔密度差的設定，等.并且在其他方面還需繼續(xù)深入，比如，其他重力梯度分量的界面恢復公式的理論推導，解算結果的詳細地球物理學解釋和綜合地震數(shù)據(jù)進行Moho面的聯(lián)合結算，等.致謝感謝審稿專家的意見！

References

Bagherbandi M, Sj?berg L E. 2012. Non-isostatic effects on crustal thickness: a study using CRUST2.0 in Fennoscandia.PhysicsoftheEarthandPlanetaryInteriors, 200-201: 37-44.

Bagherbandi M, Eshagh M. 2012. Crustal thickness recovery using an isostatic model and GOCE data.Earth,PlanetsandSpace, 64(11): 1053-1057.

Braitenberg C, Zadro M, Fang J, et al. 2000. The gravity and isostatic Moho undulations in Qinghai-Tibet plateau.J.Geodyn., 30(5): 489-505.Ebbing J, Bouman J, Gradmann S, et al. 2011. Use of GOCE gravity gradient data for lithospheric modeling—a case study for the NE Atlantic margin.∥ AGU Fall Meeting Abstracts. AGU.

Fang J. 1999. Global crustal and lithospheric thickness inversed by using satellite gravity data.CrustalDeformationandEarthquake(in Chinese), 19(1): 26-31.Grombein T, Seitz K, Heck B, et al. 2010. Untersuchungen zur effizienten Berechnung topographischer Effekte auf den Gradiententensor am Fallbeispiel der Satellitengradiometriemission GOCE. KIT Scientific Publishing. Heiskanen W A, Moritz H. 1967. Physical geodesy.BulletinGéodésique, 86(1): 491-492.Hirt C, Kuhn M, Featherstone W E, et al. 2012. Topographic/isostatic evaluation of new-generation GOCE gravity field models.J.Geophys.Res., 117(B5): B05407.

Ke X P. 2006. Crustal structures of Tibetan plateau from the research of gravity forward and inversion [Ph. D. thesis] (in Chinese). Wuhan: Institute of Geodesy and Geophysics, Chinese Academy of Sciences.Laske G, Masters G, Ma Z T, et al. 2013. Update on crust1. 0—a 1-degree global model of Earth′s crust.∥ EGU General Assembly. Vienna, Austria.

Liang Q, Chao C, Li Y G. 2014. 3-D inversion of gravity data in spherical coordinates with application to the GRAIL data.J.Geophys.Res., 119(6): 1359-1373.

Luo Z C. 1996. The theory and method for the determination of earth gravity field from satellite gravity gradient data [Ph. D. thesis] (in Chinese). Wuhan: Wuhan University of Surveying and Mapping Technology.

Mayer-Gürr T, Rieser D, H?ck E, et al. 2012. The new combined satellite only model GOCO03s.∥ Abstract, GGHS2012. Venice.Mortiz H. 1980. Geodetic reference system 1980 (GRS80).BulletinGéodésique, 54(3): 395-405.Moritz H. 1990. The Figure of the Earth: Theoretical Geodesy and the Earth′s Interior. Karlsruhe: Wichmann.

Novák P, Grafarend E W. 2006. The effect of topographical and atmospheric masses on spaceborne gravimetric and gradiometric data.Stud.Geophys.Geod., 50(4): 549-582.

Oldenburg D W. 1974. The inversion and interpretation of gravity anomalies.Geophysics, 39(4): 526-536.

Parker R L. 1973. The rapid calculation of potential anomalies.Geophys.J.Int., 31(4): 447-455.

Reguzzoni M, Sampietro D, Sansò F. 2013. Global Moho from the combination of the CRUST2.0 model and GOCE data.Geophys.J.Int., 195(1): 222-237.Shin Y H, Xu H Z, Braitenberg C, et al. 2007. Moho undulations beneath Tibet from GRACE-integrated gravity data.Geophys.J.Int., 170(3): 971-985.

Sj?berg L E. 2009. Solving Vening Meinesz-Moritz inverse problem in isostasy.Geophys.J.Int., 179(3): 1527-1536.

Sj?berg L E. 2013. On the isostatic gravity anomaly and disturbance and their applications to Vening Meinesz-Moritz gravimetric inverse problem.Geophys.J.Int., 193(3): 1277-1282. Tenzer R, Hamayun K, Vajda P. 2009. Global maps of the crust 2.0 crustal components stripped gravity disturbances.J.Geophys.Res., 114(B5): B05408.

Tenzer R, Bagherbandi M. 2012. Reformulation of the Vening-Meinesz Moritz inverse problem of isostasy for isostatic gravity disturbances.Int.J.Geosci., 3(5): 918-929.

Tsoulis D. 2001. A comparison between the Airy/Heiskanen and the Pratt/Hayford isostatic models for the computation of potential harmonic coefficients.J.Geod., 74(9): 637-643.

Wang H S, Chen X, Yang H Z. 1993. An iterative method for inversion of deep large-scale single density interface by using gravity anomaly data.ChineseJ.Geophys. (ActaGeophysicaSinica) (in Chinese), 36(5): 643-650.

Wild F, Heck B. 2008. Topographic and isostatic reductions for use in satellite gravity gradiometry.∥ VI Hotine-Marussi Symposium on Theoretical and Computational Geodesy. Berlin Heidelberg: Springer, 49-55.

Zhu L Z. 2007. Gradient modelling with gravity and DEM [Ph. D. thesis]. Columbus: The Ohio State University.

附中文參考文獻

方劍. 1999. 利用衛(wèi)星重力資料反演地殼及巖石圈厚度. 地殼形變與地震, 19(1): 26-31.

柯小平. 2006. 青藏高原地殼結構的重力正反演研究[博士論文]. 武漢: 中國科學院測量與地球物理研究所.

羅志才. 1996. 利用衛(wèi)星重力梯度數(shù)據(jù)確定地球重力場的理論和方法[博士論文]. 武漢: 武漢測繪科技大學.

汪漢勝, 陳雪, 楊洪之. 1993. 深部大尺度單一密度界面重力異常迭代反演. 地球物理學報, 36(5): 643-650.

(本文編輯何燕)

基金項目國家重大科學儀器設備開發(fā)專項(2011YQ120045), 國家自然科學基金(41304023)資助.

作者簡介葉周潤，男，1984年生，講師，主要從事重力梯度數(shù)據(jù)處理研究.E-mail：yezhourun329@hotmail.com

doi:10.6038/cjg20160207 中圖分類號P223

收稿日期2014-10-30，2015-12-29收修定稿

Formula of Moho inversion in the spectral domain using vertical gravity gradient and its application

YE Zhou-Run1,2, LIU Lin-Tao2, LIANG Xing-Hui2, LANG Jun-Jian2

1SchoolofCivilEngineering,InstituteofGeomaticsEngineering,HefeiUniversityofTechnology,Hefei230009,China2StateKeyLaboratoryofGeodesyandEarth′sDynamics，InstituteofGeodesyandGeophysics，ChineseAcademyofSciences,Wuhan430077,China

AbstractIn classical gravity research, the determination of Moho is mainly based on the Vening Meinesz-Moritz (VMM) method. This approach assumes that terrain and seawater can be considered as the load on the elastic crust. The Moho depth or crustal thickness can be obtained by solving the Fredholm integral equation which requires a form of expression for gravity gradient observations. However, there is no such accurate expression available now.

葉周潤, 柳林濤, 梁星輝等. 2016. 垂直重力梯度反演Moho面的頻譜域公式及其應用.地球物理學報，59(2):476-483,doi:10.6038/cjg20160207.

Ye Z R, Liu L T, Liang X H, et al. 2016. Formula of Moho inversion in the spectral domain using vertical gravity gradient and its application.ChineseJ.Geophys. (in Chinese),59(2):476-483,doi:10.6038/cjg20160207.