CN101359049A - 一种遥感影像融合方法 - Google Patents

Abstract

本发明涉及一种基于多尺度通用传感器模型和Contourlet域隐马尔可夫树模型(CHMT，Contourlet－domain Hidden Markov Tree)相结合的遥感影像融合方法，通过结合Contourlet变换和隐马尔可夫树模型，在多尺度通用传感器模型的框架下，对融合的Contourlet系数进行最大后验概率估计，建立了一种基于CHMT模型的遥感影像融合方法。本发明首次将CHMT模型与多尺度通用传感器模型相结合用于遥感影像融合研究，并成功地将待融合图像进行主成份变换后的第一主成份用于真实场景的先验概率估计，在提高融合影像空间细节信息的同时，有效地保持多光谱影像的光谱信息。

Description

一种遥感影像融合方法

技术领域

本发明涉及一种遥感影像融合方法，特别涉及一种基于多尺度通用传感器模型与Contourlet域隐马尔可夫树(CHMT，Contourlet-domain Hidden MarkovTree)模型相结合的遥感影像融合方法。

背景技术

早期的图像融合主要以图像视觉增强为目的，采用代数运算法、彩色空间法等。目前常用的图像融合技术有IHS融合、Brovey融合、分量替换法、小波变换法等。随着新的数学理论和计算智能理论等被应用到图像融合领域，如小波包、Curvelet、Contourlet、模糊数学和数学形态学等的使用，很大程度上推进了图像融合理论和技术的发展。然而，这些方法的理论解释大多都停留在单纯的图像处理层面，没有一个通用传感器模型的融合方法作为理论依据。

CHMT模型是一种非常重要而有效的统计建模数学工具，它将Contourlet变换(即“轮廓波”变换)和隐马尔可夫模型(HMM，Hidden Markov Model)的各自特性集于一体，在图像处理和计算机视觉的各个方面，比如去噪，图像分割等等，起到了很好的效果。但目前还没有其在遥感图像融合领域应用的报道。

发明内容

本发明目的在于提供一种遥感影像融合方法；本发明的另一个目的在于提供一种基于多尺度通用传感器模型与CHMT模型相结合的遥感影像融合方法。

本发明目的是通过如下技术方案实现的：

一种遥感影像融合方法，基于多尺度通用传感器模型与CHMT模型相结合，得融合后的影像。

本发明遥感影像融合方法包括如下步骤：建立Contourlet域多尺度通用传感器模型，计算得到由传感器观测的图像经Contourlet变换后在相应尺度某方向子波段上对应位置处的观测值；对传感器的观测值应用主成分变换，获取第一主成分分量，作为对实际场景的先验估计，然后对其利用CHMT模型进行训练，获得条件概率，最后，计算得在Contourlet分解的某尺度某方向波段上实际场景的最大后验估计值；采用最小二乘法求得待融合影像的传感器系数，然后根据实际场景的最大后验估计值，求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即得融合后的影像。

本发明遥感影像融合方法具体包括如下步骤：

A、建立Contourlet域通用传感器模型，其可以表示为：

           o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，i＝1，2，…，q为传感器的数量，也可认为是待融合的波段数，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置，l_k是实际场景L经Contourlet变换后在k处的值，o_i，k是由传感器i观测的图像o_i经Contourlet变换后在k处的值，β_i，k和α_i，k分别是经Contourlet变换后第i个传感器在k处的增益和偏差，n_i，k表示经Contourlet变换后的噪声，是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；将公式(1)式写成向量的形式为：

                 o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

B、在通用传感器模型的理论下，根据公式(2)的成像模型，已知观测值O_k的情况下，真实场景l_k的后验概率可以表示为：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)式，求得融合后的Contourlet系数表达式；接下来采用最小二乘方法要对α_k，m和β_k，m进行估计；分别得到β_k，m和α_k，m的估计值和

最后，将

和

代入公式(7)即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

上述本发明遥感影像融合方法中A步骤优选为：

对于客观存在的一个场景L，可以使用不同的传感器对其进行观测并成像；在一定的观测条件之下，利用某个特定的传感器，可以得到对场景L的一个观测图像O；可以用一个线性模型来表示由传感器i得到的观测图像O_i与实际场景L之间的关系：

                     O_i＝B_iL+A_i+N_i  公式(8)

其中，i＝1，2，…，q为传感器的数量，B_i是第i个传感器的增益，A_i是第i个传感器的偏差，N_i是第i个传感器的噪声，是与真实场景L无关的高斯白噪声；在成像过程中，针对某一特定的传感器，其增益参数和偏差参数是相对稳定的，因此，该模型中的参数以及噪声的分布特性都假定在空间位置上是缓慢变化的；

为了应用多分辨塔式融合的方法，需要将空间域的通用成像模型转换到需要的变换域(如小波域、Contourlet域等)；由于Contourlet变换为线性变换，因此，在空间域的成像模型经变换到Contourlet域之后，变换系数之间的这种线性关系保持不变，即：

             o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置；n_i，k表示经Contourlet变换后的噪声，依然是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；为描述方便，将公式(1)写成向量的形式：

                   o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目，也可认为是待融合的波段数；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

假定噪声是与真实场景L无关的高斯白噪声，因此，经过Contourlet变换后，可以认为噪声依然是一种高斯白噪声，即：

$n_k~N(0,{\mathrm{\Σ}}_{n_k})$ 公式(9)

其中， ${\mathrm{\Σ}}_{n_k}=\mathrm{diag}({\mathrm{\σ}}_{n_{1,k}}^2,{\mathrm{\σ}}_{n_{2,k}}^2,...,{\mathrm{\σ}}_{n_{q,k}}^2).$

上述本发明遥感影像融合方法中B步骤优选为：

根据成像模型，利用最大后验概率方法求解不同传感器影像融合后的影像的问题可以表述为在已知观测值o_k的情况下，求真实场景值l_k的最大后验概率

根据贝叶斯公式，在o_k的条件下l_k的后验概率可以表示为如下形式：

$p\left(l_k\right|o_k)=\frac{p\left(o_k\right|l_k)p\left(l_k\right)}{p\left(o_k\right)}$ 公式(10)

然后，通过求解真实场景经Contourlet变换后l_k的概率密度p(l_k)、传感器图像经Contourlet变换后的条件概率密度p(o_k|l_k)和传感器图像经Contourlet变换后的概率密度p(o_k)，即可求得在o_k的条件下l_k的后验概率p(l_k|o_k)；

用高斯混合模型来描述经Contourlet变换后l_k的概率密度p(l_k)，即：

$p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)f\left(l_k\right|m)$ 公式(11)

其中，p_k(m)是在位置k处的系数是状态m的概率，M代表状态数；f(l_k|m)是真实场景经Contourlet变换后的条件概率，可以使用高斯分布来进行描述：

$f\left(l_k\right|m)=\frac{1}{\sqrt{2{\mathrm{\π\σ}}_{k,m}^2}}\exp (-\frac{{(l_k-{\mathrm{\μ}}_{k,m})}^2}{2{\mathrm{\σ}}_{k,m}^2})=g(l_k;u_{k,m},{\mathrm{\σ}}_{k,m}^2)$ 公式(12)

u_k，m代表在状态m时的均值，σ_k，m ²代表状态m时的方差；那么，可以得到条件概率密度：

$p\left(o_k\right|l_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|{\mathrm{\Σ}}_{n_k}\right|}^{1/2}}\exp [-\frac{1}{2}{(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})]$ 公式(13)

由此可得边缘概率密度p(o_k)为：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(14)

其中，

U_k，m＝E(o_k)＝β_k，mμ_k，m+α_k，m       公式(15)

$C_{k,m}=E\left[(o_k-U_{k,m}){(o_k-U_{k,m})}^T\right]={\mathrm{\σ}}_{k,m}^2{\mathrm{\β}}_{k,m}{\mathrm{\β}}_{k,m}^T+{\mathrm{\Σ}}_{n_k}$ 公式(16)

最后，将公式(11)、(13)和(14)代入公式(10)整理后，可以得到后验概率p(l_k|o_k)：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2{\mathrm{\π\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\mathrm{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)。

上述本发明遥感影像融合方法中C步骤优选为：

对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来要对α_k，m和β_k，m进行估计；在这里，采用最小二乘方法；该方法假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、

和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\text{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将和公式(17)和(18)代入公式(7)，即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

本发明通过结合Contourlet变换和隐马尔可夫树模型，在多尺度通用传感器模型的框架下，对融合后的Contourlet系数进行最大后验概率估计，建立了一种基于多尺度通用传感器模型的遥感影像融合方法。本发明遥感影像融合方法首次将CHMT模型与多尺度通用传感器模型相结合用于融合研究，并成功地将待融合图像进行主成份变换后的第一主成份用于真实场景的先验概率估计，结果表明，本发明遥感影像融合方法在提高融合影像空间细节信息的同时，可以有效地保持多光谱影像的光谱信息，进一步推动了遥感影像融合理论和技术的发展。

附图说明

图1是本发明所述方法与常用方法的融合效果比较图；

图2是本发明所述方法的流程图。

其中，图1(a)为AIRSAR C波段HH极化影像；图1(b)为SPOT4多光谱影像(波段组合为3、2和1)；图1(c)为使用主成分融合法(PCA)的结果；图1(d)为IHS融合法(HIS)的结果；图1(e)为模极大值小波融合法(MWF)的结果；图1(f)为基于多尺度通用传感器模型融合法(CHMT)的结果。

下述实施例用于进一步说明本发明但不限于本发明。

实验例1本发明方法与其他方法的融合效果比较实验

将本发明方法(实施例1)与其他方法(PCA方法、IHS方法、MWF方法)进行融合效果比较实验(见图1)：从目视效果比较来看，4种方法在提高融合结果清晰度的基础之上，均能够在一定程度上保持SPOT XS影像的光谱信息。在细节保留方面，PCA方法和MWF方法的效果比较突出，但PCA方法得到的结果的颜色与SPOT XS影像相比有些失真，且失真的程度比其他两种方法严重。IHS方法的结果目视效果不佳，并在图像上出现了一些“黑洞”。CHMT方法能够较好地保持光谱信息，并提高了融合结果的清晰度。

性能分析：对于遥感图像融合效果的评价，应综合考虑空间细节信息的增强与光谱信息的保持。所以，一般应综合利用两类统计参数来进行分析与评价：一类反映空间细节信息，如信息熵(Entropy)和清晰度(Average Gradient)；另一类反映光谱信息，如扭曲程度(Spectral Distortion)与相关系数(Correlation Coefficient)。

表1列出了AIRSAR C-HH影像与SPOT XS多光谱影像融合结果的统计结果。总体来看，与目视比较的结果相一致，各个方法均提高了融合结果的清晰度。均值方面，各种方法的差别不大；各种方法的熵和平均梯度(A.G.)均有较大提高，表明各种方法都提高了融合结果的清晰度，其中MWF方法的平均梯度提高最为显著。在光谱信息保持方面，通过观察光谱扭曲(S.D.)和相关系数(C.C.)发现，PCA方法的第1波段与SPOT XS影像的第1波段保持的较好，相关系数达到0.9，但第2、3波段却均为负相关，表明PCA方法融合的结果严重偏离了多光谱影像，这一点与目视比较的结果一致。IHS方法不论是清晰度方面还是光谱保持方面，均不如MWF和CHMT方法。比较MWF方法和CHMT方法，结果发现，MWF方法在提高融合结果的清晰度方面优于CHMT方法(MWF的平均梯度高于CHMT方法的平均梯度)，而CHMT方法的光谱扭曲和相关系数均好于MWF方法，表明CHMT方法保持光谱信息的能力优于MWF方法。

           表1 融合结果的性能评价指标比较

下述实施例均能实现上述实验例的效果。

具体实施方式

实施例1：

对SAR与可见光影像进行融合，采用AIRSAR C波段HH极化影像作为SAR影像，其距离向分辨率约为6米，方位向分辨率约为9米，如图1(a)；采用SPOT4 XS影像的3、2、1波段组合作为多光谱影像，其空间分辨率为20米，如图1(b)；

I、将AIRSAR影像重采样为9米×9米；将SPOT4 XS影像与AIRSAR影像配准后，也重采样为9米×9米；

II、将AIRSAR影像与SPOT4 XS影像的4个波段的数据应用主成份变换，获取第一主成分分量PC₁；

III、对第一主成份PC₁利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)和u_k，m；

IV、在Contourlet分解的某尺度方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)

对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(1)，求得融合后的Contourlet系数表达式；

V、接下来要采用最小二乘法对α_k，m和β_k，m进行估计，具体做法是：假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\text{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将公式(17)和(18)代入公式(7)，即可求得融合后的Contourlet系数值；

VI、对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例2：本发明遥感影像融合方法

A、建立Contourlet域通用传感器模型，其可以表示为：

           o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，i＝1，2，…，q为传感器的数量，也可认为是待融合的波段数，k代表位置，即经变换后的Contourlet系数在某一尺度其对应的某一方向子波段中的位置，l_k是实际场景L经Contourlet变换后在k处的值，o_i，k是由传感器i观测的图像O_i经Contourlet变换后在k处的值，β_i，k和α_i，k分别是经Contourlet变换后第i个传感器在k处的增益和偏差，n_i，k表示经Contourlet变换后的噪声，是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；将公式(1)写成向量的形式为：

                o_k＝β_kl_k+α_k+n_k  公式(2)

各字母代表的含义如公式(3)所示；其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

B、在通用传感器模型的理论下，根据公式(2)的成像模型，已知观测值o_k的情况下，真实场景l_k的后验概率可以表示为：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量PC₁，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来采用最小二乘方法对α_k，m和β_k，m进行估计，分别得到β_k，m和α_k，m的估计值

和

最后，将

和

代入公式(7)即可求得融合后的Contourlet系数值；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例3：本发明遥感影像融合方法

A、对于客观存在的一个场景L，使用不同的传感器对其进行观测并成像；在一定的观测条件之下，利用某个特定的传感器，可以得到对场景L的一个观测图像O；可以用一个线性模型来表示由传感器i得到的观测图像O_i与实际场景L之间的关系：

                   O_i＝B_iL+A_i+N_i  公式(8)

其中，i＝1，2，…，q为传感器的数量，B_i是第i个传感器的增益，A_i是第i个传感器的偏差，N_i是第i个传感器的噪声，是与真实场景L无关的高斯白噪声；在成像过程中，针对某一特定的传感器，其增益参数和偏差参数是相对稳定的，因此，该模型中的参数以及噪声的分布特性都假定在空间位置上是缓慢变化的；

为了应用多分辨塔式融合的方法，需要将空间域的通用成像模型转换到需要的变换域(如小波域、Contourlet域等)；由于Contourlet变换为线性变换，因此，在空间域的成像模型经变换到Contourlet域之后，变换系数之间的这种线性关系保持不变，即：

            o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置；n_i，k表示经Contourlet变换后的噪声，依然是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；为描述方便，将公式(1)写成向量的形式：

            o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目，也可认为是待融合的波段数；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

假定噪声是与真实场景L无关的高斯白噪声，因此，经过Contourlet变换后，可以认为噪声依然是一种高斯白噪声，即：

$n_k~N(0,{\mathrm{\Σ}}_{n_k})$ 公式(9)

其中， ${\mathrm{\Σ}}_{n_k}=\mathrm{diag}({\mathrm{\σ}}_{n_{1,k}}^2,{\mathrm{\σ}}_{n_{2,k}}^2,...,{\mathrm{\σ}}_{n_{q,k}}^2);$

B、在通用传感器模型的理论下，根据公式(2)的成像模型，已知观测值o_k的情况下，真实场景l_k的后验概率可以表示为：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量PC₁，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来采用最小二乘方法对α_k，m和β_k，m进行估计，分别得到β_k，m和α_k，m的估计值和

最后，将

和

代入公式(7)即可求得融合后的Contourlet系数值；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例4：本发明遥感影像融合方法

A、建立Contourlet域通用传感器模型，其可以表示为：

           o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，i＝1，2，…，q为传感器的数量，也可认为是待融合的波段数，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置，l_k是实际场景L经Contourlet变换后在k处的值，o_i，k是由传感器i观测的图像o_i经Contourlet变换后在k处的值，β_i，k和α_i，k分别是经Contourlet变换后第i个传感器在k处的增益和偏差，n_i，k表示经Contourlet变换后的噪声，是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；将公式(1)式写成向量的形式为：

                o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

B、根据成像模型，利用最大后验概率方法求解不同传感器影像融合后的影像的问题可以表述为在已知观测值o_k的情况下，求真实场景值l_k的最大后验概率

根据贝叶斯公式，在o_k的条件下l_k的后验概率可以表示为如下形式：

$p\left(l_k\right|o_k)=\frac{p\left(o_k\right|l_k)p\left(l_k\right)}{p\left(o_k\right)}$ 公式(10)

然后，通过求解真实场景经Contourlet变换后l_k的概率密度p(l_k)、传感器图像经Contourlet变换后的条件概率密度p(o_k|l_k)和传感器图像经Contourlet变换后的概率密度p(o_k)，即可求得在o_k的条件下l_k的后验概率p(l_k|o_k)；

用高斯混合模型来描述经Contourlet变换后l_k的概率密度p(l_k)，即：

$p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)f\left(l_k\right|m)$ 公式(11)

其中，p_k(m)是在位置k处的系数是状态m的概率，M代表状态数；f(l_k|m)是真实场景经Contourlet变换后的条件概率，可以使用高斯分布来进行描述：

$f\left(l_k\right|m)=\frac{1}{\sqrt{2{\mathrm{\π\σ}}_{k,m}^2}}\exp (-\frac{{(l_k-{\mathrm{\μ}}_{k,m})}^2}{2{\mathrm{\σ}}_{k,m}^2})=g(l_k;u_{k,m},{\mathrm{\σ}}_{k,m}^2)$ 公式(12)

u_k，m代表在状态m时的均值，σ_k，m ²代表状态m时的方差；那么，可以得到条件概率密度：

$p\left(o_k\right|l_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|{\mathrm{\Σ}}_{n_k}\right|}^{1/2}}\exp [-\frac{1}{2}{(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})]$ 公式(13)

由此可得边缘概率密度p(o_k)为：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(14)

其中，

U_k，m＝E(o_k)＝β_k，mμ_k，m+α_k，m      公式(15)

$C_{k,m}=E\left[(o_k-U_{k,m}){(o_k-U_{k,m})}^T\right]={\mathrm{\σ}}_{k,m}^2{\mathrm{\β}}_{k,m}{\mathrm{\β}}_{k,m}^T+{\mathrm{\Σ}}_{n_k}$ 公式(16)

最后，将公式(11)、(13)和(14)代入公式(10)整理后，可以得到后验概率p(l_k|o_k)：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来采用最小二乘方法要对α_k，m和β_k，m进行估计；分别得到β_k，m和α_k，m的估计值

和

最后，将

和

代入公式(7)即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例5：本发明遥感影像融合方法

A、建立Contourlet域通用传感器模型，其可以表示为：

          o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，i＝1，2，…，q为传感器的数量，也可认为是待融合的波段数，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置，l_k是实际场景L经Contourlet变换后在k处的值，o_i，k是由传感器i观测的图像o_i经Contourlet变换后在k处的值，β_i，k和α_i，k分别是经Contourlet变换后第i个传感器在k处的增益和偏差，n_i，k表示经Contourlet变换后的噪声，是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；将公式(1)式写成向量的形式为：

                    o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

B、在通用传感器模型的理论下，根据公式(2)的成像模型，已知观测值O_k的情况下，真实场景l_k的后验概率可以表示为：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来要对α_k，m和β_k，m进行估计；在这里，采用最小二乘方法；该方法假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、

和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\text{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将和公式(17)和(18)代入公式(7)，即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例6：本发明遥感影像融合方法

A、对于客观存在的一个场景L，使用不同的传感器对其进行观测并成像；在一定的观测条件之下，利用某个特定的传感器，可以得到对场景L的一个观测图像O；可以用一个线性模型来表示由传感器i得到的观测图像O_i与实际场景L之间的关系：

                O_i＝B_iL+A_i+N_i  公式(8)

其中，i＝1，2，…，q为传感器的数量，B_i是第i个传感器的增益，A_i是第i个传感器的偏差，N_i是第i个传感器的噪声，是与真实场景L无关的高斯白噪声；在成像过程中，针对某一特定的传感器，其增益参数和偏差参数是相对稳定的，因此，该模型中的参数以及噪声的分布特性都假定在空间位置上是缓慢变化的；

为了应用多分辨塔式融合的方法，需要将空间域的通用成像模型转换到需要的变换域(如小波域、Contourlet域等)；由于Contourlet变换为线性变换，因此，在空间域的成像模型经变换到Contourlet域之后，变换系数之间的这种线性关系保持不变，即：

           o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置；n_i，k表示经Contourlet变换后的噪声，依然是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；为描述方便，将公式(1)写成向量的形式：

              o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目，是待融合的波段数；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

假定噪声是与真实场景L无关的高斯白噪声，因此，经过Contourlet变换后，可以认为噪声依然是一种高斯白噪声，即：

$n_k~N(0,{\mathrm{\Σ}}_{n_k})$ 公式(9)

其中， ${\mathrm{\Σ}}_{n_k}=\mathrm{diag}({\mathrm{\σ}}_{n_{1,k}}^2,{\mathrm{\σ}}_{n_{2,k}}^2,...,{\mathrm{\σ}}_{n_{q,k}}^2);$

B、根据成像模型，利用最大后验概率方法求解不同传感器影像融合后的影像的问题可以表述为在已知观测值o_k的情况下，求真实场景值l_k的最大后验概率

根据贝叶斯公式，在o_k的条件下l_k的后验概率可以表示为如下形式：

$p\left(l_k\right|o_k)=\frac{p\left(o_k\right|l_k)p\left(l_k\right)}{p\left(o_k\right)}$ 公式(10)

然后，通过求解真实场景经Contourlet变换后l_k的概率密度p(l_k)、传感器图像经Contourlet变换后的条件概率密度p(o_k|l_k)和传感器图像经Contourlet变换后的概率密度p(o_k)，即可求得在o_k的条件下l_k的后验概率p(l_k|o_k)；

用高斯混合模型来描述经Contourlet变换后l_k的概率密度p(l_k)，即：

$p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)f\left(l_k\right|m)$ 公式(11)

其中，p_k(m)是在位置k处的系数是状态m的概率，M代表状态数；f(l_k|m)是真实场景经Contourlet变换后的条件概率，可以使用高斯分布来进行描述：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(12)

u_k，m代表在状态m时的均值，σ_k，m ²代表状态m时的方差；那么，可以得到条件概率密度：

$p\left(o_k\right|l_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|{\mathrm{\Σ}}_{n_k}\right|}^{1/2}}\exp [-\frac{1}{2}{(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})]$ 公式(13)

由此可得边缘概率密度p(o_k)为：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(14)

其中，

U_k，m＝E(o_k)＝β_k，mμ_k，m+α_k，m      公式(15)

$C_{k,m}=E\left[(o_k-U_{k,m}){(o_k-U_{k,m})}^T\right]={\mathrm{\σ}}_{k,m}^2{\mathrm{\β}}_{k,m}{\mathrm{\β}}_{k,m}^T+{\mathrm{\Σ}}_{n_k}$ 公式(16)

最后，将公式(11)、(13)和(14)代入公式(10)整理后，可以得到后验概率p(l_k|o_k)：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)式，求得融合后的Contourlet系数表达式；接下来采用最小二乘方法要对α_k，m和β_k，m进行估计；分别得到β_k，m和α_k，m的估计值

和

最后，将

和代入公式(7)即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例7：本发明遥感影像融合方法

A、对于客观存在的一个场景L，使用不同的传感器对其进行观测并成像；在一定的观测条件之下，利用某个特定的传感器，可以得到对场景L的一个观测图像O；可以用一个线性模型来表示由传感器i得到的观测图像O_i与实际场景L之间的关系：

                 O_i＝B_iL+A_i+N_i  公式(8)

其中，i＝1，2，…，q为传感器的数量，B_i是第i个传感器的增益，A_i是第i个传感器的偏差，N_i是第i个传感器的噪声，是与真实场景L无关的高斯白噪声；在成像过程中，针对某一特定的传感器，其增益参数和偏差参数是相对稳定的，因此，该模型中的参数以及噪声的分布特性都假定在空间位置上是缓慢变化的；

为了应用多分辨塔式融合的方法，需要将空间域的通用成像模型转换到需要的变换域(如小波域、Contourlet域等)；由于Contourlet变换为线性变换，因此，在空间域的成像模型经变换到Contourlet域之后，变换系数之间的这种线性关系保持不变，即：

            o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置；n_i，k表示经Contourlet变换后的噪声，依然是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；为描述方便，将公式(1)写成向量的形式：

            o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目，是待融合的波段数；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

假定噪声是与真实场景L无关的高斯白噪声，因此，经过Contourlet变换后，可以认为噪声依然是一种高斯白噪声，即：

$n_k~N(0,{\mathrm{\Σ}}_{n_k})$ 公式(9)

其中， ${\mathrm{\Σ}}_{n_k}=\mathrm{diag}({\mathrm{\σ}}_{n_{1,k}}^2,{\mathrm{\σ}}_{n_{2,k}}^2,...,{\mathrm{\σ}}_{n_{q,k}}^2);$

B、在通用传感器模型的理论下，根据公式(2)的成像模型，已知观测值O_k的情况下，真实场景l_k的后验概率可以表示为：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来要对α_k，m和β_k，m进行估计；在这里，采用最小二乘方法；该方法假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、

和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\text{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将和公式(17)和(18)代入公式(7)，即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例8：本发明遥感影像融合方法

A、建立Contourlet域通用传感器模型，其可以表示为：

          o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，i＝1，2，…，q为传感器的数量，是待融合的波段数，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置，l_k是实际场景L经Contourlet变换后在k处的值，o_i，k是由传感器i观测的图像o_i经Contourlet变换后在k处的值，β_i，k和α_i，k分别是经Contourlet变换后第i个传感器在k处的增益和偏差，n_i，k表示经Contourlet变换后的噪声，是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；将公式(1)式写成向量的形式为：

                  o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

B、根据成像模型，利用最大后验概率方法求解不同传感器影像融合后的影像的问题可以表述为在已知观测值o_k的情况下，求真实场景值l_k的最大后验概率

根据贝叶斯公式，在o_k的条件下l_k的后验概率可以表示为如下形式：

$p\left(l_k\right|o_k)=\frac{p\left(o_k\right|l_k)p\left(l_k\right)}{p\left(o_k\right)}$ 公式(10)

然后，通过求解真实场景经Contourlet变换后l_k的概率密度p(l_k)、传感器图像经Contourlet变换后的条件概率密度p(o_k|l_k)和传感器图像经Contourlet变换后的概率密度p(o_k)，即可求得在o_k的条件下l_k的后验概率p(l_k|o_k)；

用高斯混合模型来描述经Contourlet变换后l_k的概率密度p(l_k)，即：

$p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)f\left(l_k\right|m)$ 公式(11)

其中，p_k(m)是在位置k处的系数是状态m的概率，M代表状态数；f(l_k|m)是真实场景经Contourlet变换后的条件概率，可以使用高斯分布来进行描述：

$f\left(l_k\right|m)=\frac{1}{\sqrt{2{\mathrm{\π\σ}}_{k,m}^2}}\exp (-\frac{{(l_k-{\mathrm{\μ}}_{k,m})}^2}{2{\mathrm{\σ}}_{k,m}^2})=g(l_k;u_{k,m},{\mathrm{\σ}}_{k,m}^2)$ 公式(12)

u_k，m代表在状态m时的均值，σ_k，m ²代表状态m时的方差；那么，可以得到条件概率密度：

$p\left(o_k\right|l_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|{\mathrm{\Σ}}_{n_k}\right|}^{1/2}}\exp [-\frac{1}{2}{(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})]$ 公式(13)

由此可得边缘概率密度p(o_k)为：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(14)

其中，

U_k，m＝E(o_k)＝β_k，mμ_k，m+α_k，m    公式(15)

$C_{k,m}=E\left[(o_k-U_{k,m}){(o_k-U_{k,m})}^T\right]={\mathrm{\σ}}_{k,m}^2{\mathrm{\β}}_{k,m}{\mathrm{\β}}_{k,m}^T+{\mathrm{\Σ}}_{n_k}$ 公式(16)

最后，将公式(11)、(13)和(14)代入公式(10)整理后，可以得到后验概率p(l_k|o_k)：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来要对α_k，m和β_k，m进行估计；在这里，采用最小二乘方法；该方法假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\text{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将和公式(17)和(18)代入公式(7)，即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

实施例9：本发明遥感影像融合方法

A、对于客观存在的一个场景L，使用不同的传感器对其进行观测并成像；在一定的观测条件之下，利用某个特定的传感器，可以得到对场景L的一个观测图像O；可以用一个线性模型来表示由传感器i得到的观测图像O_i与实际场景L之间的关系：

                    O_i＝B_iL+A_i+N_i  公式(8)

其中，i＝1，2，…，q为传感器的数量，B_i是第i个传感器的增益，A_i是第i个传感器的偏差，N_i是第i个传感器的噪声，是与真实场景L无关的高斯白噪声；在成像过程中，针对某一特定的传感器，其增益参数和偏差参数是相对稳定的，因此，该模型中的参数以及噪声的分布特性都假定在空间位置上是缓慢变化的；

为了应用多分辨塔式融合的方法，需要将空间域的通用成像模型转换到需要的变换域(如小波域、Contourlet域等)；由于Contourlet变换为线性变换，因此，在空间域的成像模型经变换到Contourlet域之后，变换系数之间的这种线性关系保持不变，即：

           o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置；n_i，k表示经Contourlet变换后的噪声，依然是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；为描述方便，将公式(1)写成向量的形式：

                o_k＝β_kl_k+α_k+n_k  公式(2)

其中，l_k＝[l_1，k l_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目，是待融合的波段数；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right]$ ${\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right]$ ${\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right]$ $n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

假定噪声是与真实场景L无关的高斯白噪声，因此，经过Contourlet变换后，可以认为噪声依然是一种高斯白噪声，即：

$n_k~N(0,{\mathrm{\Σ}}_{n_k})$ 公式(9)

其中， ${\mathrm{\Σ}}_{n_k}=\mathrm{diag}({\mathrm{\σ}}_{n_{1,k}}^2,{\mathrm{\σ}}_{n_{2,k}}^2,...,{\mathrm{\σ}}_{n_{q,k}}^2);$

B、根据成像模型，利用最大后验概率方法求解不同传感器影像融合后的影像的问题可以表述为在已知观测值o_k的情况下，求真实场景值l_k的最大后验概率

根据贝叶斯公式，在o_k的条件下l_k的后验概率可以表示为如下形式：

$p\left(l_k\right|o_k)=\frac{p\left(o_k\right|l_k)p\left(l_k\right)}{p\left(o_k\right)}$ 公式(10)

然后，通过求解真实场景经Contourlet变换后l_k的概率密度p(l_k)、传感器图像经Contourlet变换后的条件概率密度p(o_k|l_k)和传感器图像经Contourlet变换后的概率密度p(o_k)，即可求得在o_k的条件下l_k的后验概率p(l_k|o_k)；

用高斯混合模型来描述经Contourlet变换后l_k的概率密度p(l_k)，即：

$p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)f\left(l_k\right|m)$ 公式(11)

其中，p_k(m)是在位置k处的系数是状态m的概率，M代表状态数；f(l_k|m)是真实场景经Contourlet变换后的条件概率，可以使用高斯分布来进行描述：

$f\left(l_k\right|m)=\frac{1}{\sqrt{2{\mathrm{\π\σ}}_{k,m}^2}}\exp (-\frac{{(l_k-{\mathrm{\μ}}_{k,m})}^2}{2{\mathrm{\σ}}_{k,m}^2})=g(l_k;u_{k,m},{\mathrm{\σ}}_{k,m}^2)$ 公式(12)

u_k，m代表在状态m时的均值，σ_k，m ²代表状态m时的方差；那么，可以得到条件概率密度：

$p\left(o_k\right|l_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|{\mathrm{\Σ}}_{n_k}\right|}^{1/2}}\exp [-\frac{1}{2}{(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})]$ 公式(13)

由此可得边缘概率密度p(o_k)为：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(14)

其中，

U_k，m＝E(o_k)＝β_k，mμ_k，m+α_k，m    公式(15)

$C_{k,m}=E\left[(o_k-U_{k,m}){(o_k-U_{k,m})}^T\right]={\mathrm{\σ}}_{k,m}^2{\mathrm{\β}}_{k,m}{\mathrm{\β}}_{k,m}^T+{\mathrm{\Σ}}_{n_k}$ 公式(16)

最后，将公式(11)、(13)和(14)代入公式(10)整理后，可以得到后验概率p(l_k|o_k)：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\text{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景l的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，可得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\text{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来要对α_k，m和β_k，m进行估计；在这里，采用最小二乘方法；该方法假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、

和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\mathrm{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将和公式(17)和(18)代入公式(7)，即可求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即可得融合后的影像。

Claims (9)

1、一种遥感影像融合方法，其特征在于基于多尺度通用传感器模型与CHMT模型相结合获得融合后的影像。

2、如权利要求1所述的一种遥感影像融合方法，其特征在于该方法包括如下步骤：

建立Contourlet域多尺度通用传感器模型，计算得到由传感器观测的图像经Contourlet变换后在相应尺度某方向子波段上对应位置处的观测值；对传感器的观测值应用主成分变换，获取第一主成分分量，作为对实际场景的先验估计，然后对其利用CHMT模型进行训练，获得条件概率，最后，计算得在Contourlet分解的某尺度某方向波段上实际场景的最大后验估计值；采用最小二乘法求得待融合影像的传感器系数，然后根据实际场景的最大后验估计值，求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即得融合后的影像。

3、如权利要求1或2所述的一种遥感影像融合方法，其特征在于该方法包括如下步骤：

A、建立Contourlet域通用传感器模型，其表示为：

o_i，k＝β_i，kl_k+α_i，k+n_i，k公式(1)

其中，i＝1，2，…，q为传感器的数量或待融合的波段数，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置，l_k是实际场景L经Contourlet变换后在k处的值，o_i，k是由传感器i观测的图像o_i经Contourlet变换后在k处的值，β_i，k和α_i，k分别是经Contourlet变换后第i个传感器在k处的增益和偏差，n_i，k表示经Contourlet变换后的噪声，是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；将公式(1)式写成向量的形式为：

o_k＝β_kl_k+α_k+n_k公式(2)

其中，l_k＝[l_1，kl_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right],{\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right],{\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right],n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

B、在通用传感器模型的理论下，根据公式(2)的成像模型，已知观测值Q_k的情况下，真实场景l_k的后验概率表示为：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{{2\mathrm{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景1的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)；

C、对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)式，求得融合后的Contourlet系数表达式；接下来采用最小二乘方法要对α_k，m和β_k，m进行估计；分别得到β_k，m和α_k，m的估计值

和

最后，将和

代入公式(7)求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即得融合后的影像。

4、如权利要求3所述的一种遥感影像融合方法，其特征在于其中所述的步骤A为：

对于客观存在的一个场景L，使用不同的传感器对其进行观测并成像；在一定的观测条件之下，利用某个特定的传感器，得到对场景L的一个观测图像O；用一个线性模型来表示由传感器i得到的观测图像O_i与实际场景L之间的关系：

O_i＝B_iL+A_i+N_i公式(8)

其中，i＝1，2，…，q为传感器的数量，B_i是第i个传感器的增益，A_i是第i个传感器的偏差，N_i是第i个传感器的噪声，是与真实场景L无关的高斯白噪声；在成像过程中，针对某一特定的传感器，其增益参数和偏差参数是相对稳定的，因此，该模型中的参数以及噪声的分布特性都假定在空间位置上是缓慢变化的；

为了应用多分辨塔式融合的方法，需要将空间域的通用成像模型转换到需要的变换域(如小波域、Contourlet域等)；由于Contourlet变换为线性变换，因此，在空间域的成像模型经变换到Contourlet域之后，变换系数之间的这种线性关系保持不变，即：

o_i，k＝β_i，kl_k+α_i，k+n_i，k  公式(1)

其中，k代表位置，即经变换后的系数在某一尺度其对应的某一方向子波段中的位置；n_i，k表示经Contourlet变换后的噪声，依然是高斯白噪声，即n_i，k～N(0，σ_i，k ²)；为描述方便，将公式(2)写成向量的形式：

o_k＝β_kl_k+α_k+n_k公式(2)

其中，l_k＝[l_1，kl_2，k…l_q，k]，代表场景L经Contourlet分解后在相应层中某方向子波段上对应位置的系数，q为传感器的数目，也可认为是待融合的波段数；

$o_k=\left[\begin{array}{c}{o}_{1,k}\\ o_{1,k}\\ ...\\ o_{q,k}\end{array}\right],{\mathrm{\β}}_k=\left[\begin{array}{c}{\mathrm{\β}}_{1,k}\\ {\mathrm{\β}}_{2,k}\\ ...\\ {\mathrm{\β}}_{q,k}\end{array}\right],{\mathrm{\α}}_k=\left[\begin{array}{c}{\mathrm{\α}}_{1,k}\\ {\mathrm{\α}}_{2,k}\\ ...\\ {\mathrm{\α}}_{q,k}\end{array}\right],n_k=\left[\begin{array}{c}{n}_{1,k}\\ n_{2,k}\\ ...\\ n_{q,k}\end{array}\right]$ 公式(3)

假定噪声是与真实场景L无关的高斯白噪声，因此，经过Contourlet变换后，噪声依然是一种高斯白噪声，即：

$n_k~N(0,{\mathrm{\Σ}}_{n_k})$ 公式(9)

其中， ${\mathrm{\Σ}}_{n_k}=\mathrm{diag}({\mathrm{\σ}}_{n_{1,k}}^2,{\mathrm{\σ}}_{n_{2,k}}^2,...,{\mathrm{\σ}}_{n_{q,k}}^2).$

5、如权利要求3所述的一种遥感影像融合方法，其特征在于其中所述的步骤B为：

根据成像模型，利用最大后验概率方法求解不同传感器影像融合后的影像的问题表述为在已知观测值o_k的情况下，求真实场景值l_k的最大后验概率

根据贝叶斯公式，在o_k的条件下l_k的后验概率表示为如下形式：

$p\left(l_k\right|o_k)=\frac{p\left(o_k\right|l_k)p\left(l_k\right)}{p\left(o_k\right)}$ 公式(10)

然后，通过求解真实场景经Contourlet变换后l_k的概率密度p(l_k)、传感器图像经Contourlet变换后的条件概率密度p(o_k|l_k)和传感器图像经Contourlet变换后的概率密度p(o_k)，即可求得在o_k的条件下l_k的后验概率p(l_k|o_k)；

用高斯混合模型来描述经Contourlet变换后l_k的概率密度p(l_k)，即：

$p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)f\left(l_k\right|m)$ 公式(11)

其中，p_k(m)是在位置k处的系数是状态m的概率，M代表状态数；f(l_k|m)是真实场景经Contourlet变换后的条件概率，使用高斯分布来进行描述：

$f\left(l_k\right|m)=\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{k,m}^2}}\exp (-\frac{{(l_k-{\mathrm{\μ}}_{k,m})}^2}{{2\mathrm{\σ}}_{k,m}^2})=g(l_k;u_{k,m},{\mathrm{\σ}}_{k,m}^2)$ 公式(12)

u_k，m代表在状态m时的均值，σ_k，m ²代表状态m时的方差；得到条件概率密度：

$p\left(o_k\right|l_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|{\mathrm{\Σ}}_{n_k}\right|}^{1/2}}\exp [-\frac{1}{2}{(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})]$ 公式(13)

由此得边缘概率密度p(o_k)为：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(14)

其中，

U_k，m＝E(o_k)＝β_k，mμ_k，m+α_k，m公式(15)

$C_{k,m}=E\left[(o_k-U_{k,m}){(o_k-U_{k,m})}^T\right]={\mathrm{\σ}}_{k,m}^2{\mathrm{\β}}_{k,m}{\mathrm{\β}}_{k,m}^T+{\mathrm{\Σ}}_{n_k}$ 公式(16)

最后，将公式(11)、(13)和(14)代入公式(10)整理后，得到后验概率p(l_k|o_k)：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\mathrm{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景1的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²最后，得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)。

6、如权利要求4所述的一种遥感影像融合方法，其特征在于其中所述的步骤B为：

根据成像模型，利用最大后验概率方法求解不同传感器影像融合后的影像的问题表述为在已知观测值o_k的情况下，求真实场景值l_k的最大后验概率

根据贝叶斯公式，在o_k的条件下l_k的后验概率表示为如下形式：

$p\left(l_k\right|o_k)=\frac{p\left(o_k\right|l_k)p\left(l_k\right)}{p\left(o_k\right)}$ 公式(10)

求真实场景经Contourlet变换后l_k的概率密度p(l_k)、传感器图像经Contourlet变换后的条件概率密度p(o_k|l_k)和传感器图像经Contourlet变换后的概率密度p(o_k)；

用高斯混合模型来描述经Contourlet变换后l_k的概率密度p(l_k)，即：

$p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)f\left(l_k\right|m)$ 公式(11)

其中，p_k(m)是在位置k处的系数是状态m的概率，M代表状态数；f(l_k|m)是真实场景经Contourlet变换后的条件概率，使用高斯分布来进行描述：

$f\left(l_k\right|m)=\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{k,m}^2}}\exp (-\frac{{(l_k-{\mathrm{\μ}}_{k,m})}^2}{{2\mathrm{\σ}}_{k,m}^2})=g(l_k;u_{k,m},{\mathrm{\σ}}_{k,m}^2)$ 公式(12)

u_k，m代表在状态m时的均值，σ_k，m ²代表状态m时的方差；得到条件概率密度：

$p\left(o_k\right|l_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|{\mathrm{\Σ}}_{n_k}\right|}^{1/2}}\exp [-\frac{1}{2}{(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\β}}_{k,m}{l}_{k,m}-{\mathrm{\α}}_{k,m})]$ 公式(13)

由此得边缘概率密度p(o_k)为：

$p\left(o_k\right)=\underset{k}{\mathrm{\Σ}}p\left(o_k\right|l_k)p\left(l_k\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{{\left(2\mathrm{\π}\right)}^{q/2}{\left|C_{k,m}\right|}^{1/2}}\exp [-\frac{{(o_k-U_{k,m})}^T{C}_{k,m}^{-1}(o_k-U_{k,m})}{2}]$ 公式(14)

其中，

U_k，m＝E(o_k)＝β_k，mμ_k，m+α_k，m公式(15)

$C_{k,m}=E\left[(o_k-U_{k,m}){(o_k-U_{k,m})}^T\right]={\mathrm{\σ}}_{k,m}^2{\mathrm{\β}}_{k,m}{\mathrm{\β}}_{k,m}^T+{\mathrm{\Σ}}_{n_k}$ 公式(16)

最后，将公式(11)、(13)和(14)代入公式(10)整理后，得到后验概率p(l_k|o_k)：

$p\left(l_k\right|o_k)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{p}_k\left(m\right)\frac{1}{\sqrt{2\mathrm{\π}{\mathrm{\σ}}_{l|o,m}^2}}\exp [-\frac{{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2}{2{\mathrm{\σ}}_{l|o,m}^2}]$ 公式(4)

其中，后验均值和方差分别为：

${\mathrm{\μ}}_{l|o,m}=E\left[l_k\right|o_k,m]={\mathrm{\σ}}_{l|a,m}^{-2}\·[{\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2]$ 公式(5)

${\mathrm{\σ}}_{l|a,m}^2=\mathrm{cov}\left(l_{k,m}\right|o_k)=E\left[{(l_{k,m}-{\mathrm{\μ}}_{l|o,m})}^2\right|o_k,m]={\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+{\mathrm{\σ}}_{k,m}^{-2}$ 公式(6)

真实场景1的估计值：首先对q个传感器的观测值应用主成分变换，获取第一主成分分量，然后对其利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)、u_k，m和σ_k，m ²，最后，得在Contourlet分解的某尺度某方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)。

7、如权利要求3所述的一种遥感影像融合方法，其特征在于其中所述的步骤C为：

对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来要对α_k，m和β_k，m进行估计；在这里，采用最小二乘方法；该方法假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、

和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\mathrm{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将和公式(17)和(18)代入公式(7)，即求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即得融合后的影像。

8、如权利要求4、5或6任一所述的一种遥感影像融合方法，其特征在于其中所述的步骤C为：

对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；接下来要对α_k，m和β_k，m进行估计；在这里，采用最小二乘方法；该方法假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、

和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\mathrm{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将和公式(17)和(18)代入公式(7)，求得融合后的Contourlet系数；对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，得融合后的影像。

9、一种遥感影像融合方法，其特征在于该方法包括如下步骤：对SAR与可见光影像进行融合，采用AIRSAR C波段HH极化影像作为SAR影像，其距离向分辨率约为6米，方位向分辨率约为9米，；采用SPOT4XS影像的3、2、1波段组合作为多光谱影像，其空间分辨率为20米；

I、将AIRSAR影像重采样为9米×9米；将SPOT4XS影像与AIRSAR影像配准后，也重采样为9米×9米；

II、将AIRSAR影像与SPOT4XS影像的4个波段的数据应用主成份变换，获取第一主成分分量PC₁；

III、对第一主成份PC₁利用CHMT模型进行训练，获得条件概率p(m|l_k，m，θ_o)和u_k，m；

IV、在Contourlet分解的某尺度方向波段上l_k的最大后验估计

为：

${\hat{l}}_{\mathrm{kMAP}}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o){\mathrm{\μ}}_{l|o,m}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}p\left(m\right|l_{k,m},{\mathrm{\θ}}_o)\·\frac{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}(o_k-{\mathrm{\α}}_{k,m})+u_{k,m}/{\mathrm{\σ}}_{k,m}^2)}{({\mathrm{\β}}_{k,m}^T{\mathrm{\Σ}}_{n_k}^{-1}{\mathrm{\β}}_{k,m}+1/{\mathrm{\σ}}_{k,m}^2)}$ 公式(7)

对于经过Contourlet变换得到的待融合图像的各个方向子波段应用公式(7)，求得融合后的Contourlet系数表达式；

V、接下来要采用最小二乘法对α_k，m和β_k，m进行估计，具体做法是：假设在一个较小的区域内R_k内(譬如5×5的窗口)，β_k，m、α_k，m、和σ_k，m ²是常量，通过假设该窗口内的观测值服从高斯分布，并利用窗口内的一阶和二阶统计量，得到β_k，m和α_k，m的最小二乘估计值：

${\widehat{\mathrm{\β}}}_{k,m}=\frac{\sqrt{\mathrm{\λ}}\·u}{{\mathrm{\σ}}_{k,m}}$ 公式(17)

${\widehat{\mathrm{\α}}}_{k,m}={\mathrm{\μ}}_a-{\widehat{\mathrm{\β}}}_{k,m}\·{\mathrm{\μ}}_{k,m}$ 公式(18)

其中，λ是修正噪声协方差矩阵 $P={\mathrm{\Σ}}_o-{\mathrm{\Σ}}_{n_k}$ 的主特征值，u是主特征值对应的主特征向量，∑_o和μ_a分别如下：

${\mathrm{\Σ}}_o=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(o_n-{\mathrm{\μ}}_a){(o_n-{\mathrm{\μ}}_a)}^T$ 公式(19)

${\mathrm{\μ}}_a=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{o}_n$ 公式(20)

式中，N代表小区域R_k内的像素个数；最后，将公式(17)和(18)代入公式(7)，求得融合后的Contourlet系数值；

Ⅵ、对所有尺度所有方向波段进行融合处理之后，应用逆Contourlet变换，即得融合后的影像。

Images