趙 莉, 劉白云, 范 兵, 王文才, 張衛(wèi)東, 劉艷云, 杜建清

(1. 甘肅省敦煌文物保護研究中心, 甘肅 敦煌 736200; 2. 中國鐵塔股份有限公司甘肅省分公司, 甘肅 蘭州 730000;3. 中國地震局地質(zhì)研究所地震動力學(xué)國家重點實驗室, 北京 100029; 4. 甘肅省地震局, 甘肅 蘭州 730000;5. 新疆眾志盛合建設(shè)工程有限公司, 新疆 烏魯木齊 830011)

0 引言

據(jù)中國地震臺網(wǎng)正式測定,2022年1月8日1時45分,在青海省海北州門源縣(37.77°N,101.26°E)發(fā)生M6.9地震,震源深度10 km(https://www.cenc.ac.cn/cenc/dzxx/396391/index.html)。此次地震震中距門源縣城54 km,距西寧市138 km。

震區(qū)位于青藏高原東北緣,是青藏高原往大陸俯沖的前緣部位,它的東北部連接阿拉善塊體,東部毗鄰鄂爾多斯塊體,新構(gòu)造較為強烈。新生代以來,印度板塊持續(xù)北向擠壓,造成青藏高原內(nèi)部地殼縮短增厚,物質(zhì)向周邊擠出[1]。青藏高原東北緣由于受到具有堅硬中上地殼的阿拉善塊體阻擋[2],發(fā)育了一系列地殼尺度的大型斷裂帶。祁連—海原斷裂帶是調(diào)節(jié)青藏高原和阿拉善地塊匯聚的邊界,托萊山斷裂與冷龍嶺斷裂分列于該大型斷裂之西、東,相互為右階斜列銜接,是斷裂帶中西段的重要組成部分[3]。GPS及其相關(guān)研究表明,托萊山斷裂與冷龍嶺斷裂滑動速率差異較大,在主震震中附近閉鎖程度高,應(yīng)力積累顯著[4]。因此,在這兩條斷裂交匯處,變形不均性增大、彈性應(yīng)變能累積,從而導(dǎo)致在該區(qū)域內(nèi)地震活動較為活躍。此次地震是繼2016年1月21日門源M6.4地震后該塊體內(nèi)發(fā)生的又一次強震。震源機制解顯示為走滑型,距離最近的斷層是托勒山北緣斷裂。

震后蘭州大學(xué)與甘肅省科技廳會同中國地震局蘭州地震研究所、中山大學(xué)和青海省地震局等10余名科研人員當天趕赴震區(qū),對本次地震地表破裂帶、發(fā)震構(gòu)造和震害等進行了現(xiàn)場考察。經(jīng)考察后發(fā)現(xiàn),本次M6.9地震地表破裂帶長度超過22 km,主要由2條破裂帶組成。其中,北側(cè)主破裂帶沿廣義海原斷裂帶中段的冷龍嶺斷裂西段分布,從東向西大致經(jīng)過硫磺溝腦分水嶺、蘭新高鐵大梁隧道、道溝至下大圈溝西支溝止,整體走向NWW,長度超過17 km,向兩端逐漸衰減;南側(cè)的次級破裂帶分布在廣義海原斷裂帶中西段的托勒山斷裂東段局部段上,大致沿大西溝、獅子崖一線分布,走向近EW,長度約5 km,構(gòu)成與主破裂帶西段左階斜列的次級分支破裂,應(yīng)為地震向西破裂的端部效應(yīng)(1)資料來源于2022年1月18日科技日報(記者頡滿斌).。

為進一步探明此次地震的發(fā)震構(gòu)造及其幾何形態(tài),本文利用雙差地震定位方法對2022年1月8日青海門源M6.9地震震后幾天的余震進行相對精確的重新定位,然后根據(jù)震區(qū)小震線性分布的疏密程度并考慮地表破裂帶的展布,選擇其中余震相對密集分布的線性矩形區(qū)域,根據(jù)這些矩形區(qū)域內(nèi)的余震震源信息,利用模擬退火及高斯-牛頓的計算機算法擬合本次地震斷層面的參數(shù)。從而基于地震學(xué)的研究方法推斷本次地震的發(fā)震構(gòu)造及在深部的延展形態(tài)。

1 理論方法

1.1 雙差精確定位方法

地震的活動性可揭示斷裂帶的幾何形態(tài)結(jié)構(gòu)[5],地震活動性與構(gòu)造活動結(jié)合是推敲地震活動構(gòu)造較為易于入手的方法,但這種方式的使用卻依靠于地震定位精確度。雙差算法原理是通過各臺站計取的2個相鄰震事件的觀測和理論走時差的殘差來判定它的相對位置,因此可抵消共同誤差產(chǎn)生的不確定因素,尤其是與速度結(jié)構(gòu)特征在橫向變化及與接收有關(guān)的影響,不必進行對臺站校正和精確估計走時就可大幅度改善其相對定位精確程度。目前這種方法已在地震構(gòu)造研究中得到比較廣泛的應(yīng)用[6-8]。

在3D模型中,相對同樣一個觀測站點k,震源i和震源j兩者之間的觀測到時與理論到時之差的方程如下[9]:

(1)

在雙差地震精確定位法中,在震源i和震源j距離較近時,可以認為震源到兩個地震臺站之間的傳播路徑是相同的,此時若走時速度結(jié)構(gòu)是已知的,則式(1)可簡化為:

(2)

(3)

WGm=Wd

(4)

式中:G為M×4N階的偏導(dǎo)數(shù)矩陣式(M為觀測雙差數(shù)目,N為事件數(shù)量);m為4N階震源參數(shù)矢量;W為加權(quán)的對角矩陣表達式;d為雙差數(shù)據(jù)矢量表達式。

在反演計算過程中,對所有余震事件經(jīng)再次定位得到的所有震源參數(shù)(即3個方向上坐標與發(fā)震時刻)加上了使其平均移動為0的約束條件:

(5)

在計算時,一般采用共軛梯度數(shù)學(xué)算法求解上述方程,可以得到帶阻尼的最小二乘解,并利用奇異值分解方法得到有關(guān)模型參數(shù)的分辨度、殘差等。

1.2 斷層面與滑動角計算方法

一般來說中強震發(fā)生后小地震經(jīng)常在斷層面及其附近發(fā)生,這時利用一些小震位置參數(shù)就可以擬合得到斷層面參數(shù)。其原理為:找尋一個平面,使得所有小震的震源位置到這個平面長度的平方和為最小。建立小震到所求斷層面豎向距離與觀測誤差兩者比值平方和為最小的目標函數(shù):

(6)

式中:E為跟ρ、φ、δ三個變量相關(guān)的三元非線性函數(shù)表達式;n表示參與擬合小震個數(shù);δ為擬合斷層面的傾角值;σi是第i個小震定位的誤差。反復(fù)搜索三個參數(shù)ρ、φ、δ的值,E最小時對應(yīng)的φ值就為擬合斷層面的走向。

盡管基于余震疏密情況可計算得到斷層面參數(shù),但滑動角對斷層面之間作用的判斷也是重要的?；诩俣ā獢鄬臃较蚝蜆?gòu)造應(yīng)力場在擬合斷層面上相互作用時剪切應(yīng)力的方向是一致的。萬永革等[10]提供了斷層面滑動角的一般計算方法:

設(shè)應(yīng)力場P軸走向是φp、俯角是δp,T軸的走向是φT、俯角是δT,則在NED的地理坐標系當中方向矢量表達式為:

T=(cosφTcosδT,sinφTcosδT,sinδT)

(7)

P=(cosφpcosδp,sinφpcosδp,sinδp)

(8)

中間軸(B軸)的方向表達式如下:

B=P×T

(9)

主應(yīng)力相對R大小計算過程如下:

(10)

式中:S1、S2、S3分別為最大主壓應(yīng)力、中間主壓應(yīng)力與最小主壓應(yīng)力。

在前述得到的斷層面走向φ、傾角δ基礎(chǔ)上,斷層面滑移方向及法向在地理坐標系下可表達為:

l=(cosφ,sinφ,0)

(11)

m=(-sinφcosδ,cosφcosδ,sinδ)

(12)

n=(sinφsinδ,-cosφsinδ,cosδ)

(13)

式中:l、m、n分別為走向、向下及法向方向。在斷層坐標系l、m和n中,走向和傾向的應(yīng)力量表達式如下:

T1=cos(T,n)cos(T,l)S1+cos(B,n)cos(B,l)S2+
cos(P,n)cos(P,l)S3

(14)

T2=cos(T,n)cos(T,m)S1+cos(B,n)cos(B,m)S2+
cos(P,n)cos(P,m)S3

(15)

此時,滑動角為:

(16)

通常情況下,上述問題的解法有局部和全局算法,但都依賴于初始解,若初始解的選取不恰當,實際解會與正確解相差較遠。萬永革等[10]發(fā)展了模擬退火的全局搜索與高斯-牛頓局部搜索相結(jié)合的算法解決了上述問題,避免了對初始解的過度依賴,既可在全局搜索出最優(yōu)解,還可估計出斷層面參數(shù)的誤差值。另外,在求得斷層面參數(shù)的基礎(chǔ)上進一步利用局部構(gòu)造應(yīng)力場參數(shù)可得到相關(guān)斷層面的滑動角。因而被廣泛應(yīng)用到斷層面參數(shù)的確定研究當中[11-17]。

2 數(shù)據(jù)

本文收集了青海及周邊地震臺網(wǎng)87個地震臺站(其中甘肅地區(qū)7個測震臺、青海地區(qū)80個測震臺)于1月8—13日記錄的余震到時數(shù)據(jù)。地震臺站的分布見圖1。在收集到的地震事件中,原始震中位于研究區(qū)(37.3°～38.0°N,100.7°～102.0°E)余震事件共867次。剔除單臺記錄及震級小于0.5的地震187次,最終選擇的680次余震事件震級分布范圍為ML0.5～5.3(圖2),其中ML0.5～0.9地震1次,ML1.0～1.9地震427次,ML2.0～2.9地震201次,ML3.0～3.9地震41次,ML4.0～4.9地震8次,ML5.0～5.9地震2次,最大余震為1月8日ML5.3地震。被選地震事件總的震相到時包含了Pg震相到時為5 667條(圖3),Sg震相到時為3 604條。

F1:阿爾金斷裂帶;F2:祁連山北緣斷裂;F3:托勒山斷裂;F4:冷龍嶺—天景山斷裂;F5:毛毛山—老虎山—海原斷裂;F6:青海南山—循化南山斷裂;F7:西秦嶺北緣斷裂;F8:東昆侖斷裂;F9:瑪多—甘德斷裂;F10:巴彥喀喇山主峰斷裂;F11:甘孜—玉樹—風(fēng)火山斷裂;F12:黃河—靈武斷裂圖1 震區(qū)主要斷裂及地震臺站分布Fig.1 Spatial distribution of main faults and stations in the earthquake area

圖2 研究區(qū)重新定位前的小地震震中分布圖Fig.2 Epicenter distribution map of small earthquakes in the study area before relocation

圖3 P波走時曲線圖Fig.3 Travel time curve of P-wave

定位結(jié)果既涉及定位方法的選擇,還與采用的初始速度模型有關(guān)。本次在雙差地震定位時所使用的速度模型(圖4)綜合考慮了青藏高原東北緣地區(qū)人工激發(fā)地震和層析成像等研究[18-19],其每層詳細信息列于表1,任意點的P波速度值用插值法提供。

表1 地震重新定位所用的速度模型Table 1 Velocity model for seismic relocation in this paper

圖4 雙差定位時使用的初始1D P波速度結(jié)構(gòu)模型 (據(jù)文獻[18-19]確定)Fig.4 1-D initial P-wave velocity model used for double diff- erence relocation (After reference[18-19])

3 結(jié)果分析

3.1 重定位結(jié)果

利用雙差地震定位法對上述收集并整理好的地震觀測資料進行再次定位。雙差定位程序在讀取震相到時時,由于不同震相分辨難易程度不同,因而讀取精度也有差異,需對不同類型震相資料作加權(quán)。一般來說,初始定位時Pg到時讀取質(zhì)量比Sg高,因此對這兩種類型震相各自賦予1與0.5的比重。考慮到計算范圍內(nèi)余震數(shù)量較多、地震臺站分布較密的特點,臺站和地震間距最大值選取為200 km、地震對之間間距最大值選取為10 km,每個地震最多同時可以與10個地震組成地震對。重新定位后丟棄了47個獨立地震,得到了633個地震事件的精確震源位置。重新定位前后地震走時均方根殘差平均值從0.389 s下降到0.016 9 s(圖5)。余震水平方向的定位誤差值在43～485 m間(圖6)。對比重新定位前(圖2)、后(圖7)余震的震源位置分布,可以看出重新定位后的地震事件相對更加集中地分布于托勒山斷裂東段與冷龍嶺斷裂西段的兩側(cè)。

圖5 定位前后地震走時殘差均方根直方圖Fig.5 Histogram of RMS of seismic travel time residual before and after location

圖6 重新定位后水平方向誤差分布圖Fig.6 Horizontal error distribution after relocation

F1:托勒山斷裂;F2:昌馬—俄博斷裂;F3:祁連山北緣斷裂;F4:冷龍嶺斷裂圖7 研究區(qū)重新定位后的小地震震中分布圖(粗青色矩形框為不同的反演區(qū)域段)Fig.7 Epicenter distribution map of small earthquakes in the study area after relocation (Thick cyan rectangular frame indicates different inversion region)

對震源深度深入研究可為確定地震孕育與發(fā)生的構(gòu)造背景環(huán)境、區(qū)域變形及構(gòu)造動力學(xué)性質(zhì)等方面提供基礎(chǔ)參數(shù)[20]。對比本次研究區(qū)地震事件在重新定位前后震源深度的頻度分布特征(圖8)可以看出,重新定位前本次地震的余震初始震源深度集中分布在5～15 km,重新定位后余震震源在0～20 km深度范圍內(nèi)偏正態(tài)分布,這與Fan等[21]對此次地震序列的震源深度分布特征研究結(jié)論是一致的。

圖8 研究區(qū)小地震重新定位前、后的震源深度統(tǒng)計圖Fig.8 Statistics of focal depths of small earthquakes in the research area before and after location

3.2 斷層面參數(shù)反演及分析

根據(jù)震區(qū)活動斷裂以及本次地震地表破裂帶展布情況,并結(jié)合重新定位后余震沿主干斷裂線性分布特征,本文沿托勒山斷裂東段與冷龍嶺斷裂西段選定了兩個矩形區(qū)域。利用上述兩個矩形區(qū)域內(nèi)的余震事件信息(圖7)進行斷層面參數(shù)擬合計算。綜合模擬退火與高斯-牛頓的算法,求得的托勒山斷裂東段與冷龍嶺斷裂西段斷層面走向值、傾角值及擬合斷層面4個頂點位置列于表2。

表2 擬合2個小段的斷層走向、傾角、滑動角和位置Table 2 Strike, dip angle, slip angle, and position of fault by fitting two small segments

從擬合得到的斷層面參數(shù)值及深部形態(tài)來看,兩個斷層面均比較陡立[圖9(c),圖10(c)],托勒山斷裂東段斷層面走向為95°,即近EW向展布,總長度約為15 km。冷龍嶺斷裂西段斷層面走向為315°,即呈NNW向,總長度約為12 km。它們共同組成了本次地震的發(fā)震構(gòu)造。

圖9 托勒山斷裂東段重新定位后小震的投影分布及距離斷層面情況(圓圈表示再次定位后余震,粗線為擬合斷層面邊界;AA′為邊界端點;SD是走向,DD是傾向,DF是余震與斷層面的之間距離)Fig.9 Projection distribution of small earthquakes in the eastern part of Toleshan fault after relocation and their distance from the fault plane

圖10 冷龍嶺斷層西段再次定位后的小地震分布情況,其余說明和圖9相同F(xiàn)ig.10 Distribution of small earthquakes in the west section of Lenglongling fault after relocation, and other descriptions are the same as those in Fig.9

由于選擇的兩個矩形區(qū)域內(nèi)余震分布相對比較均勻,因此擬合得到的斷層面走向誤差也比較小。另外被選擇的地震事件基本沿斷層走向方向?qū)ΨQ分布,因而得到的余震與擬合的斷層面距離分布頻度[圖9(d)、圖10(d)]也大致對稱,說明絕大多數(shù)余震是分布在所擬合斷層面附近的,且基本以斷層面為幾何中心分別向左右兩個方向相對對稱分布。

3.3 滑動角計算結(jié)果與分析

在研究構(gòu)造應(yīng)力場與震源機制相互間關(guān)系時,不僅需要P軸走向φP和俯角δP,T軸走向φT和俯角δT,此外還需表示主應(yīng)力相對大小的量R。卜玉菲等[22]將甘肅省及附近地區(qū)分割為1°×1°子區(qū),利用這些地區(qū)M3.5以上地震震源機制方面資料反演平均應(yīng)力場參數(shù),得到了上述地區(qū)的應(yīng)力場方向及相對大小的信息。本文從他們的研究結(jié)果中選取能代表本次門源地震震源區(qū)構(gòu)造應(yīng)力場參數(shù)和R值的資料(表3)。

表3 門源地區(qū)的局部應(yīng)力場參數(shù)Table 3 Regional stress field of Menyuan zone

假設(shè)傾伏角誤差值為10°,T軸和P軸方位誤差值為5°,在前文擬合兩段斷層面走向、傾角及其標準差的基礎(chǔ)上,采用表3中的局部應(yīng)力場參數(shù)值及其標準差求解它們的滑動角,得到了每個分段斷層面滑動角(表2)。從結(jié)果來看,托勒山斷裂東段斷層面的滑動角為2.2°,冷龍嶺斷裂西段斷層面的滑動角為9.6°,表明它們的活動性質(zhì)均為左旋走滑。

將擬合得到的斷層面參數(shù)與震后考察的地表破裂帶展布特征做對比,可以看出本次擬合的斷層面與地表破裂帶分布的位置基本一致。從擬合的每段斷層面長度與其相應(yīng)的地表破裂帶長度來看,擬合的冷龍嶺斷裂西段斷層面長度與對應(yīng)的地表破裂帶長度比較一致,擬合的托勒山斷裂東段斷層面長度比相應(yīng)的地表破裂帶長度長,對于兩者相比延伸長度相差較大的原因還需要進一步分析確定。此外,本次擬合的兩段斷層面與現(xiàn)場地質(zhì)考察的發(fā)震斷層滑動性質(zhì)相比也是一致的。

4 討論

大震發(fā)震構(gòu)造研究除傳統(tǒng)方法外,還可以采用余震分布獨立給出發(fā)震構(gòu)造的幾何參數(shù)與滑動性質(zhì)。近些年來,隨著地震測站的密度增加及觀測技術(shù)的進步,地震觀測震相到時拾取越來越精確?；谛≌鹂臻g分布特征分析或者在此基礎(chǔ)上通過某種反演方式來推斷相關(guān)的地質(zhì)構(gòu)造等越來越引起地球科學(xué)家的關(guān)注?；诖舜伍T源M6.9震區(qū)小震資料信息,綜合高斯-牛頓局部搜索法與模擬退火全局搜索法擬合的斷層面與地質(zhì)考察吻合比較好,說明基于余震分布推測大震發(fā)震構(gòu)造是可行的。

從余震重新定位前、后震源深度變化可以看出,重新定位前研究區(qū)內(nèi)余震的震源深度集中分布于5～15 km,重新定位后10～15 km深度范圍內(nèi)的地震頻度明顯減少,總體上表現(xiàn)在0～20 km深度范圍內(nèi)偏正態(tài)分布,且具有震源深度較淺的特征。同時,擬合的斷層面整體埋深深度也在上述深度范圍內(nèi)。這或許是由于青藏塊體的快速抬升和在本區(qū)向NE推擠的過程中遇到阿拉善塊體和鄂爾多斯塊體的阻擋,在這一地區(qū)形成了弧形斷裂組成的斷裂束,而這些斷裂束的新活動性強、埋深淺,故在活動過程中形成了一系列震源深度較淺的地震。

對于擬合的托勒山斷裂東段斷層面長度與對應(yīng)的地表破裂帶長度相比差距較大的原因,可能是因為擬合該斷層面時地震事件的選取考慮了此段整體的余震分布,而通過地質(zhì)考察的地表破裂帶在該斷裂段兩端末端反映并不清晰。除文中開展的此類研究外,還有諸如震源區(qū)地殼速度結(jié)構(gòu)與發(fā)震構(gòu)造之間的關(guān)系等科學(xué)問題有待深入研究。

5 結(jié)論

文中利用雙差定位法對本次門源M6.9地震的余震進行了重新定位,并基于重新定位后的震源信息,利用模擬退火與高斯-牛頓相結(jié)合的算法進行斷層面擬合計算,獲得的主要結(jié)論如下:

(1) 重新定位后大部分余震基本沿托勒山斷裂東段與冷龍嶺斷裂西段分布,這可能意味著本次地震是上述兩條斷裂末端相互擠壓共同破裂形成的。其擠壓力源可能是在NEE向區(qū)域主壓應(yīng)力作用下,祁連塊體在運動過程中受到了東側(cè)鄂爾多斯和北側(cè)阿拉善剛性塊體的阻擋[23],其內(nèi)部的冷龍嶺斷裂西段與托勒山斷裂東段相連的尾端區(qū)域在長期的擠壓下形成應(yīng)力集中區(qū)。當應(yīng)力積累到一定程度后集中釋放時就會造成大震的發(fā)生。

重新定位后的余震分布與地表破裂帶與范莉蘋等[21]對本次余震重新定位等方面的研究結(jié)果是一致的。

(2) 重新定位前余震震源深度主要集中分布于5～15 km范圍內(nèi),重新定位后變化為在0～20 km的深度范圍偏正態(tài)分布。

(3) 擬合得到的托勒山斷裂東段斷層面走向為95°,即近EW向,其長度約為15 km;冷龍嶺斷裂西段斷層面走向為315°,即呈NWW向,其長度約為12 km;從傾角來看,上述兩斷層面均為陡立狀態(tài),符合純走滑斷層的形態(tài)。

(4) 滑動角反演表明,兩個斷層面的運動性質(zhì)均為左旋走滑。

致謝：本研究使用了萬永革研究員提供的斷層面計算程序,在此表示衷心感謝!