李春果，王宏偉，溫瑞智，強(qiáng)生銀，任葉飛

（1.中國地震局工程力學(xué)研究所地震工程與工程振動(dòng)重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150080；2.地震災(zāi)害防治應(yīng)急管理部重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150080）

引言

據(jù)日本氣象廳（JMA）測(cè)定，2022年3月16日23：36：29（UTC+9）日本本州東岸近海發(fā)生7.4級(jí)地震，震中位于141.623°E、37.697°N，震源深度為57 km，是2011年3月11日9.1級(jí)東日本大地震后日本東部地區(qū)最大震級(jí)地震。據(jù)日本防災(zāi)科學(xué)技術(shù)研究所（NIED）震度分布，地震最大震度達(dá)到6強(qiáng)，造成了強(qiáng)烈的地面運(yùn)動(dòng)、嚴(yán)重的工程結(jié)構(gòu)破壞及人員傷亡。此次地震的發(fā)震區(qū)位于歐亞板塊、太平洋板塊、北美板塊與菲律賓海板塊的交界，該區(qū)域板塊間相互作用復(fù)雜，太平洋板塊和菲律賓板塊向西俯沖形成狹長的板間地震帶，歐亞板塊沿邊緣東邊界向東移動(dòng)，形成雙向?qū)_式匯聚，擁有活躍的淺源和深源強(qiáng)震活動(dòng)［1－4］（見圖1）。地震發(fā)生后，密集的地震動(dòng)觀測(cè)臺(tái)網(wǎng)（K－NET和KiK－net）觀測(cè)到強(qiáng)烈的地震動(dòng)及豐富的長周期面波，其中K－NET臺(tái)網(wǎng)的MYG002臺(tái)站獲得了峰值地面加速度（PGA）最大的記錄，其東西、南北和垂直 向 的PGA分 別 為800.353、609.915、414.211 cm/s2。

圖1 日本浪江町7.4級(jí)地震震中及其3.0≤MJMA≤6.0余震（截至2022年3月21日）Fig.1 Epicenters of the M7.4 Namie,Japan earthquake and its aftershocks with 3.0≤MJMA≤6.0(up to March 21,2022)

日本東部海域地區(qū)強(qiáng)震活動(dòng)十分活躍，研究人員開展了大量的地震波動(dòng)數(shù)值模擬工作，F(xiàn)urumura等［5］、Takemura等［6］分別模擬了不同地震中日本關(guān)東盆地產(chǎn)生的長周期地震動(dòng)，前者分析了盆地內(nèi)地震波聚焦效應(yīng)引起的顯著的盆地放大，后者則重點(diǎn)討論了震源參數(shù)對(duì)長周期地震動(dòng)模擬的影響；Nakamura［7］等探討了表面地形、海水層對(duì)地震動(dòng)數(shù)值模擬的影響；Petukhin等［8］分析了包含海水層的三維速度結(jié)構(gòu)模型對(duì)大阪盆地模擬地震動(dòng)的影響，解釋了地震面波在海洋沉積層中的產(chǎn)生過程，強(qiáng)調(diào)了地震波動(dòng)模擬中考慮海水層的必要性。

文中利用日本7.4級(jí)地震的震源運(yùn)動(dòng)學(xué)破裂模型初步反演結(jié)果以及考慮地表地形、海水、板塊邊界等的三維速度結(jié)構(gòu)精細(xì)化模型，基于交錯(cuò)網(wǎng)格有限差分法，模擬大尺度區(qū)域內(nèi)地震波的產(chǎn)生和傳播過程，對(duì)比強(qiáng)震動(dòng)臺(tái)站的觀測(cè)記錄，對(duì)數(shù)值模擬結(jié)果進(jìn)行檢驗(yàn)，重點(diǎn)討論了復(fù)雜地形、海水層對(duì)地震動(dòng)傳播的影響。

1 數(shù)值模擬方法

有限差分法是一種求解微分方程數(shù)值解的近似方法，其主要原理是對(duì)微分方程中的微分項(xiàng)進(jìn)行直接差分近似，從而將微分方程轉(zhuǎn)化為代數(shù)方程組求解。交錯(cuò)網(wǎng)格有限差分法是運(yùn)用較為廣泛的方法之一［9－10］，該方法能夠處理包括起伏地表和海水平面在內(nèi)的復(fù)雜地形問題，采用局部差分算子增加運(yùn)算速度，同時(shí)減少數(shù)值頻散現(xiàn)象，提高計(jì)算精度［11－12］，結(jié)合高階有限差分能夠適當(dāng)加大網(wǎng)格間距，增加運(yùn)算效率［13－14］。文中采用交錯(cuò)網(wǎng)格有限差分法的開源程序OpenSWPC［15］（https：//tktmyd.github.io/OpenSWPC/）模擬地震波的產(chǎn)生和傳播過程，在空間上選擇四階精度，時(shí)間上為二階精度［10］。該開源程序求解了以速度－應(yīng)力分量表示的三維連續(xù)介質(zhì)力學(xué)運(yùn)動(dòng)方程，

式中：vi是第i個(gè)復(fù)合體中彈性運(yùn)動(dòng)的粒子速度：ρ是密度；σij是應(yīng)力張量的i，j個(gè)分量；fi是 體力的第i個(gè)分量，采用廣義齊納體（Generalized Zener Body）模型的本構(gòu)方程表示寬頻范圍的恒定Q值，將內(nèi)存變量納入本構(gòu)方程以表示粘彈性介質(zhì)，采用一階Crank-Nicolson法［16］顯式求解原始隱式方程。

采用有限差分法的OpenSWPC程序相較于有限元法、譜元法等程序求解形式簡(jiǎn)單，占用內(nèi)存小，計(jì)算效率高，且整合了模擬過程必要的前處理和后處理程序，并包含了日本全部地區(qū)精細(xì)化三維速度結(jié)構(gòu)模型。輸入相應(yīng)的地震震源或有限斷層參數(shù)，程序?qū)⒆詣?dòng)分配并行計(jì)算機(jī)內(nèi)存生成非均勻介質(zhì)模型的交錯(cuò)網(wǎng)格，進(jìn)行地震波傳播的數(shù)值模擬。

2 震源破裂模型及三維速度結(jié)構(gòu)模型

2.1 震源破裂模型

震后多個(gè)機(jī)構(gòu)給出了本次地震的震源機(jī)制解和破裂過程反演結(jié)果，文中采用了中國地震局地球物理研究所張旭等（https：//www.cea－igp.ac.cn/kydt/278892.html）利用遠(yuǎn)場(chǎng)體波數(shù)據(jù)初步反演的此次地震的震源運(yùn)動(dòng)學(xué)破裂模型。震中為37.702°N、141.587°E，震源深度63.1 km，矩震級(jí)Mw7.4，地震矩M0為1.696×1020Nm，破裂面走向和傾角分別為184°和40°，破裂面沿走向的長度和沿傾向的寬度為110 km和70 km，子斷層尺寸為10 km×10 km。斷層面的破裂滑動(dòng)主要集中于起始破裂點(diǎn)附近，破裂持續(xù)時(shí)間約15 s左右，破裂速度VR=3.0 km/s，其SSW側(cè)的破裂滑動(dòng)更明顯，最大滑動(dòng)量約3 m，該破裂模型充分描述了滑動(dòng)在破裂面上的不均勻分布，滑動(dòng)主要集中于破裂起始點(diǎn)附近區(qū)域，該區(qū)域可看作一個(gè)凹凸體，文中也根據(jù)Somerville等［12］的方法對(duì)破裂面進(jìn)行了裁剪，裁剪了邊緣基本不存在破裂滑動(dòng)區(qū)域（圖2虛線所示）。震源時(shí)間函數(shù)為單周期的Kupper小波函數(shù)［18］，根據(jù)Ekstr?m等［19］給出的震源上升時(shí)間與地震矩經(jīng)驗(yàn)關(guān)系確定震源上升時(shí)間TR，根據(jù)上升時(shí)間得到震源輻射頻帶為1/2TR～2/TR。

圖2 運(yùn)動(dòng)學(xué)破裂模型的震源滑動(dòng)分布Fig.2 Source-slip distribution of kinematic rupture model

2.2 三維速度結(jié)構(gòu)模型

地震波傳播特性受介質(zhì)性質(zhì)的影響，特別是對(duì)地震動(dòng)低頻成分具有顯著影響，因此精細(xì)的三維速度結(jié)構(gòu)模型是準(zhǔn)確模擬研究區(qū)域內(nèi)的地震波傳播過程的基礎(chǔ)。文中模擬了400 km×500 km×100 km范圍內(nèi)的地震波傳播過程，使用橫向墨卡托投影將模型區(qū)域轉(zhuǎn)換為笛卡爾坐標(biāo)系統(tǒng)以進(jìn)行球形測(cè)量［20］，中心坐標(biāo)為140.5°E、37.5°N。地表和海底地形起伏數(shù)據(jù)來自于美國國家海洋和大氣管理局發(fā)布的全球地表地形ETOPO1模型［21］，陸地自由表面以上的網(wǎng)格為空氣層，波速為0，密度為1 kg/m3，為考慮海水層對(duì)地震波傳播的影響，在海域?qū)⒑５字梁Ｆ矫娴姆秶O(shè)定為彈性海水層，海水S波波速為零，P波波速VP=1.5 km/s，密度為103kg/m3，剛度為0。三維地殼速度結(jié)構(gòu)模型采用日本綜合速度結(jié)構(gòu)模型［22］，該模型使用23個(gè)等速層表示三維非均質(zhì)地下結(jié)構(gòu)以及沉積盆地的復(fù)雜基底形狀，包括不規(guī)則基底地形和海底測(cè)深、低速沉積層分布、康拉德不連續(xù)面和莫霍面深度、太平洋板塊邊界、洋－殼與洋－幔邊界等速度不連續(xù)界面等，能夠用于評(píng)價(jià)日本大地震長周期強(qiáng)地震動(dòng)。研究區(qū)域的S波速度模型如圖3所示，圖中可見，速度結(jié)構(gòu)橫向不均勻、板塊交界處的顯著差異等特點(diǎn)，由于板塊交界，計(jì)算區(qū)域內(nèi)西南部速度結(jié)構(gòu)較為復(fù)雜，存在非均勻的殼內(nèi)低速體，即分割不同板塊的局部低速層。文中建立了包含復(fù)雜地表地形、海水層、板塊邊界等的交錯(cuò)網(wǎng)格，模型邊界選擇完美匹配層PML厚度為5 km（20個(gè)網(wǎng)格）［23－24］作為吸收邊界條件。

圖3 三維S波速度結(jié)構(gòu)模型透視圖Fig.3 Perspective view of model with three-dimensional S-wave velocity structure

為保證數(shù)值模擬精度，根據(jù)有限差分法的穩(wěn)定性條件，三維模型的空間離散網(wǎng)格單元尺寸設(shè)定為250 m×250 m×250 m，計(jì)算區(qū)域內(nèi)共劃分了14.4億個(gè)單元，時(shí)間步長設(shè)定為0.015 s，模擬了250 s內(nèi)的地震波傳播過程。沉積層最小剪切波速為500 m/s，根據(jù)網(wǎng)格尺寸、最小波速等，模擬的最高頻率fmax為0.8 Hz。本研究利用了國家超級(jí)計(jì)算天津中心天河二號(hào)超算，通過MPI并行來提高計(jì)算效率，實(shí)現(xiàn)大尺度區(qū)域內(nèi)的地震波傳播的數(shù)值模擬。

3 討論

3.1 速度波場(chǎng)快照

圖4顯示了不同時(shí)刻的速度波場(chǎng)快照。20 s時(shí)，地震波沿?cái)鄬与p向輻射，P波波峰到達(dá)斷層北部的沿海地區(qū)，地震波傳播相對(duì)較為規(guī)則，主要由震源破裂模型控制，破裂滑動(dòng)向東側(cè)滑沖，未表現(xiàn)出明顯的方向性效應(yīng)；35～50 s左右，強(qiáng)烈的P波及S波抵達(dá)福島東部及宮城縣地區(qū)，產(chǎn)生強(qiáng)烈地面運(yùn)動(dòng)，在震中東南部海域出現(xiàn)更大的速度峰值。60～80 s，地震波繼續(xù)傳播至關(guān)東、新潟地區(qū)，受地表地形的影響，波場(chǎng)呈邊緣不規(guī)則的圓形，并存在局部起伏地形引起的幅值放大，體波抵達(dá)后在關(guān)東、新潟地區(qū)產(chǎn)生了長周期面波，東南部海底出現(xiàn)了破碎的散射波，關(guān)東地區(qū)低速沉積層引起波形滯后現(xiàn)象，且增大了地震波的幅值，傾斜入射角使得體波轉(zhuǎn)換為面波，體波與面波的相互干涉在沉積邊緣產(chǎn)生局部放大效應(yīng)；100 s時(shí)，關(guān)東、新潟地區(qū)的面波波長變短、傳播速度下降，同時(shí)觀測(cè)到關(guān)東、新潟地區(qū)速度場(chǎng)的持續(xù)放大，在模型邊界吸收層未發(fā)現(xiàn)顯著的人工邊界反射波。模擬波場(chǎng)快照較好地體現(xiàn)了地震波的傳播特征。

圖4 模擬速度波場(chǎng)快照Fig.4 Snapshots of the simulated velocity wavefiled

東南部海底散射波的存在表明海底高速介質(zhì)層間反射和低速軟流層對(duì)入射地震波的吸收折射作用［25－27］；Noguchi等［28－29］分別研究了瑞利波在海洋區(qū)域的傳播及轉(zhuǎn)換，瑞利波對(duì)海水層的厚度和海底地形起伏等結(jié)構(gòu)非均質(zhì)性較為敏感，說明了日本海溝及海水層對(duì)海底瑞利波傳播的影響；關(guān)東、新潟等地區(qū)的滯后性放大現(xiàn)象表明深厚沉積場(chǎng)地對(duì)地震波傳播的影響，長周期面波在盆地內(nèi)低速沉積物中持續(xù)傳播，造成局部速度波場(chǎng)的振幅放大和持續(xù)時(shí)間明顯增長［30－31］。

3.2 速度波形對(duì)比

圖5對(duì)比了沿海地區(qū)固定軌跡AA′（大致平行于斷層跡線）上部分臺(tái)站的模擬與觀測(cè)記錄的三分量速度時(shí)程，其中對(duì)觀測(cè)記錄依次進(jìn)行零線校正、記錄波形首位加余弦窗并補(bǔ)零、巴特沃斯非因果帶通濾波（0.01～0.8 Hz）處理。觀測(cè)記錄速度波形中存在一個(gè)主要波峰，這與斷層上發(fā)生了一次主破裂相對(duì)應(yīng)，模擬記錄速度波形也顯示了一個(gè)主要波峰。震中距較大的臺(tái)站的模擬與觀測(cè)速度波形具有較為一致的振幅與震相，但近斷層處的臺(tái)站（例如MYG017、FKS004、FKS001等）模擬記錄幅值略偏小。近斷層地震動(dòng)高頻成分豐富，并受發(fā)震斷層的規(guī)模、破裂特點(diǎn)等因素影響，表現(xiàn)為近斷層強(qiáng)震動(dòng)的集中性等特征［32－33］，從而使近斷層觀測(cè)速度波形幅值增加；震中距較遠(yuǎn)的關(guān)東及千葉地區(qū)臺(tái)站速度波形出現(xiàn)與地形相關(guān)的較為復(fù)雜的面波持時(shí)和幅值的放大，圖4中同樣存在局部放大現(xiàn)象，表明該地區(qū)地震波在深厚沉積場(chǎng)地發(fā)生復(fù)雜的散射和波形轉(zhuǎn)換。

圖5 （續(xù)）Fig.5（Continued）

圖5 觀測(cè)與模擬記錄速度波形對(duì)比Fig.5 Comparisons of the observed and the simulated velocity waveforms

3.3 峰值參數(shù)對(duì)比

圖6對(duì)比了K－NET臺(tái)網(wǎng)觀測(cè)與模擬的峰值地面速度PGV，二者具有相似的隨斷層距的衰減趨勢(shì)，在近斷層區(qū)域模擬值普遍偏小，如東西向近斷層臺(tái)站（RJB<100 km）模擬記錄PGV約2～8 cm/s，而觀測(cè)記錄PGV約6～20 cm/s，圖7對(duì)比了觀測(cè)和模擬記錄PGV的空間分布，模擬記錄PGV空間分布與觀測(cè)記錄較為一致，在關(guān)東、新潟地區(qū)均存在明顯的地震動(dòng)放大。

圖6 觀測(cè)與模擬記錄的PGV隨距離衰減Fig.6 Attenuation of PGVs versus distance for the observed and simulated ground motions

圖7 觀測(cè)記錄和模擬記錄PGV空間分布Fig.7 Spatial distributions in terms of PGVs for the observed and simulated ground motions

3.4 傅里葉幅值譜對(duì)比

圖8為部分近場(chǎng)及遠(yuǎn)場(chǎng)臺(tái)站速度記錄的傅里葉幅值譜，其中近斷層臺(tái)站（MYG015、FKS004）模擬值與觀測(cè)值具有一致的峰值頻帶（分別約為0.4、0.3 Hz）且三分量譜值較為相近，而遠(yuǎn)場(chǎng)臺(tái)站（CHB007、CHB017）擁有較寬峰值頻帶（約0.15～0.7 Hz）的同時(shí)，水平向模擬譜值略小于觀測(cè)結(jié)果。無論是近場(chǎng)或遠(yuǎn)場(chǎng)臺(tái)站，大于0.1 Hz的范圍內(nèi)，速度記錄與模擬結(jié)果的譜值具有較好的一致性，部分臺(tái)站在小于0.1 Hz的超低頻域內(nèi)，模擬值與觀測(cè)值存在一定的偏差，這可能與記錄的有效頻帶有關(guān)。

圖8 部分臺(tái)站觀測(cè)和模擬的速度記錄的傅里葉幅值譜對(duì)比Fig.8 Comparisons between the Fourier velocity spectra of the observed and simulated recordings in selected stations

模擬與觀測(cè)記錄總體具有較好的一致性，而針對(duì)模擬記錄偏小的情況，我們分析可能有以下幾點(diǎn)原因：

（1）模擬中采用的數(shù)值網(wǎng)格尺寸有限，為保證滿足運(yùn)算的波長條件（波長內(nèi)包含5～10個(gè)離散網(wǎng)格），運(yùn)算中截止速度Vcut=500 m/s，即陸地中波速小于500 m/s的松軟表層場(chǎng)地均按500 m/s的堅(jiān)硬場(chǎng)地進(jìn)行替代，因此并未考慮近地表更軟弱的土層對(duì)地震動(dòng)的影響。而在海底未進(jìn)行此項(xiàng)替換，因而較好的模擬了由沉積物產(chǎn)生的豐富的面波。

（2）數(shù)值模擬選擇完美匹配層吸收邊界條件，假設(shè)吸收邊界區(qū)域介質(zhì)為完全彈性，此方法在介質(zhì)模型中波速比較大的界面處極偶爾會(huì)出現(xiàn)不穩(wěn)定性［18，30］。為保證結(jié)果的穩(wěn)定，運(yùn)算中忽略速度結(jié)構(gòu)中的波速比過大界面處的低速層，從而使PGV的空間分布中模擬結(jié)果偏小。

（3）低估體波速度振幅可能與考慮了海水層的三維速度模型有關(guān)，海水層相較于剛性基巖，作為地下結(jié)構(gòu)中的低速體能夠吸收一部分地震波，從而降低陸地上的地震波振幅。Maeda［34］曾模擬了2011年日本東北311地震的三維地震動(dòng)響應(yīng)，并對(duì)比了模型中是否含有海水層的模擬結(jié)果，發(fā)現(xiàn)含海水層的結(jié)果較之不含海水層的結(jié)果偏小。

4 結(jié)論

文中利用2022年3月16日日本7.4級(jí)地震的震源運(yùn)動(dòng)學(xué)破裂模型初步反演結(jié)果，采用綜合考慮了地形起伏、海水層、板塊邊界等的三維速度結(jié)構(gòu)精細(xì)化模型，基于交錯(cuò)網(wǎng)格有限差分方法，模擬了地震波的傳播過程，得到以下主要結(jié)論：

（1）計(jì)算區(qū)域內(nèi)的地震動(dòng)顯示出明顯的區(qū)域變化特征。模擬速度波場(chǎng)快照中存在局部地形引起的放大效應(yīng)；同時(shí)在含有沉積層的關(guān)東、新潟等地區(qū)呈現(xiàn)顯著的盆地效應(yīng)，出現(xiàn)較為復(fù)雜的長周期面波及放大的地震波振幅；海底持時(shí)較長的散射波的存在則顯示了海底高速介質(zhì)與地表介質(zhì)的層間反射和低速軟流層對(duì)入射地震波的吸收折射作用，體現(xiàn)了計(jì)算區(qū)域內(nèi)的地表地形起伏、深厚沉積層和海底結(jié)構(gòu)等復(fù)雜地形對(duì)地震波傳播的影響。

（2）數(shù)值模擬結(jié)果與K－NET臺(tái)站的實(shí)際觀測(cè)記錄具有較為一致的幅值與震相。通過對(duì)比模擬記錄與觀測(cè)記錄PGV的空間分布，模擬記錄具有相似的空間分布特征但峰值速度略小，可歸因于計(jì)算模型中的海水層一定程度上降低了模擬的地震波振幅，突出了海水層作為速度模型中的低速層對(duì)地震波的吸收作用。

在地震波三維數(shù)值模擬過程中，地下速度結(jié)構(gòu)的精度及網(wǎng)格劃分的疏密程度對(duì)波的傳播過程影響較為顯著，文中所采用的三維速度模型已相對(duì)比較精細(xì)，但模擬結(jié)果相對(duì)觀測(cè)記錄仍偏小。可知除速度模型以外的影響因素，包括斷層分布等震源參數(shù)的設(shè)置和由于有限差分法數(shù)值運(yùn)算的局限性而對(duì)三維速度模型進(jìn)行的處理等，都將影響數(shù)值模擬結(jié)果的準(zhǔn)確度，需進(jìn)一步進(jìn)行分析，以得到更加符合實(shí)際的模擬結(jié)果。

致謝：感謝東京大學(xué)地震研究所Takashi Furumura教授對(duì)數(shù)值模擬過程的指導(dǎo)。中國地震局地球物理研究所張旭提供日本本州近海M7.4地震的震源破裂過程反演結(jié)果。K－NET臺(tái)站的強(qiáng)震動(dòng)觀測(cè)數(shù)據(jù)來自日本科學(xué)技術(shù)研究所（https：//www.kyoshin.bosai.go.jp/kyoshin/search/index_en.html，NIED）。