張獻(xiàn)兵, 張偉, 陳曉非*, 蔡永恩
1 北京大學(xué)地球與空間科學(xué)學(xué)院地球物理系, 北京 1008712 南方科技大學(xué)地球與空間科學(xué)系, 廣東深圳 518055
北京地區(qū)位于華北平原地震帶中,這個(gè)地區(qū)及其鄰區(qū)歷史上有記錄的6級以上強(qiáng)震曾發(fā)生多次,例如,公元294年延慶東6級地震、1057年大興南6.8級地震、1484年延慶6.8級地震、1536年通縣南6級地震、1665年通縣6.5和1730年京西6.5級地震(謝毓壽和蔡美彪,1983;顧功敘,1984).在第五代中國地震動區(qū)劃圖中,北京地區(qū)50年超越概率10%的地震動峰值加速度大部分被設(shè)定為0.2g,對應(yīng)的地震烈度為Ⅷ度(中華人民共和國國家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局和中國國家標(biāo)準(zhǔn)化管理委員會,2015).
焦青等(2005)根據(jù)八寶山—黃莊—高麗營斷裂第四紀(jì)地質(zhì)構(gòu)造活動特征和首都圈地區(qū)跨斷層位移觀測資料,得到的北京地區(qū)區(qū)域應(yīng)力狀態(tài)和斷層活動特性,推測首都圈及其鄰區(qū)即將進(jìn)入新的地震活躍期.武敏捷等(2013)綜合利用b值、地區(qū)歷史和現(xiàn)今地震活動、斷裂帶分段地震活動參數(shù)等信息,研究了首都圈現(xiàn)今地震活動特征,提出未來延懷盆地具有發(fā)生中強(qiáng)地震的危險(xiǎn).秦向輝等(2014)等利用水壓致裂法地應(yīng)力測量和地應(yīng)力變化監(jiān)測臺獲得的北京地區(qū)地殼淺表層地應(yīng)力,結(jié)果表明如果摩擦系數(shù)弱化到0.4,黃莊—髙麗營斷裂附近測點(diǎn)應(yīng)力狀態(tài)將趨向失穩(wěn)滑動.張群偉和朱守彪(2019)計(jì)算了1820年以來華北主要地震的庫侖應(yīng)力變化,認(rèn)為黃莊—高麗營斷裂帶的值變大,發(fā)生地震的風(fēng)險(xiǎn)相對提高.近年來,人們通過綜合地球物理與地質(zhì)方法開展了北京地區(qū)活動斷裂研究,可控源音頻大地電磁測深、淺層人工地震、高密度電阻率法、鉆孔剖面方法均顯示黃莊—高麗營斷裂、南口—孫河斷裂和順義斷裂等具有第四紀(jì)分段活動特征(徐錫偉等,2015;張磊等,2017a,b;周永恒等,2021).最近的歷史地震和地震目錄研究表明,京津冀地區(qū)在未來一段時(shí)間仍處于地震活躍期,存在發(fā)生強(qiáng)震的風(fēng)險(xiǎn)(胡慧文等,2021).
上述研究表明北京地區(qū)未來存在發(fā)生破壞性地震的可能性,迫切需要對這個(gè)地區(qū)強(qiáng)震一旦發(fā)生可能引起的強(qiáng)地面運(yùn)動進(jìn)行研究,為抗震設(shè)防提供參考.高孟潭等(2002)用有限差分法模擬了1679年三河—平谷8級地震斷層破裂對北京地區(qū)的強(qiáng)地面運(yùn)動的影響.他們發(fā)現(xiàn)北京盆地的放大效應(yīng)與北京市區(qū)地表土層的放大效應(yīng)具有同等水平,并且提出為了全面考慮盆地效應(yīng),需要對更多的設(shè)定地震進(jìn)行研究.Ding等(2004)用橫向不均勻非彈性介質(zhì)中合成地震圖的計(jì)算方法,模擬了北京強(qiáng)地面運(yùn)動,發(fā)現(xiàn)厚的第三紀(jì)和第四紀(jì)沉積地層對強(qiáng)地面運(yùn)動的峰值有放大作用.潘波等(2006)利用Graves的三維有限差分程序,計(jì)算了1679年三河—平谷8級地震在北京地區(qū)引起的地面峰值速度,認(rèn)為有必要開展對各條發(fā)震斷層的詳細(xì)研究,以獲得某一斷層發(fā)震后強(qiáng)地面運(yùn)動值的分布.劉博研等(2007)模擬三河—平谷地震北京地區(qū)強(qiáng)地面運(yùn)動,給出了寬頻帶的三分量復(fù)合地震記錄,得到北京地區(qū)強(qiáng)地面地表水平運(yùn)動峰值加速度在2 m·s-2左右.周紅和陳曉非(2008)利用新局域離散波數(shù)法模擬了1998年張北地震引起的北京西部凹陷區(qū)二維強(qiáng)地面運(yùn)動,結(jié)果表明薄的第四紀(jì)沉積蓋層放大效應(yīng)遠(yuǎn)大于厚的第三紀(jì)蓋層.楊勇等(2008)基于地震活動性模型,采用地震危險(xiǎn)性分析方法分別計(jì)算華北地區(qū)地震動峰值加速度和各大城市的危險(xiǎn)性曲線,結(jié)果表明相比華北區(qū)域內(nèi)其他城市,唐山、太原和北京等地區(qū)的潛在地震危險(xiǎn)性更高.付長華等(2015)利用有限差分?jǐn)?shù)值算法,設(shè)定震源為三河—平谷地震、頤和園地震和通縣地震,研究了北京盆地對3~10 s周期地震動加速度反應(yīng)譜的放大系數(shù),探討了高層建筑地震致災(zāi)風(fēng)險(xiǎn)與盆地結(jié)構(gòu)的關(guān)系.廖力(2017)利用譜元法,地形數(shù)據(jù)采用etop1′,分析對比首都圈4.7級地震的地形效應(yīng),結(jié)果表明三維模型模擬的信號在0.02~0.5 Hz范圍內(nèi)接近真實(shí)臺站記錄,由于地形的影響,峰值加速度(PGA)分布模式較為復(fù)雜.巴振寧等(2022)采用譜元法模擬三河—平谷地震北京地區(qū)強(qiáng)地面運(yùn)動,雙側(cè)破裂方式的結(jié)果表明北京中心城區(qū)峰值加速度和速度分別達(dá)到0.12 g和1.0 m·s-1.
以上研究主要關(guān)注北京地區(qū)以外強(qiáng)震對北京地區(qū)地面運(yùn)動的影響,著重震源機(jī)制、場地效應(yīng)等對強(qiáng)地面運(yùn)動的影響.對于此區(qū)主要斷裂發(fā)生的幾次歷史強(qiáng)震,宏觀地震烈度的估計(jì)主要是來自野外考察和歷史地震的研究.基于北京已有的活動斷裂研究,本文將首先探討北京地區(qū)主要潛在活動斷裂的地震安全度,然后選擇其中的黃莊—高麗營斷裂(北段),設(shè)定一個(gè)與斷層尺度有關(guān)的地震大小,探討這個(gè)地震引起的強(qiáng)地面運(yùn)動特征和強(qiáng)度分布細(xì)節(jié),為防震減災(zāi)提供參考依據(jù).
斷裂是活動構(gòu)造的一種類型.活動構(gòu)造是指晚更新世(距今10~12萬年)以來一直在活動,未來一定時(shí)期內(nèi)仍可能活動的斷裂(鄧起東,1991).北京地區(qū)位于沉積平原,主要斷裂大多數(shù)都是隱伏的正斷層.由于不同作者使用的研究方法和關(guān)注點(diǎn)不同,文獻(xiàn)上給出的這個(gè)地區(qū)的斷裂帶長度、產(chǎn)狀和數(shù)目,均存在差異.本文以圖1給出的斷裂帶作為研究地震危險(xiǎn)性的主要斷裂帶.北京地區(qū)活動斷裂主要有:(1)南口—孫河斷裂,長度約50~80 km,走向北西,斷裂北西段和南東段分別傾向南西和北東,傾角70°~80°(車兆宏,1994;梁亞南,2019); (2)黃莊—高麗營斷裂,長約110~132 km(南口—孫河斷裂以北~40 km),走向北東—北東東,傾向南東—南東東,傾角50°~80°(張磊等,2017a;周永恒等,2021);(3)順義—良鄉(xiāng)斷裂,走向北東25°~30°,傾向南東,全長~100 km.以南口—孫河斷裂為界,北段稱為順義—孫河斷裂(35~40 km),走向北東25°~40°,傾向南東,傾角60°~80°;南段稱為孫河—良鄉(xiāng)斷裂(~55 km),走向北東,傾向北西,傾角60°(張世民等,2005;焦青和邱澤華,2006);(4)南苑—通縣斷裂長約110 km,北段長約20 km,走向北東,傾向北西,南段長90 km,走向北東,傾向北西,斷裂帶總體傾角為50°~75°(焦青和邱澤華,2006).
圖1 研究區(qū)域內(nèi)北京地區(qū)主要活動斷裂
研究區(qū)域如圖1所示,在歷史上曾經(jīng)發(fā)生多次6級以上地震.在南口—孫河斷裂帶,據(jù)探槽古地震研究,大約4000~10000年前曾經(jīng)發(fā)生過3次7級以上古地震(江娃利等,2001),在黃莊—高麗營斷裂附近的清河隱伏活動斷層上,曾發(fā)生1730年圓明園6.5級地震(116.6°E,39.96°N)(環(huán)文林等,1996),南苑—通縣斷裂帶曾發(fā)生1057年大興MS6.75地震和1665年通州MS6.5地震.自1730年以來北京地區(qū)沒有發(fā)生過大于等于6級的地震(江宜和殷娜,2021).
構(gòu)造地震的發(fā)生取決于斷裂強(qiáng)度及其所處的應(yīng)力狀態(tài).應(yīng)力狀態(tài)主要與斷層的幾何形態(tài)、巖石力學(xué)性質(zhì)、構(gòu)造應(yīng)力場等有關(guān).目前可以通過水壓致裂和鉆孔應(yīng)變測量方法,獲得地表淺層水平方向的應(yīng)力.對于震源深處的應(yīng)力狀態(tài)還不能直接測量,只能通過一些間接方法進(jìn)行估計(jì).不考慮地下流體、溫度和化學(xué)等因素對斷層應(yīng)力場的影響,對于給定的斷層幾何形態(tài)和區(qū)域構(gòu)造應(yīng)力場,王仁等(1981)假設(shè)斷層破裂遵從理想塑性材料的庫侖—摩爾破裂準(zhǔn)則,提出用有限元方法探討斷層的構(gòu)造應(yīng)力狀態(tài),并利用這個(gè)應(yīng)力狀態(tài)和斷層的摩擦系數(shù),定義了一個(gè)斷層能否發(fā)生地震的“地震安全度”:
(1)
其中τn、σn和μ分別為斷層面上的剪應(yīng)力、正應(yīng)力(拉張為正)和內(nèi)摩擦系數(shù),μσn稱為斷層摩擦強(qiáng)度.從式(1)可見,安全度G值越大,斷層越安全,地震的危險(xiǎn)性越小,安全度G等于零,即斷層的剪應(yīng)力達(dá)到了破裂極限,意味著發(fā)生地震.
本文利用王仁等提出的地震安全度的概念和有限元方法,計(jì)算北京地區(qū)主要地震斷裂的地震安全度.根據(jù)已有的北京活斷層研究和地震活動性,選取南口—孫河斷裂(西段和中段)、黃莊—高麗營斷裂(北段和中段)、良鄉(xiāng)—順義斷裂(北段)、南苑—通縣斷裂(北段)和南口—山前斷裂為可能發(fā)生強(qiáng)震的斷層(圖1).由于這些斷裂都是高傾角的隱伏斷裂,現(xiàn)今具有走滑性質(zhì),本文把實(shí)際問題簡化為平面應(yīng)變問題,其邊值問題可以歸結(jié)為
σyy|y=W=σ1sin2θ+σ2cos2θ,
τyx|y=W=(σ2-σ1)cosθsinθ,
ux|x=0=0,
τxy|x=0=0,
uy|y=0=0,
τyx|y=0=0,
(2′)
其中x和y分別為模型的空間坐標(biāo),它們的正方向分別指向東和北,模型長為L,寬為W,σxx、σyy和τxy分別為模型內(nèi)一點(diǎn)處x、y方向的正應(yīng)力和剪應(yīng)力,εxx、εyy和γxy分別為正應(yīng)變和剪應(yīng)變,σ1和σ2分別為區(qū)域構(gòu)造應(yīng)力場的最大和最小主應(yīng)力,θ為最大主應(yīng)力方向與x軸的夾角,ux和uy分別為x和y方向的位移分量,E和ν為分別楊氏模量和泊松比.區(qū)域構(gòu)造應(yīng)力場的最大和最小主應(yīng)力分別為-200 MPa和-600 MPa,最小主應(yīng)力方向?yàn)楸睎|70°(王仁等,1981).為了方便比較,斷裂帶內(nèi)的摩擦系數(shù)均取為0.6.根據(jù)模型區(qū)域內(nèi)地震波平均速度(縱波為4214 m·s-1、橫波為2200 m·s-1)和平均密度(2375 kg·m-3),得到楊氏模量和泊松比分別為30.2 GPa和0.32,斷層內(nèi)的楊氏模量比斷層外小一個(gè)數(shù)量級.
有限元計(jì)算結(jié)果如圖2所示.在上述區(qū)域應(yīng)力構(gòu)造環(huán)境下,黃莊—高麗營斷裂(北段)的地震安全度相對最低,平均為0.25,其次為南口—孫河斷裂(西段)、順義—良鄉(xiāng)斷裂(北段)、南口—孫河斷裂(中段)、黃莊—高麗營斷裂(中段)和南口—山前斷裂,南苑—通縣斷裂(北段)的安全度相對最高.
圖2 研究區(qū)域內(nèi)北京地區(qū)主要斷裂帶地震安全度G
本文選擇黃莊—高麗營斷裂(北段)設(shè)定地震,模擬其破裂引起的北京地區(qū)的強(qiáng)地面運(yùn)動.設(shè)定地震震級根據(jù)龍鋒等(2006)的發(fā)震斷層長度L與地震面波震級MS的統(tǒng)計(jì)關(guān)系:MS=3.821+1.860 lgL(剩余標(biāo)準(zhǔn)差為0.317),取黃莊—高麗營斷裂(北段)的長度為45 km,得到最大可能面波震級為7.2.
本文數(shù)值模擬區(qū)域?yàn)闁|經(jīng)115.99°至116.80°,北緯39.66°至40.34°,面積大約為90 km×75.6 km(圖1).
本文綜合中國地震局地球物理勘探中心提供的速度模型,程捷提供的第四系土壤層速度模型,以及北京地區(qū)的地震淺層勘探結(jié)果、地質(zhì)資料、鉆孔資料、地層波速試驗(yàn)資料(上述資料均由《北京市活斷層探測與地震危險(xiǎn)性評價(jià)》項(xiàng)目提供)構(gòu)建三維速度結(jié)構(gòu)模型(參數(shù)見表1).各層底界面的鳥瞰圖如圖3所示.圖3c為研究區(qū)內(nèi)第四系底界面,第四系沉積層厚度分布并不均勻,西部較淺,東部較深,西北部昌平凹陷和東北部順義凹陷深達(dá)幾百米,其中東北部地區(qū)部分最深點(diǎn)達(dá)760 m左右,位于順義天竺.
表1 北京地區(qū)7層三維介質(zhì)模型物理參數(shù)
圖3 速度結(jié)構(gòu)模型各層底界面鳥瞰圖(a)Hq3+4, (b)Hq2, (c)Hq1, (d)Hn, (e)Hg, (f)Hc, (g)Hm.
設(shè)定地震的黃莊—高麗營斷裂在北京地震臺網(wǎng)建立之后沒有可用于反演這些斷層詳細(xì)破裂過程的地震記錄.對于未知破裂過程的斷層一旦破裂引起的地震危害性分析,存在三種可能的處理方案:一是已知初應(yīng)力分布的情況下,可以采用斷層破裂的動力學(xué)模擬方法獲得破裂過程(Zhang et al.,2014,2017);二是采用“借用”的方案,使用其他地區(qū)斷層的破裂反演結(jié)果,作為設(shè)定地震的可能破裂過程;三是采用Somerville等(1999)提出的特征斷層(characterized source model)作為預(yù)設(shè)斷層破裂過程.本文采用特征斷層方案作為預(yù)設(shè)斷層破裂過程.特征斷層基于斷層的凹凸體(Asperity)模型, 是對斷層不均勻破裂過程中主要特征的概括和近似:通過不均勻分布在斷層上的凹凸體反映斷層破裂的不均勻性;凹凸體上的滑動量大于斷層其他區(qū)域的滑動量;在單個(gè)凹凸體內(nèi)部,通常假定具有相同的滑動量.
Somerville等(1999)通過對大量地震的統(tǒng)計(jì)分析,獲得了凹凸體模型參數(shù)之間的經(jīng)驗(yàn)規(guī)律,參數(shù)包括:斷層破裂區(qū)域的面積S、地震矩M0、破裂面上的平均應(yīng)力降Δσc、斷層面上凹凸體區(qū)域的面積Sa、凹凸體應(yīng)力降Δσa、凹凸體的數(shù)目N及其在斷層上的位置、凹凸體平均位錯(cuò)量Da、凹凸體外區(qū)域及其應(yīng)力σb、斷層滑動速度、震源時(shí)間函數(shù)以及起始破裂點(diǎn)的位置等.建立特征斷層凹凸體破裂模型時(shí),本文主要參考以下經(jīng)驗(yàn)規(guī)律:
(a) 震級同地震矩的經(jīng)驗(yàn)公式(沙海軍和呂悅軍,2018;Hanks and Kanamori,1979)
MW=2/3 lgM0-10.7,
(3)
(b) Asperity的總面積同斷層面積的比值
Sa/S=0.22,
(4)
(c) Asperity的數(shù)目N同活動斷層分段成正比,并且不同的N時(shí),Asperity上的平均位錯(cuò)Da同總平均位錯(cuò)D具有關(guān)系
N=1∶Da/D=2.3,
N=2∶Da/D=2.0,
N=3∶Da/D=1.8,
(5)
Asperity的位置由斷層的幾何形態(tài)確定.
(d) 采用了Somerville(1993)中的上升時(shí)間同地震矩的經(jīng)驗(yàn)關(guān)系數(shù)值模擬使用的斷層和震源具體參數(shù)如圖4和表2所示.
表2 黃莊—高麗營斷裂(北段)設(shè)定地震震源參數(shù)
圖4 黃莊—高麗營斷裂(北段)設(shè)定地震震源模型(a) 斷裂地面投影; (b) 斷裂空間位置; (c) 斷裂面凹凸體(黃色塊)位置和破裂時(shí)間.
(6)
地震斷層破裂引起的強(qiáng)地面運(yùn)動的計(jì)算,數(shù)學(xué)上歸結(jié)為求解一個(gè)半無限空間內(nèi)一個(gè)移動破裂面源的波動方程問題,由于介質(zhì)的不均勻性、斷層幾何的復(fù)雜性和地表的起伏,一般只能采用數(shù)值方法求解.本文利用張偉等(張偉,2006;Zhang and Chen,2006;祝賀君等,2009;Zhang et al.,2012;Sun et al.,2017)開發(fā)的能夠模擬復(fù)雜地形波場的曲線網(wǎng)格有限差分方法模擬黃莊—高麗營斷裂(北段)設(shè)定地震的近場長周期強(qiáng)地面運(yùn)動.
數(shù)值模型的網(wǎng)格數(shù)目為900×756×214;空間步長均為100 m,當(dāng)深度大于10 km后深度步長變?yōu)?00 m;時(shí)間步長為0.005 s,時(shí)間步數(shù)為10000,總時(shí)間為50 s.
黃莊—高麗營斷裂(北段)MS7.2設(shè)定地震在北京地區(qū)地表不同時(shí)刻的速度波場(南北向分量)快照如圖5所示.從圖中可以看到,在8 s前地表速度峰值不大,之后凹凸體破裂導(dǎo)致的強(qiáng)振動開始在地表顯著;9 s時(shí)在斷層?xùn)|北側(cè)出現(xiàn)較大面積向北運(yùn)動;由于地下結(jié)構(gòu)的復(fù)雜性,地震動的傳播表現(xiàn)出較為復(fù)雜的特征,如16 s的波場快照顯示波場在較厚的沉積內(nèi)出現(xiàn)較強(qiáng)的反復(fù)震蕩,持續(xù)時(shí)間較長.
圖5 黃莊—高麗營(北段)設(shè)定地震的速度波場傳播快照圖
圖6—8分別是數(shù)值模擬得到的地表速度、位移和加速度峰值分布.由圖可見,峰值速度(PGV)大于0.3 m·s-1的區(qū)域主要集中在斷裂帶附近;最大峰值速度出現(xiàn)在斷層的中部偏東側(cè),其值可達(dá)1 m·s-1;PGV>0.7 m·s-1的區(qū)域集中在沿?cái)嗔褨|側(cè)的8~9 km寬的條帶內(nèi);除天安門附近和第四紀(jì)地層較厚的地區(qū)(例如昌平凹陷)外,不在斷層附近地區(qū)的峰值速度基本小于0.3 m·s-1.峰值位移(PGD)分布不同于PGV,在沿?cái)鄬觾蓚?cè)分布較寬;最大峰值主要集中在斷層的東北端,其值達(dá)0.9 m.
圖6 黃莊—高麗營斷裂(北段)設(shè)定地震的地面運(yùn)動峰值速度
圖7 黃莊—高麗營斷裂(北段)設(shè)定地震的地面運(yùn)動峰值位移
圖8 黃莊—高麗營斷裂(北段)設(shè)定地震的地面運(yùn)動峰值加速度
由于本文模擬計(jì)算主要為長周期(低頻)運(yùn)動,PGA不適合直接用于地震動強(qiáng)度評估,因此采用PGV來計(jì)算地震烈度.計(jì)算公式使用的是《中國地震烈度表》(國家市場監(jiān)督管理總局和國家標(biāo)準(zhǔn)化管理委員會,2020)式(7):
IV=3.00 lg(PGV)+9.77,
(7)
其中PGV為峰值速度.模擬計(jì)算得到的地震烈度分布見圖9.
圖9 黃莊—高麗營斷裂帶(北段)設(shè)定地震的烈度分布圖
強(qiáng)地面運(yùn)動數(shù)值模擬結(jié)果表明,對于設(shè)定地震為MS7.2的黃莊—高麗營斷層(北段)如果發(fā)生地震,斷層附近20 km內(nèi)的地區(qū)烈度大都將在Ⅷ度以上;在斷層附近北東方向烈度超過Ⅸ度;在斷層西南端,Ⅷ度以上地區(qū)沿?cái)鄬游髂戏较蜓由炝?0 km左右;遠(yuǎn)離斷層20 km的區(qū)域烈度大都在Ⅷ度以下;第四紀(jì)較厚的區(qū)域(如昌平凹陷)烈度比周邊較高.
本文采用的數(shù)值方法能夠模擬地震動的震源、傳播途徑、場地放大等因素對強(qiáng)地面運(yùn)動的影響,可以評價(jià)潛在地震事件對目標(biāo)區(qū)的地震動烈度.本文得到的結(jié)果依賴于使用的斷層幾何和震源模型以及速度結(jié)構(gòu),不同的模型會得到不同的結(jié)果.由于近場地震動的大小會受到斷層和場地效應(yīng)的影響,地震動不可能有統(tǒng)一的圖像.這就表明,只能對具體情況進(jìn)行具體分析,即使是同一個(gè)地區(qū),不同斷裂帶產(chǎn)生的強(qiáng)地面運(yùn)動也會不同.此外,地震動的峰值只是烈度區(qū)劃參數(shù)之一,應(yīng)結(jié)合其他參數(shù),不適合單獨(dú)使用.為了避免近場非規(guī)律性最大值的影響,本文的設(shè)定地震事件作為低超越概率地震,相當(dāng)或超出區(qū)劃圖中規(guī)定的大震或安評對象的重要工程中使用的最不利地震的概率水平.
致謝本文早期得到鄧起東院士、徐錫偉研究員、盧造勛研究員、楊主恩研究員的指導(dǎo)和大力支持.趙伯明教授對本文的工作給予了熱情指導(dǎo)與幫助.北京地震局領(lǐng)導(dǎo)對本研究給予了大力支持.兩位審稿專家和編輯給予本文寶貴的修改意見.在此一并表示衷心的感謝.