基于独立成分分析和互相关技术的水平管弹状流参数测量方法.pdf

本发明提供基于独立成分分析和互相关技术的水平管弹状流参数测量方法，它是在由结构相同的双截面电学传感器阵列、电子测量电路和成像计算机构成的电学层析成像系统中，引入独立成分分析和互相关技术实现水平管两相流弹状流特征参数测量方法，即，使用独立成分分析方法分析原始测量数据，发现峭度最大的独立成分与弹状流中液弹的变化信息具有很好的对应关系，进一步利用互相关技术，对上、下游测量截面获得的峭度最大的独立成分消除循环平稳特性后进行互相关计算，获得弹状流中液弹的平均速度，利用液弹持续的时间乘以液弹的平均速度获得平均的液弹长度。这种方法利用电学层析成像系统所获得的测量值无需结合其他测量数据获得水平管弹状流的速度和长度，节省了测量的成本，拓宽了电学层析成像系统在两相流测量中的应用范围，有利于精确测量水平管弹状流参数。

1、  一种基于独立成分分析和互相关技术的水平管弹状流参数测量方法，该测量方法是在由结构相同的双截面电学传感器阵列、电子测量电路和成像计算机构成的电学层析成像系统中，对水平管两相流弹状流测量，该测量方法包含有以下步骤：
(1)双截面电学层析成像系统中，若每一截面有N个测量电极，分别将由双截面电学层析成像系统获得的水平管道弹状流的原始完全测量数据转换成N个电极上的时间序列矩阵T1_N×M、T2_N×M，其中，T1_N×M为上游测量截面的时间序列矩阵，T2_N×M为下游测量截面的时间序列矩阵，下标N为每一截面有N个电极，下标M为每个电极上的采样数为M，记TI_N×M，I＝1，2表示上述时间序列矩阵，I＝1时，即为上述T1_N×M，I＝2时，即为上述T2_N×M；
(2)分别对所述时间序列矩阵T1_N×M、T2_N×M进行预处理，包括中心化和白化，获得白化后的矩阵Z_i，下标i＝1，2，下标i＝1时，即Z₁，表示由上游测量截面时间序列矩阵T1_N×M经预处理获得的白化后的矩阵，下标i＝2时，即Z₂，表示由下游测量截面时间序列矩阵T2_N×M经预处理获得的白化后的矩阵；
(3)将所述白化后的矩阵Z_i进行变换，使输出矩阵Y_i＝W_iZ_i的各行向量之间相互独立，即求解混合矩阵W_i，进而得到独立成分集合Y_i，Y_i为N×M的矩阵，Y_i中每一行对应于一个独立成分y_k，下标i＝1，2分别对应上、下游测量截面的变量，下标i＝1时，对应上游测量截面的变量，下标i＝2时，对应下游测量截面的变量，Y_i代表独立成分集合，下标i＝1时，即Y₁，代表上游测量截面的独立成分集合，下标i＝2时，即Y₂，代表下游测量截面的独立成分集合，W_i代表解混合矩阵，下标i＝1时，即W₁，代表上游测量截面的解混合矩阵，下标i＝2时，即W₂，代表下游测量截面的解混合矩阵；
(4)分别计算所述独立成分集合Y₁、Y₂中各独立成分y_k的峭度，分别提取Y₁、Y₂中峭度最大的独立成分；
(5)分别消除双截面峭度最大的独立成分的循环平稳特性，得到特征向量C1、C2；
(6)对所述特征向量C1、C2进行互相关计算，获得相关速度即是液弹的平均速度；
(7)利用所述液弹的平均速度乘以液弹持续的时间，即可获得平均的液弹长度。

2、  根据权利要求1所述的基于独立成分分析和互相关技术的水平管弹状流参数测量方法，其特征是：所述步骤(1)中时间序列矩阵，是指对每个电极按照采样顺序得到行向量，测量电极数为N个，时间序列矩阵的行数则为N。

3、  根据权利要求1所述的基于独立成分分析和互相关技术的水平管弹状流参数测量方法，其特征是：所述步骤(2)预处理中的中心化，是指以时间序列矩阵TI_N×M的某一行向量TI_k×M为对象，计算该行元素的平均值E[TI_k×M]，将行向量TI_k×M的元素都减去该行元素的平均值E[TI_k×M]，获得向量x_k＝TI_k×M-E[TI_k×M]，下标k从1依次取到N，获得矩阵X_i＝[x₁，x₂，…，x_N]′，下标i＝1，2，分别对应上、下游测量截面对应的变量，即X₁对应上游测量截面的中心化后的矩阵，X₂对应下游测量截面的中心化后的矩阵；
所述述步骤(2)预处理中的白化，是去除信号x_k间的相关性，并得到白化后的矩阵Z_i，使得Z_i＝V_iX_i，其中，V_i为白化矩阵，X_i为中心化后的矩阵，Vi=Di-1/2EiT,]]>D_i为矩阵X_i的协方差矩阵E{X_iX_i^T}的特征值组成的对角阵，E_i为E{X_iX_i^T}的特征向量组成的正交阵，下标i＝1，2，分别对应上、下游测量截面的变量，i＝1时对应上游测量截面的变量，i＝2时对应下游测量截面的变量。

4、  根据权利要求1所述的基于独立成分分析和互相关技术的水平管弹状流参数测量方法，其特征是：所述步骤(3)中求解混合矩阵W_i，具体如下：
(1)初始化解混矩阵W_i(随机数)；
(2)计算Y_i＝W_iZ_i和vi=11+exp(-Yi);]]>
(3)计算权值增量ΔWi=μ∂H(vi)∂Wi=μ[I+(1-2vi)YiT]Wi;]]>
(4)更新权值W_i＝W_i+ΔW_i；
(5)计算性能函数Ei=Σi=1n(Σj=1n|pkj|maxl|pkl|-1)+Σj=1n(Σi=1n|pkj|maxl|plj|-1),]]>其中：
P_i＝(p_kj)＝W_iA_i，Ai=(WiTVi)-1,]]>判断是否达到收敛条件E_i＜ε，是，结束；否则，回到所述步骤(2)；下标i＝1，2，分别对应上、下游测量截面的变量，i＝1时对应上游测量截面的变量，i＝2时对应下游测量截面的变量；
所述步骤(3)中独立成分集合，是指：由于Z_i＝V_iX_i，Y_i＝W_iZ_i，所以Xi=(WiTVi)-1Yi,]]>令(WiTVi)-1=Ai,]]>Y_i＝S_i，故有X_i＝A_iS_i，矩阵S_i体现了原始数据的稀疏表达，矩阵S_i每一行称作一个独立成分，Y_i为S_i的估计，下标i＝1，2，分别对应上、下游测量截面的变量，i＝1时对应上游测量截面的变量，i＝2时对应下游测量截面的变量。

5、  根据权利要求1所述的基于独立成分分析和互相关技术的水平管弹状流参数测量方法，其特征是：所述步骤(4)求取各独立成分y_k的峭度，具体如下：
β=1T∫0T[yk-uyk]4dt,]]>K=βσyk4,]]>其中为信号y_k的标准差，uyk=0]]>为独立成分y_k的均值。

6、  根据权利要求1所述的基于独立成分分析和互相关技术的水平管弹状流参数测量方法，其特征是：所述步骤(5)中消除独立成分的循环平稳特性，是以T为周期的相关函数在一个周期T内作时间平均，在相邻激励模式下，就是对前述步骤(4)获得的独立成分每(N-1)个数据做一个平均，获得新的特征向量。

7、  根据权利要求1所述的基于独立成分分析和互相关技术的水平管弹状流参数测量方法，其特征是：所述步骤(6)中互相关计算，包括但不限于用基于互相关方法的时间延迟估计求滞后时间，S₁₂(f)是C1、C2的互谱，S₁₁(f)，S₂₂(f)分别是C1、C2的自谱，则相干函数的平方为：
C12(f)=|S12(f)|2S11(f)S22(f)]]>
极大似然窗为：
W(f)=1|S12(f)|C12(f)(1-C12(f))]]>
于是加窗的互相关函数R₁₂(τ)就等于W(f)S₁₂(f)的傅立叶反变换，R₁₂(τ)在τ＝D处存在一个峰值，延迟D的初始估计值d是R₁₂(τ)的峰值位置，所对应的延迟时间为τ=Dfc,]]>f_c为采样率；
所述步骤(6)中相关速度，为L为两电学传感器阵列截面的轴向间距。

8、  根据权利要求1所述的基于独立成分分析和互相关技术的水平管弹状流参数测量方法，其特征是：所述步骤(7)液弹持续时间，由液弹脉冲持续的数据点数/采集率获得，液弹的平均速度乘以液弹持续的时间，即为平均的液弹长度。

基于独立成分分析和互相关技术的水平管弹状流参数测量方法
技术领域
本发明属于两相流测量技术领域，涉及过程流动成像研究领域的信息技术，特别是一种基于独立成分分析和互相关技术的水平管弹状流参数测量方法。
背景技术
弹状流(slug flow)是多相流(尤其是油气混输管道)中一种最常见的流型，由于其流动的间歇性，引起管道中持液率和压力的急剧波动，并使运行在该流型下的管道相应承受间歇性应力冲击。同时，离开管道末端的大液塞会引起下游设备中的液位剧烈波动。为保证管线和下游处理设备的最优设计和安全操作，必须能够掌握、测量和预报弹状流特征参数的大小与变化规律。弹状流的特征参数强烈依赖于管道的几何尺寸和流体的物性参数，由于实验室试验条件和实际工业管道几何与运行条件的差别，实验室模拟试验的结果应经过现场试验数据的检验和修正才有应用价值。因此，如何从现场工业管道上获得有价值的试验数据，成为弹状流流动特性研究的关键。
电学层析成像技术是基于电特性敏感机理的过程层析成像技术，它的物理基础是基于不同的媒质具有不同的电特性(电导率/介电系数/复导纳)，判断出敏感场内物体的电特性分布便可推知该场中媒质的分布情况。它运用传感器阵列形成旋转的空间敏感场，对包含多种不同介质的容器或管道从不同的观测视角进行扫描，当场内介质的电特性分布发生变化时，场内电流密度、电势分布也随之变化，从而导致场域边界上的测量数据发生变化，边界测量数据的变化信息对应电特性的变化信息，运用相应的成像算法，便可重建出实际对象的物质分布，并可对物场分布实现可视化。
目前电学层析成像技术在两相流测量的应用，主要可以概括为以下几个方面：
1、提供被测流体在流体管道某一横截面的实时图像，用于识别与判断两相流的流型及确定相间界面；
2、用于确定两相流体各相组分在容器或管道中的局部浓度分布，进一步处理分析可得到各相组分的总浓度；
3、确定离散相的颗粒尺寸和运动轨迹等微尺度参数；
4、将电学层析成像技术与相关技术结合可以获得两相流的速度信息，通过图像重建可以获得截面分相含率信息，二者结合可实现对两相流体的总流量、各分相流量的测量及流体在管道横截面上速度分布的实时测量。
其中与相关技术的结合由于可以获得两相流的速度信息，拓宽了电学层析成像的应用范围。
目前文献中，通过像素相关的方法可以获得测量截面速度场的分布，如Wang M等人2006年发表于《粒子与粒子系统表征》(Particle & Particle Systems Characterization)第23卷，第321-329页，题为《应用电阻层析成像可视化漩涡流中气泡速度分布》(Visualisation of bubbly velocity distribution in a swirling flow using electrical resistancetomography)的文章，如G.Vilar等人2008年发表于《化学工程学报》(Chemical EngineeringJournal)，题为《应用电阻抗层析成像在线分析振荡折流板反应器中分布结构》(On lineanalysis of structure of dispersions in an oscillatory baffled reactor using electrical impedancetomography)等；但是由于无法将截面的每个象素与实际的多相流组分一一对应，所以无法准确获得相关速度的物理意义，也无法获得弹状流中液弹的平均速度和液弹长度。
目前文献中提及利用电学层析成像技术对弹状流进行测量，需要结合更多的其他数据，比如视频图像，如K.L.Ostrowski等2000年发表于《化学工程学报》(ChemicalEngineering Journal)的第77卷，第43-50页，题为《电容层析成像在水平管气动输送器流型在线、离线分析中的应用》(Application of capacitance electrical tomography for on-lineand off-line analysis of flow pattern in horizontal pipeline of pneumatic conveyer)。但是具体的应用电学层析成像技术测量弹状流中液弹的速度和长度，尤其是只利用电学层析成像系统所获得的测量值不结合其他测量数据的方法尚未有报道。
发明内容
本发明的目的在于提出一种基于独立成分分析和互相关技术的水平管弹状流参数测量方法，有利于精确测量水平管弹状流参数，节省了测量成本。
为实现上述目的，本发明采用的技术方案是提供一种基于独立成分分析和互相关技术的水平管弹状流参数测量方法，该测量方法是在由结构相同的双截面电学传感器阵列、电子测量电路和成像计算机构成的电学层析成像系统中，对水平管两相流弹状流测量，该测量方法包含有以下步骤：
(1)双截面电学层析成像系统中，若每一截面有N个测量电极，分别将由双截面电学层析成像系统获得的水平管道弹状流的原始完全测量数据转换成N个电极上的时间序列矩阵T1_N×M、T2_N×M，其中，T1_N×M为上游测量截面的时间序列矩阵，T2_N×M为下游测量截面的时间序列矩阵，下标N为每一截面有N个电极，下标M为每个电极上的采样数为M，记TI_N×M，I＝1，2表示上述时间序列矩阵，I＝1时，即为上述T1_N×M，I＝2时，即为上述T2_N×M；
(2)分别对所述时间序列矩阵T1_N×M、T2_N×M进行预处理，包括中心化和白化，获得白化后的矩阵Z_i，下标i＝1，2，下标i＝1时，即Z₁，表示由上游测量截面时间序列矩阵T1_N×M经预处理获得的白化后的矩阵，下标i＝2时，即Z₂，表示由下游测量截面时间序列矩阵T2_N×M经预处理获得的白化后的矩阵；
(3)将所述白化后的矩阵Z_i进行变换，使输出矩阵Y_i＝W_iZ_i的各行向量之间相互独立，即求解混合矩阵W_i，进而得到独立成分集合Y_i，Y_i为N×M的矩阵，Y_i中每一行对应于一个独立成分y_k，下标i＝1，2分别对应上、下游测量截面的变量，下标i＝1时，对应上游测量截面的变量，下标i＝2时，对应下游测量截面的变量，Y_i代表独立成分集合，下标i＝1时，即Y₁，代表上游测量截面的独立成分集合，下标i＝2时，即Y₂，代表下游测量截面的独立成分集合，W_i代表解混合矩阵，下标i＝1时，即W₁，代表上游测量截面的解混合矩阵，下标i＝2时，即W₂，代表下游测量截面的解混合矩阵；
(4)分别计算所述独立成分集合Y₁、Y₂中各独立成分y_k的峭度，分别提取Y₁、Y₂中峭度最大的独立成分；
(5)分别消除双截面峭度最大的独立成分的循环平稳特性，得到特征向量C1、C2；
(6)对所述特征向量C1、C2进行互相关计算，获得相关速度即是液弹的平均速度；
(7)利用所述液弹的平均速度乘以液弹持续的时间，即可获得平均的液弹长度。
本发明的有益效果是该方法拓宽了电学层析成像在两相流测量领域的应用范围，无需结合其他途经获得的更多的测量数据，只利用电学层析成像自身的数据即可获得水平管弹状流液弹的速度和液弹的长度这两个重要的弹状流特征参数，节省了测量成本。该方法通过对电学阵列传感器提供的原始测量数据进行独立成分分析，利用峭度最大的独立成分与弹状流中液弹的变化信息具有很好的对应关系，结合互相关技术，将上、下游测量截面获得原始测量数据进行独立成分分析，提取峭度最大的独立成分，消除循环平稳特性后进行相关计算，可以获得弹状流中液弹的平均速度，平均的液弹长度可以通过液弹持续的时间乘以液弹的平均速度获得。
附图说明
图1为本发明的基于独立成分分析和互相关技术的水平管弹状流参数测量方法的流程框图；
图2为本发明的电极、电阻层析成像系统圆形管道传感器截面分布示意图；
图3为本发明的电极、电阻层析成像系统圆形管道传感器轴向分布示意图；
图4为本发明的由上游测量截面利用独立成分分析方法估计获得的八个独立成分；
图5为本发明的上游测量截面和下游截面获取的液弹信息；
图6(a)为本发明的消除两截面峭度最大的独立成分的循环平稳特性得到特征向量C1、C2；
图6(b)为本发明的局部区域内的放大对比；
图7(a)、(b)、(c)、(d)、(e)分别为估算出的图6(b)中A、B、C、D、E处的液弹平均速度。
图中：
1、有机玻璃管    2、电极
具体实施方式
结合附图和实施例对本发明的基于独立成分分析和互相关技术的水平管弹状流参数测量方法加以说明。
本发明的基于独立成分分析和互相关技术的水平管弹状流参数测量方法引入独立成分分析和互相关技术实现水平管两相流弹状流特征参数测量方法，使用独立成分分析方法分析原始测量数据，利用峭度最大的独立成分与弹状流中液弹的变化信息具有很好的对应关系，进一步利用互相关技术，对上、下游测量截面获得的峭度最大的独立成分消除循环平稳特性后进行互相关计算，获得弹状流中液弹的平均速度，利用液弹持续的时间乘以液弹的平均速度获得平均的液弹长度而进行检测。
图2所示为一电学层析成像之一电阻层析成像系统的圆形单截面阵列电阻传感器，八个电极均匀分布在有机玻璃管道内壁。图3所示，在直径为D的有机玻璃管的内壁上安装有两个具有相同电极分布结构的阵列电阻传感器，间隔距离为L。
如图1所示，本发明的基于独立成分分析和互相关技术的水平管弹状流参数测量方法实施例包括如下具体步骤：
1、双截面电学层析成像系统中，每一截面有八个测量电极，分别将由双截面电学层析成像系统获得的水平管道弹状流的原始完全测量数据转换成八个电极上的时间序列矩阵T1_8×M、T2_8×M，其中，T1_8×M为上游测量截面的时间序列矩阵，T2_8×M为下游测量截面的时间序列矩阵，下标M为每个电极上的采样数为M；
所述时间序列矩阵，是指对每个电极按照采样顺序得到的行向量，测量电极数为八个，时间序列矩阵的行数为八行。
2、分别对所述时间序列矩阵T1_8×M、T2_8×M进行预处理，包括中心化和白化，获得白化后的矩阵Z₁、Z₂，Z₁为由上游测量截面时间序列矩阵T1_8×M经预处理获得的白化后的矩阵，Z₂为由下游测量截面时间序列矩阵T2_8×M经预处理获得的白化后的矩阵；
所述中心化，是指以上游测量截面时间序列矩阵T1_8×M的某一行向量T1_k×M为对象，计算该行元素的平均值E[T1_k×M]，将行向量T1_k×M的元素都减去该行元素的平均值E[T1_k×M]，获得向量x_k＝T1_k×M-E[T1_k×M]，下标k从1依次取到8，获得上游测量截面的中心化后的矩阵X₁＝[x₁，x₂，...，x₈]′，按同样方法，以下游测量截面时间序列矩阵T2_8×M的某一行向量T2_k×M为对象，计算该行元素的平均值E[T2_k×M]，将行向量T2_k×M的元素都减去该行元素的平均值E[T2_k×M]，获得向量x_k＝T2_k×M-E[T2_k×M]，下标k从1依次取到8，获得下游测量截面的中心化后的矩阵X₂＝[x₁，x₂，...，x₈]′；
所述白化，是去除信号x_k间的相关性，并得到白化后的矩阵Z₁、Z₂，使得Z₁＝V₁X₁、Z₂＝V₂X₂，其中，V₁为上游测量截面的白化矩阵，X₁为上游测量截面中心化后的矩阵，V₁＝D₁^-1/2E₁^T，D₁为矩阵X₁的协方差矩阵E{X₁X₁^T}的特征值组成的对角阵，E₁为E{X₁X₁^T}的特征向量组成的正交阵，V₂为下游测量截面的白化矩阵，X₂为下游测量截面中心化后的矩阵，V₂＝D₂^-1/2E₂^T，D₂为矩阵X₂的协方差矩阵E{X₂X₂^T}的特征值组成的对角阵，E₂为E{X₂X₂^T}的特征向量组成的正交阵。
3.将所述白化后的矩阵Z_i进行变换，使输出矩阵Y_i＝W_iZ_i的各行向量之间相互独立，即求解混合矩阵W_i，进而得到独立成分集合Y_i，Y_i为8×M的矩阵，Y_i中每一行对应于一个独立成分y_k，下标i＝1，2分别对应上、下游测量截面的变量；
所述求解混合矩阵W_i，具体如下：
(1)初始化解混矩阵W_i，该矩阵为8×8的矩阵，且矩阵里的各元素为任一随机实数；
(2)计算Y_i＝W_iZ_i和vi=11+exp(-Yi);]]>
(3)计算权值增量ΔWi=μ∂H(vi)∂Wi=μ[I+(1-2vi)YiT]Wi;]]>
(4)更新权值W_i＝W_i+ΔW_i；
(5)计算性能函数Ei=Σi=1n(Σj=1n|pkj|maxl|pkl|-1)+Σj=1n(Σi=1n|pkj|maxl|plj|-1),]]>其中：
P_i＝(p_kj)＝W_iA_i，Ai=(WiTVi)-1,]]>判断是否达到收敛条件E_i＜ε，是，结束；否则，回到所述步骤(2)；下标i＝1，2，分别对应上、下游测量截面的变量，i＝1时对应上游测量截面的变量，i＝2时对应下游测量截面的变量；
所述独立成分集合，是指：由于Z_i＝V_iX_i，Y_i＝W_iZ_i，所以Xi=(WiTVi)-1Yi,]]>令(WiTVi)-1=Ai,]]>Y_i＝S_i，故存X_i＝A_iS_i，矩阵S_i体现了原始数据的稀疏表达，矩阵S_i每一行称作一个独立成分，Y_i为S_i的估计，下标i＝1，2，分别对应上、下游测量截面的变量，i＝1时对应上游测量截面的变量，i＝2时对应下游测量截面的变量。
图4为由第一个测量截面利用上述独立成分分析方法估计获得的八个独立成分，横坐标为时间，单位为秒。ICi为第i个独立成分(i＝1，...，8)。
4、分别计算所述独立成分集合Y₁、Y₂中各独立成分y_k的峭度，分别提取Y₁、Y₂中峭度最大的独立成分；
所述求取各独立成分y_k的峭度，具体如下：
β=1T∫0T[yk-uyk]4dt,]]>K=βσyk4,]]>
其中：为信号y_k的标准差，uyk=0]]>为独立成分y_k的均值。表1为图4中各独立分量的峭度。峭度最大的独立成分代表了弹状流中液弹的变化。图5为上游测量截面和下游截面获取的液弹信息，即两截面峭度最大的独立成分。
表1