基于頻率-速度域多重信號分類的面波高分辨率頻散成像方法

伍敦仕,孫成禹,林美言

(中國石油大學(華東)地球科學與技術學院,山東青島266580)

基于頻率-速度域多重信號分類的面波高分辨率頻散成像方法

伍敦仕,孫成禹,林美言

(中國石油大學(華東)地球科學與技術學院,山東青島266580)

準確的相速度頻散圖像是主動源面波勘探方法反演近地表橫波速度的基礎。提出了一種基于頻率-速度域多重信號分類(multiple signal classification in frequency-velocity domain)的面波高分辨率頻散成像方法(簡稱fv-MUSIC方法)。該方法首先將傳統(tǒng)的頻率-波數(shù)域波束形成器改造成頻率-速度域形式,然后引入多重信號分類算法將空間譜相關矩陣分解為信號子空間和噪聲子空間兩個部分,最后利用噪聲子空間部分生成最終的面波頻散圖像。理論數(shù)據(jù)和實際數(shù)據(jù)應用結果表明,該方法能產生較高精度的相速度圖像,并且使用方便,計算效率高,尤其當接收排列較短時,該方法依然能保持較高的相速度分辨率,有利于提高主動源面波方法的橫向速度分辨能力。

多重信號分類;面波;相速度;高分辨率;頻散

在地震勘探領域,轉換波靜校正需要準確的近地表橫波速度[1-4]。在地震工程領域,一般需要地表到地下30m的平均橫波速度進行場地分類[5]。但長期以來獲得近地表橫波速度的手段并不多,主要有基于淺層SH折射波的旅行時解釋或層析反演和鉆井資料分析兩種,當存在下伏低速層等復雜情況時,前者會導致錯誤的解釋結果,后者由于施工成本較高,井的分布密度十分有限。雖然在地震勘探領域一般將面波視為噪聲將其去除[6-7],但這種做法正在改變,因為面波特有的頻散特征可以用來獲得近地表信息[1]。多道面波分析技術(multichannel analysis of surface waves,MASW)是20世紀90年代發(fā)展起來的一種基于面波的近地表無損探測技術,它通過反演面波相速度頻散曲線來獲得近地表橫波速度甚至品質因子。由于施工方便,不受常規(guī)折射波中存在的隱蔽層問題困擾,這項技術目前發(fā)展十分迅速[8-10]。

MASW技術的關鍵步驟在于獲得準確的面波相速度頻散圖像,進而根據(jù)能量峰值拾取頻散曲線和反演。F-K變換法是最早用于主動源面波頻散曲線成像的方法。將時-空域面波記錄進行二維傅里葉變換得到頻率-波數(shù)譜,通過拾取能量峰值和簡單的變換關系可得到相速度頻散曲線[11-14]。MCMECHAN等[15]最早提出了頻率-慢度變換法(ω-p變換法),并用其對勘探地震記錄中的面波進行了頻散成像。相移法是目前提取主動源面波頻散曲線應用最為廣泛的一項技術,已經成為許多面波處理商業(yè)軟件的標準模塊[16-18]。DAL等[19]曾對F-K變換法、ω-p變換法和相移法進行比較,結果顯示相移法的總體性能優(yōu)于前二者。LUO等[20]將高分辨率線性拉東變換(high-resolution linear Radon transform,HR-LRT)用于面波相速度頻散成像,其相速度分辨率遠高于F-K變換法、ω-p變換法以及相移法。

多重信號分類(MUSIC)是一種基于矩陣特征結構分解的超高分辨率空間譜估計方法[21],用于估計信號到達方向。該方法于20世紀90代末即被引入地震勘探領域,用于獲得高分辨率疊加速度譜[22-24]。相對于F-K等變換類方法,MUSIC方法在主動源面波頻散成像方面的應用研究還較少。IRANPOUR等[25]率先用MUSIC方法獲得Scholte波的高分辨率頻率-波數(shù)譜。隨后,WINSBORROW等[26]將其應用于海底Love波的頻散成像。BOIERO等[27]進一步將該方法應用于小尺度三維地震工程面波的頻散特征分析。MISBAH等[28]則將該方法應用于面波衰減系數(shù)的提取。但是上述研究均在頻率-波數(shù)域內進行,需要進行插值才能獲得對應的頻率-相速度頻散圖像。尤其需要指出的是,同為高分辨率方法,MUSIC法與HR-LRT法在獲取高分辨率面波頻散圖像方面存在哪些差異尚未見報道。本文提出了一種直接在頻率-速度域進行多重信號分類的高分辨率面波頻散成像方法(fv-MUSIC法),該方法無需額外的插值轉換即可獲得高分辨率的頻率-相速度圖像。此外,通過理論面波數(shù)據(jù)、工程面波數(shù)據(jù)和地震勘探中的含面波數(shù)據(jù)對比了fv-MUSIC法、相移法和HR-LRT法的頻散成像效果,驗證了fv-MUSIC方法的有效性和獨特優(yōu)勢。

1 fv -MUSIC方法原理

根據(jù)典型的主動源面波采集方式(圖1),在地表采用錘擊或落重法激發(fā)主動源面波(地震勘探中多采用爆炸震源),假設共有N個檢波器等間距均勻分布接收(也可以采用不等距分布),道距為dx,炮點與第一道距離為xS,則第i道的炮檢距為xi=xS+(i-1)dx。傳統(tǒng)的頻率域波束形成器形式為[29]:

(1)

式中:

(2)

為水平波數(shù)形式的平面波導向矢量(steering vector),H代表復共軛轉置;R(ω)為空間譜相關矩陣,

(3)

為了直接得到頻率-速度域高分辨率面波頻散圖像,利用水平波數(shù)與相速度之間的關系k=ω/v,我們將波束形成器形式改寫為:

(4)

其中,

(5)

為對應相速度v和角頻率ω的平面波導向矢量。

因為空間譜相關矩陣R(ω)是一艾爾米特對稱正定矩陣,所以它的特征向量通常與R(ω)正交,特征值一般都為正值[21]。將其逆矩陣進行特征分解,有:

(6)

式中:N為特征值個數(shù),等于接收道數(shù);λi為第i個特征值;vi為對應λi的特征向量。這些特征值與特征向量可以進一步分為對應信號和噪聲的兩個子空間部分,即:

(7)

其中較大的Ns個特征值對應信號子空間,其余特征值對應噪聲子空間。MUSIC方法只利用其中的噪

聲子空間部分進行最大分辨率空間譜估計,即:

(8)

將(8)式與平面波導向矢量(公式(5))相結合,得到最終的fv-MUSIC譜估計結果:

(9)

這里的vT表示角頻率ω的掃描相速度或實驗相速度。對每個頻率都采用一系列的相速度vT代入(9)式,即可獲得最終的頻散相速度圖像。

圖1 MASW面波采集

2 理論模型測試

首先合成一個簡單的頻散波記錄,用于檢驗在理想條件下fv-MUSIC方法的頻散成像效果,并將其與相移法、HR-LRT法進行對比。震源函數(shù)為主頻20Hz的雷克子波,假設理論相速度頻散關系為:

(10)

式中:v0代表高頻極限時的相速度,取300m/s;Δv為速度變化梯度,取500m/s;σ=30。由(10)式得到的理論相速度頻散曲線見圖2,通過(10)式和炮檢距xi可以得到每道的相位延遲,經過傅里葉反變換即得到最終的理論頻散波地震記錄(圖3)。該記錄接收道數(shù)100,道距2m,炮點與第一道距離10m,時間采樣間隔2ms。記錄中只有單純的頻散波,可以近似看成基階模式瑞雷波或勒夫波,代表橫波速度隨深度增加而遞增的典型地層條件下接收到的面波記錄。

分別采用相移法、fv-MUSIC法和HR-LRT法獲得的頻散圖像見圖4。對比發(fā)現(xiàn),fv-MUSIC法分辨率最高,幾乎達到和理論頻散曲線相同的效果;HR-LRT法效果次之;相移法效果不夠理想,除了相速度能量軸最粗之外,還存在許多能量較弱的虛假相速度。在處理實際面波資料時,有可能將這些虛假相速度能量軸誤判為某階面波模式,得到錯誤的反演結果。一種好的頻散成像方法應該在突出真實面波相速度的同時盡可能壓制上述假象。在計算效率方面,由于相移法是在頻率域進行簡單的傾斜疊加,所以速度最快,用時約7s。fv-MUSIC法由于對每個頻率都進行矩陣特征分解,所以計算量增加,用時約35s。HR-LRT法由于在每個頻率都需要進行反演問題的求解,計算效率最低,用時約338s,并且正則化參數(shù)的選擇直接影響HR-LRT法頻散成像的效果和計算效率,需要通過多次試驗來獲得滿意的頻散圖像。

圖2 理論頻散曲線

將圖3所示頻散波地震記錄加入一定程度的均勻分布隨機噪聲,使其信噪比為1,得到如圖5a所示含噪記錄,用以考察3種頻散成像方法的抗噪性。圖5b 到圖5d為3種方法頻散成像結果,可見由于噪聲的影響,3種方法高頻部分(>50Hz)的相速度聚焦都不成功,低頻部分(<10Hz)也出現(xiàn)了明顯的抖動。也就是說,從成像頻帶寬度來看,3種方法的抗噪性基本一致。因為MASW方法基于地層橫向速度不變或緩變這一假設條件,所以得到的相速度頻散曲線反映的是接收排列范圍內地下的平均效應。為了得到二維近地表橫波速度剖面,一般采用滾動式面波激發(fā)采集方式得到一系列頻散圖像,提取頻散曲線,反演得到一組一維橫波速度,通過橫向插值構建出二維橫波速度分布圖。所以,如果能用更少的接收道數(shù)得到合理的頻散圖像,將有利于提高MASW方法的橫向分辨能力。以下選擇圖5a中第41至第50道(共10道)記錄進行頻散分析,比較3種方法的成像效果(圖6)。

圖3 簡單頻散波地震記錄

對比圖6b至圖6d可見,記錄道數(shù)減少后,相移法得到的頻散圖像的相速度分辨率大幅下降,HR-LRT法得到的頻散圖像的相速度分辨率下降也十分明顯。這是因為,作為傾斜疊加類算法,相移法和HR-LRT法的分辨率直接受接收排列長度影響,接收排列越長,分辨率越高,反之越低。fv-MUSIC法的相速度分辨能力依然很高(圖6d),只是由于記錄道數(shù)減少,其相速度能量軸的光滑度相比圖5d有所下降,速度能量軸出現(xiàn)了抖動。對比圖6b,圖6c和圖6d可見,fv-MUSIC法成像效果優(yōu)于相移法和HR-LRT法。

圖4 理論模型3種方法頻散成像效果a 相移法; b HR-LRT法; c fv -MUSIC法

圖5 頻散波記錄信噪比為1時3種方法頻散成像效果a 含噪聲頻散波記錄; b 相移法頻散成像效果; c HR-LRT法頻散成像效果; d fv -MUSIC法頻散成像效果

圖6 中間10道記錄及3種方法頻散成像效果a 第41至第50道信號; b 相移法頻散成像效果; c HR-LRT法頻散成像效果; d fv -MUSIC法頻散成像效果

3 實際資料應用試驗

3.1 工程面波記錄

選擇一炮采集于美國Texas A & M 大學(TAMU)國家地學技術試驗場的MASW工程面波記錄(圖7),62道接收,道距0.6096m,時間采樣間隔0.7813ms,炮點與第一接收道距離3.048m。該數(shù)據(jù)由美國工程地球物理協(xié)會于2011年公開向外界發(fā)布,用于近地表面波科研人員測試處理與反演算法的研究[30]。

圖7 MASW工程面波記錄

由圖7可見,該試驗場淺地表存在橫向速度變化,記錄第10道、第22道附近都能看到明顯的背向散射面波,這會對后續(xù)頻散曲線的提取產生一定影響。圖8為3種頻散成像方法處理得到的頻散圖像,可以看到基階模式與第一高階模式瑞雷波能量很強,特別是第一高階模式瑞雷波從50Hz一直延續(xù)到80Hz,非常適合于后續(xù)基階與高階模式的瑞雷波聯(lián)合反演。從3種方法的成像品質來看,相移法分辨率最低,并且頻散圖像中存在大量背景噪聲干擾;fv-MUSIC法和HR-LRT法得到的頻散圖像分辨率不相上下,都明顯高于相移法,并且背景噪聲得到了較好壓制,主要瑞雷波模式得以很好凸顯。HR-LRT法得到的頻散圖像雖然分辨率很高,但對有效面波模式造成了壓制(圖8b),第一高階模式瑞雷波的能量強度與連續(xù)性明顯不如相移法和fv-MUSIC法。

圖8 工程面波記錄3種方法頻散成像效果a 相移法; b HR-LRT法; c fv -MUSIC法

3.2 地震勘探中的含面波記錄

選擇中國西部某工區(qū)一個含面波地震勘探反射波記錄進一步對比了3種頻散成像方法的效果(圖9)。該記錄接收道距不均勻,但總體保持在10m左右。記錄中面波十分發(fā)育,線性同相軸連續(xù)并清晰可辨,沒有空間假頻現(xiàn)象,多模式特征十分明顯,在2.2s和3.1s附近能看到明顯的反射波同相軸。記錄中面波的線性同相軸特征和反射波同相軸的雙曲線特征都表明,該工區(qū)地表比較平坦,地下無明顯傾斜構造,可以看作水平層狀介質，適合采用面波法進行處理。

圖9 含面波地震勘探反射波記錄

圖10給出了3種頻散成像方法得到的面波頻散

圖10 實際資料3種方法頻散成像效果a 相移法; b HR-LRT法; c fv -MUSIC法

圖像。通過比較可以發(fā)現(xiàn),相移法得到的圖像中背景噪聲比較嚴重,但可以清晰識別3～10Hz之間的基階模式瑞雷波和8～17Hz之間的第一高階模式瑞雷波。實際上,第二高階模式與第三高階模式瑞雷波也能識別,但是能量較弱,受噪聲污染比較嚴重。fv-MUSIC法得到的頻散圖像在低于3Hz的部分成像效果看起來不如其它兩種方法,但是理論模型測試(參見第2節(jié))表明,頻率太低的部分因為噪聲的影響相速度可靠性較低,不應該進行提取并參與到反演中。因此,3種方法都能夠可靠拾取的基階模式相速度頻率范圍是3～10Hz。在圖10c中,fv-MUSIC方法對第一高階模式可以識別到20Hz,而其它兩種方法只能識別到17Hz左右(圖10a,圖10b),這進一步說明3種方法中fv-MUSIC方法存在優(yōu)勢。

4 結束語

本文基于多重信號分類算法提出了一種直接在頻率-速度域進行面波高分辨率頻散成像的方法,即fv-MUSIC方法。通過與目前廣泛使用的相移法和HR-LRT法進行比較,獲得以下認識:

1) 對于不含噪聲的理論面波數(shù)據(jù),fv-MUSIC方法的頻散成像分辨率可以接近理論頻散曲線的程度；

2) 存在隨機噪聲干擾時,fv-MUSIC方法的成像頻帶范圍和相移法、HR-LRT法基本一致,所以,fv-MUSIC方法的抗噪性與相移法和HR-LRT法相當；

3) 采用短排列接收時,相移法和HR-LRT法的相速度分辨率顯著下降,而fv-MUSIC法依然能保持較高分辨率；

4) 在對實際面波資料進行頻散成像時,fv-MUSIC方法能在壓制背景干擾的同時不降低高階模式面波的成像質量。

需要指出的是,F-K變換法和HR-LRT法由于存在正反變換,在面波提取和模式分離方面依然具有很重要的作用。而fv-MUSIC法屬于單純的空間譜估計技術,在這些方面的作用有限。

[1] 孟小紅,郭良輝.利用地震瑞利波速度反演求取P-SV波橫波靜校正量[J].石油地球物理勘探,2007,42(4):448-453 MENG X H,GUO L H.Using velocity inversion of seismic Rayleigh wave to compute s-wave statics of P-SV wave[J].Oil Geophysical Prospecting,2007,42(4):448-453

[2] MARI J L.Estimation of static corrections for shear-wave profiling using the dispersion properties of Love waves[J].Geophysics,1984,49(8):1169-1179

[3] BOIERO D,MARSDEN P,ESAULOV V,et al.Building a near-surface velocity model in the South Ghadames basin:surface-wave inversion to solve complex statics[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:1811-1815

[4] PIATTI C,SOCCO L V,BOIERO D,et al.Constrained 1D joint inversion of seismic surface waves and P-refraction travel times[J].Geophysical Prospecting,2013,61(S1):77-93

[5] 夏江海,高玲利,潘雨迪,等.高頻面波方法的若干新進展[J].地球物理學報,2015,58(8):2591-2605 XIA J H,GAO L L,PAN Y D,et al.New findings in high-frequency surface wave method[J].Chinese Journal of Geophysics,2015,58(8):2591-2605

[6] 徐穎,劉晨,呂秋玲,等.多域組合去噪技術在塔中奧陶系低信噪比資料處理中的應用[J].石油物探,2015,54(2):172-179 XU Y,LIU C,LV Q L,et al.Application of multi-domain composite denoising technology for the processing of Ordo-vician low SNR seismic data in Tazhong area[J].Geophysical Prospecting for Petroleum,2015,54(2):172-179

[7] 王德營,李振春,王姣.S域時頻濾波分析與改進[J].石油物探,2015,54(4):396-403 WANG D Y,LI Z C,WANG J.The analysis and improvement on time-frequency filtering in S-transform domain[J].Geophysical Prospecting for Petroleum,2015,54(4):396-403

[8] PARK C B,MILLER R D,XIA J.Multichannel analysis of surface waves[J].Geophysics,1999,64(3):800-808

[9] SOCCO L V,FOTI S,BOIERO D.Surface-wave analysis for building near-surface velocity models—established approaches and new perspectives[J].Geophysics,2010,75(5):A83-A102

[10] FOTI S,PAROLAI S,ALBARELLO D,et al.Application of surface-wave methods for seismic site characterization[J].Surveys in Geophysics,2011,32(6):777-825

[11] GABRIELS P,SNIEDER R,NOLET G.In situ measurements of shear-wave velocity in sediments with higher-mode Rayleigh waves[J].Geophysical Prospecting,1987,35(2):187-196

[12] TSELENTIS G,DELIS G.Rapid assessment of S-wave profiles from the inversion of multichannel surface wave dispersion data[J].Annali Di Geofisica,1998,41(1):1-15

[13] FOTI S,LANCELLOTTA R,SAMBUELLI L,et al.Notes on f-k analysis of surface waves[J].Annali Di Geofisica,2000,43(6):1199-1209

[14] 李子偉,劉學偉.空間假頻對頻率-波數(shù)域瑞利面波頻散曲線反演的影響[J].石油地球物理勘探,2013,48(3):395-402 LI Z W,LIU X W.Effects of spatial aliasing on frequency-wave number domain inversion of Rayleigh-wave dispersion curve[J].Oil Geophysical Prospecting,2013,48(3):395-402

[15] MCMECHAN G A,YEDLIN M J.Analysis of dispersive waves by wave field transformation[J].Geophysics,1981,46(6):869-874

[16] PARK C B,MILLER R D,XIA J.Imaging dispersion curves of surface waves on multi-channel records[J].Expanded Abstracts of 68thAnnual Internat SEG Mtg,1998:1377-1380

[17] RYDEN N,PARK C B,ULRIKSEN P,et al.Multimodal approach to seismic pavement testing[J].Journal of Geotechnical and Geoenvironmental Engineering,2004,130(6):636-645

[18] OBANDO E A,PARK C B,RYDEN N,et al.Phase-scanning approach to correct time-shift inaccuracies in the surface-wave walk-away method[J].Soil Dynamics and Earthquake Engineering,2010,30(12):1528-1539

[19] DAL M G,PIPAN M,FORTE E,et al.Determination of Rayleigh-wave dispersion curves for near surface applications in unconsolidated sediments[J].Expanded Abstracts of 73rdAnnual Internat SEG Mtg,2003:1247-1250

[20] LUO Y,XIA J,MILLER R D,et al.Rayleigh-wave dispersive energy imaging using a high-resolution linear Radon transform[J].Pure and Applied Geophysics,2008,165(5):903-922

[21] SCHMIDT R O.Multiple emitter location and signal parameter estimation[J].IEEE Transactions on Antennas and Propagation,1986,34(3):276-280

[22] BIONDI B L,KOSTOV C.High-resolution velocity spectra using eigenstructure methods[J].Geophysics,1989,54(7):832-342

[23] KIRLIN R L.The relationship between semblance and eigenstructure velocity estimators[J].Geophysics,1992,57(8):1027-1033

[24] BARROS T,LOPES R,TYGEL M.Implementation aspects of eigendecomposition-based high-resolution velocity spectra[J].Geophysical Prospecting,2015,63(1):99-115

[25] IRANPOUR K E,MUYZERT E,GRION S.Local velocity analysis by parametric wavenumber estimation in seismic (fk-MUSIC)[J].Expanded Abstracts of 64thEAGE Annual International Conference and Exhibition,2002:171-174

[26] WINSBORROW G,HUWS D G,MUYZERT E.Acquisition and inversion of Love wave data to measure the lateral variability of geoacoustic properties of marine sediments[J].Journal of Applied Geophysics,2003,54(1/2):71-84

[27] BOIERO D,BERGAMO P,REGE R B,et al.Estimating surface-wave dispersion curves from 3D seismic acquisition schemes:part 1,1D models[J].Geophysics,2011,76(6):G85-G93

[28] MISBAH A S,STROBBIA C L.Joint estimation of modal attenuation and velocity from multichannel surface wave data[J].Geophysics,2014,79(3):EN25-EN38

[29] JOHNSON D H,DUDGEON D E.Array signal processing[M].Upper Saddle River:PTRPrentice-Hall Inc,1993:349-392

[30] TRAN K T,HILTUNEN D R.One-dimensional inversion of full waveforms using a genetic algorithm[J].Journal of Environmental & Engineering Geophysics,2012,17(4):197-213

(編輯：戴春秋)

High resolution dispersion imaging of surface waves based on multiple signal classification in frequency-velocity domain

WU Dunshi,SUN Chengyu,LIN Meiyan

(SchoolofGeosciences,ChinaUniversityofPetroleum,Qingdao266580,China)

Accurate phase velocity dispersion image is the basis of surface wave inversion for near-surface shear wave velocity.We present a high resolution dispersion imaging method of surface wave named multiple signal classification in frequency-velocity domain (fv-MUSIC) based on multiple signal classification in frequency velocity domain.Firstly,a beam-former is constructed in frequency velocity domain instead of the conventional frequency wavenumber domain.Secondly,the MUSIC algorithm is introduced to decompose the spatial-spectral cross-correlation matrix into signal and noise subspaces.At last,the ultimate surface wave dispersion image is produced using the noise subspace.Theoretical and real surface wave data processing cases indicate that it has the advantages of convenient operation,less calculation time and high resolution phase velocity images.Especially when the length of receiving array is shorter,thefv-MUSIC method still remains higher phase velocity resolution.It helps to improve the lateral shear wave velocity resolution of the active source surface wave method.

multiple signal classification,surface wave,phase velocity,high resolution,dispersion

2016-02-01;改回日期：2016-03-19。

伍敦仕(1983—),男,博士在讀,主要從事地震面波的正演、處理及反演研究。

國家自然科學基金項目(41374123)資助。

P631

A

1000-1441(2017)01-0141-09

10.3969/j.issn.1000-1441.2017.01.016

This research is financially supported by the National Natural Science Foundation of China (Grant No.41374123).