Application of generative remote sensing image super-resolution reconstruction technology in crop phenological extraction

研究区位于辽宁南部辽河出海口处，春季和夏季云层降水量较高，秋冬降水量较低。辽河平原是重要的东北黑土地重点保护区域，辽中南地区的“鱼米之乡”，也是中国最重要粮食主产区之一。在研究区选取了大面积的水田和旱田作为物候提取试验区域，试验区域示意图见图1。

Figure 1.  Schematic diagram of the study area

收集2023年在40.84°~40.92°N，121.98°~122.07°E拍摄的全部无云、无噪声的Sentinel-2影像，作为原始时序影像。利用欧洲航天局提供的 Sen2cor 插件对下载的影像进行辐射校正和大气校正，生成 L2A级数据，用于全年的时序监测与物候提取应用。在水田区域共获取有效时序影像57幅，平均时间分辨率为6.40 d，旱田区域共获取有效时序影像55幅，平均时间分辨率为6.64 d。每幅图像均包含了物候提取最必要的4个颜色通道，即红、绿、蓝和近红外波段。

①超分辨率网络训练数据集。在训练数据集的制作方面，首先将研究区不同季节的20幅Sentinel-2影像经过双三次下采样为40 m，制作4倍超分辨率重建的基础数据集，用于大规模预训练。同时在训练数据中加入了79组Sentinel-2和Landsat-8相同日期拍摄的影像对，制作异源影像数据并基于这些数据精细化训练生成式神经网络。②超分辨率影像数据生成方法。在数据集上训练超分辨率重建网络，并对网络的图像重建能力进行定量评价。为了适应时序影像的超分辨率需求，设计了一个简洁、轻量化的PGT-GAN模型(图2)，该模型兼具目前主流网络的生成性能和简洁性。它的生成器由深度特征提取和分辨率升级两部分组成，首先将分辨率影像数据输入到由残差模块(residual block)和残差STB(residual swin transformer block)构成的深度特征提取模块中，而后经多组由卷积层、像素归一化层和LRelu层逐步进行分辨率升级。判别器部分则沿用PGGAN的模式，经卷积层和LRelu层逐级下采样，最后归一化并计算损失函数。所提出的PGT-GAN方法使用Tensorflow 2.8.0框架实现。在训练过程中应用Adam优化器，初始学习率为0.002，批处理大小为4。通过多次测试和比较，将绝对误差损失(L1 loss)设置为200，深度(depth)设置为64。研究中使用上述设置对该模型进行了500个训练轮次(epoch)的训练。所有的实验均是在配备了NVIDIA RTX P5000 GPU的Ubuntu 20.04系统上进行。

Figure 2.  PGT-GAN network structure

选取最广泛应用的2个图像评价指标：峰值信噪比(peak signal to noise ratio, PSNR)和结构相似性指数(structural similarity，SSIM)^[19]，用于评估单幅遥感影像超分辨率重建效果。PSNR数学定义如式(1)所示：

式(1)中：R_PSNR指PSNR；C_MAX为图像中颜色的最大值；E_MS为均方误差。

SSIM是基于人类视觉提取图像中结构化信息的假设，是基于亮度、对比度和结构来衡量幅图像相似度的指标。计算如式所示：

式(2)中：R_SSIM指SSIM；x, y分别代表两幅对比影像中相对应的像素；l(x, y)用于衡量亮度差异；c(x, y)用于衡量对比度差异即像素的标准差；s(x, y)用于衡量结构差异，α、β、γ分别表示不同特征在SSIM衡量中的占比。

将提出的超分辨率时间序列重构方法与插值填补法和原始数据提取的物候信息在多个平滑条件下进行鲁棒性和实用性对比。这些平滑模型已被证明具有有效性，包括滑动平滑法(Moving)、局部加权平滑(locally weighted scatterplot smootling, LOWESS)^[20]、稳健局部加权回归(RLOWESS)^[21]、Savitzky-Golay滤波法^[22]。在目前物候提取研究常使用的插值方法中，高斯过程回归(Gaussian process regression，GPR)方法^[23]是填补时间序列影像空白的首选方法之一。此外，它也是诸多插值方法中唯一能提供相关不确定性的方法。本研究采用GPR和常用的K近邻回归(K-nearest neighbors regression，KNNR) 2种虚拟填补时序影像剩余缺失值的方法作为对比研究。

本研究的技术路线包括4个主要阶段：数据预处理、模型训练、时序重建和物候提取。数据预处理阶段被设计用于生成训练所需的输入数据和目标，并将测试图像预处理为网络的输入格式，即异源影像的影像对。在第二步模型训练中，采用2.1.2中生成的训练数据训练生成式深度学习模型，将Landsat-8影像输入到生成式神经网络，重建为与Sentinel-2光谱类似的10 m分辨率影像，并与原始10 m分辨率Sentinel-2影像构成更加密集的时间序列。在物候提取阶段，本研究从重建后的影像中计算植被指数，并提取多个物候指标进行比对分析。

重建性能分析以2023年11月4日的Sentinel-2 (10 m)和Landsat-8 (40 m)异源影像为例，以验证本研究提出方法的分辨率重建结果。如表1所示：通过和ESRGAN、pix2pix和SwinIR定量对比可以看出，本研究提出的PGT-GAN超分辨率重建网络和目前主流超分辨率网络相比略有优势，相较于SwinIR， SSIM和PSNR分别提升0.004和0.13。

通过图像的定性对比(图3)可以看出：PGT-GAN网络具有较强的影像还原能力，重建的影像纹理细节饱满，与原图色差很小，可以用于后续时序重建及物候提取。

Figure 3.  Comparison of results of super-resolution reconstruction methods　　　　

使用3.1中重建的遥感数据对原始影像数据进行填补，水田和旱田分别由每年57、55幅影像分别提升为每年64、62幅，平均观测回访周期分别由6.40、6.63 d提升至5.70、5.88 d，并进行基于重建数据的物候提取。为验证该方法提取物候的鲁棒性和实用性，分别采用4种不同的平滑方法，使用NDVI提取水田和旱田的多个物候指标，分别为：每年的季初日期(start of the season, SOS)、季末日期(end of the season, EOS)、振幅(amplitude)，生长周期长度(length)、最大值(max value)和峰值日期(max day)。图4~5为水田Moving平滑法和Savitzky-Golay滤波法的可视化提取结果，图6~7为旱田LOWESS和Savitzky-Golay滤波法的可视化提取结果。

Figure 4.  Various phenological indexes of paddy field extracted by Moving filter

Figure 5.  Various phenological indexes of paddy field extracted by Savitzky Golay filter

Figure 6.  Various phenological indexes in dry field extracted by LOWESS filter

Figure 7.  Various phenological indexes in dry field extracted by Savitzky-Golay filter

由图4~7可见：在相同的平滑方法下，重建后的时序数据和原始数据在各个物候指标的制图上有着明显的差异。水田和旱田不论是通过平滑方法还是Savitzky Golay滤波法提取的SOS和EOS都出现了多处差异，Moving方法提取的峰值日期制图差异最为明显，原始数据的峰值日期更为靠后。重建后的影像可以在很大程度上减少提取SOS和EOS时产生的异常值点，采用Savitzky Golay滤波法提取的EOS更为靠前，在最大值和振幅方面都有明显差异，原始数据相较于重建后数据振幅偏大，这表示原始数据在提取物候时稳定性欠佳。综上，无论是Sentinel-2原始影像还是重建后的影像，相较于MODIS数据提取的物候结果都更为精细。MODIS虽具有很高的时间分辨率，但如图8所示，仅适用于广泛的大面积检测。

Figure 8.  Schematic diagram of the phenological indicators extracted from MODIS data

表2~3分别表示各方法提取的水田、旱田的SOS、EOS、生长周期、峰值日期、最大值和振幅平均值。由于计算的是每个像素的平均值，因此即便是很小的数值差距也表明整张图巨大变化。通过定量对比可以发现：相较于插值方法对数据缺失部分进行的填补，本研究提出的填补方法在各个平滑方法的对比中，与原始的结果都存在一定差异。SOS、EOS的差异分别在±4、±2 d，物候期的长度差异在±3 d。

由图9时间序列曲线可见：相较于原始时间序列，本研究方法的数据更为密集，红框标注部分突出展示了本研究提出的时序重建方法的优越性，即可以在试验区生成有效的植被指数数据，明显且有效地填补时间序列中的缺失。这对于作物生育期遥感影像缺失的导致无法精准检测的情况尤为重要。尤其是在作物生长周期的主要拐点处，重建方法可以填补因云雾、噪声等因素影响而缺失的数据空白，这进一步导致了生长周期峰值迁移、斜率的上升以及物候指标日期的变化，这与图5~8展现出的结果相一致。

Figure 9.  Comparison of time series curves between original data and reconstructed data

精准的作物物候期提取是卫星遥感数据下游应用的关键任务^[24]。时间序列中的连续间隙会阻碍准确的物候提取^[25]，这在高分辨率时序影像中的表现更为明显。在原始地表反射率数据的年时间序列中，在作物的生育期，入海口处常伴随大量云雨天气，遥感影像在年度时间序列中经历了较大的、连续的和重复的缺失区域，这会导致物候指标的检索复杂化^[26]。原始影像中检索到的SOS和EOS中都存在大量的间隙，时间序列曲线也波动较大。在使用超分辨率方法重建这些空白影像后，作物生育期影像缺失的问题得到了缓解。

在今后的研究中可以选择生长季影像缺失更为严重的试验区进行探讨，以更广泛地验证其在物候提取中的作用。另外本研究仅验证了单幅遥感分辨率重建技术在卫星影像时间序列重建和物候提取方面的可行性以及提取结果和原始数据、MODIS数据提取结果之间的差异，更精准的数据验证需要结合PhenoCam^[27]等地面数据进行更进一步的探究。

本研究提出了一种基于生成式遥感影像重建技术的物候提取新方法，通过GAN网络将异源影像分辨率升级并实现了光谱迁移，生成有效的图像数据来填补空白以重建中分辨率卫星Sentinel-2数据的时间序列。并提出了一种轻量化的PGT-GAN网络，它提取多个深度特征，以捕获分辨率重建时的光谱和空间信息，这种方法在SSIM和PSNR指标上表现良好。这表明该方法能有效重建缺失区域的纹理和细节信息，很大程度地保持光谱范围的一致性。本研究提出的方法可以通过一次完整训练实现全年影像较高精度的分辨率重建，且不同于时空融合，并不需要与研究区精确对应的高分辨率参考影像进行影像重建。在Sentinel-2/Landsat-8真实数据集上对重建后的时序数据进行了实验，将结果与原始物候提取结果、插值提取结果和MODIS提取结果进行了对比分析，结果表明：本研究方法在各组结果中都展现了可靠的性能，显著提升了10 m分辨率数据的时间分辨率。提出的结合遥感影像分辨率重建的时间序列重建方法为产生高质量的地表数据提供了一种可行的方法，生成的更密集的时间序列数据不但可以检测到作物不同物候阶段的精准变化，也可以为更高时空分辨率陆地表面物候提取提供支持。