一种基于叠层矩阵分解的低频电磁特性仿真方法

摘要

本发明公开了一种基于叠层矩阵分解的低频电磁特性仿真方法，该方法首先采用增量型电场积分方程方法来消除传统电场积分方程方法面临的低频崩溃问题；然后针对增量型电场积分方程的系数矩阵，采用低频多层快速多极子技术来构建其叠层矩阵表达式；进而采用一种重压缩方法对该叠层矩阵进一步压缩，去除冗余信息；最后采用叠层矩阵格式算法对压缩后的叠层矩阵进行上下三角分解，该叠层矩阵分解能够将计算复杂度和内存消耗分别降低到几乎线性，不仅为迭代解法构造了一种预条件技术，而且还构造了一种直接解法。本发明所述方法具有求解速度快、内存消耗低和精度可控的特点，可用于各种低频电磁特性分析。

说明书

技术领域

本发明涉及电磁仿真技术领域，特别是一种基于叠层矩阵分解的低频电磁特性仿真方法。

背景技术

随着微电子技术和制造工艺的发展，在计算机芯片封装设计、高速集成电路设计、天线设计、目标电磁散射特性分析、电磁兼容分析等领域中，出现了越来越多的具有复杂精细结构的目标。这些目标的尺寸远小于波长或几何建模剖分尺寸远小于波长，均属于低频问题的范畴。对于低频问题的分析，准静态方法随着工作频率的升高的误差将越来越大。因此，电磁数值仿真方法在低频问题的分析中发挥着越来越重要的作用。表面积分方程方法是一类具有代表性的电磁数值仿真方法，它只需离散目标表面，因此未知量少、求解效率高。然而，传统的表面积分方程方法在分析低频问题时会出现数值不稳定的情况，让低频电磁仿真无法获得精确解，这一现象被称为“低频崩溃”。许多方法致力于解决低频崩溃问题，例如：基函数Helmholtz分解方法和Calderón预条件方法等。

其中，增量型电场积分方程方法(Augmented Electric Field IntegralEquation,A-EFIE)是一种极其有效的消除低频崩溃的方法(Z.G.Qian and W.C.Chew.Anaugmented EFIE for high speed interconnect analysis.Microwave and OpticalTechnology Letter,vol.50,no.10,pp.2658-2662,2008)。A-EFIE通过在电场积分方程中引入电流连续性方程，对电场积分方程中不平衡的谱分支进行了归一化，从而有效解决了低频崩溃问题。尽管A-EFIE分离了电流矢量和电荷标量，A-EFIE系统矩阵仍然是高度稠密的。迭代解法是求解A-EFIE系统矩阵方程的常用方法，原因在于，迭代解法的主要操作是矩阵-矢量乘，其计算复杂度低且能够通过引入快速算法(如：快速多极子算法)进行加速。然而，迭代解法在面临病态的系统矩阵时，迭代收敛所需步数增多；预条件技术能够在一定程度上提高迭代解法的收敛速度，其中，约束(Constraint)预条件是低频电磁仿真中最为常用的一种预条件技术(Z.G.Qian and W.C.Chew.Fast full-wave surface integralequation solver for multiscale structure modeling.IEEE Transactions onAntennas and Propagation,vol.50,pp.3594-3601,2009)。此外，迭代解法在求解多激励问题时，对不同的右边向量均需重新计算，上述情况都会导致迭代解法效率低下。与迭代解法不同，直接解法不仅不存在收敛速度的问题，而且对于多右边向量问题无需重复计算。然而，由于直接解法巨大的计算消耗，采用直接解法来求解大规模低频问题通常被认为是不可行的。因此，发展一种针对低频电磁问题的高效求解方法具有重要的应用价值。

发明内容

本发明所要解决的技术问题是克服现有技术的不足而提供一种基于叠层矩阵分解的低频电磁特性仿真方法，该方法不仅消除了低频崩溃现象，而且具有求解速度快、内存消耗低和精度可调的特点。

本发明为解决上述技术问题采用以下技术方案：

根据本发明提出的一种基于叠层矩阵分解的低频电磁特性仿真方法，包括以下步骤：

第1步，针对目标表面建立电场积分方程，并引入电流连续性方程，采用三角形单元离散目标表面，用RWG基函数展开表面电流，用冲击基函数展开表面电荷，生成增量型电场积分方程的系统矩阵方程；

第2步，利用矩阵收缩方法或者缩减电荷自由度方法对该系统矩阵方程进行修正，消除电荷中性原理在极低频率时造成的矩阵奇异性；

第3步，建立增量型电场积分方程系统矩阵的叠层矩阵；

第4步，采用叠层矩阵格式的运算法则对该叠层矩阵进行上下三角分解；

第5步，将叠层矩阵的上下三角分解与前后向回代结合，采用预条件的迭代解法或者直接解法进行求解，获得目标表面的电流及电荷分布，通过计算提取所需的电磁特性参数。

作为本发明所述的一种基于叠层矩阵分解的低频电磁特性仿真方法进一步优化方案，第3步中，建立增量型电场积分方程系统矩阵的叠层矩阵的步骤如下：

第(1)步，分别针对定义在棱边上的RWG基函数簇和定义在面片上的冲击基函数簇，采用几何递归二分的方法，构造电流簇树和电荷簇树；

第(2)步，电流簇树和电荷簇树相互作用，构造电流-电流块树、电流-电荷块树、电荷-电流块树和电荷-电荷块树，并引入容许条件来划分容许块和非容许块，其中，非容许块以满阵形式表示，容许块以低秩分解矩阵形式表示；

第(3)步，将系统矩阵中各子矩阵的非零元素分别填入对应的块树中，其中，对于电流子矩阵和电荷子矩阵，采用矩量法构造非容许块，采用低频多层快速多极子技术构造容许块，并利用基于QR分解和奇异值分解的重压缩方法进一步对容许块进行压缩；对于其余子矩阵，将其非零元素对应填入非容许块或容许块，生成系统矩阵的叠层矩阵。

作为本发明所述的一种基于叠层矩阵分解的低频电磁特性仿真方法进一步优化方案，第(3)步中，采用低频多层快速多极子技术构造容许块具体包括，提取低频多层快速多极子的聚合、转移和配置矩阵，将容许块表示成低秩分解矩阵的形式，并通过调节多极子数目来控制低秩分解的精度。

作为本发明所述的一种基于叠层矩阵分解的低频电磁特性仿真方法进一步优化方案，第4步中，叠层矩阵格式的运算法则中涉及低秩分解矩阵运算，通过调节低秩分解矩阵运算中的截断精度来控制叠层矩阵上下三角分解的精度，高精度的分解用于构造直接解法，低精度的分解用于构造预条件来加速迭代解法的收敛，且不同精度的分解具有不同的加速效果。

作为本发明所述的一种基于叠层矩阵分解的低频电磁特性仿真方法进一步优化方案，叠层矩阵上下三角分解的精度通过分解前后矩阵之差的二范数和分解前矩阵的二范数之比来定义，高精度指达到实际仿真精度要求的精度，低精度指低于实际仿真精度要求的精度。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

(1)本发明提出的基于叠层矩阵分解的低频电磁特性仿真方法不仅为低频电磁问题的迭代解法提供一种新型预条件技术，而且为低频电磁问题的求解提供一种直接解法，显著提高现有低频电磁特性仿真方法的性能；

(2)与现有流行的约束预条件技术相比，本发明提出的基于叠层矩阵分解的预条件技术能够显著提高迭代收敛速度，有效解决在系统矩阵性态差时迭代收敛慢甚至不收敛的问题；

(3)本发明提出的基于叠层矩阵分解的直接解法能够将计算复杂度降低到O(Nlog²N)，将内存消耗降低到O(NlogN)(其中，N为未知量的数目)，为低频电磁问题的直接求解提供一种可行途径。

附图说明

图1是基于叠层矩阵分解的低频电磁特性仿真方法流程图。

图2是递归二分定界盒子方法示意图。

图3是电流簇树示意图。

图4是两个簇树相互作用示意图。

图5是电流-电流块树示意图。

图6a是低频单层快速多极子算法转移过程示意图。

图6b是低频多层快速多极子算法转移过程示意图。

图7是飞机目标模型及电流分布图。

图8是不同预条件技术的迭代收敛曲线图。

图9是叠层矩阵分解的计算复杂度曲线图。

图10是叠层矩阵分解的内存消耗曲线图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

图1给出了本发明所述方法的主要流程，下面结合图1对本发明所述方法的步骤作进一步详细描述：

第1步，针对目标表面建立电场积分方程，

其中，i为根号负一，ω为角频率，μ_r为相对磁导率，ε_r为相对介电常数，和分别表示源点和观察点的位置矢量，为格林函数，表示观察点处的表面电流，S'表示目标表面，▽'为梯度算子，表示入射电场矢量。然后，引入描述表面电流和表面电荷之间关系的电流连续性方程，

其中，表示观察点处的表面电荷。

采用三角形单元离散目标表面，用RWG基函数来展开表面电流，用冲击基函数来展开表面电荷，采用伽辽金方法测试后，生成增量型电场积分方程的系统矩阵方程：

其中，表示单位阵，上标T表示矩阵转置，k₀表示自由空间波数，c₀表示光速，η₀表示自由空间波阻抗，j_m和ρ_m分别表示未知的电流系数和电荷系数，矩阵的元素的元素和右边向量的元素b_m具有如下形式：

其中，下标m和n分别表示矩阵的行编号和列编号，和分别表示第m个和第n个RWG基函数，h_m(·)和h_n(·)分别表示第m个和第n个冲击基函数，b_m为入射波和第m个RWG基函数作用生成的右边向量元素，S表示目标的积分表面。矩阵是一个稀疏的关联矩阵，描述了单元和棱边之间的相互关系，其矩阵元素的表达式如下：

第2步，利用矩阵收缩方法或者缩减电荷自由度方法对该系统矩阵方程进行修正，消除电荷中性原理在极低频率时造成的矩阵奇异性。

对于矩阵收缩方法，首先将系统矩阵方程变换为：

令方程(8)的系数矩阵为进而将修正为

其中，trace表示矩阵的迹，e表示单元离散生成的内边数目，p表示单元离散生成的三角形面片数目。；是一个归一化的向量：

其中，的前e个元素是0，后p个元素是这样，基于矩阵收缩方法的增量型电场积分方程的系统可写为：

对于缩减电荷自由度方法，只需在目标的每个独立部分中去除一个电荷未知量，即在每个独立部分让表示电荷的面片数量从p缩减为p-1，从而消除矩阵的奇异性。

第3步，建立增量型电场积分方程系统矩阵的叠层矩阵表达式。进一步地，该过程分三步进行：

第(1)步，构造电流簇树和电荷簇树。簇树T_I一般通过递归地细分一组有限指标簇I＝{1,2,…,N}来构造。对于电流簇树，I表示定义在棱边自由度上的RWG基函数簇I_j；对于电荷簇树，I表示定义在面片自由度上的冲击基函数簇I_ρ。采用基于定界盒子的递归二分方法来构造具有二叉树结构的电流簇树和电荷簇树以构造电流簇树为例，如图2所示，首先以定界盒子包围整个目标表面，然后不断沿坐标轴三个方向对定界盒子进行递归二分，这样包含在定界盒子中的RWG基函数簇I_j也被递归二分，这一过程直至达到终止条件时结束。终止条件定义为，定界盒子达到预设尺寸或者定界盒子中的指标个数小于预设值。假定目标在三角形离散后生成25条棱边，生成的电流簇树如图3所示。类似地，电荷簇树也通过这种方式来构造，不同的是递归二分的不是棱边簇，而是面片簇。

第(2)步，构造电流-电流块树、电流-电荷块树、电荷-电流块树和电荷-电荷块树。块树T_I×J由两个簇树T_I和T_J相互作用生成。在伽辽金方法中，T_I×J由两个相同的簇树T_I和T_J生成，其中，T_I可视作原始基函数簇树，T_J可视作测试基函数簇树。T_I×J通过递归细分块I×J来构造，细分过程直至块t×s∈T_I×J满足如下容许条件：

min{diam(B_t),diam(B_s)}≤ηdist(B_t,B_s)>

其中，t表示簇树T_I中的任意簇，s表示簇树T_J中的任意簇，B_t表示包围t的最小定界盒子，B_s表示包围s的最小定界盒子，diam和dist分别表示定界盒子的欧几里得直径和距离；常数η＞0控制容许条件的严格程度，η越小，容许条件越严格。

以构造电流-电流块树为例，原始基函数簇的电流簇树和表示测试基函数簇的电流簇树层层相互作用，如图4所示，从而生成电流-电流块树如图5所示。

容许条件将块树中的所有块划分成两类：容许块和非容许块。满足容许条件的块称为容许块，对于任一容许块具有如下低秩分解的矩阵形式：

其中，和均为满阵，上标p,q和k表示矩阵维度，符号表示复数域。容许块之外的其余块称为非容许块，以满阵形式表示。如图4所示，在中，白色块为容许块，黑色块为非容许块。

类似地，电流-电荷块树电荷-电流块树和电荷-电荷块树均可通过相应的簇树之间的相互作用并选用适当的容许条件来构造。需要指出的是，和是方树，即行簇树和列簇树相同；但和不是方树，原因在于行簇树和列簇树不同。

第(3)步，将系统矩阵中各子矩阵的非零元素分别填入对应的块树中，生成系统矩阵的叠层矩阵表达式。根据(9)式，基于矩阵收缩方法的增量型电场积分方程的系统矩阵包含五个子矩阵：和因此，需要分别构造它们的叠层矩阵表达式和来构造出的叠层矩阵表达式

对于和的构造，需要将中的元素填入电流-电流块树中，中的元素填入电荷-电荷块树中。其中，以满阵形式存储的非容许块采用矩量法直接构造，以低秩分解矩阵形式存储的容许块则采用低频多层快速多极子算法(Low-frequency MultilevelFast Multipole Algorithm,LF-MLFMA)来构造。在簇树层面，传统LF-MLFMA基于八叉树结构，而叠层矩阵算法通常基于二叉树结构，因此需要将八叉树结构转化为二叉树结构。由于八叉树的一层相当于二叉树的三层，因此，当一个n层的八叉树建立时，就相当于建立了一个3n层的二叉树。在块树层面，传统LF-MLFMA和叠层矩阵中，容许块和非容许块的划分方式必须一致，这要求容许条件的选取必须一致。根据传统LF-MLFMA的近远场划分方式可知，通过选取在(12)式所示的容许条件中选取η＝1能够实现LF-MLFMA和叠层矩阵块树的统一。这样，便可采用LF-MLFMA来构造和了。

三维LF-MLFMA的核心转移因子具有如下形式：

其中，父层转移因子，表示子层配置因子，表示子层转移因子，表示子层聚合因子。L_i为多极子模式数。这里，其中P表示多极子的数目，用来控制展开形式和转移因子的精度。(14)式可写成矩阵形式如下：

其中，和分别表示父层转移矩阵、子层配置矩阵、子层转移矩阵和子层聚合矩阵。考虑多层的情况，如图6a、图6b所示，每一层的矩阵形式如下：

第l层：

第l-1层：

对于满足容许条件(12)的两个簇t和s，根据(15)式集合其中所有指标j∈t和i∈s,即可构造出具有如下低秩分解形式的容许块

其中，为配置矩阵，为聚合矩阵，为转移矩阵，下标#t、#s和r表示矩阵的维度。这里，符号“#”用来表示簇里的指标个数。对于中的容许块，#t表示观察组(簇)t中RWG基函数的个数，#s则表示源组(簇)s中RWG基函数的个数；r＝3×(P+1)²，其中(P+1)²表示L_i的总数，常数项“3”描述的是电流矢量的三个方向。对于中的容许块，#t和#s分别表示簇t和s中冲击基函数的个数，由于电荷是标量，因此有r＝(P+1)²。通过调节多极子的数目P，可以控制低秩分解的精度。

尽管LF-MLFMA构造的容许块已经具有(16)式所示的低秩分解矩阵形式，该低秩分解矩阵仍然存在冗余信息，可被进一步压缩，生成更为紧凑的叠层矩阵表达式。下面采用一种基于QR分解和奇异值分解的重压缩方法对容许块进一步压缩，具体过程如下：

1.对执行QR分解，生成矩阵Q_t和R_t：(这里符号表示复数域，上标#t和r表示矩阵维度)；

2.对执行QR分解，生成矩阵Q_s和R_s：(这里符号表示复数域，上标#s和r表示矩阵维度)；

3.执行的操作，并对进行奇异值分解，生成矩阵U，∑和V：其中∑为对角阵；

4.从矩阵∑的对角元素中提取部分元素构造出矩阵这里∑₁₁,∑₂₂…,∑_kk分别表示矩阵∑的第1,2…,k个对角元素，且∑_(k+1)(k+1)≤ε_rec∑₁₁＜∑_kk，其中，参数ε_rec表示相对截断误差，用来控制重压缩的精度；

5.分别从U和V中提取前k列元素获得矩阵和即：其中，U_k和V_k分别表示U和V的第k列元素；

6.获得其中，和

对于的构造，由于是稀疏矩阵，因此只需将中的非零元素填入电荷-电流块树中。中的非零元素均可被填入非容许块中，而不会被填入容许块中，原因在于，满足容许条件(12)的两个簇t和s必然离开一定的距离，而根据(7)式，离开一定距离的i∈t,j∈s生成的矩阵元素因此，可由无损地表达。只是的转置。此外，通过将电荷-电荷块树的主对角块的主对角线元素赋值为1即可构造出来。

对于的构造，只需将中的非零元素填入电荷-电荷块树中，这里，表示提取向量中与电荷相关的部分向量，即：与的构造过程不同,中的容许块极易生成，这是因为天然地具有低秩分解矩阵的形式，换句话说，中的任一容许块均可直接从中提取，并具有秩k＝1的低秩分解矩阵形式：

表示向量中与簇t相关的部分向量，其维度为#t×1；表示向量中与簇s相关的部分向量，其维度为#s×1；下标#t和#s分别表示簇t和s中基函数的个数。需要指出的是，当分析目标包含多个独立物体时，k的取值等于独立物体的数目。中的非容许块通过直接计算可方便地获得。

至此，系统矩阵的五个子矩阵的叠层矩阵表达式和均构造完成，从而获得了系统矩阵的叠层矩阵表达式

需要说明的是，上述第(1)步到第(3)步构造过程针对的是矩阵收缩方法。如果采用缩减电荷自由度的方法，只需在第(1)步中针对目标的每个独立物体去掉一个面片自由度即可，而无需加入修正矩阵

第4步，采用叠层矩阵格式的运算法则对该叠层矩阵进行上下三角分解。

根据(9)式，可表示为一个2×2的分块矩阵形式：

其中，及这里，表示中与电荷相关的子矩阵。

首先，通过计算来获得这里，采用叠层矩阵格式的乘法取代传统矩阵乘法来完成的操作。此外，采用叠层矩阵格式的加法取代传统矩阵加法来完成的操作。然后，从(18)式所示分块矩阵的上下三角分解出发，

这里，和分别表示对进行上下三角分解后获得的矩阵块，通过层层递归执行如下步骤来完成叠层矩阵的上下三角分解：

1.对执行上下三角分解：获得和

2.求解下三角方程：获得

3.求解上三角方程：获得

4.计算然后对其执行上下三角分解：获得和

上述步骤中涉及到的矩阵加法和乘法均采用叠层矩阵格式的加法和乘法来完成。叠层矩阵格式的加法和乘法中涉及低秩分解矩阵的运算，采用基于QR分解和奇异值分解的截断算子来定义，即：和引入自适应截断误差ε_t来控制截断算子的精度，ε_t的含义是用低秩分解矩阵与满阵之间的误差。叠层矩阵上下三角分解的精度也由ε_t来调控。叠层矩阵上下三角分解在低频电磁仿真时具有O(Nlog²N)计算复杂度和O(NlogN)内存消耗，其中，N为未知量的数目。完成叠层矩阵的上下三角分解之后，获得叠层矩阵格式的上下三角因子，结合叠层矩阵格式的前后向回代，即可用于完成对增量型电场积分方程的系统矩阵方程的求解。叠层矩阵格式的前后向回代在低频电磁仿真时的计算复杂度为O(NlogN)。

第5步，将叠层矩阵的上下三角分解与前后向回代结合，采用预条件的迭代解法或者直接解法进行求解。将系统矩阵方程简写作：

其中，表示待求的未知系数向量，表示右边向量。

采用预条件的迭代解法时，将方程(20)左右两端同时乘以一个预条件矩阵让具有比更小的条件数。叠层矩阵的上下三角分解用作预条件时，由于叠层矩阵分解之后，和均已获得，只需执行叠层矩阵格式的前后向回代即可完成一次的矩阵-向量乘操作，参与到迭代求解方程(20)的每一次矩阵-向量乘操作中。

采用直接解法时，由于方程(20)可写为：

通过执行一次前后向回代，即可获得方程的解。

实际中，通过调节截断误差ε_t可以调控叠层矩阵分解的精度。高精度(达到工程需求的精度)的叠层矩阵分解可用作直接解法，低精度(未达到工程需求的精度)的叠层矩阵分解可用作迭代解法的预条件。叠层矩阵分解精度越高，预条件的效果越好，但同时构造时间和内存消耗也越大。因此，选用何种精度的叠层矩阵分解作为预条件，以及选用直接解法还是预条件的迭代解法，都需要考虑多种因素进行权衡。一般地，直接解法更适合于多右边向量问题的求解，迭代解法更适合于单右边向量问题的求解。

求解出电流和电荷系数后，即可获得目标表面的电流及电荷分布，通过进一步计算即可提取出所需的电磁特性参数，可广泛用于芯片封装、集成电路、电磁散射、天线设计及电磁兼容等领域中的低频电磁特性仿真。

实施例1

采用本发明所述方法对一个具有复杂外形的飞机目标的低频电磁散射特性进行分析。该飞机目标的电尺寸为0.16λ×0.11λ×0.04λ，其中，λ表示入射电磁波的波长，如图7所示。目标受10MHz的入射平面波照射。采用三角形单元离散，生成40704条棱边和27136个面片，最小棱边长度仅为0.910×10^-4λ。按照本发明所述步骤，构建基于叠层矩阵分解的预条件。采用低频多层快速多极子算法的层数取为5层，构造的叠层矩阵层数为18层，多极子个数P取为5，容许块重压缩精度ε_rec取为10^-3。采用叠层矩阵分解预条件，来加速基于广义最小余量法(GMRES)的迭代解法的收敛，迭代收敛误差设置为10^-6。叠层矩阵分解预条件中，分别测试截断精度ε_t为10^-2和10^-3两种情况。图8给出了迭代收敛曲线，并与常用的约束预条件进行了对比，可以看出，本发明所述方法能够显著提高迭代收敛速度。表1给出求解效率的对比，可以看出，约束预条件尽管构造时间短，但预条件效果不理想，收敛速度慢，因此求解时间长；而本发明所述方法则显著提高了求解效率。此外，从表1不难看出，对于单激励问题，低截断精度的预条件求解效率高；但是对于多激励问题(如：单站雷达散射截面分析等)，高截断精度的预条件甚至直接解法则更有效，因为叠层矩阵分解只需构造一次，而迭代过程则需进行多次。图9和图10表明，本发明所述叠层矩阵分解方法能够将计算复杂度和内存消耗分别降低到O(Nlog²N)和O(NlogN)。

表1几种预条件技术之间的性能比较

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替代，都应当视为属于本发明的保护范围。