pyDMS一种热红外的超分辨率方法-源码解读与疑难点记录

前言

pyDMS是利用决策树方法将低分辨率影像升采样为高分辨率影像的python package。由Radoslaw Guzinski and contributors在2019年编写。原算法为2012 Gao Feng在Remote sensing上的文章。
本文目的：做项目时查到了这个包，有很多文献都使用了该方法，但无奈本人英语奇懒，原文算法看了三遍仍有很多细节不清楚，因此翻出源码并结合文章来一点一点扣细节，对疑惑和认为重要的点进行记录。如有勘误和表述不恰之处，欢迎指正。
原文链接：

[Gao2012] Gao, F., Kustas, W. P., & Anderson, M. C. (2012). A Data Mining Approach for Sharpening Thermal Satellite Imagery over Land. Remote Sensing, 4(11), 3287–3319. https://doi.org/10.3390/rs4113287
[Guzinski2019] Guzinski, R., & Nieto, H. (2019). Evaluating the feasibility of using Sentinel-2 and Sentinel-3 satellites for high-resolution evapotranspiration estimations. Remote Sensing of Environment, 221, 157–172. https://doi.org/10.1016/j.rse.2018.11.019

项目源码链接：pyDMS

要点1：

文章提到TsHARP有个前提，假设T-NDVI（热红外与NDVI的关系）在不同空间分辨率下是恒定的，因此低分辨的也可以用于高分辨图像；
反射率的聚合（应该是降采样，原文是aggregation）是线性的，但温度的聚合不是线性的（参考斯忒藩-玻尔兹曼定律，辐射出射度与温度的四次方成正比），因此混合像元的T-R（热红外与反射率）与T-NDVI可能随着空间分辨率变化。当然，均质像元的T-R关系是可以迁用的。因此引出一个问题——怎么确定均质像元呢？
文中提出了变异系数（Kustas et al.）。文中给出寻找的均质像元样本定义如下：如果在高分辨率上反射率相似（similar），则在低分辨率上这个差异是趋于0的，这类像元定义为“均质像元”（homogeneous pixels）。
变异系数计算如下：

这个公式是我认为是有误解的（最起码我第一次读的时候是懵逼的），准确来说，这里的$\sigma 和\mu$是以低分辨率像元大小为窗口，求取对应范围内高分辨率像元的均值和方差。
低分辨率由多少行列，就有多少个$\sigma 和\mu$。高分辨率有多少个波段，就有几层的$\sigma 和\mu$，然后对每个像元求变异系数Cv。
该部分源码如下：

def _resampleHighResToLowRes(bandData_HR, ySize_LR, yRes_LR, yRes_HR, xSize_LR, xRes_LR, xRes_HR,
                             gt_HR, gt_LR):
    aggregatedMean = np.zeros((ySize_LR, xSize_LR))
    aggregatedStd = np.zeros((ySize_LR, xSize_LR))
    for yPix_LR in range(ySize_LR):
        yPos_LR_min = gt_LR[3] + yPix_LR*yRes_LR
        yPix_HR_min = int(round(max(0, gt_HR[3] - yPos_LR_min) / yRes_HR))
        yPix_HR_max = int(round(max(0, gt_HR[3] - (yPos_LR_min + yRes_LR)) / yRes_HR))
        for xPix_LR in range(xSize_LR):
            xPos_LR_min = gt_LR[0] + xPix_LR*xRes_LR
            xPix_HR_min = int(round(max(0, xPos_LR_min - gt_HR[0]) / xRes_HR))
            xPix_HR_max = int(round(max(0, xPos_LR_min + xRes_LR - gt_HR[0]) / xRes_HR))
            aggregatedMean[yPix_LR, xPix_LR] =\
                np.nanmean(bandData_HR[yPix_HR_min:yPix_HR_max, xPix_HR_min:xPix_HR_max])
            aggregatedStd[yPix_LR, xPix_LR] =\
                np.nanstd(bandData_HR[yPix_HR_min:yPix_HR_max, xPix_HR_min:xPix_HR_max])

    return aggregatedMean, aggregatedStd
-------------------------------------------------------------------------------------------
resMean, resStd = utils.resampleHighResToLowRes(scene_HR, subsetScene_LR)
resMean[resMean == 0] = 0.000001
resCV = np.sum(resStd/resMean, 2) / resMean.shape[2]
resCV[np.isnan(resCV)] = 1000

要点2：

全局模型与局部模型结合的结果输出。
文中用了两种策略来保障输出结果：

• 第一种是扩大采样窗口，实际扩大一倍的采样数量，从而保证相邻移动窗口共享部分样本。
  目的是为了消除边界效应。局部模型对局部的效果更好，但容易造成边界效应（窗口与窗口之间存在分割线），如下图所示。
  
• 第二种是使用全局和局部的加权结果作为最终输出。
  每次运行时，会在缓存中保留局部和全局的结果，然后分别计算每个像元的局部和全局结果的残差（计算流程和源码见下）及权重，对残差结果进行高斯滤波（目的也是为了消除边界效应）；最后的输出结果是逐像元的全局与局部的加权平均和。
  源码如下：

# windowed weight
 ww = (1/windowedResidual)**2/((1/windowedResidual)**2 + (1/fullResidual)**2)
 # full weight
 fw = 1 - ww
 outData = outWindowData*ww + outFullData*fw

---

残差计算：
首先提取高分辨率与低分辨率重叠区域，然后判断低分辨率是否存在QA数据，如果存在将非good pixls设为nan值。
然后判断是否为温度数据，如果是需要转为辐亮度数据（因为辐射出射度与温度的四次方成正比，所以温度的降采样不是线性关系）源码是直接四次幂。
最后对锐化的高分辨率图像进行均值降采样，并用低分辨率影像减去高分辨率影像。

 # Then resample high res scene to low res pixel size
        if self.disaggregatingTemperature:
            # When working with tempratures they should be converted to
            # radiance values before aggregating to be physically accurate.
            radianceScene = utils.saveImg(downscaledScene.GetRasterBand(1).ReadAsArray()**4,  # 注意这里，四次幂
                                          downscaledScene.GetGeoTransform(),
                                          downscaledScene.GetProjection(),
                                          "MEM",
                                          noDataValue=np.nan)
            resMean, _ = utils.resampleHighResToLowRes(radianceScene,
                                                       subsetScene_LR)
            # Find the residual (difference) between the two)
            residual_LR = data_LR**4 - resMean[:, :, 0]  # 还有这里，都是四次幂
        else:
            resMean, _ = utils.resampleHighResToLowRes(downscaledScene,
                                                       subsetScene_LR)
            # Find the residual (difference) between the two
            residual_LR = data_LR - resMean[:, :, 0]
        # Smooth the residual and resample to high resolution
        residual = utils.binomialSmoother(residual_LR)  # 高斯滤波

权重计算公式如下：
c表示这个数据来自低分辨率影像，i应该是1或者2，代表某一像元的残差数量（局部或全局）,r代表残差，w为某一像元在局部模型或全局模型下的权重。