CO2地质封存中四维多分量地震监测的转换波AVO计算方法

摘要

CO2地质封存中四维多分量地震监测的转换波AVO计算方法，模拟CO2注入与地质封存的不同阶段，储层孔隙压力、CO2饱和度及裂缝密度变化情况下，不同入射角和不同方位所对应的PS1转换波和PS2转换波AVO；针对储层为各向同性介质与各向异性介质，在注入CO2的不同阶段，在储层孔隙压力和混合流体饱和度变化情况下，储层弹性参数变化及AVO变化；针对HTI介质模型，模拟注入井压力增大可能导致的裂缝张开，及生产井压力减小导致的裂缝闭合，针对孔隙度与裂缝密度变化；分析裂缝密度变化导致的注入CO2前后各向异性介质中，四维转换波AVO及方位AVO的变化；有助于标定和解释四维转换波地震资料，确定CO2驱油效果、CO2在地下的封存状态与分布范围和评价CO2地质封存的安全性。

说明书

技术领域

本发明专利涉及CO₂驱油、地质封存中四维多分量地震监测技术领域，具体涉及CO₂地质封存中四维多分量地震监测的转换波AVO计算方法，用于模拟CO₂注入地下不同阶段， 四维多分量地震监测中转换波AVO地震响应及四维转换波差异响应。

背景技术

自从挪威科技大学的Landro等(2003)开展四维地震的PP与PS波AVO联合反演砂岩 储层压力与流体饱和度研究取得成功后，四维多分量转换波地震在国际上的油气勘探中得 到越来越多的应用。这主要是因为单独利用常规纵波进行四维地震监测，难以区别储层流 体饱和度变化与压力变化产生的振幅差异，因此难以通过地震振幅差异判断那些振幅差异 异常是流体饱和度产生的、那些是压力变化产生的？由于转换横波对于流体不敏感，所以四 维转换波的差异主要是压力产生的。这样，通过PP与PS的联合解释，我们可以有效区别压力 变化与饱和度变化在储层的分布范围。

国际上在CO₂-EOR(CO₂EnhancedOilRecovery)与地质封存领域，利用四维多分量 地震资料来监测压力、流体饱和度或裂缝变化已成为一种趋势。加拿大Weyburn油田CO₂地 质封存项目开展了世界上第一次完善的四维多分量地震监测，到目前为止在其Phase1A工 区采集了5次三维三分量(3D3C)地震资料。在Weyburn油田的CO₂注入区Phase1B工区，科罗 拉多矿院的ThomasL.Davis(2003)教授采集了三次开展四维3D9C地震，用于监测CO₂地质 封存过程中的裂缝分布，及CO₂泄漏风险评估。此前，ThomasL.Davis(2002)教授在美国新 墨西哥州VacuumField开展了四维3D9C监测注水驱油的研究，在科罗拉多州非常规气田 RulisonField进行盖层的四维多分量监测，和进行四维多分量谱分解研究。并在 Louisiana的DelhiField运用四维多分量技术开展残余油饱和度与盖层完整性研究，特别 是CO₂注入过程中横波分裂如何用于监测盖层应力状态。Bale(2005，2013)等在EOR监测重 油过程中，采用四维多分量中横波分裂研究应力状态，Grossman(2013)研究了加拿大 Alberta省西北PouceCoupeField油田，非常规油气藏多分量四维地震中横波分裂与油气 产量之间的关系。研究表明岩石(包括储层)所受水平最大应力方向对应产生的裂缝方向。 我们知道一般顺裂缝方向速度快而垂直裂缝方向速度慢，横波遇到裂缝可以发生分裂。 Dillen(2000)发现了应力各向异性也可以导致横波分裂。而国内目前还没有开展四维多分 量地震资料的实际监测。

目前，国际上的四维多分量相关研究大多考虑是基于油、水两相介质，并且不考虑 储层开发过程中注入井、生产井压力的变化。即便在CO₂地质封存与CO₂-EOR四维多分量监测 中，也没有分析压力变化对储层弹性参数的影响。即对于各向同性介质，无论注入CO₂压力 和生产井压力多大，储层弹性参数在注入CO₂的不同阶段都是相同的。对于各向异性HTI介 质，则不考虑注入压力及生产井压力变化对于裂缝密度的影响。

本发明考虑了注入CO₂过程中，注入井压力孔隙高和生产井孔隙压力低导致的各 向同性介质储层弹性参数变化。同时考虑了注入井压与生产井压力变化对于垂向裂隙介质 (HTI)的影响，即裂缝密度在不同压力下的变化。并分别依据各向同性介质和各向异性介质 (HTI)理论，进行注入CO₂不同阶段四维转换波地震AVO响应的模拟，使得模拟结果更符合实 际注入状况。为准确的四维多分量地震解释奠定基础。

AVO为AmplitudeVersusOffset的简写。

发明内容

为克服上述现有技术的不足，本发明的目的在于提供，CO₂地质封存中四维多分量 地震监测的转换波AVO计算方法，CO₂注入地下进行驱油和地质封存的不同阶段，模拟四维 多分量转换波地震响应的方法。可以从转换波地震资料正演的角度解释四维转换波差异产 生的原因，本方法分别考虑储层为各向同性均匀介质和各向异性HTI介质两种情况，针对各 向同性的孔隙介质，本方法基于Gassman流体替换理论，在流体替换过程中同时考虑孔隙压 力和CO₂饱和度变化；运用岩石物理测试获得的干岩石弹性模量与切变模量随压力变化关 系，对干岩石体变模量和剪切模量做校正。

对于各向异性HTI介质，在流体替换过程中加入储层裂缝密度随CO₂注入压力变化 的关系，即考虑流体替换过程中储层裂缝随压力变化的张开或闭合；然后，引入Thomson弱 各向异性理论，量化CO₂注入不同阶段各向异性参数变化规律；最后，采用Rüger(1998)对平 面波入射条件下各向异性HTI介质中的简化公式，计算CO₂注入不同阶段，及不同入射角和 方位角下PS1波和PS2波AVO的变化。

为了实现上述目的，本发明采用的技术方案是：CO₂地质封存中四维多分量地震监 测的转换波AVO计算方法，包括以下步骤：

1)从测井曲线中，得到CO₂注入前的纵波速度、横波速度、密度、孔隙度参数的测井 曲线，结合阿尔奇公式得到含水饱和度和含油饱和度，从岩心的岩石物理实验中获取干岩 石弹性模量随压力变化关系；

2)以储层岩石物理参数和测井曲线为基础，运用Batzle-Wang公式计算未发生流 体替换时，未注入CO₂，储层条件下油气水混合流体的体变模量和密度，然后通过Gassmann 方程，求取出储层段干岩石随深度变化的岩石物理参数曲线；

3)将岩心干岩石的岩石物理测试数据，对步骤2)通过测井曲线求出的干岩石随深 度变化的岩石物理参数曲线，进行校正，并带入干岩石弹性模量随压力变化关系，获得储层 段随压力变化的干岩石弹性模量曲线，即干岩石弹性模量随深度变化曲线；

4)对于各向同性的孔隙介质，以Gassmann流体替换理论为基础对储层做注入CO₂后的流体替换，并将随压力变化的干岩石体变模量和剪切模量带入Gassmann方程，计算孔 隙压力和油气、水、CO₂混合流体饱和度变化，及不同温度、矿化度、气油比、API值条件下的 储层弹性参数；

5)对于各向同性介质，在上述流体替换的基础上，以精确的Zoeppritz方程或者 Zoeppritz近似公式计算转换波AVO；

6)对于各向异性介质，采用Tod理论计算压力变化造成的裂缝密度变化及岩石孔 隙度变化，及模拟注入CO₂前原始储层压力、注入CO₂过程中，注入CO₂饱和度增大、压力增大 及生产井压力小条件下的压力变化导致的裂缝密度及孔隙度变化，再将孔隙度变化带入 Gassmann公式，加入混合流体饱和度变化参数后，计算CO₂注入前、后注入井及生产井的储 层弹性物性参数；

7)以Thomsen弱各向异性理论为基础，求取储层各向异性参数以及沿着裂缝方向 的纵、横波速度；

对于各向异性HTI介质，以Rüger简化的方位AVO理论为基础，在Thomson各向异性 参数的基础上计算方位AVO下的各向异性参数，以新的各向异性参数为基础，求取不同入射 角、方位角下的PS1波和PS2波两层介质模型AVO曲线，模拟不同CO₂注入压力和混合流体饱 和度下的PS1波和PS2波AVO曲线；

8)求取不同注入压力条件下、CO₂流体替换前、后井模型AVO反射系数，并以褶积模 型为基础分别制作井模型的各向同性介质和各向异性(HTI)介质转换波人工合成地震记 录；对于各向同性介质模型，将测井资料的深度域采样间隔如0.1米或0.125米看作各向同 性介质，从测井曲线起点开始，以Zoeppritz方程计算不同入射角度的平面波反射系数，然 后与子波褶积获得转换波人工合成地震记录；对于各向异性HTI介质，在储层顶面按照Rü ger简化的方位AVO计算公式，计算不同方位的AVO反射系数；在储层之外，仍然采样各向同 性介质方法计算转换波AVO反射系数，并计算注入CO₂前后转换波人工合成地震记录的差 异；进而模拟不同注入CO₂阶段的转换波地震响应及四维地震转换波差异地震响应，为准确 解释CO₂地质封存中四维转换波信息奠定基础。

本发明建立的计算四维多分量转换波反射系数与AVO方法，其关键在于分别建立 了四维转换波AVO响应与孔隙压力、流体饱和度、裂缝等的联系，使得模拟的转换波地震响 应更符合实际。

本发明的有益效果是：

对于各向同性介质，这里指孔隙介质，本发明通过对Gassman流体替换理论的改 进，模拟CO₂注入地下进行EOR和地质封存过程中，储层孔隙压力和混合流体饱和度变化条 件下的四维地震转换波AVO响应。

对于各向异性HTI介质，通过模拟CO₂注入地下进行EOR和地质封存过程中，除了储 层流体变化外，储层裂缝的张开和闭合造成的裂缝密度变化产生的不同方位、不同注入阶 段转换波AVO的变化，为利用四维转换波地震资料解释和评价CO₂-EOR效果与地质封存的安 全性提供依据。

附图说明

图1加拿大Weyburn油田CO₂封存区1020423井储层部分原始测井曲线。

图2岩石物理实验结果，岩心干岩石物性参数随压力的变化趋势。

图3储层段干岩石随深度变化的岩石物理参数曲线。

图4运用岩心校正后储层段干岩石随深度变化的岩石物理参数曲线。

图5为各向同性均匀介质下，一定孔隙压力时，不同CO₂饱和度(SCO₂)的PS转换波 AVA曲线图；其中图5(a)为孔隙压力为9MPa时，不同CO₂饱和度时转换波的AVA曲线图；图5 (b)为孔隙压力为15MPa时，不同CO₂饱和度时转换波的AVA曲线图；图5(c)为孔隙压力为 21MPa时，不同CO₂饱和度时转换波的AVA曲线图。

图6为各向同性均匀介质下，一定CO₂饱和度(SCO₂)时，不同孔隙压力的PS转换波 AVA曲线；其中图6(a)为CO₂饱和度为0时，不同孔隙压力时的AVA曲线图；其中图6(b)为CO₂饱 和度为20％时，不同孔隙压力时的AVA曲线图；其中图6(c)为CO₂饱和度为40％时，不同孔隙 压力时的AVA曲线图。

图7为随孔隙压力和CO₂饱和度变化时裂缝密度的等值线图。

图8(a)为裂缝密度随孔隙压力的变化曲线图；图8(b)为储层总孔隙度随孔隙压力 的变化图；图8(c)为储层总孔隙度随裂缝密度的变化关系图。

图9为各向异性(HTI)介质模型中，随孔隙压力和CO₂饱和度变化时饱和岩石体变 模量(单位GPa)的等值线图。

图10为各向异性(HTI)介质模型中，随孔隙压力和CO₂饱和度变化时饱和岩石剪切 模量(单位GPa)的等值线图。

图11为各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度(SCO₂)下，随方位角 变换的PS1转换波的AVA曲线图；其中图11(a)为孔隙压力为9MPa，CO₂饱和度为0％时，方位 角分别为0°、30°、60°、90°时PS1波的AVA曲线图；图11(b)为孔隙压力为15MPa，CO₂饱和度为 20％时，方位角分别为0°、30°、60°、90°时PS1波的AVA曲线图；图11(c)为孔隙压力为21MPa， CO₂饱和度为40％时，方位角分别为0°、30°、60°、90°时PS1波的AVA曲线图。

图12为各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度(SCO₂)下，随入射角 和方位角变化时PS1转换波反射系数的等值线图；其中图12(a)为孔隙压力为9MPa，CO₂饱和 度为0％时，PS1波反射系数随方位角和入射角的变化图；图12(b)为孔隙压力为15MPa，CO₂饱和度为20％时，PS1波反射系数随方位角和入射角的变化图；图12(c)为孔隙压力为 21MPa，CO₂饱和度为40％时，PS1波反射系数随方位角和入射角的变化图。

图13为各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度(SCO₂)下，随方位角 变换的PS2转换波的AVA曲线图；其中图13(a)为孔隙压力为9MPa，CO₂饱和度为0％时，方位 角分别为0°、30°、60°、90°时PS2波的AVA曲线图；图13(b)为孔隙压力为15MPa，CO₂饱和度为 20％时，方位角分别为0°、30°、60°、90°时PS2波的AVA曲线图；图13(c)为孔隙压力为21MPa， CO₂饱和度为40％时，方位角分别为0°、30°、60°、90°时PS2波的AVA曲线图。

图14为各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度(SCO₂)下，随入射角 和方位角变化时PS2转换波反射系数的等值线图；其中图14(a)为孔隙压力为9MPa，CO₂饱和 度为0％时，PS2波反射系数随方位角和入射角的变化图；图14(b)为孔隙压力为15MPa，CO₂饱和度为20％时，

PS2波反射系数随方位角和入射角的变化图；图14(c)为孔隙压力为21MPa，CO₂饱 和度为40％时，PS2波反射系数随方位角和入射角的变化图。

图15为各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度下，随孔隙压力和 CO₂饱和度变化时PS1转换波反射系数等值线图。

图16为各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度下，随孔隙压力和 CO₂饱和度变化时PS2转换波反射系数的等值线图。

图17为各向异性(HTI)介质模型中，一定孔隙压力下，不同CO₂饱和度(SCO₂)时PS1 转换波和PS2转换波的AVA曲线图；其中图17(a)为孔隙压力为9MPa时，不同CO₂饱和度时PS1 波和PS2波反射系数的AVA曲线图；图17(b)为孔隙压力为15MPa时，不同CO₂饱和度时PS1波 和PS2波反射系数的AVA曲线图；图17(c)为孔隙压力为21MPa时，不同CO₂饱和度时PS1波和 PS2波反射系数的AVA曲线图。

图18为各向异性(HTI)介质模型中，一定CO₂饱和度(SCO₂)下，不同孔隙压力时PS1 转换波和PS2转换波的AVA曲线图；其中图18(a)为CO₂饱和度为0％时，不同孔隙压力下PS1 波和PS2波反射系数的AVA曲线图；图18(b)为CO₂饱和度为20％时，不同孔隙压力下PS1波和 PS2波反射系数的AVA曲线图；图18(c)为CO₂饱和度为40％时，不同孔隙压力下PS1波和PS2 波反射系数的AVA曲线图。

图19为各向同性均匀介质下，CO₂饱和度为0％时，不同孔隙压力下入射角度道集 转换波人工合成地震记录及其差异。

图20为各向同性均匀介质下，孔隙压力为15MPa，不同CO₂饱和度下入射角度道集 转换波人工合成地震记录及其差异。

图21为各向异性(HTI)介质模型中，方位角为30°时，CO₂饱和度为0％时，不同孔隙 压力下入射角度道集PS1波人工合成地震记录及其差异。

图22为各向异性(HTI)介质模型中，方位角为30°时，孔隙压力为15MPa，不同CO₂饱 和度下入射角度道集的PS1波人工合成地震记录及其差异。

图23为各向异性(HTI)介质转换波人工合成地震记录与Weyburn油田实际径向转 换波叠加偏移资料的对比。

具体实施方式

下面结合附图对本发明作进一步详细说明。

CO₂地质封存的四维地震监测中转换波AVO的计算方法，包括下列步骤：

第一，对工区测井曲线进行整理，得到CO₂注入前的纵波速度、横波速度、密度、孔 隙度等参数的测井曲线，结合阿尔奇公式得到含水饱和度曲线，如图1是加拿大Weyburn油 田1020423井整理后储层部位测井曲线。通过对储层进行岩心采样，并对岩心做岩石物理实 验，获得岩心干岩石的岩石物理参数，表1和表2分别为Weyburn油田储层岩心的岩石物理参 数，图2是岩石物理实验过程中岩心干岩石物性参数随压力的变化情况；

图1加拿大Weyburn油田1020423井整理后储层部位测井曲线，从左到右分别为纵 波速度、横波速度、密度、孔隙度、密度曲线。储层段深度为1398m～1423m，储层段可以分别 为Marly层和Vuggy层，其中Marly层的深度为1398m～1405m，Vuggy层深度为1405m～1423m。

图2岩心岩石物理实验结果，测试了岩心的岩石物理参数随着压力变化的情况，从 左到右分别为纵波速度、横波速度、体变模量、剪切模量、泊松比、扬氏模量随着孔隙压力从 7MPa～24MPa变化时的变化规律。

表1为两层介质模型原始参数(孔隙压力＝15MPa,上覆压力32MPa)。

表1

表2为储层原始物性参数。

表2

名称 物性参数数值 温度T(℃) 56 含盐度(％) 8.5 原始裂缝密度Yc 0.03 原始裂缝孔隙度fic 0.000041448 孔喉比ar 0.00033 储层岩石骨架的体变模量K_m83GPa 储层岩石骨架的剪切模量μ_m48GPa 储层岩石颗粒密度den_grain2.87g/cm³

表3为Weyburn油田干岩石物性参数。

表3

名称 干岩石物性参数数值 纵波速度(Vp_dry) 3.531519km/s 横波速度(Vs_dry) 2.130550km/s 密度(den_dry) 2.035582g/cm³体变模量(K_dry)13.037017GPa 剪切模量(μ_dry) 9.24GPa

表4为计算的随孔隙压力、流体饱和度变化的储层物性参数。

表4

第二，以第一步获得的储层岩石物性参数和测井曲线为基础，包括纵波、横波、密 度、孔隙度、含水饱和度，运用Batzle-Wang公式计算未发生流体替换时(未注入CO₂时)，储 层条件下油气水混合流体的体变模量和密度，然后通过Gassmann方程，获得储层段干岩石 随深度变化的岩石物理参数曲线，岩石物理参数包括纵波速度、横波速度、密度、体变模量、 剪切模量；

在此引入科罗拉多矿院Brown(2002)在Weyburn油田岩石物理测试获得的，关于体 变模量剪切模量与孔隙压力之间的关系：

K_DRY＝1.731*10^-4*(P_d-P_f)²-1.325*10^-2*(P_d-P_f)²+0.3542×(P_d-P_f)+13.11(1)

μ_DRY-1.157*10^-4*(P_d-P_f)³-9.828*10^-2*(P_d-P_f)²+0.2989*(P_d-P_f)+8.78(2)

式中K_DRY为压力校正后干岩石的体变模量，μ_DRY为压力校正后的干岩石剪切模量，P_d是储层围压，P_f是孔隙压力，运用式(1)和(2)校正干岩石体变模量和切变模量。图3描述了 测井曲线求出的干岩石物理参数随深度变化的曲线；

图3描述了由测井曲线求出的储层段干岩石物理参数随深度变化的曲线，从左到 右分别为储层段干岩石的纵波速度、横波速度、体变模量、剪切模量、泊松比、扬氏模量曲 线。

第三，将第一步岩心求取干岩石的岩石物理参数，对第二步通过测井曲线求出的 干岩石随深度变化的岩石物理参数曲线进行校正(如图4)，并带入干岩石弹性模量随压力 变化关系，获得储层段随压力变化的干岩石弹性模量曲线，即干岩石弹性模量随深度变化 曲线；图4左边两列为饱和岩石纵横波速度，其中虚线为原始测井数据经流体替换后所得的 新的测井数据，实线为岩石物理测试数据校正后的饱和岩石的纵横波速度，右侧四列岩石 物理参数曲线中虚线为从测井数据干岩石弹性模量随深度变化曲线，实线为运用岩心求取 干岩石的岩石物理测试数据校正后的干岩石弹性模量以及饱和岩石纵横波速度随深度变 化曲线；

第四，在各向同性的孔隙介质中，不考虑流体替换过程中储层孔隙度的变化，在上 述计算结果的基础上，重复运用Batzle-Wang公式和Gassman方程，求取流体替换过程随孔 隙压力和流体饱和度变化的混合流体体变模量，混合流体密度，进而求取饱和岩石的体变 模量、剪切模量，求取饱和岩石的纵波速度、横波速度、密度、体变模量、剪切模量等新的岩 石物理参数，具体计算过程见式(3)～(7)。

$K_{sub}=K_{DRY}+\frac{{(1-\frac{K_{DRY}}{K_m})}^2}{\frac{porosity}{K_{fmix}}+\frac{1-porosity}{K_m}-\frac{K_{DRY}}{{K_m}^2}}---\left(3\right)$

den_sub＝dengrain*(1-porosity)+den_fmix*porosity(4)

mu_sub＝μ_DRY(5)

${\mathrm{Vp}}_{sub}=\sqrt{\frac{K_{sub}+\frac{4}{3}{\mathrm{mu}}_{sub}}{{\mathrm{den}}_{sub}}}---\left(6\right)$

${\mathrm{Vs}}_{sub}=\sqrt{\frac{{\mathrm{mu}}_{sub}}{{\mathrm{den}}_{sub}}}---\left(7\right)$

在式(3)～(7)中，K_sub是饱和岩石的体变模量，K_DRY是压力校正后干岩石的体变模 量，K_m是储层岩石骨架的体变模量，porosity是储层总孔隙度，K_fmix是混合流体体变模量， den_sub是饱和岩石密度，dengrain是储层岩石颗粒密度，den_fmix是储层流体密度，mu_sub是饱 和岩石的剪切模量，μ_DRY是压力校正后干岩石的剪切模量，Vp_sub是饱和岩石纵波速度，V_sub是 饱和岩石横波速度。

同时，计算孔隙压力和油气、水、CO₂混合流体饱和度变化，及不同温度、矿化度、气 油比、API值条件下的储层弹性参数；

第五，以精确的Zoeppritz方程为理论基础，利用第四步所计算的储层岩石物理参 数，在此基础上计算流体替换过程中不同储层孔隙压力和储层CO₂饱和度下的转换波AVA曲 线，如图5和图6；

如图5所示，模拟了当孔隙压力分别为9、15、21MPa时，即对应Weyburn油田CO₂注入 区生产井、原始地层和注入井压力，CO₂饱和度为0％、20％、40％时孔隙介质中PS转换波的 AVA曲线。图5(a)为孔隙压力为9MPa时，不同CO₂饱和度时转换波的AVA曲线图；图5(b)为孔 隙压力为15MPa时，不同CO₂饱和度时转换波的AVA曲线图；图5(c)为孔隙压力为21MPa时，不 同CO₂饱和度时转换波的AVA曲线图。从图5中可以看到，不同孔隙压力和CO₂饱和度构成的流 体模型所对应的AVA曲线的变化趋势一致。在入射角较小时，PS转换波反射系数随着入射角 的增加而增加，当入射角达到一定值后，反射系数达到最大值，并且所有模型的PS转换波反 射系数在相同的入射角达到最大值。比较相同压力下不同CO₂饱和度时的PS转换波AVA曲 线，我们可以看到当压力一定时，在相同的入射角下，PS转换波的反射系数基本不随CO₂饱 和度的变化而变化，即在孔隙介质模型中，CO₂饱和度对于PS转换波反射系数影响很小；

如图6，描述了在CO₂饱和度一定时，不同孔隙压力下的PS转换波的AVA曲线。其中 图6(a)为CO₂饱和度为0时，不同孔隙压力时的AVA曲线图；其中图6(b)为CO₂饱和度为20％ 时，不同孔隙压力时的AVA曲线图；其中图6(c)为CO₂饱和度为40％时，不同孔隙压力时的 AVA曲线图。横向比较图6中3幅图，在相同CO₂饱和度下和一定入射角时，PS转换波反射系数 随着孔隙压力的增加而增加，孔隙压力越大，反射系数增加越快；但对与不同CO₂饱和度下 PS转换波反射系数变化趋势大致相同(见图5)，这表明在孔隙介质的两层模型中，CO₂饱和 度的变化对于PS转换波的反射系数影响有限，孔隙压力的变化在PS转换波的反射系数的变 化过程中占主导地位；

第六，对于各向异性介质，采用Tod理论计算压力变化造成的裂缝密度变化及岩石 孔隙度变化，及模拟注入CO₂前原始储层压力、注入CO₂过程中，注入CO₂饱和度增大、压力增 大及生产井压力小条件下的压力变化导致的裂缝密度及孔隙度变化，再将孔隙度变化带入 Gassmann公式，加入混合流体饱和度变化参数后，计算CO₂注入前、后注入井及生产井的储 层弹性物性参数；

各向异性(HTI)介质模型是在各向同性介质模型的基础上构建的，最主要的区别 在于在第三步计算干岩石物性参数后，在流体替换过程中，引入Tod关于裂缝密度随压力变 化关系的理论，考虑流体替换过程中随着储层孔隙压力的变化而导致裂缝的张开和闭合。 图7描述流体替换过程中裂缝密度随孔隙压力和CO₂饱和度的变化情况，图8描述了流体替 换过程中裂缝密度和岩石孔隙度随孔隙压力的变化情况；

图7描述流体替换过程中裂缝密度随孔隙压力和CO₂饱和度的变化情况，其中孔隙 压力变化范围是7MPa～24MPa，CO₂饱和度变化范围是0％～100％，图中曲线为不同孔隙压 力和CO₂饱和度下储层裂缝密度的等值线图，从图中可看出在孔隙压力保持一定时，随着CO₂饱和度的变化，储层孔隙裂缝密度始终保持不变，而当CO₂饱和度一定时，随着孔隙压力增 加，储层裂缝密度不断增加，并且孔隙压力越大，在相同的压力梯度下裂缝密度变化越小。

图8(a)为裂缝密度随孔隙压力的变化曲线图；图8(b)为储层总孔隙度随孔隙压力 的变化图；图8(c)为储层总孔隙度随裂缝密度的变化关系图。图8(a)中随着孔隙压力不断 增加，储层裂缝密度不断变大，裂缝密度与孔隙压力呈一次函数关系变化；图8(b)中随着孔 隙压力不断增加，储层总孔隙度不断变大，储层总孔隙度与孔隙压力也呈一次函数变化；图 8(c)中储层总孔隙度与储层裂缝密度也呈一次函数变化。

Tod(2002)的公式表示为：

$Y_{cf}=Y_{c0}\{1-\sqrt{\frac{P_d-P_f}{P_d}}-\frac{1}{2}{e}^{-C_r(P_d-P_f)}*\sqrt{\frac{\pi }{C_r{P}_d}}*erf[\sqrt{C_r(P_d-P_f)}\left]\right\}---\left(8\right)$

其中，P_d是储层围压，P_f是孔隙压力，Y_c0是储层原始裂缝密度，π是圆周率，C_r是中间 变量，可表示为

其中，v_DRY是压力校正后干岩石的泊松比，μ_DRY是压力校正后干岩石的剪切模量，ar 为储层裂缝孔喉比；

式(8)所求Y_cf为压力变化过程中裂缝密度，换算成裂缝所代表的孔隙度fic_f：

${\mathrm{fic}}_f=\frac{4}{3}*\pi *ar*Y_{cf}---\left(10\right)$

重复运用Batzle-Wang公式和Gassman方程，求取流体替换过程随孔隙压力和流体 饱和度变化的混合流体体变模量，混合流体密度，进而求取饱和岩石的体变模量(见图9)、 剪切模量(见图10)，依照岩石物理参数之间的关系，求取饱和岩石的纵波速度、横波速度、 密度、体变模量、剪切模量，以及考虑裂缝变化后的岩石总孔隙度(见图8)，表4为流体替换 后不同孔隙压力和流体组成下物性的参数变化情况；

图9描述了饱和岩石的体变模量随储层孔隙压力和储层CO₂饱和度的变化情况。其 中，储层孔隙压力从7MPa～24MPa变化，CO₂饱和度从0％～100％变化，图中曲线为饱和岩石 的体变模量的等值线。从图中可看到，当CO₂饱和度一定时，随着储层孔隙压力的增加，饱和 岩石的体变模量逐渐减少，但减小幅度很小，相反的，当孔隙压力不变时，在CO₂饱和度较小 时(20％以下)，随着CO₂饱和度的增加，饱和岩石的体变模量迅速减小，当CO₂饱和度超过 40％后，随着CO₂饱和度的增加，饱和岩石的体变模量逐渐减小。

图10描述了饱和岩石的剪切模量随储层孔隙压力和储层CO₂饱和度的变化情况。 其中，储层孔隙压力从7MPa～24MPa变化，CO₂饱和度从0％～100％变化，图中曲线为饱和岩 石的剪切模量的等值线。从图中可看到，当孔隙压力一定时，随着储层CO₂饱和度的增加，饱 和岩石的剪切模量始终不发生变化，而当储层CO₂饱和度不变时，随着储层孔隙压力的增 加，饱和岩石的剪切模量不断减小，并且随着空隙压力越来越大，饱和岩石的剪切模量的变 化幅度越来越小。

第七，引入Thomson弱各向异性介质理论，结合第六步流体替换过程产生的各项岩 石物理参数，计算储层各向异性参数Epsilon、Gamma、Delta，同时，结合Thomson各向异性参 数，求取沿着裂缝方向的纵波速度和横波速度；

沿裂缝方向的纵波速度a0可表示为：

其中，Vp_sub为流体替换后储层的饱和岩石的纵波速度，im是Thomson各向异性参数 Epsilon与流体替换过程中考虑孔隙度变化后的泊松比v_f的函数，可表示为：

im＝2*Epsilon*(1-v_f)*(1-v_f)*(1-v_f)(12)

沿裂缝方向的横波速度b0可表示为：

式(12)中，im是Thomson各向异性参数Epsilon与流体替换过程中考虑孔隙度变化 后的泊松比v_f的函数；

式(13)中，Vs_sub为流体替换后储层的饱和岩石的横波速度，Gamma是Thomson各向 异性参数。

引入Rüger弱各向异性介质下方位AVO理论，结合流体替换后饱和岩石的纵波速 度、横波速度、密度，各向异性参数Epsilon、Gamma、Delta以及沿着裂缝方向的纵波速度和 横波速度，计算PS1波和PS2波方位AVO；

Rüger新的方位各向异性参数可表示为：

$Epsilon\_v=\frac{-1*Epsilon}{1+2*Epsilon}---\left(14\right)$

$Gamma\_v=\frac{-1*Gamma}{1+2*Gamma}---\left(15\right)$

$Delta\_v=\frac{Delta*Epsilon*(1+\frac{Epsilon}{f})}{(1+2*Epsilon)*(1+2*\frac{Epsilon}{f})}$其中，$f=1-{\left(\frac{b0}{a0}\right)}^2---\left(16\right)$

式(14)、(15)、(16)中Epsilon、Gamma、Delta为Thomson各向异性参数，a0、b0分别 为沿着裂缝方向的纵、横波速度；

依照Rüger方位AVO理论，PS1波和PS2波反射系数R_ps1和R_ps2可表示为：

R_ps1＝R_psa*cos(fry)+R_psb*sin(fry)(17)

R_ps2＝R_psa*sin(fry)+R_psb*cos(fry)(18)

式(17)和式(18)中R_psa和R_psb为横波沿着裂缝方向以及垂直裂缝方向震动时的反 射系数，fry为横波极化角；

$\begin{array}{c}{R}_{psa}=V_1*\sin \left(\frac{i*\pi }{180}\right)*\cos \left(\frac{i*\pi }{180}\right)+V_2*\frac{\sin \left(\frac{i*\pi }{180}\right)}{\cos \left(j\right)}+V_3*{\sin }^3\left(\frac{i*\pi }{180}\right)*\cos \left(\frac{i*\pi }{180}\right)+\\ V_4*\frac{{\sin }^3\left(\frac{i*\pi }{180}\right)}{\cos \left(j\right)}+V_5*\frac{{\sin }^3\left(\frac{i*\pi }{180}\right)}{\cos \left(j\right)}\end{array}---\left(19\right)$

$\begin{array}{c}{R}_{psb}=H_1*\sin \left(\frac{i*\pi }{180}\right)+H_2*\frac{\sin \left(\frac{i*\pi }{180}\right)*\cos \left(\frac{i*\pi }{180}\right)}{\cos \left(j\right)}+H_3*{\sin }^3\left(\frac{i*\pi }{180}\right)+H_4*\\ \frac{{\sin }^3\left(\frac{i*\pi }{180}\right)*\cos \left(\frac{i*\pi }{180}\right)}{\cos \left(j\right)}\end{array}---\left(19\right)$

$\left(20\right)$

$fry={\sin }^{-1}\left(\frac{-1*\sin \left(sita\right)}{\sqrt{1-\sin \left(j\right)*\sin \left(j\right)*\cos \left(sita\right)*\cos \left(sita\right)}}\right)---\left(21\right)$

其中，V₁、V₂、V₃、V₄、V₅、H₁、H₂、H₃、H₄均为Rüger方位各向异性参数以及反射面上下层 岩石的物性参数的函数；

运用式(19)和式(20)所求反射系数，绘制不同孔隙压力和流体饱和度下的PS1波 和PS2波AVO曲线，其结果如图11～图18；

图11～图14描述了在一定孔隙压力和CO₂饱和度下，PS1波和PS2波反射系数随入 射角和方位角的变化情况。图11描述了各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱 和度(SCO₂)下，随方位角变换的PS1转换波的AVA曲线图；其中图11(a)为孔隙压力为9MPa， CO₂饱和度为0％时，方位角分别为0°、30°、60°、90°时PS1波的AVA曲线图；图11(b)为孔隙压 力为15MPa，CO₂饱和度为20％时，方位角分别为0°、30°、60°、90°时PS1波的AVA曲线图；图11 (c)为孔隙压力为21MPa，CO₂饱和度为40％时，方位角分别为0°、30°、60°、90°时PS1波的AVA 曲线图。

图12描述了各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度(SCO₂)下，随入 射角和方位角变化时PS1转换波反射系数的等值线图；其中图12(a)为孔隙压力为9MPa，CO₂饱和度为0％时，PS1波反射系数随方位角和入射角的变化图；图12(b)为孔隙压力为15MPa， CO₂饱和度为20％时，PS1波反射系数随方位角和入射角的变化图；图12(c)为孔隙压力为 21MPa，CO₂饱和度为40％时，PS1波反射系数随方位角和入射角的变化图。

图13描述了各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度(SCO₂)下，随方 位角变换的PS2转换波的AVA曲线图；其中图13(a)为孔隙压力为9MPa，CO₂饱和度为0％时， 方位角分别为0°、30°、60°、90°时PS2波的AVA曲线图；图13(b)为孔隙压力为15MPa，CO₂饱和 度为20％时，方位角分别为0°、30°、60°、90°时PS2波的AVA曲线图；图13(c)为孔隙压力为 21MPa，CO₂饱和度为40％时，方位角分别为0°、30°、60°、90°时PS2波的AVA曲线图。

图14描述了各向异性(HTI)介质模型中，不同孔隙压力和CO₂饱和度(SCO₂)下，随入 射角和方位角变化时PS2转换波反射系数的等值线图；其中图14(a)为孔隙压力为9MPa，CO₂饱和度为0％时，PS2波反射系数随方位角和入射角的变化图；图14(b)为孔隙压力为15MPa， CO₂饱和度为20％时，PS2波反射系数随方位角和入射角的变化图；图14(c)为孔隙压力为 21MPa，CO₂饱和度为40％时，PS2波反射系数随方位角和入射角的变化图。

从图11～图14中可看出，对于不同压力和CO₂饱和度，PS1波和PS2波反射系数随入 射角和方位角的变化趋势大致相同，对于PS1波，在方位角为90°时，其反射系数恒为零；在 任意孔隙压力与CO₂饱和度下，当方位角一定时(不为90°)，PS1波反射系数随着入射角的增 加先减小再增大，PS1波反射系数的最小值所对应的入射角随孔隙压力与饱和度变化；当入 射角较小时，PS1波反射系数随方位角的增大而增加，当入射角增加到一定程度后，PS1波反 射系数随方位角增大而减小，同时，在入射角较小，且当在入射角一定时，PS1波反射系数随 方位角的增加而减小。对于PS2波，在方位角为0°时，其反射系数恒为零；在任意孔隙压力与 CO₂饱和度下，当方位角一定时(不为0°)，PS2波反射系数随着入射角的增加先增大再减小， PS2波反射系数的最大值所对应的入射角随孔隙压力与饱和度变化；当入射角较小时，PS2 波反射系数随方位角的增大而增加，当入射角增加到一定程度后，PS2波反射系数随方位角 增大而减小，同时，在入射角较小，且当在入射角一定时，PS2波反射系数随方位角的增加而 增加。

图15和图16，分别描述了当方位角为30°，纵波入射角为15°时，随着孔隙压力和 CO₂饱和度的变化，PS1波和PS2波反射系数的变化情况，为了描述方便，将PS1波和PS2波的 反射系数校正到正极性，从图中可以看到，PS1波和PS2波的反射系数的变化趋势近似相同， 但PS2波的反射系数整体小于PS1波，在孔隙压力保持不变时，随着CO₂饱和度的增加，PS1波 和PS2波的反射系数有着略微减小的趋势，但变化幅度很小，几乎为一常量；而当CO₂饱和度 一定时，随着孔隙压力的增加，PS1波和PS2的反射系数有着明显的增加。

图17是一定孔隙压力下，不同CO₂饱和度时PS1波和PS2波的AVA曲线。在CO₂饱和度 和孔隙压力一定时，PS1波和PS2波的反射系数均随着入射角的增加先增大再减小，并且在 相同的入射角达到最大值，在相同的入射角下，PS1波反射系数总是大于PS2波反射系数(均 为正极性时)，并且在二者达到最大值时，PS1波反射系数与PS2波反射系数之间的差异达到 最大；同时，从图中可以看到，当孔隙压力不变时，不同CO₂饱和度下PS1波和PS波的AVA曲线 基本重合，并且随着孔隙压力增大，不同CO₂饱和度下的PS1波和PS2波反射系数重合性越 好。图18是一定CO₂饱和度(SCO₂)下，不同孔隙压力时PS1波和PS2波的AVA曲线。描述了当方 位角为30°时，在CO₂饱和度(0％、20％、40％)一定时，不同孔隙压力(9MPa、15MPa、21MPa)下 的PS1波和PS2波的AVA曲线，同图17，为了描述方便，将PS1波与PS2波极性均调整为正极性。 横向比较图18中3幅图，在相同CO₂饱和度下和一定入射角时，PS1波与PS2波反射系数均随 着孔隙压力的增加而增加，在不同孔隙压力下PS1波和PS2波的反射系数相对差异较大；但 在不同CO₂饱和度下，不同压力之间的PS1波和PS2波反射系数的差异基本相同，PS1波与PS2 波反射系数变化趋势大致相同(见图17)。

第八，求取不同注入压力条件下、CO₂流体替换前后井模型PS转换波AVO反射系数， 并以褶积模型为基础分别制作井模型的各向同性介质和各向异性(HTI)介质转换波人工合 成地震记录。对于各向同性介质模型，将测井资料的深度域采样间隔如0.1米或0.125米看 作各向同性介质，从测井曲线起点，以Zoeppritz方程计算不同入射角度平面波转换波反射 系数，然后褶积获得合成地震记录。对于各向异性(HTI)介质，在储层顶面按照Rüger简化的 方位AVO计算公式计算不同方位的PS转换波AVO反射系数，而储层之外则仍然采样各向同性 介质方法计算PS转换波AVO反射系数，进而模拟不同注入CO₂阶段的转换波地震响应及四维 地震转换波差异地震响应，为准确解释CO₂地质封存中四维转换波信息奠定基础。

图19描述了各向同性介质模型中，CO₂饱和度为0％时，不同孔隙压力(9MPa、 15MPa、21MPa)下的入射角度道集人工合成地震记录，左侧三列分别为纵波、横波以及密度 测井曲线，中间三列为一定CO₂饱和度下不同孔隙压力时的转换波人工合成地震记录，右边 三列为孔隙压力下转换波人工合成地震记录的差异。从图19中可以看出，在相同的CO₂饱和 度下，在储层上部(2040ms以上)，不同孔隙压力下的转换波人工合成地震记录几乎没有差 异，不同射线参数下的差异都非常接近0，而从储层段开始，不同孔隙压力下的转换波人工 合成地震记录的差异开始出现，差异从转换波人工合成地震记录的储层部位就开始出现， 随着深度的增加，差异越来越明显，横向上看，孔隙压力变化幅度较大时(9MPa与21MPa之间 的差异)，差异相对更加明显，另外，在相同的压力梯度下，15MPa～21MPa之间的差异比9MPa ～15MPa之间的差异大。

图20描述了各向同性介质模型下，孔隙压力为15MPa时，不同CO₂饱和度下的入射 角度道集转换波合成地震记录及其差异，图左侧三列测井曲线分别为纵波、横波、密度曲 线，中间五列为孔隙压力为15MPa时，不同CO₂饱和度下的入射角度道集转换波人工合成地 震记录，右边四列为不同饱和度下转换波人工合成地震记录的差异。从图中可以看出，在相 同的储层压力下，在储层上部(2040ms以上)，不同CO₂饱和度下的转换波人工合成地震记录 几乎没有差异，不同射线参数下的差异都非常接近0，而从储层段开始，不同CO₂饱和度下的 转换波人工合成地震记录的差异开始出现，随着深度的增加，特别是强反射轴附近，差异越 来越明显，并且不同CO₂饱和度所对应的转换波人工合成地震记录之间的差异都相近。

结合图19和图20，可以看到，在各向同性介质的井模型中，孔隙压力的变化对转换 波人工合成地震记录的影响更大，是地震记录发生变化的主要因素。

在各向异性(HTI)介质模型AVO的基础上，我们以井模型和褶积模型为基础，构建 各向异性(HTI)介质条件下，PS1波和PS2波的人工合成地震记录，结果如图21和图22。图21 描述了各向异性(HTI)介质模型中，方位角为30°时，CO₂饱和度分别为0％时，不同孔隙压力 (9MPa、15MPa、21MPa)下PS1波的人工合成地震记录，图21左侧三列测井曲线分别为纵波、横 波、密度曲线，中间三列为不同孔隙压力和CO₂饱和度下的入射角度道集人工合成地震记 录，右边三列为不同饱和度下人工合成地震记录的差异。

从图21中可以看出，对于PS1波之间的差异，在相同的CO₂饱和度下，在储层上部 (2040ms以上)，不同孔隙压力下的入射角度道集人工合成地震记录几乎没有差异，不同射 线参数下的差异都非常接近0，而从储层段开始，不同孔隙压力下的人工合成地震记录的差 异开始出现，但差异都不明显，横向上看，孔隙压力变化幅度较大时(9MPa与21MPa之间的差 异)，差异相对更加明显，另外，在相同的压力梯度下，15MPa～21MPa之间的差异比9MPa～ 21MPa之间的差异大。

图22是各向异性(HTI)介质模型中，当方位角为30°，孔隙压力为15MPa时，不同CO₂饱和度下入射角度道集的PS2波人工合成地震记录及其差异，图左侧三列测井曲线分别为 纵波、横波、密度曲线，中间五列为孔隙压力为15MPa时，不同CO₂饱和度下的入射角度道集 人工合成地震记录，右边四列为不同饱和度下人工合成地震记录的差异。从图中可以看出， 在相同的储层压力下，在储层上部(即测井曲线1200米～1398米)，不同CO₂饱和度下的人工 合成地震记录几乎没有差异，不同射线参数下的差异都非常接近0，而从储层段开始，不同 CO₂饱和度下的人工合成地震记录的差异开始出现，而从储层段开始，不同CO₂饱和度下的人 工合成地震记录的差异开始出现，并且差异主要发生在合成地震记录中储层下部的强反射 层位，储层段的差异相对不太明显，并且不同CO₂饱和度所对应的人工合成地震记录之间的 差异都相近。图23是人工合成地震记录与实际转换波叠偏资料的对照结果，左边三列测井 曲线分别为纵波、横波、密度曲线，两列入射角度道集转换波人工合成地震记录分别为孔隙 压力为15MPa，CO₂饱和度为0％以及孔隙压力为15MPa，CO₂饱和度为10％的PS1波人工合成地 震记录，实际资料分别为1999年和2001年井附近的实际转换波叠加剖面地震资料，依照标 定结果，工区储层大致在2050ms左右，如图中所示Marly和Vuggy储层，转换波人工合成地震 记录与实际转换波叠偏资料显示了比较好的匹配，最右边两列分别为入射角度道集转换波 合成地震记录的差异与实际转换波资料的差异，合成地震记录的差异与图20中差异类似， 在储层上部几乎没有差异，在储层部位开始出现明显差异。对比两年转换波实际地震资料， 两年资料在剖面上均可划分出明显的储层，并且与合成地震记录能进行较好的标定，但因 两年实际资料的采集处理差异，两年实际转换波地震资料的差异稍显凌乱，但仍可依照井 资料做出标定。