基于多頻段Pds波振幅的地幔過渡帶間斷面橫波速度結(jié)構(gòu)反演方法及應(yīng)用

薛 靜，于欣鑫，李 君

(1.成都理工大學(xué) 地球物理學(xué)院，成都 610059；2.中國地震局第二監(jiān)測中心，西安 710054)

位于地下大約410 km和660 km深度的地幔過渡帶間斷面(以下簡稱410間斷面與660間斷面)是地球內(nèi)部重要的相變界面，分別對應(yīng)為橄欖石—瓦茲利石(α→β)相變和林伍德石—鈣鈦礦/鐵方鎂石(γ→pv+mw)相變[1]。由于相變過程受溫度和壓力條件影響較大，地幔過渡帶間斷面的深度變化可用于推測過渡帶內(nèi)溫度的空間變化[2-9]；得益于高溫高壓巖石物理實驗和相關(guān)理論的發(fā)展，地震學(xué)觀測所得到的跨越間斷面的波速和密度變化可用于約束間斷面附近的物質(zhì)組分[10-15]，而研究地幔過渡帶內(nèi)溫度和物質(zhì)組分，對于探討板塊俯沖、克拉通破壞、地幔柱上升等地幔動力學(xué)過程[16]，進而理解地表構(gòu)造活動和巖石化學(xué)性質(zhì)有著重要意義[17-18]。

在地幔過渡帶中，跨越間斷面的地震波速度變化的幅值主要取決于間斷面兩側(cè)礦物相的彈性性質(zhì)以及相變礦物的體積含量。由于相變涉及的兩種礦物相可以在一定的溫壓范圍內(nèi)共存，相變界面并非表現(xiàn)為一階速度間斷面，所產(chǎn)生的地震波速變化也會發(fā)生在一定的深度范圍內(nèi)[19]，這個深度范圍稱為間斷面的銳度(sharpness)或厚度(thickness).410間斷面的厚度一般隨鐵含量和水含量的增加而增大。比如，在鐵鎂的相對含量為1∶9時，干燥條件下橄欖石和瓦茲利石可以在7 km的范圍內(nèi)共存；而在水含量為0.4%時，這一深度范圍可增加至11 km.660間斷面處的情形更為復(fù)雜：在660間斷面處，除林伍德石相變外，還存在石榴石-鈣鈦礦/鐵方鎂石的相變。相較于快速發(fā)生相變的林伍德石，具有更高鋁、硅含量的石榴石相變可在數(shù)十千米的范圍內(nèi)進行[19]。這兩種相變之間的關(guān)系與地幔過渡帶的溫度有關(guān)。當(dāng)?shù)蒯＿^渡帶溫度較低時，660間斷面以林伍德石的相變?yōu)橹?，間斷面較尖銳，因此高頻和低頻地震波都能觀測到。當(dāng)?shù)蒯＿^渡帶溫度較高時，石榴石的相變成為主導(dǎo)，石榴石和相變產(chǎn)物鈣鈦礦可以在很大的深度范圍內(nèi)共存，使得660間斷面的地震波速度在發(fā)生很小的躍變之后在較寬的深度范圍內(nèi)以較大的梯度增加[20]。

地幔巖石樣品大多來自深度小于200 km的冷卻巖石圈，所以想要獲取地幔過渡帶深度的巖石樣本從而測量其物理性質(zhì)和礦物組分十分困難，更為可行的手段是結(jié)合地震學(xué)觀測與礦物物理模型推斷過渡帶的相關(guān)性質(zhì)。產(chǎn)生于間斷面的反射或轉(zhuǎn)換震相的強弱反映了間斷面上下的波速和密度變化信息，用于探測地幔過渡帶間斷面的方法主要有PP前驅(qū)波和SS前驅(qū)波，Pds接收函數(shù)以及P’P’前驅(qū)波。來自660間斷面的信號在疊加的SS前驅(qū)波中十分清晰[21]，但SS前驅(qū)波是長周期信號，對660間斷面銳度的分辨率有限。局部區(qū)域的PP前驅(qū)波和接收函數(shù)研究能夠探測到660間斷面并非簡單的單個界面，來自660間斷面的信號可能表現(xiàn)為具有明顯不對稱性的波形或分裂成兩個波形，推測在640 km到720 km之間可能存在兩個界面[9，22-26]。短周期的P’P’前驅(qū)波能夠觀測到連續(xù)的660間斷面[27]，推測間斷面是比較尖銳的，但它在全球的疊加PP前驅(qū)波中并未出現(xiàn)，這可能與石榴石的相變有關(guān)[22]。由此可見，高頻地震波僅對尖銳的速度變化敏感，長周期地震波雖能分辨速度梯度但分辨率有限[28-29]，要研究660間斷面的厚度和速度變化，僅靠單一頻段的地震信號是非常困難的。本文利用接收函數(shù)方法分離來自地幔過渡帶間斷面的Pds轉(zhuǎn)換波(見圖1)，提取不同頻段的Pds轉(zhuǎn)換波與直達P波之間的振幅比(Pds/P)，使用馬爾科夫鏈蒙特卡洛方法同時反演跨越間斷面的橫波速度變化和間斷面的厚度。首先通過理論測試驗證該方法的可行性，繼而將其應(yīng)用于中國東部地區(qū)的實際數(shù)據(jù)并得到該地區(qū)的660間斷面厚度。結(jié)合高溫高壓礦物相變實驗和礦物物理模型計算，則可進一步約束過渡帶附近的地幔礦物體積組分。

(a) 來自深度為d的間斷面的Pds轉(zhuǎn)換波與直達P波，圖中藍色三角形表示接收臺站，實線表示P波而虛線表示S波; (b) 間斷面橫波速度結(jié)構(gòu)示意圖，h表示間斷面厚度，ΔvS表示跨越間斷面的橫波速度變化百分比圖1 Pds轉(zhuǎn)換波射線路徑及間斷面橫波速度結(jié)構(gòu)示意圖Fig.1 Schematics showing the Pds ray path and the shear velocity structure near the discontinuity

1 Pds轉(zhuǎn)換波振幅反演方法

遠震P波在傳播過程中遇到間斷面會有一部分能量轉(zhuǎn)化為SV波，稱為Pds轉(zhuǎn)換波，其中d代表間斷面深度，即(410，660). Pds波與P波之間的到時差取決于間斷面的深度和間斷面之上的地震波速結(jié)構(gòu)，二者之間的振幅比則對跨越間斷面的波速變化和間斷面厚度非常敏感[23，28]。本文使用接收函數(shù)法[30-32]提取Pds轉(zhuǎn)換波并得到多頻段Pds波振幅。

1.1 正演問題

正演問題中，需要由一維速度模型生成理論接收函數(shù)。本文中，使用Receiver Function Package[32-33]正演計算層狀介質(zhì)中的垂向與徑向響應(yīng)，將其旋轉(zhuǎn)至L與Q分量后，以L分量為母波形，分別將L分量與Q分量與之作反卷積，得到L接收函數(shù)與Q接收函數(shù)。進而，將二者濾波至同一頻段，提取Q分量接收函數(shù)上間斷面對應(yīng)的Pds轉(zhuǎn)換波，并用L分量上P波振幅歸一，得到Pds振幅(見圖2和圖3).

圖2顯示了具有不同橫波速度躍變(ΔvS)的一階間斷面(以660為例)產(chǎn)生的Pds接收函數(shù)波形及多頻段Pds/P振幅比。在3個模型中，間斷面厚度均為0 km，速度躍變分別為4.25%、6.25%和8.25%，地幔中S波和P波的平均品質(zhì)因子均設(shè)置為Qβ=225及Qα=500.以0.025～0.200 Hz頻段為例，3個模型中以P波振幅歸一的Q分量接收函數(shù)波形顯示，P660s和直達P波的振幅比分別為0.045，0.034和0.024，這表明在同一頻段，不同速度躍變的間斷面所對應(yīng)的Pds波形主要表現(xiàn)為振幅不同：速度躍變越大，振幅越大(圖2(b)).固定高通截止頻率為0.025 Hz，并使低通截止頻率在0.1～0.8 Hz之間取不同值，用巴特沃斯濾波器將接收函數(shù)濾波后，提取不同頻段的Pds/P振幅比(圖2(c)).在不同頻段，速度躍變越大的模型對應(yīng)的Pds/P振幅越大，而3個模型所對應(yīng)的Pds/P振幅比隨頻率的變化規(guī)律大致相同。

(a) 速度模型，其中紅色實線為IASP91參考模型(ΔvS=6.25%)，黑色和藍色實線分別為ΔvS=4.25%和ΔvS=8.25%的模型; (b) 3個模型對應(yīng)的接收函數(shù)波形(用P波歸一并濾波至0.025～0.200 Hz); (c) P660s/P振幅比隨頻率的變化，橫坐標(biāo)表示低通截止頻率，高通截止頻率固定為0.025 Hz圖2 具有不同橫波速度躍變的一階間斷面產(chǎn)生的P660s轉(zhuǎn)換波振幅(射線參數(shù)p=0.055 3 s/km)Fig.2 Amplitude of the P660s converted wave in models with different shear velocity jump across the 660 km discontinuity (ray parameter p=0.055 3 s/km)

圖3顯示了具有不同厚度的間斷面(同樣以660為例)產(chǎn)生的Pds接收函數(shù)波形及多頻段Pds/P振幅比。在3個模型中，間斷面處的速度躍變均為6.25%，分別發(fā)生于0 km、10 km和20 km的深度范圍，地幔中S波和P波的平均品質(zhì)因子均設(shè)置為Qβ=225及Qα=500.同樣以0.025～0.200 Hz頻段為例，3個模型中以P波振幅歸一的Q分量接收函數(shù)波形顯示，P660s和直達P波的振幅比分別為0.034、0.031和0.027，這表明在同一頻段，間斷面厚度越大，Pds轉(zhuǎn)換波振幅越小且波形更寬(圖3(b))；同時，隨著間斷面厚度的增加，Pds/P振幅比隨著頻率升高的衰減越快。

(a) 橫波速度模型，其中黑色線為IASP91參考模型(間斷面厚度h=0 km)，紅色實線和藍色實線分別為維持ΔvS=6.25%不變而間斷面厚度為10 km和20 km的模型; (b) 3個模型對應(yīng)的接收函數(shù)波形(用P波歸一并濾波至0.025～0.200 Hz); (c) P660s/P振幅比隨頻率的變化圖3 具有不同厚度的間斷面產(chǎn)生的P660s轉(zhuǎn)換波振幅(射線參數(shù)p=0.055 3 s/km)Fig.3 Amplitude of the P660s converted wave in models with different 660 thickness (ray parameter p=0.055 3 s/km)

兩組正演表明，跨越間斷面的橫波速度變化主要影響Pds/P振幅比大小，而間斷面的厚度同時影響Pds/P的振幅比大小及其隨頻率的變化規(guī)律。因此，可以利用多頻段振幅信息同時約束跨越間斷面的橫波速度變化ΔvS和間斷面厚度h.

1.2 反演問題

根據(jù)多頻段Pds/P振幅比反演跨越間斷面的橫波速度變化ΔvS和間斷面厚度h是一個非線性問題，使用馬爾科夫鏈蒙特卡洛方法(MCMC)[34]估計兩個參數(shù)的概率分布。假設(shè)觀測數(shù)據(jù)的不確定性均服從高斯分布，則似然函數(shù)滿足：

式中：di是由各頻段的Pds/P振幅比組成的數(shù)據(jù)向量，σi是相應(yīng)的不確定性，m=[ΔvS，h]為模型參數(shù)向量，g(m)則表示由模型參數(shù)得到Pds/P振幅比的非線性正演過程。假定ΔvS與h先驗分布分別服從[0，10%]和[0，30]km上的均勻分布，在410 000步的迭代過程中，每一步在兩個參數(shù)上分別增加標(biāo)準差為0.001和1、均值為0的高斯擾動以得到一組新的模型參數(shù)，根據(jù)Metropolis-Hastings采樣算法[34]接受或者拒絕新的模型參數(shù)，去掉剛開始的10 000組模型參數(shù)后，在所有接受的模型中每隔100步取一組形成最終的4 000組作為反演所得的參數(shù)分布。

1.3 理論測試

以厚度10 km、橫波速度變化為6%的間斷面產(chǎn)生的多頻段Pds/P振幅比加上5%的不確定性作為觀測數(shù)據(jù)，利用上述MCMC方法進行反演，得到的參數(shù)分布如圖4所示。反演結(jié)果顯示，間斷面厚度與橫波速度變化大致服從正態(tài)分布，間斷面的厚度為(9.69±1.59)km，而跨越間斷面的橫波速度變化為(5.96±0.20)%，均與已知的模型參數(shù)接近，并且反演所得的速度模型所對應(yīng)的Pds波振幅也在觀測值的誤差范圍內(nèi)，可見通過該反演方法能有效得到間斷面的厚度與跨越間斷面的橫波速度變化。

(a-d) 反演所得模型參數(shù)ΔvS和h的概率分布(對角線)和參數(shù)散點圖(非對角線); (e) 觀測數(shù)據(jù)(即由輸入模型正演所得的Pds/P振幅比，紅色方塊，誤差棒表示2σ不確定性)與一系列接受的模型參數(shù)預(yù)測的Pds/P振幅比(黑色虛線)圖4 反演方法理論測試Fig.4 Synthetic test

表1 反演結(jié)果與輸入?yún)?shù)的對比Table 1 Comparison between inverted and input model parameters

2 實際應(yīng)用

2.1 數(shù)據(jù)

從中國國家測震臺網(wǎng)數(shù)據(jù)備份中心[35]收集了2011-2018年間中國地震局固定臺網(wǎng)共970個固定臺站記錄到的震級6級以上、震中距30°～95°、震源深度小于100 km的地震(圖5(a))的三分量波形數(shù)據(jù)，對原始波形進行去儀器響應(yīng)轉(zhuǎn)換為位移記錄后，挑選出垂向分量P波信噪比大于3的數(shù)據(jù)(此處信噪比定義為以P波理論到時為中心的40 s時窗內(nèi)最大振幅與此之前40 s時窗內(nèi)最大振幅之比)，將E、N、Z三分量旋轉(zhuǎn)到L、Q、T分量，截取L分量上的P波段(P波之前20 s到P波之后50 s)作為母波形，將完整的Q與L分量(P波之前20 s到P波之后200 s)作為子波形，使用Wiener濾波法[31-32]提取L分量與Q分量接收函數(shù)。對提取的接收函數(shù)進行質(zhì)量控制，只保留Q分量接收函數(shù)最大幅值不超過L分量上P波幅值30%的數(shù)據(jù)作進一步分析。

2.2 接收函數(shù)疊加

選擇P660s入射點位于以(110°E，30°N)為圓心、半徑5°的圓形區(qū)域內(nèi)(圖5(b))的接收函數(shù)，按震中距疊加，疊加的Q分量接收函數(shù)顯示出清晰的來自地幔過渡帶間斷面的轉(zhuǎn)換波P410s與P660s(圖5(c)).為增大信噪比，將不同震中距的接收函數(shù)按照P660s的到時對齊進一步疊加?？紤]到不同震中距P660s波射線參數(shù)(即入射角)不同，根據(jù)P波透射系數(shù)將不同震中距的P660s振幅校正到70°并按照一維參考模型IASP91[36]計算所得P660s-P到時差對齊后疊加，得到一條接收函數(shù)(圖5(c)).

(a) 本文中所用到密集覆蓋研究區(qū)域的970個固定臺站(藍色三角)和861個地震事件(紅色圓點)分布; (b) 所選研究區(qū)域內(nèi)P660s在660間斷面的入射點分布，A、B、C為3條東西向測線，其上紅色方框表示了測線上3個小區(qū)域的中心; (c) Q分量接收函數(shù)疊加。左側(cè)為按震中距1°間隔疊加的Q分量接收函數(shù)(濾波至0.025～0.200 Hz頻段)，黑色虛線標(biāo)記了P410s與P660s的理論到時，上方為每個震中距內(nèi)疊加的接收函數(shù)數(shù)量，右側(cè)為將60°～95°范圍內(nèi)所有接收函數(shù)按照P660s理論到時對齊后疊加形成的一條接收函數(shù)圖5 臺站、事件、P660s入射點分布以及疊加的接收函數(shù)Fig.5 Map showing the distribution of stations, events, and piercing points of P660s waves

2.3 多頻段Pds/P振幅比反演

將疊加的接收函數(shù)濾波至不同頻段，測量Pds/P振幅比，并用Bootstrap方法[37]給出其不確定性，得到多頻段Pds/P振幅比觀測數(shù)據(jù)(圖6(a)).然而，研究區(qū)域內(nèi)間斷面深度的橫向變化以及間斷面上覆地幔波速結(jié)構(gòu)的不均勻性會導(dǎo)致從不同位置入射的P660s波與直達P波之間的走時差與一維速度模型的預(yù)測值不符，從而造成接收函數(shù)上P660s波形不能對齊，疊加所得的振幅小于實際振幅?；谝痪S參考速度模型IASP91[36]進行共轉(zhuǎn)換點疊加，建立了研究區(qū)內(nèi)660間斷面的深度模型(圖6(b))，其中660間斷面的深度范圍為640 km～690 km.根據(jù)每一條接收函數(shù)P660s入射點處的間斷面深度，計算其實際到時與理論到時的差，構(gòu)建相同頻率的脈沖函數(shù)，求出理想對齊情況下疊加振幅與非對齊情況下的振幅之比，作為振幅校正系數(shù)，得到校正后的Pds/P振幅比(圖6(a)).

將校正后的振幅作為觀測數(shù)據(jù)，應(yīng)用1.2節(jié)所述方法進行反演，得到了研究區(qū)域內(nèi)660間斷面上下的橫波速度差異和間斷面的厚度的后驗分布(圖7).結(jié)果顯示，中國東部地區(qū)660間斷面較為尖銳，厚度為(1.66±1.18)km，而其間斷面上下的橫波速度躍變?yōu)?6.5±0.1)%.

進一步將研究區(qū)域細分為若干半徑2°的圓形區(qū)域(圖5(b))，將P660s入射點位于各圓形區(qū)域內(nèi)的數(shù)據(jù)分別疊加，以上述方法反演各圓形區(qū)域內(nèi)跨越660間斷面的橫波速度變化和間斷面厚度，從而得到間斷面性質(zhì)的空間變化。圖8顯示了3條測線上，各個小圓形區(qū)域內(nèi)的Pds/P振幅比以及反演所得的ΔvS與h.在研究區(qū)內(nèi)不同位置，P660s/P振幅比幅值及其隨頻率的變化趨勢不盡相同，而整個區(qū)域的疊加振幅約為各小區(qū)域的平均值。由于小區(qū)域內(nèi)數(shù)據(jù)量變少，觀測誤差增大，反演所得的速度變化和間斷面厚度誤差均更大。考慮到不確定性，ΔvS的空間變化不大，在區(qū)域均值上下波動。各小區(qū)域間斷面厚度雖存在一定空間差異，但因其較大誤差，差異并不顯著。各小區(qū)域的間斷面厚度均高于整個區(qū)域疊加振幅所反演得到的間斷面厚度，這是因為從圖8(a)中可知，不同變化規(guī)律的Pds振幅一經(jīng)平均，其隨頻率的變化趨于平緩，而由前文中的正演模擬可知，變化平緩的轉(zhuǎn)換波振幅反映了較為尖銳的間斷面。

(a)整個區(qū)域內(nèi)P660s的疊加振幅(紅色方框)以及各個小圓形區(qū)域內(nèi)的P660s振幅(灰色)；(b)測線A、B、C上反演所得ΔvS，灰色條帶表示了正文圖7中用整個區(qū)域平均振幅反演所得的范圍；(c)測線A、B、C上反演所得間斷面厚度h，灰色條帶表示了正文圖7中用整個區(qū)域平均振幅反演所得的范圍圖8 研究區(qū)域為660間斷面性質(zhì)的空間變化Fig.8 Spatial variation of the 660 discontinuity property in the study area

3 討論

地幔過渡帶間斷面上下地震波速度變化，以及該間斷面是尖銳的界面還是具有一定厚度的漸變區(qū)域，是關(guān)系到地幔礦物組分以及水含量等的重要性質(zhì)。本文提出利用多頻段Pds/P振幅比，對間斷面厚度和跨越間斷面橫波速度變化進行概率反演。正演模型及理論測試表明，在存在觀測誤差的條件下，該方法能有效給出跨越間斷面的橫波速度變化和間斷面厚度的概率分布。注意到Pds轉(zhuǎn)換波對縱波速度變化和密度變化均不如對橫波速度變化敏感(圖9)，因此在反演中并未對縱波和密度的結(jié)構(gòu)進行特別設(shè)置。

黑色實線表示固定ΔvP與Δρ為5%而改變ΔvS時的透射系數(shù)；紅色實線表示固定ΔvS與Δρ為5%而改變ΔvP時的透射系數(shù)；藍色實線表示固定ΔvP與ΔvS為5%而改變Δρ時的透射系數(shù)圖9 Pds透射系數(shù)隨跨越間斷面的橫波速度躍變、縱波速度躍變以及密度躍變的變化Fig.9 Pds transmission coefficient as a function of model parameters ΔvP, ΔvS, and Δρ

以660間斷面為例，收集了中國東部地區(qū)固定臺網(wǎng)記錄的波形數(shù)據(jù)，利用接收函數(shù)方法分離了來自該區(qū)域660間斷面的Pds轉(zhuǎn)換波，并反演了660間斷面處的橫波速度變化和間斷面厚度。結(jié)果顯示，該區(qū)域660間斷面非常尖銳，在不到3 km的深度范圍內(nèi)橫波速度增加了約6.5%，與前人利用三重震相得到的結(jié)果較接近[38]，反映了該區(qū)域地幔過渡帶底部斜方輝橄巖富集的特點。盡管地幔巖石礦物組分的混合方式不同，計算的理論地球速度模型中，660間斷面上下的橫波速度差隨斜方輝橄欖巖含量的減少而減小[10-15]，進一步說明了本研究所得較大的速度躍變是由斜方輝橄巖富集引起的。另外，由于尖晶石相變發(fā)生在數(shù)十千米的深度范圍內(nèi)，660間斷面可能表現(xiàn)為一個較尖銳的速度間斷連接一個較大的速度梯度，要更準確地反演這種結(jié)構(gòu)，需要在未來的研究中加入對此更為敏感的反射波進行聯(lián)合約束。將地震學(xué)的研究結(jié)果結(jié)合高溫高壓礦物相變實驗和礦物物理模型計算，可更好地約束過渡帶附近的地幔礦物體積組分。

本文所選擇的研究區(qū)域橫向展布范圍約1 000 km，因此所反演所得的660間斷面處的橫波速度結(jié)構(gòu)反映了近千公里橫向范圍內(nèi)的平均性質(zhì)。在如此大的空間范圍內(nèi)，由間斷面地形起伏和波速異常結(jié)構(gòu)所引起的地震波振幅的聚焦和發(fā)散等效應(yīng)[39]可以認為被平均掉了。將研究區(qū)域進一步細分后，跨越間斷面的橫波速度變化并未顯示出明顯的空間差異，與整個區(qū)域的均值接近。間斷面厚度存在一定的空間差異，但可能是小區(qū)域內(nèi)的小尺度結(jié)構(gòu)影響了Pds振幅的變化規(guī)律，而這種影響在大區(qū)域范圍的數(shù)據(jù)疊加中被平均掉了。

4 結(jié)束語

接收函數(shù)方法是地幔間斷面結(jié)構(gòu)研究中常用的一種地震學(xué)方法，通過接收函數(shù)提取的轉(zhuǎn)換波和多次波的走時信息用于約束間斷面的深度。在大部分一維參考模型中，地幔過渡帶間斷面被表示為尖銳的一階間斷面。寬頻帶地震觀測資料的增加使得提取多頻段的Pds轉(zhuǎn)換波振幅信息約束間斷面的厚度和速度變化等復(fù)雜信息成為可能。本文提出利用多頻段的Pds/P振幅比同時反演跨越間斷面的速度變化和間斷面的厚度，為地幔過渡帶間斷面性質(zhì)的研究提供了可行的思路。在今后的工作中，我們將對大量實際觀測資料進行深入分析和處理，將研究方法進一步拓展到間斷面的多次反射波，以更好地約束間斷面附近的波速和密度結(jié)構(gòu)，探索地幔過渡帶的物質(zhì)組分和動力學(xué)過程。