随着地震勘探采集技术的不断发展, 采集到的多波多分量地震数据越来越多。传统多分量地震数据处理方法是将水平分量当作转换波进行成像, 而垂直分量则是采用纵波的成像方法来处理。此种方法没有考虑地震波的矢量特征, 且成像精度主要取决于波场分离的精度, 因而当波场分离不彻底时, 残存的非本型能量波会严重影响成像质量。因此, 研究基于弹性波的成像算法具有重要意义。

作为处理强横向变速区域复杂构造成像的关键技术, 叠前深度偏移在近些年取得了快速发展。弹性波偏移算法依据理论的不同主要分为两类: 弹性波Kirchhoff偏移和弹性波逆时偏移。弹性波Kirchhoff偏移基于射线理论, 具有较高的计算效率和灵活性, 同时能够面向目标成像。PAO等^[1]首先推导了非均匀介质中的弹性波Kirchhoff积分。KUO等^[2]基于格林函数的远场近似, 实现了各向同性介质中的弹性波Kirchhoff偏移。SENA等^[3]将弹性波Kirchhoff偏移推广到各向异性介质。DRUZHININ^[4]采用极化滤波进行波形分离, 实现了解耦的弹性波Kirchhoff偏移。DU等^[5]和杜启振等^[6]研究了角度域弹性波Kirchhoff偏移与速度分析方法。弹性波逆时偏移基于全波方程, 对弹性波正反向延拓波场采用相应的成像条件进行成像, 不受成像角度的限制, 具有很高的成像精度。CHANG等^[7-8]首先利用有限差分将逆时偏移应用于多分量地震数据, 并将其推广到三维介质。YAN等^[9]实现了角度域各向异性弹性波逆时偏移。弹性波逆时偏移不仅计算量和存储量较大, 而且对速度场的精确性要求也较高; 弹性波Kirchhoff偏移计算效率较高, 但其不能处理多次波至和陡倾角问题, 成像精度相对较低。

高斯束偏移作为改进的Kirchhoff偏移, 不仅继承了射线类偏移灵活、高效的特性, 而且克服了Kirchhoff偏移固有的缺陷(焦散区、阴影区等问题)^[10-14]。POPOV等^[15]将高斯束和逆时偏移相结合, 提出了采用高斯束构建格林函数的逆时偏移算法。黄建平等^[16]在Kirchhoff积分的基础上, 实现了格林函数高斯束逆时偏移。张凯等^[17]和肖建恩等^[18]引入基于弹性参数的各向异性射线追踪体系将高斯束逆时偏移推广到各向异性介质。毕丽飞等^[19]将高斯束逆时偏移的思想应用于弹性波成像, 研究了弹性多波高斯束逆时偏移。

各向异性广泛存在于地下介质中, 地震波在各向异性介质中传播主要表现为速度频散与振幅衰减。随着采集到的宽方位角和长偏移距数据越来越多, 各向异性对偏移成像的影响变得不可忽略, 因此在地震偏移成像时, 需要根据观测系统的特性, 合理考虑速度各向异性的影响。本文基于ZHU等^[20-21]的研究, 发展了一种基于相速度的各向异性弹性波射线追踪算法, 并将该算法应用于弹性波高斯束逆时偏移, 实现了弹性波各向异性高斯束逆时偏移。模型试算证明了本文方法的成像效果及计算效率优势。

1 基本原理

弹性波各向异性高斯束逆时偏移首先通过各向异性弹性波射线追踪来计算弹性波高斯束, 然后采用弹性波高斯束来表征格林函数, 进而实现正向和反向波场的延拓, 最后计算正反向延拓波场的互相关获取成像值。

1.1 各向异性弹性波射线追踪

弹性波高斯束的计算主要通过各向异性弹性波射线追踪来实现。PP波射线追踪原理和声波一样; PS波则是在震源处利用P波进行射线追踪, 检波点处则需要用S波进行射线追踪。

1.1.1 运动学射线追踪

运动学射线追踪主要用于求解中心射线的路径和走时信息。ČERVENÝ^[22]推导了基于弹性参数的各向异性运动学射线追踪方程组, 但推导过程需在每一步射线追踪时对Christoffel方程的特征值进行求解, 计算量巨大。为了解决这个问题, ZHU等^[20-21]基于相速度对各向异性运动学射线追踪方程组进行了重新定义。本文将基于相速度的各向异性追踪体系^[20-21]推广到各向异性弹性波:

$ \left\{\begin{array}{l} \frac{\mathrm{d} x_{j}}{\mathrm{~d} \tau}=V_{\mathrm{P}_{j}} \\ \frac{\mathrm{d} p_{P_{j}}}{\mathrm{~d} \tau}=-\frac{1}{v_{\mathrm{P}}} \frac{\partial v_{\mathrm{P}}}{\partial x_{j}} \\ \frac{\mathrm{d} x_{k}}{\mathrm{~d} \tau}=V_{\mathrm{S}_{k}} \\ \frac{\mathrm{d} p_{S_{k}}}{\mathrm{~d} \tau}=-\frac{1}{v_{\mathrm{S}}} \frac{\partial v_{\mathrm{S}}}{\partial x_{k}} \end{array}\right. $

(1)

式中: V_P和V_S分别为P波和S波的群速度; j, k分别代表x和z方向, 其中, x, z分别代表水平方向和垂直方向; τ代表中心射线旅行时; v_P和v_S分别为P波和S波的相速度。

$ \left\{\begin{array}{l} v_{\mathrm{P}}=v_{\mathrm{P} 0}\left(1+\delta \sin ^{2} \theta \cos ^{2} \theta+\varepsilon \sin ^{4} \theta\right) \\ v_{\mathrm{S}}=v_{\mathrm{S} 0}\left(1+\sigma \sin ^{2} \theta \cos ^{2} \theta\right) \end{array}\right. $

(2)

式中: ε和δ为Thomosen参数; v_P0和v_S0分别代表P波和S波的垂向速度; θ为相速度角; σ=v_S0²/v_P0²(ε－δ)。

1.1.2 动力学射线追踪

在获得中心射线走时和路径后, 需要通过动力学射线追踪来计算高斯束的振幅和相位信息。对于各向异性介质, 射线中心坐标系不再正交, 此时需要引入一个沿射线路径的权值变量来消除非正交性带来的误差。HANYGA^[23]推导了基于弹性参数的各向异性的动力学射线追踪方程组, 但其求解比较复杂。ZHU等^[20-21]推导了基于相速度的各向异性动力学射线追踪方程, 本文在此基础上给出射线中心坐标系(y_M, y_N, τ)下的各向异性弹性波动力学射线追踪方程:

$ \left\{\begin{array}{l} \frac{\mathrm{d} P_{P M}}{\mathrm{~d} \tau}=-A_{M N} Q_{P N}-B_{M N} P_{P N} \\ \frac{\mathrm{d} Q_{P N}}{\mathrm{~d} \tau}=C_{M N} Q_{P N}+D_{M N} P_{P N} \end{array}\right. $

(3)

$ \left\{\begin{array}{l} \frac{\mathrm{d} P_{S N}}{\mathrm{~d} \tau}=-A_{M N}^{\prime} Q_{S N}-B_{M N}^{\prime} P_{S N} \\ \frac{\mathrm{d} Q_{S M}}{\mathrm{~d} \tau}=C_{M N}^{\prime} Q_{S N}+D_{M N}^{\prime} P_{S N} \end{array}\right. $

(4)

式中: P和Q为各向异性动力学射线追踪参数; A_MN, A′_MN, B_MN, B′_MN, C_MN, C′_MN, D_MN, D′_MN为相关系数。

$ \left\{\begin{array}{l} A_{M N}=\frac{1}{v_{P}} \frac{\partial^{2} v_{P}}{\partial y_{M} \partial y_{N}}, B_{M N}=\frac{\partial^{2} \ln v_{P}}{\partial y_{M} \partial q_{N}} \\ C_{M N}=\frac{\partial^{2} \ln v_{P}}{\partial y_{N} \partial q_{M}}, D_{M N}=\frac{\partial V_{P N}}{\partial q_{M}} \end{array}\right. $

(5)

$ \left\{\begin{array}{l} A_{M N}^{\prime}=\frac{1}{v_{S}} \frac{\partial^{2} v_{S}}{\partial y_{M} \partial_{N}}, B_{M N}^{\prime}=\frac{\partial^{2} \ln v_{S}}{\partial y_{M} \partial_{N}} \\ C_{M N}^{\prime}=\frac{\partial^{2} \ln v_{S}}{\partial y_{N} \partial q_{M}}, D_{M N}^{\prime}=\frac{\partial V_{S N}}{\partial q_{M}} \end{array}\right. $

(6)

式中: V_PN和V_SN分别为P波和S波群速度在射线中心坐标系下的分量; q_M=∂τ/∂y_M, q_N=∂τ/∂y_N。

1.2 正反向弹性波场延拓

ČERVENÝ等^[24]给出了二维抛物线波动方程的求解过程, 在此基础上可以得到P波和S波的高斯束表达式:

$ \begin{gathered} \boldsymbol{u}_{\mathrm{P}}(s, n, \omega, t)=\frac{\varPsi_{\mathrm{P}}}{\sqrt{v_{\mathrm{P}}(s) \rho(s) Q(s)}}\left[\boldsymbol{t}+\frac{n P(s) v_{\mathrm{P}}(s)}{Q(s)} \boldsymbol{n}\right] \times \\ \exp \left\{-\mathrm{i} \omega\left[t-\int\limits_{0}^{s} \frac{\mathrm{d} s}{v_{\mathrm{P}}(s)}-\frac{1}{2} \frac{P(s)}{Q(s)} n^{2}\right]\right\} \end{gathered} $

(7)

$ \begin{gathered} \boldsymbol{u}_{S}(s, n, \omega, t)=\frac{\varPsi_{\mathrm{S}}}{\sqrt{v_{\mathrm{S}}(s) \rho(s) Q(s)}}\left[\boldsymbol{n}-\frac{n P(s) v_{\mathrm{S}}(s)}{Q(s)} \boldsymbol{t}\right] \times \\ \exp \left\{-\mathrm{i} \omega\left[t-\int\limits_{0}^{s} \frac{\mathrm{d} s}{v_{\mathrm{S}}(s)}-\frac{1}{2} \frac{P(s)}{Q(s)} n^{2}\right]\right\} \end{gathered} $

(8)

式中: s为任一点在射线上投影点到初始点的距离; n为任一点到射线的法线距离; ω为圆频率; t为旅行时; P(s)和Q(s)为各向异性动力学射线追踪参数; n和t为射线中心坐标系的基矢量; ρ为密度; Ψ_P和Ψ_S为权系数。

$ \left\{\begin{array}{l} \varPsi_{\mathrm{P}}=\frac{\mathrm{i}}{4 {\rm{ \mathsf{ π} }}\left[v_{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{s}}\right)\right]^{2}} \sqrt{\frac{\omega_{\mathrm{r}} w_{0}^{2}}{\rho\left(\boldsymbol{x}_{\mathrm{s}}\right)}} \\ \varPsi_{\mathrm{S}}=\frac{\mathrm{i}}{4 {\rm{ \mathsf{ π} }}\left[v_{\mathrm{S}}(\boldsymbol{x} \mathrm{s})\right]^{2}} \sqrt{\frac{\omega_{\mathrm{r}} w_{0}^{2}}{\rho\left(\boldsymbol{x}_{\mathrm{s}}\right)}} \end{array}\right. $

(9)

式中: ω_r为参考频率; w₀为高斯束初始宽度; x_s为震源位置。

对于各向异性弹性波高斯束逆时偏移, 通过P波高斯束位移矢量来表征震源波场, 可以得到震源处正向延拓波场:

$ \begin{gathered} U_{F}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{s}} ; t_{0}\right)=\frac{\mathrm{i}}{4 v_{\mathrm{P}}^{2}\left(\boldsymbol{x}_{\mathrm{s}}\right)} \mathrm{Re} \int\limits_{0}^{\infty} \mathrm{d} \omega \mathrm{e}^{-\mathrm{i} \omega t} f(\omega) \sqrt{\frac{\omega_{\mathrm{r}} w_{0}^{2}}{\rho\left(\boldsymbol{x}_{\mathrm{s}}\right)}} \times\\ \int \frac{\mathrm{d} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{s}}\right)}{p^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{s}}\right)} A^{\mathrm{P*}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{s}}\right) \exp \left[-\mathrm{i} \omega T^{\mathrm{P*}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{s}}\right)\right] \end{gathered} $

(10)

式中: A和T分别为高斯束的复值振幅和走时; p_x, p_z为慢度; 符号*表示共轭。

根据PAO等^[1]和毕丽飞等^[19]的研究, 接收点x_r处的反向延拓波场P波和S波位移波场可以表示为:

$ \begin{aligned} u_{m}^{\mathrm{P}}&\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}} ; t_{0}\right)=-\int\limits_{t_{0}}^{\mathrm{T}} \mathrm{d} t \frac{\sqrt{\omega_{\mathrm{r}} w_{0}^{2}}}{4 {\rm{ \mathsf{ π} }}} \int_{\partial \varOmega}\left[u_{x}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right) \cdot\right. \\ &\left.W_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+u_{z}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right) W_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)\right] \mathrm{d} x{ }_{\mathrm{r}} \int \boldsymbol{e}_{m}^{\mathrm{P}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right)\cdot \\ &A^{\mathrm{P} *}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right) \exp \left[-\mathrm{i} \omega T^{\mathrm{P} *}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right)\right] \frac{\mathrm{d} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)}{p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)} \end{aligned} $

(11)

$ \begin{aligned} u_{m}^{\mathrm{S}}&\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}} ; t_{0}\right)=-\int\limits_{t_{0}}^{\mathrm{T}} \mathrm{d} t \frac{\sqrt{\omega_{\mathrm{r}} w_{0}^{2}}}{4 {\rm{ \mathsf{ π} }}} \int_{\partial \varOmega}\left[u_{x}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right)\cdot\right. \\ &\left.W_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+u_{z}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right) W_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)\right] \mathrm{d} x_{\mathrm{r}} \int \boldsymbol{e}_{m}^{\mathrm{S}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right) \cdot\\ &A^{\mathrm{S} *}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right) \exp \left[-\mathrm{i} \omega T^{\mathrm{S} *}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right)\right] \frac{\mathrm{d} p_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)}{p_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)} \end{aligned} $

(12)

式中: u_x(x_r; t)和u_z(x_r; t)分别为在x_r处接收到的地震记录水平分量和垂直分量; e_m^P(x; x_r)和e_m^S(x; x_r)|_{m=1, 2}分别为P波和S波的极性矢量; W_x^P(x_r), W_z^P(x_r), W_x^S(x_r), W_z^S(x_r)为相关权系数。

$ \begin{gathered} W_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)=C_{15} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{55} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+ \\ C_{55} p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{35} p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) \\ W_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)=C_{13} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{35} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+ \\ C_{35} p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{33} p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right) \\ W_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)=C_{15} p_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{55} p_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+ \\ C_{55} p_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{35} p_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) \\ W_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)=C_{13} p_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{35} p_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+ \\ C_{35} p_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+C_{33} p_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) e_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right) \end{gathered} $

(13)

式中: C为弹性参数。

1.3 弹性波各向异性成像公式

根据反射波成像准则, 成像值可以通过震源波场和接收点反向延拓波场的互相关来求取, PP波和PS波成像公式分别为:

$ \begin{gathered} I^{\mathrm{PP}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{s}}\right)=-\frac{\mathrm{i} \omega_{\mathrm{r}} w_{0}^{2}}{16 {\rm{ \mathsf{ π} }} v_{\mathrm{P}}^{2}\left(\boldsymbol{x}_{\mathrm{s}}\right)} \int \mathrm{d} t_{0} \int_{\partial \varOmega} \sqrt{\frac{\rho\left(\boldsymbol{x}_{\mathrm{r}}\right)}{\rho\left(\boldsymbol{x}_{\mathrm{s}}\right)}} \mathrm{d} \boldsymbol{x}_{\mathrm{r}} \cdot \\ \iint \frac{\mathrm{d} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{s}}\right) \mathrm{d} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)}{p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{s}}\right) p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)} u_{\mathrm{GB} z}^{\mathrm{P}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{s}}\right) u_{\mathrm{GB} z}^{\mathrm{P}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right) \cdot \\ {\left[u_{x}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right) W_{1}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+u_{z}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right) W_{2}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{r}}\right)\right]} \end{gathered} $

(14)

$ \begin{gathered} I^{\mathrm{PS}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{s}}\right)=-\frac{\mathrm{i} \omega_{\mathrm{r}} w_{0}^{2}}{16 {\rm{ \mathsf{ π} }} v_{\mathrm{P}}^{2}\left(\boldsymbol{x}_{\mathrm{s}}\right)} \int \mathrm{d} t_{0} \int_{\partial \varOmega} \sqrt{\frac{\rho\left(\boldsymbol{x}_{\mathrm{r}}\right)}{\rho\left(\boldsymbol{x}_{\mathrm{s}}\right)}} \mathrm{d} \boldsymbol{x}_{\mathrm{r}} \cdot\\ \iint \frac{\mathrm{d} p_{x}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{s}}\right) \mathrm{d} p_{x}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)}{p_{z}^{\mathrm{P}}\left(\boldsymbol{x}_{\mathrm{s}}\right) p_{z}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)} \operatorname{sgn} \alpha u_{\mathrm{GB} z}^{\mathrm{P}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{s}}\right) u_{\mathrm{GB} x}^{\mathrm{S}}\left(\boldsymbol{x} ; \boldsymbol{x}_{\mathrm{r}}\right) \cdot\\ \left[u_{x}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right) W_{1}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)+u_{z}\left(\boldsymbol{x}_{\mathrm{r}} ; t\right) W_{2}^{\mathrm{S}}\left(\boldsymbol{x}_{\mathrm{r}}\right)\right] \end{gathered} $

(15)

式中: I^PP(x; x_s)和I^PS(x; x_s)分别为PP波和PS波的单炮成像值; u_GBz^P(x; x_s)和u_GBx^S(x; x_r)分别为P波高斯束和S波高斯束; sgnα为符号函数, 用于校正转换波成像中的极性反转问题。

1.4 技术实现流程

各向异性弹性波高斯束逆时偏移具体实现流程包括以下5步:

1) 读入速度场和炮记录以及相关各向异性参数场;

2) 通过各向异性弹性波运动学和动力学射线追踪求取中心射线的运动学和动力学信息;

3) 在震源处利用(10)式计算正向传播波场, 同时在接收点处利用(11)式或(12)式构建反向延拓波场;

4) 利用(14)式或(15)式计算震源处正向传播波场和接收点处反向延拓波场的互相关, 得到单炮成像值;

5) 所有炮计算完成后, 将单炮成像值进行叠加, 得到最终的成像结果。

2 模型试算 2.1 复杂构造模型

采用VTI复杂构造模型测试本文方法的可行性和有效性。模型由一个平层、两个断层和底部两个起伏界面组成, 速度场和各向异性参数场如图 1所示。速度场模型横向宽度为16 010 m, 纵向深度为4 000 m, 网格间距均为10 m。本文通过各向异性高斯束正演方法来合成地震记录, 震源子波为雷克子波, 主频为25 Hz。合成地震记录共281炮, 每炮201道接收, 炮间距为50 m, 道间距为10 m。时间采样点数为4 001, 采样间隔为1 ms。图 2为单炮记录的X分量和Z分量, 可以看出, 在X分量中存在极性反转的情况。

采用图 1所示复杂构造模型进行试算。图 3a和图 3b分别是采用各向同性PP波和PS波进行偏移成像的结果, 可以看出, 由于忽略了各向异性的影响, 偏移剖面中反射界面不能准确归位, 而且存在大量的成像噪声, 同相轴的聚焦性和连续性较差。而在采用本文方法对PP波和PS波进行偏移成像的结果(图 3c, 图 3d)中, 对各向异性构造进行了精确的刻画, 成像噪声也得到了压制, 同相轴具有更好的连续性和聚焦性, 剖面整体信噪比更高。对比图 3c和图 3d可以看出, PS波的成像范围更宽, 这是由于PS波传播速度较慢, 传播角度大, 包含更为丰富的地下介质信息。

2.2 断层模型

采用断层模型测试本文方法在各向异性TTI介质中的适用性。模型网格大小为1 001×301, 网格间距均为10 m。模型包含两个平层和两个断层, 图 4是纵、横波速度场和各向异性参数场。合成地震数据由各向异性高斯束正演算法获得, 震源子波为雷克子波, 主频为25 Hz。地震记录共201炮, 每炮201道接收, 炮间距和道间距分别为40 m和10 m。时间采样长度为2.5 s, 采样间隔为1 ms, 单炮记录的X分量和Z分量如图 5所示。

模型试算结果如图 6所示, 对比可以看出, 采用各向同性算法处理各向异性数据会导致结果中存在大量的成像噪声和绕射波, 成像位置也发生了偏离, 同时断层和底部平层同相轴的能量聚焦性较差(如图 6a和图 6b中的红线、箭头和红框所示), 剖面信噪比较低。而在基于弹性参数的各向异性算法(图 6c和图 6d)和本文方法的PP波和PS波偏移结果(图 6e和图 6f)中, 构造位置准确, 反射层得到了清晰的刻画, 而且同相轴的聚焦性和连续性更好, 成像质量得到明显提高。

在相同的软件和硬件配置下, 各向同性算法、基于弹性参数各向异性算法和本文算法的平均单炮计算时间分别为384.6, 612.8, 468.5 s, 可以看出, 在不影响成像精度的前提下, 相比于传统基于弹性参数的各向异性算法, 本文方法在计算效率方面具有较大优势。

3 结论与认识

弹性波高斯束逆时偏移采用弹性波动力学高斯束来表征格林函数, 进而实现弹性波场的延拓, 最后通过正反向延拓波场的互相关来获取成像值。不仅继承了弹性波高斯束偏移的高效性和灵活性, 而且具有接近弹性波逆时偏移的成像精度, 能够实现面向目标成像。本文引入基于相速度的各向异性弹性波射线追踪算法, 实现了弹性波各向异性高斯束逆时偏移。模型试算表明, 本文方法能够对各向异性构造进行准确成像, 而且相比于传统基于弹性参数的各向异性算法, 在不影响成像精度的前提下, 在计算效率方面更有优势。