2022年3月16日日本近海7.4級地震地震動場三維有限差分模擬

李春果，王宏偉，溫瑞智，強生銀，任葉飛

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

引言

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

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

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

1 數(shù)值模擬方法

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

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

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

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

2.1 震源破裂模型

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

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

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

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

3 討論

3.1 速度波場快照

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

圖4 模擬速度波場快照Fig.4 Snapshots of the simulated velocity wavefiled

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

3.2 速度波形對比

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

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

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

3.3 峰值參數(shù)對比

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

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

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

3.4 傅里葉幅值譜對比

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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