[发明专利]一种基于OpenMP的频率域高斯束格林函数计算方法

权利要求书

1.一种基于OpenMP的频率域高斯束格林函数计算方法，其特征在于，包括以下步骤：

步骤1：生成频率域格林函数，具体步骤如下所示：

(1)申请构体数组Cells(Nz,Nx),用于存储单个方向高斯束对应的格林函数，申请格林函数数组Green(Nz,Nx,Nω,Nshot)，用于存储格林函数，z和x指空间点二维坐标，Nz,Nx分别为z方向和x方向上的网格点数，ω为频率，Nω指频率的循环次数，Nshot指炮点的循环次数，shot指炮点；

(2)进行炮点循环，即对炮点shot进行Nshot次循环，对此循环利用OpenMP并行处理，私有化结构体数组Cells(Nz,Nx)及涉及的循环相关变量和临时变量,共享格林函数数组Green(Nz,Nx,Nω,Nshot)；

(3)进行射线方向循环，即对射线高斯束进行Np个方向的循环，对某一个方向Ip的射线高斯束进行时空域射线高斯束计算，形成该方向的Cells(Nz,Nx)，并对射线高斯束有值位置做好标记；

(4)进行频率循环，即对频率ω进行Nω次循环，对某一频率Iω的有标记的射线高斯束进行计算，得到频率域格林函数相关信息；

(5)进行模型空间循环，模型空间大小为Nx*Nz，对模型空间内第Ix和Iz的网格点，进行格林函数的积分计算，即得到Green(Iz,Ix,Iω,Ishot)，其中I＝(1，N)；z和x指空间点二维坐标；

步骤2：利用格林函数Green(Nz,Nx,Nω,Nshot)实现地震波Born正演的应用即反偏移，具体步骤如下所示：

(1)申请单道数据数组Trace(Lt)，申请数据体数组Dat(Lt,Nreceiver,Nshot)，Lt指离散地震记录时间t长度，Nshot和Nreceiver分别为炮点和炮检点的个数，shot为炮点，receiver为地震波检波器；

(2)进行炮点循环，即对编号范围从1到Nshot多炮地震数据进行计算，共计算Nshot次，利用OpenMP对此循环进行并行处理，私有Trace(Lt)及涉及的循环相关变量和临时变量,共享数据体数组Dat(Lt,Nreceiver,Nshot)；Lt指离散地震记录时间t长度，shot指炮点，receiver为地震波检波器，Nshot和Nreceiver分别为炮点和炮检点的个数；

(3)进行检波点循环，即对任意检测地震波的检波器Ireceiver进行Nreceiver次循环；

(4)进行频率循环，即对地震波有效频带范围内任一频率Iω进行Nω次循环；

(5)进行模型空间循环，模型空间大小为Nx*Nz，对地震波速度模型空间内第Ix和Iz的网格点利用已知格林函数Green(Iz,Ix,Iω,Ishot)进行积分，得到Trace(Iω)；

(6)对Trace(Iω)数据进行反傅里叶变换得到Trace(Lt)并赋值给Dat(Lt,Ireceiver,Ishot)，至此完成了对格林函数Green(Iz,Ix,Iω,Ishot)的地震波Born正演。

2.根据权利要求1所述的一种基于OpenMP的频率域高斯束格林函数计算方法，其特征在于，步骤1所述频率域格林函数为：

其中，x表示成像点坐标，x_s表示炮点坐标，表示炮点局部平面波的初始波矢量方向，ω为地震波频率，u_GB(x,x_s,p^s,ω)表示笛卡尔坐标系中一个高斯束，单个高斯束可以在射线中心坐标系由动力学射线追踪求取，i表示复数。

3.根据权利要求2所述的一种基于OpenMP的频率域高斯束格林函数计算方法，其特征在于，所述高斯束表示如下：

其中，s为高斯束射线中心坐标系沿射线方向坐标，q^T＝(q₁,q₂)为垂直射线方向坐标，v(s)和τ(s)表示沿射线方向的速度和时间，ω为地震波频率，P(s)和Q(s)为沿射线方向的两个动力学射线追踪参数，Q(s₀)为Q(s)的初值，v(s₀)为速度v(s)的初值，该射线中心坐标系的高斯束可以通过坐标变换变成笛卡尔坐标系下的高斯束。

4.根据权利要求1所述的一种基于OpenMP的频率域高斯束格林函数计算方法，其特征在于，所述地震波Born正演即反偏移公式表示如下：

D_Born(x_r,x_s,ω)＝2ω²∫_Dm(x)G(x,x_s,ω)G(x,x_r,ω)dx

上式中，m(x)和D_Born(x_r,x_s,ω)分别表示地下模型空间x的像和地表炮点位于(x_r,x_s)的地震波Born正演数据，x_r表示检波点坐标，x_s表示炮点坐标，z_r表示地震波检波点深度，其中，G(x,x_s,ω)和G(x,x_r,ω)分别是炮点和检波点的格林函数，ω为地震波频率，G^*(x,x_r,w)是G(x,x_r,ω)的共轭，D_s(x_r,x_s,ω)指记录的频率域地震数据，D表示数据。