CN111535802B - 泥浆脉冲信号处理方法 - Google Patents

泥浆脉冲信号处理方法 Download PDF

Info

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
Application number
CN202010381133.3A
Other languages
English (en)
Other versions
CN111535802A (zh
Inventor
鄢志丹
陈果
王俊飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN202010381133.3A priority Critical patent/CN111535802B/zh
Publication of CN111535802A publication Critical patent/CN111535802A/zh
Application granted granted Critical
Publication of CN111535802B publication Critical patent/CN111535802B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

本发明涉及一种泥浆脉冲信号处理方法，在充分依据随钻泥浆脉冲信号产生、传输及加噪特性的基础上，将单压力传感器采集的泥浆脉冲信号依次进行泵噪声状态空间模型构建、卡尔曼滤波、去除随机噪声、去除基线漂移和设置自适应阈值处理，获得明显的有效传输脉冲信号。本发明方法简单可靠，能够快速、准确、可靠识别有效脉冲信号，实现井下数据高效传输，去噪效果好，具有很大的实用价值。

Description

泥浆脉冲信号处理方法
技术领域
本发明属于石油天然气勘探技术领域，涉及石油天然气随钻测量技术，具体地说，涉及了一种基于单压力传感器的泥浆脉冲信号处理方法。
背景技术
泥浆脉冲传输系统以其低成本和可靠性作为国内外最常用的随钻信号传输系统，仍具有数据传输效率低和信号提取困难等问题。稳定且高传输效率的脉冲编码技术与良好的脉冲信号处理算法是随钻传输系统的前沿研究热点。泥浆脉冲信号在传输过程中会受到泥浆泵噪声、反射噪声和一系列随机噪声(如：井底机械振动、钻柱失稳、泥浆摩擦等引起的噪声)的影响，当其传至地面时，已经极其微弱，信号提取和处理工作变得十分困难，因此如何从幅值较大，带宽宽广的复杂噪声背景中检测出有用信号成为了泥浆脉冲传输技术中急需解决的关键问题。特别是，由于泥浆泵噪声幅值大、谐波多、频带宽，对有效脉冲信号的干扰尤为剧烈，研究切实可行的泵噪声抑制方法显得十分必要。
公开号为CN 102900430 A的中国专利公开了一种钻井液连续压力波信号的泵压干扰消除方法，采用两个压力传感器安装在直管道上并相距一段距离，直管道一端与钻井液泵连接，另一端与井口连接，反映井下随钻测量信息的井下钻井液压力信号由井口通过直管道传输到地面并被安装在直管道上的压力传感器接收，将两个压力传感器的测量信号相减得到延迟差动检测信号，由于泵压力干扰的传输方向与井下钻井液压力信号传输方向相反，因此钻井液泵产生的压力干扰在这一过程中被消除，延迟差动检测信号中已无泵压干扰成分，而延迟差动检测信号中包含的井下钻井液压力信号通过基于时域差分方程或基于傅里叶变换的信号重构方法得以恢复，从而达到消除信号中的泵压干扰影响，提高信号信噪比的目的。该专利采用双压力传感器，具有积极意义，但在高压管道上多一个传感器势必增加安全隐患。此外，该专利研究对象仅限于仿真波形，无法得知其实际效果。
公开号为CN 107083957 A的中国专利申请公开了一种钻井液随钻信号的泵冲干扰消除方法及系统，该专利申请利用信号稀疏分离方法对随钻信号中的泵冲干扰进行分离，以获得目标脉冲信号。该专利申请针对泵冲干扰和目标脉冲信号在离散余弦变换(DCT)域和时域上的稀疏性不同的特征，通过适当的分解方法，将泵冲干扰和脉冲信号分离出来，从而达到消除泵冲干扰的目的，但不能准确反应传输过程中噪声的实时波动情况，去噪的实时性相对较差。
公开号为CN 104265278 A的中国专利公开了一种利用回音抵消技术消除随钻测井中的泵冲噪声的方法，该专利公开了一种采用自适应滤波技术，将泵冲传感器采集的信号作为自适应滤波器参考值，泥浆泵噪声作为期望值，并将滤波器输出与原始信号反向相加以消除泵冲噪声对信号的影响。该专利在常规单压力传感器的泥浆脉冲检测的基础上，增加了一个泵冲传感器来测量泥浆泵工作状态，有其特别之处，但其实际效果究竟增强多少，却不得而知；此外，该案并未说明如何即时获得泥浆泵噪声以实现实时地自适应滤波，实时性相对较差。
公开号为CN 105041303 A的中国专利公开了一种钻井液随钻数据传输系统的泵冲干扰信号消除方法，该专利通过获取泵冲干扰信号基频，并基于泵冲干扰信号基频，采用自适应傅立叶级数合成的方式合成泵冲干扰信号和对接收的井下数据信号同步，在此基础上从经同步后的接收的井下数据信号中消除合成后的所述泵冲干扰信号。该案所提方法可有效去除泵冲干扰信号，但若泵冲频率与采样点处的泵噪声频率存在细微差别时，其去噪效果可能会变差。
由上可知，现有去除泵冲干扰信号的方法，虽然可以对泵冲噪声进行消除，但仍存在实时性和去噪效果相对较差等问题。
发明内容
本发明针对现有技术存在的上述问题，提供一种泥浆脉冲信号处理方法，能够有效快速、准确、可靠识别脉冲信号，实现井下数据高效传输。
为了达到上述目的，本发明提供了一种泥浆脉冲信号处理方法，其具体步骤为：
根据泵噪声通用模型构建泵噪声状态空间模型；
采用卡尔曼滤波算法对泵噪声状态空间模型进行计算，获得泵噪声状态向量估计值
Figure GDA0004088707120000031
通过泵噪声状态向量估计值
Figure GDA0004088707120000032
对采集泥浆脉冲信号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)。
优选的，依据泥浆泵工作机理，所述泵噪声通用模型为：
Figure GDA0004088707120000033
式中，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，其中各元素表示为:
Figure GDA0004088707120000041
则时不变线性模型的状态方程表示为：
Figure GDA0004088707120000042
式中，F表示状态转移矩阵，为时不变的常值矩阵；w(n)为过程噪声；
Figure GDA0004088707120000043
x(n+1)为泥浆脉冲传输系统的状态值；
时不变线性模型的测量方程表示为：
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，其中各元素表示为:
Figure GDA0004088707120000051
则时变线性模型的状态方程表示为：
Figure GDA0004088707120000052
式中，I为2阶单位矩阵；
时变线性模型的测量方程表示为：
Figure GDA0004088707120000053
式中，H(n)表示观测转移矩阵，为时变矩阵；
构建非线性模型的方法为：
设含K阶次谐波的泵噪声的状态向量为x(n)＝[x1(n) x2(n)…x2K-1(n)x2K(n)]T，其中各元素表示为:
Figure GDA0004088707120000061
则非线性模型的状态方程表示为：
x(n+1)＝x(n)+w(n)                     (9)
非线性模型的测量方程表示为：
Figure GDA0004088707120000062
式中，h(*)表示非线性模型测量方程的非线性函数。
优选的，采用标准卡尔曼滤波算法计算泵噪声状态向量估计值
Figure GDA0004088707120000063
其具体方法为：
将时不变线性模型或时变线性模型的状态转移矩阵F和观测转移矩阵H带入以下公式循环递推，即可获得泵噪声状态向量估计
Figure GDA0004088707120000064
Figure GDA0004088707120000065
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)
Figure GDA0004088707120000066
P(n)＝[I-K(n)H(n)]P′(n)                  (15)
式中，Q(n-1)为过程噪声w(n)的协方差矩阵；R(n)为测量噪声v(n)的协方差矩阵；
Figure GDA0004088707120000067
为一步的预测状态向量；
Figure GDA0004088707120000068
为上一时刻的状态向量估计值；P′(n)为一步的预测均方误差；P(n-1)为上一时刻均方误差估计值；K(n)为卡尔曼增益；
Figure GDA0004088707120000071
为当前时刻泵噪声的最佳状态向量估计值；P(n)为当前时刻的均方误差估计值；
采用时不变线性模型获得的重构泵噪声信号B(n)表示为：
Figure GDA0004088707120000072
采用时变线性模型获得的重构泵噪声信号B(n)表示为：
Figure GDA0004088707120000073
优选的，采用扩展卡尔曼滤波算法计算泵噪声状态向量估计值
Figure GDA0004088707120000074
其具体方法为：
对非线性模型测量方程的非线性函数h(*)就泵噪声状态向量估计值
Figure GDA0004088707120000075
做一阶泰勒展开，得到：
Figure GDA0004088707120000076
令
Figure GDA0004088707120000078
Figure GDA0004088707120000079
则非线性模型的测量方程改写为：
z(n)＝H′(n)x(n)+γ2(n)+v(n)                (21)
线性化后的观测转移矩阵H′(n)为h(*)的雅克比矩阵，表示为：
Figure GDA00040887071200000710
将H′(n)作为H(n)带入公式(11)至公式(15)进行循环递推即可获得各时刻泵噪声的最佳状态向量估计值
Figure GDA00040887071200000711
进而得到重构泵噪声信号B(n)表示为：
Figure GDA0004088707120000081
优选的，采用无迹卡尔曼滤波算法计算泵噪声状态向量估计值
Figure GDA00040887071200000812
其具体方法为：
由公式(11)和公式(12)得到一步预测状态向量
Figure GDA0004088707120000082
和一步预测均方误差P′(n)；
在
Figure GDA0004088707120000083
附近取4K+1个采样点形成集合σ，σ＝{χi′(n),i∈{0,1,2,3,...4K}}，其中
Figure GDA0004088707120000084
式中，χi′(n)为n时刻的σ采样点值，r＝1,2,…,2K；λ＝α2(2K+κ)-2K；α为σ采样点散布程度因数，α＝10-4～1；κ为缩放参数，保证
Figure GDA0004088707120000085
为半正定矩阵，κ＝0～(3-2K)；
Figure GDA0004088707120000086
为χi′(n)的均值权重，
Figure GDA00040887071200000813
为χi′(n)的方差权重；β为σ采样点值的分布因数；
将χi′(n)代入非线性模型测量方程的非线性函数h(*)，得到观测更新之后的观测向量z′(n)及其均值
Figure GDA0004088707120000087
zi′(n)＝h(n,χi′(n))                    (25)
Figure GDA0004088707120000088
滤波更新为：
Figure GDA0004088707120000089
Figure GDA00040887071200000810
Figure GDA00040887071200000811
Figure GDA0004088707120000091
Figure GDA0004088707120000095
式中，
Figure GDA0004088707120000096
为观测向量z′(n)的误差协方差；
Figure GDA0004088707120000092
为预测状态向量
Figure GDA0004088707120000093
与观测向量z′(n)的误差互协方差；
通过上述公式(27)-(31)循环递推，即可泵噪声状态向量估计值
Figure GDA0004088707120000094
通过公式(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，则总采样点数为
Figure GDA0004088707120000101
各等级幅值的概率为pi＝ni/N；用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类，即e0＝{0,1,…,T}，e1＝{T+1,T+2,…,L-1}，则泥浆脉冲信号E(n)的幅值的均值μ为：
Figure GDA0004088707120000102
e0的均值μ0和e1的均值μ1分别为：
Figure GDA0004088707120000103
Figure GDA0004088707120000104
式中，
Figure GDA0004088707120000105
e0和e1的类间方差
Figure GDA0004088707120000106
为：
Figure GDA0004088707120000111
依照上述流程，使T在[0,L-1]范围内依次取值并计算相应的
Figure GDA0004088707120000112
当
Figure GDA0004088707120000113
最大时所对应的T值即为泥浆脉冲信号E(n)中相应时间段自适应阈值。
与现有技术相比，本发明的优点和积极效果在于：
本发明提供的泥浆脉冲信号处理方法，在充分依据随钻泥浆脉冲信号产生、传输及加噪特性的基础上，将单压力传感器采集的泥浆脉冲信号依次进行泵噪声状态空间模型构建、卡尔曼滤波、去除随机噪声、去除基线漂移和设置自适应阈值处理，获得明显的有效传输脉冲信号。本发明方法简单可靠，能够快速、准确、可靠识别有效脉冲信号，实现井下数据高效传输，去噪效果好，具有很大的实用价值。
附图说明
图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、针对不同的变量参数，根据泵噪声通用模型构建时不变线性模型、时变线性模型和非线性模型三种泵噪声状态空间模型。
具体地，依据泥浆泵工作机理，所述泵噪声通用模型为：
Figure GDA0004088707120000121
式中，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，其中各元素表示为:
Figure GDA0004088707120000122
则时不变线性模型的状态方程表示为：
Figure GDA0004088707120000131
式中，F表示状态转移矩阵，为时不变的常值矩阵；w(n)为过程噪声；
Figure GDA0004088707120000132
x(n+1)为泥浆脉冲传输系统的状态值；
时不变线性模型的测量方程表示为：
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，其中各元素表示为:
Figure GDA0004088707120000133
则时变线性模型的状态方程表示为：
Figure GDA0004088707120000134
式中，I为2阶单位矩阵；
时变线性模型的测量方程表示为：
Figure GDA0004088707120000141
式中，H(n)表示观测转移矩阵，为时变矩阵；
具体地，构建非线性模型的方法为：
设含K阶次谐波的泵噪声的状态向量为x(n)＝[x1(n)x2(n)…x2K-1(n)x2K(n)]T，其中各元素表示为:
Figure GDA0004088707120000142
则非线性模型的状态方程表示为：
x(n+1)＝x(n)+w(n)                     (9)
非线性模型的测量方程表示为：
Figure GDA0004088707120000143
式中，h(*)表示非线性模型测量方程的非线性函数。
S2、采用卡尔曼滤波算法对泵噪声状态空间模型进行计算，获得泵噪声状态向量估计值
Figure GDA0004088707120000144
通过泵噪声状态向量估计值
Figure GDA0004088707120000145
对采集泥浆脉冲信号A(n)中的泵噪声信号进行重构，获得重构泵噪声信号B(n)，则去除泵噪声后的泥浆脉冲信号C(n)＝A(n)-B(n)。卡尔曼滤波算法包括面向线性模型的标准卡尔曼滤波算法，面向非线性模型的扩展卡尔曼算法和无迹卡尔曼滤波算法。泥浆脉冲信号C(n)包括随机噪声W(n)和有效脉冲信号S(n)。
具体地，采用标准卡尔曼滤波算法计算泵噪声状态向量估计值
Figure GDA0004088707120000151
其具体方法为：
将时不变线性模型或时变线性模型的状态转移矩阵F和观测转移矩阵H带入以下公式循环递推，即可获得泵噪声状态向量估计
Figure GDA0004088707120000152
Figure GDA0004088707120000153
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)
Figure GDA0004088707120000154
P(n)＝[I-K(n)H(n)]P′(n)                  (15)
式中，Q(n-1)为过程噪声w(n)的协方差矩阵；R(n)为测量噪声v(n)的协方差矩阵；
Figure GDA0004088707120000155
为一步的预测状态向量，即状态向量的一步预测值，也就是状态向量的先验估计；
Figure GDA0004088707120000156
为上一时刻的状态向量估计值；P′(n)为一步的预测均方误差；P(n-1)为上一时刻均方误差估计值，即均方误差的一步预测值，也就是均方误差的先验估计；K(n)为卡尔曼增益；
Figure GDA0004088707120000157
为当前时刻泵噪声的最佳状态向量估计值；P(n)为当前时刻的均方误差估计值。
需要说明的是，过程噪声w(n)的协方差矩阵Q(n)、测量噪声v(n)的协方差矩阵R(n)，以及初始状态向量
Figure GDA0004088707120000158
和初始均方误差P(0)可根据实际噪声情况灵活调整。
采用时不变线性模型获得的重构泵噪声信号B(n)表示为：
Figure GDA0004088707120000159
采用时变线性模型获得的重构泵噪声信号B(n)表示为：
Figure GDA00040887071200001510
。
具体地，采用扩展卡尔曼滤波算法计算泵噪声状态向量估计值
Figure GDA0004088707120000162
其具体方法为：
对非线性模型测量方程的非线性函数h(*)就泵噪声状态向量估计值
Figure GDA0004088707120000163
做一阶泰勒展开，得到：
Figure GDA00040887071200001612
令
Figure GDA0004088707120000165
Figure GDA0004088707120000166
则非线性模型的测量方程改写为：
z(n)＝H′(n)x(n)+γ2(n)+v(n)                (21)
线性化后的观测转移矩阵H′(n)为h(*)的雅克比矩阵，表示为：
Figure GDA0004088707120000167
将H′(n)作为H(n)带入公式(11)至公式(15)进行循环递推即可获得各时刻泵噪声的最佳状态向量估计值
Figure GDA0004088707120000168
进而得到重构泵噪声信号B(n)表示为：
Figure GDA0004088707120000169
具体地，采用无迹卡尔曼滤波算法计算泵噪声状态向量估计值
Figure GDA00040887071200001610
其具体方法为：
由公式(11)和公式(12)得到一步预测状态向量
Figure GDA00040887071200001611
和一步预测均方误差P′(n)；
在
Figure GDA0004088707120000171
附近取4K+1个采样点形成集合σ，σ＝{χi′(n)，i∈{0，1,2,3,...4K}}，其中
Figure GDA0004088707120000172
式中，χi′(n)为n时刻的σ采样点值，r＝1,2,…,2K；λ＝α2(2K+κ)-2K；α为σ采样点散布程度因数，α＝10-4～1；κ为缩放参数，保证
Figure GDA0004088707120000173
为半正定矩阵，κ＝0～(3-2K)；
Figure GDA0004088707120000174
为χi′(n)的均值权重，
Figure GDA0004088707120000175
为χi′(n)的方差权重；β为σ采样点值的分布因数。需要说明的是，对于高斯分布的状态向量，β＝2最佳。
将χi′(n)代入非线性模型测量方程的非线性函数h(*)，得到观测更新之后的观测向量z′(n)及其均值
Figure GDA0004088707120000176
zi′(n)＝h(n,χi′(n))                    (25)
Figure GDA0004088707120000177
滤波更新为：
Figure GDA0004088707120000178
Figure GDA0004088707120000179
Figure GDA00040887071200001710
Figure GDA00040887071200001711
Figure GDA00040887071200001712
式中，
Figure GDA00040887071200001713
为观测向量z′(n)的误差协方差；
Figure GDA00040887071200001714
为预测状态向量
Figure GDA00040887071200001715
与观测向量z′(n)的误差互协方差；
通过上述公式(27)-(31)循环递推，即可泵噪声状态向量估计值
Figure GDA0004088707120000181
通过公式(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，则总采样点数为
Figure GDA0004088707120000201
各等级幅值的概率为pi＝ni/N；用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类，即e0＝{0,1,…,T}，e1＝{T+1,T+2,…,L-1}，则泥浆脉冲信号E(n)的幅值的均值μ为：
Figure GDA0004088707120000202
e0的均值μ0和e1的均值μ1分别为：
Figure GDA0004088707120000203
Figure GDA0004088707120000204
式中，
Figure GDA0004088707120000205
e0和e1的类间方差
Figure GDA0004088707120000206
为：
Figure GDA0004088707120000207
依照上述流程，使T在[0,L-1]范围内依次取值并计算相应的
Figure GDA0004088707120000208
当
Figure GDA0004088707120000209
最大时所对应的T值即为泥浆脉冲信号E(n)中相应时间段自适应阈值。
本发明上述方法在充分依据随钻泥浆脉冲信号产生、传输及加噪特性的基础上，将单压力传感器采集的泥浆压力信号依次进行泵噪声状态空间模型构建、卡尔曼滤波、抑制随机噪声、消除基线漂移和设置自适应阈值等处理，可以获得明显的有效传输脉冲信号。该方法简便可靠，具有很大的实用价值。
以下结合具体地实施例对本发明上述方法作出进一步说明。
参见图2，图2所示为地面立管处连续采集的某段100s含泥浆脉冲信号的原始压力信号A(n)，首先对含泥浆脉冲信号的原始压力信号A(n)进行FFT变换，获得原始泥浆压力信号频谱FA(n)(如图3)。根据该频谱特征采用本发明所述泥浆脉冲从处理方法对上述原始压力信号A(n)进行处理，得到波形特征明显的泥浆脉冲信号，用于实现解码及后续数据处理。其具体步骤为：
S1、依据泥浆泵工作机理，泵噪声通用模型表示为：
Figure GDA0004088707120000211
式中，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，其中各元素表示为:
Figure GDA0004088707120000221
则泵噪声状态空间模型的状态方程表示为：
Figure GDA0004088707120000222
式中，F表示状态转移矩阵，为时不变的常值矩阵；w(n)为过程噪声；
Figure GDA0004088707120000223
x(n+1)为泥浆脉冲传输系统的状态值；
泵噪声状态空间模型的测量方程表示为：
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；泵噪声状态空间模型为线性时不变模型。
S2、采用标准卡尔曼滤波算法对泵噪声状态空间模型进行计算，获得泵噪声状态向量估计值
Figure GDA0004088707120000224
其具体方法为：
将泵噪声状态空间模型的状态转移矩阵F和观测转移矩阵H带入以下公式循环递推，即可获得泵噪声状态向量估计
Figure GDA0004088707120000225
Figure GDA0004088707120000226
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)
Figure GDA0004088707120000231
P(n)＝[I-K(n)H(n)]P′(n)                  (15)
式中，Q(n-1)为过程噪声w(n)的协方差矩阵，R(n)为测量噪声v(n)的协方差矩阵，
Figure GDA0004088707120000232
为一步的预测状态向量；
Figure GDA0004088707120000233
为上一时刻的状态向量估计值；P′(n)为一步的预测均方误差；P(n-1)为上一时刻均方误差估计值；K(n)为卡尔曼增益；
Figure GDA0004088707120000234
为当前时刻泵噪声的最佳状态向量估计值；P(n)为当前时刻的均方误差估计值。
本实施例中，过程噪声w(n)和测量噪声v(n)设为均值为0、方差分别是0.026和1的高斯白噪声；初始状态向量
Figure GDA0004088707120000235
为2K×1维的零向量，初始均方误差P(0)为2K×2K维的单位矩阵。
通过迭代获得的泵噪声状态向量估计
Figure GDA0004088707120000236
对泵噪声信号进行重构，获得重构泵噪声信号B(n)(如图4)表示为：
Figure GDA0004088707120000237
从含泥浆脉冲信号的原始压力信号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，则总采样点数为
Figure GDA0004088707120000241
各等级幅值的概率为pi＝ni/N；用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类，即e0＝{0,1,…,T}，e1＝{T+1,T+2,…,L-1}，则泥浆脉冲信号E(n)的幅值的均值μ为：
Figure GDA0004088707120000251
e0的均值μ0为：
Figure GDA0004088707120000252
式中，
Figure GDA0004088707120000253
e1的均值μ1为：
Figure GDA0004088707120000254
式中，
Figure GDA0004088707120000255
e0和e1的类间方差
Figure GDA0004088707120000256
为：
Figure GDA0004088707120000257
依照上述流程，使T在[0,L-1]范围内依次取值并计算相应的
Figure GDA0004088707120000258
当
Figure GDA0004088707120000259
最大时所对应的T值即为泥浆脉冲信号E(n)中相应时间段自适应阈值。
将泥浆脉冲信号E(n)中小于相应时间段自适应阈值T的采样点值置零，得到泥浆脉冲信号信号F(n)，如图8所示。
上述实施例用来解释本发明，而不是对本发明进行限制，在本发明的精神和权利要求的保护范围内，对本发明做出的任何修改和改变，都落入本发明的保护范围。

Claims (8)

1.一种泥浆脉冲信号处理方法，其特征在于，其具体步骤为：
根据泵噪声通用模型构建泵噪声状态空间模型；依据泥浆泵工作机理，所述泵噪声通用模型为：
Figure FDA0004088707100000011
式中，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，其中各元素表示为:
Figure FDA0004088707100000012
则时不变线性模型的状态方程表示为：
Figure FDA0004088707100000013
式中，F表示状态转移矩阵，为时不变的常值矩阵；w(n)为过程噪声；
Figure FDA0004088707100000021
x(n+1)为泥浆脉冲传输系统的状态值；时不变线性模型的测量方程表示为：
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，其中各元素表示为:
Figure FDA0004088707100000022
则时变线性模型的状态方程表示为：
Figure FDA0004088707100000023
式中，I为2阶单位矩阵；
时变线性模型的测量方程表示为：
Figure FDA0004088707100000024
式中，H(n)表示观测转移矩阵，为时变矩阵；
构建非线性模型的方法为：
设含K阶次谐波的泵噪声的状态向量为x(n)＝[x1(n) x2(n) … x2K-1(n) x2K(n)]T，其中各元素表示为:
Figure FDA0004088707100000031
则非线性模型的状态方程表示为：
x(n+1)＝x(n)+w(n)                     (9)
非线性模型的测量方程表示为：
Figure FDA0004088707100000032
式中，h(*)表示非线性模型测量方程的非线性函数；
采用卡尔曼滤波算法对泵噪声状态空间模型进行计算，获得泵噪声状态向量估计值
Figure FDA0004088707100000033
通过泵噪声状态向量估计值
Figure FDA0004088707100000034
对采集泥浆脉冲信号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所述的泥浆脉冲信号处理方法，其特征在于，采用标准卡尔曼滤波算法计算泵噪声状态向量估计值
Figure FDA0004088707100000041
其具体方法为：将时不变线性模型或时变线性模型的状态转移矩阵F和观测转移矩阵H带入以下公式循环递推，即可获得泵噪声状态向量估计
Figure FDA0004088707100000042
Figure FDA0004088707100000043
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)
Figure FDA0004088707100000044
P(n)＝[I-K(n)H(n)]P′(n)                  (15)
式中，Q(n-1)为过程噪声w(n)的协方差矩阵；R(n)为测量噪声v(n)的协方差矩阵；
Figure FDA0004088707100000045
为一步的预测状态向量；
Figure FDA0004088707100000046
为上一时刻的状态向量估计值；P′(n)为一步的预测均方误差；P(n-1)为上一时刻均方误差估计值；K(n)为卡尔曼增益；
Figure FDA0004088707100000047
为当前时刻泵噪声的最佳状态向量估计值；P(n)为当前时刻的均方误差估计值；
采用时不变线性模型获得的重构泵噪声信号B(n)表示为：
Figure FDA0004088707100000048
采用时变线性模型获得的重构泵噪声信号B(n)表示为：
Figure FDA0004088707100000049
3.如权利要求1所述的泥浆脉冲信号处理方法，其特征在于，采用扩展卡尔曼滤波算法计算泵噪声状态向量估计值
Figure FDA00040887071000000410
其具体方法为：对非线性模型测量方程的非线性函数h(*)就泵噪声状态向量估计值
Figure FDA00040887071000000411
做一阶泰勒展开，得到：
Figure FDA0004088707100000051
令
Figure FDA0004088707100000052
Figure FDA0004088707100000053
则非线性模型的测量方程改写为：
z(n)＝H′(n)x(n)+γ2(n)+v(n)                (21)
线性化后的观测转移矩阵H′(n)为h(*)的雅克比矩阵，表示为：
Figure FDA0004088707100000054
将H′(n)作为H(n)带入公式(11)至公式(15)进行循环递推即可获得各时刻泵噪声的最佳状态向量估计值
Figure FDA0004088707100000055
进而得到重构泵噪声信号B(n)表示为：
Figure FDA0004088707100000056
4.如权利要求1所述的泥浆脉冲信号处理方法，其特征在于，采用无迹卡尔曼滤波算法计算泵噪声状态向量估计值
Figure FDA0004088707100000057
其具体方法为：由公式(11)和公式(12)得到一步预测状态向量
Figure FDA0004088707100000058
和一步预测均方误差P′(n)；
在
Figure FDA0004088707100000059
附近取4K+1个采样点形成集合σ，σ＝{χi′(n),i∈{0,1,2,3,...4K}}，其中
Figure FDA00040887071000000510
式中，χi′(n)为n时刻的σ采样点值，r＝1,2,…,2K；λ＝α2(2K+κ)-2K；
α为σ采样点散布程度因数，α＝10-4～1；κ为缩放参数，保证
Figure FDA0004088707100000061
为半正定矩阵，κ＝0～(3-2K)；
Figure FDA0004088707100000062
为χi′(n)的均值权重，Ψi c为χi′(n)的方差权重；β为σ采样点值的分布因数；
将χi′(n)代入非线性模型测量方程的非线性函数h(*)，得到观测更新之后的观测向量z′(n)及其均值
Figure FDA0004088707100000063
zi′(n)＝h(n,χi′(n))                    (25)
Figure FDA0004088707100000064
滤波更新为：
Figure FDA0004088707100000065
Figure FDA0004088707100000066
Figure FDA0004088707100000067
Figure FDA0004088707100000068
Figure FDA0004088707100000069
式中，
Figure FDA00040887071000000610
为观测向量z′(n)的误差协方差；
Figure FDA00040887071000000611
为预测状态向量
Figure FDA00040887071000000612
与观测向量z′(n)的误差互协方差；
通过上述公式(27)-(31)循环递推，即可泵噪声状态向量估计值
Figure FDA00040887071000000613
通过公式(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，则总采样点数为
Figure FDA0004088707100000081
各等级幅值的概率为pi＝ni/N；用整数T将泥浆脉冲信号E(n)中的采样点按幅值等级高低划分为e0和e1两类，即e0＝{0,1,…,T}，e1＝{T+1,T+2,…,L-1}，则泥浆脉冲信号E(n)的幅值的均值μ为：
Figure FDA0004088707100000082
e0的均值μ0和e1的均值μ1分别为：
Figure FDA0004088707100000083
Figure FDA0004088707100000084
式中，
Figure FDA0004088707100000085
e0和e1的类间方差
Figure FDA0004088707100000086
为：
Figure FDA0004088707100000087
依照上述流程，使T在[0,L-1]范围内依次取值并计算相应的
Figure FDA0004088707100000088
当
Figure FDA0004088707100000089
最大时所对应的T值即为泥浆脉冲信号E(n)中相应时间段自适应阈值。
CN202010381133.3A 2020-05-08 2020-05-08 泥浆脉冲信号处理方法 Active CN111535802B (zh)

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 (6)

* Cited by examiner, † Cited by third party
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
CN112418174A (zh) * 2020-12-08 2021-02-26 中国石油天然气集团有限公司 随钻泥浆随机噪声的去除方法
CN112554875A (zh) * 2020-12-08 2021-03-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 中石化经纬有限公司 一种随钻测量连续波信号去噪方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
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 浙江大学 应用经验模态分解和粒子群优化算法的泵噪消除方法

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
US9007232B2 (en) Mud pulse telemetry noise reduction method
CN107083957B (zh) 钻井液随钻信号的泵冲干扰消除方法及系统
CN105738948B (zh) 一种基于小波变换的微地震数据降噪方法
US20150168573A1 (en) Geologic quality factor inversion method
CN104849756B (zh) 一种提高地震数据分辨率增强有效弱信号能量的方法
CN103645507B (zh) 地震记录的处理方法
CN109031422A (zh) 一种基于CEEMDAN与Savitzky-Golay滤波的地震信号噪声抑制方法
CN112835103B (zh) 自适应去鬼波与宽频准零相位反褶积联合处理方法及系统
CN107179550B (zh) 一种数据驱动的地震信号零相位反褶积方法
CN109143331B (zh) 地震子波提取方法
US8942330B2 (en) Interference reduction method for downhole telemetry systems
CN108508489B (zh) 一种基于波形微变化匹配的地震反演方法
CN105041303B (zh) 钻井液随钻数据传输系统的泵冲干扰信号消除方法
Chen et al. Study on model-based pump noise suppression method of mud pulse signal
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
WO2021020985A1 (en) A method and system for monitoring a wellbore object using a reflected pressure signal
CN113341463B (zh) 一种叠前地震资料非平稳盲反褶积方法及相关组件
CN110749923A (zh) 一种基于范数方程提高分辨率的反褶积方法
CN112554875A (zh) 随钻泥浆正脉冲信号的处理方法
DONG et al. Fast implementation technique of adaptive Kalman filtering deconvolution via dyadic wavelet transform
Li et al. Random-noise attenuation by an amplitude-preserved time-frequency peak based on empirical wavelet transform predictive filtering

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