一種提高振動(dòng)源掃描算法信噪比的方法

劉雙慶，謝 靜，王熠熙

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

0 引言

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

1 振動(dòng)源掃描（SSA）模型

1.1 常用的SSA模型

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

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

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

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

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

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

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

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

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

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

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

2 振動(dòng)源事件的SSA 分析

2.1 有源地震震中定位

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

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

2.2 無(wú)源噪聲定位

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

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

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

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

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

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

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

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

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

3 認(rèn)識(shí)及討論

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

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

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

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

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

［3］ 李文軍，李麗，陳棋福.用震源掃描算法（SSA）研究列車(chē)源的運(yùn)動(dòng)［J］.地球物理學(xué)報(bào)，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.