低渗透油藏CO2驱试井解释方法

苏玉亮，蔡明玉，孟凡坤，范理尧，李 蕾

（1.中国石油大学（华东）石油工程学院，山东青岛 266580；2.中国石油大学胜利学院，山东东营 257061）

受制于低渗透油藏水驱注入能力，CO2驱已逐渐成为提高低渗透油藏采收率的主要方法之一，具有广泛的应用前景［1-4］。在CO2注入过程中开展试井测试，可了解油藏性质及井的状况，对指导后续生产具有重要意义［5-6］。目前，针对CO2驱，提出的多为复合油藏试井模型，TANG等率先揭示了CO2驱流体分布特征，将油藏划分为CO2区、CO2-原油过渡区及未波及原油区［7］；随后，AMBASTHA 等建立了三区复合试井解释模型，分析试井曲线变化规律［8］，ISSAKA 等假定过渡区流体性质为指数变化，给出了变性质过渡区复合模型半解析解［9］。在中国，王敬瑶运用传统三区复合模型，分析某CO2驱区块试井数据，但未考虑CO2驱流体分布特征［10］；朱建伟等基于两区复合模型，考虑流体性质变化所产生压降，但试井曲线无明显变化规律［11］；苏玉亮等充分考虑了CO2驱替过程中流体幂律型、非连续型等变化方式，建立了三区复合模型，分析其对试井曲线的影响［12-14］。近年来，众多学者又考虑应力敏感、水平井、垂直裂缝直井等复杂条件，建立多区复合模型，研究试井曲线变化特征［15-17］。

针对目前研究集中于试井曲线分析，而缺乏对CO2驱试井解释方法研究的现状，考虑启动压力梯度、应力敏感等低渗透油藏特性以及流体性质变化，建立改进复合区试井模型，并综合SPSA 自动搜索优化算法，最终形成CO2驱试井解释方法流程。

1 CO2驱试井模型建立与求解

1.1 物理模型

由封闭边界低渗透油藏CO2驱三区复合模型剖面（图1）可见，一区为驱替波及区（CO2区）；二区为CO2-原油过渡区（CO2+原油轻组分），该区流体性质与流体组分及油藏压力、温度等紧密相关，难以用特定函数统一表征；三区代表驱替未波及区（原油区）。温度（304.14K）和临界压力（7.382MPa），故注入CO2实际处于超临界状态，其性质类似于液体，因此可简化模型为单相渗流。同时考虑低渗透油藏非达西渗流的特点，在模型中加入表征低渗透特点的

一般地层条件下的温度和压力均超过CO2临界启动压力梯度及地层压力敏感性。模型的基本假设条件如下：①地层为水平、等厚且恒温的封闭地层，油藏原始地层压力在各处相等。②注入井CO2注入量恒定。③忽略重力和毛管力的影响，考虑表皮污染和井筒存储效应。④一区、三区中流体均微可压缩，且流体性质均一。⑤二区中流体黏度及压缩系数随半径呈连续（幂律、指数）、非连续等多种变化形式。⑥整个流动过程为单相渗流，受流体性质差别影响，各区启动压力梯度不同，均服从非达西渗流规律，渗流速度表达式为：

⑦考虑CO2注入过程中地层压敏效应，其渗透率变化可表示为：

1.2 数学模型

根据物理模型的假设，推导渗流数学模型并进行无因次化后可得：

模型初始条件为：

（3）式中，在R1D＜rD＜R2D的条件下，V，W，Q的定义式分别为：

当黏度及压缩系数随半径成解析函数形式变化时，对于给定的rD可以直接计算出该半径处的V，W，Q值；对于黏度及压缩系数随半径成不规则变化时，对于给定的rD可以采用插值法（如样条插值）计算出该半径处的V，W，Q值。

（3）—（8）式中其他各无因次参数定义式为：

1.3 模型求解

借鉴文献［14］中所采用的求解方法，首先对渗流方程及边界条件（（1）—（4）式）进行对数坐标的转化，并对主控方程（（3）式）进行离散可得：

设一区、二区及三区网格数分别为N1，N2，N3，则对于一区1＜i＜N1：

对于三区N2＜i＜N3：

内边界（i=1）离散为：

其中：

一区-二区（i=N1）、二区-三区（i=N2）交界面及外边界（i=N3）离散方程分别为：

对（27）—（30），（32）—（34）式进行统一整理，可得到线性代数方程组：

利用Newton-Raphson方法构造牛顿迭代矩阵：

即可求解关于压力的非线性方程组：

2 SPSA搜索优化算法

SPSA（Simultaneous Perturbation Stochastic Approximation）方法最先由SPALL［18］于1992年提出，与传统的有限差分近似法（梯度法）不同，该方法通过设置一定的随机变量，对所有控制变量进行同步扰动，以此来获得扰动梯度，进行损失函数等性能指标函数的计算。尽管扰动梯度是随机的，但通过选择适当的随机变量，可以保证对于最小化问题来讲恒为下山方向，且期望值为真实梯度。算法的计算流程包括以下6 个步骤［19］：①模型参数初始化及选择。设定计数标记k=0。选择初始估计值θ0，根据迭代系数αk和ck计算表达式（αk=a/（A+k+1），ck=c/（k+1）γ），确定非负系数a，c，A，α和γ。②随机扰动向量生成。运用Bernoulli±1 分布生成p维扰动向量Δk，为提高早期迭代稳定性，采用多次计算后的平均梯度。③目标函数估计。基于目前值附近的同步扰动，根据步骤①和②中设定的ck与Δk，计算两者目标函数值。④梯度的近似计算。生成同步扰动变化梯度，其表达式为：

gk()中p个变量反映中所有元素同步扰动；⑤参数迭代更新。对进行更新的表达式为：

⑥迭代终止条件判断。返回到第②步，运用k+1 代替k。当相邻2 次迭代步骤间计算结果差值较小，满足设定阈值，或达到最大迭代步数时终止运算。

3 试井解释与自动拟合方法

对于低渗透油藏CO2驱的试井解释，商业试井软件不易考虑其低渗透及CO2驱的特性。因此，可构建低渗透油藏CO2驱试井解释自动拟合方法，其基本思路为：先通过常规试井解释方法得到初始参数值，再通过自动搜索算法结合改进的模型进行拟合矫正。其具体流程如下：①拟定初始解释参数，可采用常规试井解释软件，得到基于均质复合区模型的解释值作为初始值。②根据矿场、实验或者数值模拟数据，给定启动压力梯度、压敏系数以及过渡区参数的变化形式。③根据上述参数重新计算压力及压力导数曲线。④计算目标函数。⑤利用SPSA 算法在变量定义域(xmin,xmax)中进行搜索，使得目标函数值最小。⑥调整搜索步长和搜索次数，使所求的x满足精度要求，并符合实际情况。⑦绘制优化后的压力及压力曲线图。

目标函数为：

目标函数采用的是对压力曲线和压力导数曲线同时进行拟合。由于目前已有的低渗透油藏CO2实际试井数据极为有限，因此这里采用文献［3］中CO2试井数据进行算例示范操作。在该实例中，因采用压力恢复试井，测试压力随时间的增长而增大，初始增大速率较小，后期逐步稳定，由于3 个区间存在界面，使得压力增大速率出现局部突变。通过商业软件采用二区复合油藏模型（表1）初步拟合参数，初步试井参数解释结果为：半径为31.3 m，表皮系数为2.72，井储系数为0.004 4，流度比为1.09，导压系数为0.37，渗透率为31.3 mD，指数I为0.7，所得初始拟合曲线见图2。

将基础数据及初始值导入自动拟合程序。假设该油藏为低渗透油藏，且注二氧化碳时启动压力梯度为0.01 MPa/m，地层压敏系数为0.1 mD/MPa。同时，假定在上述二区模型中存在过渡带。由初步试井参数解释结果可以看出，两区流度比值相差较小，因此黏度的变化在此算例中忽略不计，仍设为恒定值。在过渡区域内压缩系数的变化可简化为呈指数形式变化（由组分模型计算归纳得到），即：此时所得优化前曲线见图2。

表1 油藏参数Table1 Reservoir parameters

图2 优化曲线前后对比Fig.2 Comparison of curves before and after optimization

由于所需搜索的各个参数在数量级以及数值上具有较大差异，因此，为了减小搜索范围，提高搜索速度，采用对各参数的变化倍数进行搜索的方法。即初始搜索值为全为1的向量。根据油藏实际，可给定参数变化倍数的上下限。最终优化所得曲线见图2，其最小二乘值变化如图3所示。

图3 优化过程中最小二乘值变化Fig.3 Change of target least-square value during optimization

根据最终优化的向量，可得出优化后的数值（表2）。由表2 可见，在该算例中，渗透率、井储系数、导压系数以及指数I的变化最为明显。除了井储系数，其余参数都与油藏低渗透特性及非均质流体有直接关系。当存在启动压力梯度及压敏效应时，理论曲线会向上抬高；而渗透率增大时，会使理论曲线相对下降。因此优化后，渗透率解释值增大。而当考虑非均质压缩系数流体时，指数I与导压系数与其直接相关，因此其参数解释值变化较大。

表2 优化后参数变化值Table2 Parameters after optimization

4 结论

考虑低渗透油藏CO2驱替过程中存在启动压力梯度以及压敏效应的特点，同时对CO2驱所导致的流体性质变化运用解析或插值函数进行描述，建立了低渗透油藏CO2驱非均质流体三区复合模型，并进行了数值求解。

基于改进的复合模型，结合SPSA 搜索优化算法，形成了一套低渗透油藏CO2驱试井解释方法，并对实际数据进行了算例计算。

参数修正解释结果表明，当考虑低渗透油藏低渗透特性时，对储层渗透率的参数解释影响较大；当考虑流体非均质分布时，与之相关的指数系数和导压系数等的参数解释值会受到较大影响。

符号解释

rw——井筒原始半径，cm；R——区域半径，m；下标1，2，3——一区、二区和三区；v——流体渗流速度，cm/s；K——地层渗透率，mD；μ——流体黏度，mPa·s；∇p——压力梯度，10-1MPa/cm；G——启动压力梯度，10-1MPa/cm；ε——压敏系数，10-1MPa/cm；p——压力，10-1MPa；下标D——无因次；r——储层径向半径，cm；S——井筒表皮系数；t——时间，s；V，W，Q——表征区域二中流体性质变化的中间变量，无因次；η12，η13——一区和二区、一区和三区导压系数比；C——井储系数，cm3/10-1MPa；pw——井底压力，10-1MPa；M12，M23——一区和二区、二区和三区流体流度比；Ct——地层压缩系数，10 MPa-1；h——地层厚度，cm；q——注入量，cm3/s；Bg——注入流体的体积系数，cm3/cm3；pi——原始地层压力，10-1MPa；ϕ——地层孔隙度；i——网格数，1＜i＜N3；j——时间步；x——对数坐标转换结果，x=lnr；F——残差矩阵；J——系数矩阵；P——压力矩阵；gk()——同步扰动变化梯度——第k次迭代的估计值；Δki——向量Δk中的第i个元素（可能是随机变量±1）；——一区外边界处的压缩系数，10 MPa-1；I——压缩系数变化指数，初始值设为0.7。