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

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

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

引言

據日本氣象廳（JMA）測定，2022年3月16日23：36：29（UTC+9）日本本州東岸近海發生7.4級地震，震中位于141.623°E、37.697°N，震源深度為57 km，是2011年3月11日9.1級東日本大地震后日本東部地區最大震級地震。據日本防災科學技術研究所（NIED）震度分布，地震最大震度達到6強，造成了強烈的地面運動、嚴重的工程結構破壞及人員傷亡。此次地震的發震區位于歐亞板塊、太平洋板塊、北美板塊與菲律賓海板塊的交界，該區域板塊間相互作用復雜，太平洋板塊和菲律賓板塊向西俯沖形成狹長的板間地震帶，歐亞板塊沿邊緣東邊界向東移動，形成雙向對沖式匯聚，擁有活躍的淺源和深源強震活動［1－4］（見圖1）。地震發生后，密集的地震動觀測臺網（K－NET和KiK－net）觀測到強烈的地震動及豐富的長周期面波，其中K－NET臺網的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)

日本東部海域地區強震活動十分活躍，研究人員開展了大量的地震波動數值模擬工作，Furumura等［5］、Takemura等［6］分別模擬了不同地震中日本關東盆地產生的長周期地震動，前者分析了盆地內地震波聚焦效應引起的顯著的盆地放大，后者則重點討論了震源參數對長周期地震動模擬的影響；Nakamura［7］等探討了表面地形、海水層對地震動數值模擬的影響；Petukhin等［8］分析了包含海水層的三維速度結構模型對大阪盆地模擬地震動的影響，解釋了地震面波在海洋沉積層中的產生過程，強調了地震波動模擬中考慮海水層的必要性。

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

1 數值模擬方法

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

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

采用有限差分法的OpenSWPC程序相較于有限元法、譜元法等程序求解形式簡單，占用內存小，計算效率高，且整合了模擬過程必要的前處理和后處理程序，并包含了日本全部地區精細化三維速度結構模型。輸入相應的地震震源或有限斷層參數，程序將自動分配并行計算機內存生成非均勻介質模型的交錯網格，進行地震波傳播的數值模擬。

2 震源破裂模型及三維速度結構模型

2.1 震源破裂模型

震后多個機構給出了本次地震的震源機制解和破裂過程反演結果，文中采用了中國地震局地球物理研究所張旭等（https：//www.cea－igp.ac.cn/kydt/278892.html）利用遠場體波數據初步反演的此次地震的震源運動學破裂模型。震中為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。斷層面的破裂滑動主要集中于起始破裂點附近，破裂持續時間約15 s左右，破裂速度VR=3.0 km/s，其SSW側的破裂滑動更明顯，最大滑動量約3 m，該破裂模型充分描述了滑動在破裂面上的不均勻分布，滑動主要集中于破裂起始點附近區域，該區域可看作一個凹凸體，文中也根據Somerville等［12］的方法對破裂面進行了裁剪，裁剪了邊緣基本不存在破裂滑動區域（圖2虛線所示）。震源時間函數為單周期的Kupper小波函數［18］，根據Ekstr?m等［19］給出的震源上升時間與地震矩經驗關系確定震源上升時間TR，根據上升時間得到震源輻射頻帶為1/2TR～2/TR。

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

2.2 三維速度結構模型

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

圖3 三維S波速度結構模型透視圖Fig.3 Perspective view of model with three-dimensional S-wave velocity structure

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

3 討論

3.1 速度波場快照

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

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

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

3.2 速度波形對比

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

圖5 （續）Fig.5（Continued）

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

3.3 峰值參數對比

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

圖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）的同時，水平向模擬譜值略小于觀測結果。無論是近場或遠場臺站，大于0.1 Hz的范圍內，速度記錄與模擬結果的譜值具有較好的一致性，部分臺站在小于0.1 Hz的超低頻域內，模擬值與觀測值存在一定的偏差，這可能與記錄的有效頻帶有關。

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

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

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

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

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

4 結論

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

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

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

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

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