海原-六盤山斷裂帶現今地震危險性的數值模擬分析

蔣鋒云，季靈運,2，趙 強

1.中國地震局第二監測中心，陜西 西安 710054；2.防災科技學院，河北 三河 065201

0 引言

地震是在區域構造應力作用下，應力在活動斷裂帶上不斷積累并達到極限狀態，而后突發失穩破裂的結果(Scholz，1990，1998; 張培震等, 2003，2013)，彈性回跳理論表明，當斷層處于震間期時，隨著構造作用力的持續加載，斷層面逐漸閉鎖，斷層兩側地殼運動在斷層面及附近不斷積累應變能，直至發生破裂。因此，對活動斷裂震間期斷層面閉鎖特征，相關區域構造變形與應變積累的研究，是地震研究不可或缺的內容。基于此，許多學者利用GPS觀測得到的地殼水平運動速度場，對海原-六盤山斷裂帶構造變形、應變積累、斷層面滑動分布及閉鎖特征進行了深入的研究。這些研究對認識海原-六盤山斷裂帶現今構造變形特征以及強震發生的危險性，理解青藏高原東北部擴展變形機理具有重要的意義。然而這些研究大多從運動學的角度出發，很少涉及地殼介質屬性和斷層力學性質參數等動力學方面的研究和討論。實際上，區域地殼形變本身受到這些參數的影響和制約，斷裂帶介質力學參數與斷層面上的滑動速率、應力狀態、區域地殼運動速度場等密切相關。如果不考慮這些力學參數對斷裂帶運動的影響，模擬結果可能和真實情況會存在出入。如Mildon et al.(2019)研究表明，庫侖預應力場的存在和斷層復雜結構對地震觸發和強震危險性研究具有不可忽視的作用。也有學者借助有限元技術，顧及地殼介質和斷裂帶力學參數對包括海原-六盤山斷裂帶在內的區域進行形變場的模擬，并討論了斷裂帶力學參數和地殼介質屬性對區域地殼形變場的影響，給出了合理的參數解釋。如He et al.(2013)利用GPS資料，借助有限元技術探討了斷裂帶摩擦系數和地殼介質流變學參數對區域地殼速度場地的影響。研究表明斷層低摩擦系數和恰當的地殼流變學參數是海原-六盤山區域現今地殼運動速度圖像控制的主要因素。石富強等(2018)同樣以GPS資料作為約束，采用三維有限元技術，模擬研究了青藏高原東北緣斷層剪切力學性能(剪切模量)對區域地殼運動速度場圖像的控制作用，進而在最優模型的基礎上分析了當前青藏高原東北緣不同斷裂的應力狀態。結果表明，不同的斷裂其剪切力學性能差異較大，對區域速度場的控制作用也不盡相同，解釋了不同學者給出斷裂摩擦性能存在差異的原因。這些力學模擬，對認識區域地殼介質屬性和斷裂帶力學性質參數具有重要的意義。需要指出的是，由于各自研究側重點的不同，上述研究將某一斷裂帶或斷裂帶的某一段的力學參數(摩擦系數或者剪切模量)由淺至深僅僅作為一個參數，并沒有考慮到沿著深度方向上斷層力學參數的變化。實際上斷層深淺部力學參數可能存在較大的差異，如按照負位錯理論，離開斷層的遠場變形主要由斷層深部無震蠕滑(剪切滑動)來控制，而斷層近場形變主要由斷層淺部斷層震間期運動來控制(部分閉鎖；Reinoza et al., 2015)。

因此，文中將以海原-六盤山斷裂帶附近GPS觀測地殼水平運動速度場作為約束，通過構建三維彈性有限元模型，借助非線性反演方法遺傳算法，參考石富強等(2018a)的文章，將斷層看做正交各向異性介質構成弱化帶，來反演海原-六盤山斷裂帶沿著斷層面平行方向上剪切模量深淺部的差異。進一步結合區域構造特征及地震活動分析海原-六盤山斷裂現今地震危險性。

1 模型構建

海原-六盤山斷裂帶位于青藏高原東北部，該區域處于青藏高原北東向擴展的前沿地帶，新生代構造變形和地震活動十分強烈。海原-六盤山斷裂為該區域內主干斷裂帶，是青藏高原東北部的地貌邊界與構造邊界，也是青藏高原塊體與阿拉善塊體、鄂爾多斯塊體的一級塊體活動邊界帶(鄭文俊等, 2019；李錦軼等，2019)。該斷裂帶自西向東由近東西走向的冷龍嶺斷裂、金強河斷裂、毛毛山-老虎山斷裂、狹義海原斷裂和近南北走向的六盤山斷裂組成。從地震活動性來看，廣義海原斷裂及周邊地區歷史強震非常活躍，曾發生過1920年海原8.5級和1927年古浪8.0級兩次8級以上地震以及多次7級左右地震。兩次8級以上地震之間長達220 km的天祝破裂空段(Gaudemer et al., 1995)有歷史地震記載以來，沒有發生過6級以上地震，存在強震發生的危險。六盤山斷裂帶中段也存在地震空區或者破裂空段，斷裂帶北段及其附近曾發生1219年固原7級地震、1306年固原7級地震和1622年固原北7級地震等；斷裂帶南段是公元600年天水—隴縣間地震的發震斷裂(王師迪等, 2018)；斷裂帶中段的隆德—隴縣之間，歷史上無強震與大地震記載，屬于地震空區或者破裂空段，長約70 km，也具有強震發生的危險性(圖1)。

圖1 研究區構造特征及歷史強震

文中主要研究海原-六盤山弧形斷裂帶斷層深淺部平行斷層面方向剪切模量的差異，因此在模型構建時，為了減少周邊次級斷裂的影響，將模型范圍限制于沿斷裂兩側100 km的范圍內，斷層的展布主要參考徐錫偉等(2017)最新活斷層研究結果，將海原-六盤山斷裂帶近似為連續折線展布。

通過對所調研的地質資料進行假設，構建海原-六盤山斷裂帶三維有限元模型：考慮海原-六盤山地區地殼速度結構(陳九輝等, 2005)，以及地殼孕震層20 km厚深度(段星北, 1997; 張國民等, 2002)，將模型從地面延伸至30 km深。彈性力學參數參考Crust1.0(Laske et al., 2013)以及嘉世旭和張先康(2008)與張洪雙等(2015)給出地殼速度結構，地殼彈性模量采用統一的值為85 GPa，泊松比為0.25。

1.1 有限元斷層模型及網格劃分

斷層的模擬是文中研究的重點。相關研究表明，可以借助斷層兩盤之間的摩擦系數(Hergert and Heidbach, 2010)或者斷裂帶的剪切模量(Fialko et al., 2002)來模擬斷層的力學性能。對應于有限元模擬中，通常采用接觸摩擦模型，或者斷裂帶介質弱化模型來模擬斷層的運動及變形。從長期構造演化來看(萬年甚至百萬年)，斷層作為弱化帶來處理是沒有問題的。同樣，接觸摩擦模型如果來研究長期構造演化最后其實也是連續的滑動變形(時間步長大于地震復發周期)。但對于研究相對短期的(幾十年)震間斷層非連續變形而言，這二者之間就會存在差異。采用介質弱化模型來模擬斷層短期變形行為，會使得斷層內垂直斷層走向力的平衡被破壞，而實際上在長期的構造作用下，斷層與周邊介質已經處于力平衡狀態。因此介質弱化模型使得模擬結果會存在誤差，而相對震間斷層短期變形而言，忽略這種差異可能會使得模擬結果的合理性和可靠性存在問題。而接觸摩擦模型，能夠很好的反映斷層兩盤變形的非連續性，從接觸摩擦理論的角度考慮，也具有明確的物理意義，更接近于實際情況。但摩擦接觸計算是個非線性問題，不僅涉及到計算量大的問題，還涉及到解的收斂問題(李玉江等, 2009)。特別是文中如果考慮斷層面深淺部摩擦系數的不同，不僅使得收斂變得困難，同時將其作為目標參數進行非線性反演，需要的計算量是非常巨大的，使得反演幾乎不可能完成。

石富強等(2018b)在介質弱化帶模型基礎上引入正交各向異性理論，用正交各向異性本構方程替代傳統弱化帶中的線彈性本構方程，實現了彈性模量與剪切模量相互獨立。即將斷層帶看作夾在塊體之間的類似層合板的特殊介質，其在垂直于斷層面方向剪切模量與周邊介質一致，使之不易于發生大的擠壓或拉張變形；而在平行斷層面方向剪切模量相比周邊介質要小，從而能夠沿著斷層面產生較大的滑動變形。這種正交各向異性模型，相比傳統的介質弱化模型，一定程度上克服了由于降低彈性模量使得斷層內垂直斷層走向的力的平衡被破壞的情況，使得模擬結果和摩擦接觸模型更為接近。同時，由于其是連續模型，又能避免摩擦接觸非線性引起的計算結果的收斂和計算量過大的問題。因此，文中采用基于正交各向異性理論的介質弱化模型來模擬斷斷層的變形行為。根據復合材料力學基礎知識，材料正交各向異性斷層本構關系如下，斷層內應力應變關系在其三個主方向上滿足關系：

(1)

其中ε=[ε11ε22ε33ε12ε23ε13]T；σ=[σ11σ22σ33σ12σ23σ13]T分別為三個主方向的應變和應力分量；K柔度矩陣。斷層介質的變形主要表現為走向和傾向的剪切行為，其他方向上的力學性質與周圍介質一致，因此正交各項異性介質本構關系的柔度矩陣可以表述為：

(2)

其中E，v分別為周邊介質彈性模量和泊松比；Gslip為斷層在剪切方向上(走向和傾向)的剪切模量，獨立于周邊介質參數E，v。為了后面反演的需要，引入一無量綱參數λ=Gslip/G0，其中G0=E/2(1+v)為周邊地殼介質剪切模量。斷層的剪切方向上的剪切模量通常小于周邊地殼介質剪切模量，所以λ通常介于0到1之間。

參照利用負位錯進行斷層面滑動分布反演時斷層面的設置情況(李強等, 2014；Li et al., 2016；郝明等, 2017)，將斷層沿走向和傾向，劃分為不同的區，不同的區具有獨立的平行斷層面方向的剪切模量。海原-六盤山斷裂斷層淺部傾角較大，為了建模方便，將其均看作直立斷層。垂直斷層方向剪切模量和周圍介質一致，而斷層泊松比相比周圍介質較大，取為0.28。考慮到文中以斷層為中心構建模型，且模型尺度較小，GPS觀測結果可能更多的反映了地殼淺部的運動，因此對模型進行格網劃分的時候對于20 km以上，深度方向上采取較為密集的網格劃分，每5 km劃分一個網格，20 km以下采用較為稀疏的網格網劃分每10 km劃分一個網格。考慮到計算量和研究目的，在網格劃分的時候，在水平面內也是在斷層附近網格較密，而離開斷層網格較稀疏。整個有限元數值模型采用SOLID186單元進行網格化，共劃分3850個單元，18690個節點(圖2)。

圖2 有限元模型及網格劃分

1.2 模型邊界約束

在數值模擬研究中，模型的邊界條件對最后的結果起著至關重要的作用，合理的邊界條件是模擬能夠取得成功的關鍵因素之一。從平面圖(圖3)上來看，文中模型邊界區域均有一定數量的GPS觀測地表水平運動速度場分布，將其附近GPS速度采用克里金差值方法內插到邊界上作為模型水平向邊界約束。GPS觀測得到的地殼運動速度，是深部地殼運動在地表的直接反應，因此認為在文中給定的模型深度上，地殼運動速度沿深度方向不發生變化。不考慮深部復雜的地球動力學作用，約束模型底部垂向位移為0，水平向不作約束，可以自由移動。參照傳統介質弱化模型斷層模量的取值方法，在正交各向異性斷層模型中，將平行斷層面方向剪切模量和周圍介質的比值的系數作為待確定參數，以模型內部斷層兩側100 km以內GPS觀測速度場作為約束，通過不斷調整比例系數，使得觀測值與擬合值的方差最小，從而得到最優的平行斷層運動方向上斷層的剪切模量的分布。GPS速度場采用Wang and Shen(2020)最新公開發表的1999—2016年長期速度場成果。目標函數如下：

圖3 研究區GPS觀測速度場(1999—2016年)

(3)

uix，vix，uiy，viy分別為x方向的模擬值和觀測值以及y方向的模擬值與觀測值。

1.3 反演方法

目前，基于限元軟件的APDL語言與MATLAB進行交互式樣訪問的技術已經非常成熟(夸克工作室, 2002)，因此，對于文中的有限元非線性反演，可以方便的利用MATLAB提供的非線性反演工具遺傳算法來實現。遺傳算法GA(Genetic Algorithm)是根據生物進化思想而啟發得出的一種全局優化算法，在本質上是一種不依賴具體問題的直接搜索方法。直接以目標函數值作為搜索信息，它僅僅使用適應度函數值來度量個體的優良程度，不涉及目標函數值求導求微分的過程，具有廣泛的適應性。同時，遺傳算法又具有群體搜索的特性，且基于概率規則，而不是確定性規則，具有獲得全局最優解的特性。因此，其在各學科和領域中得到了廣泛的應用。當然，它也有自身的缺點，比如編碼不規范、搜索效率低，計算量大，過早收斂等問題。文中將通過調整遺傳算法多個控制參數，擴大種群數和增加進化代數來獲得全局最優解。通過反演得到的每一代的目標函數的最佳值和平均值圖像(圖4)可以看出經過多達300多次進化，目標函數最小值趨于穩定，且和平均值得差異很小，說明找到了全局最優解。

圖4 目標函數收斂過程

2 模擬結果可靠性評價分析

對于模擬結果的檢驗一般追求的是模擬結果與觀測結果一致，但實際上模擬結果和實際觀測結果不可能完全一致，為了定量的給出模擬結果的可靠程度，石富強等(2018a)對He et al.(2013)的檢驗方法做了一定的補充之后，認為只要當模擬值與觀測值之差位于觀測結果誤差橢圓之內時，模擬結果較為可靠，若超出誤差橢圓，則可信度降低。將第i個觀測點的觀測結果與模擬結果差矢量的模與其對應方位誤差橢圓極半徑之比δi，作為第i個觀測點的評價參數，當δi≤1時，i點的模擬結果落在該點觀測值誤差橢圓之類。具體計算公式如下：

(4)

(5)

(6)

從模擬結果和實際觀測結果的對比來看(圖5)，除了靠近斷裂個別點之外，模擬結果和實際觀測結果在速度大小和方向上基本一致。反映了在周緣穩定阿拉善地塊、鄂爾多斯地塊阻擋下，地殼運動沿著海原六盤山弧形斷裂帶呈順時針旋轉特征，構造上則表現為從海原斷裂帶的左旋走滑向六盤山構造帶的逆沖推覆轉換。從對模擬結果的殘差統計分析上來看(圖6)，接近90%站點的殘差在1 mm以內，且絕大多數對于位于觀測值誤差橢圓之內，整個區域測點δi整體較小，算術平均值小于1，表明模擬結果與觀測結果整體上具有較好的一致性。

圖5 模擬速度場和GPS觀測結果的對比

圖中柱狀灰度深淺代表文中公式(4)中δi的大小

應變作為微小量，不受參考框架的影響，能更好的反映出模擬結果觀測結果之間是否具有一致性。因此，分別利用GPS站點觀測值和模擬值，借助Shen et al.(2015)計算應變的程序，分別使用GPS站點上的模擬值和觀測值，計算了研究區最大剪應變率和主應變率場(圖7)。結果顯示，模擬結果計算的應變率場空間分布特征和觀測結果計算的應變率場基本上一致，只是在量值上略小，可能與模型沒有考慮中下地殼介質粘彈性效應有關。從長期變形來看，介質的粘彈性效應對區域形變場的影響不可忽視，但由于文中模擬的是較短時間震間期地殼運動，且模型尺度較小，模擬殘差在可以接受的范圍內，因此可以近似認為模擬結果是可靠的。從實際觀測值計算的最大剪應變率和主應變率場來看，研究區西段，應變高值區主要集中在毛毛山斷裂、老虎山斷裂和香山-天景山斷裂之間，說明香山-天景山斷裂和金強河斷裂、毛毛山斷裂、老虎山斷裂一起承擔鄂爾多斯地塊和青藏地塊之間在該區域的差異變形。模擬結果計算雖然也顯示了這一特征，但不論是主壓應變率的強度，還是最大剪應變率強度，相比實際觀測計算結果都要小，可能與模型沒有考慮香山-天景山等次級斷裂有關。而考慮多條斷裂更復雜的模型可能會對模擬結果得到改進，但同時會增加反演參數和計算量。狹義海原斷裂帶用模擬結果和實際觀測結果分別計算應變率的差異要比金強河斷裂、毛毛山斷裂、老虎山斷裂小得多，其高值區主要集中于斷裂帶附近，反映狹義海原斷裂帶兩側地殼差異運動明顯。而六盤山斷裂帶利用模擬結果和觀測結果分別計算的應變率場，在空間分布特征、量值上都有較好的一致性。從應變率圖上來看，六盤山斷裂帶在其北段以東平涼附近有一明顯的局部應變高值區，且主張應變和斷裂近乎垂直，經過仔細分析認為，這一現象可能是由于位于平涼市內單個GPS站速度(GSPL)相對周圍測站速度較大引起，與斷裂的整體運動相關性不大。總體來看，六盤山斷裂附近一定范圍內應變不明顯，可能反映六盤山斷裂和汶川地震之前的龍門山斷裂所處的情況類似，由于斷裂高度閉鎖，近場幾乎沒有相對變形。從主應變率來看，整個區域明顯受到近北東向主壓應變率和與垂直方向主張應變率的控制，且自西向東，主應變率方向總體上存在順時針旋轉特征。從主應變率和斷裂的夾角來看，毛毛山斷裂、老虎山斷裂和狹義海原斷裂帶整體上以左旋運動為主，而六盤斷裂帶整體上，特別是相對斷裂中遠場，垂直斷裂顯示輕微的逆沖特征，和地質學上的認識較為一致。

圖中色卡對應顏色表示最大剪應變率大小，十字叉表示主壓、主張應變率大小和方向

3 討論

在一個地震復發周期中，地震發生時，斷層在短時間內發生快速滑動破裂，斷層內阻礙斷層運動的強度快速降低，即對應平行斷層剪切模量快速減小，和周圍介質剪切模量的比值也快速減小。隨著地震離逝時間的增加，斷層面逐漸愈合，斷層力學屬性逐漸增高，即對應平行斷層走向剪切模量逐漸增大，和周圍介質剪切模量的比值也逐漸增大，直到下次地震來臨之前，斷層力學屬性又增加到較強的狀態，即對應平行斷層走向剪切模量幾乎和周圍介質相同，比值接近于1。因此，沿著斷層面方向剪切模量和周邊介質剪切模量的比值，一定程度上能夠反映斷層的閉鎖特征。即比值越小表明斷層越處于蠕滑狀態，越容易滑動，越不容易積累應變能，比值越大表明斷層更可能處于粘結狀態，越不容易滑動，更容易積累應變能。因此可以用沿著斷層面方向剪切模量和周邊介質剪切模量的比值，即上文中所述的λ大小來表示斷層面的閉鎖程度。

從反演結果來看(圖8)，六盤山斷裂中南段閉鎖程度較高，比例系數接近于1(s8和s9段)，深度最深達到30 km(s9段)，而其北段(s7段)和海原斷裂交匯段，僅僅在淺部(0～5 km)存在一定程度的閉鎖。這和利用負位錯方法反演得到的六盤山斷裂帶斷層面閉鎖結果基本一致(李強等，2014; Li et al., 2016；郝明等，2017)，也和地震地質資料研究得出的六盤山斷裂帶中南段處于強震破裂空段具有一致性(M7專項工作組，2012)。整個狹義海原斷裂帶(s4—s6段)閉鎖程度均較低，僅在淺部存在一定程度的閉鎖，比例系數均在0.4以下，可能反映了1925年海原8.5級地震之后，整個斷裂仍然處于震后調整狀態。西段金強河斷裂、毛毛山斷裂、老虎山斷裂淺部閉鎖程度較低，而在5～20 km存在一定程度的閉鎖，表明該斷裂可能存在淺部蠕滑，而較深的位置存在能量積累，具有強震發生的危險性，這和相關研究的結論也基本吻合。

色塊代表每個分塊平行斷層面的剪切模量和周圍介質剪切模量的比值，白色圓圈代表距斷層兩側各10 km范圍內小震(3.0級以上)精定位結果在斷層面上的投影

為了將數值模擬結果和小震精定位結果進行對比，利用中國地震局地球物理研究所房立華提供的小震重新定位目錄，將1980以來發生在研究區內海原-六盤山斷裂兩側各10 km范圍內3.0級以上地震投影到斷層面并疊加到數值模擬結果之上。結果顯示，反演結果和地震活動沿著斷層面空間分布特征具有較好的對應關系：六盤山斷裂帶中南段和北段淺部地震活動較為稀疏，和反演得到的斷層高度閉鎖相對應；海原斷裂帶整體地震活動稀疏，雖然反演結果顯示斷層閉鎖程度較低，但可能與海原斷裂帶1925年8.5級地震使得整條斷裂發生貫通性破裂有關(Gaudemer et al., 1995)；金強河斷裂、毛毛山斷裂(s1，s2段)5～20 km地震活動也較為稀疏，反演顯示該區域具有一定閉鎖程度，而老虎山斷裂(s3段)由地表至深度25 km處，均顯示地震活動較為密集，和其5～20 km存在一定程度的閉鎖區相矛盾，根據老虎山斷裂近場InSAR最新成果(Jolivet et al., 2012)，認為可能是反演所使用的GPS觀測站點較為稀疏，空間分辨率不夠所致。

雖然從模擬效果來看，文中的模擬結果在可接受的范圍內，但從有限元模型構建的合理性、可靠性來看，仍然存在一些需要改進和深入的問題。①考慮更復雜的斷層模型、地殼介質的粘彈性效應和初始應力場等方面的問題，會使模擬結果更好，也會提高模擬結果的科學性和合理性；②對于斷層的模擬，選擇什么樣的模型與所研究的問題有關，如對于地質長時間尺度的模擬而言，傳統的介質弱化模型是合理的，而對于短期的同震、震間或者震后形變而言，可能接觸摩擦模型或基于正交各項異性的介質弱化模型更為合理；③僅以GPS觀測的地殼水平運動速度場作為約束來反演斷層剪切模量，缺乏InSAR、水準等垂向資料的約束，如果加入這些約束資料可能會增加反演可信度；④將斷層統一看做直立斷層，對于垂直于斷層方向上傾滑運動(拉張或擠壓分量)的模擬可能不夠，特別是六盤山斷裂，和龍門山斷裂一樣具有鏟形特征，深部傾角較小。

此外，需要說明的是，相比基于塊體運動的負位錯模型而言，文中不需要考慮塊劃分對模擬結果的影響，但受限于反演方法帶來的計算量和計算效率的問題，使得斷層網格劃分不宜過密，反演結果無法像負位錯反演一樣做到更高的分辨率，因此，文中反演得到的斷層面網格劃分較粗，反演結果不夠精細。

4 結論

文中基于正交各向異性理論的介質弱化模型來模擬斷層的變形行為，將平行斷層面方向的剪切模量和周圍介質剪切模量的比值作為反演參數，以海原-六盤山斷裂附近現今GPS觀測地殼水平運動速度場作為約束，通過構建三維有限元模型，采用遺傳算法，反演了海原-六盤山斷裂平行斷層面的剪切模量分布。進一步結合研區域構造特征和斷層附近地震精定位結果分析了海原-六盤山斷裂帶現今強震危險性。結果表明：六盤山斷裂帶中南段閉鎖程度較高，沿斷層面3.0級以上地震動活動稀疏，未來發生強震的危險性較高，而其北段僅在淺部有一定程度閉鎖，深部3.0級以上地震活動相比中南段更為密集，未來發生強震的可能性較低；狹義海原斷裂帶整體閉鎖程度不高，3.0級以上地震活動較為稀疏，反映在1920年海原8.5級地震之后，該斷裂仍然處于震后調整狀態，推測距下一次強震發生仍較為遙遠；西段金強河斷裂、毛毛山斷裂和老虎山斷裂在淺部存在一定程度的蠕滑，3.0級以上地震活動較為密集，但在5～20 km處具有一定的閉鎖，并可與金強河斷裂帶、毛毛山斷裂帶的小震活動稀疏區域相對應，推測未來具有發生強震的背景。而老虎山斷裂，雖然反演結果顯示在5～20 km同樣具有一定程度的閉鎖，但沿斷層面地震由地表至深部地震活動較為密集，可能處于貫通性蠕滑狀態，未來發生強震的可能性不大。

致謝：中國地震局房立華研究員提供了研究區小震精定位結果，在此表示感謝。