融合点云残差网络和考虑滑移的流场和压力场预测方法

1.本发明属于智能油气田开发技术领域，具体涉及一种融合点云残差网络和考虑滑移的流场和压力场预测方法。


背景技术：

2.油气资源事关国家社会经济可持续发展和能源安全，随着常规油气资源的开发进入中后期，我国石油对外依存度逐年上升，我国页岩油储量丰富，是未来增储上产的主要阵地。页岩油藏孔隙介质中发育着丰富的有机孔隙和无机孔隙，页岩油在这两种不同孔隙内的运移规律以及赋存状态均不相同，纳米级孔隙的壁面矿物属性对页岩油流动有极大的影响不同矿物壁面具有不同的润湿性，流体在壁面的吸附程度也不同。有机质孔隙表面为干酪根，油相润湿；无机质孔隙表面为黏土矿物，水相润湿。此外，在不同矿物相孔隙内的流动具有不同的滑移条件，纳米孔隙中的有机质壁面滑移对流体流动的影响很大，无机质的壁面滑移几乎可以忽略，常规基于达西定律的渗流数学模型不再适用。开发技术政策制定、开发方案优化、提高采收率技术实施的基本前提是探索出页岩储层中流体的渗流机理，以及明确致密油有效动用条件与流动参数。并且通过岩心实验来实现揭示致密多孔隙的渗流汇率存在耗时长、成本高等问题。
3.到目前为止，针对页岩油流动模拟有一系列方法，如分子动力学模拟(mds)、孔隙网络建模(pnm)和直接数值模拟(dns)、岩心分析等。mds需要大量计算，其模拟规模太小。因此，直接使用mds在孔隙尺度上模拟多孔介质中的流动是不现实的。pnm是一种具有孔隙和喉道的简单几何形状，用于表示多孔介质中的复杂孔隙空间，其优点是反映数字岩心的拓扑和连通性，并减少模拟计算然而，由于孔隙与喉道空间的简化，pnm进行流动模拟存在一定的误差。对于dns虽然求解精度较高，但实际流体动力学问题的传统数值模拟需要在特定模型参数的网格上求解偏微分方程(pde)系统。求解非线性耦合偏微分方程的求解过程通常需要构造大型矩阵，这使得迭代过程计算复杂。
4.深度学习技术近年来随着机器学习和人工智能的诞生和发展逐渐兴起，具有自我组织学习输出特征，在学习中不断适应问题的特性，能够十分稳定并且准确地去分析处理大规模的数据集，深度学习在非线性计算、模式识别、数据挖掘等领域都被大量应用，其中的卷积神经网络方法也越来越广泛地被用于石油行业中。因此，可以利用卷积神经网络来解决流动模拟需要求解大量的耦合pde问题，利用卷积网络可以得到多孔介质的几何边界，在网络中控制偏微分方程的残差，将边界条件等加入损失函数中。通过一定的稀疏点来进行训练，利用卷积网络求解逆问题，从而得到整个多孔介质考虑滑移流动的速度场和压力场。


技术实现要素：

5.为了解决上述问题，本发明提出了一种融合点云残差网络和考虑滑移的流场和压力场预测方法，通过将点云残差网络中的损失函数加入流动方程，选取一定的稀疏点作为
训练集，通过稀疏点求解逆问题，得到多孔介质的压力场、速度场，以此来实现页岩油多孔介质考虑滑移的流动精准预测。
6.本发明的技术方案如下：
7.一种融合点云残差网络和滑移的流场和压力场预测方法，具体包括如下步骤：
8.步骤1、构建三维数字岩心，并在数字岩心内随机生成孔隙结构，数字岩心的孔隙结构内充填流体，其余部分均设置为骨架；
9.步骤2、对构建的三维数字岩心进行竖直切片得到二维多孔介质，对二维多孔介质进行图像分割，得到二维多孔介质孔隙空间；
10.步骤3、利用多物理场仿真模拟软件分别求取二维多孔介质孔隙空间所对应的速度矢量，基于二维多孔介质孔隙空间及其所对应的速度矢量，构建二维多孔介质数据集；
11.步骤4、构建点云残差网络，并耦合流动方程构建损失函数；
12.步骤5、数据归一化处理；
13.步骤6、采用归一化后的数据对点云残差网络模型进行训练，设置学习率、参数epoch、batch，训练结束后，采用训练完成的模型对二维多孔介质流场和压力场进行预测，当预测结果表明当前模型性能良好时，输出该模型，将该模型用于真实场景下流场和压力场的预测。
14.进一步地，步骤1中，孔隙结构包括有机质孔隙和无机质孔隙，针对三维数字岩心划分有机质孔隙与无机质孔隙边界。
15.进一步地，步骤2中，采用image j软件对二维多孔介质片图像上的有机质、无机质和孔隙进行标注后，再将二维岩心切片的几何单位更改为纳米级别，获得边长为100nm的二维多孔介质孔隙空间。
16.进一步地，步骤3的具体过程如下：
17.步骤3.1、将二维多孔介质内的流体设置为稳态流动的不可压缩牛顿流体，流体流速不随时间而变化，则ns方程简化为等式(1)、(2)所示；设置二维多孔介质的边界条件以及有机孔隙与无机质壁面滑移条件，考虑有机孔隙壁面滑移条件如式(3)所示；设置二维多孔介质内流体的密度和粘度，将二维多孔介质的一侧设置为流体入口，将多孔介质的另一侧设置为流体出口，设置流体的入口压力和出口压力，入口与出口的初始压力条件如(4)和(5)；
[0018][0019][0020][0021]
p
in
＝c1ꢀꢀ
(4)
[0022]
p
out
＝c2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0023]
其中，μ为流体的黏度，单位为pa
·
s；表示梯度算子；为对矢量求散度；u为速度矢量，单位为m/s；p为流体的压力，单位为pa；u|
wall
为壁面的速度；b为有机孔隙壁面滑移长度，单位为m；n为壁面的法向单位向量；p
in
和p
out
分别为入口和出口的压力，单位为pa；c1和c2均为设定的压力常数；
[0024]
步骤3.2、分别针对二维多孔介质孔隙空间，基于数字岩心孔隙空间构建二维多孔
介质，利用多物理场仿真模拟软件模拟得到流体在二维多孔介质内的流动过程，获取各多孔介质稳定状态下的流场分布，获取流场流动模拟过程中各数据采集点的横坐标、纵坐标、压力、速度水平分量u和速度竖直分量v，进而构建二维多孔介质数据集；
[0025]
步骤3.3、在构建的二维多孔介质数据集内选取部分点作为数据采集点，构建训练集和测试集；
[0026]
为了保证训练集对整个平面的均匀取样，首先设置速度水平分量u和速度竖直分量v的大小阈值分别为u《10-20
,v《10-20
得到无机质边界，并标记无机质边界区域点，统计无机质边界点个数为m2，无机质壁面为无滑移条件；自定义设置点密度k，然后设置x1,y1，具体公式如下式，
[0027][0028]
其中，k为自定义的点密度，i为一个循环，i∈[0,k)且为整数；
[0029]
通过公式(6)和(7)得到k个x1和k个y1，对得到的x1和y1相互组合得到坐标(x1,y1)，由此得到k2个点，统计这k2个点在孔隙空间及数据集的点个数记为m3，并标记为内部区域点；由上述操作无机质壁面边界区域点与内部区域点作为数据采集点，构建训练和测试所用的数据集，且数据采集点的个数满足m1＝m2+m3。
[0030]
进一步地，步骤4的具体过程如下：
[0031]
步骤4.1、构建的点云残差网络共有15个一维卷积模块conv1d、1个最大池化模块maxpooling和1个拼接模块concat；
[0032]
点云残差网络以离散的采样点为输入，由于输入数据中的采样点没有先后顺序，因此首先通过两个一维卷积模块和一个点云残差模块进行校正；然后进行特征转换，通过3个一维卷积模块和3个点云残差模块，1个最大池化层maxpooling，1个concat；之后，经由4个一维卷积层模块和3个点云残差模块，其中这3个点云残差模块都是在卷积输出时与前一个卷积层的输出进行跳跃连接；在经过最后一个一维卷积层的输出后获得每个采样点在速度水平分量u、速度竖直分量v和压力p上的预测值，卷积模块激活函数均为tanh；
[0033]
步骤4.2、结合流动方程构建损失函数；损失函数包括设置数据采集点的连续性方程损失r1、数据采集点速度水平分量u的ns方程损失r2、数据采集点速度竖直分量v的ns方程损失r3、数据采集点中无机质边界区域点的速度竖直分量v的均方损失r4、数据采集点中内部区域点u和v的均方损失r5、数据采集点的压力p的均方损失r6，具体公式如下：
[0034][0035][0036]
[0037][0038][0039][0040]
loss＝100
×
r1+r2+r3+r4+r5+r6ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0041]
其中，表示微分算子，其中m1表示数据采集点的个数，m2表示数据采集点中无机质边界区域点的个数，m3表示数据采集点中内部区域点的个数，uk、vk、pk表示速度水平分量u、速度竖直分量v和压力p的预测值，xk表示点的横坐标，yk表示点的纵坐标，u
t
、v
t
、p
t
表示速度水平分量u、速度竖直分量v和压力p的真实值；loss为总损失函数，对r1、r2、r3、r4、r5、r6进行自适应参数调整。
[0042]
进一步地，步骤5的具体过程如下：
[0043]
对数据集中所有的u,v,p的值归一化到[-1,1]，使用下一个10的最高幂作为规范化常数，归一化方法如下式(15)-(16)所示，对于任何一组大小为l的数据，都能找到一个整数a满足以下不等式的最小整数；对a＝a1的速度场所有的u、v数据集进行规范化处理，a＝a2的压力场p规范化处理；对空间坐标x、y数据归一化的参数a3＝a1，对于黏度μ标准化的a4＝a2；
[0044][0045][0046]
其中，αi表示归一化前的数据，l表示数据的总量，a表示归一化参数，表示归一化后的数据。
[0047]
本发明所带来的有益技术效果：
[0048]
本发明构建了一种含有机质分布的岩心，建立了点云残差网络并耦合流动方程等作为损失函数。通过对稀疏数据集的训练以及pinn的自动微分，从而得到整个岩心的流场和压力场，可以大幅度减少计算开销。因此本发明对于油气渗流以及计算流体力学考虑滑移流动的模拟有着至关重要的作用。本发明点云残差网络预测精度更高。
附图说明
[0049]
图1为本发明融合点云残差网络和考虑滑移的流场和压力场预测方法的流程图；
[0050]
图2为本发明含有机质与有机质的二维多孔介质图；
[0051]
图3为本发明点云残差网络的结构图；
[0052]
图4为本发明采用模型进行测试时速度水平分量u的预测值和真实值的对比图；
[0053]
图5为本发明采用模型进行测试时速度竖直分量v的预测值和真实值的对比图；
[0054]
图6为本发明采用模型进行测试时压力p的预测值和真实值的对比图；
[0055]
图7为本发明距入口50nm的数据点的速度水平分量u的预测值和真实值的对比图；
[0056]
图8为本发明距入口50nm的数据点的速度竖直分量v的预测值和真实值的对比图；
[0057]
图9为本发明距入口50nm的数据点的压力p的预测值和真实值的对比图。
具体实施方式
[0058]
下面结合附图以及具体实施方式对本发明作进一步详细说明：
[0059]
点云指输入是一组点的空间坐标，将网格顶点作为点云。残差网络resnet是经典卷积神经网络中的一种结构，使用一种称为跳跃连接的技术。该技术跨过了几个层次，直接连接到输出。残差网络是经典卷积神经网络中的一种结构，使用一种称为跳跃连接的技术。该技术跨过了几个层次，直接连接到输出。残差网络不是学习底层映射的层，而是允许网络适应残差映射，残差网络的跨层连接能够有效的避免梯度消失或精度饱和等问题。残差网络能够有效的避免梯度消失或精度饱和等问题。点云残差网络能够学习空间位置与速度场、压力场之间的关系。
[0060]
滑移是指在固壁上(边界上)的切向速度不为零，即产生滑移速度。一般来讲，按照固壁与流体之间的滑移速度，可以将边界条件分成三类，即无滑移、部分滑移和完全滑移。无滑移指的是，流体在固壁上(边界上)的速度为0，即无滑移速度；部分滑移指的是，边界切向速度成一定梯度；完全滑移指的是，流体在固壁上切向速度与流体在流场中沿固壁切向速度相等。
[0061]
如图1所示，一种融合点云残差网络和考虑滑移的流场和压力场预测方法，具体步骤包括：
[0062]
步骤1、构建三维数字岩心并在数字岩心内随机生成孔隙结构，孔隙结构包括有机质孔隙和无机质孔隙，数字岩心的孔隙结构内充填流体，其余部分均设置为骨架。针对三维数字岩心划分有机质孔隙与无机质孔隙边界。
[0063]
步骤2、对构建的三维数字岩心进行竖直切片得到二维多孔介质，对二维多孔介质进行图像分割，得到二维多孔介质孔隙空间。具体过程如下：
[0064]
本实施例中，对三维数字岩心进行竖直切片得到如图2所示的二维多孔介质。利用image j软件对二维多孔介质上的有机质、无机质和孔隙进行标注后，深黑色代表无机质，浅黑色代表有机质。再将二维多孔介质几何单位更改为纳米级别，获得边长为100nm的二维多孔介质孔隙空间。
[0065]
步骤3、利用多物理场仿真模拟软件分别求取孔隙空间所对应的速度矢量(包括垂直和水平速度矢量)，基于二维多孔介质孔隙空间及其所对应的速度矢量，构建二维多孔介质数据集。具体过程如下：
[0066]
步骤3.1、将二维多孔介质内的流体设置为稳态流动的不可压缩牛顿流体，流体流速不随时间而变化，则ns方程可以简化为等式(1)、(2)所示；设置二维多孔介质的边界条件以及有机孔隙壁面滑移条件，考虑有机孔隙壁面滑移条件如式(3)所示；设置二维多孔介质孔隙流体的密度和粘度，将二维多孔介质孔隙的一侧设置为流体入口，将孔隙的另一侧设置为流体出口，设置流体的入口压力和出口压力，入口与出口的初始压力条件如(4)和(5)；
[0067][0068]
[0069][0070]
p
in
＝c1ꢀꢀ
(4)
[0071]
p
out
＝c2ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0072]
其中，μ为流体的黏度，单位为pa
·
s；表示梯度算子；为对矢量求散度；u为速度矢量，
[0073]
单位为m/s；p为流体的压力，单位为pa；u|
wall
为壁面的速度；b为有机孔隙壁面滑移长度，单位为m；n为壁面的法向单位向量；p
in
和p
out
分别为入口和出口的压力，单位为pa；c1和c2均为设定的压力常数。
[0074]
步骤3.2、分别针对二维多孔介质孔隙空间，基于数字岩心孔隙空间构建二维多孔介质，利用多物理场仿真模拟软件模拟得到流体在二维多孔介质孔隙内的流动过程，获取各孔隙区域稳定状态下的流场分布，获取流场流动模拟过程中各数据采集点的横坐标、纵坐标、压力、速度水平分量u和速度竖直分量v，进而构建二维多孔介质数据集。
[0075]
步骤3.3、在构建的二维多孔介质数据集内选取部分点作为数据采集点，构建训练集和测试集。
[0076]
为了保证训练集对整个平面的均匀取样，首先设置速度水平分量u和速度竖直分量v的大小阈值分别为u《10-20
,v《10-20
得到无机质边界，并将点标记为蓝色，蓝色点为无机质边界区域点，统计无机质边界点个数为m2，无机质壁面为无滑移条件。自定义设置点密度k，然后设置x1,y1，具体公式如下式，
[0077]
x1＝50
×
(i
÷
k+0.001)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(6)
[0078]
y1＝50
×
(i
÷
k+0.001) (7)
[0079]
其中，k为自定义的点密度，i为一个循环，i∈[0,k)且为整数。
[0080]
通过公式(6)和(7)可以得到k个x1和k个y1，对得到的x1和y1相互组合得到坐标(x1,y1)，由此可得到k2个点，统计这k2个点在孔隙空间及数据集的点个数记为m3，并标记为橙色，橙色点记为内部区域点。由上述操作无机质壁面蓝色点与内部橙色点作为数据采集点，构建训练和测试所用的数据集。且数据采集点的个数满足m1＝m2+m3。
[0081]
本实施例设置入口初始压力为0.05pa，出口初始压力为0pa，流体密度和黏度分别为1000kg/m3、0.001pa
·
s，有机质的滑移长度为130nm，无机质的滑移长度为0nm。自定义点密度k为30。经过采集及计算后，可以得到速度矢量u的水平分量u的速度场图、竖直分量v的速度场图、压力场图。然后图中选取共594个点作为数据采集点为训练集，剩余的点即为测试集。
[0082]
步骤4、构建点云残差网络，并耦合流动方程构建损失函数。具体过程为：
[0083]
步骤4.1、如图2所示，构建的点云残差网络共有15个一维卷积模块conv1d、1个最大池化模块maxpooling和1个拼接模块concat，图中方框中的数字第一个代表通道数，第二个代表数据的序列长度(如trans3中的64
×
64表示该64个通道数与64个序列长度)。
[0084]
点云残差网络以离散的采样点为输入，由于输入数据中的采样点没有先后顺序，因此首先通过两个一维卷积模块(trans1、trans3)和一个点云残差模块trans2(在卷积输出时与上一个卷积层trans2进行跳跃连接)进行校正，解决旋转不变性问题。然后进行特征转换，主要通过3个一维卷积模块(global1、global2、gloal3)和3个点云残差模块
(global1_1、global2_2、glabal3_3)，1个最大池化层maxpooling，1个concat(进行数据拼接映射到高维空间)。global1_1、global2_2在卷积输出时与前一个卷积层的输出进行跳跃连接，而glabal3_3是在concat前卷积层trans3的输出进行跳跃连接。这些特征转换的操作将二维平面中的每一个点都映射到高维空间，由此得到的每一个二维坐标点都由一个1024维的向量表示，然后通过最大池化减少冗余信息，解决坐标点的无序性问题。之后，经由4个一维卷积层(layer1、layer2、layer3、layer4)模块和3个点云残差模块(layer1_1、layer2_1、layer3_1)，其中这3个点云残差模块都是在卷积输出时与前一个卷积层的输出进行跳跃连接。在经过最后一个一维卷积层layer4的输出后可获得每个采样点在速度水平分量u、速度竖直分量v和压力p上的预测值，卷积模块激活函数均为tanh。
[0085]
步骤4.2、结合流动方程构建损失函数。损失函数包括设置数据采集点的连续性方程损失r1、数据采集点速度水平分量u的ns方程损失r2、数据采集点速度竖直分量v的ns方程损失r3、数据采集点中无机质边界区域点的速度竖直分量v的均方损失r4、数据采集点中内部区域点u和v的均方损失r5、数据采集点的压力p的均方损失r6，具体公式如下：
[0086][0087][0088][0089][0090][0091][0092]
loss＝100
×
r1+r2+r3+r4+r5+r6ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0093]
其中，表示微分算子，其中m1表示数据采集点的个数，m2表示数据采集点中无机质边界区域点的个数，m3表示数据采集点中内部区域点的个数，uk、vk、pk表示速度水平分量u、速度竖直分量v和压力p的预测值，xk表示点的横坐标，yk表示点的纵坐标，u
t
、v
t
、p
t
表示速度水平分量u、速度竖直分量v和压力p的真实值。
[0094]
为了避免速度场和压力场都收敛到0，对总损失函数loss进行了式(14)的操作，对r1、r2、r3、r4、r5、r6进行了自适应参数调整。上述公式(8)-(14)实现了pinn的自动微分。
[0095]
步骤5、数据归一化处理；具体过程如下：
[0096]
为了更好地输入、输出并训练点云残差网络，对数据集中所有的u,v,p的值归一化到[-1,1]，使用下一个10的最高幂作为规范化常数，归一化方法如下式(15)-(16)所示，对于任何一组大小为l的数据，都能找到一个整数a满足以下不等式的最小整数。对于本发明
如果对a＝a1的速度场所有的u、v数据集进行规范化处理，a＝a2的压力场p规范化处理。为了保持连续性方程标准化的一致性，则对空间坐标x、y数据归一化的参数a3＝a1，对于黏度μ标准化的a4＝a2。
[0097][0098][0099]
其中，αi表示归一化前的数据，l表示数据的总量，a表示归一化参数，表示归一化后的数据。
[0100]
上述操作主要是把u,v,p的值归一化到[-1,1]，同时保证公式(8)(9)(10)数量级的一致性，对u和v除以相同的数，保证(8)分子分母数量级的一致性a＝a1，对u和p除以相同的数保证(9)减号前后两个式子数量级的一致性a4＝a2。
[0101]
本实施例中，归一化处理的过程如下：对a1＝-6的速度场所有的u、v数据集进行规范化处理，a2＝-1的压力场p规范化处理。为了保持连续性方程标准化的一致性，则对空间坐标x、y数据归一化的参数a3＝a1＝-6，对于黏度μ标准化的a4＝a2＝-1。
[0102]
步骤6、采用归一化后的数据对点云残差网络模型进行训练，设置学习率、参数epoch、batch，训练结束后，采用训练完成的模型对二维多孔介质孔隙流场的速度水平分量u、速度竖直分量v和压力p进行预测，当预测结果表明当前模型性能良好时，输出该模型，将该模型用于真实场景下流场和压力场的预测。
[0103]
本发明实施例中，训练时，学习率设置为0.003，epoch即训练轮数为60000，batch为1。采用测试集进行模型测试时，预测结果如图4-9所示。图4-6分别为u、v、p的预测值与真实值对比图，图7-9分别为距出口50nm的数据点u、v、p的预测值与真实值对比图。从图可以看出，与压力场p相比，预测速度场u、v经历了更高水平的误差，这是因为压力场在x方向上简单地线性减小-而速度模式更复杂。此外，速度矢量涉及质量和动量平衡，而压力仅在动量方程(2)中起作用。整体上u、v、p的相对误差分别为13.2％、14.2％、0.012％，预测精度较良好。
[0104]
当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。