基于波譜分析方法的波速比計(jì)算軟件

高 亮 湯曜煒 楊 選

（中國(guó)廣州510070 廣東省地震局）

0 引言

傳統(tǒng)測(cè)定波速比的方法為和達(dá)法，有單臺(tái)和達(dá)法和多臺(tái)和達(dá)法。實(shí)際中較多采用多臺(tái)和達(dá)法。多臺(tái)和達(dá)法利用一組臺(tái)站記錄的地震事件的P 波走時(shí)tP和S 波與P 波走時(shí)差tS-P來測(cè)定波速比，該方法存在以下問題：①受時(shí)間精度的影響；②震相識(shí)別存在誤差；③得到的波速比是各臺(tái)站至震源傳播路徑的平均值，并非震源處的波速比；④受臺(tái)站布局的影響，臺(tái)站分布不均勻時(shí)得到的只是局部區(qū)域的波速比。

在對(duì)中強(qiáng)地震前后地震序列的研究中，地震學(xué)家發(fā)現(xiàn)波速比往往會(huì)呈現(xiàn)出“下降—恢復(fù)”的異常變化。孕震理論對(duì)以上現(xiàn)象解釋為，在地震孕育過程中，特別是在其后期，孕震區(qū)內(nèi)巖石的大量微裂隙會(huì)發(fā)生空間排列的變化，介質(zhì)的孔隙度和各向異性程度等隨之發(fā)生改變，以致地震波的傳播速度受到影響。由于和達(dá)法求取波速比存在缺陷，部分學(xué)者傾向利用地震波波譜分析求取震源處的波速比。這種方法的優(yōu)點(diǎn)是不受時(shí)間精度和震相識(shí)別誤差的影響。在地震預(yù)報(bào)實(shí)踐中，二者之間可以彼此印證，互相補(bǔ)充。目前國(guó)內(nèi)外有不少學(xué)者進(jìn)行過相關(guān)研究，均取得了較好的效果。許健生等（1998）對(duì)天祝—古浪5.4 級(jí)地震的研究表明，主震后地震的縱、橫波位移譜零頻極限的對(duì)數(shù)比值從1.2 左右下降到1.0 左右；張永久等（2004）通過縱、橫波位移譜零頻極限比值的計(jì)算，得到雅江6.0 級(jí)地震序列前后波速比的異常變化；陳麗娟等（2015）計(jì)算了縱橫波位移譜零頻極限比值與拐角頻率比值，得到蘆山MS7.0 地震序列前后波速比異常變化。利用拐角頻率比與零頻極限比值計(jì)算波速比的過程較為復(fù)雜，涉及地震儀器響應(yīng)計(jì)算、縱橫地震波譜分析、縱橫波拐角頻率與零頻極限計(jì)算等。當(dāng)前，研發(fā)一款使用方便、界面友好的軟件以進(jìn)行波譜分析下的波速比計(jì)算是有意義的。筆者在MATLAB 平臺(tái)下，開發(fā)了一款利用縱橫波拐角頻率與位移譜零頻極限值求解波速比的軟件。本軟件用于計(jì)算震源處的波速比，集成了上述計(jì)算功能，且具備繪圖功能，在操作界面設(shè)有對(duì)應(yīng)控件，界面友好，操作簡(jiǎn)便。

1 波速比求解理論

1.1 地震波形譜分析

要計(jì)算拐角頻率、位移譜零頻極限值，波譜分析是計(jì)算基礎(chǔ)。波譜分析在軟件中采用平移窗譜方法（Chael，1987）。首先，將選取的P、S波波形分成若干包含256個(gè)采樣點(diǎn)的小段，并使相鄰小段之間有50%的重疊。對(duì)于采樣率為50 Hz 的地震記錄來說，每個(gè)小段的時(shí)間長(zhǎng)度是 5.1 s。然后，在每一小段波形的起始和末尾加5% COS 邊瓣，通過快速傅里葉變換得到每個(gè)小段的傅里葉譜。最后，對(duì)每個(gè)小段的傅里葉譜進(jìn)行儀器響應(yīng)校正，并通過式（1）得到位移譜振幅a(f)。

其中，vi(f)是經(jīng)過去儀器響應(yīng)之后的第i小段傅里葉速度譜，f指傅里葉變換的各個(gè)頻率點(diǎn)，T為選取P、S 波波形總持續(xù)時(shí)間，n為選取的P、S 波波形所拆分的小段數(shù)量，t為包含256 采樣點(diǎn)小段時(shí)間長(zhǎng)度，t=5.1 s。

對(duì)于噪聲扣除，根據(jù)P 波初動(dòng)前256 個(gè)采樣點(diǎn)的噪聲記錄，通過式（2）可以計(jì)算歸一到與信號(hào)相同持續(xù)時(shí)間的噪聲位移譜振幅。

1.2 位移譜拐角頻率與零頻極限值理論計(jì)算方法

在頻率域?qū)Φ卣鸩ㄟM(jìn)行波譜分析，可得到拐角頻率、位移譜零頻極限值以及高頻衰減斜率。對(duì)于中小地震，震源譜符合Brune 圓盤模型（Brune，1970），在忽略非彈性衰減的情況下，震源位移譜U(f)可表示為

其中，U0為位移譜零頻極限值，f為頻率，fc為拐角頻率。當(dāng)U0和fc確定時(shí)，即可得到理論震源位移譜。在實(shí)踐中，需要根據(jù)已知地震記錄位移譜來反演U0和fc。一般采用擬合法計(jì)算U0和fc（蔡杏輝，2015）。

根據(jù)式（4），在低頻段（f1—f2），U(f)≈U0，則位移譜平方的積分為

根據(jù)式（5）—（6），先在低頻段（f1—f2）積分計(jì)算SD，進(jìn)而求取U0。

1.3 縱橫波波速比計(jì)算

1.3.1 利用縱橫波拐角頻率比值求解。圓盤面位錯(cuò)源產(chǎn)生的地震P 波和S 波遠(yuǎn)場(chǎng)位移的頻譜可以表示為

其中，ρ是介質(zhì)密度，m0是地震矩，α與β分別是P 波與S 波速度，R為震源距，ω為圓頻率分別為P 波、S 波輻射圖形因子，F(xiàn)α(ω)與Fβ(ω)為縱橫波頻譜形狀的函數(shù)，可以統(tǒng)一寫作

式中，c代表α或β，S為圓盤位錯(cuò)面積，為任意一點(diǎn)Ω的位錯(cuò)元dΣ的平面極坐標(biāo)，Δu(Ω)表示位錯(cuò)元在圓盤位錯(cuò)面上的靜力錯(cuò)距，是整個(gè)位錯(cuò)面上的平均值，G(ω)是震源時(shí)間函數(shù)G(t)的頻譜，G(t)是階躍函數(shù)，v表示圓盤中心向外破裂速度。陳運(yùn)泰（1976）給出式（10）的積分結(jié)果，進(jìn)而給出式（7）與式（8）的值，討論了UP(ω)與US(ω)的高頻和低頻漸進(jìn)行為，并給出2 條漸進(jìn)線的交點(diǎn)所對(duì)應(yīng)的拐角頻率，則有

其中，θ為觀測(cè)點(diǎn)在球極坐標(biāo)系(r,θ,φ)的θ項(xiàng)，r0為圓盤形剪切面的破裂半徑，當(dāng)考慮破裂速度v無限大的情況時(shí)，可近似當(dāng)作同時(shí)破裂情形處理，此時(shí)拐角頻率的表達(dá)式可以為對(duì)于在震源球上隨機(jī)分布的觀測(cè)點(diǎn)，則地震拐角頻率期望值可近似為

則由式（13）可得

由于ω=2πf，則利用縱橫波拐角頻率比值，可以計(jì)算得到震源處縱橫波波速比。

1.3.2 利用縱橫波位移譜零頻極限比值求解。位移譜曲線在低頻部分平坦，當(dāng)ω→0 時(shí)，由式（7）—（8）推導(dǎo)得到縱橫波零頻極限值為

比較式（15）—（16）可以得到

由式（17）可知，縱橫波波速比與縱橫波零頻極限比值及輻射圖形因子有關(guān)，但對(duì)于同一序列、同一地區(qū)的地震，可以近似認(rèn)為輻射圖形因子比值為常數(shù)。計(jì)算中，輻射圖形因子可取所有觀測(cè)點(diǎn)輻射圖形因子平均值，一般來說，P 波平均輻射圖形因子為0.52，S 波平均輻射圖形因子為0.63（Aki K et al，1980）。

2 軟件計(jì)算流程

軟件計(jì)算流程為：①導(dǎo)入原始地震編目文件；②導(dǎo)入地震臺(tái)站儀器參數(shù)數(shù)據(jù)，計(jì)算儀器響應(yīng)數(shù)據(jù)；③設(shè)置編目文件的相對(duì)工作目錄，自動(dòng)導(dǎo)入編目對(duì)應(yīng)原始地震的各臺(tái)波形文件；④批量計(jì)算P 波、S 波位移譜、拐角頻率與零頻極限值；⑤由拐角頻率與零頻極限值計(jì)算波速比；⑥循環(huán)計(jì)算完成某地震序列所有地震波速比后，繪制其空間與時(shí)間分布圖。

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

計(jì)算流程中有3 個(gè)基本數(shù)據(jù)需要導(dǎo)入，分別是地震編目數(shù)據(jù)、波形數(shù)據(jù)、臺(tái)站儀器參數(shù)數(shù)據(jù)。以上數(shù)據(jù)均可以從JOPNES 地震速報(bào)分析系統(tǒng)中獲取，地震編目數(shù)據(jù)是目前通用的.phase 數(shù)據(jù)，波形數(shù)據(jù)是JOPNES 地震速報(bào)分析系統(tǒng)截取導(dǎo)出的.sac 數(shù)據(jù)文件。.sac 數(shù)據(jù)格式廣泛用于地震學(xué)分析軟件，也是地震波形數(shù)據(jù)的通用格式有ASCII 碼和二進(jìn)制碼2種格式，本軟件數(shù)據(jù)導(dǎo)入采用.sac 二進(jìn)制碼格式。地震臺(tái)站地震計(jì)參數(shù)也是JOPNES 地震速報(bào)分析系統(tǒng)導(dǎo)出的.ascII 文件中的.par 文件，包含對(duì)應(yīng)地震記錄的地震計(jì)靈敏度、傳遞函數(shù)零極點(diǎn)等信息。導(dǎo)入相關(guān)儀器參數(shù)信息后，可繪制儀器幅頻響應(yīng)圖并直觀顯示，能在對(duì)地震波譜進(jìn)行儀器響應(yīng)校正之前判斷儀器響應(yīng)是否準(zhǔn)確。

2.2 地震P 波、S 波位移譜求取

P、S 波的位移譜求取是計(jì)算拐角頻率與零頻極限值的基礎(chǔ)，需要在軟件中解決以下3個(gè)問題。

2.2.1 P、S 波時(shí)間長(zhǎng)度的確定。P 波時(shí)間長(zhǎng)度為P 波初動(dòng)至S 波初動(dòng)之間波形的時(shí)間段。S 波時(shí)間長(zhǎng)度為從S 波開始到包含90% S 波能量的時(shí)間段。劉杰等（2003）通過回歸分析1 900 條地震記錄，得到S 波截取時(shí)間長(zhǎng)度與Pg、Sg 波走時(shí)差為線性關(guān)系，并擬合出具體線性表達(dá)式。本軟件采用劉杰等（2003）的結(jié)果。在軟件中，通過波形顯示界面拾取Pg、Sg 波走時(shí)差，確定S 波時(shí)間長(zhǎng)度。

2.2.2 扣除儀器響應(yīng)。通過儀器參數(shù)文件計(jì)算儀器的幅頻響應(yīng)，用地震記錄速度譜除以各頻點(diǎn)的儀器幅頻響應(yīng)。

2.2.3 去除背景噪聲。選取Pg 波或Pn 波到時(shí)前256 個(gè)點(diǎn)的時(shí)間段波形，計(jì)算背景噪聲位移譜，根據(jù)式（3），得到所需地震波位移譜。

2.3 波速比的計(jì)算

一個(gè)地震一般會(huì)被多個(gè)地震觀測(cè)臺(tái)站記錄到，可以使用每個(gè)臺(tái)站的波形資料分別計(jì)算拐角頻率fc與零頻極限值U0，再計(jì)算多臺(tái)平均值，最后將平均拐角頻率與平均零頻極限值代入式（14）或式（17）算出波速比。為消除個(gè)別臺(tái)站異常高值對(duì)平均值的影響，在計(jì)算平均值時(shí)采用式（18）—（19）。

式中，xi為各臺(tái)站的拐角頻率fc或零頻極限值U0，N為臺(tái)站數(shù)。Δx為誤差因子，其意義為當(dāng)x以對(duì)數(shù)坐標(biāo)作圖時(shí)的標(biāo)準(zhǔn)差。

3 軟件使用說明

3.1 軟件界面

軟件主界面（圖1）包含基本參數(shù)輸入、計(jì)算方法選擇、結(jié)果顯示?；緟?shù)輸入部分可導(dǎo)入地理底圖、地震編目、儀器參數(shù)等。計(jì)算方法選擇部分可以點(diǎn)選多臺(tái)和達(dá)法或者波譜分析方法來計(jì)算波速比。結(jié)果顯示部分將波速比結(jié)果用等值線繪制在界面中間地理底圖上，正下方坐標(biāo)圖顯示波速比隨時(shí)間變化曲線。在導(dǎo)入各臺(tái)儀器參數(shù)文件之后，軟件界面右上角坐標(biāo)圖可以顯示不同臺(tái)站儀器幅頻響應(yīng)曲線。

圖1 程序界面Fig.1 The graphic interface of the program

如圖2 所示，在扣除儀器響應(yīng)之前，可以觀察地震計(jì)幅頻響應(yīng)曲線以避免錯(cuò)誤。

圖2 地震計(jì)分量?jī)x器響應(yīng)曲線Fig.2 The amplitude frequency response curve of a seismograph

3.2 軟件操作流程

首先，導(dǎo)入基礎(chǔ)地理數(shù)據(jù)，例如研究區(qū)域的省級(jí)區(qū)域底圖，圖2 中顯示的是廣東省地理底圖，還可以選擇添加各類地理信息，如臺(tái)站、縣名、構(gòu)造、河流等，這些地理要素需要事先打包在軟件包里。然后，導(dǎo)入儀器參數(shù)，計(jì)算儀器響應(yīng)。最后，導(dǎo)入地震編目數(shù)據(jù).phase 文件及與地震編目對(duì)應(yīng)觀測(cè)臺(tái)站原始波形記錄.sac 文件。如圖3 所示，數(shù)據(jù)導(dǎo)入完成后，軟件界面中間顯示震中分布，右側(cè)表格顯示地震目錄相關(guān)信息。

圖3 顯示地震信息并繪制震中分布Fig.3 Show the seismic information and plot seismic distribution map

當(dāng)選擇某經(jīng)緯度區(qū)域內(nèi)的地震序列進(jìn)行計(jì)算時(shí)，對(duì)于使用不同臺(tái)站地震波原始記錄計(jì)算的該地震序列的拐角頻率值及零頻極限值結(jié)果，軟件會(huì)自動(dòng)保存圖像，如圖4 所示。

圖4 GD.201810281418 地震事件DNB 臺(tái)地震記錄拐角頻率與零頻極限值結(jié)果Fig.4 The corner frequency and zero spectrum limit value of GD.201810281418 earthquake at DNB station

在地震序列所中有地震的計(jì)算完成后，每個(gè)地震的拐角頻率與零頻極限值以表格形式顯示。

4 實(shí)例分析

4.1 實(shí)例計(jì)算

以2018 年10 月至2020 年2 月廣東陽江區(qū)域19 個(gè)ML2.0 以上地震為例，分別計(jì)算P 波、S 波拐角頻率、零頻極限值及對(duì)應(yīng)波速比，計(jì)算結(jié)果在軟件中彈出窗口以表格顯示，如圖5 所示。使用2 種方法計(jì)算的該序列波速比隨時(shí)間變化曲線如圖6 所示。

圖5 地震序列P 波S 波拐角頻率、零頻極限值及對(duì)應(yīng)波速比結(jié)果Fig.5 The P,S wave corner frequency,zero spectrum limit value,and relative wave velocity ratio result of the earthquake sequence

圖6 拐角頻率法波速比變化曲線與零頻極限值比值法波速比變化曲線Fig.6 The wave velocity ratio curve by P,S wave corner frequency method,and zero spectrum limit value method

為了方便的比較2 種不同方法計(jì)算的波速比隨時(shí)間變化趨勢(shì)，可通過本軟件將計(jì)算結(jié)果歸一化于[0，1]區(qū)間內(nèi)，如圖7 所示。

圖7 歸一化后2 種方法計(jì)算的波速比變化趨勢(shì)對(duì)比曲線Fig.7 The comparative curves of the normalized velocity ratio by two methods

4.2 軟件計(jì)算結(jié)果分析

陳麗娟等（2015）利用拐角頻率比值法得到波速比平均值約1.27，零頻極限比值法得到波速比平均值為1.67。由零頻極限比值法所得波速比較拐角頻率法更接近泊松比。本軟件對(duì)上述地震序列19 例地震計(jì)算結(jié)果顯示，利用拐角頻率比值所得波速比平均值約為1.20，零頻極限比值法所得波速比為1.767 6，與陳麗娟等（2015）的計(jì)算結(jié)果較為一致。且本軟件用零頻極限比值法所得波速比更接近泊松比數(shù)值1.73。對(duì)于2 種結(jié)果的差異，陳麗娟等（2015）解釋為，在零頻極限值與拐角頻率值的反演過程中，對(duì)震源譜的低頻水平段f1、f2和衰減段f3有個(gè)初始的預(yù)設(shè)限定，若預(yù)初始值不恰當(dāng)，可能會(huì)造成某些地震記錄拐角頻率計(jì)算誤差較大，而零頻極限值受初始值影響相對(duì)較小。另外，拐角頻率與破裂過程有關(guān)，且隨臺(tái)站與斷層面之間的夾角變化，而零頻極限值與頻率、破裂過程無關(guān)，所以零頻極限比值法計(jì)算結(jié)果更為穩(wěn)定。

由圖7 可見，零頻極限比值法與拐角頻率比值法所得波速比曲線隨時(shí)間的變化趨勢(shì)有較好的一致性。為了檢查2 種結(jié)果之間的相關(guān)性，利用MATLAB 中自帶corrcoef 函數(shù)，計(jì)算得到零頻極限比值法與拐角頻率比值法的相關(guān)系數(shù)，數(shù)值為0.73，可以認(rèn)為，2 種結(jié)果具有較強(qiáng)的相關(guān)性，2 種不同波譜分析計(jì)算方法能夠相互印證。

中強(qiáng)地震前地震波速比值出現(xiàn)異常。諸多研究表明，地震前縱橫波速比表現(xiàn)為“下降—恢復(fù)”的異常形態(tài)。由本軟件計(jì)算結(jié)果圖可見，特別是在整個(gè)地震序列最大地震ML3.5 地震前和第二大地震ML3.3 地震前，波速比表現(xiàn)出“下降—恢復(fù)”的異常變化，與前人所得結(jié)果一致。

5 結(jié)束語

通過地震波譜分析求取震源處的波速比，其優(yōu)點(diǎn)是不受時(shí)間精度和震相識(shí)別誤差的影響。在地震預(yù)報(bào)實(shí)踐中，可以與常用的多臺(tái)和達(dá)法互相補(bǔ)充。筆者以MATLAB 為開發(fā)平臺(tái)，編制了基于P、S 波拐角頻率與零頻極限值比值的波速比計(jì)算軟件，將導(dǎo)入波形數(shù)據(jù)、去儀器響應(yīng)、求取位移譜、計(jì)算拐角頻率與零頻極限值、結(jié)果可視化及比較集成一體，可用于某個(gè)區(qū)域范圍內(nèi)波速比計(jì)算、大震序列前后震源區(qū)波速比追蹤、判斷未來地震序列發(fā)展趨勢(shì)、判別大震之后強(qiáng)余震的發(fā)生等工作。在實(shí)際應(yīng)用中，為了保證拐角頻率與零頻極限值的計(jì)算準(zhǔn)確性，需要對(duì)研究區(qū)域的地震序列位移譜水平段f1、f2和衰減段f3做預(yù)先判定，設(shè)置較恰當(dāng)?shù)某跏贾?，以期獲得較好的反演結(jié)果。

該軟件使得繁瑣的頻譜計(jì)算過程簡(jiǎn)化，且界面友好，操作簡(jiǎn)單，便于推廣使用。