简介：本资源为2022年全国大学生数学建模竞赛本科组B题的完整论文PDF，聚焦无人机纯方位无源定位问题，适合备战数模竞赛的本科生及指导教师参考。论文围绕圆形与锥形两种编队，系统建立被动接收信号无人机的定位模型与位置调整模型，涵盖余弦定理反解、三圆共点判定、枚举法确定发射信号机数量、步长遍历调整等核心方法，并给出问题一三个小问与问题二的完整求解过程。资源包内仅含1个PDF文件，约2.63MB，内容为获奖论文全文，含摘要、问题重述、模型假设、符号说明及模型建立与求解等标准章节，结构完整便于对照学习。该论文曾获河北省一等奖，目前已有8039人学习下载，读者可从中获取赛题分析思路、MATLAB编程实现要点及论文写作规范，适合作为数模竞赛入门与进阶的参考范例。
1. 纯方位无源定位到底难在哪：从一道竞赛题说起
2022年全国大学生数学建模竞赛本科组B题，题目背景是无人机集群在编队飞行中的纯方位无源定位问题。简单说，就是若干架无人机只靠测量彼此之间的方位角（角度信息），不依赖GPS、不主动发射信号，推算出整个编队中每架无人机的位置。这个问题的核心难点在于：角度信息是非线性的，而且存在测量误差，少量观测点要反推多维位置，天然存在多解性和病态性。
我第一次接触这类问题时，直觉是“角度交会嘛，两条射线交点就是目标”。但实际动手才发现，当无人机自身位置也有偏差、角度测量有噪声时，最小二乘的解会剧烈震荡，迭代不收敛是家常便饭。这道题真正考察的是：如何用极坐标几何建立观测方程、如何用最小二乘或优化方法求解、如何分析误差传播和发射方案的最优性。适合有线性代数、最优化基础，想把这套方法用到实际定位场景的读者。
2. 从极坐标到观测方程：把几何关系翻译成可解的数学形式
2.1 坐标系选择与变量定义
纯方位定位的第一步不是写代码，而是把“谁测谁、测什么角”定义清楚。常见做法是建立全局平面直角坐标系，设编号为 FY00 的无人机为基准（位置已知或作为参考原点），其余无人机编号 FY01 到 FY0n。每架无人机搭载一个能测量相对方位角的传感器，测得的是目标无人机相对于自身朝向的角度。
这里有个容易翻车的点：方位角是相对于正北还是相对于机体轴？题目里通常给的是相对于某条参考线的夹角。我一般会统一约定：所有角度以全局坐标系 x 轴正方向为起点，逆时针为正。如果原始数据是相对机体朝向的，必须先加上机体航向角做一次坐标旋转，否则后面方程全错。
设第 i 架无人机位置为 $(x_i, y_i)$，它观测第 j 架无人机得到的方位角为 $\theta_{ij}$，则观测方程是：
$$\tan(\theta_{ij}) = \frac{y_j - y_i}{x_j - x_i}$$
这个方程看起来简单，但它对未知量是非线性的，而且当 $x_j - x_i$ 接近零时会数值爆炸。实际处理时不会直接用 tan，而是转成：
$$(y_j - y_i)\cos\theta_{ij} - (x_j - x_i)\sin\theta_{ij} = 0$$
这样避免了除零，也把方程写成了关于未知坐标的线性形式——前提是角度已知。这个变形是整道题的关键突破口。
2.2 观测方程组的矩阵组装
假设基准机 FY00 和另外两架位置已知的无人机作为发射源，它们对未知无人机 FY0k 各测一个方位角，那么对每架未知机就有若干条上述线性方程。把所有未知坐标排成一个列向量 $\mathbf{p} = [x_1, y_1, x_2, y_2, \dots]^T$，每条观测方程贡献一行，组装成：
$$\mathbf{A}\mathbf{p} = \mathbf{b}$$
其中 $\mathbf{A}$ 的每一行由 $\sin\theta$ 和 $\cos\theta$ 组成，$\mathbf{b}$ 由已知发射源坐标和角度项组成。下面是一段组装矩阵的 Python 代码：
import numpy as np def build_observation_matrix(known_pos, unknown_ids, angles): """ known_pos: dict, {发射源编号: (x, y)} unknown_ids: list, 未知无人机编号列表 angles: dict, {(发射源编号, 目标编号): 方位角(弧度)} 返回 A, b，使得 A @ p = b，p 为未知坐标展平向量 """ n_unknown = len(unknown_ids) idx_map = {uid: i for i, uid in enumerate(unknown_ids)} rows = [] rhs = [] for (src_id, tgt_id), theta in angles.items(): if tgt_id not in idx_map: continue xs, ys = known_pos[src_id] row = np.zeros(2 * n_unknown) col = 2 * idx_map[tgt_id] # 方程: (y_j - y_i)cos - (x_j - x_i)sin = 0 # 展开后未知项系数: x_j 系数 -sin, y_j 系数 cos row[col] = -np.sin(theta) # x_j 系数 row[col + 1] = np.cos(theta) # y_j 系数 rows.append(row) # 已知项移到右边: x_i*sin - y_i*cos rhs.append(xs * np.sin(theta) - ys * np.cos(theta)) A = np.array(rows) b = np.array(rhs) return A, b逻辑说明：每条观测方程里，未知的是目标坐标 $(x_j, y_j)$，已知的是发射源坐标和角度。把未知项系数放进矩阵 $\mathbf{A}$，已知项移到等号右边放进 $\mathbf{b}$。参数上，angles的键是（发射源，目标）二元组，值必须统一为弧度制；如果原始数据是角度制，先乘np.pi/180。idx_map保证未知量在向量里的排列顺序和矩阵列一一对应，这是后面求解不出错的前提。
2.3 最小二乘求解与秩条件判断
方程组装好之后，直接调用最小二乘：
def solve_positions(A, b): p, residuals, rank, sv = np.linalg.lstsq(A, b, rcond=None) if rank < A.shape[1]: print(f"警告: 矩阵秩亏, rank={rank}, 未知量数={A.shape[1]}, 解不唯一") return p, rank这里必须检查rank。纯方位定位有一个经典结论：如果所有发射源和目标共线，或者观测几何退化，矩阵会秩亏，解不唯一。实际题目里，三架发射源如果近似在一条直线上，对远处目标的定位精度会急剧下降。我一般会额外算一下条件数np.linalg.cond(A)，超过 $10^6$ 就要警惕数值不稳定。
3. 迭代优化与误差分析：当最小二乘不够用时怎么办
3.1 为什么线性化之后还要做非线性优化
上一章的线性化有个隐含假设：角度是精确已知的。但实际测量角度有噪声，而且当我们把非线性方程强行转成线性形式时，噪声的传播方式变了——角度误差 $\delta\theta$ 在方程里被放大为与距离相关的项。当目标距离发射源较远时，同样的角度误差对应更大的位置误差，这就是纯方位定位的“距离放大效应”。
所以更稳的做法是：用线性最小二乘的解作为初值，再做一次非线性最小二乘，直接最小化角度残差：
$$\min \sum_{i,j} \left( \theta_{ij}^{\text{measured}} - \theta_{ij}^{\text{computed}}(\mathbf{p}) \right)^2$$
其中 $\theta_{ij}^{\text{computed}}$ 由当前估计位置反算：
from scipy.optimize import least_squares def angle_residuals(p_flat, known_pos, unknown_ids, angles): idx_map = {uid: i for i, uid in enumerate(unknown_ids)} res = [] for (src_id, tgt_id), theta_meas in angles.items(): if tgt_id not in idx_map: continue xs, ys = known_pos[src_id] col = 2 * idx_map[tgt_id] xj, yj = p_flat[col], p_flat[col + 1] theta_calc = np.arctan2(yj - ys, xj - xs) # 角度残差归一化到 [-pi, pi] d = theta_meas - theta_calc d = (d + np.pi) % (2 * np.pi) - np.pi res.append(d) return np.array(res) def refine_positions(p_init, known_pos, unknown_ids, angles): result = least_squares( angle_residuals, p_init, args=(known_pos, unknown_ids, angles), method='lm', max_nfev=2000 ) return result.x, result.cost, result.optimality参数说明：method='lm'适合残差数大于未知量数的场景；如果残差数接近未知量数，改用'trf'更稳。max_nfev控制最大函数评估次数，太小会提前终止，太大会浪费时间，一般 1000 到 5000 够用。optimality反映梯度范数，如果它不接近零，说明没收敛到局部极小，需要换初值重试。
3.2 误差传播的量化：蒙特卡洛比公式更实用
理论上可以用雅可比矩阵做一阶误差传播，但纯方位问题非线性强，一阶近似在角度误差较大时不准。我一般直接用蒙特卡洛：给角度加高斯噪声，重复求解 N 次，统计位置估计的均值和标准差。
def monte_carlo_error(known_pos, unknown_ids, angles_true, sigma_deg=1.0, n_trials=500): sigma = np.deg2rad(sigma_deg) errors = [] for _ in range(n_trials): angles_noisy = {k: v + np.random.normal(0, sigma) for k, v in angles_true.items()} A, b = build_observation_matrix(known_pos, unknown_ids, angles_noisy) p0, _ = solve_positions(A, b) p_refined, _, _ = refine_positions(p0, known_pos, unknown_ids, angles_noisy) errors.append(p_refined) errors = np.array(errors) return errors.mean(axis=0), errors.std(axis=0)这段代码的价值在于：它能告诉你当前发射源布局下，哪些无人机的位置估计最不稳定。如果某架机的标准差是其他的好几倍，说明它相对发射源的几何构型不好，需要调整发射方案。
3.3 发射方案优化：枚举加指标筛选
题目通常会问“最少需要几架发射源”“哪几架发射最合适”。我的做法是枚举所有可能的发射源组合，对每种组合算蒙特卡洛标准差，选标准差最小且发射源数量最少的方案。
| 发射源数量 | 可选组合数（从n架中选） | 平均定位标准差 | 是否满足精度 |
|---|---|---|---|
| 2 | C(n,2) | 较大，部分目标发散 | 否 |
| 3 | C(n,3) | 中等 | 视几何而定 |
| 4 | C(n,4) | 较小 | 通常满足 |
实际枚举时，n 一般不超过 10，组合数可控。关键是评价指标要统一：我一般用所有未知机位置标准差的均方根作为总指标，再单独看最差的那架，避免“平均好但个别发散”的方案被选中。
4. 避坑与排查：纯方位定位里最容易翻车的五个地方
4.1 角度单位混用导致结果整体偏移
现象：求解出来的位置整体旋转了一个固定角度，或者完全离谱。原因：原始数据里角度是度，代码里当弧度用了，或者arctan2返回弧度但和角度制数据直接相减。解决：在数据入口处统一转换，写一个to_radian函数，所有角度进模型前必须过这一道。我习惯在angles字典构建时就全部转成弧度，后面不再碰单位。
4.2 角度残差没有归一化导致迭代震荡
现象：非线性优化迭代过程中残差忽大忽小，不收敛。原因：两个角度相减可能得到 $350°$ 这种值，实际只差 $10°$。解决：残差统一用(d + pi) % (2*pi) - pi归一化到 $[-\pi, \pi]$，上面代码里已经处理。这个坑我在第一次做角度相关优化时踩过，调了一晚上才反应过来。
4.3 发射源共线导致矩阵秩亏
现象：lstsq返回的解数值巨大，或者rank小于未知量数。原因：所有发射源近似在一条直线上，观测几何退化。解决：组装矩阵后先检查rank和条件数，如果秩亏，换发射源组合。实际题目里，如果三架发射源夹角小于 15°，定位精度就会明显下降。
4.4 初值太差导致非线性优化收敛到局部极小
现象：线性解已经不错，但非线性优化后反而变差。原因：least_squares默认初值敏感，如果线性解本身偏差大，可能收敛到错误的局部极小。解决：用多组初值（线性解、线性解加随机扰动）分别优化，取 cost 最小的。我一般跑 5 到 10 组初值，多花几秒但稳很多。
4.5 蒙特卡洛次数太少导致标准差估计不可信
现象：重复跑两次蒙特卡洛，标准差差很多。原因：试验次数不够，统计量还没稳定。解决：先用 100 次粗跑看趋势，正式结果至少 500 次。如果标准差本身很大，说明方案确实不稳定，加次数也救不了，得改发射源布局。
5. 进阶技巧：用几何精度因子快速筛选发射方案
前面枚举加蒙特卡洛虽然可靠，但组合多的时候计算量大。实际比赛中时间有限，我一般先用几何精度因子（GDOP）做快速筛选，再对前几名做精细蒙特卡洛。
GDOP 的核心思想是：定位误差不仅取决于测量误差，还取决于观测几何。对纯方位定位，可以构造雅可比矩阵 $\mathbf{J}$，其中每一行是角度对未知坐标的偏导：
$$\frac{\partial \theta_{ij}}{\partial x_j} = \frac{-(y_j - y_i)}{r_{ij}^2}, \quad \frac{\partial \theta_{ij}}{\partial y_j} = \frac{x_j - x_i}{r_{ij}^2}$$
其中 $r_{ij}^2 = (x_j - x_i)^2 + (y_j - y_i)^2$。组装雅可比后，GDOP 定义为：
$$\text{GDOP} = \sqrt{\text{trace}((\mathbf{J}^T\mathbf{J})^{-1})}$$
GDOP 越小，几何构型越好。下面是一段计算代码：
def compute_gdop(known_pos, unknown_ids, positions, angles): """ 在给定位置估计处计算 GDOP positions: dict, {编号: (x, y)} 包含已知和未知 """ idx_map = {uid: i for i, uid in enumerate(unknown_ids)} n_unknown = len(unknown_ids) J = [] for (src_id, tgt_id), _ in angles.items(): if tgt_id not in idx_map: continue xs, ys = known_pos[src_id] xj, yj = positions[tgt_id] dx, dy = xj - xs, yj - ys r2 = dx*dx + dy*dy row = np.zeros(2 * n_unknown) col = 2 * idx_map[tgt_id] row[col] = -dy / r2 # d(theta)/d(xj) row[col + 1] = dx / r2 # d(theta)/d(yj) J.append(row) J = np.array(J) try: cov = np.linalg.inv(J.T @ J) gdop = np.sqrt(np.trace(cov)) except np.linalg.LinAlgError: gdop = np.inf return gdop参数说明：positions需要包含所有无人机的当前估计位置，已知机用真实坐标，未知机用上一步的解。r2是距离平方，不能为零，否则说明目标和发射源重合，实际不会出现但要防御性判断。np.linalg.inv可能失败，用try-except兜住，返回inf表示该方案不可用。
用 GDOP 筛选的流程是：对每个候选发射源组合，先用线性最小二乘快速解一遍得到粗略位置，再算 GDOP，按 GDOP 从小到大排序，取前 3 到 5 个组合做完整蒙特卡洛。这样能把计算量降一个数量级，同时不丢好方案。
我自己的习惯是：GDOP 只用来排序，最终决策一定看蒙特卡洛的实际标准差。因为 GDOP 是一阶近似，在角度误差大或者目标距离远的时候，排序可能和实际有出入。有一次我完全信 GDOP 选了一个组合，结果蒙特卡洛跑出来最差那架机的标准差比第二名大了 40%，从那以后就改成“GDOP 筛、蒙特卡洛定”。
最后一个实用技巧：如果题目允许无人机之间互相通信，可以把所有观测方程集中起来做全局最小二乘，而不是每架机单独解。全局解的精度通常更高，因为信息共享了。代价是矩阵变大，但 n 在几十架以内，现代计算机完全扛得住。希望帮到你。
本文还有配套的精品资源，点击获取
