CN106842326A - 无横波速度测井时砂泥互层co2地质封存时移地震正演模拟方法 - Google Patents
无横波速度测井时砂泥互层co2地质封存时移地震正演模拟方法 Download PDFInfo
- Publication number
- CN106842326A CN106842326A CN201510886011.9A CN201510886011A CN106842326A CN 106842326 A CN106842326 A CN 106842326A CN 201510886011 A CN201510886011 A CN 201510886011A CN 106842326 A CN106842326 A CN 106842326A
- Authority
- CN
- China
- Prior art keywords
- rock
- modulus
- shear
- wave velocity
- dry
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明的无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟方法包括:采集数据;分层提取数据;利用测井资料计算砂岩矿物组成曲线,标定CO2注入之后压力条件下的纵、横波速度,预测横波速度;进行混合流体替换注入之后的纵、横波速度曲线;保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射系数;与Ricker子波褶积,得到合成地震记录;利用注入CO2前后及注入不同阶段、不同主频的时移正演人工合成地震记录相减,获得差异正演地震响应。本方法针对复杂的砂泥薄互层的特点,较好地解决了CO2驱油与地质封存过程中不同储层特征的混合流体饱和度替换问题和预测出CO2注入阶段压力变化的纵、横波速度。
Description
技术领域
本发明涉及温室气体减排技术、CO2地质封存的四维地震监测与CO2驱油技术监测技术、油气田开发的地球物理监测技术领域。具体涉及一种利用Digby公式进行不同压力下纵、横波速度的预测,利用Gassmann理论进行混合流体替换和利用保留薄层信息的非均匀采样方法,进行井资料的深时转换,然后进行时移地震响应模型和差异地震模型制作的方法。
背景技术
在利用测井资料制作时移地震正演模型的过程中,需要用到纵波速度、横波速度以及密度。然而由于测量成本原因和早期技术的限制,多数情况下研究工区是没有横波速度测井资料的,特别是在CO2地质封存区。因此很多的地球物理学者都在研究横波速度的预测方法,有采用经验公式,也有基于岩石物理理论的。
Wyllie等人在1958年和1963年陆续提出了充满盐水的孔隙介质的孔隙度与速度之间的经验关系:1/V=(1-φ)/Vma+φ/Vf1,其中V为岩石的整体速度,Vma为岩石基质的速度,Vf1为孔隙流体的速度,φ为孔隙度。这个公式还通常可以写成层间旅行时的表达式:Δt=(1-φ)Δtma+φΔtf1,其中Δt代表整个岩层的旅行时,Δtma为基质的旅行时,而Δtf1为孔隙流体的旅行时,而Wyllie的这个时间平均公式还包括许多假设和限制,如:这个方程要用于孔隙流体是盐水的情形,用于深度小于2700米的岩石,而且这个岩石的胶结程度和固结程度要很好,并且孔隙度为中等。Gardner于1974年给出了不同的岩性之间的速度与密度的关系,其中它的平均变换式为ρ=0.23V0.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公式成为目前地震正反演研究使用最多的公式,也是我们采用的反射系数计算公式之一。
而在进行CO2地质封存的过程中,时移地震(或者四维地震)是监测CO2是否安全封存于地下的一种必不可少的方法。而在进行时移地震资料解释的时候,需要利用测井资料制作时移人工合成地震记录来标定层位,并确定两次地震的差异特征,如监测地震数据(Monitor)与本底(Baseline)地震数据的差异,进而解释差异地震信息。而上述早期的方法无论是纵、横波速度的计算方法,还是近似反射系数的计算方法,都不适合CO2注入阶段的特点。
CO2捕集、利用与地质封存技术(Carbon Capture,Utlization and Storage简称CCUS)是目前国际上快速降低温室气体排放最有效的技术。CO2地质封存中的时移(四维)地震监测技术是确定CO2是否被安全地被封存在地下的关键技术。而开展砂泥薄互储层,在CO2注入情况下的时移地震正演模拟,是监测CO2地质封存不同阶段CO2在地下长期封存安全性、泄漏风险以及确定CO2驱油效果的关键。只有通过准确的模拟CO2不同注入阶段的时移地震响应,才能有效建立正确的时移地震解释模型,为后续的时移地震定量解释奠定基础。
在目前国内、外开展的CO2地质封存项目中,主要是利用老油田、废弃油田的已有油井开展CO2的注入。这种方法与选择新的地区开展碳封存相比,新的地区选址需要钻探新井和开展大量的基础地质研究。因此在老油田开展碳封存,具有节省钻井、测井、地震采集等成本的显著优势,还可以利用CO2驱油获得额外的收益,因此是国际上首选的CO2地质封存方法,这种方法在国外的加拿大Weyburn油田、挪威Sleipner气田等及中国的胜利油田、陕西延长石油集团靖边油田、吉林油田等成功地开展。在老油田、废弃油田开展CO2地质封存项目过程中,需要利用两次以上的地震资料对比,即注入CO2前(Baseline),和注入CO2后(Monitor)及注入后不同时期(Monitors)采集的地震资料,来确定地下CO2的分布,确定CO2在地下是否安全的被封存?同时可以检验CO2注入地下后驱油的效果。
为了标定注入CO2后不同阶段采集的地震监测数据(Monitor),需要利用现有的井资料(注入CO2前的测井数据),制作CO2注入后不同阶段的时移(四维)人工合成地震记录。或者说要制作随CO2注入后,注入井点储层的压力增大、CO2饱和度增大状况下的合成地震记录。此外,工业上为了获得最大的CO2驱油效果,一般将生产井点的压力控制为较小的压力。这样在储层中会形成注入井点、油层内被CO2驱替部分、油层未被CO2波及(sweep)部分和生产井点等储层内部不同的CO2饱和度、地层压力。但是测井工作一般是在油田勘探开发初期或者在CO2注入前进行的,现有的技术还难以开展注入CO2的二次测井。加上注入CO2后是不宜进行钻井取心和测井的,一般油田早期的探井虽有岩心资料,但是往往距离CO2地质封存的开发区域较远,而无法用探井的岩心来标定CO2注入区内的测井资料。
这种情况下的CO2地质封存区域的井资料存在下面几个问题:
(1)没有横波速度测井;
(2)没有岩心资料开展岩石物理测试而无法标定测井资料;
(3)测井所得到的纵波速度是在原始地层压力下,不是注入CO2状态、生
产井压力状态下的纵波速度。
在CO2注入过程中,随着CO2的注入,储层的流体成分和饱和度和压力发生了变化。因此需要将注入前的纵、横波速度曲线校正成为含CO2和随压力变化的纵、横波速度。即将注入前地层压力状态和没有CO2影响的纵、横波测井曲线,校正为注入CO2井和生产井压力状态和饱和度下的测井曲线。这样所计算得到的地震响应模型才与实际CO2注入情况吻合,才能用来标定和解释注入CO2后采集的地震监测资料。而在进行反射系数的计算过程中,可以采用Zoeppritz精确公式也可以采用Shuey等的近似公式计算反射系数,然后与不同频率、相位的Ricker子波褶积,得到时移地震的地震响应。两次不同时移地震人工合成地震记录的差异,即为差异人工合成地震记录。
同时,在储层为砂泥薄互层的情况下,CO2注入的层位为砂岩层,发生变化的为砂岩层,这其中既包括油层,也包括盐水层。因此在进行储层流体替换的纵、横波速度计算之前,还需根据测井和录井解释结果,将CO2注入段内的薄含油、盐水砂岩层提取出来单独计算。这都是前述方法中未考虑到的问题。
发明内容
本发明的目的在于针对现有技术的不足,提供一种为CO2地质封存的时移地震监测,提供可靠的地震正演模拟技术的无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟方法。
该方法是在无横波速度测井资料情况下,考虑了储层为砂泥薄互层情况下,如何进行多套薄储层、薄互层砂体的混合流体替换?如何进行不同压力和混合流体饱和度(CO2、油、盐水)条件下的纵、横波速度曲线的计算?然后如何对测井资料进行保留薄层的深时转换,再制作时移地震响应和差异地震响应的方法。
本发明采用如下的技术方案:
本发明包括以下步骤:
1)采集数据:采集测井资料,包括随深度变化的岩石的孔隙度φ,体积密度ρ,纵波速度Vpmeasured,油水饱和度以及岩石各个组分的含量。同时整理资料,确定各种流体的体变模量,各种岩石骨架的体变模量,岩石骨架的切变模量,差异压力P,配位数Cp',变形之前接触区域的半径a与颗粒的半径R;
2)分层提取数据:根据详细的测井解释结果,将储层中的砂岩层和泥岩层分离。根据整理的测井资料将砂岩层的随深度变化的岩石的孔隙度φ,体积密度ρ,纵波速度Vpmeasured,油水饱和度以及岩石各个组分的含量提取出来。
3)由于砂岩的组分基本上都是由石英,长石以及其它矿物组成,因此需要利用测井资料计算砂岩矿物组成,再利用Hill平均值法对矿物进行混合得到混合矿物的体变模量Kma和切变模量μma;利用Wood方程来计算砂岩层CO2、油、盐水混合流体的体变模量Kf;
4)在CO2注入区域没有岩心或者岩石物理资料的情况下,从整个研究区域或者CO2注入区外围区域收集岩石物理测试资料,得到不同压力下的纵横波速度。根据取心的深度将其归位到相应的测井资料上,通过对比相同深度的岩心和测井资料,得到在相同测井深度下,岩性相同或者相近的岩心。利用得到岩心的纵、横波速度标定CO2注入之后压力条件下的纵、横波速度。
5)预测横波速度:计算纵、横波速度的公式如下,
(1)
(2)
其中ρ为岩石的体积密度,Ksat为饱和岩石的切变模量,μdry为干岩石的切变模量,Vp为纵波速度,Vs为横波速度。
将含有加权系数W的配位数公式代入Digby公式,可以得到干岩石的体变模量Kdry和切变模量μdry的表达式,然后将Digby公式所表示的干岩石的体变模量和切变模量代入Gassmann方程得到饱和岩石体变模量Ksat,而饱和岩石的切变模量与干岩石的切变模量μdry相等,然后将饱和岩石的体变模量和干岩石的切变模量代入公式(1),得到含有未知数W的纵波速度的表达式,作为预测纵波速度,减去实测纵波速度为0,得到一个只有未知数W的方程,公式如下
|Vpmeasured-Vppredicted(W)|→min (3)
在该公式中Vpmeasured为实测的纵波速度,φ为孔隙度,ρ为岩石的体积密度,这些均可从测井资料获得;Kma和μma为骨架的体变模量和切变模量,Kf为混合流体的体变模量,均可可以从步骤2)获得;a/R,b/R为岩石变形前后接触区域的半径和颗粒半径的比值,前者从步骤1)获得,后者从Digby方程计算获得;W为加权系数,为未知量;该方程只有W一个未知量,解方程得到W。将W代入Digby方程可得到干岩石的切变模量μdry,然后代入公式(2)得到横波速度。
6)进行混合流体替换计算注入之后的纵、横波速度曲线:
依据测井资料获得的渗透特性资料,按照Gassmann方程,对不同砂岩储层进行不同层位、不同混合流体饱和度(CO2、油和盐水)的流体替换。
根据岩石物理测试得到的注入CO2之后压力条件下的纵波速度,利用公式(4)得到其加权系数,然后将该岩心归位到该深度的测井资料中,那么所得到的加权系数W为该深度岩石在注入CO2之后的加权系数。其它深度的加权系数,可利用公式(5)得到:
Wafter=W*Wcore/Wlog (5)
其中Wafter为注入CO2之后的配位数,Wcore为计算得到岩心的配位数,Wlog为岩心所对应深度岩石的配位数;得到注入CO2之后岩石的配位数之后,将其代入Digby方程,得到注入CO2之后压力下干岩石的体变模量和切变模量;然后利用Wood方程计算注入CO2之后混合流体的体变模量,然后将注入CO2之后压力下干岩石的体变模量和切变模量以及流体的体变模量代入Gassmann方程得到饱和岩石的体变模量,然后饱和岩石的体变模量和干岩石的切变模量代入公式(1)和(2)得到注入CO2之后的纵横波速度。
7)保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射系数:
在得到注入CO2前后压力饱和度下的纵、横波速度之后,首先对得到的数据进行深时转换,将纵、横波速度以及密度从深度域转换到时间域,采用非均匀采样,或者采用更高采样率,将深度域薄层信息全部转换到时间域。这样,时间域的纵、横波速度以及密度曲线将保留薄互储层与泥岩夹层的全部信息,而不会漏掉薄储层与夹层信息,然后再利用精确的Zoeppritz方程或者其近似公式(如Shuey近似公式)计算平面波反射系数序列;
8)与Ricker子波褶积,得到合成地震记录:在得到时间域井资料的反射系数序列后,将其与不同主频的Ricker子波进行褶积,得到不同子波主频、垂向上不同层位具有不同混合流体(CO2、油和盐水)饱和度及地层压力变化条件下的人工合成地震记录;
9)差异正演地震响应获得:利用注入CO2前后及注入不同阶段、不同主频的时移正演人工合成地震记录相减,获得差异正演地震响应。差异地震响应可用于研究CO2地质封存时移地震监测资料的其他属性信息,可以用来标定时移地震资料的层位,是时移地震解释的基础。
上述方案还包括:
所述步骤3)中,Hill平均值法为:
其中,Mv是Voigt上限,MR是Reuss下限,fi and Mi分别是第i中组分的体积分数以及模量(切变模量或者体变模量);
所述步骤3)以及步骤6)中的Wood方程为:
Kf=1/(∑(Si/kfi)) (i=1,2,3...) (9)
其中Kf为混合流体的体变模量,Si和kfi分别为第i种流体的饱和度和体变模量。
所述步骤5)中,配位数公式为:
Cp=W(11.759e1-φ-12.748) (10)
其中W为加权系数,φ为孔隙度。
所述步骤5)、6)中,Digby公式为:
其中b可以表示为,
而d满足公式(4),
介质的泊松比可以用如下公式计算
其中Kdry与μdry分别为干岩石的体变模量和切变模量;v与μma分别为岩石颗粒的泊松比与切变模量;φ为孔隙度;Cp为配位数;P为差异压力;α为变形之前接触区域的半径,b为变形之后接触区域的半径,R为颗粒的半径。νx为某种介质的泊松比,比如骨架的泊松比,干岩石的泊松比等等,Kx和μx为某种介质的体变模量和切变模量,如果求的是骨架的泊松比,那么Kx和μx则为骨架的体变模量和切变模量。
所述步骤5)、6)中,Gassmann方程为:
μsat=μdry (17)
其中Ksat和μsat分别为饱和岩石的体变模量和切变模量,Kdry与μdry分别为干岩石的体变模量和切变模量,Kma为骨架的体变模量,φ为孔隙度,Kf为流体的体变模量。
所述步骤7)中,Zoeppritz方程为:
其中Rpp、Rps、Tpp和分Tps别是纵波反射系数,横波反射系数,纵波透射系数,横波透射系数;α1、β1、ρ1、α2、β2和ρ2分别为分界面两侧纵、横波速度及介质密度。
所述步骤7)中,Ricker子波公式为:
g(t)=[1-(2πft)2]exp[-(πft)2] (19)
其中f为地震子波主频;t为时间。
本发明的有益效果是:本发明考虑到在储层为多套薄砂层、砂泥薄互层的情况下,发生流体替换的层位为砂岩层,首先将砂岩层提取出来,根据各套储层的渗透特征,进行混合流体(CO2、油、盐水)替换的计算。即储层的渗透特性不同,注入CO2进行流体替换的多少不同;含油砂岩进行的是CO2、油和盐水三相混合流体替换,而储层段内的不含油砂岩,即盐水层,进行的是CO2和盐水两相混合流体替换。其次,本发明考虑到了在CO2注入过程中,压力发生了变化,提出了随压力以及混合流体饱和度变化的纵、横波速度的计算方法。利用注入CO2前的测井曲线,预测注入CO2后饱和度、压力变化下的纵、横速度曲线。这样,使得到的纵、横波速度更加符合实际,才能正确计算注入CO2后的人工合成地震记录。为了多套薄层、薄互层砂岩信息从测井资料的深度域转换到时间域,采用了非均匀采样方法和高采样率的方法,进行测井资料的深时转换。而在进行时移地震响应模型的制作过程中,使用Zoeppritz方程或其近似公式(如Shuey公式)进行反射系数的计算。这样时移地震响应正演模型同时考虑了CO2注入层的压力以及饱和度的变化及薄互层的特点,可以有效进行时移地震资料的正演模拟。
附图说明
图1胜利油田某CO2注入区,预测的纵波速度随压力、CO2饱和度变化图。其中,淡实线为差异压力14.6MPa,黑线为差异压力20MPa。实线、虚线和点线依次代表CO2饱和度为0、20%和40%.
图2利用胜利油田某CO2注入区附近实际采集的纵、横波速度测井资料,将储层的砂岩层提取出来进行原始地层压力下横波速度预测的结果,其中实线为实测横波速度,虚线为预测横波速度。黑点为横波速度预测误差。误差很小,说明预测效果的可靠性,和方法可以用于预测相邻地质区域的横波速度。
图3对图2砂岩储层原始地层压力下预测的横波速度,进行不同CO2注入压力的横波速度预测结果。
图4原始地层压力下,实际采集胜利油田某CO2注入区的深度域纵、横波速度、密度、反射系数序列反射系数,与采用保留薄层信息的深时转换方法得到的时间域的反射系数。其中REI代表反射系数序列。注意,在时间域测井曲线及反射系数序列中里,采用非均匀采样和高采样率的深时转换后,深度域的3234-3436m,3290-3291m,3311-薄砂层,全部保留在时间域170ms,190ms,以及210ms位置。即时间域反射系数保留了更多的薄储层信息。
图5是CO2注入之前压力下的不同饱和度的地震响应,每个图中从左至右依次为纵波速度、横波速度、密度、人工合成地震记录。Angle为平面地震波入射角度,differ为差异地震记录。(b)(c)(d)中最后一栏的differ图,为该孔隙压力饱和度下的人工合成地震记录,与(a)也就是与原始未注入CO2的情况下合成地震记录的差异,即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa,CO2饱和度为0%的合成地震记录;(b)为利用预测横波速度计算的差异压力为20MPa,CO2饱和度为0%的合成地震记录;(c)为利用预测纵、横波速度计算的差异压力为20MPa,CO2饱和度为10%的合成地震记录;(d)为利用预测纵、横波速度计算的差异压力为20MPa,CO2饱和度为40%的合成地震记录。
图6是CO2注入之后的压力下的不同饱和度的地震响应。(b)(c)(d)差异压力为14.6MPa下的合成地震记录,每个图中从左至右依次为纵波速度、横波速度、密度以及人工合成地震记录。Angle为平面地震波入射角度,differ为差异地震记录。(b)(c)(d)图中最后一栏的differ图,为该孔隙压力饱和度下的人工合成地震记录,与(a)也就是与原始未注入CO2的情况下的差异,即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa,CO2饱和度为0%的合成地震记录;(b)为利用预测纵、横波速度计算的差异压力为14.6MPa,CO2饱和度为0%的合成地震记录;(c)为利用预测纵、横波速度计算的差异压力为14.6MPa,CO2饱和度为10%的合成地震记录;(d)为利用预测纵、横波速度计算的差异压力为14.6MPa,CO2饱和度为40%的合成地震记录。
具体实施方式
下面结合附图对本发明做详细描述;本发明公式中的*表示乘号。
本发明针对复杂储层即砂泥岩互层的储层,提出一种CO2地质封存中,时移地震正演模型的制作方法。随着的CO2注入,储层中的流体发生变化,那么储层的压力以及流体饱和也发生了变化,需要特别注意的是,发生这些变化的层位为注入的砂岩层。因此,本发明的做法首先需要将砂岩层提取出来,然后进行流体替换以及压力变化的纵横波速度的预测。得到不同压力饱和度下的纵、横波速度以及密度之后,对其进行深时转换,利用精确的Zoeppritz方程进行反射系数的计算,然后将得到的反射系数与相应频率的Ricker褶积,得到时移地震响应。得到结果更加符合实际情况。
本发明的步骤如下:
1)采集数据:采集测井资料,包括随深度变化的岩石的孔隙度φ,体积密度ρ,纵波速度Vpmeasured,油水的饱和度以及岩石各个组分的含量。同时整理资料,确定各种流体的体变模量,各种岩石骨架的体变模量,岩石骨架的切变模量,差异压力P,配位数Cp',变形之前接触区域的半径a与颗粒的半径R;
2)分层提取数据:根据详细的测井解释结果,将储层中的砂岩层和泥岩层分离。根据整理的测井资料将砂岩层的随深度变化的岩石的孔隙度φ,体积密度ρ,纵波速度Vpmeasured,油水饱和度以及岩石各个组分的含量提取出来。这里我们以胜利油田某CO2地质封存区为例,提取得到的砂岩层数据如下表,
同时,由于上述砂岩层岩石的组分基本上都是由砂岩、泥岩以及灰岩组成,因此需要利用Hill平均值法对矿物进行混合得到混合矿物的体变模量Kma和切变模量μma;Hill平均值法为:
其中,Mv是Voigt上限,MR是Reuss下限,fi and Mi分别是第i中组分的体积分数以及模量(切变模量或者体变模量);
利用Wood方程来计算砂岩层混合流体的体变模量Kf,需要注意的是,CO2注入之前的流体为油和水两相,CO2注入之后的混合流体则为油、水以及CO2三相;Wood方程为:
Kf=1/(∑(Si/kfi)) (i=1,2,3...) (9)
其中Kf为混合流体的体变模量,Si和kfi分别为第i种流体的饱和度和体变模量。
3)在整个区域或者外围区域收集岩石物理测试资料,得到不同压力下的纵、横波速度。根据取心的深度将其归位到测井资料上,通过对比相同深度的岩心和测井资料,得到在相同测井深度下,岩性相同或者相近的岩心。利用得到岩心的纵、横波速度标定CO2注入之后压力条件下的纵横波速度。
4)预测横波速度:计算纵横波速度的公式如下,
其中ρ为岩石的体积密度,Ksat为饱和岩石的切变模量,μdry为干岩石的切变模量,Vp为纵波速度,Vs为横波速度。
将含有加权系数W的配位数公式代入Digby公式,可以得到干岩石的体变模量Kdry和切变模量μdry的表达式,
其中配位数的表达式为:
Cp=W(11.759e1-φ-12.748) (10)
其中W为加权系数,φ为孔隙度。
Digby公式为:
其中b可以表示为,
而d满足公式(4),
介质的泊松比可以用如下公式计算
其中Kdry与μdry分别为干岩石的体变模量和切变模量;v与μma分别为岩石颗粒的泊松比与切变模量;φ为孔隙度;Cp为配位数;P为差异压力;α为变形之前接触区域的半径,b为变形之后接触区域的半径,R为颗粒的半径。νx为某种介质的泊松比,比如骨架的泊松比,干岩石的泊松比等等,Kx和μx为某种介质的体变模量和切变模量,如果求的是骨架的泊松比,那么Kx和μx则为骨架的体变模量和切变模量。
5)然后将Digby公式所表示的干岩石的体变模量和切变模量代入Gassmann方程得到饱和岩石体变模量Ksat,而饱和岩石的切变模量与干岩石的切变模量μdry相等,Gassmann方程为:
μsat=μdry (17)
其中Ksat和μsat分别为饱和岩石的体变模量和切变模量,Kdry与μdry分别为干岩石的体变模量和切变模量,Kma为骨架的体变模量,φ为孔隙度,Kf为流体的体变模量。
然后将饱和岩石的体变模量和干岩石的切变模量代入公式(1),得到含有未知数W的纵波速度的表达式,作为预测纵波速度,减去实测纵波速度为0,得到一个只有未知数W的方程,公式如下
|Vpmeasured-Vppredicted(W)|→min (3)
在该公式中Vpmeasured为实测的纵波速度,φ为孔隙度,ρ为岩石的体积密度,这些均可从测井资料获得;Kma和μma为骨架的体变模量和切变模量,Kf为混合流体的体变模量,均可可以从步骤2)获得;a/R,b/R为岩石变形前后接触区域的半径和颗粒半径的比值,前者从步骤1)获得,后者从Digby方程计算获得;W为加权系数,为未知量;该方程只有W一个未知量,解方程得到W。将W代入Digby方程可得到干岩石的切变模量μdry,然后代入公式(2)得到横波速度。得到的横波速度的预测结果如图1。其中实线为实测横波速度,虚线为预测横波速度。
6)进行流体替换计算注入之后的纵、横波速度:
依据测井资料获得的渗透特性资料,按照Gassmann方程,对不同砂岩储层进行不同混合流体饱和度的(CO2、油、盐水)的流体替换。
根据岩石物理测试得到的注入CO2之后压力条件下的纵波速度,利用公式(4)得到其加权系数,然后将该岩心归位到该深度的测井资料中,那么所得到的加权系数W为该深度岩石在注入CO2之后的加权系数。其它深度的加权系数,可利用公式(5)得到:
Wafter=W*Wcore/Wlog (5)
其中Wafter为注入CO2之后的配位数,Wcore为计算得到岩心的配位数,Wlog为岩心所对应深度岩石的配位数;得到注入CO2之后岩石的配位数之后,将其代入Digby方程,得到注入CO2之后压力下干岩石的体变模量和切变模量;然后利用Wood方程计算注入CO2之后混合流体的体变模量,然后将注入CO2之后压力下干岩石的体变模量和切变模量以及流体的体变模量代入Gassmann方程得到饱和岩石的体变模量,然后饱和岩石的体变模量和干岩石的切变模量代入公式(1)和(2)得到注入CO2之后的纵横波速度。计算得到的不同压力下的横波速度见图2。
7)保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射系数:
在得到注入CO2前后压力饱和度下的纵、横波速度之后,首先对得到的数据进行深时转换,将纵、横波速度以及密度从深度域转换到时间域,采用非均匀采样,或者采用更高采样率,将深度域薄层信息全部转换到时间域。这样,时间域的纵、横波速度以及密度曲线将保留薄互储层与泥岩夹层的全部信息,然后在利用利用精确的Zoeppritz方程或者其近似公式计算平面波反射系数序列;其中Zoeppritz方程为:
其中Rpp、Rps、Tpp和分Tps别是纵波反射系数,横波反射系数,纵波透射系数,横波透射系数;α1、β1、ρ1、α2、β2和ρ2分别为分界面两侧纵、横波速度及介质密度。
8)与Ricker子波褶积,得到合成地震记录:在得到时间域井资料的反射系数序列后,将其与不同主频的Ricker子波进行褶积,得到不同子波主频、垂向上不同层位具有不同混合流体(CO2、油和盐水)饱和度及地层压力变化条件下的人工合成地震记录;
9)差异正演地震响应获得:利用注入CO2前后及注入不同阶段、不同主频的人工合成地震记录相减,获得差异正演地震响应,并用于研究CO2地质封存时移地震监测资料的其他属性信息。可以用来标定时移地震资料的层位。
其中Ricker子波公式为:
g(t)=[1-(2πft)2]exp[-(πft)2]
时移人工合成地震记录与差异地震记录如图5,图6。其中图5为CO2注入之前的压力下的不同饱和度的地震响应。每个图中从左至右依次为纵波速度、横波速度、密度、人工合成地震记录。angle为平面地震波入射角度,differ为差异地震记录。(b)(c)(d)图中最后一栏differ图为该孔隙压力饱和度下的人工合成地震记录,与(a)也就是与原始未注入CO2的情况下的差异,即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa,CO2饱和度为0%的合成地震记录;(b)为利用预测横波速度计算的差异压力为20MPa,CO2饱和度为0%的合成地震记录;(c)为利用预测纵横波速度计算的差异压力为20MPa,CO2饱和度为10%的合成地震记录;(d)为利用预测纵横波速度计算的差异压力为20MPa,CO2饱和度为40%的合成地震记录。
图6为CO2注入之后的压力下的不同饱和度的地震响应。其中(b)(c)(d)为差异压力为14.6MPa下的合成地震记录,每个图中从左至右依次为纵波速度、横波速度、密度以及合成地震记录,angle为平面地震波入射角度,differ为差异地震记录。(b)(c)(d)图中最后一栏的differ图为该孔隙压力饱和度下的人工合成地震记录,与(a)也就是与原始未注入CO2的情况下的差异,即差异地震记录。其中(a)为利用实际测井资料即差异压力为20MPa,CO2饱和度为0%的合成地震记录;(b)为利用预测纵、横波速度计算的差异压力为14.6MPa,CO2饱和度为0%的合成地震记录;(c)为利用预测纵、横波速度计算的差异压力为14.6MPa,CO2饱和度为10%的合成地震记录;(d)为利用预测纵、横波速度计算的差异压力为14.6MPa,CO2饱和度为40%的合成地震记录。
从图中我们可以看到,随着CO2的注入,也就是差异压力的减小,振幅差异越来越明显,而在相同的压力下,随着CO2饱和度的增大,振幅差异也越来越明显。说明本发明能够很好的模拟发生流体替换前后的地震响应与差异地震响应。
Claims (4)
1.无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟方法,其特征包括:
1)采集数据:采集测井资料包括随深度变化的岩石的孔隙度φ,体积密度ρ,纵波速度Vpmeasured,岩石各个组分的含量以及油水的饱和度,确定各种流体的体变模量,各种岩石骨架的体变模量,岩石骨架的切变模量,差异压力P,配位数Cp',变形之前接触区域的半径a与颗粒的半径R;
2)分层提取数据:根据测井解释结果,将储层中的砂岩层和泥岩层分离,将砂岩层的随深度变化的岩石的孔隙度φ,体积密度ρ,纵波速度Vpmeasured,油水饱和度以及岩石各个组分的含量提取出来;
3)利用测井资料计算砂岩矿物组成曲线,再利用Hill平均值法对矿物进行混合得到混合矿物的体变模量Kma和切变模量μma;利用Wood方程来计算砂岩层CO2、油、盐水混合流体的体变模量Kf;
4)利用CO2注入区域或者CO2注入区外围区域收集岩石物理测试资料,得到不同压力下的纵横波速度,根据取心的深度将其归位到相应的测井资料上,通过对比相同深度的岩心和测井资料,得到在相同测井深度下,岩性相同或者相近的岩心;利用得到岩心的纵、横波速度标定CO2注入之后压力条件下的纵、横波速度;
5)预测横波速度:计算纵、横波速度的公式如下,
其中ρ为岩石的体积密度,Ksat为饱和岩石的切变模量,μdry为干岩石的切变模量,Vp为纵波速度,Vs为横波速度;
将含有加权系数W的配位数公式代入Digby公式,得到干岩石的体变模量Kdry和切变模量μdry的表达式,然后将Digby公式所表示的干岩石的体变模量和切变模量代入Gassmann方程得到饱和岩石体变模量Ksat,而饱和岩石的切变模量与干岩石的切变模量μdry相等,然后将饱和岩石的体变模量和干岩石的切变模量代入公式(1),得到含有未知数W的纵波速度的表达式,作为预测纵波速度,减去实测纵波速度为0,得到一个只有未知数W的方程,公式如下
|Vpmeasured-Vppredicted(W)|→min (3)
在该公式中Vpmeasured为实测的纵波速度,φ为孔隙度,ρ为岩石的体积密度;Kma和μma为骨架的体变模量和切变模量,Kf为混合流体的体变模量;a/R,b/R为岩石变形前后接触区域的半径和颗粒半径的比值,后者从Digby方程计算获得;W为加权系数,为未知量;
该方程只有W一个未知量,解方程得到W,将W代入Digby方程得到干岩石的切变模量μdry,然后代入公式(2)得到横波速度;
6)进行混合流体替换计算注入之后的纵、横波速度曲线:
依据测井资料获得的渗透特性资料,按照Gassmann方程,对不同砂岩储层进行包括CO2、油和盐水的不同混合流体饱和度的流体替换;
根据岩石物理测试得到的注入CO2之后压力条件下的纵波速度,利用公式(4)得到其加权系数,然后将该岩心归位到该深度的测井资料中,那么所得到的加权系数W为该深度岩石在注入CO2之后的加权系数;其它深度的加权系数,可利用公式(5)得到:
Wafter=W*Wcore/Wlog (5)
其中Wafter为注入CO2之后的配位数,Wcore为计算得到岩心的配位数,Wlog为岩心所对应深度岩石的配位数;得到注入CO2之后岩石的配位数之后,将其代入Digby方程,得到注入CO2之后压力下干岩石的体变模量和切变模量;利用Wood方程计算注入CO2之后混合流体的体变模量,将注入CO2之后压力下干岩石的体变模量和切变模量以及流体的体变模量代入Gassmann方程得到饱和岩石的体变模量,然后饱和岩石的体变模量和干岩石的切变模量代入公式(1)和(2)得到注入CO2之后的纵横波速度;
7)保留薄层信息的测井资料时深转换以及利用Zoeppritz方程计算反射系数:
在得到注入CO2前后压力饱和度下的纵、横波速度之后,首先对得到的数据进行深时转换,将纵、横波速度以及密度从深度域转换到时间域,采用非均匀采样或者采用更高采样率,将深度域薄层信息全部转换到时间域;时间域的纵、横波速度以及密度曲线将保留薄互储层与泥岩夹层的全部信息,然后再利用精确的Zoeppritz方程或者其近似公式计算平面波反射系数序列;
8)与Ricker子波褶积,得到合成地震记录:
在得到时间域井资料的反射系数序列后,将其与不同主频的Ricker子波进行褶积,得到不同子波主频、垂向上不同层位具有不同混合流体饱和度及地层压力变化条件下的人工合成地震记录;
9)差异正演地震响应获得:利用注入CO2前后及注入不同阶段、不同主频的时移正演人工合成地震记录相减,获得差异正演地震响应。
2.根据权利要求1所述的无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟方法,其特征在于,所述步骤5)中,配位数公式为:
Cp=W(11.759e1-φ-12.748) (10)
其中W为加权系数,φ为孔隙度。
所述步骤5)、6)中,Digby公式为:
其中b可以表示为,
而d满足公式(4),
介质的泊松比可以用如下公式计算
其中Kdry与μdry分别为干岩石的体变模量和切变模量;v与μma分别为岩石颗粒的泊松比与切变模量;φ为孔隙度;Cp为配位数;P为差异压力;α为变形之前接触区域的半径,b为变形之后接触区域的半径,R为颗粒的半径。νx为某种介质的泊松比,比如骨架的泊松比,干岩石的泊松比等等,Kx和μx为某种介质的体变模量和切变模量,如果求的是骨架的泊松比,那么Kx和μx则为骨架的体变模量和切变模量。
3.根据权利要求1或2所述的无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟方法,其特征在于,所述步骤3)中,Hill平均值法为:
其中,Mv是Voigt上限,MR是Reuss下限,fi和Mi分别是第i中组分的体积分数以及切变模量或者体变模量;
所述步骤3)以及步骤6)中的Wood方程为:
Kf=1/(∑(Si/kfi)) (i=1,2,3...) (9)
其中Kf为混合流体的体变模量,Si和kfi分别为第i种流体的饱和度和体变模量;所述步骤5)、6)中,Gassmann方程为:
μsat=μdry (17)
其中Ksat和μsat分别为饱和岩石的体变模量和切变模量,Kdry与μdry分别为干岩石的体变模量和切变模量,Kma为骨架的体变模量,φ为孔隙度,Kf为流体的体变模量;
所述步骤7)中,Zoeppritz方程为:
其中Rpp、Rps、Tpp和分Tps别是纵波反射系数、横波反射系数、纵波透射系数、横波透射系数;α1、β1、ρ1、α2、β2和ρ2分别为分界面两侧纵、横波速度及介质密度。
4.根据权利要求3所述的无横波速度测井时砂泥互层CO2地质封存时移地震正演模拟方法,其特征在于,所述步骤8)中,Ricker子波公式为:
g(t)=[1-(2πft)2]exp[-(πft)2] (19)
其中f为地震子波主频;t为时间。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510886011.9A CN106842326B (zh) | 2015-12-04 | 2015-12-04 | 无横波速度测井时砂泥互层co2地质封存时移地震正演模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510886011.9A CN106842326B (zh) | 2015-12-04 | 2015-12-04 | 无横波速度测井时砂泥互层co2地质封存时移地震正演模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106842326A true CN106842326A (zh) | 2017-06-13 |
CN106842326B CN106842326B (zh) | 2020-10-13 |
Family
ID=59150555
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510886011.9A Active CN106842326B (zh) | 2015-12-04 | 2015-12-04 | 无横波速度测井时砂泥互层co2地质封存时移地震正演模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106842326B (zh) |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110441817A (zh) * | 2019-06-27 | 2019-11-12 | 西北大学 | 孔隙介质中基于4d pp-ps波反演ccus压力和饱和度变化的方法 |
CN111025396A (zh) * | 2020-01-06 | 2020-04-17 | 中国石油化工股份有限公司 | 基于人工智能算法的油藏物性参数地震预测方法 |
CN111323817A (zh) * | 2020-04-15 | 2020-06-23 | 中国矿业大学(北京) | 基于深度学习的二氧化碳封存监测方法及装置 |
CN112255688A (zh) * | 2020-10-27 | 2021-01-22 | 中国海洋石油集团有限公司 | 一种基于岩石物理理论的三维地震反演地层压力的方法 |
CN117491502A (zh) * | 2023-11-02 | 2024-02-02 | 陕西地矿创新研究院有限公司 | 一种不同co2饱和度岩石物理测试装置及实验方法 |
CN118011497A (zh) * | 2024-02-06 | 2024-05-10 | 中国海洋大学 | 基于时移地震属性差异识别海洋碳封存变化的方法及介质 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2002075363A1 (en) * | 2001-03-15 | 2002-09-26 | Shell Internationale Research Maatschappij B.V. | Time-lapse seismic surveying at sea |
US20100312480A1 (en) * | 2009-04-24 | 2010-12-09 | Hansteen Fredrik | Method for monitoring fluid flow in a multi-layered system |
CN102508294A (zh) * | 2011-10-20 | 2012-06-20 | 西北大学 | 一种利用时移地震勘探资料进行差异avo分析的方法 |
CN103576195A (zh) * | 2013-10-28 | 2014-02-12 | 西北大学 | 一种随压力变化的裂隙介质横波速度预测方法 |
CN103576196A (zh) * | 2013-10-28 | 2014-02-12 | 西北大学 | 一种随压力变化的孔隙介质横波速度预测方法 |
CN104678438A (zh) * | 2015-03-27 | 2015-06-03 | 西北大学 | 一种co2地质封存中四维地震资料co2分布预测的方法 |
CN104749617A (zh) * | 2013-12-26 | 2015-07-01 | 中国石油化工股份有限公司 | 一种多尺度裂缝储层正演模型建立方法 |
-
2015
- 2015-12-04 CN CN201510886011.9A patent/CN106842326B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2002075363A1 (en) * | 2001-03-15 | 2002-09-26 | Shell Internationale Research Maatschappij B.V. | Time-lapse seismic surveying at sea |
US20100312480A1 (en) * | 2009-04-24 | 2010-12-09 | Hansteen Fredrik | Method for monitoring fluid flow in a multi-layered system |
CN102508294A (zh) * | 2011-10-20 | 2012-06-20 | 西北大学 | 一种利用时移地震勘探资料进行差异avo分析的方法 |
CN103576195A (zh) * | 2013-10-28 | 2014-02-12 | 西北大学 | 一种随压力变化的裂隙介质横波速度预测方法 |
CN103576196A (zh) * | 2013-10-28 | 2014-02-12 | 西北大学 | 一种随压力变化的孔隙介质横波速度预测方法 |
CN104749617A (zh) * | 2013-12-26 | 2015-07-01 | 中国石油化工股份有限公司 | 一种多尺度裂缝储层正演模型建立方法 |
CN104678438A (zh) * | 2015-03-27 | 2015-06-03 | 西北大学 | 一种co2地质封存中四维地震资料co2分布预测的方法 |
Non-Patent Citations (5)
Title |
---|
JINFENG MA ET AL.: "AVO modeling of pressure-saturation effects in Weyburn CO2 sequestration", 《THE LEADING EDGE》 * |
李军等: "CO2驱正演模拟研究", 《中国地球科学联合学术年会2014》 * |
杨扬等: "CO2地质封存四维多分量地震监测技术进展", 《地球科学进展》 * |
由荣军等: "CO2地质储存的地震监测", 《物探与化探》 * |
郭建强等: "《二氧化碳地质储存技术方法概论》", 31 January 2015 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110441817A (zh) * | 2019-06-27 | 2019-11-12 | 西北大学 | 孔隙介质中基于4d pp-ps波反演ccus压力和饱和度变化的方法 |
CN111025396A (zh) * | 2020-01-06 | 2020-04-17 | 中国石油化工股份有限公司 | 基于人工智能算法的油藏物性参数地震预测方法 |
CN111025396B (zh) * | 2020-01-06 | 2021-11-05 | 中国石油化工股份有限公司 | 基于人工智能算法的油藏物性参数地震预测方法 |
CN111323817A (zh) * | 2020-04-15 | 2020-06-23 | 中国矿业大学(北京) | 基于深度学习的二氧化碳封存监测方法及装置 |
CN112255688A (zh) * | 2020-10-27 | 2021-01-22 | 中国海洋石油集团有限公司 | 一种基于岩石物理理论的三维地震反演地层压力的方法 |
CN112255688B (zh) * | 2020-10-27 | 2022-08-02 | 中国海洋石油集团有限公司 | 一种基于岩石物理理论的三维地震反演地层压力的方法 |
CN117491502A (zh) * | 2023-11-02 | 2024-02-02 | 陕西地矿创新研究院有限公司 | 一种不同co2饱和度岩石物理测试装置及实验方法 |
CN118011497A (zh) * | 2024-02-06 | 2024-05-10 | 中国海洋大学 | 基于时移地震属性差异识别海洋碳封存变化的方法及介质 |
Also Published As
Publication number | Publication date |
---|---|
CN106842326B (zh) | 2020-10-13 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Eastwood et al. | Seismic monitoring of steam-based recovery of bitumen | |
CN106842326A (zh) | 无横波速度测井时砂泥互层co2地质封存时移地震正演模拟方法 | |
White | Monitoring CO 2 storage during EOR at the Weyburn-Midale Field | |
CN106368691B (zh) | 基于岩石物理地震信息三维异常孔隙压力预测方法 | |
AU2011320352B2 (en) | Model based inversion of seismic response for determining formation properties | |
CN103544361A (zh) | 一种油气田开发中co2地质封存潜力的评价方法 | |
CN104678438B (zh) | 一种co2地质封存中四维地震资料co2分布预测的方法 | |
CN107807412B (zh) | 一种瓦斯地质溯源重构的方法 | |
Moscariello | Exploring for geo-energy resources in the Geneva Basin (Western Switzerland): opportunities and challenges | |
Barajas-Olalde et al. | Joint impedance and facies inversion of time-lapse seismic data for improving monitoring of CO2 incidentally stored from CO2 EOR | |
Kulikowski et al. | The cooper–eromanga petroleum province, australia | |
Utley | Time-lapse PP seismic to characterize stimulation and production behavior in the Niobrara and Codell Reservoirs Wattenberg Field, Colorado, US | |
Cao et al. | Pressure evolution and hydrocarbon migration-accumulation in the Moliqing fault depression, Yitong basin, Northeast China | |
Jin et al. | Application of CO2 injection monitoring techniques for CO2 EOR and associated geologic storage | |
Neil et al. | Monitoring Oil/Water Fronts by Direct Measurement (includes associated papers 23943 and 24103) | |
Ojo et al. | Quantitative modeling of the architecture and connectivity properties of reservoirs in ‘Royal’Field, Niger-Delta | |
Jansen | Seismic investigation of wrench faulting and fracturing at Rulison Field, Colorado | |
Marquez et al. | Improved reservoir characterization of a mature field through an integrated multidisciplinary approach. LL-04 reservoir, Tia Juana Field, Venezuela | |
Ehinola et al. | Seismic attributes mapping and 3D static modeling of reservoirs within “OYA” field, offshore Depobelt, Niger delta sedimentary basin, Nigeria | |
Airen et al. | Reservoir Hydrocarbon Volumetric Analysis of Sapele Deep Field, Niger Delta, Southern Nigeria | |
Gogri | Investigation and real-time monitoring for safe waste-water disposal with a focus on Arbuckle Group, Oklahoma | |
Morris et al. | Fault compartmentalization in a mature clastic reservoir: An example from elk hills field, California | |
Gargouri | Multicomponent 3D seismic interpretation of the Marcellus shale Bradford county, Pennsylvania | |
Riestenberg et al. | SOSRA Modeling-Based MVA Recommendations (Deliverable 9.1. a) | |
Soleimani et al. | Petroleum reservoir simulation, Ramin Oil Field, Zagros, Iran |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |