山西南部區域應變場演化與地震關系的數值模擬

李宏偉，劉瑞春，陳永前

(1.山西省地震局，山西 太原 030021；2.太原大陸裂谷動力學國家野外科學觀測研究站，山西 太原 030025)

0 引言

山西中南部至晉陜交界地區歷史上曾發生過多次強震，2014年先后又出現了多項準同步前兆異常，部分異常測項于2015年下半年又同步轉折變化，且與山西南部GPS區域應變參數轉折變化較為一致。從區域應變特征來看，最小主應變于2014年5月出現偏離背景狀態的加速負值變化，表明區域壓性應變在持續增強，與基線的縮短狀態也較為吻合。在此背景下，最大主應變和第一剪應變在2015年5月至12月同時出現短期增大-恢復的變化過程，之后發生了2016年3月12日運城鹽湖ML4.8地震；另外，ML≥3.0地震視應力空間圖像顯示山西南部區域高值異常仍在持續[1-2]。因此需要開展基于地球物理場觀測的具有物理意義的地震預測工作，為山西南部地震趨勢判斷提供科學依據。

地震的孕育是區域應力積累到逐漸加強的過程，當斷裂積累的應變能超過斷裂所能承受的臨界閾值時，就會以地震的形式突然釋放應變能。研究者在實踐中也不斷驗證了這一結論：地震的發生與區域應變能的積累程度有關[3-6]。近年來，隨著GPS觀測精度以及測點密度的提高，逐漸被應用到數值模型的邊界約束條件中。從1999年到現在，先后有劉杰等[7]、馮向東等[8]、劉峽等[9-11]對華北地區進行了三維有限元數值模擬。其中，2007年劉峽[9]利用研究區域的GPS速度場作為模型的邊界約束條件，分析了模型內部的GPS速度場與模擬結果的一致性，反復調整模型材料參數，使兩者盡可能誤差最小，最終得到該區域的應力場特征。2013年劉峽等[11]采用“塊體加載”的思想，模擬了汶川地震前后山西斷陷帶的應力場空間格局，結果顯示，汶川地震后，整個山西斷陷帶呈現構造擠壓增強態勢。隨后，山西斷陷帶于2010年相繼發生河津4.8級地震、大同4.6級地震以及陽曲4.6級地震。2017年李平恩等[12]利用三維粘彈性模型模擬了太原盆地應力場演化與強震的關系，表明盆地的地震活動性受區域應力場控制。目前，山西斷陷帶的地震視應力計算結果顯示，山西斷陷帶南部存在高視應力區[1-2]，因此，該文利用ANSYS有限元軟件，結合研究區地質構造以及斷層幾何分布建立山西南部三維有限元模型，以2011-2016年GPS速度場作為邊界加載條件，分析山西南部近期應變場空間分布特征。

1 山西斷陷帶南部地質構造背景及模型的構建

1.1 研究區地質構造背景

山西斷陷帶位于鄂爾多斯塊體東緣，整體上呈“S”型展布，長約1 200 km，寬約幾十公里，屬于拉張型斷陷。山西斷陷帶南部主要包括臨汾盆地和運城盆地兩個斷陷盆地。其中臨汾盆地由洪洞凹陷和臨汾凹陷組成，受控斷裂主要有蘇堡斷裂、羅云山斷裂以及大陽斷裂。1303年和1695年先后發生在洪洞凹陷和臨汾凹陷的2個歷史大地震，使臨汾盆地劃為高烈度區。運城盆地的主要活動斷裂呈NEE向，西南長邊走向NEE至近EW向，往東北逐漸變為NE、NNE向,盆地由兩凹一壘組成，即侯馬凹陷、運城凹陷和峨嵋臺地。侯馬凹陷受呂梁山山前龍門山斷裂控制，單側斷陷，新生界厚約1 200 m；運城凹陷受中條山北麓斷裂控制，單側斷陷，新生界達5 500 m[13-16]。研究區內主要活動斷裂空間分布及斷裂幾何特征如圖1和表1所示。

圖1 研究區地質構造簡圖Fig.1 Geological structure of the study area

斷裂編號斷層名稱走向傾向傾角(°)運動性質長度kmF1蘇堡斷裂NWWS高角度正斷層30F2羅云山山前斷裂NNESE80°以上右旋兼正斷145F3韓城斷裂NE30°NEE60°左右正斷層100F4峨眉臺地北緣斷裂NE70°N70°～80°正斷層156F5峨眉臺地南緣斷裂NESE70°～80°正斷層50F6中條山斷裂NE-NEENW58°～75°正斷層137

1.2 模型的構建

基于研究區斷裂分布以及GPS站點位置，在經度109°E～113°E、緯度34°N～37°N區域內建立山西斷陷南部三維有限元模型，垂向上包括上地殼、中地殼和下地殼，建模時采用高斯投影變換將研究區內的大地坐標轉化為平面坐標。研究區斷裂帶由很多細小的斷層系、斷層束構成，空間分布繁雜，因此在構建過程中做一定的簡化，即按照斷裂帶的主體走向，將屬于同一斷裂帶的斷層系處理為一條斷層，將斷層作為相對容易變形的軟弱化帶處理，同時適當增加斷層及其附近區域的網格密度，確保之后的模擬計算具有比較好的精度和穩定度。根據研究區地殼速度結構的反演結果，綜合設定上地殼為15 km、中地殼為12 km、下地殼為12 km，斷層帶寬度約5 km左右，斷層深度延伸至中地殼底部。采用有限元靜力學分析方法，假設研究區地殼為彈性結構，具體的材料參數如表2所示，利用ANSYS軟件單元庫中solid185劃分單元，最終建立模型的節點數為375 143個，單元數為2 379 224個(見圖2)。

表2 研究區彈性材料參數表Table 2 Elastic material parameter of the study area

圖2 研究區三維有限元模型Fig.2 Three-dimensional finite element model of the study area

1.3 模型邊界約束條件

一般認為，在幾年時間尺度范圍內，巖石的力學性質主要以彈性為主，區域應力的大小和方向也認為是基本穩定，因此可以將某一時段的GPS速度場作為模型邊界約束條件，得到該時期研究區內應變場增量的空間分布圖。目前由于GPS點位分布較稀疏，同時也存在一些速度不連續的點位，模型邊界約束條件需要通過節點周圍的GPS點進行插值計算，因此需要將誤差較大、速度不連續的點位進行剔除，插值計算過程中盡量避免斷層附近及斷陷帶內部的點，選擇較為穩定的GPS站點。假定速度場在地殼的垂直方向上不隨深度的變化而改變，將經插值計算得到的邊界約束條件加載到模型周緣，同時將模型底部的垂直方向固定，水平方向設為自由面，模型表面不做任何約束，計算時將模型視為無自重彈性體。本次研究采用中國地震局第一監測中心提供的中國大陸GPS水平形變速度場，參考基準為區域無旋轉基準。原始觀測數據來源為“中國地殼運動觀測網絡”和“中國大陸構造環境監測網絡”項目產出的連續和流動觀測數據以及IGS提供的連續觀測數據，第9頁圖3為研究區2011-2016年GPS水平形變速度場。由圖3可知，2011-2016年研究區內變形速率較小，整體運動方向SEE向。該速度場所反應的變形性質與地質研究結果所給出的山西地震帶NW-SE向拉張，NE-SW向壓縮的應力場狀態較為一致。

圖3 研究區內水平形變速度場Fig.3 Horizontal deformation velocity field in the study area

2 計算結果分析

本文所述的應變場為該計算時段的應變場增量，不代表各斷裂實際的應變積累程度，定義應變場結果的正值表示拉張變形，負值表示擠壓變形，剪切應變為水平面方向的剪切應變。

通過GPS速度場給定位移約束，得到2011-2016年時段的主壓形變場和剪切形變場圖像(見圖4)，結果顯示，該時段整體形變量較弱，除斷裂端部、斷裂轉折處應變量約為10-7量級，斷層整體的變化量級約10-8，表明汶川地震對山西斷陷帶的影響逐漸減弱，但區域形變場仍處于弱擠壓形變狀態。在中條山斷裂附近和羅云山山前斷裂附近有相對較高的應變增量集中區，同時剪切形變也具有一致性，羅云山山前斷裂附近的臨汾跨斷層水準也顯示2013年之后斷層張性活動明顯減弱，逐漸轉為壓性變化(見圖5)，與該時段數值模擬計算結果較為一致。2016年鹽湖ML4.8地震發生在中條山斷裂應變增量集中區附近，表明應力場的擾動對地震有促進作用，因此通過分析應力場的相對變化來研究形變場分布與地震的發生具有一定的意義。

圖4 研究區應變場分布圖Fig.4 Distribution of strain field in the study area

圖5 臨汾跨斷層4-1水準測線觀測值Fig.5 Observation of the 4-1 leveling line of Linfen cross-fault

3 結論

在收集研究區斷層幾何參數以及地質材料參數的基礎上，運用ANSYS軟件構建三維地質模型，將GPS觀測結果作為模型邊界約束條件，計算得到研究區的主壓應變場以及水平面方向的剪切應變場空間分布圖像，同時結合2016年鹽湖ML4.8地震進行了驗證分析，初步得到以下結論：

(1) 數值模擬結果顯示，2011-2016年山西斷陷南部區域以擠壓形變為主，跨斷層水準測量結果也在該時段出現張性活動減弱狀態，分析可能受區域應力場調整的影響。

(2) 2016年鹽湖ML4.8地震發生在形變場累積集中區附近，表明利用形變場的相對變化來研究、預測形變場的擾動與地震之間的關系有一定的意義。

(3) 斷層的幾何形態會對應變場空間分布有一定的影響，尤其是在斷層端部、轉折處以及斷裂交匯部位。

(4) 模型計算的應變場表示的是區域應變場的應變增量，斷裂危險性的大小主要取決于斷裂絕對應變能的積累程度，因此將某一期的GPS速度場作為模型的約束條件僅能對本時段各斷裂應變積累的相對量作一定的比較。

衷心感謝中國地震局第一監測中心和中國地震局地震研究所為本次研究提供GPS觀測數據。