CN111535802B - 泥浆脉冲信号处理方法 - Google Patents
泥浆脉冲信号处理方法 Download PDFInfo
- Publication number
- CN111535802B CN111535802B CN202010381133.3A CN202010381133A CN111535802B CN 111535802 B CN111535802 B CN 111535802B CN 202010381133 A CN202010381133 A CN 202010381133A CN 111535802 B CN111535802 B CN 111535802B
- Authority
- CN
- China
- Prior art keywords
- pulse signal
- mud pulse
- value
- noise
- pump noise
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Measuring Volume Flow (AREA)
Abstract
本发明涉及一种泥浆脉冲信号处理方法,在充分依据随钻泥浆脉冲信号产生、传输及加噪特性的基础上,将单压力传感器采集的泥浆脉冲信号依次进行泵噪声状态空间模型构建、卡尔曼滤波、去除随机噪声、去除基线漂移和设置自适应阈值处理,获得明显的有效传输脉冲信号。本发明方法简单可靠,能够快速、准确、可靠识别有效脉冲信号,实现井下数据高效传输,去噪效果好,具有很大的实用价值。
Description
技术领域
本发明属于石油天然气勘探技术领域,涉及石油天然气随钻测量技术,具体地说,涉及了一种基于单压力传感器的泥浆脉冲信号处理方法。
背景技术
泥浆脉冲传输***以其低成本和可靠性作为国内外最常用的随钻信号传输***,仍具有数据传输效率低和信号提取困难等问题。稳定且高传输效率的脉冲编码技术与良好的脉冲信号处理算法是随钻传输***的前沿研究热点。泥浆脉冲信号在传输过程中会受到泥浆泵噪声、反射噪声和一系列随机噪声(如:井底机械振动、钻柱失稳、泥浆摩擦等引起的噪声)的影响,当其传至地面时,已经极其微弱,信号提取和处理工作变得十分困难,因此如何从幅值较大,带宽宽广的复杂噪声背景中检测出有用信号成为了泥浆脉冲传输技术中急需解决的关键问题。特别是,由于泥浆泵噪声幅值大、谐波多、频带宽,对有效脉冲信号的干扰尤为剧烈,研究切实可行的泵噪声抑制方法显得十分必要。
公开号为CN 102900430 A的中国专利公开了一种钻井液连续压力波信号的泵压干扰消除方法,采用两个压力传感器安装在直管道上并相距一段距离,直管道一端与钻井液泵连接,另一端与井口连接,反映井下随钻测量信息的井下钻井液压力信号由井口通过直管道传输到地面并被安装在直管道上的压力传感器接收,将两个压力传感器的测量信号相减得到延迟差动检测信号,由于泵压力干扰的传输方向与井下钻井液压力信号传输方向相反,因此钻井液泵产生的压力干扰在这一过程中被消除,延迟差动检测信号中已无泵压干扰成分,而延迟差动检测信号中包含的井下钻井液压力信号通过基于时域差分方程或基于傅里叶变换的信号重构方法得以恢复,从而达到消除信号中的泵压干扰影响,提高信号信噪比的目的。该专利采用双压力传感器,具有积极意义,但在高压管道上多一个传感器势必增加安全隐患。此外,该专利研究对象仅限于仿真波形,无法得知其实际效果。
公开号为CN 107083957 A的中国专利申请公开了一种钻井液随钻信号的泵冲干扰消除方法及***,该专利申请利用信号稀疏分离方法对随钻信号中的泵冲干扰进行分离,以获得目标脉冲信号。该专利申请针对泵冲干扰和目标脉冲信号在离散余弦变换(DCT)域和时域上的稀疏性不同的特征,通过适当的分解方法,将泵冲干扰和脉冲信号分离出来,从而达到消除泵冲干扰的目的,但不能准确反应传输过程中噪声的实时波动情况,去噪的实时性相对较差。
公开号为CN 104265278 A的中国专利公开了一种利用回音抵消技术消除随钻测井中的泵冲噪声的方法,该专利公开了一种采用自适应滤波技术,将泵冲传感器采集的信号作为自适应滤波器参考值,泥浆泵噪声作为期望值,并将滤波器输出与原始信号反向相加以消除泵冲噪声对信号的影响。该专利在常规单压力传感器的泥浆脉冲检测的基础上,增加了一个泵冲传感器来测量泥浆泵工作状态,有其特别之处,但其实际效果究竟增强多少,却不得而知;此外,该案并未说明如何即时获得泥浆泵噪声以实现实时地自适应滤波,实时性相对较差。
公开号为CN 105041303 A的中国专利公开了一种钻井液随钻数据传输***的泵冲干扰信号消除方法,该专利通过获取泵冲干扰信号基频,并基于泵冲干扰信号基频,采用自适应傅立叶级数合成的方式合成泵冲干扰信号和对接收的井下数据信号同步,在此基础上从经同步后的接收的井下数据信号中消除合成后的所述泵冲干扰信号。该案所提方法可有效去除泵冲干扰信号,但若泵冲频率与采样点处的泵噪声频率存在细微差别时,其去噪效果可能会变差。
由上可知,现有去除泵冲干扰信号的方法,虽然可以对泵冲噪声进行消除,但仍存在实时性和去噪效果相对较差等问题。
发明内容
本发明针对现有技术存在的上述问题,提供一种泥浆脉冲信号处理方法,能够有效快速、准确、可靠识别脉冲信号,实现井下数据高效传输。
为了达到上述目的,本发明提供了一种泥浆脉冲信号处理方法,其具体步骤为:
根据泵噪声通用模型构建泵噪声状态空间模型;
采用卡尔曼滤波算法对泵噪声状态空间模型进行计算,获得泵噪声状态向量估计值通过泵噪声状态向量估计值对采集泥浆脉冲信号A(n)中的泵噪声信号进行重构,获得重构泵噪声信号B(n),则去除泵噪声后的泥浆脉冲信号C(n)=A(n)-B(n);
对泥浆脉冲信号C(n)进行随机噪声去除,得到滤波后泥浆脉冲信号D(n);
对泥浆脉冲信号D(n)去除基线漂移得到泥浆脉冲信号E(n);
按设定时间间隔分段对泥浆脉冲信号E(n)进行自适应阈值划分,获得各时间段的自适应阈值T,将泥浆脉冲信号E(n)中小于相应时间段自适应阈值T的采样点值置零,大于相应时间段自适应阈值T的采样点值不变,得到阈值划分后的泥浆脉冲信号F(n)。
优选的,依据泥浆泵工作机理,所述泵噪声通用模型为:
式中,p(n)为泵噪声,n为连续采样点,取值为1,2,…,N,N为采样点总个数;k为谐波阶次,取值为1,2,…,K,K为谐波总个数;Mk(n)为n时刻第k次谐波的幅值;θk(n)为n时刻第k次谐波的相位;ω是基波角频率,kω是k次谐波的角频率,ω=2πf0/fs,f0为基波频率,fs为采样频率。
优选的,针对不同的变量参数,根据泵噪声通用模型构建时不变线性模型、时变线性模型和非线性模型三种泵噪声状态空间模型;其中:
构建时不变线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n)x2(n)…x2K-1(n)x2K(n)]T,其中各元素表示为:
则时不变线性模型的状态方程表示为:
时不变线性模型的测量方程表示为:
z(n)=Hx(n)+v(n)=[1 0 … 1 0]1×2Kx(n)+v(n) (4)式中,H表示观测转移矩阵,为时不变的常值矩阵;v(n)为测量噪声;z(n)为泥浆脉冲传输***的观测值;
构建时变线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n) x2(n)…x2K-1(n)x2K(n)]T,其中各元素表示为:
则时变线性模型的状态方程表示为:
式中,I为2阶单位矩阵;
时变线性模型的测量方程表示为:
式中,H(n)表示观测转移矩阵,为时变矩阵;
构建非线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n) x2(n)…x2K-1(n)x2K(n)]T,其中各元素表示为:
则非线性模型的状态方程表示为:
x(n+1)=x(n)+w(n) (9)
非线性模型的测量方程表示为:
式中,h(*)表示非线性模型测量方程的非线性函数。
P′(n)=FP(n-1)FT+Q(n-1) (12)
K(n)=P′(n)HT(n)[H(n)P′(n)HT(n)+R(n)]-1 (13)
P(n)=[I-K(n)H(n)]P′(n) (15)
式中,Q(n-1)为过程噪声w(n)的协方差矩阵;R(n)为测量噪声v(n)的协方差矩阵;为一步的预测状态向量;为上一时刻的状态向量估计值;P′(n)为一步的预测均方误差;P(n-1)为上一时刻均方误差估计值;K(n)为卡尔曼增益;为当前时刻泵噪声的最佳状态向量估计值;P(n)为当前时刻的均方误差估计值;
采用时不变线性模型获得的重构泵噪声信号B(n)表示为:
采用时变线性模型获得的重构泵噪声信号B(n)表示为:
令
则非线性模型的测量方程改写为:
z(n)=H′(n)x(n)+γ2(n)+v(n) (21)
线性化后的观测转移矩阵H′(n)为h(*)的雅克比矩阵,表示为:
式中,χi′(n)为n时刻的σ采样点值,r=1,2,…,2K;λ=α2(2K+κ)-2K;α为σ采样点散布程度因数,α=10-4~1;κ为缩放参数,保证为半正定矩阵,κ=0~(3-2K);为χi′(n)的均值权重,为χi′(n)的方差权重;β为σ采样点值的分布因数;
zi′(n)=h(n,χi′(n)) (25)
滤波更新为:
通过公式(23)对泵噪声信号进行重构,获得重构泵噪声信号B(n)。
优选的,采用FIR数字低通滤波器进行随机噪声去除,所述FIR数字低通滤波器的截止频率大于井下脉冲发生器所产生脉冲的频率,所述FIR数字低通滤波器采用窗函数法进行随机噪声去除。
优选的,采用小波阈值去噪法进行随机噪声去除,其具体步骤为:
选择小波基,并依据信号采样率和有效正泥浆脉冲信号所在频段确定分解的层数U;
对泥浆脉冲信号C(n)进行变换获得一组小波系数;
选择阈值和阈值函数,对小波系数进行阈值划分,并依据阈值函数对小波系数进行处理;
将处理后的小波系数进行小波反变换,从而获得滤波后泥浆脉冲信号D(n)。
优选的,采用中值滤波法、或最小二乘拟合算法、或小波变换法对泥浆脉冲信号D(n)去除基线漂移得到泥浆脉冲信号E(n),其中:
采用中值滤波法去除基线漂移的步骤为:设定一个中值滤波窗口长度y,将泥浆脉冲信号D(n)前后各扩展y/2个数得到泥浆脉冲信号D1(m);对泥浆脉冲信号D1(m)进行中值滤波medfilt(D1(m),y)得到泥浆脉冲信号D2(m);删除泥浆脉冲信号D2(m)的前后各y/2个数得到泥浆脉冲信号D3(n);将泥浆脉冲信号D3(n)作为基线,并在泥浆脉冲信号D(n)中减去泥浆脉冲信号D3(n),得到去除基线漂移后的泥浆脉冲信号E(n);
采用最小二乘拟合算法去除基线漂移的步骤为:泥浆脉冲信号D(n)按设定时间间隔分段分别进行多项式最小二乘拟合,获得各时间段的拟合曲线;将各时间段拟合曲线形成的数组作为基线,并在泥浆脉冲信号D(n)中减去,得到去除基线漂移后的泥浆脉冲信号E(n);
采用小波变换法去除基线漂移的步骤为:选定小波基和分解层数对泥浆脉冲信号D(n)进行小波变换,将低频段小波系数置零,再对处理后小波系数进行小波反变换即可得到去除基线漂移后的泥浆脉冲信号E(n)。
优选的,采用最大类间方差法进行自适应阈值划分,其具体步骤为:将一段离散信号按幅值大小分为0~(L-1)级,L为幅值的个数;设幅值为i的采样点数为n,则总采样点数为各等级幅值的概率为pi=ni/N;用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类,即e0={0,1,…,T},e1={T+1,T+2,…,L-1},则泥浆脉冲信号E(n)的幅值的均值μ为:
e0的均值μ0和e1的均值μ1分别为:
与现有技术相比,本发明的优点和积极效果在于:
本发明提供的泥浆脉冲信号处理方法,在充分依据随钻泥浆脉冲信号产生、传输及加噪特性的基础上,将单压力传感器采集的泥浆脉冲信号依次进行泵噪声状态空间模型构建、卡尔曼滤波、去除随机噪声、去除基线漂移和设置自适应阈值处理,获得明显的有效传输脉冲信号。本发明方法简单可靠,能够快速、准确、可靠识别有效脉冲信号,实现井下数据高效传输,去噪效果好,具有很大的实用价值。
附图说明
图1为本发明实施例所述泥浆脉冲信号处理方法的流程图;
图2为本发明实施例含泥浆脉冲信号的某原始压力信号示意图;
图3为本发明实施例原始泥浆压力数据频谱图;
图4为本发明实施例重构泵噪声信号示意图;
图5为本发明实施例去除重构泵噪声后的泥浆脉冲信号和去除随机噪声后泥浆脉冲信号示意图;
图6为本发明实施例去除基线漂移后的泥浆脉冲信号示意图;
图7为本发明实施例自适应阈值示意图;
图8为本发明实施例低于阈值数据置零后的泥浆脉冲信号示意图。
图中,A(n)、含泥浆脉冲信号的原始压力信号,FA(n)、原始泥浆压力信号频谱,B(n)、重构泵噪声信号,C(n)、原始泥浆压力去除重构泵噪声后信号,D(n)、去除随机噪声后的泥浆脉冲信号,E(n)、去除基漂后的泥浆脉冲信号E(n),T、自适应阈值,F(n)、低于阈值数据置零后的泥浆脉冲信号。
具体实施方式
下面,通过示例性的实施方式对本发明进行具体描述。然而应当理解,在没有进一步叙述的情况下,一个实施方式中的元件、结构和特征也可以有益地结合到其他实施方式中。
参见图1,本发明提供了一种泥浆脉冲信号处理方法,其具体步骤为:
S1、针对不同的变量参数,根据泵噪声通用模型构建时不变线性模型、时变线性模型和非线性模型三种泵噪声状态空间模型。
具体地,依据泥浆泵工作机理,所述泵噪声通用模型为:
式中,p(n)为泵噪声,n为连续采样点,取值为1,2,…,N,N为采样点总个数;k为谐波阶次,取值为1,2,…,K,K为谐波总个数;Mk(n)为n时刻第k次谐波的幅值;θk(n)为n时刻第k次谐波的相位;ω是基波角频率,kω是k次谐波的角频率,ω=2πf0/fs,f0为基波频率,fs为采样频率。
具体地,构建时不变线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n)x2(n)…x2K-1(n)x2K(n)]T,其中各元素表示为:
则时不变线性模型的状态方程表示为:
时不变线性模型的测量方程表示为:
z(n)=Hx(n)+v(n)=[1 0 … 1 0]1×2Kx(n)+v(n) (4)(4)
式中,H表示观测转移矩阵,为时不变的常值矩阵;v(n)为测量噪声;z(n)为泥浆脉冲传输***的观测值;
具体地,构建时变线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n) x2(n)…x2K-1(n)x2K(n)]T,其中各元素表示为:
则时变线性模型的状态方程表示为:
式中,I为2阶单位矩阵;
时变线性模型的测量方程表示为:
式中,H(n)表示观测转移矩阵,为时变矩阵;
具体地,构建非线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n)x2(n)…x2K-1(n)x2K(n)]T,其中各元素表示为:
则非线性模型的状态方程表示为:
x(n+1)=x(n)+w(n) (9)
非线性模型的测量方程表示为:
式中,h(*)表示非线性模型测量方程的非线性函数。
S2、采用卡尔曼滤波算法对泵噪声状态空间模型进行计算,获得泵噪声状态向量估计值通过泵噪声状态向量估计值对采集泥浆脉冲信号A(n)中的泵噪声信号进行重构,获得重构泵噪声信号B(n),则去除泵噪声后的泥浆脉冲信号C(n)=A(n)-B(n)。卡尔曼滤波算法包括面向线性模型的标准卡尔曼滤波算法,面向非线性模型的扩展卡尔曼算法和无迹卡尔曼滤波算法。泥浆脉冲信号C(n)包括随机噪声W(n)和有效脉冲信号S(n)。
P′(n)=FP(n-1)FT+Q(n-1) (12)
K(n)=P′(n)HT(n)[H(n)P′(n)HT(n)+R(n)]-1 (13)
P(n)=[I-K(n)H(n)]P′(n) (15)
式中,Q(n-1)为过程噪声w(n)的协方差矩阵;R(n)为测量噪声v(n)的协方差矩阵;为一步的预测状态向量,即状态向量的一步预测值,也就是状态向量的先验估计;为上一时刻的状态向量估计值;P′(n)为一步的预测均方误差;P(n-1)为上一时刻均方误差估计值,即均方误差的一步预测值,也就是均方误差的先验估计;K(n)为卡尔曼增益;为当前时刻泵噪声的最佳状态向量估计值;P(n)为当前时刻的均方误差估计值。
采用时不变线性模型获得的重构泵噪声信号B(n)表示为:
采用时变线性模型获得的重构泵噪声信号B(n)表示为:
令
则非线性模型的测量方程改写为:
z(n)=H′(n)x(n)+γ2(n)+v(n) (21)
线性化后的观测转移矩阵H′(n)为h(*)的雅克比矩阵,表示为:
式中,χi′(n)为n时刻的σ采样点值,r=1,2,…,2K;λ=α2(2K+κ)-2K;α为σ采样点散布程度因数,α=10-4~1;κ为缩放参数,保证为半正定矩阵,κ=0~(3-2K);为χi′(n)的均值权重,为χi′(n)的方差权重;β为σ采样点值的分布因数。需要说明的是,对于高斯分布的状态向量,β=2最佳。
zi′(n)=h(n,χi′(n)) (25)
滤波更新为:
通过公式(23)对泵噪声信号进行重构,获得重构泵噪声信号B(n)。
S3、对泥浆脉冲信号C(n)进行随机噪声去除,得到滤波后泥浆脉冲信号D(n)。
在一具体实施方式中,采用FIR数字低通滤波器进行随机噪声去除,所述FIR数字低通滤波器的截止频率大于井下脉冲发生器所产生脉冲的频率,所述FIR数字低通滤波器采用窗函数法进行随机噪声去除。具体地,FIR数字低通滤波器采用矩形窗、三角窗、汉宁窗、海明窗、布莱克曼窗和凯泽窗中的任意一种窗函数法进行随机噪声去除。需要说明的是,所述FIR数字低通滤波器的窗函数和阶数随实际信号可灵活调整,以得到较佳的消噪信号输出。
在另一具体实施方式中,采用小波阈值去噪法进行随机噪声去除,其具体步骤为:
选择小波基,并依据信号采样率和有效正泥浆脉冲信号所在频段确定分解的层数U;
对泥浆脉冲信号C(n)进行变换获得一组小波系数;
选择阈值和阈值函数,对小波系数进行阈值划分,并依据阈值函数对小波系数进行处理;
利用小波第U级的低频系数和第1级至第U级的高频系数的特性,将处理后的小波系数进行小波反变换,从而获得滤波后泥浆脉冲信号D(n)。
具体地,小波阈值去噪法可选用软阈值函数、硬阈值函数、双曲线阈值函数、指数阈值函数和幂阈值函数等。
需要说明的是,所述小波阈值去噪法的小波基与分解层数随实际信号可灵活调整,以得到较佳的消噪信号输出。
S4、对泥浆脉冲信号D(n)去除基线漂移得到泥浆脉冲信号E(n)。
在一具体实时方式中,采用中值滤波法去除基线漂移,其具体步骤为:设定一个中值滤波窗口长度y,将泥浆脉冲信号D(n)前后各扩展y/2个数得到泥浆脉冲信号D1(m);对泥浆脉冲信号D1(m)进行中值滤波medfilt(D1(m),y)得到泥浆脉冲信号D2(m);删除泥浆脉冲信号D2(m)的前后各y/2个数得到泥浆脉冲信号D3(n);将泥浆脉冲信号D3(n)作为基线,并在泥浆脉冲信号D(n)中减去泥浆脉冲信号D3(n),得到去除基线漂移后的泥浆脉冲信号E(n)。
具体地,所述中值滤波窗口长度为立管压力信号采样频率的3-5倍,且为偶数;所述数组D1(m)的前y/2个数为0,后y/2个数与D(n)的最后一个数相等。
在另一具体实施方式中,采用最小二乘拟合算法去除基线漂移的步骤为:泥浆脉冲信号D(n)按设定时间间隔分段分别进行多项式最小二乘拟合,获得各时间段的拟合曲线;将各时间段拟合曲线形成的数组作为基线,并在泥浆脉冲信号D(n)中减去,得到去除基线漂移后的泥浆脉冲信号E(n)。
具体地,所述时间间隔通常为5-20秒,最小二乘拟合多项式为一元三次、或四次、或五次多项式,具体根据实际情况进行选择。
在又一具体实施方式中,采用小波变换法去除基线漂移的步骤为:选定小波基和分解层数对泥浆脉冲信号D(n)进行小波变换,将低频段小波系数置零,再对处理后小波系数进行小波反变换即可得到去除基线漂移后的泥浆脉冲信号E(n)。
具体地,所述小波变换算法的小波基与分解层数随实际信号的采样率和有效正脉冲所在频段进行调整,分解后的低频系数中不包含有效正脉冲成分。
S5、按设定时间间隔分段对泥浆脉冲信号E(n)进行自适应阈值划分,获得各时间段的自适应阈值T,将泥浆脉冲信号E(n)中小于相应时间段自适应阈值T的采样点值置零,大于相应时间段自适应阈值T的采样点值不变,得到阈值划分后的泥浆脉冲信号F(n)。需要说明的是,所述时间间隔可依据泥浆脉冲信号E(n)波形质量和处理时间限制灵活调整。
具体地,采用最大类间方差法进行自适应阈值划分,其具体步骤为:将一段离散信号按幅值大小分为0~(L-1)级,L为幅值的个数;设幅值为i的采样点数为n,则总采样点数为各等级幅值的概率为pi=ni/N;用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类,即e0={0,1,…,T},e1={T+1,T+2,…,L-1},则泥浆脉冲信号E(n)的幅值的均值μ为:
e0的均值μ0和e1的均值μ1分别为:
本发明上述方法在充分依据随钻泥浆脉冲信号产生、传输及加噪特性的基础上,将单压力传感器采集的泥浆压力信号依次进行泵噪声状态空间模型构建、卡尔曼滤波、抑制随机噪声、消除基线漂移和设置自适应阈值等处理,可以获得明显的有效传输脉冲信号。该方法简便可靠,具有很大的实用价值。
以下结合具体地实施例对本发明上述方法作出进一步说明。
参见图2,图2所示为地面立管处连续采集的某段100s含泥浆脉冲信号的原始压力信号A(n),首先对含泥浆脉冲信号的原始压力信号A(n)进行FFT变换,获得原始泥浆压力信号频谱FA(n)(如图3)。根据该频谱特征采用本发明所述泥浆脉冲从处理方法对上述原始压力信号A(n)进行处理,得到波形特征明显的泥浆脉冲信号,用于实现解码及后续数据处理。其具体步骤为:
S1、依据泥浆泵工作机理,泵噪声通用模型表示为:
式中,p(n)为泵噪声,n为连续采样点,取值为1,2,…,N,N为采样点总个数;k为谐波阶次,取值为1,2,…,K,K为谐波总个数;Mk(n)为n时刻第k次谐波的幅值;θk(n)为n时刻第k次谐波的相位;ω是基波角频率,kω是k次谐波的角频率,ω=2πf0/fs,f0为基波频率,fs为采样频率。
根据泵噪声通用模型构造泵噪声状态空间模型,其步骤为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n)x2(n)…x2K-1(n)x2K(n)]T,其中各元素表示为:
则泵噪声状态空间模型的状态方程表示为:
泵噪声状态空间模型的测量方程表示为:
z(n)=Hx(n)+v(n)=[1 0 … 1 0]1×2Kx(n)+v(n) (4)(4)
式中,H表示观测转移矩阵,为时不变的常值矩阵;v(n)为测量噪声;z(n)为泥浆脉冲传输***的观测值
本实施例中,谐波阶次K为12,基波频率f0为1.664Hz,采样频率fs为200Hz;泵噪声状态空间模型为线性时不变模型。
P′(n)=FP(n-1)FT+Q(n-1) (12)
K(n)=P′(n)HT(n)[H(n)P′(n)HT(n)+R(n)]-1 (13)
P(n)=[I-K(n)H(n)]P′(n) (15)
式中,Q(n-1)为过程噪声w(n)的协方差矩阵,R(n)为测量噪声v(n)的协方差矩阵,为一步的预测状态向量;为上一时刻的状态向量估计值;P′(n)为一步的预测均方误差;P(n-1)为上一时刻均方误差估计值;K(n)为卡尔曼增益;为当前时刻泵噪声的最佳状态向量估计值;P(n)为当前时刻的均方误差估计值。
从含泥浆脉冲信号的原始压力信号A(n)中去除重构泵噪声信号B(n),得到原始泥浆压力去除重构泵噪声后泥浆脉冲信号C(n),C(n)=A(n)-B(n),如图5所示。其中,泥浆脉冲信号C(n)包括随机噪声W(n)和有效脉冲信号S(n)。
S3、采用小波阈值去噪法对泥浆脉冲信号C(n)进行随机噪声去除,得到滤波后泥浆脉冲信号D(n)。其具体步骤为:
选择小波基为Sym 8,并依据信号采样率和有效正泥浆脉冲信号所在频段确定分解的层数为7层;
对泥浆脉冲信号C(n)进行变换获得一组小波系数;
对小波系数进行阈值划分,并依据幂阈值函数对小波系数进行处理;
利用小波第7级的低频系数和第1级至第7级的高频系数的特性,将处理后的小波系数进行小波反变换,从而获得滤波后泥浆脉冲信号D(n),如图5所示。
S4、采用小波变换法对泥浆脉冲信号D(n)去除基线漂移得到泥浆脉冲信号E(n)。其具体步骤为:
设定中值滤波窗口长度y为900,将泥浆脉冲信号D(n)前后各扩展y/2个数得到泥浆脉冲信号D1(m);然后对泥浆脉冲信号D1(m)进行中值滤波medfilt(D1(m),y)得到泥浆脉冲信号D2(m);最后删除泥浆脉冲信号D2(m)的前后各y/2个数得到泥浆脉冲信号D3(n),将泥浆脉冲信号D3(n)作为基线,并在泥浆脉冲信号D3(n)D(n)中减去泥浆脉冲信号D3(n),得到去除基漂后的泥浆脉冲信号E(n),如图6所示。
本实施例中,立管压力信号采样频率为200Hz,所述泥浆脉冲信号D1(m)的前450个数为0,后450个数与泥浆脉冲信号D(n)的最后一个数相等。
S5、对于去除基线漂移后的泥浆脉冲信号E(n),按时间间隔10s分段采用最大类间方差法进行自适应阈值划分,获得各时间段的自适应划分阈值T,如图7所示。
采用最大类间方差法进行自适应阈值划分具体步骤为:
将一段离散信号按幅值大小分为0~(L-1)级,L为幅值的个数;设幅值为i的采样点数为n,则总采样点数为各等级幅值的概率为pi=ni/N;用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类,即e0={0,1,…,T},e1={T+1,T+2,…,L-1},则泥浆脉冲信号E(n)的幅值的均值μ为:
e0的均值μ0为:
将泥浆脉冲信号E(n)中小于相应时间段自适应阈值T的采样点值置零,得到泥浆脉冲信号信号F(n),如图8所示。
上述实施例用来解释本发明,而不是对本发明进行限制,在本发明的精神和权利要求的保护范围内,对本发明做出的任何修改和改变,都落入本发明的保护范围。
Claims (8)
1.一种泥浆脉冲信号处理方法,其特征在于,其具体步骤为:
根据泵噪声通用模型构建泵噪声状态空间模型;依据泥浆泵工作机理,所述泵噪声通用模型为:
式中,p(n)为泵噪声,n为连续采样点,取值为1,2,…,N,N为采样点总个数;k为谐波阶次,取值为1,2,…,K,K为谐波总个数;
Mk(n)为n时刻第k次谐波的幅值;θk(n)为n时刻第k次谐波的相位;
ω是基波角频率,kω是k次谐波的角频率,ω=2πf0/fs,f0为基波频率,fs为采样频率;
针对不同的变量参数,根据泵噪声通用模型构建时不变线性模型、时变线性模型和非线性模型三种泵噪声状态空间模型;其中:
构建时不变线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n)x2(n)…x2K-1(n)
x2K(n)]T,其中各元素表示为:
则时不变线性模型的状态方程表示为:
z(n)=Hx(n)+v(n)=[1 0 … 1 0]1×2Kx(n)+v(n) (4)
式中,H表示观测转移矩阵,为时不变的常值矩阵;v(n)为测量噪声;z(n)为泥浆脉冲传输***的观测值;
构建时变线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n) x2(n) … x2K-1(n) x2K(n)]T,其中各元素表示为:
则时变线性模型的状态方程表示为:
式中,I为2阶单位矩阵;
时变线性模型的测量方程表示为:
式中,H(n)表示观测转移矩阵,为时变矩阵;
构建非线性模型的方法为:
设含K阶次谐波的泵噪声的状态向量为x(n)=[x1(n) x2(n) … x2K-1(n) x2K(n)]T,其中各元素表示为:
则非线性模型的状态方程表示为:
x(n+1)=x(n)+w(n) (9)
非线性模型的测量方程表示为:
式中,h(*)表示非线性模型测量方程的非线性函数;
采用卡尔曼滤波算法对泵噪声状态空间模型进行计算,获得泵噪声状态向量估计值通过泵噪声状态向量估计值对采集泥浆脉冲信号A(n)中的泵噪声信号进行重构,获得重构泵噪声信号B(n),则去除泵噪声后的泥浆脉冲信号C(n)=A(n)-B(n);
对泥浆脉冲信号C(n)进行随机噪声去除,得到滤波后泥浆脉冲信号D(n);
对泥浆脉冲信号D(n)去除基线漂移得到泥浆脉冲信号E(n);
按设定时间间隔分段对泥浆脉冲信号E(n)进行自适应阈值划分,获得各时间段的自适应阈值T,将泥浆脉冲信号E(n)中小于相应时间段自适应阈值T的采样点值置零,大于相应时间段自适应阈值T的采样点值不变,得到阈值划分后的泥浆脉冲信号F(n)。
2.如权利要求1所述的泥浆脉冲信号处理方法,其特征在于,采用标准卡尔曼滤波算法计算泵噪声状态向量估计值其具体方法为:将时不变线性模型或时变线性模型的状态转移矩阵F和观测转移矩阵H带入以下公式循环递推,即可获得泵噪声状态向量估计
P′(n)=FP(n-1)FT+Q(n-1) (12)
K(n)=P′(n)HT(n)[H(n)P′(n)HT(n)+R(n)]-1 (13)
P(n)=[I-K(n)H(n)]P′(n) (15)
式中,Q(n-1)为过程噪声w(n)的协方差矩阵;R(n)为测量噪声v(n)的协方差矩阵;为一步的预测状态向量;为上一时刻的状态向量估计值;P′(n)为一步的预测均方误差;P(n-1)为上一时刻均方误差估计值;K(n)为卡尔曼增益;为当前时刻泵噪声的最佳状态向量估计值;P(n)为当前时刻的均方误差估计值;
采用时不变线性模型获得的重构泵噪声信号B(n)表示为:
采用时变线性模型获得的重构泵噪声信号B(n)表示为:
式中,χi′(n)为n时刻的σ采样点值,r=1,2,…,2K;λ=α2(2K+κ)-2K;
zi′(n)=h(n,χi′(n)) (25)
滤波更新为:
通过公式(23)对泵噪声信号进行重构,获得重构泵噪声信号B(n)。
5.如权利要求1所述的泥浆脉冲信号处理方法,其特征在于,采用FIR数字低通滤波器进行随机噪声去除,所述FIR数字低通滤波器的截止频率大于井下脉冲发生器所产生脉冲的频率,所述FIR数字低通滤波器采用窗函数法进行随机噪声去除。
6.如权利要求1所述的泥浆脉冲信号处理方法,其特征在于,采用小波阈值去噪法进行随机噪声去除,其具体步骤为:
选择小波基,并依据信号采样率和有效正泥浆脉冲信号所在频段确定分解的层数U;
对泥浆脉冲信号C(n)进行变换获得一组小波系数;
选择阈值和阈值函数,对小波系数进行阈值划分,并依据阈值函数对小波系数进行处理;
将处理后的小波系数进行小波反变换,从而获得滤波后泥浆脉冲信号D(n)。
7.如权利要求1所述的泥浆脉冲信号处理方法,其特征在于,采用中值滤波法、或最小二乘拟合算法、或小波变换法对泥浆脉冲信号D(n)去除基线漂移得到泥浆脉冲信号E(n),其中:
采用中值滤波法去除基线漂移的步骤为:设定一个中值滤波窗口长度y,将泥浆脉冲信号D(n)前后各扩展y/2个数得到泥浆脉冲信号D1(m);对泥浆脉冲信号D1(m)进行中值滤波medfilt(D1(m),y)得到泥浆脉冲信号D2(m);删除泥浆脉冲信号D2(m)的前后各y/2个数得到泥浆脉冲信号D3(n);将泥浆脉冲信号D3(n)作为基线,并在泥浆脉冲信号D(n)中减去泥浆脉冲信号D3(n),得到去除基线漂移后的泥浆脉冲信号E(n);
采用最小二乘拟合算法去除基线漂移的步骤为:泥浆脉冲信号D(n)按设定时间间隔分段分别进行多项式最小二乘拟合,获得各时间段的拟合曲线;将各时间段拟合曲线形成的数组作为基线,并在泥浆脉冲信号D(n)中减去,得到去除基线漂移后的泥浆脉冲信号E(n);
采用小波变换法去除基线漂移的步骤为:选定小波基和分解层数对泥浆脉冲信号D(n)进行小波变换,将低频段小波系数置零,再对处理后小波系数进行小波反变换即可得到去除基线漂移后的泥浆脉冲信号E(n)。
8.如权利要求1所述的泥浆脉冲信号处理方法,其特征在于,采用最大类间方差法进行自适应阈值划分,其具体步骤为:将一段离散信号按幅值大小分为0~(L-1)级,L为幅值的个数;设幅值为i的采样点数为n,则总采样点数为各等级幅值的概率为pi=ni/N;用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类,即e0={0,1,…,T},e1={T+1,T+2,…,L-1},则泥浆脉冲信号E(n)的幅值的均值μ为:
e0的均值μ0和e1的均值μ1分别为:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010381133.3A CN111535802B (zh) | 2020-05-08 | 2020-05-08 | 泥浆脉冲信号处理方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010381133.3A CN111535802B (zh) | 2020-05-08 | 2020-05-08 | 泥浆脉冲信号处理方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111535802A CN111535802A (zh) | 2020-08-14 |
CN111535802B true CN111535802B (zh) | 2023-04-14 |
Family
ID=71971679
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010381133.3A Active CN111535802B (zh) | 2020-05-08 | 2020-05-08 | 泥浆脉冲信号处理方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111535802B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP3752214A1 (en) * | 2018-02-16 | 2020-12-23 | Gambro Lundia AB | Filtering a pressure signal from a medical apparatus |
CN112554875A (zh) * | 2020-12-08 | 2021-03-26 | 中国石油天然气集团有限公司 | 随钻泥浆正脉冲信号的处理方法 |
CN112418174A (zh) * | 2020-12-08 | 2021-02-26 | 中国石油天然气集团有限公司 | 随钻泥浆随机噪声的去除方法 |
CN113806945B (zh) * | 2021-07-22 | 2024-03-22 | 中国石油大学(北京) | 控压钻采井控过程中的地层反演方法、装置及存储介质 |
CN116906313B (zh) * | 2023-09-12 | 2023-11-28 | 山东亿昌装配式建筑科技有限公司 | 一种建筑节能给排水智能控制*** |
CN116955941B (zh) * | 2023-09-21 | 2023-12-19 | 中石化经纬有限公司 | 一种随钻测量连续波信号去噪方法 |
CN117888889B (zh) * | 2024-03-14 | 2024-06-07 | 德州联合石油科技股份有限公司 | 泥浆脉冲信号的处理方法及其装置、电子设备 |
Family Cites Families (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US8860583B2 (en) * | 2008-04-03 | 2014-10-14 | Baker Hughes Incorporated | Mud channel characterization over depth |
CN103670380B (zh) * | 2013-12-18 | 2017-01-25 | 贝兹维仪器(苏州)有限公司 | 一种井下泥浆脉冲信号的产生装置 |
CN106158301A (zh) * | 2015-04-10 | 2016-11-23 | 刘会灯 | 配电变压器运行噪声有源降噪*** |
CN105227503B (zh) * | 2015-09-08 | 2019-01-18 | 北京航空航天大学 | 一种基于无线随钻测量***的井下信源信道联合编码方法 |
CN106301289A (zh) * | 2016-08-03 | 2017-01-04 | 广东工业大学 | 利用自适应滤波算法消除泥浆脉冲信号中的泵冲噪声的方法 |
CN106437689B (zh) * | 2016-09-13 | 2019-04-09 | 中国石油大学(华东) | 一种随钻泥浆正脉冲信号的处理方法 |
CN109522802B (zh) * | 2018-10-17 | 2022-05-24 | 浙江大学 | 应用经验模态分解和粒子群优化算法的泵噪消除方法 |
-
2020
- 2020-05-08 CN CN202010381133.3A patent/CN111535802B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN111535802A (zh) | 2020-08-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111535802B (zh) | 泥浆脉冲信号处理方法 | |
Gaci | A new ensemble empirical mode decomposition (EEMD) denoising method for seismic signals | |
Wu et al. | Noise attenuation for 2-D seismic data by radial-trace time-frequency peak filtering | |
CN107083957B (zh) | 钻井液随钻信号的泵冲干扰消除方法及*** | |
CN105738948B (zh) | 一种基于小波变换的微地震数据降噪方法 | |
CN103487835B (zh) | 一种基于模型约束的多分辨率波阻抗反演方法 | |
US20150168573A1 (en) | Geologic quality factor inversion method | |
US20130154845A1 (en) | Mud Pulse Telemetry Noise Reduction Method | |
CN104849756B (zh) | 一种提高地震数据分辨率增强有效弱信号能量的方法 | |
CN112835103B (zh) | 自适应去鬼波与宽频准零相位反褶积联合处理方法及*** | |
CN103645507B (zh) | 地震记录的处理方法 | |
CN109143331B (zh) | 地震子波提取方法 | |
US8942330B2 (en) | Interference reduction method for downhole telemetry systems | |
CN105041303B (zh) | 钻井液随钻数据传输***的泵冲干扰信号消除方法 | |
Chen et al. | Study on model-based pump noise suppression method of mud pulse signal | |
CN108508489B (zh) | 一种基于波形微变化匹配的地震反演方法 | |
Chen et al. | MWD drilling mud signal de-noising and signal extraction research based on the pulse-code information | |
Li et al. | A new comprehensive filtering model for pump shut-in water hammer pressure wave signals during hydraulic fracturing | |
CN113341463B (zh) | 一种叠前地震资料非平稳盲反褶积方法及相关组件 | |
CN110749923A (zh) | 一种基于范数方程提高分辨率的反褶积方法 | |
CN112554875A (zh) | 随钻泥浆正脉冲信号的处理方法 | |
DONG et al. | Fast implementation technique of adaptive Kalman filtering deconvolution via dyadic wavelet transform | |
CN106249284B (zh) | 基于q值差反射的衰减信号分解方法 | |
Li et al. | Random-noise attenuation by an amplitude-preserved time-frequency peak based on empirical wavelet transform predictive filtering | |
CN114428283B (zh) | 薄储层及地震弱反射层的地震识别方法及电子设备 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |