秦 峰，闫正和，唐圣来，罗睿乔，龚 斌
（1.中海石油（中国）有限公司深圳分公司，广东 深圳 518000；
2.中国地质大学（武汉）资源学院，湖北 武汉 430074）
数值模拟技术目前已是大型油气藏开发决策常用的辅助分析手段之一［1-5］，然而由于控制方程的非线性性质和地质描述的多尺度特性，油气藏模型包含了复杂的渗流机理和巨量网格，当模型应用于优化、不确定性量化和数据同化时可能需要运行数千次模拟计算［6-9］，这样模型的计算成本极高。应用基于深度学习的代理模型进行油气藏模拟是油气藏仿真研究的一个新方向［10-20］，代理模型比高保真油气藏模型的模拟运算速度大幅提升，因此计算成本更加低廉，在实际应用中可替代高保真油气藏模型，实现大型复杂油气藏开发方案的快速迭代。
目前国外对构建地下流体流动模拟的代理模型已进行了广泛的研究，其中基于图像处理网络架构的代理模型发展较快，这些构建代理模型的方法已被应用于地质参数化、不确定性量化和代理/降阶建模等方面。CANCHUMUNI 等将基于图像处理的网络架构应用于地质建模，使用变分自动编码器（VAE，Variational Auto Encoder）将随机化的低维隐藏变量生成新的地质实现［21-22］。在流体流动模拟方面，WANG 等使用长短期记忆（LSTM，Long-Short Term Memory）循环神经网络模型来逼近由本征正交分解（POD，Proper-Orthogonal-Decomposition）方法所建立的流动动力学模型［23-24］。GONZALEZ 等用自编码器（AE，Auto-Encoder）架构代替了POD 步骤进行低维表示［25］。TANG 等引入一种基于深度学习的代理模型R-U-Net［26］，该模型同时包含了卷积和递归神经网络，在井控条件不变的情况下用于预测不同地质模型的流量响应，在油气藏辅助历史拟合过程中应用较好。TEMIRCHEV 等构造了一个与降阶模型相似的模型，用VAE 来表示低维的储层状态［27］。模型结合线性回归、LSTM 或门控循环单元（GRU）进行了实际油藏动态模拟测试，但结果显示在不考虑井控条件的情况下，模型结果的相对误差在某些验证方案中相对较大。上述模型只适用于油气藏内边界条件不变的情况，即没有考虑井控条件的变化，而实际油气藏开采过程中井控条件一般都会随时间发生变化，因此上述模型对油气藏实际流动模拟问题适应性较差。研究前期许多模型都是端到端的纯数据驱动方法，没有考虑流体流动的基本控制方程。为了提高模型性能，最新进展中已有许多方法将物理约束纳入到深度神经网络模型构建过程中。MO 等将VAE 模型扩展到处理多相流问题，并通过引入额外的物理约束进一步提高性能［28］。RAISSI 等提出了使用密集连接前馈神经网络的物理知识深度学习框架（被称为物理信息神经网络或PINN）［29］。在PINN中，将与控制偏微分方程相关联的残差函数引入神经网络的损失函数中，并利用控制方程的残差构建了用于不确定性量化的代理模型。
笔者采用嵌入式控制框架（E2C，Embed to Con⁃trol）模型构建流动控制方程的降阶模型，并对模型架构进行了调整优化，将编码的卷积输出结果加入到解码的卷积输入中，这样解码器可以获得更多的油气藏属性空间信息；
同时将油气藏井控条件作为物理约束加入到模型中，因此该模型适用于油气藏边界条件可变的情况，更加符合油气藏实际流动模拟状态。
达西定律描述了在孔隙介质中流体的宏观流动方式，因此，采用基于达西定律的常规油气水三相流动模型。从各组分的质量守恒出发，结合达西定律并忽略毛管压力的影响，导出非混相油气水流动的控制方程为：
其中：
（1）式中，ϕ和K表示地质特征，Sj和p为求解变量。一般采用有限体积方法对（1）式求解，并针对每个网格块计算它们的解。在本次研究中，使用特雷西能源科技有限公司的油气藏数值模拟器Cloud［30］来构建所有训练数据。
WATTER 等提出了E2C 模型，在机器人规划系统的背景下将直接感觉数据（图像）和时变控制作为输入来预测系统状态的演化［31］。E2C 模型将VAE 模型和生成模型相结合，VAE 模型是将系统状态投影到低维子空间的推理模型，生成模型是对预测结果进行重构生成的模型和线性转化模型。后者基于时变控制输入近似低维状态的演化。借鉴E2C模型的设计思路，开发一个包含物理约束（井控条件）的嵌入式控制框架的降阶模型，主要包括3个处理部分：将系统变量从高维空间投射到低维子空间（这里称为隐藏空间）的编码器、低维近似系统动力学的线性转化模型、将解映射回高维（全阶）空间的解码器或生成模型。该模型与WATTER 所设计模型的区别在于对模型架构进行了调整优化，将编码的卷积输出结果加入到解码的卷积输入中，这样解码器可以获得更多的油气藏属性空间信息；
同时将油气藏井控条件作为物理约束加入到模型中，更适用于油气藏实际流动模拟。
2.1 编码器
编码器提供了全序状态变量的低维表示。在训练时编码器可简单地表示为：
E2C 模型中编码器结构如图1 所示。这里使用的编码器网络由3 个卷积块、3 个残差块和1 个全连接层堆积组成。在训练过程中，压力和饱和度被输入到编码器网络，隐藏状态的特征矩阵zt被生成。
图1 编码器结构Fig.1 Structure of encoder
2.2 线性转化模型
线性转化模型根据给定的控制将隐藏变量从一个时间步发展到下一个时间步，用于描述模型由于外界作用而发生变化的过程，文中主要用于描述模型由于井控条件的变化导致压力场、饱和度场的变化过程。线性转化模型（图2）由3 个转化块和2个全连接层组成。线性转化模型遵循一个Densebatchnorm-ReLU 架构，这被认为是全连接网络的标准处理步骤。该模型的输入包括：隐藏状态的特征矩阵zt、井控条件ut+1和预测的时间间隔Δt。模型输出为下一时刻隐藏状态的特征矩阵ẑt+1。
图2 线性转化模型结构Fig.2 Structure of linear converter
2.3 解码器
解码器的结构类似于编码器，其结构如图3 所示，具体可以表示为：
解码器结构如图3所示。将线性转化模型的输出zt+1与编码器的卷积输出Et结合作为解码器的输入，预测高维状态的输出，即油气藏的压力场及饱和度场。解码器结构与编码器结构类似，只是顺序是颠倒的，由1 个全连接层、3 个残差块、4 个逆卷积块和1个卷积层堆叠组成。密集层将一个低维的潜在向量转换为特征映射的堆积，这些特征图先通过残差块和逆卷积块转化得到多个特征图，再通过卷积层将多个特征映射转换为压力和饱和度场作为最终输出。
图3 解码器结构Fig.3 Structure of decoder
2.4 模型实现及训练细节
E2C模型的输入数据为数值模拟算例中包含的压力场、饱和度场及时序类型的井控条件数据。对输入数据的预处理主要是对数据进行归一化处理，可以提高学习过程和输出质量［32］。饱和度S∈[0,1]不需要归一化，井控条件的压力数据需要归一化。标准压力p0（包括网格压力和井底流压）的表达式为：
E2C模型的输出数据为压力场、饱和度场，为了便于实际应用，在计算出压力场、饱和度场结果的基础上，利用皮斯曼方程计算单井的产量。关于多少个训练样本可以使模型达到最佳训练效果，即使模型达到精准度与训练时间的平衡还没有进一步的深入研究，文中使用了数值模拟结果的500 个算例作为训练数据，模型训练后达到了预期效果。
3.1 区块概况及算例
对实际气藏进行三维模拟以测试E2C 模型效果。番禺35-1气田位于珠江口盆地流花区块，区域上位于白云凹陷北坡，北部紧靠番禺低隆起。气田所在海域水深约为195～202 m，整体为受岩性控制的层状边水砂岩气藏。
利用特雷西能源科技有限公司的油气藏数值模拟器Cloud 建立气水模拟的算例，作为E2C 模型的训练集。具体算例的网格模型（图4）包含3 口生产井，其中2 口直井A2 和A3，一口水平井A1H，实际生产为衰竭式开发。网格块的尺寸为50 m×50 m×0.5 m，总网格数为154×81×124 共1 546 776。模型渗透率的算术平均值为25.4 mD，孔隙度平均值为0.147。相对渗透率的表达式为：
图4 番禺35-1气田算例的网格模型Fig.4 Grid model of PY35-1 Gas Field case
（6）式和（7）式中，K0ro为1.0，K0rw为0.7，Sor为0.3，Swr为0.1，a为3.6，b为1.5。其他流体参数包括：ρg为0.738 kg/m3，ρw为1 000 kg/m3，μg为0.02 mPa·s，μw为0.91 mPa·s，毛管压力的影响可以忽略。
模型采用了来自于研究区1 400 d 的实际生产数据。模拟时间间隔对于训练和测试运行都是相同的。生产井采用变井底流压控制，井底流压数据按实际输入。
3.2 应用效果
通过对比E2C 模型与Cloud 模拟器模拟结果的差异及模拟用时，评价E2C 模型的实际应用效果。在相同的生产时间节点，对比E2C 模型与Cloud 模拟器模拟的产量曲线、饱和度场、压力场，并定量分析其差异。
产量曲线对比 分析实际工区3口生产井的井底流压曲线（图5），其井底流压为0~32 MPa，Cloud模拟器和E2C 模型在相同的井底流压曲线即井控条件的约束下进行运算，E2C 模型在训练阶段的日产气量曲线整体与Cloud 模拟器差别很小（图5），日产气量曲线的整体误差绝对值为0~0.5 Sm3/d，日产水量曲线的整体误差绝对值为0~1.2 Sm3/d，E2C 模型训练效果理想。预测阶段Cloud 模拟器和E2C 模型都在井底流压为25 MPa 条件下定压生产，E2C 模型在预测阶段误差稍微增加，日产气量曲线的整体误差绝对值为0~0.95 Sm3/d，日产水量曲线的整体误差绝对值为0~3.6 Sm3/d。由此可见，E2C 模型模拟的生产曲线较Cloud 模拟器差别不大，曲线形态与Cloud模拟器匹配较好。
图5 番禺35-1气田3口生产井的日产量曲线对比Fig.5 Comparison of daily production curves of three producing wells in PY35-1 Gas Field
饱和度场对比 以番禺35-1 气田生产1 400 d的饱和度场（图6）为例说明。由于Cloud 模拟器模拟结果（图6a）和E2C 模型模拟结果（图6b）从视觉上看是相似的，因此对两者之间的差异进行了定量描述（图6c），可以看出偏差基本为-0.05～0.01，平均相对误差小于5%，饱和度场模拟效果较好。
图6 番禺35-1气田饱和度场对比及误差分析Fig.6 Comparison of saturation fields of PY35-1 Gas Field and error analysis
压力场对比 以番禺35-1 气田生产1 400 d 的压力场（图7）为例说明。由于Cloud 模拟器模拟结果（图7a）和E2C 模型模拟结果（图7b）从视觉上看是相似的，对两者之间的差异进行了定量描述（图7c），E2C 模型局部位置的压力场与Cloud 模拟器模拟的压力场存在差异，且Cloud 模拟器计算的压力场更加光滑连续，分析原因主要是基于物理渗流过程的Cloud 模拟器可以再现渗流过程，基于数据驱动的E2C 模型对于某些渗流过程的细节特征难以提取与保留，这也是未来E2C 模型的改进方向。与Cloud模拟器模拟结果相比，E2C模型模拟的压力场平均相对误差为8%，效果较好。
图7 番禺35-1气田压力场对比及误差分析Fig.7 Comparison of pressure fields of PY35-1 Gas Field and error analysis
运行时间对比 在Intel Xeon ES-2670 双CPU（24 核）的计算节点上，Cloud 模拟100 个测试算例，每个耗时60 s，运行用时共计6 000 s。在分配了8 GB 内存的Tesla V100 GPU 节点上，使用E2C 模型可以在大约1.25 s 内评估完100 个测试用例。这种条件下进行比较可以得到6 000/1.25=4 800倍的加速，但是这个比较没有考虑GPU 和CPU 的差别。在相同的CPU（Intel Xeon ES-2670）条件下，E2C 模型的100 次测试算例运行时间大约需要16 s，这相当于6 000/16=375倍的加速，虽然比GPU的加速要小，但仍然非常可观。
采用一种基于深度学习的嵌入式控制框架模型模拟油气藏地下三维属性场的变化，探索了深度学习模型在油气藏模拟方面的应用。利用E2C 模型对实际气藏进行三维模拟，评价了深度学习模型在真实三维气藏条件下的性能和适用性。算例结果显示：在准确度上，E2C模型可以很好地模拟气藏在不同的井控条件下压力场、饱和度场的变化，与Cloud 模拟器结果对比误差较小；
在时间性能上，E2C 模型模拟耗时是传统方法的1/375。综上所述，E2C模型在优化和历史拟合方面比传统数值模拟方法更具潜力。
研究中实际应用的E2C 模型预测耗时不包括前期数值模拟模型准备、批量计算、代理模型训练等耗时，如果综合考虑全部耗时，E2C模型的时间优势仍然比较明显：第一，前期模型的准备、批量计算、代理模型训练等工作可以在机器空闲时进行，即线下进行；
训练好的模型由于具备极快的求解速度，可以在一些智能油藏的场景，比如生产方案实时优化、实时决策等，无延迟完成，即线上完成。第二，通过相对有限的训练算例得到的代理模型，具备在更广泛参数空间采样后快速求解的能力，因此在运用智能优化算法（如粒子群算法）时，能够避免成百上千优化算例的直接求解，这部分时间的节省也非常显著。模型泛化性方面，本次实际应用的E2C 模型主要针对具体的一个气藏，目前的研究进展表明，代理模型的实际应用在现阶段还不具备泛化性，需要针对具体的地质模型和井位部署、在一定边界的生产制度下完成替代训练。目前正在探索具有泛化能力的网络模型路线，但需要建立极其海量的训练集（不同区块、不同地质模型的多个实现，以及不同生产制度下的多个算例），经过超大型的训练方能完成。
符号解释
