本发明涉及一种油水两相流含水率自适应估计方法，包括：将来自电导传感器上、下游的测量数据进行时间配准；基于自适应更新的状态估计；判定两相流是否处于水连续的状态；对传感器的含率测量值进行卡尔曼状态估计和更新。本发明能够提升含水率的估计精度。
1.一种油水两相流含水率自适应估计方法，利用阵列式电导传感器上下游电极测量数据进行估计，包括如下步骤：
1)将来自电导传感器上、下游的测量数据进行时间配准：
将取自上游和下游测量单元的含率测量值z1(n)和z2(n)看作来自各态历经的平稳随机过程的两个样本函数，其中n＝1,2,3,…N为测量时间内的数据点，两传感器间距为L，设渡越时间为τ0，对应离散系统中上、下游测量数据平移组数n0；计算z1(n)和z2(n)在离散条件下的互相关函数 公式如下：Rz1z2(n)的最大值出现在n0处，即互相关函数峰值对应流体从上游测量单元到下游测量单元的渡越时间，通过传感器上游测量单元获取的前、后时刻的实测系统状态向量确定状态转移矩阵，反映 上游测量所得的系统状态矩阵的转化情况：其中，z1(k-1)和z1(k)是取自上游测量单元，分别为测量区域内前后时刻的含率测量值；
利用n0结合电导传感器采样频率，对上下游测量数据进行时间配准，如公式(3)所示：z2(k)＝z'2(k+n0) (3)
其中，z'2为下游传感器单元采集数据，z2是配准后的数据；
将由传感器上游测量单元确定的状态转移矩阵用于时间配准后的下游含水率状态更新中：其中，F2(k)为传感器下游测量单元状态转移矩阵，F1(k)为上游测量单元的状态转移矩阵；
传感器测量数据时间配准后，依据传感器上下游具体数据对F(k)进行实时的自适应更新，如公式(5)所示：
将电导传感器下游测量单元k时刻含率测量值z2(k)带入公式(6)，计算k时刻的量测噪声w(k)，将量测噪声w(k)带入公式(8)，计算流态判别系数q(k)；
z2(k)-H(k)x(kk-1)＝w(k) (6)S＝H(k)p(k|k-1)+R(k) (7)q(k)＝wT(k)S-1w(k) (8)公式(6)中x(k|k-1)是k时刻系统状态预测值，估计方式如公式(9)所示；H(k)是量测矩阵，取值为1；公式(7)中，R(k)是k时刻量测噪声的方差，可取传感器含率测量数据的标准差作为恒值代入；p(k|k-1)是传感器k-1时刻对k时刻估计的协方差，计算方式如公式(10)所示；流态判别系数q(k)在卡尔曼估计过程中针对具体对象衡量环境信息的门限值，当q过大，代表流动状况处于非水相连续的状态，此处，q的限值取0.118。
4)对传感器的含率测量值进行卡尔曼状态估计和更新：
x(k|k-1)＝F(k)x(k-1|k-1) (9)p(k|k-1)＝F(k-1)p(k-1|k-1)FT(k-1)+Q(k-1) (10)根据量测信息z(k)修正估计值：
x(k|k)＝x(k|k-1)+K(k)×[z(k)-H(k)x(k|k-1)] (11)p(k|k)＝[I-K(k)H(k)]p(k|k-1) (12)其中，K(k)是传感器系统在k时刻的卡尔曼增益：
K(k)＝p(k|k-1)HT(k)×[H(k)p(k|k-1)HT(k)+R(k)]-1 (13)其中，p(k-1|k-1)是k-1时刻的协方差，其迭代初值取含水率测量值的方差；z(k)是传感器在k时刻对含水率的测量值，x(k|k-1)是传感器在k时刻对状态x(k)的预测值，x(k|k)是传感器k时刻对状态x(k)的最小均方误差估计，p(k|k)是对k时刻估计x(k|k)的协方差，Q(k-1)为k-1系统过程噪声的方差，可取恒值为0.01，I为单位矩阵，其维数等同于状态向量。
油水两相流含水率自适应估计方法
技术领域
[0001] 本发明属于流体测量技术领域，利用电导传感器与自适应卡尔曼算法测量含水率信息，用于提高两相流含水率估计精度。技术背景
[0002] 两相流广泛存在于日常生活与工业生产过程中，如食品加工，生物工程，化工产业，冶金工业和石油产业等行业。与一般单相流相比，两相流的流动状态更为复杂，建模更加困难，因此两相流的在线过程参数检测一直是学术界和工业界关注的难点和重点。特别在石油工业中，油水两相流含水率对油井监测、生产预测和开采计划制定具有重要意义。
[0003] 目前两相流含水率的测量手段包括快关阀测量方法、射线方法、微波方法、Wire-mesh方法、超声方法和电阻抗方法等。其中，电学方法具有测量速度快，敏感性强，传感器结构简单，易于安装，成本低廉等优点。两相流是一种随机的动态流动过程，流动状态复杂多变，传感器的测量信号有随机误差，且只能反映系统的外部特性。为获得多相流含水率的准确信息，需要根据实测数据对流动状态进行估计。状态估计对认知和分析测量系统具有重要的应用价值，常用的估计理论包括最小二乘估计、线性最小方差估计和卡尔曼估计等。不同的估计方法性能指标不同，使用的信息不同，其适用范围也不同。最小二乘估计通常应用于随机向量和常值向量的估计，算法结构简单但精度有限；线性最小方差估计以均方误差最小作为优化参数，适用于平稳过程，其测量信息越丰富，估计精度越高，计算量越大，而且会随着估计过程的累积，造成计算负担较重，不适用于复杂的长时间的估计过程。卡尔曼估计采用递推的形式，通过一个状态模型和一个观测模型的迭代，利用上一时刻的估计值和这一时刻的观测值，对这一刻的状态进行估计。
发明内容
[0004] 本发明提供一种油水两相流含水率自适应估计方法，该方法利用阵列式电导传感器上下游电极测量数据，实时更新卡尔曼估计迭代过程的相关参数，进行自适应状态转换，提高含水率的估计精度。本发明的技术方法如下：
[0005] 一种油水两相流含水率自适应估计方法，利用阵列式电导传感器上下游电极测量数据进行估计，包括如下步骤：
[0006] 1)将来自电导传感器上、下游的测量数据进行时间配准：
[0007] 将取自上游和下游测量单元的含率测量值z1(n)和z2(n)看作来自各态历经的平稳随机过程的两个样本函数，其中n＝1,2,3,…N为测量时间内的数据点，两传感器间距为L，设渡越时间为τ0，对应离散系统中上、下游测量数据平移组数n0；计算z1(n)和z2(n)在离散条件下的互相关函数 公式如下：
[0008]
[0009] 的最大值出现在n0处，即互相关函数峰值对应流体从上游测量单元到下游测量单元的渡越时间，通过传感器上游测量单元获取的前、后时刻的实测系统状态向量确定状态转移矩阵，反应上游测量所得的系统状态矩阵的转化情况：
[0010]
[0011] 其中，z1(k-1)和z1(k)是取自上游测量单元，分别为测量区域内前后时刻的含率测量值；
[0012] 利用n0结合电导传感器采样频率，对上下游测量数据进行时间配准，如公式(3)所示：
[0013] z2(k)＝z'2(k+n0) (3)
[0014] 其中，z'2为下游传感器单元采集数据，z2是配准后的数据；
[0015] 2)基于自适应更新的状态估计：
[0016] 将由传感器上游测量单元确定的状态转移矩阵用于时间配准后的下游含水率状态更新中：
[0017]
[0018] 其中，F2(k)为传感器下游测量单元状态转移矩阵，F1(k)为上游测量单元的状态转移矩阵；
[0019] 传感器测量数据时间配准后，依据传感器上下游具体数据对F(k)进行实时的自适应更新，如公式(5)所示：
[0020]
[0021] 3)判定两相流是否处于水连续的状态：
[0022] 将电导传感器下游测量单元k时刻含率测量值z2(k)带入公式(6)，计算k时刻的量测噪声w(k)，将量测噪声w(k)带入公式(8)，计算流态判别系数q(k)；
[0023] z2(k)-H(k)x(k|k-1)＝w(k) (6)
[0024] S＝H(k)p(k|k-1)+R(k) (7)
[0025] q(k)＝wT(k)S-1w(k) (8)
[0026] 公式(6)中x(k|k-1)是k时刻系统状态预测值，估计方式如公式(9)所示；H(k)是量测矩阵，取值为1；公式(7)中，R(k)是k时刻量测噪声的方差，可取传感器含率测量数据的标准差作为恒值代入；p(k|k-1)是传感器k-1时刻对k时刻估计的协方差，计算方式如公式(10)所示；流态判别系数q(k)在卡尔曼估计过程中针对具体对象衡量环境信息的门限值，当q过大，代表流动状况处于非水相连续的状态，此处，q的限值取0.118。
[0027] 4)对传感器的含率测量值进行卡尔曼状态估计和更新：
[0028] x(k|k-1)＝F(k)x(k-1|k-1) (9)
[0029] p(k|k-1)＝F(k-1)p(k-1|k-1)FT(k-1)+Q(k-1) (10)
[0030] 根据量测信息z(k)修正估计值：
[0031] x(k|k)＝x(k|k-1)+K(k)×[z(k)-H(k)x(k|k-1)] (11)
[0032] p(k|k)＝[I-K(k)H(k)]p(k|k-1) (12)
[0033] 其中，K(k)是传感器系统在k时刻的卡尔曼增益：
[0034] K(k)＝p(k|k-1)HT(k)×[H(k)p(k|k-1)HT(k)+R(k)]-1 (13)
[0035] 其中，p(k-1|k-1)是k-1时刻的协方差，其迭代初值取含水率测量值的方差；z(k)是传感器在k时刻对含水率的测量值，x(k|k-1)是传感器在k时刻对状态x(k)的预测值，x(k|k)是传感器k时刻对状态x(k)的最小均方误差估计，p(k|k)是对k时刻估计x(k|k)的协方差，Q(k-1)为k-1系统过程噪声的方差，可取恒值为0.01，I为单位矩阵，其维数等同于状态向量。
[0036] 本发明的实质性特点是：针对多相流流动过程，采用电导传感器上下游测量所得数据，在测量精度难以提高的情况下，利用可以对被测对象估计过程中的状态转移矩阵进行实时更新和变化的特点，改进卡尔曼估计算法的更新和迭代过程，进行自适应卡尔曼估计。利用“不确定度补偿”提高估计精度。本发明的有益效果及优点如下：
[0037] 1、该方法结合传感器具体结构，充分利用传感器上游和下游测量数据；
[0038] 2、该方法相对于传感器直接测量，提高了估计精度；
[0039] 3、该方法在计算量较小，可以实现实时在线估计。
附图说明
[0040] 以下附图描述了本发明所选择的实施例，均为示例性附图而非穷举或限制性，以电导传感器测量油水两相流为例，其中：
[0041] 图1本发明的测量方法电导传感器测量结构示意图，(a)为正视图，(b)为侧视图，(c)为立体图；
[0042] 图2本发明的在线自适应卡尔曼估计算法计算步骤。
具体实施方式
[0043] 当信号频率在10kHz到100kHz之间时，油水两相流的负载特性体现为电阻性，即电极对之间的电压与分布其间的流体阻抗成正比。在油水两相流流体可以形成回路时，油水两相流的阻抗与相分布有对应关系，因此可以通过测量电极对之间的电压来测量其相含率。
[0044] 下面结合说明书附图说明本发明的计算方法。
[0045] 图1为传感器测量空间示意图。电导传感器实物图如图1所示，详细流程如图2所示，由六个沿轴向排列的环形电极组成，其中1、6为激励电极，2、3、4、5为测量电极。采用电流激励电压测量的方式，通过向电极1中注入恒定电流并将电极6接地，可在电极1和6中间形成稳定的电学测量场，分别测量电极对2-3、电极对3-4以及电极对4-5之间的电势差V可获得其间的两相流平均电阻率。
[0046] 电导传感器为6环结构电导环传感器阵列以轴向排列的6个环状金属电极嵌入管路中的形式构成，管径D为50mm，相邻2个环状电极的距离分别为：1D、0.4D、1.2D、0.4D、1D。在此，依照图2所示流动方向，选取电极2和电极3之间的区域为上游测量单元，电极4和电极
5之间的区域为下游测量单元，两个测量单元完全相同。两侧环形电极为激励电极，连接峰峰值2mA，20kHz交流方波电流信号，中间电极为测量电极，采集电势差。管道中流体含水率不同时，流体的混合电导率会发生变化，因此每一组电势差反映的是该两环形电极之间的测量区域内部的介质分布。但是，当管道内连续相不导电，且为分散流时，电导传感器无法形成电流回路，测量失效。测量实验过程中，定义无量纲电导，
[0047]
[0048] 其中，V0是管道充满水时的电势差测量值，Vi是管道内流过两相流时的电势差测量值。水平管道油水两相流的流型主要分为分层流和分散流。对于分层流，V*即为相含率。对于分散流如油包水或水包油的混合流态，根据Maxwell理论推导，含水率与无量纲电导有如下关系：
[0049]
[0050] 其中，hw代表含水率，作为估计过程中的量测向量z。
[0051] 将来自电导传感器上、下游的测量数据进行时间配准，具体方式是将来自上游、下游的含水率数据进行离散形式下的互相关计算：
[0052]
[0053] 上游和下游测量单元的含率测量值分别为z1(n)和z2(n)，其中n＝1,2,3,…N为测量时间内的数据点，z1(n)和z2(n)的互相关函数为Rz1z2(n)。Rz1z2(n)的最大值出现在n0处，对应离散系统中上、下游测量数据平移组数n0。
[0054] 利用n0结合电导传感器采样频率，对上下游测量数据进行时间配准，使配准后下游测量单元含水率所代表的系统状态向量与上游系统状态向量一一对应：
[0055] z2(k)＝z'2(k+n0) (18)
[0056] 数据配准后，判定两相流是否处于水连续的状态。
[0057] 将电导传感器下游测量单元k时刻含率测量值z2(k)带入公式子(19)，计算k时刻的量测噪声w(k)，将量测噪声w(k)带入公式(21)，计算流态判别系数q(k)；
[0058] z2(k)-H(k)x(k|k-1)＝w(k) (19)
[0059] S＝H(k)p(k|k-1)+R(k) (20)
[0060] q(k)＝wT(k)S-1w(k) (21)
[0061] 公式(19)中x(k|k-1)是k时刻系统状态预测值，H(k)取值为1；公式(20)中，R(k)是k时刻量测噪声的方差，可取传感器含率测量数据的标准差作为恒值代入；p(k|k-1)是传感器k-1时刻对k时刻估计的协方差；流态判别系数q(k)在卡尔曼估计过程中针对具体对象衡量环境信息的门限值；将所有q(k)取均值作为判别系数q，当q过大，代表流动状况处于非水相连续的状态，不适用本发明的内容，本发明中q的限值取0.118。
[0062] 图2为本例中，在自适应卡尔曼估计方法的计算步骤。
[0063] 1)本发明在估计过程中传感器系统的量测矩阵H(k)取值为1。
[0064] 2)由传感器上游测量单元前后时刻的含水率测量值获得某一时刻的状态转移矩阵。
[0065] 传感器测量数据时间配准后，依据传感器上下游具体数据对状态转移矩阵F(k)进行实时的自适应更新：
[0066]
[0067] 3)上游测量单元和下游测量单元同时开始数据采集，而上游测量单元和下游测量单元相同，所以在忽略相间速度滑脱的前提下，可以将下游测量单元所测得含水率数值视作上游流态在时间和空间上的平移。利用油、水相混合流速结合上游和下游测量单元之间的距离，计算得出状态平移的实验数据采集组数，进而进行上游和下游测量单元数据在时间上的配准。
[0068] 4)在整个过程中量测矩阵H(k)始终取1的前提下，F(k)和F(k-1)的数值实时更新，将此F(k)和F(k-1)带入下游测量单元的卡尔曼估计和更新迭代过程中，以上游测量数据前后时刻含水率状态的转移情况影响下游测量数据的状态估计和更新过程，其状态估计过程如式(23)所示，此时，状态估计过程如下：
[0069]
[0070] 其中，x2(k|k-1)表示下游状态的预测，x1(k-1)和x1(k)分别表示前后时刻上游的实际状态。
[0071] 且下游估计的协方差可表示为：
[0072]
[0073] 其中，Q(k-1)是k-1时刻过程噪声的方差。此后，按照如下基本的形式进行卡尔曼估计：
[0074] x(k|k-1)＝F(k)x(k-1|k-1) (25)
[0075] p(k|k-1)＝F(k-1)p(k-1|k-1)FT(k-1)+Q(k-1) (26)
[0076] x(k|k)＝x(k|k-1)+K(k)×[z(k)-H(k)x(k|k-1)] (27)
[0077] p(k|k)＝[I-K(k)H(k)]p(k|k-1) (28)
[0078] K(k)＝p(k|k-1)HT(k)×[H(k)p(k|k-1)HT(k)+R(k)]-1 (29)
[0079] 5)在完成每组含水率条件下所有实验数据的迭代和滤波之后，对所有时刻的自适应估计值进行平均，作为该组实验最终的含率估计值。
