基于海洋溢油溯源数值模拟方法的长岛保护区污染风险分析

陈勤思，胡松，魏永亮

（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 结论

本文基于海洋溢油溯源数值模拟方法，评估月平均统计条件下长岛保护区海域整体污染危害风险，分析可能对保护区造成危害的来源区域及其溢油对保护区造成污染影响的概率等信息。模拟结果表明，①受到长岛保护区海域强潮流作用和主导风向的显著影响，高概率区域的分布常年表现为东-西走向的分布态势，并在夏季向东南方向延伸，最远可延伸至山东烟台市北岸；在冬季向西北方向延伸，最远可延伸至渤海中部和北部。②高概率区域的区域面积表现为春夏季节较小、秋冬季节较大的特征。溢油污染的最大风险出现在冬季。