基于多域特征提取的气液两相流流型识别

2023-11-03张立峰

计量学报 2023年10期

张立峰, 王智

(华北电力大学自动化系，河北保定 071003)

1 引言

气液两相流系统作为一个高度随机、非线性的瞬时动态系统,广泛存在于自然环境与工业生产中,如火力发电、核电、油气运输、航空航天等领域[1,2]。两相界面不稳定的热量、动量及质量传递过程使得流动机理复杂,难以用精确的数学模型描述。流型是两相流的流态特征之一,准确实现流型识别对定量塑造两相流流动、深入理解流动机理及工业设施的安全生产具有重要意义,然而目前尚未见合适的流动模型,对工艺参数的处理方法也不完善,因此,急需探索一种鲁棒性强、准确率高的流型识别框架。流型测量方法大致可分为射线法[3]、光学法[4]、超声波法[5]、差压法[6]、电学法等,其中电学法可以获得管道截面的相分布状态,且具有成本低、灵活性高等优点,被认为是一种有前途的多相流检测技术。电学法主要包括电阻层析成像[7,27](electrical resistance tomography, ERT)、电容层析成像[8](electrical capacitance tomography, ECT)及电磁层析成像[9](electromagnetic tomography, EMT),其中ERT在气液两相流检测中得到了广泛应用,为研究流体行为提供了一种非侵入、非放射性的检测方法,有助于研究不同尺度与极端条件下的动态流动行为。ERT系统通过多点测量提供了较为全面的流动信息,使用特定的算法可以实现截面电导率分布的实时重建[10]。若能利用ERT系统获取流动数据,进而即可通过数据处理、特征提取与分类模型实现流型识别。

对于传感器测量值序列的处理方法有功率谱分析[11]、概率密度分析[12]、多尺度熵[13]、复杂网络[14]、时频分析[15]、统计分析[16]等,但其中大部分方法仅仅从时域或频域单空间中认识测量数据,尽管时频分析法兼顾了时域与频域空间,但其忽视了ERT系统可以重建截面相分布的优点,因此,为充分利用ERT系统测量信息,本文将重建图像作为测量数据的补充,在多域空间中处理ERT系统采集的流动信息。ERT系统测量值序列蕴含着介质分布变化的时域信息,重建图像提供了介质分布在空间域的直观表达,而变换域有助于提取更加全面的图像信息。离散傅里叶变换(discrete fourier transform, DFT)常被用于实现图像由空域到频域的变换,但DFT在计算中涉及到了复数运算,作为频域的扩展,Walsh变换[17]为在列率域中研究频率提供了理论依据,且具有计算过程简单、能量集中的优点。然而,提取的高维特征往往包含着冗余信息,采用合适的方法实现特征简约对提高计算效率、提升识别准确率有重要作用。

流形学习[18]方法已成功用于寻找高维空间中数据流形结构的低维嵌入映射,从而实现维数简约,广泛用于生物医学[19]、故障诊断[20]、图像处理[21]等领域。流形学习方法可分为线性与非线性两类,线性方法如独立成分分析 (independent component analysis, ICA)、主成分分析 (principal components analysis, PCA)等,非线性方法有局部线性嵌入[22](locally Linear Embedding, LLE)、等度量映射 (isometric Mapping, Isomap)、统一流形逼近与投影[23](uniform manifold approximation and projection, UMAP)等。其中,UMAP基于k-近邻思想,通过计算成对点在高维空间与低维映射空间中的距离度量并使用随机梯度下降法(stochastic gradient descent, SGD)优化结果实现降维,具有保留全局结构、计算快速及非线性降维等优点。

为实现垂直上升管道气液两相流流型识别,本文提出了一种多域特征的流型识别方案:1)在内径50 mm的垂直上升管道上进行动态试验操作,并利用ERT系统监测两相流流动以获取气液两相介质分布变化信息;2)为融合ERT传感器测量数据与ERT重建图像蕴含的流动分布信息,提出了一种多域特征处理方案,从时域、空域及列率域提取并量化流型特征。使用Walsh变换代替时频变换将图像由空域转换至列率域简化计算过程;3)为实现特征简约,使用UMAP算法对特征进行降维处理,以获得低维空间中具有代表性的特征向量,并与ICA、LLE及Isomap降维方法进行对比,验证UMAP的有效性;4)搭建支持向量机(support vector machine, SVM)实现流型辨识,并进一步研究低维映射空间维度对辨识性能的影响。对上述方案进行实验验证,结果表明,该流型识别框架具有较高的稳定性与识别效果。

2 ERT系统原理与两相流实验

2.1 电阻层析成像

电阻层析成像基于不同介质(相)电导率不同的原理。通过1对或多对电极向被测区域施加激励以形成敏感场,当场内介质分布改变时,电导率分布发生变化从而引起电势分布变化,这种改变可以通过测量边界电极输出所反应。本文所用16电极ERT传感器结构如图1所示,管道内径50 mm,圆形电极直径5 mm,均匀排列在管壁上。

图1 ERT传感器结构Fig.1 Structure of ERT sensor

ERT系统工作模式如图2所示,采用电压激励电流检测的单电极激励、测量传感策略。激励信号为频率256 kHz、峰-峰值4 V的正弦交流电压。ERT系统工作时,首先向电极E1施加激励,依次测量其余电极输出,然后将激励施加到电极E2,仍依次测量其余电极输出,重复上述过程直到完成所有测量。1次完整扫描可以获得16×15=240个测量数据,ERT系统数据采集速度为120 帧/s。

ERT图像重建是由测量值到电导率分布的反演过程,可以描述为式(1),其中I表示测量值,W为灵敏度矩阵,反应电导率分布g与边界测量值的映射关系。

I=Wg

(1)

线性反投影[24](linear back projection, LBP) 算法由于其快速性已广泛用于电阻层析成像的在线实时重建,可以表达为式(2),

(2)

式中uI中所有元素为1。

2.2 两相流实验

在华北电力大学成像与智能感测实验室的气液两相流实验装置上进行实验,实验装置示意图如图3所示。气、水输送管道由不锈钢管与透明有机玻璃管组成,其中玻璃管安装在透明直管段,用于观测流型,记录流动图像。测试段安装16电极ERT传感器并连接ERT系统,通过计算机保存数据。

图3 实验装置Fig.3 Experimental device

实验介质为空气与水。空气通过空气压缩机压缩后存储在储气罐内,并依靠稳压阀将气体出口压强固定在0.5 MPa。存储在水箱内的水由水泵抽出,与气体混合后继续向上运动,在测试段完成数据采集后返回水箱构成水循环回路,同时气体排放至环境中。分别使用电磁流量计与浮子流量计监测气液两相流量,准确度为±0.5%与±1.5%。

参考Taitel[25]流型图改变气液流量,同时通过透明直管段观测流动图像,利用流动结构的形态学特征确定流型。气体以小气泡形式离散的分布在液相中为泡状流;随着气相流速的增加,气泡间的碰撞行为加剧,小气泡聚合成直径稍大的气泡,偶尔出现直径比管径略小的“Taylor”气泡,系统由泡状流向段塞流过渡;进一步增大气相流速,气泡以近于管道内径大小的大型弹状气塞形式出现,气塞与液塞交替发生并具有一定的伪周期性趋势,流动系统发展为段塞流;当气塞表面液膜变得不稳定时,长气塞逐渐破碎,流动系统展现出了一定的震荡特性,发展为严重段塞流。

使用ERT系统采集123种不同工况点流动数据,其中泡状流13组,泡状-段塞过渡流型20组,段塞流42组,严重段塞流48组,每种工况点采集10次,每次采集10 s,流动图像分别如图4(a)～图4(d)所示。

图4 气液两相流流型Fig.4 Flow pattern of gas-liquid two-phase flow

为获得流动中两相分布的实时变化情况,依据式(2)使用LBP算法反演测量数据,图5展示了试验中基于该算法的ERT系统重建图像,图中像素灰度值较高部分代表气相,灰度值较低部分表示液相。

图5 LBP算法重建图像Fig.5 Image reconstruction based on LBP algorithm

3 流型识别框架

基于ERT系统采集的流动数据,采用多域特征提取的流型识别方案,如图6所示。该方案从测量值序列与LBP重建图像两个角度分析,由时域、空域、列率域提取量化特征,将数据集划分为训练集与测试集两类,为实现特征简约,使用UMAP算法对训练集高维特征做降维处理以获取低维流形嵌入空间,并使用训练集降维数据训练SVM分类器,最后使用测试集数据验证分类效果。

图6 流型识别框架Fig.6 Structure of flow pattern recognition

3.1 Walsh-Hadamard变换

Walsh变换[17]可以将函数映射到正交、归一化的Walsh空间,实现了图像由空域到列率域的转换。不同于Fourier变换中频率的概念,列率是一种广义的频率,其代表信号单位时间内跨越零点次数的1/2。Walsh函数系的取值仅有-1与+1,由于其运算只涉及实数的加减而没有复数的乘法,因此相对于Fourier变换更为快速简单。

对于二维离散图像信号G(x,y),x∈[0,M-1],y∈[0,N-1],其Walsh变换为:

(3)

式中:X为图像G的列率谱;u∈[0,M-1];v∈[0,N-1];φ(x,y,u,v)为Walsh基函数。

当M=2m,N=2n时,可采用具有简单递推关系的Hadamard排列作为基函数。即

(4)

式中bi(x) 为x倒序自然二进制码的第i位,仅可取0或1。

对于2k阶Hadamard排列可以表示为:

(5)

式中⊗表示克罗内克积。

3.2 特征量化

对ERT系统采集的高维测量值序列降维处理,将每1帧数据压缩至一维以简化分析:

(6)

式中:Iij为气液两相流时第i次激励下第j次测量的电流值;Iij0为Iij在管道内充满水时的测量数据;则I可表示为1次完整扫描下测量电极对气液两相介质分布变化的平均响应。

进一步提取时间序列时域特征并量化。使用均值tμ表征数据集中趋势,使用标准差tσ表征测量数据离散程度,使用偏度tskew表征数据分布的对称性,使用峰度tkurt表征数据分布的平滑、陡峭程度。

(7)

(8)

(9)

(10)

式中:T为时间序列长度;I(i)为第i次激励的电流值。

(11)

式中Gi为LBP算法重建出的第i幅图像。

提取纹理均值Sμ表示纹理分布的全局均匀程度,使用纹理方差Sσ表示纹理分布局部变化的剧烈程度,使用图能量SE表示纹理分布的粗粒程度。

(12)

(13)

(14)

式中g为像素点个数。

为更全面地反应图像的纹理特征,进一步将其转化至列率域。由于LBP算法得到的图像为仅有812个像素点的不规则圆形,因此在进行Walsh-Hadamard变换之前需对其补0填充至32×32的形状,记为G′。提取列率谱系数矩阵X′的一阶矩Cμ、L1范数CL1与L2范数CL2描述纹理全局特征。

X′=Wal(G′)

(15)

(16)

(17)

(18)

式中:Wal表示Walsh-Hadamard变换;j=1,2,…,32;λ为X′TX′的最大特征值。综上则可得到每个样本的特征向量F=[tμ,tσ,tskew,tkurt,Sμ,Sσ,SE,Cμ,CL1,CL2]。

3.3 UMAP

UMAP算法[23]基于黎曼几何与代数拓扑理论,常用于高维特征的降维与可视化,其能在降维过程中最大程度保留数据原始特征,且运算快速。下面为该方法的主要步骤描述。

给定RM维空间中含有n个样本的向量集X={xi},i=1,2,…,n,寻找每个样本点xi的k个近邻点集合ηi。对于每个点xi,寻找其最近邻点与距离:

ρi=min(d(xi,xj)|xj∈ηi,d(xi,xj)>0)

(19)

(20)

式中:ρi为xi与第1个最近邻点的距离;d(xi,xj) 为点xi与xj之间的距离;σi为xi最近邻点的半径。进一步便可得到无向加权图G,其邻接矩阵B定义为:

B=A+AT-A°AT

(21)

式中°表示哈达玛积。A中元素由相应的有向图的权重给出:

(22)

式中w((xi,xj)) 为点xi与xj的连接权重。

UMAP在低维空间中使用力定向图布局算法,通过在迭代中向边或顶点上施加作用力与反作用力实现,与模拟退火算法中所使用的方法相似。设yi与yj是高维空间中点xi与xj的低维表示,则yi与xj间的作用力定义为:

(23)

式中a、b为超参数。

排斥力定义为:

(1-w((xi-xj)))×(yi-yj)

(24)

式中ε为一个极小数,默认为0.001。

3.4 支持向量机

SVM[26]在解决小样本、非线性问题方面具有优势,其思想为在数据空间中构造最优超平面实现分类。采用成对分类法构建多分类模型,对于含有L个类别的数据集,该策略在每两个类别之间构建一个二分类SVM进而实现多分类,则共需训练L(L-1)/2个二分类器。

给定含有n个样本的训练集X={(x1,y1),(x2,y2),…,(xn,yn)},xi∈Rk,i=1,2,…,n,k为特征空间维度,yn={1,2,…,L},表示样本类别。

对于一个二分类模型,使用高斯核函数时判决函数为:

(25)

(26)

(27)

式中:αi为拉格朗日乘数;b为阈值;σ为可选参数。

4 结果分析

将本文所提流型识别模型应用于辨识垂直上升管道的气液两相流流型,将样本特征归一化至[0,1]之间作为模型输入,模型输出为代表流型的标签值,即泡状流、泡状-段塞过渡流型、段塞流及严重段塞流的标签值分别为1、2、3、4。设置训练集与测试集比例为6∶4。使用分类准确率定量评价模型识别效果:

(28)

式中:SACC,i为类别i的分类准确率;TP,i为类别i分类正确的样本个数;Ti为类别i的测试集样本总数。

为验证UMAP算法有效性,使用LLE、Isomap、ICA作为对比,同时研究了降维维度对训练集分类准确的影响,如表1所示。此时LLE算法邻居数量设置为15,UMAP与Isomap算法的邻居数量均为5。

表1 不同分类模型在训练集上的分类准确率Tab.1 Classification accuracy of different classification models on the training set (%)

由表1可以发现UMAP-SVM的分类性能在训练集上展现出最优效果,且降维空间的维度对分类结果影响较小。观察其他分类模型可以发现保留的流形空间维度并不是越高越好,过高的维度可能会带来无效特征,这不仅不利于降低模型计算量,甚至会导致分类准确率略有下降。为确定UMAP算法的最佳近邻数量,并保证其有一定的鲁棒性,将流形空间设置为2维,通过改变近邻点数量观察UMAP-SVM分类器在测试集与训练集上的分类效果,结果如图7所示。

由图7可以发现,邻居数量对UMAP-SVM分类器在训练集上的分类性能影响较小,但会对模型泛化能力产生较大影响,在邻居数量为7时模型鲁棒性达到最佳效果,此时模型对测试集数据集的分类结果如表2所示,混淆矩阵如图8所示。

表2 UMAP-SVM框架测试集分类准确率Tab.2 Classification accuracy of UMAP-SVM framework on test set

图8 UMAP-SVM辨识结果Fig.8 Identification results of UMAP-SVM

由表2可以发现UMAP-SVM流型识别模型在每种流型上的分类准确率均可达到94.8%以上,平均分类准确率可达95.7%,拥有较好的辨识性能,且具有较高的稳定性。观察图8可以发现,由于间歇流系统稳定性较差,湍流动能相对泡状流更强,参数变化范围更广,易于其他流型产生交叉重叠,因此对间歇流的识别准确率相对于泡状流要低,特别在严重段塞流时,液膜不稳定使得气塞破碎现象相对明显,系统具有一定的震荡特性,使得测量参数与其他流型的重叠范围更广,因此其余流型样本更易错分为严重段塞流,这种现象可以通过增加训练集样本数量,建立更加完备的样本集来减弱,同时使用合适的去噪方法降低测量值序列噪声也对提高模型辨识性能有利。