CN103020654B - 基于样本产生与核局部特征融合的sar图像仿生识别方法 - Google Patents

Abstract

本发明提供了一种基于样本产生与核局部特征融合的SAR图像仿生识别方法，属于图像处理技术及SAR目标识别领域。该方法先构建了超完备训练样本集，进行训练得到几何流形，然后再对待识别样本进行识别，具体对各样本先采用K-SVD字典学习方法进行图像去噪，再用目标质心方法实现目标区域提取；然后结合局部相位量化（LPQ）与Gabor滤波方法分别进行特征提取，并进行特征融合，最后通过高维几何流形的覆盖进行分类，运用仿生模式进行识别。本发明方法的图像相干噪声抑制效果较为显著，能有效地提取SAR图像特征，解决了SAR图像因姿态角变化引起提取的特征不稳定的问题，识别精度高，具有很好的鲁棒性。

Description

基于样本产生与核局部特征融合的SAR图像仿生识别方法

技术领域

本发明属于图像处理技术及SAR目标识别领域，特别涉及一种基于样本产生与核局部特征融合的SAR图像仿生识别方法，可应用于军事和民用的目标识别中。

背景技术

合成孔径雷达(SyntheticApertureRadar，简称SAR)由于其能够全天候、全天时获取地面战场数据，在现代战场感知、对地军事打击等领域有良好的应用前景，使得基于SAR图像的自动目标识别(AutomaticTargetRecognition，简称ATR)受到了越来越多研究学者的重视，并且已经逐渐成为国内外的研究热点之一。在SAR目标识别的研究中，目前国内外许多研究机构大多都是基于由美国DARPA/AFRLMSTAR项目组提供的运动、静止目标探测和识别（MovingandStationaryTargetAcquisitionandRecognition,简称MSTAR）数据库进行的，MSTAR数据库对SAR图像目标识别的研究工作发挥了至关重要的作用。该数据库中所成的SAR图像包含幅度、相位和极化等丰富的特征信息，同时在SAR图像中还存在大量的乘性噪声。另一影响识别的重要因素是目标姿态角的变化导致难以提取有效稳定的特征。

因此，如何消除或减少噪声对目标识别的影响，以及如何有效地克服目标姿态角变化，以便提取目标旋转不变的稳定特征，成为了提高目标识别精度的关键技术之一。目前，对于消除噪声的干扰方法，通常采用同态滤波的方法，在对数域内进行噪声消除，然后再用指数变换生成原图像，但这一方法对于抑制SAR图像的相干噪声难以取得理想的方法。而对于目标姿态角变化的难点问题，通常采用对目标进行姿态角按区间分类的方法，构造多个分类器进行分类。这一解决方法需要已知目标姿态角或者估计目标姿态角，而且需要构造多个分类器，因此，复杂度较高。另一方面，寻找符合SAR目标特征特性的有效分类器，这也是SAR目标图像识别的难点问题之一。

发明内容

基于上述已有技术的缺陷与不足，本发明提出了一种样本产生与核局部特征融合的SAR图像识别方法。本发明方法可有效出去SAR图像中相干噪声干扰，并克服SAR图像姿态角变化引起特征提取不稳定的问题，从而提高SAR图像目标识别的正确识别率的目的。

本发明提供的一种基于样本产生与核局部特征融合的SAR图像仿生识别方法，其具体步骤如下：

步骤1：构建超完备训练样本集；对输入每一类SAR的原始训练图像分别采用样本旋转的方法产生虚拟训练样本图，由真实的SAR的原始训练图像和所产生的所有虚拟训练样本图构成最终的超完备训练样本集；

步骤2：采用K-SVD字典学习算法对输入的每个样本图像，进行噪声处理；

步骤3：对噪声处理后的各样本图像，采用质心方法提取目标识别区域；

步骤4：采用高斯核函数将所有经步骤3输出的样本图像映射到高维特征空间；

步骤5：采用二维Gabor滤波器对步骤4得到的样本图像进行特征提取，得到样本图像的特征图；

步骤6：采用局部相位量化方法对步骤4得到的样本图像进行特征提取，得到样本图像的特征图；

步骤7：将步骤5和步骤6得到的同一样本图像的特征图进行加权融合，得到各样本图像的融合特征图；

步骤8：根据超完备训练样本集中的各样本图像的融合特征图，构建超完备训练样本集中所有样本图像的仿生流形；所构建的仿生流形优选采用”超香肠”神经元模型；

步骤9：对待识别的样本图像输入步骤2中，经过步骤2至步骤7的处理，得到待识别的样本图像的融合特征图，然后采用仿生模式识别方法进行SAR目标图像识别。

所述的步骤2具体是：设Y为输入的样本集中一幅要进行噪声处理的原始样本图像，将样本图像Y分割成N块，则有y_i表示样本图像Y分割后的第i块，则通过K-SVD字典学习算法，找出一个最优字典D和最优稀疏分解系数矩阵X表示原始样本图像Y，设随机初始化的字典D⁽⁰⁾为N×K的矩阵，D⁽⁰⁾为当前的目标训练字典，令迭代次数J初始为1，然后进行如下步骤：

步骤2.1：优化每个图像小块y_i对应的稀疏分解局部系数得到更新的稀疏分解系数矩阵X。保持当前的目标训练字典D^(J-1)不变，采用追踪算法循环计算最优稀疏分解局部系数优化公式如下：

$\underset{x_i}{\min }\left\{{\left|\right|y_i-{\mathrm{Dx}}_i^r\left|\right|}_2^2\right\}\mathrm{subjectto}{\∀}_i,{\left|\right|x_i^r\left|\right|}_0\≤T_0,$ (i＝1,2,…,N)

subjectto表示约束条件，为稀疏分解系数矩阵X的第i列；为任意第i列，||·||₀为向量的0范数，T₀为稀疏度控制阈值，F表示F范数，D为当前的目标训练字典D^(J-1)。

步骤2.2：更新目标训练字典。保持步骤2.1优化得到的图像小块y_i的稀疏分解局部系数保持不变，对图像小块y_i进行循环学习，更新当前的目标训练字典D^(J-1)的第k列原子d_k，k＝1,2,3,…,K，具体过程如下：

步骤2.2.1：获取当前目标训练字典D^(J-1)的原子d_k；字典D⁽⁰⁾的第k列原子d_k的第i行元素d_k(i)为： 表示稀疏分解系数矩阵X的第k行第i列元素；

步骤2.2.2：确定误差矩阵E_k：其中，d_j为字典D的第j列原子，为稀疏分解系数矩阵X的第j行；

步骤2.2.3：将误差矩阵E_k乘以矩阵R_k得到新的误差矩阵 其中R_k为N×|d_k|的矩阵，|d_k|为d_k的模；

步骤2.2.4：对误差矩阵进行奇异值分解：Δ为所要求取的奇异值矩阵，U和V^T为两个正交矩阵，用正交矩阵U的第一列u₁更新第k列原子d_k：d_k＝u₁；

步骤2.2.5：重复步骤2.2.2至步骤2.2.4，对目标训练字典D^(J-1)中的每列原子都进行更新，更新完毕后得到新的目标训练字典D^(J)，目标训练字典D^(J)；

步骤2.3：更新迭代次数J＝J+1；然后判断迭代次数J是否达到设定的最大迭代次数，若是，停止学习，输出目标优化字典D^(J-1)以及相对应的稀疏分解系数矩阵X；若否，判断当前得到的目标优化字典D^(J-1)是否满足总优化公式的约束条件，若是，输出目标优化字典D^(J-1)以及相对应的稀疏分解系数矩阵X，否则转步骤2.1继续执行；

总优化公式为： (i＝1,2,…,N)；F表示F范数，D为当前得到的目标优化字典D^(J-1)；

步骤2.4：对样本图像Y进行去除噪声处理，得到新图像Y"＝D^(J-1)X。

所述的步骤6的具体实现方法是：

步骤6.1：进行短时傅里叶变换，对于样本图像中的每一个像素点(x,y)，获取该像素点在大小为M×M的邻域N_x上的局部相位信息F(u',v')：

$F(u^{\′},v^{\′})={\mathrm{\Σ}}_{y\∈N_x}f(x-y)e^{-j2{\mathrm{\πu}}^{\′T}y}=w_{u^{\′}}^T{f}_x$

(u',v')表示对坐标(x，y)进行短时傅里叶变换后的坐标，u'表示频率，u'^T表示u'的转置，f(x-y)表示样本图像函数f(x)的卷积，w_u'是二维离散傅里变换在频率u'的基本向量，是w_u′的转置；f_x是包含N_x里面所有M²个像素点的向量，M为整数；

设a是一个标量频率，对于四个频率u'₁＝[a,0]^T,u'₂＝[0,a]^T,u'₃＝[a,a]^T,u'₄＝[a,-a]^T对应的四个局部相位信息分别为：F(u'₁,v'),F(u'₂,v'),F(u'₃,v′),F(u'₄,v')，则四个局部相位信息组成矩阵为： $F_x^c=[F({u^{\′}}_1,v^{\′}),F({u^{\′}}_2,v^{\′}),F({u^{\′}}_3,v^{\′}),F({u^{\′}}_4,v^{\′})],$ 进一步根据集合可得到转换系数矩阵F_x：其中，Re{*},Im{*}分别表示一个复数的实部和虚部；转换系数矩阵F_x对应的8*M²转换矩阵W是： $W=\left[\mathrm{Re}\right\{W_{{u^{\′}}_1},W_{{u^{\′}}_2},W_{{u^{\′}}_3},W_{{u^{\′}}_4}\},\mathrm{Im}\{W_{{u^{\′}}_1},W_{{u^{\′}}_2},W_{{u^{\′}}_3},W_{{u^{\′}}_4}\left\}\right],$ 其则有F_x＝Wf_x；

步骤6.2：获得转换系数矩阵F_x的协方差矩阵，设样本图像函数f(x)是一阶马尔可夫过程，临近像素值之间的相关系数是ρ，样本方差是σ²，设σ²＝1，则两个不同像素位置x_I和x_J的协方差σ_I,J表示为： ${\mathrm{\σ}}_{I,J}={\mathrm{\ρ}}^{\left|\right|x_I-x_J\left|\right|},$ ||*||表示L₂模；

M个样本在领域N_x的协方差矩阵C为：

因此，获得转换系数矩阵F_x的协方差矩阵D_x为：D_x＝WCW^T；

步骤6.3：对协方差矩阵D_x进行奇异值分解，得到后奇异值正交矩阵V^T，然后使用白变换将转换系数矩阵F_x达到高斯独立分布，得到矩阵G_x=V^TF_x；获取所有像素点对应的矩阵G_x，然后对得到的所有矩阵G_x进行量化，量化公式为：其中，q_j表示量化值，g_i表示矩阵G_x中的第j个元素；

用0-255之间的二进制编码整数值b来代表量化的系数：

最后，生成一个所有像素点的二进制编码整数值的直方图。

本发明方法的优点与积极效果在于：

（1）本发明方法的图像相干噪声抑制效果较为显著，从而保证了有效地提取SAR图像特征。本发明通过K-SVD字典学习的方法，进行图像噪声抑制，其抑制效果较为显著，获得了高信噪比的SAR图像，解决了相干噪声对于SAR图像识别的干扰问题，从而保证了有效地提取SAR图像特征；（2）本发明提出采用Gabor和LPQ方法进行姿态角不变特征的提取及加权融合，能够有效解决SAR图像因姿态角变化引起提取的特征不稳定的问题；（3）识别精度高。首次将仿生模式识别方法应用于SAR目标识别中，通过实验证明，获得了较高的识别率，并具有很好的鲁棒性；（4）本发明提出了旋转样本产生的方法。采用真实SAR样本进行旋转产生了虚拟SAR样本的方法，提高了同类样本的流形连续性，提高了仿生模式几何流形覆盖体的构造精度，从而提高了识别率。

附图说明

图1是本发明的SAR图像识别方法的实现流程图；

图2是本发明采用的三类军事目标光学图像与SAR图像示意图；

图3是本发明提出的旋转样本产生方法实现流程图；

图4是本发明采用不同尺度的局部相位量化方法进行特征提取的效果示意图；

图5是本发明采用Gabor滤波方法进行特征提取的效果示意图；

图6是本发明仿生模式识别构建及判决示意图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明提供的SAR图像识别方法，先采用K-SVD字典学习方法进行图像去噪、再用目标质心方法实现目标区域提取；然后结合局部相位量化（LocalPatternQuantization,LPQ）与Gabor滤波方法进行特征提取并进行特征融合，最后运用仿生模式识别方法进行高维几何流形的覆盖进行分类。如图1所示，为本发明提供的SAR图像识别方法，下面对各步骤进行详细说明。

步骤1：构建超完备训练样本集。

本发明实施例中采用的三类军事目标光学图像与SAR图像如图2所示，图2所示的三类图像为T72，BTR70和BMP2。对输入每一类的SAR的原始训练图像分别采用样本旋转的方法产生虚拟样本，由真实的SAR的原始训练样本和虚拟的训练样本构成最终的超完备训练样本集，具体如图3所示，通过样本旋转方法产生虚拟训练样本：

步骤1.1：确定旋转角度θ₁。θ₁一般选取能被360整除的数值。

步骤1.2：计算每一类真实训练样本的虚拟样本数为：(360°/θ₁)-1。

步骤1.3：每旋转一次角度θ₁，得到对应旋转角度下的虚拟训练样本图；如图3所示，旋转角度θ₁为60度，得到5个虚拟训练样本图：SAR虚拟样本1～5。

步骤2：采用K-SVD字典学习算法对输入的超完备训练样本集或识别样本集中的每个样本图像进行噪声处理。本发明方法首先输入超完备训练样本集进行处理，构建超香肠神经元链，在构建好之后再输入识别样本进行SAR图像仿生识别。

K-SVD字典学习算法能够有效地去除噪声，获得更高信噪比的图像，其目的在于找出一个最优的字典D和最优的稀疏分解系数X来表示原始样本数据。采用K-SVD字典学习算法进行噪声处理的具体过程如下：

在K-SVD字典学习算法中，总优化公式为：

$\underset{D,X}{\min }\left\{{\left|\right|Y-\mathrm{DX}\left|\right|}_F^2\right\}\mathrm{subjectto}{\∀}_i,{\left|\right|x_i^r\left|\right|}_0\≤T_0,$ (i＝1,2,…,N)

Y为输入的一幅原始样本图像，对输入的样本图像Y分割成N小块，则有y_i表示原始图像分割后的第i块，Y是一个N行×N列的矩阵。D为目标训练字典，X表示稀疏分解系数矩阵，为稀疏分解系数矩阵X的第i列，subjectto表示约束条件，为任意i，||·||₀为向量的0范数，T₀为稀疏度控制阈值，为一整数，由用户根据需要设置。F表示F范数。

设初始化的随机规范字典为D⁽⁰⁾，D⁽⁰⁾为N行×K列的矩阵，为当前目标训练字典。初始迭代次数J＝1。

步骤2.1：优化每个图像小块y_i对应的稀疏分解局部系数得到更新的稀疏分解系数矩阵X。保持目标训练字典D^(J-1)不变，采用任意追踪算法循环计算原始样本图像Y分割小块y_i对应的最优稀疏分解局部系数其优化公式如下：

$\underset{x_i}{\min }\left\{{\left|\right|y_i-{\mathrm{Dx}}_i^r\left|\right|}_2^2\right\}\mathrm{subjectto}{\∀}_i,{\left|\right|x_i^r\left|\right|}_0\≤T_0,$ (i＝1,2,…,N)

公式中的D为当前的目标训练字典D^(J-1)。

步骤2.2：更新目标训练字典。保持由步骤2.1求解的图像小块y_i的稀疏分解局部系数不变，对图像小块y_i进行循环学习，更新当前的目标训练字典D^(J-1)的第k列原子d_k，k＝1,2,3,…,K，从而得到优化的原子集合，由优化后的原子集合构成与步骤2.1更新的稀疏分解系数矩阵X对应的当前新的目标训练字典D^(J)。由D⁽⁰⁾优化得到的新的当前目标训练字典为D⁽¹⁾。具体由目标训练字典的更新过程如下：

步骤2.2.1：获取当前目标训练字典D^(J-1)的原子d_k。

定义初始化随机规范字典D⁽⁰⁾的原子d_k各元素为：d_k(i)表示原子d_k中的第i行元素，表示稀疏分解系数矩阵X的第k行的第i列元素。

步骤2.2.2：确定误差矩阵E_k：其中，d_j为字典D的第j列原子，为稀疏分解系数矩阵X的第j行，E_k表示不使用D的第k列原子进行稀疏分解所产生的误差矩阵。

步骤2.2.3：将E_k乘以矩阵R_k得到新的误差矩阵 其中R_k为N×|d_k|的矩阵，|d_k|为d_k的模；

步骤2.2.4：对误差矩阵进行奇异值分解：Δ为所要求取的奇异值矩阵，U和V^T为两个正交矩阵，用正交矩阵U的第一列u₁更新目标训练字典D的第k列原子d_k，即d_k＝u₁,Δ(1,1)表示分解后奇异值矩阵中最大的奇异值，v₁表示给奇异值Δ(1,1)对应的特征向量，为由最大奇异值重构后的特征向量。

步骤2.2.5：重复步骤2.2.2至步骤2.2.4，对目标训练字典D^(J-1)中的每列原子进行更新处理，更新完毕后得到新的字典D^(J+1)。

步骤2.3：更新迭代次数J＝J+1；然后判断迭代次数J是否达到设定的最大迭代次数，若是，停止学习，输出目标优化字典D^(J-1)以及相对应的稀疏分解系数矩阵X，否则，判断当前得到的目标优化字典D^(J-1)是否满足总优化公式的约束条件，若是，输出目标优化字典D^(J-1)以及相对应的稀疏分解系数矩阵X，则转步骤2.1继续执行，继续进行字典学习。

步骤2.4：对样本图像Y进行去除噪声处理，得到新图像Y″＝D^(J-1)X。

步骤3：采用质心方法对去除噪声后的训练样本或识别样本进行目标识别区域提取。即提取样本图像的感兴趣区域(RegionofInterest，ROI)，其方法为：以感兴趣目标区域的质心(x₀,y₀)为中心，分别向四周扩展出一个L×L的矩阵，L为整数。目标的质心位置(x₀,y₀)可用一阶矩和零阶矩的比值表示：

$x_0=\frac{m_{10}}{m_{00}},y_0=\frac{m_{01}}{m_{00}}$

m_pq表示p+q阶原点矩，定义为：m_pq＝∑_x∑_yx^py^qf(x,y)

其中，(x,y)表示样本图像的像素坐标，而f(x,y)表示图像分割后的目标ROI区域上的象素幅度值，经过质心方法提取后，输出的具备目标识别区域的样本图像大小为L×L。

步骤4：对所有经过步骤3处理后的训练样本图像或识别样本图像，采用高斯核函数将样本图像映射到高维特征空间，得到特征空间中的训练样本集或识别样本集。所采用的高斯核函数的计算公式为：

k_G＝exp(-|Y_u-Y_v|²/(2*σ_G ²))

其中，Y_u、Y_v表示样本集中不同的两个样本图像，k_G为核映射结果，σ_G为高斯核函数的宽度参数。

步骤5：采用二维Gabor滤波器，对经步骤4处理后得到的训练样本图像或识别样本图像进行特征提取。如图4所示，为经过Gabor滤波后不同方向不同尺度下的特征图。Gabor滤波器的多尺度及方向性滤波能够较好地克服由于SAR图像的方位角变化带来的识别困难，最大限度地改善样本在特征空间的分布，为后续的识别工作提供极大的便利，二维Gabor滤波G(u,v)定义为：

G(u,v)＝∫∫f(x,y)g(u-x,v-y)dxdy

(u,v)表示对坐标(x，y)进行Gabor空域变化的坐标，g(u-x,v-y)通过代入二维Gabor滤波器核函数的得到，二维Gabor滤波器核函数具体形如下：

$g(x,y,\mathrm{\λ},{\mathrm{\δ}}_x,{\mathrm{\σ}}_y)=\exp [-\frac{1}{2}\{\frac{R_1^2}{{\mathrm{\σ}}_x^2}+\frac{R_2^2}{{\mathrm{\σ}}_y^2}\left\}\right]\exp \left\{i\frac{{2\mathrm{\πR}}_1}{\mathrm{\λ}}\right\}$

其中，参数R₁=xcosθ_k+ysinθ，参数R₂=-xsinθ_k+ycosθ_k，θ表示Gabor函数的方位角，代表一个振荡函数，实部为余弦函数，g(x,y,λ，σ_x,σ_y)虚部为正弦函数。是高斯函数，σ_x，σ_y分别是Gauss函数包络线沿x轴，y轴的标准方差，λ代表Gauss窗的波长,θ_k代表Gauss窗的振荡方向，通过它可以控制采样的方向。将两者相乘实际上是通过加窗限制了振荡函数的范围，使其只在局部有效。由高斯函数局部性可知，二维Gabor滤波器实际上是在(x，y)处抽取特征。

步骤6：采用局部相位量化（LPQ）方法对SAR训练样本或识别样本再次进行特征提取，如图5所示，（a）表示原图，（b）～（f）表示对（a）分别在LPQ子窗口大小为13、23、33、43和53下提取得到的特征图。具体步骤6特征提取的步骤为：

步骤6.1：进行短时傅里叶变换(STFT)。局部相位量化（LPQ）方法基于傅里叶相位谱的模糊不变性，它利用局部相位信息，首先对经步骤4处理后的样本图像上的每一个像素点，获取该像素点M×M邻域N_x上的局部相位信息，其计算公式如下：

$F(u^{\′},v^{\′})={\mathrm{\Σ}}_{y\∈N_x}f(x-y)e^{-j2{\mathrm{\πu}}^{\′T}y}=w_{u^{\′}}^T{f}_x$

F(u',v')表示经过短时傅立叶变换后得到的局部相位信息，(x,y)表示样本图像中像素点的坐标，(u',v')表示对样本图像中的坐标(x，y)进行短时傅里叶变换后的坐标，u'表示频率，u′^T表示u'的转置，f(x-y)表示样本图像f(x)的卷积；w_u'是二维离散傅里变换（DFT）在频率u'的基本向量，是w_u'的转置；f_x是包含N_x里面所有M²个图像样本的向量。

M为整数，本发明中对于分辨率为48×48的感兴趣目标区域，M的优选值为21，23或25，根据优选值进行STFT处理后，能得到有效的特征图像。

在LPQ方法里仅需考虑四个复值系数，对应四个二维频率u'₁～u'₄：

u'₁＝[a,0]^T,u'₂＝[0,a]^T,u'₃＝[a,a]^T,u'₄＝[a,-a]^T

其中a是一个标量频率。

设四个频率u'₁～u'₄对应的局部相位信息为F(u'₁,v'),F(u'₂,v'),F(u'₃,v′),F(u'₄,v′)，则四个局部相位信息组成矩阵为： $F_x^c=[F({u^{\′}}_1,v^{\′}),F({u^{\′}}_2,v^{\′}),F({u^{\′}}_3,v^{\′}),F({u^{\′}}_4,v^{\′})],$ 则根据集合可得到转换系数矩阵F_x：其中，Re{*},Im{*}分别表示一个复数的实部和虚部。

矩阵F_x对应的8*M²转换矩阵W是： $W=\left[\mathrm{Re}\right\{W_{{u^{\′}}_1},W_{{u^{\′}}_2},W_{{u^{\′}}_3},W_{{u^{\′}}_4}\},\mathrm{Im}\{W_{{u^{\′}}_1},W_{{u^{\′}}_2},W_{{u^{\′}}_3},W_{{u^{\′}}_4}\left\}\right],$ 其中， 分别表示二维离散傅里变换在频率u'₁,u'₁,u'₁,u'₄的转换基矩阵。

则有F_x＝Wf_x。

步骤6.2：系数统计分析，获得转换系数矩阵F_x的协方差矩阵。设图像函数f(x)是一阶马尔可夫过程，临近像素值之间的相关系数是ρ，样本方差是σ²，为不失一般性，设σ²＝1。因此两个相邻像素点x_i和x_j的协方差σ_i,j表示为：

${\mathrm{\σ}}_{i,j}={\mathrm{\ρ}}^{\left|\right|x_i-x_j\left|\right|},i\≠j,$

||*||表示L₂模，在领域N_x内M²个像素点的协方差矩阵可表达成：

因此，获得转换系数矩阵F_x的协方差矩阵D_x为：D_x＝WCW^T。

步骤6.3：解相关与量化。在量化之前系数是解相关的，可以证明，如果量化的样本是统计独立的，那么该信息能够被最大限度的保存在标量量化中。使用高斯白变换将转换系数矩阵F_x达到高斯独立分布，得到的G_x：G_x=V^TF_x，V^T是由矩阵D_x进行奇异值分解得出来的规范正交矩阵，D_x＝U∑V^T，U和V^T分别称为前奇异值正交矩阵和后奇异值正交矩阵，∑为分解的奇异值矩阵。

获取所有像素点的高斯独立分布矩阵G_x，例如：对于像素点x∈{x₁,x₂,…,x_N}，分别得到一个对应的矩阵G_x，然后对得到的所有矩阵G_x进行量化，使用一个简单的量化规则对结果向量进行量化，量化值q_j如下：

g_i是G_x第j个元素。量化的系数用0-255之间的二进制编码整数值b来代表：

$b={\mathrm{\Σ}}_{j=1}^8{q}_j{2}^{j-1}$

最后，生成一个所有像素点编码整数值的直方图，而该直方图就是后续分类中使用的一个256维特征向量。

步骤7：特征融合。针对每个样本图像，将步骤5采用Gabor滤波器得到的特征图和步骤6采用LPQ方法得到的特征图，进行加权融合，形成对应的样本图像的融合特征图。

步骤8：根据超完备训练样本集中的各样本图像得到的融合的特征图，构建超完备训练样本集中所有样本图像的仿生流形，如图6所示，圆圈和三角代表不同的两类样本集，图6示意了不同类的样本在高维空间中构建出来的仿生流形。仿生流形的构建方法如下：

步骤8a，在同类样本具有连续性的前提下，对同类样本在高维空间分布的流形进行维数分析；若某类样本在高维空间分布的轨迹是线性的，则为一维流形；若其分布轨迹是面，则为二维流形，以此类推。

步骤8b，合理选择和排列训练样本。

步骤8c，由训练样本流形维数及经过合理选择和排列后的训练样本确定覆盖该类样本子空间的近似几何形体B；

步骤8d，根据步骤8c所确定的近似几何形体B来构造该类封闭、有界的子空间；

步骤8e，重复步骤8c，直至所有的训练样本都完成子空间的构造。

优选地，本发明所有训练样本的仿生流形采用”超香肠”神经元模型，其构建具体步骤如下：

步骤8.1：设在某一类训练样本集P中有n个样本图像P＝{ρ₁,ρ₂,…ρ_n}，先找出这n个样本中距离最近的两个样本p_i和p_j，即ρ(p_i,p_j)≤ρ(p_x,p_y)；其中：p_x和p_y为训练样本集P中的任意两个样本，ρ(·)代表两个样本间的欧式距离函数。

步骤8.2：将两个样本p_i和p_j用一条直线连接起来，构成一个“超香肠”神经元，依次计算P中余下的n-2个点到该“超香肠”神经元的距离，若存在某一样本到该“超香肠”神经元的距离小于半径阈值r₀，则认为该样本是坏点并剔除出去。阈值r₀的最优值要通过实验来选取。

步骤8.3：确定覆盖该类样本子空间的近似几何形体B。剔除坏点后，根据以下公式来确定覆盖该类样本子空间的近似几何形体B：

$f_{\mathrm{HSN}}\left(x\right)=\mathrm{sgn}(2^{\frac{d^2(x,\stackrel{\‾}{x_1{x}_2})}{r_0^2}}-0.5)$

其中， $d^2(x,\stackrel{\&OverBar;}{x_1{x}_2})=\left\{\begin{array}{c}{\left|\right|x-x_1\left|\right|}^2,q(x,x_1,x_2)<0\\ {\left|\right|x-x_2\left|\right|}^2,q(x,x_1,x_2)>\left|\right|x-x_2\left|\right|\\ {\left|\right|x-x_1\left|\right|}^2-q^2(x,x_1,x_2),\mathrm{otherwise}\end{array}\right.$

$q(x,x_1,x_2)=<{x-x}_1,\frac{x_2-x_1}{\left|\right|x_2-x_1\left|\right|}>$

f_HSN(x)表示椭圆神经元函数，x为训练样本集P中的某个样本，x₁,x₂为步骤8.1中在训练样本集P中找到的距离最近的两个样本，即x₁,x₂也就是p_i和p_j。此处，x，x₁，x₂都表示融合后的特征点。d²表示样本x到样本x₁与样本x₂所构成的线段的高维距离。q函数用于判断样本x位于样本x₁与样本x₂所构成的椭圆神经元的左侧、右侧还是上方。

步骤8.4：根据近似几何形体B构造该类的子空间，从训练样本集Y中选取j个样本，形成新的样本集Y，训练样本集Y和训练样本集P为同一类，其表达式如下：

B＝{x|x＝Y_i'(i＝0,1,2,...,j-1),ρ(Y_i',Y_i+1′)≤ξ≤ρ(Y_i-1′,Y_i+1′))},

其中，ξ为一个无穷小的数值，Y_i'表示样本集Y’中的第i个样本。

用j个神经元近似地覆盖样本子空间：

第i个覆盖神经元P_i为：

P_i＝{x|＝ρ(x,y)≤r₀,y∈B_i,x∈Rⁿ}

B_i＝{x|x＝αY_i'+(1-α)Y_i+1′,α∈(0,1)}

α为0到1范围内的参数，Rⁿ表示n维实数空间，B_i表示第i个覆盖神经元的近似几何体。

则覆盖所有的j个神经元的子空间P_a′为：

$P_a^{\&prime;}={\&cup;}_{i=0}^{j-1}{P}_i$

步骤8.5：重复步骤8.2、步骤8.3和步骤8.4，直至该类内所有样本都已学习完毕。

步骤9：对识别样本进行识别。将识别样本图像输入步骤2中，开始进行处理，经过步骤2～步骤7的处理，得到识别样本图像的融合特征图，对得到的融合特征图进行识别。

本发明采用仿生模式识别方法进行SAR目标图像识别，本发明采用仿生模式识别方法进行SAR目标图像识别，为了简化计算过程和降低计算复杂度，首先计算识别样本中样本到各个单形的距离，而点到覆盖流形的距离定义为点到各单形距离中的最小值，其具体步骤如下：

首先，采用点到单形距离的算法计算出每个样本到各“超香肠”链的距离，对于第i类待识别样本中的第j个样本p_ij确定该样本到各“超香肠”链的距离d_ij。

其次，找出d_ij中的距离最小值d_min，即该样本p_ij到各“超香肠”链距离的最小值。设此时样本p_ij与第I个“超香肠”链距离最近。

最后，对于待识别样本p_ij若该距离最小值d_min小于等于半径阈值r₀，则认为此样本p_ij属于第I类样本，再若i＝I，则认为此次识别正确；若i≠I，则认为此次识别错误。若该最小值大于半径阈值r₀，则认为此次错误，待识别样本p_ij被拒绝识别。对于拒识样本，若该距离最小值d_min大于小于半径阈值r₀，则认为此次识别错误，若该距离最小值d_min大于半径阈值r₀，则认为此次正确拒绝。

Claims (9)

1.一种基于样本产生与核局部特征融合的SAR图像仿生识别方法，其特征在于，通过以下步骤实现：

步骤1：构建超完备训练样本集，具体是：对输入的每一类SAR的原始训练图像分别采用样本旋转的方法产生虚拟训练样本图，由SAR的原始训练图像和所产生的虚拟训练样本图构成最终的超完备训练样本集，将超完备训练样本集输入步骤2中；

步骤2：对样本集中的每个样本图像，采用K-SVD字典学习算法进行噪声处理；

步骤3：对噪声处理后的各样本图像，采用质心方法提取目标识别区域；

步骤4：采用高斯核函数将所有经步骤3输出的样本图像映射到高维特征空间；

步骤5：采用二维Gabor滤波器对步骤4得到的样本图像进行特征提取，得到样本图像的特征图；

步骤6：采用局部相位量化方法对步骤4得到的样本图像进行特征提取，得到样本图像的特征图；

步骤7：将步骤5和步骤6得到的同一样本图像的特征图进行加权融合，得到融合的特征图；

步骤8：根据超完备训练样本集中的各样本图像的融合特征图，构建超完备训练样本集中所有样本图像的仿生流形；所构建的仿生流形采用”超香肠”神经元模型；

步骤9：对待识别样本集中的各样本图像输入步骤2中，经过步骤2至步骤7的处理，得到待测试的样本图像的融合特征图，然后采用仿生模式识别方法进行SAR目标图像识别。

2.根据权利要求1所述的SAR图像仿生识别方法，其特征在于，步骤1中所述的采用样本旋转的方法产生虚拟训练样本图，具体是：

步骤1.1：确定旋转角度θ₁；

步骤1.2：确定每一类SAR的原始训练图像的虚拟训练样本图的数量为：(360°/θ₁)-1；

步骤1.3：将SAR的原始训练图像每旋转一次角度θ₁，得到对应旋转角度下的虚拟训练样本图。

3.根据权利要求1所述的SAR图像仿生识别方法，其特征在于，所述的步骤2具体是：设Y为输入的样本集中一幅原始样本图像，将样本图像Y分割成N块，则有y_i表示样本图像Y分割后的第i块，则通过K-SVD字典学习算法，寻找最优的目标训练字典D和最优的稀疏分解系数矩阵X来表示样本图像Y，设随机初始化的字典D⁽⁰⁾为N×K的矩阵，D⁽⁰⁾为当前的目标训练字典，令迭代次数J初始为1，然后进行下面步骤：

步骤2.1：优化每个图像小块y_i对应的稀疏分解局部系数得到更新的稀疏分解系数矩阵X；具体优化稀疏分解局部系数的方法是：保持当前的目标训练字典D^(J-1)不变，采用追踪算法循环计算最优稀疏分解局部系数优化公式如下：

$\underset{x_i}{\min }\left\{\right||y_i-{\mathrm{Dx}}_i^r|{|}_2^2\},\mathrm{su}bject\mathrm{to},{\&ForAll;}_i,|\left|x_i^r\right|{|}_0\&le;T_0,(i=1,2,...,N)$

subjectto表示约束条件，表示任意i，表示稀疏分解系数矩阵X的第i列，||·||₀为向量的0范数，T₀为稀疏度控制阈值，D为当前的目标训练字典D^(J-1)；

步骤2.2：更新目标训练字典，具体是：保持步骤2.1优化得到的图像小块y_i的稀疏分解局部系数保持不变，对图像小块y_i进行循环学习，更新当前的目标训练字典D^(J-1)的第k列原子d_k，k＝1,2,3,…,K，具体过程为：

步骤2.2.1：获取当前目标训练字典D^(J-1)的原子d_k；

字典D⁽⁰⁾的第k列原子d_k的第i行元素d_k(i)为：表示稀疏分解系数矩阵X的第k行第i列元素；

步骤2.2.2：确定误差矩阵E_k：其中，d_j为字典D的第j列原子，为稀疏分解系数矩阵X的第j行向量；

步骤2.2.3：将误差矩阵E_k乘以矩阵R_k得到新的误差矩阵其中R_k为N×|d_k|的矩阵，|d_k|为d_k的模；

步骤2.2.4：对误差矩阵进行奇异值分解：Δ为所要求取的奇异值矩阵，U和V^T为两个正交矩阵，用正交矩阵U的第一列原子u₁更新第k列原子d_k：d_k＝u₁；

步骤2.2.5：重复步骤2.2.2至步骤2.2.4，对目标训练字典D^(J-1)中的每列原子都进行更新，更新完毕后得到新的目标训练字典D^(J)；

步骤2.3：更新迭代次数J＝J+1；然后判断迭代次数J是否达到设定的最大迭代次数，若是，停止学习，输出目标优化字典D^(J-1)以及相对应的稀疏分解系数矩阵X；若否，判断当前得到的目标优化字典D^(J-1)是否满足总优化公式的约束条件，若是，输出目标优化字典D^(J-1)以及相对应的稀疏分解系数矩阵X，否则转步骤2.1继续执行；

总优化公式为： $\underset{D,X}{\min }\left\{\right||Y-DX|{|}_F^2\},\mathrm{su}bject\mathrm{to},{\&ForAll;}_i,|\left|x_i^r\right|{|}_0\&le;T_0,(i=1,2,...,N);$ F表示F范数，D为当前得到的目标优化字典D^(J-1)；

步骤2.4：对样本图像Y进行去除噪声处理，得到新图像Y"＝D^(J-1)X。

4.根据权利要求1所述的SAR图像仿生识别方法，其特征在于，所述的步骤3具体实现方法是：以样本图像中感兴趣目标区域的质心(x₀,y₀)为中心，分别向四周扩展出一个L×L的矩阵；目标的质心位置(x₀,y₀)表示为：

m_pq表示p+q阶原点矩，p＝0,1；q＝0,1；m_pq＝∑_x∑_yx^py^qf(x,y)，其中，f(x,y)表示感兴趣目标区域上的象素幅度值，(x,y)表示图像的像素坐标；采用质心方法提取的目标识别区域为L×L大小的图像。

5.根据权利要求1所述的SAR图像仿生识别方法，其特征在于，所述步骤4采用的高斯核函数为：k_G＝exp(-|Y_u-Y_v|²/(2*σ_G ²))，k_G为核映射结果，Y_u、Y_v表示同一样本集中不同的两个样本图像；σ_G为高斯核函数的宽度参数。

6.根据权利要求1所述的SAR图像仿生识别方法，其特征在于，所述的步骤5中采用的二维Gabor滤波器G(u,v)为：

G(u,v)＝∫∫f(x,y)g(u-x,v-y)dxdy

其中，f(x,y)表示感兴趣目标区域上的象素幅度值，(x,y)表示样本图像的像素坐标，(u,v)表示对坐标(x,y)进行Gabor空域变化后的坐标，g(u-x,v-y)代入二维Gabor滤波器核函数的得到，二维Gabor滤波器核函数为：

$g(x,y,\mathrm{\&lambda;},{\mathrm{\&sigma;}}_x,{\mathrm{\&sigma;}}_y)=\exp \&lsqb;-\frac{1}{2}\{\frac{R_1^2}{{\mathrm{\&sigma;}}_x^2}+\frac{R_2^2}{{\mathrm{\&sigma;}}_y^2}\}\&rsqb;\exp \left\{i\frac{2{\mathrm{\&pi;R}}_1}{\mathrm{\&lambda;}}\right\}$

其中，参数R₁＝xcosθ_k+ysinθ，参数R₂＝-xsinθ_k+ycosθ_k，表示一个振荡函数，实部为余弦函数，虚部为正弦函数；是高斯函数，σ_x和σ_y分别是该高斯函数包络线沿x轴和y轴的标准方差，λ表示高斯窗的波长，θ_k表示高斯窗的振荡方向，θ表示Gabor函数的方位角。

7.根据权利要求1所述的SAR图像仿生识别方法，其特征在于，所述的步骤6的具体实现步骤是：

步骤6.1：进行短时傅里叶变换，对于样本图像中的每一个像素点(x,y)，获取该像素点在大小为M×M的邻域N_x上的局部相位信息F(u',v')：

$F(u^{\&prime;},v^{\&prime;})={\mathrm{\&Sigma;}}_{y\&Element;N_x}f(x-y)e^{-j2{{\mathrm{\&pi;u}}^{\&prime;}}^{T}y}=w_{u^{\&prime;}}^T{f}_x$

(u',v')表示对坐标(x,y)进行短时傅里叶变换后的坐标，u'表示频率，u^'T表示u'的转置，f(x-y)表示样本图像函数f(x)的卷积，w_u'是二维离散傅里变换在频率u'的基本向量，是w_u'的转置；f_x是包含N_x里面所有M²个像素点的向量，M为整数；

设a是一个标量频率，对于四个频率u'₁＝[a,0]^T,u'₂＝[0,a]^T,u'₃＝[a,a]^T,u'₄＝[a,-a]^T对应的四个局部相位信息分别为：F(u'₁,v'),F(u'₂,v'),F(u'₃,v'),F(u'₄,v')，则四个局部相位信息组成矩阵为： $F_x^c=\&lsqb;F({u^{\&prime;}}_1,v^{\&prime;}),F({u^{\&prime;}}_2,v^{\&prime;}),F({u^{\&prime;}}_3,v^{\&prime;}),F({u^{\&prime;}}_4,v^{\&prime;})\&rsqb;,$ 进一步根据集合可得到转换系数矩阵F_x：其中，Re{*},Im{*}分别表示一个复数的实部和虚部；转换系数矩阵F_x对应的8*M²转换矩阵W是： $W=\&lsqb;\mathrm{Re}\{W_{{u^{\&prime;}}_1},W_{{u^{\&prime;}}_2},W_{{u^{\&prime;}}_3},W_{{u^{\&prime;}}_4}\},\mathrm{Im}\{W_{{u^{\&prime;}}_1},W_{{u^{\&prime;}}_2},W_{{u^{\&prime;}}_3},W_{{u^{\&prime;}}_4}\}\&rsqb;,$ 则有F_x＝Wf_x；

步骤6.2：获得转换系数矩阵F_x的协方差矩阵，设样本图像函数f(x)是一阶马尔可夫过程，临近像素值之间的相关系数是ρ，样本方差是σ²，设σ²＝1，则在领域N_x内的两个相邻像素点x_i和x_j的协方差σ_i,j表示为：||*||表示L₂模；

则在领域N_x内M²个像素点的协方差矩阵C为：

因此，获得转换系数矩阵F_x的协方差矩阵D_x为：D_x＝WCW^T；

步骤6.3：对协方差矩阵D_x进行奇异值分解，得到后奇异值正交矩阵V^T，然后使用白变换将转换系数矩阵F_x达到高斯独立分布，得到矩阵G_x＝V^TF_x；获取所有像素点对应的矩阵G_x，然后对得到的所有矩阵G_x进行量化，量化公式为：其中，q_j表示量化值，g_j表示矩阵G_x中的第j个元素；

用0-255之间的二进制编码整数值b来代表量化的系数：

最后，生成一个样本图像中所有像素点对应的二进制编码整数值b的直方图，该直方图就是样本图像的特征图。

8.根据权利要求1所述的SAR图像仿生识别方法，其特征在于，步骤8中所述的仿生流形为”超香肠”神经元模型，具体构建方法是：

步骤8.1：设在某一类训练样本集P中有n个样本P＝{p₁,p₂,…p_n}，找出n个样本中距离最近的两个样本p_i和p_j；

步骤8.2：将样本p_i和p_j用一条直线连接起来，构成一个“超香肠”神经元，依次计算余下的n-2个样本到该“超香肠”神经元的距离，若存在某一样本到该“超香肠”神经元的距离小于半径阈值r₀，则剔除该样本；

步骤8.3：根据以下公式确定覆盖该类样本子空间的近似几何形体B：

$f_{HSN}\left(x\right)=\mathrm{sgn}(2^{\frac{d^2(x,\stackrel{\&OverBar;}{x_1{x}_2})}{r_0^2}}-0.5)$

其中， $d^2(x,\stackrel{\&OverBar;}{x_1{x}_2})=\left\{\begin{array}{c}\left|\right|x-x_1|{|}^2,q(x,x_1,x_2)<0\\ \left|\right|x-x_2|{|}^2,q(x,x_1,x_2)>||x-x_2||\\ \left|\right|x-x_1|{|}^2-q^2(x,x_1,x_2),otherwise\end{array}\right.$

$q(x,x_1,x_2)=<x-x_1,\frac{x_2-x_1}{\left|\right|x_2-x_1\left|\right|}>$

f_HSN(x)表示椭圆神经元函数，x为训练样本集P中的某个样本，x₁,x₂为步骤8.1中在训练样本集P中找到的距离最近的两个样本p_i和p_j；d²表示样本x到样本x₁与样本x₂所构成的线段的高维距离；q函数用于判断样本x位于样本x₁与样本x₂所构成的椭圆神经元的左侧、右侧还是上方；

步骤8.4：根据近似几何形体B构造该类样本的子空间，具体是：

从与训练样本集P为同一类的训练样本集Y中选取j个样本，形成新的样本Y’，其表达式如下：

$B=\left\{x\right|x=Y_i^{\&prime;}(i=0,1,2,...,j-1),\mathrm{\&rho;}(Y_i^{\&prime;},Y_{i+1}^{\&prime;})\&le;\mathrm{\&xi;}\&le;\mathrm{\&rho;}(Y_{i-1}^{\&prime;},Y_{i+1}^{\&prime;})\left)\right\},Y^{\&prime;}\&Subset;Y$

其中，ξ为一个无穷小的数值，Y_i'表示样本集Y’中的第i个样本，ρ(·)代表两个样本间的欧式距离函数；

第i个覆盖神经元P_i为：P_i＝{x|x＝ρ(x,y)≤r₀,y∈B_i,x∈Rⁿ}；

B_i表示第i个覆盖神经元的近似几何体：B_i＝{x|x＝αY'+(1-α)Y′_i+1,α∈(0,1)}；α为0到1范围内的参数，Rⁿ表示n维实数空间；

则覆盖所有的j个神经元的子空间P′_α为：

步骤8.5：重复步骤8.2、步骤8.3和步骤8.4，直至该类内所有样本都已学习完毕。

9.根据权利要求8所述的SAR图像仿生识别方法，其特征在于，步骤9所述的采用仿生模式识别方法进行SAR目标图像识别，具体方法是：首先，对于第i类待识别样本中的第j个样本p_ij，确定该样本到各“超香肠”链的距离，并找到其中距离最小值d_min，设样本p_ij距离第I个“超香肠”链最近；然后，判断距离最小值d_min是否小于等于半径阈值r₀，若是，则样本p_ij属于第I类样本，进一步判断i是否等于I，若是，则识别正确，否则，识别错误，若距离最小值d_min大于半径阈值r₀，则样本p_ij被拒绝识别，对于被拒绝识别的样本，再次判断距离最小值d_min是否大于半径阈值r₀，若是，则正确拒绝，否则，识别错误。