✅作者简介：热爱科研的Matlab仿真开发者，擅长数据处理、建模仿真、程序设计、完整代码获取、论文复现及科研仿真。
🍎 往期回顾关注个人主页：Matlab科研工作室
👇 关注我领取海量matlab电子书和数学建模资料
🍊个人信条：格物致知,完整Matlab代码获取及仿真咨询内容私信
🔥 内容介绍
一、引言：山地环境下无人机路径规划的技术痛点与 APF 算法价值
无人机三维路径规划是实现自主飞行、任务执行的核心技术，在山地搜救、资源勘探、电力巡检等场景中具有不可替代的作用。复杂山地环境以高程突变、地形崎岖、障碍物密集（如山峰、峡谷、树林）、信号遮挡为典型特征，对路径规划提出三大核心要求：安全性（避开障碍与地形碰撞）、经济性（路径长度与能耗最优）、平滑性（满足无人机机动约束）。
传统路径规划算法（如 A*、RRT）在二维平面或简单三维场景中表现优异，但面对山地地形的连续高程约束与动态障碍感知需求时，存在搜索效率低、路径冗余度高、鲁棒性不足等问题。人工势场算法（Artificial Potential Field, APF）凭借模型简洁、实时性强、局部优化能力突出的优势，成为复杂环境下三维路径规划的优选方案 —— 其通过模拟 “引力 - 斥力” 的势场作用，使无人机在目标点引力牵引下向终点移动，同时在障碍物斥力作用下避开危险区域，无需复杂的全局搜索即可快速生成可行路径。
本文将聚焦复杂山地模型的建模方法，深入解析 APF 算法的三维扩展与改进策略，通过仿真验证其在山地环境中的适用性与优越性。
二、复杂山地模型的构建与关键约束分析
（一）山地地形建模技术
复杂山地模型的精准构建是路径规划的基础，需兼顾地形真实性与计算效率，常用方法包括：
数字高程模型（DEM）离散建模：基于实测高程数据（如 GPS、激光雷达点云），通过栅格化处理生成 DEM 模型，每个栅格单元存储对应位置的海拔高度。设栅格坐标系为（X, Y, Z），其中 Z 为高程值，通过双线性插值算法可获取任意坐标（x, y）的连续高程，满足路径规划的精细化需求。
地形特征提取与简化：从 DEM 模型中提取山峰（局部高程极大值）、峡谷（局部高程极小值）、陡坡（高程梯度＞45°）等关键特征，采用四叉树分块策略对平缓区域进行栅格简化，对复杂区域加密栅格，平衡建模精度与计算速度。
障碍物融合建模：将山地中的树林、建筑物、输电塔等静态障碍物，以 “圆柱模型” 或 “不规则多面体模型” 嵌入 DEM 中，明确障碍物的空间坐标范围与危险半径（如树林的碰撞风险半径设为 5m）。
（二）无人机路径规划的核心约束
地形约束：无人机飞行高度需高于地形高程一定安全裕度（通常 5-10m），避免与地面或凸起地形碰撞；同时需限制最大爬升 / 俯冲角（一般≤30°），匹配无人机动力性能。
动力学约束：路径曲率需满足无人机最小转弯半径（如多旋翼无人机最小转弯半径≥2m），速度与加速度需在无人机额定范围内（如最大飞行速度≤20m/s，最大加速度≤5m/s²）。
任务约束：根据实际场景需求，可能需加入航点遍历（如勘探点、巡检点）、避禁飞区等约束条件。
三、APF 算法的三维扩展与山地环境适配改进
（一）传统 APF 算法的核心原理
APF 算法的核心思想是在规划空间中构建虚拟势场：
引力势场：目标点对无人机产生引力，引力大小与无人机到目标点的距离成正比，引导无人机向终点移动；
斥力势场：障碍物对无人机产生斥力，斥力大小与无人机到障碍物的距离成反比，距离越近斥力越强，避免碰撞。
二维平面中，引力 F_att 与斥力 F_rep 的计算公式如下：
Image
⛳️ 运行结果
Image
📣 部分代码
function PlotSolution(BestSol ,model)
%% 绘图函数
subplot(2,2,1)
mesh( model.x_data , model.y_data , model.z_data ); hold on
colorbar; box on ,
set(gcf,'Color',[1 1 1]);
% set(gcf, 'unit' , 'centimeters','position' , [2 2 30 15]);
h3= plot3( BestSol.sol.xx , BestSol.sol.yy , BestSol.sol.zz , '-r'); hold on
% temp = 10^-2 ;
h1 = plot3( model.xs, model.ys, model.zs, 'o' , 'MarkerEdgeColor','r', ...
'MarkerFaceColor','r'); hold on
h2 = plot3( model.xt , model.yt , model.zt, '^' , 'MarkerEdgeColor','r', ...
'MarkerFaceColor','r'); hold on
if isfield( model ,'Barrier')
for ind = 1: model.Num_Barrier
[X,Y,Z] = cylinder( model.Barrier(ind, 3) ,100);
h4= surf(X+model.Barrier(ind, 1),Y+model.Barrier(ind, 2), model.zmin+Z*( max(model.z_data(:))- model.zmin) ) ; hold on
set(h4, 'edgecolor','m','facecolor', 'm') ;
end
end
if ~isfield( model ,'Barrier')
legend( [ h1 , h2 , h3] , '起点' , '终点' , '线路', 'Location','southoutside' , 'Orientation','horizontal')
else
legend( [ h1 , h2 , h3 h4] , '起点' , '终点' , '线路' , '无法通行区域' , 'Location','southoutside', 'Orientation','horizontal')
end
xlabel('x / km','fontsize',10 ,'fontname','Times new roman');
ylabel('y / km','fontsize',10 ,'fontname','Times new roman');
zlabel('z / km','fontsize',10 ,'fontname','Times new roman');
set(gca, 'xlim' , [ model.xmin model.xmax]) ;
set(gca, 'ylim' , [ model.ymin model.ymax]) ;
axis tight
axis normal
hold off
set(gca, 'cameraposition', [-337 , 0 , 14.2])
%%
subplot(2,2,2)
mesh( model.x_data , model.y_data , model.z_data ); hold on
view(2)
% contour( model.x_data , model.y_data , model.z_data ); hold on
colorbar; box on ,
set(gcf,'Color',[1 1 1]);
% set(gcf, 'unit' , 'centimeters','position' , [2 2 30 15]);
h3= plot3( BestSol.sol.xx , BestSol.sol.yy , BestSol.sol.zz , '-r'); hold on
% temp = 10^-2 ;
h1 = plot3( model.xs, model.ys, model.zs, 'o' , 'MarkerEdgeColor','b', ...
'MarkerFaceColor','b'); hold on
h2 = plot3( model.xt , model.yt , model.zt, '^' , 'MarkerEdgeColor','g', ...
'MarkerFaceColor','g'); hold on
% legend( [ h1 , h2 , h3] , '起点' , '终点' , '线路', 'Location','southoutside' , 'Orientation','horizontal')
if isfield( model ,'Barrier')
for ind = 1: model.Num_Barrier
[X,Y,Z] = cylinder( model.Barrier(ind, 3) ,100);
h4= surf(X+model.Barrier(ind, 1),Y+model.Barrier(ind, 2), model.zmin+Z*( max(model.z_data(:))- model.zmin) ) ; hold on
set(h4, 'edgecolor','m','facecolor', 'm') ;
end
end
if ~isfield( model ,'Barrier')
legend( [ h1 , h2 , h3] , '起点' , '终点' , '线路', 'Location','southoutside' , 'Orientation','horizontal')
else
legend( [ h1 , h2 , h3 h4] , '起点' , '终点' , '线路' , '无法通行区域' , 'Location','southoutside', 'Orientation','horizontal')
end
xlabel('x / km','fontsize',10 ,'fontname','Times new roman');
ylabel('y / km','fontsize',10 ,'fontname','Times new roman');
zlabel('z / km','fontsize',10 ,'fontname','Times new roman');
set(gca, 'xlim' , [ model.xmin model.xmax]) ;
set(gca, 'ylim' , [ model.ymin model.ymax]) ;
axis tight
axis normal
hold off
%%
subplot(2,2,3)
% mesh( model.x_data , model.y_data , model.z_data ); hold on
% view(2)
contour( model.x_data , model.y_data , model.z_data ); hold on
colorbar; box on ,
set(gcf,'Color',[1 1 1]);
% set(gcf, 'unit' , 'centimeters','position' , [2 2 30 15]);
h3= plot3( BestSol.sol.xx , BestSol.sol.yy , BestSol.sol.zz , '-r'); hold on
% temp = 10^-2 ;
h1 = plot3( model.xs, model.ys, model.zs, 'o' , 'MarkerEdgeColor','b', ...
'MarkerFaceColor','b'); hold on
h2 = plot3( model.xt , model.yt , model.zt, '^' , 'MarkerEdgeColor','g', ...
'MarkerFaceColor','g'); hold on
% legend( [ h1 , h2 , h3] , '起点' , '终点' , '线路', 'Location','southoutside' , 'Orientation','horizontal')
if isfield( model ,'Barrier')
for ind = 1: model.Num_Barrier
[X,Y,Z] = cylinder( model.Barrier(ind, 3) ,100);
h4= surf(X+model.Barrier(ind, 1),Y+model.Barrier(ind, 2), model.zmin+Z*( max(model.z_data(:))- model.zmin) ) ; hold on
set(h4, 'edgecolor','m','facecolor', 'm') ;
end
end
if ~isfield( model ,'Barrier')
legend( [ h1 , h2 , h3] , '起点' , '终点' , '线路', 'Location','southoutside' , 'Orientation','horizontal')
else
legend( [ h1 , h2 , h3 h4] , '起点' , '终点' , '线路' , '无法通行区域' , 'Location','southoutside', 'Orientation','horizontal')
end
xlabel('x / km','fontsize',10 ,'fontname','Times new roman');
ylabel('y / km','fontsize',10 ,'fontname','Times new roman');
zlabel('z / km','fontsize',10 ,'fontname','Times new roman');
set(gca, 'xlim' , [ model.xmin model.xmax]) ;
set(gca, 'ylim' , [ model.ymin model.ymax]) ;
axis tight
axis normal
hold off
%% 收敛曲线
subplot(2,2,4)
% semilogy( BestSol.BestCost ,'LineWidth',2);
plot( BestSol.BestCost ,'LineWidth',2);
xlabel('迭代次数');
ylabel('目标函数');
grid on;
set(gca,'XLim',[0 BestSol.MaxIt]);%X轴的数据显示范围
🔗 参考文献
🎈 部分理论引用网络文献，若有侵权联系博主删除
🏆团队擅长辅导定制多种科研领域MATLAB仿真，助力科研梦：
