基于分阶段拟合的多重介质试井解释方法

A Stage Wise Fitting-Based Well-Test Interpretation Method for Multiple Media

摘要: 由于缝洞型碳酸盐岩油藏显著的非均质性与多重渗流机理的并存，使得常规模型化的试井解释模板一直面临着试井解释模板有限的局限性与拟合失配等问题。针对此问题，本文提出了一种基于分阶段拟合的多重介质试井解释方法。该方法首先在理论上建立了包含基质、裂缝与洞穴多相介质以及多类边界条件的数学模型，并根据数学模型绘制出相应的特征响应曲线，建立响应曲线库；随后对现场试井曲线按流动阶段进行分段划分，并对各流动阶段分别与响应曲线库中相对应的理论模板进行参数化拟合；最终通过组合优化将拟合后的模板按时间序列重构为完整试井响应曲线，并在此基础上开展储层判断与参数反演。根据数值算例表明，该方法在曲线拟合精度和局部响应的判识能力方面均优于基于整体固定模板的传统方法。该方法可有效降低固有模板选择带来的局限性，提高对缝洞型碳酸盐岩复杂渗流行为的识别与定量估计能力。

Abstract: The pronounced heterogeneity of fracture-vug carbonate reservoirs and the coexistence of multiple flow mechanisms have long limited conventional template-based well-test interpretation, causing template constraints and poor fits. To address this problem, this paper proposes a staged-fitting, multi-medium well-test interpretation method. First, a theoretical mathematical model is established that incorporates matrix, fracture and cavity media and multiple types of boundary conditions, and the corresponding characteristic response curves are generated to form a response-curve library. Field well-test curves are then segmented according to flow stages, and each segment is parametrically fitted to the corresponding theoretical template from the library. Finally, by means of combinatorial optimization the fitted templates are reconstructed in time sequence to produce a complete well-test response curve, on the basis of which reservoir classification and parameter inversion are carried out. Numerical examples show that this method outperforms traditional fixed-template approaches in both curve-fitting accuracy and the ability to identify local responses. The proposed approach effectively reduces the limitations associated with fixed template selection and improves the recognition and quantitative estimation of the complex flow behaviours in fracture-vug carbonate reservoirs.

1. 引言

缝洞型碳酸盐岩储层具有复杂的介质结构，既包含微溶孔、裂缝，也存在米级至十米级的大型溶洞，储集空间尺度跨越广泛。这种多尺度的孔–缝–洞结构导致储层整体非均质性较强、孔隙连通性差[1]。此外，储层中流体运动通常不满足传统达西定律，在不同条件下可能出现高速达西流、地下自由流或管流等多种流动形式。这使得缝洞型储层的渗流行为非常复杂，给动态参数解释和储层表征带来巨大挑战。

传统的试井解释方法主要是基于双对数图版或类型曲线拟合，即将井底压力和时间曲线与一系列预设的典型特征曲线拟合，以识别油藏类型并反演储层参数[2]。该类方法在均质或近均质储层中已被广泛验证并取得较好效果，但在多重介质缝洞型储层应用中存在许多不足：一方面，现有模板数量与类型有限，多数模板基于单一或双介质假设，难以覆盖多重介质并存的多模态响应；另一方面，匹配过程较多依赖经验判断，自动化程度较低，且在存在大型空洞或显著非达西流动时，基于等效介质的简化处理将引入明显偏差，从而降低解释精度并增大不确定性[3]。

基于上述背景，本文提出并验证了一种针对缝洞型碳酸盐岩的分阶段拟合试井解释策略。该策略的核心在于：在理论上建立了包含基质、裂缝与洞穴多相介质以及多类边界条件的数学模型，并根据数学模型绘制出相应的特征响应曲线，建立响应曲线库，将压力–时间响应按流动阶段进行分段识别，并将各分段与响应曲线库中的理论模板逐段拟合，从而实现对不同流动模式的分离与定量化。具体流程包括：首先建立覆盖基质、裂缝与溶洞多重介质以及常见边界条件的理论试井曲线库；接着对现场试井曲线按不同的相应阶段进行曲线划分并分别与理论试井曲线库相对应阶段进行参数化拟合；再通过组合拟合度较高的分段理论试井曲线将拟合得到的曲线重构为完整的响应曲线；最后在完整相应曲线基础上开展油藏流动阶段划分和参数反演，并以数值模型对结果进行校核。该方法通过阶段性解耦不同渗流模式，显著降低了参数耦合与多解现象，增强了解释过程的可解释性。

2. 研究现状

Warren与Root提出的双孔介质理论被视为天然裂缝储层试井解释的基础，该理论将储层视为高导通裂缝与低渗透基质两类连续体，提出了储量比与形状因子等关键概念，并基于线性扩散方程构建了类型曲线分析框架[4]。该解析框架在识别裂缝主导的中后期流动特征、基质与裂缝间互流速率方面被广泛采用；但在存在大尺度空洞、显著非线性流动或复杂几何边界的情形下，其等效化假设的适用性受到限制，因而在缝洞并存的碳酸盐岩体系中需谨慎应用并辅以数值校核。

为弥补经典双孔介质解析解在传输细节上的近似，后续研究对形状因子和裂缝传输项进行了系统的数值化扩展[5] [6]。基于有限差分和有限体积等数值方法的研究表明，等效形状因子对几何形状与边界条件敏感，且矩阵–裂缝传输常表现出时间尺度依赖性。因此在具有复杂构造或高度非均质的储层中，形状因子不宜被简化为常数，应作为拟合或校准对象；在工程反演中建议将其作为待估参数，以降低解析近似引入的偏差。

为刻画基质–裂缝–溶洞并存的多模态体系，三重介质模型被提出并逐步发展为重要的理论工具[7] [8]。该模型在连续介质框架下引入第三类孔隙体，并通过多组互渗传输项描述各连续体间的质量交换，从而在压力响应上产生多阶段特征。尽管三重介质模型在辨识多阶流动阶段方面具优势，但其参数数量显著增加且相互耦合，工程应用通常需结合分阶段拟合、先验约束或多级优化以获得稳定可靠的参数估计。

近年综述系统评估了碳酸盐岩中裂缝、洞室的尺度分类、典型渗流机理以及现有试井模型的适用边界[9] [10]。当溶洞尺度达到或超过井筒影响范围的某一临界值时，连续等效模型的近似性将被削弱；此外，重力驱动、局部湍流及孔隙–洞室界面处的非线性效应，均可能在特定时间窗内改变导压曲线的斜率与拐点特征。基于此，应以“尺度对照 + 多源证据驱动”的原则选择模型，即依据测井、成像与钻井记录等现场资料判断是否需要对洞室进行显式建模或采用Darcy-Stokes类耦合方案以提高物理一致性。

近年深度学习与全局优化方法被引入试井解释领域，用以实现曲线自动分类与快速参数估计[11] [12]。这些方法在合成数据和部分现场案例中显示出较高的识别效率与较低的计算延时，适合用于大批量数据的初筛与快速估值；但其对标注样本的依赖性强、对训练集分布敏感且物理约束较弱，导致跨场景泛化能力与工程可解释性存在不足。故在工程流程中，建议将数据驱动方法作为先验识别或初始估计工具，并与物理模型或先验约束相结合以提高可靠性。

3. 多重介质的数学表述与求解框架

3.1. 理论数学模型

油藏渗流模型通常由控制方程、初始条件与边界条件共同构成，用以刻画不同构造介质间的质量守恒与动量传输。对于多重介质体系，三相连续体在达西流动假设和线性互渗律下耦合，显式引入井筒效应和边界约束。为了保持表达的一致性与工程可复现性，本文在统一符号与量纲体系下列出完整的控制方程组及相应的初始–边界条件，随后将对这些方程进行无因次化与拉普拉斯域求解。

井筒储集阶段：

$\Delta P=\frac{qB}{24C}t$ (1)

$\lg \Delta P=\lg t+\lg \frac{qB}{24C}$ (2)

径向流动：

$P_{wf}\left(t\right)=P_i-\frac{2.121qB\mu }{hK}\left(\lg \frac{Kt}{\varphi C_t\mu r_w^2}-2.0923+0.8686S\right)$ (3)

无限导流性垂直裂缝切割井筒：

$\lg \Delta P=\frac{1}{2}\lg t+\lg \left(\frac{0.1959qB}{h{x}_f}\sqrt{\frac{\mu }{\varphi C_{t}K}}\right)$ (4)

有限导流性垂直裂缝：

$\lg \Delta P=\frac{1}{4}\lg t+\lg \frac{1.1054\mu qB}{h\sqrt{K_f\omega \sqrt[4]{\varphi C_{t}K}}}$ (5)

双线性流动：

$\lg \Delta P=\frac{1}{4}\lg t+\lg \left(\frac{1.1054q\mu B}{h\sqrt{K_f\omega }\sqrt[4]{\varphi \mu C_{t}K}}\right)$ (6)

球形流动：

$\Delta P=\frac{0.933q\mu B}{r_{SPH}K}-8.833\frac{q\mu B}{K}\sqrt{\frac{\varphi \mu C_t}{K}}{t}^{-\frac{1}{2}}$ (7)

$r_{SPH}$ -球形流动等价半径。

情况1：长期行为($t\to \infty$ )

当时间t较大时，第二项$t^{-1/2}$ 趋近于0，可忽略。此时公式近似为：

$\Delta P\approx \frac{0.933q\mu B}{r_{SPH}K}$ (8)

对两边取自然对数：

$\lg \left(\Delta P\right)=\lg \left(\frac{0.933q\mu B}{r_{SPH}K}\right)$ (9)

情况2：短期行为($t\to 0$ )

当时间t较小时，第二项主导，第一项可忽略。公式近似为：

$\Delta P\approx -8.833\frac{q\mu B}{K}\sqrt{\frac{\varphi \mu C_t}{K}}{t}^{-\frac{1}{2}}$ (10)

取绝对值后取自然对数(假设$\Delta p$ 为压力降，符号由物理意义决定)：

$\lg \left(\left|\Delta P\right|\right)\approx -\frac{1}{2}v+\lg \left(8.833\frac{q\mu B}{K}\sqrt{\frac{\varphi \mu C_t}{K}}\right)$ (11)

长期：稳态流动(与时间无关)。

短期：瞬态扩散(与$t^{-1/2}$ 相关)。

无限作用径向流动：

$P_{wf}=-\frac{2.121q\mu B}{Kh}\lg t+\left[P_i-\frac{2.121q\mu B}{Kh}\left(\lg \frac{K}{\varphi \mu C_t{r}_w^2}-2.0923+0.8686S\right)\right]$ (12)

均质油藏：

$P_{wf}\left(t\right)=P_i-\frac{2.121qB\mu }{hK}\left(\lg \frac{Kt}{\varphi C_t\mu r_w^2}-2.0923+0.8686S\right)$ (13)

断层边界：

$P=P_i+\frac{0.921\mu qB}{hK}\left[Ei\left(-\frac{r_1^2}{0.0144\eta t}\right)+Ei\left(-\frac{r_2^2}{0.0144\eta t}\right)\right]$ (14)

恒压边界：

$P_{wf}\left(t\right)=P_i-\frac{2.121qB\mu }{hK}\left[\lg \frac{8.085\times {10}^{-3}\eta t}{r_w^2}-\lg \frac{8.085\times {10}^{-3}\eta t}{{\left(2L\right)}^2}\right]$ (15)

层间串流：

$P_D=\frac{K_1{h}_1+K_2{h}_2}{1.842qB\mu }\Delta P$ (16)

3.2. 拉普拉斯变换

在油藏渗流问题的研究中，控制方程通常为非稳态偏微分方程，直接在时域中求解不仅过程复杂，而且在处理边界条件时容易出现困难。为此，常采用拉普拉斯变换将时间域问题转化到复频域中。通过该方法，时间变量t上的一阶或高阶导数可以转化为代数形式，从而显著简化了数学模型的求解过程。

设$f\left(t\right)$ 为定义在$t\ge 0$ 上的时域函数，其拉普拉斯变换$F\left(s\right)$ 定义为：

$F\left(s\right)={\displaystyle {\int }_0^{\infty }f\left(t\right)\cdot {\text{e}}^{-st}\text{d}t}$ (17)

其中，s为复频域变量。利用该变换，可以将时间导数项转化为代数式，例如：

$\frac{L\text{d}f\left(t\right)}{\text{d}t}=sF\left(s\right)-f\left(0\right)$ (18)

这使得控制方程中的非稳态项由时间微分转化为代数运算，极大简化了解析求解难度。

在缝洞型油藏的三重介质模型中，裂缝、基质和溶洞系统的控制方程均包含随时间变化的压力导数项。经过拉普拉斯变换处理后，各介质系统方程可统一转化为拉普拉斯域下的代数方程，进而得到井底压力与井筒径向分布的解析解。此时，径向渗流问题的空间部分通常转化为修正贝塞尔方程，其解为修正第一类与第二类贝塞尔函数$I_0\left(\alpha r\right)$ 、$K_0\left(\alpha r\right)$ 的线性组合。由此可见，拉普拉斯变换在井筒问题中不仅简化了时间项处理，还为求解径向分布提供了明确的数学形式。最终，得到的井底压力及井筒压力导数均以拉普拉斯域表达。为了恢复时域解，还需采用数值反演方法。

3.3. Stehfest数值反演

由于最终结果需要恢复到时域，因此必须对拉普拉斯域解进行数值反演。本文采用Stehfest数值反演方法，其基本公式为：

$f\left(t\right)\approx \frac{\ln 2}{t}\ast {\displaystyle {\sum }_{i=1}^N{V}_i}\cdot F\left(\frac{i\cdot \ln 2}{t}\right)$ (19)

其中：

$F\left(s\right)$ 为拉普拉斯域函数；

$V_i$ 为Stehfest系数；

N为反演项数，一般取N = 8或10。

需要注意的是，N的取值过小会导致反演精度不足，而N过大则可能因舍入误差而引起数值不稳定。在极小或极大时间段内，Stehfest方法可能出现振荡或发散现象，此时可以通过调整N或采用多精度运算进行误差控制。

3.4. 贝塞尔函数的数值计算

在井筒径向流问题中，拉普拉斯变换后的解常包含修正贝塞尔函数I₀、I₁、K₀、K₁等。

当$\alpha r$ 较小时，可采用幂级数展开，例如：

$I_0\left(x\right)\approx 1+\left(\frac{x^2}{4}\right)+\left(\frac{x^4}{64}\right)+\cdots$ (20)

当$\alpha r$ 较大时，可采用渐近展开公式，避免数值溢出。

各阶函数$I_n$ 、$k_n$ 之间还满足递推关系，可用于提高数值计算效率与稳定性。

在实际计算中，需根据参数范围选择合适的计算方法，并通过设定截断条件来控制误差。

对于上述模型的非稳态渗流方程，常采用拉普拉斯变换将时间导数转化为代数形式，从而得到井筒问题在拉普拉斯域(s)下的解。径向二阶空间微分方程可化为修正贝塞尔方程，其通解为：

$P\left(r,s\right)=A\ast I_0\left(\alpha r\right)+B\ast K_0\left(\alpha r\right)$ (21)

其中：

$I_0$ 为第一类零阶修正贝塞尔函数；

$K_0$ 为第二类零阶修正贝塞尔函数；

$\alpha$ 为与渗流参数相关的常数；

A、B为由边界条件确定的常数。

通过拉普拉斯域解法与Stehfest数值反演，可得到井底压力和压力导数随时间的无因次曲线，并绘制对数–对数类型曲线。该曲线可分阶段识别流动机理(井筒效应、裂缝流、溶洞窜流、基质窜流、边界效应等)，并用于参数反演和油藏类型识别。

4. 实例计算

假设单溶洞型储层，半径$r_v$ ，外边界$r_e$ ；流体压缩性线性，可压缩系数$c_t$ ；井口恒流q，外边界恒压$P_i$ 。

控制方程：

$\frac{\partial P\left(r,t\right)}{\partial t}=\frac{k}{\varphi \mu c_t}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial p}{\partial r}\right)$ (22)

初始条件：

$P\left(r,0\right)=P_i$ (23)

边界条件：

$\{\begin{array}{l}-2\pi r_{v}hk{\frac{\partial P}{\partial r}|}_{r=r_v}=q\mu \left(恒流井\right)\\ P\left(r_e,t\right)=P_i\left(外�界恒�\right)\end{array}$ (24)

无因次化：

$P_D=\frac{2\pi kh}{q\mu }\left(P_i-P\right)$ (25)

$r_D=\frac{r}{r_v},t_D=\frac{kt}{\varphi \mu c_t{r}_v^2},r_{eD}=\frac{r_e}{r_v}$ (26)

$\frac{\partial P_D}{\partial t_D}=\frac{1}{r_D}\frac{\partial }{\partial r_D}\left(r_D\frac{\partial P_D}{\partial r_D}\right)$ (27)

对时间$t_D$ 做拉普拉斯变换：

${\overline{P}}_D\left(r_D,s\right)={\displaystyle {\int }_0^{\infty }{\text{e}}^{-s{t}_D}{P}_D\left(r_D,t_D\right)\text{d}{t}_D}$ (28)

得到ODE：

$\frac{1}{r_D}\frac{\text{d}}{\text{d}{r}_D}\left(r_D\frac{\text{d}{\overline{P}}_D}{\text{d}{r}_D}\right)-s{\overline{P}}_D=0$ (29)

得拉普拉斯域解析解：

${\overline{P}}_D\left(r_D,s\right)=A{I}_0\left(\sqrt{s}{r}_D\right)+B{K}_0\left(\sqrt{s}{r}_D\right)$ (30)

$I_0$ ，$K_0$ 为修正贝塞尔函数。

系数A，B由边界条件确定。

恒流井(径向流)：

${r_D\frac{\text{d}{\overline{P}}_D}{\text{d}{r}_D}|}_{r_D=1}=-\frac{1}{s}$ (31)

外边界恒压：

${\overline{P}}_D\left(r_{eD},s\right)=0$ (32)

解得：

${\overline{P}}_D\left(r_D,s\right)=\frac{K_0\left(\sqrt{s}{r}_D\right)-\frac{K_0\left(\sqrt{s}{r}_{eD}\right)}{I_0\left(\sqrt{s}{r}_{eD}\right)}{I}_0\left(\sqrt{s}{r}_D\right)}{s\left[\sqrt{s}\left(K_1\left(\sqrt{s}\right)+\frac{K_0\left(\sqrt{s}{r}_{eD}\right)}{I_0\left(\sqrt{s}{r}_{eD}\right)}{I}_1\left(\sqrt{s}\right)\right)\right]}$ (33)

使用Stehfest数值反演法：

$P_D\left(r_D,t_D\right)=\frac{\ln 2}{t_D}{\displaystyle \sum }_{i=1}^N{V}_i{\overline{P}}_D\left(r_D,s=\frac{i\ln 2}{t_D}\right)$ (34)

$P\left(r,t\right)=P_i-\frac{q\mu }{2\pi kh}{P}_D\left(r_D,t_D\right)$ (35)

5. 影响因素分析

图1展示了井筒储集系数对瞬态压力曲线的影响。在所有曲线在早期阶段基本重合并呈现斜率接近1的直线段，随后随着$C_W$ 的增大，该线性段明显延长，中期特征峰值被延后显现且数值被削弱，谷值位置亦相应右移。较大$C_W$ 条件下的曲线整体在早期表现为相对下移，并在中期压差峰值减小、谷值出现更晚而幅度减弱；而$C_W$ 较小时峰、谷出现提前且数值更接近真实地层响应。这表明井筒储存主要通过掩盖早期地层特征，延迟峰谷的出现时间并降低其幅值，从而影响参数反演的准确性。因此在试井解释时，应优先校正井筒储存效应，否则会低估近井区损伤或高估渗透率。

Figure 1. The effect of wellbore storage coefficient on the curve

图1. 井筒储集系数对曲线影响

图2反映了表皮系数对压力曲线的中期响应特征影响。结果表明，表皮系数增大时，曲线整体向上平移，导致中期压差峰值与谷值的数值显著升高，而峰谷出现的时间位置变化不大。与基准相比，正表皮系数条件下，中期的峰高增加、谷更深，且峰–谷幅差明显增大；而表皮系数减小时，曲线整体下移，峰、谷压差值降低，幅差减小。由于井口附加压降近似与S成正比，该效应在中期尤为突出。因此表皮系数主要控制曲线的纵向位移和中期幅值大小，是解释近井区损伤与完井效率的关键参数。

Figure 2. Effect of skin factor on curve

图2. 表皮系数对曲线影响

图3表明渗透率对曲线的影响。随着K的增大，曲线整体下移，即中期与晚期的压差峰值和谷值显著减小，同时峰、谷的出现时间提前；而当K减小时，曲线整体抬升，峰谷压差增大并延后出现。定量关系上，稳态或准稳态近似下${\Delta }_p\propto \frac{1}{K}$ ，特征时间$t^{*}\propto \frac{1}{K}$ ，因此渗透率提高不仅使峰谷幅值降低，还使其出现更早。图中不同颜色曲线在早期即显著分离，说明K对压力传播速率和数值的影响极为敏感。进一步观察中期谷值深浅可发现，当渗透率较高时，谷值较浅且回升迅速，而低渗透率条件下谷值更深且恢复缓慢，提示裂缝–基质交换或连通性在中期有明显作用。综合来看，渗透率是控制峰谷时间与幅值的核心参数，在试井拟合中应通过峰谷的压差数值与时间位移同时约束渗透率大小。

Figure 3. The effect of permeability on the curve

图3. 渗透率对曲线影响

Figure 4. The influence of $P_i$ on the curve

图4. $P_i$ 对曲线影响

图4展示了$P_i$ 对瞬态曲线的影响规律。随着$P_i$ 增大，曲线整体下移，即中期与晚期的压差值减小，表现为峰、谷幅值降低但峰—谷比值保持近似不变；而当$P_i$ 降低时，曲线整体抬升，峰、谷数值同步增加。因此$P_i$ 的变化不明显改变峰、谷出现的时间，但对其数值幅度有成比例的影响。具体而言，较高$P_i$ 井表现为较低峰谷压差，说明产能优越；较低$P_i$ 则对应于更高峰谷压差，暗示产能不足或近井区连通性差。由此可见，$P_i$ 主要决定曲线整体水平与产能表征，反演时可通过峰谷绝对压差幅度直接估算$P_i$ ，现场应用上则应结合完井与井筒条件，必要时通过增大有效流通截面积来提升$P_i$ ，从而在试井曲线上表现为整体幅值下降。

6. 结论

本文提出了一种基于分阶段拟合的多重介质试井解释方法，并在理论建模与数值算例中对此方法进行了验证。该方法在数学上建立了多介质耦合模型，构建了类型响应曲线库，通过按流动阶段分段识别并分别拟合相应模板，再将分段结果时序重构，能够有效解耦多模态流动特征，显著提高整体拟合精度及对局部响应的识别能力。敏感性分析表明，不同物理参数呈现明显的时间尺度分离：井筒储集主要影响早期响应，表皮效应主导中期特征，渗透率决定峰谷出现时刻与幅值；基于此的分阶段拟合可针对性约束主导参数，增强参数的物理可解释性。方法流程的模块化设计便于与类型曲线库、分段识别器及组合优化器集成，从而在工程解释平台中与测井和成像资料联动以提升现场快速诊断与反演的可靠性。但该方法的适用性受限于类型曲线库的覆盖度与现场数据质量，建议在后续工作中开展多井与大规模现场验证、引入贝叶斯/蒙特卡洛等不确定性量化方法，并探索与数据驱动技术的耦合以兼顾计算效率与物理约束。

基金项目

重庆科技大学研究生创新计划项目(项目批准号：YKJCX2420138)。