閔德權(quán)，江可鑒，劉 茹

(大連海事大學(xué) 交通運(yùn)輸工程學(xué)院，遼寧 大連 116026)

0 引 言

隨著全球氣候持續(xù)變暖，北極海冰的覆蓋面積逐漸縮小，北極航線的通航問題得到了更多的關(guān)注與研究?？傮w來看，北極航線可分成東北航道、西北航道以及穿極航線。北極航線東北航道的通航期從7月中下旬開始，在10月中下旬結(jié)束[1]。但相較傳統(tǒng)航線，其具有航線時間和航線成本上的優(yōu)勢。

針對船舶的航速優(yōu)化，需先對航線進(jìn)行總體分析。ZHENG Binglei等[2]通過對北極航線上船舶航行的航速、干擾以及延遲等因素進(jìn)行分析，建立以航線運(yùn)營成本最小為目標(biāo)的優(yōu)化模型；關(guān)曉光等[3]在對北極航線的東北航線進(jìn)行航次的收益和效率分析的基礎(chǔ)上，建立環(huán)繞模式集裝箱航線優(yōu)化模型；J.H.NAM等[4]基于海冰和環(huán)境等數(shù)據(jù)，建立交互式的仿真系統(tǒng)，從而求解出最佳的北極航線。在上述研究的基礎(chǔ)上，可對船舶的航速做進(jìn)一步分析。A.AFONIN等[5]利用包含的不同船舶有關(guān)信息的地理信息系統(tǒng)(GIS)，并結(jié)合冰情數(shù)據(jù)，統(tǒng)計處理冰區(qū)航行船舶的航速以識別航速減少的趨勢；FU Shanshan等[6]建立概率模型以預(yù)測船舶在北極水域的困阻概率，并使用貝葉斯信念網(wǎng)絡(luò)(bayesian belief networks)模型預(yù)測船舶的航速；李振福等[7]基于冰區(qū)氣象的條件，分析海況對船舶的干擾力和冰阻力等因素，建立以不同冰情航速優(yōu)化模型。

基于上述研究，筆者以北極東北航線為基礎(chǔ)，建立通航期冰情條件下的混合整數(shù)雙目標(biāo)船舶航速優(yōu)化模型。使用線性加權(quán)和離散化的方法對模型進(jìn)行處理，選擇MATLAB軟件對模型進(jìn)行求解，以期得出不同船舶的航速優(yōu)化方案。

1 航線通航期冰情分析

北極東北航線的可通航時間主要集中于7月到10月[8]，將上述4個月作為通航期，并以通航期的冰情數(shù)據(jù)為基礎(chǔ)對船舶航速進(jìn)行優(yōu)化分析。航線的冰情包括海冰厚度、海冰密集度等數(shù)據(jù)，目前可使用的數(shù)據(jù)庫有：美國國家冰雪數(shù)據(jù)中心(National Snow & Ice Data Center)數(shù)據(jù)庫、德國不來梅大學(xué)(The University of Bremen)的海冰數(shù)據(jù)庫等。

選擇德國不來梅大學(xué)海冰數(shù)據(jù)庫中2020年7月至10月數(shù)據(jù)，對航線經(jīng)過海域的海冰密集度和海冰厚度數(shù)據(jù)進(jìn)行分析。

1.1 通航期海冰密集度數(shù)據(jù)分析

不來梅大學(xué)海冰數(shù)據(jù)庫所提供的海冰密集度數(shù)據(jù)以HDF、netCDF等數(shù)據(jù)格式進(jìn)行存儲[9]，可使用相應(yīng)軟件進(jìn)行讀取。

在航線的涉及海域中，渤海在冬季部分海域開始結(jié)冰，到次年的3月海冰完全消失[10]。巴倫支海西南部存在常年無冰區(qū)，且在夏季期間出現(xiàn)無冰情況[11]。挪威海和北海在夏季可通航，亦不做冰情數(shù)據(jù)分析。在現(xiàn)有海冰密集度數(shù)據(jù)中，白令海7至10月海冰密集度均值在1%以下，其海冰密集度僅作為計算參數(shù)進(jìn)行使用。因此設(shè)定楚科奇海、東西伯利亞海、拉普捷夫海以及喀拉海作為冰區(qū)海域進(jìn)行分析，其通航期月均海冰密集度數(shù)據(jù)如圖1。

從圖1中海域的海冰密集度數(shù)據(jù)可以看出，楚科奇海在2020年的4個月的月均冰密集度呈現(xiàn)出較高態(tài)勢，最高冰密集度接近60%。拉普捷夫海和喀拉海的月均冰密集度較低，均在10%以下。

圖1 冰區(qū)海域的海冰密集度

從總體來看，上述海域的冰密集度數(shù)據(jù)在8月或9月較低，但在7月或10月的波動較為劇烈。除喀拉海以外，其他3個海域在通航期4個月的變化趨勢基本一致。

1.2 通航期冰厚度數(shù)據(jù)分析

使用不來梅大學(xué)海冰數(shù)據(jù)庫的薄海冰厚度(the thickness of thin sea ice, SIT)庫所記錄的北極海冰數(shù)據(jù)[12]，其中部分?jǐn)?shù)據(jù)存儲于內(nèi)部資料庫。

選擇2020年7至10月的北極海冰月均厚度數(shù)據(jù)作為后期模型求解的參考數(shù)據(jù)，不再做具體海域的數(shù)值分析，其數(shù)值如表1。

表1 海冰的月均厚度

2 船舶阻力模型構(gòu)建

在冰區(qū)航行的船舶不僅要考慮水的作用效應(yīng)，還要考慮海冰對船舶的影響。因此，分別對船舶的靜水阻力和冰阻力進(jìn)行模型構(gòu)建。對阻力模型進(jìn)行以下設(shè)定：①對靜水阻力和冰阻力區(qū)分分析；②不考慮船舶附屬體帶來的附體阻力。

2.1 靜水阻力模型

當(dāng)船舶在無冰的靜水中航行時，其所受的阻力(即裸船體阻力)f0的計算如式(1)[13]：

(1)

式中：C0為裸船體阻力的阻力系數(shù)；ρ為海水密度；v為船舶的航行速度；S為船舶的濕表面積。

船舶的濕表面積計算如式(2)[14]:

(2)

式中：Lw為船舶的水線長度；d為船舶的吃水；B為船舶的型寬；Cb為方形系數(shù)。

依據(jù)文獻(xiàn)[13]，船舶在無冰的靜水中航行所受到的阻力f0可分為摩擦阻力f1、粘壓阻力f2以及興波阻力f3：

f0=f1+f2+f3

(3)

考慮到船體表面的粗糙度，計算船舶摩擦阻力時需額外加入粗糙度補(bǔ)貼系數(shù)CΔ[13]。

CΔ可由船長L進(jìn)行選取，其具體值如表2[7]，設(shè)定表中船長不計邊界值。

表2 粗糙度補(bǔ)貼系數(shù)選取

2.1.1 摩擦阻力

摩擦阻力計算如式(4)：

(4)

依據(jù)文獻(xiàn)[13]，船舶的摩擦阻力系數(shù)可以使用1957年國際船模實驗池會議所提出的計算公式：

(5)

式中：Re為雷諾數(shù)；γ為海水的運(yùn)動粘性系數(shù)。

(6)

2.1.2 粘壓阻力

依據(jù)文獻(xiàn)[13]，粘壓阻力系數(shù)可使用巴甫米爾提出的近似公式進(jìn)行計算。粘壓阻力f2和粘壓阻力系數(shù)C2的計算如式(7)、式(8)：

(7)

(8)

式中：Am為船舶中橫剖面的水線下面積，后簡稱為中橫剖面積。Am可通過式(9)進(jìn)行計算，Cm為中橫剖面系數(shù)。

Am=CmBd

(9)

Lr為船舶的去流段長度，為了避免嚴(yán)重的船舶的航行漩渦效應(yīng)，選取去流段長度為[15]：

(10)

2.1.3 興波阻力

通過對Michell積分進(jìn)行簡化，可對興波阻力進(jìn)行計算[16]。具體公式如式(11)：

(11)

式中：λ=secθ，θ為凱爾文波系(Kelvin wave system)相對于x軸的角；ω1與ω2與船舶表面函數(shù)的縱向梯度有關(guān)。

簡化后的興波阻力公式也難以直接計算出船舶受到的興波阻力。而文獻(xiàn)[17]提供的測試程序可進(jìn)行簡化測算，該程序通過計算Michell積分來估計Wigley船型的波阻力，且提供的阻力數(shù)據(jù)可作為計算興波阻力的基礎(chǔ)。

首先，在上述測試程序中輸入相應(yīng)的船舶參數(shù)，運(yùn)行程序得到多航速對應(yīng)的阻力數(shù)據(jù)。其次，對結(jié)果數(shù)據(jù)進(jìn)行MATLAB的傅里葉級數(shù)擬合，得到多參數(shù)的興波阻力和船舶航速的擬合函數(shù)，其擬合函數(shù)如式(12)：

(12)

式中：a0、aj、bj以及z為對應(yīng)的擬合參數(shù)。

2.2 冰阻力模型

設(shè)定船舶在通航期航行所受到冰阻力為浮冰(pack ice)阻力fe，可由式(13)進(jìn)行計算[18]。

fe=CeρegBhv2Wn

(13)

式中：ρe為海冰密度；h為海冰厚度；W為海冰密集度；Ce浮冰阻力系數(shù)；n為海冰密集度的冪次數(shù)；g為重力加速度。依據(jù)文獻(xiàn)[18]，使用n=2對于計算具有極小的誤差。

依據(jù)文獻(xiàn)[18]，加拿大海洋技術(shù)研究所經(jīng)過多年試驗得出lnCe和lnFr呈線性關(guān)系，Ce和Fr的關(guān)系如式(14)。

Ce=4.4Fr-0.8267

(14)

式中：Fr為浮冰的弗勞德數(shù)，計算如式(15)：

(15)

3 航速優(yōu)化模型構(gòu)建

3.1 問題概述

在傳統(tǒng)航線的運(yùn)營中，船舶通常有固定的掛靠港口[19]，且在港口間的貨運(yùn)量相對確定。船公司可以根據(jù)自身的經(jīng)營狀況，對航線上運(yùn)營的多種船舶進(jìn)行合理調(diào)配，以期得到最大的經(jīng)濟(jì)效益。

在北極航線中，船舶一般只在通航期進(jìn)行航行，難以在航線內(nèi)進(jìn)行多航段和多船組合的調(diào)配。因此，在現(xiàn)有船舶條件下，建立北極航線單船單航次的雙目標(biāo)航速優(yōu)化模型，在船舶的港口使費(fèi)、固定成本及燃油成本在內(nèi)的運(yùn)營成本最低，以及船舶的航行時間和停泊時間在內(nèi)的運(yùn)營時間達(dá)到最低的條件下，并求解船舶的最佳航速。依據(jù)第2節(jié)中的船舶阻力公式，可對船舶克服阻力進(jìn)行航行的做功進(jìn)行轉(zhuǎn)化求出船舶的燃油成本，進(jìn)而確定船舶運(yùn)營成本。同時對航線進(jìn)行航段劃分，使不同船舶的航速在不同航段達(dá)到最優(yōu)。對船舶的航行過程做出3點(diǎn)假設(shè)：①航線上僅配置單個船舶，且船舶在航線的每個航段勻速航行；②不考慮航線極端氣象狀況，選取理想航行狀態(tài)。③不使用額外破冰船進(jìn)行輔助航行。

3.2 規(guī)劃模型

第1個目標(biāo)函數(shù)如式(16)，表示船舶完成航線航行的運(yùn)營成本最小化。式(16)由船舶的港口使費(fèi)、船舶航行的固定成本、船舶的燃油成本3個部分構(gòu)成。

(16)

第2個目標(biāo)函數(shù)如式(17)，表示船舶完成航線航行的運(yùn)營時間最小化。式(17)由船舶航行時間、船舶停泊時間兩部分組成。

(17)

(18)

使用燃油的燃燒熱值對船舶做功進(jìn)行燃油量的換算，從而求出燃油成本。船舶的燃機(jī)和傳動機(jī)構(gòu)復(fù)雜，難以計算其實際效率，因此不計船舶的能量損耗。

(19)

式中：DL為輕質(zhì)燃油單位價格；DH為重質(zhì)燃油單位價格；ε1為輕質(zhì)燃油的燃燒熱值；ε2為重質(zhì)燃油的燃燒熱值；vki為船舶k在i航段的航速。

(20)

船舶于不同航段的航行時間如式(21)，分為3項：第1項為冰區(qū)航段的船舶航行時間；第2項為ECA航段的船舶航行時間；第3項為普通航段的船舶航行時間。

(21)

(22)

航段的航速限制范圍如式(23)：

(23)

式(24)表示航線只能使用單個船舶k。

(24)

4 算例的求解與分析

4.1 算法設(shè)計

建立的航速優(yōu)化模型為雙目標(biāo)模型，兩個目標(biāo)函數(shù)均為線性函數(shù)，且存在離散化的參數(shù)條件。求解上述模型可使用線性加權(quán)法或啟發(fā)式算法。啟發(fā)式算法的求解速度較快，但求解結(jié)果通常為近似最優(yōu)解。因此，針對雙目標(biāo)航速優(yōu)化模型，使用線性加權(quán)法進(jìn)行離散化求解。求解軟件為MATLAB，求解步驟為：

1)使用線性加權(quán)法進(jìn)行重要度賦權(quán)，將多目標(biāo)問題轉(zhuǎn)化為單目標(biāo)混合整數(shù)規(guī)劃問題。設(shè)置相應(yīng)權(quán)值?1和?2，構(gòu)造新目標(biāo)函數(shù)Y3：

minY3=?1Y1+?2Y2

(25)

2)根據(jù)多船舶和多航段條件進(jìn)行目標(biāo)問題拆分，使決策變量xk失效，目標(biāo)問題從而離散化為若干個一般的線性規(guī)劃子問題。

3)使用MATLAB求解上述若干子問題，得出對應(yīng)的最優(yōu)解。

4.2 船舶信息

選擇4艘冰級運(yùn)輸船作為算例，船舶具體類型分別為：雜貨運(yùn)輸船、集裝箱運(yùn)輸船、原油運(yùn)輸船、液化天然氣運(yùn)輸船。上述船舶基于俄羅斯船級社標(biāo)準(zhǔn)的冰級分別為：Arc 4、Arc 4、Arc 7、Arc 7。船舶具體數(shù)據(jù)如表3。

表3 船舶數(shù)據(jù)

基于現(xiàn)有海冰密集度和海冰厚度數(shù)據(jù)，設(shè)定中度冰情(冰密集度在30%～60%和冰厚度在30 cm內(nèi))時船舶應(yīng)按照預(yù)定航速的80%進(jìn)行航行[4]。有冰情數(shù)據(jù)但低于上述值時，設(shè)為輕度冰情。船舶航速的80%為冰區(qū)航段中度冰情時的安全航速范圍。

4.3 航線的航段劃分

算例航線選擇天津港為起始掛靠港，鹿特丹港為終止掛靠港。并以符拉迪沃斯托克、彼得羅巴甫洛夫斯克、普羅維杰尼亞、佩韋克、季克西、迪克森以及摩爾曼斯克的港口為基準(zhǔn)港進(jìn)行航線長度的測算。設(shè)定船舶物資儲備充足，可滿足船舶完成航行的需求，船舶不掛靠基準(zhǔn)港。設(shè)定船舶在起始港裝載貨物，在終止港卸載貨物。整條航線按照航行海域可劃分11個航段，航線內(nèi)有1個起始港、1個終止港、7個途經(jīng)港。

對航線的航段劃分做出以下設(shè)定：①對部分基準(zhǔn)港區(qū)間內(nèi)的非冰區(qū)海域進(jìn)行整體航段劃分;②可收集到冰情數(shù)據(jù)(海冰密集度、海冰厚度)的海域跨越多基準(zhǔn)港區(qū)間時：基于離散化的方法，對上述涉及海域的航段進(jìn)行分屬基準(zhǔn)港區(qū)間的子航段劃分。同屬一個海域的子航段的冰情參數(shù)和海水參數(shù)取值一致。為避免基準(zhǔn)港區(qū)間的干擾，選取同海域的第一個子航段的航速作為航速分析的基礎(chǔ)。航線的總體航段劃分如表4。

表4 航線的總體航段劃分

4.4 優(yōu)化結(jié)果分析

針對船舶的燃油成本，重質(zhì)燃油的參考價格為430美元/t，輕質(zhì)燃油的參考價格為530美元/t。在起始港，選取包括船舶停泊和引航(移泊)費(fèi)用、貨物的港務(wù)費(fèi)在內(nèi)的港口使費(fèi)，終止港視具體規(guī)定進(jìn)行選取。依據(jù)既定條件對模型進(jìn)行求解，可得出算例船舶在多航段中的不同月的最優(yōu)航速，航速單位為kn。

表5為4艘船舶航速的優(yōu)化結(jié)果。從表5可以看出，在非冰區(qū)航段，4艘船舶的航速波動較小，其航速均值分別為14.88、14.46、13.39、16.12 kn。而在冰區(qū)航段，相較于非冰區(qū)航段其航速有所下降，船舶在中度冰情航段航速均值分別為：12.72、12.42、11.86、13.91 kn。

表5 船舶的最優(yōu)航速

圖2為4艘船舶在非冰區(qū)航段和冰區(qū)航段通航期不同月份的航速對比，冰區(qū)航段出現(xiàn)中度冰情顯示有降幅航速的均值。圖3為航速優(yōu)化后的單航次船舶運(yùn)營成本和運(yùn)營時間。

圖2 船舶航速對比

由圖2可知：在非冰區(qū)航段，船舶的航速無明顯變化。而在冰區(qū)的中度冰情航段，船舶的航速有一定的降幅。依據(jù)表5的數(shù)據(jù)，在中度冰情航段，4艘船舶的航速均值相較非冰區(qū)航段的降幅最大達(dá)到14.53%、14.12%、11.41%、13.67%。由圖3可知，通航期冰情變化對船舶的運(yùn)營成本以及運(yùn)營時間的影響較小。因此，筆者所構(gòu)建的模型能夠較好的反映船舶在北極航線通航期的實際航行狀態(tài)。

圖3 單航次運(yùn)營成本和運(yùn)營時間

5 結(jié) 論

筆者對北極東北航線的通航期冰情數(shù)據(jù)進(jìn)行分析，考慮船舶的阻力因素，建立以船舶單航次的運(yùn)營成本和運(yùn)營時間最低為目標(biāo)的航速優(yōu)化模型。其次，使用線性加權(quán)和離散化的方法處理所建立的模型。最后，使用MATLAB軟件對模型進(jìn)行求解，得出多種船舶不同航段的最優(yōu)航速。從航速結(jié)果可以看出，在非冰區(qū)航段的船舶航速變化較小，而在中度冰情的條件下船舶航速有一定的降幅。相較非冰區(qū)航段，4艘算例船舶在中度冰情航段的航速均值最大降幅達(dá)到14.53%、14.12%、11.41%、13.67%。