CN103400004A - 基于多区域模型矩量法的介质粗糙面电磁散射仿真方法 - Google Patents

Abstract

本发明公开了一种基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，包括：输入介质粗糙面功率谱密度函数和粗糙面参数，用蒙特卡洛法得到粗糙面上N个离散面元位置坐标；将粗糙面划分为主区域与非主区域；分别获取这两个区域面元上的等效电流；得到该粗糙面的电磁散射系数，实现介质粗糙面电磁散射的仿真。本发明对粗糙面进行区域划分，仅在对散射结果影响较大的主区域，运用矩量法建立边界积分方程并获取主区域等效电流，在非主区域，近似获取等效电流，在保证一定精度的前提下减少未知量个数，降低内存，提高计算效率，解决了实际大区域介质粗糙面电磁散射仿真的技术问题。用于国防和民用领域的粗糙面雷达散射特性研究。

Description

基于多区域模型矩量法的介质粗糙面电磁散射仿真方法

技术领域

本发明属于雷达电磁仿真技术领域，主要涉及电磁散射数值仿真，具体是一种基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，用于获取介质粗糙面的电磁散射系数。

背景技术

随着雷达技术的快速发展，粗糙面的电磁散射研究在理论分析与实际应用中具有重要的意义。当机载或者星载雷达对粗糙面等背景中的目标进行电磁探测时，雷达的回波信号中包含了粗糙面背景的电磁散射信号，通过分析这些回波信号，进而可以得出粗糙面以及探测目标的某些特性。因而，粗糙面雷达散射特性的研究在国防领域和民用领域具有显著的学术价值和广泛的应用前景。

在过去的几十年中，许多电磁仿真技术被学者提出用以处理地海面的电磁散射问题，大致分为近似方法和数值方法。近似方法的优点是消耗内存低、分析速度快，然而近似方法都是基于特定的物理近似，其精度往往较低并且不具有通用性。相比于近似方法，数值方法能保持较高的仿真精度。矩量法作为一种典型的数值方法，由于其计算结果的高精度得到广泛的应用。但该方法也存在一定缺陷，运用传统矩量法对粗糙面的电磁散射进行仿真时，首先要对粗糙面的表面轮廓进行离散，然后在整个介质粗糙面上建立边界积分方程，并对任意两个离散面元间的相互作用进行计算，最后求解一个Ax＝b的矩阵方程，其中A称为阻抗矩阵，如果粗糙表面被离散成N份，则A为一个N×N大小的矩阵，A中的每个元素代表两个离散面元之间的相互作用；b称为电压向量，是一个长度为N的向量，其中元素与入射波信息有关；x为所要求的未知向量，即粗糙表面上每个面元上的等效电流。可以看出一旦计算对象的表面积增大，粗糙面被剖分的份数N随之增大，矩量法计算时将产生大量的未知量，阻抗矩阵的规模将会以N²的数量级增长，内存消耗也会以N²的数量级增长，同时求解矩阵方程将会耗费大量机时。而对于实际的粗糙面，通常区域较大，采用单纯的矩量法对其散射特性进行仿真时，由于其内存消耗大，分析速度慢，普通个人计算机的配置已不能满足仿真需求，故该方法很难应用到实际大区域粗糙面的电磁散射仿真中。

发明内容

本发明的目的在于克服已有技术的不足，针对大区域粗糙面的电磁散射仿真问题，提供了一种基于多区域模型的介质粗糙面矩量法电磁仿真方法，在保证一定精度的前提下减少了未知量个数，降低了仿真的内存需求，提高了仿真效率。

实现本发明目的的技术方案，包括如下步骤：

(1)输入介质粗糙面的功率谱密度函数和粗糙面参数，通过蒙特卡洛方法得到粗糙面上N个离散面元的位置坐标r_-N/2+1＝(x_-N/2+1，z_-N/2+1)，…，r_-1＝(x_-1，z_-1)，r₀＝(x₀，z₀)，r₁＝(x₁，z₁)，…，r_N/2＝(x_N/2，Z_N/2)，即对介质粗糙面进行建模，也就是生成介质粗糙面。

(2)使用锥形波入射到该仿真的粗糙面上，将入射波强度大于e^-1/4的面元划分进主区域，入射波强度小于e^-1/4的面元划分进非主区域。

(3)在主区域中，根据电磁场中的边界条件，获得介质粗糙面主区域的矩量法边界积分方程，并将该方程转化为矩阵方程，求解该矩阵方程，得到介质粗糙面每个主区域面元上的等效电流。

(4)在非主区域的每个面元上，仅考虑入射磁场和主区域上等效电流在其表面的感应，得到每个非主区域面元上的等效电流。

(5)利用已得到的介质粗糙面每个面元上的等效电流，获取整个介质粗糙面的电磁散射系数，完成该介质粗糙面的电磁散射仿真。

本发明的实现还在于＝对介质粗糙面进行区域划分，包括如下步骤：

(2a)采用锥形波照射粗糙面，锥形波的入射电场E^inc(r)和入射磁场H^inc(r)分别表示如下：

$E^{\mathrm{inc}}\left(r\right)=\hat{y}\exp [-j{k}_a(x\sin {\mathrm{\θ}}_i-z\cos {\mathrm{\θ}}_i)(1+w\left(r\right))]\×\exp (-\frac{{(x+z\tan {\mathrm{\θ}}_i)}^2}{g^2})$

$H^{\mathrm{inc}}\left(r\right)=\hat{k}\×E^{\mathrm{inc}}\left(r\right)/Z_a$

其中 $w\left(r\right)=\exp [2{(x-\mathrm{yctg}{\mathrm{\θ}}_i)}^2/g^2-1]/{\left(k_{a}g\sin {\mathrm{\θ}}_i\right)}^2,$

是空间波数的单位矢量，Z_a为上半空间波阻抗，μ_a为上半空间相对磁介电常数，ε_a为上半空间相对介电常数，k_a为上半空间中的波数，r＝(x，z)是面元的位置坐标，g是锥形波因子，θ_i为入射角。

(2b)依次读入步骤(1)中生成介质粗糙面上的每个面元的位置坐标r_-N/2+1＝(x_-N/2+1，z_-N/2+1)，…，r_-1＝(x_-1，z_-1)，r₀＝(x₀，z₀)，r₁＝(x₁，z₁)，…，r_N/2＝(x_N/2，Z_N/2)，将面元位置坐标带入步骤(2a)中，依次得到每个面元上的入射电场E^inc(r)和入射磁场H^inc(r)，并得到每个面元上电场的幅值|E^inc(r)|。

(2c)设置门限值为e^-1/4，将入射波电场幅值|E^inc(r)|大于门限值的面元划分进主区域，入射波电场幅值|E^inc(r)|小于门限值的面元划分进非主区域。

(2d)将位置坐标r的下标小于0的主区域面元划分为1l区域，由这些面元组成的区域用符号S_1l表示，将位置坐标r的下标大于等于0的主区域面元划分进1r区域，这些面元组成的区域用符号S_1r表示。

(2e)与步骤(2d)类似，将位置坐标r的下标小于0的非主区域面元划分为2l区域，由这些面元组成的区域用符号S_2l表示，将位置坐标r的下标大于等于0的非主区域面元划分进2r区域，这些面元组成的区域用符号S_2r表示。

本步骤进行了介质粗糙面的区域划分，即引入多区域模型，以便下一步对不同区域进行不同处理，该多区域模型将入射波较强的面元划入主区域，由于该区域对电磁散射仿真结果贡献较大，在后续步骤中仅在该区域运用矩量法建立边界积分方程并求解主区域等效电流；入射波较弱的面元划入非主区域，该区域对电磁仿真结果影响相对较小，在后续步骤中该区域的等效电流将近似求解。

本发明的实现还在于：获得每个主区域面元上的等效电流，包括如下步骤：

(3a)在主区域建立边界积分方程如下：

J可以代指主区域上每个面元上的等效电流，在主区域的每个面元上，该方程均成立，上式中，下标a和b分别表示粗糙面上方和下方，

是对应面元单位法向量，Z_b为下半空间波阻抗，μ_b为下半空间相对磁介电常数，ε_b为下半空间相对介电常数，算子

和

分别定义为：

${\stackrel{\‾}{L}}_a\left(X\right)=-j{k}_a\underset{S}{\∫}[X+\frac{1}{k_a^2}\▿(\▿\·X)]G_a(r,r^{\′})d{s}^{\′}$

${\stackrel{\‾}{K}}_a\left(X\right)=-\underset{S}{\∫}X\×\▿G_a(r,r^{\′})d{s}^{\′}$

${\stackrel{\‾}{L}}_b\left(X\right)=-j{k}_b\underset{S}{\∫}[X+\frac{1}{k_b^2}\▿(\▿\·X)]G_b(r,r^{\′})d{s}^{\′}$

${\stackrel{\‾}{K}}_b\left(X\right)=-\underset{S}{\∫}X\×\▿G_b(r,r^{\′})d{s}^{\′}$

表示对整个介质粗糙面进行积分，G_a(r，r′)为粗糙面上半空间中的格林函数，G_b(r，r′)为粗糙面下半空间中的格林函数。

(3b)采用狄拉克函数作为基函数，用伽辽金方法，将上式转化为矩阵方程：

$\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\·\left[\begin{array}{c}{J}_{1l}\\ J_{1r}\end{array}\right]=\begin{array}{}\left[\begin{array}{c}{b}_1\\ b_2\end{array}\right]\end{array}$

式中，J_1l表示一个向量，其中的每个元素对应S_1l中每个面元上的等效电流，J_1r表示一个向量，其中的每个元素对应S_1r中每个面元上的等效电流，A是子矩阵元素，b为入射波的列矩阵。

(3c)采用LU矩阵分解方法求解步骤(3b)中的矩阵方程，得到主区域的等效电流分布J_1l，J_1r。

本步骤在步骤(2)划分区域的基础上，得到S_1l和S_1r这两个主区域面元上的等效电流J_1l，J_1r，该步骤得到的主区域等效电流可用于下一步求解非主区域的等效电流，在最后进行整个介质粗糙面散射系数仿真时，也需用到本步骤所得的主区域等效电流。

传统矩量法需在整个介质粗糙面上建立积分方程，未知量个数为N，在将整个介质粗糙面的积分方程转化为矩阵方程后，内存需求为N²的倍数；而本步骤中的矩量法是基于多区域模型的，即只在主区域建立积分方程，未知量个数为N_d，在将主区域的积分方程转化为矩阵方程后，内存需求为

的倍数。

本发明的实现还在于：获得每个非主区域面元上的等效电流，包括如下步骤：

(4a)根据步骤(2b)得到的区域S_2l中任意一个面元上的入射磁场

通过下式：

$J_{2l}^{\mathrm{inc}}\left(r\right)=2\hat{n}\×H_{2l}^{\mathrm{inc}}\left(r\right)$

求得区域S_2l中任意一个面元上由入射场感应产生的等效电流

其中，r为该面元的位置坐标，

为面元的单位法向矢量。

(4b)根据步骤(3c)得到S_1l中每一个面元上的电流分布J_1l，S_2l中任意一个面元上的感应得到的等效电流

(r)可通过下式：

$J_{2l}^{S_{1l}}\left(r\right)=-2\hat{n}\×\underset{S_{1l}}{\mathrm{\Σ}}{J}_{1l}(r\′)\×\▿G_a(r,r\′)\·\mathrm{\Δx}$

依次将S_1l中每个面元上的等效电流对S_2l中该面元的作用进行叠加，其中r′对应区域S_1l中个面元的位置坐标，G_a(r，r′)为粗糙面上半空间中的格林函数，

为该面元的单位法向矢量。

(4c)对S_2l每个面元执行步骤(4a)，(4b)，得到S_2l每个面元上由入射场感应产生的等效电流

(r)和J_1l在S_2l上感应产生的等效电流

(r)，将这两项相加，得到S_2l中每个面元上的等效电流J_2l(r)，即 $J_{2l}\left(r\right)=J_{2l}^{\mathrm{inc}}\left(r\right)+J_{2l}^{S_{1l}}\left(r\right);$

(4d)根据步骤(2b)可得到区域S_2r中任意一个面元上的入射磁场

(r)，通过下式：

$J_{2r}^{\mathrm{inc}}\left(r\right)=2\hat{n}\×H_{2r}^{\mathrm{inc}}\left(r\right)$

求得区域S_2r中任意一个面元上由入射场感应产生的等效电流

(r)，其中，r为该面元的位置坐标，

为面元的单位法向矢量。

(4e)根据步骤(3c)得到S_1r中每一个面元上的电流分布J_1r，S_2r中任意一个面元上的感应得到的等效电流

(r)可通过下式：

$J_{2r}^{S_{1r}}\left(r\right)=-2\hat{n}\underset{S_{1r}}{\mathrm{\Σ}}{J}_{1r}(r\′)\×\▿G_a(r,r\′)\·\mathrm{\Δx}$

依次将S_1r中每个面元上的等效电流对S_2r中该面元的作用进行叠加，其中r′对应区域S_1r中个面元的位置坐标。

(4f)对S_2r每个面元执行步骤(4d)，(4e)，得到S_2r每个面元上由入射场感应产生的等效电流(r)和J_1r在S_2r上感应产生的等效电流

(r)，将这两项相加，得到S_2r中每个面元上的等效电流J_2r(r)，即

本步骤在步骤(2)划分区域的基础上，得到S_2l和S_2r这两个主区域面元上的等效电流J_2l，J_2r，与传统矩量法不同，该区域没有通过求解矩阵方程求解，而是根据步骤(2b)得到的入射磁场和步骤(3c)得到的主区域等效电流J_1l，J_1r直接求解，减小了内存需求和计算量。

与现有技术相比，本发明的技术优势：

第一，由于使用了多区域模型，仅在该区域运用矩量法建立边界积分方程，并将该积分方程转化为矩阵方程，假设整个介质粗糙面共有N个面元，主区域含有N_d个面元，则传统矩量法的未知量个数为N，内存需求为N²的倍数，而本发明只含有N_d个未知量，内存需求为

的倍数，故本发明可降低计算机内存需求。

第二，由于采用了多区域模型，对于求解的矩阵方程Ax＝b，传统矩量法中A为一个N×N大小的矩阵，而本发明中A的大小仅为N_d×N_d，矩阵方程规模减小，求解矩阵方程所需的时间随之减小，计算效率随之提高。

附图说明

图1是本发明的流程图；

图2是介质粗糙面轮廓图；

图3是粗糙面主区域和非主区域划分示意图；

图4是θ_i＝30°时本发明和传统矩量法仿真得到的粗糙面等效电流对比曲线图；

图5是θ_i＝30°时本发明和传统矩量法仿真得到介质粗糙面散射系数对比曲线图；

图6是θ_i＝70°时本发明和传统矩量法仿真得到的粗糙面等效电流对比曲线图；

图7是θ_i＝70°时本发明和传统矩量法仿真得到介质粗糙面散射系数对比曲线图。

具体实施方式

下面结合附图和实施例对本发明再做说明。

实施例1

本发明是一种基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，当机载或者星载雷达对粗糙面等背景中的目标进行电磁探测时，雷达的回波信号中包含了粗糙面背景的电磁散射信号，对这些回波信号分析，可得出粗糙面以及探测目标的某些特性。因此，粗糙面的电磁散射特性研究具有重要意义，本发明即是通过实验数据得到粗糙面的参数，输入这些参数建立仿真的介质粗糙面，再运用基于多区域模型矩量法对介质粗糙面电磁散射进行仿真，本发明得到的仿真数据可用于后续的粗糙面散射特性分析研究。本发明所需的实验设备是一台具有Fortran语言编译环境的计算机。

本发明的基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，参照图1，本发明的具体实现步骤如下：

(1)输入介质粗糙面的功率谱密度函数和粗糙面参数，通过蒙特卡洛方法得到粗糙面上N个离散面元的位置坐标r_-N/2+1＝(x_-N/2+1，z_-N/2+1)，…，r_-1＝(x_-1，z_-1)，r₀＝(x₀，z₀)，r₁＝(x₁，z₁)，…，r_N/2＝(x_N/2，Z_N/2)，即对介质粗糙面进行建模，也就是生成介质粗糙面。

(1a)通过实验获得所要仿真介质粗糙面的均方根高度、相关长度的粗糙面参数和功率谱密度函数，将其作为输入参数，在输入参数的作用下，选取一系列振幅独立的高斯随机变量的高斯谐波。

(1b)输入入射波的频率f，以Δx＝3.0×10⁸/(10f)为采样间隔，取N个离散采样点x_-N/2+1，…，x_-1，x₀，x₁，…，x_N/2，令x₀＝0，则x_-1＝-Δx，x₁＝Δx，…x_-N/2+1＝(-N/2+1)·Δx，x_N/2＝N/2·Δx。

(1c)将这N个离散采样点作为粗糙面的横坐标，在离散点处通过蒙特卡洛方法将大量的振幅独立的高斯谐波，即步骤(1a)所选取的一系列振幅独立的高斯随机变量的高斯谐波进行叠加，然后作逆傅里叶变换得到每个采样点对应的纵坐标z_-N/2+1，…，z_-1，z₀，z₁，…，Z_N/2，生成所要仿真的介质粗糙面。该仿真的介质粗糙面是具有位置坐标的轮廓面，参见图2。

以上步骤实现了介质粗糙面的轮廓仿真，得到了每个面元的位置坐标，这些位置坐标为步骤(2)的区域划分提供依据。

(2)使用锥形波入射到该仿真的粗糙面上，将入射波强度大于e^-1/4的面元划分进主区域，入射波强度小于e^-1/4的面元划分进非主区域。

(2a)采用锥形波照射粗糙面，锥形波的入射电场E^inc(r)和入射磁场H^inc(r)分别表示如下：

$E^{\mathrm{inc}}\left(r\right)=\hat{y}\exp [-{\mathrm{jk}}_a(x\sin {\mathrm{\θ}}_i-z\cos {\mathrm{\θ}}_i)(1+w\left(r\right))]\×\exp (-\frac{(x+z\tan {\mathrm{\θ}}_i{)}^2}{g^2})$

H^inc(r)＝k×E^inc(r)/Z_a

其中 $w\left(r\right)=\exp [2(x-\mathrm{yctg}{\mathrm{\θ}}_i{)}^2/g^2-1]/(k_{a}g\sin {\mathrm{\θ}}_i{)}^2,$ θ_i为入射角，

是空间波数的单位矢量，

Z_a为上半空间波阻抗，

μ_a为上半空间相对磁介电常数，ε_a为上半空间相对介电常数，k_a为上半空间中的波数，

r＝(x，z)是面元的位置坐标，g是锥形入射波因子。

(2b)依次读入步骤(1)中生成介质粗糙面上的每个面元的位置坐标r_-N/2+1＝(x_-N/2+1，z_-N/2+1)，…，r_-1＝(x_-1，z_-1)，r₀＝(x₀，z₀)，r₁＝(x₁，z₁)，…，r_N/2＝(x_N/2，Z_N/2)，将面元位置坐标带入步骤(2a)中，依次可得到每个面元上的入射电场E^inc(r)和入射磁场H^inc(r)，并且求得每个面元上电场的幅值|E^inc(r)|。

(2c)设置门限值为e^-1/4，将入射波电场幅值|E^inc(r)|大于门限值的面元划分为主区域，入射波电场幅值|E^inc(r)|小于门限值的面元划分为非主区域。

(2d)将位置坐标r的下标小于0的主区域面元划分为1l区域，由这些面元组成的区域用符号S_1l表示，将位置坐标r的下标大于等于0的主区域面元划分进1r区域，这些面元组成的区域用符号S_1r表示。

(2e)与步骤(2d)类似，将位置坐标r的下标小于0的非主区域面元划分为2l区域，由这些面元组成的区域用符号S_2l表示，将位置坐标r的下标大于等于0的非主区域面元划分进2r区域，这些面元组成的区域用符号S_2r表示。

本发明在介质粗糙面轮廓模拟的基础上进行了介质粗糙面的区域划分，即引入多区域模型，以便下一步对不同区域进行不同处理，该多区域模型将入射波较强的面元划入主区域，由于该区域对电磁散射仿真结果贡献较大，在后续步骤中仅在该区域运用矩量法建立边界积分方程并求解主区域等效电流；入射波较弱的面元划入非主区域，该区域对电磁仿真结果影响相对较小，在后续步骤中该区域的等效电流将近似求解。

值得注意的是，根据(2a)中入射电场的表达形式，在整个介质粗糙面上，位于介质粗糙面最中间的面元入射电场幅值最大，面元越靠近边缘，电场幅值越小，所以区域S_1l，S_1r位于粗糙面靠中间的部分，S_2l与S_2r分别位于粗糙面的边缘部分，如图3所示。

(3)在主区域中，根据电磁场中的边界条件，获得介质粗糙面主区域的矩量法边界积分方程，并将该方程转化为矩阵方程，求解该矩阵方程，得到介质粗糙面每个主区域面元上的等效电流。

(3a)在主区域建立边界积分方程如下：

J可以代指主区域上每个面元上的等效电流，在主区域的每个面元上，该方程均成立，上式中，下标a和b分别表示粗糙面上方和下方，

是对应面元单位法向量，Z_b为下半空间波阻抗，μ_b为下半空间相对磁介电常数，ε_b为下半空间相对介电常数算子和

分别定义为：

${\stackrel{\‾}{L}}_a\left(X\right)\text{=-}{\mathrm{jk}}_a\underset{S}{\∫}[X+\frac{1}{k_a^2}\▿(\▿\·X)]G_a(r,r\′)\mathrm{ds}\′$

${\stackrel{\‾}{K}}_a\left(X\right)=-\underset{S}{\∫}X\×\▿G_a(r,r\′)\mathrm{ds}\′$

${\stackrel{\‾}{L}}_b\left(X\right)=-{\mathrm{jk}}_b\underset{S}{\∫}[X+\frac{1}{k_b^2}\▿(\▿\·X)]G_b(r,r\′)\mathrm{ds}\′$

${\stackrel{\‾}{K}}_b\left(X\right)=-\underset{S}{\text{\∫}}X\×\▿G_b(r,r\′)\mathrm{ds}\′$

表示对整个介质粗糙面进行积分，G_a(r，r′)为粗糙面上半空间中的格林函数，其中j为虚数单位，

为0阶第二类汉克尔函数，G_b(r，r′)其中为粗糙面下半空间中的格林函数，r，r′分别是作为场点和源点的面元位置坐标。

(3b)采用狄拉克函数作为基函数，使用伽辽金方法，将上式转化为矩阵方程：

$\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\·\begin{array}{c}\left[\begin{array}{c}{J}_{1l}\\ J_{1r}\end{array}\right]\end{array}=\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right]$

式中，J_1l表示一个向量，其中的每个元素对应S_1l区域中每个面元上的等效电流，J_1r表示一个向量，其中的每个元素对应S_1r区域中每个面元上的等效电流，A是子矩阵元素，B为入射波的列矩阵。

(3c)采用LU矩阵分解方法求解步骤(3b)中的矩阵方程，得到主区域的等效电流分布J_1l，J_1r。

步骤(3b)中积分方程转化为矩阵方程的步骤是传统矩量法中的经典转化方法，其中，基函数的定义和伽辽金方法的具体内容可参见R.F.Harrington所著，王尔杰译的《计算电磁场的矩量法》，步骤(3c)中的LU矩阵分解方法为一种经典求解矩阵方程的方法，详细内容可见程云鹏所著的《矩阵论》。

由于本发明使用了多区域模型，仅在该区域运用矩量法建立边界积分方程，并将该积分方程转化为矩阵方程，假设整个介质粗糙面共有N个面元，主区域含有N_d个面元，则传统矩量法的未知量个数为N，内存需求为N²的倍数，而本发明只含有N_d个未知量，内存需求为

的倍数，故本发明降低了计算机内存需求，进而加快了仿真速度。

(4)在非主区域的每个面元上，本发明仅考虑入射磁场和主区域上等效电流在其表面的感应，得到每个非主区域面元上的等效电流J。

(4a)根据步骤(2b)得到的区域S_2l中任意一个面元上的入射磁场

通过下式：

$J_{2l}^{\mathrm{inc}}\left(r\right)=2\hat{n}\×H_{2l}^{\mathrm{inc}}\left(r\right)$

求得区域S_2l中任意一个面元上由入射场感应产生的等效电流

其中，r为该面元的位置坐标，为面元的单位法向矢量。

(4b)根据步骤(3c)得到S_1l中每一个面元上的电流分布J_1l，S_2l中任意一个面元上的感应得到的等效电流可通过下式：

$J_{2l}^{S_{1l}}\left(r\right)=-2\hat{n}\×\underset{S_{1l}}{\mathrm{\Σ}}{J}_{1l}(r\′)\×\▿G_a(r,r\′)\·\mathrm{\Δx}$

依次将S_1l中每个面元上的等效电流对S_2l中该面元的作用进行叠加，其中r′对应区域S_1l中个面元的位置坐标，G_a(r，r′)为粗糙面上半空间中的格林函数，

为该面元的单位法向矢量。

(4c)对S_2l每个面元执行步骤(4a)，(4b)，得到S_2l每个面元上由入射场感应产生的等效电流和j_1l在S_2l上感应产生的等效电流

(r)，将这两项相加，得到S_2l中每个面元上的等效电流J_2l(r)，即

(4d)根据步骤(2b)可得到区域S_2r中任意一个面元上的入射磁场

，通过下式：

$J_{2r}^{\mathrm{inc}}\left(r\right)=2\hat{n}\×H_{2r}^{\mathrm{inc}}\left(r\right)$

求得区域S_2r中任意一个面元上由入射场感应产生的等效电流

，其中，r为该面元的位置坐标，

为面元的单位法向矢量。

(4e)根据步骤(3c)得到S_1r中每一个面元上的电流分布x_1r，S_2r中任意一个面元上的感应得到的等效电流

可通过下式：

$J_{2r}^{S_{1l}}\left(r\right)=-2\hat{n}\×\underset{S_{1r}}{\mathrm{\Σ}}{J}_{1r}(r\′)\×\▿G_a(r,r\′)\·\mathrm{\Δx}$

依次将S_1r中每个面元上的等效电流对S_2r中该面元的作用进行叠加，其中r′对应区域S_1r中个面元的位置坐标。

(4f)对S_2r每个面元执行步骤(4d)，(4e)，得到S_2r每个面元上由入射场感应产生的等效电流

和J_1r在S_2r上感应产生的等效电流

将这两项相加，得到S_2r中每个面元上的等效电流J_2r(r)，即

(5)利用己得到的介质粗糙面每个面元上的等效电流，获取整个介质粗糙面的电磁散射系数，完成该介质粗糙面的电磁散射仿真。

由于采用了多区域模型，矩阵方程规模减小，求解矩阵方程所需的时间随之减小，计算效率随之提高。

实施例2

基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，同实施例1，其中，步骤(5)的具体实现方法如下：

(5a)由表面等效电流分布得到远区散射场E^sca，具体步骤如下式：

$E^{\mathrm{sca}}=Z_a\underset{p=1}{\overset{2}{\mathrm{\Σ}}}\underset{S_{\mathrm{pl}}}{\∫}{J}_{\mathrm{pl}}{G}_a(r,r^{\′}){\mathrm{ds}}^{\′}+Z_a\underset{p=1}{\overset{2}{\mathrm{\Σ}}}\underset{S_{\mathrm{pr}}}{\∫}{J}_{\mathrm{pr}}{G}_a(r,r\′)\mathrm{ds}\′$

p分别取1，2，使区域S_1lS_1rS_2l S_2r中每个面元等效电流辐射产生的远区场叠加。

(5b)根据双站雷达散射系数方程，得到介质粗糙面的双站雷达散射系数σ(θ)：

$\mathrm{\σ}\left(\mathrm{\θ}\right)=10\lg \left(2\mathrm{\πr}{\left|\frac{E^{\mathrm{sca}}}{E^{\mathrm{inc}}}\right|}^2\right)$

其中，θ是雷达的散射角度，r是仿真位置离原点的距离，E^sca是散射电场，E^inc表示入射电场。

实施例3

基于多区域模型矩量法的介质粗糙面电磁散射仿真方法同实施例1-2。

仿真条件：

介质粗糙面的粗糙面参数选取如下：粗糙面均方根高度δ＝0.02m，相关长度为l＝0.1m，粗糙面的谱密度函数选取高斯谱函数，

仿真中使用的雷达入射波频率为f＝3GHz，介质粗糙面的面元离散间隔为Δx＝0.01m，介质粗糙面上离散采样点数N＝1024，粗糙面的总长度为L＝NΔx＝10.24m，S_1l区域的长度为1.28m，S_1r区域长度为1.28m，S_2l区域长度为3.84m，S_2r区域的长度为3.84m，锥形波因子g＝L/4，粗糙面上方介质为真空，即μ_a＝1，ε_a＝(1，0)，下方介质的相对磁介电常数μ_b＝1，相对介电常数为ε_b＝(9.8，1.2)。入射波的角度θ_i＝30°。

仿真实验在CPU为Intel(R)Core(TM)i3，主频2.93GHz，可用内存为3.45GB的Windows7系统上用Compaq Visual Fortran6软件编程完成。

仿真内容及分析：

在执行完步骤(4)后可得到θ_i＝30°时整个介质粗糙面上每个面元的等效电流，即粗糙面的等效电流分布曲线图，如图4所示，图4将本发明仿真得到的等效电流分布曲线与传统矩量法得到的等效电流分布曲线进行了对比，图4中可看出在主区域上，本发明与传统矩量法得到的结果非常吻合，因为在主区域上，本发明仍使用矩量法建立积分方程并求解出等效电流分布，所以仿真结果具有精确性，在非主区域部分，本发明仿真得到的等效电流与传统矩量法相比存在一定的误差，这是因为本方法在非主区域的等效电流是近似得到的。本发明将有限的计算资源用于仿真对散射系数贡献较大的主区域等效电流，非主区域等效电流由于近似获得，计算资源消耗非常小。而传统矩量法在整个粗糙面上都采用精确解法，虽然结果准确，但计算资源消耗较大，仿真效率低下，尤其在粗糙面区域增大时，粗糙面上离散面元个数增多，未知量个数增多，内存需求与计算时间都随之增长。采用本发明可有效降低计算资源的消耗，并保证一定精度，在仿真实际大区域粗糙面时具有很大优势。

图4的结果证明了本发明在介质粗糙面表面电流分布的仿真中具有较好的的准确性。

实施例4

基于多区域模型矩量法的介质粗糙面电磁散射仿真方法同实施例1-2。本实施例的仿真条件同实施例3。

本实施例在实施例3得到介质粗糙面等效电流分布的基础上继续执行步骤(5)，得到介质粗糙面的电磁散射系数随散射角变化的曲线。

在图5中，将本发明仿真得到的介质粗糙面散射系数曲线和传统矩量法的仿真结果进行了对比，从图5可看出在-60°～+80°的散射角范围内本发明的仿真结果与传统矩量法结果吻合，而在其它散射角区域，由本发明仿真的雷达散射系数存在误差。该仿真实验证明了在-60°～+80°散射角范围内，本发明可有效准确地对雷达散射系数进行仿真。表1给出了本实施例的仿真条件下本发明与传统矩量法未知量个数和仿真时间对比。可看出本发明处理的未知量个数大大减少，仿真所需时间也大幅减少，变为传统矩量法的5.28％。由此可见，本发明提高了仿真效率。为实际的大区域粗糙面电磁散射提供了一种工程化的、快速有效的仿真方法。

表1θ_i＝30°时本发明与传统矩量法未知量个数与仿真时间对比

仿真方法	未知量个数	仿真时间(s)
			传统矩量法	1024	153.733
本发明	256	8.1184

实施例5

基于多区域模型矩量法的介质粗糙面电磁散射仿真方法同实施例1-2。

仿真条件：

介质粗糙面的粗糙面参数选取如下：粗糙面均方根高度δ＝0.02m，相关长度为l＝0.1m，粗糙面的谱密度函数选取高斯谱函数，

仿真中使用的雷达入射波频率为f＝3GHz，介质粗糙面的面元离散间隔为Δx＝0.0lm，介质粗糙面上离散采样点数N＝2048，粗糙面的总长度为L＝NΔx＝20.48m，S_1l区域的长度为2.56m，S_1r区域长度为2.56m，S_2l区域长度为7.68m，S_2r区域的长度为7.68m，锥形波因子g＝L/4，粗糙面上方介质为真空，即μ_a＝1，ε_a＝(1，0)，下方介质的相对磁相介电常数μ_b＝1，相对介电常数为ε_b＝(9.8，1.2)。入射波的角度θ_i＝70°。

仿真实验在CPU为Intel(R)Core(TM)i3，主频2.93GHz，可用内存为3.45GB的Windows7系统上用Compaq Nisual Fortran6软件编程完成。

仿真内容及分析：

本实施例实施方法如实施例1-2，在执行到完步骤(4)后可得到θ_i＝70°时整个介质粗糙面的等效电流分布曲线图，如图6所示，图6将θ_i＝70°时本发明仿真得到的等效电流分布曲线与传统矩量法得到的等效电流分布曲线进行了对比，图6中可看出在主区域上，本发明与传统矩量法得到的结果吻合很好，因为在主区域上，本发明仍使用矩量法建立积分方程并求解出等效电流分布，所以仿真结果具有精确性，在非主区域部分，本发明仿真得到的等效电流与传统矩量法相比存在一定的误差，这是因为本方法在非主区域的等效电流是近似得到的。每台计算机由于配置有限，本发明将有限的计算资源用于仿真对散射系数贡献较大的主区域等效电流，非主区域等效电流由于近似获得，计算资源消耗非常小。而传统矩量法在整个粗糙面上都采用精确解法，虽然结果准确，但计算资源消耗较大，仿真效率低下，尤其在粗糙面区域增大时，粗糙面上离散面元个数增多，未知量个数增多，内存需求与计算时间都随之增长。采用本发明可有效降低计算资源的消耗，并保证一定精度，在仿真实际大区域粗糙面时具有很大优势。

图6的结果证明了本发明在介质粗糙面表面电流分布的仿真中具有较好的的准确性。

实施例6

基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，同实施例1-2，本实施例的仿真条件与实施例5相同。

本实施例在实施例5得到θ_i＝70°时介质粗糙面等效电流分布的基础上继续执行步骤(5)，得到入射角θ_i＝70°时介质粗糙面的电磁散射系数随散射角变化的曲线。

从图7可看出在-60°～+80°的散射角范围内本发明的仿真结果与传统矩量法结果吻合，而在其它散射角区域，由本发明仿真的雷达散射系数存在误差。该仿真实验证明了在-60°～+80°散射角范围内，本发明可有效准确地对雷达散射系数进行仿真。同时，表2给出了本实施例的仿真条件下本发明与传统矩量法未知量个数和仿真时间对比。可看出本发明处理的未知量个数大大减少，仿真所需时间也大幅减少，变为传统矩量法的8.31％。由此可见，本发明提高了仿真效率。为实际的大区域粗糙面电磁散射提供了一种工程化的、快速有效的仿真方法。

表2θ_i＝70°时本发明与传统矩量法未知量个数与仿真时间对比

仿真方法	未知量个数	仿真时间(s)
			传统矩量法	2048	1405.634
本发明	512	116.883

以上描述仅是本发明在入射角为θ_i＝30°和θ_i＝70°时的具体实例，实际对于0°～90°范围内的任意入射角，本发明均能在-60°～+80°的散射角范围内获得准确的雷达散射系数°

以上描述仅是本发明的几个具体实例，对于本领域内的专业人员来说，在了解了本发明内容和原理后，可以进行形式上和细节上的各种修正和改变，但是这些基于本发明思想的修正和改变仍在本发明的权利要求保护范围之内。

简而言之，本发明的基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，与传统矩量法相比减小了未知量的个数，降低内存需求，提高计算效率，又在一定范围内上保证了其精确性。其实现步骤为(1)根据介质粗糙面的功率谱密度和粗糙面参数，通过蒙特卡洛方法产生仿真的介质粗糙面，并对仿真的粗糙面进行面元离散；(2)使用锥形波入射，根据每个面元上入射波的强度的大小将面元分别划分进主区域和非主区域；(3)在主区域中，根据电磁场中的边界条件，获得介质粗糙面离散主区域的矩量法边界积分方程，并将该方程转化为电磁耦合矩阵，求解该矩阵方程，得到介质粗糙面每个主区域面元上的等效电流；(4)在非主区域的每个面元上，仅考虑入射磁场和主区域上等效电流产生的散射电磁场在其表面的感应，得到每个非主区域面元上的等效电流；(5)利用以得到的介质粗糙面每个面元上的电磁流分布，获取整个介质粗糙面的电磁散射系数。本发明解决了实际大区域介质粗糙面电磁散射仿真的技术问题。用于国防和民用领域的粗糙面雷达散射特性研究。

Claims (4)

1.一种基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，包括如下步骤：

(1)输入介质粗糙面的功率谱密度函数和粗糙面参数，通过蒙特卡洛方法得到粗糙面上N个离散面元的位置坐标r_-N/2+1＝(x_-N/2+1，z_-N/2+1)，…，r_-1＝(x_-1，z_-1)，r₀＝(x₀，z₀)，r₁＝(x₁，z₁)，…，r_n/2＝(x_N/2，Z_N/2)；

(2)使用锥形波入射到该仿真的粗糙面上，将入射波强度大于e^-1/4的面元划分进主区域，入射波强度小于e^-1/4的面元划分进非主区域；

(3)在主区域中，根据电磁场中的边界条件，获得介质粗糙面主区域的矩量法边界积分方程，并将该方程转化为矩阵方程，求解该矩阵方程，得到介质粗糙面每个主区域面元上的等效电流；

(4)在非主区域的每个面元上，仅考虑入射磁场和主区域上等效电流在其表面的感应，得到每个非主区域面元上的等效电流；

(5)利用已得到的介质粗糙面每个面元上的等效电流，获取整个介质粗糙面的电磁散射系数。

2.根据权利要求1所述的基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，其特征在于：对介质粗糙面进行区域划分，包括如下步骤：

(2a)采用锥形波照射粗糙面，锥形波的入射电场E^inc(r)和入射磁场H^inc(r)分别表示如下：

$E^{\mathrm{inc}}\left(r\right)=\hat{y}\exp [-j{k}_0(x\sin {\mathrm{\θ}}_i-z\cos {\mathrm{\θ}}_i)(1+w\left(r\right))]\×\exp (-\frac{{(x+z\tan {\mathrm{\θ}}_i)}^2}{g^2})$

$H^{\mathrm{inc}}\left(r\right)=\hat{k}\×E^{\mathrm{inc}}\left(r\right)/Z_a$

其中 $w\left(r\right)=\exp [2{(x-\mathrm{yctg}{\mathrm{\θ}}_i)}^2/g^2-1]/{\left(\mathrm{kg}\sin {\mathrm{\θ}}_i\right)}^2,$

是空间波数的单位矢量，Z_a为上半空间波阻抗，

μ_a为上半空间相对磁介电常数，ε_a为上半空间相对介电常数，r＝(x，z)是面元的位置坐标，g是锥形入射波因子，θ_i为入射角；

(2b)依次读入步骤1中生成介质粗糙面上的每个面元的位置坐标r_-N/2+1＝(x_-N/2+1，z_-N/2+1)，…，r_-1＝(x_-1，z_-1)，r₀＝(x₀，z₀)，r₁＝(x₁，z₁)，…，r_N/2＝(x_N/2，Z_N/2)，将面元位置坐标带入步骤(2a)中，依次可得到每个面元上的入射电场E^inc(r)和入射磁场H^inc(r)，并且求得每个面元上电场的幅值|E^inc(r)|；

(2c)设置门限值为e^-1/4，将入射波电场幅值|E^inc(r)|大于门限值的面元划分为主区域，入射波电场幅值|E^inc(r)|小于门限值的面元划分为非主区域；

(2d)将位置坐标r的下标小于0的主区域面元划分为1l区域，由这些面元组成的区域用符号S_1l表示，将位置坐标r的下标大于等于0的主区域面元划分进1r区域，这些面元组成的区域用符号S_1r表示；

(2e)与步骤(2d)类似，将位置坐标r的下标小于0的非主区域面元划分为2l区域，由这些面元组成的区域用符号S_2l表示，将位置坐标r的下标大于等于0的非主区域面元划分进2r区域，这些面元组成的区域用符号S_2r表示。

3.根据权利要求2所述的基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，其特征在于：获得每个主区域面元上的等效电流J，包括如下步骤：

(3a)在主区域建立边界积分方程如下：

J可以代指主区域上每个面元上的等效电流，在主区域的每个面元上，该方程均成立，上式中，下标a和b分别表示粗糙面上方和下方，

是对应面元单位法向量，Z_a为上半空间波阻抗，Z_b为下半空间波阻抗，算子和分别定义为：

${\stackrel{\‾}{L}}_a\left(X\right)=-j{k}_a\underset{S}{\∫}[X+\frac{1}{k_a^2}\▿(\▿\·X)]G_a(r,r^{\′})d{s}^{\′}$

${\stackrel{\‾}{K}}_a\left(X\right)=-\underset{S}{\∫}X\×\▿G_a(r,r^{\′})d{s}^{\′}$

${\stackrel{\‾}{L}}_b\left(X\right)=-j{k}_b\underset{S}{\∫}[X+\frac{1}{k_b^2}\▿(\▿\·X)]G_b(r,r^{\′})d{s}^{\′}$

${\stackrel{\‾}{K}}_b\left(X\right)=-\underset{S}{\∫}X\×\▿G_b(r,r^{\′})d{s}^{\′}$

表示对整个介质粗糙面进行积分，G_a(r，r′)为粗糙面上半空间中的格林函数，G_b(r，r′)为粗糙面下半空间中的格林函数；

(3b)采用狄拉克函数作为基函数，采用伽辽金方法，将上式转化为矩阵方程：

$\left[\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\·\left[\begin{array}{c}{J}_{1l}\\ J_{1r}\end{array}\right]=\left[\begin{array}{c}{B}_1\\ B_2\end{array}\right]$

式中，J_1l表示一个向量，其中的每个元素对应S_1l区域中每个面元上的等效电流，J_1r表示一个向量，其中的每个元素对应S_1r区域中每个面元上的等效电流，A是子矩阵元素，B为入射波的列矩阵；

(3c)采用LU矩阵分解方法求解步骤(3b)中的矩阵方程，得到主区域的等效电流分布J_1l，J_1r。

4.根据权利要求3所述的基于多区域模型矩量法的介质粗糙面电磁散射仿真方法，其特征在于：获得每个非主区域面元上的等效电流J，包括如下步骤：

(4a)根据步骤(2b)得到的区域S_2l中任意一个面元上的入射磁场

通过下式：

$J_{2l}^{\mathrm{inc}}\left(r\right)=2\hat{n}\×H_{2l}^{\mathrm{inc}}\left(r\right)$

求得区域S_2l中任意一个面元上由入射场感应产生的等效电流

其中，r为该面元的位置坐标，

为面元的单位法向矢量；

(4b)根据步骤(3c)得到S_1l中每一个面元上的电流分布J_1l，S₂₁中任意一个面元上的感应得到的等效电流(r)可通过下式：

$J_{2l}^{S_{1l}}\left(r\right)=-2\hat{n}\×\underset{S_{1l}}{\mathrm{\Σ}}{J}_{1l}\left(r^{\′}\right)\×\▿G_a(r,r^{\′})\·\mathrm{\Δx}$

依次将S_1l中每个面元上的等效电流对S_2l中该面元的作用进行叠加，其中r′对应区域S_1l中个面元的位置坐标，G_a(r，r′)为粗糙面上半空间中的格林函数，

为该面元的单位法向矢量；

(4c)对S_2l每个面元执行步骤(4a)，(4b)，得到S_2l每个面元上由入射场感应产生的等效电流(r)和J_1l在S_2l上感应产生的等效电流

将这两项相加，得到S_2l中每个面元上的等效电流J_2l(r)，即 $J_{2l}\left(r\right)=J_{2l}^{\mathrm{inc}}\left(r\right)+J_{2l}^{S_{1l}}\left(r\right);$

(4d)根据步骤(2b)可得到区域S_2r中任意一个面元上的入射磁场

通过下式：

$J_{2r}^{\mathrm{inc}}\left(r\right)=2\hat{n}\×H_{2r}^{\mathrm{inc}}\left(r\right)$

求得区域S_2r中任意一个面元上由入射场感应产生的等效电流

其中，r为该面元的位置坐标，

为面元的单位法向矢量；

(4e)根据步骤(3c)得到S_1r中每一个面元上的电流分布J_1r，S_2r中任意一个面元上的感应得到的等效电流可通过下式：

$J_{2r}^{S_{1r}}\left(r\right)=-2\hat{n}\×\underset{S_{1r}}{\mathrm{\Σ}}{J}_{1r}\left(r^{\′}\right)\×\▿G_a(r,r^{\′})\·\mathrm{\Δx}$

依次将S_1l中每个面元上的等效电流对S_2l中该面元的作用进行叠加，其中r′对应区域S_1l中个面元的位置坐标，G_a(r，r′)为粗糙面上半空间中的格林函数，为该面元的单位法向矢量；

(4f)对S_2r每个面元执行步骤(4d)，(4e)，得到S_2r每个面元上由入射场感应产生的等效电流

和J_1r在S_2r上感应产生的等效电流

将这两项相加，得到S_2r中每个面元上的等效电流J_2r(r)，即 $J_{2r}\left(r\right)=J_{2r}^{\mathrm{inc}}\left(r\right)+J_{2r}^{S_{1r}}\left(r\right).$

Images