多模态动态核主成分分析的气液两相流状态监测

董峰 李昭 李凌涵 张淑美

气液两相流是一种气体和液体同时存在且具有明确分界面的流体流动形式,广泛存在于自然界及能源、动力、石油、化工、冶金、制药等众多工业领域中,具有流动结构复杂、过程状态波动、过程变量多及状态的变化具有随机性、非线性等特点.由于气液两相流流动状态随环境、各分相流体性质及相间的相互作用而改变,其过程参数的检测和流动状态的监测相比单相流具有更大的挑战[1-2].如能及时、准确地获得流动过程状态信息,对流动过程的产生、发展及转化进行分析,实现流动状态的表征和描述,可为进一步实现气液两相流流动过程的控制提供关键参数,对促进多相流问题研究的进一步发展,保障实际生产过程的安全稳定运行具有重要意义.

在气液两相流流动过程中,流速、相含率、压力或压力降等过程参数及相介质的分布反映了其流动过程特性.其中,相含率可以通过电学法[3-4]、射线法[5]、微波法[6]、快关阀法[7]等获取;流速可以通过差压法[8]、超声多普勒方法[9]等测量;截面阵列式电阻传感器可以获取流体不同位置处的相分布信息.在获取相关流动参数及相分布特性等基础上,采用小波变换[10]、经验模态分解、希尔伯特-黄变换[11]等时频分析方法,可以进一步揭示气液两相流流动过程状态的演化特性,为流动状态的判别和过程监控提供了有益的借鉴.

通过多传感器获取同类传感器的多模式信息以及异类传感器的多检测角度信息,可以更加完整地表征气液两相流的流动状态[1].根据不同传感器特点进行原始信号预处理,获得的多维时间序列蕴含了流动状态的变化特性.通过历史数据建立状态监测模型,获取实时数据进行状态监测,为实现气液两相流复杂流动状态的监测和分析提供了可行性方法.

基于数据驱动的多元统计监测方法如主成分分析(Principal component analysis,PCA)[12]、独立成分分析(Independent component analysis,ICA)[13]、偏最小二乘(Partial least squares,PLS)[14]等可以简化数据结构,解决多变量的耦合问题,在实际过程得到广泛应用[15-16].在气液两相流研究领域,PCA、ICA 等方法多用于流型识别和流体测量中传感器信号的特征提取.Shaban等[17]利用归一化后的差压信号建立PCA-ICA 模型,实现了垂直气液两相流流量测量;李凯锋等[18]利用PCA 提取电导率信息主成分,结合K-均值聚类算法实现了气水两相流流型辨识;Dong等[19]对电导环信号的多域特征采用等距特征映射、PCA 等方法得到三维向量实现了流型样本间结构的可视化.

PCA 作为实际工业过程最为成熟和常用的过程监测方法,可以从相互关联的高维数据中提取有用信息,实现数据的有效降维.Li等[20]将PCA 方法和支持向量数据描述结合起来用于制冷系统的故障诊断;Aouabdi等[21]结合多尺度熵和PCA,通过分析电机电流信号监测齿轮退化.PCA 假设过程变量是线性的且在时间上静态独立,仅提取了变量间的静态互相关,但实际的工业过程数据往往具有非线性、动态性等特点[22-23].Ku等[24]通过增加迟滞数据对数据矩阵进行动态拓展,提出了一种动态主成分分析方法(Dynamic PCA,DPCA).DPCA 提取变量在时间上的自相关和动态互相关关系,但仍未考虑到过程数据的非线性特性.基于以上问题,Choi等[25]提出了基于动态核主成分分析(Dynamic kernel PCA,DKPCA)的非线性动态过程监测方法,Zhang等[26]将特征向量选择与DKPCA 相结合提高了故障检测效率.

DKPCA 在工业领域中主要用于监测正常工况是否出现异常故障,时滞数据扩展适合描述过程的动态特性,核函数有利于捕获过程的非线性特性,其方法本身同时考虑了过程数据的动态性和非线性,符合气液两相流流动过程特点.在应用于气液两相流流动状态监测时,更应注重对流动状态的描述和分析[27-28].

为实现气液两相流动过程状态的监测,以实际流动过程中出现频率较高的泡状流、塞状流和弹状流3 种典型流动状态为研究对象,利用不同流动状态下多传感器测试数据所呈现的流动特性差异性,提出一种针对气液两相流状态监测的多模态动态核主成分分析方法(Multiple dynamic kernel PCA,MDKPCA).通过对多种流动状态采用历史数据建立多模态监测模型,实现对典型流动状态的判别及过渡状态的监测.

1 多模态动态核主成分分析监测

1.1 动态核主成分分析方法

DKPCA 将过程的动态特性和非线性特性融入模型,将输入空间X∈Rn×m拓展到具有l个时滞的增广矩阵X(l)以提取过程数据的动态特性[29],即

构造从输入空间到特征空间的非线性映射,以提取过程数据的非线性特性[30],在线性特征空间的协方差可以表示为

式中,Φ(·)表示将输入空间的输入向量映射到特征空间F的非线性映射函数.

通过求取特征空间的特征值可以对应得到特征空间的线性主元

式中,λ≥0,ν∈F{0},将λ值降序排列,最大的λ对应的ν为特征空间F中的第一个主元,最小的λ对应的ν为最后一个主元.

输入向量x(l)的主元成分t通过在特征空间F将Φ(·)映射到νk(k=1,2,···,p)而被提取

为避免直接计算非线性映射,引入点积形式的核函数k(x,y)=〈Φ(x),Φ(y)〉.

1.2 气液两相流多模态建模

气液两相流中不同的典型流动状态具有不同的流动特性,体现在过程参数的差异.因此,可采取多模态建模[31],即分别对每一种流动状态利用历史数据采用DKPCA 建模,根据实时数据进行状态监测.以实际流动过程中出现频率较高的泡状流、塞状流和弹状流3 种典型流动状态为研究对象,MDKPCA 建模步骤为:

4)将步骤2)中的径向基核函数k(x,y)代入式(4)中,分别提取3 种模态数据的非线性主成分:

5)计算T2和SPE 统计量及统计量控制限,统计量计算式为

式中,t为保留p个主元的主元得分向量,对角矩阵S=diag{λ1,···,λp}由建模数据协方差矩阵与保留主元数相对应的p个特征值所构成,表示保留所有主元的得分向量的平方和,表示保留p个主元的得分向量的平方和.

根据方差累计贡献率(Cumulative percent variance,CPV)选择主元数量.经检验,气液两相流的过程数据不符合高斯分布,因此选用核密度估计法(Kernel density estimation,KDE)计算T2及SPE 统计量控制限和SPEq,lim.二者的置信度水平均选择为α=0.99.

1.3 气液两相流状态监测

气液两相流流动过程作为动态系统,状态与状态间并非完全独立,当前的状态很大程度上依赖于过去的状态.单一时刻的过程参数在一定程度上并不能准确表征流动状态的特性:对于塞状流流动状态,一个完整的流动周期应包括由细小气泡聚集而成的长塞型气泡以及紧随其后的细小气泡,二者交替出现,间断地沿管壁顶部流动;对于弹状流流动状态,在一个流动周期内,含大量液体和细小气泡的液弹和含有大量气体的气弹交替出现.当采样点取到塞状流的细小气泡段及弹状流的液弹段,并不能表征其完整的流动过程特性.因此,在监测时应考虑到之前时刻的数据值,采用滑动窗(Moving window)技术,将当前采样时刻数据与前l个时刻的采样数据构成一个窗口宽度为l+1 的新数据样本.

监测步骤为:

1)获取新的实时数据xnew,按照所建立的不同状态模型的均值和方差进行标准化处理,得到经标准化处理后的xq,new.

2)将xq,new与前l个标准化处理后的数据纳入滑动窗,作为当前时刻采集到的数据样本new=[xq,l,xq,l-1,···,xq,new]T.

3)利用建模步骤3)中的Kq和1N,计算核向量kt∈R1×N并进行均值中心化处理:

4)利用建模步骤4)的式(6),令xq(l)=xq,new,提取非线性主成分tq,new.

如图1 所示为气液两相流流动状态MDKPCA建模与监测原理图.

图1 气液两相流MDKPCA 建模与监测原理图Fig.1 MDKPCA modeling and monitoring schematic diagram of gas-liquid two-phase flow

2 流动过程数据采集、预处理与建模

MDKPCA 方法用于气液两相流状态监测的可行性和有效性通过实验进行验证.实验产生不同工况下的多种流动状态,通过多传感器采集流动过程数据并进行预处理,进而进行MDKPCA 建模.

2.1 气液两相流实验

气液两相流实验于气液两相流水平环管实验装置上完成,实验装置如图2 所示.

图2 气液两相流水平环管实验装置Fig.2 Experimental apparatus of horizontal loop for gas-liquid two-phase flow

实验所用的管道为内径50 mm 的不锈钢管,流体分别为自来水(密度为998 kg/m3,动力粘度为1.01×10-3Pa · s)和干燥空气(密度为1.2 kg/m3),在管道入口处安装的混合器进行混合.流体进入混合器之前,使用标准单相流量计计量各相流量,由入口至出口距离为16.6 m.为获得充分发展的流动状态,多传感器测试管段安装在位于入口处下游约12 m.测试管段为透明有机玻璃管,可观察、并采用高速摄像机记录流动过程状态.流动过程的工况条件由测试管段上的压力和温度仪表进行监测和记录.

在气液两相流流动过程中,管道内各相流速、相分布、相含率及压力、温度等过程参数可表现流动状态的特性.为了多方面获取流动过程信息,选择用于数据获取的多传感器有:截面阵列式电阻、连续波超声多普勒、电容、电导环以及压力计.测试管段多传感器布置结构如图3 所示.各传感器所获数据表征的流体特性信息及相关的过程参数如表1所示.

图3 测试管段多传感器结构Fig.3 Structure of multiple sensors in test section

表1 气液两相流测量多传感器Table 1 Multiple sensors for gas-liquid two-phase flow measurement

实验包含7 组、共42 个实验条件测试点.每组水流量固定不变,通过调节空气流量由低至高变化,产生不同的流动状态,包括泡状流、塞状流、弹状流3 种典型流动状态和从泡状流到弹状流的过渡状态.水的流量依次设置为2,5,7,9,11,13和14 m3/h,标况(0 摄氏度,1 个大气压)下空气流量为1.5,2.5,5,7,19和36 m3/h.图4 所示为实验点在Mandhane 气液两相流流型图中的分布[32].图中实线为流型转换边界,不同图标点表示实验中出现的典型流动状态.

图4 气液两相流实验点分布Fig.4 Distribution of experimental points for gas-liquid two-phase flow

2.2 多传感器数据预处理

1)截面阵列式电阻

截面阵列式电阻传感器的工作原理为通过激励电极在测量区域建立敏感场,利用流体流动状态改变引起敏感场域内的电导率改变,导致测量边界电压发生变化.16 个电极沿着管壁等距离排列,采用相邻激励、相邻测量模式,采样频率为140 帧/s,每帧包含208 个边界电压测量值,对测量边界电压数据进行预处理[33]

式中,Vij0表示管道内充满水时第i个激励电极下的第j个边界电压测量值,Vij表示两相流动时第i个激励电极下的第j个边界电压测量值,VRi表示两相流动时第i个激励电极下13 个边界电压的平均值.

由于气液两相流流动结构复杂,截面阵列式电阻的测量数据反映了流体流动过程的二维介质分布信息,因此1 帧的数据并不能完整地表征流动过程的状态.根据测试系统中其他传感器的测试频率,需进行测试数据的时间配准.采用每次移动1 帧的滑动窗方法,对连续42帧VRi数据的幅值进行特征提取,并实现数据降维.提取的特征值包括均值、方差SD、偏度SK和峭度KI等用于表征流动气、液相介质的分布信息

式中,xi为时间序列VRi中的采样点,为xi的样本均值,n为时间序列VRi的采样点数量.

2)连续波超声多普勒

连续波超声多普勒传感器由一对压电陶瓷超声波换能器构成.根据管径及流速范围,选用直径9 mm 的压电晶片,安装倾斜角为55°.超声波发射器以1 MHz 的频率向被测流体中发射连续声波,超声波接收器接收被测流体中经离散相液滴调制的声波.根据超声多普勒效应,多普勒频移大小直接取决于离散相液滴的流速.连续波超声多普勒采样频率为50 kHz,采集的电压数据为时间序列,采用短时傅里叶变换(Short-time Fourier transform,STFT)将时域信号转变为频域信号,并求取多普勒速度[34]

式中,Sd(fd)为多普勒频移的能量谱,fd为多普勒频移组分,f0为超声换能器发出的超声波频率,c为超声波在流体中的纵波波速,θ为超声波声束方向与水平方向的夹角.

3)电导环

电导环传感器由轴向排列的6 个环形金属电极组成,采样频率为2 kHz,其中,E、G 为激励电极,M 为测量电极,相邻电极之间的距离为:E -M1和M4 -G 为50 mm,M1 -M2和M3 -M4 为20 mm,M2 -M3 为60 mm.采用电流激励电压测量的方式,通过激励电极对E-G,向管道内施加20 kHz的方波激励电流形成电场.连续相导电时(水为连续相),可以通过测量电极对之间的电压表征连续相(水)含率信息[34]

式中,Vw为满水测量值,Vmeas为实际测量值,Vn为归一化电压值.

4)电容

电容传感器由一对凹面金属极板C1和C2 构成,极板间混合流体的介电常数由各相的介电常数和相分布决定,因此测量电容与相含率存在对应关系.当连续相不导电时(气为连续相),可通过电容值计算离散相(水)含率.为消除管壁对检测结果带来的影响,采用相对电容变化量RCD表征离散相(水)含率信息[35]

式中,Cm为流体测量电容值,Ct为管道充满导电相时的电容值,CI为管道充满非导电相时的电容值,Vm为流体测量电压值,Vt为管道充满导电相时的电压值,VI为管道充满非导电相时的电压值.

压力计测量管道内压力,主要反映了流动过程的工况信息.由于环境温度变化幅度很小且对气液两相流流动状态影响较小,故未考虑温度参数.

2.3 过程参数建模

各传感器每个原始数据集的采样时间均为10 s,但由于采样频率不同,因此通过重采样对齐多传感器样本数据.经重采样后,每组的数据为Y(1 000×8)矩阵.以此为基础,对泡状流、塞状流、弹状流3 种典型流动状态进行MDKPCA 建模,并分别获取各典型流动状态下的T2 及SPE 统计量控制限.考虑到实际过程以及实验所获取的各流动状态的采样数量,采用的建模数据:泡状流为5 000 个采样点,塞状流为4 000 个采样点,弹状流为10 000 个采样点.由于气液两相流流动过程的状态需要一段时间内结构的变化进行描述,因此在分析过程状态的动态特性时,需采用一定时长的测试数据.根据实际气液两相流流动过程状态变化的特点,取采样的时滞l=29.即以时间序列长度为30 (对应时间长度为0.3 s)的采样点为一个样本点用于建模和测试.

3 气液两相流流动状态监测

3.1 典型流动状态判别

以气液两相流在实际流动过程中最为常见的泡状流、塞状流、弹状流等典型流动状态为研究对象,3 种典型流动状态如图5 所示.泡状流中,气泡非均匀分散于管道的中上部;塞状流中,部分细小气泡聚集成长气塞,气塞尾部伴随细小气泡在管道顶部流动;弹状流中,夹带气泡的液弹与较大速度的冲击性气弹快速交替.

图5 3 种典型流动状态Fig.5 Three typical flow states

分别采用多模态的主成分分析(Multiple PCA,MPCA)、多模态的动态主成分分析方法(Multiple dynamic PCA,MDPCA)和MDKPCA 方法进行状态判别,结果分别如图6～8 所示.

图6 中,由于采用的MPCA 方法仅提取变量间的静态互相关信息,组成数据样本的时滞l=0,样本数据为:1～1 000 采样点为泡状流状态,1 001～2 000 采样点为塞状流状态,2 001～3 000 采样点为弹状流状态.MPCA 状态判别结果表明:泡状流在3 种监测模型下除个别点外均未超限,即泡状流符合泡状流、塞状流及弹状流3 种流动状态模型的特性;此时,会将泡状流误判为塞状流或者弹状流.这一现象符合3 种流动状态的流动过程规律:塞状流和弹状流的完整周期包括细小气泡段.MPCA 方法只针对单一时刻的状态进行独立的线性分析,无法描述流体流动过程中前后状态间的关联性和状态变化的非线性,导致泡状流的误判.而塞状流和弹状流整体上判别效果好于泡状流,可以比较直观地区分,但在非塞状流和弹状流对应的监测模型下,仍有大量采样点在控制限下,存在较高的误判率.

图6 MPCA 监测模型典型状态判别Fig.6 MPCA-based state identification for typical states

图7 中,采用MDPCA 方法进行3 种流动状态的判别时,组成数据样本的时滞l=29,样本数据为:1～971 样本点为泡状流状态,972～1 942 样本点为塞状流状态,1 943～2 913 样本点为弹状流状态.MDPCA 状态判别结果表明:泡状流在3 种监测模型下均未超限,但相比于MPCA 方法,泡状流样本点在非泡状流监测模型下,均比模型对应的塞状流或弹状流数据整体更接近控制限.

图7 MDPCA 监测模型典型状态判别Fig.7 MDPCA-based state identification for typical states

图8 中,采用MDKPCA 方法对3 种流动状态进行判别时,组成数据样本的时滞l=29,样本数据为:1～971 样本点为泡状流状态,972～1 942 样本点为塞状流状态,1 943～2 913 样本点为弹状流状态.MDKPCA 状态判别结果表明:泡状流在塞状流监测模型下,SPE 统计量完全超限;在弹状流监测模型下,虽仍有部分样本点低于控制限下,但相比于弹状流数据,T2和SPE 统计量整体更接近控制限.塞状流和弹状流数据的误判率进一步降低.相比于前两种方法,3 种典型流动状态的判别效果有了很大的提升.

图8 MDKPCA 监测模型典型状态判别Fig.8 MDKPCA-based state identification for typical states

为进一步讨论不同典型流动状态对MPCA、MDPCA和MDKPCA 等3 种方法所建立监测模型判别效果的影响,对T2和SPE 统计量进一步分析.定义模型对各典型流动状态的符合度为:观测的典型流动状态被判别为所用模型对应流动状态的比率.定义判别正确率β为:3 种监测模型对观测的典型流动状态均正确判别的比率.

分别采用MPCA、MDPCA和MDKPCA 方法得到各监测模型对所观测到的3 种典型流动状态的判别结果如表2 所示.

表2 所示数据中,MDPCA 方法对塞状流和弹状流的判别效果明显优于MPCA 方法,将塞状流和弹状流在非对应监测模型下误判为其他流动状态的比例明显降低,判别正确率明显提高.表明了MDPCA 方法通过提取状态变化的动态自相关和互相关信息,把流体流动过程近似为一种线性的动态过程处理后,对过程数据中包含的动态特性敏感性增强.但MDPCA 方法难以将泡状流和其他两种流动状态的流动特性区分,观测的泡状流状态在塞状流和弹状流监测模型下仍存在100%和76.93%的误判率.分析原因,主要是MDPCA 方法未考虑流动过程的非线性变化及测试数据中耦合的这种非线性变化特点.

表2 观测流动状态下3 种方法所建立监测模型判别结果(%)Table 2 Identification results of monitoring models in three methods under observation of flow states (%)

MDKPCA 方法在MDPCA 方法的基础上,进一步将测试数据中耦合的非线性变化特点,通过核方法映射到高维空间进行线性化处理,提取塞状流和弹状流在流动过程中液滴和气泡不断聚并、破碎的动力学演化过程的细节特征,用于区分不具有上述演化过程特征的泡状流状态.因此,MDKPCA方法与MPCA和MDPCA 方法相比,建立的监测模型对模型非对应的典型流动状态误判率明显降低,3 种典型流动状态判别正确率明显提高.MDKPCA 方法可以更好地描述气液两相流不同流动状态的流动过程特性,进行状态监测时效果更佳.

3.2 过渡状态监测

在气液两相流流动过程中,从一种稳定的流动状态演化到另一种稳定的流动状态之间存在过渡状态,以稳定的泡状流过渡到稳定的弹状流这一完整的流动过程为研究对象,由高速摄像机记录观测到的过程状态如图9 所示.分别采用MPCA、MDPCA、MDKPCA 方法进行状态监测,监测结果分别如图10～12 所示.

图9 泡状流过渡到弹状流的过程状态Fig.9 Transition from bubble flow to slug flow

图10 中,采用MPCA 方法进行状态监测,1～2 000 个采样点为稳定的泡状流,2 001～4 000 个采样点为过渡状态,4 001～6 000 个采样点为稳定的弹状流.MPCA 监测结果表明:在泡状流监测模型下,只有T2统计量可以体现从稳定的泡状流过渡到稳定的弹状流状态转变的过程;在弹状流监测模型下,两个统计量值在整个过程中均位于对应的控制限下,完全无法区分流动状态和实现状态监测.

图10 MPCA 监测模型过渡状态监测Fig.10 MPCA-based monitoring for transitions

图11 中,采用MDPCA 方法进行状态监测,1～1 942 个样本点为稳定的泡状流,1 943～3 884 个样本点为过渡状态,3 885～5 826 个样本点为稳定的弹状流.MDPCA 监测结果表明:在泡状流监测模型下,两个统计量均可以体现从稳定的泡状流过渡到稳定的弹状流状态转变的过程;在弹状流监测模型下,虽然两个统计量都可以反映整个过渡过程的变化趋势,但不能准确区分泡状流稳定状态.

图11 MDPCA 监测模型过渡状态监测Fig.11 MDPCA-based monitoring for transitions

图12 中,采用MDKPCA 方法进行状态监测,1～1 942 个样本点为稳定的泡状流,1 943～3 884 个样本点为过渡状态,3 885～5 826 个样本点为稳定的弹状流.MDKPCA 监测结果表明:在泡状流监测模型下,两个统计量均直观地体现了整个流动过程的状态变化;在弹状流监测模型下,SPE 统计量可以准确区分泡状流,实现对稳定和过渡状态的准确监测.当处于过渡状态时,统计量值的移动和波动反映了在泡状流过渡到弹状流过程中,气相流量不断增大,气相和液相呈现出剧烈、复杂的相互作用,揭示了气液流动过程的转化并非是瞬态的,而是过去状态的不断累积和演化.

图12 MDKPCA 监测模型过渡状态监测Fig.12 MDKPCA-based monitoring for transitions

4 结束语

针对水平气液两相流流动状态监测问题,从其状态的变化具有时变性、非线性、随机性等特点出发,在多传感器获取流动状态测试数据的基础上,通过采用动态自相关及互相关方法提取过程的动态特性,核方法提取非线性特性,结合主成分分析的MDKPCA 方法,建立不同典型流动状态的监测模型,并以T2和SPE 统计量作为监测指标,采用滑动窗技术保障气液两相流流动过程信息的连续性,结合所建立的多模态模型对不同的典型流动状态进行判别,并进一步实现流动过渡状态监测.

与采用MPCA和MDPCA 方法相比,MDKPCA方法在前者的基础上,通过进一步提取流动过程的动态和非线性特征,能够更全面地描述不同流动状态的过程特性.气液两相流实验装置测试数据的处理结果验证了MDKPCA 方法对典型流动状态判别的准确性和对过渡状态监测的有效性.