CN105354798A - 基于几何先验和分布相似性测度的sar图像去噪方法 - Google Patents

Abstract

本发明公开了一种基于几何先验和分布相似性测度的SAR图像去噪方法。主要解决现有技术中无法有效调节SAR图像滤波中细节信息保持和相干斑抑制的矛盾。其实现步骤如下：(1)提取SAR图像的素描图；(2)按照素描图将合成孔径雷达SAR图像中的像素分为具有方向特性的像素和不具有方向特性的像素；(3)估计具有方向特性像素的像素值；(4)估计不具有方向特性像素的像素值；(5)对于步骤(3)和步骤(4)中的重叠像素采用基于统计加权的方式对估计值进行融合。本发明具有较好的相干斑抑制效果，可用于合成孔径雷达SAR图像的去噪。

Description

基于几何先验和分布相似性测度的SAR图像去噪方法

技术领域

本发明属于图像处理技术领域，更进一步涉及合成孔径雷达SAR图像滤波处理技术领域中的一种基于几何先验和分布相似性测度的合成孔径雷达图像(SyntheticApertureRadar，SAR)的去噪方法。本发明可用于对合成孔径雷达SAR图像的相干斑进行抑制。

背景技术

合成孔径雷达图像是采用主动发射电磁波并接收电磁回波的方式成像的。然而，由于反向电磁回波之间常常伴随有相干现象，所获得的合成孔径雷达SAR图像中存在有大量的斑点信息。虽然这些相干斑信息在一定程度上体现了合成孔径雷达SAR图像成像场景中目标的物理几何特性，但是也给合成孔径雷达SAR图像的理解和解译造成了巨大的困难。因此，相干斑抑制是对合成孔径雷达SAR图像进行理解和解译的关键技术之一。

在图像的噪声抑制中，频域和空域是比较常见的用于对图像滤波方法进行分类的准则。其中，频域滤波方法是通过将信号变换到频域，利用真实信号和噪声在频域具有不同的频谱特性这一特点实现对真实信号的估计。空域滤波方法则是以局部邻域像素之间的统计相关性为基础，实现对当前像素真实值的估计。

HongxiaoFeng，BiaoHou，MaoguoGong在其发表的论文“SARImageDespecklingBasedonLocalHomogeneous-RegionSegmentationbyUsingPixel-RelativityMeasurement.IEEETransactionsonGeoscienceandRemoteSensing,2011,49(7)”中提出了一种使用像素间的相关性测度来构建块的相似性测度进行局部同质区域搜索的合成孔径雷达SAR图像去噪方法。该方法使用比值距离来度量两个含噪图像块之间的距离，通过核函数得到局部匀质区域，然后使用极大似然准则来估计每个匀质区域中的真实信号得到去噪图像。但是该方法仍然存在的不足之处是，没有考虑合成孔径雷达SAR图像中不同图像块内像素之间的统计相关性，从而导致去噪结果中边缘细节信息的模糊泛化。

西安电子科技大学在其申请的专利“基于素描图和核选择的SAR图像降斑方法”(专利申请号201310093148.X，公开号CN103377465A)中公开了一种基于素描图和核选择的合成孔径雷达SAR图像降斑方法。该方法根据合成孔径雷达SAR图像的素描图对合成孔径雷达SAR图像进行分类，利用所划分的区域的几何结构特性构造具有方向特性的各向异性核映射函数，对每一个像素以块相似性测度来搜索局部最大同质区域，并用该同质区域的极大似然值作为当前像素的估计值来实现对整个合成孔径雷达SAR图像的相干斑抑制。但是该方法仍然存在的不足之处是，没有考虑到SAR图像边缘结构信息的稀疏性，仍采用基于局部同质区域搜索的滤波方法来估计信号的真实值，最终导致细节信息的保持不好，尤其是点目标的保持。

发明内容

本发明的目的在于克服现有技术的不足，提出一种基于几何先验和分布相似性测度的SAR图像去噪方法。克服了当前合成孔径雷达SAR图像相干斑滤波方法中无法同时实现相干斑的有效抑制和细节信息的保持问题，提高了合成孔径雷达SAR图像信号的估计精度。同时，对基于同质区域搜索的滤波方法，提出了一种更有效的测度，改善了同质区域搜索的质量，实现了对真实信号的有效估计。

为实现上述目的，本发明具体实现步骤包括如下：

(1)素描化图像：

采用合成孔径雷达SAR素描模型，获得输入合成孔径雷达SAR图像的素描图；

(2)像素分类：

以素描图中素描线所包含的素描点的位置作为参考位置，将输入合成孔径雷达SAR图像中位置与参考位置相同的像素点作为具有方向特性的像素，将输入合成孔径雷达SAR图像中剩余的像素作为不具有方向特性的像素；

(3)采用基于几何结构块的非局部滤波方法，估计每个具有方向特性像素的像素值；

(4)估计不具有方向特性像素的像素值：

(4a)采用基于局部线性最小均方误差准则的估计方法，对于每个不具有方向特性的像素，在以该像素为中心的3×3个像素大小的邻域窗口内分别计算其初始估计值和初始等效视数；

(4b)从不具有方向特性的像素中选取任意一个未选像素，以该像素为中心，15个像素为边长构造正方形窗口；

(4c)采用区域生长方法，对正方形窗口中的所有像素进行生长和标记，将标记后的所有像素构成局部邻域，所述的区域生长方法中的相似性测度按照下式计算得到：

$Sim(E_p,E_q)={\[\frac{(2N_q-1)\·N_p}{(2N_p-1)\·N_q}\·r^2\]}^{N_p-0.5}\·{\[\frac{(N_p\·r^2+N_q)(2N_q-1)}{2N_q\·(N_p+N_q-1)}\]}^{1-N_p-N_q}$

其中，Sim(·)表示正方形窗口中种子点像素p的估计值与待生长像素q的估计值之间的相似性值，E_p和E_q分别表示正方形窗口中种子点像素p和待生长像素q的估计值，N_p和N_q分别表示正方形窗口中种子点像素p和待生长像素q的等效视数，r表示正方形窗口中种子点像素p的估计值E_p与待生长像素q的估计值E_q之间的比值，r＝E_p/E_q；

(4d)采用基于局部线性最小均方误差准则的估计方法，利用局部邻域内包含的所有像素，更新正方形窗口中心像素的估计值和等效视数；

(4e)采用相似性测度公式，利用正方形窗口中未标记像素的初始估计值和初始等效视数以及正方形窗口中心像素更新后的估计值和等效视数重新计算未标记像素与中心像素之间的相似性值，将相似性值满足区域生长阈值的像素加入到局部邻域并标记，得到扩充后的局部邻域；

(4f)采用极大似然准则，利用扩充后的局部邻域中的所有像素，估计正方形窗口中心像素的像素值；

(4g)判断所有不具有方向特性的像素中是否存在未选像素，若是，执行步骤(4b)，否则，执行步骤(5)；

(5)估计重叠像素的像素值：

(5a)采用下式表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数：

$p_s\left(s\right)=\frac{2L^L}{\mathrm{\Γ}\left(L\right)}{s}^{2L-1}\exp (-{\mathrm{Ls}}^2),s\≥0$

其中，p_s(·)表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数，s表示输入合成孔径雷达SAR图像的相干斑噪声，L表示输入合成孔径雷达SAR图像的视数，Γ(·)表示伽玛函数，exp(·)表示指数函数；

(5b)采用基于统计分布的加权融合方式，计算重叠像素的融合估计值：

$D=\frac{1}{Z}\underset{i}{\Σ}{p}_s\left(\frac{A}{D_i}\right)\·D_i,i\∈\{0,1\}$

其中，D表示重叠像素的融合估计值，Z表示加权融合的归一化因子，∑表示求和操作，i表示像素估计值的编号，A表示重叠像素在输入合成孔径雷达SAR图像中的灰度值，D₀表示重叠像素在步骤(3)得到的估计值，D₁表示重叠像素在步骤(4)得到的估计值，p_s(·)表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数。

本发明与现有的技术相比具有以下优点：

第一，本发明利用所估计的局部等效视数构建了一种基于局部统计分布相似性的测度准则实现局部同质区域的搜索，克服了现有技术没有考虑合成孔径雷达SAR图像块内像素之间的统计相关性，导致了去噪结果中边缘细节信息模糊泛化的问题。采用本发明能够较好地保持图像的边缘细节信息，取得了较为理想的去噪效果。

第二，本发明利用素描图对合成孔径雷达SAR图像中结构信息的表示，将基于几何结构块的非局部滤波方法和基于统计测度的自适应邻域滤波方法有效地结合在一起，克服了现有技术没有考虑合成孔径雷达SAR图像边缘结构稀疏性的缺点，导致了去噪结果中图像的细节信息，尤其是点目标信息的丢失。采用本发明能够较好地保持合成孔径雷达SAR图像的细节信息，提高了图像的去噪效果。

附图说明

图1是本发明的流程图；

图2是计算不具有方向特性像素的像素值的流程图；

图3是分别采用本发明和现有技术的SK-LHRS滤波方法、迭代PPB滤波方法对TerraSAR-X波段1米分辨率的Nordlingerties图像的去噪结果对比图；

图4是分别采用本发明和现有技术的SK-LHRS滤波方法、迭代PPB滤波方法对DRASARX波段3米分辨率的Bedfordshire图像的去噪结果对比图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照附图1，本发明的具体步骤如下。

步骤1，素描化图像。

采用合成孔径雷达SAR素描模型，提出基于几何加权边线检测和多算子融合的SAR图像素描追踪算法，获得输入合成孔径雷达SAR图像的素描图。

所述的合成孔径雷达SAR素描模型及SAR图像素描追踪算法参见JieWu等人于2014年发表在IEEETransactionsonGeoscienceandRemoteSensing期刊上的文章《LocalmaximalhomogeneousregionsearchforSARspecklereductionwithsketch-basedgeometricalkernelfunction》。

合成孔径雷达SAR素描模型包含以下内容：

定义合成孔径雷达SAR图像的可素描性，将合成孔径雷达SAR图像分为可素描部分和不可素描部分，即亮度结构可辨识部分和不可辨识部分。

构造基于几何加权的边线检测算子，提取合成孔径雷达SAR图像中具有不同亮度结构的几何特征。

利用匹配追踪算法实现SAR图像素描图的提取。

基于几何加权边线检测和多算子融合的SAR图像素描追踪算法的具体步骤如下。

构造具有不同方向和尺度的边、线模板，并利用模板的方向和尺度信息构造各向异性高斯函数来计算该模板中每一点的加权系数，其中尺度个数取值为3～5，方向个数取值为18。

按照下式，计算模板不同区域对应在合成孔径雷达SAR图像中像素的均值μ和方差ν：

$\mathrm{\μ}=\frac{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g{A}_g}{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g}$

$v=\frac{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g{(A_g-\mathrm{\μ})}^2}{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g}$

其中，Ω表示模板中的某一区域，g表示区域Ω中点的位置，μ和ν分别表示区域Ω对应在合成孔径雷达SAR图像中像素的均值和方差，∑表示求和操作，w_g表示区域Ω中位置g处的权重系数，w_g的取值范围为0～1，A_g表示区域Ω中位置g对应在合成孔径雷达SAR图像中的像素值。

按照下式，计算合成孔径雷达SAR图像中每个像素对比值算子的响应值：

$R=1-min\{\frac{{\mathrm{\μ}}_a}{{\mathrm{\μ}}_b},\frac{{\mathrm{\μ}}_b}{{\mathrm{\μ}}_a}\}$

其中，R表示合成孔径雷达SAR图像中每个像素对比值算子的响应值，min{·}表示求最小值操作，a和b分别表示模板不同区域的编号，μ_a和μ_b分别表示区域a和区域b对应在合成孔径雷达SAR图像中像素的均值。

按照下式，计算合成孔径雷达SAR图像中每个像素对相关性算子的响应值：

$C=\sqrt{\frac{1}{1+2\·\frac{v_a^2+v_b^2}{{({\mathrm{\μ}}_a+{\mathrm{\μ}}_b)}^2}}}$

其中，C表示表示合成孔径雷达SAR图像中每个像素对相关性算子的响应值，a和b分别表示模板不同区域的编号，v_a和v_b分别表示区域a和区域b对应在合成孔径雷达SAR图像中像素的方差，μ_a和μ_b分别表示区域a和区域b对应在合成孔径雷达SAR图像中像素的均值，表示平方根操作。

按照下式，融合合成孔径雷达SAR图像中像素对比值算子的响应值和合成孔径雷达SAR图像中像素对相关性算子的响应值，计算合成孔径雷达SAR图像中每个像素对各个模板的响应值：

$F=\sqrt{\frac{R^2+C^2}{2}}$

其中，F表示合成孔径雷达SAR图像中每个像素对各个模板的响应值，R和C分别表示合成孔径雷达SAR图像中像素对比值算子和合成孔径雷达SAR图像中像素对相关性算子的响应值，表示平方根操作。

选择具有最大响应值的模板作为合成孔径雷达SAR图像中像素的模板，并将最大响应值作为该像素的强度，将具有最大响应值的模板的方向作为该像素的方向，获得合成孔径雷达SAR图像的边线响应图和方向图。

利用合成孔径雷达SAR图像中每个像素所选择的模板，获得合成孔径雷达SAR图像基于模板的梯度图。

按照下式，将归一化到[0,1]的边线响应图与归一化到[0,1]的基于模板的梯度图进行融合，得到最终的强度图：

$I=\frac{xy}{1-x-y+2xy}$

其中，I表示强度图中的强度值，x表示边线响应图中的值，y表示梯度图中的值。

利用方向图，采用非极大值抑制方法对该强度图进行检测，得到建议草图。

选取建议草图中具有最大强度的像素，将建议草图中与该像素连通的像素连接形成建议线段，得到建议素描图。

利用合成孔径雷达SAR图像的统计分布，计算建议素描图中素描线的编码长度增益CLG：

$CLG=\underset{t}{\overset{m}{\Σ}}\[\frac{A_t^2}{A_{t,0}^2}+ln\left(A_{t,0}^2\right)-\frac{A_t^2}{A_{t,1}^2}-ln\left(A_{t,1}^2\right)\]$

其中，CLG表示建议素描图中素描线的编码长度增益，∑表示求和操作，m表示当前素描线邻域中像素的个数，t表示当前素描线邻域中像素的编号，A_t表示当前素描线邻域中第t个像素的观测值，A_t,0表示在当前素描线不能表示结构信息的假设下，该素描线邻域中第t个像素的估计值，ln(·)表示对数函数，A_t,1表示在当前素描线能够表示结构信息的假设下，该素描线邻域中第t个像素的估计值。

设定阈值T，T的取值范围为5～50，选择CLG>T的建议素描线作为最终素描图中的素描线，获得输入合成孔径雷达SAR图像对应的素描图。

由于素描线是由多个素描线段组成，因此该素描图可以看成是由素描线段组成的图形，其中每条素描线段由具有方向特性的素描点组成，每个素描点的方向为其所在素描线段的方向。

步骤2，像素分类。

从素描图中选取具有方向特性的素描点即素描图中素描线所包含的素描点，将其位置作为参考位置，将输入合成孔径雷达SAR图像中位置与参考位置相同的像素点作为具有方向特性的像素，将输入合成孔径雷达SAR图像中剩余的像素作为不具有方向特性的像素。

步骤3，采用基于几何结构块的非局部滤波方法，估计每个具有方向特性像素的像素值。

其具体估计方法，参见西安电子科技大学的专利申请“基于线状奇异性信息的SAR图像相干斑抑制方法”(授权专利号：ZL201110095598.3)，步骤如下。

将素描图中的每一条线段作为处理单位，以线段上的每一点为中心，沿着该线段的方向构建目标矩形块B₀。

在具有方向特性的像素所组成的区域中，搜索与目标矩形块B₀方向相同的相似矩形块B_j。

按照下式，计算目标矩形块B₀中所有像素的估计值：

$B_0^{\′}=\underset{j}{\Σ}{w}_{0,j}{B}_j$

其中，B′₀表示目标矩形块B₀中所有像素的估计值，j表示相似矩形块B_j的编号，∑表示求和操作，w_0,j表示目标矩形块B₀与相似矩形块B_j之间的相似性权重，B_j表示与目标矩形块B₀方向相同的相似矩形块。

步骤4，估计不具有方向特性像素的像素值。

参照附图2，对本步骤做进一步的描述。

第1步，采用基于局部线性最小均方误差准则的估计方法，对于每个不具有方向特性的像素，在以该像素为中心的3×3个像素大小的邻域窗口内分别计算其初始估计值和初始等效视数。

基于局部线性最小均方误差准则的估计方法的步骤如下：

按照下式，计算待估计像素的估计值：

$E=\frac{1}{n}\underset{k}{\overset{n}{\Σ}}{A}_k+\mathrm{\ξ}(A-\frac{1}{n}\underset{k}{\overset{n}{\Σ}}{A}_k)$

其中，E表示待估计像素的估计值，n表示邻域中所有像素的个数，k表示邻域中像素的编号，∑表示求和操作，A_k表示邻域中像素k在输入合成孔径雷达SAR图像中的灰度值，ξ表示衡量邻域内异质性程度的系数，ξ的取值范围为0～1，A表示待估计像素在输入合成孔径雷达SAR图像中的灰度值。

按照下式，计算等效视数的估计值：

$N=\frac{n}{(n-1){\mathrm{\ξ}}^2+1}\·L$

其中，N表示等效视数的估计值，n表示邻域中所有像素的个数，ξ表示衡量邻域内异质性程度的系数，ξ的取值范围为0～1，L表示输入合成孔径雷达SAR图像的视数。

第2步，从不具有方向特性的像素中选取任意一个未选像素，以该像素为中心，15个像素为边长构造正方形窗口。

第3步，采用区域生长方法，对正方形窗口中的所有像素进行生长和标记，将标记后的所有像素构成局部邻域，所述的区域生长方法中的相似性测度按照下式计算得到：

$Sim(E_p,E_q)={\[\frac{(2N_q-1)\·N_p}{(2N_p-1)\·N_q}\·r^2\]}^{N_p-0.5}\·{\[\frac{(N_p\·r^2+N_q)(2N_q-1)}{2N_q\·(N_p+N_q-1)}\]}^{1-N_p-N_q}$

其中，Sim(·)表示正方形窗口中种子点像素p的估计值与待生长像素q的估计值之间的相似性值，E_p和E_q分别表示正方形窗口中种子点像素p和待生长像素q的估计值，N_p和N_q分别表示正方形窗口中种子点像素p和待生长像素q的等效视数，r表示正方形窗口中种子点像素p的估计值E_p与待生长像素q的估计值E_q之间的比值，r＝E_p/E_q。

采用区域生长方法得到局部邻域的步骤如下：

选取正方形窗口的中心像素作为种子点像素。

设定区域生长的阈值为0.1。

选取种子点像素邻域内的任一像素作为待生长像素。

采用相似性测度公式，计算正方形窗口中种子点像素的初始估计值与待生长像素的初始估计值之间的相似性值。

判断正方形窗口中种子点像素的初始估计值与待生长像素的初始估计值之间的相似性值是否满足区域生长阈值，若满足，则将待生长像素加入种子点像素所在的区域并标记。

在已加入种子点像素所在区域的像素的邻域中，寻找是否存在包含在正方形窗口内但未访问的像素，若存在，将该像素作为下一个待生长像素，计算正方形窗口中种子点像素的初始估计值与待生长像素的初始估计值之间的相似性值并与区域生长的阈值进行判断，若满足区域生长的阈值，则将待生长像素加入种子点像素所在的区域并标记，重复本操作直到正方形窗口中不存在满足条件的待生长像素。

将标记后的所有像素构成局部邻域。

第4步，采用基于局部线性最小均方误差准则的估计方法，利用局部邻域内包含的所有像素，更新正方形窗口中心像素的估计值和等效视数。

第5步，采用相似性测度公式，利用正方形窗口中未标记像素的初始估计值和初始等效视数以及正方形窗口中心像素更新后的估计值和等效视数重新计算未标记像素与中心像素之间的相似性值，将相似性值满足区域生长阈值的像素加入到局部邻域并标记，得到扩充后的局部邻域。

未标记像素是指在区域生长过程中曾选为待生长像素但不包含在局部邻域内的像素。

第6步，采用极大似然准则，利用扩充后的局部邻域中的所有像素，估计正方形窗口中心像素的像素值。

第7步，判断所有不具有方向特性的像素中是否存在未选像素，若是，执行第2步，否则，执行步骤5。

步骤5，估计重叠像素的像素值。

采用下式表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数：

$p_s\left(s\right)=\frac{2L^L}{\mathrm{\Γ}\left(L\right)}{s}^{2L-1}\exp (-{\mathrm{Ls}}^2),s\≥0$

其中，p_s(·)表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数，s表示输入合成孔径雷达SAR图像的相干斑噪声，L表示输入合成孔径雷达SAR图像的视数，Γ(·)表示伽玛函数，exp(·)表示指数函数。

采用基于统计分布的加权融合方式，计算重叠像素的融合估计值：

$D=\frac{1}{Z}\underset{i}{\Σ}{p}_s\left(\frac{A}{D_i}\right)\·D_i,i\∈\{0,1\}$

其中，D表示重叠像素的融合估计值，Z表示加权融合的归一化因子，∑表示求和操作，i表示像素估计值的编号，A表示重叠像素在输入合成孔径雷达SAR图像中的灰度值，D₀表示重叠像素在步骤3得到的估计值，D₁表示重叠像素在步骤4得到的估计值，p_s(·)表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数。

重叠像素是指包含在目标矩形块内的不具有方向特性的像素。

下面结合仿真图对本发明做进一步的描述。

1.仿真条件：

图像来源为DRASARX波段3米分辨率的Bedfordshire图像(记为Field)和TerraSAR-X波段1米分辨率的Nordlingerties图像(记为Nord)。

2.仿真内容：

在以上仿真条件下，使用两种不同的对比方法对Field图像和Nord图像进行去噪：

对比方法1为SK-LHRS滤波方法，对Field图像和Nord图像进行去噪；

对比方法2为PPB滤波方法，对Field图像和Nord图像进行去噪；

两种对比方法和本发明方法的去噪视觉效果如图3和图4所示，其中图3(a)为来源于DRASARX波段3米分辨率的Bedfordshire图像(记为Field)，图3(b)为SK-LHRS滤波方法的去噪结果图，图3(c)为PPB滤波方法迭代25次的去噪结果图，图3(d)为本发明的去噪结果图，图4(a)为来源于TerraSAR-X波段1米分辨率的Nordlingerties图像(记为Nord)，图4(b)为SK-LHRS滤波方法的去噪结果图，图4(c)为PPB滤波方法迭代25次的去噪结果图，图4(d)为本发明的去噪结果图。

本发明与现有的SK-LHRS滤波方法和PPB滤波方法的合成孔径雷达SAR图像相干斑抑制评价指标：滤波前后图像的均值比RoM和比值图的等效视数SL，垂直方向的边缘保持指数EPI(V)、水平方向的边缘保持指数EPI(H)以及对于图4(a)所示的区域A和区域B的均值保持和相干斑抑制指数MPSSI(A)和MPSSI(B)的值，如表1和表2所示。

表1SAR图像相干斑抑制评价指标对比表

表2均值保持和相干斑抑制指数对比表

	MPSSI(A)	MPSSI(B)
			SK-LHRS	0.0104	0.0062
PPB	0.0086	0.0062
			本发明方法	0.0076	0.0057
理想值	0.0000	0.0000

3.仿真结果分析：

从图3和图4可以看出，本发明相比SK-LHRS滤波方法和PPB滤波方法的优点在于，将基于几何结构块的非局部滤波方法和基于同质区域搜索的滤波方法有机的结合在一起，在对相干斑抑制的同时，有效地保持了SAR图像中结构信息。对比图3(b)和图3(d)，可以看出本发明相比于SK-LHRS滤波方法具有很好的细节保持特性。而对比图4(c)和图4(d)，可以看出PPB滤波方法不仅增强了合成孔径雷达SAR图像的细节信息，也在同质区域引入了一些人工痕迹，同时，该方法也把细节信息附近的噪声增强了出来。相比之下，我们的算法对于噪声的抑制和细节信息的保持具有很好的性能。

从表1和表2中可以看出，本发明算法在抑制斑点噪声的同时，能够很好地保持图像散射特性和几何结构信息的特性。相对于SK-LHRS滤波方法和PPB滤波方法，本发明在图像奇异信息保持和合成孔径雷达SAR图像相干斑抑制之间做到了较好的平衡。

综上所述，本发明实现了合成孔径雷达SAR图像斑点噪声抑制中结构信息保持与相干斑抑制之间的平衡，获得了合成孔径雷达SAR图像良好的去噪效果。

Claims (7)

1.一种基于几何先验和分布相似性测度的SAR图像去噪方法，包括如下步骤：

(1)素描化图像：

采用合成孔径雷达SAR素描模型，获得输入合成孔径雷达SAR图像的素描图；

(2)像素分类：

以素描图中素描线所包含的素描点的位置作为参考位置，将输入合成孔径雷达SAR图像中位置与参考位置相同的像素点作为具有方向特性的像素，将输入合成孔径雷达SAR图像中剩余的像素作为不具有方向特性的像素；

(3)采用基于几何结构块的非局部滤波方法，估计每个具有方向特性像素的像素值；

(4)估计不具有方向特性像素的像素值：

(4a)采用基于局部线性最小均方误差准则的估计方法，对于每个不具有方向特性的像素，在以该像素为中心的3×3个像素大小的邻域窗口内分别计算其初始估计值和初始等效视数；

(4b)从不具有方向特性的像素中选取任意一个未选像素，以该像素为中心，15个像素为边长构造正方形窗口；

(4c)采用区域生长方法，对正方形窗口中的所有像素进行生长和标记，将标记后的所有像素构成局部邻域，所述的区域生长方法中的相似性测度按照下式计算得到：

$Sim(E_p,E_q)={\[\frac{(2N_q-1)\·N_p}{(2N_p-1)\·N_q}\·r^2\]}^{N_p-0.5}\·{\[\frac{(N_p\·r^2+N_q)(2N_q-1)}{2N_q\·(N_p+N_q-1)}\]}^{1-N_p-N_q}$

其中，Sim(·)表示正方形窗口中种子点像素p的估计值与待生长像素q的估计值之间的相似性值，E_p和E_q分别表示正方形窗口中种子点像素p和待生长像素q的估计值，N_p和N_q分别表示正方形窗口中种子点像素p和待生长像素q的等效视数，r表示正方形窗口中种子点像素p的估计值E_p与待生长像素q的估计值E_q之间的比值，r＝E_p/E_q；

(4d)采用基于局部线性最小均方误差准则的估计方法，利用局部邻域内包含的所有像素，更新正方形窗口中心像素的估计值和等效视数；

(4e)采用相似性测度公式，利用正方形窗口中未标记像素的初始估计值和初始等效视数以及正方形窗口中心像素更新后的估计值和等效视数重新计算未标记像素与中心像素之间的相似性值，将相似性值满足区域生长阈值的像素加入到局部邻域并标记，得到扩充后的局部邻域；

(4f)采用极大似然准则，利用扩充后的局部邻域中的所有像素，估计正方形窗口中心像素的像素值；

(4g)判断所有不具有方向特性的像素中是否存在未选像素，若是，执行步骤(4b)，否则，执行步骤(5)；

(5)估计重叠像素的像素值：

(5a)采用下式表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数：

$p_s\left(s\right)=\frac{2L^L}{\mathrm{\Γ}\left(L\right)}{s}^{2L-1}\exp (-{\mathrm{Ls}}^2),s\≥0$

其中，p_s(·)表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数，s表示输入合成孔径雷达SAR图像的相干斑噪声，L表示输入合成孔径雷达SAR图像的视数，Γ(·)表示伽玛函数，exp(·)表示指数函数；

(5b)采用基于统计分布的加权融合方式，计算重叠像素的融合估计值：

$D=\frac{1}{Z}\underset{i}{\Σ}{p}_s\left(\frac{A}{D_i}\right)\·D_i,i\∈\{0,1\}$

其中，D表示重叠像素的融合估计值，Z表示加权融合的归一化因子，∑表示求和操作，i表示像素估计值的编号，A表示重叠像素在输入合成孔径雷达SAR图像中的灰度值，D₀表示重叠像素在步骤(3)得到的估计值，D₁表示重叠像素在步骤(4)得到的估计值，p_s(·)表示输入合成孔径雷达SAR图像中相干斑噪声的概率密度函数。

2.根据权利要求1所述的基于几何先验和分布相似性测度的SAR图像去噪方法，其特征在于：步骤(1)中所述的采用合成孔径雷达SAR素描模型，获得输入合成孔径雷达SAR图像的素描图的具体步骤如下：

第1步，构造具有不同方向和尺度的边、线模板，并利用模板的方向和尺度信息构造各向异性高斯函数来计算该模板中每一点的加权系数，其中尺度个数取值为3～5，方向个数取值为18；

第2步，按照下式，计算模板不同区域对应在合成孔径雷达SAR图像中像素的均值μ和方差ν：

$\mathrm{\μ}=\frac{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g{A}_g}{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g}$

$v=\frac{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g{(A_g-\mathrm{\μ})}^2}{\underset{g\∈\mathrm{\Ω}}{\Σ}{w}_g}$

其中，Ω表示模板中的某一区域，g表示区域Ω中点的位置，μ和ν分别表示区域Ω对应在合成孔径雷达SAR图像中像素的均值和方差，∑表示求和操作，w_g表示区域Ω中位置g处的权重系数，w_g的取值范围为0～1，A_g表示区域Ω中位置g对应在合成孔径雷达SAR图像中的像素值；

第3步，按照下式，计算合成孔径雷达SAR图像中每个像素对比值算子的响应值：

$R=1-min\{\frac{{\mathrm{\μ}}_a}{{\mathrm{\μ}}_b},\frac{{\mathrm{\μ}}_b}{{\mathrm{\μ}}_a}\}$

其中，R表示合成孔径雷达SAR图像中每个像素对比值算子的响应值，min{·}表示求最小值操作，a和b分别表示模板不同区域的编号，μ_a和μ_b分别表示区域a和区域b对应在合成孔径雷达SAR图像中像素的均值；

第4步，按照下式，计算合成孔径雷达SAR图像中每个像素对相关性算子的响应值：

$C=\sqrt{\frac{1}{1+2\·\frac{v_a^2+v_b^2}{{({\mathrm{\μ}}_a+{\mathrm{\μ}}_b)}^2}}}$

其中，C表示合成孔径雷达SAR图像中每个像素对相关性算子的响应值，a和b分别表示模板不同区域的编号，ν_a和ν_b分别表示区域a和区域b对应在合成孔径雷达SAR图像中像素的方差，μ_a和μ_b分别表示区域a和区域b对应在合成孔径雷达SAR图像中像素的均值，表示平方根操作；

第5步，按照下式，融合合成孔径雷达SAR图像中像素对比值算子的响应值和合成孔径雷达SAR图像中像素对相关性算子的响应值，计算合成孔径雷达SAR图像中每个像素对各个模板的响应值：

$F=\sqrt{\frac{R^2+C^2}{2}}$

其中，F表示合成孔径雷达SAR图像中每个像素对各个模板的响应值，R和C分别表示合成孔径雷达SAR图像中像素对比值算子和合成孔径雷达SAR图像中像素对相关性算子的响应值，表示平方根操作；

第6步，选择具有最大响应值的模板作为合成孔径雷达SAR图像中像素的模板，并将最大响应值作为该像素的强度，将具有最大响应值的模板的方向作为该像素的方向，获得合成孔径雷达SAR图像的边线响应图和方向图；

第7步，利用合成孔径雷达SAR图像中每个像素所选择的模板，获得合成孔径雷达SAR图像基于模板的梯度图；

第8步，按照下式，将归一化到[0,1]的边线响应图与归一化到[0,1]的基于模板的梯度图进行融合，得到最终的强度图：

$I=\frac{xy}{1-x-y+2xy}$

其中，I表示强度图中的强度值，x表示边线响应图中的值，y表示梯度图中的值；

第9步，利用方向图，采用非极大值抑制方法对强度图进行检测，得到建议草图；

第10步，选取建议草图中具有最大强度的像素，将建议草图中与该像素连通的像素连接形成建议线段，得到建议素描图；

第11步，利用合成孔径雷达SAR图像的统计分布，计算建议素描图中素描线的编码长度增益CLG：

$CLG=\underset{t}{\overset{m}{\Σ}}\[\frac{A_t^2}{A_{t,0}^2}+ln\left(A_{t,0}^2\right)-\frac{A_t^2}{A_{t,1}^2}-ln\left(A_{t,1}^2\right)\]$

其中，CLG表示建议素描图中素描线的编码长度增益，∑表示求和操作，m表示当前素描线邻域中像素的个数，t表示当前素描线邻域中像素的编号，A_t表示当前素描线邻域中第t个像素的观测值，A_t,0表示在当前素描线不能表示结构信息的假设下，该素描线邻域中第t个像素的估计值，ln(·)表示对数函数，A_t,1表示在当前素描线能够表示结构信息的假设下，该素描线邻域中第t个像素的估计值；

第12步，设定阈值T，T的取值范围为5～50，选择CLG>T的建议素描线作为最终素描图中的素描线，获得输入合成孔径雷达SAR图像对应的素描图。

3.根据权利要求1所述的基于几何先验和分布相似性测度的SAR图像去噪方法，其特征在于：步骤(3)所述的基于几何结构块的非局部滤波方法，计算所有具有方向特性像素的估计值的具体步骤如下：

第1步，将素描图中的每一条线段作为处理单位，以线段上的每一点为中心，沿着该线段的方向构建目标矩形块；

第2步，在具有方向特性的像素所组成的区域中，搜索与目标矩形块方向相同的相似矩形块；

第3步，按照下式，计算目标矩形块中所有像素的估计值：

${B_0}^{\′}=\underset{j}{\Σ}{w}_{0,j}{B}_j$

其中，B′₀表示目标矩形块B₀中所有像素的估计值，j表示相似矩形块的编号，∑表示求和操作，w_0,j表示目标矩形块B₀与相似矩形块B_j之间的相似性权重，B_j表示与目标矩形块B₀方向相同的相似矩形块。

4.根据权利要求1所述的基于几何先验和分布相似性测度的SAR图像去噪方法，其特征在于：步骤(4a)、步骤(4d)所述基于局部线性最小均方误差准则的估计方法的具体步骤如下：

第1步，按照下式，计算待估计像素的估计值：

$E=\frac{1}{n}\underset{k}{\overset{n}{\Σ}}{A}_k+\mathrm{\ξ}(A-\frac{1}{n}\underset{k}{\overset{n}{\Σ}}{A}_k)$

其中，E表示待估计像素的估计值，n表示邻域中所有像素的个数，k表示邻域中像素的编号，∑表示求和操作，A_k表示邻域中像素k在输入合成孔径雷达SAR图像中的灰度值，ξ表示衡量邻域内异质性程度的系数，ξ的取值范围为0～1，A表示待估计像素在输入合成孔径雷达SAR图像中的灰度值；

第2步，按照下式，计算等效视数的估计值：

$N=\frac{n}{(n-1){\mathrm{\ξ}}^2+1}\·L$

其中，N表示等效视数的估计值，n表示邻域中所有像素的个数，ξ表示衡量邻域内异质性程度的系数，ξ的取值范围为0～1，L表示输入合成孔径雷达SAR图像的视数。

5.根据权利要求1所述的基于几何先验和分布相似性测度的SAR图像去噪方法，其特征在于：步骤(4c)所述采用区域生长方法，得到局部邻域的具体步骤如下：

第1步，选取正方形窗口的中心像素作为种子点像素；

第2步，设定区域生长的阈值为0.1；

第3步，选取种子点像素邻域内的任一像素作为待生长像素；

第4步，采用相似性测度公式，计算正方形窗口中种子点像素的初始估计值与待生长像素的初始估计值之间的相似性值；

第5步，判断正方形窗口中种子点像素的初始估计值与待生长像素的初始估计值之间的相似性值是否满足区域生长阈值，若满足，则将待生长像素加入种子点像素所在的区域并标记；

第6步，在已加入种子点像素所在区域的像素的邻域中，寻找是否存在包含在正方形窗口内但未访问的像素，若存在，将该像素作为下一个待生长像素，执行第4步；否则，执行第7步；

第7步，将标记后的所有像素构成局部邻域。

6.根据权利要求1所述的基于几何先验和分布相似性测度的SAR图像去噪方法，其特征在于：步骤(4e)所述的未标记像素是指在区域生长过程中曾选为待生长像素但不包含在局部邻域内的像素。

7.根据权利要求1所述的基于几何先验和分布相似性测度的SAR图像去噪方法，其特征在于：步骤(5)所述的重叠像素是指包含在目标矩形块内的不具有方向特性的像素。