1.   引　言

在卫星遥感器的发射过程中，周围环境的剧烈变化将导致其物理性能发生改变，在轨后，随着时间的推移其性能也会出现不同程度的衰减，因此需要对遥感器重新进行在轨定标以保证其观测数据的科学准确性^[1]。

目前，太阳反射波段（RSB）卫星载荷的在轨定标精度在10⁻²量级^[2]，与观测气候变化所需的10⁻³量级的精度水平还存在一定差距^[3]，考虑到进一步提高现有定标方法精度的难度较大，各国均将发展空基溯源标准计划列入日程，通过在太空建立辐射基准并采用交叉定标的方法将辐射基准进行传递，从而提高在轨定标精度水平。国际上，已经有NASA主导的CLARREO计划和NPL主导的THRUTH计划^[4-6]。我国近年来也开展了国产高精度在轨辐射基准的相关研制计划^[4-8]。

交叉定标方法是将在轨辐射基准传递至其他卫星的重要方法之一，是一种利用已有的定标精度较高的卫星作为参考对其他卫星进行定标的方法，其不受定标场景的限制，定标成本低，定标频次高^[1,9]。目前国际上关于交叉定标的研究已经取得了诸多成果^[1,9]，针对红外波段（IR）已经建立起基于METOP IASI和AQUA AIRS的交叉定标系统^[9]，在RSB波段普遍以定标手段丰富且定标精度较高的AQUA MODIS作为辐射基准，传递至其他卫星，但AQUA MODIS已在轨运行了16年，目前，S-NPP VIIRS已正成为RSB波段交叉定标基准的新选择^[10]。综合分析可知目前RSB波段交叉定标能够达到的平均精度水平在5%左右^[1,2,11]。

一般的交叉定标方法是将经过近一致观测筛选匹配后的数据点对代入回归模型中，然后使用普通最小二乘法（Ordinary Least Square，OLS）得出定标系数。而根据实际的交叉定标流程知，各项匹配环节对基准辐射值的影响有差异，使得各数据点对的质量不同，此时，若继续使用OLS回归定标系数将会降低结果的有效性^[12-13]。针对以上问题，本文提出一种基于不确定度分析的交叉定标改进方法，利用加权最小二乘法（Weighted Least Square，WLS）回归定标系数，从而提高定标结果的准确性。

2.   方　法

本文将RSB高光谱成像仪作为参考仪器对宽通道成像仪进行交叉定标，分析交叉定标过程中的辐射误差传递源，通过蒙特卡罗法量化计算各数据点对中辐射基准值的不确定度，从而解算出权重系数矩阵，利用WLS方法回归定标系数，并通过与待定标仪器官方定标系数进行比对来检验WLS方法的有效性。

2.1   辐射基准替代载荷的选取

因现阶段的基准载荷尚在研制中，缺乏可用数据，为提供具体的算法实现以及验证本文所提出WLS方法的有效性，本文通过比较现有的常见RSB波段高光谱载荷参数（表1），拟选用与基准载荷预定指标（光谱分辨率优于10 nm、光谱范围在380~2 350 nm内）最接近的EO-1卫星上搭载的HYPERION作为参考载荷，实现具体的定标算法。

表  1  与基准载荷指标相近的高光谱成像仪参数

Table  1.  Hyperspectral imager′s parameters similar to the radiance standard remote sensor index

载荷名称	光谱覆盖 范围/nm	光谱分辨率/ nm	空间分辨率/ m	幅宽/ km	通道数
HICO	380~960	5.7	100	—	—
HYPERION	400~2 500	10	30	7.7	242
LAC	930~1 580	2~6	250	185	256
HIS	450~950	4	100	50	128

下载: 导出CSV 

| 显示表格

由于TERRA MODIS与HYPERION的轨道交点较多，易于满足交叉定标的时空匹配条件，且其定标精度高，定标系数可以用于验证本文算法，同时产品数据易获取。因此，本文利用HYPERION作为参考仪器对TERRA MODIS 1~7通道进行定标。选取2012年7月4日起连续10天内HYPERION和MODIS的观测数据作为基础数据集^[14-15]，通过交叉定标各匹配环节获取满足近一致观测的数据点对，从而建立起MODIS 1~7通道的回归模型。

2.2   定标回归模型

交叉定标流程如图1所示，通过时间、空间、观测几何、光谱匹配环节选取多个满足近一致观测的数据点对，建立起回归模型，通过线性回归获得待定标仪器的定标系数。值得注意的是，在使用高光谱遥感器作为参考仪器对宽通道仪器进行定标时，通常直接采用光谱卷积的方式对其他宽通道式成像仪进行光谱匹配。

图  1  交叉定标流程图

Figure  1.  Flowchart of cross-calibration approach

下载: 全尺寸图片 幻灯片

设筛选后满足近一致观测的数据点对共N个，将来自参考仪器和待定标仪器的数据分别用N×1的列向量L_ref和DN_cal表示，$\left({{L_{{\rm ref}\;i}},D{N_{{\rm cal}\;i}}} \right)$代表一个数据点对，回归模型可用矩阵表示为：

其中，β₀和β₁为定标系数，将定标系数用2×1的列向量β表示。

一般采用OLS对定标系数矩阵进行求解，得出待定标遥感器的定标系数。OLS通过寻找系数矩阵使得残差平方和最小来对模型参数进行估计^[13]，可以得出β的解为：

考虑交叉定标回归模型，由于定标过程中各项匹配环节对不同数据点对中基准辐射值的影响不同，使得获取的各数据点对中辐射基准值的不确定度存在差异，即数据点对的质量存在差异，此时若继续使用OLS来回归定标系数将会降低定标系数的有效性，影响定标精度，此外，由于交叉定标时各数据点对的获取过程相互独立，各数据点对相互之间没有影响，因此本文采用WLS法来计算定标系数，提高定标系数的有效性。

2.3   加权最小二乘法

加权最小二乘法通过给获取的每个交叉定标数据点对加上权重系数提高回归结果的有效性，该权重系数可以反映该定标数据点对的质量以及对定标曲线回归的影响大小，数据点对对应的权重系数一般用数据点对中自变量不确定度平方的倒数表示^[12-13]。结合交叉定标回归模型，权重系数用数据点对中辐射基准值的不确定度平方的倒数表示，即获取的数据点对中辐射基准值的不确定度越大，说明数据点对的质量越差，对应权重系数越小，则该数据点对对定标系数回归的影响就越小，第i个定标数据点对对应的权重系数w_i可表示为$1/{u^2}\left({{L_{{\rm ref}\;i}}} \right)$。WLS法通过最小化广义残差平方和求出模型系数^[12-13]，可表示为：

式中，W为权重系数矩阵，结合定标模型，权重系数矩阵可以表示为：

在交叉定标时，要使用WLS计算定标系数，则需要计算定标数据点对中辐射基准值的不确定度，具体的计算方法如下文。

2.4   交叉定标误差传递源分析以及不确定度量化方法

在交叉定标数据点对的筛选过程中，L_ref的误差主要来源于以下几点：参考仪器自身观测的不确定度、时相匹配中过境时间的差异、空间匹配时空间分辨率不一致和空间位置匹配不当产生的误差、地表方向性即观测几何匹配不当引起的误差、光谱匹配时光谱响应函数的衰减和波长漂移等引起的误差^[2,16]。

在待定标仪器单个像元观测视场内，将与待定标仪器满足时间，空间位置，且观测角度一致的像元输出值设为L₀，经过近一致观测后配准的数据点对中辐射基准值表示为L_ref，并将各环节不确定度对最终获取的辐射基准值的影响看作乘性影响，因此，可在L₀基础上添加各项乘性因子（参考仪器自身K_ref、时间匹配K_time、空间匹配K_space、观测几何匹配K_angle、光谱匹配K_spectral），则L_ref表示为：

则根据不确定度的表达式^[17]，辐射基准值的不确定度可以表示为：

式（6）中，$u\left({{K_{\rm ref}}} \right)$、$u\left({{K_{\rm time}}} \right)$、$u\left({{K_{\rm space}}} \right)$、$u\left({{K_{\rm angle}}} \right)$、$u\left({{K_{\rm spectral}}} \right)$为对应各项K因子的不确定度，采用蒙特卡罗法（MCM）计算各K因子的不确定度。MCM法在原有输入数据的基础上经过添加扰动项进行重复试验得出输出量的概率分布后，给定置信概率，再根据概率分布进行统计得出置信区间，从而得出输出量的不确定度^[17]。MCM法中各项K因子不确定度的具体计算方法描述如下。

参考仪器自身观测不确定度引起的K_ref不确定度，采用MCM法结合参考仪器自身观测值的不确定度，在参考仪器观测值的基础上加入适当的扰动，输入交叉定标流程，重复试验，最终得出K_ref的分布，给定置信率，从而得出K_ref的不确定度$u\left({{K_{\rm ref}}} \right)$。

对于时相匹配不当引起的K_time因子的变化，主要考虑当地气象条件变化的影响，利用大气辐射传输模型，结合过境时间差内当地气象参数的变化，利用MCM法得出K_time的分布，结合置信率给出置信区间，从而得出其不确定度$u\left({{K_{\rm time}}} \right)$。

空间位置配准不当则采用滑动窗口法^[16]进行分析，使得两卫星图像在原有的配准基础上发生相对位移，分析辐亮度值的变化，利用MCM法得出K_space分布及其不确定度$u\left({{K_{\rm space}}} \right)$。

观测几何匹配不当主要受地表方向的影响，结合当地的BRDF模型和实际的观测几何角度的变化，输入辐射传输模型，利用MCM方法得出K_angle的分布及其不确定度$u\left({{K_{\rm angle}}} \right)$。

光谱匹配主要考虑遥感器光谱响应的变化，通过在已知光谱响应函数的基础上适当地加入扰动，得出K_spectral的分布及其不确定度$u\left({{K_{\rm spectral}}} \right)$。

在得出各项K因子的不确定度后，利用公式（6）得到辐射基准值的合成不确定度，针对每个数据点对都计算其辐射基准值的不确定度，从而得出权重系数矩阵作为WLS法回归模型系数。

3.   验证及比较分析

以HYPERION作为参考仪器对TERRA MODIS的1~7通道进行定标，验证WLS方法的有效性。为得出WLS法的权重系数矩阵，需要量化定标数据点对中辐射基准值的不确定度。以2012年7月10日卫星过境Libya4定标场地经纬度坐标（28.55°N，23.67°E）对应的数据点对为例，展示单个数据点对中辐射基准值的不确定度的获取过程。图2（左）、2（右）分别为HYPERION和MODIS过境Liyba4场地时的影像。

图  2  2012年7月10日HYPERION（左）和MODIS（右）过境Liyba4场地时的影像

Figure  2.  Images of HYPERION(left) and MODIS(right) when they passed the Liyba4 site on 10th July, 2012

下载: 全尺寸图片 幻灯片

在前文的不确定度分析方法中，关于时间和观测几何匹配引入的不确定度计算需要获取地面配合测量的气象参数，并将其输入至大气传输模型，本文通过选取晴朗无云场景的观测图像以及设置尽量小的过境时间差和观测几何角度差阈值，减小时相匹配和观测几何匹配对数据点对中基准辐射值的影响。对于所有的数据点对中的辐射基准值，忽略了时相和观测几何匹配引入的误差，使K_time和K_angle的不确定度近似为零。

已知HYPERION自身观测相对不确定度为3.4%^[11]；认为MODIS和HYPERION的空间位置匹配误差不超过单个MODIS像元大小，即空间匹配的误差在1 km内；MODIS在轨光谱响应的波长漂移在−2~2 nm范围内，响应度变化在5%范围内^[18]。（对获取的所有数据点对均使用以上条件）采用蒙特卡罗法分别计算K_ref、K_space和K_spectral的分布，求出置信率为95%时对应的置信区间，求得各因子对应的不确定度，忽略时间和观测几何引入的不确定度，计算合成不确定度。针对MODIS CH01通道Libya4定标场获取的数据点对的辐射基准值，计算其不确定度，获得的各因子的分布如图3（a）、3（b）、3（c）所示。

图  3  MCM法各项K因子概率分布图

Figure  3.  The probability density distribution of K-factors of each term of MCM

下载: 全尺寸图片 幻灯片

在该数据点对不确定度的量化过程中，由空间匹配不当引入的不确定度占比最大，K_space的相对不确定度为1.87%，其次是遥感器自身观测精度的影响，K_ref的相对不确定度为1.58%，由光谱匹配不当引入的不确定度占比最小，K_spectral的相对不确定度为0.067%。图3（c）可以看出，由光谱匹配不当引起的K_spectral概率分布呈偏斜状，这是由MODIS CH1的光谱响应函数和HYPERION观测输出的高光谱辐亮度曲线线型共同作用导致的。针对不同场景获得的高光谱辐亮度曲线不同，在光谱匹配不确定度分析时获得的K_spectral概率分布图也不尽相同。最终计算得出MODIS CH01通道在Libya4定标场获取的数据点对中辐射基准值总的相对不确定度为2.45%。

对所有的数据点对重复上述计算过程，得出每个数据点对中辐射基准值的不确定度，从而得出权重系数矩阵，利用WLS和OLS分别计算HYPERION对MODIS CH 1~7通道的定标系数β_WLS和β_OLS，并与MODIS产品数据中的定标系数β进行比对。为评价两种方法计算出的定标系数的定标效果，本文利用WLS和OLS定标系数对同一DN值序列进行定标并与官方定标系数的定标结果进行比对，并分别计算定标结果与官方定标结果的平均相对误差ε以及均方根误差RMSE。设DN值序列共包含N个点，对DN值序列中第i个DN值采用官方系数定标后的辐亮度值为${L_{0\;i}}$，经计算，获取的系数定标后的辐亮度值为L_i，则相对差异计算公式为：

平均相对误差计算公式为：

均方根误差计算公式为：

以HYPERION对MODIS CH01通道定标为例，所采用的WLS和OLS的回归曲线如图4（彩图见期刊电子版）所示，WLS和OLS计算的定标系数分别为（−0.325 3，0.026 1）、（−1.902 5，0.026 1），官方定标系数为（0.0，0.027 2），可以看出WLS的定标系数要更接近官方定标系数，同时可以发现两种方法回归得到的定标系数的差异主要体现在截距β₀上。将WLS和OLS得到的定标系数对同一DN值序列进行定标并计算其与官方定标结果的相对差异，如图5所示。从两种回归方法获得的定标系数与官方定标系数差异结果来看，采用WLS计算的定标系数的定标相对误差在低端明显减少，最大由9.12%降低到4.86%，平均相对误差由5.89%降低到4.30%，均方根误差由6.905 W/(m²·sr·μm)降低至5.455 W/(m²·sr·μm)，说明WLS方法有效提高了该通道的定标曲线拟合精度。

图  4  HYPERION对MODIS CH01定标数据点对以及OLS和WLS的回归定标线

Figure  4.  The data points of the MODIS CH01 and regression results of WLS and OLS

下载: 全尺寸图片 幻灯片

图  5  WLS和OLS计算的回归系数对同一DN值序列的定标结果与官方定标结果的相对误差

Figure  5.  The relative error between the regressed coefficients by WLS and OLS calibrating a DN sequence with the official calibrated coefficients

下载: 全尺寸图片 幻灯片

HYPERION对MODIS CH1~7通道的WLS和OLS计算定标系数结果为β_WLS和β_OLS、官方定标系数β、定标评价参数最大相对误差${\epsilon_{{\rm{max}}}}$、平均相对误差ε以及均方根误差RMSE的计算结果如表2所示。

表  2  WLS和OLS对MODIS CH1~7通道的定标系数回归结果

Table  2.  The calibration result of MODIS CH1~7 by WLS and OLS

MODIS CH	1	2	3	4	5	6	7
β0	0.000 0	0.000 0	0.000 0	0.000 0	0.000 0	0.000 0	0.000 0
β1	0.027 2	0.010 2	0.036 1	0.024 9	0.005 9	0.002 7	0.000 8
βOLS 0	−1.902 5	−1.951 6	−4.023 8	−0.900 5	−2.885 8	−0.746 6	−0.496 0
βOLS 1	0.026 1	0.009 6	0.031 8	0.023 4	0.005 3	0.002 4	0.000 8
βWLS 0	−0.325 3	−0.288 3	−5.891 5	−2.501 9	−2.379 2	−0.730 2	−0.364 4
βWLS 1	0.026 1	0.009 6	0.032 5	0.024 0	0.005 3	0.002 4	0.000 8
${{{\epsilon}}_{{{\max}}\;{\rm{OLS}}}}$	0.091 2	0.110 1	0.215 4	0.062 9	0.213 1	0.142 8	0.229 2
${{{\epsilon}}_{{{\max}}\;{\rm{WLS}}}}$	0.048 6	0.065 4	0.238 6	0.042 0	0.160 9	0.106 0	0.170 2
εOLS	0.058 9	0.076 1	0.164 7	0.068 9	0.149 9	0.126 5	0.071 6
εWLS	0.043 0	0.060 4	0.164 3	0.064 9	0.138 8	0.119 4	0.054 4
RMSEOLS	6.904 6	7.364 8	16.573 1	6.934 6	8.633 2	3.893 7	0.466 3
RMSEWLS	5.455 2	6.254 1	16.095 9	5.896 0	8.061 9	3.667 8	0.359 7

下载: 导出CSV 

| 显示表格

从表2数据中可以看出，除MODIS CH3通道WLS和OLS的定标评价参数相近外，其余通道WLS的定标效果均要好于OLS的定标效果。不同的MODIS通道在进行光谱匹配后获取的定标数据点对不同，各项K因子的不确定度计算结果也存在差异，因此WLS方法获取的定标结果也不尽相同。对于MODIS的1、2、4、5、6、7通道，WLS方法回归得到的定标系数与OLS相比，其定标平均相对误差的降幅在0.5%~1.5%范围内，在观测低亮度目标时，其定标相对误差降幅在3%~5%范围内。针对MODIS CH3通道WLS方法定标效果不理想的分析原因如下：MODIS CH3通道带宽较窄（19 nm左右），仅对应2个高光谱通道，相对于其它通道可对应4~9个高光谱通道而言，CH3通道采用光谱卷积得出辐射基准值的准确性较差，各项K因子的不确定度计算准确度下降权重系数计算不理想，导致最终WLS回归结果准确性下降。同时，可以发现采用WLS和OLS得到的回归系数与官方系数相比，斜率差异不明显，主要体现在截距的差异上。亮度较低的目标一般存在较大的不确定度，而WLS方法的引进可以有效修正低亮度目标对截距的干扰，使定标结果更接近真实值。还需要指出的是，HYPERION载荷自身的定标精度低于MODIS载荷的定标精度，因此，WLS计算所得的MODIS定标系数与官方定标系数还存在一定差异，待我国自研的高精度太阳反射波段高光谱基准载荷在轨后，可使用其作为参考仪器向MODIS传递溯源基准，将有望消减该差异，从而更准确地检验WLS方法的优势。

4.   结　论

本文提出的基于不确定度分析的加权最小二乘法的交叉定标方法，在得出近一致观测的数据点对后，利用蒙特卡罗法对各项匹配环节对辐射基准值的影响进行量化，针对每个数据点对都可算出其辐射基准值的不确定度，利用不确定度平方的倒数作为权重系数，采用WLS方法回归了定标系数。以HYPERION作为参考基准载荷，对MODISCH1~7通道的定标结果表明，1、2、4、5、6、7通道采用WLS法回归的定标系数的评价参数要优于OLS法，其中定标系数的定标平均相对误差的降幅在0.5%~1.5%范围内，在观测低亮度目标时的定标相对误差降幅在3%~5%范围内，说明WLS方法有效地提高了定标系数的精度，验证了基于不确定度的WLS方法的有效性。

在后续的工作中将进一步评估通道带宽较窄，对应的高光谱通道有限时，基于高光谱地物谱型库的光谱匹配方法对应的不确定度，同时在计算不确定度时进一步分析时间匹配和观测几何匹配引入的不确定度，从而进一步提高权重系数和定标数据点对的计算精度，提升WLS的回归效果。

致谢：本文由国家重点研发计划课题（2018YFB0504601）以及国家自然基金委员会面上项目（41271373）资助。感谢ladsweb.modaps.eosdis.nasa.gov提供MODIS数据及earthexplorer.usgs.gov提供HYPERION数据。