速度是地震成像的关键, 在当前单参数全波形反演中占主导地位。为了适应复杂介质, 提高地震成像和储层预测精度, 降低反演结果解释的多解性, 需要利用多种参数提供更加充分信息用以确认地下地层结构和内部充填物质特征^[1-2]。对于复杂油气藏, 仅利用基于走时信息的射线类反演方法已不能满足当前对精细地震勘探的要求。近十年全波形反演技术发展迅速, 同时我国高性能计算集群计算能力不断提高, 基于波动方程的多参数全波形反演可以充分利用地震波的运动学和动力学特征对复杂油气藏进行精细刻画, 进而为更加准确的储层预测奠定了基础^[3-4]。

实际地球介质是复杂的, 利用各种复杂介质地球物理模型能够更加逼近真实的地下介质, 但基于复杂介质的波动方程地震波数值模拟导致全波形反演计算效率大幅度降低。因此, 当前仍然以基于声介质的纵波全波形反演为主^[5-6]。

现有的波动方程通常采用纵波速度和密度两个参数对声介质进行描述, 因而当提及声波方程多参数全波形反演时, 通常是纵波速度、密度双参数反演^[7-9]。由于阻抗在储层预测中占据更加重要的地位, 综合利用速度、密度、阻抗三参数信息可进一步提高储层预测精度。针对声波方程无法直接反演三参数问题, 本文提出了一种基于阻抗-速度方程和速度-密度方程的多参数全波形反演方法:首先利用阻抗-速度方程实现低波数速度模型和阻抗模型的全波形反演, 通过方程转换再利用速度-密度方程实现高波数速度反演和密度反演。最后利用SEG的Overthrust推覆体模型数据进行测试, 验证了方法的有效性。

1 方法原理 1.1 三参数反演原理

全波形反演目标函数的非线性源自波动方程的非线性, 需要从波场产生机制进行深入理论剖析, 进而为多参数全波形反演提供理论指导。对于声介质可以利用纵波速度、密度、阻抗、模量4个参数进行描述。何兵红等^[10]研究了6种参数化模式下不同参数的辐射模式, 以及每种模式的反演策略。参数之间的耦合性一直是多参数反演的难点问题, 反演的参数个数越多, 参数之间的耦合性越强, 一定程度上降低了反演的精度和可靠性。声波介质中描述储层特征最基本的弹性参数是速度、密度和阻抗, 其它参数如拉梅参数、体积模量等都可以通过这3个参数转换得到。现有的声波方程最多只能用2个参数来表示, 三参数同时反演不可行。因此本文需要从双参数全波形反演的特征出发, 研究一种可行的基于波动方程转换的反演方法。

基于速度、密度和阻抗可形成3种参数化模式, 如图 1, 利用3种参数化模式通过方程转换可形成三参数全波形反演。方程转换有多种选择方案:可以将密度作为中间过渡参数形成基于速度-密度和密度-阻抗方程的三参数全波形反演(称为方案一); 可以将阻抗作为中间过渡参数形成基于速度-阻抗和阻抗-密度方程的三参数全波形反演(称为方案二); 亦可形成将速度作为中间过渡参数的基于密度-速度和速度-阻抗方程的三参数全波形反演(称为方案三)。同时考虑2种方程先后反演顺序的变化, 总共有6种方案可选。下面从速度-密度、密度-阻抗和阻抗-速度的辐射模式出发(图 1), 研究几种方案的可行性。这3种模式下各参数的敏感核函数和辐射模式可以参考文献[10]和文献[11], 具体推导过程本文不再陈述。

对于方案一, 密度在速度-密度模式和密度-阻抗中分别影响小角度和大角度地震数据, 可以先开展密度-阻抗反演, 得到阻抗和低波数的密度, 然后进行速度-密度反演得到速度和高波数密度。但密度反演的稳定性一直是难以解决的问题, 以密度作为2种模式的衔接存在很大的隐患, 且该模式下无法得到可靠的速度和阻抗。

对于方案二, 阻抗应该按照从低波数到高波数的多尺度策略进行反演。在密度-阻抗模式中, 阻抗和密度在大角度耦合严重, 难以得到较为准确的低波数阻抗。在速度-阻抗模式中, 阻抗主要影响小角度地震数据, 对应的反演中得到的是高波数的阻抗。所以, 该种方案需要首先解决低波数阻抗反演难题, 并不是理想的实施方案。

对于方案三, 速度是地震勘探中最关键的因素, 以速度作为多种双参数全波形反演之间的过渡和衔接是保证基于波动方程转换反演方法的基本思想。对于阻抗-速度参数化模式, 在小角度主要是阻抗扰动对扰动波场的贡献, 在大角度主要是速度扰动对扰动波场的贡献, 是一种较为理想的多参数全波形反演模式。所以, 本文采用这种方案进行多参数全波形反演。

下面分别对阻抗-速度和速度-密度2种模式下的双参数反演进行简单介绍。

1.2 阻抗-速度模式下双参数反演

目前, 常用的方程是时间域变密度声波方程, 即速度-密度方程:

$ \begin{array}{*{20}{l}} {\frac{1}{{\rho (\mathit{\boldsymbol{x}}){v^2}(\mathit{\boldsymbol{x}})}}\frac{\partial }{{\partial {t^2}}}\mathit{\boldsymbol{P}}(\mathit{\boldsymbol{x}},t) - \boldsymbol{\nabla} \cdot \left[ {\frac{1}{{\rho (\mathit{\boldsymbol{x}})}}\boldsymbol{\nabla} \mathit{\boldsymbol{P}}(\mathit{\boldsymbol{x}},t)} \right]}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = f(\mathit{\boldsymbol{x}},t)\delta (\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_s})} \end{array} $

(1)

式中:v(x)为速度; ρ(x)为密度; P(x, t)为压力场; f(x, t)为时间域震源; x_s为震源位置坐标; δ(x－x_s)表示将震源放置于坐标点x_s处。一阶速度-应力方程实用性更强, 为了反演得到阻抗值, 我们建立了与阻抗相关的自伴随时间域一阶声波波动方程, 并引入卷积完全匹配层(CPML)边界条件:

$ \left\{ \begin{array}{l} \frac{1}{{I_{\rm{P}}^2}}\frac{{\partial p}}{{\partial t}} = \frac{{{v_{\rm{P}}}}}{{{I_{\rm{P}}}}}\left[ {\left( {\frac{1}{{{\chi _x}}}\frac{\partial }{{\partial x}} + {\varphi _x}} \right){v_x} + \left( {\frac{1}{{{\chi _y}}}\frac{\partial }{{\partial y}} + {\varphi _y}} \right) \cdot } \right.\\ \left. {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {v_y} + \left( {\frac{1}{{{\chi _z}}}\frac{\partial }{{\partial z}} + {\varphi _z}} \right){v_z}} \right]\\ \frac{{\partial {v_x}}}{{\partial t}} = \frac{{{v_{\rm{P}}}}}{{{I_{\rm{P}}}}}\left( {\frac{1}{{{\chi _x}}}\frac{\partial }{{\partial x}} + {\varphi _x}} \right)p\\ \frac{{\partial {v_y}}}{{\partial t}} = \frac{{{v_{\rm{P}}}}}{{{I_{\rm{P}}}}}\left( {\frac{1}{{{\chi _y}}}\frac{\partial }{{\partial y}} + {\varphi _y}} \right)p\\ \frac{{\partial {v_z}}}{{\partial t}} = \frac{{{v_{\rm{P}}}}}{{{I_{\rm{P}}}}}\left( {\frac{1}{{{\chi _z}}}\frac{\partial }{{\partial z}} + {\varphi _z}} \right)p \end{array} \right. $

(2)

式中:p为压力场; I_P为纵波阻抗; v_P为地震波传播的纵波速度; v_x, v_y, v_z为x, y, z三个方向的质点震动速度; χ_i, d_i和α_i表示在i(i=x, y, z)方向阻尼项; φ_i表示在i方向的递归函数。

$ \begin{array}{*{20}{l}} {{\varphi _i}(n) = {\varphi _i}(n - 1){\rm{exp}}\left[ { - \left( {{\alpha _i} + \frac{{{d_i}}}{{{\chi _i}}}} \right)\delta t} \right] + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left\{ {{\rm{exp}}\left[ { - \left( {{\alpha _i} + \frac{{{d_i}}}{{{\chi _i}}}} \right)\delta t} \right] - 1} \right\}\frac{{{d_i}}}{{{\chi _i}({\alpha _i}{\chi _i} + {d_i})}}{\partial _i}} \end{array} $

(3)

式中:n代表迭代次数。

根据模拟数据和观测数据波场误差的二范数形式建立对应的目标泛函:

$ {E_1} = \frac{1}{2}|{F_{\rm{R}}}[{I_{\rm{P}}},{v_{\rm{P}}}] - {d_{{\rm{obs}}}}{|^2} $

(4)

式中:F[I_P, v_P]表示采用阻抗-速度方程得到模拟波场; F_R[I_P, v_P]表示在检波点处接收到的地震记录; d_obs为观测数据。采用伴随状态法可以得到阻抗-速度模式下关于速度和阻抗的梯度:

$ {\rm{gra}}{{\rm{d}}_v}{E_1}[{I_{\rm{P}}},{v_{\rm{P}}}] = - \langle q,{\rm{D}}{F_v}[{I_{\rm{P}}},{v_{\rm{P}}}]\rangle $

(5a)

$ {\rm{gra}}{{\rm{d}}_{{I_{\rm{P}}}}}{E_1}[{I_{\rm{P}}},{v_{\rm{P}}}] = - \langle q,{\rm{D}}{F_{{I_{\rm{P}}}}}[{I_{\rm{P}}},{v_{\rm{P}}}]\rangle $

(5b)

式中:q为伴随波场, 是波场残差的反传波场; 符号〈·〉表示内积; DF_v[I_P, v_P]、DF_IP[I_P, v_P]分别表示速度-阻抗方程对速度和阻抗的一阶偏导数。

$ {\rm{D}}{F_v}[{I_{\rm{P}}},{v_{\rm{P}}}] = - \frac{1}{{{I_{\rm{P}}}v_{\rm{P}}^2}}\frac{{\partial p}}{{\partial t}} - \frac{{{I_P}}}{{v_P^2}}\left( {\frac{{\partial {v_x}}}{{\partial t}} + \frac{{\partial {v_y}}}{{\partial t}} + \frac{{\partial {v_z}}}{{\partial t}}} \right) $

(6a)

$ {\rm{D}}{F_{{I_{\rm{p}}}}}[{I_{\rm{P}}},{v_{\rm{P}}}] = - \frac{1}{{I_{\rm{P}}^2{v_{\rm{P}}}}}\frac{{\partial p}}{{\partial t}} + \frac{1}{{{v_{\rm{P}}}}}\left( {\frac{{\partial {v_x}}}{{\partial t}} + \frac{{\partial {v_y}}}{{\partial t}} + \frac{{\partial {v_z}}}{{\partial t}}} \right) $

(6b)

则可进一步得到更新后的速度和阻抗:

$ {v_{1,k + 1}} = {v_{1,k}} + \alpha {\rm{gra}}{{\rm{d}}_v}{E_1}[{I_{\rm{P}}},{v_{\rm{P}}}] $

(7a)

$ {I_{{\rm{P}},k + 1}} = {I_{{\rm{P}},k}} + \alpha {\rm{gra}}{{\rm{d}}_{{I_{\rm{P}}}}}{E_1}[{I_{\rm{P}}},{v_{\rm{P}}}] $

(7b)

式中:v_{1, k}中下标1表示在阻抗-速度模式中的速度; k代表迭代次数; α为迭代步长。

在该模式下我们可以通过中大角度的地震数据得到大尺度的速度模型, 利用中小角度的地震数据获得具有较高分辨率的阻抗模型。

1.3 速度-密度模式下双参数反演

将前述由阻抗-速度模式下反演得到的低波数速度模型作为本节中速度的初始模型, 实现反演参数的衔接。仍采用包含CPML边界的自伴随的速度-密度一阶速度-应力方程进行波场模拟:

$ \left\{ {\begin{array}{*{20}{l}} {\frac{1}{{v_{\rm{P}}^2{\rho ^2}}}\frac{{\partial p}}{{\partial t}} = \frac{1}{\rho }\left[ {\left( {\frac{1}{{{\chi _x}}}\frac{\partial }{{\partial x}} + {\varphi _x}} \right){v_x} + \left( {\frac{1}{{{\chi _y}}}\frac{\partial }{{\partial y}} + {\varphi _y}} \right) \cdot } \right.}\\ {\quad \quad \quad \quad {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left. {{v_y} + \left( {\frac{1}{{{\chi _z}}}\frac{\partial }{{\partial z}} + {\varphi _z}} \right){v_z}} \right]}\\ {\frac{{\partial {v_x}}}{{\partial t}} = \frac{1}{\rho }\left( {\frac{1}{{{\chi _x}}}\frac{\partial }{{\partial x}} + {\varphi _x}} \right)p}\\ {\frac{{\partial {v_y}}}{{\partial t}} = \frac{1}{\rho }\left( {\frac{1}{{{\chi _y}}}\frac{\partial }{{\partial y}} + {\varphi _y}} \right)p}\\ {\frac{{\partial {v_z}}}{{\partial t}} = \frac{1}{\rho }\left( {\frac{1}{{{\chi _z}}}\frac{\partial }{{\partial z}} + {\varphi _z}} \right)p} \end{array}} \right. $

(8)

同样, 根据模拟数据和观测数据波场误差最小建立(9)式目标泛函:

$ {E_2} = \frac{1}{2}|{F_{\rm{R}}}[{v_{\rm{P}}},\rho ] - {d_{{\rm{obs}}}}{|^2} $

(9)

式中:F_R[v_P, ρ]表示模拟波场F[v_P, ρ]在检波点处接收到的波场。同样采用伴随状态法可得到速度-密度模式下速度和密度的更新梯度:

$ {\rm{gra}}{{\rm{d}}_v}{E_2}[{v_{\rm{P}}},\rho ] = - \langle q,{\rm{D}}{F_v}[{v_{\rm{P}}},\rho ]\rangle $

(10a)

$ {\rm{gra}}{{\rm{d}}_\rho }{E_2}[{v_{\rm{P}}},\rho ] = - \langle q,{\rm{D}}{F_\rho }[{v_{\rm{P}}},\rho ]\rangle $

(10b)

式中:DF_v[v_P, ρ], DF_ρ[v_P, ρ]分别表示速度-密度方程对速度和密度的一阶偏导数。

$ {\rm{D}}{F_v}[{v_{\rm{P}}},\rho ] = - \frac{2}{{v_{\rm{p}}^3\rho }}\frac{{\partial p}}{{\partial t}} $

(11a)

$ {\rm{D}}{F_\rho }[{I_{\rm{p}}},{v_{\rm{P}}}] = - \frac{1}{{v_{\rm{P}}^2{\rho ^2}}}\frac{{\partial p}}{{\partial t}} + \left( {\frac{{\partial {v_x}}}{{\partial t}} + \frac{{\partial {v_y}}}{{\partial t}} + \frac{{\partial {v_z}}}{{\partial t}}} \right) $

(11b)

则更新后的速度和密度分别为:

$ {v_{2,k + 1}} = {v_{2,k}} + \beta {\rm{gra}}{{\rm{d}}_v}{E_2}[{v_{\rm{P}}},\rho ] $

(12a)

$ {\rho _{k + 1}} = {\rho _k} + \beta {\rm{gra}}{{\rm{d}}_\rho }{E_2}[{v_{\rm{P}}},\rho ] $

(12b)

式中:v_{2, k}中下标2表示在速度-密度模式中的速度; β为迭代步长。

在阻抗-速度模式的基础上利用速度-密度模式进一步反演得到速度中的高波数成分, 同时以阻抗-速度模式下反演得到的低波数的速度作为初始速度模型, 提高了速度和密度反演的稳定性。此外, 只对2种参数化模式下只进行速度反演时, 可以实现从低波数速度到高波数速度的多尺度反演^[11]。

2 模型测试

选择构造复杂和速度横向变化剧烈的Overthrust模型数据进行速度、密度、阻抗三参数全波形反演。图 2至图 7分别给出了二维速度、密度和阻抗的真实模型和初始模型。为了提高计算效率, 我们对原始模型(二维线801个CDP点, 深度187采样点)进行重采样, 新的模型数据有200个CDP点, CDP间距都是20 m, 深度方向186个点, 采样间隔为15 m。初始模型是通过对真实模型进行平滑得到。

模型深度500 m、CDP100左右处存在一个高速层, 速度为4 540 m/s, 其上部地层速度为3 300 m/s、下部地层速度3 500 m/s, 厚度约100 m, 我们标记该层为高速层a。在其正下方深度为770 m处存在另一厚度约为40 m的更薄的高速层, 速度约为4 800 m/s, 标记该层为高速层b。对于高速层a可分辨的地震数据频率最低为5.7 Hz, 对于高速层b可分辨的地震数据频率最低为15 Hz, 该分辨频率是在地震波垂直入射的条件下基于可分辨厚度为λ/4(λ为地震波长)得到的。

首先利用主频为10 Hz地震数据进行阻抗-速度以及速度-密度反演, 结果见图 8和图 9。在阻抗-速度模式中速度主要利用大偏移距数据进行低波数的速度更新, 阻抗则可以利用小偏移距数据进行阻抗更新, 明显看出阻抗模型的分辨率远高于速度模型, 两个高速层在速度模型上没有体现出来(图 9)。在深度250~500 m, CDP130附近的断层构造在阻抗中得到了较好的刻画(图 8)。

进一步, 利用主频为10 Hz的地震数据进行速度-密度反演。在速度-密度模式中由于地震数据中的小角度地震数据参与了速度反演, 浅层的断层构造更加清晰, 高速层b得到了进一步反演, 而高速层a变化并不明显(图 10)。但是全波形反演得到的密度反演结果与真实构造特征相差甚远, 密度反演不稳定(图 11)。同时基于第一步阻抗-速度反演的结果, 利用阻抗、速度、密度三者之间的关系经过数学运算得到密度模型(图 12), 由于采用除法运算, 密度更不稳定。

针对主频为10 Hz的地震数据分辨率不足问题, 本文利用主频为30 Hz的地震数据开展阻抗-速度以及速度-密度反演, 结果见图 13和图 14。此时在阻抗模型中2个高速层和浅层断层都得到非常清晰的刻画(图 13), 速度虽然仍以低波数为主, 但主要的构造特征得以充分体现(图 14)。进一步进行速度-密度反演, 结果见图 15和图 16。此时速度模型中2个高速层更加接近真实模型, 构造起伏形态与真实模型完全一致, 浅层的断层边界清晰, 断面干脆, 满足了精细地震勘探对速度的要求(图 15)。但是当速度和密度同时反演时, 高频时密度的反演更加不稳定, 特别是2个高速层的影响导致密度反演结果完全失真(图 16)。在速度准确的情况下, 通过速度-密度方程可得到较为可靠的密度反演结果(图 17)。在该种条件下, 本文在利用速度-密度模式进行反演时, 采用顺序反演策略, 首先得到精度较高的速度模型(图 15), 并以此作为初始模型进行密度反演。同时在密度反演中通过加入加德纳经验公式对密度值的阈值进行约束, 从而得到更加稳健的密度模型(图 18)。图 19对比了CDP100处约束前、后密度反演的结果, 可以明显看出, 经过约束后的密度反演更加稳定。

3 结论与认识

本文提出了一种三参数全波形反演方法, 首先利用阻抗-速度方程实现低波数速度模型和阻抗模型的全波形反演, 然后通过方程转换再利用速度-密度方程实现高波数速度反演和密度反演。该方法以低波数到高波数速度为衔接, 将2种双参数全波形反演结合实现声波方程下的三参数全波形反演。将三参数全波形反演分解为2个双参数全波形反演, 体现出两方面的优势。

1) 常规声波方程大多采用密度固定只用速度表示的单参数方程或双参数方程, 方程中参数个数小于需要反演参数个数, 常规的三参数同时反演不可实现, 本方法通过方程转换的形式实现了基于声波方程的三参数直接反演;

2) 将三参数反演分解为双参数反演, 弱化了反演过程中多参数之间的耦合性, 通过递进式的反演策略提高了反演结果的可靠性。同时将阻抗-速度模式中得到的低波数的速度作为速度-密度模式下的初始速度模型, 进一步提高了速度和密度反演的稳定性。

基于复杂的Overthrust推覆体模型对主频10 Hz和30 Hz地震数据进行反演, 得出的速度和阻抗构造起伏形态与真实模型吻合, 特别是断层边界刻画非常清晰、断面干脆, 满足了精细地震勘探对速度和阻抗的要求。主频越高反演结果分辨率越高, 对波动方程的稳定性要求也提高, 计算量增大, 在实际资料反演中可根据目的层厚度确定合适的主频。此外本文通过加入加德纳经验公式对密度的阈值进行约束提高了密度反演的稳定性, 但对于实际数据已知的信息与真实地下情况相差甚远, 经验公式约束方法的实用效果会远不如模型数据, 这也是今后需要解决的问题之一。