江苏省位于116°18′E~121°57′E、30°45′N~35°20′N，是长江三角洲经济带的重要组成部分，人均GDP、综合竞争力、地区发展与民生指数均居中国各省第一。

本文选择江苏省2015年1月-12月共12期的VCM月度NPP/VIIRS影像作为实验数据，数据来自美国国家海洋和大气管理局（National Oceanic and Atmospheric Administration，NOAA）。原始VCM数据中杂散光、闪电光已经被剔除，但是火光、船只以及其他临时性灯光未被剔除；此外，月度数据还存在负值像元。在数据预处理过程中，首先将负像元值赋值为0，再进行临时灯光去除。考虑到如果一个国家或地区没有自然灾害等重大事件发生，那么该区域的灯光亮度在连续的月份间不应出现剧烈变化。鉴于此，本文借助箱线图中稳健统计量来识别临时灯光像元。提取每个像元1月-12月有效像元值（非0值）构成时间序列$x$，其中$t$月份稳定灯光像元值$x_t$（1≤$t$≤12）满足：

式中，$Q_{\mathrm{1}}$为时间序列$x$的下四分位数；$Q_{\mathrm{3}}$为时间序列$x$的上四分位数。

当像元值$x_t$位于箱线图最大上限或最小下限以外时，认为该像元为数据缺失或异常像元，将其赋值为0，作为插值对象。经数据预处理后，江苏省2015年VCM月度影像待插值像元占总像元的比例见表 1，其中6月份待插值像元占比最高，达到40.84%。

本文插值方法的基本思路是：以时、空关系互相作为约束条件，将时序变化一致性较好的像元数据作为空间维度插值的参考，将空间关系一致性较好的月度数据作为时序维度插值的参考；利用参考值及其与待插值的差值进行像元值的初步估计，进一步根据时、空一致性关系，借助卷积运算求得待插值像元的时空插值结果。为了简便起见，本文将以时序变化一致性为依据，利用空间邻近像元作为参考值进行插值的方法简称为空间维度插值；以空间关系一致性为依据，利用时序邻近像元作为参考值进行插值的方法简称时间维度插值；将二者插值进行一步加权混合的方法称为综合时空插值，简称时空插值。

本文方法的实现过程借助卷积运算完成。卷积运算是一种在空间域中对图像进行邻域检测的算法，可用于图像平滑、信号处理和转换边缘检测等，其本质是通过图像邻域像素的加权求和得到卷积值，使用的加权矩阵称为卷积模板或卷积核。为了得到不同维度的插值结果，首先从时间和空间维度分别评估邻域数据集与包含待插值数据的数据集的一致性，并将该一致性指标的标准化值作为插值结果的可靠性权重系数；然后构建不同的卷积核计算时间和空间维度插值，进而将不同维度插值结果加权融合，得到缺失数据的最终插值结果。时空插值方法技术路线如图 1所示。

图  1  时空插值技术路线

Figure 1.  Flowchart of Spatio-temporal Interpolation

数据插值即从有限已知数据推测出未知数据值，选择合适的插值参考是插值过程至关重要的步骤。待插值数据集与参考数据集的一致性越高，插值结果精度越高。为此，构建一致性指标评价待插值数据集与参考数据集的时空一致性。设$\mathit{l}=[l_{\mathrm{1}},l_{\mathrm{2}}\dots l_n]$和$\mathit{L}=[L_{\mathrm{1}},L_{\mathrm{2}}\dots L_n]$为长度相等的两个向量，两个向量具有可比的数据尺度和单位，$\mathit{\Delta }=[{\Delta }_{\mathrm{1}},{\Delta }_{\mathrm{2}}\dots {\Delta }_n]$为向量$\mathit{l}$和L的元素差值向量，即${\Delta }_m=l_m-L_m$，$l_m$、$L_m$和${\Delta }_m$分别代表向量$\mathit{l}\mathbf{、}\mathit{L}\mathrm{和}\mathit{\Delta }$中第$m$个元素值。$\mathit{\Delta }$中元素的标准差代表向量元素的离散程度，标准差越小，元素${\Delta }_m$越趋近同一值，则向量$\mathit{l}$和L元素的值呈现同步变化规律。基于此，本文构建一致性指标$\mathrm{S}\mathrm{i}\mathrm{m}\left(\mathit{l},\mathit{L}\right)$，评价向量$\mathit{l}$和L的一致性，计算式为：

式中，$\overline{\Delta }$为向量$\mathit{\Delta }$中元素的平均值；$n$为向量$\mathit{\Delta }$的长度；$\epsilon$为一个正的无穷小量（避免分母为零的情况）；$\mathrm{S}\mathrm{i}\mathrm{m}\left(\mathit{l},\mathit{L}\right)$为$\mathit{l}$与$\mathit{L}$向量的一致性程度，其值越大，$\mathit{l}$与$\mathit{L}$向量的同步变化一致性越高，反之越低。

空间维度插值是根据数据缺失像元所在时间序列与邻域像元时间序列的同步性规律选择空间插值的参考像元，利用估算的待插值像元与参考像元的差值，完成待插值像元值的估计，结果如图 2（a）所示。

图  2  不同插值方法示意图

Figure 2.  Schematic Diagram of Different Interpolation Methods

设${\mathit{S}}_{i,j}$为$i$行$j$列的时间序列向量，$S_{i,j}^t$代表$t$月份$i$行$j$列的像元值，${\mathit{S}}_{i,j+\mathrm{2}}$为${\mathit{S}}_{i,j}$邻域内某一时间序列向量，$S_{i,j+\mathrm{2}}^t$代表${\mathit{S}}_{i,j+\mathrm{2}}$中与$S_{i,j}^t$处于同一月份的邻域像元值，当向量${\mathit{S}}_{i,j}$或${\mathit{S}}_{i,j+\mathrm{2}}$中存在其他缺失数据时，去除另一向量中的对应数据。假设$S_{i,j}^t$为缺失数据，计算向量${\mathit{S}}_{i,j}$和${\mathit{S}}_{i,j+\mathrm{2}}$向量的一致性指标$\mathrm{S}\mathrm{i}\mathrm{m}({\mathit{S}}_{i,j},{\mathit{S}}_{i,j+\mathrm{2}}$），同时求得向量${\mathit{S}}_{i,j}$和${\mathit{S}}_{i,j+\mathrm{2}}$中对应元素差值的均值${\overline{\Delta }}_{i,j+\mathrm{2}}$；以${\mathit{S}}_{i,j+\mathrm{2}}$作为插值参考时，数据缺失像元$S_{i,j}^t$的估计值计算式为：

式中，${\widehat{S}}_{i,j}^t$为数据缺失像元的估计值；$S_{i,j+\mathrm{2}}^t$为向量${\mathit{S}}_{i,j+\mathrm{2}}$中与$S_{i,j}^t$处于同一个时间节点相邻像元值；${\overline{\Delta }}_{i,j+\mathrm{2}}$为向量${\mathit{S}}_{i,j}$和${\mathit{S}}_{i,j+\mathrm{2}}$中对应元素差值平均值。

将数据缺失像元某个邻域（以5×5为例）内时间序列作为插值参考，计算待插值向量${\mathit{S}}_{\mathit{i},\mathit{j}}$与插值参考${\mathit{S}}_{i,j+\mathrm{2}}$的一致性指标$\mathrm{S}\mathrm{i}\mathrm{m}\left({\mathit{S}}_{i,j},{\mathit{S}}_{i,j+\mathrm{2}}\right)$，并将其标准化至[0, 1]，计算式为：

式中，$\mathrm{S}\mathrm{i}\mathrm{m}\mathrm{\text{'}}\left({\mathit{S}}_{i,j},{\mathit{S}}_{i,j+\mathrm{2}}\right)$代表一致性指标$\mathrm{S}\mathrm{i}\mathrm{m}\left({\mathit{S}}_{i,j},{\mathit{S}}_{i,j+\mathrm{2}}\right)$的标准化结果，取值范围为[0, 1]；$\mathrm{M}\mathrm{a}\mathrm{x}$、$\mathrm{M}\mathrm{i}\mathrm{n}$为数据缺失像元5×5邻域内时空一致性的最大值、最小值。

借助卷积运算构建空间维度插值矩阵$\widehat{\mathit{S}}$和空间维度卷积核$\mathit{K}$，计算式分别为：

式（5）中，${\widehat{S}}_{i,j+\mathrm{1}}^t$是以$S_{i,j+\mathrm{1}}$作为插值参考时，数据缺失像元的估计值。进而得到数据缺失像元的空间维度插值$R_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}}$，计算式为：

时间维度插值是根据数据缺失像元及其空间邻域数据所在的数据集与其他月份相同位置像元的一致性作为依据，选取时序参考数据进行插值，结果如图 2（b）所示。

${\mathit{S}}^t$为t月份包含缺失数据的空间数据集，$S_{i,j}^t$代表数据缺失像元值，${\mathit{S}}^k$为k月份与${\mathit{S}}^t$行列号相同的数据集，k≠t，$S_{i,j}^k$代表${\mathit{S}}^k$中与$S_{i,j}^t$处于不同时间节点但行列号相同的像元值。当向量${\mathit{S}}^k$或${\mathit{S}}^t$中存在其他缺失数据时，去除另一向量中对应数据，计算两个数据集的一致性指标$\mathrm{S}\mathrm{i}\mathrm{m}({\mathit{S}}^t,{\mathit{S}}^k$），将向量${\mathit{S}}^t$和${\mathit{S}}^k$按照先行后列顺序排列组成一维向量，同时求得${\mathit{S}}^t$和${\mathit{S}}^k$中对应元素差值的平均值${\overline{\Delta }}_k$；${\widehat{S}}_{i,j}^t$代表以${\mathit{S}}^{\mathit{k}}$作为插值参考时数据缺失像元$S_{i,j}^t$的估计值，计算式为：

类似地，以其他月份空间数据作为插值参考，分别计算待插值向量${\mathit{S}}^t$与插值参考的一致性指标，并将其标准化至[0, 1]，借助卷积运算构建时间维度插值矩阵${\widehat{\mathit{S}}}^{\mathit{\text{'}}}$和时间维度卷积核$\mathit{K}\mathit{\text{'}}$，计算式分别为：

式中，${\widehat{S}}_{i,j}^k$代表以${\mathit{S}}^{\mathit{k}}$作为插值参考时数据缺失像元$S_{i,j}^t$的估计值；$\mathrm{S}\mathrm{i}{\mathrm{m}}^{\mathrm{\text{'}}}\left({\mathit{S}}^t,{\mathit{S}}^k\right)$为一致性指标$\mathrm{S}\mathrm{i}\mathrm{m}\left({\mathit{S}}^t,{\mathit{S}}^k\right)$的标准化结果。进而得到数据缺失像元的时间维度插值$R_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}$，计算式为：

将数据缺失像元的空间维度插值$R_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}}$和时间维度插值$R_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}$融合，得到综合时空插值$R_{\mathrm{s}\mathrm{t}}$，计算式为：

式中，$q_{\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e}}$和$q_{\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{e}}$分别为空间和时间维度插值时参与运算的初步估计值的一致性指标之和。

图 3为江苏省2015年6月份不同插值方法影像插值结果，鉴于可作为参考影像的数据较少，同时便于与已有研究进行对比分析，本文选择VCMSL作为参考影像，将已有研究中插值精度较高的三次Hermite插值^[13]、空间维度插值和时间维度插值与时空插值方法进行对比，从插值结果的灯光亮度总和、绝对误差、相对误差以及与参考影像逐像元差值的整体方差等多个方面评价时空插值算法插值精度。

图  3  江苏省2015年6月份插值结果

Figure 3.  Interpolation Results of Jiangsu Province in June 2015

选择插值结果灯光值总和与参考影像进行对比，结果如图 4所示。从图 4可以看出，Hermite插值在1、6、7、11月份插值结果偏高，其中数据缺失最严重的6月份插值结果偏高现象最为明显，三次Hermite插值在5、6、7月份结果分别为163.12×10⁴ nW·cm^-2·sr^-1、154.60×10⁴ nW·cm^-2·sr^-1和146.52×10⁴ nW·cm^-2·sr^-1，6月份插值灯光亮度总和接近5月、7月份灯光亮度总和的平均值。分析其原因在于三次Hermite插值仅考虑了像元值的时序关系，忽略了NPP/VIIRS数据包含的重要空间信息，当插值点与两个插值端点距离相同时，插值结果接近两个插值端点的均值，因此出现6月份插值结果灯光总和接近相邻月份灯光总和均值现象，插值结果灯光亮度总和较参考影像偏高明显。空间维度插值以时序变化一致性较好的像元的数据作为插值的参考，顾及到像元的空间关联性，仍无法满足数据大范围缺失的插值要求，12个月份的插值结果中有11个月份低于同时期参考影像，6月份灯光总和偏低最明显。时间维度插值以空间关系一致性较好的月度数据作为插值参考，更侧重像元的时间趋势性，插值精度有一定提高，但6月份灯光亮度总和仍有较大偏差。时空插值以时、空关系互相作为约束条件，插值过程综合考虑了NPP/VIIRS数据时间趋势性和空间关联性，在4种插值方法中插值结果精度最高，单一月份插值结果灯光亮度总和未出现较大偏差，较同时期参考影像灯光总和最接近。

图  4  插值结果灯光辐射亮度值总和

Figure 4.  Sum Radiance Value of Interpolation Results

图 5为不同插值方法的插值结果与参考数据的误差对比。由图 5（a）可以看出，Hermite插值结果的绝对误差连续12个月均为最大，其中6月绝对误差为19.94×10⁴ nW·cm^-2·sr^-1，为全年绝对误差最大月份，12月灯光总和出现较大负偏差，绝对误差为-12.52×10⁴ nW·cm^-2·sr^-1。空间维度插值的绝对误差有11个月份为负偏差，整体插值结果偏低。时间维度插值的绝对误差在1、3、7月份均小于时空插值，但其他月份均大于时空插值。时空插值6月份灯光亮度总和的绝对误差仅为2.31×10⁴ nW·cm^-2·sr^-1，12个月份中有9个月份的绝对误差值在4种插值方法中最小。从图 5（b）可以看出，三次Hermite插值、空间维度插值和时间维度插值的最大相对误差均出现在6月份，最大相对误差分别为14.81%、9.42%和6.34%，而时空插值6月份相对误差仅1.04%，全年最大月份相对误差为4.85%，插值结果精度最高。此外，三次Hermite插值、空间维度插值和时间维度插值的相对误差在不同月份出现较大波动，而时空插值的相对误差较平稳，未出现较大波动；与其他3种插值方法相比，时空插值方法稳定性更好，插值结果不会因季节变化出现较大偏差。

图  5  插值结果与参考数据的误差对比

Figure 5.  Error Comparison Between Interpolation Results and Reference Data

为了验证插值结果与参考影像的整体吻合程度，将不同月份插值结果与同时期参考影像逐像元作差，计算逐像元差值的方差，结果如图 6所示。从图 6中可以看到，4种方法的插值结果与参考影像逐像元差值的整体方差较为接近，整体方差最大值均出现在6、7月份，其中时空插值方法逐像元差值的整体方差有10个月均为最小，表明插值结果影像与参考影像越相近；与其他3种插值方法相比，时空插值方法的插值结果与参考影像整体吻合度最高。

图  6  插值结果与参考影像逐像元差值的方差

Figure 6.  Variance of Differences Between Interpolation Results with Reference Data

本文以数据缺失最严重的6月份为例，选择图 3中区域1和区域2作为两个局部区域，进行不同方法的局部插值精度对比，结果见表 2。

由表 2可知，相较于Hermite插值，时空插值在区域1的灯光总和相对误差由25.01%降低至13.87%，与参考影像逐像元差值方差由2.69 nW·cm^-2·sr^-1降低至2.55 nW·cm^-2·sr^-1；区域2的灯光总和相对误差由28.24%降低至18.52%，逐像元差值方差由5.48 nW·cm^-2·sr^-1降低至4.86 nW·cm^-2·sr^-1。总体来看，时空插值的局部插值精度最高，时间维度插值次之，空间维度插值效果最差。

通过求均值的方式将月度影像合成年度影像，与NOAA提供的VCM_ORM（VIIRS cloud mask-outlier removed）年度影像进行对比，结果见表 3。

由表 3可知，时空插值合成的年度影像的灯光总和、像元平均值与年度影像更接近，灯光总和相对误差低至2.20%。4种插值方法的像元最大值在270 nW·cm^-2·sr^-1左右，与年度影像相比偏差较大，查看原始影像数据发现，VCM_ORM年度影像中像元值大于282.87 nW·cm^-2·sr^-1仅有1个像元，且该像元5月份辐射亮度值为1 155.14 nW·cm^-2·sr^-1，数据预处理时，该像元值被识别为临时灯光像元，影像插值时，该像元被重新赋值，因此该像元插值结果小于VCM_ORM年度影像像元值。此外，Hermite插值时所有插值结果均在已有数据范围内不会出现插值结果越界现象，而时间维度和空间维度插值可能出现插值结果为负值，当插值结果为负值时，认为该插值结果无效，将插值结果及其权重赋值为0，因此时空插值方法插值结果未出现负像元值。通过逐像元差值的整体方差对比发现，时空插值结果与年度影像插值整体方差最小，与年度影像吻合度最高。

针对NPP/VIIRS月度影像数据缺失问题，本文以2015年江苏省NPP/VIIRS月度影像插值为例，将时、空关系互相作为约束条件，分别进行了空间维度插值、时间维度插值和时空插值实验，得到的主要结论如下：

1）空间维度插值虽然顾及到像元的空间关联性，仍无法满足数据大范围缺失的插值要求，插值结果整体偏低。

2）时间维度插值考虑到像元的时间趋势性，插值精度较空间维度插值有一定提高，但部分月份插值结果有较大偏差。

3）相较于三次Hermite插值，时空插值获得的月度数据灯光亮度总和最大相对误差由14.81%降至4.85%，年度影像灯光亮度总和相对误差由2.44%降至2.20%，逐像元差值整体方差由1.25降至1.20。

总的来看，本文中时空插值方法在插值过程中同时顾及到NPP/VIIRS数据时间平稳性和空间结构稳定性，插值精度明显提高，且对待插值月份前后时序数据没有严格要求，更具广泛性。

值得注意的是，邻域参考数据范围的大小对时空插值精度和运算效率都有较大影响。通过多次实验发现，选择数据缺失像元3×3邻域像元作为插值参考时，可作为插值参考的数据集较少，插值结果较差；选择5×5邻域像元为参考时，可以得到较好的插值结果；如果继续扩大参考像元邻域范围至7×7，将导致计算量显著增加，算法效率降低，因此，本文以数据缺失像元5×5邻域内数据集作为插值参考。虽然本文中综合时空插值精度更高，但是该方法也存在一定不足之处，即在插值起始位置必须找到有效的参考数据。如果某一月份数据全部缺失，则需要扩大插值范围，从可以找到有效参考位置进行影像插值，然后提取所需范围影像，由此可能会带来计算量的增加，对时间成本有较高的要求。