CN105046025B - 一种核磁共振多相流测量中各相分离的方法 - Google Patents

一种核磁共振多相流测量中各相分离的方法 Download PDF

Info

Publication number
CN105046025B
CN105046025B CN201510545805.9A CN201510545805A CN105046025B CN 105046025 B CN105046025 B CN 105046025B CN 201510545805 A CN201510545805 A CN 201510545805A CN 105046025 B CN105046025 B CN 105046025B
Authority
CN
China
Prior art keywords
mrow
signal
imf
msub
oil
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.)
Expired - Fee Related
Application number
CN201510545805.9A
Other languages
English (en)
Other versions
CN105046025A (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.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
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 Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to CN201510545805.9A priority Critical patent/CN105046025B/zh
Publication of CN105046025A publication Critical patent/CN105046025A/zh
Application granted granted Critical
Publication of CN105046025B publication Critical patent/CN105046025B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30Assessment of water resources

Abstract

本发明公开了一种核磁共振测量多相流分相含率的信号分离方法。多相流核磁共振测量信号是一维数组,如式(1)所示,该方法能够分离出与油、气、水分相含率相关的核磁共振衰减系数A、A、A。首先将用核磁共振方法测量到的包含油气水三相物质的自由衰减信号(FID)进行经验模态分解(EMD),得到若干个本征模态函数(imf(t));然后imf(t)从低阶往高阶依次叠加,分别计算与测量信号的相关系数,找到使相关系数产生跳变的imf(t),将该imf(t)至最后一个imf(t)累加获得重构信号,达到滤波的目的;最后将滤波后的信号进行分段拟合,得到各相参数的估计值。其显著效果是:本发明滤波效果相较于小波滤波法、最小均方根滤波法,滤波效果好,自适应性强,经过分段拟合后的参数估计误差小。

Description

一种核磁共振多相流测量中各相分离的方法
技术领域
本发明属于多相流测量领域,涉及到一种核磁共振多相流测量中各相分离的方法。
背景技术
多相流指的是气体、固体、液体三相或其中两相混合的复杂流体,在能源、环保、石油等工程领域大量存在。在石油工业中,从油井开采出来的石油是多组分的油、气、水混合多相流,是最复杂的多相流之一。油气水多相流检测是对油、气、水三种物质的含量和速度进行在线精确计量。目前多相流检测方法主要有射线法、电学层析成像法、核磁共振法等。核磁共振技术(Nuclear Magnetic Resonance,NMR)由于其穿透性与高精度等特点,在多相流检测中逐渐显示出优势。在利用核磁共振法进行多相流检测时,检测到的核磁共振自由感应衰减信号(Free Induction Decay,FID)包含油、气、水含量等重要信息。有效去除FID信号中的噪声是提高多相流检测精度的根本措施。
经验模态分解(Empirical Mode Decomposition,EMD)方法是HUANG等人提出的一种信号分析方法。它依据数据自身的时间尺度特征进行信号分解,自适应强,能处理非线性、非平稳数据。
基于EMD的滤波算法的关键是找到包含原信号较多的分量,剔除含噪声信号较多的分量,然后进行重构。对此国内外学者进行大量研究。Boudraa等人在《EMD-Based SignalFiltering》提出基于连续均方差的EMD滤波算法,即找到imf能量的全局极小值位置作为原信号占主导的分量与噪声占主导分量的临界点,但当某些原信号占主导的分量的能量低于噪声占主导分量时,分选准则存在偏差。林丽等人在《基于相关系数的EMD改进算法》中提出了基于相关系数的EMD滤波算法,设置相关系数序列中最大值的十分之一为阀值,认为相关系数小于该值的imf分量是伪分量,应该剔除。但是当噪声信号集中在某个imf分量,该imf分量的相关系数值也会很大,会被错误判断为真分量,导致没有被剔除。贾瑞生等人在《基于经验模态分解及独立成分分析的微震信号降噪方法》提出基于相关系数最小值的分选方法,但是容易陷入局部最优解,导致分选偏差。
发明内容
本发明的目的是针对现有算法的不足,提供一种核磁共振多相流测量中各相分离的方法。
本发明的目的是通过以下技术方案实现的。
一种核磁共振测量油气水多相流分相含率的信号分离方法包括以下步骤:
1)用核磁共振方法测量到的包含油气水三相含率信息的自由衰减信号,对自由衰减信号进行经验模态分解,得到m个本征模态函数imf(t);
2)将本征模态函数imf(t)从低阶往高阶依次叠加,分别计算与测量信号的相关系数,找到使相关系数产生跳变的本征模态函数imf(t),将该本征模态函数imf(t)至最后一个本征模态函数imf(t)累加获得重构信号,达到滤波的目的;
3)将滤波后的信号进行分段拟合,得到油气水各相参数的估计值。
所述的步骤1)包括:
2.1)核磁共振方法测量的自由衰减信号用公式(1)描述:
y(t)=A·exp(-t/T)+A·exp(-t/T)+A·exp(-t/T)+n(t) (1)
其中A、A、A分别为油气水各相自由衰减信号的幅值,代表油气水各相的含量,油气水多相流分相含率的信号分离方法的目的是计算出A、A、A的估计值;T、T、T分别为油气水各相自由衰减信号的横向弛豫时间,油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s,T、T、T在核磁共振测量中用以区分各相物质;n(t)为测量中的噪声信号,是需要剔除的信号;
2.2)将测量的自由衰减信号y(t)按如下步骤进行经验模态分解:
2.2.1)分别找到测量的自由衰减信号y(t)的所有极大值点和极小值点;然后通过三次样条函数,将所有极大值点拟合出包络线e+(t),将所有极小值点拟合出包络线e-(t),按公式计算均值包络线m1(t);
2.2.2)将测量信号y(t)减去m1(t)得到一个去掉低频信号的新信号
2.2.3)重复步骤2.2.1)至步骤2.2.2),假设经过k次之后满足对于整个信号,极值点的个数和过零点的个数相差不超过1;并且在任意时刻,局部极大值和极小值的包络线的平均值是零;则测量信号y(t)的一阶imf分量为;
2.2.4)将测量信号y(t)减去imf1(t),获得一个去掉高频分量的残余分量res1(t);
res1(t)=y(t)-imf1(t) (5)
2.2.5)再对残余分量res1(t)重复步骤2.2.1)至步骤2.2.4),得到第二个分量imf2(t),如此进行多次分解一直到残余分量resn(t)是常量或单调函数时,经验模态分解过程结束,得到:
2.3)将残余分量resn(t)当做最后一个imf分量,imfn+1(t)=resn(t)得到:
其中m=n+1。
所述的步骤2)包括:
3.1)在用核磁共振技术测量多相流时,对多数包含噪声的测量信号,其主要信号的能量集中在低频段,大部分噪声分布在高频段;因此存在某个imfk(t)分量,使得对于该分量之后的imf(t)中信号占主要部分,而其前k个imf(t)中噪声占主要部分;
3.2)imf(t)的信号重构函数C(j)
噪声和信号都与测量信号存在相关性,每个imf分量与测量信号y(t)也存在一定相关性,随着分量从低阶到高阶依次相加重构,重构后的信号与测量信号的相关性依次增大;对于信噪比大于0dB的测量信号,信号与测量信号的相关性比噪声与测量信号的相关性更大;因此在加入imfk+1(t)分量之后,重构信号与测量信号的相关系数会有一个跳变;通过步骤1)将测量信号y(t)分解为m个imf(t)分量,分别将分量从低阶到高阶依次相加重构得到m个函数,C(j),j=1,2,…,m;
3.3)C(j)和y(t)的相关系数corr(j)计算
分别计算C(j)和输入信号的相关系数corr(j),corr(j)的绝对值取值在0到1之间,越接近1表示相关性越强;计算相关系数的公式为:
其中分别为数组C(j)、y(t)的平均值,N为采样点数;
3.4)k值计算
过多次仿真实验,确定相关系数产生跳变的阀值为5,即当加入imfk+1(t)分量后的相关系数绝对值的增量是加入imfk(t)分量后的相关系数绝对值增量的5倍,k即为所求的值;
3.5)滤波后的信号y2(t)
用核磁共振技术测量多相流的自由衰减信号,如式(1)所示,其有效信号的能量集中在低频段,大部分噪声分布在高频段;据此,存在某个imfk(t)分量,使得对于该分量之后的imfk+1(t)、、、imfm(t)中信号占主要部分,imf1(t)、、、imfk(t)中噪声占主要部分,剔除这部分imf(t),就达到了滤波的效果;即信号y(t)滤波后的信号为y2(t):
所述的步骤3)包括:
4.1)经过滤波后的自由衰减信号y2(t)是一个数组,需要通过指数函数拟合的方法得到表达式,完成油气水三相物质的自由衰减信号分离,采用最小二乘法和模拟退火优化算法对滤波后的自由衰减信号进行拟合,步骤2)获得滤波后的自由衰减信号为:
y2(t)=A·exp(-t/T)+A·exp(-t/T)+A·exp(-t/T) (12)
将y2(t)分为t<0.3s和t>0.3s两段;
4.2)对于指数衰减函数A·exp(-t/T),当t>5T时,幅值衰减为0.0067A,近似为0;油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s因此当t>0.3s,气体的自由衰减信号认为衰减至0,则原三相自由衰减信号测量信号则变成油水两相,简化成两相拟合问题,结果得到:
其中分别为A、A、T、T的参数估计值;
4.3)由步骤4.2)中已经求出油、水两相参数的估计值,因此截取测量信号中t<0.3s的数据,直接拟合求出气体单项参数的估计值,得到:
其中分别为A、T的参数估计值;结果中代表着油气水各相含量的估计值。
本发明对比已有的技术具有以下创新点:
1.利用相关系数跳变找到临界imf分量
2.利用拟合算法对多相流各相FID信号进行分离
本发明对比已有技术具有以下显著优点:
1.基于EMD滤波信噪比高(SNR),均方根误差(MSE)小
2.经过分段拟合后参数估计误差小
附图说明
图1是各相分离算法的流程图;
图2是测量信号y(t);
图3是滤波效果图,其中y2(t)是滤波后的信号,是未加噪声的信号;
图4是三种滤波算法效果对比。
具体实施方式
一种核磁共振测量油气水多相流分相含率的信号分离方法包括以下步骤:
1)用核磁共振方法测量到的包含油气水三相含率信息的自由衰减信号,对自由衰减信号进行经验模态分解,得到m个本征模态函数imf(t);
2)将本征模态函数imf(t)从低阶往高阶依次叠加,分别计算与测量信号的相关系数,找到使相关系数产生跳变的本征模态函数imf(t),将该本征模态函数imf(t)至最后一个本征模态函数imf(t)累加获得重构信号,达到滤波的目的;
3)将滤波后的信号进行分段拟合,得到油气水各相参数的估计值。
所述的步骤1)包括:
2.1)核磁共振方法测量的自由衰减信号用公式(1)描述:
y(t)=A·exp(-t/T)+A·exp(-t/T)+A·exp(-t/T)+n(t) (1)
其中A、A、A分别为油气水各相自由衰减信号的幅值,代表油气水各相的含量,油气水多相流分相含率的信号分离方法的目的是计算出A、A、A的估计值;T、T、T分别为油气水各相自由衰减信号的横向弛豫时间,油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s,T、T、T在核磁共振测量中用以区分各相物质;n(t)为测量中的噪声信号,是需要剔除的信号;
2.2)将测量的自由衰减信号y(t)按如下步骤进行经验模态分解:
2.2.1)分别找到测量的自由衰减信号y(t)的所有极大值点和极小值点;然后通过三次样条函数,将所有极大值点拟合出包络线e+(t),将所有极小值点拟合出包络线e-(t),按公式计算均值包络线m1(t);
2.2.2)将测量信号y(t)减去m1(t)得到一个去掉低频信号的新信号
2.2.3)重复步骤2.2.1)至步骤2.2.2),假设经过k次之后满足对于整个信号,极值点的个数和过零点的个数相差不超过1;并且在任意时刻,局部极大值和极小值的包络线的平均值是零;则测量信号y(t)的一阶imf分量为;
2.2.4)将测量信号y(t)减去imf1(t),获得一个去掉高频分量的残余分量res1(t);
res1(t)=y(t)-imf1(t) (5)
2.2.5)再对残余分量res1(t)重复步骤2.2.1)至步骤2.2.4),得到第二个分量imf2(t),如此进行多次分解一直到残余分量resn(t)是常量或单调函数时,经验模态分解过程结束,得到:
2.3)将残余分量resn(t)当做最后一个imf分量,imfn+1(t)=resn(t)得到:
其中m=n+1。
所述的步骤2)包括:
3.1)在用核磁共振技术测量多相流时,对多数包含噪声的测量信号,其主要信号的能量集中在低频段,大部分噪声分布在高频段;因此存在某个imfk(t)分量,使得对于该分量之后的imf(t)中信号占主要部分,而其前k个imf(t)中噪声占主要部分;
3.2)imf(t)的信号重构函数C(j)
噪声和信号都与测量信号存在相关性,每个imf分量与测量信号y(t)也存在一定相关性,随着分量从低阶到高阶依次相加重构,重构后的信号与测量信号的相关性依次增大;对于信噪比大于0dB的测量信号,信号与测量信号的相关性比噪声与测量信号的相关性更大;因此在加入imfk+1(t)分量之后,重构信号与测量信号的相关系数会有一个跳变;通过步骤1)将测量信号y(t)分解为m个imf(t)分量,分别将分量从低阶到高阶依次相加重构得到m个函数,C(j),j=1,2,…,m;
3.3)C(j)和y(t)的相关系数corr(j)计算
分别计算C(j)和输入信号的相关系数corr(j),corr(j)的绝对值取值在0到1之间,越接近1表示相关性越强;计算相关系数的公式为:
其中分别为数组C(j)、y(t)的平均值,N为采样点数;
3.4)k值计算
过多次仿真实验,确定相关系数产生跳变的阀值为5,即当加入imfk+1(t)分量后的相关系数绝对值的增量是加入imfk(t)分量后的相关系数绝对值增量的5倍,k即为所求的值;
3.5)滤波后的信号y2(t)
用核磁共振技术测量多相流的自由衰减信号,如式(1)所示,其有效信号的能量集中在低频段,大部分噪声分布在高频段;据此,存在某个imfk(t)分量,使得对于该分量之后的imfk+1(t)、、、imfm(t)中信号占主要部分,imf1(t)、、、imfk(t)中噪声占主要部分,剔除这部分imf(t),就达到了滤波的效果;即信号y(t)滤波后的信号为y2(t):
所述的步骤3)包括:
4.1)经过滤波后的自由衰减信号y2(t)是一个数组,需要通过指数函数拟合的方法得到表达式,完成油气水三相物质的自由衰减信号分离,采用最小二乘法和模拟退火优化算法对滤波后的自由衰减信号进行拟合,步骤2)获得滤波后的自由衰减信号为:
y2(t)=A·exp(-t/T)+A·exp(-t/T)+A·exp(-t/T) (12)
将y2(t)分为t<0.3s和t>0.3s两段;
4.2)对于指数衰减函数A·exp(-t/T),当t>5T时,幅值衰减为0.0067A,近似为0;油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s因此当t>0.3s,气体的自由衰减信号认为衰减至0,则原三相自由衰减信号测量信号则变成油水两相,简化成两相拟合问题,结果得到:
其中分别为A、A、T、T的参数估计值;
4.3)由步骤4.2)中已经求出油、水两相参数的估计值,因此截取测量信号中t<0.3s的数据,直接拟合求出气体单项参数的估计值,得到:
其中分别为A、T的参数估计值;结果中代表着油气水各相含量的估计值。
下面结合附图和具体实施例对本发明作进一步的详细说明。
实施例:
如图1所示,本实施例采用MATLAB仿真来验证算法的效果,设采样频率为3288Hz,采样时间长度为1s。油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s。
假设FID信号y(t)为:
y(t)=10*exp(-t/0.8)+2*exp(-t/0.04)+n(t)+5*exp(-t/0.4)+n(t) (15)
其中n(t)是信噪比为16dB的白噪声,y(t)如图2所示。y(t)按如下步骤进行EMD分解:
分别找到测量信号y(t)的所有极大值点和极小值点;然后通过三次样条函数,将所有极大值点拟合出包络线e+(t),将所有极小值点拟合出包络线e-(t),按公式计算均值包络线m1(t);
将测量信号y(t)减去m1(t)得到一个去掉低频信号的新信号
重复上述步骤,假设经过k次之后满足对于整个信号,极值点的个数和过零点的个数相差不超过1;并且在任意时刻,局部极大值和极小值的包络线的平均值是零;则测量信号y(t)的一阶imf分量为;
将测量信号y(t)减去imf1(t),获得一个去掉高频分量的残余分量res1(t);
res1(t)=y(t)-imf1(t) (19)
再对res1(t)重复上述步骤,得到第二个分量imf2(t),如此进行多次分解一直到残余分量resn(t)是常量或单调函数时,EMD分解过程结束,得到:
为了方便处理将残余分量res7(t)当做最后一个imf(t)分量,imf8(t)=res7(t)得到:
在用核磁共振技术测量多相流时,对多数包含噪声的测量信号,其主要信号的能量集中在低频段,大部分噪声分布在高频段;因此存在某个imfk(t)分量,使得对于该分量之后的imf(t)中信号占主要部分,而其前k个imf(t)中噪声占主要部分。
噪声和信号都与测量信号存在相关性,每个imf分量与测量信号y(t)也存在一定相关性,随着分量从低阶到高阶依次相加重构,重构后的信号与测量信号的相关性依次增大;对于信噪比大于0dB的测量信号,信号与测量信号的相关性比噪声与测量信号的相关性更大;因此在加入imfk+1(t)分量之后,重构信号与测量信号的相关系数会有一个跳变;通过EMD方法,将测量信号y(t)分解为8个imf(t)分量,分别将分量从低阶到高阶依次相加重构得到8个函数,C(j),j=1,2,…,8。
分别计算C(j)和输入信号的相关系数corr(j),corr(j)的绝对值取值在0到1之间,越接近1表示相关性越强;计算相关系数的公式为:
其中分别为数组C(j)、y(t)的平均值,N为采样点数;
计算得到的相关系数如表1所示:
表1皮尔森相关系数
过多次仿真实验,确定相关系数产生跳变的阀值为5,即当加入imfk+1(t)分量后的相关系数绝对值的增量是加入imfk(t)分量后的相关系数绝对值增量的5倍,k即为所求的值,计算得到k=7。
用核磁共振技术测量多相流的FID信号,其有效信号的能量集中在低频段,大部分噪声分布在高频段;据此,存在某个imfk(t)分量,使得对于该分量之后的imfk+1(t)、、、imfm(t)中信号占主要部分,imf1(t)、、、imfk(t)中噪声占主要部分,剔除这部分imf(t),就达到了滤波的效果;即信号y(t)滤波后的信号为y2(t),如图2所示。
y2(t)=imf8(t)=res7(t) (25)
为了验证本文得到k=7的正确性,将经过EMD分解后的7个本征模态函数和一个残余分量从后往前依次累加,重构成8个信号。并采用信噪比(Signal To NoiseRatio,SNR)和均方根误差(Mean Square Error,MSE)两个指标来进行衡量和比较。其中,信噪比越大说明滤波效果越好,而均方根误差越小滤波效果越好,如表2所示。
表2重构信号对比
从表2可以看出由res7(t)分量重构获得的信号的SNR最高,并且MSE最小,信号滤波效果最佳。
表3符号对应的表达式
为了验证该滤波算法的优越性,现对测量信号分别利用最小均方(Least MeanSquare,LMS)滤波法、小波滤波法进行滤波处理,然后对比三种方法的滤波效果,这里也采用SNR和MSE作为评价指标。
表4三种滤波算法的效果对比
从表4可知,相比于小波滤波法和LMS法滤波,采用EMD滤波法得到的SNR最高,而MSE最小,即滤波效果最好。
经过滤波后的FID信号y2(t)是一个数组,需要通过指数函数拟合的方法得到表达式,完成油、气、水三相物质的FID信号分离,采用最小二乘法和模拟退火优化算法对滤波后的FID信号进行拟合,获得滤波后的FID信号为:
y2(t)=A·exp(-t/T)+A·exp(-t/T)+A·exp(-t/T) (26)
由于三相FID含有6个参数,当直接进行拟合时,拟合结果如表5所示。
表5拟合结果
从上述结果可以看出三相FID信号拟合的效果不理想,误差较大。因此本发明采用分段拟合的方法来改进拟合效果。
对于指数衰减函数A·exp(-t/T),当t>5T时,幅值衰减为0.0067A,可以近似为0;油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s因此当t>0.3s,气体的FID信号可以认为衰减至0,则原三相FID测量信号则变成油水两相,简化成两相拟合问题,结果得到:
其中分别为A、A、T、T的参数估计值,如表6所示;
由上述步骤求出油、水两相参数的估计值,截取测量信号中t<0.3s的数据,直接拟合求出气体单项参数的估计值,得到:
其中分别为A、T的参数估计值,如表6所示。
表6分段拟合的结果
由表5和表6对比可知,经过分段拟合,参数估计值的绝对误差减小到10%以内,结果中代表着油气水各相含量的估计值,是本发明所求的参数。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉该技术的人在本发明所揭露的技术范围内,可理解想到的变换和替换,都应涵盖在本发明的包含范围之内,本发明请求保护的技术内容已经全部记载在权利要求书中。

Claims (1)

1.一种核磁共振测量油气水多相流分相含率的信号分离方法,其特征在于,包括以下步骤:
1)用核磁共振方法测量到的包含油气水三相含率信息的自由衰减信号,对自由衰减信号进行经验模态分解,得到m个本征模态函数imf(t),具体步骤如下:
1.1)核磁共振方法测量的自由衰减信号用公式(1)描述:
y(t)=A·exp(-t/T)+A·exp(-t/T)+A·exp(-t/T)+n(t) (1)
其中A、A、A分别为油气水各相自由衰减信号的幅值,代表油气水各相的含量,油气水多相流分相含率的信号分离方法的目的是计算出A、A、A的估计值;T、T、T分别为油气水各相自由衰减信号的横向弛豫时间,油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s,T、T、T在核磁共振测量中用以区分各相物质;n(t)为测量中的噪声信号,是需要剔除的信号;
1.2)将测量的自由衰减信号y(t)按如下步骤进行经验模态分解:
1.2.1)分别找到测量的自由衰减信号y(t)的所有极大值点和极小值点;然后通过三次样条函数,将所有极大值点拟合出包络线e+(t),将所有极小值点拟合出包络线e-(t),按公式计算均值包络线m1(t);
<mrow> <msub> <mi>m</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <msub> <mi>e</mi> <mo>+</mo> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>e</mi> <mo>-</mo> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> <mn>2</mn> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
1.2.2)将测量信号y(t)减去m1(t)得到一个去掉低频信号的新信号
<mrow> <msubsup> <mi>h</mi> <mn>1</mn> <mn>1</mn> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>y</mi> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>m</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
1.2.3)重复步骤1.2.1)至步骤1.2.2),假设经过k次之后满足对于整个信号,极值点的个数和过零点的个数相差不超过1;并且在任意时刻,局部极大值和极小值的包络线的平均值是零;则测量信号y(t)的一阶imf分量为;
<mrow> <msub> <mi>imf</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mi>h</mi> <mn>1</mn> <mi>k</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
1.2.4)将测量信号y(t)减去imf1(t),获得一个去掉高频分量的残余分量res1(t);
res1(t)=y(t)-imf1(t) (5)
1.2.5)再对残余分量res1(t)重复步骤1.2.1)至步骤1.2.4),得到第二个分量imf2(t),如此进行多次分解一直到残余分量resn(t)是常量或单调函数时,经验模态分解过程结束,得到:
<mrow> <msub> <mi>y</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>imf</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>res</mi> <mi>n</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
1.3)将残余分量resn(t)当做最后一个imf分量,imfn+1(t)=resn(t)得到:
<mrow> <msub> <mi>y</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mi>imf</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
其中m=n+1;
2)将本征模态函数从低阶往高阶依次叠加,分别计算与测量信号的相关系数,找到使相关系数产生跳变的本征模态函数,将该本征模态函数至最后一个本征模态函数累加获得重构信号,达到滤波的目的,具体步骤如下:
2.1)在用核磁共振技术测量多相流时,对多数包含噪声的测量信号,其主要信号的能量集中在低频段,大部分噪声分布在高频段;因此存在某个imfk(t)分量,使得对于该分量之后的imf(t)中信号占主要部分,而其前k个imf(t)中噪声占主要部分;
2.2)imf(t)的信号重构函数C(j)
噪声和信号都与测量信号存在相关性,每个imf分量与测量信号y(t)也存在一定相关性,随着分量从低阶到高阶依次相加重构,重构后的信号与测量信号的相关性依次增大;对于信噪比大于0dB的测量信号,信号与测量信号的相关性比噪声与测量信号的相关性更大;因此在加入imfk+1(t)分量之后,重构信号与测量信号的相关系数会有一个跳变;通过步骤1)将测量信号y(t)分解为m个imf(t)分量,分别将分量从低阶到高阶依次相加重构得到m个函数,C(j),j=1,2,…,m;
<mrow> <mtable> <mtr> <mtd> <mrow> <mi>C</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>j</mi> </munderover> <msub> <mi>imf</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>j</mi> <mo>&amp;le;</mo> <mi>m</mi> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
2.3)C(j)和y(t)的相关系数corr(j)计算
分别计算C(j)和输入信号的相关系数corr(j),corr(j)的绝对值取值在0到1之间,越接近1表示相关性越强;计算相关系数的公式为:
<mrow> <mi>c</mi> <mi>o</mi> <mi>r</mi> <mi>r</mi> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mrow> <mo>(</mo> <mi>C</mi> <msub> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>C</mi> <mo>&amp;OverBar;</mo> </mover> <mo>(</mo> <mi>j</mi> <mo>)</mo> <mo>)</mo> </mrow> <mo>&amp;CenterDot;</mo> <mrow> <mo>(</mo> <mi>y</mi> <msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>y</mi> <mo>&amp;OverBar;</mo> </mover> <mo>(</mo> <mi>t</mi> <mo>)</mo> <mo>)</mo> </mrow> </mrow> <msqrt> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <mrow> <mo>(</mo> <mi>C</mi> <msub> <mrow> <mo>(</mo> <mi>j</mi> <mo>)</mo> </mrow> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>C</mi> <mo>&amp;OverBar;</mo> </mover> <mo>(</mo> <mi>j</mi> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>&amp;CenterDot;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msup> <mrow> <mo>(</mo> <mi>y</mi> <msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>y</mi> <mo>&amp;OverBar;</mo> </mover> <mo>(</mo> <mi>t</mi> <mo>)</mo> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
其中分别为数组C(j)、y(t)的平均值,N为采样点数;
2.4)k值计算
经过多次仿真实验,确定相关系数产生跳变的阀值为5,即当加入imfk+1(t)分量后的相关系数绝对值的增量是加入imfk(t)分量后的相关系数绝对值增量的5倍,k即为所求的值;
<mrow> <mfrac> <mrow> <mi>c</mi> <mi>o</mi> <mi>r</mi> <mi>r</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>2</mn> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mi>o</mi> <mi>r</mi> <mi>r</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow> <mrow> <mi>c</mi> <mi>o</mi> <mi>r</mi> <mi>r</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mi>o</mi> <mi>r</mi> <mi>r</mi> <mrow> <mo>(</mo> <mi>k</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>&gt;</mo> <mn>5</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
2.5)滤波后的信号y2(t)
用核磁共振技术测量多相流的自由衰减信号,如式(1)所示,其有效信号的能量集中在低频段,大部分噪声分布在高频段;据此,存在某个imfk(t)分量,使得对于该分量之后的imfk+1(t)、、、imfm(t)中信号占主要部分,imf1(t)、、、imfk(t)中噪声占主要部分,剔除这部分imf(t),就达到了滤波的效果;即信号y(t)滤波后的信号为y2(t):
<mrow> <msub> <mi>y</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msub> <mi>imf</mi> <mi>i</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
3)将滤波后的信号进行分段拟合,得到油气水各相参数的估计值,具体步骤如下:
3.1)经过滤波后的自由衰减信号y2(t)是一个数组,需要通过指数函数拟合的方法得到表达式,完成油气水三相物质的自由衰减信号分离,采用最小二乘法和模拟退火优化算法对滤波后的自由衰减信号进行拟合,步骤2)获得滤波后的自由衰减信号为:
y2(t)=A·exp(-t/T)+A·exp(-t/T)+A·exp(-t/T) (12)
将y2(t)分为t<0.3s和t>0.3s两段;
3.2)对于指数衰减函数A·exp(-t/T),当t>5T时,幅值衰减为0.0067A,近似为0;油田三相流中T为0.3s~1.0s,T为0.03~0.06s,T为0.01s~0.5s,因此当t>0.3s,气体的自由衰减信号认为衰减至0,则原三相自由衰减信号测量信号则变成油水两相,简化成两相拟合问题,结果得到:
其中分别为A、A、T、T的参数估计值;
3.3)由步骤3.2)中已经求出油、水两相参数的估计值,因此截取测量信号中t<0.3s的数据,直接拟合求出气体单项参数的估计值,得到:
其中分别为A、T的参数估计值;结果中代表着油气水各相含量的估计值。
CN201510545805.9A 2015-08-31 2015-08-31 一种核磁共振多相流测量中各相分离的方法 Expired - Fee Related CN105046025B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510545805.9A CN105046025B (zh) 2015-08-31 2015-08-31 一种核磁共振多相流测量中各相分离的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510545805.9A CN105046025B (zh) 2015-08-31 2015-08-31 一种核磁共振多相流测量中各相分离的方法

Publications (2)

Publication Number Publication Date
CN105046025A CN105046025A (zh) 2015-11-11
CN105046025B true CN105046025B (zh) 2017-11-17

Family

ID=54452566

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510545805.9A Expired - Fee Related CN105046025B (zh) 2015-08-31 2015-08-31 一种核磁共振多相流测量中各相分离的方法

Country Status (1)

Country Link
CN (1) CN105046025B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109946384A (zh) * 2019-04-08 2019-06-28 山东大学 一种基于rapid层析成像技术的信号获取过程优化方法

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105928701B (zh) * 2016-04-29 2017-07-18 石家庄铁道大学 基于相关分析的emd过程中有效imf的判定方法
CN106154344A (zh) * 2016-08-01 2016-11-23 湖南文理学院 一种基于组合滤波的大地电磁信号去噪方法
CN107766793A (zh) * 2017-09-20 2018-03-06 天津大学 基于混合型方法的mems陀螺仪信号去噪处理方法
CN108680212B (zh) * 2018-05-18 2020-01-07 中国石油天然气股份有限公司 多相流磁共振流量计刻度装置及其含水率、流速刻度方法
CN111934647A (zh) * 2020-07-14 2020-11-13 广东技术师范大学 一种emd-lms混合滤波方法及滤波器

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6046587A (en) * 1997-06-24 2000-04-04 Southwest Research Institute Measurement of flow fractions, flow velocities, and flow rates of a multiphase fluid using NMR sensing
CN102297896A (zh) * 2011-05-23 2011-12-28 浙江大学 一种多相流设备中雾沫夹带量的检测方法及装置

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6046587A (en) * 1997-06-24 2000-04-04 Southwest Research Institute Measurement of flow fractions, flow velocities, and flow rates of a multiphase fluid using NMR sensing
CN102297896A (zh) * 2011-05-23 2011-12-28 浙江大学 一种多相流设备中雾沫夹带量的检测方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
EMD与ICA算法在三相流流量软测量中的应用;于莉娜等;《PROCESS AUTOMATION INSTRUMENTATION》;20140228;第35卷(第2期);第69-72页 *
多相流相分率测量技术研究进展;吕宇玲等;《管道技术与设备》;20021230(第5期);第10-12页 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109946384A (zh) * 2019-04-08 2019-06-28 山东大学 一种基于rapid层析成像技术的信号获取过程优化方法

Also Published As

Publication number Publication date
CN105046025A (zh) 2015-11-11

Similar Documents

Publication Publication Date Title
CN105046025B (zh) 一种核磁共振多相流测量中各相分离的方法
CN103956756B (zh) 一种电力系统低频振荡模态辨识方法
CN104268883B (zh) 一种基于边缘检测的时频谱曲线提取方法
CN106096242A (zh) 一种基于改进emd分解的尾水管压力脉动综合评价方法
CN104278989B (zh) 一种获取低孔低渗储层饱和度指数的方法
CN103630742A (zh) 一种动态信号参数的获取方法
CN104165925B (zh) 随机共振的离心式压缩机半开式叶轮裂纹故障检测方法
Zhang et al. Energy operator demodulating of optimal resonance components for the compound faults diagnosis of gearboxes
CN104748961A (zh) 基于svd分解降噪和相关性eemd熵特征的齿轮故障诊断方法
CN109164489A (zh) 一种基于vmd与tk能量算子的地震流体预测方法
CN104990854B (zh) 确定束缚水饱和度的方法及装置
CN103245832A (zh) 基于快速s变换的谐波时频特性参数估计方法及分析仪
CN106441547A (zh) 一种变压器振动监测方法及装置
CN107957566A (zh) 基于频率选择奇异谱分析的磁共振测深信号提取方法
CN106198482A (zh) 基于拉曼光谱的检测保健品中是否添加有西药的方法
CN104297791A (zh) 一种基于地震优势频率的反演方法及系统
CN104678288B (zh) 基于信息熵和小波变换的开关电流电路故障字典获取方法
Chen et al. Improved VMD-FRFT based on initial center frequency for early fault diagnosis of rolling element bearing
CN109884697A (zh) 基于完备总体经验模态分解的砂砾岩沉积相地震预测方法
CN102359815A (zh) 一种基于小波分形组合的爆破振动信号特征提取方法
CN106092893A (zh) 一种光谱判别分析的波长优选方法
CN106405237A (zh) 一种应用于多通道电力系统信号中的低频振荡模式识别的分析方法
CN103487837A (zh) 拟饱含水核磁共振自旋回波信号的分解与合成方法
CN104820145A (zh) 用于对锁相放大器进行测试的测试仪及其测试方法
Ding et al. A fault feature extraction method of motor bearing using improved LCD

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171117

Termination date: 20200831