一种基于薄板近似的超松弛迭代叠前偏移成像方法

本发明属于地震散射波偏移成像，具体涉及利用超松弛迭代法在薄板内求取lippmann–schwinger积分方程的格林函数解及其法向导数，特别涉及一种基于薄板近似的超松弛迭代叠前偏移成像方法。

背景技术：

1、随着地震工区从简单区到复杂区的扩展，尤其是反射地震技术的发展和亟需解决更为复杂地质问题的提出，当勘探目标地质体构造复杂、断层发育、地层倾角较陡或岩性横向突变、不同尺度的非均匀地质体共生时，会形成极为复杂的、多种波组相互干涉的地震波场。只针对反射波和绕射波不能对复杂区域进行精细成像，然而散射波也携带着与复杂构造和复杂岩性有关的几何和物理信息，如何能在利用常规反射波实现成像的同时，利用散射波实现复杂构造偏移成像就显得十分重要。

2、为了实现对散射波的偏移成像，散射波场的延拓问题转化为针对green函数的求解问题。利用扰动理论，可将描述散射问题的helmholtz方程转化成为第二类fredholm型积分方程，即lippmann schwinger方程；简记为ls方程。从而，散射波数值模拟问题就转化成为了ls方程的求解问题。从而，散射波数值模拟问题就转化成为了ls方程的求解问题。相对于其他方法：传递(传播)矩阵法、扰(摄)动法、有限单元法、几何射线法等，积分方程法因其本身具有半解析特征，且理论研究和公式推导简单、容易实现而成为散射波分析的有力工具。

3、现有技术公开了一种lippmann-schwinger积分方程广义超松弛迭代解法及其收敛特性，将广义超松弛迭代法应用于地震散波场正演模拟中，在网格间距、时间采样间隔较大的情况下可以得到与有限差分相当的数值模拟结果。现有技术还公开了一种利用离散与fft快速褶积的散射地震波并行计算方法，对原ls方程进行改写，构造等价ls方程并用法进行离散，根据格林函数矩阵的toeplitz性质，通过快速傅氏变换来加速矩阵向量乘积，达到降低存储和计算复杂的目的。然而，该方法计算效率依然低下，如何提高计算效率是一个亟需解决的问题。到目前为止，还没有将超松弛迭代解法求解散射波场方法应用到叠前偏移成像的波场延拓中去。

技术实现思路

1、本发明的目的就在于提供一种基于薄板近似的超松弛迭代叠前偏移成像方法，以解决偏移成像过程中，针对广义超松弛迭代法在解lippermann-schwinger(ls)方程时，离散系数矩阵满秩导致占用存储空间大、计算效率低的问题。

2、本发明的目的是通过以下技术方案实现的：

3、首先，对速度模型进行薄板划分；其次，在第一个薄板内利用震源子波和叠前地震记录分别进行求取正传偏移green函数的一阶超松弛迭代解gd(x,z,w)和反向的偏移green函数的一阶超松弛迭代解的法向导数b(x,z,w)；然后，将计算结果作为下一个薄板的边界条件，直至模型底部；最后，将gd(x,z,w)和b(x,z,w)做互相关成像，获得最终的偏移成像结果。

4、一种基于薄板近似的超松弛迭代叠前偏移成像方法，包括以下步骤：

5、a、读取光滑后的二维速度模型，横向网格点数为x，纵向长度为z，模型每个位置的速度为vxz，将速度模型沿着z方向划分成一个个薄板，共计m个；

6、b、设置震源参数，这里采用delta脉冲，主频为f hz，第一炮的位置为(xa,za)，其中，a代表炮数，a＝1,2,…,n，道间距h m，共计n炮，检波器r个，道间距s m，采样间隔为δt，采样长度为n，网格间距分别为dx、dz，每个薄板厚度为dz；

7、c、在第一个薄板内，利用delta脉冲在震源位置激发，作为green函数的初值g(x,zj,t)，t为传播时间，(x,zj＝1)代表空间的位置，zj为第j个薄板的深度(j＝1,2…,m)，这里的zj为初值取炮点位置(x0,z0)；

8、d、利用傅里叶变换将green函数g(x,zj,t)从时间-空间域变成频率波数域具体公式如下：

9、

10、其中，i为虚数单位，kx为横向波数，w为角频率，数值上等于2πf；

11、e、基于薄板近似理论并利用一阶超松弛迭代法求解薄板内的green函数值，即前向散射green函数场：

12、

13、其中，为傅里叶反变换。k0(zj)＝w/cj，cj为zj深度薄板内的最小速度，为介质的速度扰动，为松弛因子，其中l＝i-a，i为单位矩阵，a＝g0(x,zj+dz；x,zj；w)o(x,z)dz，其中式(2)为频率波数域和频率空间域交替进行计算，频率波数域做相移处理，频率空间域做时移校正；

14、f、将步骤e的计算结果g(x,zj+dz,w)作为下一个薄板的边界条件，重复步骤e，获得不同薄板内的前向散射green函数值；

15、g、当计算到模型底部时，循环停止，保存每个深度的不同频率的green函数值(即为正传延拓场gd(x,z,w))；

16、h、利用该震源处获得的去除直达波的地震记录d(x,z,t)进行傅里叶变换，从时间空间域变换到频率波数域最为第一个薄板的初值；

17、i、利用一阶超松弛迭代法求解薄板内反向的green函数值：

18、

19、这里由于求解使上行波，i前面取正号，由于反向递推，时间往时间反向递推，也是取负号，因此，反向递推公式与正向递推公式形式上相同，但物理意义不同，式(3)也是在频率波数域和频率空间域交替计算，

20、j、求取薄板内的green函数的法向导数：

21、

22、其中，表示法向导数，sinθ＝kx/k0；

23、k、将步骤i的计算结果d(x,zj+dz,w)，作为下一个薄板的边界条件，重复步骤九-十，获得不同薄板内的反向散射green函数值d(x,zj+dz,w)及其法向导数值b(x,zj,w)；

24、l、当计算到模型底部时，循环停止，保存每个深度的不同频率的green函数的导数值b(x,z,w)；

25、m、进行互相关成像，获得最终的偏移成像结果i1(x,z)，计算公式如下：

26、i1(x,z)＝∫dωgd(x,z,w)b(x,z,w)               (5)

27、n、将剩余的n-1炮分别按步骤b-步骤m进行循环，获得每一个单炮偏移成像结果i2(x,z),i3(x,z),…,in(x,z)；

28、o、将所有单炮偏移成像结果进行叠加获得最终的偏移成像结果i(x,z)：

29、i(x,z)＝i1(x,z)+i2(x,z)+…+in(x,z)           (6)。

30、进一步地，所述步骤e，具体为：

31、e1、将速度模型沿着z轴方向划分成m个薄板；

32、e2、当在第一个薄板内，给定detla函数作为初值将其进行傅里叶变换，从时间域转换到频率域；

33、e3、将lippmann–schwinger积分方程改写成：

34、g(x,zj,w)＝g0(xa,za,w)+ag(x,zj,w)           (7)

35、其中，其中为第一个薄板内的介质背景green函数，表达式为g0(xa,za,w)为detla函数的初值，同时将上式写成如下形式：

36、(i-a)g(x,zj,w)＝g0(xa,za,w)              (8)

37、e4、式(8)的一阶迭代解形式为：

38、

39、其中，β为超松弛因子；i-l＝a；

40、e5、为了得到精确的解，求得取剩余向量r1＝g(x,za,w)-βlg(x,za,w)的极小值，取β使||r1||＝||g(x,za,w)-βlg(x,za,w)||在hilbert空间中为极小，则有：

41、

42、其中，g0＝g(x,za,w)，为g0复共轭。满足上式的必要条件为：

43、

44、得到：

45、

46、e6、将l＝i-a代入到(9)式，得到：

47、

48、e7、利用薄板近似理论利用下面公式将结果传递到下一个薄板：

49、

50、其中，为傅里叶反变换。k0(zj)＝w/cj，cj为zj深度薄板内的最小速度，

51、

52、为介质的速度扰动，为松弛因子，其中l＝i-a，i为单位矩阵，a＝g0(x,zj+dz；x,zj；w)o(x,z)dz，其中式(2)为频率波数域和频率空间域交替进行计算，频率波数域做相移处理，频率空间域做时移校正。

53、与现有技术相比，本发明的有益效果是：

54、1、本发明将超松弛迭代解法求解散射波场方法应用到叠前偏移成像的波场延拓中去，克服了传统的born散射级数受弱散射假设限制，且只适用于短程传播的问题；

55、2、借鉴薄板近似的思想，将模型进行划分，将全局的松弛因子转换成局部松弛因子，减少存储空间，从而大大减少了计算量，在不影响精度的情况下提高了计算效率。

56、3、对基于地震记录的反传递推场求取法向导数(类似对信号进行高通滤波)，提高了偏移成像精度。