基于隨機振動模型的強地面運動數值模擬及其在廣西靈山M6地震中的應用

陳剛 李茂峰 李克華 唐勇 張清

馬桂芳2) 申文豪3,4) 姜文亮3)

1)中國南方電網有限責任公司超高壓輸電公司百色局, 廣西百色 533000

2)廣西壯族自治區地震局, 南寧 530022

3)應急管理部國家自然災害防治研究院, 北京 100085

4)上海佘山地球物理國家野外科學觀測研究站, 上海 201602

0 引言

地震災害的根源主要是斷層突然運動產生地表位移,激發地震波造成地表強烈振動,即強地面運動。強地面運動可直接導致建、構筑物的倒塌、破壞,也可引起大量的崩塌、滑坡、滾石、地裂縫、砂土液化等,從而間接導致基礎設施等遭受嚴重破壞。強地面運動研究的基礎是強震記錄數據,但是由于地震的隨機性和不確定性,強震記錄往往集中在某一些地區,地理環境和地質條件的不同使得許多數據的應用價值受到限制。為了在強震數據匱乏地區獲得滿足特定地理地質條件的地震數據,需要開展強地面運動模擬研究工作。

廣西是我國少震、弱震地區,其所在的華南地塊是中國大陸一個相對穩定的構造單元,地震活動性較弱(黃強強等,2022)。雖然廣西地震活動性水平較低,但是在歷史上卻發生過華南內陸地區自有地震記載以來最大的地震——1936年靈山M6地震。1936年4月1日9時31分,廣西靈山縣平山圩東羅陽山發生M6強烈地震,震中位于靈山縣東北約20km的羅陽山西北麓一帶,震中烈度達Ⅸ度或Ⅸ度強(李冰溯等,2018)。該地震造成101人死亡,263人受傷,8000余間房屋倒塌破壞(陳恩民等,1984; 李偉琦,1989; 任鎮寰等,1996)。對于缺乏儀器數據的歷史地震,進行強地面運動模擬是評估該地區潛在地震危害性的重要手段。作為廣西乃至華南內陸有記錄以來最大的歷史地震,開展靈山M6地震的強地面運動模擬,進而補充本地區強震記錄數據庫,具有重要的科學價值和工程意義。但是由于歷史久遠和儀器記錄匱乏等原因,自發震至21世紀初70余年的時間里,靈山地震的研究程度一直較低(李細光等,2017b)。近年來,隨著技術手段進步,相關學者對靈山地震開展了一系列研究,獲得了一批突破性成果,已基本厘清此次地震的震級、震中位置、極震區等級范圍和發震構造等關鍵信息(周本剛等,2008; 何軍等,2012; 李細光等,2017a、2017b; 李冰溯等,2018; 李細光等,2018)。因此,針對靈山地震開展強地面運動模擬的條件已經成熟。

本研究使用隨機振動有限斷層法對1936年廣西靈山M6地震進行了強地面運動模擬,在收集地震地質資料的基礎上,確定發震斷層的走向、傾角、長度、寬度、震級、應力降等參數,選取適當的傳播路徑參數和場地效應參數,將上述參數輸入隨機有限斷層模型進行強地面運動模擬,并將模擬結果與地震烈度野外調查數據和地震動衰減關系進行對比,驗證模擬結果的可靠性。本研究最終獲得1936年廣西靈山M6地震區域及極震區峰值加速度(PGA)、峰值速度(PGV)及烈度的空間分布,相關結果可為工程抗震設計、政府決策和地震災害風險評估提供依據。

1 方法與原理

1.1 隨機振動有限斷層方法

隨機振動法是由Hanks等(1981)和Boore(1983)于20世紀80年代提出的一種用隨機振動來模擬高頻地震波的方法。這種方法的本質是把遠場地震動看作有限帶寬、有限持時的高斯白噪聲,并把影響地震動參數的一系列因素(包括震源、路徑、場地)在頻率域里表達為簡單的乘積形式,然后通過窗函數調整噪聲的傅里葉譜,最后通過傅里葉反變換得到地震動參數的時程曲線(Boore,1983、2003)。對于一次大地震,合成地震動的有效方法是基于小震模擬。在隨機有限斷層模型中,斷層破裂面被分成N個大小相等的矩形子源,每個子源即為一個點源,每個子源滿足ω平方衰減規律,通過設定震源模式和破裂傳播速度,可以得到子源破裂的時間順序,然后根據子源與場地的幾何關系,計算每個子源對場地的影響。所有子源在觀測點引起的地震動在時域中以適當的延遲時間疊加,便可獲得場地的地震動時程a(t),即

(1)

其中,NL、NW分別為斷層沿長度和寬度劃分的子斷層數,i和j分別代表沿斷層走向和傾向方向的斷層指標,Δtij為子源到場地的滯后時間,aij(t)為子源利用隨機點源方法得到的在觀測點的地震動時程。針對有限斷層模型采用部分子斷層多次破裂的方式來保證地震矩守恒的不合理之處,Motazedian等(2005)提出了“動力學拐角頻率”的概念,認為子斷層拐角頻率隨著破裂的傳播而減小。每個子斷層動力學拐角頻率表達式為

fcij(t)=4.9×106β(Δσ/M0ave)1/3×NR(t)-1/3×S

(2)

其中,fcij(t)為第ij個子斷層的動力學拐角頻率,t為第ij個子源被觸發的時刻,NR(t)為在時刻t已破裂子斷層的累計數,M0ave=M0/N為子斷層的平均地震矩,S為表示子斷層輻射強度的一個常數。為了使斷層總輻射能守恒,引入一個標度因子H。動力學拐角頻率對子斷層總輻射能的遞減影響得到標度因子H的補償,總輻射能保持不變,這樣在保證輻射能守恒的前提下避免了子斷層多次破裂這一非物理現象。至此,隨機振動有限斷層方法的理論框架得到最終完善。

1.2 場地條件對地震動模擬結果的影響

(3)

圖1 研究區內地形(a)及 分布(b)

2 模型參數

2.1 發震構造

發震斷層的幾何形態(包括走向、傾角、長度和寬度等)是隨機振動有限斷層模型的基本輸入參數。靈山M6地震所在的華南地區由于歷史強震較少,人類改造活動頻繁且降水較多,地震遺址保存非常困難,因此很長一段時間內關于此次地震的發震斷層存在爭議。陳國達(1939)最早研究靈山M6地震時,認為靈山地震發震斷層為NE走向的靈山—大容山斷裂,而茹錦文等(1985)則認為發震斷層應歸屬于NW向的昆侖關斷裂南延部分。前人對該地震發震斷層的研究大多從極震區震害和地裂縫的空間分布情況進行推測(陳國達,1939; 茹錦文等,1985; 李偉琦,1989; 任鎮寰等,1996),缺乏直接的地質學證據。2015年研究人員首次在靈山斷裂北段發現靈山地震2條地表破裂帶遺跡(李細光等,2017a、2017b),地表破裂帶沿線開挖的探槽清晰揭露出斷層對上覆土層的錯動,這為靈山斷裂北段晚更新世以來的活動提供了直接證據。探槽顯示的最新一次地震事件距今約80年,與靈山地震發震時間吻合,證明了靈山斷裂北段即為1936年靈山M6地震的發震斷裂。另外,李偉琦(1989)和任鎮寰等(1996)根據極震區的等震線形狀和低烈度區的長軸方向,推測靈山地震可能為ENE向斷裂與NNW向斷裂共軛破裂的結果,其中ENE向構造為控震構造。李冰溯等(2018)綜合極震區震害及地質學研究成果,認為靈山地震主發震構造為NE-ENE走向的靈山斷裂,羅陽山南麓的泗州斷裂同時發生共軛破裂,造成局部烈度增強。因此,本文在開展模擬時將發震斷層設定為2條,即靈山斷裂北段和泗州斷裂(圖2),2條斷裂的幾何參數見表1。

表1 模型計算輸入參數

注: 斷層數據來自國家地震科學數據中心活斷層數據分中心[注] https://www.activefault-datacenter.cn/。

2.2 斷層滑動模型

斷層滑動模型即斷層面上的滑動量靜態分布是進行強地面運動模擬的必要參數,一般通過震后地表形變數據或地震波數據反演獲得(申文豪等,2019)。對于靈山M6地震這類沒有觀測數據的歷史地震,已很難通過反演方式獲得滑動模型。除了反演模型之外,地震學者還根據地震波觀測總結出斷層滑動模型的空間分布規律,提出了通過正演方式生成斷層滑動分布的方法(王海云,2004)。本文采用應用廣泛的k平方模型正演獲得斷層面滑動分布。Gallovi等(2004)提出的k平方模型二維波數譜為

(4)

(5)

圖3 正演獲得的靈山斷裂(a)和泗州斷裂(b)斷層面滑動分布模型

2.3 震源及路徑傳播參數

根據現有文獻及地質資料基本可以確定1936年靈山M6地震發震斷層的破裂長度、走向、傾角、埋深、震中位置等幾何參數(表1)。除此之外,為描述斷層破裂和傳播特征,還需要確定地震矩、應力降、kappa值、品質因子Q模型、幾何衰減模型等參數。龍政強等(2015)選取廣西數字地震臺網記錄到的2007—2013年廣西主要斷裂帶ML≥2.5地震事件,采用Atkinson方法反演得到介質品質因子Q=366.3f0.47。根據Tsang等(2010),Kappa值κ可以由經驗關系式確定,即

κ=0.145-0.12ln(Vuc)

(6)

其中,Vuc為上地殼4km深度內平均橫波速度。根據廣西中西部地區背景噪聲成像結果,靈山地區Vuc約為3.25km/s(黃強強等,2022),因此Kappa值計算結果為0.0036。楊曉瑜等(2021)利用中國國家地震臺網336個固定地震臺站記錄的遠震波形資料,通過P波接收函數分析估算了中國華南地區的地殼厚度,其結果顯示靈山地震震源區地殼厚度約為28km,根據三段式幾何衰減模型的定義(Atkinson等,1992),該地區幾何衰減模型可以表示為

(7)

其中,R為震中距。龍政強等(2015)研究顯示,廣西地區10MPa及其以上的高應力降地震主要分布在桂東南,且歷史上發生過5.0級以上地震的震源區及其附近仍具有高應力背景。2010年8月17日發生的靈山ML3.6地震其應力降為11.54MPa,該地震位于1936年靈山M6地震和1958年9月25日靈山M5地震震中附近,這說明靈山地震區具有高應力背景,本研究將應力降的值取為8.0MPa。

3 結果分析與討論

靈山M6地震雖然缺乏儀器記錄,但具有豐富的烈度調查數據。本研究針對宏觀區域(調查烈度覆蓋區域及其擴展區域)和極震區(Ⅷ度和Ⅸ度區)分別劃分不同的計算區域,并設置不同尺度的計算網格,利用隨機振動有限斷層模型計算每個網格點的PGA和PGV,通過經驗轉換關系(丁寶榮等,2017)分別轉化為烈度,最后取兩種烈度的加權平均為最終烈度值,并將其與野外調查結果進行比較。丁寶榮等(2017)給出的PGA、PGV與烈度MMI轉換關系為

(8)

3.1 宏觀烈度模擬

宏觀烈度模擬區域劃定在106.0°E～112.0°E、20.5°N～25.5°N范圍內,區域網格大小為6.0°×5.0°,經度及緯度方向網格精度分別為0.06°和0.05°。為更好地展示強地面運動參數衰減趨勢,繪制了lg(PGA)和lg(PGV)空間分布圖,如圖4(a)、4(b)所示。從圖中可以看到模擬的強地震動參數具有以下特征:①PGA最大值在900gal左右,PGV最大值達100cm/s,轉化為烈度均已超過Ⅸ度; ②PGA和PGV均隨距斷層距離增加而迅速衰減,且PGA的衰減趨勢更快; ③PGA和PGV均受到場地條件的影響,在平原等沉積層較厚區域存在不同程度的放大效應。

圖4 基于隨機振動模型模擬得到的靈山M6地震宏觀區域強地震動參數分布

3.2 極震區模擬

為更好地研究雙斷層共軛破裂造成的強地面運動特點,本研究針對極震區開展了網格尺度更精細的數值模擬,模擬區域劃為109.1°E～109.7°E、22.2°N～22.8°N,區域網格大小為0.6°×0.6°,經度及緯度方向網格精度為0.006°。極震區地震動模擬結果如圖5所示,圖中白色實線為野外調查Ⅷ度和Ⅸ度區等值線,黑色實線為模擬烈度Ⅷ度和Ⅸ度區等值線。國家地震局全國地震烈度區劃編圖組(1979)匯編的《中國地震等烈度線圖集》中,靈山M6地震的震中烈度最大為Ⅸ度,并且形狀為一長軸方向呈NE向的橢圓形?！吨袊卣鹉夸洝?顧功敘,1983)和《中國近代地震目錄》(中國地震局震害防御司,1999)均采用這一結果。此后,陳恩民等(1984)、李偉琦(1989)以及任鎮寰等(1996)對極震區烈度進行了重新調查和評定,結果均顯示極震區的烈度分布有NEE和NWW兩個優勢方向,呈“T”型分布。由圖5可以看到,無論是PGA、PGV分布,還是烈度分布,本研究的模擬結果較好地反映了極震區烈度分布的這種“T”型特征,尤其是在Ⅸ度區。這一結果從模擬的角度驗證了模型參數的合理性,也證實了靈山地震極震區分布是由靈山斷裂及泗州斷裂共軛破裂造成的。兩斷層發生共軛破裂的控制因素可能是羅陽山地區最大有效力矩準則約束下的以NWW-SEE向為主壓應力方向的力學系統(張沛全等,2012)。圖5(c)還顯示出與調查烈度不同的是,模擬烈度Ⅷ度區在北部及西南部有兩處明顯的烈度異常區,而這兩處烈度異常區均為河流沖積平原,北部為郁江平原區,西南部為欽江平原區。這兩處平原區沉積層較厚,水平較低,從而導致明顯的場地放大效應。

3.3 衰減關系對比

除了與調查烈度進行空間分布對比之外,本研究還將宏觀區域和極震區模擬得到的地震動參數與衰減關系NGA(Boore等,2007)進行對比。圖6(a)、6(b)分別展示了宏觀區域模擬PGA、PGV與NGA的衰減趨勢,圖6(c)、6(d)分別展示了極震區內模擬PGA、PGV與NGA的衰減趨勢,圖中黑色實線為NGA衰減關系曲線,黑色虛線為95%置信區間曲線。從圖6 可以看到,模擬結果與NGA結果的衰減趨勢大體一致,在極震區內模擬結果與NGA曲線符合程度較高,但在震中距大于100km時,地震動模擬結果整體略高于NGA預測值。圖7 展示了本研究模擬的烈度與王繼等(2008)針對本地區烈度衰減公式的對比,其中藍色實線為長軸衰減公式,藍色虛線為短軸衰減公式。從圖7可以看到,模擬烈度與理論公式在不同震中距尺度上具有相似的衰減特征。在震中距小于5km的范圍內,基于烈度衰減公式預測的最大烈度值為Ⅷ度強,顯然低于模擬烈度和實際調查烈度。造成這一差異原因可能是王繼等(2008)所使用的地震事件中超過6.0級的地震樣本數較少,從而可能會低估大震級事件極震區內的烈度水平。在震中距大于200km處,模擬烈度則低于基于烈度衰減公式預測的烈度水平?？紤]到在遠場區域地震動主要以面波形式出現,而品質因子Q對其有著重要影響,在本研究過程中只采用了單一的Q衰減模型。事實上,許多研究已經表明Q值是依賴于震中距變化的(Tsang等,2010),因此,在今后工作中需要建立隨距離變化的Q值模型來改善模擬結果。

注: (a)和(b)分別為宏觀區域模擬PGA和PGV與NGA對比; (c)和(d)分別為極震區內模擬PGA和PGV與NGA對比。

圖7 模擬得到的烈度值與華南地區烈度衰減關系(王繼等,2008)對比

4 結論與展望

本研究收集整理了1936年靈山M6地震的地質資料、活動斷層探測最新成果等,利用隨機振動有限斷層模型計算了宏觀區域及極震區內網格點的峰值加速度及峰值速度等參數,并且加入了淺層橫波速度結構對模擬結果的影響,最終得到此次地震的地震動分布,分析了地震動特征。本研究將模擬結果與歷史地震烈度調查數據及地震動衰減關系進行對比,結果顯示模擬結果與調查烈度值在整體特征、極震區的分布等方面均符合較好。

由于廣西地區中強地震少發,缺乏足夠的強地面運動記錄回歸擬合得到本地區地震動衰減關系,而模擬得到的強震數據可以補充強震數據的不足,為建立廣西地區地震動衰減關系曲線提供有益補充。另外,由于隨機振動有限斷層模型未考慮地震波在介質中的傳播過程,不需要像復合震源模型、有限差分法或者譜元法等方法耗費大量時間用于地質體網格剖分和格林函數計算,具備時效性強的突出特點,這一點正是地震烈度速報所需要的。因此,未來隨機振動有限斷層模型可以應用于本地區地震烈度速報,為震后災情快速評估和應急救援提供有效幫助。