引言

近年来，地球轨道上分布的废弃卫星、火箭残骸、航天器碎片等人造物体大量增加，严重影响了航天器在轨运行的安全。地基光电探测系统在中高轨的空间目标观测中占据绝对优势^[1-4]。而杂散光是影响地基光电探测系统成像质量的重要因素，若不能得到有效抑制，轻则降低成像对比度、目标信噪比，严重时甚至会将探测目标淹没^[5-6]。杂散光是光学系统中辐射于探测器表面的非成像光线，以及经非正常路径传播到达探测器表面的成像光线^[7-9]。杂散光按照其来源可分为三类：外部杂散光、内部杂散光和成像杂散光。针对可见光探测系统，由系统外辐射源直接照射或经过系统内构件多次反射、散射到达探测器而造成的外部杂散光起主要作用，这类辐射源包括太阳光、月光、地球表面和大气的散射、漫射光等。

当前，针对外部杂散光抑制的研究可分为以下两类^[10]：一是从光学系统的角度，通过对杂散光来源进行分析并改进光学系统的设计，从根本上抑制杂散光；二是从图像角度，根据图像中杂散光的特征，利用退化图像恢复出清晰图像。其中，针对后者的研究主要应用于抑制遥感图像中的杂散光噪声。针对我国FY-2气象卫星扫描辐射计可见光通道图像存在的明显杂散光，文献[11]在分析了杂散光产生机理的基础上，进行分析建模，通过学习各通道地球圆盘外的杂散光特征并将其推广到圆盘内，从而估计出全视场范围的杂散光分布，最终实现杂散光的去除；文献[12]根据光机建模确定了杂散光分布并建立漏光模板，模拟了杂散光的叠加过程，并通过反卷积实现了卫星图像中杂散光的消除；此外，文献[13]针对SJ-9A卫星光学遥感图像中的条纹杂散光噪声，通过统计杂散光在图像上的分布特征建立特征函数，并采用分块自适应算法实现了条纹噪声的去除。然而，上述方法存在一定的不足：1） 通过改进光学系统设计抑制杂散光的方法一般在设备设计阶段完成，对于已投入使用的光学设备，重新对其结构进行设计、改进难以实现；2） 由于地基光电探测设备布置于地面，拍摄到的图像与遥感图像有所区别。同时，观测过程中受气象条件以及大气影响，使得观测图像中杂散光噪声的分布不断变化，现有图像处理方法往往难以适用。

本文针对地基可见光观测图像中由于大气干扰而产生的杂散光噪声问题，提出了一种基于图像亮度先验的杂散光噪声去除方法。杂散光去除后图像整体亮度均匀，同时目标信杂比（signal-to-clutter ratio, SCR）得到大幅提升，相较于现有星图预处理算法，取得了更好的实验效果。

1   星图亮度先验

本文通过对大量受杂散光干扰的实测数据进行统计分析，对杂散光噪声的空间分布等图像特征展开了相应的研究，发现实际拍摄星图中的杂散光噪声呈现以下复杂特性：1） 降低区域内目标信噪比，在噪声强度不大时，目标特性不明显，清晰度降低；在噪声集中的区域，目标与噪声混杂、难以分辨；2） 在图像中呈现不均匀弥散特性，杂散光在图像中往往集中于一个区域，其强度沿中心向外逐渐减弱；3） 星图中像素的亮度值随着杂散光噪声强度的变化而急剧变化。

以某张地基可见光观测设备拍摄的实际星图为例，见图1。

图  1  星图特性及亮度先验示意图

Figure  1.  Schematic diagram of star maps characteristics and intensity prior

下载: 全尺寸图片 幻灯片

图1（a）为原始图像，由于大气干扰的原因，图中存在明显的杂散光噪声，严重降低了区域内目标的信噪比。图1（b）和图1（c）展示了对图1（a）所示图像按行和列统计像素灰度均值的结果，可以看到，杂散光在星图中的空间分布并不均匀，呈现出由中心向外弥散的特性。为进一步说明星图亮度与杂散光噪声强度的关系，在图中选取噪声干扰程度不同的3个区域，分析其亮度值变化如下：图1（d）展示了无杂散光噪声的清晰区域，其亮度值较低；图1（e）展示了在有一定杂散光噪声干扰的区域，由于杂散光的存在，亮度值显著增加；图1（f）为杂散光严重干扰的区域，图像中的暗弱目标几乎被杂散光噪声覆盖，亮度值较图1（e）进一步增加。通过以上观察可以得出，杂散光噪声的强度与星图的亮度呈现出有规则的变化。下面通过分析星图成像过程来进一步验证这一规律。

在无杂散光干扰的理想情况下，空间目标反射来自入射光源的能量，到达成像设备时几乎没有能量损失，成像系统接收到的目标反射的入射能量聚焦到图像平面上，在没有杂散光影响时，空间目标成像清晰，亮度分布均匀。而在有杂散光干扰的情况下，如图2所示，由于地球大气的存在，目标的反射光在穿过大气的过程中，会受到空气中的各种分子以及悬浮粒子吸收、散射的影响：一方面，由于目标反射光的衰减，其在图像中的亮度值降低；另一方面，由于包括月光等在内的大气光散射后增强了天空背景的亮度，且大气散射对成像设备的影响是可以不断累加的，在大多数情况下大气散射模型对成像的影响占主导作用。因此，实际场景中大气散射严重的区域对应于图像中杂散光噪声集中的区域，呈现出亮度高的特点。

图  2  地基观测设备的杂散光形成过程

Figure  2.  Stray light formation process of ground-based observation equipment

下载: 全尺寸图片 幻灯片

根据以上分析，杂散光噪声强度与大气深度、星图亮度呈正相关，即：

式中：$s(x)$表示图像中像素点$x$处的杂散光噪声强度；$d(x)$表示该点处大气的深度，其值越大，表示大气干扰程度越高；$i(x)$表示该点的亮度值。将（1）式所描述的统计规律称为亮度先验。值得指出的是，（1）式只是一个观察到的直观结果，并非定量的关系表达式，下一节将具体分析并建立描述$d(x)$与$i(x)$之间关系的数学表达式。

通过以上对星图特性及杂散光形成过程的分析，根据McCartney提出的大气散射模型^[14]，建立星图在杂散光干扰下的退化模型如下：

式中：$x$表示星图上像素点的位置坐标；$I(x)$表示实际拍摄到的受杂散光噪声干扰的星图，也即场景反射光经过衰减后到达观测设备的光线强度；$t(x)$表示大气透射率，反映了光线穿透大气的能力，其值越大，表明从场景表面的反射光穿透大气到达观测设备视场的数量越多；$J(x)$表示待恢复的无杂散光干扰图像，也即空间目标直接反射光强；$A$为大气光值，表示大气干扰最严重区域的光线强度；$\;\beta $是散射系数，表示入射光线的散射程度，在均匀介质中，可以将其看作常数；$d(x)$表示大气深度，值越大，表明大气散射干扰的程度越高。式中只有$I(x)$为已知，杂散光去除的实质就是求出$t(x)$和$A$的值，代入（2）式中解得$J(x)$。

2   基于亮度先验的星图杂散光噪声去除算法

2.1   基于图像亮度的大气深度恢复

根据所描述的大气散射模型可知，恢复无杂散光干扰的清晰图像需要获得大气光值和大气透射率，而大气深度对这两者的估计至关重要。根据文献[15]，假设大气为均匀介质，而将散射系数$\beta $视为常数，此时，若能获得大气深度$d(x)$，便可估计出大气透射率$t(x)$。同时，根据（2）～（3）式可知，当$d(x) \to \infty $时，有：

即：大气深度趋于无穷大的像素点的亮度可以代表大气光值$A$。

因此，去除杂散光的任务进一步转化为大气深度的恢复，根据亮度先验，可以通过亮度值从单幅图像中直接恢复出大气的深度信息，从而实现无杂散光干扰清晰图像的求解。

2.1.1   线性模型建立

根据观察到的亮度先验，亮度与大气深度呈正相关，因此，可建立线性模型：

式中：$x$是图像中的像素点；$d(x)$是该点的大气深度；$i(x)$是杂散光干扰星图中对应的亮度值；${\theta _0}$、${\theta _1}$为待求解的线性系数；$\varepsilon (x)$是随机变量，表示模型的随机误差，假设$\varepsilon (x) \sim N(0,{\sigma ^2})$，由正态分布的性质可知：

2.1.2   数据集构建

由于图像的亮度值$i$是已知的，因此，使用有监督学习的方法确定模型中的变量${\theta _0}$、${\theta _1}$和$\sigma $。参考文献[16]的方法，收集大量的清晰图像，并根据星图的退化模型，使用清晰图像产生合成大气深度图和相应的杂散光干扰图像，以获得足够的训练样本。

2.1.3   学习策略

首先建立关于大气深度的联合概率密度函数如（7）式：

式中：$n$表示训练的杂散光干扰星图中的像素总数；$d\left( {{x_n}} \right)$表示第$n$个像素点的大气深度；$L$表示似然性。假设各场景点的随机误差相互独立，即$p\left( {{\varepsilon _1}, \cdots ,{\varepsilon _{\rm{n}}}} \right) = \prod\limits_{k = 1, \cdots ,{\rm{n}}} p \left( {{\varepsilon _k}} \right)$，则可得似然函数：

将（6）式代入上式，可得：

式中${g_k}$表示像素点$k$处真实的大气深度。似然函数取得最大值表示相应的模型参数能够使观测数据出现的概率最大，因此，问题转化为寻找最优的参数${\theta _0}$、${\theta _1}$和$\sigma $使得似然函数L最大化。为了便于计算，对（9）式两端同取对数，然后对$\ln L$做极大似然估计，可得：

为了求解该问题，取$\ln L$关于$\sigma $的偏导数，并令其为0，有：

根据（11）式，变量${\sigma ^2}$的最大似然估计为

使用梯度下降法对线性系数${\theta _0}$、${\theta _1}$进行求解，分别取$\ln L$对${\theta _0}$、${\theta _1}$的偏导数，得到：

通过大量训练样本学习该线性模型，得到的最佳系数分别为${\theta _0} = 0.121\;779$，${\theta _1} = 0.959\;710$,$\sigma = $$ 0.041\;337$。确定了相关系数，即可实现图像亮度与大气深度之间的线性变换。

2.1.4   深度图恢复及优化

根据大气深度与星图亮度之间建立的关系模型，以及学习出的最优系数，可以根据（5）式实现亮度$i$到大气深度$d$的直接转换，通过这种方式直接获取的深度图如图3（a）所示。可以看出，该模型会将星图中亮度值较高的恒星当作杂散光噪声，导致估计出的大气深度不准确。虽然对原始深度图进行最小值滤波可以解决该问题，但如图3（b）所示，处理后的深度图中会出现块效应。因此，为了抑制产生的块效应以及得到更精确的深度图，本文使用导向滤波^[17]对求得的深度图进行优化。

图  3  优化后的深度图（局部放大）

Figure  3.  Optimized depth map (partial zoom)

下载: 全尺寸图片 幻灯片

导向滤波可以平滑图像细节和保持图像的边缘信息，本文采用原始星图I作为引导图，根据导向滤波局部线性模型假设，输出的优化后的深度图q可表示为

式中：${\omega _k}$为引导图I中像素点$k$为中心的邻域，其窗口半径为r；$j$表示邻域${\omega _k}$中的像素点；系数${a_k}$、${b_k}$在邻域${w_k}$中为常数。根据（15）式，在局部区域，输出图像可以捕获与引导图相似的细节信息，其捕获细节信息的能力与半径r的大小密切相关，经大量实验测试，本文取$r = 15$。

导向滤波通过最小化输出图像q与输入图像p之间的差异来寻求最优的系数${a_k}$、${b_k}$，即在邻域${w_k}$中，最小化（16）式所描述的目标函数：

式中$\lambda $是正则化参数，用来防止${a_k}$过大，其对优化结果不敏感，经大量实验测试，本文取$\lambda = {10^{ - 3}}$。

本文通过（5）式对深度图进行粗估计，然后使用最小值滤波去除原始深度图中估计错误的亮目标，最后利用导向滤波抑制最小值滤波产生的块效应并对深度图进行精确化，从而有效改善杂散光的抑制效果。导向滤波后的深度图如图3（c）所示。

2.2   大气光值估计

根据上一节的分析，大气光值可以根据输入图像恢复的深度图进行估计。由于深度图中像素点的亮度与输入图像中杂散光噪声的强度呈正相关，因此可以通过深度图获取杂散光噪声最强的区域，以此来估计大气光值。具体而言，选取深度图中亮度最高的前0.1%像素点，并在这些点对应的输入图像中，选择亮度最高的像素点作为大气光值。根据上述方法，只要能够准确估计出深度图，大气光值A便可落到图像中杂散光噪声干扰区域。如图4中红色点表示的区域，即为根据上述方法得到的像素点对应到原始图像中的像素位置。

图  4  大气光值在原图中的位置（局部放大）

Figure  4.  Position of atmospheric light value in original image (partial zoom)

下载: 全尺寸图片 幻灯片

2.3   清晰图恢复

在已知大气深度d和大气光值A后，通过（3）式即可估计出大气透射率t，并通过（2）式恢复无杂散光干扰的图像$J(x)$，即：

经大量实验测试，当$t(x)$的值在0.3～1.0之间时，对透射率图中的噪声具有较好的抑制效果。本文用于恢复无杂散光清晰图像的最终函数可表示为

$J(x)$即为求解的无杂散光噪声干扰星图。

3   实验结果与分析

为验证本文方法对实际图像预处理的性能，对大视场地基光学望远镜拍摄的星图进行杂散光去除实验，每幅图像的尺寸为4096×4096像素。选取5个受不同程度杂散光影响的区域，通过统计算法处理前后星图序列中具有针对性的空间目标区域的变化情况来评价本文方法的有效性。

3.1   算法评价指标

由于本文算法作为后续空间目标检测的预处理算法，最终目的是为了提高空间目标检测、跟踪的精度。因此，在评价算法时，除算法处理前后的主观视觉效果外，主要关注目标区域的背景抑制情况及信杂比变化。

采用空间目标检测领域常用的SCR增益作为目标增强效果的客观评价指标。具体而言，选择序列图像中5个观测得到的空间目标，统计序列中每帧图像处理前后，以各目标为中心的15×15像素窗口区域内的信杂比增益。信杂比增益${R_{{\rm{SCRgain}}}}$的定义如（19）式：

式中：${R_{{\rm{SCRout}}}}$表示算法处理后图像的信杂比；${R_{{\rm{SCRin}}}}$表示原始图像的信杂比。${R_{{\rm{SCRout}}}}$和${R_{{\rm{SCRin}}}}$均可以通过（20）式计算：

式中：${R_{{\rm{SC}}}}$表示信杂比；${m_{\rm{t}}}$为目标平均灰度；${m_{\rm{b}}}$为背景平均灰度；${\sigma _{\rm{b}}}$为背景标准差。

采用背景抑制因子（background suppression factor, BSF）作为背景抑制效果的客观评价指标，为了更好地反映星图的背景变化情况，选择序列图像中5个观测得到的空间目标，统计序列中每帧图像处理前后，以各目标为中心的205×205像素窗口区域内的背景抑制因子。背景抑制因子${F_{{\rm{BS}}}}$的定义如（21）式：

式中：${\sigma _{{\rm{in}}}}$表示原始图像的标准差；${\sigma _{{\rm{out}}}}$表示算法处理后图像的标准差。

3.2   算法实现与对比

利用所提出的算法对星图序列进行处理，图5（a）给出了包含空间目标的原始图像局部放大图，图5（b）给出了本文算法处理后的结果，图中方框标示的是空间目标。图6更直观地展示了本文算法的有效性，给出了星图序列处理前后分别以5个观测目标为中心、15×15像素区域内SCR的变化情况。可以看出，本文算法提高了整个序列中各目标的SCR。值得指出的是：1） 同一图像序列中的不同目标处于图像的不同区域，受到的杂散光干扰差异极大，因此目标的信杂比有所不同；2） 针对单组图像序列中的同一目标，本文算法增强后的目标SCR提升效果均有所不同，这是由原图目标运动过程中复杂的背景变化造成的。当空间目标由杂散光干扰较大的区域运动到杂散光干扰较小的区域时，由于背景噪声的大幅下降，其局部区域内的信杂比会有所提升，反之则有所降低。

图  5  本文算法处理前后效果对比

Figure  5.  Comparison of effects before and after algorithm processing

下载: 全尺寸图片 幻灯片

图  6  本文算法处理前后星图序列中5个目标的信杂比变化

Figure  6.  Changes in SCR of 5 targets in star maps sequence before and after algorithm processing

下载: 全尺寸图片 幻灯片

表1进一步说明了所提出算法的目标增强和背景抑制效果，对本文算法处理星图序列前后的信杂比、信杂比增益以及BSF进行统计；同时，与其他5种典型算法的处理结果进行对比，包括max-mean^[18]、TDLMS^[19]、ILCM^[20]、Top-hat和median-filter。其中，max-mean、TDLMS、top-hat是经典的目标增强算法，ILCM是基于人类视觉系统的算法，median-filter是经典的去噪算法。表1中SCR和BSF的值是星图序列中每帧结果的平均值，其中，信杂比增益的定义如（19）式所示，BSF的定义如（21）式所示。从表1可以看出，对含有杂散光噪声的星图而言，在信杂比提升方面，经典的max-mean、TDLMS算法以及Top-hat算法效果较差，ILCM和median-filter算法对于目标4和目标5有一定效果，而本文算法效果最为突出；在背景抑制方面，median-filter算法和本文算法效果较好。总体而言，本文算法不仅大幅度提升了目标的信杂比，也有效地抑制了背景噪声；同时，对于不同背景下的目标具有较高的鲁棒性，相比于传统去噪算法及空间目标检测预处理算法，具有较为明显的优势。

表  1  不同算法去除星图序列杂散光的定量评价结果

Table  1.  Quantitative evaluation results of different algorithms for removing stray light from star maps sequence

目标 序号	原图 SCR		算法 名称	处理 后SCR	SCR 增益	BSF
1	max	9.267 8	max-mean	4.279 5	−4.853 3	1.423 7
			TDLMS	6.908 0	−0.731 3	1.259 4
	mean	7.513 3	ILCM	8.265 7	0.832 5	0.683 4
			Top-hat	6.518 9	−1.298 6	1.165 3
	min	6.083 8	median-filter	9.959 7	2.518 4	1.721 2
			Ours	17.736 2	7.391 0	2.126 1
2	max	8.840 5	max-mean	3.420 0	−6.191 1	1.739 0
			TDLMS	6.703 9	−0.374 3	1.194 0
	mean	6.992 1	ILCM	6.642 8	−0.465 8	0.651 1
			Top-hat	7.027 3	0.039 8	1.009 8
	min	5.682 0	median-filter	6.909 5	−0.098 5	2.896 6
			Ours	29.674 9	12.005 3	2.922 6
3	max	6.242 8	max-mean	3.479 5	−1.607 8	1.571 9
			TDLMS	3.810 0	−1.003 0	1.848 4
	mean	4.227 3	ILCM	6.301 8	3.561 7	1.113 3
			Top-hat	4.050 5	−0.380 7	1.099 5
	min	3.118 9	median-filter	4.094 7	−0.495 0	1.955 4
			Ours	24.548 9	15.046 2	2.142 0
4	max	2.924 9	max-mean	2.860 4	2.933 7	1.375 3
			TDLMS	1.981 5	−0.300 0	1.650 8
	mean	2.051 7	ILCM	5.066 8	7.846 4	1.355 5
			Top-hat	1.800 4	−1.125 9	1.069 6
	min	1.341 2	median-filter	3.698 3	4.959 2	2.411 0
			Ours	8.350 7	11.811 0	1.920 1
5	max	3.890 9	max-mean	2.820 8	1.439 8	1.685 1
			TDLMS	2.232 3	−0.216 5	1.366 8
	mean	2.297 4	ILCM	5.157 0	5.970 6	1.140 0
			Top-hat	2.111 8	−0.684 3	1.038 7
	min	1.248 1	median-filter	3.532 8	3.195 1	2.609 2
			Ours	16.570 0	15.528 6	3.268 0

下载: 导出CSV 

| 显示表格

4   结论

本文提出了一种基于亮度先验的星图杂散光噪声去除方法。在对大量地基可见光观测设备获取的实际星图进行统计分析的基础上，建立了星图在杂散光干扰下的退化模型。根据亮度先验，从单幅星图中可直接恢复出杂散光强度信息，从而实现无杂散光干扰清晰图像的求解。与现有算法相比，对于受不同程度杂散光干扰的目标，本文方法在背景抑制和目标信杂比提升上均获得了更好的实验效果。其中，针对序列星图中信杂比为2.05以上的空间目标，处理后SCR增益和BSF分别能够达到7.39和1.92以上。本文方法利用杂散光在图像中的分布特征来解决星图杂散光干扰问题，相比于改进光学系统设计的方法成本更低、更易于实现，且可作为地基光电望远镜设备部署后的补救措施，具有广泛的应用价值。