【论文】罗山贵,等：基于混合优化算法和深度神经网络模型结合的致密砂岩气藏裂缝参数优化

致密砂岩气作为一种重要的非常规天然气资源， 近年来产量不断增长，已成为我国天然气增储上产的重要增长极^[1-2]。水平井分段压裂是致密气藏的主要开发方式^[3]，合理的裂缝参数设计是气藏有效开发的关键^[4]。然而，储层低孔隙度、低渗透和强非均质性增加了裂缝参数选择的难度^[5]，传统的裂缝参数设计依赖于现场工程师经验和正交实验等手段，难以得到这些条件下的最佳裂缝参数。因此，石油行业逐步发展了基于油藏数值模拟和智能优化算法的自动生产优化方法^[6-8]，将裂缝参数的设计转化为一个高维、非线性、多约束的数学最优化问题，以提升裂缝参数优化设计的效率和效益。求解此类问题的优化算法可分为局部搜索算法和全局搜索算法^[9-10]。全局搜索算法在不需要梯度信息下即可搜索全局最优解，因而在非常规储层裂缝参数优化研究中备受青睐^[11-13]。但是，以遗传算法（Genetic Algorithm, GA）等群智能优化算法为代表的全局搜索方法，存在收敛速度慢、计算资源需求高的问题。尤其是在非常规油气藏裂缝参数优化中，目标函数的评估涉及数值模拟运行， 计算成本更高^[14]。多算法混合优化是提高优化算法寻优速度和求解精度的常用策略。大多数混合优化策略基于不同群智能搜索算法的混合^[15]，而全局搜索和局部搜索的混合主要是在全局搜索结束后，再使用局部搜索方法精确寻优^[16]。

除了提升优化算法的效率，加快目标函数的求解也是提高优化问题求解效率的常用方法^[17]。随着人工智能的发展，基于机器学习（Machine Learning, ML） 建立代理模型以能解决石油行业中面临的各种昂贵、耗时的计算问题而得到广泛的应用和发展^[18-20]。相比于传统机器学习方法，深度神经网络（Deep Neural Network, DNN）具有特征提取能力强、能够处理大规模数据和适应不同类型数据等优点，是使用频率最高的机器学习模型^[21]。在油气藏生产优化领域，机器学习代理模型可用于减少优化过程中数值模拟运行次数。机器学习代理辅助优化分为非自适应学习和自适应学习。自适应学习仅选择少量样本点建立初始代理模型，随后根据优化过程不断采集新的样本点迭代更新代理模型，所需油藏数值模拟次数更少，同时具有更高的计算效率和预测精度^[22-23]。新采样点的获取即自适应学习策略是影响代理模型精度的关键。常用自适应策略以最有希望解（当前优化阶段获得的最佳解） 和最不确定解（整个可行域空间中代理模型预测精度最低的解）作为新采样点。在选择最不确定解时，现有研究基于相同的样本数据主要采取两种方法建立代理模型：①建立多个目标函数的代理模型，然后以多个模型预测值的方差最大的点作为最不确定解；②同时建立目标函数代理模型和不确定性代理模型^[24-25]。然而，这些过程均涉及多个代理模型的构建和参数调整过程，增加了优化问题的复杂程度。

针对上述问题，笔者建立了基于混合优化算法和自适应DNN 代理模型加速的致密砂岩气藏裂缝参数优化方法。为充分结合全局搜索算法的广泛探索能力与局部搜索算法的快速寻优能力，混合优化算法采用GA 和贝叶斯自适应直接搜索（Bayesian Adaptive Direct Search, BADS）之间循环迭代的混合优化策略。为减少数值模拟运行次数和提高代理模型精度，在自适应学习过程中提出了以最大平均距离点作为最不确定解的自适应学习策略。最后，以实际致密砂岩气藏裂缝参数优化问题测试了所提出的优化方法的先进性。

1 　混合优化算法建立

1.1 　遗传算法

GA 是一种基于自然界生物进化机制的全局元启发式优化方法^[26]，其能够在不需要梯度信息的情况下处理各种类型变量、非线性目标函数以及约束条件。在GA 中，通过编码组成初始种群后，遗传操作的任务就是对种群的个体按照它们对环境的适应度施加一定的操作，从而实现优胜劣汰的进化过程。GA 主要包括染色体编码、种群初始化、适应度函数设计、选择操作、交叉操作和变异操作，而选择、交叉和变异是遗传算法具备解决复杂优化问题和各种优点的关键。

为确保最优解不会因随机的选择、交叉和变异操作而丢失，GA 常采用精英策略。Eiben 等^[27] 用马尔可夫链理论证明了采用精英策略的GA 可以弥补标准选择算子导致的统计误差，保证得到全局最优解。精英策略意味着每次迭代时，都会按照一定的比例从当前种群中选择若干个最优的个体直接保留到下一代种群中。由于最优个体被保留，GA 搜索方向更可能指向全局最优解，从而提高收敛速度。

1.2 　贝叶斯自适应直接搜索

BADS 算法是一种将网格自适应直接搜索（Mesh Adaptive Direct Search, MADS）和贝叶斯优化（Bayesian Optimization）相结合的局部优化方法^[28]。该算法在求解缺乏梯度信息以及非解析或有噪声的目标函数上，表现出较好的性能。BADS 算法通过局部的贝叶斯优化过程（即搜索过程）和系统的、较慢的网格搜索（即轮询过程）之间的相互交替，更为有效地进行空间探索，得到更为精确的代理模型。在BADS 算法的搜索阶段，高斯过程（Gaussian Process, GP）被拟合到当前迭代步数下评估点的子集， 根据低置信限策略，即在不确定性较高的搜索区域（高GP 不确定性）与概率较大的区域（低GP 平均值） 之间权衡，迭代选择新的评估点。与传统的贝叶斯优化不同的是，当搜索过程不断失败时，会进入轮询阶段。该阶段每一次在一个方向前进步数，对该网格上的点进行评估，直到发现改进或尝试了所有方向。如果成功，步长加倍，否则减半。其间，BADS 收集目标函数局部形状的信息，以便为下一次搜索建立更好的代理。

1.3 　混合优化算法

在混合优化算法建立前，首先以多峰非线性Rastrigin 函数^[29] 测试GA 和BADS 算法性能。Rastrigin 函数在全局最小值周围存在多个规则分布的局部最小值，且较多局部最小值与全局最小值接近。对于二维Rastrigin 函数， 在x=(0，0) 处有唯一最小值f(x)=0。如图1 所示，当初始点为（－10，10）时， BADS 算法经过86 次目标函数评价后陷入局部最优解。在更好的初始点（－0.1，0.1）时，BADS 算法经过43 次目标函数评价后获得全局最优解。可见， BADS 等局部搜索算法寻优精度受初始点影响较大， 从优良的初始点出发，BADS 算法更易跳出局部最优解。GA 等群智能优化算法不受初始点的影响，但在本算例中，GA 需经过1 463 次目标函数评价才能收敛到全局最小值。

由上述可知，BADS 算法具有收敛速度快的优点，但优化结果受初始点影响大且易陷入局部最优解。GA 具备较好的全局搜索能力，然而收敛速度较慢、所需目标函数评价次数多。结合两种算法各自的优点，笔者建立了GA 和BADS 循环交替的混合优化算法，其优化步骤如下：

步骤1 ：初始化GA 和BADS 算法参数。

步骤2 ：基于GA 执行k 代的全局搜索。

步骤3 ：若满足截止条件，则结束优化；否则， 以GA 的优化结果作为初始点执行BADS 局部搜索。

步骤4 ：若满足截止条件，则结束优化；否则， 以BADS 算法的优化结果替换种群中最差的个体后转到步骤2。

在上述算法步骤中，连续k 次GA 迭代和一轮BADS 局部搜索的过程称为混合优化算法的一个迭代轮次。其中k 的选择要保证初次GA 优化时目标函数有较大的改进，但是优化结果又未达到全局最优。在该算法中，全局搜索算法为局部搜索算法提供好的初始点，局部搜索算法为全局搜索提供更优良的个体。如图2 所示，针对Rastrigin 函数，GA+BADS 混合优化经过380 次目标函数评价收敛到全局最优解，相比于GA 优化效率提高3.85 倍。

对于裂缝参数优化等实际工程优化问题，目标函数评价依赖于数值模拟等耗时较长的手段。因此，在相同的计算资源下，GA 的时间复杂度最高， BADS 算法的时间复杂度最低，而GA+BADS 混合算法的时间复杂度介于两者之间。

2 　自适应深度神经网络代理加速

2.1 　深度神经网络模型

DNN 的多隐藏层结构极大提高了其提取数据特征的能力，能够更高效地拟合输入变量和输出变量之间复杂的非线性映射关系^[30]。全连接神经网络是最常用的DNN 架构，其具有典型的分层结构，包括一个输入层、多个隐藏层和一个输出层（图3）。每个隐藏层由多个神经元构成，每个神经元与其相邻层中的所有神经元相互连接，同一层中的神经元之间没有连接。隐藏层的数量决定了网络的深度，每层神经元的数量决定了网络的宽度。

神经元是DNN 最基本的计算单元。神经元的计算过程包含线性叠加与非线性激活两部分（如图3 的红色框线内所示）。线性叠加部分是全部输入与各自权重的乘积之和，并与偏置相加。非线性激活函数的引入是神经网络具备非线性映射关系学习能力的关键。以线性整流函数（ReLU）及其变体（LeakyReLU、ELU 等）为代表的非饱和激活函数减轻了模型训练过程中可能出现的梯度消失问题，是最常用的激活函数^[31-32]。

在DNN 训练时，笔者采用Adam（Adaptive Moment Estimation）算法更新网络的所有权重和偏置，以最小化预定义的损失函数（Loss Function）。Adam 是一种结合了动量优化和RMSProp（Root Mean Square Propagation）算法的自适应学习率随机优化算法， 兼具动量优化跟踪过去梯度的指数衰减平均值和RMSProp 跟踪过去平方梯度的指数衰减平均值的特点。Adam 算法具有容易实现、计算效率高和收敛速度快等优点，在大规模深度学习领域备受青睐，其参数更新规则如式（1）～（5）所示。为加快训练速度和提高模型泛化能力，梯度类优化器常与小批量学习（Mini-Batch Learning）结合。小批量梯度下降类算法将总训练集随机划分为多个互斥的小批量样本集，每次基于单个小批量样本执行梯度下降过程，直到遍历整个样本集，即完成一个轮次的训练。经过多个轮次的遍历，最终达到收敛条件。

对于回归任务，常用的损失函数为均方误差（Mean Squared Error, MSE）。为了避免模型的过拟合问题，可采用L2 正则化（又称为岭回归）技术，通过在损失函数中添加一个等于权重向量的L2 范数的平方之半的正则项，迫使学习算法不仅拟合数据，同时使模型权重尽可能小。当添加L2 正则化项后，损失函数如下所示：

批量归一化（Batch Normalization, BN）也是DNN 训练常用的优化技术。批量归一化对神经网络所有隐藏层的每个输入零中心并归一化，然后每层使用两个新的参数向量缩放和偏移其结果，以解决在训练过程中，中间层数据分布发生改变的情况。为了使输入零中心并归一化，该算法需要评估当前小批次上的输入的均值和标准差，其计算步骤如式（8）～（11）所示。在批量归一化操作中，缩放参数（γ）和偏移参数（β） 是神经网络训练时的待学习参数。在模型训练阶段， 基于小批量样本数据计算各层的输入均值向量和输入标准差向量。在模型预测阶段，基于整个训练数据集计算输入均值向量和输入标准差向量。批量归一化可以降低模型对权重初始化和激活函数的敏感性，同时可以使用更大的学习率，从而加快学习过程。

2.2 　自适应学习策略

在机器学习代理辅助的裂缝参数优化中，常以最有希望的解和最不确定解作为新的采样点来更新代理模型。但现有自适应学习策略中对于最不确定解的选择过程过于复杂。因此，笔者提出了以最大平均距离点作为最不确定解的自适应学习策略。如图4 所示，基于拉丁超立方采样^[33] 获得多个新拟采样点， 求得这些拟采样点到所有数值模拟评价点的平均距离。其中，平均距离最大的采样点即为最大平均距离点。

在代理模型更新前，数值模拟评价的真实样本点已经全部用于构建代理模型，能很好地预测这些真实样本点及其周围点的目标函数值。但是在离真实样本点较远的区域，模型在训练阶段所能学到的信息有限，因而这些区域预测精度较低，不确定性较大。而最大平均距离点代表着在采样空间（优化问题的可行域空间）中离所有真实样本点最远的区域，因此可作为代理模型预测不确定性最大的点。

由图4 可知，最不确定解为新采样点中平均距离最大所对应的样本点，最大平均距离计算公式如下：

进一步，笔者建立了基于混合优化算法和自适应DNN 加速的裂缝参数优化流程，其具体执行步骤如下。

步骤1 ：首先确定优化问题的目标函数和优化算法的初始化参数，并基于GA+BADS 混合优化算法和油藏数值模拟执行一定数目的优化过程。同时根据优化问题的规模和复杂程度，可采用拉丁超立方采样等方法获得额外的样本数据点并通过数值模拟评价目标函数值。将此过程的所有数据作为初始样本集。

步骤2 ：基于初始样本集，建立初始DNN 代理模型。

步骤3 ：基于DNN 代理模型和GA+BADS 混合优化算法继续执行优化过程。

步骤4 ：DNN 代理模型优化结果精度计算。以步骤3 所求近似最优解在代理模型和数值模拟上预测值的相对误差作为代理模型的精度。若该相对误差小于阈值，则满足精度条件执行步骤6 ；若不满足则执行步骤5。

步骤5：采用最大平均距离方法获得最不确定解， 基于数值模拟评价最有希望解和最不确定性解的真实目标函数值，同时也可利用拉丁超立方采样获得少量样本点并评价目标函数作为补充。这些少量的补充采样点与最有希望解和最不确定解共同用于更新DNN 代理模型（即采用与初始DNN 模型相同的样本划分比例和超参数再训练一次DNN 模型），随后转到步骤3。

步骤6 ：以代理模型优化的近似最优解作为初始点，采用BADS 局部搜索算法和油藏数值模拟继续寻优，直到收敛。

与以往的自适应代理辅助油藏生产优化方法相比，笔者提出的流程主要包含4 点不同：①优化流程的初始阶段并不是拉丁超立方采样，而是先执行数代的基于数值模拟的优化过程，而拉丁超立方采样仅根据实际问题的复杂程度用于辅助采样；②最不确定解的求取避免了常用的多代理模型构建的复杂过程，仅依据多个新采样点的最大平均距离获取； ③在代理模型更新时，不仅使用了最有希望解和最不确定解，同时也可采用拉丁超立方采样获得少量样本点共同用于更新代理模型；④虽然代理模型的预测精度随优化过程逐渐提高，但是在最优解附近可能难以实现准确的目标函数评价。因此，在代理模型优化结束后，基于局部搜索算法和油藏数值模拟继续执行优化，以提高最终寻优精度。

3 　优化方法验证与应用

3.1 　单段裂缝参数优化

为验证笔者建立的混合优化算法在实际工程问题中的性能，基于四川盆地天府气田中侏罗统沙溪庙组致密砂岩气藏^[34-35] 和Eclipse 油藏模拟器建立了2D 非均质属性模型（图5）。数值模拟模型孔隙度在0.03 ～ 0.12，平均0.08 ；渗透率在0.004 ～ 0.310 mD， 平均0.050 mD。模型大小为300 m×600 m×26 m， 平面网格尺寸为5 m×5 m，总网格数目为7 200。假设此模型包含3 条等间距分布的水力裂缝。水力裂缝垂直于水平井轨迹，因此可采用局部网格加密（Local Grid Refinement, LGR）表征水力裂缝。模型的详细参数如表1 所示。模拟所用气水相渗曲线和应力敏感曲线如图6 所示。

针对上述数值模拟模型，在总水力裂缝长度和水力裂缝导流能力固定的条件下，优化变量为3 条水力裂缝的长度，优化目标为15 a 的累计产气量（Cumulative Gas Production, CGP）。在水力裂缝长度优化过程中，优化变量还需要满足一定边界约束条件。考虑到优化算法通常是寻求最佳自变量以使目标函数最小，因此，本节的优化问题可描述如下：

首先遍历所有优化变量的组合，得到本节优化问题的全局最优解。对于上述优化问题，3 条水力裂缝的长度共有2 792 种组合方案。通过数值模拟得到所有2 792 组方案的15 a 累计产气量。其中， 最小累计产气量为0.239 4×10⁸ m³，最大累计产气量为0.252 9×10⁸ m³。相比于等水力裂缝长度的设计（累计产气量为0.243 5×10⁸ m³），最优水力裂缝长度方案下单个压裂段产量可提高3.86%。

本节优化问题较简单，在GA+BADS 混合优化算法求解时，GA 执行一次迭代就转到BADS 局部搜索求解阶段。其中，GA 的初始化参数如表2 所示， BADS 初始网格步长为1，最小网格步长为0.01。如图7 所示，仅使用GA 优化时，需要156 次目标函数评价才能找到全局最优解。与GA 相比，GA+BADS 混合优化算法所需目标函数评价次数降低18%（128 次）。可以发现，GA+BADS 混合优化算法与GA 的寻优精度相同，但是寻优速度明显提高。如图8 所示， 对于本节的单压裂段裂缝参数优化，最优设计倾向于两边布长缝、中间布短缝的方案。

3.2 　单井裂缝参数优化

本节以全井裂缝参数优化为例，进一步对比优化算法和优化策略的性能。数值模拟模型大小为 1 000 m×600 m×26 m，基质网格大小5 m×5 m，总网格数24 000 个。模型孔隙度在0.030 ～ 0.124，平均0.082；渗透率在0.004 ～ 0.368 mD，平均0.06 mD （图9）。假设水平井位于模型中间，水平段长800 m， 等缝间距分布16 条水力裂缝。气水相对渗透率曲线和应力敏感曲线等与上一节一致。所不同的是，本节以净现值（Net Present Value, NPV）为优化目标，同时优化各条水力裂缝的缝长和裂缝导流能力。因此， 本节的优化问题如下所示：

当水平井位置固定时，钻井成本和地面建设成本基本固定。参考天府气田致密砂岩气藏，钻井成本取3 000 万元，地面建设成本为1 500 万元。而总压裂成本包括压裂液和支撑剂的材料费以及压裂服务费，总压裂成本计算如下：

假设每个压裂段包含4 条裂缝，则压裂段数N=4。而压裂液和支撑剂总购买费用与水力裂缝性质有关。基于天府气田致密砂岩气藏现场资料拟合得到压裂液购买成本与水力缝长和缝高的近似计算公式如下：

当支撑剂类型确定后，初始支撑裂缝渗透率可近似恒定，因此由支撑裂缝导流能力即可计算支撑裂缝缝宽。假设支撑剂在水力裂缝中均匀铺置，支撑剂购买成本可由下式计算^[37] ：

由于优化变量增加（32 个），GA 种群大小取100，最大迭代次数为20，BADS 算法最小网格步长取0.1，其余初始化参数与3.1 节一致。在混合优化算法中，执行5 次GA 迭代后进入BADS 局部搜索。首先，测试GA+BADS 混合优化算法在高维裂缝参数优化问题上的适应性。如图10 所示，GA 经过 2 005 次目标函数评价后最大NPV 为1 799 万元。而GA+BADS 仅需要655 次目标函数评价就获得了与GA 相同的NPV，最终经过1 076 次目标函数评价收敛于1 930 万元的NPV。对于高维优化问题，相比于仅采用GA 优化，GA+BADS 混合优化算法的收敛速度和寻优精度都显著提升（在约1/2 的GA 总目标函数评价次数下明显提高了NPV）。

随后，进一步测试DNN 代理加速优化的效果。在该方法中，首先执行一轮次的基于混合优化算法和数值模拟求解目标函数的优化过程。随后建立初始DNN 代理模型，在后续优化过程中，则采用代理模型求解目标函数值。其中，以最有希望解上DNN 与数值模拟预测值相对误差作为代理模型精度判断条件。当两者相对误差小于1×10^－3，则结束基于代理模型的优化。在采用最大平均距离方法求取最不确定解时，利用拉丁超立方采样获得了500 个可行解。在代理模型更新时，除了最有希望解和最不确定解， 还利用拉丁超立方采样获得了额外的5 个样本点。因此，每次代理模型更新时，需要7 次基于数值模拟的目标函数评价。这些新的样本点和原有样本共同构成代理模型更新时的样本数据集。

将第一轮次优化过程收集的样本数据按照9∶1 的比例随机划分得到训练集和验证集。DNN 输入参数为所有裂缝的缝长和导流能力，输出为NPV。模型以均方误差为损失函数，并采用表3 所示超参数训练500 轮次，以训练过程验证集上均方误差最小时的网络权重和偏置作为最终模型。在代理模型更新时， 仅样本数据的规模逐渐增大，样本划分比例和超参数均保持不变。如图11 所示，初始DNN 共迭代2 500 次（每轮次迭代5 次），在第1 703 次迭代时验证集上损失和均方根误差最小。此时，训练集均方根误差为30.24，验证集均方根误差为35.93。

DNN 训练仅需要数秒钟，因而仍然以数值模拟目标函数评价次数作为优化效率对比指标。如图12-a 所示，使用GA+BADS 和DNN 代理模型加速的优化流程经过812 次目标函数评价收敛于1 930 万元的NPV。其中，用于样本收集的首轮优化经过566 次目标函数值计算获得1 820.2 万元的NPV。自适应DNN 优化阶段再经过50 次目标函数评价获得了最大1 924 万元的NPV。受训练样本随机划分和总训练样本数目较少的影响，自适应DNN 优化阶段NPV 呈波动变化。因而以此阶段最大NPV 对应的优化解作为下一阶段局部搜索的起点。最终的BADS 局部搜索阶段经过196 次目标函数评价，NPV 由1 924 万元提高至1 930 万元。

不同方法对比可知，GA 花费了最多的数值模拟次数且NPV 最低（图12-b）。与GA 相比，GA+BADS 混合优化和DNN 代理加速优化在更少的目标函数评价次数下NPV 均明显增加。相比于GA+BADS 混合优化算法，DNN 代理加速优化可减少24.54% 的目标函数评价次数。

本研究聚焦于混合优化算法的建立和代理模型的加速过程，采用了截断法和罚函数法处理边界条件和线性等式约束。对于更复杂的约束条件，可能需要过滤器法^[38] 或增广拉格朗日函数法来进一步提高算法的求解效率。此外，对于致密砂岩气等非常规油气藏的开发优化，可以考虑将裂缝扩展作为优化过程的一部分，进而实现对压裂施工参数的直接优化。

4 　结论

1）BADS 局部搜索算法具有收敛速度快的优势， 然而，对于多模态的复杂优化问题，算法性能受初始点影响较大且容易陷入局部最优解；GA 具有良好的全局探索能力，但收敛速度慢，所需目标函数评价次数高；所建立的GA+BADS 混合优化算法表现出了显著优于GA 的收敛速度和寻优精度。

2）针对低维致密气藏裂缝参数优化问题，在收敛于相同最优解时，GA+BADS 混合优化算法可降低18% 的数值模拟运算次数；针对高维裂缝参数优化问题，GA+BADS 混合优化算法在约1/2 的GA 总数值模拟次数下提高了131 万元的NPV。

3）在高维致密砂岩气藏裂缝参数优化问题上， 基于GA+BADS 混合优化算法和自适应DNN 代理加速的优化流程性能最佳。相比于GA+BADS 混合优化算法，在获得相同NPV 时，自适应DNN 代理加速优化可减少24.54% 的数值模拟运行次数。