-
雷达目标的电磁散射特性仿真可为目标识别、武器制导和引信控制等武器装备研制与验证试验提供关键的数据支撑[1]。以导弹类目标为例,在靶场开展飞行试验之前,需获取目标的电磁散射特性数据,以检验目标是否满足考核指标要求,或为试验方案设计中的配试装备布站提供参考依据。以往的靶场飞行验证试验方案设计中,目标电磁散射特性数据主要依赖工业部门提供,或依据工程积累得到的经验公式[2]进行估算。工业部门提供的数据精度较高但计算周期较长,当目标需通过多次调整其散射特性以满足试验考核要求时,会由于计算量巨大而对试验进度造成严重影响。近年来,随着靶场飞行试验贴近实战化程度的提升,如何高效地进行目标电磁散射特性的快速估算,已成为一个亟待解决的重要课题。
众所周知,任何雷达目标散射特性的定量获取从原理上讲,最终都可归结为求解满足各种边界条件的麦克斯韦 (Maxwell) 方程组[3-5]的解。然而当前的仿真分析软件,无论是基于微分方程方法、积分方程方法或高频近似方法,都需要工程操作人员具有较强的电磁数理基础,以完成雷达目标的复杂几何建模和网格剖分,不利于电磁仿真软件的应用推广。同时,即使是采用高频近似方法,无论是基于面元离散[6]还是基于NUBUS曲面建模[7],子域求和的积分策略使得电大尺寸 (数十至数千波长) 目标的单站散射计算仍然需要耗费一定的计算时间,导致大口径、高分辨率的一维或多维散射特征提取仍然非常困难,难以满足方案设计和论证等工程应用的实时计算需求。
针对靶场飞行试验方案设计中对导弹目标电磁散射特性的实时计算需求,本文将常规导弹目标分解为平板、圆柱、圆台 (或圆锥) 和椭球冠等基本几何体散射单元,首先采用物理光学方法[8]和边缘绕射理论[9]对各基本单元的电磁散射特性进行解析计算,再通过参数化组合建模、散射场解析计算和矢量场叠加的方法,快速仿真得到导弹目标的电磁散射数据。该方法可实现任意导弹类目标基于合理参数的快速组合建模与高效电磁仿真计算。无需复杂的几何建模和网格离散,设计的仿真软件可近乎实时地仿真获取导弹目标的RCS、一维距离像和二维ISAR像等电磁散射特性。
-
当目标的电尺寸较大时,散射变成了一种局部现象,散射体各部分之间的相互影响很小,这样就可把整个散射体作为多个独立散射单元的集合来处理。对于导弹目标,可以将其近似为平板、圆柱、圆台 (或圆锥) 和椭球冠等基本几何体散射单元的空间组合。因此整个弹体的电磁散射可等效为这些组合单元的电磁散射叠加。
金属目标的散射场由其表面感应电流的二次辐射引起,该散射贡献可由物理光学法近似计算。物理光学法 ((physical optical, PO) 的数理基础是斯特拉顿—朱兰成 (Stratton-Chu) 散射场积分方程[10]。即由金属表面的磁场边界条件:
$$ \hat{n}\times \boldsymbol{H}(\boldsymbol{r})=\hat{n}\times \left[ \boldsymbol{H}_{{}}^{\text{inc}}(\boldsymbol{r})+\boldsymbol{H}_{{}}^{\text{sca}}(\boldsymbol{r}) \right]=\boldsymbol{J}(\boldsymbol{r}) $$ (1) 以及场-源积分表达式:
$$ \boldsymbol{H}_{{}}^{\text{sca}}(\boldsymbol{r})=\int_{s}{\nabla G(\boldsymbol{r},{\boldsymbol{r}}')\times \boldsymbol{J}({\boldsymbol{r}}')\text{d}{\boldsymbol{r}}'} $$ (2) 可以建立金属表面感应电流满足的磁场积分方程 (MFIE) 为:
$$ \frac{1}{2}\boldsymbol{J}(\boldsymbol{r})-\hat{n}\times PV\int_{s}{\nabla G(\boldsymbol{r},{\boldsymbol{r}}')\times \boldsymbol{J}({\boldsymbol{r}}')\text{d}{\boldsymbol{r}}'}=\hat{n}\times \boldsymbol{H}_{{}}^{\text{inc}}(\boldsymbol{r}) $$ (3) 式中,左端第二项表示目标表面除场点以外其他各处源点r'(r'≠r) 对场点r的电磁耦合贡献;G (r, r') 为均匀背景空间的格林函数。在入射波能够直接到达的目标区域 (亮区),根据高频场的局部性原理,可忽略目标其他部分的耦合贡献,即近似认为式 (3) 中的积分贡献为零,则式 (3) 可写为:
$$ \boldsymbol{J}(\boldsymbol{r})=2\hat{n}\times \boldsymbol{H}_{{}}^{\text{inc}}(\boldsymbol{r})\text{ }\boldsymbol{r}\in 亮区 $$ (4) 式中,Hinc为入射波的磁场分量,可表示为:
$$ \boldsymbol{H}_{{}}^{\text{inc}}(\boldsymbol{r})=({{\boldsymbol{\hat{k}}}_{i}}\times {{\boldsymbol{\hat{e}}}_{i}})\frac{{{E}_{0}}}{\eta }{{\text{e}}^{-\text{j}k{{{\hat{k}}}_{i}}\cdot r}} $$ (5) 式中,${{\boldsymbol{\hat{k}}}_{i}}$为电波传播方向矢量;${{\boldsymbol{\hat{e}}}_{i}}$为入射波极化方向矢量。
而在入射波不能直接照射的阴影区 (暗区),其他各源点耦合贡献和入射波近乎抵消,使得暗区的电流强度很小。基于此,物理光学法直接假定在物体暗区电流密度为零:
$$ \boldsymbol{J}(\boldsymbol{r})=0\text{ }\boldsymbol{r}\in 暗区 $$ (6) 因此,由式 (4) 和式 (6),可根据入射场独立地近似确定目标表面感应电流分布。可见,用PO方法计算物体的散射场时,必须对物体进行面元的遮挡判别。对于关注的导弹目标,亮、暗区的判别则显得比较简单,在基本散射单元反射贡献求解和参数组合建模中再分别进行说明。
近似得到目标表面的感应电流后,由场-源积分可求解得到目标的散射场。其远区散射电磁场表达式为:
$$ \begin{array}{l} {\boldsymbol{E}^{{\rm{sca}}}}(\boldsymbol{r},{{\boldsymbol{\hat k}}_s}) = - {\rm{j}}k\eta \int_s {\boldsymbol{\bar G}(\boldsymbol{r},\boldsymbol{r}') \cdot \boldsymbol{J}(\boldsymbol{r}'){\rm{d}}\boldsymbol{r}'} \approx \\ - {\rm{j}}k{E_0}\frac{{{{\rm{e}}^{ - {\rm{j}}k{\rm{r}}}}}}{{2{\rm{\pi }}\boldsymbol{r}}}\int_S {\boldsymbol{\hat n}' \times ({{\boldsymbol{\hat k}}_i} \times {{\boldsymbol{\hat e}}_i}){{\rm{e}}^{{\rm{j}}2k{{\hat k}_s} \cdot r'}}{\rm{d}}s'} \end{array} $$ (7) 式中,${\hat k_s}$为散射方向矢量。${\hat k_s} = - {\hat k_i}$即对应单站散射。
为了便于散射贡献的量化分析,可用复雷达散射截面 (复RCS)[11]表示散射场的主极化分量:
$$ \sqrt \sigma = \mathop {\lim }\limits_{r \to \infty } 2\sqrt {\rm{\pi }} \boldsymbol{r}\frac{{{\boldsymbol{E}^{{\rm{sca}}}} \cdot {{\boldsymbol{\hat e}}_i}}}{{{E_0}}}{{\rm{e}}^{{\rm{j}}kr}} $$ (8) 将式 (7) 代入式 (8),得:
$$ \sqrt \sigma = - {\rm{j}}\frac{{kI}}{\pi } $$ (9) 式中,积分核I定义为:
$$ I = \int_S {\boldsymbol{\hat n}'} \cdot {\boldsymbol{\hat k}_s}{{\rm{e}}^{{\rm{j}}2k{{\hat k}_s} \cdot r'}}{\rm{d}}S' $$ (10) 复RCS体现了散射的相位和接收机极化的影响,是一个复数量,在计算中可保持目标上各部件散射场之间的相位关系,正确地表示干涉作用。因此不同散射单元的复RCS可以直接相加,以表示不同部件总的散射效果。
目标总的RCS可以表示为:
$$ \sigma = {\left| {\sum\limits_n {{{\sqrt \sigma }_n}} } \right|^2} $$ (11) 式中,∑表示对所有组合单元的散射贡献进行求和叠加。
-
由上节的分析可知,基于物理光学近似,雷达目标的散射场可由式 (10) 进行计算。本节对平板、圆柱、圆台和椭球等基本散射单元的积分计算进行解析公式推导。式 (10) 中,${\hat k_s}$为散射波方向矢量,对于$\left( {\theta ,\varphi } \right)$方向的雷达来波,单站情形的散射方向矢量可表示为:
$$ {\boldsymbol{\hat k}_s} = (\sin \theta \cos \phi ,\sin \theta \sin \phi ,\cos \theta ) $$ (12) 可见,积分核仅与散射单元的外形 ($\boldsymbol{\hat n}'$,r'和dS') 密切相关。
-
如图 1所示,平行于XOY平面的矩形平板长和宽分别为Lx和Ly,平板中心位于点 (0, 0, z0) 处。矩形平板上的源点和法向矢量可分别表示为:$r' = (x',y',{z_0})$,$\boldsymbol{\hat n}' = (0,0,1)$,结合式 (10),散射积分可表示为:
$$ \begin{array}{l} I = \cos \theta \int_{ - \frac{{{L_x}}}{2}}^{\frac{{{L_x}}}{2}} {\int_{ - \frac{{{L_y}}}{2}}^{\frac{{{L_y}}}{2}} {{{\rm{e}}^{{\rm{j}}2k(x'\sin \theta \cos \phi + y'\sin \theta \sin \phi + {z_0}\cos \theta )}}{\rm{d}}x'{\rm{d}}y} } = \\ {{\rm{e}}^{{\rm{j2}}k{z_{\rm{0}}}\cos \theta }}\frac{{\sin (k\sin \theta \cos \phi {L_x})}}{{k\sin \theta \cos \phi }}\frac{{\sin (k\sin \theta \sin \phi {L_y})}}{{k\sin \theta \sin \phi }}\cos \theta \end{array} $$ (13) 当分母为零时,可用洛必达 (L'Hospital) 法则[12]求极限化简式 (13),求得特殊角度时的散射结果。
-
如图 2所示,平行于XOY平面的圆形平板半径为a,平板中心位于点 (0, 0, z0) 处。圆形平板上的源点和法向矢量可分别表示为$r' = (\rho '\cos \phi ',\rho '\sin \phi ',{z_0})$,$n' = (0,0,1)$,结合式 (10),散射积分可表示为:
$$ \begin{array}{l} I = \cos \theta \int_{{\rm{ }}0}^{{\rm{ }}a} {\int_{{\rm{ }} - {\rm{\pi }}}^{{\rm{ \pi }}} {\rho '{{\rm{e}}^{{\rm{j}}2k(\sin \theta \cos \phi \rho '\cos \phi ' + \sin \theta \sin \phi \rho '\sin + {z_0}\cos \theta )}}{\rm{d}}\rho '{\rm{d}}\phi '} } = \\ \;\;\;\;\;\;\cos \theta {{\rm{e}}^{{\rm{j}}2k{z_0}\cos \theta }}\int_{{\rm{ }}0}^{{\rm{ }}a} {\int_{{\rm{ }} - {\rm{\pi }}}^{{\rm{ \pi }}} {\rho '{{\rm{e}}^{{\rm{j}}2k\sin \theta \rho '\cos (\phi ' - \phi )}}{\rm{d}}\rho '{\rm{d}}\phi '} } \end{array} $$ (14) 为了求取式 (14) 中的内层积分,利用贝塞尔函数恒等式:
$$ {{\rm{e}}^{{\rm{j}}2k\sin \theta \rho '\cos (\phi ' - \phi )}}{\rm{ = }}\sum\limits_{m = - \infty }^\infty {{J_m}(2k\rho '\sin \theta ){{\rm{e}}^{{\rm{j}}m(\phi ' + \phi )}}} $$ (15) 将式 (15) 代入式 (14),得到:
$$ \begin{array}{c} I = \cos \theta {{\rm{e}}^{{\rm{j}}2k{z_0}\cos \theta }}\sum\limits_{m = - \infty }^\infty {\int_{{\rm{ }}0}^{{\rm{ }}a} {\rho '{J_m}(2k\rho '\sin \theta ){\rm{d}}\rho '} } \times \\ \int_{{\rm{ }} - {\rm{\pi }}}^{{\rm{ \pi }}} {{{\rm{e}}^{{\rm{j}}m(\phi ' + \phi )}}{\rm{d}}\phi '} = 2{\rm{\pi cos}}\theta {{\rm{e}}^{{\rm{j}}2k{z_0}\cos \theta }}\int_{{\rm{ }}0}^{{\rm{ }}a} {\rho '{J_0}(2k\rho '\sin \theta ){\rm{d}}\rho '} = \\ 2{\rm{\pi }}{a^{\rm{2}}}{\rm{cos}}\theta {{\rm{e}}^{{\rm{j}}2k{z_0}\cos \theta }}\frac{{{J_1}(2ka\sin \theta )}}{{2ka\sin \theta }} \end{array} $$ (16) 当$\sin \theta = 0$时,$\frac{{{J}_{1}}(2ka\sin \theta )}{2ka\sin \theta }\text{ }=\frac{1}{2}$,此时有:
$$ I = {\rm{\pi }}{a^2}\cos \theta {{\rm{e}}^{{\rm{j}}2k{z_0}\cos \theta }} $$ (17) -
如图 3所示,圆柱面半径为a,高度为l,中心位于点 (0, 0, z0) 处。圆柱面上的源点和法向矢量可分别表示为$r' = (a\cos \phi ',a\sin \phi ',z')$,$\boldsymbol{\hat n}' = (\cos \phi ',$ $\sin \phi ',0)$,结合式 (10),散射积分可表示为:
$$ \begin{array}{c} I = \int_{{\rm{ }}S} {(\sin \theta \cos \phi \cos \phi ' + \sin \theta \sin \phi \sin \phi ')} \times \\ {{\rm{e}}^{{\rm{j}}2k[a(\sin \theta \cos \phi \cos \phi ' + \sin \theta \sin \phi \sin \phi ') + \cos \theta z']}}{\rm{d}}S' = \\ a\sin \theta \int_{{\rm{ }}{z_0} - \frac{l}{2}}^{{\rm{ }}{z_0} + \frac{l}{2}} {{{\rm{e}}^{{\rm{j}}2kz'\cos \theta }}} {\rm{d}}z'\int_{{\rm{ }} - \frac{{\rm{\pi }}}{2}}^{{\rm{ }}\frac{{\rm{\pi }}}{2}} {\cos (\phi ' - \phi ){{\rm{e}}^{{\rm{j}}2ka\sin \theta \cos (\phi ' - \phi )}}{\rm{d}}\phi '} \end{array} $$ (18) 其中,关于z'的积分可解析求得:
$$ {I_z} = {{\rm{e}}^{{\rm{j}}2k{z_0}\cos \theta }}\frac{{\sin (kl\cos \theta )}}{{k\cos \theta }} $$ (19) 关于$\phi '$的积分,因对于所有方向来波,柱面总是只有一半被照亮,所以$\phi $向的积分范围限定为半个周期。该积分可以由鞍点法进行近似计算。
鞍点法[13]首项近似公式为:
$$ \int_{ - \infty }^\infty {f(t){{\rm{e}}^{{\rm{j}}h(t)}}{\rm{d}}t \approx f({t_0}){{\rm{e}}^{{\rm{j}}h({t_0})}}\sqrt {\frac{{2{\rm{\pi }}}}{{ - {\rm{j}}h''({t_0})}}} } $$ (20) 其中t0为驻相点,即满足h'(t0)=0。
$$ \int_{ - \frac{{\rm{\pi }}}{2}}^{\frac{{\rm{\pi }}}{2}} {\cos (\phi ' - \phi ){{\rm{e}}^{{\rm{j2kasin}}\theta {\rm{cos(}}\phi ' - \phi )}}{\rm{d}}\phi ' \approx {{\rm{e}}^{{\rm{j}}2ka\sin \theta }}\sqrt {\frac{{\rm{\pi }}}{{{\rm{j}}ka\sin \theta }}} } $$ (21) 因此,式 (10) 可最终表达为:
$$ I = \frac{l}{k}{{\rm{e}}^{{\rm{j}}2k{z_0}\cos \theta }}\frac{{\sin (kl\cos \theta )}}{{kl\cos \theta }}{{\rm{e}}^{{\rm{j}}2ka\sin \theta }}\sqrt { - {\rm{j\pi }}ka\sin \theta } $$ (22) 当${\cos \theta }= 0 $时,$\frac{{\sin (kl\cos \theta )}}{{kl\cos \theta }} = 1$,此时:
$$ I = \frac{l}{k}\sqrt { - {\rm{j\pi }}ka\sin \theta } {{\rm{e}}^{{\rm{j}}2ka\sin \theta }} $$ (23) -
如图 4所示,圆台 (锥) 底面半径为a2,底面中心位于点 (0, 0, z2) 处,圆台顶面半径为a1,顶面中心位于点 (0, 0, z1) 处;圆锥顶点位于点 (0, 0, z0) 处,半锥角$\alpha = {\rm{arctg}}({a_1}/|{z_0} - {z_1}|)$,其中$\alpha \in (0,\pi /2)$。
圆台 (锥) 面上的源点r'用柱坐标变量可表示$r' = (\rho '\cos \varphi ',\rho '\sin \varphi ',z')$。由图 4中的几何关系可知,$\rho ' = ({z_0} - z'){\rm{tg}}\alpha $。r'处的法向矢量可表示为:$\boldsymbol{\hat n}' = (\cos \alpha \cos \varphi ',\cos \alpha \sin \varphi ',\sin \alpha )$。结合式 (10),散射积分可表示为:
$$ \begin{array}{c} I = \int_{{\rm{ }}{z_2}}^{{\rm{ }}{z_1}} {\int_{{\rm{ }}{\phi _1}}^{{\rm{ }}{\phi _2}} {({z_0} - z'){\rm{tg}}\alpha [\sin \theta \cos \alpha \cos (\phi - \phi ') + } } \\ \cos \theta \sin \alpha ]{{\rm{e}}^{{\rm{j}}2k({z_0} - z')\sin \theta tg\alpha \cos (\phi ' - \phi )}}{{\rm{e}}^{{\rm{j}}2kz'\cos \theta }}{\rm{d}}\varphi '{\rm{d}}z'/\cos \alpha = \\ \frac{{{\rm{tg}}\alpha }}{{{\rm{cos}}\alpha }}\int_{{\rm{ }}{z_2}}^{{\rm{ }}{z_1}} {({z_0} - z'){{\rm{e}}^{{\rm{j}}2kz'\cos \theta }}} \times \\ \int_{{\rm{ }}{\phi _1}}^{{\rm{ }}{\phi _2}} {(\sin \theta \cos \alpha \cos (\phi ' - \phi ) + \cos \theta \sin \alpha )} \times \\ {{\rm{e}}^{{\rm{j}}2k({z_0} - z')\sin \theta {\rm{tg}}\alpha \cos (\phi ' - \phi )}}{\rm{d}}\varphi '{\rm{d}}z' \end{array} $$ (24) 当θ < α时,圆台 (锥) 面全部为亮区,式 (24) 内层积分范围$[{\phi _1},{\phi _2}] = [ - \pi ,\pi ]$,内层关于$\phi '$的积分可由$\phi ' - \phi = 0$和$\phi ' - \phi = {\rm{\pi }}$两个驻相点的贡献近似计算;当α < θ < π-α时,圆台 (锥) 面仅一半为亮区,式 (24) 中内层积分范围$[{\phi _1},{\phi _2}] = [ - {\rm{\pi /2,\pi /2]}}$,关于$\phi '$的积分由$\phi ' - \phi = 0$单个驻相点的贡献近似计算;当$\theta = 0$时,圆台 (锥) 面均未被照亮,散射贡献被忽略。近似得到内层积分的表达式后,式 (24) 简化为关于z'的线积分,可用数值积分 (如辛普森积分)[14]快速计算。值得注意的是,当θ=0时,式 (24) 中的积分核变得非常简单,散射积分可以直接解析求得。
-
如图 5所示,z1, z2分别为椭球冠下底面和顶点的z轴位置,a1为下底面半径,a2为z向轴半径。
定义横向轴半径为a,则椭球面方程为:
$$ f(x',y',z') = \frac{{{{x'}^2}}}{{{a^2}}} + \frac{{{{y'}^2}}}{{{a^2}}} + \frac{{{{\left( {z' - {z_0}} \right)}^2}}}{{{a_2}}} = 1 $$ (25) 椭球面上的源点可以用球坐标变量表示为:
$$ \boldsymbol{r}' = (a\sin \theta '\cos \varphi ',a\sin \theta '\sin \varphi ',{a_2}\cos \theta ' + {z_0}) $$ (26) 由图中几何关系有:
$$ \left\{ \begin{array}{l} {a_1} = a\sin \alpha \\ {z_1} - ({z_2} - {a_2}) = {a_2}\cos \alpha \end{array} \right. $$ (27) 根据式 (27),可以由a1, a2, z1, z2求得椭球横向轴半径。因此,椭球冠形状可以由a1, a2, z1, z24个参量唯一确定。
因椭球面上任意一点处的外法向矢量与椭球函数f的梯度方向一致,r'处的法向矢量可表示为:
$$ \boldsymbol{\hat n}' = \frac{{\nabla f}}{{|\nabla f|}} $$ (28) 其中
$$ \nabla f{\rm{ = }}2\left( {\frac{1}{a}\sin \theta '\cos \varphi ',\frac{1}{a}\sin \theta '\sin \varphi ',\frac{1}{{a_2^{}}}\cos \theta '} \right) $$ (29) $$ \left| {\nabla f} \right| = 2\sqrt {\frac{1}{{{a^2}}}{{\sin }^2}\theta ' + \frac{1}{{a_2^2}}{{\cos }^2}\theta '} $$ (30) $$ \begin{array}{c} \boldsymbol{\hat n}' \cdot {{\hat k}_s} = \frac{1}{{\sqrt {\frac{1}{{{a^2}}}{{\sin }^2}\theta ' + \frac{1}{{a_2^2}}{{\cos }^2}\theta '} }} \times \\ \left[ {\frac{1}{a}\sin \theta \sin \theta '\cos (\varphi ' - \varphi ) + \frac{1}{{{a_2}}}\cos \theta \cos \theta '} \right] \end{array} $$ (31) $$ {\boldsymbol{\hat k}_s} \cdot \boldsymbol{r}' = a\sin \theta \sin \theta '\cos (\varphi ' - \varphi ) + \cos \theta ({a_2}\cos \theta ' + {z_0}) $$ (32) 将式 (31)、式 (32) 代入式 (10),得[15]:
$$ \begin{array}{l} I{\rm{ = }}\int_{{\rm{ }}0}^{{\rm{ }}\alpha } {\int_{{\rm{ }}{\phi _1}}^{{\rm{ }}{\phi _2}} {a\sin \theta '\frac{{\sqrt {a_2^2{{\sin }^2}\theta ' + {a^2}{{\cos }^2}\theta '} }}{{\sqrt {\frac{1}{{{a^2}}}{{\sin }^2}\theta ' + \frac{1}{{a_2^2}}{{\cos }^2}\theta '} }}} } \times \\ \left. {{\rm{ }}\left[ {\frac{1}{a}\sin \theta \sin \theta '\cos (\varphi ' - \varphi )} \right. + \frac{1}{{{a_2}}}\cos \theta \cos \theta '} \right] \times \\ {{\rm{e}}^{{\rm{j}}2k(a\sin \theta \sin \theta '\cos (\varphi ' - \varphi ) + {a_2}\cos \theta \cos \theta ' + {z_0}\cos \theta )}}{\rm{d}}\varphi '{\rm{d}}\theta ' \end{array} $$ (33) 与圆台面的分区讨论和计算类似,可由鞍点法对式 (33) 中关于$\phi '$的内层积分进行化简,仅剩下关于$\theta '$的线积分,然后通过辛普森积分实现快速计算。
-
目标的高频散射主要来自于光滑曲面的镜面反射,以及表面不连续处的边缘或尖端绕射。物理光学方法能较为精确地描述光滑曲面的镜面反射贡献[16]。用平板、圆柱、圆台和椭球冠等基本散射单元组合成弹体目标时,基本单元结合处的边缘会产生绕射现象,在某些角度可能成为目标RCS的重要组成部分,需要单独考虑其贡献。
用于目标绕射分析的方法很多,常用的方法包括几何绕射理论 (GTD)[17]、一致性几何绕射理论 (UTD)[18]、等效边缘电磁流 (EEC) 法[9]、增量长度绕射系数 (ILDC)[19]等,为了更好地与物理光学 (PO) 法进行结合,本文选择EEC法计算边缘绕射贡献。
EEC采用等效电磁流方法解决边缘绕射问题,其原理是:当任何有限电磁流分布的远区绕射场通过一个辐射积分来求和时,将得到一个有限的结果,若能找到这种适当的分布,则可避免几何绕射理论的发散问题。其基本方法是假设在环绕表面奇异性 (边缘回路) 的各点处存在线电流Ie和线磁流Im,并用远场辐射积分的形式表示绕射场:
$$ {\boldsymbol{E}^{{\rm{EEC}}}}(r,\boldsymbol{\hat k}) = \frac{{{\rm{j}}k{{\rm{e}}^{ - {\rm{j}}kr}}}}{{4{\rm{\pi }}r}}\int_l {[\eta {I_e}\boldsymbol{\hat k} \times (\boldsymbol{\hat k} \times \boldsymbol{\hat l}) + {I_m}(\boldsymbol{\hat k} \times \boldsymbol{\hat l})]{{\rm{e}}^{{\rm{j}}k\hat k \cdot r'}}{\rm{d}}r'} $$ (34) $$ \boldsymbol{H}_{}^{{\rm{EEC}}}(r,\boldsymbol{\hat k}) = - \frac{{{\rm{j}}k{{\rm{e}}^{ - {\rm{j}}kr}}}}{{4{\rm{\pi }}r}}\int_l {\left[ {{I_e}(\boldsymbol{\hat k} \times \boldsymbol{\hat l}) - \frac{1}{\eta }{I_m}\boldsymbol{\hat k} \times (\boldsymbol{\hat k} \times \boldsymbol{\hat l})} \right]{{\rm{e}}^{{\rm{j}}k\hat k \cdot r'}}{\rm{d}}r'} $$ (35) 式中,$\boldsymbol{\hat l}$为沿该边缘回路方向的单位矢量;$\boldsymbol{\hat k} = {\boldsymbol{\hat k}_s}$。
本文采用物理绕射理论-等效边缘电流法 (PTD-EEC)[20]计算目标边缘绕射作用,采用三角变换的方法去除虚假奇异点[21],再将求出的绕射贡献对物理光学法计算所得散射场进行修正,从而得到更为精确的RCS数据。
-
在雷达波束的照射下目标所表现出来的回波特性是雷达进行目标探测与识别的根本特征信息。常用的雷达目标特征信息包括目标的雷达散射截面 (RCS)、一维距离像、二维ISAR像,以及角闪烁和多普勒信息等[23]。
对目标进行一维和多维成像时,需要计算一定带宽或角度范围内各个采样点处目标的单站散射场信息。以二维ISAR成像为例,假定要求的成像口径范围为[−Xmax, Xmax] [−Ymax, Ymax],两个维度的分辨率分别为Δx和Δy,则对雷达信号的带宽和姿态角变化范围有如下要求[24]:
$$ B = \frac{c}{{2\Delta x}} $$ (36) $$ \Omega = \frac{{{\lambda _c}}}{{2\Delta y}} $$ (37) 同时,频域采样间隔和角度采样间隔需满足:
$$ \Delta f = \frac{c}{{2{X_{\max }}}} $$ (38) $$ \Delta \phi = \frac{{{\lambda _c}}}{{2{Y_{\max }}}} $$ (39) 求得各个采样点处的散射场数据后,对ES(f) 或ES(f, φ) 作傅里叶逆变换即可得到目标的一维距离像和二维ISAR像。
-
以某型导弹为例,如图 7所示,将该导弹分解为钝头锥、圆柱和圆台3个基本电磁散射计算单元。其中,钝头锥高1.1 m,底面半径为1 m,由一个圆台 (a1=1 m, a2=0.2 m, z1=9 m, z2=10 m) 和椭球冠 (a1=0.2 m, a2=2 m, z1=10 m, z2=10.1 m) 相切组合而成;圆柱半径为1 m,高6 m;圆台上下底面半径分别为1 m和2 m,高3 m。利用参数化组合建模及计算方法对其RCS、一维距离像和二维ISAR图像进行快速仿真计算。
RCS仿真频率f=1 GHz,雷达入射波的方位角为0°,俯仰角从0°开始以0.5°为步进增加至180°,通过本文方法快速计算得到的单站RCS与多层快速多极子 (MLFMA) 精确数值仿真方法[25]的计算结果对比如图 8所示。由图可见,两种计算方法得到的目标随入射角 (俯仰角) 变化所表现出的散射特性的趋势是一致的。该算例中,即使基于24个核心的高性能工作站进行计算,MLFMA算法耗时2.6 h才能得到图中所示的一条单站RCS曲线,占用计算机内存33.5 GB。而本文开发的组合参数建模方法在普通的笔记本电脑上单核运算仅耗时2 s,占用内存250 MB,可满足实时计算的应用需求。
采用相对均方根误差对两种计算结果的差异进行衡量。对于该导弹模型,参数化组合 (本文方法) 计算结果与MLFMA结果的相对均方根误差仅为0.5%。可见对于导弹目标,参数化组合建模计算得到的RCS具有非常高的精度:
$$ {\rm{Relative RMS Error}} = \sqrt {\frac{{\sum\limits_{i = 1}^N {{{\left| {{\sigma _i} - \sigma _i^{{\rm{ref}}}} \right|}^2}} }}{{\sum\limits_{i = 1}^N {{{\left| {\sigma _i^{{\rm{ref}}}} \right|}^2}} }}} $$ (40) RCS曲线随俯仰角变化在51.5°、71.5°和90°共3个角度出现峰值,这和导弹模型的钝头圆锥侧面法线、截头圆锥侧面法线,以及圆柱体法线方向保持一致,在这些角度刚好有较强的镜面反射。在0°和180°出现的峰值刚好和导弹顶部和底部的镜面反射保持一致。
为了得到导弹的一维距离像信息,设置轴向观测范围为15 m,距离分辨率为0.25 m,雷达波中心频率为1 GHz,根据式 (36) 和式 (38) 计算得到观测频率范围为[0.7, 1.3]GHz,频率采样间隔Δf=1 MHz,共需61个采样频点。图 9给出了分别采用本文方法和MLFMA算法计算得到的导弹模型在轴线上 ($\theta = 0^\circ $) 的一维距离像对比。
雷达回波信号峰值主要由不连续结构处的感应电流引起。因此雷达回波信号标识的距离像信息与导弹物理尺寸能够非常准确地对应。由图可见,虽然各个散射峰的幅度略有差异,但两种分析方法都能准确地标识导弹模型的距离像信息。同时观测到在MLFMA的距离像中多了一个与爬行波对应的散射峰。这是因为MLFMA通过数值方法精确求解麦克斯韦方程,可模拟物体各部件间的电磁耦合和精确描述物体表面的感应电流分布;而本文方法的算法基础是物理光学近似,无法模拟爬行波的贡献。
进一步对导弹模型进行二维ISAR像分析。设置径向 (距离维) 观测范围为15 m,距离分辨率为0.25 m,横向 (方位维) 观测范围为10 m,空间分辨率为0.2 m,雷达信号中心频率为1 GHz,中心观测角度为$\theta = 0^\circ $。由式 (36)~式 (39) 计算得到频率范围为[0.7, 1.3]GHz,频率采样间隔Δf=10 MHz,角度范围为[-21.49°, 21.49°],角度采样间隔$\theta = 0.86 ^\circ $,对应61个频率采样点和51个角度采样点,需3 111次单站散射计算。若基于MLFMA精确数值算法,在24个核心的高性能工作站上计算需要耗时约30 h,才能完成该ISAR图像提取。而基于本文提出的参数化组合建模和快速数值仿真方法,在普通的PC机上运行时间不到30 s就可完成仿真分析。由图 10可见,ISAR图像中的散射热点分布与导弹的物理结构几乎完全一致。
通过以上的应用案例分析不难看出,参数化组合建模及计算方法基本能够实现导弹目标电磁散射特性的实时计算,仿真精度满足工程应用需求。这为总体方案设计时,对调整外形参数后目标散射特性快速估算提供了有效的实现途径,并为试验方案验证提供了高效的解决手段。
Compositional Modeling and Electromagnetic Scattering Characteristics Extracting for Missile Targets
-
摘要: 该文将常规导弹目标分解为平板、圆柱、圆台和球冠等基本几何体散射单元,通过参数化组合建模、散射场解析计算和矢量场叠加的方法,快速仿真得到导弹目标的电磁散射数据,以支持其RCS、一维距离像和二维ISAR图像等电磁散射特征提取。计算结果表明,该方法基本能够满足目标电磁散射实时、高精度计算需求,开发的仿真软件为导弹目标基于方案验证需求的电磁散射特性快速计算与评估提供了有力支持。Abstract: In this paper, conventional missile targets are divided into basic scattering cells, such as disks, cylinders, cones and ellipsoids. By carrying on the parameterized compositional modeling, analytical calculation of integrals, and superposition the field contribution of each unit, we can obtain the scattering fields of missile-like objects in time. Consequently, the electromagnetic scattering characteristics such as RCS, the range profile and ISAR images are available. The numerical results demonstrate the proposed method can meet the real-time requirement as well as the high precision. Numerical software is developed and it gives a best support to the RCS calculation and design evaluation of missile targets in engineering application.
-
[1] 聂在平.目标与环境电磁散射特性建模-理论、方法与实现[M].北京:国防工业出版社, 2009. NIE Zai-ping. Modeling of electromagnetic scattering characteristics of target and environments-theory, method and realization[M]. Beijing:National Defense Industry Press, 2009. [2] 庄钊文, 袁乃昌, 莫锦军, 等.军用目标雷达散射截面预估与测量[M].北京:科学出版社, 2007. ZHUANG Zhao-wen, YUAN Nai-chang, MO Jin-jun, et al. Prediction and measurement of radar cross section for military targets[M]. Beijing:Science Press, 2007. [3] 傅君眉, 冯恩信.高等电磁场理论[M].西安:西安交通大学出版社, 2000. FU Jun-mei, FENG En-xin. The advanced electromagnetic field theory[M]. Xi'an:Xi'an Jiaotong University Press, 2000. [4] 方宙奇, 孟敏.电磁场数值方法[M].成都:电子科技大学出版社, 2012. FANG Zhou-qi, MENG Min. Numerical methods for electromagnetic fields[M]. Chengdu:University of Electronic Science and Technology of China Press, 2012. [5] 盛新庆.计算电磁学要论[M].北京:中国科学技术大学出版社, 2008. SHENG Xin-qing. Computational electromagnetics[M]. Beijing:University of Science and Technology of China Press, 2008. [6] 张京国, 梁晓庚.基于物理光学法和面元法的目标近场RCS计算[J].探测与控制学报, 2008, 30(6):42-45. http://www.cnki.com.cn/Article/CJFDTOTAL-XDYX200806011.htm ZHANG Jing-guo, LIANG Xiao-geng. Calculating near-field RCS of targets based on physical-optics method and panel method[J]. Journal of Detection & Control, 2008, 30(6):42-45. http://www.cnki.com.cn/Article/CJFDTOTAL-XDYX200806011.htm [7] 严俊. 参数曲面上的高频渐进方法研究与应用[D]. 成都: 电子科技大学, 2015. YAN Jun. Research and application of parametric surface based high frequency asymptotic methods[D]. Chengdu:University of Electronic Science and Technology of China, 2015. [8] KNOTT E F. A progressing of high-frequency RCS prediction techniques[J]. Proc IEEE, 1985, 73(2):252-264. doi: 10.1109/PROC.1985.13137 [9] MICHAELI A. Equivalent edge currents for arbitrary aspects of observation[J]. IEEE Transactions on Antennas and Propagation, 1984, 32(3):252-258. doi: 10.1109/TAP.1984.1143303 [10] STRATTON J A. Electromagnetic theory[M]. New York:McGraw-Hill, 1941. [11] KNOTT E F. Radar cross section measurements[M]. New York:Springer Science & Business Media, 2012. [12] 马振华.现代应用数学手册[M].北京:清华大学出版社, 2005. MA Zhen-hua. Handbook of applied mathematics applied mathematics handbook[M]. Beijing:Tsinghua University Press, 2005. [13] CHEW W C. 非均匀介质中的场与波[M]. 聂在平, 柳清伙, 译. 北京: 电子工业出版社, 1992. CHEW W C. Waves and fields in inhomogeneous media[M]. Translated by NIE Zai-ping, LIU Qing-huo. Beijing:Publishing House of Electronics Industry, 1992. [14] 《数学手册》编写组.数学手册[M].北京:高等教育出版社, 1979. "Mathematics Handbook" Compilation Group. Mathematics handbook[M]. Beijing:Higher Education Press, 1979. [15] 徐宏枢.三轴椭球面的面积计算[J].渝州大学学报, 1998, 15(1):57-60. http://www.cnki.com.cn/Article/CJFDTOTAL-YZZK801.011.htm XU Hong-shu. The calculation of the area of the 3-axes ellipsoid[J]. Journal of Yuzhou University.1998, 15(1):57-60. http://www.cnki.com.cn/Article/CJFDTOTAL-YZZK801.011.htm [16] CONDE O M, PEREZ J, CATEDRA M F. Stationary phase method application for the analysis of radiation of complex 3-D conducting structures[J]. IEEE Transactions on Antennas and Propagation, 2001, 49(5):724-731. doi: 10.1109/8.929626 [17] KELLER J. Geometrical-theory of diffraction[M]. New York:McGraw-Hill, 1958. [18] KOUYOUMJIAN R G, PATHAK P H. A uniform geometrical theory of diffraction for an edge in a perfectly conducting surface[J]. Proceedings of the IEEE, 1974, 62(11):1448-1461. doi: 10.1109/PROC.1974.9651 [19] SHORE R A, YAGHJIAN A D. Application of incremental length diffraction coefficients to calculate the pattern effects of the rim and surface cracks of a reflector antenna[J]. IEEE Transactions on Antennas & Propagation, 1993, 41(1):1-11. [20] UFIMTSEV P Y. Method of edge waves in the physical theory of diffraction[M]. Moscow:Soviet Radio, 1962. [21] 崔索民, 吴振森. PTD和PTDEEC算式的奇异点处理[J].电波科学学报, 1997, 4(12):369-374. http://www.cnki.com.cn/Article/CJFDTOTAL-DBKX199704003.htm CUI Suo-min, WU Zhen-sen. Treatment of singularities in PTD and PTDEEC[J]. Chinese Journal of Radio Science, 1997, 4(12):369-374. http://www.cnki.com.cn/Article/CJFDTOTAL-DBKX199704003.htm [22] Mathworks. 面向科学计算的MATLAB和Simulink[EB/OL]. [2016-09-10]. http://cn.mathworks.com. Mathworks. Model-Based Design with MATLABand Simulink[EB/OL].[2016-09-10]. http://cn.mathworks.com. [23] 黄培康, 殷红成, 许小剑.雷达目标特性[M].北京:电子工业出版社, 2005. HUANG Pei-kang, YIN Hong-cheng, XU Xiao-jian. Characteristics of radar targets[M]. Beijing:Publishing House of Electronics Industry, 2005. [24] OZDEMIR C. Inverse synthetic aperture radar imaging with MATLAB algorithms[M]. Manhattan, America:John Wiley & Sons, Inc, 2012. [25] 胡俊. 复杂目标矢量电磁散射的高效方法-快速多极子方法及其应用[D]. 成都: 电子科技大学, 2000. HU Jun. The efficient method for vector electromagnetic scattering from complex object-fast multipole method and its application[D]. Chengdu:University of Electronic Science and Technology of China, 2000.