矩量法在电磁散射中的应用

1、矩量法在电磁散射中的应用一矩量法在电磁散射问题中的应用电磁散射问题是电磁学中的一个重要研究领域，研究电磁波的散射机理以及计算其散射场强的大小与分布，具有十分重要的实际意义。矩量法作为一种有效的数值计算方法在其中有着广泛的应用。但作为一种计算方法它也有着自己的缺陷，为了解决这些问题，人们提出了各种方案，矩量法在这个过程中也获得了很大的发展。MoM(Method of Moments)原本是一种近似求解线性算子方程的方法，通过它可以将算子方程转化为一矩阵方程，进而通过求解此矩阵方程得到最终的近似解。MoM最早是由两位数学家L. V. Kantorovich和V. IKrylov提出的，后来由KKM

2、ei引入计算电磁学，最终被RF. Harryington在其著作计算电磁场中的矩量法中加以系统描述。利用矩量法求解电磁问题的主要优点是：它严格地计算了各个子系统间的互耦，而算法本身又从根本上保证了误差系统总体最小而不产生数值色散。如今MoM被广泛应用于计算电磁学中，虽然它不能处理电大尺寸目标的电磁问题，但基于MoM的各种加速方法仍受到极大重视，如多层快速多极子方法MLMFA等。电磁散射问题是电磁学中的一个重要研究领域，研究电磁波的散射机理以及计算其散射场强的大小与分布，具有十分重要的实际意义。在实际生活中，遇到的散射目标往往不仅具有复杂的几何形状，而且构成的材料也各不相同。因此对复杂目标的电磁

3、散射特性进行快速、高效的分析，具有重要的理论意义和实用价值。电磁散射问题只有在相对简单的情况下才可以用严格的解析法来求解，比如对极少数形状规则的物体。对于电大物体，可以用高频近似方法，例如几何光学法(GO)、物理光学法(PO)、几何绕射理论(GTD)、物理绕射理论(PTD)、一致性几何绕射理论(UTD)、复射线法(CT)等来求解散射场。反之，对于电小物体，可以用准静态场来进行分析。介乎这两者之间的物体，一般采用数值方法。目前分析电磁散射问题的数值方法主要有微分方程法和积分方程法。微分方程法有有限差分法(FDM)、时域有限差分法(FDTD)、频域有限差分法(FDFD)、时域平面波法(PWTD)、

4、时域多分辨分析法(MRTD)、有限元法(FEM)、传输线矩阵法(TLM)等，积分方程法有表面积分方程法(SIEM)、矩量法(MOM)、边界元法(BEM)、体积分方程法(VIEM)、快速多极子法(FMM)、时域积分方程法(IETD)等。这些方法各有优缺点，有的是为了避免矩阵求逆，有的是为了加快收敛，有的是为了提高精度，还有的是为了减少贮存等。它们被广泛应用于求解复杂的工程电磁场问题中。 应用微分方程法求解电磁散射问题时，由于散射体的外空间为无限大，需要人为设置截断边界使求解区域有限，这种截断边界的引入会导致非物理的反射现象。矩量法是一种将连续方程离散化成代数方程组的方法，经常被看作数值“精确解”

5、。它既适用于电磁场积分方程又适用于微分方程，由于已经有有效的数值计算方法求解微分方程，所以目前矩量法大都用来求解积分方程。由于此方法能解决边界比较复杂的一些问题，因而得到了广泛的应用。需要注意的是，虽然矩量法中求解阻抗矩阵的表达式较为简单，但其计算工作量很大，对于以积分方程为基础的离散方程，其系数矩阵通常为满矩阵，所有元素都需大量的数值计算。尤其随着目标电尺寸的增大，矩量法得到的系数矩阵将迅速增大，这将给计算机内存和CPU带来沉重的负担。为了克服这些困难，人们对传统矩量法进行了一些改进，提出了一些快速、有效的方法，如快速多极子方法(FMM)和多层快速多极子方法(MLFMM)，阻抗矩阵局部化(I

6、ML)方法和压缩或稀疏化阻抗矩阵的小波分解法，(基于快速Fourier变换的CGFFT法、稀疏矩阵规则网格(SMCG)法和自适应积分法(AIM)，来降低计算机内存和计算量的需求。在这些快速分析方法中，需要的计算量和内存分别降为D(NlogW)和0(N)，N为未知变量数。但是，这些改进的方法仍然受传统矩量法离散尺寸的限制。采用整域基函数代替分域基函数可以降低矩量法系数矩阵的维数。然而，在绝大多数情况下，难以找到合适的整域基函数。为此，近年来，人们又相继开展了一些基于部分域(子域或块)概念来降低矩阵维数的研究，如多层矩量法(MMM)、子域多层法(SMA)、合成基函数(SBF)法等。这些方法通过对问

7、题的部分域进行分析来构造宏基函数，宏基函数的域比传统矩量法的大，因此可以降低未知变量数。这几种方法是通过迭代的方式递归地修正互耦项来改进解的收敛性。特征基函数法(CBFM)是近两年提出来的一种新的求解电磁散射问题的有效方法。二矩量法的基本原理矩量法是将算子方程转化为矩阵方程，然后求解该矩阵方程的方法，由于在求解方程过程中，需要计算广义矩量，所以我们称这种方法称为矩量法，其实质是内域基加权余量法。(一)离散化过程已知算子方程：程已知，唯一确定。设（是系数，是基函数）。其中不一定是完备的，越大，不一定越好，只有完备，越大越接近真值。方程离散为： 此过程的主要目的在于将算子方程转为代数方程，具体步骤

8、如下：(1) 在算子的定义域内适当的选择一组基函数，他们应该是线性无关的。(2) 将未知函数表示为该组基函数的线性组合，并取有限项近似，即: (2-1)(3) 将公式（2-1）利用算子的线性，可将算子方程化为代数方程，即: (2-2)于是，待求的问题化解为求解的系数的问题。(二)取样检验过程为了使的近似函数之间的误差极小，必须进行抽样检验，在抽样点上使加权平均误差为零，从而确定未知系数，这一过程的基本步骤如下：(1) 在算子的值域内取适当的选择一组加权函数，他们也是线性无关的。(2) 将与式（2-2）取内积进行抽样检验，因为要确定个未知数，需要进行次抽样检验，则 （2-3） (3) 利用算子的

9、线性和内积的性质，将（2-3）式转化为矩阵方程，即: （选配过程） （2-4）将它写成矩阵形式： (2-5)于是，求解代数方程问题转化为求解矩阵方程的问题。(三)矩阵求逆一旦求得了矩阵方程，通过常规的矩阵求逆或求解线性方程组，就可以得到矩阵方程的解。如果是非奇异的，则它的逆矩阵 是存在的，则矩阵方程的解为： （2-6）将求得的展开系数带入到（2-1）式中，就可以得到原来算子方程的近似解为： (2-7) 以上所述是矩量法求解算子方程的基本过程，在矩量法的所有应用中，通常都要遵循这个统一的过程。由于和之间存在近似，故算子方程的左端近似值与右端精确值之间存在如下关系： （2-8）三基函数与权函数的选

10、择矩量法的求解过程是简单的，求解步骤是统一的，应用起来比较方便。求解精度取决于很多因素，例如离散化的成都、基函数和权函数的选取和矩阵的求解过程等，其中基函数与权函数的选择尤其重要，可以直接影响到求解的精度。基函数可以分为全域基和分域基，权函数可以分为权函数、分域权和点匹配，它们之间的不同组合便形成不同的方法。下面简单的对此进行介绍。1、 全域基函数：全域基函数是指算子定义域内的全域上存在一组基函数。他们应该满足边界条件。在矩量法求解的离散化过程中，选择全域基函数作为展开函数，实质上将未知函数表示为全域存在的若干个离散化基函数的线性组合。2、 分域基函数：分域基函数不是在算子定义域的全域上存在的

11、，而仅仅存在于算子定义域的各个分域上的函数。选择分域基函数作为未知函数的展开函数，在矩量法求解的离散化过程中是一种区域离散，即未知函数表示为各个分域上存在的函数线性组合。3、 全域权：在算子的值域内选择一组权函数，它是一组在值域内的全域上存在的权函数。在此方法中，如果选择，即权函数等于基函数，则称此方法为伽略金法。4、 点匹配:全域基、全域权的伽略金发是一种常用的求解方法。但是，如果算子本身是复杂的积分算子，或者选择了比较复杂的基函数，由于内积运算本身又是积分运算，从而使矩阵元素和的形成十分困难。这时为了简化，我们可以利用函数的选择性将函数作为权函数将使内积计算得以简化。四算例问题：计算带电导

12、体板的电容。正方形板，边长为，位于。空间中任意一点电位，其中。在板上为常数，有，其中电荷分布为待求，电容为：。计算当时方板的归一化电容。1 选取分域基：方板各边等分，其中个子域。定义：，其中在上为1，否则为0 。将带入中，有 ， 其中。是上单位振幅的均匀电荷密度在的中心处产生的电位。因此有 物体的电容式其各小块电容的总和加上每一对小块间的互电容。 2 取 3 的计算：由本身面上的单位电荷密度在其中心处产生的电位（自作用单元）而互作用单元：经过计算模拟可以发现当取值越大时，结果将趋近于一个定值。当为20时，计算得到归一化电容为40.29PF，非常接近理论计算值40PF，误差为0.725%，在误差

13、允许范围内。程序：#pragma onceclass Martixprivate:const int M;const int N;bool flag;double *array;public:Martix(double *_array, int _M, int _N) :M(_M), N(_N)array = _array;flag = false;Martix(int _M, int _N) :M(_M), N(_N)flag = true;array = new double *M;for (int i = 0; i < M; i+)arrayi = new doubleN;for

14、(int i = 0; i < M; i+)for (int j = 0; j < N; j+)arrayij = 0;Martix Gauss();void show();void swap(int i, int j);int Max(int i);Martix operator+(Martix &T);Martix operator-(Martix &T);Martix operator*(Martix &T);Martix operator*(double t);friend Martix operator*(double t, const Marti

15、x &T)Martix sum(T.M, T.N);for (int i = 0; i < T.M; i+)for (int j = 0; j < T.N; j+)sum.arrayij = T.arrayij * t;return sum;Martix Copy();double *Getarray()return this->array;int GetM()return this->M;int GetN()return this->N;bool GetFlag()return this->flag;Martix();#include "

16、Martix.h"#include<iomanip>#include<iostream>using namespace std;void Martix:show()for (int i = 0; i < M; i+)cout << fixed << right;for (int j = 0; j < N; j+)cout << setw(10) << setprecision(6) << arrayij << " "cout << endl;Ma

17、rtix Martix:Gauss()Martix sum(M, N);for (int i = 0; i < N; i+)sum.arrayii = 1;for (int i = 0; i < N; i+)if (abs(arrayii) <= 0.00001&&i<M - 1)if (abs(arrayMax(i)i) <= 0.00001)continue;int count = Max(i);swap(i, count);sum.swap(i, count);else if (abs(arrayii) < 0.00001&&a

18、mp;i = N - 1)cout << "没有逆矩阵" << endl;exit(0);/*cout << "变换前" << endl;this->show();*/double temp = arrayii;for (int j = 0; j < N; j+)arrayij = arrayij / temp;sum.arrayij = sum.arrayij / temp;for (int k = i + 1; k < N; k+)double temp = arrayki;if (a

19、bs(temp) >= 0.00001)for (int j = 0; j < N; j+)arraykj = arraykj - temp*arrayij;sum.arraykj = sum.arraykj - temp*sum.arrayij;double t = abs(array00);for (int i = 1; i < N; i+)t = t*abs(this->arrayii);if (abs(t) <= 0.0001)cout << "没有逆矩阵" << endl;exit(0);for (int i

20、= N - 2; i >= 0; i-)for (int j = 0; j < N; j+)for (int k = i + 1; k < N; k+)sum.arrayij = sum.arrayij - arrayik * sum.arraykj;return sum;void Martix:swap(int i, int j)double temp = 0;for (int k = 0; k < N; k+)temp = arrayik;arrayik = arrayjk;arrayjk = temp;Martix Martix:Copy()Martix sum(

21、this->M, this->N);for (int i = 0; i < M; i+)for (int j = 0; j < N; j+)sum.arrayij = this->arrayij;return sum;Martix Martix:operator+(Martix &T)Martix sum(T.M, T.N);if (M != T.M | N != T.N)cout << "错误" << endl;for (int i = 0; i < M; i+)for (int j = 0; j <

22、; N; j+)sum.arrayij = this->arrayij + T.arrayij;return sum;Martix Martix:operator-(Martix &T)Martix sum(T.M, T.N);if (M != T.M | N != T.N)cout << "错误" << endl;for (int i = 0; i < M; i+)for (int j = 0; j < N; j+)sum.arrayij = this->arrayij - T.arrayij;return sum;

23、Martix Martix:operator*(Martix &T)Martix sum(T.M, T.N);if (N != T.M)cout << "错误" << endl;for (int i = 0; i < M; i+)for (int j = 0; j < T.N; j+)for (int k = 0; k < N; k+)sum.arrayij = sum.arrayij + this->arrayik * T.arraykj;return sum;Martix Martix:operator*(doub

24、le t)Martix sum(this->M, this->N);for (int i = 0; i < this->M; i+)for (int j = 0; j < this->N; j+)sum.arrayij = this->arrayij * t;return sum;int Martix:Max(int i)int count = i;double temp = abs(arrayii);for (int k = i + 1; k < N; k+)if (temp < abs(arrayki)temp = arrayki;wh

25、ile (temp != abs(arraycounti)count+;return count;Martix:Martix()#include<iostream>#include<cmath>#include<iomanip>#include"Martix.h"using namespace std;const double pi = 3.1415926;const double eps = 1 / (36 * pi)*pow(10, -9);double a = 2.0;int N = 40;double b = a / N;int M = N*N;double ds = 4*b*b;void clean(Martix &T)for (int i = 0; i < T