-
随着对空间资产和空间相关技术的日益依赖, 各国越来越重视空间势态感知技术的发展. 对于远距离的空间目标, 通常在探测器的像面上只占据几个像素; 仅凭这几个像素的灰度信息, 很难有效地获取空间目标的特征[1,2]. 能代表物体固有属性差异的红外光谱信息为空间目标的识别与特征提取带来了很大的希望[3,4].
美国波音公司Skinner等利用宽带阵列光谱系统(BASS)对人造地球卫星进行了长期的观测, 并对观测的红外光谱进行了分析. 2007年, 他们采集了三个地球同步卫星的红外光谱数据, 用普朗克函数进行了拟合, 提取了等效温度和等效面积[5]. 2009年, 他们用相同的方法提取了在轨标定球的等效温度和等效面积, 并与在轨标定球的物性进行比较, 从而验证了该方法的有效性[6]. 2014年, 他们又对死亡卫星和活跃卫星进行了研究, 根据等效温度和等效面积随时间的变化, 能很好地对它们进行区别[7,8]. 可以看出, 这些研究都是围绕目标的等效温度和等效面积开展. 但由于卫星红外光谱受卫星的表面温度、形状和探测器的方位等多种因素的影响, 以及对测量卫星的物性信息的缺乏, 外场实验无法解释等效温度和等效面积的物理意义, 并且无法根据卫星的特点进一步挖掘新特征. 通过准确地建立基于地基探测的卫星热红外光谱模型, 再对模拟的红外光谱数据进行分析, 是一种可行的研究方法.
本文在已有的卫星热红外光谱模型的基础上[9−15], 考虑了太阳辐射、地球辐射、地基探测器可探测卫星的范围、卫星各面对探测器的可见情况、大气衰减等因素的影响, 更加真实地建立了卫星热红外光谱模型. 以风云三号卫星为例, 在观测时序上仿真了地面探测器接收到该卫星的3—14
${\text{μ}}{\rm m}$ 红外光谱辐照度; 时序上分析了卫星温度场和卫星各面对探测器可见情况对红外光谱的影响; 用普朗克公式拟合该卫星的红外光谱数据, 得到等效温度和等效面积, 解释了它们的物理意义. 在卫星温度场分析的过程中, 发现帆板温度和本体温度有较大的温度差, 从而进一步实现了卫星本体和帆板的温度和面积的分离, 并实现新特征的提取. -
将卫星视为具有漫发射和漫反射的灰体, 其表面由若干温度互不相同的面元构成, 面元自身辐射由面元表面温度和发射率决定. 则面元自身的红外光谱辐射出射度表示为
式中c1为第一辐射常量, 值为3.742 × 10–16 W·m2; c2为第二辐射常量, 值为1.4388 × 10–2 m·K;
$\lambda $ 为波长,单位为${\text{μ}}{\rm m} $ ;Ti为面元表面温度, 由卫星温度场数值仿真得出 [12];${\varepsilon _i}$ 为面元的发射率. -
太阳通常可以视为一个温度为5900 K的黑体, 由普朗克公式可以计算出它的光谱辐射出射度:
式中Ts为太阳温度. 由于日地距离很远, 可以认为太阳在空间上的光谱辐射是均匀分布的, 则光谱辐射强度为
式中Rs为太阳半径, Rs = 6.9599 × 108 m; 卫星处太阳光谱辐照度为
式中Dse为日地距离, Dse = 1.4968 × 1011 m; 卫星某面元直接反射太阳辐射的光谱辐射出射度为
${\rho _{\rm{s}}}$ 为面元对太阳辐射的反射率;${\theta _{\rm{i}}}$ 为面元法线与太阳入射方向的夹角, 由坐标变换可得出[16], 若该角大于90°, 即太阳光不能照到面元上, 则$\cos {\theta _{\rm{i}}}$ = 0. -
地球辐射对于低轨道和中轨道卫星的辐射的影响是不可忽略的. 面元反射的地球辐射主要包括地球自身的红外辐射和地球反射太阳的红外辐射. 计算地球自身红外辐射时, 可以将地球视为一个温度为280 K的黑体. 由普朗克公式可得出地球的光谱辐射出射度:
式中Te为地球温度. 地球自身红外辐射在卫星处产生的光谱辐照度表示为
式中Re是地球的半径, Re = 6370 km, De是卫星离地面的高度. 地球反射辐射在卫星处产生的光谱辐照度表示为
${\rho _{{\rm{es}}}}$ 是地球对太阳辐射的平均反射率,${\rho _{{\rm{es}}}}$ = 0.35. 面元反射地球辐射的光谱辐射出射度表示为${\rho _{\rm{e}}}$ 是面元对地球自身辐射的反射率;${\theta _{\rm{e}}}$ 为面元中心指向地心方向与面元法线之间的夹角, 若该角大于90°, 地球辐射不能照到面元上, 则$\cos {\theta _{\rm{e}}}$ = 0. -
地面观察卫星的可见范围受仰角的限制, 地面观测点与卫星之间的视线方向在当地的仰角应大于5°. 但一些天文观测站由于其海拔高度较高, 在仰角为负的时候仍能观测到卫星. 例如, 在毛伊岛的空间观测系统, 其在海拔高度为3076 m的高山上, 能观察到最低仰角为–58°的空间目标[5]. 这里, 地基观测设备的最低仰角设为–30°. 地面站覆盖区是以地面观测点P为中心的可观区, 星下点在此圈内的卫星都是可观测的. 如图1所示, 该区是以P点为中心, 满足仰角E为给定值, 星下点B相对于P点的分布圈.
如图1, 从地面站观察卫星的仰角是在含观察点P、地心O和卫星S的平面内, 卫星视线方向与观察点P水平线之间的夹角E为仰角. 由几何关系可得, 在平面OPS内, 斜距
$\rho $ 和仰角E为:式中Re为地球半径, r为OS的距离,
$\psi $ 角为卫星星下点B与观察点P之间的地心夹角. 由球面三角形PNPB有式中L为P的地心纬度,
$\theta $ 为P相对于B的经度,$\varphi $ 为B的地心纬度. -
卫星在地基探测器的可见覆盖区域时, 卫星可被观察到. 卫星表面面元与探测器入瞳面的几何位置关系如图2所示. 其中, o,
$o'$ 分别为面元和探测器入瞳面的中心, ni为太阳光入射方向, nr为面元的法线方向, nd为探测器入瞳面的法线方向, ne为面元指向地心方向;${\theta _{\rm{i}}}$ 为太阳光入射方向与面元法线方向之间的夹角,${\theta _{\rm{r}}}$ 为面元法线与两面中心连线的夹角,${\theta _{\rm{d}}}$ 为探测器入瞳面法线与两面中心连线的夹角,${\theta _{\rm{e}}}$ 为面元指向地心方向与面元法线之间的夹角.卫星某面元的光谱辐射出射度可由(2), (6), (10)式得出
通常面元的面积相对于卫星与探测器之间的距离小很多, 因此可以将面元视为点源, 根据朗伯余弦定律可得面元在探测器
${\theta _{\rm{r}}}$ 方向的红外光谱辐射强度:式中Ai为面元的面积.
由面元与探测器入瞳面之间的几何关系, 探测立体角
$\varOmega $ 可写成:式中D为探测器入瞳面直径,
$\rho $ 为探测器到面元的距离即斜距. 由(14), (15)式可得出面元在探测器入瞳处红外光谱辐照度:若
${\theta _{\rm{r}}} + {\theta _{\rm{d}}}$ (面元法线方向与探测器入瞳法线反方向的夹角)大于90°, 即该面元对探测器不可见, 则Est,i = 0.卫星在探测器入瞳处红外光谱辐照度:
式中n表示卫星表面被分化成的面元数.
-
大气环境下, 地基探测器对空间目标进行红外辐射测量时, 目标辐射在到达探测器的传输过程中受到大气的衰减, 3—5
${\text{μ}}{\rm m}$ 和8—14${\text{μ}}{\rm m}$ 为两个红外大气窗口, 其中4.2—4.4${\text{μ}}{\rm m}$ 波段是一个极强的吸收带, 几乎全部吸收; 同时, 大气背景辐射也叠加到目标辐射上一同到达探测器. 则实际探测器入瞳处接收到的红外辐射:下标i表示第i个波段,
${\tau _i}$ 为第i个波段的大气透过率, Li表示第i个波段的大气背景辐射. 为了得到每个波段精确的卫星红外辐射Est(${\lambda _i}$ ), 需要进行背景扣除和大气衰减修正. 地基红外探测系统可对空间目标和周围天空背景进行测量, 目标与背景在入瞳上的辐照度分别为:相减得
目前工程上常采用扫描振镜的背景实时扣除方法和望远镜的快速斩波方法实现实时背景扣除, 其背景扣除误差在±1%之内[3,17]. 在式中还存在大气透过率
${\tau _i}$ 影响目标辐射, 利用红外标准星可以修正大气衰减, 得到准确的卫星红外辐射Est(${\lambda _i}$ )[18−20]. 目前, 红外标准星光谱数据库中有620颗星, 覆盖全天空[18]. 再对目标进行测量后, 转动光学系统将其对准离目标最近的红外标准星, 用上述背景扣除方法得到红外标准星的红外辐射:式中Eir(
$\lambda _i$ )为光谱数据库中的已知值, 由天基探测器测量得到, 不受大气影响. 可得第i个波段的大气透过率:该大气透过率的不确定度在3%以内[18]. 由此可得经大气修正后卫星在探测器入瞳处的红外辐照度:
其不确定度在4%以内.
-
卫星由许多不同面积和表面温度的部分组成. 因为卫星热红外光谱主要由卫星自身热发射引起, 所以可以忽略卫星反射的太阳和地球光谱. 将卫星的各个组成部分视为灰体. 则可得出它们在探测器处简化的红外光谱辐照度公式:
式中m表示提取卫星m个不同面积和表面温度的组成部分, 这里面积
$A_i'$ 为发射投影面积.将这个公式去拟合实际的红外光谱数据, 便是一个最优化问题. 通常的做法是用最小二乘的方法来构建目标函数:
式中N为波段数. 最小化目标函数J, 便可得到最优的面积和温度. 有很多最优化算法可以解决该问题, 这里可以采用简单且实用的拟牛顿方法.
-
对风云三号卫星于春分9:50至10:05时间段到达地面探测器的光谱辐照度进行了仿真计算, 并对光谱辐照度进行了分析.
-
我们将卫星简化为简单的几何结构模型, 如图3所示. 太阳帆板和卫星本体都视为规则的长方体. 其表面视为漫反射和漫发射的灰体. 对该模型建立坐标系, 以本体中心为坐标原点, z轴指向地心方向, x轴指向卫星运行方向, y轴由右手法则确定. 风云三号卫星 (FY-3) 的物性参数和轨道参数分别见表1和表2.
-
地面探测器与卫星降交点同经度, 纬度为北纬26°, 入瞳直径为1.2 m. 对风云三号卫星于春分时9:50到10:05期间(卫星在探测器的探测范围内)进行了仿真计算, 利用(10)和(11)式求出倾角和斜距, 如图4(a)所示. 卫星各表面法线方向与探测器方向的夹角余弦即
$\cos {\theta _{\rm{r}}}$ 如图4(b)所示, 风云三号卫星各时段下的表面稳定温度场如图4(c)所示. 由(13)和(17)式可计算求得卫星在探测器处的光谱辐照度, 并加入大气修正后4%的不确定度, 如图4(d)所示. -
当n = 1时, 用(25)式去拟合红外光谱数据得到的是卫星整体的等效温度和等效面积, 如图5所示; n = 2时, 用(25)式去拟合红外光谱数据得到两组不同的温度和面积, 如图6所示.
-
如图4(a)所示, 在观测期间卫星对地面探测器仰角和斜距变化非常大. 在短短的15 min里, 仰角的变化为–10°—70°; 斜距变化为1000—5000 km. 这导致探测器的信号强度会有几十倍的差距, 增大了地基探测器对极地卫星的观测难度. 地基探测器在对卫星观测期间, 由卫星运动引起的对探测器可见情况的改变是影响卫星红外光谱数据的主要因素. 如图4(b)和图4(c)所示, 在9:50—10:05期间, 卫星的温度场较为稳定, 而卫星各面对探测器的可见情况都有明显变化, 这导致了图4(d)中8—14
${\text{μ}}{\rm m}$ 的光谱曲线的差异. 对比图4(d)中模拟的风云三号卫星在探测器上的红外光谱辐照度和BASS系统实测的地球同步卫星红外光谱, 它们的曲线大体一致, 且数据的不确定度也大致相同, 这说明了模型的准确性.图5(a)显示等效温度很接近太阳帆板的温度, 温差仅在15 K左右; 这是由于在普朗克函数中温度在指数项上, 温度在幅值和曲线形状上都对卫星红外光谱有很大的影响, 对于帆板这种温度较高且有一定面积的部件在卫星红外光谱中占了主导地位. 图5(b)显示等效面积与卫星投影面积具有较明显的差距, 但可以看出等效面积在时序上的变化与卫星投影面积是一致的. 卫星中有的温度较低且有一定面积的部分, 由于它的温度与等效温度有较大差距, 导致计算的面积小于它的投影面积, 因此等效面积与卫星投影面积有较明显的差异. 而各低温部分的面积减小的比率大致相同, 因此等效面积能表征卫星投影面积的变化.
由卫星红外光谱反演卫星各组成部分的温度和面积等特征是一个不适定问题, 即反演的解不是稳定惟一的[3]. 从图4(c)卫星的温度场可以看出, 卫星本体和帆板的温差较大, 这使得它们的光谱差异也有所增大. 因此有望将它们分离, 分别得到帆板和本体的面积和温度. 通过设置合理的迭代初始值, 高温设为等效温度, 低温设为低于等效温度100 K, 它们的面积设为等效面积的一半; 由此可得到n = 2时稳定惟一的解. 结果如图6所示, 其拟合出的较高温度大致能和太阳帆板温度重合, 同时, 该温度对应的面积与帆板的面积也相一致, 整体的面积与卫星投影面积之间的差距也有所减小. 这有理由相信拟合得到的较低温度能表征卫星本体温度, 该温度对应的面积能表征本体面积.
-
本文考虑了地基探测器可探测卫星的范围、太阳光照、卫星各面对探测器的可见情况等因素的影响, 建立了一个更准确的基于地基探测的卫星红外光谱模型. 仿真结果发现, 对于风云三号这样的三轴稳定极地卫星, 仅有十几分钟的时间能被地基探测器探测到, 而且在这期间其对地面探测器仰角和斜距变化非常大. 在各种影响因素中, 由卫星运动引起的对探测器可见情况的改变是影响卫星红外光谱数据的主要因素. 对卫星红外光谱数据的进行时序分析, 解释了等效温度和等效面积的物理意义. 认识到卫星的等效温度更接近太阳帆板的温度, 等效面积能表征卫星投影面积的变化. 利用帆板温度和本体温度有较大的温度差, 进一步实现了卫星本体和帆板的温度和面积的分离. 从而实现了利用卫星红外光谱数据获得卫星具体信息的可能, 这对地基探测器对卫星的监测与识别有重要意义.
基于地基观测的时序卫星红外光谱建模与分析
Modeling and analyzing of time-resolved satellite infrared spectrum based on ground-based detector
-
摘要: 针对地基观测的卫星红外光谱受复杂因素的影响和外场试验对测量卫星物性信息的缺乏, 无法解释卫星红外光谱反演出特征的有效性和具体物理意义的问题, 提出了一种基于地基观测的卫星热红外光谱的建模和分析的方法. 首先, 考虑了太阳辐射、地球辐射、卫星各面对探测器的可见情况、地基探测器可探测卫星的范围、大气衰减等因素的影响, 更加准确地建立卫星热红外光谱模型. 然后, 以风云三号卫星为例, 利用该模型计算了在观测时序上卫星在地基探测器入瞳上的3—14
${\text{μ}}{\rm m}$ 红外光谱辐照度; 分析了影响卫星红外光谱变化的主要因素. 最后, 利用普朗克函数拟合卫星红外光谱, 提取出特征与卫星的物性比较, 并对其进行分析. 结果表明: 在各种影响因素中, 由卫星运动引起的对探测器可见情况的改变是影响卫星红外光谱数据的主要因素. 等效温度和等效面积物理含义能被有效地解释, 等效温度接近于太阳帆板的温度, 温差仅在15 K左右, 等效面积能表征卫星投影面积的变化; 发现利用帆板和本体有较大的温差, 能实现帆板和本体的分离, 并实现新特征的提取.Abstract: Satellite infrared spectra based on ground-based detector are affected by complex factors such as satellite surface temperature, solar radiation, observation angle, etc, whose change cannot be detected in external field experiment. Therefore, it is impossible to analyze what are the main factors that affect the satellite infrared spectra. At the same time, due to the lack of physical information about the satellites through the external field experiment, the validity and physical significance of retrieving features from satellite infrared spectrum cannot be explained. In view of the above problem, a method to model and analyze satellite thermal infrared spectra based on ground-based detector is proposed. It is a feasible research method to accurately establish satellite thermal infrared spectrum model based on the ground detection, then to analyze the simulated infrared spectrum data. Firstly, considering the solar radiation, earth radiation, detectable range of the satellite on the detector, observation angle, atmospheric attenuation, etc., the satellite thermal infrared spectrum model can be established more accurately. Then, taking FY-3 satellite for example, the physical and orbital parameters of the satellite are set up, and the 3−14${\text{μ}}{\rm m}$ infrared irradiance of the satellite on the pupil of the detector is calculated by using the model. Meanwhile, the main factors affecting the infrared spectrum of the satellite are analyzed. Finally, the equivalent temperature and equivalent area are extracted by fitting the satellite infrared spectrum with the Planck formula. And they are compared with the physical properties of the satellite. The results show that among the various factors, with the satellite’ movement, the change of the visible state of the satellite induced by the satellite’s movement is the main factor that affects the satellite infrared spectrum. The physical meanings of the equivalent temperature and equivalent area can also be explained effectively. The equivalent temperature is close to the temperatures of the solar panels, and their temperature difference is only about 15 K. The change of equivalent area is consistent with that of the satellite projected area. Moreover, it is also found that there is a large temperature difference between the solar panels and the body, which makes their infrared spectra obviously different. Therefore, it is hopeful to obtain the areas and temperatures of the solar panels and the body respectively. This research can make up for the shortcomings of the external field experiments and promote the monitoring and recognizing of satellites by ground-based infrared detectors.-
Key words:
- satellite infrared spectrum /
- ground-based detector /
- Planck formula /
- feature extraction .
-
-
图 4 在观测期间, (a) 风云三号卫星对地面探测器的倾角和斜距, (b)卫星各面法线与探测器连线的夹角余弦, (c)风云三号卫星的模拟温度场, (d)风云三号卫星在探测器上的红外光谱辐照度和BASS系统实测的地球同步卫星红外光谱[3]
Figure 4. During the observation period, (a) the elevation angle and range of the FY-3 satellite to ground-based detector, (b) the angle cosine between the normal of satellite’s side and the direction of detector, (c) the simulated temperature field of the FY-3 satellite, (d) the infrared spectral irradiance of the FY-3 satellite on the detector and the infrared spectral irradiance of geosynchronous satellite measured by BASS.
图 6 在观测期间, (a)较高温度和卫星帆板温度的比较, (b)较高温度对应的面积和帆板面积的比较, 整体面积和卫星对探测器投影面积的比较
Figure 6. During the observation period, (a) the comparison of the higher temperature of n = 2 and solar panel temperature, (b) the comparison of the area corresponding to higher temperature and the area of solar panel, the comparison of the sum of the areas of n = 2 and the projection area of the satellite to the detector.
表 1 风云三号卫星的物性参数
Table 1. Physical parameters of the FY-3 satellite.
部件名称 几何尺寸/mm 材料 发射率 吸收率 卫星本体 4460 × 2200 × 3790 F46聚酯薄膜 0.81 0.1 太阳帆板 40 × 7800 × 3790 背面 SR107白漆 0.87 0.17 正面 太阳电池 0.86 0.9 侧面 有机黑漆 0.88 0.93 表 2 风云三号卫星的轨道参数
Table 2. Orbital parameters of the FY-3 satellite.
轨道半长轴/km 偏心率 倾角/(°) 升交点赤经/(°) 降交点地方时 周期/min 7207 0.001 98.5 150 10:00 102 -
[1] 孙成明, 赵飞, 袁艳 2015 物理学报 64 034202 doi: 10.7498/aps.64.034202 Sun C M, Zhao F, Yuan Y 2015 Acta Phys. Sin. 64 034202 doi: 10.7498/aps.64.034202 [2] Fulcoly D, Kalamaroff K, Chun F 2012 J. Spacecr. Rockets 49 76 doi: 10.2514/1.A32002 [3] Lynch D, Russell R, Gutierrez D, et al. 2006 Proc. Adv. Maui Opt. Space Surveillance Technol. Conf. Hawaii, September 10–14, 2006 [4] 徐灿, 张雅声, 赵阳生, 李鹏 2017 光谱学与光谱分析 37 672 doi: http://www.gpxygpfx.com/CN/abstract/abstract8994.shtml Xu C, Zhang Y S, Zhao Y S, Li P 2017 Spectrosc. Spect. Anal. 37 672 doi: http://www.gpxygpfx.com/CN/abstract/abstract8994.shtml [5] Skinner M, Payne T, Russell R, et al. 2007 Proc. Adv. Maui Opt. Space Surveillance Technol. Conf. Hawaii, September 12–15, 2007 [6] Skinner M, Russell R, Rudy R, et al. 2009 Proceedings of the International Astronomical Congress 60th Meeting Daejeon, Republic of Korea, October 12–16, 2009 p1791 [7] Skinner M, Russell R, Kelecy T, et al. 2014 Acta Astronaut. 105 1 [8] Skinner M, Russell R, Kelecy T, et al. 2014 Proceedings of the International Astronomical Congress 65th Meeting Toronto, Canada, September 9, 2014 p1188 [9] 张伟清 2006 博士学位论文 (南京: 南京理工大学) Zhang W Q 2006 Ph.D. Dissertation (Nanjing: Nanjing University of Science and Technology)(in Chinese) [10] 张永阳 2007 硕士学位论文 (南京: 南京理工大学) Zhang Y Y 2007 M.S. Thesis (Nanjing: Nanjing University of Science and Technology)(in Chinese) [11] 马伟 2011 博士学位论文 (南京: 南京理工大学) Ma W 2011 Ph.D. Dissertation (Nanjing: Nanjing University of Science and Technology) (in Chinese) [12] 孙成明, 袁艳, 张修宝 2010 物理学报 59 7523 doi: 10.7498/aps.59.7523 Sun C M, Yuan Y, Zhang X B 2010 Acta Phys. Sin. 59 7523 doi: 10.7498/aps.59.7523 [13] 王雨飞, 李强, 廖胜 2011 红外与激光工程 40 2085 doi: 10.3969/j.issn.1007-2276.2011.11.004 Wang Y, Li Q, Liao S 2011 Infrar. Laser Eng. 40 2085 doi: 10.3969/j.issn.1007-2276.2011.11.004 [14] 汪洪源, 陈赟 2016 红外与激光工程 45 504002 doi: http://journal02.magtech.org.cn/Jwk_irla/CN/Y2016/V45/I5/504002 Wang H Y, Chen Y 2016 Infrar. Laser Eng. 45 504002 doi: http://journal02.magtech.org.cn/Jwk_irla/CN/Y2016/V45/I5/504002 [15] 李文豪, 刘朝晖, 穆猷, 等 2017 红外与激光工程 46 604003 Li W H, Liu Z H, Mu Y, et al. 2017 Infrar. Laser Eng. 46 604003 [16] 谈和平, 夏新林, 刘林华, 阮立明 2006 红外辐射特性与传输的数值计(哈尔滨: 哈尔滨工业大学出版社) 第 378页 Tan H P, Xia X L, Liu L H, Ruan L M 2006 Numerical Calculation of Infrared Radiation Characteristics and Transmission (Harbin: Harbin Institute of Technology Press) p378 (in Chinese) [17] 王先起, 廖胜, 沈忙作, 黄建明 2005 光电工程 32 5 doi: 10.3969/j.issn.1003-501X.2005.06.002 Wang X Q, Liao S, Shen M Z, Huang J M 2005 Opto-Electronic Engineering 32 5 doi: 10.3969/j.issn.1003-501X.2005.06.002 [18] 殷丽梅, 刘俊池, 王建立, 等 2014 光子学报 43 1204004 doi: http://www.photon.ac.cn/CN/abstract/abstract20550.shtml Yin L M, Liu J C, Wang J L, et al. 2014 Acta Phot. Sin. 43 1204004 doi: http://www.photon.ac.cn/CN/abstract/abstract20550.shtml [19] 刘莹奇, 刘祥意 2014 光学学报 34 0512003 doi: http://www.opticsjournal.net/Articles/abstract?aid=OJ140418000164FbHeKh Liu Y, Liu X 2014 Acta Opt. Sin. 34 0512003 doi: http://www.opticsjournal.net/Articles/abstract?aid=OJ140418000164FbHeKh [20] Skinner M, Russell R, Kelecy T, et al. 2012 Acta Astronaut. 80 154 doi: 10.1016/j.actaastro.2012.04.044 -