接收函數、面波與重力聯合約束地殼厚度與波速比＊

任志遠李永華強正陽石磊

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

2） 中國北京 100081 中國地震局震源物理重點實驗室

引 言

地殼厚度和波速比是描述地殼結構與物質組成的重要參數，其中：地殼厚度不僅可為地震波成像研究中地殼改正、重力補償與模擬等提供基本參數，還可為地殼形成演化及其動力學過程研究提供重要約束；地殼波速比則可為地殼物質組成的確定提供約束（Christensen，Fountain，1975）.

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

重力法也是研究地殼結構和成分的重要工具，該方法具有較高的水平分辨率（Guo，Gao，2018； Zhaoet al，2018）.完整的布格重力異常包括莫霍界面起伏引起的莫霍重力異常和地殼內密度分布不均勻引起的地殼重力異常（Blakely，1995；王謙身，2003），該異常不僅與地殼厚度和密度有關，還與地殼波速比κ有關.為此，Lowry和Pérez-Gussinyé （2011）提出了利用接收函數與重力聯合反演地殼厚度和波速比的方法，并將其用于美國西部地區地殼厚度和波速比的確定.Shi等（2018）對 Lowry 和 Pérez-Gussinyé （2011）提出的聯合估計技術進行了簡化和改進，使其易于操作，并具有較高的效率和精度.在不考慮地熱影響的情況下，利用接收函數與完全布格重力異常的聯合約束，確定了地殼厚度和vP/vS.該方法的優點在于，即使接收函數多次波不清晰，也可以很好地確定地殼厚度和波速比（Christensen，1996；Lowry，Pérez-Gussinyé，2011）.

上述的接收函數H-κ疊加方法和重力與接收函數聯合約束方法，都需要給定一個初始的地殼平均P波速度.當地殼平均P波速度變化時，地殼厚度和波速比也會發生變化（Ammonet al，1990）.為了減少給定 P 波速度所產生的誤差，Ma和 Zhou （2007）將面波頻散引入接收函數H-κ疊加方法中，通過同時擬合瑞雷波群速度頻散和接收函數轉換波、多次波震相到時來估計地殼厚度、vP/vS和地殼平均vP.不足之處在于，如果接收函數的多次波不清晰，該方法也無法給出可信的地殼厚度和vP/vS.

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

1 方法

1.1 接收函數H-κ疊加方法

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

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

假設vP為常數，則波速比κ＝vP/vS隨vS的變化而變化.對于給定的H和κ，相應的轉換波和多次波振幅疊加譜為（Zhu，Kanamori，2000）：

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

本研究中設置H和κ的范圍分別為20——60 km和1.50——2.00，步長分別為1 km和0.01.根據式（1）——（4）執行掃描和疊加，得到H-κ疊加譜，再得到地殼厚度H和vP/vS的最優估計.

1.2 完全布格重力異常的似然估計

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

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

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

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

其似然函數為

將其取對數，并令 μ和 σ2的一階導數為0，則

從而可得

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

1.3 面波頻散似然估計

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

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

1.4 聯合約束

利用接收函數、重力和面波頻散數據聯合約束的流程如下：

1） 計算接收函數H-κ疊加譜，將最優地殼厚度和vP/vS估計作為下一步重力估計的初始值.對于采用接收函數H-κ疊加方法不能得到地殼厚度和vP/vS的地震臺站，采用頻域濾波技術（Guoet al，2013）和密度界面反演技術（Oldenburg，1974）估計莫霍面深度.將高程值加上莫霍面深度，即可得地殼厚度的初始值，然后由接收函數H-κ疊加譜導出對應的vP/vS初始值.

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

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

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

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

6） 將接收函數H-κ疊加譜與重力反演的H-κ似然譜點乘，得到接收函數和重力聯合譜.

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

8） 設定vP范圍為 6.0——6.5 km/s，通過改變初始vP（步長為 0.02 km/s），根據步驟 7）得到不同vP所對應的面波H-κ似然譜，將步驟6）得到的聯合譜最優估計值與面波H-κ似然譜相比較.當依據前述兩個譜分別確定的最優κ值最接近時，所對應的vP為最優vP.

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

求取聯合譜，并拾取最終的H值和κ值.

1.5 不確定性估計

本研究將最大振幅標準誤差范圍內的解定義為聯合估計解，參照Eaton等（2006）評估接收函數疊加結果誤差的公式，擬采用下式來確定解的不確定性：

式中ηR和ηG分別為接收函數高斯濾波器和重力似然濾波器的寬度，NR和NG分別為用于疊加的接收函數和重力觀測的數量.當所有解在H-κ網格空間上趨于高斯分布時，產生最大振幅的解不一定位于高斯分布的中心.因此解的相關不確定性不能用解的總體均值和標準差來描述.為了計算聯合約束方法中的不確定度，我們取高斯正態分布函數下15.9%——84.1%范圍內即68.2%的范圍用于計算平均值，類似于高斯總體分布的一個標準偏差，而最大振幅代表最優解（Delphet al，2019）.

2 理論模型測試

2.1 理論模型及地震、重力數據合成

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

圖1 用于正演測試的兩個速度模型（a） 簡單速度模型（模型Ⅰ）；（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）合成的接收函數（上）和瑞雷波群速度頻散 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

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

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

圖4 基于模型（圖3）正演得到的布格重力異常（a） 莫霍面起伏引起的重力異常；（b） 地殼內密度不均勻引起的重力異常；（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 簡單地殼模型數據測試

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

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

圖5 基于模型Ⅰ（圖1a）得到的地殼厚度H-波速比κ約束圖（圖中白線代表最優值68%的置信區間）（a） 接收函數 H-κ 疊加譜；（b，c） 初始 P 波速度為 6.1 和 6.2 km/s時的面波 H-κ 似然譜；（d） 重力 H-κ 似然譜；（e，f） 初始P波速度為6.1和6.2 km/s時，接收函數與重力聯合約束譜；（g） 接收函數、面波和重力聯合約束譜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 含低速層地殼模型數據測試

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

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

圖6 基于模型Ⅱ（圖1b）得到的H-κ約束圖（圖中白線代表最優值68%的置信區間）（a） 接收函數 H-κ 疊加譜；（b） 面波 H-κ 似然譜；（c） 重力 H-κ 似然譜；（d） 初始 P 波速度為 6.1 km/s時，接收函數與重力聯合約束譜；（e） 接收函數、面波和重力聯合約束譜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 實際數據測試

我國華南地區位于秦嶺——大別造山帶以南、青藏高原以東，該地區的地質構造演變過程復雜，研究該地區的地殼厚度和波速比可為研究區域地殼變形和動力學機制提供重要依據.這方面已經受到關注并開展了一系列研究（鄧陽凡等，2011；Zhouet al，2012；Tenget al，2013；Huanget al，2015；Shenet al，2016；Guoet al，2019）.

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

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

為證實本文研究方法的有效性，通過對比前人研究結果（Huanget al，2015；Guoet al，2019；Luoet al，2019；楊曉瑜，李永華，2021），挑選前人接收函數H-κ疊加結果差異較大的兩個臺站（HB_NZH和HB_YDU）進行深入分析.

3.1 HB_NZH臺站結果

圖8a為本文計算得到的HB_NZH臺站的226個接收函數波形記錄，圖8b為其觀測頻散圖（Liet al，2013）.

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

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

圖8 臺站HB_NZH數據資料及計算得到的H-κ約束圖（圖中白線代表最優值68%的置信區間）（a） HB_NZH 臺站觀測接收函數；（b） 實測瑞雷波群速度頻散 U （Li et al，2013）；（c） 接收函數 H-κ 疊加譜；（d） 面波 H-κ 似然譜；（e） 重力 H-κ 似然譜；（f） 聯合約束譜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臺站結果

圖9a給出了從HB_YDU臺站計算的接收函數中選取的133個接收函數，圖9b為HB_YDU 臺站的觀測頻散圖（Liet al，2013）.

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

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

圖9 臺站HB_YDU數據資料及計算得到的H-κ約束圖（圖中白線代表最優值68%的置信區間）（a） HB_YDU 臺站實測接收函數；（b） 實測瑞雷波群速度頻散 U （Li et al，2013）；（c） 接收函數 H-κ疊加譜；（d） 面波頻散 H-κ 似然譜；（e） 重力 H-κ 似然譜；（f） 聯合約束譜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 討論與結論

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

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

感謝審稿專家和編輯提出的寶貴修改意見.