魯甸地震發震斷層自發破裂模擬

廖力 李平恩 劉盼 楊建思

摘要：從彈性動力學方程出發模擬了魯甸地震在包谷垴—小河斷裂地震自發破裂過程，探討影響魯甸地震破裂的因素。研究結果表明：魯甸地震的左旋走滑的震源機制以及震級主要是受背景應力場的影響，斷層滑移分布受到斷層幾何結構、水平應力場方向及相對大小的影響，非平面復雜斷層幾何結構是導致魯甸地震復雜的滑動位移分布的原因。

關鍵詞：魯甸地震;自發破裂;曲線有限差分方法;包谷垴—小河斷裂

中圖分類號：P315.332文獻標識碼：A文章編號：1000-0666（2019）01-0024-09

0引言

2014年8月3日云南魯甸發生6.5級地震，震源深度約12km。此次地震震中位于包谷垴—小河斷裂及昭通—魯甸斷裂附近（27.1°N，103.3°E），是青藏高原巴顏喀拉塊體及其鄰近地區最新地震活躍期內發生的一次中等強度地震（鄧起東等，2014）。張勇等（2014）利用地震波形數據反演了魯甸地震破裂過程，結果顯示是一次共軛破裂地震事件，約70%的地震能量釋放在近SN向的包谷垴—小河斷裂上;主震震源機制表明此次地震為一次高傾角左旋走滑地震。張振國等（2014）依據張勇等（2014）反演的破裂過程模擬了魯甸地震的強地面運動并進行了烈度預測。現場地質調查發現魯甸地震的地表破裂均在包谷垴—小河斷裂帶（李西等，2014），顯示此次地震的發震斷裂為包谷垴—小河斷裂，余震精定位（房立華等，2014）及無人機航測技術分析（于江等，2018）的結果也支持這一論斷。因此，本研究擬從彈性動力學方程出發模擬魯甸地震在包谷垴—小河斷裂的地震自發破裂過程，探討分析影響魯甸地震滑動位移分布的主要因素。

1研究方法與模型構建

1.1數值方法

本文利用曲線有限差分方法模擬包谷垴—小河斷裂（圖1）上的魯甸地震的動力學自發破裂過程。該方法采用分裂節點描述斷層間斷面，網

格沿著斷層以及自由表面等不規則界面劃分。通過重新定義曲面斷層分裂節點物理量（Dayetal，2005;Dalguer，Day，2007），曲線坐標系中的彈性動力學方程可以與摩擦準則約束的物理變量牽引力結合起來。該方法既保持了傳統有限差分方法的特點，如計算效率高、容易施加等，又在模擬復雜斷層方面有較高靈活性，可以用來模擬現實中常見的復雜幾何斷層，如傾斜、彎曲、跳躍等斷層（張偉，2006;張振國，2014）。

1.2斷層結構

包谷垴—小河斷裂東起昭通—魯甸斷裂龍頭山段，向西延伸至蓮峰斷裂，平均走向N30°W（李克昌等，1981）。震后地質調查表明，該斷裂斷層產狀為60°∠60°（李西等，2014），斷層長度約40km，圖1給出了包谷垴—小河斷裂的不規則地表跡線示意圖。

現場地質調查只能給出斷層在地表的出露，對于其深部構造則不能給出更多約束信息。斷層的三維復雜幾何結構一般由觀測資料反演獲得，特別是近場強震記錄、空間觀測技術等資料。綜合這些因素，本文構造基礎斷層三維結構，即斷層在地表的走向參考現場地質調查給出的地表跡線進行約束，深部則按照傾角為90°的斷層向下延拓，整個斷層模型沿走向方向長約60km，沿傾向方向寬20km。

現有的觀測資料表明，斷層兩側巖體的接觸面并非簡單平面，而是存在不同空間尺度的粗糙表面，數值計算表明粗糙表面會對自發破裂以及地震波輻射產生重要影響。因此筆者在基礎三維斷層面結構基礎之上附加自相似擾動，擾動尺度約為1km。

地震，尤其是破裂面積較大的地震，其發生一般是由于一定區域達到破裂條件開始破裂并向外傳播。數值模擬地震的破裂過程也需要一個成核區，在成核區內人為設置一個高于斷層強度的剪切應力值，使之處于起始破裂狀態，之后引起斷層面上其他區域的自發破裂。根據地震精定位結果（房立華等，2014），魯甸地震的起始破裂區域發生在斷層模型沿走向約20km處，參考Xie等（2015）的震源機制解及張勇等（2014）的震源破裂過程及精定位結果，震源斷層參數設定為斷層走向165°，傾角90°，震源矩心深度3km，地震自發破裂深度設定為6km，震源成核區在龍頭山附近。

1.3應力狀態與動力學參數

地質觀測資料以及研究區域的中小地震震源機制反演結果顯示研究區域的斷層以走滑斷層為主。基于云南地區地震震源機制分析，孫業君等（2017）給出了該區域內三軸構造應力分布與相對大小。根據其研究結果，本文中包谷垴—小河斷裂帶的應力比設定為：

式中：σH，σh分別是最大、最小水平主壓應力;σv為垂直主壓應力，各應力基本關系表達式為：

式中：ρ和ρw分別表示上覆巖體和水的密度;h表示質點到地表的深度;g是重力常數。斷層及周圍巖體在構造應力的驅動下發生形變及移動，當形變積累的應力能力超過最大強度時，斷層發生相對錯動，釋放應變能，產生地震波輻射。本文選取大尺度范圍內的空間觀測結果作為數值模擬中的應力場方向，并設定最大水平主壓應力σH的方位角為330°，其他兩個主壓應力軸與σH相互垂直，三者構成右手系。圖1給出了2個水平應力的方向和相對大小，其中箭頭對分別表示最大水平主壓應力σH和最小水平主壓應力σh，σv垂直于水平面。

本文用簡單速度弱化準則（Slip-weakeninglaw）作為摩擦準則（Ida，1972;Andrews，1976），臨界滑移距離Dc=0.4m，在小于臨界滑移距離時，靜摩擦系數逐漸減小直至臨界滑移距離轉變為動摩擦系數，也反映了斷層面滑動過程中應力釋放過程（圖2）。劉博研和史保平（2011）利用原地應力測量數據（郭啟良等，2009）和簡單傾斜斷層模型分析汶川地震動力學參數，得到靜摩擦系數為0.5，動摩擦系數為0.2，由于鉆井深度在400m左右，該研究僅限于淺部地震有一定參考意義。本文設定靜摩擦系數μs為0.4，動摩擦系數μd為0.24。

1.4速度結構

包谷垴—小河斷裂地處青藏高原東南緣（鄧起東等，2002），位于印度板塊與歐亞板塊碰撞帶的東部，速度結構復雜，對地震破裂過程以及強地面運動產生重要調制作用。簡單的速度模型不能反映真實地震發生時對近場造成的地面運動和地震災害。本文采用的是區域三維速度結果（鄭晨等，2016），對研究區域介質可以提供較好的約束。

2研究結果

包谷垴—小河斷裂上魯甸地震的模擬主要是非平面上斷層的動力學模擬，即地震觸發之后的傳播、擴展到愈合等一系列的自發破裂過程。為了精細地反映出地震的破裂過程，動力學模擬采用60m網格。

圖3模擬了魯甸地震后15s內的自發破裂過程，本文截取的是從破裂開始到擴展，直至破裂接觸地表6個比較典型的位移分布結果。從圖中可以看出，破裂在斷層上的傳播具有較強的非對稱性，破裂從成核區（五角星）開始向地表傳播，在地表形成了破裂。圖4顯示地震矩能量釋放主要集中在第1～10s，矩震級MW=6.56。自發破裂模擬的震源時間函數和滑動位移分布與張勇等（2014）利用近震、遠震數據聯合反演的結果也較為接近。

從圖5可以看出，魯甸地震自發破裂向深部傳播過程中，由于受到趨于靜水圍壓的三軸應力的作用，斷層深部的有效剪切應力趨于零，滑動在傳播到斷層深部邊界之前就逐漸停止而愈合。破裂在斷層走向兩側是正常傳播，地震的最終滑移量分布（圖5）顯示較強的非均勻性。最主要的是滑動位移從成核區向NW方向傳播直至地表，顯示出較強的左旋走滑性質，這與區域背景應力場的性質關系較為密切。

3討論

3.1斷層幾何彎曲的影響

本文采用的斷層信息是從房立華等（2014）研究圖中取點，以斷層N30°W（李克昌等，1981）插值得到，斷層傾角以60°從地表向下拓展。因此，模型中的斷層面幾乎是一個平面，而實際斷層面常為不規則幾何彎曲結構，從而導致斷層面上更為復雜的應力狀態。為研究斷層幾何彎曲對自發破裂的影響，將模型中斷層向ES方向延拓50km，按照斷層的平均走向在包谷垴—小河斷裂段設置一些幾何彎曲，自發破裂的成核區設置在整個斷層的中部（龍頭山段），保持震源深度、模型的背景應力場等參數設置不變模擬魯甸地震自發破裂過程。從圖6結果可以看到，魯甸地震自發破裂的方式和設置模型較為一致，都是從成核區向NW方向擴展，但破裂釋放的能量相對小一點，震級達到約MW6.35。相比于前述模型能量釋放比較均勻，在此模型中破裂能量釋放更快更集中，破裂的時間主要集中于0～5s，這種自發破裂的特點可能與成核區的斷層幾何彎曲結構有關。圖7則顯示了魯甸地震自發破裂在斷層走向和傾向上的最終滑移量，可以看到，斷層的最終滑移都是發生在走向上，這也顯示了魯甸地震左旋走滑的特點。

3.2背景應力場的影響

3.2.1最大主壓應力方向的影響

孫業君等（2017）給出了該區域內最大主壓應力的方向，但因其采用近于水平主應力中量值較大的主應力方向近似，因此存在最大可能數十度的誤差。因此，本文模擬了將水平最大主壓應力順時針和逆時針分別偏轉20°，即最大主壓應力方向分別為310°，350°的情形，模擬結果如圖8所示。從模擬結果可以發現在最大主壓應力向逆時針旋轉（從310°到350°）的過程中，自發破裂強度逐漸變大，地震震級逐漸變大，在地表的滑動位移也更大。這表明本文模型中逆時針旋轉的構造應力更有利于斷層上的整體破裂。

3.2.2三軸應力場大小相對值的影響

魯甸地震震中區域位于以巧家為中心的三軸應力場大小相對值R從高向低變化的過渡帶（孫業君等，2017），為研究三軸應力場大小相對值對魯甸地震自發破裂的影響，本文還模擬了R值設定為0.8，0.6和0.4的滑動分布。R值越大，在走滑斷層背景應力場表示水平向2個主應力大小越接近。從圖9可以看出：當R=0.8時，無論在地表還是斷層深部，滑動位移相比其他兩種情況的模擬結果都更大，R=0.6時的滑動位移比R=0.4時稍大，震級也稍大。以上說明水平向2個主應力大小越接近時，地震破裂釋放的能量越大。

3.2.3背景應力場的影響

將模型設定的背景應力場按以下方式進行設置，應力比為：

各應力基本關系表達式為：

在此應力場設置下，斷層性質表現為正斷層，其他參數設置保持不變，模擬魯甸地震后0～8s的自發破裂過程，斷層傾向滑動速率結果如圖10所示。從圖中可以看出，魯甸地震自發破裂從震源成核區出發，向震源兩側擴展傳播，直到遇到人工強邊界而停止。因為斷層幾何結構比較均勻，因此破裂顯示出比較明顯的均勻分布特性。斷層破裂在此應力場設置下與前述2個模型的破裂方式完全不同，破裂主要是在斷層傾向上發生，顯示出比較明顯的正斷層性質，斷層的最終滑移分布見圖11。由此可見背景應力場是影響自發破裂方式的重要因素。

3.3滑動摩擦系數的影響

將模型的背景應力場按3.2.3節的方式進行設置，但將動摩擦系數在斷層走向80km以南區域設為0.24，以北區域設為0.27，研究動摩擦系數對滑動位移分布的影響。由模型采用的滑動弱化準則可知，高的動摩擦系數使斷層面上的剪應力更高，使應力降更小，從而使最終斷層滑移更小。將成核區設在沿斷層走向90km處，模擬結果（圖12）與3.2節成核區2側均勻的位移分布不同，成核區以南的斷層滑移更大，而成核區以北因為較大的動摩擦系數導致較小的應力降，使最終滑動位移較小。

4結論

本文以魯甸地震為研究對象，考慮區域三維速度結構差異、區域背景應力場、地震自發破裂斷裂帶幾何結構、建立能反映地表起伏的斷層自發破裂有限差分模型。再以現場地質調查的滑動破裂結果為約束條件模擬魯甸地震自發破裂過程，得到以下結論：

（1）本文模擬的魯甸地震自發破裂過程是基于云南區域的背景應力場，魯甸斷層的幾何資料以及最新三維速度結構建立模型進行模擬，以最新震源機制反演結果的矩心深度以及震源參數設定自發破裂的成核區，模擬的震源時間函數與張勇等（2014）的反演結果較為接近，自發破裂模擬的地表破裂約6km，與魯甸地震后地質調查發現的地表破裂也較為吻合。

（2）研究區域的背景應力場三軸主應力相對大小決定了魯甸地震發震震級、走滑型的震源機制，主應力的方向變化會影響滑動位移分布。

（3）背景應力場對地震自發破裂起決定性作用，而不同的應力場參數包括水平主應力方向的改變以及相對大小的改變，會影響地震破裂的傳播、滑動位錯的分布以及地震波的輻射。

（4）斷層的幾何結構會改變斷層面上應力場的大小，從而控制斷層自發破裂，特別是發震斷層的幾何結構、參數等存在較強的非均勻性時，不同的起始破裂區域往往能導致相差較大的地震。地震起始破裂區域的位置對地震破裂形態有重要作用，可能對相應強地面運動規律與地震災害分布也有重要影響。

（5）在滑動弱化準則控制下的動摩擦系數會影響斷層破裂引起的滑動位移分布，動摩擦系數越小的斷層單元最終滑動位錯越大。

感謝南方科技大學張振國助理教授、中國科技大學張文強博士對自發破裂代碼修改的建議和意見。

參考文獻：

鄧起東，程紹平，馬冀，等.2014.青藏高原地震活動特征及當前地震活動形勢[J].地球物理學報，57（7）：2025-2042.

鄧起東，張培震，冉勇康，等.2002.中國活動構造基本特征[J].中國科學，32（12）：1020-1030.

房立華，吳建平，王未來，等.2014.云南魯甸MS6.5地震余震重定位及其發震構造[J].地震地質，36（4）：1173-1185.

郭啟良，王成虎，馬洪生，等.2009.汶川MS8.0大震前后的水壓致裂原地應力測量[J].地球物理學報，52（5）：1395-1401.

李克昌，候學英，趙維城，等.1981.滇東北地區地震地質特征[J].地震研究，4（1）：57-63.

李西，張建國，謝英情，等.2014.魯甸MS6.5地震地表破壞及其與構造的關系[J].地震地質，36（4）：1280-1291.

劉博研，史保平.2011.MS8.0汶川地震斷層的應力狀態以及對余震危險性的影響[J].地球物理學報，54（4）：1002-1009.

孫業君，趙小艷，黃耘，等.2017.云南地區震源機制及應力場特征[J].地震地質，39（2）：390-407.

于江，張彥琪，李西，等.2018.無人機航測技術在2014年魯甸MS6.5地震震區活動構造調查中的應用[J].地震研究，41（2）：166-172.

張偉.2006.含起伏地形的三維非均勻介質中地震波傳播的有限差分算法以及在強地面運動模擬中的應用[D].北京：北京大學.

張勇，許力生，陳運泰，等.2014.2014年8月3日云南魯甸MW6.1（MS6.5）地震破裂過程[J].地球物理學報，57（9）：3052-3059.

張振國，孫耀充，徐建寬，等.2014.2014年8月3日云南魯甸地震強地面運動初步模擬及烈度預測[J].地球物理學報，57（9）：3038-3041.

張振國.2014.三維非平面斷層破裂動力學研究[D].合肥：中國科學技術大學.

鄭晨，丁志峰，宋曉東.2016.利用面波頻散與接收函數聯合反演青藏高原東南緣地殼上地幔速度結構[J].地球物理學報，59（9）：3223-3236.

AbstractInthisstudy，wesimulatetheprocessofspontaneousruptureoftheLudianearthquakeatBaogunaoXiaohefault，anddiscusstheinfluencingfactorsoftheLudianearthquakerupture.TheresultsshowthatthefocalmechanismoftheleftlateralstrikeslipandthemagnitudeoftheLudianearthquakearemainlyaffectedbythebackgroundstressfield.Thedistributionoffaultslipisaffectedbythechangeoffaultgeometry，directionandrelativesizeofhorizontalstressfield.NonplanarcomplexfaultgeometrystructurecausesthecomplexslidingdisplacementdistributionoftheLudianearthquake.

Keywords：Ludianearthquake;spontaneousrupture;curvefinitedifferencemethod;BaogunaoXiaohefault