0 引言

地质建模技术兴起于上世纪八十年代，是油田中后期开展油藏描述的有效技术手段，也是从油田勘探、评价延伸至油田开发阶段的桥梁。随着计算机软、硬件技术不断革新，地质建模功能越来越强大，融合了层序地层学、构造地质学、储层沉积学、石油地质学、油藏工程学、地质统计学等多个领域，成为集井、震、藏、地质等多学科领域于一体的综合技术平台^[1]，在油藏描述领域扮演着“航母”角色。地质建模主要分为构造建模和属性建模，其中，属性模型建立在构造模型基础上，可分为离散型模型和连续型模型。离散型模型以沉积相模型和裂缝模型为主；连续型模型包括孔隙度模型、渗透率模型、饱和度模型、净毛比模型等，主要为后期油藏数值模拟提供静态数据模型。由此可见，地质建模是油田开发中后期不可缺少的技术手段之一^[2]。

目前，中国东部油田多数已进入高含水或特高含水时期。以渤海油田为首的海上油田亦逐渐进入油田开发中后期^[3]，随着油田开发程度不断提高，对储层的认识亦不断深入。以渤海湾盆地SZ36-1油田为例，油田位于辽东湾下辽河坳陷辽西低凸起中段，自下而上发育孔店组、沙河街组、东营组、馆陶组、明化镇组和平原组等地层。主要含油气目的层为古近系东营组下段，开发层系为东营组东二下段Ⅰ、Ⅱ油组，又细分为14个小层，属于大型河流三角洲沉积复合体。储层岩性以中—细粒岩屑长石砂岩、长石石英细砂岩和细粉砂岩为主，具有胶结疏松、变异系数高、非均质性强的特征，属于中高孔—中高渗储层。受复杂地质构造作用，储层连续性和连通性相对较差，油藏呈现复杂的油水系统，严重制约了油田剩余油挖潜，同时对地质建模提出了更高的技术需求^[4]。为此，一方面需要更加系统地融合地球物理、石油地质、油藏工程等领域信息作为建模原始约束条件，另一方面需借助多学科技术手段克服模型多解性和不确定性，降低错误概率，建立更高精度、更为可信的油藏地质模型，才能为油田开发提供有效的指导。

众所周知，地质建模技术多以地质统计学为核心理论^[5]，前人在储层构型建模、相控建模、多手段一体化建模等方面进行了深入探讨，并取得了颇丰的认识^[6-8]。但是，在基于多学科地质建模方面的认识仍存在一定的局限性^[9]，具体体现在三个方面。①基于单一变差函数建模精度低。由于变差函数拟合的多解性，导致井间预测模拟结果存在很大的不确定性。该方法对井距大小敏感性较强，在密井区具有较高的精确度，在大井距甚至无井区域难以分析储层属性参数变差函数结果，造成属性模型的精度、可靠性大大降低。②相控建模模式化严重。受沉积微相展布规律和已有沉积模式认识的约束，相控物性模拟结果往往与沉积微相边界趋势一致。该方法融入了地质思想，同时还受研究人员主观因素的影响，模式化间距严重。③地震约束地质建模不够规范^[10-13]。目前所报道的地震约束建模技术多样化，特别是地震约束建模的原理和流程不够明确。

针对以上问题，本文以渤海湾盆地SZ36-1油田为例，试图建立一套基于地震反演和沉积微相约束地质建模的方法。在进行相控建模之前，需分微相拟合试验变差函数，确定不同方向的变程，加入相模型作为约束条件，以地质规律控制属性模拟，表征和预测储层各沉积微相内部物性的变化特征，重点阐述地震约束地质建模原理和转化过程，以期有效指导油田下一步勘探部署。

1 地震约束建模原理

通常，地震数据约束建模采用协克里金法^[14]。笔者借鉴了前人研究成果，基于地震数据及测井数据的空间相关性，利用协克里金法对储层物性进行较为准确的描述，为了简化说明，仅使用两个变量{Z₁(x)，Z₂(x)}构成协同区域化变量，其二阶平稳假设如下：

(1) 每一个变量的期望存在且为常数，即

$ E\left[Z_{k}(x)\right]=m_{k} \quad k=1,2 $

(1)

式中E表示期望；

(2) 每一个变量的空间协方差存在且为空间滞后h的函数，与绝对空间位置无关，即

$ \operatorname{cov}\left[Z_{k}(x+h), Z_{k}(x)\right]=C_{k k}(h) \quad k=1,2 $

(2)

式中cov表示协方差；

(3) 两个变量的交叉协方差函数存在且为空间滞后h的函数，与绝对空间位置无关，即

$ \operatorname{cov}\left[Z_{k}(x), Z_{k^{\prime}}(x+h)\right]=C_{k k^{\prime}}(h) \quad k, k^{\prime}=1,2 $

(3)

交叉协方差中k和k′的顺序不能颠倒。

内蕴假设中使用变量在一定空间滞后上的增量的期望、变异函数和交叉变异函数。

二阶平稳性假设下，单一区域化变量具有关系γ(h)=C(0)-C(h)。相应地，交叉变异函数和交叉协方差函数具有下列转换关系

$ \gamma_{k k^{\prime}}(h)=C_{k k^{\prime}}(0)-\frac{1}{2}\left[C_{k k^{\prime}}(h)+C_{k^{\prime} k}(h)\right] $

(4)

假设测井数据为区域化变量Z₂(x)，观测值个数为N₂；地震数据为区域变化量Z₁(x)，为辅助变量，观测值个数为N₁，且N₂ < N₁。综合利用Z₂(x)和Z₁(x)的观测值对x₀位置的Z₂(x)进行估计，协同克里金估计公式为

$ Z_{2}^{*}\left(x_{0}\right)=\sum\limits_{i=1}^{N_{1}} \lambda_{1 i} Z_{1}\left(x_{1 i}\right)+\sum\limits_{j=1}^{N_{2}} \lambda_{2 j} Z_{2}\left(x_{2 j}\right) $

(5)

式中：Z^*(x₀)是在待估位置x₀的估计值；Z(x_i)是区域化变量Z(x)在位置x_i的观测值；λ_i是分配给Z(x_i)的权重。

协同克里金估计方法的无偏估计数学表达式为

$ \begin{array}{l} E\left[Z_{2}^{*}\left(x_{0}\right)-Z_{2}\left(x_{0}\right)\right] \\ \quad=E\left[\sum\limits_{i=1}^{N_{1}} \lambda_{1} Z_{1}\left(x_{1 i}\right)+\sum\limits_{j=1}^{N_{2}} \lambda_{2 j} Z_{2}\left(x_{2 j}\right)-Z_{2}\left(x_{0}\right)\right] \\ \quad=m_{1} \sum\limits_{i=1}^{N_{1}} \lambda_{1 i}+m_{2}\left(\sum\limits_{j=1}^{N_{2}} \lambda_{2 j}-1\right)=0 \end{array} $

(6)

式中m₁、m₂分别是变量Z₁(x)和Z₂(x)的期望。

为了保证无偏估计对于任意m₁和m₂都成立，上式第一个变量权重之和应为0，即$\sum\limits_{i = 1}^{{N_1}} {{\lambda _{1i}}} = 0$。第二个变量的权重之和为1，即$\sum\limits_{j = 1}^{{N_2}} {{\lambda _{2j}}} = 1$。

协同克里金估计方法的最优估计要求估计方差最小，var[Z^*(x₀)-Z(x₁)]=min。根据估计公式并结合无偏估计条件表达式，进一步化简估计方差表达式，获得

$ \begin{aligned} \operatorname{var}[&\left.Z_{2}^{*}\left(x_{0}\right)-Z_{2}\left(x_{0}\right)\right] \\ &=E\left[Z_{2}^{*}\left(x_{0}\right)-Z_{2}\left(x_{0}\right)\right]-E\left[Z_{2}^{*}\left(x_{0}\right)-Z_{2}\left(x_{0}\right)\right]^{2} \\ &=E\left[Z_{2}^{*}\left(x_{0}\right)-Z_{2}\left(x_{0}\right)\right]^{2} \\ &=E\left[\sum\limits_{i=1}^{N_{1}} \lambda_{1 i} Z_{1}\left(x_{1 i}\right)+\sum\limits_{j=1}^{N_{2}} \lambda_{2 j} Z_{2}\left(x_{2 j}\right)-Z_{2}\left(x_{0}\right)\right]^{2} \end{aligned} $

(7)

引入两个拉格朗日乘数μ₁和μ₂，将无偏估计条件和估计方差极小值问题转化为无条件表达式的极值问题求解，获得协同克里金方程组

$ \left\{\begin{array}{l} \sum\limits_{i=1}^{N_{1}} \lambda_{1 i} \gamma_{11}\left(x_{1 i}-x_{p p}\right)+\sum\limits_{j=1}^{N_{2}} \lambda_{2 j} \gamma_{21}\left(x_{2 j}-x_{p p}\right)+\mu_{1}=\gamma_{21}\left(x_{0}-x_{p p}\right) \quad p=1,2, \cdots, N_{1} \\ \sum\limits_{i=1}^{N_{1}} \lambda_{1 i} \gamma_{21}\left(x_{1 i}-x_{p p}\right)+\sum\limits_{j=1}^{N_{2}} \lambda_{2 j} \gamma_{22}\left(x_{2 j}-x_{p p}\right)+\mu_{2}=\gamma_{22}\left(x_{0}-x_{p p}\right) \quad p=1,2, \cdots, N_{2} \\ \sum\limits_{i=1}^{N_{1}} \lambda_{1 i}=0 \\ \sum\limits_{j=1}^{N_{2}} \lambda_{2 j}=0 \end{array}\right. $

(8)

式中γ₂₁为变量2与变量1的交叉变异函数；γ₂₂为变量2的交叉变异函数。

将上述方程组求解出的权重系数λ_1i(i=1，2，…，N₁)、λ_2j(j=1，2，…，N₂)和两个拉格朗日乘数μ₁、μ₂，代入协同克里金法估计公式进行待估计值点的估值。同时，代入估计方差公式，获得简化的协同克里金估计方差

$ \begin{aligned} \sigma_{\mathrm{CK}}^{2}=& \sum\limits_{i=1}^{N_{1}} \lambda_{1 i} \gamma_{21}\left(x_{1 i}-x_{0}\right)+\\ & \sum\limits_{j=1}^{N_{2}} \lambda_{2 j} \gamma_{22}\left(x_{2 j}-x_{0}\right)+\mu_{2} \end{aligned} $

(9)

2 地震反演约束建模过程

地震约束建模的地球物理本质是以地震属性或地震反演等为载体^[15]，与储层物性参数建立关联，借助地震资料横向连续性和钻井资料纵向高分辨率的优势，结合地质资料，开展井震联合储层研究，实现地震约束油藏建模，达到提升建模精度的目的。图 1为基于相控和地震反演约束地质建模技术流程，相控建模是基于地质资料的“硬”约束，地震约束建模主要是基于地震资料的“软”约束，二者结合构成了双“引擎”属性模型^[16]。

一般而言，波阻抗反演是地震约束建模的首选^[17-18]。以孔隙度建模为例，原始孔隙度数据与声波曲线关系密切，受岩心录井、取心费用、实验周期等因素影响，孔隙度数据多由声波曲线(AC或DT)建立图板计算获得^[19]。另一方面，波阻抗曲线通常是由声波曲线与加德纳方程(Gardner Equation)计算获得，即波阻抗曲线与声波曲线存在函数关系^[20]，从而间接把声波曲线与孔隙度建立了关联。实现过程如下：①数据预处理。将反演数据体以SGY格式输入建模软件中，针对不同反射系数进行重采样和特征重构，实现不同数据体在空间和分辨率上的统一，使之与地质建模的网格和测网相适应；②时深转换。时深转换的目的是将时间域的地震反演数据转换成深度域数据体，结合地震合成记录的时深关系，建立空间变速度场，实现地震反演数据体的深度域转化，完成井—震数据空间域匹配，为下一步地震反演约束建模奠定基础；③空间物性概率分析。地震反演约束建模由物性概率曲线实现^[21]，分析研究区河道、坝主体、坝缘、滩砂等微相的孔隙度与波阻抗数据相关关系(图 2)，调整不同小层在不同岩相的函数分布，逐层建立不同岩相下孔隙度与波阻抗的函数关系，从而实现地震反演参与构造建模的“软”约束^[22]。在权重分配方面，引入影响比概念，前文提及的λ₁和λ₂即分别表示地震数据和测井数据的影响权重。

3 技术优势与可信度分析 3.1 技术优势

油藏地质建模一般基于地质统计学或相控约束获得^[23]，基于相控约束的地质建模方法具备融入地质概念^[24]并借助钻井资料纵向高分辨率优势，依据砂体分布有效约束地质模型。然而，针对海上油田，受平台、海域等客观条件限制，各油田井网密度分布大多不均衡，沉积相的绘制主观性较强，从而增大了井间模型多解性的概率^[25]。地震反演约束建模技术恰恰能弥补这一缺陷。以SZ36-1油田东二下段Ⅰ、Ⅱ层为例，基于沉积相约束的相控模型在一定程度上克服了模型随机性，由图 3a、图 3b可见孔隙度空间分布与砂体沉积分布关系密切，受西北物源沉积影响，北部砂体发育区带亦是孔隙度较大区域，南部河道砂体不发育，为大面积滩砂微相，对应孔隙度多为10%~15%。不难发现，高孔砂体分布继承了河道微相的“痕迹”，尤其是北部孔隙度大于30%的高孔砂严重模式化。原因在于，基于井资料的属性建模结果在趋势上与沉积微相具有较好的一致性，但划分沉积相人为主观性较强，井控“痕迹”明显，缺乏横向预测性。而井震约束比单独利用井资料方式综合性更强。基于双约束下的孔隙度模型，在某种程度上消除了硬约束的“痕迹”，高孔砂的分布体不仅符合沉积微相发育趋势，在稀井区具有一定的预测性(图 3c、图 3d)^[26]，经钻井证实，南部个别井区发育数个高孔甜点，特别是D02井区单井钻遇砂岩厚度73.10m，平均孔隙度32.21%，日产油119.94m³/d。由此可见，地震反演约束建模对克服模型多解性具有一定的技术优势。

3.2 可信度评价

模型可信度是评价建模方法可行性、合理性的关键指标^[27]，直接影响建模方法的推广。按井网的排列开展选择性抽稀井，未参与井虚拟为后验井^[28-30]，主要采用两种途径进行评价：

(1) 采用抽稀井网方法评价可信度。对应反演成果分别建立未抽稀井网和抽稀井网地质模型，并对二者进行定性对比，完成基于地震反演约束地震建模的可信度评价^[31]。本文提及的后验井为未参与地震反演和地质建模的预留井，以此提高模型的精度和可信度。图 4为剖面效果对比，可见，无论是否抽稀井网，地震反演剖面的效果差别不大(图 4a、图 4b)，说明地震反演预测较合理。如盲井M17在地震反演参与前后的8号、9号层内均预测一薄砂。总体上盲井参与地震反演后，井间砂体边界更清晰，符合率较高；同理，从孔隙度模型与地震反演剖面可以看出，孔隙度模型剖面总体趋势与地震反演剖面一致，说明地震反演约束的作用较为明显。对比抽稀前后孔隙度模型可见(图 4c、图 4d)，二者差别不明显，基于抽稀井网的孔隙度模型纵向分辨率较高，纵向上可识别5~9套砂体。以M16井为例，抽稀孔隙度模型，井间高孔厚砂体符合较好，孔隙度曲线响应较明显，如6号层和8号层高孔厚砂前后趋势吻合度较高，低孔薄砂体符合一般，砂体越薄，则吻合效果越差。

(2) 利用后验井纵向孔隙度对被抽稀模型开展误差定量分析^[32-33](图 5)可知，后验井平均符合率为92.51%，其中，孔隙度误差低于10%储层多为滩坝薄层砂，该类储层为非含油主力砂体，砂体厚度较薄，发育规模较小，分析原因可能是受地震分辨率限制，薄层砂预测精度低所致。值得一提的是，针对后验井孔隙度误差大于10%这一问题，分析原因如下：一是当砂泥互层且泥岩过薄，含油气砂岩与泥岩产生了调谐效应^[34]，致使阻抗差异小而难以识别，从而降低了地震分辨率；二是误差分析方法的局限性，即孔隙度预测结果为连续数据体，误差分析通常采用均值法计算各小层孔隙度，而钻井孔隙度为离散型点数据，二者难免产生误差。综上得出，地震约束建模具有一定的应用价值。

4 结论

(1) 以SZ36-1油田为例详细阐述了“双控”地质建模原理和实现过程。基于贝叶斯原理，分析了主变量、协变量与被模拟数据的数理关系；通过协克里金法，介绍了实现地震与测井约束地质建模的组合关系。从三个层次介绍了地震反演约束建模的实现过程——数据预处理、时深转换和空间物性概率分析。结合渤海湾盆地SZ36-1油田地质资料，完成了基于相控的“硬”约束和地震反演“软”约束孔隙度模型的建立。

(2)“双控”地质建模技术具有较好的技术优势。结合相控建模技术和地震约束建模技术的效果较好，消除了硬约束的“痕迹”，高孔砂的分布体不仅与沉积微相发育趋势相符合，而且具有一定的预测性。采用定性与定量相结合技术手段，综合评价地质模型的可信度。模型效果不仅符合地质规律的认识，而且后验井平均符合率可达92.51%，说明本文方法建模技术可信度较高。