接收函數(shù)與重力聯(lián)合估計(jì)的技術(shù)與應(yīng)用研究進(jìn)展

唐晗晗，郭良輝,2*

1 中國地質(zhì)大學(xué)（北京）地球物理與信息技術(shù)學(xué)院，北京 100083

2 地質(zhì)過程與礦產(chǎn)資源國家重點(diǎn)實(shí)驗(yàn)室（中國地質(zhì)大學(xué)），北京 100083

0 引 言

地殼厚度和泊松比是研究地殼結(jié)構(gòu)和物質(zhì)組成的重要參數(shù). 其中，地殼厚度可以描述莫霍界面起伏特征和地殼增厚或減薄狀態(tài)，巖礦石泊松比表示在單向應(yīng)力作用下橫向應(yīng)變與縱向應(yīng)變比值的絕對值，可以描述地殼物質(zhì)的長英質(zhì)或鐵鎂質(zhì)礦物含量及可能的熔融狀態(tài)（Watanabe et al., 1993;Christensen, 1996）. 在以往的地質(zhì)、地球物理研究中，眾多學(xué)者利用這兩個(gè)參數(shù)量化地殼結(jié)構(gòu)、物質(zhì)成分和變形特征，分析地殼增厚或減薄模式（Dugda et al., 2005; Ji et al., 2009; He et al., 2014;Wang et al., 2016; Cheng et al., 2021）、地殼均衡狀態(tài)（Wang et al., 2017; Guo et al., 2019），區(qū)分構(gòu)造單元或識別構(gòu)造塊體之間的縫合帶（Lombardi et al., 2008; Guo et al., 2019），結(jié)合地?zé)釘?shù)據(jù)描述地殼水化狀態(tài)（Lowry and Pérez-Gussinyé, 2011; Ma and Lowry, 2017）等. 特別地，泊松比在描述物質(zhì)成分及差異方面有重要的作用，獲取精確的地殼泊松比分布將對認(rèn)識地球構(gòu)造演化具有重要意義.

遠(yuǎn)震接收函數(shù)H-κ疊加方法（Zhu and Kanamori, 2000）是獲取地殼厚度和波速比的重要手段. 該方法避免了因拾取震相到時(shí)不準(zhǔn)確而引起的誤差，被廣泛應(yīng)用于全球多個(gè)地區(qū)（Lombardi et al., 2008; Ji et al., 2009; Bashir et al., 2011; Liu and Niu, 2011; He et al., 2014; Tao et al., 2014; Zhang et al., 2018）. 然而，實(shí)際應(yīng)用中部分臺站下方往往地殼結(jié)構(gòu)比較復(fù)雜（如構(gòu)造界面傾斜、殼內(nèi)低速層等），以及存在地震波形記錄信噪比低等問題，接收函數(shù)H-κ疊加方法拾取地殼厚度和波速比參數(shù)的不確定性增大. 因此，需要引入其它數(shù)據(jù)的約束，改善接收函數(shù)H-κ疊加譜錯(cuò)亂的問題，降低估計(jì)的不確定性.

Lowry和Gussinye（2011）首先提出了接收函數(shù)與重力聯(lián)合估計(jì)地殼厚度和波速比參數(shù)的方法，該方法認(rèn)為布格重力異常包含莫霍界面起伏、殼內(nèi)密度不均勻異常體和殼幔熱異常引起局部物質(zhì)膨脹這三部分引起的重力異常. 主要思路為：（1）利用遠(yuǎn)震接收函數(shù)H-κ疊加法計(jì)算目標(biāo)臺站及其附近臺站（以目標(biāo)臺站為窗口中心）的地殼厚度和波速比分布；（2）根據(jù)接收函數(shù)得到的地殼厚度和波速比正演得到區(qū)域理論重力異常，與實(shí)測重力異常作差并利用似然比（likelihood ratio method, Beck and Arnold, 1977）的方法計(jì)算重力H-κ似然譜；（3）根據(jù)克里金插值原理計(jì)算最優(yōu)插值概率得到“OI”（optimal interpolations）似然譜，并將其與重力H-κ似然譜、接收函數(shù)H-κ譜相乘并作歸一化處理，譜峰值對應(yīng)聯(lián)合估計(jì)的最優(yōu)解. 該方法綜合考慮了地?zé)嵋蛩匾约爸車_站數(shù)據(jù)分布結(jié)果的影響，更精確地對臺站下方的地殼厚度和波速比進(jìn)行了估計(jì)，并成功應(yīng)用到美國西部地區(qū)，得到了該地區(qū)更可靠的地殼厚度和波速比分布. 隨后，該方法經(jīng)過國內(nèi)外學(xué)者多次的發(fā)展和改進(jìn)，逐步形成了一套相對完備的接收函數(shù)和重力的聯(lián)合估計(jì)技術(shù).

本文將對接收函數(shù)與重力聯(lián)合估計(jì)的技術(shù)原理與應(yīng)用進(jìn)展進(jìn)行概述. 首先，分別簡述單一接收函數(shù)、重力似然函數(shù)估計(jì)的方法原理及兩者聯(lián)合估計(jì)技術(shù)的算法流程；再者，介紹接收函數(shù)與重力聯(lián)合估計(jì)的技術(shù)進(jìn)展，分析方法存在的不確定性及其對應(yīng)的優(yōu)化和改進(jìn)；然后，介紹該方法的應(yīng)用實(shí)踐；最后，對當(dāng)前接收函數(shù)與重力聯(lián)合估計(jì)方法進(jìn)行總結(jié)，并對其未來研究方向提出展望.

1 接收函數(shù)與重力聯(lián)合估計(jì)的原理

1.1 接收函數(shù)H-κ疊加法

接收函數(shù)是去除震源、地震波傳播路徑以及儀器響應(yīng)等因素后的時(shí)間序列，它主要記錄了地震臺站下方地殼和上地幔速度間斷面所產(chǎn)生的轉(zhuǎn)換波（如圖1，Ps）及其多次波（PpPs, PpSs+PsPs）的信息（Langston, 1979）. 震中距在 30°～90° 之間的遠(yuǎn)震 P 波可視為垂直入射，其在莫霍面上的轉(zhuǎn)換點(diǎn)對應(yīng)臺站正下方，因此，可以利用接收函數(shù)記錄的震相計(jì)算臺站正下方的地殼厚度與波速比.

圖1 遠(yuǎn)震P波接收函數(shù)示意圖Fig. 1 Sketch map of teleseismic P-wave receiver functions

Zhu 和 Kanamori（2000）提出的接收函數(shù)H-κ疊加法是目前獲取地殼厚度和波速比的使用最廣泛的方法，在理想的單層地殼模型下，不同震相的走時(shí)差與地殼厚度H、波速（橫波VS，縱波VP）關(guān)系為：

式中，P表示射線參數(shù)，通過不同數(shù)值的（H,κ）組合，計(jì)算每個(gè)震相的到時(shí)差及其對應(yīng)的接收函數(shù)振幅r（t），賦予不同的權(quán)重再求和：

即可得到H-κ疊加圖（如圖2a所示）. 其中w1、w2、w3為各項(xiàng)權(quán)重系數(shù)，和為1. 疊加譜值最大點(diǎn)對應(yīng)的（H,κ）即為最優(yōu)解.

該方法通過疊加的方式減少數(shù)據(jù)中非相干噪聲的影響，也避免了人工拾取到時(shí)引起的誤差，具備處理大量地震事件數(shù)據(jù)等優(yōu)點(diǎn). 試驗(yàn)表明，在單層水平均勻地殼模型下，接收函數(shù)H-κ疊加結(jié)果受相關(guān)輸入性參數(shù)（如初始P波速度、震相的權(quán)重分配）的影響較小，準(zhǔn)確度較高（Ogden et al., 2019; 查小惠等, 2020）. 但在實(shí)際情況中，當(dāng)臺站下方地殼結(jié)構(gòu)復(fù)雜（如殼內(nèi)低速層、界面傾斜等問題的影響），應(yīng)用接收函數(shù)H-κ疊加法估計(jì)結(jié)果的可靠性降低.

1.2 重力似然估計(jì)法

Shi等（2018）在Lowry和Gussinye（2011）的方法基礎(chǔ)上，考慮到小面積情況下深源熱異常偏弱可以忽略，從而從理論重力異常的正演部分和似然函數(shù)的計(jì)算部分對重力模塊的似然估計(jì)法進(jìn)行了簡化. 這里以Shi等（2018）改進(jìn)后的重力似然估計(jì)法為例說明其基本原理.

在不考慮深部熱異常影響的情況下，布格重力異常主要包括莫霍界面起伏引起的異常gMoho和殼內(nèi)密度不均勻體引起的重力異常gCrust（Blakely,1995），其中，地殼厚度可以描述莫霍面的起伏變化，波速比則可以描述地殼內(nèi)部物質(zhì)成分差異，與密度參數(shù)有一定關(guān)聯(lián)，其各個(gè)部分的正演公式可表示為（Shi et al., 2018）：

式中，ρMoho為莫霍面上下密度差，?ρCrust/?k為地殼密度相對于波速比的變化率，D=H-E為莫霍面深度，H為地殼厚度，E為地形高度，D為研究區(qū)莫霍面深度的平均值，κ為研究區(qū)波速比的平均值，F(xiàn){}和F-1{}分別表示傅里葉變換和反傅里葉變換，G為萬有引力常數(shù)，f為波數(shù). 利用這一關(guān)系，就可以將接收函數(shù)和重力結(jié)合起來.

最大似然估計(jì)法是建立在最大似然原理基礎(chǔ)上的一種統(tǒng)計(jì)方法，在作參數(shù)估計(jì)的應(yīng)用中比較常用. 重力似然估計(jì)是根據(jù)理論重力異常gModel與實(shí)測重力異常gObs之間的偏差Δgres計(jì)算其似然值，進(jìn)而分析其正演模型的準(zhǔn)確程度. 假設(shè)實(shí)測重力異常和理論重力異常的殘差服從高斯分布，則其似然函數(shù)可表示為均值μ 和方差σ2的函數(shù)：

據(jù)此，通過設(shè)定同接收函數(shù)H-κ疊加譜一樣大小和步長的H值和κ值，進(jìn)行掃描并計(jì)算每個(gè)點(diǎn)的似然函數(shù)值，就可以得到重力H-κ似然譜（如圖2b所示），其最大值即對應(yīng)重力估計(jì)的地殼厚度和波速比參數(shù). 通常重力H-κ似然譜的譜條帶與接收函數(shù)H-κ疊加譜的譜條帶呈正交或交叉關(guān)系.

圖2 （a）接收函數(shù)H-κ疊加譜；（b）重力H-κ似然譜；（c）聯(lián)合H-κ疊加譜Fig. 2 （a） H-κ map from the receiver-function H-κ stacking; （b） H-κ likelihood map from the gravity inversion；（c） H-κ map from the joint inversion

1.3 接收函數(shù)與重力聯(lián)合估計(jì)算法與流程

在接收函數(shù)與重力聯(lián)合估計(jì)的算法流程中，重力模塊的計(jì)算是在一個(gè)以目標(biāo)臺站為中心的窗口內(nèi)進(jìn)行的，聯(lián)合估計(jì)算法流程圖見圖3（本節(jié)流程對應(yīng)實(shí)線流程框），具體步驟為：

圖3 接收函數(shù)與重力聯(lián)合估計(jì)算法流程圖. 實(shí)線流程框?qū)?yīng)接收函數(shù)結(jié)果作初始數(shù)據(jù)，將虛線流程框替換星號框（*）即對應(yīng)重力結(jié)果作初始數(shù)據(jù)Fig. 3 Flow chart of the joint estimation by receiver function and gravity data. Shapes in solid line refer to the algorithm when receiver function results are initial data.Substitute the shape marked with an asterisk （*） for shapes in dashed line, which refers to the algorithm when gravity results are initial data

（1）準(zhǔn)備初始模型數(shù)據(jù). 利用常規(guī)接收函數(shù)H-κ疊加方法獲取研究區(qū)內(nèi)所有臺站的地殼厚度和波速比分布，并將地殼厚度減去地形得到莫霍面深度數(shù)據(jù)，作為重力異常正演的初始模型數(shù)據(jù).

（2）計(jì)算密度參數(shù)ρMoho和?ρCrust/?k. 準(zhǔn)備研究區(qū)的布格重力異常數(shù)據(jù)，并將研究區(qū)所有臺站的地殼厚度和波速比結(jié)果插值成與重力數(shù)據(jù)坐標(biāo)相對應(yīng)的網(wǎng)格數(shù)據(jù)，根據(jù)公式（3）和（4）采用線性回歸算法計(jì)算ρMoho和?ρCrust/?k的數(shù)值.

（3）計(jì)算重力H-κ似然譜. 將窗口內(nèi)所有臺站的接收函數(shù)結(jié)果網(wǎng)格化，結(jié)合第（2）步得到的ρMoho和?ρCrust/?k密度參數(shù)分別進(jìn)行正演[公式（3）和（4）]，求和即可得到窗口內(nèi)的理論重力異常.同時(shí)，截取該窗口內(nèi)的實(shí)測重力異常，兩者作差，代入似然函數(shù)[公式（3）]計(jì)算其似然值；然后，掃描同接收函數(shù)H-κ疊加譜一致的H、κ范圍和步長，計(jì)算每組（H,κ）對應(yīng)的似然值，最終生成重力H-κ似然譜.

（4）獲取目標(biāo)臺站下方（H,κ）的聯(lián)合估計(jì)值. 將計(jì)算得到的重力H-κ似然譜與接收函數(shù)H-κ疊加譜相乘，并作歸一化處理，得到聯(lián)合估計(jì)譜（圖2c），譜值最大點(diǎn)所對應(yīng)的坐標(biāo)（H,κ）即為聯(lián)合估計(jì)的最優(yōu)解.

（5）將獲取的目標(biāo)臺站結(jié)果（H,κ）對初始模型數(shù)據(jù)進(jìn)行更新，重復(fù)上述步驟計(jì)算下一個(gè)臺站的聯(lián)合估計(jì)值，直到遍歷研究區(qū)內(nèi)所有的臺站.

1.4 質(zhì)量評價(jià)

接收函數(shù)H-κ疊加譜圖提供了評估方差最大值的置信區(qū)間. 目前單一和聯(lián)合估計(jì)技術(shù)均采用常規(guī)的誤差分析方法（Zhu and Kanamori, 2000），對s（H,κ） 震相譜進(jìn)行泰勒展開并略去高階項(xiàng)，得到地殼厚度和波速比的方差Ea√ton等（2006）在此基礎(chǔ)上將振幅譜方差σs改進(jìn)為，使得方差的大小會隨疊加地震事件數(shù)目n的增加而降低，更符合疊加規(guī)律.

2 接收函數(shù)與重力聯(lián)合估計(jì)的技術(shù)進(jìn)展

自從Lowry和Pérez-Gussinyé（2011） 提出接收函數(shù)與重力聯(lián)合估計(jì)的方法以來，多位學(xué)者在該方法的發(fā)展與應(yīng)用實(shí)踐中對算法做了分析和優(yōu)化改進(jìn)，主要針對：初始模型的構(gòu)建、相關(guān)密度參數(shù)的估計(jì)、地幔熱異常的影響、地殼結(jié)構(gòu)復(fù)雜（如沉積層的存在）情況下的算法等，促進(jìn)技術(shù)實(shí)用化.

2.1 初始地殼模型優(yōu)化構(gòu)建

聯(lián)合反演對初始模型數(shù)據(jù)（地殼厚度H和波速比κ）具有依賴性，國內(nèi)外學(xué)者一般使用接收函數(shù)H-κ疊加法得到的地殼厚度和波速比（H,κ）作為重力異常正演的初始模型數(shù)據(jù). 但在實(shí)際應(yīng)用中，常規(guī)接收函數(shù)H-κ疊加法的結(jié)果往往連續(xù)性較差，模型網(wǎng)格化后易形成局部異常值，導(dǎo)致重力異常的擬合效果不佳，不符合殘差呈高斯分布的假設(shè). Shi等（2018）對接收函數(shù)估計(jì)結(jié)果異常的臺站，研究采用接收函數(shù)與重力聯(lián)合建模的辦法，優(yōu)化選取臺站的初始地殼厚度和波速比參數(shù). 我們將這一初始地殼模型構(gòu)建辦法推廣到全區(qū)，這是因?yàn)樵摲椒ǖ玫降某跏嫉貧つＰ洼^為平滑，有利于后續(xù)重力異常擬合. 這種初始模型優(yōu)化構(gòu)建算法的基本原理如下：

對布格重力異常進(jìn)行分離，得到反映莫霍界面起伏引起的區(qū)域異常部分，然后應(yīng)用頻率域密度界面反演方法（Oldenburg, 1974）對該部分重力異常進(jìn)行反演，得到莫霍面深度數(shù)據(jù)，此時(shí)再加上地形值，就得到新的初始地殼厚度HGra. 將重力反演得到的臺站地殼厚度HGra投影到臺站的接收函數(shù)H-κ疊加譜圖上，此時(shí)最大值對應(yīng)的κ值作為新的初始地殼波速比κGra. 按照這種辦法可獲取研究區(qū)內(nèi)所有臺站的初始地殼厚度和波速比參數(shù)（HGra, κGra），實(shí)現(xiàn)全區(qū)的初始地殼模型優(yōu)化構(gòu)建. 本小節(jié)對應(yīng)流程圖（圖3）的虛線流程框部分（將星號框替換為虛線流程框，即為重力結(jié)果作初始模型數(shù)據(jù)的完整聯(lián)合估計(jì)流程）.

2.2 滑動窗口參數(shù)優(yōu)化選取

馬亞偉（2017）對該方法進(jìn)行參數(shù)測試，評估計(jì)算窗口對聯(lián)合估計(jì)結(jié)果的影響，得出以下認(rèn)識：計(jì)算窗口的網(wǎng)格間距對聯(lián)合反演結(jié)果的影響不大，但為提高運(yùn)算速度，建議窗口內(nèi)網(wǎng)格間距N與臺站平均間距M的關(guān)系為M/4＜N＜M；不同的窗口大小會對地殼厚度和波速比的結(jié)果產(chǎn)生一定誤差，誤差范圍分別在1 km和0.02之內(nèi)，為保證頻率域計(jì)算的穩(wěn)定性，應(yīng)保證結(jié)果在合理范圍內(nèi)適當(dāng)增加窗口大小. 同時(shí)，對計(jì)算窗口的大小進(jìn)行說明：臺站分布較密的區(qū)域?qū)?yīng)的計(jì)算窗口較小，臺站稀疏的區(qū)域?qū)?yīng)的計(jì)算窗口較大，須確保每個(gè)窗口內(nèi)至少有5個(gè)地震臺站. 陳國江（2020）在此基礎(chǔ)上對實(shí)際臺站的窗口大小進(jìn)行調(diào)試，測試臺站為NM.BHS和LN.GAX，計(jì)算窗口的網(wǎng)格間距為20 km，固定ρMoho和?ρCrust/?k均為400 kg/m3，分別進(jìn)行了步長為40 km、變化范圍從200～600 km的窗口測試.結(jié)果顯示窗口太大或太小均會導(dǎo)致聯(lián)合反演拾取的結(jié)果不穩(wěn)定（誤差在±1 km和0.02之間），并認(rèn)為在實(shí)際數(shù)據(jù)處理過程中，使計(jì)算結(jié)果穩(wěn)定的窗口大小應(yīng)設(shè)置在300～450 km范圍之間.

2.3 地幔熱異常壓制

本文在1.2節(jié)中介紹了Shi等（2018）不考慮地幔地?zé)嵋蛩兀瑢χ亓Ξ惓５恼莺蛿M合作出簡化，并采用重力最大似然估計(jì)法代替似然比（Beck and Arnold, 1977），完成了重力H-κ似然譜的求算. 當(dāng)研究區(qū)域較小，深部地幔熱變化引起的重力異常較小且平緩，因此，該方法對小面積研究區(qū)是適用的.

對于較大面積的研究區(qū)，Guo等（2019）研究使用重力垂直梯度異常gz代替布格重力異常進(jìn)行聯(lián)合估計(jì)算法改進(jìn)，此時(shí)公式（3）、（4）變?yōu)椋?/p>

梯度異常的優(yōu)點(diǎn)是能夠削弱深部背景場，且突出淺部場源的影響，因此一定程度上可以壓制深部地幔地?zé)嶙兓鸬闹亓Ξ惓?，使得接收函?shù)與重力聯(lián)合估計(jì)方法在較大面積研究區(qū)得以推廣應(yīng)用，具有更廣的適用性.

2.4 密度參數(shù)優(yōu)化估計(jì)

為了測試密度參數(shù)ρMoho和?ρCrust/?k對聯(lián)合估計(jì)值的影響，陳國江（2020）利用實(shí)測臺站數(shù)據(jù)（NM.BHS和LN.GAX）對這兩個(gè)參數(shù)進(jìn)行測試，結(jié)果顯示，莫霍面上下密度差ρMoho對聯(lián)合反演結(jié)果的影響較小，ρMoho在100～600 g/cm3范圍內(nèi)變化，聯(lián)合反演的估計(jì)值基本保持不變；地殼平均密度與波速比的偏導(dǎo)數(shù)?ρCrust/?k在100～1 600 g/cm3范圍內(nèi)變化時(shí)，兩個(gè)臺站的結(jié)果展示出不同的變化規(guī)律，NM.BHS不受該參數(shù)的影響，而LN.GAX在200 g/cm3以內(nèi)變化對聯(lián)合估計(jì)結(jié)果（H,κ）產(chǎn)生影響，大于200 g/cm3之后結(jié)果基本保持不變. 因此，?ρCrust/?k參數(shù)對不同臺站結(jié)果的影響不同，具體規(guī)律還須進(jìn)一步明確.

密度參數(shù)?ρCrust/?k表示地殼平均密度與波速比的偏導(dǎo)數(shù)，可以用來指示地下主要巖石類型的密度與波速比的相關(guān)變化關(guān)系. Lowry和Pérez-Gussinyé（2011）利用重力異常反演計(jì)算得到美國西部的?ρ/?κ為460 kg/cm3，Zhang等（2020）根據(jù)相關(guān)系數(shù)法計(jì)算得到的美國西部?ρ/?κ為1 180 kg/cm3，中東部為300 kg/cm3. 可知，該參數(shù)的不確定性較大，使用一個(gè)常數(shù)代表一個(gè)地區(qū)是不合適的，需要計(jì)算?ρ/?κ的二維分布. 為此，我們對密度參數(shù)?ρCrust/?κ的估計(jì)做了理論分析和算法改進(jìn)，以適應(yīng)研究區(qū)不同構(gòu)造單元參數(shù)值差異的實(shí)際情況.

Ji 等（2009）在實(shí)驗(yàn)室下測量主要巖石類型的平均密度和波速比，發(fā)現(xiàn)兩者存在相關(guān)變化關(guān)系.我們研究發(fā)現(xiàn)密度參數(shù)?ρCrust/?κ的絕對值對應(yīng)這一變化趨勢線的斜率大小、及?ρCrust/?κ的正負(fù)指示密度和波速比兩個(gè)參數(shù)的正負(fù)相關(guān)關(guān)系. 因此，可以利用地殼視密度和平均波速比的分布計(jì)算?ρCrust/?κ的橫向變化值. 我們定義其絕對值?ρ/?κ=，其中Δρx和Δρy指地殼視密度在水平x、y方向的變化（Δκx和Δκy同理），如果某一點(diǎn)的視密度與波速比同時(shí)高于或低于區(qū)域平均值（如2 670 kg/m3和 1.75），該點(diǎn)?ρ/?κ取為正，反之，?ρ/?κ取負(fù). 為了保證?ρ/?κ參數(shù)分布的平滑性，可以采取計(jì)算周圍一定范圍內(nèi)的平均值并進(jìn)行滑動窗口的方式.

通常，密度參數(shù)?ρCrust/?κ的變化范圍在-1 600～1 600 kg/cm3范圍之間（零值除外），在盆地地區(qū)以負(fù)值為主，指示較高的密度背景值；在山脈地區(qū)以正值為主，指示較低的密度背景值.

2.5 地殼結(jié)構(gòu)復(fù)雜情況下的算法優(yōu)化

本文在1.1節(jié)提到，在單層、水平均勻的地殼模型假設(shè)下，接收函數(shù)H-κ疊加結(jié)果受初始P波速度數(shù)值的影響較?。∣gden et al., 2019; 查小惠等,2020）. 在實(shí)際情況中，部分臺站下方不滿足這一理想假設(shè)，呈多層或介質(zhì)不均勻等復(fù)雜結(jié)構(gòu)，此時(shí)，地殼初始P波速度的錯(cuò)誤估計(jì)會對接收函數(shù)的疊加結(jié)果產(chǎn)生擾動，進(jìn)而對聯(lián)合估計(jì)結(jié)果產(chǎn)生影響. 為了減少地殼初始P波速度帶來的誤差，任志遠(yuǎn)和李永華（2020）在接收函數(shù)和重力聯(lián)合估計(jì)的基礎(chǔ)上加入面波頻散的約束，獲取臺站下方地殼平均P波速度、地殼厚度與波速比值.

此外，模型測試和實(shí)際應(yīng)用發(fā)現(xiàn)，在沉積層較厚的地區(qū)，聯(lián)合估計(jì)的效果不佳，因此，聯(lián)合估計(jì)技術(shù)還應(yīng)對考慮沉積層的影響. 目前，前人在接收函數(shù)的沉積層校正方面提出了一些方法（Yeck et al., 2013; Tao et al., 2014; Yu et al., 2015），我們把沉積層校正技術(shù)（Yu et al., 2015）與接收函數(shù)和重力的聯(lián)合估計(jì)技術(shù)結(jié)合起來，用以削弱沉積層的影響，此時(shí)，重力異常增加了沉積層底界面起伏引起的重力異常gSed，對應(yīng)的正演公式為：

式中，ρSed表示沉積層界面上下密度差，D0和D0表示沉積層深度及其平均值. 此時(shí)，聯(lián)合估計(jì)得到的地殼厚度和波速比結(jié)果指示沉積層下方的結(jié)晶基底結(jié)構(gòu)，因此，考慮沉積層的影響有助于獲取臺站下方真實(shí)的地殼結(jié)構(gòu).

3 接收函數(shù)與重力聯(lián)合估計(jì)的應(yīng)用進(jìn)展

接收函數(shù)與重力聯(lián)合估計(jì)技術(shù)可以獲取更可靠的地殼厚度和泊松比（以及波速比），為區(qū)域地殼變形、物質(zhì)組成和動力學(xué)機(jī)制提供重要依據(jù)，目前已被應(yīng)用于北美大陸和中國部分地區(qū).

3.1 美國大陸

巖石類型、溫度和是否存在結(jié)合水是地殼受力變形和發(fā)生破裂的重要影響因素（Lowry and Gussinyé, 2011）. Lowry和Pérez-Gussinyé （2011）利用接收函數(shù)、重力和地表熱流聯(lián)合估算了美國西部的地殼厚度和波速比分布，結(jié)果發(fā)現(xiàn)，波速比低值和巖石圈高溫二者與科迪勒拉山脈的變形具有驚人的相關(guān)性，利用波速比對溫度不敏感、但對石英豐度非常敏感的特性，分析認(rèn)為，大陸巖石中最弱的礦物——地殼石英的豐度可能是大陸的溫度和變形的決定性參數(shù)，韌性應(yīng)變首先位于相對較弱的石英質(zhì)地殼內(nèi)，然后促進(jìn)變暖、水化和進(jìn)一步弱化，形成一種反饋機(jī)制，進(jìn)一步解釋構(gòu)造活化和威爾遜旋回（Wilson tectonic cycle）. Ma和Lowry （2017）將接收函數(shù)與重力聯(lián)合估計(jì)技術(shù)應(yīng)用在美國大陸地區(qū)，分析認(rèn)為水化作用降低了地殼波速比和密度，并在地殼淺層釋放熱量，而在最底層吸收熱量，進(jìn)一步分析認(rèn)為科迪勒拉山脈高熱流可能反映了部分地殼水化. Zhang等（2020）利用聯(lián)合估計(jì)技術(shù)中得到的理論重力異常與實(shí)測重力異常之間的偏差，反演并發(fā)現(xiàn)了美國中部密蘇里州的布盧姆菲爾德附近相鄰臺站間的深成巖體（bloomfield pluton），該結(jié)果與Ravat等（1987）使用重、磁數(shù)據(jù)構(gòu)建的復(fù)雜密度異常體相一致. 由此可見，聯(lián)合估計(jì)技術(shù)獲取的重力異常殘差還可以幫助識別相鄰臺站之間的密度異常體.

3.2 青藏高原東南緣

李永東等（2017）通過接收函數(shù)與重力聯(lián)合估計(jì)技術(shù)估算了青藏高原東南緣的地殼厚度和波速比分布，該地區(qū)的地殼厚度整體位于30～60 km之間，由南向北逐漸增厚；地殼平均波速比在1.70～1.86之間變化，相比單一接收函數(shù)的結(jié)果，聯(lián)合反演后的結(jié)果分布更為均勻，四川盆地及攀西裂谷地區(qū)為明顯的波速比高值異常，滇南地區(qū)呈現(xiàn)出相對的低值異常. 結(jié)合水對地殼常見礦物波速比的影響，李永東等（2017）認(rèn)為青藏高原東南緣的低波速比、低強(qiáng)度特征與俯沖的巖石圈脫水及上覆巖石圈的水合作用有關(guān). 從礦物學(xué)角度分析，水的出現(xiàn)會促使地殼中石英、長石等相對較輕的礦物增長，輝石、石榴子石的含量降低，從而降低地殼波速比和密度數(shù)值，在重力均衡的作用下地表必然會被抬升，這一過程可能在青藏高原的隆升及變形中扮演著重要角色.

3.3 青藏高原東北緣

石磊等（2020）利用青藏高原東北緣及鄰區(qū)布設(shè)的密集流動地震臺陣得到的遠(yuǎn)震事件數(shù)據(jù)，結(jié)合區(qū)域重力異常資料，應(yīng)用接收函數(shù)與重力梯度聯(lián)合估計(jì)方法，獲得了青藏高原東北緣及鄰區(qū)的地殼厚度和波速比參數(shù)值及分布特征，結(jié)果顯示：該地區(qū)地殼厚度位于39～65 km之間，由西南向東北逐漸變??；地殼平均波速比整體呈低值分布，大致在1.71～1.78之間變化. 在此基礎(chǔ)上計(jì)算得到的艾里均衡深度顯示，在松潘—甘孜和西秦嶺地區(qū)下方出現(xiàn)了最大9.3 km的偏差，這兩個(gè)地區(qū)的莫霍面深度與高程的相關(guān)性較弱，說明其均衡補(bǔ)償?shù)呢暙I(xiàn)不僅來自地殼，也來自巖石圈深部的地幔. 同時(shí)，地殼厚度與波速比、地殼厚度與密度之間存在較強(qiáng)的負(fù)相關(guān)關(guān)系，支持上地殼均勻增厚或下地殼分層的假設(shè).

3.4 華南地區(qū)

Guo等（2019）在華南大陸地區(qū)應(yīng)用接收函數(shù)與重力聯(lián)合反演技術(shù)，獲取了該地區(qū)地殼厚度和泊松比的分布. 結(jié)果顯示：該地區(qū)地殼厚度位于26.1～46.5 km之間，整體呈西北厚東南薄的變化特征，沿著宜昌—吉首—百色一線呈現(xiàn)南北向的地殼厚度梯級帶，與重力梯級帶大致吻合；其均衡地殼厚度表明，揚(yáng)子塊體的地殼超過了均衡深度，補(bǔ)償過剩，江南造山帶和華夏大部分地殼未達(dá)到均衡狀態(tài)，補(bǔ)償不足. 地殼泊松比顯示，大部分地區(qū)的數(shù)值在0.20～0.31之間，西北部四川盆地及周邊地區(qū)泊松比大于0.28，推測其前寒武紀(jì)基底內(nèi)可能存在鎂鐵質(zhì)中基性火山巖；東部沿?；鹕絽^(qū)較高，在0.26～0.29之間；揚(yáng)子塊體和華夏地塊之間的江南造山帶為明顯的低值條帶，大致在0.20～0.24之間，其低值異常的原因可歸結(jié)為江南造山帶逆沖推覆構(gòu)造和南華裂谷系沉積. 同時(shí)，該低值條帶可能指示板塊碰撞拼合的位置，Guo等（2019）根據(jù)其低值分布，結(jié)合相關(guān)地質(zhì)、地球物理資料推測揚(yáng)子與華夏兩塊體之間的縫合帶位于石臺—九江—吉首—百色一線和紹興—江山—萍鄉(xiāng)—永州—貴港—北海一線之間.

3.5 華東地區(qū)

我們將接收函數(shù)與重力聯(lián)合估計(jì)方法進(jìn)一步應(yīng)用在中國大陸華東地區(qū)，結(jié)果顯示，華東地區(qū)的地殼厚度位于26.6～62.4 km之間，整體呈現(xiàn)由西向東逐漸變薄的趨勢，與地形呈鏡像關(guān)系. 地殼厚度的梯級變化帶與重力梯級帶存在相似性和差異性，兩者在大興安嶺—太行山—武陵山一線展布相似，在秦嶺—大巴和雪峰山一帶偏離較大. 該地區(qū)地殼平均泊松比范圍為0.18～0.31，高泊松比主要集中在沉積盆地及火山地區(qū)，四川盆地和鄂爾多斯盆地的基底顯示出與正磁異常相對應(yīng)的高泊松比條帶，結(jié)合巖石物理數(shù)據(jù)推測與殼內(nèi)中基性巖體的侵入有關(guān)；松遼盆地北部呈高泊松比異常，推測與新生代以來多期次的巖漿噴發(fā)有關(guān). 低泊松比主要位于華北克拉通北部，渭河地塹、江南造山帶沿線附近，與構(gòu)造塊體相對應(yīng)，利用泊松比低值可以指示構(gòu)造板塊的邊界位置的這一特點(diǎn)，我們結(jié)合前人地質(zhì)資料，推測中亞造山帶東段增生帶的延伸位置主要位于白云礦區(qū)—化德—赤峰—遼源—和龍一線和蘇尼特右旗—林西—通遼—吉林—汪清一線之間.

以上研究表明，接收函數(shù)和重力聯(lián)合估計(jì)地殼厚度和波速比的研究工作得到了較快發(fā)展，隨著方法的不斷改進(jìn)，應(yīng)用地區(qū)涉及了中國大陸和美國大陸，且獲得了更可靠的地殼結(jié)構(gòu)，能夠?yàn)閰^(qū)域地殼結(jié)構(gòu)、物質(zhì)組成及演化過程提供準(zhǔn)確的地球物理依據(jù).

4 結(jié)語與展望

接收函數(shù)與重力聯(lián)合估計(jì)方法經(jīng)過不斷發(fā)展和技術(shù)改進(jìn)，能夠適用于地殼結(jié)構(gòu)較為復(fù)雜以及地?zé)釘?shù)據(jù)短缺的地區(qū)，具有較廣的應(yīng)用前景. 該方法可以獲得臺站下方較為可靠的地殼厚度和波速比參數(shù)，為我們了解地殼結(jié)構(gòu)、物質(zhì)組成及其動力學(xué)演化等方面提供地球物理依據(jù).

然而，接收函數(shù)與重力聯(lián)合估計(jì)技術(shù)目前仍然存在一些問題，比如，在臺站分布稀少的地區(qū)，重力異常的擬合殘差較大，以及地殼各向異性結(jié)構(gòu)和地震事件方位角的變化對接收函數(shù)產(chǎn)生的影響是否可以通過重力約束作用完全消除等，是目前待解決的問題. 同時(shí)，該方法目前的誤差評估僅是對解的分布區(qū)間進(jìn)行分析，還需要從更多角度（如初始數(shù)據(jù)、相關(guān)密度參數(shù)和計(jì)算窗口引起的擾動）對聯(lián)合估計(jì)技術(shù)誤差作全面評估.

該方法可以獲得較為可靠的地殼厚度和泊松比分布，考慮到泊松比在描述地殼物質(zhì)組成方面具有重要價(jià)值，今后可以與地質(zhì)、巖石物理資料相結(jié)合，推進(jìn)區(qū)域巖性的識別與分析，更好地發(fā)揮其作用.

致謝

感謝中國地震局地球物理研究所李永華研究員、石磊副研究員的有益討論和寶貴建議. 衷心感謝兩位匿名評審專家提出的寶貴意見.