蔣 禮

(1.中國地質(zhì)大學(xué)地球物理與空間信息學(xué)院,武漢 430074;2.華北水利水電學(xué)院數(shù)學(xué)與信息科學(xué)學(xué)院,鄭州 450011)

1 引 言

傳統(tǒng)面波的頻散分析方法主要分為譜分析法(SASW)[1]和多道分析法(MASW)[2]兩類。雖然SASW具有原理簡單、運(yùn)算速度快等優(yōu)點(diǎn),但由于下列兩個主要原因的影響導(dǎo)致該算法的準(zhǔn)確性遠(yuǎn)低于MASW。第一,無論是用時差還是相差來提取群速度或相速度,由于時差和相差都是位于分母,故而任何細(xì)微誤差都會對求取的速度造成很大的影響;第二,譜分析所得到的相位并非該信號的真實相位,故而存在相位解纏的難題。

本文首先將希爾伯特黃變換應(yīng)用于時頻分析,從而獲得比傳統(tǒng)傅里葉譜更為精確的希爾伯特譜,進(jìn)而得到較為準(zhǔn)確的群速度信息;然后,依據(jù)群速度和相速度的關(guān)系,提出在群速度的基礎(chǔ)上求取相速度的新算法,該算法可以有效地解決相位解纏難題,并能最終獲得較為準(zhǔn)確的相速度信息。

2 原理簡介

通常有3類方法進(jìn)行時頻分析:基于時頻窗類的短時傅里葉變換、gabor變換和小波變換[3];基于相關(guān)性研究的Cohen類時頻分布[3];基于瞬時頻率的希爾伯特譜分析[3,4]。希爾伯特譜不受時頻不確定原理的限制,故該譜比傅里葉譜有著更精確的時頻刻畫能力;而且希爾伯特譜也不必顧慮類似于Cohen類時頻分布交叉項、邊緣條件等核干擾的問題,故該譜比Cohen分布更接近于該信號的真實物理頻率(即傅里葉頻率)。

2.1 希爾伯特黃變換

2.1.1 瞬時頻率

[3]

設(shè)某實信號 s(t)=A(t)ejφ(t),其中 A(t)表示瞬時幅度,φ(t)表示瞬時相位,其解析信號的傅里葉變換結(jié)果為S(ω),則其平均頻率為

因此瞬時頻率為

由式(1)可知,瞬時頻率即為某一時刻的平均頻率,如果該信號為單分量信號,則 φ/(t)即為真實的物理頻率(即傅里葉頻率)。

2.1.2 經(jīng)驗?zāi)B(tài)分解(EMD)

將一個多分量信號分解為一系列單分量信號之和的過程稱之為經(jīng)驗?zāi)B(tài)分解[5,6],每一個單分量信號稱為一個固有模態(tài)函數(shù)(IMF)。

對實信號s(t)進(jìn)行經(jīng)驗?zāi)B(tài)分解可得:

2.1.3 希爾伯特譜

對于各單分量信號,利用希爾伯特變換(HT)可以得到基于瞬時頻率的希爾伯特譜。希爾伯特譜H(t,ω)主要有兩類,即時間-頻率-能量譜He(t,ω)和時間-頻率-相位譜Hp(t,ω),從這兩類譜中可以獲得時間、頻率、能量和相位的分布關(guān)系。

希爾伯特黃變換的完整算法流程如下:

2.2 頻散分析

以地震波信號為例,設(shè)震源的單分量信號[1]為s(t)=A(t)ejφ(t),道間距(即相鄰檢波器的距離)為Δx,某道的該分量信號為

相鄰道的該分量變?yōu)?/p>

式中,tg為群傳播的時差,tp為相傳播的時差,則群速度 Vg=Δx/(tg2-tg1),相速度Vp=Δx/(tp2-tp1)。

2.2.1 群速度

時頻能量譜He(t,ω)為傳統(tǒng)的時頻分析譜,從中可以獲得各分量信號的能量分布情況和群速度計算結(jié)果[7]。

tg為群傳播時差,即能量傳播的時差。根據(jù)時頻譜He(t,ω),某單分量信號IMFn的能量傳播時差為

而該單分量信號的傅里葉頻率即平均頻率為

該單分量信號的振幅計算公式為

式中,f0、f1分別為該單分量信號瞬時頻率的起、止頻率,Tmean=1/fmean為該單分量信號的周期。

2.2.2 相速度

從時頻相位譜Hp(t,ω)中可以獲得各分量信號在任意時刻的瞬時相位,再根據(jù)相差計算公式可以求取各分量信號的相速度。

(1)同步相差[1,8],即同時在兩道檢波器獲取同一分量檢測信號的相位差。

式中,Δφ為同步相位差,Δx為道間距。

由于低頻信號周期較長,因此前后兩道同時檢測到該信號的相位值(相位值的取值范圍在±2π內(nèi))依然有可能位于同一周期內(nèi),此時并不需要進(jìn)行相位解纏;而中頻和高頻信號因為周期較短,所以前后兩道同時檢測到該信號的相位值很有可能已不在同一周期內(nèi),此時必須進(jìn)行相位解纏,才能獲得真實的相位差。

(2)異步相差,即相鄰兩道檢波器在不同時刻所得同一分量信號的相位差。

式中,Δφ為異步相位差。

假設(shè)沿著波的傳播方向,先到達(dá)的檢波器記錄下某分量信號某一時刻的相位,后到達(dá)的檢波器記錄下該分量信號略微延遲一下的相位,只要檢測時差控制得當(dāng),任何頻率的相位差都可控制在同一周期內(nèi),所以使用恰當(dāng)時差的異步相差完全可以避免相位解纏。

(3)校正相差,即利用群速度傳播時差代替估計時差的異步相差。

雖然使用恰當(dāng)時差的異步相差可以完全避免相位解纏,但測量時差的選取只能進(jìn)行估計嘗試。一旦該時差的選擇出現(xiàn)少許偏差,由于中頻、高頻信號周期較短,此時依然會面臨相位解纏;如果時差選擇反了,有可能低頻信號也要進(jìn)行相位解纏,而中高頻信號將面臨著更復(fù)雜的相位解纏。對于地球物理勘探常用的瑞利面波信號而言,各頻率信號的群速度和相速度的變化規(guī)律基本一致(隨著頻率的增加而減小),故各分量的群傳播時差和相傳播時差的區(qū)別通常不會太大。

校正相差為

式中,k為波數(shù)。若Vp與Vg相近,則Δφ可控制在±2π內(nèi),此時可以避免相位解纏,將式(8)計算結(jié)果直接代入式(7)即可求出相速度。校正相差可以有效解決估計時差問題,并且可以在一定程度上避免相位解纏。

3 仿真實例

3.1 模型建立

設(shè)某個瑞利波源是由3個單分量信號組合而成,各分量的具體參數(shù)詳見表1。從該表中可以看出該模擬波低頻分量能量大、高頻分量能量小,各分量信號的群速度和相速度隨著頻率的增加而減小,且群速度小于相速度,這些特征十分符合地震瑞利面波特性。

表1 瑞利波源分量參數(shù)列表Table 1 Rayleigh wave source component parameter list

兩個檢波器分別距離震源12 m和18 m,檢波時間為1.023 s,采樣時間間隔為1ms。12 m處檢波器所得信號如圖1(a)所示,18 m處的信號如圖1(b)所示。從圖中可以看出,信號所占時寬隨著傳播距離增加正在慢慢變大,各分量信號在傳播過程中因速度不同正在互相分離,這是典型的頻散效果。對這兩個信號進(jìn)行希爾伯特黃變換,并繪制出相應(yīng)的時間-頻率-能量譜(如圖 1(c)、(d)所示)和時間-頻率-相位譜(如圖1(e)、(f)所示)。

圖1 希爾伯特黃變換譜分析圖Fig.1 Spectrum analysis diagram of the Hilbert-Huang transform

3.2 希爾伯特譜分析

3.2.1 相位譜和能量譜

時間-頻率-能量譜簡稱能量譜,從圖1(c)和圖1(d)中可以看出該譜橫軸為時間軸,縱軸為頻率軸,色標(biāo)代表著多點(diǎn)平滑幅值也即能量。該圖可以看作是傳統(tǒng)的時頻分析圖,只不過此圖的頻率為瞬時頻率;時間-頻率-相位譜簡稱相位譜,該譜與能量譜類似(見圖1(e)和圖1(f)),橫軸為時間軸,縱軸為頻率軸,但色標(biāo)代表著相位。

從相位譜和能量譜的對比可以看出,不是所有的瞬時相位都有意義。對于某一單分量信號而言,有能量的相位或者能量大于某一設(shè)定閾值的相位,可看作是該分量信號的真實相位;而其余的相位可看作是能量泄漏在該頻段產(chǎn)生的干擾相位,這部分相位沒有真實的物理意義。對于同一信號使用不同的EMD方法可能會出現(xiàn)不同的干擾相位,但真實相位始終相同。

無論從能量譜還是相位譜中我們都可以清楚地看到3個單分量信號的信息,各分量的瞬時頻率非常接近于真實的物理頻率,起止時間、持續(xù)時間和能量分布也基本正確。為了分析相位譜的正確性,將圖1(e)的各分量單獨(dú)提出,并與理論值進(jìn)行比較,如圖2所示。圖2(a)、(b)、(c)為三分量信號的時域波形,圖2(d)、(e)、(f)為根據(jù)理論計算繪制出的真實時頻相位圖,圖2(g)、(h)、(i)為圖1(e)的各分量單獨(dú)的時頻相位圖。從圖2中各相應(yīng)分量逐一對比可以發(fā)現(xiàn),相位譜的瞬時相位十分準(zhǔn)確,頻率雖有所誤差但基本在正確頻率值附近擾動。

圖2 單分量信號的相位譜分析圖Fig.2 Phase spectrum analysis diagram of the single-component signal

3.2.2 相位解纏

由于我們所實測的相位值在±2π內(nèi),如果同一分量信號在兩個相鄰檢波器所測量的相位并不是位于同一周期內(nèi)的,那么此時的相位之差并不能反映真正物理意義上的相差,這種情況下必須要進(jìn)行相位解纏。

以第二分量信號為例,將圖1(e)和圖1(f)中的該分量信號提取出來,即以該分量信號能量中心為中心點(diǎn),截取相鄰兩道該分量信號相同大小的相位譜圖(如圖3所示)。從圖3中可以發(fā)現(xiàn),兩幅圖所共同覆蓋的時間內(nèi)任意一個時刻信號的相位差,除了實測相位之差外還應(yīng)該加上一個周期的補(bǔ)償。那么對于該分量信號而言,真實相位差等于實測相位差加2π。

3.3 結(jié)果分析

3.3.1 能量分布及群速度

利用能量譜和式(3)～(5)可以獲得各分量信號的頻率、振幅以及群速度。表2所示為平均頻率與能量分布。

表2 平均頻率與能量分布Table 2 Average frequency and energy distribution

圖3 相位解纏圖Fig.3 Phase unwrapping diagram

從表2中可以看出,各分量信號的頻率值非常準(zhǔn)確,但各分量的振幅值誤差較大,且絕大多數(shù)振幅值都偏小,這是由于EMD分解過程中能量泄漏所造成的。群速度計算結(jié)果及其誤差如表3所示。

表3 群速度計算結(jié)果及其誤差Table 3 Group velocity and error

從表3中可以看出,無論是時間還是群速度的計算結(jié)果都非常準(zhǔn)確。

3.3.2 相速度

從圖1(c)和圖1(d)中可以看出,各分量信號都經(jīng)過0.3 s,即無論是在12 m處還是18 m處,0.3 s時的信號包含各分量的信息。依據(jù)同步相差的原理可以選擇12 m處和18 m處在0.3 s時的瞬時相位計算各分量的相速度。

根據(jù)公式(6)可以計算出各分量的相速度,從相位譜中可以看出,低頻分量不用進(jìn)行相位解纏;中頻分量和高頻分量則需要在所求相差的基礎(chǔ)上增補(bǔ)一個周期2π和3個周期6π,才能求出真實相差。計算結(jié)果如表4所示。

表4 同步相差計算相速度Table 4 Phase velocity with synchronous phase difference

從圖1(d)中可以看出,各分量信號都經(jīng)過0.35 s,即在18 m處0.35 s時的信號包含各分量的信息。依據(jù)異步相差的原理可以選擇12 m處在0.3 s時的瞬時相位和18 m處在0.35 s時的瞬時相位計算各分量的相速度。

由于檢測時差選擇合適,低、中、高頻分量信號皆不用進(jìn)行相位解纏。依據(jù)公式(7)的計算結(jié)果如表5所示。

表5 異步相差計算相速度Table 5 Phase velocity with asynchronous phase difference

利用表3各分量的群速度求出該分量在相鄰檢波器間的傳播時差,并將該時差代入公式(7)即可求出該分量的相速度,計算結(jié)果見表6。

表6 校正相差計算相速度Table 6 Phase velocity with emendation phase difference

從表6中可以看出,低頻、中頻分量有效地避免了相位解纏。但高頻分量周期太小,且群速度計算結(jié)果誤差較大(見表3),故而實測相差與真實相差出現(xiàn)偏差。此時,實測相差需補(bǔ)上一個周期2π才能獲得真實相差。

4 結(jié) 論

運(yùn)用希爾伯特黃變換可以獲得信號的時頻-相位的精確分布信息,而異步相差和校正相差可以有效避免相位解纏,將上述方法用于進(jìn)行信號的頻散分析,不僅算法簡單而且結(jié)果準(zhǔn)確。但該方法部分細(xì)節(jié)依然有待改進(jìn):

(1)因為精度有限,不同時刻的相差不一定完全一致;而且瞬時頻率在真實頻率附近起伏,故不同時刻求取的相速度值會略有區(qū)別,所以確定某一分量的相速度時,可以對多個時間點(diǎn)的相速度進(jìn)行平均處理;

(2)如果某信號含有的單分量信號數(shù)目較多,直接采用EMD可能效果不會很理想,此時可以采用經(jīng)驗?zāi)B(tài)頻率分解[9,10](EMFD)的方法來獲得IMF,即首先對原信號進(jìn)行分頻段窄帶濾波,然后對每段濾波結(jié)果再單獨(dú)進(jìn)行希爾伯特黃變換。

[1] Dong-Soo Kim,Hyung-Choon Park.Determination of dispersive phase velocities for SASW method using harmonic wavelet transform[J].Soil Dynamics and Earthquake Engineering,2002,22(8):675-684.

[2] Jianghai Xia,Richard D Miller,Choon B Park,et al.Comparing shear-wave velocity profiles from MASW with borehole measurements in unconsolidated sediments[J].Journal of Environmental and Engineering Geophysics,2000,5(3):1-13.

[3] 科恩.時頻分析:理論與應(yīng)用[M].白居憲,譯.西安:西安交通大學(xué)出版社,1998.LeonCohen.Time-Frequency Analysis:Theory and Applications[M].Translated by BAI Ju-xian.Xi′an:Xi′an Jiaotong University Press,1998.(in Chinese)

[4] Norden E Huang,Zhaohua Wu.A Review on Hilbert-Huang Transform:Method an Its Applications to Geophysical Studies[J].Reviews of Geophysics,2008,46(2):1-23.

[5] Norden E Huang,Zheng Shen,Steven R Long.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proceedings of the Royal Society A:Mathematical Physical and Engineering Sciences,1998,454(1971):903-995.

[6] Norden Huang,Zhengshen,Steven R Long.A new view of nonlinear water waves:The Hilbert Spectrum[J].Annual Review of Fluid Mechanics,1999,31(1):417-457.

[7] Chau-Huei Chen,Cheng-Pling Li,Ta-Liang Teng.Surface-Wave Dispersion Measurements Using Hilbert-Huang Transform[J].Terrestrial,Atmospheric and OceanSciences,2002,13(2):171-184.

[8] 李白基,師潔珊,宋子安,等.地震面波的數(shù)字計算[J].地球物理學(xué)報,1977,20(4):283-298.LI Bai-ji,SHI Jie-shan,SONG Zi-an,et al.Digital Processing for Seismic Surface Wave Dispersion[J].Acta Geophysica Sinica,1977,20(4):283-298.(in Chinese)

[9] LI Jiang,XU Yixian.Analysis of dispersion of phase velocities using empirical mode frequency decomposition on SASW[C]//Proceedings of Near-Surface Geophysics and Geohazards.Chengdu:Science Press,2010:225-231.

[10] 程乾生,武連文.時間序列的經(jīng)驗?zāi)B(tài)頻率分解EMFD[J].數(shù)學(xué)的實踐與認(rèn)識,2005,36(5):151-153.CHENG Qian-sheng,WU Lian-wen.The Empirical Mode Frequency Decomposition for Time Series Analysis[J].Mathematics in Practice and Theory,2005,36(5):151-153.(in Chinese)