0 引言

煤层气作为非常规油气的重要组成部分，具有重要的研究价值^[1-3]。煤层气开发是集多种技术于一体的复杂系统工程，而井位选择是煤层气勘探开发的基础，涉及到地质构造条件、煤层发育特征、煤层气物性、施工工程条件等一系列因素，是国内外研究的热点^[4-5]。国内研究人员多依据油藏地质情况制定若干井网方案，通过油藏数值模拟进行井位优化研究，此方法对于研究人员的经验要求较高，但油藏地质由于不确定性因素过多，此方法不易得到最优解^[6-9]。

国外研究人员主要通过优化算法来筛选最优解，首先选定一个目标函数（净现值或采收率等），然后以算法为手段，调用数值模拟软件优化方案，迭代循环得到最佳方案。目前，广泛应用于井位优化的无梯度优化算法多是一些全局优化算法，比如遗传算法（GA）、粒子群算法（PSO）、模拟退火算法（SA）、蚁群算法（ACO）等。其中，模拟退火算法、粒子群算法等在井位优化中的应用越来越广泛。1995年Beckner等^[10]将模拟退火算法引入到井位优化研究中；2002年Norrena等^[11]采用模拟退火算法并考虑经济约束进行井位优化研究；2010年On‐ wunalu等^[12]将粒子群算法引入到井位优化研究中；2011年Onwunalu等^[13]采用粒子群算法进行规则井网优化研究。国内姜瑞忠等^[14-15]分别于2014年和2018年提出用遗传算法和新型遗传算法进行井位优化研究。

本次煤层气井位优化研究以沁水盆地沁端区块煤层气开发为例，基于粒子群算法的劣势和模拟退火算法的优势，提出一种混合粒子群算法，其以净现值为目标函数，单井控制面积和井位为变量，采用优化理论结合油藏数值模拟的方法，优选出净现值最大的单井控制面积和井位，并将优化结果与常规矩形布井方式进行对比，以期为煤层气开发井位优化提供新方法。

1 混合粒子群算法

粒子群算法是由Kennedy等^[16]于1995年提出的一种智能优化算法。它类似于遗传算法，也是从随机解出发，通过适应度评价解的优劣，迭代寻找最优解。其优点是概念比较清晰，过程简单，涉及的计算参数比较少，所需粒子数不多，比较容易实现，但是它也有缺陷，比如容易早熟，收敛速度慢，容易陷入局部最优等。对粒子群算法的改进大致归为2类：一类是在基本的算法基础上通过增加某些策略改进其基本算法的性能；另一类则是将粒子群算法与其他的优化算法混合而形成一种混合算法。混合运用多种智能算法可以取长补短，加快收敛速度的同时有效地避免早熟问题。

模拟退火算法的思想最早是由Metropolis等^[17]于1953年提出来的，1983年Kirkpatrick等^[18]将其应用于组合优化。模拟退火算法是基于蒙特卡罗迭代求解策略的一种随机寻优算法，其出发点是基于固体物质的退火过程与组合优化问题之间的相似性，它具有跳出局部最优和搜索精度高的优点。

基于此，笔者提出一种混合粒子群算法，该算法结合粒子群算法和模拟退火算法，将模拟退火思想引入粒子群算法中，增强了粒子群算法跳出局部最优的能力^[19-21]，并通过matlab编程来实现。混合粒子群算法按照模拟退火接受准则允许目标函数依概率在有限范围内变坏^[22]。

本文将混合粒子群算法运用到井位优化过程的主要思路为：

（1）粒子群初始化。设定粒子数N_s、维数D、惯性权重ω、最大迭代次数k_max、学习因子c₁/c₂，随机产生N_s个初始井位X_i⁰和初始速度V_i⁰，迭代次数k = 0。

（2）模拟退火初始化。设定初始温度T_i⁰、允许变坏的范围Fval和退温系数α。

While（迭代次数k＜k_max）do

（3）调用Eclipse数值模拟器计算试验区的累计产气量与累计产水量，进而计算每个粒子的适应度fit_i^k。

（4）根据粒子适应度fit_i^k，选出种群的个体极值P_{i, d}^k，全局极值P_{g, d}^k。

（5）按式（1）更新速度V_{i, d}^k+1，并把它限制在最大速度| V_max |内

$ \begin{array}{l} V_{i, d}^{k + 1} = wV_{i, d}^k + {c_1}r_1^i\left( {P_{i, d}^k - X_{i, d}^k} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;{c_2}r_2^i\left( {P_{g, d}^k - X_{i, d}^k} \right) \end{array} $

(1)

式中：r₁^k，r₁^k均为0~1的随机数。

（6）按式（2）更新井位V_{i, d}^k+1，并把它限制在最大井位| X_max |内

$ X_{i, d}^{k + 1} = X_{i, d}^k + X_{i, d}^{k + 1} $

(2)

（7）计算井位变化引起的适应度的变化量Δf=fit_i^k+1-fit_i^k，若Δf＜Fval，则接受新井位X_{i, d}^k+1及新速度V_{i, d}^k+1，转到步骤（8）；否则，计算接受概率p_i^k[式（3）]。若p_i^k＞rand_i^k(0, 1)，接受新井位X_{i, d}^k+1及新速度V_{i, d}^k+1，转到步骤（8）；若p_i^k＜rand_i^k(0, 1)，保持原井位及速度，转到步骤（8）

$ p_i^k = \exp \left( {\frac{{\Delta f}}{{T_i^k}}} \right) $

(3)

（8）k=k+1，更新温度T_i^k。

End

（9）输出最优井位X和最优净现值P_{g, d}。

2 地质模型建立 2.1 区块简介

沁端区块位于沁水盆地南部，隶属于山西省沁水县，区内由老到新依次发育奥陶系、石炭系、二叠系、三叠系和第四系。本区构造形态总体为一走向北北东、倾向北西西的单斜构造，在此基础上发育了一系列近南北—北北东向的宽缓褶曲，形成区内地层的波状起伏，岩层倾角一般不超过15°，个别地段受构造影响地层倾角变化大，断层不发育。

区块内主要含煤地层为上石炭统太原组（C₃ t）和下二叠统山西组（P₁ s）。含煤17层，煤层总厚度为14.67 m，其中主要可采煤层2层，分别为3#煤层和15#煤层，3#煤层厚度为6.05~6.60 m，平均为6.24 m；15#煤层厚度为3.8~4.5 m，平均为4.14 m。开发方式多采用直井压裂合采，单井稳定日产气为1 000~2 000 m³。

2.2 模型建立

基于Petrel地质建模软件，利用生产井测井资料，建立区块的地质模型（图 1），纵向上分为3个小层（包括隔层），x，y方向上的网格步长为20 m，模型的总结点数为241 029（图 2）。根据测井解释的成果，使用克里金插值方法建立区块的煤层气含量、孔隙度、渗透率和净毛比等属性模型。

本区块煤层气藏整体面积大，若选取整个气藏进行井位优化的数模工作，计算量大，时间长。因此，依据各层的物性分布和储量丰度特征，选取1.2 km×1.2 km的试验区进行井位优化数值模拟工作。模型网格划分为61×61×3，x，y方向上的网格步长均为20 m，所选模型原始地质储量为3.52亿m³。生产制度为直井压裂，定压生产（0.5 MPa），总共生产5 000 d。

3 基于混合粒子群算法的井位优化 3.1 井位优化目标函数（净现值）

煤层气开发的井网部署主要考虑3个方面的因素：地质因素、开发因素与经济效益。对于渗透率、煤层气含量、孔隙度等地质因素以及井间干扰、开发层系、采气速度等开发因素，可以通过油藏数值模拟进行研究，而经济效益决定该煤层气区块能否进行商业开发或者采取增产措施，这就要求合理的井位不仅须要满足采气速度的要求，更要在经济上可行。因此，本次井位优化研究选用净现值为目标函数^[23]。

净现值是指在一个项目的整个建设和生产服务年限内各时间段的净现金流量按照设定的折现率折成现值后求和所得到的值。净现值越大，项目收益越大。其表达式为

$ NPV = \int_{{t_0}}^{{t_t}} {\left[ {{C_{{\rm{in }}}}(t) - {C_{{\rm{out }}}}(t)} \right]} {(1 + r)^{ - t}}{\rm{d}}t $

(4)

其中

$ {C_{{\rm{in }}}}(t) = {Q_{\rm{g}}}(t)Gas\;Price - {Q_{\rm{w}}}(t)Water\;OPEX $

(5)

$ {C_{{\rm{out }}}}(t) = OPEX + CAPEX + TAX $

(6)

式中：NPV为净现值，元；t₀为项目的开始年份；t_f为项目的结束年份；C_in（t）为t时刻的现金流入量，主要是销售收入，元；C_out（t）为t时刻的现金流出量，主要是投资、成本及销售税金，元；r为折现率，指将未来有限期预期收益折算成现值的比率；Q_g（t）为第t年的累计产气量，m³；Q_w（t）为第t年的累计产水量，m³；Gas Price为煤层气价格，元/m³；Water OPEX为污水处理价格，元/m³；OPEX为操作成本，主要包括人员工资、水电费、运输费用、维护费、设备更新费等，元；TAX为税费，元；CAPEX为基建费用，主要包括钻井、完井、固井、射孔、压裂、酸化等措施费用以及所需物资的费用，元。

油田开发中采用净现值方法计算利润时，须要将开发时间划分为若干个时间间隔，因此，净现值的目标函数变为离散型

$ \begin{array}{l} NPV = \sum\limits_1^t {\left[ {{Q_{\rm{g}}}(t){{ Gas \;Price }} - {Q_{{\rm{w }}}}(t)} \right.} {{ Water }}\;OPEX - \\ \;\;\;\;\;\;\;\;\;\;\;\;OPEX - CAPEX - TAX]{(1 + r)^{ - t}} \end{array} $

(7)

计算净现值用的参数如表 1所列。

3.2 变量及约束条件

井位坐标X（I，J）为井位优化研究的变量，而且在地质模型中，井位位于网格内。假定井数为N_w，则问题的维数为D= 2 N_w。井位坐标X如下：

$ X = \left[ {\begin{array}{*{20}{c}} {{I_{1, 1}}}&{{J_{1, 2}}}&{{I_{1, 3}}}&{{J_{1, 4}}}& \cdots &{{I_{1, d}}}&{{J_{1, d + 1}}}& \cdots &{{I_{1, D - 3}}}&{{J_{1, D - 2}}}&{{I_{1, D - 1}}}&{{J_{1, D}}}\\ {{I_{2, 1}}}&{{J_{2, 2}}}&{{I_{2, 3}}}&{{J_{2, 4}}}& \cdots &{{I_{2, d}}}&{{J_{2, d + 1}}}& \cdots &{{I_{2, D - 3}}}&{{J_{2, D - 2}}}&{{I_{2, D - 1}}}&{{J_{2, D}}}\\ {}&{}&{}&{}& \vdots &{}&{}&{}&{}&{}&{}&{}\\ {{I_{i, 1}}}&{{J_{i, 2}}}&{{I_{i, 3}}}&{{J_{i, 4}}}& \cdots &{{I_{i, d}}}&{{J_{i, d + 1}}}& \cdots &{{I_{i, D - 3}}}&{{J_{i, D - 2}}}&{{I_{i, D - 1}}}&{{J_{i, D}}}\\ {}&{}&{}&{}& \vdots &{}&{}&{}&{}&{}&{}&{}\\ {{I_{Ns - 1, 1}}}&{{J_{Ns - 1, 2}}}&{{I_{Ns - 1, 3}}}&{{J_{Ns - 1, 4}}}& \cdots &{{I_{Ns - 1, d}}}&{{J_{Ns - 1, d + 1}}}& \cdots &{{I_{Ns - 1, D - 3}}}&{{J_{Ns - 1, D - 2}}}&{{I_{Ns - 1, D - 1}}}&{{J_{Ns - 1, D}}}\\ {{I_{Ns, 1}}}&{{J_{Ns, 2}}}&{{I_{Ns, 3}}}&{{J_{Ns, 4}}}& \cdots &{{I_{Ns, d}}}&{{J_{Ns, d + 1}}}& \cdots &{{I_{Ns, D - 3}}}&{{J_{Ns, D - 2}}}&{{I_{Ns, D - 1}}}&{{J_{Ns, D}}} \end{array}} \right] $

(8)

由于网格在x和y方向是有限个的，因此，约束条件为：

（1）1＜I_{i, j}＜61，I_{i, j}为第（j + 1）/2口井在x方向的网格位置，整数。

（2）1＜J_{i, j}＜61，J_{i, j}为第j/2口井在y方向的网格位置，整数。

3.3 井位优化流程

基于建立的地质模型及目标函数，将混合粒子群应用到井位优化中。基于混合粒子群算法的井位优化流程如图 3所示。

井位优化过程中混合粒子群算法中惯性权重ω、最大迭代次数k_max、学习因子c₁/c₂、粒子速度V、允许变坏的范围Fval和退温系数α等基本参数取值如表 2所列。

4 井位优化方法应用 4.1 可行性研究

在本次可行性研究中，通过研究一口生产井的最佳井位来证明本方法的优势。针对本模型，共有61×61= 3 721个位置，不同位置的含气量与渗透率不同。因此，通过3 721次穷举计算可以获得最佳井位W1（50，52）（参见图 2）与全区的净现值分布（图 4）。

对比图 2（a）与图 2（b）可以发现，与含气量分布相比，渗透率分布对煤层气开发的效果影响更大。分析其原因为渗透率决定井的采气速度，而含气量决定井的生产时间，只有当采气速度和生产时间达到一个最优组合时才能使单井效益最大化。

分析图 4可知，净现值面有若干个极大值点，这是由于储层的非均质性引起的，最优井位为W1（50，52）。图 5为利用混合粒子群算法进行井位优化研究的净现值图，分析发现混合粒子群算法可以快速收敛（Ns k= 270次）到最优井位W1（50，52），获得最大净现值。

可行性研究结果表明，基于混合粒子群算法的井位优化算法较穷举法计算量大幅度降低，可以准确确定最优井位，用于实际井位优化。

4.2 单井控制面积对开发效果的影响

单井控制面积是一口井所控制的开发面积的大小^[24]。煤层气主要采用矩形、菱形井网等，单井控制面积为0.09~0.64 km²。单井控制面积涉及气田开发指标和经济效益的评价，是煤层气气田开发的重要参数，其大小与井型和井间距大小有关，它取决于储层的性质以及煤层气开发的规模。

针对试验区，利用混合粒子群井位优化算法，用数值模拟技术研究不同的单井控制面积（0.06，0.08，0.10，0.20，0.30，0.40，0.50，0.80，1.44 km²）对开发效果的影响。图 6为不同单井控制面积下的最优井位分布图，分析发现不同单井控制面积下生产井数不同，最优井位会发生变化，且位于渗透率高的位置，无明显井网形式。

图 7为不同单井控制面积下的生产指标曲线，单井控制面积越大，区块的净现值越小，不存在明显的拐点；单井控制面积对气井单井净现值影响较大，存在明显的拐点。当单井控制面积大于0.2 km²时，单井净现值随其增大而减小，原因是煤层气井通过排水降压达到解吸开采，单井控制面积越大，导致井间干扰越小，压力降低幅度越小，不利于煤层气降压开采。因此，为获得最大的区块开发收益及单井效益最大化，选择最佳单井控制面积为0.2 km²。

4.3 最优单井控制面积下的井位优选

采用单井控制面积为0.2 km²进行沁水盆地煤层气田井位优化设计，当种群数量设置为30、迭代次数设置为100时，混合粒子群算法得到的净现值为34 115 452元，大于常规矩形井网的净现值30 312 502元，净现值增加幅度为12.55%。

从图 8还可以看出混合粒子群算法寻优能力较强，能够快速定位最佳井位。

图 9为混合粒子群算法和常规矩形布井的井位分布图。从图 9可以看出，最优井位与含气量及渗透率分布密切相关，与含气量分布相比，渗透率分布对煤层气井位选取的影响较大。最优井位多分布在渗透率高的部位，而渗透率低的部位少；含气量大的部位布井较少，原因是此处渗透率低，无法满足一定的采气速度。最优井位是含气量与渗透率的最优组合，此处生产井具有丰富的物质基础与合理的采气速度。

对比图 9（a）与图 9（b）还可以发现，混合粒子群算法优选出的井位与常规矩形井网类似，分布比较均匀，这与本区块非均质程度较小有关。

5 结论

（1）井位优化方法可行性研究验证了基于混合粒子群算法的井位优化方法较穷举法有较大优势，该方法能够快速收敛到最优井位，减少了计算量。

（2）基于混合粒子群算法的井位优化方法克服了常规井网部署的经验依赖性，能够确定储层非均质条件下的最优井位，最大限度地提高了煤层气田的经济效益。

（3）选取不同的单井控制面积（0.06~1.44 km²）进行沁水盆地煤层气井位优化研究，得到最优单井控制面积为0.2 km²，对于最优单井控制面积，混合粒子群算法得到的净现值比常规矩形井网净现值增加12.55%。

（4）最优井位分布与含气量、渗透率均密切相关，其中渗透率影响尤为重要，最优井位是含气量与渗透率的最优组合，既要有丰富的物质基础，又要有合理的采气速度。