基于BEMD的SAR影像与多光谱影像的融合方法

摘要

本发明一种SAR影像与多光谱影像的融合方法，包括：对多光谱影像进行IHS变换得到I、H、S分量；对SAR影像和I分量分别进行BEMD变换得到SAR影像和I分量的IMF分量和残余分量；调整SAR影像和I分量的IMF分量和残余分量，使二者的IMF分量个数保持一致，并将多余的IMF分量加入对应的残余分量；分别将调整后的SAR影像和I分量的IMF分量与IMF分量，以及残余分量与残余分量进行融合；根据融合结果进行BEMD逆变换，得到融合后的I分量；对融合后的I分量、H分量和S分量进行IHS逆变换，得到融合后的影像。本发明的基于BEMD的SAR影像与多光谱影像的融合方法保存了较好的纹理信息，能够对于同类地物的具体细分提供可靠保证，提高了对非线性、非平稳信号的分析及处理能力。

说明书

技术领域

本发明涉及影像处理领域，尤其涉及一种基于BEMD的SAR影像与多光谱影像的融合方法。

背景技术

随着传感器技术、遥感技术等的迅猛发展，多平台、多层面、多时相、多传感器、多光谱和多分辨率的遥感影像构成了多区域的多源遥感信息。不同的传感器具有不同的成像机理且工作于不同的波长范围，因此所获得的影像往往反映地物不同方面的特征。通过单一传感器所提供的信息可能是不全面、不一致、甚至是不准确的。为了更有效的利用日益复杂的多源遥感影像信息，与之相关的影像融合技术应运而生。

合成孔径雷达(Synthetic Aperture Radar,SAR)是一种微波传感器，其通过主动发射电磁波并接收返回的电磁波信息实现遥感成像。相对于光学传感器来说，合成孔径雷达具有全天时、全天候、高分辨率的成像能力，并且其作用距离远，成像范围广，穿透能力强，解决了光学遥感传感器受时间及天气的影响问题。

目前，常用的融合方法有：IHS变换、PCA变换、Brovey变换以及小波变换，其中，IHS变换、PCA变换和Brovey变换通常会导致比较严重的光谱扭曲，而小波变换则解决了光谱扭曲的问题。然而，上述方法主要用于多光谱影像及全色影像的融合，在本文提出的对光学影像与高分辨率SAR影像进行融合时，采用上述方法会丢失大量纹理信息，在对细微地物进行分类提取时有很大难度。

下面，对上述的遥感影像融合领域的现有技术做一个较为详细的描述，以便更好地理解本发明。

传统的遥感影像融合一般采用光学全色影像与多光谱影像进行融合，传统的融合方法包括：IHS变换、主成分变换(PCA)、比值融合法、主成分分析法(K-L变换)、加权融合法、小波变换法等。

以IHS变换为例：将由光学、热红外和雷达(微波)等方式得到的不同波段遥感数据合成的一个对物体颜色属性进行描述的RGB颜色空间变换到以亮度(Intensity)、色度H(Hue)、饱和度S(Saturation)来描述影像的IHS彩色空间。

IHS变换的融合算法如下：

正变公式：

$\left(\begin{array}{c}I\\ v_1\\ v_2\end{array}\right)=\left(\begin{array}{ccc}\frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}\\ \frac{1}{\sqrt{6}}& \frac{1}{\sqrt{6}}& -\frac{2}{\sqrt{6}}\\ \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0\end{array}\right)*\left(\begin{array}{c}R\\ G\\ B\end{array}\right)$

$S=\sqrt{v_1^2+v_2^2}$

H＝arctg(v₁/v₂)

式中：I表示亮度，H表示色度，S表示饱和度，v1、v2是为了计算I、H引入的中间变量。

反换公式：

$\left(\begin{array}{c}R\\ G\\ B\end{array}\right)=\left(\begin{array}{ccc}\frac{1}{\sqrt{3}}& \frac{1}{\sqrt{6}}& \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{6}}& -\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{6}}& 0\end{array}\right)*\left(\begin{array}{c}I\\ v_1\\ v_2\end{array}\right)$

以小波变换为例：小波变换融合的过程如下：

小波变换定义为：

变换核函数为：

其中，为一个基本小波，又称母小波，或者是小波基，它是以t＝0为为中心的带通函数，且时域平均值

小波变换在于，以空间大尺度处理低频信息，以空间小尺度处理高频信息。

目前，SAR影像与光学影像的融合方法中采用的多尺度分析方法主要包括：二维离散小波变换、平移不变离散小波变换以及表示各向异性的曲波变换等。本领域技术人员应当理解，小波变换在本质上是一种线性变换，对非线性，非平稳信号的处理能力有限；而且，小波变换是非适应的，其融合效果取决于小波基的选择。

因此，需要一种将高分辨率SAR影像与光学影像进行融合的方法，该方法能够克服现有方法在处理非线性、非平稳信号上较为乏力的缺陷，并且同时能够保持多源影像中的光谱信息、边缘信息及纹理信息。

本说明书中的部分关键术语的缩略语的意义在此被示出，包括：EMD一维经验模态分解(Empirical Mode Decomposition,EMD)；二维经验模态分解(BidimensionalEmpirical Mode Decomposition,BEMD)；本征模函数(Intrinsic Mode Function,IMF)；合成孔径雷达((SyntheticApertureRadar,SAR)。

发明内容

本发明的目的是提供一种能够克服上述缺陷的合成孔径雷达SAR影像与多光谱影像的融合方法。

在第一方面，本发明提供了一种合成孔径雷达SAR影像与多光谱影像的融合方法，其特征在于：对多光谱影像进行IHS变换，得到所述多光谱影像的I分量、H分量和S分量；对所述SAR影像和所述I分量分别进行二维经验模态分解BEMD变换，得到所述SAR影像的本征模函数IMF分量和残余分量，以及所述I分量的IMF分量和残余分量；对所述SAR影像的IMF分量和残余分量，以及所述I分量的IMF分量和残余分量进行调整，使所述SAR影像的IMF分量与所述I分量的IMF分量个数保持一致，并将多余的IMF分量加入对应的残余分量；分别将所述调整后的SAR影像的IMF分量与I分量的IMF分量，以及调整后的SAR影像的残余分量与I分量的残余分量进行融合；根据融合结果进行BEMD逆变换，得到融合后的I分量；对融合后的I分量、所述H分量和所述S分量进行IHS逆变换，得到融合后的影像。

优选地，所述BEMD变换的步骤包括：对待进行BEMD变换的影像进行初始化；识别初始化后的影像中的极值点；在所述极值点满足预定阈值的情况下，通过对所述极值点进行拟合，得到影像的平均包络面；根据所述平均包络面提取影像的局部趋势；在所述局部趋势满足预定的IMF条件的情况下，将所述局部趋势判定为IMF分量；根据所述IMF分量计算残余分量；在所述残余分量满足预定单调条件的情况下，将所述残余分量判定为最终的残余分量。

优选地，在所述局部趋势不满足IMF判定条件的情况下，用所述局部趋势代替初始化后的影像，继续寻找其中的极值点；在所述残余分量不满足单调条件的情况下，用所述残余分量代替初始化后的影像，继续寻找其中的极值点。

优选地，所述IMF判定条件根据在其中识别极值点的影像中的极值点与零点之间的关系以及所述平均包络面的值决定；所述单调条件通过两个相继判定出来的IMF分量计算得到。

优选地，在对多光谱影像进行IHS变换的步骤之前，还包括：对所述SAR影像进行去噪；通过重采样统一去噪后的SAR影像和所述多光谱影像的像素分辨率；对重采样之后的SAR影像和多光谱影像进行影像配准。

优选地，在所述融合的步骤前，还包括：对所述调整后的SAR影像的IMF分量与I分量的IMF分量进行拉普拉斯高通滤波，并且对调整后的SAR影像的残余分量与I分量的残余分量进行拉普拉斯低通滤波；基于所述滤波的结果分别计算所述SAR影像中的IMF分量和残余分量中各像素的权重；其中，所述融合基于所述权重进行。

本发明的基于BEMD的SAR影像与多光谱影像的融合方法保存了较好的纹理信息，能够对于同类地物的具体细分提供可靠保证，提高了对非线性、非平稳信号的分析及处理能力。

附图说明

图1为根据本发明实施例的基于BEMD的SAR与多光谱影像的融合方法的流程图；

图2为图1中的BEMD变换的流程图；

图3为基于小波变换的桂林地区的多源遥感影像的融合结果图；

图4为基于本发明的融合方法的桂林地区的多源遥感影像的融合结果图。

具体实施方式

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

针对现有的遥感影像融合变换方法在处理非线性、非平稳信号上较为乏力的缺陷，本发明提出一种新的遥感影像融合变换方法。BEMD变换适用于非线性、非平稳信号的分析，具有较高信噪比；另外，BEMD依据信号自身的时间初读特征实现信号分解，最终融合效果不依赖任何基函数的选择，具有完全的自适应性。利用BEMD变换方法，本发明实现了SAR影像与高分光学影像的融合。

首先，将本发明所采用的影像融合方法笼统地描述如下几个步骤。

1)对多源遥感影像融合的对象影像进行预处理。

例如，可以选取以TerraSAR-X卫星获取的SAR影像和GF-1相机获取的多光谱影像进行影像融合。优选地，为保证融合效果，分别先对两种影像数据进行地理定标、几何校正、影像配准、滤波去噪等预处理。由于两种影像具有不同的空间分辨率，在进行影像配准之前需统一影像的像素分辨率；另外，由于SAR影像在获取时存在入射角，需在预处理过程中消除叠影。

2)针对影像融合过程中的计算循环等问题，合理地选取影像像素和影像窗口的大小。影像像素大小可以选取8的倍数，优选地以正方形为最佳；且影像不宜过大，例如可以选取影像窗口大小为4000*4000。

3)对多光谱影像和SAR影像进行BEMD分解。

多光谱影像含有三个波段R、G、B，而IHS彩色空间三分量具有相对独立性，能够更好的体现色彩信息。因此，对多光谱影像先进行IHS变换，后进行BEMD分解，即对其三个波段分别进行分解。在实验中，可以定义分解得到四层分量，分别为三个IMF分量和一个残余分量。第一IMF分量为最高频信息，第二和第三分量的频率信息依次降低，残余分量主要为亮度信息。

SAR影像只有一个波段，即灰度波段，那么也对SAR影像的灰度波段进行BEMD分解。

4)将分解后的SAR影像灰度波段与分解后的多光谱影像I分量进行融合，然后与H、S分量合并再转换为R、G、B，最后得到融合影像。

下面，详细描述根据本发明的基于BEMD高分辨率的SAR影像与光学影像的融合方法。

图1为根据本发明实施例的基于BEMD的SAR影像与多光谱影像的融合方法的流程图。

第一步：影像预处理，其过程主要包括影像重采样，SAR影像去噪，影像配准。

首先，可以通过下采样统一两幅影像的像素分辨率。

其次，由于SAR影像易受噪声影响，从而导致影像的非线性、非平稳信号随影像空间分辨率提高而愈加明显，固有噪声也随之更加突出。因此，需要在SAR影像采样前先对其进行去噪处理，例如可以采用改进Goldstein滤波对SAR影像进行降噪处理。

最后，针对去噪、采样后的SAR影像和采样后的多光谱影像，采集同名像点对进行影像配准。

第二步：将SAR影像与多光谱影像进行融合，其具体操作步骤如下：

1、对多光谱影像进行IHS变换(即图1中的彩色空间转换)，得到多光谱影像的I、H、S分量。

$\left(\begin{array}{c}I\\ v_1\\ v_2\end{array}\right)=\left(\begin{array}{ccc}\frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}\\ \frac{1}{\sqrt{6}}& \frac{1}{\sqrt{6}}& -\frac{2}{\sqrt{6}}\\ \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0\end{array}\right)*\left(\begin{array}{c}R\\ G\\ B\end{array}\right)$

$S=\sqrt{v_1^2+v_2^2}$

H＝arctg(v₁/v₂)

注：v1、v2为中间变量。

2、对预处理后的SAR影像和变换后的多光谱影像I分量分别进行BEMD变换。

通过分别的BEMD变换，得到SAR影像的IMF分量和残余分量R，并且得到多光谱影像的IMF分量和残余分量R。经BEMD分解后，SAR影像可表示为多光谱影像I分量可表示为其中，I、R分别为IMF分量与残余分量；M、N分别为SAR与多光谱影像的I分量可分解的IMF的个数。

3、对SAR影像和多光谱影像的影像分解量进行处理，使二者具有相同的分解层数。

由于BEMD具有完全自适应性，SAR影像与多光谱影像的I分量通过BEMD分解后得到的IMF个数可能不同，从而导致在后续的融合过程中的不一致。因此，需要对影像的分解量进行处理，使待融合影像的分解层数一致。具体地，可以在BEMD变换中设定分解后将得到的IMF个数，使得两种影像最终得到的IMF个数相同。例如，设n＝min{M,N}，以n＝M为例，多光谱影像的I分量的第n个IMF之后的IMFs分量将被加入到残余分量中，使其IMF个数为n。此时，多光谱影像的I分量可以表示为：

$F^{MI}=\underset{n=1}{\overset{N}{\Sigma }}{I}_n^{MI}+R_N^{MI}$

其中，

4、对SAR影像和多光谱影像通过BEMD变换得到的分量进行拉普拉斯滤波，以检测影像中的边缘信息。例如，对SAR影像和多光谱影像的IMF分量进行拉普拉斯高通滤波，对二者的残余分量进行拉普拉斯低通滤波。影像中的边缘信息主要对应高频信息，即通过BEMD分别得到的IMF分量，相应地，其中的低频信息为残余分量。

5、基于拉普拉斯算子确定的像素窗口，检测窗口区域内是否平滑。通过检测，判断上一步骤滤波后是否得到相同频率的信息。如判断为平滑，则表明高通滤波得到的都是高频信息，低通滤波得到的都是低频信息。换言之，平滑表明滤波后得到的信息在频率上没有跳跃性差异，即没有边缘。

6、在判断像素窗口内区域平滑的情况下，表示频率信息相同或相似，则可推断像素窗口内为同一地物，那么计算其平均权重即可；否则需要分别计算权重。例如，可以根据基于区域特征的融合规则来计算SAR影像各分量中各像素的权重，并将SAR影像的IMF分量与残余分量的权重分别记作α_n、β。

7、将IMF分量及残余分量分别进行不同尺度的融合，其融合公式为如下：

$I_n^{fuse}={\alpha }_n{I}_n^{SF}+(1-{\alpha }_n)I_n^{MI}$

R^fuse=βR^SF+(1-β)R^MI

6、进行BEMD逆变换，将融合结果作为新的I分量I^fuse，其中逆变换公式如下：

$I^{fuse}=\underset{k=1}{\overset{n}{\Sigma }}{I}_k^{fuse}+R^{fuse}$

7、将新的I分量I^fuse与多光谱影像的H、S分量进行彩色空间逆变换，即将I^fuseH^MS^M逆变换到RGB空间，得到融合后的影像。

图2是图1中的BEMD变换的流程图。

在步骤201，对待进行BEMD变换的影像进行初始化。

在步骤202，识别初始化后的影像的极值点，具体地，即寻找其所有的局部极大值点和极小值点。

在步骤203，判断极值点的个数是否满足预定条件。如果满足，则流程进行到204；否则流程进行到步骤207。

在步骤204，选择有效的插值算法分别拟合信号的极大值点与极小值点，得到上下包络面，并计算得到平均包络面。

在步骤205，根据平均包络面提取影像的局部趋势；

在步骤206，判断该局部趋势是否满足IMF条件。若满足，则流程进行到步骤207；否则，用该局部趋势代替步骤202中的初始化的影像，流程回到步骤202。

IMF判定条件为：

1)对于整个数据集，|极值点-零点|≤1；

2)对于数据集上任一点，局部极大值与局部极小值确定的平均包络面值为零。

在步骤207，将该局部趋势设定为IMF分量。

在步骤208，根据该IMF分量计算残余分量。

在步骤209，判断步骤207中计算得到的残余分量是否满足单调条件。如满足，则流程结束，将该残余分量判定为最终输出的残余分量；否则用该残余分量代替步骤202中的初始化影像，流程回到步骤202。

残余分量是否单调的判定标准为：

1)不包含任何IMF；或者

2)残余分量小于SD，0.2<SD<0.3，

通常，SD设置的范围为0.2-0.3，即当SD满足0.2＜SD＜0.3时，筛选过程结束。

图3为基于小波变换的桂林地区的多源遥感影像的融合结果图；图4为基于本发明的融合方法的桂林地区的多源遥感影像的融合结果图。通过图3和图4的比较可以明显看出，根据本发明的融合方法保存了较好的纹理信息，能够对于同类地物的具体细分提供可靠保证，提高了对非线性、非平稳信号的分析及处理能力。

本发明基于BEMD变换对多源遥感影像进行融合，作为遥感影像分类的基础，能够保证在处理线性或者非线性，平稳或者非平稳信号的同时，保证多源遥感影像所包含的光谱信息、边缘信息以及纹理信息不丢失。传统SAR与光学影像融合方法主要采用二维离散小波变换、平移不变小波变换等线性变换，这些方法在很大程度上放弃了非线性及非平稳信号的处理。本发明提供的融合方案更适用于非线性、非平稳信号的分析，具有较高的信噪比。另外，BEMD变换优于小波变换在于它不同于小波变换在对信号处理时依赖于小波基的选择，它是依赖于信号自身时间尺度特征来实现信号分解，具有完全自适应性。

专业人员应该还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。

结合本文中所公开的实施例描述的方法或算法的步骤可以用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明。所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围。凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。