程勋龙,康国华,邱钰桓,武俊峰,吴佳奇,尹一蓁
(南京航空航天大学 航天学院,江苏 南京 211101)
激光雷达(Light Detection and Ranging,Lidar)因其波束窄,抗干扰能力强[1]等优点被应用于空间非合作目标的识别和预警[2-3]。通过对区域内目标跟踪、搜索以及测量,激光雷达可实时获取航天器目标的三维点云信息[4],尤其是在无光照条件下,激光雷达依然可以工作[5-6],因此成为当前空间安全防御的主要威胁。针对激光雷达建立激光雷达散射特性模型,将对后续航天器的激光雷达隐身性能评估起到重要作用。
对于目标激光雷达散射特性的研究,哈尔滨工业大学李程程[7]通过电磁散射解析理论分析不同入射波长激光、目标表面粗糙度对激光散射特性的影响,在此基础上建立三维目标的激光雷达截面(Laser Radar Cross Section,LRCS)计算方法,虽然此方法分析了不同粗糙度材质的特性,但是其光束模型采用了能量分布均匀的点光源与实际情况存在一定差异。
宫彦军等[8]采用差异积分法克服统一积分法中存在的舍掉误差,分析了观测方向的方位角和天顶角对LRCS 的影响,以及三锥体的几何参数对后向LRCS的影响,但是其在计算朗伯三锥体的LRCS 时缺少对于多材质复杂表面情况的分析。西安电子科技大学徐付昌[9]基于GRECO 算法,实现了不同光斑下的复杂目标单、双站的LRCS 计算,同时还考虑了局部照射下光斑无法完全覆盖的情况。张昊等[10]提出了一种图形电磁计算方法(Graphical Electromagnetic Computing,GRECO)以及双向反射分布函数(Bidirectional Reflectance Distribution Function,BRDF)模型的复杂目标LRCS 计算软件的一体化开发方法,实现了对复杂目标的LRCS 快速求解。李良超[11]提出了计算覆盖多种材料复杂表面的激光雷达散射特性可视化方法,使用颜色标签建立含材料信息的目标模型,并且通过OPENGL 消隐解决遮挡问题。孙鹏举等[12]在一般计算方法的基础上,重点针对镜面反射目标,提出了一种简便的LRCS 工程估算方法,并应用于激光角反射器目标和“猫眼”目标,获得相应的LRCS估算公式。JIN 等[13]提出了一种快速可视化的LRCS求解办法,其利用颜色标签包含材料信息,并利用五参数经验模型快速求解目标LRCS。西安电子科技大学杨昭[14]通过蒙特卡洛光线追迹算法计算目标在激光波段的LRCS 散射特性,分析了二次散射对计算结果的影响。上述文献虽然实现了复杂目标的LRCS 求解,但是都忽视了激光高斯能量分布特性,会导致实验结果与理论值存在一定的误差。
针对上述情况,曲卫东[15]分析了光斑模式为高斯能量分布与均匀分布之间存在的误差,并且提出将消光布遮蔽法与光斑分析法相结合,以提高LRCS 测量精度。西安电子科技大学罗龙[16]研究了均匀照射情况与非均匀高斯光束照射的差异,并提出补偿参数R 的概念,及其与目标形状尺寸的关系;但是其研究目标为简单几何体,并未对复杂表面目标进行分析。王明军等[17]分析了激光光束的高斯能量分布特性,并研究了在波长1.06 μm 激光束入射下,粗糙球体目标散射场量的二阶统计特征。上述工作都考虑了激光的高斯能量分布特性对于激光雷达散射特性产生的影响,但都是针对简单表面目标进行分析。
综上所述,以往对于激光雷达的目标散射特性分析通常采用简单可见光模型求解,未充分考虑激光特性;并且对于非连续、多材质表面的影响也未考虑充分,因此建模往往不够准确,散射特性分析和实测相差较大。本文基于激光特性和复杂航天器表面精确模型,提出一套更接近真实情况的LRCS 仿真算法,以期为航天器激光隐身的设计分析提供一个更为准确的路径。
在以往激光雷达散射特性建模分析中,将光束能量分布近似为均匀分布,如图1(a)所示,把其能量特性类似于可见光;然而激光雷达光束的强度在截面上的分布呈现高斯分布[18],其能量分布如图1(b)所示;因此等效成均匀能量分布最后会使计算结果不够精确。同时对比两者可以看出高斯能量分布的分布趋势主要为中间较高,并且逐渐向四周递减,而均匀分布则在同一截面上各点能量大小相同。
图1 均匀能量分布、高斯能量分布及Lidar 坐标系Fig.1 Schematic diagram of uniform and Guassian energy distributions and the Lidar coordinate system
针对这种情况,采用基于高斯能量分布的激光雷达光束模型用于目标特性分析。其坐标系如图1(b)所示,z轴正方向为激光光束的传播方向,ω(z)为与传播轴线相交于z点的高斯光束等相位面上的光斑半径。选取基模高斯球面波,在空间中(x,y,z)点的电矢量表达式[19]为
式中:A0为常数因子;ω0为基模高斯光束的束腰半径,m;f为高斯光束的共焦参数,m;R(z)为与传播轴线相交于z点的高斯光束等相位面的曲率半径,m;λ为激光波长,m;L为激光谐振腔腔长,m。
激光雷达散射截面其定义式[20]为
式中:σ为激光雷达散射截面,m2;R为目标与探测器之间的距离,m;Ii为入射到目标上的能流密度,W/ m2;Ir为反射到探测器上的能流密度,W/m2。
由式(7)可知,对于航天器LRCS 的求解,关键在于Ir。基于有限元的思想,通常航天器表面采用三角面元的划分,因此Ir的求解就转化为单个三角面元LRCS 在表面上的积分。
以激光雷达为原点建立坐标系,如图1(b)所示,首先设一小面元dS,入射功率为dp,入射角为θ,根据朗伯余弦定理,求解出反射能流密度的表达式为
式中:ρB为后向反射率;R为目标与探测器之间的距离,m。
则对于面积为Striangle的三角面元其反射回波能流密度为
单位三角面元的入射功率dp为
由式(8)和式(10)可得,三角面元的反射回波能流密度为
式中:x1、x2为三角面元在x轴上的积分下限、上限;y1、y2为关于x的积分下限表达式和积分上限表达式。
考虑高斯能量分布的情况,其在空间中的电场分布由式(1)可得,将其带入式(11)得到考虑高斯能量分布情况下的三角面元反射回波能流密度为
已知入射能流密度表达式[21]为
式中:Pi为入射功率大小,W。
此时将式(12)~式(14)带入式(5)可以求出三角面元的激光雷达散射截面表达式
式中:σtriangle为单个三角面元的LRCS 值,m2。
复杂表面航天器表面结构复杂、材料属性不单一,如将其简单地视为朗伯体,得到的LRCS 结果往往不符合实际情况[22]。针对材料属性多样,且物体后向散射率ρB无法直接求解的问题,通过使用BRDF 统计模型获取不同材料在不同探测条件下的后向散射率,以此改进复杂表面航天器的激光雷达散射特性模型。单站激光雷达BRDF 模型如图2 所示,采用基于菲涅尔反射现象的改进Phong 模型[23],其公式如式(16)所示:
图2 单站激光雷达BRDF 计算模型Fig.2 BRDF calculation model of single-station Lidar
式中:ρd为材质的漫反射系数;ρs为材质的镜面反射系数;α为镜像指数;cosθi为引入的修正漫反射项,用以调节镜面的反射强度;β为观测方向与镜面反射方向的夹角,β=min{π/2,β};a>0,以调节菲涅尔现象的强度;b>0 用以调节镜面反射分量的增降速度。
常用的航天器表面材质上述模型参数见表1。
表1 3 种常见航天器表面材质的Phong 模型参数Tab.1 Phong model parameters of three common spacecraft surface materials
通过航天器面元划分,可以将复杂表面分解成形状和材质单一的面元,从而解决了外形复杂、材质多样化的问题。但划分以后,沿着激光入射方向不可避免地存在面元遮挡情况。
面元遮挡分为两种方式:自遮挡和互遮挡。自遮挡即表面未被激光直接照射,如图3(a)所示,表示激光雷达照射方向的矢量i与表示面元法向的向量n1之间角度超过90°,即i·n1≤0。由此可得式(17)。
图3 面元自遮挡以及互遮挡情况Fig.3 Diagram of self-occlusion and mutual occlusion of face elements
式中:Sselfjug为面元自遮挡判断状态;当Sselfjug=1 时表示不存在自遮挡;Sselfjug=0 时表示存在自遮挡。
另外一种即为互遮挡,需要判断是否存在2 个面元之间的遮挡情况。如图3(b)所示。从图中可以看出面元A与面元B在激光雷达方向上的投影存在遮挡。通过选取常见的Z-buffer 消隐算法[24]实现面元的遮挡判断:首先以激光雷达照射方向i为z轴建立坐标系,并将A、B两面元投影至xoy平面内,同时储存两面元的深度信息即zA和zB值。其遮挡判断的具体表达如式(18)所示。
式中:Sbuffer为面元剔除状态;So为两面元的重叠面积,当两面元的重叠面积So>0 时,表示两面元存在重叠;zA和zB分别为A、B两面元的深度值。当Sbuffer=A时表示面元A被剔除,Sbuffer=B表示面元B被剔除,Sbuffer=0 表示没有面元被剔除。
大部分情况我们利用Z-buffer 消隐算法能够对面元遮挡情况进行判断,但是也会遇到误遮挡的情况,例如2 个面元只有极小部分产生遮挡或因计算机计算精度产生误判断,这些情况会导致最终LRCS 计算结果的精度降低。因此对传统Z-buffer算法进行适当的改进,添加面积容差的判断,提高LRCS 的计算精度。
以下主要总结了3 种误遮挡的情况,分别为两面元异面、两面元共线以及两面元共点,如图4所示。
图4 面元误遮挡情况Fig.4 Schematic diagram of false occlusion of surface elements
1)两面元异面表示2 个面元不存在交点,在这种情况下,会出现两面元仅由极小部分遮挡而被判定为互遮挡的情况,最终导致面元被舍弃,进而损失一部分精度,如图4(a)所示。
2)两面元共线表示两面元存在一条公共边,由于计算精度不够的原因,导致两面元存在极小的重叠面积,而两面元实际并没有遮挡,此种情况也会导致LRCS 计算结果的不准确,如图4(b)所示。
3)两面元共点表示两面元仅存在一个公共点,与上述情况(2)类似,因为计算精度导致误差,误认为两面元存在遮挡,如图4(c)所示。
针对上述3 种情况产生的误遮挡,提出了一种加入面元容差判断的算法,其具体算法流程如图5所示。图5 流程表明,首先需要计算两面元公共顶点数量N,判断两面元是否存在公共点:
图5 误遮挡判断流程Fig.5 Flow chart of error occlusion judgment
1)当两面元不存在公共点即N=0 时,则计算两面元的重叠面积So,两面元的投影面积S1和S2,如果So<min(S1,S2)/10,则认为两面元不存在遮挡,反之认为两面元存在遮挡。
2)当两面元存在公共点即N≠0 时,则计算其中一个面元在另一个面元内部的顶点的数量I,如果I=0 表明两面元不存在包含关系,接下来判断如果存在So<10-10,则认为两面元不存在遮挡,反之认为两面元存在遮挡。如果I≠0,则说明两面元存在包含关系,认为两面元存在遮挡。
经过对两面元位置关系判断、面元容差的设置后,消除了计算误差导致的两面元遮挡判断错误的情况,以及两面元仅存在极少部分重叠就判断为遮挡的情况。
同时在激光雷达存在部分照射的情况,如图6(a)所示。如果对部分照射的面元直接全部舍弃,将会降低计算结果的精度,所以本文主要采用对部分照射面元进行切割的方法,保留面元被照射到的部分。具体操作方法为:当三角面元ABC出现部分照射的面元时,首先计算三角面元与光斑边缘的交点J、K,将面ABC沿交线JK划分成两部分,并判断三角面元在光斑内部顶点的数量,当存在2 个顶点在光斑内部时,保留梯形面元JBCK,并将梯形面元JBCK沿对角线划分成2 个三角面元JBC和KJC,如图6(b)所示。当存在一个顶点时,则保留面元JAK,如图6(c)所示。此方法通过对部分照射面元进行切割提升了算法的精度。
图6 面元部分照射情况及面元切割Fig.6 Diagram of part irradiation and cutting of surface elements
综合改进后的Z-buffer 消隐算法,基于高斯光束的复杂表面航天器的激光雷达散射特性的求解步骤如图7 所示。
图7 考虑高斯能量分布的目标LRCS 计算流程Fig.7 Flow chart of target LRCS calculation considering the Gaussian energy distribution
具体步骤如下:
1)3D 模型面元划分:首先利用3Dmax 建模软件建立复杂目标3D 模型,并且通过优先面元法对3D 模型进行三角面元划分,确认材质属性,随后将目标3D 模型以obj 格式导出。
2)激光雷达参数设置:确认激光雷达的波长、激光谐振腔腔长,照射距离等参数。
3)坐标系转化:将obj 格式模型中的面元坐标以及顶点坐标转化到激光雷达坐标系下。
4)有效面元判断:利用矢量法以及Z-Buffer 消隐算法判断每次激光光束照射到的有效面元,同时对部分照射的面元进行面元切割。
5)目标整体LRCS 计算:结合BRDF 模型,计算当前光束照射下的所有有效面元的LRCS,并对整个光斑内的三角面元的LRCS 进行累加,获得整体的LRCS。
通过构造有解析解的标准体对前面提出的方法进行分析和验证。在考虑高斯能量分布的情况下,当激光雷达照射在朗伯体圆形平面上时,其LRCS 值存在解析解。因此选择朗伯体圆形薄片验证本算法的准确性。
当目标为朗伯体表面时,其半球反射率为ρ2π,BRDF为ρ2π/π,入射角为ϕ,探测距离为z,光斑半径为ω(z),则可以得到朗伯体圆形平面的LRCS 表达式为[25]
为了验证面元遮挡判断算法的准确性,构造了一个含内外层的圆形薄片,其3D 模型如图8(a)所示;圆形薄片三维重建图如图8(b)所示:红色部分为外层,蓝色部分为内层;外层半径R,内层半径R’,外层薄片厚度H,内层薄片高度H’。
图8 圆形薄片几何示意及三维重建Fig.8 Geometric diagram and 3D reconstruction of circular wafer
调节光斑大小ω(z),使其刚好等于当前圆薄片的半径R,仿真参数分别为:波长λ=1.064×10-6m,激光谐振腔腔长L=1 m,半球反射率ρ2π=0.5,入射角ϕ=0,探测距离z随光斑半径大小改变。假设外层高度H=0.1 m,光斑半径取值范围R=[0.1,1],每次取值间隔0.1 m,从而形成10 个外层圆形薄片;在这外层圆形薄片下嵌套同一大小的内层薄片,其半径假设R’=0.08 m,高度H’=0.005 m。添加Z-buffer 消隐算法前后的LRCS 计算值如图9所示。
图9 不同半径圆形薄片Z-buffer 互遮挡消隐判断前后LRCS 对比曲线及误差曲线Fig.9 Comparison and error curves of LRCS before and after Z-buffer mutual occlusion blanking judgment on circular slices with different radii
对比图9(a)可以看到,针对不同半径圆形薄片,将采取Z-buffer 互遮挡消隐判断后的LRCS 值与仅有外层圆形薄片的LRCS 值相比较,可以看出两者曲线吻合度较好,由图9(b)可以看出,未采取Z-buffer 互遮挡消隐判断的双层圆形薄片与仅有外层圆形薄片之间的LRCS 计算结果存在较大误差,且在当双层圆形薄片内外两层的面积相差越小时,误差越大,证明了消隐算法的有效性和必要性。
同时对比消隐后的LRCS 仿真值与圆形薄片理论值,以及均匀光束理论值以验证消隐算法的准确性,对比结果如图10 所示。
图10 不同半径圆形薄片的LRCS 仿真算法与解析法的对比Fig.10 Comparison of results for circular slices with different radii obtained by the LRCS simulation algorithm and the analytical method
图10(a)中黑色曲线为解析解,即可认为是真值,红色曲线为考虑高斯能量分布下的LRCS 仿真结果,蓝色曲线为考虑均匀能量分布下的LRCS 仿真结果,可以看出考虑高斯能量分布下LRCS 计算方法的结果与理论值基本吻合。由图10(b)可以看出考虑高斯能量分布下的LRCS 计算结果误差均值为0.56%,而分析造成这种误差的原因主要在于利用三角面元拟合圆面时,其边缘仍是不断趋近于圆的多边形,会造成一部分面积的缺失,所以在计算LRCS 的时候,也会产生一定的计算误差。同时对比考虑均匀能量分布下的LRCS 计算结果则与理论值相差角度,误差均值为15.13%,可以得出在考虑高斯光束能量分布下的LRCS 计算精度提升了14.58%。
选取常见的航天器3D 模型进行复杂表面航天器的LRCS 计算,航天器主要由包覆金聚酰亚胺的卫星本体,以及由单结砷化镓电池板组成太阳帆和铝制金属卫星天线3 部分组成,如图11(a)所示。仿真Lidar 采用单站模式进行LRCS 探测仿真。基于运动相对性,把探测器绕卫星进行LRCS 探测,等效为Lidar 不动,而卫星在转动,一共有2 种模式:模式1:Lidar 在Z方向,航天器绕X轴旋转一周,如图11(b)所示;模式2:Lidar 在X方向,航天器绕Z轴旋转一周,如图11(c)所示。
图11 航天器本体及Lidar 探测Fig.11 Schematic diagram of spacecraft body and Lidar detection
假设Lidar 光斑半径ω(z)=0.7 m,能够覆盖整个航天器;基于改进Z-buffer 消隐算法后的LRCS计算结果如图12 所示。
图12 Lidar 扫描下卫星不同姿态的LRCS 周向Fig.12 LRCS circumferential map of spacecraft with different attitudes under Lidar scanning
由图12 中的仿真曲线分析可知:
1)当Lidar 入射角为0°、90°、180°、360°时LRCS值较大,且呈现对称分布。因为在这4 个角度方向上,卫星外型特征呈现对称分布,并且是激光反射最大的方向。
2)将图12(a)和图12(b)合成就可以得到模型的LRCS 全向图。从数据分析可以看出,当Lidar 直射通信天线、卫星底部、太阳翼时LRCS 最大;在顺着太阳翼方向探测时,LRCS 最小,因为这些方向下反射面积最大或者最小,整个LRCS 趋势与实际情况相吻合。
本文提出了一种改进的复杂航天器激光雷达散射特性建模方法,以提高模拟的精度。
1)首先用更准确的高斯能量分布模型替代均匀能量分布模型建立激光光束模型,减少了由光源特性不准带来的误差。
2)而针对多材质、非连续表面带来散射特性难模拟的问题,采取有限元思想,通过对单一材质和形状的面元辐射特性的积分予以解决。其中激光光束照射下,面元遮挡导致的散射特性误差的问题,则通过改进Z-Buffer 消隐算法进一步消除。
3)最后仿真分析表明,本文提出的这套完整算法精度提高了14.58%;而对典型复杂表面航天器在不同激光雷达入射角度下的LRCS 分析,得到了与实际情况相一致的结论。
4)不过由于本文工作主要是和理论解析解进行对比分析,还缺乏实际验证,因此下一步将采用真实的激光雷达和卫星模型来进一步验证;但本文工作对于后续航天器的激光雷达隐身设计提供了有价值的参考模型。