一种海面溢油预测方法与流程

本发明属于海面油污扩散漂移预测技术领域，特别是一种误差小、准确度高的海面溢油预测方法。

背景技术：

当海面发生溢油事故后，为对油污进行有效清理和进行污染预防，提前判断溢油的运动轨迹和归宿至关重要。然而，因于环境条件多变，随着时间推移，对其运动轨迹和归宿难度增加。

为了预测溢油发生后随着时间和环境条件改变后其运动轨迹和归宿，更好地为海上清污防治的资源调度和配置提供技术支持和决策依据，国内外学者对溢油行为的数值模拟做了大量研究，并建立了相关的溢油预测模型。其大致可分为油膜扩展模式、对流扩散模式和“油粒子”模式3类。

扩展模式中Fay的经典三阶段扩展理论、Blokker扩展模式和刘肖孔公式受到广泛重视。

对流扩散模式包含了漂移和离散过程，漂移模型有Webb等人建立的美国海军Navy模型、Williams等人建立的SEADOCK模型以及Delaware模型等；离散模型大多采用蒙特卡罗方法。

自从Johansen、Elhot等提出了油粒子概念后，“油粒子”模型得到了长足发展，常用确定性方法模拟平流过程，用随机方法模拟对流扩散过程，油粒子模型在油膜扩展中的研究应用较少。

然而，由于对各过程间的相互影响与联系研究较少，导致溢油预测存在较大误差。

技术实现要素：

本发明的目的在于提供一种海面溢油预测方法，误差小、准确度高。

实现本发明目的的技术解决方案为：一种海面溢油预测方法，包括如下步骤：

(10)溢油图像获取：灰度化实时海洋图像，并对灰度化后的实时海洋图像进行阈值分割和二值化处理，得到初始溢油图像；

(20)海面溢油检测：根据初始溢油图像，并依据实时海洋图像的的分辨率和经纬度，得到包括溢油区域方位与面积的溢油检测图像；

(30)海面溢油预测：以溢油检测图像为基础，根据海面溢油行为预测综合模型，计算得到包括预测溢油区域方位与面积的溢油预测图像。

本发明与现有技术相比，其显著优点为：

误差小、准确度高：传统的溢油预测模型大多是对扩展、漂移和风化各变化模块独立的计算模型研究，没有系统地将各运动模块组合，实现对溢油发生后整个行为的预测计算。且较少的考虑各模块之间的相互影响与联系，预测结果可能会产生较大的误差。本发明总结的计算模型中将各变化模块组合成整体，并考虑各模块之间的相互影响与联系，增大预测结果的准确性；

下面结合附图和具体实施方式对本发明作进一步的详细描述。

附图说明

图1为本发明海面溢油预测方法的主流程图。

图2为图1中海面溢油预测步骤的流程图。

图3为图2中溢油风化预测步骤的流程图。

图4为不同油品的密度、API以及摩尔体积。

图5为实验室水槽实验结果与模拟结果的比较示意图，包括初始面积、结束面积、运动长度和运动时间的比较。

图6(a)为水槽实验初始溢油红外图像，图6(b)为水槽实验溢油发生10s后的红外图像。

图7(a)为MATLAB模拟溢油检测图像，图7(b)为MATLAB模拟溢油预测图像。

图8(a)为溢油区域图像，图8(b)为经过灰度化、阈值分割和二值化处理的溢油区域图像。

图9为渤海地图，用作溢油背景海洋图像。

图10为初始溢油检测图像。

图11为溢油预测图像

图12(a)为以渤海为背景的初始溢油检测图像，图12(b)为以渤海为背景的溢油预测图像。

图13为基于MATLAB GUI的交互式预测模拟平台示意图，包括参数设置、检测结果显示、预测结果显示和显示模式按键。

图14为原油在西南风、风速5米/秒、洋流方向向东、洋流速度为4米/秒、不同时间下的预测轨迹图像，(a)～(d)分别表示溢油发生后1小时、3小时、5小时以及10小时的预测图像。

图15原油在西南风、风速10米/秒、洋流速度5米/秒，不同洋流方向的环境条件下溢出后1小时的预测图像，(a)～(d)分别表示洋流方向向东、向南、向西以及向北的预测图像。

图16不同品种的油在西南风、风速为10米/秒、洋流速度为5米/秒、洋流方向朝东的环境条件下溢出后1小时的预测图像，(a)～(d)分别表示原油、车用汽油、轻柴油以及沥青的预测图像。

图17(a)表示原油溢出后1小时，在西南风、风速为3米/秒、洋流速度为4米/秒、洋流方向向东的环境条件下的预测图像，

图17(b)表示原油溢出后3小时，在西风、风速为5米/秒、洋流速度为4米/秒、洋流方向向东的环境条件下的预测图像，

图17(c)表示原油溢出后5小时，在西北风、风速为8米/秒、洋流速度为5米/秒、洋流方向向东的环境条件下的预测图像，

图17(d)表示原油溢出后10小时，在北风、风速为10米/秒、洋流速度为5米/秒、洋流方向向东的环境条件下的预测图像。

具体实施方式

如图1所示，本发明海面溢油预测方法，包括如下步骤：

(10)溢油图像获取：灰度化实时海洋图像，并对灰度化后的实时海洋图像进行阈 值分割和二值化处理，得到初始溢油图像；

读入溢油区域的图像，如图8(a)所示，将溢油图像转换成灰度图像，利用最大类间方差法对灰度图像进行阈值分割处理。最大类间方差法(OTSU)是1979年由日本学者Nobuyuki Otsu提出的一种自适应确定阈值的方法。依据图像灰度值的不同，将其分为背景和目标两部分。背景和目标的类间方差越大，表示组成图像的两部分区别越大，若将目标错分为背景或背景错分为目标则导致两部分区别减小。因此，类间方差最大时分割错误的概率最小。对于图像I(x，y)，计算原理如下所示：

$p_1=\frac{n_1}{m\×n}---\left(1\right)$

$p_2=\frac{n_2}{m\×n}---\left(2\right)$

n₁+n₂＝m×n (3)

p₁+p₂＝1 (4)

a＝p₁a₁+p₂a₂ (5)

S＝p₁(a₁-a)²+p₂(a₂-a)² (6)

将式(5)代入式(6)得到等价公式：

S＝p₁p₂(a₂-a₁)² (7)

使用遍历方法得到类间方差最大时的阈值，完成图像的阈值分割。其中Th表示目标和背景的分割阈值；p₁表示目标的像素数占总图像的比例，a₁为其平均灰度；p₂表示背景像素数占总图像的比例，a₂为其平均灰度；a表示总图像的平均灰度，S表示类间方差。图像大小记作m×n，n₁表示图像中灰度值小于Th的像素数，n₂表示灰度值大于Th的像素数。对阈值分割后的图像进行二值化处理，将灰度值小于阈值Th的像素点的灰度值设置为0，灰度值大于阈值Th的像素点的灰度值设置为255，如图8(b)所示，便 于后期的溢油面积和方位的检测计算。

(20)海面溢油检测：根据初始溢油图像，并依据实时海洋图像的的分辨率和经纬度，得到包括溢油区域方位与面积的溢油检测图像；

利用Google Earth软件查询中国渤海地图，并按一定比例尺截图作为本次研究的背景，如图9所示。依据图9的分辨率，做出相同分辨率大小的纯黑背景图片(比例尺和位置信息与图9也相同)，然后将油污图像缩小放至纯黑背景中，如图8所示，以图8为初始溢油图像对其进行检测和预测。对图像8进行扫描，并根据图中的比例尺、位置信息和分辨率计算油污的面积和方位。

(30)海面溢油预测：以溢油检测图像为基础，根据海面溢油行为预测综合模型，计算得到包括预测溢油区域方位与面积的溢油预测图像。

以初始溢油检测图像为基础，根据上述海面溢油行为预测综合模型进行计算，得到溢油预测图像，并计算预测溢油区域的面积和方位。利用上述预测模型对溢油进行预测计算，得到预测图像和相关面积、方位，如图11所示。将油污图像与背景图像叠加，形成更为直观的油污检测图像，如图12所示。

如图2所示，所述(30)海面溢油预测步骤包括：

(31)溢油扩展预测：溢油自身的扩展主要受油膜自身的组成成分和性质影响，因此溢出后风化引起的性质变化对扩展有一定影响。此外，表面风速对扩展的影响也不可忽略。扩展过程可根据Fay的三段扩展改进理论。

扩展各阶段油膜的长轴l、短轴r随时间的变化按下式计算：

重力扩展阶段：

r₁＝K₁(ΔgV)^1/4t^1/2， (8)

l₁＝r₁+cu_f^ηt^ε， (9)

粘性扩展阶段：

$r_2=K_2{\left(\frac{{\mathrm{\ΔgV}}^2}{\sqrt{v_w}}\right)}^{1/6}{t}^{1/4},---\left(10\right)$

l₂＝r₂+cu_f^ηt^ε， (11)

表面张力扩展阶段：

$r_3=K_3{\[\frac{{\mathrm{\δ}}^2}{{{\mathrm{\ρ}}_W}^2{v}_w}\]}^{1/4}{t}^{3/4},---\left(12\right)$

l₃＝r₃+cu_f^ηt^ε， (13)

式中：

Δ＝1-ρ₀/ρ_w， (14)

ρ₀、ρ_w分别表示油、水的密度，；g表示重力加速度，V表示溢油体积，；ν_w为水的运动粘滞系数，表征液体反抗形变的能力，是液体在流动中产生能量损失的主要原因，可由公式

ν_w＝0.01775/(1+0.0337T+0.000221T²)， (15)

求得，T表示水温，T＝10℃时，ν_w＝1.307×10^-6，δ为净表面张力系数，取0.0308；，t表示溢油时间，K₁、K₂、K₃表示各阶段扩展系数， K₁＝2.28，K₂＝2.90，K₃＝3.20；u_f表示风速，c表示经验常数，取c＝0.03；η＝4/3；ε＝3/4；

(32)溢油漂移预测：漂移扩散主要由平流过程和紊动扩散过程两部分组成。平流过程是表面流和风力作用的结果，其中最为主要的影响因素为洋流速度及其方向。

漂移速度由下式计算得到：

其中，

式中，是油膜迁移速度矢量，是表面洋流速度，α_c是表面洋流漂移系数，取1.0；a_w是风漂流系数，a_w＝3.5%；是风漂流速度矢量，由下式给出，

式中，u₁₀是海面以上10m处的风速，α是风向角，α′是柯式偏转角，取15°；

当紊动扩散为各向异性时,油粒子在水平方向上的产生的紊动扩散速度为：

式中，R表示均值是0，标准差是1的正态分布随机数，K_x、K_y分别表示x、y方向上的扩散系数，取值10～50m²/S；

(33)溢油风化预测：综合溢油蒸发、乳化、密度，预测溢油风化；

油膜的风化过程，包括蒸发、氧化、乳化、溶解、生物降解和沉降。其中蒸发和乳化对溢油的性质影响较大，不可忽略。

如图3所示，所述(33)溢油风化预测步骤包括：

(331)蒸发系数计算：蒸发是溢油风化过程中最主要的部分，受溢油的密度、API和摩尔体积等因素影响。根据实际情况选择Mackay等人提出的解析法计算蒸发系数。

按如下解析法计算蒸发系数，

$F_V=\[\ln P_0+\ln ({\mathrm{BK}}_{E}t+\frac{1}{P_0})\]/B,---\left(21\right)$

K_E＝K_MAV_M/GTV₀， (22)

K_M＝0.0025u₁₀^0.78， (23)

式中，F_V是蒸发系数，t是时间，A是油膜的面积，V_M是摩尔体积，取值150×10^-6～600×10^-6m³/mol，G是气体常数，取值8.206×10^-5，T是油的表面温度，约等于大气温度T_E，V₀是溢油的初始体积，P₀是初始挥发气压，关系式如下：

$\ln P_0=10.6(1-\frac{T_0}{T_E}),---\left(23\right)$

式中，T₀是初始沸点，T_E＝283K时，B、T₀的值可由如下公式计算：

B＝1158.9API^-1.1435， (24)

T₀＝542.6-30.275API+1.565API²-0.03439API³+0.0002604API⁴，(25)

式中，API表示美国石油协会采用的比重；

不同油品的密度、API以及摩尔体积的具体数值如图4所示。

(332)含水率计算：乳化受表面风速的影响，Mackay等提出用含水率来标志乳化程度，

按下式计算标志乳化程度的含水率：

$Y_W=\frac{1}{K_B}(1-e^{-K_A{K}_B{(1+u_f)}^2}t),---\left(26\right)$

式中，Y_W为含水率，表示乳化物中的含水量(％)，

K_A＝4.5×10^-6， (27)

$K_B=\frac{1}{Y_W^F}\≈1.25,---\left(28\right)$

表示最终含水量，取0.8，u_f为表面风速，t为溢油时间；

(333)溢油密度计算：溢油密度由于风化过程而不断变化，主要考虑乳化和蒸发的影响。综合两者的影响。

溢油密度如下式，

ρ＝(1-Y_W)[(0.6ρ₀-0.34)F_V+ρ₀]+Y_W·ρ_w， (29)

式中，ρ表示蒸发后油的密度，ρ₀表示油的初始密度，ρ_w表示海水密度。

(334)体积变化计算：蒸发系数F_V为蒸发量与油总量的比值。

蒸发系数对体积的影响由下式计算，

V＝(1-F_V)^t·V₀， (30)

(34)位移计算：设油膜某一质点在t_i时刻的坐标为S(t_i)，t_i+1时刻的坐标为S(t_i+1)，则有

为扩展、漂移矢量在Δt时段内位移的合成，即

为扩展在Δt时段内的位移矢量，即

为漂移在Δt时段内的位移矢量，即

蒸发和乳化会引起溢油的密度和体积的变化，从而对油膜扩展的长轴与短轴造成一定的影响，为计算准确，将公式(29)、(30)代入扩展公式(8)～(13)计算扩展的长轴与短轴，再代入公式(33)计算扩展的位移矢量。

为便于直观观察和简化操作，利用MATLAB GUI设计交互式界面，如图13所示，界面顶部为参数设置栏，根据影响因素分为扩展项和漂移项和时间参数。扩展项包括溢油种类，可通过下拉菜单进行选择，已有种类包括原油、车用汽油、轻柴油和沥青四种，不同油种的密度、API和摩尔体积均有所不同，具体数值如图4所示。漂移项包括风速、风向、洋流速度、洋流方向。风速和洋流速度可在编辑框中直接输入数值。风向和洋流方向可分别通过下拉菜单选择。其中风向包括南风、西南风、西风、西北风、北风、东北风、东风和东南风。洋流方向包括向北、向东北、向东、向东南、向南、向西南、向西和向西北八个方向。界面中部为图像显示区域，左边为溢油的初始图像，右边为溢油的预测图像，两幅图像下方分别为溢油初始面积和方位、预测溢油面积和方位。右下角为动态显示和静态显示按钮。选择动态显示时，只需设置溢油品种，其他条件随时间的 改变而改变，可实现溢油轨迹的动态模拟。选择静态显示则需要设置溢油品种、风速、风向、洋流速度、洋流方向和时间，实现特定条件下的溢油预测。

下面结合仿真实例对本发明做进一步的说明。

为了验证本发明预测模型的准确性，我们采用了水槽实验验证溢油的扩散漂移。实验以原油作为溢油样品，密度、API以及摩尔体积的具体数值见图4。水的密度为1000kg/m³，水温为25℃，则由公式

v_w＝0.01775/((1+0.0337T_W+0.000221T_W^2))

可得v_w为0.00896cm²/s。水槽宽30cm，长4m，取其流速均匀的一段59cm，水的平均流速约6cm/s，风速为0。利用红外摄像头采集水面溢油的运动轨迹，将采集到的初始溢油图像、溢油发生10s后的溢油图像作为实测结果，与相同条件参数下的模拟结果对比分析，如图5所示，红外实测结果与相应模型计算结果的对比，溢油扩展面积与漂移距离数值均较为接近。图6(a)为实测的初始溢油红外图像，图6(b)为10s之后的红外溢油图像，图7(a)为溢油检测图像，图7(b)为以图7(a)初始溢油图像通过模型计算得到的溢油预测图像。由图可知，实测结果与模型计算结果相符，验证了溢油预测模型的准确性。

为了证明本发明预测方法的简单、直观和快速，我们在拥有2.1GHz的CPU和4G内存的PC上利用MATLAB R2014a软件制作的预测模拟仿真平台完成了多组预测仿真实验。图13为基于MATLAB GUI的交互式预测模拟平台示意图，界面顶部为参数设置栏，根据影响因素分为扩展项和漂移项和时间参数。扩展项包括溢油种类，可通过下拉菜单进行选择，已有种类包括原油、车用汽油、轻柴油和沥青四种。漂移项包括风速、风向、洋流速度、洋流方向。风速和洋流速度可在编辑框中直接输入数值。风向和洋流方向可分别通过下拉菜单选择。界面中部为图像显示区域，左边为溢油的初始图像，右边为溢油的预测图像，两幅图像下方分别为溢油初始面积和方位、预测溢油面积和方位。右下角为动态显示和静态显示按钮。选择动态显示时，只需设置溢油品种，其他条件随时间的改变而改变，可实现溢油轨迹的动态模拟。选择静态显示则需要设置溢油品种、风速、风向、洋流速度、洋流方向和时间，实现特定条件下的溢油预测。

图14(a)为原油在西南风、风速5米/秒、洋流方向向东、洋流速度为4米/秒、溢油1小时后的预测轨迹图像；图14(b)为原油在西南风、风速5米/秒、洋流方向向东、洋流速度为4米/秒、溢油3小时后的预测轨迹图像；图14(c)为原油在西南风、风速5米/秒、洋流方向向东、洋流速度为4米/秒、溢油5小时后的预测轨迹图像；图14(d)为原油在西南风、风速5米/秒、洋流方向向东、洋流速度为4米/秒、溢油10小时后的预测轨迹图像。由图可知，溢油发生后，随着时间的推移，油膜的形状、大小和位置不断发生变化。短时间内，油膜迅速扩展，面积增长快速；后期由于厚度减小和蒸发等因素，面积增长减缓，直至油膜厚度减小到一定值扩展停止。

图15(a)为原油在西南风、风速10米/秒、洋流速度5米/秒，洋流方向向东的环境条件下溢出后1小时的预测图像；图15(b)为原油在西南风、风速10米/秒、洋流速度5米/秒，洋流方向向南的环境条件下溢出后1小时的预测图像；图15(c)为原油在西南风、风速10米/秒、洋流速度5米/秒，洋流方向向西的环境条件下溢出后1小时的预测图像；图15(d)为原油在西南风、风速10米/秒、洋流速度5米/秒，洋流方向向北的环境条件下溢出后1小时的预测图像。由图所见，溢油的漂移轨迹随着洋流方向的改变而产生明显的变化，可见洋流对预测溢油漂移起重要作用。

图16(a)为原油在西南风、风速为10米/秒、洋流速度为5米/秒、洋流方向朝东的环境条件下溢出后1小时的预测图像；图16(b)为车用汽油在西南风、风速为10米/秒、洋流速度为5米/秒、洋流方向朝东的环境条件下溢出后1小时的预测图像；图16(c)为轻柴油在西南风、风速为10米/秒、洋流速度为5米/秒、洋流方向朝东的环境条件下溢出后1小时的预测图像；图16(d)为沥青在西南风、风速为10米/秒、洋流速度为5米/秒、洋流方向朝东的环境条件下溢出后1小时的预测图像。从图可知，油膜扩展受溢油种类影响较大，溢油密度越小，扩展速度越快，反之，密度越大，扩展速度越慢。

图17(a)表示原油溢出后1小时，在西南风、风速为3米/秒、洋流速度为4米/秒、洋流方向向东的环境条件下的预测图像，图17(b)表示原油溢出后3小时，在西风、风速为5米/秒、洋流速度为4米/秒、洋流方向向东的环境条件下的预测图像，图17(c)表示原油溢出后5小时，在西北风、风速为8米/秒、洋流速度为5米/秒、洋流方向向东的环境条件下的预测图像，图17(d)表示原油溢出后10小时，在北风、风速为10米/秒、洋流速度为5米/秒、洋流方向向东的环境条件下的预测图像。由图17可得，风速对油 膜扩展有一定影响，风速越大，扩展越快；油膜扩展过程主要发生在前期，后期扩展作用减弱至无扩展；漂移主要受洋流运动的影响，海面风速和风向的影响较小。

图14～图16均为静态显示模式下的模拟仿真图像，参数设置完毕后，点击静态显示，等待10s左右时间即可获得溢油预测的仿真图像和面积、经纬度信息，操作简单，结果显示直观，仿真快速。