張 赫, 閆建鑫, 邢江豪

(大連海事大學(xué)交通運輸工程學(xué)院， 大連 116026)

隨著世界經(jīng)濟和技術(shù)的快速發(fā)展，各種類型船舶的數(shù)量隨之持續(xù)增加，船舶之間交通沖突和擁擠的現(xiàn)象愈發(fā)嚴(yán)重。例如海上船舶在行駛到工程船舶的作業(yè)區(qū)域時，大量的工程船舶根據(jù)作業(yè)地點的不同需要多次穿越航道，這與往來正常航行的其他船舶形成了一個復(fù)雜的航線網(wǎng)絡(luò)，船舶的航行受到極大的影響。在這種情況下，需要對船舶的行駛路線進(jìn)行合理規(guī)劃，而誘導(dǎo)技術(shù)是實現(xiàn)路徑優(yōu)化的有效方法之一。

路徑誘導(dǎo)技術(shù)自20世紀(jì)產(chǎn)生到現(xiàn)在已經(jīng)幾十年的歷史，相關(guān)的資料也比較完善，目前更多的應(yīng)用于車輛行駛的誘導(dǎo)。對于汽車的路徑誘導(dǎo)中外相關(guān)領(lǐng)域的專家采用不同方式和方法進(jìn)行了大量研究[1-4]，而針對于船舶的誘導(dǎo)研究較少，僅有喬春福[5]研究了如何根據(jù)船位序列判斷船舶上行和下行的算法，從而為過往船舶的上下行助航提供誘導(dǎo)幫助。從文獻(xiàn)中可以看出，目前對船舶的航線誘導(dǎo)還在發(fā)展之中，僅有的研究也只是考慮了行駛水域的特點結(jié)合地圖進(jìn)行上下行的助航，沒有針對于某一船舶在復(fù)雜的水域內(nèi)行駛的航線誘導(dǎo)。而陸上汽車的路徑誘導(dǎo)問題經(jīng)過多年的發(fā)展已經(jīng)趨于成熟，相關(guān)理論也相對完善，這對實現(xiàn)船舶的航線誘導(dǎo)具有指導(dǎo)意義?，F(xiàn)在考慮與陸上汽車誘導(dǎo)所不同的氣象水文條件影響的基礎(chǔ)上結(jié)合船舶航行的特點將誘導(dǎo)技術(shù)運用到船舶的航線規(guī)劃之中，以期為船舶航行提供科學(xué)依據(jù)。

1 海況對船舶航行的影響

船舶在海上航行時，空氣和水對船舶行駛具有很大的影響，為了便于研究船舶在不同的風(fēng)浪條件下的行駛速度，將船舶阻力分為基本阻力、空氣阻力、洶濤阻力三部分。

1.1 基本阻力

船舶基本阻力是船舶在靜水中航行受到的阻力，由摩擦阻力和壓阻力(剩余阻力)組成。結(jié)合傅汝德相似定律和雷諾相似定律可得船舶靜水行駛基本阻力R0的計算表達(dá)式[6]：

(1)

CR=Cf+ΔCf+Cr

(2)

式中：CR為基本阻力系數(shù)，ρ0為流體密度，V0為船舶在靜水中的行駛速度；S為船舶濕面積；Cf為船舶摩擦阻力系數(shù)；ΔCf為粗糙度補貼系數(shù)；Cr為船舶剩余阻力系數(shù)。

摩擦阻力系數(shù)Cf用“1957年國際船模實驗池實船-船模換算公式”[式(3)]進(jìn)行求解, 粗糙度補貼系數(shù)ΔCf如表1所示；剩余阻力系數(shù)Cr利用Lap-kelller圖譜進(jìn)行查表計算獲得[7]。

(3)

(4)

式中：Re為雷諾數(shù)；V為船速；Lwl為水線長，近似于船長L；ν為水運動黏性系數(shù)。

表1 粗糙度補貼系數(shù)Table 1 Roughness subsidy coefficients

1.2 空氣阻力

船舶在海上航行時，其水面以上的部分在風(fēng)的作用下會產(chǎn)生空氣阻力影響行駛速度，產(chǎn)生船舶艏艉方向的風(fēng)壓力Fw可按式(5)計算[8]：

(5)

(6)

(7)

(8)

lgY=αw+βwlgX

(9)

式中：ρa為空氣密度(1.226 kg/m3)；Cw為艏艉方向上風(fēng)力系數(shù)；Va為相對風(fēng)速；Aa為水線以上船體正面受風(fēng)面積；Ba為水線以上船體側(cè)面受風(fēng)面積；Vzf為真風(fēng)風(fēng)速；αa為真風(fēng)方向角，左舷為正；θ為風(fēng)舷角；Ca為風(fēng)力系數(shù)；α為風(fēng)力角，是風(fēng)壓力與船舶艏艉線之間的夾角；Y為船舶受風(fēng)面積(包括Aa、Ba)；X表示船舶載重量或總質(zhì)量，t；αw、βw為回歸系數(shù)，其數(shù)值如表2所示[9]。

表2 不同船型船舶的受風(fēng)面積回歸系數(shù)Table 2 Regression coefficient of windage area of different types of ships

風(fēng)力角α和風(fēng)力系數(shù)Ca可根據(jù)巖井聰公式[式(10)和式(11)]進(jìn)行計算[6-8]：

(10)

Ca=1.325-0.05cos(2α)-0.35cos(4α)-0.175cos(6α)

(11)

1.3 洶濤阻力

洶濤阻力分為規(guī)則波干擾和不規(guī)則波干擾兩種，通常只考慮規(guī)則波的影響。對于規(guī)則波干擾力的計算比較復(fù)雜，目前主要依靠船模實驗進(jìn)行，Daidola在船模實驗的基礎(chǔ)上提出了洶濤阻力的計算模型為[6-7]

(12)

(13)

1.4 船舶自由行駛速度

假定船舶在靜水中行駛的推進(jìn)功率和風(fēng)浪中行駛的推進(jìn)功率保持不變，則船舶在風(fēng)浪中行駛減少的速度用來抵消風(fēng)浪阻力的影響[6]，從而得到船舶在風(fēng)浪下自由航行速度Vf。

(14)

由于船舶在不同路段行駛的位置和方向并不一致這就導(dǎo)致相對風(fēng)速的大小、方向發(fā)生變化，因此船舶在不同路段行駛時自由航行速度Vf需要根據(jù)行駛方向和位置分別求解。

2 路段行程時間計算

在對船舶進(jìn)行航線誘導(dǎo)時，優(yōu)化目標(biāo)通常有距離、時間和費用三種，通過綜合考慮，現(xiàn)以航線用時最短作為優(yōu)化目標(biāo)。

各路段行程時間分為路段行駛時間和航線網(wǎng)絡(luò)交叉點延誤時間兩個部分。在沒有航線交叉的路段內(nèi)，路段行程時間就等于路段行駛時間。船舶在路段的行駛時間與路段船舶流量密切相關(guān)，根據(jù)風(fēng)浪對船舶航行的影響對美國聯(lián)邦公路局函數(shù)即BPR函數(shù)進(jìn)行改進(jìn)，建立路段行駛時間的表達(dá)式：

(15)

式(15)中:tx為兩節(jié)點之間路段行駛時間；ψ為風(fēng)浪對航速的影響系數(shù)，ψ=Vf/V0；t0為兩節(jié)點之間交通量為零時無風(fēng)浪影響下的路段行駛時間；q為路段船舶交通量；c為航路實際通行能力；γ1、γ2為相關(guān)系數(shù)。

航線交叉點平均延誤時間td根據(jù)排隊論原理進(jìn)行求解。將船舶看作“顧客”，交叉點視為“服務(wù)窗”，船舶在交叉點的排隊穿越過程可視為一個排隊問題。船舶到達(dá)交叉點時，若“服務(wù)窗”被占用，則進(jìn)行排隊等待，最大排隊長度n為航道長度和單位船舶所需的航道長度Ld的比值。對船舶到達(dá)規(guī)律進(jìn)行統(tǒng)計分析，船舶的到達(dá)規(guī)律可近似服從參數(shù)為λ的泊松分布，船舶通過交叉點的服務(wù)時間可看作服從參數(shù)為μ的負(fù)指數(shù)分布(到達(dá)率λ=q,μ可統(tǒng)計獲得)[10]。因此交叉點的船舶排隊模型為M/M/1/n/∞模型。

(16)

綜上，ij路段總行程時間tij為

tij=tx+td

(17)

3 模型構(gòu)建

3.1 影響因素分析

3.1.1 施工水域影響

在海上航行的船舶可能會遇到海上施工區(qū)域，影響船舶航行。由于施工區(qū)域內(nèi)有大量工程船舶進(jìn)行運輸作業(yè)，使得施工區(qū)內(nèi)航線密集，交叉增多。

3.1.2 海風(fēng)動態(tài)變化

船舶位置不同，船舶受到風(fēng)力大小也不一致，為此建立船舶風(fēng)速與位置的對應(yīng)關(guān)系公式Vzf=φ1Vfe，其中φ1為與位置相關(guān)的風(fēng)力作用系數(shù)(0≤φ1≤1)，Vfe為基本風(fēng)速。

3.1.3 波浪動態(tài)變化

波浪振幅的變化與海風(fēng)息息相關(guān)，并且影響洶濤阻力的大小。為此建立波浪振幅與位置的對應(yīng)關(guān)系公式Wa=φ2W0，其中φ2為與位置相關(guān)的波浪作用系數(shù)(0≤φ2≤1)，W0為基本波浪振幅。

3.1.4 沉點的影響

對船舶航行產(chǎn)生阻礙的物體所在的位置稱為沉點。在施工區(qū)域內(nèi)，進(jìn)行挖泥、疏浚作業(yè)的船舶通常移動速度較慢可視為沉點，若沉點存在于路段中相應(yīng)路阻為無窮大。淺灘、島嶼等存在的位置也是沉點，為禁止行駛區(qū)域需繞開行駛。

3.2 模型建立

在對船舶航線誘導(dǎo)時，將航線網(wǎng)絡(luò)抽象成圖論網(wǎng)絡(luò)G=(P,E,T)，其中P是節(jié)點集合表示航道節(jié)點；E是弧的集合表示路段；T={tij(Ti)|(i,j)∈E}是各弧的權(quán)值，表示在此路段的行程時間tij(Ti)為船舶在Ti時刻從節(jié)點i航行到節(jié)點j所需花費的時間。設(shè)施工區(qū)k(k=1,2,…，l)的節(jié)點集合為Pk存在的時間為[Tks,Tke]，從起點到達(dá)施工區(qū)k的時間為tk，非施工區(qū)節(jié)點結(jié)合為P0。則建立船舶動態(tài)時變網(wǎng)絡(luò)下的船舶航線誘導(dǎo)優(yōu)化模型為

(18)

s.t.

(19)

(20)

(21)

(22)

(23)

k=1,2,…,l

(24)

(25)

Vzf=φ1Vfe

(26)

Wa=φ2W0

(27)

式(18)為優(yōu)化目標(biāo)；式(19)表示路段ij在時刻被路徑使用則等于1，否則等于0；式(20)表示初始點只能是出發(fā)點不能是到達(dá)點，終點只能是到達(dá)點不能是出發(fā)點，中間節(jié)點的到達(dá)和出發(fā)次數(shù)相等；式(21)表示所有路段最多只能被選擇一次；式(22)為施工區(qū)內(nèi)沉點對行程時間的影響；式(23)、式(24)表示船舶行駛到施工區(qū)域位置時若在施工區(qū)域存在的時間范圍內(nèi)則將區(qū)域內(nèi)節(jié)點加入到航線網(wǎng)絡(luò)之中，否則保留原節(jié)點；式(25)表示船舶各路段自由行駛速度的計算；式(26)為風(fēng)速與位置的對應(yīng)關(guān)系式；式(27)為波浪振幅與位置的對應(yīng)關(guān)系式；行程時間計算過程中的其他參數(shù)的計算公式約束第1節(jié)與第2節(jié)已經(jīng)說明，此處不再列出。

對于各路段船舶流量進(jìn)行實時的更新是非常困難的，通常將時間根據(jù)實際劃分為多個時段，并假定在同一個時段航路網(wǎng)各路段的船舶流量不發(fā)生改變。

4 蟻群算法求解

(28)

螞蟻釋放在各節(jié)點之間的信息素隨著時間的流逝也會逐漸消失，設(shè)參數(shù)ε3(0<ε3<1)代表信息素的揮發(fā)系數(shù)。當(dāng)m只螞蟻都完成了一次循環(huán)后，需要對各節(jié)點間t時刻的τij(t)進(jìn)行迭代更新，公式為

(29)

(30)

式(30)中：Q為螞蟻在一次循環(huán)中所散發(fā)的信息素總量；Lk為第k只螞蟻在本次循環(huán)中所花費的總時間。

針對于本文航線誘導(dǎo)問題，以船舶出發(fā)點為螞蟻路徑起點o，以船舶目標(biāo)點為螞蟻路徑尋找的終點，以螞蟻選擇的節(jié)點集合表示船舶路徑的可行解。假設(shè)初始時刻T=T0,則用蟻群算法求解的基本過程如下：

(1)對螞蟻總數(shù)目m、信息素影響因子ε1、啟發(fā)函數(shù)影響因子ε2、信息揮發(fā)因子ε3、迭代次數(shù)t、最大迭代次數(shù)tmax等參數(shù)進(jìn)行設(shè)置,令初始時刻m只螞蟻所在節(jié)點x=o。

(2)根據(jù)生成的施工區(qū)域和各節(jié)點的坐標(biāo)，計算各節(jié)點間的行駛方向和風(fēng)浪的相對大小和方向。

(3)根據(jù)氣象水文條件和船舶本身狀況計算船舶在各路段的自由行駛速度Vf，并據(jù)此求出對應(yīng)的風(fēng)浪影響系數(shù)。

(4)螞蟻根據(jù)與x相連的可行解節(jié)點在T時刻的船舶流量和風(fēng)浪影響系數(shù)等計算從節(jié)點x到可行節(jié)點各自的行程時間，然后按照式(28)依概率選擇下一節(jié)點。

(5)令T=T+Txr(r為選擇的下一節(jié)點，Txr為路段xr行程時間)、x=r，返回(4)繼續(xù)選擇下一節(jié)點直到到達(dá)目標(biāo)點。

(6)螞蟻k(k=1,2,…,m)到達(dá)目標(biāo)點則按照實際用時記錄此路徑，依次完成所有m只螞蟻循環(huán)，記錄此次迭代中的最優(yōu)路徑及用時，并按照式(29)對信息素進(jìn)行更新。

(7)令t=t+1,x=o,T=T0若t≤tmax轉(zhuǎn)向(4),進(jìn)行下一次迭代，否則轉(zhuǎn)向(8);

(8)輸出最優(yōu)解。

5 算例

在此給出如圖1所示的從煙臺到天津的簡易模擬航線路網(wǎng)(箭頭表示只能單向通行)，陰影區(qū)域為沉點存在的區(qū)域，用MATLAB編程求解。利用隨機算法在航線網(wǎng)絡(luò)中隨機生成3個長度為6～8 km的互不影響矩形施工區(qū)域，施工區(qū)域存在的時間為[T0+Xk,T0+Yk](0≤Xk≤Yk≤24 h)。在各施工區(qū)域內(nèi)隨機生成一個施工船沉點，若沉點在航線路段當(dāng)中，則該路段行駛時間為無窮大，以15 min為一個時段。

圖1 模擬航線網(wǎng)Fig.1 Analog ship route network

設(shè)置起點的標(biāo)號為0，終點的標(biāo)號為1 000，1號航線的各點標(biāo)號為101、102、103、104、105，同理為其他各航線進(jìn)行標(biāo)號。圖1中在1、3、5航線內(nèi)生成的施工區(qū)航線網(wǎng)絡(luò)及相應(yīng)標(biāo)號分別如圖2(a)～圖2(c)所示。

圖2 海上施工區(qū)域航線網(wǎng)絡(luò)Fig.2 Route network in offshore construction area

現(xiàn)選擇某型號散貨船進(jìn)行計算，該船舶船長180 m，船寬32.2 m，載重量為48 103 t，滿載排水量54 385 t，吃水11.6 m，靜水營運航速為6.5 m/s，真風(fēng)方向為西偏北30°，波長取0.8倍船長。

在用蟻群算法求解時，取m=30,ε1=1,ε2=1,ε3=0.1,tmax=200, 求得時變網(wǎng)絡(luò)中各代最短用時和平均用時隨迭代次數(shù)的變化曲線如圖3所示。

圖3 時變網(wǎng)絡(luò)各代最短用時和平均用時對比Fig.3 Comparison of shortest and average time of each generation in time-varying networks

最終通過MTALAB編程求得動態(tài)網(wǎng)絡(luò)的最佳路徑為0→501→502→503→71→81→91→92→93→94→95→504→505→1 000，用時86 937 s。在相同情況下，求得1、2、3、4航線的最短用時分別為90 958、91 704、91 818、92 504 s。最佳路徑比4航線路徑少用時5 567 s,減少用時6.4%。

6 結(jié)論

通過建立船舶的航線誘導(dǎo)模型利用MATLAB軟件進(jìn)行算例仿真驗證，說明了本文提出的船舶航線誘導(dǎo)方法能為船舶在復(fù)雜航線網(wǎng)絡(luò)下行駛提供航行路線決策幫助，可以得出以下結(jié)論。

(1)減少總用時時間。通過算例將求解的路線與其他路線的最優(yōu)結(jié)果進(jìn)行對比，結(jié)果表明在此次算例中最佳航線比其他航線最優(yōu)結(jié)果最高可減少用時6.4%，證明了船舶動態(tài)航線誘導(dǎo)的更優(yōu)性和模型的可行性。

(2)提高模型適用性。充分考慮了海上船舶航行的特點，在建模時考慮施工水域的影響、海風(fēng)和波浪的動態(tài)變化影響以及沉點的影響，對路段行駛時間求解時對BPR模型進(jìn)行改進(jìn)并結(jié)合排隊論的理論，使得模型比以往更加適用于船舶。

當(dāng)前模型與實際情況還有些區(qū)別，缺乏對突發(fā)情況的考慮，風(fēng)浪對船舶的影響的計算也還不夠精確，這是今后改進(jìn)的方向。