CN106872530A - 两相流含水率自适应串联估计方法 - Google Patents

Abstract

本发明涉及一种两相流含水率自适应串联估计方法，利用上、下游不同传感器的测量数据进行估计，包括如下步骤：对上游传感器和下游传感器测量数据进行时间配准；状态估计算初始化；对上游传感器含率测量值进行线性卡尔曼状态估计和更新，获得对状态向量的预测和协方差矩阵的预测；3)利用下游传感器测量数据，计算两相流过程状态转换矩阵：对下游含率测量值进行基于无迹卡尔曼滤波非线性状态估计。

Description

两相流含水率自适应串联估计方法

技术领域

本发明属于流体测量技术领域，利用自适应无迹卡尔曼算法提升相含率传感器对两相流含水率估计准确度。

技术背景

两相流广泛存在于日常生活与工业生产过程中，如食品加工，生物工程，化工产业，冶金工业和石油产业等行业。与一般单相流相比，两相流的流动状态更为复杂，建模更加困难，因此两相流的在线过程参数检测一直是学术界和工业界关注的难点和重点。特别在石油工业中，油水两相流含水率对油井监测、生产预测和开采计划制定具有重要意义。

目前两相流含水率的测量手段包括快关阀测量方法、射线方法、微波方法、Wire-mesh方法、超声方法和电阻抗方法等。其中，电学方法具有测量速度快，敏感性强，传感器结构简单，易于安装，成本低廉等优点。两相流是一种随机的动态流动过程，流动状态复杂多变，传感器的测量信号有随机误差，且只能反映系统的外部特性。为获得多相流含水率的准确信息，需要根据实测数据对流动状态进行估计，以状态估计算法提升相含率传感器的测量精度。状态估计对认知和分析测量系统具有重要的应用价值，常用的估计理论包括最小二乘估计、线性最小方差估计和卡尔曼估计等。不同的估计方法性能指标不同，使用的信息不同，其适用范围也不同。最小二乘估计通常应用于随机向量和常值向量的估计，算法结构简单但精度有限；线性最小方差估计以均方误差最小作为优化参数，适用于平稳过程，其测量信息越丰富，估计精度越高，计算量越大，而且会随着估计过程的累积，造成计算负担较重，不适用于复杂的长时间的估计过程。卡尔曼估计采用递推的形式，通过一个状态模型和一个观测模型的迭代，利用上一时刻的估计值和这一时刻的观测值，对这一刻的状态进行估计。而基于不同结构的卡尔曼算法，卡尔曼估计又衍生出基本的线性卡尔曼估计，扩展卡尔曼估计，基于无迹卡尔曼滤波的非线性卡尔曼估计等估计方式，分别适用于不同对象。

发明内容

本发明提供一种两相流含水率自适应估计方法，该方法协同利用管道上、下游相含率传感器测量数据，上、下游相含率传感器可为同类传感器或数据采样特征匹配的异类传感器。对上游传感器测量数据进行线性卡尔曼滤波，计算其实时的状态转移矩阵，将该矩阵作为下游传感器状态转移矩阵的先验信息，代入下游传感器测试数据的无迹卡尔曼估计迭代过程中，实现含水率估计的自适应状态转换，同时调整无迹卡尔曼算法中Sigma点集结构，提高含水率的估计精度。本发明的技术方法如下：

一种两相流含水率自适应串联估计方法，利用上、下游不同传感器的测量数据进行估计，包括如下步骤：

1)对上游传感器和下游传感器测量数据进行时间配准：

获取上游传感器和下游传感器所获相含率时间序列z₁(j)和z₂(j)，j＝1,2,3,…N为测量时间内的数据点，下标1和2分别代表上游传感器和下游传感器，计算渡越时间，利用渡越时间结合数据采样频率，对上游传感器和下游传感器测量数据进行时间配准；

2)状态估计算初始化：

其中，是上游状态向量的估计初值，是上游协方差的迭代初值，z_1,k是上游传感器在k时刻的测量值，且 是下游状态向量的估计初值；

3)对上游传感器含率测量值进行线性卡尔曼状态估计和更新，获得对状态向量的预测和协方差矩阵的预测；

在估计过程中量测矩阵H取值均为1，状态向量x₁代表含水率，k时刻的系统模型为：

x_1,k+1＝F_1,kx_1,k+v_1,k (4)

z_1,k＝H_1,kx_1,k+w_1,k (5)

其中，F_1,k是k时刻的状态转移矩阵，x_1,k+1和x_1,k分别是上游在k+1时刻和k时刻的状态向量，H_1,k是k时刻的量测矩阵，z_1,k是量测向量，v_1,k和w_1,k分别表示系统在k时刻的过程噪声和量测噪声；则状态向量的预测的预测和协方差矩阵的预测为：

根据上游量测信息z_1,k修正估计值：

其中，K_1,k是上游传感器系统在k时刻的卡尔曼增益：

其中，H_1,k是k时刻上游传感器量测矩阵，是上游传感器在k时刻对状态的预测值，是上游传感器k时刻对状态x_1,k的最小均方误差估计，是k时刻估计的协方差，R_1,k为上游滤波过程在k时刻系统量测噪声的协方差，可取值为0.01至0.2；Q_1,k-1为上游滤波过程在k-1时刻系统过程噪声的协方差，可取值为0.01至0.09，I为单位矩阵，其维数等同于状态向量的维数；

3)利用下游传感器测量数据，计算两相流过程状态转换矩阵：

将传感器上游测量单元确定的状态转移过程加以记录，用于时间配准后下游含水率状态估计中：

其中，F_2,k为下游测量单元的状态转移矩阵，k＝n₀,n₀+1,…,N；

4)对下游含率测量值进行基于无迹卡尔曼滤波非线性状态估计：

构造Sigma点集，其基本构造方式如式(12)所示：

其中，L是状态向量的维度，是下游估计过程k-1时刻的状态估计，p_2,k-1是下游估计过程k-1时刻的全局协方差，γ与λ可由下式计算：

其中，λ为规模因子，α是决定Sigma点集中点的分散程度扩展因子，取值为0到1之间；κ是缩放系数，设为0，计算状态预测权系数W_i ^(m)和协方差预测权系数W_i ^(c)：

其中，β先验分布因子，通常设置为2，i为矩阵在Sigma点集中对应的元素序号，下游状态估计过程的状态预测方式如下：

估计过程中的协方差矩阵分为全局协方差矩阵和互协方差矩阵，更新方式分别如下：

其中，p_{2,k k-1}是对全局协方差的预测，Q_2,k-1是下游过程噪声的协方差，H_2,k是下游估计过程在k时刻的量测矩阵；

其中，是下游传感器测量值在k时刻的协方差，是k时刻互协方差矩阵，R_2,k为下游估计过程中的系统量测噪声的协方差，取值为0.01至0.09；K_2,k是k时刻下游估计过程的卡尔曼增益，按照式(19)方式进行更新；

基于无迹卡尔曼的非线性状态估计方式为：

其中，H_2,k是下游估计过程在k时刻的量测矩阵，是对含水率的最终估计结果，其均值即为测量时间段内平均相含率的最优估计。

本发明的实质性特点是：针对多相流流动过程，采用上、下游相含率传感器测量所得相含率进行估计。利用测量对象的上、下游含率状态并行测量的特点，将估计过程分为一次线性滤波和一次非线性滤波，对上游传感器测量数据进行线性卡尔曼滤波，计算其实时的状态转移矩阵，将该矩阵作为下游传感器状态转移矩阵的先验信息，对下游传感器测试数据的进行基于无迹卡尔曼滤波的非线性状态估计，随估计过程改变Sigma点集而进行自适应卡尔曼估计。本发明的有益效果及优点如下：

1、充分挖掘两相流上、下游相含率传感器测量数据之间的内在联系，实现同类传感器测试信息的有机融合；

2、增加了对于两相流对象状态转换的先验信息；

3、利用不同传感器测量的互补性，提高了估计精度；

4、计算量较小，精度较好，可以实现在线估计。

附图说明

以下附图描述了本发明所选择的实施例，均为示例性附图而非穷举或限制性，以电导传感器测量油水两相流为例，其中：

图1本发明实施例采用的电导传感器结构示意图。

图2本发明的在线自适应卡尔曼估计算法计算步骤。

具体实施方式

以电导传感器为具体实施例，不局限于电导传感器，其它类型的传感器同样适用。

当信号频率在10kHz到100kHz之间时，油水两相流的负载特性体现为电阻性，即电极对之间的电压与分布其间的流体阻抗成正比。在油水两相流流体可以形成回路时，油水两相流的阻抗与相分布有对应关系，因此可以通过测量电极对之间的电压来测量其相含率。

如图1所示，由六个沿轴向排列的环形电极组成，其中1、6为激励电极，2、3、4、5为测量电极。采用电流激励电压测量的方式，通过向电极1中注入恒定电流并将电极6接地，可在电极1和6中间形成稳定的电学测量场，分别测量电极对2-3、电极对3-4以及电极对4-5之间的电势差V可获得其间的两相流平均电阻率。详细流程如图2所示，

电导传感器为6环结构电导环传感器阵列以轴向排列的6个环状金属电极嵌入管路中的形式构成，管径D为50mm，相邻2个环状电极的距离分别为：1D、0.4D、1.2D、0.4D、1D。在此，依照图2所示流动方向，选取电极2和电极3之间的区域为上游测量单元，电极4和电极5之间的区域为下游测量单元，两个测量单元完全相同。两侧环形电极为激励电极，连接峰峰值2mA，20kHz交流方波电流信号，中间电极为测量电极，采集电势差。管道中流体含水率不同时，流体的混合电导率会发生变化，因此每一组电势差反映的是该两环形电极之间的测量区域内部的介质分布。但是，当管道内连续相不导电，且为分散流时，电导传感器无法形成电流回路，测量失效。测量实验过程中，定义无量纲电导，

其中，V₀是管道充满水时的电势差测量值，V_i是管道内流过两相流时的电势差测量值。水平管道油水两相流的流型主要分为分层流和分散流。对于分层流，V^*即为相含率。对于如油包水或水包油的分散流态，根据Maxwell理论推导，含水率与无量纲电导有如下关系：

其中，h_w代表含水率，作为估计过程中的量测向量z。

将来自电导传感器上、下游的测量数据进行时间配准，具体方式是将来自上游、下游的含水率数据进行离散形式下的互相关计算：

上游和下游测量单元的含率测量值分别为z₁(j)和z₂(j)，其中j＝1,2,3,…N为测量时间内的数据点，z₁(j)和z₂(j)的互相关函数为 的最大值出现在n₀处，对应离散系统中上、下游测量数据平移组数n₀。z_1,k与是对同一被测对象在不同时刻和不同位置处的观测。

图2为本例自适应卡尔曼估计方法的计算流程。

1)本发明在估计过程中传感器系统的量测矩阵H取值为1。

2)本算法先后进行两次基于卡尔曼滤波的状态估计。上游数据的线性估计过程获取对被测流体相含率的初始估计，获得被测对象状态转移矩阵，将其作为先验信息代入下游数据的非线性估计过程中：

3)在上游状态转移序列的基础上，构造非线性滤波状态算法中的Sigma点集。其目的主要有两点：一是增加油水两相流在随机过程角度上的描述，减轻由于实验系统采样频率限制带来的信息损失；二是通过调整无迹卡尔曼滤波过程中的Sigma点集构成方式，一定程度上抵消实验过程中非高斯噪声对含率测量结果的影响:

本专利提出的Sigma点集是围绕在状态估计值周围的一定数量的采样点，用以进行非线性变换后逼近状态真值，Sigma点集基本构成方式为而为了应对存在的非高斯噪声，本专利另提出多种Sigma点集的构造方式：

加性噪声模型：

减性噪声模型：

加性噪声模型和减性噪声模型分别处理噪声成分主要为加性噪声或减性噪声的动态系统。除此之外，还可以根据线性滤波的一次估计状态进行判断，实时更新估计过程中每次迭代的Sigma点集的构造方式：

或

其中，式(27)是自适应的减性噪声模型，式(28)是自适应混合噪声模型。L是状态向量的维度，是下游状态估计在k-1时刻的状态估计，p_2，k-1是k-1时刻的全局协方差，参数γ的数值在式(1)中已给出。

4)构造Sigma点集后，对下游含水率测量值进行非线性滤波。

将Sigma点集中的点进行加权平均，进行状态预测，其权值分配方式如式(29)所示，其中i为矩阵在Sigma点集中对应的元素序号。

状态预测的方式为：

基于无迹卡尔曼的估计过程中，存在全局的误差协方差矩阵、测量值的协方差矩阵以及测量值与状态估计的互协方差矩阵，它们都随着预测和估计的过程中更新，更新方式如下：

其中，R_2,k为下游估计过程中的系统量测噪声的协方差，K_2,k为下游测量值估计过程中的非线性卡尔曼增益，其更新方式为：

则含水率状态估计的更新为

其中，是对含水率的最终估计结果。

5)在完成每组含水率条件下所有实验数据的滤波和估计之后，对所有时刻的自适应估计值进行平均，其表示为作为该测试时间段内平均相含率的最优估计值。

Claims (1)

1.一种两相流含水率自适应串联估计方法，利用上、下游不同传感器的测量数据进行估计，包括如下步骤：

1)对上游传感器和下游传感器测量数据进行时间配准：

获取上游传感器和下游传感器所获相含率时间序列z₁(j)和z₂(j)，j＝1,2,3,…N为测量时间内的数据点，下标1和2分别代表上游传感器和下游传感器，计算渡越时间，利用渡越时间结合数据采样频率，对上游传感器和下游传感器测量数据进行时间配准；

2)状态估计算初始化：

${\hat{x}}_{1,0}=\frac{1}{N}\underset{k=1}{\overset{N}{\Σ}}{z}_{1,k}---\left(1\right)$

${\hat{p}}_{1,0}=\underset{k=1}{\overset{N}{\Σ}}\sqrt{\frac{z_1\left(k\right)-{\stackrel{\‾}{z}}_1}{N-1}}---\left(2\right)$

${\hat{x}}_{2,0}=\frac{1}{N}\underset{k=1}{\overset{N}{\Σ}}{z}_{2,k}---\left(3\right)$

其中，是上游状态向量的估计初值，是上游协方差的迭代初值，z_1,k是上游传感器在k时刻的测量值，且 是下游状态向量的估计初值；

3)对上游传感器含率测量值进行线性卡尔曼状态估计和更新，获得对状态向量的预测和协方差矩阵的预测；

在估计过程中量测矩阵H取值均为1，状态向量x₁代表含水率，k时刻的系统模型为：

x_1,k+1＝F_1,kx_1,k+v_1,k (4)

z_1,k＝H_1,kx_1,k+w_1,k (5)

其中，F_1,k是k时刻的状态转移矩阵，x_1,k+1和x_1,k分别是上游在k+1时刻和k时刻的状态向量，H_1,k是k时刻的量测矩阵，z_1,k是量测向量，v_1,k和w_1,k分别表示系统在k时刻的过程噪声和量测噪声；则状态向量的预测的预测和协方差矩阵的预测为：

${\hat{x}}_{1,k|k-1}=F_{1,k}{\hat{x}}_{1,k-1}---\left(6\right)$

${\hat{p}}_{1,k|k-1}=F_{1,k-1}{\hat{p}}_{1,k-1}{F_{1,k-1}}^T+Q_{1,k-1}---\left(7\right)$

根据上游量测信息z_1,k修正估计值：

${\hat{p}}_{1,k}=\[I-K_{1,k}{H}_{1,k}\]{\hat{p}}_{1,k|k-1}---\left(9\right)$

其中，K_1,k是上游传感器系统在k时刻的卡尔曼增益：

$K_{1,k}={\hat{p}}_{1,k|k-1}{H_{1,k}}^T{\[H_{1,k}{\hat{p}}_{1,k|k-1}{H_{1,k}}^T+R_{1,k}\]}^{-1}---\left(10\right)$

其中，H_1,k是k时刻上游传感器量测矩阵，是上游传感器在k时刻对状态的预测值，是上游传感器k时刻对状态x_1,k的最小均方误差估计，是k时刻估计的协方差，R_1,k为上游滤波过程在k时刻系统量测噪声的协方差，可取值为0.01至0.2；Q_1,k-1为上游滤波过程在k-1时刻系统过程噪声的协方差，可取值为0.01至0.09，I为单位矩阵，其维数等同于状态向量的维数；

3)利用下游传感器测量数据，计算两相流过程状态转换矩阵：

将传感器上游测量单元确定的状态转移过程加以记录，用于时间配准后下游含水率状态估计中：

$F_{2,k}=\frac{{\hat{x}}_{1,k-n_0}}{{\hat{x}}_{1,k-n_0-1}}---\left(11\right)$

其中，F_2,k为下游测量单元的状态转移矩阵，k＝n₀,n₀+1,…,N；

4)对下游含率测量值进行基于无迹卡尔曼滤波非线性状态估计：

构造Sigma点集，其基本构造方式如式(12)所示：

其中，L是状态向量的维度，是下游估计过程k-1时刻的状态估计，p_2,k-1是下游估计过程k-1时刻的全局协方差，γ与λ可由下式计算：

$\left\{\begin{array}{c}\mathrm{\γ}=\sqrt{L+\mathrm{\λ}}\\ \mathrm{\λ}={\mathrm{\α}}^2(L+\mathrm{\κ})-L\end{array}\right.---\left(13\right)$

其中，λ为规模因子，α是决定Sigma点集中点的分散程度扩展因子，取值为0到1之间；κ是缩放系数，设为0，计算状态预测权系数W_i ^(m)和协方差预测权系数W_i ^(c)：

$\left\{\begin{array}{c}{W}_0^{\left(m\right)}=\mathrm{\λ}/(L+\mathrm{\λ})\\ W_0^{\left(c\right)}=\mathrm{\λ}/(L+\mathrm{\λ})+(1-{\mathrm{\α}}^2+\mathrm{\β})\\ W_i^{\left(m\right)}=W_i^{\left(c\right)}=1/\[2(L+\mathrm{\λ})\]\end{array}\right.---\left(14\right)$

其中，β先验分布因子，通常设置为2，i为矩阵在Sigma点集中对应的元素序号，下游状态估计过程的状态预测方式如下：

$\left\{\begin{array}{c}{\mathrm{\χ}}_{i,k|k-1}^{*}=f\left({\mathrm{\χ}}_{i,k-1}\right)\\ {\hat{x}}_{2,k|k-1}=\underset{i=0}{\overset{2L}{\mathrm{\Σ}}}{W}_i^{\left(m\right)}{\mathrm{\χ}}_{i,k|k-1}^{*}\end{array}\right.---\left(15\right)$

$f\left({\mathrm{\χ}}_{i,k-1}\right)=F_{2,k}{\mathrm{\χ}}_{i,k-1}=\frac{x_{1,k-n_0}}{x_{1,k-n_0-1}}{\mathrm{\χ}}_{i,k-1}---\left(16\right)$

估计过程中的协方差矩阵分为全局协方差矩阵和互协方差矩阵，更新方式分别如下：

$\left\{\begin{array}{c}{p}_{2,k|k-1}=\underset{i=0}{\overset{2L}{\mathrm{\Σ}}}{W}_i^{\left(c\right)}\left({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{2,k|k-1}\right){\left({\mathrm{\χ}}_{i,k|k-1}^{*}-{\hat{x}}_{2,k|k-1}\right)}^T+H_{2,k-1}{\hat{x}}_{2,k-1}{Q}_{2,k-1}{\left(H_{2,k-1}{\hat{x}}_{2,k-1}\right)}^T\\ P_{2,k}=P_{2,k|k-1}-K_{2,k}{P}_{{\hat{z}}_k{\hat{z}}_k}{K}_{2,k}^T\end{array}\right.---\left(17\right)$

其中，p_2,k|k-1是对全局协方差的预测，Q_2,k-1是下游过程噪声的协方差，H_2,k是下游估计过程在k时刻的量测矩阵；

$\left\{\begin{array}{c}{P}_{{\hat{x}}_k{\hat{z}}_k}={\hat{p}}_{2,k|k-1}{H}_{2,k}^T\\ P_{{\hat{z}}_k{\hat{z}}_k}=H_{2,k}{\hat{p}}_{2,k|k-1}{H}_{2,k}^T+R_{2,k}\end{array}\right.---\left(18\right)$

其中，是下游传感器测量值在k时刻的协方差，是k时刻互协方差矩阵，R_2,k为下游估计过程中的系统量测噪声的协方差，取值为0.01至0.09；K_2,k是k时刻下游估计过程的卡尔曼增益，按照式(19)方式进行更新；

$K_{2,k}=P_{{\hat{x}}_k{\hat{z}}_k}{P}_{{\hat{z}}_k{\hat{z}}_k}^{-1}---\left(19\right)$

基于无迹卡尔曼的非线性状态估计方式为：

$\left\{\begin{array}{c}{\hat{z}}_{2,k|k-1}=\underset{i=0}{\overset{2L}{\mathrm{\Σ}}}{H}_{2,k}{\hat{x}}_{2,k|k-1}\\ {\hat{x}}_{2,k}={\hat{x}}_{2,k|k-1}+K_{2,k}(z_{2,k}-{\hat{z}}_{2,k|k-1})\end{array}\right.---\left(20\right)$

其中，H_2,k是下游估计过程在k时刻的量测矩阵，k＝n₀,n₀+1,…,N，是对含水率的最终估计结果，其均值即为测量时间段内平均相含率的最优估计。