CN102521815A - 影像快速融合系统及快速融合方法 - Google Patents

Abstract

本发明提供了一种遥感影像快速融合系统及实现方法，融合系统包括：图像预处理模块；图像IHS并行融合模块；图像PCA并行融合模块；图像小波并行融合模块；图像融合效果并行评估模块，所述图像预处理模块和图像融合效果评估模块均与所述图像IHS并行融合模块、图像PCA并行融合模块和图像小波并行融合模块相连。融合方法采用并行融合和评价。本发明提供了一种简单方便又高效的遥感影像快速融合实现方法，为灾害预防及救灾、军事等领域的遥感应用提供了高效的遥感影像快速融合方法与实现系统，解决了大数据量遥感影像融合速度慢的问题。

Description

影像快速融合系统及快速融合方法

技术领域

本发明涉及信息处理技术领域，特别涉及图像并行融合方法，具体为遥感影像快速融合系统及实现方法。

背景技术

多传感器图像融合属于多传感器信息融合的范畴，是指按照一定的算法，将不同传感器获得的同一景物的图像或同一传感器在不同时刻获得的同一景物的图像合成到一幅满足给定要求的图像中。单一传感器由于受由光的能量和衍射决定的分辨极限、成像系统的调制传递函数、信噪比三个方面的限制，要同时获得光谱、空间和时间的高分辨率是很难的。多传感器图像融合技术由于可以有效的利用多幅图像提供的互补信息和冗余信息，因此融合后的图像对场景的描述比任何单一源图像都更全面、精确。一般而言，使用多传感器信息融合技术具有以下优点：(1)可提高系统的可靠性和鲁棒性；(2)可扩展空间和时间上的观测范围；(3)可提高信息的精确程度和可信度；(4)可提高对目标物的监测和识别性能；(5)可降低对系统的冗余投资。

遥感影像数据融合分为三级：像元级融合、特征级融合、决策级融合。像元级融合是一种低水平的融合，它是直接在采集到的原始数据层上进行的融合，通常用于多源图像复合、图像分析和理解等。特征级融合是先对原始遥感影像信息进行特征提取，然后对特征进行综合分析和处理，融合结果能最大限度地给出决策分析所需要的特征信息。决策级融合是一种高层次的融合，它是在上述像元级和特征级融合所提供的各类特征信息的基础上对图像信息进行识别、分类或目标检测，并在获得有关区域决策信息后，再对所获得的专题图像进行融合处理，它的融合结果直接为指挥、控制、决策系统提供依据。

图像融合是多传感器数据融合的一个重要分支，而多传感图像融合技术最早是被应用于遥感图像的分析和处理中。1979年，Dail等人首先把雷达图像和Landsat-MSS图像的复合图像应用于地质解译，其处理过程可以看作是简单的图像融合。1981年，Laner和Todd进行了Landsat-RBV和MSS图像数据的融合试验。到80年代中后期，图像融合技术开始引起人们的关注，陆续有人将图像融合技术应用遥感多光谱图像的分析和处理，90年代以后，随着多颗遥感卫星JERS-1、ERS-1、Radarsat等的发射升空，遥感影像融合技术成为遥感图像处理和分析的研究热点之一。在融合方法上，有一些经典的算法，如各种比值和加权乘法、Brovey变换法、高通滤波法、IHS变换法(也叫HIS变换法)、小波变换法、PCA(主成分分析)法。近年来，小波变换融合方法已用于多传感器影像数据的融合。Jorge Nunez等基于多分辨率小波分解系数加到低分辨率全色图像的小波分解系数中去，通过使用这种方法对SPOT和LANDSAT(TM)图像进行融合实验，结果表明该技术在保持光谱和空间信息明显优于IHS和LHS方法融合结果。IHS变换是遥感影像融合中较为常用的方法，但有其缺陷。为此，许多研究者对IHS变换融合提出了多种改进算法，哈斯巴干、贾永红、王智均等结合小波变换对I分量进行了改进，最大限度地保留了原多光谱影像的光谱信息。李军等将小波的多分辨率分析与IHS变换相结合，提出了叠加融合的新方法，给出了黑白航空影像与TM影像，SAR影像与TM影像的融合结果。它先对高分辨率影像进行小波分解，然后将得到的各小波面叠加到多光谱影像经IHS变换后的强度I影像中。根据遥感影像融合相关文献，小波变换用于遥感影像融合算法主要有二进小波变换、多进制小波变换、多分辨率小波变换、小波包变换等。影像融合模式有小波变换对IHS融合法的改进，和基于小波变换的加权融合等。王相海等对遥感图像经Contourlet变换后的高频子带系数分布的方向特征进行统计分析，发现遥感图像经Contourlet变换后高频系数的分布具有较强的方向区域特征，在此基础上，提出一种基于Contourlet系数方向区域相关性的遥感图像融合算法。实验结果表明，该算法在提高融合图像空间分辨率的同时能够更好地保留原始多光谱图像的光谱信息。该文针对超高分辨率的全色光图像和多光谱图像的融合，所提出的融合算法在对多光谱数据进行对应分析的基础上，利用冗余小波变换提取出全色光图像的空间细节信息并将其融入到成分空间，相比现有同类融合方法，该方法能够在提高空间分辨率的同时更好地保持光谱特性，有效地减少了色彩失真的现象。柴勇等提出了B-Spline塔型方向滤波器组，由B-Spline金字塔构造误差金字塔，得到与Laplacian金字塔类似的误差金字塔分解结构，然后再采用方向滤波器组进行滤波，得到B-Spline塔型方向滤波器组。仿真表明，新算法能得到比Contourlet变换更好的融合效果，信息损失少。胡根生等研究了一种基于支持向量值轮廓波变换的遥感影像融合方法，首先将支持向量值滤波器和方向滤波器组进行合理搭配，构造支持向量值轮廓波变换，该变换具有平移不变、泛化能力好、捕捉奇异性能强等特性。然后利用该变换对多源遥感影像进行多尺度、多方向、多分辨率分解，在不同的分解水平上利用基于区域能量和基于区域的轮廓波对比度方法进行融合。实验结果表明利用该变换进行遥感影像融合时，能在不降低空间分辨率的情况下有效保留源影像的光谱信息。

图像数据融合技术作为一门很有优势的技术，成为国际技术研究和应用的热点。纵观国内外当前此领域的研究，多源遥感数据融合的面临一些主要问题为：(1)缺乏统一的数据融合模型，特别是数学模型。发展的方向和难点：建立一种相对统一的融合数学模型，以简化现有算法的繁杂性。(2)缺乏对数据融合结果的有效评价手段。现有的评价标准都过于简单，缺乏灵活性，无法适应技术发展的需要。(3)传统的多源数据关联和融合算法，不仅需要较多的先验知识，而且当融合源增加以及传感器或被观测地物目标数目增多时，在计算上出现N-P完全复杂性问题，表现出难以克服的计算组合爆炸现象。(4)由于平台(如卫星、飞机)的机动和不稳定性，多源信息的时空配准和误差补偿已成为影像遥感信息融合性能的主要问题之一。(5)大数据量遥感影像融合速度慢，难以适应一些遥感领域应用需要。多源遥感数据融合的发展方向为：(1)建立多源遥感数据融合的基本理论，对兼有稳健性和准确性的融合算法研究，并深入探讨更为复杂的融合算法。(2)另一方面是将已有的融合方法工程化和商品化，开发能够提供多种复杂融合算法的处理硬件，以便在数据获取的同时就实现实时的融合。(3)建立信息融合系统的设计和评估方法，尤其是针对不同的应用领域和融合算法，构建相应的定性和定量的评价指标。重点要侧重于单因素的评价体系和综合评价体系兼顾并重。(4)遥感影像与GIS数据库基于特征级的融合方法研究。对遥感图像与GIS数据库的大量背景数据进行叠加、分析，可大大提高GIS中的模式识别能力和可信度。(5)遥感数据与专家系统的决策级结合研究。尽管建立一种相对统一的融合数学模型是多源遥感影像融合的发展目标，但是针对不同应用领域或对象而采用不同的融合算法和模型的方法的现状将仍旧在一定的时期存在。一定程度上仍需要紧密结合实际应用背景开展应用研究。

多源遥感数据融合有着非常重要的实际应用意义。一方面，多源遥感图像所提供的信息具有冗余性、互补性和合作性，多源遥感图像数据融合不仅扩大了各数据的应用范围，而且提高了分析精度，应用效果和实用价值。另一方面，在现代遥感技术的快速发展过程中，各种对地观测卫星所提供的各种遥感信息数据量呈海量程度增长，越来越多的遥感图像数据的出现使人们对数据的处理和分析面临更多的困难和挑战。遥感图像数据的处理要远远滞后于遥感图像数据的获取，如何实现对海量数据的实时处理，如何充分利用如此多的遥感图像数据，以满足人们对高质量图像的迫切需求以及对地物观测数据智能化处理的需求，应用多源遥感图像数据融合以综合利用多源遥感图像进行信息提取和分析是人们对遥感数据的使用的迫切要求。而有效的遥感影像数据融合算法都较为复杂，要充分发挥各种对地观测卫星的作用，结合遥感影像数据融合技术的优势，研究快速的遥感影像数据融合方法显得日益迫切。

由目前出现的遥感影像数据融合方法及在遥感影像数据融合的实际研究可知，遥感影像数据融合较为复杂，计算量大，相应的耗时也较长，对计算机的速度和内存等资源都提出了很高的要求。例如，目前串行PCA融合计算过程为：将两幅配准的全色影像和多光谱影像调整为相同图像尺寸，将多光谱影像的R、G、B三个波段数据以一个矩阵输入，设多光谱影像的尺寸为高H，宽W，大小为Size＝H*W，则输入矩阵IM的大小为3*Size，

其中协方差矩阵计算如下，协方差矩阵可表示为

$\left[\begin{array}{ccc}{\mathrm{\σ}}_{11}& {\mathrm{\σ}}_{12}& {\mathrm{\σ}}_{13}\\ {\mathrm{\σ}}_{21}& {\mathrm{\σ}}_{22}& {\mathrm{\σ}}_{23}\\ {\mathrm{\σ}}_{31}& {\mathrm{\σ}}_{32}& {\mathrm{\σ}}_{33}\end{array}\right]$

${\mathrm{\σ}}_{\mathrm{xy}}=\mathrm{COV}(X,Y)=\frac{\underset{i=1}{\overset{\mathrm{Size}}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})}{\mathrm{Size}-1}(x=\mathrm{1,2,3};y=\mathrm{1,2,3})$

X_i，Y_i分别代表第i个像素R、G、B波段中的其中的一个值。通过对串行PCA融合的分析获知，其中的协方差矩阵计算量很大。

虽然计算机发展很快，但计算机的速度和内存等的扩展是有限的，计算机性能的提高需要花较大的代价，因此用其它办法研究遥感影像数据的快速融合对各种对地观测卫星的应用有着非常重要。

在高性能计算中，并行集群计算系统具有较高的性价比和良好的可扩展性，可以满足不同规模的大型计算问题，因此受到日益重视。目前HPC(HighPerformance Computing)的发展方向是集群系统。集群系统是利用高速通用网络将一组工作站或PC机，按某种结构连接起来，在并行程序设计以及可视化人机交互集成开发环境支持下，统一调度，协调处理，实现高效并行处理的系统。集群系统相对于MPP有性能/价格比高的优势，具有良好的可扩展性和容错能力，易于维护和管理。目前世界上的超级计算机500强很多基于这种结构。

由于并行计算的优势，国内外正在研究将并行计算用于图像处理之中。国外对图像处理的并行算法的研究较多，如美国北卡罗莱纳大学、田纳西大学、卡内基-梅隆大学、美国能源部等都做了大量研究，也提供了许多图像并行处理算法，西澳大利亚大学的Thomas

教授等人还出版了有关并行图象处理的专著，比较详细地讨论一般图像处理算法的并行设计方法。在国内，如陈国良教授等著的《并行算法的设计与分析》、《并行算法实践》等著作里面的许多算法思想已被应用于并行图像处理之中。黄国满和郭建峰讨论了分布式并行遥感图像处理中的数据划分问题，通过对不同数据划分方式的净通讯量和实际通讯量的分析，得出了在三类图像处理中数据划分的最佳方式，他们还进一步分析了非平均分配的数据划分问题。李军和李德仁分析和探讨了分布式遥感图像处理中的主要研究问题和若干关键技术，包括分布式处理模型、支撑环境、网络环境下的遥感数据存储管理、分布式并行处理中数据划分、功能分解以及负载平衡与任务分配等。池天河等基于并行计算研究了洪水灾害快速评估系统，并用研究成果对鄱阳湖地区1998年特大洪水遥感数据进行了实验，取得了令人满意的效果。吕捷等基于MPI，论述了图像处理算法的并行实现方法；牟胜梅和刘光明分析了遥感图像预处理算法的可并行性，提取了算法要素，提出了适用于大数据量遥感图像实时预处理的分布并行处理系统的结构设计方案。中科院计算所王贞松领导的小组实现了星载SAR数据成像的并行和实时处理。周非等针对ScanSAR成像机制的特殊性，提出了一种基于多节点-多进程的异步并行模式的星载ScanSAR并行处理的技术。实验结果表明，相对于传统的并行处理方法，该方法的成像时间得到了较好的改善，具有一定的优越性。

总之，当前国内外关于并行程序和算法设计、图像并行处理的研究正在大力进行，但在并行遥感影像融合方面的研究成果还也较少。研究有效的并行遥感影像融合方法，开发有效的并行遥感影像融合系统并探讨其应用显得日益重要。

为此，为了更好地发挥地遥感的作用，满足遥感在灾害预防及救灾、军事等邻域的实际应用需要，需要找到遥感影像融合的高效方法，开发有效的并行遥感影像融合系统，以实现大数据量遥感影像快速融合。

发明内容

为了克服上述所述的遥感融合现有技术中的不足，本发明提供了一种影像快速融合系统及实现方法。

本发明的影像快速融合系统，该系统基于并行集群计算系统环境，包括：

图像预处理模块；

图像IHS并行融合模块；

图像PCA并行融合模块；

图像小波并行融合模块；

图像融合效果并行评估模块，

所述图像预处理模块和图像融合效果评估模块均与所述图像IHS并行融合模块、图像PCA并行融合模块和图像小波并行融合模块相连。

各模块可以为并行集群处理系统中的多个结点(计算机或服务器)，实现并行处理。

本发明还提供了一种影像快速融合方法，包括如下步骤：

(1)图像预处理：对原始多光谱图像和全色图像进行配准并调整为大小相同的图像；

(2)图像并行融合，包括：

a、将预处理后的多光谱图像和全色图像进行IHS并行融合，保存融合结果，并将融合结果与原始多光谱图像进行图像IHS并行融合性能评估；

b、将预处理后的多光谱图像和全色图像进行PCA并行融合，保存融合结果，并将融合结果与原始多光谱图像进行图像PCA并行融合性能评估；

c、将预处理后的多光谱图像和全色图像进行小波并行融合，保存融合结果，并将融合结果与原始多光谱图像进行图像小波并行融合性能评估。

所述IHS并行融合，其并行融合处理步骤为：

1)根据进程数对多光谱图像和全色图像数据进行分块；

2)对结点上的小块多光谱图像和全色图像数据进行IHS融合；

3)将结点上的融合后的分块融合结果回收，组合形成融合图像；

其中步骤2)中的IHS融合方法为IHS融合中的三角变换融合，三角变换计算公式为：

$I=\frac{1}{3}(R+G+B)---\left(1\right)$

当B＝Min时， $H=\frac{(G-B)}{3(I-B)},S=1-\frac{B}{I}---\left(2\right)$

当R＝Min时， $H=\frac{(B-R)}{3(I-R)},S=1-\frac{R}{I}---\left(3\right)$

当G＝Min时， $H=\frac{(R-G)}{3(I-G)},S=1-\frac{G}{I}---\left(4\right)$

所述图像PCA并行融合的主要步骤为：

1)根据进程数进行图像数据分配，以此分别并行计算卫星多光谱遥感影像的全局R、G、B分量的均值；

2)根据下式并行计算多光谱影像的协方差矩阵；

${\mathrm{\σ}}_{\mathrm{xy}}=\mathrm{COV}(X,Y)=\frac{\underset{i=1}{\overset{\mathrm{Size}}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})}{\mathrm{Size}-1}$

$=[\underset{i=1}{\overset{\mathrm{Size}/n}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})+\underset{\mathrm{Size}/n+1}{\overset{2*\mathrm{Size}/n}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})+\·\·\·\underset{\mathrm{Size}*(n-1)/n+1}{\overset{\mathrm{Size}}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})]/(\mathrm{Size}/1)$

X_i，Y_i分别代表第i个像素R、G、B波段中的其中的一个值；

3)协方差矩阵特征值与特征向量计算，按特征值绝对值大小进行特征向量排序；

4)采用并行计算进行多光谱影像数据的PCA变换，并将变换后的第一主成分为卫星全色遥感图像替换；

5)采用并行计算进行PCA反变换；

6)形成融合图像。

所述并行小波影像融合算法主要步骤为：

1)对配准后的多光谱影像MImage进行并行计算IHS变换，得到分量I_M、H_M、S_M，对分量I_M和全色图像进行并行计算两次塔式小波变换；

2)对分量I_M和全色图像进行小波变换后的所得的相应小波系数按特征匹配度规则进行并行融合；

其中特征匹配规则为：

$E(x,y,k)=\frac{1}{M*N}\underset{x,y=0}{\overset{M*N-1}{\mathrm{\Σ}}}{\left|c(x,y,k)\right|}^2$

$F_{\mathrm{PM}}(x,y,k)=\frac{2\underset{M}{\mathrm{\Σ}}\underset{N}{\mathrm{\Σ}}{C}_P(x+M,y+N,k)C_M(x+M,y+N,k)}{E_P(x,y,k)E_M(x,y,k)}$

C_F(x，y，k)＝W_PC_P(x，y，k)+W_MC_M(x，y，k)

当F_PM≤T时，有Wmin＝0，Wmax＝1

若E_P＞E_M，W_P＝Wmax，W_M＝Wmin，否则W_M＝Wmax，W_P＝Wmin，

当F_PM＞T时，有Wmin＝1/2-(1-F_PM)/2(1-T)，Wmax＝1-Wmin

若E_P＞E_M，W_P＝Wmax，W_M＝Wmin，否则W_M＝Wmax，W_P＝Wmin；

3)将上述融合后的小波系数进行并行小波逆变换，得到融合量I_F

4)将量I_F、H_M、S_M进行并行反IHS变换，即获得融合后的图像。

所述并行融合效果评估采用标准差、偏差指数、平均梯度和相关系数四个评价因子进行并行计算，对融合结果进行并行定量分析。

所述标准差并行计算步骤为：

(1)根据进程数对融合后的图像进行数据分块；

(2)对结点上的小块融合图像进行像素各波段值求和，然后对这些结点上的和回收，并求总和，同时求平均，求得融合图像的各波段均值；

(3)对结点上的小块融合图像各像素与融合图像全局均值求差平方和，将所有结点上的差平方和求总和，总和与(M*N-1)的商求平方根，即为标准差；

标准差公式为：

$\mathrm{std}=\sqrt{\frac{1}{\mathrm{MN}-1}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(I(i,j)-\stackrel{\‾}{I})}^2}$

式中，I(i，j)、F(i，j)分别代表融合后的图像在像素(i，j)处的值，A(i，j)代表原始多光谱图像在像素(i，j)处的值，M、N为融合后的图像和原始图像的高、宽。

所述偏差并行计算步骤为：

(1)根据进程数对融合后的图像及原始多光谱图像进行数据分块；

(2)对结点上的小块融合图像及原始多光谱图像按公式下式求商，并对各结点上的商求和；

求商公式为 $\mathrm{DI}=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{|F(i,j)-A(i,j)|}{A(i,j)}$

(4)将所有结点上的和回收求总和，总和与M*N的商，即为偏差指数。

所述平均梯度的并行计算步骤为：

(1)根据进程数进行对融合后的图像进行数据分块；每次分块需从下一块中多取一行图像数据，

(2)对结点上的小块融合图像各像素按公式求平方和，并对平方和求和，

(3)将所结点上的和回收求总和，总和与M*N的商，即为平均梯度，

平均梯度计算公式为：

$\left\{\begin{array}{c}{\mathrm{\Δ}}_{x}f(i,j)=I(i,j)-I(i-1,j)\\ {\mathrm{\Δ}}_{y}f(i,j)=I(i,j)-I(i,j-1)\end{array}\right.$

$\▿\stackrel{\‾}{g}=\frac{1}{(M-1)(N-1)}\sqrt{\underset{i=2}{\overset{M}{\mathrm{\Σ}}}\underset{j=2}{\overset{N}{\mathrm{\Σ}}}({\mathrm{\Δ}}_x{f}^2(i,j)+{\mathrm{\Δ}}_y{f}^2(i,j))}$

所述相关系数的并行计算步骤为：

(1)根据进程数对融合后的图像及原始多光谱图像进行图像数据分配，以此分别并行计算对融合后的图像及原始多光谱图像的全局R、G、B分量的均值；

(2)根据相关系数公式并行计算相关系数各分子和分母中两部分的值；

(3)对各结点分子和分母两部分的值分别总体求和，得到相关系数总的分子和分母两部分的值；

(4)对总的分子和分母值求商，得到所求的相关系数；

相关系数公式为

$\mathrm{corr}(F,A)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(F(i,j)-\stackrel{\‾}{F})(A(i,j)-\stackrel{\‾}{A})}{\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(F(i,j)-\stackrel{\‾}{F})}^2\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(A(i,j)-\stackrel{\‾}{A})}^2}}$

本发明的有益效果为：依赖于并行集群计算机系统的并行处理能力，对多光谱影像和全色影像进行图像IHS并行融合处理和并行评估、图像PCA并行融合处理和并行评估、图像小波并行融合处理和并行评估，将计算工作分配给不同结点并行处理，实现了快速融合评估，解决了大数据量影像融合速度慢的问题，可用于遥感影像融合，为灾害预防及救灾、军事等领域的遥感应用提供了保障。对于其他类型的高分辨率图像和多波段彩色图像，同样可以使用本发明的方法来实现。

附图说明

图1为本发明的影像快速融合系统示意图；

图2为本发明的遥感影像快速融合方法流程图；

图3为本发明的遥感影像快速融合方法的并行IHS融合流程图；

图4为本发明的并行PCA融合流程图；

图5为本发明的并行小波融合流程图；。

图6为标准差并行计算流程图；

图7为平均梯度并行计算流程图；

图8为偏差指数并行计算流程图。

具体实施方式

下面结合附图，以遥感影像快速融合为例，对本发明的实现方法及系统作进一步的说明。

一、卫星遥感影像快速融合系统

如图1所示，本发明的影像快速融合系统，为并行集群计算系统，包括：

图像预处理模块，为并行集群计算系统中结点计算机，主要用于多光谱图像和全象图像的配准和调整，即调整为大小相同的图像，当然也可以进行除噪、几何校正等预处理；

图像IHS并行融合模块，为并行集群计算系统中多个结点计算机，以进行图像的IHS并行融合处理，得到IHS融合图像；

图像PCA并行融合模块，为并行集群计算系统中多个结点计算机，以进行图像的PCA并行融合处理，得到PCA融合图像；

图像小波并行融合模块，为并行集群计算系统中多个结点计算机，以进行图像的小波并行融合处理，得到小波融合图像；

图像融合效果并行评估模块，为并行集群计算系统中多个结点计算机，对融合图像进行并行融合性能评估。

所述图像预处理模块和图像融合效果评估模块均与所述图像IHS并行融合模块、图像PCA并行融合模块和图像小波并行融合模块通过网络通信器件相连。

二、卫星遥感影像快速融合方法

基于并行集群计算系统的卫星遥感影像快速融合方法具体描述如下：

1、利用图像预处理模块，对遥感卫星获得的多光谱图像和全色图像进行几何校正、去噪和配准，并调整到大小一致，以便为图像融合作准备。

2、利用图像融合模块，对多光谱图像和全色图像进行融合处理，采用了IHS并行融合、PCA并行融合和小波并行融合。

IHS并行融合，在提高融合图像信息含量以及与原始图像相关性方面具有良好的性能。

PCA并行融合，在保持图像的清晰度方面有优势，光谱信息损失比IHS方法稍好。后续应用若需要图像有更好的光谱特性时，PCA变换是较IHS变换更好的选择。

小波并行融合，在保持图像光谱特性方面比较有优势。

由于图像融合过程中，数据处理和计算量大，本发明采用将图像分块，以充分利用并行集群计算系统多结点计算机并行处理的高速性能，进行并行处理，提高融合处理速度。

融合处理分图像IHS并行融合模块、图像PCA并行融合模块、图像小波并行融合模块；如图2所示，并行融合步骤为：

步骤1：用户准备多光谱影像和高分辨全色影像进行预处理，得到图2中配准并调整好的多光谱影像和高分辨全色影像，以为融合处理作准备。

步骤2：配置好所用计算机，选择进程数，向三个融合模块中均输入上述两幅对应的图像，进行图像IHS并行融合处理、图像PCA并行融合处理和图像小波并行融合处理，分别得到IHS并行融合结果、PCA并行融合结果和小波并行融合结果并保存。融合结果再通过图像融合性能并行评估模块并行评估，分别得到IHS并行融合结果并行评估结果、图像PCA并行融合结果并行评估结果、小波并行融合结果并行评估结果。融合处理和评估均采用并行处理，因此使得本发明影像融合处理速度快，以解决大数据量遥感影像融合速度慢的问题。

各并行融合过程分别介绍如下：

1、图像IHS并行融合，如图3所示，步骤为：

1)根据进程数对多光谱遥感图像和高分辨全色图像数据进行分块，获得与进程结点对应的各进程结点分块多光谱影像、各进程结点分块高分辨全色影像；

2)在各结点上，计算机对各结点对应的分块多光谱影像与高分辨全色影像通过IHS变换法进行融合，通过结点上的小块多光谱遥感图像和全色图像数据进行IHS融合，在各结点上获得各结点对应的分块多光谱影像与高分辨全色影像进行融合所得的分块融合结果；

3)将各结点上的融合后的分块融合结果回收收集，组合形成融合结果图像；

其中步骤2)中的IHS融合过程为IHS融合中的三角变换融合，三角变换计算公式为：

$I=\frac{1}{3}(R+G+B)---\left(1\right)$

$H=\frac{(G-B)}{3(I-B)},S=1-\frac{B}{I}$

当B＝Min时，(2)

当R＝Min时， $H=\frac{(B-R)}{3(I-R)},S=1-\frac{R}{I}---\left(3\right)$

当G＝Min时， $H=\frac{(R-G)}{3(I-G)},S=1-\frac{G}{I}---\left(4\right)$

。

2、如图4所示，图像PCA并行融合的主要步骤为：

1)根据进程数对多光谱影像进行图像数据分配，构建多光谱影像输入矩阵IM，分别并行计算卫星多光谱遥感影像的全局R、G、B分量的均值；

2)根据下式并行计算多光谱影像输入矩阵I的协方差矩阵IX；

${\mathrm{\σ}}_{\mathrm{xy}}=\mathrm{COV}(X,Y)=\frac{\underset{i=1}{\overset{\mathrm{Size}}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})}{\mathrm{Size}-1}$

$=[\underset{i=1}{\overset{\mathrm{Size}/n}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})+\underset{\mathrm{Size}/n+1}{\overset{2*\mathrm{Size}/n}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})+\·\·\·\underset{\mathrm{Size}*(n-1)/n+1}{\overset{\mathrm{Size}}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})]/(\mathrm{Size}/1)$

X_i，Y_i分别代表第i个像素R、G、B波段中的其中的一个值；

3)协方差矩阵IX的特征值m与特征向量V计算，按特征值m的绝对值大小进行特征向量V排序，得到Vp；

4)采用并行计算进行多光谱影像数据的PCA变换，并将变换后的第一主成分为卫星全色遥感图像替换，

具体为：采用并行计算进行多光谱影像数据的PCA变换，即将Vp′*IM，得到主成份矩阵P，其中对应于最大特征值的一行P1为第一主成份；用经过拉伸后的全色影像数据替换P1，得到新的矩阵Pnew矩阵；

5)采用并行计算进行PCA反变换，变换方法为Vp*Pnew，得到融合结果矩阵IMF；

6)形成融合图像。

3、图像小波并行融合的主要步骤如图5所示，为：

1)对配准后的多光谱影像根据进程数进行并行IHS变换，得到分量I_M、H_M、S_M，对分量I_M和高分辨全色图像进行并行计算两次塔式小波变换；

2)对分量I_M和全色高分辨图像进行小波变换后的所得的相应小波分解系数按特征匹配度规则进行并行融合；

特征匹配度规则

$E(x,y,k)=\frac{1}{M*N}\underset{x,y=0}{\overset{M*N-1}{\mathrm{\Σ}}}{\left|c(x,y,k)\right|}^2$

$F_{\mathrm{PM}}(x,y,k)=\frac{2\underset{M}{\mathrm{\Σ}}\underset{N}{\mathrm{\Σ}}{C}_P(x+M,y+N,k)C_M(x+M,y+N,k)}{E_P(x,y,k)E_M(x,y,k)}$

C_F(x，y，k)＝W_PC_P(x，y，k)+W_MC_M(x，y，k)

当F_PM≤T时，有Wmin＝0，Wmax＝1

若E_P＞E_M，W_P＝Wmax，W_M＝Wmin，否则W_M＝Wmax，W_P＝Wmin，

当F_PM＞T时，有Wmin＝1/2-(1-F_PM)/2(1-T)，Wmax＝1-Wmin

若E_P＞E_M，W_P＝Wmax，W_M＝Wmin，否则W_M＝Wmax，W_P＝Wmin。

3)将上述融合后的小波系数进行并行小波逆变换，得到融合分量I_F；

4)将分量I_F、H_M、S_M进行并行IHS反变换，即获得融合后的图像。

在图像融合效果并行评估模块，进行并行融合结果的并行性能评估，包括标准差、偏差指数、平均梯度和相关系数并行计算，定量分析。标准差反映了图像中每个像素的灰度值相对平均灰度值的离散程度，标准差越大，说明图像灰度级分布分散，图像反差大，可以看出更多信息。偏差指数，即融合图像和低分辨率图像差值的绝对值与低分辨率图像值的比值，偏差指数越小，融合结果越好。平均梯度反映图像中微小细节信息反差与纹理变化特征，也反映了图像的清晰程度，平均梯度越大，表明图像的细节信息越清晰。相关系数反映了图像X和Y光谱特征的相关程度，通过比较融合前后图像的相关系数，可以看出光谱信息改变程度。

在下述融合的效果评估各公式中，I(i，j)、F(i，j)分别代表融合后的图像在像素(i，j)处的值，A(i，j)代表原始多光谱图像在像素(i，j)处的值，M、N为融合后的图像和原始图像的高、宽。

如图5所示，并行标准差，其计算步骤为：

(1)根据选择进程数对融合结果影像(融合后的图像)进行数据分块；

(2)并行计算融合结果图像各波段均值：

具体为对结点上的小块融合图像进行像素各波段(如R、G、B)值求和，然后对这些结点上的和回收，并求总和，同时求平均，求得融合图像的各波段均值(如R_mean、G_mean、B_mean)。

(3)并行计算融合结果图像的标准差，具体为对结点上的小块融合图像各像素与融合图像全局均值求差平方和，将所有结点上的差平方和求总和，总和与(M*N-1)的商求平方根，即为标准差。

标准差公式

$\mathrm{std}=\sqrt{\frac{1}{\mathrm{MN}-1}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(I(i,j)-\stackrel{\‾}{I})}^2}$

如图6所示，并行偏差指数计算步骤为：

(1)根据进程数对融合结果影像(融合后的图像)及原始多光谱图像进行数据分块；

(2)对结点上的小块融合图像及原始多光谱图像按偏差指数公式求商，并对各结点上的商求和，

(3)将所结点上的和回收求总和，总和与M*N的商，即为偏差指数。

偏差指数公式

$\mathrm{DI}=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{|F(i,j)-A(i,j)|}{A(i,j)}$

如图7所示，并行平均梯度计算步骤为：

(1)根据进程数对融合后的图像进行数据分块；每次分块需从下一块中多取一行图像数据，以便于从平均梯度公式计算平均梯度，

(2)对结点上的小块融合图像各像素按平均梯度公式求平方和，并对平方和求和，

(3)将所结点上的和回收求总和，总和的平方根与M*N的商，即为平均梯度。

平均梯度公式

$\left\{\begin{array}{c}{\mathrm{\Δ}}_{x}f(i,j)=I(i,j)-I(i-1,j)\\ {\mathrm{\Δ}}_{y}f(i,j)=I(i,j)-I(i,j-1)\end{array}\right.$

$\▿\stackrel{\‾}{g}=\frac{1}{(M-1)(N-1)}\sqrt{\underset{i=2}{\overset{M}{\mathrm{\Σ}}}\underset{j=2}{\overset{N}{\mathrm{\Σ}}}({\mathrm{\Δ}}_x{f}^2(i,j)+{\mathrm{\Δ}}_y{f}^2(i,j))}$

如图8所示，并行相关系数，其计算步骤为：

(1)根据进程数对融合结果图像(融合后的图像)及原始多光谱图像进行图像数据分配，以此分别并行计算对融合后的图像及原始多光谱图像的全局R、G、B分量的均值；

(2)根据相关系数公式并行计算相关系数各分子和分母的值；

(3)对各结点分子和分母的值分别总体求和，得到相关系数总的分子和分母值；

(4)对总的分子和分母值求商，得到所求的相关系数；

相关系数公式

$\mathrm{corr}(F,A)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(F(i,j)-\stackrel{\‾}{F})(A(i,j)-\stackrel{\‾}{A})}{\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(F(i,j)-\stackrel{\‾}{F})}^2\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(A(i,j)-\stackrel{\‾}{A})}^2}}$

以上举例的形式对本发明在遥感图像快速融合中的应用进行了详细的描述，但是本发明并不局限于遥感图像，对于其他类型的高分辨率图像和多波段彩色图像，同样可以使用本发明的方法来实现。

综上所述，依照本发明的并行遥感图像融合实现方法及系统能够进行遥感图像快速融合的操作，从而解决了大规模遥感图像融合速度慢的问题；另外本发明中的图像融合效果并行评估模块，可快速对遥感图像融合的效果进行快速评估。

Claims (10)

1.一种影像快速融合系统，该系统基于并行集群计算系统环境，包括：

图像预处理模块；

图像IHS并行融合模块；

图像PCA并行融合模块；

图像小波并行融合模块；

图像融合效果并行评估模块，

所述图像预处理模块和图像融合效果并行评估模块均与所述图像IHS并行融合模块、图像PCA并行融合模块和图像小波并行融合模块相连。

2.一种影像快速融合方法，包括如下步骤：

(1)图像预处理：对原始多光谱图像和全色图像进行配准并调整为大小相同的图像；

(2)图像并行融合，包括：

a、图像IHS并行融合：将预处理后的多光谱图像和全色图像进行IHS并行融合处理，并保存融合结果，再将融合结果与原始多光谱图像进行图像IHS并行融合性能评估；

b、图像PCA并行融合：将预处理后的多光谱图像和全色图像进行PCA并行融合处理，保存融合结果，再将融合结果与原始多光谱图像进行图像PCA并行融合性能评估；

c、图像小波并行融合：将预处理后的多光谱图像和全色图像进行小波并行融合，保存融合结果，并将融合结果与原始多光谱图像进行图像小波并行融合性能评估。

3.如权利要求2所述的影像快速融合方法，其特征在于：所述图像IHS并行融合主要步骤为：

1)根据进程数对多光谱遥感图像和全色图像数据进行分块；

2)对结点上的小块多光谱遥感图像和全色图像数据进行IHS融合；

3)将结点上的融合后的分块融合结果回收，组合形成融合图像；

其中步骤2)中的IHS融合过程为IHS融合中的三角变换融合，三角变换计算公式为：

$I=\frac{1}{3}(R+G+B)---\left(1\right)$

当B＝Min时， $H=\frac{(G-B)}{3(I-B)},S=1-\frac{B}{I}---\left(2\right)$

当R＝Min时， $H=\frac{(B-R)}{3(I-R)},S=1-\frac{R}{I}---\left(3\right)$

当G＝Min时， $H=\frac{(R-G)}{3(I-G)},S=1-\frac{G}{I}---\left(4\right)$

4.如权利要求2所述的影像快速融合方法，其特征在于：所述图像PCA并行融合的主要步骤为：

1)根据进程数进行图像数据分配，以此分别并行计算卫星多光谱遥感影像的全局R、G、B分量的均值；

2)根据下式并行计算多光谱影像的协方差矩阵；

${\mathrm{\σ}}_{\mathrm{xy}}=\mathrm{COV}(X,Y)=\frac{\underset{i=1}{\overset{\mathrm{Size}}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})}{\mathrm{Size}-1}$

$=[\underset{i=1}{\overset{\mathrm{Size}/n}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})+\underset{\mathrm{Size}/n+1}{\overset{2*\mathrm{Size}/n}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})+\·\·\·\underset{\mathrm{Size}*(n-1)/n+1}{\overset{\mathrm{Size}}{\mathrm{\Σ}}}(X_i-\stackrel{\‾}{X})(Y_i-\stackrel{\‾}{Y})]/(\mathrm{Size}/1)$

X_i，Y_i分别代表第i个像素R、G、B波段中的其中的一个值；

3)协方差矩阵特征值与特征向量计算，按特征值绝对值大小进行特征向量排序；

4)采用并行计算进行多光谱影像数据的PCA变换，并将变换后的第一主成分为全色图像替换；

5)采用并行计算进行PCA反变换；

6)形成融合图像。

5.如权利要求2所述的影像快速融合方法，其特征在于：所述图像小波并行融合的主要步骤为：

1)对配准后的多光谱影像进行并行计算IHS变换，得到分量I_M、H_M、S_M，对分量I_M和全色图像进行并行计算两次塔式小波变换；

2)对分量I_M和全色图像进行小波变换后的所得的相应小波系数按特征匹配度规则进行并行融合；

其中特征匹配规则为：

$E(x,y,k)=\frac{1}{M*N}\underset{x,y=0}{\overset{M*N-1}{\mathrm{\Σ}}}{\left|c(x,y,k)\right|}^2$

$F_{\mathrm{PM}}(x,y,k)=\frac{2\underset{M}{\mathrm{\Σ}}\underset{N}{\mathrm{\Σ}}{C}_P(x+M,y+N,k)C_M(x+M,y+N,k)}{E_P(x,y,k)E_M(x,y,k)}$

C_F(x，y，k)＝W_PC_P(x，y，k)+W_MC_M(x，y，k)

当F_PM≤T时，有Wmin＝0，Wmax＝1

若E_P＞EM，W_P＝Wmax，W_M＝Wmin，否则W_M＝Wmax，W_P＝Wmin，

当F_PM＞T时，有Wmin＝1/2-(1-F_PM)/2(1-T)，Wmax＝1-Wmin

若E_P＞E_M，W_P＝Wmax，W_M＝Wmin，否则W_M＝Wmax，W_P＝Wmin；

3)将上述融合后的小波系数进行并行小波逆变换，得到融合量I_F；

4)将量I_F、H_M、S_M进行并行反IHS变换，即获得融合后的图像。

6.如权利要求2所述的影像快速融合方法，其特征在于：所述并行融合效果评估采用标准差、偏差指数、平均梯度和相关系数四个评价因子进行并行计算，对融合结果进行并行定量分析。

7.如权利要求6所述的影像快速融合方法，其特征在于：所述标准差并行计算步骤为：

(1)根据进程数对融合后的图像进行数据分块；

(2)对结点上的小块融合图像进行像素各波段值求和，然后对这些结点上的和回收，并求总和，同时求平均，求得融合图像的各波段均值；

(3)对结点上的小块融合图像各像素与融合图像全局均值求差平方和，将所结点上的差平方和求总和，总和与(M*N-1)的商求平方根，即为标准差；

标准差公式为：

$\mathrm{std}=\sqrt{\frac{1}{\mathrm{MN}-1}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(I(i,j)-\stackrel{\‾}{I})}^2}$

式中，I(i，j)、F(i，j)分别代表融合后的图像在像素(i，j)处的值，A(i，j)代表原始多光谱图像在像素(i，j)处的值，M、N为融合后的图像和原始图像的高、宽。

8.如权利要求6所述的影像快速融合方法，其特征在于：所述偏差并行计算步骤为：

(1)根据进程数对融合后的图像及原始多光谱图像进行数据分块；

(2)对结点上的小块融合图像及原始多光谱图像按公式下式求商，并对各结点上的商求和；

求商公式为 $\mathrm{DI}=\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{|F(i,j)-A(i,j)|}{A(i,j)}$

(3)将所有结点上的和回收求总和，总和与M*N的商，即为偏差指数。

9.如权利要求6所述的影像快速融合方法，其特征在于所述平均梯度的并行计算步骤为：

(1)根据进程数进行对融合后的图像进行数据分块；每次分块需从下一块中多取一行图像数据，

(2)对结点上的小块融合图像各像素按公式求平方和，并对平方和求和，

(3)将所结点上的和回收求总和，总和与M*N的商，即为平均梯度，平均梯度计算公式为：

$\left\{\begin{array}{c}{\mathrm{\Δ}}_{x}f(i,j)=I(i,j)-I(i-1,j)\\ {\mathrm{\Δ}}_{y}f(i,j)=I(i,j)-I(i,j-1)\end{array}\right.$

$\▿\stackrel{\‾}{g}=\frac{1}{(M-1)(N-1)}\sqrt{\underset{i=2}{\overset{M}{\mathrm{\Σ}}}\underset{j=2}{\overset{N}{\mathrm{\Σ}}}({\mathrm{\Δ}}_x{f}^2(i,j)+{\mathrm{\Δ}}_y{f}^2(i,j))}$

10.根据权利要求6所述的影像快速融合方法，其特征在于所述相关系数的并行计算步骤为：

(1)根据进程数对融合后的图像及原始多光谱图像进行图像数据分配，以此分别并行计算融合后的图像及原始多光谱图像的全局R、G、B分量的均值；

(2)根据相关系数公式并行计算相关系数各分子和分母中两部分的值；

(3)对各结点分子和分母两部分的值分别总体求和，得到相关系数总的分子和分母两部分的值；

(4)对总的分子和分母值求商，得到所求的相关系数；

相关系数公式为

$\mathrm{corr}(F,A)=\frac{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}(F(i,j)-\stackrel{\‾}{F})(A(i,j)-\stackrel{\‾}{A})}{\sqrt{\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(F(i,j)-\stackrel{\‾}{F})}^2\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{(A(i,j)-\stackrel{\‾}{A})}^2}}.$

Images