基于MATLAB的場地效應及路徑衰減聯(lián)合反演軟件研制

1、PAGE PAGE 15場地響應及路徑衰減聯(lián)合反演方法軟件實現(xiàn)及應用朱新運 摘要：地震波衰減及場地響應研究是地震學研究領域重要的前沿課題。基于譜比法，聯(lián)合反演求解場地響應及地震波衰減參數(shù),這種算法理論上是可行的，據(jù)此，基于MATLAB開發(fā)平臺開發(fā)研制了地震波衰減參數(shù)與臺站場地響應聯(lián)合反演軟件。本文介紹軟件設計的思路、理論原理，方法實現(xiàn)及資料處理過程，簡述軟件模塊功能及使用方法。該軟件計算過程穩(wěn)定、可靠，界面設計友好，可用于日常地震資料的分析，也可作為地震科學研究的工具。關鍵詞 場地響應 衰減參數(shù) 聯(lián)合反演 軟件 1 引言在對耦合了路徑、場地、儀器、地震輻射的方向性等多種響應成分的地震記錄（陳運

2、泰等，2000）進行必要限定的情況下計算波譜比，可以扣除耦合的幾何擴散、儀器響應、地震輻射方向性，但是場地響應與非均勻性、非彈性路徑衰減則難以簡單地扣除的。為了獲得較為準確的衰減系數(shù)，Chun等(1987)提出的“逆向雙臺譜比法”。該方法是計算Lg波衰減參數(shù)比較可靠的方法（Chung,2002），它通過兩臺站連線兩端延長線位置上發(fā)生的成對地震的記錄譜比計算衰減參數(shù)。由于這種方法在資料選擇上要求比較苛刻，影響了該方法的推廣應用。Chung等（2005）對“逆向雙臺譜比法”進行了改進，改進后的“逆向雙臺譜比法”不再嚴格要求地震對以及臺站對處在一條直線上，只要求地震對及臺站對之間的距離滿足一定的條件

3、即可。這一改進擴大了地震資料的選擇范圍，但對資料的應用仍然受到一定的限制。場地響應描述的是近地表地層介質對地震地面運動的放大作用，研究臺站場地響應（Bonilla et al, 1997; Riepl et al, 1998;Langston, 1979; Nath et al, 2000,2001）的方法很多，但基本上是通過路徑衰減項簡化求得場地響應的。根據(jù)Shearer（1999）的研究結果，場地響應值與近地表地層介質的介質密度及地震波速度有關，與地震波傳播距離無關，而地震衰減參數(shù)是地震波路徑衰減的描述,僅與距離有關，由于地震Lg波的衰減對地震輻射模式不敏感(Mitra et al，200

4、6)，因而使用Lg波作為基本資料進行相關研究時，臺站方位可以忽略。據(jù)此，朱新運等（朱新運等，2007，ZHU Xin_Yun,et al,2007;2012;2013;2014）以地震Lg波臺站記錄譜為觀測數(shù)據(jù)，在資料選擇上，使用Ms5級以上地震，為保證通過譜比消去震源因子時不致引入較大誤差，限定臺站對地震方位角變化不大，同時為保證衰減參數(shù)穩(wěn)定，限定地震波穿過構造相對單一的區(qū)域，利用路徑衰減與臺站震中距有關、而場地響應與震中距無關這一特點，通過聯(lián)合反演將場地響應與路徑衰減互相分開，得到了較好的結果。本文中作者對場地響應與路徑衰減參數(shù)聯(lián)合反演方法論證，同時介紹了在MATLAB平臺下研制的路徑衰減

5、與場地響應聯(lián)合反演軟件，及使用該軟件研究了華北、華東等區(qū)域的地震波衰減及場地響應參數(shù)成果。軟件設計構件合理、界面美觀、計算功能完善、使用方便，可以作為地震波相關研究的實用工具。2 理論與方法臺站記錄的地震波的頻譜可以表示為（Hasegawa,1985）： （1）式中，表示頻率；表示震源因子項；表示震源輻射圖型因子；表示震中距；表示儀器效應；表示場地響應；表示幾何擴散因子， （2）是表示幾何擴散特性的參數(shù)，對于面波波，體波；表示介質衰減項，它可以表示為： （3）式中，是衰減系數(shù)， （4）在式（4）中，表示波速度；。假定參與運算的地震數(shù)為，需要評估的臺站數(shù)為，兩個臺站（第臺與第臺，）對同一個地震（

6、第個地震）的記錄譜的比可表示為： （5）在上式中，對于同一地震，震源因子項之比可以消去；震源輻射圖型因子與震源機制解有關，對較大地震選取資料時可考慮震中距離與臺站間距離相比很大，以至于震源對臺站之間方位角變化很小，則：，對較小地震則可看作點源，震源輻射因子影響不大，從而，在此情況下震源輻射圖型因子項之比可以消去；儀器效應與是已知的；是臺站記錄譜；由式（2）給出。將式（3）代入式（5）整理后得： （6） 對式（6）兩邊取自然對數(shù)可以得： （7）式中， （8） (9) （10）令： （11） （12）則： （13）選取資料時考慮研究時間段內，在研究區(qū)即地震波射線穿過的區(qū)域內，沒有發(fā)生大的構造活動，

7、衰減參數(shù)保持恒定；臺站場地響應僅與頻率有關，兩臺站對地震效應的幅度的比值僅與頻率有關；在參與評估的臺站中，對應任何頻率總有一個臺站的效應是理想的，即：對應該頻率該臺站效應為理想基巖場地響應，效應幅度為一個單位。首先以第一個臺為參考臺站，則： （14）式中，各角標意義與前述各式意義相同，寫成矩陣形式，則為：(14)(15)在式（15）中， 參與評估臺站數(shù)，分析頻率個，頻率代碼（相關角標），則對應的有個未知數(shù)，不同頻率下各臺站場地響應對參考臺站場地響應比值有，共個未知數(shù)，方程組的個數(shù)為。解方程（15）可得相應未知數(shù)，其中對應衰減部分為，利用（4）式，由可得到值；其余部分為場地響應比值數(shù)據(jù)，參考臺站

8、對自身在不同頻率下的比值均為1，對不同分析頻率在各臺站場地響應比值數(shù)據(jù)中逐一搜索比值最小值，以該最小值為相應頻率在理想基巖場地響應的幅度單位，把該頻率其余臺站對應結果與相應最小值對比，得到該頻率在相應臺站的最終效應幅度，進而得到所有臺站對所有分析頻率的場地響應。3 軟件實現(xiàn)軟件設計包括文件操作、參數(shù)設置及獲取、各個計算過程實現(xiàn)等模塊,程序設計依據(jù)MATLAB平臺特征，數(shù)據(jù)運算盡可能矩陣化，共有代碼盡可能函數(shù)化。主窗口及窗口控件初始化直接使用GUI生成，程序控件一般屬性，包括控件類別、控件在主窗口中的位置等通過GUI模板設計后進行相應賦值，主界面見圖1。 圖1 程序主界面（群速度高端太小，沒有截

9、取到合理數(shù)據(jù)段、終止計算）3.1 參數(shù)設置軟件對各關鍵控制參數(shù)均提供了數(shù)據(jù)接口，其中高、低頻率及頻率間隔數(shù)限定軟件計算從低頻到高頻以等間隔分隔的各頻率結果，考慮到在付氏變換過程中記錄譜的突跳等影響數(shù)據(jù)穩(wěn)定，軟件設計平滑功能，在平滑窗鍵入相關參數(shù)，在數(shù)據(jù)計算過程中進行數(shù)據(jù)平滑。群速度參數(shù)用于控制地震波有效成分截取，該參數(shù)包括低端速度和高端速度，軟件根據(jù)該參數(shù)截取對應時間窗的地震波形，為了有效截取波形，在準備地震文件時應提供較為準確的地震震源位置。幾何擴散參數(shù)用于對波形進行幾何擴散校正，一般面波校正參數(shù)為0.5，體波校正參數(shù)為1。截取有效波形后加窗進行付氏變換，軟件提供漢寧窗、漢明窗等六種窗函數(shù)及

10、窗寬大小編輯窗，可設計計算過程中使用的具有相關屬性的窗函數(shù)。3.2 結果顯示及處理結果處理區(qū)包括了圖像查看、結果存儲等控件，圖像查看包括過程圖像查看和結果圖像查看，過程圖像查看相當于計算過程的回放；通過查看關鍵步驟產生數(shù)據(jù)圖像及結果圖像，可以對整個計算過程及計算結果是否合理有一個全面了解。在“波、譜圖”區(qū)域，用戶可通過鍵入地震記錄條數(shù)并選擇“波形圖”或“波譜圖”查看本條地震記錄波形及在速度窗控制下波形被截取或波譜計算情況，在“衰減顯示”下方選擇衰減分向查看該分向衰減圖形（圖2）或選擇“三向綜合”查看所有三個分向計算的衰減結果，在“綜合場地響應”下方選擇分向名，繪圖區(qū)繪制在同一坐標系下各臺場地響

11、應，在“單臺場地響應”下選擇臺站代碼，繪圖區(qū)繪制該臺站場地響應結果。結果存儲包括數(shù)據(jù)結果和圖像結果存儲兩部分，“當前圖像存儲”按鈕完成當前繪圖區(qū)圖像存儲，為便于把多臺場地響應結果繪制在同一頁面，該軟件專門設計了“多臺場地響應”存儲區(qū)，鍵入單頁面縱、橫參數(shù)，點擊“存儲”按鈕完成多臺站場地響應圖件結果存儲。圖2、軟件顯示圖形實例（1）.N_S向衰減結果4 算例：首都圈區(qū)域Lg波衰減及場地響應研究4.1 資料條件首都圈數(shù)字地震臺網(wǎng)由107個地震臺站組成，均為3分量數(shù)字地震儀，采樣率為50Hz。其中44個寬頻帶臺站的儀器為在20s20 Hz速度平坦的FBS一3A地震儀；58個短周期臺站的儀器為在120

12、 Hz速度平坦的FSS-3地震儀，其余5個為在120s20 Hz之間平坦的CTS一1甚寬頻帶地震儀。甚寬帶和寬頻帶地震儀大多分布在山區(qū)的地表臺站，短周期地震儀大多分布在平原地區(qū)，5臺安裝在地表，其余均為安裝在第四紀黃土層中的井下。上述儀器的性能如圖3所示。圖3 首都圈數(shù)字地震臺網(wǎng)地震儀器性能首都圈地區(qū)地質構造復雜，尤其是山區(qū)和平原區(qū)地殼厚度、地質構造差異較大，東南為凹陷盆地，而西北部為山區(qū)?？紤]到首都圈地區(qū)地質構造的復雜性，我們根據(jù)地震射線傳播路徑將評估衰減參數(shù)的資料分成相互獨立的四組進行評估；第一組：西北部山區(qū)臺站對同區(qū)域地震的記錄；第二組：東南部凹陷盆地臺站對同區(qū)域地震記錄；第三組：西北部

13、山區(qū)臺站對東南部凹陷盆地地震的記錄；第四組：東南部凹陷盆地臺站對西北部山區(qū)地震的記錄。其中，第一、第二兩組的地震射線分別在各自相同的區(qū)域傳播，第三、四組均沿相反路徑穿過地殼陡變區(qū)域。在2006年01月至2008年底記錄的404次震級ML1.75.4的地震中，作者挑選了信噪比高的1000多條波形，涉及的地震臺達81個，地震射線復蓋了首都圈西北山區(qū)和東南盆地，地震、地震臺及地震的射線分布情況如圖2所示。圖4 地震臺、地震射線及地震分布圖。 (a)隆起區(qū)臺站對同區(qū)域地震的記錄射線（粗虛線）東南部凹陷盆地臺站對同區(qū)域地震記錄（細實線）;(b)西北部山區(qū)臺站對東南部凹陷盆地區(qū)域地震的記錄射線（粗虛線），

14、東南部凹陷盆地臺站對西北部山區(qū)地震記錄（細實線）4.2 資料處理在本研究中，Lg波段的截取采用了固定速度窗法。利用清晰的地震波震相到時求得g波的速度上界3.65km/s，取g波速度下界為2.6km/s (Murphy et al,1982; Cappillo et al, 1985)，以該速度窗截取波形（圖5），對截取的g波形段補零擴展至最近的2的整數(shù)次方倍，為防止譜泄漏，兩端加2的余弦窗進行富利葉變換、然后做儀器校正、幾何衰減校正，并以3個單位頻率步長進行平滑（圖5），使用數(shù)據(jù)段信號譜大于噪聲譜（圖6）。根據(jù)式（14）進行擬合（圖7）。圖5 Lg波段的截取（a）及波譜圖（b）實例圖6 信號與

15、噪聲譜比較圖7 數(shù)據(jù)擬合實例4.3 數(shù)據(jù)結果4.3.1 衰減系數(shù)從同一組地震射線中三個分向計算結果可以發(fā)現(xiàn)在6Hz以上總趨勢發(fā)生拐彎現(xiàn)象，且在三個分向出現(xiàn)明顯差別（圖8）。這個結果與Shin等(1987)及Ral等（2008）的研究結果一致，Shin等(1987)在研究Lg波衰減時發(fā)現(xiàn)在7Hz以上，Lg波的衰減發(fā)生了較大變化，并采取了一定的辦法進行校正，可是并沒有獲得理想結果，Shin等(1987)對地震波衰減不一致的情況做了解釋，認為Lg波的高頻部分受到了Sn及Pn波的尾波干擾，因此建議研究Lg波衰減可在7Hz以下或者挑選震中距小于600km的地震的記錄。Ral等（2008）在對墨西哥東北索

16、諾拉州地震波衰減研究中也出現(xiàn)6或7Hz以上衰減不一致的情況，但沒有進一步解釋。為此，在本研究中我們僅計算16Hz之間衰減結果并擬合相關關系。圖8 地震波1Hz至10Hz的衰減關系第一組射線包括西北部山區(qū)臺站對同區(qū)域地震的86條記錄圖（No.1），第二組射線包括東南部凹陷盆地臺站對同區(qū)域地震的1036條記錄圖（No.2），第三組射線包括西北部山區(qū)臺站對東南部凹陷盆地地震的記錄182條（No.3），第四組射線包括東南部凹陷盆地臺站對西北部山區(qū)地震記錄196條（No.4）。四組衰減關系中第一組最大，第四組最小，其余兩組居中（圖9僅畫出垂直向）。圖9 地震波衰減關系(No.1、No.2、No.3、No

17、.4分別代表第一、第二、第三、第四組數(shù)據(jù))4.3.2 場地響應根據(jù)資料情況，首都圈區(qū)域東北地震資料較少，無法滿足數(shù)據(jù)運算條件。在西北山區(qū)及東南凹陷盆地區(qū)域部分臺站前后儀器響應不一致甚至差別較大，有些在資料選擇時間段沒有足夠地震記錄或有較多地震記錄但在選擇Lg波段時資料被剔除等原因，實際上共有81個臺站得到穩(wěn)定的場地響應。在參與評估的臺站中，山區(qū)臺站基底均為硬基底類地表臺；平原區(qū)選擇的臺站中有8個臺站為地表基巖臺，l0個臺站為井下泥巖臺，其余臺站均位于井下黃土沉積層中?；鶐r的場地響應多數(shù)在16 Hz頻帶內比較平坦，部分在頻率高端有抬升（56Hz）。在低頻段（12Hz），黃土沉積的場地響應普遍高于

18、硬基底的基巖場地（圖10）。場地響應描述的是近地表地層介質對地震地面運動的放大作用，是衡量臺站附近地下介質對地震波的阻尼影響程度的物理量。根據(jù)Shearer（Shearer，1999）有關的研究結果，松軟的沉積場地對地震動的放大作用較大；堅硬的基巖場地則反之，本研究結果與此結論一致。圖10 場地響應圖。實線表示垂直向場地響應數(shù)據(jù)結果，點線代表南北向結果，虛線表示東西向結果5 討論與結論使用該軟件除了華北區(qū)域衰減及臺站參數(shù)進行研究之外，也處理了華東區(qū)域地震波衰減及場地響應參數(shù)，并對這些結果與其他研究者在全球不同區(qū)域使用不同方法所獲得的區(qū)域Lg波衰減參數(shù)，筆者對這些搜集并與之進行了比較（圖11）。

19、本研究數(shù)據(jù)結果較好的反映了研究該區(qū)域地震活動水平。圖11 不同區(qū)域Lg波衰減參數(shù)比較本研究以譜比法為基礎，利用路徑衰減與臺站震中距有關、而場地響應與震中距無關這一特點，通過數(shù)據(jù)計算實現(xiàn)衰減參數(shù)項與場地響應項分離。軟件設計以MATLAB為平臺開發(fā)，設計上充分考慮開發(fā)平臺特征，發(fā)揮其優(yōu)勢，以節(jié)省內存、運行速度快、流程可控、軟件可讀性強及操作簡便為基本要求，功能上集數(shù)據(jù)預處理、參數(shù)計算、結果處理于一體，應用可視化操作，界面友好、功能完善，數(shù)據(jù)輸入、輸出及參數(shù)配置操作方便實用。參考文獻陳運泰，吳忠良，王培德等.數(shù)字地震學北京：地震出版社：114朱新運，于俊誼，2008，基于MATLAB的小震震源參數(shù)計

20、算軟件研制，2008，西北地震學報，30（6）：380384.朱新運，陳運泰，用Lg波資料反演場地響應與地震波衰減參數(shù)，2007，地震學報，29（6）：569580.朱新運，2013，數(shù)字地震波分析系統(tǒng)，浙江大學出版社。朱新運,2016,華北盆地Lg 波衰減及場地響應特征, 地球科學， 41(12):2109-2117Borcherdt R D. 1970. Effects of local geology on ground motion near San Francisco Bay. Bull Seism Soc Amer, 60:2961.Brune J N.1970.Tectonic

21、stress and the spectra of seismic shear waves from earthquakeJ.J Geophys Res,75(26)：4 9975 009Brune J N.1971.CorrectionFJ.J Geophys Res.76：5002Chun K-K, Gogdon F W, Richard J K, et al.1987. A novel technique for measuring Lg attenuation-results from eastern Canada between 1 to10. Bull Seism Soc Amer

22、,77:398419.Garccf D, Singh S K, Herriz M, et al.,2004.Inslab Earthquakes of Central Mexico: Source spectra, and stress drop.Bull Seism Soc Amer.,94:782802.Hartzell S H. 1992. Site response estimation from earthquake data. Bull Seism Soc Amer, 82:23082327.Hasegawa H S. 1985. Attenuation of Lg waves i

23、n the Canada Shield. Bull Seism Soc Amer, 75:15691582.Herrmann R B, Kijko A.1983.Modeling some empirical vertical component Lg relation, DARPA/AFGL Seismic Research Symposium,6-8 May,1985,U.S. Air Force Academy, Colorado Springs, ColoradoKing J L, Tucker B E. 1984, Observed variations of earthquakes

24、 motion across a sediment-filled valley. Bull Seism Soc Amer., 74:137152.Langston C A. 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J Geop-hys Res. 84:47494762.Moya C A，Aguirre J，Irikura K2000Inversion of source parameters and site effects from strong ground

25、 motion records using genetic algorithms-JBull Seism Soc Amer,90(4)：977992Nath S K, Biswas N N, Dravinski M, et al. 2001. Determination of S-wave site response in Anchorage, Alaska in the 19 Hz frequency band, Pure Appl Geophys. 159:26732698.Nath S K, Sengupta P, Sengupta S, et al. 2000. Site response estimation using strong motion network