基于UKF的海上天然气井数据驱动软测量方法

王 丹 康 琦 杨居衡 宫 敬 张 奇

1. 中国石油大学（北京）经济管理学院 2. 中国石油经济技术研究院 3. 深圳清华大学研究院

4. 中国石油国际事业有限公司 5.油气管道输送安全国家工程实验室·中国石油大学（北京）

随着海洋天然气资源的开发，海上天然气生产系统的运行监测面临着严峻挑战：多相流量计安装和维护成本高，用于系统运行分析的单井流量往往无法直接测量，只能得到与之有函数关系的温度或压力等变量的观测数据，而观测数据也受到各种误差的干扰。以上因素导致生产系统运行数据获取难度大、成本高、可靠性低，制约着海上气田的安全生产[1-3]。为了解决工程中的此类问题，研究者们利用辨识建模法开发了多种计算速度快、在线适用性强的数据驱动模型，包括多项式模型和机器学习模型等[4-12]，表征系统状态变量间的函数关系，用于工业系统关键状态变量的软测量，作为实体仪表的参考或备用。由于数据驱动模型的估计准确度有限，部分研究者基于最优估计数据融合理论[13-14]，将滤波器引入软测量，利用观测数据实时校正模型模拟结果，充分利用工业系统多方面信息，得到更加贴近实际工况的估计结果。

滤波器是普遍应用的、以最优估计为基础的动态数据融合方法，计算速度快，适用于动态系统的实时估计。滤波器自1960年产生以来，从仅适用于线性系统的卡尔曼滤波器[15]发展出了扩展卡尔曼滤波器[16]、集合卡尔曼滤波器[17]、无迹卡尔曼滤波器（Unscented Kalman Filter, UKF）[18-20]等一系列适用于非线性系统的滤波器[13,21]。多年来，滤波应用范围也从航天导航系统发展到电子、油气生产及运输等工业系统。Aguirre[22]等将滤波器与电子振荡器的多项式数据驱动模型结合，用于改进电子振荡器的模拟信号。才建[23]、王永红[24]、Lorentzen[25-27]等将滤波器应用于输油和气—液两相流管道流量和压力的估计，以及采油过程中井筒压力和流量的估计。Teixeira[4]等将滤波器与多项式模型、多层感知器模型组成的数据驱动模型库相结合，提升井底压力软测量效果；
并利用模型冗余建立融合滤波器，进一步改善井底压力估计结果。在油气生产的数据驱动软测量领域，目前研究者大多仅依赖数据驱动模型进行软测量，对滤波器与数据驱动模型的结合，以及融合滤波器的应用仍然较少，可见，软测量准确度研究待进一步提高。

针对天然气生产系统，引入UKF创新建立了数据驱动的软测量方案，用于改善单井流量和井口压力的估计准确度。①建立外源输入非线性自回归（Nonlinear Auto-Regressive with Exogenous Inputs，NARX）模型库；
②利用每个NARX模型构建单个UKF；
③利用NARX模型的冗余，创新建立两种融合UKF。通过实际生产数据，验证了滤波技术与数据驱动模型结合的优势，比较了单个UKF和融合UKF的准确性与计算成本，推荐适合工程应用的最佳软测量方案，从而提高流动管理系统的软测量功能模块，提高其获取可靠运行数据的能力，以期推动智慧油气田建设。

1.1 系统介绍和算法框架

以某海上天然气两井生产系统（图1，两井编号分别为Ⅰ、Ⅱ）为例建立基于融合滤波的数据驱动软测量算法。该两井生产系统位于某水下气田，于2013年12月底投产。两支路生产的天然气—凝析液经管汇进入海管，连入配有分离器的处理平台。系统沿线分布多个配有压力和温度传感器的观测点，包括井底、井口、油嘴下游和平台分离器，用数字1、2、3、4分别代表以上位置。单井未设置多相流量计，因此，单井流量不可观测。

图1 某海上天然气两井生产系统拓扑结构图

结合系统数据驱动模型和仪表测量值两方面信息，兼顾模型和测量值的不确定度，基于数据融合理论，引入滤波技术建立“三步”数据驱动动态软测量算法，用于单井流量和井口压力的在线校正估算。所建软测量算法是过程辨识和状态估计技术的结合（图2），其中，重构井筒状态空间方程和构造合适的滤波器是两个关键点。对于前者，对各井筒建立NARX模型库用以重构井筒状态空间方程。对于后者，引入基于点集传播的非线性滤波器以适应非线性天然气生产系统。

图2 基于滤波器的数据驱动动态估计算法框架图

1.2 NARX模型库建立

建立NARX模型用以捕捉天然气—凝析液在井筒中流动的非线性动态特性，描述当期输出与先行输出、先行输入之间的动态约束关系[28-29]，为构建状态空间方程和滤波器奠定基础。离散时间下NARX模型的通式表示为：

式中g表示NARX模型；
R表示实数空间；
RN×ζ表示N行ζ列的实数矩阵，其他同理；
z(k)、u(k)、e(k)分别表示某时步下的输出向量、输入向量和误差向量；
z、u、e分别表示各时步下输入向量、输出向量、误差向量组成的矩阵；
ζ表示输出变量的个数；
μ表示输入变量的个数；
ρ=[ρ1, …,ρζ]表示输出时滞向量；
τ=[τ1, …,τμ]表示输入时滞向量，每个输出（输入）变量的最大时滞不同Θ表示参数向量；
k表示离散时步（k=1, 2, …,N，其中N表示采样点个数）。

针对每个井筒，以单井流量和井口压力为输出变量（ζ=2），以系统沿线的可观测变量为输入变量，基于时间序列分析[30-31]，综合动态和稳态训练样本（图1），采用正交最小二乘回归[28-29]和深度学习[9,10]等黑箱辨识技术，以及双目标最小二乘[5,32]和迁移学习等参数校正技术，建立了一系列相互独立的黑箱和灰箱NARX模型（gj，j=1, 2, …,m，其中m为模型个数），组成NARX模型库。模型库中包括：以多项式（Polynomial，POLY）为基础的POLY-NARX黑箱动态模型（g1）、POLY-NARX灰箱动态模型（g2），以深度前馈网络（Dense Neural Network，DNN）为基础的DNN-NARX黑箱动态模型（g3、g4），DNNNARX灰箱动态模型（g5、g6）。其中，g3和g5以支路和平台可观测变量为输入变量，g4和g6仅以支路可观测变量为输入变量。以上NARX模型相关建模过程和模型参数设置采用Wang等[33]研究成果。

1.3 UKF理论基础

对于非线性的天然气生产系统，经典卡尔曼滤波器（Kalman Filter，KF）很多条件和性质不再成立[13,21]。1990年后出现的Sigma点卡尔曼滤波器（Sigma Point Kalman Filter，SPKF）[21]适用于非线性系统。SPKF通过确定性采样法得到Sigma点集，用Sigma点集非线性映射后的期望和方差近似状态向量（X）概率分布。SPKF具有以下优势：①逼近概率分布比逼近任意非线性函数更容易实现，避免对非线性函数进行解析求导；
②对均值和方差的逼近精度提高到二阶及以上；
③计算量较小，滤波稳定性较高[13]。

Julier和Uhlmann于1994年提出的UKF是典型的SPKF方法，利用历史状态向量的后验分布，采用无迹变换法（Unscented Transform，UT）生成规模为2nX+1（nX为X维度）、对称分布的Sigma点集[20-21]，并以此进行非线性传播。

1.3.1 状态空间方程和滤波步骤

系统的非线性状态空间方程是构建UKF的基础，由非线性状态转移方程和观测方程组成[式（2）]，前者表示不同时步状态向量间的递推关系，后者表示状态向量和观测向量间的关系。

式中X(k)、U(k－1)、Y(k)表示某时步下的状态向量、输入向量、观测向量；
nX、nU、nY分别表示状态变量、输入变量和观测变量的个数；
C(k)表示某时步下的观测矩阵；
W(k－1)、V(k)分别表示某时步下的系统模型噪声向量、观测噪声向量，假设是均值为0的白噪。

滤波是一种基于最小化均方误差[13]的“两步”递推估计算法。第一步是一步预测，也叫时间更新，根据状态向量前一时步的分布信息Px(k－1) }得到其一步预测值，即当前时步的先验估计。第二步是测量更新，以滤波增益为权重，融合观测向量的测量新息与状态向量的先验估计，得到状态向量当前时步的后验估计。确定状态向量的初始分布Px(0)}后，按以上两步不断进行递推估算，可得到每个时步状态向量的后验估计。用式（3）概括UKF的滤波步骤，具体计算过程详见文献[21]。

式中Forecast表示一步预测阶段；
DataAssimilation表示数据融合，也就是测量校正阶段；
Q(k－1)、R(k)分别表示某时步下系统模型噪声和观测噪声的协方差阵，是系统模型和观测数据的不确定度；
(k)表示某时步下观测向量的测量值；
分别表示Sigma点集、观测向量点集的期望；
Px(k)表示状态向量的误差协方差阵；
Py(k)表示观测向量的误差协方差阵；
Pxy(k)表示状态向量和观测向量的误差协方差阵；
S(k)表示滤波增益，S(k)=Pxy(k)P-y1(k)；
某时步下Y(k)的新息为，即观测向量的估计残差；
k|k－1代表k－1时步到k时步的一步预测值或先验估计。

1.3.2 噪声统计特性设置

滤波算法是利用状态空间方程和白噪声激励的统计特性形成。在白噪激励的假设下，Q(k)=diag[qi2(k)]qi2(k)表示系统模型对每个状态变量一步预测的误差方差，通过模型模拟结果得到。ro2(k)表示每个观测变量的噪声方差，通过仪表的不确定度得到。

1.4 基于融合UKF的数据驱动动态估计

基于所建NARX模型库重构井筒状态空间方程，以此构建单个UKF。进而利用NARX模型的冗余，建立融合UKF。

1.4.1 天然气生产井状态空间方程重构

以天然气生产井为估计单元，对每个井筒构建基于NARX模型的状态空间方程。

1.4.1.1 状态转移方程

状态转移方程的一步预测算子（f）通过时滞的NARX模型（gj）变形得到[4,21]，gj的输入向量、输出向量与f的输入向量、状态向量之间具有明确的对应关系[式（4）]。k时步下，U(k)对应u(k)张成的空间，X(k)对应z(k)张成的空间。z=[G, p2]，则X=[G, p2]；
不同状态转移方程对应的gj不同，相应的输入向量也不同。

式中wi(i=1, 2, …, nX)表示系统模型对每个状态变量一步预测的误差。

1.4.1.2 观测方程

为便于驱动滤波器，取状态向量中的可观测变量组成观测向量[22,25]，因此Y=[z2]=[p2]。Y(k)=[z2(k), …, z2(k－ ρ2+1)]= [p2(k), …, p2(k－ ρ2+1)]，C(k)=[0，I]，其中0是全零列，I是单位阵。

1.4.1.3 初始分布与噪声协方差阵

取各状态变量样本集中第1个时刻的值作为状态向量的初始期望¯(0)。状态向量的初始误差协方差阵 Px(0)=diag[σi2(0) ]，σi2(0)=10－5。

Qk中，qi2(k)根据各模型对输出变量模拟误差的方差赋值[4]。R(k)中ro2(k)=0.003 6。

1.4.2 融合UKF建立

为发挥多模型优势，提高软测量的估计性能，利用所建NARX模型库中模型的冗余，构造了两种融合UKF。第一种融合滤波器“先融合后滤波”，简称为FF1（图3）；
第二种融合滤波器“先滤波后融合”[4]，简称为FF2（图3）。

图3 融合滤波器FF1和FF2的计算流程图

1.4.2.1 建立FF1

针对每个井筒构建1个UKF，该UKF的状态转移方程是模型库中多个模型的组合。模型库中择优选择多个模型并行计算，对所得多个模型模拟值进行加权融合得到一步预测值，在此基础上进行测量校正得到后验估计值。各时步下每个模型的相应权重根据各自的模拟误差赋值，误差较小的模型对应较大的权重。

1.4.2.2 建立FF2

每个模型分别构建1个UKF，通过多个模型得到多个状态空间方程互不相同的独立UKF，组成滤波库再进行融合。择优选择多个独立UKF并行计算，对所得多个UKF的后验估计结果进行加权融合，得到融合滤波值。各时步下每个滤波器的相应权重根据各自的新息赋值，新息较小的滤波器对应较大的权重。

以上述海上天然气两井生产系统（图1）中Ⅱ井的生产数据为例，从可行性、准确度和计算速度等方面对所建数据驱动估计算法进行验证。根据前期研究结果，灰箱模型的全局准确性和稳定性高于相应黑箱模型[33]，因此，直接基于灰箱模型g2、g5、g6构建3个相互独立的UKF（表1），分别命名为UKF2、UKF5、UKF6。对以上灰箱模型进一步择优，用于建立融合滤波器FF1和FF2。在此基础上比较单个UKF和FF1、FF2的估计结果与计算时间。

表1 Ⅱ井灰箱模型g2、g5、g6和滤波器UKF2、UKF5、UKF6设置表

利用生产系统的常规工艺参数和生产数据进行算法测试。常规工艺参数包括井身结构和流体组分，生产数据由集散控制系统采集的仪表观测数据和两相流模拟软件OLGA的模拟数据组成，共涵盖了2014年1—4月的5 000组动态数据，并从中提取出10组稳态生产数据（图4）。其中，仪表观测数据包括T1、p1、T2、p2、T3、p3、T4、p4、V；
OLGA 的 模拟数据主要包括G。以2014年1—4月p2的仪表测量值和G的OLGA模拟值作为算法估算结果的对照，为评判估算结果的准确性提供依据。受限于现场生产条件，各生产井不配备流量计，因此目前难以通过现场生产数据来验证OLGA计算结果。OLGA是工业界公认的多相流模拟专业软件，现有条件下相对可行的手段是利用OLGA的流量模拟结果来验证算法的流量估算结果。

图4 Ⅱ井2014年1—4月动态和稳态生产数据图

使用均方根误差（Root Mean Square Error，RMSE）、平均绝对误差（Mean Absolute Error，MAE）、均方误差（Mean Square Error，MSE）及决定系数R2作为算法准确度的评价指标[34]。

2.1 不同NARX模型所建UKF对比

将UKF2、UKF5、UKF6对Ⅱ井流量（G）和井口压力（p2）的估算结果（图5）、估算准确度指标（表2）与相应NARX模型g2、g5、g6的结果一一对比。与g2、g5、g6的模拟结果相比，具有测量校正功能的UKF2、UKF5、UKF6对G和p2的估算准确度明显提高，R2较大且RMSE、MRE、MAE较小。其中，UKF5和UKF6的R2达到0.9以上，而g5、g6的R2约为0.7、0.5；
与g6相比，UKF6的RMSE、MRE、MAE减小50%以上。以上结果说明，通过滤波器反馈井口压力测量值能够有效校正模型对单井流量和井口压力的模拟结果。

表2 Ⅱ井模型g2、g5、g6和滤波器UKF2、UKF5、UKF6的计算准确度指标统计表

图5 Ⅱ井模型g2、g5、g6和滤波器UKF2、UKF5、UKF6对单井流量和井口压力的计算结果图

对不同模型构建的UKF2、UKF5、UKF6进行对比可见，通过DNN-NARX灰箱模型构建的UKF5、UKF6的估算准确度高于POLY-NARX灰箱模型构建的UKF2。进一步比较UKF5与UKF6发现，UKF5的输入向量中包含支路变量和平台变量，而UKF6的输入向量仅包含支路变量。虽然UKF6的输入变量个数较少，但与UKF5的估算误差相差不大。因此，DNN-NARX模型构造的UKF具有较高的准确性和一定的稳定性，在工程实际中，如果部分观测点的仪表缺失导致输入变量个数减少，仍能保持较高的估算准确度。因此，选择DNN-NARX灰箱模型作为融合滤波器的首选模型。

2.2 单个UKF与融合UKF对比

基于DNN-NARX灰箱模型g5、g6构建融合滤波器FF1、FF2，比较单个滤波器UKF5、UKF6和融合滤波器FF1、FF2对G和p2的估算准确度（图6、表3）与计算速度。

表3 单个滤波器UKF5、UKF6和融合滤波器FF1、FF2的计算准确度指标统计表

图6 Ⅱ井单个滤波器UKF5、UKF6与融合滤波器FF1、FF2对单井流量和井口压力的计算结果图

根据结果可知，两个融合滤波器中，FF1的R2较高且RMSE、MAE、MRE等估算误差明显较小。FF1与单个滤波器中估算效果较好的UKF5相比，前者R2较高，且MAE和MRE较小。而FF2未体现绝对优势，其估算准确度介于UKF5和UKF6之间。因此，融合滤波器FF1的准确度最高。

对单个UKF和融合UKF的计算时间进行对比。计算程序使用Java语言进行编写，在Windows 7 64-bit操作系统平台上搭建JDK开发环境，测试环境为4核PC，处理器性能为Intel（R）Core（TM）CPU i5-5200 @ 2.20 GHz。UKF5、UKF6、FF1、FF2 的单步计算时间分别为8 ms、8 ms、10 ms、15 ms。可以看出，FF1和FF2比单个UKF计算耗时，但均能满足工程中在线估计的需求。同时发现，FF2比FF1耗时长，这是因为，FF1是单个滤波器运行，而FF2是多滤波器同时并行，且每个滤波器涉及Sigma点集的并行计算，计算量较大。可以推断，随着FF1和FF2所含模型个数增加，与单个UKF计算时间的差距会加大。

综上所述，工程应用中首先推荐FF1作为数据驱动软测量方案。FF1能够综合不同模型的优势，整体估算准确度最高，计算量和耗时较少。其次推荐DNN-NARX灰箱模型构造的单个UKF作为数据驱动软测量方案，该方案也具有较高的估算准确度和较快的计算速度。

1）基于UKF软测量方案的准确度高于纯模型软测量方案的准确度，所建软测量算法准确度较高，计算量小，拓展性和在线适用性强，能够提高流动管理系统获取可靠运行数据的能力。

2）所建包含POLY-NARX和DNN-NARX的模型库中，DNN-NARX灰箱模型构造的UKF准确性和稳定性较高。

3）单个UKF与FF1、FF2相比，“先融合后滤波”的FF1估计准确度较高，计算速度与单个UKF相当，推荐作为首选软测量方案。

猜你喜欢 准确度滤波器滤波 基于HP滤波与ARIMA-GARCH模型的柱塞泵泄漏量预测农业工程学报(2022年10期)2022-08-22基于改进自适应中值滤波的图像降噪方法*舰船电子工程(2022年6期)2022-08-02浅谈有源滤波器分析及仿真客联(2022年4期)2022-07-06基于多模谐振器的超宽带滤波器设计成都信息工程大学学报(2022年2期)2022-06-14影响重力式自动装料衡器准确度的因素分析科海故事博览·下旬刊(2022年4期)2022-05-07基于非下采样剪切波变换与引导滤波结合的遥感图像增强计算机应用(2016年10期)2017-05-12论提高装备故障预测准确度的方法途径价值工程(2016年32期)2016-12-20Word中“邮件合并”功能及应用考试周刊(2016年34期)2016-05-28FIR滤波器线性相位特性的研究智能计算机与应用(2016年1期)2016-03-02对GB 17167实施过程中衡器准确度要求问题的探讨中小企业管理与科技·中旬刊(2014年10期)2015-02-03