地 球 物 理 学 进 展
Vol.23 No.52008年10月(页码:1526～1531)
PRO GRESS IN GEOP H YSICS
Oct. 2008
基于流体饱和多孔隙介质波动方程
的时间推移地震反演
贺 英, 韩 波, 陈 勇
(哈尔滨工业大学数学系,哈尔滨150001)
摘 要 针对二维流体饱和多孔隙介质波动方程,设计了自适应同伦法和小波2自适应同伦法,在此基础上,将局域化思想引入时间推移地震反演,并使用小波2自适应同伦法和自适应同伦法反演两次监测数据,提出了一种适用于时间推移地震勘探的局域化反演算法.为了验证新方法是可行的,文中把它应用于同一地点不同时间的观测模型,通过数值实验表明了这种方法的有效性.
关键词 时间推移地震,局域化反演算法,小波2自适应同伦法,流体饱和多孔隙介质
中图分类号 P631 文献标识码 A 文章编号 100422903(2008)0521526206
Time 2lapse seismic data inversion based on the
w ave equ ation fluid 2saturated porous media
H E Y ing , HAN Bo , CH EN Y ong
(Depart ment of M at hematics ,Harbin I nstit ute of Technolog y ,Harbin 150001,China )
Abstract In this paper ,we propose a localized strategy using the wavelet 2homotopy method and the adaptive 2homoto 2py method as the interpretation tool for time 2lapse seismic inversion.In order to testify the new method ′s feasibility ,we apply the method to wave equations of fluid 2saturated porous media.The numerical results show that the new method is effective.
K eyw ords time 2lapse seismic ,localized inversion algorithm ,wavelet 2adaptive homotopy method ,fluid 2saturated por 2ous media
收稿日期 2007211210; 修回日期 2008202220.
基金项目 国家自然科学基金项目(40774056)和哈尔滨工业大学优秀青年教师培养计划资助项目(HITQNJ S.2008.032)联合资助.作者简介 贺英,女,1978年生,博士研究生,主要从事时间推移地震小波算法研究.(E 2mail :happybirdzhp @sohu.com )
0 引 言
近年来,随着油田勘探程度的加深,可供发现的
新油田减少,而对石油的需求却不断增加,这使得用于油藏管理的时间推移地震监测(Time 2lap se seis 2mic monitoring )技术越来越受到人们的关注.时间推移地震是指每间隔一定的时间对地质模型进行一次观测,对不同时间的观测数据进行处理,通过与初始观测数据相比较,来确定油藏随时间变化的情况,并综合利用岩石物理学,地质学和油藏工程资料,及时地对油藏进行动态监测,快速做出油藏评价,调整开发方案,对油田进行有效的开发,提高采收率.油藏动态监测的思想是由Nur.A [1]等人首次提出的,
他们用实验的方法观察到岩石中压力变化会导致地震速度变化,进而推测地震技术可以用来监测地层流体的动态变化.
在油田开采过程中,储层流体的变化会引起地震响应的变化.而时间推移地震反演的主要目的是利用地震响应接收数据的变化来识别储层参数的变化,并根据数据重复性高的特点,充分合理地利用前后勘探信息,提高解释流体流动情况的精度.就时间推移地震反演研究而言,现已取得了大量的成果:在国外,Abubakar [2],Gluck [3]分别探讨了时移地震波阻抗反演的可行性.Landro 等[4]给出了计算压力与饱和度的显式表达式,并通过对实际时间推移地震数据反演验证了压力与饱和度变化情况.Vaso 等[5]
在射线基方法的基础上,通过时间推移地震振幅变化反演了储层渗透率的变化.Veire等[6]针对时间推移多数据反演问题,提出Bayesian评估方法,估计了饱和度和压力的变化.在国内,李景叶等[7]依据油田实际情况,确定了适合计算疏松砂岩纵、横波速度和密度的岩石物理模型,并推导了时移地震AVO 反演计算方法.陈小宏等[8]在现有研究基础上提出基于岩石物理模型和混合优化算法的时移地震反演压力、饱和度变化方法并通过油藏参数的反演结果证明非线性反演方法不受油藏参数变化范围的影响,优于线性化反演方法.陈勇等[9]针对三维声波方程设计了时间推移地震的逐次递归算法.在稠油热采地震监测中,王新红[10]应用了弹性波阻抗反演方法.
尽管时间推移地震反演工作已取得了一定进展,但在算法方面的研究还不够完善,尤其是基于更符合流体在地层中的流动情况的Biot理论模型的时间推移地震反演.鉴于此,本文针对基于Biot理论的流体饱和多孔隙介质波动方程模型,构造了一种抗噪能力强,收敛范围广,程序易于实现的小波2自适应同伦法,并把它应用于时间推移地震前一次地震勘探模型反演.然后在此基础上,为了进一步减少反演计算量、提高计算速度,把局域化反演思想引入时间推移地震后一次勘探模型反演问题中,并使用自适应同伦法求解.最后针对理想的地质勘探模型给出数值模拟结果.
1 二维流体饱和多孔隙介质波动方程反演模型
在油藏开发过程中,储层流体变化所引起的地震响应变化是储层岩石弹性性质变化的反映.为此,对油藏监测而言,首先必须明确开采过程中油藏参数(孔隙度、渗透率、饱和度与压力、温度)的变化是如何影响储层岩石的弹性性质的.油藏变化利用孔隙度、压力、温度、饱和度等参数表示,利用地震表征油藏参数的变化需要建立油藏参数与地震参数的关系,即岩石物理模型,它是连接油藏参数变化与地震响应变化的桥梁.Biot[10,11]理论指出地下介质由流体和孔隙中充满流体的固体骨架组成,这一理论更加准确地描述了储层流体在地下介质中实际存在的情况,因此本文考虑使用流体饱和多孔隙介质波动方程作为时间推移地震反演模型.
忽略流体黏滞性的影响,描述地震波传播的二维流体饱和多孔隙介质模型以下控制方程为:
2
5
5xμ
5u x
5x+
5
5zμ
5u x
5z+μ
5u z
5x
+
5
5xλ
5u x
5x+λ
5u z
5z+
5
5xαM
5ωx
5x+αM
5ωz
5z
=ρ u x+ρ
f
ωx-f1,(1) 5
5xμ
5u x
5z+μ
5u z
5x+2
5
5zμ
5u z
5z
+
5
5zλ
5u x
x+λ
5u z
5z+
5
5zαM
5ωx
x+αM
5ωz
5z
=ρ u z+ρf ωz-f2,(2)
5
5xαM
5u x
5x+αM
5u z
5z+M
5ωx
5x+M
5ωz
5z
=ρfüx+m¨
ωx-f1,(3)
5
5zαM
5u x
5x+αM
5u z
5z+M
5ωx
5x+M
5ωz
5z
=ρ
f u z+m
ωz-f2,(4)其中u x(x,z,t),u z(x,z,t),ωx(x,z,t),ωz(x,z,t)为地震响应位移函数,f1(x,z,t),f2(x,z,t)是震源函数,x,z分别表示水平方向和垂直方向的坐标.参数μ,λ,α,M,m,ρ都与反演参数孔隙度β有关,关系如下:
k b=λb+
2
3
μ,α=1-k b
k s
,ρ=(1-β)ρs+βρf,
M=k s
α+βk s
k f
-1
,λ=λb+α2M,
其中k b,λb,k s,k f,ρs,ρf是流体饱和多孔隙介质中固体骨架及流体参数.
边界条件为:
5u x(x,z,t)
5x x=0=
5u x(x,z,t)
5x x=L=0, 5u x(x,z,t)
5z z=H=0,(5) 5u z(x,z,t)
5x x=0=
5u z(x,z,t)
5x x=L=0, 5u z(x,z,t)
5z z=H=0,(6) 5ωx(x,z,t)
x x=0=
5ωx(x,z,t)
x x=L=0, 5ωx(x,z,t)
5z z=H=0,(7) 5ωz(x,z,t)
5x x=0=
5ωz(x,z,t)
5x x=L=0, 5ωz(x,z,t)
5z z=H=0,(8)其中L,H为求解区域Ω的边界.
初始条件为:
7251
地 球 物 理 学 进 展23卷
u x t =0
=0,u z t =0
=0,ωx t =0
=0,ωz t =0
=0,
(9)
u x
t =1
=0,u z
t =1
=0,ωx
t =1
=0,ωz
t =1
=0.
(10)
如果以地震记录作为附加条件
u x (x ,0,t )=f (x ,t ),t ∈[0,T ],
(11)其中T 表示接收数据的时间.那么,(1～11)就构成了确定孔隙度β的二维流体饱和多孔隙介质波动方程反问题,即在油藏开采过程中,可以通过反演油藏储层的孔隙度来动态监测储层流体的流动情况.
在时间推移地震油藏监测中,用来反演地层参数的地震响应测量数据一般应由勘探实际地质模型的接收数据给出.为了简单起见,本文使用有限差分法把(1～10)进行差分离散计算得到反演地层参数的地震响应函数,并结合附加条件(11)将原反演问题
F (β
)=y ,(12)转化成求解极小化泛函
min β∈
Ω
‖F (β)-y ‖2
,(13)
其中y 是由已知的地震记录,即附加条件f k i (i =1…
L Δx ,k =1…T
Δt
)进行排列所组成的向量,
Δx ,Δt 表示进行差分离散时选取的空间和时间步长,F 是由模型方程(1)～(11)定义的一个非线性向量值函数F :
β→y .2 小波2自适应同伦法
引入局域化反演思想之前,对第一次地震勘探模型,即基础模型进行较为精确的反演计算是十分重要的,它能更好地识别油藏储层位置,从而准确地划定监测模型求解区域.针对求解基础模型反问题,本文设计了小波2自适应同伦法,这种方法不仅大范围收敛而且计算精度高.2.1 自适应同伦法
考虑求解反问题(12),定义非线性算子T 如下:
T (β):=F (β)-y ,β∈Ω,(13)则反问题就归结为求解非线性算子方程T (β
)=0.(14)这是一个非线性的、不适定的算子方程,考虑构造自适应同伦法进行求解.
构造如下同伦正则化方法
H (β,α)=(1-α)‖F (β)-y ‖2+α‖β-β0‖2
,
(15)
其中β0
是反演参数β的初值,参数α∈[0,1]既是同伦参数又是正则化参数.当α=1时,方程H (β,1)=0
的解是已知初值β0
.而当α=0时,方程H (β,0)=0的解就是所求反问题(12)的解.
划分[0,1]为1=α0>α1>…>αN =0.对于某个αk ,顺序地用某种迭代法求解H (β,αk )=0.由于α=1
时方程的解β0
已知,故可用前一个方程的近似解作为后一个方程的初值.假设第k 个方程的近似解βk
已求出,为求解第(k +1)个方程,可考虑使用阻尼高斯2牛顿法.在实际应用中,并不需要计算每一步的精确解,所以在使用阻尼高斯-牛顿法求解方程H (β,αk )=0时,只需迭代一步即可.迭代公式如下:
βk+1=βk -[(1-αk )F ′(βk )3F ′(βk )+αk I ]
-1 ・[(1-αk )F ′(βk )3(F (βk )-y )+αk
βk ], k =0,1,…,N -1
(16)把得到的第N 个方程近似解βk 作为阻尼高斯2牛顿法的初值,再利用阻尼高斯2牛顿迭代公式βk+1=βk -[F ′(βk )3F ′(βk )]-1・[F ′
(βk )3(F (βk )-y )],k =N ,N +1…(17)进行校正,就可以最终求出真解的近似值.这样公式
(16)和(17)一起称为同伦反演方法的迭代公式.
同伦反演方法克服了传统数值迭代法容易陷入局部收敛的不足,且对初始猜测的选取没有严格的,是求解非线性算子方程的一种快速稳定的大范围收敛方法.从上面的观察计算过程可以看出,同伦参数αk 的选择是至关重要的,它直接影响同伦路径的跟踪速度和精度.但是,在以往的文献中同伦参数大多等距的选取αk =
k
N
,k =1,2…,N.在反演过程中,这样选取αk 会导致某些极小值点被跳过,或
者在不存在极小值点区域仍然用等距步长进行路径跟踪.为此,在本文中,进一步修正了同伦参数的迭代步长,使它随求解过程自适应地选取,称修正了同伦参数的同伦反演方法为自适应同伦法,具体算法的步骤如下:
步骤1:选取反演参数的初值β0
,阈值ε和同伦参数初始步长h .
步骤2:令α1=α0-h ,α2=α1-h ,k =2.
步骤3:使用阻尼高斯2牛顿法求解方程H (β,
α1)=0的解β1
.
步骤4:当αk ≥
0时,执行下面循环使用阻尼高斯2牛顿法求解方程H (β,αk )=0的解βk ,令κ=‖βk -βk -1‖/‖βk -1-βk -2
‖,如果κ<ε,那么令h =2h ,否则令h =h 2
,令k =k +1,αk =
8
251
α
k-1-h.
2.2 小波2自适应同伦法
自适应同伦法虽具有大收敛范围和对初始值选取要求不高的优势,但当反问题存在多解或目标函数存在大量局部极小点,该方法的全局收敛性还需要进一步探讨,解决这类问题可采用小波多尺度分解的方法.
小波多尺度反演是利用小波分解算法把原始反问题分解在不同尺度上,根据不同尺度上目标函数的特征,逐步搜索全局极小点的一种反演方法,由于在大尺度上目标函数呈现较强的凸性和较少的局部极值点,可以先在大尺度上迭代反演,得到一个比较好的反演参数估计.将这个估计作为新的初值,在小一级的尺度上迭代反演,其反演相当于在全局极小点附近进行迭代修正,故能搜索到这一尺度上的“全局极小点”.依此类推,当尺度参数最终降至目标函数的原始尺度时,对应尺度上的“全局极小点”就可以视为真正意义上的全局极小点.从上面的分析可以看出,小波多尺度反演必须依赖某种反演方法才能实现反问题的数值求解.鉴于此,本文把小波多尺度反演和自适应同伦法结合,设计出一种新的反演方法,称为小波2自适应同伦法.小波2自适应同伦法的优点在于:在大尺度上,反演稳定且反演结果不受初值的影响,在一定程度上能避免其后的反演受局部极小困扰,加快收敛速度,进一步提高了计算效率.
不失一般性,将使用小波的3层分解来具体描述小波2自适应同伦法的求解过程.
步骤1:给定反演参数的初值β0,模型参数β3,震源函数f和监测数据y.
步骤2:把震源函数f和监测数据y利用紧支撑正交小波“DB4”作二维Mallat分解,即把反问题(12)分解为F1(β)=y1和F2(β)=y2.
步骤3:在最大尺度空间V2中,以β0为反演参数初值,使用自适应2同伦方法求解反问题F2(β)= y2,得到解的近似值β2.
步骤4:在小一级尺度空间V1中,以β2为反演参数初值,这时初值离真值比较接近,使用高斯2牛顿法求解反问题F1(β)=y1,得到解的近似值β1.
步骤5:在原问题所在空间V0中,以β1为反演参数初值,使用高斯2牛顿法求解反问题F(β)=y,得到解的近似值β,并视其为真正意义上的全局极小点.
从小波2自适应同伦法的具体实现过程可以看出:当把震源函数和监测数据利用紧支撑正交小波基分解后,随着尺度的增加,反问题(12)的复杂程度不断降低,表现出极值点少且距离较远的良好特征.在这种情况下,先在大尺度上使用自适应同伦法,得到孔隙度的估计,除此之外的其它尺度上都使用高斯2牛顿法求解,并把上一级尺度上的估计值作为下一级尺度的初值.不断细化尺度空间,直到找到满足原始反问题精度要求的全局极小点.
3 时间推移地震局域化反演算法
对于时间推移地震而言,在一定时间间隔内需要对油藏进行至少两次地震监测,为了区分两次勘探数据,把基础模型和监测模型的观测数据记为y1,y2,反演地层参数记为β1,β2.
此时,时间推移地震的基础模型反演可以描述为极小化问题(Ⅰ)
min
β
1
∈Ω
‖F(β1)-y1‖2,
其中Ω表示基础模型求解区域.针对极小化问题(Ⅰ),本文使用在第二部分设计的小波2自适应同伦法求解.因为小波2自适应同伦法是一种大范围收敛且计算精度较高的反演方法,所以通过对反演的结果成像,可以描述出油藏储层的近似位置.然后根据时间推移地震成像理论:非储层岩石的地震响应特征具有复测不变性,即在每次地震勘探中,非储层岩石的地震响应特征不随时间而变化,仅当油藏本身由于开发效应而导致地震响应特征的变化,针对监测模型反演问题,引入局域化反演思想,选定一个相对较小的区域Ω1(Ω1<Ω)作为新的求解区域,储层包含在Ω1内.此时,时间推移地震的监测模型反演可以描述为极小化问题(Ⅱ)
min
β
2
∈Ω
1
‖F(β2)-y2‖2,
区别于极小化问题(Ⅰ),监测模型求解区域较小,这时选用自适应同伦法求解极小化问题(Ⅱ)即可.引入局域化反演思想的主要优点是缩小求解计算区域,减少反演过程的计算量,提高计算精度,进而得到储层流体流动情况的精细刻画.
4 数值模拟
为了检验局域化反演算法的可靠性,本文设计了一个二维时间推移地震数据反演模型.图1表示参考模型的孔隙度β分布情况,油藏部分剖面图为一个椭圆形,孔隙度β初始值取为0.15,震源函数选取为正弦函数,其中求解区域Ω=[0m,1000m]
9251
地 球 物 理 学 进 展23卷
×[0m ,1000m ],震源位置在(500m ,500m )处.时间推移地震主要是比较前后两次地震勘探结果成像
的差异.在前面的讨论中提到,第一次地震勘探均在β∈Ω的求解区域内进行.而对于时间推移地震而言,除去第一次勘探以外,其余勘探的地质参数β仅仅在勘探区域Ω的一个局部区域Ω1内,也就是储层周围发生变化
.
图1 参考模型孔隙度分布
Fig.1 The reference model of porosity
在使用局域化反演算法具体计算时,首先使用
有限差分法计算反演地层参数的地震响应测量数据,空间步长取10m ,并在获得的数据上加3%的噪声.再以“DB4”作为紧支撑正交小波基,使用小波2
自适应同伦法求解min β∈
Ω
‖F (β)-y ‖2
来进行基础模型勘探反演,储层流体反演结果如图2.然后根据储
层流体的大概位置确定局域化反演算法的新求解区域Ω1=[0m ,550m ]×[300m ,550m ].假设经过一段时间的油藏开采,监测模型的椭圆型油藏有一部分孔隙度β从0.15变化到0.25,分布情况如图3.针对监测模型反演,选用自适应同伦法求解
min β∈Ω
1
‖F (β)-y ‖2
,反演结果如图4.比较两次地震勘探数据的反演结果的剖面图,我们可以看到,局域
化反演方法可以较好的反演出储层参数,它为动态监测储层流体的流动情况提供了解释依据.为了验证算法的效率,分别使用自适应同伦法和局域化反演算法进行监测模型的反演.取定反演参数的相对误差为1%,两种反演算法运算所需要的时间分别为675s 和430s.从运算时间可以看出,使用结合了局域化反演思想的反演方法,算法效率有了明显提高.数值模拟实验是在Inter Core Dou 1.6GHz/512m 下运行的
.
351
5 结 论
本文在时间推移地震成像理论假设基础上,设计了一种适合时间推移地震勘探的局域化反演算法,并以更符合实际地质情况的Biot流体饱和多孔隙介质波动方程模型为例,进行了储层参数孔隙度的反演研究.
数值结果表明,把局部化思想引入时间推移地震油藏监测数据反演中是可行的,通过使用小波2自适应同伦法对基础模型反演成像,可以初步地划定一个包含油藏储层的区域,进而在相对较小的求解区域内使用自适应同伦法进行监测模型反演,局域化反演算法不仅提高了计算的精度,还减少计算量,节省了计算时间.
本文所提出的方法仅仅在理想的地质模型上较好的实现了时间推移地震油藏监测储层参数反演,但对于实际应用还存在一些问题,如本文仅考虑了二维模型且地质模型求解区域相对实际问题的广度和深度还远远不够,这些将在以后的研究中加以探讨和改进.
参 考 文 献(References):
[1] Nur A,Tosaya C,V2Thanh D.Seismic monitoring of t hermal
enhanced oil recovery processes[J].SEG Technical Program Expanded Abstract s,1984,337～340.
[2] Abubakar A,Van den Berg P M,Fokkema J T.A feasibility
study on nonlinear inversion of time2lapse seismic data[J].
70t h Ann.Internet SEG.Expanded Abstract s,2001,16～1667.
[3] G luck S,Deschizeaux B,Mignot A,et al.Time2lapse imped2
ance inversion of post2stack seismic data[J].70t h Ann.Inter2 net SEG.Expanded Abstract s,2000,1509～1512.[4] Landro M.Discrimination between pressure and fluid satura2
tion changes from time2lapse seismic data.Geophysics[J].
2001,66:836～844.
[5] Vasco D W,Datta2Gupta A,Behrens R,et al.Seismic ima2
ging of reservoir flow properties:Time2lapse amplitude chan2 ges.Geophysics[J].2004,69(6):1425～1442.
[6] Veire H H,Borgos H G,Landro M.Stochastic of pressure
and sat uration changes from time2lapse multi component data [J].Geophysics,2007,55(6):805～818.
[7] 李景叶,陈小宏,金龙,等.时移地震AVO反演在油藏定量解
释中的应用[J].石油学报.2005,26(3):68～73.
Li Jingye,Chen Xiaohong,Jilong,et al.Application of time2 lapse seismic AVO inversion to quantitative interpretation of reservoir[J].Acta Pet rolei Sinica,2005,26(3):68～73. [8] 陈小宏,赵伟,马继涛,等.时移地震非线性反演压力与饱和度
变化研究[J].地球物理学进展,2006,21(3):830～836.
Chen X H,Zhao W,Ma J T,et al.Pressure and saturation change nonlinear inversion met hod using time2lapse seismic da2 ta[J].Progress in Geophysics(in Chinese),2006,21(3):830～836.
[9] 陈勇,韩波.时间推移地震反演的连续模型与算法[J].地球物
理学报,2006,49(4):11～1168.
Chen Y,Han B.Continous models and algorit hms for time2 lapse seismic inversion[J].Chinese J.Geophys.,2006,49
(4):11～1168.
[10] 王新红.弹性波阻抗反演在稠油热采地震监测中的应用[J].
地球物理学进展,2007,22(1):186～191.
Wang X H.Application of elastic impendence inversion in
seismic monitoring of t he heavy2oil t hermal recovery[J].Pro2
gress in Geophysics(in Chinese),2007,22(1):186～191..
[11] Biot M A.Theory of propagation of elastic waves in a fluid2
saturated porous solid:low2frequency range[J].J.Acoust.
Soc.Amer.,1956a,28:168～178.
[12] Biot M A.Theory of propagation of elastic waves in a fluid2
saturated porous solid:higher2frequency range[J].J.
Acoust.Soc.Amer.,1956b,28:179～191.
1351下载本文
