基于主成分分析和支持向量机的油气水多相流流型识别方法

摘要

一种基于主成分分析和支持向量机的油气水多相流流型识别方法，其特点是：利用水平管路、取压环、取压管、压力、差压变送器、数据采集卡和计算机对不同流型的流动参数信号进行采集；将采集到的压力和不同取压间距的压差波动信号进行小波包去除噪声预处理，再对上述信号进行经验模式分解，得到的信号固有模态函数，将这些固有模态函数构成一个特征矩阵X，利用主成分分析得到流型的特征向量；再将这些特征作为支持向量机的输入样本，应用支持向量机完成从特征空间到流型空间的映射，最终实现流型识别。所采用的特征向量融合了多个传感器的信息，能够全面反映流型的信息，所应用的支持向量机同神经网络分类器相比，具有学习速度快、分类能力强等特点。

说明书

技术领域

本发明涉及多相流测量技术领域，尤其涉及一种基于主成分分析和支持向量机的油气水多相流流型识别方法。

背景技术

传统的流型识别方法有两大类：一类是采用实验方法做出流型图；另一类是根据对流型转变机理的分析得到的转变准则关系式，利用现场的流动参数来确定流型。由于在生产实际中流动参数往往是需要解决的问题，因此传统的方法无法得到广泛应用。

流型的现代测量方法从工作原理上可分为直接测量和间接测量方法。常用的直接测量方法有目测法、高速摄影法等，该方法具有一定的主观性，对于多个观察者可能得出不同的结论。间接测量方法是利用测量的仪器设备测量的压力、压差和含气率等波动信号，通过对信号进行分析、提取特征，结合神经网络等模式识别技术来识别流型。该方法有两个关键问题：一是流型特征向量的提取；二是网络模型的选取，其中特征提取尤为关键。现有的方法主要是利用单一的压力或压差波动信号来提取信号的单一特征，导致提取的特征无法全面反映流型的信息，致使流型识别率不高。

发明内容

本发明的目的在于克服上述现有技术的缺陷，提出一种流型识别准确性、可靠性高，通用性强的基于主成分分析和支持向量机的油气水多相流流型识别方法。

本发明的目的是由以下技术方案来实现的：一种基于主成分分析和支持向量机的油气水多相流流型识别方法，其特征在于，它包括以下步骤：

(1)对不同流型的流动参数信号进行采集：利用取压间距分别为5倍、10倍和15倍管径D；采样频率为256Hz；采样点数为1024点；1个压力变送器和3个不同取压间距的差压变送器对流动参数信号进行实时采集；

(2)去除噪声预处理：利用小波包对压力和压差波动信号进行去除噪声处理。

(a)对1个压力波动和3个不同取压间距的压差波动信号S₁、S₂、S₃和S₄按照下面的递归式(1)进行4层小波包分解：

$>\left(\begin{array}{}>>>>u>>2>n>>>>(>t>)>>=>>2>\underset{}{\underset{}{}}\end{array}\right)$

式中h(k)为高通滤波器组；g(k)为低通滤波器组；u₀(t)＝(t)，为尺度函数；u₁(t)＝ψ(t)，为小波函数；k为分解的尺度。h(k)和g(k)满足正交关系：

                g(k)＝(-1)^kh(1-k)                              (2)

(b)对分解的各个频带分别设定阈值，采用如下的浮动阈值形式

$>>>t>k>>=>>21>n>>(>n>)>>>\sigma >/>>n>>->->->>(>3>)>>>s>$

其中σ为噪声强度，用中位数MAD(υ_k)/0.6745(υ_k为小波包分解系数)来估计，n为信号的长度。

(c)用非线性函数ζ＝sgn(υ_k)(|υ_k|-t_k)作用于信号的小波包分解系数υ_k，得到系数υ_k，即将浮动阈值以上的系数减去阈值后并保留下来。

(d)由系数υ_k重建原信号，从而达到去除噪声，得到去除噪声后的信号为S₁′、S₂′、S₃′和S₄′。

(3)在对去除噪声后的1个压力和3个压差波动信号S₁′、S₂′、S₃′、S₄′进行经验模式分解，提取流型的特征向量。EMD方法从本质上讲是对一个信号进行平稳化处理，其结果是将信号中不同尺度的波动或趋势逐级分解出来，产生一系列具有不同特征尺度的数据序列，每一个序列称为一个固有模态函数IMF。IMF分量必须满足下面两个条件：一是其极点数和零点数相同(或最多相差一个)，二是其上下包络线关于时间轴局部对称。

EMD方法步骤如下：假设任何信号都由不同的IMF组成，每个IMF可以是线性的，也可以是非线性的，这样任何一个信号就可以分解为有限个IMF之和，则IMF可以按以下方法“筛选”(Sift)：

(a)确定信号x(t)的所有局部极值点，将所有极大值点用三次样条线连接起来形成上包络线，将所有极小值点用三次样条线连接起来形成下包络线，这两条包络线包络了所有的信号数据。

(b)将两条包络线的均值记为μ₁，求出

                 y₁(t)＝x(t)-μ₁                     (4)

(c)判断y₁(t)是否为IMF，若y₁(t)不满足IMF条件，则将y₁(t)作为原始数据，此时，记y₁(t)＝c₁(t)，则c₁(t)为信号x₁(t)的第一个IMF分量，它代表信号x₁(t)中最高频率的分量；

(d)将c₁(t)从x₁(t)中分离出来，即得到一个去掉高频分量的差值信号r₁(t)，即有

                 r₁(t)＝x₁(t)-c₁(t)                  (5)

将r₁(t)作为原始数据，重复步骤(1)、(2)和(3)，得到第二个IMF分量c₂(t)，重复n次，得到n个IMF分量。这样就有

$>\left(\begin{array}{}>>>>r>1>>->>c>2>>=>>r>2>>>>>>M>>>>>>r>>n>->1>>>->>c>n>>=>>r>n>>>\end{array}\right)$

当c_n(t)或r_n(t)满足给定的终止条件(通常使r_n(t)成为一个单调函数)时，循环结束，由(2)式和(3)式可得到：

$>>x>>(>t>)>>=>>\Sigma >>i>=>1>>nover>>>c>i>>>(>t>)>>+>>r>n>>>(>t>)>>->->->>(>7>)>>>s>$

式中，r_n(t)为残余函数，代表信号的平均趋势。而各IMF分量c₁(t)，c₂(t)，Λ，c_n(t)分别包含了信号从高到低不同频率段的成分，每一频率段所包含的频率成分都是不同的，且随信号本身的变化而变化。

(e)特征向量的提取。1个压力波动信号和3个压差波动信号的n个IMF分量分别记为c₁₁(t)，c₁₂(t)，Λ，c_1n1(t)、c₂₁(t)，c₂₂(t)，Λ，c_2n2(t)、c₃₁(t)，c₃₂(t)，Λ，c_3n3(t)、c₄₁(t)，c₄₂(t)，Λ，c_4n4(t)。将这些IMF分量组成一个特这个矩阵X，即

$>>>X>>p>\times >n>>>=>=>\left(\begin{array}{}>>>>x>11>>>>>x>12>>>>\Lambda >>>>x>>1>p>>>>>>>>x>21>>>>>x>22>>>>\Lambda >>>>x>>2>p>>>>>>>M>>>M>>>M>>>M>>>>>>x>>n>1>>>>>>x>>n>1>>>>>\Lambda >>>>x>np>>>\end{array},\left(\begin{array}{}>>>>c>11>>>>>>M>>>>>>c>>1>>n>1>>>>>>>>>c>21>>>>>>M>>>>>>c>>2>>n>2>>>>>>>>>c>31>>>>>>M>>>>>>c>>3>>n>3>>>>>>>>>c>41>>>>>>M>>>>>>c>>4>>n>4>>>>>\end{array}\right)\right)$

式中：X为p×n的特征矩阵，p为采样点数，n＝n₁+n₂+n₃+n₄为4个原始信号的所有IMF分量的个数，x_i(i＝1，2，Λ，n)为p维向量。

(f)对矩阵X用下式变换为相关矩阵，即主成分矩阵R，并计算出相关矩阵R的m个特征值。

①将原始变量进行标准化

各指标量纲往往不同，分析时，不同的量纲和数量级会引出新的问题，故要标准化，使标准化的变量平均值为0，方差为1，即 $>sup>>x>ij>*sup>>=>>>>x>ij>>->ver>>x>‾>>i>>>>var>>(>>x>ij>>)>>>>>s>i＝1，2，\Lambda ，n；j＝1，2，\Lambda ，p (9)$

式中：x_ij^*为标准化后的观测量，x_i为平均值，为标准差。

②计算X的相关矩阵R

标准化后，相关阵和协方差阵相同，且R的子分量可表示为

$>>>r>ij>>=>>1>n>>>\Sigma >>i>=>1>>nover>sup>>x>ti>*sup>sup>>x>ij>*sup>>->->->>(>10>)>>>s>$

式中：x_ti^*，x_ij^*分别为特征矩阵中第t行第i列与第j列元素。

③求相关矩阵R的特征值及特征向量

计算出各特征值(λ₁＞λ₂＞Λ＞λ_n)和各特征向量(u₁，u₂，Λ，u_n)。

④选择主成分

计算各主成分的贡献率及累计贡献率，主成分的个数视具体问题而定，本发明取累计贡献率为90％所对应的m个变量作为主成分。

⑤计算得分

根据主成分表达式计算各主成分得分。在主成分分析中， $>>>\mu >j>>=>>\lambda >j>>/>>\Sigma >>i>=>1>>mover>>>\lambda >i>>>s>为主成分的贡献率，比值反映前k个主成分所代表的原始数据信息的百分数。令$

                 T＝[μ₁，μ₂，Λ，μ_m]               (11)

则T为不同流型的特征向量。

(4)支持向量机分类器。采用支持向量机来实现流型的识别。水平管内油气水多相流流型有4种，即分层流、间歇流、泡状流和环状流。根据上述分类情况，采用“一对多”的策略，构造4个二值分类器，对每一个流型构造一个支持向量机，4个支持向量机的输出是一个4维向量，每一个分量表示样本是否对应该流型，训练好的支持向量机可实现对流型的识别。

本发明基于主成分分析和支持向量机的油气水多相流流型识别方法的优点在于：

1.对1个压力和3个差压变送器的信号进行融合，提取的特征能全面反映流型的信息，因此，识别准确性和可靠性高。

2.采用支持向量机来实现流型的识别，计算速度快，有利于流型的在线识别。

附图说明

图1是基于主成分分析和支持向量机的油气水多相流流型识别方法流程图。

图2是基于主成分分析和支持向量机的油气水多相流流型信号采集装置的结构示意图。

具体实施方式

下面利用附图和实施例对本发明基于主成分分析和支持向量机的油气水多相流流型识别方法作进一步说明。

如图2所示，基于主成分分析和支持向量机的油气水多相流流型信号采集装置，具有多相流水平管路16，在多相流水平管路16上套置的取压环1a、2b、3c、4d、5e、6f、7g，取压环1a与取压管8a连接，取压环2b、3c与取压管8b连接，取压环4d、5e与取压管8c连接，取压环6f、7g与取压管8d连接，取压管8a与压力变送器9a连接，取压管8b、8c、8d分别与差压变送器10b、11c、12d连接，压力变送器9a、差压变送器10b、11c、12d均与数据采集卡13电连接，数据采集卡13与计算机14电连接。所述取压环2b、3c的间距为5倍管径D，取压环4d、5e的间距为10倍管径D，取压环6f、7g的间距为15倍管径D。其中取压环1a距离入口300mm，取压环1a与取压环2b相隔200mm，取压环3c与取压环4d相隔200mm，取压环5e与取压环6f相隔200mm。压力变送器9a、差压变送器10b、11c、12d型号为PD-23的输出标准信号与型号为IMP3595的数据采集卡13相连，数据采集卡13与计算机14相连，由计算机14完成压差信号的数据处理，进而实现对油气水多相流流型信号的采集。计算机14还可连接打印机15。基于主成分分析和支持向量机的油气水多相流流型信号采集装置的软件程序是依据自动检测技术及计算机数据处理技术编制的，为本领域技术人员所熟悉的技术。

如图1所示，本发明基于主成分分析和支持向量机的油气水多相流流型识别方法包括下述步骤：

(1)对不同流型的流动参数信号进行采集：在基于主成分分析和支持向量机的油气水多相流流型信号采集装置中，利用取压间距分别为5倍、10倍和15倍管径D；采样频率为256Hz；采样点数为1024点，采样时间为16s；用1个压力变送器9a和3个不同取压间距的差压变送器10b、11c、12d对流动参数信号进行实时采集；

(2)应用小波理论对采集的压差波动信号进行去除噪声预处理：(a)对1个压力波动和3个不同取压间距的压差波动信号S₁、S₂、S₃和S₄按照下面的递归式(1)进行4层小波包分解：

$>\left(\begin{array}{}>>>>u>>2>n>>>>(>t>)>>=>>2>\underset{}{\underset{}{}}\end{array}\right)$

式中h(k)为高通滤波器组；g(k)为低通滤波器组；u₀(t)＝(t)，为尺度函数；u₁(t)＝ψ(t)，为小波函数；k为分解的尺度。h(k)和g(k)满足正交关系：

                 g(k)＝(-1)^kh(1-k)                     (2)

(b)对分解的各个频带分别设定阈值，采用如下的浮动阈值形式

$>>>t>k>>=>>21>n>>(>n>)>>>\sigma >/>>n>>->->->>(>3>)>>>s>$

其中σ为噪声强度，用中位数MAD(υ_k)/0.6745(υ_k为小波包分解系数)来估计，n为信号的长度。

(c)用非线性函数ζ＝sgn(υ_k)(|υ_k|-t_k)作用于信号的小波包分解系数υ_k，得到系数υ_k，即将浮动阈值以上的系数减去阈值后并保留下来。

(d)由系数υ_k重建原信号，从而达到去除噪声的目的。

(3)提取流型的特征：运用EMD对采集到的压差波动信号x(t)按如下步骤进行分解：

(a)确定信号x(t)的所有局部极值点，将所有极大值点用三次样条线连接起来形成上包络线，将所有极小值点用三次样条线连接起来形成下包络线，这两条包络线包络了所有的信号数据。

(b)将两条包络线的均值记为μ₁，求出

                 y₁(t)＝x(t)-μ₁                 (4)

(c)判断y₁(t)是否为IMF，若y₁(t)不满足IMF条件，则将y₁(t)作为原始数据，此时，记y₁(t)＝c₁(t)，则c₁(t)为信号x₁(t)的第一个IMF分量，它代表信号x₁(t)中最高频率的分量；

(d)将c₁(t)从x₁(t)中分离出来，即得到一个去掉高频分量的差值信号r₁(t)，即有

                 r₁(t)＝x₁(t)-c₁(t)               (5)

将r₁(t)作为原始数据，重复步骤(1)、(2)和(3)，得到第二个IMF分量c₂(t)，重复n次，得到n个IMF分量。这样就有

$>\left(\begin{array}{}>>>>r>1>>->>c>2>>=>>r>2>>>>>>M>>>>>>r>>n>->1>>>->>c>n>>=>>r>n>>>\end{array}\right)$

当c_n(t)或r_n(t)满足给定的终止条件(通常使r_n(t)成为一个单调函数)时，循环结束，由(2)式和(3)式可得到：

$>>x>>(>t>)>>=>>\Sigma >>i>=>1>>nover>>>c>i>>>(>t>)>>+>>r>n>>>(>t>)>>->->->>(>7>)>>>s>$

式中，r_n(t)为残余函数，代表信号的平均趋势。而各IMF分量c₁(t)，c₂(t)，Λ，c_n(t)分别包含了信号从高到低不同频率段的成分，每一频率段所包含的频率成分都是不同的，且随信号本身的变化而变化。

(e)特征向量的提取。1个压力波动信号和3个压差波动信号的n个IMF分量分别记为c₁₁(t)，c₁₂(t)，Λ，c_1n1(t)、c₂₁(t)，c₂₂(t)，Λ，c_2n2(t)、c₃₁(t)，c₃₂(t)，Λ，c_3n3(t)、c₄₁(t)，c₄₂(t)，Λ，c_4n4(t)。将这些IMF分量组成一个特这个矩阵X，即

$>>>X>>p>\times >n>>>=>=>\left(\begin{array}{}>>>>x>11>>>>>x>12>>>>\Lambda >>>>x>>1>p>>>>>>>>x>21>>>>>x>22>>>>\Lambda >>>>x>>2>p>>>>>>>M>>>M>>>M>>>M>>>>>>x>>n>1>>>>>>x>>n>1>>>>>\Lambda >>>>x>np>>>\end{array},\left(\begin{array}{}>>>>c>11>>>>>>M>>>>>>c>>1>>n>1>>>>>>>>>c>21>>>>>>M>>>>>>c>>2>>n>2>>>>>>>>>c>31>>>>>>M>>>>>>c>>3>>n>3>>>>>>>>>c>41>>>>>>M>>>>>>c>>4>>n>4>>>>>\end{array}\right)\right)$

式中：X为p×n的特征矩阵，p为采样点数，n＝n₁+n₂+n₃+n₄为4个原始信号的所有IMF分量的个数。

(f)对矩阵X用下式变换为相关矩阵，即主成分矩阵R，并计算出相关矩阵R的m个特征值。

①将原始变量进行标准化

各指标量纲往往不同，分析时，不同的量纲和数量级会引出新的问题，故要标准化，使标准化的变量平均值为0，方差为1，即

$>sup>>x>ij>*sup>>=>>>>x>ij>>->ver>>x>‾>>i>>>>var>>(>>x>ij>>)>>>>>s>i＝1，2，\Lambda ，n；j＝1，2，\Lambda ，p (9)$

式中：x_ij^*为标准化后的观测量，x_i为平均值，为标准差。

②计算X的相关矩阵R

标准化后，相关阵和协方差阵相同，且R的子分量可表示为

$>>>r>ij>>=>>1>n>>>\Sigma >>i>=>1>>nover>sup>>x>ti>*sup>sup>>x>ij>*sup>>->->->>(>10>)>>>s>$

式中：x_ti^*，x_ij^*分别为特征矩阵中第t行第i列与第j列元素。

③求相关矩阵R的特征值及特征向量

计算出各特征值(λ₁＞λ₂＞Λ＞λ_n)和各特征向量(u₁，u₂，Λ，u_n)。

④选择主成分

计算各主成分的贡献率及累计贡献率，主成分的个数视具体问题而定，本发明取累计贡献率为90％所对应的m个变量作为主成分。

⑤计算得分

在主成分分析中， $>>>\mu >j>>=>>\lambda >j>>/>>\Sigma >>i>=>1>>mover>>>\lambda >i>>>s>为主成分的贡献率，比值反映前k个主成分所代表的原始数据信息的百分数。令$

               T＝[μ₁，μ₂，Λ，μ_m]              (11)

则T为不同流型的特征向量。

(4)采用支持向量机来实现流型的识别：水平管内油气水多相流流型有4种，即分层流、间歇流、泡状流和环状流。根据上述分类情况，采用“一对多”的策略，构造4个二值分类器，分别表示为SVM1、SVM2、SVM3、SVM4。如图2所示。对每一个流型构造一个支持向量机，每个支持向量机的输出y_i∈{-1，+1}，(i＝1，2，3，4)，则4个支持向量机的输出是一个4维向量，每一个分量表示样本是否对应该流型，(1，-1，-1，-1)代表分层流，(-1，1，-1，-1)代表间歇流，(-1，-1，1，-1)代表泡状流，(-1，-1，-1，1)代表环状流。用120个(4种流型各30个样本)对支持向量机进行训练，通过比较，确定本发明的核函数为径向基核函数，即 $>>K>>(>x>,>y>)>>=>>e>>->>>|>|>x>->y>|>|>>2>>/>>>2>\sigma >>2>>>>,>>s>式中\sigma 取0.05。训练好的支持向量机可实现对流型的识别。再用200个样本进行测试，结果如表1所示(部分测试样本)。$

                                                           表1

                                              测试样本  输出  流型  0.47561  0.29038  0.11616  0.08258  0.02529  0.00434  0.00473  (-1，-1，-1，1)  环状流  0.48657  0.28142  0.10764  0.07964  0.02637  0.00538  0.00456  (-1，-1，-1，1)  环状流  0.48841  0.27038  0.11524  0.08357  0.02735  0.00469  0.00504  (-1，-1，-1，1)  环状流  0.57743  0.31873  0.06160  0.01991  0.00848  0.00563  0.00469  (-1，-1，1，-1)  泡状流  0.57628  0.31772  0.06036  0.01814  0.01584  0.00433  0.00297  (-1，-1，1，-1)  泡状流  0.57628  0.31772  0.06036  0.01814  0.01584  0.00433  0.00297  (-1，-1，1，-1)  泡状流  0.69828  0.11918  0.09776  0.06454  0.01191  0.004421  0.00229  (-1，1，-1，-1)  间歇流  0.69828  0.11918  0.09776  0.06454  0.01191  0.004421  0.00229  (-1，1，-1，-1)  间歇流  0.68550  0.23403  0.05095  0.02419  0.00263  0.001564  0.00115  (-1，1，-1，-1)  间歇流  0.47509  0.28644  0.13284  0.07738  0.01721  0.005168  0.00364  (1，-1，-1，-1)  分层流  0.42308  0.21099  0.16566  0.10569  0.074812  0.01078  0.00551  (1，-1，-1，-1)  分层流  0.41356  0.23654  0.15324  0.12341  0.06244  0.01025  0.00631  (1，-1，-1，-1)  分层流