蔣 禮

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

1 引 言

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

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

2 原理簡(jiǎn)介

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

2.1 希爾伯特黃變換

2.1.1 瞬時(shí)頻率

[3]

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

因此瞬時(shí)頻率為

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

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

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

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

2.1.3 希爾伯特譜

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

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

2.2 頻散分析

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

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

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

2.2.1 群速度

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

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

而該單分量信號(hào)的傅里葉頻率即平均頻率為

該單分量信號(hào)的振幅計(jì)算公式為

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

2.2.2 相速度

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

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

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

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

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

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

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

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

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

校正相差為

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

3 仿真實(shí)例

3.1 模型建立

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

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

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

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

3.2 希爾伯特譜分析

3.2.1 相位譜和能量譜

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

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

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

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

3.2.2 相位解纏

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

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

3.3 結(jié)果分析

3.3.1 能量分布及群速度

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

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

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

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

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

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

3.3.2 相速度

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

(2)如果某信號(hào)含有的單分量信號(hào)數(shù)目較多,直接采用EMD可能效果不會(huì)很理想,此時(shí)可以采用經(jīng)驗(yàn)?zāi)B(tài)頻率分解[9,10](EMFD)的方法來獲得IMF,即首先對(duì)原信號(hào)進(jìn)行分頻段窄帶濾波,然后對(duì)每段濾波結(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] 科恩.時(shí)頻分析:理論與應(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ì)算[J].地球物理學(xué)報(bào),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] 程乾生,武連文.時(shí)間序列的經(jīng)驗(yàn)?zāi)B(tài)頻率分解EMFD[J].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí),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)