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

師涵博 張元生

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

關(guān)鍵詞：門源地震; 強(qiáng)地震動模擬; 有限斷層法; 逐步迭代射線追蹤法

中圖分類號： P315.3+3 ?????文獻(xiàn)標(biāo)志碼：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 引言

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

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

此方法的合理性和穩(wěn)定性已在文獻(xiàn)［8］中有詳細(xì)論述。本文在此基礎(chǔ)上進(jìn)行了射線路徑計算的優(yōu)化以及參數(shù)完備性的補(bǔ)充與討論，使此方法的原理更為清晰，方便后續(xù)的研究和改進(jìn)，可為準(zhǔn)確模擬地震災(zāi)害分布以及防震減災(zāi)工作提供一定的理論支持與科學(xué)依據(jù)。

1 理論與方法

1.1 三維速度結(jié)構(gòu)下的射線追蹤法

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

在設(shè)立的三維速度結(jié)構(gòu)的坐標(biāo)系O-XYZ中，O、A、B、C分別表示三維速度結(jié)構(gòu)與地表水平面的四個交點(diǎn)，Z軸垂直向下。S和R分別為隨機(jī)選取的震源地表投影點(diǎn)（震中）和接收點(diǎn)，如圖2所示。

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

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分別代表兩點(diǎn)間經(jīng)、緯度之差的絕對值。

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

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

當(dāng)射線交點(diǎn)位于橫向軸時，路徑的偏移變化如圖3所示。

從震源點(diǎn)出發(fā)，順序選取路徑上的三個點(diǎn)P1、P2、P3進(jìn)行迭代，P2為中間點(diǎn)，經(jīng)過一次迭代后，射線路徑由圖3中的虛線校正為實線，即中間點(diǎn)由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）

同理，當(dāng)交點(diǎn)位于縱向軸時，只需將x、z的坐標(biāo)調(diào)換即可，具體如下：

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

當(dāng)交點(diǎn)位于兩軸交點(diǎn)時，需區(qū)分上、下偏移兩種情況，分別計算出從點(diǎn)P2到P21和P22所對應(yīng)的偏移量，如圖4所示。

結(jié)合Snell折射定律，當(dāng)式（5）、（6）聯(lián)立無實數(shù)解時，求解式（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）

在迭代計算過程中，首先計算整條射線路徑在實數(shù)域內(nèi)的偏移量Δx和Δz;再對其檢查是否存在新的或者可以合并的交點(diǎn);完成優(yōu)化后，進(jìn)行下一輪的迭代。在本研究中，我們規(guī)定當(dāng)整條射線路徑迭代15次或者最大偏移量低于5 m 時，認(rèn)為射線路徑已達(dá)到要求，不再執(zhí)行計算。

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

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

1.2.1 傅里葉幅值譜

傅里葉幅值譜的計算需在頻率域?qū)⒄鹪醋VS（M0，f）、傳播路徑衰減P（R，f）、場地效應(yīng)G（f）、地震動類型因子I（f）=（2πf）n（n取值0、1、2分別表示位移、速度和加速度）進(jìn)行連乘，獲取最終的幅值譜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為表達(dá)輻射方向性差異的常數(shù);f表示頻率;MOij、fcij和Hij分別表示第（i，j˙）個子源的地震矩、動拐角頻率以及高頻標(biāo)度因子。其中，C的表達(dá)式可寫為：

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

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

RP=415 （12）

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

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

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

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

M0ij=M0DijD （15）

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

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

fcij=f01+M0ijM0ave13 （16）

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

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

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

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

式（17）、（18）中：Δσ表示應(yīng)力降，單位為0.1 MPa;β表示震源區(qū)介質(zhì)的剪切波速。

為了減小高頻地震動的低估，需引入標(biāo)度因子Hij進(jìn)行補(bǔ)償［8］，如式（19）所示：

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

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

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

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

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

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

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

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

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

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

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

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

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

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）可以得到不同類型的波在三個方向上的傅里葉幅值譜，將同方向的幅值譜疊加并與相應(yīng)的相位譜結(jié)合，即可得到該方向上的地震波時程數(shù)據(jù)，下文將圍繞相位譜進(jìn)行相關(guān)討論。

1.2.2 相位譜

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

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

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

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表示從震源點(diǎn)指向觀測點(diǎn)的矢量在i方向的方向余弦;ρ為密度;r表示射線路徑長度;τ為上升時間;α、β分別代表P、S波波速;δij=1，i=j0，i≠j。

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

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

則對應(yīng)的滑動時間函數(shù)為：

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

滑動時間函數(shù)的積分可表示為：

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

函數(shù)Ht，rc的表達(dá)式需分兩種情況討論：

當(dāng)rc<τ時：

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

當(dāng)rc≥τ時：

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

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

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

2 參數(shù)的選取

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

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

圖5中不同顏色代表的P波速度值與右下角的色標(biāo)尺相對應(yīng)，大致范圍在4～6.5 km/s，深度范圍在地下0～20 km?？梢钥闯?，隨著深度的增加，此速度結(jié)構(gòu)的P波速度值逐漸增大，但橫向不均勻性相對較小。在本研究中，視地下介質(zhì)為泊松介質(zhì)，因此P波波速α與S波波速β的關(guān)系應(yīng)為α=3β。此外，本文斷層模型采用王衛(wèi)民（私人通信）所做的門源地震斷層破裂模型，斷層面上的靜態(tài)滑動分布如圖6所示。

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

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

3 結(jié)果分析

根據(jù)上述的理論方法和斷層破裂模型可計算得到門源地區(qū)的PGA的分布。對于II的求取，本文采用王海云等提出的經(jīng)驗關(guān)系式［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為擬合優(yōu)度指數(shù)，其值越接近1，相關(guān)度越高。根據(jù)同震效應(yīng)的影響范圍，本研究選取的計算區(qū)域為（36.8°N～38.8°N，100°E～102.5°E）、單元網(wǎng)格尺寸為 0.1°×0.1°，共計26×21=546個地表接收點(diǎn)，采用克里金插值中的立方插值法，繪制出PGA和II的分布結(jié)果如圖7所示。

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

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

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

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

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

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

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

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

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

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

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

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

4 討論與結(jié)論

4.1 討論

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

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

4.2 結(jié)論

本文以門源地震為例，進(jìn)行了強(qiáng)地震動數(shù)值模擬的理論應(yīng)用研究，主要獲得以下結(jié)論：

（1） 基于前人的研究成果，使用解析格林函數(shù)位移解析相位譜代替隨機(jī)有限斷層法中的隨機(jī)相位譜，并與逐步迭代射線追蹤法相結(jié)合，對前人模擬強(qiáng)地震動的方法進(jìn)行了優(yōu)化。同時對相關(guān)公式進(jìn)行了詳細(xì)論述并成功模擬得到研究區(qū)PGA和II的分布，證明了所用方法的可行性。

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

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

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

參考文獻(xiàn)（References）

［1］ 左可楨，羅鈞，趙翠萍，等.青海門源地區(qū)地震活動時空分布特征和孕震環(huán)境研究［J］.地球物理學(xué)報，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］.地球物理學(xué)報，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］ 王遼，謝虹，袁道陽，等.結(jié)合野外考察的2022年門源MS6.9地震地表破裂帶的高分七號影像特征［J］.地震地質(zhì)，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］ 張元生.三維分塊模型速度隨機(jī)分布逐次迭代射線追蹤方法［C］// 中國地球物理學(xué)會第十四屆學(xué)術(shù)年會論文集.北京：中國地球物理學(xué)會，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ī)分布逐步迭代射線追蹤方法［J］.地球物理學(xué)報，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］ 孫曉丹.強(qiáng)地震動場估計中若干問題的研究［D］.哈爾濱：哈爾濱工業(yè)大學(xué)，2010.

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

［9］ 俞偉哲.基于層狀彈性介質(zhì)的全波地震波場數(shù)值模擬［D］.青島：中國海洋大學(xué)，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］ 張元生，李清河，徐果明.聯(lián)合利用走時與波形反演技術(shù)研究地殼三維速度結(jié)構(gòu)（Ⅰ）：理論與方法［J］.西北地震學(xué)報，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］ 溫瑞智，任葉飛，王宏偉，等.強(qiáng)震動記錄分析與應(yīng)用：蘆山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］ 周連慶，趙翠萍，修濟(jì)剛，等.利用天然地震研究地殼Q值的方法和進(jìn)展［J］.國際地震動態(tài)，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］ 萬永革.地震學(xué)導(dǎo)論［M］.北京：科學(xué)出版社，2016.

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

［19］ 稂子平，俞瑞芳，肖亮，等.局部場地地震動高頻衰減系數(shù)估計模型［J］.地震學(xué)報，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方法反演門源地區(qū)數(shù)字地震臺站場地響應(yīng)［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］.地球物理學(xué)報，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］ 王海云，李強(qiáng).震后近斷層震動圖的快速產(chǎn)出研究：以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］ 盧育霞，王良，魏來，等.利用場地表征參數(shù)研究岷漳地區(qū)地震動相對變化趨勢［J］.防災(zāi)減災(zāi)工程學(xué)報，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］ 史大成，溫瑞智，杜春清.區(qū)域性場地vS30及峰值加速度放大系數(shù)估算方法［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］.地球物理學(xué)報，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.

（本文編輯：張向紅）

基金項目：甘肅省科技重大專項“甘肅重點(diǎn)地區(qū)地震預(yù)測預(yù)警新方法與新技術(shù)”（21ZD4FA011）

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

通信作者：張元生（1965-），男，研究員，碩士生導(dǎo)師，主要從事地球內(nèi)部結(jié)構(gòu)及動力學(xué)環(huán)境研究。E-mail：zhangys@gsdzj.gov.cn。

師涵博，張元生.2022年1月8日青海門源6.9級強(qiáng)地震動模擬研究［J］.地震工程學(xué)報，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