CN104931504A - 一种基于非负矩阵分解的空间碎片材料分析方法 - Google Patents

Abstract

一种基于非负矩阵分解的空间碎片材料分析方法，包括如下步骤：步骤1，碎片光谱数据采集；步骤2，光谱数据预处理；(2.1)截去光谱两端高噪声波段；(2.2)阈值化：(2.3)例外点去除；步骤3，非负矩阵分解：给定非负矩阵V，寻找非负矩阵因子W和H，使得V≈WH；为评估近似性能，给出代价函数及其相应迭代更新规则，计算广义K_L散度：步骤4，碎片材料谱辨认；通过选择不同的初始化W和H值，得到碎片材料谱的候选者，利用已知的地面材料谱，由步骤4计算候选材料谱与地面材料谱的广义K_L散度，选取散度阈值D₀，最后确定小于该散度阈值的候选谱为碎片材料谱的组成。本发明有效适用于深空环境、准确性较好。

Description

一种基于非负矩阵分解的空间碎片材料分析方法

技术领域

本发明涉及载人航天实施空间监测技术领域，尤其是一种空间碎片材料分析方法。

背景技术

自1957年10月前苏联发射首颗人造地球卫星以来，人类几十年的空间探索活动产生了大量的空间碎片，对人类航天活动的安全造成严重威胁，成为空间环境的主要污染源，并在一定程度上对航天活动的正常开展产生了影响。

随着中国的经济发展以及国家安全的需要，航天空间活动将愈来愈多，并在不久的将来，对应用卫星的需求更可能大增。这些应用卫星与国民经济关系密切，一旦受损，社会影响、经济影响巨大，甚至危及国家安全。而这些卫星运行区域大都分布在低轨道，处于空间碎片密集区域，受碰撞损伤的威胁很高。同时，持续开展载人航天活动、建立永久的有人值守轨道空间站，将会是中国航天事业发展的必然趋势。可以预期，对载人航天实施空间监测、预警也会成为中国航天不可回避的事实，在这样的背景下，对空间碎片的观测技术进行研究，具有重要的现实意义和研究价值。

传统空间碎片测量以位置信息测量为主，包括碎片的三维位置坐标、速度、加速度等参数,可衍生出各类地球轨道参数。为了提高碎片的监测预警能力,对测量系统除了要求获得其位置信息之外,更需要获得碎片的特征信息，如碎片的形状、体积、表面材料参数等特征信息，它对空间监测、预警尤为重要。

人造天体碎片的大小、形状、材料类别的确定，对于空间环境监测和预警至关重要。光谱测量技术是天体碎片分析的一种重要方法，人造天体碎片一般本身并不发光，其亮度来自太阳光的反射,仅由亮度变化不足以分辨碎片的材料组成，分析其光谱特征成为我们辨认碎片类别的主要手段。

在地基(地面观测站)条件下，光谱仪所获取光谱由碎片材料成分、太阳光谱、地球大气吸收谱，以及测量过程中产生的噪声等因素综合决定，因而由光谱推断碎片种类并没有一个解析解。

空间碎片的测量有雷达测量和光学测量，雷达测量可以克服天气、太阳和及地影的影响，能全天候全天时工作，但由于雷达测量时，其反射回波信号的强度与距离的四次方成反比，因而雷达测量比较适合于低地球轨道的小碎片；另外雷达测量需要发射信号，属于主动探测形式，这在某些情形下可能不合适。

而对于无源光学测量，信号反射强度与物体的距离的平方成反比，因而能探测高轨道碎片，另外它仅接收碎片对太阳光的反射，属于被动探测形式，可较好的适合于某些特殊场合。

基于光谱数据对碎片材料进行分析，所观测的碎片光谱是多种材料光谱的混合，而材料光谱的混合机制仍很复杂，它与碎片的姿势、材料组成、光照条件、观测角度等因素有关，甚至真空环境、辐照强度、材料老化都会改变混合谱形态，这使得由混合谱计算材料谱及混合比例成为一个典型的病态问题，并没有一个解析解。

发明内容

为了克服已有空间碎片材料分析方法的无法适用于深空环境、准确性较差的不足，本发明提供一种有效适用于深空环境、准确性较好的基于非负矩阵分解的空间碎片材料分析方法。

本发明解决其技术问题所采用的技术方案是：

一种基于非负矩阵分解的空间碎片材料分析方法，包括如下步骤：

步骤1，碎片光谱数据采集；

利用地基望远镜和无狭缝光栅，采集碎片光谱数据；

步骤2，光谱数据预处理；

(2.1)截去光谱两端高噪声波段：

将所测得的可见光波段[403-779]纳米，截去两端高噪声波段，剩下波长范围为[479-729]纳米；

(2.2)阈值化：

采用阈值化处理：

$f\left(x\right)=\left\{\begin{array}{cc}\mathrm{Th}& \mathrm{if\; x}>\mathrm{Th}\\ 0& \mathrm{if\; x}<0\\ x& \mathrm{otherwise}\end{array}\right.$

其中，f(x)为阈值化函数，选取阈值Th，若光谱数据x大于该阈值，则设置该项数据为所选择的阈值Th；若光谱数据x小于0，则设置该值为0；其它情形，保持不变；

(2.3)例外点去除

取同一圈次中观测所得的光谱数据，计算各光谱间的欧氏距离，并求得欧氏距离的均值和标准差δ；计算每条光谱与同一圈次中其它光谱间的欧氏距离，并求其均值i为光谱序号，若则将第i条光谱作为例外点去除；

步骤3，非负矩阵分解；

给定非负矩阵(m行，n列矩阵)，寻找非负矩阵因子(m行，r列矩阵)和(r行，n列矩阵)，使得V≈WH；其中表示矩阵元素为非负实数，r为分解得到的基矢量个数，

为评估近似性能，给出代价函数及其相应迭代更新规则：

广义K_L散度： $D\left(V\right|\left|\left(\mathrm{WH}\right)\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}(V_{\mathrm{ij}}\log \frac{V_{\mathrm{ij}}}{{\left(\mathrm{WH}\right)}_{\mathrm{ij}}}-V_{\mathrm{ij}}+{\left(\mathrm{WH}\right)}_{\mathrm{ij}});$

$H_{\mathrm{\&alpha;\&mu;}}\&LeftArrow;H_{\mathrm{\&alpha;\&mu;}}\frac{{\mathrm{\&Sigma;}}_{i=1}^m{W}_{\mathrm{i\&alpha;}}{V}_{\mathrm{i\&mu;}}/{\left(\mathrm{WH}\right)}_{\mathrm{i\&mu;}}}{{\mathrm{\&Sigma;}}_{k=1}^r{W}_{\mathrm{k\&alpha;}}},W_{\mathrm{i\&alpha;}}\&LeftArrow;W_{\mathrm{i\&alpha;}}\frac{{\mathrm{\&Sigma;}}_{\mathrm{\&mu;}=1}^n{H}_{\mathrm{\&alpha;\&mu;}}{V}_{\mathrm{i\&mu;}}/{\left(\mathrm{WH}\right)}_{\mathrm{i\&mu;}}}{{\mathrm{\&Sigma;}}_{v=1}^n{H}_{\mathrm{\&alpha;v}}}$

其中α,μ,i,j,k,ν为矩阵元素的下标，Σ表示求和；

步骤4，碎片材料谱辨认；

通过选择不同的初始化W和H值，得到碎片材料谱的候选者，利用已知的地面材料谱，由步骤4计算候选材料谱与地面材料谱的广义K_L散度，选取散度阈值D₀，最后确定小于该散度阈值的候选谱为碎片材料谱的组成，记为W₀。

进一步，所述方法还包括以下步骤：步骤5，碎片材料丰度计算；

利用最小非负二乘法，其中，表示计算矢量各元素的平方和，V_i为第i条观测光谱矢量，W₀为步骤4所确定的碎片材料光谱，可求得材料的丰度比例矢量X，各分量即为对应材料的丰度比例。

本发明的技术构思为：采用线性混合的模型对材料光谱模型建模，并基于非负矩阵分解(NMF)，采用多次迭代方式求得材料谱候选者，再用广义K_L散度确定材料组成；最后用非负最小二乘算法得到材料丰度比例；本方法可克服深空环境雷达难以准确测量导致空间碎片类型难以辩析的困难，对于低轨碎片类型的辨认也提供了一种有效的辅助分析手段。

本发明的有益效果主要表现在：利用地基望远镜获取碎片的光谱数据，然后利用该算法辨别碎片材料组成，最后给出计算材料组成丰度方法；本发明相比雷达等其他方法的优点是：较适宜于深空环境，因为雷达难以测量远距离对象，同时对低轨碎片类型的确定提供了一种有效的辅助手段。

附图说明

图1是本发明的方法流程图。

图2是本发明的实施例的预处理前后比较结果示意图，其中，(a)为预处理前的示意图，(b)为预处理后的示意图。

图3是本发明的实施例的采集的地面材料谱图。

图4是本发明的实施例的以铝(Aluminum)材料为主的混合谱图。

图5是本发明的实施例的两种材料混合丰度的示意图。

图6是本发明的实施例的NMF估计的材料谱与真实混合材料谱比较图，其中，(a)表示铝(Aluminum)，(b)表示聚酯薄膜(Mylar)。

图7是本发明的实施例的NMF估计的两种材料丰度相对误差示意图，其中，(a)表示铝(Aluminum)，(b)表示聚酯薄膜(Mylar)。

具体实施方式

下面结合附图对本发明作进一步描述。

参照图1～图7，一种基于非负矩阵分解的空间碎片材料分析方法，分为两个过程：碎片材料组成分析过程和碎片材料丰度计算过程；其流程图见图1：

包括如下步骤：

步骤1，碎片光谱数据采集；

利用地基望远镜和无狭缝光栅，采集碎片光谱数据；

步骤2，光谱数据预处理；

(1)截去光谱两端高噪声波段：

将所测得的可见光波段[403-779]纳米，截去两端高噪声波段，剩下波长范围为[479-729]纳米；

(2)阈值化：

由于测量精度及大气湍流等因素影响，所获取光谱数据中存在某些误差，为此采用如下阈值化处理：

$f\left(x\right)=\left\{\begin{array}{cc}\mathrm{Th}& \mathrm{if\; x}>\mathrm{Th}\\ 0& \mathrm{if\; x}<0\\ x& \mathrm{otherwise}\end{array}\right.$

其中，f(x)为阈值化函数，选取阈值Th，若光谱数据x大于该阈值，则设置该项数据为所选择的阈值Th；若光谱数据x小于0，则设置该值为0；其它情形，保持不变；

(3)例外点去除

取同一圈次中观测所得的光谱数据，计算各光谱间的欧氏距离，求得欧氏距离的均值和标准差δ；计算每条光谱与同一圈次中其它光谱间的欧氏距离，并求其均值i为光谱序号，若则将第i条光谱作为例外点去除。

步骤3，非负矩阵分解(NMF)；

给定非负矩阵(m行，n列矩阵)，寻找非负矩阵因子(m行，r列矩阵)和(r行，n列矩阵)，使得V≈WH；其中表示矩阵元素为非负实数，r为分解得到的基矢量个数，

为评估近似性能，给出代价函数及其相应迭代更新规则：

广义K_L散度： $D\left(V\right|\left|\left(\mathrm{WH}\right)\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}(V_{\mathrm{ij}}\log \frac{V_{\mathrm{ij}}}{{\left(\mathrm{WH}\right)}_{\mathrm{ij}}}-V_{\mathrm{ij}}+{\left(\mathrm{WH}\right)}_{\mathrm{ij}});$

$H_{\mathrm{\&alpha;\&mu;}}\&LeftArrow;H_{\mathrm{\&alpha;\&mu;}}\frac{{\mathrm{\&Sigma;}}_{i=1}^m{W}_{\mathrm{i\&alpha;}}{V}_{\mathrm{i\&mu;}}/{\left(\mathrm{WH}\right)}_{\mathrm{i\&mu;}}}{{\mathrm{\&Sigma;}}_{k=1}^r{W}_{\mathrm{k\&alpha;}}},W_{\mathrm{i\&alpha;}}\&LeftArrow;W_{\mathrm{i\&alpha;}}\frac{{\mathrm{\&Sigma;}}_{\mathrm{\&mu;}=1}^n{H}_{\mathrm{\&alpha;\&mu;}}{V}_{\mathrm{i\&mu;}}/{\left(\mathrm{WH}\right)}_{\mathrm{i\&mu;}}}{{\mathrm{\&Sigma;}}_{v=1}^n{H}_{\mathrm{\&alpha;v}}}$

其中α,μ,i,j,k,ν为矩阵元素的下标，Σ表示求和；

步骤4，碎片材料谱辨认；

通过选择不同的初始化W和H值，得到碎片材料谱的候选者，利用已知的地面材料谱，由公式(4)计算候选材料谱与地面材料谱的广义K_L散度，选取一定的阈值D₀，最后确定小于该阈值的候选谱为碎片材料谱的组成，记为W₀。

步骤5，碎片材料丰度计算；

利用最小非负二乘法，其中V为观测光谱，W₀为步骤4所确定的碎片材料光谱，可求得每种材料的丰度比例X。

实例：一种基于非负矩阵分解的空间碎片材料分析方法，过程如下：

(1).数据预处理，如图2所示：

(2).地面材料谱的收集，分别为铝(Aluminum)、聚酯薄膜(Mylar)、太阳板(Solar Cell)、白漆(White paint)，其光谱图见图3，其波长范围为[403--779]纳米。

(3).材料谱混合

以‘Aluminum’和‘Mylar’按一定比例混合，再加1％的高斯噪声，取混合数100，得到混合谱见图4所示，两者的丰度图见图5。

(4).NMF分解

由NMF所估计的材料组成见图6。

(5).由最小非负二乘发计算材料组成，估计误差由图7表示。其中Aluminum的丰度相对误差平均为3.4556％，Mylar的丰度相对误差平均为4.0675％。

本发明基于光谱数据对碎片材料进行分析，所观测的碎片光谱是多种材料光谱的混合，而材料光谱的混合机制仍很复杂，它与碎片的姿势、材料组成、光照条件、观测角度等因素有关，甚至真空环境、辐照强度、材料老化都会改变混合谱形态，这使得由混合谱计算材料谱及混合比例成为一个典型的病态问题，并没有一个解析解，本发明采用线性混合的模型对材料光谱模型建模，并基于非负矩阵分解(NMF)，采用多次迭代方式求得材料谱候选者，再用广义K_L散度确定材料组成，最后用非负最小二乘算法得到材料丰度比例；本发明可克服深空环境雷达难以准确测量导致空间碎片类型难以辩析的困难，对于低轨碎片类型的辨认，本方法也提供了一种有效的辅助分析手段，实验结果证明了本发明的可行性。

Claims (2)

1.一种基于非负矩阵分解的空间碎片材料分析方法，其特征在于：所述空间碎片材料分析方法包括如下步骤：

步骤1，碎片光谱数据采集；

利用地基望远镜和无狭缝光栅，采集碎片光谱数据；

步骤2，光谱数据预处理；

(2.1)截去光谱两端高噪声波段：

将所测得的可见光波段[403-779]纳米，截去两端高噪声波段，剩下波长范围为[479-729]纳米；

(2.2)阈值化：

采用阈值化处理：

$f\left(x\right)=\left\{\begin{array}{cc}\mathrm{Th}& \mathrm{ifx}>\mathrm{Th}\\ 0& \mathrm{ifx}<0\\ x& \mathrm{otherwise}\end{array}\right.$

其中，f(x)为阈值化函数，选取阈值Th，若光谱数据x大于该阈值，则设置该项数据为所选择的阈值Th；若光谱数据x小于0，则设置该值为0；其它情形，保持不变；

(2.3)例外点去除

取同一圈次中观测所得的光谱数据，计算各光谱间的欧氏距离，求得欧氏距离的均值和标准差δ；计算每条光谱与同一圈次中其它光谱间的欧氏距离，并求其均值i为光谱序号，若则将第i条光谱作为例外点去除；

步骤3，非负矩阵分解；

给定非负矩阵V为m行n列矩阵，寻找非负矩阵因子W和H，使得V≈WH；其中，表示矩阵元素为非负实数， W为m行r列矩阵，H为r行n列矩阵，r为分解得到的基矢量个数，

为评估近似性能，给出代价函数及其相应迭代更新规则：

广义K_L散度： $D\left(V\right|\left|\left(\mathrm{WH}\right)\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}(V_{\mathrm{ij}}\log \frac{V_{\mathrm{ij}}}{{\left(\mathrm{WH}\right)}_{\mathrm{ij}}}-V_{\mathrm{ij}}+{\left(\mathrm{WH}\right)}_{\mathrm{ij}});$

$H_{\mathrm{\&alpha;\&mu;}}\&LeftArrow;H_{\mathrm{\&alpha;\&mu;}}\frac{{\mathrm{\&Sigma;}}_{i=1}^m{W}_{\mathrm{i\&alpha;}}{V}_{\mathrm{i\&mu;}}/{\left(\mathrm{WH}\right)}_{\mathrm{i\&mu;}}}{{\mathrm{\&Sigma;}}_{k=1}^r{W}_{\mathrm{k\&alpha;}}},W_{\mathrm{i\&alpha;}}\&LeftArrow;W_{\mathrm{i\&alpha;}}\frac{{\mathrm{\&Sigma;}}_{\mathrm{\&mu;}=1}^n{H}_{\mathrm{\&alpha;\&mu;}}{V}_{\mathrm{i\&mu;}}/{\left(\mathrm{WH}\right)}_{\mathrm{i\&mu;}}}{{\mathrm{\&Sigma;}}_{v=1}^n{H}_{\mathrm{\&alpha;v}}}$

其中α,μ,i,j,k,v为矩阵元素的下标，Σ表示求和；

步骤4，碎片材料谱辨认；

通过选择不同的初始化W和H值，得到碎片材料谱的候选者，利用已知的地面材料谱，由步骤4计算候选材料谱与地面材料谱的广义K_L散度，选取散度阈值D₀，最后确定小于该散度阈值的候选谱为碎片材料谱的组成，记为W₀。

2.如权利要求1所述的一种基于非负矩阵分解的空间碎片材料分析方法，其特征在于：所述方法还包括以下步骤：

步骤5，碎片材料丰度计算；

利用最小非负二乘法，其中，表示计算矢量各元素的平方和，V_i为第i条观测光谱矢量，W₀为步骤4所确定的碎片材料光谱，可求得材料的丰度比例矢量X，各分量即为对应材料的丰度比例。