企业官网建设流程全解析

简介：围绕2022年全国大学生数学建模竞赛本科组B题“无人机纯方位无源定位”，这份河北省一等奖论文PDF面向参加数模竞赛的本科生及指导教师，系统讨论圆形与锥形编队下的定位建模与位置调整。内容从三机发射信号、被动接收无人机定位切入，借助二维直角坐标系、余弦定理与“三圆共点”分析轨迹交点，并讨论额外增加发射信号无人机的必要性；位置调整部分给出以0.01米为步长的遍历思路，锥形编队则扩展至更复杂的相对位置约束。文中还说明MATLAB在符号推导、数值求解和轨迹可视化中的用法，便于读者理解模型落地过程。资源包仅含1个PDF文件，大小约2.63MB，结构紧凑，适合集中研读。已有8041人学习，可作为赛题复盘、模型对照与论文写作参考。

1. 无人机纯方位无源定位这道题，真正难的是把概念对上号

无人机纯方位无源定位，第一次看会觉得几何味太重：圆形编队十架机、锥形编队十五架机、三组方位角、余弦定理，像是把高中解三角形搬进了竞赛。但真正动手的人会发现，卡住人的不是推公式，而是把一堆角度假信息和编号、坐标、迭代步长对应起来——0 号机固定在圆心，九架机均匀分布在半径 100 m 的圆周上，每架机只能收到几个夹角，剩下的全靠算。这份 2022 年数学建模竞赛本科组 B 题的完整解法加四份 MATLAB 程序（定位模型、额外机数量确定、圆形调整、锥形调整），解决的正是这条从建模到枚举再到逐机调整的链路。适合两类人：正在准备竞赛、需要一份能跑通的参考实现的同学，以及想把「定角对定边」这类几何约束真正写进代码的工程师。下面按「模型怎么立、枚举怎么做、调整怎么迭代、哪里会出错」的顺序拆。

2. 方向角转坐标：三组余弦定理联立与定角轨迹的解析式

这一章的定位很明确：把题目给出的「方位角」翻译成可求解的方程。很多人一上来就想套最小二乘，其实这题的第一步是纯几何——三架位置已知的发射机、一架位置未知的接收机，只要接收机能报出三组夹角，位置就是唯一确定的。先把这层关系立住，后面所有代码才有意义。

2.1 坐标系怎么摆：0 号机固定在原点的理由

圆形编队里，0 号机恒在圆心，1～9 号机沿圆周均匀分布，相邻两架间隔 40°。题面没说 1 号机在哪个方位，这时候有个自由度可以利用：直接把 1 号机放在 x 轴正半轴上，0 号机放原点。这么摆的好处是后面所有点的极角都变成(n-1)*40°，不用额外引入全局旋转角，仿真复现时也少一个待定参数。

代价是真实飞行时 1 号机未必真的在正东方向。所以建模阶段算出来的坐标是「编队自身坐标系」下的相对位置，实际输出给飞控时要再乘一个旋转矩阵回全局系。这一步在论文里没写，但只要是落地就躲不开，我在代码里习惯留一个R_global变量占位。

坐标生成没什么技巧，一行循环就够：

%% 圆形编队准确坐标生成（编号 0~9，下标 1~10） r = 100; % 圆周半径，题目给定 100 m theta = 0:40:320; % 九架圆周机方位角，相邻间隔 40° P = [0 0; r*cosd(theta)' r*sind(theta)']; % 10×2，第 n+1 行是 n 号机 for n = 0:9 fprintf('FY%02d (%9.3f, %9.3f)\n', n, P(n+1,1), P(n+1,2)); end

这里最关键的是索引映射：MATLAB 数组从 1 开始，题目编号从 0 开始，P(n+1,:)就是 n 号机。别看这一行简单，我见过不止一个人把P(i,:)当成 i 号机，结果整套角度表整体错位一位，最后所有匹配都对不上还查不出原因。后面第 5 章会专门讲这个坑。

参数说明：r是编队半径，题目里 100 m，做单位无量纲验证时可以设成 1，看角度是否不变；theta的顺序决定了编号顺序，一旦你把它改成40:-40:-320，编号就从顺时针变成逆时针，整张角度表全反。

2.2 三个方位角反解位置：三组余弦定理联立

设接收机为 k，三架发射机分别是 0 号（圆心）、i 号、j 号。接收机测到的三个夹角为 α₁、α₂、α₃，分别对应三条连线两两之间的张角。记 d₀ = |k−O|、dᵢ = |k−i|、dⱼ = |k−j|，那么在三个三角形 O-k-i、O-k-j、i-k-j 里各写一次余弦定理：

function [d0, di, dj] = solve_range(alpha, i, j, r) % alpha = [a1 a2 a3]，单位：度 % i、j 为两架圆周发射机的编号（1~9），0 号机自动作为第三架 syms d0 di dj positive aij = 2*r*sind(20*abs(i-j)); % 圆周上 i、j 两机的弦长，圆心角 |i-j|*40° eq1 = (d0^2 + di^2 - r^2 )/(2*d0*di) == cosd(alpha(1)); % 三角形 O-k-i eq2 = (d0^2 + dj^2 - r^2 )/(2*d0*dj) == cosd(alpha(2)); % 三角形 O-k-j eq3 = (di^2 + dj^2 - aij^2)/(2*di*dj) == cosd(alpha(3)); % 三角形 i-k-j sol = vpasolve([eq1; eq2; eq3], [d0; di; dj], [r, 1.5*r, 1.5*r]); d0 = double(sol.d0); di = double(sol.di); dj = double(sol.dj); end

三个未知数、三个方程，方程数刚好够。逻辑上，这三条边一旦确定，接收机到三个已知点的距离就全定了，位置就是三个圆的交点——这就是原题里说的「三圆共点」。而余弦定理只是把「夹角已知」翻译成「距离可求」的工具。

参数说明几个容易忽略的地方。aij里的20是圆心角|i-j|*40的一半，因为弦长公式是 2r·sin(θ/2)。vpasolve的第三个参数是初值向量，全给r会让某些钝角情况收敛到不合理解；我一般按几何直觉给[r, 1.5r, 1.5r]，接收机到圆心距离不会超过2r，到圆周机距离一般也在r附近。

得到三条边之后，还要把距离转成坐标。这里用一次正弦定理比再解一遍方程省事：在三角形 O-i-k 里，∠Oik可以用asin(di·sinα₁ / r)求出来，然后用180 − α₁ − ∠Oik得到∠iOk，接收机的极角就是θᵢ − ∠iOk。角度正负号这一步要对着图看，接收机落在 i 号机内侧还是外侧，得用第二组角度判一次方向。

2.3 定角轨迹：钝角画出「眼睛」，锐角画出「葫芦」

论文里最出彩的一段，是把「一个动点到两定点的夹角恒定」这个约束解成了显式轨迹。结论是：若两定点间距离为 L、张角为 α，那么动点轨迹是两段关于两点连线对称的圆弧，圆心在连线的中垂线上，偏离连线中点的距离 h = (L/2)/tanα，圆弧半径 R = L/(2·sinα)。

当 α 是钝角时，tanα 为负，两个圆心靠近连线内侧，两段弧在弦的两端收得很紧，画出来像一只「眼睛」；当 α 是锐角时，两个圆心退到连线外侧较远的地方，两段弧向外鼓出，形状像个「葫芦」。这个形状差异不是美术问题——它直接决定了「两架发射机够不够」这个结论。

function draw_locus(A, B, alpha_deg) % A、B 为两架发射机坐标，alpha_deg 为接收机处的张角 L = norm(B - A); R = L / (2*sind(alpha_deg)); % 轨迹圆半径 h = (L/2) / tand(alpha_deg); % 圆心到弦中点的有向距离 M = (A + B) / 2; nv = [-(B(2)-A(2)), B(1)-A(1)] / L; % 弦的单位法向量 C1 = M + h*nv; C2 = M - h*nv; % 上下两个圆心 tA = atan2(A(2)-C1(2), A(1)-C1(1)); tB = atan2(B(2)-C1(2), B(1)-C1(1)); t = linspace(tA, tB, 200); % 只画对应的一段弧 plot(C1(1)+R*cos(t), C1(2)+R*sin(t), 'b'); hold on; end

这段代码只画了上半个圆心对应的弧，另外半边把C1换成C2再调用一次。实际建模时它不是用来算数的，而是用来做几何判断和画图汇报的——论文里那几张「眼睛」「葫芦」的插图，本质就是这个函数跑出来的。

真正要记住的结论是：一条轨迹是一条封闭曲线，一条曲线跟另一条曲线最多交两点，所以两组方位角只能把位置压到两个对称候选，还差一个约束。这就是为什么 2 架发射机不够、3 架刚好。第 3 章会用枚举把这个结论压实。

3. 额外一架就够：56 组枚举表与平方和最小匹配

问题一第二小问问的是「除 0 号和 1 号外，还需要几架发射信号」。思路听起来是组合数学，其实是上一章的几何结论加一次穷举验证。这一章把「为什么是 1 架」讲清楚，再把那张角度对照表用代码复现出来。

3.1 为什么两架发射机不够、三架刚好

0 号固定在圆心、1 号固定在圆周，这两架发射时，任意一架接收机只能得到「0 号与 1 号之间的张角」一个信息。按 2.3 的结论，这一个角对应一条封闭的定角轨迹，接收机的真实位置在这条曲线上，但具体在哪一点完全不知道，一条曲线上的点有无穷多个。

再加一架 m 号发射机，接收机就能多拿到「0 号与 m 之间」和「1 号与 m 之间」两个角。这两条新轨迹和旧轨迹一起，最多交出唯一一个公共交点——也就是原题里的「三圆共点」。反过来看，如果只加 0 架，也就是只有 0 号和 1 号发射，接收机永远只有一条轨迹，位置不可解。

这里有个容易被质疑的点：三架发射机是否一定唯一？理论上如果三架机共线、或者接收机恰好落在某条边的延长线上，退化成两条平行轨迹，交点仍然不唯一。所以「一架就够」的完整表述应该是「一架就够，且在一般情况下唯一」。第 5 章会把退化情况单独拿出来讲。

顺带说一句，论文里还专门证了「定角轨迹为什么是圆弧」：设两个定点为 (−a,0) 和 (a,0)，张角恒为 α，把 tan 展开整理后能得到一个标准圆方程 x² + (y ± h)² = R²，上下对称两条弧拼成那个封闭曲线。这个推导的代数细节不重要，重要的是它给了你一个可以直接调用的解析式，而不是让你去解隐函数。

3.2 56 组组合的枚举：把角度对照表用代码生成

确定了「0 号、1 号 + 再加一架 m」这个结构后，问题变成：m 可以是 2～9 号中的哪一架？接收机 k 可以是 1～9 号中除 m 外的哪一架？一共 8×7 = 56 种有效组合，论文里那张角度对照表就是把每种组合下的三个角都算出来。

手工算 56 组显然不现实，用代码就三行核心逻辑：

r = 100; theta = 0:40:320; P = [0 0; r*cosd(theta)' r*sind(theta)']; % 第 n+1 行 = n 号机 % 计算：站在 k 号机位置，看向 i 号和 j 号机所成的夹角 angle_at = @(k, i, j) acosd( dot(P(i+1,:)-P(k+1,:), P(j+1,:)-P(k+1,:)) / ... (norm(P(i+1,:)-P(k+1,:)) * norm(P(j+1,:)-P(k+1,:))) ); Tab = nan(10, 10, 3); % Tab(m, k, :) = 0、1、m 发射时，k 收到的三个角 for k = 2:9 for m = 2:9 if m == k, continue; end Tab(m, k, 1) = angle_at(k, 0, 1); % 0 号与 1 号的张角 Tab(m, k, 2) = angle_at(k, 0, m); % 0 号与 m 号的张角 Tab(m, k, 3) = angle_at(k, 1, m); % 1 号与 m 号的张角 end end

angle_at用的是向量点积公式，acosd直接返回角度制，避免弧度转换漏乘 180/π。这个匿名函数是整个流程里复用率最高的东西，后面锥形编队还要用同一套逻辑，只是把坐标矩阵换成锥形网格。

Tab是一个三维数组，第一维是增派的发射机 m，第二维是接收机 k，第三维存三个角。为什么这么摆？因为实际使用时的查询方向是「已知接收机 k、已知三个角度，反查 m」，正好是沿着第一维扫。如果摆反了，后面每次查询都要转置一次，跑批量数据时性能会差一个量级。

参数说明：k从 2 开始是因为 1 号机是固定发射机，不需要给它自己算接收角；if m == k那句跳过的是「接收机不能同时当发射机」的物理约束，漏掉这一句会在表里留下自反的无效组合。

3.3 平方和最小匹配：从实测角反查该增派哪一架

真实场景里，接收机报上来的角度是有偏差的，不可能是表格里的精确值。所以匹配不能用等号判断，得用距离度量。论文用的是三个角偏差的平方和：

beta = [b1, b2, b3]; % 接收机实测的三个方位角（度） best = inf; sel = -1; for m = 2:9 if m == k, continue; end ref = squeeze(Tab(m, k, :))'; % 该组合下的三个理论角 s = sum((beta - ref).^2); % 平方和代价 if s < best best = s; sel = m; end end fprintf('接收机 FY%02d 需要增派 FY%02d，代价 %.4f\n', k, sel, best);

逻辑说白了就是最近邻搜索：把「实测三角度」当成一个三维点，56 个理论组合就是 56 个参考点，找欧氏距离最小的那个。用平方和而不用绝对值和，是为了让大偏差被放大，避免某一项偏得离谱、另外两项偏得少就被平均过去。

参数上要注意两点。一是squeeze的使用，三维数组取单个(m,k)切片出来是 1×1×3，不 squeeze 的话beta - ref会做隐式扩展，得到 1×1×3 而不是 1×3，sum的结果数量级就不对了。二是当k已经等于当前遍历的m时必须跳过，否则会算到「接收机给自己发信号」的组合，代价虽然算得出来但物理上无意义。

还有一个隐含假设：不同组合下的理论角度差异要足够大，匹配才有意义。相邻两架机之间隔 40°，反映到张角上通常在几度到十几度的量级，只要测量误差在 1° 以内，最近邻基本不会认错。但如果误差到了 10° 量级，就会出现模糊匹配，这时候要么提高测角精度，要么按 2.2 的思路直接解方程而不查表。

4. 两种编队的逐机调整：0.01 米步进与三角形内角和补齐

建模和匹配解决的是「我在哪」，这一章解决的是「怎么飞过去」。圆形编队和锥形编队的调整逻辑用的是同一个内核，区别在于发射机怎么选、退化情况怎么绕开。整个流程是串行的：定位一架、调整一架、把它加入「位置无偏差集合」、再定位下一架。

4.1 圆形编队：0.01 米步长的极坐标遍历

题目给的条件是 0 号和 1 号初始位置已经在准确位置上，并且假设飞行过程中两者相对关系不变。这意味着每轮迭代里，这两架机始终可以作为无偏差的发射源。

对 2 号机来说，它能收到 0 号和 1 号两个方位信息，但只有两个角，位置是 2.3 里那条封闭轨迹上的任一点——所以它并不知道自己该往哪边挪。论文的处理很直接：设定步长 0.01 m，以当前位置为起点，向四面八方「试走」一小步，每一步都用定位模型重新算一遍自己到准确位置的距离，取距离最小的方向走。

step = 0.01; % 单次调整步长，单位 m pos = [x0, y0]; % 当前实际位置 tgt = [xt, yt]; % 由定位模型算出的准确位置 iter = 0; while norm(pos - tgt) > step && iter < 2e5 bestPos = pos; bestD = norm(pos - tgt); for d = 0:5:355 % 方位遍历，5° 一档 cand = pos + step * [cosd(d), sind(d)]; dc = norm(cand - tgt); if dc < bestD bestD = dc; bestPos = cand; end end pos = bestPos; iter = iter + 1; end fprintf('收敛用 %d 步，残差 %.4f m\n', iter, norm(pos - tgt));

这段代码有两个细节值得说。第一，遍历用 5° 一档而不是 1° 一档，是因为步长本身只有 0.01 m，5° 的方向误差带来的横向偏差约 0.01×(1−cos5°) ≈ 4e-6 m，远小于收敛阈值，没必要把内层循环撑到 360 次。第二，iter上限必须有，我遇到过角度表错位导致模型算出的tgt本身就是错的，循环永远收敛不了，while直接跑成死循环。

参数说明：step取 0.01 m 是题目建议值，也是精度与迭代次数的折中——步长减半，迭代次数大致翻倍；tgt不是真值，是模型算出来的估计值，所以最终残差反映的是「模型自洽程度」，不是「绝对定位误差」。

这里必须诚实指出一个问题：这段遍历是离线算法。它假装接收机知道自己应该到哪，而真实场景里接收机只知道三个角度，不知道目标坐标。工程上正确的做法是把「到目标的距离」换成「实测三角度与逐渐逼近目标时的理论角度之间的残差」，然后用梯度下降或牛顿法调整。论文用遍历法是因为题目允许「调整时间不计」，而且求解精度可控。真要做实物，别照搬。

4.2 逐机锁定：定位一架、固化一架、再定位下一架

圆形编队的调整不是一次性把所有飞机都算一遍，而是串行的。第一轮：0、1 号已知无偏差，用它们去定位 2 号机，然后按 4.1 的方法把 2 号机挪到准确位置（题目假设调整时间不计，可以认为挪到即到位）。第二轮：0、1、2 号现在都无偏差，用这三架去定位 3 号机，调整到位。依此类推，直到 9 号机。

为什么必须串行？因为 0 号和 1 号只能提供一个方位信息，只有两个发射源的时候，接收机还是那条封闭轨迹上的任一点。必须等 2 号机挪到位，第三架无偏差发射机才存在，三圆共点的条件才成立。

这个思路用伪代码表达就是一个集合不断扩大的过程：

表里没有列出每轮具体用哪三架发射，是因为选择不唯一——只要保证三架不共线、且接收机不在任何一条边的延长线上，解就是唯一的。实际操作中我倾向于固定用 {0, 1, k−1} 这个组合，因为 0 号和 1 号的位置永远不变，只有第三架在变，几何退化的情况容易预先枚举排除。

有一处论文里没明说但影响很大：每一轮调整完之后，接收机的位置变成了「准确位置」，但它的角度测量误差并没有消失。所以轮次越靠后，参与定位的「无偏差」飞机越多，但每一架都只是「在容差内被认为无偏差」。这个误差是会累积的，实际操作中每调整完一架，都建议回代一次，看它与 0 号机的距离是否仍然保持在 r ± 容差 之内。

4.3 锥形编队：两发一收，用内角和补出第三个角

问题二把队形换成锥形，十五架机排成等边三角形网格，相邻间距相等。这里的麻烦在于：锥形编队不像圆形编队有天然的中心机可以当作固定参考点，必须人为假设两架机（论文里假设的是 1 号和 5 号）位置无偏差，才能启动迭代。

先构造网格。按行从上到下、每行从左到右编号：

L = 100; % 相邻两机间距 h = L * sind(60); % 行间竖直距离，等边三角形的行高 rowX = { 0, ... [-L/2, L/2], ... [-L, 0, L], ... [-1.5*L, -0.5*L, 0.5*L, 1.5*L], ... [-2*L, -L, 0, L, 2*L] }; Pc = []; for rr = 1:numel(rowX) y = -(rr-1)*h; for xx = rowX{rr} Pc = [Pc; xx, y]; % 第 n 行 = n 号机的准确位置 end end

rowX用元胞数组是因为每行的机数不同，硬编码成一个二维矩阵会带来大量 NaN 填充，后续索引反而更容易出错。h = L·sin60°是等边三角形的高，行间距和列间距不相等，这一点是最容易被忽略的——很多人直接按 L 同时作为行距和列距，画出来的网格是斜的。

锥形编队定位 2 号机的思路，本质是利用三角形内角和。1 号和 5 号已知，2 号待求，那么在三角形 1-2-5 里：先用 1、2 号发射、5 号接收，得到一个张角 A；再用 2、5 号发射、1 号接收，得到另一个张角 B；那么 2 号处的角 C = 180° − A − B。

% 锥形编队的角度计算，Pc 的第 n 行是 n 号机 ang_c = @(k, i, j) acosd( dot(Pc(i,:)-Pc(k,:), Pc(j,:)-Pc(k,:)) / ... (norm(Pc(i,:)-Pc(k,:)) * norm(Pc(j,:)-Pc(k,:))) ); A = ang_c(5, 1, 2); % 5 号机处，1 与 2 的张角 B = ang_c(1, 2, 5); % 1 号机处，2 与 5 的张角 C = 180 - A - B; % 2 号机处的角，由三角形内角和补齐 L15 = norm(Pc(1,:) - Pc(5,:)); % 已知边 L12 = L15 * sind(A) / sind(C); % 正弦定理求 1-2 边长 dir15 = (Pc(5,:) - Pc(1,:)) / L15; % 1 指向 5 的单位向量 dir12 = [dir15(1)*cosd(B) - dir15(2)*sind(B), ... dir15(1)*sind(B) + dir15(2)*cosd(B)]; % 绕 1 号旋转 B P2_est = Pc(1,:) + L12 * dir12; % 2 号机估计位置

这段代码里最值得看的是dir12那一行：用二维旋转矩阵把 1→5 的方向旋转 B 角，得到 1→2 的方向。方向和正负号取决于锥形网格的编号顺序，一旦你改了编号规则，旋转角要跟着取反。这是整个锥形部分最容易翻车的地方，比公式本身更容易出错。

两架发射机只能给出两个角，这在圆形编队里是不够的，但在锥形编队里够用，原因是三角形内角和提供了一个额外约束。这是问题二相对问题一的核心技巧，也是这道题设计得比较巧的地方。

迭代流程分三轮。第一轮用 1、5 号定位 2 号，调整到位；第二轮 1、2、5 号都无偏差，去定位 3、6、8、10、12、15 号；第三轮选 2、6、8 号作发射机，把剩下的 4、7、9、11、13、14 号定位出来。为什么第二轮定位不了那六架？因为它们恰好落在三角形 1-2-5 某条边的延长线上，张角退化成 0° 或 180°，两个角的信息量塌缩，方程病态。这是避坑章要展开的第三类问题。

5. 避坑排查：跑这套 MATLAB 最容易翻车的五个地方

这套代码我在不同机器、不同 MATLAB 版本上重启过好几轮，出错的地方高度集中。下面五条按「现象 → 原因 → 解决」写，都是实际遇到过的，不是理论假设。

5.1vpasolve返回复数或负数边长

现象：调用solve_range之后，d0或di出来是负数，甚至是带虚部的复数，后面算坐标时报错或画到奇怪的位置。

原因：余弦定理联立的方程组在几何上可能有不止一组代数解，其中包含「边长为负」的伪解。符号求解器不带物理约束，会把这些解一并返回。另外，syms d0 di dj positive只声明了变量为正，但对vpasolve的行为约束有限，某些版本会忽略这个假设。

解决：一是给初值而不是给空，vpasolve(..., [r, 1.5*r, 1.5*r])，让它从一个物理上合理的区域开始搜；二是拿到解后强制过滤，d0 = double(sol.d0); if ~isreal(d0) || d0 <= 0, error('无物理解'); end；三是更稳的做法，直接用 2.2 里的边长范围做上下界，把vpasolve换成带区间的vpasolve(eq, var, [lb ub])，把求解限制在[0.1r, 2r]之内。

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