CN105842298A - 两相流含水率的自适应估计方法 - Google Patents

Abstract

本发明涉及一种两相流含水率自适应卡尔曼估计融合方法，包括：计算能够直接反映含水率的传感器测量值作为含水率突变的判断标准；卡尔曼估计的工作状态为稳定状态和突变状态两个工作状态，需要设定阈值与含水率突变判断参数进行比较，对含水率是否突变做出判断，以切换卡尔曼估计的工作状态；对每个传感器进行卡尔曼估计更新；根据量测信息修正估计值；在稳定状态和突变状态下，应用不同的算法，对自适应权重做归一化计算；计算融合后的含水率。本发明能够提升分含水率的估计精度。

Description

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

技术领域

本发明属于流体测量技术领域，利用传感器测量得到的含水率冗余信息进行自适应卡尔曼估计融合，用于提高两相流含水率测量精度。

技术背景

两相流广泛存在于日常生活与工业生产过程中，尤其体现在食品加工，生物工程，化工产业，冶金工业和石油产业等行业，包括油水两相流和气液两相流等。与一般单相流相比，两相流的流动状态更为复杂，因此两相流的在线过程参数检测一直以来都是学者和工业产业界关注的重点。对于石油产业而言，油水两相流含水率对油井监测、生产预测和开采计划制定具有重要意义。

目前针对两相流含水率测量，可以使用快关阀方法、射线方法、微波方法、Wire-mesh方法、超声方法和电阻抗方法等。相比于其他方法，电学方法具有测量速度快，对含水率变化敏感，传感器结构简单，易于安装，价格低廉等优点。电学含水率测量方法包括电导法、电容法、层析成像法等。单一传感器对含水率测量有一定限制，例如，针对油水两相流，当使用电导传感器和电容传感器联合测量时，电导传感器无法测量油相为连续相的分散流，电容传感器对层状流测量分辨率较低等。在含水率较高(大于30％)且流型不是层状流时，电导传感器和电容传感器均可获得含水率测量值，因此可利用双传感器的冗余数据进行融合，改善测量精度。

两相流是一种随机的动态流动过程，流动状态复杂多变，传感器的测量信号又有无法避免的随机误差，其测量信息并不准确，而且只能反映系统的外部特性。因此，很难根据传感器的测量信号直接对多相流动状态进行精确的建模，而对状态的描述都要基于大量的近似假设，因此为得到多相流的含水率信息，需要根据测量数据对动态系统的内部状态进行估计。

状态估计对认知和分析测量系统具有重要的应用价值，常用的估计理论包括最小二乘估计、线性最小方差估计和卡尔曼估计等。不同的估计方法性能指标不同，使用的信息不同，其适用范围也不同。最小二乘估计通常应用于常值向量和随机向量的估计，精度不高，算法简单；线性最小方差估计以均方误差最小作为优化参数，适用于平稳过程的估计，其测量信息越丰富，估计精度越高，计算量越大，而且会随着估计过程的累积，使得计算量越来越大，不适用于复杂的长时间的估计过程。针对动态非平稳的复杂的多相流流动过程，卡尔曼估计采用递推的形式，通过一个状态模型和一个观测模型的迭代，利用上一时刻的估计值和这一时刻的观测值，来对这一时刻的状态进行估计。卡尔曼估计中，单个传感器(k+1)时刻的状态模型为：x(k+1)＝F(k)x(k)+v(k)，式中x(k)和x(k+1)分别表示传感器在k和(k+1)时刻被测对象的状态矢量；v(k)为状态过程的高斯白噪声；F(k)是状态转移矩阵。传感器k时刻的测量模型为z(k)＝h(k)x(k)+w(k)，式中，z(k)是传感器的测量值；h(k)是传感器的测量矩阵；w(k)是传感器的测量噪声。卡尔曼估计计算量较小，可实时估计系统状态，适用于动态系统。

多相流作为非线性的动态系统，流态复杂，含水率会随着流态的不同而产生变化，难以构建理论模型，因此有必要对卡尔曼估计模型的参数进行选择和优化，以使其能够对实际的多相流状态进行估计。当流态稳定时，含水率较为稳定，可以利用普通的卡尔曼估计方法进行估计，当流态发生突变时，含水率也会随之发生突变，此时需要对卡尔曼估计方法进行改进，以便能够更快速地适应流态突变。

针对多相流的随机动态过程，不同的估计算法具有不同的问题，如低估计精度、低响应速度等。因此，利用自适应加权融合算法对不同估计算法的估计值进行融合，以避免各自算法的问题，从而获得更好的结果。单个传感器的估计结果，集中式估计结果，分布式估计结果均可以用于进行自适应加权融合，利用各种方法的优势，提高含水率输出的精度。

发明内容

本发明针对两相流含水率测量，特别是石油工业的采油过程中经常发生的含水率突变状态，需要估计算法能够对突变状态进行高精度的快速响应，提出一种基于卡尔曼滤波估计的含水率动态估计算法，能够适应测量过程中的含水率突变状态，并利用自适应加权方法，提升分含水率的估计精度。本发明的具体技术方法如下：

一种两相流含水率自适应卡尔曼估计融合方法，设传感器为H个，k时刻第j个传感器测量所得含水率为z_j(k)，包括：

1)计算能够直接反映含水率的传感器测量值作为含水率突变的判断标准J(k)：

$J\left(k\right)=|z\left(k\right)-\frac{z(k-1)+z(k-2)}{2}|---\left(1\right)$

式中，是k时刻多种传感器含水率测量值的平均值；

2)卡尔曼估计的工作状态为稳定状态和突变状态两个工作状态，需要设定阈值与含水率突变判断参数J(k)进行比较，对含水率是否突变做出判断，以切换卡尔曼估计的工作状态：

$Q=\left\{\begin{array}{cc}{Q}_1& ifJ\left(k\right)<\mathrm{\&theta;}\\ Q_2& ifJ\left(k\right)>\mathrm{\&theta;}\end{array}\right.---\left(2\right)$

式中，Q是过程噪声协方差矩阵，θ是阈值，在计算中为一固定值，其典型取值范围为：在需要对油水两相流含水率变化较为敏感的条件下，可在0％<θ≤5％范围内取一固定值；在对油水两相流含水率变化较为不敏感的条件下，可在5％<θ≤50％范围内取一固定值。Q₁和Q₂分别是稳定状态和突变状态时的估计得到的过程噪声协方差，将J(k)与θ的比较作为判断标准以切换卡尔曼估计器的工作状态：当J(k)<θ时，判断卡尔曼估计器工作于稳定状态；当J(k)>θ时，判断卡尔曼估计器工作于突变状态；

3)在判断工作状态后，利用不同的Q值，对每个传感器使用如下方程进行卡尔曼估计更新：

根据k时刻状态推出(k+1)时刻状态：

${\hat{x}}_j(k+1|k)=F(k+1){\hat{x}}_j\left(k\right|k)---\left(3\right)$

p_j(k+1|k)＝F(k)p_j(k|k)F^T(k)+Q (4)

根据量测信息修正估计值：

其中，矩阵K_j(k+1)是第j个传感器(k+1)时刻的卡尔曼增益：

式中，是k时刻第j个传感器对被测状态x(k+1)的预测值，F(k+1)是状态转移矩阵，是k时刻第j个传感器对被测状态x(k)的最小均方误差估计值，p_j(k+1|k)是k时刻第j个传感器对(k+1)时刻的估计的协方差，p_j(k|k)是k时刻第j个传感器对k时刻估计的协方差，z_j(k+1)是(k+1)时刻第j个传感器含水率测量值，h_j(k+1)是(k+1)时刻第j个传感器的测量矩阵，I是单位矩阵，R是测量噪声协方差矩阵；

4)当J(k)<θ时，分别用不同传感器的测量数据进行卡尔曼集中式估计，利用不同传感器的估计值进行卡尔曼分布式估计，分别获得k时刻含水率的集中式估计值和分布式估计值；

当J(k)>θ时，将不同传感器的测量数据进行卡尔曼集中式估计，获得k时刻含水率的集中式估计值；

5)当J(k)<θ时，使用各个传感器前一时刻卡尔曼算法、集中式估计算法和分布式估计算法得出的估计值计算自适应权重；当J(k)>θ时，使用各个传感器的前一时刻卡尔曼估计算法和集中式估计算法得出的估计值计算自适应权重。

对于第q种算法，k时刻自适应权重为W_q(k)，其计算方法为：

$W_q\left(k\right)=\sqrt{{\mathrm{\&sigma;}}_q^2+{\left(x_q\right(k-1)-\stackrel{\&OverBar;}{x}(k-1\left)\right)}^2}---\left(8\right)$

式中，σ_q是第q种算法的“不确定度补偿”参数，x_q(k-1)是(k-1)时刻第q种算法的估计值，是(k-1)时刻所使用的所有估计算法估计值的平均值，m是用于计算的估计算法的个数。对于单个传感器的卡尔曼估计，σ_q使用该传感器的满量程测量精度σ_j；对于集中式估计算法和分布式估计算法的“不确定度补偿”参数σ_q，使用如下方程计算：

${\mathrm{\&sigma;}}_q=\sqrt{\underset{j=1}{\overset{H}{\&Sigma;}}{\mathrm{\&sigma;}}_j^2}---\left(9\right)$

在稳定状态和突变状态下，应用不同的算法，对k时刻自适应权重W_q(k)做归一化计算：

$W_q^{*}\left(k\right)=\frac{W_q\left(k\right)}{{\mathrm{\&Sigma;}}_{q=1}^m{W}_q\left(k\right)}---\left(10\right)$

其中是在k时刻第q种算法的归一化权重。

6)在稳定状态与突变状态，均使用如下方程计算融合后的含水率：

$X\left(k\right)=\underset{q=1}{\overset{m}{\&Sigma;}}{x}_q\left(k\right)\&times;W_q^{*}\left(k\right)---\left(11\right)$

X(k)为k时刻融合后的含水率。

作为优选实施方式，在需要对油水两相流含水率变化较为敏感的条件下，可在0％<θ≤5％范围内取一固定值；在对油水两相流含水率变化较为不敏感的条件下，可在5％<θ≤50％范围内取一固定值。测量噪声采用零均值高斯白噪声，计算该白噪声的协方差矩阵作为R；测得的信号局部方差最小值认为是噪声方差，最大值认为是信号方差，这两个方差的比值认为是信噪比，使用当前测量k时刻之前的检测数据作为局部数据，根据信噪比重新构造随机噪声信号，并对随机噪声归一化，计算稳定状态测量噪声方差，取其最大值作为Q₁值；计算突变状态测量噪声方差，取其最大值作为Q₂值；状态转移矩阵F(k)和传感器的测量矩阵h(k)均取值为1。

本发明的实质性特点是：针对多相流流动过程，采用多传感器测量所得数据，在测量精度难以提高的情况下，利用可以对被测对象突变状态有快速响应的改进卡尔曼估计算法的估计值，集中式估计和分布式估计的估计值，进行自适应卡尔曼融合估计。利用“不确定度补偿”参数，对不同估计精度的估计方法的自适应权重进行补偿，提高估计精度。本发明的有益效果及优点如下：

1、该方法对被测对象的突变状态有更快速的响应；

2、该方法相对于传感器直接测量，提高了估计精度；

3、该方法在计算量较小，可以实现实时在线估计。

附图说明

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

图1本发明的测量方法电容电导传感器测量空间示意图；

图2本发明的测量方法整体计算步骤；

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

具体实施方式

下面结合说明书附图说明本发明的计算方法。

本例以电容传感器和电导传感器两种测量方法对油水两相流进行含水率测量。图1为本发明的测量方法中，示例的电容电导传感器测量空间示意图。电导传感器为四环结构，两侧环形电极为激励电极，连接峰峰值2mA，20kHz交流方波电流信号，中间两环形电极为测量电极，采集电势差。管道中流体含水率不同时，流体的混合电导率会发生变化，测得的电势差可以直接反映管道中流体的含水率。但是，当管道内连续相不导电，且为分散流时，电导传感器无法形成电流回路，测量失效。电容传感器为对壁结构，置于电导传感器两测量电极之间。一个极板为激励电极，连接峰峰值20V，1MHz的交流正弦电压信号，另一个极板为测量电极，利用与被测电容值成正比的测量电路，采集电压值。管道中流体含水率不同时，流体的混合介电常数会发生变化，测得的电压值可以直接反映管道中流体的含水率。但是，当管道内流体为层状流或近似层状流时，测量敏感度降低，测量失效。电容电导传感器测量敏感空间相同，采集频率相同，可以较为方便地用于数据融合。

图2为本例中，含水率的整体计算步骤，包括利用互补信息和提取用于融合的冗余信息。在k时刻，电容传感器和电导传感器均可以得到关于含水率的测量值：(1)当电导传感器含水率模型计算值z₁(k)大于100％时，说明非导电相为连续相，采用电容传感器含水率计算值z₂(k)；(2)当电容传感器含水率模型计算值z₂(k)与标定时的满管水测量值z_w满足90％×z_w≤z₂(k)≤110％×z_w时，说明管道内为层状流或近似层状流，采用电导传感器含水率计算值。除此两种情况外，两种传感器所获含水率，均可用于数据融合，利用在线自适应卡尔曼估计融合方法进行关于含水率的融合计算。

图3为本例中，在线自适应卡尔曼估计融合方法的计算步骤。

1)针对多相流流动，特别是石油工业的采油过程中经常发生的含水率突变状态，需要估计算法能够对突变状态进行高精度的快速响应。对于电容传感器和电导传感器都可以计算得到的有效的含水率值，其两者的结果进行平均，依据式(1)计算含水率突变判断参数。

2)利用阈值θ和含水率突变判断参数J(k)对含水率突变进行判断，依据式(2)将J(k)与阈值θ的比较作为判断标准以切换卡尔曼估计器的工作状态：当J(k)<θ时，令卡尔曼估计器工作于稳定状态；当J(k)>θ时，令卡尔曼估计器工作于突变状态。

在需要对油水两相流含水率变化较为敏感的条件下，可在0％<θ≤5％范围内取一固定值；在对油水两相流含水率变化较为不敏感的条件下，可在5％<θ≤50％范围内取一固定值。本例中，令阈值θ为3％。Q₁和Q₂是稳定状态和突变状态时的估计得到的过程噪声协方差。

3)当J(k)<θ时，认为多相流动处于稳定状态，计算稳定状态电容传感器和电导传感器各自的卡尔曼增益，通过式(3)～(7)，利用Q₁进行卡尔曼估计。当J(k)>θ时，认为多相流动处于突变状态，计算稳定状态电容传感器和电导传感器各自的卡尔曼增益，通过式(3)～(7)，利用Q₂进行卡尔曼估计。卡尔曼估计器的测量噪声R采用零均值高斯白噪声。测得的信号局部方差最小值认为是噪声方差，最大值认为是信号方差，这两个方差的比值作为信噪比。使用当前测量k时刻之前1秒的检测数据作为局部数据(计算初始值所需的之前1秒数据可以由标定的历史数据得到)，重新构造噪声信号，计算出稳定状态测量噪声方差Q₁，Q₁计算值为6～10之间，采用最大值10；采用同样方法计算突变状态测量噪声方差Q₂，Q₂计算值小于0.1，采用最大值0.1。当电导传感器和电容传感器均正常工作时，状态转移矩阵F(k)取值为1。电导传感器和电容传感器测量得到的含水率为真实含水率和量测噪声的叠加，传感器的测量矩阵h(k)取值为1。第j个传感器k时刻的状态估计值使用初始时刻的测量值作为初值。

4)a)对于处于稳定状态的多相流动过程，分别利用不同传感器的测量数据进行卡尔曼集中式估计和分布式估计，分别获得k时刻含水率的集中式估计值和分布式估计值。

使用如下方程进行卡尔曼集中式估计：

${\hat{X}}_{cf}\left(k\right|k-1)=F_{cf}(k-1){\hat{X}}_{cf}(k-1|k-1)---\left(13\right)$

式中，角标cf表示集中式估计。是k时刻被测状态x(k)的集中式全局估计值，是(k-1)时刻对被测状态x(k)的集中式全局预测值，K_j(k|k)是k时刻第j个传感器的卡尔曼增益，z_j(k)是k时刻第j个传感器含水率测量值，h_j(k)是第k时刻j个传感器的测量矩阵，F_cf(k-1)是(k-1)时刻的集中式转移矩阵。

使用如下方程进行卡尔曼分布式估计：

${\hat{X}}_{df}\left(k\right|k-1)=F_{df}\left(k\right){\hat{X}}_{df}(k-1|k-1)---\left(15\right)$

P_df(k|k-1)＝F_df(k-1)P_df(k-1|k-1)F_df ^T(k-1)+Q (17)

式中，角标df表示分布式估计。是k时刻被测状态x(k)的分布式全局估计，是(k-1)时刻被测状态x(k)的分布式全局预测，P_df(k|k)是k时刻分布式全局估计协方差，p_j(k|k)是第j个传感器k时刻对k时刻估计的协方差，p_j(k|k-1)是第j个传感器(k-1)时刻对k时刻的估计的协方差，是第j个传感器(k-1)时刻对被测状态x(k)的预测值，F_df(k)是k时刻的分布式转移矩阵，P_df(k|k-1)是(k-1)时刻对k时刻的分布式全局估计协方差。

本例中，利用得到的电容传感器和电导传感器k时刻含水率测量值进行计算，k时刻的分布式转移矩阵F_df(k)与集中式转移矩阵F_cf(k)取值均为1。

b)对于处于突变状态的多相流动过程，利用不同传感器的测量数据进行卡尔曼集中式估计，获得k时刻含水率的集中式估计值。本例利用得到的电容传感器和电导传感器k时刻测量值，依据式(12)～(13)进行集中式估计，k时刻的集中式转移矩阵F_cf(k)取值为1。突变状态发生后的一段时间内(本例为突变发生后3秒)，均使用估计器的突变工作状态进行计算，之后再使用估计器的稳定状态进行计算。

5)计算各个算法的自适应权重，其中“不确定度补偿”参数σ用于调节各个算法的自适应权重，补偿参数根据各个算法的估计精度选取，以补偿不同算法由于精度不同而造成的可信度不同。在计算自适应权重时，各个传感器的“不确定度补偿”参数σ_q，使用其满量程测量精度σ_j，集中式估计值和分布式估计值的“不确定度补偿”参数使用式(9)计算。在本例中，传感器个数H取2(电容传感器、电导传感器)。利用计算得到集中式估计值和分布式估计值的“不确定度补偿”参数进行进一步自适应权重计算。例如：电容传感器估计“不确定度补偿”参数为5％，电导传感器估计“不确定度补偿”参数为5％，集中式估计“不确定度补偿”参数为7.07％，分布式估计“不确定度补偿”参数为7.07％。

利用式(8)计算自适应权重，利用式(10)进行自适应权重归一化计算，得到最终的归一化权重。对于多相流动稳定状态，利用在(k-1)时刻得到的电容传感器估计值、电导传感器估计值、集中式估计值、分布式估计值，计算四种算法k时刻各自的归一化权重；对于多相流动突变状态，利用在(k-1)时刻得到的电容传感器估计值、电导传感器估计值、集中式估计值，计算三种算法k时刻各自的归一化权重。初始时刻各种算法采用的估计值，均为电容传感器和电导传感器测量值的平均值。

6)对于多相流动的稳定状态和突变状态，均使用式(11)计算在线自适应加权卡尔曼估计融合后的含水率。对于稳定状态，使用各个传感器的卡尔曼估计值、卡尔曼集中式估计值和卡尔曼分布式估计值及其相应的归一化权重计算融合后的含水率；对于突变状态，使用各个传感器的卡尔曼估计值和卡尔曼集中式估计值及其相应的归一化权重计算融合后的含水率。

Claims (1)

1.一种油水两相流含水率自适应卡尔曼估计融合方法，设传感器为H个，k时刻第j个传感器测量所得含水率为z_j(k)，包括：

1)计算能够直接反映含水率的传感器测量值作为含水率突变的判断标准J(k)：

$J\left(k\right)=|z\left(k\right)-\frac{z(k-1)+z(k-2)}{2}|---\left(1\right)$

式中，是k时刻多种传感器含水率测量值的平均值；

2)卡尔曼估计的工作状态为稳定状态和突变状态两个工作状态，需要设定阈值与含水率突变判断参数J(k)进行比较，对含水率是否突变做出判断，以切换卡尔曼估计的工作状态：

$Q=\left\{\begin{array}{ccc}{Q}_1& if& J\left(k\right)<\mathrm{\&theta;}\\ Q_2& if& J\left(k\right)>\mathrm{\&theta;}\end{array}\right.---\left(2\right)$

式中，Q是过程噪声协方差矩阵，θ是阈值，在计算中为一固定值，其典型取值范围为：在对油水两相流含水率变化较为敏感的实验条件下，可在0％<θ≤5％范围内取一固定值；在油水两相流含水率变化较为不敏感的工程条件下，可在5％<θ≤50％范围内取一固定值。Q₁和Q₂分别是稳定状态和突变状态时的估计得到的过程噪声协方差，将J(k)与θ的比较作为判断标准以切换卡尔曼估计器的工作状态：当J(k)<θ时，判断卡尔曼估计器工作于稳定状态；当J(k)>θ时，判断卡尔曼估计器工作于突变状态；

3)在判断工作状态后，利用不同的Q值，对每个传感器使用如下方程进行卡尔曼估计更新：

根据k时刻状态推出(k+1)时刻状态：

${\hat{x}}_j(k+1|k)=F(k+1){\hat{x}}_j\left(k\right|k)---\left(3\right)$

p_j(k+1|k)＝F(k)p_j(k|k)F^T(k)+Q (4)

根据量测信息修正估计值：

其中，矩阵K_j(k+1)是第j个传感器(k+1)时刻的卡尔曼增益：

式中，是k时刻第j个传感器对被测状态x(k+1)的预测值，F(k+1)是状态转移矩阵，是k时刻第j个传感器对被测状态x(k)的最小均方误差估计值，p_j(k+1k)是k时刻第j个传感器对(k+1)时刻的估计的协方差，p_j(k|k)是k时刻第j个传感器对k时刻估计的协方差，z_j(k+1)是(k+1)时刻第j个传感器含水率测量值，h_j(k+1)是(k+1)时刻第j个传感器的测量矩阵，I是单位矩阵，R是测量噪声协方差矩阵；

4)当J(k)<θ时，分别用不同传感器的测量数据进行卡尔曼集中式估计，利用不同传感器的估计值进行卡尔曼分布式估计，分别获得k时刻含水率的集中式估计值和分布式估计值；

当J(k)>θ时，将不同传感器的测量数据进行卡尔曼集中式估计，获得k时刻含水率的集中式估计值；

5)当J(k)<θ时，使用各个传感器前一时刻卡尔曼算法、集中式估计算法和分布式估计算法得出的估计值计算自适应权重；当J(k)>θ时，使用各个传感器的前一时刻卡尔曼估计算法和集中式估计算法得出的估计值计算自适应权重。

对于第q种算法，k时刻自适应权重为W_q(k)，其计算方法为：

式中，σ_q是第q种算法的“不确定度补偿”参数，x_q(k-1)是(k-1)时刻第q种算法的估计值，是(k-1)时刻所使用的所有估计算法估计值的平均值，m是用于计算的估计算法的个数。对于单个传感器的卡尔曼估计，σ_q使用该传感器的满量程测量精度σ_j；对于集中式估计算法和分布式估计算法，使用如下方程计算：

${\mathrm{\&sigma;}}_q=\sqrt{\underset{j=1}{\overset{H}{\&Sigma;}}{\mathrm{\&sigma;}}_j^2}---\left(9\right)$

在稳定状态和突变状态下，应用不同的算法，对k时刻自适应权重W_q(k)做归一化计算：

$W_q^{*}\left(k\right)=\frac{W_q\left(k\right)}{{\mathrm{\&Sigma;}}_{q=1}^m{W}_q\left(k\right)}---\left(10\right)$

其中是在k时刻第q种算法的归一化权重。

6)在稳定状态与突变状态，均使用如下方程计算融合后的含水率：

$X\left(k\right)=\underset{q=1}{\overset{m}{\&Sigma;}}{x}_q\left(k\right)\&times;W_q^{*}\left(k\right)---\left(11\right)$

X(k)为k时刻融合后的含水率。