CN104391323A - 一种利用反射波信息反演速度场中低波数成分的方法 - Google Patents

Abstract

利用反射波信息反演速度场中低波数成分的方法，包括：1)建立初始速度模型，并进行深度偏移得到反射系数；2)利用反偏移算子从背景速度场中获取反射波信息；3)利用动态图像校正方法计算反射波旅行时差；4)计算伴随源并求取反射波路径，将伴随源沿反射波路径反传，得到多炮梯度；5)确定迭代步长并更新速度场；6)重复1)至5)步，通过判定迭代收敛条件最终得到反演结果。本发明能够利用平滑速度场进行反偏移得到模拟反射波，从而构建反射波路径来反演背景速度场，为全波形反演提供很好的初始速度模型，所求取的反射波旅行时差舍弃了人工拾取过程，避免了陷入局部极值，提高了反演的稳定性，可实现地震勘探中自动化速度建模过程。

Description

一种利用反射波信息反演速度场中低波数成分的方法

技术领域

本发明属于油气勘探地震资料处理领域，具体涉及一种能有效的利用反射波信息来反演速度场中低波数成分的方法。 

背景技术

目前，层析反演主要是在折射波(初至波)基础上做的，主要目的是为了解决近地表的静校正问题，难以做到中深层速度建模。而全波形反演(FWI)是一种基于波动方程的波形匹配反演方法，其波形信息包括走时信息及振幅信息。但是应用中存在一些问题：地震数据缺乏低频及大偏移距反射波；地震波振幅信息跟模型低波数成分呈较强的非线性；模拟地震数据中地震波振幅具有不确定性。总之，常规的波形反演利用折射波数据取得一些应用效果，恢复一些浅层信息，对中深层的反演效果较差。对于透射波和折射波，其具有较大的入射角孔径，可以用来恢复速度模型中的长波长分量；但是对于反射波，由于较小的反射角孔径，通常用于恢复速度场短波长分量。当初始模型与真实模型偏离较大时FWI通常会陷入局部极值，目标泛函很难收敛到全局极小值。 

发明内容

本发明的目的是提供一种利用反射波信息反演速度场中低波数成分的方法，以克服现有技术的不足。 

本发明的技术构思是： 

由于基于波动方程的旅行时反演兼顾了射线反演与波形反演的优点，将波动方程旅行时反演应用到反射波，能够实现中深层低频信息的反演。本发明利用偏移与反偏移算子从背景速度场中获取反射波信息，构建了反射波路径，进而建立波动方程反射层析的梯度公式。通过优化合成地震记录和实际数据之间反射波走时残差，实现了利用反射波信息来反演速度场中低波数成分的过程。相比振幅信息，旅行时反演提供了更加可靠的时间信息。因此波动方程反射旅行时反演，降低了其对波动方程的依赖程度，为提高中深层全波形反演精度奠定基础。 

一种利用反射波信息反演速度场中低波数成分的方法，包括： 

0.获得预处理后的地震观测反射波场p₀； 

1.建立初始深度域地震速度模型(速度场)，并利用该速度模型进行深度偏移，得到偏移剖面，在此选择逆时偏移方法进行成像，得到反射系数； 

2.根据步骤1得到的反射系数以及深度域地震速度模型，利用反偏移方法进行反偏移，得到模拟反射波场； 

所述的反偏移方法是通过入射波场与成像值的互相关产生所需要的反射波，反偏移过程如下： 

$(\frac{1}{v^2}\frac{{\∂}^2}{{\∂t}^2}-{\▿}^2)q=s---\left(1\right)$

$(\frac{1}{v^2}\frac{{\∂}^2}{{\∂t}^2}-{\▿}^2)q=I.q---\left(2\right)$

其中：p是反射波场，q是背景波场，v是平滑的背景速度场，s是震源函数，t是地震波传播旅行时间，I是叠加成像值，通常也称为反射系数； 

(2)式左边项中的反射波场p是由右边项作为震源产生的；(1)式与(2)式两个解耦方程得到反射波场p，从而降低了传统全波形反演的周波跳跃现象； 

其特征在于还包括以下步骤： 

3.利用动态图像校正方法计算模拟反射波场与观测反射波场的走时残差： 

两道地震记录f(i)与g(i)，分别对应模拟波场和观测反射波场的一道数据，两道地震合成记录之间存在时移函数s(i)，定义时差匹配面板： 

e[i,l]＝(f[i]-g[i+l])²  (3) 

其中：i是时间采样点序号，l是延迟变量,e是两道之间的匹配误差； 

将提取准确时移量u[0:N-1]的问题分以下两步来实现： 

通过对时差匹配面板(3)的双向平滑处理和对其反向追踪准确时移量两个步骤来实现两道记录f(i)与g(i)的最佳匹配； 

双向平滑如下： 

${\tilde{e}}_f[0,l]=e[0,l],$

${\tilde{e}}_f[i,l]=e[i,l]+\min \left\{\begin{array}{c}{\tilde{e}}_f[i-1,l-1]\\ {\tilde{e}}_f[i-1,l]\\ {\tilde{e}}_f[i-1,l+1]\end{array}\right.,i=\mathrm{1,2},...,N-1---\left(4\right)$

${\tilde{e}}_r[N-1,l]=e[N-1,l],$

${\tilde{e}}_r[i,l]=e[i,l]+\min \left\{\begin{array}{c}{\tilde{e}}_f[i+1,l-1]\\ {\tilde{e}}_r[i+1,l]\\ {\tilde{e}}_r[i+1,l+1]\end{array}\right.,i=N-2,N-3...,0---\left(5\right)$

$\tilde{e}[i.l]={\tilde{e}}_f[i.l]+{\tilde{e}}_r[i,l]-e[i,l].---\left(6\right)$

其中：i是时间采样点序号，l是延迟变量，N是时间采样点数，是正向平滑后的时差匹配面板，是反向平滑后的时差匹配面板，是双向平滑后的时差匹配面板； 

用双向平滑处理，从后往前追踪准确时移量，首先计算第N个准确时移量u[N-1]，逐次向前计算。 

$u[N-1]=\underset{l}{\arg \min }{\tilde{e}}_f[N-1,l],$

$u[i-1]=\underset{l\∈\left\{u\right[i]-1,u[i],u[i]+1\}}{\arg \min }{\tilde{e}}_f[i-1,l],i=N-1,N-2,...,1---\left(7\right)$

采用动态图像校正方法，即通过优化下列非线性问题来求取模拟反射波场p与观测反射波场p0的多道准确时移量τ： 

$\mathrm{\τ}=\underset{l}{\arg \min }D\left(l\right)---\left(8\right)$

其中 $D\left(l\right)=\frac{1}{2}\∫{\mathrm{dx}}_s{\mathrm{dx}}_r\mathrm{dt}{(p(x_r,t;x_s)-p_0(x_r,t+l(x_r,t;x_s);x_s))}^2---\left(9\right)$

满足线性约束条件： $\left|\frac{\∂\mathrm{\τ}}{\∂t}\right|{\≤\mathrm{\σ}}_t,\left|\frac{\∂\mathrm{\τ}}{\∂x_r}\right|{\≤\mathrm{\σ}}_r,\left|\frac{\∂\mathrm{\τ}}{\∂x_s}\right|{\≤\mathrm{\σ}}_s---\left(10\right)$

其中：l是延迟变量，τ是两个反射波场间存在的多道准确时移量，t,x_r,x_s分别代表时间方向、检波点方向、激发点方向； 

约束条件(10)式控制着估计的时差在时间t、接收点x_r、激发点x_s三个方向上的变化率，进而保证了其平滑度； 

4.求取反射波路径： 

对于散射波，建立介质散射模型，该模型分为光滑的背景模型与模型扰动的和 

m＝m₀+δm  (11) 

其中：m是地球介质参数，这里指速度，m₀代表背景速度模型，δm代表速度扰动； 

由此可以导出基于一阶Born近似的扰动场 

$\mathrm{\δG}(x_r,\mathrm{\ω};x_s)={\mathrm{\ω}}^2\underset{V}{\∫}G(x^{\′},\mathrm{\ω};x_s)\mathrm{\δm}\left(x^{\′}\right)G(x_r,\mathrm{\ω};x^{\′}){\mathrm{dx}}^{\′}---\left(12\right)$

其中，x_s,x_r分别为炮点和检波点位置，δG是基于一阶Born近似的扰动场，ω是圆频率，V代表地球介质体，δm是速度扰动，G代表格林函数； 

利用下列扰动波场δG对背景模型参数m₀的Fréchet导数计算出反射波路径 

$\frac{\∂\mathrm{\δG}}{\∂m_0\left(x^{\′}\right)}{|}_{\mathrm{\δm}}={\mathrm{\ω}}^2(\mathrm{\δG}(x_r,\mathrm{\ω};x^{\′})G(x^{\′},\mathrm{\ω};x_s)+\mathrm{\δG}(x^{\′},\mathrm{\ω};x_s)G(x_r,\mathrm{\ω};x^{\′}))---\left(13\right)$

其中，x_s,x_r分别为炮点和检波点位置，ω是圆频率，δm是速度扰动，G代表格林函 数，δG(x′,ω；x_s)，δG(x_r,ω；x′)分别是一阶Born近似下，由模型扰动引起的反射波场； 

5.根据步骤3求取的多道准确时移量及其步骤4得到的反射波路径计算伴随源；基于多道准确时移量构建如下目标函数 

$E_T=\frac{1}{2}\underset{s}{\mathrm{\Σ}}\underset{r}{\mathrm{\Σ}}{\left|\right|\mathrm{\τ}(x_r,t;x_s)\left|\right|}^2---\left(14\right)$

其中：x_s,x_r分别为炮点和检波点位置，多道准确时移量τ是两个反射波场间存在的时差，E_T是走时残差目标函数； 

该目标泛函对背景模型参数m₀的梯度 

$\frac{\∂E_T}{\∂m_0}=\underset{s}{\mathrm{\Σ}}\underset{r}{\mathrm{\Σ}}\left(\frac{{\∂\mathrm{\τ}}^T}{\∂m_0}\right)\mathrm{\τ}---\left(15\right)$

其中：m₀是背景模型参数，这里指速度；s,r代表炮点和检波点；T代表矩阵的转置； 

根据波动方程旅行时反演中求取伴随源的方法，可以得到梯度公式 

$\frac{\∂E_T}{\∂m_0}=\underset{s}{\mathrm{\Σ}}\underset{r}{\mathrm{\Σ}\∫}\mathrm{dt}\left(\frac{\∂p(x_r,t;x_s)}{\∂m_0}\right)\·f^{\mathrm{adj}}---\left(16\right)$

其中，f^adj称为伴随源，具有如下形式 

$f^{\mathrm{adj}}=\frac{{\stackrel{\·}{p}}_0(x_r,t+\mathrm{\τ};x_s)\mathrm{\τ}(x_r,t+\mathrm{\τ};x_s)}{\∂F/\∂\mathrm{\σ}}---\left(17\right)$

其中：p₀是观测波场，τ是两个反射波场间存在的时差； 

(16)式中的即为(13)式所求取的反射波路径 

$\frac{\&PartialD;p(x_r,t;x_s)}{\&PartialD;m_0}=<\mathrm{\&delta;G}(x_r,t;x),\frac{2}{v^3}\stackrel{\&CenterDot;\&CenterDot;}{G}(x,t;x_s)>+<G(x_r,t;x),\frac{2}{v^3}\mathrm{\&delta;}\stackrel{\&CenterDot;\&CenterDot;}{G}(x,t;x_s)>---\left(18\right)$

其中：x_s,x_r分别为炮点和检波点位置，代表格林函数对时间的二阶导数，δG(x,t；x_s)，δG(x_r,t；x)分别是一阶Born近似下，由模型扰动引起的反射波场； 

6.联合(16)、(17)、(18)两式便得到了波动方程反射旅行时层析的梯度公式； 

7.计算迭代步长 

即按照常规方法计算步骤6得到的梯度公式的步长； 

8.根据梯度方向和步长计算速度模型中低波数成分的更新量，更新步骤1建立的速度模型并检查迭代收敛条件，如果满足收敛条件则停止迭代并输出速度模型，否则继续上述步骤1-7，直到满足收敛条件为止，并输出结果。 

上述步骤3中，选择 ${\mathrm{\&sigma;}}_t=25\%\frac{\mathrm{\&Delta;t}}{\mathrm{\&Delta;t}},{\mathrm{\&sigma;}}_r=100\%\frac{\mathrm{\&Delta;t}}{\mathrm{\&Delta;}{x}_r},{\mathrm{\&sigma;}}_s=\&infin;\frac{\mathrm{\&Delta;t}}{\mathrm{\&Delta;}{x}_s};$ Δt，Δx_r，Δx_s分别为三个方向的采样间隔。 

上述步骤7的步长为固定的小步长，所述的固定的小步长是以上述背景速度场v的数值除以步骤6得到的梯度后乘以1％。 

上述步骤8所述的收敛条件是：相邻的前后两次迭代的更新量小于1‰即可认为收敛。 

发明的效果 

波动方程反射旅行时反演方法在反演方面，有着其他方法不具备的优势，其具体优势和特点表现在以下几个方面： 

一、该方法超越了以前地震勘探中仅仅利用折射波(或初至波)进行层析静校正的局限，能够有效利用反射波来实现层析反演。一定程度上解决了地震勘探中中深层速度建模问题，其反演的深度与地震资料的观测孔径有关，孔径越大，其反演的深度越大。 

二、该方法能够利用平滑速度场进行反偏移得到模拟反射波，从而可以构建反射波路径来反演背景速度场。普通的全波形反演对初始速度模型的要求较高，该方法反演的中深层速度模型能够为全波形反演能够提供一个很好的初始速度模型。 

三、动态图像校正方法求取反射波旅行时差舍弃了人工拾取过程以及互相关方法，以免其陷入局部极值。该方法可以较为准确地估计复杂介质模型的反射旅行时差，而且时差纵横向变化比较平缓，能够提高反演的稳定性。 

四、该方法作为一种完整的反演体系，可以实现地震勘探中自动化速度建模过程。较人工交互拾取速度分析能够体现其价值所在，提高了速度建模的效率。 

附图说明

图1是本发明的流程图。 

图2为两道合成地震记录(a)f(i)；(b)g(i)。 

图3为误差匹配面板：(a)未平滑；(b)正向平滑；(c)反向平滑；(d)双向平滑。 

图4为提取的时差曲线。 

图5为本发明反射波路径示意图。 

图6为本发明平层模型试算示意图：(a)真实速度模型；(b)初始速度模型；(c)最终迭代的速度模型。 

图7为本发明平层模型试算示意图：(a)第51炮真实模拟炮记录；(b)迭代第一次反偏移的炮记录；(c)迭代最后一次反偏移的炮记录；(d)第一次模拟记录与真实记录的多道准确时移量(单位秒)；(e)第一次的伴随源；(f)最后一次的多道准确时移量(单位秒)。 

图8为本发明平层模型试算中反演速度场不同水平位置处的速度曲线对比图：(a)Distance＝2.0km处；(b)Distance＝4.0km处；(c)Distance＝5.5km处。 

图9为本发明sigsbee2b模型试算示意图：(a)真实速度模型；(b)初始速度模型；(c)以(b)为初始速度模型FWI的结果；(d)波动方程反射层析的结果(WERTI)；(e)以(d)为初始速度模型FWI的结果(WERTI+FWI)。 

图10为本发明sigsbee2b模型试算示意图：(a)、(b)分别为由图9(c)(e)所模拟的炮记录与真实记录的多道准确时移量(单位毫秒)；(c)图9(a)(c)(e)中Distance＝2.9km处纵向速度对比曲线。 

具体实施方式

如图1所示，一种利用反射波信息反演速度场中低波数成分的方法，包括： 

0.获得预处理后的地震观测反射波场p₀； 

1.建立初始深度域地震速度模型(速度场)，并利用该速度模型进行深度偏移，得到偏移剖面，在此选择逆时偏移方法进行成像，得到反射系数； 

2.根据步骤1得到的反射系数以及深度域地震速度模型，利用反偏移方法进行反偏移，得到模拟反射波场； 

所述的反偏移方法是通过入射波场与成像值的互相关产生所需要的反射波，反偏移过程如下： 

$(\frac{1}{v^2}\frac{{\&PartialD;}^2}{{\&PartialD;t}^2}-{\&dtri;}^2)q=s---\left(1\right)$

$(\frac{1}{v^2}\frac{{\&PartialD;}^2}{{\&PartialD;t}^2}-{\&dtri;}^2)q=I.q---\left(2\right)$

其中：p是反射波场，q是背景波场，v是平滑的背景速度场，s是震源函数，t是地震波传播旅行时间，I是叠加成像值，通常也称为反射系数； 

(2)式左边项中的反射波场p是由右边项作为震源产生的；(1)式与(2)式两个解耦方程得到反射波场p，从而降低了传统全波形反演的周波跳跃现象； 

其特征在于还包括以下步骤： 

3.利用动态图像校正方法计算模拟反射波场与观测反射波场的走时残差： 

动态时间校正(Dynamic Time Warping)是一种经典的用来估计两个时间序列时差的方法，该方法解决了在线性约束条件下的非线性最优化问题，其估计的时差是一种全局最优解。 

如图2所示，假设存在两道地震合成记录f(i)与g(i)(分别对应模拟波场和观测反射波场的一道数据)，其采样点N＝501，采样间隔0.01s。两道地震合成记录之间存在时移函数 s(i)，定义时差匹配面板： 

e[i,l]＝(f[i]-g[i+l])²  (3) 

其中：i是时间采样点序号，l是延迟变量,e是两道之间的匹配误差； 

图3(a)显示了上述两道合成记录的时差匹配面板，当延迟变量l近似等于时移函数s时，匹配误差达到了最小值，接近于零；要提取的准确时移量u[0:N-1]便是使时差匹配面板e达到最小值时所对应的延迟变量l，但是，直接在该时差匹配面板上求取准确时移量u[0:N-1]存在一定困难；故将该最优化问题分两步来实现： 

通过对时差匹配面板(3)的平滑处理和对其反向追踪准确时移量两个步骤来实现两道记录f(i)与g(i)的最佳匹配； 

平滑处理有正向平滑、反向平滑以及双向平滑。 

正向平滑： 

${\tilde{e}}_f[0,l]=e[0,l],$

${\tilde{e}}_f[i,l]=e[i,l]+\min \left\{\begin{array}{c}{\tilde{e}}_f[i-1,l-1]\\ {\tilde{e}}_f[i-1,l]\\ {\tilde{e}}_f[i-1,l+1]\end{array}\right.,i=\mathrm{1,2},...,N-1---\left(4\right)$

反向平滑： 

${\tilde{e}}_r[N-1,l]=e[N-1,l],$

${\tilde{e}}_r[i,l]=e[i,l]+\min \left\{\begin{array}{c}{\tilde{e}}_f[i+1,l-1]\\ {\tilde{e}}_r[i+1,l]\\ {\tilde{e}}_r[i+1,l+1]\end{array}\right.,i=N-2,N-3...,0---\left(5\right)$

双向平滑： 

$\tilde{e}[i.l]={\tilde{e}}_f[i.l]+{\tilde{e}}_r[i,l]-e[i,l].---\left(6\right)$

其中：i是时间采样点序号，l是延迟变量，N是时间采样点数，是正向平滑后的时差匹配面板，是反向平滑后的时差匹配面板，是双向平滑后的时差匹配面板； 

图3(b)、(c)、(d)分别展示了对图3(a)进行正向平滑、反向平滑和双向平滑处理后的效果。正向平滑处理后误差面板从前往后平滑程度依次增大，反向平滑处理效果反之。而经过双向平滑处理后，误差匹配面板整体得到了很好地平滑。这为准确求取时差奠定了很好的铺垫，同时也提高了其抗噪能力。 

在进行反向追踪准确时移量u[0:N-1]时，追踪的方向与平滑的方向是相反的， 

如果正向平滑处理后，就应该从后往前追踪准确时移量。首先计算第N个准确时移量u[N-1]，逐次向前计算。 

$u[N-1]=\underset{l}{\arg \min }{\tilde{e}}_f[N-1,l],$

$u[i-1]=\underset{l\&Element;\left\{u\right[i]-1,u[i],u[i]+1\}}{\arg \min }{\tilde{e}}_f[i-1,l],i=N-1,N-2,...,1---\left(7\right)$

如果是反向平滑后从前往后进行追踪。 

一般用双向平滑处理，可以任选一个方向进行追踪准确时移量。该方法的精度较高，能够很好地刻画记录的时移量。在此我们采用双向平滑处理,即(6)式。图4显示了用上述方法计算得到的时差曲线，可以看出，该方法的精度较高，能够很好地刻画其准确时移量。 

在估计两道记录f(i)与g(i)间的准确时移量时，当时移量随着时间或空间变化剧烈，或者图像含有噪音时，动态时间校正方法存在一定困难。由此需将动态时间校正法进行扩展，采用动态图像校正方法，动态图像校正方法通过优化下列非线性问题来求取模拟反射波场p与观测反射波场p₀的多道准确时移量τ： 

$\mathrm{\&tau;}=\underset{l}{\arg \min }D\left(l\right)---\left(8\right)$

其中 $D\left(l\right)=\frac{1}{2}\&Integral;{\mathrm{dx}}_s{\mathrm{dx}}_r\mathrm{dt}{(p(x_r,t;x_s)-p_0(x_r,t+l(x_r,t;x_s);x_s))}^2---\left(9\right)$

满足线性约束条件： $\left|\frac{\&PartialD;\mathrm{\&tau;}}{\&PartialD;t}\right|{\&le;\mathrm{\&sigma;}}_t,\left|\frac{\&PartialD;\mathrm{\&tau;}}{\&PartialD;x_r}\right|{\&le;\mathrm{\&sigma;}}_r,\left|\frac{\&PartialD;\mathrm{\&tau;}}{\&PartialD;x_s}\right|{\&le;\mathrm{\&sigma;}}_s---\left(10\right)$

其中：l是延迟变量，τ是两个反射波场间存在的多道准确时移量，t,x_r,x_s分别代表时间方向、检波点方向、激发点方向； 

约束条件(10)式控制着估计的时差在时间t、接收点x_r、激发点x_s三个方向上的变化率，进而保证了其平滑度。本文中，选择 ${\mathrm{\&sigma;}}_t=25\%\frac{\mathrm{\&Delta;t}}{\mathrm{\&Delta;t}},{\mathrm{\&sigma;}}_r=100\%\frac{\mathrm{\&Delta;t}}{\mathrm{\&Delta;}{x}_r},{\mathrm{\&sigma;}}_s=\&infin;\frac{\mathrm{\&Delta;t}}{\mathrm{\&Delta;}{x}_s};$ Δt，Δx_r，Δx_s分别为三个方向的采样间隔；σ_s为无穷大意味着不同位置的反射波场的时差可以是相互独立的； 

4求取反射波路径： 

散射是由于地球介质的不均匀性引起的波路径变化的现象，当地震波到达这些不均匀构造后，相互作用产生的波称为散射波。基于Born近似可以将散射波分为一次散射波、二次散射波等。其产生机理同样可以用惠更斯原理予以解释。 

对于散射波，建立介质散射模型，在该模型中，可以理解为将地球介质分解为光滑的背景模型与模型扰动的和 

m＝m₀+δm  (11) 

其中：m是地球介质参数，这里指速度，m₀代表背景速度模型，δm代表速度扰动； 

由此可以导出基于一阶Born近似的扰动场 

$\mathrm{\&delta;G}(x_r,\mathrm{\&omega;};x_s)={\mathrm{\&omega;}}^2\underset{V}{\&Integral;}G(x^{\&prime;},\mathrm{\&omega;};x_s)\mathrm{\&delta;m}\left(x^{\&prime;}\right)G(x_r,\mathrm{\&omega;};x^{\&prime;}){\mathrm{dx}}^{\&prime;}---\left(12\right)$

其中，x_s,x_r分别为炮点和检波点位置，δG是基于一阶Born近似的扰动场，ω是圆频率，V代表地球介质体，δm是速度扰动，G代表格林函数。 

Xu et al.(2012)推导了扰动波场δG对背景模型参数m₀的Fréchet导数 

$\frac{\&PartialD;\mathrm{\&delta;G}}{\&PartialD;m_0\left(x^{\&prime;}\right)}{|}_{\mathrm{\&delta;m}}={\mathrm{\&omega;}}^2(\mathrm{\&delta;G}(x_r,\mathrm{\&omega;};x^{\&prime;})G(x^{\&prime;},\mathrm{\&omega;};x_s)+\mathrm{\&delta;G}(x^{\&prime;},\mathrm{\&omega;};x_s)G(x_r,\mathrm{\&omega;};x^{\&prime;}))---\left(13\right)$

其中，x_s,x_r分别为炮点和检波点位置，ω是圆频率，δm是速度扰动，G代表格林函数。δG(x′,ω；x_s)，δG(x_r,ω；x′)分别是一阶Born近似下，由模型扰动引起的反射波场； 

可以看出，该梯度计算不仅与背景模型有关，而且还与模型扰动有关；而反射波场的获取可以通过波动方程反偏移来计算，由(13)式计算出的结果便是反射波路径，如图5所示；可以看出，该式进行的是反射波层析成像。 

5根据步骤3求取的多道准确时移量及其步骤4得到的反射波路径计算伴随源；基于反射波走时残差(旅行时差)构建如下目标函数 

$E_T=\frac{1}{2}\underset{s}{\mathrm{\&Sigma;}}\underset{r}{\mathrm{\&Sigma;}}{\left|\right|\mathrm{\&tau;}(x_r,t;x_s)\left|\right|}^2---\left(14\right)$

其中：x_s,x_r分别为炮点和检波点位置，多道准确时移量τ是两个反射波场间存在的时差，E_T是走时残差目标函数。 

该目标泛函对背景模型参数m₀的梯度 

$\frac{\&PartialD;E_T}{\&PartialD;m_0}=\underset{s}{\mathrm{\&Sigma;}}\underset{r}{\mathrm{\&Sigma;}}\left(\frac{{\&PartialD;\mathrm{\&tau;}}^T}{\&PartialD;m_0}\right)\mathrm{\&tau;}---\left(15\right)$

其中：m₀是背景模型参数，这里指速度；s,r代表炮点和检波点；T代表矩阵的转置。 

根据Luo Y(1991)波动方程旅行时反演中求取伴随源的方法，可以得到梯度公式 

$\frac{\&PartialD;E_T}{\&PartialD;m_0}=\underset{s}{\mathrm{\&Sigma;}}\underset{r}{\mathrm{\&Sigma;}\&Integral;}\mathrm{dt}\left(\frac{\&PartialD;p(x_r,t;x_s)}{\&PartialD;m_0}\right)\&CenterDot;f^{\mathrm{adj}}---\left(16\right)$

其中，f^adj称为伴随源，具有如下形式 

$f^{\mathrm{adj}}=\frac{{\stackrel{\&CenterDot;}{p}}_0(x_r,t+\mathrm{\&tau;};x_s)\mathrm{\&tau;}(x_r,t+\mathrm{\&tau;};x_s)}{\&PartialD;F/\&PartialD;\mathrm{\&sigma;}}---\left(17\right)$

其中：p₀是观测波场，τ是两个反射波场间存在的时差。 

(16)式中的即为(13)式所求取的反射波路径 

$\frac{\&PartialD;p(x_r,t;x_s)}{\&PartialD;m_0}=<\mathrm{\&delta;G}(x_r,t;x),\frac{2}{v^3}\stackrel{\&CenterDot;\&CenterDot;}{G}(x,t;x_s)>+<G(x_r,t;x),\frac{2}{v^3}\mathrm{\&delta;}\stackrel{\&CenterDot;\&CenterDot;}{G}(x,t;x_s)>---\left(18\right)$

其中：x_s,x_r分别为炮点和检波点位置，代表格林函数对时间的二阶导数，δG(x,t；x_s)，δG(x_r,t；x)分别是一阶Born近似下，由模型扰动引起的反射波场； 

6.联合(16)、(17)、(18)两式便得到了波动方程反射旅行时层析的梯度公式； 

7.计算迭代步长 

即按照常规方法计算步骤6得到的梯度公式的步长； 

(为了减小计算量，可以是固定的小步长，如以上述背景速度场v的数值除以步骤6得到的梯度后乘以1％) 

8.根据梯度方向和步长计算速度模型中低波数成分的更新量，更新步骤1建立的速度模型并检查迭代收敛条件，如果满足收敛条件则停止迭代并输出速度模型，否则继续上述步骤1-7，直到满足收敛条件为止，并输出结果。所述的收敛条件如：相邻的前后两次迭代的更新量小于1‰即可认为收敛，即可输出结果。 

实施实例 

1)平层模型测试。如图6(a)所示，该速度模型纵横向采样点分别为300、800，纵横向采样间隔均为10m。速度从上到下依次为2000m/s、2100m/s、2300m/s、2500m/s。观测系统为地面放炮地面接收，放炮炮数为101炮，炮间距50m，第一炮激发点在1.5km处，每炮300道接收，道间距10m。震源子波为主频25Hz的Ricker子波，炮记录最大接收时间为2.5s，采样率1ms。 

图6(b)为反演的初始速度模型——2000m/s的均匀速度场。图6(c)是迭代20次后得到的反演速度模型，可以看出，横向上2km～6km间的速度场得到了很好的恢复，两侧受照明不足影响，反演结果不准确。我们选取中间炮第51炮进行分析，图7(a)为真实速度模型的炮记录，图7(b)为用初始速度模型反偏移得到的炮记录，图7(d)为(a)、(b)的多道准确时移量，时移量单位为秒。可以看出，第一个界面的多道准确时移量基本为零，第二、三个反射界面存在很大时移量，图7(e)为计算得到的初次迭代伴随源。图7(c)为用最终反演的速度模型反偏移得到的炮记录，可以看出，反偏移后的噪音也有所减少，我们再对其求取多道准确时移量，如图7(f)所示，时移量单位同样为秒。显然，相比初始模型，其时移量很大程度上减小了。图8(a)、(b)、(c)分别是横向位置为2.0km、4.0km和5.5km处初始速度、真实速度与20次反演迭代速度曲线对比，可以看出速度场界面信息以及低波数成分基本都得到了很好的恢复，从而验证了其精度。 

2)sigsbee2b模型试算。如图9(a)所示，我们截取sigsbee2b模型的一部分进行测试。速度模型纵横向采样点分别为300、700，纵横向采样间隔分别为8m、10m。观测系统 为地面放炮地面接收，放炮炮数为101炮，炮间距40m，第一炮激发点在1.5km处，每炮300道接收，道间距10m。震源子波为主频15Hz的Ricker子波，炮记录最大接收时间为2.0s，采样率1ms。 

图9(b)为常梯度速度模型，在其基础上直接进行FWI，得到反演结果如图9(c)所示，可以看到浅层反演效果较好，中深层由于背景速度不准确，FWI反演结果陷入局部极小值。我们将(b)图作为初始速度模型首先进行波动方程反射旅行时层析，图9(d)是其层析的结果，可以看到速度场低波数成分基本得到了恢复。我们将(d)作为初始速度模型进行FWI，反演结果如图9(e)所示。与图9(c)图比较，从圈内可以看到采用本文方法更新的背景速度，FWI中深层反演结果得到有效改善。图10(a)、(b)分别显示了用图9(c)、(e)模拟的炮记录与真实记录的多道准确时移量，时移量单位为毫秒。可以看出，在浅层，多道准确时移量很大程度上减小了。由此也论证了普通的全波形反演对速度场高波数成分比较敏感，而引入反射波路径后，利用波动方程反射旅行时层析可以较好地恢复速度场中低波数成分，为后续的全波形反演能够提供一个较好的初始速度场。图10(c)显示了图9(a)、(c)、(e)中Distance＝2.9km处纵向速度对比曲线，可以看到WERTI+FWI的纵向速度曲线更加接近于真实速度曲线，具有更高的反演精度。 

Claims (4)

1.一种利用反射波信息反演速度场中低波数成分的方法，包括：

0.获得预处理后的地震观测反射波场p₀；

1.建立初始深度域地震速度模型(速度场)，并利用该速度模型进行深度偏移，得到偏移剖面，在此选择逆时偏移方法进行成像，得到反射系数；

2.根据步骤1得到的反射系数以及深度域地震速度模型，利用反偏移方法进行反偏移，得到模拟反射波场；

所述的反偏移方法是通过入射波场与成像值的互相关产生所需要的反射波，反偏移过程如下：

$(\frac{1}{v^2}\frac{{\&PartialD;}^2}{{\&PartialD;t}^2}-{\&dtri;}^2)q=s---\left(1\right)$

$(\frac{1}{v^2}\frac{{\&PartialD;}^2}{{\&PartialD;t}^2}-{\&dtri;}^2)p=I.q---\left(2\right)$

其中：p是反射波场，q是背景波场，v是平滑的背景速度场，s是震源函数，t是地震波传播旅行时间，I是叠加成像值，通常也称为反射系数；

(2)式左边项中的反射波场p是由右边项作为震源产生的；(1)式与(2)式两个解耦方程得到反射波场p，从而降低了传统全波形反演的周波跳跃现象；

其特征在于还包括以下步骤：

3.利用动态图像校正方法计算模拟反射波场与观测反射波场的走时残差：

两道地震记录f(i)与g(i)，分别对应模拟波场和观测反射波场的一道数据，两道地震合成记录之间存在时移函数s(i)，定义时差匹配面板：

e[i,l]＝(f[i]-g[i+l])²  (3)

其中：i是时间采样点序号，l是延迟变量,e是两道之间的匹配误差；

将提取准确时移量u[0:N-1]的问题分以下两步来实现：

通过对时差匹配面板(3)的双向平滑处理和对其反向追踪准确时移量两个步骤来实现两道记录f(i)与g(i)的最佳匹配；

双向平滑如下：

$\begin{array}{c}{\tilde{e}}_f[0,l]=e[0,l],\\ {\tilde{e}}_f[i,l]e[i,l]+\min \left\{\begin{array}{c}{\tilde{e}}_f[i-1,l-1]\\ {\tilde{e}}_f[i-1,l]\\ {\tilde{e}}_f[i-1,l+1]\end{array},i=\mathrm{1,2},...,N-1\right.\end{array}---\left(4\right)$

$\begin{array}{c}{\tilde{e}}_r[N-1,l]=e[N-1,l],\\ {\tilde{e}}_r[i,l]=e[i,l]+\min \begin{array}{}\end{array}\left\{\begin{array}{c}{\tilde{e}}_r[i+1,l-1]\\ {\tilde{e}}_r[i+1,l]\\ {\tilde{e}}_r[i+1\mathrm{ml}+1]\end{array}\right.,i=N-2,N-3,...,0\end{array}---\left(5\right)$

$\tilde{e}[i,l]={\tilde{e}}_f[i,l]+{\tilde{e}}_r[i,l]-e[i,l].---\left(6\right)$

其中：i是时间采样点序号，l是延迟变量，N是时间采样点数，是正向平滑后的时差匹配面板，是反向平滑后的时差匹配面板，是双向平滑后的时差匹配面板；

用双向平滑处理，从后往前追踪准确时移量，首先计算第N个准确时移量u[N-1]，逐次向前计算。

$\begin{array}{c}u[N-1]=\underset{l}{\arg \min }{\tilde{e}}_f[N-1,l],\\ u[i-1]=\underset{l\&Element;\left\{u\right[i]-1,u[i],u[i]+1\}}{\arg \min }{\tilde{e}}_f[i-1,l],i=N-1,N-2,...,1\end{array}---\left(7\right)$

采用动态图像校正方法，即通过优化下列非线性问题来求取模拟反射波场p与观测反射波场p₀的多道准确时移量τ：

$\mathrm{\&tau;}=\underset{l}{\arg \min }D\left(l\right)---\left(8\right)$

其中 $D\left(l\right)=\frac{1}{2}{\&Integral;\mathrm{dx}}_s{\mathrm{dx}}_r\mathrm{dt}{(p(x_r,t;x_s)-p_0(x_r,t+l(x_r,t;x_s);x_s))}^2---\left(9\right)$

满足线性约束条件： $\left|\frac{\&PartialD;\mathrm{\&tau;}}{\&PartialD;t}\right|\&le;{\mathrm{\&sigma;}}_t,\left|\frac{\&PartialD;\mathrm{\&tau;}}{{\&PartialD;x}_r}\right|\&le;{\mathrm{\&sigma;}}_r,\left|\frac{\&PartialD;\mathrm{\&tau;}}{{\&PartialD;x}_s}\right|\&le;{\mathrm{\&sigma;}}_s---\left(10\right)$

其中：l是延迟变量，τ是两个反射波场间存在的多道准确时移量，t,x_r,x_s分别代表时间方向、检波点方向、激发点方向；

约束条件(10)式控制着估计的时差在时间t、接收点x_r、激发点x_s三个方向上的变化率，进而保证了其平滑度；

4.求取反射波路径：

对于散射波，建立介质散射模型，该模型分为光滑的背景模型与模型扰动的和

m＝m₀+δm  (11)

其中：m是地球介质参数，这里指速度，m₀代表背景速度模型，δm代表速度扰动；

由此可以导出基于一阶Born近似的扰动场

$\mathrm{\&delta;G}(x_r;\mathrm{\&omega;};x_s)={\mathrm{\&omega;}}^2\underset{V}{\&Integral;}G(x^{\&prime;},\mathrm{\&omega;};x_s)\mathrm{\&delta;m}\left(x^{\&prime;}\right)G(x_r,\mathrm{\&omega;};x^{\&prime;}){\mathrm{dx}}^{\&prime;}---\left(12\right)$

其中，x_s,x_r分别为炮点和检波点位置，δG是基于一阶Born近似的扰动场，ω是圆频率，V代表地球介质体，δm是速度扰动，G代表格林函数；

利用下列扰动波场δG对背景模型参数m₀的Fréchet导数计算出反射波路径

$\frac{\&PartialD;\mathrm{\&delta;G}}{{\&PartialD;m}_0\left(x^{\&prime;}\right)}{|}_{\mathrm{\&delta;m}}={\mathrm{\&omega;}}^2(\mathrm{\&delta;G}(x_r,\mathrm{\&omega;};x^{\&prime;})G(x^{\&prime;},\mathrm{\&omega;};x_s)+\mathrm{\&delta;G}(x^{\&prime;},\mathrm{\&omega;};x_s)G(x_r,\mathrm{\&omega;};x^{\&prime;}))---\left(13\right)$

其中，x_s,x_r分别为炮点和检波点位置，ω是圆频率，δm是速度扰动，G代表格林函数，δG(x′,ω；x_s)，δG(x_r,ω；x′)分别是一阶Born近似下，由模型扰动引起的反射波场；

5.根据步骤3求取的多道准确时移量及其步骤4得到的反射波路径计算伴随源；基于多道准确时移量构建如下目标函数

$E_T=\frac{1}{2}\underset{s}{\mathrm{\&Sigma;}}\underset{r}{\mathrm{\&Sigma;}}{\left|\right|\mathrm{\&tau;}(x_r,t;x_s)\left|\right|}^2---\left(14\right)$

其中：x_s,x_r分别为炮点和检波点位置，多道准确时移量τ是两个反射波场间存在的时差，E_T是走时残差目标函数；

该目标泛函对背景模型参数m₀的梯度

$\frac{{\&PartialD;E}_T}{{\&PartialD;m}_0}=\underset{s}{\mathrm{\&Sigma;}}\underset{r}{\mathrm{\&Sigma;}}\left(\frac{{\&PartialD;\mathrm{\&tau;}}^T}{{\&PartialD;m}_0}\right)\mathrm{\&tau;}---\left(15\right)$

其中：m₀是背景模型参数，这里指速度；s,r代表炮点和检波点；T代表矩阵的转置；

根据波动方程旅行时反演中求取伴随源的方法，可以得到梯度公式

$\frac{{\&PartialD;E}_T}{{\&PartialD;m}_0}=\underset{s}{\mathrm{\&Sigma;}}\underset{r}{\mathrm{\&Sigma;}}\&Integral;\mathrm{dt}\left(\frac{\&PartialD;p(x_r,t;x_s)}{{\&PartialD;m}_0}\right)\&CenterDot;f^{\mathrm{adj}}---\left(16\right)$

其中，f^adj称为伴随源，具有如下形式

$f^{\mathrm{adj}}=\frac{{\stackrel{\&CenterDot;}{p}}_0(x_r,t+\mathrm{\&tau;};x_s)\mathrm{\&tau;}(x_r,t+\mathrm{\&tau;};x_s)}{\&PartialD;F/\&PartialD;\mathrm{\&tau;}}---\left(17\right)$

其中：p₀是观测波场，τ是两个反射波场间存在的时差；

(16)式中的即为(13)式所求取的反射波路径

$\frac{\&PartialD;p(x_r,t;x_s)}{{\&PartialD;m}_0}=<\mathrm{\&delta;G}(x_r,t;x),\frac{2}{v^3}\stackrel{\&CenterDot;\&CenterDot;}{G}(x,t;x_s)>+<G(x_r,t;x),\frac{2}{v^3}\mathrm{\&delta;}\stackrel{\&CenterDot;\&CenterDot;}{G}(x,t;x_s)>---\left(18\right)$

其中：x_s,x_r分别为炮点和检波点位置，代表格林函数对时间的二阶导数，δG(x,t；x_s)，δG(x_r,t；x)分别是一阶Born近似下，由模型扰动引起的反射波场；

6.联合(16)、(17)、(18)两式便得到了波动方程反射旅行时层析的梯度公式；

7.计算迭代步长

即按照常规方法计算步骤6得到的梯度公式的步长；

8.根据梯度方向和步长计算速度模型中低波数成分的更新量，更新步骤1建立的速度模型并检查迭代收敛条件，如果满足收敛条件则停止迭代并输出速度模型，否则继续上述步骤1-7，直到满足收敛条件为止，并输出结果。

2.如权利要求1所述的利用反射波信息反演速度场中低波数成分的方法，其特征在于上述步骤3中，选择 Δt，Δx_r，Δx_s分别为三个方向的采样间隔。

3.如权利要求1所述的利用反射波信息反演速度场中低波数成分的方法，其特征在于上述步骤7的步长为固定的小步长，所述的固定的小步长是以上述背景速度场v的数值除以步骤6得到的梯度后乘以1％。

4.如权利要求1所述的利用反射波信息反演速度场中低波数成分的方法，其特征在于上述步骤8所述的收敛条件是：相邻的前后两次迭代的更新量小于1‰即可认为收敛。