基于卷积曲面的动态实时星图模拟

闫劲云， 刘慧， 赵伟强， 江洁

(1. 中国计量科学研究院, 北京 100029； 2. 北京航空航天大学仪器科学与光电工程学院, 北京 100083)

星敏感器是航天器中最重要的姿态敏感设备，在设计、调试和测试等环节中，需要提供输入激励以完善和评估其各项功能和性能，星模拟器是实现这一任务的最佳途径[1-2]。星模拟器包括单星模拟器、静态多星模拟器、动态多星模拟器[3-4]。动态多星模拟器可模拟全天区任意指向的星图，然后由光学系统将星图输出为适合星敏感器接收的星光信号[5]。

随着星敏感器的动态性能和数据更新率等各种指标的不断提高，对多星模拟器的要求也随之提高[6-7]。星敏感器的动态性能是指当其随着载体运动时，星点成像发生拖尾的情况下，保持正常姿态输出的能力。为了测试动态性能，要求多星模拟器可以模拟拖尾的星点，目前多星模拟器对于动态拖尾星点模拟的算法研究还较少[5,8]。当前星敏感器的数据更新率普遍可达10 Hz，一些星敏感器可达30 Hz[6]，多星模拟器要求动态刷新显示变化的星空，其刷新率不低于星敏感器的数据更新率，通常大于30 Hz，因此对算法的实时性也提出了较高要求[9]。

多星模拟器算法的核心任务是高精度实时星图生成，主要包括导航星搜索、星点坐标变换、星点光斑弥散及像素离散化等步骤。导航星搜索是指根据星敏感器的视轴指向在全天区范围内搜索当前视场内的星点信息。星点坐标变换将搜索到的星点信息变换为模拟星图中的图像坐标和灰度信息。星点光斑弥散模拟了镜头的弥散作用，将星点变换为星点光斑[10]。像素离散化过程将连续的星点光斑变为离散的数字图像。通过星表分区[3,5]和重新组织构建星表[11]等方法，可以解决搜索导航星耗时较高的问题，本文对此不再赘述。静态下的星点坐标变换模型较为成熟，但在动态条件下，星点在模拟星图中的坐标构成一条轨迹。刘海波[8]和刘朝山[12]等推导了星点在像面的速度与星敏感器角速度的关系，但未给出星点轨迹曲线的表达式。Sun等[7]在推导过程中使用了一阶近似，因此仅适用于曝光时间较短的情况。本文在此基础上，给出了完备的星点轨迹模型。对于模拟星点光斑弥散过程，Sun等[7]利用运动模糊图像中的概念，将动态拖尾星点在频域中表示为镜头弥散函数和运动退化函数的乘积，但运动退化函数在时域中的物理含义仍有待进一步明确。Yan等[6]利用线扩散函数的概念，得出了拖尾星点光斑的解析表达式模型，但仅适用于星点在像面做匀速直线运动的情况。本文针对这些问题，建立了拖尾星点的卷积曲面模型，更准确地描述了动态拖尾星点的成像过程。对于像素离散化过程，由于涉及对光斑曲面的二重积分，普遍存在计算复杂、实时性较差的问题[13]。Yan等[6]的方法仅适用于匀速直线运动的情况，刘朝山等[12]的方法假定星点扫过每个像素的时间一致，因此不适合非匀速运动的情况。本文在卷积曲面模型的基础上提出了一种快速像素离散算法，适用于星点做任意运动的情况，并且大大加快仿真速度，以满足动态多星模拟器的实时星图仿真要求。

1 动态拖尾星点模拟

1.1 星点运动轨迹

当星敏感器随载体发生角位移时，在曝光时间内，星点在星敏感器像面移动成像，星点中心的运动轨迹为一条曲线段。

设星表中恒星i在天球坐标系下的映像为Vi。在t时刻，此恒星在星敏感器坐标系下的映像为Wi(t)，则Vi和Wi(t)可分别表示为

(1)

(2)

式中：αi和δi分别为赤经和赤纬；(xi(t),yi(t))为星点在星敏感器像面的投影坐标；F为星敏感器的镜头焦距。

Wi(t)和Vi是同一个矢量在不同坐标系下的映像，二者关系为

Wi(t)=R(t)Vi

(3)

式中：R(t)为星敏感器的姿态矩阵。

设星敏感器相对于天球坐标系的旋转角速度为ω=[ωx,ωy,ωz]T，式(3)对时间t求导得[8]

-[ω×]Wi(t)

(4)

(5)

若曝光起始时刻t0星点在像面坐标为(xi0,yi0)，则动态下星点在像面的运动轨迹C为

(6)

1.2 基于卷积曲面的光斑弥散模型

星敏感器为了获得亚像素级的定位精度，无穷远处的星光经过星敏感器镜头在图像传感器像面并非汇聚于单个像素，而是成像为一个弥散光斑。静态下星点光斑能量分布可表示为[6-7]

f(x,y)=ΦTh(x-xc,y-yc)

(7)

式中：Ф为星光被镜头收集并入射到像面的辐射通量；T为曝光时间；(xc,yc)表示星点光斑中心位置；h(x,y)为镜头的点扩散函数，通常以二维高斯函数来描述[6]。

(8)

(9)

式中：E0为0等星在地球大气外表面处的辐照度[1]，E0= 2.96×10-14W/mm2；m为入射星光的星等；D为星敏感器的光学镜头孔径；τ为镜头的透射率；ρ为镜头的弥散半径。

动态情况下，星点在像面发生移动，星点光斑的能量分布构成一个二维卷积曲面[14-15]，如图1所示。图1中，星点的轨迹曲线为C，r为轨迹曲线 上的点，r=(xc(t),yc(t))，其参数方程见式(6)。

图1 拖尾星点光斑成像过程Fig.1 Imaging process of smeared star spot

p为像平面内的任意一点，p=(x,y)。h(p)为星敏感器光学成像系统的点扩散函数，h(|p-r|)表示此函数的中心位于点r。记权重函数g(r)表示星点经过轨迹曲线上的弧微元|dr|所需的时间，则有

(10)

在点r附近取弧微元|dr|，星点扫过此微元期间的能量为Φg(r)|dr|。显然，当星点的速度越快，星光累积的能量越少；反之，当星点的速度越慢，能量累积越多。由于光学镜头的弥散作用，此星点能量按照点扩散函数描述的方式弥散到了像面，能量分布可表示为Φg(r)h(|p-r|)|dr|，那么沿曲线积分后可得星点在像面的总能量分布为

Φg(p)⊗h(p)

(11)

式中：⊗为卷积运算符。星点能量分布为一卷积曲面，即f(p)为点扩散函数h(p)与权重函数g(r)的卷积。拖尾星点光斑的卷积曲面模型具有清晰的物理含义，但不便于运算，为此进一步导出动态下星点能量分布积分模型为

(12)

式(11)和式(12)所表示的模型精确描述了动态条件下的拖尾星点成像过程。

1.3 卷积曲面的快速像素离散算法

入射到像面的光斑能量，在像元内转化为电荷信号，再经过读出和放大，转换为数字图像。那么位于第i行和第j列的像素的图像灰度值为[1,5-8,12-13]

(13)

式中：K为转换系数，是图像传感器的量子效率和转换增益的乘积。图像传感器对光斑能量的作用不仅仅体现为离散采样，还包括积分作用，此过程称为像素离散。显然，像素灰度值为一个三重积分表达式，因此仿真一幅星图需要对每个星点光斑范围内的每一个像素计算三重积分，运算的时间开销较大，不能满足动态多星模拟器实时性的要求。

由拖尾星点的卷积曲面模型的物理含义出发，提出一种快速算法，以解决运算量较大的问题。具体思路为：将曝光时间T等间隔地分为N段，每段时间间隔ΔT=T/N。相应的，轨迹曲线也被分为了N段。当ΔT足够短时，地轨迹上的每一段可近似为一点，此时的星点光斑成像近似为静态下的情况。静态下对二维高斯函数的二重积分可推导得其原函数，将耗时的数值积分转化为函数求值。最后将这N点沿轨迹累加，即可得到

(14)

(15)

式中：(xck,yck)为第k时间段时星点光斑中心在轨迹上的坐标；erf(x)为误差函数。式(15)将复杂的三重积分化简为函数求值、乘法和加法运算，大大提高了拖尾星点的仿真速度。

2 仿真及分析

仿真计算机的配置如下。CPU：Inter Core i7，内存8 GB，仿真环境：MATLAB 2015b。仿真星敏感器的主要参数为：视场20°×20°，镜头直径27 mm，焦距37.92 mm，像素尺寸5.5 μm × 5.5 μm，图像分辨率2 048像素×2 048像素。采用SAO星表，选取亮于6.0等星的星作为仿真星表。

2.1 精度对比

本节首先对比3种星点仿真算法的精度，再研究对轨迹划分到何种精度，本文算法的精度可以达到满意结果。式(13)表示的三重积分法是未优化的原始算法，速度慢，但精度高，得到的星点光斑可认为是真值。本文算法的精度依赖于对积分划分的细度，随着N增大，精度逐渐逼近真值。第3种算法为运动模糊图像卷积法，拖尾星点由原始星点图像与运动模糊核函数图像卷积而得[7,10]。定义星点图像灰度误差为

ferr=|E(j,i)-I(j,i)|

(16)

进一步定义最大灰度误差为星点区域内误差绝对值最大的灰度值；定义灰度RMS误差为星点区域内每一个像素的误差灰度值的均方差。

设[ωx,ωy,ωz]=[5(°)/s, 5(°)/s,0]，曝光时间T=20 ms，镜头的弥散半径为0.7像素，星点初始位置为像面中心，则由式(6)可得，星点拖尾长度为14.33像素。采用三重积分法得到的拖尾星点如图2(a)所示；采用本文算法，将星点轨迹划分为143段，仿真得到的拖尾星点如图2(b)所示；采用图像卷积法得到拖尾星点如图2(c)所示；后两者的星点图像灰度误差分别如图2(d)、(e)所示。为便于比较，图中星点灰度值在最大值处进行了归一化。

对比图2(a)、(b)可知，本文仿真结果与真值较为一致。由图2(d)可知，仅起始点和终点附近存在一定误差，星点拖尾中间部位的误差可忽略，其灰度RMS误差为0.005 4。由图2(a)、(c)可知，采用图像卷积法得到的拖尾星点，其星点能量弥散程度较大，图2(e)所示的星点误差也较大，其灰度RMS误差为0.102 0。

由式(14)可知，当轨迹划分越细，则仿真精度越高。本文算法在不同划分间隔下，星点灰度误差如表1所示。

由表1可知，随着轨迹划分间隔增加，星点灰度误差逐渐增大。由于成像过程中存在的背景噪声及其他噪声，星点图像灰度存在起伏，若仿真星点最大误差小于噪声引起的灰度起伏，则仿真误差对最终结果可忽略。若以最大灰度误差为限，轨迹划分间隔为0.05像素时，仿真星点精度满足要求；若以RMS误差为限，轨迹划分为0.10像素即可满足要求。

图2 拖尾星点光斑仿真精度对比Fig.2 Simulation accuracy comparison of smeared star spots

划分间隔/像素0.020.050.100.20最大灰度误差0.00530.01290.02600.0525灰度RMS误差0.00110.00270.00540.0109

2.2 仿真速度对比

令多星模拟器产生1 000个随机视轴指向，分别利用三重积分法、运动模糊图像卷积法和本文算法仿真生成拖尾星图，统计每幅仿真星图的仿真时间，3种算法仿真耗时的直方图统计如图3所示。

由图3可知，三重积分法仿真一幅星图所需时间为几秒，运动模糊图像卷积法所需时间为几百毫秒，本文算法耗时为几十毫秒，3种算法仿真一幅星图的平均耗时分别为3 030、309和22 ms。显然，本文算法的仿真速度大大优于前两种算法。

图3(c)所示的个别星图的仿真耗时大于33 ms，不满足多星模拟器刷新率的要求。为此，绘制模拟星图仿真耗时与星点数量的关系图(见图4)，以进一步分析超时原因。

由图4可知，星图仿真耗时与星点数量是线性关系。星点数量较多时，仿真耗时也越多。对于单个拖尾星点，平均仿真耗时为0.45 ms，那么为保证单幅星图仿真时间小于33 ms，只要星图中 星点个数小于70即可。为此，对于模拟星图中星点数量大于70颗星的情况，根据星点的亮度信息，筛选其中较亮的70颗星，其余亮度较暗的星予以剔除。这样的处理对实际使用影响较小，原因在于星敏感器为获得较高的姿态精度，仅采用较亮的若干颗星进行姿态解算。较暗的星形成的星点光斑图像的信噪比低，因此对其定位的精度也较低，应予以剔除。增加对多余暗星剔除的步骤，得到的仿真耗时直方图如图5所示。

图3 不同算法仿真耗时直方图Fig.3 Histogram of simulation time consumption for different algorithm

图4 星图仿真耗时与星图中星点数量的关系Fig.4 Star map simulation time consumption vs. star number of star image

由图5可知，剔除多余暗星后，本文的快速星点仿真算法可满足动态多星模拟器的实时星图仿真的要求。

图5 剔除暗星后仿真耗时直方图Fig. 5 Histogram of simulation time consumption after dark star is excluded

3 结 论

针对动态多星模拟器对拖尾星点仿真难以满足实时性的问题，本文首先建立了动态拖尾星点的卷积曲面成像模型，然后提出了一种基于卷积曲面的快速像素离散算法。得到的仿真星点误差低于图像噪声水平，仿真精度高；同时仿真一幅星图耗时小于33 ms，比现有算法速度提高14倍，满足了多星模拟器刷新率30 Hz的要求。