天然气井筒生产运行工况动态估计方法及装置

1.本发明涉及流体计算领域，尤指一种天然气井筒生产运行工况动态估计方法及装置。


背景技术：

2.随着海洋油气资源的开发，基于气-液两相流混输技术的海上天然气生产系统的运行监测面临以下难题。第一，高昂的安装和维护成本使得多相流量计安装不足，单井产量不可观测。第二，对于温度、压力等可观测变量，其传感器采集的实时生产数据会受到显著误差和随机误差的干扰，使测量数据违背基本控制方程；这些误差通常来自于仪表本身的不确定性、洋流的扰动以及人工校准的失误等等。第三，海底环境恶劣，仪表面临更高的失效风险，尤其是井底的温度、压力传感器。因此，很难直接获取有效的、可信度高的生产运行信息。
3.基于数据驱动的天然气井筒生产运行工况动态估计方法就是针对上述问题所提出的一种新方法。结合生产系统沿线观测点的温度、压力和变量间的函数关系，估计不可测的流量和沿线节点压力，并修正观测值，以获取准确的运行信息。作为一种软测量技术对于油气生产具有重要意义，基于数据驱动的天然气井筒生产运行工况动态估计方法，不仅能作为一种状态监测工具，降低测量成本，便于进行在线安全分析；更契合流动安全保障和数字油气田等研究热点，为开发流动管理系统和建设数字油气田准备高质量的数据库。


技术实现要素：

4.本发明目的在于提供一种天然气井筒生产运行工况动态估计方法及装置，以帮助工程人员及时掌握难于通过实体设备监测的井筒生产运行工况。
5.为达上述目的，本发明提供了一种天然气井筒生产运行工况动态估计方法，所述方法包含：获取生产数据，根据所述生产数据利用学习算法构建的辨识技术对天然气井筒进行灰箱系统辨识，建立非线性自回归灰箱的动态辨识模型库，并构建井筒状态空间方程；根据所述动态辨识模型库中的辨识模型构建状态转移方程，并通过观测矩阵构建观测方程，通过所述状态转移方程和所述观测方程重构所述井筒状态空间方程；通过重构后的井筒状态空间方程结合无迹卡尔曼滤波-集合卡尔曼滤波的数据融合技术，融合所述辨识模型的估计值与井筒观测点的测量值；根据融合后的所述估计值和所述测量值融合所述无迹卡尔曼滤波器-集合卡尔曼滤波器的估计结果，以过程辨识和状态估计相结合的方式对天然气井筒生产运行工况进行动态估计，获得目标估计值。
6.在上述天然气井筒生产运行工况动态估计方法中，优选的，根据所述生产数据利用学习算法构建的辨识技术对天然气井筒进行灰箱系统辨识，建立非线性自回归灰箱的动态辨识模型库包含：根据所述生产数据，在未引入各层规范化和丢弃技术的基础上利用传统梯度下降法和sigmoid激活函数对神经网络模型进行有监督预训练，得到传统神经网络模型；以泄漏修正线性单元作为激活函数，根据所述生产数据通过随机梯度下降法配合各
层规范化和丢弃技术对预设神经网络模型进行有监督预训练，得到深度前馈网络模型；以修正线性单元作为激活函数，加入卷积层、池化层、全连接层后，根据所述生产数据通过随机梯度下降法对神经网络模型进行有监督预训练，得到卷积神经网络模型；根据传统神经网络模型、深度前馈网络模型和卷积神经网络模型获得描述井筒非线性动态特性的所述动态辨识模型库。
7.在上述天然气井筒生产运行工况动态估计方法中，优选的，建立非线性自回归灰箱的动态辨识模型库还包含：根据稳态样本，通过迁移学习调整预训练所得辨识模型。
8.在上述天然气井筒生产运行工况动态估计方法中，优选的，根据融合后的所述估计值和所述测量值融合所述无迹卡尔曼滤波器-集合卡尔曼滤波器的估计结果，以过程辨识和状态估计相结合的方式对天然气井筒生产运行工况进行动态估计包含：采用无迹卡尔曼滤波-集合卡尔曼滤波分别构建独立滤波器；于所述独立滤波器中以动态辨识模型作为状态转移算子，以状态量和观测量的关系矩阵作为观测算子，通过递推估计法融合所述辨识模型的估计值与井筒观测点的测量值；根据观测点的测量值在线校正模型估计值，获得每个时步状态量的后验估计值；根据所述后验估计值对无迹卡尔曼滤波器-集合卡尔曼滤波器并行计算得到的后验估计结果，利用预设权重进行加权融合，得到融合滤波值。
9.在上述天然气井筒生产运行工况动态估计方法中，优选的，通过递推估计法融合所述辨识模型的估计值与井筒观测点的测量值还包含：在无迹卡尔曼滤波器中，以上一时步状态量的后验估计作为当前时步状态量的初始分布；根据所述初始分布通过确定性采样法构造状态量的sigma点集；通过状态转移算子，对所述sigma点集进行非线性映射，得到sigma点集的先验估计；通过观测算子，对所述sigma点集再次进行非线性映射，得到观测点集的先验估计；通过sigma点集和观测点集的先验估计，更新状态量误差协方差阵和观测量误差协方差阵的先验估计；在所述集合卡尔曼滤波器中，以上一时步状态量的后验估计作为当前时步状态量的初始分布；根据所述初始分布通过蒙特卡洛随机抽样法构造状态量的背景集合；通过状态转移算子，对所述背景集合进行非线性映射得到背景集合的先验估计；通过观测算子对所述背景集合再次进行非线性映射，得到观测点集的先验估计；通过背景集合和观测点集的先验估计，更新状态量误差协方差阵和观测量误差协方差阵的先验估计。
10.在上述天然气井筒生产运行工况动态估计方法中，优选的，根据观测点的测量值在线校正模型估计值还包含：通过无迹卡尔曼滤波预测过程中所述误差协方差阵的先验估计，更新无迹卡尔曼滤波对应的滤波增益；通过滤波增益对预设测量值和sigma 点集的先验估计值进行加权融合，获得状态量和状态量误差协方差阵的后验估计，并作为下一时步状态量的初始分布；通过集合卡尔曼滤波预测过程中所述误差协方差阵的先验估计，更新集合卡尔曼滤波对应的滤波增益；通过滤波增益对预设测量值和背景集合的先验估计值进行加权融合，并作为下一时步状态量的初始分布。
11.本发明还提供一种天然气井筒生产运行工况动态估计装置，所述装置包含：辨识模块、重构模块和估计模块；所述辨识模块用于获取生产数据，根据所述生产数据利用学习算法构建的辨识技术对天然气井筒进行灰箱系统辨识，建立非线性自回归灰箱的动态辨识模型库，并构建井筒状态空间方程；所述重构模块用于根据所述动态辨识模型库中的辨识模型构建状态转移方程，并通过观测矩阵构建观测方程，通过所述状态转移方程和所述观
测方程重构所述井筒状态空间方程；所述估计模块用于通过重构后的井筒状态空间方程结合无迹卡尔曼滤波-集合卡尔曼滤波的数据融合技术，融合所述辨识模型的估计值与井筒观测点的测量值；根据融合后的所述估计值和所述测量值融合所述无迹卡尔曼滤波器-集合卡尔曼滤波器的估计结果，以过程辨识和状态估计相结合的方式对天然气井筒生产运行工况进行动态估计，获得目标估计值。
12.在上述天然气井筒生产运行工况动态估计装置中，优选的，所述构建模块包含：根据所述生产数据，在未引入各层规范化和丢弃技术的基础上利用传统梯度下降法和 sigmoid激活函数对神经网络模型进行有监督预训练，得到传统神经网络模型；以泄漏修正线性单元作为激活函数，根据所述生产数据通过随机梯度下降法配合各层规范化和丢弃技术对预设神经网络模型进行有监督预训练，得到深度前馈网络模型；以修正线性单元作为激活函数，加入卷积层、池化层、全连接层后，根据所述生产数据通过随机梯度下降法对神经网络模型进行有监督预训练，得到卷积神经网络模型；根据传统神经网络模型、深度前馈网络模型和卷积神经网络模型获得描述井筒非线性动态特性的所述动态辨识模型库。
13.本发明还提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述方法。
14.本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有执行上述方法的计算机程序。
15.本发明的有益技术效果在于：以天然气生产井筒为对象进行过程辨识，利用历史采样数据和软件模拟值，采用传统机器学习mlp-narx、深度学习dnn-narx和卷积神经网络cnn-narx辨识模型，建立黑箱/灰箱模型库，描述井筒状态变量之间的动态约束关系，建立天然气井筒生产运行工况状态空间方程，构造ukf和enkf，以过程辨识和状态估计相结合的方式，用于井筒流量和压力的估计。利用该计算，可以替代实体流量计和传感器工作，并满足生产要求。该技术需要生产历史数据和在线数据，特别适用于采用dcs的气田生产系统。该技术具有较强的弹性和适应性，只要具有足够的基础参数和实时的工况数据，不仅可以同时预测多口天然气井筒生产运行工况，且不需要改造生产系统即可直接应用该技术。
附图说明
16.此处所说明的附图用来提供对本发明的进一步理解，构成本技术的一部分，并不构成对本发明的限定。在附图中：
17.图1为本发明一实施例所提供的天然气井筒生产运行工况动态估计方法的流程示意图；
18.图2为本发明一实施例所提供的天然气井筒生产运行工况动态估计方法的逻辑原理图；
19.图3为本发明一实施例所提供的辨识模型的构建流程示意图；
20.图4a为本发明一实施例所提供的mlp-narx和dnn-narx模型结构示意图；
21.图4b为本发明一实施例所提供的cnn-narx模型结构示意图；
22.图5为本发明一实施例所提供的估计值计算流程示意图；
23.图6为本发明一实施例所提供的天然气井筒生产系统拓扑结构示意图；
24.图7a至图7d为本发明一实施例所提供的训练数据和测试数据示意图；
25.图8为本发明一实施例所提供的迁移学习逻辑示意图；
26.图9为本发明一实施例所提供的ukf递推估计流程示意图；
27.图10为本发明一实施例所提供的天然气井筒生产运行工况动态估计装置的结构示意图；
28.图11为本发明一实施例所提供的电子设备的结构示意图。
具体实施方式
29.以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。
30.另外，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。
31.请参考图1所示，本发明所提供的天然气井筒生产运行工况动态估计方法，具体包含：
32.s101获取生产数据，根据所述生产数据利用学习算法构建的辨识技术对天然气井筒进行灰箱系统辨识，建立非线性自回归灰箱的动态辨识模型库，并构建井筒状态空间方程；
33.s102根据所述动态辨识模型库中的辨识模型构建状态转移方程，并通过观测矩阵构建观测方程，通过所述状态转移方程和所述观测方程重构所述井筒状态空间方程；
34.s103通过重构后的井筒状态空间方程结合无迹卡尔曼滤波-集合卡尔曼滤波的数据融合技术，融合所述辨识模型的估计值与井筒观测点的测量值；
35.s104根据融合后的所述估计值和所述测量值融合所述无迹卡尔曼滤波器-集合卡尔曼滤波器的估计结果，以过程辨识和状态估计相结合的方式对天然气井筒生产运行工况进行动态估计，获得目标估计值。
36.由此，结合实际生产数据和变量间的约束关系，可实现对不可测天然气井筒流量和压力的闭环在线估计，有助于工程人员及时掌握难于通过实体设备监测的井筒生产运行工况，可节省实体流量计和传感器装置的费用，降低测量成本，便于在线安全分析，更可为油气田数字智能化提供高质量的数据资源。在实际工作中，上述实施例中主要是基于现场实际数据，通过“传统神经网络(mlp)-非线性自回归(narx) 辨识模型(mlp-narx)”、“深度前馈网络(dnn)-非线性自回归(narx)辨识模型(dnn-narx模型)”、“卷积神经网络(cnn)-非线性自回归(narx)辨识模型(cnn-narx)”构建井筒生产流量、压力等信息的状态空间，结合无迹卡尔曼滤波(ukf)和集合卡尔曼滤波(enkf)等数据融合技术，以过程辨识和状态估计相结合的方式(图2)，实现对天然气井筒生产运行工况的动态估计。具体包括如下步骤：(1)基于生产数据，利用“mlp-narx模型”、“dnn-narx模型”和“cnn-narx 模型”对天然气井筒进行系统辨识，构建井筒状态空间方程；(2)在井筒状态空间基础上，构造“ukf”和“enkf”相结合的融合滤波器，通过“一步预测”后的“测量校正”两步实现不断递推估计，可应用于天然气矿藏开采行业。
37.请参考图3所示，在本发明一实施例中，根据所述生产数据利用深度前馈网络的非线性自回归的辨识模型对天然气井筒进行系统辨识包含：
38.s301根据所述生产数据，在未引入各层规范化和丢弃技术的基础上利用传统梯度下降法和sigmoid激活函数对神经网络模型进行有监督预训练，得到传统神经网络模型；
39.s302以泄漏修正线性单元作为激活函数，根据所述生产数据通过随机梯度下降法配合各层规范化和丢弃技术对预设神经网络模型进行有监督预训练，得到深度前馈网络模型；
40.s303以修正线性单元作为激活函数，加入卷积层、池化层、全连接层后，根据所述生产数据通过随机梯度下降法对神经网络模型进行有监督预训练，得到卷积神经网络模型；
41.s304根据传统神经网络模型、深度前馈网络模型和卷积神经网络模型获得描述井筒非线性动态特性的所述动态辨识模型库。
42.其中，建立非线性自回归灰箱的动态辨识模型库还包含：根据稳态样本，通过迁移学习调整预训练所得辨识模型。
43.具体的，在实际工作中基于生产数据，利用传统机器学习、深度学习等辨识技术对天然气井筒进行系统辨识，构建包含“mlp-narx模型”、“dnn-narx模型”、“cnn-narx模型”在内的灰箱动态辨识模型库，构建井筒状态空间方程，如图4a 至图4b所示；
44.mlp-narx模型的辨识第一步：确定输入量和输出量、时滞向量、选取动态训练样本和稳态训练样本、使用传统梯度下降法和sigmoid激活函数、不使用各层规范化和“丢弃”(dropout)技术，对模型进行有监督预训练，以确定mlp-narx模型超参数(网络层数、各层节点数、学习率、动量因子)和参数(各节点权重ω和偏移量b)，得到mlp-narx黑箱动态模型，描述井筒非线性动态特性。
45.dnn-narx模型的辨识第一步：确定输入量和输出量、时滞向量、选取动态训练样本和稳态训练样本、使用小批量随机梯度下降法(mini-batch sgd)、选用“泄漏修正线性单元”(leakyrelu)作为激活函数、辅以各层规范化和“丢弃”(dropout) 技术，对模型进行有监督预训练，以确定dnn-narx模型超参数(网络层数、各层节点数、学习率、丢弃率、动量因子、批量大小等)和参数(各节点的权重ω和偏移量b)，得到dnn-narx黑箱动态模型，描述井筒非线性动态特性。
46.cnn-narx模型辨识第一步：确定输入量和输出量、确定二维输入矩阵的维度(时滞向量和输入变量个数)、选取动态训练样本和稳态训练样本、选用“修正线性单元”(relu)作为激活函数、辅以批量规范化和“丢弃”(dropout)技术，对模型进行有监督预训练，以确定cnn-narx模型超参数(卷积层和池化层的层数，卷积层的卷积核个数、卷积核扫描尺寸、移动步长，池化层的池化方法、池化层过滤器尺寸和步长，学习率、动量因子、丢弃率、批量个数等)和全连接层参数(各节点的权重ω和偏移量b)，得到cnn-narx黑箱动态模型，描述井筒非线性动态特性。
47.mlp-narx模型、dnn-narx模型和cnn-narx模型的辨识第二步：利用稳态小样本，基于迁移学习对mlp-narx、dnn-narx和cnn-narx黑箱模型分别进行微调，建立相应的mlp-narx、dnn-narx和cnn-narx灰箱动态模型，扩大模型适用范围。
48.请参考图5所示，通过重构后的井筒状态空间方程结合无迹卡尔曼滤波、集合卡尔
曼滤波，融合所述辨识模型估计值和观测点测量值，融合所述无迹卡尔曼滤波器和集合卡尔曼滤波器的估计值，以过程辨识和状态估计相结合的方式对天然气井筒生产运行工况进行动态估计获得估计值包含：
49.s501采用无迹卡尔曼滤波、集合卡尔曼滤波分别构建独立滤波器；
50.s502于所述独立滤波器中以动态辨识模型作为状态转移算子，以状态量和观测量的关系矩阵作为观测算子，通过递推估计法融合所述辨识模型的估计值与井筒观测点的测量值；
51.s503根据观测点的测量值在线校正模型估计值，获得每个时步状态量的后验估计值；
52.s504根据所述后验估计值对无迹卡尔曼滤波器-集合卡尔曼滤波器并行计算得到的后验估计结果，利用预设权重进行加权融合，得到融合滤波值。
53.在上述实施例中，所述预测处理包含：在所述ukf中，以上一时步状态量的后验估计作为当前时步状态量的初始分布，根据所述初始分布，通过确定性采样法构造状态量的sigma点集；通过状态转移算子，对所述sigma点集进行非线性映射，得到 sigma点集的先验估计；通过观测算子，对所述sigma点集再次进行非线性映射，得到观测点集的先验估计；通过sigma点集和观测点集的先验估计，更新状态量误差协方差阵和观测量误差协方差阵的先验估计。在所述enkf中，以上一时步状态量的后验估计作为当前时步状态量的初始分布，根据所述初始分布，通过蒙特卡洛随机抽样法构造状态量的背景集合；通过状态转移算子，对所述背景集合进行非线性映射，得到背景集合的先验估计；通过观测算子，对所述背景集合再次进行非线性映射，得到观测点集的先验估计；通过背景集合和观测点集的先验估计，更新状态量误差协方差阵和观测量误差协方差阵的先验估计。
54.在上述实施例中，所述测量校正处理还包含：通过ukf一步预测过程中所述误差协方差阵的先验估计，更新ukf对应的滤波增益，通过滤波增益对预设测量值和 sigma点集的先验估计值进行加权融合，获得状态量和状态量误差协方差阵的后验估计，实现所述辨识模型估计值的在线校正，并作为下一时步状态量的初始分布。通过enkf一步预测过程中所述误差协方差阵的先验估计，更新enkf对应的滤波增益，通过滤波增益对预设测量值和背景集合的先验估计值进行加权融合，对模型估计值进行在线校正，获得状态量和状态量误差协方差阵的后验估计，实现所述辨识模型估计值的在线校正，并作为下一时步状态量的初始分布。
55.具体的，在实际工作中，上述实施例原理在于在井筒状态空间基础上，构造ukf 和enkf，通过“一步预测”后的“测量校正”两步实现不断递推估计。1、在上述 mlp-narx、dnn-narx和cnn-narx模型的辨识基础上，构建状态转移方程，并通过观测矩阵构建观测方程，重构天然气井筒生产状态空间方程。每个模型分别构建一个状态转移方程，与相应观测方程结合，重构相互独立的天然气井筒生产状态空间方程。
56.2、初始化状态量的分布并按以下两步不断进行递推估计，便可得到每个时步的状态量后验估计值。ukf和enkf滤波步骤均用下式概括：
[0057][0058]
在上式中，“forcast”代表一步预测阶段，“dataassimilation”代表测量校正阶段；k是离散时步；下标“k”代表k时步的后验估计，下标“k|k-1”代表k-1时步到k时步的一步预测值或先验估计；x是状态向量，y是观测向量；k是滤波增益；是观测量的测量值；p
x
是状态量的误差协方差阵。式中，分别为sigma点集、观测量点集的期望；是k时步的新息，即观测量的估计残差；py为观测量的误差协方差阵；p
xy
为状态量和观测量的误差协方差阵。
[0059]
在无迹卡尔曼滤波器中，一步预测包含根据上一时步的后验估计通过确定性采样利用无迹(ut)变换来构造sigma点集。通过状态转移算子对sigma 点集进行非线性映射，得到sigma点集的一步预测期望(先验估计)。通过观测算子对一步预测后的sigma点集再次映射，得到观测量点集的预测期望同时更新误差协方差阵p
x,k|k-1
、p
y,k
、p
xy,k
。测量校正包含引入测量值校正sigma点集的先验估计得到修正后的状态量后验估计和误差协方差阵p
x,k
。测量校正以kk为权重，对一步预测值和测量值做加权融合。kk基于系统模型和测量值的不确定度得到。
[0060]
在集合卡尔曼滤波中，一步预测过程，也称为预报过程，包含根据上一时步的后验估计通过蒙特克罗随机抽样构造状态量背景集合。通过状态转移算子对背景集合进行非线性映射，得到背景集合的一步预测期望(先验估计)。通过观测算子对一步预测后的背景集合再次映射，得到观测量点集的预测期望k|k-1，同时更新误差协方差阵p
x,k|k-1
、p
y,k
、p
xy,k
。测量校正包含引入测量值校正背景集合的先验估计得到修正后的状态量后验估计和误差协方差阵p
x,k
。测量校正以kk为权重，对一步预测值和测量值做加权融合。kk基于系统模型和测量值的不确定度得到。
[0061]
为便于更清楚的理解本发明所提供的天然气井筒生产运行工况动态估计方法的具体应用方式，以下结合海上系统对上述各实施例做进一步说明，本领域相关技术人员当可知，该说明仅为便于理解上述各实施例，并不对其做任何限制。
[0062]
为了方便描述算法并展示中间结果，结合实际海上系统进行算法的建立。该生产系统包含两个生产井(系统拓扑如图6)，为了方便描述，两井标号分别为1#、2#。对于各支路，流量用gt表示，kg/s，油嘴开度用vo表示，井底观测点处的压力、温度分别用pt1、tm1表示，单位分别为mpa、℃，井口观测点处的压力、温度分别用pt2、tm2表示，油嘴后观测点处的压力、温度分别用pt3、tm3表示。平台观测点处的压力、温度分别用pt4、tm4表示。
[0063]
利用天然气生产系统的常规工艺参数和生产数据，以井筒为单位进行黑箱/灰箱模型的辨识和测试。生产数据包含dcs系统采样的观测数据和模拟数据，用于训练模型，并作为模型测试的参照。其中，各观测点的温度、压力数据来自现场采样，而单井流量来自软
件模拟。在软件模拟中，以实际天然气生产井的工艺参数和生产数据为基础进行建模，在井底、井口压力边界条件下计算得到单井流量。常规工艺参数，包括组分、井身结构等。训练数据和测试数据见图7a至图7d，涵盖各观测点的温度、压力采样数据和相应的流量模拟数据。
[0064]
一、基于生产数据，根据所述生产数据，利用传统机器学习、深度学习等辨识技术，对天然气井筒进行灰箱系统辨识，建立包含mlp-narx、dnn-narx和 cnn-narx等多种机器学习模型的非线性自回归的灰箱动态辨识模型库，用以构建井筒状态空间方程；
[0065]
(1)、mlp-narx、dnn-narx和cnn-narx模型描述了输入量与输出量之间的动态约束关系。本发明输入为温度、压力或油嘴开度等观测量，输出为单井流量和压力等。离散时间下非线性单输入单输出(siso)或多输入多输出(mimo)系统可用抽象形式，下式表述：
[0066][0067][0068]
式中，k是离散时步；zk是输出向量；uk是输入向量；ek是误差向量；n是采样点个数；nz、nu、ne分别是输出变量和输入变量的个数；[z
k-1
,
…
,z
k-dz
]为历史输出时间序列，也称为自回归 (ar)项，代表了当前输出和历史输出的自相关关系；[u
k-1
,
…
,u
k-du
]为历史输入时间序列，也称为外源性输入(xi)项，代表了当前输出和历史输入的互相关关系；dz、 du分别是ar项、xi项的最大时滞；ar项中，[1,
…
,dz]是时滞向量，不一定连续取值，[1,
…
,du]同理。ω、b是dnn的参数，分别代表节点权重和偏移量。
[0069]
综上，mlp-narx、dnn-narx和cnn-narx模型描述了ar项、xi项和当前时刻输出zk之间的动态关系。可包含所有输出量，这是因为神经网络类模型的数据挖掘能力和稳定性更高。
[0070]
(2)、建立黑箱动态模型。
[0071]
mlp-narx模型的辨识第一步：确定输入量和输出量、时滞向量、选取动态训练样本和稳态训练样本、使用传统梯度下降法和sigmoid激活函数、不使用各层规范化和“丢弃”(dropout)技术，对模型进行有监督预训练，以确定mlp-narx模型超参数(网络层数、各层节点数、学习率、动量因子)和参数(各节点权重ω和偏移量b)，得到mlp-narx黑箱动态模型，描述井筒非线性动态特性。
[0072]
dnn-narx模型的辨识第一步：确定输入量和输出量、时滞向量、选取动态训练样本和稳态训练样本、使用小批量随机梯度下降法(mini-batch sgd)、选用“泄漏修正线性单元”(leakyrelu)作为激活函数、辅以各层规范化和“丢弃”(dropout) 技术，对模型进行有监督预训练，确定网络模型超参数(网络层数、各层节点数、学习率、丢弃率)和参数(各节点的权重ω和偏移量b)，得到dnn-narx黑箱动态模型，描述井筒非线性动态特性。
[0073]
cnn-narx模型辨识第一步：确定输入量和输出量、确定二维输入矩阵的维度 (时滞向量和输入变量个数)、选取动态训练样本和稳态训练样本、选用“修正线性单元”(relu)
作为激活函数、辅以批量规范化和“丢弃”(dropout)技术，对模型进行有监督预训练，以确定cnn-narx模型超参数(卷积层和池化层的层数，卷积层的卷积核个数、卷积核扫描尺寸、移动步长，池化层的池化方法、池化层过滤器尺寸和步长，学习率、动量因子、丢弃率、批量个数等)和全连接层参数(各节点的权重ω和偏移量b)，得到cnn-narx黑箱动态模型，描述井筒非线性动态特性。
[0074]
确定输入量和输出量，对井筒分别建立3个黑箱模型：mlp-narx黑箱模型， dnn-narx黑箱模型和cnn-narx黑箱模型，分别命名为g1、g2、g3。g1、g2、 g3的输入量包含所有可观测的支路和平台变量(pt2除外)。
[0075]
对于1#，g1、g2、g3的输入量为u＝{vo,pt1,tm1,tm2,pt3,tm3,pt4,tm4}，输出量为z＝{gt,pt2}。
[0076]
对于2#，g1、g2、g3的输入量为u＝{vo,pt1,tm1,tm2,pt3,tm3,pt4,tm4}，输出量为z＝{gt,pt2}。
[0077]
选取动态训练样本{u,z}和稳态训练样本样本点个数分别用n、ns表示。取图7a至图7d中前50％的生产数据作为动态大样本，用于训练mlp-narx、 dnn-narx和cnn-narx黑箱模型。在一定分布下抽样得到稳态小样本用于微调所述三种黑箱模型的参数。模型测试样本为图7a至图7d中所有生产数据。通过有监督预训练和模型测试，确定所述三种黑箱模型g1、g2、g3的超参数和参数。 g1超参数包括网络层数、各层节点数、学习率、动量因子等，参数包括各节点的权重ω和偏移量b；g2超参数包括网络层数、各层节点数、学习率、动量因子、批量大小、丢弃率、学习次数等，参数包括各节点的权重ω和偏移量b；g3超参数包括卷积层和池化层的层数，卷积层的卷积核个数、卷积核扫描尺寸、移动步长，池化层的池化方法、池化层过滤器尺寸和步长，学习率、丢弃率、批量大小等，参数包括全连接层各节点的权重ω和偏移量b。
[0078]
确定时滞向量：对输出序列的随机成分进行偏自相关分析，对输出-输入序列的随机成分进行互相关分析，分别得到偏自相关函数图(pacf)和互相关函数图(ccf)，然后根据pacf和ccf，确定ar项和xi项的时滞向量。各井输出量的pacf在2 阶或5阶后不再显著，ar项的dz取2或5，ar项的时滞向量为1：dz。xi项的时滞向量取ccf峰值对应时滞的并集。dnn和mlp的输入是自回归项和外源输入项组成的向量；cnn每个批量的输入是二维矩阵，矩阵的行是时滞向量，矩阵的列是各个不考虑时滞的外源输入变量和自回归变量。
[0079]
训练mlp-narx黑箱模型时，优化方法使用传统梯度下降法；参数初始化方法使用xavier初始化法，令各层的激活值和状态梯度的方差在传播过程中的方差保持一致；激活函数使用sigmoid函数：不使用各层规范化和丢弃技术；学习率0.001，动量因子0.9，批量大小100，学习次数300；在模型训练和测试过程中，使用均方根误差(rmse)、平均绝对误差(mae)、均方误差(mse)和相关系数平方(r2) 作为模型训练效果的评估指标。训练dnn-narx黑箱模型时，优化算法使用小批量随机梯度下降法(mini-batch sgd)，防止产生局部极值；参数初始化方法使用xavier 初始化法；激活函数使用“泄漏修正线性单元”(leakyrelu)；各层规范化使用第一类规范化(l1：对点规范化，针对各神经元的输出进行规范化)或第二类规范化(l2：对边规范化，针对连接神经元的参数进行规范化)；学习率0.001；动量因子0.9；批量大小100；应用“丢弃”(dropout)技术，也就是输入层和隐藏层的神经元以一定的概率被随机舍弃，能够降低神经网络的复杂度，防止过拟合，丢弃率是0.4～0.6；学习次数
为300；在模型训练和测试过程中所用误差评估指标同mlp-narx模型。
[0080]
训练cnn-narx黑箱模型时，初始化方法使用xavier初始化方法；激活函数使用relu函数；池化方法使用最大池化法；学习率0.01，动量因子0.9，批量大小200；应用“丢弃”(dropout)技术，丢弃率是0.4～0.6；学习次数为100；在模型训练和测试过程中所用误差指标同mlp-narx模型。
[0081]
(3)、建立灰箱动态模型。应用迁移学习技术引入稳态小样本，分别对 mlp-narx、dnn-narx和cnn-narx黑箱模型的参数进行微调，获得适用范围更广的mlp-narx、dnn-narx和cnn-narx灰箱模型，实现模型训练过程中动态大样本和稳态小样本的结合。
[0082]
在迁移学习时(图8)，需要冻结浅层的结构和参数，利用少量有标注的新样本，对较深层的结构或参数进行训练、微调，保证模型适应新样本的同时仍适应旧样本。
[0083]
进行微调时，根据实际需要，可以仅改变微调层的参数，也可以同时改变微调层的节点数和参数，即保持mlp-narx、dnn-narx和cnn-narx黑箱模型g1、g2、 g3的结构不变，仅改变微调层的参数来获得相应的灰箱模型g4、g5、g6。
[0084]
二、在上述mlp-narx、dnn-narx和cnn-narx模型的辨识基础上，重构天然气井筒生产状态空间方程，以“衰减记忆的ukf”ukf和enkf，通过“一步预测”后的“测量校正”两步实现不断递推估计。
[0085]
(1)在上述mlp-narx、dnn-narx和cnn-narx模型的辨识基础上，将加性噪声ukf的非线性状态空间方程，通过构建的状态转移方程，观测矩阵构建的观测方程，重构天然气井筒生产状态空间方程。
[0086][0087]
式中，z包含两个输出量z＝{gt， pt2}，nz＝2，因此x＝{gt，pt2}；状态向量中包含观测向量，用于驱动滤波器，即 x＝{gt,pt2}，y＝{pt2}；输入向量u＝{vo,pt1,tm1,tm2,pt3,tm3,pt4,tm4}。ck由0、 1构成，式中，0是全零列，i是单位阵。
[0088]
(2)初始化状态量的分布并按以下两步不断进行递推估计，便可得到每个时步的状态量后验估计值(图9)。以上滤波步骤用下式概括：
[0089][0090][0091]
式中，“forcast”代表一步预测阶段，“dataassimilation”代表测量校正阶段。k是
离散时步；下标“k”代表k时步的后验估计，下标“k|k-1”代表k-1时步到k时步的一步预测值或先验估计；x是状态向量，y是观测向量；k是滤波增益；是观测量的测量值；p
x
是状态量的误差协方差阵。式中，分别为sigma点集、观测量点集的期望；是k时步的新息，即观测量的估计残差；py为观测量的误差协方差阵；p
xy
为状态量和观测量的误差协方差阵。
[0092]
(a)无迹卡尔曼滤波器的一步预测包含根据上一时步的后验估计通过确定性采样利用无迹(ut)变换来构造sigma点集。通过状态转移算子对sigma 点集进行非线性映射，得到sigma点集的一步预测期望(先验估计)。通过观测算子对一步预测后的sigma点集再次映射，得到观测量点集的预测期望同时更新误差协方差阵p
x,k|k-1
、p
y,k
、p
xy,k
。
[0093][0094][0095][0096][0097][0098][0099][0100][0101]
无迹卡尔曼滤波器的测量校正过程包含，引入测量值校正sigma点集的先验估计得到修正后的状态量后验估计和误差协方差阵p
x,k
。测量校正以kk为权重，对一步预测值和测量值做加权融合。kk基于系统模型和测量值的不确定度得到。
[0102][0103][0104][0105]
式中，分别为sigma点集、观测量点集的期望；是k时步的新息，即观测量的估计残差；w
i(m)
、w
i(c)
分别为计算点集的一步预测均值和误差协方差时各点的权
重；py为观测量的误差协方差阵；p
xy
为状态量和观测量的误差协方差阵； n、κ、α、β是ukf的可调参数。
[0106]
(b)在集合卡尔曼滤波中，一步预测过程，也称为预报过程，包含根据上一时步的后验估计通过蒙特克罗随机抽样构造状态量背景集合。通过状态转移算子对背景集合进行非线性映射，得到背景集合的一步预测期望(先验估计)。通过观测算子对一步预测后的背景集合再次映射，得到观测量点集的预测期望k|k-1，同时更新误差协方差阵p
x,k|k-1
、p
y,k
、p
xy,k
。背景集合中点的个数n
bg
设为100，enkf 的计算式如下：
[0107][0108][0109][0110][0111][0112][0113][0114]wi(m)
＝w
i(c)
＝1/n
bg i＝0l n
bg
[0115]
测量校正过程，也称作预报过程，包含引入测量值校正背景集合的先验估计得到修正后的状态量后验估计和误差协方差阵p
x,k
。测量校正以kk为权重，对一步预测值和测量值做加权融合。kk基于系统模型和测量值的不确定度得到。
[0116][0117][0118][0119]
式中，分别为背景集合、观测量点集的期望；是k时步的新息，即观测量的估计残差；w
i(m)
、w
i(c)
分别为计算点集的一步预测均值和误差协方差时各点的权重；py为观测量的误差协方差阵；p
xy
为状态量和观测量的误差协方差阵；
[0120]
由此，利用本发明所提供的天然气井筒生产运行工况动态估计方法可对单个或者多个天然气井筒生产运行工况的在线估计。以该方法为基础编写的软件可以实现在线实时监测生产气井的井口生产运行工况，其优势在于：数据驱动估计算法不仅能够进行动态估计，而且计算量小、计算速度快，工程适用性更强；克服了传统的计量设备技术既难于应用于水下和沙漠深处的生产系统，即便应用也难于维护的难题。仅需要在工程实际中需要不断利用新的生产数据更新黑箱/灰箱模型和相应的滤波器，因此该方法应用成本也相对减
少很多。
[0121]
请参考图10所示，本发明还提供一种天然气井筒生产运行工况动态估计装置，所述装置包含：辨识模块、重构模块和估计模块；所述辨识模块用于获取生产数据，根据所述生产数据利用学习算法构建的辨识技术对天然气井筒进行灰箱系统辨识，建立非线性自回归灰箱的动态辨识模型库，并构建井筒状态空间方程；所述重构模块用于根据所述动态辨识模型库中的辨识模型构建状态转移方程，并通过观测矩阵构建观测方程，通过所述状态转移方程和所述观测方程重构所述井筒状态空间方程；所述估计模块用于通过重构后的井筒状态空间方程结合无迹卡尔曼滤波-集合卡尔曼滤波的数据融合技术，融合所述辨识模型的估计值与井筒观测点的测量值；根据融合后的所述估计值和所述测量值融合所述无迹卡尔曼滤波器-集合卡尔曼滤波器的估计结果，以过程辨识和状态估计相结合的方式对天然气井筒生产运行工况进行动态估计，获得目标估计值。其中，所述辨识模块具体包含：根据所述生产数据，在未引入各层规范化和丢弃技术的基础上利用传统梯度下降法和sigmoid激活函数对神经网络模型进行有监督预训练，得到传统神经网络模型；以泄漏修正线性单元作为激活函数，根据所述生产数据通过随机梯度下降法配合各层规范化和丢弃技术对预设神经网络模型进行有监督预训练，得到深度前馈网络模型；以修正线性单元作为激活函数，加入卷积层、池化层、全连接层后，根据所述生产数据通过随机梯度下降法对神经网络模型进行有监督预训练，得到卷积神经网络模型；根据传统神经网络模型、深度前馈网络模型和卷积神经网络模型获得描述井筒非线性动态特性的所述动态辨识模型库。其中，获得描述井筒非线性动态特性的所述辨识模型还包含：根据稳态样本，通过迁移学习调整所述预训练所得辨识模型。
[0122]
在上述实施例中，各模块的具体应用原理及流程已在前述实施例中详细说明，在此就不再一一详述。
[0123]
本发明的有益技术效果在于：以天然气生产井筒为对象进行过程辨识，利用历史采样数据和软件模拟值，采用mlp-narx、dnn-narx和cnn-narx辨识模型，建立黑箱/灰箱模型库，描述井筒状态变量之间的动态约束关系，建立天然气井筒生产运行工况状态空间方程，构造ukf和enkf相结合的融合滤波器，以过程辨识和状态估计相结合的方式，用于井筒流量和压力的估计。利用该计算，可以替代实体流量计和传感器工作，并满足生产要求。该技术需要生产历史数据和在线数据，特别适用于采用dcs的气田生产系统。该技术具有较强的弹性和适应性，只要具有足够的基础参数和实时的工况数据，不仅可以同时预测多口天然气井筒生产运行工况，且不需要改造生产系统即可直接应用该技术。
[0124]
本发明还提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述方法。
[0125]
本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有执行上述方法的计算机程序。
[0126]
如图11所示，该电子设备600还可以包括：通信模块110、输入单元120、音频处理单元130、显示器160、电源170。值得注意的是，电子设备600也并不是必须要包括图11中所示的所有部件；此外，电子设备600还可以包括图11中没有示出的部件，可以参考现有技术。
[0127]
如图11所示，中央处理器100有时也称为控制器或操作控件，可以包括微处理器或其他处理器装置和/或逻辑装置，该中央处理器100接收输入并控制电子设备600 的各个部
件的操作。
[0128]
其中，存储器140，例如可以是缓存器、闪存、硬驱、可移动介质、易失性存储器、非易失性存储器或其它合适装置中的一种或更多种。可储存上述与失败有关的信息，此外还可存储执行有关信息的程序。并且中央处理器100可执行该存储器140 存储的该程序，以实现信息存储或处理等。
[0129]
输入单元120向中央处理器100提供输入。该输入单元120例如为按键或触摸输入装置。电源170用于向电子设备600提供电力。显示器160用于进行图像和文字等显示对象的显示。该显示器例如可为lcd显示器，但并不限于此。
[0130]
该存储器140可以是固态存储器，例如，只读存储器(rom)、随机存取存储器 (ram)、sim卡等。还可以是这样的存储器，其即使在断电时也保存信息，可被选择性地擦除且设有更多数据，该存储器的示例有时被称为eprom等。存储器140还可以是某种其它类型的装置。存储器140包括缓冲存储器141(有时被称为缓冲器)。存储器140可以包括应用/功能存储部142，该应用/功能存储部142用于存储应用程序和功能程序或用于通过中央处理器100执行电子设备600的操作的流程。
[0131]
存储器140还可以包括数据存储部143，该数据存储部143用于存储数据，例如联系人、数字数据、图片、声音和/或任何其他由电子设备使用的数据。存储器140 的驱动程序存储部144可以包括电子设备的用于通信功能和/或用于执行电子设备的其他功能(如消息传送应用、通讯录应用等)的各种驱动程序。
[0132]
通信模块110即为经由天线111发送和接收信号的发送机/接收机110。通信模块 (发送机/接收机)110耦合到中央处理器100，以提供输入信号和接收输出信号，这可以和常规移动通信终端的情况相同。
[0133]
基于不同的通信技术，在同一电子设备中，可以设置有多个通信模块110，如蜂窝网络模块、蓝牙模块和/或无线局域网模块等。通信模块(发送机/接收机)110还经由音频处理器130耦合到扬声器131和麦克风132，以经由扬声器131提供音频输出，并接收来自麦克风132的音频输入，从而实现通常的电信功能。音频处理器130 可以包括任何合适的缓冲器、解码器、放大器等。另外，音频处理器130还耦合到中央处理器100，从而使得可以通过麦克风132能够在本机上录音，且使得可以通过扬声器131来播放本机上存储的声音。
[0134]
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等) 上实施的计算机程序产品的形式。
[0135]
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0136]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特
定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0137]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0138]
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。