用于油井動(dòng)液面檢測(cè)的頻率估算方法
【專利摘要】本發(fā)明公開(kāi)了一種用于油井動(dòng)液面檢測(cè)的頻率估算方法，包括以下步驟：采集油井動(dòng)液面的聲場(chǎng)信號(hào)，得到采樣信號(hào)x(n)；對(duì)采樣信號(hào)x(n)進(jìn)行加窗處理，得到加窗函數(shù)xw(n)；將加窗函數(shù)xw(n)進(jìn)行離散傅里葉變換，得到頻譜Xw(k)；從頻譜Xw(k)中尋找幅值最大的頻點(diǎn)、最大頻譜幅值Xw(l)和第二大頻譜幅值Xw(l±1)，其中幅值最大的頻點(diǎn)記為第一次迭代時(shí)的頻率初始值l1；通過(guò)插值計(jì)算真實(shí)頻率值λ0；按照公式計(jì)算信號(hào)的頻率fo。有益效果：本算法運(yùn)用范圍廣，適用于所有選擇的窗函數(shù)；計(jì)算誤差小，受噪聲的影響小，具有較好的一致性；計(jì)算過(guò)程簡(jiǎn)潔，無(wú)需對(duì)數(shù)據(jù)進(jìn)行預(yù)先計(jì)算，也不需要存儲(chǔ)窗函數(shù)的相關(guān)參數(shù)信息。
【專利說(shuō)明】
用于油井動(dòng)液面檢測(cè)的頻率估算方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明設(shè)及油井動(dòng)液面測(cè)量中信號(hào)頻率的檢測(cè)技術(shù)領(lǐng)域，具體的說(shuō)是一種用于油 井動(dòng)液面檢測(cè)的頻率估算方法。
【背景技術(shù)】
[0002] 目前已提出的基于管柱聲場(chǎng)特性的油井動(dòng)液面檢測(cè)方法，其原理是在井口發(fā)出聲 波激發(fā)井下套管內(nèi)空氣柱的共振，利用檢測(cè)到的空氣柱聲波共振頻率來(lái)估算空氣柱長(zhǎng)度， 進(jìn)而測(cè)量出空氣柱下面的動(dòng)液面深度。其數(shù)學(xué)模型為
[0003]
[0004] fn為油套環(huán)空內(nèi)空氣柱的第n階共振頻率，C為井下聲波的傳播速度，d為井下管道 的直徑。由此數(shù)學(xué)模型可知，在聲速一定的情況下，只要準(zhǔn)確測(cè)量井下空氣柱的共振頻率fn 即可得估算出動(dòng)液面深度。因此，準(zhǔn)確高效的估算出井下空氣柱的共振頻率是測(cè)量動(dòng)液面 深度的關(guān)鍵。
[0005] 但是由于事先并不知道聲波信號(hào)的頻率，很難做到對(duì)信號(hào)的整周期采樣，因此對(duì) 采集到信號(hào)的傅里葉變換一般都會(huì)出現(xiàn)頻譜泄漏和柵欄效應(yīng)，導(dǎo)致對(duì)信號(hào)頻率的估算存在 很大的偏差，人們常采用加窗、插值迭代法來(lái)克服缺陷，但是對(duì)于不同的窗函數(shù)的插值公式 不一樣，沒(méi)有適用于所用窗函數(shù)的通用算法，不能滿足人們的需求，而且目前運(yùn)些算法對(duì)譜 線錯(cuò)誤定位敏感，在大噪聲條件下誤差較大。

【發(fā)明內(nèi)容】

[0006] 針對(duì)上述問(wèn)題，本發(fā)明提供了一種用于油井動(dòng)液面檢測(cè)的頻率估算方法，來(lái)計(jì)算 油井動(dòng)液面的聲場(chǎng)信號(hào)頻率，來(lái)解決聲場(chǎng)信號(hào)頻率計(jì)算過(guò)程復(fù)雜、計(jì)算誤差大的缺陷。
[0007] 為達(dá)到上述目的，本發(fā)明采用的具體技術(shù)方案如下：
[000引一種用于油井動(dòng)液面檢測(cè)的頻率估算方法，包括W下步驟：
[0009] Sl .要隹細(xì)#動(dòng)漏而的畝揚(yáng)倍耳.浩幸IlM占長(zhǎng)底的采樣信號(hào)x(n):
[0010]
[00川其中A。表示幅值，f。表示信號(hào)的頻率，f S表示采樣頻率，與嗦示相位，n表示采樣點(diǎn) 序數(shù)，M為對(duì)信號(hào)的總的采樣點(diǎn)數(shù)；
[001^ S2:對(duì)步驟Sl得到的采樣信號(hào)x(n)進(jìn)行jj幡處理；
[0013] 構(gòu)造 K點(diǎn)長(zhǎng)度的窗函數(shù)w(n)，對(duì)采樣信號(hào)x(n)進(jìn)行加窗，得到函數(shù)：
[0014] xw(n)=x(n)w(n)；
[0015] S3:將步驟S2得到的加窗函數(shù)xw(n)進(jìn)行離散傅里葉變換，得到頻譜XwA):
[0016] S4:從步驟S3所得的頻譜Xw化）中尋找幅值最大的頻點(diǎn)W及當(dāng)前最大頻譜幅值Xw (1)和第二大頻譜幅值Xw(l±l);其中幅值最大的頻點(diǎn)記為第一次迭代時(shí)的頻率初始值h;
[0017] S5:計(jì)算真實(shí)頻率值入0;
[0018] 根據(jù)公式lm+l = lm+Sm計(jì)算第m+1次迭代時(shí)頻率估算值，且第m+1次迭代時(shí)頻率估算 值和第m個(gè)頻率估算值偏差
其中：
[0019]
，q為矩形窗與選擇的窗函數(shù)的 扇形損失比，IXw(Im) I為第m次迭代時(shí)頻率估算值對(duì)應(yīng)的頻譜幅值，I X'w(lm-l) I為IXw(Im-I) 、I XwQm-O . 5) I W及I Xw(lm+0.5) I之間的最大值，Xw(Im-I)為第m-1次迭代時(shí)頻率估算值對(duì)應(yīng) 的頻譜幅值，Xw( Im-O . 5)為第m次迭代時(shí)頻率估算值減0.5個(gè)步長(zhǎng)所在頻點(diǎn)的幅值，Xw( Im+ 0.5)為第m次迭代時(shí)頻率估算值加0.5個(gè)步長(zhǎng)所在頻點(diǎn)的幅值，當(dāng)I I <T時(shí)，T為預(yù)設(shè) 闊值，迭代結(jié)束，并設(shè)定真實(shí)頻率值A(chǔ)o= lm+1;
[0020] S6:將步驟S5得到的真實(shí)頻率值值A(chǔ)o帶入公;ES
，計(jì)算得到信號(hào)的頻率f。。
[0021 ]進(jìn)一步描述，當(dāng)所加窗函數(shù)為矩形窗時(shí)，對(duì)其進(jìn)行離散傅里葉變換得到：
[0022]
[0023] 再進(jìn)一步描述，當(dāng)所加窗函數(shù)為=角窗時(shí)，對(duì)其進(jìn)行離散傅里葉變換得到：
[0024]
[0025] 當(dāng)所加窗函數(shù)為漢寧窗時(shí)，對(duì)其進(jìn)行離散傅單葉變換得到：
[0026]
。
[0027] 再進(jìn)一步描述，步驟S5中矩形窗與選擇的窗函數(shù)的扇形損失比q為
，其中Sk 矩形窗的扇形損失，SLc為選擇的窗函數(shù)的扇形損失，且窗函數(shù)的扇形損失為：
其中W(0.5)為信號(hào)最差非整周期采樣的情形時(shí)候的最大譜線幅值，W(O)為信號(hào)整周期采樣 時(shí)最大譜線幅值。
[0028] 本發(fā)明的有益效果:本算法運(yùn)用范圍廣，適用于所有的窗函數(shù);計(jì)算誤差小，可W 克服噪聲對(duì)譜線造成的影響;計(jì)算過(guò)程簡(jiǎn)潔，無(wú)需對(duì)數(shù)據(jù)進(jìn)行預(yù)先計(jì)算，也不需要存儲(chǔ)選擇 窗函數(shù)的相關(guān)參數(shù)信息。
【附圖說(shuō)明】
[0029] 圖1是本發(fā)明的聲場(chǎng)信號(hào)頻率計(jì)算流程圖；
[0030] 圖2是常用窗函數(shù)計(jì)算信號(hào)頻率系統(tǒng)誤差圖；
[0031 ]圖3是隨機(jī)噪聲對(duì)不同信號(hào)頻率算法的影響效果圖。
【具體實(shí)施方式】
[0032] 下面結(jié)合附圖對(duì)本發(fā)明的【具體實(shí)施方式】W及工作原理作進(jìn)一步詳細(xì)說(shuō)明。
[0033] 從圖1可W看出，一種用于油井動(dòng)液面檢測(cè)的頻率估算方法，包括W下步驟：
[0034] SI:采集油井動(dòng)液面的聲場(chǎng)信號(hào)，得到M點(diǎn)長(zhǎng)度的采樣信號(hào)x(n):
[0035]
[0036] 其中A。表示幅值，f。表示信號(hào)的頻率，fs表示采樣頻率，卿嗦示相位，n表示采樣點(diǎn) 序數(shù)，M為對(duì)信號(hào)的總的采樣點(diǎn)數(shù)；
[0037] S2:對(duì)步驟Sl得到的采樣信號(hào)x(n)進(jìn)行加窗處理；
[0038] 構(gòu)造 K點(diǎn)長(zhǎng)度的窗函數(shù)w(n)，對(duì)采樣信號(hào)x(n)進(jìn)行加窗，得到函數(shù):xw(n)=x(n)w (n)；
[0039] 在本實(shí)施例中，所加窗函數(shù)為矩形窗，則其表達(dá)式為:w(n) = 1 (n = 0，1，2…K)，則； xw(n) = x(n)(n=l ,2,…K)
[0040] S3:將步驟S2得到的加窗函數(shù)xw(n)進(jìn)行離散傅里葉變換，得到頻譜X^k):
[0041 ]
(I)
[0042] S4:從步驟S3所得的頻譜Xw化）中尋找幅值最大的頻點(diǎn)W及當(dāng)前最大頻譜幅值Xw (1)和第二大頻譜幅值XwQ ± 1);其中幅值最大的頻點(diǎn)記為第一次迭代時(shí)的頻率初始值h;
[0043] S5:計(jì)算真實(shí)頻率值入0:
[0044] 在實(shí)際應(yīng)用中，由于不知道信號(hào)的實(shí)際頻率，因此對(duì)信號(hào)的截?cái)嗪茈y是整周期，那 么進(jìn)行DFT變換就會(huì)出現(xiàn)頻譜泄露，即歸一化的真實(shí)頻率值A(chǔ)o幾乎都是出現(xiàn)在DFT變換后的 兩個(gè)譜線間，因此真實(shí)頻率值A(chǔ)o假設(shè)表示為：
[0045] A〇=l"+5" (2)
[0046] 其中Im和Sm(-0.5 < 5 < 0.5)分別為整數(shù)和小數(shù);Sm為選擇的窗函數(shù)第m+1次迭代時(shí) 頻率估算值和第m個(gè)頻率估算值偏差，Sm具體計(jì)算方法為：
[0047] 將公式（2)帶入公式（1)可得到：
此處 假設(shè)5<A0<N/2-5,則：|W(21m+Sm) I << |W(-Sm(-Sm) I，最大頻率幅值的表達(dá)式為：
[004引
：3)
[0049] 第二大的頻譜幅估表試式為：
[(K)加 ]
(4)
[0051]令最大與第二大頻譜幅值之比為：
[0 化 2]
(5)
[0053] 對(duì)于任意窗函數(shù)的歸一化頻譜區(qū)間[-0.5,0.引內(nèi)有：
[0054] [Wc化)]q 3 Wr化）
[0055] 其中Wc化）、Wr化)表示任意窗函數(shù)和矩形窗的頻譜，即[Wc化）]q與Wr化)在歸一化頻 譜區(qū)間[-0.5,0.5]內(nèi)形狀幾乎一樣，將任意的窗函數(shù)逼近矩形窗，則矩形窗對(duì)應(yīng)的最大與
第二大來(lái)耐並腦估方hk At .
[0化6]
[0化7] 其中q為矩形窗與選擇的窗函數(shù)的扇形損失比，且
，其中Sk矩形窗的扇形 損失，SLc為選擇的窗函數(shù)的扇形損失，且窗函數(shù)的扇形損失為：
.其中W(0.5) 為信號(hào)最差非整周期采樣的情形時(shí)候的最大譜線幅值，W(O)為信號(hào)整周期采樣時(shí)最大譜線 幅值；
[005引對(duì)于矩形窗，即啊(11) = 1(11 = 0，1，2.'郵，對(duì)其進(jìn)行離散傅里葉變換(0尸1')可得到：
當(dāng) |k|/N<<l，則： C6)[0060] 結(jié)合公式(3) (4) (5) (6)可得到矩形窗函數(shù)存在W下關(guān)系：
[0063] 在本實(shí)施例中，選擇所加的窗函數(shù)為矩形窗函數(shù)，則矩形窗與選擇的窗函數(shù)的扇 形損失比
[0061]
[0062]
[0064] 根據(jù)步驟S4尋找的第一次迭代時(shí)的頻率初始值IiW及當(dāng)前最大頻譜幅值Xw(Ii)和 第二大頻譜幅值Xw( h ± I)，可計(jì)算出第2次迭代時(shí)頻率估算值和第1次迭代時(shí)頻率估算值偏 差
> 則第2次迭代時(shí)頻率估算值為l2 = li+Si， 若滿足Sl= I I2-I1 I <T時(shí)，T為預(yù)設(shè)闊值，此時(shí)第2次迭代時(shí)頻率估算值b為真實(shí)頻率值入0;
[00化]若Si= I h-h I >T，則迭代繼續(xù)，按照公式
開(kāi)算第m+1次迭代時(shí)頻 率估算值和第m次迭代時(shí)頻率估算值偏差;其中，為了計(jì)算第m次迭代時(shí)最大與第二大頻譜 幅值之比Orm，就要要尋找最大頻譜幅值和第二大頻譜幅值；
[0066] 首先比較第m-1次迭代時(shí)頻率估算值對(duì)應(yīng)的頻譜幅值Xw(U-I)、第m次迭代時(shí)頻率 估算值減0.5個(gè)步長(zhǎng)所在頻點(diǎn)的幅值XwQm-0.5) W及第m次迭代時(shí)頻率估算值加0.5個(gè)步長(zhǎng) 所在頻點(diǎn)的幅值U lm+0.5 )大小，S者中的最大值設(shè)為I X ' W( Im-I ) I ;然后比較I X ' W( Im-I ) I與 當(dāng)?shù)趍次迭代時(shí)頻率估算值對(duì)應(yīng)的頻譜幅值IXw(Im) I的大小，若Ix'w(lm-l) I < IXw(Im) I，則為 乂^1。）|為最大頻譜幅值，片^1。-1)|為第二大頻譜幅值;若|乂\(1。-1)|>|乂^1。）|，|乂\ (Im-I) I最大頻譜幅值，I X'w(lm-l) I為第二大頻譜幅值；
[0067] 即當(dāng) Ix'w(lm-I) I < IXwQn {m > 2) ? 當(dāng) |X'w(lm-l)| > IXw(Im)I 時(shí)： ：引；按照公式lm+l = lm+Sm計(jì)算下一 次迭代時(shí)頻率估算值，當(dāng)滿足|lm+l-l"|<T，此時(shí)得到的第m+1個(gè)頻率估算值為真實(shí)頻率值
入0;
[0068] S6:將步驟S5得到的真實(shí)頻率值A(chǔ)o帶入公式 ，計(jì)算得到信號(hào)的頻率f。。
[0069] 計(jì)算系統(tǒng)誤差：
[0070] 設(shè)置采樣信號(hào)x(n)的幅值A(chǔ)o設(shè)為1，采樣頻率fs為1024化，數(shù)據(jù)總長(zhǎng)度M為1024點(diǎn)， 故頻率分辨率為IHz，對(duì)每個(gè)窗函數(shù)的離散傅里葉變換譜線k，在區(qū)間比-0.5，k+0.5)內(nèi)W步 距0.05進(jìn)行頻率掃描。對(duì)于每一個(gè)步距頻率，相位在區(qū)間（-31，JT ]范圍內(nèi)變化，步距為V72。 即對(duì)于每個(gè)窗函數(shù)的離散傅里葉變換譜線產(chǎn)生2880個(gè)測(cè)試信號(hào)，其中與理論信號(hào)的最大絕 對(duì)誤差作為該譜線k的系統(tǒng)誤差。
[0071] 從圖2可W看出，不同窗函數(shù)在該方法下計(jì)算出的信號(hào)的頻率f。的系統(tǒng)誤差，其中 窗函數(shù)包括矩形窗（Rectangle Window)、海明窗化amming Window)、^角窗（Triangle Window)、漢寧窗（化nning Window)、布萊克曼窗（Blackman Window),從圖2可W看出，入0滿 足A〇>5和A0<N/2-5系統(tǒng)誤差是非常小的，可被忽略，對(duì)于所有窗函數(shù)，最大頻率誤差約為 1(T 3,運(yùn)足W滿足大多數(shù)工程應(yīng)用，同時(shí)表明該算法可適用于不同的窗函數(shù)。
[0072] 檢測(cè)隨機(jī)噪聲對(duì)不同信號(hào)頻率算法的影響：
[0073] 設(shè)置采樣信號(hào)x(n)的幅值A(chǔ)O設(shè)為1，采樣頻率fs為1024化，數(shù)據(jù)總長(zhǎng)度M為1024點(diǎn)， 故頻率分辨率為IHz，為了使得譜線定位錯(cuò)誤的幾率較高，將信噪比設(shè)為-5dB，為了充分驗(yàn) 證譜線定位錯(cuò)誤對(duì)本算法的影響仿真信號(hào)的頻率在[255.5,256.5]區(qū)間內(nèi)逐漸增加，步距 為0.025,仿真信號(hào)的相位在[-31，31)范圍內(nèi)隨機(jī)選取。每個(gè)頻率均獨(dú)立產(chǎn)生50000個(gè)帶有加 性高斯白噪聲干擾的實(shí)例，觀察不同算法的方差;從圖3可W看出，無(wú)論是加矩形窗還是漢 寧窗，本算法(ACIIpDFT)都具有最小的方差，且方差不隨頻率偏差S改變而改變，具有較好 的一致性。因此，本算法不但具有較好的噪聲性且不會(huì)受到譜線定位錯(cuò)誤的影響。
【主權(quán)項(xiàng)】
1. 一種用于油井動(dòng)液面檢測(cè)的頻率估算方法，其特征在于包括以下步驟： S1:采集油井動(dòng)液面的聲場(chǎng)信號(hào)，得到Μ點(diǎn)長(zhǎng)度的采樣信號(hào)x(n):其中A。表示幅值，f。表示信號(hào)的頻率，:^表示采樣頻率，_表示相位，η表示采樣點(diǎn)序數(shù)， Μ為對(duì)信號(hào)的總的采樣點(diǎn)數(shù)； S2:對(duì)步驟S1得到的采樣信號(hào)χ(η)進(jìn)行加窗處理； 構(gòu)造 Κ點(diǎn)長(zhǎng)度的窗函數(shù)w(n)，對(duì)采樣信號(hào)χ(η)進(jìn)行加窗，得到函數(shù):xw(n)=x(n)w(n); S3:將步驟S2得到的加窗函數(shù)^(1〇進(jìn)行離散傅里葉變換，得到頻譜Xw(k); S4:從步驟S3所得的頻譜Xw(k)中尋找幅值最大的頻點(diǎn)以及當(dāng)前最大頻譜幅值Xw( 1)和 第二大頻譜幅值Xw(l ± 1);其中幅值最大的頻點(diǎn)記為第一次迭代時(shí)的頻率初始值li; S5:計(jì)算真實(shí)頻率值λ〇; 根據(jù)公式lm+1 = lm+Sm計(jì)算第m+1次迭代時(shí)頻率估算值，且第m+1次迭代時(shí)頻率估算值和第m為矩形窗與選擇的窗函數(shù)的扇形損失比，I XW(U I為第m次迭代時(shí)頻率估算值對(duì)應(yīng)的頻譜幅 值，|X'w(lm-1) | 為 |Xw(lm-1) |、|Xw(lm-0.5) I 以及 |Xw(lm+0.5) I 之間的最大值，Xw(lm-i)為第m-1次迭代時(shí)頻率估算值對(duì)應(yīng)的頻譜幅值，Xw(lf 〇. 5)為第m次迭代時(shí)頻率估算值減0.5個(gè)步長(zhǎng) 所在頻點(diǎn)的幅值，Xw( U+0.5)為第m次迭代時(shí)頻率估算值加0.5個(gè)步長(zhǎng)所在頻點(diǎn)的幅值，當(dāng) | <τ時(shí)，τ為預(yù)設(shè)閾值，迭代結(jié)束，并設(shè)定真實(shí)頻率值λ〇 = lm+1; S6:將步驟S5得到的真實(shí)頻率值λ〇帶入公式，計(jì)算得到信號(hào)的頻率f。。2. 根據(jù)權(quán)利要求1所述用于油井動(dòng)液面檢測(cè)的頻率估算方法，其特征在于：當(dāng)所加窗函 數(shù)為矩形窗時(shí)，對(duì)其進(jìn)行離散傅里葉變換為得到： 3. 根據(jù)權(quán)利要求1所述用于油井動(dòng)液面檢測(cè)的頻率估算方法，其特征在于：當(dāng)所加窗函 數(shù)為三角窗時(shí)，對(duì)其進(jìn)行離散傅里葉變換得到：4. 根據(jù)權(quán)利要求1所述用于油井動(dòng)液面檢測(cè)的頻率估算方法，其特征在于：當(dāng)所加窗函 數(shù)為漢寧窗時(shí)，對(duì)其進(jìn)行離散傅里葉變換得到：5.根據(jù)權(quán)利要求1所述用于油井動(dòng)液面檢測(cè)的頻率估算方法，其特征在于:步驟S5中矩 形窗與選擇的窗函數(shù)的扇形損失比q為：其中SLr矩形窗的扇形損失，SL。為選擇的 窗函數(shù)的扇形損失，且窗函數(shù)的扇形損失為，其中W(0.5)為信號(hào)最差非整周期 采樣的情形時(shí)候的最大譜線幅值，W(0)為信號(hào)整周期采樣時(shí)最大譜線幅值。
【文檔編號(hào)】E21B47/047GK105822289SQ201610177653
【公開(kāi)日】2016年8月3日
【申請(qǐng)日】2016年3月25日
【發(fā)明人】羅久飛, 周偉, 周盼, 李曉亮, 李太福, 易軍, 張?jiān)獫? 胡剛
【申請(qǐng)人】重慶科技學(xué)院