1 引言

空间目标的光学散射特性是外部因素和内部因素共同影响的结果。其中外部因素指太阳、目标和地面观测站在空间中的相对位置;内部因素主要指目标自身的特征,包括目标的几何结构、目标形状和姿态以及表面材质等。对于光学监视设备而言,在对目标进行光学探测和跟踪的过程中,若目标的轨道信息已知,就可以确定探测过程中太阳、目标和观测站三者之间的相对位置关系,此时空间目标的光学散射特性是空间目标自身特征的体现,而描述空间目标的光学散射特性是开展目标识别和状态判定的前提^[1]。总结国内外研究现状,大部分均是对空间目标的外形结构和表面材质进行了简化,其中二面元模型是最具代表性和采纳广泛的模型^[2-5]。该模型的其中一个面元代表卫星上的太阳能帆板,另一个面元代表卫星本体,要求帆板面元近似于平整平面并具有镜面和朗伯反射特性,本体为具有特定反射率的朗伯体。文献[ 6]对卫星褶皱表面的光学散射特性进行了建模研究,指出卫星本体褶皱表面的存在可以将卫星本体的表面近似为朗伯材质,当卫星本体无镜反现象时,印证了二元模型模拟结果的正确性。

实际观测中,地球同步轨道(GEO)卫星与地球低轨(LEO)卫星均出现了特殊的闪光现象。文献[ 7-8]在实测数据中观测到了GEO卫星的帆板闪光现象,分析了卫星帆板的运动角速度,但未得到闪光的规律。文献[ 9]基于卫星帆板双向反射分布函数及卫星运动特性,量化研究了当GEO卫星帆板追踪太阳时引起季节性闪光现象的特点和规律。与GEO卫星不同的是,LEO卫星具有轨道周期短、相对地面运动速度快的特点,当卫星经过测站上空时,地面望远镜每次可观测的时间较短(约为几分钟至十几分钟不等),太阳和观测站的空间几何关系较GEO卫星更为复杂,导致其闪光现象的出现存在随机性,闪光时间较短(约为几十秒)且与季节无关。卫星帆板追踪太阳和天光地影时受到观测条件的限制,导致地面观测站不能捕捉到LEO卫星帆板的闪光现象,说明LEO卫星的闪光现象由卫星本体上的镜反部件所致。

二面元模型只能反映卫星本体的朗伯散射特性,无法描述卫星本体包含镜反现象的光学散射特性。针对上述问题,提出一种描述空间目标的光学散射特性的多元模型,该模型在二元模型的基础上增加镜反部件,用于描述卫星本体含有镜反部件的卫星光学散射特性。

2 卫星镜反材质的散射特性

为了实现对卫星镜反部件的光学特性进行反演和描述,首先对空间目标常用的镜面反射材质的光学散射特性进行分析。采用本课题组改进的Phong模型^[10],对卫星表面材质的光学散射特性进行双向反射分布函数(BRDF)^[11]的计算,表达式为

fr=ρdcosθi+ρscosαβcosaθiexp[-b(1-cosβ)1/π],(1)

式中:ρ_d和ρ_s分别为材质的漫反射系数和镜面反射系数;α为镜向指数;β为反射矢量偏离镜面反射矢量的角度;a>0,用以调节菲涅耳现象强度的系数;b>0,用以调节镜面反射分量增降速度的系数;θ_i为入射天顶角,引入cos θ_i来修正漫反射项,用以调节镜面反射的强度。前期,本课题组在400~900 nm的波长范围内测量了大量卫星表面常用材质的BRDF值,并基于这些数据反演了卫星表面材料的BRDF模型参数,所用材质对应的模型参数如表1所示。

将裸露在卫星本体表面的阳极化镜片作为分析对象,以30°入射角为例,在镜反平面内考察BRDF值随着出射角度的变化规律,结果如图1所示。

图 1. BRDF与出射角度的变化曲线

Fig. 1. Variation curve of BRDF and exit angle

下载图片 查看所有图片

从图1可以看到,阳极化镜片的镜反特性显著。研究镜反部件的光学特性,除了BRDF幅度外,还需要对不同入射角度下的镜反角度宽度进行统计分析,结果如图2所示。

图 2. 镜面反射角度宽度与入射角度的变化曲线

Fig. 2. Variation curve of specular reflection angle width and incident angle

下载图片 查看所有图片

从图2可以看到,当入射角度较小时(偏离镜面的角度),镜面反射角度宽度较小,随着入射角度的增大,镜面反射角度宽度逐渐变大,镜面反射角度宽度的均值约为5.278°。需要指出的是,实验以朗伯材质BRDF值为0.31作为参考,BRDF值大于0.31判定为镜反,5.278°是BRDF从镜反峰值向两侧分别下降到0.31时对应的角度宽度。

以边长为5 cm的阳极化镜片为例,在太阳光以30°入射的条件下,距离测站1500 km轨道处不同角度下的星等值,结果如图3所示。

图 3. 30°入射角度下的星等曲线

Fig. 3. Magnitude curve under incident angle of 30°

下载图片 查看所有图片

从图3可以看到,将边长为5 cm的阳极化镜片放在距离测站1500 km的轨道上,在镜面反射角度(30°)外,目标亮度约为12星等,在镜反角度内,下降到6星等左右,此时目标迅速变亮,地面观测站可观测到约6星等的亮度。从镜面反射外到镜面反射内,在测站观测到的目标亮度提高了近600倍。虽然其镜面反射角度宽度很窄,但能量很集中,且大部分亮度被面元镜反。低轨卫星的运行速度比测站快,导致观测的几何关系变化快,当测站处于镜反部件的反射方向附近时,测站将会观测到目标短暂的镜反现象。根据光的反射原理,可确定镜面反射部件的指向。

3 多元法描述卫星的光学散射特性

多元法模拟卫星的光学散射特性分三步展开:首先,构建包含卫星本体、帆板和镜反部件的三维卫星模型;其次确定参数,一是通过拟合光度的曲线来确定卫星本体的反射率,二是通过空间观测的几何关系来确定镜反部件的指向和面积;最后将多元模型导入光度计算程序中,可对含有镜反现象的卫星光学散射特性进行描述,理论分析过程如下。

3.1 构建卫星的多元模型

多元模型包含三个要素:第一个要素为方形的卫星本体,卫星本体的表面视为朗伯材质,反射率待定,在其表面可安装具有朗伯材质特性的天线,表面包覆白漆;第二个要素为卫星帆板,在其表面包覆空间目标帆板常用的材质,当进行卫星光度数据计算时卫星帆板追踪太阳;第三个要素为镜面反射部件,将其安装在目标本体的下方,镜反部件的表面包覆镜面反射材质的阳极化镜片,镜反部件的指向待定,初始时设定很小的面积。建立的卫星3D模型如图4所示。

图 4. 3D卫星模型示意图

Fig. 4. Schematic of 3D satellite model

下载图片 查看所有图片

其中body-x、body-y和body-z分别表示卫星的本体轴,实验将卫星姿态角度设为零,此时卫星质心轨道的坐标系与卫星本体的坐标系重合,body-x指向目标速度方向,body-z对地定向,body-y垂直于body-x和body-z,三者满足右手定律。采用常用的3DS MAX建模软件来构建3D模型,对采用3ds格式导出的文件进行目标光度曲线的计算。计算模型采用基于OpenGL拾取技术来实现快速获取可见面元信息^[12],以计算目标的OCS(Optical Cross Section)值。其中,星等与OCS^[13]的关系为

m=-26.74-2.5·lg(xOCS/R2),(2>)

式中:m为星等;R为探测器与卫星的距离。3D卫星模型创建后,对卫星本体的反射率、镜反部件的面积和指向确定展开研究。

3.2 本体反射率的计算

在不考虑镜反现象的条件下,逐步调整卫星本体的反射率,使得计算的卫星光度曲线和实测光度曲线的平均误差最小。在整体的变化趋势上拟合实测光度曲线,调整卫星本体反射率的方法如下。

采用两个相同尺寸的卫星本体,其反射率分别为ρ₁和ρ₂,由(2)式可知,两者的星等之差为

Δm=2.5·lg(ρ1/ρ2)。(3)

已知ρ₂的卫星星等值,与之相差Δm的ρ₁应为

ρ1=ρ210Δm/2.5。(4)

卫星的光度曲线仿真计算中,对于Δm需采用两条曲线的平均误差来逐步调整卫星本体的反射率。在此过程中,可结合地基望远镜来探测卫星的STK(Satellite Tool Kit)场景,探测卫星的过程中,太阳照射和地基观测卫星本体的不同表面,调整不同表面的反射率,使得仿真计算的光度曲线与实测光度曲线的整体趋势和幅值保持一致。

3.3 镜反部件的指向确定

在卫星质心轨道的坐标系下,当实测数据出现镜反现象时,太阳和测站的坐标分别为(x_S,y_S,z_S)和(x_F,y_F,z_F),其单位矢量为N_S和N_F。根据光的反射原理,镜反部件的法向矢量为N_F+N_S,将其单位化为N_M=(x_M,y_M,z_M)。N_M的方向即是当前时刻卫星本体中镜反部件的指向,采用方位角和俯仰角来描述矢量方向,β_M为法向偏离body-z轴正向的角度,α_M为法向在本体坐标系xOy平面内与body-x轴正向的夹角,如图5所示。

图 5. 镜反部件的反射示意图

Fig. 5. Reflection diagram of mirror reflection component

下载图片 查看所有图片

因方位角与俯仰角不存在耦合关系,所以可用来单独描述矢量绕单一坐标轴的转动情况,但无法描述一个矢量转向另一矢量的过程。在3DS MAX和STK中均采用欧拉角法^[14]来描述目标的转动过程,为了方便在3DS MAX软件中建模,采用欧拉角法来描述初始矢量N₀=(x₀,y₀,z₀)转向N_M的过程,转动后的欧拉角为目标的法向。

欧拉角法的原理是转动坐标系,按照XYZ的顺序旋转坐标系。则N₀与N_M的关系可表示为

NM=RN0,(5)

式中:R=R(θ_z)R(θ_y)R(θ_x),其中R(θ_x)、R(θ_y)、R(θ_z)分别表示坐标轴X、Y、Z依次旋转θ_x、θ_y、θ_z的旋转矩阵。可得N_M为

xM=y0(cosθxsinθz+cosθzsinθxsinθy)+z0(sinθxsinθz-cosθxcosθzsinθy)+x0cosθycosθzyM=y0(cosθxcosθz-sinθxsinθysinθz)+z0(cosθzsinθx+cosθxsinθysinθz)-x0cosθysinθzzM=x0sinθy+z0cosθxcosθy-y0cosθysinθx。(6)

将本体z轴的正向作为N₀的方向,即N₀=(0,0,1),则N_M为

xM=sinθxsinθz-cosθxcosθzsinθyyM=cosθzsinθx+cosθxsinθysinθzzM=cosθxcosθy。(7)

求解(7)式的欧拉角θ_x、θ_y和θ_z,将本体的坐标系作为参考坐标系,依次旋转坐标轴即可将N₀转向N_M。当不转动参考坐标系时,某一矢量方向在参考坐标系中转向另一矢量方向,这是一个绕任意轴旋转的问题。矢量绕任意轴旋转与旋转坐标轴的过程相反,N₀和N_M均处于参考坐标系中,N₀在参考坐标系中分别绕本体坐标轴z、y、x旋转-θ_z、-θ_y、-θ_x即可转至N_M。由于N₀=(0,0,1),此时-θ_z=0,则N_M为

xM=-cosθxsinθyyM=sinθxzM=cosθxcosθy。(8)

求解(8)式非线性方程组的部分解,可得

θy=-arctan(xM/zM)θx=arcsinyM,-90°≤θy,θx≤90°。(9)

由(9)式可知,在该部分解下,上述矢量旋转过程只能实现N₀在卫星本体坐标系下I、II、III、IV卦限内的转动。实际的地基望远镜观测空间目标的光照条件:空间目标被太阳直射且地面测站处于地球全影或半影区,使得N₀和N_M只在卫星本体坐标系下I、II、III、IV卦限内转动,而(9)式的解满足实际地基的观测需求。

为了方便在3DS MAX软件中构建卫星模型,根据STK卫星本体坐标系向3DS MAX坐标系的过渡矩阵,将卫星本体坐标系下的坐标在3DS MAX坐标系中进行表示。建模过程中,模型的旋转过程是在3DS MAX软件中完成的,依据(5)~(9)式推导矢量之间的旋转过程,在3DS MAX坐标系中,N₀依次绕3DS MAX坐标系Y轴、X轴旋转θ_y、θ_x,即可转至矢量N_M,其中

θy=arctan(xM/zM)θx=arcsinyM,-90°≤θy,θx≤90°。(10)

3.4 镜反部件的面积确定

在镜反的时间范围内,认为卫星星等的变化是由镜反部件引起的。在镜反的时间范围外,镜反部件不会引起卫星星等的变化,并且在短时间内由卫星本体引起的卫星星等变化较小,这与图3镜反角度较窄的结论一致。对于相同材质的两个镜反平整平面,面积分别为A₁和A₂,将其法向转至镜反时刻的法向。由OCS计算公式及其与星等转化关系可知,两者在镜反时刻的星等之差为

Δm1=2.5·lg(A1/A2)。(11)

对已知面元A₂的星等值,则面元A₁的面积为

A1=A210Δm1/2.5。(12)

与本体反射率不同的是,这里Δm₁指镜反部件在镜反时刻的星等差值,即镜反处两条星等曲线的峰值之差。需要指出的是,当镜反部件的面积确定时,可用镜反的强度和范围来描述镜反部件,镜反强度由镜反部件的面积来确定,镜反范围由镜反材质的BRDF镜反系数和漫反系数两个主要参数来确定。镜反材质的属性虽不完全对应目标真实镜反材质的属性,但可对目标镜反这一现象的描述。

确定本体反射率、镜反部件指向和面积后,将模型导入OCS计算程序即可拟合实测数据,进而实现对包含镜反现象的空间目标光学散射特性的描述。

4 实测验证

采用丽江国家天文台(25.6°N,101.1°E)获取的某低轨卫星光度数据对所提模型进行验证。根据实测数据镜反点出现的时刻并结合观测几何关系,对目标的镜反部件指向进行提取,根据3.2~3.4节分析给出的卫星本体反射率与面积确定方法在卫星本体上添加镜反部件。

表2为目标出现的镜反时刻、镜反法向单位方向矢量、镜反部件的尺寸和在3DS MAX软件中由N₀转向N_M所需的角度。

在不加入镜反部件的条件下,整体的拟合结果和实测数据如图6所示。

图 6. 无镜反部件的仿真计算结果和实测数据对比

Fig. 6. Comparison of simulation calculation results andmeasured data of mirrorless components

下载图片 查看所有图片

无镜反部件的仿真计算结果和实测数据的绝对误差如图7所示,其中平均误差为0.4404星等。从图7可以看到,目标光度曲线出现两次明显的镜反现象,说明目标包含两块镜反部件。在3DS MAX软件中建立平整平面,其法向矢量为N₀,按照表2的角度旋转平面,使得平面法向与在镜反时刻目标的法向保持一致。添加镜反部件的仿真计算结果和实测数据对比如图8所示。

图 7. 无镜反部件的绝对误差

Fig. 7. Absolute error of mirrorless components

下载图片 查看所有图片

图 8. 添加镜反部件的仿真计算结果和实测数据对比

Fig. 8. Comparison of simulation calculation results and measured data of adding mirror reflection components

下载图片 查看所有图片

图 9. 添加镜反部件的绝对误差

Fig. 9. Absolute error of adding mirror reflection components

下载图片 查看所有图片

从图8可以看到,在实测数据中出现两个镜反时刻点,仿真数据在相应的时刻点均出现镜反,验证镜反部件指向确定的正确性,表明提出的多元模型描述卫星的光学散射特性与二元模型相比更为全面,能够更具体地反映卫星个体独有的光学散射特性。含有镜反部件的仿真计算结果和实测数据的绝对误差如图9所示,其中平均误差为0.3294星等。相比不含镜反部件,含有镜反部件的精度提高了0.1110星等。

产生误差的原因:一方面是对卫星本体进行朗伯近似,为了提高精度,在后续的反演工作中,需考虑卫星表面包裹的真实材料和表面褶皱的情况;另一方面,实验采用的阳极化镜片的BRDF参数与卫星实际的镜反部件之间有差异,由此导致镜反部分并非与实测数据完全重合,尤其表现在镜反范围内,目标的镜反部件真实BRDF参数需进一步通过光谱解混等手段来反演;当望远镜在观测刚开始和快结束时,观测条件不稳定会导致观测数据本身的均方差比较大,从而产生较大的拟合误差。

5 结论

卫星闪光是卫星个体特有的现象,能够反映卫星本身独有的特征,其是研究卫星运动和工作状态的重要依据之一。为此,提出卫星光学散射特性多元模型,在二元模型的基础上增加对卫星本体镜反部件要素的光学特性描述。首先,将目标本体视为朗伯材质,通过拟合实测光度数据来确定其反射率;其次,结合卫星的运动特性解析镜反部件的指向,计算镜反部件的面积;最后,根据镜反部件的指向和面积的计算结果,在卫星本体表面包覆镜反部件,实现目标出现镜反时卫星光学散射特性的模拟,以更为全面地反映卫星细节的光学散射特性。工程应用中,初始模型可采用二元模型,根据实际观测数据来调整卫星散射特性描述模型的元数,可实现对任意目标光学散射特性的模拟。