矩量法由于积分方程自动满足辐射条件，不但适用于任何形状的复杂目标，而且数值计算精度与效率高，同时避免了色散误差传播积累和网格截断误差等影响，目前已经在天线辐射、复杂散射体的电磁散射和电磁兼容等实际工程领域获得了广泛的应用。但是，在运用该方法求解开域问题，如散射和辐射问题时，将积分方程离散后，得到的阻抗矩阵是满秩的，所消耗的计算机内存和计算复杂度分别为 $o({N^2})$ 和 $o({N^3})$ ，其中N为积分方程离散后的自由度，即未知量的个数^[1]。随着N的增加，会使得填充阻抗矩阵中元素所需的计算量增大，导致计算机的内存需求增加也是非常高的，这就大大限制了矩量法对于电尺寸较大目标的应用^[2]。为了解决该问题，近年来学者们提出了很多快速算法：其中一类是基于物理意义的快速算法，包括快速多级子算法(multilevel fast multipole algorithm，MLFMA)、基于快速傅里叶变换的方法(fast Fourier transform，FFT)，如自适应积分法(adaptive integral method，AIM)。这类算法在内存和计算复杂度方面具有优势，但是这些元素都要依赖于选择合适的元素来等效代替组的物理性质，如多级子、等效表面或者卷积基/测试基函数^[3]。比如MLFMA要求将格林函数用加法定理展开，当格林函数具有复杂形式时，如分层媒质中的格林函数，是很难做到的。

另一类是基于代数意义的快速算法被称为基于矩阵压缩的算法，如多层矩阵分解方法(multilevel matrix decomposition algorithm，MLMDA)、多层QR分解、多层UV分解、快速双多层正交化(multilevel Gram Schmidt，MGS)分解、自适应交叉算法(adaptive cross approximation，ACA)算法等^[4-5]。这类算法只依赖于对矩阵的数学操作，不涉及积分核函数和网格的投影操作。但大多数的算法分解需要知道完整的子矩阵信息，这样就导致计算量的增大。例如，多层QR分解需要 $o({N^2})$ 的计算量。而ACA最大的优势在于它不依赖于格林函数和积分方程，其分解过程只需要原矩阵的部分行和列，不需要提前知道要分解的子矩阵的全部信息，以行、列的外积来近似代替原矩阵，在分解的过程中自适应地得到矩阵的秩。ACA算法是针对稠密矩阵提出来的，其数学原理是利用线性相关性来对低秩矩阵进行压缩。ACA算法由Bebendorf于2000年首次提出^[6]，并被应用到多种数值算法中。2005年由K.Z.Zhao等^[7]引入到矩量法中，用于分析电磁辐射和散射问题。但由于传统ACA算法所产生的矩阵并不是酉矩阵，进而导致矩阵压缩不彻底。鉴于此，文中提出结合QR分解和奇异值分解(singular value decomposition，SVD)产生一种新的快速算法ACA-SVD来对其进一步压缩，并采用MATLAB和C语言混合编程来实现具体的计算仿真。

1 理论分析 1.1 矩量法原理

对于三维理想导体目标电磁散射问题，在入射场Eⁱ的照射下，根据等效原理，满足理想导体表面切向电场为零的边界条件，可得到：

$- {E^s}(r){|_{\tan }} = {E^i}(r){|_{\tan }}$

(1)

式中r为表面S上的场向量，当平面波照射到理想导体上时，导体表面会产生感应电流，感应电流可用RWG基函数 $\{ {f_n}(r)\} $ 展开，采用伽略金法，可将式(1)转化为矩阵方程：

$\mathit{\boldsymbol{ZI}} = \mathit{\boldsymbol{V}}$

式中：Z表示阻抗矩阵，I是待求电流向量，V是激励矢量。

文中以矩量法为基础，运用文献[8]所述的方法，得到阻抗矩阵元素 ${\mathit{\boldsymbol{Z}}_{mn}}$ 的表达式为

$\begin{aligned}{\mathit{\boldsymbol{Z}}_{mn}} = & {l_m}[{\rm{j}}w(A_{mn}^ + \cdot \frac{{\rho _m^{c + }}}{2} + A_{mn}^ - \cdot \frac{{\rho _m^{c - }}}{2}) + \varPhi _{mn}^ - - \varPhi _{mn}^ + ]\\{{V}_m} = & {l_m}(E_m^ + \cdot \frac{{\rho _m^{c + }}}{2} + E_m^ - \cdot \frac{{\rho _m^{c - }}}{2})\\{A}_{mn}^ \pm = & \frac{\mu }{{4{\rm{\pi }}}}\int_S {{f_n}(r')\frac{{{{\rm{e}}^{ - {\rm{j}}kR}}{m^ \pm }}}{{R_m^ \pm }})} {\rm{d}}S'\\\varPhi _{mn}^ \pm = & - \frac{\mu }{{4{\rm{\pi j}}w\varepsilon }}\int_S {{\nabla _S}^\prime \cdot {f_n}(r')} \frac{{{{\rm{e}}^{ - jkR}}{m^ \pm }}}{{R_m^ \pm }}{\rm d}S'\end{aligned}$

式中：l_m表示第m对三角面元的公共边长度； $\rho _m^{c + }$ 、 $\rho _m^{c - }$ 分别表示正、负三角形的自由顶点到中心的坐标矢量； $E_m^ \pm $ 是三角形的中心点处的电场强度；V_m为激励矢量元素； $\varPhi _{mn}^ \pm $ 为标量磁位函数； ${A}_{mn}^ \pm $ 为矢量磁位函数。

1.2 再压缩自适应交叉近似算法 1.2.1 自适应交叉近似算法

由于矩量法通过基函数和权函数来表达场点和源点间的相互作用，从而形成阻抗矩阵。根据格林函数表达式，两点距离越近，其相互作用所对应的元素值就相对越大；两点距离越远，其相互作用所对应的元素值就相对越小^[9]。自适应交叉近似算法就是利用远场组所形成的阻抗矩阵具有低秩的特点，对阻抗矩阵元素进行分解，如式(2)^[10]：

${\tilde{Z}^{m \times n}}{\rm{ = }}{\mathit{\boldsymbol{U}}^{m \times k}}{\mathit{\boldsymbol{V}}^{k \times n}}{\rm{ = }}\sum\limits_{i = 1}^k {u_i^{m \times 1}v_i^{1 \times n}} $

(2)

式中：k为矩阵 ${\mathit{\boldsymbol{Z}}^{m \times n}}$ 的有效秩， ${\mathit{\boldsymbol{U}}^{m \times k}}$ 和 ${\mathit{\boldsymbol{V}}^{k \times n}}$ 为压缩后得到的满秩的矩阵。

ACA算法选取的秩由式(3)所示：

$\parallel {\mathit{\boldsymbol{R}}^{m \times n}}\parallel = \parallel {\mathit{\boldsymbol{Z}}^{m \times n}} - {\tilde{Z}^{m \times n}}\parallel \leqslant \varepsilon \parallel {\mathit{\boldsymbol{Z}}^{m \times n}}\parallel $

(3)

式中：ε为误差迭代门限，R为误差矩阵， $\parallel \cdot \parallel $ 是矩阵的Frobenus范数。在实际的电磁问题中，一般的秩 $k{\rm{ < min(}}m{\rm{,}}n{\rm{)}}$ ，因此，ACA算法只需要保存 ${\mathit{\boldsymbol{U}}^{m \times k}}$ 和 ${\mathit{\boldsymbol{V}}^{k \times n}}$ 这2个维度较小的矩阵即可，这样矩阵的存储量就由 $m \times n$ 降到 $k \times (m + n)$ ^[11]，计算复杂度为 $o({k^2}(m + n))$ 。

ACA算法流程分为初始化部分和第k次的迭代^[7]。

初始化部分：

1) 初始化第一个行索引I₁=1，初始化 $\tilde{Z}{\rm{ = 0}}$ ；

2) 初始化误差矩阵的第一行： $\tilde{R}\left( {{I_1},:} \right)\! =\! \mathit{\boldsymbol{Z}}\left( {{I_1},:} \right)$ ；

3) 在第一行搜索最大值确定第一个列标号J₁： $\left| {\tilde{R}\left( {{I_{\rm{1}}},{J_{\rm{1}}}} \right)} \right| = {\rm{max}}j\left| {\tilde{R}\left( {{I_{\rm{1}}},{J_{\rm{1}}}} \right)} \right|$ ；

4) V矩阵第一行表示为 ${\mathit{\boldsymbol{V}}_1} = \tilde{R}\left( {{I_{\rm{1}}},:} \right)/\tilde{R}\left( {{I_{\rm{1}}},{J_{\rm{1}}}} \right)$ ；

5) 初始化误差矩阵第一列 $\tilde{R}\left( {:,{J_{\rm{1}}}} \right) = \mathit{\boldsymbol{Z}}\left( {:,{J_{\rm{1}}}} \right)$ ；

6) U的第一列为 ${\mathit{\boldsymbol{U}}_{\rm{1}}} = \tilde{R}\left( {:,{J_1}} \right)$ ；

7) $||{\mathit{\boldsymbol{Z}}^{\left( {\rm{1}} \right)}}|{|^2} = ||{\mathit{\boldsymbol{Z}}^{\left( 0 \right)}}|{|^2} + ||{u_1}|{|^2}||{v_1}|{|^2}$ ；

8) 搜索第1列最大值确定第2个行标号I₂： $\left| {\mathit{\boldsymbol{R}}\left( {{I_2},{J_1}} \right)} \right| = \mathop {\max }\limits_j \left( {\left| {\tilde{R}\left( {i,{J_1}} \right)} \right|} \right)$ ， $i \ne {I_{\rm{1}}}$ 。

第k次迭代：

1) 更新误差矩阵的 ${{I}_k}$ 行： $R({I_k},:) = \left| {Z({I_k},:) - } \right.$ $\left. {\displaystyle\sum\nolimits_{l = 1}^{k - 1} {({u_l}} {)_{{I_k}}}{v_l}} \right|$ ；

2) 搜索第I_k行中最大值，确定第k列的标号J_k。使得： $\left| {\left. {\mathit{\boldsymbol{R}}({I_k},{J_k})} \right|} \right. = \mathop {\max }\limits_j (\left| {\left. {\mathit{\boldsymbol{R}}({I_k},j)} \right|} \right.),j \ne {J_1},{J_2}, \cdot \cdot \cdot ,{J_k}$ ；

3) V矩阵第k行 ${\mathit{\boldsymbol{v}}_k} = \mathit{\boldsymbol{R}}({I_k},:)/\mathit{\boldsymbol{R}}({I_k},{J_k})$ ；

4) 更新误差矩阵第J_k列 $R(:,{J_k}) = {Z}(:,{J_k}) - $ $\displaystyle\sum\nolimits_{l = 1}^{k - 1} {({u_l}} {)_{{J_k}}}{v_l}$ ；

5) 更新U矩阵的每一列 ${\mathit{\boldsymbol{u}}_k} = \mathit{\boldsymbol{R}}(:,{J_k})$ ；

6) $|{{Z}^{\left( k \right)}}|{|^2} \!=\! ||{{Z}^{\left( {k - 1} \right)}}|{|^2} \!+\! 2\displaystyle\sum\nolimits_{l = 1}^{k - 1} {\left| {{u}_j^{\rm{T}}{{u}_k}} \right| \cdot \left| {{v}_j^{\rm{T}}{{v}_k}} \right|} \!+\! ||{{u}_k}|{|^2}||{{v}_k}|{|^2}$

7) 检查收敛性：如果 $||{\mathit{\boldsymbol{u}}_k}||||{\mathit{\boldsymbol{v}}_k}|| \leqslant \varepsilon ||{\mathit{\boldsymbol{Z}}^{\left( k \right)}}||$ ，则迭代结束，否则开始寻找下一行；

8) 寻找下一个行标号 ${I_{k + 1}}$ ： $\left| {\left. {R({I_{k + 1}},{J_k})} \right|} \right. = \mathop {\max }\limits_i $ $\left| {\left. {R(i,{J_k})} \right|} \right.$ ， $i \ne {I_1},{I_2}, \cdot \cdot \cdot ,{I_k}$ 。

在上述迭代过程中，步骤7)的终止条件是根据秩r的选取由式(3)所得。

1.2.2 再压缩技术

由于经过ACA算法得到的2个酉矩阵 ${\mathit{\boldsymbol{U}}^{m \times k}}$ 和 ${\mathit{\boldsymbol{V}}^{n \times k}}$ 并不是正交的，所以接下来采用再压缩技术对解后的低秩矩阵中存在冗余的矩阵进一步压缩，从而降低内存，提高计算速度。算法的具体原理如下^[7]：

首先需要对经过ACA算法得到的 ${\mathit{\boldsymbol{U}}^{m \times k}}$ 和 ${\mathit{\boldsymbol{V}}^{n \times k}}$ 进行QR分解：

$\begin{aligned}& {\mathit{\boldsymbol{U}}^{m \times k}} = \mathit{\boldsymbol{Q}}_U^{m \times k}\mathit{\boldsymbol{R}}_U^{k \times k}\\& {\mathit{\boldsymbol{V}}^{n \times k}} = \mathit{\boldsymbol{Q}}_V^{n \times k}\mathit{\boldsymbol{R}}_V^{k \times k}\end{aligned}$

然后对经过QR分解后得到的2个上三角阵的乘积进行SVD分解：

$\mathit{\boldsymbol{R}}_U^{k \times k}{(\mathit{\boldsymbol{R}}_V^{k \times k})^{\rm{T}}} = {\hat{U}^{k \times r}}\hat{S}{\hat{V}^{r \times k}}$

这个部分采用的是压缩的SVD分解算法，在所有的特征值中满足 ${\sigma _r} < \varepsilon {\sigma _1}$ 的奇异值将会被舍去，只保留高奇异值及所对应的奇异分量来实现对矩阵的低秩近似。其中σ_r表示 $\mathit{\boldsymbol{R}}_U^{k \times k}{(\mathit{\boldsymbol{R}}_V^{k \times k})^{\rm{T}}}$ 中的第r个特征值，一般来说此处选择的阈值是ACA算法选取阈值的10倍^[2]。经过上述操作，由ACA算法得到的2个矩阵 ${\mathit{\boldsymbol{U}}^{m \times k}}$ 和 ${\mathit{\boldsymbol{V}}^{n \times k}}$ 便转化为 ${\hat{U}^{m \times r}}$ 和 ${\hat{V}^{r \times n}}$ ，其中r<k，从而在ACA算法的基础上进一步降低了有效秩，降低了内存，节省了计算时间。

算法的流程图如图1所示。

1.3 近场预处理技术

对于式(2)阻抗矩阵元素Z_mn的计算，采用MATLAB和C语言的混合编程来实现。MATLAB具有强大的数值计算功能，但运行速度较慢。而C语言可对操作系统和应用程序以及硬件进行直接操作，优于其他解释型高级语言。因此，将MATLAB与C语言混合编程可以提高阻抗矩阵的计算效率。在完成阻抗矩阵填充之后，求解矩阵方程时采用广义最小余量法(generalized minimal residual，GMRES)，这是一种动态迭代解法，也称为Krylov子空间迭代解法，求解时可将计算量降为 $o({N^2})$ 。对于迭代问题的求解，计算时的收敛速度很大程度上取决于系数矩阵的谱特性。而预处理技术的引入可以改善矩阵的谱特性，从而达到提高迭代算法收敛性的目的。这里采用ILU预处理技术改善矩阵的性态，并结合Krylov子空间迭代法来求解矩阵方程。

2 数值算例与结果

本文算例程序均基于MATLAB2010b仿真平台，应用具有良好稳定性的GMRES算法进行迭代求解，收敛精度设置为10^–3，所有算例在Intel Pentium G2020 2.9 GHz CPU、1.95 GB内存的台式计算机上运行。

算例1　分析在平面波入射频率f=300 MHz，入射角度为(θ, ϕ) = (0°, 0°)，双站RCS散射波的观察角度为(θ_s,ϕ_s) = (0°~180°, 0°)。球体表面按照λ/10标准剖分，剖分得到的未知量个数为3 072，其中λ为自由空间的波长的条件下。分别采用传统矩量法、结合ACA的矩量法以及结合ACA-SVD的矩量法进行求解，得到的GMRES迭代求解球体的迭代次数和迭代误差曲线如图2所示。根据计算结果，比较ACA阈值和SVD阈值选取不同值时内存的压缩效果。

通过图2可以看出，ACA算法迭代误差门限选用ε=10^–3或是ε=10^–6，求解过程均在51次时到达收敛。下面给出了2个算例的内存消耗，如表1所示。

通过表1可以看出，ACA算法迭代误差门限并不是越小越好。这是因为门限越小，所需要的内存越多，算得越慢。

表2给出了ACA阈值和SVD阈值选取不同阈值时内存压缩效果的比较。

表2中的相对误差是指 $\left| {\displaystyle\frac{{{J_{{\rm{com}}}}{\rm{ - }}{J_{{\rm{uncom}}}}}}{{{J_{{\rm{uncom}}}}}}} \right|$ ，其中 ${J_{{\rm{com}}}}$ 指的是压缩后计算的电流密度， ${J_{{\rm{uncom}}}}$ 是指未经压缩矩量法计算的电流密度值。由表2可以看出，经过压缩和未经过压缩得到的电流密度相差很小，并且当SVD阈值是ACA阈值的10倍时，压缩的效果最好。

算例2　分析理想导体球的双站RCS，平面波入射频率f=300 MHz，入射角度为(θ, ϕ) = (0°, 0°)，双站RCS散射波的观察角度为(θ_s, ϕ_s)=(0°∼180°, 0°)。球体表面按照λ/10标准剖分，剖分得到的未知量个数为3 072，其中λ为自由空间的波长。分别采用传统矩量法、结合ACA的矩量法以及结合ACA-SVD的矩量法进行求解，其中，ACA算法迭代误差门限分别为ε=10^–3和ε=10^–6的条件下。SVD分解所选取的误差分别为ε=10^–2和ε=10^–5，求解得到的RCS对比曲线如图3所示。从图中可以看出2种方法结果吻合的良好。

由图3可见，分别采用结合了ACA压缩算法的矩量法以及结合了ACA-SVD压缩算法的矩量法求解双站RCS，结果与Mie级数解吻合良好，说明ACA和ACA-SVD算法在对矩阵进行压缩的同时并未损失矩量法的计算精度。

图4显示了PEC球体的剖分细度不一样时，产生不同未知量的未知数，此时运算过程所占内存的变化趋势。

通过图4可以看出，矩量法的存储量的增加趋势为 $o({N^2})$ ，自适应交叉近似算法对内存的要求为 $o({N^{4/3}}\log (N))$ ^[7]。

通过以上算例可以得知，ACA和ACA-SVD算法的运用，在不影响MOM计算精度的前提下，可以大大减少计算的存储空间，提高计算速度。例如，对于PEC球体，在ε=10^-3的条件下，相比于传统矩量法，ACA和ACA-SVD算法的引入，分别可以减少59.25%的存储空间和78.10%的存储空间，从而可以加速后续的矩阵向量乘的计算。

3 结论

文中应用ACA和ACA-SVD算法求解导体目标双站RCS。在远场组作用区，利用该方法压缩阻抗矩阵，减少矩阵填充时间、节省内存；在近场作用区，采用传统矩量法精确填充阻抗矩阵元素来提高计算效率。通过数值计算结果可以看出：

1) ACA算法可以有效地压缩阻抗矩阵，加速矩阵填充和矩阵向量乘。在保证一定计算精度的前提下，相比于传统的ACA算法，ACA-SVD的方法可以更快速地进行阻抗矩阵的填充，压缩效率更高，并通过仿真结果发现矩阵的规模越大，剖分的未知量数目越多，压缩效果越明显。

2) 验证了ACA-SVD算法的优越性。在基于RWG基函数的矩量法中加入ACA-SVD算法，大大降低计算的存储量和复杂度，提升矩量法的计算效率。

3) 下一步可以尝试将ACA-SVD算法与其他阻抗矩阵填充算法混合，进一步提高算法的性能，节省计算机内存和计算时间。该方法对研究复杂目标电磁散射特型有一定的参考价值。