一种利用Nyström离散与FFT快速褶积的散射地震波并行计算方法

徐杨杨， 孙建国 ， 商耀达

吉林大学地球探测科学与技术学院， 长春 130026

0 引言

类似于光学散射和电磁散射问题的数值模拟正演方法，常规的声波散射数值模拟方法也分为两大类：基于微分方程理论的有限差分法(Finite Difference)、有限元法(Finite Element)等和基于积分方程理论的边界元法、矩量法、Nyström法等(Fu et al.,1997;Fu and Wu,2000; Hsiao et al.,2000；刘立彬等，2020；邱长凯等，2020；徐杨杨等，2021).微分方程类的数值算法具有容易实现、计算速度快等优点，但是网格剖分需要在全空间上进行，对于较小的异常体也会生成规模非常大的离散矩阵，并且容易受到数值频散的影响.而积分方程类的数值算法具有半解析性质、只在散射异常体上进行网格离散即可.因此，对小型散射异常体的散射波场数值模拟，积分方程法生成的离散矩阵规模小，降低了计算和存储的要求，其计算速度和精度要明显优于基于微分方程的数值算法.对地震散射问题中的大尺度复杂地质条件而言，由于积分方程法形成的线性方程组系数矩阵的密实性，并且离散矩阵元素的计算过程复杂，计算耗时、存储量大，相比于计算电磁学，基于积分方程的数值方法在地球物理中的应用受到了非常大的限制(王德智，2015).

针对积分方程法处理散射问题时所存在的问题，国内外学者在积分方程快速算法方面作了大量研究工作.由于求解积分方程数值法生成的大型线性方程组非常耗时，有学者提出采用近似解法来避免解方程组.Zhdanov和Fang(1996)、Zhdanov等(2000)提出了一种拟线性的近似解法，认为散射异常场是背景场与反射系数的乘积，反射系数是最小二乘问题的解，该方法在电磁散射正反演中取得较好的应用效果.Sun(2005)根据计算电磁学的拟解析线性近似技术，将电反射率张量的概念引入直流电场积分方程数值法中，给出了拟线性级数的表达形式，提供了一种新的直流电场的正反演方法.孙建国(2006)借鉴电磁散射问题中的拟线性近似技术，通过假设散射波场与背景波场之间的线性关系，建立了声波散射的拟解析近似表达式.拟解析近似是半解析解，只需计算数值积分，该方法可以有效的避开了复杂介质情况下大型线性方程组的求解.

殷长春等(2018)利用插值技术将准线性近似算法应用于三维情况下航空电磁数据的数值模拟计算，在稀疏网格上计算散射场和点反射系数，然后进行线性插值，最后根据拟线性近似技术来求解电磁异常体的散射波场.另一类方法利用格林函数的Toeplitz性质，通过快速傅氏变换(Fast Fourier Transform，简称FFT)来加速线性方程组的迭代类求解方法中的矩阵向量乘积，称为积分方程快速傅里叶变换法(Integral Equation Fast Fourier Transform，简称IE-FFT)(Peters and Volakis,1988; Volakis and Barkeshli,1991)，达到降低存储和计算复杂的目的，如共轭梯度快速傅里叶变换(Conjugate Gradient Fast Fourier Transform，简称CG-FFT).Osnabrugge 等(2016)在背景波数中引入虚部分量，得到了带阻尼因子的背景格林函数.通过加入预条件因子进一步加速了迭代的收敛性，同时利用FFT实现了循环卷积的快速算法.该方法在光学散射数值模拟中表现出了精度高、计算效率快的优点.

此外，传统的求解积分方程的矩量法计算复杂，离散后的线性方程组系数矩阵每个元素都需要计算内积积分，计算量大，效率低.虽然快速多极子方法(Fast Multiple Method，简称FMM)和多层快速多极子方法(Multievel Fast Multiple Algorithm，简称MLFMA)(Rokhlin,1990; Çakir,2006; Çakr,2009；崔晓娜和刘洪，2020)能够有效地降低存储和计算复杂度，但对积分核依赖性强，且实现过程较为复杂.

除了上述方法之外，在计算数学中还有一种解积分方程的方法，称为Nyström法.该方法利用适当的数值求积公式来近似积分算子(Nystöm离散)，可以避免生成系数矩阵元素时的大量积分计算，相较于传统的矩量法计算成本较低.如果积分方程的核函数只由Green函数组成，则基于数值积分公式的Nyström法所生成系数矩阵只对积分核进行简单计算并且具有Toeplitz性质，计算复杂度低且易实现并行化.对于奇异性的处理，可以采用阶数较高的插值基函数、奇异性减法等技术均可处理奇异性(Canino et al.，1998).

然而，在地震散射问题中出现的L-S方程的积分核除了Green函数之外还包含有速度扰动项.因此，为了使得积分方程的核函数只包含Green函数，首先将L-S方程改写为一个与原积分方程等价的方程(等价L-S方程)，并对该方程进行逐点归一化，进而得到在离散后其系数矩阵具有Toeplitz性质的归一化等价L-S方程.然后，对归一化等价L-S方程进行Nyström离散，从而得到其系数矩阵具有Toeplitz性质的线性方程组，极大地降低了存储量.同时，为了提高空间褶积的计算效率，将归一化等价L-S方程变换到波数域，在其傅里叶谱空间中实现矩阵和向量之间的乘积，进而实现对线性方程组的快速求解.通过简单的分析可知，上述两个措施可将内存和计算复杂度由O(N×N)分别降为O(N)和O(NlogN)(N=NxNz，Nx、Nz分别为x、z方向上的网格剖分数量).此外，因数值模拟在频率域进行，各频率下的波场值计算是独立的.因此，为进一步提高计算效率，设计了基于MPI+OpenMP的进程级和线程级混合并行的实现模式，在有效降低其内存消耗的基础上最大限度地提高积分方程法散射波场数值模拟的计算效率.

1 等价与归一化等价L-S方程

利用扰动理论，可将描述散射问题的Helmholtz方程转化成为Lippmann-Schwinger(L-S)积分方程：

U(r)=Ub(r)+Us(r)

(1)

由于格林函数G(r,r′)中的Hankel函数在r=r′点包含In(k0|r-r′|)项而奇异，考虑Cauchy主值意义下的收敛性，(1)式右端的积分项在点r∈Ω上一致连续(Fu et al.,1997)，即L-S方程具有弱奇异性.应用Fredholm择一定理，可知L-S方程存在唯一解.

(2)

(3)

这是一个积分核为Green函数的积分方程.因此，如果首先对方程(2)进行某种方式的数值离散，然后再进行形如方程(3)的逐点归一化处理，则可得到其系数矩阵具有Toeplitz性质的线性方程组.

2 等价L-S方程的Nyström离散

令Pn是插值算子，{rj,j=1,2,…,n}⊂Ω是插值基点，{lj(r′),j=j=1,2,…,n}是插值基函数，满足lj(r′i)=δij.于是：

(4)

定义弱奇异积分算子：

(5)

则其近似算子满足：

(6)

其中权函数：

(7)

(8)

(9)

定义常元插值算子：

(10)

其中，Mj为剖分单元Vj的中心.此时，需要计算弱奇异积分：

(11)

从而得到权系数.由于当r∈Vj且r=r′时G(r,r′)奇异，所以必须要对(11)进行奇异性消除处理，具体见刘宁和孙建国(2007)或张金会和孙建国(2009) .

3 线性方程组的逐次超松弛迭代

Lippmann-Schwinger积分方程的常元Nyström近似生成的系数矩阵为稠密矩阵，并且每一个矩阵元都必须计算格林函数在剖分元上的弱奇异积分.当散射体规模较大时，网格剖分数量也很大，在二维情况下，系数矩阵由N2个元素组成(N=NxNz，Nx、Nz分别为x、z方向上的网格剖分数量).采用高斯消元法或LU分解等直接法求解时需要存储整个矩阵及其逆矩阵，内存需求为O(N2)并且矩阵求逆计算复杂度为O(N3).而迭代法无需存储逆矩阵，其复杂度下降为O(N2).常用的迭代算法有Gauss-Seisel迭代、超松弛迭代(SOR)、共轭梯度法、广义最小残差法以及其他Krylov子空间法(Kleinman et al., 1990a,b; Kleinman and Van Den Berg,1991; Yu and Fu, 2014).本文采用逐次超松弛迭代法(Kleinman and Van Den Berg, 1991)来求解L-S方程离散生成的线性系统.

(12)

(13)

式中，βn为复逐次超松弛因子.根据Kleinman和Van Den Berg(1991)的超松弛迭代在光学散射中的应用经验，通过在迭代过程中最小化残差‖rn‖ 来得到逐次复松弛因子βn，即取βn使得 ‖rn‖ 在Hilbert空间中取极小：

(14)

4 矩阵向量乘积的FFT实现

(15a)

(15b)

其中：

(15c)

对于原L-S积分方程，核函数为α(r′)G(r,r′).由于速度扰动函数在每个面元上的值不同，导致带速度扰动权的积分核经数值离散后的系数矩阵不具有Toeplitz性质.因此，有必要对L-S方程进行改写，使得积分核仅包含Green函数，保证离散后的系数矩阵具有Toeplitz性质.

任给一n阶循环矩阵C=circ(c0,c1,…,cn)可以被傅里叶矩阵F对角化(Davis，2013)，即:

(16)

(17)

对于一个n阶对称Toeolitz矩阵T，可以通过将它嵌到2n阶循环矩阵C中，来计算矩阵向量乘积，即:

(18)

从而有:

(19)

根据式(19)可以得到矩阵向量乘积Tx.

(20)

其中:

0≤i≤Nx-1，

(21)

5 频率域散射波场正演的并行化

MPI并行标准库通过消息传递来实现程序的并行化，是一种进程级的并行工具.它具有效率高、可移植性强等优势，是目前实际并行计算中普遍采用的并行编程环境.OpenMP作为一种在内存共享基础上实现多个线程间的并行计算编程接口，通过简单的高级语言指令即可实现多线程并行化计算(Quinn，2004).Nyström法求解Lippmann-Schwinger积分方程是在频率域中实现的，各频率片之间是独立计算的，可以通过MPI并行处理技术实现频率波场计算并行化.此外，当地下散射异常体规模较大时，剖分网格数量将会很大，求解单频率下的线性方程组时通常含有大量稠密矩阵向量乘法运算相关的for循环，在单个节点上使用OpenMP实现for循环的并行化可进一步提高效率.基于MPI+OpenMP的混合并行模式可以减少MPI的通信开销，充分利用资源，得到更高的加速比，在有效降低其内存消耗的基础上最大限度地提高积分方程法散射波场数值模拟的计算效率.

5.1 基于MPI的频率并行

基于MPI并行实现频率片并行化的算法设计过程中，为获得单炮所有频率下的波场值，可以采用主从式并行模式，主进程进行消息发送与接收以及其他输入输出任务，从进程负责所有频率的波场值计算.具体并行实现为主进程进行各从进程频率数和偏移量的计算并将消息发送给从进程，当从进程计算完相应频率数的波场值后，主进程接收各从进程发送的相应频率数、偏移量和波场值，最终得到所有接收道和所有频率波场值生成的二维数组.基于MPI并行的散射波场数值模拟实现步骤如下：

(1)调用MPI_Init函数,并生成通信域，参数初始化；

(2)根据集群中的节点数量，调用MPI_Comm_size函数来确定MPI_COMM_WORLD中的总进程数量(一个节点执行一个进程)，并调用MPI_Comm_rank确定进程编号；

(3)主进程根据设定的总进程数计算各从进程需要计算的频率数和频率偏移量并发送给相应的从进程；

(4)各从进程接收所需计算的频率数和频率偏移量，并计算相应频率数下的所有道的波场值，并将计算结果发送回主进程；

(5)主进程接收从进程发送的频率数、频率偏移量和计算结果，对所有频率下的波场值进行逆傅氏变换得到时域单炮记录；

(6)调用MPI_Finalize，结束并行环境.

5.2 基于OpenMP的节点内for循环并行

L-S方程的Nyström近似迭代解法实现过程中涉及以for循环形式出现的稠密矩阵向量乘法运算，而线程级并行的OpenMP并不需要进行复杂的线程创建、同步、销毁等工作.基于OpenMP的节点内循环并行化实现步骤如下：

(1)按照单个计算节点上的CPU核的数目，调用omp_set_num_threads设定线程数；

(2)在无数据相关性的for循环外部调用#pragma omp parallel for语句，将该循环并行化；

(3)调用private(〈variable list〉)和shared(〈variable list〉)设置并行域内的私有变量和共享变量.

6 数值试验

下面通过数值试验来测试基于Nyström法的IE-FFT并行算法在解决地震散射波场正演模拟中的有效性以及计算效率.首先设计了球形单体散射模型，通过与直接解法得到的数值结果进行对比，分析IE-FFT算法在存储和计算效率上的优势.然后用更为复杂的SEG/EAGE的盐丘模型对快速算法进行进一步测试.数值试验环境：Linux操作系统，包含11个节点，每个节点2个CPU，每个CPU包含8个核(CPU型号：Intel(R) Xeon(R) CPU E5-2620 0 @ 2.00 GHz )，使用基于MPI+OpenMP的并行算法对散射模型进行单炮数值模拟.

6.1 球体散射模型

球状散射体速度模型为均匀的背景介质中含有一半径200 m的高速球状散射体(如图1).模型大小为3000 m×3000 m，球散射体速度为3000 m×3000 m，背景速度2000 m·s-1，震源为主频15 Hz的雷克子波，震源位置rs=(1500 m,10 m)，检波器置于z=10 m的界面上，空间采样间隔为dx=dz=5 m，时间采样间隔为4 ms，记录时间2.4 s.

图1 球体散射模型Fig.1 Sphere scattering model

基于MPI+OpenMP，开辟11个进程进行基于MPI的频率并行计算，每个进程开6个线程进行基于OpenMP的for循环并行计算，分别采用直接法和IE-FFT迭代算法计算球体散射模型正演结果.对于单球体散射速度模型，仅需要在散射异常体上进行离散.两种方法得到的单频波场(f=15 Hz)和单炮记录结果如图2和图3，IE-FFT迭代快速算法与直接法求解得出的数值结果具有较好的一致性.为了更细致的分析，分别从两种算法得到的单炮记录中随机抽取2道进行对比(如图4)，结果显示出，IE-FFT算法得到的数值结果在同相轴能量上与直接求解法结果有极小差别，量级在误差允许范围之内，两种方法同相轴相位对应一致.

图2 球体散射模型直接法与IE-FFT迭代法15Hz单频波场数值结果(a) 直接法单频波场实部； (b) 直接法单频波场虚部； (c) IE-FFT迭代法单频波场实部； (d) IE-FFT迭代法单频波场虚部.Fig.2 Numerical results of 15Hz single-frequency wave field using the direct method and IE-FFT iterative method on the sphere scattering model(a) The real part of the single-frequency wave field using direct method； (b) The imaginary part of the single-frequency wave field using direct method； (c) The real part of the single-frequency wave field using the IE-FFT iterative method； (d) The imaginary part of the single-frequency wave field using the IE-FFT iterative method.

图3 球体散射模型直接法(a) 与IE-FFT迭代法; (b) 数值结果.Fig.3 Numerical results using the direct method (a) and IE-FFT iterative method (b) on the sphere scattering model

图4 球体散射模型直接法(实线)与IE-FFT迭代法(点划线)单道对比(a) 第40道单道对比； (b) 第150道弹道对比.Fig.4 Comparison between the direct method (solid line) and the IE-FFT iterative method (dotted line) for the sphere scattering model(a) Comparison of trace 40； (b) Comparison of trace 150.

球体散射模型直接法与IE-FFT迭代法的CPU时间与内存消耗对比见表1.数值实现过程中，直接法耗费CPU时间为4624.27 s，每个进程内存消耗为515 M，IE-FFT法迭代50次耗费CPU计算时间为396.37 s，每个进程内存消耗为78 M.IE-FFT迭代算法的数值结果可以较好地吻合直接法结果，而计算速度是直接法的11倍，并且应用格林函数的Toeplitz性质降低了系数矩阵的内存需求，内存消耗仅约为直接法的1/7.

6.2 SEG/EAGE盐丘模型

当地下介质中包含较大尺度的复杂散射体且速度扰动更强时，积分方程法应用于此类介质模型时需要在全空间中进行离散计算.而直接法因其离散后的矩阵为NxNz阶方阵，普通的微机难以获取相应的栈内存来存储系数矩阵且计算效率不满足实际应用需求.为了验证IE-FFT迭代算法在复杂强散射介质中的适应性，选取SEG/GAGE盐丘速度模型.模型参数如图5所示，模型大小为7800 m×2090 m，背景速度1500 m·s-1，震源子波为主频15 Hz的雷克子波，源位置rs=(3800 m,10 m)，检波器置于z=10 m的界面上，迭代法空间采样间隔为dx=dz=10 m，时间采样间隔为4 ms，记录时间4 s.

图5 SEG/EAGE盐丘模型Fig.5 SEG/EAGE salt dome model

IE-FFT快速迭代算法得到的单频波场(f=15 Hz)和单炮地震记录如图6和图7所示.从单炮地震记录上可以清晰分辨盐丘顶端以及岩体右侧横向起伏变化剧烈处的散射同相轴，对强散射体的识别能力较强，并且模型底部测试层的同相轴清晰、能量强，并未受到高速异常体的屏蔽影响.在直接法无法求解大型线性方程组的情况下，IE-FFT快速迭代法迭代100次耗费CPU计算时间和每个进程内存消耗为为9466.54 s和321 M(见表1).

图6 IE-FFT迭代算法15HZ单频波场数值结果(a) 单频波场实部； (b) 单频波场虚部.Fig.6 Numerical results of 15Hz single-frequency wave field using the IE-FFT iterative method on the SEG/EAGE salt dome model(a) The real part of the single-frequency wave field； (b) The imaginary part of the single-frequency wave field.

图7 IE-FFT迭代算法散射波场数值结果Fig.7 Numerical result using the IE-FFT iterative method on the SEG/EAGE salt dome model

表1 直接法与IE-FFT迭代法CPU时间与内存消耗Table 1 CPU time and memory consumption of the direct method and the IE-FFT iterative method

7 结论

传统上，对L-S方程的数值求解一般采用矩量法.此时，相应线性方程组的系数矩阵元素都是以内积积分的形式出现的.虽然可以利用数值积分法计算这些内积积分且可以获得较高的精度，但是其计算过程复杂、对计算机资源要求高且效率低.与此相反，基于数值积分公式的Nyström离散可使系数矩阵只包含对积分核(Green函数)的积分计算.因此，从计算效率方面考虑，对于大尺度地质条件下的地震散射波场计算问题，Nyström法会更加适合.事实上，通过对L-S方程的变形处理，可使得最终形成的线性方程组的系数矩阵具有Toeplitz性质，进而可将存储量由对N阶矩阵直接进行存储的O(N×N)降为O(N).此外，利用FFT和IFFT计算离散空间褶积，可以加速线性方程组求解迭代计算过程中的矩矢乘积，使其计算复杂度由O(N×N)降为O(NlogN).再有，根据各单频波场在计算上相互独立的特点，采用了基于MPI+OpenMP的进程级和线程级相结合的并行化计算方案.数值结果表明，相较于传统的积分方程数值算法，利用Nyström离散改写后的等价L-S方程且同时采用FFT计算离散褶积和并行计算的数值求解方案可以得到与直接求解法相当的数值结果，在保证计算精度的同时，内存和CPU耗时则显著降低，有效地克服了传统数值算法在求解L-S方程时系数矩阵存储和线性方程组求解效率低的问题.