2022年1月8日青海門源6.9級強地震動模擬研究

師涵博 張元生

摘要：為滿足數字化時代應急處置的新要求，基于隨機有限斷層法的發展，利用青海門源地區地下三維速度結構的模型和vS30數據，通過逐步迭代射線追蹤法、格林函數位移解析相位譜和有限斷層法的聯合計算，以青海門源6.9級地震為例，獲得具有地表土層放大效應的強地震動模擬數據，進而繪制出研究區內PGA和烈度分布模擬圖；并與實際臺站記錄的PGA和實地調查烈度結果相比，其烈度區劃范圍基本一致，同時也驗證此聯合計算方法可用于未來地震災害的快速評估，并為災后應急救援提供參考。

關鍵詞：門源地震; 強地震動模擬; 有限斷層法; 逐步迭代射線追蹤法

中圖分類號： P315.3+3 ?????文獻標志碼：A?? 文章編號： 1000-0844（2024）03-0680-12

DOI：10.20000/j.1000-0844.20230920001

Strong ground motion simulation of the Qinghai M6.9 earthquake on January 8， 2022

SHI Hanbo， ZHANG Yuansheng

（Lanzhou Institute of Seismology， CEA， Lanzhou 730000， Gansu， China）

Abstract：?To address the modern requirements of emergency response efforts in the digital age， we developed an approach based on the stochastic finite fault method. Specifically， utilizing an underground three-dimensional velocity structure model and vS30 data from the Menyuan region of the Qinghai Province， we applied the step-by-step iterative ray-tracing method， phase spectrum of Green's function displacement analytical solution， and finite fault method to record strong ground motion simulation data incorporating ground surface amplification effects. Furthermore， considering the Menyuan M6.9 earthquake that struck the Qinghai Province on January 8， 2022， we created simulation maps depicting the peak ground acceleration （PGA） and intensity distributions within the study area. Subsequently， comparing the simulated results with actual PGA records from monitoring stations and field survey intensity data revealed consistent intensity zonation scopes. Moreover， these results validated the utility of the proposed method for rapid assessments of future earthquake disasters， offering valuable insights for postdisaster emergency rescue efforts.

Keywords：Menyuan earthquake; strong ground motion simulation; finite fault method; step-by-step iterative ray-tracing method

0 引言

根據地震臺網中心測定，北京時間2022年1月8日1時45分27秒，青海省海北藏族自治州門源回族自治縣（37.77°N，101.26°E）發生MS6.9地震（下文簡稱為門源地震），震源深度10 km，這是門源地區自1986年8月26日MS6.4地震和2016年1月21日MS6.4地震之后發生的又一次強震活動［1］。此次門源地震的發震構造為帶有少量逆斷分量的左旋走滑斷層，位于青藏高原東北緣冷龍嶺斷裂、托萊山斷裂和肅南—祁連斷裂（俄堡段）的階區部位，破裂帶呈NWW—SEE走向［2］，區域構造較為復雜，如圖1所示。其中，冷龍嶺斷裂長約120 km，具有明顯的線性特征，屬于全新世活躍的左旋走滑斷裂；而托萊山斷裂則位于河西走廊南部邊緣，走向NNW，全長280 km，是一條晚更新世期間活動的左旋走滑斷裂［3］。此外，祁連山斷裂帶是青藏高原東北部最為活躍的斷裂帶之一，尤其是海原斷裂帶，在歷史上曾發生過多次MW7～8的破壞性地震。在過去約35年間，門源地區已發生了3次MW>6的地震，而且這些地震的主要發震斷層均位于海原斷裂帶中冷龍嶺斷裂西段［4］。頻繁的地震活動對抗震救災的時效性提出了更高的要求，因此快速、準確地模擬地震所產生的強地面運動及其地表響應具有重要意義。本文擬以門源MS6.9地震為例，進行強地面運動的模擬及烈度評估，并在此基礎上展開相關的討論和分析。

在進行地震波場模擬時，射線追蹤是一個繞不開的話題，相較于波動方程法可能產生的頻散以及積分方程法計算的復雜性，其原理簡單清晰，無需近似等價模型的建立與調整，對波場衰減的模擬更具有簡便性和合理性。本研究首先基于張元生［5］、高爾根等［6］提出的逐步迭代射線追蹤法進行射線路徑的計算，該方法除了追蹤路徑結果符合射線追蹤要求外，還具有計算速度快、精度高且可以克服射線路徑非唯一性的特點;接著將此方法與Xin等［7］提供的中國大陸巖石圈三維速度結構模型相結合，同時運用孫曉丹［8］基于隨機有限斷層法改進的格林函數位移解析法和王衛民（私人通信）利用遠震波形數據反演的斷層破裂模型，在頻率域考慮震源、傳播路徑及場地效應等影響，計算出相應的傅里葉幅值譜和格林函數位移解析相位譜，并將兩者結合，得到相應的地震時程數據，進而得到研究區域的地震動峰值加速度（Seismic Peak Ground Acceleration，PGA）及儀器烈度（Instrument Intensity，II）分布圖。

此方法的合理性和穩定性已在文獻［8］中有詳細論述。本文在此基礎上進行了射線路徑計算的優化以及參數完備性的補充與討論，使此方法的原理更為清晰，方便后續的研究和改進，可為準確模擬地震災害分布以及防震減災工作提供一定的理論支持與科學依據。

1 理論與方法

1.1 三維速度結構下的射線追蹤法

隨著地震學以及各種觀測手段的不斷發展，人們對地下介質結構及其性質有了更加清晰的認識，使得準確、快速地模擬強地震動成為可能。本文旨在模擬近場范圍內直達波產生的地震動分布，因此相較于效率不高的波動方程法［9］，采用逐步迭代射線追蹤法［5-6］為基礎，在三維不均勻塊狀速度結構中進行射線追蹤的難點主要有：當介質橫、縱向都不均勻時，對于射線上某一點的調整較為復雜;除了兩種軸交點及兩種網格交點的基本情況外，還需要考慮各種特殊情況，如算得的解是否可用、射線垂直或水平入射的參數確定，以及尋找最優解等問題。下文給出相關原理的簡單推導。

在設立的三維速度結構的坐標系O-XYZ中，O、A、B、C分別表示三維速度結構與地表水平面的四個交點，Z軸垂直向下。S和R分別為隨機選取的震源地表投影點（震中）和接收點，如圖2所示。

為了簡化計算，上述點的坐標都以經、緯度來表示。還需明確X、Y、Z軸分別被速度所劃分的刻度值，均以km為單位，從0開始，到對應的速度結構邊界結束。在開始計算之前，需要確保所有數據都在同一單位下，對于經、緯度向距離的轉換，本文采用Haversine公式，相較于Great-circle distance公式使用大量余弦函數會導致較大的舍入誤差，而其采用正弦函數，即使距離很小，也能保證足夠的有效數字：

D=2R·arctan sin2（dlat）+sin2（dlon）·cos（lat1）·cos（lat2） 1-sin2dlat2-sin2dlon2·cos（lat1）·cos（lat2） （1）

式中：D表示S（lon1，lat1）和R（lon2，lat2）之間的距離;R為地球平均半徑，單位都為km;dlon和dlat分別代表兩點間經、緯度之差的絕對值。

在速度結構確定后，根據Koch［10］、張元生等［11］的研究結果可知：進行三維射線追蹤時，如果忽略橫向折射大約5%的橫向速度擾動，走時結果誤差在0.1 s左右。本研究選取的速度結構橫向不均勻性較小，因此這一誤差是完全可以接受的;據此，在三維速度結構中，本文參考二維速度隨機分布逐步迭代射線追蹤法［6］，自動截取震源點與接收點所在的二維速度結構平面進行不均勻介質下射線路徑的計算，可大大減少計算量和計算時間。此方法的原理基于Snell折射定律，如式（2）所示：

sinθ1v1=sinθ2v2 （2）

當射線交點位于橫向軸時，路徑的偏移變化如圖3所示。

從震源點出發，順序選取路徑上的三個點P1、P2、P3進行迭代，P2為中間點，經過一次迭代后，射線路徑由圖3中的虛線校正為實線，即中間點由P2移動至P21，計算出偏移量Δx，即可得到校正后的射線路徑，如式（3）所示：

v1（x3-x2-Δx）（x3-x2-Δx）2+（z3-z2）2=?? v2（x2-x1+Δx）（x2-x1+Δx）2+（z2-z1）2（3）

同理，當交點位于縱向軸時，只需將x、z的坐標調換即可，具體如下：

v1（z2-z3+Δz）（z2-z3+Δz）2+（x3-x2）2=?? v2（z1-z2-Δz）（z1-z2-Δz）2+（x2-x1）2（4）

當交點位于兩軸交點時，需區分上、下偏移兩種情況，分別計算出從點P2到P21和P22所對應的偏移量，如圖4所示。

結合Snell折射定律，當式（5）、（6）聯立無實數解時，求解式（7）、（8）即可。其公式如下：

v1（Δz）Δx2+Δz2=v2（z1-z2-Δz）（z1-z2-Δz）2+（x2-x1）2（5）

v3（Δx）Δx2+Δz2=v2（x3-x2-Δx）（x3-x2-Δx）2+（z2-z3）2 （6）

v1（Δx）Δx2+Δz2=v21（x2-x1-Δx）（x2-x1-Δx）2+（z1-z2）2（7）

v3（Δz）Δx2+Δz2=v21（z2-z3-Δz）（z2-z3-Δz）2+（x3-x2）2 （8）

在迭代計算過程中，首先計算整條射線路徑在實數域內的偏移量Δx和Δz;再對其檢查是否存在新的或者可以合并的交點;完成優化后，進行下一輪的迭代。在本研究中，我們規定當整條射線路徑迭代15次或者最大偏移量低于5 m 時，認為射線路徑已達到要求，不再執行計算。

1.2 基于解析格林函數相位的有限斷層法

對射線路徑優化后，強地震動模擬的準確性可有一定提升。除此之外，強地震動場模擬的另一重要內容是準確計算傅里葉幅值譜與相位譜。在以往的研究中，很多使用隨機合成法得到的平均相位譜，使得各點地震動時程間幾乎不存在相關性，地震動場的空間特征無法表達。因此，本研究參考孫曉丹［8］提出的基于解析格林函數相位的有限斷層法來進行計算，其主要用格林函數位移解析相位譜代替了之前的平均白噪聲隨機相位譜，并充分論證了該方法的可行性與一致性。此方法還具有計算時間短、結果穩定等優點。具體實現過程如下：

1.2.1 傅里葉幅值譜

傅里葉幅值譜的計算需在頻率域將震源譜S（M0，f）、傳播路徑衰減P（R，f）、場地效應G（f）、地震動類型因子I（f）=（2πf）n（n取值0、1、2分別表示位移、速度和加速度）進行連乘，獲取最終的幅值譜FA（M0，f，R），如式（9）所示：

FA（M0，f，R）=S（M0，f）·P（R，f）·G（f）·I（f）（9）

其中，震源模型采用動拐角頻率模型［12］，Sij（M0，f）為子源震源譜，可表示為：

Sij（M0，f）=CMOijHij1+ffcij2 （10）

式中：C為表達輻射方向性差異的常數;f表示頻率;MOij、fcij和Hij分別表示第（i，j˙）個子源的地震矩、動拐角頻率以及高頻標度因子。其中，C的表達式可寫為：

C=RθφFV4πρc3 （11）

式中：F為自由表面放大，可由視入射角計算得到;V為地震波水平分量含剪切波能量的比例，通常取為0.707;ρ表示震源處的介質密度;c為P波或S波的速度;Rθφ為點源輻射模式，隨方位角呈對稱性變化［13］。參考Onishi等［14］的研究，其可由滑動角及斷層的傾角確定，P、SV及SH波對應的表達式如下：

RP=415 （12）

RSV=12sin2λ1415+13sin2（2δ）+cos2λ415+23cos2δ（13）

RSH=1223cos2λ（1+sin2δ）+13sin2λ［1+cos2（2δ）］（14）

式中：λ代表子斷層滑動角;δ則代表子斷層的傾角。

為使每個子斷層只觸發一次地震并且遵循地震矩守恒原則，子源地震矩應為各個子斷層滑動量的加權平均值，如式（15）所示：

M0ij=M0DijD （15）

式中：M0為總地震矩［15］，并M0=10（1.5MW+9.1），單位為N·m;Dij、D分別表示第（i，j）個子源的滑動量以及所有子源滑動量的總和。

關于式（10）中fcij的計算方法較多，在此采用孫曉丹［8］改進的公式：

fcij=f01+M0ijM0ave13 （16）

式中：M0ave表示平均地震矩，M0ave=M0/N，N為子斷層總數;f0為根據總地震矩估算的拐角頻率（定拐角頻率），如式（17）所示：

f0=4.9×106βΔσM013 （17）

同理，斷層面上第（i，j）個子源所對應的定拐角頻率可表示為：

f0ij=4.9×106βΔσM0ij13 （18）

式（17）、（18）中：Δσ表示應力降，單位為0.1 MPa;β表示震源區介質的剪切波速。

為了減小高頻地震動的低估，需引入標度因子Hij進行補償［8］，如式（19）所示：

Hij=∫f1+ff0ij22df∫f1+ffcij22df （19）

此外，地震波在地殼介質中的傳播非常復雜，傳播路徑衰減的影響Pij（Rij，f）通常考慮兩個方面：一是幾何擴散GS（Rij）;另一個是非彈性衰減D（Rij，f）。其中的關系可以表示為Pij（Rij，f）=GS（Rij）·D（Rij，f），即幾何擴散僅與距離有關，但其有多種表達形式。本文采用Atkinson等［16］研究中給出的三段式幾何擴散曲線：

GS（Rij）=R-1ij，RijR02（20）

式中：R01=1.5H，R02=2.5H。H代表研究區內的平均地殼厚度。非彈性衰減D（Rij，f）表示由于地球介質對波能的吸收及耗散產生的非彈性衰減，可以寫為：

D（Rij，f）=exp-πfRijQ（f）v （21）

式中：v表示P波或S波的波速;Q（f）為P波或S波的品質因子：Q（f）=Q0fn，Q0（頻率為1 Hz時的Q值） 的大小與介質均勻程度相關聯。根據周連慶等［17］的統計結果，中國大陸的Q值范圍在116～585，n取值范圍在0.3～1.0。對于泊松介質，P波的品質因子Qα和S波的品質因子Qβ關系［18］可表示為：Qα=9Qβ/4，將對應的v和Q（f）代入式（21）即可得到P波和S波對應的非彈性衰減。

式（9）中，場地效應項G（f）可表示為近地表幅值放大因子A（f）與高頻衰減項K（f）的乘積：

G（f）=A（f）·K（f） （22）

式中：K（f）=exp（-πfk0），f為頻率，k0為特定場地與路徑無關的高頻衰減系數。本研究采用稂子平等［19］擬合的高頻衰減系數和地表以下30 m深度內的平均剪切波速度（The Time-Averaged Shear-wave Velocity to 30 m Depth，vS30）之間的關系表達式：k0=-0.042 8×lg（vS30）+0.153 3，由此加入地表土層的影響;A（f）的計算過于復雜，為了簡便，本文采用Boore等［20］計算得出的北美基巖放大系數A′（f），得到最終的G（f）=A′（f）K（f）。

應用上述公式，即可求出不同子源和接收點對應的傅里葉幅值譜。為了與實際記錄有更好的對應，還需將傅里葉幅值譜按視入射角分解為北南（NS）向、東西（EW）向以及垂直（UD）向進行計算，所需P、SV波的真入射角可由本文1.1節中的射線路徑坐標點求出，根據萬永革［18］所做的相關研究，真入射角和視入射角之間存在著如下聯系：

siniP=αβsinlP2， （P波） （23）

tanlS=2cotiPcot2iS， （SV波） （24）

式中：lP、lS分別表示P、SV波的視入射角;iP、iS分別為入射P、SV波的真入射角，α和β分別代表P波和S波的波速，SH波在自由表面的位移為入射SH波的兩倍［18］。結合接收點的方位角（AZ）以及自由表面放大F，P、SV和SH波幅值譜在北東下坐標中的投影可以表示為：

PNS=sin （lP）cos（AZ） （25）

PEW=sin （lP）sin（AZ） （26）

PUD=-cos （lP） （27）

SVNS=cos （lS）cos （AZ+π） （28）

SVEW=cos （lS）sin（AZ+π） （29）

SVUD=-sin （lS） （30）

SHNS=2cosAZ+π2 （31）

SHEW=2sinAZ+π2 （32）

利用式（25）～（32）可以得到不同類型的波在三個方向上的傅里葉幅值譜，將同方向的幅值譜疊加并與相應的相位譜結合，即可得到該方向上的地震波時程數據，下文將圍繞相位譜進行相關討論。

1.2.2 相位譜

根據孫曉丹［8］的研究成果，可使用無限均勻介質中的格林函數位移解析相位譜進行計算，結合位移表示定理，點源引起空間中任意一點的彈性位移可表示為：

Un（x，t）=Mpq（ξ，τ）·Gnp，q（x，t;ξ，τ）（33）

式中：Gnp，q（x，t;ξ，τ）為Green函數，表示在x位置t時刻觀測到的一對τ時刻在p取向上方向相反、在q取向上分開且作用在ξ上的力，在n方向上的位移響應;Mpq（ξ，τ）代表地震矩張量，根據角動量守恒，可以用9對（p=1，2，3;q=1，2，3）力偶組合來表示;Un（x，t）為引起的地表位移響應，表示在距離震源x位置t時刻由震源引起n方向的位移量。因此，式（33）可進一步表示為：

Mpq·Gnp，q =15γnγpγq-3γnδpq-3γpδnq-3γqδnp4πρ·1r4∫rβrατMpq（t-τ）dτ+?6γnγpγq-γnδpq-γpδnq-γqδnp4πρα2·1r2Mpqt-rα- 6γnγpγq-γnδpq-γpδnq-2γqδnp4πρβ2·1r2Mpqt-rβ+?γηγpγq4πρα3·1rpqt-rα-γnγp-δnp4πρβ3·1rrqpqt-rβ（34）

式中：γi表示從震源點指向觀測點的矢量在i方向的方向余弦;ρ為密度;r表示射線路徑長度;τ為上升時間;α、β分別代表P、S波波速;δij=1，i=j0，i≠j。

式（34）可以分為5項：第一項為近場項，與滑動時間函數的積分有關;第二、三項為P、S波中場項，與滑動時間函數直接相關;第四、五項為P、S波遠場項，與滑動速度函數有關。本研究選用鐘形函數作為滑動速度函數：

（t）=0，t<01τ-1τcos2πτt，0≤t≤τ0，t>τ（35）

則對應的滑動時間函數為：

S（t）=0，t<0tτ-12πsin2πτt，0≤t≤τ1，t>τ（36）

滑動時間函數的積分可表示為：

∫rβrατMpq（t-τ）dτ=MpqHt，rβ-Ht，rα（37）

函數Ht，rc的表達式需分兩種情況討論：

當rc<τ時：

Ht，rc=0，td<0t3d6τ-τ4πtd+τ28π2sin2πτtd，0≤td

當rc≥τ時：

Ht，rc=0，td<0t3d6τ-τ4πtd+τ28π2sin2πτtd，0≤td<τt3d6τ-（td-τ）36τ+τ24π2，τ≤td

式中：td為破裂從起始點至接收點的時間延遲，td=t-tr，tr為斷層子源破裂時間。本文以地震學北東下（XYZ）坐標系為基準，求取子源矩張量Mpq，其可由子源地震矩MOij、走向角φ、傾角δ以及滑動角λ來確定［18］。

綜上即可得到含有高頻成分的地表地震動時程Un（x，t），進行傅里葉變換后，便可提取出不同子源與接收點之間所對應的相位譜。

2 參數的選取

本文選用開源且兼容性較強的Python［21］編程語言開發代碼，借助CPU并行計算實現強地震動的快速模擬。在計算之前，需要輸入的參數包括：三維速度結構與地面4個交點（本文1.1節中的O、A、B、C點）的經度和緯度、矩震級、平均應力降（MPa）、需要計算波形的時間范圍（起始時間，終止時間，時間間隔）、研究區內S波品質因子Q（f）中的Q0和n、震源處的介質平均密度（g/cm3），以及研究區地殼厚度（km）等參數。

本文選用門源地區的平均地殼厚度為55 km;根據Xu等［22］使用強震數據反演此次門源地震的相關參數獲取震源處的平均介質密度為2.72 g/cm3;矩震級MW6.6以及研究區各點的vS30數據均從美國地質調查局（United States Geological Survey，USGS）網站獲取;趙燕杰［23］基于臺網地震記錄，利用Atkinson方法反演得到門源地區S波非彈性衰減為Q（f）=137.0f0.834 2;鄭瑞等［24］基于反演模型估算的地震平均應力降約為5.9 MPa;研究區域地下三維速度結構來源于Xin等［7］提供精度為0.5°×0.5°的中國大陸巖石圈速度結構模型USTClitho1.0，進行插值后如圖5所示。

圖5中不同顏色代表的P波速度值與右下角的色標尺相對應，大致范圍在4～6.5 km/s，深度范圍在地下0～20 km。可以看出，隨著深度的增加，此速度結構的P波速度值逐漸增大，但橫向不均勻性相對較小。在本研究中，視地下介質為泊松介質，因此P波波速α與S波波速β的關系應為α=3β。此外，本文斷層模型采用王衛民（私人通信）所做的門源地震斷層破裂模型，斷層面上的靜態滑動分布如圖6所示。

從圖6中可以看出，此模型共有19×9=171個子斷層，分辨率為2 km×2 km，深度范圍在地下1～17 km，可以適用于圖5所示的速度結構范圍;此斷層結構較為復雜，斷層面上不同顏色代表的滑動量大小與左下角的色標尺相對應，其中最大滑移量為470 cm，靠近地表位置的滑移量分布集中且相對較大;子斷層傾角均為87°，但走向角的變化較為顯著，按斷層走向方向可分成5組，依次為：76°、106°、104°、121°和131°，計算地震動時所需的子斷層參數和部分數據如表1所列。

將表1中的參數與本文1.2.2節的內容相結合，可計算得到無限均勻空間中格林函數的位移解析解;然后對其采用離散數據求導，可得到相應的加速度時程;運用快速傅里葉變換，分別提取出位移和加速度的相位譜;與本文1.2.1節對應的傅里葉位移幅值譜和加速度幅值譜相結合，經過傅里葉逆變換即可得到位移與加速度的時程數據。對于同一個接收點，為了避免出現頻譜泄露，將每個斷層子源對其合成的地震動時程在時域進行疊加，最終得到高頻段的地震動時程波形數據。需注意的是，此方法合成的地震記錄與實際地震記錄之間有差異是正常的，因為使用平均隨機相位譜或格林函數相位譜法合成結果代表的是某一地震引起地震動的平均水平。從統計意義上來說，一次地震發生時的地面運動僅是一個樣本的作用結果，并不能代表多次地震引起的地震動的平均水平［8］。

3 結果分析

根據上述的理論方法和斷層破裂模型可計算得到門源地區的PGA的分布。對于II的求取，本文采用王海云等提出的經驗關系式［25］：

lgPGV=-0.397 7+0.691 4lgPGA，（R=0.863 0）（40）

II=2.259 5+1.344 9lgPGA+1.800 8lgPGV，（R=0.991 7）（41）

式中：PGV代表地震動峰值速度;R為擬合優度指數，其值越接近1，相關度越高。根據同震效應的影響范圍，本研究選取的計算區域為（36.8°N～38.8°N，100°E～102.5°E）、單元網格尺寸為 0.1°×0.1°，共計26×21=546個地表接收點，采用克里金插值中的立方插值法，繪制出PGA和II的分布結果如圖7所示。

圖7（a）為數值模擬得到的研究區PGA等值線圖，右側色標尺的不同顏色對應左圖中不同的地震動峰值加速度值。從中可以看出：PGA高于380 cm/s2的范圍較小，且靠近發震斷層的區域有較為明顯的斷層走向分布趨勢;隨著震源距的增加，在PGA大于40 cm/s2的區域，地震動逐漸衰減且趨于斷層四周大致呈橢圓狀均勻分布;整體來看，研究區PGA在震中NE向有局部小尖端、ES向有較為明顯的區域突出。結合上文的論述可以推斷，這種異常分布可能是由斷層滑動分布和區域速度結構導致的射線路徑變化共同引起的。

圖7（b）中不同II值的分布與右側色標尺的顏色相對應。與圖7（a）相比，在靠近發震斷層位置的區域，兩者分布趨勢較為相似;隨著震中距的增加，烈度分布更趨于橢圓狀，除震中ES向外，基本沒有明顯的異常凸起或凹陷。模擬結果顯示，此次門源MS6.9地震的最大烈度為Ⅸ度，且越靠近震中，烈度梯度的變化越大。

在得到數值模擬結果后，還需要將模擬結果與實際觀測結果進行對比，以檢驗方法的可靠性并分析差異產生的原因。因此，本文采用此次門源地震實際臺站記錄的PGA，以及實地調查烈度等震線與圖7的模擬結果進行比對。

將此次門源地震在甘肅省境內的臺站PGA記錄，結合部分模擬PGA數據（虛擬臺站）使用GMT6［26］內置的surface插值后，呈現出的分布趨勢如圖8所示。

從圖8中可以看出，大部分實際臺站位于震中的NE方向，其中H、L分別表示PGA的區域最高值和最低值。由于數據有限，無實際合站的區域使用PGA模擬值進行了填補，因此填補區并不具有對比分析的意義，我們將重點對圖7（a）和圖8右上角的實際合站密集區進行研究。

通過將圖7（a）和圖8震中NE方向等值線所示數值比較可得：圖8中實際臺站位置處的PGA值基本與圖7（a）一致，具體數值如表2所列。圖7（a）NE向的PGA凸起以及右側邊緣處的異常在圖8中對應位置也有所體現，PGA分布趨勢和數值基本都有很好的對應。初步分析可得，本文的模擬結果可以與實際記錄有相似的趨勢分布，且地下速度結構的影響不容忽視。

通過表2數據分析可得，PGA臺站值與模擬值基本一致：兩者比值最大為1.22，最小為0.75。越靠近震中，模擬的PGA值相對偏小，但從具體比值來看，差異都在可接受范圍內;除此之外，根據盧育霞等［27］、史大成等［28］對PGA放大效應與vS30的關系研究可知：地震動峰值加速度的放大程度主要與近地表vS30相關，vS30越大地震動峰值加速度的放大越小。對比分析表2中的后兩列數據可知：vS30的變化對臺站值與模擬值的比值影響并不大，這說明本文所用方法在不同的地表土層放大作用下，可以得到相對準確的地震動峰值加速度。將上述結論與圖7（a）和圖8所示的PGA分布結合分析可得，本文模擬得到的研究區PGA結果具備一定的合理性和可信度。

本研究將圖7（b）與2022年1月11日應急管理部發布的烈度圖（下文簡稱為實際烈度圖）進行了對比分析（圖9）。從中可以看出，模擬烈度與實地調查的烈度范圍基本一致。但在Ⅸ度和Ⅷ度區，兩者在震中NE向略有突出，Ⅶ度和Ⅵ度區相對較為吻合。在震中NW方位，模擬值與實地調查值的差別隨著烈度的減小不斷增大。

為了進一步分析實際烈度與模擬烈度差異形成的原因，本文根據來源于USGS的數據繪制了研究區內的vS30分布，如圖10所示。

從圖10可以看出，研究區內vS30的變化較為劇烈，但高低值分布的區域性較為明顯。高值區基本呈NW—SE方向分布，與圖7的趨勢分布較為一致。根據前人的研究結論［27-28］，vS30值越大的地方，烈度模擬值與實際值的差異越小，且隨著震中距的增加，這種差異越發明顯;圖9中Ⅸ度和Ⅷ度在震中NE向的突出，模擬烈度與實地調查烈度在NW向有差異的區域，基本都位于圖10中的高值區，結合上述研究結論，兩者間不應有如此大的差異，所以應該考慮是否受到了地形因素的影響。

結合圖1來看，震中NW向的構造相對更為復雜，而山谷、河谷地形對地震烈度也有一定的放大效應［29］，因此除了vS30之外，本文烈度結果的差異應該還受到了地形和地區房屋結構的共同影響。另外，本文僅計算了直達P、S波在地表的地震動響應，待后續代碼算法改進及相關數據完善后，將結合多個震例對此進行深入討論研究。

將圖7（b）和圖9結合來看，Ⅷ度以上（PGA≥380 cm/s2）的區域主要沿NWW—SEE向延伸，并與實際烈度范圍大致相同，也與斷層的主要走向方向較為一致。余下的區域大致呈橢圓狀分布，除NW向外，在震源的其余方位均與實地調查烈度區域有很好的對應。在考慮了vS30的影響后，這種差異大概率是由于地形的影響而導致的，除此之外，在實地烈度調查中，人為因素的影響也不容忽視。結合地形、vS30以及PGA和烈度的實際觀測數據可知，本文無論是PGA還是烈度的模擬結果，均與實際記錄和調查有著較好的對應。這說明本研究所采用的模擬方法可以適用于地震發生后快速估計烈度區劃范圍，為震后災害評估和應急救援提供參考，同時也表明地形構造的分布特征對地震PGA和烈度的分布會有較大影響。

4 討論與結論

4.1 討論

綜上所述，本文通過將逐步迭代射線追蹤法與格林函數位移解析相位譜的有限斷層法相結合，以2022年1月8日青海門源MS6.9地震為例，從基本的公式和理論出發，詳細論述了模擬高頻地震動時程的過程以及相關參數的選取等問題，結合青海門源地區地下三維速度結構［7］與遠震數據反演得到的斷層破裂模型，計算了P、S直達波在該區域地表引起的地震動，并將得到的PGA和II分布與實際數據進行了對比，結果如圖7～9所示。

從PGA對比結果來看，模擬結果與實際記錄有相似的趨勢分布，數值上的差異也在合理范圍之內。從烈度對比結果來看，兩者在高烈度區較為吻合，主要沿NWW—SEE向延伸;結合vS30的分布來看，低烈度區由于地形的影響在NW向有較大的差異。對于強地震動模擬來說，其結果由地形、斷層模型、速度結構以及震源模型共同決定，參數多而復雜，采用不同的模型和參數可能會有不同的結果。隨著后續代碼算法與房屋及地形等數據的優化和充實，可對此進一步研究以提高此方法的適用性，使模擬結果與實際情況更加吻合。

4.2 結論

本文以門源地震為例，進行了強地震動數值模擬的理論應用研究，主要獲得以下結論：

（1） 基于前人的研究成果，使用解析格林函數位移解析相位譜代替隨機有限斷層法中的隨機相位譜，并與逐步迭代射線追蹤法相結合，對前人模擬強地震動的方法進行了優化。同時對相關公式進行了詳細論述并成功模擬得到研究區PGA和II的分布，證明了所用方法的可行性。

（2） 通過將模擬的PGA與實際記錄對比分析，結果表明：本文所用的方法可以很好地模擬出與實際類似的趨勢分布;結合表2中的數據，兩者間的數值差異也在可接受范圍之內，充分證明了此方法獲取地震動模擬結果的合理性。

（3） 將圖7（b）與圖9對比分析可得：模擬烈度與實地調查烈度結果基本一致，存在差異的區域分別位于震中NE和NW向;結合圖1和圖10分析，推斷地形對此有一定的影響。盡管如此，相關結果仍可表明本文所用的模擬方法能夠較好地應用于震后烈度范圍的快速估算。

致謝：本文斷層滑動模型數據來源于王衛民研究員，圖件采用Python［21］和GMT6軟件［26］繪制而成，審稿專家為本文的完善提供了建設性修改意見，使本文描述更為完善，特此致謝！

參考文獻（References）

［1］ 左可楨，羅鈞，趙翠萍，等.青海門源地區地震活動時空分布特征和孕震環境研究［J］.地球物理學報，2023，66（4）：1460-1480.

ZUO Kezhen，LUO Jun，ZHAO Cuiping，et al.Spatiotemporal distribution characteristics of seismicity and seismogenic environment in the Menyuan area，Qinghai Province［J］.Chinese Journal of Geophysics，2023，66（4）：1460-1480.

［2］ 萬永革，黃少華，王福昌，等.2022年門源地震序列揭示的斷層幾何形狀及滑動特性［J］.地球物理學報，2023，66（7）：2796-2810.

WAN Yongge，HUANG Shaohua，WANG Fuchang，et al.Fault geometry and slip characteristics revealed by the 2022 Menyuan earthquake sequence［J］.Chinese Journal of Geophysics，2023，66（7）：2796-2810.

［3］ 王遼，謝虹，袁道陽，等.結合野外考察的2022年門源MS6.9地震地表破裂帶的高分七號影像特征［J］.地震地質，2023，45（2）：401-421.

WANG Liao，XIE Hong，YUAN Daoyang，et al.The surface rupture characteristics based on the GF-7 images interpretation and the field investigation of the 2022 Menyuan MS6.9 earthquake［J］.Seismology and Geology，2023，45（2）：401-421.

［4］ WANG J Y，DING L，HE J K，et al.Research of seismogenic structures of the 2016 and 2022 Menyuan earthquakes，in the northeastern Tibetan Plateau［J］.Remote Sensing，2023，15（3）：742.

［5］ 張元生.三維分塊模型速度隨機分布逐次迭代射線追蹤方法［C］// 中國地球物理學會第十四屆學術年會論文集.北京：中國地球物理學會，1998：318.ZHANG Yuansheng.Three-dimensional chunking model velocity random distribution by iterative ray tracing method［C］//Proceedings of the Fourteenth Annual Academic Conference of the Chinese Geophysical Society.Beijing：Chinese Geophysical Society，1998：318.

［6］ 高爾根，徐果明.二維速度隨機分布逐步迭代射線追蹤方法［J］.地球物理學報，1996，39（增刊1）：302-308.

GAO Ergen，XU Guoming.A new kind of step by step iterative ray-tracing method［J］.Chinese Journal of Geophysics，1996，39（Suppl01）：302-308.

［7］ XIN H L，ZHANG H J，KANG M，et al.High-resolution lithospheric velocity structure of continental China by double-difference seismic travel-time tomography［J］.Seismological Research Letters，2019，90（1）：229-241.

［8］ 孫曉丹.強地震動場估計中若干問題的研究［D］.哈爾濱：哈爾濱工業大學，2010.

SUN Xiaodan.Some issues on estimation of strong ground motion field［D］.Harbin：Harbin Institute of Technology，2010.

［9］ 俞偉哲.基于層狀彈性介質的全波地震波場數值模擬［D］.青島：中國海洋大學，2010.

YU Weizhe.Full wave seismic wave field numerical simulation based on the layered elastic medium［D］.Qingdao：Ocean University of China，2010.

［10］ KOCH M.A numerical study on the determination of the 3-D structure of the lithosphere by linear and non-linear inversion of teleseismic travel times［J］.Geophysical Journal International，1985，80（1）：73-93.

［11］ 張元生，李清河，徐果明.聯合利用走時與波形反演技術研究地殼三維速度結構（Ⅰ）：理論與方法［J］.西北地震學報，1998，20（2）：8-15.

ZHANG Yuansheng，LI Qinghe，XU Guoming.Combined inversion technique to study 3-D crustal velocity structure by using seismic wave travel-time and wave form （Ⅰ）：theory and method［J］.China Earthquake Engineering Journal，1998，20（2）：8-15.

［12］ MOTAZEDIAN D.Stochastic finite-fault modeling based on a dynamic corner frequency［J］.The Bulletin of the Seismological Society of America，2005，95（3）：995-1010.

［13］ 溫瑞智，任葉飛，王宏偉，等.強震動記錄分析與應用：蘆山MS7.0地震為例［M］.北京：地震出版社，2017：158.

WEN Ruizhi，REN Yefei，WANG Hongwei，et al.Analysis and application of strong earthquake records：taking Lushan earthquake as an example［M］.Beijing：Seismological Press，2017：158.

［14］ ONISHI Y，HORIKE M.The extended stochastic simulation method for close-fault earthquake motion prediction and comments for its application to the hybrid method［J］.Journal of Structural and Construction Engineering （Transactions of AIJ），2004，69（586）：37-44.

［15］ HANKS T C，KANAMORI H.A moment magnitude scale［J］.Journal of Geophysical Research：Solid Earth，1979，84（B5）：2348-2350.

［16］ ATKINSON G M，MEREU R F.The shape of ground motion attenuation curves in southeastern Canada［J］.Bulletin of the Seismological Society of America，1992，82（5）：2014-2031.

［17］ 周連慶，趙翠萍，修濟剛，等.利用天然地震研究地殼Q值的方法和進展［J］.國際地震動態，2008，38（2）：1-11.

ZHOU Lianqing，ZHAO Cuiping，XIU Jigang，et al.Methods and developments of research on crustal Q value by using earthquakes［J］.Recent Developments in World Seismology，2008，38（2）：1-11.

［18］ 萬永革.地震學導論［M］.北京：科學出版社，2016.

WAN Yongge.Introduction to seismology［M］.Beijing：Science Press，2016.

［19］ 稂子平，俞瑞芳，肖亮，等.局部場地地震動高頻衰減系數估計模型［J］.地震學報，2023，45（5）：919-928.

LANG Ziping，YU Ruifang，XIAO Liang，et al.An estimation model of high frequency attenuation coefficient of ground motion for local site［J］.Acta Seismologica Sinica，2023，45（5）：919-928.

［20］ BOORE D M，JOYNER W B.Site amplifications for generic rock sites［J］.Bulletin of the Seismological Society of America，1997，87（2）：327-341.

［21］ Python Language Reference［EB/OL］.Python Software Foundation，2021.[2023-09-23].Available：https：//docs.python.org/3/.

［22］ XU C Y，ZHANG Y，WANG R J，et al.Application of MEMS data to fast inversion of rupture process：tests with recordings from the IRREEW network［J］.Seismological Research Letters，2023，94（4）：1821-1835.

［23］ 趙燕杰.用Moya方法反演門源地區數字地震臺站場地響應［J］.地震研究，2016，39（增刊1）：101-106，134.

ZHAO Yanjie.Site response of digital seismic stations in Menyuan area inversed by Moya method［J］.Journal of Seismological Research，2016，39（Suppl01）：101-106，134.

［24］ 鄭瑞，王琪，鄒蓉，等.2022年青海門源MW6.6地震InSAR同震形變場與震源特征［J］.地球物理學報，2023，66（8）：3218-3229.

ZHENG Rui，WANG Qi，ZOU Rong，et al.InSAR coseismic deformation monitoring and source characteristics of the 2022 Qinghai Menyuan MW6.6 earthquake［J］.Chinese Journal of Geophysics，2023，66（8）：3218-3229.

［25］ 王海云，李強.震后近斷層震動圖的快速產出研究：以2022年1月8日青海門源地震為例［J］.世界地震工程，2022，38（2）：1-9.

WANG Haiyun，LI Qiang.Study on rapid generation of near-fault shakemaps after an earthquake：a case of Menyuan earthquake on January 8，2022，Qinhai Province［J］.World Earthquake Engineering，2022，38（2）：1-9.

［26］ WESSEL P，SMITH W H F.New，improved version of generic mapping tools released［J］.Eos，Transactions American Geophysical Union，1998，79（47）：579.

［27］ 盧育霞，王良，魏來，等.利用場地表征參數研究岷漳地區地震動相對變化趨勢［J］.防災減災工程學報，2018，38（2）：359-366.

LU Yuxia，WANG Liang，WEI Lai，et al.Study on relative variability of ground-motion in Min-Zhang Region by using representative site parameters［J］.Journal of Disaster Prevention and Mitigation Engineering，2018，38（2）：359-366.

［28］ 史大成，溫瑞智，杜春清.區域性場地vS30及峰值加速度放大系數估算方法［J］.地震工程與工程振動，2012，32（4）：40-46.

SHI Dacheng，WEN Ruizhi，DU Chunqing.Study on regional site vS30 and PGA amplification factor［J］.Journal of Earthquake Engineering and Engineering Vibration，2012，32（4）：40-46.

［29］ 李平，劉紅帥，薄景山，等.汶川MS8.0地震河谷地形對漢源縣城高烈度異常的影響［J］.地球物理學報，2016，59（1）：174-184.

LI Ping，LIU Hongshuai，BO Jingshan，et al.Effects of river valley topography on anomalously high intensity in the Hanyuan town during the Wenchuan MS8.0 earthquake［J］.Chinese Journal of Geophysics，2016，59（1）：174-184.

（本文編輯：張向紅）

基金項目：甘肅省科技重大專項“甘肅重點地區地震預測預警新方法與新技術”（21ZD4FA011）

第一作者簡介：師涵博（1999-），男，碩士研究生，研究方向：固體地球物理學。E-mail：1598139042@qq.com。

通信作者：張元生（1965-），男，研究員，碩士生導師，主要從事地球內部結構及動力學環境研究。E-mail：zhangys@gsdzj.gov.cn。

師涵博，張元生.2022年1月8日青海門源6.9級強地震動模擬研究［J］.地震工程學報，2024，46（3）：680-691.DOI：10.20000/j.1000-0844.20230920001

SHI Hanbo，ZHANG Yuansheng.Strong ground motion simulation of the Qinghai M6.9 earthquake on January 8， 2022［J］.China Earthquake Engineering Journal，2024，46（3）：680-691.DOI：10.20000/j.1000-0844.20230920001