基于有限斷層震源模型的北京地區(qū)強地震動模擬*
——以1679年三河—平谷8級地震為例

巴振寧，趙靖軒，梁建文，張郁山，張玉潔

(1.天津大學(xué) 中國地震局地震工程綜合模擬與城鄉(xiāng)抗震韌性重點實驗室，天津 300350；2.天津大學(xué) 土木工程系，天津 300350；3.中國地震災(zāi)害防御中心，北京 100029)

0 引言

近年來全球發(fā)生的破壞性地震，如1994年美國Northridge地震、1995年日本Kobe地震、1999年土耳其Izmit地震、1999年中國臺灣Chi-Chi地震、2003年伊朗Bam地震，均在近斷層區(qū)域造成了大量的人員傷亡和建筑結(jié)構(gòu)破壞(趙由佳等，2018)。其中Kobe地震造成了6 400多人死亡和接近1 000億美元的經(jīng)濟損失，工業(yè)重鎮(zhèn)神戶和大阪遭受慘重破壞(Pitarka，1998)。近斷層區(qū)域的嚴重震害促使研究人員開始重視這一區(qū)域地震動的研究，相關(guān)學(xué)者采用有限元(Hisada，1998)、有限差分(Rodgers，2019)、邊界元(劉中憲等，2021)、譜元(Komatitsch，2012；Stupazzini，2009)等數(shù)值方法模擬了近斷層地震動作用下復(fù)雜場地地面運動，試圖預(yù)測未來發(fā)震斷層附近的強地震動分布，為城市規(guī)劃和結(jié)構(gòu)抗震以及震害預(yù)防提供參考依據(jù)(李祥秀等，2021)。

北京作為中國的政治、經(jīng)濟和文化中心，其地位的重要性毋庸置疑。另一方面，北京坐落于典型盆地構(gòu)造的范圍內(nèi)，具有發(fā)生大地震的構(gòu)造背景。震害研究表明，沉積地形對強地震動有顯著的放大效應(yīng)(高文學(xué)，馬瑾，1993)，從而對自振周期較長的高層建筑、大跨橋梁等建筑物造成比較嚴重的破壞(巴振寧等，2020)，給北京地區(qū)的經(jīng)濟發(fā)展和人民生命財產(chǎn)安全帶來潛在威脅。

1679年三河—平谷8級地震是北京地區(qū)有歷史記錄以來最大的一次地震。據(jù)史料記載，此次地震由于震級較大、波及范圍廣、余震持續(xù)時間長，造成了大量的人員傷亡和財產(chǎn)損失(劉中憲等，2019)。其發(fā)震斷裂(夏墊斷裂)距北京市中心最近距離不到50 km，對北京地區(qū)建筑結(jié)構(gòu)產(chǎn)生很大威脅(余中元等，2020)。眾多學(xué)者對1679年三河—平谷地震進行了大量的模擬研究：高孟潭等(2002)采用有限差分方法，模擬的震源最高頻率為1.6 Hz，震級為8.0，并著重研究了北京盆地對地震動的放大效應(yīng)；朱耿尚(2014)采用有限差分法，模擬的震源最高頻率為1 Hz，震級為7.7，并研究了沉積層分布對強地面運動的影響；付長華(2012)采用有限差分法，通過改變斷層的埋深、傾向和滑動角等因素模擬了震級為7.5時的3～10 s長周期地震動，并研究了不同斷層震源參數(shù)對盆地放大效應(yīng)影響的差異；潘波等(2009)采用有限元法模擬得到的震源最高頻率為1 Hz，震級為8.0，并發(fā)現(xiàn)三河—平谷地區(qū)出現(xiàn)比較明顯的近斷層效應(yīng)和盆地效應(yīng)。

但值得注意的是，1679年三河—平谷地震作為歷史大震缺乏地震記錄，難以進行震源破裂的反演，斷層面上破裂細節(jié)尚不清楚(劉博研等，2007)。在模擬近斷層地震動時，由于距離發(fā)震斷層比較近，地面運動強烈依賴于發(fā)震斷層面的位錯發(fā)展過程、滑動方向和破裂速度等斷層面上的破裂細節(jié)，斷層面上破裂方式的改變影響了整個斷層面上的破裂過程，從而對近斷層地震動作用下地面運動分布有著不可忽視的影響(劉啟方，2005；張效亮等，2018)。但在進行北京地區(qū)地震動模擬研究時，尚未針對斷層破裂方式這一重要的震源條件對地面運動分布造成的影響進行分析，因此可能對強地面運動的空間分布以及峰值大小估計不足，從而影響城市規(guī)劃和重大工程抗震設(shè)防。

基于以上分析，本文選取北京地區(qū)(39°30′～40°30′N，116°10′～ 117°20′E)作為研究對象，建立包含起伏地形的地下介質(zhì)三維速度結(jié)構(gòu)模型，并參考已有研究給出的斷層位置、面積和傾角等震源基本參數(shù)建立確定性有限斷層震源模型，開展了最高頻率為2 Hz的確定性物理模擬。筆者在國家超算中心“天河一號”超級計算機上采用譜元法模擬了1679年三河—平谷8級地震作用下北京地區(qū)強地面運動分布特征，深入研究不同斷層破裂方式(單側(cè)破裂和雙側(cè)破裂)對北京地區(qū)強地震動的影響和可能引起的地震危險性分布特征，以期為近場復(fù)雜場地地震動估計和工程抗震設(shè)計等提供借鑒和參考。

1 確定性物理模型建立與參數(shù)設(shè)置

確定性物理模型的近斷層地震動模擬精度和有效性取決于精細的三維速度結(jié)構(gòu)模型(包含復(fù)雜場地條件)和合理的震源模型(張冬麗等，2009，劉啟方等，2008)。為此，本文建立了含起伏地形的地下介質(zhì)三維速度結(jié)構(gòu)模型，通過地面高程數(shù)據(jù)和地下勘探資料，確定起伏地表及其下地層的空間分布。震源模型采用運動學(xué)震源模型，相關(guān)參數(shù)包括斷層的形狀、大小、位置等全局參數(shù)，以及斷層面上的位錯分布、破裂速度和破裂傳播方式等局部參數(shù)。

1.1 運動學(xué)有限斷層震源模型

運動學(xué)有限斷層震源模型考慮為多個動態(tài)子源(點震源)地震反應(yīng)疊加。建立有限斷層震源模型時主要涉及斷層面子源的劃分以及斷層震源參數(shù)的設(shè)定。主要思路是將斷層面上每一塊區(qū)域假定為由一個地震矩點源控制的子源。子源群中存在一個最開始破裂的起始點，其余每個子源的破裂起始時間由與起始破裂點的距離、傳播方式以及時間過程確定。通過將所有子源的矩張量大小和對應(yīng)的時間函數(shù)進行疊加，模擬整個斷層面上的破裂過程。依據(jù)上述思路，筆者基于Matlab針對SPECFEM3D譜元程序開發(fā)了運動學(xué)有限斷層震源模型程序，解決了譜元模擬中有限斷層震源的輸入問題，Matlab程序流程見圖1。

圖1 有限斷層震源設(shè)置主要步驟Fig.1 The main steps of finite fault source settings

針對三河—平谷地震建立確定性震源模型時，筆者首先根據(jù)歷史考察資料進行全局參數(shù)的設(shè)置(表1)，然后使用高度約為1 km的三角形劃分斷層面，共劃分為2 846個子源斷層，進而進行局部震源參數(shù)設(shè)置，震源模型參考付長華(2012)提出的震源模型，斷層面上設(shè)置兩個凹凸體，其相對位置如圖2所示。破裂面總面積為1 280 km，取22%為凹凸體的面積，取最大凹凸體面積為200 km，次級凹凸體面積為76 km。凹凸體中的平均滑移量為7 m，背景域的滑移量為2.3 m。并針對單、雙側(cè)兩種破裂方式，共設(shè)置了3個震源模型，分別把震源破裂的起始點設(shè)置于斷層面的中心和斷層面左右邊界的中心震源。

表1 斷層震源參數(shù)設(shè)置Tab.1 Source parameter settings of the fault

(a)斷層面雙側(cè)破裂模型

(b)斷層面單側(cè)破裂模型(破裂起始點位于斷層左側(cè))

(c)斷層面單側(cè)破裂模型(破裂起始點位于斷層右側(cè))圖2 不同破裂方式對應(yīng)的運動學(xué)有限斷層模型Fig.2 Kinematic finite fault models corresponding to different fracture modes

1.2 三維速度結(jié)構(gòu)模型

本文借鑒付長華等(2015)選擇的區(qū)域范圍(39°30′～40°30′N，116°10′～117°20′E)作為研究對象，根據(jù)地面高程數(shù)據(jù)、物探及鉆孔資料，通過Csimoft公司開發(fā)的前處理軟件Trelis建立含起伏地形的北京地區(qū)精細三維速度結(jié)構(gòu)模型，模型東西長約120 km，南北寬約100 km，沿縱向深度約為40 km。模型包含地殼和上地幔頂部，被6個速度分層界面分割(圖3a)，各個分層介質(zhì)的物理參數(shù)見表2。

表2 北京地區(qū)各層介質(zhì)的物理參數(shù)Tab.2 Physical parameters of the strata in Beijing area

建模時首先將地面高程數(shù)據(jù)和各個分層界面控制點的經(jīng)緯度轉(zhuǎn)化為平面直角坐標，然后依次將各個界面的控制點導(dǎo)入到Trelis建模軟件中并將控制點掃掠生成各個分層界面，最后自上而下疊加所有分層界面，建立包含起伏地形的三維物理模型，如圖3b所示?；谧V元法模擬時需要通過劃分網(wǎng)格對整體三維模型進行離散，若將網(wǎng)格能傳播的最大頻率擴大一倍，則需要將網(wǎng)格間距尺寸減半，并且時間步長也隨之減半，這就導(dǎo)致對目標地區(qū)的地震動模擬的計算量和計算時間顯著增加。同時，為考慮地層實際情況，在進行三維模型建模時選擇包含地殼和上地幔頂部等6個速度分層界面，模型的每個分層界面均為起伏面(圖3)，這也為三維整體模型進行網(wǎng)格離散帶來一定的困難。本文通過進一步減小網(wǎng)格尺寸，將設(shè)定模型網(wǎng)格能模擬的最大頻率從1 Hz擴展為2 Hz。譜元法中為保證結(jié)果的精確可靠，要求最短波長中至少包含5個GLL積分點，根據(jù)實際地質(zhì)介質(zhì)參數(shù)，模型自地表至3 km深度處網(wǎng)格大小約為300 m，3 km至底部40 km深度處網(wǎng)格大小約為900 m，劃分網(wǎng)格總數(shù)約270萬，GLL節(jié)點數(shù)量約為2.23億。

圖3 北京地區(qū)地質(zhì)構(gòu)造界面(a)和三維速度結(jié)構(gòu)模型(b)Fig.3 Geological interfaces(a)and 3D velocity structure model(b)in Beijing area

2 模擬結(jié)果及分析

計算時依次給出不同破裂方式下速度波場快照，峰值速度()分布，加速度、速度和位移的空間分布情況。計算任務(wù)基于國家超算中心“天河一號”計算平臺，使用譜元程序SPECFEM3D，共調(diào)用200個進程計算，時間步距取為0.002 5 s，模擬90 s內(nèi)的地震波傳播，每種工況計算時間約4 h。

2.1 速度波場快照及其分析

圖4給出了3種震源破裂方式下地面東西方向速度波場快照，圖中五角星位置代表震源的起始破裂點，箭頭指向為斷層破裂方向。3種破裂方式波場快照圖均從第5 s開始，并每隔10 s輸出相應(yīng)的波場快照結(jié)果。

從圖4可以看出，不同破裂方式對應(yīng)的地震波傳播過程存在明顯差異。雙側(cè)破裂的破裂起始點位于斷層中間，波場表現(xiàn)為地震波自破裂起始位置向兩端擴散的傳播過程；單側(cè)破裂的破裂起始點分別位于斷層端部，波場相應(yīng)表現(xiàn)出沿單方向傳播的形態(tài)。對于3種震源破裂方式，震源破裂前方的地震動顯著高于破裂后方，近斷層的方向性效應(yīng)十分明顯。上述速度波場結(jié)果一定程度上證明了模擬結(jié)果的正確性和合理性。

圖4 雙側(cè)破裂(a)、單側(cè)破裂(左)(b)、單側(cè)破裂(右)(c)3種震源破裂方式下速度波場快照Fig.4 Snapshots of velocity wave field in bilateral rupture(a)，unilateral rupture(left)(b)， and unilateral rupture(right)(c)modes

2.2 PGV云圖分布及其分析

圖5給出了不同破裂方式下地面水平方向分布圖。圖中顯示3種情況下震中區(qū)域水平向處于2.3～2.5 m/s，地面分布情況與潘波等(2009)計算的北京地區(qū)的水平分布較為相符。強地面運動峰值分布的總體特征呈現(xiàn)“南強北弱，東強西弱”的趨勢，斷層跡線附近出現(xiàn)條帶狀集中分布特征。強地震動集中區(qū)的形狀受到破裂方式的影響出現(xiàn)明顯差別，雙側(cè)破裂引起的強地面運動水平分量在斷層兩側(cè)的分布形態(tài)更為狹長，單側(cè)破裂下分布形態(tài)相對更為集中。與單側(cè)破裂相比，雙側(cè)破裂產(chǎn)生的強地面運動分布范圍更廣，模擬出的高烈度區(qū)范圍明顯擴大，尤其北京市區(qū)、大興地區(qū)在雙側(cè)破裂情況下，水平可達到1.0 m/s，明顯高于單側(cè)破裂情況。

圖5 雙側(cè)破裂(a)、單側(cè)破裂(左)(b)、單側(cè)破裂(右)(c)3種震源破裂方式下PGV分布Fig.5 PGV distribution in bilateral(a)，unilateral rupture(left)(b)， and unilateral rupture(right)(c)fault modes

上述現(xiàn)象的出現(xiàn)主要由于不同斷層破裂方式下，破裂起始點和破裂方向的改變產(chǎn)生的斷層破裂方向性效應(yīng)與深厚沉積地形對地震動相互作用機制的結(jié)果。由于北京地區(qū)第四系主要發(fā)育于平谷、大廠和順義等幾個凹陷處，其中順義和大廠凹陷第四系厚度為800 m，平谷凹陷的第四系厚度也超過400 m，不同破裂方式下地震動傳播過程中受到沉積地形的影響，產(chǎn)生了地震能量遷移的現(xiàn)象，導(dǎo)致地震動響應(yīng)被明顯放大，強地震動的分布范圍以及地震動峰值大小出現(xiàn)明顯差異。

2.3 地表地震動時程及其分析

為進一步研究不同震源破裂方式對地震動空間分布的影響，選取北京市區(qū)、大廠等11個觀測點，給出了不同震源破裂方式下各個觀測點水平方向和豎直方向的加速度、速度和位移結(jié)果(圖6)以及各個觀測點對應(yīng)的峰值結(jié)果(表3)。

表3 不同破裂方式下不同觀測點的PGA、PGV、PGDTab.3 PGA，PGV and PGD values corresponding to fault-rupture modes at different observational sites

2.3.1 近斷層地震動效應(yīng)

本文模擬結(jié)果很好地反應(yīng)了近斷層地震動特性，以雙側(cè)破裂情況為例，各個觀測點地震動時程結(jié)果(圖6)表明：①近斷層集中性效應(yīng)顯著，強地震動主要集中在沿斷層兩側(cè)一個狹窄的范圍內(nèi)，其中平谷地區(qū)水平和垂直向最大分別為0.38和0.39 g，三河地區(qū)水平和豎直向最大分別為0.37和0.35 g，隨著斷層距的增加，觀測點的迅速衰減，昌平地區(qū)的約為0.02 g；②近斷層速度大脈沖現(xiàn)象顯著，靠近斷層位置的三河、平谷等地區(qū)均出現(xiàn)了速度峰值大、持時短的速度脈沖現(xiàn)象，其中平谷地區(qū)水平和豎直向的達到1.98和1.55 m/s，三河地區(qū)水平和豎直方向的達到1.63和1.53 m/s；③斷層破裂造成地面的永久變形現(xiàn)象顯著，在靠近斷層的位置變形也相對明顯，以平谷地區(qū)為例，模擬中水平達到1.03 m，永久位移達到0.45 m。

2.3.2 不同破裂方式對地震動時程的影響

從圖6還發(fā)現(xiàn)斷層破裂方式顯著影響了觀測點的地震動時程結(jié)果，尤其是在靠近斷層觀測點出現(xiàn)了明顯的差異。以大廠觀測點的地震動時程結(jié)果為例，單側(cè)破裂(破裂起始點位于斷層右側(cè))時峰值出現(xiàn)的時刻最早并且水平和豎直向分別為0.34和0.33 g，高于雙側(cè)破裂時的0.32和0.31 g以及單側(cè)破裂(破裂起始點位于斷層左側(cè))的0.31和0.28 g；對比不同破裂方式的時程結(jié)果發(fā)現(xiàn)，斷層破裂方式對速度時程結(jié)果的影響最為明顯，以大廠觀測點處水平和豎直向速度為例(圖7)，破裂起始點位于斷層右側(cè)的單側(cè)破裂情況下，大廠地區(qū)的水平和垂直向分別為2.10和1.79 m/s，地震動持續(xù)時間為45 s，地震動峰值和持時均明顯高于雙側(cè)破裂情況和破裂起始點位于左側(cè)的單側(cè)破裂情況。上述現(xiàn)象在斷層距較近的觀測點，如三河、平谷和通縣等地比較明顯，在斷層距較遠的位置如懷柔、昌平等地不同破裂方式產(chǎn)生的地震動時程則差異不大。

圖6 不同觀測點對應(yīng)的加速度(a)、速度(b)、位移(c)時程結(jié)果Fig.6 Time histories of acceleration(a)，velocity(b)，and displacement(c) corresponding to observational points

圖7 不同破裂方式下大廠觀測點水平(a)、豎直(b)方向速度Fig.7 Velocities in horizontal(a)and vertical(b)direction in fault rupture modes at Dachang observational site

出現(xiàn)上述情況的原因可能是：①震源破裂方式的不同導(dǎo)致各個觀測點的震中距不同。當(dāng)震源的破裂速度一定時，觀測點震中距越近，地震波傳播至各個觀測點的時間越早，出現(xiàn)峰值時刻的時間越早。從圖6觀察發(fā)現(xiàn)，震中距越小的觀測點如三河、平谷和大廠等，在不同破裂方式下峰值時刻差異明顯；②震源破裂方式影響了地震動峰值大小和地震動持時，由于震中距的差異，導(dǎo)致近場上覆土層對于近斷層地震動的衰減效果不同。當(dāng)觀測點距離震源中心較近時，松軟土體對地震動的衰減效果不明顯，導(dǎo)致該破裂方式下地震動峰值較高。特別是當(dāng)區(qū)域內(nèi)存在較厚沉積地形時，由于峰值時刻出現(xiàn)時間早、峰值較大且沉積地形對地震動的聚集放大效應(yīng)，導(dǎo)致地震動衰減的時間明顯延長。通過上述分析可以判定，斷層震源的破裂方式與近斷層區(qū)域內(nèi)強地面運動密切相關(guān)。

3 結(jié)論

本文采用譜元法對1679年三河—平谷8級地震進行了模擬，通過改變斷層面的不同破裂方式(即單側(cè)破裂和雙側(cè)破裂)，得到了不同破裂方式下地面的速度波場快照、分布圖和地震動的空間分布結(jié)果，通過對這些結(jié)果進行分析得出以下主要結(jié)論：

(1)北京地區(qū)強地面運動分布體現(xiàn)了條帶狀集中性分布特征、方向性效應(yīng)以及永久位移等近斷層地震動特性。北京地區(qū)強地面運動峰值速度分布呈現(xiàn)“南強北弱，東強西弱”的趨勢。三河—平谷地區(qū)最高可達0.37 g，大廠地區(qū)最高可達0.32 g，表明平谷、大廠等地區(qū)為震害分布的主要區(qū)域。

(2)斷層破裂方式顯著影響了北京地區(qū)的地面運動分布情況，根據(jù)分布圖發(fā)現(xiàn)雙側(cè)破裂比單側(cè)破裂影響的范圍更廣，北京地區(qū)遭受地震災(zāi)害的威脅最大。

(3)不同斷層破裂方式對近斷層區(qū)域內(nèi)的地震動時程結(jié)果也存在較大影響。以大廠觀測點為例，由于震中距不同，在不同破裂方式下地震動時程表現(xiàn)出的峰值時刻、峰值大小和地震動持時明顯不同，尤其當(dāng)觀測點處存在沉積地形時該現(xiàn)象更為明顯。

本文通過精細化地質(zhì)模型建模以及不同震源模型的輸入，模擬得到了北京地區(qū)的強地面運動分布，上述模擬結(jié)果可以為北京地區(qū)震害預(yù)測預(yù)防和建設(shè)規(guī)劃提供一定的指導(dǎo)和借鑒。

感謝中國地震局地球物理研究所付長華博士提供的北京地區(qū)地質(zhì)勘探資料，感謝國家超級計算天津中心“天河一號”為本文計算提供的幫助。