[发明专利]一种提高高陡构造速度反演效率的全波形联合反演方法

权利要求书

1.一种提高高陡构造速度反演效率的全波形联合反演方法，其特征在于：按照如下步骤进行：

步骤1：输入初始速度场的速度、野外观测炮记录及震源子波，并建立观测系统；

步骤2：使用下式表示的成像公式获得成像结果：

$I=\underset{x_s}{\Σ}\underset{t}{\Σ}uw/u^2$

其中I表示成像结果，u和w分别表示正传波场和反传波场，w＝L^*(R^*(p_obs))，L^*为反传算子，R^*表示将数据残差空间扩展到整个空间，p_obs为野外观测炮记录；

步骤3：将步骤2中的成像结果作为反射系数；

步骤4：应用线性正演模拟和线性波场反传，分别计算两种棱柱波prism1和prism2的梯度方向；

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)$

其中g(v)_prism1和g(v)_prism2分别表示prism1的梯度方向和prism2的梯度方向，v为速度场的速度，w'表示野外观测炮记录与数值模拟的地震记录残差的反传波场，x_s为震源位置，t为时间；

步骤5：将步骤4中所得的两个梯度方向相加即可求得棱柱波形反演的梯度方向；

$g{\left(v\right)}_{prism}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}+\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)$

其中g(v)_prism表示棱柱波形反演的梯度方向；

步骤6：用线性搜索方法或者抛物拟合方法求取更新步长α；

步骤7：利用求取的梯度方向和更新步长更新速度；

v_k＝v_k-1+α_kg(v_k-1)_prism

其中，下标k表示迭代的次数；

步骤8：利用常规全波形反演方法更新速度；

v_k+1＝v_k+α_kg(v_k)

其中，g(v_k)为第k次迭代常规全波形反演的梯度方向；

步骤9：判断v_k+1与v_k之差是否满足误差条件；

若：判断结果是v_k+1与v_k之差不满足误差条件，则执行步骤2；

或判断结果是v_k+1与v_k之差满足误差条件，则执行步骤10；

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

2.根据权利要求1所述的提高高陡构造速度反演效率的全波形联合反演方法，其特征在于：在步骤4中，具体包括

步骤4.1：定义目标函数：

$E\left(v\right)=\frac{1}{2}\underset{x_r}{\Σ}\underset{x_s}{\Σ}\underset{t}{\Σ}{\left(Ru\right(t,x_r,x_s)-p{\left(t,x_r,x_s\right)}_{obs})}^2---\left(1\right);$

其中，Ru(t,x_r,x_s)和p(t,x_r,x_s)_obs分别为数值模拟的地震记录与野外观测炮记录，x_s和x_r分别表示炮点和检波点，t为时间；

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

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\frac{1}{2}\mathrm{\δ}\[(Ru-p_{obs})\·(Ru-p_{obs})\]\\ =\[\mathrm{\δ}(Ru-p_{obs})\]\·(Ru-p_{obs})\\ =\left(R\mathrm{\δ}u\right)\·(Ru-p_{obs})\end{array}---\left(2\right);$

步骤4.3：定义二维声波波动方程：

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u=s(t,x_s)---\left(3\right);$

将方程(3)变形：

$\frac{1}{{(v+\mathrm{\δ}v)}^2}\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(4\right);$

其中，s(t,x_s)表示震源项，表示Laplace算子；

对方程(4)进行Taylor展开：

$(\frac{1}{v^2}-\frac{2\mathrm{\δ}v}{v^3})\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(5\right);$

对公式(5)展开略去高阶项:

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u+\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u-\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}=s(t,x_s)---\left(6\right);$

将方程(4)和方程(6)相减,得到：

$\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u=\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}---\left(7\right);$

进一步可得：

$\mathrm{\δ}u=L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)---\left(8\right);$

其中，L表示正传播过程；

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

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\left(R\right(L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\left)\right)\·(Ru-p_{obs})\\ =\left(L\right(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\left)\right)\·R^{*}(Ru-p_{obs})\\ =\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\·L^{*}\left(R^{*}\right(Ru-p_{obs}\left)\right)\end{array}---\left(9\right);$

其中，R^*(Ru-p_obs)表示将数据残差空间扩展到整个空间，L^*(R^*(Ru-p_obs))表示剩余波场逆时传播；

步骤4.5：令数据残差为：

Δp＝Ru-p_obs(10)；

步骤4.6：将正传波场u和数据残差Δp＝Ru-p_obs分解为如下两部分：

u(t,x_r,x_s)＝u₁(t,x_r,x_s)+u₂(t,x_r,x_s)(11)；

Δp(t,x_r,x_s)＝Δp₁(t,x_r,x_s)+Δp₂(t,x_r,x_s)(12)；

其中，u₁(t,x_r,x_s)和u₂(t,x_r,x_s)分别表示一阶反射波场和二阶棱柱波波场；Δp₁(t,x_r,x_s)和Δp₂(t,x_r,x_s)分别表示一阶反射波场的数据残差和二阶棱柱波场的数据残差；

步骤4.7：将方程(10)、(11)、(12)代入方程(9)：

步骤4.8：求取目标函数对速度模型的梯度：

$\begin{array}{c}g\left(v\right)=\frac{\mathrm{\δ}E\left(v\right)}{\mathrm{\δ}v}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1\left)\right)\\ +\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2\left)\right)+\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1\left)\right)\end{array}---\left(14\right);$

令g(v)_reflect，g(v)_prism1和g(v)_prism2分别表示一次反射波的梯度，prism1的梯度和prism2的梯度，则

$g{\left(v\right)}_{reflect}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1\left)\right)---\left(15\right);$

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1\left)\right)---\left(16\right);$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2\left)\right)---\left(17\right);$

将一阶反射波场u₁(t,x_r,x_s)由正传波场u表示，二阶棱柱波波场u₂(t,x_r,x_s)由线性正演模拟B(u)表示，一阶反射波场的数据残差的L^*(R^*(Δp₁))由反传波场w'表示，二阶棱柱波波场的数据残差的L^*(R^*(Δp₂))由反传波场w'的线性波长反传B^*(w')表示，做相应的变量替换，则

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}---\left(18\right);$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)---\left(19\right).$