CN102044072A - 基于统计模型的sar图像融合处理方法 - Google Patents

Abstract

本发明公开了一种基于统计模型的SAR图像融合处理方法，属于图像处理领域，首先提取输入图像的边缘特征，以边缘特征为基础，完成局部分析窗口的自适应修正处理，确保局部分析窗口内为同质目标；然后，结合输入图像的统计模型，利用局部分析窗口数据估计输入图像的模型参数；最后，结合输入图像的统计模型，通过最优估计理论完成图像的融合处理。本发明将非线性方程组得最优化估计处理转换为线性方程组的最优化估计处理，简化融合处理过程，从SAR图像和光学图像的统计特性出发，综合考虑SAR图像斑点噪声及图像中的边缘特征对融合处理的影响，有效改善了SAR图像与光学图像及SAR图像之间的融合处理效果。

Description

基于统计模型的SAR图像融合处理方法

技术领域

本发明属于图像处理领域，涉及一种图像融合处理方法，具体涉及一种基于统计模型的SAR图像融合处理方法。

背景技术

随着现代遥感技术的飞速发展，光学、雷达、红外和激光等传感器系统得到了广泛的应用，获取的对地遥感数据(多分辨率、多时相、多光谱)越来越多。由于成像机理和技术条件的限制，任何单一传感器所获取的数据都无法全面反映目标特征，而利用多传感器对目标进行协同探测，获取对目标多视角、多波段、多时相的信息，充分发挥各传感器自身的优势，相互弥补各自的不足，通过精细的融合处理能够得到更为全面、更为准确的目标信息，形成对目标完整一致的信息描述。

多传感器图像融合是指综合两个或多个图像的信息，获得对同一场景或目标更为准确、全面和可靠的图像，使之更加适合于人眼感知或计算机后续处理。遥感图像融合处理的基本流程如图1所示，在开展遥感图像融合处理时，首先对输入图像采取某种变换处理将输入图像转换到某一处理域内，该处理域可以是频域、小波域、甚至是空域。接着，依据融合处理的目的选取图像的某一特性作为融合的参考，进行局部特性统计并根据局部统计的结果确定融合系数。然后，根据提取的融合系数来进行图像合成处理，获取处理域内的合成图像，最后，通过图像逆变换处理得到融合后的图像。

融合处理的核心在于如何确定融合系数来实现图像的合成处理。在传统的处理方法中通常采用局部方差、邻域熵、显著性水平、空间频率等指标来确定融合系数，然而这些指标受噪声的影响很大。考虑到SAR图像中具有强烈的斑点噪声，使得这些常规的指标无法满足SAR图像的融合处理要求。

下面介绍SAR图像与光学图像的统计特性：

1、SAR图像的统计特性；

SAR图像是地面目标平均后向散射系数的反映。对于包含大量点散射体的面目标而言，当点散射体的个数多到可以应用中心极限定理时，认为面目标所对应检波包络电压的同相分量V_I和正交分量V_Q相互独立，并服从高斯分布N(0，σ²)，其中标准差σ由面目标中各点散射体的散射特性决定。相应的可知检波包络电压的功率值

服从指数分布

则：

$f\left(I\right)=\frac{1}{2{\mathrm{\σ}}^2}{e}^{-I/2{\mathrm{\σ}}^2},I\≥0---\left(1\right)$

式中：I表示检波包络电压的功率值；σ为概率分布标准差。

在(1)式中作变量代换：

I＝2σ²·S                (2)

则，

S是严格单调可微函数，S的概率密度函数为：

f(S)＝e^-S        S≥0        (3)

此时，随机变量S满足数学期望和方差均为1的指数分布e(0，1)，即E(S)＝D(S)＝1。

其中，2σ²是反映地物真实RCS的物理量，记R＝2σ²，此时，对于代表每个像素强度的随机变量I，可以将它等效为R与一个呈e(0，1)分布的随机变量S的乘积：

I＝R·S                        (4)

式(4)为斑点噪声乘性模型。其中，R在通常情况下被直接称为目标的真实RCS(忽略接收机热噪声和需要辐射校正的各因素)，S为代表了斑点噪声的随机变量，R与S相互独立。

在SAR图像的实际应用中，为了抑制斑点噪声的影响常常采用多视处理，即对目标的多个独立功率样本进行平均叠加。由于叠加针对于功率值，没有考虑相位信息，也称为非相干叠加。由概率论的知识可知，I_L服从L阶的伽马分布(即自由度为2L的χ²分布)

$f\left(I_L\right)=\frac{L^L}{\mathrm{\Γ}\left(L\right)}\·\frac{1}{R}\·{\left(\frac{I_L}{R}\right)}^{L-1}\·\exp \{-\frac{{\mathrm{LI}}_L}{R}\},I_L\≥0---\left(5\right)$

且：E(I_L)＝R，D(I)＝R²/L。

式中：L表示处理视数；Γ(L)表示L阶伽马分布。

令S_L表示多视功率图像中的斑点噪声，作变量代换I_L＝R·S_L，得到S_L也服从L阶的伽马分布：

$f\left(S_L\right)=\frac{L^L}{\mathrm{\Γ}\left(L\right)}\·{S_L}^{L-1}\·\exp \{-{\mathrm{LS}}_L\},S_L\≥0---\left(6\right)$

且E(S_L)＝1、D(S_L)＝1/L。

可见，在SAR多视功率图像中，也满足斑点噪声的乘性模型。为了表述方便，多视图像的功率和斑点噪声仍然用I和S表示。在多视图像的有关表达式中，令L＝1就得到单视图像的表达式。

2、光学图像的统计特性；

与SAR图像中乘性噪声模型相比，光学图像通常采用加性噪声模型，即光学图像认为是理想场景通过传感器系统传递函数模糊后降采样得到的结果，因此，其数学模型通常表示为：

$a(x,y)=\underset{m=-M/2}{\overset{M/2}{\mathrm{\Σ}}}\underset{n=-N/2}{\overset{N/2}{\mathrm{\Σ}}}h(m,n)R(x-m,y-n)+\mathrm{\ϵ}(x,y)---\left(7\right)$

其中：a(·)表示传感器所获得的实际图像；h(·)表示传感器的系统传递函数；ε(·)图像中的加性噪声；x，y表示图像中各像素的位置。

从上式可以看出，传感器获取的遥感图像是理想场景与传感器系统传递函数卷积处理后叠加噪声的结果。若理想化系统传递函数为冲激响应函数，同时假设引入的噪声为高斯白噪声，光学图像的数学模型可简化为：

a(x，y)＝β(x，y)R(x，y)+ε(x，y)            (8)

其中：β(x，y)表式不同位置的系统增益函数。

发明内容

本发明的目的是为了解决斑点噪声对SAR图像融合处理存在影响的问题，提出一种基于统计模型的SAR图像融合处理方法，从SAR图像和光学图像的统计特性出发来实现SAR图像与光学图像及SAR图像之间的融合处理，综合考虑边缘特征对图像融合处理的影响，根据局部分析窗口内的边缘特征，自适应的修正局部分析窗口，进而提高融合处理的效果。

本发明为一种基于统计模型的SAR图像融合处理方法，首先提取输入图像的边缘特征，以边缘特征为基础，完成局部分析窗口的自适应修正处理，确保局部分析窗口内为同质目标；然后，结合输入图像的统计模型，利用局部分析窗口内的图像数据估计输入图像的模型参数；最后，结合输入图像的统计模型，通过最优化估计处理完成图像的融合处理。从整体上来看，该处理流程主要包含提取输入图像的边缘特征、局部分析窗口自适应修正、统计模型参数及图像合成处理四个部分，包括如下几个步骤：

步骤一：提取输入图像的边缘特征；

采用Canny边缘检测算法进行光学图像边缘轮廓特征提取。该算法首先将输入图像与二维高斯函数进行卷积处理，利用高斯函数的低通特性完成对输入图像的滤波处理，接着对滤波处理后的图像实施微分处理，分别提取水平方向和垂直方向上的处理结果G_y和G_x，得到梯度的大小和方向。

$G=\sqrt{G_x^2+G_y^2}---\left(9\right)$

$\mathrm{\θ}=\mathrm{Arc}\tan \frac{G_y}{G_x}---\left(10\right)$

式中：G和θ分别表示梯度的大小与方向。

对微分处理后的图像进行非极大值抑制处理，以梯度方向为参考，判断该像素点梯度值是否为沿梯度方向上相邻像素的最大梯度值，若该梯度值为最大时，表明该像素点的梯度值为局部极大值，保留该像素的梯度值；若该梯度值不是最大时，表明该像素点的梯度值并非局部极大值，将梯度值设置为0。在完成非极大值抑制处理后，进行双阈值处理。统计微分处理后图像的均值与方差，将均值与两倍方差的和设定为高阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制后图像的高阈值边缘图像；将均值设定为低阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制后图像的低阈值边缘图像。然后以高阈值边缘图像为种子，通过在对应低阈值边缘图像中进行搜索，获得完整的光学图像边缘特征。

采用Touzi边缘检测算法进行边缘轮廓提取。该算法首先利用四个方向上的掩模模版，根据式(11)分别计算水平方向、垂直方向、右倾斜45°和左倾斜45°四个方向上均值的比值γ_0°，γ_45°，γ_90°，γ_135°。

$\mathrm{\γ}=\max \{\frac{{\mathrm{\μ}}_1}{{\mathrm{\μ}}_2},\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}\},(\mathrm{\γ}\≥1)---\left(11\right)$

式中：μ₁表示左侧或上侧区域均值；μ₂表示右侧或下侧区域均值；γ表示分区后两部分较大均值与较小均值的比值。

选择四个方向上最大的比值作为当前分析点的均值比值，并记录最大的比值所对应的方向，作为当前分析点的梯度方向，如公式(12)、(13)所示：

γ_value＝Max(γ_0°，γ_45°，γ_90°，γ_135°)            (12)

γ_direction＝direction(γ_value)                        (13)

式中：γ_value表示当前分析点的均值比值；γ_direction表示当前分析点的梯度方向。

对均值比值处理结果进行非极大值抑制处理，遍历γ_value，根据每一个点的梯度方向γ_direction，对相邻像素均值比值的大小进行比较，以判断当前分析点的均值比值是否为沿梯度方向上的局部极大值，若沿该像素点梯度方向上前后两个像素的均值比值中的任何一个均值比值大于当前像素的均值比值，表明该像素点的均值比值并不是梯度方向上的极大值，将该像素点处的均值比值设置为0；若沿该像素点梯度方向上前后两个像素的均值比值均小于当前像素的均值比值，表明该像素点的均值比值是梯度方向上的极大值，保留该像素点处的均值比值。最后得到像素均值比值的非极大值抑制处理结果。在完成非极大值抑制处理后，对输入图像进行双阈值处理，统计微分处理后图像的均值与方差，将均值与两倍方差的和设定为高阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制后图像的高阈值边缘图像；将均值设定为低阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制后图像的低阈值边缘图像。然后以高阈值边缘图像为种子，通过在对应低阈值边缘图像中进行搜索，获得完整的SAR图像边缘特征。

步骤二：局部分析窗口自适应修正；

在完成边缘特征提取后，利用局部边缘轮廓与当前分析点之间的相互关系来完成对局部分析窗口的自适应修正处理，进而确保局部分析窗口内不存在边缘轮廓特征，整个局部分析窗口内为同质目标。以输入图像的当前分析点为中心，以设定的分析窗口尺寸为半径，截取局部圆形区域作为局部分析窗口。检测局部分析窗口内是否存在边缘轮廓，若各输入图像的对应局部分析窗口内均不存在边缘轮廓时，认为当前分析点处于均匀区域，直接利用原局部分析窗口作为当前分析点的局部分析窗口；若当前分析窗口内存在边缘特征，且当前分析点在所有输入图像中均不位于边缘特征上时，判断当前分析点与边缘轮廓之间的相互位置关系，提取出局部分析窗口内包含当前分析点的边缘轮廓的内接局部均匀区域作为新的局部分析窗口；若当前分析点在部分或全部输入图像中位于边缘轮廓上时，忽略当前分析点不位于边缘轮廓上的输入图像，仅仅考虑当前分析点位于边缘轮廓上的输入图像。根据当前分析点处的局部方差进行加权平均处理，得到当前分析点的融合处理结果；

$\widehat{\mathrm{\α}}(x,y)=\frac{\underset{k}{\mathrm{\Σ}}{\mathrm{\σ}}_k^2{\mathrm{\α}}_k(x,y)}{\underset{k}{\mathrm{\Σ}}{\mathrm{\σ}}_k^2}---\left(14\right)$

式中：

表示融合后的图像；

表示输入图像k的局部方差；α_k(x，y)表示输入图像k当前分析点的像素值。

步骤三：统计模型参数；

假设存在理想场景R，光学传感器和SAR传感器均是对理想场景的部分反映，其反映的内容分别为β_oR和β_sR，根据背景技术中所介绍的SAR图像与光学图像的统计模型可知，各输入图像满足如下的方程组：

$\left\{\begin{array}{c}{a}_{o1}(x,y)={\mathrm{\β}}_{o1}(x,y)R(x,y)+{\mathrm{\ϵ}}_{o1}(x,y)\\ .\\ .\\ .\\ a_{\mathrm{on}}(x,y)={\mathrm{\β}}_{\mathrm{on}}(x,y)R(x,y)+{\mathrm{\ϵ}}_{\mathrm{on}}(x,y)\\ a_{s1}(x,y)={\mathrm{\β}}_{s1}(x,y)R(x,y)\·{\mathrm{\ϵ}}_{s1}(x,y)\\ .\\ .\\ .\\ a_{\mathrm{sm}}(x,y)={\mathrm{\β}}_{\mathrm{sm}}(x,y)R(x,y)\·{\mathrm{\ϵ}}_{\mathrm{sm}}(x,y)\end{array}\right.---\left(15\right)$

式中：a_oi(x，y)表示第i个光学传感器获取的遥感图像在(x，y)处的像素值，a_sj(x，y)表示第j个SAR传感器获取的遥感图像在(x，y)处的像素值，β_oi(x，y)表示第i个光学传感器在(x，y)处的增益因子，β_sj(x，y)表示第j个SAR传感器在(x，y)处的增益因子，ε_oi(x，y)表示第i个光学传感器遥感图像中在(x，y)处的噪声，ε_sj(x，y)表示第i个SAR传感器遥感图像中在(x，y)处的噪声，光学图像中的噪声满足高斯分布，SAR图像中的噪声满足伽马分布；i∈[1，n]，j∈[1，m]，n表示光学图像的数量，m表示SAR图像的数量，R表示理想场景，R(x，y)表示(x，y)处的理想场景，E(·)表示数学期望。。

采用Lee降斑滤波器技术，对SAR图像的统计模型进行一阶泰勒展开，将方程组(15)将简化为：

A＝HR+W                (16)

式中：A＝(a_o1(x，y)，…，a_on(x，y)，a_s1(x，y)，…，a_sm(x，y))^T；

H＝(β_o1(x，y)，…，β_on(x，y)，β_s1(x，y)，…，β_sm(x，y))^T；

$W=\left(\begin{array}{c}{\mathrm{\ϵ}}_{o1}(x,y)\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_{\mathrm{on}}(x,y)\\ E\left({\mathrm{\β}}_{s1}(x,y)R(x,y)\right)\·{\mathrm{\ϵ}}_{s1}(x,y)\\ .\\ .\\ .\\ E\left({\mathrm{\β}}_{\mathrm{sm}}(x,y)R(x,y)\right)\·{\mathrm{\ϵ}}_{\mathrm{sm}}(x,y)\end{array}\right)$

在确定局部分析窗口后，依据输入图像的统计模型及局部分析窗口，统计局部图像的模型参数，所涉及模型参数主要包括：光学图像的噪声方差、SAR图像中的噪声方差、光学传感器的增益因子和SAR传感器的增益因子等因素的估计处理。

对于光学图像中噪声方差的估计需根据输入图像的特性采用不同的估计方法，如果图像以平坦区域为主，对于平坦区域近似认为局部方差主要由噪声所引起，利用局部方差直方图中的峰值所对应的方差来近似图像中的噪声方差

如果图像不以平坦区域为主但存在一些平坦区域，由于边缘和纹理的影响，其局部方差直方图的峰值并不能很好的体现图像的噪声方差，则人为的选取平坦区域，通过统计平坦区域的局部方差来近似图像中的噪声方差

如果图像中不存在平坦区域或者不希望人工进行参与，直接将局部方差最小的区域近似为平坦区域，利用统计局部方差来近似输入图像的噪声方差

对于SAR图像中噪声方差的估计需要考虑SAR图像斑点噪声乘性模型的特性，根据式(16)给出的一阶泰勒近似展开后SAR图像中的噪声表达形式，结合SAR图像的统计特性可知，SAR图像中的局部噪声方差为：

${\mathrm{\σ}}_{\mathrm{\ϵs}}^2=E{\left(a_s(x,y)\right)}^2/L---\left(17\right)$

式中：

表示SAR图像的局部噪声方差；E(·)表示数学期望；a_s(x，y)表示输入的SAR图像；L为SAR图像的等效视数。

假设局部分析窗口内传感器增益因子保持不变，针对方程(16)两侧分别进行求方差处理后有：

$C_A=D_H^2{H}^{T}H+{\mathrm{\Σ}}_W---\left(18\right)$

式中：C_A表示输入图像的协方差矩阵；H传感器增益向量；∑_W噪声的协方差矩阵；

表示理想场景的方差。

同时，根据局部统计可以得到输入图像的局部方差为：

${\mathrm{\Σ}}_A=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(A_n-\stackrel{\‾}{A}){(A_n-\stackrel{\‾}{A})}^T---\left(19\right)$

其中：

表示局部均值；A_n表示输入图像。

两者之间的均方误差E可以表示为：

$E={\left|\right|{\mathrm{\Σ}}_A-D_R^2{H}^{T}H-{\mathrm{\Σ}}_W\left|\right|}^2---\left(20\right)$

其中：||·||表示模值运算。

为了使根据模型估算的传感器图像局部方差和根据实际图像估算的局部方差间具有最小均方误差，将上式对传感器增益因子H进行求导，并使得求导式为零，可以得到如下方程：

$({\mathrm{\Σ}}_A-{\mathrm{\Σ}}_W)H^T=D_R^2{H}^T{\mathrm{HH}}^T---\left(21\right)$

通过式(21)可以得到传感器增益因子H为：

$H=\frac{{\mathrm{\λ}}^{0.5}}{D_R}U---\left(22\right)$

其中：U和λ分别表示协方差矩阵(∑_A-∑_W)所对应的特征向量和特征值。

若引入约束条件||H||＝1，则式(22)可简化为：

H＝U                        (23)

步骤四：图像合成处理；

在确定模型参数的基础上，结合式(16)给出的图像统计模型，可以利用最小均方误差估计理论完成对理想场景的最优化估计，实现输入图像的融合处理，此时，融合处理后的图像可表示为：

$R=E\left(R\right)+{(D_R^{-1}+H^T{\mathrm{\Σ}}_w^{-1}H)}^{-1}{H}^T{\mathrm{\Σ}}_w^{-1}(A-H\·E\left(R\right))---\left(24\right)$

将步骤三中估算得到的模型参数带入到式(24)中，得到输入场景的最优估计结果，实现SAR图像与光学图像间的融合处理。

本发明的优点在于：

(1)本发明从SAR图像和光学图像的统计特性出发，综合考虑SAR图像斑点噪声及图像中的边缘特征对融合处理的影响，有效改善了SAR图像与光学图像及SAR图像之间的融合处理效果；

(2)本发明利用Lee降斑滤波器技术，将SAR图像与光学图像的融合处理从非线性方程组的最优化估计处理转换为线性方程组的最优化估计处理，简化融合处理过程；

(3)本发明利用局部边缘轮廓特征完成局部分析窗口的自适应修正处理，进而保证局部分析窗口内为同质目标，提高模型参数估计的精度，进而提升融合处理效果。

附图说明

图1是本发明背景技术中遥感图像融合处理的方法流程图；

图2是本发明的方法流程图；

图3a是本发明步骤一中局部分析窗口水平方向划分示意图；

图3b是本发明步骤一中局部分析窗口垂直方向划分示意图；

图3c是本发明步骤一中局部分析窗口右倾斜45°划分示意图；

图3d是本发明步骤一中局部分析窗口左倾斜45°划分示意图；

图4是本发明实施例中输入的光学图像；

图5是本发明实施例中输入的SAR图像；

图6a是本发明实施例中基于最大值准则的融合处理结果；

图6b是本发明实施例中基于局部方差准则的融合处理结果；

图6c是本发明实施例中基于显著性水平准则的融合处理结果；

图6d是本发明施例中采用本发明的融合处理结果；

图7a是本发明施例中第一组基于最大值准则融合处理结果的局部放大图；

图7b是本发明施例中第一组基于局部方差准则融合处理结果的局部放大图；

图7c是本发明施例中第一组基于显著性水平准则融合处理结果的局部放大图；

图7d是本发明施例中第一组采用本发明的融合处理结果的局部放大图；

图8a是本发明施例中第二组基于最大值准则融合处理结果的局部放大图；

图8b是本发明施例中第二组基于局部方差准则融合处理结果的局部放大图；

图8c是本发明施例中第二组基于显著性水平准则融合处理结果的局部放大图；

图8d是本发明施例中第二组采用本发明的融合处理结果的局部放大图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于统计模型的SAR图像融合处理方法，流程如图2所示，包括如下几个步骤：

步骤一：提取输入图像的边缘特征；

边缘轮廓特征是图像中的重要特征，也是图像融合处理中需要保留的重要特性，然而，边缘特征的存在会导致局图分析窗口的统计特性不满足式(6)和式(8)中所给出的遥感图像统计模型，进而影响遥感图像融合处理效果。同时，当前分析点位于图像边缘时，局部分析窗口跨越两类介质，无法保证局部分析窗口内的传感器增益因子保持不变，也将影响遥感图像融合处理效果。为了减小边缘特征对遥感图像融合处理的影响，并更好的保留输入图像中的边缘轮廓特征，需要预先提取图像中的边缘特征，并以此为基础来修正局部分析窗口，进而提高融合处理的效果。

所述的输入图像为SAR图像或者光学图像；

(1)当输入图像为光学图像时，采用Canny边缘检测算法进行边缘轮廓提取，具体为：

①首先将输入图像与二维高斯函数作卷积处理，利用高斯函数的低通特性来完成对输入图像的滤波处理。

②对滤波处理后的图像实施微分处理，分别提取水平方向和垂直方向上的处理结果G_y和G_x，得到梯度的大小和方向：

$G=\sqrt{G_x^2+G_y^2}$

$\mathrm{\θ}=\mathrm{Arc}\tan \frac{G_y}{G_x}$

式中：G和θ分别表示梯度的大小与方向。

③对微分处理后的图像进行非极大值抑制处理，即以梯度方向为参考，判断该像素点梯度值是否为沿梯度方向上的最大梯度值，若该梯度值为最大时，表明该像素点的梯度值为局部极大值，保留该像素的梯度值；若该梯度值不是最大时，表明该像素点的梯度值并非局部极大值，将梯度值设置为0。

④在完成非极大值抑制处理后，进行双阈值处理。统计微分处理后图像的均值与方差，将均值与两倍方差的和设定为高阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制后图像的高阈值边缘图像；将均值设定为低阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制后图像的低阈值边缘图像。然后以高阈值边缘图像为种子，通过在对应低阈值边缘图像中进行搜索，获得完整的光学图像边缘特征。

(2)当输入图像为SAR图像时，采用Touzi边缘检测算法进行边缘轮廓提取，具体为：

①以SAR图像中的当前分析点为中心，采用如图3所示的掩模模版，图中标有“中”字的点为当前分析点，将局部分析窗口分别按照水平方向、垂直方向、右倾斜45°和左倾斜45°的方式划分为两个区域，分别统计两个区域的均值大小，计算水平方向、垂直方向、右倾斜45°和左倾斜45°四个方向上均值的比值γ_0°，γ_45°，γ_90°，γ_135°。

$\mathrm{\γ}=\max \{\frac{{\mathrm{\μ}}_1}{{\mathrm{\μ}}_2},\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}\},(\mathrm{\γ}\≥1)$

式中：μ₁表示左侧或上侧区域均值；μ₂表示右侧或下侧区域均值；γ表示分区后两部分较大均值与较小均值的比值。

选择四个方向上最大的比值作为当前分析点的均值比值，并记录最大的比值所对应的方向，作为当前分析点的梯度方向。

γ_value＝Max(γ_0°，γ_45°，γ_90°，γ_135°)

γ_direction＝direction(γ_value)

式中：γ_value表示当前分析点的均值比值；γ_direction表示当前分析点的梯度方向。

②对均值比值处理结果进行非极大值抑制处理，遍历γ_value，根据每一个点的梯度方向γ_direction，对相邻像素均值比值的大小进行比较，以判断当前分析点的均值比值是否为沿梯度方向上的局部极大值，若沿该像素点梯度方向上前后两个像素的均值比值中的任何一个均值比值大于当前像素的均值比值，表明该像素点的均值比值并不是梯度方向上的极大值，将该像素点处的均值比值设置为0；若沿该像素点梯度方向上前后两个像素的均值比值均小于当前像素的均值比值，表明该像素点的均值比值是梯度方向上的极大值，保留该像素点处的均值比值。最后得到像素均值比值的非极大值抑制处理结果；

③对非极大值抑制处理结果进行双阈值处理。统计均值比值处理后图像的均值与方差，将均值与两倍方差的和设定为高阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制处理后图像的高阈值边缘图像；将均值设定为低阈值，对非极大值抑制处理后图像进行阈值处理，将幅度值低于阈值的幅度值设定为0，幅度值高于阈值的幅度值保留，获得非极大值抑制处理后图像的低阈值边缘图像。然后以高阈值边缘图像为种子，通过在对应低阈值边缘图像中进行搜索，获得完整的SAR图像边缘特征。

步骤二：局部分析窗口自适应修正；

在完成边缘特征提取后，利用局部边缘轮廓与当前分析点之间的相互关系来完成对局部分析窗口的自适应修正处理，进而确保局部分析窗口内不存在边缘轮廓特征，整个局部分析窗口内为同质目标。局部分析窗口的自适应修正处理具体包括以下几个步骤：

1)以输入图像的当前分析点为中心，以设定的分析窗口尺寸为半径，截取局部圆形区域作为局部分析窗口；

2)根据步骤一得到的图像边缘特征，检测局部分析窗口内是否存在边缘轮廓，当各输入图像的对应局部分析窗口内均不存在边缘轮廓时，认为当前分析点处于均匀区域，直接利用原局部分析窗口作为当前分析点的局部分析窗口；

3)当局部分析窗口内存在边缘轮廓，且当前分析点不位于边缘轮廓上时，判断当前待分析点与边缘轮廓之间的相互位置关系，提取出局部分析窗口内包含当前分析点的边缘轮廓的内接局部均匀区域作为新的局部分析窗口；

4)当局部分析窗口内存在边缘轮廓特征时，且当前分析点位于部分或全部输入图像的边缘轮廓上时，忽略当前分析点不位于边缘轮廓上的输入图像，仅仅考虑当前分析点位于边缘轮廓上的输入图像。根据当前分析点处的局部方差进行加权平均处理，得到当前分析点的融合处理结果；

$\widehat{\mathrm{\α}}(x,y)=\frac{\underset{k}{\mathrm{\Σ}}{\mathrm{\σ}}_k^2{\mathrm{\α}}_k(x,y)}{\underset{k}{\mathrm{\Σ}}{\mathrm{\σ}}_k^2}$

式中：

表示融合后的图像；

表示输入图像k的局部方差；α_k(x，y)表示输入图像k的当前分析点的像素值。

步骤三：统计融合参数

在确定局部分析窗口后，依据输入图像的统计模型及局部分析窗口，统计局部图像的模型参数。具体处理过程包含以下几个步骤：

1)采用Lee降斑滤波器技术，简化输入图像的统计模型；

假设存在理想场景R，光学传感器和SAR传感器均是对理想场景的部分反映，其反映的内容分别为β_oR和β_sR，各输入图像满足如下的方程组：

$\left\{\begin{array}{c}{a}_{o1}(x,y)={\mathrm{\β}}_{o1}(x,y)R(x,y)+{\mathrm{\ϵ}}_{o1}(x,y)\\ .\\ .\\ .\\ a_{\mathrm{on}}(x,y)={\mathrm{\β}}_{\mathrm{on}}(x,y)R(x,y)+{\mathrm{\ϵ}}_{\mathrm{on}}(x,y)\\ a_{s1}(x,y)={\mathrm{\β}}_{s1}(x,y)R(x,y)\·{\mathrm{\ϵ}}_{s1}(x,y)\\ .\\ .\\ .\\ a_{\mathrm{sm}}(x,y)={\mathrm{\β}}_{\mathrm{sm}}(x,y)R(x,y)\·{\mathrm{\ϵ}}_{\mathrm{sm}}(x,y)\end{array}\right.$

其中：a_oi(x，y)表示第i个光学传感器获取的遥感图像在(x，y)处的像素值，a_sj(x，y)表示第j个SAR传感器获取的遥感图像在(x，y)处的像素值，β_oi(x，y)表示第i个光学传感器在(x，y)处的增益因子，β_sj(x，y)表示第j个SAR传感器在(x，y)处的增益因子，ε_oi(x，y)表示第i个光学传感器遥感图像中在(x，y)处的噪声，ε_sj(x，y)表示第i个SAR传感器遥感图像中在(x，y)处的噪声，光学图像中的噪声满足高斯分布，SAR图像中的噪声满足伽马分布；i∈[1，n]，j∈[1，m]，n表示光学图像的数量，m表示SAR图像的数量，R表示理想场景，R(x，y)表示(x，y)处的理想场景，E(·)表示数学期望。

采用Lee降斑滤波器技术，对SAR图像的统计模型进行一阶泰勒展开，方程组(简化输入图像的统计模型)可简化为：

A＝HR+W

其中：A＝(a_o1(x，y)，…，a_on(x，y)，a_s1(x，y)，…，a_sm(x，y))^T；

H＝(β_o1(x，y)，…，β_on(x，y)，β_s1(x，y)，…，β_sm(x，y))^T；

$W=\left(\begin{array}{c}{\mathrm{\ϵ}}_{o1}(x,y)\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_{\mathrm{on}}(x,y)\\ E\left({\mathrm{\β}}_{s1}(x,y)R(x,y)\right)\·{\mathrm{\ϵ}}_{s1}(x,y)\\ .\\ .\\ .\\ E\left({\mathrm{\β}}_{\mathrm{sm}}(x,y)R(x,y)\right)\·{\mathrm{\ϵ}}_{\mathrm{sm}}(x,y)\end{array}\right).$

此时，为了确定输入图像的统计模型，需要利用局部分析窗口内的图像数据估算模型参数，所涉及的参数包括光学图像的噪声方差、SAR图像中的噪声方差、光学传感器的增益因子和SAR传感器的增益因子等，对应各参数的估计方法如下：

2)光学图像中的噪声方差估计；

光学图像中噪声方差的估计方法根据输入图像的特性采用不同的估计方法：如果图像以平坦区域为主，对于平坦区域近似认为局部方差主要由噪声所引起，利用局部方差直方图中的峰值所对应的方差来近似图像中的噪声方差

如果图像不以平坦区域为主但存在一些平坦区域，由于边缘和纹理的影响，其局部方差直方图的峰值并不能很好的体现图像的噪声方差，则人为的选取平坦区域，通过统计平坦区域的局部方差来近似图像中的噪声方差

如果图像中不存在平坦区域或者不希望人工参与，直接将局部方差最小的区域近似为平坦区域，利用统计局部方差来近似输入图像的噪声方差

3)SAR图像中的噪声方差估计；

与光学图像不同，由于斑点噪声是以乘性噪声模式引入的，局部噪声强度与局部能量有关，因此，必须逐点计算SAR图像中各点的噪声强度。结合SAR图像统计模型的一阶泰勒近似展开模型可知，SAR图像中的局部噪声方差为：

${\mathrm{\σ}}_{\mathrm{\ϵs}}^2=E{\left(a_s(x,y)\right)}^2/L$

其中：表示SAR图像的局部噪声方差；E(·)表示数学期望；a_s(x，y)表示输入的SAR图像；L为SAR图像的等效视数。

4)传感器增益因子的估算；

本发明利用根据模型估算的传感器图像局部方差和根据实际图像估算的局部方差间具有最小均方误差为准则来对传感器增益因子进行估算。

假设局部分析窗口内传感器增益因子保持不变，针对理想场景估计方程式的两侧分别进行求方差处理后有：

$C_A=D_R^2{H}^{T}H+{\mathrm{\Σ}}_W$

式中：C_A表示输入图像的协方差矩阵；H传感器增益向量；∑_W噪声的协方差矩阵；

表示理想场景的方差。

同时，根据局部统计可以得到输入图像的局部方差为：

${\mathrm{\Σ}}_A=\frac{1}{N}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}(A_n-\stackrel{\‾}{A}){(A_n-\stackrel{\‾}{A})}^T$

其中：

表示局部均值；A_n表示输入图像。

两者之间的均方误差E可以表示为：

$E={\left|\right|{\mathrm{\Σ}}_A-D_R^2{H}^{T}H-{\mathrm{\Σ}}_W\left|\right|}^2$

其中：||·||表示模值运算。

为了获取最小均方误差，将上式对传感器增益因子H进行求导，并使得求导式为零，可以得到如下方程：

$({\mathrm{\Σ}}_A-{\mathrm{\Σ}}_W)H^T=D_R^2{H}^T{\mathrm{HH}}^T$

进而可以得到传感器增益因子H为：

$H=\frac{{\mathrm{\λ}}^{0.5}}{D_R}U$

其中：U和λ分别表示协方差矩阵(∑_A-∑_W)所对应的特征向量和特征值。

若引入约束条件||H||＝1，传感器增益因子H可简化为：

H＝U

步骤四：图像合成处理；

在确定模型参数的基础上，利用最小均方误差估计理论完成对理想场景的最优化估计，实现输入图像的融合处理，此时，融合处理后的图像可表示为：

$R=E\left(R\right)+{(D_R^{-1}+H^T{\mathrm{\Σ}}_w^{-1}H)}^{-1}{H}^T{\mathrm{\Σ}}_w^{-1}(A-H\·E\left(R\right))$

将估算得到的传感器增益因子H、光学图像的噪声方差

和SAR图像的噪声方差

带入到理想场景估计方程中，得到输入场景的最优估计结果，实现SAR图像与光学图像间的融合处理。

实施例：

输入的光学图像如图4所示，输入的SAR图像如图5所示，图6a～6d给出了采用不同融合处理准则的处理结果，其中：图6a为采用最大值准则的融合处理结果，图6b为采用局部方差准则的融合处理结果，图6c为采用局部显著度准则的融合处理结果，图6d为采用本发明的融合处理结果。由于图像中斑点噪声的影响导致前三种处理算法的处理结果主要体现光学图像的特性，而SAR图像不但没有为融合结果中注入互补信息，反而导致融合处理结果变得更加模糊。本发明有效考虑斑点噪声的影响，将SAR图像中的互补信息注入到融合处理结果中，进而获得更为完成的场景信息。

为了更好的对比各种处理算法所获得的融合处理效果，图7给出了一组局部特征放大结果，其中图7a为采用最大值准则的融合处理结果的局部放大图，图7b为采用局部方差准则的融合处理结果的局部放大图，图7c为采用局部显著度准则的融合处理结果的局部放大图，图7d为采用本发明的融合处理结果的局部放大图。如图所示，在前面三种算法的处理结果中，由于SAR图像中斑点噪声的影响导致融合图像中具有较强的噪声，而本发明的处理结果噪声水平明显低于前三种融合处理结果，且同时包含SAR图像和光学图像的特性，如融合处理结果中包含了SAR图像中所显示的机场边上的铁丝网，而该信息在光学图像中很难被发现。

图8给出了另一组局部特征放大结果，其中图8a为采用最大值准则的融合处理结果的局部放大图，图8b为采用局部方差准则的融合处理结果的局部放大图，图8c为采用局部显著度准则的融合处理结果的局部放大图，图8d为采用本发明的融合处理结果的局部放大图。如图所示，SAR图像的穿透特性使得局部特征更为清晰。

为了进一步衡量不同融合处理算法的融合处理效果，引入峰值信噪比和相关系数这两个指标来描述融合处理效果，其中：峰值信噪比反映了融合图像中的噪声大小，通过对比不同融合处理算法所得到融合图像中的峰值信噪比，能够比较各种融合算法对噪声的抑制能力，其值越大表明融合处理算法的噪声抑制能力越强。相关系数反映了两幅图像的相关程度，通过对比融合前后图像的相关系数，能够评价融合处理算法对图像纹理特征的保持能力。表1给出了针对上述4种处理算法所获得融合处理结果的评估结果。对比四种处理算法可见，本发明融合处理算法不仅具有最强的噪声抑制能力，且信息保持能力也优于其余各种融合处理算法。

表1融合效果的评估结果

	最大值准则	方差准则	显著度准则	本报告算法
					峰值信噪比(dB)	38.5	38.0	38.3	40.5
相关系数	0.651	0.622	0.651	0.654

Claims (2)

1.基于统计模型的SAR图像融合处理方法，其特征在于，包括以下几个步骤：

步骤一：提取输入图像的边缘特征；

所述的输入图像为SAR图像或者光学图像；

(1)当输入图像为光学图像时，采用Canny边缘检测算法进行边缘轮廓提取；首先将输入图像f(x，y)与二维高斯函数G(x，y)作卷积处理，利用高斯函数的低通特性完成对输入图像的滤波处理；对滤波处理后的图像实施微分处理，分别提取水平方向和垂直方向上的处理结果G_y和G_x，得到梯度的大小和方向；对微分处理后的图像进行非极大值抑制处理，判断像素点梯度值是否为沿梯度方向上的极大梯度值，若梯度值为极大值时，保留梯度值；若梯度值不是极大值时，将梯度值设置为0；对非极大值抑制处理后的边缘图像进行双阈值处理，分别获得高阈值和低阈值边缘图像，以高阈值边缘图像为种子，在对应低阈值边缘图像中进行搜索，获得完整的光学图像边缘特征；

(2)当输入图像为SAR图像时，采用Touzi边缘检测算法进行边缘轮廓提取；首先以SAR图像中的当前分析点为中心，分别提取水平方向、垂直方向、右倾斜45°和左倾斜45°四个方向上均值的比值，并以四个方向上比值的最大值作为当前分析点的梯度值，所对应的方向为梯度方向；对梯度处理结果进行非极值抑制处理，根据每一个点的梯度方向，对相邻像素均值比值大小进行比较，若沿梯度方向前后两个像素的均值比值有一个大于当前像素的均值比值，表明当前分析点并非局部极大值，将对应的均值比值设置为0；若沿梯度方向前后两个像素的均值比值均小于当前像素的均值比值，表明当前分析点为局部极大值，保留对应的均值比值；对非极大值抑制处理结果进行双阈值处理，分别获得高阈值和低阈值边缘图像，然后以高阈值边缘图像为种子，在对应低阈值图像中进行搜索，获得完整的SAR图像边缘特征；

步骤二：局部分析窗口自适应修正；

以输入图像的当前分析点为中心，以设定的分析窗口尺寸为半径，截取局部圆形区域作为局部分析窗口；检测局部分析窗口内是否存在边缘轮廓，若各输入图像的对应局部分析窗口内均不存在边缘轮廓时，认为当前分析点处于均匀区域，直接利用原局部分析窗口作为当前分析点的局部分析窗口；若当前局部分析窗口内存在边缘特征，且当前分析点在所有输入图像中均不位于边缘特征上时，判断当前分析点与边缘轮廓之间的相互位置关系，提取出局部分析窗口内包含当前分析点的边缘轮廓的内接局部均匀区域作为新的局部分析窗口；若当前分析点在部分或全部输入图像中位于边缘轮廓上时，仅考虑当前分析点位于边缘轮廓上的输入图像，根据当前分析点处的局部方差进行加权平均处理，得到当前分析点的融合处理结果；

步骤三：统计融合参数；

采用Lee降斑滤波器技术，简化输入图像的统计模型，修正后图像模型涉及参数包括光学图像的噪声方差、SAR图像中的噪声方差、光学传感器的增益因子和SAR传感器的增益因子；

光学图像中噪声方差的估计方法根据输入图像的特性采用不同的估计方法，如果图像以平坦区域为主，利用局部方差直方图中的峰值所对应的方差来近似图像中的噪声方差

如果图像不以平坦区域为主但存在一些平坦区域，人为的选取平坦区域，通过统计平坦区域的局部方差来近似图像中的噪声方差

如果图像中不存在平坦区域或者不希望人工参与，直接将局部方差最小值近似输入图像的噪声方差

SAR图像中噪声方差的估计方法根据斑点噪声乘性模型的特点，利用局部图像均值的平方与图像等效视数的比值来估算得到局部SAR图像的噪声方差；

以根据模型估算的传感器图像局部方差和根据实际图像估算的局部方差间具有最小均方误差为准则，将根据两种方式得到局部方差的差值对输入传感器增益因子进行求偏导处理，并令求导式为零来计算输入传感器的增益因子；

步骤四：图像合成处理；

在确定模型参数的基础上，输入图像的融合处理过程可转换为根据输入图像完成对理想场景的最优化估计过程，依据输入的图像组采用最小均方误差估计得到理想场景的最优化估计结果，完成对输入图像的融合处理。

2.根据权利要求1所述的基于统计模型的SAR图像融合处理方法，其特征在于：所述步骤三中，简化输入图像的统计模型具体为：

A＝HR+W

其中：A＝(a_o1(x，y)，…，a_on(x，y)，a_s1(x，y)，…，a_sm(x，y))^T；

H＝(β_o1(x，y)，…，β_on(x，y)，β_s1(x，y)，…，β_sm(x，y))^T；

$W=\left(\begin{array}{c}{\mathrm{\ϵ}}_{o1}(x,y)\\ .\\ .\\ .\\ {\mathrm{\ϵ}}_{\mathrm{on}}(x,y)\\ E\left({\mathrm{\β}}_{s1}(x,y)R(x,y)\right)\·{\mathrm{\ϵ}}_{s1}(x,y)\\ .\\ .\\ .\\ E\left({\mathrm{\β}}_{\mathrm{sm}}(x,y)R(x,y)\right)\·{\mathrm{\ϵ}}_{\mathrm{sm}}(x,y)\end{array}\right)$

其中：x，y表示图像中各像素的位置，a_oi(x，y)表示第i个光学传感器获取的遥感图像在(x，y)处的像素值，a_sj(x，y)表示第j个SAR传感器获取的遥感图像在(x，y)处的像素值，β_oi(x，y)表示第i个光学传感器在(x，y)处的增益因子，β_sj(x，y)表示第j个SAR传感器在(x，y)处的增益因子，ε_oi(x，y)表示第i个光学传感器遥感图像中在(x，y)处的噪声，ε_sj(x，y)表示第i个SAR传感器遥感图像中在(x，y)处的噪声，光学图像中的噪声满足高斯分布，SAR图像中的噪声满足伽马分布；i∈[1，n]，j∈[1，m]，n表示光学图像的数量，m表示SAR图像的数量，R表示理想场景，R(x，y)表示(x，y)处的理想场景，E(·)表示数学期望。

Images