青藏高原東北緣瑞利面波成像：射線理論與程函方程結果的比較

郝識杰，黃周傳，王良書，徐鳴潔，米 寧，于大勇

南京大學 地球科學與工程學院，南京 210023

1 引言

青藏高原是印度板塊和歐亞板塊碰撞形成的年輕的高原（圖1；Yin and Harrison, 2000）。對于青藏高原的隆升機制存在熱烈的討論，前人提出許多端元模型，如地殼流模型（Clark and Royden,2000; Royden et al., 2008），剛性塊體擠出模型（Tapponnier et al., 2001），巖石圈連續(xù)變形模型（Molnar and Tapponnier, 1975）等。獲取精細的巖石圈地震波速度結構，對于討論青藏高原的構造演化過程，具有重要的意義，廣泛布設的密集地震觀測臺陣為獲取更高精度的地球內(nèi)部結構提供了可能。

“中國地震科學臺陣南北地震帶北段”觀測項目（ChinArray 二期），在青藏高原東北緣及其相鄰區(qū)域布設了675臺流動地震儀。面波速度成像是利用密集臺陣數(shù)據(jù)，對巖石圈結構進行約束的有效方法。傳統(tǒng)的面波成像研究多基于射線假設，即假設地震臺接收到的面波信號是沿震源與臺站之間的大圓弧路徑傳播（Barmin et al., 2001）。在實際情況下，由于巖石圈波速結構存在橫向變化，面波的傳播方向會偏離大圓弧路徑，導致基于射線假設的方法出現(xiàn)誤差。Foster等（2014）利用類似聚束分析的方法分析了USArray臺陣記錄的地震事件的入射角度，發(fā)現(xiàn)入射角度的偏差可以達到15°，由此造成的相速度誤差可以達到4%。Lin等（2009）提出了基于程函方程的成像方法，該方法首先計算面波的相位走時場，然后基于程函方程計算相速度的大小和方向。在該方法中，面波速度的方向從地震記錄計算得來，理論上可以減少因面波傳播偏離大圓弧路徑導致的誤差。在USArray數(shù)據(jù)集，基于程函方程的瑞利波相速度成像結果和Foster等（2014）矯正入射角度之后的成像結果有很好的一致性（Jin and Gaherty, 2015）。這一結果表明基于程函方程的成像方法可以有效地修正面波傳播方向偏離大圓弧路徑導致的誤差。

在青藏高原東北緣，前人已經(jīng)開展了許多基于面波的工作。如：基于射線的面波成像（潘佳鐵等, 2017），基于程函方程的面波成像（鐘世軍等, 2017），面波反演三維S波速度結構（Li et al.,2017），面波與接收函數(shù)聯(lián)合反演三維S波速度結構（Wang et al., 2017）。對于這些基于面波的研究工作而言，面波速度成像都是重要的基礎。在這些研究中，面波速度成像的結果有相似的總體特征（青藏高原東北緣存在低速異常，鄂爾多斯塊體西部上地殼存在低速異常，中下地殼和巖石圈表現(xiàn)出高速異常），在細節(jié)上有所區(qū)別（銀川河套地塹在地殼是否存在低速異常）。由于這些研究使用了同一個數(shù)據(jù)集，所以成像結果的區(qū)別應該與使用的成像方法有關。

本研究分別使用基于射線的成像方法（Debayle and Sambridge, 2004）和基于程函方程的成像方法（Jin and Gaherty, 2015）計算青藏高原東北緣的相速度分布。通過對兩種方法原理的闡述和兩種方法計算結果的對比，討論前人研究中存在爭議的速度結構細節(jié)。

圖1 青藏高原及周邊地區(qū)地形圖（紅色矩形標識研究區(qū)域）Fig.1 Topography map of the Tibetan Plateau and adjacent region (red rectangle denotes the study region)

2 數(shù)據(jù)和方法

2.1 數(shù)據(jù)

本研究中使用的地震數(shù)據(jù)來自ChinArray二期項目，包含675臺寬頻地震儀（圖2b），平均臺間距約40 km，記錄時間從2013年9月到2015年5月，觀測的采樣頻率為100 Hz。兩種成像方法使用了相同的天然地震瑞利波數(shù)據(jù)。為了挑選較為清晰的基階面波波形，筆者選擇了震中距在10°～160°范圍內(nèi)、矩震級大于5.5、震源深度小于50 km的地震事件，共351個（圖2a）。

2.2 計算臺站之間的相位走時差

瑞利波相速度成像的第一步是利用同一個地震事件的記錄，計算兩個臺站之間的瑞利波相位走時差。為了客觀地比較兩種成像方法，使用同一套流程計算相位走時差（Jin and Gaherty, 2015），然后分別用于兩種成像方法。首先進行預處理，通過群延遲時間計算窗口函數(shù)Ws，分離出基階面波的信號。然后通過如下公式計算互相關函數(shù)C(t):

其中S1和S2分別代表兩個臺站記錄的地震波形。對計算得到的互相關函數(shù)進行窄帶通濾波，通過擬合如下公式得到相位走時差等參數(shù)（Gee and Jordan, 1992）：

等式左側(cè)表示對互相關函數(shù)用濾波器Fi進行窄帶通濾波，并用窗口函數(shù)Wc突出主要能量。等式右側(cè)Ga表示高斯函數(shù)，用來擬合的五個參數(shù)分別為A、σ、ω、tg和tp，分別表示幅值參數(shù)，高斯窗口半寬度，窄帶通濾波的中心頻率，群延遲時間和相位延遲時間。利用擬合得到的相位走時差，可以進行下一步的成像工作。因為擬合的曲線是余弦曲線，擬合時可能產(chǎn)生誤差，導致結果跳過一個或多個波長的范圍。筆者選擇間距200 km以內(nèi)的臺站對計算相位走時，這樣面波波列較短，擬合更準確。

2.3 基于射線的面波相速度成像方法

如果震源和兩個臺站近似處于同一個大圓弧路徑上，那么兩個臺站的相位走時差可以用于計算兩個臺站之間的平均相速度。在上一步計算的相位走時差中，挑選符合條件的數(shù)據(jù)，計算了臺站對之間的平均相速度，射線密度如圖3。在得到射線上的平均相速度之后，按照Debayle和Sambridge（2004）提出的面波成像流程進行相速度成像，得到二維的相速度結構?；谏渚€理論，平均速度可以表述為如下公式：

圖2 （a）使用的地震事件分布圖（b）本研究使用的青藏高原東北緣臺站分布圖（據(jù)Taylor and Yin, 2009修改）Fig.2 (a)Seismic events used in this research (b)Stations used in this research

圖3 基于射線的成像方法中不同周期射線密度分布圖Fig.3 The density of ray paths for tomography based on ray theory

其中T表示周期，c(T,r?)是面波的傳播速度，r?表示地球表面的位置矢量，p表示面波的傳播路徑，cp(T)表示路徑上的平均相速度，Lp表示路徑的長度。面波傳播路徑近似為震源和臺站之間的大圓弧路徑，這一假設成立的條件是面波速度橫向不均一性足夠小，對射線路徑造成的擾動小于想要達到的分辨率。通過阻尼線性反演，使模型正演得到的兩點之間平均相速度逼近實際觀測的臺站之間的平均相速度。

2.4 基于程函方程的面波相速度成像方法

臺陣記錄到的，地震事件激發(fā)的單一頻率的面波波場，可以由程函方程表示：

其中δτ為兩個臺站之間的相位走時差，是慢度矢量，表示兩個臺站之間的路徑。同樣地，通過阻尼線性反演，使模型正演的相位走時差逼近觀測的相位走時差。由于基于程函方程的成像方法計算的是慢度（速度）矢量，所以任意兩個臺站對之間的相位走時差數(shù)據(jù)都能參與反演，從而更充分地利用數(shù)據(jù)。

3 結果

瑞利波相速度主要受巖石圈的P波速度，S波速度和密度影響，對S波速度變化敏感性最強，所以可以用來約束巖石圈S波速度結構。總體上，周期較長的相速度對更深部速度結構變化更敏感，速度結構變化對特定周期的相速度的影響可以通過敏感核量化地表示?；贚ITHO1.0模型（Pasyanos et al., 2014）和ak135模型（Kennett et al., 1995）設置模型，地殼分為三層，莫霍面深度約為50 km，計算了20 s、30 s、60 s的瑞利波相速度敏感核，作為分析討論成像結果的依據(jù)（圖4）。從敏感核可以看出，20 s和30 s的瑞利波相速度主要反映中下地殼的結構，60 s的結果主要反映上地幔頂部的結構。

兩種方法在20 s、30 s、60 s的成像結果如圖5所示，二者速度異常特征存在很好的一致性，且與主要構造單元的分布相關性較好。在青藏高原東北緣，兩種方法在20 s、30 s、60 s的結果都表現(xiàn)為低速異常，且低速異常分布范圍較為吻合。在20 s和30 s的結果中，松潘甘孜塊體和祁連塊體所在位置存在低速異常，且兩塊低速異常體互相分開。這一細節(jié)特征在兩種方法的結果中都有體現(xiàn)，根據(jù)敏感核的計算結果，主要反映中下地殼的S波速度結構。用背景噪聲面波反演得到的三維S波速度模型中，中下地殼松潘甘孜塊體和祁連塊體的低速異常表現(xiàn)為互相分離的兩塊（Bao et al., 2013; Li et al., 2014）。在60 s的相速度圖像中，青藏高原東北緣又整體呈現(xiàn)為低速異常體。

圖4 (a)計算敏感核使用的速度模型 (b)瑞利波相速度敏感核Fig.4 (a)Velocity model used in calculation of the kernel (b) sensitivity kernel of rayleigh wave phase velocity

在鄂爾多斯塊體，相速度圖像表現(xiàn)出顯著的高速異常，這一特征在兩種方法中都有體現(xiàn)，反映了鄂爾多斯塊體作為穩(wěn)定克拉通的特征。在鄂爾多斯塊體西北部的銀川河套地塹，兩種方法的結果在20 s、30 s都表現(xiàn)出低速異常，在60 s基于程函方程的結果有低速異常，基于射線的結果射線較稀，約束效果不理想。基于程函方程的成像結果中，低速異常體的分布與銀川河套地塹的輪廓更為吻合。

4 討論

4.1 兩種成像方法的對比

分別使用基于射線的成像方法和基于程函方程的成像方法得到了青藏高原東北緣的瑞利波相速度成像結果，兩種方法使用了同一數(shù)據(jù)集計算的相位走時差?；谏渚€的成像方法只使用了與地震事件近似處于同一大圓弧的臺站對之間的相位走時差，而基于程函方程的成像方法不受這一限制，利用數(shù)據(jù)更加充分。且后者同時反演面波速度的大小和方向，可以避免因偏離大圓弧路徑傳播導致的誤差。從結果上看，兩種方法取得了相似的相速度圖像，但有一些細節(jié)區(qū)別。為了對比兩種方法計算的瑞利波相速度差異，筆者用射線成像方法的結果減去程函方程成像方法的結果（圖6）。兩種方法得到的速度結構主要在研究區(qū)邊緣及速度異常體邊緣存在差異。筆者推測這是由于基于射線的成像利用的數(shù)據(jù)有限，對異常體形態(tài)約束不如程函方程方法充分，且在研究區(qū)邊緣射線密度較小，表現(xiàn)出接近于平均速度的相速度特征。因此筆者認為對于研究區(qū)邊緣的細節(jié)結構，基于程函方程的成像結果更加可靠。

4.2 青藏高原東北緣的低速異常結構

面波速度主要受S波速度影響，S波低速異?？赡苡捎诮橘|(zhì)溫度高、力學性質(zhì)弱或包含流體等因素造成。此次研究的成像結果中，20 s和30 s瑞利波相速度在松潘甘孜塊體和祁連塊體表現(xiàn)出顯著低速異常，且兩塊低速異常體互不相連。結合相速度敏感核，這反映了該區(qū)域中下地殼存在S波速度低速異常。

地殼流模型是解釋青藏高原變形機制的一個端元模型，即青藏高原中下地殼較為軟弱，在地形高程梯度形成的壓強差作用下，由高原中心向高原邊緣發(fā)生塑性流動（Clark and Royden, 2000;Royden et al., 2008）。地殼流模型解釋了青藏高原內(nèi)部平坦的地形特征和青藏高原邊緣陡峭和平緩兩種類型的地形變化（Clark and Royden, 2000）。青藏高原通過地殼流機制進行隆升和擴張的前提之一是存在軟弱、可以流動的中下地殼。發(fā)生部分熔融可能改變中下地殼的流變性質(zhì)，導致其可以流動。部分熔融會產(chǎn)生地球物理手段可以觀測的異?，F(xiàn)象。在松潘甘孜塊體，接收函數(shù)和面波聯(lián)合反演顯示中下地殼存在較高的P波S波波速比（Ye et al., 2017），電阻率反演顯示中下地殼存在低電阻率層（Wang et al., 2014），Lg波Q值成像反映出低Q值（Zhao et al., 2013）。這些地球物理研究結果支持松潘甘孜塊體的中下地殼存在部分熔融。本研究結果中的松潘甘孜塊體相速度低速異常反映了中下地殼S波低速異常，這些S波低速異?？赡芘c中下地殼發(fā)生部分熔融有關。由于存在部分熔融，在重力作用下，松潘甘孜塊體有產(chǎn)生地殼流的條件。同時，各向異性研究也在地殼觀測到與地形梯度平行的快波方向，符合地殼流模型的預測（Huang et al., 2017）。

圖5 青藏高原東北緣瑞利波相速度分布圖Fig.5 Rayleigh wave phase velocity maps of the Northeastern Tibetan Plateau

圖6 射線成像方法與程函方程成像方法結果的對比Fig.6 Comparison between phase velocities which are calculated based on ray theory and Eikonal equation

成像結果在祁連塊體表現(xiàn)出和松潘甘孜塊體相似的相速度低速異常，但是在祁連塊體支持中下地殼存在部分熔融的證據(jù)較弱，該塊體的縱橫波速比低于松潘甘孜地體（Tian and Zhang, 2013），但電阻率高于松潘甘孜地體（Xiao et al., 2013）。祁連塊體可能發(fā)育逆沖推覆斷層，地殼增厚，成為含水流體運移的通道（Xiao et al., 2012），造成低速異常。

4.3 銀川河套地塹的巖石圈結構

銀川河套地塹位于鄂爾多斯塊體的西北部，在印度板塊和歐亞板塊碰撞之后形成（Zhang et al., 1998）。基于前文的討論，筆者認為基于程函方程的成像能更好地約束這一區(qū)域。用該方法得到的相速度分布圖在20 s、30 s和60 s，在銀川河套地塹表現(xiàn)出低速異常，且分布范圍與地塹的輪廓較為一致，表明本區(qū)域中下地殼和上地幔頂部存在S波速度低速異常。體波成像也支持上地幔存在低速異常（Huang and Zhao, 2006）。接收函數(shù)研究顯示銀川河套地塹下方巖石圈厚度（約80 km）比相鄰的鄂爾多斯塊體（約200 km）?。–hen et al., 2009）。上地幔頂部的低速異?？赡苁怯捎谲浟魅ξ镔|(zhì)上涌造成的（Tian et al.,2011）。太平洋板塊俯沖到華北克拉通下方，可能在上地幔產(chǎn)生了地幔楔（Zhao et al., 2007）。軟流圈上涌可能與板塊脫水和地幔楔中的拐角流有關（Huang and Zhao, 2006）。中下地殼存在低速異常表明軟流圈上涌的影響到達了地殼，同時表現(xiàn)出低電阻率（Dong et al., 2014）和高的縱橫波速比（Tian et al., 2011）。

5 總結

本研究使用ChinArray二期共675臺寬頻流動地震儀記錄的遠震瑞利波數(shù)據(jù)，分別使用基于射線的面波成像方法和基于程函方程的面波成像方法計算了瑞利波相速度分布圖，揭示了地殼和上地幔頂部的結構橫向變化特征。并通過兩種方法的對比，討論了前人研究中存在爭議的細部結構特征。

在青藏高原東北緣，祁連塊體和松潘甘孜塊體在20 s、30 s、60 s存在低速異常。這一速度特征在兩種方法取得的結果中都顯著存在，反映了地殼較厚且存在S波速度低速異常。在20 s、30 s的相速度分布圖中，祁連塊體和松潘甘孜塊體的低速異常相互分開，在60 s的結果中，兩塊低速異常體不再顯著分開。這一細節(jié)特征在兩種方法的計算結果中較為一致。支持兩個塊體存在不同的變形機制，松潘甘孜塊體可能存在地殼流，祁連塊體可能通過逆沖推覆斷層發(fā)生增厚。在20 s、30 s、60 s，鄂爾多斯塊體存在顯著的高速異常，反映了穩(wěn)定克拉通的特征。

在銀川河套地塹，基于程函方程的成像結果在20 s和30 s顯示出顯著的低速異常，且分布范圍與地塹的范圍較為一致。基于射線的成像也顯示出低速異常，但是低速體分布范圍與構造體吻合較差。這可能由于基于程函方程的成像方法更充分地利用了臺站之間地相位走時差數(shù)據(jù)，對較小的速度異常體約束更加準確。

致謝：感謝“中國地震科學探測臺陣數(shù)據(jù)中心”為本文提供觀測資料。本研究由自然科學基金項目（41674044；41674049）以及中國地震科學臺陣探測項目（DQJB16A0306）資助，黃周傳受南京大學“登峰人才支持計劃”資助。文中多數(shù)圖像是由GMT繪制，在此表示衷心感謝。