張 赫, 閆建鑫, 邢江豪

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

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

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

1 海況對(duì)船舶航行的影響

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

1.1 基本阻力

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

(1)

CR=Cf+ΔCf+Cr

(2)

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

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

(3)

(4)

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

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

1.2 空氣阻力

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

(5)

(6)

(7)

(8)

lgY=αw+βwlgX

(9)

式中：ρa(bǔ)為空氣密度(1.226 kg/m3)；Cw為艏艉方向上風(fēng)力系數(shù)；Va為相對(duì)風(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)行計(jì)算[6-8]：

(10)

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

(11)

1.3 洶濤阻力

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

(12)

(13)

1.4 船舶自由行駛速度

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

(14)

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

2 路段行程時(shí)間計(jì)算

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

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

(15)

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

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

(16)

綜上，ij路段總行程時(shí)間tij為

tij=tx+td

(17)

3 模型構(gòu)建

3.1 影響因素分析

3.1.1 施工水域影響

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

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

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

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

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

3.1.4 沉點(diǎn)的影響

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

3.2 模型建立

在對(duì)船舶航線誘導(dǎo)時(shí)，將航線網(wǎng)絡(luò)抽象成圖論網(wǎng)絡(luò)G=(P,E,T)，其中P是節(jié)點(diǎn)集合表示航道節(jié)點(diǎn)；E是弧的集合表示路段；T={tij(Ti)|(i,j)∈E}是各弧的權(quán)值，表示在此路段的行程時(shí)間tij(Ti)為船舶在Ti時(shí)刻從節(jié)點(diǎn)i航行到節(jié)點(diǎn)j所需花費(fèi)的時(shí)間。設(shè)施工區(qū)k(k=1,2,…，l)的節(jié)點(diǎn)集合為Pk存在的時(shí)間為[Tks,Tke]，從起點(diǎn)到達(dá)施工區(qū)k的時(shí)間為tk，非施工區(qū)節(jié)點(diǎn)結(jié)合為P0。則建立船舶動(dòng)態(tài)時(shí)變網(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在時(shí)刻被路徑使用則等于1，否則等于0；式(20)表示初始點(diǎn)只能是出發(fā)點(diǎn)不能是到達(dá)點(diǎn)，終點(diǎn)只能是到達(dá)點(diǎn)不能是出發(fā)點(diǎn)，中間節(jié)點(diǎn)的到達(dá)和出發(fā)次數(shù)相等；式(21)表示所有路段最多只能被選擇一次；式(22)為施工區(qū)內(nèi)沉點(diǎn)對(duì)行程時(shí)間的影響；式(23)、式(24)表示船舶行駛到施工區(qū)域位置時(shí)若在施工區(qū)域存在的時(shí)間范圍內(nèi)則將區(qū)域內(nèi)節(jié)點(diǎn)加入到航線網(wǎng)絡(luò)之中，否則保留原節(jié)點(diǎn)；式(25)表示船舶各路段自由行駛速度的計(jì)算；式(26)為風(fēng)速與位置的對(duì)應(yīng)關(guān)系式；式(27)為波浪振幅與位置的對(duì)應(yīng)關(guān)系式；行程時(shí)間計(jì)算過(guò)程中的其他參數(shù)的計(jì)算公式約束第1節(jié)與第2節(jié)已經(jīng)說(shuō)明，此處不再列出。

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

4 蟻群算法求解

(28)

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

(29)

(30)

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

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

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

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

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

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

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

(6)螞蟻k(k=1,2,…,m)到達(dá)目標(biāo)點(diǎn)則按照實(shí)際用時(shí)記錄此路徑，依次完成所有m只螞蟻循環(huán)，記錄此次迭代中的最優(yōu)路徑及用時(shí)，并按照式(29)對(duì)信息素進(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所示的從煙臺(tái)到天津的簡(jiǎn)易模擬航線路網(wǎng)(箭頭表示只能單向通行)，陰影區(qū)域?yàn)槌咙c(diǎn)存在的區(qū)域，用MATLAB編程求解。利用隨機(jī)算法在航線網(wǎng)絡(luò)中隨機(jī)生成3個(gè)長(zhǎng)度為6～8 km的互不影響矩形施工區(qū)域，施工區(qū)域存在的時(shí)間為[T0+Xk,T0+Yk](0≤Xk≤Yk≤24 h)。在各施工區(qū)域內(nèi)隨機(jī)生成一個(gè)施工船沉點(diǎn)，若沉點(diǎn)在航線路段當(dāng)中，則該路段行駛時(shí)間為無(wú)窮大，以15 min為一個(gè)時(shí)段。

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

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

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

現(xiàn)選擇某型號(hào)散貨船進(jìn)行計(jì)算，該船舶船長(zhǎng)180 m，船寬32.2 m，載重量為48 103 t，滿載排水量54 385 t，吃水11.6 m，靜水營(yíng)運(yùn)航速為6.5 m/s，真風(fēng)方向?yàn)槲髌?0°，波長(zhǎng)取0.8倍船長(zhǎng)。

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

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

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

6 結(jié)論

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

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

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

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