CN106066941A - 一种基于cbfm和smw算法的电磁散射快速分析方法 - Google Patents

Abstract

本发明公开一种基于特征基函数法(CBFM)和SMW算法的快速分析目标电磁散射的方法，该方法能快速、直接地计算电大目标的电磁散射特性。首先，利用二叉树结构对离散目标的基函数进行分块，在最底层的每一块生成特征基函数(CBF)；然后使用SMW算法将CBFM中的缩减阻抗矩阵近似表示成一系列块对角矩阵的乘积，由于每一个块对角矩阵都具有特殊的结构，所以这些块对角矩阵的逆矩阵可以用SMW公式快速求出；最后求得目标的电磁散射特性。该方法用特征基函数法(CBFM)显著减少了未知量的数目，且使用基于SMW公式的直接求解算法进一步加速了计算速度。该方法不仅易于实现和适合并行计算，并且具有核独立特性，因此适合处理大部分的电磁散射问题。

Description

一种基于CBFM和SMW算法的电磁散射快速分析方法

技术领域：

本发明涉及一种快速直接分析电大尺寸导体目标电磁散射的方法，尤其涉及一种基于特征基函数法(CBFM)和SMW算法的直接求解目标电磁散射特性的混合方法。

背景技术：

矩量法(Method of Moments,MoM)将电磁积分方程转化成矩阵方程，来求解电磁散射问题。而电大目标的电磁散射一直受到国内外学者的广泛关注。有两种方法来求解矩阵方程：迭代求解和直接求解。相比而言，直接求解没有收敛问题且可以求解单站RCS。然而，对于电大目标，传统的直接求解方法比如LU分解，还是很消耗内存和时间。

特征基函数法(Characteristic Basis Function Method,CBFM)可以通过在一定区域定义特征基函数(CBF)，实现未知量的个数的降低，从而降低所占内存和消耗时间。然而随着目标电尺寸的平方增长，特征基函数法中用来直接求解的缩减阻抗矩阵还是非常大。

一种解决方法是多层特征基函数法(Multilevel CBFM,MLCBFM)，其通过递归地使用单层CBFM实现更大的压缩率。本发明提出另外一种基于特征基函数(CharacteristicBasis Function,CBF)和Sherman-Morrison-Woodbury(SMW)公式的方法来加快直接求解速度和节省内存。其中SMW算法将特征基函数法压缩后的缩减阻抗矩阵转化为一系列块对角矩阵相乘的形式，这些块对角矩阵具有特殊的结构，可以通过SMW公式快速得到逆矩阵。值得一提的是其它的一些快速直接求解算法比如多尺度压缩快分解(MSCBD)和多尺度LU分解(MSLU)同样可以加快缩减矩阵的直接求解速度。SMW算法主要的优势是易于实现，而且适合并行计算。

本发明公开一种将CBFM和SMW算法优点结合一起的快速直接求解方法。

发明内容：

发明目的：本发明解决的是快速分析导体目标电磁散射的问题，本发明提出了一种基于CBFM和SMW算法的电磁散射快速分析方法。该方法在利用二叉树结构对目标进行分块之后，在最底层生成特征基函数，降低未知量的数目，将阻抗矩阵转化为缩减阻抗矩阵；接着用SMW算法快速求解缩减矩阵的逆。该方法能够在阻抗矩阵压缩后再进一步降低求逆矩阵的复杂度，从而进一步加快求解速度。

为了达到上述目的，本发明的技术方案实现的基本步骤如下：

第1步：对导体目标的表面用三角形面片进行离散，在剖分得到的三角形网格的公共边上构造Rao-Wilton-Glisso(RWG)基函数；根据导体目标表面边界条件建立表面积分方程，用定义的RWG基函数对表面积分方程进行离散；

第2步：根据基函数的质心的位置，用L层二叉树将基函数进行多层分块，L为正整数；

第3步：在最底层的每块上生成特征基函数；

第4步：根据多层分块，将缩减阻抗矩阵的非主对角块对用自适应交叉近似(ACA)算法压缩；

第5步：利用SMW算法将缩减阻抗矩阵表示成L+1个块对角矩阵相乘的形式，并求解最终的电流系数；

第6步：根据第5步所得结果，解出目标远场雷达散射截面RCS的值。

与现有技术相比，本发明的优势在于：利用特征基函数法，在二叉树分组后的底层生成特征基函数(CBF)，并在此基础上利用SMW算法，有效的进一步加快了求解速度。该方法在利用特征基函数法有效降低阻抗矩阵维数的情况下，又进一步利用SMW直接求解算法，将求解一个大矩阵的逆转化为求解多个块对角矩阵的逆，从而减少了计算时间和内存需求，有效地提升了电磁散射问题的直接求解效率。

附图说明：

图1是本发明基函数二叉树分块示意图。

图2是本发明三层SMW算法阻抗矩阵分块示意图。

图3是本发明三层SMW算法阻抗矩阵近似转化对角块矩阵乘积的示意图。

图4是本发明方法的算例结果示意图。

具体实施方案：

下面结合附图对技术方案的实施作进一步的详细描述：

第1步：针对导体目标的表面用三角形面片进行离散，三角形面片的平均边长为0.1λ，λ为入射平面波的波长。接着在导体目标表面建立的表面积分方程为，

$\hat{n}\left(r\right)\×\hat{n}\left(r\right)\×\left({\mathrm{j\ω\μ}}_0{\∫}_{s}G\right(r,r^{\′}\left)J\right(r^{\′}){\mathrm{ds}}^{\′}-\▿\frac{1}{{\mathrm{j\ω\ϵ}}_0}{\∫}_{s}G(r,r^{\′}){\▿}_s^{\′}\·J(r^{\′}\left){\mathrm{ds}}^{\′}\right)=\hat{n}\left(r\right)\×\hat{n}\left(r\right)\×E^i\left(r\right)---\left(1\right)$

其中，为单位矢量，j为虚数单位，μ₀为自由空间磁导率，ε₀为自由空间电导率，为梯度算子，为散度算子，J(r′)为导体目标上任一源点r′处的感应电流，G(r,r′)＝e^-jk|r-r′|/(4π|r-r′|)为自由空间格林函数，Eⁱ(r)是入射平面波电场，r是任一场点位置矢量，r′是任一源点位置矢量。

用定义的RWG基函数对表面积分方程进行离散，导体表面上的感应电流可以近似展开为

$J\left(r\right)\≈\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{I}_n{f}_n\left(r\right)---\left(2\right)$

其中，f_n(r)表示第n个RWG基函数，I_n表示与第n个RWG基函数所对应的系数，N是RWG基函数的数目；

第2步：根据基函数的质心的位置，用L层二叉树将基函数进行多层分块，L为正整数。如附图1所示，第1层划分，根据基函数X轴坐标将其分为两组，且每组基函数的数目相同；第2层划分，根据基函数Y轴坐标将第一次划分得到的两个子组分别再均匀分成两组，这样就得到四组基函数；第3层划分，根据基函数Z轴坐标进行。第4层与第1层相似，按照X轴坐标划分，依次类推，经过L层划分后，整个目标的基函数被分为2^L组，且每个组含有相同数目的基函数。第L层包含2^L个块。

第3步：在最底层的每块上生成激励无关的特征基函数，这些特征基函数通过平面波普的照射并结合奇异值分解提取得到。第i块第n个特征基函数可以表示为

$F_{i,n}\left(r\right)=\underset{p=1}{\overset{N_i}{\mathrm{\Σ}}}{J}_i(p,n)f_{i,p}\left(r\right)---\left(3\right)$

其中，f_i,p(r)表示第i块的第p个RWG基函数，N_i表示第i块中包含的RWG基函数的数目，J_i(p,n)是联系第i个分块第n个CBF和第i个分块第p个RWG基函数的线性表出系数，每一列对应一个CBF。假设第i块共生产K_i个特征基函数(CBF)。

第4步：首先，根据最底层(第L层)的分组情况，整个缩减阻抗矩阵方程可以表示为

其中，是第i个和第j个分块的CBF之间的阻抗矩阵，即缩减的子阻抗矩阵，

$Z_{ij}^{CBF}=J_i^T{Z}_{ij}{J}_j---\left(5\right)$

这里，和分别是第i个和第j个分块的特征基函数矩阵，是J_i的转置。是第i个和第j个分块中RWG基函数之间的阻抗矩阵。N_i和N_j分别表示第i块和第j块包含的RWG基函数数目，K_i和K_j分别表示第i块和第j块包含的CBF数目。B＝2^L。为第i块所有特征基函数对应的电压向量

$E_i^{CBF}=J_i^T{E}_i^{RWG}---\left(6\right)$

其中，为第i块所有RWG基函数对应的电压向量。为第i块特征基函数对应的待求电流系数。

然后，根据多层二叉树分块情况，缩减矩阵可以划分成如附图2所示的结构，这里附图2给出的是根据3层二叉树结构划分的结果。然后利用自适应交叉近似(ACA)算法对各层中具有低秩特性的矩阵块进行压缩分解，这些具有低秩特性的矩阵块就是附图2中那些非对角块(互阻抗矩阵)。

第5步：如附图3所示，利用SMW算法将缩减阻抗矩阵Z近似表示成L+1个块对角矩阵相乘的形式，

Z≈Z_LZ_L-1…Z₀ (7)

其中，Z_l(l＝0,1,2,…,L)为在第l层的块对角矩阵。并利用(8)式求出所需的电流系数，

$I\≈Z_0^{-1}{Z}_1^{-1}...Z_{L-2}^{-1}{Z}_{L-1}^{-1}{Z}_L^{-1}E---\left(8\right)$

其中，E为所有CBF对应的电压向量。

SMW算法的基本原理是：

SMW算法是一个多层算法，为了更清楚说明其原理，先描述1层的SMW算法。在1层算法中，缩减矩阵方程按照第1层的分块情况可以划分成如下形式

$\left(\begin{array}{cc}{Z}_{11}& Z_{12}\\ Z_{21}& Z_{22}\end{array}\right)\left(\begin{array}{c}{I}_1\\ I_2\end{array}\right)=\left(\begin{array}{c}{E}_1\\ E_2\end{array}\right)---\left(9\right)$

Z₁₁和Z₂₂分别为第1个和第2个分块的所有CBF之间的阻抗矩阵，Z₁₂和Z₂₁为第1块和第2块CBF之间的阻抗矩阵，这四个矩阵由最底层的缩减矩阵聚集而成。E₁和E₂为第1块和第2块所有CBF对应的电压向量，I₁和I₂为第1块和第2块CBF对应的待求电流系数。

在第4步中Z₁₂和Z₂₁已经被ACA算法压缩为

Z₁₂≈U₁₂V₁₂ (10)

Z₂₁≈U₂₁V₂₁ (11)

其中，U₁₂和U₂₁的大小为(N^CBF/2)×r,V₁₂和V₂₁的大小为r×(N^CBF/2)，r是Z₁₂和Z₂₁的有效秩，其远小于N^CBF。

将(10)式和(11)式代入(9)式且同时左乘[Z₁₁,0；0,Z₂₂]的逆，

${\left[\begin{array}{cc}{Z}_{11}& 0\\ 0& Z_{22}\end{array}\right]}^{-1}\left[\begin{array}{cc}{Z}_{11}& U_{12}{V}_{12}\\ U_{21}{V}_{21}& Z_{22}\end{array}\right]\left[\begin{array}{c}{I}_1\\ I_2\end{array}\right]={\left[\begin{array}{cc}{Z}_{11}& 0\\ 0& Z_{22}\end{array}\right]}^{-1}\left[\begin{array}{c}{E}_1\\ E_2\end{array}\right]---\left(12\right)$

其中，0表示零矩阵。之后得到

$\left[\begin{array}{cc}1& {\hat{U}}_{12}{V}_{12}\\ {\hat{U}}_{21}{V}_{21}& 1\end{array}\right]\left[\begin{array}{c}{I}_1\\ I_2\end{array}\right]=\left[\begin{array}{c}{\hat{E}}_1\\ {\hat{E}}_2\end{array}\right]---\left(13\right)$

其中，1表示单位矩阵，

${\hat{U}}_{12}=Z_{11}^{-1}{U}_{12},{\hat{E}}_1=Z_{11}^{-1}{E}_1,---\left(14\right)$

${\hat{U}}_{21}=Z_{22}^{-1}{U}_{21},{\hat{E}}_2=Z_{22}^{-1}{E}_2.---\left(15\right)$

矩阵可以写成，

$\left[\begin{array}{cc}1& {\hat{U}}_{12}{V}_{12}\\ {\hat{U}}_{21}{V}_{21}& 1\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]+\left[\begin{array}{cc}0& {\hat{U}}_{12}\\ {\hat{U}}_{21}& 0\end{array}\right]\left[\begin{array}{cc}{V}_{21}& 0\\ 0& V_{12}\end{array}\right]---\left(16\right)$

根据Sherman-Morrison-Woodbury公式，(16)式的逆矩阵等于

${\left[\begin{array}{cc}1& {\hat{U}}_{12}{V}_{12}\\ {\hat{U}}_{21}{V}_{21}& 1\end{array}\right]}^{-1}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]-\left[\begin{array}{cc}0& {\hat{U}}_{12}\\ {\hat{U}}_{21}& 0\end{array}\right]{\left[\begin{array}{cc}1& V_{21}{\hat{U}}_{12}\\ V_{12}{\hat{U}}_{21}& 1\end{array}\right]}^{-1}\left[\begin{array}{cc}{V}_{21}& 0\\ 0& V_{12}\end{array}\right]---\left(17\right)$

将(17)式代入(13)式，得出电流系数为

$\left[\begin{array}{c}{I}_1\\ I_2\end{array}\right]=\left[\begin{array}{c}{\hat{E}}_1\\ {\hat{E}}_2\end{array}\right]-\left[\begin{array}{cc}0& {\hat{U}}_{12}\\ {\hat{U}}_{21}& 0\end{array}\right]{\left[\begin{array}{cc}1& V_{21}{\hat{U}}_{12}\\ V_{12}{\hat{U}}_{21}& 1\end{array}\right]}^{-1}\left[\begin{array}{c}{V}_{21}{\hat{E}}_1\\ V_{12}{\hat{E}}_2\end{array}\right]---\left(18\right)$

可以看出，使用SMW算法，我们只需要计算一个2r×2r矩阵的逆，而不是原先的N^CBF×N^CBF矩阵的逆。

(18)式中需要求解以下等式，

$\left[\begin{array}{cc}1& P_{12}\\ P_{21}& 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]---\left(19\right)$

其中

$P_{12}=V_{21}{\hat{U}}_{12},y_1=V_{21}{\hat{E}}_1,---\left(20\right)$

$P_{21}=V_{12}{\hat{U}}_{21},y_2=V_{12}{\hat{E}}_2---\left(21\right)$

通过LU分解，(19)式转化为

$\left[\begin{array}{cc}1& 0\\ P_{21}& 1\end{array}\right]\left[\begin{array}{cc}1& P_{12}\\ 0& 1-P_{21}{P}_{12}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]---\left(22\right)$

从而(19)式的解为

x₂＝(1-P₂₁P₁₂)^-1(y₂-P₂₁y₁) (23)

x₁＝y₁-P₁₂x₂ (24)

将(23)式和(24)式代入(18)式，可以得出待求的电流系数解为

$I_1={\hat{E}}_1-{\hat{U}}_{12}{x}_2---\left(25\right)$

$I_2={\hat{E}}_2-{\hat{U}}_{21}{x}_1---\left(26\right)$

将1层SMW算法推广到多层SMW算法，其步骤如下：

步骤1：生成底层块对角矩阵

其中，为第i块中所有CBF的自阻抗矩阵，它的逆可以通过LU分解得到。那么Z_L的逆由所有的的逆组合而成。

步骤2：取在第L-1层的主对角矩阵块组成矩阵Z_L-1。其中Z_L-1的每一个对角矩阵块有与(16)式相同的形式，因而其逆矩阵可以用1层SMW算法快速求出。

步骤3：按步骤2的方法，依次求出Z_L-2至Z₀。其中Z_l(l＝L-1,L-2,…,0)是在第l层的块对角矩阵。其中所需的逆矩阵可以通过1层SMW算法快速求出。

步骤4：利用求得的Z_L的逆和Z_L-1、Z_L-2…Z₀根据(8)式求出电流系数I。

第6步：利用感应电流系数I解出剩余目标的远场雷达散射截面RCS，表示为：

$RCS=\underset{r\→\mathrm{\∞}}{\lim }4{\mathrm{\πr}}^2\frac{{\left|E^s\right|}^2}{{\left|E^i\right|}^2}---\left(28\right)$

其中，E^s为远场散射场，Eⁱ为入射场。

下面以一具体实例对本发明方法作进一步说明：

本发明以一个金属球的散射问题为研究对象加以详细论述，金属球的半径为3m。入射平面波的工作频率为300MHz，入射方向为方向，入射波的电场方向为下面按照技术方案的过程实现该目标电磁散射问题的高效求解。

首先根据第1～3步，对金属球进行离散得到40851个RWG基函数。用8层二叉树对目标进行分块，最底层每块平均包含159个RWG基函数。在最底层每一块中生成特征基函数，共得到10965个特征基函数。

然后根据第4～6步，使用SMW算法快速计算出缩减阻抗矩阵的逆矩阵，计算导体目标的感应电流系数，最终求出导体目标的RCS。

最终解出金属球的远场雷达散射截面如附图4。从附图4可以看出，用本发明提出的方法对金属球的计算结果与Mie级数解析解吻合良好。然后我们将本发明中的方法与传统的CBFM-ACA进行计算效率和内容需求的比较。使用CBFM-ACA求解的时间和内存分别为1001s和931MB，而使用本发明提出的方法只需要323s的时间和395MB的内存。在传统的CBFM-ACA方法中，使用的传统的LU分解求解缩减矩阵方程。如果使用传统的矩量法(MoM)求解该目标，预计需要12.4GB的内存。可以看出本发明提出的快速算法的效率显著比传统的CBFM-ACA和MoM高。值得说明的是，当计算电尺寸更大的导体目标时，本发明实现的时间和内存的缩减会变得更加明显。

Claims (4)

1.一种基于CBFM和SMW算法的电磁散射快速分析方法，其特征在于，步骤如下：

第1步：对导体目标的表面用三角形面片进行离散，在剖分得到的三角形网格的公共边上构造RWG基函数；根据导体目标表面边界条件建立表面积分方程，用定义的RWG基函数对表面积分方程进行离散；

第2步：根据基函数的质心的位置，用L层二叉树将基函数进行多层分块，L为正整数；

第3步：在最底层的每块上生成特征基函数；

第4步：根据多层分块，将缩减阻抗矩阵的非主对角块对用自适应交叉近似算法压缩；

第5步：利用SMW算法将缩减阻抗矩阵表示成L+1个块对角矩阵相乘的形式，并求解最终的电流系数；

第6步：根据第5步所得结果，解出目标远场雷达散射截面RCS的值。

2.根据权利要求1所述的一种基于CBFM和SMW算法的电磁散射快速分析方法，其特征在于，第2步中，根据基函数的质心的位置，用L层二叉树将基函数进行多层分块，L为正整数；第1层划分，根据基函数X轴坐标将其分为两组，且每组基函数的数目相同；第2层划分，根据基函数Y轴坐标将第1次划分得到的两个子组分别再均匀分成两组，这样就得到四组基函数；第3层划分，根据基函数Z轴坐标进行；第4层与第1层相似，按照X轴坐标划分，依次类推，经过L层划分后，整个目标的基函数被分为2^L组，且每个组含有相同数目的基函数；第L层包含2^L个块。

3.根据权利要求1所述的一种基于CBFM和SMW算法的电磁散射快速分析方法，其特征在于，第4步中，首先，根据最底层的分组情况，整个缩减阻抗矩阵方程可以表示为

其中，是第i个和第j个分块的CBF之间的阻抗矩阵，即缩减的子阻抗矩阵，

$Z_{ij}^{CBF}=J_i^T{Z}_{ij}{J}_j---\left(5\right)$

这里，和分别是第i个和第j个分块的特征基函数矩阵，是J_i的转置；是第i个和第j个分块中RWG基函数之间的阻抗矩阵；N_i和N_j分别表示第i块和第j块包含的RWG基函数数目，K_i和K_j分别表示第i块和第j块包含的CBF数目，B＝2^L；为第i块所有特征基函数对应的电压向量

$E_i^{CBF}=J_i^T{E}_i^{RWG}---\left(6\right)$

其中，为第i块所有RWG基函数对应的电压向量，为第i块特征基函数对应的待求电流系数；

然后利用自适应交叉近似(ACA)算法对各层中具有低秩特性的矩阵块进行压缩分解。

4.根据权利要求1所述的一种基于CBFM和SMW算法的电磁散射快速分析方法，其特征在于，第5步，利用SMW算法将缩减阻抗矩阵Z近似表示成L+1个块对角矩阵相乘的形式，

Z≈Z_LZ_L-1…Z₀ (7)

其中，Z_l(l＝0,1,2,…,L)为在第l层的块对角矩阵；并利用(8)式求出所需的电流系数，

$I\≈Z_0^{-1}{Z}_1^{-1}...Z_{L-2}^{-1}{Z}_{L-1}^{-1}{Z}_L^{-1}E---\left(8\right)$

其中，E为所有CBF对应的电压向量；

其步骤如下：

步骤1：生成底层块对角矩阵

其中，为第i块中所有CBF的自阻抗矩阵，它的逆可以通过LU分解得到；那么Z_L的逆由所有的的逆组合而成；

步骤2：取在第L-1层的主对角矩阵块组成矩阵Z_L-1；其中Z_L-1的每一个对角矩阵块有与(16)式相同的形式，因而其逆矩阵可以用1层SMW算法快速求出；

步骤3：按步骤2的方法，依次求出Z_L-2至Z₀；其中Z_l(l＝L-1,L-2,…,0)是在第l层的块对角矩阵；其中所需的逆矩阵可以通过1层SMW算法快速求出；

步骤4：利用求得的Z_L的逆和Z_L-1、Z_L-2…Z₀根据(8)式求出电流系数。