接收函數(shù)、面波與重力聯(lián)合約束地殼厚度與波速比＊

任志遠(yuǎn)李永華強(qiáng)正陽(yáng)石磊

1） 中國(guó)北京 100081 中國(guó)地震局地球物理研究所

2） 中國(guó)北京 100081 中國(guó)地震局震源物理重點(diǎn)實(shí)驗(yàn)室

引 言

地殼厚度和波速比是描述地殼結(jié)構(gòu)與物質(zhì)組成的重要參數(shù)，其中：地殼厚度不僅可為地震波成像研究中地殼改正、重力補(bǔ)償與模擬等提供基本參數(shù)，還可為地殼形成演化及其動(dòng)力學(xué)過(guò)程研究提供重要約束；地殼波速比則可為地殼物質(zhì)組成的確定提供約束（Christensen，F(xiàn)ountain，1975）.

確定地殼厚度的方法較多，包括接收函數(shù)、面波、重力等方法.接收函數(shù)是由遠(yuǎn)震水平分量與垂直分量之間的反褶積產(chǎn)生的時(shí)間序列，由于該方法消除了震源和傳播路徑的影響，更適于獲取地殼和上地幔不連續(xù)界面的信息，是獲取地殼和地幔結(jié)構(gòu)最常用的工具之一（Burdick，Langston，1977；Vinnik，1977；吳慶舉，曾融生，1998）.Zhu和 Kanamori （2000）提出了接收函數(shù)H-κ疊加技術(shù)來(lái)估計(jì)地殼厚度H和地殼平均波速比κ，該方法根據(jù)地殼參數(shù)對(duì)H值和κ值求取理論到時(shí)，并依據(jù)該理論到時(shí)所對(duì)應(yīng)接收函數(shù)的轉(zhuǎn)換波和多次波的振幅大小來(lái)確定最優(yōu)地殼厚度和波速比.該方法不需要人工拾取地震震相到時(shí)，且操作簡(jiǎn)單，被廣泛應(yīng)用于各種構(gòu)造區(qū)域的地殼厚度和波速比κ的計(jì)算中（Julià，Mejía，2004； Wanget al，2010）.然而由于地震臺(tái)站下方復(fù)雜地殼結(jié)構(gòu)的影響，當(dāng)多次波震相無(wú)法被清晰識(shí)別時(shí)，則很難得到可靠的地殼結(jié)構(gòu)（Zhu，Kanamori，2000； Liet al，2014）.

重力法也是研究地殼結(jié)構(gòu)和成分的重要工具，該方法具有較高的水平分辨率（Guo，Gao，2018； Zhaoet al，2018）.完整的布格重力異常包括莫霍界面起伏引起的莫霍重力異常和地殼內(nèi)密度分布不均勻引起的地殼重力異常（Blakely，1995；王謙身，2003），該異常不僅與地殼厚度和密度有關(guān)，還與地殼波速比κ有關(guān).為此，Lowry和Pérez-Gussinyé （2011）提出了利用接收函數(shù)與重力聯(lián)合反演地殼厚度和波速比的方法，并將其用于美國(guó)西部地區(qū)地殼厚度和波速比的確定.Shi等（2018）對(duì) Lowry 和 Pérez-Gussinyé （2011）提出的聯(lián)合估計(jì)技術(shù)進(jìn)行了簡(jiǎn)化和改進(jìn)，使其易于操作，并具有較高的效率和精度.在不考慮地?zé)嵊绊懙那闆r下，利用接收函數(shù)與完全布格重力異常的聯(lián)合約束，確定了地殼厚度和vP/vS.該方法的優(yōu)點(diǎn)在于，即使接收函數(shù)多次波不清晰，也可以很好地確定地殼厚度和波速比（Christensen，1996；Lowry，Pérez-Gussinyé，2011）.

上述的接收函數(shù)H-κ疊加方法和重力與接收函數(shù)聯(lián)合約束方法，都需要給定一個(gè)初始的地殼平均P波速度.當(dāng)?shù)貧て骄鵓波速度變化時(shí)，地殼厚度和波速比也會(huì)發(fā)生變化（Ammonet al，1990）.為了減少給定 P 波速度所產(chǎn)生的誤差，Ma和 Zhou （2007）將面波頻散引入接收函數(shù)H-κ疊加方法中，通過(guò)同時(shí)擬合瑞雷波群速度頻散和接收函數(shù)轉(zhuǎn)換波、多次波震相到時(shí)來(lái)估計(jì)地殼厚度、vP/vS和地殼平均vP.不足之處在于，如果接收函數(shù)的多次波不清晰，該方法也無(wú)法給出可信的地殼厚度和vP/vS.

本研究試圖借鑒上述接收函數(shù)與重力、接收函數(shù)與面波頻散聯(lián)合方法各自的優(yōu)勢(shì)，利用接收函數(shù)、重力和面波頻散信息聯(lián)合約束地殼厚度、地殼波速比和平均P波速度，并采用理論和真實(shí)觀測(cè)資料進(jìn)行測(cè)試，以證明該方法的可行性和有效性.

1 方法

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

接收函數(shù)波形記錄中，在直達(dá)P波之后是莫霍面的Ps轉(zhuǎn)換波和多次波（PpPs和PsPs＋PpSs），它們相對(duì)于直達(dá) P 波的到時(shí)可以表達(dá)為 （Zhu，Kanamori，2000）：

式中H為地殼厚度，vP和vS分別為P波和S波速度，p為入射波射線參數(shù).

假設(shè)vP為常數(shù)，則波速比κ＝vP/vS隨vS的變化而變化.對(duì)于給定的H和κ，相應(yīng)的轉(zhuǎn)換波和多次波振幅疊加譜為（Zhu，Kanamori，2000）：

式中S（H，κ）為H-κ疊加譜，w1，w2和w3為加權(quán)系數(shù)，r（ti） （i＝1，2，3）為轉(zhuǎn)換波和多次波的波形振幅.

本研究中設(shè)置H和κ的范圍分別為20——60 km和1.50——2.00，步長(zhǎng)分別為1 km和0.01.根據(jù)式（1）——（4）執(zhí)行掃描和疊加，得到H-κ疊加譜，再得到地殼厚度H和vP/vS的最優(yōu)估計(jì).

1.2 完全布格重力異常的似然估計(jì)

Lowry和Pérez-Gussinyé （2011）認(rèn)為布格重力異常ΔgB由莫霍面起伏引起的重力異常?gMolo、地殼密度分布不均勻引起的地殼重力異常Δgcrust和熱流引起的重力異常組成.在重力似然反演方法中引入滑動(dòng)窗口技術(shù)，在如此小的區(qū)域內(nèi)，地溫的變化通常較小，熱流引起的重力影響可以忽略不計(jì)（Shiet al，2018； Guoet al，2019），因此 ΔgB的表達(dá)式為

莫霍面重力異常與地殼厚度的關(guān)系式、地殼重力異常與地殼波速比的關(guān)系式分別如下（Lowry，Pérez-Gussinyé，2011）：

式中：ΔρMoho為莫霍面上下密度差，可以描述地殼密度相對(duì)于波速比的變化率；D為莫霍面深度，D＝H－E，其中H為地殼厚度，E為高程；D為研究區(qū)莫霍面深度的平均值， κ為研究區(qū)波速比κ的平均值；F{·}和F?1{·}分別表示傅里葉變換和反傅里葉變換；G為萬(wàn)有引力常數(shù)，f為波數(shù)，c為地殼厚度變化的校正因子.

若已知若干點(diǎn)的地殼厚度、波速比和實(shí)測(cè)布格重力異常，通過(guò)線性回歸方法可以得到式（6）和（7）中的 ΔρMoho和 ? ρcrust/?κ；將其重新代入式（5）——（7），則可正演得到理論重力異常；再將實(shí)測(cè)布格重力異常與理論布格重力異常作差求得噪聲異常.

其似然函數(shù)為

將其取對(duì)數(shù)，并令 μ和 σ2的一階導(dǎo)數(shù)為0，則

從而可得

這樣，采用似然估計(jì)法 ［ 式（12）和（13） ］ 分別計(jì)算出μ和σ2，再根據(jù)式（9）計(jì)算似然函數(shù)值.以接收函數(shù)H-κ疊加相同的范圍和步長(zhǎng)選取H和κ值，根據(jù)式（5）——（9）進(jìn)行掃描并計(jì)算似然函數(shù)值，則得到重力H-κ似然譜.

1.3 面波頻散似然估計(jì)

觀測(cè)頻散和合成頻散之間的最小二乘擬合m（H，κ）由理論面波頻散與觀測(cè)頻散差的均方根和得到.Ma和Zhou （2007）定義了一個(gè)適度函數(shù)來(lái)得到面波H-κ似然譜：

式中mmax和mmin是H-κ平面內(nèi)m（H，κ）的最大值和最小值.在計(jì)算m（H，κ）時(shí)，H值和κ值的選擇與接收函數(shù)H-κ疊加方法相同，進(jìn)行網(wǎng)格搜索時(shí)固定vP不變，vS隨κ值變化.

1.4 聯(lián)合約束

利用接收函數(shù)、重力和面波頻散數(shù)據(jù)聯(lián)合約束的流程如下：

1） 計(jì)算接收函數(shù)H-κ疊加譜，將最優(yōu)地殼厚度和vP/vS估計(jì)作為下一步重力估計(jì)的初始值.對(duì)于采用接收函數(shù)H-κ疊加方法不能得到地殼厚度和vP/vS的地震臺(tái)站，采用頻域?yàn)V波技術(shù)（Guoet al，2013）和密度界面反演技術(shù)（Oldenburg，1974）估計(jì)莫霍面深度.將高程值加上莫霍面深度，即可得地殼厚度的初始值，然后由接收函數(shù)H-κ疊加譜導(dǎo)出對(duì)應(yīng)的vP/vS初始值.

2） 計(jì)算重力估計(jì)的第i個(gè)臺(tái)站的ΔρMoho和?ρcrust/?κ參數(shù).設(shè)定滑動(dòng)窗口的大小，搜索所有地殼厚度和vP/vS以及分布在該觀測(cè)站中心窗口內(nèi)的重力異常，然后經(jīng)網(wǎng)格化后，在一個(gè)規(guī)則的網(wǎng)格中獲取地殼厚度、vP/vS和觀測(cè)的重力異常數(shù)據(jù)，將這三個(gè)網(wǎng)格數(shù)據(jù)帶入式（5）——（7），利用線性回歸算法求解 ΔρMoho和?ρcrust/?κ.

3） 將求得的 ΔρMoho和?ρcrust/?κ重新代入式（5）——（7），求出模型理論重力異常.

4） 將觀測(cè)重力異常與模型理論重力異常作差求得噪聲異常，采用似然估計(jì)法的式（12）和（13）分別計(jì)算出μ和σ2，再根據(jù)式（9）計(jì)算似然函數(shù)值.

5） 以與接收函數(shù)H-κ疊加相同的范圍和步長(zhǎng)選取H和κ值，重復(fù)步驟3）和4），形成重力反演的H-κ似然譜.

6） 將接收函數(shù)H-κ疊加譜與重力反演的H-κ似然譜點(diǎn)乘，得到接收函數(shù)和重力聯(lián)合譜.

7） 采用面波頻散似然估計(jì)方法（Ma，Zhou，2007），以接收函數(shù)H-κ疊加的范圍和步長(zhǎng)選取H和κ值，通過(guò)改變初始κ值改變vS，得到理論面波頻散；求取其與觀測(cè)頻散差的均方根，得到面波擬合的H-κ似然譜.

8） 設(shè)定vP范圍為 6.0——6.5 km/s，通過(guò)改變初始vP（步長(zhǎng)為 0.02 km/s），根據(jù)步驟 7）得到不同vP所對(duì)應(yīng)的面波H-κ似然譜，將步驟6）得到的聯(lián)合譜最優(yōu)估計(jì)值與面波H-κ似然譜相比較.當(dāng)依據(jù)前述兩個(gè)譜分別確定的最優(yōu)κ值最接近時(shí)，所對(duì)應(yīng)的vP為最優(yōu)vP.

9） 設(shè)定最終vP，分別求取接收函數(shù)H-κ疊加譜、面波H-κ似然譜和重力H-κ似然譜，依據(jù)

求取聯(lián)合譜，并拾取最終的H值和κ值.

1.5 不確定性估計(jì)

本研究將最大振幅標(biāo)準(zhǔn)誤差范圍內(nèi)的解定義為聯(lián)合估計(jì)解，參照Eaton等（2006）評(píng)估接收函數(shù)疊加結(jié)果誤差的公式，擬采用下式來(lái)確定解的不確定性：

式中ηR和ηG分別為接收函數(shù)高斯濾波器和重力似然濾波器的寬度，NR和NG分別為用于疊加的接收函數(shù)和重力觀測(cè)的數(shù)量.當(dāng)所有解在H-κ網(wǎng)格空間上趨于高斯分布時(shí)，產(chǎn)生最大振幅的解不一定位于高斯分布的中心.因此解的相關(guān)不確定性不能用解的總體均值和標(biāo)準(zhǔn)差來(lái)描述.為了計(jì)算聯(lián)合約束方法中的不確定度，我們?nèi)「咚拐龖B(tài)分布函數(shù)下15.9%——84.1%范圍內(nèi)即68.2%的范圍用于計(jì)算平均值，類(lèi)似于高斯總體分布的一個(gè)標(biāo)準(zhǔn)偏差，而最大振幅代表最優(yōu)解（Delphet al，2019）.

2 理論模型測(cè)試

2.1 理論模型及地震、重力數(shù)據(jù)合成

為檢測(cè)該方法的可行性，本文首先建立了兩種不同的地殼速度模型（圖1），并基于這兩種模型采用Hrftn96和Surf96正演模擬程序（Herrmann，2013）合成理論接收函數(shù)和瑞雷波群速度頻散（圖2）.

圖1 用于正演測(cè)試的兩個(gè)速度模型（a） 簡(jiǎn)單速度模型（模型Ⅰ）；（b） 含低速層的速度模型（模型Ⅱ）Fig.1 Two velocity models used for forward testing（a） Simple velocity model；（b） Velocity model including a low velocity layer

圖2 基于圖1 中的模型Ⅰ（a）和Ⅱ（b）合成的接收函數(shù)（上）和瑞雷波群速度頻散 U （下）Fig.2 Synthetic receiver function （upper panels） and Rayleigh wave group velocity dispersion U（lower panels） based on the models Ⅰ （a） and Ⅱ （b） shown in Fig.1

為了對(duì)理論重力異常進(jìn)行正演，我們以待反演臺(tái)站為中心構(gòu)建一個(gè)半窗口為150 km的地殼厚度和波速比模型（圖3a，b），其中數(shù)據(jù)網(wǎng)格的間距為50 km，地殼厚度為36.2——43.1 km，vP/vS范圍為1.59——2.05.本文所用的聯(lián)合估計(jì)方法假設(shè)布格重力異常是由莫霍面起伏和殼內(nèi)密度不均勻共同引起.假設(shè)其觀測(cè)面的高程為0 km，莫霍面深度等于地殼厚度.地殼地幔密度差ΔρMoho為0.5 g/cm3，地殼密度與地殼波速比之比的偏導(dǎo)數(shù)?ρcrust/?κ為0.25.圖4為基于理論模型計(jì)算得到的布格重力異常.

圖3 用于理論重力異常正演而構(gòu)建的地殼厚度（a）和波速比（b）分布模型Fig.3 The crustal thickness （a） and vP/vS ratio （b） used for synthetic gravity anomaly

圖4 基于模型（圖3）正演得到的布格重力異常（a） 莫霍面起伏引起的重力異常；（b） 地殼內(nèi)密度不均勻引起的重力異常；（c） 理論合成的布格重力異常Fig.4 Synthetic Bouguer gravity anomalies for the model in Fig.3（a） The modeled gravity anomalies associated with undulation Moho interface；（b） The modeled gravity anomalies associated with crustal heterogeneous density distribution； （c） The modeled Bouguer gravity anomalies with 5% Gaussian noise

2.2 簡(jiǎn)單地殼模型數(shù)據(jù)測(cè)試

假定由地殼和上地幔組成的簡(jiǎn)單模型（圖1a），中心臺(tái)站下方的地殼厚度和vP/vS分別為40 km 和 1.75，在莫霍面上方和下方的vP分別為 6.10 km/s和 8.15 km/s，圖2a展示了由該模型合成的中心臺(tái)站的接收函數(shù).在這個(gè)簡(jiǎn)單模型中，接收函數(shù)的多次波震相清晰.然而，在實(shí)際應(yīng)用中，多次波的震相往往難以識(shí)別，Ps波震相是唯一比較明顯的震相.

為了模擬多次波不易識(shí)別這一現(xiàn)實(shí)情況，我們假設(shè)只存在Ps波震相，即在式（4）中設(shè)置加權(quán)系數(shù)w1為1，w2和w3均為0.通過(guò)展開(kāi)接收函數(shù)H-κ疊加，得到單極值條帶的H-κ疊加譜（圖5a）；對(duì)重力異常進(jìn)行似然估計(jì)算法，得到重力H-κ似然譜（圖5d）；將重力H-κ似然譜與接收函數(shù)H-κ疊加譜相乘，得到重力與接收函數(shù)聯(lián)合約束H-κ反演譜（圖5e）；根據(jù)面波頻散似然估計(jì)方法得到面波H-κ似然譜（圖5b），根據(jù)初始P波速度的變化范圍6.0——6.5 km/s和步長(zhǎng)0.02 km/s對(duì)比聯(lián)合譜與面波頻散H-κ似然譜，兩圖極值的最大值κ的坐標(biāo)最接近時(shí)，取此時(shí)的P波速度為H-κ疊加方法的初始P波速度；將H-κ反演譜與面波頻散H-κ似然譜點(diǎn)乘，得到H-κ聯(lián)合譜（圖5g）.由此可以得出，中心臺(tái)站的地殼厚度、vP/vS比值和初始P波速度的最優(yōu)估計(jì)值分別為（40±1.62） km，（1.75±0.032）和 6.1 km/s，與理論地殼模型相同.

圖5 基于模型Ⅰ（圖1a）得到的地殼厚度H-波速比κ約束圖（圖中白線代表最優(yōu)值68%的置信區(qū)間）（a） 接收函數(shù) H-κ 疊加譜；（b，c） 初始 P 波速度為 6.1 和 6.2 km/s時(shí)的面波 H-κ 似然譜；（d） 重力 H-κ 似然譜；（e，f） 初始P波速度為6.1和6.2 km/s時(shí)，接收函數(shù)與重力聯(lián)合約束譜；（g） 接收函數(shù)、面波和重力聯(lián)合約束譜Fig.5 H-κ stacking map based on model Ⅰ（Fig.1a） where white lines delineate 68% confidence interval（a） The receiver function H-κ stacking map；（b，c） The surface wave H-κ likelihood map with initial vP＝6.1 and 6.2 km/s；（d） The gravity H-κ likelihood map； （e，f） Normalized SRSG with vP＝6.1 km/s and vP＝6.2 km/s；（g） Joint H-κ stacking map

2.3 含低速層地殼模型數(shù)據(jù)測(cè)試

假設(shè)簡(jiǎn)單地殼模型中存在一個(gè)10 km厚的殼內(nèi)低速層（圖1b），各層的密度分別為2.65，2.7，2.8和 3.3 g/cm3，P 波速度分別為 5.60，5.20，6.80和 8.15 km/s，圖2b所示為由該模型合成的中心臺(tái)站的接收函數(shù)（高斯系數(shù)為2.5）以及由該模型合成的面波群速度頻散.

在式（4）中設(shè)置w1＝0.6，w2＝0.3，w3＝0.1，中心臺(tái)站所產(chǎn)生的接收函數(shù)H-κ疊加譜如圖6a所示.由于地殼內(nèi)部存在低速層，接收函數(shù)H-κ疊加結(jié)果存在兩個(gè)極值條帶，難以估計(jì)地殼厚度和vP/vS.本研究利用重力、接收函數(shù)和面波頻散聯(lián)合約束生成聯(lián)合H-κ疊加譜（圖6d），據(jù)此給出的地殼厚度和vP/vS最優(yōu)估計(jì)值為H＝（40±0.78） km，vP/vS＝1.76±0.028，vP＝6.1 km/s，這與本文構(gòu)建的速度結(jié)構(gòu)模型相符.

圖6 基于模型Ⅱ（圖1b）得到的H-κ約束圖（圖中白線代表最優(yōu)值68%的置信區(qū)間）（a） 接收函數(shù) H-κ 疊加譜；（b） 面波 H-κ 似然譜；（c） 重力 H-κ 似然譜；（d） 初始 P 波速度為 6.1 km/s時(shí)，接收函數(shù)與重力聯(lián)合約束譜；（e） 接收函數(shù)、面波和重力聯(lián)合約束譜Fig.6 H-κ stacking map based on model Ⅱ（Fig.1b） where white lines represent 68% confidence interval（a） Receiver function H-κ stacking map；（b） Surface wave H-κ likelihood map with vP＝6.1 km/s；（c） Gravity H-κ likelihood map；（d） Normalized SRSG with vP＝6.1 km/s；（e） Joint H-κ stacking map

3 實(shí)際數(shù)據(jù)測(cè)試

我國(guó)華南地區(qū)位于秦嶺——大別造山帶以南、青藏高原以東，該地區(qū)的地質(zhì)構(gòu)造演變過(guò)程復(fù)雜，研究該地區(qū)的地殼厚度和波速比可為研究區(qū)域地殼變形和動(dòng)力學(xué)機(jī)制提供重要依據(jù).這方面已經(jīng)受到關(guān)注并開(kāi)展了一系列研究（鄧陽(yáng)凡等，2011；Zhouet al，2012；Tenget al，2013；Huanget al，2015；Shenet al，2016；Guoet al，2019）.

本文使用楊曉瑜和李永華（2021）的接收函數(shù)資料用于實(shí)際數(shù)據(jù)分析.該研究從國(guó)家測(cè)震臺(tái)網(wǎng)數(shù)據(jù)備份中心（鄭秀芬等，2009）下載了中國(guó)國(guó)家地震數(shù)字臺(tái)網(wǎng)在華南地區(qū)的336個(gè)固定臺(tái)站在2009年1月至2018年2月期間記錄的遠(yuǎn)震波形數(shù)據(jù)，并從中挑選震中距在30°——90°之間、M＞5.5且具有清晰P波初至和高信噪比的遠(yuǎn)震資料用于P波接收函數(shù)計(jì)算.面波頻散數(shù)據(jù)源自Li等（2013），該研究利用中國(guó)國(guó)家地震數(shù)字臺(tái)網(wǎng)及GEOSCOPE，K-NET和KZ-NET等計(jì)劃在東亞大陸周邊地區(qū)布設(shè)的184個(gè)臺(tái)站所記錄的區(qū)域地震波形數(shù)據(jù)，采用單臺(tái)法提取了面波群速度頻散，并構(gòu)建了周期為10——145 s的瑞雷波群速度頻散分布圖，檢測(cè)板測(cè)試結(jié)果顯示其多數(shù)周期的頻散圖分辨率為2°.重力數(shù)據(jù)選取世界重力數(shù)據(jù)庫(kù)WGM2012中華南地區(qū)的完整布格重力異常數(shù)據(jù)（Balminoet al，2012），其網(wǎng)格間距為 2″ （圖7）.

圖7 研究區(qū)布格重力異常Fig.7 The complete Bouguer gravity anomalies in study area

為證實(shí)本文研究方法的有效性，通過(guò)對(duì)比前人研究結(jié)果（Huanget al，2015；Guoet al，2019；Luoet al，2019；楊曉瑜，李永華，2021），挑選前人接收函數(shù)H-κ疊加結(jié)果差異較大的兩個(gè)臺(tái)站（HB_NZH和HB_YDU）進(jìn)行深入分析.

3.1 HB_NZH臺(tái)站結(jié)果

圖8a為本文計(jì)算得到的HB_NZH臺(tái)站的226個(gè)接收函數(shù)波形記錄，圖8b為其觀測(cè)頻散圖（Liet al，2013）.

假定地殼厚度H和κ分別在20——60 km和1.5——2.0之間變化，依據(jù)前述方法分別計(jì)算接收函數(shù)H-κ疊加譜（圖8c）、面波H-κ似然譜（圖8d）和重力H-κ似然譜（圖8e）.在接收函數(shù)H-κ疊加過(guò)程中，當(dāng)設(shè)定w1＝0.6，w2＝0.3，w3＝0.1時(shí)，臺(tái)站接收函數(shù)H-κ疊加譜呈單極值條帶（圖8c）.由此可見(jiàn)，僅依據(jù)接收函數(shù)得到的地殼厚度和波速比具有較大的誤差.

圖8f展示了本文利用改進(jìn)的聯(lián)合約束方法生成的聯(lián)合H-κ疊加譜.從該結(jié)果來(lái)看，該臺(tái)站下方的地殼厚度和波速比最優(yōu)估計(jì)分別為（36±1.05） km和1.75±0.036，平均P波速度為 6.34 km/s，這與前人（Heet al，2013；Guoet al，2019）利用接收函數(shù)H-κ疊加方法以及接收函數(shù)與重力聯(lián)合反演方法得到的地殼厚度和波速結(jié)果（Guo，Gao，2018）較為一致，但是本文利用聯(lián)合方法給出的誤差估計(jì)明顯要小.

圖8 臺(tái)站HB_NZH數(shù)據(jù)資料及計(jì)算得到的H-κ約束圖（圖中白線代表最優(yōu)值68%的置信區(qū)間）（a） HB_NZH 臺(tái)站觀測(cè)接收函數(shù)；（b） 實(shí)測(cè)瑞雷波群速度頻散 U （Li et al，2013）；（c） 接收函數(shù) H-κ 疊加譜；（d） 面波 H-κ 似然譜；（e） 重力 H-κ 似然譜；（f） 聯(lián)合約束譜Fig.8 Data and H-κ stacking map of station HB_NZH where white lines in Figs.（c）–（f）represent the 68 percent confidence interval（a） Observed receiver functions of station HB_NZH；（b） Observed dispersions U （Li et al，2013）；（c） Receiver function H-κ stacking map；（d） Surface wave H-κ likelihood map；（e） Gravity H-κ likelihood map；（f） Joint H-κ stacking map

3.2 HB_YDU臺(tái)站結(jié)果

圖9a給出了從HB_YDU臺(tái)站計(jì)算的接收函數(shù)中選取的133個(gè)接收函數(shù)，圖9b為HB_YDU 臺(tái)站的觀測(cè)頻散圖（Liet al，2013）.

與HB_NZH臺(tái)站相同，假定地殼厚度H和κ分別在20——60 km和1.5——2.0之間變化，依據(jù)前述方法分別計(jì)算了接收函數(shù)H-κ疊加譜（圖9c）、面波頻散似然譜（圖9d）和重力似然譜（圖9e）.在接收函數(shù)H-κ疊加過(guò)程中，當(dāng)設(shè)定w1＝0.6，w2＝0.3，w3＝0.1時(shí)，臺(tái)站接收函數(shù)H-κ疊加譜呈單極值條帶（圖9c）.由此可見(jiàn)，僅僅依據(jù)接收函數(shù)得到的地殼厚度和波速比具有較大的誤差.

圖9f展示了我們利用改進(jìn)的聯(lián)合約束方法生成的聯(lián)合H-κ疊加譜.從該結(jié)果來(lái)看，該臺(tái)站下方的地殼厚度和波速比最優(yōu)估計(jì)分別為（40±1.34） km和1.67±0.027，平均P波速度為6.38 km/s，這與Luo等（2019）利用接收函數(shù)與PmP走時(shí)共同約束給出的結(jié)果（H＝40.4 km，κ＝1.67，vP＝6.34 km/s）一致，但與前人利用接收函數(shù)、接收函數(shù)與重力聯(lián)合反演方法得到的地殼厚度與波速比結(jié)果相差較大.例如：Huang等（2015）利用接收函數(shù)H-κ疊加方法得到的結(jié)果約為33.5 km和1.73，楊曉瑜和李永華（2021）利用接收函數(shù)H-κ疊加方法得到的結(jié)果為32.5 km和1.75，而Guo等（2019）利用接收函數(shù)與重力聯(lián)合反演方法得到的結(jié)果為32.5 km和1.75.事實(shí)上，該臺(tái)站的接收函數(shù)多次波并不清晰，因此，依靠接收函數(shù)或接收函數(shù)和重力聯(lián)合地約束很難準(zhǔn)確約束地殼厚度和波速比.

圖9 臺(tái)站HB_YDU數(shù)據(jù)資料及計(jì)算得到的H-κ約束圖（圖中白線代表最優(yōu)值68%的置信區(qū)間）（a） HB_YDU 臺(tái)站實(shí)測(cè)接收函數(shù)；（b） 實(shí)測(cè)瑞雷波群速度頻散 U （Li et al，2013）；（c） 接收函數(shù) H-κ疊加譜；（d） 面波頻散 H-κ 似然譜；（e） 重力 H-κ 似然譜；（f） 聯(lián)合約束譜Fig.9 Data and H-κ stacking map of station HB_YDU where white lines in Figs.（c）?（f）represent 68 percent confidence interval（a） Observed receiver functions；（b） Observed dispersions U （Li et al，2013）；（c） Receiver function H-κ stacking map； （d） Surface wave H-κ likelihood map；（e） Gravity H-κ likelihood map；（f） Joint H-κ stacking map

4 討論與結(jié)論

由于接收函數(shù)對(duì)P波速度敏感，因此，利用接收函數(shù)H-κ疊加方法（Zhu，Kanamori，2000）、重力和接收函數(shù)聯(lián)合約束方法（Shiet al，2018； Lowry，Pérez-Gussinyé，2011）計(jì)算地殼厚度和波速度比時(shí)，都需要給定一個(gè)初始的地殼平均P波速度.

接收函數(shù)和面波頻散聯(lián)合約束方法（Ma，Zhou，2007）可以用于同時(shí)估計(jì)地殼厚度、vP/vS和地殼平均vP，但其缺點(diǎn)是當(dāng)接收函數(shù)多次波不清晰時(shí)，估計(jì)的地殼結(jié)構(gòu)參數(shù)存在較大誤差.本文提出了一種利用接收函數(shù)、重力數(shù)據(jù)及面波頻散聯(lián)合約束地殼厚度、波速比和平均P波速度的新方法.通過(guò)對(duì)合成數(shù)據(jù)和實(shí)際數(shù)據(jù)的測(cè)試，結(jié)合以往研究，驗(yàn)證了該方法的可行性，證明了重力和面波約束的介入可以減少H-κ疊加結(jié)果的不確定性，提高估計(jì)的精度，并得到可靠的平均P波速度.精準(zhǔn)的地殼參數(shù)可為研究地殼物質(zhì)的組成、地殼形成演化及其動(dòng)力學(xué)過(guò)程提供有利條件.此外，當(dāng)研究區(qū)存在較厚的低速沉積層時(shí)，會(huì)導(dǎo)致接收函數(shù)多次波和轉(zhuǎn)換波的延時(shí)效應(yīng)，該方法計(jì)算得到的地殼厚度較真實(shí)值大.因此，該方法仍需進(jìn)一步優(yōu)化.

感謝審稿專(zhuān)家和編輯提出的寶貴修改意見(jiàn).