基于解混策略的时空融合模型稳定性分析

背景及科学问题

地表温度（Land Surface Temperature，LST）是地表能量和水平衡的重要组成部分，能够提供从局部到全球表面平衡状态的时空变化信息^[1]^。目前，地表温度已被广泛应用于蒸散发估算、城市热岛监测、局部气候变化、植被监测、农业监测以及森林火灾监测等领域^[2–5]^。然而，由于地表热辐射能量较低，加之卫星传感器受限于成本和硬件技术，单一卫星传感器获取的遥感数据存在“时空矛盾”，即难以同时满足高空间分辨率和高时间分辨率的需求^[6,7]^，这在很大程度上限制了地表温度数据的应用。

针对如何获取高时空分辨率LST这一问题，众多学者提出了不同的LST降尺度方法。这些方法按照算法理论可分为基于尺度因子的降尺度方法和基于时空融合的降尺度方法。其中，时空融合方法因其能够有效降低对辅助数据的依赖，且成本低、便捷高效，被广泛应用于地表温度降尺度^[8]^。根据实现原理，时空融合方法又可进一步分为基于权重函数的方法、基于学习的方法以及基于解混的方法。

基于权重函数的方法是假定地表反射率的时间变化在空间尺度上的一致性，通过权重函数组合所有输入影像的信息来预测目标影像高空间分辨率的像素值。Gao等人^[9]^最早提出时空自适应反射融合模型（Spatial and Temporal Adaptive Reflectance Fusion Model，STARFM），该模型假定粗分辨率影像中的像素为纯净像元，引入移动窗口和相似像元的概念，从光谱差异性、时间差异性以及空间距离来计算权重。对于复杂地貌中土地类型发生变化的区域，STARFM方法无法准确预测出未知时相真实的地物类型，导致融合结果产生较大误差。针对STARFM方法不足，Hilker等^[10]^提出了一种新的时空自适应反射率变化检测算法（Spatial Temporal Adaptive Algorithm for Mapping Reflectance Change，STAARCH），使用缨帽变换检测反射率变化，提高从低分辨率影像上检测出土地覆盖时空变化的能力。但STAARCH 中的变化检测仅适用于植被地表，因此Zhu等^[11]^提出的增强型时空自适应反射融合模型（Enhanced Spatial and Temporal Adaptive Reference Fusion Model，ESTARFM），引入转换系数和线性光谱解混理论，用光谱相关系数来代替空间距离，保留更多的空间细节的同时提高了在复杂地貌中的预测精度。针对ESTARFM 方法预测有形状变化的地物覆被随时间变化的准确性，Zhao等^[12]^人提出了鲁棒自适应时空融合模型（Robust Adaptive Spatial and Temporal Fusion Model ，RASTFM），由基于非局部线性回归的加权平均模块和基于非局部线性回归的图像超分辨率模块组成，分别用于预测地表非形状变化（包括物候变化和没有形状变化的土地覆盖变化）和形状变化（具有形状变化的土地覆盖变化）。STARFM系列模型参数简单，计算效率高，但在地物覆盖复杂、不同时相地物类型变化的预测仍具有局限性。

基于学习的方法通过机器学习算法用非线性的方式对高、低分辨率影像对建模，预测出高时空分辨率影像。Huang等^[13]^提出基于稀疏表示的时空反射融合模型（SParse-Representation-Based SpatioTemporal Reflectance Fusion Model，SPSTFM），利用字典对学习建立低空间分辨率影像和高空间分辨率影像的对应反射率变化关系，并通过时间加权预测目标高分辨率影像。针对字典对中扰动问题，Wu等^[14]^提出误差约束正则化的半耦合字典学习（Error-Bound-Regularized Semi-Coupled Dictionary Learning，EBSCDL），利用误差约束的正则化方法解决字典扰动，构建优化的半耦合字典解决低空间分辨率影像和高分辨率影像之间差异性，提高融合精度。然而上述基于稀疏表示的融合方法需要人为设计字典基元，在算法实现过程中，字典学习、稀疏编码和影像重建等步骤是分离的，提升算法的不稳定性和复杂性。对此，Song等^[15]^提出深度卷积神经网络时空融合模型（Spatiotemporal Fusion Method Based on Deep Convolutional Neural Networks，STFDCNN），结合卷积神经网络、非线性映射模型设计了一个双重卷积神经网络，实现了自动提取影像特征并提高预测精度，有效解决稀疏表示方法类的缺陷。基于学习的方法可以捕获到更多的地表空间细节，适用于异质区域等，但该类方法过度依赖训练样本和模型参数，具有较高的时间复杂性。

基于解混的方法很好地弥补了上述方法的不足。此类方法基于线性光谱混合理论，从高空间分辨率影像中提取地物类别和丰度，以分解低空间分辨率像元并获得类别的光谱值^[16]^。其实现过程通常包括端元选取、丰度计算和解混等步骤。最早由Zhukov等^[17]^提出的解混方法假设同类地物反射率相等，在移动窗口内对低空间分辨率影像进行解混，并将结果分配给未知时相的高空间分辨率影像。然而，该方法在解混过程中可能产生较大误差，且其假设在地物类型随时间变化的区域并不成立。对此，Zurita-Milla等^[18]^在线性解混过程中引入约束条件，用以处理解混结果中的负值和异常值问题；Wu等^[19]^提出了时空数据融合模型（STDFA），通过引入时间变化信息提高预测结果的准确性；针对STDFA算法未考虑影像端元反射率的空间异质性，Zhang等^[20]^结合多尺度分割算法和ISODATA算法生成分类图，并利用移动窗口方法对低空间分辨率影像进行解混，最终引入时间权重概念预测未知时相影像；为了提升模型在地物类型变化区域的预测准确性，Huang等^[21]^提出了基于解混的时空反射融合模型（U-STFM），该模型假设地物像元反射率在同质变化单元（homogeneous change regions，HCRs）上的变化率相等，并利用多尺度分割算法生成的超像素进行解混，从而提高了不同时相遥感影像在地表类型空间变化中的融合效果。

由于温度解混的复杂性和病态矩阵解混的准确性，解混模型稳定性受到很大程度的影响，其主要体现在以下三个方面：首先是HCRs分类方式对模型稳定性影响，由于不同空间分辨率像素之间的温度变化具有差异性，影响从高分辨率影像提取HCRs的准确性。HCRs分类方式主要有聚类方法和面向对象分类方法，面向对象分类方法是以影像分割后的对象为分类单元来提取影像光谱信息、纹理信息以及几何信息，聚类方法是将空间内数据按照某种相似度聚集为不同类别；其次是解混次数对模型稳定性影响，解混类方法常采用最小二乘算法来求解线性混合方程组，是一个病态矩阵求解问题，解混次数的叠加也对模型的稳定性有一定的影响。最后是HCRs数量对模型稳定性影响，HCRs数量越多越能细致刻画地表变化，同时也增加了方程解的不确定性。因此，本文从HCRs分类方式、解混次数以及HCRs数量三方面来分析解混解混策略对模型稳定性影响。

本文以粤港澳大湾区部分区域为研究区，选取Landsat7 LST产品和MODIS LST产品，基于U-STFM算法从HCRs的分类方式、解混次数以及HCRs的数量三方面分析对模型稳定性的影响，然后将改进后的方法与经典的STARFM、ESTARFM时空融合方法对比，验证该方法在LST数据中预测的有效性和稳定性。

研究方法介绍

具体实验设计如图2所示：

本文采用Landsat7 LST和MODIS LST数据从HCRs分类方式、解混次数以及HCRs数量三方面来分析解混模型的稳定性，并将改进后的模型与经典模型对比LST融合后的精度。在U-STFM方法上，首先采用ISODATA算法[26]、K-Means算法[27]以及多尺度分割算法[28]分别获取HCRs分类图，选取2001年11月20日为预测日期，排列组合共12组实验日期对，通过对比基于HCRs不同分类方式得到的融合结果分析其对模型稳定性的影响；其次基于空间像元分解的时空融合方法原理[17]，在分析HCRs分类方式结果的基础上，控制变量HCRs的类别数，分别对MODIS单期影像、MODIS差分影像、MODIS影像变化趋势比率进行解混，选取表1中的8个预测日期排列组合共169组影像，通过对比基于解混次数得到的融合结果分析其对模型稳定性的影响；然后在上述分析的结果基础上，同样选取表1中的8个预测日期排列组合共169组影像，通过对比不同HCRs数量得到的融合结果分析其对模型稳定性的影响。最后选取最佳解混策略后的模型与两个经典模型对比分析，测试该模型在LST数据融合效果。其文中各模型融合结果绝对误差分布图和各模型影像对比图均来源于20000914-20011120-20011222日期组。

U-STFM时空融合方法

本次实验选取解混方法中最具有代表性的U-STFM算法，将预测日期当天（tp）的MODIS LST数据以及前后两期（t0、te）的MODIS LST、Landsat7 LST数据作为该模型的输入，以预测目标日期Landsat7 LST数据。该模型基于线性光谱混合模型，忽略不同传感器数据之间的误差，其低空间分辨率像素等于HCRs光谱值与其覆盖率（fractional）的线性组合。解混过程计算如下：

精度评价

本文采用定性和定量评价指标来评价模型预测地表温度的效果。通过对比分析预测影像和真实影像可视化效果，实现对模型融合精度的定性评估。其中定量评价采用峰值信噪比PSNR（Peak Signal to Noise Ratio）、相关系数CC（Correlation Coefficient）、均方根误差RMSE（Root-Mean-Square Error）以及平均绝对误差MAE（Mean Absolute Error）。PSNR是一种全参考的影像质量评价指标，PSNR越高表示影像融合效果越好；CC、RMSE以及MAE表示预测影像和参考影像之间的相关性和差异性，CC值有效范围在区间[-1,1]之间，越接近1表示融合结果越好，而RMSE、MAE值越小代表融合影像质量越高。上述定量评价指标均通过调用scikit-learn模块对应的评价指标函数来实现。

试验结果

HCRs分类方式对模型稳定性影响

HCRs选取是解混模型第一步，其划分的准确性对模型融合效果至关重要。因此本文选取多尺度分割算法、K-means算法以及ISODATA算法来划分HCRs，其中多尺度分割算法是一种自下而上的区域合并技术，通过识别遥感影像的光谱、形状等特征差异，从空间上划分出不同的区域；k-mean是基于距离划分的聚类算法[29]，以数据对象之间距离作为相似度准则，通过迭代更新将数据划分为不同的簇；ISODATA是动态聚类分析算法，在K-means算法的基础上，增加“合并”和“分裂”操作，预先设定聚类中心和聚类数，通过定义相似度准则函数将全部样本调整完后再重新计算样本均值作为新的聚类中心[30]。

图3所示分别为HCRs三种分类方法得出的融合结果的定量评价结果。可以看出，ISODATA算法相对于K-means算法和多尺度分割算法更能准确划分HCRs，更适合捕捉地物时间序列上的LST变化，其融合结果具有较高的PNSR以及较低的RMSE和MAE；三种分类算法融合结果CC均在85%以上，基于多尺度分割算法的预测影像与参考影像相关性最高，聚类算法对HCRs划分受数据源质量的影响较大。图4所示为基于三种HCRs分类算法的融合结果与参考影像进行差值处理并取绝对值，以更直观的方式来分析不同分类算法对模型精度的影像。从空间分布来看，基于三种HCRs分类算法的融合结果的绝对误差大部分分布在5K以内；其中基于ISODATA算法（图4（c））融合结果的绝对误差普遍小于2K，表现出较好的结果。总体来看，对于时刻变化的地表温度数据，ISODATA算法更能有效的刻画出LST在时空变化中的HCRs，因此最终选取其作为HCRs的分类方式。

解混次数对模型稳定性影响

解混方程是解混类模型的核心组成部分，而HCRs、丰度以及MODIS像元组成的解混矩阵是一种病态矩阵，其解受异常值影响较大。因此本节将探讨解混次数的累加对模型稳定性的影响，采用ISODATA算法划分HCRs基础上将最小二乘算法求解线性光谱解混应用到影像跨尺度融合的不同位置，选取RMSE指标定量分析对融合结果精度的影响，结果如图5所示。图中解混1次是指直接对MODIS影像变化趋势比率进行解混，解混2次是指分别对俩组MODIS差分影像解混，解混3次是指分别对三期MODIS单期影像解混。从各评价指标箱线图可以明显看出解混3次时模型融合效果明显低于其他两种解混次数；基于解混1次和基于解混2次的模型融合效果相似，从箱线图中值来看，基于解混1次模型融合结果的CC和RMSE较高、MAE和PNSR较低，但整体上基于解混1次模型融合结果偏差较小，更稳定。这是因为随着解混方程次数的累加其解的奇异值也会随之增加，而解混2次前MODIS影像做了差分处理，削弱了奇异值的数量，因此会和解混1次时的效果相似。图6为基于上述三种解混策略得到的融合影像与参考影像的绝对误差分布图，可以看出解混1次和解混2次时模型预测误差大部分都在在5k以下，其中2k-5k解混2次时模型预测误差分布相对较多，而解混3次时模型预测误差大部分分布在5k-20k，模型误差较大。上述结果进一步揭示了解混算法使用次数越多，模型误差越大，也反映出基于最小二乘的解混算法的不稳定性。

HCRs类别个数对模型稳定性影响

理论上来讲HCRs划分越细就会获得越多的细节，但当HCRs数量大于MODIS像元数量时，会导致解混方程解不稳定，而HCRs数量也会影响方程组解的准确性，因此选取适当HCRs数量对模型的性能至关重要。本节在采用ISODATA算法划分HCRs和解混1次的基础上，探讨HCRs类别数对融合结果的影响。图7所示为169组模型在不同HCRs类别数上的预测影像评价指标中位数散点图，可以看出图中MAE与HCRs类别数量有很强的相关性，随着HCRs类别数量增加MAE呈现下降趋势；PSNR与RMSE在HCRs类别数为104以内时分别呈现上升、下降趋势，之后逐渐趋于稳定；CC值均在0.92以上，在 HCRs类别数为100以内时值相对较高，之后有下降趋势。图5为选取基于HCRs类别10、类别104和类别207时模型融合LST结果绝对误差的可视化结果，可以看出基于这三种类别的融合影像与参考影像的绝对误差普遍分布在5K以内，其中绝对误差在2K-5K分布最多的为基于HCRs10类模型，基于 HCRs104类模型以及HCRs207类模型总体分布较相似，前者较集中。总体上来看当HCRs类别数量在104时模型融合效果较好。

对比其他模型

基于上述实验结果，先对Landsat7 LST数据采用ISODATA算法进行聚类，调参后得到HCRs104类分类图，之后直接对MODIS影像变化趋势比率进行解混，然后基于HCRs时间变化率一致性假设预测出未知高空间分辨率地表温度。为了将原模型HCRs划分方式区分开，我们将其记为iso_USTFM。图9是iso_USTFM（图9（c））与两个经典的融合模型STARFM、ESTARFM预测生成的LST分布图（图9（a）、图9（b））以及参考影像Landsat7 LST分布图（图9（d）），图中红色表示高温地物，由房屋、公路等非植被区域控制；绿色表示低温地物，由草地、森林等植被区域控制。从图中可以看出，由iso_USTFM预测生成的LST分布图与Landsat7的LST产品分布图最为接近，且展现更多的纹理信息，但在城镇区错误的预测出LST高温值，而ESTARFM却相反，将与植被区域接近的非植被区域错误的预测为低温值；STARFM预测的LST与Landsat7的LST产品误差最大，LST分布出现大量与现实不符的低温和高温极值。图10为三种模型预测生成的LST与Landsat7的LST产品的散点图，从数据分布来看，三组的数据点都接近对角线，其中（图10（a））有大量的异常值，而（图10（c））相关性最高。从评价指标来看iso_USTFM预测出的LST与Landsat7的LST产品决定系数为0.874，RMSE为1.0948K，均优于STARFM、ESTARFM预测生成的LST，融合效果较好。综上可以得出，iso_USTFM方法与STARFM方法、ESTARFM方法预测生成的地表温度相比，预测效果最好。

核心结论

本文从HCRs分类方式、解混次数以及HCRs数量三方面分析了解混策略对模型稳定性的影响。首先，将USTFM方法应用于生成高时空分辨率LST，并基于LST特性分析模型中的解混策略；然后，应用iso_USTFM方法融合Landsat7 LST产品和MODIS LST产品，生成具有MODIS时间分辨率的30m空间分辨率LST影像；最后，与经典的时空融合模型STARFM和ESTARFM对比预测出的LST影像效果，得出以下结论：

（1）ISODATA聚类算法相较于多尺度分割算法和K-means算法，更能有效刻画LST在时空变化中的HCRs。

（2）解混算法使用次数越少，模型预测误差越小，稳定性越高。

（3）随着HCRs类别数量的增加，模型预测效果逐渐提升，并在类别数为104时趋于稳定。

（4）在地表温度融合方面，iso_USTFM算法生成的预测影像纹理清晰、精度较高，相较于STARFM和ESTARFM算法更适合预测LST。

此外，该解混类融合方法受限于遥感影像自身分辨率的限制。当两个传感器影像的空间分辨率差距较大时，无法详细预测高分辨率影像像元中的空间细节。iso_USTFM中HCRs占比通过其在MODIS像元中的面积确定，这种简单的线性解混可能无法有效刻画LST的变化，且其基于最小二乘算法的解混过程易受异常值影响，导致模型解不准。因此，未来研究将考虑结合非线性解混策略和深度学习方法优化模型，以提升模型的稳定性和鲁棒性。

‍