青藏高原東北緣基于W2度量的全波形成像

董興朋, 楊頂輝, 蒙偉娟

1 中國地震局地質研究所, 北京 100029

2 清華大學數學科學系, 北京 100084

0 引言

青藏高原東北緣位于鄂爾多斯塊體、秦祁地塊、四川盆地及其與高原內部各塊體的交匯部位(圖1),是青藏高原向歐亞大陸內部擴展的前緣過渡帶和印度—歐亞大陸碰撞作用由近南北向向北東和東方向轉換的重要區域(Zhang et al., 2004; 劉啟民等,2014).該區地質構造復雜,晚新生代構造變形顯著且強震活動頻繁,當前仍處于地殼增厚和巖石圈變形的初期階段(王椿鏞等,2016).因此青藏高原東北緣是研究高原橫向擴展及其與周圍地塊相互作用機制等地球動力學問題的關鍵區域,同時也是理解強震孕育及發生機制的理想場所.前人已經在青藏高原東北緣開展了廣泛的、多學科的地球物理研究工作.GPS觀測結果表明,該區正在遭受顯著的NE向的擠壓(Zhang et al., 2004; Gan et al., 2007).深地震反射、寬頻帶地震學、大地電磁測深以及重力異常研究揭示了地殼增厚、殼內普遍存在低速高導層等特點,加深了人們對青藏高原東北緣深部構造和變形機制的認識(王椿鏞等, 2003; 趙國澤等,2004; Pavlis et al., 2012; 李孟洋等,2022).

當前,青藏高原東北緣地表抬升和橫向擴展機制的研究仍具有爭議性.Tapponnier等(1982)認為,該區域的橫向擴展主要體現為沿青藏高原內部的大型走滑斷裂(如昆侖斷裂、阿爾金斷裂等)的東向擠出,并推測地殼與巖石圈地幔間可能存在強耦合現象.Clark和Royden(2000)基于數值模擬結果,提出高原深部地殼的弱物質因重力驅動在地質年代尺度上可發生水平向塑性流動,從而導致高原的隆升和橫向擴展,且殼幔由于低黏滯中下地殼的存在而可能導致該區域殼幔解耦或至少部分解耦.隨著地震層析成像和大地電磁測深研究的不斷深入,證實該區域下方存在大范圍殼內低波速區和高導層,地殼流模型被廣泛用來解釋高原的橫向擴展(Clark et al., 2005).該模型為青藏高原東緣的低地形梯度以及在剛性四川盆地阻擋下,龍門山斷裂帶所形成的高海拔、陡峭地形提供了合理的解釋.高分辨率深部結構成像是檢驗上述不同端元模型的基礎,特別是在研究青藏高原東北緣區域中下地殼低速層的空間分布、異常幅度及其連通性等方面.

地震全波形反演理論由Tarantola于1984年提出,最初在勘探地震學領域得到應用.近年來,得益于計算機性能的顯著提升和現代數值計算方法(如譜元法)的進步,實現計算區域乃至全球尺度的三維實際地球模型的理論地震圖已經成為可能(Komatitsch and Tromp, 2002).不同于基于高頻近似的傳統射線類成像方法,全波形成像采用數值方法求解地震波方程,不僅充分考慮了實際地震波的有限頻率特性,也包含了在傳播過程中的散射、衍射和波前愈合等多種復雜波場效應,實現了高精度地震波形數據的有效計算.該方法同時利用地震波的相位和振幅信息約束地下結構,并涵蓋射線路徑之外介質對波場擾動的影響,在同等條件下其成像結果分辨率顯著優于傳統的射線類成像方法(Virieux and Operto, 2009).Tape等(2009)首次利用伴隨成像方法,揭示了美國南加州地殼結構的顯著非均質性,與初始模型相比,局部地區速度擾動可達30%;其后,Fichtner等(2009)和Zhu等(2012)應用該方法獲得了澳大利亞和歐洲地區高分辨率殼幔速度結構圖像.一些學者也將該方法應用到中國大陸及周邊地區殼幔結構研究中,例如Chen等(2015)和Xiao等(2020)分別利用該方法獲得了東亞大陸和青藏高原及其周邊區域殼幔速度結構.對于區域尺度的研究,Dong等(2020)、Dong和Yang(2020)利用全波形成像方法,獲得了青藏高原東北緣和華北地區高分辨率巖石圈速度結構.值得注意的是,盡管上述研究均依賴數值方法求解地震波動方程,其反演中所采用的目標函數主要還是互相關走時信息.這種處理方法(也稱為擬波形成像方法,黃雪源等,2021)部分是因為走時信息與地球內部結構之間存在近似線性關系,有利于后續的穩定迭代計算;同時,由于觀測數據常受到噪聲的干擾和場地效應及儀器響應因素的影響,獲取地震波形的絕對振幅值在實踐中頗具挑戰.不過,這種處理意味著可能會犧牲一定的分辨率,因為振幅信息在結構約束方面的潛力并未得到充分利用.

傳統的L2度量目標函數具有很強的局限性,對數據的擾動和微小形變都會產生較大誤差,原因在于基于L2范數的目標函數對數據高度非凸、存在很多局部極小值點,因此降低它的非凸性是一個長期存在的難題.在地震全波形反演研究中,這種局限性常表現為難以克服周期跳躍(cycle-skipping)問題,導致反演過程容易陷入局部極小值,降低了后續地質解釋的準確性和可信度(Virieux and Operto, 2009).近年來,最優傳輸理論快速發展,由其引入的Wasserstein距離可以衡量兩個概率分布之間的相異程度.與傳統L2范數的“點對點”振幅對比不同,Wasserstein距離還考慮了觀測波形和合成波形之間的相位差異,使得該度量更加精確(Villani,2003).Engquist和Froese(2014)、Yang等(2018)證明了Wasserstein度量在處理地震信號時移和振幅擾動方面表現出良好的凸性,數值算例表明采用Wasserstein距離構建目標泛函,能為全波形反演提供更加穩定且精確的目標函數,顯著減少局部極小值,一致地收斂到全局最優解,提高了成像結果的可靠性.此外,該度量對地震信號中的噪聲不敏感,特別適合于實際地震觀測數據的波形反演.然而,Engquist和Froese(2014)、Yang等(2018)研究工作仍局限于勘探地震學領域,隨后,Wu等(2018)將該度量應用于天然地震精定位研究.Dong等(2019)首次將基于二次Wasserstein度量的全波形反演方法擴展到區域尺度地震成像,并引入多窗Wasserstein距離這一新技術,進一步壓制了實際地震記錄中噪聲對反演結果的影響,獲得了川滇地區可靠的殼幔結構圖像.

本文研究利用2009年1月至2022年12月期間位于青藏高原東北緣區域的114個寬頻帶地震臺站記錄到的50個區域地震事件的波形數據(圖1),采用最新發展的基于W2度量的全波形成像方法,獲得了該區域可靠的、高分辨率的殼幔S波速度、徑向各向異性和泊松比結構.根據成像結果,重點探討了青藏高原東北緣中下地殼流分布范圍、高原與周邊各塊體之間相互作用機制等科學問題.

1 方法和數據

1.1 基于W2度量的地震全波形成像方法

對于單分量地震波形數據,在進行重編碼和歸一化處理后,可視為一維的概率密度函數,如果地震波信號的持續時間為T0,則二次Wasserstein距離定義如下(Yang et al., 2018; Dong et al., 2019)：

(1)

式中f和g分別代表合成波形和觀測波形,T(t)為兩個波形記錄之間的最優傳輸映射,其表達式為：

T(t)=G-1(F(t)),

(2)

其中,G和F分別為合成波形f和觀測波形g的累積分布函數：

(3)

在地震波形反演中,目標泛函的Fréchet梯度通常采用伴隨狀態法進行求解,即由震源激發、正向傳播的地震波場與臺站處激發、反向傳播的伴隨波場之間做互相關得到.因此,伴隨震源推導和計算方法非常重要,因為它是計算目標函數梯度的前提.根據Yang等(2018)和Dong等(2019)的研究,伴隨震源的數學表達式如下：

°F(t))](t-G-1°F(t)),

(4)

其中Γ表示非零元素為1的上三角矩陣.

1.2 基于W2度量的全波形反演處理流程

從數學上講,全波形反演是一個帶偏微分方程(PDE)約束的優化問題,其目標是通過對初始模型進行逐步修正來最小化誤差函數(Tarantola,1984).首先,參考前人的研究成果(例如射線層析成像)或相關的地質先驗信息,構建初始模型,并將該模型進行網格剖分,結合相應地震事件的震源機制解和震源子波函數,進行地震波場的“正傳”模擬,從而生成相應的合成地震波形;其次,采用Wasserstein度量標準,對合成波形與實際觀測波形進行定量對比,從而獲得相應的伴隨震源,并在臺站處開展地震波場“反傳”模擬,然后將“正傳”和“反傳”模擬獲得的波場位移做互相關運算得到Fréchet敏感核,進而得到誤差梯度;再次,利用梯度類的優化算法(如共軛梯度法等),基于計算出的誤差梯度對當前模型進行更新;模型更新后,再次進行正演模擬,并與觀測數據進行比對,形成一個迭代的反演循環.該循環將持續進行,直至模型理論波形與實際觀測數據之間的誤差降至可接受水平或滿足其他預設的終止條件.

1.3 基于W2度量的青藏高原東北緣全波形成像

本研究一共50個區域地震事件參與反演,其矩震級(MW)介于4.7～6.7之間,該震級范圍既能確保觀測數據具備較高的信噪比,同時又避免因震級過大和破裂過程復雜而引發的點源近似失效問題.

圖1 青藏高原東北緣地質構造圖左上角子圖展示了青藏高原東北緣在中國大陸及鄰區的位置,而右上角子圖為本研究中使用的地震事件及臺站分布情況. ALB：阿拉善塊體;OB：鄂爾多斯塊體;QQB：秦祁地塊;QO：秦嶺造山帶;SGB：松潘—甘孜地塊;SCB：四川盆地; KLF：昆侖斷裂;XSHF：鮮水河斷裂;LMSF：龍門山斷裂;LRBF：龍日壩斷裂.

反演涉及4個參數：壓縮波速度VP、水平極化SH波速度VSH、垂直極化SV波速度VSV及密度ρ.對于單個選定的波形窗口,基于W2度量的誤差函數i表達式如下：

(5)

其中T0和T1分別代表所選定波形窗口的起始與結束時間.誤差函數的變分可寫成上述4個參數的Fréchet敏感核與擾動量的體積分(Zhu et al., 2015)：

+Kρδlnρ)d3x,

(6)

式中KVP,KVSV,KVSH和Kρ分別為壓縮波、SV波、SH波和密度的敏感核.本研究采用譜元法程序SES3D(Fichtner, 2009)來模擬地震波的傳播并進行Fréchet導數的數值計算.

本次反演的地震波頻帶為10～100 s,其初始模型基于FWEA18(Tao et al., 2018).波形窗口的選擇通過自動選取走時窗程序FLEXWIN實現(Maggi et al., 2009),并在該目標波形窗口內計算相應的誤差和伴隨震源.此外,針對反演過程中地震事件和臺站分布不均帶來的問題,本研究采用了一種地理加權策略(Ruan et al., 2019),其核心是通過提升稀疏地震事件(或臺站)分布區域的權重,并降低密集地震事件(或臺站)分布區域的權重,從而更有效地平衡數據覆蓋率,優化反演過程.圖2展示了本研究所用地震事件和臺站的權重分布情況.在反演時,所有選定波形窗口的誤差函數可寫成：

(7)

圖2 地震事件(a)和觀測臺站(b)權重分布圖

圖3 歸一化誤差函數隨迭代次數變化

1.4 目標模型分辨率測試

在全波形反演分辨率測試中,計算資源與評估準確性是兩個必須考慮的關鍵因素.雖然傳統的棋盤格式檢測板被廣泛用來評估成像質量和可靠性,但由于全波形成像所需的巨大計算量,使該方法存在困難,因為它的計算需求與實際地下結構反演相當(Zhu et al., 2012; Chen et al., 2015).與此同時,盡管Fichtner和Trampert(2011a,b)提出的點擴散函數測試方法(point-spread function test)可以評估最終模型的成像質量,該方法與目標模型分辨率之間的關系不夠直觀.考慮到目標函數的梯度(或稱Fréchet敏感核)可以清晰地揭示模型的更新方向(Tromp et al., 2005),采用梯度的空間分布評估模型分辨率尤為合適,特別是在當前全波形成像主要基于梯度類方法迭代更新的背景下.

在研究區中心,深度為35 km的位置,設置一個0.75°×0.75°的低SV速度區,速度擾動的異常幅度為5%(圖5a).結果顯示,盡管SV波的Fréchet敏感核分布與速度異常區存在一定偏差,但其主要特征卻得到了較好的恢復,并且VSH與VSV、VP之間的相互影響較小,表明該位置處的反演結果是可靠的(圖5b—d).為進一步驗證結果的可靠性,我們在上地幔85 km深度進行了額外的分辨率測試,設置參數與35 km深度相同.圖6展示的測試結果證明我們的反演結果在此深度同樣具有較高的可靠性.

圖4 反演前后初始模型(a)與最終模型(b)的W2誤差分布直方圖

2 結果與討論

圖6 青藏高原東北緣模型中心深度85 km處的分辨率測試(a) 真實SV波擾動; (b) P波Fréchet敏感核; (c) SH波Fréchet敏感核; (d) SV波Fréchet敏感核. 顏色棒單位為1.0×10-18 s2·m-4.

圖7 青藏高原東北緣不同深度剖面的各向同性S波速度分布子圖(a)中的A-A′和B-B′代表兩條垂向切片的位置.

Clark和Royden(2000)的地殼流模型指出,由青藏高原內部流出的地殼弱物質在遇到剛性的四川盆地時,在東北方向可能存在兩個主要的流動路徑.其一是繞過四川盆地,向東流入秦嶺造山帶;其二是從松潘—甘孜地塊內部經昆侖斷裂向北流動進入秦祁地塊,并最終可能抵達阿拉善塊體的南邊界.然而,根據我們的成像結果,秦嶺造山帶在中下地殼尺度表現出明顯的高S波速特征(圖7b和c),這似乎表明該地區的物質黏度相對較高,與地殼流模型中所預期的低黏度狀況不符.此外,Wei等(2017)基于瑞雷面波層析成像結果也揭示東秦嶺下方巖石圈具有高波速特征.因此,我們認為高原內部的中下地殼弱物質并未流入秦嶺造山帶.

圖8 沿垂向剖面A-A′的地形高程(a)、S波速度及其擾動(b)、徑向各向異性(c)和泊松比分布(d)黑色曲線代表Moho面,對應模型中VS=4.1 km·s-1深度.HV：高速異常;PRA：正徑向各向異性;NRA：負徑向各向異性.

圖8展示了沿剖面A-A′的地形高程、S波速度及其擾動、徑向各向異性以及泊松比分布.基于反演得到的模型,選取VS=4.1 km·s-1作為參考波速來估算剖面下方Moho面的深度,這是地震層析成像中常用的一種估算Moho面埋深的手段.Moho面作為殼幔結構分界面,其上下兩側S波速度存在明顯變化,S波在向下穿過Moho面時,其速度由約3.8 km·s-1激增到約4.4 km·s-1.因此,我們取它們之間的平均值4.1 km·s-1作為Moho的分界線.值得注意的是,不同學者選取的Moho面參考波速值可能存在細微差別(Bao et al., 2015).從松潘—甘孜塊體經秦祁地塊到阿拉善塊體,Moho面埋深逐漸變淺,這意味著從高原延伸到外部塊體時,地殼厚度逐漸減薄(劉啟民等,2014).在松潘—甘孜塊體,中下地殼明顯表現出低S波速特征,并且伴有顯著的正徑向各向異性(VSH>VSV)及高泊松比異常.研究表明,巖石的泊松比會隨地殼物質部分熔融程度的增加而增大(Watanabe, 1993),而中下地殼的各向異性與礦物如黑云母的晶格優勢排列方向有關(Shapiro et al., 2004),這可能意味著該區域的介質強度較弱,存在部分熔融并可能在重力驅動下發生水平向運移.此特征一直延伸至秦祁地塊內部,暗示地殼流可能已經跨過昆侖斷裂進入秦祁地塊.但在秦祁地塊中部,中下地殼則呈現較高的S波速和負徑向各向異性(VSH

圖9 沿垂向剖面B-B′的地形高程(a)、S波速度及其擾動(b)、徑向各向異性(c)和泊松比分布(d)黑色曲線代表Moho面,對應模型中VS=4.1 km·s-1深度.PRA：正徑向各向異性;NRA：負徑向各向異性.

圖9為沿剖面B-B′的地形高程、S波速度及其擾動、徑向各向異性以及泊松比分布,用于深入探討青藏高原與四川盆地間的相互作用機制.在松潘—甘孜地塊內,位于龍日壩斷裂以西的區域,Moho界面埋深約為60 km.此處的中下地殼呈現出明顯的低速、正徑向各向異性(VSH>VSV)以及高泊松比的特征,為地殼流模型提供了有力的支持.然而,在剖面東部尤其是跨越龍日壩斷裂之后,Moho界面開始明顯抬升.在巖石圈尺度,青藏高原塊體整體表現為明顯的低速特征,代表其巖石圈強度較弱,而四川盆地則為顯著的高速,反映了其作為一個剛性塊體的特性.在印度板塊持續推擠下,低強度的高原巖石圈與高強度的四川盆地相互作用,造成兩者交界區域陡峭的地形變化(張風雪等,2018).值得注意的是,四川盆地的巖石圈基底(圖9b,深度100～150 km)展現出向西延伸至龍日壩斷裂附近的趨勢.在徑向各向異性分布上(圖9c),龍日壩斷裂作為一個明顯的分界線,其西側的地殼和上地幔頂部都呈現出明顯的正徑向各向異性,而東側的中下地殼部分則只有弱的正徑向各向異性,上地幔頂部和地殼底部的徑向各向異性均為負值.這些結果表明,當地殼流從高原內部流出并越過龍日壩斷裂后,其活動顯著減弱.同時,龍日壩斷裂東側上地幔頂部和地殼底部以垂向變形為主,這可能是由揚子塊體巖石圈基底向西楔入高原內部所致.大地電磁和剪切波分裂研究顯示龍門山次級塊體可能構成一個相對獨立的地質單元(Zhao et al., 2012; Gao et al., 2019).在地殼尺度上,我們的成像結果與上述研究相吻合,表現為龍門山區域下方的S波速度同四川盆地和龍日壩西側的松潘—甘孜塊體存在差別.然而,需要指出的是,上述大地電磁研究的成像深度有限,主要局限在地殼尺度,而剪切波分裂研究獲得的各向異性結構則缺乏深度方向上的分辨率.在巖石圈尺度上,本研究揭示高速揚子塊體巖石圈基底已向西楔入高原內部并抵達龍日壩斷裂附近,因此,我們認為龍日壩斷裂是青藏高原和揚子塊體間的一條重要的構造邊界,它代表了揚子塊體向高原內部延伸的最西緣(Guo et al., 2015).

3 結論

全波形成像技術為揭示地球內部結構的高分辨率特征提供了有效手段,該方法克服了傳統射線類成像技術的諸多局限性(例如無限高頻假設、多路徑效應、分辨率等),顯著提升了對復雜地質體的解析能力.然而,當采用傳統的基于L2范數的目標函數進行處理時,全波形反演結果容易受到周期跳躍問題的影響,導致反演收斂到局部極小值.Wasserstein度量源于最優傳輸理論,它能同時考慮兩個概率分布之間的相位和振幅差異,有效地克服了周期跳躍問題,從而提高了成像結果的可靠性.

本研究采用基于W2度量的全波形成像方法,利用青藏高原東北緣的密集地震臺站和區域地震事件數據,獲得了該區域地殼和上地幔的高分辨率結構.顯著的低S波速、強正徑向各向異性和高泊松比異常暗示松潘—甘孜地塊中下地殼存在水平運動的弱物質流,該地殼流向北跨越昆侖斷裂進入秦祁地塊,并終止于秦祁地塊中部.然而,在中下地殼深度,秦嶺造山帶卻表現出顯著的高速異常,這表明高原內部中下地殼流并未延伸至此.值得注意的是,龍日壩斷裂兩側巖石圈存在明顯的結構差異,表明揚子塊體巖石圈基底已向西楔入高原內部并抵達該斷裂附近.

致謝謹以此文祝賀滕吉文先生90華誕暨從事地球物理工作70年.感謝中國地震局地球物理研究所地震科學國際數據中心(http：∥www.esdc.ac.cn)提供的觀測波形數據.

附錄A 初始模型沿A-A′剖面的深部結構

附圖1 初始模型FWEA18沿垂向剖面A-A′的地形高程(a)、S波速度及其擾動(b)、徑向各向異性(c)和泊松比分布(d)黑色曲線代表Moho面,對應模型中VS=4.1 km·s-1深度.

為說明反演前后最終模型相比初始模型的提升,本研究展示了初始模型FWEA18沿剖面A-A′下方的S波速度及其擾動、徑向各向異性和泊松比的分布圖像(附圖1),并與反演后的結果進行了詳盡的對比分析(圖8).相較于初始模型,反演后S波速度及其擾動圖案呈現更為豐富的細節特征,Moho界面表現出更高的平滑度,其深度從高原內部向邊緣逐漸抬升,這一趨勢與本區域地質構造特征相吻合.最顯著的提升體現在徑向各向異性結果上,在初始模型中,剖面A-A′下方整體表現為大面積的正值(VSH>VSV),而反演后的模型揭示了更加精細的變形特征,特別是在秦祁地塊中部地殼及上地幔頂部,觀察到顯著的負徑向各向異性(VSH