简介：该资源面向石油工程领域从事气藏开发与试井分析的研究人员和工程师，围绕裂缝性气藏分支水平井的不稳定试井模型展开，重点解决井底压力动态预测与流动阶段识别问题。压缩包内仅含1个PDF文件，约921KB，集中呈现论文复现所需的完整数学推导、Python代码实现及逐段解释，涵盖Laplace变换解析求解、Stehfest数值反演、典型曲线绘制与参数敏感性分析等核心环节。读者可据此复现9个流动阶段的划分逻辑，理解分支长度、分支间距对第二拟径向流段的影响规律，并借助代码框架完成模型构建、压力动态模拟与现场数据验证。内容还讨论了未考虑应力敏感效应、假设均匀分支流量分布等模型局限，为后续耦合井筒管流、引入应力敏感效应等研究提供切入点。目前已有70人学习，适合具备一定试井理论基础、希望深入掌握分支水平井压力动态分析方法的读者参考。
1. 裂缝性气藏分支水平井试井模型：从论文公式到可运行代码的完整复现路径
裂缝性气藏的分支水平井试井解释，是近两年油气藏工程方向被反复搜索的硬骨头。原因很直接：碳酸盐岩、页岩这类双重介质气藏里，天然裂缝和基质两套系统渗流特征差异巨大，再叠加分支水平井的多段井筒结构，压力响应曲线会出现多个径向流段和过渡段，常规直井试井图版根本套不上。我最初接触这个方向时，对着论文里的 Laplace 空间解和 Stehfest 数值反演公式看了整整一周，公式能看懂，但落到代码上压力导数曲线就是不对——要么早期续流段消失，要么晚期拟径向流斜率不是 0.5。后来把每个中间变量拆开单独验证，才把整条链路跑通。这篇笔记就是把我踩过的坑和最终能复现的代码路径完整讲清楚，适合做气藏工程数值模拟、试井解释方向的研究生和一线工程师，也适合需要快速把论文模型转成可调参工具的人。
2. 双重介质 + 分支水平井：模型到底在算什么
2.1 物理模型拆解：三区两系统一井筒
裂缝性气藏的核心特征是双重介质：基质岩块是主要储集空间，天然裂缝是主要渗流通道。当分支水平井钻入这类储层，压力波传播会经历几个阶段：早期井筒储集和表皮效应主导，压力导数出现单位斜率段；随后裂缝系统先响应，压力波沿裂缝快速传播，出现裂缝径向流；接着基质向裂缝发生窜流，压力导数出现下凹的“凹子”；窜流结束后，整个双重介质系统一起响应，进入总系统径向流；最后边界效应出现，压力导数上翘或下掉。
分支水平井的加入让问题复杂一层：每个分支都相当于一个有限导流或无限导流的内边界，分支之间的干扰、分支与主井筒的连接方式、分支长度和夹角，都会改变压力响应形态。论文里通常用源函数法或边界元法建立半解析模型，在 Laplace 空间求解，再用数值反演回到实空间。
我一般会把整个模型拆成三个模块来理解：储层模型（双重介质渗流方程）、井筒模型（分支水平井内边界条件）、耦合求解（Laplace 空间叠加 + 数值反演）。这三个模块任何一个参数设错，最终曲线都会翻车。
2.2 为什么必须用 Laplace 变换 + Stehfest 反演
直接时域求解双重介质分支水平井的渗流方程，解析解几乎不可能拿到，数值解又面临网格划分和计算量的问题。Laplace 变换把时间导数项消掉，偏微分方程变成常微分方程或代数方程，求解难度大幅下降。得到 Laplace 空间的压力解后，再用 Stehfest 数值反演算法回到实空间。
Stehfest 反演的公式不复杂：
p(t) = (ln2/t) * Σ Vi * P_Laplace(si) si = i * ln2 / t其中 Vi 是反演系数，N 通常取 8 到 12。N 太小精度不够，N 太大数值振荡。我实测下来 N=10 对大多数试井曲线够用，N=12 在晚期边界效应段更稳，但计算量翻倍。
常见做法是：先把 Laplace 空间解写成函数，再写一个 Stehfest 反演函数，最后用双对数坐标画压力曲线和压力导数曲线。压力导数用数值差分计算，时间步长取对数均匀分布。
注意：Stehfest 反演对晚期数据敏感，如果 Laplace 空间解在 s 很小时有数值误差，反演出来的晚期曲线会剧烈振荡。建议在反演前检查 s 从 1e-6 到 1e2 范围内解是否光滑。
2.3 分支水平井内边界条件的三种处理方式
论文里对分支水平井的处理通常有三种：无限导流假设、有限导流假设、均匀流量假设。无限导流假设最简单，认为井筒内压力损失为零，整个分支水平井段压力相同；有限导流考虑井筒内摩阻和加速度项，更接近实际但方程复杂；均匀流量假设认为流量沿分支均匀分布，适合早期段分析。
我建议复现时先从无限导流入手，把双重介质部分跑通，再逐步加入有限导流修正。很多论文的“新模型”其实就是在无限导流基础上加了分支干扰项或非均质性修正，核心求解框架没变。
3. 从零写代码：Laplace 空间解 + Stehfest 反演的完整实现
3.1 双重介质储层的 Laplace 空间压力解
双重介质模型常用 Warren-Root 模型描述，基质和裂缝之间的窜流用形状因子控制。在 Laplace 空间，无限大双重介质储层的点源解可以写成：
import numpy as np from scipy.special import kv # 第二类修正贝塞尔函数 def laplace_pressure_double_porosity(s, rD, omega, lam): """ 双重介质无限大储层 Laplace 空间点源解 s: Laplace 变量 rD: 无量纲半径 omega: 储容比 (裂缝系统储容 / 总储容) lam: 窜流系数 """ # 双重介质综合参数 f_s = omega + (1 - omega) * lam / (s + lam) # 修正贝塞尔函数参数 u = np.sqrt(s * f_s) # 点源解（无限大外边界） p_bar = kv(0, u * rD) / s return p_bar这段代码是整个模型的地基。omega控制裂缝系统储容占比，典型值 0.01 到 0.1；lam控制基质向裂缝窜流速度，典型值 1e-6 到 1e-3。f_s是双重介质的综合压缩系数项，当lam很大时f_s趋近 1，退化为单重介质；当lam很小时窜流段明显。
kv(0, u*rD)是第二类零阶修正贝塞尔函数，描述径向渗流。如果外边界不是无限大，比如圆形封闭或定压边界，需要换成对应的边界条件解，通常用贝塞尔函数的组合形式。
3.2 分支水平井的叠加：分段源 + 镜像反映
分支水平井可以看作多个线段源的组合。每个分支分成若干段，每段用一个点源或线源近似，然后在 Laplace 空间叠加。如果考虑井筒储集和表皮，用杜哈美原理叠加：
def laplace_pressure_branch_well(s, rD, omega, lam, CD, S, branch_segments): """ 分支水平井 Laplace 空间压力解（无限导流 + 井储 + 表皮） CD: 无量纲井筒储集系数 S: 表皮系数 branch_segments: 各分支分段坐标列表 """ # 无限导流假设下，各段压力相同，总流量为1 # 先计算无井储无表皮时的压力响应 p_bar_no_wellbore = 0.0 for seg in branch_segments: rD_seg = seg['rD'] q_seg = seg['q'] # 各段流量分配 p_bar_no_wellbore += q_seg * laplace_pressure_double_porosity(s, rD_seg, omega, lam) # 叠加井筒储集和表皮 p_bar = p_bar_no_wellbore / (1 + s * CD * (s * p_bar_no_wellbore + S)) return p_bar这里的关键是branch_segments的构建：每个分支按长度等分或不等分，每段到观测点的无量纲距离rD不同，流量分配q_seg在无限导流假设下通过求解线性方程组得到，保证各段压力相等。实际写代码时，我会先用均匀流量试算，确认曲线形态合理后再切换到无限导流求解。
井筒储集和表皮的叠加公式是标准形式，CD典型值 1e-3 到 1e2，S典型值 -2 到 10。注意CD太大会把早期段完全掩盖，双对数曲线上看不到裂缝径向流特征。
3.3 Stehfest 数值反演：系数表 + 向量化计算
Stehfest 反演的系数需要预先算好。N 取偶数，系数公式涉及阶乘和求和，直接算容易溢出，建议用查表或对数域计算：
def stehfest_coefficients(N): """计算 Stehfest 反演系数""" V = np.zeros(N) for i in range(1, N+1): j_min = (i + 1) // 2 j_max = min(i, N // 2) sum_val = 0.0 for j in range(j_min, j_max + 1): numerator = j ** (N // 2) denominator = (N // 2 - j) # 阶乘项用对数域计算避免溢出 log_fact = 0.0 for k in range(1, N // 2 - j + 1): log_fact += np.log(k) # 组合数项 log_comb = 0.0 for k in range(1, j + 1): log_comb += np.log(k) for k in range(1, i - j + 1): log_comb -= np.log(k) for k in range(1, 2*j - i + 1): log_comb -= np.log(k) sum_val += np.exp(log_comb + (N//2 - j) * np.log(j) - log_fact) V[i-1] = (-1) ** (i + N//2) * sum_val return V def stehfest_inversion(p_bar_func, t, N=10): """Stehfest 数值反演""" V = stehfest_coefficients(N) ln2 = np.log(2) p = np.zeros_like(t) for idx, ti in enumerate(t): s_i = np.arange(1, N+1) * ln2 / ti p_sum = 0.0 for i in range(N): p_sum += V[i] * p_bar_func(s_i[i]) p[idx] = (ln2 / ti) * p_sum return p这段代码里stehfest_coefficients用对数域计算阶乘和组合数，避免 N=12 时整数溢出。stehfest_inversion对每个时间点计算 N 个 Laplace 变量下的压力值，加权求和。
实际跑的时候，p_bar_func就是前面定义的laplace_pressure_branch_well的偏函数。时间t取对数均匀分布，比如np.logspace(-4, 4, 100)。
提示：N=10 时系数表算一次就够，不用每个时间点重算。把
V缓存起来，反演速度能快一个数量级。
3.4 压力导数计算与双对数曲线绘制
压力导数用数值差分，时间步长取对数均匀后，差分公式要小心：
def pressure_derivative(t, p): """计算压力导数 t * dp/dt，用于双对数诊断""" ln_t = np.log(t) ln_p = np.log(p) # 中心差分 dlnp_dlnt = np.gradient(ln_p, ln_t) return dlnp_dlnt # 绘制双对数曲线 import matplotlib.pyplot as plt t = np.logspace(-4, 4, 100) p = stehfest_inversion(lambda s: laplace_pressure_branch_well(s, ...), t, N=10) dp = pressure_derivative(t, p) plt.loglog(t, p, 'b-', label='压力') plt.loglog(t, dp, 'r--', label='压力导数') plt.xlabel('无量纲时间 tD') plt.ylabel('无量纲压力 pD') plt.legend() plt.grid(True, which='both') plt.show()压力导数曲线是试井解释的灵魂。早期单位斜率段对应井筒储集，裂缝径向流段导数水平，窜流段导数下凹，总系统径向流段导数再次水平，边界段导数上翘或下掉。如果曲线形态不对，先检查omega和lam是否在合理范围，再检查 Stehfest 的 N 是否合适。
4. 参数敏感性分析与曲线形态诊断
4.1 窜流系数 lam 和储容比 omega 怎么调
lam和omega是双重介质模型最核心的两个参数。omega决定窜流凹子的深度，omega越小凹子越深；lam决定凹子出现的时间，lam越小凹子出现越晚。
我一般这样调：先固定omega=0.05，lam从 1e-6 扫到 1e-3，看凹子位置移动；再固定lam=1e-5，omega从 0.01 扫到 0.2，看凹子深度变化。扫参时用循环批量计算，画在一张图上对比。
omega_list = [0.01, 0.05, 0.1] lam_list = [1e-6, 1e-5, 1e-4] for omega in omega_list: for lam in lam_list: p = stehfest_inversion( lambda s: laplace_pressure_branch_well(s, rD=1.0, omega=omega, lam=lam, CD=1e-3, S=0, branch_segments=segs), t, N=10 ) plt.loglog(t, pressure_derivative(t, p), label=f'ω={omega}, λ={lam}') plt.legend()实际气藏中，omega通常小于 0.1，因为裂缝系统储容远小于基质。如果拟合出来omega大于 0.3，要么是模型选错了，要么是早期数据质量有问题。
4.2 分支参数对压力响应的影响
分支水平井的分支数量、分支长度、分支夹角都会影响压力曲线。分支越多，早期径向流段越不明显，因为各分支的干扰叠加会掩盖单个分支的特征。分支长度越长，总系统径向流出现越晚。
我做过一组对比：单分支、双分支、三分支，其他参数相同。单分支曲线在早期有明显的裂缝径向流段，双分支开始模糊，三分支几乎看不到。这说明如果现场数据早期段特征不明显，不一定全是储层问题，也可能是井身结构复杂导致的。
分支夹角的影响相对小，除非夹角小于 30 度，分支之间干扰显著增强。实际建模时，如果夹角数据不确定，可以先按 90 度算，再做敏感性分析。
4.3 井筒储集 CD 和表皮 S 的剥离
井筒储集和表皮主要影响早期段。CD越大，早期单位斜率段越长；S越大，早期压力曲线整体上移。
剥离方法是：先在双对数曲线上找到单位斜率段结束的时间，反推CD；再用早期压力值与理论值的偏差估算S。如果CD和S耦合太强，固定一个调另一个，直到早期段拟合满意。
常见错误是把CD设得太大，导致裂缝径向流段被完全掩盖，然后误以为模型不对。我一般先设CD=1e-3，拟合后再逐步放大。
5. 避坑与排查：复现论文模型时最容易翻车的五个地方
5.1 压力导数曲线晚期振荡
现象：双对数曲线晚期段压力导数剧烈振荡，无法判断边界效应。
原因：Stehfest 反演的 N 太大，或者 Laplace 空间解在 s 很小时有数值误差。
解决：先把 N 降到 8 试，如果振荡消失，说明是反演问题；如果还在，检查laplace_pressure_double_porosity里kv(0, u*rD)在u*rD很小时的数值稳定性，必要时用渐近展开替代。
5.2 窜流凹子不出现
现象：压力导数曲线上看不到双重介质的特征凹子。
原因：lam太大，窜流段和总系统径向流段合并；或者omega太接近 1，双重介质退化为单重介质。
解决：把lam降到 1e-6 量级，omega降到 0.05 以下，重新计算。如果还不出现，检查f_s的表达式是否正确。
5.3 分支干扰导致早期段异常
现象：早期压力导数斜率不是 1，也不是 0.5，而是介于两者之间。
原因：分支之间的干扰在早期就显现，多个分支的井筒储集叠加。
解决：这是物理现象，不一定是错误。但如果斜率异常到无法解释，检查分支分段是否太粗，增加分段数看曲线是否收敛。
5.4 无量纲化参数不一致
现象：曲线形态对，但数值和论文对不上。
原因：无量纲定义不同。有的论文用rD = r/rw，有的用rD = r/L，tD的定义也五花八门。
解决：复现前先把论文的无量纲定义抄下来，和自己的代码逐项对照。我一般会在代码里加注释，标明每个变量的定义来源。
5.5 Stehfest 系数计算溢出
现象：N=12 时系数计算出 NaN 或 Inf。
原因：阶乘和组合数直接计算溢出。
解决：用对数域计算，或者直接查表。网上有现成的 Stehfest 系数表，N=8 到 16 都有，复制过来用最省事。
6. 进阶技巧：用拟合残差自动调参 + 现场数据验证
模型跑通后，下一步是用现场数据验证。手动调参效率太低，我一般用scipy.optimize.minimize做自动拟合。目标函数是压力曲线和压力导数曲线的加权残差：
from scipy.optimize import minimize def objective(params, t_obs, p_obs, dp_obs): omega, lam, CD, S = params p_calc = stehfest_inversion( lambda s: laplace_pressure_branch_well(s, rD=1.0, omega=omega, lam=lam, CD=CD, S=S, branch_segments=segs), t_obs, N=10 ) dp_calc = pressure_derivative(t_obs, p_calc) # 压力和压力导数加权残差 res_p = np.sum((np.log(p_calc) - np.log(p_obs)) ** 2) res_dp = np.sum((np.log(dp_calc) - np.log(dp_obs)) ** 2) return res_p + 0.5 * res_dp # 初始猜测 x0 = [0.05, 1e-5, 1e-3, 0] bounds = [(0.001, 0.5), (1e-8, 1e-2), (1e-5, 1e2), (-5, 20)] result = minimize(objective, x0, args=(t_obs, p_obs, dp_obs), bounds=bounds, method='L-BFGS-B')权重 0.5 是经验值，压力导数对参数更敏感，但噪声也更大。如果现场数据噪声大，把权重降到 0.2 到 0.3。
自动拟合的坑在于初值选择。omega和lam的初值如果偏离太远，优化器容易陷入局部最优。我一般先用网格搜索粗调，再用minimize精调。网格搜索时omega取 0.01 到 0.2，lam取 1e-7 到 1e-3，各取 10 个点，跑 100 组，选残差最小的作为初值。
现场数据验证时，还要注意单位换算。现场压力单位是 MPa，时间单位是小时，代码里是无量纲量，换算系数必须统一。我吃过一次亏：压力换算系数少乘了一个 0.1，拟合出来的omega偏大 10 倍，曲线形态看着对，参数完全没意义。
最后一个习惯：每次拟合完，把拟合曲线和实测曲线叠在一起看，不要只看残差数值。残差小不代表拟合好，有时候曲线整体偏移但残差也小，那是参数补偿的结果。我一般会检查三个点：早期单位斜率段、窜流凹子底部、晚期径向流段，这三个点都对上了，参数才可信。
希望帮到你。
本文还有配套的精品资源，点击获取
