CN104700379A - 一种基于多尺度形态成分分析的遥感图像融合方法 - Google Patents

Abstract

本发明公开了一种基于多尺度形态成分分析的遥感图像融合方法，属于信号处理与遥感图像处理交叉领域。在不同尺度下，分别对高分辨率遥感图像和多光谱遥感图像进行形态成分分析，采用迭代收缩法进行稀疏分解，把待融合的目标图像分成多个尺度的纹理成分和卡通成分，去除高分辨率图像中的卡通成分和噪声成分、多光谱图像中的纹理成分和噪声成分，保留高分辨率图像的有效尺度纹理成分、多光谱图像中的卡通成分，并对二者进行稀疏重建，获得融合图像。基于多尺度形态成分分析的遥感图像融合方法能对高分辨率遥感图像和多光谱遥感图像进行有效地融合，相比现有的融合方法具有更高的空间分辨率和更低的光谱失真，并且比现有的稀疏重建法有速率有很大的提升。

Description

一种基于多尺度形态成分分析的遥感图像融合方法

技术领域：

本发明属于信号分析、遥感图像处理交叉领域，是一种基于多尺度形态成分分析的遥感图像融合方法。

背景技术：

随着计算机技术、航空航天技术以及遥感技术的飞速发展，图像资源的获取方式(多传感器、多平台)日益丰富，获得的遥感图像也呈现多空间分辨率、多光谱和多时相等特点。目前，遥感图像存在的问题有：①单一遥感图像往往只包含一部分地物特征，而不能反映出全部内容；②与单源数据相比，多源数据之间是互补的，同时存在一定量的冗余；③由于多源遥感数据信息量大、冗余度高，其利用率低。为了能充分利用多源遥感数据，发挥各传感器优势的数据，将不同类型遥感图像数据进行融合是一条行之有效的途径。遥感图像的融合可为决策提取更丰富、更有用、更可靠的信息。

遥感图像融合目的就是使融合的结果包含最大的信息量，在保证光谱失真率尽可能小的同时，能最大限度提高空间分辨率。目前已有很多的遥感图像融合方法，如IHS变换法、Brovey法、主成分分析法、小波变换法、稀疏重建法以及这些方法的组合，这些方法从不同的层面对多源遥感数据去除了部分冗余，取得了一定的融合效果。然而，现有的这些方法仍存在不同程度的弊端：

1)IHS变换法、Brovey法和主成分分析法的融合结果存在较大的光谱失真，且空间分辨率有较大的提升空间；(参考对比文件1、对比文件2、对比文件3)

2)小波变换分析法会引人较多的人为噪声，且存在较大的光谱失真；(参考对比文件4)

3)基于基追踪和匹配追踪的稀疏重建法计算量巨大，仿真实验常采用小图像进行测试，对真实的大幅遥感图像难以投入实际应用。(参考对比文件5)

形态成分分析(MCA,Morphological Component Analysis)是一种新型的图像分析方法，它能集合多种基的优点来更加稀疏的描述和分解图像，取得了比离散余弦变换、小波变换等变换域分析法和基追踪、匹配追踪等稀疏重建法更好的效果(参考对比文件6)，并且形态成分分析的稀疏基寻优实现是基于迭代收缩算法，有较好的执行效率。然而，现有形态成分分析法是在单一尺度下进行的，这对地物成分复杂的遥感图像分析和分解不利。因此，本发明首次提出在多个尺度下利用形态成分分解遥感图像，选择高分辨率图像的最有效尺度纹理成分与多光谱图像的卡通成分(分段平滑成分)，然后采用稀疏重建的方法进行融合。对比文件1：M.Choi.A new intensity-hue-saturation fusion approach to image fusion with a tradeoff parameter.IEEE Transactions on Geoscience and Remote Sensing,2006,44(6):1672–1682.

对比文件2：周会珍，吴双，毛德发等.多卫星传感器数据的Brovey融合改进方法.遥感学报，2012,16(2):343–360

对比文件3：张庆河,邹峥嵘,余加勇.遥感图像像素级融合方法比较研究.测绘工程,2008,17(4):35–42.对比文件4：G.Pajares,J.Cruz.A wavelet-based image fusion tutorial.Pattern Recognition,2004,37(9):1855–1872.

对比文件5：尹雯,李元祥,周则明,刘世前.基于稀疏表示的遥感图像融合方法.光学学报,2013,33(4):31–38.

对比文件6：李映,张艳宁,许星.基于信号稀疏表示的形态成分分析:进展和展望.电子学报,2009,37(1):146–152

发明内容：

本发明提出了一种基于多尺度形态成分分析的遥感图像融合方法，联合曲波变换基Φ₁和局部离散余弦变换基Φ₂作为稀疏分解字典Φ＝[Φ₁,Φ₂]，把遥感图像按字典系数阈值从不同的尺度分解成纹理成分和卡通成分，选取高分辨率遥感图像的前4个尺度的纹理成分和多光谱遥感图像的卡通(分段平滑)成分进行稀疏重建，重建结果即为所得融合图像。具体步骤为：

1)对高分辨率遥感图像(全色，SAR等)进行多尺度形态成分分解，保留纹理细节成分去掉

卡通成分和最精细的纹理成分(噪声)；

2)对多光谱遥感图像进行单一尺度形态成分分解，保留卡通成分去掉纹理成分(含噪声)；

3)合成高分辨图像纹理成分和多光谱图像卡通成分，进行重建，即重建图像

$I_{\mathrm{Rec}}=({\mathrm{\α}}_{\mathrm{MS}}^A+{\mathrm{\α}}_{\mathrm{HR}}^B)\mathrm{\Φ}---\left(1\right)$

4)重建图像I_Rec即为融合结果。

附图说明：

图1TerraSAR-X高分辨率遥感图像多尺度分解得到的8个尺度纹理成分。

图2各种融合方法下的TerraSAR-X高分辨率遥感图像与SPOT5多光谱图像的融合结果。

图3各种融合方法下的Landsat ETM+高分辨率遥感图像与中巴资源卫星多光谱图像的融合结果。

具体实施方式：

曲波变换基能有效的捕捉图像卡通成分，局部离散余弦变换基能有效的捕捉图像的纹理成分。所以，联合曲波变换基Φ₁和局部离散余弦变换基Φ₂作为形态成分分析的分解字典Φ＝[Φ₁,Φ₂]，分别对高分辨率遥感图像和多光谱遥感图像进行多尺度形态成分分解。舍弃高分辨率遥感图像的含噪声内容的尺度部分，保留其他尺度分解下的纹理成分；舍弃多光谱图像的纹理成分，保留TV约束下多光谱图像的卡通(分段平滑)成分。把保留的高分辨率遥感图像的纹理成分和多光谱图像的分段平滑成分作为新图像的纹理成分和卡通成分进行重建，获得融合结果。

其中高分辨率遥感图像和多光谱图像分解为纹理成分和卡通成分的过程如下：

设第i尺度下的分解字典其中，是第i尺度下的曲波基字典，是第i尺度下的局部离散余弦变换基字典。则在第i尺度下的高分辨率图像I_HR和多光谱图像I_MS可分解表示为

$I_{\mathrm{HR}}=\underset{i=1}{\overset{T}{\mathrm{\Σ}}}(I_{\mathrm{HR}}^{\mathrm{Ai}}+I_{\mathrm{HR}}^{\mathrm{Bi}})=\underset{i=1}{\overset{T}{\mathrm{\Σ}}}({\mathrm{\α}}_{\mathrm{HR}}^{\mathrm{Ai}}{\mathrm{\Φ}}_1^i+{\mathrm{\α}}_{\mathrm{HR}}^{\mathrm{Bi}}{\mathrm{\Φ}}_2^i)---\left(2\right)$

$I_{\mathrm{MS}}=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}(I_{\mathrm{MS}}^{\mathrm{Ai}}+I_{\mathrm{MS}}^{\mathrm{Bi}})=\underset{i=1}{\overset{K}{\mathrm{\Σ}}}({\mathrm{\α}}_{\mathrm{MS}}^{\mathrm{Ai}}{\mathrm{\Φ}}_1^i+{\mathrm{\α}}_{\mathrm{MS}}^{\mathrm{Bi}}{\mathrm{\Φ}}_2^i)---\left(3\right)$

其中，和表示第i尺度下高分辨率图像I_HR的卡通成分和纹理成分，和表示第i尺度下高分辨率图像I_HR的卡通成分和纹理成分的稀疏分解系数，和表示第i尺度下多光谱图像I_MS的卡通成分和纹理成分，和表示第i尺度下高分辨率图像I_HR的卡通成分和纹理成分的稀疏分解系数，T和K表示分解的尺度数。

卡通成分采用了全变差(TV,Total Variation)惩罚处理，使这部分更加适合分段平滑模型，提高重构质量。在字典Φⁱ下对输入图像进行稀疏分解，即求解：

$\underset{I_1,I_2}{\min }\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{\left|\right|{\left({\mathrm{\Φ}}_k^i\right)}^T{I}_k\left|\right|}_1+\mathrm{\γ}{\left|\right|I_A\left|\right|}_{\mathrm{TV}},$ 满足 ${\left|\right|I-\underset{k=1}{\overset{2}{\mathrm{\Σ}}}{I}_k\left|\right|}_2\≤\mathrm{\σ}---\left(4\right)$ 其中，表示对I_k做相应的正向变换(如曲波变换或局部离散余弦变换)，得到变换系数α。I_A表示卡通(分段平滑)部分，γ是TV正则化参数，卡通部分的惩罚项TV为：

${\left|\right|I_A\left|\right|}_{\mathrm{TV}}=\underset{i=1}{\overset{N_1}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_2}{\mathrm{\Σ}}}\sqrt{{[I_A(i+1,j)-I_A(i,j)]}^2+{[I_A(i,j+1)-I_A(i,j)]}^2}---\left(5\right)$

其中，I_A(i,j)表示卡通部分I_A在位置(i,j)处的值。

通过曲波变换基和局部离散余弦变换基的系数大小来控制分解尺度，在某尺度下对输入图像进行稀疏分解(为了简化，分解字典不再带上标i表示尺度)，求解式(4)的迭代收缩算法步骤如下：

Step1：输入图像I，分解字典Φ＝[Φ₁,Φ₂]，迭代次数N_iter，γ，停止阈值λ_min。

Step2：初始化： 

①初始解： $I_k^{\left(0\right)}=0,\∀k.$

②初始残差：r⁽⁰⁾＝I

③初始阈值： $k^{*}={\max }_k{\left|\right|{\mathrm{\Φ}}_k^{T}I\left|\right|}_{\∞},$ 设置 ${\mathrm{\λ}}^{\left(0\right)}={\max }_{k\≠k^{*}}{\left|\right|{\mathrm{\Φ}}_k^{I}I\left|\right|}_{\∞}$

Step3：主迭代过程：

从t＝1到N_iter

从k＝1到2

①计算边缘残差： $r_k^{\left(t\right)}=r^{(t-1)}+I_k$

②根据阈值更新第k个成分系数： ${\mathrm{\α}}_k^{\left(t\right)}={\mathrm{TH}}_{{\mathrm{\λ}}^{(t-1)}}\left({\mathrm{\Φ}}_k^T{r}_k^{\left(t\right)}\right)=\max (1-\frac{{\mathrm{\λ}}^{(t-1)}}{{\mathrm{\Φ}}_k^T{r}_k^{\left(t\right)}},0)$

③更新第k个成分：

④如果k是卡通部分且γ≠0，应用TV对进行约束

更新残差： $r_k^{(t+1)}=I-{\mathrm{\Σ}}_{k=1}^2{I}_k$

更新阈值：λ^(t)＝λ^(t-1)-λ_min/N_iter

如果λ^(t)≤λ_min，停止

停止迭代

Step4：输出：形态成分和稀疏系数α^A、α^B，即对应的卡通成分和纹理成分。

重建图像时，分别选取高分辨率图像和多光谱图像形态成分对应的纹理成分和卡通成分按照式(1)进行重建，重建结果即为所得的融合图像。

仿真结果： 

设置形态成分分析的分解算法参数为：N_iter＝300，λ_min＝10^-6，γ＝1。

1)高分辨率遥感图像的多尺度形态成分分析的纹理成分

高分辨率遥感图像选自德国TerraSAR-X卫星于2008年拍摄的珠三角地区的空间分辨率为1m的SAR图像(1024×1024pixels)，实验对8个尺度下的纹理成分进行了提取，第1尺度到第8尺度的系数阈值参数分别为：0.0625、0.125、0.1875、0.25、0.375、0.5、0.625、0.75。实验结果如图1所示，可见，第5尺度及第5尺度以下的纹理成分微弱，大部分为噪声。实验中的第5尺度的字典中的系数阈值>0.25。

2)高分辨遥感图像与多光谱遥感图像的融合

融合结果的客观指标计算方法和度量标准如下：

①峰值信噪比(PSNR,Peak of Signal to Noise Ratio)。表示融合结果中的信号和噪声的功率比，计算方法如式(6)所示。

$\mathrm{PSNR}=10\×\log \left(\frac{{255}^2}{\mathrm{MSE}}\right)---\left(6\right)$

其中，I_F(i,j)与I_R(i,j)分别是融合图像与参考图像在(i,j)处的像素值，Framesize表示图像的大小。PSNR越大，说明噪声含量越小。

②相关系数(CC,Correlation Coefficients)。反映融合结果与参考图像之间的相关程度，计算方法如(7)式所示。

$\mathrm{CC}(I_F,I_R)=\left|\frac{{\mathrm{\Σ}}_{i=1}^M{\mathrm{\Σ}}_{j=1}^N(I_F(i,j)-\stackrel{\‾}{I_F})(I_R(i,j)-\stackrel{\‾}{I_R})}{\sqrt{{\mathrm{\Σ}}_{i=1}^M{\mathrm{\Σ}}_{j=1}^N{(I_F(i,j)-\stackrel{\‾}{I_F})}^2{(I_R(i,j)-\stackrel{\‾}{I_R})}^2}}\right|---\left(7\right)$

其中，与分别是融合图像与参考图像的像素平均灰度值。CC值越大，说明两幅图像间相关程度越高，即越相似。

③光谱角(SAM,Spectral Angle Mapper)。表示融合图像与参考图像之间的光谱扭曲程度，计算方法如(8)式所示。

$\mathrm{SAM}(v,\hat{v})=\arccos \left(\frac{\⟨v,\hat{v}\⟩}{{\left|\right|v\left|\right|}_2{\left|\right|\hat{v}\left|\right|}_2}\right)---\left(8\right)$

其中，向量v和都是n维的，v＝{v₁,v₂,…,v_n}表示融合图像n个波段像素点的集合，

表示对应的参考图像n个波段像素点的集合。若光谱角值等于0，则表示融合后的光谱没有扭曲误差。④相对整体维数综合误差(ERGAS,Relative Dimensionless Global Error in Synthesis)。表示融合图像与参

考图像之间的光谱扭曲程度，计算方法如(9)式所示。

$\mathrm{ERGAS}=100\frac{h}{l}\sqrt{\frac{1}{k}\underset{b=1}{\overset{K}{\mathrm{\Σ}}}{\left(\frac{\mathrm{RMSE}\left(b\right)}{\mathrm{\μ}\left(b\right)}\right)}^2}---\left(9\right)$

其中，h/l表示高分辨率图像和多光谱图像像素值的大小比，K表示波段个数，μ(b)表示第b个波段的均值；RMSE(b)表示第b个波段图像与参考图像间的均方根误差。融合光谱质量越高，ERGAS越小，理想值为0。

⑤信息熵(IE,Information Entropy)。对融合图像总信息量的描述的指标为信息熵，计算方法如式(10)所示。

$\mathrm{IE}=-\underset{i=0}{\overset{L-1}{\mathrm{\Σ}}}{p}_i{\log }_2^{p_i}---\left(10\right)$

其中，p＝{p₀,p₁,…,p_i…,p_L-1}表示一幅图像的灰度分布，p_i为像素值为i的像素数与图像总像素数之比。一般来说，IE越大越好，有的融合方法降噪效果比较好，也可能会使融合图像的信息熵比较小，因此，该指标可以辅助做遥感图像融合方法性能评价。

仿真1：TerraSAR-X高分辨率遥感图像与SPOT5多光谱遥感图像之间的融合

本组实验数据的高分辨率遥感图像依然选用德国TerraSAR-X卫星于2008年获取的珠三角地区的空间分辨率为1m的雷达图像，如图2(b)所示(1024×1024pixels)，多光谱遥感图像选自SPOT5号卫星拍摄的相应地区的分辨率为10米的多光谱图像，是由B3、B2、B1三个波段合成，如图2(a)所示(1024×1024pixels)。

根据1)中的实验结果，选取高分辨率遥感图像的第1尺度至第4尺度的纹理成分与多光谱的卡通成分进行融合。采用IHS、Brovey、主成分分析法、小波变换法、稀疏重建法和多尺度形态成分分析法融合的结果分别如图2(c–h)所示。表1给出了不同融合方法融合结果的统计参数对比，参考图像为待融合的多光谱遥感图像。

表1仿真1中的不同融合方法融合结果的评价参数对比

仿真2：Landsat ETM+高分辨率全色图像和中巴资源卫星(CRERS)多光谱图像之间的融合

本组实验采用的数据是Landsat ETM+高分辨率全色图像(1024×1024pixels，分辨率19.5米)和中巴资源卫星(CRERS)多光谱图像(1024×1024pixels，分辨率15米)，位于广东省珠海市斗门区，实验区地 物涵盖了水体、林地、建筑用地、农业用地和未利用地等多种类型。其中，多光谱图像的显示是由B3、B2、B1三个波段组成的假彩色图像，如图3(a)所示，Landsat ETM+高分辨率全色图像如图3(b)所示。

根据1)中的实验结果，尺度选择与仿真1中相同，依然选取高分辨率遥感图像的第1至第4尺度的纹理成分与多光谱的卡通成分进行融合。采用IHS、Brovey、主成分分析法、小波变换法、稀疏重建法和多尺度形态成分分析法融合的结果分别如图3(c–h)所示。表2给出了不同融合方法融合结果的统计参数对比，参考图像为待融合的多光谱遥感图像。

表2仿真2中的不同融合方法融合结果的评价参数对比

由仿真1和仿真2的实验结果，从目视效果(图2和图3)和客观指标统计(表1和表2)两个方面分别对原始遥感图像和融合后图像进行对比和分析，可以得出如下结论：

(1)目视来看，相比IHS、Brovey、主成分分析法、小波变换法，基于多尺度形态成分分析的遥感图像融合方法得到融合图像的颜色(光谱信息)与原始多光谱图像最相似，这说明多尺度形态成分分析法的光谱失真最小；相比IHS、主成分分析法、稀疏重建法，基于多尺度形态成分分析法有更高的空间分辨率；

(2)由于目前还没有一种有效的遥感图像融合的客观度量指标，通常做法是采用多指标综合评测，由表1和表2可见，基于多尺度形态成分法的大部分评价参数占优；算法执行速率(CPU Time)比稀疏重建法有很大的提升。

Claims (7)

1.基于多尺度形态成分分析的遥感图像融合方法，其特征在于：联合曲波变换基和局部离散余弦变换基作为分解字典，对高分辨率的遥感图像按系数阈值进行多个尺度形态成分分解，得到纹理成分和卡通成分，选取有效尺度的纹理成分；对多光谱遥感图像进行单一尺度的形态成分分解，得到纹理成分和卡通成分，选取卡通成分；对选取的高分辨率遥感图像的有效纹理成分和多光谱图像的卡通进行稀疏重建，得到融合结果。

2.如权利要求1所述的基于多尺度形态成分分析的遥感图像融合方法，其中形态成分分析特征在于：采用迭代收缩算法把图像分解成纹理成分和卡通成分。

3.如权利要求1所述的基于多尺度形态成分分析的遥感图像融合方法，其中纹理成分特征在于：局部离散余弦变换基捕捉的成分。

4.如权利要求1所述的基于多尺度形态成分分析的遥感图像融合方法，其中卡通成分特征在于：图像中的分段平滑成分，即曲波变换基捕捉的成分。

5.如权利要求1所述的基于多尺度形态成分分析的遥感图像融合方法，其中选取有效尺度的纹理成分特征在于：选取第1到第4尺度的纹理成分。

6.如权利要求1所述的基于多尺度形态成分分析的遥感图像融合方法，其中高分辨率的遥感图像特征在于：高空间分辨率、单一波段的遥感图像。

7.如权利要求1所述的基于多尺度形态成分分析的遥感图像融合方法，其中多光谱遥感图像特征在于：同一传感器、多个波段的遥感图像。