反舰导弹分段弹道建模与MATLAB复现：从四段曲线到三维参数方程

简介：第八届MathorCup数学建模竞赛优秀论文C3289.pdf聚焦陆基导弹打击航母的数学建模问题，面向数学建模参赛者、军事运筹与弹道仿真爱好者，完整呈现特等奖论文的建模思路与算法设计。论文从静态轨道模型入手，针对发射段、中段、末段分段建立数学表达式，并利用MATLAB绘制轨道曲线；随后扩展至动态目标场景，引入三维坐标和参数方程，通过Excel生成散点并经SPSS绘制三维轨迹；最后综合风速、气压、控制延迟等因素，分析系统误差、随机误差并计算命中率。资源仅含1个PDF文件，大小约511KB，内容为完整赛题论文，包含问题重述、模型假设、分段参数方程、图表与结果展示，可直接用于备赛参考、论文研读和建模方法复盘。目前已有106人学习下载，适合希望提升复杂系统建模与数值仿真能力的读者。

1. 这份数学建模论文的价值：把一条弹道拆成四段曲线

拿到这份第八届 MathorCup 竞赛优秀论文 PDF 时，我原以为里面就是一堆密密麻麻的公式截图，翻完才发现它把"反舰导弹轨道设计"这件事拆成了可以照着复现的三步：先对静止目标做分段弹道建模，再对运动目标换成三维参数方程，最后落到误差与命中率计算。最打动我的是它的分段思路——发射段、中段两个过程、末段，每一段受力不同，必须单独建立表达式再拼接，而不是用一条曲线硬套全段。这份资源适合两类人：准备数学建模竞赛的学生，以及做目标跟踪、飞行器轨迹仿真相关工作的从业者。我去年带学生复现类似题目时就吃过亏：直接追求单一连续弹道方程，结果发射段的抛物线形态和末段的平抛形态根本合不到一起，这份论文恰好给出了最省事的解法。

2. 把静态弹道拆成四段：受力分析与分段模型

2.1 为什么不能一条公式走到底

弹道建模最常见的翻车点，是试图用一个连续函数描述导弹全程。实际上导弹在飞行不同阶段的受力结构差异极大：发射段导弹竖直初速度 500 m/s 向上，同时施加一个水平指向目标的牵引力，整体做类斜抛运动，轨迹是抛物线；进入中段后水平牵引力消失，只保留竖直方向的牵引力调整高度，水平速度保持不变；末段则完全撤去推力，导弹依靠水平初速度做平抛运动，借助重力完成末段下压。

这三段受力不同，运动学方程的形式也完全不同——发射段是 v-t 和抛物线消参，中段是两个竖直分运动的拼接，末段是典型的平抛。硬用一个高阶多项式去拟合整条弹道，不仅物理含义说不清，参数回归也会很勉强。这份论文的聪明之处在于承认分段事实，把总轨迹写成带定义域的分段函数，再用"曲线连续光滑"作为拼接约束。这是做物理建模的第一原则：先看受力再看方程，而不是先选函数形式。

2.2 关键参数是怎么定出来的

要复现这个模型，先把几个关键数字的来源理清楚。以下是论文里通过运动学关系反推出来的核心参数表，建议复现时逐一核对：

这些参数不是凭空设定的，而是通过"已知末段时间、末段最大速度、重力加速度、发射段初速度"四个条件联立反推得到的。复现时最容易出错的是末段水平速度：论文里用动能定理 mg·h3 = 0.5·m·(v末² − v水平²)，代入 v末 = 1000 m/s、h3 = 1960 m，解得 v水平 ≈ 980.6 m/s。注意这个 980.6 代表末段水平初速度，它和发射段末端的水平速度在严格物理意义上应当连续衔接，这一点我在后面避坑章会专门展开。

2.3 MATLAB 复现分段轨道：匿名函数与拼接绘制

拿到分段表达式后，最直接的做法是在 MATLAB 里用匿名函数逐段定义，再按定义域拼接绘制。问题 1 的总轨道分段函数可以用下面这段代码实现：

% 分段弹道曲线定义，x 单位米，y 单位米 f1 = @(x) 224.06 * sqrt(x) - 0.9839 * x; % 发射段：0 <= x < 13003 f2 = @(x) -1.0275e-4 * x.^2 + 0.2672 * x + 11017.7235; % 中段第一过程：13003 <= x < 16038 f3 = @(x) 8.334e-6 * x.^2 - 6.25175e-3 * x + 13661.9198; % 中段第二过程：16038 <= x < 359206 f4 = @(x) -5.0958e-11 * x.^3 + 3.66088e-5 * x.^2 - 3.65546e-1 * x + 655546; % 末段 x1 = linspace(0, 13003, 200); x2 = linspace(13003, 16038, 200); x3 = linspace(16038, 359206, 200); x4 = linspace(359206, 378818, 200); hold on; plot(x1, f1(x1), 'b-', 'LineWidth', 1.5); plot(x2, f2(x2), 'r-', 'LineWidth', 1.5); plot(x3, f3(x3), 'g-', 'LineWidth', 1.5); plot(x4, f4(x4), 'm-', 'LineWidth', 1.5); xlabel('水平距离 x / m'); ylabel('高度 y / m'); legend('发射段', '中段第一过程', '中段第二过程', '末段'); grid on;

代码逻辑上，四段函数依次对应发射段抛物线、中段两段曲线和末段下压曲线。linspace的第三个参数控制每段的采样点数，200 个点足够让曲线光滑，又不至于让文件过大。绘制后有一个值得做的检查：把分界点 x = 13003、16038、359206 分别代入相邻两个函数，看高度差是否接近 0——这就是"连续光滑"约束的数值验证。如果某个分界点出现明显台阶，说明系数抄错了或者定义域写错了。

3. 动态目标模型：从二维平面到三维参数方程的改进

3.1 思路一为什么被舍弃

问题 2 要求打运动中的航母，航速 32 节，航向正南。最先想到的思路往往是这样：既然航母速度和方向已知，那就先算出导弹飞行 t 秒后航母在哪，直接把导弹打击终点设为那个位置。这个思路看起来顺理成章，实际却有一个致命缺陷——飞行时间 t 本身就是未知量，需要联立求解；更关键的是，指挥中心无法在导弹发射前准确预知 t 时刻航母的最终位置，因为航母的实时坐标要由卫星和无人机持续提供，打击过程是动态修正的，不是发射前一次解算完的。

论文里也验证了这个思路的方程复杂度：水平距离 L 需要写成 L = √(a + 16.4622t)² + b² 的形式，其中 a 是经线方向间距、b 是纬线方向间距，再与发射段时间 51.02 s、末段时间 20 s、中段时间一起联立解 t。这个模型即便数学上可解，也违背了"动态制导"的实际工作方式，所以被舍弃。这个判断对做工程仿真很有启发：模型不是越复杂越好，而是要和实际系统的信息流匹配。

3.2 三维坐标系与发射段参数方程

改进模型的关键一步是升维。论文以导弹初始位置为原点，导弹初始位置与航母初始位置连线方向为 Y 轴，水平面上从 Y 轴顺时针旋转 90° 为 X 轴，垂直水平面向上为 Z 轴，建立右手坐标系。这样把航母的正南运动分解到 X 和 Y 两个方向：导弹初始位置与航母初始位置的连线与正南方向的夹角记为 θ，计算得 sinθ = 0.83673391，cosθ = 0.24672125。

初始速度分解为：Z 方向分速度 500 m/s（与问题 1 相同），正南方向分速度 v = 16.46 m/s（32 节换算）。发射段三个方向的参数方程如下，代码里直接按时间 t 生成采样点：

% 发射段参数方程，t 单位秒，x/y/z 单位米，0 <= t <= 51.02 t1 = linspace(0, 51.02, 300); x1 = 13.7745 * t1; % X 方向：v·sinθ 匀速 y1 = 4.98 * t1.^2 + 4.06158 * t1; % Y 方向：初速 v·cosθ + 水平牵引加速 z1 = -4.9 * t1.^2 + 500 * t1; % Z 方向：初速 500 向上减速 % 发射段终点 fprintf('发射段终点坐标: (%.3f, %.3f, %.3f)\n', ... x1(end), y1(end), z1(end)); % 输出：发射段终点坐标: (702.824, 13167.278, 12754.935)

注意第三行 Z 方向的表达式，−4.9·t² 是重力加速度一半的负值，500·t 是竖直初速度的效果，两者叠加正是竖直上抛运动。Y 方向的 4.98·t² 来自水平牵引力加速度的一半，4.06158·t 是正南分速度贡献。三个方向分别独立写成参数方程，比二维的 y = f(x) 形式信息量更大，也为后续中段和末段的拼接提供了统一的时间基准。

3.3 中段到末段：时间平移与三维散点可视化

中段第一过程从发射段终点开始，但参数方程在推导时是以"中段起点为时间零点"写的，绘图前需要做时间平移。论文给出的核心方程中，X 方向依旧是匀速 13.7745·t，Y 方向变为 984.66158·t − 37067.071，Z 方向变为 −9.88·t² + 1008.1552·t − 12962.9372，定义域为 51.02 < t ≤ 54.11476。这里把平移后的常数项直接做进了系数里，复现时最容易抄错的是那三个常数项——它们是由发射段终点坐标代入后合并出来的。

由于动态模型是三维曲线，MATLAB 里可以用plot3直接看三维轨迹：

% 中段第一过程参数方程，t 单位秒，51.02 < t <= 54.11476 t2 = linspace(51.02, 54.11476, 300); x2 = 13.7745 * t2; y2 = 984.66158 * t2 - 37067.071; z2 = -9.88 * t2.^2 + 1008.1552 * t2 - 12962.9372; plot3(x1, y1, z1, 'b-', 'LineWidth', 1.5); hold on; plot3(x2, y2, z2, 'r-', 'LineWidth', 1.5); xlabel('X / m'); ylabel('Y / m'); zlabel('Z / m'); grid on;

论文里采用的方式是用 Excel 按时间步长生成一系列三维散点，再导入 SPSS 绘制三维散点图，这样做的优势是可视化直观，能清楚看到轨迹从发射段过渡到中段时 X、Y、Z 三个方向的变化趋势。用 MATLAB 复现时，plot3更省事，但要注意把每段的采样时间接上，否则三维图上会看到轨迹在中段起点处断掉。我的习惯是用linspace让两段时间轴首尾重叠一个点，这样图上的曲线是连续闭合的。

4. 复现避坑记录：单位、假设与数据自洽性

4.1 单位换算翻车：32 节到底是多少米每秒

现象：把 32 节直接当 32 m/s 代入参数方程，绘制出的弹道整体偏离目标位置，Y 方向位移明显偏大。

原因：节是航海速度单位，1 节 = 1.852 km/h，32 节换算成国际单位需要先乘 1.852 得 59.264 km/h，再除以 3.6 才是 16.46 m/s。直接代入 32 等于把速度放大了近一倍。

解决：建模前先做单位统一表，所有速度一律换算为 m/s，所有距离一律换算为米。论文里给出了换算结果 v = 16.46 m/s，但如果你拿到的是原始题面数据，这一步必须自己算一遍，并且换算后用 sinθ、cosθ 分量回代验证发射段终点的 Y 坐标是否落在合理区间。

4.2 假设与抛物线形态的矛盾

现象：忽略空气阻力后，模型里发射段轨迹是抛物线，但有人尝试加入空气阻力系数后画出完全不同的曲线，和论文图对不上。

原因：抛物线运动的前提是合力大小和方向恒定。空气阻力总是与速度反向，大小随速度变化，一旦计入，合力就不再恒定，轨迹不可能保持抛物线形态。

解决：这个模型从一开始就建立在"忽略空气阻力"的假设上，论文里明确写了这一点。复现时不要中途私自加阻力项，否则需要重推所有分段函数。如果要做阻力敏感性分析，正确的做法是在论文原始模型基础上单独建立一个带阻力的对照模型，而不是在原模型里混入阻力项。

4.3 分段点处速度和位置的连续性问题

现象：把四段分段函数分别求导后，发现 x = 13003 处左右两侧斜率差异很大；更进一步，发射段末端的水平速度约 508 m/s，而末段反推出的水平初速度约 980.6 m/s，差距明显。

原因：论文的建模在位置坐标上保证了分段函数连续，但"连续光滑"要求一阶导数也连续，当分段边界处的受力状态突然改变时，速度跳变难以避免；末段 980.6 m/s 是由动能定理结合 1000 m/s 末速度反推出来的，它和发射段末端的水平速度不是同一个来源。

解决：复现时把速度连续性作为单独的校验项。检查方法是分别对相邻两段求导，比较分界点处的导数值。若出现明显跳变，说明中间可能缺失了一个加速或减速过程，需要补一段变力过程，或者修正牵引力作用时间。这个校验能帮你发现论文里没有展开说明的过渡态，也算是对原始模型的合理补全。

4.4 延迟时间量的量纲陷阱

现象：命中率计算时，把控制系统延迟时间 Δtd 代入位移公式，得到的结果单位对不上，数值大得离谱。

原因：延迟时间的单位是秒，位移 = 速度 × 时间，如果速度用了 km/h、时间用了秒，结果就变成千米级错误；如果延迟时间按毫秒计算而没有除以 1000，也会放大三个数量级。

解决：统一用米和秒计算，先写出位移公式 Δs = v × Δtd，再检查量纲。比如水平速度 980.6 m/s、延迟 0.1 s，位移约 98 m，这个量级和航母 335 m 的舰长相比才有意义。我在复现时会把延迟时间换算成秒之后单独存一个变量，不在公式里反复写数字。

4.5 中段时间基准错位

现象：绘制三维轨迹时，发射段曲线正常，但从中段起点开始整体平移，明显错位。

原因：中段参数方程的推导以中段起点为时间零点，直接把它放进以发射时刻为 0 的时间轴上，所有坐标都会平移 51.02 秒对应的量。

解决：建模时标注每一段的时间基准。发射段是 0 ≤ t ≤ 51.02，中段表达式里的 t 需要整体平移，实际绘图坐标是 t_实际 − 51.02。论文给出的系数里把平移常数合并进去了，看起来简洁，但复现时容易忽略它隐含的时间偏移。建议在代码注释里写明每一段的时间基准和定义域，避免混用。

5. 验证整体弹道的关键技巧：边界点三查

拿到这套模型后，指望一次性画对并不现实。我的经验是复现任何分段弹道模型，都要做一次"边界点三查"：查位置连续、查速度连续、查终点命中。下面这段代码把三查集中在一个脚本里，建议每次改完参数都跑一遍：

% 边界点三查：位置、速度、终点 % 位置连续：相邻函数在分界点处的值应近似相等 x_bounds = [13003, 16038, 359206, 378818]; tol = 1e-3; % 容差，单位米 for i = 1:length(x_bounds)-1 y_left = feval(f_cell{i}, x_bounds(i)); y_right = feval(f_cell{i+1}, x_bounds(i)); fprintf('分界点 x=%.0f 高度差: %.4f m\n', x_bounds(i), abs(y_left - y_right)); end % 速度连续：对比相邻函数在分界点的导数值 syms x; df1 = diff(f1(x), x); df2 = diff(f2(x), x); v1 = double(subs(df1, x, 13003)); v2 = double(subs(df2, x, 13003)); fprintf('发射段/中段斜率差值: %.4f\n', abs(v1 - v2)); % 终点命中：末段终点坐标应接近目标位置 target_y = 378818; % 目标水平距离 endpoint = f4(target_y); fprintf('末段终点高度: %.2f m（应接近 0）\n', endpoint);

三查的逻辑很简单：位置连续是对"轨道曲线连续光滑"的直接验证，速度连续能暴露过渡态缺失，终点命中则确认整条弹道确实落在了目标航母坐标上。我在第一次复现这份论文时就因为只画图、没查边界，在中段第一过程的常数项里抄错了一个数字，导致 13003 米处分界点高度差了将近 3 千米——如果当时先跑一遍三查脚本，这个错误一眼就能看出来。从那以后，我每次做这类分段轨迹复现，都强制把三查脚本放在代码最前面，先过校验再谈画图。这套习惯帮我省下的排查时间，远超过写校验脚本本身那十分钟。希望帮到你。

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