★ FEATURED ARTICLE 简介面向备战全国大学生数学建模竞赛、希望复盘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%足以让重叠率从合格区间掉出去。第二个口径是坡度的正负号。同一个坡面从东往西推和从西往东推左右半宽的表达式是对调的。这份代码包的做法是α 恒为正约定“测线右侧变深”递推方向固定从深水侧出发。这个约定写进了函数注释里改坡度段的时候只需要改数值不需要改符号。坡度 α左侧半宽 / m右侧半宽 / m水平覆盖宽度 / m0°173.2173.2346.41.5°173.0174.1347.15°171.0183.5354.510°166.6215.4382.030°——分母为零无解表里的 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目录 / 文件作用是否需要改geometry.py半宽、覆盖宽度、重叠率三个纯函数一般不动q1_single_line.py单条测线的覆盖宽度与面积计算只改命令参数q2_layout.py固定坡度下的测线间距递推只改命令参数q3_slope.py分段坡面的测线布置坡面分段表在脚本顶部q4_grid.py读附件网格数据、按等深线布线需要指定数据文件路径plot_bands.py把测线位置画成覆盖条带图出图用不参与计算results/数值结果 csv 与图片输出自动生成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测线编号、x 坐标、深度、左半宽、右半宽、与前一条的重叠率重叠率列应全部在 [0.10, 0.20]q3_lines.csv同上额外带坡度段编号段与段的交线处重叠率单独检查q4_lines.csv测线编号、x 坐标、该线上最大深度、估计条数最大深度取该列最大值不是平均值q2_bands.png覆盖条带俯视图目视检查有无纵向白缝q4_bands.png实测网格上的布线图结合等深线看测线走向是否合理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第 i1 条测线在它左边左边界是 x_{i1} − left_{i1}。重叠率定义为overlap [ (x_i right_i) − (x_{i1} − left_{i1}) ] / (left_i right_i)注意 left_{i1} 和 right_{i1} 又依赖 x_{i1} 处的深度也就是说方程两边都有未知量没有显式解。这一步用二分法最稳不需要求导收敛区间也好定。# 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, overlap0.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, xtol1e-6)residual里x - left是下一条测线的左边界edge_right - (x - left)是重叠长度除以本条覆盖宽度得到重叠率。函数关于 x 单调x 越小越往左挪左边界越靠左重叠长度越大残差越大。单调性成立二分法不会卡在假根上。lo和hi的初值取W_prev的 1 倍和 0.02 倍作为间距上下界。1 倍覆盖宽度对应重叠率为 00.02 倍对应几乎完全重叠真实解一定落在这个区间内。外面的for循环是给极浅水段准备的那里 W 很小深度函数的斜率又可能让残差在初值区间内不变号向外扩几次就能兜住。xtol1e-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, headerNone).to_numpy(dtypefloat) return raw def build_lines(grid, dx, theta_deg, overlap0.10): 按列扫描每列取最大深度逐列递推测线位置 col_max grid.max(axis0) # 每条候选测线上的最深点 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(axis0)里的axis0表示沿南北方向压缩留下每条东西向位置上的最大深度方向搞反了会得到“每条测线上最深处”的另一个维度条数直接差一倍。这个错误在图上不容易发现因为两种布线看起来都挺整齐只有把总长度和参考值对比才会露馅。正式版代码里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 行索引是南北向列索引是东西向但代码里按axis1取最大值相当于沿东西方向压缩得到的“每条测线最深处”变成了另一个方向的极值。解决先用两个已知点校核。在附件里找一行一列的手工值打印grid.max(axis0)[k]和肉眼从表里挑出的该列最大值比一次。对得上再往下跑比跑完再回头查省事得多。这一条我吃过两次亏第一次是把行列搞反第二次是单位没换算两次的表象都是“条数不对”但排查路径完全不同。6. 交付前的最后一次校核把测线位置回代成覆盖率曲线布线跑完、图也画出来很容易就点保存交卷了。我现在的习惯是加一道回代把最终测线位置当成输入重新算一遍每相邻两条的覆盖区间和重叠率和求解时输出的值对一遍两次必须一致到 1e-6。def audit(xs, depth_at, alpha_at, theta_deg, overlap_lo0.10, overlap_hi0.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} 与 {i1} 条之间有 {edge_l - edge_r:.1f} m 空隙) elif not (overlap_lo - 1e-6 ov overlap_hi 1e-6): issues.append(f第 {i} 与 {i1} 条重叠率 {ov:.4f} 越界) return issues这个函数只干一件事拿最终结果反算不复用求解过程中的任何中间变量。求解和校核共用同一个函数是好事共用同一份中间结果就不是了——中间结果一旦错两边会错得一模一样校核完全失效。回代之外还有三张表我每次都会对一遍。第一张是测线条数与总长度确认总长度等于各段间距之和差超过 1 m 就说明边界那条测线的处理有偏移。第二张是重叠率的极值分布正常结果是大量贴近 0.10、少数几条因为边界截断略微偏高如果出现一堆 0.18 以上说明某段坡度接错了。第三张是深度范围测线上的最小深度和最大深度应该和海域实际深度范围吻合出现负深度或者超出量程的值基本是单位换算或者行列方向的问题。这三张表加起来不到二十行代码但能挡住我见过的大部分低级错误。从那以后我每次交结果前不管多赶都强制把audit跑一遍、把这三张表打印出来看一眼——省下的不是时间是那种交完才发现重叠率整列偏掉的后悔药。压缩包里的脚本按第 3 章的顺序跑一遍再把这一节的两个校核接上去基本就能复现出一份能对得上、改得动、讲得清的结果。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
