用于油井动液面检测的频率估算方法
【专利摘要】本发明公开了一种用于油井动液面检测的频率估算方法，包括以下步骤：采集油井动液面的声场信号，得到采样信号x(n)；对采样信号x(n)进行加窗处理，得到加窗函数xw(n)；将加窗函数xw(n)进行离散傅里叶变换，得到频谱Xw(k)；从频谱Xw(k)中寻找幅值最大的频点、最大频谱幅值Xw(l)和第二大频谱幅值Xw(l±1)，其中幅值最大的频点记为第一次迭代时的频率初始值l1；通过插值计算真实频率值λ0；按照公式计算信号的频率fo。有益效果：本算法运用范围广，适用于所有选择的窗函数；计算误差小，受噪声的影响小，具有较好的一致性；计算过程简洁，无需对数据进行预先计算，也不需要存储窗函数的相关参数信息。
【专利说明】
用于油井动液面检测的频率估算方法
技术领域
[0001] 本发明设及油井动液面测量中信号频率的检测技术领域，具体的说是一种用于油 井动液面检测的频率估算方法。
【背景技术】
[0002] 目前已提出的基于管柱声场特性的油井动液面检测方法，其原理是在井口发出声 波激发井下套管内空气柱的共振，利用检测到的空气柱声波共振频率来估算空气柱长度， 进而测量出空气柱下面的动液面深度。其数学模型为
[0003]
[0004] fn为油套环空内空气柱的第n阶共振频率，C为井下声波的传播速度，d为井下管道 的直径。由此数学模型可知，在声速一定的情况下，只要准确测量井下空气柱的共振频率fn 即可得估算出动液面深度。因此，准确高效的估算出井下空气柱的共振频率是测量动液面 深度的关键。
[0005] 但是由于事先并不知道声波信号的频率，很难做到对信号的整周期采样，因此对 采集到信号的傅里叶变换一般都会出现频谱泄漏和栅栏效应，导致对信号频率的估算存在 很大的偏差，人们常采用加窗、插值迭代法来克服缺陷，但是对于不同的窗函数的插值公式 不一样，没有适用于所用窗函数的通用算法，不能满足人们的需求，而且目前运些算法对谱 线错误定位敏感，在大噪声条件下误差较大。
【发明内容】
[0006] 针对上述问题，本发明提供了一种用于油井动液面检测的频率估算方法，来计算 油井动液面的声场信号频率，来解决声场信号频率计算过程复杂、计算误差大的缺陷。
[0007] 为达到上述目的，本发明采用的具体技术方案如下：
[000引一种用于油井动液面检测的频率估算方法，包括W下步骤：
[0009] Sl .要隹細#动漏而的亩扬倍耳.浩幸IlM占长底的采样信号x(n):
[0010]
[00川其中A。表示幅值，f。表示信号的频率，f S表示采样频率，與嗦示相位，n表示采样点 序数，M为对信号的总的采样点数；
[001^ S2:对步骤Sl得到的采样信号x(n)进行jj幡处理；
[0013] 构造 K点长度的窗函数w(n)，对采样信号x(n)进行加窗，得到函数：
[0014] xw(n)=x(n)w(n)；
[0015] S3:将步骤S2得到的加窗函数xw(n)进行离散傅里叶变换，得到频谱XwA):
[0016] S4:从步骤S3所得的频谱Xw化）中寻找幅值最大的频点W及当前最大频谱幅值Xw (1)和第二大频谱幅值Xw(l±l);其中幅值最大的频点记为第一次迭代时的频率初始值h;
[0017] S5:计算真实频率值入0;
[0018] 根据公式lm+l = lm+Sm计算第m+1次迭代时频率估算值，且第m+1次迭代时频率估算 值和第m个频率估算值偏差
其中：
[0019]
，q为矩形窗与选择的窗函数的 扇形损失比，IXw(Im) I为第m次迭代时频率估算值对应的频谱幅值，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次迭代时频率估算值对应 的频谱幅值，Xw( Im-O . 5)为第m次迭代时频率估算值减0.5个步长所在频点的幅值，Xw( Im+ 0.5)为第m次迭代时频率估算值加0.5个步长所在频点的幅值，当I I <T时，T为预设 阔值，迭代结束，并设定真实频率值Ao= lm+1;
[0020] S6:将步骤S5得到的真实频率值值Ao带入公;ES
，计算得到信号的频率f。。
[0021 ]进一步描述，当所加窗函数为矩形窗时，对其进行离散傅里叶变换得到：
[0022]
[0023] 再进一步描述，当所加窗函数为=角窗时，对其进行离散傅里叶变换得到：
[0024]
[0025] 当所加窗函数为汉宁窗时，对其进行离散傅单叶变换得到：
[0026]
。
[0027] 再进一步描述，步骤S5中矩形窗与选择的窗函数的扇形损失比q为
，其中Sk 矩形窗的扇形损失，SLc为选择的窗函数的扇形损失，且窗函数的扇形损失为：
其中W(0.5)为信号最差非整周期采样的情形时候的最大谱线幅值，W(O)为信号整周期采样 时最大谱线幅值。
[0028] 本发明的有益效果:本算法运用范围广，适用于所有的窗函数;计算误差小，可W 克服噪声对谱线造成的影响;计算过程简洁，无需对数据进行预先计算，也不需要存储选择 窗函数的相关参数信息。
【附图说明】
[0029] 图1是本发明的声场信号频率计算流程图；
[0030] 图2是常用窗函数计算信号频率系统误差图；
[0031 ]图3是随机噪声对不同信号频率算法的影响效果图。
【具体实施方式】
[0032] 下面结合附图对本发明的【具体实施方式】W及工作原理作进一步详细说明。
[0033] 从图1可W看出，一种用于油井动液面检测的频率估算方法，包括W下步骤：
[0034] SI:采集油井动液面的声场信号，得到M点长度的采样信号x(n):
[0035]
[0036] 其中A。表示幅值，f。表示信号的频率，fs表示采样频率，卿嗦示相位，n表示采样点 序数，M为对信号的总的采样点数；
[0037] S2:对步骤Sl得到的采样信号x(n)进行加窗处理；
[0038] 构造 K点长度的窗函数w(n)，对采样信号x(n)进行加窗，得到函数:xw(n)=x(n)w (n)；
[0039] 在本实施例中，所加窗函数为矩形窗，则其表达式为:w(n) = 1 (n = 0，1，2…K)，则； xw(n) = x(n)(n=l ,2,…K)
[0040] S3:将步骤S2得到的加窗函数xw(n)进行离散傅里叶变换，得到频谱X^k):
[0041 ]
(I)
[0042] S4:从步骤S3所得的频谱Xw化）中寻找幅值最大的频点W及当前最大频谱幅值Xw (1)和第二大频谱幅值XwQ ± 1);其中幅值最大的频点记为第一次迭代时的频率初始值h;
[0043] S5:计算真实频率值入0:
[0044] 在实际应用中，由于不知道信号的实际频率，因此对信号的截断很难是整周期，那 么进行DFT变换就会出现频谱泄露，即归一化的真实频率值Ao几乎都是出现在DFT变换后的 两个谱线间，因此真实频率值Ao假设表示为：
[0045] A〇=l"+5" (2)
[0046] 其中Im和Sm(-0.5 < 5 < 0.5)分别为整数和小数;Sm为选择的窗函数第m+1次迭代时 频率估算值和第m个频率估算值偏差，Sm具体计算方法为：
[0047] 将公式（2)带入公式（1)可得到：
此处 假设5<A0<N/2-5,则：|W(21m+Sm) I << |W(-Sm(-Sm) I，最大频率幅值的表达式为：
[004引
：3)
[0049] 第二大的频谱幅估表试式为：
[(K)加 ]
(4)
[0051]令最大与第二大频谱幅值之比为：
[0 化 2]
(5)
[0053] 对于任意窗函数的归一化频谱区间[-0.5,0.引内有：
[0054] [Wc化)]q 3 Wr化）
[0055] 其中Wc化）、Wr化)表示任意窗函数和矩形窗的频谱，即[Wc化）]q与Wr化)在归一化频 谱区间[-0.5,0.5]内形状几乎一样，将任意的窗函数逼近矩形窗，则矩形窗对应的最大与
第二大来耐並脑估方hk At .
[0化6]
[0化7] 其中q为矩形窗与选择的窗函数的扇形损失比，且
，其中Sk矩形窗的扇形 损失，SLc为选择的窗函数的扇形损失，且窗函数的扇形损失为：
.其中W(0.5) 为信号最差非整周期采样的情形时候的最大谱线幅值，W(O)为信号整周期采样时最大谱线 幅值；
[005引对于矩形窗，即啊(11) = 1(11 = 0，1，2.'邮，对其进行离散傅里叶变换(0尸1')可得到：
当 |k|/N<<l，则： C6)[0060] 结合公式(3) (4) (5) (6)可得到矩形窗函数存在W下关系：
[0063] 在本实施例中，选择所加的窗函数为矩形窗函数，则矩形窗与选择的窗函数的扇 形损失比
[0061]
[0062]
[0064] 根据步骤S4寻找的第一次迭代时的频率初始值IiW及当前最大频谱幅值Xw(Ii)和 第二大频谱幅值Xw( h ± I)，可计算出第2次迭代时频率估算值和第1次迭代时频率估算值偏 差
> 则第2次迭代时频率估算值为l2 = li+Si， 若满足Sl= I I2-I1 I <T时，T为预设阔值，此时第2次迭代时频率估算值b为真实频率值入0;
[00化]若Si= I h-h I >T，则迭代继续，按照公式
开算第m+1次迭代时频 率估算值和第m次迭代时频率估算值偏差;其中，为了计算第m次迭代时最大与第二大频谱 幅值之比Orm，就要要寻找最大频谱幅值和第二大频谱幅值；
[0066] 首先比较第m-1次迭代时频率估算值对应的频谱幅值Xw(U-I)、第m次迭代时频率 估算值减0.5个步长所在频点的幅值XwQm-0.5) W及第m次迭代时频率估算值加0.5个步长 所在频点的幅值U lm+0.5 )大小，S者中的最大值设为I X ' W( Im-I ) I ;然后比较I X ' W( Im-I ) I与 当第m次迭代时频率估算值对应的频谱幅值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] 即当 Ix'w(lm-I) I < IXwQn {m > 2) ? 当 |X'w(lm-l)| > IXw(Im)I 时： ：引；按照公式lm+l = lm+Sm计算下一 次迭代时频率估算值，当满足|lm+l-l"|<T，此时得到的第m+1个频率估算值为真实频率值
入0;
[0068] S6:将步骤S5得到的真实频率值Ao带入公式 ，计算得到信号的频率f。。
[0069] 计算系统误差：
[0070] 设置采样信号x(n)的幅值Ao设为1，采样频率fs为1024化，数据总长度M为1024点， 故频率分辨率为IHz，对每个窗函数的离散傅里叶变换谱线k，在区间比-0.5，k+0.5)内W步 距0.05进行频率扫描。对于每一个步距频率，相位在区间（-31，JT ]范围内变化，步距为V72。 即对于每个窗函数的离散傅里叶变换谱线产生2880个测试信号，其中与理论信号的最大绝 对误差作为该谱线k的系统误差。
[0071] 从图2可W看出，不同窗函数在该方法下计算出的信号的频率f。的系统误差，其中 窗函数包括矩形窗（Rectangle Window)、海明窗化amming Window)、^角窗（Triangle Window)、汉宁窗（化nning Window)、布莱克曼窗（Blackman Window),从图2可W看出，入0满 足A〇>5和A0<N/2-5系统误差是非常小的，可被忽略，对于所有窗函数，最大频率误差约为 1(T 3,运足W满足大多数工程应用，同时表明该算法可适用于不同的窗函数。
[0072] 检测随机噪声对不同信号频率算法的影响：
[0073] 设置采样信号x(n)的幅值AO设为1，采样频率fs为1024化，数据总长度M为1024点， 故频率分辨率为IHz，为了使得谱线定位错误的几率较高，将信噪比设为-5dB，为了充分验 证谱线定位错误对本算法的影响仿真信号的频率在[255.5,256.5]区间内逐渐增加，步距 为0.025,仿真信号的相位在[-31，31)范围内随机选取。每个频率均独立产生50000个带有加 性高斯白噪声干扰的实例，观察不同算法的方差;从图3可W看出，无论是加矩形窗还是汉 宁窗，本算法(ACIIpDFT)都具有最小的方差，且方差不随频率偏差S改变而改变，具有较好 的一致性。因此，本算法不但具有较好的噪声性且不会受到谱线定位错误的影响。
【主权项】
1. 一种用于油井动液面检测的频率估算方法，其特征在于包括以下步骤： S1:采集油井动液面的声场信号，得到Μ点长度的采样信号x(n):其中A。表示幅值，f。表示信号的频率，:^表示采样频率，_表示相位，η表示采样点序数， Μ为对信号的总的采样点数； S2:对步骤S1得到的采样信号χ(η)进行加窗处理； 构造 Κ点长度的窗函数w(n)，对采样信号χ(η)进行加窗，得到函数:xw(n)=x(n)w(n); S3:将步骤S2得到的加窗函数^(1〇进行离散傅里叶变换，得到频谱Xw(k); S4:从步骤S3所得的频谱Xw(k)中寻找幅值最大的频点以及当前最大频谱幅值Xw( 1)和 第二大频谱幅值Xw(l ± 1);其中幅值最大的频点记为第一次迭代时的频率初始值li; S5:计算真实频率值λ〇; 根据公式lm+1 = lm+Sm计算第m+1次迭代时频率估算值，且第m+1次迭代时频率估算值和第m为矩形窗与选择的窗函数的扇形损失比，I XW(U I为第m次迭代时频率估算值对应的频谱幅 值，|X'w(lm-1) | 为 |Xw(lm-1) |、|Xw(lm-0.5) I 以及 |Xw(lm+0.5) I 之间的最大值，Xw(lm-i)为第m-1次迭代时频率估算值对应的频谱幅值，Xw(lf 〇. 5)为第m次迭代时频率估算值减0.5个步长 所在频点的幅值，Xw( U+0.5)为第m次迭代时频率估算值加0.5个步长所在频点的幅值，当 | <τ时，τ为预设阈值，迭代结束，并设定真实频率值λ〇 = lm+1; S6:将步骤S5得到的真实频率值λ〇带入公式，计算得到信号的频率f。。2. 根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于：当所加窗函 数为矩形窗时，对其进行离散傅里叶变换为得到： 3. 根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于：当所加窗函 数为三角窗时，对其进行离散傅里叶变换得到：4. 根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于：当所加窗函 数为汉宁窗时，对其进行离散傅里叶变换得到：5.根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于:步骤S5中矩 形窗与选择的窗函数的扇形损失比q为：其中SLr矩形窗的扇形损失，SL。为选择的 窗函数的扇形损失，且窗函数的扇形损失为，其中W(0.5)为信号最差非整周期 采样的情形时候的最大谱线幅值，W(0)为信号整周期采样时最大谱线幅值。
【文档编号】E21B47/047GK105822289SQ201610177653
【公开日】2016年8月3日
【申请日】2016年3月25日
【发明人】罗久飞, 周伟, 周盼, 李晓亮, 李太福, 易军, 张元涛, 胡刚
【申请人】重庆科技学院