青藏高原東南緣現今地殼形變與多尺度應變率場特征

陳 斌 瞿 偉 張 勤 王慶良 郝 明

1 長安大學地質工程與測繪學院，西安市雁塔路126號，710054 2 中國地震局第二監測中心，西安市西影路316號，710054

青藏高原東南緣是我國現今地殼構造形變與地震活動非常強烈的地區之一[1-2]，區域內發育多條深大活動斷裂[3-5]，曾發生過多次中強地震，其中危害最嚴重的是2008年汶川MS8.0強震[6-7]。開展青藏高原東南緣現今地殼構造形變，特別是強震活動影響下的地殼形變與應變場變化研究，對于深入認知青藏高原現今地殼構造活動特征與地震危險性評估具有重要意義。

大范圍尺度下的高精度GNSS觀測資料為研究地殼構造活動與動力學機制提供了有力的基礎參考信息。相關學者基于青藏高原東南緣GNSS監測資料，通過計算庫侖應力變化獲取汶川地震影響范圍[8-9]，并分別利用塊體形變模型[10]、位錯模型[11]、塊體模型結合最小二乘配置模型[12]及塊體應變模型[13]等探討汶川強震對周邊塊體運動特征與龍門山斷裂的影響，為認知汶川強震同震形變及震后短期(1～2 a內)效應對青藏塊體東南緣地殼活動影響提供了重要參考。

本文利用青藏高原東南緣在汶川地震前8 a(1999～2007年)與震后6 a(2011～2017年)的高精度GNSS監測資料，通過構建區域地殼運動與應變多尺度球面小波模型，獲得研究區現今地殼形變與多尺度應變率場特征，分析汶川強震前后不同空間尺度下研究區應變率場差異性特征，并進一步深入探討在何種尺度范圍下模型計算結果是合理的及在何種尺度下模型計算結果最能反映出區域構造形變的真實情況。本文成果可為研究青藏高原東南緣構造活動及其演化規律提供有價值的參考信息。

1 青藏高原東南緣現今地殼運動及變化特征

本文GNSS監測數據來自中國地殼運動觀測網絡與中國大陸構造環境監測網絡。數據處理流程概述如下：采用GAMIT軟件單日松弛解處理基線[12]，采用IGS提供的sp3精密衛星星歷削弱軌道誤差，采用模型改正固體潮和海潮對GNSS臺站的影響，引入中國大陸及周邊IGS站進行聯合解算獲得基線單日松弛解h文件，并將h文件代入GLOBK軟件進行統一平差計算[13]，獲得ITRF2008框架下速度場[14]，并進一步計算獲得研究區域相對于歐亞板塊的速度場(圖1)。需要說明的是，為顧及7.0級以上強震對地殼速度場的影響，采用相應同震模型對2001年可可西里MW8.1[14]、2011年日本宮城MW9.0[15]與2013年蘆山MW7.0地震[16]進行處理。

圖1 青藏高原東南緣GNSS水平速度場Fig.1 GNSS horizontal velocity field in the southeastern margin of the Qinghai-Tibet plateau

圖1(a)和1(b)均顯示，青藏高原東南緣地殼運動整體呈現出順時針旋轉的特征，并在龍門山、鮮水河、安寧河、則木河與小江斷裂帶處形成明顯的速度差異梯度帶，且該速度梯度帶西側GNSS地殼運動速度場量值顯著大于東側。研究區2011～2017年GNSS速度場量值整體大于1999～2007年，1999～2007年研究區EW向平均速率約為6.7 mm/a，NS向平均速率約為-5.8 mm/a；2011～2017年研究區EW向平均速率約為8.0 mm/a，NS向平均速率約為-6.3 mm/a。圖1(a)和1(b)最顯著的差異是汶川強震后龍門山斷裂帶附近西北側地殼運動GNSS速度場量值顯著增大，1999～2007年速率平均值約為8.7 mm/a，2011～2017年速率平均值約為14.3 mm/a；龍門山斷裂帶東南側地殼運動速度量值則變化不大。2011～2017年羌塘、巴顏喀拉與川滇地塊呈現出加速向南東運移并推動華南塊體的趨勢。上述結果表明，汶川強震后研究區內構造塊體地殼運動差異性特征明顯增強。

2 地殼運動與應變多尺度球面小波模型

球面小波技術具有對空間和頻率局部化的功能，將其用于分析地殼運動時，不同的分解尺度會反映出不同空間范圍內的地殼運動和形變特征信息，其中大空間范圍反映的是區域的整體特征，小空間范圍反映的是局部細節特征[17]。本文選用高斯函數差(difference of Gaussians, DOG)構建的球面DOG小波作為球面小波模型基函數[18]，其表達式為：

(1)

(2)

地殼運動監測中，GNSS站點通常分布不均，因此需根據站點疏密程度選擇相應的尺度因子[20-21]。基于多尺度球面小波模型構建原理，根據研究區域范圍的球面積得到區域半徑L，最小尺度因子下球面小波基函數涵蓋的空間范圍應小于2L，不同尺度因子q下球面小波的基函數涵蓋的空間范圍如表1所示。根據本文研究區域GNSS站點分布范圍，可計算得出L≈1 400 km，因此取最小尺度因子q=3；根據滿足建立球面小波模型的基本條件，即最大尺度因子下球面小波模型的每個小波點涵蓋的空間范圍內應至少包括3個GNSS測站點，可確定最大尺度因子q=8。圖2為根據本文研究區域內GNSS站點的空間分布情況確定的不同尺度因子q的情況。

表1 不同尺度因子下球面小波的涵蓋范圍

圖2 GNSS站點空間分布確定的小波模型分解的尺度因子Fig.2 Scale factors of wavelet model decompositionbased on spatial distribution of GNSS stations

相對于歐亞板塊下的區域GNSS速度場包含區域整體旋轉運動，這一整體旋轉運動可通過歐拉旋轉矢量進行模擬[22]。本文在進行地殼應變解算前，先將研究區域GNSS水平速度場進行歐拉旋轉擬合，得到歐拉旋轉角速率為0.276°/Ma，歐拉極為(14.603°S, 89.170°W)，再將兩期GNSS水平速度場減去該歐拉旋轉矢量，得到研究區域無旋轉GNSS水平速度場，并將實測GNSS速度場與模型模擬值進行對比分析(圖3)。從圖3(a)和3(d)可看出，無旋轉GNSS水平速度場可較好地反映出區域地殼相對運動特征，羌塘、川滇塊體保持順時針旋轉運動趨勢，華南塊體整體表現出北北西向運動趨勢，巴顏喀拉塊體整體運動趨勢一致性較差；安寧河-則木河-小江斷裂帶兩側地殼運動GNSS速度場量值在2011～2017年明顯增強；龍門山斷裂兩側地殼運動GNSS速度場在2011～2017年不僅在量值上有顯著增大趨勢，且在運動方向上也表現出明顯的差異性。圖3(b)、3(c)、3(e)、3(f)分別展示了東向與北向速度分量模擬值與實測值殘差統計直方圖。可以看出，絕大部分殘差趨于0且近似服從正態分布，1999～2007年東向與北向模擬值與觀測值之差標準差分別僅為0.63 mm/a和0.64 mm/a，2011～2017年分別僅為0.53 mm/a和0.50 mm/a, 均小于各時間段GNSS速度場觀測誤差。

圖3 青藏高原東南緣無旋轉GNSS水平速度場及模擬值與觀測值殘差直方圖(q=3～7)Fig.3 GNSS horizontal velocity field without rotation in the southeastern margin of the Qinghai-Tibet plateau and histograms of residuals between observation and simulation(q=3-7)

此外，為防止構建的多尺度球面小波模型在確定球面小波函數空間分布時出現過擬合現象，基于研究區域GNSS監測站點分布分別計算1999～2007年與2011～2017年多尺度球面小波模型東向與北向最優正則化參數(圖4)。

圖4 多尺度球面小波模型東向、北向正則化參數Fig.4 The eastward and northward regularization parameters of the multi-scale spherical wavelet model

進一步將無旋轉GNSS水平速度場作為解算球面小波模型應變參數的輸入數據，采用球坐標系計算水平向應變張量D：

(3)

式中，vθ=-vN，vN為速度場北向分量；vλ=vE，vE為速度場東向分量；r為地球半徑；θ為余緯；λ為經度；D12=D21表示矩陣元素對稱。由于在球坐標系中式(3)以南向為正，而GNSS水平速度場常以北向為正，因此通過坐標旋轉分別得到北向、北東向、東北向和東向應變分量εNN、εNE、εEN、εEE，其中，εNE=εEN，進一步可計算出最小主應變率、最大主應變率、最小主應變率方位角、面膨脹率和最大剪應變率[2]。

3 青藏高原東南緣現今地殼多尺度應變率場變化特征

為揭示青藏高原東南緣在汶川強震前后不同尺度下的應變場時空演化特征，利用球面小波模型分別計算獲得研究區在1999～2007年和2011～2017年多尺度應變特征參數。

3.1 主應變率和面膨脹率變化特征

圖5為研究區1999～2007年與2011～2017年在尺度因子q=3～5、6、7、8時的主應變率與面膨脹率分布。可以看出，不同尺度因子下研究區應變特征表現出較顯著的差異，不同空間與尺度范圍內應變積累與變化特征也反映出青藏高原東南緣復雜的地殼構造運動與形變特性。

圖5(a)和5(b)為尺度因子q=3～5時的計算結果，主要反映的是區域整體應變場特征。可以看出，研究區主壓應變方向整體以NW-SE向為主，在川滇塊體南部主壓應變率呈近EW向；羌塘、巴顏喀拉與川滇塊體內部主應變率量值顯著大于華南塊體內部。對比1999～2007年與2011～2017年的結果看出，羌塘與川滇塊體內部呈現出面膨脹特征，而面壓縮特征主要位于巴顏喀拉塊體內，不同的是，2011～2017年巴顏喀拉塊體與華南塊體交界的龍門山斷裂帶處也表現出一定的面壓縮特征。

圖5(c)和5(d)為尺度因子q=6時的計算結果。可以看出，q=6時應變率結果較q=3～5時體現出更多的細節差異，最顯著的特征為1999～2007年主壓應變率主要位于鮮水河(與龍門山斷裂交界處)-安寧河-則木河-小江斷裂帶處，最大量值約-21×10-9/a，方向以NW-SE向為主，且該區域也表現出顯著的面壓縮狀態；2011～2017年主壓應變率雖仍然沿鮮水河-安寧河-則木河-小江斷裂帶處展布，但主壓應變最顯著的區域則位于龍門山斷裂帶及其附近區域，最大量值約-29×10-9/a，同時龍門山斷裂帶也是面壓縮最為顯著的區域，量值高達約30×10-9/a。

圖5(e)和5(f)為尺度因子q=7時的計算結果。可以看出，研究區域應變場特征較q=6時反映出更豐富的局部細節特征。2個時間段內主應變率和面膨脹率整體分布特征較為類似，地殼面膨脹和壓縮區域整體表現出成對輪流出現的特點，但各部分形變特征并不完全一致。該尺度下的應變特征主要反映研究區域內局部連續形變特性，進一步展示出青藏塊體東南緣構造形變的廣泛分布和復雜性。1999～2007年和2011～2017年主應變率與面膨脹率最顯著的差異特征仍主要體現在龍門山斷裂帶附近，1999～2007年龍門山斷裂帶北西側主要表現為輕微面膨脹，東南側則表現出微弱面壓縮狀態；而2011～2017年龍門山斷裂帶兩側均表現出顯著的NW-SE向主壓應變與面壓縮狀態。

圖5(g)和5(h)為尺度因子q=8時的計算結果。可以看出，其所反映的區域主應變和面膨脹分布與尺度因子q=3～5、6、7時差異較大，有些地方甚至出現相反的結果。如當q=8時，2011～2017年龍門山斷裂帶處并未呈現出較顯著的面壓縮(壓應變)特征，反而呈現出一定的面膨脹特征。而根據區域構造活動特性，2008年汶川地震后，在龍門山斷裂帶附近應呈現出較顯著的面壓縮(壓應變)特征。據此，本文舍棄尺度因子q=8，保留合理的最大尺度因子q=7。

從上述分析可看出，當尺度因子q=6時，即可較好地揭示出汶川強震后龍門山斷裂帶附近顯著的主壓應變與面壓縮高值現象。從圖5(c)也可發現，研究區域的主應變中面膨脹和面壓縮區的高值區主要沿川滇塊體的東西邊界交替分布，這一交替區域也是地殼運動速度場顯著差異帶，易造成區域應變積累。而汶川地震后巨大的應力釋放打破了區域原有的應力平衡，在應力再平衡的過程中，研究區面膨脹和壓縮區高值區轉換為沿龍門山-安寧河-則木河-小江斷裂帶對稱分布(圖5(d))，表明在尺度因子q=6時，研究區地殼形變受邊界斷裂帶控制。同時，圖5(e)、5(f)中主應變率與面膨脹率正負值輪流出現的特征也精細地反映出研究區復雜的地殼運動形變特性與地質構造背景、地貌特征；面膨脹與面壓縮相間區域反映出研究區地殼巖性不均勻的特征，可能與區域內廣泛發育的波浪形褶皺有關，反映出了局部區域應力平衡的力學機制。

進一步，為充分發揮出多尺度球面小波模型在表達區域復雜地殼形變特征中的優勢，更加清晰地綜合揭示出區域地殼形變大尺度(整體)與局部形變特征，分別計算尺度因子q=3～5、3～6、3～7下球面小波模型模擬GNSS速度場與實測GNSS速度場殘差的RMSE量值(表2)。表2(單位mm/a)顯示，當尺度因子q=3～7時，模型模擬與實測速度場RMSE量值最小，可較好地綜合揭示出區域地殼大尺度(整體)與局部形變特征，即最佳尺度因子為q=3～7。

表2 不同尺度因子下模型模擬和實測GNSS速度場差值RMSE

圖6為q=3～7時研究區主應變率與面膨脹率分布結果。可以看出，2個時間段中研究區均呈現出以龍門山、安寧河、則木河、小江斷裂帶為界，主應變率總體表現出西強東弱的趨勢；甘孜-玉樹斷裂南段、鮮水河、安寧河、則木河、小金河-麗江、紅河、金沙江和小江斷裂帶處是主應變率與面膨脹率高值區；甘孜-玉樹斷裂南段、鮮水河、安寧河、則木河、金沙江斷裂帶處均呈現出NW-SE向擠壓應變與NE-SW向拉張應變的特征；小金河-麗江、紅河、小江斷裂帶處呈現出近EW向擠壓應變與近SN向拉張應變的特征。

圖6 尺度因子q=3～7下青藏高原東南緣主應變率與面膨脹率分布Fig.6 Principal strain and dilatation rate with scalefactor q=3-7 in the southeastern margin of the Qinghai-Tibet plateau

圖6(a)顯示，1999～2007年研究區內最大擠壓應變率高值區位于甘孜-玉樹斷裂南段、鮮水河、安寧河、則木河和小江斷裂帶處，量值約為-70×10-9/a。小金河-麗江與紅河斷裂帶處擠壓應變率平均量值約為-25×10-9/a。面壓縮特征也主要沿上述幾條主干斷裂帶展布。巴顏喀拉和華南塊體內主應變率與面膨脹率數值相對均較小。龍門山斷裂帶整體主應變率量值較小，在龍門山斷裂帶最南端與鮮水河及安寧河斷裂帶交界處呈現出一定NW-SE向擠壓應變與面壓縮狀態。圖6(b)顯示，2011～2017年研究區內整體主應變率與面膨脹率量值較1999～2007年均有所增加，特別是沿著甘孜-玉樹、鮮水河、安寧河、則木河、龍門山、小江、小金河-麗江與紅河斷裂帶處尤為突出。2011～2017年研究區域內最大擠壓應變率高值區位于龍門山斷裂帶及其附近區域，斷裂帶區NW-SE向擠壓應變率量值高達約-91×10-9/a，該區域也是整個研究區域內面壓縮最大量值區，面壓縮量值高達-72×10-9/a。甘孜-玉樹斷裂南段、鮮水河、安寧河、則木河和小江斷裂帶處呈現出更顯著的擠壓應變特征，擠壓應變量值平均約為-45×10-9/a。巴顏喀拉塊體東部與華南塊體西部也呈現出較強地震前顯著的擠壓應變特征。

3.2 最大剪應變率變化特征及與地震活動性關系

圖7為研究區1999～2007年和2011～2017年在尺度因子q=3～5、6、7時的最大剪應變率分布特征。可以看出，研究區最大剪應變率在不同尺度因子下雖呈現出細節差異，但最大剪應變率高值區幾乎都沿主干斷裂展布。

圖7 不同尺度因子下青藏高原東南緣最大剪應變率分布Fig.7 Maximum shear strain rate with differentscale factors in the southeastern margin of theQinghai-Tibet plateau

圖7(a)和7(b)為研究區在尺度因子q=3～5時的計算結果，反映的是大尺度范圍內區域最大剪應變率整體分布特征。可以看出，最大剪應變率高值區在整個研究區內未呈現出較多的細節性特征，但在1999～2007年區域內最大剪應變率最大值位于鮮水河斷裂帶北東段與甘孜-玉樹斷裂帶交界處，量值約26×10-9/a；而2011～2017年最大值則位于鮮水河與安寧河以及龍門山斷裂帶南段的交界處，量值約29×10-9/a。

圖7(c)和7(d)為研究區在尺度因子q=6時的計算結果。可以看出，最大剪應變率較q=3～5時的結果有所不同，最顯著的特征為相較于1999～2007年，2011～2017年時段龍門山斷裂帶呈現出明顯的最大剪應變率高值區，量值約17×10-9/a；在巴顏喀拉塊體東部以及華南塊體西部靠近龍門山斷裂帶處也呈現出較明顯的最大剪應變率高值區，量值約15×10-9/a。上述特征表明，汶川強震后龍門山斷裂帶及其附近區域構造剪切運動顯著增強。

圖7(e)和7(f)為研究區在尺度因子q=7時的計算結果。可以看出，1999～2007年與2011～2017年研究區內主干斷裂帶剪切活動在汶川強震后有顯著的增強趨勢。

圖8為青藏高原東南緣1999～2007年與2011～2017年在尺度因子q=3～7時的最大剪應變率分布特征。可以看出，研究區最大剪應變率高值區主要分布于甘孜-玉樹斷裂南段、鮮水河、安寧河、則木河、小江、小金河-麗江與紅河斷裂帶及其附近區域。表明主干斷裂及其附近區域是青藏塊體東南緣地殼構造活動最為劇烈的區域。同時，對比2個時間段研究區發生的MS>4.0地震空間分布特征也可看出，上述主干斷裂及其附近區域也是地震活動的高頻帶。

圖8 尺度因子q=3～7下青藏高原東南緣最大剪應變率分布及與地震活動疊加圖Fig.8 Maximum shear strain rate with scale factor q=3-7 and its superposition with seismicity in the southeastern margin of the Qinghai-Tibet plateau

如圖8(a)所示，1999～2007年研究區最大剪應變率最大值位于甘孜-玉樹斷裂帶與鮮水河斷裂帶交界處，量值約79×10-9/a；除研究區域內幾條主干斷裂處呈現出最大剪應變率高值區，在巴顏喀拉地塊中部也呈現出部分高量值區，但在龍門山斷裂帶及其附近區域最大剪應變量值很小。如圖8(b)所示，2011～2017年研究區主干斷裂及其附近區域最大剪應變率量值大小和范圍整體均較1999～2007年有所增加。雖然該時間段內最大剪應變最大值仍位于甘孜-玉樹斷裂帶與鮮水河斷裂帶交界處及鮮水河斷裂帶南段，量值高達72×10-9/a，但相較于1999～2007年該時間段最顯著的差異特征是龍門山斷裂帶及其附近呈現出最大剪應變率高值區，量值達約50×10-9/a。同時，羌塘、巴顏喀拉、川滇甚至華南地塊西部均呈現出較1999～2007年更顯著的地殼活動性。甘孜-玉樹斷裂南段、鮮水河、安寧河、則木河、小江、小金河-麗江與紅河斷裂帶活動性進一步增強，此外，小金河-麗江與紅河斷裂帶南部區域地殼活動性也進一步增強。從地震空間活動性來看，汶川強震對青藏高原東南緣區域內主干斷裂帶及其附近區域的地震活動性影響最為顯著。

相比于1999～2007年，2011～2017年研究區沿著鮮水河、安寧河、則木河、小江、小金河-麗江與金沙江斷裂帶及其附近MS>4.0地震明顯增多，特別是龍門山斷裂帶及其附近區域地震沿斷裂走向密集展布，且在華南塊體西部與龍門山、安寧河與則木河斷裂帶交界處地震活動也明顯增多。2個時間段地震活動的顯著差異性進一步揭示出汶川強震對于研究區，特別是震中附近地殼構造活動性的顯著影響。

圖9為1999～2007年與2011～2017年尺度因子q=3～7下青藏高原東南緣的面膨脹率和最大剪應變率差值分布。從圖9(a)看出，在距離汶川地震震中較近的龍門山、安寧河、則木河斷裂帶及華南塊體西側呈現顯著的面壓縮狀態，反映汶川強震后巴顏喀拉地塊與川滇塊體加速向南東運移并推動華南塊體的運動趨勢；在其他主干斷裂帶附近也呈現出面應變率的高值區。圖9(b)也顯示出整個研究區內的主干斷裂帶附近是最大剪切應變率的高值區所在，但最顯著的特征還是在汶川地震震中的龍門山斷裂帶與鮮水河斷裂南段呈現出的明顯高值區。圖9所示的這些顯著的應變率差異特征進一步表明，汶川強震對于研究區，特別是龍門山斷裂帶附近地殼構造活動性具有顯著影響。

圖9 尺度因子q=3～7下青藏高原東南緣1999～2007年和2011～2017年面膨脹率和最大剪切應變率差值Fig.9 Differences of dilatation rate and maximum shear strain rate with scale factor q=3-7 from 1999 to 2007 and from 2011 to 2017 in the southeastern margin of the Qinghai-Tibet plateau

4 結 語

1)汶川強震前后研究區地殼運動形變具有一定的繼承性發展特征，地殼運動整體呈現出向南東運移同時自身作順時針旋轉運動的特征，在龍門山、安寧河、則木河與小江斷裂帶處形成明顯的速度差異梯度帶，且高應變率值也主要位于主干構造斷裂及其附近區域。

2)不同尺度因子能夠揭示出研究區域在不同空間范圍內的地殼應變場積累特征。本文依據滿足建立球面小波模型的基本條件，首先初步定量確定出最小與最大尺度因子，再根據不同尺度因子計算結果與區域構造活動特性的符合程度確定出合理的最大尺度因子，進一步通過計算不同組合尺度因子下模型模擬與實測地殼速度場殘差的RMSE大小，最終確定能夠較好地綜合揭示出區域地殼形變大尺度(整體)與局部形變特征的最佳組合尺度因子。結果表明，對于青藏高原東南緣區域，合理的最大尺度因子q=7；當尺度因子q=6時，可較好地揭示出區域整體構造活動特性，即清晰地揭示出汶川強震后龍門山斷裂帶處呈現出的主壓應變率、面壓縮率與最大剪應變率高值特征；當尺度因子q=3～7時(最佳組合尺度因子)，可較好地綜合揭示出區域地殼大尺度(整體)形變與局部形變特征。

需要說明的是，基于GNSS觀測獲得的應變率結果可較好地反映出區域現今構造應力場特征，但本文聚焦于分析青藏塊體東南緣在不同尺度上對汶川強震發生前后的較長期響應特征，即使細致地劃分尺度因子，對研究區域實際地殼形變復雜性的表征還是有限；且考慮到區域構造與形變是一個長期的地球動力學過程，GNSS觀測也僅反映了現今地殼形變活動。此外，本文主要依據地殼活動程度分析地震活動與最大剪切應變率的空間對應關系，但地震的發生和發展極為復雜，特別是對于汶川強震震后效應的細致研究，還需結合震后粘彈性松弛等機理模型進行深入分析。

致謝：感謝野外測量工作者獲取GNSS觀測數據，感謝中國地震局第二監測中心提供高精度GNSS數據。