基于程函方程與三維速度模型的中國東北地區(qū)地幔過渡帶接收函數(shù)研究

張炎， 鈕鳳林， 寧杰遠(yuǎn),4*

1 北京大學(xué)地球與空間科學(xué)學(xué)院， 北京 100871 2 美國萊斯大學(xué)地球、環(huán)境與行星科學(xué)系， 休斯頓 77005 3 中國石油大學(xué)石油資源與勘探國家重點實驗室和非常規(guī)天然氣研究所， 北京 102249 4 河北紅山巨厚沉積與地震災(zāi)害國家野外科學(xué)觀測研究站， 北京 100871

0 引言

中國東北地區(qū)位于歐亞板塊東緣，中朝塊體和西伯利亞地臺之間，西北部與中亞造山帶相連，西南部與華北克拉通相接，東臨日本海.受構(gòu)造運動及俯沖板塊的影響，該地區(qū)地質(zhì)歷史復(fù)雜，經(jīng)歷了蒙古—鄂霍次克洋的張開與閉合，古亞洲洋的關(guān)閉，阿穆爾板塊與西伯利亞板塊的碰撞，太平洋板塊俯沖等地質(zhì)歷史事件，以盆地-山脈相間的形式為主，形成了三個主要的組成部分：中部的松遼盆地，西部的大興安嶺以及東南部的長白山山脈(杜曉娟等，2009；Zhou and Wilde，2013).早白堊時期，東北地區(qū)經(jīng)歷了大規(guī)模的地殼拉張運動，開始形成大量的斷陷盆地，并伴隨火山活動.晚白堊紀(jì)以來，該地區(qū)發(fā)育了大量的火山活動(Liu et al.，2001；Chen et al.，2007).中國東北地區(qū)的新生代火山巖主要是堿性玄武巖，火山活動以零星偶發(fā)為主(Fan and Hooper，1991；Liu et al.，2001)，廣泛分布于盆地邊緣與山脈中，包括長白山火山群、五大連池火山群、阿爾山火山群和阿巴嘎火山群等(劉嘉麒，1999)(圖1a).研究該區(qū)域的地幔結(jié)構(gòu)，可為該地區(qū)上下地幔之間的動力學(xué)關(guān)系，俯沖板塊的現(xiàn)狀和火山起源等科學(xué)問題的研究提供重要線索.

對于該地區(qū)新生代火山活動，普遍認(rèn)為與太平洋板塊的俯沖有一定的關(guān)系，但對板內(nèi)火山活動起源的詳細(xì)機(jī)制，仍然存在著爭議.通過不同的觀測資料和研究方法，前人給出了該地區(qū)新生代火山形成原因的幾種觀點.全球和區(qū)域?qū)游龀上窠Y(jié)果顯示(如Fukao et al.，1992；Bijwaard et al.，1998；Zhao，2004；Huang and Zhao，2006)，中國東北地區(qū)下方，地幔過渡帶中存在近水平延伸約1500 km的地震波高速異常體，被解釋為在地幔過渡帶中滯留的俯沖板塊；很多研究者(如Lei and Zhao, 2005，2006; Huang and Zhao, 2006；Zhao et al.，2007，2009；Duan et al.，2009；Li and Van Der Hilst，2010；Wei et al.，2012，2015；Lei et al.，2013)認(rèn)為長白山火山與滯留在地幔過渡帶內(nèi)的太平洋板塊脫水造成的物質(zhì)在“大地幔楔”中上涌有關(guān)；Tang等(2014)則認(rèn)為太平洋滯留板塊存在“空缺”， 來自下地幔的熱物質(zhì)沿著“空缺”上涌是長白山火山的深部動力學(xué)機(jī)制，該觀點也受到一些研究結(jié)果的支持(Liu et al.，2015；Zhang et al.，2016)；另外，Zou等(2008)認(rèn)為俯沖板塊在地幔過渡帶內(nèi)的停滯堆積，可能驅(qū)動軟流圈上涌至淺部導(dǎo)致減壓熔融；Zhang等(2016)認(rèn)為五大連池和阿爾山火山的活動可能是地幔巖石圈的拆沉造成軟流圈物質(zhì)上涌所引起的，同時指出在阿巴嘎火山群附近可能存在小型地幔柱.

地幔內(nèi)部在全球尺度存在兩個地震波速度和密度發(fā)生快速躍變的間斷面，其全球平均深度約為410 km和660 km，因此被稱為410-km和660-km間斷面.兩個間斷面間的地幔被稱為是上下地幔的地幔過渡層或地幔過渡帶.通常認(rèn)為410-km間斷面是由地幔中礦物橄欖石到瓦茲利石，660-km間斷面是由林伍德石到布里奇曼石與鎂方鐵石的高溫高壓相變所致(Katsura and Ito，1989；Helffrich and Wood，2001；周春銀等，2010).由于與410-km和660-km兩個地幔間斷面相關(guān)的相變分別具有正的和負(fù)的克拉伯龍斜率(Ito and Takahashi，1989)，在溫度低于周圍地幔的區(qū)域，410-km間斷面深度將會變淺而660-km間斷面深度則會變深，形成增厚的地幔過渡帶；在溫度高于周圍地幔的區(qū)域，情況則正好相反.由于地幔物質(zhì)相變對溫度的依賴性，可以通過測量這兩個間斷面深度起伏及地幔過渡帶厚度(410-km和660-km間斷面的深度差)的橫向變化，約束地幔過渡帶內(nèi)的低溫(如板塊俯沖區(qū)域)與高溫(如熱物質(zhì)上涌區(qū)域)結(jié)構(gòu)的分布.

地震學(xué)探測是認(rèn)識地球深部結(jié)構(gòu)、物質(zhì)組成和動力學(xué)過程的重要手段.對于地球內(nèi)部間斷面結(jié)構(gòu)的研究，主要依靠間斷面處體波的折射、反射以及震相轉(zhuǎn)換等信息，其中，接收函數(shù)是獲取間斷面深度最有效的方法之一.前人在中國東北地區(qū)的接收函數(shù)研究工作中(如Ai et al.，2003, 2008；Li and Yuan，2003；Liu et al.，2015；Tian et al.，2016；Zhang et al.，2016；楊凡等，2021)，由于使用的方法和參考速度模型不同，導(dǎo)致各項研究獲得的間斷面深度與過渡帶厚度結(jié)果存在一定差異.例如Li和Yuan(2003)與Ai 等(2003)發(fā)現(xiàn)在長白山火山下方660-km間斷面有20～30 km的下陷，而在其西側(cè)下陷程度明顯減??；Liu等(2015)得到了相似的結(jié)果，同時在東北地區(qū)北部得到了大范圍660-km間斷面略微下陷區(qū)和長白山火山群以西660-km間斷面局部抬升的成像結(jié)果；Zhang等(2016)同樣發(fā)現(xiàn)了長白山火山下方和西側(cè)660-km間斷面下陷和局部抬升的結(jié)構(gòu)特征，但在北部沒有發(fā)現(xiàn)大范圍的660-km間斷面下陷跡象.在進(jìn)行410-km和660-km間斷面接收函數(shù)成像時，通常采用共轉(zhuǎn)換點(Common Conversion Point，CCP)疊加方法(Dueker and Sheehan，1997；Zhu，2000；Gilbert et al.，2003).為了獲得高質(zhì)量的地震圖像，則需要一個相對準(zhǔn)確的三維速度參考模型，采用一種既高效又準(zhǔn)確的方法計算地震波在三維模型中的走時.目前，大多數(shù)研究采用一維射線路徑走時校正算法，即首先利用一個一維參考速度模型進(jìn)行射線追蹤，計算一維射線路徑，然后計算三維速度模型在一維路徑中走時擾動的總和來獲取近似的三維走時(Ai et al.，2008； Liu et al., 2015； Zhang et al., 2016).在后文中，我們將簡稱該近似三維走時計算方法為射線追蹤校正法.當(dāng)參考模型中存在較強(qiáng)速度異常時，三維射線路徑的影響就會變得十分明顯，從而導(dǎo)致射線追蹤校正法計算走時的準(zhǔn)確度大為降低.同時，射線追蹤校正法需要對每個接收函數(shù)進(jìn)行射線追蹤與走時校正計算，計算效率較低.

本文采用了一種利用快速行進(jìn)法(Fast Marching Method，F(xiàn)M法)求解程函方程的方法，以提高計算三維模型中地震波傳播時間的準(zhǔn)確度和效率(Rawlinson and Sambridge，2004a，2004b；De Kool et al.，2006；Guan and Niu，2018).FM法是采用有限差分法求解程函方程，追蹤走時變化(Sethian and Popovici，1999).Guan和Niu(2018)優(yōu)化了利用FM法計算P-S轉(zhuǎn)換波相對走時的參數(shù)設(shè)置，并通過對一維與三維合成模型的數(shù)值模擬計算，驗證了該方法的高效性.該方法考慮了三維路徑效應(yīng)的影響，提高了走時計算的準(zhǔn)確度；在計算過程中先求出臺站下方目標(biāo)區(qū)域內(nèi)P波與S波走時場，再根據(jù)臺站與轉(zhuǎn)換點的地理位置進(jìn)行插值運算，求取P-S轉(zhuǎn)換波相對走時，從而提高了走時計算效率.Cheng等(2016)在Kirchhoff疊前偏移成像中，使用了FM法計算轉(zhuǎn)換波相對走時，取得了較好的結(jié)果.該方法橫向分辨率高，對復(fù)雜地下結(jié)構(gòu)成像較好，但是受臺站分布與計算效率較低的影響，不太適用于臺站分布不均的較大尺度區(qū)域研究.相比之下，本文所用的基于FM法的CCP疊加方法兼顧了計算效率與準(zhǔn)確性，更適合于大尺度區(qū)域地幔過渡帶研究.

本文中，我們首先利用國際合作項目NECESSArray以及中國國家地震臺網(wǎng)記錄到的實際波形數(shù)據(jù)，對該方法在大尺度研究區(qū)域中計算的可行性進(jìn)行了評估，然后利用Tao等(2018)全波形反演獲得的東亞地區(qū)高分辨率三維速度模型作為參考模型，對中國東北地區(qū)地幔過渡帶結(jié)構(gòu)進(jìn)行了接收函數(shù)成像，最后依據(jù)本文獲取的成像結(jié)果，探討了太平洋板塊在地幔過渡帶的存在狀態(tài)，并對板內(nèi)火山的起源進(jìn)行了嘗試性解釋.

1 數(shù)據(jù)與方法

我們采用了Liu等(2015)的接收函數(shù)數(shù)據(jù)，該數(shù)據(jù)是通過國際合作在中國東北地區(qū)布設(shè)的臨時地震儀臺站組成的大型地震臺陣(NECESSArray)和中國國家地震臺網(wǎng)，共255個臺站觀測的遠(yuǎn)震記錄計算得到(圖1a)，遠(yuǎn)震地震事件數(shù)為788個(圖1b)，挑選出接收函數(shù)總計45505條.臺站較為均勻地分布在東經(jīng)115°—135°，北緯40°—49°的研究區(qū)域內(nèi)，觀測時間為2009年9月至2011年8月.三維速度參考模型采用的是Tao等(2018)通過三維全波形反演方法獲得的東亞地區(qū)P波和S波速度模型(FWEA18).

1.1 一維模型相對走時計算

CCP疊加成像方法將轉(zhuǎn)換震相振幅按照射線路徑反向投影到其真實的空間位置，利用疊加提高信噪比來獲得間斷面的空間分布特征.首先通過選取適當(dāng)?shù)膮⒖妓俣饶Ｐ陀嬎悴煌疃鹊腜-S轉(zhuǎn)換波相對走時，再將觀測到的時間序列中各個采樣點轉(zhuǎn)換到相應(yīng)的轉(zhuǎn)換深度，從而得到轉(zhuǎn)換震相與轉(zhuǎn)換點位置的對應(yīng)關(guān)系，完成時深轉(zhuǎn)換.接收函數(shù)上各個點對應(yīng)的振幅可以視為在這個轉(zhuǎn)換位置上生成的轉(zhuǎn)換波振幅.最后將相同轉(zhuǎn)換點附近的接收函數(shù)對應(yīng)的轉(zhuǎn)換震相的振幅進(jìn)行疊加，就可以得到該轉(zhuǎn)換點的轉(zhuǎn)換波振幅.

參考速度模型的選取會對CCP疊加結(jié)果產(chǎn)生較大的影響.參考速度模型更接近真實的地幔，計算得到的轉(zhuǎn)換波相對走時則更準(zhǔn)確，成像結(jié)果也更可靠.對于一條接收函數(shù)，將從震源到臺站的直達(dá)P波的走時記作tP0，P-S轉(zhuǎn)換點的深度記作d，從震源到該轉(zhuǎn)換點的P波和從該轉(zhuǎn)換點到臺站的S波的走時分別記作tP和tS，則該P-S轉(zhuǎn)換波相對直達(dá)P波的相對走時ΔtPds可以表示如下：

ΔtPds=tP+tS-tP0.

(1)

本研究采用了De Kool等(2006)通過快速行進(jìn)法計算走時的程序(程序包名稱為fm3d)，該程序采用了高效且無條件穩(wěn)定的快速行進(jìn)算法數(shù)值求解球面程函方程(Sethian，1996；Popovici and Sethian，2002).以一維IASP91速度模型為例，對于任一地震事件E到任一臺站S的ΔtPds計算過程如下.

首先，選擇一個范圍合適的計算區(qū)域，并選擇合適的水平向與垂直向間距劃分出走時場的三維網(wǎng)格點.對于地震事件E，運行fm3d代碼，計算遠(yuǎn)震P波的走時場，并將三維網(wǎng)格中每個網(wǎng)格點對應(yīng)的走時儲存起來(圖2a).直達(dá)P波走時tP0則可以由地表離臺站S位置最近的4個網(wǎng)格點儲存的走時，通過插值得到.對于轉(zhuǎn)換深度d，通過一維射線追蹤得到該轉(zhuǎn)換點的地理位置，在P波走時場內(nèi)找到對應(yīng)轉(zhuǎn)換點附近8個網(wǎng)格點儲存的走時信息，插值得到tP(圖2b).然后，通過相同的計算過程，可以求出S波的走時場，并將從震源到臺站的S波走時記為tS0，從震源到轉(zhuǎn)換深度d的S波走時記為tSd.通過(tS0-tSd)可以計算得到約等于公式(1)里的tS的值(圖2b).

圖2 FM法計算相對走時示意圖(a)、(b)與(c)、(d)分別表示在一維與三維速度模型中計算轉(zhuǎn)換波相對走時的過程. 紅色五角星表示計算區(qū)域的地表中心點；綠色倒三角代表任一臺站位置；紅色與紫色圓點分別代表計算出的P波與S波走時計算區(qū)域的節(jié)點.(a)Δl與Δh分別表示走時場內(nèi)水平方向與垂直方向網(wǎng)格間隔；黑色箭頭表示走時場外沿一維IASP91模型入射到走時場邊界的初始遠(yuǎn)震體波.(c)黑色箭頭表示以任一臺站為虛擬源產(chǎn)生的S波.(b)、(d)黑色實線代表地表與任一轉(zhuǎn)換深度；藍(lán)色射線代表直達(dá)P波；灰色方形代表任一轉(zhuǎn)換點位置；紅色和紫色射線分別代表直達(dá)轉(zhuǎn)換點的P波和S波；紫色虛線代表轉(zhuǎn)換S波；紅色和紫色虛線方框分別代表離轉(zhuǎn)換點最近的P波和S波走時場網(wǎng)格范圍.Fig.2 Illustration of the FM method for calculating relative traveltimes(a), (b) and (c), (d) show the process of calculating the relative traveltimes of the converted waves in the 1-D and 3-D velocity models, respectively. The red pentagrams represent the central receiving points on the surface; the green inverted triangles represent the location of any station; the red and purple dots represent the calculated P- and S-wave traveltime field nodes, respectively. (a) Δl and Δh denote the horizontal and vertical intervals in the traveltime field, respectively; the black arrows indicate the initial teleseismic waves incident to the traveltime field boundary along the IASP91 model outside the field. (c) The black arrows represent the S-wave generated by any station as a virtual source. (b), (d) The black solid lines represent the surface and any conversion depth; the blue lines indicate the direct P-wave; the gray square means the location of any conversion point; the red and purple lines represent the direct P- and S-wave at the conversion point, respectively; the purple dashed line represents the converted S-wave; the red and purple dashed boxes represent the nearest P- and S-wave traveltime field grids to the conversion point, respectively.

最后，通過P波和S波走時場以及臺站和轉(zhuǎn)換點地理位置，由插值得到每個接收函數(shù)在各個轉(zhuǎn)換深度P-S轉(zhuǎn)換波相對直達(dá)P波的相對走時ΔtPds.對于每個接收函數(shù)，用射線追蹤法(一維模型中射線追蹤法不需要進(jìn)行三維速度校正)得到的ΔtPds減去FM法得到的ΔtPds，獲得兩種方法轉(zhuǎn)換波相對走時之間的差值ΔtErrors.在一維速度模型中，F(xiàn)M法與射線追蹤法計算走時的結(jié)果本質(zhì)上應(yīng)該一致，因此可以通過ΔtErrors的大小來判斷FM法結(jié)果的可靠性并確定計算參數(shù).

對于走時場計算區(qū)域范圍的參數(shù)設(shè)置，主要包括三維體外部大小(水平方向L與垂直方向H)和內(nèi)部網(wǎng)格點間距(水平方向Δl與垂直方向Δh)兩個部分(圖3c和圖2a).Guan 和Niu (2018)對網(wǎng)格點內(nèi)部間距，用了6組不同的參數(shù)進(jìn)行了測試，發(fā)現(xiàn)在水平方向上間隔0.05°(Δl)，深度間隔5 km(Δh)的情況下，求出的ΔtErrors相對較小.經(jīng)過確認(rèn)其正確性，本研究的計算中，對網(wǎng)格點內(nèi)部間距沿用了其參數(shù)設(shè)置.

由于Guan和Niu(2018)的研究區(qū)域較小，沒有對計算區(qū)域外部范圍進(jìn)行討論，而是直接選取水平方向邊長為10°，垂直方向800 km的走時場范圍進(jìn)行計算.然而，當(dāng)研究區(qū)域水平范圍較大時，如果選取一個可以包含整個研究區(qū)域在內(nèi)的走時場進(jìn)行計算，會非常耗時.因此，本研究中將整個研究區(qū)域分為若干子區(qū)域進(jìn)行計算(圖3a).

圖3 子區(qū)域劃分與走時場范圍示意圖藍(lán)色三角形表示臺站分布位置；紅色五角星表示中心接收點的位置.(a)黑色長方形表示研究區(qū)域范圍；紅色虛線將該區(qū)域劃分為15個子區(qū)域.(b)5個不同大小的走時場范圍俯視圖. 邊長4°的正方形代表研究區(qū)域的任一子區(qū)域.(c) 任一走時場范圍立體圖. L和H分別表示走時場水平與垂直范圍.Fig.3 Illustration of sub-area division and traveltime field rangeBlue triangles indicate the distribution of stations; red pentagrams denote the location of the central receiving point. (a) Black rectangle shows the study area; the red dashed lines divide the area into 15 sub-areas. (b) The top view of 5 different sizes of the traveltime field. The square with 4° represent any sub-area. (c) 3-D view of any traveltime field. L and H indicate the horizontal and vertical ranges of the traveltime field, respectively.

在選擇走時場外部范圍時，考慮到要將地幔過渡帶包括在內(nèi)，因此垂直方向的深度定為800 km(圖3c中H).如果選擇的計算區(qū)域水平范圍太小，可能導(dǎo)致射線路徑從側(cè)面入射，在660-km之上傳入；反之，如果水平范圍太大，會因為三維速度模型空間范圍不夠而產(chǎn)生假象；同時，水平范圍太大或太小也可能導(dǎo)致原本震中距在30°～90°的接收函數(shù)超出震中距閾值范圍，受到其他震相的影響.本文選取了8°、10°、12°、14°和16°(圖3c中L對應(yīng)的值)共5個不同大小的走時場計算區(qū)域水平范圍進(jìn)行測試(圖3b)，在每個計算區(qū)域內(nèi)，通過FM法和射線追蹤法分別計算了45505條接收函數(shù)P410s(410-km間斷面P-S轉(zhuǎn)換波)和P660s(660-km間斷面P-S轉(zhuǎn)換波)的ΔtErrors.

1.2 三維模型相對走時計算

本研究中三維速度參考模型選取的是Tao等(2018)全波形反演得到的東亞地區(qū)高分辨率P和S波三維速度模型，通過插值的方法，將三維模型中的地震波波速賦值到計算區(qū)內(nèi)各個網(wǎng)格點上.在FM法的計算過程中，可以將三維速度模型看作特殊的只有一層的一維模型.因此，tP的計算與1.1節(jié)中介紹的方法相同，而tS的計算過程存在一定的差異.

首先選擇一個范圍合適的三維計算區(qū)域，從震源到該區(qū)域邊界采用一維IASP91模型計算.在三維計算區(qū)域內(nèi)，采用1.1節(jié)中介紹的方法，計算出tP0與tP(圖2a、d).然后，對于任一臺站，將該臺站視為一個虛擬源，在整個走時場三維體內(nèi)傳播S波，計算出S波的走時場，并將三維網(wǎng)格中每個網(wǎng)格點的走時儲存起來(圖2c).對于轉(zhuǎn)換深度d，在S波走時場內(nèi)找到對應(yīng)轉(zhuǎn)換點附近8個網(wǎng)格點儲存的走時信息，插值得到從臺站到該轉(zhuǎn)換點的S波走時，由互易原理可知該走時即為tS.最后，通過P波和S波走時場以及臺站和轉(zhuǎn)換點地理位置，由插值得到每個接收函數(shù)在各個轉(zhuǎn)換深度，P-S轉(zhuǎn)換波相對直達(dá)P波的相對走時ΔtPds(圖2d).為了提高計算速度，接收函數(shù)的轉(zhuǎn)換點選取為一維射線追蹤得到的地理位置.

2 結(jié)果

2.1 基于一維模型的ΔtErrors和CCP疊加成像

利用45505條接收函數(shù)中的P410s和P660s轉(zhuǎn)換震相，在3個不同水平范圍的計算區(qū)域中(8°、12°和16°)，通過射線追蹤法和FM法得到了相對走時差ΔtErrors隨震中距的分布情況，并通過柱狀圖進(jìn)行了統(tǒng)計(圖4).

從圖4a—c可以看出，當(dāng)計算區(qū)域水平范圍為8°時，4729個接收函數(shù)P410s相對走時差大于0.1 s或小于-0.1 s.當(dāng)水平范圍大于或等于12°時，幾乎全部的相對走時差都在±0.1 s之間，僅有2個(12°)大于0.1 s.

圖4d—f的結(jié)果顯示，P660s相對走時差在各個計算區(qū)域水平范圍內(nèi)，大于0.1 s或小于-0.1 s的數(shù)量都占有一定的比例.其中，計算區(qū)域水平范圍為8°時計算出的相對走時差，約有14.7%超出了±0.1 s的范圍，最大值達(dá)到了約7.5 s.隨著計算區(qū)域水平范圍的增加，超出±0.1 s的比例逐漸降低，水平范圍12°時，比例最低，只有約2.5%.而隨著走時場水平范圍繼續(xù)增加，比例又開始有所增長.同時，隨著水平范圍的增加，計算每個走時場的運算時間以及儲存走時場數(shù)據(jù)占用的空間也隨之增長.

圖4 一維模型中P410s和P660s相對走時差計算結(jié)果(a)、(b)、(c) 分別表示走時場水平范圍為8°、12°和16°時，一維模型中基于射線追蹤法和FM法得到的P410s相對走時之差分布的統(tǒng)計結(jié)果; (d)、(e)、(f) 分別表示走時場水平范圍為8°、12°和16°時，一維模型中基于射線追蹤法和FM法得到的P660s相對走時之差分布的統(tǒng)計結(jié)果.Fig.4 Relative traveltime differences of P410s and P660s in the 1-D model(a), (b), (c) represent the statistical results about the distribution of the relative traveltime difference of P410s obtained by ray tracing method and FM method in the 1-D model when the horizontal range of traveltime field is 8°, 12° and 16°, respectively, while (d), (e), (f) represent the corresponding statistical results of P660s.

因此，本研究在一維模型中，走時場計算區(qū)域最終選擇水平范圍12°，深度范圍800 km，并使用與Liu等(2015)相同的CCP疊加方法進(jìn)行成像.通過對比發(fā)現(xiàn)，基于FM法的CCP疊加成像結(jié)果與基于射線追蹤法得到的結(jié)果在相同緯度基本一致.

2.2 基于三維模型的ΔtErrors和CCP疊加成像

基于一維模型的ΔtErrors結(jié)果，本文在基于三維模型的計算中，選擇計算區(qū)域水平范圍12°的走時場進(jìn)行計算.數(shù)值檢驗結(jié)果表明，基于三維速度模型的相對走時差ΔtErrors的分布情況與一維速度模型時類似.P410s相對走時差在各個計算區(qū)域水平范圍內(nèi)，±0.5 s之間的數(shù)量占據(jù)的比例都超過了90%，其中使用12°的水平范圍計算結(jié)果占據(jù)的比例最高，達(dá)到了約91.9%.同樣，在12°的計算區(qū)域水平范圍時，得出的P660s相對走時差在±0.5 s之間的比例最高，達(dá)到約89.3%.

圖5展示了410-km和660-km間斷面的P-S轉(zhuǎn)換點的地理位置分布，同時通過色標(biāo)表示出共享相同轉(zhuǎn)換點附近參與疊加的接收函數(shù)數(shù)量.CCP疊加過程中劃分的網(wǎng)格大小，深度間隔為1 km，水平方向0.1°×0.1°；圖5所示轉(zhuǎn)換點密度，是高度為1 km，半徑1°的圓柱形疊加體內(nèi)參與疊加的轉(zhuǎn)換點數(shù)量(具體方法參見Liu et al.，2015).可以看出研究區(qū)域中絕大部分區(qū)域疊加的接收函數(shù)數(shù)量都超過了200條，靠近中部的大部分區(qū)域普遍在1000條以上，最多達(dá)到2000條左右，表明本研究最終獲得的成像結(jié)果擁有較高的可信度.

圖5 P410s和P660s轉(zhuǎn)換點位置分布圖(a)、(b)分別展示了P410s和P660s轉(zhuǎn)換點位置的分布情況. 紅色圓點表示轉(zhuǎn)換點的地理位置，背景顏色的深淺代表轉(zhuǎn)換點的密集程度，即高度為1 km，半徑1°的圓柱形疊加體內(nèi)參與疊加的轉(zhuǎn)換點數(shù)量.Fig.5 Distribution of conversion point locations for P410s and P660s(a) and (b) show the distribution of conversion point locations for P410s and P660s, respectively. The red dots indicate the geographical location of the conversion points, and the shade of the background color represents the density of the conversion points, which is the number of conversion points involved in stacking within a cylinder with 1 km height and 1° radius.

最終的成像區(qū)域，即是圖5展示的接收函數(shù)疊加數(shù)量大于200條的區(qū)域.基于FM法計算出的P-S轉(zhuǎn)換波相對走時，將接收函數(shù)進(jìn)行CCP疊加，獲得了該區(qū)域410-km和660-km間斷面的起伏形態(tài)，并通過兩個間斷面深度之差，得到了過渡帶厚度的橫向變化情況(圖6a、b、c).同時，本文展示了使用相同接收函數(shù)和三維速度模型的情況下，基于射線追蹤校正法的CCP疊加成像結(jié)果(圖6d、e、f).東北角(圖6a、d中A區(qū)域)，410-km間斷面轉(zhuǎn)換波波峰并不明顯，接收函數(shù)數(shù)量也不是特別多，疊加結(jié)果不太穩(wěn)定，推測該區(qū)域間斷面深度的結(jié)果并不可靠，因此不在本文討論范圍內(nèi).從圖6中可以看出，通過對比，使用相同的三維參考模型的情況下，兩種方法得到的間斷面起伏形態(tài)與過渡帶厚度的成像結(jié)果基本一致，但是由于三維路徑效應(yīng)的影響，導(dǎo)致在范圍、深度與厚度方面存在一定的差別.

FM法得到的410-km間斷面深度在402～426 km之間變化(圖6a)，660-km間斷面深度在646～702 km之間變化(圖6b)，過渡帶厚度在232～282 km之間變化(圖6c).660-km間斷面最明顯的下陷區(qū)位于東經(jīng)128°—131°之間，靠近長白山火山群，呈近南北向展布，集中于北緯40°—46°之間(圖6b中F區(qū)域)，該下陷區(qū)域的位置和走向與Wadati-Benioff帶600 km深度的位置對應(yīng)較好.在研究區(qū)域北部和西南部(圖6b中G、H區(qū)域)，同樣可以觀察到一定范圍的660-km間斷面下陷區(qū)，其中G區(qū)域的660-km間斷面下陷程度較深，而H區(qū)域的660-km間斷面下陷較淺.下陷區(qū)F西側(cè)約200 km處，660-km間斷面存在顯著的抬升，范圍在東經(jīng)122°—126°，北緯40°—45°之間(圖6b中I區(qū)域).另外，在東部(圖6b中J區(qū)域)與阿爾山火山附近區(qū)域，660-km間斷面也有明顯抬升.

圖6 地幔過渡帶成像結(jié)果左列和右列分別展示了基于FM法和射線追蹤校正法的CCP疊加成像結(jié)果. 其中(a)、(d)與(b)、(e)分別表示410-km間斷面和 6 60-km間斷面深度的橫向變化；(c)、(f)表示過渡帶厚度的橫向變化. 圖中注釋與圖1a相同.Fig.6 Mantle transition zone imaging resultsThe left and right columns show the imaging results of CCP stacking based on FM method and ray tracing correction method, respectively. (a), (d) and (b), (e) show the lateral variation of the 410-km and 660-km discontinuity, respectively; (c) and (f) show the lateral variation of the transition zone thickness. The figures are annotated as in Fig.1a.

410-km間斷面較為明顯的下陷區(qū)位于松遼盆地(圖6a中B區(qū)域)和長白山火山附近區(qū)域，其中松遼盆地410-km間斷面下陷大致集中在北緯41°—46°和東經(jīng)122°—126°之間，呈北東-南西走向，長約900 km，寬約200 km，下陷程度約為12 km；長白山火山西南和東北方向的下陷區(qū)，程度約10 km.大興安嶺西南與東北部 (即阿巴嘎火山和阿爾山火山附近)也可以觀察到410-km間斷面在小范圍內(nèi)的微弱下陷.在大興安嶺西南和中部，以及研究區(qū)域北部(圖6a中C、D、E區(qū)域)，存在幾處獨立的410-km間斷面微升區(qū)域.

射線追蹤校正法得到的410-km間斷面深度在402～428 km之間變化(圖6d)，660-km間斷面深度在650～704 km之間變化(圖6e)，過渡帶厚度在232～284 km之間變化(圖6f).與FM法的結(jié)果相比，射線追蹤校正法得到的660-km間斷面起伏形態(tài)，在長白山火山下方的下陷程度更深一些；而在北部，下陷的范圍與程度略小一些，同時整體上660-km間斷面的抬升程度更小(圖6e).410-km間斷面的整體抬升相對于FM法略小，而松遼盆地410-km間斷面的下陷更明顯(圖6d).這種差異主要是FM法計算轉(zhuǎn)換波相對走時的時候，直接在三維速度模型中進(jìn)行計算，減小了三維路徑效應(yīng)的影響，使得相對走時的計算結(jié)果更加接近真實情況.

與間斷面深度的對比結(jié)果類似，從圖中可以看出，兩種方法得到的過渡帶厚度的橫向變化基本一致，而在增厚與減薄的程度上略有不同(圖6c、f).由于410-km間斷面深度整體起伏不大，過渡帶厚度受660-km間斷面深度影響較大，在長白山火山周圍、研究區(qū)域北部、大興安嶺中部與西南部等區(qū)域，分別發(fā)現(xiàn)了增厚的過渡帶；在松遼盆地周邊、阿爾山火山周圍和研究區(qū)域東部則出現(xiàn)了減薄的過渡帶.

3 討論

如前文所述，在使用FM法計算三維速度模型中轉(zhuǎn)換波相對走時的時候，沿用了在一維速度模型中通過射線追蹤法得到的轉(zhuǎn)換點地理位置(以下稱其為“一維轉(zhuǎn)換點”).圖7a繪制出了在三維模型中，射線追蹤校正法所使用的一維射線路徑(紅色與藍(lán)色實線)和FM法直接在三維速度模型中使用的射線路徑(紅色與藍(lán)色虛線).由于路徑上存在速度異常區(qū)域，雖然轉(zhuǎn)換點的地理位置保持不變，但是FM法得到的直達(dá)P波與P-S轉(zhuǎn)換波的射線路徑與一維速度模型中得到的射線路徑相比，存在一定的差異，差異的大小與射線路徑經(jīng)過的速度異常程度和范圍有關(guān).因此，總體而言，在相同的三維速度模型中，射線路徑上速度異常區(qū)的存在，對最終得到的轉(zhuǎn)換波相對走時存在明顯的影響.

圖7 (a)三維S波速度模型剖面中直達(dá)P波與轉(zhuǎn)換P410s波射線路徑. 紅色和藍(lán)色實線代表在一維速度模型中通過射線追蹤分別求出的直達(dá)P波和轉(zhuǎn)換P410s波射線路徑；紅色和藍(lán)色虛線代表在三維速度模型中分別求出的直達(dá)P波和轉(zhuǎn)換P410s波射線路徑. 黑色虛線表示410-km和660-km深度位置. (b)展示了基于FM法，使用三維速度模型中的轉(zhuǎn)換點地理位置求出的P410s相對走時，與使用一維速度模型中轉(zhuǎn)換點地理位置求出的P410s相對走時之差的統(tǒng)計結(jié)果Fig.7 (a) Direct P-wave and converted P410s ray paths in the 3-D S-wave velocity model profile. The red and blue solid lines represent the direct P-wave and converted P410s ray paths derived by ray tracing in the 1-D velocity model; the red and blue dashed lines represent the direct P-wave and converted P410s ray paths in the 3-D velocity model. The black dashed lines indicate the positions of 410-km and 660-km. (b) shows the statistical results of the difference between the relative traveltime of P410s based on FM method by using the geographic location of the converted points in the 3-D and 1-D velocity models

本文對FM法使用一維轉(zhuǎn)換點，在三維速度模型中計算轉(zhuǎn)換波相對走時會對結(jié)果造成的影響進(jìn)行了討論.從前文所提到的45505條接收函數(shù)中隨機(jī)抽取2000條，首先使用FM法，沿用一維轉(zhuǎn)換點，在三維速度模型中計算P410s的相對走時.然后直接在三維速度模型中搜索轉(zhuǎn)換點地理位置(以下稱其為“三維轉(zhuǎn)換點”)，并基于三維速度模型計算P410s的相對走時.由于實際當(dāng)中，一維轉(zhuǎn)換點與三維轉(zhuǎn)換點距離不會太遠(yuǎn)，本文搜索三維轉(zhuǎn)換點的方法是，在三維模型中，以一維轉(zhuǎn)換點為中心，在邊長3°的正方形范圍內(nèi)，水平間距每0.001°取一個點，假設(shè)其為三維轉(zhuǎn)換點，計算P410s的相對走時，則這些假設(shè)的三維轉(zhuǎn)換點中P410s相對走時最小的，即對應(yīng)真實的三維轉(zhuǎn)換點.最后將FM法使用三維轉(zhuǎn)換點求出的P410s相對走時，與FM法使用一維轉(zhuǎn)換點求出的相對走時相減，取其絕對值為兩個轉(zhuǎn)換點相對走時之差.由于FM法使用三維轉(zhuǎn)換點計算出的相對走時更加接近實際情況，因此相對走時差越小，說明FM法使用一維轉(zhuǎn)換點計算相對走時的結(jié)果越準(zhǔn)確.

本文隨機(jī)抽取的2000條接收函數(shù)，基于FM法計算相對走時，使用一維轉(zhuǎn)換點與使用三維轉(zhuǎn)換點得到結(jié)果的差別，大部分小于0.1 s(約86%)，最大值不超過0.6 s，平均差異值約為0.05 s(圖7b)，即使用一維轉(zhuǎn)換點得到的相對走時與真實情況十分接近，因此不需要再計算三維速度模型中轉(zhuǎn)換點的位置，而是沿用一維轉(zhuǎn)換點即可獲得接近真實的結(jié)果.

圖8展示了東經(jīng)126.6°，北緯41.0°的CCP疊加成像結(jié)果，藍(lán)色實線表示在計算相對走時的過程中，只使用一維速度模型而不進(jìn)行三維校正的方法 (以下簡稱為1D法).從圖中可以看出，基于FM法(紅色實線)與射線追蹤校正法(綠色實線)，P410s和P660s的疊加轉(zhuǎn)換振幅接近，其中基于FM法的結(jié)果振幅略大一些，同時兩種方法均比基于1D法的振幅大.基于FM法得到的振幅與1D法相比，在整個區(qū)域內(nèi)的振幅平均值比例為AP410s_FM/AP410s_1D=1.368，AP660s_FM/AP660s_1D=1.127(AP410s_FM表示基于FM法的P410s振幅平均值，AP410s_1D表示基于1D法的P410s振幅平均值；AP660s_FM表示基于FM法的P660s振幅平均值，AP660s_1D表示基于1D法的P660s振幅平均值).由此可以看出，選擇準(zhǔn)確的參考速度模型并考慮三維路徑效應(yīng)的影響對于結(jié)果的可靠性至關(guān)重要.FWEA18模型通過全波形反演，結(jié)合多種體波震相和面波信息，同時約束了P波和S波速度結(jié)構(gòu)，具有較高的分辨率.我們認(rèn)為基于該模型的FM法，能夠使得接收函數(shù)的CCP疊加成像結(jié)果更加可信.

圖8 CCP疊加結(jié)果對比圖Fig.8 Comparison of CCP stacking imaging

結(jié)果顯示，在整個研究區(qū)域中，410-km和660-km間斷面的平均深度與過渡帶的平均厚度分別為414 km、667 km和253 km，與Liu 等(2015)和Zhang 等(2016)的研究結(jié)果相比，410-km間斷面深度的結(jié)果基本一致，660-km間斷面深度和過渡帶厚度存在一定的差異，比前人結(jié)果減小了約10 km.圖6a中顯示的410-km間斷面抬升與下陷的區(qū)域，與Tao 等(2018)的三維模型S波速度的異常區(qū)有較好的對應(yīng).其中C、D和E區(qū)域410-km間斷面的抬升與Zhang 等(2016)的成像結(jié)果基本一致，由于這三個區(qū)域離俯沖板塊的作用區(qū)域較遠(yuǎn)，同時結(jié)合新生代火山分布與全波形反演結(jié)果，本文支持Zhang等(2016)的觀點，認(rèn)為該區(qū)域410-km間斷面的抬升，主要的原因可能是由于巖石圈拆沉，致使低溫物質(zhì)下沉，從而造成410-km附近溫度降低，相變位置變淺.

結(jié)合Tao等(2018)的全波形反演結(jié)果，長白山—五大連池連線一帶660-km間斷面下陷的形態(tài)與范圍和太平洋俯沖板塊在過渡帶底部的滯留有較好的對應(yīng)關(guān)系，其中南部F區(qū)域?qū)?yīng)于日本海俯沖帶在地幔過渡帶的滯留區(qū)；北部G區(qū)域?qū)?yīng)于千島俯沖帶在地幔過渡帶的滯留區(qū)；H區(qū)域的660-km間斷面下陷則可能與巖石圈物質(zhì)拆沉有關(guān).660-km間斷面最明顯的下陷區(qū)域(圖6b中F區(qū)域)，下陷幅度約為20～40 km，與Wadati-Benioff帶的600 km輪廓線形態(tài)基本一致，位于深源地震帶的延伸方向上.這一結(jié)果與Liu 等(2015)和Zhang 等(2016)的研究結(jié)果相似，但是在范圍與深度方面有一定的差異.更深的660-km間斷面相變位置通常歸因于較低的過渡帶溫度，根據(jù)實驗室測試結(jié)果，660-km間斷面的克拉伯龍斜率變化范圍在-1.3～-2.8 MPa/K之間(Ito and Takahashi，1989；Fei et al.，2004)，以-2.0 MPa/K為例計算，取圖6b中F區(qū)域660-km間斷面平均下陷深度30 km進(jìn)行計算，則該區(qū)域660-km間斷面的下陷對應(yīng)大約-500 K的溫度變化.Tao 等(2018)全波形反演結(jié)果顯示在該區(qū)域存在P波高速異常區(qū)，與本文中下陷范圍有較好的對應(yīng)，異常值約為2%，考慮到地震波速對溫度的敏感度約為-0.43%/100K(Cammarano et al.，2003)，則該高速異常對應(yīng)的溫度變化約為-465 K，與本文估算的結(jié)果大致一致.圖6b中F區(qū)域的下陷范圍向西延伸到約東經(jīng)128°，橫向范圍約300 km，大致說明了俯沖的太平洋板塊在地幔過渡帶的滯留范圍，與前人(如Zhao，2004；Huang and Zhao，2006)顯示的滯留板塊近水平向西延伸約1500 km的結(jié)果并不一致.

沿北緯42°由西向東，在深度剖面中(圖9c)可以較為清晰地識別出P410s和P660s轉(zhuǎn)換波，從而判斷出兩個間斷面大致的深度分布.410-km間斷面在東經(jīng)116.4°附近有小幅度抬升(圖9c中紅色橢圓位置)，在松遼盆地與長白山火山下方區(qū)域明顯下陷(圖9c中藍(lán)色橢圓位置)；660-km間斷面在大興安嶺與長白山火山附近區(qū)域有明顯下陷(圖9c中綠色橢圓位置)，在松遼盆地下方區(qū)域顯著抬升(圖9c中黃色橢圓位置).深度剖面的結(jié)果與前文所描述的間斷面深度異常情況有較好的對應(yīng)關(guān)系，同時更加直觀地顯示出了部分異常所處的位置、范圍以及間斷面起伏的程度.

圖9 (a)地表地形起伏示意圖. 淺藍(lán)色表示海平面以下；紅色三角形表示長白山火山群. (b)紅色與藍(lán)色虛線分別表示對應(yīng)地理位置410-km和660-km間斷面參與疊加的接收函數(shù)數(shù)量；黑色虛線的位置表示疊加數(shù)量為200條.(c)沿北緯42°的CCP疊加深度剖面. 黑色虛線分別表示410-km和660-km深度位置；紅色表示轉(zhuǎn)換波振幅大于0的部分Fig.9 (a) Topographic relief map. The light blue color indicates below sea level; the red triangle indicates the Changbaishan Volcanic zones. (b) The red and blue dashed lines represent the number of receiver function involved at the corresponding geographic locations of 410-km and 660-km discontinuity, respectively; the position of the black dashed line indicates the number of receiver function involved in stacking is 200. (c) The depth profile along 42°N. The black dashed lines indicate the position of 410 km and 660 km, respectively; the red color represents the part of the conversion wave with amplitude greater than 0

研究區(qū)域內(nèi)660-km間斷面明顯抬升的位置分別位于圖6b中的I、J和阿爾山火山附近區(qū)域，與全波形層析成像結(jié)果S波的低速異常區(qū)有較好的對應(yīng)，說明這三處位置660-km間斷面附近可能有高溫異常存在.其中最顯著的抬升區(qū)I，位于俯沖板塊前端西側(cè)，其南北向狹長的形狀與F區(qū)域類似，空間上近似平行，抬升幅度約為10 km，推測可能與太平洋俯沖板塊的破壞或下沉有關(guān)，從而導(dǎo)致局部下地幔高溫物質(zhì)上涌.同時長白山火山的巖漿活動，也有可能是熱物質(zhì)沿俯沖板塊前端破裂缺口上涌形成的.上涌的熱物質(zhì)受到淺層結(jié)構(gòu)的影響，在410-km間斷面附近形成不同的熱物質(zhì)上升區(qū)，造成了松遼盆地及長白山火山410-km間斷面附近溫度上升、深度加深，為該區(qū)域高溫物質(zhì)來源與長白山火山起源問題提供了約束.圖6b中與F和I區(qū)域空間上近似平行的J區(qū)域660-km間斷面抬升，則可能是由于板塊后端破壞所造成的高溫物質(zhì)上涌，致使660-km間斷面附近溫度升高，相變位置變淺.阿爾山火山附近區(qū)域660-km間斷面的抬升，也可能是受到深部熱流作用而形成的.由于阿爾山火山與J區(qū)域均位于成像區(qū)邊緣，其圖像需要補(bǔ)充更多數(shù)據(jù)，結(jié)合更大范圍的成像結(jié)果進(jìn)行確認(rèn).

本文主要討論了溫度與地幔過渡帶速度間斷面深度之間的關(guān)系.410-km和660-km間斷面的深度除了受溫度影響以外，還可能受到水含量和化學(xué)成分等因素的影響，從而使得實際情況變得更加復(fù)雜.許多礦物物理學(xué)實驗與地球動力學(xué)模型表明，古老冷卻的海洋板塊中含水礦物能夠?qū)⑺畮氲蒯＿^渡帶內(nèi)，停滯于過渡帶內(nèi)的海洋板塊可以將水釋放到周圍的地幔中，從而對410-km和660-km間斷面深度產(chǎn)生很大影響(如Thompson，1992；Helffrich，2000；Higo et al.，2001；Harte，2010).一些研究表明，水的存在對地幔過渡帶內(nèi)間斷面深度的影響有著類似于低溫的作用，尤其是過渡帶底部，水的存在對660-km間斷面下陷的影響甚至比溫度降低還要大，俯沖板塊中2.0 wt%的水含量可能會引起660-km間斷面下陷約15 km (Litasov et al.，2005；Helffrich，2000；Pearson et al.，2014).也有實驗證據(jù)表明化學(xué)成分對間斷面深度的影響，比如，含鐵量對410-km間斷面深度的影響與含水量相似，這意味著在富含鐵元素的區(qū)域，410-km間斷面深度將變得更淺(Katsura et al.，2004).我們會在未來的研究工作中對這部分內(nèi)容進(jìn)行更加詳細(xì)的分析與討論.

4 結(jié)論

本文使用了中國東北地區(qū)NECESSArray記錄到的遠(yuǎn)震事件接收函數(shù)，采用快速行進(jìn)法計算P-S轉(zhuǎn)換波相對走時，通過與射線追蹤校正法得到的相對走時以及成像結(jié)果的對比，說明了FM法在大尺度區(qū)域研究中的可行性.基于東亞地區(qū)高精度三維速度模型，通過FM法得到了中國東北地區(qū)地幔過渡帶的三維結(jié)構(gòu).成像結(jié)果顯示，在長白山—五大連池連線一帶，有明顯的660-km間斷面南北向下陷和過渡帶增厚區(qū)，東西方向延伸約300 km.同時在該區(qū)域東西兩側(cè)，發(fā)現(xiàn)了在空間上與下陷區(qū)近似平行延伸的660-km間斷面抬升與過渡帶減薄區(qū)，其中660-km間斷面最顯著抬升區(qū)位于長白山西側(cè)約200 km.表明日本海俯沖帶在地幔過渡帶中近水平滯留的延伸距離有限，而且可能造成地幔熱物質(zhì)沿其所造成的破裂缺口上涌，提供了長白山等新生代火山活動的巖漿源.而松遼盆地西側(cè)660-km間斷面的下陷，反映的更有可能是該地區(qū)巖石圈物質(zhì)拆沉的結(jié)果.研究區(qū)域內(nèi)410-km間斷面顯示出的所有下陷區(qū)，與地表新生代火山活動區(qū)/拉伸區(qū)有很好的對應(yīng)關(guān)系，表明這些地表構(gòu)造與深部熱物質(zhì)上涌有關(guān)，同時也表明地表構(gòu)造可能對深部熱物質(zhì)上涌通道有塑造作用.

致謝感謝NECESSArray國際合作項目和中國地震局地球物理研究所“中國地震臺網(wǎng)數(shù)據(jù)管理中心”提供的波形數(shù)據(jù).感謝中國地震局地球物理研究所柳正博士提供的接收函數(shù)數(shù)據(jù).感謝匿名審稿人提出的很好的建設(shè)性意見.本研究工作得到北京大學(xué)高性能計算校級公共平臺支持.