地震動非參數(shù)化譜反演可靠性分析

王宏偉 ，任葉飛 ，溫瑞智 ，周 影

（1.中國地震局工程力學研究所, 黑龍江 哈爾濱 150080；2.中國地震局地震工程與工程振動重點實驗室, 黑龍江哈爾濱 150080）

地震動譜反演是頻域上同時分離地震動震源、傳播路徑和場地效應的高效方法，該方法普遍應用于地震震源特征[1-4]、地殼介質(zhì)衰減特性[5-6]、場地效應[7-8]、地震動不確定性[9]等的研究，這些研究對深入理解地震震源物理、探究地震發(fā)震機理、可靠預測地震動及地震危險性有著重要貢獻.

根據(jù)傳播路徑衰減表示方式的差異，譜反演方法可分為參數(shù)化和非參數(shù)化兩種.Castro等[10]首次提出的非參數(shù)化方法將傳播路徑衰減約束為距離的單調(diào)遞減平滑曲線，既體現(xiàn)了傳播路徑衰減的復雜性，又描述了地震波在傳播中的不斷衰減，解決了參數(shù)化方法利用經(jīng)驗衰減模型會導致品質(zhì)因子反演結(jié)果為極大值甚至負值，違背地震波衰減損耗的物理事實.之后非參數(shù)化方法開始廣泛應用于地震活躍地區(qū)的地震波衰減特征及地震震源特征等的研究，例如日本[2]、新西蘭[11]、歐洲[3,9,12]、中國四川地區(qū)[13]等，相比參數(shù)化方法，近年來非參數(shù)化方法的應用更普遍.

盡管非參數(shù)化譜反演的應用十分普遍，但反演結(jié)果的可靠性需要謹慎處理諸多關(guān)鍵環(huán)節(jié)來確保.譜反演必須考慮震源與場地項之間的權(quán)衡，一般通過選擇合適的參考場地來實現(xiàn)，參考場地的選取顯然會影響反演結(jié)果的可靠性，目前尚無研究針對如何選擇合理的參考場地給出合理建議.根據(jù)反演步驟的差異，非參數(shù)化譜反演具體分為兩步法[10]和單步法[2]，兩種方法的反演結(jié)果是否存在明顯差異并不清楚，也未有研究評價應優(yōu)先選用哪種方法.根據(jù)理論震源模型[14]，震源譜以地震矩和拐角頻率表示，因此震源譜的反演結(jié)果常用于地震震源參數(shù)的估計及震源物理特性研究，地震矩和拐角頻率之間的權(quán)衡可能會影響震源參數(shù)估計的可靠性[15]，評估震源參數(shù)之間的權(quán)衡對其估計的影響十分必要.

針對影響地震動非參數(shù)化譜反演可靠性的關(guān)鍵環(huán)節(jié)，包括參考場地的確定、單步法或兩步法的選擇、震源參數(shù)估計，本文以2016—2017年意大利中部地震序列的地震動非參數(shù)化譜反演作為典型實例，歐洲強震動數(shù)據(jù)中心（engineering strong-motion，ESM）發(fā)布了地震序列八千余組人工處理的高質(zhì)量的地震動記錄，本研究限定記錄峰值地面加速度（PGA）上限（100 cm/s2）、震源距上限（100 km）及臺站（地震記錄數(shù)下限）10 （10），選用了142個臺站收集的78次震級（3.5～6.1）級地震的5 135組觀測記錄用于譜反演研究，截取記錄S波并在其前后增加余弦邊瓣窗口消除截斷誤差，S波到達時刻和結(jié)束時刻分別通過人工識別和能量法[3]確定，最后以兩水平方向S波傅里葉譜均方根代表觀測記錄水平向S波傅里葉譜.本研究選用大量強震動記錄可確保足夠的數(shù)據(jù)冗余，避免因記錄數(shù)過少造成反演結(jié)果的不可靠，通過系統(tǒng)分析反演結(jié)果，提出確保反演可靠性的上述關(guān)鍵環(huán)節(jié)的處理建議，可靠反演結(jié)果有助于地震震源物理的準確認識及地震動和地震危險性的有效預測.

1 非參數(shù)化譜反演方法

第i個地震（震級Mi）中第j個臺站（震源距Rij）觀測到的地震動水平向S波的傅里葉幅值譜（Oij）可表示為

式中：f為頻率；Si(f,Mi)為第i個地震的震源譜；Aij(f,Rij)為傳播路徑衰減；Gj(f)為第j個臺站的場地效應.

兩步法首先從S波傅里葉譜中分離傳播路徑衰減項（第1步），即

式中：Mai(f)為與第i次地震震級大小相關(guān)的量.

之后從傳播路徑衰減修正的S波傅里葉譜中同時分離震源和場地項（第2步），即

單步法直接求解式（1），從S波傅里葉譜中同時分離震源、傳播路徑衰減和場地項，求解時需約束傳播路徑衰減及震源和場地項的權(quán)衡.

2 參考場地選擇

參考場地一般為場地放大恒為1的無風化出露堅硬巖石場地[16]或平均場地放大為1的多個巖石場地[10]或至少一個放大效應已知的場地[17].實際應用中找到理想的出露堅硬巖石場地作為參考場地十分困難，多數(shù)研究采用場地放大效應恒為1的一個或多個巖石類場地作為參考場地[3,9,11].

一般考慮下列條件選取參考巖石場地：1）場地地質(zhì)資料，例如：場地類別、場地等效剪切波速（例如地下30 m土層的等效剪切波速Vs30）等；2）水平/垂直（H/V）譜比曲線，巖石場地的譜比曲線在整個頻段上平坦且幅值較小.為保證反演結(jié)果的穩(wěn)定性和可靠性，Bindi等[17]規(guī)定盡可能選取觀測記錄足夠多且記錄到震級和方位角范圍足夠大的巖石臺站.Pacor等[3]則假定所有臺站均為巖石場地且平均場地放大恒為1，利用譜反演給出所有臺站的場地放大效應，再選取整個頻段上場地放大小于2的場地作為參考巖石場地.Oth等[11]首先對參考巖石場地進行初選，通過震源譜和場地效應的譜反演試算結(jié)果從初選中確定.鑒于理想巖石場地較難獲取，上述這些研究提供了選取參考巖石場地的一系列方法，但這些方法均是不得已而為之，參考場地選取直接影響譜反演的可靠性.

根據(jù)Eurocode 8規(guī)范，本研究的142個觀測臺站中A類場地（巖石場地或其他接近巖石地質(zhì)構(gòu)造的場地，場地表面軟弱層不大于 5 m 且Vs30＞800 m/s）臺站共33個，其中7個為流動觀測臺站，A類場地臺站作為參考場地的備選，因流動觀測臺站建設要求較低，受環(huán)境影響較大，7個流動觀測臺站不納入備選之列.26個備選臺站中，分別取記錄數(shù)≥ 50組、[40～50）組、[30～40）組、[20～30）組、＜ 20 組的臺站分別為7、5、2、4和8個，有詳細鉆孔數(shù)據(jù)的3個（LSS、AQP 和 PSC，Vs30分別為 1 091、836、1 000 m/s），只有LSS臺站的觀測記錄充足（66組），AQP和PSC臺站觀測記錄較少（分別為15組和14組）.

進一步縮小備選范圍，僅以記錄數(shù)不少于50組（73、111、155、105、60、58、157）的 7個臺站（LSS、TOD、FIAM、GNU、MNF、SPM和 MGAB）作為備選，若最終證明均不合適，可適當放寬對備選臺站記錄數(shù)的要求.分別利用 7 個臺站的 PGA≤ 50 cm/s2的觀測記錄計算0.25～20.00 Hz的H/V譜比曲線（見圖1）.FIAM臺站缺乏地形地貌描述且其H/V譜比曲線在1 Hz附近出現(xiàn)較為明顯的峰值；SPM臺站位于陡坡邊緣，H/V譜比曲線幅值大且6.00 Hz附近出現(xiàn)明顯峰值；GNU和MGAB臺站的H/V譜比曲線幅值較大，說明這兩個臺站存在明顯場地放大效應.上述4個臺站均不適合作為參考巖石場地.LSS臺站的觀測記錄方位角范圍廣（大約100°）、記錄震源距范圍較大（＞30 km）、H/V譜比曲線平坦且幅值??；TOD和MNF臺站的H/V譜比曲線平坦且幅值小，此3個臺站可作為初選的參考巖石場地.

圖1 7 個臺站的 H/V 譜比曲線Fig.1 The horizontal-to-vertical spectral ratios at seven stations

本研究以H/V譜比作為3個初選參考場地的場地放大效應，反演給出2016—2017年意大利中部地震序列觀測記錄的傳播路徑衰減、震源譜和場地效應，以LSS、TOD和MNF分別作為參考場地得到的6次典型地震的加速度震源譜如圖2所示，6次地震（1#～6#）為 20 160 824 013 632_ML6.0（2016 年8月 24日 01：36：32發(fā)生的ML6.0級地震，余同）、20 160 825 031 716_ML4.5、 20 160 824 232 205_ML3.8、20 161 026 191 806_ML5.9、 20 161 031 070 545_ML4.1、20 161 103 003 501_ML4.8，震源譜大體上符合ω-2理論震源模型[14]，但不同參考場地得到的震源譜的形狀及幅值差異明顯.LSS和MNF臺站作為參考場地的震源譜形狀相似，LSS臺站作為參考場地的震源譜的幅值更小；TOD臺站作為參考場地的震源譜的形狀明顯不同，低頻段（0.50～2.00 Hz）幅值更大，高頻段（＞10.00 Hz）衰減更快.上述結(jié)果清楚地展示了參考場地選取對反演結(jié)果的顯著影響.

圖2 6個典型地震的加速度震源譜反演結(jié)果Fig.2 The inverted acceleration source spectra of six typical earthquakes

ω-2模型[14]和高頻衰減模型[18]表示反演震源譜，即

式中：C為常數(shù)，定義見文獻[4]；M0為地震矩；fc為拐角頻率.

網(wǎng)格搜索方法估計地震的震源參數(shù)（M0、fc）及高頻衰減參數(shù)κ，根據(jù)M0計算矩震級Mw，利用M0、fc計算地震應力降Δσ[14].

LSS、MNF和TOD臺站分別作為參考場地時，利用反演震源譜得到的Mw估計值與INGV正式發(fā)布的地方震級ML的差值如圖3（a）所示，整體上Mw估計值依次增大，LSS和MNF臺站作為參考場地的Mw估計值接近，明顯低于TOD臺站作為參考場地的結(jié)果.ML＜ 5.3 時，隨震級增大，Mw估計值與ML之差整體上呈下降趨勢，Mw估計值從明顯高于ML逐漸變化為低于或接近ML；ML＞5.3 時，Mw估計值普遍低于ML，且隨ML增大，Mw與ML的差異逐漸縮小.圖3（a）同時比較了TDMT和MTS提供的Mw與ML之差.TDMT和MTS提供的Mw值普遍低于ML，LSS、MNF和TOD臺站作為參考場地的Mw估計與TDMT、MTS提供的Mw存在一定偏差，尤其是小震的矩震級估計偏大，這可能與距離相關(guān)的能量法截取的小震記錄S波窗口過大有關(guān)，但TOD作為參考場地的Mw估計與TDMT、MTS提供的Mw的偏差最顯著.

圖3（b）、（c）給出了 LSS、MNF 和 TOD 臺站作為參考場地的fc和Δσ估計值，參考場地不同時fc和Δσ估計值差異明顯，LSS作為參考場地的fc估計值明顯低于MNF和TOD作為參考場地的結(jié)果，相應地LSS作為參考場地的Δσ估計值也最小.參考場地不同時M0-fc和Δσ-Mw的變化趨勢較為一致，Mw＜ 5.00時，Δσ估計值隨Mw的增大而增大，LSS、MNF和TOD作為參考場地的Δσ估計值為0.11～2.90、0.32～9.97 和 0.55～13.45 MPa；Mw≥5.00時，Δσ估計值隨Mw的變化不明顯，LSS、MNF和TOD作為參考場地的Δσ為5.40～12.19、12.96～26.50和28.57～78.84 MPa，明顯高于全球板間地震平均應力降（3.31 MPa）[19]及 Bindi等[17]估計的 1997—1998年Umbria-Marche地震序列（本研究選用地震序列的西北部）中15次Mw=4.20～6.00余震的應力降（圖3（c））.Pacor等[3]利用非參數(shù)化譜反演確定的2009年L’Aquila地震序列（本研究地震序列的南部）應力降為 0.25～25.00 MPa，與本研究選用LSS、MNF臺站作為參考場地得到的Δσ估計值的范圍較為一致（圖3（c））.

LSS、MNF和TOD臺站作為參考場地的κ估計值依次增大（圖3（d）），κ估計值近似服從正態(tài)分布，LSS和MNF作為參考場地的κ估計值較為接近，均值分別為 0.017 s和 0.023 s，明顯低于 TOD作為參考場地的κ估計值（均值 0.056 s）.κ估計值反映了參考場地中未予考慮的可能由于近地表巖層風化等引起的高頻衰減[20].一般認為巖層風化越嚴重，其κ值越大.Oth等[2]利用日本井下鉆孔記錄反演得到的井下鉆孔臺站的κ均值為0.015 s，LSS和MNF臺站作為參考場地的κ均值與日本井下鉆孔臺站κ均值接近，低于一般巖石場地的κ值（0.020～0.040 s），說明 LSS 和 MNF 臺站近地表巖層風化不嚴重，而TOD作為參考場地的κ均值較高則說明該臺站近地表巖層風化可能較為嚴重.

圖3 地震震源參數(shù) （ Mw, fc and Δσ ）及高頻衰減參數(shù) （ κ ）估計值Fig.3 Estimated source parameters (Mw, fc and Δσ) and high-frequency decay parameters (κ)

通過分析分別選用LSS、MNF和TOD臺站作為參考場地的震源譜、震源參數(shù)及高頻衰減參數(shù)的試算結(jié)果，LSS和MNF臺站均可作為參考巖石場地的合理選擇.我們建議在選取參考場地時，首先綜合考慮場地地質(zhì)資料、觀測記錄H/V譜比曲線幅值和形狀初步選擇參考巖石場地，再分別以初選作為參考場地進行譜反演，根據(jù)震源譜、震源參數(shù)、高頻衰減參數(shù)的試算結(jié)果，最終選取合適的參考場地.

3 單步法與兩步法

本研究分別采用兩步法和單步法完成2016—2017年意大利中部地震序列地震動記錄的非參數(shù)化譜反演，均以LSS臺站作為參考場地.兩種方法分別確定的0.25～20.00 Hz傳播路徑衰減差異明顯（圖4（a）），差異大小與震源距、頻率顯著相關(guān)（圖4（b）、（c））.近場（＜ 20 km）或低頻 S 波（＜ 2.00 Hz）的傳播路徑衰減較為一致，高頻 S 波（＞5.00 Hz）傳播路徑衰減的差異隨震源距增大迅速擴大，高頻S波（＞5.00 Hz）在遠場（＞60 km）的傳播路徑衰減的差異顯著，兩者相差可達5倍.與單步法相比，兩步法確定的傳播路徑衰減更強烈（數(shù)值更?。?采用幾何擴散和非彈性衰減簡單表示傳播路徑衰減曲線，利用奇異值分解方法給出兩種方法確定的傳播路徑衰減的幾何擴散指數(shù)n和品質(zhì)因子Q，兩步法和單步法的n值較為一致，分別為 -1.01和 -1.07，品質(zhì)因子的差異明顯，分別為 96.45f0.59和 212.45f0.61，說明兩種方法路徑衰減差異主要為頻率相關(guān)的非彈性衰減差異.

圖4 基于兩步法和單步法反演路徑衰減的差異Fig.4 Differences between path attenuations inverted by the two-step and one-step methods

利用單步法和兩步法分別確定的6次典型地震的加速度震源譜如圖5所示，兩步法的震源譜幅值普遍高于單步法，大于2.00 Hz的高頻段尤為明顯.基于反演震源譜，利用網(wǎng)格搜索法估計地震震源參數(shù)（Mw、fc）及高頻衰減參數(shù)（κ）（見圖6（a）～（c）），兩種方法Mw估計值相近，多數(shù)地震Mw估計值相差不超過0.02；兩種方法fc估計值差異較小，單步法fc估計值整體上略小于兩步法，多數(shù)地震單步法fc估計值約為兩步法的88%～98 %，我們還注意到fc估計值的差異與fc大小有關(guān)，fc估計值的差異近似隨fc增大而增大；兩種方法κ估計值差異并不明顯，單步法κ估計值普遍高于兩步法，多數(shù)地震單步法κ估計值約為兩步法的1.0倍～1.5倍，只有極少數(shù)κ估計值＜ 0.01 s的地震，單步法κ估計值約為兩步法的2.0倍～6.0倍.利用Mw和fc估計值計算得到Δσ（如圖6（d））和輻射地震能量Es（如圖6（e）），兩種方法Δσ、Es估計值差異不明顯，單步法Δσ和Es估計值均小于兩步法，且大多數(shù)地震單步法Δσ和Es估計值約為兩步法的70%～90%.

圖5 單步法和兩步法反演給出的6次典型地震的加速度震源譜Fig.5 Acceleration source spectra for six typical events derived from the one-step and two-step methods

圖6 基于兩步法和單步法的震源參數(shù)（Mw、fc、Δσ、Es）及高頻衰減參數(shù)（κ）估計值的差異Fig.6 Differences in estimated source parameters (Mw, fc, Δσ and Es) and high-frequency decay parameter (κ) by the one-step and two-step methods

非參數(shù)化譜反演的單步法和兩步法給出的遠場高頻路徑衰減差異顯著，高頻震源譜有一定差異，Mw、fc、Δσ、Es及κ的差異在可接受范圍.兩步法將觀測記錄的震源和場地項假定為地震大小相關(guān)的量Ma，包含了該地震所有觀測記錄的平均震源譜和平均場地效應，單個觀測記錄的實際Ma與平均Ma存在偏差，每個距離分段的傳播路徑衰減會吸收該距離分段的所有觀測記錄的Ma偏差（主要來自場地項），如果該距離分段所有觀測記錄的Ma偏差均值明顯偏離0，此時兩步法傳播路徑衰減項出現(xiàn)偏差，推斷這是兩步法和單步法傳播路徑衰減偏差明顯的根本原因.由于高頻路徑衰減差異更顯著，推斷Ma在高頻段的變化性更顯著，可能與高頻場地效應強變化性有關(guān).為避免兩步法路徑衰減被場地項干擾，我們建議盡可能采用單步法，但只研究地震震源特征時，由于震源參數(shù)影響較小，也可用兩步法.

4 震源參數(shù)求解中的權(quán)衡問題

常用網(wǎng)格搜索估計震源譜的特征參數(shù)（Mw、fc、κ），從一系列特征參數(shù)構(gòu)建的理論震源譜中，利用網(wǎng)格搜索尋找最接近反演震源譜的一組特征參數(shù)表示的理論震源譜，作為震源譜特征參數(shù)的最佳估計.估計Mw、fc和κ時，這些參數(shù)之間可能存在的權(quán)衡會影響結(jié)果估計的可靠性.

為評估這些參數(shù)之間的權(quán)衡，本研究利用一系列Mw、fc和κ先驗值，根據(jù)式（4）、（5）構(gòu)建光滑的先驗震源加速度譜（352個），Mw取值 3.75、4.25、5.00 和 6.00，fc取值 0.3、0.5、0.8、1.0、1.5、2.0、3.0、4.0 Hz，κ取值為 0～0.05 s范圍內(nèi)以 0.005 s等間隔遞增的一系列數(shù)值.利用6次典型地震的理論震源譜與反演震源譜的殘差修正先驗震源譜，最終構(gòu)建了2 112個先驗震源譜.網(wǎng)格搜索給出先驗震源譜的Mw、fc和κ估計值，Mw估計值與先驗值之差（m）、κ估計值與先驗值之差（e）隨fc估計值與先驗值之比（l）的變化見圖7.

圖7 震源譜特征參數(shù)之間的權(quán)衡Fig.7 Trade-offs among those parameters describing source spectra

fc先驗值≤ 2.00 Hz時，fc估計值普遍偏小，而Mw估計值普遍偏大，fc估計值偏小的程度越顯著則Mw估計值偏大的程度也越顯著，這說明fc估計值越小則Mw估計值越大，兩者之間存在一定的權(quán)衡.fc先驗值為 0.50～2.00 Hz 時，fc估計值約為先驗值的80%～100%，Mw估計值與先驗值相差小于 0.1，Mw、fc估計值與先驗值的差異均較小，說明Mw和fc之間的權(quán)衡對Mw和fc估計值的影響有限.fc先驗值大于2.00 Hz時，fc估計值偏大而Mw估計值略有偏小，不論fc估計值偏大的程度如何Mw估計值偏小的程度近似不變，說明兩者之間幾乎不存在權(quán)衡.fc先驗值≤ 2.00 Hz 時，fc和κ估計值均偏小，fc估計值越小κ估計值也越小；fc先驗值＞2.00 Hz，fc和κ估計值均偏大，fc估計值越大κ估計值也越大，說明fc與κ之間存在權(quán)衡.由于κ主要控制震源譜的高頻段，為避免fc與κ之間權(quán)衡的影響，確保Mw、fc的可靠估計，應盡量利用中低頻段的震源譜估計Mw和fc.

針對6次典型地震，κ為最佳估計值，Mw與fc的組合構(gòu)建的理論震源譜與反演震源譜的擬合方差見圖8.圖中：紅色實線表示擬合方差最小值對應的fc值.擬合方差的外包絡線均呈開口向上拋物線，fc對應于最小擬合方差值，fc＞0.50 Hz 的 4 次地震，外包絡拋物線在fc最佳估計附近谷值明顯，容易識別；fc≤ 0.50 Hz的 2 次地震，外包絡拋物線在fc最佳估計值附近十分平坦，谷值不明顯.結(jié)果說明，fc＞0.50 Hz時，Mw和fc之間的權(quán)衡對Mw、fc估計影響很小，可不予考慮，fc≤ 0.50 Hz時，需考慮Mw和fc之間的權(quán)衡，一般情況下fc較小地震的震級較大，其矩震級可直接采用震源機制解的結(jié)果以處理Mw和fc之間的權(quán)衡，確保fc的可靠估計.

圖8 擬合方差隨 fc 變化Fig.8 Fitting variance vs.fc

5 結(jié) 論

1）反演結(jié)果與參考場地的選取密切相關(guān)，建議首先綜合考慮臺站觀測記錄數(shù)、場地地質(zhì)資料、H/V譜比曲線等初選參考巖石場地，再以初選作為參考場地，根據(jù)震源譜、震源參數(shù)、高頻衰減參數(shù)試算結(jié)果，最終選取合適的參考場地.

2）單步法和兩步法的遠場高頻路徑衰減差異顯著，高頻震源譜有一定差異，Mw、fc、Δσ、Es及κ差異較小.為避免兩步法路徑衰減被場地項干擾，建議采用單步法；如只研究地震震源特征，亦可采用兩步法.

3）網(wǎng)格搜索估計震源參數(shù)（Mw、fc）和高頻衰減參數(shù)（κ）時，fc和κ之間存在權(quán)衡，為了盡可能避免此影響，建議利用中低頻段的震源譜估計Mw和fc確保其可靠估計.fc＜ 2.00 Hz 時Mw和fc之間存在一定權(quán)衡，但fc＞0.50 Hz 時，此權(quán)衡對Mw、fc估計影響很小，fc較小時地震的震級較大，建議利用Mw處理Mw和fc之間的權(quán)衡，提高fc估計的可靠性.

致謝：中國地震局工程力學研究所基本科研業(yè)務費專項資助（2018B03）.