劉雙慶，謝 靜，王熠熙

（天津市地震局，天津 300201）

0 引言

對產(chǎn)生信號的信號源進(jìn)行有效定位，是被動式觀測系統(tǒng)研究的主要課題。當(dāng)信號源周邊散布多個監(jiān)測儀器，由信號源同一時刻發(fā)出的信號到達(dá)各儀器存在較確切的初至?xí)r間時，可以利用相應(yīng)的走時方程進(jìn)行聯(lián)合定位。在地震科學(xué)中，對地震震源位置進(jìn)行定位即是其中的一個例子。傳統(tǒng)的地震定位，需要較清晰的震相初至到時及明確的震相屬性（Pg、Sg、Pn、P、PkP、PKIKP 等等），以便利用恰當(dāng)?shù)淖邥r方程［1］。對于震相屬性不明確，或初至到時不清晰的振動信號，對其產(chǎn)生源進(jìn)行定位則較困難。為此，一些學(xué)者提出了反向投影算法。這種算法基于如下假定：從振動源發(fā)出的到達(dá)各接收點(diǎn)的信號，其時間序列上存在先后到達(dá)順序。如果消除振動源到達(dá)各接收點(diǎn)的傳播走時，則時間序列上的信號起始點(diǎn)將排列得很整齊。因此再將這些消除走時差異的信號疊加，將有較高的信噪比。這個高能量的信號即為各接收點(diǎn)反算出的震源點(diǎn)某時刻發(fā)出的振動。如果假定的震源點(diǎn)非真實(shí)震源位置，則到達(dá)各接收點(diǎn)的走時與實(shí)際走時不一致，消減這些走時后，各時間序列上的信號起始點(diǎn)排列仍不整齊，疊加后的信噪比較低。因此，在真實(shí)震源附近區(qū)域?qū)δ硶r刻反算出的信噪比進(jìn)行成像，可以獲得真實(shí)震源最可能的位置，即具有最高信噪比的點(diǎn)位。2011年日本MS9.0地震破裂過程反演中，不同頻率的振動源強(qiáng)度分布特征研究就使用了反投影算法［2］。SSA（source－scanning algorithm）算法則是反投影算法中的一種計算實(shí)現(xiàn)模型，它引入“亮點(diǎn)函數(shù)”進(jìn)行刻畫［3］。由于在近場條件及較高頻成分中，隨著震中距離的改變，從同一震源發(fā)出的信號將出現(xiàn)拉長現(xiàn)象，直接矢量疊加效果較差。利用絕對值處理有明顯改善，但是對振動型信號進(jìn)行直接絕對值處理，會引入偽高頻成分。其他學(xué)者還針對振動拉長現(xiàn)象進(jìn)行非等間隔重采樣處理，但改進(jìn)效果不是很明顯［4］。本文將利用局部極值特性，對地震數(shù)據(jù)構(gòu)建較好的包絡(luò)線，以該包絡(luò)線作為SSA 計算的預(yù)處理數(shù)據(jù)，討論這種處理方式的優(yōu)點(diǎn)，并用白噪聲檢驗(yàn)其穩(wěn)定性及其他特性。

1 振動源掃描（SSA）模型

1.1 常用的SSA模型

SSA 模型由Kao和Shan于2004年提出，是一種通過波形對震源分布進(jìn)行定位成像的方法［5］。假設(shè)一個振動事件由N 個臺站記錄到（如圖1 中的A、B、C三個臺站），首先歸一化每個臺站的各通道記錄值un，然后計算某個空間點(diǎn)η在τ時刻的亮度函數(shù)：

圖1 振動源掃描算法原理示意圖（據(jù)文獻(xiàn)［3］修改）空間上的點(diǎn)η在τ時刻的“亮度”可以由相應(yīng)臺站的理論到時（即τ加上相應(yīng)的走時taη、tbη、tcη）計算給出，如果某處“亮度”大（圖中的實(shí)六角星），表明這里有事件發(fā)生，如果“亮度”小（空六角星），表示這里沒有事件發(fā)生

式（1）中：un為歸一化的地震記錄，在本文中取垂直分量的直達(dá)震相波列進(jìn)行分析；tηn為 從 點(diǎn)η到臺站n計算的某個最大能量震相走時，如果所有的最大振幅都是由點(diǎn)η和時間τ產(chǎn)生，那么Bright（η，τ）＝1，即圖1中紅色脈沖都被剛好選擇上。如果脈沖信號都沒被選上，則Bright（η，τ）＝0，即圖1中藍(lán)色空心圈圈定位置。通常情況下，Bright（η，τ）既不等于1也不等于0。另外，由于區(qū)域速度模型跟真實(shí)地殼結(jié)構(gòu)很難完全一樣，為了抑制理論走時誤差的影響，需對求和的N點(diǎn)數(shù)據(jù)進(jìn)行高斯窗加權(quán)處理。

1.2 本文的改進(jìn)模型

地震儀器記錄的地震波列一般呈正負(fù)振蕩衰減。按式（1）直接將小于0值的數(shù)據(jù)進(jìn)行x軸對稱翻轉(zhuǎn)，則原在0點(diǎn)附近近似光滑的波形將不再可導(dǎo)，從而人為引入高頻成分（圖2b）。為此，本文將不直接對負(fù)值區(qū)振幅取絕對值，而是先擬合數(shù)據(jù)的包絡(luò)線。最常用的包絡(luò)線是對數(shù)據(jù)先進(jìn)行等間隔處理，然后提取各間隔區(qū)內(nèi)的最大或最小值作為該區(qū)間的代表值，連接這些值所獲取的曲線作為原始數(shù)據(jù)的包絡(luò)線。由于地震信號是一種非平穩(wěn)時間序列過程，對間隔長度的不同選擇，會對結(jié)果有較大影響，而且等間隔選擇也不能突出局部非平穩(wěn)性。鑒于這種情形，本文提出利用光滑后的局部極值特性進(jìn)行包絡(luò)線追蹤擬合。此外在一般情況下，地震波列還含有環(huán)境噪聲及地震儀器熱噪聲等隨機(jī)干擾，因此，操作過程分以下4步。

（1）按分析的震相（或波列）主要周期段，先進(jìn)行ButterWorth帶通濾波，獲取相對光滑的波形；

（2）對光滑后的數(shù)據(jù)求導(dǎo)，獲取導(dǎo)數(shù)為0的時間軸點(diǎn)；

（3）按奇數(shù)或偶數(shù)腳標(biāo)提取導(dǎo)數(shù)為0的時間點(diǎn)及相應(yīng)的波形數(shù)據(jù)；

（4）對這些數(shù)據(jù)進(jìn)行分段3次埃爾米特插值，獲取數(shù)據(jù)包絡(luò)線。其中兩端無0值導(dǎo)數(shù)的數(shù)據(jù)段按濾波后的數(shù)據(jù)進(jìn)行補(bǔ)充（圖2a）。

經(jīng)過上述處理后，再利用式（1）進(jìn)行SSA 分析（本文所有算例都使用了高斯窗加權(quán)處理）。通過以上方式，可以較好地捕捉到原數(shù)據(jù)的局部極值，同時降低高頻成分的引入。

圖2 對原始波形進(jìn)行兩種不同處理方式的結(jié)果對比

2 振動源事件的SSA 分析

2.1 有源地震震中定位

2014年9月6日18時37分在河北延懷盆地中發(fā)生了一次MS4.3 級地震，其震中定位結(jié)果為（115.433°E，40.278°N，深度約14km）。首都圈臺網(wǎng)中大部分臺站對這次地震都有較高的信噪比記錄。本節(jié)利用該地震對本文的方法進(jìn)行詳細(xì)對比。其中，圖3a是直接采用絕對值方法給出的結(jié)果；圖3b為本文引入的方法的結(jié)果；圖3c為直接使用原始波形進(jìn)行SSA 掃描的結(jié)果，圖中的三角形為參與掃描計算的臺站。計算結(jié)果清晰顯示，直接使用原始波形，亮度函數(shù)值有正負(fù)性，數(shù)值很低（0.001的量級）；用絕對值的方法，亮點(diǎn)函數(shù)的最大值約0.26；而本文所用方法提高到0.41左右。亮度函數(shù)值的空間分布上，本文給出的結(jié)果其邊緣特征更顯著。作者還對掃描求和的長度N 進(jìn)行調(diào)整，其計算的結(jié)果相類似。但如果N 很大，比如大于振動時間長度的2倍，則分辨率較低；N 更長時，分辨率更低。N 宜取分析的臺站中，95%振動能量持續(xù)的平均時間。

圖3 對原始數(shù)據(jù)進(jìn)行不同處理后獲取的SSA 掃描圖像

2.2 無源噪聲定位

當(dāng)存在有源信號時，只要信號發(fā)生時刻較準(zhǔn)確，并且傳播的介質(zhì)速度模型比較合理，則利用絕對值法和本文的方法進(jìn)行SSA 分析一般都能獲取到信號源的有效位置。但是如果振動源產(chǎn)生的待研究信號的發(fā)生時刻未知，信號的信噪比也未明確給出時，對該信號的產(chǎn)生源進(jìn)行定位則非常困難。因此，評價一個觀測系統(tǒng)若干臺儀器記錄的信號中是否存在可探測的同源信息是一個棘手的難題。在此利用無源白噪聲及SSA 算法進(jìn)行分析，對其部分特征進(jìn)行表達(dá)。無源噪聲基礎(chǔ)數(shù)據(jù)采用Matlab軟件自帶的rand函數(shù)進(jìn)行生成［6］。所使用的仍是圖3 中所示的臺站，但每個臺站的偽隨機(jī)數(shù)生成種子并不相同，以保證各道的隨機(jī)數(shù)列幾乎完全互不相關(guān)。圖4給出了直接利用均勻分布的偽隨機(jī)分布噪聲（圖4a）、該噪聲絕對值（圖4b）、對該噪聲處理的本文算法（圖4c，掃描長度擴(kuò)大為2倍）；以及利用噪聲絕對值法，改變掃描長度擴(kuò)大為原長度的2 倍、4 倍、10倍的效果（圖4d～圖4f）；改變掃描起始時間點(diǎn)10 s、20s后的效果（圖4g、圖4h）。其中N 為1 201點(diǎn)，如果采樣率取100sps，則對應(yīng)12s長度。

圖4 利用均勻分布的噪聲對不同掃描長度、起止時間、不同預(yù)處理方法的結(jié)果對比

通過圖3和圖4的結(jié)果對比，可以獲得以下幾點(diǎn)主要認(rèn)識：

（1）當(dāng)不存在發(fā)生源時，不同設(shè)置下所計算的亮度函數(shù)最大最小值差異都很小，一般在0.02左右。而有源存在時，亮度函數(shù)差異在0.2以上，兩者相差一個量級。

（2）噪聲信號計算出的亮度函數(shù)值在0.5左右。

（3）有源信號存在時，直接用絕對值處理的SSA 方法和本文包絡(luò)線的SSA 方法，兩種亮度函數(shù)的結(jié)果分布形態(tài)很相似。而噪聲信號，兩者的結(jié)果差異較大（圖4c、圖4d）。

（4）掃描長度N 的變化對噪聲信號的結(jié)果影響更顯著，當(dāng)N 擴(kuò)大10倍時，幾乎顯示不出N 未變化時的亮度函數(shù)分布特征。

（5）噪聲起始點(diǎn)延遲或提前，亮度函數(shù)的空間分布近似以臺站中心為匯聚或發(fā)散點(diǎn)進(jìn)行匯聚或向外傳播。其原因是延后掃描的起始點(diǎn)時間，相當(dāng)于將噪聲數(shù)據(jù)時間軸提前，因此亮度函數(shù)將向臺站中心匯聚（圖4d、圖4g、圖4h）。而對于有源信號，延后掃描點(diǎn)將導(dǎo)致只有部分臺站能接收到由該源發(fā)出的波列的后半部分，而其他臺站漏記，因此掃描出的亮點(diǎn)函數(shù)將向這些能接收到部分波列的臺站中心匯聚（圖3b、圖3e、圖3f）。即有源信號延遲掃描，亮度函數(shù)高值區(qū)呈有向性偏移。而隨機(jī)信號以匯聚或發(fā)散為主要特征。

（6）掃描起始點(diǎn)偏移，對隨機(jī)信號的亮點(diǎn)函數(shù)值影響不大，仍保持在0.5左右，而有源信號信噪比顯著降低。

3 認(rèn)識及討論

本文利用局部極值構(gòu)建地震儀記錄的原始數(shù)據(jù)包絡(luò)線，在此基礎(chǔ)上對振動源掃描算法SSA 的信噪比進(jìn)行分析。經(jīng)地震事件掃描檢驗(yàn)，本文給出的方法確實(shí)能提高SSA 的信噪比，最大值提高約1.6倍，而且亮度函數(shù)空間分布的邊緣特征更清晰，有利于確定振動源位置。本文的SSA 掃描結(jié)果的高值區(qū)與地震編目的震中位置很吻合。經(jīng)過多次測試顯示：式（1）中的求和長度N 宜取分析的臺站中，95%振動能量持續(xù)的平均時間。N 過短可能導(dǎo)致遠(yuǎn)臺未接收到振動的主要能量，N 過長導(dǎo)致非源振動能量進(jìn)入，且分辨率降低。

本文還對無源均勻噪聲進(jìn)行SSA 掃描，以研究其與有源信號的區(qū)別。計算顯示，無源噪聲掃描的亮度函數(shù)值在0.5附近，而且最大最小值差異很小，約0.02左右。而有源信號一般差異為0.2，二者相差一個量級。對有源信號掃描起點(diǎn)變動后，亮度函數(shù)的高值區(qū)呈方向性偏移，主要偏于能接收到該源部分信號的臺站（網(wǎng)）中心。而無源信號掃描起點(diǎn)變動后，亮度函數(shù)表現(xiàn)為以臺站（網(wǎng)）為中心而進(jìn)行的發(fā)散或聚斂。而且有源信號的亮點(diǎn)函數(shù)值顯著下降，但隨機(jī)信號的亮度函數(shù)值未出現(xiàn)明顯變化。隨著掃描長度N 的增長，隨機(jī)信號的亮度函數(shù)空間分布特征持續(xù)性較差。

當(dāng)然，引起SSA 算法的信噪比變化的因素很多，其與信號源的自身特征也有很大的關(guān)系。本文從幾個主要參數(shù)進(jìn)行對比分析，雖然較清晰地展示了本文提出的方法的優(yōu)點(diǎn)，以及對無源隨機(jī)信號判別的一些特征，但對不同信噪比的有源信號的分析仍未進(jìn)行展開，這是后續(xù)需要進(jìn)一步強(qiáng)化的內(nèi)容。

［1］ 國家地震局地球物理研究所.近震分析［M］.北京：地震出版社，1978：3－8.

［2］ 姚華建.用基于時間域和頻率域的反投影算法研究大地震的破裂過程［R］.北京：中國地震局地震預(yù)測研究所學(xué)術(shù)報告，2013.

［3］ 李文軍，李麗，陳棋福.用震源掃描算法（SSA）研究列車源的運(yùn)動［J］.地球物理學(xué)報，2008，51（4）：1146－1151.

［4］ 張瑞紅，林大超，喬蘭.一種改進(jìn)的震源掃描算法微震定位［J］.煤礦安全，2011，42（5）：48－51.

［5］ Honn Kao，Shao－ju Shan.The source－scanning algorithm：mapping the distribution of seismic sources in time and space［J］.Geophys.J.omt.，2004，157：589－594.

［6］ Cleve B.Moler.Numerical Computing with MATLAB［M］.America：Society for Industrial and Applied Mathematics，2004：255－264.