茅遠哲 呂國軍 孟立朋

1 河北省地震局，石家莊市槐中路262號，050021 2 河北省地震動力學重點實驗室，河北省三河市學院街465號，065201

地震的孕育和發(fā)生是孕震體應力應變能不斷積累、進入臨界狀態(tài)并最終失穩(wěn)的力學過程，地殼巖石破裂過程中的應力分布、應力擾動與斷裂構造明顯相關[1]。因此，要探討強震的遷移規(guī)律并開展地震危險性分析，最根本的途徑是研究斷裂帶的應力狀態(tài)及其動態(tài)演化過程[2]。

本文首先基于邯鄲地區(qū)活動斷裂地表、淺層、中層、深層探測和小震重新定位等地球物理探測結果以及前人研究成果，建立該區(qū)域三維粘彈性有限元模擬的斷層模型；然后根據(jù)歷史地震和古地震研究結果，加入運動學特征邊界條件，模擬地震發(fā)生在目標區(qū)活動斷裂上引起的應力應變場變化；最后計算1830年磁縣地震引起的周圍斷裂面和滑動方向上產(chǎn)生的庫侖破裂應力變化，給出研究區(qū)活動斷裂地震危險性綜合評價結果。

1 計算方法與模型

1.1 目標區(qū)內(nèi)主要斷裂分布

邯鄲地區(qū)主要斷裂呈北北東向和北西西向分布(圖1，藍色線段表示斷裂，紅色實心圓為磁縣地震)，北北東向的斷裂自西向東為紫山西斷裂、太行山山前斷裂和邯東斷裂(邯鄲隱伏斷裂)，北西西向斷裂自北向南為永年斷裂和磁縣斷裂。這些斷裂近似呈兩兩正交，控制著研究區(qū)內(nèi)的構造活動。研究區(qū)主要活斷層探測結果見表1。

圖1 目標區(qū)主要構造背景

表1 邯鄲地區(qū)主要斷裂要素

1.2 幾何模型的建立及介質(zhì)參數(shù)確定

1.2.1 幾何模型

根據(jù)地表和深部斷裂展布資料，建立目標區(qū)幾何模型(圖2)。為防止邊界出現(xiàn)突變，要求模型區(qū)域比目標區(qū)稍大，模型范圍以紫山西斷裂、太行山山前斷裂、邯東斷裂(邯鄲隱伏斷裂)、永年斷裂和磁縣斷裂為界，每層劃分為6個塊體。因目標區(qū)內(nèi)地殼厚度起伏不大，設計模型每層水平展布，模型各層厚度見表2。

圖2 根據(jù)目標區(qū)斷裂分布得到的平面幾何模型

表2 模型中不同深度塊體的橫波波速

由淺中層地震勘探獲得了主要斷裂的傾角，但在地殼深部這些主要斷裂的展布情況不明，因此在模型中將主要斷裂均設為近直立斷裂。由于目標區(qū)主要在模型中心部分，因此在模型中部對網(wǎng)格進行加密劃分，分成10 563個節(jié)點、8 736個單元(圖3，圖中藍色線段表示斷裂)。

1.2.2 模型介質(zhì)楊氏模量

利用有限元程序TEKTON進行模擬計算[3]，該程序使用的計算方法為根據(jù)巖石上施加的應力控制巖石破裂的機制，將莫爾圓實際半徑(O1-O3)/2與破裂時半徑之間的比值定義巖石極限破裂值Pf，可表示為：

式中，σ1與σ3分別為被施加的最大主應力與最小主應力。計算需要確定模型介質(zhì)的楊氏模量和粘滯性系數(shù)等參數(shù)。

利用深部探測結果來計算模型內(nèi)各不同深度塊體的楊氏模量。由于寬頻帶地震臺陣探測范圍比模型區(qū)域要小，主要在東西向跨太行山山前斷裂帶和邯東斷裂(邯鄲隱伏斷裂)、南北向跨永年斷裂和磁縣斷裂的區(qū)域內(nèi)布設；同時模型的主要關注區(qū)域在模型的中心部分。因此，用圖3中A1、A2、A5、A6塊體的波速代表整個塊體的波速(有一定程度的簡化)。根據(jù)寬頻帶臺陣深部探測結果獲取不同深度模型塊體的橫波波速，見表2。

圖3 有限元模型結構和單元劃分

地震波速與介質(zhì)彈性常數(shù)的關系為：

(1)

E=2μ(1+ν)

(2)

式中，E為楊氏模量，VS為S波速度，ρ為巖石密度，μ為剪切模量，ν為泊松比。根據(jù)地殼巖石的常用量，ν設定為0.25；根據(jù)文獻[4]，ρ取2.7×109g/m3。由式(1)可用VS和ρ求得μ，再根據(jù)式(2)求得E，由此得到有限元模型中不同深度塊體的楊氏模量(表3)。

表3 有限元模型不同深度塊體的楊氏模量

1.2.3 模型介質(zhì)有效粘滯性系數(shù)

周永勝等[5]在總結近30 a來高溫高壓流變實驗資料的基礎上，應用流變數(shù)據(jù)結合地震震源深度分布，對華北地殼流變性質(zhì)進行了研究，得到華北地區(qū)巖石圈強度剖面模型中巖石圈各層的流變強度數(shù)據(jù)。本文根據(jù)其研究成果得出30 km深度范圍內(nèi)的有效粘滯性系數(shù)(表4)。根據(jù)GPS觀測的中國大陸地殼水平位移場計算出的最大剪應變速率值，目標區(qū)的剪切應變率在1.268×10-15/s左右，但對于模型中30～80 km深度范圍，由于不確定其應變速率是否與地表一致，因此應變速率采用平均值1×10-15/s。

表4 有限元模型不同深度的有效粘滯性系數(shù)

1.2.4 模型邊界條件

影響目標區(qū)動力學特征的主要因素為巖石圈介質(zhì)特性和運動學特征，因此模擬時將目標區(qū)的運動學特征作為邊界條件代入有限元模型進行計算。

假設模型底部垂直于地表方向位移固定、地表自由，由于模型只到80 km深度，根據(jù)GPS位移場計算出模型4個邊界節(jié)點上的速度矢量(圖4)，其中目標區(qū)內(nèi)部斷裂標注較多，其GPS點位所示方向與區(qū)外一致，呈SE向。

圖4 華北地區(qū)GPS速度場

2 粘彈性有限元模型的建立和地震危險性分析

利用得到的斷裂活動模型、歷史強震數(shù)據(jù)以及區(qū)內(nèi)巖石圈介質(zhì)分布，可以建立粘彈性有限元模型，模擬地震發(fā)生在目標區(qū)內(nèi)活動斷裂上引起的應力應變場變化，計算活動斷裂的斷裂面和滑動方向上產(chǎn)生的庫侖破裂應力變化，以分析目標區(qū)未來強震危險性。

1830年磁縣斷裂西段發(fā)生7.5級強震(表5)，利用斷層有限元計算中的劈節(jié)點技術[3]，在模擬的彈性步中按表5數(shù)據(jù)施加在磁縣斷裂上。

表5 1830年磁縣地震主要參數(shù)

研究區(qū)內(nèi)紫山西斷裂和磁縣斷裂中、西段現(xiàn)今時有小震發(fā)生，據(jù)此設計模型：磁縣強震發(fā)生，并且紫山西斷裂和磁縣斷裂中、西段小震不斷。

1830年磁縣地震距今已有190 a，而華北地區(qū)強震復發(fā)周期約在千年尺度，據(jù)此模擬的時間過程設計為：1)發(fā)震；2)震后幾天的調(diào)整；3)震后100 a的時間過程；4)震后1 000 a的時間過程，設計模擬的時間步為：1 d×2+2 a×2+10 a×10+100 a×3+1 000 a×3。

參照文獻[6-7]的方法，取μ′= 0.4。數(shù)值實驗表明，改變此值對計算得到的庫侖破裂應力變化的空間分布影響不大，但對其大小有一定的影響。將剪切應力變化投影到假定的滑動方向，與假定的滑動方向一致取正，反之取負。假定的滑動方向取自后續(xù)破裂事件的震源機制。

3 模擬結果與討論

根據(jù)以上模型設計，得到模型的模擬結果。首先根據(jù)時間尺度的不同，選取0 a、4 a、204 a、3 404 a的時間節(jié)點，提取時間步t0(彈性步)、t4(4 a)、t15(204 a)和t20(3 404 a)的位移場進行分析(圖5、6)?？梢钥闯觯趴h地震發(fā)生時(圖5)，由于存在水平和垂向位錯，使得震中附近的水平位移場發(fā)生扭轉(zhuǎn)，相對于遠離震中的位置，目標區(qū)內(nèi)震中的位移量明顯大于其他地區(qū)。

圖5 磁縣地震時和震后4 a模型的位移場

震后4 a，目標區(qū)內(nèi)不同位置位移量的差異減小，但區(qū)內(nèi)位移方向分布沒有明顯變化，說明強震引起的局部位移異常雖然在慢慢減弱，但是還未消除。震后204 a和3 404 a的位移圖像(圖6)表明，整個區(qū)域內(nèi)位移比較均勻，已經(jīng)看不到局部位移方向和大小的異常，說明磁縣地震的影響已經(jīng)被完全吸收。

圖6 磁縣地震震后204 a和3 404 a模型的位移場

圖7(單位MPa)為1830年磁縣地震引起的周圍斷裂面和滑動方向上的庫侖破裂應力變化。圖中藍色至綠色表示斷裂上庫侖應力降低，地震危險性降低；紅色至黃色表示斷裂上庫侖應力增加，地震危險性增加。由于紫山西斷裂和磁縣斷裂中、西段深部小震沿斷裂分布[8]，確定為2條現(xiàn)今活動斷裂，其庫侖應力的變化對于地震危險性評價至關重要。模型顯示，磁縣強震發(fā)生后，磁縣斷裂中、西段庫侖應力加強，在震中附近位置庫侖應力降低。對于紫山西斷裂，在斷裂南端，特別是與磁縣斷裂交接處庫侖應力加強；在斷裂中段，庫侖應力減弱；在斷裂中偏北的一小段斷裂上，庫侖應力增強明顯；在斷裂北段，庫侖應力增量值較小。說明在磁縣斷裂中、西段兩端與2條北東向斷裂交接的部位，地震危險性值得關注。

圖7 磁縣地震在周圍活動斷裂的斷裂面和滑動方向上產(chǎn)生的庫侖破裂應力變化

4 結 語

研究區(qū)內(nèi)，磁縣斷裂中、西段兩端與2條北東向斷裂(紫山西斷裂和太行山山前斷裂)交接處庫侖應力有所加強，未來地震危險性值得關注。

致謝：感謝中國地震局地質(zhì)研究所陶瑋副研究員為本研究提供巨大幫助。