pointcloud_generation.m 点云生成将cfar检测到的点从极坐标系转直角坐标系。functionpointCloudpointcloud_generation(detections)% detections: CFAR检测后的[range, doppler, angle, rcs]detectPointsNumsize(detections,1);pointCloudzeros(detectPointsNum,5);%[x y zvpower]fori1:detectPointsNum rangedetections(i,1);% 距离估计 velocitydetections(i,2);% 速度估计 azdetections(i,3);% 角度估计 % 转笛卡尔坐标 xrange * sind(az);yrange * cosd(az);z0;% 当前假设2D线性阵列垂直高度设为0【可扩展为俯仰】 pointCloud(i,:)[x, y, z, velocity, detections(i,4)];%[x y zvpower]end % % 绘制点云图 % X_pcpointCloud(:,1);% 纵向距离(m)% Y_pcpointCloud(:,2);% 横向距离(m)% Z_pcpointCloud(:,3);% 高度(m)% V_pcpointCloud(:,4);% 径向速度(m/s)% P_pcpointCloud(:,5);% 反射功率(dB或线性)% %2. 功率映射为点的大小(防止功率数值太大或太小导致点消失或撑满屏幕)% % 使用归一化动态映射将点的大小限制在20到150之间 %ifmax(P_pc)~min(P_pc)% marker_size20130*(P_pc - min(P_pc))/(max(P_pc)- min(P_pc));%else% marker_size40* ones(size(P_pc));% 若功率全等固定大小 % end % % 绘制点云图 % figure(Color,[111]);% % 使用 scatter3 绘制4D点云 % % X, Y, Z 为空间坐标marker_size 决定大小V_pc 决定颜色 % scatter3(X_pc,Y_pc, Z_pc, marker_size, V_pc,filled,MarkerEdgeColor,[0.30.30.3]);% %4. 绘制雷达原点方便对比相对位置 % hold on;% plot3(0,0,0,rp,MarkerSize,15,MarkerFaceColor,r);% text(0,0,0.5, 雷达(0,0,0),Color,r,FontWeight,bold);% %5. 润色图表 % grid on;box on;% colormap(jet);% 蓝冷色代表速度小/靠近红暖色代表速度大/远离 % cbcolorbar;% ylabel(cb,径向速度 v (m/s),FontSize,11);% % 坐标轴标签(严格对照雷达直角坐标系)% xlabel(纵向距离 X (米),FontSize,11);% ylabel(横向距离 Y (米),FontSize,11);% zlabel(高度 Z (米),FontSize,11);% title(雷达 4D 空间点云成像结果 (大小:功率, 颜色:速度),FontSize,12,FontWeight,bold);% % 保持几何比例1:1防止飞机被压扁或者拉长 % axis equal;view(0,90);end中间结果对于仿真的单个飞机目标该目标质心的[距离, 目标中心与雷达角度, 航向速度, 航向角][60, 10, -4, 225]绘制的点云图如下和仿真飞机目标位置与速度基本对应上了。聚类% 点的聚类 % 输入1聚类算法1DBSCAN2K-Means3X-Means % 输入2数据集 data(N x5)% 输入3聚类半径 EPS % 输入4最小数据量 minPoints % 输出1聚类数据 data_clusterfunctiondata_clusterpoints_clustering(clusterMethod, data, eps, minpoints)% data: N x5矩阵[x, y, z, v, power]ifnargin2eps2.5;% 最大搜索半径 【根据飞机大小设置不然单个飞机会出现多个簇】 minpoints10;% 最小点数 end飞机的姿态航向角直接决定了最终点云在直角坐标系里的空间疏密和延展方向。当飞机朝着雷达迎头直飞时机头挡住了机身而后掠翼的边缘与雷达射线的夹角极其锐利此时雷达漏掉大批点迹点云在空间上会缩水成一个极其狭窄的“一字形”或几个孤立的强点如发动机、机头。可根据飞机姿态调整聚类参数。聚类半径侧飞/斜切进场点云饱满此时 eps 不需要太大迎头对冲/背道而驰点迹可能断层调大eps最小聚集点数侧飞/斜切进场抬高最小聚集点数迎头对冲/背道而驰减小聚集点数如何知道飞机是否侧飞呢1.根据cfar检测结果每个点的径向速度已知如果飞机是笔直迎头飞来机身所有点的径向速度几乎一致。飞机侧飞时不同部位的径向速度就会产生明显差异。所以可计算检测点的径向速度标准差作为评估指标越大就表明飞机斜切越厉害。不过这种计算标准差仅适用于单个目标。2.可以利用前面帧卡尔曼滤波器估计到的飞机运动轨迹和航速判断飞机是否侧飞做自适应聚类。%【除了用距离也可以加上多普勒速度进行聚类】 featuresdata(:,1:3);% 提取用于聚类的特征维度[x,y,z]方便做 norm 距离计算ifclusterMethod1% DBSCAN %1-----调用MATLAB内置的高效函数dbscan % idx 返回每个点对应的簇标签。如果点被判定为噪声/离群点idx 会被赋予-1% idxdbscan(features, eps, minpoints);%[dataNum,1]%2-----调用实现的my_dbscan函数 idxmy_dbscan(features, eps, minpoints);% 剔除噪声点离群点 outliers_mask(idx-1);data(outliers_mask,:)[];idx(outliers_mask)[];data_cluster[data, idx];%[dataNum,6]最后一列 簇标签idx聚类ID对于仿真的单个飞机目标该目标质心的[距离, 目标中心与雷达角度, 航向速度, 航向角][60, 10, -4, 225]matlab自带dbscan中间结果my_dbscan中间结果elseif clusterMethod2% K-Means % 局限性说明K-Means 无法识别噪声点且必须输入固定的 K 值 % 在实际多目标雷达管线中通常不建议直接用标准 K-Means % 假设场上有3个目标 fixed_K3;ifsize(features,1)fixed_K idxones(size(features,1),1);% 点数不够分塞进一类else% 调用matlab内置的 kmeansReplicates,3表示尝试3次初始化取最优防止陷入局部极小值 idxkmeans(features, fixed_K,Replicates,3,Display,off);end data_cluster[data, idx];不建议直接用kmeans。elseif clusterMethod3% X-Means % 自适应K值从 K_min 开始通过 BIC 评分自动分裂直到 K_max k_min1;k_max8;% 限制单帧最大可能的目标数 idxrun_xmeans(features, k_min, k_max, targetnum);data_cluster[data, idx];end plot_cluster_birdseye(data_cluster, clusterMethod);end子函数部分自实现的dbscan算法用于了解算法过程。functionidxmy_dbscan(X, epsilon, minPts)% 输入: % X - N x2或 N x3的点云/点迹坐标矩阵 % epsilon - 邻域半径距离门限 % minPts - 成为核心点所需的最小邻居数含自身 % 输出: % idx - N x1向量记录每个点所属的聚类ID0代表噪声/野点 numPointssize(X,1);idxzeros(numPoints,1);% 初始化0 代表未分类Unvisited / Noise clusterId0;% 聚类标签计数器 % 步骤1利用 MATLAB 矩阵矩阵秒算“两两欧氏距离矩阵”(N x N)% pdist2Pairwise Distance Matrix成对距离矩阵计算 比双重for循环快几百倍 % Dpdist2(X, Y)矩阵D的行数等于X的点数列数等于Y的点数 distMatrixpdist2(X, X);%[N, N]对称矩阵对角线全0自己到自己的距离是0 % 步骤2主循环遍历每一个点fori1:numPoints % 如果这个点已经被某个群吞并了直接跳过ifidx(i)~0continue;end % 寻找当前点的所有邻居索引 neighborsfind(distMatrix(i,:)epsilon);% 步骤3密度检查iflength(neighbors)minPts % 邻居太少暂时标记为噪声0等别人来吞并它 idx(i)0;else% 邻居足够恭喜你成为核心点开辟一个新帮派 clusterIdclusterId 1;idx(i)clusterId;% 开始“顺藤摸瓜”把邻居放入一个动态队列中去扩展 queueneighbors;queue(queuei)[];% 队列中移出自己while~isempty(queue)currentPointqueue(1);% 取出队列第一个点 queue(1)[];% 出队 % 如果这个点之前被误判为噪声现在把它收编进当前帮派ifidx(currentPoint)0idx(currentPoint)clusterId;end % 如果这个点已经被其他正常帮派收编了跳过它 % 边界点的归属采用先到先得的原则允许微小的边界模糊比强行去合并两个独立的物理目标安全得多ifidx(currentPoint)~clusterIdidx(currentPoint)0continue;end % 关键扩展检查这个邻居是不是也是一个核心点 currentNeighborsfind(distMatrix(currentPoint,:)epsilon);% 如果它也是核心点有带小弟的能力iflength(currentNeighbors)minPts % 它的所有未分类小弟都要加入我们的考察队列fork1:length(currentNeighbors)neighborPointcurrentNeighbors(k);ifidx(neighborPoint)0% 未访问过或曾是噪声 idx(neighborPoint)clusterId;% 现场收编 queue[queue, neighborPoint];%#okAGROW % 加入扩展队列end end end end %while队列扩展结束 end end %for主循环结束 endxmeans聚类算法functionidxrun_xmeans(X, k_min, k_max)[N, M]size(X);% N个点M个维度(3)ifNk_min idxones(N,1);% 塞进一类return;end % 初始状态从 k_min 开始做初始 K-means current_Kk_min;idxkmeans(X, current_K,Replicates,2,Display,off);% 开始迭代尝试分裂whilecurrent_Kk_max changedfalse;new_idxidx;offset0;% 用于修正分裂后的标签偏移fori1:current_K % 抠出当前第 i 个簇的所有点 cluster_mask(idxi);X_subX(cluster_mask,:);N_subsize(X_sub,1);% 如果这个簇点数太少少于4个点就没必要再分裂了ifN_sub4new_idx(cluster_mask)i offset;continue;end % 【物理几何尺寸叫停机制】 % 如果一个簇在 X距离、Y横向距离、Z高度方向上的最大跨度已经小于你雷达的物理分辨率 % 或者一个合理的物理阈值比如0.3米那它在物理上绝对已经是一个不可分割的“单体”了。 % 计算这个簇在 X, Y, Z 三个物理维度上的最大跨度几何尺寸 cluster_spanmax(X_sub,[],1)- min(X_sub,[],1);% 如果它在任何一个方向的跨度都小于0.6米说明已经聚成一团完美点迹了 % 直接强行跳过不准它参与任何 BIC 计算和分列ifall(cluster_span0.6)% 这个值经过几次调试得到一开始设置的0.4还是会将单个点目标分割 new_idx(cluster_mask)i offset;continue;% 强行保命不许切 end % 计算不分裂时的原始 BIC 分数 bic_parentcalculate_bic(X_sub,1);% 尝试对这个簇进行2-Means 强行一分为二 try sub_idxkmeans(X_sub,2,Replicates,2,Display,off);% 计算分裂后的新 BIC 分数 bic_childcalculate_bic(X_sub, sub_idx);catch bic_child-Inf;% 如果 kmeans 报错放弃分裂 end % 判断标准如果分裂后的 BIC 分数更高说明切开更好ifbic_childbic_parent % 接受分裂将原先的标签 i变成两个互不干扰的新标签 sub_idx_adjustedsub_idx;sub_idx_adjusted(sub_idx1)i offset;sub_idx_adjusted(sub_idx2)current_K offset 1;new_idx(cluster_mask)sub_idx_adjusted;offsetoffset 1;% 增加了一个新簇 changedtrue;else% 拒绝分裂保持原样仅做标签修正 new_idx(cluster_mask)i offset;end end % 更新总体局势 idxnew_idx;current_Kcurrent_K offset;% 如果遍历了所有簇没有任何一个簇想继续分裂或者达到了最大 K 限制就退出if~changed||current_Kk_maxbreak;end end % 规整标签使其连续例如把1,3,7变成1,2,3[~, ~, idx]unique(idx);end %% 数学工具BIC(贝叶斯信息准则)计算器functionbiccalculate_bic(X, idx)% 该函数专门用于评估在 X 数据集上当前的聚类划分 idx 分数是多少[N, M]size(X);Klength(unique(idx));% 计算全场方差(估计最大似然)variance0;fori1:K cluster_mask(idxi);X_subX(cluster_mask,:);ifsize(X_sub,1)1centroidmean(X_sub,1);variancevariance sum(sum((X_sub-centroid).^2));end end variancevariance /(N - K);ifvarianceeps bic-Inf;% 方差为0说明重合属于无效划分return;end % 计算对数似然度(Log-Likelihood)log_likelihood0;fori1:K cluster_mask(idxi);N_nsum(cluster_mask);ifN_n0% X-Means 论文中的标准极大似然估计公式 log_likelihoodlog_likelihood N_n * log(N_n)- N_n * log(N)... -(N_n * M /2)* log(2* pi * variance)... -(N_n -1)/2;end end % 计算自由参数个数(自由度)% 每个中心有 M 维坐标加上一个方差项 num_parametersK *(M 1);% BIC 终极公式得分似然度 - 参数复杂度惩罚项 % 惩罚项能有效防止算法无节制地分裂成无数个极小的簇 penalty_factor2.5;biclog_likelihood - penalty_factor*(num_parameters /2)* log(N);endfunctionplot_cluster_birdseye(data_cluster, clusterMethod)% plot_cluster_birdseye - 绘制雷达点云聚类的二维鸟瞰图 % % 输入参数: % data_cluster:N x6矩阵[x, y, z, v, power, cluster_id]% 最后一列必须是聚类算法输出的有效簇标签 idx % clusterMethod: 标量1DBSCAN2K-Means3X-Means(仅用于标题字符串切换)% ---1. 安全防呆检查 ---ifisempty(data_cluster)warning(输入的数据集为空取消图表绘制);return;endifsize(data_cluster,2)6error(输入数据格式不合理矩阵至少需要6列 [x, y, z, v, power, cluster_id]。);end % ---2. 提取聚类标签并去重 --- idxdata_cluster(:,6);unique_clustersunique(idx);num_clusterslength(unique_clusters);% ---3. 核心绘图引擎 ---if~isempty(unique_clusters)% 创建独立画布设置合适的分辨率和宽高比 figure(Color,w,Position,[200,200,800,650]);hold on;box on;grid on;% 定义不同独立目标的图形样式库防止目标太多导致标记重复 marker_library{o,s,^,d,v,,,p,h};forii1:num_clusters cluster_idunique_clusters(ii);% 筛选出当前目标集群的所有散射点迹 cluster_mask(data_cluster(:,6)cluster_id);cluster_datadata_cluster(cluster_mask,:);x_valcluster_data(:,1);y_valcluster_data(:,2);v_valcluster_data(:,4);% 提取第四列的物理速度 % 根据当前的循环序数自适应选取图形形状 m_stylemarker_library{mod(ii-1, length(marker_library))1};% 使用 scatter 绘制二维投影通过 v_val 实现动态颜色映射 scatter(x_val, y_val,45, v_val, m_style,filled,...LineWidth,1,MarkerEdgeColor,k,...DisplayName,[目标集群 , num2str(cluster_id)]);end % ---4. 视角锁定与几何形状保真对账 --- view(0,90);% 强制锁定正上方绝对俯视视角 axis equal;% 强行保持一比一物理横纵比防止飞机等扩展目标形状被拉伸变形 % ---5. 色彩轴(Colorbar)渲染配置 --- ccolorbar;c.Label.String雷达径向速度 (m/s);c.Label.FontSize11;colormap(jet);% jet 渲染器蓝色靠近负速红色远离正速 % 动态界定速度色彩轴的显示范围防止全场单一速度时 colorbar 报错 v_minmin(data_cluster(:,4));v_maxmax(data_cluster(:,4));ifv_max ~v_min caxis([v_min -1, v_max 1]);end % ---6. 标签与图例张贴 --- method_names{DBSCAN,K-Means,X-Means};ifclusterMethod1clusterMethod3method_strmethod_names{clusterMethod};elsemethod_str未知算法;end title(sprintf(%s 聚类鸟瞰图 (颜色表示速度), method_str),FontSize,12);xlabel(X 轴 (横向距离 m),FontSize,11);ylabel(Y 轴 (纵向距离 m),FontSize,11);legend(Location,best);elsewarning(数据集中未检测到任何有效分类标签请检查前级聚类算法);end endx-means是k-means的高级版。自适应分类数通过比较每个簇的分裂前后评分确定最佳分类数。关键点我的仿真点总是被划分为好几类所以加了一个物理几何叫停机制就是当两个簇之间的跨度小于指定距离时可通过雷达分辨率确定或者自己尝试几次确定那它就是一个单体。分裂前后的BIC评分机制是贝叶斯和最大似然估计加上参数复杂度惩罚项。特点有个物理几何尺寸叫停机制。这里我尝试到est_siz 9.0 时才聚类出实际的1个目标。