雷达目标识别最常用的聚类算法之一——密度聚类DBSCAN
当雷达检测目标点的数据达到一定程度,基于原始检测点目标的跟踪算法耗时耗力,且容易导致目标分裂,所以常见的处理方式是在建立目标的关联航迹之前,先对点目标进行聚类预处理,比如200个目标检测点,经过DBSCAN聚类之后,产生50个左右的聚类目标,在50个聚类目标的基础上再去进行航迹的关联跟踪,这样做有几点好处:
1) 大大减少了目标跟踪的复杂度和时间;
2) 对大目标的分裂问题友好;
3) 通过聚类的结果对单一的虚假点进行过滤,减少虚假航迹;
4) 通过类内的点的特征与分布,可以帮助估计目标的特征与运动形态;
本文提供了DBSCAN的matlab代码实现,但在具体应用的过程中,需要根据实际雷达检测点的质量、精度、分辨率等信号的指标特性,动态调整DBSCAN的代码逻辑或者聚类门限,以下代码仅供大家学习和研究。
%----------------DBSCAN聚类算法-----------------
% 基于密度聚类,不需要设定类的数量
% 优势:首先,它根本不需要pe-set数量的簇。 它还将异常值识别为噪声,而不像mean-shift,即使数据点非常不同,
% 它们也会将它们引入群集中。 另外,它能够很好地找到任意大小和任意形状的簇。
% 缺点:当簇的密度不同时,DBSCAN的性能不如其他组织。 这是因为当密度变化时,用于识别邻近点的距离阈值ε和
% minPoints的设置将随着群集而变化。 对于非常高维的数据也会出现这种缺点,因为距离阈值ε再次难以估计。
%-----------------------------------------------------
function [IDX, isnoise]=DBSCAN(X,epsilon,MinPts)
C = 0;
n = size(X,1);
IDX = zeros(n,1);
D = pdist2(X,X); % 计算向量之间的距离
visited = false(n,1);
isnoise = false(n,1);
for i = 1:n
if ~visited(i)
visited(i) = true;
Neighbors = RegionQuery(i); % 找到第i个点epsilon范围内的点
if numel(Neighbors) < MinPts % numel()返回满足条件的元素个数
% X(i,:) is NOISE
isnoise(i) = true;
else
C=C+1;
ExpandCluster(i,Neighbors,C);
end
end
end
function ExpandCluster(i,Neighbors,C)
IDX(i) = C;
k = 1;
while true
j = Neighbors(k);
if ~visited(j)
visited(j) = true;
Neighbors2 = RegionQuery(j);
if numel(Neighbors2) >= MinPts
Neighbors=[Neighbors Neighbors2]; %#ok
end
end
if IDX(j) == 0
IDX(j) = C;
end
k = k + 1;
if k > numel(Neighbors)
break;
end
end
end
function Neighbors = RegionQuery(i)
Neighbors = find(D(i,:) <= epsilon);
end
end
clc;
clear;
close all;
%% test Data
% 这里使用的iris数据的一部分,由于第3维和第4为数据数据区分度好,因此用3、4维数据测试
% Iris-setosa
X1 = [5.1,3.5,1.4,0.2;
4.9,3.0,1.4,0.2;
4.7,3.2,1.3,0.2;
4.6,3.1,1.5,0.2;
5.1,3.7,1.5,0.4;
4.6,3.6,1.0,0.2;
5.1,3.3,1.7,0.5;
5.0,3.6,1.4,0.2;
5.4,3.9,1.7,0.4;
4.6,3.4,1.4,0.3;
5.0,3.4,1.5,0.2;
4.4,2.9,1.4,0.2;
4.9,3.1,1.5,0.1;
5.4,3.7,1.5,0.2;
4.8,3.4,1.6,0.2;
4.8,3.0,1.4,0.1;
4.3,3.0,1.1,0.1;
5.8,4.0,1.2,0.2;
5.7,4.4,1.5,0.4;
5.4,3.9,1.3,0.4;
5.1,3.5,1.4,0.3;
5.7,3.8,1.7,0.3;
5.1,3.8,1.5,0.3;
5.4,3.4,1.7,0.2;
6.4,3.2,4.5,1.5;
6.9,3.1,4.9,1.5;
5.5,2.3,4.0,1.3;
6.5,2.8,4.6,1.5;
5.7,2.8,4.5,1.3;
6.3,3.3,4.7,1.6;
4.9,2.4,3.3,1.0;
4.9,2.4,3.3,1.0;
6.6,2.9,4.6,1.3;
5.2,2.7,3.9,1.4;
5.0,2.0,3.5,1.0;
5.9,3.0,4.2,1.5;
6.0,2.2,4.0,1.0];
X = X1(:,3:4);
A = X;
%% 计算每个点的最邻近距离,并存储在Kdist中
% knnsearch K近邻算法是找到K个最近的邻居。% KNN k distance graph, to determine the epsilon
numData = size(A,1);
Kdist = zeros(numData,1);
[IDX,Dist] = knnsearch(A(2:numData,:),A(1,:));% IDX 返回A(2:numData,:)中离A(1,:)最近的向量索引值,Dist表示对应的距离
Kdist(1) = Dist;
for i = 2:numData
[IDX,Dist] = knnsearch(A([1:i-1,i+1:numData],:),A(i,:));
Kdist(i) = Dist;
end
[sortKdist,sortKdistIdx] = sort(Kdist,'descend');
distX = [1:numData]';
plot(distX,sortKdist,'r+-','LineWidth',2);
set(gcf,'position',[1000 340 350 350]);
grid on;
%% Run DBSCAN Clustering Algorithm
epsilon = 0.2;
MinPts = 3;
IDX1=DBSCAN(A,epsilon,MinPts);
%% Plot Results
figure;
PlotClusterinResult(A, IDX1);
title(['DBSCAN Clustering (\epsilon = ' num2str(epsilon) ', MinPts = ' num2str(MinPts) ')']);
set(gcf,'position',[30 -10 500 500]);
%% Run DBSCAN Clustering Algorithm
epsilon= 0.25 ;
MinPts= 3 ;
IDX2=DBSCAN(A,epsilon,MinPts);
%% Plot Results
figure;
PlotClusterinResult(A, IDX2);
title(['DBSCAN Clustering (\epsilon = ' num2str(epsilon) ', MinPts = ' num2str(MinPts) ')']);
set(gcf,'position',[530 -10 500 500]);
%% Run DBSCAN Clustering Algorithm
epsilon= 0.5 ;
MinPts= 3 ;
IDX3=DBSCAN(A,epsilon,MinPts);
%% Plot Results
figure;
PlotClusterinResult(A, IDX3);
title(['DBSCAN Clustering (\epsilon = ' num2str(epsilon) ', MinPts = ' num2str(MinPts) ')']);
set(gcf,'position',[30 380 500 500]);
%结果显示子函数,需另外创建.m文件保存
% Copyright (c) 2015, Yarpiz (www.yarpiz.com)
% All rights reserved. Please read the "license.txt" for license terms.
%
% Project Code: YPML110
% Project Title: Implementation of DBSCAN Clustering in MATLAB
% Publisher: Yarpiz (www.yarpiz.com)
%
% Developer: S. Mostapha Kalami Heris (Member of Yarpiz Team)
%
% Contact Info: sm.kalami@gmail.com, info@yarpiz.com
function PlotClusterinResult(X, IDX)
k=max(IDX);
Colors=hsv(k);
Legends = {};
for i=0:k
Xi=X(IDX==i,:);
if i~=0
Style = 'x';
MarkerSize = 8;
Color = Colors(i,:);
Legends{end+1} = ['Cluster #' num2str(i)];
else
Style = 'o';
MarkerSize = 6;
Color = [0 0 0];
if ~isempty(Xi)
Legends{end+1} = 'Noise';
end
end
if ~isempty(Xi)
plot(Xi(:,1),Xi(:,2),Style,'MarkerSize',MarkerSize,'Color',Color);
end
hold on;
end
hold off;
axis equal;
grid on;
legend(Legends);
legend('Location', 'NorthEastOutside');
end