一种基于辅助坐标系的起伏地表波形反演方法

摘要

本发明公开了一种基于辅助坐标系的起伏地表波形反演方法，属于石油地球物理勘探技术领域，本发明首先采用辅助坐标系下的早至波波形反演方法对曲网格区域进行反演，以得到更新的近地表速度，再采用辅助坐标系下的全波形反演全局速度场进行速度更新，以克服近地表速度不准确对深部速度场反演的影响；采用时间域多尺度反演方法从低频到高频进行速度反演，以克服波形反演方法对初始速度模型的依赖性。本发明能够很好地反演剧烈起伏地表的速度场，为高精度成像方法提供准确的偏移速度场。

说明书

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种基于辅助坐标系的起伏地表波形反演 方法。

背景技术

剧烈起伏地表对地震资料采集、处理带来了严重的影响，地球物理研究人员发展了一系 列方法来克服这个问题。目前，针对起伏地表的处理主要有两种策略：一是对表层波场进行 校正，二是基于起伏地表进行深度偏移成像。但复杂条件下静校正量难以准确计算，且静校 正无法彻底消除地表起伏对地震波场造成的扭曲，因此直接基于起伏地表的深度偏移成像方 法逐渐成为研究的热点。

由于全波形反演高精度以及高分辨率的特点，使其成为速度建模的一种有力工具，逐渐 成为研究的热点。全波形反演是一个非线性数据拟合的过程，通过减少观测数据与预测数据 的之间的差值来更新参数模型，这个过程以迭代的方式重复下去，直到数据差值足够小为止。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种基于辅助坐标系的起伏地表波 形反演方法，设计合理，具有良好的推广价值。

为了实现上述目的，本发明采用如下技术方案：

一种基于辅助坐标系的起伏地表波形反演方法，按照如下步骤进行：

步骤1：输入初始全局速度场、常规叠前炮记录、起伏高程及震源子波，并建立观测系 统；

步骤2：根据初始全局速度场及起伏高程进行网格剖分，将近地表附近的网格剖分成曲 网格，深层的网格剖分成矩形网格；

步骤3：将初始全局速度场变换到辅助坐标系下的矩形网格，采用下式所示的变换格式：

$\left(\begin{array}{c}x(\xi ,\eta )=\xi \\ z(\xi ,\eta )=\frac{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}{{\eta }_{i-1}\left(\xi \right)-{\eta }_i\left(\xi \right)}(\eta -{\eta }_i)+z_i\left(\xi \right)\end{array}\right);$

其中，x和z表示笛卡尔坐标系下的横纵坐标；ξ和η表示辅助坐标系下的横纵坐标；z_i-1(ξ) 和z_i(ξ)是笛卡尔坐标系下第i层顶界面、底界面的高程，定义最深层为零；η_i-1(ξ)和η_i(ξ)是辅 助坐标系第i层顶界面、底界面的纵向采样点数，定义最深层为零；

步骤4：对常规叠前炮记录划分时窗，在辅助坐标系下对近地表附近的曲网格区域速度 场应用早至波波形反演，更新近地表的速度场，早至波波形反演的梯度方向如下：

$g\left(v\right)=\frac{2}{v^3}\underset{x_s}{\Sigma }\underset{t=0}{\overset{T_e}{\Sigma }}\frac{\partial p}{\partial t}·p^{*};$

其中，g为梯度；x_s为炮点坐标；v为速度场的速度；p为声压；p^*是残差波场的反向传 播；t为时间；T_e表示早至波时窗；

步骤5：判断使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮 记录之差是否满足误差条件；

若：判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前 炮记录之差满足误差条件，则近地表速度场更新完成，然后执行步骤6；

或判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮 记录之差不满足误差条件，则执行步骤4；

步骤6：将步骤1中输入的常规叠前炮记录分解成不同主频的多尺度炮记录，采用如下 式所示的分解公式：

$f\left(\omega \right)=\frac{W_t\left(\omega \right)W_o^{*}\left(\omega \right)}{|W_o\left(\omega \right){|}^2+{ϵ}^2}$

其中，f是维纳滤波器；W_o表示初始炮记录子波；W_t是生成的炮记录子波；ω是角频 率；ε为一个小数；*表示共轭转置；

步骤7：将更新过近地表速度的全局速度场作为初始速度场，在辅助坐标系下应用全波 形反演从低频到高频更新全局速度场，全波形反演的梯度方向如下：

$g\left(v\right)=\frac{2}{v^3}\underset{x_s}{\Sigma }\underset{t=0}{\overset{T_{\max }}{\Sigma }}\underset{f=f_1}{\overset{f_{\max }}{\Sigma }}\frac{\partial p}{\partial t}·p_f^{*}$

其中，f表示主频；f₁和f_max为多尺度分解的最低频率和最高频率；T_max表示常规叠前 炮记录的最大记录时间；表示主频是f的残差波场反传；

步骤8：判断使用全局速度场正演模拟的炮记录与常规叠前炮记录之差是否满足误差条 件；

若：判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差满足误差条件， 则全局速度场更新完成，然后执行步骤9；

或判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差不满足误差条 件，则执行步骤7；

步骤9：将更新完成的全局速度场反变换到笛卡尔坐标系下，反变换公式如下：

$\left(\begin{array}{c}\xi (x,z)=x\\ \eta (x,z)=\frac{{\eta }_{i-1}\left(\xi \right)-{\eta }_i\left(\xi \right)}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}(z-z_i)+{\eta }_i\left(\xi \right)\end{array}\right);$

步骤10：输出反演的速度场。

优选地，在步骤4中，具体包括

步骤4.1：定义目标函数：

$E=\frac{1}{2}\underset{x_s}{\Sigma }\underset{x_r}{\Sigma }\underset{t}{\Sigma }{\left(Ru\right(t,x_r,x_s)-d_{obs}(t,x_r,x_s\left)\right)}^2---\left(1\right);$

其中，u(t,x_r,x_s)代表模拟波场u＝(u,w,p)^T，其中T表示转置；R为限定算子； d_obs(t,x_r,x_s)是常规叠前炮记录；x_s和x_r表示震源点和检波点的位置坐标；t表示时间；E 为目标函数值；

步骤4.2：将目标函数变分，得到变分表达式：

$\begin{array}{c}\delta E\left(v\right)=\frac{1}{2}\delta <(Ru-d_{obs}),(Ru-d_{obs})>\\ =<\delta (Ru-d_{obs}),(Ru-d_{obs})>=<\left(R\delta u\right),(Ru-d_{obs})>\end{array}---\left(2\right);$

步骤4.3：定义变换格式：

$\left(\begin{array}{c}x(\xi ,\eta )=\xi \\ z(\xi ,\eta )=\frac{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}{{\eta }_{i-1}\left(\xi \right)-{\eta }_i\left(\xi \right)}(\eta -{\eta }_i)+z_i\left(\xi \right)\end{array}\right)---\left(3\right);$

通过变换格式(3)与锁链法则得到下面的映射公式：

$\{\begin{array}{c}\frac{\partial \xi }{\partial x}=1\\ \frac{\partial \xi }{\partial z}=0\\ \frac{\partial \eta }{\partial z}=\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\\ \frac{\partial \eta }{\partial x}=\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi }\end{array}---\left(4\right);$

通过映射公式(4)，得到辅助坐标系下的一阶方程：

$\left(\begin{array}{c}\rho \frac{\partial u}{\partial t}-\frac{\partial p}{\partial \xi }-\frac{\partial p}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })=0\\ \rho \frac{\partial w}{\partial t}-\frac{\partial p}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)=0\\ \frac{1}{v^2}\frac{\partial p}{\partial t}-\rho [\frac{\partial u}{\partial \xi }+\frac{\partial u}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\\ \frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })+\frac{\partial w}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)]=f\end{array}\right)---\left(5\right);$

其中，p是声压；u和w分别是水平方向和垂直方向的质点速度；v是介质速度；f表示 震源项；

速度的扰动δv可以引起地震波场的扰动δu，δu＝(δu,δw,δp)^T，将v+δv→u+δu 代入一阶方程(5)后与一阶方程(5)相减得到如下方程：

$\left(\begin{array}{c}\rho \frac{\partial \delta u}{\partial t}-\frac{\partial \delta p}{\partial \xi }-\frac{\partial \delta p}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })=0\\ \rho \frac{\partial \delta w}{\partial t}-\frac{\partial \delta p}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)=0\\ \frac{1}{v^2}\frac{\partial \delta p}{\partial t}-\rho [\frac{\partial \delta u}{\partial \xi }+\frac{\partial \delta u}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\\ \frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })+\frac{\partial \delta w}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)]=\frac{2\delta v}{v^2}\frac{\partial p}{\partial t}\end{array}\right)---\left(6\right);$

进一步可得：

$\delta u=L\left(\frac{2\delta v}{v^3}\frac{\partial p}{\partial t}\right)---\left(7\right);$

其中，L代表正演算子；

步骤4.4：将方程(7)代入变分表达式(2)，可得：

$\delta E\left(v\right)=<R\left(L\right(\frac{2\delta v}{v^3}\frac{\partial p}{\partial t}\left)\right),(Ru-d_{obs})>=<\left(\frac{2\delta v}{v^3}\frac{\partial p}{\partial t}\right),L^{*}{R}^{*}(Ru-d_{obs})>---\left(8\right);$

其中，L^*R^*(Ru-d_obs)表示波场的逆时传播；

利用伴随状态法，L^*R^*(Ru-d_obs)可由下式求得：

$\left(\begin{array}{c}\rho \frac{\partial u^{*}}{\partial t}-\frac{\partial p^{*}}{\partial \xi }-\frac{\partial p^{*}}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}+\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })=0\\ \rho \frac{\partial w^{*}}{\partial t}-\frac{\partial p^{*}}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)=0\\ \frac{1}{v^2}\frac{\partial p^{*}}{\partial t}-\rho [\frac{\partial u^{*}}{\partial \xi }+\frac{\partial u^{*}}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\\ \frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })+\frac{\partial w^{*}}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)]=p-p_{obs}\end{array}\right)---\left(9\right);$

则L^*R^*(Ru-d_obs)＝p^*；

步骤4.5：求取目标函数对速度模型的梯度，得到早至波波形反演的梯度方向：

$g\left(v\right)=\frac{\delta E\left(v\right)}{\delta v}=\frac{2}{v^3}\underset{x_s}{\Sigma }\underset{t=0}{\overset{T_e}{\Sigma }}\frac{\partial p}{\partial t}·p^{*}---\left(10\right);$

其中，g为梯度；x_s为炮点坐标；v为速度场的速度；p为声压；p^*是残差波场的反向传 播；t为时间；T_e表示早至波时窗。

本发明所带来的有益技术效果：

本发明提出了一种基于辅助坐标系的起伏地表波形反演方法，与现有技术相比，一种基 于辅助坐标系的起伏地表波形反演方法，采用一阶方程变密度方程克服传统常密度二阶方程 在处理密度变化比较大地区的速度反演不准确的缺点；同时在起伏地表处采用曲网格剖分， 在深部区域采用矩形网格剖分，并统一变换到辅助坐标系下的矩形网格进行计算，克服起伏 地表对地震波传播的影响；首先采用辅助坐标系下的早至波波形反演方法对曲网格区域进行 反演，以得到更新的近地表速度，再采用辅助坐标系下的全波形反演全局速度场进行速度更 新，以克服近地表速度不准确对深部速度场反演的影响；采用时间域多尺度反演方法从低频 到高频进行速度反演，以克服波形反演方法对初始速度模型的依赖性。本发明能够很好地反 演剧烈起伏地表的速度场，为高精度成像方法提供准确的偏移速度场。

附图说明

图1为本发明的一种基于辅助坐标系的起伏地表波形反演方法的流程图。

图2为本发明使用的真实起伏地表速度模型。

图3为变换到辅助坐标系下的速度模型。

图4为采用统一曲网格变换到辅助坐标系下的速度模型。

图5为初始全局速度场。

图6为输入的常规叠前炮记录。

图7为采用统一曲网格正演模拟得到的炮记录。

图8为多尺度分解得到的5Hz叠前炮记录。

图9为多尺度分解得到的15Hz叠前炮记录。

图10为采用本发明方法得到的第一尺度反演速度场(主频为5Hz)。

图11为采用本发明方法得到的第二尺度反演速度场(主频为15Hz)。

图12为采用本发明方法得到的第三尺度反演速度场(主频为25Hz)。

图13为采用单尺度波形反演方法得到的速度场。

图14为采用统一曲网格得到的波形反演速度场。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

一种基于辅助坐标系的起伏地表波形反演方法的流程图(如图1所示)，包括如下步骤：

输入初始全局速度场、常规叠前炮记录、起伏高程及震源子波，并建立观测系统；

根据初始全局速度场及起伏高程进行网格剖分，近地表附近的网格剖分成曲网格，深层 的网格剖分成矩形网格；

将初始全局速度场变换到辅助坐标系下的矩形网格；

对常规叠前炮记录划分时窗，在辅助坐标系下对近地表附近的曲网格区域速度场应用早 至波波形反演，更新近地表的速度场，判断使用近地表速度场正演模拟的炮记录与划分时窗 的常规叠前炮记录的差是否满足误差条件，如果满足则近地表速度场更新完成，如果不满足 再次应用早至波波形反演更新近地表速度场；

将输入的常规叠前炮记录分解成不同主频的多尺度炮记录，将更新过近地表速度的全局 速度场作为初始速度场，在辅助坐标系下应用全波形反演从低频到高频更新全局速度场，判 断使用全局速度场正演模拟的炮记录与常规叠前炮记录的差是否满足误差条件，如果满足则 全局速度场更新完成，如果不满足再次应用全波形波形反演更新全局速度场；

将更新完成的全局速度场变换到笛卡尔坐标系下，并输出反演的速度场。

具体步骤为：

步骤1：输入初始全局速度场、常规叠前炮记录、起伏高程及震源子波，并建立观测系 统；

步骤2：根据初始全局速度场及起伏高程进行网格剖分，将近地表附近的网格剖分成曲 网格，深层的网格剖分成矩形网格；

步骤3：将初始全局速度场变换到辅助坐标系下的矩形网格，采用下式所示的变换格式：

$\left(\begin{array}{c}x(\xi ,\eta )=\xi \\ z(\xi ,\eta )=\frac{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}{{\eta }_{i-1}\left(\xi \right)-{\eta }_i\left(\xi \right)}(\eta -{\eta }_i)+z_i\left(\xi \right)\end{array}\right)---\left(1\right);$

其中，x和z表示笛卡尔坐标系下的横纵坐标；ξ和η表示辅助坐标系下的横纵坐标；z_i-1(ξ) 和z_i(ξ)是笛卡尔坐标系下第i层顶界面、底界面的高程，定义最深层为零；η_i-1(ξ)和η_i(ξ)是辅 助坐标系第i层顶界面、底界面的纵向采样点数，定义最深层为零。

通过变换格式(1)与锁链法则得到下面的映射公式：

$\left(\begin{array}{c}\frac{\partial \xi }{\partial x}=1\\ \frac{\partial \xi }{\partial z}=0\\ \frac{\partial \eta }{\partial z}=\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\\ \frac{\partial \eta }{\partial x}=\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi }\end{array}\right)---\left(2\right);$

步骤4：对常规叠前炮记录划分时窗，在辅助坐标系下对近地表附近的曲网格区域速度 场应用早至波波形反演，更新近地表的速度场，早至波波形反演的梯度方向如下：

通过映射公式(2)，得到辅助坐标系下的一阶方程：

$\left(\begin{array}{c}\rho \frac{\partial u}{\partial t}-\frac{\partial p}{\partial \xi }-\frac{\partial p}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })=0\\ \rho \frac{\partial w}{\partial t}-\frac{\partial p}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)=0\\ \frac{1}{v^2}\frac{\partial p}{\partial t}-\rho [\frac{\partial u}{\partial \xi }+\frac{\partial u}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\\ \frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })+\frac{\partial w}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)]=f\end{array}\right)---\left(3\right);$

其中，p是声压；u和w分别是水平方向和垂直方向的质点速度；v是介质速度；f表示 震源项。

对于波形反演需要构建目标函数。

本发明采用常规叠前炮记录与模拟炮记录残差的L2模作为目标函数：

$E=\frac{1}{2}\underset{x_s}{\Sigma }\underset{x_r}{\Sigma }\underset{t}{\Sigma }{\left(Ru\right(t,x_r,x_s)-d_{obs}(t,x_r,x_s\left)\right)}^2---\left(4\right);$

其中，u(t,x_r,x_s)代表模拟波场u＝(u,w,p)^T，其中T表示转置；R为限定算子； d_obs(t,x_r,x_s)是常规叠前炮记录；x_s和x_r表示震源点和检波点的位置坐标；t表示时间；E 为目标函数值；

将目标函数变分，得到变分表达式：

$\begin{array}{c}\delta E\left(v\right)=\frac{1}{2}\delta <(Ru-d_{obs}),(Ru-d_{obs})>\\ =<\delta (Ru-d_{obs}),(Ru-d_{obs})>=<\left(R\delta u\right),(Ru-d_{obs})>\end{array}---\left(5\right);$

速度的扰动δv可以引起地震波场的扰动δu，δu＝(δu,δw,δp)^T，将v+δv→u+δu代 入辅助坐标系下的一阶方程(3)中并与之相减得到如下方程：

$\left(\begin{array}{c}\rho \frac{\partial \delta u}{\partial t}-\frac{\partial \delta p}{\partial \xi }-\frac{\partial \delta p}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })=0\\ \rho \frac{\partial \delta w}{\partial t}-\frac{\partial \delta p}{\partial \eta }\left(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)=0\\ \frac{1}{v^2}\frac{\partial \delta p}{\partial t}-\rho [\frac{\partial \delta u}{\partial \xi }+\frac{\partial \delta u}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\\ \frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })+\frac{\partial \delta w}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)]=\frac{2\delta v}{v^2}\frac{\partial p}{\partial t}\end{array}\right)---\left(6\right);$

从方程(6)，我们可以得到其中L代表正演算子，方程(5)进一步 可表示为：

$\delta E\left(v\right)=<R\left(L\right(\frac{2\delta v}{v^3}\frac{\partial p}{\partial t}\left)\right),(Ru-d_{obs})>=<\left(\frac{2\delta v}{v^3}\frac{\partial p}{\partial t}\right),L^{*}{R}^{*}(Ru-d_{obs})>---\left(7\right);$

其中，L^*R^*(Ru-d_obs)表示波场的逆时传播；

利用伴随状态法，L^*R^*(Ru-d_obs)可由下式求得：

$\left(\begin{array}{c}\rho \frac{\partial u^{*}}{\partial t}-\frac{\partial p^{*}}{\partial \xi }-\frac{\partial p^{*}}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })=0\\ \rho \frac{\partial w^{*}}{\partial t}-\frac{\partial p^{*}}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)=0\\ \frac{1}{v^2}\frac{\partial p^{*}}{\partial t}-\rho [\frac{\partial u^{*}}{\partial \xi }+\frac{\partial u^{*}}{\partial \eta }(\frac{{\eta }_{i-1}-\eta }{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_i\left(\xi \right)}{\partial \xi }+\\ \frac{\eta -{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}·\frac{\partial z_{i-1}\left(\xi \right)}{\partial \xi })+\frac{\partial w^{*}}{\partial \eta }\left(\frac{{\eta }_{i-1}-{\eta }_i}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}\right)]=p-p_{obs}\end{array}\right)---\left(8\right);$

则L^*R^*(Ru-d_obs)＝p^*；

早至波波形反演的梯度方向如下：

$g\left(v\right)=\frac{\delta E\left(v\right)}{\delta v}=\frac{2}{v^3}\underset{x_s}{\Sigma }\underset{t=0}{\overset{T_e}{\Sigma }}\frac{\partial p}{\partial t}·p^{*}---\left(9\right);$

其中，g为梯度；x_s为炮点坐标；v为速度场的速度；p为声压；p^*是残差波场的反向传 播；t为时间；T_e表示早至波时窗；

步骤5：判断使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮 记录之差是否满足误差条件；

若：判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前 炮记录之差满足误差条件，则近地表速度场更新完成，然后执行步骤6；

或判断结果是使用近地表速度场正演模拟的炮记录与步骤4中划分时窗后的常规叠前炮 记录之差不满足误差条件，则执行步骤4；

步骤6：将步骤1中输入的常规叠前炮记录分解成不同主频的多尺度炮记录，采用如下 式所示的分解公式：

$f\left(\omega \right)=\frac{W_t\left(\omega \right)W_o^{*}\left(\omega \right)}{|W_o\left(\omega \right){|}^2+{ϵ}^2}---\left(10\right);$

其中，f是维纳滤波器；W_o表示初始炮记录子波；W_t是生成的炮记录子波；ω是角频 率；ε为一个小数；*表示共轭转置；

步骤7：将更新过近地表速度的全局速度场作为初始速度场，在辅助坐标系下应用全波 形反演从低频到高频更新全局速度场，全波形反演的梯度方向如下：

$g\left(v\right)=\frac{2}{v^3}\underset{x_s}{\Sigma }\underset{t=0}{\overset{T_{\max }}{\Sigma }}\underset{f=f_1}{\overset{f_{\max }}{\Sigma }}\frac{\partial p}{\partial t}·p_f^{*}---\left(11\right);$

其中，f表示主频；f₁和f_max为多尺度分解的最低频率和最高频率；T_max表示常规叠前 炮记录的最大记录时间；表示主频是f的残差波场反传；

步骤8：判断使用全局速度场正演模拟的炮记录与常规叠前炮记录之差是否满足误差条 件；

若：判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差满足误差条件， 则全局速度场更新完成，然后执行步骤9；

或判断结果是使用全局速度场正演模拟的炮记录与常规叠前炮记录之差不满足误差条 件，则执行步骤7；

步骤9：将更新完成的全局速度场反变换到笛卡尔坐标系下，反变换公式如下：

$\left(\begin{array}{c}\xi (x,z)=x\\ \eta (x,z)=\frac{{\eta }_{i-1}\left(\xi \right)-{\eta }_i\left(\xi \right)}{z_{i-1}\left(\xi \right)-z_i\left(\xi \right)}(z-z_i)+{\eta }_i\left(\xi \right)\end{array}\right)---\left(12\right);$

步骤10：输出反演的速度场。

本发明的一种基于辅助坐标系的起伏地表波形反演方法，能够很好地反演剧烈起伏地表 的速度场，为高精度成像方法提供准确的偏移速度场。

应用实验

本发明一种基于辅助坐标系的起伏地表波形反演方法，应用于复杂起伏地表模型数据。 图2为本发明使用的真实起伏地表速度模型；图3为变换到辅助坐标系下速度模型；图4为 采用统一曲网格变换到辅助坐标系下的速度模型；应用本发明的网格剖分策略对近地表附近 的网格剖分成曲网格，深层的网格剖分成矩形网格，再变换到辅助坐标系下的速度模型(图 3)与应用传统统一曲网格剖分，再变换到辅助坐标系下的速度模型(图4)相比，本发明的 结果更好地保存了深部构造的原始形态。图5为输入初始全局速度场，图6为输入的常规叠 前炮记录，也是采用本发明网格剖分策略正演模拟得到的炮记录，与采用统一曲网格剖分得 到的叠前炮记录(图7)对比，本发明方法得到的炮记录同相轴更加清晰，而且虚假的绕射 波和散射波信息更少，因为波形反演的基础是正演模拟，通过正演模拟的结果对比也能体现 本发明的优势。将输入的常规叠前炮记录分解成不同主频的多尺度炮记录，图8为主频为5Hz 的炮记录，图9为主频为15Hz的炮记录。在辅助坐标系下应用全波形反演从低频到高频更新 全局速度场，图10为采用本发明方法得到的第一尺度反演速度场(主频为5Hz)，图11采用 本发明方法得到的第二尺度反演速度场(主频为15Hz)，图12采用本发明方法得到的第三尺 度反演速度场(主频为25Hz)，从三图的对比可以看出，采用第一尺度得到的反演速度场， 低频信息得到了很好地恢复，随着反演频率逐渐升高，高频信息也得到了很好地改善，反演 结果非常接近于真实速度模型(图2)，作为对比，这里给出采用单尺度波形反演方法得到的 速度场(图13)，因为初始速度模型(图5)极不准确，所以很难得到较好的反演结果，层位 信息与速度信息都没有很好地反演出来。本发明的反演结果(图12)相比于采用统一曲网格 得到的波形反演速度场(图14)也更加准确。

为此本发明提供了一种基于辅助坐标系的起伏地表波形反演方法，能够很好地反演剧烈 起伏地表的速度场，为高精度成像方法提供准确的偏移速度场。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的 技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护 范围。