崔康靖 鄭元洲 陳國成 胡衛(wèi)東

(武漢理工大學(xué)航運(yùn)學(xué)院1) 武漢 430063) (內(nèi)河航運(yùn)技術(shù)湖北省重點(diǎn)實(shí)驗(yàn)室2) 武漢 430063)(湄洲灣港引航站3) 莆田 351100)

0 引 言

隨著經(jīng)濟(jì)全球化，海運(yùn)量已占據(jù)全球商品貿(mào)易的80%以上[1]，而由氣象因素導(dǎo)致的海難事故時有發(fā)生．在氣象航線及其優(yōu)化方面國內(nèi)外學(xué)者已有眾多研究．Park等[2]通過分階段來優(yōu)化油耗問題，在航速恒定下用A*算法求解轉(zhuǎn)向點(diǎn)位置，結(jié)合幾何規(guī)劃法重新分配各航路航速．Zaccone等[3]提出了一種基于三維動態(tài)規(guī)劃的船舶航線優(yōu)化方法，通過分配航路點(diǎn)尋求最優(yōu)航線．張進(jìn)峰等[4]針對我國近海船舶進(jìn)行了避臺航線優(yōu)化設(shè)計，基于動態(tài)規(guī)劃算法確定各航路點(diǎn)．在航線中船舶各航段的航速優(yōu)化方面，李明峰等[5]以浪高為約束條件構(gòu)建氣象航線模型，引入懲罰函數(shù)應(yīng)用于遺傳進(jìn)化算法進(jìn)行優(yōu)化．劉洋等[6-7]擬定危險風(fēng)區(qū)為障航物，提出一種基于向量的多步滾動窗口優(yōu)化算法進(jìn)行航線設(shè)計．

上述研究在模擬海上氣象時考慮的氣象條件都比較單一，應(yīng)用算法在規(guī)避惡劣氣象區(qū)域時也不易設(shè)置合適的安全距離．針對大風(fēng)浪問題，文中提出了一種結(jié)合人工勢場法和模擬退火算法的優(yōu)化方法．該方法在考慮風(fēng)速和浪高的雙重氣象因素下對各航段航速影響，并以航程和航行時間為優(yōu)化目標(biāo)，利用人工勢場法能較好的避讓大風(fēng)浪區(qū)，并通過模擬退火算法替代原有的負(fù)梯度下降法對相應(yīng)的航路點(diǎn)位置進(jìn)行調(diào)整，增加可選航路點(diǎn)位置來優(yōu)化人工勢場法航線，以尋求最優(yōu)氣象航線．

1 船舶氣象航線設(shè)計的數(shù)學(xué)模型

1.1 航線模型定義

在考慮船舶航線設(shè)計時設(shè)定為整條航線由多個航路點(diǎn)組成，即將整條航線分為多個航路段，每段航路以恒向線方式表示，包含該段航路端點(diǎn)的航行速度、經(jīng)緯度位置和航向角等信息．

設(shè)定航線中有n個航路點(diǎn)，用向量R表示．各個航路點(diǎn)的位置坐標(biāo)為

R={(Lat1,Long1),…,(Lati,Longi),…,

(Latn,Longn)}

式中：Lati為第i個航路點(diǎn)坐標(biāo)的緯度；Longi為第i個航路點(diǎn)坐標(biāo)經(jīng)度．

船舶在每段航路間將以恒向方式航行，航線中各航段的航向變量φ為

φ=[φ1,…,φi,…,φn]

(1)

式中：φi為第i段航路的航向z，由第i個和第i+1個航路點(diǎn)經(jīng)緯度坐標(biāo)決定，φi∈[0°,360°]，以正北方向?yàn)?°，順時針方向遞增．

V0為船舶在靜水中航行速度，而在實(shí)際航行中船舶容易受到風(fēng)浪流等外界因素的影響而導(dǎo)致速度發(fā)生變化，因此用向量V表示各航路段航速，即

V=[V1,…,Vi,…,Vn]

(2)

式中：Vi為第i段航路的航速．

1.2 風(fēng)浪對船舶速度的影響

船舶在海上航行時，受到風(fēng)、浪、涌等外界環(huán)境因素影響，船舶行駛受到的阻力增加，在主機(jī)輸出功率不變時，船舶實(shí)際航速低于靜水時航速，這種現(xiàn)象稱為船舶失速．采用文獻(xiàn)[8]以最小二乘法為基礎(chǔ)迭代計算船舶失速公式.

Va=V0-(1.08h-0.126qh+2.77×10-3×

Fcosα)(1-2.33×10-7DV0)

(3)

式中：Va為船舶在風(fēng)浪中實(shí)際速度，kn；V0為所給定的船舶靜水速度，kn;F為風(fēng)速大小，m/s；D為船舶實(shí)際排水量，t;h為有義波高，m;q為船舶航向與浪向間的夾角，rad；α為船舶航向與風(fēng)向之間的相對夾角，rad．

1.3 大風(fēng)浪區(qū)的界定

在海上航行時，不同船舶所能應(yīng)對的風(fēng)浪情況并不相同．根據(jù)船舶的船型、噸位、結(jié)構(gòu)強(qiáng)度、操作性能，以及裝載貨物等，風(fēng)浪對船舶造成的影響會有很大差異．因此，在大風(fēng)浪區(qū)的具體界定時，應(yīng)當(dāng)考慮該船的現(xiàn)實(shí)情況，只要某片海域的風(fēng)浪環(huán)境超出該船的承受能力，便能將其定義為大風(fēng)浪區(qū)域．船舶抗風(fēng)浪等級的公式為[9]

(4)

h=2(M0/g+K1DLw1-K2W0Lw1-

(5)

式中：V10為船舶在改裝載狀態(tài)下可安全航行的最大風(fēng)速；D為船舶排水量；lq為最小傾覆力臂；Af為船舶受風(fēng)面積；K為穩(wěn)性衡準(zhǔn)數(shù)；Zf為風(fēng)作用力臂；h為船舶船體強(qiáng)度所能承受的最大浪高；M0為船舶的船中彎矩；g為重力加速度；K1為平水排水量力矩系數(shù)；Lw1為船舶水線長；K2為空船重量力矩系數(shù)；W0為空船重量；∑PiXi為所有載荷對船中的力矩；K3為波形排水量力矩修正系數(shù)；B為船寬．

2 動態(tài)航線規(guī)劃模型

2.1 目標(biāo)函數(shù)

1) 航線航程 根據(jù)前文中航線模型的設(shè)定，航線總航程應(yīng)為各個航路段相加之和，為

(6)

式中：L為總航線航程；Li為各分段航程．

在求航線航程時，需要根據(jù)坐標(biāo)點(diǎn)的經(jīng)緯度信息求實(shí)際航行距離，由前文介紹在各航路段中航線表示為恒向線，選用墨卡托算法計算航程，為

(7)

L=(Lat2-Lat1)×secφ

(8)

式中：Lat1,Long1分別為第一個點(diǎn)的緯、經(jīng)度坐標(biāo)；Lat2,Long2分別為第二個點(diǎn)的緯、經(jīng)度坐標(biāo)；φ為恒向線方向；L為兩點(diǎn)間航程；e為地球偏心率．

2) 航線航時 從起點(diǎn)到終點(diǎn)的總航行時間可由各航路段航行時間求和得到，為

(9)

式中：Ttotal為總的航行時間；Vi為第i段的實(shí)際航速．

2.2 人工勢場算法

文中在模型仿真中采用的是人工勢場算法，其原理為人通過引力場和斥力場共同作用使物體能夠沿梯度下降最快的方向運(yùn)動，其中目標(biāo)點(diǎn)對物體產(chǎn)生引力，引導(dǎo)物體朝向其運(yùn)動．障礙物對物體產(chǎn)生斥力，避免物體與之發(fā)生碰撞．物體在路徑上每一點(diǎn)所受的合力等于該點(diǎn)所有斥力和引力之和．具體的算法模型如下.

常見的引力場：

(10)

式中：ξ為引力尺度因子；ρ(q,qgoal)為物體當(dāng)前狀態(tài)與目標(biāo)的距離．

引力即為引力場對距離的導(dǎo)數(shù)：

Fatt(q)=-Uatt(q)=ρ(qgoal-q)

(11)

斥力場：

(12)

式中：η為斥力尺度因子；ρ(q,qobs)為物體與障礙物之間的距離；ρ0為每個障礙物的影響半徑．

斥力為斥力場的導(dǎo)數(shù)：

(13)

引力和斥力尺度因子的取值對人工勢場算法航路規(guī)劃結(jié)果影響較大，如果取值不當(dāng)，會導(dǎo)致航路不可行或航路點(diǎn)冗余[10]．取ξ=10且保持不變，η分別取5、8、10、12、15、18時人工勢場算法規(guī)劃結(jié)果見圖1．

圖1 η取不同值時對航線規(guī)劃的影響

由圖1可知，當(dāng)η取5、8、10時，會出現(xiàn)航路途徑大風(fēng)浪障礙區(qū)的情形，規(guī)劃航線不可取；而相比η=12，當(dāng)η取15和18時，發(fā)現(xiàn)航路點(diǎn)有明顯冗余現(xiàn)象，因此本文參數(shù)選擇取ξ=10，η=12，此時人工勢場算法得到的規(guī)劃航路最佳．

2.3 模擬退火算法

經(jīng)由人工勢場法得到航線并非是最優(yōu)航線，船舶駛于航線中途的航路點(diǎn)時，一方面可能會面臨引力和斥力的合力矢量為零陷入局部最小而終止路徑規(guī)劃．另一方面人工勢場法可能受引、斥力增益系數(shù)的影響造成迭代步伐過大，引起抖動現(xiàn)象，所以需采用模擬退火算法進(jìn)行動態(tài)航線優(yōu)化．

在人工勢場法航線中加入模擬退火算法，可以通過增加隨機(jī)擾動值調(diào)整新航路點(diǎn)，避免陷入局部最小值，還可以搜索大風(fēng)浪區(qū)附近的航路點(diǎn)位置，替代原有的負(fù)梯度下降受力控制的方法，避免人工勢場法中參數(shù)系數(shù)對航路點(diǎn)選取的影響，可有效緩解抖動現(xiàn)象．

每次經(jīng)過調(diào)整之后的航路點(diǎn)組合便構(gòu)成了一條新航線，可以利用目標(biāo)函數(shù)值比較以及劣航線的接受準(zhǔn)則進(jìn)行處理．將新航線的目標(biāo)函數(shù)值設(shè)為f(A)，當(dāng)前最優(yōu)航線目標(biāo)函數(shù)值設(shè)為f(B)．若f(A)

Metropolis準(zhǔn)則判斷具體方法為

設(shè)當(dāng)前航線為xi，新航線為xj，其目標(biāo)函數(shù)值分別記為f(xi)和f(xj)，新航線被接受可以作為最優(yōu)航線的概率Pt為

(14)

在判斷過程中，系統(tǒng)將產(chǎn)生一個隨機(jī)數(shù)ε∈(0,1)，并與Pt進(jìn)行比較，若Pt>ε，則將新航線設(shè)為當(dāng)前最優(yōu)航線；否則，舍棄新航線，繼續(xù)使用當(dāng)前航線作為最優(yōu)航線．從式(14)可見，溫度越低，接受概率Pt的值也就越小，即劣航線被接受的概率會越來越小，當(dāng)退火到最小值時，劣航線基本被全部淘汰．

經(jīng)過多次實(shí)驗(yàn)，選用了以下參數(shù)：初始溫度t0=100；衰減系數(shù)α=0.7；終止溫度tf=10；迭代次數(shù)L=500．

模擬退火算法的具體流程為

步驟1初始參數(shù)設(shè)置 初始溫度t0，衰減系數(shù)α，終止溫度tf，迭代總次數(shù)L．

步驟2初始值設(shè)置 初始航線x0，當(dāng)前航線xcurrent，新航線為xnew，最優(yōu)航線為xbest，初始目標(biāo)值設(shè)為E0，當(dāng)前航線航程為Ecurrent，新航線航程為Enew，最優(yōu)航線航程為Ebest，當(dāng)前迭代次數(shù)l=0．

步驟3調(diào)整航路點(diǎn) 生成隨機(jī)數(shù)n和n個2維擾動變量ΔLat,ΔLon，隨機(jī)選擇當(dāng)前航線上n個航路點(diǎn)(除始點(diǎn)和終點(diǎn))，令Lat′=Lat+ΔLat,Lon′=Lon+ΔLon,得到新航線xnew，計算其目標(biāo)值Enew，迭代次數(shù)l=l+1．

步驟4計算新航線xnew，與最優(yōu)航線xbest之間的航程差值ΔE，ΔE=Enew-Ecurrent，如果ΔE<0，那么令xbest=xcurrent=xnew，Ebest=Ecurrent=Enew，轉(zhuǎn)至步驟6．否則轉(zhuǎn)至步驟5．

步驟5令ΔE=Enew-Ecurrent，若exp(-ΔE/ti)>ε，則xcurrent=xnew，Ecurrent=Enew，其中ε為(0，1)內(nèi)的隨機(jī)數(shù)；否則xcurrent和Ecurrent的值不變．

步驟6若l≥L，令ti+1=αti、l=0，轉(zhuǎn)至步驟7；否則返回步驟3．

步驟7若ti+1≤tf，算法終止；否則返回步驟3．

3 仿真結(jié)果及分析

3.1 氣象資料獲取和處理

研究采用歐洲中期天氣預(yù)報中心(european center for medium-range weather forecasts，ECMWF)提供的氣象數(shù)據(jù)．ECMWF主要提供中長期天氣預(yù)報如實(shí)時的天氣、海況信息來為遠(yuǎn)洋航行船舶提供服務(wù)．由于本文在模型中主要考慮風(fēng)浪對船舶的影響，所以采用ECMWF提供海平面10 m高處的風(fēng)速與風(fēng)向數(shù)據(jù)．氣象數(shù)據(jù)文件格式采用網(wǎng)絡(luò)通用數(shù)據(jù)格式NetCDF(Network Common Data Format)，數(shù)據(jù)更新間隔為6 h，數(shù)據(jù)精度的網(wǎng)格大小為1°×1°．利用Matlab可以將氣象數(shù)據(jù)可視化，更直觀顯示風(fēng)浪分布情況[11]．圖2為2019年6月太平洋區(qū)域第8日12:00時刻的風(fēng)速數(shù)據(jù)信息．圖3為2019年6月太平洋區(qū)域第8日中12:00時刻的浪高數(shù)據(jù)信息．

圖2 風(fēng)速與風(fēng)向數(shù)據(jù)信息

圖3 浪高數(shù)據(jù)信息

3.2 結(jié)果及分析

選用大海域太平洋作為航行水域，起點(diǎn)為日本橫濱港(35°N,141°E)，終點(diǎn)為美國舊金山港(37°N,123°W)．仿真選用的船舶為標(biāo)準(zhǔn)排水量23 740 t的集裝箱船S175，船舶靜水速度設(shè)為15 kn，浪高與風(fēng)級合并考慮，僅將船舶與風(fēng)浪的遭遇角作為海況考慮[12]．選用2019年5月30日和6月8日的大風(fēng)浪數(shù)據(jù)在仿真中展示效果圖(見圖4)，航線為大圓航線．

圖4 航行途徑大風(fēng)浪區(qū)示意圖

大圓航線為球面上兩點(diǎn)之間距離最短的航線，即在理想環(huán)境下船舶走大圓航線將是最優(yōu)選擇，但在現(xiàn)實(shí)中船舶并不容易以大圓航線航行，如考慮危險氣象環(huán)境時，航線必將有所調(diào)整．因此本文在進(jìn)行航線規(guī)劃時將以大圓航線為基準(zhǔn)并利用人工勢場算法進(jìn)行動態(tài)調(diào)整，在經(jīng)過人工勢場法避障之后的航線其航程里數(shù)并不是最優(yōu)解，因此再利用模擬退火算法對該航線進(jìn)行優(yōu)化，優(yōu)化結(jié)果見圖5，幾種航線的航程與航時對比見表1．

圖5 優(yōu)化航線圖

表1 三種航線結(jié)果對比

由圖5和表1可知：單純的人工勢場算法所規(guī)劃航線可以很好的避開大風(fēng)浪區(qū)域，且相較大圓航線船舶平均航速有所提高，但航行距離和航行時間有明顯增加．在人工-模擬法航線中則有較好的改進(jìn)，雖總航程相較大圓航線有所增加，但總航時明顯降低而且平均航速也有明顯提高．

4 結(jié) 束 語

文中提出了一種結(jié)合人工勢場法和模擬退火算法的氣象航線優(yōu)化方法，并用Matlab進(jìn)行模擬仿真，在考慮風(fēng)速和浪高氣象因素下，調(diào)整風(fēng)浪中各航段航速，優(yōu)化縮短航程和航時．該方法能根據(jù)船舶自身條件設(shè)置安全距離，有效避過大風(fēng)浪區(qū)，人工-模擬退火算法對航線優(yōu)化效果良好，可以顯著提高船舶平均航速，縮短航行時間．

文中的氣象因素僅考慮了風(fēng)速和浪高，為更準(zhǔn)確模擬海上環(huán)境，后續(xù)研究中還應(yīng)考慮云霧、暴雨等更多海上氣象要素的影響，并從經(jīng)濟(jì)層面考慮，構(gòu)建如油耗等模型以進(jìn)一步完善優(yōu)化算法．