简介本资源是一套面向石油地质工程师、测井数据处理初学者及高校地球物理专业学生的测井综合实践工具包聚焦测井数据处理、岩性识别与解释核心能力培养。包内共263个文件以87个C源码cpp和86个头文件h构成主体程序框架支撑测井曲线生成、岩性分类算法实现与解释逻辑封装辅以63幅BMP格式岩性/曲线图示、8个ICO图标及配套资源脚本bat、rc、ini等便于可视化分析与工程集成整体压缩包仅871KB轻量易部署。已有305人下载学习适合需快速掌握从原始测井信号预处理、深度校正、多参数联合岩性判识到典型曲线打印输出全流程的实践者。资源包含完整可编译项目结构含dsp/dsw工程文件、帮助文档生成脚本MakeHelp.bat及多类岩石特征图库如ROCK.BMP、GZZ.BMP为理解测井解释中电阻率-声波-自然伽马三参数协同判别砂岩/泥岩/碳酸盐岩提供代码级支撑。1. 测井数据处理不是“调个模型就出岩性”而是用物理约束统计建模把电阻率、声波、密度曲线翻译成地质语言在油田现场工程师常遇到这样的困境同一口井不同软件输出的岩性分类结果差异显著——砂岩段被标成泥岩含气层被误判为水层。问题不在算法多先进而在于测井曲线本身是间接响应它不直接测量岩性而是记录岩石孔隙中流体与矿物对电磁波、声波的响应。station_95_测井数据处理_测井_岩性曲线_测井解释_测井岩性分类_这个标题指向的是一套以物理驱动为锚点、以数据驱动为增强的闭环流程从原始测井曲线预处理开始构建能反映矿物组合与孔隙结构的岩性敏感曲线如Vsh、PHIE、SW再通过多参数联合判别实现岩性分类。它适合两类人一是刚接触测井解释的地球物理/地质工程师需要可复现的本地化处理链二是已有解释经验但希望将传统交会图法与现代分类器如随机森林、XGBoost融合落地的技术人员。本方案不依赖商业软件许可证所有步骤基于开源工具链Python OpenCV Scikit-learn LASIO重点解决三个真实痛点曲线深度对齐偏差导致的岩性跳变、低信噪比段分类置信度不可靠、以及碳酸盐岩与碎屑岩混合段的矿物解耦困难。2. 用LASIO读取并校正测井曲线深度对齐、坏值剔除与标准化三步不可省测井数据质量直接决定后续所有分析的上限。原始LAS文件常存在深度采样不均、仪器漂移、接箍干扰等问题若直接输入模型分类结果会在井段交界处出现突兀跳变。常见做法是先做物理层清洗再进入统计建模。2.1 用LASIO加载LAS文件并检查深度基准一致性import lasio import numpy as np import pandas as pd # 加载LAS文件注意编码部分老版本用latin-1 las lasio.read(well_A.las, encodingutf-8) # 检查深度索引是否为单调递增且等间距 depth las.index print(f深度范围: {depth.min():.2f} ~ {depth.max():.2f} m) print(f采样点数: {len(depth)}) print(f平均采样间隔: {np.mean(np.diff(depth)):.4f} m) # 检查关键曲线是否存在 required_curves [GR, RT, AC, DEN, CNL] missing_curves [c for c in required_curves if c not in las.keys()] if missing_curves: print(f缺失关键曲线: {missing_curves})提示lasio.read()默认按DEPTH通道作为索引但部分文件使用DEPT或TVD。若报错KeyError: DEPTH需先用las.curves查看实际深度通道名并用las.set_depth_unit(M)统一单位。2.2 基于滑动窗口的坏值检测与插值修复测井曲线中的尖峰spike和平台段flatline会严重干扰后续计算。我们采用双阈值滑动窗口法对每个采样点计算其前后10个点的标准差σ若当前值偏离局部均值超过3σ且该点连续不变超过5个采样点则判定为坏值。def clean_curve(curve_data, depth, window_size10, spike_sigma3, flat_length5): cleaned curve_data.copy() n len(curve_data) for i in range(window_size, n - window_size): window curve_data[i-window_size:iwindow_size1] local_mean np.nanmean(window) local_std np.nanstd(window) # 判定尖峰偏离局部均值 3σ 且非NaN if not np.isnan(curve_data[i]) and abs(curve_data[i] - local_mean) spike_sigma * local_std: # 同时检查是否为平台段连续相同值 if i flat_length n: is_flat np.all(curve_data[i:iflat_length] curve_data[i]) if is_flat: # 平台段用前后窗口均值插值 left_mean np.nanmean(curve_data[max(0,i-5):i]) right_mean np.nanmean(curve_data[i1:min(n,i6)]) cleaned[i] (left_mean right_mean) / 2 else: # 尖峰用线性插值 cleaned[i] np.interp(depth[i], [depth[i-1], depth[i1]], [curve_data[i-1], curve_data[i1]]) return cleaned # 应用清洗以GR曲线为例 gr_raw las[GR] gr_clean clean_curve(gr_raw, depth)2.2.1 参数说明与调优逻辑window_size10对应约0.5米窗口按0.05m采样率太小易受噪声干扰太大则无法捕捉局部异常spike_sigma3符合高斯分布3σ原则对碳酸盐岩中天然放射性异常如钾长石富集保留容忍flat_length5对应0.25米足以过滤仪器停顿导致的平台又不会误删致密灰岩段的真实低GR平台。2.3 多井深度对齐用井壁微电阻率图像FMI或自然伽马GR特征点匹配单井处理后需跨井对比但不同次测井深度系统存在系统偏差可达0.3米。最可靠方法是基于地层标志层对齐。若无FMI图像可用GR曲线的峰值点如煤层、火山灰夹层作为锚点from scipy.signal import find_peaks # 提取GR峰值点要求高度100API宽度3个采样点 peaks, _ find_peaks(gr_clean, height100, distance3) # 获取峰值深度与幅度 anchor_depths depth[peaks] anchor_values gr_clean[peaks] # 假设参考井well_ref的锚点已知计算偏移量 # ref_anchors np.array([1250.3, 1387.6, 1522.1]) # 示例 # offset np.median(anchor_depths - ref_anchors) # 全局偏移 # las.df().index las.df().index - offset # 校正深度索引注意深度对齐必须在生成岩性曲线前完成。若先计算孔隙度再对齐会导致PHIE与AC曲线在深度上错位使Sw计算完全失效。3. 构建岩性敏感曲线从原始测井到Vsh、PHIE、SW的物理公式推导与代码实现岩性分类的根基不是原始曲线而是能表征地质意义的派生参数。station_95流程中Vsh泥质含量、PHIE有效孔隙度、SW含水饱和度是三大核心中间变量其计算必须严格遵循测井物理模型而非黑箱拟合。3.1 泥质含量Vsh用自然伽马GR与声波时差AC双参数交叉验证单一GR法在碳酸盐岩中失效因灰岩本身GR低故引入AC曲线泥岩声波时差高90μs/ft灰岩低50μs/ft。我们采用加权平均法融合两者def calculate_vsh(gr_clean, ac_clean, gr_min20, gr_max120, ac_min45, ac_max110): gr_min/gr_max: 纯砂岩/纯泥岩的GR值需根据区域标定 ac_min/ac_max: 纯砂岩/纯泥岩的AC值同上 # GR法线性刻度 vsh_gr np.clip((gr_clean - gr_min) / (gr_max - gr_min), 0, 1) # AC法同样线性刻度 vsh_ac np.clip((ac_clean - ac_min) / (ac_max - ac_min), 0, 1) # 加权融合AC在致密段更稳定GR在渗透层更敏感 weight_ac 0.7 * (1 - vsh_gr) 0.3 # Vsh越低AC权重越高 vsh_final weight_ac * vsh_ac (1 - weight_ac) * vsh_gr return vsh_final vsh calculate_vsh(gr_clean, las[AC])3.1.1 区域标定参数表需根据本区岩心数据修正岩性类型GR典型值(API)AC典型值(μs/ft)标定建议纯砂岩20–3545–55取试油成功层段的最低GR/AC纯泥岩90–13095–120取大段厚泥岩段的最高GR/AC灰岩5–1540–48若存在单独建立灰岩Vsh模型3.2 有效孔隙度PHIE密度-中子交会法消除岩性影响密度DEN与中子CNL曲线对孔隙流体响应相反但对骨架矿物响应一致。二者交会可消除岩性影响得到真实孔隙度def calculate_phie(den_clean, cnl_clean, den_ma2.65, den_fl1.0, # 砂岩骨架/流体密度 cnl_ma0.0, cnl_fl1.0): # 砂岩骨架/流体中子值 密度-中子交会法PHIE计算假设砂岩骨架 实际应用中需根据岩性切换ma参数如灰岩den_ma2.71, cnl_ma0.0 # 密度孔隙度 phie_den (den_ma - den_clean) / (den_ma - den_fl) # 中子孔隙度 phie_cnl (cnl_fl - cnl_clean) / (cnl_fl - cnl_ma) # 交会法取二者平均但当差异0.05时优先信密度因中子易受氯离子干扰 diff np.abs(phie_den - phie_cnl) phie_final np.where(diff 0.05, phie_den, (phie_den phie_cnl) / 2) return np.clip(phie_final, 0.01, 0.35) # 物理约束孔隙度0.01~0.35 phie calculate_phie(las[DEN], las[CNL])提示若井区存在大量白云岩需修改den_ma2.87, cnl_ma0.05否则PHIE系统性偏低。此参数必须由岩心分析数据标定不可套用邻井。3.3 含水饱和度SW阿尔奇公式与Simandoux公式的工程选择阿尔奇公式Archie仅适用于干净砂岩而Simandoux公式引入Vsh修正项更适合含泥质储层def calculate_sw(rt_clean, phie, vsh, a1.0, m2.0, n2.0, rw0.1, rsh5.0): Simandoux公式1/SW^n (a/PHIE^m) * (RT/RW) * (1-Vsh) (Vsh/RSH) * (1/PHIE^m) * (RT/RW) 参数说明 a,m,n: 阿尔奇参数本区岩心实验标定 rw: 地层水电阻率查SP曲线或实验室数据 rsh: 泥岩电阻率取Vsh0.7段的RT均值 term1 (a / (phie ** m)) * (rt_clean / rw) * (1 - vsh) term2 (vsh / rsh) * (1 / (phie ** m)) * (rt_clean / rw) sw_power 1 / (term1 term2) sw sw_power ** (1/n) return np.clip(sw, 0.15, 1.0) # 设定最小可动水饱和度0.15 sw calculate_sw(las[RT], phie, vsh)3.3.1 关键参数获取路径rw从自然电位SP曲线基线偏移量反算或取邻近水层实测值rsh在Vsh0.7的纯泥岩段取RT曲线的中位数a,m,n必须用本区岩心孔渗数据与测井RT数据回归禁用教科书默认值。4. 测井岩性分类从交会图规则到XGBoost模型的渐进式建模策略有了Vsh、PHIE、SW等物理参数岩性分类就从“经验猜”变为“证据链推理”。station_95流程强调先建规则、再调模型用交会图划定粗粒度岩性区间再用机器学习在边界模糊区精细区分。4.1 三参数交会图Vsh-PHIE-SW构成的岩性立方体将三维参数空间划分为8个卦限每个卦限对应一种主导岩性import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 构建岩性标签数组初始全为Unknown litho_labels np.full(len(vsh), Unknown, dtypeobject) # 定义规则示例基于某盆地标准 litho_labels[(vsh 0.15) (phie 0.18) (sw 0.4)] Clean_Sand litho_labels[(vsh 0.6) (phie 0.08)] Shale litho_labels[(vsh 0.2) (phie 0.06) (sw 0.8)] Tight_Carbonate litho_labels[(vsh 0.3) (phie 0.12)] Sandy_Shale # 可视化交会图 fig plt.figure(figsize(10, 8)) ax fig.add_subplot(111, projection3d) scatter ax.scatter(vsh, phie, sw, c[{Clean_Sand:0,Shale:1,Tight_Carbonate:2,Sandy_Shale:3}.get(l,4) for l in litho_labels], cmaptab10, s1) ax.set_xlabel(Vsh); ax.set_ylabel(PHIE); ax.set_zlabel(SW) plt.show()4.1.1 规则制定的地质依据Clean_Sand低泥质高孔隙低含水 → 典型主力产层Shale高泥质低孔隙 → 非储层但可作盖层Tight_Carbonate低泥质低孔隙高含水 → 白云岩致密层需酸压Sandy_Shale中高泥质中等孔隙 → 水平井靶窗优选区兼顾产能与可压性。4.2 XGBoost模型训练用岩心描述数据监督学习边界模糊区交会图无法区分“泥质粉砂岩”与“粉砂质泥岩”此时需引入岩心描述数据lithology description作为标签from xgboost import XGBClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report # 特征矩阵Vsh, PHIE, SW, GR, RT, AC6维 X np.column_stack([vsh, phie, sw, gr_clean, las[RT], las[AC]]) # 标签需提前准备岩心点深度对应的岩性编码如1Sand, 2Shale, 3Carbonate... # y_core load_core_labels(depth, core_log_file) # 此函数需自行实现 # 划分训练集仅用有岩心的深度点 X_train, X_test, y_train, y_test train_test_split( X[~np.isnan(y_core)], y_core[~np.isnan(y_core)], test_size0.3, random_state42 ) # 训练XGBoost关键参数控制过拟合 model XGBClassifier( n_estimators200, max_depth5, # 防止单棵树过深 learning_rate0.1, # 小步长提升稳定性 subsample0.8, # 行采样防过拟合 colsample_bytree0.8, # 列采样增加泛化 random_state42 ) model.fit(X_train, y_train) # 预测全井段 y_pred_full model.predict(X)4.2.1 模型评估与可信度阈值设定# 输出预测概率用于设定置信度阈值 y_proba model.predict_proba(X) max_proba np.max(y_proba, axis1) # 仅当最高概率0.85时采纳模型结果否则回退到交会图规则 final_litho np.where(max_proba 0.85, y_pred_full, litho_labels_encoded)注意XGBoost的predict_proba输出的是各岩性类别的概率估计非绝对确定性。实践中若最大概率0.7应标记为“Uncertain”触发人工复核而非强行分类。5. 岩性曲线可视化与解释验证用深度轨迹图误差热力图定位分类薄弱区最终输出的岩性曲线必须能被地质师快速验证。station_95流程强制要求双视图输出左侧为深度-岩性轨迹图直观右侧为分类误差热力图可追溯。5.1 绘制标准测井解释图Track Plotimport matplotlib.patches as mpatches def plot_litho_track(depth, gr_clean, rt_clean, phie, vsh, litho_labels, litho_colors{Clean_Sand:yellow,Shale:gray,Tight_Carbonate:blue,Sandy_Shale:brown}): fig, axes plt.subplots(1, 4, figsize(12, 10), shareyTrue) # GR Track axes[0].plot(gr_clean, depth, g, labelGR) axes[0].set_xlim(0, 150) axes[0].set_xlabel(GR (API)) axes[0].grid(True) # RT Track axes[1].semilogx(rt_clean, depth, r, labelRT) axes[1].set_xlim(0.2, 2000) axes[1].set_xlabel(RT (ohm.m)) axes[1].grid(True) # PHIE Vsh Track axes[2].plot(phie, depth, b, labelPHIE) axes[2].fill_betweenx(depth, 0, vsh, color0.7, alpha0.5, labelVsh) axes[2].set_xlim(0, 0.3) axes[2].set_xlabel(PHIE/Vsh) axes[2].legend() axes[2].grid(True) # Lithology Track柱状图 litho_numeric np.array([list(litho_colors.keys()).index(l) if l in litho_colors else 0 for l in litho_labels]) axes[3].fill_betweenx(depth, 0, 1, wherenp.isin(litho_labels, list(litho_colors.keys())), facecolor[litho_colors.get(l, white) for l in litho_labels], alpha0.8) axes[3].set_xlim(0, 1) axes[3].set_xlabel(Lithology) axes[3].set_yticks([]) # 隐藏y轴刻度保持整洁 # 添加图例 legend_elements [mpatches.Patch(facecolorcolor, labellith) for lith, color in litho_colors.items()] axes[3].legend(handleslegend_elements, loccenter left, bbox_to_anchor(1, 0.5)) plt.tight_layout() plt.show() plot_litho_track(depth, gr_clean, las[RT], phie, vsh, litho_labels)5.2 生成分类误差热力图定位模型失效的深度段误差热力图的核心是对比模型预测与交会图规则的分歧点这些点往往是地质复杂区如断层破碎带、流体突变层# 计算规则法与模型法的差异0一致1不一致 rule_vs_model (litho_labels ! [list(litho_colors.keys())[i] for i in y_pred_full]) # 转换为滑动窗口误差率每10米统计一次不一致比例 window_size_m 10 window_points int(window_size_m / np.mean(np.diff(depth))) error_rate [] for i in range(0, len(depth), window_points): window rule_vs_model[i:iwindow_points] error_rate.append(np.mean(window) if len(window) 0 else 0) # 绘制热力图 plt.figure(figsize(10, 6)) plt.imshow([error_rate], cmapRdYlBu_r, aspectauto, extent[0, len(error_rate), depth.min(), depth.max()]) plt.colorbar(labelClassification Disagreement Rate) plt.xlabel(Window Index (10m each)) plt.ylabel(Depth (m)) plt.title(Model vs Rule Disagreement Heatmap) plt.show()5.2.1 误差热力图的地质解读指南红色热点误差率0.6立即检查该深度段的岩心照片或FMI图像常对应断层角砾岩、沥青充填缝洞或钻井液侵入带连续黄色条带误差率0.3~0.6提示该层段岩性过渡渐变需调整交会图边界或增加训练样本全蓝区域误差率0.1模型与规则高度一致可放心用于批量处理。关键技巧在误差热力图中叠加试油结论如“1250–1255m日产油32t”若高产层恰好位于红色热点区说明此处岩性分类虽难但恰恰是优质储层——这正是station_95流程的价值不追求全局准确率而聚焦于识别“高价值不确定性”。本文还有配套的精品资源点击获取 测井岩性分类：物理建模与XGBoost融合的开源实现
简介本资源是一套面向石油地质工程师、测井数据处理初学者及高校地球物理专业学生的测井综合实践工具包聚焦测井数据处理、岩性识别与解释核心能力培养。包内共263个文件以87个C源码cpp和86个头文件h构成主体程序框架支撑测井曲线生成、岩性分类算法实现与解释逻辑封装辅以63幅BMP格式岩性/曲线图示、8个ICO图标及配套资源脚本bat、rc、ini等便于可视化分析与工程集成整体压缩包仅871KB轻量易部署。已有305人下载学习适合需快速掌握从原始测井信号预处理、深度校正、多参数联合岩性判识到典型曲线打印输出全流程的实践者。资源包含完整可编译项目结构含dsp/dsw工程文件、帮助文档生成脚本MakeHelp.bat及多类岩石特征图库如ROCK.BMP、GZZ.BMP为理解测井解释中电阻率-声波-自然伽马三参数协同判别砂岩/泥岩/碳酸盐岩提供代码级支撑。1. 测井数据处理不是“调个模型就出岩性”而是用物理约束统计建模把电阻率、声波、密度曲线翻译成地质语言在油田现场工程师常遇到这样的困境同一口井不同软件输出的岩性分类结果差异显著——砂岩段被标成泥岩含气层被误判为水层。问题不在算法多先进而在于测井曲线本身是间接响应它不直接测量岩性而是记录岩石孔隙中流体与矿物对电磁波、声波的响应。station_95_测井数据处理_测井_岩性曲线_测井解释_测井岩性分类_这个标题指向的是一套以物理驱动为锚点、以数据驱动为增强的闭环流程从原始测井曲线预处理开始构建能反映矿物组合与孔隙结构的岩性敏感曲线如Vsh、PHIE、SW再通过多参数联合判别实现岩性分类。它适合两类人一是刚接触测井解释的地球物理/地质工程师需要可复现的本地化处理链二是已有解释经验但希望将传统交会图法与现代分类器如随机森林、XGBoost融合落地的技术人员。本方案不依赖商业软件许可证所有步骤基于开源工具链Python OpenCV Scikit-learn LASIO重点解决三个真实痛点曲线深度对齐偏差导致的岩性跳变、低信噪比段分类置信度不可靠、以及碳酸盐岩与碎屑岩混合段的矿物解耦困难。2. 用LASIO读取并校正测井曲线深度对齐、坏值剔除与标准化三步不可省测井数据质量直接决定后续所有分析的上限。原始LAS文件常存在深度采样不均、仪器漂移、接箍干扰等问题若直接输入模型分类结果会在井段交界处出现突兀跳变。常见做法是先做物理层清洗再进入统计建模。2.1 用LASIO加载LAS文件并检查深度基准一致性import lasio import numpy as np import pandas as pd # 加载LAS文件注意编码部分老版本用latin-1 las lasio.read(well_A.las, encodingutf-8) # 检查深度索引是否为单调递增且等间距 depth las.index print(f深度范围: {depth.min():.2f} ~ {depth.max():.2f} m) print(f采样点数: {len(depth)}) print(f平均采样间隔: {np.mean(np.diff(depth)):.4f} m) # 检查关键曲线是否存在 required_curves [GR, RT, AC, DEN, CNL] missing_curves [c for c in required_curves if c not in las.keys()] if missing_curves: print(f缺失关键曲线: {missing_curves})提示lasio.read()默认按DEPTH通道作为索引但部分文件使用DEPT或TVD。若报错KeyError: DEPTH需先用las.curves查看实际深度通道名并用las.set_depth_unit(M)统一单位。2.2 基于滑动窗口的坏值检测与插值修复测井曲线中的尖峰spike和平台段flatline会严重干扰后续计算。我们采用双阈值滑动窗口法对每个采样点计算其前后10个点的标准差σ若当前值偏离局部均值超过3σ且该点连续不变超过5个采样点则判定为坏值。def clean_curve(curve_data, depth, window_size10, spike_sigma3, flat_length5): cleaned curve_data.copy() n len(curve_data) for i in range(window_size, n - window_size): window curve_data[i-window_size:iwindow_size1] local_mean np.nanmean(window) local_std np.nanstd(window) # 判定尖峰偏离局部均值 3σ 且非NaN if not np.isnan(curve_data[i]) and abs(curve_data[i] - local_mean) spike_sigma * local_std: # 同时检查是否为平台段连续相同值 if i flat_length n: is_flat np.all(curve_data[i:iflat_length] curve_data[i]) if is_flat: # 平台段用前后窗口均值插值 left_mean np.nanmean(curve_data[max(0,i-5):i]) right_mean np.nanmean(curve_data[i1:min(n,i6)]) cleaned[i] (left_mean right_mean) / 2 else: # 尖峰用线性插值 cleaned[i] np.interp(depth[i], [depth[i-1], depth[i1]], [curve_data[i-1], curve_data[i1]]) return cleaned # 应用清洗以GR曲线为例 gr_raw las[GR] gr_clean clean_curve(gr_raw, depth)2.2.1 参数说明与调优逻辑window_size10对应约0.5米窗口按0.05m采样率太小易受噪声干扰太大则无法捕捉局部异常spike_sigma3符合高斯分布3σ原则对碳酸盐岩中天然放射性异常如钾长石富集保留容忍flat_length5对应0.25米足以过滤仪器停顿导致的平台又不会误删致密灰岩段的真实低GR平台。2.3 多井深度对齐用井壁微电阻率图像FMI或自然伽马GR特征点匹配单井处理后需跨井对比但不同次测井深度系统存在系统偏差可达0.3米。最可靠方法是基于地层标志层对齐。若无FMI图像可用GR曲线的峰值点如煤层、火山灰夹层作为锚点from scipy.signal import find_peaks # 提取GR峰值点要求高度100API宽度3个采样点 peaks, _ find_peaks(gr_clean, height100, distance3) # 获取峰值深度与幅度 anchor_depths depth[peaks] anchor_values gr_clean[peaks] # 假设参考井well_ref的锚点已知计算偏移量 # ref_anchors np.array([1250.3, 1387.6, 1522.1]) # 示例 # offset np.median(anchor_depths - ref_anchors) # 全局偏移 # las.df().index las.df().index - offset # 校正深度索引注意深度对齐必须在生成岩性曲线前完成。若先计算孔隙度再对齐会导致PHIE与AC曲线在深度上错位使Sw计算完全失效。3. 构建岩性敏感曲线从原始测井到Vsh、PHIE、SW的物理公式推导与代码实现岩性分类的根基不是原始曲线而是能表征地质意义的派生参数。station_95流程中Vsh泥质含量、PHIE有效孔隙度、SW含水饱和度是三大核心中间变量其计算必须严格遵循测井物理模型而非黑箱拟合。3.1 泥质含量Vsh用自然伽马GR与声波时差AC双参数交叉验证单一GR法在碳酸盐岩中失效因灰岩本身GR低故引入AC曲线泥岩声波时差高90μs/ft灰岩低50μs/ft。我们采用加权平均法融合两者def calculate_vsh(gr_clean, ac_clean, gr_min20, gr_max120, ac_min45, ac_max110): gr_min/gr_max: 纯砂岩/纯泥岩的GR值需根据区域标定 ac_min/ac_max: 纯砂岩/纯泥岩的AC值同上 # GR法线性刻度 vsh_gr np.clip((gr_clean - gr_min) / (gr_max - gr_min), 0, 1) # AC法同样线性刻度 vsh_ac np.clip((ac_clean - ac_min) / (ac_max - ac_min), 0, 1) # 加权融合AC在致密段更稳定GR在渗透层更敏感 weight_ac 0.7 * (1 - vsh_gr) 0.3 # Vsh越低AC权重越高 vsh_final weight_ac * vsh_ac (1 - weight_ac) * vsh_gr return vsh_final vsh calculate_vsh(gr_clean, las[AC])3.1.1 区域标定参数表需根据本区岩心数据修正岩性类型GR典型值(API)AC典型值(μs/ft)标定建议纯砂岩20–3545–55取试油成功层段的最低GR/AC纯泥岩90–13095–120取大段厚泥岩段的最高GR/AC灰岩5–1540–48若存在单独建立灰岩Vsh模型3.2 有效孔隙度PHIE密度-中子交会法消除岩性影响密度DEN与中子CNL曲线对孔隙流体响应相反但对骨架矿物响应一致。二者交会可消除岩性影响得到真实孔隙度def calculate_phie(den_clean, cnl_clean, den_ma2.65, den_fl1.0, # 砂岩骨架/流体密度 cnl_ma0.0, cnl_fl1.0): # 砂岩骨架/流体中子值 密度-中子交会法PHIE计算假设砂岩骨架 实际应用中需根据岩性切换ma参数如灰岩den_ma2.71, cnl_ma0.0 # 密度孔隙度 phie_den (den_ma - den_clean) / (den_ma - den_fl) # 中子孔隙度 phie_cnl (cnl_fl - cnl_clean) / (cnl_fl - cnl_ma) # 交会法取二者平均但当差异0.05时优先信密度因中子易受氯离子干扰 diff np.abs(phie_den - phie_cnl) phie_final np.where(diff 0.05, phie_den, (phie_den phie_cnl) / 2) return np.clip(phie_final, 0.01, 0.35) # 物理约束孔隙度0.01~0.35 phie calculate_phie(las[DEN], las[CNL])提示若井区存在大量白云岩需修改den_ma2.87, cnl_ma0.05否则PHIE系统性偏低。此参数必须由岩心分析数据标定不可套用邻井。3.3 含水饱和度SW阿尔奇公式与Simandoux公式的工程选择阿尔奇公式Archie仅适用于干净砂岩而Simandoux公式引入Vsh修正项更适合含泥质储层def calculate_sw(rt_clean, phie, vsh, a1.0, m2.0, n2.0, rw0.1, rsh5.0): Simandoux公式1/SW^n (a/PHIE^m) * (RT/RW) * (1-Vsh) (Vsh/RSH) * (1/PHIE^m) * (RT/RW) 参数说明 a,m,n: 阿尔奇参数本区岩心实验标定 rw: 地层水电阻率查SP曲线或实验室数据 rsh: 泥岩电阻率取Vsh0.7段的RT均值 term1 (a / (phie ** m)) * (rt_clean / rw) * (1 - vsh) term2 (vsh / rsh) * (1 / (phie ** m)) * (rt_clean / rw) sw_power 1 / (term1 term2) sw sw_power ** (1/n) return np.clip(sw, 0.15, 1.0) # 设定最小可动水饱和度0.15 sw calculate_sw(las[RT], phie, vsh)3.3.1 关键参数获取路径rw从自然电位SP曲线基线偏移量反算或取邻近水层实测值rsh在Vsh0.7的纯泥岩段取RT曲线的中位数a,m,n必须用本区岩心孔渗数据与测井RT数据回归禁用教科书默认值。4. 测井岩性分类从交会图规则到XGBoost模型的渐进式建模策略有了Vsh、PHIE、SW等物理参数岩性分类就从“经验猜”变为“证据链推理”。station_95流程强调先建规则、再调模型用交会图划定粗粒度岩性区间再用机器学习在边界模糊区精细区分。4.1 三参数交会图Vsh-PHIE-SW构成的岩性立方体将三维参数空间划分为8个卦限每个卦限对应一种主导岩性import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 构建岩性标签数组初始全为Unknown litho_labels np.full(len(vsh), Unknown, dtypeobject) # 定义规则示例基于某盆地标准 litho_labels[(vsh 0.15) (phie 0.18) (sw 0.4)] Clean_Sand litho_labels[(vsh 0.6) (phie 0.08)] Shale litho_labels[(vsh 0.2) (phie 0.06) (sw 0.8)] Tight_Carbonate litho_labels[(vsh 0.3) (phie 0.12)] Sandy_Shale # 可视化交会图 fig plt.figure(figsize(10, 8)) ax fig.add_subplot(111, projection3d) scatter ax.scatter(vsh, phie, sw, c[{Clean_Sand:0,Shale:1,Tight_Carbonate:2,Sandy_Shale:3}.get(l,4) for l in litho_labels], cmaptab10, s1) ax.set_xlabel(Vsh); ax.set_ylabel(PHIE); ax.set_zlabel(SW) plt.show()4.1.1 规则制定的地质依据Clean_Sand低泥质高孔隙低含水 → 典型主力产层Shale高泥质低孔隙 → 非储层但可作盖层Tight_Carbonate低泥质低孔隙高含水 → 白云岩致密层需酸压Sandy_Shale中高泥质中等孔隙 → 水平井靶窗优选区兼顾产能与可压性。4.2 XGBoost模型训练用岩心描述数据监督学习边界模糊区交会图无法区分“泥质粉砂岩”与“粉砂质泥岩”此时需引入岩心描述数据lithology description作为标签from xgboost import XGBClassifier from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report # 特征矩阵Vsh, PHIE, SW, GR, RT, AC6维 X np.column_stack([vsh, phie, sw, gr_clean, las[RT], las[AC]]) # 标签需提前准备岩心点深度对应的岩性编码如1Sand, 2Shale, 3Carbonate... # y_core load_core_labels(depth, core_log_file) # 此函数需自行实现 # 划分训练集仅用有岩心的深度点 X_train, X_test, y_train, y_test train_test_split( X[~np.isnan(y_core)], y_core[~np.isnan(y_core)], test_size0.3, random_state42 ) # 训练XGBoost关键参数控制过拟合 model XGBClassifier( n_estimators200, max_depth5, # 防止单棵树过深 learning_rate0.1, # 小步长提升稳定性 subsample0.8, # 行采样防过拟合 colsample_bytree0.8, # 列采样增加泛化 random_state42 ) model.fit(X_train, y_train) # 预测全井段 y_pred_full model.predict(X)4.2.1 模型评估与可信度阈值设定# 输出预测概率用于设定置信度阈值 y_proba model.predict_proba(X) max_proba np.max(y_proba, axis1) # 仅当最高概率0.85时采纳模型结果否则回退到交会图规则 final_litho np.where(max_proba 0.85, y_pred_full, litho_labels_encoded)注意XGBoost的predict_proba输出的是各岩性类别的概率估计非绝对确定性。实践中若最大概率0.7应标记为“Uncertain”触发人工复核而非强行分类。5. 岩性曲线可视化与解释验证用深度轨迹图误差热力图定位分类薄弱区最终输出的岩性曲线必须能被地质师快速验证。station_95流程强制要求双视图输出左侧为深度-岩性轨迹图直观右侧为分类误差热力图可追溯。5.1 绘制标准测井解释图Track Plotimport matplotlib.patches as mpatches def plot_litho_track(depth, gr_clean, rt_clean, phie, vsh, litho_labels, litho_colors{Clean_Sand:yellow,Shale:gray,Tight_Carbonate:blue,Sandy_Shale:brown}): fig, axes plt.subplots(1, 4, figsize(12, 10), shareyTrue) # GR Track axes[0].plot(gr_clean, depth, g, labelGR) axes[0].set_xlim(0, 150) axes[0].set_xlabel(GR (API)) axes[0].grid(True) # RT Track axes[1].semilogx(rt_clean, depth, r, labelRT) axes[1].set_xlim(0.2, 2000) axes[1].set_xlabel(RT (ohm.m)) axes[1].grid(True) # PHIE Vsh Track axes[2].plot(phie, depth, b, labelPHIE) axes[2].fill_betweenx(depth, 0, vsh, color0.7, alpha0.5, labelVsh) axes[2].set_xlim(0, 0.3) axes[2].set_xlabel(PHIE/Vsh) axes[2].legend() axes[2].grid(True) # Lithology Track柱状图 litho_numeric np.array([list(litho_colors.keys()).index(l) if l in litho_colors else 0 for l in litho_labels]) axes[3].fill_betweenx(depth, 0, 1, wherenp.isin(litho_labels, list(litho_colors.keys())), facecolor[litho_colors.get(l, white) for l in litho_labels], alpha0.8) axes[3].set_xlim(0, 1) axes[3].set_xlabel(Lithology) axes[3].set_yticks([]) # 隐藏y轴刻度保持整洁 # 添加图例 legend_elements [mpatches.Patch(facecolorcolor, labellith) for lith, color in litho_colors.items()] axes[3].legend(handleslegend_elements, loccenter left, bbox_to_anchor(1, 0.5)) plt.tight_layout() plt.show() plot_litho_track(depth, gr_clean, las[RT], phie, vsh, litho_labels)5.2 生成分类误差热力图定位模型失效的深度段误差热力图的核心是对比模型预测与交会图规则的分歧点这些点往往是地质复杂区如断层破碎带、流体突变层# 计算规则法与模型法的差异0一致1不一致 rule_vs_model (litho_labels ! [list(litho_colors.keys())[i] for i in y_pred_full]) # 转换为滑动窗口误差率每10米统计一次不一致比例 window_size_m 10 window_points int(window_size_m / np.mean(np.diff(depth))) error_rate [] for i in range(0, len(depth), window_points): window rule_vs_model[i:iwindow_points] error_rate.append(np.mean(window) if len(window) 0 else 0) # 绘制热力图 plt.figure(figsize(10, 6)) plt.imshow([error_rate], cmapRdYlBu_r, aspectauto, extent[0, len(error_rate), depth.min(), depth.max()]) plt.colorbar(labelClassification Disagreement Rate) plt.xlabel(Window Index (10m each)) plt.ylabel(Depth (m)) plt.title(Model vs Rule Disagreement Heatmap) plt.show()5.2.1 误差热力图的地质解读指南红色热点误差率0.6立即检查该深度段的岩心照片或FMI图像常对应断层角砾岩、沥青充填缝洞或钻井液侵入带连续黄色条带误差率0.3~0.6提示该层段岩性过渡渐变需调整交会图边界或增加训练样本全蓝区域误差率0.1模型与规则高度一致可放心用于批量处理。关键技巧在误差热力图中叠加试油结论如“1250–1255m日产油32t”若高产层恰好位于红色热点区说明此处岩性分类虽难但恰恰是优质储层——这正是station_95流程的价值不追求全局准确率而聚焦于识别“高价值不确定性”。本文还有配套的精品资源点击获取