CN105822289A - 用于油井动液面检测的频率估算方法 - Google Patents

Abstract

本发明公开了一种用于油井动液面检测的频率估算方法，包括以下步骤：采集油井动液面的声场信号，得到采样信号x(n)；对采样信号x(n)进行加窗处理，得到加窗函数x_w(n)；将加窗函数x_w(n)进行离散傅里叶变换，得到频谱X_w(k)；从频谱X_w(k)中寻找幅值最大的频点、最大频谱幅值X_w(l)和第二大频谱幅值X_w(l±1)，其中幅值最大的频点记为第一次迭代时的频率初始值l₁；通过插值计算真实频率值λ₀；按照公式计算信号的频率f_o。有益效果：本算法运用范围广，适用于所有选择的窗函数；计算误差小，受噪声的影响小，具有较好的一致性；计算过程简洁，无需对数据进行预先计算，也不需要存储窗函数的相关参数信息。

Description

用于油井动液面检测的频率估算方法

技术领域

本发明涉及油井动液面测量中信号频率的检测技术领域，具体的说是一种用于油井动液面检测的频率估算方法。

背景技术

目前已提出的基于管柱声场特性的油井动液面检测方法，其原理是在井口发出声波激发井下套管内空气柱的共振，利用检测到的空气柱声波共振频率来估算空气柱长度，进而测量出空气柱下面的动液面深度。其数学模型为

$f_n=\frac{(2n-1)c}{4l}-0.3d,(n=1,2,3,...)$

f_n为油套环空内空气柱的第n阶共振频率，c为井下声波的传播速度，d为井下管道的直径。由此数学模型可知，在声速一定的情况下，只要准确测量井下空气柱的共振频率f_n即可得估算出动液面深度。因此，准确高效的估算出井下空气柱的共振频率是测量动液面深度的关键。

但是由于事先并不知道声波信号的频率，很难做到对信号的整周期采样，因此对采集到信号的傅里叶变换一般都会出现频谱泄漏和栅栏效应，导致对信号频率的估算存在很大的偏差，人们常采用加窗、插值迭代法来克服缺陷，但是对于不同的窗函数的插值公式不一样，没有适用于所用窗函数的通用算法，不能满足人们的需求，而且目前这些算法对谱线错误定位敏感，在大噪声条件下误差较大。

发明内容

针对上述问题，本发明提供了一种用于油井动液面检测的频率估算方法，来计算油井动液面的声场信号频率，来解决声场信号频率计算过程复杂、计算误差大的缺陷。

为达到上述目的，本发明采用的具体技术方案如下：

一种用于油井动液面检测的频率估算方法，包括以下步骤：

S1：采集油井动液面的声场信号，得到M点长度的采样信号x(n)：

其中A_o表示幅值，f_o表示信号的频率，f_s表示采样频率，表示相位，n表示采样点序数，M为对信号的总的采样点数；

S2：对步骤S1得到的采样信号x(n)进行加窗处理；

构造K点长度的窗函数w(n)，对采样信号x(n)进行加窗，得到函数：

x_w(n)＝x(n)w(n)；

S3：将步骤S2得到的加窗函数x_w(n)进行离散傅里叶变换，得到频谱X_w(k)：

S4：从步骤S3所得的频谱X_w(k)中寻找幅值最大的频点以及当前最大频谱幅值X_w(l)和第二大频谱幅值X_w(l±1)；其中幅值最大的频点记为第一次迭代时的频率初始值l₁；

S5：计算真实频率值λ₀；

根据公式l_m+1＝l_m+δ_m计算第m+1次迭代时频率估算值，且第m+1次迭代时频率估算值和第m个频率估算值偏差其中：

q为矩形窗与选择的窗函数的扇形损失比，|X_w(l_m)|为第m次迭代时频率估算值对应的频谱幅值，|X'_w(l_m-1)|为|X_w(l_m-1)|、|X_w(l_m-0.5)|以及|X_w(l_m+0.5)|之间的最大值，X_w(l_m-1)为第m-1次迭代时频率估算值对应的频谱幅值，X_w(l_m-0.5)为第m次迭代时频率估算值减0.5个步长所在频点的幅值，X_w(l_m+0.5)为第m次迭代时频率估算值加0.5个步长所在频点的幅值，当|l_m+1-l_m|＜τ时，τ为预设阈值，迭代结束，并设定真实频率值λ₀＝l_m+1；

S6：将步骤S5得到的真实频率值值λ₀带入公式计算得到信号的频率f_o。

进一步描述，当所加窗函数为矩形窗时，对其进行离散傅里叶变换得到：

$W\left(k\right)=\frac{sin\mathrm{\π}k}{\mathrm{\π}k}\·e^{-j\mathrm{\π}k}.$

再进一步描述，当所加窗函数为三角窗时，对其进行离散傅里叶变换得到：

$W\left(k\right)=\frac{2}{N}{e}^{-j\frac{N-2}{N}\·k\mathrm{\π}}{\[\frac{sin\frac{k\mathrm{\π}}{2}}{\sin \frac{k\mathrm{\π}}{N}}\]}^2.$

当所加窗函数为汉宁窗时，对其进行离散傅里叶变换得到：

$W\left(k\right)=\sin \mathrm{\π}k\·e^{-jk\mathrm{\π}}\[\frac{1}{4N}\frac{\sin \frac{2\mathrm{\π}k}{N}}{{\left(\sin \frac{\mathrm{\π}k}{N}\right)}^2}-\frac{1}{4N}\frac{\sin \frac{2\mathrm{\π}k}{N}}{\sin \frac{\mathrm{\π}(k-1)}{N}\·\sin \frac{\mathrm{\π}(k+1)}{N}}\].$

再进一步描述，步骤S5中矩形窗与选择的窗函数的扇形损失比q为：其中SL_r矩形窗的扇形损失，SL_c为选择的窗函数的扇形损失，且窗函数的扇形损失为：其中W(0.5)为信号最差非整周期采样的情形时候的最大谱线幅值，W(0)为信号整周期采样时最大谱线幅值。

本发明的有益效果：本算法运用范围广，适用于所有的窗函数；计算误差小，可以克服噪声对谱线造成的影响；计算过程简洁，无需对数据进行预先计算，也不需要存储选择窗函数的相关参数信息。

附图说明

图1是本发明的声场信号频率计算流程图；

图2是常用窗函数计算信号频率系统误差图；

图3是随机噪声对不同信号频率算法的影响效果图。

具体实施方式

下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。

从图1可以看出，一种用于油井动液面检测的频率估算方法，包括以下步骤：

S1：采集油井动液面的声场信号，得到M点长度的采样信号x(n)：

其中A_o表示幅值，f_o表示信号的频率，f_s表示采样频率，表示相位，n表示采样点序数，M为对信号的总的采样点数；

S2：对步骤S1得到的采样信号x(n)进行加窗处理；

构造K点长度的窗函数w(n)，对采样信号x(n)进行加窗，得到函数：x_w(n)＝x(n)w(n)；

在本实施例中，所加窗函数为矩形窗，则其表达式为：w(n)＝1(n＝0,1,2…K)，则；x_w(n)＝x(n)(n＝1,2,…K)

S3：将步骤S2得到的加窗函数x_w(n)进行离散傅里叶变换，得到频谱X_w(k)：

$\begin{array}{c}{X}_w\left(k\right)=\underset{n=0}{\overset{N-1}{\Σ}}{x}_w\left(n\right)e^{-j\frac{2\mathrm{\π}}{N}nk}\\ \begin{array}{cc}=\frac{A_0}{2}{e}^{{\mathrm{j\φ}}_0}W(k-{\mathrm{\λ}}_0)+\frac{A_0}{2}{e}^{-{\mathrm{j\φ}}_0}W(k+{\mathrm{\λ}}_0)& 0\≤k\≤N-1\end{array}\end{array}---\left(1\right)$

S4：从步骤S3所得的频谱X_w(k)中寻找幅值最大的频点以及当前最大频谱幅值X_w(l)和第二大频谱幅值X_w(l±1)；其中幅值最大的频点记为第一次迭代时的频率初始值l₁；

S5：计算真实频率值λ₀：

在实际应用中，由于不知道信号的实际频率，因此对信号的截断很难是整周期，那么进行DFT变换就会出现频谱泄露，即归一化的真实频率值λ₀几乎都是出现在DFT变换后的两个谱线间，因此真实频率值λ₀假设表示为：

λ₀＝l_m+δ_m(2)

其中l_m和δ_m(-0.5≤δ≤0.5)分别为整数和小数；δ_m为选择的窗函数第m+1次迭代时频率估算值和第m个频率估算值偏差，δ_m具体计算方法为：

将公式(2)带入公式(1)可得到：此处假设5＜λ₀＜N/2-5，则：|W(2l_m+δ_m)|＜＜|W(-δ_m(-δ_m)|，最大频率幅值的表达式为：

$\left|X_w\left(l_m\right)\right|\≈\frac{A_0}{2}\left|W(-{\mathrm{\δ}}_m)\right|---\left(3\right)$

第二大的频谱幅值表达式为：

$\left|X_w(l_m\±1)\right|\≈\frac{A_0}{2}\left|W(-{\mathrm{\δ}}_m\±1)\right|---\left(4\right)$

令最大与第二大频谱幅值之比为：

${\mathrm{\α}}_m=\frac{\left|X_w(l_m\±1)\right|}{\left|X_w\left(l_m\right)\right|}---\left(5\right)$

对于任意窗函数的归一化频谱区间[-0.5,0.5]内有：

[W_c(k)]^q≈W_r(k)

其中W_c(k)、W_r(k)表示任意窗函数和矩形窗的频谱，即[W_c(k)]^q与W_r(k)在归一化频谱区间[-0.5,0.5]内形状几乎一样,将任意的窗函数逼近矩形窗，则矩形窗对应的最大与第二大频谱幅值之比为：

${\mathrm{\α}}_{rm}={{\mathrm{\α}}_{cm}}^q=\frac{|X_w(l_m\±1){|}^q}{|X_w\left(l_m\right){|}^q}$

其中q为矩形窗与选择的窗函数的扇形损失比，且：其中SL_r矩形窗的扇形损失，SL_c为选择的窗函数的扇形损失，且窗函数的扇形损失为：其中W(0.5)为信号最差非整周期采样的情形时候的最大谱线幅值，W(0)为信号整周期采样时最大谱线幅值；

对于矩形窗，即w_r(n)＝1(n＝0,1,2…N)，对其进行离散傅里叶变换(DFT)可得到：当|k|/N＜＜1，则：

$W_r\left(k\right)=\frac{sin\mathrm{\π}k}{\mathrm{\π}k}\·e^{-j\mathrm{\π}k}---\left(6\right)$

结合公式(3)(4)(5)(6)可得到矩形窗函数存在以下关系：

${\mathrm{\α}}_{rm}\≈\frac{\left|X_{wr}(l_m\±1)\right|}{\left|X_{wr}\left(l_m\right)\right|}=\frac{\left|W_r(-{\mathrm{\δ}}_m\±1)\right|}{\left|W_r(-{\mathrm{\δ}}_m)\right|}=\frac{\left|\frac{sin\mathrm{\π}(-{\mathrm{\δ}}_m\±1)}{\mathrm{\π}(-{\mathrm{\δ}}_m\±1)}\right|}{\left|\frac{sin\mathrm{\π}(-{\mathrm{\δ}}_m)}{\mathrm{\π}(-{\mathrm{\δ}}_m)}\right|}\frac{\left|{\mathrm{\δ}}_m\right|}{|-{\mathrm{\δ}}_m\±1|}\≈\frac{\left|{\mathrm{\δ}}_m\right|}{|-{\mathrm{\δ}}_m\±1|}$

则频率偏差

在本实施例中，选择所加的窗函数为矩形窗函数，则矩形窗与选择的窗函数的扇形损失比则频率偏差

根据步骤S4寻找的第一次迭代时的频率初始值l₁以及当前最大频谱幅值X_w(l₁)和第二大频谱幅值X_w(l₁±1)，可计算出第2次迭代时频率估算值和第1次迭代时频率估算值偏差其中则第2次迭代时频率估算值为l₂＝l₁+δ₁，若满足δ₁＝|l₂-l₁|＜τ时，τ为预设阈值，此时第2次迭代时频率估算值l₂为真实频率值λ₀；

若δ₁＝|l₂-l₁|＞τ，则迭代继续，按照公式计算第m+1次迭代时频率估算值和第m次迭代时频率估算值偏差；其中，为了计算第m次迭代时最大与第二大频谱幅值之比α_rm，就要要寻找最大频谱幅值和第二大频谱幅值；

首先比较第m-1次迭代时频率估算值对应的频谱幅值X_w(l_m-1)、第m次迭代时频率估算值减0.5个步长所在频点的幅值X_w(l_m-0.5)以及第m次迭代时频率估算值加0.5个步长所在频点的幅值X_w(l_m+0.5)大小，三者中的最大值设为|X'_w(l_m-1)|；然后比较|X'_w(l_m-1)|与当第m次迭代时频率估算值对应的频谱幅值|X_w(l_m)|的大小，若|X'_w(l_m-1)|＜|X_w(l_m)|，则为|X_w(l_m)|为最大频谱幅值，|X'_w(l_m-1)|为第二大频谱幅值；若|X'_w(l_m-1)|≥|X_w(l_m)|，|X'_w(l_m-1)|最大频谱幅值，|X'_w(l_m-1)|为第二大频谱幅值；

即当|X'_w(l_m-1)|＜|X_w(l_m)|时，

当|X'_w(l_m-1)|≥|X_w(l_m)|时，按照公式l_m+1＝l_m+δ_m计算下一次迭代时频率估算值，当满足|l_m+1-l_m|＜τ，此时得到的第m+1个频率估算值为真实频率值λ₀；

S6：将步骤S5得到的真实频率值λ₀带入公式计算得到信号的频率f_o。

计算系统误差：

设置采样信号x(n)的幅值A₀设为1，采样频率f_s为1024Hz，数据总长度M为1024点，故频率分辨率为1Hz，对每个窗函数的离散傅里叶变换谱线k，在区间[k–0.5，k+0.5)内以步距0.05进行频率扫描。对于每一个步距频率，相位在区间(-π，π]范围内变化，步距为π/72。即对于每个窗函数的离散傅里叶变换谱线产生2880个测试信号，其中与理论信号的最大绝对误差作为该谱线k的系统误差。

从图2可以看出，不同窗函数在该方法下计算出的信号的频率f_o的系统误差，其中窗函数包括矩形窗(RectangleWindow)、海明窗(HammingWindow)、三角窗(TriangleWindow)、汉宁窗(HanningWindow)、布莱克曼窗(BlackmanWindow),从图2可以看出，λ₀满足λ₀＞5和λ₀＜N/2-5系统误差是非常小的，可被忽略，对于所有窗函数，最大频率误差约为10^-3，这足以满足大多数工程应用，同时表明该算法可适用于不同的窗函数。

检测随机噪声对不同信号频率算法的影响：

设置采样信号x(n)的幅值A0设为1，采样频率f_s为1024Hz，数据总长度M为1024点，故频率分辨率为1Hz，为了使得谱线定位错误的几率较高，将信噪比设为–5dB，为了充分验证谱线定位错误对本算法的影响仿真信号的频率在[255.5，256.5]区间内逐渐增加，步距为0.025，仿真信号的相位在[-π,π)范围内随机选取。每个频率均独立产生50000个带有加性高斯白噪声干扰的实例，观察不同算法的方差；从图3可以看出，无论是加矩形窗还是汉宁窗，本算法(ACIIpDFT)都具有最小的方差，且方差不随频率偏差δ改变而改变，具有较好的一致性。因此，本算法不但具有较好的噪声性且不会受到谱线定位错误的影响。

Claims (5)

1.一种用于油井动液面检测的频率估算方法，其特征在于包括以下步骤：

S1：采集油井动液面的声场信号，得到M点长度的采样信号x(n)：

其中A_o表示幅值，f_o表示信号的频率，f_s表示采样频率，表示相位，n表示采样点序数，M为对信号的总的采样点数；

S2：对步骤S1得到的采样信号x(n)进行加窗处理；

构造K点长度的窗函数w(n)，对采样信号x(n)进行加窗，得到函数：x_w(n)＝x(n)w(n)；

S3：将步骤S2得到的加窗函数x_w(n)进行离散傅里叶变换，得到频谱X_w(k)；

S4：从步骤S3所得的频谱X_w(k)中寻找幅值最大的频点以及当前最大频谱幅值X_w(l)和第二大频谱幅值X_w(l±1)；其中幅值最大的频点记为第一次迭代时的频率初始值l₁；

S5：计算真实频率值λ₀；

根据公式l_m+1＝l_m+δ_m计算第m+1次迭代时频率估算值，且第m+1次迭代时频率估算值和第m个频率估算值偏差其中：q为矩形窗与选择的窗函数的扇形损失比，|X_w(l_m)|为第m次迭代时频率估算值对应的频谱幅值，|X'_w(l_m-1)|为|X_w(l_m-1)|、|X_w(l_m-0.5)|以及|X_w(l_m+0.5)|之间的最大值，X_w(l_m-1)为第m-1次迭代时频率估算值对应的频谱幅值，X_w(l_m-0.5)为第m次迭代时频率估算值减0.5个步长所在频点的幅值，X_w(l_m+0.5)为第m次迭代时频率估算值加0.5个步长所在频点的幅值，当|l_m+1-l_m|＜τ时，τ为预设阈值，迭代结束，并设定真实频率值λ₀＝l_m+1；

S6：将步骤S5得到的真实频率值λ₀带入公式计算得到信号的频率f_o。

2.根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于：当所加窗函数为矩形窗时，对其进行离散傅里叶变换为得到：

$W\left(k\right)=\frac{sin\mathrm{\π}k}{\mathrm{\π}k}\·e^{-j\mathrm{\π}k}.$

3.根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于：当所加窗函数为三角窗时，对其进行离散傅里叶变换得到：

$W\left(k\right)=\frac{2}{N}{e}^{-j\frac{N-2}{N}\·k\mathrm{\π}}{\[\frac{sin\frac{k\mathrm{\π}}{2}}{sin\frac{k\mathrm{\π}}{N}}\]}^2.$

4.根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于：当所加窗函数为汉宁窗时，对其进行离散傅里叶变换得到：

$W\left(k\right)=sin\mathrm{\π}k\·e^{-jk\mathrm{\π}}\[\frac{1}{4N}\frac{sin\frac{2\mathrm{\π}k}{N}}{{\left(sin\frac{\mathrm{\π}k}{N}\right)}^2}-\frac{1}{4N}\frac{sin\frac{2\mathrm{\π}k}{N}}{sin\frac{\mathrm{\π}(k-1)}{N}\·sin\frac{\mathrm{\π}(k+1)}{N}}\].$

5.根据权利要求1所述用于油井动液面检测的频率估算方法，其特征在于：步骤S5中矩形窗与选择的窗函数的扇形损失比q为：其中SL_r矩形窗的扇形损失，SL_c为选择的窗函数的扇形损失，且窗函数的扇形损失为：其中W(0.5)为信号最差非整周期采样的情形时候的最大谱线幅值，W(0)为信号整周期采样时最大谱线幅值。