張 赫, 閆建鑫, 邢江豪

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

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

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

1 海況對船舶航行的影響

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

1.1 基本阻力

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

(1)

CR=Cf+ΔCf+Cr

(2)

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

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

(3)

(4)

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

表1 粗糙度補貼系數Table 1 Roughness subsidy coefficients

1.2 空氣阻力

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

(5)

(6)

(7)

(8)

lgY=αw+βwlgX

(9)

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

表2 不同船型船舶的受風面積回歸系數Table 2 Regression coefficient of windage area of different types of ships

風力角α和風力系數Ca可根據巖井聰公式[式(10)和式(11)]進行計算[6-8]：

(10)

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

(11)

1.3 洶濤阻力

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

(12)

(13)

1.4 船舶自由行駛速度

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

(14)

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

2 路段行程時間計算

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

各路段行程時間分為路段行駛時間和航線網絡交叉點延誤時間兩個部分。在沒有航線交叉的路段內，路段行程時間就等于路段行駛時間。船舶在路段的行駛時間與路段船舶流量密切相關，根據風浪對船舶航行的影響對美國聯邦公路局函數即BPR函數進行改進，建立路段行駛時間的表達式：

(15)

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

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

(16)

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

tij=tx+td

(17)

3 模型構建

3.1 影響因素分析

3.1.1 施工水域影響

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

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

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

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

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

3.1.4 沉點的影響

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

3.2 模型建立

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

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

4 蟻群算法求解

(28)

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

(29)

(30)

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

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

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

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

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

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

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

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

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

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

5 算例

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

圖1 模擬航線網Fig.1 Analog ship route network

設置起點的標號為0，終點的標號為1 000，1號航線的各點標號為101、102、103、104、105，同理為其他各航線進行標號。圖1中在1、3、5航線內生成的施工區(qū)航線網絡及相應標號分別如圖2(a)～圖2(c)所示。

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

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

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

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

最終通過MTALAB編程求得動態(tài)網絡的最佳路徑為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 結論

通過建立船舶的航線誘導模型利用MATLAB軟件進行算例仿真驗證，說明了本文提出的船舶航線誘導方法能為船舶在復雜航線網絡下行駛提供航行路線決策幫助，可以得出以下結論。

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

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

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