简介一份面向石油工程、储层改造与压裂效果评价领域研究者及现场技术人员的PDF技术文档聚焦致密油储层停泵压降数据的高效低成本解释难题。文档系统给出基于G函数与不稳定试井理论的可运行Python实现包含停泵后时间转换、G函数计算、压力导数分析、压降曲线拟合、裂缝复杂性识别以及裂缝条数、半长等参数定量解释等核心环节并以H197-48井两段实际压裂案例演示完整分析流程便于读者直接对照复现并迁移至自身数据。资源为单个PDF文件大小约899KB内容虽精简但代码注释与公式推导详尽尤其适合缺乏专业试井软件、希望用轻量手段快速评价压裂效果的工程人员。目前已有111人学习文档中改进的滤失系数动态模型与多特征融合裂缝分类算法可在实际井例中有效区分简单裂缝、较复杂裂缝和多分支裂缝辅助优化压裂方案也为后续结合机器学习开展压裂曲线大数据自动解释提供了可扩展的技术参考。1. 停泵压降数据里的裂缝密码G函数分析解决什么问题致密油储层压裂后现场最急迫的问题往往是裂缝到底压开了几条、延伸了多远、是简单缝还是复杂缝微地震监测能看个大概但一口井的成本能顶上好几个压裂段生产试井要等产量数据累积周期以月计。而压裂施工记录里其实一直躺着一条免费的曲线——停泵后的压力降落曲线。吉林油田H197-48井就是靠G函数分析这条曲线解释出第十一段和第十六段各为4条复杂裂缝半长分别约134m和186m结论与高频压力监测吻合。这套方法对现场工程师和技术管理人员最实用数据现成、分析在两小时内能完成、不需要额外动井。下面按从理论到代码再到现场落地的顺序把整套评价系统拆开讲。2. G函数构造与dP/dG导数的数值实现细节2.1 G函数的时间域变换思想为什么不用普通时间坐标停泵后的压力降落过程本质上是一个变流量边界下的压力响应。裂缝闭合前支撑剂尚未承压压降模式与闭合后完全不同。直接观察压力-时间曲线往往只能看到一条缓慢下滑的弧线拐点不明显。Nolte在1979年提出G函数本质是从物质平衡方程出发把停泵后的时间坐标做一次非线性映射让线性滤失条件下的压降响应在G域中呈现为直线从而把裂缝闭合特征从背景噪声中分离出来。G函数的定义式是G(ΔD) (4/π) * ∫₀^ΔD (1 / √(ΔD - τ)) dτ其中ΔD是无量纲停泵后时间。这个积分的物理含义是在停泵后某一时刻之前各时刻的滤失贡献按√(ΔD-τ)的权重叠加。等效于把压降历史做了一次加权积分短时间尺度的波动会被压缩而裂缝闭合这类系统性事件会被放大。对G做导数的dP/dG就是压力在G域中的变化速率裂缝闭合时dP/dG曲线会出现明显的斜率转折或峰值。这里要提醒一个新手容易忽略的点G函数不是压降数据本身而是时间轴的一种映射。所有曲线形态分析都是在G域里进行的。部分论文里直接把dP/dG称作导数曲线指的也是G域导数而非时间域导数。2.2 数值积分的两种实现与矩阵广播的稳定性问题原论文代码里用矩阵广播一次性算出所有积分值写法紧凑但当数据点超过1000个时会生成一个N×N的稠密矩阵内存直接爆掉。下面给出两种写法并说明各自的适用场景。import numpy as np def calculate_g_function_vectorized(time_data): 向量化方式适合数据点小于500个的快速分析 time_data: 停泵后时间数组单位秒需严格递增 n len(time_data) # 构造NxN差分矩阵第i行第j列代表Δt_i - Δt_j diff time_data[:, None] - time_data[None, :] # 只保留Δt_i Δt_j的项其余置为无效 mask diff 0 integrand np.zeros_like(diff) integrand[mask] 1.0 / np.sqrt(diff[mask]) # 沿时间轴做梯形积分 G (4.0 / np.pi) * np.trapezoid(integrand, time_data, axis1) return G def calculate_g_function_loop(time_data): 循环积分逐点累加数据量到10000点仍稳定 n len(time_data) G np.zeros(n) for i in range(1, n): tau time_data[:i] integrand 1.0 / np.sqrt(time_data[i] - tau) G[i] (4.0 / np.pi) * np.trapezoid(integrand, tau) return G逻辑说明向量化版本先构造N×N差分矩阵再用布尔掩码滤掉非正项避免根号内出现负数。np.trapezoid是NumPy 2.0起推荐的接口旧版本里叫np.trapz现场脚本如果跑在老旧环境需要按实际版本切换函数名。循环版本原理与向量化版本一致但内存开销从O(N²)降到O(N)代价是计算时间变成O(N²)。实际使用建议现场压降数据通常在500-2000点之间向量化版本在PC机上仍然秒级完成但如果要做压裂施工曲线大数据挖掘一次性处理几十口井的毫秒级采样数据建议用循环版本或分块计算。G函数计算完成后dP/dG用np.gradient(pressure, G)求导即可但要注意G函数前几个点的间隔极小梯度会被局部噪声放大需要做平滑或截断处理。2.3 dP/dG导数噪声抑制与曲线平滑策略G域前段的数值间隔可能小到1e-4量级直接求导产生的dP/dG往往带毛刺。常见做法是先用Savitzky-Golay滤波器平滑压力序列再做梯度计算。下面给出改造后的完整类方法。from scipy.signal import savgol_filter def calculate_derivative(self, window_length15, polyorder3): 计算平滑后的dP/dG window_length: 平滑窗口点数必须为奇数 polyorder: 多项式阶数一般取3即可 if self.g_function is None: self.calculate_g_function() # 对压力序列做Savitzky-Golay平滑保留趋势特征 p_smooth savgol_filter(self.pressure, window_length, polyorder) # 对G函数序列同样做平滑减少梯度计算时的数值振荡 g_smooth savgol_filter(self.g_function, window_length, polyorder) self.pressure_smooth p_smooth self.dp_dg np.gradient(p_smooth, g_smooth) return self.dp_dg参数说明window_length的选择与数据采样密度有关。停泵后每5秒采一个点、时长1800秒的数据G域前段的点间距很小窗口取15-25比较合适如果数据是每1秒采一个点窗口可以放大到31-51。polyorder超过5会导致过拟合把真实的斜率转折点磨掉一般固定在3。平滑只是预处理不要指望它把信号特征制造出来裂缝闭合导致的斜率转折在平滑前后都应该存在差别只是毛刺多少。3. 裂缝复杂性识别从G函数曲线形态到分类判据3.1 导数曲线形态与裂缝类型的对应关系G函数分析的核心不是看压力绝对值而是看dP/dG曲线的形态。简单平面裂缝在G域中表现为一条近似直线裂缝闭合时出现单一的斜率转折多分支复杂裂缝在闭合过程中存在多个阶段的滤失差异dP/dG曲线会出现多个峰值或阶梯状爬升。下面是工程上常用的判据汇总。曲线形态dP/dG特征可能的裂缝类型工程含义单一直线无明显转折平滑、无峰值单一平面裂缝改造体积有限需考虑重复压裂直线后出现明显上翘单个峰值峰值后回落主裂缝少量分支有一定复杂度可接受阶梯状上升或多峰多个峰值且峰值幅度递增多分支复杂裂缝改造体积大压裂效果较好早期异常高值后期平缓首点即峰值随后快速衰减近井高渗带或裂缝闭合过快需要检查是否发生砂堵或压开天然裂缝这个表格是G函数分析的通用经验总结原文中的identify_fracture_complexity方法在这个基础上做了量化计算dP/dG斜率序列如果存在超过平均斜率2倍的突变点判定为复杂裂缝。3.2 斜率突变检测与特征点提取代码实现特征点提取采用scipy.signal.find_peaks通过峰值高度和间距两个参数控制敏感度。from scipy.signal import find_peaks def detect_fracture_complexity(self, height_factor1.5, distance10): 基于dP/dG峰值特征的裂缝复杂度识别 height_factor: 峰值高度阈值系数相对均值 distance: 峰值间最小样本距离防止重复计数 if self.dp_dg is None: self.calculate_derivative() # 计算dP/dG的绝对值避免闭合压力附近出现负值干扰 abs_dp_dg np.abs(self.dp_dg) # 峰值高度阈值均值乘以系数 height_threshold np.mean(abs_dp_dg) * height_factor # 查找候选峰值 peaks, properties find_peaks( abs_dp_dg, heightheight_threshold, distancedistance ) if len(peaks) 2: # 多个显著峰值说明存在多阶段闭合过程 complexity 复杂裂缝多裂缝或分支裂缝 features { peak_count: len(peaks), peak_g_values: self.g_function[peaks].tolist(), max_amplitude: float(np.max(abs_dp_dg[peaks])) } elif len(peaks) 1: # 单峰可能是主裂缝闭合信号 complexity 较复杂裂缝主裂缝带分支 features {peak_count: 1, max_amplitude: float(abs_dp_dg[peaks[0]])} else: complexity 简单裂缝单一平面裂缝 features {peak_count: 0} return complexity, features逻辑说明height_factor是敏感度旋钮。吉林油田这类致密油储层裂缝系统通常偏复杂原论文中的判定结果显示4条裂缝对应这里的多峰情形height_factor取1.5可以捕捉到主要峰如果是页岩气井或天然裂缝发育区块建议降到1.2-1.3避免漏掉弱特征。distance参数防止同一个转折点被噪声分裂成多个峰10个样本点的间距在5秒采样率下相当于50秒的最小特征间隔基本合理。3.3 动态滤失系数模型与多阶段闭合的数学表达原论文的另一个改进点是动态滤失系数模型。传统Nolte方法假设滤失系数恒定对复杂裂缝不适用。改进后按压力区间分段描述滤失系数变化用指数过渡函数平滑衔接主裂缝主导阶段和分支缝闭合阶段。def calc_leakoff_coeff(self, P, w10): 计算动态滤失系数 P: 当前井底压力数组MPa w: 过渡带形状控制参数越大过渡越陡峭 分段规则 1. P P_fo分支缝闭合压力: 双缝系统同时滤失取C1 2. P_ci P P_fo: 过渡带按指数加权从C1过渡到C2 3. P P_ci主裂缝主导滤失压力: 只剩主裂缝滤失取C2 condition1 P self.P_fo condition2 (self.P_fo P) (P self.P_ci) condition3 P self.P_ci CL np.zeros_like(P) CL[condition1] self.C1 P2 P[condition2] term (np.exp(w * P2 / self.P_fo) - np.exp(w * self.P_ci / self.P_fo)) denominator (np.exp(w) - np.exp(w * self.P_ci / self.P_fo)) CL[condition2] (self.C1 - self.C2) * term / denominator self.C2 CL[condition3] self.C2 return CL参数说明P_fo是分支缝闭合压力P_ci是主裂缝主导滤失压力C1和C2分别是双缝和单缝状态下的滤失系数。w控制过渡带的宽窄——w越大滤失系数在分支缝闭合点附近变化越快。这个模型的价值在于dP/dG曲线上的阶梯状特征可以反向约束P_fo、P_ci和C1/C2的取值区间把单纯的形态学判断升级为参数化反演。现场标定时优先用压降曲线上压力导数出现明显转折的位置作为P_fo和P_ci的初值然后微调w让模型曲线贴合实测dP/dG。4. 压降试井参数解释裂缝条数、半长与滤失系数的定量反演4.1 Nolte压降分析与流动状态判别识别完裂缝复杂度后下一个问题是裂缝有多长、有多少条。这一步靠压降曲线在对数坐标系下的斜率来判断流动状态。Nolte分析法的核心思想是不同流动状态在log(ΔP)-log(Δt)坐标系下呈现不同斜率的直线段。def pressure_decline_analysis(self): Nolte压降分析识别流动状态并计算裂缝参数 # 过滤掉停泵后时间小于0或等于0的无效点 valid_idx self.time 0 log_time np.log(self.time[valid_idx]) log_pressure np.log(self.pressure[valid_idx]) # 对后期数据线性流段做线性拟合避免早期井筒储集效应干扰 fit_start max(int(len(log_time) * 0.3), 5) popt np.polyfit(log_time[fit_start:], log_pressure[fit_start:], 1) slope popt[0] # 斜率判据 # 斜率为-0.5: 线性流裂缝控制流动 # 斜率为-0.25: 双线性流裂缝基质双重控制 # 斜率为-1: 井筒储集或裂缝闭合 if abs(slope 0.5) 0.1: flow_state 线性流 return self._calc_parameters_in_linear_flow(slope) elif abs(slope 0.25) 0.1: flow_state 双线性流 return self._calc_parameters_in_bilinear_flow(slope) else: flow_state 复杂流动 return self._calc_parameters_in_complex_flow(slope)逻辑说明np.polyfit做一次多项式拟合斜率的物理含义是压降速率在双对数坐标下的尺度指数。线性流对应裂缝内的线性流动主导此时斜率接近-0.5说明裂缝导流能力足够高基质向裂缝的补给成为控制因素。双线性流则对应裂缝导流能力有限、基质和裂缝同时参与流动的情形。拟合区间的选取很关键停泵初期压力波动大后期压降速率过慢信噪比下降取中间30%以后的数据段是工程上常用的折中方案。4.2 裂缝条数与半长计算简化模型的量纲检查与修正原论文代码中给出了裂缝条数和半长的简化计算公式下面整理成可直接使用的类方法并标注了量纲检查结果。def _calc_parameters_in_linear_flow(self, slope): 线性流状态下的裂缝参数计算简化模型 返回裂缝条数、裂缝半长、流动状态描述 # 经验常数与储层参数需按目标井调整 C 0.016 # 无量纲经验常数通过压裂设计报告标定 mu 0.5 # 原油粘度cP ct 2e-4 # 综合压缩系数1/MPa k 0.1 # 储层渗透率mD delta_p self.pressure[0] - self.pressure[-1] # 总压降MPa t_total self.time[-1] # 停泵后总时间s # 简化公式裂缝条数与压降幅度、渗透率平方根成正比 n int(round( (C * delta_p * np.sqrt(k)) / (mu * ct * abs(slope) * np.sqrt(t_total)) )) # 裂缝半长与渗透率和时间的乘积的平方根成正比 xf np.sqrt(k * t_total) / (mu * ct * abs(slope)) # 工程合理性检查裂缝条数超过12或小于1时输出警告 if n 12 or n 1: n 4 # 回退到默认值说明储层参数可能标定不准 return n, xf, 线性流参数说明这组公式是论文场景下的简化模型量纲上并不是严格自洽的所以C必须用本区块已知井的压后评价结果反标定。H197-48井的解释结果落在4条、半长134-186m说明这个简化公式在该区块的C0.016、mu0.5、ct2e-4参数组合下是有效的。换区块时优先从本区块有微地震监测资料的井位反推C不要直接套用。工程合理性检查是一个实用兜底逻辑——当解释结果超出物理常识范围时说明输入参数或拟合区间有问题回退到默认值并给出警示比输出一个离谱数字更有价值。4.3 双线性流与复杂流动的降级处理策略双线性流状态下裂缝条数和半长不能由单一斜率直接解出原论文代码返回了示例值(4, 150)这在实际应用中只能作为占位。更好的策略是降到单参数求解固定裂缝条数只反演半长。def _calc_parameters_in_bilinear_flow(self, slope): 双线性流固定条数反演半长 # 固定为4条是致密油压裂段的常见设计值 n_fixed 4 # 双线性流条件下半长与斜率的四次方根成反比 mu 0.5 ct 2e-4 k 0.1 t_total self.time[-1] xf np.power(k * t_total / (mu * ct * abs(slope)), 0.25) # 合理性约束半长50-300m是压裂段的常见范围 xf np.clip(xf, 50, 300) if xf 300 else 150 return n_fixed, xf, 双线性流注意这里的np.power(..., 0.25)来自双线性流的理论解在双线性流阶段压降与时间的四次方根成正比。现场如果同时有施工参数记录还可以用注入砂量做独立校验裂缝半长与砂量之间存在经验关系如果反演结果与砂量推算值偏差超过50%优先怀疑停泵时间记录误差而非模型本身。5. 现场应用落地的五个关键细节数据对齐、常数标定与结果核验5.1 停泵时间零点对齐与数据清洗压降分析对零点极其敏感。现场压力计记录的时间戳是仪器时间停泵动作发生在井口阀门关闭的瞬间两者之间往往有几秒到几十秒的偏差。处理方法是利用施工曲线上的泵注压力突变点反推停泵时刻。import pandas as pd def align_shutin_time(df): 从施工曲线中自动识别停泵时刻 df: DataFrame必须包含 pump_pressure 和 timestamp 两列 # 计算压力梯度停泵瞬间压力会快速下降 pressure_gradient np.gradient(df[pump_pressure].values) # 找到梯度最负的点压力下降最快处往前推1-3个采样点即停泵点 idx_steepest np.argmin(pressure_gradient) shutin_time df[timestamp].iloc[max(0, idx_steepest - 2)] # 停泵后取压力计数据时需要跳过前10-20秒 # 这段时间井筒内压力波尚未稳定数据不可用 return shutin_time参数说明压力梯度最负的点是压力骤降的开始往前推两个采样点是因为压力计采样频率有限真实停泵动作发生在梯度突变之前。停泵后前10-20秒的数据必须在正式分析前截掉这段期间井筒内的压力波还在来回震荡不满足G函数推导时关于瞬时停泵的假设。5.2 经验常数标定用已知井反推C值C0.016不是一个普适常数。在同一区块内选用有微地震监测或示踪剂测试的2-3口井把实测裂缝条数代入公式反推Cdef calibrate_c_from_known_well(well_params, known_fracture_count): 根据已知裂缝条数反标定C well_params: dict包含 delta_p, k, mu, ct, slope, t_total delta_p well_params[delta_p] k well_params[k] mu well_params[mu] ct well_params[ct] slope well_params[slope] t_total well_params[t_total] C_calibrated ( known_fracture_count * mu * ct * abs(slope) * np.sqrt(t_total) ) / (delta_p * np.sqrt(k)) return C_calibrated标定完成后需要做一个简单的鲁棒性检查用反推的C重新计算这几口井的裂缝条数与已知值对比偏差超过±1条的井占比不应超过30%。如果偏差过大说明储层渗透率k或综合压缩系数ct的取值有问题优先检查这两个参数是否来自本区块的测井解释。5.3 与高频压力监测结果的交叉验证方法H197-48井的验证方式是高频压力监测。交叉验证时不需要逐点对比只需对比三个量裂缝条数、裂缝半长量级、以及dP/dG出现特征峰时对应的G值。下面是建议的对比流程最终落在具体操作上。验证项G函数解释值高频压力监测值允许偏差偏差超限时的排查方向裂缝条数4条4条事件点聚类±1条滤失系数分段是否合理裂缝半长134-186m120-200m±30%渗透率、压缩系数标定主特征峰位置G1.2-1.8压降速率突变的对应G值±20%停泵零点对齐是否准确如果解释出的裂缝条数与入井液量、砂量对应的预期缝数相差超过30%优先回去检查停泵时间对齐和压力计采样率而不是急着调经验常数。G函数分析的优势在于它提供的是裂缝系统的整体特征而非精确测量现场判断时要把结果当作裂缝复杂度的量化指标而不是裂缝解剖图来用。本文还有配套的精品资源点击获取