青藏高原東北緣構造應力場及動力學特征

孟文, 郭祥云, 李永華, 韓立波, 張重遠

1 中國地震局地球物理研究所, 北京 100081 2 中國地震局震源物理重點實驗室, 北京 100081 3 中國地質科學院地質力學研究所, 北京 100081 4 自然資源部活動構造與地質安全重點實驗室, 北京 100081

0 引言

新生代以來，印度板塊與歐亞板塊的碰撞和持續匯聚作用導致喜馬拉雅造山帶和青藏高原的快速崛起與不斷向外擴展(Molnar and Tapponnier, 1975).前人圍繞青藏高原深淺地質過程和變形模式提出了多種動力學模型，主要包括大陸側向逃逸模型(Tapponnier et al., 1982, 2001)、地殼連續變形模型(England and Houseman, 1986; Dewey et al., 1988)和下地殼流模型(Royden et al., 1997, 2008; Clark and Royden, 2000).側向逃逸模型認為青藏高原內部塊體沿深大斷裂發生側向擠出，其變形主要局限在塊體的邊界斷裂帶上；地殼連續變形模型以持續的擠壓動力、均勻的黏性物質和垂直連貫變形為特征，認為地表的隆升主要依靠地殼均勻的黏性增厚；下地殼流模型以上地殼高強度和無明顯縮短，以及上下地殼發生解耦為特征.各種變形機制和動力學假說都獲得了來自地震地質、地球物理與形變測量等學科的證據(Shen et al., 2001a; Chen et al., 2004; 王新勝等, 2013; Wang et al., 2016; Liu et al., 2014; Sun et al., 2012)，有關青藏高原隆升和動力學機制的爭論仍在繼續.實際上，整個青藏高原大陸變形很難用某一種變形模型去解釋，不同地區或塊體可能以某一種變形方式為主，同時兼具其他變形方式(Replumaz and Tapponnier, 2003; 許志琴等, 2010; Liu et al., 2014; 李海兵等, 2021)，高原形成機制與構造模式以及變形方式主次仍存在爭議的一個重要原因就是應力應變狀態沒有得到很好的限定(潘正洋等, 2020).

青藏高原東北緣是青藏高原擴張的前沿地帶.青藏高原快速隆升并向北東方向推擠，在其東北緣受到阿拉善塊體和鄂爾多斯塊體的阻擋，導致區內地殼發生縮短、抬升和左旋剪切，構造變形強烈且復雜(黃汲清等, 1980; 臧紹先等, 1992; 鄧起東等, 2002; 張培震等, 2003)，發育北西西向和北北西向兩組主導性活動構造，包括以左旋走滑為主的阿爾金斷裂帶和海原斷裂帶，以逆沖-左旋走滑為主的祁連山北緣斷裂帶、西秦嶺北緣斷裂帶和東昆侖斷裂帶，以左旋逆沖為主的六盤山斷裂，以及以右旋走滑為主的日月山斷裂帶和鄂拉山斷裂帶等邊界斷裂帶(圖1，Yin and Harrison, 2000; Tapponnier et al., 2001; 虢順民等, 2000; 張培震等, 2013; 鄧起東等, 2002).這些斷裂直接制約和影響了東北緣晚第四紀活動構造變形，形成了以擠壓推覆構造、次級剪切構造、剪切壓扁構造和弧形擠出構造為主的4種構造變形方式，且由活動斷裂控制的菱形塊體旋轉運動明顯(Duvall and Clark, 2010; 虢順民等, 2000; 袁道陽等, 2004)，形成了青藏高原東北緣強烈而獨特的構造環境，是我國西部擠壓構造向東部拉張構造轉換的過渡區域(Lasserre et al., 1999)，也是新生代構造變形最為強烈的地區之一.歷史上發生過1920年海原8.5級地震、1927年古浪8.0級地震、1932年祁連山北西端昌馬7.6級地震和2001年昆侖山口西8.1級地震.

地殼構造變形和地震的孕育均與地殼應力狀態密切相關，地震是在區域構造應力場作用下，巖體沿某一構造面發生破裂的結果(Scholz, 2002)，因而地殼構造應力場研究能夠幫助深入理解青藏高原東北緣的構造變形及孕震環境.地殼構造應力場的研究方法諸多，除水壓致裂法、應力解除法等可直接獲得應力狀態的原地應力測量方法外，通過地質學、地球物理學等相關方法反演也可提取應力信息(Zoback, 2007).地震震源機制不僅可以確定地震斷層類型，也反映了區域構造應力狀態，是研究區域構造應力場，特別是地殼深部應力狀態的有效方法之一(Gephart and Forsyth, 1984; Hardebeck and Michael, 2006; Wan et al., 2016； Li et al., 2018; Tian et al., 2019).地震各向異性與應力應變場密切相關，也是揭示巖石圈變形最直接的有效手段(Boness and Zoback, 2006; Gao et al., 2011; Crampin and Peacock, 2005; 高原等, 2018).利用近震剪切波分裂開展的上地殼各向異研究(張輝等, 2012; 郭桂紅等, 2015; Hu et al., 2020)認為青藏高原東北緣塊體內部地殼各向異性與研究區主壓應力有關，而深大斷裂帶附近的各向異性則與斷裂活動變形有關.部分學者基于XKS波分裂研究認為青藏高原東北緣殼幔變形具有強烈的耦合性，呈現垂直連貫變形的特征(常利軍等, 2016; Chang et al., 2017; Wang et al., 2016).也有學者利用接收函數確定地殼各向異性特征的分析方法，認為地殼各向異性主要來源于中下地殼，且受控于構造應力場及下地殼介質逃逸作用，同時揭示了存在局部殼幔解耦的可能(邵若潼等, 2019).利用震源機制解反演應力場的研究同樣認為東北緣部分區域殼幔以及上下地殼變形可能是解耦的(Han et al., 2019).前人通過地震各向異性或震源機制研究對青藏高原東北緣變形和構造應力場分布特征取得了一定認識，但關于研究區殼幔變形方式仍然存在爭議，并且缺乏聯合不同研究方法和研究成果的綜合性分析.

本文利用CAP全波形反演方法及HASH初動極性和振幅比方法，系統反演了32°N—41°N，96°E—107°E區域內ML≥3.0地震震源機制解781組，同時收集了Global CMT給出的MW≥4.6中強地震震源機制解96組，并利用STASI(阻尼區域應力場反演)算法對該區的構造應力場分布進行分析，同時聯合基于地應力實測獲得的地殼淺表層地應力狀態、基于地震各向異性獲得的地殼和上地幔變形特征以及基于GPS觀測獲得的地殼淺表層形變特征，討論了地殼淺表層、上地殼、中下地殼和上地幔的垂向變形特征，以期能夠對青藏高原東北緣的變形特征和動力學機制獲得較為深入的認識.

1 數據資料與計算方法

1.1 數據資料

本研究收集了研究區及周邊103個寬頻帶固定臺站2008年10月—2019年6月記錄的ML≥3.0以上地震波形數據(臺站分布見圖1b).為了保證研究結果的可靠性，僅從中選擇信噪比較高的地震波形記錄進行分析.同時，本研究補充采用了全球矩心矩張量目錄(GCMT)中1976年1月—2020年12月96組MW≥4.6地震事件的震源機制解數據，共同參與應力場反演.

圖1 青藏高原東北緣及區域構造背景圖(a) 青藏高原東北緣構造背景. 淺灰色圓形符號代表本文收集并進行反演的1976—2020年地震事件(ML≥3.0); 黑色曲線代表全新世活動斷裂, 灰色曲線代表更新世活動斷裂(鄧起東等, 2002); 紫色虛線表示圖12中地形剖面位置. (b) 青藏高原東北緣區域構造背景. 紫紅色三角符號代表本研究使用的寬頻帶固定臺站，淺紫色粗實線代表塊體邊界.Fig.1 Tectonic settings in and around the NE margin of the Tibetan Plateau(a) Geological structure in the NE margin of the Tibetan Plateau. Light gray circles denote local earthquakes with magnitudes equal to or greater than 3.0 occurred between 1976 and 2020 collected and inverted in this study; The black curves denote Holocene active faults, and the gray curves denote Pleistocene active faults (Deng et al., 2002); The purple dashed line indicates the position of the terrain profiles in Fig.12. (b) Regional tectonic settings in the NE margin of the Tibetan Plateau. The fuchsia triangular symbols represent 103 permanent seismic stations deployed in the NE margin of the Tibetan Plateau and the light purple thick solid lines represent block boundaries.

1.2 震源機制解計算

1.2.1 CAP全波形反演法及數據處理

本研究中ML≥3.5地震事件的震源機制解采用“裁剪-粘貼”(Cut-And-Paste，CAP)全波形反演方法(Zhao and Helmberger, 1994; Zhu and Helmberger, 1996).其主要思想是將寬頻帶地震記錄分成Pnl波(包含P波及其部分后續波形)和面波波段兩個部分進行反演.通過分別對Pnl和面波兩部分波形賦予不同的權重，計算其理論合成波形和真實觀測波形的互相關系數及誤差函數，并運用網格搜索方法搜索出理論合成波形與真實觀測波形差異最小時的地震震源機制解及最佳震源矩心深度.CAP方法具有對速度模型和臺站分布依賴性較小、計算結果相對穩定等特點.

采用CAP方法反演之前，需要對實測波形數據和理論波形數據做如下準備工作：首先對原始臺站波形數據進行預處理，主要進行去儀器響應并校正參考時間為發震時刻；而后將N-E-U坐標系下的三分量波形旋轉至R-T-Z坐標系下，并將單位轉換為cm·s-1；利用F-K頻率-波速法計算格林函數時(Zhu and Rivera, 2002)，采樣間隔為0.05 s，采樣點數為1024個.在震源機制反演過程中，對Pnl和面波設置不同的帶通濾波范圍，Pnl波為0.02～0.15 Hz，面波為0.02～0.12 Hz，這樣可以盡量減少背景噪聲對波形的干擾，亦能壓制地下小尺度結構散射造成的影響(陳偉文等，2012).為確保反演結果的可靠性，選取震中距500 km以內、位置分布均勻、三分量波形完整并且信噪比較高的臺站參與反演，去除理論波形與實際波形的擬合相關系數較低的波形.最終獲得研究區域287組3.5級以上地震事件的震源機制解.

1.2.2 HASH反演法及數據處理

對于研究區內一些地震震級較小(3.0≤ML<3.5)的地震事件，本文利用基于P波初動極性和S/P振幅比的HASH方法(Hardebeck and Shearer, 2002, 2003)開展了震源機制解反演.HASH方法的基本原理是基于不同方位、不同震中距的臺站的地震波形記錄，利用從震源向上射出的直達P波(Pg)和直達S波(Sg)引起的地動位移振幅比求解震源機制.其主要優點是可在一定程度上克服震源位置、速度模型和極性觀測誤差等，對于中小震地震有良好的適用性.在求解過程中為保證結果的可靠性，數據處理時僅選擇震中距200 km以內的、信噪比大于3的數字波形資料.利用SAC軟件人工讀取了每個臺站垂直向的初至波極性及S波到時，對于振幅的度量，將原始波形旋轉至R-T-Z坐標系下.P波振幅的度量：第一個半周期的峰值，量取Z和R方向的振幅，S波振幅的度量：三分量Sg波到時2 s內的最大值.在震源機制反演過程中，臺站分布范圍也是十分關鍵.在數據處理過程中，每個地震事件的相鄰臺站最大空隙角小于150°；離源角小于90°.HASH方法還從斷層不確定度、臺站分布比、振幅比的數量等參數對斷層面解的結果進行質量評價.在后文的討論中只采用了A、B、C三類的結果，共計494組.

1.2.3 速度模型的選取

現有結果表明，青藏高原存在著明顯的橫向不均勻性(李敏娟等, 2018; 許英才等, 2018).如果在研究區域采用單一的速度模型計算離源角，可能會造成較大的計算誤差.HASH方法能夠最多采用10個分區的速度模型，因此本文基于可靠的地殼速度結構測量，采用了多速度模型(青海省地震局日常定位用速度模型; 王周元, 1984; 張少泉等, 1985; 趙珠和張潤生, 1987; 周民都等, 2012; 金春華等, 2015)參與離源角計算并反演震源機制解(圖2).CAP方法在震源機制反演過程中，允許體波和面波進行時間平移，因此反演結果對速度模型的依賴性較小(韋生吉等，2009)，本文直接采用每個地震震中位置處Crust2.0一維速度模型計算格林函數(Bassin et al., 2000).

圖2 HASH反演法采用的地殼速度結構模型Fig.2 Velocity model used by HASH inversion method in this study

1.3 應力場反演方法

本研究采用Hardebeck和Michael(2006)提出的、基于震源機制解的MSATSI程序(Martínez-Garzón et al., 2014)進行構造應力場反演.該方法主要是基于線性反演方法，與網格搜索方法相比反演速度更快，能夠最大程度上擬合數據，消除人為影響.主要工作原理是首先將震源機制解按照分布區域劃分到若干個網格，然后采用阻尼最小二乘法進行反演，得到每個網格內的應力張量，并在最小二乘反演的過程中通過加入最優阻尼系數進行約束，以消除傳統應力反演方法中人為因素帶來的應力旋轉，最大限度減少相鄰區域間的應力差異.反演過程中，阻尼系數的確定是通過設定一系列阻尼值，得到應力場反演模型長度與數據擬合誤差關系的折中曲線.曲線拐點即為最佳的阻尼系數值，低于該值，模型復雜程度得到提高，但對反演殘差的改善作用不大，且觀測值與理論值匹配度急劇降低；高于該值，隨著模型的簡化，反演誤差會急劇增加(Hardebeck and Michael, 2006; 王曉山等, 2015; 郭祥云等, 2017).本研究將研究區劃分為1°×1°的網格進行構造應力場反演.根據雙力偶震源模型，存在兩個相同可能性的斷層面，即發震節面和輔助節面，因此MSATSI程序通過多次重采樣的方法，隨機選擇其中的一個節面進行反演，統計在一定置信度下的反演結果，從而確定最優解和不確定范圍.本研究反演過程中，在95%的置信區間內對原始數據進行2000次bootstrap重采樣.

1.4 R值反演方法

應力張量是對稱張量，有6個獨立參數，完整的應力張量需要三個主應力的量值及P、T和B軸方向，但是確定絕對應力大小是非常困難的(萬永革等, 2006).通常采用三個主應力的比例關系歸一化應力大小，引入應力形因子R表征相對應力大小(Gephart and Forsyth, 1984; 萬永革等, 2008)：

(1)

其中σ1、σ2、σ3分別表示最大、中間、最小主應力.當R值為0.5時，三個主應力成等差排列，σ1和σ3明確；當R值接近1時，σ2與σ3表現為一致的張應力狀態，并在與穩定的壓應力軸σ1垂直的平面內自由旋轉呈雙軸拉張狀態，兩軸不易區分；當R值接近0時，σ1與σ2表現為一致的壓應力狀態，并在與穩定的張應力軸σ3垂直的平面內自由旋轉呈雙軸壓縮狀態，兩軸不易區分(Guiraud et al., 1989; Wan et al., 2016; 劉兆才等, 2019).

2 結果與分析

2.1 反演計算結果

圖3給出了2012年6月3日18時49分青海共和MS4.2地震CAP反演結果，圖3a為在最佳擬合深度處各臺站的理論波形與實際波形擬合圖，其中約78%的分量擬合相關系數均高于80%，表明理論波形與實際波形具有較好的擬合關系.圖3b為該地震CAP波形反演誤差隨深度的變化圖，具有最小誤差的最佳擬合深度為6 km.

圖3 利用CAP方法反演2012年6月3日18時49分青海共和MS4.2級地震示例(a) 地震震源機制和理論波形與觀測波形擬合圖. 紅線代表理論地震圖, 黑線代表觀測地震圖, 波形下方數字分別代表理論地震圖相對觀測地震圖的偏移時間(單位: s)及二者的相關系數(%), 臺站下方數字分別代表震中距(單位: km)和方位角(單位: °); (b) 地震震源機制解反演殘差隨深度分布圖.Fig.3 Example of the focal mechanism determination by the cut-and-paste (CAP) method for an earthquake that occurred at 18∶49 on 3 June 2012(a) The corresponding waveform comparisons of synthetics (red) and records (black) at 8 stations. The two numbers under each segment are the time shift in second (upper) between the synthetic and record (positive means a delayed record) and the waveform correlation coefficient (lower, 100 indicates the best fit). Epicentral distances (in km) and azimuths (in degree) are given below the station codes; (b) Variation of fitting error with depth during the focal mechanism inversion of the earthquake occurred at 18∶49 on 3 June 2012.

為了驗證震源機制解反演結果的可靠性，選取研究區部分震級大于3.8級的中小地震，同時采用CAP方法和HASH方法反演其震源機制解，并對結果進行對比；選取5.0級以上中強地震，將CAP反演結果與GCMT結果進行對比，結果如表1和圖4a所示.不同方法的反演結果具有高度的一致性，說明采用的方法和反演結果可靠.

表1 不同方法求取震源機制結果比較Table 1 Comparison of different methods to obtain focal mechanism solutions

圖4 (a) 本研究獲得的781個地震事件的震源機制解分布, 其中綠色填充的膨脹區表示走滑型, 藍色填充的膨脹區表示逆斷型和逆走滑型(統稱為逆斷型), 紅色填充的膨脹區表示正斷型和正走滑型(統稱為正斷型), 灰色矩形內展示了分別利用不同方法得到的同一地震事件的震源機制解; (b) 從GCMT收集的96個地震事件的震源機制解分布, 沙灘球顏色標注同(a); (c) 本研究獲得的及從GCMT收集的共計877個地震事件的震源機制解P軸水平分量投影分布; (d) 877個地震事件的震源機制解T軸水平分量投影分布. 其他標注符號同圖1Fig.4 (a) 781 FMSs determined in this study. the focal mechanisms with green, blue, and red filled in expansion quadrants represent strike-slip, thrust, and normal fault earthquakes, respectively; the grey boxes show the comparison of the same events with GCMT, CAP and HASH solutions, respectively; (b) Same as (a) but for 96 FMSs collected from GCMT catalogue; (c) Horizontal projections of the P axes (blue line segments) of the 877 FMSs, including 781 events in this study and 96 events from the GCMT catalogue; (d) Same as (c) but for T axes (red line segments). The other labelling are the same as those in Fig.1

2.2 震源機制解總體特征

基于本文獲得的781個震源機制解，參照世界應力圖WSM(World Stress Map)的劃分原則(Zoback and Healy, 1992; Heidbach et al., 2018)，根據震源機制解P、T、B軸傾伏角的大小，將震源機制解類型劃分為6種，分別為正斷型(NF)、正走滑型(NS)、走滑型(SS)、逆走滑型(TS)、逆斷型(TF)和不確定型(U)，如表2所示.結果表明：研究區內走滑型地震為433個，占55.44%，正走滑型和正斷型地震為177個，占22.66%，逆走滑型和逆斷型地震為171個，占21.90%.研究區整體以走滑型地震為主，其次為正斷型和具有一定走滑分量的正斷型地震.

表2 震源機制解分類表Table 2 Categories of tectonic stress regime for focal mechanism

圖4a給出了本研究獲得的781個地震的震源機制解，圖4b為收集的GCMT中96個地震的震源機制解，圖4c和4d分別為本研究反演獲得的以及收集的共計877個震源機制解的主壓應力P軸和主張應力T軸的水平投影分布.整體上P軸呈NNE向～ENE向分布，受控于印度板塊北東向歐亞板塊持續俯沖和青藏高原物質東向擴散作用.P軸、T軸大多傾伏角較小(其中，P軸62.26%<30°，79.13%<45°；T軸52.57%<30°，74.91%<45°)，體現了近水平的擠壓應力作用，且分布特征與區內走滑型斷裂活動性質相對應，揭示了研究區所處的青藏高原東北緣塊體間強烈的剪切作用方式.

2.3 構造應力場特征

利用MSATSI程序對本研究獲得的781組震源機制解和收集的96組震源機制解數據進行應力場反演，最佳阻尼值取為1.1(圖5)，網格內震源機制解的最小個數為5個，得到各應力單元應力張量反演結果及最佳擬合構造應力張量(圖6、圖7、表3)，在此過程中，剔除了不確定性較大的數據.反演結果表明，研究區最大主壓應力軸σ1方向以NE向為主，且由巴顏喀拉塊體西南部向外呈扇形輻射，至祁連造山帶逆時針偏轉為NNE向，至秦嶺造山帶順時針偏轉為ENE向—ESE向，與鄂拉山斷裂、日月山斷裂、祁連山北緣斷裂、海源斷裂、東昆侖斷裂、西秦嶺北緣斷裂等邊界斷裂活動性質一致，且與青藏高原東北緣地形陡變帶跡線近垂直，表明構造變形和地殼隆升與構造應力場關系密切.最小主壓應力軸σ3方向以NW向—NNW向為主.σ1軸和σ3軸傾伏角較小，大部分近水平，表明研究區整體處于近水平的擠壓或剪切應力環境，局部最大主壓應力軸近垂直，處于拉張應力環境.中間主應力軸σ2傾伏角較大，多數近垂直.應力形因子R表征相對應力大小，為了更直觀的獲得研究區相對應力分布特征，采用插值法給出了連續的R值分布(圖7、圖9).揭示了研究區內R值具有明顯的空間差異分布特征，表明塊體間相互作用強烈.R值較小區域表明中間主應力值與主壓應力值較為接近，中間主應力軸表現為壓應力特征，R值較大區域表明中間主應力值與主張應力值較為接近，中間主應力軸表現為張應力特征.

圖5 模型長度與數據擬合誤差之間的折中曲線圖空心圓旁所標數字是阻尼參數, “+”符號表示所選擇的最佳阻尼系數.Fig.5 Trade-off between model length and misfit calculated from the corresponding damping parametersDamping values are shown beside each point and the best selected damping parameter is shown with a cross.

圖6 主應力軸投影分布圖黑色圓點表示主應力軸σ1、σ2、σ3; 深灰色扇形區表示主應力軸方向玫瑰圖; 等值線云圖表示主應力軸分布密度, 由冷色調藍到暖色調紅表示主應力軸投影分布密度由低到高.Fig.6 Stereonet map of principal stress axesThe black dots indicate the principal stress axes σ1, σ2, and σ3; The dark gray sector represents rose diagram of the principal stress direction; The contour map shows the distribution density of the principal stress axes, from cool blue to warm red, indicating the distribution density of the principal stress axes from low to high.

圖7 應力場反演結果紅色、綠色和藍色圓點分別表示95%置信度下最大主壓應力軸σ1、σ2和σ3的不確定范圍; 黑色實線加號表示最優解; 底圖顏色表示反演的應力形因子R值分布; 灰色曲線表示活動斷裂; 藍色曲線表示河流; 網格內數字代表該網格內參與反演的震源機制解個數.Fig.7 Stereomap of stress field inversion resultsThe coloured dots represent the stress axes in the 95 percent confidence interval by bootstrap, while the red, green and blue dots denote the orientations of the σ1, σ2 and σ3, respectively; the black crosses denote the best solutions for the principal stress axes in each grid; the color of base map represents distribution of the inverted stress shape factor R value; the gray curves represent active faults; the blue curves represent rivers; the number in the box at each grid shows the number of FMSs assigned to that grid in the stress inversion.

MSATSI程序采用Lund和Townend(2007)的方法計算了各應力單元的水平最大主應力SHmax方向，并基于WSM分類標準對其應力狀態進行了劃分(Zoback and Healy, 1992; Heidbach et al., 2018).筆者基于由MSATSI程序反演得到的最佳擬合構造應力張量，同樣采用Lund和Townend(2007)的計算方法，得到各應力單元水平最大主應力最優方向，結果如圖8所示.SHmax方向與最大主壓應力軸基本一致，同樣表現為由研究區西南部向外呈扇形擴展.研究區內應力結構類型以走滑型為主，祁連造山帶和柴達木盆地西部區域、隴西塊體東部區域和秦嶺造山帶則以逆斷型或無法明確應力結構類型為主，這些區域地震活動較為頻發，且地震震源機制類型多樣(圖4a)，地殼同時存在擠壓、剪切和拉張變形.此外，從《中國及鄰區現代構造應力場圖》中搜集研究區內基于水壓致裂和應力解除等原位實測方法獲得的地殼淺表層(<2 km)SHmax方向(謝富仁和崔效鋒, 2015)，并按照本研究應力單元劃分標準，將地殼淺表層SHmax方向進行分區平均，同樣繪制到圖8中.在研究區中部區域，地殼淺表層SHmax方向總體以NE向為優勢方位，同時存在近NS向和近EW向分布.西南部以NE向為主，東南部以NW向為主，自西向東呈現順時針偏轉，與震源深度SHmax方向分布特征一致.地殼淺表層實測SHmax方向與所在應力單元內震源層SHmax方向近平行或呈一定角度斜交，兩者夾角大多數在30°以內，表明地殼淺表層應力場與震源深度應力場具有較好的一致性，與深部運動即青藏高原物質東北向擴散相耦合.但同時地殼淺表層應力場影響因素較多，受地形地貌、構造、巖性等影響易形成局部應力場，因而與深部應力場存在一定的差異.在研究區東南部，即柴達木—隴西塊體、巴顏喀拉塊體與秦嶺造山帶、四川盆地北緣過渡區域，地殼淺表層實測SHmax方向由NE向偏轉為NW向，且與震源深度SHmax方向夾角較大，該區地殼深淺構造應力場以及變形特征可能存在較大的差異性.

3 討論

3.1 區域構造應力場與動力學過程

據中國構造應力分區圖(謝富仁和崔效鋒, 2015)，研究區及周邊區域被劃分為祁連山—河西走廊應力區、海原—六盤山應力區、柴達木盆地應力區、西秦嶺應力區和巴顏喀拉應力區，其中，除祁連山—河西走廊應力區為擠壓區外，其余應力區均為剪切區(圖9).該應力分區主要考慮塊體和斷裂相互作用的影響，將應力狀態具有較好一致性的區域劃分為同一分區.本研究表明，構造應力在空間上的非均勻性分布特征較為顯著，在現有應力分區基礎上，可進行次一級或局部應力分區.根據應力形因子R值的定義(式(1))，圖9中暖色調區域表示雙軸拉張區，冷色調區域表示雙軸壓縮區.雙軸壓縮區以祁連造山帶，隴西塊體東部，以及由鄂拉山斷裂帶、西秦嶺北緣斷裂帶、貴德斷裂和東昆侖斷裂帶圍限區域為主；雙軸拉張區以柴達木盆地、秦嶺造山帶、隴西塊體西部和巴顏喀拉塊體為主.雙軸拉張區內主應力軸方向以σ1明確為特征，表現為NE-SW向壓縮和NW-SE向拉張，且擠壓應力大于拉張應力.雙軸壓縮區內主應力軸方向以σ3明確為特征，表現為NE-SW向擠壓和NW-SE向拉張，拉張應力可能遠小于擠壓應力，也可能與擠壓應力相差不大.Wang和Shen(2020)基于1991至2016年GPS觀測結果，得到了中國大陸地殼形變場.其研究結果表明，在祁連山地區以NE-SW向縮短和NW-SE擴展為特征，并且縮短速率明顯高于擴展速率，而在東昆侖斷裂帶顯示出不低于20×10-9/yr的左旋應變率.本文研究揭示的應力機制與基于GPS觀測揭示的地殼形變特征基本吻合，表明研究區應力場與應變率場基本耦合，地殼形變主要受控于背景構造應力場(圖9).綜合區內R值、主應力分布特征，表明研究區總體主要受控于印歐板塊碰撞造成的NE向擠壓，并在這種擠壓背景下，存在高原物質向北東推擠受周邊穩定剛性塊體阻擋后由重力位能垂向差異產生的拉張力，并向外擴展，最終形成研究區內現今走滑為主的應力結構類型(圖8).同時區內構造受控于區域應力場，吸收拉張和剪切構造變形，在地質演化時期內存在由擠壓向走滑構造轉變的過程，如鄂拉山斷裂在第四紀初期由擠壓逆沖轉換為右行走滑(袁道陽，2004).

表3 1°×1°網格劃分反演得到的應力場參數Table 3 The inverted stress field parameters on grid of 1°×1°

續表3

圖8 (a) 水平最大主應力SHmax分布. 中心為圓形符號短線表示震源機制解反演結果，其中綠色、藍色、紅色和黑色短線分別代表走滑型、逆斷型、正斷型和過渡型應力狀態; 中心為菱形符號短線表示由水壓致裂和應力解除等原位實測方法得到的最大水平主應力方向(謝富仁和崔效鋒, 2015), 應力單元左下角或左上角數字角度代表該單元內震源深度SHmax方向與地殼淺表層SHmax方向的夾角; 32°N—33°N緯度行及96°E—97°E經度列左上角圓圈內數字代表應力單元坐標, 緯度行對應y坐標, 經度列對應x坐標. (b) 震源深度SHmax方向與地殼淺表層SHmax方向夾角統計圖Fig.8 (a) Distribution of the maximum horizontal stress orientations. Short lines with a circular symbol in the center represent SHmax orientations obtained by inversion of focal mechanism solutions, and the green, blue, red and black symbols represent strike-slip, reverse, normal and transitional stress states, respectively; Short lines with a diamond symbol in the center represent SHmax orientations obtained by in-situ measurement methods such as hydraulic fracturing and stress relief (Xie and Cui, 2015); The numbers in the lower left or upper left corner of the stress unit represent the angles between SHmax orientation of the focal depths and those of the shallow crust; The number in the circle in the upper left corner of 32°N—33°N latitude row and 96°E—97°E longitude column stand for the coordinate of the stress unit, the latitude row corresponds to the y coordinate, and the longitude column corresponds to the x coordinate. (b) Statistical diagram of the angle between SHmax orientation of the focal depths and those of the shallow crust

圖9 青藏高原東北緣動力模型圖彩色底圖為由震源機制解反演得到的R值分布, 其中暖色調區域表示雙軸拉張區, 冷色調區域表示雙軸壓縮區; 應力分區引自中國構造應力分區圖(謝富仁和崔效鋒, 2015); 紫紅色互相垂直箭頭代表基于GPS得到的主應變率及其方向(Wang and Shen, 2020); 淺灰色寬箭頭代表塊體運動方向.Fig.9 The dynamic model of the NE margin of the Tibetan PlateauColor base map is the distribution of R value obtained by inversion of focal mechanism solutions, the warm color area represents the biaxial stretching area, and the cool color area represents the biaxial compression area. Stress zones is quoted from China's tectonic stress field map (Xie and Cui, 2015). The fuchsia mutually perpendicular arrows represent the principal strain rate and its direction based on GPS observations (Wang and Shen, 2020); The light gray wide arrows represent the direction of movement of the blocks.

圖10 地殼不同深度應力場分布特征(a,c,e) 應力場反演結果(95%置信區間), 其他符號代表意義同圖7; (b,d,f) 水平最大主應力分布(95%置信區間, 應力狀態分類采用WSM標準: 綠色代表走滑型, 藍色代表逆斷型, 紅色代表正斷型, 黑色代表過渡型), FD: 震源深度.Fig.10 Stress field of different depths of the crust(a,c,e) Stereomap for the principal stress field (within 95 per cent confidence region) of the shallow and deep crust (see Fig.7 for the other labelling); (b,d,f) Maximum horizontal stress orientations (SHmax, within 95 per cent confidence region with stress regime categories according to the WSM standard: Blue, green, red and black lines denote thrust, strike-slip, normal and mixed environment, respectively) in different depths of the crust; FD: Focal depth.

圖11 青藏高原東北緣構造應力場及各向異性分布特征(a) 不同深度地殼構造應力場分布, 虛線矩形區指示(b)范圍; (b) 青藏高原東北緣東南部區域深淺地殼構造應力場分布(<2 km, ≤10 km, 和>20 km), 其中2km以淺SHmax方向由原位地應力實測獲得, 震源深度SHmax方向由震源機制解反演獲得, 背景圖為該區20 km深度地殼橫波速度(Shen et al., 2016); (c) 近震直達S波分裂(Hu et al., 2020; 錢旗偉等, 2017; 張輝等, 2012); (d) Ps波分裂(Xu et al., 2018; Wang et al., 2016; 邵若潼等, 2019); (e) XKS波分裂(Chang et al., 2017); (f) 青藏高原東北緣深部動力學機制示意圖, 修自Clark and Royden (2000).Fig.11 Stress field and seismic anisotropy of the NE margin of the Tibetan Plateau(a) Characteristics of stress fields in different crustal depths, the dotted rectangular area indicates the range of (b); (b) Stress fields of the shallow and deep crust (<2 km, ≤10 km, and >20 km) in southeastern of the NE margin of the Tibetan Plateau, SHmax within 2 km are obtained by in situ stress measurement methods, SHmax in seismogenic depths are obtained by focal mechanism inversion, the background colours denote the average S-wave velocity at depth of 20km (Shen et al., 2016). (c) Local S wave splitting (Hu et al., 2020; Qian et al., 2017; Zhang et al., 2012); (d) Ps wave splitting (Xu et al., 2018; Wang et al., 2016; Shao et al., 2019); (e) XKS wave splitting (Chang et al., 2017); (f) Deep dynamic mechanism model of the NE margin of the Tibetan Plateau, modified from Clark and Royden (2000).

利用本研究CAP反演法得到的最佳震源矩心深度，給出了不同地殼深度(≤10 km, 10～20 km和>20 km)的應力場分布特征(圖10, 圖11a).由于按地殼深度分層后地震活動空間分布差異較大，在該部分研究中每個網格采用的震源機制至少為1個，盡管部分網格數據量相對較少，但程序采用阻尼最小二乘法反演同時得到每個網格點的應力張量并進行平滑，依然可以得到較為可靠的結果(Menke, 1989; 劉兆才等, 2019).結果表明不同深度的SHmax方位均以NE向為主，總體上地殼構造應力場分布具有垂向一致性，且與GPS速度場總體特征一致，受控于印歐板塊持續碰撞導致的青藏高原快速隆升并向北東方向推擠作用.Han等(2019)通過對比深(>20 km)淺(<20 km)地殼應力場，得到與本文較為一致的應力場分布特征.圖11a明顯揭示出SHmax方向由祁連造山帶和柴達木—隴西塊體向秦嶺造山帶呈現由NE→ENE→ESE向偏轉的特征，大致與地形差異邊界或地塊碰撞邊界垂直.已有研究表明，地震震源深度的垂向分布受到巖石圈熱流變性質的控制，地震主要發生于巖石圈的脆性域內(Ranalli and Murphy, 1987; 張國民等, 2002)，因此震源深度可以反映地殼巖石圈脆性變形層的厚度.從圖11a可以看出，地震多集中于20 km以淺，中下地殼多以塑性為主，青藏高原東北緣構造應力場分布特征似乎是深部軟弱物質北東向擠出并向周緣擴散作用在淺部的響應.在柴達木—隴西塊體和巴顏喀拉塊體東南部與秦嶺造山帶和四川盆地的過渡區域，中下地殼應力場與上地殼應力場存在差異，深部地殼變形可能與上地殼變形解耦.考慮到震源深度的不確定性，為了更為清楚的表示深淺地殼應力場分布特征，圖11b給出了10 km以淺和20 km以深由震源機制解反演得到的SHmax方向，同時疊加了2 km以淺由原位地應力測量方法實測得到的SHmax方向.可以看出，隴西塊體10 km以淺地殼深度SHmax方向耦合性較好，而秦嶺造山帶和巴顏喀拉塊體東南部不同地殼深度SHmax方向差異較大(圖11b中虛線橢圓范圍內).2 km以淺SHmax方向以WNW-ESE為主，10 km以淺SHmax方向以WSW-ENE為主，而20 km以深SHmax方向同時具有WNW-ESE和WSW-ENE兩種優勢方位.Han等(2019)認為深部地殼SHmax方向一般垂直于低到高的橫波速度邊界.圖11b顯示，在橫波速度高低過渡區域地殼SHmax方向出現明顯的解耦現象，深部SHmax方向耦合于下地殼流方向(Clark and Royden, 2000; Zhao et al., 2021)，與青藏高原東南緣共同構成下地殼弱物質東向流出的兩個分支通道(Clark and Royden, 2000; Sun et al., 2012).

3.2 殼幔變形特征與變形機制

青藏高原東北緣近震各向異性研究表明上地殼各向異性同時受到區域應力與斷裂構造的雙重約束(圖11c).快波偏振以NE向為主，與SHmax方向一致，同時耦合于地表形變場(張輝等, 2012; 郭桂紅等, 2015; 錢旗偉等, 2017; Wang and Shen, 2020).該區與阿拉善塊體交接區受到NE向強烈擠壓作用，特別是祁連造山帶受到南北兩側阿拉善塊體和柴達木塊體的俯沖以雙向逆沖形式擠出，匯聚速率達16 mm·a-1(許志琴等, 1999)，表明該區經歷著NE向地殼縮短.阿拉善塊體、祁連塊體、柴達木—隴西塊體交匯區內，同時表現出與斷裂構造一致的NW向優勢偏振方向，且與上地幔XKS偏振方向和板塊絕對運動方向(APM)一致(圖11c, e)，殼幔變形方式同時存在差異性和關聯性，與塊體交匯區復雜的深部動力學機制相關.柴達木—隴西塊體東南部和秦嶺造山帶快波偏振方向則以WNW向—NW向為主，與上地幔XKS偏振方向以及中下地殼SHmax方向近似平行(張輝等, 2012; 郭桂紅等, 2015; Chang et al., 2017)，與上地殼SHmax方向呈大角度斜交，因此主要受斷裂構造影響.

基于遠震接收函數數據的研究揭示了青藏高原東北緣中下地殼顯著的各向異性特征，其快剪切波方向整體以WNW-ESE或NW-SE向為主(圖11d)，大致與斷裂構造及造山帶走向一致且垂直于水平最大主壓應力方向，從而將青藏高原東北緣中下地殼各向異性歸因于擠壓碰撞相關的純剪切引起的角閃石晶格優勢排列(Wang et al., 2016; Xu et al., 2018; 邵若潼等, 2019).但在共和盆地和青海湖盆地交匯區域，以及鄂爾多斯塊體西緣則呈現以NE向—ENE向為主的快波偏振方向(圖11d橢圓形區域)，與SHmax方向、GPS速度場、上地殼各向異性等均具有較好的一致性，而與近NW向分布的XKS偏振方向呈大角度斜交或近于垂直(圖11f)，表明殼幔具有不同的變形機制.在秦嶺造山帶及其周邊區域，Ps快波方向與中下地殼SHmax方向一致，而與上地殼SHmax方向解耦，地殼通道流的存在可能是造成該區深淺地殼應力場解耦以及中下地殼各向異性的主要原因(Han et al., 2019; Zhao et al., 2021).

相對于青藏高原東北緣大多小于0.4 s的地殼各向異性時間延遲(張輝等, 2012; 郭桂紅等, 2015; 錢旗偉等, 2017; Wang and Shen, 2020; Wang et al., 2016; Xu et al., 2018; 邵若潼等, 2019)，XKS波分裂結果估算的時間延遲平均為1.3 s，因此主要反映了上地幔各向異性的特征(常利軍等, 2021).Silver(1996)認為大陸上地幔各向異性受構造應力作用引起的巖石圈變形影響顯著，構造活動區快波方向通常與大型走滑斷裂和造山帶構造方向一致.XKS波分裂研究表明青藏高原東北緣上地?？觳ǚ较蛞訬W-SE為主(圖11e)，與區域走滑斷裂和造山帶構造走向一致，揭示了強烈的NE-SW擠壓應力環境，與由震源機制解揭示的NE-SW向主壓應力一致.而在秦嶺造山帶及周邊區域，XKS偏振方向偏轉為E-W或WNW-ESE向，且與中下地殼SHmax方向一致，前人根據該區較大的XKS慢波時間延遲和較薄的巖石圈，認為秦嶺造山帶是青藏高原巖石圈物質的擠出通道，同時也是軟流圈物質東流的地幔流通道(王瓊等, 2013; Chang et al., 2017; 常利軍等, 2021).根據地殼應力場和上地幔各向異性特征，不能排除該區存在殼幔解耦變形的可能.

3.3 青藏高原東北緣動力學特征

綜合以上研究，本文初步提出了青藏高原東北緣巖石圈多種變形機制共存的模式(圖11f, 圖12).高原物質向北東推移過程中受到阿拉善和鄂爾多斯兩個穩定剛性塊體的阻擋作用，總體上沿NW-SE向發生了伸展變形，形成了一系列NW-SE方向的造山帶和斷裂帶，在上地幔產生了NW-SE方向的快剪切波(王瓊等, 2013; Chang et al., 2017; 常利軍等, 2021).在阿拉善塊體與柴達木盆地之間的祁連造山帶北部，地殼各向異性也以NW-SE向分布的快波偏振方向為主要特征(張輝等, 2012; 錢旗偉等, 2017; Wang et al., 2016; Xu et al., 2018)，SHmax方向則與快波偏振方向近垂直，體現了以巖石圈連續增厚變形模式為主.在柴達木盆地東南部地區以及鄂爾多斯塊體西緣，地殼剪切波快波方向以NE-SW向為主，各向異性來源于由SHmax作用導致的沿應力方向定向排列的微裂隙(張輝等, 2012; 錢旗偉等, 2017; Wang et al., 2016; Xu et al., 2018)，地殼和地幔的變形機制可能不同.柴達木盆地強度低于阿拉善塊體和鄂爾多斯塊體，對下地殼弱物質流沒有形成有效的阻擋，物質流可能繼續向北東方向推移直至受到阿拉善塊體阻擋，繼而產生伸展變形以及地殼增厚.鄂爾多斯塊體西緣是構造活動強烈的青藏高原向穩定的鄂爾多斯塊體的過渡地帶，深部運動和構造復雜，殼?？赡艽嬖诮怦瞵F象(許英才等, 2019).柴達木—隴西塊體和巴顏喀拉塊體東南部與秦嶺造山帶和四川盆地過渡區域殼幔剪切波快波偏振方向均以近EW向或WNW向為主(張輝等, 2012; 錢旗偉等, 2017; Hu et al., 2020; Wang et al., 2016; Xu et al., 2018; 邵若潼等, 2019; Chang et al., 2017; Han et al., 2019)，中下地殼SHmax方向與殼幔剪切波快波偏振方向一致(圖8, 圖11a,b)，耦合于下地殼流方向，與青藏高原東南緣共同構成下地殼弱物質東向流出的兩個分支通道(Clark and Royden, 2000; Sun et al., 2012).上地殼SHmax方向與快波偏振方向斜交，深淺地殼SHmax方向存在一定差異，深部地殼變形與淺部發生解耦.

總體而言，青藏高原東北緣殼幔變形機制和深部動力學特征極為復雜，很難用單一模式解釋，可能同時存在多種變形方式.根據地殼構造應力場分布特征并結合殼幔各向異性，青藏高原東北緣可能以巖石圈連續增厚變形與下地殼通道流的共同作用模式為主，共同影響并導致地殼增厚和高原隆升(圖12).青藏高原北東向持續推擠，向北受到具有古老克拉通性質的阿拉善塊體阻擋，發生地殼均勻增厚導致的地殼縮短，從而形成陡峭的地形邊緣，如圖12所示A-A′剖面在約50 km范圍內地形陡變近3000 m.同時，受青藏高原北東向擠出作用以及上地幔物質加熱作用，巴顏喀拉塊體深部地殼物質發生部分熔融并產生北東向擴展的下地殼通道流(Liu et al., 2014; Zhao et al., 2021)，導致下地殼增厚并伴隨著地殼流的向外擴展形成較為平緩的地形，如圖12所示B-B′和C-C′剖面.而D-D′剖面位于兩種變形方式的過渡區域，其地形起伏特征介于兩種典型地形之間.

圖12 (a) 青藏高原東北緣深部動力學機制三維示意圖, 藍色箭頭代表高原物質逃逸方向, 紫色箭頭代表下地殼變形方向, 白色實線代表(b)中地形剖面; (b) 青藏高原東北緣典型地形剖面, 剖面位置如(a)及圖1所示Fig.12 (a) 3D perspectives of deep dynamic mechanism model of the NE margin of the Tibetan Plateau, the blue arrows indicate orientation of material escape, the purple arrows indicate the lower crustal deformation, the white solid lines represent the terrain profile in (b); (b) Typical topographic profiles for the NE margin of the Tibetan Plateau, the locations of profiles are shown in (a) and Fig.1

4 結論

本研究利用CAP全波形反演方法及HASH初動極性和振幅比方法，系統反演了研究區內ML≥3.0地震震源機制解781個，同時收集了MW≥4.6中強地震震源機制解96個，并使用MSTASI程序進行了應力場反演，獲得了研究區應力張量、R值及水平最大主應力優勢方向，詳細分析了該地區震源機制解和構造應力場特征，最后根據計算和分析結果，聯合地表GPS觀測和殼幔各向異性研究，綜合討論了青藏高原東北緣變形機制和動力學特征.分析結果顯示：

(1)研究區震源機制解類型以走滑型為主，P軸、T軸大多傾伏角較小，最大主壓應力軸σ1方向以NE向為主，且由巴顏喀拉塊體西南部向外呈扇形輻射，至祁連造山帶逆時針偏轉為NNE向，至秦嶺造山帶順時針偏轉為ENE向—ESE向，分布特征與區內走滑型邊界斷裂活動性質相對應.σ1軸和σ3軸傾伏角較小，大多近水平，表明研究區整體處于近水平的擠壓或剪切應力環境.SHmax方向與最大主壓應力軸分布特征基本一致，應力結構類型以走滑型為主.

(2)青藏高原東北緣最大主壓應力方向與最大主應變方向具有很好的一致性，應力機制與基于GPS觀測揭示的地殼形變特征基本吻合，表明研究區應力場與應變場基本耦合，地表形變主要受控于區域構造應力場.大部分地區不同地殼深度的SHmax方向一致性較好，而柴達木—隴西塊體和巴顏喀拉塊體東南部與秦嶺造山帶和四川盆地的過渡區域，深淺地殼SHmax方向具有差異性，深部地殼變形與淺部可能發生解耦.地殼構造應力場可能是深部軟弱物質北東向擠出并向周緣擴散作用在淺部的響應.

(3)青藏高原東北緣殼幔變形機制和深部動力學特征極為復雜，很難用單一模式解釋.阿拉善塊體與柴達木盆地之間的祁連造山帶北部地區，殼幔變形垂向連貫，體現了以巖石圈連續增厚變形模式為主.柴達木盆地東南部及鄂爾多斯塊體西緣，地殼各向異性主要來源于由SHmax作用導致的沿應力方向定向排列的微裂隙，殼幔變形機制可能不同.青藏高原東北緣東南部區域，深部地殼變形可能與淺部發生解耦，至鄂爾多斯塊體與四川盆地之間的秦嶺造山帶區域，中下地殼SHmax方向與殼幔剪切波快波偏振方向一致，并耦合于下地殼流方向.青藏高原東北緣可能以巖石圈連續增厚變形與下地殼通道流的共同作用模式為主.