CN103278848A - 基于mpi并行预条件迭代的地震成像正演方法 - Google Patents

Abstract

本发明公开了一种基于MPI并行预条件迭代的地震成像正演方法，具体为：在待勘测区域中激发人工地震波，并由频率域波动方程模型-ΔU(x,y,z)-k²U(x,y,z)=G(x,y,z)来刻画，其中

为拉普拉斯算子，U为待勘测区域的连续空间压力场，k为波数，G为连续的震源项；正演过程中，采用并行预条件迭代法快速有效地求解上述方程，即求解出勘测区域内的压力场，模拟地震波的传播，从而为地下层介质图像的反演提供坚实的基础。该法能够利用高性能多核计算机进行快速计算，并在相同规模下对内存的需求较低。

Description

基于MPI并行预条件迭代的地震成像正演方法

技术领域

本发明涉及地球物理勘探中的地震波数值模拟领域，更具体地，涉及一种基于MPI并行预条件迭代的地震成像正演方法。

背景技术

在地球物理、航天航空、海洋技术等工程技术领域中，频率域波动方程用来刻画声波、地震波等波的传播。如在地震勘察和油气勘探中,有关地下结构的地震信息十分重要.地下结构的精确地震信息有助于寻找可能的地下层油藏。这种地球的地下图像信息，是通过收集人工地震波从不同物理属性的地下层界面反射回地面的频次得到的。在这个过程中频率域波动方程正演（即求解频率域波动方程）是一项关键基础，频率域波动方程的求解最终归结为大规模线性系统的求解。当前求解频率域波动方程的一般直接法和迭代法。直接法一般基于高斯消去法求解，而迭代法常常采用预条件处理技术通过反复的矩阵向量积来实现求解，这两种求解方法都是通过计算机实现。

频率域波动方程离散后的线性系统规模非常大，而现有的直接法和预条件迭代法存在一些问题。直接法的缺点就是受到了计算机内存的限制。虽然当前超级计算机发展迅速，但直接法需要的计算机内存十分巨大。当处理的数据规模很大时，直接法难以承受，即使是并行的直接法。迭代法通过反复的矩阵向量积来实现求解，能够有效减少计算机内存使用量，当前使用的迭代法如最有效的多重网格预条件Krylov子空间方法是串行的方式。相对于并行算法，串行算法的求解速度是很慢的。在当前高性能超级计算机迅速发展的形势下，串行算法日益不能满足超大规模计算的实际需求。

发明内容

本发明为克服上述现有技术中在获取地下图像时，对频率域波动方程正演时采用直接法内存需求大以及采用串行迭代法速度慢的不足，提出一种基于MPI并行预条件迭代的地震成像正演方法，该法能快速模拟波传播从而快速反演地下层图像，并在同等规模下对计算机的内存需求较低。

为了实现上述目的，本发明的技术方案为：

一种基于MPI并行预条件迭代的地震成像正演方法，在待勘测区域中人工激发地震波，通过频率域波动方程模型-ΔU(x,y,z)-k²U(x,y,z)=G(x,y,z)来模拟波传播；其中

为拉普拉斯算子，U为待勘测区域的连续空间压力场，k为波数，G为连续的震源项；对频率域波动方程进行并行预条件正演求解，求解出勘测区域内的压力场，再根据压力场进行地震信息反演，获取地球的地下层图像；所述对频率域波动方程进行并行预条件正演包括以下步骤：

S1.采用差分法对频率域波动方程进行离散，得到线性系统Au=g；其中系数矩阵A为稀疏的；右端项向量g是震源项的离散，未知向量u为压力场；

S2.对线性系统Au=g进行预条件处理，得到等价的线性系统；

AM^-1v=g，Mu=v，               （1）

其中M是基于移位拉普拉斯算子离散得到的预条件子，和A具有同样的稀疏结构；

S3.对待勘测区域进行并行划分；

S4.对步骤S3划分的子区域进行内外迭代，从而求解出线性系统（1）在各子区域内的局部解，最后将这些局部解合并，得到频率域波动方程的整体解，整体解即为离散的压力场u。

为了更好的更快地利用超级计算机进行频率域波动方程(Helmholtz方程)的正演，为地震波油气勘探等工业技术领域提供重要的应用基础。

更进一步的，所述步骤S1采用的差分法为基于频散极小化的有限差分法。

具体的，采用基于频散极小化的27点有限差分法离散频率域波动方程，得到线性系统Au=g，其中系数矩阵A为稀疏的，并呈27对角的形状；人工地震波k为输入，另一输入数据为G=δ(x-x₁,y-y₁,z-z₁)，即G在计算区域内的某点(x₁,y₁,z₁)为1其余为0，(x₁,y₁,z₁)可以在区域内任选。采用基于频散极小化的27点有限差分法离散频率域波动方程的离散误差小。

更进一步的，所述步骤S2中M采用基于频散极小化差分法离散获得。

更进一步的，所述步骤S3对待勘测区域进行并行网格划分的具体方式为：

S31.在三维空间的每个维度上进行2ⁿ分割，n>0为自然数，设三维空间中x,y,z方向对应的剖分分别为p_x，p_y，p_z，则待勘测区域被分成Np=p_x·p_y·p_z子区域，即待勘测区域中数据被分划成Np个子区域，所述Np个子区域采用Np个进程分别处理；

S32.确定子区域的边界，待勘测区域的计算剖分分别为Nx，Ny，Nz，即网格规模为N=Nx·Ny·Nz，其中Nx，Ny，Nz取为2的幂次方加1；在界定分化后的子区域边界时进行多重网格的并行实施；

S33.多重网格FMG的重要组成元素有各层的网格算子A^h和网格间的转移算子，即细层到粗层的限制算子和粗层到细层的延拓算子；在进行多重网格预条件Bi-CGStab法求解时，首先要构造其组成元素，令A^h表示细网格算子，

表示细网格到粗网格的限制算子，

表示粗网格到细网格的延拓算子，则粗网格算子为

采用全粗化的多重网格，即在设三维空间的x,y,z三个方向都进行粗化，并利用全加权限制算子和黑盒子延拓算子，即延拓算子的构造是基于离散的稀疏矩阵；

S34.各子区域间的数据通信，数据通信通过非阻塞通信MPI_Isend和MPI_Irecv来接收和发送，各子区域共享的数据通过MPI_Allreduce收集。

更进一步的，所述步骤S4对格子区域进行内外迭代为：采用Krylov子空间方法作为外迭代来求解线性系统（1），并用瀑布型多重网格FMG作为内迭代来求解子系统Mu=v，作为子系统。

更进一步的，所述采用双共轭梯度法Bi-CGStab作为外迭代来求解线性系统（1）。

更进一步的，所述瀑布型多重网格FMG作为内迭代来求解子系统Mu=v时只迭代一次。因为Mu=v不需要精确求解，虽然越精确其预条件效果越好，但是付出的代价也越多。

与现有技术相比，本发明技术方案的有益效果是：采用本法能快速准确地实现地震波正演从而加速地下层反演成像，且对计算机的内存需要降低。

附图说明

图1为频率域波动方程正演流程图。

图2为本发明实施例中对待勘测区域进行Np=8的区域划分示意图。

图3为两重网格并行划分的x-y平面视图。

图4-7为四种不同震源位置示意图。

具体实施方式

附图仅用于示例性说明，不能理解为对本发明的限制。

下面结合附图和实施例对本发明的技术方案做进一步的说明。

一种基于MPI并行预条件迭代的地震成像正演方法，在待勘测区域中激发人工地震波，并由频率域波动方程模型-ΔU(x,y,z)-k²U(x,y,z)=G(x,y,z)来刻画，其中为拉普拉斯算子，U为待勘测区域的连续空间压力场，k为波数，G为连续的震源项；对频率域波动方程正演，求解出勘测区域内的压力场，模拟波传播，再根据压力场实现地震信息反演，获取地球的地下层图像；在本实施例中，如图1，频率域波动方程正演流程图包括：

步骤1：用基于频散极小化的27点有限差分方案来离散得到线性系统Au=g.其中系数矩阵A为稀疏的，并呈27对角的形状.右端项向量g是离散震源项，而向量u是要求的未知量，即离散压力场。由于针对的是实际工业应用，这里的数据全是三维的，并且规模巨大。这里要输入的数据为k，本发明中提供k=60和k=300两个情况。另一输入数据为G=δ(x-x₁,y-y₁,z-z₁)，即G在计算区域内的某点(x₁,y₁,z₁)为1其余为0，(x₁,y₁,z₁)可以在区域内任选。通过27点差分离散得到的系数矩阵A由于规模巨大、数据非常复杂，在这里不提供其具体数据。要输出的结果为离散压力场u，对于后面的计算结果，由于u的数据量非常大，也不一一显示其数据，只给出其实部的图像显示。

步骤2：对线性系统Au=g进行预条件处理，得到等价的线性系统：

AM^-1v=g，Mu=v，              （1）

其中，M是基于移位拉普拉斯算子离散得到的预条件子，和A具有同样的稀疏结构，其离散也是采用同样的频散极小化差分方法。

步骤3：对待勘测区域进行并行划分；具体为：

步骤31：为了建立并行的多重网格预条件Bi-CGStab法，首先对待勘测区域进行并行分划。由于实际待勘测区域为矩形，并行分划为在每个维度上进行2ⁿ分割，n>0为自然数。设x,y,z方向对应的剖分分别为p_x，p_y，p_z，则待勘测区域被分成Np=p_x·p_y·p_z部分，即数据被分划成Np个部分。一个进程处理一个部分的数据，这时有Np个进程可以同时处理数据。图2展示了p_x=p_y=p_z=2，即Np=8时的区域分划,8个子区域按照秩rank=0到rank=7排序。

步骤32：确定子区域的边界。令待勘测区域的计算剖分分别为Nx，Ny，Nz，即网格规模(未知量个数)为N=Nx·Ny·Nz；为表述方便这里Nx，Ny，Nz取为2的幂次方加1.在界定分化后的子区域边界时要考虑多重网格的并行实施。

以两重网格为例，粗网格的计算剖分分别为

$\mathrm{nz}=\frac{\mathrm{Nz}-1}{2},$ 设[·]表示取整，令 $\mathrm{lz}=\left[\frac{\mathrm{rank}}{\mathrm{px}\·\mathrm{py}}\right],$ $\mathrm{ly}=\left[\frac{\mathrm{rank}-\mathrm{lz}\·\mathrm{px}\·\mathrm{py}}{\mathrm{px}}\right],$ lx=[rank-ly·px-lz·px·py]，则对应秩为rank的子区域的细网格x,y,z方向边界取为：

$x:\mathrm{lx}\·\mathrm{nx}\→(\mathrm{lx}+1)\·\mathrm{nx}+2\·(1-[\frac{\mathrm{lx}+1}{\mathrm{px}}\left]\right),$

$y:\mathrm{ly}\·\mathrm{ny}\→(\mathrm{ly}+1)\·\mathrm{ny}+2\·(1-[\frac{\mathrm{ly}+1}{\mathrm{py}}\left]\right),$

$z:\mathrm{lz}\·\mathrm{nz}\→(\mathrm{lz}+1)\·\mathrm{nz}+2\·(1-[\frac{\mathrm{lz}+1}{\mathrm{pz}}\left]\right).$

对应秩为rank的子区域的粗网格x,y,z方向边界取为：

$x:\mathrm{lx}\·\mathrm{nx}\→(\mathrm{lx}+1)\·\mathrm{nx}+1-\left[\frac{\mathrm{lx}+1}{\mathrm{px}}\right],$

$y:\mathrm{ly}\·\mathrm{ny}\→(\mathrm{ly}+1)\·\mathrm{ny}+1-\left[\frac{\mathrm{ly}+1}{\mathrm{py}}\right],$

$z:\mathrm{lz}\·\mathrm{nz}\→(\mathrm{lz}+1)\·\mathrm{nz}+1-\left[\frac{\mathrm{lz}+1}{\mathrm{pz}}\right].$

这样，细网格子区域间的交叠部分有3个网格点，而粗网格子区域间的交叠部分为2个网格点。这样的分划可以保证多重网格的顺利实施。以p_x=p_y=p_z=2，Nx=Ny=Nz=17为例，对应rank=0的子区域细网格边界为x:0→6，y:0→6，z:0→6，而粗网格边界为x:0→3，y:0→3，z:0→3。图3展示了x-y平面对应rank=0，1，2，3的子区域边界。其中，数字0-8是细网格的标号，数字0’-4’是粗网格的标号，阴影部分为子区域的交叠部分。可见对应rank=0的子区域有49个细网格点和16个粗网格点。

步骤33：多重网格FMG的求解性能直接决定着Bi-CGStab的整体求解效果，至关重要。多重网格的重要组成元素有各层的网格算子A^h和网格间的转移算子，即细层到粗层的限制算子和粗层到细层的延拓算子。在进行多重网格预条件Bi-CGStab法求解时，首先要构造其组成元素。令A^h表示细网格算子，

表示细网格到粗网格的限制算子，

表示粗网格到细网格的延拓算子，则粗网格算子为

在这里采用全粗化的多重网格，即在设x,y,z三个方向都进行粗化。同时，利用全加权限制算子和黑盒子延拓算子，即延拓算子的构造是基于离散的稀疏矩阵。

步骤34：各子区域间的数据通讯。各子区域间虽有交叠的部分，但每个子区域都有人工边界，即不是物理区域的边界。边界上的值全为假值，这些假值会影响靠近人工边界处的内部值，从而造成错误的计算结果。为了结果的正确性;需要在各相邻的子区域间进行数据通信，互通有无。数据通信通过非阻塞通信MPI_Isend和MPI_Irecv来接收和发送。各自区域共享的数据通过MPI_Allreduce实现收集。

步骤4：在GPU集群上执行多重网格FMG预条件的Bi-CGStab迭代法，求方程在各子区域内的局部解，最后将这些局部解合并，得到频率域波动方程的整体解，即离散压力场u。

本实施例是采用广东省计算科学重点实验室的曙光多核计算集群（GPUCluster，峰值运算能力136T）上完成。通过以上方案计算的结果如图4-7所示，其中计算的波数k=60，离散的网格规模(未知量个数)为N=129×129×129.表格1给出了对应网格规模为N=513×513×513，约为1.35亿时的MPI并行计算加速比，此时波数k=300，实际计算区域为3776m×3776m×53776m.图4-7显示了计算结果的正确性，而表1显示了良好的并行加速比，表明本发明方案的优越性。

核数	计算时间(单位：秒)	加速比
			16	6487.8	---
32	3344.4	1.94
			64	1769.5	1.89
128	972.0	1.82
			256	562.2	1.73
512	355.8	1.58

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (7)

1.一种基于MPI并行预条件迭代的地震成像正演方法，其特征在于，在待勘测区域中激发人工地震波，通过频率域波动方程模型-ΔU(x,y,z)-k²U(x,y,z)=G(x,y,z)来模拟波传播；其中为拉普拉斯算子，U为待勘测区域的连续空间压力场，k为波数，G为连续的震源项；对频率域波动方程进行并行预条件正演求解，求解出勘测区域内的压力场，模拟地震波传播，再根据压力场实现地震信息反演，获取地球的地下层图像；所述对频率域波动方程进行并行预条件正演包括以下步骤：

S1.采用差分法对频率域波动方程进行离散，得到线性系统Au=g；其中系数矩阵A为稀疏的；右端项向量g是震源项的离散，未知向量u为压力场；

S2.对线性系统Au=g进行预条件处理，得到等价的线性系统；

AM^-1v=g，Mu=v，            （1）

其中M是基于移位拉普拉斯算子离散得到的预条件子，和A具有同样的稀疏结构；

S3.对待勘测区域进行并行划分；

S4.对步骤S3划分的子区域进行内外迭代，从而求解出线性系统（1）在各子区域内的局部解，最后将这些局部解合并，得到频率域波动方程的整体解，整体解即为离散的压力场u。

2.根据权利要求1所述的基于MPI并行预条件迭代的地震成像正演方法，其特征在于，所述步骤S1采用的差分法为基于频散极小化的有限差分法。

3.根据权利要求1所述的基于MPI并行预条件迭代的地震成像正演方法，其特征在于，所述步骤S2中M采用基于频散极小化差分法离散获得。

4.根据权利要求1、2或3所述的基于MPI并行预条件迭代的地震成像正演方法，其特征在于，所述步骤S3对待勘测区域进行并行网格划分的具体方式为：

S31.在三维空间的每个维度上进行2ⁿ分割，n>0为自然数，设三维空间中x,y,z方向对应的剖分分别为p_x，p_y，p_z，则待勘测区域被分成Np=p_x·p_y·p_z子区域，即待勘测区域中数据被分划成Np个子区域，所述Np个子区域采用Np个进程分别处理；

S32.确定子区域的边界，待勘测区域的计算剖分分别为Nx，Ny，Nz，即网格规模为N=Nx·Ny·Nz，其中Nx，Ny，Nz取为2的幂次方加1；在界定分化后的子区域边界时进行多重网格的并行实施；

S33.多重网格FMG的重要组成元素有各层的网格算子A^h和网格间的转移算子，即细层到粗层的限制算子和粗层到细层的延拓算子；在进行多重网格预条件Bi-CGStab法求解时，首先要构造其组成元素，令A^h表示细网格算子，

表示细网格到粗网格的限制算子，

表示粗网格到细网格的延拓算子，则粗网格算子为

采用全粗化的多重网格，即在设三维空间的x,y,z三个方向都进行粗化，并利用全加权限制算子和黑盒子延拓算子，即延拓算子的构造是基于离散的稀疏矩阵；

S34.各子区域间的数据通信，数据通信通过非阻塞通信MPI_Isend和MPI_Irecv来接收和发送，各子区域共享的数据通过MPI_Allreduce收集。

5.根据权利要求4所述的基于MPI并行预条件迭代的地震成像正演方法，其特征在于，所述步骤S4对格子区域进行内外迭代为：采用Krylov子空间方法作为外迭代来求解线性系统（1），并用瀑布型多重网格FMG作为内迭代来求解子系统Mu=v，作为子系统。

6.根据权利要求5所述的基于MPI并行预条件迭代的地震成像正演方法，其特征在于，所述采用双共轭梯度法Bi-CGStab作为外迭代来求解线性系统（1）。

7.根据权利要求5所述的基于MPI并行预条件迭代的地震成像正演方法，其特征在于，所述瀑布型多重网格FMG作为内迭代来求解子系统Mu=v时只迭代一次。

Images