蒙特卡羅方法模擬公路自行車運動員的功率及速度分布

摘 要 近些年來，公路自行車賽引起眾多關注，如何使運動員在比賽中獲得更好的成績成為眾多研究的目標。為給出自行車運動員在比賽過程中功率分配方案，本文首先建立了公路自行車賽車手的功率曲線模型，然后分析對比了不同賽車選手的特點;再通過分析賽車手的在不同性質路段的動力學模型，獲得運動員經過具體賽段的時間和速度隨功率的變化。我們運用蒙特卡羅方法得到了不同賽車手在不同路段最佳的功率和速度分配曲線，并將模型應用于2021公路自行車UCI世錦賽，得到不同類型運動員具體功率和速度曲線在路程中的分布。通過這一模型我們可以得知，賽車手在賽程中盡力保持勻速狀態并在擅長的路段進行一定的加速可以得到最佳的比賽結果。最后我們分析了該方法得到的結果的穩定性。

關鍵詞 蒙特卡羅方法;自行車公路賽;自行車運動員功率模型

近幾年，公路自行車賽的風潮席卷全球，引起眾多騎行愛好者和專業騎手的關注，也吸引了許多研究者的注意。越來越多的研究者研究如何使運動員在比賽中獲得更好的比賽成績。

眾多課題組現階段的研究主要是局限于基于參賽運動員個人統計數據的能量管理[1]、營養需求、體重、骨密度[2]等生理性因素對自行車選手比賽發揮的影響，也有對適合自行車性能的最佳海拔高度的研究[3，4]，而缺乏將自行車動力學和運動員的生理性因素相結合的針對自行車個人計時賽這一類高強度運動的研究[5]。

20世紀20年代希爾首先提出臨界功率，后來勞埃爾等人在希爾的研究基礎上，逐步建立功率關于運動時間的嚴謹數學模型[6]。基于此，我們可以建立出運動員功率輸出的大致模型。

密爾沃基心臟研究所提出個人公路自行車計時賽的成績與選手在整個過程中的功率分配密切相關[7]，克瑞斯·阿吉里斯指出不同的比賽項目、賽況條件以及運動員的自身差異性都會對功率分配的結果產生影響[8-10]。但是目前，關于如何更好地在比賽過程中分配功率的分析是十分缺乏的，需要更多的研究來確定比賽的策略，從而使不同運動員在不同的比賽中取得更好的成績。

在實際訓練中，教練員及參賽選手通常是通過反復試驗和經驗依賴來確定最佳的比賽方案，這樣做往往會耗費很多時間和精力而最終效果不佳。因此，使用計算機來模擬不同因素對比賽結果造成的影響，并進行優化是一個簡便優質的方法，并且目前已經有一些課題組進行了類似的研究。2014年，通過對能量供給方和能量需求方進行數學分析將人體運動建模運用于自行車運動[11]。2020年，有團隊模擬了彎道和路況對自行車個人計時賽成績的影響[12]。但是目前關于比賽的模擬結果只是針對理想狀態下功率分布對時間的影響，而忽略了對于賽道的天氣因素以及地形的影響。

在前人的研究基礎上，我們的工作是將運動員的功率輸出模型運用到包含不同地形的具體的自行車賽道上，并通過蒙特卡洛方法求解運動員在最佳功率分配情況。在地形的分析中，我們將地形主要分為平地、彎道和坡道，通過受力分析結合運動員當前的功率得到運動員的運動速度，最終得到運動員總體的運動時間。

由于實際比賽過程是變化的，受多種因素的約束，因此我們將主要影響因素設為變量并且加以條件限制，使用蒙特卡羅方法實現動態優化，而優化目標是找到最短完成時間。

1 計算理論

1.1 運動員個人輸出功率模型

我們首先要建立運動員功率輸出的模型。運動員的功率輸出主要體現在耗氧速率上，與他們的個人身體素質又有著極大關系。

在20世紀20年代希爾首先提出了臨界功率，定義為“長時間可維持的運動強度”，到現在人們認為是“高強度和極高強度運動的臨界點”。臨界功率與個人耗氧的最大速率VO2max 顯著相關。勞埃爾基于希爾的觀測數據建立了跑步者的功率輸出P（t）和比賽總持續時間T 之間的關系[10]

P（t）=（A/T）+R （1）

其中，A 和R 分別表示厭氧代謝的能力和有氧代謝的能量釋放率。然而這個模型，沒有考慮到可持續的有氧能量輸出會逐漸減少，不能準確描述跑步距離超過10000m 時的功率輸出情況。在后續的發展中在勞爾的理論基礎上導出了有氧運動情況下功率隨時間和路程的變化，并且可以表示超過7分鐘左右的運動，即

P =[S/t（ 1-e-t/20） ]+1／t∫t0［BMR +B（1-e-x/30）］dx （2）

其中，S 表示的是無氧運動中運動員實際可以消耗的功率（單位：J/kg），BMR 代表基礎代謝功率，而B 等于瞬時消耗功率峰值和BMR 的差值。該模型很好地預測了60m～42.195km 范圍內的自行車世界紀錄，其中對男性的預測誤差為0.91%，對女性為1.3%。

莫爾頓以數學為工具導出了另一個功率時間關系的模型，建立了一個由磷利用、厭氧糖酵解和有氧動力組成的三組分系統，表示出了生物能量過程的發生，有氧運動中能量供應的延遲。莫爾頓的三組份模型最大可持續有氧運動功率也隨著運動持續時間的延長而下降，但不會最終下降到零，而是維持在氧氣消耗最大速率VO2max 的84%左右。

莫爾頓對一個普通男性建立出如下方程表示該男性的平均輸出功率和運動持續時間的關系[13]

臨界功率并不是無限小，而是約為208W;而當實際功率較大遠遠超過臨界功率時，持續運動時間將很快減小。臨界功率又由于個人身體素質不同而不同，運動員由于長期的科學訓練，臨界功率比普通人更高，而自行車運動員通常為250～350W[14]。

莫爾頓的模型為后來的研究者們打下了扎實的基礎，后來的Omni功率持續時間模型就是基于此發展得來的，如式（4）（5）。Omni功率持續時間模型很好地描述了較長時間內選手的最大輸出功率隨時間的變化，其中考慮了總能量的消耗w'、最大沖刺功率Pmax 和持續時間t 等因素[15，16]。

其中，CP 為臨界功率，W'為在臨界功率以上消耗的能量，Pmax 為功率峰值，CPTTF 為臨界功率下有氧運動持續時間。該模型可以在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]

Ff =Crmg （6）

其中，Cr 為輪胎與地面的滑動摩擦系數，取決于路面和輪胎的材質，在公路自行車比賽中，我們考慮輪胎的材質相同，因此主要與騎行路面的情況有關。

騎行時，前面的空氣由于被壓縮對運動系統產生壓力，兩側的空氣由于相對運動產生，共同構成空氣阻力Fw 。空氣阻力的表達式為

Fw =1／2CdAρV2 （7）

其中，Cd 是空氣阻力系數，A 是最大迎風橫截面積，ρ 是空氣密度，V 為車體相對空氣的速度。

除了上述空氣阻力外，由于非對稱流體的動力效應，車輪旋轉也會產生空氣阻力Fd 。Greenwell研究表明，旋轉阻力的大小取決于車輪的大小及外形，不會隨車輪的轉速發生明顯改變，同時由于人體的作用，后輪所受的阻力也會減少。假定后輪旋轉阻力減少25%，由上述分析，獲得前后車輪轉動時受到的空氣阻力為

其中，Cw 是自行車車輪的空氣阻力系數，r是車輪半徑（前后輪相同），ρ 是空氣密度，V 是車體相對空氣的速度。

根據牛頓第二定律，公路自行車在騎行過程中的動力學方程為

其中，x 代表位移，x'、x″分別代表系統的速度和加速度。If ，Ir 分別代表前輪和后輪的轉動慣量。

在此基礎上，考慮坡角為θ 賽道上的情況。

如圖3所示，此時需要將重力分解為沿運動方向和垂直于運動方向。沿運動方向的重力直接影響運動，垂直運動方向通過改變支持力進而改變摩擦力的大小來影響運動[19]。上坡時，θ 為正，重力分量為阻力;下坡時，θ 為負，重力分量為動力。將支持力表示為F'g ，運動方程為

利用以上的動力學方程，通過MATLAB 利用龍格庫塔算法進行求解。在已知距離和路況的情況下，通過求解該方程可以得到以一定功率通過路段的時間以及過程中的速度和加速度。

1.2.2 自行車彎道動力學模型

假定車身與后輪始終保持在一條直線上，轉彎方向由前輪決定，如圖4所示。

對于縱向方向（x 方向），由向心運動的動力學方程

x″=φ'y'+ax （13）

其中，ax 為沿x 方向由于人的蹬踏產生的加速度。

其中，Fyf 、Fyr 分別為前后輪的橫向作用力，我們將其簡化為簡單的線性輪胎模型，即表達式如下

同樣地，我們使用MATLAB進行求解，獲得以一定功率通過該路段的時間以及該過程速度與加速度的大小。

1.3 程序思路

我們使用MATLAB 軟件編寫程序，并采用蒙特卡羅方法對該問題進行了優化求解。

由以上分析，我們將該模型轉化為求解最小時間的優化問題，并將人的功率曲線以及賽道的性質作為限制條件。為了計算模擬的方便，我們將比賽的總路段劃分為多個足夠小的階段，由于單個小段的路程足夠小，因此可以假定在單個階段的路況是一致的，同樣的，可以假定運動員的功率在這一段是一個定值。同時，在主程序的模擬中，我們暫時忽略了風速對運動員成績的影響，具體原因和解決方案將在后文中討論。

1.3.1 求解時間

將整個賽道分成多個足夠小的階段，因此每個階段只包含一種路況，并將每個階段替換為不同條件下的運動方程。利用MATLAB 中的龍格庫塔算法，用ode45函數求解上述運動微分方程，獲得通過每個路段的速度分布和時間，ode45函數是MATLAB中提供的求常微分方程數值解的函數，表示采用的是四階五階龍格庫塔算法。將每個路段的時間求和，得到總時間。在這個過程中，每個階段的最終速度將作為下一階段的初始速度。

1.3.2 功率選擇

本程序主要采用蒙特卡羅方法對功率進行選擇和優化。由于假設運動員在每小段的功率保持一致，因此只需要隨機確定在該路段的功率即可。在本文中，每段賽道的功率是在功率曲線的限制范圍內使用random 函數隨機選擇，random 函數是生成隨機數的函數，其返回值是一個設定范圍內的隨機數。將此功率作為求解運動微分方程的已知條件，并由此求得該賽道總時間，即最終成績。通過多次迭代，選擇總時間最少的功率組合，并輸出該組功率組合作為最終結果。

1.3.3 功率曲線的改進

最后，我們基于上述結果對功率曲線進行了修正。累積輸出功率的大小會對后續的功率曲線產生影響，如果長時間高功率輸出，將意味著騎手需要額外的時間在較低的功率水平下恢復，這將影響后續的功率分配，通過計算已消耗的能量來對后續的功率曲線進行修正。

2 數值結果

為了設定具體的賽道數值，我們將模型運用到2021公路自行車UCI世錦賽的男子精英賽和女子精英賽中，賽程的路線如圖5所示，不同路段的坡度、長度、彎道的轉角將作為程序的初始條件，本文依據路況將其劃分為36段，具體的路線信息可以在官網獲得。

由于程序求解的方程中有很多參數如轉動慣量、前后輪半徑等，與本文解決的主要問題即功率分配的相關性較弱，因此只取相關參數的大致數值，在實際應用的過程中可以進一步準確設定，基于此，我們得到以下的結果。

2.1 世錦賽不同類型選手模擬結果

對比兩類騎行選手的模擬結果，首先是爬山類型選手，善于在山地進行加速且山地的能量消耗對于后續運動的影響更小;還有一類計時賽專業選手，他們在各種地形有更綜合的優勢。

得到爬山類型選手在世錦賽路段上功率和速度的曲線如圖7所示。

由于計時賽專業選手和爬山類型選手設置的是在不同路況的P 反映曲線和恢復情況不同，所以他們較為舒適的功率曲線應該是一致的。

從整體來看，計時賽專業選手較爬山類型選手的總體功率變化較大。但爬山類型選手的總體速度大于計時賽專業選手。

爬山類型選手在賽程的后半段能夠保持較大的功率輸出，但是計時賽專業選手功率減小很快，這也展示出爬山類型選手擁有長耐力的特性。同時計時賽專業選手在前半段路況較好的時候，選擇較大的功率輸出，也體現出他的經驗。

計時賽專業選手在世錦賽路段上功率和速度的曲線如圖8所示。

沖刺階段，計時賽專業選手提高功率輸出，但是爬山類型選手的輸出功率反而減小，可能是由于爬山類型選手前期長時間功率較高的輸出導致體力下降，而恢復程度不如計時賽專業選手。

在較為復雜的路段，如轉彎、上下坡的時候，計時賽專業選手傾向于通過多次減小輸出功率來維持速度的變化，但是爬山類型選手該作用不顯著，通過藍色的實際輸出功率曲線也可以看到，計時賽專業選手比爬山類型選手更擅長控制自身的功率變化，以適應地形的改變。

2.2 世錦賽不同性別選手模擬結果

在功率的圖中，淺色是騎手能夠維持長時間穩定輸出的功率曲線，而深色代表實際選擇的可以達到最好成績的功率曲線。

從圖9和圖10中我們可以看到，整體上，男性的功率比女性的功率要高，速度也相應地要高。這說明男性能夠在比賽過程中產生更大的速度，成績也會更好。但是，男性功率的下降幅度較女性的大，這說明男性在耐力等方面不如女性。且男性遇到彎道、山坡等路況時輸出功率的表現也不如女性[21]，功率圖中顯示女性的斜率要略大于男性，且持續的路程更長，這代表著女性的沖刺路線長度長于男性，且沖刺的爆發力也越大。這是因為在更大功率的輸出后，男性的恢復能力不如女性。而在研究中表明，人體的質量可以用脂肪質量（fatmass，FM）與瘦體質量（lean mass，LM）這兩個量來描述，其中瘦體指數（lean mass index，LMI）與個體運動的相關性更大。女性由于體重和體脂率的影響，會表現得更輕盈更有耐力，用瘦體指數LMI來衡量女性個體的運動情況會更加準確。[22]

結合路線發現，男性在遇到急彎等復雜路況時更偏向于加大功率來維持自身速度的穩定，而女性的功率表現則相對平穩，同時女性的速度變化也更小。

2.3 程序性質分析

2.3.1 針對風速的敏感性分析

天氣會影響自行車計時賽運動員的運動狀態和速度，斯溫使用普羅斯佩羅等人的數學模型預測[23]，自行車計時賽的最佳比賽策略是讓功率隨坡度和風速變化。斯溫在實際中運用了他的模型來驗證其優越性，當平均功率在丘陵和多風地區增加到相同程度時，它可以在丘陵地區提高速度并節省時間，但是有更多的方案比斯溫的模型能節省更多的時間。在平均功率為289W 的40km 標準賽試跑中，如果功率以接近10%的速率變化，最終可節省約26秒的時間。所以如果我們忽略風速對功率改變的影響，也可以得到最佳功率分配曲線。

馬丁等人導出了自行車行駛速度的模型[24]，涉及的因變量包括空氣阻力、風速、風力、車輪軸承摩擦、勢能和動力學以及機械效率。由此產生的表達式是

P ={V2aVG1/2P CDA +Fw + ——— 空氣阻力VGGRRmTgcos arctan GR +——— 滾動阻力VG 91+8.7v0 VG10-3 +——— 車輪軸承的摩擦損失VGGRmTgsin arctan GR——— 勢能變化1/2 mT +1/r2 V2Gf +V2Gi / ti-tf /Ec——— 能量轉換

研究結果證實，在自行車計時比賽中，如果騎手的力量與風向平行變化，可以縮短比賽時間。也就是說，如果騎手可以通過改變功率來補償風速的變化，從而節省時間。然而，在實際的比賽線路中，由于賽道的方向可能會發生變化，因此行進方向不一定與風速相同，運動員對風的感知以及運動員在運動過程中的反應時間和反應動作很難模仿。而在實驗中發現風速與節省的時間呈線性關系，因此，我們不將風速這個因素納入模型的建立中，但可以在分析結果時考慮風速的影響。

2.3.2 針對運動員實際功率分配的穩定性分析

當運動員判斷失誤時，我們以爬山類型選手為例，在2021公路自行車UCI世錦賽的項目上。我們假設，當遇到山地路段時，爬山類型選手錯誤地使用力量向前移動的可能性有20%。我們假設在這一點上，他將以爬山時能達到的最大功率前進，然后我們使用之前在1.1中建立的個人功率輸出模型，我們可以得到P-x 和V-x 曲線來預測登山者之后的表現。

比較最佳情況下爬山類型選手的輸出功率和速度分布，我們可以看到，爬山類型選手在誤判情況下的最大速度遠高于其他情況下的速度，但誤判花費的總時間比我們模型預測的最佳情況高19.12%。雖然爬山類型選手可以在如此長的距離內恢復到合適的功率輸出情況，但由于對高功率輸出的錯誤判斷，之后的速度會均勻下降，從而影響后續的比賽進度，因此，正如我們在上文計算得到的結果那樣，更均勻地分配功率輸出和速度可以在比賽中獲得更好的表現。

3 結語

3.1 總結

為了幫助騎手在具體線路中制定一個功率分配方案，我們建立了三個模型，即騎手功率模型、自行車動力學模型和個人計時賽的功率分配模型，并將這些模型應用于一條包含不同路況的實際比賽線路，得到了最佳方案，驗證了模型的有效性。基于以上結果，我們通過爬山類型選手和計時賽專業選手在不同路段的表現來分析他們之間的差異。最后，還分析了模型對天氣的敏感性以及實際比賽中當選手不按計算結果分配功率時導致的花費的時間偏差。該模型可作為實際自行車計時賽選手功率分配的參考，并作為教練和運動員規劃初步的路線的依據。

3.2 展望

為了更好地模擬實際線路，獲得更具實際意義的結果，我們在后續工作中可以考慮開發一個程序來識別路況，該程序將自動讀取相關參數，再將其代入模型中進行計算。另外，我們將查閱更多數據，通過添加輪胎材質等可變參數對進行功率優化模型進行進一步的完善，從而實現更符合個人需求的模型定制。

參 考 文 獻

[1] LOUCKS A B， KIENS B， WRIGHT H H. Energy availabilityin athletes[J]. Journal of Sports Sciences， 2011，29（sup1）： S7-S15.

[2] BARRETO G H， OLIVEIRA L， SAITO T， et al. TheRoad to the UCI： Profiling a professional Brazilian femaleroad cycling team[J]. 2020. http：//dx.doi.org/10.31236/osf.io/wnpjk

[3] OLDS T. The optimal altitude for cycling performance： amathematical model[J]. Excel， 1992（8）： 155-155.

[4] POWERS SK， WILLIAMS J. Exercise-induced hypoxaemiain highly trained athletes[J]. Sports Med. 1987， 4： 46-53.

[5] KOLBE T， DENNIS S C， SELLEY E， et al. The relationshipbetween critical power and running performance[J].Journal of Sports Sciences， 1995， 13（3）： 265-269.

[6] PERONNET F， THIBAULT G. Mathematical analysis ofrunning performance and world running records[J]. Journalof Applied Physiology， 1989， 67（1）： 453-465.

[7] FOSTER D C， SCHRAGER M， SNYDER A C， et al. Pacingstrategy and Athletic performance[J]. Sports Medicine， 1994，17（2）： 77.

[8] van INGEN SCHENAU G J， De KONING J J， De GROOTG. The distribution of anaerobic energy in 1000 and 4000metre cycling bouts[J]. International Journal of SportsMedicine， 1992， 13（6）： 447-451.

[9] MARINO F E. Anticipatory regulation and avoidance of catastropheduring exercise-induced hyperthermia[J]. ComparativeBiochemistry and Physiology Part B Biochemistryand Molecular Biology， 2004， 139（4）： 561-569.

[10] FOSTER C， SNYDER A C， THOMPSON N N， et al.Effect of pacing strategy on cycle time trial performance[J]. Medicine and Science in Sports and Exercise， 1993，25（3）： 383-388.

[11] OLDS T. Modelling Human locomotion applications to cycling[J]. Sports Med. 2001， 31（7）： 497-509.

[12] ZIGNOLI A. Influence of corners and road conditions oncycling individual time trial performance and ‘optimal’ pacingstrategy： A simulation study[J]. Proceedings of theInstitution of Mechanical Engineers， Part P： Journal ofSports Engineering and Technology， 2021， 235（3）：227-236.

[13] MORTON R H. A three component model of humanbioenergetics[J]. Math Biology. 1986， 24（4）：451-466.

[14] LEO P， SPRAGG J， PODLOGAR T， et al.Power profilingand the power-duration relationship in cycling： a narrativereview[J]. European Journal of Applied Physiology，2022， 122（2）： 301-316.

[15] PETHICK J， WINTER S L， BURNLEY M. Physiologicalevidence that the critical torque is a phase transition， not athreshold[J]. Medicine amp; Science in Sports amp; Exercise，2020， 52.

[16] VOGT S， SCHUMACHER Y， ROECKER K， et al. PowerOutput during the Tour de France[J]. Sports Med. 2007，28（9）： 756-761.

[17] 蔡睿，陳亮，何申杰，等.公路自行車的動力學分析及仿真研究[J].中國體育科技， 2014， 50（1）： 4.

CAI R， CHEN L， HE S J， et al. Analysis on the dynamicsand simulation of road cycling[J]. China Sport Science andTechnology， 2014， 50（1）： 125-128.

[18] KIM K U， SHIN B S. Modeling motion resistance of rigidwheels[J]. Journal of Terramechanics， 1985， 22（4）： 225-236.

[19] MUELLER A. Modern robotics： Mechanics， planning，and control[bookshelf][J]. IEEE Control Systems Magazine，2019， 39（6）： 100-102.

[20] XU Y， DAI Y， ZHANG H. Strategy arrangement of roadcycling individual time trial based on topology optimizationalgorithm[C]//Proceedings of the 6th International Conferenceon Medical and Health Informatics. 2022： 242-246.

[21] HAAKONSSEN E C， BARRAS M， BURKE L M， et al.Body composition in female road and track endurance cyclists：Normative values and typical changes in female roadand track endurance cyclists[J]. European Journal of SportScience， 2015：1-9.

[22] HOPKINS W G， ANDERSON M E， LEE H， et al. TheLmi： A new practical index of lean mass for monitoringathletes[J]. Medicine amp; Science in Sports amp; Exercise，2002， 34（5）： S237.

[23] ATKINSON G， PEACOCK O， PASSFIELD L. Variableversus constant power strategies during cycling time-trials：Prediction of time savings using an up-to-date mathematicalmodel[J]. Journal of Sports Sciences， 2007， 25（9）： 1001-1009.

[24] James C. Martin， et al. Validation of a mathematical modelfor road cycling power[J]. Journal of Applied Biomechanics，1998， 14（3）.。