测井岩性分类：物理建模与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", encoding='utf-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_size=10, spike_sigma=3, flat_length=5): cleaned = curve_data.copy() n = len(curve_data) for i in range(window_size, n - window_size): window = curve_data[i-window_size:i+window_size+1] 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:i+flat_length] == curve_data[i]) if is_flat: # 平台段用前后窗口均值插值 left_mean = np.nanmean(curve_data[max(0,i-5):i]) right_mean = np.nanmean(curve_data[i+1:min(n,i+6)]) cleaned[i] = (left_mean + right_mean) / 2 else: # 尖峰用线性插值 cleaned[i] = np.interp(depth[i], [depth[i-1], depth[i+1]], [curve_data[i-1], curve_data[i+1]]) return cleaned # 应用清洗（以GR曲线为例） gr_raw = las['GR'] gr_clean = clean_curve(gr_raw, depth)

2.2.1 参数说明与调优逻辑

window_size=10：对应约0.5米窗口（按0.05m采样率），太小易受噪声干扰，太大则无法捕捉局部异常；
spike_sigma=3：符合高斯分布3σ原则，对碳酸盐岩中天然放射性异常（如钾长石富集）保留容忍；
flat_length=5：对应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, height=100, distance=3) # 获取峰值深度与幅度 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_min=20, gr_max=120, ac_min=45, ac_max=110): """ 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–35	45–55	取试油成功层段的最低GR/AC
纯泥岩	90–130	95–120	取大段厚泥岩段的最高GR/AC
灰岩	5–15	40–48	若存在，单独建立灰岩Vsh模型

3.2 有效孔隙度PHIE：密度-中子交会法消除岩性影响

密度（DEN）与中子（CNL）曲线对孔隙流体响应相反，但对骨架矿物响应一致。二者交会可消除岩性影响，得到真实孔隙度：

def calculate_phie(den_clean, cnl_clean, den_ma=2.65, den_fl=1.0, # 砂岩骨架/流体密度 cnl_ma=0.0, cnl_fl=1.0): # 砂岩骨架/流体中子值 """ 密度-中子交会法PHIE计算（假设砂岩骨架） 实际应用中需根据岩性切换ma参数（如灰岩：den_ma=2.71, cnl_ma=0.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_ma=2.87, cnl_ma=0.05，否则PHIE系统性偏低。此参数必须由岩心分析数据标定，不可套用邻井。

3.3 含水饱和度SW：阿尔奇公式与Simandoux公式的工程选择

阿尔奇公式（Archie）仅适用于干净砂岩，而Simandoux公式引入Vsh修正项，更适合含泥质储层：

def calculate_sw(rt_clean, phie, vsh, a=1.0, m=2.0, n=2.0, rw=0.1, rsh=5.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: 泥岩电阻率（取Vsh>0.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：在Vsh>0.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', dtype=object) # 定义规则（示例：基于某盆地标准） 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, projection='3d') 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], cmap='tab10', s=1) 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, AC（6维） X = np.column_stack([vsh, phie, sw, gr_clean, las['RT'], las['AC']]) # 标签：需提前准备岩心点深度对应的岩性编码（如1=Sand, 2=Shale, 3=Carbonate...） # 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_size=0.3, random_state=42 ) # 训练XGBoost（关键参数：控制过拟合） model = XGBClassifier( n_estimators=200, max_depth=5, # 防止单棵树过深 learning_rate=0.1, # 小步长提升稳定性 subsample=0.8, # 行采样防过拟合 colsample_bytree=0.8, # 列采样增加泛化 random_state=42 ) 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, axis=1) # 仅当最高概率>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 Plot）

import 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), sharey=True) # GR Track axes[0].plot(gr_clean, depth, 'g', label='GR') 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', label='RT') 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', label='PHIE') axes[2].fill_betweenx(depth, 0, vsh, color='0.7', alpha=0.5, label='Vsh') 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, where=np.isin(litho_labels, list(litho_colors.keys())), facecolor=[litho_colors.get(l, 'white') for l in litho_labels], alpha=0.8) axes[3].set_xlim(0, 1) axes[3].set_xlabel('Lithology') axes[3].set_yticks([]) # 隐藏y轴刻度，保持整洁 # 添加图例 legend_elements = [mpatches.Patch(facecolor=color, label=lith) for lith, color in litho_colors.items()] axes[3].legend(handles=legend_elements, loc='center 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:i+window_points] error_rate.append(np.mean(window) if len(window) > 0 else 0) # 绘制热力图 plt.figure(figsize=(10, 6)) plt.imshow([error_rate], cmap='RdYlBu_r', aspect='auto', extent=[0, len(error_rate), depth.min(), depth.max()]) plt.colorbar(label='Classification Disagreement Rate') plt.xlabel('Window Index (10m each)') plt.ylabel('Depth (m)') plt.title('Model vs Rule Disagreement Heatmap') plt.show()