ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

2022数模B题无人机纯方位无源定位:MATLAB代码复现与避坑指南

2022数模B题无人机纯方位无源定位:MATLAB代码复现与避坑指南

简介:本资源为2022年全国大学生数学建模竞赛本科组B题的完整论文PDF,聚焦无人机纯方位无源定位问题,适合备战数模竞赛的本科生及指导教师参考。论文围绕圆形与锥形两种编队,系统建立被动接收信号无人机的定位模型与位置调整模型,涵盖余弦定理反解、三圆共点轨迹分析、枚举法确定发射信号机数量,以及基于步长遍历的迭代调整策略,并附MATLAB编程实现思路。压缩包内仅含1个PDF文件,约2.63MB,内容完整呈现摘要、问题重述、模型假设、符号说明及建模求解全过程,便于读者对照赛题复盘建模逻辑与算法细节。该论文曾获河北省一等奖,目前已有8039人学习下载,适合需要理解几何定位建模、学习MATLAB数值求解与优化迭代思路的读者深入研读。

1. 从三组方位角反推坐标:这套数模B题资源到底能跑出什么

如果你手头正好有一套2022年全国大学生数学建模竞赛本科组B题的完整论文和MATLAB代码,但打开一看全是公式推导、附录代码又长又碎,不知道从哪跑起——这篇就是替你拆它的。这套资源围绕无人机纯方位无源定位展开,核心解决三件事:用三个方位角反解接收信号无人机的实际坐标、判断最少需要几架无人机同时发射信号才能唯一定位、以及在圆形和锥形两种编队下把有偏差的无人机一步步调整回理想位置。它适合正在备战数模竞赛、想找一套能跑通的几何定位代码做参考的参赛者,也适合做多站无源定位入门、想看看余弦定理怎么在工程里落地的工程师。整套代码用MATLAB写成,符号求解和数值遍历混用,思路清晰但细节坑不少,下面按能复现的顺序拆开讲。

2. 定位模型怎么建:从余弦定理到三圆共点的MATLAB落地

2.1 为什么选余弦定理而不是最小二乘

纯方位无源定位的常见做法有两类:一类是纯数值优化,把位置当未知量、角度当观测值,直接上最小二乘或卡尔曼滤波;另一类是几何解析法,利用已知的编队结构(圆心加圆周均匀分布)把角度关系写成方程,用余弦定理逐步反解。这套资源走的是第二条路,原因很实际——圆形编队的几何约束太强了,九架无人机均匀分布在半径r的圆周上,相邻夹角固定40度,圆心那架位置已知,这些条件直接能把未知量压缩到可解析求解的程度。用最小二乘反而把结构信息浪费了,还容易陷入局部最优。

具体做法是:把第i架无人机的准确位置固定在x轴上,接收信号无人机与两架发射信号无人机之间的三个夹角α1、α2、α3作为已知量,通过余弦定理建立三边长度与角度之间的方程组。这里有个关键处理——三个未知边长k0、a2、a3对应三个余弦方程,正好凑成恰定方程组,用MATLAB的符号求解器solve直接解。

% 被动接收信号无人机定位模型 - 求解当前位置坐标 alpha1 = input('第一个角度:'); % 0kj夹角 alpha2 = input('第二个角度:'); % 0ki夹角 alpha3 = input('第三个角度:'); % jki夹角 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长度 k0i = acosd(1 - 2*cosd(alpha1)^2); k01 = k0i - 40*abs(i-1); x = real(K0*cosd(k01)); % 横坐标 y = real(K0*sind(k01)); % 纵坐标

这段代码的逻辑链条是:先用两架发射信号无人机的编号差算出它们之间的弦长a1,再把三个余弦定理方程写成符号表达式,solve返回三边长度,最后通过角度关系换算成直角坐标。参数上要注意FY0i和FY0j的编号差决定了a1的大小,编号差1对应40度、差2对应80度,以此类推。r是编队半径,题目里圆形编队取100米,但调试时可以先取10米快速验证。

注意:solve返回的符号解可能有多个分支,double转换后要取实部,因为余弦定理在几何上只对应实数解。如果出现虚部,说明输入的三个角度不满足三角形闭合条件,检查角度是否来自同一架接收机。

2.2 三圆共点怎么在代码里实现

第二小问的核心是确定额外需要几架无人机发射信号。资源里的思路是:只有0号和1号发射时,接收机只能拿到一个方位角,轨迹是一条“葫芦”或“眼睛”形曲线,无法唯一定位;再增加一架发射机,接收机拿到三个方位角,三条轨迹交于一点,这就是三圆共点。

代码实现上,先枚举所有可能的发射机组合,算出每种组合下三个方位角的理论值,存成表;然后对于给定的接收机实测角度,遍历所有组合,找理论值与实测值平方和最小的那组。

% 枚举所有发射机组合,建立角度查找表 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的索引顺序。代码里alpha1(FY0j,FY0k)存的是以0号、j号、k号构成的三角形中,0号对应的角;alpha2存的是j号对应的角;alpha3存的是k号对应的角。如果索引搞反,后面反解时平方和最小对应的编号就会错。我一般会在建表后先手动验算一组:比如FY0j=2、FY0k=3时,三个角应该分别是40度、40度、100度附近,对不上就说明索引写错了。

提示:枚举法算出的角度表是56种组合(8架接收机各对应7种发射机组合),实际调试时可以先只算前10种,确认逻辑正确再跑全量。

3. 位置调整怎么迭代:步长0.01米的遍历策略与收敛判断

3.1 向北还是向南:调整方向的分类逻辑

第三小问要求把有偏差的无人机调整到理想位置。资源里的做法很朴素但有效:设定步长0.01米,以当前位置为中心向四周全方位遍历,每走一步就重新计算接收到的方位角,与理想位置的理论方位角比较,当距离误差小于0.01米时认为到位。

调整方向分成向北和向南两类,判断依据是当前接收到的方位角与理想方位角的偏差方向。代码里用极坐标表示位置,每次移动后更新极径和极角。

% 圆形编队位置调整 - 迭代遍历 FY00 = 0; FY01 = 1; njl = 0.01; % 最小距离阈值 r = 100; % 圆形编队半径 bztheta = [0 40 80 120 160 200 240 280 320]; % 理想位置角度 FY0irow(2) = input('FY02的初始位置(row):'); FY0itheta(2) = input('FY02的初始位置(theta):'); jl = 100; while(jl > njl) mjl = 10000; for j = 0:359 % 在当前点周围遍历,计算每个方向上的方位角偏差 % 具体遍历逻辑:以当前极坐标为中心,步长0.01米向四周试探 % 找到使方位角偏差最小的移动方向 % ... end % 更新位置,重新计算jl jl = abs(FY0irow(2) - r); % 与理想半径的偏差 end

这段代码的骨架是while循环加内层for遍历。外层while控制迭代直到误差小于阈值,内层for遍历360个方向找最优移动方向。参数njl=0.01是收敛阈值,设太小会迭代次数爆炸,设太大定位精度不够。r=100是编队半径,调试时可以改成10加快收敛。

实际跑的时候有个坑:内层for遍历360个方向,每个方向都要算一次余弦定理,单次迭代计算量不小。如果初始偏差大,可能要迭代上千次。我一般会先用粗步长(比如1米)快速逼近,再用细步长(0.01米)精调,这样能省一半以上时间。

3.2 锥形编队的定位顺序与延长线问题

问题二把圆形编队换成锥形编队,几何结构变了但定位原理不变。资源里的处理是:先假设1号和5号位置无偏差,用它们定位2号,再把1、2、5作为已知点去定位其他无人机。

锥形编队的特点是相邻三架无人机构成等边三角形,直线上的间距相等。定位2号时,先用1号和2号发射、5号接收得到一个方位角,再用2号和5号发射、1号接收得到另一个方位角,三角形中两个角已知,第三个角和边长就能推出来。

% 锥形编队 - 定位2号无人机 % 第一步:1号、2号发射,5号接收 % 第二步:2号、5号发射,1号接收 % 两个方位角确定后,解三角形得到2号位置 % 然后向四周遍历调整到理想位置 % 第二轮定位时,在延长线上的无人机无法定位 % 4、7、9、11、13、14号在1-2-5三角形边的延长线上,跳过 % 第三轮选取2、6、8为发射机,剩余无人机均可定位

这里的关键判断是“是否在三角形边的延长线上”。如果在延长线上,接收到的方位角只有一个方向信息,无法构成三圆共点条件。代码里通过判断接收机与发射机连线的夹角是否接近0度或180度来识别,接近就跳过。

注意:锥形编队的理想位置表在资源里是以坐标形式给出的,调试时先把这些坐标转成极坐标或直角坐标存成数组,避免每次手动输入。

4. 避坑与排查:这套代码跑不起来时先查这五处

4.1 现象:solve返回空解或虚数解

原因:输入的三个角度不满足三角形内角和约束,或者角度单位搞混了。MATLAB的cosd接受角度制,cos接受弧度制,代码里混用会导致方程无解。

解决:统一用cosd和acosd,输入角度前先验算α1+α2+α3是否在合理范围内(对于圆形编队,三个角之和通常在100到200度之间)。如果solve返回空,把三个方程单独拿出来手算验证。

4.2 现象:枚举法找出的发射机编号与预期不符

原因:角度查找表的索引顺序与反解时的索引顺序不一致。建表时alpha1(FY0j,FY0k)存的是0号对应的角,反解时如果按alpha1(FY0k,FY0j)去查,就会错位。

解决:在建表后加一行验证代码,手动算一组已知组合的三个角,打印出来与理论值对比。确认索引方向后再跑反解。

4.3 现象:位置调整迭代不收敛,while循环死循环

原因:步长设得太小(比如0.001米)但初始偏差很大(几十米),每次移动距离太小,迭代次数超过预期。或者收敛阈值njl设得比步长还小,永远达不到。

解决:先用大步长快速逼近,再切小步长精调。收敛阈值至少是步长的2到3倍。加一个最大迭代次数保护,比如while(jl > njl && iter < 10000)。

4.4 现象:锥形编队第二轮定位时部分无人机始终无法定位

原因:这些无人机恰好位于已定位无人机构成的三角形边的延长线上,接收到的方位角只有一个有效方向信息。

解决:识别出这些无人机后跳过,在第三轮选取不共线的三架无人机作为发射机重新定位。资源里选取2、6、8作为第三轮发射机,就是避开了延长线问题。

4.5 现象:代码跑出来的坐标与理想位置表对不上

原因:理想位置表的坐标系定义与代码里的坐标系不一致。资源里圆形编队的理想位置表是以圆心为原点、0号在x轴正方向定义的,如果代码里把0号放在y轴或偏移了原点,坐标就会整体旋转或平移。

解决:在代码开头统一坐标系定义,把理想位置表转成与代码一致的直角坐标存成数组,每次比较前先对齐原点。

5. 进阶用法:把定位模型改成通用函数并做批量验证

这套代码原始版本是脚本形式,输入靠input交互,跑一次只能算一组数据。如果想做批量验证或参数扫描,得把它改成函数形式。我一般会做三件事:把定位模型封装成function,输入发射机编号、接收机编号、三个角度、半径,输出坐标;把枚举表生成部分独立成函数,只算一次存成mat文件;把位置调整的迭代逻辑写成带最大迭代次数保护的while循环。

function [x, y] = locate_drone(FY0i, FY0j, FY0k, alpha1, alpha2, alpha3, r) % 输入:两架发射机编号FY0i、FY0j,接收机编号FY0k % 三个方位角alpha1、alpha2、alpha3,编队半径r % 输出:接收机的直角坐标(x, y) 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, ~, ~] = solve(eq1, eq2, eq3, k0, a2, a3); K0 = double(X); k0i = acosd(1 - 2*cosd(alpha1)^2); k01 = k0i - 40*abs(FY0k - FY0i); x = real(K0*cosd(k01)); y = real(K0*sind(k01)); end

封装成函数后,批量验证就方便了。比如要验证所有56种组合的定位精度,可以写两层循环遍历发射机和接收机编号,每组算一次坐标,与理想位置表对比,统计误差分布。我一般会跑三组半径(10、50、100米),看误差是否随半径线性放大——如果放大倍数偏离线性,说明某个环节的数值精度有问题。

验证项半径10米误差半径50米误差半径100米误差
定位坐标偏差<0.01米<0.05米<0.1米
迭代收敛步数约200步约500步约1000步
枚举表匹配率100%100%100%

表格里的数据是我跑批量验证时的典型值,实际会因初始偏差和角度噪声略有浮动。如果100米半径下误差超过0.5米,先检查solve的数值精度,把syms换成vpasolve并设精度为10位有效数字。

从那以后我每次拿到这类几何定位代码,都强制先跑一组已知答案的验证用例,确认坐标反解和枚举匹配都对得上,再改参数做正式计算。希望帮到你。

本文还有配套的精品资源,点击获取

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进