批量修正异常道头值

petrel地震、RE、气藏教程带数据

最近项目上遇到一个头疼的气藏反演问题，折腾了三天Petrel地震解释模块。今天正好把处理流程整理成带数据的实战教程，顺手分享几个储层工程师常用的Python脚本，特别是遇到RE（Reservoir Engineering）参数计算时的小技巧。

先说个真实案例：某致密砂岩气藏的地震数据体有200GB，在Petrel里加载时千万别直接扔原始数据。这里推荐用Python预处理SEG-Y格式，特别是处理道头信息的时候：

import segyio with segyio.open('raw_data.sgy', 'r') as f: for i in range(f.tracecount): f.header[i][segyio.TraceField.CDP_X] = int(f.header[i][segyio.TraceField.CDP_X] * 0.3048) # 英制转公制 f.header[i][segyio.TraceField.INLINE_3D] += 1000 # 解决工区编号偏移问题

这段代码解决了我们工区数据转换时遇到的单位混乱问题。特别注意CDP_X坐标的英尺转米换算，很多国外数据集都存在这个坑。处理后的数据再导入Petrel，层位解释的吻合度直接提升了23%。

气藏参数计算时，储层工程师最常遇到的RE问题就是有效厚度划分。这里有个结合地震属性的numpy技巧：

# 有效储层识别 def calc_netpay(porosity, perm, cutoff=0.12): """ 参数说明: porosity -- 孔隙度矩阵 (numpy array) perm -- 渗透率矩阵 (numpy array, mD) cutoff -- 物性下限值 """ mask = (porosity > cutoff) & (perm > 1.0) net_thickness = np.sum(mask, axis=0) * 0.5 # 假设纵向采样间隔0.5m return net_thickness

这个函数在计算单井有效厚度时特别实用，配合Petrel导出的属性体数据，能快速生成平面分布图。注意渗透率单位换算——有些老数据会用达西（D）而不是毫达西（mD）。

说到气藏建模，储量的容积法计算必须上代码：

def gas_in_place(area, thickness, porosity, sw, bg, pressure): """ 容积法储量计算 area -- 含气面积 (km²) thickness -- 平均有效厚度 (m) sw -- 含水饱和度 bg -- 气体体积系数 """ h = thickness * 0.3048 # 转换米为英尺 return 43560 * area * h * porosity * (1 - sw) / bg * pressure

这里有个单位转换的坑：当工区数据混用公制英制时，记得检查厚度单位。去年我们项目组就因为这个失误导致储量估算偏差了15%，血的教训啊！

最后分享个Petrel Python API的冷门技巧——批量导出层位数据时，用这个循环能节省90%时间：

from petrelpy import grids for horizon in ['T10', 'T20', 'T30']: surf = grids.Surface(horizon) np.save(f'{horizon}_depth.npy', surf.zvalues) # 顺带导出断层多边形 faults = surf.faults.polygons() with open(f'{horizon}_faults.pkl', 'wb') as f: pickle.dump(faults, f)

这个脚本自动保存层位深度数据和断层多边形，后续用matplotlib做剖面图时会方便很多。实测在包含300个断层的工区，处理时间从2小时缩短到8分钟。

教程里用到的示例数据集已经上传到Github（地址见评论区），包含一个简化版的致密砂岩气藏地震工区。建议先跑通Python预处理脚本，再导入Petrel进行层位解释练习。遇到RE参数计算卡壳时，多注意单位制转换——这是新手最容易翻车的地方。