侯晴宇 巩晋南 樊志鹏 王一惠

(哈尔滨工业大学空间光学工程研究中心,哈尔滨 150001)

1 引 言

在天基光学空间目标态势感知领域,目前的研究侧重于远距离点目标检测、目标轨迹的预测与确认.而在远距离成像无目标结构、纹理信息时,利用时序点信号对空间目标的在轨工作状态及基本物理属性参量进行反演具有重要的意义,可有效支撑空间目标状态判别以及在轨维护的相关决策.

天基光学空间目标态势感知的前提为目标的特性认知.目前,在轨目标光学成像特性的研究分为实测以及仿真两类方法.由于在轨测量的成本较高,得到的测量数据相对有限,因此在该领域多采用地面材料特性测量与成像仿真结合的方式实现在轨目标成像特性研究.近些年来,对于目标材料表面可见光波段双向反射分布函数(BRDF)测量以及光度信号、图像仿真的研究较为深入,给出了多种空间目标材料BRDF的实验室测量方法[1,2]以及建模方法[3],并建立了基于BRDF的空间目标图像仿真方法[4−7].

空间目标的可见光波段光度信号由目标反射太阳光得到,它是目标外形结构、尺寸、姿态指向、表面材料等物理属性参量的函数[8,9].由于耦合量多,基于光度信号的参量估计的难度较大,一般情况下需要给定先验知识.在目标表面材料反射特性模型已知的条件下,国内外开展了一些基于地基光学探测系统对于空间目标姿态、外形的估计研究.Calef等[10]在假设目标姿态和指向已知的前提下,利用时序光度和热辐射数据反演目标的三维外形.Hinks等[11]分析了姿态变化与光度信号变化之间的关系,以及利用时序光度信号推导姿态变化的可行性.Wetterer等[12,13]从滤波角度给出了姿态估计的流程.文献[14]跟踪了国内外基于光度数据反演目标特征信息的最新动态,分析了目前基于地基平台进行光度数据姿态反演的主要方法及特点.文献[15]研究了国外几种基于点目标测量信息的空间目标识别方法的基本原理和实现途径.文献[16]基于时序光度信号匹配的思想开展了卫星形状反演的研究.文献[17]基于地面实验测量获得了目标光谱特性,反向提取目标的材料、大小和状态等特征,具有一定的参考价值.

在上述文献中,针对空间目标的相关状态及特性反演研究多以目标材料的光学特性为已知输入.但是,在实际应用中,非合作目标的本体材料组成未知,参数反演难度较大.因此,本文以含太阳翼的三轴稳定卫星为研究对象,基于双面模型及BRDF的多级融合表征模型,建立了融合帆板在轨动态特性的光学特性宏观表征模型,以在轨测量的时序光度信号为输入,建立优化方法估计模型参数,完成了对宏观表征模型的重构;时序光度信号仿真以及参数反演仿真验证了本文所提方法的有效性.

2 可见光波段时序光度建模

2.1 总体方案

为了体现空间目标在轨过程中帆板与本体运动特性上的差异,本文在文献[4—7]的基础上进行了相应完善,融入了帆板的对日指向计算模型.在轨空间目标的可见光波段时序光度特性建模流程如图1所示.首先,根据空间目标的几何结构模型进行面元划分,确定本体面元的法向量及位置等相关属性;然后,根据空间目标的三轴稳定姿态,确定该时刻的空间目标照明与观测矢量.考虑帆板对日指向,确定帆板面元的法向量及位置等相关属性;再根据照明与观测矢量以及面元的相对位置信息,确定既能够被光照又能够被观测的有效面元.结合有效面元的BRDF模型及太阳的辐射特性计算面元的反射特性,结合观测星的位置确定光瞳处接收到的总的光度特性.

图1 在轨空间目标时序光度特性建模流程Fig.1. Modeling process of sequence photometric characteristics of on-orbit space object.

2.2 时序光度信号模型

空间目标成像特性建模需要定义下列坐标系:卫星本体坐标系Rb(OXY Z)、帆板坐标系Rp(OXpYpZp)、相机坐标系R′(OX′Y′Z′)及探测器坐标系R′′(OMN),如图2所示.

各个坐标系之间的矢量变换满足:

式中RY(β)为绕Y(Yp)轴转动β角时的旋转矩阵;ξ,ψ,ζ分别为相机坐标系坐标方向在本体坐标系下的单位矢量;T0是相机坐标系原点在本体坐标系下的坐标;M,N分别为探测器的行列数;a为探测器像元尺寸;f为光学系统焦距.

图2 坐标系之间的关系示意图Fig.2.Relationship diagram among coordinate systems.

卫星在轨运行过程中,反射能量绝大部分为太阳辐射,因此本文研究的卫星反射特性主要是指卫星反射直接入射太阳光的特性.设本体坐标系下空间目标第j个面中第k个面元在入瞳方向的辐射强度为反射光谱强度时波长和时间的函数,表示为

式中Esun(λ)为太阳辐照度;Aj,k为第j个面中第k个面元面积;面元矢量nj,k,oj,k为本体坐标系下第j个面中第k个面元-观测星入瞳中心的单位矢量;s为本体坐标系下空间目标-太阳的单位矢量;fj(λ,nj,k,oj,k,s)表示第k个面元为第j种材料时的BRDF,〈·〉为点积运算.其中,oj,k,s,fj(λ,nj,k,oj,k,s)都为时变量.

探测系统入瞳接收到的面元辐照度为

光学系统入瞳接收到的目标辐通量为

式中Aaperture为光学系统的入瞳直径.

光学系统入瞳接收到的目标辐通量为

考虑到远距离点目标成像过程,面元-观测星光学系统入瞳中心的单位矢量基本一致,记为o.(7)式简化为

可以看出,空间目标光度信号与面元在本体坐标系下的法向量nj,k以及照明矢量s和观测矢量o有关.

2.3 帆板的对日指向计算模型

对于三轴稳定卫星,本体面元的相关属性如位置、法向量在本体坐标系下保持不变.但对于太阳帆板面元,由于保持对日指向运动,帆板与本体之间存在实时相对转动,因此需要单独计算太阳帆板面元在本体坐标系下的法向量.

由于太阳距地球距离远,入射太阳光可近似为平行光.对于三轴稳定卫星,一般具有较大的太阳能帆板,从卫星姿态控制的角度考虑,帆板通常为单自由度旋转,即帆板绕卫星本体OY轴旋转,则帆板坐标系下照明矢量sp表示为

式中,s为本体坐标系中的太阳光入射方向矢量,s=[sx,sy,sz].

为确保太阳翼获得最大能量,需保证帆板法线与太阳光入射矢量之间的夹角最小.为保证该条件,如图2,在帆板坐标系内,sp的Xp分量为零,即

因此,本体坐标系下帆板的法向量表示为

3 光度模型的在轨重构方法

3.1 空间目标的等效双面模型

对于正常工作状态下的对地定向三轴稳定卫星,对其进行同轨道面或近轨道面观测时,目标本体和观测系统的相对姿态基本不发生变化,使得本体被观测到的表面时序不变,假设该表面为平面(如图3所示),并且本体和帆板不发生相互遮挡,可将观测模型简化为双面模型,表示为

式中

其中Φp(t,λ),Φb(t,λ)分别为帆板和本体的光通量;Ap为帆板的面积;Ab为本体中被观测表面的面积.

图3 对三轴稳定卫星的在轨观测过程Fig.3.On orbit observation process of three axis stabilized satellite.

在轨观测过程中,光度信号Φ(t,λ),Esun(λ),Aaperture由光学成像系统获得,为已知量;T0以及ζ,o,nb与观测几何相关,可基于空间定轨、定位方法获得,可认为是已知量;np可以根据帆板对日指向的约束计算得到.Apfp(λ,np,o,s)和Abfb(λ,nb,o,s),即本体和帆板的面积-BRDF乘积,为未知量.

3.2 材料表面BRDF的多级融合表征模型

对于太阳帆板,其表面材料为电池片,电池片的材料组分相差不大,其BRDF可以通过地面测量得到[3].对于空间目标本体,由于卫星功能差异,可能存在多种包覆不同表面材料的构件,因此其BRDF未知,这里采用多级融合表征模型,该模型表示为

式中bm为与β′m(nb,o,s)对应的反射率,b={bm},bm需要满足0≤bm≤1,bm≤1.表示为

m=0时,表示漫反射;m逐渐增大,表示镜面反射的增强.在这里,提出应用Cook-Torrance BRDF作为参数化模型,即

式中α为n和b之间的夹角,b=(s+o)/2.

因为该模型考虑了面元微观分布特性,并且对漫反射以及不同角分布的镜面反射的表征能力较强,其形状由µm调节,µm表征了镜反峰的宽度,可将其设定为0.04≤µm≤0.5,µm按照一定的间隔进行取值.此时,(14)式表示为

式中Q(t,λ)=〈np·s〉〈np·o〉·Apfpp(λ,np,o,s)+,则未知量为面积与反射率的乘积fpp及{Abbm}m=0,1,···,M.

3.3 光学特性模型的在轨重构方法

将以上Q(t,λ)式转化为矩阵形式有

式中Q为T×1维向量,

βp为T×1向量,其中向量元素=p();βb为T×M矩阵,其中矩阵元素[为,

将(19)式简化为

(20)式的误差向量可以用矩阵向量形式表示:r(B′)=Q−βB′,由于B′中元素的物理意义为反射率面积的乘积,所以对B′中的元素做如下约束:b′≥ 0.

特征矩阵β为帆板和本体的BRDF信息矩阵,当T＞M时,即特征矩阵β行数大于列数时,只需适当选取B′使误差r(B′)在2范数意义下最小,即‖r(B′)‖=‖Q−βB′‖最小,即可对B′进行求解,即完成了空间目标复杂材料表面的物性参数反演以及宏观表征模型的重构.之所以称之为宏观表征模型,是因为该模型利用BRDF的基函数重构得到,是微观BRDF模型的宏观光度体现.

4 仿真验证

4.1 时序光度信号的仿真

4.1.1 仿真输入

卫星本体尺寸为1.2 m×1.2 m×1.2 m,帆板尺寸为6.8 m×1.2 m×0.1 m;上圆柱的直径为0.1 m,高0.3 m;下圆柱的直径为1 m,高0.5 m;与下圆柱相连接的半球高0.3 m.卫星的结构如图4所示.

图4 卫星结构图 (a)卫星结构简图;(b)面元划分后的结果Fig.4.Structure of satellites:(a)Satellite structure;(b)result of surface division.

仿真计算过程中,卫星本体表面包覆黄色热控材料,帆板正面贴满太阳能电池片,帆板背面涂有机黑漆,本体上下两个组件涂有机白漆.对于有机黑漆(反射率0.04)和有机白漆(反射率0.9),其反射特性遵从朗伯漫反射定律.黄色热控材料和电池片具有较强的镜反射特性,利用双向反射分布函数进行描述.对于黄色热控材料及太阳电池片,其BRDF的测量及建模方法见文献[5],采用改进的Sun模型对太阳电池片的BRDF进行描述,能够较好地拟合测量数据,如图5所示.图5给出了入射角为30°,45°,60°时测量得到的双向反射分布函数与反射角的关系.可以看出,两种材料表现出了很强的镜反射特性,反射能量主要集中在镜面反射方向±10°的范围.

图5 (网刊彩色)黄色热控材料与太阳电池片的BRDF测量与建模结果 (a)入射角30°;(b)入射角45°;(c)入射角60°Fig.5.(color online)BRDF measurement and modeling results of yellow thermal control materials and solar cells:(a)Incident angle of 30°;(b)incident angle of 45°;(c)incident angle of 60°.

观测卫星探测系统参数如表1所列.

选择两种轨道参数作为仿真参数,表2表示近轨,表3表示观测星与目标星同轨.

表1 观测星探测系统参数Table 1.Parameters of observing satellite detection system.

表2 目标卫星与观测卫星的轨道参数1Table 2.Orbital elements 1 of the target satellite and observation satellite.

表3 目标卫星与观测卫星的轨道参数2Table 3.Orbital elements 2 of the target satellite and observation satellite.

4.1.2 仿真结果

对于表2的轨道参数进行仿真得到的空间目标时序光度信号如图6所示,本体和帆板的反射峰值均在38帧,帆板能量占较大比例.

针对表3中的轨道参数仿真得到的目标整体时序光度信号如图7(a)所示,图中信号有两个峰值,分别为卫星帆板和本体的反射峰值.帆板的镜面反射峰值在30帧,本体反射峰值点在42帧;此时,光度信号呈现了多峰特性.目标本体和帆板的时序光度信号如图7(b)和图7(c)所示.

4.2 光学特性宏观模型的重构

4.2.1 轨道参数1对应的重构结果

针对上述卫星模型、卫星轨道参数1、探测器参数仿真所得到的目标时序光度曲线,进行基于双面法的多级BRDF反演,结果列于表4.根据反演的结果,分别对本体、帆板、目标整体的光度信号进行重构,结果如图8所示.

由图8可见,观测星与目标星在同一轨道时,时序光度信号呈现单峰特性,重构与实测数据整体拟合效果很好,但帆板、本体重构误差较大,主要原因为帆板和本体信号发生耦合,本文方法对混合信号曲线中各组分的特征辨识度低,无法精确分离和提取各个信号,故在同轨观测情况下,算法误差较大.

表4 卫星轨道1对应的参数反演结果Table 4.Parameters inversion results on orbit 1.

4.2.2 轨道参数2对应的重构结果

针对上述卫星模型、卫星轨道参数2、探测器参数仿真所得到的目标时序光度曲线,进行基于双面法的多级BRDF反演,结果列于表5.根据反演的结果,分别对本体、帆板、目标整体的光度信号进行重构,结果如图9所示.

图6 轨道1目标时序光度信号 (a)整体时序光度信号;(b)本体时序光度信号;(c)帆板时序光度信号Fig.6.Space object photometric sequence signal on orbit 1:(a)Photometric signal of space object;(b)photometric signal of body;(c)photometric signal of panel.

表5 卫星轨道参数2对应的参数反演结果Table 5.Parameters inversion result on orbit 2.

图7 轨道2目标时序光度信号 (a)整体时序光度信号;(b)本体时序光度信号;(c)帆板时序光度信号Fig.7.Space object photometric sequence signal on orbit 2:(a)Photometric signal of space object;(b)photometric signal of body;(c)photometric signal of panel.

由图9可见,观测星和目标星处于相近轨道时,时序光度信号为多峰特性,分别体现了本体和帆板分别占优时的特性,对于二者的信号解混起到了决定性的作用,此时针对二者的参数反演及信号重构精度较高.

4.2.3 重构结果的对比分析

重构精度用确定系数(R-square)来表征,确定系数是通过数据的变化来表示两组数据的相关性,取值范围在0—1之间,越接近1,表明数据相关性越好,对原始数据的解释能力越强;确定系数表达式如下:

图8 (网刊彩色)轨道1对应的时序光度信号重构结果(a)目标整体重构结果;(b)本体重构结果;(c)帆板重构结果Fig.8.(color online)Reconstructing photometric sequence signal for orbit 1:(a)Reconstructing result of space object;(b)reconstructing result of body;(c)reconstructing result of panel.

其中SSE为拟合数据与原始数据对应点的误差的平方和,SST为原始数据和均值之差的平方和.

经过计算轨道参数1对应的目标卫星整体、帆板、本体的R-square分别为0.99913,0.97883,0.86196.轨道参数2对应的目标卫星整体、帆板、本体的R-square分别为0.99753,0.99978,0.97304.

对于轨道1,帆板和本体的确定系数均低于轨道2的反演结果,但是目标整体光度信号的确定系数却高于轨道2的反演结果.说明该反演方法对于目标星和观测星在同一轨道的情况下,帆板本体信号的分离有一定误差,但对目标整体光学模型具有很强的描述能力.

图9 (网刊彩色)轨道2对应的时序光度信号重构结果(a)目标整体重构结果;(b)本体重构结果;(c)帆板重构结果Fig.9.(color online)Reconstructing photometric sequence signal for orbit 2:(a)Reconstructing result of space object;(b)reconstructing result of body;(c)reconstructing result of panel.

对于轨道2,帆板的重构精度最高,几乎与原始数据重合.重构精度高是因为帆板的主要材料为太阳能电池片,且没有其他构件,因此光学特性显著,可以进行精确的反演与重构.本体的重构精度较帆板略低,由图4卫星结构可知,卫星本体除了一个立方体构件,还有其他构件.本文反演采用双面法理论,即将本体、帆板各视为一个面,所以本体的其他构件对本体重构的精度有一定影响,但重构误差在可允许的范围内.

综上,采用基于双面法的多级BRDF融合模型对空间目标光学特性宏观模型具有很强的描述能力,当目标信号有多峰特性的条件下,能够进行精准的卫星帆板、本体信号分离;在无多峰特性时,能够较为精确地重现目标整体的光度信号.

5 结 论

针对在轨非合作空间目标的光学特性反演问题,基于在轨观测过程的双面假设以及BRDF多级表征,提出了基于时序光度信号分析的光学特性模型参数反演及模型重构方法.利用在轨空间目标光度模型进行了时序光度信号的仿真,结果表明提出的宏观表征模型的在轨重构方法能够实现近轨条件下针对本体、帆板97%以上的信号重构精度.该方法可为天基平台空间目标光学态势感知提供一种解决途径,并可为空间目标姿态、形状反演提供技术支撑.

[1]Yuan Y,Sun C M,Zhang X B 2010Acta Phys.Sin.59 2097(in Chinese)[袁艳,孙成明,张修宝 2010物理学报59 2097]

[2]Wang H Y,Zhang W,Dong A T 2013Measurement46 3654

[3]Hou X Y,Zhi X Y,Zhang H L,Zhang W 2014Proc.SPIE9299 929914

[4]Yuan Y,Sun C M,Zhang X B,Zhao H J,Wang Q 2010Acta Opt.Sin.30 2748(in Chinese)[袁艳,孙成明,张修宝,赵慧洁,王潜2010光学学报30 2748]

[5]Wang H Y,Zhang W 2012J.Mod.Opt.59 547

[6]Wang H Y,Zhang W,Dong A T 2012Appl.Opt.51 7810

[7]Wang H Y,Zhang W,Wang F G 2012Sci.China:Technol.Sci.55 982

[8]Linares R,Shoemaker M,Walker A,Mehta P M,Palmer D M,Thompson D C,Koller J,Crassidis J L 2013Proceedings of the Advanced Maui Optical and Space Surveillance Technologies ConferenceMaui,Hawaii,United States,September 10–13,2013 p1889

[9]Coughlin J 2014Proceedings of the Advanced Maui Optical and Space Surveillance Technologies ConferenceMaui,Hawaii,United States,September 9–12,2014

[10]Calef B,Africano J,Birge B,Hall D,Kervin P 2006Proc.SPIE6307 6307E

[11]Hinks J C,Linares R,Crassidis J L 2013AIAA Guidance,Navigation,and Control(GNC)ConferenceBoston,MA,United States,August 19–22,2013 p1

[12]Wetterer C J,Jah M 2009J.Guid.Control.Dynam.32 1648

[13]Holzinger M J,Alfriend K T,Wetterer C J,Luu K K,Sabol C,Hamada K 2014J.Guid.Control.Dynam.37 921

[14]Gou X R,Du X P,Liu H 2016Laser Opt.Progress.53 1(in Chinese)[苟瑞新,杜小平,刘浩 2016激光与光电子学进展53 1]

[15]Zhuang X X,Ruan N J,Zhao S S 2016Infrared Laser Eng.45 201(in Chinese)[庄绪霞,阮宁娟,赵思思 2016红外与激光工程45 201]

[16]Du X P,Liu H,Chen H,Gou X R,Du M J 2016Acta Opt.Sin.36 251(in Chinese)[杜小平,刘浩,陈杭,苟瑞新,杜明江2016光学学报36 251]

[17]Sun C M,Zhao F,Yuan Y 2015Acta Phys.Sin.64 034202(in Chinese)[孙成明,赵飞,袁艳 2015物理学报64 034202]