CN106842326B - 无横波速度测井时砂泥互层co2地质封存时移地震正演模拟方法 - Google Patents

Abstract

本发明的无横波速度测井时砂泥互层CO₂地质封存时移地震正演模拟方法包括：采集数据；分层提取数据；利用测井资料计算砂岩矿物组成曲线，标定CO₂注入之后压力条件下的纵、横波速度，预测横波速度；进行混合流体替换注入之后的纵、横波速度曲线；保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射系数；与Ricker子波褶积，得到合成地震记录；利用注入CO₂前后及注入不同阶段、不同主频的时移正演人工合成地震记录相减，获得差异正演地震响应。本方法针对复杂的砂泥薄互层的特点，较好地解决了CO₂驱油与地质封存过程中不同储层特征的混合流体饱和度替换问题和预测出CO₂注入阶段压力变化的纵、横波速度。

Description

无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟 方法

技术领域

本发明涉及温室气体减排技术、CO₂地质封存的四维地震监测与CO₂驱油技术监测技术、油气田开发的地球物理监测技术领域。具体涉及一种利用Digby公式进行不同压力下纵、横波速度的预测，利用Gassmann理论进行混合流体替换和利用保留薄层信息的非均匀采样方法，进行井资料的深时转换，然后进行时移地震响应模型和差异地震模型制作的方法。

背景技术

在利用测井资料制作时移地震正演模型的过程中，需要用到纵波速度、横波速度以及密度。然而由于测量成本原因和早期技术的限制，多数情况下研究工区是没有横波速度测井资料的，特别是在CO₂地质封存区。因此很多的地球物理学者都在研究横波速度的预测方法，有采用经验公式，也有基于岩石物理理论的。

Wyllie等人在1958年和1963年陆续提出了充满盐水的孔隙介质的孔隙度与速度之间的经验关系：1/V＝(1-φ)/V_ma+φ/V_f1，其中V为岩石的整体速度，V_ma为岩石基质的速度，V_f1为孔隙流体的速度，φ为孔隙度。这个公式还通常可以写成层间旅行时的表达式：Δt＝(1-φ)Δt_ma+φΔt_f1，其中Δt代表整个岩层的旅行时，Δt_ma为基质的旅行时，而Δt_f1为孔隙流体的旅行时，而Wyllie的这个时间平均公式还包括许多假设和限制，如：这个方程要用于孔隙流体是盐水的情形，用于深度小于2700米的岩石，而且这个岩石的胶结程度和固结程度要很好，并且孔隙度为中等。Gardner于1974年给出了不同的岩性之间的速度与密度的关系，其中它的平均变换式为ρ＝0.23V^0.25，这个平均变换式是对所有岩性的速度与密度关系的最佳拟合，它适合于所有岩性。比较著名的经验公式包括Castagna等人在 1985年提出的著名的泥岩线为Vs＝0.862Vp-1.172。

随着岩石物理理论的完善，基于岩石物理理论的横波速度预测方法逐渐成为研究的主流。许多地球物理学家都喜欢用Gassmann方程进行横波速度的预测，这是因为Gassmann方程的大部分参数如颗粒的体变模量Kma与切变模量μ_ma等都是很好获得的，因此很多方法的给出无论是对于砂岩还是对碳酸盐都是基于 Gassmann方程的。不过在Gassmann方程中干岩石的弹性模量是个很难解决的问题，因此很多的地球物理学家给出了干岩石弹性模量的计算方法，知道了干岩石的弹性模量，岩石的纵横波速度就很好获得。Xu和White将Kuster和

于1974年建立的理论与差分有效介质理论结合，进行岩石弹性模量的计算，具体表现为利用孔隙纵横比来表征砂泥成分之间的关系。

无论用何种方法，在得到了岩石的纵、横波速度以及密度之后，需要进行反射系数的计算，然后利用褶积模型来计算地震响应。Zoeppritz方程给出的反射系数形式复杂、不易进行数值计算，并且物理意义不直观。许多学者对 Zoeppritz方程的解析进行了不同形式的简化。1955年，Koefoed研究发现泊松比的改变对P波(纵波)反射系数随入射角变化有显著影响。Koefoed的这个发现是在人们意识到流体对弹性参数有影响之后才变得有实用性。1961年 Bortfeld详细论述了平面纵波反射系数近似计算方法，在假设反射/透射界面两侧岩性参数变化较小时，给出了第一个区分流体和固体的简化公式；Richards& Frasier在1976年研究弹性波的散射问题时，在假设相邻地层介质弹性参数变化较小的情况下，从位移和应力连续角度对Zoeppritz方程进行了近似，给出了较为简单直观且精度较好的反射和透射系数的近似表达式。1985年，Shuey 在Aki和Richards近似公式基础上，给出了用不同角度项表示的突出泊松比的反射系数近似表达形式。Shuey公式成为目前地震正反演研究使用最多的公式，也是我们采用的反射系数计算公式之一。

而在进行CO₂地质封存的过程中，时移地震(或者四维地震)是监测CO₂是否安全封存于地下的一种必不可少的方法。而在进行时移地震资料解释的时候，需要利用测井资料制作时移人工合成地震记录来标定层位，并确定两次地震的差异特征，如监测地震数据(Monitor)与本底(Baseline)地震数据的差异，进而解释差异地震信息。而上述早期的方法无论是纵、横波速度的计算方法，还是近似反射系数的计算方法，都不适合CO₂注入阶段的特点。

CO₂捕集、利用与地质封存技术(Carbon Capture,Utlization and Storage简称CCUS)是目前国际上快速降低温室气体排放最有效的技术。CO₂地质封存中的时移(四维)地震监测技术是确定CO₂是否被安全地被封存在地下的关键技术。而开展砂泥薄互储层，在CO₂注入情况下的时移地震正演模拟，是监测CO₂地质封存不同阶段CO₂在地下长期封存安全性、泄漏风险以及确定CO₂驱油效果的关键。只有通过准确的模拟CO₂不同注入阶段的时移地震响应，才能有效建立正确的时移地震解释模型，为后续的时移地震定量解释奠定基础。

在目前国内、外开展的CO₂地质封存项目中，主要是利用老油田、废弃油田的已有油井开展CO₂的注入。这种方法与选择新的地区开展碳封存相比，新的地区选址需要钻探新井和开展大量的基础地质研究。因此在老油田开展碳封存，具有节省钻井、测井、地震采集等成本的显著优势，还可以利用CO₂驱油获得额外的收益，因此是国际上首选的CO₂地质封存方法，这种方法在国外的加拿大Weyburn油田、挪威Sleipner气田等及中国的胜利油田、陕西延长石油集团靖边油田、吉林油田等成功地开展。在老油田、废弃油田开展CO₂地质封存项目过程中，需要利用两次以上的地震资料对比，即注入CO₂前 (Baseline)，和注入CO₂后(Monitor)及注入后不同时期(Monitors)采集的地震资料，来确定地下CO₂的分布，确定CO₂在地下是否安全的被封存？同时可以检验CO₂注入地下后驱油的效果。

为了标定注入CO₂后不同阶段采集的地震监测数据(Monitor)，需要利用现有的井资料(注入CO₂前的测井数据)，制作CO₂注入后不同阶段的时移(四维)人工合成地震记录。或者说要制作随CO₂注入后，注入井点储层的压力增大、CO₂饱和度增大状况下的合成地震记录。此外，工业上为了获得最大的CO₂驱油效果，一般将生产井点的压力控制为较小的压力。这样在储层中会形成注入井点、油层内被CO₂驱替部分、油层未被CO₂波及(sweep)部分和生产井点等储层内部不同的CO₂饱和度、地层压力。但是测井工作一般是在油田勘探开发初期或者在CO₂注入前进行的，现有的技术还难以开展注入CO₂的二次测井。加上注入CO₂后是不宜进行钻井取心和测井的，一般油田早期的探井虽有岩心资料，但是往往距离CO₂地质封存的开发区域较远，而无法用探井的岩心来标定CO₂注入区内的测井资料。

这种情况下的CO₂地质封存区域的井资料存在下面几个问题：

(1)没有横波速度测井；

(2)没有岩心资料开展岩石物理测试而无法标定测井资料；

(3)测井所得到的纵波速度是在原始地层压力下，不是注入CO₂状态、生

产井压力状态下的纵波速度。

在CO₂注入过程中，随着CO₂的注入，储层的流体成分和饱和度和压力发生了变化。因此需要将注入前的纵、横波速度曲线校正成为含CO₂和随压力变化的纵、横波速度。即将注入前地层压力状态和没有CO₂影响的纵、横波测井曲线，校正为注入CO₂井和生产井压力状态和饱和度下的测井曲线。这样所计算得到的地震响应模型才与实际CO₂注入情况吻合，才能用来标定和解释注入CO₂后采集的地震监测资料。而在进行反射系数的计算过程中，可以采用Zoeppritz精确公式也可以采用Shuey等的近似公式计算反射系数，然后与不同频率、相位的Ricker 子波褶积，得到时移地震的地震响应。两次不同时移地震人工合成地震记录的差异，即为差异人工合成地震记录。

同时，在储层为砂泥薄互层的情况下，CO₂注入的层位为砂岩层，发生变化的为砂岩层，这其中既包括油层，也包括盐水层。因此在进行储层流体替换的纵、横波速度计算之前，还需根据测井和录井解释结果，将CO₂注入段内的薄含油、盐水砂岩层提取出来单独计算。这都是前述方法中未考虑到的问题。

发明内容

本发明的目的在于针对现有技术的不足,提供一种为CO₂地质封存的时移地震监测，提供可靠的地震正演模拟技术的无横波速度测井时砂泥互层CO₂地质封存时移地震正演模拟方法。

该方法是在无横波速度测井资料情况下，考虑了储层为砂泥薄互层情况下，如何进行多套薄储层、薄互层砂体的混合流体替换？如何进行不同压力和混合流体饱和度(CO₂、油、盐水)条件下的纵、横波速度曲线的计算？然后如何对测井资料进行保留薄层的深时转换，再制作时移地震响应和差异地震响应的方法。

本发明采用如下的技术方案：

本发明包括以下步骤：

1)采集数据:采集测井资料，包括随深度变化的岩石的孔隙度φ，体积密度ρ，纵波速度Vp_measured，油水饱和度以及岩石各个组分的含量。同时整理资料，确定各种流体的体变模量，各种岩石骨架的体变模量，岩石骨架的切变模量，差异压力P，配位数C_p'，变形之前接触区域的半径a与颗粒的半径R；

2)分层提取数据：根据详细的测井解释结果，将储层中的砂岩层和泥岩层分离。根据整理的测井资料将砂岩层的随深度变化的岩石的孔隙度φ，体积密度ρ，纵波速度Vp_measured，油水饱和度以及岩石各个组分的含量提取出来。

3)由于砂岩的组分基本上都是由石英，长石以及其它矿物组成，因此需要利用测井资料计算砂岩矿物组成，再利用Hill平均值法对矿物进行混合得到混合矿物的体变模量K_ma和切变模量μ_ma；利用Wood方程来计算砂岩层CO₂、油、盐水混合流体的体变模量K_f；

4)在CO₂注入区域没有岩心或者岩石物理资料的情况下，从整个研究区域或者CO₂注入区外围区域收集岩石物理测试资料，得到不同压力下的纵横波速度。根据取心的深度将其归位到相应的测井资料上，通过对比相同深度的岩心和测井资料，得到在相同测井深度下，岩性相同或者相近的岩心。利用得到岩心的纵、横波速度标定CO₂注入之后压力条件下的纵、横波速度。

5)预测横波速度：计算纵、横波速度的公式如下，

(1)

(2)

其中ρ为岩石的体积密度，K_sat为饱和岩石的切变模量，μ_dry为干岩石的切变模量，V_p为纵波速度，V_s为横波速度。

将含有加权系数W的配位数公式代入Digby公式，可以得到干岩石的体变模量K_dry和切变模量μ_dry的表达式，然后将Digby公式所表示的干岩石的体变模量和切变模量代入Gassmann方程得到饱和岩石体变模量K_sat，而饱和岩石的切变模量与干岩石的切变模量μ_dry相等，然后将饱和岩石的体变模量和干岩石的切变模量代入公式(1)，得到含有未知数W的纵波速度的表达式，作为预测纵波速度，减去实测纵波速度为0，得到一个只有未知数W的方程，公式如下

|Vp_measured-Vp_predicted(W)|→min (3)

在该公式中Vp_measured为实测的纵波速度，φ为孔隙度，ρ为岩石的体积密度，这些均可从测井资料获得；K_ma和μ_ma为骨架的体变模量和切变模量，K_f为混合流体的体变模量，均可可以从步骤2)获得；a/R，b/R为岩石变形前后接触区域的半径和颗粒半径的比值，前者从步骤1)获得，后者从Digby方程计算获得；W 为加权系数，为未知量；该方程只有W一个未知量，解方程得到W。将W代入 Digby方程可得到干岩石的切变模量μ_dry，然后代入公式(2)得到横波速度。

6)进行混合流体替换计算注入之后的纵、横波速度曲线：

依据测井资料获得的渗透特性资料，按照Gassmann方程，对不同砂岩储层进行不同层位、不同混合流体饱和度(CO₂、油和盐水)的流体替换。

根据岩石物理测试得到的注入CO₂之后压力条件下的纵波速度，利用公式(4) 得到其加权系数，然后将该岩心归位到该深度的测井资料中，那么所得到的加权系数W为该深度岩石在注入CO₂之后的加权系数。其它深度的加权系数，可利用公式(5)得到：

W_after＝W*W_core/W_log (5)

其中W_after为注入CO₂之后的配位数，W_core为计算得到岩心的配位数，W_log为岩心所对应深度岩石的配位数；得到注入CO₂之后岩石的配位数之后，将其代入 Digby方程，得到注入CO₂之后压力下干岩石的体变模量和切变模量；然后利用 Wood方程计算注入CO₂之后混合流体的体变模量，然后将注入CO₂之后压力下干岩石的体变模量和切变模量以及流体的体变模量代入Gassmann方程得到饱和岩石的体变模量，然后饱和岩石的体变模量和干岩石的切变模量代入公式(1)和(2)得到注入CO₂之后的纵横波速度。

7)保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射系数：

在得到注入CO₂前后压力饱和度下的纵、横波速度之后，首先对得到的数据进行深时转换，将纵、横波速度以及密度从深度域转换到时间域，采用非均匀采样，或者采用更高采样率，将深度域薄层信息全部转换到时间域。这样，时间域的纵、横波速度以及密度曲线将保留薄互储层与泥岩夹层的全部信息，而不会漏掉薄储层与夹层信息，然后再利用精确的Zoeppritz方程或者其近似公式(如 Shuey近似公式)计算平面波反射系数序列；

8)与Ricker子波褶积，得到合成地震记录：在得到时间域井资料的反射系数序列后，将其与不同主频的Ricker子波进行褶积，得到不同子波主频、垂向上不同层位具有不同混合流体(CO₂、油和盐水)饱和度及地层压力变化条件下的人工合成地震记录；

9)差异正演地震响应获得：利用注入CO₂前后及注入不同阶段、不同主频的时移正演人工合成地震记录相减，获得差异正演地震响应。差异地震响应可用于研究CO₂地质封存时移地震监测资料的其他属性信息，可以用来标定时移地震资料的层位，是时移地震解释的基础。

上述方案还包括：

所述步骤3)中，Hill平均值法为：

其中，M_v是Voigt上限,M_R是Reuss下限,f_i and M_i分别是第i中组分的体积分数以及模量(切变模量或者体变模量)；

所述步骤3)以及步骤6)中的Wood方程为：

K_f＝1/(∑(S_i/k_fi)) (i＝1,2,3...) (9)

其中K_f为混合流体的体变模量，S_i和k_fi分别为第i种流体的饱和度和体变模量。

所述步骤5)中，配位数公式为：

C_p＝W(11.759e^1-φ-12.748) (10)

其中W为加权系数，φ为孔隙度。

所述步骤5)、6)中，Digby公式为：

其中b可以表示为,

而d满足公式(4)，

介质的泊松比可以用如下公式计算

其中K_dry与μ_dry分别为干岩石的体变模量和切变模量；v与μ_ma分别为岩石颗粒的泊松比与切变模量；φ为孔隙度；C_p为配位数；P为差异压力；α为变形之前接触区域的半径，b为变形之后接触区域的半径，R为颗粒的半径。ν_x为某种介质的泊松比，比如骨架的泊松比，干岩石的泊松比等等，K_x和μ_x为某种介质的体变模量和切变模量，如果求的是骨架的泊松比，那么K_x和μ_x则为骨架的体变模量和切变模量。

所述步骤5)、6)中，Gassmann方程为：

μ_sat＝μ_dry (17)

其中K_sat和μ_sat分别为饱和岩石的体变模量和切变模量，K_dry与μ_dry分别为干岩石的体变模量和切变模量，K_ma为骨架的体变模量，φ为孔隙度，K_f为流体的体变模量。

所述步骤7)中，Zoeppritz方程为：

其中R_pp、R_ps、T_pp和分T_ps别是纵波反射系数，横波反射系数，纵波透射系数，横波透射系数；α₁、β₁、ρ₁、α₂、β₂和ρ₂分别为分界面两侧纵、横波速度及介质密度。

所述步骤7)中，Ricker子波公式为：

g(t)＝[1-(2πft)²]exp[-(πft)²] (19)

其中f为地震子波主频；t为时间。

本发明的有益效果是：本发明考虑到在储层为多套薄砂层、砂泥薄互层的情况下，发生流体替换的层位为砂岩层，首先将砂岩层提取出来，根据各套储层的渗透特征，进行混合流体(CO₂、油、盐水)替换的计算。即储层的渗透特性不同，注入CO₂进行流体替换的多少不同；含油砂岩进行的是CO₂、油和盐水三相混合流体替换，而储层段内的不含油砂岩，即盐水层，进行的是CO₂和盐水两相混合流体替换。其次，本发明考虑到了在CO₂注入过程中，压力发生了变化，提出了随压力以及混合流体饱和度变化的纵、横波速度的计算方法。利用注入CO₂前的测井曲线，预测注入CO₂后饱和度、压力变化下的纵、横速度曲线。这样，使得到的纵、横波速度更加符合实际，才能正确计算注入CO₂后的人工合成地震记录。为了多套薄层、薄互层砂岩信息从测井资料的深度域转换到时间域，采用了非均匀采样方法和高采样率的方法，进行测井资料的深时转换。而在进行时移地震响应模型的制作过程中，使用Zoeppritz方程或其近似公式(如Shuey公式)进行反射系数的计算。这样时移地震响应正演模型同时考虑了CO₂注入层的压力以及饱和度的变化及薄互层的特点，可以有效进行时移地震资料的正演模拟。

附图说明

图1胜利油田某CO₂注入区，预测的纵波速度随压力、CO₂饱和度变化图。其中，淡实线为差异压力14.6MPa，黑线为差异压力20MPa。实线、虚线和点线依次代表CO₂饱和度为0、20％和40％.

图2利用胜利油田某CO₂注入区附近实际采集的纵、横波速度测井资料，将储层的砂岩层提取出来进行原始地层压力下横波速度预测的结果，其中实线为实测横波速度，虚线为预测横波速度。黑点为横波速度预测误差。误差很小，说明预测效果的可靠性，和方法可以用于预测相邻地质区域的横波速度。

图3对图2砂岩储层原始地层压力下预测的横波速度，进行不同CO₂注入压力的横波速度预测结果。

图4原始地层压力下，实际采集胜利油田某CO₂注入区的深度域纵、横波速度、密度、反射系数序列反射系数，与采用保留薄层信息的深时转换方法得到的时间域的反射系数。其中REI代表反射系数序列。注意，在时间域测井曲线及反射系数序列中里，采用非均匀采样和高采样率的深时转换后，深度域的3234- 3436m,3290-3291m,3311-薄砂层，全部保留在时间域170ms，190ms,以及210ms 位置。即时间域反射系数保留了更多的薄储层信息。

图5是CO₂注入之前压力下的不同饱和度的地震响应，每个图中从左至右依次为纵波速度、横波速度、密度、人工合成地震记录。Angle为平面地震波入射角度，differ为差异地震记录。(b)(c)(d)中最后一栏的differ图，为该孔隙压力饱和度下的人工合成地震记录，与(a)也就是与原始未注入CO₂的情况下合成地震记录的差异，即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa，CO₂饱和度为0％的合成地震记录；(b)为利用预测横波速度计算的差异压力为20MPa，CO₂饱和度为0％的合成地震记录；(c)为利用预测纵、横波速度计算的差异压力为20MPa，CO₂饱和度为10％的合成地震记录；(d)为利用预测纵、横波速度计算的差异压力为20MPa，CO₂饱和度为40％的合成地震记录。

图6是CO₂注入之后的压力下的不同饱和度的地震响应。(b)(c)(d)差异压力为14.6MPa下的合成地震记录，每个图中从左至右依次为纵波速度、横波速度、密度以及人工合成地震记录。Angle为平面地震波入射角度，differ为差异地震记录。(b)(c)(d)图中最后一栏的differ图，为该孔隙压力饱和度下的人工合成地震记录，与(a)也就是与原始未注入CO₂的情况下的差异，即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa，CO₂饱和度为0％的合成地震记录；(b)为利用预测纵、横波速度计算的差异压力为14.6MPa，CO₂饱和度为0％的合成地震记录；(c)为利用预测纵、横波速度计算的差异压力为14.6MPa，CO₂饱和度为10％的合成地震记录；(d)为利用预测纵、横波速度计算的差异压力为14.6MPa，CO₂饱和度为40％的合成地震记录。

具体实施方式

下面结合附图对本发明做详细描述；本发明公式中的*表示乘号。

本发明针对复杂储层即砂泥岩互层的储层，提出一种CO₂地质封存中，时移地震正演模型的制作方法。随着的CO₂注入，储层中的流体发生变化，那么储层的压力以及流体饱和也发生了变化，需要特别注意的是，发生这些变化的层位为注入的砂岩层。因此，本发明的做法首先需要将砂岩层提取出来，然后进行流体替换以及压力变化的纵横波速度的预测。得到不同压力饱和度下的纵、横波速度以及密度之后，对其进行深时转换，利用精确的Zoeppritz方程进行反射系数的计算，然后将得到的反射系数与相应频率的Ricker褶积，得到时移地震响应。得到结果更加符合实际情况。

本发明的步骤如下：

1)采集数据:采集测井资料，包括随深度变化的岩石的孔隙度φ，体积密度ρ，纵波速度Vp_measured，油水的饱和度以及岩石各个组分的含量。同时整理资料，确定各种流体的体变模量，各种岩石骨架的体变模量，岩石骨架的切变模量，差异压力P，配位数C_p'，变形之前接触区域的半径a与颗粒的半径R；

2)分层提取数据：根据详细的测井解释结果，将储层中的砂岩层和泥岩层分离。根据整理的测井资料将砂岩层的随深度变化的岩石的孔隙度φ，体积密度ρ，纵波速度Vp_measured，油水饱和度以及岩石各个组分的含量提取出来。这里我们以胜利油田某CO₂地质封存区为例，提取得到的砂岩层数据如下表，

同时，由于上述砂岩层岩石的组分基本上都是由砂岩、泥岩以及灰岩组成，因此需要利用Hill平均值法对矿物进行混合得到混合矿物的体变模量K_ma和切变模量μ_ma；Hill平均值法为：

其中，M_v是Voigt上限,M_R是Reuss下限,f_i and M_i分别是第i中组分的体积分数以及模量(切变模量或者体变模量)；

利用Wood方程来计算砂岩层混合流体的体变模量K_f，需要注意的是，CO₂注入之前的流体为油和水两相，CO₂注入之后的混合流体则为油、水以及CO₂三相； Wood方程为：

K_f＝1/(∑(S_i/k_fi)) (i＝1,2,3...) (9)

其中K_f为混合流体的体变模量，S_i和k_fi分别为第i种流体的饱和度和体变模量。

3)在整个区域或者外围区域收集岩石物理测试资料，得到不同压力下的纵、横波速度。根据取心的深度将其归位到测井资料上，通过对比相同深度的岩心和测井资料，得到在相同测井深度下，岩性相同或者相近的岩心。利用得到岩心的纵、横波速度标定CO₂注入之后压力条件下的纵横波速度。

4)预测横波速度：计算纵横波速度的公式如下，

其中ρ为岩石的体积密度，K_sat为饱和岩石的切变模量，μ_dry为干岩石的切变模量，V_p为纵波速度，V_s为横波速度。

将含有加权系数W的配位数公式代入Digby公式，可以得到干岩石的体变模量K_dry和切变模量μ_dry的表达式，

其中配位数的表达式为：

C_p＝W(11.759e^1-φ-12.748) (10)

其中W为加权系数，φ为孔隙度。

Digby公式为：

其中b可以表示为,

而d满足公式(4)，

介质的泊松比可以用如下公式计算

其中K_dry与μ_dry分别为干岩石的体变模量和切变模量；v与μ_ma分别为岩石颗粒的泊松比与切变模量；φ为孔隙度；C_p为配位数；P为差异压力；α为变形之前接触区域的半径，b为变形之后接触区域的半径，R为颗粒的半径。ν_x为某种介质的泊松比，比如骨架的泊松比，干岩石的泊松比等等，K_x和μ_x为某种介质的体变模量和切变模量，如果求的是骨架的泊松比，那么K_x和μ_x则为骨架的体变模量和切变模量。

5)然后将Digby公式所表示的干岩石的体变模量和切变模量代入 Gassmann方程得到饱和岩石体变模量K_sat，而饱和岩石的切变模量与干岩石的切变模量μ_dry相等，Gassmann方程为：

μ_sat＝μ_dry (17)

其中K_sat和μ_sat分别为饱和岩石的体变模量和切变模量，K_dry与μ_dry分别为干岩石的体变模量和切变模量，K_ma为骨架的体变模量，φ为孔隙度，K_f为流体的体变模量。

然后将饱和岩石的体变模量和干岩石的切变模量代入公式(1)，得到含有未知数W的纵波速度的表达式，作为预测纵波速度，减去实测纵波速度为0，得到一个只有未知数W的方程，公式如下

|Vp_measured-Vp_predicted(W)|→min (3)

在该公式中Vp_measured为实测的纵波速度，φ为孔隙度，ρ为岩石的体积密度，这些均可从测井资料获得；K_ma和μ_ma为骨架的体变模量和切变模量，K_f为混合流体的体变模量，均可可以从步骤2)获得；a/R，b/R为岩石变形前后接触区域的半径和颗粒半径的比值，前者从步骤1)获得，后者从Digby方程计算获得；W 为加权系数，为未知量；该方程只有W一个未知量，解方程得到W。将W代入 Digby方程可得到干岩石的切变模量μ_dry，然后代入公式(2)得到横波速度。得到的横波速度的预测结果如图1。其中实线为实测横波速度，虚线为预测横波速度。

6)进行流体替换计算注入之后的纵、横波速度：

依据测井资料获得的渗透特性资料，按照Gassmann方程，对不同砂岩储层进行不同混合流体饱和度的(CO₂、油、盐水)的流体替换。

根据岩石物理测试得到的注入CO₂之后压力条件下的纵波速度，利用公式(4) 得到其加权系数，然后将该岩心归位到该深度的测井资料中，那么所得到的加权系数W为该深度岩石在注入CO₂之后的加权系数。其它深度的加权系数，可利用公式(5)得到：

W_after＝W*W_core/W_log (5)

其中W_after为注入CO₂之后的配位数，W_core为计算得到岩心的配位数，W_log为岩心所对应深度岩石的配位数；得到注入CO₂之后岩石的配位数之后，将其代入 Digby方程，得到注入CO₂之后压力下干岩石的体变模量和切变模量；然后利用 Wood方程计算注入CO₂之后混合流体的体变模量，然后将注入CO₂之后压力下干岩石的体变模量和切变模量以及流体的体变模量代入Gassmann方程得到饱和岩石的体变模量，然后饱和岩石的体变模量和干岩石的切变模量代入公式(1)和 (2)得到注入CO₂之后的纵横波速度。计算得到的不同压力下的横波速度见图 2。

7)保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射系数：

在得到注入CO₂前后压力饱和度下的纵、横波速度之后，首先对得到的数据进行深时转换，将纵、横波速度以及密度从深度域转换到时间域，采用非均匀采样，或者采用更高采样率，将深度域薄层信息全部转换到时间域。这样，时间域的纵、横波速度以及密度曲线将保留薄互储层与泥岩夹层的全部信息，然后在利用利用精确的Zoeppritz方程或者其近似公式计算平面波反射系数序列；其中Zoeppritz方程为：

其中R_pp、R_ps、T_pp和分T_ps别是纵波反射系数，横波反射系数，纵波透射系数，横波透射系数；α₁、β₁、ρ₁、α₂、β₂和ρ₂分别为分界面两侧纵、横波速度及介质密度。

8)与Ricker子波褶积，得到合成地震记录：在得到时间域井资料的反射系数序列后，将其与不同主频的Ricker子波进行褶积，得到不同子波主频、垂向上不同层位具有不同混合流体(CO₂、油和盐水)饱和度及地层压力变化条件下的人工合成地震记录；

9)差异正演地震响应获得：利用注入CO₂前后及注入不同阶段、不同主频的人工合成地震记录相减，获得差异正演地震响应，并用于研究CO₂地质封存时移地震监测资料的其他属性信息。可以用来标定时移地震资料的层位。

其中Ricker子波公式为：

g(t)＝[1-(2πft)²]exp[-(πft)²]

时移人工合成地震记录与差异地震记录如图5，图6。其中图5为CO₂注入之前的压力下的不同饱和度的地震响应。每个图中从左至右依次为纵波速度、横波速度、密度、人工合成地震记录。angle为平面地震波入射角度，differ为差异地震记录。(b)(c)(d)图中最后一栏differ图为该孔隙压力饱和度下的人工合成地震记录，与(a)也就是与原始未注入CO₂的情况下的差异，即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa，CO₂饱和度为0％的合成地震记录；(b)为利用预测横波速度计算的差异压力为20MPa，CO₂饱和度为0％的合成地震记录；(c)为利用预测纵横波速度计算的差异压力为20MPa，CO2饱和度为10％的合成地震记录；(d)为利用预测纵横波速度计算的差异压力为20MPa， CO₂饱和度为40％的合成地震记录。

图6为CO₂注入之后的压力下的不同饱和度的地震响应。其中(b)(c)(d) 为差异压力为14.6MPa下的合成地震记录，每个图中从左至右依次为纵波速度、横波速度、密度以及合成地震记录，angle为平面地震波入射角度，differ为差异地震记录。(b)(c)(d)图中最后一栏的differ图为该孔隙压力饱和度下的人工合成地震记录，与(a)也就是与原始未注入CO₂的情况下的差异，即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa，CO₂饱和度为0％的合成地震记录；(b)为利用预测纵、横波速度计算的差异压力为14.6MPa，CO₂饱和度为0％的合成地震记录；(c)为利用预测纵、横波速度计算的差异压力为14.6MPa，CO₂饱和度为10％的合成地震记录；(d)为利用预测纵、横波速度计算的差异压力为14.6MPa，CO₂饱和度为40％的合成地震记录。

从图中我们可以看到，随着CO₂的注入，也就是差异压力的减小，振幅差异越来越明显，而在相同的压力下，随着CO₂饱和度的增大，振幅差异也越来越明显。说明本发明能够很好的模拟发生流体替换前后的地震响应与差异地震响应。

Claims (4)

1.无横波速度测井时砂泥互层CO₂地质封存时移地震正演模拟方法，其特征包括：

1)采集数据:采集测井资料包括随深度变化的岩石的孔隙度φ，体积密度ρ，实测纵波速度Vp_measured，岩石各个组分的含量以及油水的饱和度，确定各种流体的体变模量，各种岩石骨架的体变模量，岩石骨架的切变模量，差异压力p，配位数C_p，变形之前接触区域的半径a与颗粒的半径R；

2)分层提取数据：根据测井解释结果，将储层中的砂岩层和泥岩层分离，将砂岩层的随深度变化的岩石的孔隙度φ，体积密度ρ，实测纵波速度Vp_measured，油水饱和度以及岩石各个组分的含量提取出来；

3)利用测井资料计算砂岩矿物组成曲线，再利用Hill平均值法对矿物进行混合得到混合矿物的体变模量和切变模量；利用Wood方程来计算砂岩层CO₂、油、盐水混合流体的体变模量K_f；

4)利用CO₂注入区域或者CO₂注入区外围区域收集岩石物理测试资料，得到不同压力下的纵、横波速度，根据取心的深度将其归位到相应的测井资料上，通过对比相同深度的岩心和测井资料，得到在相同测井深度下，岩性相同或者相近的岩心，利用得到岩心的纵、横波速度标定CO₂注入之后压力条件下的纵、横波速度；

5)预测横波速度：计算纵、横波速度的公式如下，

其中ρ为体积密度，K_sat为饱和岩石的体变模量，μ_dry为干岩石的切变模量，V_p为纵波速度，V_s为横波速度；

将含有加权系数W的配位数公式代入Digby公式，可以得到干岩石的体变模量K_dry和切变模量μ_dry的表达式，然后将Digby公式所表示的干岩石的体变模量和切变模量代入Gassmann方程得到饱和岩石体变模量K_sat，而饱和岩石的切变模量与干岩石的切变模量μ_dry相等，然后将饱和岩石的体变模量和干岩石的切变模量代入公式(1)，得到含有未知数W的纵波速度的表达式，作为预测纵波速度Vp_predicted，减去实测纵波速度为0，得到一个只有未知数W的方程，公式如下：

|Vp_measured-Vp_predicted(W)|→0 (3)

在该公式中Vp_measured为实测的纵波速度，W为加权系数，e为自然常数，φ为孔隙度，μ_ma为岩石骨架的切变模量，b为岩石变形后接触区域的半径，π为圆周率，ν为岩石骨架的泊松比，R为颗粒的半径，a为岩石变形前接触区域的半径，K_ma为岩石骨架的体变模量，ρ为体积密度；其中Vp_measured，φ，ν以及ρ均可从测井资料获得；K_ma和μ_ma为骨架的体变模量和切变模量，K_f为混合流体的体变模量，均可以从步骤2)获得；a/R，b/R为岩石变形前后接触区域的半径和颗粒半径的比值，前者从步骤1)获得，后者从Digby方程计算获得；W为加权系数，为未知量；

该方程只有W一个未知量，解方程得到W，将W代入Digby方程可得到干岩石的切变模量μ_dry，然后代入公式(2)得到横波速度；

6)进行混合流体替换计算注入之后的纵、横波速度曲线：

依据测井资料获得的渗透特性资料，按照Gassmann方程，对不同砂岩储层进行不同混合流体饱和度CO₂、油和盐水的流体替换；

根据岩石物理测试得到的注入CO₂之后压力条件下的纵波速度，利用公式(4)得到其加权系数，然后将该岩心归位到该深度的测井资料中，那么所得到的加权系数W为该深度岩石在注入CO₂之后的加权系数；其它深度的加权系数，可利用公式(5)得到：

W_after＝W*W_core/W_log (5)

其中W_after为注入CO₂之后的加权系数，W_core为计算得到岩心的加权系数，W_log为岩心所对应深度岩石的加权系数；得到注入CO₂之后岩石的加权系数之后，将其代入Digby方程，得到注入CO₂之后压力下干岩石的体变模量和切变模量；然后利用Wood方程计算注入CO₂之后混合流体的体变模量，然后将注入CO₂之后压力下干岩石的体变模量和切变模量以及流体的体变模量代入Gassmann方程得到饱和岩石的体变模量，然后饱和岩石的体变模量和干岩石的切变模量代入公式(1)和(2)得到注入CO₂之后的纵横波速度；

7)保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射数：

在得到注入CO₂前后压力饱和度下的纵、横波速度之后，首先对得到的数据进行深时转换，将纵、横波速度以及密度从深度域转换到时间域，采用非均匀采样或者采用更高采样率，将深度域薄层信息全部转换到时间域，时间域的纵、横波速度以及密度曲线将保留薄互储层与泥岩夹层的全部信息，然后再利用精确的Zoeppritz方程或者其近似公式计算平面波反射系数序列；

8)与Ricker子波褶积，得到合成地震记录：

在得到时间域井资料的反射系数序列后，将其与不同主频的Ricker子波进行褶积，得到不同子波主频、垂向上不同层位具有不同混合流体饱和度及地层压力变化条件下的人工合成地震记录；

9)差异正演地震响应获得：利用注入CO₂前后、不同主频的时移正演人工合成地震记录相减，获得差异正演地震响应。

2.根据权利要求1所述的无横波速度测井时砂泥互层CO₂地质封存时移地震正演模拟方法，其特征在于，所述步骤5)中，配位数公式为：

C_p＝W(11.759e^1-φ-12.748) (10)

其中W为加权系数，φ为孔隙度；

所述步骤5)、6)中，Digby公式为：

其中b可以表示为,

而d满足公式(14)，

介质的泊松比可以用如下公式计算

其中K_dry与μ_dry分别为干岩石的体变模量和切变模量；ν与μ_ma分别为岩石骨架的泊松比与切变模量；φ为孔隙度；C_p为配位数；p为差异压力；a为变形之前接触区域的半径，b为变形之后接触区域的半径，R为颗粒的半径；ν_x为某种介质的泊松比，K_x和μ_x为某种介质的体变模量和切变模量。

3.根据权利要求1或2所述的无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟方法，其特征在于，所述步骤2)中，Hill平均值法为：

其中，M_v是Voigt上限,M_R是Reuss下限,f_i和M_i分别是第i中组分的体积分数以及切变模量或者体变模量；

所述步骤3)以及步骤6)中的Wood方程为：

K_f＝1/(∑(S_j/k_j)) j＝1,2,3... (9)

其中K_f为混合流体的体变模量，S_j和K_fj分别为第j种流体的饱和度和体变模量；所述步骤5)、6)中，Gassmann方程为：

μ_sat＝μ_dry (17)

其中K_sat和μ_sat分别为饱和岩石的体变模量和切变模量，K_dry与μ_dry分别为干岩石的体变模量和切变模量，K_ma为骨架的体变模量，φ为孔隙度，K_f为混合流体的体变模量；

所述步骤7)中，Zoeppritz方程为：

其中R_pp、R_ps、T_pp和分T_ps别是纵波反射系数，横波反射系数，纵波透射系数，横波透射系数；α₁、β₁、ρ₁、α₂、β₂和ρ₂分别为分界面两侧纵、横波速度及介质密度。

4.根据权利要求1所述的无横波速度测井时砂泥互层CO₂地质封存时移地震正演模拟方法，其特征在于，所述步骤8)中，Ricker子波公式为：

g(t)＝[1-(2πft)²]exp[-(πft)²] (19)

其中f为地震子波主频；t为时间。

Images