湯曜瑋，高 亮，鄧遠立，陳 瑛，龔 萱，曾曉燕

（廣東省地震局，廣州 510070）

0 引言

一般來說，固定臺站均選址在地質(zhì)條件較好的基巖之上，可以認為臺站的場地響應(yīng)近似取1，但是也有研究結(jié)果表明，基巖場地并不是完全地?zé)o放大作用。而在臺網(wǎng)的加密過程中，為了臺網(wǎng)布局的需要，有些臺站的場地響應(yīng)也并不是很好，在大震應(yīng)急中，也有很多流動臺站并不是布設(shè)在基巖上，而是在軟土層上。對于松散土層，由于其密度和波速相對較低，造成介質(zhì)的阻抗較低（介質(zhì)的阻抗等于密度和波速的乘積），而地震波的振幅和阻抗的平方根成反比，因此布設(shè)在松散土層上的流動臺相比與理想基巖，其地震計振幅會更大，基于此，在地震研究中，考慮對應(yīng)臺站的場地響應(yīng)，研究場地響應(yīng)計算方法，開發(fā)了與之對應(yīng)的操作簡單、實用化的場地響應(yīng)計算軟件，有助于基層地震工作者用于新建固定臺站、流動臺場地響應(yīng)測算或其他工程地震研究方面，具有一定的應(yīng)用價值。

Matlab 是由MathWorks 公司推出的一款強大的數(shù)值計算軟件，是數(shù)值計算、符號運算和圖形處理等多種功能的強有力實現(xiàn)工具[3]。近年來Matlab這一強大的科學(xué)計算軟件包已經(jīng)得到業(yè)界的廣泛認可，并已深入到了各行各業(yè)的眾多學(xué)科、在各大公司、科研機構(gòu)、大學(xué)校園得到了日益普及與廣泛應(yīng)用。其具有以下優(yōu)點：①代碼簡捷直觀，以矩陣為運算基礎(chǔ)，在計算方面避免使用者花費大量時間去編寫不必要的小程序；②數(shù)值計算方面，Matlab基本涵蓋了所有的數(shù)值計算方法，從基本的求極限、解微分方程，數(shù)字信號處理FFT，到有限差分、線性擬合，再到最新的最優(yōu)化技術(shù)諸如神經(jīng)網(wǎng)絡(luò)、遺傳算法等等，使復(fù)雜的數(shù)值計算變得容易；③數(shù)據(jù)可視化方面，Matlab提供了多種基本圖件處理函數(shù)，極大的方便繪制出高質(zhì)量的工程圖譜；④Matlab提供了強大的編程接口，支持Matlab與其他編程語言進行數(shù)據(jù)交換。自地震觀測從模擬觀測步入數(shù)字化觀測以來，地震儀器拾取的地面震動信號本質(zhì)上也是一種隨機數(shù)字信號，而Matlab編程平臺自帶的數(shù)字信號處理工具箱對于數(shù)字信號的各類型常規(guī)處理諸如卷積、濾波、加窗、FFT、相關(guān)性分析等內(nèi)置了多種函數(shù)[4]，為本軟件的研制也提供極大的編程方便，節(jié)省了開發(fā)時間。

1 譜比法場地響應(yīng)計算原理及數(shù)據(jù)處理方法

1.1 譜比法場地響應(yīng)計算原理

設(shè)VS為地表處地震動垂直分量的地震波振幅譜，VB為基底處地震動垂直分量的地震波振幅譜，HS為地表處水平分量的地震波振幅，HB為基底處的水平分量的地震波振幅譜，則經(jīng)驗傳遞函數(shù)[5]為

即井下（基巖處）地震動的水平分量和垂直分量大致相等，垂直分量經(jīng)過軟土層的放大效應(yīng)遠遠小于水平分量的放大效應(yīng)，基本沒有被放大，在此基礎(chǔ)上，經(jīng)驗傳遞函數(shù)可以簡化為[6]譜比率S=SN，于是有

根據(jù)上式（6）用地表記錄到的水平向地震波振幅譜和垂直向地震波振幅譜比就可以計算得到各臺站的場地響應(yīng)，根據(jù)以往文獻，除了利用地震事件計算場地響應(yīng)，同樣也可以利用地震計記錄到的地脈動來計算場地響應(yīng)。在本軟件中，可以對地震事件和地脈動分別進行計算求取對應(yīng)臺站場地響應(yīng)。

1.2 數(shù)據(jù)處理方法

當采用臺站記錄到的地震事件來計算場地響應(yīng)時，在對地震記錄數(shù)據(jù)進行譜分析時，由于不同震中距的臺站記錄到的同一地震事件，其地震記錄S 波持續(xù)時間并不相同，如果對記錄S 波直接進行FFT 轉(zhuǎn)換求取幅度譜，其頻率譜間隔疏密必定不同，為了得到相同間隔的頻譜，在軟件中，我們采用經(jīng)典平移窗譜[7]的方法。即把所取的S 波窗分割成若干個包含256個數(shù)據(jù)點的小段，并使相鄰數(shù)據(jù)段有50%的重疊間隔，以獲得穩(wěn)定的波形記錄的傅里葉譜，同時為了防止切割波形記錄產(chǎn)生頻譜泄露，在每一小段波形的起始和末尾加以5%的COS 邊瓣。對于采樣率為100 Hz 的地震記錄來說，每個小段的時間長度是2.56 s，通過快速傅里葉變換得到每個小段的傅里葉譜。然后對每個小段的傅里葉譜進行儀器響應(yīng)校正，這樣通過下式[8]

我們可以得到整個“S 窗”內(nèi)觀測信號的位移譜振幅，式中v2i( f )是經(jīng)過校正掉儀器響應(yīng)之后的S波波形記錄分割后第i小段的傅里葉快速變換譜，T 為“S 窗”的持續(xù)時間，該S 窗包含了n 個時間長度為t 的256 個采樣點的小段。最后然后通過內(nèi)插得到對數(shù)頻率為0.0， 0.05， 0.1， 0.15， 0.2，…，1.3 的共計27 個頻率點的位移譜振幅值。

在以上進行譜分析之前，由于地震臺站地震記錄中均可能出現(xiàn)直流偏移，這些偏移均可能包含在地震事件記錄中或者使地脈動波形記錄不是真正的地脈動。在實際計算中均需要進行直流分量剔除。

角礫狀構(gòu)造：早期形成的黃鐵礦被錯動成棱角分明，大小不等的角礫分布在石英脈中，有的角礫被次生鐵錳礦物膠結(jié)。

2 軟件設(shè)計思路及部分實現(xiàn)方法

根據(jù)上述譜比法場地響應(yīng)計算原理以及平移窗譜法數(shù)據(jù)處理方法，軟件總體設(shè)計思路可分為以下幾個步驟:①導(dǎo)入原始地震事件波形數(shù)據(jù)或地脈動波形；②導(dǎo)入對應(yīng)地震臺站儀器參數(shù)數(shù)據(jù)；③導(dǎo)入地震事件震相數(shù)據(jù)，根據(jù)P 波與S 波走時差經(jīng)驗估算S波持續(xù)長度，對于地脈動記錄則設(shè)置統(tǒng)一時間長度；④計算某批次地震各臺三分量位移振幅譜；⑤合成水平向振幅譜比垂直向得到某批次記錄各臺場地響應(yīng)值；⑥對多批次地震記錄或地脈記錄求得場地響應(yīng)求均值。

2.1 數(shù)據(jù)導(dǎo)入

軟件導(dǎo)入的地震記錄數(shù)據(jù)為JOPNES 地震速報分析系統(tǒng)截取的導(dǎo)出的sac 數(shù)據(jù)文件，sac 數(shù)據(jù)格式廣泛用于地震學(xué)分析軟件中，sac 格式也是地震波形數(shù)據(jù)的通用格式。sac 數(shù)據(jù)格式有ASCII 碼和二進制兩種形式，本軟件數(shù)據(jù)導(dǎo)入采用sac 二進制碼格式。軟件實現(xiàn)了對多臺地震事件記錄的批量讀取。在批量讀取地震事件記錄的同時，也批量讀入對應(yīng)地震臺站地震計參數(shù)信息，其地震計參數(shù)信息是JOPNES 地震速報分析系統(tǒng)截取的導(dǎo)出的ascii 文件中的par 文件，其中包含了對應(yīng)地震記錄的地震計靈敏度、傳遞函數(shù)零極點等信息。另外，對于地震事件，為了準確確定地震記錄S波流逝事件，軟件也直接批量讀入JOPNES 地震速報分析系統(tǒng)產(chǎn)出的震相phase文件。

2.2 地震S波窗長確定

精確的S波持續(xù)窗長是地震記錄位移譜分析的基礎(chǔ)，在軟件中導(dǎo)入的phase 文件均是對應(yīng)地震記錄不同臺站各類波到時的編目結(jié)果，精度可靠。在導(dǎo)入phase 文件后，對于S 波持續(xù)時間長度的確定，根據(jù)劉杰等指出用于計算位移譜的S波時間長度應(yīng)為為從S 波開始到包含90%S 波能量的時間段[9]，通過回歸分析1900 條地震記錄得到S 波截取時間長度與Pg、Sg 波走時差為線性關(guān)系，并擬合出具體線性表達式，如下式所示：

其中Ttotal是S波截取持續(xù)時間，TSg是S波到達時間，TPg是P波到達時間，在軟件圖形化波形顯示區(qū)域，根據(jù)導(dǎo)入的phase文件編目震相和走時經(jīng)驗公式其S波窗長是可靠的。

2.3 場地響應(yīng)值求取

對所計算臺站根據(jù)其每次地震記錄三分向利用公式（7）計算地震波位移振幅譜，然后將水平EW 與NS兩分向利用公式（9）做合成，垂直向VS（f）即為地震計垂直向位移振幅譜。

再將（9）與（10）求比值，即得該次地震事件不同頻點的譜比法場地響應(yīng)值。對于某個臺站，可以計算多次不同地震事件的場地響應(yīng)值，最后將多批次地震事件場地響應(yīng)統(tǒng)計求均值即為該臺站場地響應(yīng)。

如果采用地脈動記錄來計算場地響應(yīng)，無需考慮地震事件中S波持續(xù)事件問題，我們將平靜時段的地脈動數(shù)據(jù)統(tǒng)一按600 s 一段截取[10]，也是每段數(shù)據(jù)疊加前一段50%記錄，分別FFT 計算每段速度譜并轉(zhuǎn)換位移振幅譜后，每段水平向譜值比垂直向譜得到每段場地響應(yīng)，最終求得該臺站場地響應(yīng)均值。

2.4 壞譜結(jié)果剔除

在計算過程中發(fā)現(xiàn)，對于批量處理地震事件時候，部分臺站由于漂移、地震記錄中斷等原因，計算得到的譜值明顯異常，為了避免異常譜值影響最終多批次場地響應(yīng)平均計算結(jié)果，在軟件中設(shè)計了譜值存儲前剔除功能。

3 軟件使用的簡單說明

3.1 軟件基本界面

軟件主界面如下圖1所示，界面分為包括以下幾個部分，參數(shù)輸入及操作流程控制部分，波形顯示部分，位移譜值結(jié)果及場地響應(yīng)結(jié)果顯示三部分，面板左邊主要是參數(shù)輸入及計算流程控制部分，中間三個坐標軸是波形記錄三分量顯示部分，左邊兩個坐標軸用于顯示位移譜值結(jié)果與場地響應(yīng)結(jié)果。參數(shù)輸入部分又分別按照地震記錄或地脈動記錄兩部不同計算方式進行分類輸入。在左邊的操作流程控制框里，依次導(dǎo)入各項參數(shù)數(shù)據(jù)后，點擊生成位移振幅譜即可在軟件后臺自動完成位移譜計算全過程，然后通過下拉框可以選擇對應(yīng)臺站，在中間及右邊觀察對應(yīng)臺站三分向波形記錄或地脈動記錄、位移譜和場地響應(yīng)結(jié)果、根據(jù)觀察情況選擇是否保存該次譜結(jié)果，在多批次的結(jié)果計算完成之后，可以存檔保存，最后點擊第三個面板框里的求取場地響應(yīng)均值完成最終計算并保存結(jié)果圖件。

3.2 軟件操作流程

分別以地震事件和地脈動記錄為數(shù)據(jù)源介紹軟件使用流程，軟件的基本功能架構(gòu)圖如下圖2所示，首先導(dǎo)入各類基礎(chǔ)數(shù)據(jù)，后臺計算，結(jié)果顯示，結(jié)果重處理。

圖1 程序界面Fig.1 Program interface

圖2 基本功能架構(gòu)圖Fig.2 the basic function diagram

當處理地震事件記錄時，點擊左側(cè)上方面板框內(nèi)批量導(dǎo)入地震事件、批量導(dǎo)入震相文件、批量導(dǎo)入儀器參數(shù)完成后，再點擊批量生成場地響應(yīng)，即可在軟件完成該批次地震所有臺站的場地響應(yīng)計算，場地響應(yīng)包含的如S波窗長確定、FFT譜計算、譜比法計算過程都在軟件內(nèi)部后臺自動計算，計算完成后，在右邊上部坐標軸顯示出該次地震所有臺站的水平與垂直位移譜線，同時在按臺站顯示的下拉菜單中會自動導(dǎo)入該次計算中的所有臺站，然后點選某各臺站，可以在右中部顯示坐標軸中顯示該臺站此次地震波形、水平及垂直向位移譜結(jié)果，如下圖3所示。

中間三坐標軸顯示德慶臺該次地震記錄，并用三條紅線顯示了PG、SG到時與SG波截斷時長。右側(cè)下部坐標軸為該次地震記錄位移譜。同時軟件會將該次地震所有計算臺站位移譜結(jié)果保存至硬盤某文件夾，文件保存內(nèi)容如下圖4 所示，含有該次地震相關(guān)信息，各臺的不同頻點位移譜值通過軟件的顯示界面，我們可以觀察求得的位移譜，如果發(fā)現(xiàn)有某臺的譜值異常，我們可以選擇剔除壞譜，不予保存至文件里。

當計算多批次地震之后，多批次的地震位移譜結(jié)果會依次保存下來，然后點選軟件面板左部最下方多批次結(jié)果求平均的控制框，讀入保存文件里的震幅譜結(jié)果，再選擇臺站既可以得到該臺多次記錄結(jié)果的場地相應(yīng)的平均值，如下圖5 所示，圖四上方坐標軸為兩次地震記錄的震幅譜，其中粗線為其震幅譜均值，下方坐標軸為其譜求平均之后，然后進行頻比求得場地相應(yīng)，其中綠線為最終場地響應(yīng)。地脈動記錄軟件操作過程與之類似，在此不在贅述。

圖3 德慶臺地震波形記錄、該次地震位移譜及場地響應(yīng)Fig.3 Seismic waveform record，displacement spectrum and site response of the earthquake recorded by Deqing station

圖4 地震記錄位移譜結(jié)果文件Fig.4 Displacement spectrum result file of seismic records

圖5 震幅位移譜均值及場地響應(yīng)結(jié)果Fig.5 The mean value of the seismic amplitude displacement spectrum and site response results

圖6 部分臺站場地響應(yīng)結(jié)果Fig.6 Site response results of some stations

4 軟件使用實例

選取廣東臺網(wǎng)2017—2018年間7次發(fā)生在省內(nèi)速報地震，所有地震震級均在ML2.5 級以上，廣東臺網(wǎng)大部分臺站以及部分省外臺站均有清晰的記錄，按照以上介紹的軟件使用流程，依次批量導(dǎo)入這7 次地震的多臺地震記錄、震相記錄及儀器參數(shù)，最后計算得到廣東大部分測震國家固定臺站的場地響應(yīng)，現(xiàn)選擇部分臺站結(jié)果顯示如下圖6所示:從結(jié)果看出大部分臺站場地響應(yīng)結(jié)果在100上下左右浮動，比較符合一般性實際臺站場地響應(yīng)結(jié)果預(yù)期。說明軟件計算場地響應(yīng)達到預(yù)期效果。

5 結(jié)論與探討

本文以Matlab為開發(fā)平臺，編制了基于譜比法臺站場地響應(yīng)計算軟件，該軟件界面友好，批量導(dǎo)入地震記錄、儀器參數(shù)以及震相記錄，在導(dǎo)入基礎(chǔ)數(shù)據(jù)后，在軟件后臺計算各臺位移譜，并用文件形式保存結(jié)果，最后將某臺多次位移譜結(jié)果求平均后計算該臺場地響應(yīng)，整個計算過程在軟件中集成一體，對于地震工作者來說，該軟件界面友好，操作簡單，可用于新建固定臺站、流動臺場地響應(yīng)測算或其他工程地震研究方面、具有一定的應(yīng)用價值。