徐文丰, 李颖晖, 裴彬彬, 宋亚南, 孙鹏飞
(1.空军工程大学航空工程学院,西安,710038; 2.95841部队,甘肃酒泉,735000)
雷达散射截面(radar cross section,RCS)是表征目标在雷达波照射下所产生的回波强度的物理量,是飞机隐身技术中最基本和关键的指标之一[1]。因此,对RCS的计算和测量具有十分重要的意义。按照测算目标的状态分类,针对RCS的研究可分为动态[2-5]和静态[6-7]2类。关于动态RCS的研究由于考虑了目标机动、紊流、天气等实际因素对静态RCS进行修正补充,有更高的工程价值。按照测算的方法分类,获取RCS的方法可分为包含微波暗室法[8],标准外场测量[9]等方法的实验测量法和计算仿真法[10-12]。实验测量法是获取目标RCS的主要方式,但其通常需要耗费大量人力、物力和时间,实验不能大量开展;计算仿真法则作为实验测量法的补充和先验知识,具有成本低、效率高等优势。
目前,现有大多数研究成果针对位置和姿态发生变化的飞机动态RCS仿真计算方法是:首先建立目标RCS关于雷达俯仰角和方位角的二维数据表格或统计模型,然后通过坐标变换,得到雷达在飞机本体系下的坐标,进而得到雷达在飞机本体系下的俯仰角和方位角,最后根据飞机本体系中的俯仰角和方位角计算当前姿态和位置下的飞机动态RCS值。这种方法简单容易实施,基本能够适应大多数情形下的飞行。近年来,国内外学者们对动态RCS的研究也基本在该方法的框架内开展,例如,文献[13]中建立了无人机关于俯仰角和方位角的二维表格,通过查表法获得不同姿态下的飞机RCS值,并基于稀疏A*算法提出一种隐身突防航迹规划方法。文献[14]对飞机动态RCS的统计模型展开研究,针对经典模型描述复杂航迹情况下的隐身飞机动态雷达散射截面统计分布特性精度不足的问题,提出一种精度更高、拟合效果更好的改进混合对数正态分布模型。文献[15]针对现有的典型统计分布模型无法满足双基地雷达隐身飞机动态电磁散射分布特性精度要求的问题,提出一种改进的混合对数正态分布模型。但这些方法无法对不同姿态位置下的RCS进行精确表征。文献[16]建立了以对数均值和对数标准差为参数的动态RCS统计模型,并基于动态RCS统计模型研究了飞机动态RCS峰值对雷达探测概率的影响,但其只考虑了雷达方位角,未考虑飞机的俯仰和滚转。文献[17~18]中给出了一种使用飞机滚转角、飞机相对于雷达的方位角和俯仰角来共同表征飞机动态RCS的方法,相对于传统方法具有更高的准确性,能够准确计算平飞条件下的飞机动态RCS,但其假定飞机轨迹在平面内运动,不适用于飞机存在俯仰的情形。
对此,本文提出了一种基于雷达坐标系和飞机本体系之间的欧拉旋转角的飞机动态RCS表征和计算方法。
式中:变换矩阵Tbr为:
Tbr=TbrxTbryTbrz
(2)
式中:
本文中定义雷达站心坐标系:X轴指向北,Y轴指向东,Z轴指向地。定义飞机本体坐标系:X轴指向机头方向,Y轴沿机体轴线指向右,Z轴垂直于XOY平面向下,二者均为右手系。
定义俯仰角θr与方位角φr如图1所示,即θr为本体系原点-雷达连线与本体系Z轴的夹角,φr为本体系原点-雷达连线在XOY平面的投影与本体系X轴的夹角。
图1 飞机本体系中的雷达俯仰角与方位角
计算俯仰角和方位角:
最后根据预先建立的关于θr、φr的二维全空域动态RCS数据库插值计算当前姿态和位置下的RCS值,也可基于FEKO[19]等电磁仿真软件对当前θr和φr条件下的RCS进行实时交互计算。
式(3)中arctan(·)函数的值域为[-π/2,π/2],而方位角需要在[-π,π]范围内取值,这显然是不能满足要求的,举例说明,考虑以下2种情形:
情形1:
[xp1,yp1,zp1]T=[100,100,-100]T
[φ1,θ1,ψ1]T=[0,0,0]T
情形2:
[xp2,yp2,zp2]T=[-100,-100,-100]T
[φ2,θ2,ψ2]T=[0,0,0]T
2种情形下飞机和雷达的相对位置关系见图2。
(a)情形1
(b)情形2
由图2可知,2种情形下的RCS值理应不同。然而按式(3)计算得到二者的方位角和俯仰角为φr1=φr2=45°,θr1=θr2=45°,故2种情形下的RCS值相同,这显然是不正确的。这是由于式(3)只能在一、二、五、六卦限内正确计算雷达与飞机的方位角导致的。
针对本节中指出的传统方法只能在有限角域内正确反映雷达方位角的问题,修正雷达方位角的计算公式如下:
这样定义的φr能够在[-π,π]内取值,其正方向与飞机欧拉偏航角ψ的正方向保持一致,使得基于雷达俯仰角和方位角表征RCS的方法能够适用于全卦限范围。
本节中将以图3所示的F-5飞机模型为例,对修正前后的动态RCS计算方法进行基于FEKO电磁学仿真的对比验证。
图3 FEKO中的F-5飞机模型。
选取飞行轨迹从东经50°0′0″,北纬60°0′0″,先后经过航路点1:东经50°0′0″,北纬60°12′0″,航路点2:东经50°18′0″,北纬60°12′0″。飞机的飞行高度保持为5 000 m,雷达的位置位于东经50°6′0″,北纬60°6′0″,海拔高度为0 m。飞行轨迹及雷达的相对位置见图4,飞行过程中飞机的速度保持在200 m/s。飞行过程中飞机的欧拉姿态角见图5,雷达在飞机本体系中的坐标变化情况见图6。
仿真中使用多层快速多极子算法(MLFMA)[20],设置雷达极化方式为水平极化,水平接收。照射频率为1 GHz。由式(1)、式(2)、式(4)计算不同姿态角与位置下的θr、φr,在FEKO中相应设置雷达相对于飞机的方位仿真得到的飞行动态RCS结果见图7。
假定修正后传统方法的FEKO电磁计算的结果为动态RCS的真实值,计算修正前的方法将会带来的相对误差为:
图4 雷达站心坐标系下的飞行轨迹(航路点1-2)
图5 飞机欧拉姿态角(航路点1-2)
图6 雷达在飞机本体系中的坐标(航路点1-2)
图7 修正前后的动态RCS对比
图8 修正前传统方法的RCS相对误差
由图6~8可见,在0~53 s和124~153 s时,雷达在飞机本体系中的坐标位于第一、五卦限,2种方法的RCS值相同,而飞机在其他时间段内处于第三、四、七、八卦限,未经修正的传统RCS计算方法由于无法清晰的分辨当前雷达相对于飞机的方位角而造成较大的误差,相对误差最大可达14 000%。
需要说明的是,这里动态RCS仿真计算可能是不够准确的,其准确性程度与在FEKO软件中选择的计算方法,以及网格划分的精度紧密相关。但这种不完全准确的方法对说明本文对动态RCS方法进行的修正是足够的,因为这种基于飞机姿态信息建立动态RCS数据库或统计模型的方法也可用于其他RCS仿真算法或飞行试验中。
尽管在第2节中对基于雷达方位角和俯仰角计算动态RCS的算法进行了修正,但其仍存在无法精确表征不同姿态下飞机RCS的问题,在本节中将对这一问题进行说明。
式(4)中定义的φr、θr,描述了在飞机本体坐标系中雷达的视线角,即雷达与本体坐标系原点的连线与本体系三轴的位置关系,然而若飞机和雷达的位置确定,φr、θr并不能唯一确定飞机的欧拉姿态角φ、θ、ψ,若将飞机以雷达视线轴为轴进行旋转,φr、θr保持不变,即一组雷达视线角φr、θr有无穷多组欧拉姿态角φ、θ、ψ与之相对应。虽然在不同的欧拉姿态角下,飞机本体系中的X、Y、Z轴与雷达视线角的相对方位保持不变,但其相对于电磁波中电场矢量和磁场矢量的相对方位发生改变。这无疑改变了飞机的RCS值。
显然,将飞机以雷达视线角为轴进行旋转的情形与保持飞机姿态不变、将雷达装置以视线角为轴进行旋转的情形相等价,这相当于改变了雷达的极化方式,根据电磁波散射理论,不同极化方式对同一目标有不同的RCS值[21]。因此,φr、θr并不能与动态RCS形成一一对应关系。
假定飞机在雷达站心系中的位置[xp,yp,zp]T已知,本体系中雷达相对于飞机本体系坐标原点的方位角和俯仰角φr、θr已知,求解飞机欧拉姿态角φ、θ、ψ。
由式(1)~(2)可得:
(sinφsinθsinφ+cosφcosφ)yp+sinφcosθzp]
(cosφsinθsinφ-sinφcosφ)yp+cosφcosθzp]
假定飞机在雷达站心系中的位置[xp,yp,zp]T=[-1 000,0,0]T,令[φ3,θ3,ψ3]T=[0°,0°,0°]T,[φ4,θ4,ψ4]T=[60°,0°,0°]T
2种情形下的[φr,θr]T均为[0°,90°]T。采用水平极化方式,2种情形下垂直于雷达视线角的截面示意图见图9(a~b)。
在FEKO中同样设置雷达频率为1 GHz,使用多层快速多极子算法进行动态RCS仿真,得到水平接收回波信号情形下的2种姿态的RCS值分别为-7.036 dBsm,-7.395 dBsm。若将[φ4,θ4,ψ4]T姿态下的雷达极化角度旋转60°,见图9(c)。
仿真计算得该情形下的RCS为-7.036 dBsm,与[φ3,θ3,ψ3]T在水平极化下对应的RCS相同。显然,飞机姿态绕雷达视线角旋转等价于电磁波极化方向的旋转,不同旋转角下的飞机有不同RCS值。
(a)[φ3,θ3,ψ3]T
(c)极化角度旋转60°时的[φ4,θ4,ψ4]T
图10 雷达站心系中的飞机俯仰角与方位角
其数学表达式如下:
雷达照射系、雷达站心系、飞机本体系三者的旋转关系如图11所示。
图11 坐标系之间的转换关系
图11中的箭头指向代表通过坐标系旋转进行的坐标系变换,其中红色,绿色,蓝色的箭头分别表示将坐标系绕其自身的X,Y,Z轴进行右手旋转,旋转的次序由箭头的次序确定,旋转的角度由箭头上标定的角度确定。雷达站心系Fr中的坐标转换至飞机本体系Fb的转换矩阵Tbr已在式(2)中给出,雷达照射系Fs中的坐标转换至飞机本体系Fb的转换矩阵为:
Tbs=TbsxTbsyTbsz
(7)
式中:
雷达照射系中Fs的坐标转换至雷达站心系Fr的转换矩阵为:
Trs=TrsxTrsy
(8)
式中:
显然有:
Tbs=TbrTrs
(9)
综上,给出单组航路飞行数据的动态RCS在线仿真计算流程见图12所示。
图12 单组航路飞行数据的动态RCS在线仿真计算流程
本节中将通过仿真实例对比第2节中的传统基于飞机本体系雷达视线角计算动态RCS的方法和第4节中提出的基于欧拉角旋转角[φx,θx,ψx]T计算动态RCS的方法。同样使用图4中给出的F-5飞机模型,极化方式选择水平极化、水平接收,雷达频率设置为1 GHz,解算方法选择精度更高的多层快速多极子算法[22]。
设置飞行轨迹为从东经50°0′0″,北纬60°0′0″,高度5 km出发,先后经过航路点1:东经50°0′0″,北纬60°6′0″,高度7 km,航路点2:东经50°0′0″,北纬60°12′0″,高度4 km,航路点3:东经50°18′0″,北纬60°12′0″,高度4 km。雷达位于东经50°6′0″,北纬60°6′0″,海拔高度0 m。速度保持为200 m/s。仿真过程中飞机与雷达的位置及其欧拉姿态角随时间的变化见图13~14。
图13 雷达站心坐标系下的飞行轨迹(航路点3-5)
图14 飞机欧拉姿态角(航路点3-5)
(10)
2种方法的动态RCS对比,以及按式(10)计算所得的传统动态RCS计算方法产生的误差分别在图15~16中给出。
图15 传统方法与本文中方法的对比
图16 传统方法的RCS相对误差
可见,当飞机处于平飞阶段时,2种方法得到的结果差距较小,而在0~20 s,60~80 s,120~160 s时间段时,飞机姿态角较大,传统的方法会产生一定的误差。在本例中,相对误差的峰值达到1 300% 。这说明本文中所提出的动态RCS计算的新方法能够在飞机处于机动阶段时大幅度提高动态RCS的准确性。
1)针对传统的基于飞机本体系内雷达方位角和俯仰角计算动态RCS方法只能在有限角域内有效的问题,本文对方位角进行了重新定义,使得传统动态RCS方法具有更大的适用范围。
2)指出了修正后的传统动态RCS方法仍存在的问题,即飞机的动态RCS值无法与飞机本体系中定义的方位俯仰角[φr,θr]T建立一一映射关系,不能精确表征不同姿态与方位下的飞机动态RCS值。
3) 提出一种新的动态RCS计算方法,定义了雷达照射坐标系,并根据雷达照射坐标系与飞机本体系之间的欧拉旋转角[φx,θx,ψx]T表征飞机的动态RCS值。给出了根据飞机的欧拉姿态角[φ,θ,ψ]T与雷达站心系中的方位俯仰角[ψr,θr]T计算[φx,θx,ψx]T的方法,可按此方法计算不同姿态和位置信息下的[φx,θx,ψx]T,进而通过插值或在线计算的方法获取当前姿态和位置下的飞机RCS。