CN103778613B - 一种窗口自适应的极化sar影像滤波方法 - Google Patents

Abstract

本发明公开了一种窗口自适应的极化SAR影像滤波方法，本发明利用大小自适应和形状自适应的窗口模板对图像进行滤波：在同质区域用大的方形窗口进行滤波，在异质区域用小的非方形窗口进行滤波，从而实现在抑制噪声的同时保持点、线、边缘等细节信息，最大程度上提高相干斑噪声抑制效果。

Description

一种窗口自适应的极化SAR影像滤波方法

技术领域

本发明属于遥感图像处理技术领域，是一种图像相干斑噪声抑制方法，涉及一种对全极化合成孔径雷达（PolarimetricSyntheticApertureRadar，PolSAR）影像中的相干斑噪声进行滤除的新方法，具体涉及一种窗口自适应的极化SAR影像滤波方法。

背景技术

全极化SAR***相对于单极化SAR可以提供更加丰富的全极化信息，因此正在逐渐代替单极化SAR在地表分类、目标识别、变化监测、参数反演等方面扮演越来越重要的角色。但是由于SAR***采用相干成像，SAR图像中不可避免的存在相干斑噪声。相干斑噪声降低了数据的信噪比，令SAR图像解译、参数反演等工作变得更加困难。为了减小相干斑噪声的影响，通常在进一步处理之前首先对SAR图像进行相干斑噪声滤波，以尽可能的抑制噪声，同时尽量保持图像的细节信息。因此，自从SAR***产生以来，SAR图像滤波算法的研究就从未停止过，并且到目前为止仍然是一个热点。目前，针对极化SAR数据进行相干斑噪声滤波的方法已经有很多，其中最为经典的是精制极化Lee滤波算法，该方法采用非方形的边缘方向窗口达到保持边缘信息的目的，并且对协方差矩阵的所有元素利用相同的参数进行滤波，以防止通道间串扰。这种防止通道间串扰的思想对极化SAR滤波算法的研究影响深远，成为极化SAR滤波的一个基本原则。然而，此滤波算法存在明显的缺陷：在同质区域有明显但是扇贝效应，在较暗的线状区域有明显的虚假亮线。因此对于有些应用如道路提取、图像分割等会造成较大影响。

针对以上缺陷，目前已经有一些改进算法。这些算法的改进思路主要有两个：1）增加同质像素选择步骤；2）令滤波窗口大小自适应。第一种思路在非方形窗口中进一步筛选同质像素，会降低噪声抑制能力，增加算法复杂度；第二种思路仅能改善虚假细线问题，但不能从根本上克服这两种缺陷。

发明内容

本发明的目的在于从根本上克服精制极化Lee滤波存在的两种缺陷，提出了一种窗口自适应的极化SAR影像滤波方法。

本发明所采用的技术方案是：一种窗口自适应的极化SAR影像滤波方法，其特征在于，包括以下步骤：

步骤1：设置参数：设置最大滤波窗口Wmax和最小滤波窗口Wmin，获得N=(Wmax–Wmin)/2+1个不同大小的窗口W；

步骤2：计算极化同质度图：对不同大小的窗口W,计算对应的极化同质度图PolH(W)，同时得到像素方向图θ(W)和像素类型图T_p(W)；

步骤3：自动检测同质度阈值：对不同窗口下的极化同质度图PolH(W)利用阈值自动检测算法计算同质度阈值T(W)；

步骤4：自适应窗口滤波，其具体实现包括以下子步骤：

步骤4.1：对当前像素P，设置当前窗口大小W=W_max；

步骤4.2：判断：

如果PolH(W)>T(W)，则选择大小为W的方形窗口，并转到步骤4.5；

否则，转到步骤4.3；

步骤4.3：判断：

如果W=W_min，则转到步骤4.4；

否则，令W=W–2，回转执行所述的步骤4.2；

步骤4.4：判断：

如果像素类型T_p(W)是线，则选择线性窗口；

否则，如果像素类型T_p(W)是边缘，则选择边缘窗口；

步骤4.5：用选择的滤波窗口对中心像素进行局部线性最小均方差LLMMSE滤波，移到下个像素，并回转执行所述的步骤4.1，直到所有像素都滤波完毕。

作为优选，步骤2中所述的对不同大小的窗口W,计算对应的极化同质度图PolH(W)，其极化同质度图PolH的计算公式为：PolH=ENL/LAE，其中，ENL是等效视数图，LAE是线和边缘检测图。

作为优选，所述的等效视数图ENL的计算公式为：

$\mathrm{ENL}=\frac{\mathrm{tr}{\left(\mathrm{\&Sigma;}\right)}^2}{<\mathrm{tr}\left(\mathrm{ZZ}\right)>-\mathrm{tr}\left(\mathrm{\&Sigma;\&Sigma;}\right)}$

其中，Z=X/L，X表示极化SAR图像相干矩阵或者协方差矩阵，L表示极化SAR图像的视数；Σ是当前估计窗口下Z的均值，tr(·)表示求矩阵的迹。

作为优选，所述的线和边缘检测图LAE的计算方法包括以下子步骤：

步骤2.1：设置模板参数l_f，w_f，d_f，和θ_f，获得N_f=π/θ_f个模板；其中，l_f是模板长度，w_f是模板宽度，d_f是区域R₁和R₂间的距离，θ_f是模板对应的角度；

步骤2.2：对N_f个不同方向θ下的模板，计算区域R₁和R₂之间的梯度D₁₂(θ)、R₁和R₃之间的梯度D₁₃(θ)，以及R₂和R₃之间的梯度D₂₃(θ)，获得3*N_f个梯度值，其中，梯度D₁₂(θ)的计算公式为：

D₁₂(θ)=2ln|X₁+X₂|-ln|X₁|-ln|X₂|+2qln2

其中，X₁、X₂分别表示方向θ下模板区域R₁和区域R₂的q*q的极化SAR图像相干矩阵或协方差矩阵；其余梯度的计算方法与梯度D₁₂(θ)的计算方法相同；

步骤2.3：求出3*N_f个梯度值中的最大值D_max及其对应的方向θ_max，判断：

如果D_max=D₁₃(θ)或者D_max=D₂₃(θ)，则转到步骤2.4；

否则，转到步骤2.5；

步骤2.4：当θ≤180°，θ=θ；当θ≥180°，θ=θ–180°；计算θ=θ_max时中心像素P₀和区域R₃之间的非相似度D₀₃(θ_max)，以及θ=θ_max+90°时中心像素P₀和区域R₃之间的非相似度D₀₃(θ_max+90°)；判断：如果D₀₃(θ_max)<D₀₃(θ_max+90°)，则最终的选择的模板方向为θ_max；否则，最终选择的模板方向为θ_max+90°；最终选择的模板对应的梯度为新的D_max，其中，非相似度为中心像素P₀和区域R₃在Span图像上的差值；

步骤2.5：获得最终的梯度D_max，即为LAE值，同时获得模板方向θ和目标类型T_p(W)，其中，如果D_max=D₁₂，则是T_p(W)边缘；否则，T_p(W)是线。

作为优选，步骤3中所述的对不同窗口下的极化同质度图PolH(W)利用阈值自动检测算法计算同质度阈值T(W)，其阈值自动检测算法为熵阈值自动检测算法，具体实现包括以下子步骤：

步骤3.1：图像增强；首先，对极化同质度图PolH进行对数变换，计算公式为：x=10*log₁₀(x)，再进行最大最小值拉伸，处理方法为：将图像像素值按从小到大的顺序排序，取第6.25%位置上的像素值为最小值min，取第93.75%位置上的像素值为最大值max，然后利用最大值最小值进行归一化处理，归一化公式为：x=(x-min)/(max-min)；

步骤3.2：将PolH图中的像素分为n个灰度级，计算每个灰度级i中像素占总像素数的比率p_i；

步骤3.3：假设灰度级s将图像分为A和B两部分，每个部分的熵定义如下：

$H_A=-\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}\frac{p_i}{P_s}\ln \frac{p_i}{P_s}$

$H_B=-\underset{i=i+s}{\overset{n}{\mathrm{\&Sigma;}}}\frac{p_i}{1-P_s}\ln \frac{p_i}{1-P_s}$

其中， $P_s=\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{p}_i,$

则A和B的熵的和为：H_sum=H_A+H_B，计算每个灰度级i对应的H_sum；

步骤3.4：求出最大的H_sum及其对应的灰度级s，则s就是最终的分割阈值。

作为优选，步骤1中所述的用选择的滤波窗口对中心像素进行局部线性最小均方差LLMMSE滤波，其中LLMMSE滤波的公式如下：

$\stackrel{\&OverBar;}{x}=\stackrel{\&OverBar;}{y}+b(y-\stackrel{\&OverBar;}{y})$

其中，表示滤波后的像素值，表示局部均值，y是中心像素值，b是权值参数，其计算公式为：

$b=\frac{\mathrm{var}\left(y\right)-{\stackrel{\&OverBar;}{y}}^2{{\mathrm{\&sigma;}}_v}^2}{(1+{{\mathrm{\&sigma;}}_v}^2)\mathrm{var}\left(y\right)}$

其中，var(y)表示局部方差，σ_v ²表示噪声方差。

本发明利用大小自适应和形状自适应的窗口模板对图像进行滤波：在同质区域用大的方形窗口进行滤波，在异质区域用小的非方形窗口进行滤波，从而实现在抑制噪声的同时保持点、线、边缘等细节信息，最大程度上提高相干斑噪声抑制效果。

本发明的创新之处在于：

1）提出新的线和边缘检测方法，能更准确的检测线状区域的方向和梯度；

2）提出极化同质度图，利用等效视数图和边缘检测图之比增强了同质区域和异质区域间的对比度，从而可以通过自动阈值检测算法确定同质区域和异质区域间的阈值；

3）利用极化同质度图将同质区域和异质区域分别处理，在同质区域用大的方形窗口滤波，从而最大程度上抑制相干斑噪声，并且避免了精制极化Lee滤波扇贝效应问题；在异质区域用小的非方形窗口滤波，从而最大程度上保持细节信息，并且克服了精制极化Lee滤波中虚假亮线的问题。

附图说明

图1：是本发明的流程图。

图2：是本发明实施例中的边缘检测模板图。

图3：是本发明实施例中以W=5为例的边缘窗口模板图。

图4：是本发明实施例中以W=3为例的梯度算子图。

图5：是本发明实施例中以W=3为例的非方形窗口图。

具体实施方式

以下将结合附图和具体实施例对本发明做进一步的阐述。

请见图1，本发明所采用的技术方案是：一种窗口自适应的极化SAR影像滤波方法，包括以下步骤：

步骤1：设置参数：设置最大滤波窗口W_max和最小滤波窗口W_min，获得N=(W_max–W_min)/2+1个不同大小的窗口W；

步骤2：计算极化同质度图：对不同大小的窗口W，计算对应的极化同质度图PolH(W)，同时得到像素方向图θ(W)和像素类型图T_p(W)；

其极化同质度图PolH的计算公式为：PolH=ENL/LAE，其中，ENL是等效视数图，LAE是线和边缘检测图。

等效视数图ENL的计算公式为：

$\mathrm{ENL}=\frac{\mathrm{tr}{\left(\mathrm{\&Sigma;}\right)}^2}{<\mathrm{tr}\left(\mathrm{ZZ}\right)>-\mathrm{tr}\left(\mathrm{\&Sigma;\&Sigma;}\right)}$

其中，Z=X/L，X表示极化SAR图像相干矩阵或者协方差矩阵，L表示极化SAR图像的视数；Σ是当前估计窗口下Z的均值，tr(·)表示求矩阵的迹。

请见图2，线和边缘检测图LAE的计算方法包括以下子步骤：

步骤2.1：设置模板参数l_f，w_f，d_f，和θ_f，获得N_f=π/θ_f个模板；其中，l_f是模板长度，w_f是模板宽度，d_f是区域R₁和R₂间的距离，θ_f是模板对应的角度；一般情况下模板为正方形，即l_f=2*w_f+d_f；本实施例设置l_f=3，w_f=1，d_f=1，θ_f=45°，则共有4个不同方向的模板。

步骤2.2：对4个不同方向θ下的模板，计算区域R₁和R₂之间的梯度D₁₂(θ)、R₁和R₃之间的梯度D₁₃(θ)，以及R₂和R₃之间的梯度D₂₃(θ)，获得12个梯度值，请见图4，为对应的梯度算子，其中，a～d为边缘梯度，e～l为线性梯度。梯度D₁₂(θ)的计算公式为：

D₁₂(θ)=2ln|X₁+X₂|-ln|X₁|-ln|X₂|+2qln2

其中，X₁、X₂分别表示方向θ下模板区域R₁和区域R₂的q*q的极化SAR图像相干矩阵或协方差矩阵；其余梯度的计算方法与梯度D₁₂(θ)的计算方法相同；

步骤2.3：求出12个梯度值中的最大值D_max及其对应的方向θ_max，判断：

如果D_max=D₁₃(θ)或者D_max=D₂₃(θ)，则转到步骤2.4；

否则，转到步骤2.5；

步骤2.4：当θ≤180°，θ=θ；当θ≥180°，θ=θ–180°；计算θ=θ_max时中心像素P₀和区域R₃之间的非相似度D₀₃(θ_max)，以及θ=θ_max+90°时中心像素P₀和区域R₃之间的非相似度D₀₃(θ_max+90°)；判断：如果D₀₃(θ_max)<D₀₃(θ_max+90°)，则最终的选择的模板方向为θ_max；否则，最终选择的模板方向为θ_max+90°；最终选择的模板对应的梯度为新的D_max，其中，非相似度为中心像素P₀和区域R₃在Span图像上的差值；

步骤2.5：获得最终的梯度D_max，即为LAE值，同时获得模板方向θ和目标类型T_p(W)，其中，如果D_max=D₁₂，则是T_p(W)边缘；否则，T_p(W)是线。

以上步骤在计算LAE的过程中会同时得到像素方向图θ(W)和像素类型图T_p(W)，这两个图会在后续的步骤中用到。

步骤3：自动检测同质度阈值：对不同窗口下的极化同质度图PolH(W)利用阈值自动检测算法计算同质度阈值T(W)；

本实施例中，其阈值自动检测算法为熵阈值自动检测算法，具体实现包括以下子步骤：

步骤3.1：图像增强；首先，对极化同质度图PolH进行对数变换，计算公式为：x=10*log₁₀(x)，再进行最大最小值拉伸，处理方法为：将图像像素值按从小到大的顺序排序，取第6.25%位置上的像素值为最小值min，取第93.75%位置上的像素值为最大值max，然后利用最大值最小值进行归一化处理，归一化公式为：x=(x-min)/(max-min)；

步骤3.2：将PolH图中的像素分为n个灰度级，计算每个灰度级i中像素占总像素数的比率p_i；本实施例中n=256；

步骤3.3：假设灰度级s将图像分为A和B两部分，每个部分的熵定义如下：

$H_A=-\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}\frac{p_i}{P_s}\ln \frac{p_i}{P_s}$

$H_B=-\underset{i=i+s}{\overset{n}{\mathrm{\&Sigma;}}}\frac{p_i}{1-P_s}\ln \frac{p_i}{1-P_s}$

其中， $P_s=\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}{p}_i,$

则A和B的熵的和为：H_sum=H_A+H_B，计算每个灰度级i对应的H_sum；

步骤3.4：求出最大的H_sum及其对应的灰度级s，则s就是最终的分割阈值。

步骤4：自适应窗口滤波，请见图2、图3和图5，其具体实现包括以下子步骤：

步骤4.1：对当前像素P，设置当前窗口大小W=W_max；

步骤4.2：判断：

如果PolH(W)>T(W)，则选择大小为W的方形窗口，并转到步骤4.5；

否则，转到步骤4.3；

步骤4.3：判断：

如果W=W_min，则转到步骤4.4；

否则，令W=W–2，回转执行所述的步骤4.2；

步骤4.4：判断：

如果像素类型T_p(W)是线，则选择线性窗口，请见图2，即R₃区域，其中l_f=W，θ_f=θ(W)；请见图5中，其中，1～8为边缘窗口，9～12为线性窗口，本实施例以W=3为例，其为9～12号窗口；

否则，如果像素类型T_p(W)是边缘，则选择边缘窗口，请见图3，即白色区域，请见图5中，其中，1～8为边缘窗口，9～12为线性窗口，本实施例以W=3为例，其为1～8号窗口；

由于每个角度对应两个边缘窗口，因此需进一步计算，这两个窗口哪个与中心像素在Span图上的差值小，就选择哪个。

步骤4.5：用选择的滤波窗口对中心像素进行局部线性最小均方差LLMMSE滤波，移到下个像素，并回转执行所述的步骤4.1，直到所有像素都滤波完毕。

其中LLMMSE滤波的公式如下：

$\stackrel{\&OverBar;}{x}=\stackrel{\&OverBar;}{y}+b(y-\stackrel{\&OverBar;}{y})$

其中，表示滤波后的像素值，表示局部均值，y是中心像素值，b是权值参数，其计算公式为

$b=\frac{\mathrm{var}\left(y\right)-{\stackrel{\&OverBar;}{y}}^2{{\mathrm{\&sigma;}}_v}^2}{(1+{{\mathrm{\&sigma;}}_v}^2)\mathrm{var}\left(y\right)}$

其中，var(y)表示局部方差，σ_v ²表示噪声方差。

本实施例中，首先用Span图像计算权重b，然后对相干矩阵或协方差矩阵的所有元素用相同的权重及滤波窗口进行滤波。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代。如模板参数θ_f可设置为其它值；梯度及非相似度可采用其它公式代替；阈值自动检测算法可以选择其它方法代替熵阈值法；等等。但这并不会超出本发明所提算法的框架，不会偏离本发明的精神，或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种窗口自适应的极化SAR影像滤波方法，其特征在于，包括以下步骤：

步骤1：设置参数：设置最大滤波窗口W_max和最小滤波窗口W_min，获得N＝(W_max–W_min)/2+1个不同大小的窗口W；

步骤2：计算极化同质度图：对不同大小的窗口W，计算对应的极化同质度图PolH(W)，同时得到像素方向图θ(W)和像素类型图T_p(W)；

其中所述的对不同大小的窗口W，计算对应的极化同质度图PolH(W)，其极化同质度图PolH的计算公式为：PolH＝ENL/LAE，其中，ENL是等效视数图，LAE是线和边缘检测图；所述的等效视数图ENL的计算公式为：

$ENL=\frac{tr{\left(\mathrm{\&Sigma;}\right)}^2}{<tr\left(ZZ\right)>-tr\left(\mathrm{\&Sigma;}\mathrm{\&Sigma;}\right)};$

其中，Z＝X/L，X表示极化SAR图像相干矩阵或者协方差矩阵，L表示极化SAR图像的视数；Σ是当前估计窗口下Z的均值，tr(·)表示求矩阵的迹；

步骤3：自动检测同质度阈值：对不同窗口下的极化同质度图PolH(W)利用阈值自动检测算法计算同质度阈值T(W)；

步骤4：自适应窗口滤波，其具体实现包括以下子步骤：

步骤4.1：对当前像素P，设置当前窗口大小W＝W_max；

步骤4.2：判断：

如果PolH(W)>T(W)，则选择大小为W的方形窗口，并转到步骤4.5；

否则，转到步骤4.3；

步骤4.3：判断：

如果W＝W_min，则转到步骤4.4；

否则，令W＝W–2，回转执行所述的步骤4.2；

步骤4.4：判断：

如果像素类型T_p(W)是线，则选择线性窗口；

否则，如果像素类型T_p(W)是边缘，则选择边缘窗口；

步骤4.5：用选择的滤波窗口对中心像素进行局部线性最小均方差LLMMSE滤波，移到下个像素，并回转执行所述的步骤4.1，直到所有像素都滤波完毕。

2.根据权利要求1所述的窗口自适应的极化SAR影像滤波方法，其特征在于：所述的线和边缘检测图LAE的计算方法包括以下子步骤：

步骤2.1：设置模板参数l_f，w_f，d_f，和θ_f，获得N_f＝π/θ_f个模板；其中，l_f是模板长度，w_f是模板宽度，d_f是区域R₁和R₂间的距离，θ_f是模板对应的角度；

步骤2.2：对N_f个不同方向θ下的模板，计算区域R₁和R₂之间的梯度D₁₂(θ)、R₁和R₃之间的梯度D₁₃(θ)，以及R₂和R₃之间的梯度D₂₃(θ)，获得3*N_f个梯度值，其中，梯度D₁₂(θ)的计算公式为：

D₁₂(θ)＝2ln|X₁+X₂|-ln|X₁|-ln|X₂|+2qln2

其中，X₁、X₂分别表示方向θ下模板区域R₁和区域R₂的q*q的极化SAR图像相干矩阵或协方差矩阵；其余梯度的计算方法与梯度D₁₂(θ)的计算方法相同；

步骤2.3：求出3*N_f个梯度值中的最大值D_max及其对应的方向θ_max，判断：

如果D_max＝D₁₃(θ)或者D_max＝D₂₃(θ)，则转到步骤2.4；

否则，转到步骤2.5；

步骤2.4：当θ≤180°，θ＝θ；当θ≥180°，θ＝θ–180°；计算θ＝θ_max时中心像素P₀和区域R₃之间的非相似度D₀₃(θ_max)，以及θ＝θ_max+90°时中心像素P₀和区域R₃之间的非相似度D₀₃(θ_max+90°)；判断：如果D₀₃(θ_max)<D₀₃(θ_max+90°)，则最终的选择的模板方向为θ_max；否则，最终选择的模板方向为θ_max+90°；最终选择的模板对应的梯度为新的D_max，其中，非相似度为中心像素P₀和区域R₃在Span图像上的差值；

步骤2.5：获得最终的梯度D_max，即为LAE值，同时获得模板方向θ和目标类型T_p(W)，其中，如果D_max＝D₁₂，则T_p(W)是边缘；否则，T_p(W)是线。

3.根据权利要求1所述的窗口自适应的极化SAR影像滤波方法，其特征在于：步骤3中所述的对不同窗口下的极化同质度图PolH(W)利用阈值自动检测算法计算同质度阈值T(W)，其阈值自动检测算法为熵阈值自动检测算法，具体实现包括以下子步骤：

步骤3.1：图像增强；首先，对极化同质度图PolH进行对数变换，计算公式为：x＝10*log₁₀(x)，再进行最大最小值拉伸，处理方法为：将图像像素值按从小到大的顺序排序，取第6.25％位置上的像素值为最小值min，取第93.75％位置上的像素值为最大值max，然后利用最大值最小值进行归一化处理，归一化公式为：x＝(x-min)/(max-min)；

步骤3.2：将PolH图中的像素分为n个灰度级，计算每个灰度级i中像素占总像素数的比率p_i；

步骤3.3：假设灰度级s将图像分为A和B两部分，每个部分的熵定义如下：

$H_A=-\underset{i=1}{\overset{s}{\mathrm{\&Sigma;}}}\frac{p_i}{P_s}ln\frac{p_i}{P_s}$

$H_B=-\underset{i=i+s}{\overset{n}{\mathrm{\&Sigma;}}}\frac{p_i}{1-P_s}ln\frac{p_i}{1-P_s}$

其中，

则A和B的熵的和为：H_sum＝H_A+H_B，计算每个灰度级i对应的H_sum；

步骤3.4：求出最大的H_sum及其对应的灰度级s，则s就是最终的分割阈值。

4.根据权利要求1所述的窗口自适应的极化SAR影像滤波方法，其特征在于：步骤4.5中所述的用选择的滤波窗口对中心像素进行局部线性最小均方差LLMMSE滤波，其中LLMMSE滤波的公式如下：

$\hat{x}=\stackrel{\&OverBar;}{y}+b(y-\stackrel{\&OverBar;}{y})$

其中，表示滤波后的像素值，表示局部均值，y是中心像素值，b是权值参数，其计算公式为：

$b=\frac{\mathrm{var}\left(y\right)-{\stackrel{\&OverBar;}{y}}^2{{\mathrm{\&sigma;}}_v}^2}{(1+{{\mathrm{\&sigma;}}_v}^2)\mathrm{var}\left(y\right)}$

其中，var(y)表示局部方差，σ_v ²表示噪声方差。