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

摘要

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

主权项

一种利用反射波信息反演速度场中低波数成分的方法，包括：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)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\∈\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{\τ}=\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{\∂p(x_r,t;x_s)}{{\∂m}_0}=\<\mathrm{\δG}(x_r,t;x),\frac{2}{v^3}\stackrel{\·\·}{G}(x,t;x_s)>+\<G(x_r,t;x),\frac{2}{v^3}\mathrm{\δ}\stackrel{\·\·}{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，直到满足收敛条件为止，并输出结果。