[发明专利]一种波动方程叠前偏移中的分层延拓成像方法

权利要求书

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)三维偏移距平面波方程：

(4v2∂2∂τ2+∂2∂x2+∂2∂y2)U~+4ω2v2(1-v2p→h24)U~=0---(3)]]>

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

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

∂U~∂τ=-iω1-v2p→h2/4U~---(4)]]>

∂U~∂τ=-iav24ω1-v2p→h2/4(∂2∂x2+∂2∂y2)1+bv24ω2(1-v2p→h2/4)(∂2∂x2+∂2∂y2)U~---(5)]]>

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

U~(x,y,τ+Δτ;ω)=U~(x,y,τ;ω)exp(-iω1-v2p→h24Δτ)---(6)]]>

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

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

{I-[(α+βx)-iηx]T}U~n+1={I-[(α+βx)+iηx]T}U~n---(7)]]>

{I-[(α+βy)-iηy]T}U~n+1={I-[(α+βy)+iηy]T}U~n---(8)]]>

其中：ηx=av2δτ8ωΔx21-v2p→h24,]]>ηy=av2δτ8ωΔy21-v2p→h24,]]>βx=bv24Δx2ω21-v2p→h24,]]>I＝(0，1，0)，T＝(-1，2，-1)。Δx和Δy表示x、y方向的差分间隔。α、a、b为差分方程系数。

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

Ip(x,y,τ)=ΣωIω(x,y,τ)=Σωreal(U~(x,y,τ;ω))---(9)]]>

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