1. 赛题回顾与核心难点拆解2022年的国赛B题题目是“无人机遂行编队飞行中的纯方位无源定位”。这个题目一出来当时很多队伍都懵了。它不像传统的优化题或者数据题那样有明确的套路可循而是把物理建模、几何计算、非线性优化和算法设计拧在了一起对参赛者的综合能力提出了很高的要求。简单来说题目描述了一个场景若干架无人机FY00-FY09组成一个圆形编队其中FY00是“主机”它只知道其他无人机相对于它的方位角也就是方向而不知道距离。我们需要做的就是仅凭这些方位角信息去推算出其他所有无人机在平面上的精确位置。这听起来有点像“只闻其声不见其人”然后要在地图上把说话的人全找出来。核心的难点非常突出这是一个典型的非线性、欠定方程组求解问题。每个方位角信息只能提供一个方向约束而我们需要求解的是每个无人机的二维坐标x, y也就是两个未知数。理论上仅凭一个方位角方程是无法确定一个点的位置的这就像你知道一个人在正东方向但不知道是100米外还是1000米外。题目巧妙的地方在于它给出了多架无人机FY01-FY09相对于FY00的方位角并且FY00自身的位置是已知的通常设为原点。此外题目还隐含了一个极强的约束所有无人机要均匀分布在一个圆周上。这个几何约束是解决整个问题的“钥匙”。很多队伍一开始就卡在了建模上不知道如何将“方位角”和“圆周均匀分布”这两个条件用数学语言表达出来。有的队伍试图用纯几何作图法去凑但无法处理有误差的情况有的队伍直接套用最小二乘法但模型建得不对导致优化结果完全偏离。所以理解并构建出正确的数学模型是攻克这道题的第一道也是最重要的一道坎。2. 核心数学模型构建从方位角到坐标约束要定位首先得把题目描述翻译成数学方程。这是整个解题思路的基石一步错步步错。2.1 方位角观测方程的建立设主机FY00位于坐标原点O(0, 0)。对于任意一架僚机FYi其坐标为(xi, yi)。题目中给出的方位角αi是指以FY00为原点从正北方向通常对应y轴正方向顺时针旋转到视线FY00指向FYi的向量所经过的角度。这里有一个关键的坐标轴定义需要统一。在数学建模中我们通常使用平面直角坐标系x轴向右y轴向上。而方位角“正北顺时针”是大地测量学中的习惯。为了简化计算一种常见且有效的处理方式是将“正北”视为y轴正方向“正东”视为x轴正方向那么方位角α就是从y轴正方向顺时针旋转到向量(xi, yi)的角度。因此向量(xi, yi)的方向与方位角αi的关系可以通过正切函数来描述tan(αi) xi / yi但这里要非常小心因为tan函数周期是π无法区分第一象限和第三象限或第二和第四象限。直接使用tan(αi) xi / yi会丢失象限信息导致方程多解。更稳健的方法是使用正弦和余弦函数sin(αi) xi / dicos(αi) yi / di其中di sqrt(xi^2 yi^2)是FY00到FYi的距离。将距离di代入我们可以得到两个等效的约束方程xi di * sin(αi)yi di * di * cos(αi)或者更常见的是将其写为比例形式消去未知的距离dixi / yi tan(αi)需结合象限判断 或者使用向量正交的概念向量(xi, yi)与方向向量(sin(αi), -cos(αi))是平行的不这里更容易出错。一个更不易出错的方法是使用点积或叉积来构建垂直关系。最推荐且不易出错的方法是利用方向向量的垂直关系。与方位角αi方向垂直的一个向量是(cos(αi), -sin(αi))你可以画图验证将方向向量顺时针旋转90度。那么位置向量(xi, yi)在方向向量(sin(αi), cos(αi))上的投影长度就是距离di而它与垂直向量的点积应为0。这引出了线性约束xi * cos(αi) - yi * sin(αi) 0这个方程完美地表达了“点(xi, yi)位于方位角为αi的射线上”这一几何事实且它是一个关于xi和yi的线性方程这大大简化了问题。我们记为观测方程1。2.2 圆形编队约束的数学表达题目要求所有无人机包括FY00均匀分布在一个半径为R的圆周上。FY00位于圆心。这是一个非常强的约束。 对于僚机FYi其坐标必须满足xi^2 yi^2 R^2我们记为约束方程2。R是未知的半径。此外“均匀分布”意味着相邻无人机之间的圆心角是相等的。对于9架僚机它们应该把360度的圆周均匀分成9份即每份40度。如果给FY01-FY09编号并确定一个起始角那么它们的理论方位角应该是已知的例如FY01在40度方向FY02在80度方向……。但题目给出的方位角αi是带有观测误差的实测值并不严格等于这些理论值。我们不能直接用理论值去反推坐标而必须将理论分布作为一个优化目标或约束。设理论方位角为θi(例如θi (i * 40)°i1,2,...,9)。那么“均匀分布”可以理解为所有无人机的实际位置应该尽可能使它们与圆心的连线方向接近这些理论角度。这可以转化为另一个优化目标最小化实际方位角与理论方位角的偏差。 实际方位角可以通过坐标计算βi atan2(xi, yi)注意atan2函数能返回正确的象限角需要转换为题目定义的“正北顺时针”角度制。 那么偏差为|βi - θi|需要考虑角度循环如355度和5度的偏差是10度不是350度。所以“均匀分布”约束不是一个硬性的等式约束而是一个软性的优化目标与观测方程一起引导优化算法找到既符合观测数据又符合编队几何形状的解。2.3 综合数学模型一个带约束的优化问题将上述分析结合起来对于每一架僚机FYi我们有以下信息来自FY00的观测数据一个带有误差的方位角测量值αi_meas。来自编队要求的几何约束它应在半径为R的圆周上且其理论方位角为θi。我们需要求解所有(xi, yi)和R。一个自然而然的建模思路是建立一个非线性最小二乘优化模型优化变量X [x1, y1, x2, y2, ..., x9, y9, R]共19个变量。目标函数由两部分构成旨在同时满足观测和编队形状。观测拟合项最小化实际位置与观测方位角之间的不一致性。利用2.1中推导的线性方程我们可以计算残差f1_i xi * cos(αi_meas) - yi * sin(αi_meas)这一项的理想值是0。将其平方和最小化意味着让所有无人机的位置尽可能满足FY00的观测。F1 Σ (f1_i)^2(i1 to 9)编队形状项最小化两个偏差。 a)半径一致性项所有无人机到原点的距离应尽可能等于R。f2_i sqrt(xi^2 yi^2) - RF2 Σ (f2_i)^2(i1 to 9) b)分布均匀性项实际方位角与理论方位角的偏差尽可能小。 计算实际方位角βi atan2(xi, yi)并转换为与αi_meas相同的度量如北东坐标系下的角度。 计算角度偏差δi min(|βi - θi|, 360 - |βi - θi|)这是处理角度循环差的标准方法。F3 Σ (δi)^2(i1 to 9)总目标函数Minimize F w1 * F1 w2 * F2 w3 * F3其中w1, w2, w3是权重系数用于平衡三项的重要性。通常观测数据项F1的权重应最大因为它是直接的测量依据。F2和F3是正则化项用于引入先验知识编队是圆形且均匀的在观测信息不足欠定时引导求解。可能的约束可以添加R 0作为边界约束。这个模型清晰地描述了问题在观测数据不准的情况下寻找一组最可能的位置使得它们既“看起来”符合FY00看到的方位又“看起来”像一个均匀的圆形编队。3. 求解算法选择与实现细节模型建好了怎么解这是一个多变量非线性优化问题。直接求解析解几乎不可能必须依赖数值优化算法。3.1 算法选型为什么是LM算法可供选择的算法很多如梯度下降法、牛顿法、高斯-牛顿法、Levenberg-Marquardt (LM) 算法等。对于本题Levenberg-Marquardt (LM) 算法是最佳选择之一。原因如下专门针对非线性最小二乘我们的目标函数是平方和形式F Σ ri(x)^2LM算法正是为此类问题设计的效率很高。鲁棒性强LM算法是高斯-牛顿法的改进版通过引入一个阻尼因子在参数更新时能在最速下降法和高斯-牛顿法之间自适应切换。当初始猜测离最优解很远时它更像梯度下降法保证收敛当接近最优解时它切换到高斯-牛顿法能快速收敛。这对于我们初始值可能给得不准的情况非常友好。成熟库支持在MATLAB、Python (SciPy) 中都有非常成熟稳定的LM算法实现如MATLAB的lsqnonlin SciPy的least_squares我们不需要自己从头编写复杂的优化代码可以专注于模型本身。注意使用这些库函数时我们需要提供残差函数而不是目标函数。即我们需要构建一个向量函数r(X) [r1, r2, ..., rM]使得F r(X)^T * r(X)。我们的f1_i,f2_i,δi就是这些残差分量。3.2 关键实现步骤与代码框架以Python SciPy为例下面勾勒一个具体的实现流程步骤1定义理论角度和观测数据import numpy as np from scipy.optimize import least_squares # 理论均匀分布角度 (假设FY01从正东方向开始逆时针排列。需根据题目图示调整) # 注意这里需要与题目中方位角定义正北顺时针进行转换。假设转换后理论角为 theta_i theta_theoretical np.deg2rad(np.array([40, 80, 120, 160, 200, 240, 280, 320, 360])) # 示例 # 实测方位角 (来自题目数据)转换为弧度制 alpha_measured np.deg2rad(np.array([...])) # 填入FY01-FY09的实测角度步骤2构建残差函数这是最核心的部分。def residuals(vars): vars: 优化变量数组 [x1, y1, x2, y2, ..., x9, y9, R] 返回残差向量 n 9 x vars[0:2*n:2] # x坐标 y vars[1:2*n:2] # y坐标 R vars[-1] # 半径 res [] # 1. 观测方程残差 (线性约束) for i in range(n): # 使用线性约束: x*cos(alpha) - y*sin(alpha) 0 res.append(x[i] * np.cos(alpha_measured[i]) - y[i] * np.sin(alpha_measured[i])) # 2. 半径约束残差 for i in range(n): dist np.sqrt(x[i]**2 y[i]**2) res.append(dist - R) # 3. 均匀分布残差 (角度偏差) for i in range(n): # 计算实际点的方位角 (注意atan2参数顺序通常atan2(y, x)返回与x轴夹角需转换) # 假设坐标系x轴东y轴北。则方位角β atan2(x, y)结果在[-pi, pi]需转换到[0, 2pi) beta np.arctan2(x[i], y[i]) # 根据你的坐标轴定义调整 if beta 0: beta 2 * np.pi # 计算与理论角度的最小差值处理360度循环 delta_angle beta - theta_theoretical[i] delta_angle np.mod(delta_angle np.pi, 2*np.pi) - np.pi # 将差值映射到[-pi, pi] res.append(delta_angle) # 残差就是角度差弧度 # 4. 可选权重可以通过在返回前对res数组不同部分乘以权重系数来实现 # weight_obs, weight_r, weight_angle 1.0, 0.5, 0.2 # res[:n] [val * weight_obs for val in res[:n]] # res[n:2*n] [val * weight_r for val in res[n:2*n]] # res[2*n:] [val * weight_angle for val in res[2*n:]] return np.array(res)步骤3提供初始猜测并求解初始值对非线性优化至关重要。一个合理的初始猜测能极大提高收敛速度和成功率。# 初始猜测假设半径R0根据理论角度和半径R0生成初始点 R_guess 100.0 # 例如100米 x0 [] y0 [] for theta in theta_theoretical: x0.append(R_guess * np.sin(theta)) # 注意sin/cos取决于你的坐标轴与角度定义关系 y0.append(R_guess * np.cos(theta)) # 将初始点列表和初始半径猜测拼接成优化变量数组 initial_guess [] for i in range(9): initial_guess.extend([x0[i], y0[i]]) initial_guess.append(R_guess) # 设置边界例如半径必须为正 bounds ([-np.inf]*len(initial_guess), [np.inf]*len(initial_guess)) # 无边界 bounds[0][-1] 0.1 # 半径R的下界为0.1 # 调用LM算法求解 result least_squares(residuals, initial_guess, boundsbounds, methodtrf, ftol1e-10, xtol1e-10, max_nfev2000) # methodtrf (Trust Region Reflective) 是SciPy中处理边界问题的稳健算法内部包含类似LM的机制。 # 提取结果 optimized_vars result.x optimized_R optimized_vars[-1] optimized_positions optimized_vars[:-1].reshape(9, 2) # 前18个是坐标步骤4结果分析与可视化求解后一定要验证结果。计算目标函数终值result.cost。绘制优化前后的位置对比图。将FY00画在原点用圆圈表示理论圆用星号表示初始猜测点用圆点表示优化后的点并用连线表示FY00到各点的视线方向应与实测方位角大致一致。计算优化后各点的实际方位角、到原点的距离与理论值、观测值对比评估定位精度。3.3 权重系数调整的经验权重w1, w2, w3的选择不是一成不变的它体现了你对不同约束的置信度。高观测权重 (w1)如果你相信FY00的方位角测量非常精确那么应赋予F1更高的权重让解更贴合观测数据即使这可能导致编队形状略有失真。高形状权重 (w2,w3)如果你认为方位角观测噪声很大而“圆形均匀编队”这个先验知识非常可靠那么可以增大w2和w3的权重。这相当于用强几何先验去“纠正”观测数据中的误差。通常的起手式可以先将所有权重设为1进行求解。观察残差项result.fun中哪一部分的残差明显偏大。如果观测残差很大说明要么模型有误要么观测数据噪声太大可能需要调整模型或降低w1。如果半径残差或角度残差很大说明优化结果严重偏离圆形可能需要检查初始值或增加形状约束的权重。在实际比赛中可以通过交叉验证的思路来调整用部分数据如7架无人机建模预测剩余2架的位置看预测误差。调整权重使预测误差最小。4. 可能遇到的陷阱与进阶讨论即使按照上述思路在实际编程求解中依然会碰到不少坑。4.1 角度循环与象限处理这是最容易出错的地方。主要体现在两个环节实测方位角输入题目给出的αi是0-360度的值。在代入cos(αi), sin(αi)计算时必须确保αi的单位是弧度并且三角函数计算正确。sin(α)和cos(α)本身不涉及循环问题。实际方位角计算在计算βi atan2(yi, xi)注意参数顺序时atan2返回的是(-π, π]范围内的弧度值。你需要将其转换到与你的理论角度θi和实测角度αi_meas相同的坐标系和范围通常是[0, 2π)。转换时要一致。角度偏差计算计算|βi - θi|时必须考虑360度的循环性。355度和5度的差是10度不是350度。正确的做法是delta np.abs(np.mod(βi - θi np.pi, 2*np.pi) - np.pi)这个公式能将任意角度差映射到[0, π]区间。踩坑实录我们队最初直接用abs(beta - theta)计算角度残差结果在角度接近0度和360度的边界处优化算法出现了不连续的剧烈跳动导致无法收敛。加上循环处理公式后问题立刻解决。4.2 初始值敏感性与多解问题非线性优化问题对初始值敏感且可能存在局部最优解。我们的问题由于观测信息不足只有方向理论上存在镜像解的可能。例如所有无人机可能位于一个半径相同的圆上但排列顺序可能是顺时针也可能是逆时针或者整体旋转一个角度。如何应对利用先验知识给出好初值题目中FY00在圆心其他无人机均匀分布在圆周上。这就是最强的先验。我们的初始猜测就基于此假设一个半径如100按理论角度放置。这通常能引导算法找到正确的全局最优解即符合我们认知的那个解。多起点优化如果担心陷入局部最优可以随机生成多组初始点在合理范围内如半径在50-200之间随机角度在理论值附近小幅随机扰动分别进行优化然后选择目标函数值最小的解作为最终结果。分析解的物理合理性优化结束后检查结果。所有无人机是否大致均匀分布到FY00的距离是否大致相等如果出现某个点特别近或特别远或者角度堆积在一起那很可能陷入了不好的局部解。4.3 模型扩展与变体思考原题是最基础的定位。但数学建模竞赛往往鼓励深入思考。这里有几个可以深化的方向用于提升论文的深度和亮点考虑观测误差模型题目说方位角有“微小误差”。这个误差是什么分布高斯白噪声有没有系统误差可以在目标函数中引入加权最小二乘假设不同方向的观测精度不同或者使用更鲁棒的损失函数如Huber损失代替平方损失以抑制可能存在的粗差outliers。时间序列与滤波如果题目给的是多个时刻的方位角序列动态定位那么问题就变成了一个跟踪问题。我们可以建立无人机的运动模型如匀速圆周运动然后使用卡尔曼滤波EKF, UKF或批处理优化滑动窗口优化来估计无人机的位置和速度。这能将不同时刻的观测信息融合得到更平滑、更准确的轨迹。这是从静态定位到动态定位的巨大飞跃。编队形状的鲁棒性不一定非要严格均匀分布。可以将其作为一个软约束并探讨当一两架无人机偏离理论位置时如故障定位算法是否依然稳健。可以引入l1范数正则化使得算法对个别异常点不敏感。三维空间扩展如果无人机不在同一高度呢方位角就包含了俯仰角。问题就变成了三维空间中的纯方位定位需要更多的几何约束或更多的观测站多架主机无人机。这可以作为一个重要的模型推广讨论。4.4 论文写作中的呈现技巧思路清晰算法有效最后要靠论文把故事讲好。模型部分一定要画出清晰的示意图标明FY00、FYi、方位角αi、半径R、理论角θi。图示能极大帮助评委理解你的建模思想。算法部分给出清晰的算法流程图伪代码说明LM算法如何与你的残差函数结合。不要只写“我们使用了lsqnonlin函数”。结果部分必须要有可视化对比图。一张图显示初始猜测、优化结果、理论圆和观测射线。另一张图可以显示残差收敛曲线证明算法有效。表格列出优化前后的坐标、距离、角度偏差对比。灵敏度分析这是拿高分的关键。分析权重系数变化对结果的影响。分析方位角测量误差增大时定位精度如何下降。分析当编队半径未知或变化时模型的适应性。这些分析能体现你对模型性能的深刻理解。模型评价与推广客观评价模型的优缺点例如对初始值敏感但加入形状约束后鲁棒性增强。简要讨论前述的扩展方向动态、三维等展示你的思考深度。最后想说的是这道题考察的远不止是编程或数学它考察的是将模糊的实际问题转化为精确数学模型的能力以及利用先验知识编队形状弥补观测信息不足的思维。很多队伍败在第一步——模型建得似是而非。一旦抓住了“线性观测方程”和“非线性几何约束”这个核心矛盾并用优化框架将其统一起来剩下的就是技术实现和细节打磨了。我们在实际求解时花了大量时间在角度坐标系的转换和残差函数的调试上一个符号错误就可能导致完全错误的结果。所以耐心、细致的验证和可视化是确保你不走弯路的最重要保障。