沈紅薇，李 環(huán)

(沈陽(yáng)理工大學(xué) 信息科學(xué)與工程學(xué)院，沈陽(yáng) 110159)

地震波初至波與scholte波的時(shí)差提取研究

沈紅薇，李 環(huán)

(沈陽(yáng)理工大學(xué) 信息科學(xué)與工程學(xué)院，沈陽(yáng) 110159)

通過(guò)對(duì)初至波與scholte波到達(dá)時(shí)刻提取的方法研究，提出計(jì)盒分形維數(shù)算法和修正后的平滑偽魏格納分布時(shí)頻分析算法分別對(duì)初至波和scholte波到達(dá)時(shí)刻進(jìn)行提取，兩者時(shí)差主要應(yīng)用于水下艦船目標(biāo)探測(cè)研究。在此背景下，通過(guò)模擬淺海海底環(huán)境進(jìn)行地震波數(shù)據(jù)采集，并拾取初至波與scholte波，獲取兩者時(shí)差，從而驗(yàn)證了所提算法的有效性。

初至波；scholte波；計(jì)盒分形維數(shù)法；修正的平滑偽魏格納分布

隨著現(xiàn)代艦船防護(hù)和隱身能力的不斷增強(qiáng)，水雷及聲吶等傳統(tǒng)探測(cè)技術(shù)對(duì)艦船目標(biāo)識(shí)別與探測(cè)難度日益加大，迫切需要某種新型的水下艦船探測(cè)及定位技術(shù)來(lái)適應(yīng)現(xiàn)代科技的發(fā)展。本文利用艦船地震波測(cè)距這種目前在水下探測(cè)技術(shù)方面鮮少用到的方法進(jìn)行目標(biāo)識(shí)別與定位。

艦船在水下行進(jìn)時(shí)會(huì)產(chǎn)生巨大的振動(dòng)，進(jìn)而產(chǎn)生噪聲信號(hào)，這種振動(dòng)波會(huì)以水體為媒介傳播到海底，引起海床振動(dòng)，并在海底地層表面以彈性波的形式向四周傳播。這種以行進(jìn)的艦船為震源所產(chǎn)生的彈性波和普通的地震波相同，被稱為艦船海底地震波。地震波又分為體波和面波兩類，體波又分為橫波和縱波[1]。在海洋地震測(cè)量中，均勻地層中地震觀測(cè)點(diǎn)最先接收到的由震源直接傳過(guò)來(lái)的地震波稱為初至波[2]。而沿著彈性介質(zhì)表面或兩個(gè)不同彈性介質(zhì)表面上傳播的地震波稱為面波，它是介質(zhì)中橫波與縱波相互疊加與干涉所產(chǎn)生的波。本文選用初至波與面波的時(shí)差進(jìn)行測(cè)距進(jìn)而實(shí)現(xiàn)水下艦艇目標(biāo)定位。初至波相對(duì)于橫波和面波起跳時(shí)間最早，到達(dá)時(shí)間最快，頻率最高，幅值最低。正因?yàn)槌踔敛ň哂羞@些特性，因而大部分初至波拾取算法中均利用這些特點(diǎn)來(lái)進(jìn)行研究。面波主要能量分布在10～80Hz頻帶內(nèi)。主要在地表傳播，其能量約占地震波總能量的百分之七十左右，波速約為1～3km/s，相比于體波波速較低。面波根據(jù)其傳播時(shí)媒介特質(zhì)的不同而有不同的面波分類，當(dāng)僅在固體表面上傳播時(shí)稱為瑞利波，在固體與固體界面上傳播時(shí)稱為斯通利波，而沿著液體與固體界面上傳播時(shí)稱為scholte波[3]。由于本文研究所得結(jié)論旨在用在海底艦船目標(biāo)探測(cè)及定位上，故重點(diǎn)研究scholte波的特性及拾取。scholte波沿著海底表平面?zhèn)鞑?，離開界面時(shí)波在這兩種介質(zhì)中逐漸消失，幅度則按指數(shù)規(guī)律進(jìn)行衰減。質(zhì)點(diǎn)在深度-距離平面內(nèi)做橢圓運(yùn)動(dòng)，沒有低頻截止頻率。且具有到達(dá)時(shí)間最晚、頻率低、幅值居中、能量最大等特點(diǎn)。scholte波所在頻帶內(nèi)40～70Hz的能量占主要地位，故scholte波的提取通常采用時(shí)頻分析的方法，觀察其能量及頻率分布來(lái)提取其到達(dá)時(shí)刻，進(jìn)而測(cè)量其距離可實(shí)現(xiàn)水下艦船目標(biāo)定位。

1 計(jì)盒維數(shù)法提取初至波原理

1.1 分形維數(shù)算法

本文采用計(jì)盒分形維數(shù)算法來(lái)提取初至波的到達(dá)時(shí)刻。分形維數(shù)由法國(guó)數(shù)學(xué)家Mandelbrot于1973年率先引入到信號(hào)處理方面進(jìn)行研究，它研究極不規(guī)律、極不連續(xù)的幾何問題。它不同于以往人們?cè)跉W式空間所研究問題的習(xí)慣，而是用極簡(jiǎn)單的方法描述復(fù)雜現(xiàn)象，用規(guī)律的方法來(lái)處理不規(guī)律的問題，它為一些非均勻型、突變型、差異性和間斷型的研究提供了一種非常有效的方法[4]。分形維數(shù)的計(jì)算方法有許多，方法不同則得到的維數(shù)值的含義也不盡相同，在關(guān)于地理問題的研究中，計(jì)盒維數(shù)法是相對(duì)用的較多的一種分形維數(shù)算法。取邊長(zhǎng)為r的多個(gè)盒子，用其將分形曲線覆蓋起來(lái)。覆蓋全部曲線所得的盒子數(shù)記為N，由此計(jì)算其分形維數(shù)及其變化趨勢(shì)。通過(guò)幾中不同的分形維數(shù)計(jì)算方法的對(duì)比，不難發(fā)現(xiàn)計(jì)盒維數(shù)法更適用于地震道信號(hào)的提取與研究。

1.2 地震數(shù)據(jù)的分形維數(shù)計(jì)算

當(dāng)初至波到達(dá)以后，地震記錄的信號(hào)變成了噪聲信號(hào)與初至波信號(hào)的疊加，其能量更強(qiáng)，振幅也較大。通過(guò)對(duì)比初至波到達(dá)前、后計(jì)盒分形維數(shù)的變化，可以較準(zhǔn)確地描繪出其對(duì)應(yīng)時(shí)間序列的能量變化趨勢(shì)，分形維數(shù)的第一個(gè)突變點(diǎn)即指出了地震波初至到達(dá)相對(duì)應(yīng)的時(shí)間位置。

定義n個(gè)邊長(zhǎng)為r的盒子覆蓋地震道曲線，即其測(cè)量長(zhǎng)度為L(zhǎng)=n×r，改變盒子邊長(zhǎng)r，則測(cè)量的分形曲線長(zhǎng)度L也隨之變化，發(fā)現(xiàn)盒子邊長(zhǎng)r與分形曲線長(zhǎng)度L存在一定的指數(shù)關(guān)系，即：

L=C×r1-D

(1)

式中,C是比例常數(shù)，D為所要求得的分形維數(shù)。

在實(shí)際算法運(yùn)用中即在某個(gè)需要處理的地震道上確定長(zhǎng)度為r的一個(gè)工作窗，將工作窗逐個(gè)沿地震道向前移動(dòng)直至覆蓋整個(gè)地震道曲線，計(jì)算每個(gè)工作窗內(nèi)地震道采樣序列的分形維數(shù)，得到一條分形維數(shù)變化軌跡趨勢(shì)圖，此軌跡第一個(gè)突變點(diǎn)即為初至波到達(dá)的時(shí)間點(diǎn)。改變窗口長(zhǎng)度r，不同的窗口長(zhǎng)度會(huì)對(duì)應(yīng)得到不同的分形維數(shù)，r越小則計(jì)算精度越高，但會(huì)有計(jì)算時(shí)間上的限制。將窗口長(zhǎng)度r和分形曲線長(zhǎng)度L繪制在同一雙對(duì)數(shù)坐標(biāo)平面上，代入式(1)，得到分形維公式如下：

logL=(1-D)logr+logC

(2)

擬合各點(diǎn)曲線的斜率和它所分析的分形維數(shù)，不難發(fā)現(xiàn)兩者間存在一定的關(guān)系，如圖1所示。

圖1 窗口長(zhǎng)度與分形曲線對(duì)數(shù)關(guān)系

即斜率k=logL/logr，故分形維數(shù)與此擬合直線的斜率關(guān)系為D=1-k，因此分形維數(shù)也可以利用分形曲線長(zhǎng)度L與窗口長(zhǎng)度r的對(duì)數(shù)比來(lái)求得。這和利用最小二乘法對(duì)式(2)來(lái)估計(jì)分形維數(shù)D的原理相同。由此分形維數(shù)理論可知，利用計(jì)盒維數(shù)算法提取初至波可行。系統(tǒng)會(huì)自動(dòng)計(jì)算出分形維數(shù)突變點(diǎn)所在位置，再根據(jù)地震道數(shù)據(jù)的采樣序列查找出突變點(diǎn)所對(duì)應(yīng)的初至波到達(dá)時(shí)刻。

2 修正的平滑偽魏格納分布提取scholte波原理

時(shí)頻分析的主要研究對(duì)象是非平穩(wěn)信號(hào)或時(shí)變信號(hào)，它能清晰地描繪出信號(hào)的時(shí)變頻譜特性及能量分布，在信號(hào)處理中有著非常重要的作用。典型的時(shí)頻分析方法有S變換、小波變換、短時(shí)傅里葉變換、魏格納分布、平滑魏格納分布等。

時(shí)頻分析是利用時(shí)間和頻率的聯(lián)合分布函數(shù)來(lái)表示某一非平穩(wěn)信號(hào)，并對(duì)此非平穩(wěn)信號(hào)進(jìn)行分析和處理的一個(gè)過(guò)程[5]。地震波信號(hào)是一種非平穩(wěn)信號(hào)，因此利用時(shí)頻分析能更準(zhǔn)確地將時(shí)間、頻率和能量清晰地表現(xiàn)出來(lái)，以此最終精確地計(jì)算出面波scholte波的到達(dá)時(shí)刻。

本文采用修正的魏格納分布來(lái)對(duì)截取后的濾波信號(hào)進(jìn)行時(shí)頻分析，觀察其頻率、時(shí)間及能量分布的情況來(lái)確定scholte波到達(dá)時(shí)刻。

魏格納分布是由魏格納在二十世紀(jì)初提出的，起初被應(yīng)用到量子力學(xué)方面，此后由Ville將它引進(jìn)到信號(hào)分析處理方面進(jìn)行算法研究。魏格納分布是完全非局部的，能更好地適用于非平穩(wěn)信號(hào)的時(shí)頻分析。

(3)

式中:W(t,ω)是信號(hào)s(t)的魏格納分布;t為時(shí)移;τ為積分變量。

根據(jù)希爾伯特變換對(duì)s(t)構(gòu)造解析信號(hào)z(t)形式。希爾伯特變換公式如下：

z(t)=s(t)+jH[s(t)]

(4)

式中H[s(t)]定義為

(5)

而在頻域上，由實(shí)信號(hào)s(t)構(gòu)造的解析信號(hào)z(t)魏格納分布為

(6)

由于魏格納分布具有良好的時(shí)頻聚焦性，有時(shí)魏格納分布在時(shí)間和頻率上把某些值置于兩個(gè)信號(hào)的中間，有時(shí)這些值又處于時(shí)頻平面與所希望得到的信號(hào)爭(zhēng)奪位置，這樣就產(chǎn)生了交叉項(xiàng)，是指不同信號(hào)分量之間的相互交叉作用所產(chǎn)生的“虛假信號(hào)”[6]。交叉項(xiàng)的幅值可以達(dá)到其信號(hào)本值的兩倍，該現(xiàn)象會(huì)嚴(yán)重影響到信號(hào)的時(shí)頻分析特征[7]。所以，為抑制這種“虛假信號(hào)”人們提出多種解決辦法，如平滑偽魏格納分布、choi-william分布等。目前應(yīng)用較為廣泛的是平滑偽魏格納分布。

信號(hào)s(t)的平滑偽魏格納分布為

(7)

式中g(shù)(u-t)、h(τ)均為窗函數(shù)，用來(lái)抑制交叉項(xiàng)，其中當(dāng)且僅當(dāng)g(0)=h(0)=1時(shí)，g(u)和h(τ)為兩個(gè)實(shí)偶窗函數(shù)，魏格納分布經(jīng)過(guò)加窗處理后即為平滑偽魏格納分布。

若將平滑偽魏格納分布適當(dāng)重新排列，則會(huì)進(jìn)一步改善其時(shí)頻分析的性能，這種特殊形式的平滑偽魏格納分布被稱作修正的平滑偽魏格納分布。其分布公式如下所示：

(8)

3 初至波與scholte波的時(shí)差提取實(shí)驗(yàn)與分析

3.1 采集系統(tǒng)的組成

由于條件限制，無(wú)法利用真正的艦船來(lái)拾取其艦船地震波進(jìn)行試驗(yàn)，故模擬了一個(gè)淺海海底環(huán)境來(lái)進(jìn)行測(cè)試。

實(shí)驗(yàn)場(chǎng)地：某水聲研究所院內(nèi)，選擇外界干擾較少時(shí)進(jìn)行。

實(shí)驗(yàn)器材：7m深水池，方形鐵板(約70cm×80cm)，18.6kg大鐵錘(用于敲擊鐵板作為震源)。

實(shí)驗(yàn)過(guò)程：將方形鐵板浮在水池上，利用大鐵錘從高于水面5m左右的空中自由落體砸在鐵板上，以此作為震源放在右側(cè)聲源口上。傳感器置于地面連接采集系統(tǒng)進(jìn)行數(shù)據(jù)采集。

采樣頻率：2048Hz。

其他硬件實(shí)驗(yàn)器材：CA-YD-189型傳感器、pxi-3342型采集卡、泛華PXI-9106型機(jī)器、YE-3821型信號(hào)調(diào)理器。

所用軟件：Lzbview2013、Matlab2012b。

3.2 數(shù)據(jù)采集和數(shù)據(jù)處理

圖2為采用泛華 PXI-9106型號(hào)機(jī)器和Labview2013軟件編譯的數(shù)據(jù)采集系統(tǒng)前面板。

圖2 數(shù)據(jù)采集程序前面板

測(cè)量幅值范圍根據(jù)采集卡固有的電壓范圍，設(shè)置為最小-5V最大+5V,采樣頻率為2048Hz。利用上述采集系統(tǒng)采集到的原始地震數(shù)據(jù)波形如圖3所示。

圖3 采集到的原始數(shù)據(jù)波形

圖4為截取其中部分有用信號(hào)所得數(shù)據(jù)波形。

圖4 截取后的數(shù)據(jù)波形

3.3 初至波到達(dá)時(shí)刻提取

利用計(jì)盒分形維數(shù)算法對(duì)截取后的部分原始地震波數(shù)據(jù)進(jìn)行初至波提取。在提取過(guò)程中改變盒子邊長(zhǎng)即窗口長(zhǎng)度r分別計(jì)算其分形維數(shù)。并通過(guò)其變化趨勢(shì)來(lái)確定最適合的窗口長(zhǎng)度，進(jìn)而更為準(zhǔn)確地提取初至波。通過(guò)實(shí)際程序中r的更改對(duì)比，當(dāng)r=20時(shí)，既是采樣頻率的整數(shù)倍，且能在計(jì)算機(jī)計(jì)算精度較高的情況下計(jì)算時(shí)間較短，因此選則窗口長(zhǎng)度r為20進(jìn)行初至波到達(dá)時(shí)刻拾取。

圖5為對(duì)截取后的數(shù)據(jù)進(jìn)行無(wú)量綱化處理后，利用計(jì)盒維數(shù)法得到的窗口長(zhǎng)度為20時(shí)地震信號(hào)分形維數(shù)趨勢(shì)圖。

由圖5可見，圖中程序自動(dòng)標(biāo)出的第一個(gè)突變點(diǎn)即為初至波到達(dá)時(shí)刻所對(duì)應(yīng)的分形維數(shù)離散點(diǎn)。根據(jù)程序計(jì)算結(jié)果如圖6所示，按其窗口長(zhǎng)度為20計(jì)算其分形維數(shù)突變點(diǎn)在第15點(diǎn)。其對(duì)應(yīng)的實(shí)際地震信號(hào)的離散點(diǎn)為15+1再乘以窗口長(zhǎng)度20即320點(diǎn)，對(duì)應(yīng)實(shí)際信號(hào)數(shù)據(jù)如圖7所示，可得其初至波到達(dá)時(shí)刻為17.055664s。

3.4 scholte波到達(dá)時(shí)刻提取

由于scholte波主要能量分布在10～80Hz之間，故對(duì)圖4的截取后的信號(hào)波進(jìn)行濾波去噪，采用階數(shù)為3階的巴特沃茲帶通濾波器進(jìn)行濾波，通帶范圍為10～80Hz，由此截取其頻帶在10～80Hz之間的頻率信號(hào)。濾波后的截取信號(hào)如圖8所示。

圖5 地震信號(hào)分形維數(shù)趨勢(shì)圖

圖6 分形維數(shù)突變處對(duì)應(yīng)的離散點(diǎn)

圖7 初至波到達(dá)時(shí)刻數(shù)據(jù)

圖8 濾波后信號(hào)

再將圖8所示濾波后的有用信號(hào)利用修正后的平滑偽魏格納分布進(jìn)行時(shí)頻分析，得到圖9所示的時(shí)-頻能量分布圖。可以準(zhǔn)確地看出10～80Hz頻帶內(nèi)信號(hào)的時(shí)間、頻率及能量密度的分布情況。顏色越深則表示能量越大。由于scholte波同其他干擾波相比能量最大，頻率居中，且已經(jīng)10～80Hz外的信號(hào)濾掉，因此將其時(shí)頻分析能量分布圖轉(zhuǎn)換視角并進(jìn)行部分放大如圖10后，可明顯觀測(cè)到在6.72s時(shí)刻scholte波到達(dá)，成功接收到scholte波。

圖9 時(shí)頻分析能量分布圖

圖10 放大后的能量分布圖

3.5 時(shí)差提取

由以上結(jié)論可知，對(duì)于所采集的實(shí)際地震信號(hào)進(jìn)行數(shù)據(jù)處理后，提取到其初至波到達(dá)時(shí)刻為17.056s，其scholte波到達(dá)時(shí)刻為17.42s，兩者時(shí)差為0.364s。在此初至波與scholte波時(shí)差提取成功的基礎(chǔ)上，可以設(shè)置多個(gè)傳感器在不同方位進(jìn)行數(shù)據(jù)采集，利用時(shí)延定位法計(jì)算每個(gè)傳感器對(duì)該數(shù)據(jù)的接收時(shí)間時(shí)差和已知的波速得到目標(biāo)到傳感器的距離進(jìn)而對(duì)目標(biāo)進(jìn)行定位。

4 結(jié)論

(1)利用計(jì)盒分形維數(shù)法能較準(zhǔn)確地提取初至波到達(dá)時(shí)刻，且效率明顯提高。

(2)利用修正后的平滑偽魏格納分布進(jìn)行scholte波的提取，能很好地抑制交叉項(xiàng)，且很好地表現(xiàn)了能量的分布情況。結(jié)合scholet波能量最強(qiáng)頻率居中等特點(diǎn)進(jìn)行時(shí)頻分析，這樣可以非常清晰地提取scholte波到達(dá)時(shí)刻，進(jìn)而得到初至波與scholte波的到達(dá)時(shí)差，以便實(shí)現(xiàn)水下目標(biāo)定位。

[1]李響,顏冰,周穗華.基于直達(dá)波與地聲界面波時(shí)延的目標(biāo)定位方法[J].探測(cè)與控制學(xué)報(bào),2010,32(6):50-53.

[2]曾富英,李敏鋒,申維.地震波初至波拾取的分形研究[J].現(xiàn)代地質(zhì),2002,16(2):209-213.

[3]張海剛.淺海甚低頻聲傳播建模與規(guī)律研究[D].哈爾濱:哈爾濱工程大學(xué),2010.

[4]錢光萍.基于分形維地震道初至拾取方法研究[D].成都:成都理工大學(xué),2001.

[5]馬見青.多分量地震自適應(yīng)極化濾波方法[D].西安:長(zhǎng)安大學(xué),2012.

[6]田琳,陳穎頻,梁華蘭.平滑偽Wigner-Ville分布在地震信號(hào)處理中的應(yīng)用[J].新疆師范大學(xué)學(xué)報(bào):自然科學(xué)版,2013,32(3):1-4.

[7]金銀燕,于鳳芹.Gabor變換與雙線性時(shí)頻分布的時(shí)頻結(jié)構(gòu)[J].計(jì)算機(jī)工程與應(yīng)用,2011,47(25):146-148.

(責(zé)任編輯：馬金發(fā))

Research on Time Difference Extraction of Seismic Preliminary Wave and Scholte Wave

SHEN Hongwei,LI Huan

(Shenyang Ligong University,Shenyang 110159,China)

By researching the method of preliminary wave and scholte wave of arrival time extraction,a box counting fractal dimension algorithm is proposed to extract the time of preliminary wave arrival and give a correction to extract the time of scholte wave arrival after the smoothed pseudo Wigner Ville distribution frequency analysis method.The time difference between the two waves is mainly used in underwater ship target survey research.In this context,the shallow seabed environment model system is taken to extract preliminary wave and scholte wave with time difference for seismic wave data collection，which verifies the algorithm validity.

preliminary wave；scholte wave;box counting fractal dimension;modified smoothed pseudo Wigner-Ville transform

2016-01-08

沈紅薇(1990—)，女，碩士研究生；通訊作者：李環(huán)(1964—)，女，教授，研究方向：擴(kuò)頻通信技術(shù)及應(yīng)用等。

1003-1251(2016)06-0082-06

TP391

A