简介这份资源面向地质勘探、地球物理与油藏工程方向的研究人员和学生聚焦岩石物理中弹性参数的理论计算与流体饱和影响分析。包内提供一组MATLAB函数围绕Gassman公式展开可用于模拟不同流体饱和状态下岩石的弹性响应并涉及Reuss下限、Voigt上限以及二者结合的估算方法帮助理解多相介质平均弹性模量的边界与修正思路。资源共6个文件全部为m脚本压缩包约3KB涵盖Gassman公式实现、上下限计算、流体敏感性参数测量等模块结构紧凑、便于直接调用与二次开发。已有366人学习下载适合具备一定岩石物理与MATLAB基础、希望将理论模型落地到地层评价与石油勘探实践中的读者参考。1. 岩石物理性质从测井曲线到油藏模型的“翻译官”钻井取心成本高、周期长一口井取心段往往只有几十米而测井曲线却是全井段连续采集的。岩石物理要解决的核心矛盾就在这里如何用连续的测井响应反推出离散岩心才能标定的储层参数。孔隙度、渗透率、饱和度、泥质含量这四个参数直接决定了一个区块有没有商业产能。做岩石物理性质评价的人日常就是在岩心刻度测井、测井标定地震、地震外推储层这条链条上反复迭代。如果你手上有常规测井曲线和少量岩心分析数据想建立孔渗饱解释模型或者想判断某个层位是气层还是水层这篇笔记会从数据准备、模型标定到实际验证把中间容易翻车的地方都摊开讲。岩石物理性质不是一门纯理论学科它更像一门手艺——参数选对了模型能外推几公里选错了连井间都推不过去。2. 岩心刻度测井孔隙度与渗透率的标定流程2.1 为什么孔隙度模型必须分岩性标定孔隙度测井的响应本质上是岩石骨架、孔隙流体和泥质三者的加权平均。砂岩骨架密度约 2.65 g/cm³灰岩约 2.71 g/cm³白云岩约 2.87 g/cm³如果用一个统一的骨架密度值去计算孔隙度在岩性变化的井段误差可以超过 5 个孔隙度单位。我见过一个区块上部砂岩段用 2.65 算孔隙度跟岩心吻合得很好到了下部含灰质段计算孔隙度比岩心高了 8 个单位直接导致储量计算虚高。常见做法是先用自然伽马或中子-密度交会图做岩性划分然后对每个岩性段分别标定骨架参数。具体步骤是整理岩心分析孔隙度数据剔除裂缝发育和破碎样品将岩心深度与测井深度对齐这一步的深度校正没做好后面所有标定都是白费按岩性分组做交会图回归出孔隙度计算公式。import numpy as np import pandas as pd from scipy import stats # 读取岩心分析数据和测井曲线数据 core_data pd.read_csv(core_porosity.csv) # 列: depth, phi_core, lithology log_data pd.read_csv(log_data.csv) # 列: depth, rhob, dt, gr, nphi # 深度对齐将岩心深度匹配到最近的测井采样点 def align_depth(core_df, log_df, tolerance0.15): 岩心深度与测井深度对齐tolerance为最大允许深度差(米) merged pd.merge_asof( core_df.sort_values(depth), log_df.sort_values(depth), ondepth, tolerancetolerance, directionnearest ) return merged.dropna() aligned align_depth(core_data, log_data) # 分岩性标定密度孔隙度公式: phi (rho_ma - rhob) / (rho_ma - rho_f) def calibrate_density_porosity(df, rho_f1.0): 对每种岩性回归骨架密度rho_ma results {} for lith in df[lithology].unique(): subset df[df[lithology] lith] # 由岩心孔隙度反推骨架密度: rho_ma rhob phi*(rho_ma - rho_f) # 线性回归: rhob rho_ma - phi*(rho_ma - rho_f) slope, intercept, r_value, p_value, std_err stats.linregress( subset[phi_core], subset[rhob] ) rho_ma_est intercept # 当phi0时的密度即为骨架密度 r_squared r_value ** 2 results[lith] { rho_ma: round(rho_ma_est, 3), r_squared: round(r_squared, 4), n_samples: len(subset) } return results calibration_results calibrate_density_porosity(aligned) for lith, params in calibration_results.items(): print(f岩性: {lith}, 骨架密度: {params[rho_ma]} g/cm³, fR²: {params[r_squared]}, 样品数: {params[n_samples]})这段代码的关键在于用线性回归反推骨架密度当岩心孔隙度趋近于零时对应的密度值就是骨架密度。tolerance参数控制深度对齐的容差一般取 0.10.2 米太大容易把不同层的数据混在一起太小则匹配不上几个点。r_squared低于 0.7 的岩性段需要检查是否混入了非均质样品。2.2 渗透率标定的三种思路与适用边界渗透率是岩石物理性质里最难准确预测的参数没有之一。孔隙度与渗透率的相关性在碎屑岩中通常只有 0.60.8碳酸盐岩中更差。常见做法有三类孔隙度-渗透率回归、Timur 公式、以及基于流动单元的分类标定。孔隙度-渗透率回归最简单在对数坐标下做线性拟合log(K) a * phi b。这种做法在砂岩中能用但误差通常在一个数量级以内。Timur 公式引入了束缚水饱和度作为第二变量K 0.136 * phi^4.4 / Swi^2适合有核磁共振测井的井。流动单元分类标定是最靠谱但也是最费事的做法先按孔喉半径或 FZI 把储层分成若干类每类单独回归。# 流动单元分类标定渗透率 def calculate_fzi(phi, k): 计算流动单元指标FZI # RQI 0.0314 * sqrt(k/phi), 单位转换后 rqi 0.0314 * np.sqrt(k / phi) fzi rqi / (phi / (1 - phi)) return fzi def classify_flow_units(fzi_values, n_units3): 按FZI值聚类划分流动单元 from sklearn.cluster import KMeans fzi_log np.log10(fzi_values.reshape(-1, 1)) kmeans KMeans(n_clustersn_units, random_state42, n_init10) labels kmeans.fit_predict(fzi_log) # 按FZI均值排序确保单元编号有意义 centers kmeans.cluster_centers_.flatten() order np.argsort(centers) label_map {old: new for new, old in enumerate(order)} return np.array([label_map[l] for l in labels]) # 对每个流动单元单独回归孔渗关系 aligned[FZI] calculate_fzi(aligned[phi_core], aligned[k_core]) aligned[flow_unit] classify_flow_units(aligned[FZI].values) for unit in sorted(aligned[flow_unit].unique()): subset aligned[aligned[flow_unit] unit] log_k np.log10(subset[k_core]) slope, intercept, r, p, se stats.linregress(subset[phi_core], log_k) print(f流动单元 {unit}: log(K) {slope:.2f}*phi {intercept:.2f}, fR² {r**2:.3f}, 样品数 {len(subset)})n_units的选择需要结合孔喉半径分布和压汞曲线确定一般 24 类。分类过少则不同流动能力的储层混在一起回归精度上不去分类过多则每类样品数不足回归系数不稳定。我一般会先看 FZI 的对数分布直方图如果呈现明显的双峰或三峰就按峰的数量定类别数。3. 饱和度模型从阿尔奇公式到低阻油层修正3.1 阿尔奇公式的参数确定与适用条件阿尔奇公式 Sw [(a * Rw) / (phi^m * Rt)]^(1/n) 是饱和度解释的基石但它的适用条件很苛刻岩石必须满足孔隙结构均匀、不含泥质、不导电矿物少、润湿性为水湿。实际储层中完全满足这些条件的很少所以参数 a、m、n 的确定就成了关键。a 和 m 通过岩电实验确定对同一岩性组的岩心样品测量不同矿化度地层水下的地层因素 F Ro/Rw 与孔隙度 phi 的关系在双对数坐标下回归得到 F a / phi^m。n 通过电阻增大率 I Rt/Ro 与含水饱和度 Sw 的关系确定I 1 / Sw^n。# 岩电实验数据回归a, m, n def fit_archie_params(porosity, formation_factor, sw, resistivity_index): 回归阿尔奇参数a, m, n # F a / phi^m log(F) log(a) - m*log(phi) log_phi np.log10(porosity) log_F np.log10(formation_factor) m_slope, log_a, r_f, p_f, se_f stats.linregress(log_phi, log_F) m -m_slope a 10 ** log_a # I 1 / Sw^n log(I) -n*log(Sw) log_sw np.log10(sw) log_I np.log10(resistivity_index) n_slope, intercept_i, r_i, p_i, se_i stats.linregress(log_sw, log_I) n -n_slope return { a: round(a, 3), m: round(m, 3), n: round(n, 3), R2_F: round(r_f**2, 4), R2_I: round(r_i**2, 4) } # 示例数据 phi_samples np.array([0.15, 0.18, 0.22, 0.25, 0.28, 0.30]) F_samples np.array([45.2, 32.1, 21.5, 16.8, 13.2, 11.5]) sw_samples np.array([1.0, 0.8, 0.6, 0.5, 0.4, 0.3]) I_samples np.array([1.0, 1.56, 2.78, 4.0, 6.25, 11.1]) params fit_archie_params(phi_samples, F_samples, sw_samples, I_samples) print(fa{params[a]}, m{params[m]}, n{params[n]}) print(fF回归R²{params[R2_F]}, I回归R²{params[R2_I]})R2_F低于 0.85 说明岩性组内孔隙结构差异大需要进一步细分岩性。m 值一般在 1.82.2 之间n 值在 1.82.5 之间偏离这个范围要检查实验数据质量。3.2 低阻油层与泥质校正的实战处理低阻油层是饱和度解释中最容易翻车的地方。现象是油层电阻率跟水层差不多用阿尔奇公式算出来含水饱和度 70% 以上但试油却产纯油。原因通常有三类泥质附加导电、微孔隙发育导致束缚水饱和度高、以及导电矿物存在。泥质校正的常见做法是采用 Waxman-Smits 模型或双水模型。Waxman-Smits 模型的核心思想是黏土阳离子交换容量 CEC 提供了额外的导电通路使得含油饱和度下的电阻率降低。实际应用中CEC 可以通过自然伽马或中子-密度交会估算也可以直接测量。# Waxman-Smits模型计算含水饱和度 def waxman_smits_sw(Rt, Rw, phi, m, n, Qv, B4.0): Waxman-Smits模型 Rt: 地层电阻率, Rw: 地层水电阻率, phi: 孔隙度 m: 胶结指数, n: 饱和度指数, Qv: 阳离子交换容量(meq/mL) B: 阳离子等效电导率, 默认4.0 (S/m)/(meq/mL) # 迭代求解Sw Sw 0.5 # 初始猜测 for _ in range(100): # 考虑黏土导电后的等效地层水电阻率 Rw_eff 1.0 / (1.0/Rw B * Qv / Sw) Sw_new ((1.0 / Rw_eff) * (1.0 / (phi ** m)) * (1.0 / Rt)) ** (1.0 / n) if abs(Sw_new - Sw) 1e-6: break Sw Sw_new return Sw # 对比阿尔奇与Waxman-Smits结果 Rt, Rw, phi, m, n 5.0, 0.1, 0.25, 2.0, 2.0 Qv 0.15 # 阳离子交换容量 Sw_archie ((1.0/Rw) * (1.0/(phi**m)) * (1.0/Rt)) ** (1.0/n) Sw_ws waxman_smits_sw(Rt, Rw, phi, m, n, Qv) print(f阿尔奇含水饱和度: {Sw_archie:.3f}) print(fWaxman-Smits含水饱和度: {Sw_ws:.3f})Qv的取值对结果影响很大每 0.05 meq/mL 的变化可以导致含水饱和度变化 510 个百分点。没有 CEC 测量数据时可以用自然伽马相对值估算Qv ≈ 0.1 * (GR - GR_min) / (GR_max - GR_min) * 孔隙度。这个估算关系是经验性的不同区块需要重新标定。4. 岩石物理性质评价的避坑与排查4.1 深度校正没做好后面全白搭现象岩心孔隙度与测井计算孔隙度交会图相关性极差R² 不到 0.3但单看数据质量都没问题。原因岩心深度和测井深度来自不同的测量系统钻具深度、电缆深度、岩心归位深度之间存在系统偏差。这个偏差在直井中可能只有 0.5 米但在斜井或水平井中可以超过 2 米。解决先用自然伽马曲线做深度校正基准因为岩心的自然伽马测量值可以和测井自然伽马直接对比。具体做法是把岩心自然伽马数据与测井自然伽马曲线做互相关找到最佳匹配深度偏移量。没有岩心伽马数据时用岩心描述中的泥质夹层或钙质夹层深度与测井曲线上的对应特征深度对齐。深度校正的精度要控制在 0.2 米以内否则薄层的标定会完全失真。4.2 用全井段数据回归忽略了岩性变化现象孔隙度模型在标定井段误差很小但应用到邻井后误差突然变大有的层位偏差超过 6 个孔隙度单位。原因标定井的岩心段恰好集中在某一岩性中回归出的骨架参数只适用于该岩性。邻井钻遇了不同的岩性组合参数就不适用了。解决回归前先按自然伽马或 Pe 曲线做岩性分组确保每个岩性组至少有 1520 个岩心样品。样品数不足时可以从区域上收集同层位的岩心数据补充但要注意不同区块的成岩作用差异。我一般会在交会图上用不同颜色标注岩性如果看到数据点明显分成几团就说明必须分组标定。4.3 渗透率回归直接外推到未取心段现象取心段渗透率回归 R² 达到 0.85但外推到未取心段后试井渗透率比测井解释渗透率低了两个数量级。原因岩心分析的渗透率是基质渗透率而试井渗透率反映的是包含裂缝和层理的综合渗透率。在裂缝发育的储层中两者可以差几个数量级。另外岩心分析通常只测量垂向渗透率而水平渗透率可能高得多。解决渗透率解释结果必须用试井或生产动态数据验证。如果试井渗透率远大于岩心渗透率说明裂缝贡献不可忽略需要用裂缝孔隙度或裂缝密度曲线做二次修正。水平井中还要考虑渗透率各向异性用水平渗透率与垂向渗透率的比值做校正。4.4 饱和度参数直接借用邻区经验值现象新区的第一口探井用邻区经验参数 a1, m2, n2 计算含水饱和度结果油层算出来 Sw 高达 65%差点把油层判成水层。原因a、m、n 受孔隙结构、润湿性、地层水矿化度影响不同区块之间可以差很多。m 值从 1.8 到 2.5 的变化在孔隙度 20% 时会导致含水饱和度变化 15 个百分点以上。解决新区第一口井必须做岩电实验没有岩心时至少要做地层水分析确定 Rw。如果连 Rw 都没有可以用自然电位曲线或水分析资料估算。在参数不确定的情况下用敏感性分析给出含水饱和度的范围而不是一个确定值。我一般会做 m 值 ±0.2、n 值 ±0.3 的敏感性分析看看含油饱和度是否跨越了油水界限。4.5 忽略泥质含量对孔隙度的影响现象含泥质砂岩段密度孔隙度与中子孔隙度差异很大密度孔隙度 25%中子孔隙度 35%两者交会图上的点偏离砂岩线。原因泥质的中子响应与石英完全不同蒙脱石的中子孔隙度可以到 40% 以上而石英骨架的中子孔隙度接近零。泥质含量每增加 10%中子孔隙度就会偏高 35 个单位。解决先用自然伽马或中子-密度交会计算泥质含量然后对孔隙度做泥质校正。常见做法是phi_corrected phi_total - Vsh * phi_sh其中 phi_sh 是纯泥岩的孔隙度一般取 0.250.35。校正后的孔隙度再用岩心数据验证。如果泥质含量超过 30%密度-中子交会法计算孔隙度的误差会显著增大此时核磁共振测井是最可靠的选择。5. 用生产数据反向验证岩石物理模型岩石物理模型建完之后最容易被忽略的一步是验证。我见过太多人把交会图 R² 做到 0.9 就收工结果模型应用到油藏数值模拟里历史拟合怎么都拟合不上。问题出在岩石物理模型只用了静态数据标定没有用动态数据检验。验证方法分三个层次。第一层是测井解释成果与试油数据对比把解释的含水饱和度与试油产水率做交会如果 Sw 低于 40% 的层试油产纯水说明饱和度模型有问题。第二层是岩心渗透率与试井渗透率对比在双对数坐标下如果数据点偏离 1:1 线超过一个数量级需要检查裂缝或各向异性影响。第三层是数值模拟历史拟合把岩石物理模型输出的孔渗饱参数输入模拟器看压力和生产历史是否吻合。# 用试油数据验证饱和度模型 def validate_sw_model(sw_log, water_cut, sw_cutoff0.4): sw_log: 测井解释含水饱和度 water_cut: 试油产水率 sw_cutoff: 含水饱和度界限 import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(8, 6)) ax.scatter(sw_log, water_cut, alpha0.6, edgecolorsk, linewidth0.5) ax.axvline(xsw_cutoff, colorr, linestyle--, labelfSw界限{sw_cutoff}) ax.axhline(y0.5, colorgray, linestyle:, label产水率50%) ax.set_xlabel(测井解释含水饱和度 Sw) ax.set_ylabel(试油产水率) ax.set_title(饱和度模型验证交会图) ax.legend() ax.grid(True, alpha0.3) plt.tight_layout() plt.savefig(sw_validation.png, dpi150) # 统计误判率 false_oil np.sum((sw_log sw_cutoff) (water_cut 0.5)) false_water np.sum((sw_log sw_cutoff) (water_cut 0.2)) total len(sw_log) print(f油层误判为水层: {false_oil}/{total} ({false_oil/total*100:.1f}%)) print(f水层误判为油层: {false_water}/{total} ({false_water/total*100:.1f}%)) return false_oil, false_water # 示例数据 sw_interpreted np.array([0.25, 0.35, 0.45, 0.55, 0.65, 0.30, 0.50, 0.70, 0.20, 0.40]) water_cut_measured np.array([0.05, 0.10, 0.30, 0.70, 0.90, 0.08, 0.55, 0.95, 0.02, 0.25]) validate_sw_model(sw_interpreted, water_cut_measured)误判率超过 15% 时模型需要重新标定。重点检查低阻油层和高阻水层这两个极端情况。低阻油层误判为水层通常是泥质校正不够高阻水层误判为油层通常是忽略了导电矿物或裂缝的影响。我自己的习惯是每建完一个岩石物理模型先不急着外推而是找三口有试油数据的井做盲测。盲测通过率超过 80% 才敢往新区推广。这个习惯让我避免了好几次重大误判有一次盲测发现某区块的 m 值应该取 2.3 而不是 2.0如果直接外推整个区块的含油饱和度会系统性偏高 12 个百分点。希望帮到你。本文还有配套的精品资源点击获取