張昊，陶寧蓉*,2，楊男

(1.上海海洋大學，工程學院，上海 201306；2.上海海洋可再生能源工程技術研究中心，上海 201306)

0 引言

隨著經濟全球化的快速發(fā)展，各國經濟越來越依賴于海上運輸，然而海上事故頻發(fā)，海上運輸安全問題一直備受社會關注[1]。其中，最受關注的當屬海上溢油事故。根據《1990 年國際油污防備、反應和合作公約》，海上溢油事故是指原油及其相關煉制品等油品進入海洋或河流的事故[2]。海上溢油事故作為一種突發(fā)事件，不僅會造成巨大經濟損失，而且溢出的油品會嚴重污染海洋環(huán)境、破壞海洋生態(tài)。據統(tǒng)計，1973—2018 年，我國沿海港口共發(fā)生3336 起船舶溢油事故，平均每年發(fā)生76 起。因此，及時合理地調度溢油事故應急物資至關重要。

目前，國內外針對海上溢油事故應急物資調度的研究較少。張莉等[3]針對事故點的需求緊迫度，研究如何優(yōu)先配送緊迫度高的事故點以降低總成本。郝國柱等[4]考慮海上溢油事故應急物資需求量和調度時間不確定，引入三角模糊函數，建立了總成本最小和調度時間最少的雙目標優(yōu)化模型。Zhang等[5]針對海上溢油事故多配送中心協(xié)同調度進行研究，提出一種基于粒子群優(yōu)化的啟發(fā)式方法。李松等[6]考慮各應急物資集散點和船舶并行作業(yè)的特點，將物資需求量轉換為當量體積，以船舶航行時間最小為目標，建立了大型海上溢油事故應急物資聯動調度模型。張聆曄等[7]統(tǒng)籌海上應急物資調度與陸上補給應急物資調度，構建兩階段應急物資動態(tài)優(yōu)化調度模型，并提出混合啟發(fā)式算法進行求解。汪強等[8]結合海上溢油事故特點，構建了應急救援中心供應充足條件下運輸物資延誤總時間最小的數學模型，并用遺傳算法進行求解。

溢油事故發(fā)生時，受到氣象條件和海況的影響，油品進入海洋后會漂移和擴散，進而導致救援需求點變動。應急救援如果不考慮需求點的漂移擴散特性，會造成救援船舶偏離需求點實際位置而影響救援方案實施。劉曉佳等[9]針對溢油的擴散性和應急物資需求的多樣性引入三角模糊數確定應急物資需求區(qū)間，以船舶運輸總時間最短為目標建立調度模型，并運用于遺傳算法對其進行求解。李晶等[10]分析了海上油膜漂移特性和海上溢油點對救援物資需求量不確定的特點，建立救援成本最小、時間最短的雙目標規(guī)劃模型并求解。然而，以上研究都沒有考慮油膜漂移擴散對應急調度路徑規(guī)劃的影響。王軍等[11]考慮單一遇險目標受海上風浪影響產生漂移，建立了水陸兩階段應急物資優(yōu)化調度模型，但是該研究在海上只有一個需求點，而海上溢油事故發(fā)生后的求援需求點通常是多個。Zhang等[5,12]考慮了多個需求點的救援，結合需求點的漂移特性提出了兩階段優(yōu)化模型，分別采用了改進粒子群算法和基于改進粒子群算法的混合啟發(fā)式算法進行求解獲得Pareto解集，但是該研究在評估環(huán)境污染時沒有涉及需求點漂移導致的油膜面積動態(tài)變化。

綜上所述，已有研究大多沒有在應急救援調度中考慮需求點受海上風浪影響而產生位置漂移的時變特性。本文在綜合考慮海上溢油事故需求點漂移特性以及油膜擴散造成的環(huán)境損失后，建立最小化調度運輸成本和環(huán)境污染損失的數學模型，并在傳統(tǒng)遺傳算法的變異操作中加入模擬退火的思想，提出改進的遺傳模擬退火混合優(yōu)化算法進行求解。

1 問題描述

當海上發(fā)生溢油事故時，相關部門首先迅速定位事故發(fā)生點，分析氣象條件、海況、溢油類型、溢油規(guī)模等信息，計算溢油擴散數據和應急物資需求量；接著協(xié)調各應急救援中心做出調度安排；最后實施應急措施，對泄漏油品等污染物進行圍控、回收、處理。如圖1所示，本文研究的問題為：在海上溢油事故發(fā)生初期，如何調度多個應急救援中心的物資，對污染區(qū)域內多個動態(tài)需求點實施求援，從而將事故控制在較小范圍，實現救援運輸成本和環(huán)境損失的最小化。基本假設如下：

圖1 海上溢油事故應急物資調度Fig.1 Emergency material scheduling for offshore oil spill accident

(1)溢油地位置、物資需求量和溢油地氣象條件等基本信息可通過衛(wèi)星技術計算獲取。

(2)考慮油膜的漂移和擴散特性，事故存在多個求援需求點，且隨著油膜的漂移擴散，需求點的位置是動態(tài)變化的。

(3)救援船舶為負責對溢油進行圍控、導流的溢油工作船舶。

(4)當應急物資運送到需求點后救援工作即刻啟動。

(5)海上沒有設置中間節(jié)點進行船舶之間的物資交付。

(6)船舶在海上航行速度保持一致。

模型中符號說明如下。

(1)集合

V——海上溢油事故應急物資調度節(jié)點集合，V=A∪B，其中，A表示需求點集合，B表示應急救援中心集合；

K——應急救援船舶集合。

(2)參數

U——溢油總量；

T0——調度開始時間；

Qk——船舶k的最大運載能力；

qi——需求點i的需求量；

vk——船舶k的航行速度；

vi——需求點i的移動速度；

c1——單位距離運輸成本；

c2——船舶調用成本；

ρw——海水密度；

ρ0——溢油密度；

β——風速與水平坐標的傾角。

(3)中間變量

S(T)——T時刻的溢油油膜面積；

θij——船舶從點i向j航行時，點i到點j之間的距離連線與j點漂移路徑的夾角；

(4)決策變量

2 油膜漂移擴散分析

2.1 溢油油膜的擴散

關于海上溢油擴散行為的研究，目前最具代表性的是Fay 提出的三階段擴展模型[13]?；谠撃Ｐ?，Leendertse[14]考慮風作為影響因素，通過溢油總量U、溢油事故等待救援時間t、海水密度ρw和溢油密度ρ0計算需求點油膜擴散速度，即

式中：Δρ=ρw-ρ0。

2.2 溢油油膜的漂移

依據溢油類型、溢油發(fā)生時間和地點、溢油量和持續(xù)時間等信息，計算油膜漂移的具體數值。本文基于Li等[15]的研究，認為油由大量質量相等的粒子構成，并受風、波浪和湍流擴散影響，則通過需求點i處的風速和洋流速度，計算需求點i的漂移速度為

2.3 需求點時變位置分析

研究海上溢油事故救援需求點位置的時變特性，使應急調度方案與需求點實際位置匹配，是提高求援響應效率的前提與保障。圖2 為需求點位置漂移示意圖。

圖2 需求點位置漂移示意圖Fig.2 Schematic diagram of demand point location drift

根據油膜的漂移速度和方向、船舶行駛速度、應急救援中心/當前船舶所在的位置等信息可以計算船舶到達下一個目標點的時間以及該點的實際位置。如圖2 所示，當救援船舶k從出發(fā)點i向目標點j行駛時，出發(fā)時刻兩點間的距離為

三角形三邊關系為

夾角θij的計算公式為

則已知需求點漂移速度和方向、船舶速度，可推出船舶所需航行時間為

2.4 環(huán)境損失評估

大多數針對石油等有害物質造成的環(huán)境損失評估研究需要大量的數據和很長時間。為降低評估成本，本文采用弗羅里達公式[9]，該公式被證明是一種用于溢油生態(tài)損害評估的快速、低成本方法。環(huán)境污染損失E為

式中：E(T)——T時刻的環(huán)境污染損失；

R——基準利率；

店主懶洋洋按了播放，迪斯科和說唱又想起。警察一直耐心得聽到了“送到派出所”，琢磨了半天，轉身對左小龍道：“沒問題啊，沒反黨啊。”

L——地理位置影響因素(沿海為8，近岸為5，離岸為1)；

M——特殊管理區(qū)域(是為2，否為1)；

A*——受影響動物棲息地的附加費用(按照生物類型，每平方米1～50元)；

PC——污染物特質(重油為8，中重油為4，輕油為1)。

由于本文研究關注油膜的漂移擴散特性，即油膜面積S隨時間變化。當攔油類物資到達需求點實施救援后，油膜漂移擴散可得到有效控制，因此本文通過救援需求點的位置信息來計算某時刻t下的油膜面積S。如圖3 所示，假設在時刻t跟蹤到事故區(qū)域內的若干救援點，則油膜面積S的具體計算方法為：(a)獲取時刻t的所有需求點位置信息；(b)根據離散點最小邊界查找法確定邊界點坐標；(c)用向量法計算油膜污染區(qū)域面積S。

圖3 最小邊界查找示意圖Fig.3 Schematic diagram of minimum boundary search

式中：a,b,…,f——確認后的邊界點坐標信息。

2.5 數學模型

2.5.1 目標函數

基于以上分析，本文模型的優(yōu)化目標為最小化應急救援運輸成本和環(huán)境損失，即

式中：a1、a2——權重系數，a1+a2=1；

Z1——應急救援運輸成本，如式(11)所示，該部分包括派出救援的船舶固定成本和油耗成本；

Z2——環(huán)境污染損失，如式(12)所示，取救援實施前后的差值；

c1——單位距離運輸成本；

c2——船舶派遣成本；

E(T0)——救援前的環(huán)境污染值；

E——所有需求點的救援工作到位后的環(huán)境污染值。

2.5.2 約束條件

式(13)表示船舶k的 載重約束。式(14)和式(15)為救援需求點距離約束，式(14)表示船舶k在需求點i時，需求點i,j之間的距離；式(15)表示船舶k由i到j的實際航行距離。式(16)和式(17)為時間約束，式(16)表示若船舶k不經過i點，則=0；式(17)表示船舶k由i到j的運輸時間為船舶k到達j點的時刻減到達i點的時刻。式(18)和式(19)是對求援船舶出發(fā)點的約束，式(18)表示每艘救援船的出發(fā)點必須是應急救援中心；式(19)表示每艘救援船僅能從多個應急救援中心的其中一個出發(fā)。式(20)約束救援船不可在應急救援中心之間往返。式(21)表示救援船在各救援需求點的流量守恒。式(22)表示每個救援需求點只需一次救援。

3 算法設計

本文模型屬于多中心多需求點且需求點漂移的路徑規(guī)劃問題，而路徑規(guī)劃問題已被證明是NP-hard問題，在實際應用中，這類問題很難在有限時間內通過精確算法求出高效解。遺傳算法是一種強魯棒性的啟發(fā)式隨機搜索算法，具有較好的全局搜索性能，但基本遺傳算法存在過早收斂于局部最優(yōu)解的現象，而模擬退火算法通過以一定概率接受“劣等解”具有較好的局部搜索能力。因此本文采用基于遺傳算法和模擬退火算法的混合算法，在遺傳變異操作中引入模擬退火思想，通過多策略變異方法對個體進行連續(xù)多次變異并以一定概率接受差解，算法流程如圖4所示。

圖4 遺傳模擬退火算法流程圖Fig.4 Flow chart of genetic and simulated annealing algorithm

3.1 編碼與解碼

本文編碼采用整數編碼方式，染色體為J+I-1(需求點數+應急物資中心點數-1)的整數序列，其中，數值1～J分別為需求點編號，而J+1～J+I-1 的數值為應急物資中心的分割碼。以8 個需求點、2 個應急物資中心為例，數值9 為2 個應急物資中心的分割碼，如圖5所示。

圖5 染色體解碼示意圖Fig.5 Schematic diagram of chromosome decoding

解碼時，從染色體的第1 個需求點7 開始累加物資需求量，假設到需求點5累加的需求量超過船舶載重，則在需求點5 前斷開，則7-3-2 由1 號應急物資中心的船舶1 求援，從需求點5 開始重復上述操作。最終解碼的救援方案為1 號應急物資中心船舶1的路徑為0-7-3-2-0；船舶2的路徑為0-5-4-0；2號應急物資中心船舶1的路徑為0-6-1-8-0。

3.2 適應度函數

各個染色體的適應度值由模型的目標函數構建，目標函數值越小，適應度越大。

式中：fi——染色體i的適應度值；

假設當前需求點序列為0-7-3-2-0，如圖6 所示，首先，獲取救援中心0和需求點7的位置以及需求點7 的漂移方向和速度，根據式(7)計算0-7 的運輸時間，記錄船舶1 行駛距離·vk；更新時刻剩余需求點的位置；基于更新后的各點位置以及下一個需求點3的漂移速度和方向，計算下一步運輸時間和距離，如圖6(b)所示，重復以上操作直到船舶返回應急救援中心，計算Z1。在救援結束后獲取所有需求點更新后的位置信息，根據最小邊界查找法獲得油膜面積，計算Z2。

圖6 船舶救援路徑示意圖Fig.6 Schematic diagram of ship rescue route

3.3 選擇操作

選擇操作采用精英保留策略和輪盤賭策略，將種群中的個體按照適應度大小排序，通過精英保留策略保留適應度最大的若干個體。剩下的個體采用輪盤賭策略，即將剩余個體適應度值與當代最優(yōu)適應度值相除得到選擇概率，進行選擇。

將輪盤賭策略選擇出的個體進行交叉、變異操作，將精英保留策略保留的個體直接插入交叉、變異之后的種群中，計算種群數量，如果數量小于種群規(guī)模N，則隨機生成新個體，確保種群規(guī)模的一致性。

3.4 交叉操作

多點交叉可以避免過早收斂，但是過多的交叉點會影響優(yōu)良基因的傳承，本文權衡利弊后采用交叉點為2 的交叉操作。如圖7 所示，交叉后找出單個個體中重復的點，對應父代X1、X2將個體非交叉區(qū)域中重復點替換為因交叉缺失的點，得到子代x1，x2。

圖7 交叉操作示意圖Fig.7 Schematic diagram of cross operation

3.5 模擬退火變異

需求點漂移影響應急救援路徑的選擇，單點變異較難跳出局部最優(yōu)，故將模擬退火思想引入遺傳變異操作，通過交換、逆轉和插入的多策略變異方法對個體進行連續(xù)多次變異，并以Pt=exp的 概 率 接 受“ 劣 等 解”，ΔF=-fi如圖8 所示，交換結構為隨機選擇染色體上的兩個位置，然后將這兩個位置上的元素進行交換；逆轉結構為隨機選擇染色體上的兩個位置，將兩個位置之間的元素進行逆序排列；插入結構為隨機選擇染色體上的兩個位置，將第1個位置后的元素插入到第2個元素的后面。

圖8 染色體交換、逆轉、插入操作示意圖Fig.8 Schematic diagram of chromosome exchange,reversal and insertion operations

4 算例驗證與結果分析

4.1 實例數據

本文基于2011 年蓬萊石油泄露事故的相關數據[7]開展實驗。該事故的應急物資運輸行動發(fā)生在事故發(fā)生6 h后。根據溢油實際信息得到溢油事故發(fā)生位置為離岸非特殊管理區(qū)域，則L=1 、M=1。根據渤海生態(tài)環(huán)境[16]，受影響動物棲息地的附加費用A*=10，溢油類型為中重油PC=4。經過測算，溢油造成的一大片油膜污染區(qū)域有30 個求援需求點需要圍油欄進行圍控。該事故溢油量為250 t，渤海海域水密度ρw為1.025 g·cm-3，渤海海域的風向為南風，海平面10 m的風速為10.6 kn，需求點受到兩種洋流影響，分別是流速1.8 kn、流向東偏南37°漂移的洋流和流速為1.2 kn、流向為北偏西24°漂移的洋流。

各需求點的當前坐標、漂移速度和救援物資需求量等信息如表1所示。為方便計算，本文采用高斯投影法將經緯度坐標轉換成為直角坐標。應急救援中心位置信息如表2 所示。船舶平均航行速度為15 m·s-1，船舶額定裝載量為500 t，單位運輸成本c1為450 元·km-1，固定運營成本c2為10000元·艘-1，權重因子a1和a2均取0.5。

表1 救援物資需求點信息Table 1 Rescue material demand point information

表2 應急救援中心信息Table 2 Information of emergency rescue center

4.2 實例驗證分析

本模型基于MATLAB R2021b 求解，在PC 端運行，操作系統(tǒng)為Window11，運行內存16 G，CPU為Inter(R)Core(TM)i5-10600kf，主頻4.1 GHz。

在算法參數選擇上，本文參照已有研究，采取逐個調整一個參數固定其他參數的方法，最終將參數確定為：種群規(guī)模N=50，最大迭代次數Gmax=400，交叉概率Pc=0.9，變異算子中交換結構概率Ps=0.2、逆轉結構概率Pr=0.5、插入結構概率Pi=0.3，模擬退火初始溫度Tsa=0.025，每次迭代冷卻系數α=0.8，外層循環(huán)最大迭代次數=400，里層循環(huán)最大迭代次數=3。

為驗證算法的有效性，引入傳統(tǒng)遺傳算法、模擬退火算法和參考文獻[5]構建的改進粒子群算法與本文算法進行比較，算法收斂情況如圖9 所示。模擬退火算法收斂緩慢，需要更多次迭代才能得到較好的解，遺傳算法較容易陷入局部最優(yōu)，改進粒子群算法雖然也能夠獲得較好的解，但較易陷入局部最優(yōu)解，本文改進后的遺傳模擬退火混合算法在保證收斂速度的前提下能夠搜索到更優(yōu)解。

圖9 算法收斂性對比Fig.9 Comparison of algorithm convergence

本文算法獲得的最優(yōu)方案如表3所示，即同時派出6 艘救援船，其中大連基地3 艘、煙臺2 艘、威海1 艘，具體救援路線如圖10 所示，環(huán)境污染損失10199241.84 元，救援物資運輸費用625251.89 元，整個溢油事故控制在約5384.68 km2的油膜污染面積內。此外，基于本文模型的求解實現了對救援點漂移位置的追蹤，圖11 為實施救援前后需求點位置差異。

圖10 最優(yōu)救援方案路線圖Fig.10 Roadmap of optimal rescue plan

圖11 救援前后需求點位置示意圖Fig.11 Schematic diagram of demand points before and after rescue

表3 應急救援配送方案Table 3 Emergency rescue distribution scheme

在此基礎上，本文分析了制定應急救援方案時考慮需求點漂移的必要性。如表4所示，不考慮漂移得出的救援方案中所定位的需求點位置存在偏差，因此在實施時增加了額外的行駛距離，同時由于救援工作延誤，導致油膜污染面積增加。與不考慮需求點漂移下獲得的救援方案相比，本文方法在航行距離上減少了9.11%、環(huán)境污染降低了41.17%。

表4 考慮/不考慮需求點漂移因素的方案比較Table 4 Scheme comparison with/without demand point drift

4.3 數值實驗分析

設計6組不同規(guī)模的算例，分別采取本文所提的混合算法和遺傳算法、模擬退火算法進行求解。算例設計主要考慮以下因素：地理數據、船舶服務的需求點數量、需求點漂移擴散速度。地理數據依據Solomon實例數據集種R101、C101的數據，需求點在問題集R101中隨機生成，而在問題集C101中采用聚類的方式產生。由于基準算例中沒有漂移速度，本文為每個需求點隨機生成-2＜vj＜2 的速度。對相同問題集分別取前25、50、100 個需求點構成小規(guī)模、中等規(guī)模和大規(guī)模的數據進行實驗。

每次實驗運行5次記錄最大值、最小值和平均值，并通過計算每組算例中4種算法的目標值偏差來評估算法全局搜索能力和算法穩(wěn)定性。如表5所示，本文算法與改進粒子群算法(IPSO)、模擬退火算法(SA)、遺傳算法(GA)求解結果的平均偏差分別為1.63%、1.76%和20.39%。在小規(guī)模問題中，本文算法、IPSO 和SA 都可獲得最優(yōu)解，隨著規(guī)模增大，本文算法的優(yōu)越性凸顯，在R101_100 算例中，IPSO、SA和GA的最優(yōu)解與本文算法相比，分別增加了0.55%、1.41%和35.19%。算法穩(wěn)定性上，在算例C101-25～R101-50 中，本文算法穩(wěn)定性突出；在算例C101-100、R101-100中，本文算法與IPSO穩(wěn)定性接近，優(yōu)于SA和GA算法。綜合比較，本文算法在保持穩(wěn)定性的基礎上在不同類型、不同規(guī)模的算例中具有更強的尋優(yōu)能力。

表5 算法求解結果對比Table 5 Comparison of algorithm solution results

4.4 靈敏度分析

本文進一步分析了船舶容量、航行速度和需求點漂移速度等因素的影響。針對每個因素，分別在原數據基礎上取-40%、-20%、20%、40%進行實驗，實驗結果如圖12～圖14所示。

由圖12 可知：船舶容量增大可以有效降低運輸成本，船舶容量減小時，船舶數量增加，頻繁從應急救援中心派遣船舶直接導致船舶總運輸成本增加；只要船舶數量充足，其容量的變化對油膜污染損失的影響并不明顯。由圖13 可知：當船舶航行速度加快時，船舶更快到達需求點，可以更有效地控制油膜污染面積，防止其擴散；由于本文算例基于渤海海域，溢油油膜實際是朝向海岸線擴散，因此出現船舶航行速度越快，救援越及時，反而需求點處在與陸地應急救援中心較遠的位置，進而運輸成本上升，所以單純增加船舶航行速度對運輸成本的影響不大。由圖14可知：當需求點漂移速度增大時，油膜擴散加快，油膜污染面積明顯增大；而漂移速度增大時，由于受到海上風浪的影響加劇，需求點之間的距離增加，進而增加了船舶總航行距離。

圖12 船舶容量靈敏度分析Fig.12 Sensitivity analysis of ship capacity

圖13 船舶航行速度靈敏度分析Fig.13 Sensitivity analysis of ship navigation speed

圖14 需求點移動速度靈敏度分析Fig.14 Sensitivity analysis of moving speed of demand point

5 總結

(1)實例結果驗證了海上溢油應急救援調度問題中考慮需求點漂移的重要性。與不考慮需求點漂移獲得的救援方案相比，本文方法在總航行距離上減少了9.11%，環(huán)境污染降低了41.17%。本文方法能夠獲得較優(yōu)的調度方案，實現更準確高效的應急救援，使險情得到及時控制，減少環(huán)境污染損失。

(2)算例分析結果表明：本文算法與IPSO、SA、GA算法求解結果的平均偏差分別為1.63%、1.76%和20.39%。在小規(guī)模問題中，本文算法、IPSO 和SA都可獲得最優(yōu)解，隨著規(guī)模增大，本文算法的優(yōu)越性凸顯，在R101_100 算例中，IPSO、SA 和GA 的最優(yōu)解與本文算法相比，分別增加了0.55%、1.41%和35.19%。在算法穩(wěn)定性上，在算例C101-25～R101-50中本文算法穩(wěn)定性突出，在算例C101-100、R101-100 中本文算法與IPSO 穩(wěn)定性接近，優(yōu)于SA 和GA 算法。綜合比較，本文算法在保持穩(wěn)定性的基礎上在不同類型、不同規(guī)模的算例中具有更強的尋優(yōu)能力。

(3)船舶容量對運輸成本的影響更靈敏，而對油膜污染面積影響不明顯。船舶航行速度對油膜污染面積影響更靈敏，而對運輸成本的影響較小。因此，提升船舶容量可以明顯降低應急救援運輸成本，提升船舶航行速度可以更好地降低環(huán)境污染。