基于波譜分析方法的波速比計算軟件

高 亮 湯曜煒 楊 選

（中國廣州510070 廣東省地震局）

0 引言

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

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

1 波速比求解理論

1.1 地震波形譜分析

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

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

對于噪聲扣除，根據(jù)P 波初動前256 個采樣點的噪聲記錄，通過式（2）可以計算歸一到與信號相同持續(xù)時間的噪聲位移譜振幅。

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

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

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

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

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

1.3 縱橫波波速比計算

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

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

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

其中，θ為觀測點在球極坐標系(r,θ,φ)的θ項，r0為圓盤形剪切面的破裂半徑，當考慮破裂速度v無限大的情況時，可近似當作同時破裂情形處理，此時拐角頻率的表達式可以為對于在震源球上隨機分布的觀測點，則地震拐角頻率期望值可近似為

則由式（13）可得

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

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

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

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

2 軟件計算流程

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

2.1 數(shù)據(jù)導入

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

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

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

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

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

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

2.3 波速比的計算

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

式中，xi為各臺站的拐角頻率fc或零頻極限值U0，N為臺站數(shù)。Δx為誤差因子，其意義為當x以對數(shù)坐標作圖時的標準差。

3 軟件使用說明

3.1 軟件界面

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

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

如圖2 所示，在扣除儀器響應之前，可以觀察地震計幅頻響應曲線以避免錯誤。

圖2 地震計分量儀器響應曲線Fig.2 The amplitude frequency response curve of a seismograph

3.2 軟件操作流程

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

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

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

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

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

4 實例分析

4.1 實例計算

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

圖5 地震序列P 波S 波拐角頻率、零頻極限值及對應波速比結果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 種不同方法計算的波速比隨時間變化趨勢，可通過本軟件將計算結果歸一化于[0，1]區(qū)間內(nèi)，如圖7 所示。

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

4.2 軟件計算結果分析

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

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

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

5 結束語

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

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