張 晨，季靈運，朱良玉，徐 晶

(中國地震局第二監(jiān)測中心，陜西 西安 710054)

0 引言

印度板塊與歐亞板塊之間的相互擠壓碰撞形成了青藏高原，這是地質歷史上最重要的造山事件。喜馬拉雅東構造結地區(qū)屬于該構造碰撞的前緣，位于青藏高原東南緣、印度大陸東北以及緬甸板塊西北角。該地區(qū)地貌復雜，構造活動劇烈，且地震頻發(fā)。1900年以來，該地區(qū)發(fā)生5.5級以上地震近百次，且大多數(shù)集中于區(qū)域內部的主要地震帶上。因此，該地區(qū)被稱為地震研究天然實驗場(楊建亞等，2017；王凱悅等，2021)。

地殼應力是地震孕育和發(fā)生的直接驅動力(Hardebeck，Michael，2006)，地殼應力場的方向也是評估斷層潛在滑移的一個重要參數(shù)，地殼應力的時空變化能夠較好地反映地震與應力之間的關系特征，因而能在地震危險性和風險評估中發(fā)揮至關重要的作用(Heidbach，Ben-Avraham，2007；Martinez-Garzon，2014；Walsh，Zoback，2016；吳微微，2020)?，F(xiàn)有推算構造應力場的方法主要包括：測井擴張破裂、實地應力測量中的水壓破裂和應力解除(Zoback，1992；Reynolds，2006)、震源機制解(徐紀人等，2008；王曉山等，2015；黃驥超等，2016；崔華偉等，2017)、斷層滑移測量(Angelier，1979)等。其中，利用震源機制解求解構造應力場是研究應力構造的主要方法(謝富仁等，1993；崔效鋒，謝富仁，1999；Hardebeck，Michael，2006；Wan，2010；Martinez-Garzon，2014；楊帆等，2019；王曉山等，2020)。在所有利用震源機制解反演應力場的方法中，阻尼應力張量反演(SATSI)是較為成熟的計算方法(Hardebeck，Michael，2006；Zhao，2013；Luo，2015；王曉山等，2015)，該方法被廣泛應用于青藏高原地區(qū)的應力場研究(徐紀人等，2008；侯強等，2018；李君等，2019)。

隨著現(xiàn)代大地測量學的發(fā)展，GPS大地測量技術已成為地震監(jiān)測以及形變測量的主要手段(王敏，沈正康，2020)，國內外學者利用GPS數(shù)據(jù)獲取了許多構造活躍區(qū)域的應變率場(Allmendinger，2007)，尤其是喜馬拉雅造山帶(Bilham，1997；Vernant，2015)、青藏高原(張培震等，2002；張清志等，2005)以及阿爾金斷裂帶(鄭文俊等，2009)。前人利用GPS觀測手段和震源機制解，獲取全球范圍的應力方向分布圖(Zoback，1992)，該結果為研究全球地殼應力應變特征提供了指導作用，但針對喜馬拉雅東構造結及其周邊地區(qū)構造應力場的精細分析尚顯不足。

因此，本文利用最新發(fā)表的GPS速度場數(shù)據(jù)，計算喜馬拉雅東構造結及其鄰區(qū)的地殼應變率場，然后收集該地區(qū)≥2.0地震的震源機制解，利用MSATAI程序包反演研究區(qū)的應力場，最后對比GPS主應變方向與震源應力場方向，分析研究區(qū)地殼深淺部構造變形特征。

1 GPS數(shù)據(jù)與應變率場

為了獲取研究區(qū)的現(xiàn)今GPS應變率場，本文使用Wang和Shen(2020)收集的1991—2016年的長期GPS觀測數(shù)據(jù)(圖1)，參考框架為歐亞大陸，采用Shen等(2015)提出的一種通過速度場計算應變的方法。該方法較為穩(wěn)健，不依賴于對數(shù)據(jù)的某種假設，且平衡了穩(wěn)定性與分辨率兩方面，計算結果更為全面可靠。計算結果如圖2所示。

圖1 喜馬拉雅東構造結及其鄰區(qū)速度場

由圖2a可知，最大主應變分布自西向東、由北至南呈現(xiàn)順時針偏轉。區(qū)域內部的大型走滑斷裂，其最大主壓應變也較高，表明其具有較為強烈的構造變形活動。在96°E以西、30°N以北的青藏高原地區(qū)，應變以拉張為主，方向由近EW向逐漸轉為NW向；而30°N以南的喜馬拉雅區(qū)域，主要以SN-NE向擠壓為主，主壓應變方向近似與板塊邊界相垂直，表明該地區(qū)受到來自印度板塊的NNE擠壓作用。在96°E以東地區(qū)，主要包括地震活動最劇烈的青藏高原東南緣區(qū)域——川滇菱形塊體，該區(qū)域受華南地塊的阻擋，主應變復雜多變，整體表現(xiàn)為繞喜馬拉雅構造結呈SN-NNW-NW向的順時針偏轉。巴顏喀拉塊體內部主壓應變方向為近EW向，川滇菱形塊體東邊界主壓應變方向從北向南逐漸由近EW向轉為近NW向。最大剪應變率高值區(qū)主要位于川滇菱形塊體東邊界，甘孜—玉樹斷裂帶、鮮水河斷裂帶、安寧河斷裂帶、則木河斷裂帶、小江斷裂帶及其周邊次級斷裂(圖2b)。本文研究區(qū)應變率場結果與前人研究得出的空間分布特性一致。

圖2 喜馬拉雅東構造結及其鄰區(qū)面應變率(a)和最大剪應變率(b)分布

2 震源機制解數(shù)據(jù)與應力反演

本文使用的數(shù)據(jù)為1976年1月1日至2020年9月25日喜馬拉雅東構造結及周邊地區(qū)(23°～34°N，90°～105°E)的7 092個震源機制解，對應震級范圍為2～7.9級。其中，420個震源機制解數(shù)據(jù)來源于GCMT(Global Centroid-Moment-Tenso)，對應震級范圍為4.7～7.9級，其余數(shù)據(jù)來源于四川地震臺，震級范圍為2～6.2級。按照Aki和Richards(1980)的分類方法對震源機制解數(shù)據(jù)進行分類(圖3)，分別是：正斷層，逆斷層，走滑斷層以及斜滑斷層，具體分類依據(jù)見表1。

圖3 喜馬拉雅東構造結及其鄰區(qū)震源機制解分布

表1 本文研究震源機制解劃分依據(jù)

通常情況下，由于研究區(qū)空間應力場分布的非均勻性，我們常將其劃分成更為細小的多個子區(qū)域，并獨立擬合出每個子區(qū)域的應力張量，但這種模型下的應力場空間變化存在難以解釋的情況。因此本文采用阻尼應力張量反演法盡可能消除應力變化產生的虛假信息，該方法能使相鄰子區(qū)域的應力張量分量之間的差異最小化，也令震源機制解與應力張量的匹配程度達到最高。鑒于該方法在空間應力場分布反演上具有極大潛力，Material-Garzon等(2014)在SATSI算法的基礎上開發(fā)了MATLAB環(huán)境中運行計算的MSATAI程序包。在進行應力張量反演時，將震源機制解數(shù)據(jù)按照0.5°×0.5°的網格進行分區(qū)，共計660個子區(qū)域，如圖4a所示。在輸入參數(shù)設置時，為了保證更多的子區(qū)域能夠反演得到應力結果，將每個網格地震最少個數(shù)設置為1，但由于每個子區(qū)域至少需要包含4個地震數(shù)據(jù)才能計算，此時，程序會根據(jù)該子區(qū)域周圍的震源機制解對該區(qū)域的應力張量反演進行約束。最后，對于該區(qū)域，本文選取最佳阻尼系數(shù)為1.2(圖4b)。

圖4 0.5°×0.5°網格地震數(shù)目統(tǒng)計(a)及模型擬合曲線圖(b)

本文的反演結果提供了3個主應力軸的方向以及反映3個主應力相對大小的參數(shù)：應力形因子(Pesicek，2010)，可表示為：

(1)

式中：、和分別是最大主壓應力、中間主壓應力以及最小主壓應力；主要反映中間應力更接近于最大主壓應力或是最小主壓應力，其取值為0～1。當<05時，求得的最大主張應力軸(最小主壓應力軸)方向可靠性較最大主壓應力軸的方向相對較高，此時，中間應力呈現(xiàn)為壓應力，與的大小表現(xiàn)為負相關，即越小，中間應力軸壓應力越明顯；當>05時，最大主壓應力軸的方向可靠性更大，中間應力呈現(xiàn)張應力，與的大小表現(xiàn)為正相關。極端情況下，=1時，最大主壓軸確定、最大主張軸不確定；=05時，兩軸均確定；=0時，最大主張軸確定、最大主壓軸不確定。

利用上述震源機制解數(shù)據(jù)，通過MSATAI程序包反演得到研究區(qū)的應力結果，如圖5a所示，藍色線段代表最大主壓應力軸，紅色線段代表最小主壓應力軸。線段的長度表示應力軸傾角的大小，線段越短傾角越大，越長則傾角越小。

圖5 震源應力場反演結果(a)、應力形因子R與最大水平應力SHmax分布結果(b)

為了更加清晰直觀地展示該區(qū)域的應力狀態(tài)，本文基于Lund和 Townend(2007)提出的最大水平應力計算方法，獲取了該區(qū)域的最大水平應力分布情況，如圖5b所示，圖中紅色直線表示最大水平應力max，最大水平應力的背景值為應力形因子。

3 研究結果

3.1 應力場反演結果

從圖5a可以看出，在99°E以西地區(qū)，為NE向，呈NWW或NW走向；在羌塘塊體內部，垂直于斷裂走向，則與斷裂走向平行；自西向東、由北向南存在NNE-NE向偏轉，軸傾角也逐漸增大，其中怒江、瀾滄江斷裂處軸傾角近垂直(Qiu，Qiao，2017)；自西向東呈現(xiàn)出NWW-NW向偏轉，且該區(qū)域的軸傾角較小，接近于水平。在緬甸弧俯沖帶區(qū)域指向NNE向，指向SEE向。

在99°E以東地區(qū)，由北至南呈現(xiàn)EW-SEE-NS的偏轉，軸傾角也存在較大的變化。巴顏喀拉塊體內部、川滇菱形塊體東邊界軸傾角較小，其西邊界軸傾角較大；表現(xiàn)為NS-NNE-EW的順時針偏轉，軸傾角較小。該區(qū)域主要是地震活動最強烈的川滇菱形塊體，應力方向在區(qū)域內存在很大差異，川滇菱形塊體東邊界主要為主壓應力方向是SEE向的走滑斷層，其西邊界則主要是主壓方向為SSE向的正斷結構，其內部西北區(qū)域為正斷型應力狀態(tài)，主要集中于金沙江斷裂帶以及理塘斷裂，龍門山斷裂則主要為逆斷型應力狀態(tài)。

從圖5b可以看出，區(qū)域內最大水平應力max方向呈現(xiàn)出較為明顯的分區(qū)特征。96°E以西地區(qū)，max主要為NNE向，且在羌塘地塊內部<0.5，表明該區(qū)域處于雙軸壓縮狀態(tài)；96°～99°E范圍內，max方向呈現(xiàn)出NNE-NE-WE的偏轉，該范圍屬于應力軸方向偏轉的過渡區(qū)域，值大部分都大于0.5，表明與雙軸拉張；99°E以東地區(qū)，max表現(xiàn)為圍繞喜馬拉雅構造結的順時針偏轉特征。本文所獲應力場分布與已有研究結果基本一致(徐志剛，2017；Zhao，2013)。

本文所反演的應力場結果與該區(qū)域內的地震類型較為吻合，可見，使用該方法能較為準確地反演出區(qū)域震源應力場狀態(tài)，這為我們分析地殼深部構造應力方向提供了一種可靠約束。

3.2 GPS主應變率與最大水平應力SHmax

GPS觀測可提供地殼淺部的應變結果，而震源機制解所反演的應力場結果則從地殼深部揭示了地球動力學過程(Pan，2020)。將兩者相結合分析，有助于認識研究區(qū)的應力狀態(tài)以及構造活動。因此，將GPS主壓應變率方向與利用震源機制解所獲取的最大水平應力max進行比較分析，如圖6所示，從而更好地認識研究區(qū)的地殼構造變形特征。由圖6可知，整體上，最大水平應力max方向與GPS主壓應變的方向幾乎一致，呈現(xiàn)出順時針偏轉的現(xiàn)象。為了精細分析圖6反映的各區(qū)域主要特征，將研究區(qū)分為A、B、C、D、E、F 六個區(qū)域。其中，A區(qū)主要包括羌塘地塊的東南部分，B區(qū)主要是印度板塊東北部的阿薩姆地區(qū)，E區(qū)則主要是滇西北區(qū)域，F(xiàn)區(qū)域主要是滇南塊體范圍，C區(qū)域范圍最小，位于羌塘塊體中部。A、E、F三個區(qū)域將川滇菱形塊體的西邊界包含于其中，D區(qū)域包含川滇菱形塊體的東邊界以及塊體內部。

圖6 GPS主壓應變率方向與最大水平應力SHmax方向對比

按照上述分區(qū)，本文統(tǒng)計了該區(qū)域兩種方向的夾角，如圖7所示。其中，角度小于10°的占40%，10°～20°占比25.1%；夾角均值在A區(qū)為39.33°，在B區(qū)為37.35°，在C區(qū)為38.8°，在E、F區(qū)分別為42.63°、28.71°。GPS主壓應變方向與最大水平應力軸方向分布方向一致的D區(qū)夾角均值為5.3°。

圖7 GPS主壓應變方向與最大水平應力SHmax方向夾角統(tǒng)計(a)、分區(qū)夾角均值(b)

由上述應力場反演與統(tǒng)計結果可知，在川滇菱形塊體東邊界以及塊體內部，最大水平應力max優(yōu)勢方向與GPS主壓應變幾乎相同，且由北向南呈現(xiàn)NS-NW的偏轉，表明該處地表淺部與上地殼深部受到了相同機制的驅動力；青藏高原向東南部擠出的作用力在川滇菱形塊體東邊界一直延續(xù)至小江斷裂。而在川滇菱形塊體西邊界，最大水平應力max優(yōu)勢方向與GPS主壓應變展現(xiàn)出不一致性，位于滇西南塊體的紅河斷裂以及瀾滄江斷裂附近區(qū)域也存在此現(xiàn)象，說明在上地殼深、淺部可能具有不相同的驅動機制。

4 討論

本文所獲取的GPS主壓應變方向與max方向分別從地殼淺部與深部提供了認識該區(qū)域構造活動的兩種方式。但應力分布可能受到局部構造環(huán)境或者其它外部因素影響，因此在對結果的分析中需要綜合考慮其它因素。

4.1 川滇菱形塊體西邊界SHmax與GPS主壓應變方向差異分析

楊曉松和馬瑾(2003)研究表明，青藏高原及其周邊地區(qū)表現(xiàn)為地殼內部解耦現(xiàn)象，且在其東南緣該現(xiàn)象更為明顯。印度板塊向歐亞大陸的俯沖作用帶來NNE向的壓力，該作用力與青藏高原東南向的逃逸擠出相互作用(圖8)，可能是造成區(qū)域A、B、E區(qū)差異較大的主要原因。

圖8 喜馬拉雅東構造結及其鄰區(qū)地殼變形模式的示意圖

常利軍等(2015)指出，在青藏高原東構造結內部及其周邊地區(qū)，測量的快波方向分別與最大左旋剪切方向和最大拉伸方向一致；通過GPS應變計算所得的地幔各向異性與測量的快波方向的一致性進一步說明該地區(qū)巖石圈變形屬于垂直連貫變形模式。然而，本文反演結果中A區(qū)GPS主壓應變與max方向存在較大差異，產生這種差異的原因可能有2個：一是震源機制解獲取的地殼應力場是地殼彈性層以上的應力方向，而常利軍等(2015)獲取的快波方向是上地幔與地殼累積效應，觀測應力方向的深部不一致；二是震源機制解的應力場方向是近50 a的結果，而快波方向通常被認為是巖石晶格長達百萬年尺度的構造作用下產生的方向，兩者在時間尺度上存在較大差異。因此，本文結果與通過快波方向得到的結果可能并不是同一種物理量。此外，該分區(qū)的震源機制解數(shù)據(jù)較少，也可能對結果有一定影響。

4.2 滇西南塊體SHmax與GPS主壓應變差異分析

F區(qū)主要包含紅河斷裂以及瀾滄江斷裂兩條主要斷裂帶。GPS觀測顯示該區(qū)地殼主要圍繞喜馬拉雅東構造結旋轉，方向為近SN向或NEE向；地球物理等多方面觀測資料顯示，緬甸微板塊下方呈現(xiàn)高度傾斜的俯沖狀態(tài)，層析成像結果顯示在印緬山脈下方有一個高速異常體向東俯沖(Pesicek，2010)，而同時緬甸微板塊整體隨印度板塊向北移動?？觳ńY果表現(xiàn)出近EW向的形變，該處軟流圈上涌，致使巖石圈減薄，SKS剪切波更多反映的是軟流圈的方向，拖拽巖石圈下部一起運動，導致上下解耦，但速度場剖面結果顯示該區(qū)域沒有成熟的殼內軟弱帶，可能的原因是巖石圈的下沉引起如軟流層拖拽力不強烈。這些使得地殼內耦合，巖石圈與軟流圈解耦(侯強等，2018)，可能是引起F區(qū)中GPS主壓應變與max方向分布雜亂的原因。

4.3 喜馬拉雅東構造結SHmax與GPS主壓應變差異分析

地震各向異性特性可以反映區(qū)域應力場、變形特征以及斷裂活動特征，有助于研究地殼深部構造以及動力學機制。王凱悅等(2021)研究表明，在活動斷裂或其附近快波偏振方向與斷裂走向相一致。在拉薩地塊東部地區(qū)，郭鐵龍和高原(2020)關于剪切波分裂研究結果表明快波偏振方向為NS或NNE，這與本文B區(qū)右上角拉薩地塊內部水平最大主壓應力方向吻合；在察隅斷裂南部，近震直達S波分裂結果表明，快波偏振優(yōu)勢方向與斷裂走向一致，為NW向(黃臣宇等，2021)，與本文B區(qū)左下角最大水平應力方向一致；在印度板塊東北角的阿薩姆地區(qū)，應力反演結果與印度板塊向歐亞板塊NE或NNE向擠壓作用方向一致；而該區(qū)域GPS觀測數(shù)據(jù)較少，這可能是導致B區(qū)max與GPS主壓應變方向存在差異的原因。

5 結論

本文主要利用喜馬拉雅東構造結及其鄰區(qū)≥2.0震源機制解計算得到max方向與GPS主壓應變方向，并對二者進行比較，以研究該區(qū)的變形模式。研究結果表明：

(1)在川滇菱形塊體東邊界以及塊體內部，最大水平應力max方向與GPS主壓應變方向基本一致，且由北向南呈現(xiàn)NS-NW向偏轉，方向變化具有一定的連續(xù)性，說明該處地表淺部與上地殼深部受到了相同機制的驅動力；青藏高原向南部擠出的作用力在川滇菱形塊體東邊界從鮮水河斷裂、安寧河斷裂、則木河斷裂一直延續(xù)至小江斷裂。

(2)在川滇菱形塊體西邊界，max方向與GPS主壓應變呈現(xiàn)出了顯著的差異性，可能由于其構造環(huán)境復雜，深、淺部存在差異運動，在上地殼深、淺部具有不相同的驅動機制；在紅河斷裂附近區(qū)域，max方向與GPS主壓應變方向也存在差異，可能是由于緬甸板塊的俯沖作用力與巖石圈的下沉引起的軟流層拖拽力等的共同結果。

感謝四川地震臺龍鋒高級工程師為本文研究提供的震源機制解數(shù)據(jù)。