CN114693005A - 基于卷积傅里叶神经网络的三维地下油藏动态预测方法 - Google Patents

基于卷积傅里叶神经网络的三维地下油藏动态预测方法 Download PDF

Info

Publication number
CN114693005A
CN114693005A CN202210603016.6A CN202210603016A CN114693005A CN 114693005 A CN114693005 A CN 114693005A CN 202210603016 A CN202210603016 A CN 202210603016A CN 114693005 A CN114693005 A CN 114693005A
Authority
CN
China
Prior art keywords
pressure
convolution
fourier
saturation
dimensional
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202210603016.6A
Other languages
English (en)
Other versions
CN114693005B (zh
Inventor
张凯
左袁德
王晓雅
张黎明
刘丕养
严侠
张华清
杨勇飞
孙海
张文娟
姚军
樊灵
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN202210603016.6A priority Critical patent/CN114693005B/zh
Publication of CN114693005A publication Critical patent/CN114693005A/zh
Application granted granted Critical
Publication of CN114693005B publication Critical patent/CN114693005B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/044—Recurrent networks, e.g. Hopfield networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/048—Activation functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G06N3/084—Backpropagation, e.g. using gradient descent
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/02—Agriculture; Fishing; Forestry; Mining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Business, Economics & Management (AREA)
  • General Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computing Systems (AREA)
  • Artificial Intelligence (AREA)
  • Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Molecular Biology (AREA)
  • Computational Linguistics (AREA)
  • Biomedical Technology (AREA)
  • Strategic Management (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Human Resources & Organizations (AREA)
  • Economics (AREA)
  • Geometry (AREA)
  • Marketing (AREA)
  • Algebra (AREA)
  • General Business, Economics & Management (AREA)
  • Tourism & Hospitality (AREA)
  • Computational Mathematics (AREA)
  • Discrete Mathematics (AREA)
  • Animal Husbandry (AREA)
  • Databases & Information Systems (AREA)
  • Mining & Mineral Resources (AREA)
  • Marine Sciences & Fisheries (AREA)
  • Primary Health Care (AREA)
  • Medical Informatics (AREA)
  • Computer Vision & Pattern Recognition (AREA)

Abstract

本发明公开了一种基于卷积傅里叶神经网络的三维地下油藏动态预测方法，属于油藏工程技术领域，利用构建的3D卷积傅里叶神经网络模型对三维地下油藏剩余油饱和度或压力进行动态预测，具体包括以下步骤：采集三维地下油藏数据；构建3D卷积傅里叶神经网络，结合LSTM同时考虑时间和空间信息；设置卷积傅里叶网络模型的超参数，训练3D卷积傅里叶网络模型；评估3D卷积傅里叶神经网络的性能；训练完成后，输出测试指标良好的3D卷积傅里叶网络模型。本发明实现了三维油藏模型的高精度快速预测，很好的针对三维油藏的时空性和物理***的特点，更符合现场实际需求。

Description

基于卷积傅里叶神经网络的三维地下油藏动态预测方法
技术领域
本发明属于油藏工程技术领域，具体涉及一种基于3D卷积傅里叶神经网络模型的三维地下油藏剩余油饱和度或压力的生产动态预测方法。
背景技术
对多孔介质中的多相流进行数值模拟对许多地质科学应用来说是至关重要。在油藏工程领域，油水两相偏微分方程可以准确描述地下流动***，模拟地下流动需要进行大量的偏微分方程求解计算，这也是数值模拟技术的基础。传统的数值模拟基于数值计算方法如有限元法，有限差分法等对偏微分方程进行求解。油藏数值模拟可以准确的计算出未来油藏的压力、饱和度、产量等等油藏状态用于生产开发动态评价。而真实的三维油藏模型网格数多，模型大且边界复杂，求解往往需要耗费大量时间和计算资源。
随着信息化的深入和硬件水平的不断发展，深度学习等人工智能方法在图像识别、目标检测等场景得到了广泛的应用。而在油藏工程领域，数据驱动的机器学习方法通过构建神经网络代理模型，为传统的模拟器提供了更快的替代方案。机器学习代理模型与数值模拟数据的映射，为传统的模拟器提供了更快速的选择，可以有效的节约预测时间。目前流行的代理模型方法大多针对于二维油藏，对于一个三维的复杂的地质模型，网格数更多，同时要考虑到重力影响，以及油井穿透多个区域。所以三维油藏预测比二维油藏更具有挑战性，更具有实用价值。现有的一些应用于三维模型的方法只考虑图像之间的局部信息，不能考虑到物理信息和时空特征，存在在精度低，泛化性能差，训练时间慢等条件限制。
发明内容
针对当前使用油藏数值模拟进行剩余油饱和度或压力预测时间过长，普通机器学习代理模型不能有效的预测三维复杂油藏等问题，本发明提出了一种基于卷积傅里叶神经网络模型三维油藏剩余油饱和度或压力的生产动态预测方法，同时考虑到三维油藏的时空性和物理信息和重力因素的影响，有效提高预测精度，适应不同的油藏条件，节省计算时间。
本发明的技术方案如下：
一种基于卷积傅里叶神经网络的三维地下油藏动态预测方法，利用构建的3D卷积傅里叶神经网络模型对三维地下油藏剩余油饱和度或压力进行动态预测，具体包括以下步骤：
步骤1、采集三维地下油藏数据，利用SGeMS地质统计工具生成符合地质特征的三维渗透率场，利用数值模拟器求解出不同渗透率对应的一段时间内的饱和度或压力作为样本库，按比例划分数据集，渗透率作为输入，饱和度或压力作为输出；
步骤2、构建3D卷积傅里叶神经网络，将3D卷积和傅里叶变换相结合，利用3D卷积算子提取图像局部信息，利用傅里叶变换提取物理信息近似微分算子，在网络中添加重力约束；利用编码解码网络实现渗透率到油藏状态的映射，结合LSTM同时考虑时间和空间信息；
步骤3、设置卷积傅里叶网络模型的超参数，在训练集下训练3D卷积傅里叶网络模型；
步骤4、使用测试集评估3D卷积傅里叶神经网络的性能，计算RMSE和量化饱和度或压力的相对误差；
步骤5、训练完成后，输出测试指标良好的3D卷积傅里叶网络模型；基于油藏安装的监测设备实时采集三维地下油藏数据，利用地质统计工具在线生成三维油藏的渗透率场，然后利用3D卷积傅里叶网络模型在线预测该油藏未来一段时间内油藏饱和度或压力场分布，通过压力和饱和度可以进一步计算产量，为制定生产策略提供直观的参考。
进一步地，步骤1的具体过程如下：
步骤1.1、采集三维地下油藏数据，使用SGeMS地质统计模型生成三维渗透率场，网格尺寸为40×40×20，生成2000个渗透率样本；
步骤1.2、利用数值模拟器使用有限元方法计算三维油藏***的油水两相流动方程得到所有网格块中的饱和度或压力结果；油藏地下油水流动***是三维的，需要考虑重力效应，由质量守恒方程（1）和相达西速度方程（2）构成；
质量守恒方程表达式如（1），
Figure 393993DEST_PATH_IMAGE001
（1）
其中，
Figure 933558DEST_PATH_IMAGE002
表示梯度算子；j表示不同的相，包括油水两相，o表示油，w表示水；
Figure 331173DEST_PATH_IMAGE003
表 示相的密度；v j 表示相达西速度；
Figure 639794DEST_PATH_IMAGE004
表示源/汇项，上标l表示井处的参数；
Figure 713930DEST_PATH_IMAGE005
表示孔隙度；S j 表示相的饱和度；t表示时间；
相达西速度方程的表达式如（2），
Figure 317955DEST_PATH_IMAGE006
（2）
其中，k rj 表示相的相对渗透率，
Figure 378315DEST_PATH_IMAGE007
表示绝对渗透率向量，
Figure 666077DEST_PATH_IMAGE008
表示相的粘度，p j 表示 相的压力，g表示重力，z表示深度；
步骤1.3、按照8：2的比例分割数据集和测试集；输入数据渗透率为(n,nx,ny,nz,1)的张量，输出数据饱和度或压力为(n,nx,ny,nz,1,T)的张量，n是样本数，nx为横向网格数，ny为纵向网格数，nz为深度方向的网格数，T为时间步数量。
进一步地，步骤2的具体过程如下：
3D卷积傅里叶神经网络，主要包含编码网络和解码网络两个部分，考虑到时间序列影响增加了处理时间的循环神经模块；输入参数是已知观测值渗透率场，输出参数是未来一段时间的饱和度或压力；
编码网络结构包含M个3D卷积傅里叶层，每一个卷积傅里叶层结构如下，
Figure 52059DEST_PATH_IMAGE009
（3）
Figure 441583DEST_PATH_IMAGE010
（4）
其中，x是最初输入参数，W是权重，
Figure 305634DEST_PATH_IMAGE011
是激活函数，使用Relu激活；v N+1(x)是卷积 傅里叶模块N+1层的输出结果，v N (x)是N层的输出结果；
Figure 447902DEST_PATH_IMAGE012
积是三维傅里叶积分算 子，用以提取物理信息；
Figure 270365DEST_PATH_IMAGE013
是快速傅里叶变换，
Figure 6239DEST_PATH_IMAGE014
是傅里叶逆变换，C是截断超过最大阈 值的傅里叶级数的操作，傅里叶变换操作能通过傅里叶级数即一系列正余弦曲线去近似数 据所符合的物理规律；Conv3D是3D卷积结果，用以提取局部信息，每一层的卷积核大小和步 幅会进行调整，卷积核的大小均为3×3×3，卷积步幅有1×1×1和2×2×2两种；因为傅里 叶变换是无限维度映射的，使用傅里叶提取的结果需要与卷积结果维度一致，相加之后不 会改变原来维度，再经过五层傅里叶层，得到的最终结果进入解码网络；
解码网络结构由反卷积层构成，解码网络的部分层的输入由编码网络对应的输出矩阵和当前输入复制拼接而成，目的是保存更多信息防止卷积过程中出现的特征丢失；解码网络是将中间结果解码到原数据空间的过程，得到的结果为某一时间步的压力值或饱和度；
考虑到时间序列的影响，将多个时间步的输入经过解码网络的中间结果通过LSTM结构，结合时间序列信息，经过LSTM处理后对每个特征图用对应的解码网络得到每一个时间步的输出结果；
LSTM在t时刻的输入有三个：细胞状态C t-1，隐层状态h t-1，当前输入特征x t ；
C t 承载着空间时间信息，计算公式如下，
Figure 50812DEST_PATH_IMAGE015
(5)
其中，f t 代表遗忘门，控制丢弃细胞状态C t-1中的某些信息；
Figure 922953DEST_PATH_IMAGE016
代表新的候选细胞 状态，i t 代表更新门，选择保留候选细胞状态
Figure 775372DEST_PATH_IMAGE016
中的某些信息；
隐层输出结果h t 由输出门o t 和上一步计算出的C t 经过tanh激活函数计算得出，如下，
Figure 608330DEST_PATH_IMAGE017
(6)
遗忘门f t ，更新门i t ，输出门o t ，以及新的候选细胞状态
Figure 814183DEST_PATH_IMAGE016
计算公式分别如下，每 种门有各自对应的权重和偏置，LSTM中的
Figure 931044DEST_PATH_IMAGE018
是sigmoid激活函数，激活结果在（0，1）之间， tanh激活结果在（-1，1）之间；
Figure 829729DEST_PATH_IMAGE019
(7)
Figure 382939DEST_PATH_IMAGE020
(8)
Figure 392484DEST_PATH_IMAGE021
(9)
Figure 895009DEST_PATH_IMAGE022
(10)
其中，W f 、W i 、W o 和W c 分别为遗忘门、更新门、输出门和新的候选细胞状态对应的权重；b f 、b i 、b o 和b c 分别为遗忘门、更新门、输出门和新的候选细胞状态对应的偏置；
输入渗透率场经过卷积傅里叶编码层得到隐层变量，再经过LSTM结构和解码网络 得到同时包含时间信息和空间信息的饱和度或压力的预测结果
Figure 699017DEST_PATH_IMAGE023
。
进一步地，步骤3的具体过程如下：
步骤3.1、设置3D卷积傅里叶网络模型的输入输出；输入参数是三维油藏的渗透率；设置输出值，预测未来一段时间内的油藏饱和度或压力分布；
步骤3.2、初始化超参数，批处理大小为8，激活函数使用ReLu函数，随机初始化权重W；超参数根据网络评估性能进行优化调整，损失函数使用MSE函数计算损失，训练压力时MSE损失函数值为：编码解码网络使用相同的超参数设置，训练在200个epoch内收敛；
饱和度的损失函数loss sat_MSE 如下：
Figure 772147DEST_PATH_IMAGE024
（11）
其中，y sat_i 与
Figure 319803DEST_PATH_IMAGE025
分别为训练集中第i个样本的输出饱和度真实值和3D卷积傅 里叶神经网络模型的饱和度计算结果，n train 为训练集样本总数；
对于钻过多个区块的直井，在深度d层处的井筒压力
Figure 411255DEST_PATH_IMAGE026
受深度z变化的结引力的 影响；最上面的射孔处的井筒压力称为井底压力，继续随深度往下的d+1层处深度的井筒压 力
Figure 386165DEST_PATH_IMAGE027
计算如方程（12），
Figure 336803DEST_PATH_IMAGE028
是深度差，
Figure 64981DEST_PATH_IMAGE029
是流体混合密度；
当3D卷积傅里叶网络模型对压力进行预测时，考虑在深度方向上压力受到重力作用的影响，压力满足方程（12），
Figure 886306DEST_PATH_IMAGE030
（12）
压力的损失函数loss press_MSE 如下：
Figure 422330DEST_PATH_IMAGE031
（13）
其中，y press_i 与
Figure 594685DEST_PATH_IMAGE032
分别为训练集中第i个样本的输出压力真实值和3D卷积傅 里叶神经网络模型的压力计算结果；
Figure 359510DEST_PATH_IMAGE033
是井所在网格处的d+1层深度的井筒压力实际 值，
Figure 300921DEST_PATH_IMAGE034
是井所在网格处的d层深度的井筒压力预测值；
Figure 742267DEST_PATH_IMAGE035
是数据约束的系数，
Figure 136339DEST_PATH_IMAGE036
是方程 约束的系数；
网络根据设定好的超参数和网络迭代结构进行正向传播；
步骤3.3、根据步骤3.2计算得到的损失误差进行反向传播，利用ADAM优化器，更新权重W，使MSE损失逐渐减小；ADAM优化算法的初始学习率lr=0.003，批量大小为8，混合权重为1000；重复这个过程，在200个epoch之后误差下降到接受范围之内，训练完成。
进一步地，步骤4的具体过程如下：
对于所有的测试样本，计算相对均方根误差RMSE，量化预测饱和度或压力的相对误差；
RMSE计算公式为：
Figure 203390DEST_PATH_IMAGE037
（14）
其中，n test 为测试集样本总数；
Figure 999308DEST_PATH_IMAGE038
为模型预测结果，包括
Figure 877134DEST_PATH_IMAGE039
与
Figure 758502DEST_PATH_IMAGE040
；y i 为数 值模拟解真实值，包括y sat_i 与y press_i ；
计算三维油藏的饱和度或压力相对误差的公式如下：
Figure 865130DEST_PATH_IMAGE041
（15）
饱和度的相对误差为
Figure 515554DEST_PATH_IMAGE042
，压力的相对误差为
Figure 298702DEST_PATH_IMAGE043
，n test 是测试样本总数，n M 是网格 总数，n t 是时间步数，
Figure 667367DEST_PATH_IMAGE044
是测试集中第i个样本第m个网格t时刻的饱和度的预测结果，
Figure 67431DEST_PATH_IMAGE045
是对应饱和度的真实值；
Figure 572362DEST_PATH_IMAGE046
是测试集中第i个样本第m个网格t时刻压力的预测结 果，
Figure 791990DEST_PATH_IMAGE047
是对应压力的真实值；利用最大压力减去最小压力的方式进行压力归一化，
Figure 647951DEST_PATH_IMAGE048
是第i个测试样本t时刻的最大压力，
Figure 830802DEST_PATH_IMAGE049
是第i个测试样本t时刻的最小压力。
进一步地，步骤5的具体过程如下：
保存性能指标评价较好的3D卷积傅里叶神经网络模型，将训练完成的权重和超参数保存下来，利用3D卷积傅里叶神经网络模型快速预测出三维油藏的未来的饱和度或压力状态；输入是三维油藏的渗透率场，输出是一段时间t内的压力和饱和度分布，输入经过一次前向传播的计算即可得到输出；
利用压力和饱和度通过peaceman模型计算出油水井产量，计算公式如下，
Figure 580452DEST_PATH_IMAGE050
（16）
其中，
Figure 846348DEST_PATH_IMAGE051
是第m个网格中井的油/水相的注入/采油速率，k m 为第m个网格的绝 对渗透率，k r, j 为k m 对应的油/水相对渗透率，S w, m 为第m个网格处的含水饱和度值，
Figure 297927DEST_PATH_IMAGE052
是 网格厚度，
Figure 940261DEST_PATH_IMAGE053
，默认
Figure 544418DEST_PATH_IMAGE054
，油藏在x与y两个方向上宽度相等，r w 是井 眼半径，
Figure 981215DEST_PATH_IMAGE055
是水或油的粘度，P m 与P l 分别是第m个网格处压力与井眼压力。
本发明所带来的有益技术效果：
本发明提供基于3D卷积傅里叶模型剩余油饱和度或压力的生产动态预测方法，利用卷积傅里叶层同时提取局部信息和物理信息，利用编码解码网络实现渗透率到饱和度或压力的转化，同时在网络中采用了LSTM结构考虑时序问题，在神经网络中同时考虑到空间信息和时间信息。在预测压力时利用物理方程在损失函数添加约束，在神经网络中考虑重力因素对压力变化的影响。本发明实现了三维油藏模型的高精度快速预测，很好的针对三维油藏的时空性和物理***的特点，更符合现场实际需求。
附图说明
图1为本发明基于3D卷积傅里叶神经网络三维油藏剩余油饱和度或压力的生产动态预测方法的流程示意图；
图2为本发明3D卷积傅里叶神经网络模型的结构示意图；
图3为本发明卷积傅里叶层的结构示意图；
图4为本发明编码解码结构和LSTM结合的示意图；
图5为本发明LSTM结构详细示意图；
图6为本发明实验中的相对渗透率曲线；
图7为本发明实验中均方误差RMSE效果图；
图8为本发明实验中生产井1和生产井2的油井产量真实值和测试值对比；
图9为本发明实验中生产井3和生产井4的油井产量真实值和测试值对比。
具体实施方式
下面结合附图以及具体实施方式对本发明作进一步详细说明：
本发明通过3D卷积傅里叶神经网络构建了神经网络代理模型，利用卷积傅里叶层同时提取局部信息和物理信息，利用编码解码网络实现渗透率到饱和度或压力的转化。该方法针对于符合现场实际的三维油藏，能够很好的替代数值模拟方法，快速的进行油藏生产动态预测。与现有的代理模型相比，本方法能更好的处理三维油藏模型，预测精度和泛化性更好。
参照图1，本发明提出了一种基于卷积傅里叶神经网络的三维地下油藏动态预测方法，利用构建的3D卷积傅里叶神经网络模型对三维地下油藏剩余油饱和度或压力进行动态预测，具体包括以下步骤：
步骤1、采集三维地下油藏数据，利用SGeMS地质统计工具生成符合地质特征的三维渗透率场，利用数值模拟器求解出不同渗透率对应的一段时间内的油藏状态即饱和度或压力作为样本库，按比例划分训练测试集数据集，渗透率作为输入，饱和度或压力作为输出。
具体过程如下：
步骤1.1、采集三维地下油藏数据，使用SGeMS地质统计模型生成三维渗透率场，网格尺寸为40×40×20，生成2000个渗透率样本。
步骤1.2、利用数值模拟器使用有限元方法计算三维油藏***的油水两相流动方程得到所有网格块中的饱和度或压力结果。油藏地下油水流动***是三维的，需要考虑到重力效应，由质量守恒方程（1）和相达西速度方程（2）构成。
质量守恒方程表达式如（1），
Figure 421555DEST_PATH_IMAGE056
（1）
其中，
Figure 336421DEST_PATH_IMAGE057
表示梯度算子；j表示不同的相，包括油水两相，o表示油，w表示水；
Figure 60664DEST_PATH_IMAGE058
表 示相的密度；v j 表示相达西速度；
Figure 668362DEST_PATH_IMAGE059
表示源/汇项（上标l表示井处的参数）；
Figure 97463DEST_PATH_IMAGE060
表示孔隙度；S j 表示相的饱和度；t表示时间；
相达西速度方程的表达式如（2），
Figure 816020DEST_PATH_IMAGE061
（2）
其中，k rj 表示相的相对渗透率，
Figure 394769DEST_PATH_IMAGE062
表示绝对渗透率向量，
Figure 173370DEST_PATH_IMAGE063
表示相的粘度，p j 表示 相的压力，g表示重力，z表示深度。
步骤1.3、按照8：2的比例分割数据集和测试集。输入数据渗透率为(n,nx,ny,nz,1)的张量，输出数据饱和度或压力为(n,nx,ny,nz,1,T)的张量，n是样本数，nx为横向网格数，ny为纵向网格数，nz为深度方向的网格数，T为时间步数量。
步骤2、构建如图2-图5所示的3D卷积傅里叶神经网络，将3D卷积和傅里叶变换相结合，利用3D卷积算子提取图像局部信息，利用傅里叶变换提取物理信息近似微分算子，在网络中添加重力约束。利用编码解码网络实现渗透率到油藏状态的映射，结合LSTM同时考虑时间和空间信息，最终输出未来某一段时间的饱和度或压力。具体过程如下：
3D卷积傅里叶神经网络，主要包含编码网络和解码网络两个部分，考虑到时间序列影响增加了处理时间的循环神经模块（RNN）。输入参数是已知观测值渗透率场。输出参数是未来一段时间的饱和度或压力。
编码网络结构包含M个3D卷积傅里叶层，如图4所示，每一个卷积傅里叶层结构如下。
Figure 588301DEST_PATH_IMAGE064
（3）
Figure 110550DEST_PATH_IMAGE065
（4）
其中，x是最初输入参数，W是权重，
Figure 278226DEST_PATH_IMAGE066
是激活函数，使用Relu激活；v N+1(x)是卷积 傅里叶模块N+1层的输出结果，v N (x)是N层的输出结果。
Figure 227727DEST_PATH_IMAGE067
积是三维傅里叶积分算 子，用以提取物理信息。
Figure 362911DEST_PATH_IMAGE068
是快速傅里叶变换，
Figure 688850DEST_PATH_IMAGE069
是傅里叶逆变换，C是截断超过最大阈 值的傅里叶级数的操作，傅里叶变换操作能通过傅里叶级数即一系列正余弦曲线去近似数 据所符合的物理规律。Conv3D是3D卷积结果，用以提取局部信息，每一层的卷积核大小和步 幅可能会进行调整，卷积核的大小均为3×3×3，卷积步幅有1×1×1和2×2×2两种。因为 傅里叶变换是无限维度映射的，在这里使用傅里叶提取的结果需要与卷积结果维度一致， 相加之后不会改变原来维度，再经过五层傅里叶层，即该层不包含公式（3）中Conv3D这一项 操作的网络层，得到的最终结果会进入解码网络。
解码网络通过反卷积变换来实现，同时结合LSTM在模拟时间序列变换上的优势。解码网络结构是由反卷积层构成，解码网络的部分层的输入是由编码网络对应的输出矩阵和当前输入复制拼接的而成，目的是保存更多信息防止卷积过程中出现的特征丢失。解码网络是将中间结果解码到原数据空间的过程，得到的结果即为某一时间步的压力值或饱和度。考虑到时间序列的影响，将多个时间步的输入经过解码网络的中间结果通过LSTM结构，结合时间序列信息，LSTM结构如图2中所示，经过LSTM处理后对每个特征图用对应的解码网络得到每一个时间步的输出结果。
LSTM的迭代结构如图5所示，LSTM在t时刻的输入有三个：细胞状态C t-1，隐层状态h t-1，当前输入特征x t 。
C t 承载着空间时间信息，计算公式如下，
Figure 976612DEST_PATH_IMAGE070
(5)
其中，f t 代表遗忘门，控制丢弃细胞状态C t-1中的某些信息；
Figure 97015DEST_PATH_IMAGE071
代表新的候选细胞 状态，i t 代表更新门，选择保留候选细胞状态
Figure 486539DEST_PATH_IMAGE072
中的某些信息；
隐层输出结果h t 由输出门o t 和上一步计算出的C t 经过tanh激活函数计算得出，如下，
Figure 350590DEST_PATH_IMAGE073
(6)
遗忘门f t ，更新门i t ，输出门o t ，以及新的候选细胞状态
Figure 492858DEST_PATH_IMAGE074
计算公式分别如下，每 种门有各自对应的权重和偏置，与网络中其他情况不同的是，LSTM中的
Figure 49742DEST_PATH_IMAGE075
是sigmoid激活 函数，激活结果在（0，1）之间，tanh激活结果在（-1，1）之间。
Figure 896868DEST_PATH_IMAGE076
(7)
Figure 423665DEST_PATH_IMAGE077
(8)
Figure 561385DEST_PATH_IMAGE078
(9)
Figure 554749DEST_PATH_IMAGE079
(10)
其中，W f 、W i 、W o 和W c 分别为遗忘门、更新门、输出门和新的候选细胞状态对应的权重；b f 、b i 、b o 和b c 分别为遗忘门、更新门、输出门和新的候选细胞状态对应的偏置；
输入渗透率场经过卷积傅里叶编码层得到隐层变量，再经过LSTM结构和解码网络 得到同时包含时间信息和空间信息的饱和度或压力的预测结果
Figure 122127DEST_PATH_IMAGE080
。
图3表示了单个时间步的具体的编码和解码结构细节，首先经过6层傅里叶卷积层，提取局部和物理信息，中间变量维度大小发生了改变，再经过五层傅里叶层，不改变维度。解码结构中，会将对应的大小相同的编码层的输出与上一解码层输出复制拼接到一起。最终中间变量会被解码得到每一个时间步的预测结果。具体形状大小数值变化情况如表1所示。
步骤3、设置卷积傅里叶网络模型的超参数，在训练集下训练3D卷积傅里叶网络模型。具体过程如下：
步骤3.1、设置3D卷积傅里叶网络模型的输入输出；输入参数是三维油藏的渗透率；设置输出值，预测未来一段时间内的油藏饱和度或压力分布；
步骤3.2、初始化超参数，批处理大小为8，激活函数使用ReLu函数，随机初始化权重W；超参数根据网络评估性能进行优化调整，损失函数使用MSE函数计算损失，训练压力时MSE损失函数值为：编码解码网络使用相同的超参数设置，训练在200个epoch内收敛；
饱和度的损失函数loss sat_MSE 如下：
Figure 718194DEST_PATH_IMAGE081
（11）
其中，y sat_i 与
Figure 710421DEST_PATH_IMAGE082
分别为训练集中第i个样本的输出饱和度真实值和3D卷积傅 里叶神经网络模型的饱和度计算结果，n train 为训练集样本总数。
对于钻过多个区块的直井，在深度d层处的井筒压力
Figure 717429DEST_PATH_IMAGE083
受深度z变化的结引力的 影响。最上面的射孔处的井筒压力称为井底压力（BHP），继续随深度往下的d+1层处深度的 井筒压力
Figure 427896DEST_PATH_IMAGE084
计算如方程（12），
Figure 562074DEST_PATH_IMAGE085
是深度差，
Figure 408807DEST_PATH_IMAGE086
是流体混合密度。
当3D卷积傅里叶网络模型对压力进行预测时，为了提高预测精度，需要考虑在深度方向上压力受到重力作用的影响，压力满足方程（12），
Figure 353760DEST_PATH_IMAGE087
（12）
压力的损失函数loss press_MSE 如下：
Figure 285944DEST_PATH_IMAGE088
（13）
其中，y press_i 与
Figure 223813DEST_PATH_IMAGE089
分别为训练集中第i个样本的输出压力真实值和3D卷积傅 里叶神经网络模型的压力计算结果；
Figure 190632DEST_PATH_IMAGE090
是井所在网格处的d+1层深度的井筒压力实际值，
Figure 745635DEST_PATH_IMAGE091
是井所在网格处的d层深度的井筒压力预测值；
本发明在压力的损失函数中将重力约束作为惩罚相，嵌入物理意义。使方程在训 练过程中会同时使输出满足数据约束和方程约束，
Figure 555328DEST_PATH_IMAGE092
是数据约束的系数，
Figure 906675DEST_PATH_IMAGE093
是方程约束 的系数，都是网络的超参数，分别初始化为0.9和0.1，在训练过程中会自动更新。
网络根据设定好的超参数和网络迭代结构进行正向传播；
步骤3.3、根据步骤3.2计算得到的损失误差进行反向传播，利用ADAM优化器，更新权重W，使MSE损失逐渐减小。ADAM优化算法的初始学习率lr=0.003，批量大小为8，混合权重为1000；重复这个过程，一般在200个epoch之后误差会下降到接受范围之内，即可训练完成。
步骤4、使用测试集来评估3D卷积傅里叶神经网络的性能，计算RMSE和量化饱和度或压力的相对误差。具体过程如下：
对于所有的测试样本，计算相对均方根误差RMSE，并且量化了预测饱和度或压力的相对误差；RMSE计算公式为：
Figure 603367DEST_PATH_IMAGE094
（14）
其中，n test 为测试集样本总数；
Figure 749177DEST_PATH_IMAGE095
为评估预测解，包括
Figure 780587DEST_PATH_IMAGE096
与
Figure 935625DEST_PATH_IMAGE097
；y i 为数值 模拟解真实值，包括y sat_i 与y press_i 。RMSE值越接近0，两参数之间的相似度越高。
计算三维油藏的饱和度或压力相对误差的公式如下：
Figure 985358DEST_PATH_IMAGE041
（15）
饱和度的相对误差为
Figure 567649DEST_PATH_IMAGE098
，压力的相对误差为
Figure 820776DEST_PATH_IMAGE099
，n test 是测试样本总数400，n M 是网 格总数40×20×20=16000，n t 是时间步数为20，
Figure 779505DEST_PATH_IMAGE100
是测试集中第i个样本第m个网格t时 刻的饱和度的预测结果，
Figure 450789DEST_PATH_IMAGE101
是对应饱和度的真实值；
Figure 203981DEST_PATH_IMAGE102
是测试集中第i个样本第m个 网格t时刻压力的预测结果，
Figure 944404DEST_PATH_IMAGE103
是对应压力的真实值。饱和度的范围在0-1之间不需要归 一化，而压力则不在这一范围内需要进行归一化，利用最大压力减去最小压力的方式进行 压力归一化，
Figure 441244DEST_PATH_IMAGE104
是第i个测试样本t时刻的最大压力，
Figure 214639DEST_PATH_IMAGE105
是第i个测试样本t时刻的最 小压力。相对误差越接近于0，说明对三维油藏的饱和度或压力的预测越精准。
步骤5、训练完成后，输出测试指标良好的3D卷积傅里叶网络模型。基于油藏安装的监测设备实时采集三维地下油藏数据，利用地质统计工具在线生成三维油藏的渗透率场，然后利用3D卷积傅里叶网络模型在线预测该油藏未来一段时间内油藏饱和度或压力场分布，通过压力和饱和度可以进一步计算产量，为制定生产策略提供直观的参考。具体过程如下：
保存性能指标评价较好的3D卷积傅里叶神经网络模型，即将训练完成的权重和超参数保存下来，利用3D卷积傅里叶神经网络模型可以快速预测出三维油藏的未来的饱和度或压力状态。输入是三维油藏的渗透率场，输出是一段时间t内的压力和饱和度分布，输入经过一次前向传播的计算即可得到输出。
利用压力和饱和度通过peaceman模型计算出油水井产量，计算公式如下，
Figure 263367DEST_PATH_IMAGE106
（16）
其中，
Figure 632031DEST_PATH_IMAGE107
是第m个网格中井的油/水相的注入/采油速率，也是要计算的注采量 结果，k m 为第m个网格的绝对渗透率，k r, j 为k m 对应的油/水相对渗透率，S w, m 为第m个网格处 的含水饱和度值，
Figure 542349DEST_PATH_IMAGE108
是网格厚度，
Figure 47280DEST_PATH_IMAGE109
，默认
Figure 266909DEST_PATH_IMAGE110
，油藏在x与y两个 方向上宽度相等，r w 是井眼半径，
Figure 857290DEST_PATH_IMAGE111
是水或油的粘度，P m 与P l 分别是第m个网格处压力与井 眼压力。根据以上预测结果，可以进一步分析油藏未来的状态，为制定生产策略提供直观的 参考。
为了证明本发明的可行性，进行了如下实验进行说明。
实验数据来源于某地区的一个三维油田区块，该模型包含40×40×20网格块，每 个网格块在下x,y,z方向的尺寸均为20m。四口产油井，两口注水井，一口注水井和两口产油 井射开油藏上面的10层，而另一个注水井和两个产油井射开油藏底部10层。所有井在定压 条件生产，产油井的井底压力设为33.5Mpa，注水井的井底压力设为31.0Mp。假设x，y，z各方 向上渗透率相等。
Figure 69834DEST_PATH_IMAGE112
初始为0.25，初始含油饱和度设为0.90，初始含水饱和度设为0.10。相 对渗透率曲线表示含水饱和度和相对渗透率之间的关系，如图6所示，分别示出了水相和油 相的相对渗透率曲线。模拟该油藏未来1000天内的油水流动状态，每段生产周期为50天，共 20个时间步。首先利用序贯高斯模拟生成样本数为2000网格数为40×40×20的随机三维渗 透场，利用数值模拟器计算油水两相偏微分方程计算一段时间的高保真剩余油饱和度或压 力场结果作为数据集，随机选取1600个样本作为训练集，其余400个样本作为测试集。网络 包含编码、解码和LSTM结构，具体结构如表1所示，卷积核大小均为3×3×3。设置网络超参 数初始值，初始化网络中的超参数（傅里叶卷积编码解码网络和LSTM网络权重和偏置，损失 函数中的贡献度系数等），使用ADAM优化算法，初始学习率lr=0.003，批量大小为8，混合权 重为1000。神经网络模型训练200个epoch。完成训练之后对神经网络模型进行测试。本实验 通过训练3D卷积傅里叶神经网络代理模型，输出50个时间步、每个时间步为30天的饱和度 场与压力场场图。
表 1 傅里叶卷积网络结构细节表
Figure 694851DEST_PATH_IMAGE113
如图7所示，经过200个epoch训练完成的3D卷积神经的测试平均均方误差RMSE为0.0651，与训练损失接近，并未产生过拟合。
表2为采用40×40×20网格、定压生产的三维油藏的饱和度场和压力预测结果，可 见，饱和度场和压力预测结果的RMSE值分别为均接近于0，拟合回归效果较好，精度较高。通 过公式（15）量化的饱和度的相对误差
Figure 85381DEST_PATH_IMAGE114
和压力的相对误差
Figure 897479DEST_PATH_IMAGE099
分别为5%和0.67%，这些相 对较低的误差性能表明神经网络对三维油藏模型的预测精度很高。傅里叶卷积结构在三维 空间上很好的提取了空间信息，同时LSTM结构很好的考虑了时序特征。
表 2 三维油藏饱和度或压力预测结果
Figure 149600DEST_PATH_IMAGE115
图8（生产井1和生产井2）和图9（生产井3和生产井4）展示了通过peaceman公式计算出的产量和真实产量对比结果，在整个生产周期内，预测结果和真实产量很接近，油井产量在750m3/day到1000m3/day范围之间拟合的较好，说明3D卷积傅里叶神经网络模型预测效果良好。
上述实验可以有效证明，通过本发明的方法，可以有效提高三维油藏的剩余油和压力预测速度。在网络中采用傅里叶卷积结构和LSTM结构相结合考虑时空性和物理意义，同时也考虑重力约束。本方法处理三维油藏有更好的效果，训练速度更快，精度更高泛化性更好，更具有现场实用价值。
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (6)

1.一种基于卷积傅里叶神经网络的三维地下油藏动态预测方法，其特征在于，利用构建的3D卷积傅里叶神经网络模型对三维地下油藏剩余油饱和度或压力进行动态预测，具体包括以下步骤：
步骤1、采集三维地下油藏数据，利用SGeMS地质统计工具生成符合地质特征的三维渗透率场，利用数值模拟器求解出不同渗透率对应的一段时间内的饱和度或压力作为样本库，按比例划分数据集，渗透率作为输入，饱和度或压力作为输出；
步骤2、构建3D卷积傅里叶神经网络，将3D卷积和傅里叶变换相结合，利用3D卷积算子提取图像局部信息，利用傅里叶变换提取物理信息近似微分算子，在网络中添加重力约束；利用编码解码网络实现渗透率到油藏状态的映射，结合LSTM同时考虑时间和空间信息；
步骤3、设置卷积傅里叶网络模型的超参数，在训练集下训练3D卷积傅里叶网络模型；
步骤4、使用测试集评估3D卷积傅里叶神经网络的性能，计算RMSE和量化饱和度或压力的相对误差；
步骤5、训练完成后，输出测试指标良好的3D卷积傅里叶网络模型；基于油藏安装的监测设备实时采集三维地下油藏数据，利用地质统计工具在线生成三维油藏的渗透率场，然后利用3D卷积傅里叶网络模型在线预测该油藏未来一段时间内油藏饱和度或压力场分布，通过压力和饱和度可以进一步计算产量，为制定生产策略提供直观的参考。
2.根据权利要求1所述基于卷积傅里叶神经网络的三维地下油藏动态预测方法，其特征在于，所述步骤1的具体过程如下：
步骤1.1、采集三维地下油藏数据，使用SGeMS地质统计模型生成三维渗透率场，网格尺寸为40×40×20，生成2000个渗透率样本；
步骤1.2、利用数值模拟器使用有限元方法计算三维油藏***的油水两相流动方程得到所有网格块中的饱和度或压力结果；油藏地下油水流动***是三维的，需要考虑重力效应，由质量守恒方程（1）和相达西速度方程（2）构成；
质量守恒方程表达式如（1），
Figure 226574DEST_PATH_IMAGE001
（1）
其中，
Figure 878135DEST_PATH_IMAGE002
表示梯度算子；j表示不同的相，包括油水两相，o表示油，w表示水；
Figure 251348DEST_PATH_IMAGE003
表示相 的密度；v j 表示相达西速度；
Figure 115399DEST_PATH_IMAGE004
表示源/汇项，上标l表示井处的参数；
Figure 664192DEST_PATH_IMAGE005
表示孔隙度；S j 表 示相的饱和度；t表示时间；
相达西速度方程的表达式如（2），
Figure 830862DEST_PATH_IMAGE006
（2）
其中，k rj 表示相的相对渗透率，
Figure 566737DEST_PATH_IMAGE007
表示绝对渗透率向量，
Figure 359112DEST_PATH_IMAGE008
表示相的粘度，p j 表示相的 压力，g表示重力，z表示深度；
步骤1.3、按照8：2的比例分割数据集和测试集；输入数据渗透率为(n,nx,ny,nz,1)的张量，输出数据饱和度或压力为(n,nx,ny,nz,1,T)的张量，n是样本数，nx为横向网格数，ny为纵向网格数，nz为深度方向的网格数，T为时间步数量。
3.根据权利要求2所述基于卷积傅里叶神经网络的三维地下油藏动态预测方法，其特征在于，所述步骤2的具体过程如下：
3D卷积傅里叶神经网络，主要包含编码网络和解码网络两个部分，考虑到时间序列影响增加了处理时间的循环神经模块；输入参数是已知观测值渗透率场，输出参数是未来一段时间的饱和度或压力；
编码网络结构包含M个3D卷积傅里叶层，每一个卷积傅里叶层结构如下，
Figure 762412DEST_PATH_IMAGE009
（3）
Figure 598518DEST_PATH_IMAGE010
（4）
其中，x是最初输入参数，W是权重，
Figure 821689DEST_PATH_IMAGE011
是激活函数，使用Relu激活；v N+1(x)是卷积傅里 叶模块N+1层的输出结果，v N (x)是N层的输出结果；
Figure 152177DEST_PATH_IMAGE012
积是三维傅里叶积分算子， 用以提取物理信息；
Figure 144403DEST_PATH_IMAGE013
是快速傅里叶变换，
Figure 308668DEST_PATH_IMAGE014
是傅里叶逆变换，C是截断超过最大阈值 的傅里叶级数的操作，傅里叶变换操作能通过傅里叶级数即一系列正余弦曲线去近似数据 所符合的物理规律；Conv3D是3D卷积结果，用以提取局部信息，每一层的卷积核大小和步幅 会进行调整，卷积核的大小均为3×3×3，卷积步幅有1×1×1和2×2×2两种；因为傅里叶 变换是无限维度映射的，使用傅里叶提取的结果需要与卷积结果维度一致，相加之后不会 改变原来维度，再经过五层傅里叶层，得到的最终结果进入解码网络；
解码网络结构由反卷积层构成，解码网络的部分层的输入由编码网络对应的输出矩阵和当前输入复制拼接而成，目的是保存更多信息防止卷积过程中出现的特征丢失；解码网络是将中间结果解码到原数据空间的过程，得到的结果为某一时间步的压力值或饱和度；
考虑到时间序列的影响，将多个时间步的输入经过解码网络的中间结果通过LSTM结构，结合时间序列信息，经过LSTM处理后对每个特征图用对应的解码网络得到每一个时间步的输出结果；
LSTM在t时刻的输入有三个：细胞状态C t-1，隐层状态h t-1，当前输入特征x t ；
C t 承载着空间时间信息，计算公式如下，
Figure 894502DEST_PATH_IMAGE015
(5)
其中，f t 代表遗忘门，控制丢弃细胞状态C t-1中的某些信息；
Figure 638467DEST_PATH_IMAGE016
代表新的候选细胞状态，i t 代表更新门，选择保留候选细胞状态
Figure 140992DEST_PATH_IMAGE016
中的某些信息；
隐层输出结果h t 由输出门o t 和上一步计算出的C t 经过tanh激活函数计算得出，如下，
Figure 945000DEST_PATH_IMAGE017
(6)
遗忘门f t ，更新门i t ，输出门o t ，以及新的候选细胞状态
Figure 343813DEST_PATH_IMAGE016
计算公式分别如下，每种门 有各自对应的权重和偏置，LSTM中的
Figure 891469DEST_PATH_IMAGE018
是sigmoid激活函数，激活结果在（0，1）之间，tanh 激活结果在（-1，1）之间；
Figure 982922DEST_PATH_IMAGE019
(7)
Figure 488990DEST_PATH_IMAGE020
(8)
Figure 908470DEST_PATH_IMAGE021
(9)
Figure 135183DEST_PATH_IMAGE022
(10)
其中，W f 、W i 、W o 和W c 分别为遗忘门、更新门、输出门和新的候选细胞状态对应的权重；b f 、b i 、b o 和b c 分别为遗忘门、更新门、输出门和新的候选细胞状态对应的偏置；
输入渗透率场经过卷积傅里叶编码层得到隐层变量，再经过LSTM结构和解码网络得到 同时包含时间信息和空间信息的饱和度或压力的预测结果
Figure 956508DEST_PATH_IMAGE023
。
4.根据权利要求3所述基于卷积傅里叶神经网络的三维地下油藏动态预测方法，其特征在于，所述步骤3的具体过程如下：
步骤3.1、设置3D卷积傅里叶网络模型的输入输出；输入参数是三维油藏的渗透率；设置输出值，预测未来一段时间内的油藏饱和度或压力分布；
步骤3.2、初始化超参数，批处理大小为8，激活函数使用ReLu函数，随机初始化权重W；超参数根据网络评估性能进行优化调整，损失函数使用MSE函数计算损失，训练压力时MSE损失函数值为：编码解码网络使用相同的超参数设置，训练在200个epoch内收敛；
饱和度的损失函数loss sat_MSE 如下：
Figure 758111DEST_PATH_IMAGE024
（11）
其中，y sat_i 与
Figure 664887DEST_PATH_IMAGE025
分别为训练集中第i个样本的输出饱和度真实值和3D卷积傅里叶 神经网络模型的饱和度计算结果，n train 为训练集样本总数；
对于钻过多个区块的直井，在深度d层处的井筒压力
Figure 928247DEST_PATH_IMAGE026
受深度z变化的结引力的影响； 最上面的射孔处的井筒压力称为井底压力，继续随深度往下的d+1层处深度的井筒压力
Figure 869658DEST_PATH_IMAGE027
计算如方程（12），
Figure 717529DEST_PATH_IMAGE028
是深度差，
Figure 236235DEST_PATH_IMAGE029
是流体混合密度；
当3D卷积傅里叶网络模型对压力进行预测时，考虑在深度方向上压力受到重力作用的影响，压力满足方程（12），
Figure 929384DEST_PATH_IMAGE030
（12）
压力的损失函数loss press_MSE 如下：
Figure 335089DEST_PATH_IMAGE031
（13）
其中，y press_i 与
Figure 88281DEST_PATH_IMAGE032
分别为训练集中第i个样本的输出压力真实值和3D卷积傅里叶 神经网络模型的压力计算结果；
Figure 94283DEST_PATH_IMAGE033
是井所在网格处的d+1层深度的井筒压力实际值，
Figure 591124DEST_PATH_IMAGE034
是井所在网格处的d层深度的井筒压力预测值；
Figure 618379DEST_PATH_IMAGE035
是数据约束的系数，
Figure 542473DEST_PATH_IMAGE036
是方程约束的系 数；
网络根据设定好的超参数和网络迭代结构进行正向传播；
步骤3.3、根据步骤3.2计算得到的损失误差进行反向传播，利用ADAM优化器，更新权重W，使MSE损失逐渐减小；ADAM优化算法的初始学习率lr=0.003，批量大小为8，混合权重为1000；重复这个过程，在200个epoch之后误差下降到接受范围之内，训练完成。
5.根据权利要求4所述基于卷积傅里叶神经网络的三维地下油藏动态预测方法，其特征在于，所述步骤4的具体过程如下：
对于所有的测试样本，计算相对均方根误差RMSE，量化预测饱和度或压力的相对误差；
RMSE计算公式为：
Figure 35771DEST_PATH_IMAGE037
（14）
其中，n test 为测试集样本总数；
Figure 70723DEST_PATH_IMAGE038
为模型预测结果，包括
Figure 106812DEST_PATH_IMAGE039
与
Figure 811594DEST_PATH_IMAGE040
；y i 为数值模 拟解真实值，包括y sat_i 与y press_i ；
计算三维油藏的饱和度或压力相对误差的公式如下：
Figure 667555DEST_PATH_IMAGE041
（15）
饱和度的相对误差为
Figure 630831DEST_PATH_IMAGE042
，压力的相对误差为
Figure 990269DEST_PATH_IMAGE043
，n test 是测试样本总数，n M 是网格总数，n t 是时间步数，
Figure 895645DEST_PATH_IMAGE044
是测试集中第i个样本第m个网格t时刻的饱和度的预测结果，
Figure 973323DEST_PATH_IMAGE045
是对应饱和度的真实值；
Figure 740291DEST_PATH_IMAGE046
是测试集中第i个样本第m个网格t时刻压力的预测结果，
Figure 219814DEST_PATH_IMAGE047
是对应压力的真实值；利用最大压力减去最小压力的方式进行压力归一化，
Figure 922190DEST_PATH_IMAGE048
是 第i个测试样本t时刻的最大压力，
Figure 362530DEST_PATH_IMAGE049
是第i个测试样本t时刻的最小压力。
6.根据权利要求5所述基于卷积傅里叶神经网络的三维地下油藏动态预测方法，其特征在于，所述步骤5的具体过程如下：
保存性能指标评价较好的3D卷积傅里叶神经网络模型，将训练完成的权重和超参数保存下来，利用3D卷积傅里叶神经网络模型快速预测出三维油藏的未来的饱和度或压力状态；输入是三维油藏的渗透率场，输出是一段时间t内的压力和饱和度分布，输入经过一次前向传播的计算即可得到输出；
利用压力和饱和度通过peaceman模型计算出油水井产量，计算公式如下，
Figure 542976DEST_PATH_IMAGE050
（16）
其中，
Figure 1639DEST_PATH_IMAGE051
是第m个网格中井的油/水相的注入/采油速率，k m 为第m个网格的绝对渗透 率，k r, j 为k m 对应的油/水相对渗透率，S w, m 为第m个网格处的含水饱和度值，
Figure 609338DEST_PATH_IMAGE052
是网格厚 度，
Figure 304018DEST_PATH_IMAGE053
，默认
Figure 22575DEST_PATH_IMAGE054
，油藏在x与y两个方向上宽度相等，r w 是井眼半 径，
Figure 601324DEST_PATH_IMAGE055
是水或油的粘度，P m 与P l 分别是第m个网格处压力与井眼压力。
CN202210603016.6A 2022-05-31 2022-05-31 基于卷积傅里叶神经网络的三维地下油藏动态预测方法 Active CN114693005B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210603016.6A CN114693005B (zh) 2022-05-31 2022-05-31 基于卷积傅里叶神经网络的三维地下油藏动态预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210603016.6A CN114693005B (zh) 2022-05-31 2022-05-31 基于卷积傅里叶神经网络的三维地下油藏动态预测方法

Publications (2)

Publication Number Publication Date
CN114693005A true CN114693005A (zh) 2022-07-01
CN114693005B CN114693005B (zh) 2022-08-26

Family

ID=82131298

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210603016.6A Active CN114693005B (zh) 2022-05-31 2022-05-31 基于卷积傅里叶神经网络的三维地下油藏动态预测方法

Country Status (1)

Country Link
CN (1) CN114693005B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115204531A (zh) * 2022-09-16 2022-10-18 中科数智能源科技(深圳)有限公司 基于傅里叶神经算子的油藏预测方法、设备及介质
CN115204530A (zh) * 2022-09-16 2022-10-18 中科数智能源科技(深圳)有限公司 基于傅里叶神经算子和循环神经网络的油藏预测方法
CN115618750A (zh) * 2022-12-16 2023-01-17 中国石油大学(华东) 一种基于耦合神经网络的地下油水渗流代理模型
CN117494617A (zh) * 2023-12-29 2024-02-02 中国石油大学(华东) 基于内嵌物理信息神经网络的二氧化碳驱油快速模拟方法

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5251286A (en) * 1992-03-16 1993-10-05 Texaco, Inc. Method for estimating formation permeability from wireline logs using neural networks
CN102061910A (zh) * 2010-06-30 2011-05-18 中国石油大学(北京) 砂岩油藏注水开发可预测物理模型建立方法
CN106291700A (zh) * 2016-09-28 2017-01-04 西安交通大学 基于同步挤压变换的地震加权平均瞬时频率提取方法
CN110967471A (zh) * 2019-11-14 2020-04-07 广东电网有限责任公司 一种变压器油中溶解气体浓度预测方法
CN111561313A (zh) * 2020-05-15 2020-08-21 中国地质大学（北京） 基于物理模型和机器学习的致密砂岩储层参数预测方法
CN112100930A (zh) * 2020-11-11 2020-12-18 中国石油大学(华东) 一种基于卷积神经网络和Eaton公式的地层孔隙压力计算方法
CN112541572A (zh) * 2020-11-25 2021-03-23 中国石油大学(华东) 基于卷积编码器-解码器网络的剩余油分布预测方法
CN113052371A (zh) * 2021-03-16 2021-06-29 中国石油大学(华东) 基于深度卷积神经网络的剩余油分布预测方法及装置
US20220101494A1 (en) * 2020-09-30 2022-03-31 Nvidia Corporation Fourier transform-based image synthesis using neural networks
CN114492213A (zh) * 2022-04-18 2022-05-13 中国石油大学(华东) 基于小波神经算子网络模型剩余油饱和度和压力预测方法
CN114492211A (zh) * 2022-04-15 2022-05-13 中国石油大学(华东) 一种基于自回归网络模型的剩余油分布预测方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5251286A (en) * 1992-03-16 1993-10-05 Texaco, Inc. Method for estimating formation permeability from wireline logs using neural networks
CN102061910A (zh) * 2010-06-30 2011-05-18 中国石油大学(北京) 砂岩油藏注水开发可预测物理模型建立方法
CN106291700A (zh) * 2016-09-28 2017-01-04 西安交通大学 基于同步挤压变换的地震加权平均瞬时频率提取方法
CN110967471A (zh) * 2019-11-14 2020-04-07 广东电网有限责任公司 一种变压器油中溶解气体浓度预测方法
CN111561313A (zh) * 2020-05-15 2020-08-21 中国地质大学（北京） 基于物理模型和机器学习的致密砂岩储层参数预测方法
US20220101494A1 (en) * 2020-09-30 2022-03-31 Nvidia Corporation Fourier transform-based image synthesis using neural networks
CN112100930A (zh) * 2020-11-11 2020-12-18 中国石油大学(华东) 一种基于卷积神经网络和Eaton公式的地层孔隙压力计算方法
CN112541572A (zh) * 2020-11-25 2021-03-23 中国石油大学(华东) 基于卷积编码器-解码器网络的剩余油分布预测方法
CN113052371A (zh) * 2021-03-16 2021-06-29 中国石油大学(华东) 基于深度卷积神经网络的剩余油分布预测方法及装置
CN114492211A (zh) * 2022-04-15 2022-05-13 中国石油大学(华东) 一种基于自回归网络模型的剩余油分布预测方法
CN114492213A (zh) * 2022-04-18 2022-05-13 中国石油大学(华东) 基于小波神经算子网络模型剩余油饱和度和压力预测方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
SOUMI CHAKI: "A Novel Preprocessing Scheme to Improve the Prediction of Sand Fraction From Seismic Attributes Using Neural Networks", 《IEEE XPLORE》 *
何欣: "灰岩裂缝储层地震识别技术研究——以四扣地区为例", 《中国优秀硕士学位论文全文数据库 (基础科学辑)》 *
包志强等: "快速的卷积神经网络算法及应用", 《计算机工程与设计》 *
孟宪伟等: "《各向异性油藏水平井变流量压力动态响应计算方法》", 《石油钻采工艺》 *
张凯等: "断块油藏高部位剩余油径向钻孔产能预测", 《油气地质与采收率》 *
王谦等: "径向基函数神经网络法致密砂岩储层相对渗透率预测与含水率计算", 《石油地球物理勘探》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115204531A (zh) * 2022-09-16 2022-10-18 中科数智能源科技(深圳)有限公司 基于傅里叶神经算子的油藏预测方法、设备及介质
CN115204530A (zh) * 2022-09-16 2022-10-18 中科数智能源科技(深圳)有限公司 基于傅里叶神经算子和循环神经网络的油藏预测方法
CN115204531B (zh) * 2022-09-16 2022-12-27 中科数智能源科技(深圳)有限公司 基于傅里叶神经算子的油藏预测方法、设备及介质
CN115618750A (zh) * 2022-12-16 2023-01-17 中国石油大学(华东) 一种基于耦合神经网络的地下油水渗流代理模型
CN117494617A (zh) * 2023-12-29 2024-02-02 中国石油大学(华东) 基于内嵌物理信息神经网络的二氧化碳驱油快速模拟方法
CN117494617B (zh) * 2023-12-29 2024-04-16 中国石油大学(华东) 基于内嵌物理信息神经网络的二氧化碳驱油快速模拟方法

Also Published As

Publication number Publication date
CN114693005B (zh) 2022-08-26

Similar Documents

Publication Publication Date Title
CN114693005B (zh) 基于卷积傅里叶神经网络的三维地下油藏动态预测方法
CN113052371B (zh) 基于深度卷积神经网络的剩余油分布预测方法及装置
CN112541572B (zh) 基于卷积编码器-解码器网络的剩余油分布预测方法
CN114492213B (zh) 基于小波神经算子网络模型剩余油饱和度和压力预测方法
CN111523713A (zh) 一种预测油田中剩余油饱和度分布的方法和装置
US20230350096A1 (en) System and method for effective hydrocarbon reservoir pressure prediction and control
Bai et al. Accelerating geostatistical modeling using geostatistics-informed machine Learning
CN106355003A (zh) 基于t分布的马尔科夫链蒙特卡洛自动历史拟合方法及***
Jiang et al. Use of multifidelity training data and transfer learning for efficient construction of subsurface flow surrogate models
Gao et al. Modeling of meandering river deltas based on the conditional generative adversarial network
CN112541256A (zh) 基于深度学习降维重构的强非均质油藏历史拟合方法
Han et al. Multiple point geostatistical simulation with adaptive filter derived from neural network for sedimentary facies classification
Shokir et al. Permeability estimation from well log responses
Chen et al. Physics-guided meta-learning method in baseflow prediction over large regions
Simonov et al. Application of machine learning technologies for rapid 3d modelling of inflow to the well in the development system
Zhang et al. Efficient surrogate modeling based on improved vision transformer neural network for history matching
Razak et al. Embedding physical flow functions into deep learning predictive models for improved production forecasting
Wen et al. Multiphase flow prediction with deep neural networks
Le et al. Physics-constrained deep learning for production forecast in tight reservoirs
Ju et al. Hydrologic simulations with artificial neural networks
Wang et al. On the Feasibility of An Ensemble Multi-Fidelity Neural Network for Fast Data Assimilation for Subsurface Flow in Porous Media
Costa et al. A new methodology to reduce uncertainty of global attributes in naturally fractured reservoirs
GB2584449A (en) Apparatus method and computer-program product for calculating a measurable geological metric
Silva et al. Generative network-based reduced-order model for prediction, data assimilation and uncertainty quantification
Azizzadeh Mehmandoust Olya et al. Hydrocarbon reservoir potential mapping through Permeability estimation by a CUDNNLSTM Deep Learning Algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant