何兵红,方伍宝,胡光辉,刘定进,孙思宇
(1.中国石油化工股份有限公司石油物探技术研究院,江苏南京211103;2.中国石油大学(北京),北京102249)
随着计算机硬件的发展和高性能计算能力的提高,基于波动理论的地震全波形反演逐步成为研究热点[1]。全波形反演充分利用了地震数据的运动学特征和动力学特征,得到较高分辨率的速度模型,从而提高地震成像精度[2-5]。随着全波形反演研究的不断深入,多参数全波形反演辅助储层预测成为新的研究方向[6-9]。目前理论研究主要聚焦于弹性介质多参数全波形反演[10-13],但是由于当前采集得到的地震数据以纵波为主,且在应用推广中弹性波全波形反演受计算能力和存储能力的限制而难以有所突破[14-16],因而有必要开展声波方程的多参数全波形反演研究。
多参数间的耦合是声波方程多参数全波形反演面临的最主要问题之一。JANNANE等[17]讨论了L2范数下目标函数与长波长、中波长和短波长扰动之间的关系。但是该研究存在一定的局限,主要问题是最大偏移距只有1760m,限制了长波长的反演。ZHOU[18]等在速度-阻抗参数化模式下先固定阻抗值,利用潜波建立低波数速度场,再固定速度值,利用反射波进行高波数阻抗反演。董良国等[19]则通过变密度声波方程讨论了不同地震数据子集对应的目标函数随物性变化的非线性特征,提出了分步骤、分尺度的多参数全波形反演方案。张广智等[20-21]依靠初始模型对速度-密度进行顺序反演。YANG等[22-23]研究了声波方程中密度-速度参数化模式下的参数耦合现象,并给出了分偏移距的反演策略。石玉梅等[24]采用频率域单程声波延拓法进行波场模拟并反演得到地层密度和体积模量,提高了气藏成像的准确性。
因多参数全波形反演中目标函数的非线性源自波动方程的非线性,故需对波场产生机制进行理论分析,进而为多参数全波形反演提供理论指导。本文首先利用格林函数积分求解波动方程,在线性近似条件下推导了地震波振幅在6种参数化模式下(模量-密度、速度-密度、阻抗-密度、阻抗-速度、模量-速度、模量-阻抗)的敏感核函数;然后研究了不同参数化模式下参数扰动的辐射模式;再从理论上分析了波场随参数扰动的角度变化特征,对比了6种参数化下模式参数耦合特性,并提出了针对不同参数化模式的多参数全波形反演弱耦合化策略;最后进行了数值实验分析。
地震正演是地震反演的基础。声波方程正演通常不考虑介质密度的变化,只采用速度来描述地下介质的物性。但实际地层中,密度也是重要的物性参数,同时介质的变化通常伴随着密度的变化。考虑时间域变密度声波方程为:
(1)
式中:v(x)为速度;ρ(x)为密度;P(x,t)为压力场;f(x,t)为时间域震源;xs为震源位置坐标,δ(x-xs)表示将震源置于坐标点xs处。
根据介质的扰动理论,总介质(m)可以分为背景介质(m0)和扰动介质(δm),其中,m=m0+δm,对应的总波场(P)可以分解为背景波场(P0)和扰动波场(δP),其中,P=P0+δP。在背景介质中产生的背景波场满足波动方程:
(2)
式中:ρ0(x),v0(x)分别为背景密度和背景速度。其中,
(3)
扰动波场满足方程:
(4)
式中:ζ(x,t)为虚震源;P(x,t)为压力场,在Born近似条件下,背景波场远远大于扰动波场,虚震源项可以近似为:
(5)
式中:P0(x,ω)为背景波场。扰动波场对应的积分解可以用对应方程的格林函数G0(r,t;x)表示:
(6)
式中:P0(x,t)为背景波场。
在频率域,卷积项对应于乘积项,各个参数统一用模型的相对变化量δm/m表示,(6)式可以表示为:
(7)
整理(7)式可以得到速度相对扰动量和密度相对扰动量引起的波场扰动:
(8)
其中,Kvρ-v(x)和Kvρ-ρ(x)分别为速度项和密度项在速度-密度参数化模式下的敏感核函数:
速度敏感核函数Kvρ-v(x)表示地球介质用速度和密度表示时速度的相对变化量产生的扰动波场的强度。同理,密度敏感核函数Kvρ-ρ(x)表示在速度-密度参数化模式下密度的相对变化量产生的扰动波场的强度。
当用不同的参数化模式描述地球介质时,各个参数的敏感核函数具有不同的特征。声波方程常用的4种参数化模式包含速度-密度、模量-密度、阻抗-密度和阻抗-速度。本文又扩展了模量-速度、阻抗-模量两种参数化模式。我们从最常用的速度-密度声波方程出发推导了6种参数化模式的敏感核函数。
当采用体积模量和密度对地球介质进行模型参数化时,在Born近似条件下,总扰动波场是体积模量和密度相对变化量引起的扰动波场的线性叠加。
(11)
在弱扰动条件下,存在如下关系:
(12)
在速度-密度和模量-密度两种参数化模式下总波场扰动量相等,故有:
(13)
因此,可以得到模量-密度参数化模式下模量和密度的敏感核函数:
依据模量-密度参数化模式下敏感核函数推导流程,可以得到其它4种参数化模式的敏感核函数。
1) 阻抗-密度参数化模式下阻抗和密度的敏感核函数:
(16)
(17)
2) 阻抗-速度参数化模式下阻抗和速度的敏感核函数:
(18)
(19)
3) 模量-速度参数化模式下模量和速度的敏感核函数为:
(20)
(21)
4) 模量-阻抗参数化模式下模量和阻抗的敏感核函数:
(22)
(23)
本部分只给出这4种参数化模式下敏感核函数的解析表达式,敏感核函数的物理意义将在下一部分结合参数辐射模式进行阐述。
经过推导得到了不同参数化模式下扰动介质对扰动场的贡献量,同时因为介质的扰动存在一定的角度特征,所以根据不同参数化模式下的敏感核函数得到各种扰动介质对波场贡献量及角度特征,并制定了多参数全波形反演策略。
在三维介质中,存在关系式:
(24)
式中:θ为入射波和散射波之间的夹角。
因此速度-密度参数化模式下密度的敏感核函数可以表示为:
(25)
结合密度的敏感核函数表达式(10),得到速度和密度的辐射模式表达形式:
(26)
同理,我们可以分别得到其它参数化模式下的辐射模式表达式:
将垂直向下的方向视作入射方向,利用辐射模式的表达式得到空间各个方向的散射波能量。在各向同性介质中,任意两个垂直方向的平面内扰动参数的辐射模式具有完全等价性。
图1中箭头表示波的传播方向,入射波和反射波夹角为θ,反射波与x轴的夹角为φ。为了清楚地表示各参数的辐射模式,我们将依次展示单个垂直剖面内(φ=135°)的各个参数的辐射情况。
图1 地震波入射反射示意
图2a中速度辐射为标准的圆形,表明速度扰动引起的波场扰动在各个方向上振幅大小一样。图2b 中,密度扰动引起的波场扰动随着反射角的增大而减小。在小角度范围内相同的速度和密度相对扰动量引起的波场扰动量相当,因此仅依靠小角度地震数据难以区分速度和密度。因大角度地震数据主要由速度扰动产生,故可以通过大角度地震数据实现速度反演。全波形反演策略为先通过大角度地震数据反演速度,再利用小角度地震数据同时反演速度和密度。
图3为速度-密度模式下,波动方程数值模拟计算得到的扰动波场。背景速度和密度分别为3500m/s,2300kg/m3,扰动点位于模型的中心,扰动速度值为背景值的80%。炮点置于地表中心位置,位于扰动点正上方(本文中出现的6种参数化模式下的扰动波场的扰动点位置与图3完全相同,且采用了同样的观测方式)。由图3可知,在速度和密度相对扰动量相同的前提下,扰动引起的小角度波场扰动量基本一致,而密度扰动引起的大角度波场扰动量逐渐变弱,这与图2中的分析结果完全一致。
图2 φ=135°时速度-密度辐射情况a 速度;b 密度
图3 速度-密度模式下扰动波场a 速度扰动引起的扰动波场;b 密度扰动引起的扰动波场
图4a中模量辐射为圆形,表明模量扰动对波场扰动的贡献量在各个方向上相等。图4b中密度扰动对波场扰动的贡献量随着角度的增加先减小后增加,大角度的地震数据并不能解决模量和密度之间的耦合。当θ为90°时密度的改变量对扰动场的贡献为零(图5),可以先利用中角度地震数据先反演模量,然后利用小角度和大角度地震数据同时反演模量和密度。
图6为阻抗-密度辐射情况,阻抗扰动对扰动场的贡献量在各个方向相等,密度扰动在θ为小角度时对扰动场的贡献很小,相同的密度相对扰动量在θ为大角度时和阻抗相对扰动量对扰动波场的贡献量相当(图7)。对于阻抗-密度模式,可以先利用小角度地震数据实现阻抗反演,再利用大角度地震数据进行阻抗和密度的同时反演。
图8为阻抗-速度辐射情况,可以看出,阻抗扰动对小角度地震数据的扰动波场贡献较大,对大角度地震数据的扰动波场贡献较小。与此相反,速度扰动对大角度地震数据的扰动波场贡献较大,而对小角度地震数据的扰动波场贡献较小。图9中数值模拟计算结果验证了图8推导结果的正确性,可以首先利用大角度地震数据进行速度反演,再利用小角度地震数据进行阻抗反演。
图10为模量-速度辐射模式,模量扰动主要引起小角度的扰动波场,在该参数化模式下,速度扰动对大角度和小角度扰动波场贡献较大,图11验证了该结论。图12为模量-阻抗辐射模式,模量扰动主要引起大角度的扰动波场,在这种参数化模式下,阻抗扰动对大角度和小角度扰动波场贡献较大,对中角度的扰动波场贡献较小。图13进一步验证了该结论。
图4 φ=135°时模量-密度辐射情况a 模量;b 密度
图5 模量-密度模式下扰动波场a 模量扰动引起的扰动波场;b 密度扰动引起的扰动波场
图6 φ=135°时阻抗-密度辐射情况a 阻抗;b 密度
图8 φ=135°时阻抗-速度辐射情况a 阻抗;b 速度
图9 阻抗-速度模式下扰动波场a 阻抗扰动引起的扰动波场;b 速度扰动引起的扰动波场
图10 φ=135°时模量-速度辐射情况a 模量;b 速度
图11 模量-速度模式下扰动波场a 模量扰动引起的扰动波场;b 速度扰动引起的扰动波场
图12 φ=135°时模量-阻抗辐射情况a 模量;b 阻抗
图13 模量-阻抗模式下扰动波场a 模量扰动引起的扰动波场;b 阻抗扰动引起的扰动波场
综合以上6种参数化模式,在理想状态下具备所有角度的扰动波场时,阻抗-速度模式下能够利用大角度地震数据和小角度地震数据分别进行速度和阻抗反演,该种模式下两种参数耦合化程度较弱,是声波方程多参数反演中较为理想的参数化模式。阻抗-密度和速度-密度辐射图具有类似的特征,可先利用大角度或者小角度地震数据进行单参数反演,然后进行两种参数的同时反演。在模量-密度、模量-速度以及模量-阻抗3种参数化模式下进行多参数反演则难以区分不同参数的作用,多解性问题较严重。
理论上认为大角度地震数据是指120°≤θ<180°的地震数据,当θ=180°时,该数据已经不是反射波而是透射波。在地面地震观测方式中,大角度总是相对的,在数值模拟计算中只能尽可能提高偏移距与目标层深度的比值,来增加最大入射角度。
本文对速度-密度、模量-密度、阻抗-密度和速度-阻抗4种参数化模式进行了数值测试及分析。图14 为真实模型和对应的初始模型,初始模型由将平滑算法作用于真实模型得到。模型设计时,保证了两种参数化模式下波动方程正演所得到的模拟数据一致。图15为单炮观测数据,我们首先将全波场数据作为全波形反演的观测数据进行速度-密度同时反演。在图16速度-密度反演中,速度和密度真实异常得到了较好的反演,但是在速度中出现了密度异常引起的速度异常,并且密度项中速度异常引起的变化与密度本身的异常大小接近,即速度-密度同时反演时出现强烈的耦合现象。
本文还对模量-密度以及阻抗-密度参数化模式进行了数值测试。在模量-密度反演中模量的高低波数信息反演结果较为理想,但是密度的反演结果稳定性较差(图17)。
图14 真实模型及初始模型对比a 真实速度模型;b 初始速度模型;c 真实密度模型;d 初始密度模型
图15 单炮观测数据
在阻抗与密度的反演中,阻抗以小尺度信息为主,密度项中的部分低波数信息得到了恢复(图18)。但总体上这两种参数化反演中参数的耦合化现象比速度-密度参数化模式下的反演结果更为严重。
相比之下,速度-阻抗反演模式下,速度和阻抗表现出非常弱的耦合。同时,也不难发现两组反演中的速度表现出不同的特征。在速度-密度参数化模式下得到的速度异常具有清晰的边界,而速度-阻抗参数化模式下速度异常的中、低波数较为丰富,较高波数成分缺失,边界模糊(图19),该现象与前面的理论推导吻合。
图16 全波形反演速度(a)全波形反演密度(b)
图17 真实模量(a)、初始模量(b)与利用全波场观测数据同时反演模量(c)和密度(d)的结果
图18 真实阻抗(a)、初始阻抗(b)与利用全波场观测数据同时反演阻抗(c)和密度(d)的结果
速度-阻抗参数化模式下,速度异常主要与大角度地震数据有关,而阻抗主要与小角度地震数据有关。为了进一步降低速度和阻抗之间的耦合特征,本文采用中角度和大角度地震数据作为全波形反演的观测数据(图20)进行速度-阻抗的同时反演。反演结果如图21所示,可以看出速度中的高波数成分进一步减少,中低波数成分更加突出。同时阻抗对速度的干扰降低,即速度和阻抗之间的耦合现象在速度反演中得到了缓解。由于缺失了小角度地震数据,相对于图19d中的全波场观测数据的阻抗反演,仅采用中、大角度地震数据进行阻抗反演不但阻抗异常未能准确的反演,而且速度对阻抗的影响变大。因而在进行分角度反演时,需对反演结果进行合理的取舍。
图19 真实阻抗(a)、初始阻抗(b)与利用全波场观测数据同时反演速度(c)和阻抗(d)的结果
图20 单炮观测数据(中角度和大角度地震数据)
图21 利用中、大角度地震观测数据同时反演速度(a)和阻抗(b)
本文从理论上推导了声波方程不同参数化模式的敏感核函数,并分析了各种模式下参数的角度特征。速度-密度参数化模式下小角度地震数据的参数耦合严重不易区分,大角度地震数据则能够较好的区分;模量-密度参数化模式下小角度和大角度地震数据参数耦合都很严重,中角度地震数据耦合现象有所减弱;阻抗-密度参数化模式下大角度地震数据耦合严重,小角度地震数据各参数容易区分;速度-阻抗参数化模式下小角度和大角度地震数据参数耦合化现象都相对较弱。同时还分析了模量-速度和模量阻抗两种参数化模式,模量-速度参数化模式下小角度地震数据参数耦合严重,模量-阻抗参数化模式下大角度地震数据参数耦合严重。综上可知在理论上速度-阻抗参数化模式是最为理想的参数化模式。
通过理论模型的数值模拟计算分析了常用的速度-密度和速度-阻抗参数化模式下多参数反演。本文采用的同时反演策略,顺序反演会将所有的波场扰动归于其中一种参数扰动,容易导致不稳定性。数值模拟计算结果表明:速度-阻抗参数化模式下速度和阻抗的耦合化现象明显减弱。相对于速度-密度参数化模式反演得到的速度,速度-阻抗参数化模式下反演速度呈现出低波数特性,这是由于该模式下主要是大角度地震数据作用于速度反演。
理论推导在实际应用中往往受限。例如实际采集中最大角度有限,特别是陆上资料采集偏移距受限程度较大。此外,地震数据的大角度和小角度在理论上也没有定量的界限,因此需要根据观测系统制定合理的反演策略。