基于全相位FFT的油套環(huán)空中聲速計算方法

仵 杰，王 軻，劉彥萍，高建申

(1. 西安石油大學(xué)電子工程學(xué)院，陜西西安710065；2. 陜西省油氣井測控技術(shù)重點實驗室，陜西西安710065)

0 引 言

動液面是抽油井在正常生產(chǎn)時，油管和套管環(huán)形空間的一個液面。在油氣井開采過程中，地層的供液能力制約著產(chǎn)量和綜合效率的最大化。動液面不僅是反映地層供液能力的一個重要指標[1]，而且是油田確定合理沉沒度的重要依據(jù)。在動液面測量中，有兩個關(guān)鍵問題：一是聲速的準確確定，二是正確識別液面回波位置(通過到達時刻計算)。關(guān)于聲速的準確確定，目前主要有基于短時平均幅度差函數(shù)(Average Magnitude Difference Function, AMDF)計算抽油井套管環(huán)隙內(nèi)聲速的方法[2]；人工神經(jīng)網(wǎng)絡(luò)方法[3]；基于短時自相關(guān)函數(shù)(Autocorrelation Function, ACF)的套管環(huán)隙內(nèi)聲速計算方法[4-5]。關(guān)于液面回波位置的識別，目前主要有低通橢圓濾波器濾波方法[6]，分型模糊控制濾波方法[7]，基于譜減算法的方法[8]。以上方法雖然取得了一定的成果，能較好地對聲速進行計算以及對動液面信號進行濾波并辨識其位置。但上述方法都存在一些不可避免的缺點，例如基于AMDF算法計算出的聲速并不是很平滑的變化，選取不同的幀對結(jié)果影響很大；譜減算法中的過減系數(shù)和譜減系數(shù)的取值會影響測量結(jié)果，如選取不當(dāng)最終可能導(dǎo)致計算回波位置時產(chǎn)生更大的誤差，準確性反而下降。為了更準確地計算出動液面深度，首先需要準確地估算油套環(huán)空中的聲波速度。針對這一問題，本文在詳細分析油氣井中聲速測量原理的基礎(chǔ)上，提出采用一種性能更優(yōu)的譜分析方法——全相位快速傅里葉變換(all-phase Fast Fourier Transform, apFFT)對接箍信號進行譜分析，然后根據(jù)譜峰位置進一步計算出聲速。全相位 FFT方法的提出主要是為了解決信號直接截斷后其頻譜的幅值、頻率、相位產(chǎn)生的誤差問題[9-11]，目前主要應(yīng)用于巖石聲發(fā)射源定位[12]、合成孔徑水聲通信運動補償[13]、非穩(wěn)態(tài)電力信號[14]、航空設(shè)備的旋轉(zhuǎn)葉片異步振動辨識[15]、感應(yīng)式磁力儀的頻譜校正[16]、時移相位差估計器[17]、數(shù)字濾波器[18]、電力系統(tǒng)相量測量裝置[19]等方面。本文對此方法的實驗驗證中，分別采用快速傅里葉變換(Fast Fourier Transform, FFT)和apFFT對不同信噪比的模擬接箍信號進行頻譜分析與比較，根據(jù)各自的頻譜計算出聲速并比較其精度。然后將兩種方法應(yīng)用于實際接箍信號的處理中，進一步驗證apFFT在聲速計算時的優(yōu)越性。

1 油氣井中聲速測量原理

圖1為油氣井動液面測量的物理模型。動液面測量儀器安裝在井口，儀器產(chǎn)生的脈沖聲源經(jīng)套管和油管之間的環(huán)空傳播到動液面，聲波接收傳感器以一定采樣率采集油套環(huán)空中、動液面和油井周圍的聲波反射信號。油井內(nèi)數(shù)千米長的油管是由許多節(jié)長度基本一致的短節(jié)油管連接而成，在兩節(jié)油管的連接處存在向外凸出的接箍，如圖1所示。采用具有一定壓力的套管氣或壓縮空氣作為發(fā)聲源來產(chǎn)生次聲波，次聲波在套管環(huán)隙內(nèi)向下傳播的過程中，遇到油管接箍和油氣界面后均會產(chǎn)生特定頻率的反射波，如果套管環(huán)隙內(nèi)聲速為一常數(shù)且聲波不衰減以及無其他干擾，則從井口接收的接箍波應(yīng)接近于一個周期信號，但由于套管環(huán)隙內(nèi)各處聲速有可能不同且聲波會衰減并存在其他干擾，因此，在動液面檢測過程中，接箍波一般為準周期信號。深度較淺的油管接箍對聲波具有明顯的反射，油管接箍是周期性存在的，那么反射波形也可視為周期性出現(xiàn)。利用這一特點，可計算波形走過一個油管所需的時間，進而可計算出聲波速度：

圖1 油氣井液位測量示意圖Fig.1 Schematic diagram of liquid level measurement in oil and gas well

2 基于全相位 FFT的動液面聲速計算方法原理

動液面聲速計算方法以譜峰位置的提取作為核心，譜的辨識度及譜峰位置的測量精度至關(guān)重要。全相位FFT法具備相位不變性和克服頻譜泄露的特性[20-22]，從很大程度上解決了傳統(tǒng)FFT法測相位時因數(shù)據(jù)截斷帶來的誤差，同時抑制了諧波的干擾。該算法的數(shù)據(jù)處理精度高、過程簡單、運算量較小，便于實時操作。

下面是對全相位FFT算法原理的簡單介紹。

假設(shè)原始數(shù)據(jù)長度為2N?1的序列x(n),n = ? N + 1 ,? N + 2 ,… , N ?1，對其進行全相位數(shù)據(jù)預(yù)處理的第一步為列出包含序列中心點x(0)的所有長度為 N的序列。 xi(n): i, n = 0,1,… ,N ?1，共形成 N個包含x(0)的N維向量：

當(dāng)f( n)和b( n)均為矩形窗時，數(shù)據(jù)處理成為無窗全相位數(shù)據(jù)預(yù)處理；當(dāng)其中一個為非矩形窗時，數(shù)據(jù)處理稱為單窗全相位預(yù)處理；當(dāng)兩個均為非矩形窗時，數(shù)據(jù)處理稱為雙窗全相位數(shù)據(jù)預(yù)處理。

然后，將全相位預(yù)處理的結(jié)果向量進行傳統(tǒng)的FFT運算，得全相位FFT(apFFT)結(jié)果：

最后，對經(jīng)過全相位處理的N點序列進行峰值搜索即可獲得相對應(yīng)的幅值和頻率。在動液面聲速計算的應(yīng)用中，需對譜峰位置進行提取，將譜峰的橫坐標數(shù)值代入式(1)計算出聲速。

3 正演模型實驗及結(jié)果分析

本文采用COMSOL軟件，在壓力聲學(xué)模塊下，以高斯脈沖作為信號源構(gòu)建理想情況下的動液面聲波反射模型，而在實際生產(chǎn)中，需將聲壓信號轉(zhuǎn)換為電壓信號，所以最后采集到的數(shù)據(jù)其波形幅值單位為V或mV。設(shè)空氣中的聲速為343 m·s-1，采樣頻率為Fs=2 000 Hz。計算接收點處總聲壓場并提取其總聲壓場數(shù)據(jù)，因此所得波形幅值單位為Pa。將其作為實驗所用的模擬信號，波形如圖2所示：

圖2 聲波法測油井模擬信號波形圖及接箍信號波形放大圖Fig.2 The waveforms of analog signal by acoustic wave method for measuring oil wells and enlarged coupling wave respectively

對上述所得接箍數(shù)據(jù)進行傳統(tǒng) FFT與全相位FFT譜分析，求取油套環(huán)隙中的聲速，并對結(jié)果進行分析對比。實驗中，對純凈的模擬接箍信號加入不同頻率、不同強度的噪聲，以測試兩種方法的抗噪性及穩(wěn)定性，并根據(jù)各自所得的頻譜計算聲速，對比其精度。

加入的噪聲信號為

式中：中、低頻正弦干擾的采樣頻率為Fs=2 000 Hz，A為隨機噪聲幅值，2N?1為數(shù)據(jù)采樣點數(shù)，n′ =1 ,2,… 2N?1，r a ndn( 2N?1,1)表示產(chǎn)生2N?1行、1列的正態(tài)分布隨機數(shù)。

對采樣點數(shù)為373的含噪接箍信號分別加入不同強度的噪聲，在同等信噪比條件下分別進行傳統(tǒng)FFT與全相位FFT譜分析，并根據(jù)各自的頻譜求取油套環(huán)隙中的聲速，計算其相對誤差。表1為兩種方法的聲速計算結(jié)果。

表1 FFT和apFFT方法對模擬數(shù)據(jù)聲速估算結(jié)果的對比Table 1 Comparison of the estimated results of acoustic velocity by FFT and apFFT spectrums

由表1可知：在無噪情況下，F(xiàn)FT法所得聲速值與apFFT所得值相同，兩者的相對誤差都很小。但是在加噪情況下，根據(jù)FFT譜 計算出的聲速會出現(xiàn)誤差很大的現(xiàn)象，而根據(jù)apFFT譜計算出的聲速仍保持很小的誤差值。由此可知apFFT法比FFT法的抗干擾能力更強，穩(wěn)定性更好，聲速計算精度更高。

將表 1中信噪比為無噪、2.307 3 dB信號的FFT與apFFT譜圖進行分析對比，結(jié)果如圖3、圖4所示。

圖3 FFT和apFFT得到的無噪接箍信號的頻譜Fig.3 The spectrums of noise-free coupling signal obtained by FFT and apFFT

圖4 含噪接箍信號及其頻譜Fig.4 The noise-containing coupling signal and its spectrums obtained by FFT and apFFT

在無噪情況下，雖然根據(jù)FFT譜和apFFT譜計算出的聲速值相同，但是兩者的頻譜形狀有較大的差異，apFFT譜更為規(guī)則，幾乎沒有譜泄露的現(xiàn)象。對比加噪信號的兩種頻譜，發(fā)現(xiàn)FFT譜已經(jīng)出現(xiàn)譜峰缺少的現(xiàn)象，這會導(dǎo)致聲速計算的誤差增大，精度大幅下降。而apFFT譜仍保持正常數(shù)目的譜峰，為準確計算聲速提供了較好的前提條件。

4 現(xiàn)場接箍反射波頻譜分析與聲速估算

由于實際中無法獲得真正的聲速值，因此估算出的聲速并沒有評價標準。如果實際數(shù)據(jù)中能夠看出接箍波的時域波形特征，則先以此進行時域估算，算出一個聲速值作為評價頻域方法性能的參考。前文講過，聲速在油套環(huán)空中還有變化的可能，那么僅用時域波形估計聲速的方法局限性較大。當(dāng)從深井采集而來的數(shù)據(jù)中接箍波形特征不明顯時，時域估算法的誤差將會很大，且聲速是否在不同深度有所變化也難以獲知。頻域法對波形特征的依賴性較小，能夠在接箍波時域波形特征不明顯時對聲速進行有效估算，而在頻域法中又將尋求精度更高的方法實現(xiàn)對聲速更為準確的估算。

現(xiàn)選取采樣頻率為 470 Hz的實際測井?dāng)?shù)據(jù)中的接箍信號進行FFT和apFFT頻譜分析。淺井信號波形及其分區(qū)如圖5所示。

圖5 淺井信號波形及其分區(qū)Fig.5 The signal waveform of shallow well and its partition

由于淺井信號接箍波段數(shù)據(jù)距發(fā)射源較近的部分時域波形特征較為明顯，隨著波的傳播，其時域波形特征逐漸不明顯，為此截取其中波形明顯的數(shù)據(jù)段進行聲速估算，估算出的聲速值作為后續(xù)研究的參考。圖5的實測信號波形中，第一區(qū)間接箍波波形特征較為明顯，可根據(jù)其時域波形估算聲速。而第二區(qū)間至第五區(qū)間接箍波衰減較嚴重且被噪聲淹沒，因此不易于聲速估算，此時需采用頻域法進行估算。實際中，如果條件惡劣或環(huán)境復(fù)雜將會存在從第一區(qū)間開始接箍波就不易識別的情況，因此根據(jù)時域波形估算聲速的方法并不通用。

檢測出第一區(qū)間接箍波峰所對應(yīng)的采樣序列號為101、128、156、184、213、242，據(jù)此求出聲波信號渡越每個接箍所對應(yīng)采樣序列號的平均值為

則淺井接箍信號第一區(qū)間時域聲速的計算過程為

本文提出采用頻域分析法來估算聲速，分別采用FFT和apFFT算法做譜分析并進一步估算聲速。下面是采用這兩種頻域法對圖5各區(qū)間中接箍波的聲速進行估算的過程，其頻譜如圖6所示。

圖6 FFT和apFFT得到的不同區(qū)間接箍波段數(shù)據(jù)頻譜圖Fig.6 The spectrums of coupling wave data in different regions obtained by FFT and apFFT

可以看出，圖6(a)中第一區(qū)間譜峰規(guī)律性較強，第二區(qū)間到第五區(qū)間譜線雜亂無章，有效峰值少，將會導(dǎo)致聲速計算誤差大。圖6(b)中第一區(qū)間到第三區(qū)間的譜峰規(guī)律性均較強，一致性較好，第四區(qū)間到第五區(qū)間才顯示出峰值很小的譜線。由此說明apFFT算法的性能優(yōu)于FFT算法。

表2為淺井信號不同區(qū)間聲速估算結(jié)果對比。表2中第二列第一個值是根據(jù)接箍信號的時域波形估算的聲速，也以此為參考來評價通過兩種譜分析方法所估算的聲速。由于從第二區(qū)間開始，接箍信號的時域波形已不易辨識，因此無法據(jù)此估算聲速。從表2的第三列和第四列可以看出，根據(jù)接箍波第一區(qū)間的譜分析，F(xiàn)FT法所得聲速值與apFFT所得值相同，第五列和第六列所示兩者的相對誤差都很小；但是在第二區(qū)間和第三區(qū)間，根據(jù)FFT譜計算出的聲速誤差很大，而根據(jù)apFFT譜計算出的聲速仍保持很小的誤差值。此外，根據(jù)apFFT譜分析法還可得知，該井中的聲速基本不變。通過對上述圖表的分析，得出一致結(jié)論：apFFT法比FFT法的抗干擾能力更強，穩(wěn)定性更好，適用于有效信號很微弱的情況，通過該法估算出的聲速精度更高。

表2 淺井信號不同區(qū)間聲速估算結(jié)果對比Table 2 Comparison of the estimated results of acoustic velocity in different regions by FFT and apFFT spectrums

深井信號及接箍波段數(shù)據(jù)的頻譜如圖7所示。

由于深井信號接箍波段的數(shù)據(jù)信號更弱，其時域波形特征更加不明顯，如果根據(jù)其時域波形進行聲速計算，則會存在較大的問題。圖7(a)所示的是來自深井的信號，靠近發(fā)射端的接箍波段信號時域波形仍能呈現(xiàn)出一定的規(guī)律性。為此，截取其中波形特征明顯的數(shù)據(jù)段進行時域的聲速估計，檢測出每個波峰所對應(yīng)的采樣序列號為 77、107、136、166、197。

圖7 深井信號及FFT和apFFT得到的接箍波段數(shù)據(jù)頻譜Fig.7 The deep well signal and the spectrums of coupling wave data obtained by FFT and apFFT

根據(jù)接箍波峰所對應(yīng)的采樣序列號求出聲波信號渡越每個接箍所對應(yīng)的采樣序列號的平均值：

則深井時域聲速計算結(jié)果為

下面采用頻譜分析法計算深井中的聲速。

頻譜分析法中，對接箍波段數(shù)據(jù)的145個采樣點做FFT和apFFT，頻譜圖中展示出其頻譜的一半即可，如圖7(c)和圖7(e)所示。

(1) FFT譜分析法

根據(jù)圖7(c)所示的頻譜圖，檢測各諧波頻率所對應(yīng)的采樣序列號為 7、11、16、20、23，并根據(jù)各諧波頻率所對應(yīng)的采樣序列號求出接箍波基頻所對應(yīng)采樣序列號的平均值為

(2) apFFT譜分析法

根據(jù)圖7(e)所示的頻譜圖，檢測各諧波頻率所對應(yīng)的采樣序列號為 7、11、15、20、25，并根據(jù)各諧波頻率所對應(yīng)的采樣序列號求出接箍波基頻所對應(yīng)采樣序列號的平均值為

由圖7可以看出，通過傳統(tǒng)FFT得到的頻譜中，由于其頻譜泄露較嚴重，導(dǎo)致有效譜峰的位置不易辨識。而通過全相位FFT得到的頻譜，幾乎沒有頻譜泄漏現(xiàn)象，譜峰更為規(guī)則且其位置易于識別，從而有利于接箍反射波頻域特征的辨識與提取。根據(jù)有限的一段接箍波時域波形計算出的聲速為300.8 m·s-1，通過apFFT譜分析法計算出的聲速為 298.73 m·s-1，而根據(jù) FFT譜計算出的聲速為273.83 m·s-1，前兩個值較為接近，再次驗證了由apFFT法計算出的聲速值更為準確。當(dāng)估算更深層井中的聲速時，時域法已不再適用，而基于apFFT的頻域法可發(fā)揮較大的作用。此外，由于聲速的大小與介質(zhì)、溫度、壓力等因素有關(guān)，溫度越低、壓力越小，聲速就越小。隨著井深的增加，壓力減小，因此深井中的聲速較淺井也有所減小。

5 結(jié) 論

本文針對聲波法測油井中接箍信號時域波形特征微弱的情況下不易計算聲速的問題，提出頻域分析法計算聲速。而在接箍信號較弱時，基于FFT的聲速計算方法誤差較大，因此采用基于全相位FFT的聲速計算方法。首先，在設(shè)定的已知聲速下，模擬一組聲波法測油井?dāng)?shù)據(jù)，采用 FFT和 apFFT兩種頻域分析法對接箍波進行頻譜分析，并根據(jù)其頻譜估算聲速，然后以設(shè)定的聲速為參考標準，對比兩種估算方法的精度。通過實驗驗證與結(jié)果分析，可知apFFT能夠有效改善傳統(tǒng)FFT的頻譜泄露問題，在抗噪性和穩(wěn)定性方面更優(yōu)，從而使基于apFFT的聲速計算方法相比基于FFT的聲速計算方法精度更高。最后，選取了兩組分別來自淺井和較深井的現(xiàn)場測井?dāng)?shù)據(jù)，對其中的接箍波段數(shù)據(jù)采用時域法以及基于FFT和apFFT的兩種頻域法進行聲速估算，根據(jù)估算結(jié)果發(fā)現(xiàn)，apFFT法所得結(jié)果與時域法所得結(jié)果更為接近，并在時域法無法估算時仍能進行有效估算，從而再次驗證了apFFT譜分析法進行聲速估算的優(yōu)越性。這也為在更深井中接箍反射信號極微弱的情況下進行聲速估算提供了一種有效的方法。