CN110455281A - 暗弱小天体光学导航特征成像模拟方法 - Google Patents

Abstract

本发明公开的暗弱小天体光学导航特征成像模拟方法，属于深空探测领域。本发明通过构建深空星图的成像CCD模型、噪声模型及星点光度模型，实现暗弱小天体探测器任务远距离接近段长曝光下拖曳星图的成像模拟。在给定的光学系统参数和背景噪声参数下，首先生成噪声背景，然后确定能够进入视场的导航星及导航星在CCD平面上的位置并计算各星点能量，最后随机生成一系列光轴指向扰动，计算由该扰动产生的星图拖曳轨迹的能量，生成拖曳星图。本发明能够模拟高保真的深空暗弱小天体探测器长曝光星图，能够在没有相关器件配套的情况下进行图像处理算法的测试及确定实际深空探测任务的光学系统参数。

Description

暗弱小天体光学导航特征成像模拟方法

技术领域

本发明涉及一种暗弱小天体光学导航特征成像模拟方法，尤其涉及一种小天体探测器远距离接近段长曝光拍摄暗弱小天体目标时的成像模拟方法，属于深空探测领域。

背景技术

强自主、高精度且实时的自主光学导航是深空探测任务中的重要技术环节。在小天体探测任务的远距离接近段，光学导航目标一般为目标小天体光心。目标小天体的星等一般都较低，需要进行长时间曝光才能够成像，在曝光时间内，热力缩涨，机械振动，姿态偏移等多种原因会引起成像系统光轴在曝光时间内发生多次扰动，导航目标在CCD平面上的成像会发生偏移，成像将不再成点状或斑状分布而是产生拖曳，实际成像为轨迹，难以提取导航目标的特征信息进行精确导航。对深空导航图像展开研究具有非常重要的意义，为了验证相关特征提取算法和地面试验的顺利展开，研究暗弱小天体光学导航特征成像模拟方法极为关键。

目前，已有一些长曝光星图的模拟方法，但主要应用在星敏感器上，详细考虑长曝光情况下噪声分布和信号能量分布的接近段目标小天体成像模拟方法还未曾出现过。

发明内容

本发明公开的暗弱小天体光学导航特征成像模拟方法要解决的技术问题是：模拟暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图，具有如下3个优点：(1)通过构建深空星图的成像CCD模型、噪声模型及星点光度模型，能够模拟高保真的深空暗弱小天体探测器长曝光星图；(2)能够在没有相关器件配套的情况下验证轨迹图像星点质心提取的图像处理算法的有效性，进行图像处理算法的测试；(3)通过模拟不同CCD硬件参数及拍摄条件下的长曝光星图，有助于确定实际深空探测任务的光学系统参数。

所述的暗弱小天体指视星等远低于背景导航恒星的小天体目标。

所述远距离接近段指探测任务目标在探测器光学系统能够成像但成像大小小于等于3个像素的阶段。

所述的长曝光拖曳图像是由成像系统光轴在曝光时间内发生多次扰动所造成的。所述的引起成像系统光轴在曝光时间内发生多次扰动的原因包括热力缩涨、机械振动、姿态偏移。

本发明的目的是通过下述技术方案实现的。

本发明公开的暗弱小天体光学导航特征成像模拟方法，通过构建深空星图的成像CCD模型、噪声模型及星点光度模型，实现暗弱小天体探测器任务远距离接近段长曝光下拖曳星图的成像模拟。在给定的光学系统参数和背景噪声参数下，首先生成噪声背景，然后确定能够进入视场的导航星及导航星在CCD平面上的位置并计算各星点能量，最后随机生成一系列光轴指向扰动，计算由该扰动产生的星图拖曳轨迹的能量，生成拖曳星图。本发明能够模拟高保真的深空暗弱小天体探测器长曝光星图，能够在没有相关器件配套的情况下进行图像处理算法的测试及确定实际深空探测任务的光学系统参数。

本发明公开的暗弱小天体光学导航特征成像模拟方法，包括如下步骤：

步骤一、给定暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图的光学系统参数。

给定暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图的光学系统参数，所述光学系统参数包括CCD相机的光圈半径A、视场角FOV、像素尺寸px/py、拍摄曝光时间T、光学效率η、量子效率Q和ADC增益K。

步骤一的具体实现方法为：

CCD相机的感应器件将入射光子转化为电子并将生成的电子贮存起来，随后，读入器件测量每个像素中的电荷量在电容器两端产生模拟电压，模数转换器将每个像素中的电压数字化，产生一个计数写入电脑磁盘中；最后将电脑磁盘中的计数转化为灰度，生成灰度图像。

每个像素点接收到的入射光子数决定该像素点处信号的大小。当光学系统入瞳处的光通量为F时，CCD处接受到的光子数P为：

其中，为平均波长，h为普朗克常数，c为光速，A为光圈半径，η为光学效率，。

入射光子数P与生成光电子C的比率为光学系统的量子效率Q：

C＝Q*P (2)

光电子数C与计数值N之间的关系由ADC的增益K所确定：

N＝K*C (3)

星载光学系统生成的图像为8位灰度图像，有256级灰度，计数N和灰度DN间的转换关系为线性关系：

DN＝a*N+b (4)

步骤二、给定暗弱小天体探测器拖曳星图的背景噪声参数，并根据给定的背景噪声参数生成噪声背景。

给定暗弱小天体探测器拖曳星图的背景噪声参数，所述背景噪声参数包括读数噪声偏置值λ_bias、暗电流噪声泊松参数λ_dark、成像背景亮度Sky_level，并根据给定的背景噪声参数生成噪声背景。

步骤二的具体实现方法为：

总噪声包括读数噪声、暗电流噪声、成像背景亮度噪声：

S_噪＝I_background+I_read+I_dark (5)

由于深空中星际尘埃、太空碎片及其他物质的影响，拍摄的拖曳星图中存在大量杂散光，背景并不是完全的黑色。模拟星图的背景由用户输入参数背景亮度Sky_level和曝光时间T所决定，服从泊松分布，背景信号强度写为：

I_background～P(Sky_level*T) (6)

读数噪声是读数器件自身的电压所带来的噪声，所述读数噪声服从关于偏置值的高斯分布，其信号强度写为：

暗电流噪声是由CCD的结构内热造成的，其在单位时间内的噪声水平与CCD温度相关，属于散粒噪声的一种，服从泊松分布，其信号强度写为。

I_dark～P(λ_dark) (8)

根据给定的背景噪声参数生成噪声背景。

步骤三、根据目标小天体所在位置确定光轴指向，检索星表确定在该指向下能够进入视场范围的导航恒星列表，并计算所述导航星在CCD平面上的位置。

步骤三的具体实现方法为：

定义光学相机的视轴方向为(α₀,δ₀)，则出现在光学相机视场内的导航目标所需满足的赤经赤纬条件如下：

其中，(α_i,δ_i)为导航恒星的赤经赤纬，FOV为视场角大小。

考虑天球坐标系和相机坐标系之间的坐标变换。由于导航星距离遥远，将相机坐标系的原点O和天球坐标系的原点O'视为重合，忽略坐标系平移所带来的误差,仅需要考虑旋转变换。导航恒星在天球坐标系下的坐标为(α_i,δ_i)，天球赤道坐标系O-UVW和相机坐标系O'-XYZ间的转换关系为：

[X Y Z]^T＝M[U V W]^T (10)

[U V W]^T＝[cosα_icosδ_i sinα_icosδ_i sinδ_i]^T (11)

所述O'-XYZ坐标系的Z轴为视轴正向，X轴和Y轴方向与CCD相机平面方向一致。

其中，j角与CCD相机安装方向有关，取j＝90°，θ为视轴与OW轴的夹角，O为O'X轴与OU轴的夹角，θ＝90°-δ₀，O＝90°+α₀。因此，M矩阵简化为：

最后，考虑相机坐标系和CCD成像平面坐标系之间的坐标变换。根据光学相机的成像模型，得到相机坐标系和CCD平面坐标系之间的几何关系：

上式中，(x,y)为导航目标在CCD成像平面投影的坐标，f为相机焦距，px,py为像素的高和宽。

同时，根据CCD相机视场角大小、传感器尺寸和焦距间的关系知：

f＝2pxNx^-1tan(FOV/2)＝2pyNy^-1tan(FOV/2) (14)

其中，Nx,Ny是CCD平面尺寸，FOV为视场角大小。

最终得导航目标在CCD成像平面投影的坐标位置的表达式：

所述导航目标在CCD成像平面投影的坐标位置即为导航星在CCD平面上的位置。

步骤四、给定PSF能量扩散半径，根据背景导航星及目标小天体的视星等调用高斯扩散公式得到各星点散焦后的能量。

步骤四的具体实现方法为：

在天文观测中，视星等与单色光流间的关系为：

m_v＝-2.5log₁₀F_v+c (16)

其中，m_v为信源视星等，F_v为信源单色光流，c为常数。

由公式(17)得出，信源的成像强度主要由信源的星等决定，信源视亮度与信号强度之间的关系为：

其中，I为某像素点上的光能，I_m为0等星所对应的光能，mag为某像素点对应的星等。

由于深空导航目标为点光源，需要考虑在CCD焦平面上造成的点扩散问题，在输出的数字灰度图上实际表现为包含若干像素点的光斑。采用Gauss点扩散函数能够分别计算出目标及导航星在CCD面的投影位置周围像素收集的能量：

上式中，(m,n)为所计算的像素点坐标，I(m,n)为所计算的像素点收集的光能，I₀为目标点光源的光能，σ_PSF为Gauss弥散半径，(x₀,y₀)为目标点光源在CCD面的投影坐标。

根据公式(18)能够分别计算各星点散焦后的能量。

步骤五、随机生成一系列光轴指向扰动，计算由该扰动产生的星图拖曳轨迹的能量，生成拖曳星图，即实现模拟暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图。

步骤五的具体实现方法为：

由于光轴扰动，在单次曝光时间T内相机光轴发生k次变化，每次扰动后的天球赤道坐标系和相机坐标系之间转换矩阵记为M_k，扰动后导航目标投影的位置坐标：

上式中，Nx,Ny是CCD平面尺寸，FOV为视场角大小，M_k1,M_k2,M_k3为M_k各行，X_i为导航恒星在天球直角坐标系下的位置。

由于两次扰动之间的时间间隔较短，故将两次光轴方向变化间的运动视为是匀速的，在任意时刻t有：

M_t＝M_k·(I+ω·Δt) (20)

其中，M_t为任意时刻t的转换矩阵，ω为光轴的方向变换角速度，Δt为时间步长。

在任意时刻t，将M_t作为天球赤道坐标系和相机坐标系间的状态转移矩阵M，按照步骤三和步骤四计算任意时刻t星点的位置和能量。将不同时刻下的星点能量叠加即得到暗弱小天体探测器拖曳星图。

还包括步骤六：将步骤五得到暗弱小天体探测器拖曳星图用于验证轨迹图像星点质心提取的图像处理算法的有效性，进行图像处理算法的测试。

步骤六具体实现方法为：

步骤6.1：对长曝光下产生的拖曳星图进行去噪预处理。

考虑到图像噪声的非高斯性及低信噪比性，采用阈值小波变化的方法进行去噪。小波去噪方法的思路是：进行小波变换后，信号主要集中在低频，而噪声集中在高频上，故给定适当的阈值除去噪声。

采用软阈值小波方法进行去噪，软阈值是指若小波变换系数比设置的阈值小则将该系数置零，若小波变换系数比该阈值大则将该系数向阈值方向收敛，小波变换系数表达形式为：

其中，ω_i,j为小波变换系数，λ为设定的阈值。

步骤6.2：通过连通域判定轨迹所在的兴趣区域(ROI)，通过兴趣轨迹选取准则选取兴趣轨迹，采用逐像素细化算法(One Pass Thinning Algorithm)对兴趣轨迹进行细化并提取出轨迹骨架。

长曝光拖曳图像模型表示式为：

I_Blur＝K·I_static+E (22)

其中，I_Blur是长曝光的拖曳星图，I_static是需要被复原的静态星图，E为噪声项，K为是每一像素的模糊转移矩阵。当星图中各像素发生的模糊是均匀时，则式(2)变换为卷积的形式，如式(3)所示：

I_Blur＝k*I_static+e (23)

其中，k被称为模糊核，与静态星图I_static的每一像素发生卷积运算使图像产生模糊，其物理含义是反映每一像素点的运动特性。

由于长曝光拖曳星图的拖曳是点状目标的运动模糊，故星点拖曳轨迹本身能够反映地曝光时间内的运动，而模糊核的本质就是像素点的运动特性，故采用星点拖曳轨迹来估计拖曳图像的模糊核。

通过连通域判定轨迹所在的兴趣区域(ROI)，通过兴趣轨迹选取准则选取兴趣轨迹。

作为优选，所述兴趣轨迹选取准则选为亮度兴趣轨迹选取准则，即选取亮度最亮的轨迹为兴趣轨迹，兴趣轨迹的选取准则为：

max average(G_{ROI_i}) (24)

其中，G_{ROI_i}是第i个兴趣区的灰度矩阵。

OPTA算法是一种基于形态学的单循环并行模板匹配算法，其优点是运算简单，计算速度很快，且能够较好地保留图像的连通性和形状特征。采用OPTA算法对兴趣轨迹进行细化并提取出轨迹骨架，具体实现方法如下：按从左到右，从上到下的方式，从最左上角的像素点开始抽取出每个灰度值不为0的像素点的邻域，将该邻域与消除模板中的任意一个消除模板匹配，若不能匹配则保留像素点的值，否则要再将该邻域与保留模板匹配，若能匹配则仍然保留像素点的值，若不能匹配则删去该像素点，即将像素点的灰度值置零。对所有像素点进行上述操作，保留下的像素值不为0的像素点即为提取出的轨迹骨架。

步骤6.3：以步骤6.2中提取出的轨迹骨架为模糊核的估计初值，通过稀疏先验信息构建图像复原算法的代价函数，对长曝光拖曳星图进行盲复原，得到长曝光拖曳星图复原图像。

复原图像问题需估计出合适的模糊核及复原图像，使得图像复原的代价函数尽量最小，代价函数形式为：

然而，复原图像问题本质上是NP(non-deterministic polynomial)问题，即由于信息的匮乏，会导致结果不唯一甚至病态。为解决星图模糊复原的问题，需要根据导航星图的实际特点引入相关先验信息建立正则惩罚项。

深空导航星图的重要特征是星点能量集中在极少数像素点中，通过步骤一去噪声预处理后，绝大部分像素点为黑色背景，故利用图像的灰度稀疏作为先验信息。同理，由于有值像素点较少，故图像的边缘像素点也就更少，所以图像的灰度梯度也是稀疏的。稀疏先验信息通过L₀范数构建，即

最终采用的星图复原代价函数的形式为：

其中，项是为保证估计出的模糊核的简单性。

将步骤6.2中提取出的轨迹骨架为模糊核的估计初值，以代价函数最小为指标估计复原图像；由估计出的复原图像更新模糊核；如此循环迭代直至达到给定的循环的终止条件，此时估计出的复原图像即为最终得到的长曝光拖曳星图的复原图像。

步骤6.4：对步骤6.3得到的复原图像进行连通域判定，并提取复原图像中各星点的灰度光心，完成长曝光拖曳星图特征提取。

对步骤6.3得到的长曝光拖曳星图复原图进行连通域判定，去除连通域小于预设像素数量的目标。

作为优选，步骤6.4中去除连通域小于预设像素数量的目标优选去除连通域小于2个像素的目标。

然后，根据导航星覆盖像素的像素值确定复原星图导航星光心所在的位置，即实现提取复原图像中各星点的灰度光心，完成长曝光拖曳星图特征提取，并验证轨迹图像星点质心提取的图像处理算法的有效性，即进行图像处理算法的测试。

作为优选，步骤6.4中计算灰度光心的方法采用灰度加权算法，灰度加权算法的公式为：

还包括步骤七：改变光学系统参数，通过步骤一至步骤五模拟不同光学系统参数下的暗弱小天体探测器拖曳星图，进而选定满足实际深空探测任务的光学系统参数，实现实际深空探测任务。

有益效果：

1、针对星载光学系统在长曝光情况下光轴受扰动的问题，本发明公开的暗弱小天体光学导航特征成像模拟方法，通过构建深空星图的成像CCD模型、噪声模型及星点光度模型，能够模拟高保真的深空暗弱小天体探测器长曝光星图。

2、本发明公开的暗弱小天体光学导航特征成像模拟方法，能够在没有相关器件配套的情况下验证轨迹图像星点质心提取的图像处理算法的有效性，进行图像处理算法的测试。

3、本发明公开的暗弱小天体光学导航特征成像模拟方法，改变光学系统参数，通过步骤一至步骤五模拟不同光学系统参数下的暗弱小天体探测器拖曳星图，进而选定满足实际深空探测任务的光学系统参数，实现实际深空探测任务。

4、本发明公开的暗弱小天体光学导航特征成像模拟方法，还公开轨迹图像星点质心提取方法，该方法把轨迹图像的特征提取问题分解为轨迹细化和运动模糊图像盲复原两个步骤，将拖曳轨迹的细化结果作为盲复原算法的迭代初值，能够显著提高盲复原算法的计算效率，从而有效地提高星图光心特征提取的速度具有提取精度高、计算速度快的优点。所述轨迹图像星点质心提取方法长曝光图像中准确提取出导航星的光心位置，能够弥补传统静态星图提取算法的不足，能够实现扰动情况完全不确知的拖曳星图特征提取。

附图说明

图1为本发明公开的暗弱小天体光学导航特征成像模拟方法流程图；

图2为本发明方法应用软件的操作界面；

图3为本发明方法实施例1中的曝光时间100ms下的星图噪声背景；

图4为本发明方法实施例1中的信源分布状况

图5为本发明方法实施例1中光学系统参数1下曝光时长为100ms时的静态星图

图6为本发明方法实施例1中光学系统参数1下曝光时长为2s时的拖曳星图

图7为本发明实施例1采用的细化算法的(a)邻域定义，(b)消除模板，(c)保留模板；

图8为本发明实施例1中的(a)原兴趣轨迹，(b)提取出的辅助细化轨迹；

图9为本发明实施例1中仿真条件1下的星图复原结果；

图10为本发明实施例1中仿真条件1下的光心提取误差；

图11为本发明实施例1中(a)光学系统参数2下曝光时长为2s时的拖曳星图，(b)光学系统参数3下曝光时长为2s时的拖曳星图。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。

实施例1：

为了验证发明的可行性，针对小天体目标的远距离接近段，选取Tycho-2星表为完备星表，极限星等为15，选取的小天体目标为Eros 433，仿真时刻为2018年7月14日00:00TT，该时刻Eros 433在天球坐标系下的赤经和赤纬位置为(0.0252,-0.3010)，视星等为14.2。

如图1所示，本实施例公开的暗弱小天体光学导航特征成像模拟方法，具体实现步骤如下：

步骤一、给定暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图的光学系统参数。

给定暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图的光学系统参数，设定光学系统参数1下：CCD相机的光圈半径A为60mm，视场角为1.0°，图像像素尺寸为512×512，像素尺寸为27μm×27μm，光学效率90％，量子效率30％和ADC增益为1e⁹。

步骤二、给定暗弱小天体探测器拖曳星图的背景噪声参数，并根据给定的背景噪声参数生成噪声背景。

给定暗弱小天体探测器拖曳星图的背景噪声参数，设定星图背景亮度为50/pixel/s，读数噪声偏置值为3e^-及暗电流噪声泊松参数为1e^-。总噪声包括读数噪声、暗电流噪声、成像背景亮度噪声：

S_噪＝I_background+I_read+I_dark (29)

由于深空中星际尘埃、太空碎片及其他物质的影响，拍摄的拖曳星图中存在大量杂散光，背景并不是完全的黑色。模拟星图的背景由用户输入参数背景亮度Sky_level和曝光时间T所决定，服从泊松分布，背景信号强度写为：

I_background～P(Sky_level*T) (30)

读数噪声是读数器件自身的电压所带来的噪声，所述读数噪声服从关于偏置值的高斯分布，其信号强度写为：

暗电流噪声是由CCD的结构内热造成的，其在单位时间内的噪声水平与CCD温度相关，属于散粒噪声的一种，服从泊松分布，其信号强度写为。

I_dark～P(λ_dark) (32)

根据给定的背景噪声参数生成噪声背景如图3所示。

步骤三、根据目标小天体所在位置确定光轴指向，检索星表确定在该指向下能够进入视场范围的导航恒星列表，并计算所述导航星在CCD平面上的位置。

当拍摄时刻Eros 433位于相机视场的中心，天球坐标系下相机光轴指向为(1°26'43.8″,-18°45'20″)。则出现在光学相机视场内的导航目标所需满足的赤经赤纬条件如下：

其中，(α_i,δ_i)为导航恒星的赤经赤纬，FOV为视场角大小。

考虑天球坐标系和相机坐标系之间的坐标变换。由于导航星距离遥远，将相机坐标系的原点O和天球坐标系的原点O'视为重合，忽略坐标系平移所带来的误差,仅需要考虑旋转变换。导航恒星在天球坐标系下的坐标为(α_i,δ_i)，天球赤道坐标系O-UVW和相机坐标系O'-XYZ间的转换关系为：

[X Y Z]^T＝M[U V W]^T (34)

[U V W]^T＝[cosα_icosδ_i sinα_icosδ_i sinδ_i]^T (35)

所述O'-XYZ坐标系的Z轴为视轴正向，X轴和Y轴方向与CCD相机平面方向一致。

其中，j角与CCD相机安装方向有关，取j＝90°，θ为视轴与OW轴的夹角，O为O'X轴与OU轴的夹角，θ＝90°-δ₀，O＝90°+α₀。因此，M矩阵简化为：

最后，考虑相机坐标系和CCD成像平面坐标系之间的坐标变换。根据光学相机的成像模型，得到相机坐标系和CCD平面坐标系之间的几何关系：

上式中，(x,y)为导航目标在CCD成像平面投影的坐标，f为相机焦距，px,py为像素的高和宽。

同时，根据CCD相机视场角大小、传感器尺寸和焦距间的关系知：

f＝2pxNx^-1tan(FOV/2)＝2pyNy^-1tan(FOV/2) (38)

其中，Nx,Ny是CCD平面尺寸，FOV为视场角大小。

最终得导航目标在CCD成像平面投影的坐标位置的表达式：

视场内的导航目标的分布如图4所示，共有27颗导航恒星，进入视场内的导航恒星的具体信息如表1所示。

表1 进入视场的导航恒星具体信息

步骤四、给定PSF能量扩散半径，根据背景导航星及目标小天体的视星等调用高斯扩散公式得到各星点散焦后的能量。

由于深空导航目标为点光源，需要考虑在CCD焦平面上造成的点扩散问题，在输出的数字灰度图上实际表现为包含若干像素点的光斑。采用Gauss点扩散函数能够分别计算出目标及导航星在CCD面的投影位置周围像素收集的能量：

上式中，(m,n)为所计算的像素点坐标，I(m,n)为所计算的像素点收集的光能，I₀为目标点光源的光能，σ_PSF为Gauss弥散半径，(x₀,y₀)为目标点光源在CCD面的投影坐标。

根据公式(18)能够分别计算各星点散焦后的能量。若不考虑光轴的转动，曝光时长为时100ms，生成的静态星图如图5所示，此时暗弱的小天体目标已经能够成像，但灰度值较低，占像素个数也较少，被噪声淹没。

步骤五、随机生成一系列光轴指向扰动，计算由该扰动产生的星图拖曳轨迹的能量，生成拖曳星图，即实现模拟暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图。

设定曝光视场为2s，在曝光时间内光轴指向发生5次偏移，偏移量的范围为[-1,1]×10^-3rad。

由于两次扰动之间的时间间隔较短，故将两次光轴方向变化间的运动视为是匀速的，在任意时刻t有：

M_t＝M_k·(I+ω·Δt) (41)

其中，M_t为任意时刻t的转换矩阵，ω为光轴的方向变换角速度，Δt为时间步长。

在任意时刻t，将M_t作为天球赤道坐标系和相机坐标系间的状态转移矩阵M，按照步骤三和步骤四计算任意时刻t星点的位置和能量。将不同时刻下的星点能量叠加即得到暗弱小天体探测器拖曳星图，生成的长曝光拖曳星图如图6所示。此时暗弱的小天体目标已经可以明显成像，但需要注意的是背景噪声也随曝光时间增加而增大，且部分导航恒星轨迹出现过曝。

步骤六：将步骤五得到暗弱小天体探测器拖曳星图用于验证轨迹图像星点质心提取的图像处理算法的有效性，进行图像处理算法的测试。

步骤6.1：对长曝光下产生的拖曳星图进行去噪预处理。

采用软阈值小波方法进行去噪，软阈值是指若小波变换系数比设置的阈值小则将该系数置零，若小波变换系数比该阈值大则将该系数向阈值方向收敛，小波变换系数表达形式为：

其中，ω_i,j为小波变换系数，λ为设定的阈值。设定软阈值λ为

其中，σ_noise为噪声方差，L和H为图像尺寸大小。

步骤6.2：通过连通域判定轨迹所在的兴趣区域(ROI)，通过兴趣轨迹选取准则选取兴趣轨迹，采用逐像素细化算法(One Pass Thinning Algorithm)对兴趣轨迹进行细化并提取出轨迹骨架。

通过连通域判定轨迹所在的兴趣区域(ROI)，通过兴趣轨迹选取准则选取兴趣轨迹。

所述兴趣轨迹选取准则选为亮度兴趣轨迹选取准则，即选取亮度最亮的轨迹为兴趣轨迹，兴趣轨迹的选取准则为：

max average(G_{ROI_i}) (44)

其中，G_{ROI_i}是第i个兴趣区的灰度矩阵。

选取最亮的迹作为兴趣轨迹，定义图像中每一个像素的邻域如图7(a)所示，采用较为通用的Chin-Wan细化算法的消除模板和保留模板进行计算，消除模板和保留模板如图7(b)和(c)所示。该星原轨迹图像和提取出的细化轨迹如图8所示，细化程序的运行时间为0.096993秒，细化后的轨迹能够较好地保留原轨迹的形状特征。

步骤6.3：以步骤6.2中提取出的轨迹骨架为模糊核的估计初值，通过稀疏先验信息构建图像复原算法的代价函数，对长曝光拖曳星图进行盲复原，得到长曝光拖曳星图复原图像。

复原图像问题需估计出合适的模糊核及复原图像，使得图像复原的代价函数尽量最小，代价函数形式为：

深空导航星图的重要特征是星点能量集中在极少数像素点中，通过步骤一去噪声预处理后，绝大部分像素点为黑色背景，故利用图像的灰度稀疏作为先验信息。同理，由于有值像素点较少，故图像的边缘像素点也就更少，所以图像的灰度梯度也是稀疏的。稀疏先验信息通过L₀范数构建，即

最终采用的星图复原代价函数的形式为：

其中，项是为保证估计出的模糊核的简单性。

将步骤二得到的细化轨迹作为模糊核初值代入到盲复原算法中，以代价函数最小为指标估计复原图像，代价函数为：

然后，由估计出的复原图像更新对模糊核的估计，同样以代价函数最小为指标，代价函数为：

如此循环迭代直至达到循环的终止条件，循环的终止条件是迭代次数或者误差界限。

本实施案例中，选取迭代次数为循环终止条件。仿真条件1参数如下：λ₁＝λ₂＝4e^-2，γ＝2，迭代次数为2次。仿真条件1下的复原星图如图9所示，盲复原程序的运行时间为10.560349秒。

步骤6.4：对步骤6.3得到的复原图像进行连通域判定，并提取复原图像中各星点的灰度光心，完成长曝光拖曳星图特征提取。

对步骤6.3得到的长曝光拖曳星图复原图进行连通域判定，去除连通域小于预设像素数量的目标。

步骤6.4中去除连通域小于预设像素数量的目标优选去除连通域小于2个像素的目标。

然后，根据导航星覆盖像素的像素值确定复原星图导航星光心所在的位置，即实现提取复原图像中各星点的灰度光心，完成长曝光拖曳星图特征提取。

步骤6.4中计算灰度光心的方法采用灰度加权算法，灰度加权算法的公式为：

仿真条件1下得到的复原星图中各导航星x轴和y轴方向的光心提取误差如图10所示，x轴方向平均误差为0.0466个像素，标准差为0.0272个像素，y轴方向平均误差为0.1445个像素，标准差为0.0423个像素，灰度光心提取精度很高，证明本发明算法具体优越性。

步骤七：改变光学系统参数，通过步骤一至步骤五模拟不同光学系统参数下的暗弱小天体探测器拖曳星图，进而选定满足实际深空探测任务的光学系统参数，实现实际深空探测任务。

设定不同光学系统参数如下。

光学系统参数2：CCD相机的光圈半径A为30mm，视场角为1.0°，图像像素尺寸为512×512，像素尺寸为27μm×27μm，光学效率90％，量子效率30％和ADC增益为1e⁹。

光学系统参数3：CCD相机的光圈半径A为150mm，视场角为1.0°，图像像素尺寸为512×512，像素尺寸为27μm×27μm，光学效率90％，量子效率30％和ADC增益为1e⁹。

设定曝光视场为2s，在曝光时间内光轴指向发生5次偏移，偏移量的范围为[-1,1]×10^-3rad。该条件下光学系统参数2下模拟的拖曳星图如图11(a)所示，光学系统参数3下模拟的拖曳星图如图11(b)所示。图11(a)中的轨迹亮度较低，多条轨迹淹没在噪声之中，信噪比较低，而图11(a)中的轨迹亮度过高，多条轨迹过曝，灰度细节丢失较多，两种情况均会对星图特征提取产生不利影响，故应优选光学系统参数1。

本发明保护范围不仅局限于实施例，实施例用于解释本发明，凡与本发明在相同原理和构思条件下的变更或修改均在本发明公开的保护范围之内。

Claims (9)

1.暗弱小天体光学导航特征成像模拟方法，其特征在于：包括如下步骤，

步骤一、给定暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图的光学系统参数；

给定暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图的光学系统参数，所述光学系统参数包括CCD相机的光圈半径A、视场角FOV、像素尺寸px/py、拍摄曝光时间T、光学效率η、量子效率Q和ADC增益K；

步骤二、给定暗弱小天体探测器拖曳星图的背景噪声参数，并根据给定的背景噪声参数生成噪声背景；

给定暗弱小天体探测器拖曳星图的背景噪声参数，所述背景噪声参数包括读数噪声偏置值λ_bias、暗电流噪声泊松参数λ_dark、成像背景亮度Sky_level，并根据给定的背景噪声参数生成噪声背景；

步骤三、根据目标小天体所在位置确定光轴指向，检索星表确定在该指向下能够进入视场范围的导航恒星列表，并计算所述导航星在CCD平面上的位置；

步骤四、给定PSF能量扩散半径，根据背景导航星及目标小天体的视星等调用高斯扩散公式得到各星点散焦后的能量；

步骤五、随机生成一系列光轴指向扰动，计算由该扰动产生的星图拖曳轨迹的能量，生成拖曳星图，即实现模拟暗弱小天体探测器任务远距离接近段长曝光下产生的拖曳星图。

2.如权利要求1所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：还包括步骤六，将步骤五得到暗弱小天体探测器拖曳星图用于验证轨迹图像星点质心提取的图像处理算法的有效性，进行图像处理算法的测试。

3.如权利要求1所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：还包括步骤七，改变光学系统参数，通过步骤一至步骤五模拟不同光学系统参数下的暗弱小天体探测器拖曳星图，进而选定满足实际深空探测任务的光学系统参数，实现实际深空探测任务。

4.如权利要求1、2或3所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：步骤一的具体实现方法为，

CCD相机的感应器件将入射光子转化为电子并将生成的电子贮存起来，随后，读入器件测量每个像素中的电荷量在电容器两端产生模拟电压，模数转换器将每个像素中的电压数字化，产生一个计数写入电脑磁盘中；最后将电脑磁盘中的计数转化为灰度，生成灰度图像；

每个像素点接收到的入射光子数决定该像素点处信号的大小；当光学系统入瞳处的光通量为F时，CCD处接受到的光子数P为：

其中，为平均波长，h为普朗克常数，c为光速，A为光圈半径，η为光学效率，；

入射光子数P与生成光电子C的比率为光学系统的量子效率Q：

C＝Q*P (2)

光电子数C与计数值N之间的关系由ADC的增益K所确定：

N＝K*C (3)

星载光学系统生成的图像为8位灰度图像，有256级灰度，计数N和灰度DN间的转换关系为线性关系：

DN＝a*N+b (4)

5.如权利要求4所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：步骤二的具体实现方法为，

总噪声包括读数噪声、暗电流噪声、成像背景亮度噪声：

S_噪＝I_background+I_read+I_dark (5)

由于深空中星际尘埃、太空碎片及其他物质的影响，拍摄的拖曳星图中存在大量杂散光，背景并不是完全的黑色；模拟星图的背景由用户输入参数背景亮度Sky_level和曝光时间T所决定，服从泊松分布，背景信号强度写为：

I_background～P(Sky_level*T) (6)

读数噪声是读数器件自身的电压所带来的噪声，所述读数噪声服从关于偏置值的高斯分布，其信号强度写为：

暗电流噪声是由CCD的结构内热造成的，其在单位时间内的噪声水平与CCD温度相关，属于散粒噪声的一种，服从泊松分布，其信号强度写为；

I_dark～P(λ_dark) (8)

根据给定的背景噪声参数生成噪声背景。

6.如权利要求5所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：步骤三的具体实现方法为，

定义光学相机的视轴方向为(α₀,δ₀)，则出现在光学相机视场内的导航目标所需满足的赤经赤纬条件如下：

其中，(α_i,δ_i)为导航恒星的赤经赤纬，FOV为视场角大小；

考虑天球坐标系和相机坐标系之间的坐标变换；由于导航星距离遥远，将相机坐标系的原点O和天球坐标系的原点O'视为重合，忽略坐标系平移所带来的误差,仅需要考虑旋转变换；导航恒星在天球坐标系下的坐标为(α_i,δ_i)，天球赤道坐标系O-UVW和相机坐标系O'-XYZ间的转换关系为：

[X Y Z]^T＝M[U V W]^T (10)

[U V W]^T＝[cosα_icosδ_i sinα_icosδ_i sinδ_i]^T (11)

所述O'-XYZ坐标系的Z轴为视轴正向，X轴和Y轴方向与CCD相机平面方向一致；

其中，j角与CCD相机安装方向有关，取j＝90°，θ为视轴与OW轴的夹角，O为O'X轴与OU轴的夹角，θ＝90°-δ₀，O＝90°+α₀；因此，M矩阵简化为：

最后，考虑相机坐标系和CCD成像平面坐标系之间的坐标变换；根据光学相机的成像模型，得到相机坐标系和CCD平面坐标系之间的几何关系：

上式中，(x,y)为导航目标在CCD成像平面投影的坐标，f为相机焦距，px,py为像素的高和宽；

同时，根据CCD相机视场角大小、传感器尺寸和焦距间的关系知：

f＝2pxNx^-1tan(FOV/2)＝2pyNy^-1tan(FOV/2) (14)

其中，Nx,Ny是CCD平面尺寸，FOV为视场角大小；

最终得导航目标在CCD成像平面投影的坐标位置的表达式：

所述导航目标在CCD成像平面投影的坐标位置即为导航星在CCD平面上的位置。

7.如权利要求6所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：步骤四的具体实现方法为，

在天文观测中，视星等与单色光流间的关系为：

m_v＝-2.5log₁₀F_v+c (16)

其中，m_v为信源视星等，F_v为信源单色光流，c为常数；

由公式(17)得出，信源的成像强度主要由信源的星等决定，信源视亮度与信号强度之间的关系为：

其中，I为某像素点上的光能，I_m为0等星所对应的光能，mag为某像素点对应的星等；

由于深空导航目标为点光源，需要考虑在CCD焦平面上造成的点扩散问题，在输出的数字灰度图上实际表现为包含若干像素点的光斑；采用Gauss点扩散函数能够分别计算出目标及导航星在CCD面的投影位置周围像素收集的能量：

上式中，(m,n)为所计算的像素点坐标，I(m,n)为所计算的像素点收集的光能，I₀为目标点光源的光能，σ_PSF为Gauss弥散半径，(x₀,y₀)为目标点光源在CCD面的投影坐标；

根据公式(18)能够分别计算各星点散焦后的能量。

8.如权利要求7所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：步骤五的具体实现方法为：

由于光轴扰动，在单次曝光时间T内相机光轴发生k次变化，每次扰动后的天球赤道坐标系和相机坐标系之间转换矩阵记为M_k，扰动后导航目标投影的位置坐标：

上式中，Nx,Ny是CCD平面尺寸，FOV为视场角大小，M_k1,M_k2,M_k3为M_k各行，X_i为导航恒星在天球直角坐标系下的位置；

由于两次扰动之间的时间间隔较短，故将两次光轴方向变化间的运动视为是匀速的，在任意时刻t有：

M_t＝M_k·(I+ω·Δt) (20)

其中，M_t为任意时刻t的转换矩阵，ω为光轴的方向变换角速度，Δt为时间步长；

在任意时刻t，将M_t作为天球赤道坐标系和相机坐标系间的状态转移矩阵M，按照步骤三和步骤四计算任意时刻t星点的位置和能量；将不同时刻下的星点能量叠加即得到暗弱小天体探测器拖曳星图。

9.如权利要求8所述的暗弱小天体光学导航特征成像模拟方法，其特征在于：步骤六具体实现方法为，

步骤6.1：对长曝光下产生的拖曳星图进行去噪预处理；

考虑到图像噪声的非高斯性及低信噪比性，采用阈值小波变化的方法进行去噪；小波去噪方法的思路是：进行小波变换后，信号主要集中在低频，而噪声集中在高频上，故给定适当的阈值除去噪声；

采用软阈值小波方法进行去噪，软阈值是指若小波变换系数比设置的阈值小则将该系数置零，若小波变换系数比该阈值大则将该系数向阈值方向收敛，小波变换系数表达形式为：

其中，ω_i,j为小波变换系数，λ为设定的阈值；

步骤6.2：通过连通域判定轨迹所在的兴趣区域(ROI)，通过兴趣轨迹选取准则选取兴趣轨迹，采用逐像素细化算法(One Pass Thinning Algorithm)对兴趣轨迹进行细化并提取出轨迹骨架；

长曝光拖曳图像模型表示式为：

I_Blur＝K·I_static+E (22)

其中，I_Blur是长曝光的拖曳星图，I_static是需要被复原的静态星图，E为噪声项，K为是每一像素的模糊转移矩阵；当星图中各像素发生的模糊是均匀时，则式(2)变换为卷积的形式，如式(3)所示：

I_Blur＝k*I_static+e (23)

其中，k被称为模糊核，与静态星图I_static的每一像素发生卷积运算使图像产生模糊，其物理含义是反映每一像素点的运动特性；

由于长曝光拖曳星图的拖曳是点状目标的运动模糊，故星点拖曳轨迹本身能够反映地曝光时间内的运动，而模糊核的本质就是像素点的运动特性，故采用星点拖曳轨迹来估计拖曳图像的模糊核；

通过连通域判定轨迹所在的兴趣区域(ROI)，通过兴趣轨迹选取准则选取兴趣轨迹；

所述兴趣轨迹选取准则选为亮度兴趣轨迹选取准则，即选取亮度最亮的轨迹为兴趣轨迹，兴趣轨迹的选取准则为：

max average(G_{ROI_i}) (24)

其中，G_{ROI_i}是第i个兴趣区的灰度矩阵；

OPTA算法是一种基于形态学的单循环并行模板匹配算法，其优点是运算简单，计算速度很快，且能够较好地保留图像的连通性和形状特征；采用OPTA算法对兴趣轨迹进行细化并提取出轨迹骨架，具体实现方法如下：按从左到右，从上到下的方式，从最左上角的像素点开始抽取出每个灰度值不为0的像素点的邻域，将该邻域与消除模板中的任意一个消除模板匹配，若不能匹配则保留像素点的值，否则要再将该邻域与保留模板匹配，若能匹配则仍然保留像素点的值，若不能匹配则删去该像素点，即将像素点的灰度值置零；对所有像素点进行上述操作，保留下的像素值不为0的像素点即为提取出的轨迹骨架；

步骤6.3：以步骤6.2中提取出的轨迹骨架为模糊核的估计初值，通过稀疏先验信息构建图像复原算法的代价函数，对长曝光拖曳星图进行盲复原，得到长曝光拖曳星图复原图像；

复原图像问题需估计出合适的模糊核及复原图像，使得图像复原的代价函数尽量最小，代价函数形式为：

然而，复原图像问题本质上是NP(non-deterministic polynomial)问题，即由于信息的匮乏，会导致结果不唯一甚至病态；为解决星图模糊复原的问题，需要根据导航星图的实际特点引入相关先验信息建立正则惩罚项；

深空导航星图的重要特征是星点能量集中在极少数像素点中，通过步骤一去噪声预处理后，绝大部分像素点为黑色背景，故利用图像的灰度稀疏作为先验信息；同理，由于有值像素点较少，故图像的边缘像素点也就更少，所以图像的灰度梯度也是稀疏的；稀疏先验信息通过L₀范数构建，即

最终采用的星图复原代价函数的形式为：

其中，项是为保证估计出的模糊核的简单性；

将步骤6.2中提取出的轨迹骨架为模糊核的估计初值，以代价函数最小为指标估计复原图像；由估计出的复原图像更新模糊核；如此循环迭代直至达到给定的循环的终止条件，此时估计出的复原图像即为最终得到的长曝光拖曳星图的复原图像；

步骤6.4：对步骤6.3得到的复原图像进行连通域判定，并提取复原图像中各星点的灰度光心，完成长曝光拖曳星图特征提取；

对步骤6.3得到的长曝光拖曳星图复原图进行连通域判定，去除连通域小于预设像素数量的目标；

步骤6.4中去除连通域小于预设像素数量的目标选去除连通域小于2个像素的目标；

然后，根据导航星覆盖像素的像素值确定复原星图导航星光心所在的位置，即实现提取复原图像中各星点的灰度光心，完成长曝光拖曳星图特征提取，并验证轨迹图像星点质心提取的图像处理算法的有效性，即进行图像处理算法的测试；

步骤6.4中计算灰度光心的方法采用灰度加权算法，灰度加权算法的公式为：