CPML在瞬变电磁法时域有限差分三维正演中的应用分析

目前瞬变电磁法常用的数值模拟方法主要有积分方程法(矩量法、边界单元法)、微分方程法(有限元法、有限差分法、多分辨分析法)和混合方法(有限体积法).上述方法中频域方法较多见，在处理地空边界采用时多采用延拓方法，地下边界采用Dirichlet边界条件(Mur, 1981；邓晓红，2001；邢涛等，2008；赵云威，2012)，这要求边界必须距离场源足够远，才能使得场量在地下边界处近似为零，导致计算空间过大，计算耗时过长.为了减少计算量，一些学者(徐凯军和李桐林, 2004；孙怀凤等，2013)在低频条件下，忽略感应电流的影响，将Maxwell波动方程简化为扩散方程，从而可以扩大时间步长，减少迭代次数和计算量.但是这种忽略电磁波早期响应的方法在背景电阻率较高(如干旱地区)时加大了探测盲区，并且模拟高阻目标体时容易导致计算发散.还有一些学者(熊彬和罗延钟，2006；杨海燕和岳建华，2008)在建模方面采用非均匀网格、粗细网格剖分方法，对重点研究区的网格剖分较细，非关心区域的网格剖分较大，可以有效减少网格数量，其缺点是模型形态复杂时三维网格剖分很困难，适合规则目标形态的地质体建模.

在计算电磁学中，对于电磁场传播的开域问题，通常采用的方法是对空间进行截断，只计算重点区域，在截断边界上场量通常不满足Dirichlet边界条件，但是可以采用合适的吸收边界，使电磁波穿过吸收边界时无反射并且迅速衰减到零，在分辨率不变的情况下可以大大减少计算空间.

吸收边界的计算精度决定了时域有限差分法的结果精度，直接限制了有限差分的发展.1981年Mur提出了吸收边界条件，Mur吸收边界条件易于实现，占用资源少，但吸收效率与电磁波的入射角度有关，在顶角区域数值反射严重.Liao等(1984)提出了廖氏吸收边界条件，吸收效率优于Mur，但它的高阶吸收边界条件不稳定.Berenge(1994, 1997)提出了完善匹配层(Perfectly Matched Layer，PML)边界条件，PML是一种由虚构的本构参量组成的特殊媒质，通过匹配条件实现电磁波无反射向外传播，与频率和入射角度无关，但BPML吸收凋落波效率较低，网格剖分变大或者迭代时间较长时，会有迭代后期反射.随后Sacks等(1995)、Gedney(1996)提出各向异性PML(UPML)吸收边界条件，使用各向异性材料而不改变电磁场形式，但是电磁波的低频成份吸收特性差.Kuzuoglu和Mittra(1996)和Mittra(2000)提出了复频率参数PML(CFS-PML)，在吸收凋落波和长时间时域计算中具有高效率，但是边界条件的卷积需要保存以往所有时间步的场值，存储空间消耗极大，不易实现.Roden和Gedney(2000)将CFS-PML的卷积项离散为递归算法，从而可以高效地执行CFS-PML，称为卷积PML(Convolution PML CPML)，作者提出的参数分布对高频波的吸收很好，但对低频的吸收效率缺乏理论支持，数值模拟表明，较小的复频移参数变化对瞬变电磁法长周期扩散的电磁波计算精度影响很大，本文对CPML参数对瞬变电磁法常见介质的吸收系数和反射系数进行了分析和推导，表明合适的参数分布可以在保证计算精度的情况下，有效缩小计算空间.

在无磁性介质中，处于坐标伸缩中的CPML频域方程为(Roden and Gedney, 2000)：

式中：Sew =κew +σpew /(αew+jωε0)，Smw=κmw+σpmw/(αew+jωμ0)，(w=x, y, z)，为坐标伸缩因子，κew、σpew、αew、κmw、σpmw、αmw为新引入的参数，σw、εw为背景介质的电导率和介电常数，σpew、σpmw为C PML区域人为添加的电导率和磁导率，加入下标p是为了区分背景介质电导率.在计算域中，取Sew=Smw=1，则(1) 式和(2) 式变为正常的Maxwell方程.在计算域与PML层的内边界以及PML各层边界上，电磁波零反射的条件为：Sew=Smw；边界两侧的介质本构参数相同；电磁波传播方向上的横向伸缩因子相同(葛德彪和闫玉波，2011).将(1) 式和(2) 式的频域方程变换到时域，并且采用YEE均匀剖分网格进行差分离散，电场分量Ex的更新方程为

其中卷积递推项ϕexy定义在Ex的y方向CPML区域，ϕexz定义在Ex的z方向CPML区域，在不同的方向只需要计算其中一个卷积递推项即可，比如在计算Ex沿y方向的传播时只计算ϕexy，电磁场其他分量及其卷积递推项类似.

目前CPML吸收边界条件在雷达散射截面、天线以及微波器件等领域应用较多，对于比例伸缩因子，Roden和Gedney(2000)提出参数分布和取值范围适合吸收高频电磁波，不适用于瞬变电磁法的低频电磁波.通过对空气介质、似稳态条件下的大地介质的吸收边界参数进行理论推导，分析了衰减系数和反射系数对误差的影响，改进了吸收边界参数的分布，适合低频电磁波的吸收，取得了很好的吸收效果.

数为了简化CPML引入的参数对吸收效率的影响，考虑各向同性介质(本构参数ε=ε0εr，σ)中，二维情况下的任意方向传播的TEz极化的平面波，场量仅存Ex、Ey和Hz，见图 1.

场量的频域表达式为

代入(1) 式(2) 式，可解出A、B为

式中

将A、B代回(7)~(9) 中，归一化场量可以表示为

式中Re指实部，Im为虚部，伸缩因子Sei、Smi满足零反射条件时G=1.第一项表示平面波，第二项控制电磁场沿x方向的衰减，第三项控制y方向的衰减，可见CPML的吸收效率与频率有关.λ=ωIm(A)可以统一表示电磁波在整个区域(计算域和CPML域)的x方向的衰减系数，λ绝对值越大，电磁波的衰减越强.由于任意方向入射的平面波都可以分解为TM波和TE波，在垂直传播方向(θ=0) 的场量不衰减，只有平行传播方向(θ=90) 的场量才有衰减，因此CPML吸收效率与入射角度无关.公式为

在空气介质(σ=0) 的CPML区域，衰减系数为

取α=0则有最大衰减，衰减与频率无关，高频和低频的衰减系数都是Z0σp(Z0为自由空间波数)，d层CPML的反射系数r为

Δi为CPML的网格大小.由于电导率σp在CPML的内部边界从零值开始，逐步向外部边界取得最大值，因此α在内部边界不能为零，为了能够有效吸收低频波，离开此边界后，α应该降低为零.

在大地介质中，满足似稳态条件下(σ≫ωε0)，有：

衰减系数可表示为

假定α=rωε0，r为待定系数，式(16) 可以表示为

在大地电导率σ一定的情况下，对各个频率而言，κ、σp越大则吸收效果越强，但是由于存在反射，κ、σp值增加会使伸缩因子Sei和Smi的差异变大，反射率也会增加，取值并不是越大越好.在分析κ、σp的分布之前，首先要考虑α参数的影响.为了有效吸收低频波，α初始值应该足够小但不能为零，从(17) 可知，r=0.01~0.02时，f(r)≈1，此时可以忽略α参数对衰减系数的影响，因此α=(0.01~0.02)ωε0.

在YEE差分网格剖分中，电场和磁场在空间上相差半个网格，比例伸缩因子Sei、Smi的参数取值也因此相差半个网格，并不严格满足Sei=Smi的无反射条件，G≠1，这使得电磁波在CPML各层分界面上有反射.假定分界面1和分界面2的比例因子分别为S1ei、S2mi，则垂直入射的电磁波在相邻两层介质分界面的反射系数为

式(18) 表明反射系数与入射角度无关，只与频率有关.在时域瞬变电磁法三维有限差分正演计算中，主要考虑空气介质(电导率为零)和大地介质(电导率高于10-4 S/m)，在频率低于106 Hz的似稳态条件下的吸收效率.忽略α的影响，在满足无反射条件时有σpew/ε0=σpmw/μ0，磁场比例因子Smw中的μ0可以替换成ε0，电场比例因子Sew和磁场比例因子Smw可以表示为

σp和κ在分界面两侧都满足各自的离散函数分布，电场比例伸缩因子Se的参数σpe、κe和磁场比例伸缩因子Se的参数σpm、κm相差半个网格.假设κ的分布函数为g，σp的分布函数为f，见图 2，x1、x2、x3、x4之间的间隔为半个网格.

分别将κ1m、κ2e、κ2m在κ1e处、将σ1pm、σ2pe、σ2pm在σ1pe处用一阶泰勒展开，代回式(18)，衰减系数化简可得：

式中s1= g’(x1)Δ/2、t1 =f’(x1)Δ/2, (20) 的简化条件是s1、t1为常数或者同分布.s1、t1越小则反射系数越小，在无耗介质和似稳态条件下，κ对低频衰减不敏感，在CPML各层上取恒定常量κ=1，此时s1=0；为减小t1，将Roden提出的σp和α的分布函数修改为

瞬变电磁中心回线法是在地表铺设不接地矩形发射线框发射脉冲电磁场，通过测量二次场响应特征来分析地下介质的电性参数和空间形态.发射电磁场一部分向空中逸散，一部分向地下传播，数值模拟需要处理无耗的空气和有耗的大地两种边界.

为了分析CPML边界条件对空气介质的吸收效率，如图 3所示，建立516×516网格的空气模型，网格大小10 m，CPML厚度8个网格，相对介电常数εr=1，电导率σ=0，发射线框大小110 m×110 m，发射电流为宽度2us的高斯脉冲.线框中心距离CPML右侧边界15个网格，采样点距离边界6个网格，这样在前500个时间步只有右侧边界有反射.作为对比，将网格扩大到1200×1200，将发射线框放在网格中心，测点与线框的位置关系不变，前500个时间步没有任何方向的反射干扰，以此作为标准值对有反射的参考点进行比对，分析CPML的反射误差.误差计算公式为(Gedney，1996)：

式中Xi为原计算模型采样点的场值，Xiref为参考模型空间对应点的场值，Xiref_max为参考模型空间对应点场值的最大值.

根据(21) 式，CPML吸收边界参数取m=3，α=10-6，u0在计算Se时取0.75，在计算Sm时取0.25.为了验证理论分析结果，κ值分别取1、5、10时，采样点处Hz的反射误差如图 4所示.

可以看出在无耗介质中，κ分布越平缓反射误差越小，这和理论分析的结果是一致的.线框与采样点的距离以及采样点与CPML边界距离不同时，也有同样的结果.为了避免CPML对发射源处场量剧烈变化的影响，空气层的厚度至少需要4~6个网格.

大地介质的视电阻率多数情况下低于104 Ω·m，在频率低于106 Hz时仍然满足似稳态条件，典型情况下，大地视电阻率取100 Ω·m.建立216×216网格的大地介质模型，模型和发射线框的参数与空气层参数相同，设置4个采样点分别距离CPML边界5、15、25、35个网格.参考模型大小为1500×1500网格，将发射线框置于网格中心，采样点与线框的位置关系不变，这样电磁波从发射再返回到各个采样点，即使没有CPML吸收层，实际走过的距离接近3000个网格(30000 m)，也能够使各个频率的电磁波充分衰减，以此模型的Hz场量作为参考标准值.在不同频率下，距离CPML吸收边界不同采样点的反射误差见图 5，图中d为采样点距离CPML边界的网格数.

建立136×136×92网格的三维地空模型(模型1)，其中地上空气层厚度6个网格，大地层厚度70个网格，其余参数与空气介质的模型参数相同.在埋深10个网格处放置21×21×5的高阻空气空腔.在地表线框中心和距离地下高阻空气层右边界5个网格处分别设置采样点P1、P2，见图 5，采样点距离CPML边界都大于35个网格，反射误差小于10-5，将此采样点数值作为标准值.为考察模型空间变小时计算误差，在空气层厚度不变的情况下，将模型缩小到96×96×72网格(模型2) 和76×76×52网格(模型3).分别计算三种模型，分析模型2、模型3的采样点相对模型1的误差，其中P2采样点的反射误差见图 6，而P1点的反射误差主要来自顶部空气层，三种模型的反射误差和图 3的结果相同，模型2、模型3与模型1之间的相对误差为0.

在模型2中，P2采样点距离CPML边界15个网格，迭代104步后与模型1的相对误差小于10-4，网格点数为原模型的1/3；模型3中，P2采样点距离CPML边界5个网格，迭代104步后误差小于10-2, 网格点数为原模型的1/6.如果将高阻空气空腔换为低阻模型，由于电磁波受低阻的屏蔽和吸引作用，反射误差小于高阻模型.

CPML的反射误差随着频率的下降而增加，在瞬变电磁法数值模拟时，大地介质中的涡流相当于向下不断扩散的烟圈电流，高频成分衰减快，后期低频成分起主要作用.式(21) 中，α的取值与最大发射频率相关，在CPML厚度较小时，α的取值分布很难有效吸收宽频电磁波，特别是晚期的低频波.如果计算的时间比较长，迭代次数多，反射误差会更加严重，因此需要增加CPML的厚度.模拟结果表明，大地视电阻率为100 Ω·m情况下，计算时长超过3 ms或者迭代次数超过105步的情况下，CPML的厚度需要增加到8~12层，误差才会小于4%.如前分析，截断边界距离地下模型边界在10~15个网格就可以满足精度要求，可以大幅度缩小计算空间.如果误差不能满足精度要求，则适当增加CPML的厚度.孙怀凤等(2013)中的均匀大地介质模型大小301×301×100，网格大小为10 m，计算了均匀半空间的感应电动势衰减曲线.作为对比，采用(21) 式的CPML参数分布，建立76×76×62的网格模型，其他参数相同，计算结果和孙怀凤等(2013)的一致，而计算空间则要小的多.

本文推导了CPML吸收边界参数在空气介质和大地介质中的吸收系数和反射系数，分析了参数分布对计算结果的影响，改进了参数分布，能够有效吸收低频电磁波.在瞬变电磁数值模拟中，影响计算效率的主要因素是网格分辨率和计算空间的大小.在网格分辨率一定的情况下，计算空间太大导致网格节点数量过大，计算耗时过长，限制了时域瞬变电磁三维数值模拟的实际应用.瞬变电磁早、晚期响应的幅值相差大于106倍，并且迭代次数多，边界截断误差会对晚期数据有很大影响，通过对CPML吸收边界在瞬变电磁法中的应用参数变化的分析和数值模拟，结果表明，采用CPML吸收边界条件可以在满足精度要求的情况下有效缩小计算空间，提高计算效率，也可以对低频条件下大地电磁数值模拟提供参考.