2001年昆侖山口西MS8.1地震破裂時空過程的統(tǒng)一模型

董仁東 姚強 施賀青 王琪

1)中國地質大學(武漢)，地球物理與空間信息學院，武漢 430074 2)中國地質大學(武漢)，地球科學學院，武漢 430074 3)中國地震局地震研究所，武漢 430071

0 引言

2001年11月14日9時26分10秒，青海昆侖山口西發(fā)生MS8.1 地震，該地震是自三分量寬頻帶數(shù)字記錄以來，青藏高原內部發(fā)生的最大地震。美國地質調查局(USGS)給出的震中位置為35.946°N，90.541°E，矩震級MW為7.9，該特大地震導致昆侖山斷層西段地表大范圍、高強度破裂。昆侖山斷裂是一條近東西走向、距離長達數(shù)千千米以上的大型左旋走滑斷層，是調節(jié)青藏高原受印度板塊擠壓而造成物質東向滑移的關鍵部位(Tapponnier et al，2001；Wang et al，2001； 王敏等，2020)。歷史上該斷裂大震頻繁，如1963年MS7.0 都蘭地震、1973年MS7.3 和1997年MS7.9 瑪尼地震(Triep et al，1997； 劉剛等，2020)分別位于該地震東西兩側，2001年昆侖山口西地震則填補了昆侖山斷裂帶西段的部分地震空區(qū)。

Xu等(2002)的地質調查結果表明，此次地震破裂長度約400km，產(chǎn)生了最大達7～8m的地表破裂。喬學軍等(2002)首次公布GPS同震觀測結果，其中距斷層最近的點位位移超過1m。任金衛(wèi)等(2005)系統(tǒng)地分析震后GPS觀測資料，給出了該地震同震變形場及震后位移時間序列。Lasserre等(2005)通過獲取橫跨斷層兩側的4幅InSAR影像數(shù)據(jù)，更為準確地圈定地震破裂及影響范圍。綜合上述觀測資料顯示，地震自西向東導致太陽湖、庫塞湖、昆侖山口3段集中破裂區(qū)。

Bouchon等(2003)利用近場區(qū)域寬頻帶數(shù)據(jù)估算地震破裂的平均速度為3.7～3.9km/s，并探測到其最大速度達 5km/s，超過當?shù)丶羟胁ㄋ俣?。Antolik等(2004)利用遠震體波反演地震破裂時空過程，揭示此次地震可分為3個子事件，初次子事件發(fā)生在拉張型的走滑地塹上，后續(xù)兩次為走滑型斷層破裂，總破裂持續(xù)時間約120s； 當破裂傳播到達主斷層時，平均破裂速度達到 3.6km/s。Robinson等(2006)利用遠場SH波資料，顯示破裂主要集中在20km以上脆性地殼內，指出最大破裂速度、最大應力降、最大滑移量之間具有一致性，且昆侖山斷層在西大灘段的分叉導致超剪切破裂迅速終止。

不同學者曾給出一些昆侖山口西地震的破裂模型。近年來的代表性成果有屠泓為等(2016)利用GPS/InSAR反演的破裂分布模型及Wen等(2009)利用遠震地震波與地表破裂資料聯(lián)合反演的破裂時空過程，但迄今尚無綜合所有大地測量及地震觀測資料的聯(lián)合反演模型。鑒于利用不同資料、不同算法構建的模型間存在一定的差異，為了深入、細致地揭示此次地震的發(fā)震過程，本研究融合全球遠震動態(tài)波形、近場GPS和InSAR靜態(tài)位移數(shù)據(jù)，基于有限斷層位錯理論和模擬退火算法開展聯(lián)合反演，構建統(tǒng)一的破裂模型； 同時利用反投影方法，不依賴斷層幾何模型和地球內部結構模型等先驗信息，直接確定地震破裂過程，分析破裂速度變化與有限斷層反演破裂空間分布的內在關系。

1 資料與方法

1.1 有限斷層反演數(shù)據(jù)

本研究使用美國地震學研究聯(lián)合會(IRIS)數(shù)據(jù)中心(1)http：//www.iris.washington.edu/wilber3震中距為30°～90°的遠震體波數(shù)據(jù)，該區(qū)域數(shù)據(jù)無其他波形干擾，依據(jù)信噪比和臺站空間分布最終選取35個P波臺站和33個SH波臺站(圖1)?？紤]到地震持續(xù)時間，本研究拾取P波、SH波初至前30s至初至后130s，共計160s的時程數(shù)據(jù)，通過去除儀器響應和帶通濾波(0.0033～0.5Hz)，獲得用于模型反演的速度波形。

注： (a)區(qū)域構造背景； 紅色五角星表示震中，黃色虛線方框為研究區(qū)域，綠色虛線方框為InSAR覆蓋區(qū)域，藍色和紅色箭頭分別表示同震位移觀測值和模擬值，插圖中灰色實線為研究區(qū)域的斷層分布，沙灘球為昆侖山口西地震和歷史地震空間分布及震源機制解，五角星為震后半年余震分布，最大震級為MW5.6，紅色實線為聯(lián)合反演采用的斷層幾何模型沿地表的跡線； (b)遠震體波臺站分布； 藍色三角形表示P波臺站、紅色三角形表示SH波臺站、黃色三角形表示P波、SH波臺站； (c)淺藍色三角形表示所選用的歐洲區(qū)域臺站； 紫色三角形表示參考站圖1 2001年昆侖山口西地震區(qū)域構造背景和臺站分布

在變形數(shù)據(jù)方面，本研究使用Lasserre等(2005)處理的4幅ERS-1/2衛(wèi)星SAR干涉影像數(shù)據(jù)，從中提取了4468個樣點(樣點最大視線向位移121cm)，完全覆蓋破裂區(qū)域； 另外，從任金衛(wèi)等(2005)公布的GPS位移數(shù)據(jù)中選取34個近場GPS站點。GPS站點雖然對破裂區(qū)域覆蓋遠不及InSAR均勻、完整，但最靠近斷層的XZ02和BS33兩點位于斷層南北兩側，均有約1m的同震位移，方向相反，清楚地顯示了本次地震導致的左旋走滑破裂。由于近場資料時效性較弱，測量時間跨度達數(shù)年，故包含有地震快速震后變形信息。

1.2 斷層幾何形態(tài)及有限斷層反演

本研究幾何模型基于Xu等(2002)的地質調查結果，并按萬永革等(2008)的方式將斷層自西向東分為3段。前兩段主要考慮彎曲的斷層幾何結構及拉張型破裂，第一段位于破裂最西端的太陽湖，長40km，走向為105°； 第二段位于布格達坂，長48km，走向為76°； 第三段為純走滑的主破裂段，長360km，走向為97°。各段的傾角均為81°，寬度均為32km，較萬永革等(2008)的處理略有不同。將各段按8km×4km網(wǎng)格大小分割為子斷層，估算每個子斷層的滑動量、滑動角、上升時間及平均破裂速度(Ji et al，2002、2003)，使用公式為

(1)

為求解各類參數(shù)，本研究采用模擬退火方法，利用前文中的遠震波形和近場變形數(shù)據(jù)，在參數(shù)空間內按一定步長搜尋最優(yōu)值，使得建模目標函數(shù)最小。目標函數(shù)由觀測值的誤差函數(shù)、權重比和模型參數(shù)約束函數(shù)構成。對于靜態(tài)位移數(shù)據(jù)，采用如下誤差函數(shù)

(2)

對于遠震體波數(shù)據(jù)，在小波域內定義如下誤差函數(shù)

errwf=el+eh

(3)

其中，el，eh分別為低、高頻波形數(shù)據(jù)的誤差函數(shù)，具體形式上與靜態(tài)位移誤差函數(shù)類似(Ji et al，2002)。

為使反演結果穩(wěn)定且滿足一定的物理條件，反演計算設定針對模型參數(shù)的約束函數(shù)，強制模型參數(shù)滿足如下條件：①相鄰子斷層上滑移幅度差值應盡量小(光滑約束Scon)；②由滑移量估算的總地震矩逼近某個特定值(震級約束Wcon)，實現(xiàn)參數(shù)解算的正則化。

綜上，動態(tài)、靜態(tài)觀測值誤差函數(shù)及模型參數(shù)約束函數(shù)的建模目標函數(shù)為

err=errwf+Wst*errst+γScon+Wcon

(4)

其中，Wst為動態(tài)波形與靜態(tài)位移數(shù)據(jù)間的權重比；γ為平滑因子，是控制子斷層滑移幅度及其分布的關鍵參數(shù)，一般通過試錯的方式來選取。本研究選取平滑因子時，重點參考地表破裂的地質考察結果，使模型在地表處的預測值與地質結果一致。

將模型的初始破裂點設定在USGS確定的震中位置，模型的總地震矩接近GCMT結果(5.9×1020Nm)。依據(jù)野外地質調查結果，模型的滑動量大小限定在0～9m范圍，搜索步長為30cm； 滑動角范圍為-30°～30°(步長2°)； 平均破裂速度變化范圍為2.0～4.0km/s(步長 0.1km/s)； 上升時間變化范圍為2～10s； 搜索間隔為1s。本研究采用Galvé等(2002)的波速模型計算格林函數(shù)，青藏高原北部的巖石圈厚度在50～65km。相較于P波，SH波對破裂速度具有更高的敏感性，在波形誤差函數(shù)中將P波與SH波的權重比設為1︰1。GPS數(shù)據(jù)點數(shù)較少，與InSAR數(shù)據(jù)點數(shù)相差約100倍，為避免反演模型過分依賴單一資料，故將位移誤差函數(shù)中GPS與InSAR權重比設為100︰1。此外，本研究將目標函數(shù)中Wst設為1，等權對待遠場和近場觀測。

1.3 反投影方法

本研究將反投影方法應用于歐洲區(qū)域的臺網(wǎng)，來求解該地震在P波高頻段的能量分布特征，使用德國地學中心(2)https：//www.gfz-potsdam.de/歐洲區(qū)域臺網(wǎng)的110個寬頻帶垂直分量數(shù)據(jù)(圖1)，將其重采樣至100Hz。臺站震中距在40°～75°之間，此段距離無其他震相對直達P波的干擾； 方位角覆蓋范圍為278°～346°。

反投影方法是將震源區(qū)網(wǎng)格化，將每個網(wǎng)格節(jié)點看作一次可能的破裂子事件，通過疊加波形的高頻信號來求解地震破裂時空分布。該方法不依賴高精度震中位置信息或其他先驗信息(例如經(jīng)驗格林函數(shù)、破裂持續(xù)時間、斷層幾何等特征)，可以用來快速地估算破裂長度、破裂持續(xù)時間及破裂速度等運動學參數(shù)(Ishii et al，2005；Wang et al，2017；Yao et al，2019)。

本研究在震源區(qū)域布置了3600個網(wǎng)格點，格點間距5km，深度30km。采用上述的震中和速度模型計算P波到時，通過波形的互相關運算(時窗長度10s，滑動長度5s，相關系數(shù)閾值0.6)，將波形對齊，并使用1s的滑動時窗疊加波形(總時窗長度200s)，頻率范圍為0.5～2Hz，最終得到該事件的時空展布特征。

2 反演結果

圖2(a)為聯(lián)合遠震和近場永久位移資料反演的昆侖山口西地震同震滑動模型。反演結果表明，盡管昆侖山口西地震具有向西發(fā)展的趨勢，但主要表現(xiàn)單側破裂。如考慮深度15km以上、幅度大于1m的破裂，其累計長度已超過400km。反演得到的地震矩大小為6.1×1020Nm，對應的矩震級為MW7.78，略大于GCMT給出的結果。

注： (a)中白色實線為破裂時間等值線(單位：s)，矢量箭頭大小和方向分別代表塊體破裂大小和方向，走向方向比例尺為深度方向0.5倍，紅色五角星為震中； (c)中紅色虛線為模型地表破裂值，黑色柱體為地質考察結果(Xu et al，2002)圖2 昆侖山口西地震同震破裂模型(a)、地震矩能量釋放函數(shù)(b)和與地質調查結果對比(c)

模型清楚地顯示：在主破裂段上3個最大破裂幅度超過6m的集中破裂區(qū)，其深度不超過10km。其中，最大破裂位于庫塞湖附近，距震中240～280km，滑動量約8m。與地表考察結果對比，在距震中40～150km的區(qū)段模型顯示的破裂幅度(小于1m)低于現(xiàn)場實測和基于光學遙感觀測的結果(Xu et al，2006；Klinger et al，2005)。

模型顯示滑移量大于2m的破裂區(qū)域一般不深于20km。值得注意的是，庫塞湖下方30km處的小塊區(qū)域，其破裂幅度達到3～4m。Mechie等(2004)推斷該區(qū)域地下20km溫度可達700℃，發(fā)生剛性破裂的可能性較小，而萬永革等(2008)和屠泓為等(2016)單獨反演GPS/InSAR資料得到的破裂模型并無此特征，因此其結果可能并不真實可靠，推測是遠震波形中包含地表反射波成分導致的假象。

昆侖山口西地震的破裂持續(xù)時間約130s。前20s發(fā)生雙側破裂，破裂傳播速度約為 2km/s。在第一段上表現(xiàn)為明顯的走滑特征，在第二段拉張型斷層上，釋放能量較小。20s后破裂向東傳播到主斷層，主破裂區(qū)內的平均破裂速度達到 3.5km/s。主破裂區(qū)上的擴展可分為3個階段：①20～55s，破裂的集中程度和地震釋放能量均較弱，平均破裂速度約 3km/s；②55～80s，矩能量釋放速率達到峰值，此過程破裂向東傳播約125km，平均速度約 5km/s，大幅度超過局部剪切波傳播速度(3.2～4.0km/s)，表現(xiàn)為明顯的超剪切破裂，該階段釋放了整個地震約70%的能量，破裂速度與破裂尺度具有一致性； ③80s后，破裂速度和能量釋放速率迅速下降，100s時矩能量釋放速率達到一較小峰值，約120s時，整個地震破裂過程基本結束。該地震發(fā)生后6個月內發(fā)生了5次強余震，參考CMT的定位結果，余震主要集中在庫塞湖東段，與反演模型中的集中破裂區(qū)域位置具有互補關系，證實了反演模型的可靠性。

聯(lián)合模型能較好地擬合遠震波形及近場變形資料。對于遠震波形資料而言(圖3)，除個別位于走向方向上的臺站對破裂速度不敏感，整體擬合效果較好，波形擬合的均方根誤差為0.6。對于GPS數(shù)據(jù)，距離斷層較近的XZ02和BS33兩點擬合效果較好，遠場擬合度下降，模擬得到的同震位移均小于觀測值。InSAR數(shù)據(jù)整體擬合效果較好(圖4)，LOS向位移平均殘差約0.05m，昆侖山口段殘差稍大(0.10～0.15m)。

注： 黑色實線為觀測值； 紅色實線為理論值； 左側數(shù)字分別表示方位角和震中距； 字母表示體波類型和臺站名； 右側數(shù)字表示時程內觀測值峰值圖3 遠震體波擬合情況

圖4 InSAR數(shù)據(jù)擬合情況

本研究同樣嘗試僅依賴大地測量資料得到反演，得到的破裂滑動分布特征與聯(lián)合反演模型類似，但反演得到的地震矩大小為7.9×1020Nm，大于聯(lián)合反演結果6.1×1020Nm。GPS和InSAR數(shù)據(jù)殘差有一定減小。

反投影結果如圖5所示。相比于有限斷層反演結果，反投影從破裂速度變化的角度詮釋了昆侖山口西地震過程的階段性特征。其中，反投影得出最大能量釋放在距震中250km處附近，與有限斷層反演顯示的最大破裂區(qū)完全一致。其區(qū)別僅在超剪切破裂階段，反投影顯示破裂速度可達到 6km/s，大于有限斷層反演得到的 5km/s的破裂速度，且矩能量釋放更為集中，破裂80s后能量釋放速率迅速下降，但破裂速度未明顯減小。

注： 紅色五角星為震中位置； 灰色圓圈為震源區(qū)域網(wǎng)格； 彩色圓圈為反投影得到的地震能量信息； 圓圈半徑大小代表能量相對大?。?顏色代表色標尺所示的破裂時間； 左上角為破裂過程在能量-時間上的投影； 右上角為地震破裂過程在距離-時間上的投影； 圓圈位置代表能量點在時間和距離上的分布； 紅色背景線為速度參考線圖5 反投影地震破裂過程

3 討論與結論

基于大地測量資料的模型反演不能反映震源時間過程信息(Lasserre et al，2005； 萬永革等，2008)，而完全利用遠震體波反演模型的空間分辨率較低(Antolik et al，2004)，難以揭示破裂的細節(jié)特征。大地測量數(shù)據(jù)對斷層淺層部位(深度小于15km)的破裂敏感，對深部破裂的分辨率逐步降低，但遠震波形對深淺破裂的分辨率基本相同。因此，融合遠震波形及近場大地測量數(shù)據(jù)建模不僅將震源時空過程統(tǒng)一表述，同時也提高了深部破裂分布分辨率。由于近場變形觀測數(shù)據(jù)的約束，聯(lián)合反演模型更少受地球結構復雜性和地震波傳播的影響，同時，遠震波形的約束還可將大地測量觀測中可能存在的震后形變信息分離出來，獲得準確的地震矩能量。萬永革等(2008)根據(jù)GPS和InSAR數(shù)據(jù)反演得到的總地震矩為9.3×1020Nm，本研究反演得到的地震矩為6.1×1020Nm，更接近于GCMT的結果，也能解釋部分GPS臺站及InSAR觀測值大于模擬值。

Wen等(2009)以地質調查結果作為約束，僅用遠震體波獲得斷層滑動分布，與本研究結果大體相似，但也有所差別。本文聯(lián)合反演結果顯示在主破裂段西端地表(91°E～91.5°E)未發(fā)生明顯破裂，而Wen等(2009)強制該部位破裂大小與地表觀測值一致。此外，萬永革等(2008)反演顯示該區(qū)域地表破裂為1～2m，相比于Lasserre等(2005)的結果(滑動分量3～4m)，萬永革等(2008)的結果由于GPS數(shù)據(jù)的加入，避免了反演結果的隨意性。Xu等(2006)地質調查結果表明在該區(qū)域地層具有不確定性，地表觀測到斷層位移可能并非實際破裂，很大程度上可能是深處破裂產(chǎn)生的強地震運動導致的地表裂縫，地表測量結果夸大了實際破裂值。

有限斷層反演結果表明，沿斷層走向155～280km之間的破裂速度為 5.0km/s，且滑移量大于5m的破裂集中在該區(qū)域。280km后破裂速度和能量釋放速率迅速減小，根據(jù)地質調查結果，斷層在此處分叉且破裂沿南部的昆侖山口斷層繼續(xù)傳播，余震主要集中在該區(qū)域，相較于主震，余震震級相對較小，結合歷史地震分布，昆侖山口東段仍具有較強的地震危險性。反投影結果的最大破裂速度約 6km/s，破裂穿過該分叉口能量釋放速率減小，但仍以超過 5km/s的速度繼續(xù)傳播到距震中320km處，與Wang等(2016)利用Hi-net得到的結果基本一致。

Vallée等(2008)通過分析破裂速度轉換帶的高頻地震輻射能量，結合地質調查資料，表明在320km處，斷層走向存在7.8°的改變，阻礙了超剪切破裂的傳播。同樣，在距震中150km處，斷層走向存在5.7°的改變，使得破裂達到超剪切速度，表明斷層的幾何結構變化對超剪切破裂具有較大控制作用。由于反投影方法相對于有限斷層反演方法使用頻率更高的信號，對速度變化具有更高的敏感性，且在有限斷層反演中，破裂速度與滑移量大小具有折中關系(Fan et al，2014)，故有限斷層反演模型未能識別出由于走向細微改變而導致的破裂速度快速變化。