基于三维数字岩心的介观尺度流弹性响应数值模拟.pdf

书书书第 卷 第期 年月地球物理学报 ，孙超，张怀，唐跟阳等 基于三维数字岩心的介观尺度流弹性响应数值模拟地球物理学报，（）：，：，犆 犺 犻 狀 犲 狊 犲犑犌 犲 狅 狆 犺 狔 狊（），（）：，：基于三维数字岩心的介观尺度流弹性响应数值模拟孙超，张怀，唐跟阳，王尚旭，贺艳晓，潘冬明中国矿业大学资源与地球科学学院，江苏徐州 中国科学院计算地球动力学重点实验室，中国科学院大学，北京 油气资源与探测国家重点实验室，中国石油大学，北京 摘要复杂介质油气藏的地震波频散和衰减是地球物理勘探中油气特征反演与解释的重要属性介观尺度处于微观孔隙尺度与宏观地震波场尺度之间，地震波穿过饱和流体岩石时激励流体发生频散和衰减，导致地震波能量损耗，被称为介观波致流其形成机理、特征和规律一直是勘探地球物理学和岩石物理学等领域共同关心的科学问题，也是研究的前沿和难点之一本文首先选取印第安纳灰岩（）作为研究对象，通过排水方式饱和样品，使用 方法研究介观流体分布特征和规律；为研究其机理，引入 孔弹性理论构建三维（）数值模型，模拟部分饱和岩石的纵波速度频散和衰减特征经对比，数值模拟结果与实验结果匹配，表明 孔弹性理论可以用来描述介观波致流的物理机理本文发展的岩石物理实验结合数值模拟的研究手段为今后研究介观流衰减和频散的特征和规律提供了理论和方法借鉴关键词 数字岩心；中等（介观）尺度流；地震波；频散；衰减 ：中图分类号 收稿日期 ，收修定稿基金项目国家自然科学基金项目（、）；国家重点研发计划（）；油气资源与探测国家重点实验室开放课题（）；徐州市基础研究计划 青年科技人才项目（）资助第一作者简介孙超，男，年生，中国矿业大学副教授 ：通讯作者张怀，男，年生，中国科学院大学教授 ：犇犱 犻 犵 犻 狋 犪 犾 犮 狅 狉 犲犪 狆 狆 犾 犻 犲 犱狋 狅狀 狌 犿 犲 狉 犻 犮 犪 犾 犾 狔狆 狉 犲 犱 犻 犮 狋 犲 犾 犪 狊 狋 犻 犮 狉 犲 狊 狆 狅 狀 狊 犲 犮 犪 狌 狊 犲 犱犫 狔犿 犲 狊 狅 狊 犮 狅 狆 犻 犮 犳 犾 狅 狑 ，犛 犮 犺 狅 狅 犾 狅 犳犚 犲 狊 狅 狌 狉 犮 犲 狊犪 狀 犱犌 犲 狅 狊 犮 犻 犲 狀 犮 犲 狊，犆 犺 犻 狀 犪犝 狀 犻 狏 犲 狉 狊 犻 狋 狔狅 犳犕 犻 狀 犻 狀 犵犪 狀 犱犜 犲 犮 犺 狀 狅 犾 狅 犵 狔，犡 狌 狕 犺 狅 狌犑 犻 犪 狀 犵 狊 狌 ，犆 犺 犻 狀 犪犓 犲 狔犔 犪 犫 狅 狉 犪 狋 狅 狉 狔狅 犳犆 狅 犿 狆 狌 狋 犪 狋 犻 狅 狀 犪 犾犌 犲 狅 犱 狔 狀 犪 犿 犻 犮 狊，犆 犺 犻 狀 犲 狊 犲犃 犮 犪 犱 犲 犿 狔狅 犳犛 犮 犻 犲 狀 犮 犲 狊，犝 狀 犻 狏 犲 狉 狊 犻 狋 狔狅 犳犆 犺 犻 狀 犲 狊 犲犃 犮 犪 犱 犲 犿 狔狅 犳犛 犮 犻 犲 狀 犮 犲 狊，犅 犲 犻 犼 犻 狀 犵 ，犆 犺 犻 狀 犪犛 狋 犪 狋 犲犓 犲 狔犔 犪 犫 狅 狉 犪 狋 狅 狉 狔狅 犳犘 犲 狋 狉 狅 犾 犲 狌 犿犚 犲 狊 狅 狌 狉 犮 犲 狊犪 狀 犱犘 狉 狅 狊 狆 犲 犮 狋 犻 狀 犵，犆 犺 犻 狀 犪犝 狀 犻 狏 犲 狉 狊 犻 狋 狔狅 犳犘 犲 狋 狉 狅 犾 犲 狌 犿，犅 犲 犻 犼 犻 狀 犵 ，犆 犺 犻 狀 犪犃 犫 狊 狋 狉 犪 犮 狋 ，期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟 ，犓 犲 狔 狑 狅 狉 犱 狊 ；引言地震波频散和衰减大致可以分为弹性频散（波场总能量守恒，如散射衰减、几何频散）和非弹性频散（如波场能量转化成热能或应变能）其中，由岩石孔隙中流体引起的非弹性频散是勘探地球物理领域的研究热点当地震波通过时，诱发岩石孔隙中流体流动产生衰减和频散，称为波致流现象（）依据孔隙流体流动尺度的大小，波致流可以分为宏观波致流、微观波致流和介观波致流（，）宏观波致流是指由波长尺度的流体压力梯度引起的流体流动，遵循 孔弹性理论（，）定量描述的固液摩擦而诱发的衰减和频散研究表明，宏观流诱发的频散及衰减通常出现在 以上，远在地震和测井勘探频带之外此外，等（）的研究表明，对于致密储层岩石，存在非达西流，可使宏观流频散和衰减出现在地震频带微观波致流是指由于孔隙形状或方向不同，致使孔隙之间出现局部流体压力梯度而诱发的流体流动，进而产生频散和衰减喷射流是微观波致流中最为熟知的一种，国内外许多学者对其进行了研究（例如，；邓继新等，；欧阳芳等，），结果表明喷射流通常发生在砂岩等粒状岩石中相邻颗粒之间的软孔和硬孔之间此外，巴晶等（，）建立了一种新的黏土喷射流双孔隙度模型，提出发生在黏土的微孔和粒间大孔之间的流体流动同样可以引起显著的衰减微观喷射流可以与 宏观流同时存在，产生双尺度频散和衰减特征（，；杨宽德等，；申义庆和杨顶辉，；，）除宏观流和微观流外，在远大于微观岩石孔隙，远小于宏观地震波长的介观尺度下，当地震波通过饱和多相流体的饱和储层时，由于压缩系数存在差异，多相不相溶流体（例如，水、油和气）受地震波激励会发生震荡式流动（，），简称介观波致流（）现有理论模型（例如，；，；，；，；，；刘炯等，；杜天玮等，；，；，；宗兆云等，）、数值模型（例如，）及实验测量结果（例如，郭梦秋等，；，；，；任舒波等，）表明，介观波致流是诱发地震波速度频散和衰减的重要原因之一在介观波致流理论中，流体分布特征是影响地震波速度频散和衰减的敏感参数与此同时，它也是复杂介质油气藏勘探中油气特征反演与解释的重要目标之一（，；，）然而，在致密薄互层型孔隙弹性地层中，介观流体分布特征的尺度介于微观和宏观之间，致使对其分布规律和特征的探测、监测及反演都异常困难亟需建立介观尺度的流体分布与宏观尺度地震波频散和衰减的量化相关性联系当前，从岩石学实验角度，主要存在两类研究方法表征其相关性联系，即隐性流体分布预测法和显性流体分布预测法隐性流体分布预测法是指忽略流体分布模式细节，以宏观参数（例如，斑块尺度）代表其特征，在孔弹性理论框架下直接推导出动态（频率变化）的混合流体体积模量，并将其代入 （）方程中，即可获得频率相关的弹性模量，进而获得频散和衰减规律（，；，；，）在这类方法中，只需了解饱和度、斑块尺度及频率等参数即可计算流体动态模量，模型使用相对简单，效率高然而，估计实际流体斑块的准确尺度困难，因此多以拟合参数形式解释实测数据，该方法在实际应用中存在较大误差地 球 物 理 学 报（）卷显性流体分布预测法是指流体分布已知，以 框架为理论基础，预测特定流体分布下的频散和衰减如周期层状和球状水气分布模型（，；，；，），以及使用随机介质模型中统计平滑概念（，；，；，）表征的具有随机分布特征的复杂流体分布模型显性流体分布预测法因考虑了流体分布的准确信息而使得预测结果更加准确，但面对复杂的物理方程，无法获得准确的解析解因此，在该领域，数值模拟方法应用更为广泛在饱和条件下，利用有限元方法求解准静态 孔弹性方程，可实现对介观流诱发的频散和衰减的预测（，；，）经验证，数值模拟法可对 的层状模型、的球状模型及随机介质模型正确求解然而，数值模拟方法同样需要流体分布作为输入条件，但实际储层中流体分布情况往往是未知的为此，基于数字岩石物理技术获得样品内各相分布特征，建立数值模型，是一种常用的替代方法，如景惠敏等（，）基于数字样品讨论了热导率变化及岩石应力应变关系等数字岩石物理技术旨在基于高分辨率岩石图像（扫描获得）建立模型，以此揭示岩石各种物理性质例如，通过数字岩石物理获得孔隙空间几何形状（，；，）、水力渗透率（，；，）、电导率（，；包宏帅等，）和弹性岩石模量等（，；，；，；，）此外，利用数字岩心技术获得实际岩石中流体分布是直观且可靠的技术，如 等（）通过 扫描结果粗糙地计算了 灰岩中的流体分布特征，林青杨等（，）和王山等（，）则在微观尺度上给出了不同饱和方式下的流体分布特征 图像虽然刻画的是微观尺度，但通过统计学方法，可以实现利用微观属性预测更高尺度的属性例如，刘洁等（，）利用等效体积单元结合统计学方法，通过岩石微观属性分析介观岩石样品的塑性特征由于众多属性均根据同一样品的 图像预测，可方便地建立各个属性之间的关系，这也是数字岩心技术最重要的优势之一尽管利用数字岩心技术可以建立表征流体分布的模型，但将其用于预测地震波频散和衰减仍然缺少有利依据最近，朱伟等（）利用数字岩心技术对非均匀介质波致流流固相对运动刻画，对含裂隙致密岩石的频散和衰减预测，均表明了数字岩心技术适用 于 频 散 和 衰 减 的 预 测此 外，等（）提出了使用 扫描图像建立流体分布，并尝试据此预测介观流频散和衰减特征，遗憾的是其流体分布的计算及预测方法相对粗糙，且其建立的数值模型只具备一维特征，预测结果与测量结果存在较大差距本文作者在 研究结果的基础上提出在 框架下构建数值模型并预测速度频散和衰减特征（，），将一维模型拓展到二维模型，预测结果与实验结果更加接近，但仍然存在不可忽略的差异结果表明，数值模拟过程中输入的二维流体分布并不能完整的表征样品中的流体分布上述理论结果与实验数据的差异性表明，以数值模型预测地震波频散和衰减仍然存在挑战为消除使用二维 图像表征流体分布带来的偏差，本文首先基于岩石物理实验 图像获得介观三维流体分布特征；引入 框架，构建数值模型，利用有限元求解，厘清介观波致流引起的速度频散和衰减物理机理，揭示介观流体空间分布与地震波频散和衰减的定量关系；建立岩石物理实验、数值模拟结果与实际测量数据之间的量化联系，并形成标签，从而为深部介观流体空间分布特征成像与反演建立理论和方法基础流体分布建模 样品的选取和制备本文主要考虑介观波致流导致的岩石弹性频散和衰减为规避宏观流及微观流等频散机制的影响，研究样品选用 灰岩，别名萨拉姆（）石灰岩前人研究结果表明，该岩石是一种纯方解石生物碎屑石灰岩，由碎屑（生物化石碎片和鲕粒）和方解石胶结物组成，粒度小于 ，样品中极少见微裂隙，该样品具备的特征是：（）地震频带内未测量到喷射流现象；（）岩石渗透率较低，预估其 流的特征频率在 左右，远超地震频带；（）未在地震频带内观测到非达西流（，）引起的频散和衰减特征，例如横波速度频散；（）由于岩石样品内的非均质性尺度远小于地震频带波长，散射衰减基本可以忽略基于上述原因，灰岩满足研究需要（，期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟 ；，；，；，；，；，；，）本文所用的 灰岩（见图）与作者前期实验研究使用的是同一样品（，），长，直径，通过称重法实测干燥样品的密度为 ，通过三重量法（，）实测岩石孔隙度为 ，通过流体驱替法获得 差压（围压减去孔压）下岩石渗透率为 ，具体物性见表图带有红色矩形扫描区域的印第安纳灰岩体分辨率为 ，蓝箭头所指示的铝材料确保了在不同饱和度下 扫描位置的一致性 ，表 犐 狀 犱 犻 犪 狀 犪石灰岩、水和空气的特性犜 犪 犫 犾 犲犘 狉 狅 狆 犲 狉 狋 犻 犲 狊狅 犳 狋 犺 犲 犐 狀 犱 犻 犪 狀 犪 犾 犻 犿 犲 狊 狋 狅 狀 犲，狑 犪 狋 犲 狉，犪 狀 犱犪 犻 狉属性 灰岩水空气孔隙度（）长度犔（）直径犇（）渗透率（）排水体积模量犓（）排水剪切模量（）密度（）黏度（）注：干燥纵波模量为 差压（即围压减去孔压）下所测 ：（，）饱和样品通过前期研究可知，以排水方式饱和 灰岩可以在地震频带内观测到频散和衰减（，），因此本文使用同样的排水方式饱和样品排水方法是指利用气将样品中的水排出的方法，与传统利用驱替排水过程不同，本文利用蒸发过程中空气进入样品实现排水，具体步骤为：（）对干燥样品进行称重，获得重量犿；（）将样品放入水槽中，通过抽真空进行负压饱和在完成抽真空过程后可通过适当加压提高饱和速度，饱和完成后，样品称重为犿，据此计算饱和流体含水质量犿犿；（）将完全饱和（含水 ）的样品在室温（）下放置一天（大约 小时），进行蒸发，通过高精度电子称称重，为犿，然后利用犿犿犿犿计算含水饱和度，发现其降低至约（见图）；（）重复步骤（）、（）和（），可将饱和度分别降低到（见图）、（见图）图使用蒸发排水法饱和样品的过程 地 球 物 理 学 报（）卷 流体分布扫描使用 微型 扫描仪（，）扫描饱和后的样品，从而获得排水饱和方式下的流体分布具体的 扫描内容为：干燥样品、完全饱和样品以及通过排水（蒸发）过程获得的部分饱和样品扫描视野中心位于样品中心 处，扫描尺寸为（直径）（长度），扫描区域见图中的红色矩形框所示范围，体素为 在扫描区域的上下边界贴有铝制标记（简称“铝标”，见图蓝色箭头所指区域），用来确保不同饱和度的样品在放入 扫描设备中时处于同一个位置此外，将碘化钾（，钾的碘化物）水溶液作为流体，以增加孔隙中空气和盐水之间的射线衰减强度的对比度为防止溶液在蒸发过程中出现结晶现象，碘化钾质量分数应低于 ，使用塑料箔包裹样品，以避免在扫描过程中出现蒸发现象 图像去噪及空间位置校正使用非局部均值边缘保持滤波器（，；，）对图像进行滤波，以降低图像噪声该滤波器根据每个体素的灰度值和特定窗口内的邻域体素的相似性计算每个体素的权重值，为相似性高的体素分配高权重然后，将窗口内的加权平均值以新的灰度值分配给每个体素使用该滤波器对干燥样品的放大区域（见图）进行滤波，得到滤波后的结果（见图），可以看到，在保持原有细节的同时，观测信号的随机噪声得到有效压制另一方面，对比干燥样品（见图）和含水饱和度 样品（见图）的 扫描图可知，尽管实验图 扫描图像去噪及空间位置校正（）干燥样品的放大区域；（）使用非局部均值边缘保持滤波器对干燥样品放大区域去噪；（）使用非局部均值边缘保持滤波器对部分饱和样品放大区域去噪；（）使用 中的非局部均值边缘保持滤波器和 空间校正模块对部分饱和样品放大区域去噪和空间校正 （）；（）；（）；（）期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟中使用“铝制标记”确保在不同含水饱和度进行 扫描时的空间位置一致，但仍然存在轻微不匹配（见图 和图，中心位置的孔隙形状不同）计算饱和度需要准确的空间位置，因此，以干燥样品 图像上的“铝制标记”作为参考，对部分饱和样品 图像 进 行 空 间 位 置 精 细 校 正，上 述 过 程 可 使 用 中的 模块完成校正结果见图，对比图 可知，部分饱和样品校正后的 图像在空间位置上与干燥样品完美匹配 饱和度计算首先，将 扫描图像的体素由 粗化成 一方面，适当的粗化使等效体积单元（）更具有代表性，如孔隙度更加均匀（见图），减少物性评估误差（，）；另一方面，粗化使数据大幅度减少，可提高计算速度图像粗化后，利用 扫描图像结果计算 样品中的含气饱和度：（）使用 中的 模块将完全饱和样品的灰度值“归一化”为干燥样本的灰度值（，）所谓“归一化”是指将不同 扫描结果中同一材料的灰度值进行统一对于 样品而言，用来包裹样品的胶套和样品上用来标记位置的铝制材料（见图）是已知材料，胶套的灰度值最低，铝制材料的灰度值最高，且二者灰度值不随流体饱和度发生变化，因此适合作为参考材料具体的“标准化”公式如下：犐 （犐犐）犐 犐 犐 犐 犐，（）其中犐 是校正后的灰度值，犐是原始灰度值，犐 是饱和样品上胶套的平均灰度值，犐 是干燥样品上标记材料铝的平均灰度值，犐 是干燥样品上胶套的平均灰度值，犐 是饱和流体样品上铝材料的平均灰度值（）“归一化”后的灰度图像作为输入，通过做差计算方法获得样品孔隙度：犐犐，（）其中犐是样品完全饱和水时的灰度值，犐是样品饱和空气（干燥样品）时的灰度值，是归一化系数，通常是水和空气 数的差值孔隙度计算结果见图 和图，可以看出其孔隙度均匀图 展示了孔隙度随切片的变化，可以看到其值在 之间浮动，表明样品的孔隙度均匀上述孔隙度变化范围与实验测量结果一致（见表），表明了计算结果可靠（）计算含气饱和度犛：犛犐犐 犐犐，（）其中犐 是水气混合岩石的灰度值（即同时包括水和空气）含气饱和度的计算结果见图 和图，其中白色为气，黑色为水可以看出，排水方法获得的流体分布中气泡呈斑块状物理场建立与求解 犅 犻 狅 狋孔弹性方程按照 扫描区域岩石尺寸构建圆柱形多孔材料模型，见图根据 （，）的研究结果，多孔材料中固体和流体应满足动量守恒方程（）和达西定律（），其频率域形式为：（）狌（）犘，（）（）犘狌（）狊犘犕 狌狊，（）其中固体位移矢量狌和流体压力犘是待求解的变量，是哈密顿算子，是散度算子，是角频率，犓犓是 系数，犓是样品的排水体积模量，犓是岩石颗粒的体积模量，（）是饱和流体的岩石密度，其中是流体密度，是基质密度，是孔隙度复密度（）犻 与岩石的孔隙度，流体密度，孔隙迂曲度，流体黏度，岩石渗透率，角频率和虚数单位犻有关犕称为流体模量，表征了流体对岩石骨架的作用是应力张量，与所求变量流体压力犘，固体位移矢量狌有关，其空间各分量表示为：犻 犼犻 犼犻 犼（狌狊 犘），（）其中犻，犼，是欧几里得空间维数，对应空间狓、狔和狕方向，拉梅常数犓，是样品排水剪切模量，犻 犼为 置换张量，犻 犼犻犼犻犼，犻 犼为岩石骨架的应变张量，与固体位移矢量狌的导数有关：犻 犼（狌犻，犼狌犼，犻），（）地 球 物 理 学 报（）卷图（）三维孔隙度分布；（）采用排水方法饱和，含水饱和度 样品内的三维流体分布；（）图（）的犢 犣剖面；（）图（）的犢 犣剖面；（）孔隙度随着犡 犢平面的变化（沿着犣轴）（）；（）；（）犢 犣 （）；（）犢 犣 （）；（）犡 犢 （犣 ）狌狊犻，犼为狌狊的犻分量在犼方向的导数在方程组（）中，参数犕定义为：犕犓犓（），（）其中流体的体积模量犓与流体饱和度相关，对于本文中空气和水混合的情况，空气饱和度犛和水饱和度犛满足犛犛依据 （）公式，水与空气的混合体积模量犓（，）表示为：犓犛犓犛犓（），（）其中，犓和犓分别为空气和水的体积模量此外，在方程组（）中，以及与流体饱和度相关：犛（犛），（）其中，和分别为空气和水的密度混合流体黏度定义为（，）：（）犛，（）其中，和分别为空气和水的黏度期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟相关参数更多信息见表表岩石物理参数符号定义犜 犪 犫 犾 犲犇 犲 犳 犻 狀 犻 狋 犻 狅 狀 狊狅 犳 狊 狔 犿 犫 狅 犾 狊狌 狊 犲 犱犳 狅 狉 狋 犺 犲狆 犲 狋 狉 狅 狆 犺 狔 狊 犻 犮 犪 犾狆 犪 狉 犪 犿 犲 狋 犲 狉 狊符号物理意义取值渗透率约 流体黏度约 孔隙度 干燥样品剪切模量 犓干燥样品体积模量犓犓基质颗粒体积模量 饱和流体的岩石密度（）基质密度 流体密度 犓流体体积模量 复密度犻 孔隙迂曲度 常数犓 系数犓犓犕流体模量犓犓（）犛含气饱和度犛含水饱和度注：表中取值范围参考了 岩石物理手册 中常见储层岩石及流体的取值范围（，），犻为虚数单位 ：（，），犻 网格化使用有限元方法求解上述方程组需要对数据进行网格化对 数据网格化的常规方法是在图像处理（如阈值分割）后，可使用 软件中 模块将立方体数据转换为其他软件可用的网格使用该方法构建网格直观、方便，但当研究目标较小时，产生的网格质量差，或者整体节点编号难于优化，使得有限元计算需要极大的内存，且有时由于研究目标几何形状纵横比较大使得剖分网格过程失败，或在实际计算过程中出现问题本文中，通过 扫描获得的流体分布可知（见图 和图），样品中的气泡相对较小，使用上述划分网格方法无法获得高质量网格，更无法进行稳定的有限元数值计算为此，本文提出一种方法：按照样品尺寸构建均匀背景模型，进行常规四面体网格划分（见图）；根据 图像获得的流体分布中的气泡位置构建网格表达式，利用其对四面体网格进行自适应精细化网格剖分（见图）这种剖分网格方法的优点在于其规避了处理复杂几何模型的过程，如未完美闭合的曲面、大纵横比流体分布等最后利用流体饱和度函数计算方程组（）中对应的参数，即可完成数值模型求解的前置条件 边界条件及求解求解孔弹性方程组（）需要边界条件，在不排水条件下，选择流体第二类边界条件，样品所有边界流量应为，其频率域形式为：狀犻（犘（），（）其中狀是边界外法向向量此外，在 笛 卡 尔 坐 标 系 下，根 据 纵 波 模 量 的定义：样品侧面边界（见图 中蓝色区域）横向位移为，即：狌犻，（）其中，犻，是欧几里得空间维数，对应狓、狔方向样品下边界为固定边界，所有固体位移分量均为，即：狌（）样品上边界（见图）加载周期震荡应力 ，对应的频率域边界条件为：（）根据纵波模量定义，犕为：犕 ，（）其中，为样品垂向（狕方向）应力，为样品轴向应变纵波衰减为：犙犕犐（犕）犚（犕），（）其中，犐和犚分别是虚数和实数算子利用有限元方法对上述物理场方程和边界条件进行求解结果与分析 误差当流体的空间分布已知时，数值模型可直接估计其等效纵波复模量一般来说，如果数值样本足够大，充分表征流体分布的影响，则可以获得与真实岩石近似的平均结果然而，受计算能力限制，通常无法选取足够大的数值样品尺寸为此，可使用 地 球 物 理 学 报（）卷图（）网格粗化；（）自适应细化网格；（）纵波边界条件位置；（）应力位置 （）；（）；（）；（）方法解决上述问题（，）在自然界中，可认为岩石中流体分布具有随机性，例如，等（）使用 自相似相关函数来模拟部分饱和情况下的流体分布因此，基于 方法预测其统计平均特征是合理的具体为，使用 谱密度函数生成一组饱和度相同，流体分布随机的样品集犖犚对于每个数值样品狀，犖犚而言，单个频率犳犿，犿，犖犉下的弹性响应狀（犳犿）是随机函数，随机函数的均值和方差代表了岩石样品集弹性响应的统计性特征数值样品的方差由（）式计算（，）：（犖犚）犖犉犖犉犿（犳犿，犖犚），（）其中，（犳犿，犖犚）犖犚犖犚狀狀（犳犿）（犳犿，犖犚），（）（犳犿，犖犚）犖犚犖犚狀狀（犳犿）（）根据 和 （）建议的标准，当样品集方差稳定时，停止 的数值模拟，此时数值样品的平均值（犳犿，犖犚）为等价复弹性模量图展示了在不同饱和下，依据 方法获得的纵波模量的方差犕，其中测试样品数目犖犚 本文目标频带为地震频带，为突出目标频带内方差的影响，选取频率范围为 到 （，）如需考虑更宽广的目标频带，最大频率不宜超过高频极限（频散上限）对应的最小频率，以防止过量加权高频极限结果导致方差量级过小频率数目在保证预测结果连续的同时不宜过多，本文取犖犉 流体分布通过 函数建立，其中关键参数是非均质性尺度，能够衡量流体分布非均质性的大小图中不同颜色的曲线代表不同非均质性尺度（，）下纵波模量的平均方差总的来说，随着测试样品数目犖犚增加，纵波模量的方差减少，表明足够数量的数值测试是减少方差的关键此外，当试样中非均质性尺度减少时，只需较少的测试样品数量即期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟图纵波模量平均方差随测试样品数量犖犚的变化（）含气饱和度；（）含气饱和度；（）含气饱和度 不同颜色分别代表，的非均质性尺度 犖犚（）；（）；（），可导致纵波模量方差稳定，即数值测试的误差相应减小对比同一非均质性尺度（如图，和 中非均质性为），随着饱和度增加，纵波模量方差减少当含气饱和度增加到 时（见图），纵波模量方差达到 ，几乎为从上述分析可知，非均质性尺度越小，数值模型的误差越小而非均质性尺度的大小具有相对性，当数值样品体积变大时，非均质性相对变小，即非均质性尺度减小此外，流体饱和度同样会引起非均质性特征发生变化那么，在一定饱和度下，多大尺寸的岩石试样才能表征天然岩石的弹性响应呢？景惠敏等（，）在对岩石导热系数的研究中针对该类似问题进行了讨论本文采用相似的研究思路，取样品体积与以非均质性尺度为直径的球体的体积的比值（以下简称体积比）作为变量，探索其影响实验结果见图，其展示了含气饱和度为，和 的情况下，体 积 比 分 别 为，时纵波模量的方差变化在不同含气饱和度下，方差分布呈现相同趋势，即方差随着体积比的增加而减小当体积比增加到一定值时，方差减小的速度变缓，并趋于常值从这些图中可以看出，体积比的临界值在 左右换言之，对于一定的非均质性尺度，要获得较小的实验误差，样品体积应大于流体分布引起的非均质性尺度代表的球体体积的 倍图纵波模量平均方差随体积比的变化（）含气饱和度；（）含气饱和度；（）含气饱和度 体积比定义为样品体积与以异质长度为直径的球体体积之比 （）；（）；（），流体压力随频率变化依据上述误差分析可知，当体积比大于 时，介观尺度样品的数值模拟与宏观尺度结果误差最小在本文中，样品的体积约为 以 扫描结果为约束，依据 密度函数估计含气饱和度为 时流体诱发的非均质性尺度，其值约为，对应体积（假设为球状）为 据此计算体积比为 ，大于临界值 因此，由 结果预测的弹性响应具有很好的代表性地 球 物 理 学 报（）卷依据上述孔弹性方程和边界条件，利用有限元求解即可获得样品内部流体压力分布特征为消除加载应力大小对结果的影响，需将流体压力犘除以加载应力（）进行归一化，图到图 中的流体压力均为归一化的流体压力图为利用三维流体分布，预测的不同频率下流体压力的分布特征，图为图中对应频率下的流体分布特征的犢 犣切面图三维归一化流体压力分布 可以看到：当频率为 时，孔压充分释放，样品内部流体压力约为（见图 及其切面图）；当频率增加到 时，部分含水区域压力未充分释放到含气区域，如样品顶部（压力值约 ，见图）及样品内部（压力值约 ，见图）均出现局部高压；当频率为 时，大部分含水区域压力未充分释放到含气区域，样品顶部和中部含水位置均为高压（压力值约 ，见图），犢 犣切面表明样品内部含水区域同样为高压（压力值约 ，见图），高压流体范围较 时有所增加；当频率为 时，样品中几乎所有含水位置压力均未释放到含气区域，水气各自保持压力（水压力值约为 ，气压力值约为 ，见图），样品内部同样如此（见图）流体压力随频率变化是引起纵波模量频散和衰减的根本原因岩石纵波模量和衰减随频率的变化关系见图 当频率为 时，岩石的纵波模量约 ，衰减约为 相较于岩石排水纵波模量 ，纵波模量几乎未发生变化，表明岩石内的流体几乎未对岩石骨架起支撑作用，流体压力较低，而较小的衰减表明流体未处于流动状态低的流体压力和流体流动性表明岩石内部流体压力已经完全释放，该结论与图 及图 的流体压力分布特征一致当频率增加到 时，岩石纵波模量增加到 ，纵波衰减增加到 纵波模量及衰减增加表明岩石内部流体压力对岩石骨架起到一定的支撑作用，即存在高流体压力，此外，高衰减表明流体流动性较大高流体压力和高流体流动性表明岩石内部流体处于流动状态，未完成压力释放，该结论与图 和图 展示的部分含水区域流体压力偏高现象一致随着频率继续增加到 和 ，岩石纵波模量分别增加到 和 ，衰减则期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟图三维归一化流体压力分布的犢 犣方向横截面 犢 犣 降低到 和 ，模量增加表明岩石内更多高流体压力出现，对岩石骨架起到更多的支撑作用，而逐渐降低的衰减峰值表明流体流动性趋弱高的流体压力和低流体流动性表明此时岩石内部流体处于“封闭”状态，流体压力未得到时间释放，对比图、图、图 和图，可以看到流体压力表现出“斑块化”分布特征 与二维（犇）数值模型的对比为对比使用 流体分布和 流体分布预测频散和衰减的差异，以图 作为流体分布输入到 孔弹性模型中，流体压力结果见图 对比 流体压力分布犢 犣剖面（见图）可知，除（见图 ）频率外，在 （见图 ）、（见图 ）和 （见图 ）时流体压力分布及数值均存在差异，总结为：模型预测的流体压力（见图）的高压范围小于 模型的预测结果（见图）；模型预测结果在含水区域的流体压力值（，见图）略高于 模型在含水区域预测的流体压力值（，见图）上述差异主要来源于两点：模型的流体分布只是一个切面，与整体的流体分布并不相同，含水饱和度可能偏高或者偏低；模型假设犢轴方向处处相同，使得样品成为一个长方体，因此与 圆柱形样品流体分布不同，造成频散和衰减曲线的差异对于本样品而言，流体分布预测所得频散曲线的特征频率（）略低于 流体分布预测结果（），意味着 模型中的流体释放达到压力均衡需要更长时间此外，模型 预 测 的 纵 波模 量 频 散曲 线的 高 频（）结果略高于 模型的预测结果，意味着 流体分布中含水饱和度略高于 流体分布 与理论模型对比模型和模型是纵波模量频散特征的上下边界，模量频散范围应位于二者之间（，）地 球 物 理 学 报（）卷图 样品中心横截面的归一化流体分布 模型 低频极限当频率充分低时，流体压力完全释放，多相流体压力均衡，可看作一种流体，其体积模量可以用 （）方程描述，其形式与方程（）所示一致，将计算出的混合流体体积模量犓代入 （）方程，获得饱和混合流体岩石的体积模量：犓犓犓犓（）犓犓犓犓，（）而纵波模量：犕犓，（）被称为 （）极限，由于上述极限的假设条件是极低频率，因此极限也称作低频极限图 中黑实线为纵波模量低频极限的预测结果，可以看到在频率 以下，低频极限与数值模拟结果是一致的，即对于本样品而言，极低频率为小于 的频率 模型 高频极限当频率充分高时，流体压力没有时间释放，各流体斑块处于“隔离”状态，每个斑块内部均可看作流体完全饱和，因此可以使用 方程（见公式（）计算各个斑块区域不排水体积模量，即饱和水体积模量犓 和饱和气体积模量犓 ，通过 （）公式定义整个样品的纵波模量：犕犛犓 犛犓 烄烆烌烎（）上述模量计算使用了 和 理论，因此被称为极限，因其成立假设条件为极高频率，因此也叫高频极限 预测结果见图 ，可以看到，当频率高于 时，数值模拟结果与预测结果相同，即对于本样品而言，极高频率是指大于 的频率期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟图（）纵波模量和（）纵波衰减随频率变化的曲线实验测量（正方形）、数值预测结果（实线和虚线）和理论模型预测结果（蓝线和黑线）五角星为超声波测量结果实验测量结果来自含水饱和度 的 灰岩，饱和方法为排水 （）（）（），（），（），饱和岩石的连续随机介质模型（，模型）对于介于 和 的中间频率，流体压力处于部分释放状态，样品的纵波模量介于低频极限和高频极限之间预测中间频率的频散和衰减的方法有很多，其中 等（）以 图像计算流 体 分 布 相 关 函 数，利 用 连 续 随 机 介 质 模 型（）预测频散和衰减特征的方法与本文相关性最高，因此将其预测结果与数值结果进行对比在连续随机介质模型中，饱和的复纵波模量为：珮犎（）犎（），（）其中犎是低频极限，犎是高频极限，犎犎犎，犖 犕犕犎，犖犕犔犎，犔犓是干燥纵波模量，犎犓是饱和水的纵波模量，犕为流体模量，如公式（）所示，犕犕是犕的归一化方差，引起纵波模量频变效应的是，（）犽?狉 犅犕（狉）犽 狉狉，（）其中犽 槡犖是慢波数，犅犕（狉）是表征岩石流体分布特征的相关函数公式（）和公式（）中的其他未说明参数参照表通过公式（）可知，频变（）是关键函数，其核心是互相关函数 等（）指出，可以通过 扫描图像结果获得流体分布的互相关函数，该函数可以用传统的 和 （）函数表征对于本文的流体分布特征，拟合 函数的相关长度是，据此计算纵波模量频散和衰减见图 中蓝色曲线可以看出，其与 数值模型中的预测结果类似，衰减特征频率一致，约 ，但衰减峰值更低，约 ，与 数值模型预测结果相比，特征频率更高，衰减峰值更小 与实验结果对比为了解 数值模型预测结果的精度，使用实测数据进行对比分析本文实测纵波模量和衰减实验数据 来 自 于 作 者 前 期 研 究 结 果（，地 球 物 理 学 报（）卷 ），该成果中使用的 灰岩与本文所用灰岩为同一块，饱和流体步骤及饱和方法（排水）与本文所述完全一致，因此可用该实验数据与各个模型预测结果进行对比实验所用测量方法为应力应变法和行波法，测量频带为 和 测量设备为巴黎高师研发的低频应力应变设备（，），基本原理是使用岩石样品的轴向应变 和径向应变，获得岩石的纵横比及纵横比相位差犙 狏：，（）犙 狏 （），（）利用岩石轴向应变、参考铝块轴向应变 及各自相位，获得岩石的杨氏模量犈 及杨氏衰减犙：犈 犈 ，（）犙 犈 （），（）其中、分别代表应变及相位，犈是杨氏模量，是岩石泊松比，犙是衰减，上标 和 分别代表铝参考样品和岩石样品，下标，分别代表轴向和径向横波衰减犙 获得方法参照 和 （）的研究结果：犙 狏狏犙（狏）犙 犙）（）狏 犙（狏）犙，（）最终得到纵波模量犘和衰减犙：犘犈（狏）（狏）（狏），（）犙 狏（狏）（狏）犙 狏（狏）（狏）（狏）犙（）测量结果见图 中红色方框数据点，纵波模量频散过渡带为 ，衰减特征频率约为 ，衰减峰值约为 ，超声获得的纵波模量为 对于超声结果，模型和 模型的预测结果与实验完全一致，数值模型和 模型预测结果略低，但差异不明显对于低频结果，当频率小于 时，各模型预测结果与实验结果均相符当频率在 之间时，对于纵波模量，模型与实验 结 果 一 致，频 散 可 达 左 右，而 及 模型预测结果整体偏低对于纵波衰减，模型预测的峰值和特征频率与实验一致，即在 左右存在 的衰减峰值，而 模型和 预测的衰减特征频率比实测结果高倍，约 左右，衰减峰值偏小，约为 通过上述结果的对比分析可知，模型预测结果与实验结果一致，较其他模型预测结果精度更高，成功解释了流体分布与频散和衰减的定量关系结论本文详细而完整的介绍了使用岩石样本的三维 图像定量估计介观流频散和衰减的方法为此，首先以排水方法获得非均质流体分布，通过三维 图像数据，将流体分布进行数字化，构建数值模型，将其输入到 孔弹性模型中，使用有限元方法计算相应流体分布下的纵波速度频散和衰减特征依据三维数值模型的预测结果得到如下结论：（）岩石样品体积越大，流体非均质性尺度越小，数值模型的预测结果越具有代表性当“体积比”大于 时，数值模型预测结果的误差最小，且达到稳定；（）流体饱和度影响数值模型的预测误差，高含水饱和度使模型预测误差更大，但令模型误差达到稳定的“体积比”基本不变；（）对比传统的二维数值模型及 模型，三维数值模型预测的纵波模量频散更高，衰减特征频率更低，与实验结果完全一致，表明精确的流体分布是预测介观流频散和衰减的关键；（）三维数值模型及实验均表明，介观流诱发的纵波模量频散可出现在地震频带，对 灰岩而言，频散幅度达 左右，显著而不可忽略总的来说，基于三维数字岩心的数值模型消除了二维数值模型预测结果与实验结果之间的差异，确认了流体分布和地震衰减之间的定量关系（孔弹性理论），为解释实验室介观尺度频散和衰减的测量结果提供了可靠依据致谢感谢匿名审稿人对本文的有益建议犚 犲 犳 犲 狉 犲 狀 犮 犲 狊 ，犌 犲 狅 狆 犺 狔 狊 犻 犮 犪 犾 犚 犲 狊 犲 犪 狉 犮 犺犔 犲 狋 狋 犲 狉 狊，（）：，：，：犌 犲 狅 狆 犺 狔 狊 犻 犮 狊，（）：期孙超等：基于三维数字岩心的介观尺度流弹性响应数值模拟 ，：，犛 犘 犈犑 狅 狌 狉 狀 犪 犾，（）：，：，犌 犲 狅 狆 犺 狔 狊 犻 犮 犪 犾犚 犲 狊 犲 犪 狉 犮 犺犔 犲 狋 狋 犲 狉 狊，（）：，犑 狅 狌 狉 狀 犪 犾 狅 犳犃 狆 狆 犾 犻 犲 犱犘 犺 狔 狊 犻 犮 狊，（）：，：犜 犺 犲犑 狅 狌 狉 狀 犪 犾 狅 犳狋 犺 犲犃 犮 狅 狌 狊 狋 犻 犮 犪 犾犛 狅 犮 犻 犲 狋 狔狅 犳犃 犿 犲 狉 犻 犮 犪，（）：，：犑 狅 狌 狉 狀 犪 犾 狅 犳犃 狆 狆 犾 犻 犲 犱犘 犺 狔 狊 犻 犮 狊，（）：，：，：犑 狅 狌 狉 狀 犪 犾 狅 犳犌 犲 狅 狆 犺 狔 狊 犻 犮 犪 犾犚 犲 狊 犲 犪 狉 犮 犺：犛 狅 犾 犻 犱 犈 犪 狉 狋 犺，（）：，：，犚 犲 狏 犻 犲 狑狅 犳犛 犮 犻 犲 狀 狋 犻 犳 犻 犮 犐 狀 狊 狋 狉 狌 犿 犲 狀 狋 狊，（）：，：，（），：，：，犌 犲 狅 狆 犺 狔 狊 犻 犮 狊，（）：，：，：犑 狅 狌 狉 狀 犪 犾 狅 犳犌 犲 狅 狆 犺 狔 狊 犻 犮 犪 犾犚 犲 狊 犲 犪 狉 犮 犺：犛 狅 犾 犻 犱犈 犪 狉 狋 犺，（）：，：，犑 狅 狌 狉 狀 犪 犾 狅 犳犃 狆 狆 犾 犻 犲 犱犘 犺 狔 狊 犻 犮 狊，（）：，：，犆 犺 犻 狀 犲 狊 犲 犑 狅 狌 狉 狀 犪 犾