毫米波段内卫星电磁散射特性的数值方法

摘要

本发明公开了一种毫米波段内卫星电磁散射特性的数值方法。针对卫星这种金属介质混合结构，仅需对卫星金属部分建立电场积分方程，介质部分不需要建立方程，最终形成性态良好的方程，便于迭代求解。仅需对卫星金属表面进行网格剖分，不需要对涂敷介质部分进行网格剖分，奇异性容易处理。由于卫星斜长的结构特点，转移因子分量的方向性强，在多层快速多级子的基础上对远场作用采用了加速方法，有效的降低计算内存需求并且节省计算时间。

说明书

技术领域

本发明属于目标电磁散射特性的快速计算技术，特别是毫米波段内卫星电磁散射特性的数值方法。 

背景技术

卫星一般主要由金属材料构成，为了达到防腐蚀、装饰以及隐身等目的，通常会采用介质材料进行涂敷（QJ2754-95）。对于卫星这种由金属介质涂覆混合结构的电磁散射问题，学者们提出了一系列的电磁数值计算方法，其中基于积分方程的矩量法（MoM）（Ding-Feng Yu，Si-Yuan He，Hai-Tao Chen，Guo-Qiang Zhu，“Research on the electromagnetic scattering of3D target coated with anisotropic medium using impedance boundary condition”，Microwave and Optical Technology Letters，2011,53(2):458–462.）得到了最广泛的应用。 

针对这种金属介质混合结构，金属部分通常被作为理想导电体（PEC）来处理,并且容易被面积分方程方法（SIE）来分析求解，然而对于介质部分，不管是使用体积分方程方法，还是基于等效原理的面积分方程方法以及他们的演变算法，相对于纯金属目标的电磁计算分析方法，这些算法都需要更多的内存和时间资源。并且由于卫星涂敷介质部分非常薄的特点，如果选用面积分方程来分析介质部分，上下表面之间的奇异性较难处理；如果采用体积分分析介质部分，会导致剖分网格畸变，最终导致矩阵方程性态差，难于求解。 

毫米波波段内的卫星属于电大尺寸目标，即便采用了多层快速多级子技术就行加速计算，仍然面临计算资源消耗巨大的问题。在多层快速多级子算法中，随着层数的增加，转移矩阵的计算量会急剧增加。两个远场组之间的转移矩阵由球汉克尔函数、勒让德级数组成，转移因子被定义在有大量角谱分量的单位球上，即使采用插值和反插值方法，整个转移过程的计算也是十分耗时的。 

发明内容

本发明的目的在于提供一种毫米波段内卫星电磁散射特性的数值方法，从而实现快速得到电磁散射特性参数的方法。 

实现本发明目的的技术解决方案为：一种毫米波段内卫星电磁散射特性的数值方法，步骤如下： 

第一步，令均匀平面波照射到一个涂有薄介质层的卫星结构上，在卫星的表面将产生表面感应面电流J_S、面电荷ρ_S、薄介质内产生体极化电流J_pol、体极化磁流M_pol和极化电荷 ρ_S,pol，根据理想导体的电场边界条件，即金属表面的总场切向分量为0，得到涂有薄介质层的卫星结构目标的电场积分方程EFIE，如下 

[E^inc(r)+E^sca(r)]_tan=0    (1)其中，下标tan表示电场的切向分量，E^inc表示照射在目标上的电磁波，E^sca表示目标在电磁波照射后产生的散场，散射场的表达形式为： 

$E^{\mathrm{sca}}=-{\mathrm{j\omega \mu }}_0\underset{S}{\int }{J}_S\left(r^{\prime }\right)\bar{G}(r,r^{\prime }){\mathrm{dS}}^{\prime }-{\mathrm{j\omega \mu }}_0\underset{V}{\int }{J}_{\mathrm{pol}}\left(r^{\prime }\right)\bar{G}(r,r^{\prime }){\mathrm{dV}}^{\prime }+\underset{V}{\int }{M}_V\left(r^{\prime }\right)\times ▿G(r,r^{\prime }){\mathrm{dV}}^{\prime }$

$=-{\mathrm{j\omega \mu }}_0\underset{S}{\int }{J}_S\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{ϵ}_0}\underset{S}{\int }{\rho }_s\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }$

      (2) 

$-{\mathrm{j\omega \mu }}_0\underset{V}{\int }{J}_{\mathrm{pol}}\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{ϵ}_0}\underset{S}{\int }{\rho }_{s,\mathrm{pol}}\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{ϵ}_0}\underset{S^{\Delta }}{\int }{\rho }_{s,\mathrm{pol}}\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\Delta \prime }$

$+\underset{V}{\int }{M}_V\left(r^{\prime }\right)\times ▿G(r,r^{\prime }){\mathrm{dV}}^{\prime }$

其中,V表示薄介质涂层的体积单元，S表示金属表面单元即薄介质涂层的下表面单元，S^Δ表示薄介质涂层的上表面单元，ω为电磁波的角频率，μ₀和ε₀分别为真空中的磁导率以及介电参数，r和r'分别为场和源的位置坐标，G(r,r')为自由空间的格林函数，表达式为： 

$G(r,r^{\prime })=\frac{e^{-\mathrm{jk}|r-r^{\prime }|}}{4\pi |r-r^{\prime }|}---\left(3\right)$

其中k是自由空间的波数； 

第二步，将第一步中建立的电场积分方程中的5个未知量转换成金属表面电流密度J_S，5个未知量分别是金属表面电流密度J_S、金属表面电荷密度ρ_S、介质体极化电流密度J_pol、介质层的上下表面极化电荷密度ρ_S,pol和体极化磁流M_pol，其具体表示形式如下， 

式（1）最终变成： 

$[{\mathrm{j\omega \mu }}_0\underset{S}{\int }{J}_s\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{\mathrm{j\omega ϵ}}_0}\underset{S}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\prime }$

$-{\mathrm{j\omega \mu }}_0\kappa \underset{V}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)\hat{n}G(r,r^{\prime }){\mathrm{dS}}^{\prime }+\frac{\kappa }{{\mathrm{j\omega ϵ}}_0}\underset{S}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{\kappa }{{\mathrm{j\omega ϵ}}_0}\underset{S^{\Delta }}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\Delta \prime }$

$+{\mathrm{j\omega \mu }}_0({\mu }_r-1)\underset{V}{\int }{\hat{n}}^{\prime }\times J_S\left(r^{\prime }\right)\times ▿G(r,r^{\prime }){\mathrm{dV}}^{\prime }{]}_{\tan }=E_{\tan }^{\mathrm{inc}}---\left(5\right)$

第三步，对方程（5）进行矩量法求解，得到矩阵方程，采用多层快速多级子技术加速矩阵求解； 

第四步，当场源之间的距离满足一定的条件时，只转移场源连线方向上的角谱分量，将转移量降到最低，距离条件为： 

r_mn＞＞3γD    (6) 

γ是控制参量，γ越小，越多的非空组可以用快速远场近似来实现，相应结果精度的损失也就越大，远区组数随着γ逐渐增大而减少； 

第五步，求解矩阵方程，得到电流系数，再根据互易定理由电流系数计算电磁散射参量。 

所述第二步中未知量转换的方法如下： 

在交界面上由极化电流的法向不连续性可得 

$(J_{\mathrm{pol},1}-J_{\mathrm{pol},2})·\hat{n}=-{\mathrm{j\omega \rho }}_S---\left(7\right)$

目标模型中存在两种交界面，第1种：1为介质，2为金属，即介质下表面，即卫星与涂敷交界面，由于金属内J_pol,2=0，因此在介质下表面上可得 

$J_{\mathrm{pol},1}·\hat{n}=-{\mathrm{j\omega \rho }}_{S,\mathrm{pol}}---\left(8\right)$

第2种：1为空气，2为介质，即介质上表面，即涂敷与空气交界面，同样由于空气中J_pol,1=0，因此在介质上表面上可得 

$J_{\mathrm{pol},2}·\hat{n}=-{\mathrm{j\omega \rho }}_{S,\mathrm{pol}}---\left(9\right)$

为了将极化电流与卫星表面上的电流联系起来，利用麦克斯韦方程可得 

$J_{\mathrm{pol}}=\mathrm{j\omega }(ϵ-{ϵ}_0)E---\left(10\right)$

在卫星表面因此式（10）可以写成 

$J_{\mathrm{pol}}={\mathrm{j\omega \kappa \rho }}_S\hat{n}---\left(11\right)$

其中$\kappa =\frac{ϵ-{ϵ}_0}{ϵ};$

将电流连续性方程 

${\rho }_S=-\frac{1}{\mathrm{j\omega }}{▿}_S·J_S---\left(12\right)$

代入式（11）得， 

$J_{\mathrm{pol}}=-\kappa (▿·J_S)\hat{n}---\left(13\right)$

结合式（8），（11）和式（12），在介质下表面可得 

${\rho }_{S,\mathrm{pol}}=-\frac{\kappa }{\mathrm{j\omega }}▿·J_S---\left(14\right)$

同理在介质上表面可得 

${\rho }_{S,\mathrm{pol}}=\frac{\kappa }{\mathrm{j\omega }}▿·J_S---\left(15\right)$

忽略法向磁场，同时切向磁场是一个常数，根据理想导体的磁场边界条件得到以下表达式： 

$M_{\mathrm{pol}}\left(r^{\prime }\right)\approx \mathrm{j\omega }(\mu -{\mu }_0)H_t\left(r^{\prime }\right)=\mathrm{j\omega }\frac{({\mu }_r-1)}{{\mu }_r}{B}_t\left(r^{\prime }\right)---\left(16\right)$

μ_r为薄介质的相对磁导率，将B_t(r')做变形： 

$B_t\left(r^{\prime }\right)=-{\hat{n}}^{\prime }\times {\hat{n}}^{\prime }\times B_t\left(r^{\prime }\right)=-\mu {\hat{n}}^{\prime }\times {\hat{n}}^{\prime }\times H_t\left(r^{\prime }\right)=-\mu {\hat{n}}^{\prime }\times J_S\left(r^{\prime }\right)---\left(17\right)$

将（17）式代入到（16）中得到M_pol用金属表面的感应面电流J_s表示的表达式如下： 

$M_{\mathrm{pol}}\left(r^{\prime }\right)=-{\mathrm{j\omega \mu }}_0({\mu }_r-1){\hat{n}}^{\prime }\times J_S\left(r^{\prime }\right)---\left(18\right)$

将（12）、（13）、（14）、（15）、（18）代入（2）中得到： 

$[{\mathrm{j\omega \mu }}_0\underset{S}{\int }{J}_s\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{\mathrm{j\omega ϵ}}_0}\underset{S}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\prime }$

$-{\mathrm{j\omega \mu }}_0\kappa \underset{V}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)\hat{n}G(r,r^{\prime }){\mathrm{dS}}^{\prime }+\frac{\kappa }{{\mathrm{j\omega ϵ}}_0}\underset{S}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{\kappa }{{\mathrm{j\omega ϵ}}_0}\underset{S^{\Delta }}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\Delta \prime }$

$+{\mathrm{j\omega \mu }}_0({\mu }_r-1)\underset{V}{\int }{\hat{n}}^{\prime }\times J_S\left(r^{\prime }\right)\times ▿G(r,r^{\prime }){\mathrm{dV}}^{\prime }{]}_{\tan }=E_{\tan }^{\mathrm{inc}}---\left(5\right)$

所述第三步中采用多层快速多级子技术加速矩阵求解的方法如下： 

步骤3.1，建立卫星的树形结构； 

步骤3.2，假设在给定的一层中，任一组m中有一个场点r_i，任一组n中有一个源点r_j，r_m、r_n分别是场点组和源点组的中心点；两个非空组中场点和源点的空间矢量记为： r_ij=r_i-r_j=r_im+r_mn+r_nj，当所分析的场点组和源点组不重合也不相邻时，满足|r_im+r_nj|<|r_mn|，自由空间的标量格林函数可在角谱空间写成： 

$\frac{e^{-\mathrm{jk}|r_i-r_j|}}{|r_i-r_j|}=\frac{\mathrm{jk}}{4\pi }\underset{S_E}{\int }{e}^{-\mathrm{jk}·(r_{\mathrm{im}}+r_{\mathrm{nj}})}{\alpha }_{\mathrm{mn}}(\hat{k}·{\hat{r}}_{\mathrm{mn}})---\left(19\right)$

整个积分是定义在单位球S_E上的，α_mn是两个组的转移因子，定义如下： 

${\alpha }_{\mathrm{mn}}=\underset{l=0}{\overset{L}{\Sigma }}{(-j)}^l(2l+1)h_l^{\left(2\right)}\left(\mathrm{kr}\right)P_l(\hat{k}·\hat{r})---\left(20\right)$

其中，j_l(kd)为球贝塞函数，为第二类球汉克尔函数，为勒让德函数L是无限级数的截断长度，L＝kD+β(kD)^1/3，D是分组的尺寸，β为精度参数，β≥2；用（θ，φ）表示单位球S_E坐标，积分点数为K_L＝2L²，其中在θ方向共采集L个点的一维高斯积分，φ方向共采集2L个点数的梯形法则积分； 

步骤3.3，当两个远场作用组中心距离和组的大小确定时，将多个转移因子构成的转移矩阵加上一个窗函数： 

$w_l=\left(\begin{array}{cc}1& l\le L/2\\ \sin \left(\frac{1}{L}\pi \right)& l>L/2\end{array}\right)---\left(21\right)$

在转移因子上乘上窗函数使得两组中心连线方向附近的转移分量的值相对于其它部分转移因子的值有陡峭的变化，省略对转移过程影响较小的角谱分量的作用。 

所述步骤3.1中建立卫星的树形结构的方法如下：将卫星放入到一个能够包围它的立方体空间中，记该空间坐标系为XYZ，该空间记为树形结构的第零层，然后对该空间三维方向上分别一分为二，形成8个子空间，记为第一层，在8个子空间的基础上继续细分，当最细层空间的尺寸小于等于设定的电波长时则停止细分；整个卫星被转到一个树状结构上，对每一层的空间进行编码建立索引，保存每一层非空组空间的个数。 

所述步骤3.1中设定的波长在0.2个波长至0.5个波长之间。 

所述步骤3.3中判定对转移过程影响较小的角谱分量的方法如下：设θ_r为两作用组连线方向与空间z轴的夹角，φ_r为两作用组连线方向向量在XOY面上投影和x轴的夹角，θ_e为设定的阈值对应的角谱分量方向与两作用组连线方向的夹角；在计算转移因子的过程中，转移因子模值的大小随着θ增大而减小，提取并计算角谱分量方向与两作用组中心连线方向夹 角不大于θ_e的角谱分量。 

本发明与现有技术相比，其显著优点：（1）方程建立简单。针对卫星这种金属介质混合结构，仅需对金属部分建立电场积分方程，介质部分不需要建立方程进行描述。（2）所需计算资源少。由于卫星斜长的结构特点，转移因子分量的方向性强。而射线追踪方法中由于只需要转移场源连线方向周围的角谱分量，相较于多层快速多级子中对单位球上的所有角谱分量进行转移，节省了转移因子的计算时间和所需的存储内存。同时当满足快速远场近似的条件下，仅需计算转移因子在场源组中心连线方向上的角谱分量，在满足计算精度的条件下进一步节省了计算资源。 

附图说明

图1是本发明金属涂覆薄介质层目标示意图，a是几何模型示意图，b是电磁参数模型结构示意图。 

图2是本发明边界极化电流的法向分量的不连续性示意图。 

图3是本发明多层快速多极子远场作用示意图。 

图4是本发明卫星模型的结构及涂覆示意图 

图5是本发明双站RCS曲线示意图。 

具体实施方式

本发明毫米波段内卫星电磁散射特性的数值方法，步骤如下： 

第一步，电磁散射积分方程的建立，即基于理想导体的电场边界条件，在金属表面的总场切向分量为0，而总场即为入射电场与散射电场之和。入射电场为已知激励，均匀平面波通常被用来作为入射电场，散射电场可以用待求的表面未知电流来表示。 

第二步，未知量之间的转换，即第一步中建立的电场积分方程中共有5个未知量，分别是金属表面电流密度、金属表面电荷密度，介质体极化电流密度，介质层的上下表面极化电荷密度。在采用矩量法来求解该方程的过程中，由于仅有一个测试项：金属表面电流密度，显然是不能求解出这个方程正确的结果，因此必须将所有的未知源都用金属表面电流密度表示； 

第三步，结合卫星结构的特点，利用射线追踪的基本原理，合理的选择出转移因子中角谱分量大的部分，将角谱分量小的部分舍弃掉。 

第四步，结合卫星结构的特点，在更高层采用快速远场近似的方法快速计算远场作用。 

第五步，矩阵方程求解以及电磁散射参数的计算。 

下面结合附图对本发明作进一步详细描述。 

结合图1，本发明基于一种分析毫米波段内卫星电磁散射特性的数值方法的数值方法，步骤如下： 

第一步，令均匀平面波照射到一个涂有薄介质层的卫星结构上，在卫星的表面将产生表面感应面电流J_S和面电荷ρ_S，薄介质内产生体极化电流J_pol、体极化磁流M_pol、极化电荷ρ_S,pol，根据理想导体的电场边界条件，即金属表面的总场切向分量为0，得到涂有薄介质层的卫星结构目标的电场积分方程（EFIE），如下 

[E^inc(r)+E^sca(r)]_tan=0    (1)其中，下标tan表示电场的切向分量，E^inc表示照射在目标上的电磁波，E^sca表示目标在电磁波照射后产生的散场，而散射场包括金属和薄介质两部分的作用，散射场的表达形式为： 

$E^{\mathrm{sca}}=-{\mathrm{j\omega \mu }}_0\underset{S}{\int }{J}_S\left(r^{\prime }\right)\bar{G}(r,r^{\prime }){\mathrm{dS}}^{\prime }-{\mathrm{j\omega \mu }}_0\underset{V}{\int }{J}_{\mathrm{pol}}\left(r^{\prime }\right)\bar{G}(r,r^{\prime }){\mathrm{dV}}^{\prime }+\underset{V}{\int }{M}_V\left(r^{\prime }\right)\times ▿G(r,r^{\prime }){\mathrm{dV}}^{\prime }$

$=-{\mathrm{j\omega \mu }}_0\underset{S}{\int }{J}_S\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{ϵ}_0}\underset{S}{\int }{\rho }_s\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }$

           (2) 

$-{\mathrm{j\omega \mu }}_0\underset{V}{\int }{J}_{\mathrm{pol}}\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{ϵ}_0}\underset{S}{\int }{\rho }_{s,\mathrm{pol}}\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{ϵ}_0}\underset{S^{\Delta }}{\int }{\rho }_{s,\mathrm{pol}}\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\Delta \prime }$

$+\underset{V}{\int }{M}_V\left(r^{\prime }\right)\times ▿G(r,r^{\prime }){\mathrm{dV}}^{\prime }$

其中V表示薄介质涂层的体积单元，S表示金属表面单元（即薄介质涂层的下表面单元），S^Δ表示薄介质涂层的上表面单元，ω为电磁波的角频率，μ₀和ε₀分别为真空中的磁导率以及介电参数，r和r'分别为场和源的位置坐标，G(r,r')为自由空间的格林函数，表达式为： 

$G(r,r^{\prime })=\frac{e^{-\mathrm{jk}|r-r^{\prime }|}}{4\pi |r-r^{\prime }|}---\left(3\right)$

其中k是自由空间的波数。 

第二步，第一步中建立的电场积分方程中共有5个未知量，分别是金属表面电流密度J_S、金属表面电荷密度ρ_S，介质体极化电流密度J_pol，介质层的上下表面极化电荷密度ρ_S,pol，体极化磁流M_pol，如果采用矩量法来求解该方程，仅有一个测试项：金属表面电流密度J_S，显然是不能求解出这个方程正确的结果，因此必须将所有的未知源都用来金属表面电流密度表示。 

边界极化电流的法向分量的不连续性如图2所示。在卫星与涂敷交界面、涂敷与空气交 界面上由极化电流的法向不连续性可得 

$(J_{\mathrm{pol},1}-J_{\mathrm{pol},2})·\hat{n}=-{\mathrm{j\omega \rho }}_S---\left(4\right)$

本发明中目标模型中存在两种交界面，第1种：1为介质，2为金属，即介质下表面，由于金属内J_pol,2=0，因此在介质下表面上可得 

$J_{\mathrm{pol},1}·\hat{n}=-{\mathrm{j\omega \rho }}_{S,\mathrm{pol}}---\left(5\right)$

第2种：1为空气，2为介质，即介质上表面，同样由于空气中J_pol,1=0，因此在介质上表面上可得 

$J_{\mathrm{pol},2}·\hat{n}=-{\mathrm{j\omega \rho }}_{S,\mathrm{pol}}---\left(6\right)$

为了将极化电流与卫星表面上的电流联系起来，利用麦克斯韦方程可得 

J_pol＝jω(ε-ε₀)E    (7)在卫星表面因此式（7）可以写成 

$J_{\mathrm{pol}}=-{\mathrm{j\omega \kappa \rho }}_S\hat{n}---\left(8\right)$

其中$\kappa =\frac{ϵ-{ϵ}_0}{ϵ}.$

将电流连续性方程 

${\rho }_S=-\frac{1}{\mathrm{j\omega }}{▿}_S·J_S---\left(9\right)$

代入式（8）得， 

$J_{\mathrm{pol}}=-\kappa (▿·J_S)\hat{n}---\left(10\right)$

由式（9）和式（10）可以看出金属上的电荷密度和介质体内的极化电流密度均已经转变成了金属面电流密度的形式表示； 

结合式（5），（8）和式（9），在介质下表面可得 

${\rho }_{S,\mathrm{pol}}=-\frac{\kappa }{\mathrm{j\omega }}▿·J_S---\left(11\right)$

同理在介质上表面即卫星涂敷与空气交界面可得 

${\rho }_{S,\mathrm{pol}}=\frac{\kappa }{\mathrm{j\omega }}▿·J_S---\left(12\right)$

考虑法向磁场是可以忽略，同时切向磁场可以认为是一个常数，所以根据理想导体的磁场边界条件可以得到以下表达式： 

$M_{\mathrm{pol}}\left(r^{\prime }\right)\approx \mathrm{j\omega }(\mu -{\mu }_0)H_t\left(r^{\prime }\right)=\mathrm{j\omega }\frac{({\mu }_r-1)}{{\mu }_r}{B}_t\left(r^{\prime }\right)---\left(13\right)$

μ_r为薄介质的相对磁导率。将B_t(r')做简单的变形： 

$B_t\left(r^{\prime }\right)=-{\hat{n}}^{\prime }\times {\hat{n}}^{\prime }\times B_t\left(r^{\prime }\right)=-\mu {\hat{n}}^{\prime }\times {\hat{n}}^{\prime }\times H_t\left(r^{\prime }\right)=-\mu {\hat{n}}^{\prime }\times J_S\left(r^{\prime }\right)$

将（14）式代入到（13）中可以得到M_pol用金属表面的感应面电流J_s表示的表达式如下： 

$M_{\mathrm{pol}}\left(r^{\prime }\right)=-{\mathrm{j\omega \mu }}_0({\mu }_r-1){\hat{n}}^{\prime }\times J_S\left(r^{\prime }\right)---\left(15\right)$

将（9）、（10）、（111）、（12）、（15）代入（2）中可以得到： 

$[{\mathrm{j\omega \mu }}_0\underset{S}{\int }{J}_s\left(r^{\prime }\right)G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{1}{{\mathrm{j\omega ϵ}}_0}\underset{S}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\prime }$

$-{\mathrm{j\omega \mu }}_0\kappa \underset{V}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)\hat{n}G(r,r^{\prime }){\mathrm{dS}}^{\prime }+\frac{\kappa }{{\mathrm{j\omega ϵ}}_0}\underset{S}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\prime }-\frac{\kappa }{{\mathrm{j\omega ϵ}}_0}\underset{S^{\Delta }}{\int }{▿}^{\prime }·J_s\left(r^{\prime }\right)▿G(r,r^{\prime }){\mathrm{dS}}^{\Delta \prime }$

$+{\mathrm{j\omega \mu }}_0({\mu }_r-1)\underset{V}{\int }{\hat{n}}^{\prime }\times J_S\left(r^{\prime }\right)\times ▿G(r,r^{\prime }){\mathrm{dV}}^{\prime }{]}_{\tan }=E_{\tan }^{\mathrm{inc}}---\left(16\right)$

第三步，对方程（16）进行矩量法求解，得到矩阵方程。对自由空间的格林函数按照加法定理展开，利用射线追踪的基本原理，结合体面积分方程的特性合理的选择出转移因子中角谱分量大的部分，将角谱分量小的部分舍弃掉。 

毫米波段卫星属于电大尺寸目标，多层快速多级子技术被用来加速矩阵求解。多层快速多级子技术首先需要建立树形结构：将卫星放入到一个能够包围它的立方体空间中，记该空间坐标系为XYZ，该空间记为树形结构的第零层，然后对该空间三维方向上分别一分为二，形成8个子空间，记为第一层，以为类推细分，直到最细层空间的尺寸小于等于一定的电波长（0.2个波长至0.5个波长之间）则停止细分。整个卫星就被转到一个树状结构上，对每一层的空间进行编码建立索引，保存每一层非空组（空间）的个数，所有的计算都在非空组中进行。 

假设在给定的一层中，组m中有一个场点r_i和组n中有一个源点r_j，r_m、r_n分别是场点组和源点组的中心点。两个非空组中场点和源点的空间矢量记为：r_ij=r_i-r_j=r_im+r_mn+r_nj，当所分析的场点组和源点组不重合也不相邻时，满足|r_im+r_nj|<|r_mn|，自由空间的标量格林函数可以在角谱空间写成： 

$\frac{e^{-\mathrm{jk}|r_i-r_j|}}{|r_i-r_j|}=\frac{\mathrm{jk}}{4\pi }\underset{S_E}{\int }{e}^{-\mathrm{jk}·(r_{\mathrm{im}}+r_{\mathrm{nj}})}{\alpha }_{\mathrm{mn}}(\hat{k}·{\hat{r}}_{\mathrm{mn}})---\left(17\right)$

整个积分是定义在单位球S_E上的，α_mn是两个组的转移因子，定义如下： 

${\alpha }_{\mathrm{mn}}=\underset{l=0}{\overset{L}{\Sigma }}{(-j)}^l(2l+1)h_l^{\left(2\right)}{\left(\mathrm{kr}\right)P}_l(\hat{k}·\hat{r})---\left(18\right)$

其中，j_l(kd)为球贝塞函数，为第二类球汉克尔函数，为勒让德函数。L是无限级数的截断长度，L＝kD+β(kD)^1/3，D是分组的尺寸，β为经验参数，本发明中取2。用（θ，φ）表示单位球坐标，积分点数为K_L＝2L²，其中在θ方向共采集L个点的一维高斯积分，φ方向共采集2L个点数的梯形法则积分。 

当两个远场作用组中心距离和组的大小确定时，转移因子的特性类似于线性天线阵列的辐射方向图特性，将转移矩阵加上一个窗函数： 

$w_l=\left(\begin{array}{cc}1& l\le L/2\\ \sin \left(\frac{1}{L}\pi \right)& l>L/2\end{array}\right)---\left(19\right)$

其中，L为球坐标下θ方向L个离散点的个数，l为每个转移因子分量的编号。具体的做法是：在转移因子上乘上窗函数的作用，使得两组中心连线方向附近的转移分量的值相对于其它部分转移因子的值有陡峭的变化，在这种情况下，可以省略对转移过程影响较小的角谱分量的作用，从而使得射线多极子获得更高的效率。同时，可以验证这种处理对计算精度影响很小。接着需要设定需要转移的角谱分量的范围。如图3所示，θ_r为两作用组连线方向与z轴的夹角，φ_r为两作用组连线方向向量在XOY面上投影和x轴的夹角，θ_e为设定的阈值对应的角谱分量方向与两作用组连线方向的夹角。在计算转移因子的过程中，转移因子角谱分量的大小随着与远场作用组中心连线方向的夹角的增大而减小。这样可以对转移因子设定阈值对圆锥区域的大小进行判定。在具体操作中，只需提取并计算角谱分量方向与两作用组中心连线方向夹角不大于θ_e的角谱分量。在算法验证时可以得到这部分角谱分量可以很好地联系两个远场作用组的作用。 

第四步，结合体面积分方程的特性，在更高层采用快速远场近似的方法快速计算远场作用。快速远场近似方法中，当场源之间的距离满足一定的条件时，我们可以只计算场源连线方向上的角谱分量，这样将转移量降到最低。本发明中采用组与组之间采用快速远场近似的判据条件为： 

r_mn＞3γD    (20) 

其中D为该层组的电尺寸，γ为结合体面积分方程特点设置的一个经验值，本发明中取 为1.5。γ越小，越多的非空组可以用快速远场近似来实现，相应结果精度的损失也就会越大，远区组数随着γ逐渐增大而减少。由此可以将快速远场近似、射线追踪和多层快速多级子方法相结合，基于远场组之间的距离可以将每一层组与组之间的作用采用三种快速方法进行处理：1快速多极子加速远场作用；2射线多极子加速远场作用；3快速远场近似加速远场作用。 

第五步，求解矩阵方程，得到电流系数，再根据互易定理由电流系数计算电磁散射参量。 

为了验证本发明方法的正确性与有效性，下面给出了计算毫米波段内空间卫星雷达散射截面的示例，并且计算结果均与传统算法进行了比较，吻合得很好。 

空间卫星的组合结构由两部分组成，中间有一个边长为2m立方体星体，在xoy上有两个长度为10.87m的太阳能电池板。星体的上下前后四面涂覆介质为（3.5，-0.008）。太阳能电池板的涂覆介质为（11.9,0.0）。入射波频率为35GHz,入射波的方向θ=0°→180°,φ=0°。剖分尺寸为0.14波长，未知量个数为125160048，最细层组尺寸：0.2波长，层数为12。从图4可以看出，加入快速远场近似和射线多极子后，与传统方法相比，RCS结果基本吻合，满足精度要求。而在表1中可以看出，加入了快速远场近似和射线多极子后，在LC=8（LC是使用快速远场近似的层数）时，程序的计算时间比原来的程序节省了49.6%；内存节省了52.8%。因此快速远场近似和射线多极子能够节省计算时间和内存。 

表1计算内存与时间统计表。