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

張炎， 鈕鳳林， 寧杰遠,4*

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

0 引言

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

對于該地區新生代火山活動，普遍認為與太平洋板塊的俯沖有一定的關系，但對板內火山活動起源的詳細機制，仍然存在著爭議.通過不同的觀測資料和研究方法，前人給出了該地區新生代火山形成原因的幾種觀點.全球和區域層析成像結果顯示(如Fukao et al.，1992；Bijwaard et al.，1998；Zhao，2004；Huang and Zhao，2006)，中國東北地區下方，地幔過渡帶中存在近水平延伸約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)認為長白山火山與滯留在地幔過渡帶內的太平洋板塊脫水造成的物質在“大地幔楔”中上涌有關；Tang等(2014)則認為太平洋滯留板塊存在“空缺”， 來自下地幔的熱物質沿著“空缺”上涌是長白山火山的深部動力學機制，該觀點也受到一些研究結果的支持(Liu et al.，2015；Zhang et al.，2016)；另外，Zou等(2008)認為俯沖板塊在地幔過渡帶內的停滯堆積，可能驅動軟流圈上涌至淺部導致減壓熔融；Zhang等(2016)認為五大連池和阿爾山火山的活動可能是地幔巖石圈的拆沉造成軟流圈物質上涌所引起的，同時指出在阿巴嘎火山群附近可能存在小型地幔柱.

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

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

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

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

1 數據與方法

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

1.1 一維模型相對走時計算

CCP疊加成像方法將轉換震相振幅按照射線路徑反向投影到其真實的空間位置，利用疊加提高信噪比來獲得間斷面的空間分布特征.首先通過選取適當的參考速度模型計算不同深度的P-S轉換波相對走時，再將觀測到的時間序列中各個采樣點轉換到相應的轉換深度，從而得到轉換震相與轉換點位置的對應關系，完成時深轉換.接收函數上各個點對應的振幅可以視為在這個轉換位置上生成的轉換波振幅.最后將相同轉換點附近的接收函數對應的轉換震相的振幅進行疊加，就可以得到該轉換點的轉換波振幅.

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

ΔtPds=tP+tS-tP0.

(1)

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

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

圖2 FM法計算相對走時示意圖(a)、(b)與(c)、(d)分別表示在一維與三維速度模型中計算轉換波相對走時的過程. 紅色五角星表示計算區域的地表中心點；綠色倒三角代表任一臺站位置；紅色與紫色圓點分別代表計算出的P波與S波走時計算區域的節點.(a)Δl與Δh分別表示走時場內水平方向與垂直方向網格間隔；黑色箭頭表示走時場外沿一維IASP91模型入射到走時場邊界的初始遠震體波.(c)黑色箭頭表示以任一臺站為虛擬源產生的S波.(b)、(d)黑色實線代表地表與任一轉換深度；藍色射線代表直達P波；灰色方形代表任一轉換點位置；紅色和紫色射線分別代表直達轉換點的P波和S波；紫色虛線代表轉換S波；紅色和紫色虛線方框分別代表離轉換點最近的P波和S波走時場網格范圍.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波走時場以及臺站和轉換點地理位置，由插值得到每個接收函數在各個轉換深度P-S轉換波相對直達P波的相對走時ΔtPds.對于每個接收函數，用射線追蹤法(一維模型中射線追蹤法不需要進行三維速度校正)得到的ΔtPds減去FM法得到的ΔtPds，獲得兩種方法轉換波相對走時之間的差值ΔtErrors.在一維速度模型中，FM法與射線追蹤法計算走時的結果本質上應該一致，因此可以通過ΔtErrors的大小來判斷FM法結果的可靠性并確定計算參數.

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

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

圖3 子區域劃分與走時場范圍示意圖藍色三角形表示臺站分布位置；紅色五角星表示中心接收點的位置.(a)黑色長方形表示研究區域范圍；紅色虛線將該區域劃分為15個子區域.(b)5個不同大小的走時場范圍俯視圖. 邊長4°的正方形代表研究區域的任一子區域.(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.

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

1.2 三維模型相對走時計算

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

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

2 結果

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

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

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

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

圖4 一維模型中P410s和P660s相對走時差計算結果(a)、(b)、(c) 分別表示走時場水平范圍為8°、12°和16°時，一維模型中基于射線追蹤法和FM法得到的P410s相對走時之差分布的統計結果; (d)、(e)、(f) 分別表示走時場水平范圍為8°、12°和16°時，一維模型中基于射線追蹤法和FM法得到的P660s相對走時之差分布的統計結果.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.

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

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

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

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

圖5 P410s和P660s轉換點位置分布圖(a)、(b)分別展示了P410s和P660s轉換點位置的分布情況. 紅色圓點表示轉換點的地理位置，背景顏色的深淺代表轉換點的密集程度，即高度為1 km，半徑1°的圓柱形疊加體內參與疊加的轉換點數量.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.

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

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

圖6 地幔過渡帶成像結果左列和右列分別展示了基于FM法和射線追蹤校正法的CCP疊加成像結果. 其中(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間斷面較為明顯的下陷區位于松遼盆地(圖6a中B區域)和長白山火山附近區域，其中松遼盆地410-km間斷面下陷大致集中在北緯41°—46°和東經122°—126°之間，呈北東-南西走向，長約900 km，寬約200 km，下陷程度約為12 km；長白山火山西南和東北方向的下陷區，程度約10 km.大興安嶺西南與東北部 (即阿巴嘎火山和阿爾山火山附近)也可以觀察到410-km間斷面在小范圍內的微弱下陷.在大興安嶺西南和中部，以及研究區域北部(圖6a中C、D、E區域)，存在幾處獨立的410-km間斷面微升區域.

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

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

3 討論

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

圖7 (a)三維S波速度模型剖面中直達P波與轉換P410s波射線路徑. 紅色和藍色實線代表在一維速度模型中通過射線追蹤分別求出的直達P波和轉換P410s波射線路徑；紅色和藍色虛線代表在三維速度模型中分別求出的直達P波和轉換P410s波射線路徑. 黑色虛線表示410-km和660-km深度位置. (b)展示了基于FM法，使用三維速度模型中的轉換點地理位置求出的P410s相對走時，與使用一維速度模型中轉換點地理位置求出的P410s相對走時之差的統計結果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法使用一維轉換點，在三維速度模型中計算轉換波相對走時會對結果造成的影響進行了討論.從前文所提到的45505條接收函數中隨機抽取2000條，首先使用FM法，沿用一維轉換點，在三維速度模型中計算P410s的相對走時.然后直接在三維速度模型中搜索轉換點地理位置(以下稱其為“三維轉換點”)，并基于三維速度模型計算P410s的相對走時.由于實際當中，一維轉換點與三維轉換點距離不會太遠，本文搜索三維轉換點的方法是，在三維模型中，以一維轉換點為中心，在邊長3°的正方形范圍內，水平間距每0.001°取一個點，假設其為三維轉換點，計算P410s的相對走時，則這些假設的三維轉換點中P410s相對走時最小的，即對應真實的三維轉換點.最后將FM法使用三維轉換點求出的P410s相對走時，與FM法使用一維轉換點求出的相對走時相減，取其絕對值為兩個轉換點相對走時之差.由于FM法使用三維轉換點計算出的相對走時更加接近實際情況，因此相對走時差越小，說明FM法使用一維轉換點計算相對走時的結果越準確.

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

圖8展示了東經126.6°，北緯41.0°的CCP疊加成像結果，藍色實線表示在計算相對走時的過程中，只使用一維速度模型而不進行三維校正的方法 (以下簡稱為1D法).從圖中可以看出，基于FM法(紅色實線)與射線追蹤校正法(綠色實線)，P410s和P660s的疊加轉換振幅接近，其中基于FM法的結果振幅略大一些，同時兩種方法均比基于1D法的振幅大.基于FM法得到的振幅與1D法相比，在整個區域內的振幅平均值比例為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振幅平均值).由此可以看出，選擇準確的參考速度模型并考慮三維路徑效應的影響對于結果的可靠性至關重要.FWEA18模型通過全波形反演，結合多種體波震相和面波信息，同時約束了P波和S波速度結構，具有較高的分辨率.我們認為基于該模型的FM法，能夠使得接收函數的CCP疊加成像結果更加可信.

圖8 CCP疊加結果對比圖Fig.8 Comparison of CCP stacking imaging

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

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

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

圖9 (a)地表地形起伏示意圖. 淺藍色表示海平面以下；紅色三角形表示長白山火山群. (b)紅色與藍色虛線分別表示對應地理位置410-km和660-km間斷面參與疊加的接收函數數量；黑色虛線的位置表示疊加數量為200條.(c)沿北緯42°的CCP疊加深度剖面. 黑色虛線分別表示410-km和660-km深度位置；紅色表示轉換波振幅大于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

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

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

4 結論

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

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