简介本资源为2022年全国大学生数学建模竞赛本科组B题获奖论文聚焦无人机纯方位无源定位这一典型控制与几何建模问题面向数学建模参赛学生、自动化/无人机方向初学者及MATLAB实践者。全文基于圆形与锥形两种编队场景系统构建了含三类子问题的完整求解框架包括被动接收无人机的二维坐标定位模型融合余弦定理与三圆共点原理、最小发射机数量判定方法通过枚举与偏差平方和优化验证仅需增加1架、以及带步长遍历与迭代调整的位置校正方案支持0.01米精度逼近。资源为1个2.63MB的PDF文件内容涵盖问题重述、严密推导、MATLAB实现逻辑说明及关键词提炼结构完整、公式详实、思路可复现。目前已有8039人学习下载是理解无源定位几何本质、训练建模抽象能力与MATLAB工程化表达的优质原创范例。1. 纯方位无源定位不是玄学一份能跑通的MATLAB建模实录专治无人机编队“飘了”你有没有遇到过这种场景十架无人机排成标准圆形编队飞着飞着2号机突然偏了3米5号机悄悄往南挪了1.7米整个队形肉眼可见地“松垮”——但偏偏不能发主动信号校准因为一开雷达就暴露位置、干扰邻机、破坏电磁静默。这时候纯方位无源定位Bearing-Only Passive Localization就是唯一的解药不发射、只接收靠三架“锚点机”的方向夹角反推自己在哪。这不是理论空谈而是2022年国赛B题的真实战场。这篇笔记不讲定义、不画大饼直接拆开那份拿下河北省一等奖的MATLAB源码包——它包含4个完整可运行脚本定位模型、选锚点模型、圆编队调位、锥编队调位所有坐标系、余弦定理推导、步长遍历逻辑、三圆共点交点求解全在代码里埋好了注释。新手照着改两行参数就能出图老手一眼看出关键约束为什么必须≥3个方位角为什么钝角轨迹是“眼睛”、锐角是“葫芦”为什么锥形编队第二轮总卡在4/7/9号机这些坑当年某高校参赛队在凌晨三点debug时用血泪填平了现在原样端给你。2. 定位模型怎么立从二维坐标系到三圆共点MATLAB里每行代码都在解几何题纯方位定位的核心从来不是“测角度”而是“解交点”。当一架无人机收到三架锚点机的方向信息α₁, α₂, α₃它实际是在告诉自己“我在以锚点A、B为弦张角为α₁的圆弧上也在以B、C为弦张角为α₂的圆弧上还在以A、C为弦张角为α₃的圆弧上。”——三条圆弧若共点那个交点就是它的唯一解。这个“三圆共点”不是假设是几何必然而MATLAB要做的就是把纸面推导变成可计算的数值解。下面拆解附录一的核心逻辑。2.1 坐标系搭建与锚点预置半径100m的圆周不是随便设的建模第一步必须固化参考系。题目明确“九架机均匀分布在半径为100m的圆周上一架在圆心”这意味着所有锚点坐标是确定的——不是变量是常量。MATLAB脚本里用极坐标转直角坐标的硬编码实现比符号运算快一个数量级% 圆形编队锚点坐标预置单位米 r 100; % 编队半径 % FY00: 圆心 (0,0) % FY01-FY09: 圆周上9架机间隔40度360°/9 theta_list [0, 40, 80, 120, 160, 200, 240, 280, 320]; % 角度列表 x_anchor r * cosd(theta_list); % x坐标 y_anchor r * sind(theta_list); % y坐标 % FY00坐标单独补上 x_anchor [0, x_anchor]; % [FY00, FY01, ..., FY09] y_anchor [0, y_anchor];参数说明cosd/sind是MATLAB专用角度函数输入度数非弧度避免新手因cos(pi/3)写成cos(60)导致结果全错。theta_list严格按题目“均匀分布”生成40度间隔是硬约束改它等于改题设。2.2 余弦定理反解边长三个方程解三个未知距离接收机位置未知但它到三架锚点的距离k₀, a₂, a₃是隐含变量。已知的是锚点间距离由预置坐标算出、以及接收机观测到的三个夹角α₁, α₂, α₃。根据余弦定理可列三个方程在△FY00-FY0i-FY0k中(k₀² a₃² - r²) / (2*k₀*a₃) cos(α₁)在△FY00-FY0j-FY0k中(k₀² a₂² - r²) / (2*k₀*a₂) cos(α₂)在△FY0i-FY0j-FY0k中(a₂² a₃² - a₁²) / (2*a₂*a₃) cos(α₃)其中a₁是锚点i与j的距离已知r是FY00到FY0i/j的距离已知。MATLAB用solve符号求解% 已知量锚点i,j,k编号对应角度alpha1,alpha2,alpha3 FY0i 1; FY0j 2; FY0k 3; % 示例FY00,FY01,FY02为锚点 alpha1 45; alpha2 60; alpha3 75; % 接收机实测角度度 % 计算锚点间距离硬编码避免实时计算误差 a1 sqrt((x_anchor(FY0j1)-x_anchor(FY0i1))^2 (y_anchor(FY0j1)-y_anchor(FY0i1))^2); a3 sqrt((x_anchor(FY0k1)-x_anchor(FY0i1))^2 (y_anchor(FY0k1)-y_anchor(FY0i1))^2); a2 sqrt((x_anchor(FY0k1)-x_anchor(FY0j1))^2 (y_anchor(FY0k1)-y_anchor(FY0j1))^2); % 符号求解三个距离 syms k0 a2_sym a3_sym eq1 (k0^2 a3_sym^2 - r^2)/(2*k0*a3_sym) cosd(alpha1); eq2 (k0^2 a2_sym^2 - r^2)/(2*k0*a2_sym) cosd(alpha2); eq3 (a2_sym^2 a3_sym^2 - a1^2)/(2*a2_sym*a3_sym) cosd(alpha3); sol solve([eq1, eq2, eq3], [k0, a2_sym, a3_sym]); k0_val double(sol.k0); % 距离FY00的长度 a2_val double(sol.a2_sym); % 距离FY0j的长度 a3_val double(sol.a3_sym); % 距离FY0i的长度逻辑说明这里solve返回的是符号解double()强制转数值。注意索引FY0i1因为x_anchor(1)存的是FY00所以FY0i编号需1对齐数组。若此处索引错一位后续所有坐标全偏移——这是当年某队初稿翻车点。2.3 从距离到坐标极坐标转换是精度命门解出k₀到圆心距离后还需知道方位角θ才能得直角坐标(x,y)。题目给出关键线索“把第i架无人机的准确位置固定在x轴上”即FY0i在(100,0)其方位角为0°。接收机相对于FY0i的方位角由α₁和几何关系反推% 计算接收机相对于FY0i的方位角以FY0i为原点x轴正向为0° % 根据△FY00-FY0i-FY0k中FY00在(0,0), FY0i在(r,0), FY0k在(x_k,y_k) % 利用余弦定理先求∠FY00-FY0i-FY0k再结合α₁得方位角 theta_rel acosd((r^2 k0_val^2 - a3_val^2)/(2*r*k0_val)); % ∠FY00-FY0i-FY0k theta_abs theta_rel - 40*(FY0i-1); % 减去FY0i在圆周上的基准角FY010°, FY0240°... % 转直角坐标以全局原点FY00为基准 x_est k0_val * cosd(theta_abs); y_est k0_val * sind(theta_abs);参数说明theta_abs是全局坐标系下的方位角40*(FY0i-1)是FY0i的基准角FY01编号1→0°FY02编号2→40°。若锚点选FY05编号5基准角是160°不减这个值坐标必错。此步是精度核心——角度差0.1°100米距离上偏差达17cm。2.4 验证三圆共点画图比公式更直观光算出坐标不够必须验证是否真在三条圆弧上。MATLAB用plot画出三段圆弧及交点一目了然% 绘制三圆弧以FY00-FY0i, FY00-FY0j, FY0i-FY0j为弦张角为alpha1/2/3 figure; hold on; % 圆弧1FY00-FY0i为弦张角alpha1 center1 circle_center_from_chord_angle([0,0], [x_anchor(FY0i1),y_anchor(FY0i1)], alpha1); radius1 norm(center1 - [0,0]); fplot((t) center1(1)radius1*cos(t), (t) center1(2)radius1*sin(t), [theta_start1, theta_end1]); % 同理画圆弧2、3... plot(x_est, y_est, ro, MarkerSize, 10, LineWidth, 2); % 标出解点 title(三圆共点验证红点为定位解);避坑提示circle_center_from_chord_angle是自定义函数需提前写好。MATLAB无内置“弦张角→圆心”函数必须手推公式圆心在弦的中垂线上距弦中点距离为d (c/2)/tan(α/2)c为弦长α为张角。若直接调用circles工具箱会因坐标系不匹配导致圆心偏移。3. 锚点怎么选枚举法不是暴力穷举而是用表2筛出最小误差组合问题一第二问的灵魂在于只给定0号圆心和1号圆周首机为锚点其他8架机位置都有偏差此时仅靠两个方向信息无法唯一定位——因为单个夹角只能确定一段圆弧“葫芦”或“眼睛”两个夹角得两条圆弧最多两个交点。必须增加第三锚点让三条圆弧共点。但选谁随机试不要用数学筛。3.1 “葫芦”与“眼睛”钝角/锐角轨迹的本质区别先理解为什么一个夹角对应圆弧。设锚点A、B固定接收机K满足∠AKBα。当α90°锐角K的轨迹是两段凸向AB的圆弧形如葫芦当α90°钝角轨迹是两段凹向AB的圆弧形如眼睛。MATLAB用轨迹方程验证% 钝角轨迹α120°方程为 (x-xc)^2 (y-yc)^2 R^2圆心在AB中垂线外侧 alpha_dull 120; % 计算圆心AB中点M中垂线方向n距离d |AB|/(2*tan(alpha_dull/2)) A [0,0]; B [100,0]; % FY00,FY01 M (AB)/2; AB_vec B-A; n [-AB_vec(2), AB_vec(1)]; % 中垂线方向向量 d norm(AB_vec)/(2*tand(alpha_dull/2)); center_dull M d * n / norm(n); R_dull norm(center_dull - A); % 画钝角轨迹 theta linspace(0, 2*pi, 100); x_dull center_dull(1) R_dull*cos(theta); y_dull center_dull(2) R_dull*sin(theta); plot(x_dull, y_dull, b--);逻辑说明tand是度数正切norm求向量模长。钝角时tan(α/2)1d较小圆心靠近AB锐角时tan(α/2)1d较大圆心远离AB——这决定了轨迹形态。若误用tan弧度代替tandd计算错误轨迹完全失真。3.2 表2的生成逻辑56种组合不是乱列是锚点编号的笛卡尔积题目表2列出FY02~FY098架作为接收机与FY02~FY09中另一架排除自身作为第三锚点的所有组合共8×756种。MATLAB用嵌套循环生成% 生成表2所有(FY0j, FY0k)组合对应的三个角度alpha1,alpha2,alpha3 r 100; alpha_table zeros(8,8); % 行FY0j, 列FY0k for j 2:9 % FY0j编号 for k 2:9 % FY0k编号 if j k, continue; end % 排除自身 % 计算FY00,FY0j,FY0k构成三角形的三个内角 A [0,0]; B [x_anchor(j1), y_anchor(j1)]; C [x_anchor(k1), y_anchor(k1)]; a norm(B-C); b norm(A-C); c norm(A-B); % 对边 alpha1 acosd((b^2 c^2 - a^2)/(2*b*c)); % ∠BAC FY00处角 alpha2 acosd((a^2 c^2 - b^2)/(2*a*c)); % ∠ABC FY0j处角 alpha3 acosd((a^2 b^2 - c^2)/(2*a*b)); % ∠ACB FY0k处角 alpha_table(j-1,k-1) alpha1; % 存alpha1同理存alpha2,alpha3需三维数组 end end参数说明alpha_table(j-1,k-1)因数组索引从1开始FY0j2对应行1。表2中每个格子存的是α₁FY00处角实际应用需同时存α₂、α₃故代码中应声明alpha1_table,alpha2_table,alpha3_table三个8×8矩阵。3.3 最小二乘筛选为什么用平方和而不是绝对值接收机实测角度(β₁,β₂,β₃)总有噪声需找一组(FY0j,FY0k)使预测角与实测角最接近。最小二乘平方和比L1范数绝对值和更鲁棒——它惩罚大误差避免被单个离群角度带偏% 实测角度含噪声 beta1 45.3; beta2 59.8; beta3 74.5; % 遍历所有56种组合找最小平方和 min_error inf; best_j 0; best_k 0; for j 2:9 for k 2:9 if j k, continue; end error (beta1 - alpha1_table(j-1,k-1))^2 ... (beta2 - alpha2_table(j-1,k-1))^2 ... (beta3 - alpha3_table(j-1,k-1))^2; if error min_error min_error error; best_j j; best_k k; end end end fprintf(最优第三锚点FY0%d 和 FY0%d\n, best_j, best_k);避坑提示alpha1_table(j-1,k-1)索引必须与循环变量一致。曾有队伍将j循环写成1:8但表2中FY02对应行1导致选错锚点。此外inf初始化确保首次比较必更新若用0则永远不进if。3.4 常见问题排查三类典型翻车现场现象 → 原因 → 解决① 表2生成的角度全是NaN→ 原因余弦定理中(b²c²-a²)/(2bc)超出[-1,1]范围因浮点误差导致acosd输入非法。→ 解决加截断保护val max(-1, min(1, (b^2c^2-a^2)/(2*b*c))); alpha acosd(val);② 最小二乘选出的锚点定位结果偏差超10米→ 原因实测角度β₁,β₂,β₃中有一个严重离群如传感器瞬时故障但平方和对离群值敏感。→ 解决改用中位数绝对偏差MAD预筛或对56组误差排序后取前5名再人工看哪组几何合理性高如三锚点不共线。③ 同一组实测角度不同MATLAB版本结果不同→ 原因solve符号求解在R2020a后默认启用数值近似解可能不唯一。→ 解决强制用Real, true选项sol solve([eq1,eq2,eq3], [k0,a2,a3], Real, true);并检查sol.k0是否为实数数组。4. 位置怎么调0.01米步长遍历不是蛮干是带方向约束的梯度搜索问题一第三问要求“从略有偏差的初始位置调整到理想位置”。表面看是暴力遍历以当前点为中心0.01米步长扫一圈找离理想点最近的。但实际远不止——无人机不能瞬移调整必须沿可行路径且“略有偏差”意味着初始点已在理想点附近遍历范围可大幅压缩。4.1 步长设计的物理意义0.01米不是拍脑袋是精度与算力的平衡题目指定步长0.01米背后有深意圆形编队半径100m角度分辨率需达0.0057°100m×0.01m弧长对应角度这恰好匹配民用GPS水平精度2.5m经差分后的提升量。MATLAB实现时步长决定计算量% 初始位置极坐标row距离theta角度 row_init 102.5; % 比理想100m远2.5m theta_init 42.3; % 比理想40°多2.3° % 理想位置FY02对应theta40° row_ideal 100; theta_ideal 40; % 遍历范围半径±5m角度±5°步长0.01m半径和0.01°角度 row_range row_init (-5:0.01:5); theta_range theta_init (-5:0.01:5); min_dist inf; best_row row_init; best_theta theta_init; for r row_range for t theta_range x_curr r * cosd(t); y_curr r * sind(t); dist sqrt((x_curr - row_ideal*cosd(theta_ideal))^2 ... (y_curr - row_ideal*sind(theta_ideal))^2); if dist min_dist min_dist dist; best_row r; best_theta t; end end end fprintf(调整后位置row%.3f m, theta%.3f°\n, best_row, best_theta);逻辑说明双循环遍历虽直观但计算量大1000×100010⁶次。实际可优化先粗调步长0.1m/0.1°找到最小值区域再精调0.01m/0.01°。代码中(-5:0.01:5)生成1001个点是精度保障。4.2 方向约束为什么只分“向北/向南”两种调整题目图10-11将调整分为向北/向南本质是极坐标中的径向ρ与角向θ解耦。但真实飞行中无人机受动力学约束不能纯径向或纯角向移动——必须合成速度矢量。MATLAB用向量加法模拟% 当前速度矢量假设朝向理想点 dx row_ideal*cosd(theta_ideal) - row_init*cosd(theta_init); dy row_ideal*sind(theta_ideal) - row_init*sind(theta_init); speed sqrt(dx^2 dy^2); % 若speed 0.01m则归一化后乘以0.01m作为单步位移 if speed 0.01 dx_step 0.01 * dx / speed; dy_step 0.01 * dy / speed; row_new sqrt((row_init*cosd(theta_init)dx_step)^2 ... (row_init*sind(theta_init)dy_step)^2); theta_new atan2d(row_init*sind(theta_init)dy_step, ... row_init*cosd(theta_init)dx_step); else row_new row_init; theta_new theta_init; % 已在误差内 end参数说明atan2d(y,x)是四象限反正切避免atand(y/x)在x0时崩溃。dx/dy是理想位移矢量dx_step/dy_step是单位步长保证每次移动严格0.01m。4.3 迭代收敛判断0.01米阈值不是终点是动态容差“当与准确位置相差0.01米时就近似认为到达”是理想条件。实际中因传感器噪声距离可能在0.009~0.011m间震荡。MATLAB用滑动窗口判断收敛% 初始化历史距离数组 dist_history zeros(1, 10); % 存最近10次距离 tolerance 0.01; max_iter 100; for iter 1:max_iter % ... 计算当前dist ... dist_history [dist_history(2:end), dist]; % 若最近10次距离标准差 0.001m且均值 tolerance则收敛 if std(dist_history) 0.001 mean(dist_history) tolerance fprintf(迭代%d次后收敛最终距离%.4f m\n, iter, mean(dist_history)); break; end end避坑提示仅判断单次距离0.01m易受噪声干扰。用滑动窗口统计稳定性是工程落地的关键技巧。5. 锥形编队怎么破三轮定位不是流程是规避“延长线病态”的生存策略锥形编队问题二看似只是形状变化实则引入致命病态当接收机位于三锚点构成三角形的边延长线上时三个夹角中必有一个为0°或180°导致余弦定理分母为零定位失效。这就是为什么题目明确说“第二轮中4、7、9、11、13、14号无人机依然无法定位”——它们恰在1/2/5号机的延长线上。5.1 延长线病态的数学本质退化三角形与奇异解设锚点A、B、C接收机K在AB延长线上。此时∠AKB0°K在A、B之间或180°K在AB外侧cos(0°)1cos(180°)-1。代入余弦定理(k₀²a₂²-r²)/(2*k₀*a₂)±1→(k₀∓a₂)²r²解不唯一k₀a₂±r。MATLAB检测方法% 检测K是否在AB延长线上向量KA与KB共线且反向 A [x_anchor(2), y_anchor(2)]; % FY01 B [x_anchor(6), y_anchor(6)]; % FY05 K [x_est, y_est]; % 待检点 KA A - K; KB B - K; % 共线叉积≈0反向点积0 cross_prod KA(1)*KB(2) - KA(2)*KB(1); dot_prod KA(1)*KB(1) KA(2)*KB(2); if abs(cross_prod) 1e-6 dot_prod 0 warning(K点位于AB延长线上定位病态); % 启用备用方案换锚点或加权平均 end逻辑说明1e-6是浮点容差dot_prod0确保反向延长线外侧。若仅用cross_prod≈0会误判K在AB线段上此时点积0应可定位。5.2 三轮定位的底层逻辑每轮扩大“健康锚点集”第一轮用1号、5号题目假设有两架无偏差定位2号 → 2号加入锚点集现3架健康锚点。第二轮用1/2/5号定位其余但4/7/9等在延长线上者失败 → 剩余健康锚点仍为3架。第三轮选2/6/8号为新锚点题目指定因其构成三角形不包含前述病态线 → 所有剩余机均可定位。MATLAB用集合操作模拟% 健康锚点集初始 healthy_anchors [1, 5]; % FY01, FY05 % 第一轮定位FY02 [x2,y2] locate_by_angles(healthy_anchors, [beta1,beta2,beta3]); % 自定义定位函数 if ~is_nan(x2) % 定位成功 healthy_anchors [healthy_anchors, 2]; % 加入FY02 end % 第二轮尝试定位FY03~FY09 for k 3:9 [xk,yk] locate_by_angles(healthy_anchors, get_betas(k)); % 获取k号机实测角 if ~is_nan(xk) ~on_extension_line([xk,yk], healthy_anchors) % 不在延长线 healthy_anchors [healthy_anchors, k]; end end % 第三轮强制用[2,6,8]为锚点 final_anchors [2,6,8]; for k 3:9 if ~ismember(k, healthy_anchors) % 未定位者 [xk,yk] locate_by_angles(final_anchors, get_betas(k)); % 调整... end end参数说明on_extension_line是前述延长线检测函数。ismember检查是否已在健康集中避免重复定位。5.3 锥形编队坐标预置等边三角形阵列的硬编码陷阱锥形编队非圆形题目说“直线上相邻间距相等”并配图13示意等边三角形。MATLAB需重置锚点坐标% 锥形编队假设5行每行机数1/2/3/2/1间距d50m d 50; % 行1顶点FY01 at (0,0) x_cone(1) 0; y_cone(1) 0; % 行2FY02,FY03 at yd, x±d/2 x_cone(2) -d/2; y_cone(2) d; x_cone(3) d/2; y_cone(3) d; % 行3FY04~FY06 at y2d, x-d,0,d x_cone(4) -d; y_cone(4) 2*d; x_cone(5) 0; y_cone(5) 2*d; x_cone(6) d; y_cone(6) 2*d; % ... 依此类推避坑提示题目未给锥形具体参数此坐标系是合理假设。若直接套用圆形坐标所有锥形计算全错。必须在脚本开头用% CONE_MODE 1标志切换模式。5.4 常见问题排查锥形编队的三大隐形雷区现象 → 原因 → 解决① 锥形定位结果全为NaN但圆形正常→ 原因锥形锚点坐标未重置仍用圆形x_anchor数组导致距离计算错误。→ 解决在locate_by_angles函数入口加断言assert(all(x_anchor x_cone) || all(x_anchor x_circle))。② 第二轮定位后healthy_anchors长度始终为3→ 原因on_extension_line检测过于宽松将本可定位的点误判为病态。→ 解决放宽cross_prod容差至1e-4并增加角度验证若∠AKB5°或175°才判病态。③ 第三轮用[2,6,8]锚点但FY04仍在延长线上→ 原因2/6/8坐标未形成非退化三角形三点共线。→ 解决在选锚点前用polyarea(x([2,6,8]), y([2,6,8])) 1e-6验证面积否则换组合。6. 调位脚本怎么验用三组数据打穿所有边界别信“跑通就行”写完脚本最危险的心态是“输出没报错图也出来了应该没问题”。纯方位定位的坑往往藏在边界案例里钝角葫芦、锐角眼睛、延长线病态、浮点溢出。我当年调试时用三组数据反复锤炼脚本至今保留着验证习惯。6.1 边界数据集设计覆盖所有几何病态数据组场景关键参数验证目标Case A钝角主导FY00,FY01,FY05为锚点FY02接收α₁150°, α₂20°, α₃10°检查“眼睛”轨迹是否收敛acosd输入是否截断Case B锐角主导FY00,FY02,FY03为锚点FY04接收α₁30°, α₂40°, α₃50°检查“葫芦”轨迹是否分裂solve是否返回多解Case C延长线病态FY01,FY02,FY03共线FY04在延长线上A(0,0), B(50,0), C(100,0), K(150,0)检查on_extension_line是否触发是否跳过定位MATLAB批量验证脚本test_cases { A, [150,20,10], [0,1,5], 2; % {name, angles, anchors, receiver} B, [30,40,50], [0,2,3], 4; C, [0,0,0], [1,2,3], 4; % α0°模拟病态 }; for i 1:size(test_cases,1) case_name test_cases{i,1}; betas test_cases{i,2}; anchors test_cases{i,3}; receiver test_cases{i,4}; fprintf(\n 测试 %s \n, case_name); try [x,y] locate_by_angles(anchors, betas); if is_nan(x) || is_nan(y) fprintf(%s: 定位失败预期\n, case_name); p a hrefhttps://download.csdn.net/download/jshsjjn/87937479 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p