CN103064110A - 一种波动方程叠前偏移中的分层延拓成像方法 - Google Patents

Abstract

本发明为一种波动方程叠前偏移中的分层延拓成像方法，先根据节点内存大小，对成像空间进行分层；后从首个单平面波数据开始，依次对每一个单平面波数据进行延拓和成像操作；直至单平面波全部各层的单频数据延拓成像完成；所有单平面波延拓成像结束输出成像结果。本发明针对常规波动方程叠前偏移方法内存需求过大及并行节点利用率不高的问题，根据计算节点可用内存空间大小自动调整偏移计算的内存需求，做到内存开支与CPU运算均衡，充分利用计算机集群的并行计算能力，实现大规模地震数据的叠前偏移成像。

Description

一种波动方程叠前偏移中的分层延拓成像方法

技术领域

本发明属于三维波动方程叠前偏移成像处理方法，针对常规方法内存需求过大及并行节点利用率不高的问题，涉及地震叠前偏移成像和高性能并行计算等，可应用于油气勘探地震资料成像处理。

背景技术

三维有限差分波动方程叠前偏移存在计算过程中内存开支过大的问题，特别是当今三维地震勘探采集数据量越来越大，内存问题成为困扰波动方程叠前偏移方法实用化的一种重要因素。假设有一个三维工区inline方向点数为Nx，crossline方向点数Ny，成像深度点数为Nτ，偏移频率空间域复数类型数据体的单频个数为Nω，那么波动方程叠前偏移主进程内存需求为4×(2Nτ+2Nω)×Nx×Ny字节，从进程内存需求为4×Nτ×Nx×Ny字节。一个满覆盖面积100km²(目前数据通常可达到1000km²以上)的地震数据波动方程偏移时，主进程所需内存空间约为9G，一个从进程所需内存空间约为4G；一个8核16G内存(主流配置)的计算节点，至多可用4核进行偏移计算，如此造成并行节点无法得到有效利用。数据量越大，上述问题越严重。当数据量为上述情况的两三倍时，目前主流配置的计算机集群已无法进行波动方程叠前偏移运算。

发明内容

本发明为了解决现有技术中存在的非分层延拓成像方法占用内容过大，无法进行操作和偏移运算的技术问题，研发了一种波动方程叠前偏移中的分层延拓成像方法。本发明针对常规波动方程叠前偏移方法内存需求过大及并行节点利用率不高的问题，根据计算节点可用内存空间大小自动调整偏移计算的内存需求，做到内存开支与CPU运算均衡，充分利用计算机集群的并行计算能力，实现大规模地震数据的叠前偏移成像。

本发明为了解决上述现有技术中存在的技术问题，采用的技术方案如下，

一种波动方程叠前偏移中的分层延拓成像方法，所述延拓成像方法先根据节点内存大小，对成像空间进行分层；后从首个单平面波数据开始，依次对每一个单平面波数据进行延拓和成像操作；逐层延拓成像操作中包括主进程和从进程，所述主进程负责读写地震数据和向从进程分发数据，从进程负责该单层单频数据的延拓成像；各单层计算结束后，从进程把延拓至本层深度后的单频数据送回到主进程；由主进程再把数据发送出去；并进行下一层的延拓成像操作，直至单平面波全部各层的单频数据延拓成像完成；所有单平面波延拓成像结束输出成像结果。

具体的，所述的一种波动方程叠前偏移中的分层延拓成像方法，如下步骤：

步骤1，读取输入输出文件名、成像范围、最大可用内存等参数

步骤2，成像空间进行分层：

设每个从进程最大可用内存为Memory_max兆，假设层厚度以100点数为最小单位，如果有整数Nt_layer满足：

4×Nt_layer×100×Nx×Ny≤Memory_max×1024×1024         (1)

4×(Nt_layer+1)×100×Nx×Ny＞Memory_max×1024×1024      (2)那么分层延拓的层厚度为Nt_layer×100。如果Nτ/(100×Nt_layer)为整数，则该值为延拓成像的层数，记为N_layer；如果Nτ/(100×Nt_layer)不为整数，则延拓成像的层数为其值取整加1。

Nx、Ny、Nτ依次是成像空间在三个坐标x，y，τ方向上的点数。Nt_layer等于层厚度除以100，若Nt_layer等于10，表示每层的层厚度为1000。N_layer是层数，假设成像深度τ方向共3000点，每层的厚度为1000个点，则层数N_layer＝3。

步骤3，单平面波P循环，逐一对各个单平面波数据进行延拓和成像；

步骤4，单个平面波数据逐层延拓成像

(1)判断：是否是第一层：

如果是第一层，主进程从磁盘读取波场数据，利用快速傅里叶变换(FFT)将数据转换到频率空间域后，再利用并行库MPI_SEND函数将单频数据

(τ＝0表示地表，ω₁表示一个频率)发送到从进程；如果不是第一层进入(2)；

(2)主进程则利用并行库MPI_RECV函数接收从进程完成层延拓后的波场数据，再利用并行库MPI_SEND函数将单频数据

(τ＝τ₁表示某延拓后的深度)发送到从进程；

(3)从进程接收主进程发送来的单频数据，按有限差分解进行波场延拓，再计算单频的成像值I_ω1(x，y，τ)；

(4)单层延拓和成像结束后，利用并行库MPI_REDUCE函数规约累加所有从进程的单频数据成像结果，并计入磁盘；

步骤5，判断是否层循环结束，即各层全部延拓成像是否结束，如果否转入步骤4，如果是转入步骤6；

步骤6，进行下一个平面波数据偏移，判断是否单面波循环结束，否转入步骤3，是转入步骤7；

步骤7，结束，输出波动方程偏移成像数据体结果。

其中，MPI_SEND是并行库MPI(Message Passing Interface)中定义的一个消息发送函数，用于将一个消息(变量、数组等)从当前进程发送至目标进程。MPI_RECV是并行库MPI(Message Passing Interface)中定义的一个消息接收函数，用于将指定源进程的一个消息接收到当前进程。MPI_REDUCE是并行库MPI(Message Passing Interface)中定义的一个全局归约操作函数，用于将所有进程中的同一数据按设定的算术运算符(加减等)归约到一个指定进程。所述在步骤4的(3)中，

有限差分解进行波场延拓和成像采用：

(31)三维偏移距平面波方程：

$(\frac{4}{v^2}\frac{{\∂}^2}{{\∂\mathrm{\τ}}^2}+\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})\tilde{U}+4\frac{{\mathrm{\ω}}^2}{v^2}(1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4})\tilde{U}=0---\left(3\right)$

其中，x和y是平面坐标，τ是时间深度坐标，ω是圆频率，

是表示平面波的射线参数，为频率空间域波场，v＝v(x，y，τ)为时间域介质层速度。

(32)方程(3)可以分裂为以下两个方程：

$\frac{\∂\tilde{U}}{\∂\mathrm{\τ}}=-\mathrm{i\ω}\sqrt{1-v^2{\stackrel{\→}{p}}_h^2/4}\tilde{U}---\left(4\right)$

$\frac{\∂\tilde{U}}{\∂\mathrm{\τ}}=-i\frac{\frac{{\mathrm{av}}^2}{4\mathrm{\ω}\sqrt{1-v^2{\stackrel{\→}{p}}_h^2/4}}(\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})}{1+\frac{{\mathrm{bv}}^2}{4{\mathrm{\ω}}^2(1-v^2{\stackrel{\→}{p}}_h^2/4)}(\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})}\tilde{U}---\left(5\right)$

(33)方程(4)的解析解为：

$\tilde{U}(x,y,\mathrm{\τ}+\mathrm{\Δ\τ};\mathrm{\ω})=\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω})\exp (-\mathrm{i\ω}\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}\mathrm{\Δ\τ})---\left(6\right)$

其中，表示向下延拓步长Δτ后的波场。

(34)方程(5)可以表示如下x、y两个方向的差分方程，采用追赶法求解本方程，即实现波场由

到延拓。

$\{I-[(\mathrm{\α}+{\mathrm{\β}}_x)-{\mathrm{i\η}}_x\left]T\right\}{\tilde{U}}^{n+1}=\{I-[(\mathrm{\α}+{\mathrm{\β}}_x)+{\mathrm{i\η}}_x\left]T\right\}{\tilde{U}}^n---\left(7\right)$

$\{I-[(\mathrm{\α}+{\mathrm{\β}}_y)-{\mathrm{i\η}}_y\left]T\right\}{\tilde{U}}^{n+1}=\{I-[(\mathrm{\α}+{\mathrm{\β}}_y)+{\mathrm{i\η}}_y\left]T\right\}{\tilde{U}}^n---\left(8\right)$

其中： ${\mathrm{\η}}_x=\frac{{\mathrm{av}}^2\mathrm{\δ\τ}}{8\mathrm{\ω\Δ}{x}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$ ${\mathrm{\η}}_y=\frac{{\mathrm{av}}^2\mathrm{\δ\τ}}{8\mathrm{\ω\Δ}{y}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$ ${\mathrm{\β}}_x=\frac{{\mathrm{bv}}^2}{4\mathrm{\Δ}{x}^2{\mathrm{\ω}}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$ I＝(0，1，0)，T＝(-1，2，-1)。Δx和Δy表示x、y方向的差分间隔。α、a、b为差分方程系数。

(35)单平面偏移成像值可由频率空间域波场表示为：

$I_p(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}{I}_{\mathrm{\ω}}(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{real}\left(\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω})\right)---\left(9\right)$

real表示取实部。

本发明有效降低了波动方程叠前偏移中的内存需求，若总延拓成像深度点数为N_all，分层延拓中的层厚度为Nt_layer×100，则主进程内存需求减小到非分层整体延拓成像的(Nt_layer×100+Nω)/(N_all+Nω)倍，从进程内存需求减小到非分层整体延拓成像的Nt_layer×100/N_all倍。通过减小从进程内存需求，可以利用节点中更多的核参与偏移计算，提高了并行节点的有效利用率。本发明与非分层整体延拓成像相比，在基本不增加额外运行时间的条件下(仅增加部分节点通信时间)，使得波动方程叠前时间偏移具备了大规模地震数据处理能力。

附图说明

图1是本发明方法的流程图；

图2是本发明实施例1的效果输出图

图3是本发明实施例2的效果输出图

将结合具体的实施方式加以说明

具体实施方式

本发明中，

(1)平面波波动方程叠前偏移方法

由时间域双平方根算子的频散关系可推导出

$\frac{4}{v^2}{k}_{\mathrm{\τ}}^2+{\stackrel{\→}{k}}_m^2-4\frac{{\mathrm{\ω}}^2}{v^2}{\cos }^2\mathrm{\γ}=0---\left(1\right)$

其中，k_τ为双程旅行时对应的波数，

γ为炮检点射线的半张角，v为常速度。上式在频率-空间域可以表示为：

$(\frac{4}{v^2}\frac{{\∂}^2}{{\∂\mathrm{\τ}}^2}+\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})\tilde{U}+4\frac{{\mathrm{\ω}}^2}{v^2}{\cos }^2\mathrm{\γ}\tilde{U}=0---\left(2\right)$

其中，v＝v(x，y，τ)， $\tilde{U}=\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω}).$

γ与平面波矢量

有如下关系式：

$\cos \mathrm{\γ}=\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}---\left(3\right)$

进而得三维偏移距平面波方程：

$(\frac{4}{v^2}\frac{{\∂}^2}{{\∂\mathrm{\τ}}^2}+\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})\tilde{U}+4\frac{{\mathrm{\ω}}^2}{v^2}(1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4})\tilde{U}=0---\left(4\right)$

对上式利用有限差分算法可进行单平面波的波场延拓。

利用零时间成像条件，单平面偏移成像值可用频率空间域波场表示为：

$I_p(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}{I}_{\mathrm{\ω}}(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{real}\left(\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω})\right)---\left(5\right)$

real表示取实部。

逐一对每一个单平面波

进行延拓和成像，可以实现平面波波动方程偏移成像。

(2)分层延拓层厚度计算方法

波动方程偏移中，主进程负责从硬盘读写数据和向从进程分发作业，内存数组包括一个Nω×Nx×Ny的复数类型的三维数组和两个Nτ×Nx×Ny的实数类型的三维数组，这两个实数数组其中一个存放单频成像数据体，另一个存放最终成像数据体；从进程负责对接收自主进程的单频数据进行延拓处理，内存数组是一个Nτ×Nx×Ny的实数类型的三维数组。

在计算机集群运行波动方程偏移时，一般主进程单独使用一个节点，可用内存空间较大，多个从进程共用一个节点，共享一个节点的内存空间，因此利用从进程最大可用内存进行分层控制。记每个从进程最大可用内存为Memory_max兆，假设层厚度以100点数为最小单位，如果有整数Nt_layer满足：

4×Nt_layer×100×Nx×Ny≤Memory_max×1024×1024              (6)

4×(Nt_layer+1)×100×Nx×Ny＞Memory_max×1024×1024          (7)

那么分层延拓的层厚度为Nt_layer×100。如果Nτ/(100×Nt_layer)为整数，则该值为延拓成像的层数，记为N_layer；如果Nτ/(100×Nt_layer)不为整数，则延拓成像的层数为其值取整加1。

如图1所述的方法包括：

步骤1，读取输入输出文件名、成像范围、最大可用内存等参数

步骤2，成像空间进行分层：

设每个从进程最大可用内存为Memory_max兆，假设层厚度以100点数为最小单位，如果有整数Nt_layer满足：

4×Nt_layer×100×Nx×Ny≤Memory_max×1024×1024               (1)

4×(Nt_layer+1)×100×Nx×Ny＞Memory_max×1024×1024           (2)

那么分层延拓的层厚度为Nt_layer×100。如果Nτ/(100×Nt_layer)为整数，则该值为延拓成像的层数，记为N_layer；如果Nτ/(100×Nt_layer)不为整数，则延拓成像的层数为其值取整加1。

Nx、Ny、Nτ依次是成像空间在三个坐标x，y，τ方向上的点数。Nt_layer等于层厚度除以100，若Nt_layer等于10，表示每层的层厚度为1000。N_layer是层数，假设成像深度τ方向共3000点，每层的厚度为1000个点，则层数N_layer＝3。

步骤3，单平面波P循环，逐一对各个单平面波数据进行延拓和成像；

步骤4，单个平面波数据逐层延拓成像

(1)判断：是否是第一层：

如果是第一层，主进程从磁盘读取波场数据，利用快速傅里叶变换(FFT)将数据转换到频率空间域后，再利用并行库MPI_SEND函数将单频数据

(τ＝0表示地表，ω₁表示一个频率)发送到从进程；如果不是第一层进入(2)；

(2)主进程则利用并行库MPI_RECV函数接收从进程完成层延拓后的波场数据，再利用并行库MPI_SEND函数将单频数据(τ＝τ₁表示某延拓后的深度)发送到从进程；

(3)从进程接收主进程发送来的单频数据，按有限差分解进行波场延拓，再计算单频的成像值I_ω1(x，y，τ)；

(4)单层延拓和成像结束后，利用并行库MPI_REDUCE函数规约累加所有从进程的单频数据成像结果，并计入磁盘；

步骤5，判断是否层循环结束，即各层全部延拓成像是否结束，如果否转入步骤4，如果是转入步骤6；

步骤6，进行下一个平面波数据偏移，判断是否单面波循环结束，否转入步骤3，是转入步骤7；

步骤7，结束，输出波动方程偏移成像数据体结果。

其中，MPI_SEND是并行库MPI(Message Passing Interface)中定义的一个消息发送函数，用于将一个消息(变量、数组等)从当前进程发送至目标进程。MPI_RECV是并行库MPI(Message Passing Interface)中定义的一个消息接收函数，用于将指定源进程的一个消息接收到当前进程。MPI_REDUCE是并行库MPI(Message Passing Interface)中定义的一个全局归约操作函数，用于将所有进程中的同一数据按设定的算术运算符(加减等)归约到一个指定进程。

所述在步骤4的(3)中，

有限差分解进行波场延拓和成像采用：

(31)三维偏移距平面波方程：

$(\frac{4}{v^2}\frac{{\∂}^2}{{\∂\mathrm{\τ}}^2}+\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})\tilde{U}+4\frac{{\mathrm{\ω}}^2}{v^2}(1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4})\tilde{U}=0---\left(3\right)$

其中，x和y是平面坐标，τ是时间深度坐标，ω是圆频率，

是表示平面波的射线参数，

为频率空间域波场，v＝v(x，y，τ)为时间域介质层速度。

(32)方程(3)可以分裂为以下两个方程：

$\frac{\∂\tilde{U}}{\∂\mathrm{\τ}}=-\mathrm{i\ω}\sqrt{1-v^2{\stackrel{\→}{p}}_h^2/4}\tilde{U}---\left(4\right)$

$\frac{\∂\tilde{U}}{\∂\mathrm{\τ}}=-i\frac{\frac{{\mathrm{av}}^2}{4\mathrm{\ω}\sqrt{1-v^2{\stackrel{\→}{p}}_h^2/4}}(\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})}{1+\frac{{\mathrm{bv}}^2}{4{\mathrm{\ω}}^2(1-v^2{\stackrel{\→}{p}}_h^2/4)}(\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})}\tilde{U}---\left(5\right)$

(33)方程(4)的解析解为：

$\tilde{U}(x,y,\mathrm{\τ}+\mathrm{\Δ\τ};\mathrm{\ω})=\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω})\exp (-\mathrm{i\ω}\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}\mathrm{\Δ\τ})---\left(6\right)$

其中，

表示

向下延拓步长Δτ后的波场。

(34)方程(5)可以表示如下x、y两个方向的差分方程，采用追赶法求解本方程，即实现波场由到

延拓。

$\{I-[(\mathrm{\α}+{\mathrm{\β}}_x)-{\mathrm{i\η}}_x\left]T\right\}{\tilde{U}}^{n+1}=\{I-[(\mathrm{\α}+{\mathrm{\β}}_x)+{\mathrm{i\η}}_x\left]T\right\}{\tilde{U}}^n---\left(7\right)$

$\{I-[(\mathrm{\α}+{\mathrm{\β}}_y)-{\mathrm{i\η}}_y\left]T\right\}{\tilde{U}}^{n+1}=\{I-[(\mathrm{\α}+{\mathrm{\β}}_y)+{\mathrm{i\η}}_y\left]T\right\}{\tilde{U}}^n---\left(8\right)$

其中： ${\mathrm{\η}}_x=\frac{{\mathrm{av}}^2\mathrm{\δ\τ}}{8\mathrm{\ω\Δ}{x}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$ ${\mathrm{\η}}_y=\frac{{\mathrm{av}}^2\mathrm{\δ\τ}}{8\mathrm{\ω\Δ}{y}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$ ${\mathrm{\β}}_x=\frac{{\mathrm{bv}}^2}{4\mathrm{\Δ}{x}^2{\mathrm{\ω}}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$

I＝(0，1，0)，T＝(-1，2，-1)。Δx和Δy表示x、y方向的差分间隔。α、a、b为差分方程系数。

(35)单平面偏移成像值可由频率空间域波场表示为：

$I_p(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}{I}_{\mathrm{\ω}}(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{real}\left(\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω})\right)---\left(9\right)$

real表示取实部。

图2和图3分别表示实施例1和实施例2

其中实施例1，南方某资料满覆盖面积为54.17km²，面元大小为25m×25m，采样时间为6s，采样间隔为2ms，采样点数为3000，CMP道集数据量为70G。采用非分层整体延拓成像方法，主进程内存需求为2457M，每个从进程的内存需求为1006M。采用本发明的分层延拓成像方法，将从进程最大可用内存空间设置为350M，偏移程序将成像深度分为3层，每层厚度为1000个采样点。分层延拓成像结果正确，图2是该资料利用本发明获得的一条线的成像剖面。

实施例2，西部某高密度资料，面元大小为7.5m×7.5m，采样时间为7s，采样间隔为2ms，采样点数为3500。对面积30km²的部分地震数据做波动方程叠前偏移处理，处理范围是360(线)×1450(CDP)，成像深度为3000个点数。采用非分层整体延拓成像方法，主进程内存需求为12.5G，每个从进程的内存需求为5.83G，由于内存需求很大，一个8核16G内存的计算节点只能用2核进行偏移计算，并行节点利用率很低。采用本发明的分层延拓成像方法，为使一个节点的8核全部参与计算，将从进程最大可用内存空间设置为1.6G，则偏移程序将成像深度分为4层，第1到第3层的厚度为800个点，第4层厚度为600个点。该资料分层延拓成像结果见图3，同样正确。

Claims (3)

1.一种波动方程叠前偏移中的分层延拓成像方法，其特征在于，所述延拓成像方法先根据节点内存大小，对成像空间进行分层；后从首个单平面波数据开始，依次对每一个单平面波数据进行延拓和成像操作；在逐层延拓成像操作中包括主进程和从进程，所述主进程负责读写地震数据和向从进程分发数据，从进程负责该单层单频数据的延拓成像；各单层计算结束后，从进程把延拓至本层深度后的单频数据送回到主进程；由主进程再把数据发送出去；并进行下一层的延拓成像操作，直至单平面波全部各层的单频数据延拓成像完成；所有单平面波延拓成像结束输出成像结果。

2.根据权利要求1所述的一种波动方程叠前偏移中的分层延拓成像方法，其特征在于，所述方法包括如下步骤：

步骤1，读取输入输出文件名、成像范围、最大可用内存参数；

步骤2，成像空间进行分层：

设每个从进程最大可用内存为Memory_max兆，假设层厚度以100点数为最小单位，如果有整数Nt_layer满足：

4×Nt_layer×100×Nx×Ny≤Memory_max×1024×1024         (1)

4×(Nt_layer+1)×100×Nx ×Ny＞Memory_max×1024×1024    (2)

那么分层延拓的层厚度为Nt_layer×100。如果Nτ/(100×Nt_layer)为整数，则该值为延拓成像的层数，记为N_layer；如果Nτ/(100×Nt_layer)不为整数，则延拓成像的层数为其值取整加1；

其中，Nx、Ny、Nτ依次是成像空间在三个坐标x，y，τ方向上的点数；Nt_layer等于层厚度除以100，若Nt_layer等于10，表示每层的层厚度为1000；N_layer是层数，假设成像深度τ方向共3000点，每层的厚度为1000个点，则层数N_layer＝3；

步骤3，单平面波P循环，逐一对各个单平面波数据进行延拓和成像；

步骤4，单个平面波数据逐层延拓成像

(1)判断：是否是第一层：

如果是第一层，主进程从磁盘读取波场数据，利用快速傅里叶变换(FFT)将数据转换到频率空间域后，再利用并行库MPI_SEND函数将单频数据

(τ＝0表示地表，ω₁表示一个频率)发送到从进程；如果不是第一层进入(2)；

(2)主进程则利用并行库MPI_RECV函数接收从进程完成层延拓后的波场数据，再利用并行库MPI_SEND函数将单频数据(τ＝τ₁表示某延拓后的深度)发送到从进程；

(3)从进程接收主进程发送来的单频数据，按有限差分解进行波场延拓，再计算单频的成像值I_ω1(x，y，τ)；

(4)单层延拓和成像结束后，利用并行库MPI_REDUCE函数规约累加所有从进程的单频数据成像结果，并计入磁盘；

步骤5，判断是否层循环结束，即各层全部延拓成像是否结束，如果否转入步骤4，如果是转入步骤6；

步骤6，进行下一个平面波数据偏移，判断是否单面波循环结束，否转入步骤3，是转入步骤7；

步骤7，结束，输出波动方程偏移成像数据体结果。

其中，MPI_SEND是并行库MPI(Message Passing Interface)中定义的一个消息发送函数，用于将一个消息(变量、数组等)从当前进程发送至目标进程；MPI_RECV是并行库MPI(Message Passing Interface)中定义的一个消息接收函数，用于将指定源进程的一个消息接收到当前进程；MPI_REDUCE是并行库MPI(Message Passing Interface)中定义的一个全局归约操作函数，用于将所有进程中的同一数据按设定的算术运算符(加减等)归约到一个指定进程。

3.根据权利要求2所述的一种波动方程叠前偏移中的分层延拓成像方法，其特征在于，所述在步骤4的(3)中，

有限差分解进行波场延拓和成像采用：

(31)三维偏移距平面波方程：

$(\frac{4}{v^2}\frac{{\∂}^2}{{\∂\mathrm{\τ}}^2}+\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})\tilde{U}+4\frac{{\mathrm{\ω}}^2}{v^2}(1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4})\tilde{U}=0---\left(3\right)$

其中，x和y是平面坐标，τ是时间深度坐标，ω是圆频率，

是表示平面波的射线参数，

为频率空间域波场，v＝v(x，y，τ)为时间域介质层速度。

(32)方程(3)可以分裂为以下两个方程：

$\frac{\∂\tilde{U}}{\∂\mathrm{\τ}}=-\mathrm{i\ω}\sqrt{1-v^2{\stackrel{\→}{p}}_h^2/4}\tilde{U}---\left(4\right)$

$\frac{\∂\tilde{U}}{\∂\mathrm{\τ}}=-i\frac{\frac{{\mathrm{av}}^2}{4\mathrm{\ω}\sqrt{1-v^2{\stackrel{\→}{p}}_h^2/4}}(\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})}{1+\frac{{\mathrm{bv}}^2}{4{\mathrm{\ω}}^2(1-v^2{\stackrel{\→}{p}}_h^2/4)}(\frac{{\∂}^2}{{\∂x}^2}+\frac{{\∂}^2}{{\∂y}^2})}\tilde{U}---\left(5\right)$

(33)方程(4)的解析解为：

$\tilde{U}(x,y,\mathrm{\τ}+\mathrm{\Δ\τ};\mathrm{\ω})=\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω})\exp (-\mathrm{i\ω}\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}\mathrm{\Δ\τ})---\left(6\right)$

其中，

表示

向下延拓步长Δτ后的波场。

(34)方程(5)可以表示如下x、y两个方向的差分方程，采用追赶法求解本方程，即实现波场由

到

的向下延拓。

$\{I-[(\mathrm{\α}+{\mathrm{\β}}_x)-{\mathrm{i\η}}_x\left]T\right\}{\tilde{U}}^{n+1}=\{I-[(\mathrm{\α}+{\mathrm{\β}}_x)+{\mathrm{i\η}}_x\left]T\right\}{\tilde{U}}^n---\left(7\right)$

$\{I-[(\mathrm{\α}+{\mathrm{\β}}_y)-{\mathrm{i\η}}_y\left]T\right\}{\tilde{U}}^{n+1}=\{I-[(\mathrm{\α}+{\mathrm{\β}}_y)+{\mathrm{i\η}}_y\left]T\right\}{\tilde{U}}^n---\left(8\right)$

其中： ${\mathrm{\η}}_x=\frac{{\mathrm{av}}^2\mathrm{\δ\τ}}{8\mathrm{\ω\Δ}{x}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$ ${\mathrm{\η}}_y=\frac{{\mathrm{av}}^2\mathrm{\δ\τ}}{8\mathrm{\ω\Δ}{y}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$ ${\mathrm{\β}}_x=\frac{{\mathrm{bv}}^2}{4\mathrm{\Δ}{x}^2{\mathrm{\ω}}^2\sqrt{1-\frac{v^2{\stackrel{\→}{p}}_h^2}{4}}},$

I＝(0，1，0)，T＝(-1，2，-1)。Δx和Δy表示x、y方向的差分间隔。α、a、b为差分方程系数。

(35)单平面偏移成像值可由频率空间域波场表示为：

$I_p(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}{I}_{\mathrm{\ω}}(x,y,\mathrm{\τ})=\underset{\mathrm{\ω}}{\mathrm{\Σ}}\mathrm{real}\left(\tilde{U}(x,y,\mathrm{\τ};\mathrm{\ω})\right)---\left(9\right)$

real表示取实部，逐一对每一个单平面波进行延拓和成像，可以实现平面波波动方程偏移成像。

Images