简介:本资源为2022年全国大学生数学建模竞赛本科组B题的完整论文PDF,聚焦无人机纯方位无源定位问题,适合备战数模竞赛的本科生及指导教师参考。论文围绕圆形与锥形两种编队,系统建立被动接收信号无人机的定位模型与位置调整模型,运用二维直角坐标系、余弦定理与“三圆共点”原理,推导发射信号无人机数量需求,并给出基于0.01米步长的遍历调整方案,全部求解过程均以MATLAB编程实现。压缩包内仅含1个PDF文件,约2.63MB,完整呈现摘要、问题重述、模型假设、符号说明及建模求解全过程,便于读者直接研读建模思路与代码逻辑。该论文曾获河北省一等奖,目前已有8039人学习下载,可为参赛者提供从模型构建到编程求解的完整参考范例。
1. 从三组方位角反推坐标:这套无人机纯方位无源定位方案到底能跑通什么
编队飞行里最怕的不是风,是电磁信号被干扰后整个队形散架。纯方位无源定位的思路很直接:无人机不主动发波,只靠接收其他已知位置无人机发来的方向信息,反推自己在哪、该往哪挪。这套 2022 年全国大学生数学建模竞赛本科组 B 题的完整解答,核心就是把这个过程拆成可计算的模型——用二维直角坐标系加余弦定理,把三个方位角变成坐标,再用步长遍历把偏差位置拉回理想位置。它适合两类人:正在备战数模竞赛、需要一套能直接跑 MATLAB 的完整代码框架的选手;以及做无人机编队控制、想理解纯方位定位最小可解条件的工程师。整套资料包含问题重述、模型建立、MATLAB 程序附录和结果分析,不是纯理论推导,附录里四个完整程序覆盖了从定位到调整的全流程。
2. 定位模型怎么建:从三个方位角到二维坐标的完整推导链
2.1 为什么选余弦定理而不是最小二乘
纯方位无源定位的常见做法有两种:一种是基于最小二乘的伪线性估计,把非线性角度方程线性化后迭代求解;另一种是直接利用几何关系构造解析解。这套方案选了后者,原因很实际——圆形编队里发射信号的无人机位置是精确已知的,接收信号的无人机与它们之间的角度关系可以用余弦定理直接串起来,不需要迭代,也不存在初值敏感的问题。
具体来说,把第 i 架无人机的准确位置固定在 x 轴上,圆形编队半径记为 r,相邻无人机间隔 40 度。三架发射信号无人机分别记为 FY0i、FY0j、FY0k,接收信号无人机为 FY0x。接收到的三个方位角 α1、α2、α3 分别对应接收机与三架发射机连线之间的夹角。在三角形中反复使用余弦定理,就能把三条边的长度解出来,再反解出接收机的坐标。
这里有个容易翻车的地方:当接收机位于发射机连线的延长线上时,角度信息退化,方程会出现多解甚至无解。这也是后面问题二里第二轮定位时 4、7、9、11、13、14 号无人机无法定位的根本原因。
2.2 被动接收信号无人机定位模型的 MATLAB 实现
附录一的程序分两段:第一段求解当前位置坐标,第二段求解准确位置坐标。先看第一段的核心逻辑。
% 被动接收信号无人机定位模型的建立 % 求解当前位置的坐标 alpha1 = input('第一个角度:'); % 0k i alpha2 = input('第二个角度:'); % 0k j alpha3 = input('第三个角度:'); % j k i FY0i = input('发射信号无人机 1:'); FY0j = input('发射信号无人机 2:'); r = input('圆形编队半径:'); % 计算两架发射信号无人机之间的距离 a1 = sqrt(2*r^2 - cosd(abs(FY0i - FY0j)*40)*2*r^2); syms k0 a2 a3 % 三个余弦定理方程,三个未知量 eq1 = (k0^2 + a3^2 - r^2)/(2*k0*a3) - cosd(alpha1); eq2 = (k0^2 + a2^2 - r^2)/(2*k0*a2) - cosd(alpha2); eq3 = (a2^2 + a3^2 - a1^2)/(2*a2*a3) - cosd(alpha3); [X, Y, Z] = solve(eq1, eq2, eq3, k0, a2, a3); K0 = double(X); % 求得 k0 长度 % 计算接收机相对于 x 轴的偏转角 k0i = acosd(1 - 2*cosd(alpha1)^2); k01 = k0i - 40*abs(i - 1); x = real(K0*cosd(k01)) % 横坐标 y = real(K0*sind(k01)) % 纵坐标这段代码的关键参数有三个:alpha1、alpha2、alpha3是接收到的三个方位角,单位是度;FY0i和FY0j是两架发射信号无人机的编号,用来计算它们之间的基准距离;r是圆形编队半径,题目里是 100 米。syms声明符号变量后,solve一次性解出三个未知边长,再用double转成数值。最后两行把极坐标转成直角坐标,real是防止符号计算产生微小虚部。
第二段程序反过来,已知接收机编号和发射机编号,直接算出三个方位角的准确值:
% 求解无人机接收到的角度信息 FY00 = 0; FY0i = input('发射信号无人机 1:'); FY0j = input('发射信号无人机 2:'); FY0k = input('接收信号无人机:'); r = input('圆形编队半径:'); % 计算三条边的长度 a1 = sqrt(2*r^2 - cosd((FY0j - FY0i)*40)*2*r^2); a3 = sqrt(2*r^2 - cosd((FY0k - FY0i)*40)*2*r^2); a2 = sqrt(2*r^2 - cosd((FY0k - FY0j)*40)*2*r^2); % 用余弦定理求三个角 alpha1 = acosd((r^2 + a3^2 - r^2)/(2*r*a3)); alpha2 = acosd((r^2 + a2^2 - r^2)/(2*r*a2)); alpha3 = acosd((a3^2 + a2^2 - a1^2)/(2*a3*a2));注意alpha1和alpha2的公式里出现了r^2 + a3^2 - r^2,看起来像冗余,实际上是因为圆心无人机 FY00 到圆周无人机的距离恒为 r,代入余弦定理后 r 平方项抵消了一部分。这是圆形编队的特殊性质,换成锥形编队就不能这么简化。
提示:运行第一段程序前,需要先确定接收机在圆周上的大致方位,否则
k01的计算会引入 40 度整数倍的偏移误差。我一般会先用第二段程序算出准确角度,再拿第一段反解验证。
3. 额外需要几架发射机:枚举法确定最小可解条件
3.1 为什么一架额外发射机就够了
问题一第二小问的核心结论是:除 0 号和 1 号外,只需要额外增加 1 架发射信号无人机,就能实现有效定位。这个结论的推导依赖一个几何事实——一个动点到两个定点的夹角为定值时,动点轨迹是上下对称的两段圆弧组成的封闭曲线。当夹角为钝角时,轨迹像一只“眼睛”;当夹角为锐角时,轨迹像一个“葫芦”。
三架发射机产生三条这样的轨迹,三条轨迹相交于一点,该点就是接收机的实际位置。如果只有两架发射机,接收机只能得到一个方位角,轨迹是一条曲线而非一个点,无法唯一确定位置。所以额外增加一架是理论下限。
但这里有个隐藏条件:三条轨迹必须确实相交于同一点。当接收机位于某两架发射机连线的延长线上时,对应的轨迹退化成直线,交点不唯一。这就是为什么在问题二的锥形编队里,第二轮定位时会有 6 架无人机无法定位。
3.2 枚举法求解额外发射机编号的 MATLAB 实现
附录二的程序用枚举法遍历所有可能的发射机组合,找出使方位角偏差平方和最小的那一组。
% 求出所有的可能性 FY00 = 0; FY0i = 1; % 两个已知的定点 r = input('圆形编队半径:'); % 遍历所有可能的发射机组合,预计算角度表 for FY0j = 2:9 for FY0k = 2:9 if FY0k == FY0j continue; end a1 = sqrt(2*r^2 - cosd((FY0j - FY0i)*40)*2*r^2); a3 = sqrt(2*r^2 - cosd((FY0k - FY0i)*40)*2*r^2); a2 = sqrt(2*r^2 - cosd((FY0k - FY0j)*40)*2*r^2); alpha1(FY0j, FY0k) = acosd((r^2 + a3^2 - r^2)/(2*r*a3)); alpha2(FY0j, FY0k) = acosd((a3^2 + a2^2 - a1^2)/(2*a3*a2)); alpha3(FY0j, FY0k) = acosd((r^2 + a2^2 - r^2)/(2*r*a2)); end end % 反解:根据接收到的角度信息找最匹配的发射机 beta1 = input('第一个角度:'); beta2 = input('第二个角度:'); beta3 = input('第三个角度:'); FY0k = input('接收信号的无人机编号:'); min_val = 1000; mini = 0; for i = 2:9 if i == FY0k continue; end m = (beta1 - alpha1(i, FY0k))^2 + ... (beta2 - alpha2(i, FY0k))^2 + ... (beta3 - alpha3(i, FY0k))^2; if m <= min_val min_val = m; mini = i; end end mini % 输出需要额外增加的发射信号无人机编号这段代码的逻辑分两步:先建一张角度查找表,把所有可能的发射机组合对应的三个方位角算出来存进alpha1、alpha2、alpha3三个矩阵;然后对于给定的接收角度beta1、beta2、beta3,遍历所有候选发射机,找平方和最小的那个。
参数说明:FY0i固定为 1,因为题目规定 0 号和 1 号已经发射信号;FY0j和FY0k从 2 到 9 遍历,跳过相同编号;m是三个角度偏差的平方和,越小说明匹配度越高。最终输出的mini就是需要额外增加的发射机编号。
注意:枚举法在 8 架无人机里找 1 架额外发射机,计算量是 C(7,1)=7 种组合,非常小。但如果问题扩展到几十架无人机,这种暴力枚举就不合适了,需要换成基于几何约束的剪枝搜索。
4. 位置调整怎么收敛:步长遍历与极坐标更新
4.1 向北还是向南:调整方向的分类逻辑
问题一第三小问要求把有偏差的无人机调整到理想位置。方案设定步长为 0.01 米,向四周全方位遍历,当与准确位置相差小于 0.01 米时认为到达理想位置。调整方向分成向北和向南两类,本质上是极坐标下角度增量的正负。
具体做法:无人机当前位置记为 B,向北移动一个微小步长后到达 C。已知 B 点的方位信息,可以算出 BC 的长度和角度变化,再用余弦定理求出 C 点相对于圆心的极坐标。重复这个过程直到偏差小于阈值。
这里有个工程上很实用的简化:调整时间不计,意味着无人机知道自己需要调整位置的瞬间就已经到位。这在实际系统里当然不成立,但在建模阶段可以大幅简化动力学部分,把注意力集中在定位精度上。
4.2 圆形编队位置调整的 MATLAB 实现
附录三的程序用 while 循环实现迭代调整。
% 从无人机接收的方位信息去定位无人机,找到无人机的位置 FY00 = 0; FY01 = 1; njl = 0.01; % 最小距离阈值 r = 100; % 圆形编队半径 bzrow = 100; % 准确位置 row bztheta = [0 40 80 120 160 200 240 280 320]; % 准确位置 theta FY0irow(2) = input('FY02 的初始位置(row):'); FY0itheta(2) = input('FY02 的初始位置(theta):'); jl = 100; while(jl > njl) mjl = 10000; for j = 0:359 % 遍历所有可能的角度,找距离最小的方向 % 此处省略具体角度计算,核心是极坐标更新 % ... end % 更新位置,缩小 jl jl = mjl; end这段代码的框架是:外层 while 控制迭代直到距离小于阈值,内层 for 遍历 0 到 359 度所有方向,找使距离最小的调整角度。njl是收敛阈值,设成 0.01 米;r是编队半径;bztheta是 9 架圆周无人机的理想角度分布,间隔 40 度。
实际运行时,每次迭代需要根据当前极坐标和方位角信息,用余弦定理算出移动一个步长后的新极坐标,再与理想位置比较。如果偏差没有缩小,说明方向选反了,需要反向遍历。
提示:步长 0.01 米在 100 米半径的编队里相当于万分之一精度,迭代次数可能上百次。如果只做仿真验证,可以先用 0.1 米步长快速收敛,最后再用 0.01 米精调。
5. 锥形编队怎么套用:从圆形到锥形的模型迁移与边界
5.1 锥形编队的几何差异与定位策略
问题二的锥形编队和圆形编队本质区别在于:圆周上相邻无人机等距,但锥形编队里相邻三架无人机构成等边三角形,直线上的间距相等。这意味着圆形编队里“相邻间隔 40 度”的简化不再成立,需要重新建立坐标系。
方案假设 1 号和 5 号无人机位置无偏差,先定位 2 号。具体做法:将 1 号、2 号作为发射机,5 号接收,得到一个方向信息;再将 2 号、5 号作为发射机,1 号接收,得到另一个方向信息。在三角形 1-2-5 中,两个角已知,第三个角随之确定,2 号的位置就解出来了。
但这里有个硬边界:当某架无人机位于 1、2、5 号构成的三角形的边的延长线上时,只能传送一个方向信息,无法准确定位。所以第二轮定位时,4、7、9、11、13、14 号无人机依然无法确定位置,只能先定位 3、6、8、10、12、15 号。第三轮选取 2、6、8 为发射机,剩下的未定位无人机才都能解算。
5.2 锥形编队位置调整的 MATLAB 实现要点
附录四的程序结构与附录三类似,但坐标系和理想位置表不同。锥形编队的理想位置表如下:
| 无人机编号 | X 坐标 | Y 坐标 |
|---|---|---|
| 1 | 100√3 | 50 |
| 2 | 50√3 | 50 |
| 3 | 0 | 0 |
| 4 | -50√3 | 50 |
| 5 | -100√3 | 100 |
| 6 | -100√3 | 0 |
| 7 | -100√3 | -100 |
| 8 | -50√3 | -50 |
| 9 | -50√3 | 50 |
| 10 | -50√3 | -150 |
| 11 | -100√3 | 200 |
| 12 | -100√3 | 100 |
| 13 | -100√3 | 0 |
| 14 | -100√3 | -100 |
| 15 | -100√3 | -200 |
调整逻辑与圆形编队一致:设定步长,向四周遍历,直到偏差小于阈值。区别在于锥形编队的理想位置不是均匀分布在圆周上,而是按等边三角形网格排列,所以极坐标更新时需要先转成直角坐标再计算距离。
% 锥形编队位置调整核心逻辑 target_x = -100*sqrt(3); % 示例:15 号机目标 X target_y = -200; % 示例:15 号机目标 Y step = 0.01; current_x = input('当前 X 坐标:'); current_y = input('当前 Y 坐标:'); while sqrt((current_x - target_x)^2 + (current_y - target_y)^2) > step % 计算当前距离 dist = sqrt((current_x - target_x)^2 + (current_y - target_y)^2); % 沿梯度方向移动一个步长 dx = (target_x - current_x) / dist * step; dy = (target_y - current_y) / dist * step; current_x = current_x + dx; current_y = current_y + dy; end这段代码比圆形编队的极坐标遍历更直接:已知目标坐标,直接沿梯度方向走步长。参数step控制精度,target_x和target_y从理想位置表里查。实际运行时,每架无人机需要先通过纯方位定位确定自己的当前位置,再查表找到目标位置,然后执行梯度下降。
注意:锥形编队里 4、7、9、11、13、14 号无人机在第二轮无法定位,必须等第三轮 2、6、8 号调整到位后才能解算。如果强行在第二轮对这些无人机执行调整,会因为位置估计错误导致越调越偏。
6. 从跑通到跑对:四个容易翻车的参数与验证习惯
第一个坑是角度单位。MATLAB 的acosd返回的是度,但cos和sin默认是弧度。附录一的程序里混用了cosd和cos,如果输入的角度没有统一,算出来的坐标会差一个 π/180 的系数。我一般会在程序开头加一行clear; clc;,然后把所有三角函数统一成带 d 的版本。
第二个坑是符号计算的虚部。solve解方程组时可能产生复数解,double转换后虚部虽然很小但不为零。real函数能去掉虚部,但如果虚部量级超过 1e-6,说明方程本身有问题,需要检查角度输入是否合理。
第三个坑是枚举法的边界。附录二里FY0j和FY0k从 2 到 9 遍历,但alpha1矩阵的索引是(FY0j, FY0k),当FY0j或FY0k等于FY0k(接收机编号)时,对应的角度值没有物理意义。程序里用continue跳过了相同编号,但没有跳过接收机编号,实际使用时需要手动排除。
第四个坑是步长与收敛速度的权衡。0.01 米步长在 100 米半径下需要约 10000 次迭代才能从 100 米外收敛到目标,每次迭代都要重新计算方位角,计算量不小。我一般会先用 1 米步长粗调,再用 0.01 米精调,总迭代次数能降到十分之一。
验证方法很简单:跑完附录一算出坐标后,代回附录二反解发射机编号,如果得到的编号与输入一致,说明定位模型自洽。再用附录三调整位置,检查最终坐标与理想位置表的偏差是否小于 0.01 米。三个程序串起来跑通,这套纯方位无源定位方案就算真正落地了。
从那以后我每次拿到新的编队模型,都强制先跑一遍“正解-反解-调整”的闭环验证,确认三个环节的数值能对上再往下做。希望帮到你。
本文还有配套的精品资源,点击获取