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

摘要

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

权利要求

1.一种基于K-SVD和稀疏表示的极化SAR图像分类方法，包括如下步骤：(1)将表示极化SAR图像的每个像素点的大小为3*3的极化相干矩阵T作为输入数据，计算每个像素点的大小为3*3的协方差矩阵C；

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

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

(4)利用K-SVD算法解如下公式： Subject to 得到训练字典D＝[D1,D2,...,DK]，Di是第i类训练样本集Yi训练得到的字典，i＝1,2,…,K；

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

(5)利用OMP算法解如下公式：min||Xi||0 Subject to Fj＝DiXi，i＝1,2,…,K，更新K个稀疏系数矩阵X＝[X1,X2,...,XK]，Xi是第i类字典对应的稀疏系数矩阵，i＝1,2,…,K，式中，Fj为特征矩阵中F的第j个列向量，即为特征向量，j＝1,2,…,N，N为待分类的总像素点数；

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

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

2.根据权利要求1所述的基于K-SVD和稀疏表示的极化SAR图像分类方法，其中，步骤(1)所述的将表示极化SAR图像的每个像素点的大小为3*3的极化相干矩阵T作为输入数据，计算每个像素点的大小为3*3的协方差矩阵C，按如下步骤进行：

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

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

C＝M*T*M’，

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

0sqrt(2)0]，sqrt(2)表示求2的平方根。

3.根据权利要求1所述的基于K-SVD和稀疏表示的极化SAR图像分类方法，其中，步骤(6)所述的利用步骤(5)求出的稀疏系数矩阵重构特征向量Fj，通过求解重构特征向量与原向量的差值计算Fj的类别，按如下步骤进行：

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

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

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

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

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

说明书

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

技术领域

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

背景技术

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

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

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

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

发明内容

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

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

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

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

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

[0011] (4)利用K-SVD算法解如下公式： Subject to 得到训练字典D＝[D1,D2,...,DK]，Di是第i类训练样本集Yi训练得到的字典，i＝1,2,…,K，[0012] 式中，X为稀疏系数矩阵，min||·||表示让·的值达到最小，Subject to表示约束条件， 表示任意第i列。||·||0表示向量的0范数， 为矩阵的2范数的平方，T0为稀疏表示中稀疏向量中非零值的个数的上限；

[0013] (5)利用OMP算法解式：min||Xi||0 Subject to Fj＝DiXi，i＝1,2,…,K，更新K个稀疏系数矩阵X＝[X1,X2,...,XK]，Xi是第i类字典对应的稀疏系数矩阵，i＝1,2,…,K，[0014] 式中，Fj为特征矩阵中F的第j个列向量，即为特征向量，j＝1,2,…,N，N为待分类的总像素点数；

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

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

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

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

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

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

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

附图说明

[0022] 图1是本发明的流程图；

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

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

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

具体实施方式

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

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

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

[0029]

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

[0031] C＝M*T*M’，

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

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

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

[0035]

[0036]

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

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

[0039]

[0040] 式中，H表示水平极化，V表示垂直极化，SHH表示水平发射水平接收的回波数据，SVV表示垂直发射垂直接收的回波数据，SHV表示水平发射垂直接收的回波数据，*表示矩阵的共轭，<·>表示按视数平均；

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

[0042]

[0043] 式中，s为表面散射分量的协方差矩阵，d为二面角散射分量的协方差矩阵，v为体散射分量的协方差矩阵，fs为平面散射分量的分解系数，fd为二面角散射分量的分解系数，fv为体散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射反射系数的比值，α被定义为α＝RghRvhRgvRvv，Rgh和Rgv分别表示地表的水平及垂直反射系数，Rvh和Rvv表示垂直墙体的水平及垂直反射系数；

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

[0045]

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

[0047] 2f)根据求解出的fs、fv、fd、α、β，按照下式求解出散射功率Pv、Pd、Ps：

[0048]

[0049] 式中，Ps表示表面散射功率，Pd表示二面角散射功率，Pv表示体散射功率；

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

[0051] SPAN＝T11+T22+T33，

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

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

[0054] 步骤四，利用K-SVD算法对K类训练样本集Y进行训练，得到K个训练字典D1,D2,...,DK，Di是第i类训练样本集Yi训练得到的字典，i＝1,2,…,K。

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

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

[0057]

[0058] 式中，dj为D的第j列原子， 为X的第j行，L为D的总列数，dz为D的第z列原子， 为X的第z行，Ez为去除第z列原子dz进行稀疏分解所产生的误差矩阵；

[0059] 4c)对变形后所得的公式 乘以矩阵Ωz，得到目标分解公式

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

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

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

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

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

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

[0066]

[0067] 6b)根据下式计算Fj与 的差值：

[0068]

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

[0070] 6c)根据下式计算Fj所属的类别：

[0071]

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

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

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

[0075] 1.实验条件和与方法：

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

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

[0078] 2.实验内容与结果分析

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

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

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

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

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

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

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

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

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

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

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

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

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