宋情洋,刘铁军,陈宜倩

(中国计量大学 计量测试工程学院,浙江 杭州 310018)

近些年,许多国内外专家在气液两相流流动规律研究领域取得丰硕的成果.特别是在流态图转变和流型识别方面.随着计算机控制理论的发展,神经网络因其操作简单,能很好地处理相关问题,成为流型检测的一个非常有用的方法.王妍芃[1]等将BP神经网络应用于气液两相流流型的识别中,BP神经网络的缺点在于其对样本的处理过程非常繁琐并且十分不灵活.陈露阳[2]等人以涡街流量计为元件诱发钝体绕流,通过管壁差压法获取气液两相流钝体绕流的尾迹波动信号,再采用集总经验模态分解对信号进行分解,利用Hilbert变换得到Hilbert边际谱,并利用最大互相关系数法对固有模态函数进行筛选,分别构造出了与体积含气率以及两相雷诺数相关的两种流型图,得到了较高流型识别率.D-Markov 模型是基于马尔科夫过程,用于复杂系统的模式识别的一种快速时间序列方法,最早由Ray在2004年提出.Chin、Ray等将该模型应用于异常检测领域中,并与传统的统计方法和神经网络方法进行对比分析,发现D-Markov模型的性能比上述两种方法要好得多[3].Khatkhate、Gupta、Ray等将D-Markov模型应用于联轴器的异常检测并取得了较好的结果[4].周迅、向馗等用D-Markov模型识别发动机零部件的疲劳裂纹的发展模式,并且发现这个模型对模式转变具有预测功能[5].

本文提出了一种基于符号动态滤波的气液两相流流型识别方法,这是一种改进数据序列系统处理信号的新方式.在对实数序列进行符号化的过程中,需要除去一些无关的特征,保留能表征流型信息的一些特征,基于这些有用特征建立符号序列模型.通过符号序列分析方法将数据序列转变为符号序列后,计算效率有了明显的提高.对符号数据分析来说,噪声的影响并不是很大,因此其鲁棒性也得到了加强.该方法是基于超声波气液两相流参数检测原理,采用一收一发式超声换能器,两相流体在管道流动过程中不断采集换能器接收端输出信号强度大小；不同流型情况下换能器输出信号强度不同,对不同流型下实时的信号强度进行统计,分析得出不同流型下信号强度的差异.将差异数据运用符号动态滤波方法处理得出流型与信号差异信号的对应关系；进而为现有的气液两相流流行识别方法提供参考[6-7].

1 符号动态滤波方法原理

声波在介质传播过程中会伴随着能量的传递,当超声波在气液两相流传播中遇到气相和液相分相界面会产生声波的反射、散射和衍射，使换能器接收到超声波幅值、相位、频率产生变化.本文通过符号动态滤波算法分析流经气液两相流后被超声波换能器接收信号幅值衰减情况，从而得到气液两相流的流型信息.对超声信号处理采用符号动态滤波法,其计算过程可以分为解析信号相空间分割和D-Markov模型两大部分[8].

1.1 解析信号相空间分割

解析信号相空间分割是符号动态滤波算法的一个非常重要的环节,它主要分为3个部分,分别为Hilbert变换[9-10]、映射和分割,最终将一维数据序列转换为符号序列.整个流程如图1.

图1 解析信号相空间分割流程图Figure 1 Flow chart of phase separation of the analytical signal

设符号集为S={s1,s2,…,sk},那么该符号集的符号熵定义为

(1)

式(1)中:Pi表示输出符号Si的概率.

为了表述符号熵随符号集大小k的变化,引入一个变量h(k),定义为

h(k)=H(k+1)-H(k).

(2)

取阈值ε,当h(k)<ε时,可得到符号集大小k的取值.

1.2 概率有限状态机

经过1.1描述的步骤后,实现了从原始的时间序列到符号序列的转换,最后获得一连串的符号序列.对于得到的符号序列,需要将其输入到D-Markov模型中进行计算输出状态转移概率矩阵,提取其左特征向量作为状态转移概率向量,从而对矩阵状态转移概率矩阵的相似性进行度量.符号处理的流程图如图2.

图2 符号序列处理流程图Figure 2 Symbol sequence processing flow chart

D-Markov模型最早运用于工程领域的一些检测中.对于给定的一串符号序列,构造一个长度为D的窗口,从最初的符号开始,将长度为D的窗口向左滑动,每次滑动窗口产生一个新的符号,将这个新产生的符号加上前一个状态的D-1个符号,构成了一个新的状态.将原状态的符号串的个数表示为N(sn-1sn-2…sn-D),将原状态加上新产生的符号后形成的新的符号串的个数表示为N(snsn-1sn-2…sn-D).

2 气液两相流流型识别实验装置设计

2.1 气液两相流流型识别装置设计

超声气液两相流流型识别实验系统流程图如图3,主要由信号发生模块、气/水源模块、测量模块和数据采集模块四个模块组成.其中,信号发生模块包括信号发生器、功率放大器,气/水源模块包括气泵、水箱、调节阀和流量计,测量模块包括透明管道和一发一收式的超声换能器,数据采集模块包括示波器、数据采集卡和PC上位机.信号发生器产生一个激励信号,经过功率放大器放大后传到超声换能器的发射端,由此生成超声信号注入到垂直管道中.同时,打开调节阀,根据调节阀的开度调整气泵和水箱分别注入到垂直管道的气液两相流量中.在垂直管道的另一端,由超声换能器的接收端接收超声波信号,在连接示波器观测波形的同时,连接数据采集卡并将接收到的数据上传至PC上位机.

图3 实验系统流程图Figure 3 Experimental system flow chart

2.2 实验测量设备及条件

本文使用的是型号为DG1022Z的信号发生器,正弦波范围1 μHz～25 MHz,正弦波幅值范围1.0 mVpp～10 Vpp.采用的功率放大器是TYH-200M中频功率放大器,可放大增益为0～10.超声换能器的型号为DYU-200K-01F,其中心频率为(200±30) kHz,量程为0.02～1 m,工作温度为-40～80 ℃,耐压值不大于0.5 MPa.示波器是tekronix生产的DPO 2024B示波器.由于本文发生的超声信号的频率为200 kHz,按照一般的选型规则,数据采集卡的采样频率应为所需采集信号频率的5～10倍,那么数据采集卡的采样频率应该要达到1 M以上,因此本文采用的数据采集卡为美国NI公司生产的USB6363 DAQ数据采集卡,16位AD分辨率,2 MS/s采样速率,10 ns时间分辨率.

2.3 实验方案设计及数据预处理

气液两相流流型识别实验主要包括产生流型、发送超声波信号以及接收采集超声信号三个部分,其中泡状流、弹状流和环状流分别进行多组实验,每组实验时长为2 min.由于纯水状态下信号并未发生衰减,所以做1组实验.

本文在垂直透明管道中设置了四种气液两相流流型,分别为纯水状态、泡状流、弹状流和环状流.纯水状态、泡状流、弹状流和环状流的衰减信号分别如图4-7.

图4 纯水状态衰减信号图Figure 4 Pure water state attenuation signal diagram

图5 泡状流衰减信号图Figure 5 Bubble flow attenuation signal diagram

图6 弹状流衰减信号图Figure 6 Slug flow attenuation signal diagram

图7 环状流衰减信号图Figure 7 Annular flow attenuation signal diagram

由图4-7可见,在不同流型下,超声换能器接收端接收到的信号强度差异比较明显.纯水状态下,超声信号并没有发生衰减;泡状流状态下,当超声信号遇到小气泡时,会产生较小程度上的衰减;弹状流状态下,当遇到弹状大气泡时,超声信号在一段时间内会产生较大程度上的衰减,遇到小气泡时,产生较小程度上的衰减;而环状流状态下,超声信号则呈现出长时间内很大程度的衰减.

3 基于符号动态滤波的识别方法

3.1 流型的特征值提取

因为符号集k的大小无关乎信号测量时间段内具体的信号波动数值,所以可以在整个实验过程中截取任意时间段、任意流型的实验数据进行分析,以此来确定k的大小.而对于阈值ε的选择既要考虑到计算准确度同时也要兼顾超声信号采集速度和噪音影响等因素.

取泡状流的一个时间段内的数据集来分析,在当h(k)<ε时,得到k的取值,本文选择ε<0.25,则得到k=5.

一般来说,阈值较小,则求得的符号集就偏大,也就是划分得更加细致,这可能会增加噪声,导致错误符号出现的概率也增大,而且计算量也会增大.如果阈值取得较大,求得的符号集就较小,也许不能较好地体现原始时间序列的所包含的关于流型的信息.所以,在选择阈值的时候,需要兼顾计算速度和准确度.

本文选用均匀分割的方式对解析信号进行分割.其中,纯水状态、泡状流、弹状流和环状流的部分数据分组结果如图8-11.

图8 纯水状态分组结果Figure 8 Pure water status grouping results

图9 泡状流分组结果Figure 9 Bubble flow grouping results

图10 弹状流分组结果Figure 10 Slug flow grouping results

图11 环状流分组结果Figure 11 Ring stream grouping results

接下来,需要将分组完成后的数据进行符号化.将0～0.2之间的数据命名为a,0.2～0.4之间的数据命名为b,0.4～0.6之间的数据命名为c,0.6～0.8之间的数据命名为d,0.8～1之间的数据命名为e.由此,得到各流型对应的符号序列.

3.2 建立流型识别模型

Ray在应用过程中发现,取D=1和D=2时,所产生的结果相近[3].所以,本文取D=1.那么状态集Q={a,b,c,d,e}.然后根据3.1小节选择的较优的符号集大小,可以构造出性能较好的D-Markov模型.选取部分实验数据作为训练样本,将训练样本输入D-Markov模型,得出可以体现每种流型信号波动情况的基准状态转移矩阵.同理,得到能代表测试样本的状态转移矩阵,与每种流型的基准状态转移矩阵进行相似性度量.相似度高说明其波动信号更接近与之相对应的流型,从而达到流型识别的目的.

3.3 相似性度量

根据状态转移矩阵,可以分别求得训练样本的左特征向量.其中纯水状态训练样本的左特征向量记为q1,将5组泡状流训练样本的左特征向量记为q21,q22,q23,q24,q25对这5个左特征向量求其欧拉距离,得到相似度,其统计如表1.

表1 泡状流的状态转移矩阵的相似度

可见,5组泡状流训练样本的状态转移矩阵之间的距离很小,分布在0.020 5～0.074 1之间.这可以说明相对应的5个状态转移矩阵之间的相似度较高.由于q21相对于其它4个左特征向量,与每个特征向量的距离最小,所以将q21作为泡状流的基准特征向量,与之对应的基准状态转移矩阵为Π21.

将5组弹状流的左特征向量记为q31,q32,q33,q34,q35分别对任意两个特征向量求欧拉距离,得到对应相似度,如表2.

表2 弹状流的状态转移矩阵的相似度

可见,5组弹状流训练样本的状态转移矩阵之间的距离较小,分布在0.054 0～0.139 2之间.这可以说明相对应的5个状态转移矩阵之间的相似度较高.由于q34相对于其它4个左特征向量,与每个特征向量的距离最小,所以将q34作为弹状流的基准特征向量,与之对应的基准状态转移矩阵为Π34.

将5组环状流训练样本的左特征向量记为q41,q42,q43,q44,q45分别对任意两个特征向量求欧拉距离,得到对应相似度,如表3.

表3 环状流的状态转移矩阵的相似度

可见,5组环状流训练样本的状态转移矩阵之间的距离较小,分布在0.037 1～0.080 1之间.这可以说明相对应的5个状态转移矩阵之间的相似度较高.由于q41相对于其它4个左特征向量,与每个特征向量的距离最小,所以将q41作为泡状流的基准特征向量,与之对应的基准状态转移矩阵为Π41.

接下来,将纯水左特征向量q1、泡状流基准特征向量q31、弹状流基准特征向量q31和环状流基准特征向量q41,进行相似度测量,分别计算它们两两之间的欧拉距离,如表4.

表4 流型状态转移矩阵的相似度

可见,4种不同流型的基准状态转移矩阵之间的欧氏距离相较于同种流型下的各组状态转移矩阵的欧氏距离要大得多.于是可以得出,在不同流型下,采集到的信号强度差异较大,而同种流型下采集到的信号强度较为相近.因此,相似性度量可以为流型识别提供一种新的思路.

4 结语

为了更加深入地认识气液两相流流动结构及相应流动参数的测量,本文搭建了一套垂直上升的气液两相流流型识别实验系统,对管径为50 mm的透明管道进行了气液两相流流型的识别实验,其中超声换能器用于获取垂直上升气液两相流的波动信号.然后运用符号动态滤波方法对实验过程采集到超声波信号分析处理,得到能表征两相流信号波动情况的状态转移矩阵,对矩阵进行相似性度量可判别不同流型的状态转移矩阵之间的差异程度,从而为以后气液两相流流型识别提供了一种可行性思路.