蒙特卡罗方法模拟公路自行车运动员的功率及速度分布

摘要

近几年，公路自行车赛的风潮席卷全球，引起众多骑行爱好者和专业骑手的关注，也吸引了许多研究者的注意。越来越多的研究者研究如何使运动员在比赛中获得更好的比赛成绩。

众多课题组现阶段的研究主要是局限于基于参赛运动员个人统计数据的能量管理^[1]、营养需求、体重、骨密度^[2]等生理性因素对自行车选手比赛发挥的影响，也有对适合自行车性能的最佳海拔高度的研究^[3,4]，而缺乏将自行车动力学和运动员的生理性因素相结合的针对自行车个人计时赛这一类高强度运动的研究^[5]。

20 世纪 20 年代希尔首先提出临界功率，后来劳埃尔等人在希尔的研究基础上，逐步建立功率关于运动时间的严谨数学模型^[6]。基于此，我们可以建立出运动员功率输出的大致模型。

密尔沃基心脏研究所提出个人公路自行车计时赛的成绩与选手在整个过程中的功率分配密切相关^[7]，克瑞斯·阿吉里斯指出不同的比赛项目、赛况条件以及运动员的自身差异性都会对功率分配的结果产生影响^[8-10]。但是目前，关于如何更好地在比赛过程中分配功率的分析是十分缺乏的，需要更多的研究来确定比赛的策略，从而使不同运动员在不同的比赛中取得更好的成绩。

在实际训练中，教练员及参赛选手通常是通过反复试验和经验依赖来确定最佳的比赛方案，这样做往往会耗费很多时间和精力而最终效果不佳。因此，使用计算机来模拟不同因素对比赛结果造成的影响，并进行优化是一个简便优质的方法，并且目前已经有一些课题组进行了类似的研究。2014 年，通过对能量供给方和能量需求方进行数学分析将人体运动建模运用于自行车运动^[11]。2020 年，有团队模拟了弯道和路况对自行车个人计时赛成绩的影响^[12]。但是目前关于比赛的模拟结果只是针对理想状态下功率分布对时间的影响，而忽略了对于赛道的天气因素以及地形的影响。

在前人的研究基础上，我们的工作是将运动员的功率输出模型运用到包含不同地形的具体的自行车赛道上，并通过蒙特卡洛方法求解运动员在最佳功率分配情况。在地形的分析中，我们将地形主要分为平地、弯道和坡道，通过受力分析结合运动员当前的功率得到运动员的运动速度，最终得到运动员总体的运动时间。

由于实际比赛过程是变化的，受多种因素的约束，因此我们将主要影响因素设为变量并且加以条件限制，使用蒙特卡罗方法实现动态优化，而优化目标是找到最短完成时间。

1 计算理论

1.1 运动员个人输出功率模型

我们首先要建立运动员功率输出的模型。运动员的功率输出主要体现在耗氧速率上，与他们的个人身体素质又有着极大关系。

在 20 世纪 20 年代希尔首先提出了临界功率，定义为“长时间可维持的运动强度”，到现在人们认为是“高强度和极高强度运动的临界点”。临界功率与个人耗氧的最大速率 V_O2max 显著相关。劳埃尔基于希尔的观测数据建立了跑步者的功率输出 P_(t)和比赛总持续时间T之间的关系^[10]       

莫尔顿以数学为工具导出了另一个功率时间关系的模型，建立了一个由磷利用、厌氧糖酵解和有氧动力组成的三组分系统，表示出了生物能量过程的发生，有氧运动中能量供应的延迟。莫尔顿的三组份模型最大可持续有氧运动功率也随着运动持续时间的延长而下降，但不会最终下降到零，而是维持在氧气消耗最大速率  V_O2max 的 84% 左右。

莫尔顿对一个普通男性建立出如下方程表示该男性的平均输出功率和运动持续时间的关系^[13]

临界功率并不是无限小，而是约为 208W；而当实际功率较大远远超过临界功率时，持续运动时间将很快减小。临界功率又由于个人身体素质不同而不同，运动员由于长期的科学训练，临界功率比普通人更高，而自行车运动员通常为 250~350W^[14]。

莫尔顿的模型为后来的研究者们打下了扎实的基础，后来的 Omni 功率持续时间模型就是基于此发展得来的，如式（4）（5）。Omni 功率持续时间模型很好地描述了较长时间内选手的最大输出功率随时间的变化，其中考虑了总能量的消耗 w′、最大冲刺功率 P_max 和持续时间t等因素^[15，16]。

其中，CP 为临界功率，W′ 为在临界功率以上消耗的能量，P_max 为功率峰值，CP_TTF 为临界功率下有氧运动持续时间。该模型可以在 250W 和 320W 的输出功率下较好地描绘出功率的变化情况，而自行车比赛中输出功率通常大于 300W，故我们采用该模型来预测自行车运动员在运动过程中的功率变化。

以骑手在环法自行车赛的骑行数据为例^[16]，可以知道在不同地形和不同持续时间下最佳输出功率是不同,功率曲线的变化也不同。平地的输出功率约为 218±21W[3.1±0.3W/kg]，在丘陵的输出功率为 228±22W[3.3±0.3W/kg], 在山地约为 234±13W[3.3±0.2W/kg])。在模型中我们分别对不同地形设置基础输出功率和变化范围，在每一次改变地形时随机生成范围内输出功率。

1.2 自行车动力学模型

1.2.1 自行车坡道动力学模型

首先，我们讨论了坡度为 0 的赛道的情况。在骑行过程中，自行车的动力来自于骑行者带动车轮轴转动，阻力则包括来自空气和地面的阻力以及自行车本身的机械摩擦力。沿自行车运动方向和垂直于地面方向建立了二维坐标系，并画出受力示意图，如图 2 所示。

图 2 中，动力 F 由骑手的蹬踏产生。由于在一个蹬踏循环中，力不是恒定的，因此动力用瞬时功率与速度之比表示。同时，人产生的能量也会在链式驱动过程中流失。研究表明，链条的传递效率与驱动力的大小呈正相关，公路自行车一般为 98.5%。^[17]

由于压力使轮胎发生弹性形变，同时轮胎和路面粗糙，车体相对于路面产生摩擦力。骑行过程中，踩动踏板转动链条带动后轮的转动，产生的静摩擦驱动自行车，而对于前轮不直接受到驱动力，而是受到来自地面的摩擦阻力，车体的摩擦力表示为^[18]

骑行时，前面的空气由于被压缩对运动系统产生压力，两侧的空气由于相对运动产生，共同构成空气阻力Fw。空气阻力的表达式为

除了上述空气阻力外，由于非对称流体的动力效应，车轮旋转也会产生空气阻力 Fd。Greenwell 研究表明，旋转阻力的大小取决于车轮的大小及外形，不会随车轮的转速发生明显改变，同时由于人体的作用，后轮所受的阻力也会减少。假定后轮旋转阻力减少 25%，由上述分析，获得前后车轮转动时受到的空气阻力为

其中，Cw 是自行车车轮的空气阻力系数，r 是车轮半径(前后轮相同)，ρ 是空气密度，V 是车体相对空气的速度。

根据牛顿第二定律，公路自行车在骑行过程中的动力学方程为

将以上分析代入方程

其中，x 代表位移，x′、x″ 分别代表系统的速度和加速度。I_f，I_r 分别代表前轮和后轮的转动惯量。

在此基础上，考虑坡角为 θ 赛道上的情况。如图 3 所示，此时需要将重力分解为沿运动方向和垂直于运动方向。沿运动方向的重力直接影响运动，垂直运动方向通过改变支持力进而改变摩擦力的大小来影响运动^[19]。上坡时，θ为正，重力分量为阻力；下坡时，θ为负，重力分量为动力。将支持力表示为F′_g，运动方程为

利用以上的动力学方程，通过 MATLAB 利用龙格库塔算法进行求解。在已知距离和路况的情况下，通过求解该方程可以得到以一定功率通过路段的时间以及过程中的速度和加速度。

1.2.2 自行车弯道动力学模型

假定车身与后轮始终保持在一条直线上，转弯方向由前轮决定，如图 4 所示。

对于纵向方向（x 方向），由向心运动的动力学方程   

同理对于横向方向

  

又有横向加速度

由牛顿第二定律以及扭矩平衡方程得

最终我们有方程组

同样地，我们使用 MATLAB 进行求解，获得以一定功率通过该路段的时间以及该过程速度与加速度的大小。

1.3 程序思路

我们使用 MATLAB 软件编写程序，并采用蒙特卡罗方法对该问题进行了优化求解。

由以上分析，我们将该模型转化为求解最小时间的优化问题，并将人的功率曲线以及赛道的性质作为限制条件。为了计算模拟的方便，我们将比赛的总路段划分为多个足够小的阶段，由于单个小段的路程足够小，因此可以假定在单个阶段的路况是一致的，同样的，可以假定运动员的功率在这一段是一个定值。同时，在主程序的模拟中，我们暂时忽略了风速对运动员成绩的影响，具体原因和解决方案将在后文中讨论。

1.3.1 求解时间

将整个赛道分成多个足够小的阶段，因此每个阶段只包含一种路况，并将每个阶段替换为不同条件下的运动方程。利用 MATLAB 中的龙格库塔算法，用ode45 函数求解上述运动微分方程，获得通过每个路段的速度分布和时间，ode45 函数是 MATLAB 中提供的求常微分方程数值解的函数，表示采用的是四阶五阶龙格库塔算法。将每个路段的时间求和，得到总时间。在这个过程中，每个阶段的最终速度将作为下一阶段的初始速度。

1.3.2 功率选择

本程序主要采用蒙特卡罗方法对功率进行选择和优化。由于假设运动员在每小段的功率保持一致，因此只需要随机确定在该路段的功率即可。在本文中，每段赛道的功率是在功率曲线的限制范围内使用 random函数随机选择，random 函数是生成随机数的函数，其返回值是一个设定范围内的随机数。将此功率作为求解运动微分方程的已知条件，并由此求得该赛道总时间，即最终成绩。通过多次迭代，选择总时间最少的功率组合，并输出该组功率组合作为最终结果。

1.3.3 功率曲线的改进

最后，我们基于上述结果对功率曲线进行了修正。累积输出功率的大小会对后续的功率曲线产生影响，如果长时间高功率输出，将意味着骑手需要额外的时间在较低的功率水平下恢复，这将影响后续的功率分配，通过计算已消耗的能量来对后续的功率曲线进行修正。

2 数值结果

为了设定具体的赛道数值，我们将模型运用到2021 公路自行车 UC I世锦赛的男子精英赛和女子精英赛中，赛程的路线如图5所示，不同路段的坡度、长度、弯道的转角将作为程序的初始条件，本文依据路况将其划分为 36 段，具体的路线信息可以在官网获得。

由于程序求解的方程中有很多参数如转动惯量、前后轮半径等，与本文解决的主要问题即功率分配的相关性较弱，因此只取相关参数的大致数值，在实际应用的过程中可以进一步准确设定，基于此，我们得到以下的结果。

2.1 世锦赛不同类型选手模拟结果

对比两类骑行选手的模拟结果，首先是爬山类型选手，善于在山地进行加速且山地的能量消耗对于后续运动的影响更小；还有一类计时赛专业选手，他们在各种地形有更综合的优势。

得到爬山类型选手在世锦赛路段上功率和速度的曲线如图7所示。

由于计时赛专业选手和爬山类型选手设置的是在不同路况的P反映曲线和恢复情况不同，所以他们较为舒适的功率曲线应该是一致的。

从整体来看，计时赛专业选手较爬山类型选手的总体功率变化较大。但爬山类型选手的总体速度大于计时赛专业选手。

爬山类型选手在赛程的后半段能够保持较大的功率输出，但是计时赛专业选手功率减小很快，这也展示出爬山类型选手拥有长耐力的特性。同时计时赛专业选手在前半段路况较好的时候，选择较大的功率输出，也体现出他的经验。

计时赛专业选手在世锦赛路段上功率和速度的曲线如图 8 所示。

冲刺阶段，计时赛专业选手提高功率输出，但是爬山类型选手的输出功率反而减小，可能是由于爬山类型选手前期长时间功率较高的输出导致体力下降，而恢复程度不如计时赛专业选手。

在较为复杂的路段，如转弯、上下坡的时候，计时赛专业选手倾向于通过多次减小输出功率来维持速度的变化，但是爬山类型选手该作用不显著，通过蓝色的实际输出功率曲线也可以看到，计时赛专业选手比爬山类型选手更擅长控制自身的功率变化，以适应地形的改变。

2.2 世锦赛不同性别选手模拟结果

在功率的图中，浅色是骑手能够维持长时间稳定输出的功率曲线，而深色代表实际选择的可以达到最好成绩的功率曲线。

从图 9 和图 10 中我们可以看到，整体上，男性的功率比女性的功率要高，速度也相应地要高。这说明男性能够在比赛过程中产生更大的速度，成绩也会更好。但是，男性功率的下降幅度较女性的大，这说明男性在耐力等方面不如女性。且男性遇到弯道、山坡等路况时输出功率的表现也不如女性^[21]，功率图中显示女性的斜率要略大于男性，且持续的路程更长，这代表着女性的冲刺路线长度长于男性，且冲刺的爆发力也越大。这是因为在更大功率的输出后，男性的恢复能力不如女性。而在研究中表明，人体的质量可以用脂肪质量（fat mass，FM）与瘦体质量（lean mass，LM）这两个量来描述，其中瘦体指数（lean mass index，LMI）与个体运动的相关性更大。女性由于体重和体脂率的影响，会表现得更轻盈更有耐力，用瘦体指数LMI来衡量女性个体的运动情况会更加准确。^[22]

结合路线发现，男性在遇到急弯等复杂路况时更偏向于加大功率来维持自身速度的稳定，而女性的功率表现则相对平稳，同时女性的速度变化也更小。

2.3 程序性质分析

2.3.1 针对风速的敏感性分析

天气会影响自行车计时赛运动员的运动状态和速度，斯温使用普罗斯佩罗等人的数学模型预测^[23]，自行车计时赛的最佳比赛策略是让功率随坡度和风速变化。斯温在实际中运用了他的模型来验证其优越性，当平均功率在丘陵和多风地区增加到相同程度时，它可以在丘陵地区提高速度并节省时间，但是有更多的方案比斯温的模型能节省更多的时间。在平均功率为 289W的 40km 标准赛试跑中，如果功率以接近 10% 的速率变化，最终可节省约 26 秒的时间。所以如果我们忽略风速对功率改变的影响，也可以得到最佳功率分配曲线。

马丁等人导出了自行车行驶速度的模型^[24]，涉及的因变量包括空气阻力、风速、风力、车轮轴承摩擦、势能和动力学以及机械效率。由此产生的表达式是

研究结果证实，在自行车计时比赛中，如果骑手的力量与风向平行变化，可以缩短比赛时间。也就是说，如果骑手可以通过改变功率来补偿风速的变化，从而节省时间。然而，在实际的比赛线路中，由于赛道的方向可能会发生变化，因此行进方向不一定与风速相同，运动员对风的感知以及运动员在运动过程中的反应时间和反应动作很难模仿。而在实验中发现风速与节省的时间呈线性关系，因此，我们不将风速这个因素纳入模型的建立中，但可以在分析结果时考虑风速的影响。

2.3.2 针对运动员实际功率分配的稳定性分析

当运动员判断失误时，我们以爬山类型选手为例，在 2021 公路自行车 UCI 世锦赛的项目上。我们假设，当遇到山地路段时，爬山类型选手错误地使用力量向前移动的可能性有 20%。我们假设在这一点上，他将以爬山时能达到的最大功率前进，然后我们使用之前在 1.1 中建立的个人功率输出模型，我们可以得到P-x 和 V-x 曲线来预测登山者之后的表现。

比较最佳情况下爬山类型选手的输出功率和速度分布，我们可以看到，爬山类型选手在误判情况下的最大速度远高于其他情况下的速度，但误判花费的总时间比我们模型预测的最佳情况高 19.12%。虽然爬山类型选手可以在如此长的距离内恢复到合适的功率输出情况，但由于对高功率输出的错误判断，之后的速度会均匀下降，从而影响后续的比赛进度，因此，正如我们在上文计算得到的结果那样，更均匀地分配功率输出和速度可以在比赛中获得更好的表现。

3 结语

3.1 总结

为了帮助骑手在具体线路中制定一个功率分配方案，我们建立了三个模型，即骑手功率模型、自行车动力学模型和个人计时赛的功率分配模型，并将这些模型应用于一条包含不同路况的实际比赛线路，得到了最佳方案，验证了模型的有效性。基于以上结果，我们通过爬山类型选手和计时赛专业选手在不同路段的表现来分析他们之间的差异。最后，还分析了模型对天气的敏感性以及实际比赛中当选手不按计算结果分配功率时导致的花费的时间偏差。该模型可作为实际自行车计时赛选手功率分配的参考，并作为教练和运动员规划初步的路线的依据。

3.2 展望

为了更好地模拟实际线路，获得更具实际意义的结果，我们在后续工作中可以考虑开发一个程序来识别路况，该程序将自动读取相关参数，再将其代入模型中进行计算。另外，我们将查阅更多数据，通过添加轮胎材质等可变参数对进行功率优化模型进行进一步的完善，从而实现更符合个人需求的模型定制。

参考文献

通信作者:  景鹏飞，北京师范大学物理学系助理研究员，pfjing@bnu.edu.cn。