中国计划执行“ 天问二号” 探测器任务, 通过一次发射实现一颗近地小行星(2016HO3)的取样返回和一颗主带彗星(133P)绕飞探测。 其中的热辐射光谱仪用于探测小行星表面温度, 以研究小天体热辐射物理参数。 热物理模型是求解小行星热物理参数, 如表面反照率、 热惯量、 表面粗糙度以及有效直径等的常用手段。 热物理模型需要考虑小行星的反照率、 热传导率、 比热容、 热发射率、 密度、 粗糙度等诸多因素, 用于模拟行星发射出的红外辐射, 并与观测到的红外辐射比较。 而获取这些输入物理量, 必须先获得行星的表面温度分布。

美国国家航空航天局(NASA)在2016年成功发射了OSIRIS(欧西里斯)探测器, 用于观测近地小行星“ 贝努” (bennu), 并采样回收。 2018年12月, “ 欧西里斯” 探测器进入环绕贝努小行星的轨道, 开始探测和采样任务。 “ 欧西里斯” 探测计划具有深刻的科学意义: 贝努小行星与地球之间最小轨道交会距离为0.002 AU, 具有潜在碰撞危险。 “ 欧西里斯” 可以对贝努进行全球温度绘图, 研究Yarkovsky和YORP效应, 预测贝努未来轨道位置。 对小行星贝努风化层表面物质采样, 研究小行星矿物质和有机物组成。 对贝努的矿物质分布和化学特性进行全球绘图, 研究其地质特征^{[1, 2, 3]}。

“ 欧西里斯-雷克斯” 热辐射光谱仪(OTES)对贝努进行热物理性质和矿物学测量, 得到全球热辐射光谱数据。 OTES是一种傅里叶变换光谱仪, 具有1 750~100 cm^-1的光谱范围, 光谱分辨率为8.66 cm^-1, 视场角6.5 mrad^[4]。 OTES的全球绘图具有40 m的空间分辨率, 局部光谱信息采集的空间分辨率为大约4 m。 使用OTES的热红外数据, 可以反演贝努的表面温度, 这对确定最佳采样区域和明确贝努小行星的热惯量和矿物质组成有重要意义^{[5, 6, 7, 8, 9]}。

基于“ 欧西里斯-雷克斯” 热辐射光谱仪(OTES)获取的高光谱数据, 采用温度-比辐射率分离算法(temperature-emissivity separation algorithm, TES), 计算贝努表面温度^[10]。

热红外遥感反演的关键在于温度与比辐射率分离。 传感器接收到的辐亮度表达式如式(1)

$I={{\varepsilon }_{i}}\frac{2h{{c}^{2}}{{\nu }^{3}}}{{{\text{e}}^{\frac{hcv}{kT}}}-1}$(1)

式(1)中, ε _i是行星物质在第i个波数处的比辐射率, ν 是电磁波波数, h是普朗克常数, k是玻尔兹曼常数。

表面温度和发射率是耦合的, 从辐射测量中求解它们是一个欠定问题, 即用N个方程(每个波段的N个辐射值)求解N+1个未知数(一个温度和N个发射率值)。 因此必须施加一些约束, 以减少需要估计的未知数的数量。 实现温度和比辐射率分离的代表性算法主要有: 光谱平滑迭代法、 温度-比辐射率分离算法、 发射率归一化法(emissivity normalization method, NEM)、 归一化发射率最大最小值差(maximum-minimum normalized emissivity difference, MMD)法、 比值法(spectral ratio method)、 α 残差法、 包络线法等^{[11, 12, 13, 14]}。 这些算法通过多波段或高光谱的热红外观测, 使用对发射率或温度无关波谱指数(temperature-independent thermal infrared spectral indices, TISI)或近似灰体假设等目标物质光谱先验约束条件, 以求解欠定方程。 这些算法都可得到表面温度, 但存在着各自的优点和缺点。

比如热对数残差模型法, 使用维恩公式近似普朗克定律, 通过对数运算分离温度和发射率, 使用辐射度值、 辐射度值对像素或波长的平均值、 测量的辐射度值对像素和波长的平均值, 求解辐射率。 它不可避免地引入了维恩公式近似普朗克定律所带来的误差, 而且误差随着温度的增加而增大。 发射率归一化法, 始终采用单一的发射率值(如0.97)进行表面温度的反演, 忽略了表面的空间变化, 降低了反演的精度。 比值法利用在同样的温度误差下, 相对发射率的精度远远大于发射率的精度这一特性, 对相对发射率最高的波段假设发射率为1, 计算表面温度。 常用其他方法计算表面温度后, 再用比值法进一步精确发射率和温度。 反演的精度受到实际发射率最大值的影响, 当实际发射率最大值远远小于1时, 算法误差不可忽略。 α 残差模型法是热对数残差模型的简化, 只使用了辐射度值、 测量的辐射度值对波长的平均值来构造只包含辐射度值的辐射率函数。 α 残差模型简化了热对数残差模型, 降低了反演的精度, 但提高了运算的效率。 光谱平滑迭代法, 在一个参考温度附近设置一系列等间隔温度, 计算每一个温度节点对应的地物发射率曲线的平滑度指数。 当平滑度指数最小时, 发射率曲线最平滑, 对应温度为最终温度结果。 因为许多物质的发射率在波长大于10 μ m时, 随波长缓慢变化, 所以灰体发射率法假设在波长大于10 μ m时目标发射率近似不变。 该方法已经被证实适用于多种天然材料。 温度独立光谱指数法利用表面温度的N次数值接近普朗克辐射亮度值, 先假设一个初始的发射率, 然后利用迭代接近表面真实温度。 这种近似关系同样引入了误差。 一种改进的发射率归一化法(adjusted emissivity normalization method, ANEM)利用土壤和植被的发射率相差较大, 通过归一化植被指数(normalized difference vegetation index, NDVI), 构建发射率最大值关于NDVI的函数式。 实际计算时, 先求解NDVI来判断像素中植被和土地的比例, 代入函数式中求发射率最大值, 从而计算温度。 这种方法考虑了不同材料在混合像元中对发射率的影响, 算法精度依赖于NDVI的精确度, 并且需要预先计算或已知表面的NDVI分布。

其中, 温度-比辐射率分离算法综合了发射率归一化法(NEM)、 MMD法、 比值法(RAT)三种算法的优点, 是目前公认精度较高的算法。 温度-比辐射率分离算法已经应用于国外多种传感器上, 比如ASTER和MODIS传感器, 均证实有良好的温度反演精度。

温度-比辐射率分离算法分为以下步骤^[15]:

(1)NEM模块

发射率归一化法的原理在于, 目标物质在某一波段的发射率达到最大值(ε _max), 接近于1, 并且发射率的变化很小。 认为物质在此波段发射率为常数(ε _max), 计算温度

${{T}_{i}}=\frac{hcv}{k\text{ln}\left( \frac{2\ \ {{\varepsilon }_{i}}\ \ {{c}^{2}}\ \ h\ \ {{\nu }^{3}}}{I}+1 \right)}$(2)

εi=εmax(3)

其中, T_i表示第i个波长通道反演的表面温度。

发射率归一化法, 采用单一的发射率值(如0.97)反演表面温度, 忽略了目标的空间变化。 同时, 目标的实际发射率最大值和该单一发射率值之间的差异, 是发射率归一化法的主要误差来源。 目标的实际发射率最大值应尽可能地包含在反演波段内。 近地小行星贝努每6年靠近地球一次, 在地望远镜观测的可见光和近红外数据表明, 贝努是一颗B类小行星。 B类小行星的矿物质构成以硅酸盐为主, 包括长石、 辉石、 橄榄石等。 粉末状的硅酸盐和陨石在实验室测量时表现出独特的光谱特征, 比如克里斯琴森特征(Christiansen feature, CF)^[16]。 克里斯琴森特征是指硅酸盐在7.5~10 μ m波段范围内, 存在某一波长, 这个波长的发射率最大、 反射率最小。 这一波长称为该硅酸盐的CF位置。 纯净矿物的CF位置取决于矿物的聚合程度, 例如长石的CF位置比辉石、 橄榄石的CF位置向短波移动。 硅酸盐的粒径大小主要影响了克里斯琴森特征的光谱对比度。 因此, 可以认为, 贝努表面物质在7.5~10 μ m波段范围内发射率达到最大值, 并使用发射率归一化法进行反演, 初步得到发射率谱和温度。

对初步得到的温度, 选择温度最大值, 认为是NEM模块反演的温度结果, 计算各波长通道的发射率。

(2)比值模块

使用比值法计算相对发射率, 保留了发射率的形状。 温度-比辐射率分离算法通过比值模块来缩小反演温度和发射率的误差。

相对发射率β _i表达式为

${{\beta }_{i}}=\frac{{{\varepsilon }_{i}}}{\frac{1}{N}\overset{N}{\mathop{\mathop{\sum }_{j=1}}}\,{{\varepsilon }_{j}}}$(4)

式(4)中, ε _i是NEM模块计算的发射率, N是波长通道数。

比值法的核心在于, 在同样的温度误差下, 相对发射率的精度远远大于发射率的精度。 因此, 先计算相对发射率β _i, 再通过MMD模块还原真实发射率ε _i。 以下, 证明在同样的温度误差下, 相对发射率的精度远远大于发射率的精度^[13]:

假设第i个波长通道的辐亮度为I_i, 波长为λ _i; 第i+1个波长通道的辐亮度为I_i+1, 波长为λ _i+1。 则两个波长通道的相对发射率H_i满足

Hi=Ii+1Ii(5)

温度真值为T, 而计算的温度为< T> , 导致相对发射率H_i的误差为δ , 则泰勒展开有:

Hi=Hi(1+δ)(6)

δ≤(T-< T> )∂Hi∂T(7)

使用维恩位移定律近似普朗克定律, 相对发射率H_i对温度T求导 ∂Hi∂T为

∂Hi∂T≈-HiC2(λi+1-λi)(< T> λi)2(8)

C2=hck(9)

(3)MMD模块

归一化发射率最大最小值差(maximum-minimum normalized emissivity difference, MMD)表示为

MMD=max(βi)-min(βi)(10)

TES算法在发射率最小值(ε _min)和MMD之间建立了一种函数关系, 来还原发射率

εmin=a+b·MMDc(11)

εi=βiβminεmin(12)

其中, β _min是β _i的最小值。

MMD模块基于天然或人造材料普遍具有一种规律: 当材料的发射率平均值越大时, 发射率光谱越平坦, 变化程度小, 即MMD值小; 反之, 当材料的发射率平均值越小时, 发射率光谱越波动, 变化程度大, 即MMD值大。 对于不同类型的目标(水、 岩石、 森林等)或不同类型的矿物质(硅酸盐、 硫化物等), 该函数关系的系数a、 b、 c不同。 这种函数关系代表了一种目标发射率特性的先验知识。

因此, 鉴于小行星的主要物质成分为硅酸盐, 在使用算法对行星“ 贝努” 进行表面温度反演之前, 必须根据实验室硅酸盐发射率特性进行MMD模块函数关系的拟合, 得到系数a、 b、 c。 之后系数a、 b、 c作为先验知识, 在算法迭代中保持不变。

(4)普朗克定律计算温度

${{T}_{i}}=\frac{hc\nu }{k\text{ln}\left( \frac{2\ \ {{\varepsilon }_{i}}\ \ {{c}^{2}}\ \ h\ \ {{\nu }^{3}}}{I}+1 \right)}$(13)

Ts=max(Ti)(14)

取计算的温度中的最大值(T_s)作为每次迭代的反演结果。 根据本次迭代得到的温度结果T_s, 计算发射率ε _i, 作为下一次迭代的输入

${{\varepsilon }_{i}}=I\frac{1}{2h\ \ {{c}^{2}}\ \ {{\nu }^{3}}}\left( {{\text{e}}^{\frac{hc\nu }{k{{T}_{s}}}}}-1 \right)$(15)

重复迭代, 直到满足迭代终止条件。 通常, 设置迭代终止条件为迭代到一定次数或相邻两次反演的温度之差满足要求。

选取2019年2月22日19时19分开始的序列, 进行温度反演。 取OTES官方产品为真值, 去除OTES产品中质量不佳、 太阳天顶角过大的样本点, 保留纬度在60° S— 60° N之间的点。 迭代终止条件为两次迭代结果的最大温差不超过1.0 K。 图6给出了OTES温度产品(横轴)和本算法反演的温度(纵轴)之间比较。 其中, 蓝色散点的横坐标和纵坐标分别表示样本点的OTES给出的温度值和本算法计算的温度值。 绿色拟合线是对OTES温度产品和本算法反演的温度结果进行最小二乘法进行函数拟合, 得到的函数关系, 非常接近于y=x。 这说明, 本算法反演的温度结果和OTES温度产品具有良好的一致性。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 OTES温度产品(横轴)和本算法反演的温度(纵轴)Fig.6 OTES temperature product (horizontal axis) and the temperature retrieved by the algorithm proposed in this paper (vertical axis)

图6中, 散点表示样本点的OTES给出的温度值和本算法计算的温度值。 绿色直线表示OTES给出的温度值和本算法计算的温度值拟合的函数关系, 用于衡量二者的一致性。

存在几个误差较大点, 偏离绿色拟合线, 表示在这几个样本点上, 本算法结果与OTES产品有较大误差。 分析这些样本点的测量辐射值和OTES记录的仪器光轴与行星表面相交情况, 分析原因可能是介于行星观测和深空定标之间, 行星未充满视场导致。 即, 这些样本点实际是深空和贝努表面构成的混合像元, 深空可以看成是真实温度为4 K的黑体。 反演误差(误差表示为: 误差平均值± 误差标准差)为(-0.366 0± 1.039 3) K。 可以看出, 本工作的温度-比辐射率分离算法具有良好的精度和普适性, 适用于小行星的温度反演。 中国“ 天问二号” 计划将抵达近地小行星2016HO3, 使用热辐射光谱仪进行温度和辐射率的全球测绘, 研究具有良好的精度和普适性的温度-比辐射率反演算法有重要意义。

使用温度-比辐射率分离算法, 进行了算法敏感性验证, 结果表明: (1)设置温度范围为115~415 K, 步长5 K, 随着温度的增加, 反演温度的均方根误差增加, 而发射率的误差大致不变。 (2)样本温度为295 K, 算法使用的波段起点为7.5 μ m, 采样间隔0.04 μ m, 终点为10~13.8 μ m, 步长0.2 μ m。 发现波长范围在7.5~13.8 μ m, 算法的MMD模块的精度最高。

使用温度-比辐射率分离算法, 对OTES辐亮度数据进行了温度和发射率的计算, 结果表明: 在剔除了大太阳高度角和高纬像素的情况下, 计算的温度的误差平均值为-0.366 0 K, 误差的标准差为1.039 3 K。 中国“ 天问二号” 在不远的将来将抵达近地小行星2016HO3, 该算法可应用于其搭载的热辐射光谱仪, 对温度进行反演。