CN104680183A - 基于散射点和k中心一类分类器的sar目标鉴别方法 - Google Patents

Abstract

本发明公开了基于散射点和K中心一类分类器的SAR目标鉴别方法，主要解决了现有技术未从雷达成像的本质揭示目标和杂波虚警的区别，鉴别总正确率低的缺点。其技术方案是：1)对SAR图像进行恒虚警检测，提取切片；2)从提取的切片中选取包含真实目标的切片构成训练样本；3)依散射点模型从训练样本中提取散射点矩阵，并进行幅度2-范数归一；4)对训练样本散射点矩阵进行K中心聚类，得到簇中心；5)计算K中心一类分类器的鉴别门限Thr；6)计算测试样本和簇中心的最小双向Hausdorff距离；7)根据该距离判断测试样本是否为目标。本发明降低了杂波虚警的虚警率，有效提高了鉴别总正确率，适用于SAR图像中车辆等具有显著强散射点分布特征的人造目标的鉴别。

Description

基于散射点和K中心一类分类器的SAR目标鉴别方法

技术领域

本发明属于雷达技术领域，涉及目标的检测与鉴别方法，可用于SAR图像中车辆等具有显著的强散射点分布特征的人造目标的鉴别。

背景技术

雷达成像技术是在上世纪50年代发展起来的，在之后的60年里突飞猛进的发展，目前，已经在军事、农林、地质、海洋、灾害、测绘等诸多方面得到广泛的应用。合成孔径雷达SAR具有全天候、全天时、分辨率高以及穿透力强等特点，成为目前对地观测和军事侦察的重要手段。

SAR图像解译是当前SAR应用的前沿课题，同时也是近几年来对地观测技术应用的一个热门领域，具有重要的研究意义和广泛的应用前景。美国林肯实验室提出了SAR图像自动目标识别的三级处理流程图并被广泛使用。该流程包含三个基本阶段：检测、鉴别、分类。目标鉴别属于目标检测和目标识别的中间阶段，是对检测阶段得到的潜在目标区域进一步区分目标和杂波，去除杂波虚警，减小目标识别阶段的计算代价。自林肯实验室提出SAR图像目标鉴别的概念后，研究者们在鉴别特征提取、鉴别器设计、性能分析、硬件实现上做了大量的研究。

鉴别首先需要解决的问题是特征提取，提取的应该是能够揭示目标和杂波虚警的本质差异，使得目标和杂波虚警在特征空间上可以分离的特征。传统的SAR图像鉴别特征，如标准偏差、分形维数、加权填充比等从纹理方面揭示目标和杂波虚警的差异；如目标直径、归一化转动惯量、峰值CFAR、均值CFAR等从尺寸方面揭示目标和杂波虚警的差异；揭示目标像素集合的空间分布特征，如空间分布、拐点特征、加速度特征等从目标像素集合的空间分布方面揭示目标和杂波虚警的差异。这些特性对目标和杂波虚警的分离是有用的，但是这些特征没有从雷达相干成像的本质揭示目标和杂波虚警的区别，其鉴别总正确率较低。另外，部分传统SAR图像鉴别特征的提取，需要先对样本进行恒虚警检测，而恒虚警检测的好坏直接影响特征的可分性。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于散射点和K中心一类分类器的SAR目标鉴别方法，提高目标鉴别正确率。

实现本发明目的的技术方案是：从雷达相干成像的本质出发，提出散射点特征，利用该特征描述目标与杂波虚警的区别；根据不同样本中提取散射点个数的不同，引入双向Hausdorff距离计算点数不同的两个点集间的距离，并用K中心一类分类器对样本进行鉴别。其具体步骤包括如下：

(1)训练步骤：

(1a)对SAR图像I进行恒虚警检测，得到二值图像C，对该二值图像进行聚类，得到疑似目标区域，以疑似目标区域的重心为几何中心，截取m×n的切片x_i，i＝1,...,N，N是从SAR图像I中提取切片的总个数；

(1b)从切片{x₁,x₂,...,x_i,...,x_N}中选取包含真实目标的切片{y₁,y₂,...,y_j,...,y_Q}构成训练样本集，其中j＝1,...,Q，Q≤N，Q是训练样本个数；

(1c)估计训练样本y_j的信号能量比，以此为门限从训练样本y_j中提取散射点矩阵A_j，并对散射点矩阵A_j的幅度进行2-范数归一，2-范数归一后的散射点矩阵记为A'_j；

(1d)根据预先给定的中心个数K，K≤Q，用K中心一类分类器，将训练样本对应的散射点集{A'₁,A'₂,...,A'_j,...,A'_Q}聚为K簇，同时保存这K簇的中心t＝1,...,K；

(1e)根据预先给定的拒判率P和K簇的中心计算K中心一类分类器的门限Thr，其中0≤P＜1；

(2)测试步骤：

(2a)估计测试切片z的信号能量比，并依此为门限提取散射点矩阵B，对散射点矩阵B的幅度进行2-范数归一，2-范数归一后的散射点矩阵记为B'；

(2b)计算测试样本的散射点矩阵B'与步骤(1e)得到的K簇的中心之间的最小双向Hausdorff距离d_B；

(2c)将测试样本的散射点矩阵B'与K簇的中心之间的最小双向Hausdorff距离d_B与步骤(1e)得到的K中心一类分类器的判决门限Thr进行比较，如果d_B≤Thr，则测试样本z为目标，否则，测试样本z为杂波虚警。

本发明从雷达成像的本质出发提取散射点特征，更好的揭示了目标和杂波虚警的区别。与现有技术相比，能够有效的提高总正确率，同时降低了杂波虚警的虚警率。

附图说明

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

图2是本发明中的恒虚警检测子流程图；

图3是本发明中的空心滑窗示意图；

图4是本发明实验使用的SAR图像；

图5是用本发明和传统方法对图4进行鉴别的结果对比图。

具体实施方式

以下结合附图对本发明技术方案及效果作进一步详细表述。

参照图1，本发明的实现包括训练阶段和测试阶段，详细步骤如下：

一、训练阶段

步骤1，对SAR图像I进行恒虚警检测，得到二值图像C。

参照图2，本步骤的具体实现如下：

(1a)假设目标的最大长度为l，雷达的距离分辨率和方位分辨率分别为ρ_r、ρ_a，则目标在SAR图像中的最大长度l_r和最大宽度l_a分别为：和选取警戒窗长度 $L>2\sqrt{l_r^2+l_a^2};$

(1b)给定杂波区宽度N_r，虚警概率P_fa；

(1c)选择双参数恒虚警检测器作为恒虚警检测器，根据高斯杂波统计模型和虚警概率P_fa计算出检测阈值T₁＝Φ^-1(1-P_fa)，其中Φ^-1(·)表示标准正态分布的反函数；

(1d)对SAR图像I中像素的幅度进行如下变化：

b＝10*log₁₀(a²)

其中a表示图像I中某一像素的幅度，b表示图像I中某一像素的幅度变化后的幅度，变化后的图像用I₁表示；

(1e)选取I₁(i₁,j₁)单元作为待检测单元，以待检测单元I₁(i₁,j₁)为中心，按照图3所示选取警戒区和杂波区Ω_c，用N_c表示杂波区域像素的个数，其中(i₁,j₁)表示待检测单元的坐标；

(1f)估计杂波区的均值方差其中c_w是杂波区域第w个单元的的幅度；

(1g)用矩阵C标记图像I中对应位置是否是目标像素，设待检测单元I₁(i₁,j₁)处的幅度为t，若则I₁(i₁,j₁)是目标像素，矩阵C中对应位置标记为1，否则，I₁(i₁,j₁)是非目标像素，矩阵C中对应位置标记为0；

(1h)对变化后图像I₁中的所有像素重复(1e)—(1g)，得到二值图像C。

步骤2，对二值图像C进行聚类，并进行切片提取。

(2a)对二值图像C进行区域合并，若C中任意两个区域重心的距离小于目标的最大长度l的一半，则合并区域，否则保留区域，区域合并后的二值图像记为C₁；

(2b)以C₁中某一区域的重心为中心，在I中对应位置截取m×n的切片x_i，遍历C₁中的所有区域，得到N个切片，其中i＝1,...,N。

步骤3，从切片{x₁,x₂,...,x_i,...,x_N}中选取包含真实目标的切片{y₁,y₂,...,y_j,...,y_Q}构成训练样本集，其中j＝1,...,Q，Q≤N，Q是训练样本个数。

步骤4，估计训练样本y_j的信号能量比，以此为门限从训练样本y_j中提取散射点矩阵A_j。

(4a)对切片y_j进行恒虚警检测，得到切片y_j的二值图像y_jb，y_jb中被置为1的为目标像素y_jt，置为0的为非目标像素y_jc；

(4b)计算切片y_j的信号能量比，其中⊙表示矩阵的Hadamard积，vec(·)表示将矩阵向量化；

(4c)对切片y_j做2维快速逆傅立叶变换并向量化得到原始信号ξ＝vec（ifft2(y_j)），初始化信号残差r＝ξ，其中ifft2(·)表示对信号进行2维快速逆傅立叶变换；

(4d)根据切片y_j的大小m×n，以及雷达的距离分辨率和方位分辨率ρ_r、ρ_a，假设距离维对应的单元个数为m，方位维对应的单元个数为n，则散射点的位置坐标范围为X＝[0:ρ_a:nρ_a]，Y＝[0:ρ_r:mρ_r]；

(4e)根据散射点模型 $E(f,\mathrm{\ψ})=\∂\·\exp \left(\frac{-\mathrm{\υ}4\mathrm{\πf}}{c}(x\cos \mathrm{\ψ}+y\sin \mathrm{\ψ})\right),$ 构造字典原子d_i'j'：

$d_{i^{\′}{j}^{\′}}=\frac{d_{i^{\′}{j}^{\′}}^{\′}}{{\left|\right|d_{i^{\′}{j}^{\′}}^{\′}\left|\right|}_2},$

$d_{i^{\′}{j}^{\′}}^{\′}=\mathrm{vec}\left(\exp \left(\frac{-\mathrm{\υ}4\mathrm{\πf}}{c}(X\left(i^{\′}\right)\cos \mathrm{\ψ}+Y\left(j^{\′}\right)\sin \mathrm{\ψ})\right)\right),$ i'＝1,...,N_x，j'＝1,...,N_y，

其中f、ψ分别是雷达的频率和方位角，exp(·)表示自然指数函数，||·||₂表示向量的2-范数，(·)^T表示转置，N_x是是散射点的位置坐标X的维数，N_y是散射点的位置坐标Y的维数，是某散射点的幅度，x、y分别是某散射点的距离维坐标和方位维坐标，υ是虚数单位，π是圆周率，c是光速，X(i')是X坐标的第i'个元素，Y(j')是Y坐标的第j'个元素；

(4f)由(4e)中的d_i'j'，构造字典 $D=[d_{\mathrm{1,1}},d_{\mathrm{1,2}},...,d_{1,N_y},...,d_{i^{\′}{j}^{\′}},...,d_{N_x,1},d_{N_x,2},...,d_{N_x,N_y}];$

(4g)用空集θ作为参数集合，空集S作为原子集合；

(4h)将字典D和信号残差r做内积，κ＝D·r，内积κ中的最大元素的序号为M，将D(:,M)加入到原子集合S，该列对应的参数(1 X(i') Y(j'))^T加入到参数集合θ，i'、j'满足i'×(N_y-1)+j'＝M，其中D(:,M)表示字典D中的第M列，X(i')是X坐标的第i'个元素，Y(j')是Y坐标的第j'个元素；

(4i)更新散射点的幅度更新重构信号更新信号残差更新重构能量比用替换参数集合θ的第一行，更新参数散射点的幅度，其中表示穆尔-彭罗斯(Moore-Penrose)逆矩阵，(·)^H表示共轭转置；

(4j)若重构能量比δ小于信号能量比R，则执行(4h)—(4i)，否则执行(4k)；

(4k)将参数集合θ赋值到散射点矩阵A_j，则 $A_j=\left[\begin{array}{ccccc}{a}_{j_1}& ...& a_{j_{i^{\′\′}}}& ...& a_{j_{N_j}}\\ x_{j_1}& ...& x_{j_{i^{\′\′}}}& ...& x_{j_{N_j}}\\ y_{j_1}& ...& y_{j_{i^{\′\′}}}& ...& y_{j_{N_j}}\end{array}\right],$ 其中分别表示散射点的幅度、距离维坐标、方位维坐标，i”＝1,...,N_j，N_j表示从y_j中提取的散射点的个数；

(4l)对j＝1,...,Q重复执行(4a)—(4k)，得到散射点集{A₁,A₂,...,A_j,...,A_Q}。

步骤5，对散射点矩阵A_j的幅度进行2-范数归一，2-范数归一后的散射点矩阵记为A'_j。

(5a)对散射点矩阵A_j的幅度i”＝1,...,N_j，按如下公式进行2-范数归一：

${\stackrel{\‾}{a}}_{j_{i^{\′\′}}}=\frac{a_{j_{i^{\′\′}}}}{{||\begin{array}{c}{\left(\begin{array}{ccccc}{a}_{j_1}& ...& a_{j_{i^{\′\′}}}& ...& a_{j_{N_j}}\end{array}\right)}^T\end{array}||}_2},$

(5b)将2-范数归一化后的散射点幅度代入散射点矩阵A_j的幅度得到归一化后的散射点矩阵A'_j为：

$A_j^{\′}=\left[\begin{array}{ccccc}{\stackrel{\‾}{a}}_{j_1}& ...& {\stackrel{\‾}{a}}_{j_{i^{\′\′}}}& ...& {\stackrel{\‾}{a}}_{j_{N_j}}\\ x_{j_1}& ...& x_{j_{i^{\′\′}}}& ...& x_{j_{N_j}}\\ y_{j_1}& ...& y_{j_{i^{\′\′}}}& ...& x_{j_{N_j}}\end{array}\right],j=1,...,Q;$

步骤6，根据预先给定的中心个数K，K≤Q，用K中心一类分类器，将训练样本对应的散射点集{A'₁,A'₂,...,A'_j,...,A'_Q}聚为K簇，并保存这K簇的中心t＝1,...,K。

(6a)将2-范数归一化后的散射点矩阵 $A_j^{\′}=\left[\begin{array}{ccccc}{\stackrel{\‾}{a}}_{j_1}& ...& {\stackrel{\‾}{a}}_{j_{i^{\′\′}}}& ...& {\stackrel{\‾}{a}}_{j_{N_j}}\\ x_{j_1}& ...& x_{j_{i^{\′\′}}}& ...& x_{j_{N_j}}\\ y_{j_1}& ...& y_{j_{i^{\′\′}}}& ...& x_{j_{N_j}}\end{array}\right]$ 分块表示为：

$A_j^{\′}=\left[\begin{array}{ccccc}{b}_{j_1}& ...& b_{j_p}& ...& b_{j_{N_j}}\end{array}\right],$

其中 $b_{j_p}={\left(\begin{array}{ccc}{\stackrel{\‾}{a}}_{j_p}& x_{j_p}& y_{j_p}\end{array}\right)}^T,$ p＝1,...,N_j，j＝1,...,Q；

(6b)计算训练样本散射点集的协方差矩阵：其中Σ_j是切片y_j的散射点矩阵A'_j的协方差矩阵，其表示如下：

${\mathrm{\Σ}}_j=\left[\begin{array}{ccc}{\mathrm{\σ}}_{j1}^2& 0& 0\\ 0& {\mathrm{\σ}}_{j2}^2& 0\\ 0& 0& {\mathrm{\σ}}_{j3}^2\end{array}\right],$

其中，分别是训练样本散射点矩阵 $A_j^{\′}=\left[\begin{array}{ccccc}{\stackrel{\‾}{a}}_{j_1}& ...& {\stackrel{\‾}{a}}_{j_{i^{\′\′}}}& ...& {\stackrel{\‾}{a}}_{j_{N_j}}\\ x_{j_1}& ...& x_{j_{i^{\′\′}}}& ...& x_{j_{N_j}}\\ y_{j_1}& ...& y_{j_{i^{\′\′}}}& ...& x_{j_{N_j}}\end{array}\right]$ 的第1行元素，第2行元素和第3行元素的方差；

(6c)随机从散射点集{A'₁,A'₂,...,A'_j,...,A'_Q}中选取K个样本作为簇中心计算代价函数 $J0=\underset{j=1,...,Q}{\max }\left(\underset{t=1,...,k}{\min }{\left|\right|A_j-A_{c_t}\left|\right|}_2^2\right);$

(6d)计算训练样本的簇标号 $L\left(j\right)=\underset{t=1,...K}{\arg \min }\mathrm{HD}(A_j^{\′},A_{c_t}),j=1,...,Q,$ 其中：表示训练样本A'_j和簇中心的双向Hausdorff距离，

$\mathrm{HD}(A_j^{\′},A_{c_t})=\max (h(A_j^{\′},A_{c_t}),h(A_{c_t},A_j^{\′}))$

$h(A_j^{\′},A_{c_t})=\frac{1}{N_j}\underset{p=1}{\overset{N_j}{\mathrm{\Σ}}}\underset{q=1,...,N_{c_t}}{\min }\left(d(b_{j_p},b_{c_{\mathrm{tq}}})\right)$

$h(A_{c_t},A_j^{\′})=\frac{1}{N_{c_t}}\underset{q=1}{\overset{N_{c_t}}{\mathrm{\Σ}}}\underset{p=1,...,N_j}{\min }\left(d(b_{c_{\mathrm{tq}}},b_{j_p})\right)$

$d(b_{j_p},b_{c_{\mathrm{tq}}})=d(b_{c_{\mathrm{tq}}},b_{j_p})=\sqrt{{(b_{j_p}-b_{c_{\mathrm{tq}}})}^T{\mathrm{\Σ}}^{-1}(b_{j_p}-b_{c_{\mathrm{tq}}})},$ (·)^-1表示矩阵的逆，表示第c_t个簇中心的散射点个数；

(6e)任意取一非中心样本备份并用代替第L(m₁)簇的中心计算代价函数若J1＜J0，则接受这种代替，将J1赋值给J0，若J1≥J0，则不接受这种代替，用备份的恢复

(6f)重复(6d)—(6e)直至簇中心不在发生变化。

(6g)保存这K簇的中心t＝1,...,K。

步骤7，根据预先给定的拒判率P和K簇的中心计算K中心一类分类器的门限Thr。

(7a)计算j＝1,...,Q，d_A(j)是第j个样本和簇中心间的距离；

(7b)根据预先给定的拒判率P和样本和簇中心间的距离，计算K中心一类分类器的门限 $\mathrm{Thr}=\left\{\begin{array}{cc}(d_A\left(\mathrm{num}\right)+d_A(\mathrm{num}+1)),& 1\&le;\mathrm{num}<Q\\ d_A\left(\mathrm{num}\right),& \mathrm{num}=Q\end{array}\right.,$ 表示向下取整，其中0≤P＜1。

二、测试阶段

步骤8，估计测试样本z的信号能量比，并依此为门限提取散射点矩阵B，对散射点矩阵B的幅度进行2-范数归一，得到2-范数归一后的散射点矩阵B'。

(8a)对待鉴别的SAR图像做恒虚警检测和聚类，得到疑似目标区域，以疑似目标区域的重心为中心从待鉴别的SAR图像中截取m×n的测试样本z；

(8b)从测试样本z提取散射点矩阵并归一，即重复步骤2—步骤5，得到测试样本z的2-范数归一化的散射点矩阵 $B^{\&prime;}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{a}}_{B_1}& ...& {\stackrel{\&OverBar;}{a}}_{B_{p^{\&prime;}}}& ...& {\stackrel{\&OverBar;}{a}}_{B_{N_B}}\\ x_{B_1}& ...& x_{B_{p^{\&prime;}}}& ...& x_{B_{N_B}}\\ y_{B_1}& ...& y_{B_{p^{\&prime;}}}& ...& y_{B_{N_B}}\end{array}\right],$ 其中p'＝1,...,N_B，N_B是从测试样本z中提取的散射点的个数。

步骤9，计算测试样本的散射点矩阵B'与K簇的中心之间的最小双向Hausdorff距离d_B。

(9a)将散射点矩阵B'分块表示为： $B^{\&prime;}=\left[\begin{array}{ccccc}{b}_{B_1}& ...& b_{B_{p^{\&prime;}}}& ...& b_{B_{N_B}}\end{array}\right],$ 其中

$b_{B_{p^{\&prime;}}}={\left[\begin{array}{ccc}{\stackrel{\&OverBar;}{a}}_{B_{p^{\&prime;}}}& x_{B_{p^{\&prime;}}}& y_{B_{p^{\&prime;}}}\end{array}\right]}^T,p^{\&prime;}=1,...,N_B;$

(9b)根据散射点矩阵B'的块和簇中心散射点矩阵的块计算单向Hausdorff距离和其计算公式如下：

$h\left(B^{\&prime;}{A}_{c_t}\right)=\frac{1}{N_B}\underset{p^{\&prime;}=1}{\overset{N_B}{\mathrm{\&Sigma;}}}\underset{q=1,...,N_{c_t}}{\min }\left(d(b_{B_{p^{\&prime;}}},b_{c_{\mathrm{tq}}})\right)$

$h(A_{c_t},B^{\&prime;})=\frac{1}{N_{c_t}}\underset{q=1}{\overset{N_{c_t}}{\mathrm{\&Sigma;}}}\underset{p^{\&prime;}=1,...,N_B}{\min }\left(d(b_{c_{\mathrm{tq}}},b_{B_{p^{\&prime;}}})\right)$

$d(b_{B_{p^{\&prime;}}},b_{c_{\mathrm{tq}}})=d(b_{c_{\mathrm{tq}}},b_{B_{p^{\&prime;}}})=\sqrt{{(b_{B_{p^{\&prime;}}}-b_{c_{\mathrm{tq}}})}^T{\mathrm{\&Sigma;}}^{-1}(b_{B_{p^{\&prime;}}}-b_{c_{\mathrm{tq}}})};$

(9c)根据单向Hausdorff距离和计算测试样本散射点矩阵B'与簇中心间的双向Hausdorff距离

$\mathrm{HD}(B^{\&prime;},A_{c_t})=\max (h(B^{\&prime;},A_{c_t}),h(A_{c_t},B^{\&prime;})),t=1,...,K;$

(9d)计算测试样本散射点矩阵B'与簇中心间的最小双向Hausdorff距离：

$d_B=\underset{t=1,...,K}{\min }\mathrm{HD}\left(B^{\&prime;}{A}_{c_t}\right).$

步骤10，将测试样本的散射点矩阵与K簇的中心之间的最小双向Hausdorff距离d_B与步骤(1e)得到的K中心一类分类器的判决门限Thr进行比较，如果d_B≤Thr，则测试样本z为目标，否则，测试样本z为杂波虚警。

下面结合附图和实验对本发明的效果做进一步的说明：

为了说明本发明的效果，将本发明与用林肯实验室提出的传统鉴别特征中的强度百分比恒虚警(PB，Percent Bright CFAR Feature)特征做鉴别特征、SVDD做一类分类器的传统鉴别方法作对比。

1.实验数据：

实验中采用的数据是美国桑尼迪实验室公开的MiniSAR数据中的两幅，如图4所示，其中图4(a)是郊区某停车场的SAR图像，图4(b)是某停机坪附近的SAR图像。

目标训练样本由对图4(a)做恒虚警检测得到的143幅切片中包含真实目标的51幅切片构成；目标测试样本由对图4(b)做恒虚警检测后到的117幅切片中包含真实目标的41幅切片构成；杂波虚警测试样本由对图4(a)做恒虚警检测得到的143幅切片中不包含真实目标的71幅切片和对图4(b)做恒虚警检测得到的117幅切片中不包含真实目标的45幅切片构成。

2.实验参数：

雷达参数：中心频率f_c＝16.8GHz，距离分辨率ρ_r＝0.1m，方位分辨率ρ_a＝0.1m。

恒虚警检测的参数：警戒窗L＝30，杂波宽度N_r＝1，虚警率P_fa＝0.01；

K中心一类分类器的参数：中心数K＝10，拒判率P＝2％。

3.实验内容与结果：

用本发明方法中的目标训练样本的散射点特征训练K中心一类分类器，然后用训练好的K中心一类分类器鉴别目标测试样本和杂波虚警测试样本；

用林肯特征中的强度百分比恒虚警特征做鉴别特征，用目标训练样本的强度百分比恒虚警特征训练SVDD一类分类器，然后用训练好的SVDD一类分类器鉴别目标测试样本和杂波虚警测试样本。鉴别结果如图5所示，其中：

图5(a)是用本发明方法对目标测试样本的鉴别结果；

图5(b)是用现有强度百分比恒虚警特征对目标测试样本的鉴别结果；

图5(c)是用本发明方法对图4(a)的杂波虚警测试样本的鉴别结果；

图5(d)是用现有强度百分比恒虚警特征对图4(a)的杂波虚警测试样本的鉴别结果；

图5(e)是用本发明方法对图4(b)的杂波虚警测试样本的鉴别结果；

图5(f)是用强度百分比恒虚警特征对图4(b)的杂波虚警测试样本的鉴别结果。

图5中的实线框表示鉴别正确的切片，虚线框表示鉴别错误的切片。

由图5可见，本发明的方法较现有方法的目标鉴别正确率高，杂波虚警率低，鉴别总正确率高。

对上述鉴别进行定量说明，如表1所示：

表1 鉴别结果

	目标正确率	杂波虚警率	总正确率
				强度百分比恒虚警特征+SVDD一类分类器	90.24％	22.41％	80.89％
本发明方法	92.68％	11.21％	89.81％

从上表中可以看出，本发明所提的方法，无论是在目标正确率，杂波虚警率和总正确率方面均有优势，特别是在杂波虚警率方面，本发明的优势更为明显，说明了本发明的有效性。

Claims (7)

1.一种基于散射点和K中心一类分类器的SAR目标鉴别方法，包括：

(1)训练步骤：

(1a)对SAR图像I进行恒虚警检测，得到二值图像C，对该二值图像进行聚类，得到疑似目标区域，以疑似目标区域的重心为几何中心，截取m×n的切片x_i，i＝1,...,N，N是从SAR图像I中提取切片的总个数；

(1b)从切片{x₁,x₂,...,x_i,...,x_N}中选取包含真实目标的切片{y₁,y₂,...,y_j,...,y_Q}构成训练样本集，其中j＝1,...,Q，Q≤N，Q是训练样本个数；

(1c)估计训练样本y_j的信号能量比，以此为门限从训练样本y_j中提取散射点矩阵A_j，并对散射点矩阵A_j的幅度进行2-范数归一，2-范数归一后的散射点矩阵记为A'_j；

(1d)根据预先给定的中心个数K，K≤Q，用K中心一类分类器，将训练样本对应的散射点集{A'₁,A'₂,...,A'_j,...,A'_Q}聚为K簇，同时保存这K簇的中心 $\{A_{c_1},A_{c_2},...,A_{c_t},...,A_{c_K}\},t=1,...,K;$

(1e)根据预先给定的拒判率P和K簇的中心计算K中心一类分类器的门限Thr，其中0≤P＜1；

(2)测试步骤：

(2a)估计测试切片z的信号能量比，并依此为门限提取散射点矩阵B，对散射点矩阵B的幅度进行2-范数归一，2-范数归一后的散射点矩阵记为B'；

(2b)计算测试样本的散射点矩阵B'与步骤(1e)得到的K簇的中心之间的最小双向Hausdorff距离d_B；

(2c)将测试样本的散射点矩阵B'与K簇的中心之间的最小双向Hausdorff距离d_B与步骤(1e)得到的K中心一类分类器的判决门限Thr进行比较，如果d_B≤Thr，则测试样本z为目标，否则，测试样本z为杂波虚警。

2.根据权利要求1所述的基于散射点和K中心一类分类器的SAR目标鉴别方法，其中步骤(1a)所述的对SAR图像I进行恒虚警检测，按如下步骤进行：

(1a1)假设目标的最大长度为l，雷达的距离分辨率和方位分辨率分别为ρ_r、ρ_a，则目标在SAR图像中的最大长度l_r和最大宽度l_a分别为：和选取警戒窗长度 $L>2\sqrt{l_r^2+l_a^2};$

(1a2)给定杂波区宽度N_r，虚警概率P_fa；

(1a3)选择双参数恒虚警检测器作为恒虚警检测器，根据高斯杂波统计模型和虚警概率P_fa计算出检测阈值T₁＝Φ^-1(1-P_fa)，其中Φ^-1(·)表示标准正态分布的反函数；

(1a4)对SAR图像I中像素的幅度进行如下变化：

b＝10*log₁₀(a²)

其中a表示图像I中某一像素的幅度，b表示图像I中某一像素的幅度变化后的幅度，变化后的图像用I₁表示；

(1a5)选取I₁(i₁,j₁)单元作为待检测单元，以待检测单元I₁(i₁,j₁)为中心，按照图3所示选取警戒区和杂波区Ω_c，用N_c表示杂波区域像素的个数，其中(i₁,j₁)表示待检测单元的坐标；

(1a6)估计杂波区的均值 ${\widehat{\mathrm{\&mu;}}}_c=\frac{1}{N_c}\underset{w=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{c}_w,$ 方差 ${\widehat{\mathrm{\&sigma;}}}_c=\frac{1}{N_c-1}\underset{w=1}{\overset{N_c}{\mathrm{\&Sigma;}}}{(c_w-{\widehat{\mathrm{\&mu;}}}_c)}^2,$ 其中c_w是杂波区域第w个单元的的幅度；

(1a7)用矩阵C标记图像I中对应位置是否是目标像素，设待检测单元I₁(i₁,j₁)处的幅度为t，若则I₁(i₁,j₁)是目标像素，矩阵C中对应位置标记为1，否则，I₁(i₁,j₁)是非目标像素，矩阵C中对应位置标记为0；

(1a8)对变化后图像I₁中的所有像素重复(1a5)—(1a7)，得到二值图像C。

3.根据权利要求1所述的基于散射点和K中心一类分类器的SAR目标鉴别方法，其中步骤(1a)所述的对二值图像C进行聚类，是对二值图像C进行区域合并，若C中任意两个区域重心的距离小于目标的最大长度l的一半，则合并区域，否则保留区域，区域合并后的二值图像记为疑似目标区域C₁。

4.根据权利要求1所述的基于散射点和K中心一类分类器的SAR目标鉴别方法，其中所述步骤(1c)中从训练样本y_j中提取散射点矩阵A_j，按如下步骤进行：

(1c1)对切片y_j进行恒虚警检测，得到切片y_j的二值图像y_jb，y_jb中被置为1的为目标像素y_jt，置为0的为非目标像素y_jc；

(1c2)计算切片y_j的信号能量比，其中⊙表示矩阵的Hadamard积，vec(·)表示将矩阵向量化；

(1c3)对切片y_j做2维快速逆傅立叶变换并向量化得到原始信号ξ＝vec(ifft2(y_j))，初始化信号残差r＝ξ，其中ifft2(·)表示对信号进行2维快速逆傅立叶变换；

(1c4)根据切片y_j的大小m×n，以及雷达的距离分辨率和方位分辨率ρ_r、ρ_a，假设距离维对应的单元个数为m，方位维对应的单元个数为n，则散射点的位置坐标范围为X＝[0:ρ_a:nρ_a]，Y＝[0:ρ_r:mρ_r]；

(1c5)根据散射点模型 $E(f,\mathrm{\&psi;})=\&PartialD;\&CenterDot;\exp \left(\frac{-\mathrm{\&upsi;}4\mathrm{\&pi;f}}{c}(x\cos \mathrm{\&psi;}+y\sin \mathrm{\&psi;})\right),$ 构造字典原子d_i'j'：

$d_{i^{\&prime;}{j}^{\&prime;}}=\frac{d_{i^{\&prime;}{j}^{\&prime;}}^{\&prime;}}{{\left|\right|d_{i^{\&prime;}{j}^{\&prime;}}^{\&prime;}\left|\right|}_2},$

$d_{i^{\&prime;}{j}^{\&prime;}}^{\&prime;}=\mathrm{vec}\left(\exp \left(\frac{-\mathrm{\&upsi;}4\mathrm{\&pi;f}}{c}(X\left(i^{\&prime;}\right)\cos \mathrm{\&psi;}+Y\left(j^{\&prime;}\right)\sin \mathrm{\&psi;})\right)\right),i^{\&prime;}=1,...,N_x,j^{\&prime;}=1,...,N_y,$

其中f、ψ分别是雷达的频率和方位角，exp(·)表示自然指数函数，||·||₂表示向量的2-范数，(·)^T表示转置，N_x是是散射点的位置坐标X的维数，N_y是散射点的位置坐标Y的维数，是某散射点的幅度，x、y分别是某散射点的距离维坐标和方位维坐标，是虚数单位，π是圆周率，c是光速，X(i')是X坐标的第i'个元素，Y(j')是Y坐标的第j'个元素；

(1c6)由(1c5)中的d_i'j'，构造字典 $D=[d_{\mathrm{1,1}},d_{\mathrm{1,2}},...,d_{1,N_y},...,d_{i^{\&prime;}{j}^{\&prime;}},...,d_{N_x,1},d_{N_x,2},...,d_{N_x,N_y}];$

(1c7)用空集θ作为参数集合，空集S作为原子集合；

(1c8)将字典D和信号残差r做内积，κ＝D·r，内积κ中的最大元素的序号为M，将D(:,M)加入到原子集合S，该列对应的参数(1 X(i') Y(j'))^T加入到参数集合θ，i'、j'满足i'×(N_y-1)+j'＝M，其中D(:,M)表示字典D中的第M列，X(i')是X坐标的第i'个元素，Y(j')是Y坐标的第j'个元素；

(1c9)更新散射点的幅度更新重构信号更新信号残差更新重构能量比用替换参数集合θ的第一行，更新参数散射点的幅度，其中表示穆尔-彭罗斯(Moore-Penrose)逆矩阵，(·)^H表示共轭转置；

(1c10)若重构能量比δ小于信号能量比R，则执行(1c8)—(1c10)，否则执行(1c11)；

(1c11)将参数集合θ赋值到散射点矩阵A_j，则 $A_j=\left[\begin{array}{ccccc}{a}_{j_1}& ...& a_{j_{i^{\&prime;\&prime;}}}& ...& a_{j_{N_j}}\\ {x_j}_1& ...& x_{j_{i^{\&prime;\&prime;}}}& ...& x_{j_{N_j}}\\ {y_j}_1& ...& y_{j_{i^{\&prime;\&prime;}}}& ...& y_{j_{N_j}}\end{array}\right],$ 其中分别表示散射点的幅度、距离维坐标、方位维坐标，i″＝1,...,N_j，N_j表示从y_j中提取的散射点的个数；

(1c12)对j＝1,...,Q重复执行(1c1)—(1c11)。

5.根据权利要求1所述的基于散射点和K中心一类分类器的SAR目标鉴别方法，其中所述步骤(1c)对散射点矩阵A_j的幅度进行2-范数归一，2-范数归一后的散射点矩阵记为A'_j，按如下步骤进行：

(1c13)对散射点矩阵A_j的幅度a_ji″，i″＝1,...,N_j，按如下公式进行2-范数归一：

${\stackrel{\&OverBar;}{a}}_{j_{i^{\&prime;\&prime;}}}=\frac{a_{j_{i^{\&prime;\&prime;}}}}{{\left|\right|{\left(\begin{array}{ccccc}{a}_{j_1}& ...& a_{j_{i^{\&prime;\&prime;}}}& ...& \end{array}{a}_{j_{N_j}}\right)}^T\left|\right|}_2},$

(1c14)将2-范数归一化后的散射点幅度代入散射点矩阵A_j的幅度得到归一化后的散射点矩阵A'_j为：

$A_j^{\&prime;}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{a}}_{j_1}& ...& {\stackrel{\&OverBar;}{a}}_{j_{i^{\&prime;\&prime;}}}& ...& {\stackrel{\&OverBar;}{a}}_{j_{N_j}}\\ x_{j_1}& ...& x_{j_{i^{\&prime;\&prime;}}}& ...& x_{j_{N_j}}\\ y_{j_1}& ...& y_{j_{i^{\&prime;\&prime;}}}& ...& x_{j_{N_j}}\end{array}\right],j=1,...,Q;$

6.根据权利要求1所述的基于散射点和K中心一类分类器的SAR目标鉴别方法，其中步骤(1d)所述的根据预先给定的中心个数K，K≤Q，用K中心一类分类器，将训练样本对应的散射点集{A'₁,A_'2,...,A'_j,...,A'_Q}聚为K簇，按如下步骤进行：

(1d1)将2-范数归一化后的的散射点矩阵 $A_j^{\&prime;}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{a}}_{j_1}& ...& {\stackrel{\&OverBar;}{a}}_{j_{i^{\&prime;\&prime;}}}& ...& {\stackrel{\&OverBar;}{a}}_{j_{N_j}}\\ x_{j_1}& ...& x_{j_{i^{\&prime;\&prime;}}}& ...& x_{j_{N_j}}\\ y_{j_1}& ...& y_{j_{i^{\&prime;\&prime;}}}& ...& x_{j_{N_j}}\end{array}\right]$ 分块表示为：

$A_j^{\&prime;}\text{=}\left[\begin{array}{ccccc}{b}_{j_1}& ...& b_{j_p}& ...& b_{j_{N_j}}\end{array}\right],$

其中 $b_{j_p}={\left(\begin{array}{ccc}{\stackrel{\&OverBar;}{a}}_{j_p}& x_{j_p}& y_{j_p}\end{array}\right)}^T,p=1,...,N_j,j=1,...,Q;$

(1d2)计算训练样本散射点集的协方差矩阵：其中Σ_j是切片y_j的散射点矩阵A'_j的协方差矩阵，其表示如下：

${\mathrm{\&Sigma;}}_j=\left[\begin{array}{ccc}{\mathrm{\&sigma;}}_{j1}^2& 0& 0\\ 0& {\mathrm{\&sigma;}}_{j2}^2& 0\\ 0& 0& {\mathrm{\&sigma;}}_{j3}^2\end{array}\right],$

其中，分别是训练样本散射点矩阵 $A_j^{\&prime;}=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{a}}_{j_1}& ...& {\stackrel{\&OverBar;}{a}}_{j_{i^{\&prime;\&prime;}}}& ...& {\stackrel{\&OverBar;}{a}}_{j_{N_j}}\\ x_{j_1}& ...& x_{j_{i^{\&prime;\&prime;}}}& ...& x_{j_{N_j}}\\ y_{j_1}& ...& y_{j_{i^{\&prime;\&prime;}}}& ...& x_{j_{N_j}}\end{array}\right]$ 的第1行元素，第2行元素和第3行元素的方差；

(1d3)随机从散射点集{A'₁,A'₂,...,A'_j,...,A'_Q}中选取K个样本作为簇中心计算代价函数 $J0=\underset{j=1,..,Q}{\max }\left(\underset{t=1,...,K}{\min }{\left|\right|A_j-A_{c_t}\left|\right|}_2^2\right);$

(1d4)计算训练样本的簇标号 $L\left(j\right)=\underset{t=1,...,K}{\arg \min }\mathrm{HD}(A_j^{\&prime;},A_{c_t}),j=1,...,Q,$ 其中：表示训练样本A'_j和簇中心的双向Hausdorff距离，

$\mathrm{HD}(A_j^{\&prime;},A_{c_t})=\max (h(A_j^{\&prime;},A_{c_t}),h(A_{c_t},A_j^{\&prime;}))$

$h(A_j^{\&prime;},A_{c_t})=\frac{1}{N_j}\underset{p=1}{\overset{N_j}{\mathrm{\&Sigma;}}}\underset{q=1,...,N_{c_t}}{\min }\left(d(b_{j_p},b_{c_{\mathrm{tq}}})\right)$

$h(A_{c_t},A_j^{\&prime;})=\frac{1}{N_{c_t}}\underset{q=1}{\overset{N_{c_t}}{\mathrm{\&Sigma;}}}\underset{p=1,...,N_j}{\min }\left(d(b_{c_{\mathrm{tq}}},b_{j_p})\right)$

$(b_{j_p},b_{c_{\mathrm{tq}}})=d(b_{c_{\mathrm{tq}}},b_{j_p})=\sqrt{{(b_{j_p}-b_{c_{\mathrm{tq}}})}^T{\mathrm{\&Sigma;}}^{-1}(b_{j_p}-b_{c_{\mathrm{tq}}})},$ (·)^-1表示矩阵的逆，表示第c_t个簇中心的散射点个数；

(1d5)任意取一非中心样本备份并用代替第L(m₁)簇的中心计算代价函数若J1＜J0，则接受这种代替，将J1赋值给J0，若J1≥J0，则不接受这种代替，用备份的恢复

(1d6)重复(1d4)—(1d5)直至簇中心不在发生变化。

7.根据权利要求1所述的基于散射点和K中心一类分类器的SAR目标鉴别方法，其中步骤(1e)所述根据预先给定的拒判率P和K簇的中心计算K中心一类分类器的门限Thr，按如下步骤进行：

(1e1)计算 $d_A\left(j\right)=\underset{t=1,...,K}{\min }\mathrm{HD}(A_j^{\&prime;},A_{c_t}),j=1,...,Q,$ d_A(j)是第j个样本和簇中心间的距离；

(1e2)根据预先给定的拒判率P和样本和簇中心间的距离，计算K中心一类分类器的门限 $\mathrm{Thr}=\left\{\begin{array}{cc}(d_A\left(\mathrm{num}\right)+d_A(\mathrm{num}+1))/2,& 1\&le;\mathrm{num}<Q\\ d_A\left(\mathrm{num}\right),& \mathrm{num}=Q\end{array}\right.,$ 表示向下取整，其中0≤P＜1。