植被指数时间序列在农业、自然资源规划等方面具有参考价值，常用于作为当地农业资源规划、作物长势监测的依据。然而目前常用的植被指数时间序列数据大多存在空间分辨率低或时间分辨率低的问题，难以满足实际需求。因此植被指数时间序列的重建方法成为解决当前实际需求的有效手段，目前常用的植被指数时间序列重建方法可以分为基于滤波的方法和基于函数的方法2类。

基于滤波的方法是在原有植被指数时间序列的基础上通过滤波的方式去除噪声影响，以达到提高数据质量的目的。如基于Savitzky-Golay（SG）滤波拟合法，这个方法是在时域内基于局域多项式最小二乘法拟合的滤波，这种滤波法在去除噪声的同时能够保证信号的形状和宽度不变，在长时间序列植被指数重建方面得到广泛应用^[1]。边金虎等^[2]同样采用SG滤波，对若尔盖高原地区的MODIS植被指数时间序列进行重构，并与傅里叶变换法进行了对比，然而仅采用SG滤波的方法在应对连续噪声影响时会产生无法有效去除的问题。基于这个问题，胡顺石等^[3]在此基础上增加了自适应参数采用自适应加权的SG滤波重构了MODIS植被指数时间序列，解决了多数植被指数序列重构方法不能有效去除连续云噪声的问题。通过SG滤波以及基于此滤波方式进行改进的重建方法是滤波重建的主流方法，也存在另外的滤波方法对植被指数时间序列进行重建。如除SG滤波外傅里叶谐波分析法也是常用的滤波方法，这个方法是将时间域的信息先变化到频率域上，再对频率域上的信息进行低通滤波处理^[4]，如王鹏新等^[5]通过傅里叶变换对研究区域主要农作物进行识别提取；然而单一的傅里叶谐波分析法无法在所有场景下都达到很好的效果，因此产生了通过加上变量权重的改变因素的移动加权谐波分析^[6]。通过滤波方式对植被指数时间序列重建是在已有植被指数时间序列的基础上对结果进行去噪处理，对原始数据的空间分辨率有很高的要求，且在时间序列存在空缺值的情况下难以对空缺值进行预测和填补。

基于函数的方法即函数拟合法。非对称高斯拟合法^[7]是一种分析时间序列数据时采用的策略，它首先识别出数据中的高点或低点，并将这些点所在的区域作为局部拟合的目标区间；接着，对这些选定的局部区域应用高斯函数进行数据拟合，高斯函数能够很好地模拟数据的分布形状；最后，将这些局部拟合的结果整合起来，形成一个统一的全局拟合曲线。这个过程需要通过一些优化技术来调整局部拟合的参数，以确保它们在整个数据集上能够平滑地连接。这种方法能够有效地结合数据的局部特性和整体趋势，但同时也需要仔细选择局部区间，并在全局拟合时进行精细的参数调整。如邵亚婷等^[8]采用非对称高斯拟合法重构MOD13Q1产品，用以分析蒙古国植被物候空间分布及年际变化趋势。双逻辑拟合法^[9]通过结合2个逻辑函数来创建一个复合模型，以拟合具有2个显著特征或过渡区域的数据。这种方法通过调整逻辑函数的参数（如斜率、中点和幅度），使得模型能够捕捉数据的复杂动态。樊辉^[10]采用双逻辑拟合法与其他插值滤波方法重构了云南省归一化植被指数（normalized differential vegetation index，NDVI）时间序列达到分析NDVI时空分布的目的。目前使用函数拟合法进行植被指数时间序列的重建更多地还是从数据本身出发，采用数学方法将已有的数据点进行拟合，却忽略了作物在实际的生长过程中本身特有的生长规律，且重建效果的好坏也在很大程度上依赖于参数选取以及作物生长期的自主划分，难以做到对局部地区的特质化参数选择以及对特定作物生长期的精准划分。Sun等^[11]采用参考曲线法重建植被指数时间序列的工作，在函数拟合的基础上通过构建当地目标作物的作物参考曲线，在包含研究区域内作物生长期特质性信息的基础上对目标作物进行植被指数时间序列的重建。因此本文采用能够综合考虑研究区内作物特质性的作物参考曲线法进行植被指数时间序列重建工作。

综上所述，目前长时间序列的植被指数重建的方法仍存不足。如在数据量较少的情况下针对某个小区域内的作物进行高时空分辨率的植被指数时间序列重建时，目前的重建方法往往会遇到重建的数据量不足、重建结果的时空分辨率较低、重建结果无法准确反映作物生长期的变化等问题。为解决目前重建方法遇到的问题且同时考虑提高目前国产卫星的应用情况，本文应用车向红教授的反射率标准化方法^[12]，以VIIRS数据为参考数据源，协同高空间分辨率的多源卫星数据（GF-1，GF-6与Sentinel-2）进行反射率标准化后作为作物参考曲线重建方法的输入数据，提高可用数据的数量达到提高重建结果质量的目的；并且使用高时间分辨率的VIIRS卫星数据构建目标作物的NDVI和增强型植被指数（enhanced vegetation index，EVI）时间序列库，以达到对区域内目标作物植被指数变化规律的精准描述，解决传统函数拟合法重建结果无法准确反映作物生长期的问题。

本文研究区域涵盖衡水市冀州区、武邑县等区县，范围为114.99°~115.92°E，37.41°~37.75°N，研究区Sentinel-2卫星影像合成真彩图如图1所示。研究区域大部分属河北省衡水市，衡水市位于河北省东南部，属于冲积平原地形，地跨115.17°~116.57°E，37.05°~38.38°N。衡水市下设有11个县级行政区，包括2个市辖区（桃城、冀州）、1个县级市（深州）、8个县（武邑、饶阳、故城、阜城、枣强、武强、安平、景县）。衡水市地势平坦，海拔在12~30 m。全市具有充沛的光照及热量，年日照小时数可达2 500 h，年平均气温约13 ℃。衡水市的主要土壤类型为潮土，潮土亚类及脱潮土的亚类占衡水市全市土地总面积的82%^[13-14]。基于衡水市的地形及光照、水土等因素，衡水市具备冬小麦的适宜生长条件，因此选择衡水市的部分区域作为本文研究区域，对于研究本文目标作物——冬小麦具有一定的区域代表性^[15-16]。衡水市的冬小麦通常会在10月上旬至中旬进行播种，次年2月下旬至3月上旬冬小麦进入返青期，4月上旬进入拔节期，5月上旬开花，6月上旬成熟^[17]。

本文使用的数据源分为4类：用于进行目标作物像元识别的作物分布数据、用于参与植被指数时间序列重建及反射率标准化工作的多卫星反射率数据、用于构建参考曲线数据的VIIRS数据及用于作为卫星反射率标准化基准的VIIRS地表反射率数据集。本文使用2021年的全球冬小麦作物分布数据，空间分辨率为10 m，选择了包含研究区域的2021年冬小麦物候期1—6月的反射率数据，其中涵盖了高分系列卫星（GF-1和GF-6）以及哨兵系列卫星（Sentinel-2）的反射率数据。本文使用数据具体信息如表1所示。

本研究主要技术路线如图2所示。首先准备了VIIRS数据、目标作物分布数据以及多种卫星反射率数据，在对数据进行必要的预处理后，对经过预处理的GF-1，GF-6，Sentinel-2数据进行反射率标准化工作。在对GF卫星以及Sentinel-2卫星进行预处理时，本文采用ENVI中的PRC Orthorectification Workflow模块进行正射纠正，为达到更好的纠正效果，在进行处理时将GF及Sentinel-2卫星的分辨率重采样至15 m；然后进行目标作物的NDVI及EVI时间序列重建，首先制作经过SG滤波的NDVI及EVI时间序列参考曲线数据库，得到研究区域的NDVI及EVI参考曲线数据库；接着对多卫星反射率数据进行NDVI及EVI的计算，同时记录对应的时间信息用于作物参考曲线法的植被指数时间序列重建；再使用研究区域的目标作物分布数据对数据进行目标作物像元的筛选及判断；最后输入至作物参考曲线法重建代码中完成研究区域内冬小麦的NDVI及EVI时间序列的重建。

由于传感器间光谱响应的差异以及校准和大气校正过程中的残留误差，地表反射率的传感器间不一致性问题依旧存在。这种不一致性很容易与地表变化相混淆，尤其是在监测植被对气候变化的细微响应方面。鉴于此，本研究采用一种贝叶斯方法来标准化GF-1，GF-6和Sentinel-2卫星的反射率数据。为了确保后续建模过程的顺利进行，本研究首先应用点扩散函数（point spread function，PSF）对GF系列卫星和Sentinel-2卫星数据进行空间退化，并通过与VNP09GA相应波段的相关性分析来确定最优的σ值作为PSF参数，该参数能够使二者的决定系数达到最大。确定PSF参数后，本研究利用最大后验概率（maximum a posteriori，MAP）估计，并采用共轭梯度下降法进行数值求解，实现了GF-1，GF-6卫星和Sentinel-2卫星的反射率图像与VNP09GA图像的融合，构建了一套统计关系模型，从而推导出标准化后的GF-1，GF-6和Sentinel-2卫星反射率图像数据。

1）使用点扩散函数进行空间退化。为了在之后在贝叶斯建模过程中能够将GF卫星及Sentinel-2卫星的高分辨数据与低分辨率的VIIRS数据匹配，本文需要对GF卫星及Sentinel-2卫星数据进行空间退化，退化为与VIIRS分辨率一致。本文采用PSF来进行空间退化，计算公式为：

式中：L_{1 000 m}（λ）为目标变量，代表波长为λ的高分辨卫星空间退化后的结果；psf（λ）为卷积核矩阵；*为卷积操作；L_{15 m}（λ）为输入的原始高分辨卫星影像；subsample为子采样操作。

然而psf（λ）的卷积核矩阵在进行每次卷积操作时需要获取对应的权重值，这个权重值受到传感器光学、图像运动等的影响，因此本文选择使用传统的高斯估计，计算公式为：

式中：psf_λ（i，j）为目标变量，代表滤波器的权重；i，j可以取1，2，3，且由于该滤波器是45×45的滤波器，因此i，j均减去23使得位于该滤波器中心；σ_λ则是需要通过对GF卫星及Sentinel-2卫星数据进行退化后与VIIRS影像（即本文的VNP09GA数据）之间的相关系数得出。σ_λ的计算公式为：

式中：R²（L_{1 000 m}（λ，σ_λ），V_{1 000 m}（λ））为L_{1 000 m}（λ，σ_λ）与V_{1 000 m}（λ）2个变量R²值；L_{1 000 m}（λ，σ_λ）为由等式（1）中由高分辨GF卫星、Sentinel-2卫星数据导出的空间退化后的低分辨率数据；V_{1 000 m}（λ）为原始VIIRS数据。

基于以上公式，本文对GF-1，GF-6，Sentinel-2卫星数据的退化后数据与VNPO9GA数据的蓝光波段、绿光波段、红光波段、近红外波段之间进行R²的计算，并从中选取最大的R²值对应的σ_λ。

2）构建统计关系模型进行反射率标准化。卫星反射率标准化工作在本文中是通过结合GF数据、Sentinel-2反射率数据的空间信息以及VIIRS反射率图像的光谱信息获取的。这些信息是通过构建GF卫星、Sentinel-2卫星反射率数据与目标图像以及VIIRS的反射率数据与目标图像的统计关系，并利用MAP求解获取的。假设一个正态分布，计算公式为：

式中：v_{15 m}（λ）为15 m分辨率标准化后的GF卫星或Sentinel-2卫星图像的一维形式，也是本文需要获取的目标图像；P（v_{1 000 m}（λ）|v_{15 m}（λ））为在给定v_{15 m}（λ）的条件下v_{1 000 m}（λ）的概率；k为向量v_{1 000 m}（λ）的长度；${U}_{\lambda }^{VIIRS}$为VIIRS的反射率不确定矩阵；S_λ为子采样操作；H_λ为PSF卷积运算。

接着继续假设v_{15 m}（λ）的卫星图像具有给定l_{15 m}（λ）的多元正态分布，计算公式为：

式中：P（v_{15 m}（λ）|l_{15 m}（λ））为在给定l_{15 m}（λ）的条件下v_{15 m}（λ）的概率；n为v_{15 m}（λ）的长度；${U}_{\lambda }^{GF}$为GF卫星的反射率不确定矩阵（同样适用于本文使用的Sentinel-2卫星数据）。

最后便是使贝叶斯最大后验概率最大，在本文中即是使式（4）及（5）的指数部分的和达到最小，将输出满足该条件的结果以获取反射率标准化后的卫星数据。将以上推导公式集成于卫星反射率标准化代码中，以实现多卫星反射率标准化工作。

本研究选用了GF-1，GF-6和Sentinel-2经过标准化处理后的15 m分辨率反射率数据。经过反射率标准化处理后，这些数据被用作作物参考曲线法的输入数据，并与目标作物的分布数据相结合，以实现精确的植被指数时间序列重建。本研究首先从参考曲线库中挑选标准参考曲线，采用最小二乘法挑选与输入数据点拟合情况最好的曲线作为最优曲线；然后，进行曲线平移，通过平移最优曲线的位置获取能够最好表达输入数据点的最优曲线位置；最后，局部调整最优作物参考曲线，使得调整后的曲线更加贴合输入数据点的实际情况。

1）制作参考曲线数据库。为构建本文使用的参考曲线数据库，本研究利用Google Earth Engine（GEE）平台获取了研究区域内冬小麦生长关键期2021年1—6月的VIIRS数据。通过数据处理及筛选无云影响的像元的操作，获取了无云条件下的卫星反射率数据，进而初步推导出了NDVI时间序列。然而，由于其他恶劣天气条件及传感器等问题的影响，初步计算得到的NDVI时间序列可能包含噪声。为了消除这些不可避免的干扰，本研究采用了一种简化的SG滤波算法对数据进行处理，从而得到了高质量且记录了各个像元点信息的NDVI和EVI时间序列，并将其输出为TIFF格式的数据集。本文使用的参考曲线库如图3所示。

2）挑选最优作物参考曲线并进行曲线局部调整。最优参考曲线能够反映不同像元中作物最有可能的生长情况，因此对不同像元挑选作物参考曲线是重建结果的保证。本文研究假设在作物的生长周期中，GF系列数据及Sentinel数据在反射率标准化后与VIIRS数据在像元尺度上具有线性关系，因此根据式（6），将参考曲线库中的曲线进行平移并逐一计算每条曲线在不同平移后位置以最小二乘法的拟合情况，选取出决定系数最大且拟合情况最好的曲线，再根据实际计算出的植被指数值进行局部调整，计算公式为：

式中：x为该景影像每年积日（day of year，DOY）；S（x）为多卫星标准化后的植被指数函数；x₀为由于像元物候期差异的时间偏移量（本文选取的偏移量范围为-30—30 d）；a，b为最小二乘拟合法得到的参数；M（x）为NDVI/EVI参考曲线函数。

本文反射率标准化结果与精度评价采用统计指标变化的形式进行定量评价。采用决定系数（R²）、均方根误差（root mean squared error，RMSE）、相对均方根误差（relative RMSE，rRMSE）3种统计参数进行评估，计算公式为：

式中：y_i为标准化图像像元值；$\stackrel{︿}{{y}_{i}}$为VIIRS图像像元值；n为图像像元数量；$\stackrel{-}{{y}_{i}}$为VIIRS像元平均值。

从评价指标来看，反射率标准化工作显著提升了R²的值、降低了RMSE和rRMSE的值。针对GF-1，GF-6卫星和Sentinel-2卫星，在蓝光，绿光，红光和近红外波段R²分别有所上升；RMSE和rRMSE分别有所下降。GF-1卫星与VIIRS地表反射率数据R²大部分提高0.05，少部分提高超过了0.1。RMSE降低，大部分RMSE降低了0.01，rRMSE降低幅度在2%左右。GF-6卫星大部分R²提高了约0.12，RMSE大部分减小了0.03，rRMSE普遍减小，幅度在3%~4%之间。Sentinel-2卫星R²整体提升约0.05，RMSE及rRMSE的降低大部分在0.001及2%左右。具体如表2—4所示。

本文重建了NDVI及EVI时间序列，重建结果如图4和图5所示。为准确描述重建结果，选取了偏差（Bias）、相对误差（relative error，RE）、RMSE和R²这4个统计指标对重建结果进行描述，RE和Bias的计算公式为：

式中：$\overline{V{I}_{i}}$为验证数据组中卫星影像计算得到的NDVI或EVI所有像元的平均值；$\stackrel{̑}{V{I}_{i}}$为重建组中卫星重建后影像像元的NDVI或EVI值；VI_i为验证数据组中卫星影像像元的NDVI或EVI值；n为样本的总数。具体重建结果精度评价表见表5。

根据表5显示，所有验证数据集的原始NDVI图像平均值与重建后的NDVI图像平均值非常接近，两者之间的最大差异仅为0.07。RMSE作为衡量重建结果与原始影像差异的指标，其数值越小，表明重建结果越接近真实情况。在本研究的验证数据中，RMSE的最大值为0.10，出现在DOY 117和127天，而最小值为0.03，出现在DOY 42和47天，其余时间点的RMSE均不超过0.1。RE是另一个描述重建准确性的指标，通常以百分比形式表示，反映了重建结果相对于真实值的偏差。在本研究的NDVI时间序列验证中，RE在绝大部分情况下小于15%，仅1景验证影像超18%，出现在DOY 28天，而最小值为3.90%，出现在DOY 125天，所有验证天数的RE大致在10%。R²值衡量了重建结果与原始影像之间的线性相关性，反映了重建结果的拟合程度。在本研究中，大多数验证数据的R²值超过了0.5，并且在DOY 47和125天时，R²值达到了0.8以上，显示出良好的拟合效果。

表6展示了EVI时间序列重建精度结果，与NDVI时间序列的重建结果相比，EVI时间序列在DOY 28，47和47天的RMSE指标上显示出更小或相当的水平。然而，从整体误差来看，EVI时间序列的重建误差仍然高于NDVI时间序列。在RE指标方面，EVI时间序列的重建结果在大多数验证期间保持在20%以下，多数情况下维持在13%左右，但部分验证日期的RE偏差超过了20%，特别是在EVI指数快速上升或下降的DOY 82，127这几天。这些日期的EVI变化幅度较大，且可用于重建的数据量较少，导致了明显的偏差。从R²来看，EVI时间序列的重建结果有5景影像的R²不低于0.44，而在DOY 47和125天时，R²值超过了0.7，显示出较好的拟合度。尽管如此，仍有部分日期的R²值较低，如DOY 82天，这些日期在RE指标评估中也表现出较大的偏差，这一现象在NDVI时间序列的重建结果中也有所体现。

综合来看，尽管EVI时间序列的重建结果在大多数验证日期中展现出较小的RMSE和RE以及较大的R²，但在EVI指数快速变化的时期，重建结果的偏差相对较大。与NDVI时间序列的重建结果相比，尽管EVI时间序列在整体上显示出较低的偏差和较高的与原始影像相关性，但由于EVI时间序列在快速变化时期的波动更大，且缺乏明显的波峰期，因此在重建效果上略逊于NDVI时间序列。

卫星反射率标准化前后冬小麦的NDVI时间序列图见图6。

图6展示了标准化前后冬小麦NDVI时间序列的差异，从中可以直观地观察到卫星反射率标准化工作后生成NDVI时间序列图的变化。图中的红色垂直虚线用以划分冬小麦的物候期开始时间。从图中可以观察到，在返青期和抽穗期附近，原始的NDVI曲线出现了一次显著的下降，同样在抽穗期DOY 125—127天，NDVI下降幅度超过了0.1的异常下降。此外，在DOY第127—147天，即抽穗期之后，NDVI出现了一次异常的上升。经过反射率的标准化处理，返青期附近的异常下降现象得到了有效纠正，NDVI曲线呈现出了平稳上升的趋势。同样，在抽穗期内超过0.1 NDVI值的异常下降也得到了显著改善，使得处于下降阶段的DOY第127天的异常NDVI值，经过修正后，回归到了一个合理的范围内。

为满足农业上对于作物生长监测及土地利用类型变化等需求，本研究选取了河北省部分区域作为研究区域，以冬小麦作为本文的目标作物，采用GF-1，GF-6，Sentinel-2数据对研究区域进行NDVI时间序列及EVI时间序列的重建。本文首先使用VIIRS的产品VNP09GA作为基础，对GF-1，GF-6，Sentinel-2数据预处理后进行数据标准化工作，并将标准化后的数据与VNP09GA数据及时进行定量的数据验证。在获取各个卫星标准化的反射率数据后，本文使用VIIRS数据结合作物分布数据制作出研究区的EVI及NDVI的参考曲线库。结合EVI参考曲线、NDVI参考曲线、冬小麦分布数据、反射率标准化后数据，采用作物参考曲线法对研究区域内的冬小麦像元进行EVI与NDVI时间序列的重建，重建出研究区域内空间分辨率为15 m、时间分辨率为1 d的高时空分辨率植被指数时间序列。基于重建结果，选取了一些日期进行精度验证，主要结论如下：

1）通过对GF-1，GF-6，Sentinel-2卫星反射率数据的标准化工作，获取了3种卫星数据在以VNP09GA数据为基础上标准化后的结果。经过精度检验，标准化后的卫星在红光波段、蓝光波段、绿光波段以及近红外波段的反射率数据整体与VNP09GA数据之间的R²提升，rRMSE及RMSE降低。该反射率标准化的方法使得不同卫星之间的反射率数据能够以某一卫星反射率为标准，提高不同卫星反射率数据与标准化卫星反射率数据之间的相关性并降低整体误差，实现不同的卫星数据能够更加可靠地联合使用的目的。

2）基于VIIRS数据，提取了河北省内2021年的NDVI参考曲线与EVI参考曲线。该参考曲线能够直观地反映出2021年河北省冬小麦NDVI，EVI的变化趋势与变化速率。以此参考曲线为基础，可实现物候信息提取、土地利用变化分析等。

3）基于NDVI和EVI参考曲线、作物分布数据、反射率标准化后卫星数据，采用作物参考曲线法重建出了研究区内15 m空间分辨率、1 d时间分辨率的NDVI时间序列及EVI时间序列结果。该重建结果携带空间信息及时间信息且具有较高的时空分辨率，能够作为基础信息为当地提供生态环境评估、农业管理、气候变化研究、灾害监测评估等服务。该重建方法没有限制输入数据的数量和质量，能够通过建立当地植被指数参考曲线库的形式，在输入数据较少的情况下，根据已有输入数据及作物参考曲线携带的信息对时间维度上缺失的指数值进行合理填补，并依据已输入的数据对植被指数时间序列的结果进行修正。能够更好地解决对某个小区域内的作物进行植被指数时间序列重建的时候，由于可用的数据量不足而无法进行重建或者无法准确代表作物生长期变化的问题。且本文将该重建方法与卫星反射率标准化方法进行结合，使得在进行植被指数时间序列重建时获得了更多的可输入数据，提高了重建结果的质量。

4）本文将反射率标准化与作物参考曲线法结合，重建了河北省部分区域的冬小麦NDVI以及EVI的高时空分辨率时间序列。然而该结果仅在单一区域进行了单个作物的植被指数时间序列重建，无法完全证明在其他区域及其他作物上同样较好的表现，因此在其他区域及作物上验证本文方法也是今后研究的方向之一。