石油物探  2018, Vol. 57 Issue (5): 705-716, 725  DOI: 10.3969/j.issn.1000-1441.2018.05.010
0
文章快速检索     高级检索

引用本文 

何兵红, 方伍宝, 胡光辉, 等. 声波方程参数化模式及多参数全波形反演去耦合化策略[J]. 石油物探, 2018, 57(5): 705-716, 725. DOI: 10.3969/j.issn.1000-1441.2018.05.010.
HE Binghong, FANG Wubao, HU Guanghui, et al. Parameterization of acoustic wave equation and strategy for multi-parameter full waveform inversion[J]. Geophysical Prospecting for Petroleum, 2018, 57(5): 705-716, 725. DOI: 10.3969/j.issn.1000-1441.2018.05.010.

基金项目

国家科技重大专项(2017YFB0202904)资助

作者简介

何兵红(1985—), 女, 博士, 主要从事基于波动理论的地震波传播及全波形反演方面的研究。Email:hbh9517@163.com

文章历史

收稿日期:2017-11-20
改回日期:2018-01-22
声波方程参数化模式及多参数全波形反演去耦合化策略
何兵红1,2 , 方伍宝1 , 胡光辉1 , 刘定进1 , 孙思宇1     
1. 中国石油化工股份有限公司石油物探技术研究院, 江苏南京211103;
2. 中国石油大学(北京), 北京102249
摘要:参数耦合化是制约多参数全波形反演应用的关键因素之一。从速度-密度方程出发, 基于介质弱扰动线性假设, 利用波动方程的格林函数积分解推导了速度-密度、模量-密度、阻抗-密度、阻抗-速度、模量-速度及模量-阻抗6种参数化模式下的敏感核函数; 研究了每种参数化模式下各参数辐射模式, 分析总结了参数化模式下参数耦合性态; 提出在阻抗-速度参数化模式下先利用大角度地震数据进行速度反演, 再利用小角度地震数据进行阻抗反演的声波方程全波形反演优化策略。通过理论模型数值实验实现了速度-密度、模量-密度、阻抗-密度以及阻抗-速度4种参数化模式下的反演, 反演结果与理论推导一致。
关键词全波形反演    耦合化    参数化    敏感核函数    声波方程    
Parameterization of acoustic wave equation and strategy for multi-parameter full waveform inversion
HE Binghong1,2, FANG Wubao1, HU Guanghui1, LIU Dingjin1, SUN Siyu1     
1. Sinopec Geophysical Research Institute, Nanjing 211103, China;
2. China University of Petroleum, Beijing 102249, China
Foundation item: This research is financially supported by the National Science and Technology Major Project of China (Grant No.2017YFB0202904)
Abstract: Parameter coupling is one of the critical factors restricting multi-parameter full waveform inversion to practical application.In accordance with the linear hypothesis of media weak disturbance, the sensitive kernel function of six parameter patterns are derived, namely, velocity-density, modulus-density, impedance-density, impedance-velocity, modulus-velocity, and modulus-impedance.Green's function integral solution for the density-velocity wave equation is used for such derivation.The radiation pattern and parametric coupling behavior for each parameter pattern are analyzed.Subsequently, a two-step optimal full waveform inversion strategy is proposed.Firstly, velocity inversion is conducted using the large-angle data.Secondly, impedance inversion is launched using the small-angle data.The four types of parameter-pattern inversion of velocity-density, modulus-density, impedance-density, and impedance-velocity were tested by synthetic numerical experiments.The inversion results were in agreement with the theoretical deduction.
Key words: full waveform inversion    parameter coupling    parameterization    sensitive kernel function    acoustic wave equation    

随着计算机硬件的发展和高性能计算能力的提高, 基于波动理论的地震全波形反演逐步成为研究热点[1]。全波形反演充分利用了地震数据的运动学特征和动力学特征, 得到较高分辨率的速度模型, 从而提高地震成像精度[2-5]。随着全波形反演研究的不断深入, 多参数全波形反演辅助储层预测成为新的研究方向[6-9]。目前理论研究主要聚焦于弹性介质多参数全波形反演[10-13], 但是由于当前采集得到的地震数据以纵波为主, 且在应用推广中弹性波全波形反演受计算能力和存储能力的限制而难以有所突破[14-16], 因而有必要开展声波方程的多参数全波形反演研究。

多参数间的耦合是声波方程多参数全波形反演面临的最主要问题之一。JANNANE等[17]讨论了L2范数下目标函数与长波长、中波长和短波长扰动之间的关系。但是该研究存在一定的局限, 主要问题是最大偏移距只有1 760 m, 限制了长波长的反演。ZHOU[18]等在速度-阻抗参数化模式下先固定阻抗值, 利用潜波建立低波数速度场, 再固定速度值, 利用反射波进行高波数阻抗反演。董良国等[19]则通过变密度声波方程讨论了不同地震数据子集对应的目标函数随物性变化的非线性特征, 提出了分步骤、分尺度的多参数全波形反演方案。张广智等[20-21]依靠初始模型对速度-密度进行顺序反演。YANG等[22-23]研究了声波方程中密度-速度参数化模式下的参数耦合现象, 并给出了分偏移距的反演策略。石玉梅等[24]采用频率域单程声波延拓法进行波场模拟并反演得到地层密度和体积模量, 提高了气藏成像的准确性。

因多参数全波形反演中目标函数的非线性源自波动方程的非线性, 故需对波场产生机制进行理论分析, 进而为多参数全波形反演提供理论指导。本文首先利用格林函数积分求解波动方程, 在线性近似条件下推导了地震波振幅在6种参数化模式下(模量-密度、速度-密度、阻抗-密度、阻抗-速度、模量-速度、模量-阻抗)的敏感核函数; 然后研究了不同参数化模式下参数扰动的辐射模式; 再从理论上分析了波场随参数扰动的角度变化特征, 对比了6种参数化下模式参数耦合特性, 并提出了针对不同参数化模式的多参数全波形反演弱耦合化策略; 最后进行了数值实验分析。

1 声波方程多参数耦合模式 1.1 速度-密度参数化模式

地震正演是地震反演的基础。声波方程正演通常不考虑介质密度的变化, 只采用速度来描述地下介质的物性。但实际地层中, 密度也是重要的物性参数, 同时介质的变化通常伴随着密度的变化。考虑时间域变密度声波方程为:

$ \begin{array}{l} \frac{1}{{\rho \left( x \right){v^2}\left( x \right)}}\frac{\partial }{{\partial {t^2}}}P\left( {x, t} \right) - \nabla \cdot \left[ {\frac{1}{{\rho \left( x \right)}}\nabla P\left( {x, t} \right)} \right] = \\ \;\;\;\;\;\;f\left( {x, t} \right)\delta (x - {x_{\rm{s}}}) \end{array} $ (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。在背景介质中产生的背景波场满足波动方程:

$ \begin{array}{l} \frac{1}{{{\rho _0}\left( x \right)v_0^2\left( x \right)}}\frac{\partial }{{\partial {t^2}}}{P_0}\left( {x, t} \right) - \nabla \cdot \\ \;\;\;\;\;\;\;\;\left[ {\frac{1}{{{\rho _0}\left( x \right)}}\nabla {P_0}\left( {x, t} \right)} \right] = f\left( {x, t} \right)\delta (x - {x_{\rm{s}}}) \end{array} $ (2)

式中:ρ0(x), v0(x)分别为背景密度和背景速度。其中,

$ \begin{array}{l} \frac{1}{{{v^2}\left( x \right)}} = \frac{1}{{v_0^2\left( x \right)}} + \delta \left[ {\frac{1}{{{v^2}\left( x \right)}}} \right]\\ \frac{1}{{\rho \left( x \right)}} = \frac{1}{{{\rho _0}\left( x \right)}} + \delta \left[ {\frac{1}{{\rho \left( x \right)}}} \right] \end{array} $ (3)

扰动波场满足方程:

$ \begin{array}{l} \left\{ {\frac{1}{{{\rho _0}\left( x \right)}}\frac{1}{{v_0^2\left( x \right)}} + \frac{1}{{{\rho _0}\left( x \right)}}\delta \left[ {\frac{1}{{{v^2}\left( x \right)}}} \right] + \delta \left[ {\frac{1}{{\rho \left( x \right)}}} \right]\frac{1}{{v_0^2\left( x \right)}}} \right\}\\ \frac{\partial }{{\partial {t^2}}}\delta P\left( {x, t} \right) - \nabla \left[ { \cdot \frac{1}{{{\rho _0}\left( x \right)}}\nabla \delta P\left( {x, t} \right)} \right] = \zeta \left( {x, t} \right) \end{array} $ (4)

式中:ζ(x, t)为虚震源; P(x, t)为压力场, 在Born近似条件下, 背景波场远远大于扰动波场, 虚震源项可以近似为:

$ \begin{array}{l} \zeta \left( {x, t} \right) = \nabla \cdot \left\{ {\delta \left[ {\frac{1}{{\rho \left( x \right)}}} \right]\nabla {P_0}\left( {x, t} \right) - \frac{1}{{{\rho _0}\left( x \right)}} \cdot } \right.\\ \left. {\delta \left[ {\frac{1}{{{v^2}\left( x \right)}}} \right] + \delta \left[ {\frac{1}{{\rho \left( x \right)}}} \right]\frac{1}{{v_0^2\left( x \right)}}} \right\}\frac{\partial }{{\partial {t^2}}}{P_0}\left( {x, t} \right) \end{array} $ (5)

式中:P0(x, ω)为背景波场。扰动波场对应的积分解可以用对应方程的格林函数G0(r, t; x)表示:

$ \begin{array}{l} \delta P\left( {x, t} \right) = {\smallint _\mathit{\Omega }} - \left\{ {\frac{1}{{{\rho _0}\left( x \right)}}\delta \left[ {\frac{1}{{{v^2}\left( x \right)}}} \right] + \delta \left[ {\frac{1}{{\rho \left( x \right)}}} \right] \cdot } \right.\\ \left. {\frac{1}{{v_0^2\left( x \right)}}} \right\}\frac{\partial }{{\partial {t^2}}}{G_0}\left( {r, t;x} \right)*{P_0}\left( {x, t} \right){\rm{d}}{x^3} - {\smallint _\mathit{\Omega }}\delta \left[ {\frac{1}{{\rho \left( x \right)}}} \right] \cdot \\ \nabla {G_0}\left( {r, t;x} \right)*\nabla {P_0}\left( {x, t} \right){\rm{d}}{x^3} \end{array} $ (6)

式中:P0(x, t)为背景波场。

在频率域, 卷积项对应于乘积项, 各个参数统一用模型的相对变化量δm/m表示, (6)式可以表示为:

$ \begin{array}{l} \delta P\left( {x, \omega } \right) = {\smallint _\mathit{\Omega }} - \left\{ {\left[ {2\frac{{\delta v}}{{v\left( x \right)}} + \frac{{\delta \rho }}{{\rho \left( x \right)}}} \right]} \right\}\frac{1}{{\rho \left( x \right)}}\frac{1}{{{v^2}\left( x \right)}} \cdot \\ {\omega ^2}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right){\rm{d}}{x^3} + {\smallint _\mathit{\Omega }}\frac{{\delta \rho }}{{\rho \left( x \right)}}\frac{1}{{\rho \left( x \right)}} \cdot \\ \nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right){\rm{d}}{x^3} \end{array} $ (7)

整理(7)式可以得到速度相对扰动量和密度相对扰动量引起的波场扰动:

$ \begin{array}{l} \delta P\left( {x, \omega } \right) = {\smallint _\mathit{\Omega }}{K_{v\rho - v}}\left( x \right)\frac{{\delta v}}{{v\left( x \right)}}{\rm{d}}{x^3} + {\smallint _\mathit{\Omega }}{K_{v\rho - \rho }}\left( x \right) \cdot \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{\delta \rho }}{{\rho \left( x \right)}}{\rm{d}}{x^3} \end{array} $ (8)

其中, K-v(x)和K-ρ(x)分别为速度项和密度项在速度-密度参数化模式下的敏感核函数:

$ {K_{v\rho - v}}\left( x \right) = - \frac{{2{\omega ^2}}}{{\rho {v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) $ (9)
$ \begin{array}{l} {K_{v\rho - \rho }}\left( x \right) = - \frac{{{\omega ^2}}}{{\rho {v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) \end{array} $ (10)

速度敏感核函数K-v(x)表示地球介质用速度和密度表示时速度的相对变化量产生的扰动波场的强度。同理, 密度敏感核函数K-ρ(x)表示在速度-密度参数化模式下密度的相对变化量产生的扰动波场的强度。

当用不同的参数化模式描述地球介质时, 各个参数的敏感核函数具有不同的特征。声波方程常用的4种参数化模式包含速度-密度、模量-密度、阻抗-密度和阻抗-速度。本文又扩展了模量-速度、阻抗-模量两种参数化模式。我们从最常用的速度-密度声波方程出发推导了6种参数化模式的敏感核函数。

1.2 模量-密度参数化模式

当采用体积模量和密度对地球介质进行模型参数化时, 在Born近似条件下, 总扰动波场是体积模量和密度相对变化量引起的扰动波场的线性叠加。

$ \delta P\left( {x, \omega } \right) = {\smallint _\mathit{\Omega }}{K_{\kappa \rho - \kappa }}\left( x \right)\frac{{\delta \kappa }}{\kappa }{\rm{d}}{x^3} + {\smallint _\mathit{\Omega }}{K_{\kappa \rho - \rho }}\left( x \right)\frac{{\delta \rho }}{\rho }{\rm{d}}{x^3} $ (11)

在弱扰动条件下, 存在如下关系:

$ \frac{{\delta \kappa }}{\kappa } \approx 2\frac{{\delta v}}{v} + \frac{{\delta \rho }}{\rho } $ (12)

在速度-密度和模量-密度两种参数化模式下总波场扰动量相等, 故有:

$ \begin{array}{l} \delta P\left( {x, \omega } \right) = {\smallint _\mathit{\Omega }}{K_{\kappa \rho - \kappa }}\left( x \right)\frac{{\delta \kappa }}{\kappa }{\rm{d}}{x^3} + {\smallint _\mathit{\Omega }}{K_{\kappa \rho - \rho }}\left( x \right)\frac{{\delta \rho }}{\rho }{\rm{d}}{x^3}\\ \;\;\;\;\;\; = {\smallint _\mathit{\Omega }}2{K_{\kappa \rho - \kappa }}\left( x \right)\frac{{\delta v}}{{v\left( x \right)}}{\rm{d}}{x^3} + {\smallint _\mathit{\Omega }}[{K_{\kappa \rho - \rho }}\left( x \right) + \\ \;\;\;\;\;\;\;\;\;{K_{\kappa \rho - \kappa }}\left( x \right)]\frac{{\delta \rho }}{\rho }{\rm{d}}{x^3} \end{array} $ (13)

因此, 可以得到模量-密度参数化模式下模量和密度的敏感核函数:

$ {K_{\kappa \rho - \kappa }}\left( x \right) = - \frac{{{\omega ^2}}}{\kappa }{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) $ (14)
$ {K_{\kappa \rho - \rho }}\left( x \right) = \frac{1}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) $ (15)
1.3 其它参数化模式

依据模量-密度参数化模式下敏感核函数推导流程, 可以得到其它4种参数化模式的敏感核函数。

1) 阻抗-密度参数化模式下阻抗和密度的敏感核函数:

$ {K_{Ip\rho - Ip}}\left( x \right) = - \frac{{2{\omega ^2}}}{{\rho {v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) $ (16)
$ \begin{array}{l} {K_{Ip\rho - \rho }}\left( x \right) = - \frac{{{\omega ^2}}}{{\rho {v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) \end{array} $ (17)

2) 阻抗-速度参数化模式下阻抗和速度的敏感核函数:

$ \begin{array}{l} {K_{Ipv - Ip}}\left( x \right) = - \frac{{{\omega ^2}}}{{\rho {v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) \end{array} $ (18)
$ \begin{array}{l} {K_{Ipv - v}}\left( x \right) = - \frac{{{\omega ^2}}}{{\rho {v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) \end{array} $ (19)

3) 模量-速度参数化模式下模量和速度的敏感核函数为:

$ \begin{array}{l} {K_{\kappa v - \kappa }}\left( x \right) = - \frac{{{\omega ^2}}}{\kappa }{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) \end{array} $ (20)
$ {K_{\kappa v - v}}\left( x \right) = - \frac{2}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) $ (21)

4) 模量-阻抗参数化模式下模量和阻抗的敏感核函数:

$ \begin{array}{l} {K_{\kappa Ip - \kappa }}\left( x \right) = - \frac{{{\omega ^2}}}{\kappa }{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) \end{array} $ (22)
$ {K_{\kappa Ip - Ip}}\left( x \right) = \frac{2}{\rho }\nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {x, \omega } \right) $ (23)

本部分只给出这4种参数化模式下敏感核函数的解析表达式, 敏感核函数的物理意义将在下一部分结合参数辐射模式进行阐述。

2 参数辐射模式及全波形反演策略

经过推导得到了不同参数化模式下扰动介质对扰动场的贡献量, 同时因为介质的扰动存在一定的角度特征, 所以根据不同参数化模式下的敏感核函数得到各种扰动介质对波场贡献量及角度特征, 并制定了多参数全波形反演策略。

在三维介质中, 存在关系式:

$ \begin{array}{l} \nabla {G_0}\left( {r, \omega ;x} \right)\nabla {P_0}\left( {r, \omega ;x} \right) = \frac{{\partial {G_0}}}{{\partial x}}\frac{{\partial {P_0}}}{{\partial x}} + \frac{{\partial {G_0}}}{{\partial y}} \cdot \\ \frac{{\partial {P_0}}}{{\partial y}} + \frac{{\partial {G_0}}}{{\partial z}}\frac{{\partial {P_0}}}{{\partial z}} = - \cos \theta \frac{{{\omega ^2}}}{{{v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {r, \omega ;x} \right) \end{array} $ (24)

式中:θ为入射波和散射波之间的夹角。

因此速度-密度参数化模式下密度的敏感核函数可以表示为:

$ {K_{v\rho - \rho }}\left( x \right) = - 2{\cos ^2}\left( {\frac{\theta }{2}} \right)\frac{{{\omega ^2}}}{{\rho {v^2}}}{G_0}\left( {r, \omega ;x} \right){P_0}\left( {x, \omega } \right) $ (25)

结合密度的敏感核函数表达式(10), 得到速度和密度的辐射模式表达形式:

$ {Q_{v - \rho }} = \left[ { - 2, - 2{{\cos }^2}\left( {\frac{\theta }{2}} \right)} \right] $ (26)

同理, 我们可以分别得到其它参数化模式下的辐射模式表达式:

$ {Q_{\kappa \rho }} = \left( { - 1, - \cos \theta } \right) $ (27)
$ {Q_{Ip - \rho }} = \left[ { - 2, - 2{{\sin }^2}\left( {\frac{\theta }{2}} \right)} \right] $ (28)
$ {Q_{Ip - v}} = \left[ { - 2{{\cos }^2}\left( {\frac{\theta }{2}} \right), - 2{{\sin }^2}\left( {\frac{\theta }{2}} \right)} \right] $ (29)
$ {Q_{\kappa - v}} = \left[ { - 2{{\cos }^2}\left( {\frac{\theta }{2}} \right), - 2\cos \theta } \right] $ (30)
$ {Q_{\kappa - Ip}} = \left[ { - 2{{\sin }^2}\left( {\frac{\theta }{2}} \right), - 2\cos \theta } \right] $ (31)

将垂直向下的方向视作入射方向, 利用辐射模式的表达式得到空间各个方向的散射波能量。在各向同性介质中, 任意两个垂直方向的平面内扰动参数的辐射模式具有完全等价性。

图 1中箭头表示波的传播方向, 入射波和反射波夹角为θ, 反射波与x轴的夹角为φ。为了清楚地表示各参数的辐射模式, 我们将依次展示单个垂直剖面内(φ=135°)的各个参数的辐射情况。

图 1 地震波入射反射示意

图 2a中速度辐射为标准的圆形, 表明速度扰动引起的波场扰动在各个方向上振幅大小一样。图 2b中, 密度扰动引起的波场扰动随着反射角的增大而减小。在小角度范围内相同的速度和密度相对扰动量引起的波场扰动量相当, 因此仅依靠小角度地震数据难以区分速度和密度。因大角度地震数据主要由速度扰动产生, 故可以通过大角度地震数据实现速度反演。全波形反演策略为先通过大角度地震数据反演速度, 再利用小角度地震数据同时反演速度和密度。

图 2 φ=135°时速度-密度辐射情况 a速度; b密度

图 3为速度-密度模式下, 波动方程数值模拟计算得到的扰动波场。背景速度和密度分别为3 500 m/s, 2 300 kg/m3, 扰动点位于模型的中心, 扰动速度值为背景值的80%。炮点置于地表中心位置, 位于扰动点正上方(本文中出现的6种参数化模式下的扰动波场的扰动点位置与图 3完全相同, 且采用了同样的观测方式)。由图 3可知, 在速度和密度相对扰动量相同的前提下, 扰动引起的小角度波场扰动量基本一致, 而密度扰动引起的大角度波场扰动量逐渐变弱, 这与图 2中的分析结果完全一致。

图 3 速度-密度模式下扰动波场 a速度扰动引起的扰动波场; b密度扰动引起的扰动波场

图 4a中模量辐射为圆形, 表明模量扰动对波场扰动的贡献量在各个方向上相等。图 4b中密度扰动对波场扰动的贡献量随着角度的增加先减小后增加, 大角度的地震数据并不能解决模量和密度之间的耦合。当θ为90°时密度的改变量对扰动场的贡献为零(图 5), 可以先利用中角度地震数据先反演模量, 然后利用小角度和大角度地震数据同时反演模量和密度。

图 4 φ=135°时模量-密度辐射情况 a模量; b密度
图 5 模量-密度模式下扰动波场 a模量扰动引起的扰动波场; b密度扰动引起的扰动波场

图 6为阻抗-密度辐射情况, 阻抗扰动对扰动场的贡献量在各个方向相等, 密度扰动在θ为小角度时对扰动场的贡献很小, 相同的密度相对扰动量在θ为大角度时和阻抗相对扰动量对扰动波场的贡献量相当(图 7)。对于阻抗-密度模式, 可以先利用小角度地震数据实现阻抗反演, 再利用大角度地震数据进行阻抗和密度的同时反演。

图 6 φ=135°时阻抗-密度辐射情况 a阻抗; b密度
图 7 阻抗-密度模式下扰动波场 a阻抗扰动引起的扰动波场; b密度扰动引起的扰动波场

图 8为阻抗-速度辐射情况, 可以看出, 阻抗扰动对小角度地震数据的扰动波场贡献较大, 对大角度地震数据的扰动波场贡献较小。与此相反, 速度扰动对大角度地震数据的扰动波场贡献较大, 而对小角度地震数据的扰动波场贡献较小。图 9中数值模拟计算结果验证了图 8推导结果的正确性, 可以首先利用大角度地震数据进行速度反演, 再利用小角度地震数据进行阻抗反演。

图 8 φ=135°时阻抗-速度辐射情况 a阻抗; b速度
图 9 阻抗-速度模式下扰动波场 a阻抗扰动引起的扰动波场; b速度扰动引起的扰动波场

图 10为模量-速度辐射模式, 模量扰动主要引起小角度的扰动波场, 在该参数化模式下, 速度扰动对大角度和小角度扰动波场贡献较大, 图 11验证了该结论。图 12为模量-阻抗辐射模式, 模量扰动主要引起大角度的扰动波场, 在这种参数化模式下, 阻抗扰动对大角度和小角度扰动波场贡献较大, 对中角度的扰动波场贡献较小。图 13进一步验证了该结论。

图 10 φ=135°时模量-速度辐射情况 a模量; b速度
图 11 模量-速度模式下扰动波场 a模量扰动引起的扰动波场; b速度扰动引起的扰动波场
图 12 φ=135°时模量-阻抗辐射情况 a模量; b阻抗
图 13 模量-阻抗模式下扰动波场 a模量扰动引起的扰动波场; b阻抗扰动引起的扰动波场

综合以上6种参数化模式, 在理想状态下具备所有角度的扰动波场时, 阻抗-速度模式下能够利用大角度地震数据和小角度地震数据分别进行速度和阻抗反演, 该种模式下两种参数耦合化程度较弱, 是声波方程多参数反演中较为理想的参数化模式。阻抗-密度和速度-密度辐射图具有类似的特征, 可先利用大角度或者小角度地震数据进行单参数反演, 然后进行两种参数的同时反演。在模量-密度、模量-速度以及模量-阻抗3种参数化模式下进行多参数反演则难以区分不同参数的作用, 多解性问题较严重。

3 数值实验分析

理论上认为大角度地震数据是指120°≤θ<180°的地震数据, 当θ=180°时, 该数据已经不是反射波而是透射波。在地面地震观测方式中, 大角度总是相对的, 在数值模拟计算中只能尽可能提高偏移距与目标层深度的比值, 来增加最大入射角度。

本文对速度-密度、模量-密度、阻抗-密度和速度-阻抗4种参数化模式进行了数值测试及分析。图 14为真实模型和对应的初始模型, 初始模型由将平滑算法作用于真实模型得到。模型设计时, 保证了两种参数化模式下波动方程正演所得到的模拟数据一致。图 15为单炮观测数据, 我们首先将全波场数据作为全波形反演的观测数据进行速度-密度同时反演。在图 16速度-密度反演中, 速度和密度真实异常得到了较好的反演, 但是在速度中出现了密度异常引起的速度异常, 并且密度项中速度异常引起的变化与密度本身的异常大小接近, 即速度-密度同时反演时出现强烈的耦合现象。

图 14 真实模型及初始模型对比 a真实速度模型; b初始速度模型; c真实密度模型; d初始密度模型
图 15 单炮观测数据
图 16 全波形反演速度(a)全波形反演密度(b)

本文还对模量-密度以及阻抗-密度参数化模式进行了数值测试。在模量-密度反演中模量的高低波数信息反演结果较为理想, 但是密度的反演结果稳定性较差(图 17)。

图 17 真实模量(a)、初始模量(b)与利用全波场观测数据同时反演模量(c)和密度(d)的结果

在阻抗与密度的反演中, 阻抗以小尺度信息为主, 密度项中的部分低波数信息得到了恢复(图 18)。但总体上这两种参数化反演中参数的耦合化现象比速度-密度参数化模式下的反演结果更为严重。

图 18 真实阻抗(a)、初始阻抗(b)与利用全波场观测数据同时反演阻抗(c)和密度(d)的结果

相比之下, 速度-阻抗反演模式下, 速度和阻抗表现出非常弱的耦合。同时, 也不难发现两组反演中的速度表现出不同的特征。在速度-密度参数化模式下得到的速度异常具有清晰的边界, 而速度-阻抗参数化模式下速度异常的中、低波数较为丰富, 较高波数成分缺失, 边界模糊(图 19), 该现象与前面的理论推导吻合。

图 19 真实阻抗(a)、初始阻抗(b)与利用全波场观测数据同时反演速度(c)和阻抗(d)的结果

速度-阻抗参数化模式下, 速度异常主要与大角度地震数据有关, 而阻抗主要与小角度地震数据有关。为了进一步降低速度和阻抗之间的耦合特征, 本文采用中角度和大角度地震数据作为全波形反演的观测数据(图 20)进行速度-阻抗的同时反演。反演结果如图 21所示, 可以看出速度中的高波数成分进一步减少, 中低波数成分更加突出。同时阻抗对速度的干扰降低, 即速度和阻抗之间的耦合现象在速度反演中得到了缓解。由于缺失了小角度地震数据, 相对于图 19d中的全波场观测数据的阻抗反演, 仅采用中、大角度地震数据进行阻抗反演不但阻抗异常未能准确的反演, 而且速度对阻抗的影响变大。因而在进行分角度反演时, 需对反演结果进行合理的取舍。

图 20 单炮观测数据(中角度和大角度地震数据)
图 21 利用中、大角度地震观测数据同时反演速度(a)和阻抗(b)
4 讨论与认识

本文从理论上推导了声波方程不同参数化模式的敏感核函数, 并分析了各种模式下参数的角度特征。速度-密度参数化模式下小角度地震数据的参数耦合严重不易区分, 大角度地震数据则能够较好的区分; 模量-密度参数化模式下小角度和大角度地震数据参数耦合都很严重, 中角度地震数据耦合现象有所减弱; 阻抗-密度参数化模式下大角度地震数据耦合严重, 小角度地震数据各参数容易区分; 速度-阻抗参数化模式下小角度和大角度地震数据参数耦合化现象都相对较弱。同时还分析了模量-速度和模量阻抗两种参数化模式, 模量-速度参数化模式下小角度地震数据参数耦合严重, 模量-阻抗参数化模式下大角度地震数据参数耦合严重。综上可知在理论上速度-阻抗参数化模式是最为理想的参数化模式。

通过理论模型的数值模拟计算分析了常用的速度-密度和速度-阻抗参数化模式下多参数反演。本文采用的同时反演策略, 顺序反演会将所有的波场扰动归于其中一种参数扰动, 容易导致不稳定性。数值模拟计算结果表明:速度-阻抗参数化模式下速度和阻抗的耦合化现象明显减弱。相对于速度-密度参数化模式反演得到的速度, 速度-阻抗参数化模式下反演速度呈现出低波数特性, 这是由于该模式下主要是大角度地震数据作用于速度反演。

理论推导在实际应用中往往受限。例如实际采集中最大角度有限, 特别是陆上资料采集偏移距受限程度较大。此外, 地震数据的大角度和小角度在理论上也没有定量的界限, 因此需要根据观测系统制定合理的反演策略。

参考文献
[1]
VIRIEUX J, OPERTO S. An overview of full-wave form inversion in exploration geophysics[J]. Geophysics, 2009, 74(6): WCC1-WCC26. DOI:10.1190/1.3238367
[2]
TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 1984, 49(8): 1259-1266. DOI:10.1190/1.1441754
[3]
PRATT R G, SHIN C, HICK G J. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion[J]. Geophysical Journal International, 1998, 133(2): 341-362. DOI:10.1046/j.1365-246X.1998.00498.x
[4]
PRATT R G, WORTHINGTON M H. Inversion theory applied to multi-source cross-hole tomography.Part 1:Acoustic wave-equation method[J]. Geophysical Prospecting, 1990, 38(3): 287-310. DOI:10.1111/gpr.1990.38.issue-3
[5]
BROSSIER R, TIVIER L M, VIRIEUX J, et al. A review of some methodological developments on fullwaveform inversion tackled in the SEISCOPE group[J]. Geophysical Prospecting for Petroleum, 2017, 56(1): 3-19.
[6]
张广智, 姜岚杰, 孙昌路, 等. 基于照明预处理的分步多参数时间域声波全波形反演方法研究[J]. 石油物探, 2017, 56(1): 31-37.
ZHANG G Z, JIANG L J, SUN C L, et al. The stepped multi-parameter FWI of acoustic media in time-domain by L-BFGS method with illumination analysis[J]. Geophysical Prospecting for Petroleum, 2017, 56(1): 31-37. DOI:10.3969/j.issn.1000-1441.2017.01.004
[7]
JEONG W, LEE H Y, MIN D J. Full waveform inversion strategy for density in the frequency domain[J]. Geophysical Journal International, 2012, 188(3): 221-242.
[8]
KÖHN D, NIL D D, KURZMANN A, et al. On the influence of model parametrization in elastic full waveform tomography[J]. Geophysical Journal International, 2012, 191(1): 325-345. DOI:10.1111/gji.2012.191.issue-1
[9]
OH J W, MIN D J. Multi-parameter full waveform inversion using Poisson's ratio for elastic media[J]. Exploration Geophysics, 2017, 48(4): 456-475.
[10]
SINGH S, SEARS T, ROBERTS M, et al. Full elastic waveform inversion:future of quantitative seismic imaging[J]. Expanded Abstracts of 78th Annual Internat SEG Mtg, 2008, 1905-1909.
[11]
BROSSIER R, OPERTO S, VIRIEUX J. Which data residual norm for robust elastic frequency-domain full waveform inversion?[J]. Geophysics, 2010, 75(3): R37-R46. DOI:10.1190/1.3379323
[12]
SEARS T J, BARTON P J, SINGH S C. Elastic full waveform inversion of multicomponent ocean-bottom cable seismic data:application to Alba Field, U.K.North Sea[J]. Geophysics, 2010, 75(6): R109-R119. DOI:10.1190/1.3484097
[13]
OPERTO S, GHOLAMI Y, PRIEUX V, et al. A guided tour of multiparameter full waveform inversion for multicomponent data:from theory to practice[J]. The Leading Edge, 2013, 32(9): 1040-1054. DOI:10.1190/tle32091040.1
[14]
RAKNES E B, ARNTSEN B. Challenges and solutions for performing 3D time-domain elastic full-waveform inversion[J]. The Leading Edge, 2017, 36(1): 88-93. DOI:10.1190/tle36010088.1
[15]
RAKNES E B, WEIBULL W. Efficient 3D elastic full-waveform inversion using wavefield reconstruction methods[J]. Geophysics, 2016, 81(2): R45-R55. DOI:10.1190/geo2015-0185.1
[16]
GROOS L, SCHÄFER M, FORBRIGER T, et al. Application of a complete workflow for 2D elastic full-waveform inversion to recorded shallow-seismic Rayleigh waves[J]. Geophysics, 2017, 82(2): R109-R1117. DOI:10.1190/geo2016-0284.1
[17]
JANNANE M, BEYDOUN W, CRASE E, et al. Wavelengths of earth structures that can be resolved from seismic reflection data[J]. Geophysics, 1989, 54(7): 906-910. DOI:10.1190/1.1442719
[18]
ZHOU W, BROSSIER R, OPERTO S, et al. Acoustic multiparameter full-waveform inversion of diving and reflected waves through a hierachical scheme[J]. Expanded Abstracts of 84th Annual Internat SEG Mtg, 2014, 1249-1253.
[19]
董良国, 迟本鑫, 陶纪霞, 等. 声波全波形反演目标函数性态[J]. 地球物理学报, 2013, 56(10): 3445-3460.
DONG L G, CHI B X, TAO J X, et al. Objective function behavior in acoustic full waveform inversion[J]. Chinese Journal of Geophysics, 2013, 56(10): 3445-3460. DOI:10.6038/cjg20131020
[20]
张广智, 孙昌路, 潘新朋, 等. 变密度声波全波形反演中密度影响因素及反演策略[J]. 吉林大学学报(地球科学版), 2016, 46(5): 1550-1560.
ZHANG G Z, SUN C L, PAN X P, et al. Influence factors and strategy of inversion for density of acoustic full waveform inversion with variable density[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(5): 1550-1560.
[21]
张广智, 姜岚杰, 孙昌路, 等. 基于照明预处理的分步多参数时间域声波全波形反演方法研究[J]. 石油物探, 2017, 56(1): 31-37.
ZHANG G Z, JIANG L J, SUN C L, et al. The stepped multi-parameter FWI of acoustic media in time-domain by LBFGS method with illumination analysis[J]. Geophysical Prospecting for Petroleum, 2017, 56(1): 31-37. DOI:10.3969/j.issn.1000-1441.2017.01.004
[22]
杨积忠, 刘玉柱, 董良国. 变密度声波方程多参数全波形反演策略[J]. 地球物理学报, 2014, 57(2): 628-643.
YANG J Z, LIU Y Z, DONG L G. A multi-parameter full waveform inversion strategy for acoustic media with variable density[J]. Chinese Journal of Geophysics, 2014, 57(2): 628-643.
[23]
YANG J Z, LIU Y Z, DONG L G. A multi-parameter full waveform inversion strategy in acoustic media[J]. Expanded Abstracts of 76thAnnual Internat EAGE Mtg, 2014, 316-323.
[24]
石玉梅, 张研, 姚逢昌, 等. 基于声学全波形反演的油气藏地震成像方法[J]. 地球物理学报, 2014, 57(2): 607-617.
SHI Y M, ZHANG Y, YAO F C, et al. Methodology of seismic imaging for hydrocarbon reservoir based on acoustic full waveform inversion[J]. Chinese Journal of Geophysics, 2014, 57(2): 607-617.