基于K-SVD和稀疏表示的极化SAR图像分类方法

摘要

本发明公开了一种基于K-SVD和稀疏表示的极化SAR图像分类方法，解决现有方法对分类类别数目限制以及极化特征信息利用不充分的问题。其实现步骤是：(1)将极化SAR的极化相干矩阵作为输入数据，计算协方差矩阵；(2)对每个像素点提取相干矩阵、协方差矩阵、Ps、Pd、Pv、H、α等，组成特征矩阵；(3)根据实际地物分布，选取训练样本，组成初始字典；(4)用K-SVD算法训练初始字典，得到训练字典；(5)将特征矩阵用训练字典表示，用OMP算法求解稀疏系数；(6)用求解出的稀疏系数重构特征矩阵，确定像素点的类别，得到最终分类结果。本发明利用极化SAR图像的极化特征，不限制分类类别数目，可用于极化SAR图像分类。

说明书

技术领域

本发明属于图像处理技术领域，涉及极化SAR图像的分类，该方法可用于极化 SAR图像目标分类与识别。

背景技术

雷达是一种可以实现全天候工作的主动探测系统，它可以穿透一定的地表，并 且可以改变发射波的频率、强度。合成孔径雷达技术(SAR)是成像雷达技术的一 种，它是利用雷达与目标的相对运动把尺寸较小的真实无线孔径用数据处理的方法 合成一个较大的等效天线孔径雷达。而极化SAR是用来测量回波信号的新型雷达， 它可以记录不同极化状态组合回波的相位差信息，大大提高了对地物的识别能力。 极化SAR图像分类是极化SAR图像解译的重要步骤，是边缘提取、目标检测和识别 的基础，可广泛应用于军事侦察、地形测绘、农作物生长监测等领域。

目前经典的极化SAR图像分类方法有：

1997年，Cloude等人提出了一种基于H/α极化分解的分类方法，该方法通过 Cloude分解得到特征参数散射熵H和散射角α，然后根据两个特征参数不同的值，将 目标分成8类。该方法的缺陷是位于类别边界特征相似的像素点会以随机的方式分配 给不同的类别并且这两个特征不足以表示所有的极化SAR信息。

1999年，Lee等人提出了一种基于H/α极化分解和复Wishart分类器的 H/α-Wishart分类方法，该方法将H/α极化分解方法得到的结果作为复Wishart分类 器的初始分类，对划分好的8个类别中的每一个像素进行重新划分，从而提高分类的 精度。该方法的缺陷是将分类类别数固定为8类，不能适应不同类别数的地物分类。 2004年，J.S.Lee等人提出了一种基于Freeman-Durden分解的极化SAR图像分类方法， 该方法首先通过Freeman分解获取表征散射体散射特征的三个特征：平面散射功率、 二面角散射功率和体散射功率，然后根据这三个特征的大小对数据进行初始划分， 然后利用Wishart分类器进行进一步精确划分。但是该方法由于Freeman分解中的多类 的划分以及合并，计算复杂度较高。

发明内容

本发明的目的在于针对上述现有技术方法的不足，提出了一种基于K-SVD和稀 疏表示的极化SAR图像分类方法。

本发明的具体实现步骤如下：

(1)将表示极化SAR图像的每个像素点的大小为3*3的极化相干矩阵T作为输 入数据，计算每个像素点的大小为3*3的协方差矩阵C；

(2)对极化SAR图像的每个像素点的极化相干矩阵T，采用克劳德cloude分解方 法得到H、α共两个散射参数，采用弗里曼-德登Freeman-Durden分解方法得到 P_s、P_d、P_v共三个散射参数，将相干矩阵T的9个元素、协方差矩阵C的9个元素、H、 α、P_s、P_d、P_v、总功率特征参数SPAN共24个特征作为每个像素的特征向量，该特 征向量的维数为24×1。用大小为24×N的特征矩阵F表示整幅极化SAR图像，N为 整幅极化SAR图像总的像素点数；

(3)根据实际地物分布，从每类点所对应的特征向量中选取一部分作为对应的 训练样本集Y＝Y₁,Y₂,...,Y_K，K为总类别数，Y_i是从第i类特征向量中选取的训练样本， i＝1,2,…,K；

(4)利用K-SVD算法解如下公式：Subject to得到训练字典D＝[D₁,D₂,...,D_K]，D_i是第i类训练样本集Y_i训练得到的字典， i＝1,2,…,K，

式中，X为稀疏系数矩阵，min||·||表示让·的值达到最小，Subject to表示约束 条件，表示任意第i列。||·||₀表示向量的0范数，为矩阵的2范数的平方，T₀为 稀疏表示中稀疏向量中非零值的个数的上限；

(5)利用OMP算法解式：min||X_i||₀ Subject to F_j＝D_iX_i，i＝1,2,…,K，更新K 个稀疏系数矩阵X＝[X₁,X₂,...,X_K]，X_i是第i类字典对应的稀疏系数矩阵， i＝1,2,…,K，

式中，F_j为特征矩阵中F的第j个列向量，即为特征向量，j＝1,2,…,N，N为待分 类的总像素点数；

(6)利用步骤(5)求出的稀疏系数矩阵重构特征向量F_j，通过求解重构特征向 量与原向量的差值确定F_j的类别；

(7)重复步骤(6)直到图像中所有的点都分类完毕，根据像素点所属的类别， 用不同颜色将不同类别的像素点表示出来。

本发明首先输入一幅极化SAR图像，提取它的极化特征，根据实际地物分布选 取训练样本，将训练样本的特征作为初始字典。其次，用K-SVD对初始字典进行训 练得到训练字典。然后，将图像用训练字典进行稀疏表示。最后，通过求解稀疏系 数进行极化SAR图像分类。本发明在不限制分类类别数目的同时，有效地提高了极 化SAR图像的分类精度，并且有效地利用了极化SAR图像的极化特征信息。与现 有的技术相比，本发明具有以下优点：

1、本发明将极化SAR图像的极化特征作为初始字典，有效地利用了极化SAR 图像的极化特征信息；

2、本发明利用K-SVD算法得到较好的训练字典，使得分类结果有更好的鲁棒 性；

3、本发明根据字典对样本的重构误差最小化原则来确定样本的类别，避免使用 一般距离测度算法带来的缺陷，得到更好的极化SAR图像分类结果；

4、仿真结果表明，本发明方法较经典的H/α极化分解的分类方法和 H/α-Wishart分类方法能更有效的对极化SAR图像进行分类。

附图说明

图1是本发明的流程图；

图2是本发明仿真所用的两幅测试图像；

图3是本发明与现有两种方法对San Francisco数据的分类实验结果对比图；

图4为本发明与现有两种方法对Flevoland数据划分的分类实验结果对比图；

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤一，将极化SAR图像的每个像素点的极化相干矩阵T作为输入数据，计算 每个像素点的协方差矩阵C。

1a)读入极化SAR图像的每个像素点的极化相干矩阵T，每个像素点的极化相 干矩阵T为一个含有9个元素的3×3的矩阵：

$T=\left(\begin{array}{ccc}{T}_{11}& T_{12}& T_{13}\\ T_{21}& T_{22}& T_{23}\\ T_{31}& T_{32}& T_{33}\end{array}\right);$

1b)通过下式计算每个像素点的协方差矩阵C：

C＝M*T*M’，

式中，M’表示计算矩阵M的转置矩阵，计算矩阵M＝[1/sqrt(2)]*m，m＝[101；10 -1；0sqrt(2)0]，sqrt(2)表示求2的平方根；

步骤二，对于极化SAR图像的每个像素点，采用克劳德cloude分解方法得到H、 α共两个散射参数，采用弗里曼-德登Freeman-Durden分解方法得到P_s、P_d、P_v共三 个散射参数；将相干矩阵T的9个元素、协方差矩阵C的9个元素、H、α、P_s、P_d、P_v、 总功率特征参数SPAN共24个特征作为每个像素的特征向量，维数为24×1；整幅 极化SAR图像组成大小为24×N的特征矩阵F，N为总像素点数。

2a)当T矩阵是一个3*3的矩阵时，其有3个特征值，则按照下式，得到极化 SAR图像Cloude分解的两个参数H和α：

$H=\underset{i=1}{\overset{3}{\Sigma }}-P_i{\log }_3{P}_i$

$\alpha =\underset{i=1}{\overset{3}{\Sigma }}{P}_i{\alpha }_i$

式中，P_i表示步骤一中的极化SAR图像相干矩阵T的第i个特征值与所有特征 值总和的比值，H表示极化SAR图像散射矩阵的散射熵，α表示极化SAR图像散射 矩阵的散射角；

2b)由步骤1b)计算得到的极化SAR图像的协方差矩阵C可表示为：

$C=\left(\begin{array}{ccc}⟨{\left|S_{\mathrm{HH}}\right|}^2⟩& \sqrt{2}⟨S_{\mathrm{HH}}{S}_{\mathrm{HV}}^{*}⟩& ⟨S_{\mathrm{HH}}{S}_{\mathrm{VV}}^{*}⟩\\ \sqrt{2}⟨S_{\mathrm{HV}}{S}_{\mathrm{HH}}^{*}⟩& 2⟨{\left|S_{\mathrm{HV}}\right|}^2⟩& \sqrt{2}⟨S_{\mathrm{HV}}{S}_{\mathrm{VV}}^{*}⟩\\ ⟨S_{\mathrm{VV}}{S}_{\mathrm{HH}}^{*}⟩& \sqrt{2}⟨S_{\mathrm{VV}}{S}_{\mathrm{HV}}^{*}⟩& ⟨{\left|S_{\mathrm{VV}}\right|}^2⟩\end{array}\right)---\left(1\right)$

式中，H表示水平极化，V表示垂直极化，S_HH表示水平发射水平接收的回波 数据，S_VV表示垂直发射垂直接收的回波数据，S_HV表示水平发射垂直接收的回波数 据，*表示矩阵的共轭，<·>表示按视数平均；

2c)将协方差矩阵C按下式分解：

$\left(\begin{array}{c}C={⟨C⟩}_s+{⟨C⟩}_d+{⟨C⟩}_v\\ =f_s\left(\begin{array}{ccc}{\left|\beta \right|}^2& 0& \beta \\ 0& 0& 0\\ {\beta }^{*}& 0& 1\end{array}\right)+f_d\left(\begin{array}{ccc}{\left|\alpha \right|}^2& 0& \alpha \\ 0& 0& 0\\ {\alpha }^{*}& 0& 1\end{array}\right)+f_v\left(\begin{array}{ccc}1& 0& 1/3\\ 0& 2/3& 0\\ 1/3& 0& 1\end{array}\right)\end{array}\right)---\left(2\right)$

式中，<C>_s为表面散射分量的协方差矩阵，<C>_d为二面角散射分量的协方差矩 阵，<C>_v为体散射分量的协方差矩阵，f_s为平面散射分量的分解系数，f_d为二面角 散射分量的分解系数，f_v为体散射分量的分解系数，β是水平发射水平接收后向散 射反射系数与垂直发射垂直接收后向散射反射系数的比值，α被定义为 α＝R_ghR_vhR_gvR_vv，R_gh和R_gv分别表示地表的水平及垂直反射系数，R_vh和R_vv表示垂 直墙体的水平及垂直反射系数；

2d)将式(1)中协方差矩阵C的元素与式(2)中矩阵的元素相对应，获得一 个具有五个未知数f_s、f_v、f_d、α、β和四个方程的方程组如下：

$\left(\begin{array}{c}⟨{\left|S_{\mathrm{HH}}\right|}^2⟩=f_s{\left|\beta \right|}^2+f_d{\left|\alpha \right|}^2+f_v\\ ⟨{\left|S_{\mathrm{VV}}\right|}^2⟩=f_s+f_d+f_v\\ ⟨{\left|S_{\mathrm{HH}}{S}_{\mathrm{VV}}^{*}\right|}^2⟩=f_s\beta +f_d\alpha +f_v/3\\ ⟨{\left|S_{\mathrm{HV}}\right|}^{2⟩}=f_v/3\end{array}\right)---\left(3\right)$

2e)计算像素点协方差矩阵C中的的值并判断正负，如果 则α＝-1，如果则β＝1，给定α或β的值后，剩 余的4个未知数可根据式(3)求解得出，其中Re(·)表示取实部；

2f)根据求解出的f_s、f_v、f_d、α、β，按照下式求解出散射功率P_v、P_d、P_s：

$\left(\begin{array}{c}{P}_v=\frac{8f_v}{3}\\ P_d=f_d(1+{\left|\alpha \right|}^2)\\ P_s=f_s(1+{\left|\beta \right|}^2)\end{array}\right)$

式中，P_s表示表面散射功率，P_d表示二面角散射功率，P_v表示体散射功率；

2g)按照下式计算总功率SPAN：

SPAN＝T11+T22+T33，

式中，T11、T22、T33为极化SAR图像相干矩阵对角线上的三个元素。

步骤三，根据实际地物分布，从每类点所对应的特征向量中选取一部分作为对 应的训练样本集Y＝Y₁,Y₂,...,Y_K，K为总类别数。Y_i是从第i类特征向量中选取的训 练样本，i＝1,2,…,K。

步骤四，利用K-SVD算法对K类训练样本集Y进行训练，得到K个训练字典 D₁,D₂,...,D_K，D_i是第i类训练样本集Y_i训练得到的字典，i＝1,2,…,K。

4a)从K-SVD算法中给出总的优化分布Subject to 其中D＝[D₁,D₂,...,D_K]，D_i是第i类训练样本集Y_i训练得到的字典， i＝1,2,…,K，X为稀疏系数矩阵，min||·||表示让·的值达到最小，Subject to表示约 束条件，表示任意第i列。||·||₀表示向量的0范数，为矩阵的2范数的平方， T₀为稀疏表示中稀疏向量中非零值的个数的上限；

4b)对总的优化公式中的进行变形得到：

${\left|\right|Y-\mathrm{DX}\left|\right|}_2^2={\left|\right|Y-\underset{j=1}{\overset{L}{\Sigma }}{d}_j{x}_j^T\left|\right|}_2^2={\left|\right|(Y-\underset{j\ne z}{\overset{L}{\Sigma }}{d}_j{x}_j^T)-d_z{x}_z^T\left|\right|}_2^2={\left|\right|E_z-d_x{x}_z^T\left|\right|}_2^2,$

式中，d_j为D的第j列原子，为X的第j行，L为D的总列数，d_z为D的 第z列原子，为X的第z行，E_z为去除第z列原子d_z进行稀疏分解所产生的误 差矩阵；

4c)对变形后所得的公式乘以矩阵Ω_z，得到目标分解公式 ${\left|\right|E_z{\Omega }_z-d_z{x}_z^T{\Omega }_z\left|\right|}_2^2={\left|\right|E_z^R{-d}_z{x}_z^R\left|\right|}_2^2,$

式中，变形误差矩阵为误差矩阵E_z的变形，是E_z去掉零输入后的收 缩结果，Ω_z的大小为P×|ω_z|，P为训练样本集Y的列数， |ω_z|为ω_z的模值，且Ω_z在(ω_z(j),j)处为1，其他地方全 为0，其中1≤j≤|ω_z|，ω_z(j)为ω_z的第j个数；

4d)用SVD方法分解变形误差矩阵式中，U表示左奇异矩 阵，V^TV表示右奇异矩阵，Δ表示奇异值矩阵；用U的第一列更新目标训练字典D 的第z列原子d_z；

4e)重复步骤4c)到步骤4d)，直到对D中的所有原子进行更新处理，得到K 个新的字典[D₁,D₂,...,D_K]；

步骤五，利用OMP算法解如下公式：min||X_i||₀ Subject to F_j＝D_iX_i， i＝1,2,…,K，更新K个稀疏系数矩阵X＝[X₁,X₂,...,X_K]，X_i是第i类字典对应的稀疏 系数矩阵，i＝1,2,…,K，式中，F_j为特征矩阵中F的第j个列向量，即为特征向量， j＝1,2,…,N，N为待分类的总像素点数。

步骤六，利用步骤五求出的稀疏系数矩阵重构特征向量F_j，通过求解重构特征 向量与原向量的差值计算F_j的类别。

6a)对于每一个像素点的特征向量F_j，根据下式重构特征向量F_j：

${\hat{F}}_j=D_i{X}_i,i=\mathrm{1,2},...,K;$

6b)根据下式计算F_j与的差值：

$r_i\left(F_j\right)={\left|\right|F_j-{\hat{F}}_j\left|\right|}_2,$

式中，||·||₂表示计算矩阵的2范数；

6c)根据下式计算F_j所属的类别：

$\mathrm{Class}\left(F_j\right)=\underset{i}{\min }{r}_i\left(F_j\right),i=1,...,K,$

式中，Class(F_j)表示F_j的所属类别，·表示·取最小值时i的值；

步骤七，重复步骤六直到极化SAR图像中所有的像素点都分类完毕，根据像素 点所属的类别，用不同颜色将不同类别的像素点表示出来。

本发明的效果可通过以下仿真进一步说明。

1.实验条件和与方法：

实验仿真环境：MATLAB 2009a，Windows XP Professional。

实验方法：分别为H/α极化分解的分类方法和H/α-Wishart分类方法和本发 明，其中前两种方法为极化SAR图像分类的经典方法。

2.实验内容与结果分析

实验内容：本发明使用的是图2所示的两组极化SAR图像数据，图2(a)为美 国San Francisco地区数据，视数为四，图2(b)为荷兰Flevoland地区的数据，视 数为四，两组数据都来源于美国宇航局喷气推进实验室的AIRSAR传感器。

仿真1，用本发明以及H/α极化分解的分类方法和H/α-Wishart分类方法对美 国San Francisco地区数据进行分类实验，结果如图3所示，其中：

图3(a)为H/α极化分解的分类方法的分类结果，分为9类；

图3(b)为H/α-Wishart分类方法的分类结果，分为9类；

图3(c)为用本发明方法的分类结果，分为3类，每类选取1000个像素点作为 训练样本。

从图3(a)来看，H/α极化分解的分类方法分类结果很不理想，很多区域并没 有被区分出来；

从图3(b)来看，H/α-Wishart分类方法的分类结果优于H/α极化分解的分 类方法，区域划分的更加细致，但图像细节保持较差；

从图3(c)来看，本发明的分类结果从视觉上看分类效果更好，其中跑马场、 高尔夫球场等区域在分类后的图中，区域一致性好于H/α极化分解的分类方法和 H/α-Wishart分类方法，金门大桥和小岛也都轮廓清楚。

仿真2，用本发明以及H/α极化分解的分类方法和H/α-Wishart分类方法对荷 兰Flevoland地区数据进行分类实验，结果如图4所示，其中：

图4(a)为H/α极化分解的分类方法的分类结果，分为9类；

图4(b)为H/α-Wishart分类方法的分类结果，分为9类；

图4(c)为用本发明方法的分类结果，分为13类，每类选取300个像素点作为 训练样本。

从图4可以看出，H/α-Wishart分类方法和H/α极化分解的分类方法由于固定 了类别数目，对此图不能很精确的划分，很多类别都被归为了一类，而本发明清晰 了地分出了各类别的轮廓，分类效果明显好于H/α-Wishart分类方法和H/α极化分 解的分类方法，区域划分更细致更准确。