混合像元分解的基础是线性光谱混合模型，在忽略不同传感器数据之间的地理配准误差和大气校正误差的条件下，混合像元的光谱值等于该像元内地物类别的光谱值与其丰度的线性组合(Settle和Drake，1993)(式(2))。然而，由于空间异质性的存在和最小二乘算法理论的需要，式(2)需要在一定大小的窗口内进行，但过大的窗口会削弱空间的异质性。同时，窗口的使用既可以消除中、高分辨率影像之间因配准误差带来的影响，又可以体现同类地物的空间差异性(Wu等，2012；Zhang等，2013)。因此，选择合适的矩形窗口进行混合像元分解是必要的。为选择一个合适窗口，设置不同大小的窗口进行混合像元分解，对各窗口下降尺度的数据与真实数据进行比较，通过定量精度指标选择合适的窗口。在一定大小窗口下，利用式(1)求得窗口内每个混合像元内各类别丰度，组成到n²行k列的丰度矩阵A(n²×k)(n²＞k)(n表示矩形窗口大小，k表示窗口内类别数)，利用最小二乘法(式(3))求解窗口内各类别的光谱值(Busetto等，2008；Wu等，2012；Zhang等，2013)。然后，把求得的值依照类别类型对应到窗口内相应像元的位置上，获得分解后的降尺度数据。

$ R\left( {i,t} \right) = \sum\limits_{c = 0}^k {{f_c}\left( {i,c} \right) \times \bar r\left( {c,t} \right) + \xi \left( {i,t} \right)} $ (2)

$ \left[ {\begin{array}{*{20}{c}} {R\left( {1,t} \right)}\\ \cdots \\ {R\left( {i,t} \right)}\\ \cdots \\ {R\left( {n,t} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{f_c}\left( {1,1} \right)}& \cdots &{{f_c}\left( {1,c} \right)}& \cdots &{{f_c}\left( {1,k} \right)}\\ \cdots &{}& \cdots &{}& \cdots \\ {{f_c}\left( {i,1} \right)}& \cdots &{{f_c}\left( {i,c} \right)}& \cdots &{{f_c}\left( {i,k} \right)}\\ \cdots &{}& \cdots &{}& \cdots \\ {{f_c}\left( {n,1} \right)}& \cdots &{{f_c}\left( {n,c} \right)}& \cdots &{{f_c}\left( {n,k} \right)} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\bar r\left( {1,t} \right)}\\ \cdots \\ {\bar r\left( {c,t} \right)}\\ \cdots \\ {\bar r\left( {k,t} \right)} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {\xi \left( {1,t} \right)}\\ \cdots \\ {\xi \left( {c,t} \right)}\\ \cdots \\ {\xi \left( {k,t} \right)} \end{array}} \right] $ (3)

约束条件： $ \sum\limits_{c = 0}^k {{f_c}\left({i，c} \right)= 1;0 < } \bar r\left({c，t} \right)< 1 $

式中，R(i，t)为t时期i位置低分辨率像元的反射率，f_c(i，c)为i位置低分辨率像元内类别c占该像元的面积比，(c，t)为t时期类别c的平均反射率，ξ(i，t)为残差，k为窗口内类别数。

CDSTARFM利用t₀和t_k两期低分辨率影像降尺度的数据和t_k时期中分辨率遥感影像预测t₀时期中分辨率遥感数据。STARFM算法实现的关键在于如何权重空间信息，且这个权重算法可以自适应于不同的区域(Gao等，2006)。Gao等人(2006)通过在一定窗口内筛选与中心像元光谱相似的相似性像元，利用相似性像元的光谱差异性，时间差异性和与中心像元的相对距离来加权计算预测期的中心像元值。由于CDSTARFM继承了STARFM算法，其实现过程也是通过加权窗口内相似性像元计算预测期的中心像元值。

窗口内中心像元的相似性像元由t_k时期中分辨率影像各波段的阈值θ_b确定，窗口内与中心像元差值的绝对值小于各波段θ_b的像元作为相似性像元。

$ {\theta _b} = \sqrt {\frac{1}{N}{{\sum\limits_{i = 1}^N {\left( {{x_i} - \mu } \right)} }^2}} /m $ (4)

式中，θ_b表示窗口内b波段的阈值，N为窗口内中分辨率像元个数，x_i为i位置的像元反射率，μ为窗口内像元反射率均值，m为类别数(Gao等，2006；Zhu等，2010)。

相似像元的权重大小W_ijk有3个指标来衡量：(1)光谱的差异性S_ijk；(2)时间的差异性T_ijk；(3)相似性像元与中心像元的相对距离D_ijk；

$ {S_{ijk}} = \left| {L\left( {{x_i},{y_i},{t_k}} \right) - M\left( {{x_i},{y_i},{t_k}} \right)} \right| $ (5)

$ {T_{ijk}} = \left| {M\left( {{x_i},{y_i},{t_k}} \right) - M\left( {{x_i},{y_i},{t_0}} \right)} \right| $ (6)

$ {D_{ijk}} = 1 + \sqrt {{{\left( {{x_{w/2}} - {x_i}} \right)}^2} + {{\left( {{y_{w/2}} - {y_j}} \right)}^2}} /A $ (7)

式中，A为表示空间距离相对于光谱差异性和时间差异性重要性的常量。A=750 m(Gao等，2006)

利用式(8)(9)加权窗口内相似像元的反射率计算预测期的中心像元反射率，生成最终的融合影像。

$ \begin{array}{l}L\left({{x_{w/2}}，{y_{w/2}}，{t_0}} \right)= \\\sum\limits_{i = 1}^w {\sum\limits_{j = 1}^w {\sum\limits_{k = 1}^w {{W_{ijk}} \times \left({M\left({{x_i}，{y_i}，{t_0}} \right)+ L\left({{x_i}，{y_i}，{t_k}} \right)- M\left({{x_i}，{y_i}，{t_k}} \right)} \right)} } } \end{array} $ (8)

$ \begin{array}{l}{W_{ijk}} = \left({\frac{1}{{{S_{ijk}} \times {T_{ijk}} \times {D_{ijk}}}}} \right)/\\\sum\limits_{i = 1}^w {\sum\limits_{j = 1}^w {\sum\limits_{k = 1}^w {\left({\frac{1}{{{S_{ijk}} \times {T_{ijk}} \times {D_{ijk}}}}} \right)} } } \end{array} $ (9)

式(5)—(9)中，w为窗口大小，L(x_i，y_j，t_k)和M(x_i，y_j，t_k)分别为t_k时期给定位置(x_i，y_j)的中分辨率数据和降尺度数据的像元值，(x_w/2，y_w/2)为窗口的中心像元，光谱差异性S_ijk值越小说明给定位置与邻近像元的光谱相似度越高，赋予较高的权重；时间差异性T_ijk值越小说明该段时间内光谱变化越小，赋予较高的权重；相对距离D_ijk值越小，赋予较高的权重(Gao等，2006)。

首先，利用2014-08-19的L and sat 8多光谱数据进行ISODATA聚类，子类合并后类别为10类，结合聚类数据和最小二乘算法分别在窗口大小w(MODIS像元)为5×5、7×7、11×11、15×15、21×21、31×31和41×41对2014-08-19的MODIS数据进行像元分解。利用分解后的降尺度数据与真实L and sat 8数据之间的相关系数r、偏差Bias、均方根误差RMSE和ERGAS(the Erreur Relative Globale Adimensionalle de Synthèse)定量指标对降尺度数据进行精度评价，选出最适的窗口大小。其次，在该最佳窗口下对2014-09-04的MODIS数据进行像元分解降尺度。最后，利用CDSTARFM算法和2014-08-19和2014-09-04两期最佳窗口下的降尺度数据以及基期2014-08-19 L and sat 8数据进行数据融合。CDSTARFM 算法的搜寻窗口w(L and sat OLI像元)设为7×7、11×11、31×31、61×61、101×101和151×151。为了便于该算法与STARFM算法比较，STARFM算法利用与CDSTARFM算法相同的参数进行数据融合。通过融合数据与真实数据之间的r、Bias、RMSE和ERGAS等指标对融合数据进行精度评价(Amorós-López等，2013；Gevaert和García-Haro，2015)。

$r = \frac{{\sum\limits_{j = 1}^N {\left( {{x_j} - \bar x} \right)\left( {{y_j} - \bar y} \right)} }}{{\sqrt {\sum\limits_{j = 1}^N {{{\left( {{x_j} - \bar x} \right)}^2}} \sum\limits_{j = 1}^N {{{\left( {{y_j} - \bar x} \right)}^1}} } }}$ (10)

$ Bais = \frac{1}{N}\sum\limits_{j = 1}^N {\left({{x_j} - {y_i}} \right)} $ (11)

$ RMSE = \sqrt {\frac{{{{\sum\limits_{j = 1}^N {\left({{x_j} - {y_i}} \right)} }^2}}}{N}} $ (12)

式中，N表示像元个数，j表示融合数据和真实数据的像元位置，x_j和y_j分别表示j位置融合数据和真实数据的像元值，$ {\bar x} $和$ {\bar y} $分别表示融合数据和真实数据的平均值。

$ ERGAS = 100\frac{H}{L}\sqrt {\frac{1}{{{N_{ban}}}}\sum\limits_{i = 1}^{{N_{ban}}} {{{\left({RMS{E_i}/{M_i}} \right)}^2}} } $ (13)

式中，H表示高分辨率影像的像元大小(30 m)，L表示低分辨影像的像元大小(480 m)，N_ban表示波段的数量，RMSE_i表示i波段融合数据与真实数据对应波段的均方根误差，M_i表示真实L and sat 8数据i波段的平均值。ERGAS值越小说明融合数据与真实数据在空间上的差异就越小，融合数据的质量就越高(Amorós-López等，2013；Gevaert和García-Haro，2015)。

像元降尺度算法、STARFM算法和CDSTARFM算法的最佳融合影像分别见图 3(c)—(e)。3种方法融合的影像与2014-08-19和2014-09-04真实L and sat 8影像进行目视比较，图 3(c)—(e)与图 3(b)相比较可以看出，3种方法对面积较大的地物(水体、居民区和较宽的道路)均可清晰地识别。但从放大的局部细节图(图 3(h)—(j))来看，在较破碎的区域，像元降尺度算法的融合影像(图 3(h))有“图斑”，不能识别较小的地物及地物内部的纹理信息，而图 3(i)(j)能够反映出细小地物和纹理特征；由于STARFM算法使用的是重采样的MODIS数据，使得同一MODIS像元内重采样的像元之间具有很强的均质性，相邻的MODIS像元之间在光谱上具有一定的差异，尤其在破碎区域其差异性更大，这是导致融合影像出现“MODIS像元边界”(图 3(i)中黑色框内所示)的主要原因，而CDSTARFM算法的融合影像可以有效消除MODIS边界。图 3(g)中黄色椭圆框是图 3(f)发生变化的地物，但3种融合方法均没有捕捉到该变化信息，这是因为用于融合的基期L and sat 8影像(图 3(f))没有体现该变化信息(Gao等，2006)，由于CDSTARFM算法大部分继承了STARFM算法，所以也没有解决该问题；但Hilker等提出的STAARCH方法通过选择合适的基期L and sat数据解决了这个问题(Hilker等，2009a)。总体来说，CDSTARFM融合影像的目视效果更接近真实影像。

点击放大

下载原图

图3 3种方法在最佳窗口融合的2014-09-04影像与真实影像的比较(NIR-Red-Green组合)

Figure 3 Comparison of reference image and fusion images predicted by the downscaling method, STARFM, and CDSTARFM on September 4, 2014 at the best window size, respectively (NIR-Red-Green combination)