中国油田大多属于陆相沉积, 其非均质性普遍较强, 在长期的注水开发过程中容易形成水窜现象, 严重影响了开发效果, 统计结果表明水驱后仍有约2/3的原油滞留于地下, 提高采收率潜力巨大^[1-2]。聚合物/表面活性剂驱是目前应用最为广泛的化学驱方式之一, 大庆油田、胜利油田等均进行了大量的矿场实践并取得了较好的开发效果, 胜利油田聚表二元驱矿场应用统计结果表明, 聚表二元驱可提高采收率9%~15% ^[3-4]。虽然化学驱可以获得较好的开发效果, 但其同时具有高投资和高风险的特点, 在当前国际油价持续低迷的形势下, 如何对注采参数进行优化从而实现降本增效是化学驱面临的主要问题^[5-6]。目前油田常用的注采参数优化方法主要为单因素分析法或正交设计法^[7]。单因素分析法需要在固定其他参数不变的情况下对单一因素进行对比, 没有考虑各因素之间的相互影响^[8]; 正交设计法对各个因素取不同的水平进行优化, 相比于单因素分析法考虑了因素间的互相干扰, 但该方法只能优化出所选水平的最优组合, 很难得到真实的最优参数组合^[9]。近年来基于无梯度优化算法的油藏生产优化技术得到了越来越多的应用, 该技术通过数值模拟结果分析得到目标函数的近似梯度并通过迭代不断逼近目标函数的最大值, 相比于传统优化技术大大提升了获得全局最优的可能性^[10]。然而目前该技术还主要局限于水驱注采量优化、新井井位选取等方面, 尚未系统建立聚表二元复合驱的注采参数优化方法^[11-13]。为此, 笔者以经济净现值为目标函数, 考虑矿场注入能力及地层条件限制建立约束条件, 建立一套二元复合驱注采参数优化模型, 并采用最优扰动近似梯度算法对优化模型进行求解和实例应用。

1 注采优化模型建立 1.1 二元复合驱数学模型

考虑聚表二元溶液的黏度、残余阻力因子、不可及孔隙体积、扩散、吸附以及表面活性剂对相对渗透率的影响等物化特征^[14-15], 建立二元复合驱三维两相(油相、水相)四组分(油、水、聚合物、表面活性剂)数学模型, 该模型与黑油模型的区别主要在于聚合物和表面活性剂物化特征参数的表征, 化学剂的质量守恒方程为

$ \begin{array}{l} \nabla \cdot \left( {\frac{{k{K_{{\rm{rw}}}}{c_{\rm{c}}}}}{{{\mu _{{\rm{wc}}}}{B_{\rm{w}}}{r_{\rm{k}}}}}\left( {\nabla {p_{\rm{w}}} - {\rho _{\rm{w}}}g\nabla D} \right)} \right) + \\ \nabla \cdot \left[ {{d_{\rm{c}}}\varphi \left( {1 - {F_{\rm{c}}}} \right){S_{\rm{w}}}\nabla {c_{\rm{c}}}} \right] + {q_{\rm{w}}}{c_{\rm{c}}} = \\ \frac{\partial }{{\partial t}}\left[ {\frac{{\varphi \left( {1 - {F_{\rm{c}}}} \right){S_{\rm{w}}}{c_{\rm{c}}}}}{{{B_{\rm{w}}}}}} \right] + \frac{{\partial \left[ {\left( {1 - {F_{\rm{c}}}} \right)\left( {1 - \varphi } \right){\rho _{\rm{b}}}{{\hat c}_{\rm{c}}}} \right]}}{{\partial t}}. \end{array} $

(1)

式中, k为绝对渗透率; K_rw为水相相对渗透率; μ_wc为聚表二元溶液的黏度, mPa·s; B_w为水相体积系数; r_k为水相渗透率下降系数; p_w为水相压力, MPa; S_w为水相饱和度; q_w为源汇项, m³/(d·m³); d_c为扩散系数, m²/s; F_c为不可及孔隙体积; c_c为化学剂质量浓度, kg/m³; $\hat c$_c为吸附质量分数。

由于表面活性剂对聚表二元溶液黏度的影响较小, 为了简化计算, 聚表二元溶液黏度表征模型仅与聚合物质量浓度有关, 其表达式为

$ {\mu _{\rm{wc}}}= {\mu _{\rm{w}}}\left[ {1 + \left( {{a_1}{c_{\rm{p}}} + {a_2}c_{\rm{p}}^2 + {a_3}c_{\rm{p}}^3} \right){{\left( {\frac{\gamma }{{{\gamma _{1/2}}}}} \right)}^{{n_{\rm{p}}} - 1}}{{\left( {\frac{{{c_{\rm{s}}}}}{{{c_{{\rm{smin}}}}}}} \right)}^{{n_{\rm{s}}}}}} \right]. $

(2)

式中, c_p为聚合物溶液质量浓度, kg/m³; μ_w为水相黏度，mPa·s; a₁、a₂和a₃为实验拟合参数; γ为剪切速率, s^-1; γ_1/2为对应于1/2零剪切速率黏度时的剪切速率, s^-1; c_s为盐离子质量浓度, kg/m³; c_smin为最小盐质量浓度, kg/m³; n_p和n_s为对应指数, 可通过实验结果得到。

表面活性剂分子可以大幅降低油水界面张力, 进而增大毛管数, 同时, 表面活性剂可以使储层岩石变得更加水湿, 束缚水饱和度增加。残余油饱和度和束缚水饱和度的表征模型可表达为

$ {S_{{\rm{or}}}} = {S_{{\rm{orh}}}} = \frac{{{S_{{\rm{orl}}}} - {S_{{\rm{orh}}}}}}{{1 + {T_1}{N_{{\rm{co}}}}}}, $

$ {S_{{\rm{wc}}}} = \min \left( {{S_{{\rm{wc}}}},{S_{{\rm{wch}}}} + \frac{{{S_{{\rm{wcl}}}} - {S_{{\rm{wch}}}}}}{{1 + {T_1}{N_{{\rm{co}}}}}}} \right). $

(3)

式中, S_orl为低毛管数、低弹性时对应的残余油饱和度; S_orh为高毛管数、高弹性时对应的残余油饱和度; T₁为实验拟合参数, MPa^-1; N_co为毛管数; S_wcl为低毛管数、低弹性时对应的束缚水饱和度; S_wch为高毛管数、高弹性时对应的束缚水饱和度。

1.2 优化目标函数及约束条件

二元复合驱通过向注入水中添加聚合物和表面活性剂来改善油藏波及效果和提高洗油效率, 从而大幅提高原油采收率。由于二元复合驱注采设备投资大且化学剂注入及产出液处理等操作成本高, 因此采用动态经济评价法评估预测二元复合驱开发方案的经济效益是其矿场实施成功与否的关键。经济净现值是动态经济评价法中最常用的评价指标, 通常指整个项目实施周期内, 各年限净现金流入量折现后求和所得到收益值, 其表达式为

$ J = \sum\limits_{i = 1}^t {\left( {C_{{\rm{in}}}^i - C_{{\rm{out}}}^i} \right){{\left( {1 + {i_{\rm{c}}}} \right)}^{ - i}}} . $

(4)

式中, J为经济净现值, 元; C_inⁱ为第i年的现金流入量, 元; C_outⁱ为第i年的现金流出量, 元; t为项目实施年限, a; i_c为折现率, %。

针对矿场实施二元复合驱问题, 年现金流入量为二元复合驱年增油量的销售收入, 表达式为

$ C_{{\rm{in}}}^i = {Q_{{\rm{oi}}}}{P_{\rm{o}}}\alpha . $

(5)

式中, Q_oi为二元复合驱年增油量, t; P_o为原油销售价格, 元/t; α为原油商品率, %。

二元复合驱年现金流出量包括注采设备投入、聚合物和表面活性剂注入费用、税费以及增油操作成本等, 其表达式为

$ C_{{\rm{out}}}^i = {Q_{{\rm{oi}}}}{C_{\rm{m}}} + {Q_{{\rm{pi}}}}{P_{\rm{p}}} + {Q_{{\rm{si}}}}{P_{\rm{s}}} + \left( {{R_{\rm{s}}} + R{P_{\rm{o}}}} \right){Q_{{\rm{oi}}}}\alpha + n{I_{\rm{s}}}. $

(6)

式中, C_m为增油操作成本, 元/t; Q_pi为聚合物年注入量, t; P_p为聚合物购买价格, 元/t; Q_si为表面活性剂年注入量, t; P_s为表面活性剂购买价格, 元/t; R_s为资源税, 元/t; R为综合税率, %; n为二元复合驱实施注采井数, 口; I_s为实施二元复合驱的单井增量投资费用, 元/口。

假设注采井增量投资仅发生在二元复合驱项目初始阶段, 将式(5)和(6)代入式(4), 则得到二元复合驱项目的增量经济净现值为

$ \begin{array}{l} J\left( \mathit{\boldsymbol{x}} \right) = \sum\limits_{i = 1}^t {\left( {{Q_{{\rm{oi}}}}{P_{\rm{o}}}\alpha - {Q_{{\rm{oi}}}}{C_{\rm{m}}} - {Q_{{\rm{pi}}}}{P_{\rm{p}}} - {Q_{{\rm{si}}}}{P_{\rm{s}}} - \left( {{R_{\rm{s}}} + } \right.} \right.} \\ \left. {\left. {R{P_{\rm{o}}}} \right){Q_{{\rm{oi}}}}\alpha } \right)/{\left( {1 + {i_{\rm{c}}}} \right)^i} - n{I_{\rm{s}}}. \end{array} $

(7)

式中, x为优化参数向量。

矿场实施二元复合驱项目时, 由于注采设备及地层条件限制, 注采参数需满足一定的约束条件, 包括注采液量平衡、化学剂总注入成本恒定、边界约束等, 其表达式为

$ f\left( {{x_i}} \right) = 0,i = 1, \cdots ,{N_u}, $

(8)

$ {x_{i{\rm{low}}}} \le {x_i} \le {x_{i{\rm{up}}}},i = 1, \cdots ,{N_u}. $

(9)

其中, 式(8)表示等式约束; 式(9)表示边界约束。

二元复合驱注采参数优化问题即在满足约束条件的前提下, 优化寻找注采参数使得经济净现值达到最大。

2 注采优化模型的求解 2.1 约束条件处理方法

本文中所建立的注采优化模型主要涉及边界约束和等式约束条件, 其中边界约束将采用对数变换进行处理, 表达式为

$ {s_i} = \ln \left( {\frac{{{x_i} - {x_{i{\rm{low}}}}}}{{{x_{i{\rm{up}}}} - {x_i}}}} \right). $

(10)

式中, s_i为对数变换后对应原始x_i的优化变量。

对于等式约束条件, 本文中采用拉格朗日乘子法进行处理, 构造拉格朗日函数:

$ L\left( {\mathit{\boldsymbol{s}},\lambda } \right) = J\left( \mathit{\boldsymbol{s}} \right) + \lambda f\left( \mathit{\boldsymbol{s}} \right). $

(11)

式中, λ为拉格朗日乘子; s=[s₁, s₂, …, s_Nu]^T为对数变换后对应原始x=[x₁, x₂, …, x_Nu]^T的优化变量。

在二元复合驱注采参数优化过程中, 实际上是采用优化算法对经过对数变换后的优化变量进行迭代计算, 优化结束后再逆变换回真实注采参数，表达式为

$ {x_i} = \frac{{{x_{i{\rm{low}}}} + \exp \left( {{s_i}} \right){x_{i{\rm{up}}}}}}{{1 + \exp \left( {{s_i}} \right)}}. $

(12)

2.2 最优近似梯度算法

由最优化理论可知, 沿目标函数对优化参数的梯度方向迭代计算即可求得目标函数的最大值及对应参数。然而, 二元复合驱注采参数优化是一类非线性复杂问题, 难以求解其真实梯度, 因此实际研究中多采用近似梯度算法。其中有限差分近似梯度精度高, 但计算量正比于优化参数个数, 难以应用于多参数优化问题; 同时扰动随机逼近算法^[16-17]同时求取目标函数对所有优化参数的近似梯度, 大幅减少了每个迭代步的近似梯度计算量, 但同时也造成了近似梯度精度的降低, 增加了计算收敛所需要的迭代次数。为此, 相关学者提出了一种最优近似梯度算法^[18-19], 在较小梯度计算量的基础上提高了梯度近似精度, 减少了迭代次数, 因此本文中选用该方法进行二元复合驱注采优化模型的求解。

假设优化参数扰动向量为γΔ, 则目标函数模拟结果向量形式下应满足:

$ L\left( {\mathit{\boldsymbol{u}} + \gamma \mathit{\boldsymbol{ \boldsymbol{\varDelta} }}} \right) = L\left( \mathit{\boldsymbol{u}} \right) + \gamma {\mathit{\boldsymbol{g}}^{\rm{T}}}{\mathit \Delta} + o\left( {{{\left\| {\gamma \mathit{\boldsymbol{ \boldsymbol{\varDelta} }}} \right\|}^2}} \right). $

(13)

式中, u=[s, λ]^T为对数变换后的优化参数向量; γ为扰动步长; Δ为扰动向量; g^T为L在u处的真实梯度。

当γ较小时, 式(13)可以简化为

$ L\left( {\mathit{\boldsymbol{u}} + \gamma \mathit{\boldsymbol{ \boldsymbol{\varDelta} }}} \right) = L\left( \mathit{\boldsymbol{u}} \right) + \gamma {\mathit{\boldsymbol{g}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}\mathit{.} $

(14)

对应优化参数改变量γΔ的目标函数改变量为

$ \Delta L = L\left( {\mathit{\boldsymbol{u}} + \gamma \mathit{\boldsymbol{ \boldsymbol{\varDelta} }}} \right) - L\left( \mathit{\boldsymbol{u}} \right) = \gamma {\mathit{\boldsymbol{g}}^{\rm{T}}}\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}\mathit{.} $

(15)

根据同时扰动随机逼近梯度的形式, 构造近似扰动梯度的一般表达式为

$ \mathit{\boldsymbol{\hat g}} = \frac{1}{c}\mathit{\boldsymbol{ \boldsymbol{\varDelta} M}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\Delta L. $

(16)

式中, $\mathit{\boldsymbol{\hat g}}$为L在u处的近似扰动梯度; c为常数; M为一个下三角阵。

L在u处的近似扰动梯度$\mathit{\boldsymbol{\hat g}}$与真实梯度g^T的夹角余弦为

$ \cos \left\langle {\mathit{\boldsymbol{g,\hat g}}} \right\rangle = \frac{{{\mathit{\boldsymbol{g}}^{\rm{T}}}\mathit{\boldsymbol{\hat g}}}}{{\left| \mathit{\boldsymbol{g}} \right|\left| {\mathit{\boldsymbol{\hat g}}} \right|}} = \frac{{{{\left\| {{\mathit{\boldsymbol{M}}^{\rm{T}}}\Delta L} \right\|}^2}}}{{c\left| \mathit{\boldsymbol{g}} \right|\left| {\mathit{\boldsymbol{ \boldsymbol{\varDelta} M}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\Delta L} \right|}}. $

(17)

由式(17)可知, 求L在u处的最优近似扰动梯度, 即求解M, 使cos〈g, $\mathit{\boldsymbol{\hat g}}$〉取最大值, 即

$ \max g\left( \mathit{\boldsymbol{M}} \right) = \frac{{{{\left\| {{\mathit{\boldsymbol{M}}^{\rm{T}}}\Delta L} \right\|}^2}}}{{c\left| {\mathit{\boldsymbol{ \boldsymbol{\varDelta} M}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\Delta L} \right|}}. $

(18)

利用最速下降法等最优化方法求解式(18)得到M, 然后代入式(16)即可得到当前迭代步下的最优近似梯度, 然后迭代更新对数变换后的优化参数:

$ \mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{u}}_{{\rm{old}}}} + \alpha \times \frac{{\mathit{\boldsymbol{\hat g}}}}{{{{\left\| {\mathit{\boldsymbol{\hat g}}} \right\|}_\infty }}}. $

(19)

在此基础上, 利用式(6)即可更新得到真实注采参数。

图 1为采用最优近似扰动梯度算法求解二元复合驱注采优化问题时的一般流程。

3 应用实例

基于Eclipse油藏数值模拟软件, 建立二元复合驱数值模拟模型, 如图 2所示。共划分为3个模拟层, 采用直角网格系统进行离散, 包含53×53×3=8 427个网格, 其中x, y方向网格尺寸为10.6 m×10.6 m。油层埋深1 261~1 294 m, 地层压力12.4 MPa, 地层温度68 ℃, 有效厚度12 m, 平均渗透率1 500×10^-3 μm², 渗透率变异系数0.7, 孔隙度34%, 地下原油黏度50 mPa·s, 原始含油饱和度72%。模型采用五点法井网, 包括4口注入井和9口生产井, 油藏整体保持注采平衡, 速度均为0.1V_P/a。油藏水驱至含水率95%时转二元复合驱开发, 当含水率上升至98%时油藏开发结束, 分别统计二元复合驱及相应水驱方案的开发效果, 计算经济净现值。其中经济评价参数包括:单井增量投资100万元/井, 原油销售收入2 500元/t, 吨油增量操作成本600元/t, 聚合物成本20 000元/t, 表面活性剂成本10 000元/t, 资源税率14元/t, 原油商品率97%, 综合税率14%, 基准收益率12%。

本文中主要对二元复合驱各注入井的聚合物质量浓度、表面活性剂质量分数和段塞尺寸进行优化, 共包括12个优化注入参数。考虑二元复合驱矿场实施特点和注采条件限制, 初始均匀注剂方案中各注入井聚合物质量浓度设置为1.5 g/L, 表面活性剂质量分数设置为0.3%, 段塞尺寸设置为0.1V_P。优化迭代计算中, 聚合物注入质量浓度上、下限分别设置为2和1 g/L, 段塞尺寸上下限分别设置为0.05V_P和0.4V_P, 矿场实际应用中表面活性剂质量分数与聚合物质量浓度需满足一定比例关系, 本次研究设定表面活性剂与聚合物质量浓度比值下限为1.8, 上限为2.0。此外, 优化过程中保证聚合物和表面活性剂总注入成本恒定。

图 3为二元复合驱经济净现值优化迭代曲线。可以看出, 经济净现值在早期迭代阶段快速增加, 但随着迭代次数增多, 二元复合驱经济净现值增速变小并经过41次迭代后最终收敛于3 772万元, 表明基于最优扰动近似梯度算法设计注入方案能够有效提高二元复合驱开发项目的经济效益。

表 1中对比了二元复合驱优化注剂方案与均匀注剂方案。可以看出, 优化后注入井Inj-1和Inj-4的聚合物质量浓度、表面活性剂质量分数和段塞尺寸相比均匀注剂方案均有一定程度增加, 这主要是因为Inj-1和Inj-4通过高渗透条带与生产井Pro-1、Pro-5、Pro-6窜通, 降低了注入水波及效果, 适当提高注剂质量浓度和注入量, 有利用扩大波及系数并提高原油开发效果。而Inj-2和Inj-3井所在区域油藏物性均质性较好, 水驱阶段已取得较好的开发效果, 因此与均匀注剂方案相比, 优化后适当降低Inj-2和Inj-3井的注剂质量浓度和注入量, 并不会严重影响井组开发效果, 而将成本较高的聚合物和表面活性剂分配至其他注入井可以提高油藏整体开发效果。

图 4为二元复合驱优化注剂方案与均匀注剂方案模拟得到的油藏含水率及累积采出程度变化曲线。可以看出, 优化后二元复合驱含水率下降漏斗加深2.53%, 累积采出程度提高1.23%, 即在聚合物和表面活性剂总注入成本不变的前提下, 通过各单井聚合物质量浓度、表面活性剂质量和段塞尺寸的优化分配改善了降水增油效果, 提高了化学剂利用效率。

图 5为采用优化注剂方案后各生产井相对于均匀注剂方案的增油量变化。图 6为二元复合驱优化注剂方案与均匀注剂方案开发结束时的剩余油饱和度分布。可以看出, 在聚合物和表面活性剂总注入成本恒定的条件下, 优化后位于弱非均质区域的Inj-2和Inj-3井注剂质量浓度和注剂总量均有所降低, 从而造成Pro-3、Pro-7和Pro-8井产油量小于均匀注剂方案; 但优化后位于高渗透条带的Inj-1和Inj-4井注剂质量浓度和注剂总量都得到了增加, 从而减轻水驱阶段造成的窜流情况, 一定程度上扩大了波及范围, 从而大幅增加Pro-1、Pro-4、Pro-5、Pro-6和Pro-9井的产油量; 因此, 在化学剂注入总成本未增加的前提下, 油藏整体产油量得到了提高, Inj-1和Inj-4井组所在区域剩余油饱和度降低。

图 7为二元复合驱优化注剂方案与均匀注剂方案模拟得到的经济净现值变化曲线。可以看出, 在聚合物和表面活性剂注入总量一定的情况下, 二元复合驱优化注剂方案经济净现值比各单井均匀注剂方案提高977.8万元, 且投资回收期比各单井均匀注剂方案缩短, 能够使二元复合驱开发方案更早获得经济效益, 减轻油田企业实施二元复合驱项目的负债压力。

4 结论

(1) 以经济净现值为优化目标函数, 考虑矿场注入能力及地层条件限制为约束条件, 建立二元复合驱注采优化模型, 符合矿场实施二元复合驱开发特点, 能够综合反映二元复合驱项目的经济、技术开发效果, 可以用来优化制定二元复合驱注采方案。

(2) 在聚合物和表面活性剂总注入成本一定的情况下, 适当增加水驱窜流较严重区域的聚合物质量浓度、表面活性剂质量分数和段塞尺寸, 并相应减少其他区域的聚合物、表面活性剂配注量, 能够有效提高高价化学剂的利用效率, 改善二元复合驱油藏整体经济和技术开发效果。