基于光谱与空间局部相关的遥感影像融合方法

摘要

本发明提供了一种基于光谱与空间局部相关的遥感影像融合方法，采用自适应分块局部优化的策略，将整幅影像进行分块，然后在每个影像块内分别建立观测模型和相关融合模型；在求解光谱与空间相关影像融合模型过程时利用非负最小二乘方法进行求解权系数矩阵，在遇到矩阵不可逆情况下进行异常的处理得到权系数矩阵；利用解求出的块权系数矩阵进行双线性内插获得每个像素的权系数，从而获得每个像素的模拟低分辨率全色影像值，将这个影像值带入观测模型就可以逐一像素的进行融合处理。本方法可以在提高原始多光谱影像空间分辨率的同时有效保持原有的光谱信息，并对现有的光学传感器，具有较强的自适应处理能力。

说明书

技术领域

本发明属于遥感影像处理数据融合技术领域，涉及一种基于光谱与空间局部相关的遥感 影像融合方法。

背景技术

随着对地观测技术的发展，大多数对地观测卫星如：Landsat-7,Spot5,IKONOS, QuickBird,GeoEye-1,和WorldView-2可以同时提供高光谱分辨率的多光谱影像和高空间分 辨率的全色影像。为了充分利用光谱和空间特性，遥感影像融合技术得到了迅速的发展。

遥感影像融合技术是指将高空间分辨率的全色影像与高光谱分辨率的多光谱影像的优 点集成在一起，生成一幅具有高空间和光谱分辨率的彩色影像的影像处理过程。

目前的遥感影像融合方法可以分为四类。以IHS、PCA、Gram-Schmidt Spectral、 HCS(Hyperspherical Color Sharpening)和偏微分替换为代表的投影替换类算法，利用全色波段 或其变体替换多光谱影像投影变换后得到的亮度分量，然后进行投影逆变换获得融合结果； 以Brovey、SVR(Synthetic Variable Ratio，合成变量比)为代表的相对光谱分布类算法，基于全 色影像与多光谱影像之间的线性组合假设，把全色影像信息按比例分布到多光谱影像的各个 波段；以Wavelets、curvelets、HPF(high-pass filtering)、multiscale toggle contrast operator、 Markov random field(MRF)models和多尺度光谱增益调制为代表的ARSIS类算法在不同尺度 上提取影像的空间信息和光谱信息进行影像融合。投影替换和相对光谱分布这两类算法能够 较好地保持空间分辨率，但都依赖于全色影像与多光谱影像的相关性，相关性越好融合结果 越好，且光谱信息的保持依赖于传感器的配置和影像信息本身；ARSIS类算法能够较好地保 持光谱特性，但是全色影像和多光谱影像空间信息的不一致容易造成空间细节信息的混叠， 导致空间分辨率的降低；混合模型算法则试图吸收前三类算法在空间信息保持或光谱信息保 持的优势，并引入Map估计、稀疏表达、压缩感知等新理论，该类算法研究处于起步阶段， 算法也比较复杂，离实用化还有一定距离。

总体而言，相对光谱分布类的代表性算法SVR算法具有明确的物理意义和理论基础， 在全色波段和多光谱波段（或部分波段）能够相互覆盖的条件下（如IKONOS、QuickBird）， 不仅能保持全色影像的空间细节，也能很好地保持多光谱影像的光谱信息。然而，对于不满 足光谱覆盖条件的影像，SVR会造成融合影像的光谱扭曲；另外，对于满足光谱覆盖条件的 影像，也可能会由于影像上地物类型分布不均衡，导致组合系数求解的精度不够，使合成的 全色波段（PAN）影像与实际不符，造成融合影像的光谱失真。如何克服偏色因素，使融合 影像在保持空间细节的同时又能很好地保持光谱特性是一个需要解决的问题。

发明内容

本发明为解决现有光学遥感传感器在影像融合过程中采用基于SVR方法存在的上述技 术问题，提供了一种基于光谱与空间局部相关的遥感影像融合方法（SSCSVR，Spectral and  Spatial Correlation-Based Synthetic Variable Ratio）。

本发明的技术方案为一种基于光谱与空间局部相关的遥感影像融合方法，包括以下步骤：

步骤1，对影像进行分块，得到多个影像块；

步骤2，对每个影像块分别求取块内的权系数矩阵，包括以下步骤，

步骤2.1，提取全色影像中每个像素的高频成分G_s；

步骤2.2，建立光谱与空间相关模型；

步骤2.3，利用非负最小二乘方法对步骤2.2中光谱与空间相关模型的权系数矩阵进行求解， 并进行异常处理；

步骤3，利用步骤2所得权系数矩阵进行双线性内插获得影像中每个像素的权系数进而获得合成的低分辨率全色影像Pan_Syn，最后对每个像素逐一进行融合，得到融合影像。

而且，步骤1中对影像进行分块采用自适应方式，包括将分块大小记为w×w，取 w＝5r+1，比值r＝l/h，其中h和l分别为原始的全色影像和要融合的多光谱影像的空间分 辨率。

而且，步骤2.1通过下式实现，

$G_s={\mathrm{Pan}}_{\mathrm{ori}}-G(x,y;\sigma )\otimes {\mathrm{Pan}}_{\mathrm{ori}}---\left(1\right)$

其中，高斯函数σ为尺度参数，Pan_ori是原始的全色影像灰度值。

而且，步骤2.2中建立光谱与空间相关模型如下式

式中，是待求的光谱成分权值，β是待求的空间相关成分的权值，MS_i是要融合的多 光谱影像中第i波段的灰度值，i的取值范围为1,2，...N，N为多光谱影像波段数；

光谱与空间相关模型的矩阵形式如下式所示

Pan_ori＝BX    (3)

其中为待求的权系数矩阵， B＝[MS₁,MS₂,...MS_N,G_s]为光谱与空间相关成分，矩阵维数为N+1。

而且，步骤2.3实现方式如下，

设为权系数矩阵X的最小二乘估计量，非负约束最小二乘原 理描述如下，

Minimize LSE＝(BX-Pan_ori)^T(BX-Pan_ori)subject to X≥0

LSE为非负约束的优化误差，X≥0表示对于所有的1≤k≤N+1都有x_k≥0，引入N+1维 未知的被动约束常量c＝[c₁,c₂,...c_N+1]^T＞0,c_k＞0，k的取值为1,2,...N+1，形成拉格朗日式J 如下

$J=\frac{1}{2}{(\mathrm{BX}-{\mathrm{Pan}}_{\mathrm{ori}})}^T(\mathrm{BX}-{\mathrm{Pan}}_{\mathrm{ori}})+\lambda (X-c)---\left(4\right)$

当推出

$B^{T}B{\hat{X}}_{\mathrm{NCLS}}-B^T{\mathrm{Pan}}_{\mathrm{ori}}+\lambda =0---\left(5\right)$

从公式(4)和(5)得出公式(6)

${\hat{X}}_{\mathrm{NLCS}}={\left(B^{T}B\right)}^{-1}{B}^T{P}_{{\mathrm{an}}_{\mathrm{ori}}}-{\left(B^{T}B\right)}^{-1}\lambda$

(6)

${=\hat{X}}_{\mathrm{LS}}-{\left(B^{T}B\right)}^{-1}\lambda$

和公式(7)

$\lambda =B^T({\mathrm{Pan}}_{\mathrm{ori}}-B{\hat{X}}_{\mathrm{NCLS}})---\left(7\right)$

其中，为非负约束优化解，拉格朗日系数λ＝(λ₁λ₂...λ_N+1)^T；

当实对称矩阵B^TB可逆时，采用快速非负约束最小二乘算法求解权系数矩阵，求解方式 为将估计量分离为活动集R和被动集P，活动集R包括估计量中所有负值和零值，被 动集P包括估计量中所有剩余的正值；快速非负约束最小二乘算法处理开始时被动集P 为空集，假设活动集R包括估计量中所有元素，通过公式(6)和(7)迭代处理调节活动集和 被动集，当活动集R为空集或者活动集R不为空集且所有拉格朗日系数λ_k，k∈R为负值时迭 代结束，得到权系数矩阵；

当实对称矩阵B^TB不可逆时，求解权系数矩阵采用下式，

其中，np为用于求解系数的像素个数，j的取值范围为1,2,...np，Pan_j是全色影像第 j个像素的灰度值，MS_i，j是多光谱影像第i波段第j个像素的灰度值。

而且，步骤3的具体实现操作如下，

利用双线性内插对步骤2.3求解出的块光谱成分权值进行内插，得到每个像素的然后 采用下式得到每个像素合成的低分辨率全色影像值Pan_Syn，

最后利用以下多光谱关系模型获得每个像素的融合值，

${\mathrm{Fused}}_i={\mathrm{MS}}_i\times \frac{{\mathrm{Pan}}_{\mathrm{ori}}}{{\mathrm{Pan}}_{\mathrm{Syn}}}$

其中，Fused_i是融合后影像中第i波段的灰度值，Pan_Syn是合成的低分辨率全色影像灰度 值。

本发明与现有技术的相比，改进主要体现在：首先，本发明建立了基于光谱与空间相关 的回归模型，使融合影像和原始多光谱影像之间的关系更加准确；其次，在建立回归模型时， 采用局部分块处理，从而更好的顾及局部的变化，使求的权值更能够适应影像的局部特征； 再有，求解回归模型时采用非负的最小二乘方法，改进了已有使用最小二乘方法的策略，使 求出回归值更有物理意义，得到的结果更加准确。可见，本发明方法利用光谱与空间相关特 性、局部优化和非负最小二乘求解的融合策略，可有效提高融合影像的质量，达到比较理想 的融合效果，光谱和空间细节保持性好；该方法的另一个大的优点就是对现有的光学传感器， 具有较高的自适应处理能力。

附图说明

图1是本发明的实施例流程图。

具体实施方式

以下结合附图和实施例详细说明本发明技术方案。

本发明技术方案可采用计算机软件技术实现自动运行流程。如图1所示，实施例的运行 流程包括以下步骤：

步骤1，将影像进行自适应分块，得到多个影像块，可记为block₁、block₂、block₃…， 便于局部处理。

步骤2，对每个影像块分别求取块内的权系数矩阵，对每个影像块处理过程包括以下步 骤，

步骤2.1，提取全色影像中每个像素的高频成分G_s；

步骤2.2，建立光谱与空间相关（SSC）模型；

步骤2.3，利用非负最小二乘方法对步骤2.2中光谱与空间相关模型的权系数矩阵进行求解， 并进行异常处理。

分别求取得到影像块block₁、block₂、block₃…相对应的权系数矩阵X₁、X₂、X₃…

步骤3，利用步骤2所得权系数矩阵进行双线性内插获得影像的每个像素的权系数获得合成的低分辨率影像Pan_Syn，最后对每个像素逐一进行融合，得到融合影像。

为便于实施参考起见，提供实施例涉及各操作的详细实现说明如下：

一、建立SVR多光谱遥感影像的关系模型；

该模型认为高分辨率全色影像和合成的低分辨率全色影像的比值反映了全色影像和多光 谱影像的细节差异，如果把这些差异按比例分布到多光谱影像的各个波段，就可以生成和全 色影像分辨率一致的高分辨率多光谱影像,从而提供了高空间分辨率融合多光谱影像与低空 间分辨率多光谱影像之间的对应关系

${\mathrm{Fused}}_i={\mathrm{MS}}_i\times \frac{{\mathrm{Pan}}_{\mathrm{ori}}}{{\mathrm{Pan}}_{\mathrm{Syn}}}---\left(9\right)$

式中，Fused_i是融合后影像中第i波段的灰度值，MS_i是要融合的多光谱影像（MS）中 第i波段的灰度值，Pan_ori是原始的全色影像(PAN)灰度值，Pan_Syn是合成的低分辨率全色影像 灰度值。Pan_Syn由下式(10)来计算得到。N为多光谱影像的波段数，i的取值范围为1,2,...N。

二、建立低分辨率全色影像合成模型；

传统SVR方法认为低分辨率全色影像能够利用多光谱影像合成：

式中，Pan_Syn是合成的全色影像灰度值，MS_i是要融合的多光谱影像中第i波段的灰度值， 是第i波段的权系数值。

根据全色影像和多光谱影像的线性组合假设，全色影像合成的权系数可通过式(11)来 解得。

Munechika等人利用城区、土壤、水、草地和树木五类地物来求取权系数

Zhang Yun对SVR的方法进行进一步研究，明确了合成全色影像的多光谱波段应该是全 色波段覆盖的波段，而不一定是多光谱影像的全部波段。同时，Zhang Yun认为求合成系数不能仅仅由指定的五类地物来求解，而应该利用影像的全部整体信息来求解。为此，Zhang Yun 的改进算法采用了直方图归一化、融合后影像直方图重新拉伸等处理，使SVR成为一类具有 独特意义的实用算法。然而对于不满足光谱覆盖条件的影像，Zhang Yun SVR会造成融合影 像的光谱扭曲。另外，对于满足光谱覆盖条件的影像，也可能会由于影像上地物类型分布的 不均衡，导致组合系数求解的精度不够，使合成的PAN影像与实际不符，造成融合影像的光 谱失真。本发明利用全部参与融合的多光谱波段参与低分辨率全色影像合成，这样求解出的 权系数值才能准确反映各参与融合的波段的信息。

三、建立基于空间与光谱相关的回归模型（SSC）；

融合中光谱的保持与多光谱波段结构有关，而空间细节的保持与波段内的空间结构有关。 计算Pan_Syn时传统的SVR只考虑了PAN与多光谱波段的关系，然而我们都知道影像中像元 之间是空间相关的，如果一个像元对融合影像有贡献，这个像元周围的像元也会对融合影像 有贡献；PAN波段覆盖范围一般较大，不仅含有多光谱影像的光谱信息而且含有多于多光谱 影像的空间细节信息。基于以上分析，建立了基于空间与光谱相关的回归模型，模型如公式 (12)。该模型能够更好地描述高空间分辨率的PAN与低分辨率的多光谱之间的关系，得到更 准确的权值，从而得到更加准确的合成PAN值

式中，β是待求的空间相关成分的权值，G_s是空间相关成分（即高频成分），上述公式(12) 的矩阵形式如下式所示

Pan_ori＝BX    (13)

其中为待求的权系数矩阵，矩阵维数为N+1， B＝[MS₁,MS₂,...MS_N,G_s]为已知的光谱和空间相关成分。

为了更好的获得高频成分G_s，我们选用了高斯函数卷积的方法。高斯函数具有低通性质， 并且可以模拟人眼视觉机理，通过不同的尺度参数σ设定，能够模拟人眼远近处所看到的图 像。高斯函数表达式为

$G(x,y;\sigma )={\mathrm{Ke}}^{-\frac{x^2+y^2}{{2\sigma }^2}}---\left(14\right)$

式中，系数K为归一化值，(x,y)是指影像中以滤波中心点为坐标系原点的坐标。将高斯函数 与影像进行卷积，可以得到低频成分的信息，从原始影像中减去低频成分，即可得到高频成 分G_s，表示卷积操作，G_s通过下式决定

$G_s={\mathrm{Pan}}_{\mathrm{ori}}-G(x,y;\sigma )\otimes {\mathrm{Pan}}_{\mathrm{ori}}---\left(15\right)$

为了减少光谱偏差，本发明σ=1，高斯核窗口kernelsize越大，G_s高频细节信息越多， 但是计算量大，反之亦然，为了取得好的融合效果，根据误差3σ规则，本发明建议设高斯 核窗口kernelsize=3σ。

四、建立自适应局部优化处理策略

在地物类型分布严重不均衡情况下，利用影像的整体信息来计算权值，会引起合成的全 色和原始全色影像光谱存在较大不同，从而引起偏色问题。考虑到影像像元的空间差异性， 局部处理可能会比全局处理具有更大的优势。所以将整个影像可分为很多影像块，利用每个 影像块里面的像素分别建立光谱与空间相关回归模型，并采用步骤5中非负最小二乘方法进 行求解，得到每个影像块里像素的光谱与空间相关权值，从而得到合成的全色PAN值。

设h和l分别为PAN和MS的空间分辨率，则r＝l/h为两者空间分辨率的比值，综合考 虑融合效果与处理效率，分块大小w×w，本发明实验中取w＝5r+1，建议w值不大于256， 不小于6。

五、利用非负最小二乘方法求解回归模型中的权系数及异常处理；

实际中，权系数矩阵X是未知的，传统SVR算法利用最小二乘理论进行求解，求解公式 如式(16)，为X的最小二乘估计量，该式不能保证求解出的系数非负，一旦出现负值就 使系数在SVR算法中失去本身的意义。

${\hat{X}}_{\mathrm{LS}}={\left(B^{T}B\right)}^{-1}{B}^T{\mathrm{Pan}}_{\mathrm{ori}}---\left(16\right)$

本发明采用非负约束最小二乘方法进行求解，可以使求解的系数非负，使系数具有物理 意义。非负约束最小二乘是一种较为成熟的优化理论,原理描述如下:

Minimize LSE＝(BX-Pan_ori)^T(BX-Pan_ori)subject to X≥0    (17)

LSE为非负约束的优化误差，X≥0表示对于所有的1≤k≤N+1都有x_k≥0。为了应用拉 格朗日方法解求优化问题，我们引入N+1维未知的被动约束常量c＝[c₁,c₂,...c_N+1]^T＞0,c_k＞0， k的取值为1,2,…N+1。通过c，我们可以形成拉格朗日式J：

$J=\frac{1}{2}{(\mathrm{BX}-{\mathrm{Pan}}_{\mathrm{ori}})}^T(\mathrm{BX}-{\mathrm{Pan}}_{\mathrm{ori}})+\lambda (X-c)---\left(18\right)$

当可以推出

$B^{T}B{\hat{X}}_{\mathrm{NCLS}}-B^T{\mathrm{Pan}}_{\mathrm{ori}}+\lambda =0---\left(19\right)$

从公式(19)可以得出公式(20)

${\hat{X}}_{\mathrm{NLCS}}={\left(B^{T}B\right)}^{-1}{B}^T{\mathrm{Pan}}_{\mathrm{ori}}-{\left(B^{T}B\right)}^{-1}\lambda$

(20)

$={\hat{X}}_{\mathrm{LS}}-{\left(B^{T}B\right)}^{-1}\lambda$

和公式(21)

$\lambda =B^T({\mathrm{Pan}}_{\mathrm{ori}}-B{\hat{X}}_{\mathrm{NCLS}})---\left(21\right)$

公式(20)和(21)可以用来解求非负约束优化解和拉格朗日系数λ＝(λ₁λ₂...λ_N+1)^T。

本发明求解上述式子采用了现有的快速非负约束最小二乘（FNNLS,Fast  non-negativity-constrained least squares）算法，该算法主要思想是将估计量分离为活动集R 和被动集P，前者包括中所有负值和零值，后者包括所有剩余的正值。FNNLS处理开始 时被动集P为空集，假设活动集R包括中所有元素，通过公式(20)和(21)迭代处理调节活 动集和被动集。具体实施时，迭代结束条件可设为活动集R为空集或者活动集R不为空集且 所有拉格朗日系数λ_k,k∈R为负值。

利用FNNLS算法求解权系数矩阵过程中可以看出需对实对称矩阵求逆，这样就要考虑 该矩阵是否可逆的问题，如果该矩阵可逆，则可以利用FNNLS方法较容易的求得准确的权 系数，反之则需进行异常情况的处理。通过实验发现实对称矩阵不可逆是由于用来求解系数 的像素灰度值变化不大，非常接近，如大片水的区域容易发生这种情况。当系数求解不出时， 我们采用公式(22)，即系数等于块内所有多光谱像素的和与所有全色像素的和的比值。

其中，np为用于求解系数的像素个数，j的取值范围为1,2,…np，Pan_j是全色影像第j个 像素的灰度值，MS_i，j是多光谱影像第i波段第j个像素的灰度值。一般用于求解系数的像素 采用影像块内所有像素，np=w×w。

六、利用双线性内插解决局部优化处理引起的分块效应问题，对每个像素逐一进行融合，得 到融合影像。

局部分块优化处理往往会产生分块效应，综合考虑融合效果和时间复杂度的影响，我们 采用双线性内插利用影像块获得的权系数，来获得每个像素的权系数将带入公式(10)得 到合成的低分辨率全色影像值Pan_Syn，最后利用多光谱关系模型(9)获得每个点的融合值，从 而得到整个融合影像。

为了检验本发明的优越性，选用了Spot-5和QuickBird影像进行实验，此外，选择 Gram-Schmidt,HCS,HPF,PCA以及UNB-Pansharp算法作为对比。

融合影像质量评价：对于融合影像的质量评价包括主观和客观指标评价两个方面，好的 融合结果不仅能保持原始多光谱影像的光谱信息，而且能够同时较好地保持全色波段的空间 信息。为了更好的客观评价融合结果，本发明采用Wald的策略，即将原始多光谱和全色影像 进行降维，然后对降维数据进行融合，原始多光谱影像作为参考数据，融合后多光谱影像与 参考影像进行比较。本发明采用了相关系数(CC)、相对整体误差（relative global dimensional  synthesis error，ERGAS）、通用质量指数(universal image quality indicator，UIQI)和光谱扭曲度 (SF)四种客观评价指标。

1)相关系数和通用质量指数：反应了融合后影像与参考影像的相似度，这两项指标越 大表明接近程度越大，即空间和光谱保持越好；

2)相对整体误差和光谱扭曲度：反应了空间和光谱的差别和扭曲程度，这两项指标越 小，表明空间细节与光谱变形越小，融合结果越好。

实验一：Spot-5卫星降维后的全色影像（空间分辨率10m），降维后的多光谱影像（空 间分辨率20m）和原始多光谱影像（空间分辨率10m），多光谱采用波段4,3,和2假彩 色合成显示。

实验二：QuickBird卫星降维后的全色影像（空间分辨率2.8m），降维后的多光谱影像 （空间分辨率11.2m）和原始多光谱影像（空间分辨率2.8m），多光谱采用波段3,2,和1 真彩色合成显示。

主观评价：

主观评价即目视视觉判断，从实验一中可以看出所有融合方法都能够较好地提高原始多 光谱影像空间分辨率，但是与原始参考多光谱影像相比，Gram-Schmidt,HCS,PCA和 UNB-Pansharp融合结果有不同程度的偏色现象，见红色矩形框区域，HPF融合结果与 本发明方法相对空间细节较为模糊；从实验二中可以看出除了本发明融合结果外，其它 方法融合结果都有不同程度偏色现象。从目视评价看来，相比以上方法，本发明的方法 在保持空间细节的同时可以很好地保护光谱信息。

客观评价：

两个实验的客观评价结果见表1和表2，评价最好的结果用黑体标注。

表1各方法融合结果的定量评价（Spot-5数据）

表2各方法融合结果的定量评价（QuickBird数据）

由实验一(表1)的实验数据可知，本发明方法具有除MIR波段之外的所有CC最大值， SF最小值，ERGAS最小值以及UIQI最大值。虽然HPF方法的MIR波段的CC值比本发明 方法大，但是在别的指标上均低于本发明方法。综合分析表1各指标数据，可以看出相比其 他方法本发明方法能够更好的保持原始影像的光谱和空间特征。

由实验二（表2）的实验数据可知，相比其他方法，PCA方法在NIR波段的各项指标都 是最差的，有的指标甚至是负值，因此，即使在R和G波段的CC和UIQI指标具有较高的 值，该方法仍然严重扭曲了原始影像的各项信息。本发明融合结果具有最优的SF和ERGAS 指标，因此，该方法能够最大程度的减少光谱与空间信息的扭曲；此外，本发明融合方法与 原始影像的CC和UIQI大于PCA以外的其它方法，表明该方法很好地保持了原始影像的光 谱与空间信息。

从实验一和实验二可知，本发明融合结果能够同时具有高的空间分辨率和丰富的光谱信 息，能够同时保持原始影像的空间与光谱特征信息。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技 术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不 会偏离本发明的精神或者超越所附权利要求书所定义的范围。