CN104504674A - 空间碎片星点提取与定位方法 - Google Patents

Abstract

本发明公开了一种空间碎片星点提取与定位方法，用于解决现有方法对相邻空间碎片目标分割困难的技术问题。技术方案是通过MSER算法在候选空间碎片区域内过构建数据结构树，提取出该区域内空间碎片像素集合并确定空间碎片星点个数。在MSER算法得到的像素点树形结构中，通过检测该候选空间碎片区域中最大的MSER，为后续拟合GMM模型生成训练样本。在基于GMM的定位过程中，通过EM算法迭代求解混合模型中高斯成分的均值与协方差，对应的均值就是各个空间碎片星点光斑的质心。该方法结合MSER和GMM模型的特点，通过GMM模型对近距离弥散光斑质心进行共同求解，解决了相邻空间碎片目标分割困难的技术问题。

Description

空间碎片星点提取与定位方法

技术领域

本发明涉及一种空间碎片星点提取与定位方法。

背景技术

文献“基于能量累加的空间目标星像质心定位，《光学精密工程》，2011，第12期，p2043-3048”提出了基于能量累积的插值求解质心方法。该方法使用插值方式降低离散采样点和实际感光区域光线照度的不一致性；利用弥散光斑能量累加自适应地确定空间碎片目标的有效窗口。由于这种自适应确定空间碎片目标有效窗口的方法仅适用于单个独立的空间碎片进行定位，当多个空间碎片由于距离较近，而使得图像上的星点光斑相互黏连或重叠时就会难以分割提取和定位。

发明内容

为了克服现有方法对相邻空间碎片目标分割困难的不足，本发明提供一种空间碎片星点提取与定位方法。该方法通过MSER算法在候选空间碎片区域内，通过构建数据结构树，提取出该区域内空间碎片像素集合并确定空间碎片星点个数，为利用高斯混合模型对星点光斑建模和利用最大期望EM算法求解是的初始化提供了基础。在MSER算法得到的像素点树形结构中，通过检测该候选空间碎片区域中最大的MSER可以为生成后续拟合GMM模型生成训练样本。在基于GMM的定位过程中，通过EM算法迭代求解混合模型中高斯成分的均值与协方差，对应的均值就是各个空间碎片星点光斑的质心。该方法结合MSER和GMM模型的特点，通过GMM模型对近距离弥散光斑质心进行共同求解，避免了在空间碎片目标的弥散光斑靠近、黏连甚至部分重叠时分割的困难。

本发明解决其技术问题所采用的技术方案是：一种空间碎片星点提取与定位方法，其特点是包括以下步骤：

步骤一、对于星空图像I，首先对针对星空图像中每个像素统计图像的灰度直方图，根据像素值按升序对图像中的像素点进行排序，去除像素灰度值最小的λ₁％的像素点，以减弱暗噪声对背景的影响，去除像素灰度值最高的λ₂％的像素点，以减弱强噪声和恒星对背景估计的影响，其中λ₁,λ₂为设定的背景阈值参数。利用剩余的(100-λ₁-λ₂)％像素点，求得背景的灰度均值μ和灰度标准差σ，使用阈值μ+κσ对图像 进行二值化处理，κ为需要设置的分割系数。二值化处理公式为：

$B(x,y)=\left\{\begin{array}{cc}1& I(x,y)>\mathrm{\μ}+\mathrm{\κ\σ}\\ 0& I(x,y)\≤\mathrm{\μ}+\mathrm{\κ\σ}\end{array}\right.---\left(1\right)$

式中，I(x,y)表示在坐标为(x,y)处的图像I中像素的灰度值，B(x,y)就是(x,y)处像素点的二值化结果，把二值化的后的图像记为B。这一步为星图的二值化粗分割，因此每个B(x,y)＝1的连通区域中可能包含多于一个的星点。

步骤二、对空间碎片区域进行粗分割。根据上一步图像二值化的结果，当B(x,y)＝1时表示该像素为空间碎片星点区域，当B(x,y)＝0时表示该像素为非有效空间碎片星点区域，即背景区域。对B中B(x,y)＝1的点进行广度优先搜索，提取到所有满足B(x,y)＝1的点的连通域，所有连通域集合为其中一个连通域表示为Ω。Ω为一个空间碎片星点区域，包含了一系列二值化结果为1的像素点坐标，并且他们相互连通，而且对于任一个Ω中的坐标为(x,y)的点，它在原始图像中的灰度值为I(x,y)。

步骤三、在提取空间碎片区域后，需要对每一个连通域Ω通过MSER算法进行局部极大灰度值和最大稳定极值区域检测。MSER检测是通过不同阈值对空间碎片区域进行二值化处理，根据不同的结果对稳定极大灰度值点和区域进行检测。对于一个星点区域Ω，当令二值化阈值为t时，公式(2)中的Ω_>t表示在空间碎片区域Ω中，像素灰度值大于t的像素点集合，当二值化阈值t不断变大的时候集合Ω_>t的元素个数会越来越少，同理t变小的时候集合Ω_>t的元素个数会变多。在MSER检测过程中通过集合中元素个数的变化速率来表示一个区域Ω_>t是否是一个稳定的空间碎片弥散光斑。公式(3)表示当分割阈值变化大小为d时，光斑变化的稳定度计算方法。这里|Ω_>t|，|Ω_>t+d|和|Ω_>t|分布表示集合Ω_>t，Ω_>t+d和Ω_>t的元素个数，d为需要设置参数，st为光斑的稳定度。

Ω_>t＝{(x,y)|I(x,y)>t,(x,y)∈Ω}   (2)

$\mathrm{st}=\frac{\left|{\mathrm{\Ω}}_{>t-d}\right|-\left|{\mathrm{\Ω}}_{>t+d}\right|}{\left|{\mathrm{\Ω}}_t\right|}---\left(3\right)$

利用MSER算法在一个空间碎片星点连通区域Ω中分割提取星点的步骤为：

1.首先对于图像中某个空间碎片星点连通区域Ω内的每个像素点灰度值按照大小排序，依次初始化他们为未访问状态再根据各个像素点的灰度值大小和他们空间位 置关系，按照步骤2、步骤3的描述把他们组织成一颗树形结构，这个树的各个节点都代表了图像中的一个像素点。

2.按照像素灰度值从大到小遍历Ω中所有的像素，针对某个(x,y)处的像素点遍历它周围相邻的8个点。假设(x′,y′)为一个与(x,y)相邻的位置，那么如果(x′,y′)处像素点未被遍历，则直接遍历下一个相邻点，如果已经被遍历那么：

a)当Root(x,y)＝Root(x′,y′)相同时，则不进行任何处理，遍历(x,y)处像素的下一个相邻点，因为(x,y)处像素点与(x′,y′)处像素点已经处于同一个子树中。这里的Root(x,y)表示(x,y)处像素点在当前树结构中对应的根像素点的位置。

b)当I(Root(x,y))＝I(Root(x′,y′))时，则把Root(x′,y′)处像素点作为子节点添加到Root(x,y)处像素点下。

c)当I(Root(x,y))＜I(Root(x′,y′))时，则把Root(x′,y′)处像素点添加为Root(x,y)处像素点的子节点。

3.标记当前(x,y)处像素点为已被遍历像素点，并重新回到步骤2直到完成对所有像素点的遍历。

4.当所有的像素点按照上面步骤遍历完毕，依照像素点灰度值大小关系的树形结构建立完成，该树中叶节点个数就等于该空间碎片区域Ω中局部极大灰度值像素点个数，也对应了区域内星点光斑数目，这个数目记为K，极值点的个数、位置将用于以下基于GMM模型和EM迭代算法的星点定位的初始化。

为了使用EM算法拟合GMM模型提供样本，需要在空间碎片区域内检测最大稳定极值区域。对于空间碎片光斑区域Ω内，每当取到一个二值化分割阈值t，总是对空间碎片区域Ω进行分割得到一个子光斑为Ω_>t，当该光斑的稳定度st小于t_mser时，即满足公式(4)时，该子光斑就是为Ω的一个MSER，这里t_mser为判断一个子光斑是否为MSER的阈值。当调整t，就得到一系列MSER，这一系列的MSER呈现树形的包含现象，即t较小时分割得到的MSER包含t较大时分割得到的MSER，而最大的MSER包含所有其他的MSER。对区域Ω提取不同t情况下的MSER之后得到MSER集合其中包含像素点最多的MSER记为C。

$\frac{\left|{\mathrm{\&Omega;}}_{>t-d}\right|-\left|{\mathrm{\&Omega;}}_{t+d}\right|}{\left|{\mathrm{\&Omega;}}_{>}t\right|}<t_{\mathrm{mser}}---\left(4\right)$

步骤四、考虑单个空间碎片弥散光斑可以被合理假设为一个高斯函数形状的的点状光斑，对于包含多个弥散光斑的空间碎片区域，利用混合高斯模型(GMM)对空间碎片区域内的多个弥散光斑进行建模。GMM在质心定位的时候需要多个像素点样本。

C为空间碎片区域Ω的MSER集合中包含最多像素点的MSER，对其范围内任一位置(x,y)处的像素点进行重采样，得到个样本点，这个样本点的集合记为S_(x,y)，为向上取整符号。那么C中所有的样本点集合记为S_C，其中n＝1,...,L表示像素点样本的序号，L表示样本点个数，即集合S_C的中的成员数目，m_n表示第n个样本像素点的位置。

$S_C=\left\{m_n\right\}=\underset{(x,y)\&Element;c}{\&cup;}{S}_{(x,y)},n=1,...,L---\left(5\right)$

获得了样本点集合S_C和星点光斑数目K，利用GMM模型和相应的EM算法对空间碎片光斑区域Ω内K个光斑同时进行定位。具体的模型描述和定位算法如下：

1.定义用于拟合空间碎片区域星点光斑的混合高斯模型。由于空间碎片区域是有多个弥散光斑组成，首先将对于空间碎片区域中的任一位置点变量(x,y)记为随机变量m＝(x,y)，它服从公式(6)描述的高斯混合分布，其中π_k、μ_k和Σ_k分别为第k个光斑的权重系数、均值和协方差矩阵，k＝1,2,...,K为GMM中高斯成分的序号，N(m|μ_k,Σ_k)表示μ_k和Σ_k为参数的随机变量m的高斯分布，高斯分布成分数即为区域内星体光斑数总数为K。利用EM算法根据样本点求解得到的K个高斯成分的参数π_k、μ_k和Σ_k，而μ_k为对应的第k个星点光斑的质心。由于EM算法为迭代优化算法，所以需要先根据以下步骤2、步骤3对各参数进行初始化。

$p\left(m\right)=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{k}N\left(m\right|{\mathrm{\&mu;}}_k,{\mathrm{\&Sigma;}}_k)---\left(6\right)$

2.初始化每个弥散光斑的中心。每个弥散光斑的质心，即K个高斯成分的均值μ_k,k＝1,2,......,K，被初始化化为MSER检测中得到的极大灰度值像素的中心位置。

3.根据每个样本点到各个弥散光斑的初始化中心的欧式距离，为每个样本点选取最近中心点为类别中心，把样本点分为K类，并且把属于第k类的样本集合记为S_Ck。样本点集合S_C被分割为K个集合，S_Ck中像素点样本被记为 上标表示其属于S_Ck集合，下标为S_Ck集合元素序号。通过公式(7)初始化所有的π_k，其中 运算符|·|表示求集合的元素个数。

${\mathrm{\&pi;}}_k=\frac{\left|S_{\mathrm{Ck}}\right|}{L}---\left(7\right)$

根据公式(8)对K个高斯分布成分的协方差矩阵进行初始化

${\mathrm{\&Sigma;}}_k=\frac{1}{\left|S_{\mathrm{Ck}}\right|}\underset{n_k=1}{\overset{\left|S_{\mathrm{Ck}}\right|}{\mathrm{\&Sigma;}}}(m_{n_k}^k-{\mathrm{\&mu;}}_k){(m_{n_k}^k-{\mathrm{\&mu;}}_k)}^T---\left(8\right)$

4.在完成初始化工作后，进行质心估计，高斯混合模型的质心使用EM算法来估计，EM算法为迭代算法，每次迭代步骤分为两步，第一步为期望估计，第二步为最大化。第一次迭代为初始化，对于第α次迭代，α≥2，具体如下：

a)通过初始化的参数计算第n个样本点m_n由第个k个高斯模型生成的概率

${\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}=\frac{{\mathrm{\&pi;}}_{k}N\left(m_n\right|{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}-1},{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}-1})}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{j}N\left(m_n\right|{\mathrm{\&mu;}}_j^{\mathrm{\&alpha;}-1},{\mathrm{\&Sigma;}}_j^{\mathrm{\&alpha;}-1})}---\left(9\right)$

估计更新参数，包括方差、均值和权值系数。具体通过公式(10)计算新的均值，通过公式(11)计算新的协方差方差矩阵，通过公式(12)得到新的权重系数。

${\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}}=\frac{1}{N_k^{\mathrm{\&alpha;}}}\underset{n=1}{\overset{\left|S_C\right|}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}{m}_n---\left(10\right)$

${\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}}=\frac{1}{N_k^{\mathrm{\&alpha;}}}\underset{n=1}{\overset{\left|S_C\right|}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}(m_n-{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}}){(m_n-{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}})}^T---\left(11\right)$

${\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}}=\frac{N_k^{\mathrm{\&alpha;}}}{L}---\left(12\right)$

其中计算方法如公式(13)所示，L为训练样本的数目。

$N_k^{\mathrm{\&alpha;}}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}---\left(13\right)$

迭代执行步骤b)，对样本点集合S_C，计算在不同参数下的对数似然函数值，当迭代到满足不等式(14)时，迭代结束。t_gmm是判断两次对数似然值是否收敛的阈值。其中 与分布表示模型参数与的似然值。

$|\ln p\left(S_C\right|{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}})-\ln p\left(S_C\right|{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}-1},{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}-1},{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}-1})|<t_{\mathrm{gmm}}---\left(14\right)$

其中的似然值计算方式如公式(15)所示

$\ln p\left(S_C\right|{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}})=\underset{n=1}{\overset{\left|S_C\right|}{\mathrm{\&Sigma;}}}\ln \left\{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}}N\left(m_n\right|{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}})\right\}---\left(15\right)$

对的计算同理。 

迭代停止后得到的K个高斯成分的均值μ_k,k＝1,...,K为空间碎片区域Ω中的K个空间碎片星点光斑的亚像素坐标定位结果。依次遍历联通区域集合中所有空间碎片区域，完成对图像I中所有空间碎片星点目标的提取与定位工作。

本发明的有益效果是：该方法通过MSER算法在候选空间碎片区域内，通过构建数据结构树，提取出该区域内空间碎片像素集合并确定空间碎片星点个数，为利用高斯混合模型对星点光斑建模和利用最大期望EM算法求解是的初始化提供了基础。在MSER算法得到的像素点树形结构中，通过检测该候选空间碎片区域中最大的MSER可以为生成后续拟合GMM模型生成训练样本。在基于GMM的定位过程中，通过EM算法迭代求解混合模型中高斯成分的均值与协方差，对应的均值就是各个空间碎片星点光斑的质心。该方法结合MSER和GMM模型的特点，通过GMM模型对近距离弥散光斑质心进行共同求解，避免了在空间碎片目标的弥散光斑靠近、黏连甚至部分重叠时分割的困难。

以下结合具体实施方式详细说明本发明。

具体实施方式

本发明空间碎片星点提取与定位方法具体步骤如下：

(a)星空图像二值化处理：

对于星空图像I，首先对针对星空图像中每个像素统计图像的灰度直方图，根据像素值按升序对图像中的像素点进行排序，去除像素灰度值最小的λ₁％的像素点，这样就能减弱暗噪声对背景的影响，去除像素灰度值最高的λ₂％的像素点，用于减弱强噪声和恒星对背景估计的影响，其中λ₁,λ₂为设定的背景阈值参数。利用剩余的(100-λ₁-λ₂)％像素点，求得背景的灰度均值μ和灰度标准差σ，使用阈值μ+κσ对图像进行二值化处理，这里κ为需要设置的分割系数。二值化处理公式为：

$B(x,y)=\left\{\begin{array}{cc}1& I(x,y)>\mathrm{\&mu;}+\mathrm{\&kappa;\&sigma;}\\ 0& I(x,y)\&le;\mathrm{\&mu;}+\mathrm{\&kappa;\&sigma;}\end{array}\right.---\left(1\right)$

式中，I(x,y)表示在坐标为(x,y)处的图像I中像素的灰度值，这样B(x,y)就是(x,y)处像素点的二值化结果，把二值化的后的图像记为B。这一步为星图的二值化粗分割， 因此每个B(x,y)＝1的连通区域中可能包含多于一个的星点。

设置参数λ₁＝20与λ₂＝30，在剩余的百分之五十的像素点中，求得该星空图像的背景均值与方差。令参数κ＝3，就可以得到该星空图像的二值化分割阈值为μ+3σ，利用该阈值与公式(1)，就可以得到二值化分割后的星空图像B。

(b)空间碎片星点区域粗分割提取：

由于后续的极大灰度值检测步骤需要利用粗分得到的空间碎片星点区域，所以需要对空间碎片区域进行粗分割。根据上一步图像二值化的结果，当B(x,y)＝1时表示该像素为空间碎片星点区域，当B(x,y)＝0时表示该像素为非有效空间碎片星点区域，即背景区域。对B中B(x,y)＝1的点进行广度优先搜索，可以提取到所有满足B(x,y)＝1的点的连通域，所有连通域集合为其中一个连通域可以表示为Ω。Ω为一个空间碎片星点区域，包含了一系列二值化结果为1的像素点坐标，并且它们相互连通，而且对于任一个Ω中的坐标为(x,y)的点，它在原始图像中的的灰度值为I(x,y)。

(c)星点区域局部极大灰度值和最大稳定极值区域检测：

在提取空间碎片区域后，需要对每一个连通域Ω通过MSER算法进行局部极大灰度值和最大稳定极值区域检测。MSER检测在这里的核心思想是通过不同阈值对空间碎片区域进行二值化处理，根据不同的结果对稳定极大灰度值点和区域进行检测。对于一个星点区域Ω，当令二值化阈值为t时，公式(2)中的Ω_>t表示在空间碎片区域Ω中，像素灰度值大于t的像素点集合，当二值化阈值t不断变大的时候集合Ω_>t的元素个数会越来越少，同理t变小的时候集合Ω_>t的元素个数会变多。在MSER检测过程中通过集合中元素个数的变化速率来表示一个区域Ω_>t是否是一个稳定的空间碎片弥散光斑。公式(3)表示当分割阈值变化大小为d时，光斑变化的稳定度计算方法。这里|Ω_>t|，|Ω_>t+d|和|Ω_>t|分布表示集合Ω_>t，Ω_>t+d和Ω_>t的元素个数，d为需要设置参数，st为光斑的稳定度。

Ω_>t＝{(x,y)|I(x,y)＞t,(x,y)∈Ω}   (2)

$\mathrm{st}=\frac{\left|{\mathrm{\&Omega;}}_{>t-d}\right|-\left|{\mathrm{\&Omega;}}_{>t+d}\right|}{\left|{\mathrm{\&Omega;}}_t\right|}---\left(3\right)$

利用MSER检测算法在一个空间碎片星点连通区域Ω中分割提取星点的步骤为：

1、首先对于图像中某个空间碎片星点连通区域Ω内的每个像素点灰度值按照大小 排序，依次初始化他们为未访问状态再根据各个像素点的灰度值大小和他们空间位置关系，按照步骤2、步骤3的描述把它们组织成一颗树形结构，这个树的各个节点都代表了图像中的一个像素点。

2、按照像素灰度值从大到小遍历Ω中所有的像素，针对某个(x,y)处的像素点遍历它周围相邻的8个点。假设(x′,y′)为一个与(x,y)相邻的位置，那么如果(x′,y′)处像素点未被遍历，则直接遍历下一个相邻点，如果已经被遍历那么：

a)当Root(x,y)＝Root(x′,y′)相同时，则不进行任何处理，遍历(x,y)处像素的下一个相邻点，因为(x,y)处像素点与(x′,y′)处像素点已经处于同一个子树中。这里的Root(x,y)表示(x,y)处像素点在当前树结构中对应的根像素点的位置。

b)当I(Root(x,y))＝I(Root(x′,y′))时，则把Root(x′,y′)处像素点作为子节点添加到Root(x,y)处像素点下。

c)当I(Root(x,y))<I(Root(x′,y′))时，则把Root(x′,y′)处像素点添加为Root(x,y)处像素点的子节点。

3、标记当前(x,y)处像素点为已被遍历像素点，并重新回到步骤2直到完成对所有像素点的遍历。

4、当所有的像素点按照上面步骤遍历完毕，依照像素点灰度值大小关系的树形结构建立完成，该树中叶节点个数就等于该空间碎片区域Ω中局部极大灰度值像素点个数，也对应了区域内星点光斑数目，这个数目记为K，极值点的个数、位置将用于以下基于GMM模型和EM迭代算法的星点定位的初始化。

为了为使用EM算法拟合GMM模型提供样本，需要在空间碎片区域内检测最大稳定极值区域。对于空间碎片光斑区域Ω内，每当取到一个二值化分割阈值t，总是可以对空间碎片区域Ω进行分割得到一个子光斑为Ω_>t，当该光斑的稳定度st小于t_mser时，即满足公式(4)时，该子光斑就是为Ω的一个MSER，这里t_mser为判断一个子光斑是否为MSER的阈值。当调整t，可以得到一系列MSER，这一系列的MSER呈现树形的包含现象，即t较小时分割得到的MSER包含t较大时分割得到的MSER，而最大的MSER包含所有其他的MSER。对区域Ω提取不同t情况下的MSER之后得到MSER集合其中包含像素点最多的MSER记为C。

$\frac{\left|{\mathrm{\&Omega;}}_{>t-d}\right|-\left|{\mathrm{\&Omega;}}_{t+d}\right|}{\left|{\mathrm{\&Omega;}}_{>}t\right|}<t_{\mathrm{mser}}---\left(4\right)$

令公式(4)中的d＝5，MSER判断阈值为t_mser＝0.3，那么所有满足公式(4)的区域都可以认为是一个MSER，这就可得到包含空间碎片区域Ω所有的MSER的集合其中包含最多像素点的MSER为C，它为后续质心定位提供采样基础。

(d)空间碎片定位：

考虑单个空间碎片弥散光斑可以被合理假设为一个高斯函数形状的的点状光斑，对于包含多个弥散光斑的空间碎片区域，利用混合高斯模型(GMM)对空间碎片区域内的多个弥散光斑进行建模。GMM在质心定位的时候需要多个像素点样本。因为图像中每个像素点的灰度值与像元接收光子数成正比，如果某个位置处像素点灰度值越大说明越接近光子分布的空间中心位置，所以使用重采样技术，对一个像素点进行重采样得到与其灰度值成正比的数目的像素点样本，这些重采样得到的像素点与原像素点位置和灰度值相同。

C为空间碎片区域Ω的MSER集合中包含最多像素点的MSER，对其范围内任一位置(x,y)处的像素点进行重采样，得到个样本点，这个样本点的集合记为S_(x,y)，为向上取整符号。那么C中所有的样本点集合记为S_C，其中n＝1,...,L表示像素点样本的序号，L表示样本点个数，即集合S_C的中的成员数目，m_n表示第n个样本像素点的位置。

$S_C=\left\{m_n\right\}=\underset{(x,y)\&Element;c}{\&cup;}{S}_{(x,y)},n=1,...,L---\left(5\right)$

获得了样本点集合S_C和星点光斑数目K，就可以利用GMM模型和相应的EM算法对空间碎片光斑区域Ω内K个光斑同时进行定位。具体的模型描述和定位算法如下：

1、定义用于拟合空间碎片区域星点光斑的混合高斯模型。由于空间碎片区域是有多个弥散光斑组成，首先将对于空间碎片区域中的任一位置点变量(x,y)记为随机变量m＝(x,y)，它服从公式(6)描述的高斯混合分布，其中π_k、μ_k和Σ_k分别为第k个光斑的权重系数、均值和协方差矩阵，k＝1,2,...,K为GMM中高斯成分的序号，N(m|μ_k,Σ_k)表示μ_k和Σ_k为参数的随机变量m的高斯分布，高斯分布成分数即为区域内星体光斑数总数为K。利用EM算法根据样本点求解得到的K个高斯成分的参数π_k、μ_k和Σ_k，而μ_k为对应的第k个星点光斑的质心。由于EM算法为迭代优化算法，所以需要先根据 以下步骤2、步骤3对各参数进行初始化。

$p\left(m\right)=\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{k}N\left(m\right|{\mathrm{\&mu;}}_k,{\mathrm{\&Sigma;}}_k)---\left(6\right)$

2、初始化每个弥散光斑的中心。每个弥散光斑的质心，即K个高斯成分的均值μ_k(k＝1,2,......,K)，被初始化化为MSER检测中得到的极大灰度值像素的中心位置。

3、根据每个样本点到各个弥散光斑的初始化中心的欧式距离，为每个样本点选取最近中心点为类别中心，可以把样本点分为K类，并且把属于第k类的样本集合记为S_Ck。样本点集合S_C被分割为K个集合，S_Ck中像素点样本被记为n_k＝1,...,|S_Ck|，上标表示其属于S_Ck集合，下标为S_Ck集合元素序号。通过公式(7)初始化所有的π_k，其中运算符|·|表示求集合的元素个数。

${\mathrm{\&pi;}}_k=\frac{\left|S_{\mathrm{Ck}}\right|}{L}---\left(7\right)$

根据公式(8)对K个高斯分布成分的协方差矩阵进行初始化

${\mathrm{\&Sigma;}}_k=\frac{1}{\left|S_{\mathrm{Ck}}\right|}\underset{n_k=1}{\overset{\left|S_{\mathrm{Ck}}\right|}{\mathrm{\&Sigma;}}}(m_{n_k}^k-{\mathrm{\&mu;}}_k){(m_{n_k}^k-{\mathrm{\&mu;}}_k)}^T---\left(8\right)$

4、在完成初始化工作后，接下来就是进行质心估计，高斯混合模型的质心一般使用EM算法来估计，EM算法为迭代算法，每次迭代步骤分为两步，第一步为期望估计，第二步为最大化。第一次迭代为初始化，对于第α次迭代，α≥2，具体如下：

b)通过初始化的参数计算第n个样本点m_n由第个k个高斯模型生成的概率

${\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}=\frac{{\mathrm{\&pi;}}_{k}N\left(m_n\right|{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}-1},{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}-1})}{\underset{j=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_{j}N\left(m_n\right|{\mathrm{\&mu;}}_j^{\mathrm{\&alpha;}-1},{\mathrm{\&Sigma;}}_j^{\mathrm{\&alpha;}-1})}---\left(9\right)$

估计更新参数，包括方差，均值和权值系数。具体通过公式(10)计算新的均值，通过公式(11)计算新的协方差方差矩阵，通过公式(12)得到新的权重系数。

${\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}}=\frac{1}{N_k^{\mathrm{\&alpha;}}}\underset{n=1}{\overset{\left|S_C\right|}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}{m}_n---\left(10\right)$

${\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}}=\frac{1}{N_k^{\mathrm{\&alpha;}}}\underset{n=1}{\overset{\left|S_C\right|}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}(m_n-{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}}){(m_n-{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}})}^T---\left(11\right)$

${\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}}=\frac{N_k^{\mathrm{\&alpha;}}}{L}---\left(12\right)$

其中计算方法如公式(13)所示，L为训练样本的数目。

$N_k^{\mathrm{\&alpha;}}=\underset{n=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_{\mathrm{nk}}^{\mathrm{\&alpha;}}---\left(13\right)$

迭代执行步骤b)，对样本点集合S_C，计算在不同参数下的对数似然函数值，当迭代到满足不等式(14)时，迭代结束。t_gmm是判断两次对数似然值是否收敛的阈值。其中 与分布表示模型参数与的似然值。

$|\ln p\left(S_C\right|{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}})-\ln p\left(S_C\right|{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}-1},{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}-1},{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}-1})|<t_{\mathrm{gmm}}---\left(14\right)$

其中的似然值计算方式如公式(15)所示

$\ln p\left(S_C\right|{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}})=\underset{n=1}{\overset{\left|S_C\right|}{\mathrm{\&Sigma;}}}\ln \left\{\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_k^{\mathrm{\&alpha;}}N\left(m_n\right|{\mathrm{\&mu;}}_k^{\mathrm{\&alpha;}},{\mathrm{\&Sigma;}}_k^{\mathrm{\&alpha;}})\right\}---\left(15\right)$

对的计算同理。

迭代停止后得到的K个高斯成分的均值μ_k(k＝1,...,K)为空间碎片区域Ω中的K个空间碎片星点光斑的亚像素坐标定位结果。依次遍历联通区域集合中所有空间碎片区域，完成对图像I中所有空间碎片星点目标的提取与定位工作。

Claims (1)

1.一种空间碎片星点提取与定位方法，其特征在于包括以下步骤：

步骤一、对于星空图像I，首先对针对星空图像中每个像素统计图像的灰度直方图，根据像素值按升序对图像中的像素点进行排序，去除像素灰度值最小的λ₁％的像素点，以减弱暗噪声对背景的影响，去除像素灰度值最高的λ₂％的像素点，以减弱强噪声和恒星对背景估计的影响，其中λ₁,λ₂为设定的背景阈值参数；利用剩余的(100-λ₁-λ₂)％像素点，求得背景的灰度均值μ和灰度标准差σ，使用阈值μ+κσ对图像进行二值化处理，κ为需要设置的分割系数；二值化处理公式为：

式中，I(x,y)表示在坐标为(x,y)处的图像I中像素的灰度值，B(x,y)就是(x,y)处像素点的二值化结果，把二值化的后的图像记为B；这一步为星图的二值化粗分割，因此每个B(x,y)＝1的连通区域中可能包含多于一个的星点；

步骤二、对空间碎片区域进行粗分割；根据上一步图像二值化的结果，当B(x,y)＝1时表示该像素为空间碎片星点区域，当B(x,y)＝0时表示该像素为非有效空间碎片星点区域，即背景区域；对B中B(x,y)＝1的点进行广度优先搜索，提取到所有满足B(x,y)＝1的点的连通域，所有连通域集合为其中一个连通域表示为Ω；Ω为一个空间碎片星点区域，包含了一系列二值化结果为1的像素点坐标，并且他们相互连通，而且对于任一个Ω中的坐标为(x,y)的点，它在原始图像中的灰度值为I(x,y)；

步骤三、在提取空间碎片区域后，需要对每一个连通域Ω通过MSER算法进行局部极大灰度值和最大稳定极值区域检测；MSER检测是通过不同阈值对空间碎片区域进行二值化处理，根据不同的结果对稳定极大灰度值点和区域进行检测；对于一个星点区域Ω，当令二值化阈值为t时，公式(2)中的Ω_>t表示在空间碎片区域Ω中，像素灰度值大于t的像素点集合，当二值化阈值t不断变大的时候集合Ω_>t的元素个数会越来越少，同理t变小的时候集合Ω_>t的元素个数会变多；在MSER检测过程中通过集合中元素个数的变化速率来表示一个区域Ω_>t是否是一个稳定的空间碎片弥散光斑；公式(3)表示当分割阈值变化大小为d时，光斑变化的稳定度计算方法；这里|Ω_>t|，|Ω_>t+d|和|Ω_>t|分布表示集合Ω_>t，Ω_>t+d和Ω_>t的元素个数，d为需要设置参数，st 为光斑的稳定度；

Ω_>t＝{(x,y)|I(x,y)>t,(x,y)∈Ω}         (2) 

利用MSER算法在一个空间碎片星点连通区域Ω中分割提取星点的步骤为：

（1）.首先对于图像中某个空间碎片星点连通区域Ω内的每个像素点灰度值按照大小排序，依次初始化他们为未访问状态再根据各个像素点的灰度值大小和他们空间位置关系，按照步骤（2）、步骤（3）的描述把他们组织成一颗树形结构，这个树的各个节点都代表了图像中的一个像素点；

（2）.按照像素灰度值从大到小遍历Ω中所有的像素，针对某个(x,y)处的像素点遍历它周围相邻的8个点；假设(x′,y′)为一个与(x,y)相邻的位置，那么如果(x′,y′)处像素点未被遍历，则直接遍历下一个相邻点，如果已经被遍历那么：

a)当Root(x,y)＝Root(x′,y′)相同时，则不进行任何处理，遍历(x,y)处像素的下一个相邻点，因为(x,y)处像素点与(x′,y′)处像素点已经处于同一个子树中；这里的Root(x,y)表示(x,y)处像素点在当前树结构中对应的根像素点的位置；

b)当I(Root(x,y))＝I(Root(x′,y′))时，则把Root(x′,y′)处像素点作为子节点添加到Root(x,y)处像素点下；

c)当I(Root(x,y))<I(Root(x′,y′))时，则把Root(x′,y′)处像素点添加为Root(x,y)处像素点的子节点；

（3）.标记当前(x,y)处像素点为已被遍历像素点，并重新回到步骤（2）直到完成对所有像素点的遍历；

（4）.当所有的像素点按照上面步骤遍历完毕，依照像素点灰度值大小关系的树形结构建立完成，该树中叶节点个数就等于该空间碎片区域Ω中局部极大灰度值像素点个数，也对应了区域内星点光斑数目，这个数目记为K，极值点的个数、位置将用于以下基于GMM模型和EM迭代算法的星点定位的初始化；

为了使用EM算法拟合GMM模型提供样本，需要在空间碎片区域内检测最大稳定极值区域；对于空间碎片光斑区域Ω内，每当取到一个二值化分割阈值t，总是对空间碎片区域Ω进行分割得到一个子光斑为Ω_>t，当该光斑的稳定度st小于t_mser时，即满 足公式(4)时，该子光斑就是为Ω的一个MSER，这里t_mser为判断一个子光斑是否为MSER的阈值；当调整t，就得到一系列MSER，这一系列的MSER呈现树形的包含现象，即t较小时分割得到的MSER包含t较大时分割得到的MSER，而最大的MSER包含所有其他的MSER；对区域Ω提取不同t情况下的MSER之后得到MSER集合其中包含像素点最多的MSER记为C；

步骤四、考虑单个空间碎片弥散光斑可以被合理假设为一个高斯函数形状的的点状光斑，对于包含多个弥散光斑的空间碎片区域，利用混合高斯模型对空间碎片区域内的多个弥散光斑进行建模；GMM在质心定位的时候需要多个像素点样本；

C为空间碎片区域Ω的MSER集合中包含最多像素点的MSER，对其范围内任一位置(x,y)处的像素点进行重采样，得到个样本点，这个样本点的集合记为S_(x,y)，为向上取整符号；那么C中所有的样本点集合记为S_C，其中n＝1,...,L表示像素点样本的序号，L表示样本点个数，即集合S_C的中的成员数目，m_n表示第n个样本像素点的位置；

获得了样本点集合S_C和星点光斑数目K，利用GMM模型和相应的EM算法对空间碎片光斑区域Ω内K个光斑同时进行定位；具体的模型描述和定位算法如下：

（1）.定义用于拟合空间碎片区域星点光斑的混合高斯模型；由于空间碎片区域是有多个弥散光斑组成，首先将对于空间碎片区域中的任一位置点变量(x,y)记为随机变量m＝(x,y)，它服从公式(6)描述的高斯混合分布，其中π_k、μ_k和Σ_k分别为第k个光斑的权重系数、均值和协方差矩阵，k＝1,2,...,K为GMM中高斯成分的序号，N(m|μ_k,Σ_k)表示μ_k和Σ_k为参数的随机变量m的高斯分布，高斯分布成分数即为区域内星体光斑数总数为K；利用EM算法根据样本点求解得到的K个高斯成分的参数π_k、μ_k和Σ_k，而μ_k为对应的第k个星点光斑的质心；由于EM算法为迭代优化算法，所以需要先根据以下步骤（2）、步骤（3）对各参数进行初始化；

（2）.初始化每个弥散光斑的中心；每个弥散光斑的质心，即K个高斯成分的均值μ_k,k＝1,2,......,K，被初始化化为MSER检测中得到的极大灰度值像素的中心位置；

（3）.根据每个样本点到各个弥散光斑的初始化中心的欧式距离，为每个样本点选取最近中心点为类别中心，把样本点分为K类，并且把属于第k类的样本集合记为S_Ck；样本点集合S_C被分割为K个集合，S_Ck中像素点样本被记为上标表示其属于S_Ck集合，下标为S_Ck集合元素序号；通过公式(7)初始化所有的π_k，其中运算符|·|表示求集合的元素个数；

根据公式(8)对K个高斯分布成分的协方差矩阵进行初始化

（4）.在完成初始化工作后，进行质心估计，高斯混合模型的质心使用EM算法来估计，EM算法为迭代算法，每次迭代步骤分为两步，第一步为期望估计，第二步为最大化；第一次迭代为初始化，对于第α次迭代，α≥2，具体如下：

a)通过初始化的参数计算第n个样本点m_n由第个k个高斯模型生成的概率

估计更新参数，包括方差、均值和权值系数；具体通过公式(10)计算新的均值，通过公式(11)计算新的协方差方差矩阵，通过公式(12)得到新的权重系数；

其中计算方法如公式(13)所示，L为训练样本的数目；

迭代执行步骤b)，对样本点集合S_C，计算在不同参数下的对数似然函数值，当迭代到满足不等式(14)时，迭代结束；t_gmm是判断两次对数似然值是否收敛的阈值；其中 与分布表示模型参数与的似然值；

其中的似然值计算方式如公式(15)所示

对的计算同理；

迭代停止后得到的K个高斯成分的均值μ_k,k＝1,...,K为空间碎片区域Ω中的K个空间碎片星点光斑的亚像素坐标定位结果；依次遍历联通区域集合中所有空间碎片区域，完成对图像I中所有空间碎片星点目标的提取与定位工作。