周顯華,沈金松,李亞曦,侯 桐,高 鑫

(中國石油大學(xué)(北京)地球物理學(xué)院,北京 102249)

陣列聲波測井儀接收的波形受井孔周圍巖石性質(zhì)的影響,其慢度和頻散特征等均與地層性質(zhì)密切相關(guān)[1]。準(zhǔn)確、有效地提取模式波中的頻散信息,能夠更好地獲取與地層性質(zhì)相關(guān)的信息。然而,實測波形數(shù)據(jù)受噪聲和模式波相互干擾等影響并不嚴格滿足數(shù)學(xué)模型假設(shè),相比于參數(shù)化頻散分析算法對理論模型誤差的敏感性,非參數(shù)頻散分析算法無需假設(shè)信號模型且具有更好的魯棒性,更適合處理復(fù)雜的實測波形數(shù)據(jù)[2]。

1898年,SCHUSTER[3]利用傅里葉變換算法進行譜分析時提出了周期圖的概念,但該算法存在主瓣能量泄露和方差較大的問題。NOLTE等[4]基于傅里葉變換理論提出了處理陣列譜相干函數(shù)的加權(quán)頻譜相干法(WSS),采用多點加權(quán)策略提高了算法精度及抗噪性,但該方法不適用于處理單一頻率下多個模式波共存的波形數(shù)據(jù)。前人對周期圖進行了大量研究,其中包括帶通濾波器組的處理算法,如最小方差無失真波束形成算法(Capon),通過設(shè)計數(shù)據(jù)相關(guān)的濾波器進行譜估計[5]。LI等[6]提出了一種自適應(yīng)有限脈沖濾波器算法,也稱為振幅相位估計法(APES),最初用于目標(biāo)范圍的特征估計與合成孔徑雷達成像。與最小方差無失真波束形成算法相比,振幅相位估計法算法對幅度估計更加準(zhǔn)確。曹正良等[7]利用模擬和實測陣列聲波測井?dāng)?shù)據(jù)比較了普羅尼(Prony)方法、同態(tài)處理方法和加權(quán)頻譜相干法3種頻散分析方法。王瑞甲等[8]將矩陣束算法與加權(quán)頻譜相干法相結(jié)合,給出了一種多道波形信號的頻散分析方案,壓制了加權(quán)頻譜相干法處理中出現(xiàn)的周期性偽解。劉航等[9-10]通過縱波的頻散分析提出縱波頻散譜處理方法,并形成基于頻散譜的碳酸鹽巖儲層評價方法。

陣列聲波測井波場分離技術(shù)用于分離全波列中的不同模式波[11]。宮昊等[12]用慢度-時間相干法與線性預(yù)測濾波法的組合從全波列中提取較弱的反射波,有效地解決了反射聲波測井中成像結(jié)果的多解性和可靠性差的問題。實測數(shù)據(jù)處理中,由于模式波相互干涉及噪聲影響等,在慢度-頻譜圖中存在模式波慢度-頻率相關(guān)曲線不連續(xù)或無法有效分辨等問題。因此,有必要將時、頻域波場分離與頻散提取方法相結(jié)合以實現(xiàn)波形數(shù)據(jù)的頻散分析。首先,采用長短時方差比值法[13]準(zhǔn)確提取初至,然后,將其作為時域濾波窗的起止時刻,并實施時域波場分離。

雖然加權(quán)頻譜相干法、振幅相位估計法和最小方差無失真波束形成法等算法的性能在以往文獻[7,14-15]中已有闡述,但幾種算法的應(yīng)用效果以及性能的綜合對比尚不全面,實測數(shù)據(jù)的頻散信息提取效果也不理想。本文利用模擬數(shù)據(jù)與實測數(shù)據(jù)詳細比較了傅里葉變換法(FTM)、加權(quán)頻譜相干算法、振幅相位估計算法和最小方差無失真波束形成算法及其前后向算法等共6種頻散信息提取算法的性能及應(yīng)用效果,并認識到將基于長短時方差比值的時、頻域濾波的波場分離算法與非參數(shù)頻散估計算法相結(jié)合,可以顯著改善頻散信息的提取效果。

1 頻散信息提取方法的對比分析

考慮由N個等間距d接收器組成的聲波陣列,它接收的波形信號是平面波與噪聲線性疊加后的信號。若以第一個陣列聲波接收器為參考,則頻率域波形信號的解析解可以表示為[16]:

(1)

式中:n=1,2,…,N表示接收器的編號;αp和sp分別為第p個模式波的復(fù)幅度及慢度;ω=2πf是圓頻率,f為頻率;en(ω)表示接收信號中的噪聲以及第n個接收器中模式波之間的相互干涉。

1.1 頻散提取算法的基本關(guān)系

振幅相位估計算法(APES)是根據(jù)噪聲協(xié)方差矩陣的特征值分布與大小自適應(yīng)地設(shè)計相應(yīng)有限脈沖濾波器,使得慢度為s的信號無失真地通過,并對其它信號及噪聲進行壓制,振幅相位估計算法及其前后向(FB)算法估計所得復(fù)振幅分別為[6]:

(2)

(3)

(4)

(5)

加權(quán)頻譜相干法包括頻率域相關(guān)及加權(quán)兩個步驟[4]。將全波列波形數(shù)據(jù)dn(t)變換到頻率域{|Dn(ω)|}(n代表接收器編號),第n個接收器與第n+1個接收器之間的關(guān)系為:

Dn+1(ω)=Dn(ω)exp(-a0d)exp(-iωsd)

(6)

式中:a0表示頻率ω下模式波的衰減系數(shù);s表示其模式波慢度;d為接收器間距。(6)式中,第一個指數(shù)部分代表衰減,第二個指數(shù)部分代表波的傳播。若定義向量D(ω)與X(ω,s)為:

(7)

(8)

式中:X(ω,s)中的元素代表各接收器相對于第一個接收器的相移。則頻率慢度相關(guān)歸一化函數(shù)為:

(9)

式中:符號*表示共軛轉(zhuǎn)置。對(9)式進行加權(quán)得到加權(quán)相干系數(shù),即

(10)

式中:W(ωj)為權(quán)函數(shù);l為加權(quán)的波譜點半徑?？墒褂酶咚购瘮?shù)作權(quán)函數(shù),在被加權(quán)的頻率點處具有最大的權(quán)重,即

(11)

式中:σ=NwΔω,Nw表示參與加權(quán)的數(shù)據(jù)點數(shù),Δω代表頻率間隔。僅單點加權(quán)(Nω=1)時,加權(quán)頻譜相干算法將退化為(9)式。

傅里葉變換法(FTM)利用傅里葉變換算法得到頻率域波形信號,消除由于傳播帶來的相位差異即可得到近似的復(fù)幅度估計。在實際數(shù)據(jù)處理中,(9)式分子部分即為FTM。

1.2 提取方法性能對比分析

1.2.1 保幅性及分辨率

根據(jù)陣列聲波波形信號的解析解(1)式,設(shè)計慢度分別為50μs/ft,60μs/ft,80μs/ft 3個模式波的模型(1ft≈30.48cm)。模擬數(shù)據(jù)處理結(jié)果如圖1a所示,黑色圓點為模式波實際慢度及振幅,星標(biāo)(曲線極值點)為算法估計的慢度及幅度;振幅相位估計算法、前后向振幅相位估計算法以及前后向最小方差無失真波束形成算法估計的80μs/ft模式波幅度相近;加權(quán)頻譜相干算法及傅里葉變換算法不能有效識別鄰近較弱模式波,弱模式主瓣被強模式旁瓣所淹沒,兩者均是基于傅里葉變換理論來估計信號所包含的慢度成分,它是周期函數(shù)截斷的結(jié)果,造成截斷邊界處存在能量泄漏,導(dǎo)致估計結(jié)果失真、分辨率下降。振幅相位估計類算法及最小方差無失真波束形成類算法以目標(biāo)慢度信號保真為前提設(shè)計的濾波器,能準(zhǔn)確地估計模式波的慢度。振幅相位估計算法及前后向振幅相位估計算法均可準(zhǔn)確估計信號幅度和相位。最小方差無失真波束形成算法估計的信號幅度普遍偏低,前后向最小方差無失真波束形成算法比僅利用前向數(shù)據(jù)的最小方差無失真波束形成算法幅度估計精度更高。

對遼河某地實測數(shù)據(jù)進行頻散分析,6種算法的實際數(shù)據(jù)處理(2.41kHz處)的最大分辨率對比結(jié)果如圖1b所示。此處僅進行分辨率對比,所以將所有掃描慢度點的幅度值均除以幅度最大值進行歸一化。最小方差無失真波束形成類頻散分析算法無旁瓣、能量集中、分辨率最佳,振幅相位估計類算法分辨率次之,傅里葉變換法及加權(quán)頻譜相干算法分辨率最低。

1.2.2 抗噪能力

聲波測井實測數(shù)據(jù)不可避免地含有噪聲。在單頻(8kHz)單模式波(80μs/ft)波形解析解(1)式信號中加入隨機高斯白噪聲,以慢度估計的相對誤差衡量算法抗噪性。模擬數(shù)據(jù)接收器總數(shù)設(shè)為13個,最小方差無失真波束形成算法、振幅相位估計算法及其前后向算法的濾波器權(quán)系數(shù)個數(shù)均為6,曲線每點均為1000次重復(fù)實驗的相對誤差平均結(jié)果。圖2顯示了非參數(shù)算法的抗噪能力對比結(jié)果。由圖2可見,最小方差無失真波束形成算法和振幅相位估計算法對應(yīng)的前后向算法抗噪性能得到提升;傅里葉變換算法與單點加權(quán)的加權(quán)頻譜相干算法具有極強的抗噪性,兩者慢度估計相對誤差隨信噪比變化曲線幾乎重合,抗噪性能相當(dāng);多頻點數(shù)據(jù)加權(quán)的頻譜相干算法增大了數(shù)據(jù)量,有效降低了噪聲的影響。

圖2 非參數(shù)算法抗噪能力對比結(jié)果

1.3 提取方法應(yīng)用效果

利用遼河某地實測單極子聲波波形數(shù)據(jù)對比6種頻散提取方法的應(yīng)用效果。每個波列共有672個時間采樣點,采樣間隔為16μs,8個接收器,振幅相位估計算法與最小方差無失真波束形成算法以及它們的前后向算法濾波器權(quán)系數(shù)個數(shù)均為4。處理結(jié)果如圖3所示。由圖3可見,傅里葉變換算法(圖3c)與加權(quán)頻譜相干算法(圖3f)的處理結(jié)果中能量較強的旁瓣會對頻散信息準(zhǔn)確提取造成干擾,但其模式波的慢度-頻率曲線連續(xù)性好,能量分布較為均衡。振幅相位估計算法的處理結(jié)果(圖3a)受噪聲等影響較大,存在如白圈中心處的能量局部異常點。前后向算法使得能量更加聚集,多模式分析方法中的前后向最小方差無失真波束形成算法(圖3e)與前后向振幅相位估計算法(圖3d)處理效果相近,單模式分析方法中的加權(quán)頻譜相干算法處理效果更好(圖3f)。6種頻散提取方法處理結(jié)果中模式波及偽解特征相似,但均未有效識別出縱波以及偽瑞利波。

圖3 頻散提取方法實測數(shù)據(jù)應(yīng)用效果對比a 振幅相位估計算法; b 最小方差無失真波束形成算法; c 傅里葉變換算法; d 前后向振幅相位估計算法; e 前后向最小方差無失真波束形成算法; f 加權(quán)頻譜相干算法

表1為數(shù)值模擬計算采用的井孔流體和地層參數(shù),所用聲源函數(shù)選擇主頻為12kHz的Ricker子波,最小源距為3.66m,13個換能接收器等間隔(0.1524m)排列。

表1 正演模型參數(shù)

衰減模型為[18]:

(12)

式中:Q為60,ωref為12kHz,頻散提取方法模擬數(shù)據(jù)應(yīng)用效果對比如圖4所示。結(jié)果中虛線分別為井內(nèi)流體慢度(203.2μs/ft)以及地層縱、橫波理論頻散曲線。

由圖4可見,傅里葉變換算法(圖4c)及加權(quán)頻譜相干算法(圖4f)由于強信號能量的泄露會嚴重干擾較弱模式的識別,所以更適用于單一模式波識別。即使在不含噪聲的情況下,由于算法分辨率及不同模式波之間的相互干涉等影響,6種非參數(shù)頻散信息提取方法也都不能同時且完整地提取所有模式波。

圖4 頻散提取方法模擬數(shù)據(jù)應(yīng)用效果對比a 振幅相位估計算法; b 最小方差無失真波束形成算法; c 傅里葉變換算法; d 前后向振幅相位估計算法; e 前后向最小方差無失真波束形成算法; f 加權(quán)頻譜相干算法

針對實際應(yīng)用中頻散信息提取效果變差的問題,我們認為,首先,應(yīng)進行基于長短時方差比值的時域濾波與頻域濾波相結(jié)合實現(xiàn)波場分離。然后,再進行頻散提取處理。波場分離將多種模式波疊加的信號分解成單一模式,降低模式波間的相互干涉及噪聲等干擾對處理效果的影響;同時,也能使單模式頻散分析方法較為準(zhǔn)確地提取縱波、橫波及偽瑞利波的頻散信息。

2 模式波分離方法

不同模式波出現(xiàn)時會表現(xiàn)出明顯的能量變化。定義反映瞬態(tài)能量變化和整體能量變化的“短時方差”及“長時方差”,兩者的比值可以突出初至到達時的能量變化,從而檢測縱、橫波波至[19](圖5)。設(shè)信號f(t)在第n點處的短時窗能量與長時窗能量的比值為[13]:

圖5 長短時方差提取波至

(13)

式中:l為短時窗的長度。

聲波全波列中,各模式波傳播速度不同,同時儀器較大源距使得時域波形分離較好,利用長短時方差比值法提取初至從而完成時域波場分離,再結(jié)合頻域分離即可較好地分離目標(biāo)模式波。波場分離使得單模式頻散分析方法能夠識別出被較強模式旁瓣淹沒的較弱模式波,圖6a為圖4對應(yīng)的地層中單極子激發(fā)所得完整波形數(shù)據(jù)縱波分離后的頻散分析結(jié)果,即

圖4f縱波分離后的結(jié)果,白色虛線為縱波理論頻散曲線。圖6b為與圖3相同的實測數(shù)據(jù)體的處理結(jié)果,即圖3f縱波分離后的結(jié)果,白色虛線為該深度的聲波時差參考值(54.8μs/ft)。處理結(jié)果對比表明,基于長短時方差比值提取模式波初至與時域濾波結(jié)合,再進行頻域濾波的波場分離不改變目標(biāo)模式波慢度及頻譜分布特征。

圖6 縱波分離后加權(quán)頻譜相干算法處理結(jié)果對比a 模擬數(shù)據(jù)分離后; b 實測數(shù)據(jù)分離后

3 頻散信息提取

基于頻散分析方法的性能及應(yīng)用效果的綜合對比可知,最小方差無失真波束形成類算法抗噪能力強、模式波慢度估計精度及分辨率高,受旁瓣影響極小,且能夠有效識別同一頻率下多個強弱不同的模式波。前后向最小方差無失真波束形成算法相對僅利用前向最小方差無失真波束形成算法綜合性能更好,所以選擇前后向最小方差無失真波束形成算法提取頻散信息較為理想。

實測數(shù)據(jù)處理時,若斯通利波能量軸連續(xù)、峰值強,可直接提取;否則,利用長短時方差比值法提取縱、橫以及斯通利波初至分離縱橫波,然后分別提取頻散信息;最后提取斯通利波及橫波慢度范圍內(nèi)能量較強的多階偽瑞利波(p-R),得到所有模式波頻散信息如圖7a所示。從圖7a中看到,提取的頻散曲線與細線的理論頻散曲線吻合較好。圖7b給出了圖3中使用的遼河某地實測數(shù)據(jù)得到的頻散曲線,圖中看到縱波、橫波和斯通利波的頻散曲線均連續(xù)性好,圖3 中不能提取的縱波模式也得到了較好的提取結(jié)果。

圖7 頻散信息提取結(jié)果a 模擬數(shù)據(jù); b 實測數(shù)據(jù)

4 結(jié)論

基于非參數(shù)估計算法能夠較好地處理模式波數(shù)量未知且慢度頻散規(guī)律復(fù)雜的實測數(shù)據(jù),利用模擬數(shù)據(jù)及實測數(shù)據(jù)綜合對比了傅里葉變換算法、加權(quán)頻譜相干算法、振幅相位估計算法、前后向振幅相位估計算法、最小方差無失真波束形成算法和前后向最小方差無失真波束形成算法共6種非參數(shù)頻散信息提取算法的保幅性、分辨率、抗噪性以及實際應(yīng)用效果,并與波場分離相結(jié)合提取了不同模式波的頻散信息,得到以下認識:

1) 加權(quán)頻譜相干算法和傅里葉變換算法具有良好的抗噪能力,但旁瓣大、保幅能力差且分辨率低;振幅相位估計類算法具有較高分辨率,振幅、相位估計優(yōu)勢明顯,但抗噪能力弱;最小方差無失真波束形成算法振幅估計值偏低,但其模式波慢度估計精度高,同時具有良好的分辨率及抗噪能力。與僅使用前向數(shù)據(jù)的算法相比,前后向算法的抗噪能力有所提高;前后向最小方差無失真波束形成算法的幅值估計精度比最小方差無失真波束形成算法的估計精度高近一倍。因此,建議在實測數(shù)據(jù)的速度頻散信息分析中選用前后向最小方差無失真波束形成算法。

2) 在處理實測數(shù)據(jù)時,振幅相位估計算法和最小方差無失真波束形成算法及其前后向4種多模式分析方法處理結(jié)果的慢度-頻譜曲線連續(xù)性差,峰值能量變化大;單模式分析方法僅能分辨能量較強的模式波,難以分辨慢度間隔較小的相鄰模式且旁瓣能量強;由于噪聲、模式波的相互干涉及算法分辨率等影響,6種頻散信息提取方法均難以同時有效且完整地分辨出所有模式波。

3) 基于長短時方差比值提取初至并結(jié)合頻域濾波的波場分離方法不改變模式波慢度及頻散特征,與頻散分析方法組合使用可以降低噪聲及模式波相互干涉的影響,提升頻散分析效果。而且,波場分離后用單模式頻散分析方法也能較好提取縱波和斯通利波等模式波的頻散信息。

由于偽瑞利波與斯通利波和橫波在時頻域均有混疊,基于長短時方差比值的時域與頻域相結(jié)合的波場分離方法仍然不能有效分離偽瑞利波,且單模式分析方法僅能識別同一頻率下較強的信號,所以加權(quán)頻譜相干算法無法從模擬數(shù)據(jù)中提取得到偽瑞利波,而前后向最小方差無失真波束形成算法由于其較高的分辨率可以提取偽瑞利波頻散信息。