基于改进BP神经网络的页岩地应力预测模型

2021-08-02尚福华王玮卿曹茂俊

计算机技术与发展 2021年7期

尚福华，王玮卿，曹茂俊

(东北石油大学计算机与信息技术学院，黑龙江大庆 163318)

0 引言

页岩气是中国能源勘探的重要领域和目标。其中，地应力分析技术是指导页岩压裂施工的核心技术[1]。目前，初始地应力数据主要靠现场实测获得，但是通常需要在现场利用水力压裂法，应力-应变恢复法等方法进行测量，这些测量方法测试周期较长，成本昂贵，测试条件有限，所以目前地应力场一般的演算方法为以实测地应力数据为基础，通过数值模拟或数值分析来求解[2]。然而，由于地层在漫长的时期中经历多次地质作用及构造演化，地质构造类型、作用程度和方向不同，导致构造应力场的大小和方向存在时间和空间的变异性，只能通过它的外部相关因素的状态和方式来把握该事物的信息[3]。

神经网络具有自学习、自组织、可以高速寻找最优解等优点，可成功拟合多数非线性连续函数。因此，针对建立工程区域内相关的参数与岩体初始地应力非线性映射关系问题，采用神经网络是一种有效的解决手段[4]。目前，已有一些学者应用BP神经网络分别对三维应力场和岩石力学参数进行了反演分析。侯连浪等[5]建立并训练了BP神经网络对页岩静弹性模量进行预测，杨志强等[6]利用金川矿区实测地应力样本建立了人工神经网络模型进行训练，获得了金川矿区地应力分布规律。然而，由于BP神经网络的优化算法为梯度下降法，目标函数通常又极其复杂，因此，存在着网络训练易过拟合、易陷入局部最优、收敛速度较慢以及隐层节点数难于确定等缺点[7]。

目前，一些传统的自适应梯度下降算法可以通过自适应地为各个参数分配不同学习率来解决该问题，但同样存在随着迭代次数增加，学习率越来越慢而导致提前收敛的问题[8]。因此，针对页岩气地应力预测模型的构建问题，该文以页岩横观各向同性为基础，针对传统的自适应梯度下降算法使用改进的梯度下降算法，并将改进后的算法应用于神经网络中，以长宁工区若干口井的声波测井资料与对应深度的实测地应力值数据为依据，反演出目标层段的应力边界条件，最后代入到数值模拟模型中进行正演计算，进行该工区初始地应力预测。

1 页岩地应力物理模型

1.1 横观各向同性物理模型

近年来，地应力计算模型的建立大多数针对各向同性储层，其中黄荣樽模型[9]的应用比较广泛：

(1)

(2)

由于页岩物质组成不均匀，节理发育不均匀，因此平行结构面方向和垂直结构面方向物理力学性质不同，呈现出横观各向同性的特征。鉴于此特性，计算时仍采用上述各向同性模型显然不够准确。目前，研究中一般采用Thiercelin模型[10]，具体如下：

(3)

(4)

式中，σv为垂直地应力；σH为水平最大主应力；σh为水平最小主应力；Pp为地层孔隙压力；α为Biot系数；ξ为渗流系数；Eh为层理面横向的杨氏模量；EV为层理面法向的杨氏模量；vh为层理面横向的泊松比；vv为层理面法向的泊松比；εH为水平最大构造应变；εh为水平最小构造应变。其中垂向主应力的数值近似等于上覆岩层压力σv，可以看作为地层密度和深度的函数，因此可用如下公式求出垂向地应力：

(5)

式中，h为地层埋藏深度(单位m)；ρ(h)为地层密度随地层深度变化的函数；g为重力加速度，9.8。

实际的地层密度ρ(h)随深度的变化规律复杂，用连续函数不易表示，所以在实际计算中可以用分段加和的方法来计算：

(6)

式中，ρi为第i段地层平均体积密度(单位kg/m3)；ΔDi为第i段地层厚度(单位m)。

上覆岩层压力计算完毕后，地层岩石的弹性参数可由钻井取芯进行室内试验获得，但由于钻井取芯的不连续性及高成本，可以通过声波测井资料、密度测井资料和自然伽马测井资料来计算获得。

1.2 改进的页岩地应力计算模型

由于大多数研究区构造复杂，页岩具有明显的层理结构，受层理分布的影响，地层自上而下倾角与倾向均有明显变化，而地层倾角使得地层的上覆岩层压力降低，增大了水平应力，使水平地应力的计算结果有一定误差。并且，由重力、构造应力共同作用下的岩石应变值在实际工程中不易获得。为了提高地应力计算精度，该文在Thiercelin模型的基础上，考虑地层倾角因素，对上述地应力计算模型进行改进：

(σv-αPp)sinφsin(β-β0)+αPp

(7)

(σv-αPp)sinφsin(β-β0)+αPp

(8)

式中，A为最大水平主应力方向构造运动(位移)系数；B为最小水平主应力方向构造运动(位移)系数；φ为地层倾角(单位°)；β为方位角(单位°)；β0为最大水平主应力方向角(单位°)；K为侧压系数。

模型基于空间三角函数应力分布关系进行计算，其中第一项首先通过构造应力系数反映水平构造剧烈程度，其次反映由于倾角作用使得上覆岩层压力对地层的压实程度的减少情况；第二项考虑井斜方位角与构造运动方向不同而产生的剪切应力；最后，公式采用上覆岩层压力产生的水平构造应力与原始水平构造应力加和的方式，代替利用地层的弹性应变求弹性应力的方式，避免使用因深度过大，取芯困难从而难以获得的应变数据。当地层倾角为0时，对地应力的影响也为0，此时模型还原为黄氏模型。

2 BP神经网络模型

在上述分析模型中，岩层岩石力学参数，包括各方向的杨氏模量、泊松比，井斜方位角等参数虽可用测井曲线资料计算获得，但由于岩石在不同方向所受的构造运动，温度影响不同，使得模拟地层构造环境过程复杂，不同水平方向的构造应力往往不同。因此模型中构造应力系数只能通过该地区的经验系数获得，结果误差较大。为了解决地层构造应力系数获得较难的问题，该文利用BP神经网络的映射功能反演某区块气藏地层的水平构造应力系数，形成最终该地区的地应力解释模型。

2.1 BP神经网络算法

BP神经网络是一种按误差逆传播算法训练的多层前馈网络。学习过程中由信号的正向传播与误差的逆向传播两个过程组成。正向传播时,各神经元从输入层开始，接收前一级输入，经激活函数计算后并输入到下一级，直至输出层。反向传播时,将输出层输出值与期望值的误差逐层返回,并不断修改权值矩阵。权值不断修改的过程,也就是网络学习的过程，此过程一直进行到网络输出的误差逐渐减少到可接受的程度或达到设定的学习次数为止。三层网络结构如图1所示，其中[X1,X2,…,Xn]为输入样本，[Y1,Y2,…,Ym]为输出样本。

图1 三层BP神经网络结构

令隐层权值矩阵为Wij，输出层权值矩阵为Tjk，则：

隐层各单元的输出按下式计算：

Oj=f(∑WijXi-θj)

(9)

输出层的输出结果为：

Yj=f(∑TjkOj-θk)

(10)

式中，θ表示神经元阈值，f表示非线性函数。随后按BP神经网络的梯度下降准则，定义误差函数ek后，进行权值更新，直到误差收敛至期望时训练结束：

(11)

Wjk=Wjk+ηHjek

(12)

式中，i=1,2,…,n；j=1,2,…,l；k=1,2,…,m。

2.2 改进的神经网络算法

虽然BP神经网络具有误差反向传递的特征，但由于学习速率固定，仍存在以下显著不足：

(1)若学习速率过快，有时会出现“锯齿形现象”，导致网络难以收敛；

(2)若学习速率过慢，即使网络成功收敛，但首先收敛速度过慢，其次因为误差平方和函数易陷入局部极小值，可能得到的只是局部最优解[11]。

鉴于常数型学习率的不足，目前为了解决此问题常采用自适应学习率算法。例如，Adagrad算法中网络每代的学习率利用初始学习率除以历代梯度平方和的均方根，根据历史梯度的变化情况调整大小达到实现自适应效果[12]，即：

(13)

其中，Gt为累积历史梯度的平方，即：

Gt=Gt-1+gt2

(14)

由于Adagrad算法的全局学习速率除以的是累加梯度的平方，到后面累加的比较大时，会导致每次的学习率逐渐减小，导致梯度更新缓慢，最终可能会导致训练中后期就停止迭代。为了克服其局限性，该文将Adagrad算法中的学习率通过对过去累积的梯度赋予不同的权重，对学习率改进的数学表达式如下：

(15)

(16)

该改进的好处在于，在充分考虑了为权重设计不同的学习率的基础上，分子上利用幂指数函数保证学习率稳定下降且计算量较小,分母上利用加权系数自适应学习速率对最新的随机梯度估计赋予了更多的权重，从而把握梯度变化的方向而加快收敛，避免震荡，此外还可以避免Adagrad算法关于考虑历史梯度的问题。

3 地应力预测模型构建

该文地应力预测模型构建包括页岩地应力计算模型的构建和神经网络的反分析两部分。神经网络的反分析即利用神经网络对初始地应力进行反演得到边界条件预测值，即模型中的参数。选择目标层段内若干点的实测应力值作为输入，通过网络反分析，得到区域的模型参数，然后建立上述地应力预测模型，利用测井曲线获得的声波波速计算岩石力学参数，最后将岩石力学参数，模型参数代入模型内进行计算，即可得到目标区域内任意一点的应力值。具体步骤如下所示。

3.1 测井资料处理解释

该文的初始数据来源为目标工区裸眼井阵列声波测井资料。首先对实际的波形数据进行了处理：根据仪器特性，选择合适的时窗长度，使得在相关系数等值线图上可有效区分各个模式波对应的相关系数峰值；根据各模式波在全波波形上的持续时间，确定时窗寻峰的开始时间和结束时间以及时窗移动的步长；通过分析目的层段地质、测井等资料，估算各模式波时差的大致范围，从而确定寻峰的最小时差、最大时差以及时差变化的步长；根据上述确定的时窗长度、时窗寻峰的开始时间、时窗寻峰的结束时间等参数计算相应的相关系数；由相关系数最大值对应的时间和时差得到模式波的波至时间以及模式波的时差，最后得出声波波速数据。

声波波速确定后，接下来确定刚度矩阵中的C11,C13,C33,C44和C66五个参数。C33和C44可通过波速资料直接求取，由于缺少斯通利波数据，C66通过MANNIE3模型[13]的横波和纵波各向异性参数之间的相关性进行计算，C11,C13用标准的ANNIE模型[14]进行计算，最终利用刚度矩阵元素根据相应岩石力学参数经验公式[15-16]计算水平、竖直向的杨氏模量、泊松比等其他岩石力学参数。

3.2 确定神经网络训练样本

按上述步骤计算后，模型中还缺少最大、最小水平主应力方向的构造运动系数，此时需建立实测应力值与模型参数的映射，从而获得模型参数，即地应力的反演。首先选择5个地应力实测点，根据提取的声波资料值计算相应的岩石力学参数，其次根据目标区块的地质分析，确定两个模型中待反演的参数A,B的上下限，并取5个均布的水平。最后将岩石力学参数条件按5×5正交与均匀设计原理进行组合，通过正演计算计算出相应条件下的地应力值，从而建立样本。

3.3 神经网络特征选取与结构设计

(17)

式中，h为隐层节点个数，m为输入层节点个数，n为输出层节点个数，α为调节常数。网络训练采用sigmoid函数；由于输出层的神经元个数与实际输出的个数相对应，取纯线性传递函数purelin根据随机抽取的方法将样本数据分为训练集和测试集。训练集用于应力值与岩石力学参数的拟合，测试集用于对网络模型计算均方误差和筛选活动神经元。应用以上样本，初始化权值、阈值矩阵、最大迭代次数，学习率，对网络进行训练，采用均方误差函数作为回归损失函数，即：

(18)

式中，N表示样本个数，y'表示神经网络预测输出样本，y表示真实输出。当误差函数下降至期望时迭代结束。

3.4 神经网络反演结果处理

将最优化构造运动参数施加在已建立的横观各向同性模型建立正分析模型上进行正分析计算，最后与所选测点的实测应力值进行对比。

实验流程如图2所示。

图2 实验流程

4 工程实例应用

针对西南油气田某区块内的龙马溪组页岩进行地应力分布预测。龙马溪组页岩为脆性层理页岩，应力敏感性系数差异较大,各向异性明显，岩体较为致密，岩石力学强度较高，页岩气含量为5.09～7.48 m3/t，并且随深度增大而增大。该地区历经的多期次构造活动对区域页岩气的保存产生了重大影响，是页岩气的有利保存区，因此具有一定的研究价值。

4.1 样本数据预处理

根据目标区域现场地质测试分析，可确定若干个测点的实测应力值与岩石力学参数，取5个均布的测点，根据岩石力学参数并反算出这些测点的模型中参数值(见表1)并组合，将岩石力学参数、模型中的参数按照正交设计表进行设计，共得到25组样本(见表2)，并根据上述公式计算得到结果(见表3)。

表1 初始岩石力学参数

表2 正交设计计划表(部分输出样本)

表3 部分输入样本

4.2 实验及结果对比

该文基于Matlab软件对三层BP神经网络进行建模。将网络的最大迭代次数设置为1 000，误差期望值设置为10-4。隐藏节点个数取输入单元数的(2n+1)倍，该文根据折中的思想，将隐含节点个数分别设置为13、15、17、19、21分别进行仿真实验，经比较后隐层选择17个节点进行网络训练，网络的误差函数收敛情况如图3所示。

图3 误差函数收敛曲线对比

由于200代之后误差未发生显著变化，则图3只保留200代前的收敛情况。图3显示，普通的BP神经网络的误差函数由于学习率固定出现了提前收敛至10-2并不再下降的情况，而改进的自适应学习率算法的网络收敛速度较快，均方误差值稳定下降，当迭代过程进行到80代之后，误差趋于稳定并达到期望值，隐含层最优中心权值已获得。

网络训练结束后，将第21～25组样本作为测试样本代入数值模拟模型进行正演计算，将计算结果作为模拟实测值，再代入网络进行反分析，即可模拟出实测值下的岩石力学参数与模型参数。反演结果如表4所示。

表4 参数反演结果

将上述得到的反演结果代入到模型中，重新进行正演计算，并进行5个测点的预测值与现场实测数据的结果对比与误差计算，结果如表5所示。

表5 不同的BP神经网络反演参数计算结果对比

表5显示，两种训练方法得到的结果是接近的，利用BP神经网络反演结果的最大相对误差为16.58，最小相对误差为0.01；利用自适应学习率网络进行反演的最大相对误差为11.4，最小相对误差为0.1。根据国内外资料显示，初始地应力值测量的允许误差达到25%～30%[17]。虽然2种方法的计算精度都能够较好地满足工程要求，但改进BP神经网络的相对误差较小(除个别测点外)，精度较高，计算收敛代数少，网络性能得到很大程度的提高。所以该改进有一定的实际意义。