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

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

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

电磁散射问题是电磁学中的一个重要研究领域，研究电磁波的散射机理以及计算其散射场强的大小与分布，具有十分重要的实际意义。矩量法作为一种有效的数值计算方法在其中有着广泛的应用。但作为一种计算方法它也有着自己的缺陷，为了解决这些问题，人们提出了各种方案，矩量法在这个过程中也获得了很大的发展。

MoM(Method of Moments)原本是一种近似求解线性算子方程的方法，通过它可以将算子方程转化为一矩阵方程，进而通过求解此矩阵方程得到最终的近似解。MoM最早是由两位数学家L. V. Kantorovich和V. I．Krylov提出的，后来由K．K．Mei引入计算电磁学，最终被R．F. Harryington在其著作《计算电磁场中的矩量法》中加以系统描述。利用矩量法求解电磁问题的主要优点是：它严格地计算了各个子系统间的互耦，而算法本身又从根本上保证了误差系统总体最小而不产生数值色散。如今MoM被广泛应用于计算电磁学中，虽然它不能处理电大尺寸目标的电磁问题，但基于MoM的各种加速方法仍受到极大重视，如多层快速多极子方法MLMFA等。

电磁散射问题是电磁学中的一个重要研究领域，研究电磁波的散射机理以及计算其散射场强的大小与分布，具有十分重要的实际意义。在实际生活中，遇到的散射目标往往不仅具有复杂的几何形状，而且构成

的材料也各不相同。因此对复杂目标的电磁散射特性进行快速、高效的分析，具有重要的理论意义和实用价值。

电磁散射问题只有在相对简单的情况下才可以用严格的解析法来求解，比如对极少数形状规则的物体。对于电大物体，可以用高频近似方法，例如几何光学法(GO)、物理光学法(PO)、几何绕射理论(GTD)、物理绕射理论(PTD)、一致性几何绕射理论(UTD)、复射线法(CT)等来求解散射场。反之，对于电小物体，可以用准静态场来进行分析。介乎这两者之间的物体，一般采用数值方法。目前分析电磁散射问题的数值方法主要有微分方程法和积分方程法。微分方程法有有限差分法(FDM)、时域有限差分法(FDTD)、频域有限差分法(FDFD)、时域平面波法(PWTD)、时域多分辨分析法(MRTD)、有限元法(FEM)、传输线矩阵法(TLM)等，积分方程法有表面积分方程法(SIEM)、矩量法(MOM)、边界元法(BEM)、体积分方程法(VIEM)、快速多极子法(FMM)、时域积分方程法(IETD)等。这些方法各有优缺点，有的是为了避免矩阵求逆，有的是为了加快收敛，有的是为了提高精度，还有的是为了减少贮存等。它们被广泛应用于求解复杂的工程电磁场问题中。应用微分方程法求解电磁散射问题时，由于散射体的外空间为无限大，需要人为设置截断边界使求解区域有限，这种截断边界的引入会导致非物理的反射现象。矩量法是一种将连续方程离散化成代数方程组的方法，经常被看作数值“精确解”。它既适用于电磁场积分方程又适用于微分方程，由于已经有有效的数值计算方法求解微分方程，所以目前矩量法大都用来求解积分方程。由于此方法能解决边界比较复杂的一些问题，因而得到了广泛的应用。需要注意的是，

虽然矩量法中求解阻抗矩阵的表达式较为简单，但其计算工作量很大，对于以积分方程为基础的离散方程，其系数矩阵通常为满矩阵，所有元素都需大量的数值计算。尤其随着目标电尺寸的增大，矩量法得到的系数矩阵将迅速增大，这将给计算机内存和CPU带来沉重的负担。为了克服这些困难，人们对传统矩量法进行了一些改进，提出了一些快速、有效的方法，如快速多极子方法(FMM)和多层快速多极子方法(MLFMM)，阻抗矩阵局部化(IML)方法和压缩或稀疏化阻抗矩阵的小波分解法，(基于快速Fourier变换的CG—FFT法、稀疏矩阵规则网格(SMCG)法和自适应积分法(AIM)，来降低计算机内存和计算量的需求。在这些快速分析方法中，需要的计算量和内存分别降为D(NlogW)和0(N)，N为未知变量数。但是，这些改进的方法仍然受传统矩量法离散尺寸的限制。采用整域基函数代替分域基函数可以降低矩量法系数矩阵的维数。然而，在绝大多数情况下，难以找到合适的整域基函数。为此，近年来，人们又相继开展了一些基于部分域(子域或块)概念来降低矩阵维数的研究，如多层矩量法(MMM)、子域多层法(SMA)、合成基函数(SBF)法等。这些方法通过对问题的部分域进行分析来构造宏基函数，宏基函数的域比传统矩量法的大，因此可以降低未知变量数。这几种方法是通过迭代的方式递归地修正互耦项来改进解的收敛性。特征基函数法(CBFM)是近两年提出来的一种新的求解电磁散射问题的有效方法。

二矩量法的基本原理

矩量法是将算子方程转化为矩阵方程，然后求解该矩阵方程的方

法，由于在求解方程过程中，需要计算广义矩量，所以我们称这种方法称为矩量法，其实质是内域基加权余量法。

(一)离散化过程

已知算子方程：

Lf g =

L 程g 已知，f 唯一确定。设n n n

f f =?∑（n ?是系数，n f 是基函数）。其

中{}n f 不一定是完备的，n 越大，不一定越好，只有{}n f 完备，n 越大越接近真值。方程离散为：

n n

n

Lf

g ?=∑

此过程的主要目的在于将算子方程转为代数方程，具体步骤如下： (1) 在算子L 的定义域内适当的选择一组基函数，123,,...n f f f f ，他们应该是线性无关的。

(2) 将未知函数()f x 表示为该组基函数的线性组合，并取有限项近似，即:()1

1

()N

n n N n n n n f x f f x f ∞

===?≈=?∑∑ (2-1)

(3) 将公式（2-1）利用算子的线性，可将算子方程化为代数方程，即:

1

()N

n n

n L f

g =?=∑ (2-2)

于是，待求()f x 的问题化解为求解n f 的系数n ?的问题。

(二)取样检验过程

为了使()f x 的近似函数()N f x 之间的误差极小，必须进行抽样检验，

在抽样点上使加权平均误差为零，从而确定未知系数n ?，这一过程的基本步骤如下：

(1) 在算子L 的值域内取适当的选择一组加权函数123,,...n w w w w ，他们也是线性无关的。

(2) 将m w 与式（2-2）取内积进行抽样检验，因为要确定N 个未知数，需要进行N 次抽样检验，则

1

(),,N

n

n m m n L f w g w =?=∑

(1,2,3...)m N = （2-3）

(3) 利用算子的线性和内积的性质，将（2-3）式转化为矩阵方程，即:

1

(),,N

n

n m m n L f w g w =?=∑

(1,2,3...)m N =（选配过程） （2-4）

将它写成矩阵形式：

[][][]mn n n l g ?= (2-5) 于是，求解代数方程问题转化为求解矩阵方程的问题。

(三)矩阵求逆

一旦求得了矩阵方程，通过常规的矩阵求逆或求解线性方程组，就可以得到矩阵方程的解。如果[]mn l 是非奇异的，则它的逆矩阵 []1

mn l -是存在的，则矩阵方程的解为：

[][][]1

n mn m l g -?= （2-6） 将求得的展开系数n ?带入到（2-1）式中，就可以得到原来算子方程的近似解为：

()1()N

N n n n f x f x ==?∑ (2-7)

以上所述是矩量法求解算子方程的基本过程，在矩量法的所有应用中，通常都要遵循这个统一的过程。由于()f x 和()N f x 之间存在近似，故算子方程的左端近似值与右端精确值之间存在如下关系：

()1(())()N

n n n x L f x g x ε==?-∑ （2-8）

三基函数与权函数的选择

矩量法的求解过程是简单的，求解步骤是统一的，应用起来比较方便。求解精度取决于很多因素，例如离散化的成都、基函数和权函数的选取和矩阵的求解过程等，其中基函数与权函数的选择尤其重要，可以直接影响到求解的精度。

基函数可以分为全域基和分域基，权函数可以分为权函数、分域权和点匹配，它们之间的不同组合便形成不同的方法。下面简单的对此进行介绍。

1、 全域基函数：全域基函数是指算子L 定义域内的全域上存在

一组基函数。他们应该满足边界条件。在矩量法求解的离散化过程中，选择全域基函数作为展开函数，实质上将未知函数表示为全域存在的若干个离散化基函数的线性组合。 2、 分域基函数：分域基函数不是在算子L 定义域的全域上存在

的，而仅仅存在于算子定义域的各个分域上的函数。选择分域基函数作为未知函数的展开函数，在矩量法求解的离散化过程中是一种区域离散，即未知函数表示为各个分域上存在的函数线性组合。

3、 全域权：在算子L 的值域内选择一组权函数123,,...m w w w w ，它

是一组在L 值域内的全域上存在的权函数。在此方法中，如果选择n n w f =，即权函数等于基函数，则称此方法为伽略金法。

4、 点匹配:全域基、全域权的伽略金发是一种常用的求解方法。

但是，如果算子本身是复杂的积分算子，或者选择了比较复杂的基函数，由于内积运算本身又是积分运算，从而使矩阵元素mn l 和m g 的形成十分困难。这时为了简化，我们可以利用

δ函数的选择性将δ函数作为权函数将使内积计算得以简

化。

四算例

问题：计算带电导体板的电容。正方形板，边长为2a ，位于0z =。空间中任意一点电位()''

''

0,(,,)4a

a

a a x y x y z dx dy

R

σφπε--=??

，其中

R =。在板上V φ=为常数，有

()

'''

'

0,4a

a

a

a

x y V dx dy

R

σπε--=??，其中电荷分布()'',x y σ

为待求，

R =，电容为：

()''''1

,a

a

a

a

q C dx dy x y V V

σ--=

=?

?。计算当4a =时方板的归一化电容

2C a

。 1． 选取分域基：方板各边N 等分，22a

b N

=

，其中2M N =个子域。定义：10

n f ?=??，其中在n S ?上为1，否则为0 。将()2

''

1,n N

n n n x y f σ==≈?∑带入

()''

''

0,4a a

a

a

x y V dx dy

R

σπε--=??中，有

2

1

n N mn

n n V l

===

?∑ ，

其中''

a

a

mn a a

l dx dy --=??。mn l 是n S ?上单位振幅的均

匀电荷密度在m S ?的中心处产生的电位。因此有

2

11,1n N

n n mn n n m n

C S l S V =-=≈??=?∑∑

物体的电容式其各小块电容的总和加上每一对小块间的互电容。

()()()

''''''

,,a a

a

a

x y L f L x y dx dy V σσ--===??

2． 取()()m m m w x x y y δδ=--

[][][][][]

[][][][]2

2

2

2

2

11111

1111mn n n mn N N N N N n n n mn M M M M l V l V C S S l V

ααα-????-???=?=??=?=???

3． mn l 的计算：

nn l ：由n S ?本身面上的单位电荷密度在其中心处产生的电位

（自作用

单元）

(0

2ln 1b

b

nn b

b

b

l dx πε--==

??

而互作用单元：

''

n

mn S l ?=

≈

??

经过计算模拟可以发现当N 取值越大时，结果将趋近于一个定值。当N 为20时，计算得到归一化电容为40.29PF ，非常接近理论计算值40PF ，误差为0.725%，在误差允许范围内。

程序：

#pragma once

class Martix

{

private:

c onst int M;

c onst int N;

b ool flag;

d oubl

e **array;

public:

M artix(double **_array, int _M, int _N) :M(_M), N(_N) {

array = _array;

flag = false;

}

M artix(int _M, int _N) :M(_M), N(_N)

{

flag = true;

array = new double *[M];

for (int i = 0; i < M; i++)array[i] = new double[N];

for (int i = 0; i < M; i++)

{

for (int j = 0; j < N; j++)

{

array[i][j] = 0;

}

}

}

M artix Gauss();

v oid show();

v oid swap(int i, int j);

i nt Max(int i);

M artix operator+(Martix &T);

M artix operator-(Martix &T);

M artix operator*(Martix &T);

M artix operator*(double t);

f riend Martix operator*(double t, const Martix &T)

{

Martix sum(T.M, T.N);

for (int i = 0; i < T.M; i++)

{

for (int j = 0; j < T.N; j++)

{

sum.array[i][j] = T.array[i][j] * t;

}

}

return sum;

}

M artix Copy();

d oubl

e **Getarray()

{

return this->array;

}

i nt GetM()

{

return this->M;

}

i nt GetN()

{

return this->N;

}

b ool GetFlag()

{

return this->flag;

}

~Martix();

};

#include "Martix.h"

#include

#include

using namespace std;

void Martix::show()

{

f or (int i = 0; i < M; i++)

{

cout << fixed << right;

for (int j = 0; j < N; j++)

{

cout << setw(10) << setprecision(6) << array[i][j] << " ";

}

cout << endl;

}

}

Martix Martix::Gauss()

{

M artix sum(M, N);

f or (int i = 0; i < N; i++)sum.array[i][i] = 1;

f or (int i = 0; i < N; i++)

{

if (abs(array[i][i]) <= 0.00001&&i

{

if (abs(array[Max(i)][i]) <= 0.00001)continue;

int count = Max(i);

swap(i, count);

sum.swap(i, count);

}

else if (abs(array[i][i]) < 0.00001&&i == N - 1)

{

cout << "没有逆矩阵" << endl;

exit(0);

}

/*cout << "变换前" << endl;

this->show();*/

double temp = array[i][i];

for (int j = 0; j < N; j++)

{

array[i][j] = array[i][j] / temp;

sum.array[i][j] = sum.array[i][j] / temp;

}

for (int k = i + 1; k < N; k++)

{

double temp = array[k][i];

if (abs(temp) >= 0.00001)

{

for (int j = 0; j < N; j++)

{

array[k][j] = array[k][j] - temp*array[i][j];

sum.array[k][j] = sum.array[k][j] - temp*sum.array[i][j];

}

}

}

}

d oubl

e t = abs(array[0][0]);

f or (int i = 1; i < N; i++)t = t*abs(this->array[i][i]);

i f (abs(t) <= 0.0001)

{

cout << "没有逆矩阵" << endl;

exit(0);

}

f or (int i = N - 2; i >= 0; i--)

{

for (int j = 0; j < N; j++)

{

for (int k = i + 1; k < N; k++)

{

sum.array[i][j] = sum.array[i][j] - array[i][k] * sum.array[k][j];

}

}

}

r eturn sum;

}

void Martix::swap(int i, int j)

{

d oubl

e temp = 0;

f or (int k = 0; k < N; k++)

{

temp = array[i][k];

array[i][k] = array[j][k];

array[j][k] = temp;

}

}

Martix Martix::Copy()

{

M artix sum(this->M, this->N);

f or (int i = 0; i < M; i++)

{

for (int j = 0; j < N; j++)

{

sum.array[i][j] = this->array[i][j];

}

}

r eturn sum;

}

Martix Martix::operator+(Martix &T)

{

M artix sum(T.M, T.N);

i f (M != T.M || N != T.N)

{

cout << "错误" << endl;

}

f or (int i = 0; i < M; i++)

{

for (int j = 0; j < N; j++)

{

sum.array[i][j] = this->array[i][j] + T.array[i][j];

}

}

r eturn sum;

}

Martix Martix::operator-(Martix &T)

{

M artix sum(T.M, T.N);

i f (M != T.M || N != T.N)

{

cout << "错误" << endl;

}

f or (int i = 0; i < M; i++)

{

for (int j = 0; j < N; j++)

{

sum.array[i][j] = this->array[i][j] - T.array[i][j];

}

}

r eturn sum;

}

Martix Martix::operator*(Martix &T)

{

M artix sum(T.M, T.N);

i f (N != T.M)

{

cout << "错误" << endl;

}

f or (int i = 0; i < M; i++)

{

for (int j = 0; j < T.N; j++)

{

for (int k = 0; k < N; k++)

{

sum.array[i][j] = sum.array[i][j] + this->array[i][k] * T.array[k][j];

}

}

}

r eturn sum;

}

Martix Martix::operator*(double t)

{

M artix sum(this->M, this->N);

f or (int i = 0; i < this->M; i++)

{

for (int j = 0; j < this->N; j++)

{

sum.array[i][j] = this->array[i][j] * t;

}

}

r eturn sum;

}

int Martix::Max(int i)

{

i nt count = i;

d oubl

e temp = abs(array[i][i]);

f or (int k = i + 1; k < N; k++)

{

if (temp < abs(array[k][i]))temp = array[k][i];

}

w hile (temp != abs(array[count][i]))

{

count++;

}

r eturn count;

}

Martix::~Martix()

{

}

#include

#include

#include

#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.GetM(); i++)delete[](T.Getarray())[i]; }

double lmn(int m, int n)

{

double l=0;

int i1=0, j1=m%N, i2=0, j2=n%N;

if (m == n)

{

l = 2 * b*0.8814 / pi ;;

}

else

{

if (m%N != 0)

{

i1 = m / N + 1;

}

else

{

i1 = m / N;

}

if (n%N != 0)

{

i2 = n / N + 1;

}

else

{

i2 = n / N;

}

double xmn = (i1-i2)*2*b;

double ymn = (j1-j2)*2*b;

l = ds / (4*pi*sqrt(xmn*xmn+ymn*ymn));

}

return l;

}

void main()

{

double **L = new double *[M];

double c = 0;

for (int i = 0; i < M; i++)L[i] = new double[M];

for (int i = 0; i < M; i++)

{

for (int j = 0; j < M; j++)

{

L[i][j] = lmn(i+1, j+1);

}

}

Martix l(L,M, M);

/*cout << "矩阵L:" << endl;

l.show();*/

Martix ln = l.Gauss();

/*cout << "逆矩阵：" << endl;

ln.show();*/

for (int i = 0; i < ln.GetM(); i++)

{

for (int j = 0; j < ln.GetN(); j++)

{

c = c + (ln.Getarray())[i][j];

}

}

c = c*ds*eps;

cout << "电容为：" << c/(2*a) << endl;

clean(l);

clean(ln);

}