2023数学建模B题多波束测线代码包：几何骨架与递推实现

简介：面向备战全国大学生数学建模竞赛、希望复盘2023年国赛B题解法的参赛者与指导教师，这份压缩包提供一套可参考的代码与结果材料，可用于理解题意拆解、模型实现和结果核对。包内共24个文件，涵盖xlsx结果表、py求解脚本、xml项目配置、png图示以及pdf/doc说明文档等，其中xlsx多用于结果汇总，py脚本承载求解过程，xml与md辅助项目配置和说明，压缩包约2.34MB，体量轻便，便于按代码、数据、文档三类快速检索。已有2084人学习下载。读者可重点查看四问对应的Python脚本、结果表格、附录数据与README说明，结合图示和文档梳理建模思路、复现计算流程并对照结果差异与改进方向，也可借鉴其数据组织与代码结构；使用时应注意遵守竞赛规则与学术规范，理解方法而非直接照搬。

1. 拆开这份 B 题代码包之前，先想清楚你要的是答案还是骨架

2023 年数学建模竞赛 B 题那道多波束测线题，题目正文读起来不长，真正让人卡住的不是建模思路，而是几何关系：覆盖宽度到底按水平投影算还是按沿坡面算，相邻测线的重叠率分母是谁，深度从换能器吃水线量还是从海平面量。这三个问题任意一个口径搞错，后面问题一到问题四的结果会整条链歪掉。这份“2023年数学建模国赛B题代码包”就是围着这套几何骨架铺开的：先封一个几何层，再往上搭测线递推、坡面分段、实测网格布线四层脚本，每一层都能单独跑、单独改参数、单独出图。适合两类人：赛前想拿真题练手、需要一份跑得通又能改的实现；或者模型已经写完，某条结果和参考值对不上，想找个对照。

2. 多波束测线的几何骨架：三条关系式决定代码能不能跑

几何层是整份代码的地基。这一层如果只是“抄个公式对上答案”，后面问题二到问题四的每一个参数都会变成玄学调参；反过来，把这三条关系式的定义域、口径和适用边界搞清楚，剩下三问本质上是同一个递推在不同坡面上的套用。我在拆包时最先看的就是几何函数那个文件，因为那里一旦有符号或角度约定问题，出图上看不出来，数值上会偏得莫名其妙。

2.1 覆盖宽度的推导：从斜距落到水平投影

设换能器位于海面附近的 O 点，向下的竖直深度为 D，海底是一个与水平面夹角 α 的斜面，换能器开角（全角）为 θ。在垂直于测线方向的剖面里建立坐标系：竖直向下为 y 轴，水平方向为 x 轴，坡面沿 x 正方向变深，则海底直线可以写成 y = D + x·tanα。

某条波束射线与竖直方向的夹角记为 φ，方向向量是 (sinφ, cosφ)，沿这条射线走 t 距离后的落点是 (t·sinφ, t·cosφ)。把它代进海底直线：

t·cosφ = D + t·sinφ·tanα t·(cosφ - sinφ·tanα) = D t = D·cosα / cos(φ + α)

这个 t 是斜距，也就是声波实际传播的距离。工程上更关心的是落点在水平方向上的位置：

x(φ) = t·sinφ = D·cosα·sinφ / cos(φ + α)

把 φ 取到两侧边界 ±θ/2，就得到左右两条边界射线的水平落点。右侧半宽 x(+θ/2)、左侧半宽 x(−θ/2)，两者相加即该条测线在水平方向上的覆盖宽度：

W_h = D·cosα·sin(θ/2)·[ 1/cos(θ/2 + α) + 1/cos(θ/2 − α) ]

沿坡面的实际覆盖宽度是 W_s = W_h / cosα。平坦海底（α = 0）时退化为 W_h = 2D·tan(θ/2)，θ = 120°、D = 100 m 对应 346.4 m，这个数值可以直接拿来校核函数写得对不对。

2.2 深度 D 的定义与两个必须统一的口径

深度 D 在题目里指的是换能器到海底的竖直距离，而不是斜距，也不是“海水深度减去吃水”这种含糊说法。代码里如果 D 取了斜距，覆盖宽度会系统性偏小，而且偏小的比例随坡度变化，在小坡度下看着还行，一到问题三的陡坡段就露馅。

第一个口径是水平投影宽度和沿坡面宽度。重叠率的分子是相邻两条测线覆盖带的重叠长度，分母是本条测线的覆盖宽度，两者必须同口径。用水平投影算间距、用沿坡面宽度算重叠率，结果会稳定偏大一个 1/cosα 的因子，1.5° 坡度下只差 0.03%，看不出问题；10° 坡度下差 1.5%，足以让重叠率从合格区间掉出去。

第二个口径是坡度的正负号。同一个坡面，从东往西推和从西往东推，左右半宽的表达式是对调的。这份代码包的做法是：α 恒为正，约定“测线右侧变深”，递推方向固定从深水侧出发。这个约定写进了函数注释里，改坡度段的时候只需要改数值，不需要改符号。

表里的 D 统一取 100 m、θ 取 120°。最后一行不是笔误：θ/2 + α = 90° 时 cos(θ/2 + α) = 0，那条边界射线与坡面平行，永远打不到海底。这条边界在后面避坑章节还会再出现一次。

2.3 二十行代码把几何层封成函数

几何层我一般写成纯函数，不读文件、不画图、不打印，输入输出全是标量，这样上面三问的脚本都能直接调用，也方便用表格里的数对拍。

# geometry.py import numpy as np def half_widths(D, theta_deg, alpha_deg): """ D : 换能器到海底的竖直深度，单位 m theta_deg : 换能器开角（全角），单位 度 alpha_deg : 海底坡度，单位 度，正号表示测线右侧变深 返回 : (左侧半宽, 右侧半宽)，均为水平投影，单位 m """ h = np.deg2rad(theta_deg) / 2.0 a = np.deg2rad(alpha_deg) # 边界射线与坡面平行或越过坡顶，几何上无交点 if h + abs(a) >= np.pi / 2 - 1e-9: raise ValueError( f"波束与坡面不相交: theta/2={np.rad2deg(h):.2f}°, alpha={alpha_deg}°" ) # 深侧（右侧）斜距更长，半宽更大；浅侧（左侧）相反 right = D * np.cos(a) * np.sin(h) / np.cos(h + a) left = D * np.cos(a) * np.sin(h) / np.cos(h - a) return left, right def coverage_width(D, theta_deg, alpha_deg): """返回 (水平投影全宽, 沿坡面全宽)""" left, right = half_widths(D, theta_deg, alpha_deg) w_h = left + right return w_h, w_h / np.cos(np.deg2rad(alpha_deg))

三个参数里，θ 直接取题目给定的开角，不做半角换算，换算放在函数内部做一次，避免每个调用点都写一遍/2。α 的符号约定是“右侧变深”，这一点必须和使用方的坐标系对齐，否则left和right会被当成两个不同的东西用。

异常分支不是防御性编程的摆设。问题三的坡面如果给出 30° 以上的坡度，这里会直接抛错而不是返回一个负的宽度值。宁可让脚本停在几何层，也不要让一个负宽度流到递推里，最后画出一张看起来正常、实际方向反了的覆盖图。

3. 代码包怎么跑：环境、目录与四个问题的执行链

拿到一份别人写的竞赛代码，最耗时间的往往不是读算法，而是弄清楚哪个脚本先跑、中间产物放在哪、哪些文件是手改的常量、哪些是从附件读的。这份包的结构我按“几何层—求解层—输出层”拆了一遍，跑通顺序基本固定，中间不需要手工改代码。

3.1 依赖环境与目录约定

环境用 Python 3.10 以上就够了，没有用到很新的语法。三个第三方库是硬依赖，绘图和 Excel 读取各多一个。

python -m venv .venv source .venv/bin/activate # Windows 下换成 .venv\Scripts\activate pip install numpy pandas scipy matplotlib openpyxl

requirements.txt在包里，但版本号锁得比较松，只写了主版本；如果本地 numpy 是 2.x，个别老写法可能报弃用警告，按提示改就行，不影响结果。用 MATLAB 的读者也能对上：geometry.py里的半宽公式可以直接搬过去，q2_layout.py里用到的求根换成fzero即可。

3.2 从问题一到问题四的执行顺序

四个脚本之间没有强耦合，除了共用几何层，各自读各自的参数，所以可以单独跑。但第一次上手建议按顺序来，因为问题一的结果是后面三问的校核基准。

# 问题一：单条测线，检查几何层是否与手工计算一致 python q1_single_line.py --theta 120 --depth 120 --alpha 1.5 # 问题二：固定坡度下递推全部测线位置，输出覆盖宽度与重叠率 python q2_layout.py --theta 120 --depth 120 --alpha 1.5 \ --overlap 0.10 --xmin 0 --xmax 4000 \ --out results/q2_lines.csv # 问题三：分段坡面，坡度表写在脚本顶部，直接跑 python q3_slope.py --theta 120 --depth 120 # 问题四：读附件网格，按最深处逐列递推 python q4_grid.py --grid data/attach.xlsx --theta 120 \ --overlap 0.10 --out results/q4_lines.csv

参数里最容易改错的是--overlap。题目要求重叠率落在 10% 到 20% 之间，很多人会顺手填 0.15 取中间值求“稳妥”，结果测线条数比取 0.10 多出三成以上，总长度自然下不来。既然目标是总测线长度最短，就该贴着下限取 0.10，再留一个极小的正余量防止浮点误差把某一段压到 9.99%。

--depth在问题二里是起始测线处的竖直深度，不是海域平均深度。这个参数和--xmin联用：脚本内部按D(x) = depth + x·tanα生成深度剖面，xmin那条测线的深度就是depth。填错了，整条测线序列会整体平移，条数可能刚好对得上，但位置全偏。

--xmax是海域的水平跨度。递推的最后一段通常不会刚好落在边界上，脚本的处理是：若下一条测线的中心越界，就停在边界内侧，并在输出里标记一个edge_flag = 1，提示这一段的实际重叠率高于下限。

3.3 输出结果怎么读

跑完results/下会有几个文件，数值文件和图片要对着看，单看一个容易漏问题。

q2_lines.csv的重叠率列是最直接的体检报告。如果出现负值，说明某两条测线的覆盖带出现了空隙，通常是递推方向或者半宽左右对调导致；如果出现 0.3 以上的值，多半是最后一条测线被边界截断，属于正常，但要确认它确实只在边界处出现。

q4_bands.png那张图我建议放大看边界。实测网格布线在东西两侧边界附近最容易出问题，因为边界列上的深度不一定单调，某几列的极大值会突然变浅，导致按最深处推出来的间距在边界处过宽。

4. 测线间距设计：递推求解、坡面分段与实测网格

问题二到问题四是同一个递推在不同数据形态下的三次套用：坡度固定、坡度分段、深度由网格给定。把递推这一步想清楚，后两个只是换掉深度函数。

4.1 均匀间距为什么会直接翻车

覆盖宽度与深度成正比，深度沿海域线性变化时，测线间距也必须跟着变。举个具体数字：θ = 120°，深度从 100 m 变到 200 m，覆盖宽度从 346 m 变到 693 m。如果按平均宽度取一个固定间距 500 m 铺下去，深水区（D = 200 m）重叠率还有 28%，浅水区（D = 100 m）间距已经超过覆盖宽度本身，中间直接是一条没扫到的缝。

这条缝在俯视条带图上表现为一条贯穿南北的白带，在数值上表现为某相邻两条的重叠率为负。血泪经验是：先跑均匀间距看现象，再换递推，能直观理解为什么这题的核心在“间距随深度变”而不是“算得快”。

4.2 递推求解：每条测线位置用二分法解出来

递推的核心不是“下一条测线挪多远”，而是解一个方程。已知第 i 条测线的位置 x_i，它的右边界是 x_i + right_i；第 i+1 条测线在它左边，左边界是 x_{i+1} − left_{i+1}。重叠率定义为：

overlap = [ (x_i + right_i) − (x_{i+1} − left_{i+1}) ] / (left_i + right_i)

注意 left_{i+1} 和 right_{i+1} 又依赖 x_{i+1} 处的深度，也就是说方程两边都有未知量，没有显式解。这一步用二分法最稳，不需要求导，收敛区间也好定。

# q2_layout.py 核心递推 import numpy as np from scipy.optimize import brentq from geometry import half_widths def next_line_x(x_prev, depth_at, alpha_deg, theta_deg, overlap=0.10): """已知上一条测线位置，解出下一条的位置（向 x 减小方向递推）""" D_prev = depth_at(x_prev) left_prev, right_prev = half_widths(D_prev, theta_deg, alpha_deg) W_prev = left_prev + right_prev edge_right = x_prev + right_prev # 上一条的右边界 def residual(x): # x 是待求的下一条测线中心位置 D = depth_at(x) left, right = half_widths(D, theta_deg, alpha_deg) W = left + right overlap_now = (edge_right - (x - left)) / W return overlap_now - overlap # 目标：正好等于下限 lo, hi = x_prev - 1.0 * W_prev, x_prev - 0.02 * W_prev # 扩大搜索区间直到两端异号，防止极浅水处 residual 没有变号 for _ in range(20): if residual(lo) * residual(hi) < 0: break lo -= 0.5 * W_prev return brentq(residual, lo, hi, xtol=1e-6)

residual里x - left是下一条测线的左边界，edge_right - (x - left)是重叠长度，除以本条覆盖宽度得到重叠率。函数关于 x 单调：x 越小（越往左挪），左边界越靠左，重叠长度越大，残差越大。单调性成立，二分法不会卡在假根上。

lo和hi的初值取W_prev的 1 倍和 0.02 倍作为间距上下界。1 倍覆盖宽度对应重叠率为 0，0.02 倍对应几乎完全重叠，真实解一定落在这个区间内。外面的for循环是给极浅水段准备的：那里 W 很小，深度函数的斜率又可能让残差在初值区间内不变号，向外扩几次就能兜住。

xtol=1e-6意味着位置精度到微米级，对米级的海域来说远超需要，但代价可以忽略，一个 20 条测线的序列跑完不到 0.1 秒。反过来，如果xtol放到 1e-2，最后一条测线的重叠率可能掉到 0.0999，判断题里就算不合格，这种亏吃过一次就记住了。

4.3 坡面分段：交界处的重叠率要单独算

问题三的坡面不是单一坡度，中间有转折。转折点两侧的坡度不同，如果直接把整段当成一个 α 递推，转折点附近的覆盖宽度会算错，表现为交界处那两条测线的重叠率突然跳到 0.25 以上或者掉到 0.05 以下。

处理方式是分段递推、交界处重新起算：在某一段内用该段的 α 递推，推到测线中心越过转折点后，把当前测线作为新段的起点，用新段的 α 继续解下一条。代码里q3_slope.py顶部有一张分段表：

# 每段：(起始 x, 坡度, 该段起点深度) SEGMENTS = [ (0.0, 1.5, 120.0), (1800.0, 3.0, 168.5), (3200.0, 1.0, 241.0), ]

第二段的起点深度不是独立参数，而是由第一段末端推出来的，写在表里是为了让脚本可以单独跑某一段做对比。改坡度时这两处要一起改，否则会出现“坡度变了但深度没接上”的断层，交界处那条测线的覆盖宽度会凭空跳一个台阶。

交界处还有一个细节：如果转折点正好落在某条测线的覆盖带中间，这条测线本身的左右半宽就分属两个坡度。严格做法是左右半宽分别用各自的坡度算，简化做法是取两者中较小的坡度。这份代码用的是前者，q3_slope.py里对跨越转折点的测线单独调用了一次half_widths两次再拼起来。差别在小坡度下不到 1 m，但在题目要求重叠率上下限都卡死的情况下，这 1 m 可能就是把 0.099 顶回 0.101 的关键。

4.4 实测网格：深度取列最大值而不是平均值

问题四给的是网格化的海域深度数据，测线沿南北向布设，每条测线对应一个固定的东西坐标 x，深度沿这条线随南北位置变化。一条测线上深度不是常数，覆盖宽度也就不是常数，这时取哪个深度作为代表，直接决定结果。

取平均深度是常见误用。平均深度偏浅，算出来的覆盖宽度偏小，相邻测线间距被压窄，条数变多，总长度变长，答案看着“更保险”，实际是浪费。正确做法是取该列的最大深度：只要最深处能被覆盖，整条测线就都能覆盖，因为覆盖宽度随深度单调增。

import pandas as pd import numpy as np def load_grid(path): """附件为纯数值网格，行 = 南北向索引，列 = 东西向索引""" raw = pd.read_excel(path, header=None).to_numpy(dtype=float) return raw def build_lines(grid, dx, theta_deg, overlap=0.10): """按列扫描，每列取最大深度，逐列递推测线位置""" col_max = grid.max(axis=0) # 每条候选测线上的最深点 xs, x = [], 0.0 cols = np.arange(grid.shape[1]) * dx # 列索引换算成东西向坐标 while x <= cols[-1]: xs.append(x) D = np.interp(x, cols, col_max) # 该位置的深度用列最大值插值 left, right = half_widths(D, theta_deg, slope_at(x)) W = left + right x += W * (1.0 - overlap) # 简化版递推，正式版仍走二分 return np.array(xs)

dx是网格间距，单位米，必须和坐标换算一致。如果附件给的是“每格 0.01 海里”这类单位，要先换算成米再进递推，混用单位会让np.interp的自变量和因变量对不上，结果是一堆看似合理但位置全错的测线。

grid.max(axis=0)里的axis=0表示沿南北方向压缩，留下每条东西向位置上的最大深度，方向搞反了会得到“每条测线上最深处”的另一个维度，条数直接差一倍。这个错误在图上不容易发现，因为两种布线看起来都挺整齐，只有把总长度和参考值对比才会露馅。

正式版代码里x += W*(1-overlap)这一行会被 4.2 的二分递推替换掉。上面的写法是简化示意，用来说明列的选取和深度的取值方式；直接用简化版跑问题四，重叠率会在边界列附近漂出区间，原因是列最大值随 x 变化不是单调的，线性近似在这里不够用。

5. 避坑与排查：这份代码包里最容易翻车的五个点

跑别人的代码，报错反而是好事，最怕的是不报错但结果不对。下面五条是我在拆这份包时按顺序踩过的，每一条都能在输出文件里找到对应的痕迹。

5.1 覆盖宽度算出负数或爆炸

现象：q3_lines.csv里某几行的左右半宽出现负值，或者数量级突然变成 1e6。

原因：坡度与半开角之和逼近 90°，cos(h + a)趋近于零。几何上那条边界射线与坡面平行甚至越过坡顶，根本没有交点，公式在这个定义域外失效。

解决：在half_widths里保留那道h + abs(a) >= pi/2的判断，别为了“让脚本能跑完”把它注释掉。真遇到超限坡度，正确做法是分段或者改用射线与折线求交，而不是硬套公式。判断阈值留 1e-9 的余量，是因为 89.9999° 时结果已经到 1e5 量级，没有工程意义。

5.2 重叠率整列偏大或偏小一个固定比例

现象：所有相邻测线的重叠率都在 0.13 到 0.14 之间，而目标写的是 0.10。

原因：间距按水平投影宽度算，重叠率却用了沿坡面宽度作分母，两个口径混用。1.5° 坡度下这个偏差约 0.03%，肉眼看不出来；坡度换到 10°，偏差 1.5%，整列就系统性偏了。

解决：全流程只用一种宽度。这份包的约定是几何层返回两个值，但递推和重叠率计算一律用水平投影，沿坡面宽度只在最终报告里作为“实际扫测面积”出现。改代码时在q2_layout.py顶部搜coverage_width，确认拿到的是第一个返回值。

5.3 递推方向反了，条带越推越窄

现象：测线间距从一端到另一端单调变小，最后几条挤在一起，重叠率超过 0.5。

原因：从浅水端往深水端推。覆盖宽度随深度增大，从浅往深推时每一条的间距都要变大，但如果深度函数写反了符号，深度越推越浅，间距就越推越小。

解决：固定从深水端出发，向浅水端递推，x单调递减。这条约定写进函数 docstring，改深度函数时只改depth_at的实现，不改递推方向。判据很简单：输出 csv 里测线间距应该单调不增。

5.4 总测线长度比参考值长两成

现象：布线看起来没问题，覆盖也完整，但总长度就是下不来。

原因：--overlap填了 0.15 取中间值。目标函数是总长度最短，重叠率只要落在 [0.10, 0.20] 内就算合格，取中间值等于主动放弃 5% 的间距。

解决：贴近下限取 0.10，同时留 1e-4 的正余量。余量的作用不是“保险”，而是抵消浮点误差：二分法解出的位置代回去算重叠率，可能得到 0.09999998，四舍五入看着是 0.10，但有些评委的检查脚本用严格小于判断，就判不合格。留这么一点，两边都不吃亏。

5.5 网格行列方向与坐标轴对不上

现象：问题四的测线条数和参考值差一倍，或者布线图旋转了 90°。

原因：pandas读进来的 DataFrame 行索引是南北向，列索引是东西向，但代码里按axis=1取最大值，相当于沿东西方向压缩，得到的“每条测线最深处”变成了另一个方向的极值。

解决：先用两个已知点校核。在附件里找一行一列的手工值，打印grid.max(axis=0)[k]，和肉眼从表里挑出的该列最大值比一次。对得上再往下跑，比跑完再回头查省事得多。这一条我吃过两次亏，第一次是把行列搞反，第二次是单位没换算，两次的表象都是“条数不对”，但排查路径完全不同。

6. 交付前的最后一次校核：把测线位置回代成覆盖率曲线

布线跑完、图也画出来，很容易就点保存交卷了。我现在的习惯是加一道回代：把最终测线位置当成输入，重新算一遍每相邻两条的覆盖区间和重叠率，和求解时输出的值对一遍，两次必须一致到 1e-6。

def audit(xs, depth_at, alpha_at, theta_deg, overlap_lo=0.10, overlap_hi=0.20): """回代校核：返回问题列表，空列表代表通过""" issues = [] for i in range(len(xs) - 1): D = depth_at(xs[i]) l1, r1 = half_widths(D, theta_deg, alpha_at(xs[i])) l2, r2 = half_widths(depth_at(xs[i + 1]), theta_deg, alpha_at(xs[i + 1])) edge_r, edge_l = xs[i] + r1, xs[i + 1] - l2 ov = (edge_r - edge_l) / (l1 + r1) if edge_r < edge_l: issues.append(f"第 {i} 与 {i+1} 条之间有 {edge_l - edge_r:.1f} m 空隙") elif not (overlap_lo - 1e-6 <= ov <= overlap_hi + 1e-6): issues.append(f"第 {i} 与 {i+1} 条重叠率 {ov:.4f} 越界") return issues

这个函数只干一件事：拿最终结果反算，不复用求解过程中的任何中间变量。求解和校核共用同一个函数是好事，共用同一份中间结果就不是了——中间结果一旦错，两边会错得一模一样，校核完全失效。

回代之外还有三张表我每次都会对一遍。第一张是测线条数与总长度，确认总长度等于各段间距之和，差超过 1 m 就说明边界那条测线的处理有偏移。第二张是重叠率的极值分布，正常结果是大量贴近 0.10、少数几条因为边界截断略微偏高，如果出现一堆 0.18 以上，说明某段坡度接错了。第三张是深度范围，测线上的最小深度和最大深度应该和海域实际深度范围吻合，出现负深度或者超出量程的值，基本是单位换算或者行列方向的问题。

这三张表加起来不到二十行代码，但能挡住我见过的大部分低级错误。从那以后我每次交结果前，不管多赶，都强制把audit跑一遍、把这三张表打印出来看一眼——省下的不是时间，是那种交完才发现重叠率整列偏掉的后悔药。压缩包里的脚本按第 3 章的顺序跑一遍，再把这一节的两个校核接上去，基本就能复现出一份能对得上、改得动、讲得清的结果。希望帮到你。

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