0 引言

Helmholtz方程是一类非常重要的椭圆型偏微分方程，它可以用来描述波传播和散射现象，因此在电磁辐射、声学、光学以及地震学方面具有重要的研究价值.

为了在有界区域数值求解Helmholtz方程，我们引入PML(Bérenger, 1994; Tsynkov and Turkel, 2001;Chen et al., 2011, 2013).可以采用有限差分方法离散该方程：Jo等(1996)提出旋转九点差分格式，相比经典的五点差分格式，它在保持同样计算精度的前提下，能减少一个波长内的采样点数，从而减少计算量；由于Jo的方法不适用于横纵不等间距网格，Chen(2012)提出了平均导数法，该方法能适应于任意横纵网格间距; 与此同时，Chen等(2013)证明了Jo的方法不适用于PML，并且指出若采用平均导数法，差分方程在PML区域与微分方程逐点相容.

采用平均导数法离散该方程之后，我们得到了一个稀疏的线性系统，由于大波数的影响，这个线性系统是不定的，并且不具有厄密性.求解线性方程组的方法，一般可以分为两类：直接法(LU分解)和迭代法.二者各有优缺点：LU分解得到的结果可以适用于多炮并行计算，然而受到计算机内存和计算量的限制，LU分解不适用于求解大规模问题；基于Krylov子空间的迭代方法，例如Bi-CGSTAB方法和GMRES方法可以适用于大规模问题，但是如果线性系统没有好的谱分布，迭代方法的收敛速度很慢甚至停滞.因此我们需要一个好的预条件算子来改善线性系统的谱分布，许多学者为此做了很多工作：Bayliss等(1983)提出拉普拉斯预条件算子，结合迭代法CGNR高效地求解了Helmholtz方程；Laird等(2002)在拉普拉斯预条件算子的基础上增加了一个实数项，构造出来了实移位拉普拉斯预条件算子；Erlangga等(2004)总结了以上两种预条件算子，提出复移位拉普拉斯预条件算子，Van Gijzen等(2007)分析了此类预条件后的算子的谱分布，给出了如何选择最优的“移位”值，并且结合Krylov子空间迭代方法加速了Helmholtz方程的收敛；Erlangga等(2006b)比较了不完全LU分解预条件方法和复位移拉普拉斯算子预条件方法，并得到结论：用多重网格方法得到复移位拉普拉斯算子的近似逆后，结合Bi-CGSTAB迭代法，是一种快速且稳健的方法.

对于预条件的Krylov子空间方法，至关重要的一步就是如何求取预条件算子的逆.获得逆算子的显式形式是没有必要的，我们在尽量减少计算量的前提下，只需求取一个近似逆.因此多重网格方法是一个很好的选择.对于横纵网格等间距的情形，多重网格中又称为各向同性问题，反之称为各向异性问题(Briggs et al., 2000).

对于各向同性问题，Erlangga等(2006a, 2008)提出V(1, 1)循环的多重网格法去近似求解预条件算子的逆，这种多重网格方法改进了De Zeeuw(1990)提出的依赖差分算子的延拓算子，提出了一种适用于Helmholtz方程的延拓算子；Riyanti等(2007)提出将Erlangga的方法与Bi-CGSTAB方法结合，取得了很好的收敛速度；Riyanti等(2007)和Plessix(2007)将上述方法应用于三维Helmholtz方程求解，证明了迭代法对于大规模问题的有效性.Chen等(2011)针对Helmholtz方程给出了一种不同于De Zeeuw(1990)和Erlangga等(2006a)的全新的依赖于差分算子的延拓算子，取得了比经典多重网格方法收敛速度更快的效果.

对于各向异性问题，曹健等(2015)、Cao等(2018)指出采用经典的多重网格方法结合Bi-CGSTAB方法，并不能获得很好的收敛速度，甚至完全不收敛；本文分析了经典多重网格的三个重要组成成分：延拓算子、光滑算子、限制算子.进一步提出了适用于各向异性问题的多重网格算法：半粗化技术、线迭代技术以及依赖差分算子的延拓算子.本文的两个数值算例说明：经典多重网格方法不适用于各向异性问题；本文的算法对于中低频率的各向异性问题，收敛速度很快.

1 Helmholtz方程的离散以及谱分析 1.1 Helmholtz问题的数学定义

求解在二维声介质中波的传播问题，相应的Helmholtz方程为：

(1)

其中，u代表单频压力场，是拉普拉斯算子，是虚数单位.α是一个实数，表征介质的阻尼性质，通常0≤α < < 1，本文设定为0.k是波数，定义为k:=2πf/v(x)，f是频率，v(x)是介质速度.g(x, ω)是源项.

为了在有界区域数值求解Helmholtz问题，我们引入PML：

(2)

其中，是角频率，σ_x是仅与x坐标有关的函数，定义为：

(3)

σ_y也可以类似定义.其中f₀是雷克子波主频.如图 1所示：L_PML是PML区域的宽度.l_x是当前点距离∂Ω₀的距离，a₀是经验常数，一般取为a₀=1.79(Chen et al., 2011).

采用平均导数法离散Helmholtz方程(Chen, 2012)，从而得到线性系统：

(4)

上述线性系统，可以通过矩阵A表示；另外一种较为简洁的表示是差分模板：在(x_m, y_n)点，差分模板记为：

(5)

写出差分模板的各个元素：

(6)

其中，h_x, h_y分别是x和y方向的采样步长；α, β, c, d, e是通过最小化相速度误差得到的优化权重系数(Chen, 2012)；分别定义为：

(7)

其中：

(8)

从差分模板可以得到阻抗矩阵的一些性质，例如，若要求矩阵对称，需满足

(9)

若忽略PML，阻抗矩阵是对称的实矩阵，根据线性代数的理论，它的特征值全部是实数；由于PML的存在，阻抗矩阵是非对称复矩阵.

1.2 预条件

阻抗矩阵A具有A=L-z₁D的形式，z₁=1－αi(本文α取为0)，L和D分别是离散的Laplacian算子和零阶项e_xe_yk².由于PML的层数相比内部区域的层数小得多，所以A是一个复数域的、带状分布的、稀疏、近对称矩阵.

为了求解上述线性系统，我们采用Krylov子空间方法，特别是Bi-CGSTAB方法.但是该方法要求线性系统有较好的谱分布，从而保证问题快速收敛.对于Helmholtz方程而言，预条件算子可以是(Erlangga et al., 2006a;Chen et al., 2011)：

(10)

至此，我们得到了预条件的Helmholtz微分方程：

(11)

同样采用平均导数法在Ω上离散，从而得到预条件矩阵M以及预条件线性系统

(12)

其中.若(β₁, β₂)=(0, 0)，则它是拉普拉斯预条件算子(Bayliss et al., 1983)；Laird和Erlangga提到的预条件算子对应于.如果β₁>0, β₂>0，则它是复移位拉普拉斯算子.本文选取β₁=1, β₂=0.5.

1.3 谱分析

首先，忽略PML，考虑零值Dirichlet边界条件，形如的差分模板离散得到的系数矩阵的特征值为：

(13)

N_x, N_y分别为x, y方向的网格点数，l=1, 2，…，N_y, m=1, 2，…，N_x.对于本文的线性系统而言：

(14)

实际数值模拟中，我们需要考虑PML，此时，阻抗矩阵是近对称复矩阵，特征值没有解析表达式.

我们采用如下的试验计算两种情况下矩阵的谱分布：设定参数：N_x=N_y=70, 均匀介质v=3000 m·s^-1，f=30Hz(或者设定波数k=0.0628)，定义横纵采样间距比 .图 2所示为由(13)式计算得到的矩阵的谱分布(未考虑PML，取γ=2)，图 3所示为数值计算(考虑PML，每边层数设定为5，取γ=4)的矩阵的谱分布.

从图 2a可以看到，A是实对称矩阵，特征值全部是实数，而且特征值在(0, 0)多有聚集，条件数较大.从图 2b可以看到，用复移位拉普拉斯算子预条件后的阻抗矩阵AM^-1，聚集在(0, 0)处的特征值被“推”至复平面，从而改善了系统的特征值分布.可以证明(Chen et al., 2011)，如果μ是D^-1L的特征值，当且仅当是AM^-1的特征值，因此图 2b和图 2d等价；并且这个分式线性变换(共形变换)能将实轴变换成一个以为圆心，R= 为半径的圆，虚部小于0的点被映射至此圆的内部，反之被映射至此圆的外部.可以从图 2c和图 2d看出，AM^-1的特征值全部位于圆心C=(0.5, 0), 半径R=0.5的圆上.

从图 3a可以看到，由于PML的影响，A的特征值不再只有实数，而且矩阵的条件数仍然较大.图 3b是AM^-1的谱.图 3c是D^-1L的特征值分布，图 3d是将图 3c中的特征值作分式线性变换得到，它与图 3b等价.可以看出，由于D^-1L的所有特征值虚部全部小于0，因此经过分式线性变换之后，这些特征值全部位于圆心C=(0.5, 0), 半径R=0.5的圆的内部，如图 3d所示.换言之，我们证明了AM^-1的特征值全部位于圆心C=(0.5, 0), 半径R=0.5的圆的内部.可以看到，复移位拉普拉斯预条件算子对于改善Helmholtz方程的谱分布是很有效果的.

2 多重网格求解M^－1

为了行文方便，我们给出一些记号，如表 1所示.

2.1 多重网格的思想

我们可以通过研究最简单的两网格法来理解多重网格的思想，假设待求解的线性方程组为：

(15)

令表示某近似解，u_l表示精确解.则的修正值为

(16)

残差为

(17)

残差与修正值之间满足残量方程：

(18)

现在需要用M_l的一个近似来求解e_l.

经典迭代方法，例如Jacobi迭代，Gauss-Seidel迭代是通过寻找下述方程的一个近似解来完成的：

(19)

是一个"较简单"的算子. 例如, 在Jacobi迭代中, 是M_l的对角线部分; 在Gauss-Seidel迭代中, 是M_l的下三角部分.第二个近似是根据(20)式得到的:

(20)

与经典迭代方法不同，多重网格是寻找“粗糙化”的算子M_l－1而不是“简单化”的算子.在网格大小为h_{l－1, x}, h_{l－1, y}的一个更粗糙的网格上(常取h_{l－1, x}=2h_{l, x}, h_{l－1, y}=2h_{l, y}，即经典的全粗化方法；也有特殊的取法，如本文采用的半粗化：h_{l－1, x}=h_{l, x}, h_{l－1, y}=2h_{l, y})，构造M_l的某个适当近似算子M_l－1.则残量方程近似为：

(21)

由于M_l－1的维数较小，该方程比原方程更易求解.为了定义粗网格上的亏损r_l－1，需要一个限制算子R_l－1，将细网格上的亏损r_l限制在粗网格上

(22)

一旦得到(22)式的一个解后，还需要一个延拓算子P_l将修正值内插到细网格Ω_l上

(23)

最后u_l近似解可以校正为

(24)

总而言之，多重网格的逻辑是：

(1) 误差可划分为高频摆动分量和低频光滑分量；

(2) 通过简单的定常迭代去消除高频摆动分量；

(3) 在粗网格上消除那些顽固的低频光滑分量.

2.2 各向异性问题中的光滑、限制以及延拓算子

最流行，而且也是最应该试验的光滑方法是Gauss-Seidel方法，因为它能产生很好的收敛速度.然而Gauss-Seidel方法的精确形式依赖于网格点所选取的排序.对于典型的等间距差分格式，最好使用红-黑排序；当沿某一维比沿另一维的耦合程度更强时，例如本文γ=2, 4的情形，沿y维比沿x维的耦合程度更强，我们可以考虑沿着y维同时将整条线松弛：

(25)

我们称之为y-line Gauss-Seidel.最邻近耦合的线松弛涉及一个三对角系统的求解，因而仍然具有很高的效率；如果选择红黑斑马线迭代，则根据m的奇偶性，对网格点涂色分类：当m为奇数，则为红色；反之为黑色.迭代可以分为两步，首先对红网格点更新：

(26)

再对黑网格点更新：

(27)

红与黑两部分迭代是分开进行的，因此适合并行.

延拓和限制算子的一个简洁的表示方法是，给出它们的符号.选择粗网格Ω_k－1上的脉冲函数u_k－1=δ_k－1∈U_k－1, 它在某个(x, y)∈Ω_k－1上取1，其他点取0.再定义Pδ_k－1，从而得到符号P.最流行的拓展算子是简单的双线性插值算子，这种算子适用于各向同性问题，它给出Ω_k上九个点(x, y), (x+h_{k, x}, y), …, (x－h_{k, x}, y－h_{k, y})上的非零值，分别为1，1/2，…，1/4.因此符号为：

对于本文的各向异性问题，我们选择的半粗化延拓算子的符号为：

关于限制算子的选择，我们通常可以选为延拓算子的伴随算子.根据算子的符号记法，限制算子和延拓算子的符号元素之间的关系为：

(28)

因此，完全加权限制算子，就是双线性插值算子的伴随.对于非均匀不连续介质问题，我们可以选择依赖差分算子的延拓算子(De Zeeuw, 1990；Erlangga et al., 2006a).本文采用Chen(2012)提出的适用于Helmholtz的延拓算子.

2.3 本文的算法框架

针对γ=2的情形，我们给出多重网格V循环求解M_lu_l=f_l的递归算法, 如表 2所示.

这里，粗网格Ω_k－1上的阻抗矩阵M_k－1，通常可以通过两种方式获得：

(ⅰ)类似于M_l，可以通过平均导数法，在Ω_k－1上离散算子得到.

(ⅱ) Galerkin方法：M_k－1=R_k－1M_kP_k.

本文选用第二种方法来获取M_k－1，它的好处在于，我们无须考虑边界条件.

3 数值算例

本节选取了两个模型，来说明本文算法的有效性.我们采用V(1, 1)多重网格，层数设定为9层.PML层数设定为每一个维度的网格点数的1%，震源设置在PML吸收层以下第二层中间的位置.收敛的准则是相对残差的二范数达到10^-6，如果迭代次数超过10000次仍未达到此准则，我们认定它不收敛.所有测试均在Intel 8核，内存16 GB的机器上，由Fortran 77程序完成.

3.1 均匀模型

首先，考虑一个速度大小为3000 m·s^-1的均匀模型，横纵采样点数均为401；为了测试本文算法对于各向异性问题的有效性，我们选定γ=1, 2, 4，从低到高选取了6个频率，分别是f=10, 20, 30, 40, 50, 60 Hz.本文的方法我们记为Ⅰ，经典的多重网格算法我们记为Ⅱ，测试结果如表 3和图 4所示，表中数据为Bi-CGSTAB的迭代次数以及CPU耗时.

从测试结果可以看出，对于各向同性问题，本文算法和经典多重网格算法一致，收敛速度较快；对于弱各向异性问题(γ=2)，经典多重网格仅在低频段缓慢收敛，中频段不收敛，而本文算法仍然能快速收敛；对于强各向异性问题(γ=4)，经典多重网格已经失效，而本文算法在中低频段仍然能有令人满意的收敛速度.

3.2 部分Marmousi模型

本节选用部分Marmousi模型，来说明本文算法对于强不均匀模型的有效性；我们考虑的模型横纵采样点数为301×301，速度模型如图 5所示，选定γ=1, 2, 4，测试了两个频率f=10, 30 Hz，结果如图 6所示.

经典的多重网格方法对于本例已经失效，从测试结果可以看到，本文的算法对于强不均匀介质的中低频段，收敛速度仍然很快；随着频率的提高，由于零特征值的影响，收敛速度变慢.

3.3 内存消耗和耗时对比

本小节测试了本文方法的内存需求以及计算速度，采用的模型是3.1节的均匀模型和3.2节的部分Marmousi模型，并将计算结果与直接分解方法作了对比，结果如表 4所示.

表中，M_LU和T_LU分别代表LU分解法的内存消耗以及CPU耗时，M_iter和T_iter分别代表本文的迭代法内存消耗以及耗时.迭代法内存消耗远小于直接法，并且在弱各向异性的情况下，中低频率的迭代问题，迭代法相比直接法也明显更省时.

4 结论与建议

Helmholtz方程的高效数值求解，不仅依赖于离散方法，更与算法的选择密切相关.本文的研究重点是：横纵不等间距采样情况下，Helmholtz方程的高效迭代解法.本文采用平均导数法离散Helmholtz方程；进一步，我们分析了原有线性系统的谱，发现如果不做预条件处理，条件数太大，特征值在0处多有聚集，导致经典的迭代方法失效.复移位拉普拉斯算子是一种高效的针对Helmholtz方程的预条件算子，根据谱分析结果，我们知道，预条件后的线性系统的特征值将全部落在圆心C=(0.5, 0)，半径R=0.5的圆内，极大地改善了系统的谱分布，从而使得Krylov子空间方法奏效.

对于预条件算子的近似逆，我们通过多重网格的方法获得.对于横纵不等间距采样的情况，多重网格中称为各向异性问题，从我们的数值试验可以看出，经典的多重网格方法对于此类问题已经失效.本文提出的方法对于中低频率的各向异性问题的求解是高效的；并且与直接法的内存消耗与耗时对比也证明了本文方法相比直接法的优势.

对于Helmholtz方程的高频迭代问题，仍然是一个难题，可能结合deflation技术(Nicolaides, 1987;Morgan, 1995; Frank and Vuik, 2001)，消除小特征值对于迭代的影响，能加速高频Helmholtz问题的收敛.

致谢  感谢审稿专家和编辑提出的宝贵意见