CN103258207A - 基于散射功率和强度联合统计的极化sar图像分类方法 - Google Patents

Abstract

本发明公开了一种基于散射功率和强度联合统计的极化SAR图像分类方法，主要解决现有技术对散射特性相似区域难以区分和分类类别数固定的问题。其实现过程是：1）采用Lee滤波器对相干矩阵T滤波；2)利用Freeman分解得到功率矩阵；3）利用特征值分解得到强度矩阵；4)对功率矩阵和强度矩阵分别进行8邻域求平均；5)选取k类匀质区域作为训练样本；6）采用EM算法对k类样本的功率矩阵和强度矩阵概率密度分布函数的参数进行估计；7)求解k类样本功率矩阵和强度矩阵的联合概率分布；8)对待分类极化SAR数据进行贝叶斯分类，得到分类结果。本发明具有对极化SAR图像分类效果显著的优点，可进一步用于极化SAR图像目标检测和目标识别。

Description

基于散射功率和强度联合统计的极化SAR图像分类方法

技术领域

本发明属于图像处理技术领域，特别是涉及极化SAR图像处理，可用于雷达图像中目标检测和目标分类与识别。

背景技术

随着信息处理技术和电子技术的发展，合成孔径雷达开始朝着高分辨、多极化、多工作模式、多波段等方向发展，致力于为人们提供更多更丰富的目标散射信息。其中多极化方向的发展是一个很重要的方向。相比于传统的单极化合成孔径雷达仅能在特定电磁波极化收发组合下的目标散射特性探测能力，多极化的合成孔径雷达无疑为目标的解译和正确的认识提供了更多更全面的信息量。更加丰富的信息量为深入地研究目标的散射特性提供了重要的依据，进而对目标的自动检测和自动识别有着非常重要的指导作用。现有的极化SAR图像分类方法大致可以分为基于散射特性的方法和基于统计特性方法，如基于Freeman分解的分类方法和复Wishart分类方法等。

其中，基于统计特性的极化SAR图像分类方法，是将极化SAR图像中具有某方面相似特征的像素点划分为一类，已经有很多成熟的分类算法被用到极化SAR图像分类中。作为一门新兴学科,目前主要的研究工作集中在散射特征的提取和适合极化SAR图像分类的分类器设计上。不同极化散射特征的提取方法得到不同的特征，不同的特征又适合于不同性能的分类器。现有的特征提取方法主要有Freeman分解、Cloude分解方法、SDH分解和Pauli基分解等；典型的分类器有复Wishart分类器和H/alpha分类器等。其中，利用Freeman分解方法得到特征对极化SAR图像进行分类时，存在对特征分布特性认知上的不足，对于散射机制较复杂的类别，难以很好地进行识别与分类；利用H/alpha分类器进行类别判定时，存在类别判决界限需要人为确定的不足，而且H/alpha和复Wishart都是基于极化SAR散射机制实现的分类，对于一些散射特性相似但不同类别的地物，该方法表现出了分类与识别混淆的缺陷。这些不足限制了其在极化SAR图像分类上的广泛应用。

发明内容

本发明的目的在于克服上述已有技术的不足，针对极化SAR分类中混合散射机制区域分类难、散射特性相似但不同类情况下易造成错分等问题，提出了一种基于散射功率和强度联合统计的极化SAR分类方法，以提高极化SAR图像分类的精度和分类方法的普适性。

为实现上述目的，本发明提出了基于散射功率和强度联合统计的极化SAR分类方法，包括如下步骤：

(1)将待分类极化SAR图像的相干矩阵T作为输入，T的大小为3×3×z，其中z是图像像素点总个数；

(2)采用Lee滤波器对相干矩阵T进行去噪，得到去噪后的相干矩阵T'；

(3)对去噪后的相干矩阵T'进行Freeman分解，得到体散射功率，取每个像素点的体散射功率构成大小为z×1的功率矩阵P_v；

(4)对去噪后的相干矩阵T'进行特征值分解，提取极化SAR图像的散射强度信息，得到每个像素点的三个特征值λ₁,λ₂,λ₃，取最大特征值λ₁构成大小为z×1的强度矩阵E；

(5)对功率矩阵P_v和强度矩阵E分别进行大小为m×n的变形，得到变形后的功率矩阵P_vT和强度矩阵E_T，m,n分别是待分类图像的行数与列数，对变形后功率矩阵P_vT和变形后的强度矩阵E_T的每个像素点分别进行8邻域求平均来代替该像素点的值，得到平滑后的功率矩阵P'_v和平滑后的强度矩阵E'；

(6)依据实际地物分布，选择不同的k类匀质区域作为训练样本，k=8；

(7)采用EM算法分别对k类样本的平滑后功率矩阵P'_v和平滑后强度矩阵E'的高斯参数进行估计，得到k类样本的平滑后功率矩阵P'_v的概率密度分布函数P(x|μ_vi,σ_vi)和平滑后强度矩阵E'的概率密度分布函数P(y|μ_λi,σ_λi)：

$P\left(x\right|{\mathrm{\μ}}_{\mathrm{vi}},{\mathrm{\σ}}_{\mathrm{vi}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{vi}}}\exp (-\frac{{(x-{\mathrm{\μ}}_{\mathrm{vi}})}^2}{2{\mathrm{\σ}}_{\mathrm{vi}}^2}),$

$P\left(y\right|{\mathrm{\μ}}_{\mathrm{\λi}},{\mathrm{\σ}}_{\mathrm{\λi}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{\λi}}}\exp (-\frac{{(y-{\mathrm{\μ}}_{\mathrm{\λi}})}^2}{2{\mathrm{\σ}}_{\mathrm{\λi}}^2}),$

其中x为平滑后功率矩阵P'_v的值，y为平滑后强度矩阵E'的值、μ_vi为第i类样本的平滑后功率矩阵P'_v的均值，σ_vi为第i类样本的平滑后功率矩阵P'_v的方差，μ_λi为第i类样本的平滑后强度矩阵E'的均值，σ_λi为第i类样本的平滑后强度矩阵E'的方差，i=1,2,...,k；

(8)利用步骤(7)求得的k类样本的概率密度分布函数，构造平滑后功率矩阵P'_v和平滑后强度矩阵E'的联合概率密度分布函数P(a_j|a_j∈C_i)：

P(a_j|a_j∈C_i)=P(x|μ_vi,σ_vi)·P(y|μ_λi,σ_λi)，

其中a_j为第j个像素点，j=1,2,...,z,C_i表示第i类样本集合,i=1,2,...,k；

(9)根据MAP原则，依据表达式：a_j∈C_i,st.p(a_j|C_i)=maxP(a_j|a_j∈C_i)，i=1,2,...,k，j=1,2,...,z，对待分类极化SAR数据进行贝叶斯分类，并输出分类结果。

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

1.本发明相对已有的极化SAR图像分类技术，从统计的角度入手，利用散射功率和散射强度的联合统计实现散射特性相似但不同类目标区域的分类，克服了H/alpha和复Wishart算法在这方面的缺陷，从而提高了分类精度。

2.本发明的分类类别数由人为确定，克服了H/alpha和复Wishart算法中分类类别数固定的局限，提高了算法的普适性。

仿真实验结果表明，本发明提出的基于散射功率和强度联合统计的极化SAR图像分类方法，能够有效地应用于极化SAR图像分类，具有普适性，可进一步应用于雷达图像目标检测与分类识别。

附图说明

图1是本发明的实现流程图；

图2是本发明仿真采用的极化SAR图像的RGB合成图及在该图上标记的8类训练样本；

图3是本发明仿真图2中8类训练样本平滑后功率矩阵实际分布和估计分布示意图；

图4是本发明仿真图2中8类训练样本平滑后强度矩阵实际分布和估计分布示意图；

图5是采用本发明对图2进行分类的结果图。

具体实施方式

参照图1，本发明的具体实施如下：

步骤1.将待分类极化SAR图像的相干矩阵T作为输入，T的大小为3×3×z，其中z是图像像素点总个数。

步骤2.采用Lee滤波器对相干矩阵T进行去噪，得到去噪后的相干矩阵T'。

极化SAR数据含有大量的噪声，分类前需要对原始数据去噪以减少噪声对分类结果的影响。极化SAR数据的噪声属于乘加性噪声，一般的去噪方法对极化SAR去噪效果很不理想，这里采用Lee滤波器对相干矩阵T进行去噪，得到去噪后的相干矩阵T'。

步骤3.对去噪后的相干矩阵T'进行Freeman分解，提取极化SAR数据的体散射功率，构成大小为z×1的功率矩阵P_v。

首先，利用去噪后的相干矩阵T'得到体散射分量的分解系数F_v：

$F_v=\frac{3}{2}{T}^{\′};$

然后，通过利用Freeman分解公式，得到功率矩阵P_v：

P_v=8F_v/3。

步骤4.对去噪后的相干矩阵T'进行特征值分解，提取极化SAR的散射强度信息，得到每个像素点的三个特征值λ₁,λ₂,λ₃，取最大特征值λ₁构成大小为z×1的强度矩阵E。

由于地物的散射机制是由多种散射特性混合而成的，对散射特性不同的区域，利用散射功率P_v可以实现分类，但是，对散射特性相似但不同类别的区域，散射功率P_v是难以实现分类的，因此需要加入强度矩阵来解决这一问题。

利用MATLAB软件中的eigs函数对去噪后的相干矩阵T'进行特征值分解，分解表达式如下：

$\left[T^{\′}\right]=\left[U\right]\left[\begin{array}{ccc}{\mathrm{\λ}}_1& 0& 0\\ 0& {\mathrm{\λ}}_2& 0\\ 0& 0& {\mathrm{\λ}}_3\end{array}\right]{\left[U\right]}^{*};$

其中*为共轭转置符，U为去噪后相干矩阵T'特征值分解的特征向量，λ₁,λ₂,λ₃为去噪后相干矩阵T'特征值分解所得的大小不同的三个特征值，且λ₁>λ₂>λ₃，取最大特征值λ₁构成大小为z×1的强度矩阵E。

步骤5.对功率矩阵P_v和强度矩阵E分别进行大小为m×n的变形，得到变形后的功率矩阵P_vT和强度矩阵E_T，m,n分别是待分类图像的行数与列数，对变形后功率矩阵P_vT和变形后的强度矩阵E_T的每个像素点分别进行8邻域求平均来代替该像素点的值，得到平滑后的功率矩阵P'_v和平滑后的强度矩阵E'。

步骤6.依据地物实际分布，选择不同的k类匀质区域作为训练样本，本实例取k=8。

步骤7.采用EM算法分别对k类样本的平滑后功率矩阵P'_v和平滑后强度矩阵E'的高斯模型参数进行估计，得到k类样本的平滑后功率矩阵P'_v的概率密度分布函数P(x|μ_vi,σ_vi)和平滑后强度矩阵E'的概率密度分布函数P(y|μ_λi,σ_λi)：

$P\left(x\right|{\mathrm{\μ}}_{\mathrm{vi}},{\mathrm{\σ}}_{\mathrm{vi}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{vi}}}\exp (-\frac{{(x-{\mathrm{\μ}}_{\mathrm{vi}})}^2}{2{\mathrm{\σ}}_{\mathrm{vi}}^2}),$

$P\left(y\right|{\mathrm{\μ}}_{\mathrm{\λi}},{\mathrm{\σ}}_{\mathrm{\λi}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{\λi}}}\exp (-\frac{{(x-{\mathrm{\μ}}_{\mathrm{\λv}})}^2}{2{\mathrm{\σ}}_{\mathrm{\λi}}^2}),$

其中x为平滑后功率矩阵P'_v的值，y为平滑后强度矩阵E'的值、μ_vi为第i类样本的平滑后功率矩阵P'_v的均值，σ_vi为第i类样本的平滑后功率矩阵P'_v的方差，μ_λi为第i类样本的平滑后强度矩阵E'的均值，σ_λi为第i类样本的平滑后强度矩阵E'的方差，i=1,2,...,k；

7a)将各类样本的平滑后功率矩阵P'_v和平滑后强度矩阵E'分别作为EM算法的输入数据；

7b)设定EM算法迭代终止条件为迭代误差小于0.00001，分别将第i类样本的平滑后功率矩阵P'_v的均值μ_vi和方差σ_vi，第i类样本的平滑后强度矩阵E'的均值μ_λi和方差σ_λi初始化为随机值，其中i=1,2,...,k；

7c)通过EM算法中的期望最大化对步骤（7b）中所述的参数μ_vi,μ_λi,σ_vi和σ_λi，进行更新，每次更新后计算迭代误差并判断是否满足迭代终止条件，若满足迭代终止条件，则迭代停止，输出参数μ_vi,μ_λi,σ_vi和σ_λi，否则继续执行步骤（7c）；

7d)将输出的参数μ_vi,μ_λi,σ_vi和σ_λi代入平滑后功率矩阵P'_v的概率密度分布函数P(x|μ_vi,σ_vi)和平滑后强度矩阵E'的概率密度分布函数P(y|μ_λi,σ_λi)后，得到样本平滑后功率矩阵P'_v的概率密度分布p_v和平滑后强度矩阵E'的概率密度分布p_λ：

$p_v=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{vi}}}\exp (-\frac{{(x-{\mathrm{\μ}}_{\mathrm{vi}})}^2}{2{\mathrm{\σ}}_{\mathrm{vi}}^2}),$

$p_{\mathrm{\λ}}=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{\λi}}}\exp (-\frac{{(y-{\mathrm{\μ}}_{\mathrm{\λi}})}^2}{2{\mathrm{\σ}}_{\mathrm{\λi}}^2}).$

步骤8.构造联合概率密度分布函数，并计算联合概率密度分布。

8a)利用步骤(7)求得的k类样本的概率密度分布函数，构造平滑后功率矩阵P'_v和平滑后强度矩阵E'的联合概率密度分布函数P(a_j|a_j∈C_i)：

P(a_j|a_j∈C_i)=P(x|μ_vi,σ_vi)·P(y|μ_λi,σ_λi)，

其中a_j为第j个像素点，j=1,2,...,z,C_i表示第i类样本集合,i=1,2,...,k；

8b)将步骤(7d)求得的k类样本平滑后功率矩阵P'_v的概率密度分布p_v和平滑后强度矩阵E'的概率密度分布p_λ代入平滑后功率矩阵P'_v和平滑后强度矩阵E'的联合概率密度分布函数P(a_ja_j∈C_i)后，得到平滑后功率矩阵P'_v和平滑后强度矩阵E'的联合概率密度分布p(a_j|a_j∈C_i)：

p(a_j|a_j∈C_i)=p_v·p_λ。

步骤9.根据最大后验概率MAP准则，依据表达式：

a_j∈C_i,st.p(a_j|C_i)=maxp(a_j|a_j∈C_i)，i=1,2,...,k，j=1,2,...,z，对待分类极化SAR数据进行贝叶斯分类，即利用步骤（8）得到的每个样本的联合概率密度分布p(a_j|a_j∈C_i)，将像素点a_j划分到概率密度分布p(a_j|a_j∈C_i)最大的那一类，如图5所示，图5中颜色条上的数字1-8表示不同的分类标签。

步骤10.输出分类结果。

本发明的效果可以通过如下仿真实验进行验证。

1.实验条件

1.1）试验中选取的极化SAR图像如图2所示，该图为荷兰Flevoland农田区域的L波段极化SAR的RGB合成图像，其中R代表|HH-VV|，G代表|HV|+|VH|，B代表|HH+VV|，图像尺寸大小为750像素×1024像素。该极化SAR图像中区域的地表覆盖物为各种农作物和水域。

2.2）实验中选取8类区域作为训练样本，8类样本分别代表不同的地物。

2.实验内容

2.1）对8类训练样本功率矩阵和强度矩阵的实际分布和估计分布进行实验，其结果如图3和4所示，其中图3和4中的实线为数据的实际分布，虚线为数据的估计分布。

从图3和4可以看出，数据实际分布和估计分布之间略有偏差，这是估计参数设置以及噪声问题造成的，但总体上拟合结果较好，这表明极化SAR数据的散射功率和散射强度分布特性可以通过仿真很好的估计出来。

2.2）将本发明与两种经典算法H/alpha和复Wishart对比，分类精度如表一所示：

表一三种算法分类精度对比表

	本发明	H/alpha	复Wishart
				样本1	0.6089	0.9400	0.8842
样本2	0.8487	0	0.9987
				样本3	0.9247	0.9355	0.8395
样本4	0.6653	0	0
				样本5	0.9767	0.9044	0.5837
样本6	0.7624	0	0
				样本7	0.5370	0	0
样本8	0.8995	0	0.5450
				平均	0..7779	0.3475	0.4814

从表一可以看出，本发明的平均分类精度比H/alpha和复Wishart高。H/alpha和复Wishart算法都是基于散射特性的分类方法，这类方法对散射特性相似但实际地物并不同类的区域容易造成错分，从而使得该区域分类精度为0。然而，本发明是基于散射功率和强度联合统计的分类算法，对于散射特性相似的区域，其散射强度是不同的，从而可以对这类地物实现分类。

2.3）利用本发明方法对极化SAR图像进行分类，结果如图5所示，从图5可以看出，本发明得到的分类结果区域一致性较好，这就表明本发明可以有效提高极化SAR图像分类的精度。

Claims (4)

1.一种基于散射功率和强度联合统计的极化SAR图像分类方法，包括如下步骤：

(1)将待分类极化SAR图像的相干矩阵T作为输入，T的大小为3×3×z，其中z是图像像素点总个数；

(2)采用Lee滤波器对相干矩阵T进行去噪，得到去噪后的相干矩阵T'；

(3)对去噪后的相干矩阵T'进行Freeman分解，得到体散射功率，取每个像素点的体散射功率构成大小为z×1的功率矩阵P_v；

(4)对去噪后的相干矩阵T'进行特征值分解，提取极化SAR图像的散射强度信息，得到每个像素点的三个特征值λ₁,λ₂,λ₃，取最大特征值λ₁构成大小为z×1的强度矩阵E；

(5)对功率矩阵P_v和强度矩阵E分别进行大小为m×n的变形，得到变形后的功率矩阵P_vT和强度矩阵E_T，m,n分别是待分类图像的行数与列数，对变形后功率矩阵P_vT和变形后的强度矩阵E_T的每个像素点分别进行8邻域求平均来代替该像素点的值，得到平滑后的功率矩阵P'_v和平滑后的强度矩阵E'；

(6)依据实际地物分布，选择不同的k类匀质区域作为训练样本，k=8；

(7)采用EM算法分别对k类样本的平滑后功率矩阵P'_v和平滑后强度矩阵E'的高斯参数进行估计，得到k类样本的平滑后功率矩阵P'_v的概率密度分布函数P(x|μ_vi,σ_vi)和平滑后强度矩阵E'的概率密度分布函数P(y|μ_λi,σ_λi)：

$P\left(x\right|{\mathrm{\μ}}_{\mathrm{vi}},{\mathrm{\σ}}_{\mathrm{vi}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{vi}}}\exp (-\frac{{(x-{\mathrm{\μ}}_{\mathrm{vi}})}^2}{2{\mathrm{\σ}}_{\mathrm{vi}}^2}),$

$P\left(y\right|{\mathrm{\μ}}_{\mathrm{\λi}},{\mathrm{\σ}}_{\mathrm{\λi}})=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{\mathrm{\λi}}}\exp (-\frac{{(y-{\mathrm{\μ}}_{\mathrm{\λi}})}^2}{2{\mathrm{\σ}}_{\mathrm{\λi}}^2}),$

其中x为平滑后功率矩阵P'_v的值，y为平滑后强度矩阵E'的值、μ_vi为第i类样本的平滑后功率矩阵P'_v的均值，σ_vi为第i类样本的平滑后功率矩阵P'_v的方差，μ_λi为第i类样本的平滑后强度矩阵E'的均值，σ_λi为第i类样本的平滑后强度矩阵E'的方差，i=1,2,...,k；

(8)利用步骤(7)求得的k类样本的概率密度分布函数，构造平滑后功率矩阵P'_v和平滑后强度矩阵E'的联合概率密度分布函数P(a_j|a_j∈C_i)：

P(a_j|a_j∈C_i)=P(x|μ_vi,σ_vi)·P(y|μ_λi,σ_λi)，

其中a_j为第j个像素点，j=1,2,...,z,C_i表示第i类样本集合,i=1,2,...,k；

(9)根据最大后验概率MAP原则，依据表达式：

a_j∈C_i,st.p(a_j|C_i)=maxP(a_j|a_j∈C_i)，i=1,2,...,k，j=1,2,...,z，对待分类极化SAR数据进行贝叶斯分类，并输出分类结果。

2.根据权利要求1所述的方法，其中步骤（3）所述的对去噪后的相干矩阵T'进行Freeman分解，利用如下公式进行：

P_v=8F_v/3，其中F_v为体散射分量的分解系数。

3.根据权利要求书1所述的方法，其中步骤（4）所述的对去噪后的相干矩阵T'进行特征值分解，是利用MATLAB软件中的eigs函数进行。

4.根据权利要求书1所述的方法，其中步骤（7）所述的采用EM算法对k类样本的平滑后功率矩阵P'_v和强度矩阵E'的高斯模型参数进行估计，按照如下步骤进行：

4a)将各类样本的平滑后功率矩阵P'_v和强度矩阵E'分别作为EM算法的输入数据；

4b)设定EM算法迭代终止条件为迭代误差小于0.00001，分别将第i类样本的平滑后功率矩阵P'_v的均值μ_vi和方差σ_vi，第i类样本的平滑后强度矩阵E'的均值μ_λi和方差σ_λi初始化为随机值，其中i=1,2,...,k；

4c)通过EM算法中的期望最大化对步骤（4b）中所述的参数μ_vi,μ_λi,σ_vi和σ_λi进行更新，每次更新后计算迭代误差并判断是否满足迭代终止条件，若满足迭代终止条件，则迭代停止，输出参数μ_vi,μ_λi,σ_vi和σ_λi，否则继续执行步骤（4c）。

Images