基于海洋溢油溯源數值模擬方法的長島保護區污染風險分析

陳勤思，胡松，魏永亮

（1.上海海洋大學 海洋生態系統與環境實驗室，上海 201306；2.上海海洋大學 海洋科學學院，上海 201306）

渤海是我國重要的石油產區，也是重要的海上運輸通道。隨著海洋油氣業和海上運輸業的發展，渤海逐漸成為溢油污染事故的高發區。渤海海域分布著眾多海上石油平臺，石油開采過程中的泄漏和船舶事故引起的溢油污染事件時有發生，給區域海洋生態環境造成損害。柯麗娜等總結了渤海海域中來自船舶和石油平臺的溢油風險源分布情況（圖1）[1]。此外，一些可能來源于石油管道滲油、沉船和非法傾倒等多種原因引發的溢油，常常因其來源不明而形成無主漂油，這類溢油的污染危害同樣不可小覷，需要對此進行風險分析從而及時采取防污染的措施。

圖1 渤海海域溢油風險源分布

山東長島國家級自然保護區（以下簡稱長島保護區）是位于渤海海峽南部的國家級自然保護區，具有重要的生態和旅游價值，是鳥類和太平洋斑海豹生存和繁殖的重要區域。但是，保護區海域時常受到來源不明的溢油污染物的影響[2-5]，不僅給當地的漁業、旅游業等產業造成巨大損失，而且危及鳥類和斑海豹等保護物種的生存。溢油可能會隨著潮流、風等到達保護區海域，也可能是海域內的船舶泄漏造成，因此很難從溢油來源進行污染預測，給風險分析帶來挑戰。

以往研究中，海洋溢油溯源數值模擬方法常常被用于解決一類污染受體已知、溢油源未知（或過多）情景的災害預防問題，溯源方法在這類情景中所需要花費的計算時間少[6]，并且可以提供更為全面的潛在溢油泄漏情景和污染物擴散軌跡信息[7-9]。不同于一些隨機和典型情景的溢油事故模擬方法[10-11]，溯源方法從污染受體位置出發開始，采用溯源軌跡對溢油的可能來源進行追蹤，判斷可能對污染受體造成危害的溢油來源區域和概率，從而對高概率區域和低概率區域執行相應等級的監測和預防措施[12]。Ciappa 等將該方法應用在地中海的海洋保護區溢油風險分析中，評估了該區域來自船舶的潛在污染風險和溢油監測的重點海域[8]。高松等基于該方法對秦皇島海域無主漂油的可能來源區域進行了估計，為溢油來源的排查和應急處置提供了幫助[3,13]。位于渤海的長島保護區海域位置特殊，受到周期性變化的強潮流和季風的共同影響，在進行溢油污染風險分析時應當盡可能全面地考慮這些狀況。

本文采用海洋溢油溯源數值模擬方法對長島保護區海域進行溢油災害風險評估和分析。研究使用FVCOM模擬的潮流數據和ECMWF提供的ERA5海面10 m 風數據，在溢油模型GNOME （General NOAA Oil Modeling Environment）中進行溯源模擬。首先，通過一次真實溢油事故的后報驗證了溯源模擬方法的有效性，然后將這套資料和方法應用于渤海長島保護區溢油災害風險的分析，評估了月平均統計條件下保護區溢油污染的來源分布和危害風險，探討了該方法的適用性，以期為我國海洋溢油風險評估相關的研究工作提供參考。

1 研究區域和數據

1.1 研究區域

渤海海峽是我國渤海的唯一出口，地理位置特殊。長島保護區海域位于渤海海峽南側，自北向南與老鐵山水道、北砣磯水道和長山水道等主要航道相鄰。長島保護區海域受強潮流影響，海流多為東-西向往復流。區域盛行風向隨季節轉換而有明顯變化，全年中春、秋、冬季的季風強度較大且持續時間較長，夏季大風天氣較少。

在長島保護區海域的三處位置周圍選取了66個位置并設置為污染受體地點，即圖2（b）中A、B和C 三處位置所示的受體點。對于A、B 和C 處的每一部分，相鄰受體點之間的距離是相等的。研究不設置潛在溢油事故的發生地點，油粒子從受體點位置釋放，隨時間逆向回溯其潛在的來源地點和粒子溯源軌跡。

圖2 研究區域

1.2 數據

海面10 m 高的風場數據來自ECMWF 提供的第五代大氣再分析產品ERA5再分析數據（https://cds.climate.copernicus.eu/#!/home）。該產品同化了衛星數據、現場觀測數據等多種數據，具有較高的精度。本文采用2000－2019年20年間的逐時風場資料，空間分辨率為1/4°×1/4°。在風險分析的研究部分，為體現保護區海域風場季節性差異和多年平均狀態，本文對數據進行了月平均的預處理。

蓬萊19-3 溢油事故的衛星遙感圖像數據來自歐空局ENVISAT 衛星ASAR（Advanced Synthetic Aperture Radar）圖像，并采用TCNNA（Texture-Classifying Neural Network Algorithm）算法[14]提取了圖像中油膜的范圍信息，與前人的研究結果較為一致[15-16]。為了驗證模型的可靠性，本文采用同一次溢油不同時刻的兩張衛星圖像，以較晚時刻的圖像中獲得的溢油時間和位置信息作為溯源模擬的初始條件，以較早時刻的圖像中獲得的溢油時間和位置信息作為溯源模擬結果的驗證。

2 分析方法和模型設置

2.1 分析方法

受體模式（receptor mode）由Galt 等[17]提出并較早應用于溢油溯源模擬中。受體模式基于拉格朗日觀點，將溢油視為以油粒子的形式運動，在溯源模擬中會根據平流方程，逆著洋流和海面風作用方向，在時間上逆向反推油粒子的來源位置和溯源軌跡。在受體模式中，溢油受海流影響的輸運部分被簡化為位移，即油粒子代表的浮油隨著水體以拉格朗日方式運動：

式中：x代表油粒的位置；為表層海流速度。根據Galt 等的描述，風的效應和溢油粒子的擴散通過額外的方法模擬，溢油輸運的最終結果是矢量疊加的效果。Yu等[18]將這個過程表述為：

式（2）中：p(t)為監測到的溢油粒子位置；p(t-1)表示僅受海流和風影響下的前一時刻的位置；和uc(t-1)表示每個時間步長中由風和海流引起的速度。式（3）中：R(t)表示隨機行走產生的距離半徑；uR()t-1為隨機行走的速度。式（4）中：St-1則表示由隨機行走造成的溢油粒子的可能影響范圍。通過式（2）和（4）可以得到溢油來源的高概率區域。

由于油粒子的隨機行走過程在溯源過程中是不可逆的，因此應根據實際情況判斷回溯模擬的參數以滿足對當地環境特征的適用性。關于溯源模擬中的隨機行走帶來的問題可以參考以往研究[6,19]。

本文根據溢油溯源軌跡分析了可能對長島保護區造成污染危害的來源區域，以及來源區域的溢油到達保護區的概率分布。首先，在長島保護區海域的66 個受體點（圖2（b））釋放了用于溯源的油粒子，采用月平均的風數據和隨時間變化的潮流數據模擬得到研究區域的溢油溯源軌跡。然后，從這些軌跡中分析可能對保護區造成污染影響的溢油來源位置，將這些信息按月合并在一起，最終得到了研究區域整體的概率圖。因此，概率圖展示了保護區周圍海域某一位置的溢油到達受體點的概率，其由穿過統計單元的油粒子軌跡數和由受體點釋放的油粒子間的比例組成[7]。本文標記了概率圖中置信度為95%的區域。

2.2 動力模型及驗證

通過FVCOM 海洋模型模擬得到了渤海海域的潮流數據。FVCOM 是基于三維海洋原始方程的非結構網格、有限體積的海洋模式[20]，結合了有限元法擬合邊界、局部加密和有限差分法計算高效等優點。本文建立了渤海三維正壓潮汐模型，模型共生成32 237 個三角元與17 173 個節點。為保證長島保護區海域網格質量，對該海域進行了網格加密，使海域網格最高空間分辨率在100 m 以內。模型開邊界由1/30°分辨率的TPXO9主要分潮的潮汐水位驅動[21-22]，運行時間為2011 年4 月1 日至7 月31 日。此外，在風險分析的研究部分采用模型模擬的5 月1 日至7 月1 日的潮流數據，從而體現出研究區域漲落潮、大小潮等潮汐狀態下溢油溯源模擬的結果。

通過事故同期的潮汐資料對模型水動力模擬結果進行驗證，結果如圖3 所示。本文驗證了2011年6月7－12日P1位置潮位資料與2011年5月18日P2位置潮流資料，潮位與潮流的模式計算結果與事故同期實測資料的比對結果良好。

圖3 FVCOM 模式計算網格及驗證

2.3 溢油模型

GNOME 是美國國家海洋和大氣管理局（National Oceanic and Atmospheric Administration，NOAA）開發的溢油建模工具，可對溢油漂移和擴散軌跡進行預測[23]。GNOME 在診斷模式（diagnostic mode）下可輸出“最佳預測軌跡”和“最小遺憾軌跡”。前者默認輸入信息完全正確，代表了對溢油去向的最佳判斷，后者在模型中包含不確定性，顯示可能受到溢油影響的更大區域。本文采用了“最小遺憾軌跡”模擬結果進行分析。

GNOME 自 早 期 版 本OSSM （On-Scene Spill Model）開始采用受體模式，并支持用于溢油溯源的“Run Backwards”模塊。在本文中，用于溯源模擬的油粒子是在圖2（b）所示的受體點位置釋放的。模型運行中設置的時間步長為30 min，每個受體點在一個時間步長釋放3 個油粒子，并在一個月中連續釋放15 日。不考慮油粒子風化作用和歸宿，單個粒子在模擬中設置的有效時長為30日，即模擬時長超過30 日的粒子會被剔除。因此，模擬過程的整體有效模擬周期為45日。

3 結果與討論

3.1 模擬結果驗證

本研究首先對蓬萊19-3 油田溢油事故中溢油真實的漂移運動情況進行了后報，其目的是驗證數據和溢油溯源模擬方法的有效性，并確定符合本文研究區域環境特征的溢油溯源模擬的適用參數。蓬萊19-3 溢油事故最早發生于2011 年6 月4日，持續泄漏時間在2個月以上，約700桶原油和2 500 桶含油泥漿泄漏到海水中。6 月19 日，兩張不同時刻的衛星圖像記錄了事故中溢油在海面的位置和分布情況，圖4（b）和圖4（c）展示了提取的溢油位置和范圍信息。

圖4 蓬萊19-3事故中不同時刻的溢油位置

以較晚時刻的圖像中獲得的溢油時間和位置信息作為模擬的初始條件（圖5（a）），溯源模擬的開始時刻為2011 年6 月19 日9 時0 分。以較早時刻的圖像中獲得的溢油時間和位置信息作為溯源模擬結果的驗證，模擬結束的時刻為2011 年6 月19 日2 時21 分。圖5（b）展示了本文采用的溢油溯源模擬參數方案在事故后報中得到的模擬結果，且與真實位置接近。結果表明本文的資料和方法可以在短期溯源模擬中得到較好的結果。

圖5 溢油事故溢油溯源模擬結果

由于污染風險分析采用溯源方法評價的研究尚不多見，對于溯源模擬結果的精度評價案例相對較少。因此，本文在表1 中給出了通過此次模擬所優化后的相關參數取值。在GNOME 中，各種輸入項的不確定性部分可以設置為一個區間內的隨機值。比如，風漂移系數默認在1%～4%中取值，并且該值會每15 min重置一次。

表1 溯源模擬中相關參數取值

3.2 溢油概率分布模擬

GNOME 以月平均的風數據和隨時間變化的潮流數據作為驅動，在圖2（b）所示的受體點位置釋放了用于溯源的油粒子，每個受體點每日的釋放頻率為144 個，一個月中連續釋放時長為15 日，每個粒子的有效模擬時長為30 日。然后，從66×144×45 條溯源軌跡中得到可能對保護區造成污染影響的溢油來源位置，這些信息按月合并在一起，最終得到了研究區域整體的概率圖。

溯源軌跡的概率圖顯示了保護區周圍海域某一位置的溢油到達保護區造成污染的概率。圖6不同顏色的區域表示了對應位置溢油的污染概率，黑色實線標記了圖中置信度為95%的區域。結果表明，①受強潮流作用影響，區域內東-西走向的往復流運動使東岸和西岸的區域周圍常年表現為高概率區域，也使保護區內島嶼的東西兩岸更容易受到溢油污染物的影響。②受季節性變化的典型風場影響，高概率區域的主要分布態勢呈現出較為明顯的季節性差異。高概率區域的分布形態在夏季向東南方向延伸，對應季節中盛行東南風，在冬季向西北方向延伸，對應季節中盛行西北風。春夏季節，保護區內島嶼更容易受到東部和南部海域溢油污染的影響，在對應季節中溢油可能在島嶼的東部和南部更早被發現，而在溢油污染事件多發的秋冬季節，保護區內島嶼更容易受到來自西北方向溢油的污染影響，季風會導致來自渤海中部及北部海域的油污更容易到達保護區，這可能是在保護區內島嶼的西海岸往往更早發現溢油的原因。③在表2 中給出了置信度為95%區域的風險分析表，從區域的分布來看，受主導風向、往復流的共同影響，在1－2 月以及9－12 月，保護區周圍的高概率區域向保護區北部延伸，分布面積較大，并在12 月達到最大。而在4－8 月保護區南部的近岸部分成為高概率區域，分布面積較小。高概率區域的覆蓋面積表現為春夏季節較小、秋冬季節較大的特征。

表2 溢油污染風險分析

圖6 月平均條件下長島保護區海域溢油溯源模擬結果的概率圖

4 結論

本文基于海洋溢油溯源數值模擬方法，評估月平均統計條件下長島保護區海域整體污染危害風險，分析可能對保護區造成危害的來源區域及其溢油對保護區造成污染影響的概率等信息。模擬結果表明，①受到長島保護區海域強潮流作用和主導風向的顯著影響，高概率區域的分布常年表現為東-西走向的分布態勢，并在夏季向東南方向延伸，最遠可延伸至山東煙臺市北岸；在冬季向西北方向延伸，最遠可延伸至渤海中部和北部。②高概率區域的區域面積表現為春夏季節較小、秋冬季節較大的特征。溢油污染的最大風險出現在冬季。