引　言

近年来，电磁散射在地球物理勘探^[1]、遥感^[2]、微带天线^[3]、微波成像^[4]、无损检测^[5]等民用和军用领域的应用愈发广泛. 在这些领域的电磁散射研究过程中，经常需要考虑一种层状介质的模型. 由于异常体的非均匀性，传统的表面积分方程^[6]是不适用的，这种情况下应用体积分方程更为合适. 解决非均匀物体电场积分方程(electric field integral equation, EFIE)最直接的方法是矩量法(method of moments, MoM)^[7]. 传统求解体积分方程的一类方法是将迭代方法与快速傅里叶变换(fast Fourier transform, FFT)结合，比较典型的有共轭梯度快速傅里叶变换(conjugate gradient fast Fourier transform, CG-FFT)方法、双共轭梯度快速傅里叶变换(biconjugate gradient fast Fourier transform, BCG-FFT)方法、稳定双共轭快速傅里叶变换(stabilized biconjugate gradient fast Fourier transform, BCGS-FFT)方法等. 由文献[8-9]中几个均匀背景下的算例可知，BCG-FFT方法比CG-FFT方法更高效. BCGS-FFT方法^[10]在求解均匀背景下大规模散射问题比BCG-FFT方法和CG-FFT方法精确度更高，效率更快^[11].

但是在求解EFIE的过程中，我们需要得到空域格林函数，这个过程一般是由索末菲积分完成的. 在索末菲积分中，尽管层状介质谱域格林函数是解析的^[12]，但是这个积分通常是震荡的、奇异的、收敛很慢甚至是不收敛的，并且用常规的积分方法是费时费力的. 因此，很有必要加速索末菲积分的收敛速度. 目前计算索末菲积分的方法主要有复平面积分法^[13]、复镜像法^[14]、近似法^[15]、最陡下降法^[16]，以及对索末菲尾部积分应用分部外推法^[17-19]等等.

本文我们在复平面中采用一种椭圆路径来避开格林函数的奇异点并用分部外推法加速索末菲尾部积分的收敛速度. 这种方法的优势是在尾部积分部分沿实轴积分时贝塞尔函数只有实参，并且不需要确定格林函数极点的精确位置. 该方法比复镜像^[20]等近似方法更精确可靠，比直接积分法^[21]速度更快，并且不用耗时去提取奇异点^[22]. 一般来讲这种积分方法是普遍适用的，但应用该方法时尽量不要让源点与场点之间的水平距离大于100个波长，否则会使贝塞尔函数在初始区间中有过多的震荡，影响最终的积分结果^[23]. 在各类分部外推法中，目前使用较多的有Euler变换^[24]、Aitken迭代变换^[25]、加权平均^[26]、Shanks变换^[27]，以及一般Levin变换^[28]结合W算法^[29]等，文献[17]中的数值测试部分对比了这些方法，因为索末菲积分的渐进现象，一般Levin变换结合W算法与加权平均算法是加速索末菲尾部积分最通用有效的方法. 本文的工作主要分为三个部分：①在复平面中使用一种简单实用的椭圆路径来避开格林函数的所有奇异点. ②将分部外推法应用于传统计算索末菲尾部积分的高斯积分法^[30]中，以此来计算多层介质并矢格林函数，填充并矢格林函数矩阵. ③将分部外推法与BCGS-FFT方法结合求解层状介质中的电磁散射问题，提高计算效率.

1   EFIE

如图1所示的分层介质的电磁散射模型中共有n层，在第m层中嵌入非均匀电介质物体. 因此，在第j层中总电场${E_j}({\boldsymbol{r}})$等于入射电场${\boldsymbol{E}}_j^{\rm{in}}({\boldsymbol{r}})$与散射电场${\boldsymbol{E}}_j^{{\rm{sc}}}({\boldsymbol{r}})$的和，即：

式中：${\boldsymbol{r}}$是空间中的位置矢量；${\boldsymbol{E}}_j^{\rm{in}}({\boldsymbol{r}})$是在j层没有异常体，仅仅有分层介质情况下的入射电场. 关于散射场${\boldsymbol{E}}_j^{{\rm{sc}}}({\boldsymbol{r}})$，我们可以用磁矢势${{\boldsymbol{A}}^{jm}}({\boldsymbol{r}})$和标量势${\phi ^{jm}}(\boldsymbol r)$来表示：

这里磁矢势${{\boldsymbol{A}}^{jm}}({\boldsymbol{r}})$是由嵌入第m层的非均匀物体产生的感应电流源${\boldsymbol{J}}$引起的，表达式为

式中：${\mu _j}$是第j层的介磁常数；$V$是异常体体积；${{\boldsymbol{G}}^{jm}}({\boldsymbol{r}},{\boldsymbol{r}}')$是传统形式的并矢格林函数^[12]；${\boldsymbol{J}}$可以表示为${\boldsymbol{J}}{\rm{ = i}}\omega \chi {\boldsymbol{D}}$，其中$\chi ({\boldsymbol{r}})=(\tilde{\varepsilon }({\boldsymbol{r}})-{\tilde{\varepsilon }}_{m})/\tilde{\varepsilon }({\boldsymbol{r}})$是异常体的对比度函数，${\tilde \varepsilon _m} = {\varepsilon _m} - {\rm{i}} {\sigma _m}/\omega$是复介电常数，${\varepsilon _m}$与${\sigma _m}$分别是第m层介电常数、磁导率，$\tilde \varepsilon ({\boldsymbol{r}}) = \varepsilon ({\boldsymbol{r}}) + \sigma ({\boldsymbol{r}})/ ({\rm{i}}\omega )$，电通量$D({\boldsymbol{r}}) = \tilde \varepsilon ({\boldsymbol{r}}){\boldsymbol{E}}({\boldsymbol{r}})$. 通过洛伦兹条件，我们可以把磁矢势${{\boldsymbol{A}}^{jm}}({\boldsymbol{r}})$和标量势${\phi ^{jm}}(\boldsymbol r)$联系起来，即：

式中，${k_j} = \omega \sqrt {{{\tilde \varepsilon }_j}{\mu _j}}$为第j层的波数.

将式(2)～(4)代入式(1)中，当${\boldsymbol{r}} \in V$(也就是m=j)时，我们可以得到EFIE：

这里我们可以引入一个$\gamma$算子，即：

可以看出，只要电通量${{\boldsymbol{D}}_m}({\boldsymbol{r}})$确定，我们就可以通过式(2)～(4)求出第j层任意位置的散射场. 我们的方案是使用Zwamborn和Berge提出的屋顶基函数弱离散EFIE方程的方法^[31]，并用BCGS代替CG方法来求电通量${{\boldsymbol{D}}_m}({\boldsymbol{r}})$.

图  1  平面分层介质中电磁散射的非均匀物体的典型几何模型

Fig.  1  A typical geometric model of an inhomogeneous object with electromagnetic scattering in a plane layered medium

下载: 全尺寸图片 幻灯片

2   分层介质格林函数的计算

在式(3)中，我们需要得到空域并矢格林函数${{\boldsymbol{G}}^{jm}}({\boldsymbol{r}},{\boldsymbol{r}}')$，${{\boldsymbol{G}}^{jm}}({\boldsymbol{r}},{\boldsymbol{r}}')$一般可被写为^[12]

谱域格林函数$\tilde G_{xx}^{jm}$的封闭形式可以由传输线理论推导得出^[2]，空域格林函数是由谱域格林函数经过x和y方向上傅里叶逆变换得到的，这与索末菲积分或者汉克尔变换是等价的，即：

式中：n为贝塞尔函数${\rm{J}}_{n}(x)$的阶数；$\;\rho =$$\sqrt {{x^2} + {y^2}}$；${k_\rho } = \sqrt {{k_x} + {k_y}}$，${k_x}$和${k_y}$为波矢量${\boldsymbol{k}}$在x和y方向的分量.

一般来说，传统方法多采用在复平面区域沿着实轴进行积分，这样可以避开贝塞尔函数复杂的参数计算. 但这类方法需要先精确定位平面波的极点，再计算剩余路径的积分值. 这种做法在只有一层结构时比较容易做到，但是在多层结构时会变得非常的耗时. 因此，本文采用一种椭圆路径来避开所有极点^[32]，当然路径的选择并不是唯一的^[33]，但是椭圆是简单并且灵活的路径，这可以使我们不用去计算极点的位置，这种路径可以最小化贝塞尔函数的振荡和表面波极点的影响^[34]. 格林函数的极点在${k_0} < {k_\rho } < {k_0}\sqrt {\max ({\varepsilon _{\rm{r}}}{\mu _{\rm{r}}})}$的范围内，格林函数在这之间沿着实轴是无损的，沿着虚轴是指数变化的. 因此，我们可以得到一条如图2的椭圆积分路径. 椭圆积分路径的中心在$({k_0} + {k_0}\sqrt {\max ({\varepsilon _{\rm{r}}}{\mu _{\rm{r}}})} \;)/2$处，椭圆的短半轴平行于虚轴，与源点和场点之间的水平距离$\,\rho$成反比关系. 这种路径避免了因贝塞尔函数自变量的虚部增长带来的指数变化引起的数值问题. 在数值积分时，可以先采用常规方法将椭圆参数化，为了保证积分的精度，积分点的数量应该随着场点与源点的水平距离$\,\rho$的增大而增多，随后采用高斯积分或其他积分方法. 传统的BCGS-FFT做法是在尾部积分部分继续使用高斯积分，每得到一小段的积分就与前面所得到的积分和来对比，当这段积分值小于积分和与截止误差的积时停止积分. 在本文中，我们在尾部积分部分使用高效计算的Levin变换^[28]分部外推算法.

图  2  索末菲积分路径

Fig.  2  Sommerfeld integral path

下载: 全尺寸图片 幻灯片

3   分部外推算法

关于各种具体的分部外推法在格林函数尾部积分段的表达形式，Michalski已经做了详细的说明^[17]. 本文只介绍用到的Levin变换和Sidi的W算法. 一般Levin变换被定义为一个模型序列，即：

式中：${\omega _n}$为残差估计系数；k为Levin变换的阶数；${c_i}$为未知系数；序列$\left\{ {{x_n}} \right\}$是插值点的辅助序列，即图2中的${\xi _{ - 1}}$，${\xi _0},\cdots$.

${S_n}$为前n项积分值的和，即：

式中，${u_i}$为第i段积分的值，${u_0}$即为图2中${\xi _{ - 1}}$与${\xi _0}$之间的积分值. 令${\omega _n} = {u_n}$，通过式(9)我们可以看出，Levin变换的模型序列里的未知量有待求极限$S$和k个系数${c_{ - 1}},{c_0}, \cdots ,{c_{k - 1}}$. 因此，只需要求出$k + 1$个前n项积分和$\left\{ {{S_n},{S_{n + 1}}, \cdots ,{S_{n + k}}} \right\}$，即可构造一个拥有$k + 1$个方程的方程组，得益于克莱姆法则，我们不需要求出未知系数${c_i}$，就可以获得所求的极限$S$：

式中，

可以看出，这种Levin变换的劣势是非递归的，即每次设置k的值时都必须从新开始计算. 为了解决这一问题，Sidi使用一种W算法使上面的算法变为递归的^[29]. W算法通过引入一个除法差分算子${\delta ^k}$：

这里的${\delta ^0}({u_n}) \equiv {u_n}$，这样式(9)可以改写为

可以发现，式(14)等号右边为$k - 1$次的关于$x_n^{ - 1}$的多项式，所以，在式(14)两边同时应用${\delta ^k}$算子，可以得到

因此，利用${\delta ^k}$算子的线性性质，我们可以把积分极限$S$表示为

初始值$S_n^{(0)} = {S_n}$，经过上面的推导可以看出，极限$S$的最佳近似值为$S_0^{(k)}$. 在数值计算过程中，使用两个可变长度的一维数组即可实现.

4   数值算例

在数值算例部分，我们模拟了三种典型的平面分层结构中嵌入异常体的电磁散射情况来验证上述方法. 随后，将本文提出的分部外推BCGS-FFT方法与传统BCGS-FFT方法在数值精度与计算时间上做了对比. 为了方便起见，将三种算例的源都设为一个单位电偶极子，极化方向为(1, 1, 1)，放置在(0, 0, −0.5)的位置. 所有BCGS的迭代截止误差和索末菲积分的截止误差都设为10⁻⁵. 所有算例的仿真都在Intel(R)Core(TM)i5-8250U的CPU以及16 G内存的工作平台上进行.

4.1   半空间模型中嵌入一个圆柱散射体

算例1：如图3，考虑一个直径$d =$$1.6\;{\rm{ m}}$、高$h = 1.6\;{\rm{ m}}$的圆柱体嵌入到半空间的下层. 半空间模型的分界面为$z = 0$平面，上层空间为真空，下层空间的介质为相对介电常数、磁导率、电导率分别为${\varepsilon }_{\rm{r}2}=4.0$、${\mu }_{\rm{r}2}=1.0$、${\sigma }_{2}{=0}.01 \;{\rm{S/m}}$的类似粘性干土物质. 圆柱散射体的电性参数设置为${\varepsilon }_{\rm{r}}=$$7.0$、${\mu }_{\rm{r}}=1.0$、$\sigma {=}$$0.001\;{\rm{S/m}}$的类似湿花岗岩物质. 圆柱散射体的中心在(0, 0, 2.2)位置.

图  3  直径d = 1.6 m，高h = 1.6 m的圆柱散射体嵌入半空间的下层示意图

Fig.  3  A cylindrical scatterer with a diameter of d = 1.6 m and a height of h = 1.6 m embedded in a semi-spatial background medium

下载: 全尺寸图片 幻灯片

本次仿真的工作频率设置为200 MHz，异常体剖分如图4所示，离散成$50 \times 50 \times 50$的立方体单元，每个正方体边长为0.04 m. 因此每个波长采样点数为37.5，设置61个接收点. 图5和图6展示了两种方法计算出的接收点处散射电场与磁场的精度对比，可以发现新方法可以达到与传统方法相同的计算精度. 但是，在计算机内存消耗相同的情况下，分部外推BCGS-FTT方法与传统BCGS-FFT方法的计算耗时分别为11238.3 s和17819.3 s，加速后的BCGS-FFT可以节省36.7%的时间.

图  4  剖分圆柱体网格区域以及设置61个接收点的位置坐标

Fig.  4  Sectional cylindrical grid area and the coordinates of 61 receiving points

下载: 全尺寸图片 幻灯片

图  5  算例1中接收点处的散射电场对比

Fig.  5  Comparison of the scattered electric field at the receiving points in example 1

下载: 全尺寸图片 幻灯片

图  6  算例1中接收点处的散射磁场对比

Fig.  6  Comparison of the scattered magnetic field at the receiving points in example 1

下载: 全尺寸图片 幻灯片

4.2   三层结构中嵌入一个球体

算例2：如图7所示，我们考虑一个直径$d = 10\;{\rm{ m}}$的球体嵌入到三层空间的中间层中. 第一层介质为真空，中间层为${\varepsilon }_{\rm{r}2}=4.0$、${\mu }_{\rm{r}2}=1.0$、${\sigma }_{2}{=0}.01\; {\rm{S/m}}$类似粘性干土的物质，第三层空间为${\varepsilon }_{\rm{r}3}=3.0$、${\mu }_{\rm{r}3}=$$1.0$、${\sigma }_{3}{=0}.01\; {\rm{S/m}}$类似干砂的物质. 球体还是类似湿花岗岩的物质，其电性参数为${\varepsilon }_{\rm{r}}=7.0$、${\mu }_{\rm{r}}=1.0$、$\sigma {=0}.001\; {\rm{S/m}}$. 第一层、第二层的底部分界面分别在${z}_{1}=0\text{、}{z}_{2}=20.0$处. 球体的中心在(0, 0, 10)处.

图  7  直径为d = 10 m的球体嵌入三层介质空间中间层示意图

Fig.  7  A spherical scatterer with a diameter of d = 10 m embedded in the middle layer of a three-layer medium space

下载: 全尺寸图片 幻灯片

本次仿真工作频率设置为20 MHz，剖分结果如图8所示，每个波长采样点数为75，设置61个接收点. 图9和图10展示了两种方法计算出的接收点处散射电场与磁场的精度对比，可以发现新方法达到了传统方法的计算精度. 但是，分部外推BCGS-FFT方法与传统的BCGS-FFT方法的计算耗时分别为9581.4 s和12065.4 s，加速后的BCGS-FFT可以节省20.6%的计算耗时.

图  8  剖分球体网格区域以及设置61个接收点的位置坐标

Fig.  8  Sectional spheriform grid area and the coordinates of 61 receiving points

下载: 全尺寸图片 幻灯片

图  9  算例2中接收点处的散射电场对比

Fig.  9  Comparison of the scattered electric field at the receiving points in example 2

下载: 全尺寸图片 幻灯片

图  10  算例2中接收点处的散射磁场对比

Fig.  10  Comparison of the scattered magnetic field at the receiving points in example 2

下载: 全尺寸图片 幻灯片

4.3   七层结构中嵌入一个双层结构正方体

算例3：如图11所示，我们考虑一个边长$l = 60\;{\rm{ m}}$的双层正方体嵌入到七层空间的第六层中. 第一层为真空，第二层到第七层的介质分别为：

双层正方体的上下层电性参数分别为：

第一层至第六层的底部分界面分别在${z}_{1}=0,$${z}_{2}=4.0,\;{z}_{3}=8.0,\;{z}_{4}=12.0,\;{z}_{5}=16.0,\;{z}_{6}=86.0$处. 双层正方体的中心在(0, 0, 51)的位置.

图  11  上层高h₁ = 30 m，下层h₂ = 30 m的双层正方体嵌入七层介质空间第六层示意图

Fig.  11  The upper level of the double cube with h₁ = 30 m and the lower level with h₂ = 30 m, and the double cube is embedded into the sixth layer of a seven-layer medium space

下载: 全尺寸图片 幻灯片

本次仿真工作频率设置为2 MHz. 分部外推BCGS-FFT算法与传统BCGS-FFT算法都将双层正方体离散成$60 \times 60 \times 60$的立方体单元，每个正方体边长为1 m，如图12所示. 因此每个波长采样点数为150，设置81个接收点. 图13和图14展示了两种方法计算出的接收点处散射电场与磁场的精度对比. 可以发现新方法与传统方法达到的计算精度相同. 但是，分部外推BCGS-FTT方法与传统的BCGS-FFT方法的计算耗时分别为22651.9 s和36128.9 s，加速后的BCGS-FFT算法可以节省37.3%的计算耗时.

图  12  剖分双层正方体网格区域以及设置81个接收点的位置坐标

Fig.  12  Divide the spheriform grid area and set the coordinates of 81 receiving points

下载: 全尺寸图片 幻灯片

图  13  算例3中接收点处的散射电场对比

Fig.  13  Comparison of the scattered electric field at the receiving points in example 3

下载: 全尺寸图片 幻灯片

图  14  算例3中接收点处的散射磁场对比

Fig.  14  Comparison of the scattered magnetic field at the receiving points in example 3

下载: 全尺寸图片 幻灯片

5   结　论

本文采用了分部外推法计算空域格林函数的尾部积分，且与BCGS-FFT结合以便加速计算非均匀散射体嵌入多层空间的电磁散射问题，并与传统BCGS-FFT方法进行比较. 通过三个不同类型异常体的算例验证了该算法的可靠性. 将分部外推BCGS-FFT算法与传统BCGS-FFT的计算结果进行了对比，验证了算法精确度. 结果表明，本文提出的分部外推BCGS-FFT算法适用于计算异常体在层状结构中的电磁散射问题. 在精度相同的情况下，三个算例分别能节省36.7%、20.6%、37.3%的计算耗时.

值得注意的是，本文公式推导仅适用于各向同性介质，并且不考虑磁介质材料. 此外，文中提出的分部外推BCGS-FFT算法仅适用于散射体嵌入单一层中，暂不适用于跨层散射体. 诸多待解决问题将成为我们今后的研究方向.