海冰是纬度较高海域发生的一种海洋灾害^[1-2], 每年冬季渤海及黄海北部均会受到海冰灾害的影响^[3]。海冰外缘线及分布面积的确定, 是研究海冰变化规律、确定海冰冰情、启动国家海冰灾害应急预案的一个重要指标。海冰的外缘线受风、浪、流的影响变化很快^[4], 快速、准确地识别海冰外缘线, 对船舶航行、海上工程等尤其重要^[5-6]。

遥感监测作为可以全天候、全天时的监测方式, 是海冰监测的有力手段^[7-10]。目前科研人员利用遥感影像进行了多方法的海冰要素提取研究, 如基于形态学^[11]、欧式距离^[12]、稀疏分解和改进马尔可夫随机场(Markov random filed, MRF)模型^[13]、神经网络^[14-15]、机器学习^[16-18]、阈值方法^[19]和深度学习^[20]等多种方法进行海冰分类、确定海冰外缘线, 并进行提取。

然而, 上述研究均对海冰边缘线提取不够精确, 并且大多研究的是局部海冰信息提取, 没有考虑业务化海冰监测对大范围或整个海区海冰信息提取的需求, 另外都没有研究提取海冰信息的具体面积和时间效率问题。一副遥感影像, 因海冰和海水反照率差异造成影像上冰水的灰度明显不同, 基于此能够进行海冰外缘线的判别、提取。将影像上的冰水进行分离、提取海冰外缘线, 在数字图像处理领域即为图像分割。图像分割的方法较多, 基于变分和偏微分方程方法的图像分割因提取信息精度高、效率高, 已成功应用在医学领域的影像提取、人脸及车牌的智能识别、道路和河道以及建筑的提取等多方面。基于变分水平集方法、无需重新初始化的快速算法进行海冰外缘线检测及提取, 目前尚没有这方面的研究报道。

本文尝试使用基于变分水平集方法和无需重新初始化的、分裂Bregman投影算法(split Bregman projection method, SBPM)进行海冰外缘线提取及海冰分布面积统计, 并提出业务化应用的解决方案。将得出的图像和量化结果与人工提取结果进行比较后发现, 在提取的海冰外缘线的精度和提取时间效率方面, 本文提出的方法具有结果精确稳定、速度快、不受人工干扰等优点。

1 研究区域和数据来源

辽东湾是我国海冰灾害的典型海域, 每年11月份开始出现海面结冰现象, 至次年3月上旬左右海冰融化。本文选择该海域2018年1月中旬至2月中旬海冰灾害最严重时期的五景遥感影像进行研究。

遥感监测海冰分为可见光和合成孔径雷达(Synthetic Aperture Radar, SAR)两种方式, 两者在监测上各有优势。为满足业务化大范围监测海冰的需要, 本文选用五景中高时空分辨率、高幅宽的光学遥感卫星监测数据, 如中分辨率成像光谱仪(Moderate-resolution Imaging Spectroradiometer, MODIS)数据、环境与灾害监测预报小卫星B星(简称HJ-1B, 下同)数据和高分辨率对地观测四号卫星(简称“高分四号”, 下同)数据进行仿真实验。数据技术指标见表 1, 具体数据情况见表 2。

2 研究方法 2.1 遥感影像预处理方法

原始的遥感影像通常无法直接分析解译, 必须经过一些预处理工作。本文主要进行了图像增强、图像镶嵌、图像裁剪等预处理工作。因陆地、云等是无关区域, 若参与分割将降低运算效率, 因此对无关区进行了不规则图像裁剪工作。遥感影像预处理工作均基于美国ESRI公司的ENVI软件(以下简称ENVI)。实验区域遥感影像图及不规则影像裁剪结果, 如图 1、图 2所示。

2.2 变分水平集Chan-Vese模型

Chan-Vese模型^[21]是基于Mumford-Shah模型^[22] (M-S模型)的著名模型。M-S模型因最小轮廓线和光滑图像的维度不一致, 导致该模型无解; Chan和Vese对M-S模型提出了一个近似解, 即两相图像分割的Chan-Vese模型^[21](C-V模型), 定义如下:

$ \begin{array}{l} \begin{array}{*{20}{c}} {\rm{Min}} \\ {u, \phi } \end{array}\left\{ {E\left( {u, \phi } \right) = {\alpha _1}\int_{\mathit{\Omega} 1} {{{\left( {{u_1} - f} \right)}^2}H\left( \phi \right){\rm{d}}x} + } \right. \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\alpha _2}\int_{\mathit{\Omega} 2} {{{\left( {{u_2} - f} \right)}^2}} \left[ {1 - H\left( \phi \right)} \right]{\rm{d}}x + \\ \left. {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\gamma \int_\mathit{\Omega} {\left| {\nabla H\left( \phi \right)} \right|{\rm{d}}x} } \right\}, \\ \end{array} $

(1)

约束条件:

$ \left| {\nabla \phi } \right| = 1. $

(2)

式(1)中, f表示输入源图像, u=(u₁, u₂), u₁、u₂分别表示图像的前景和背景的图像灰度均值, ϕ表示水平集函数, H(ϕ)为Heaviside函数, 见式(3), α₁、α₂及γ表示惩罚参数, $\left| {\nabla H\left( \phi \right)} \right|$表示长度项。式(2)为式(1)的约束条件。

$ H\left( \phi \right) = \left\{ {\begin{array}{*{20}{c}} {1, }&{\phi \geqslant 0} \\ {0, }&{\phi < 0} \end{array}} \right.. $

(3)

对式(1)求解, 可以利用交替优化方法, 即ϕ固定, 分别优化u₁和u₂:

$ {u_1} = \frac{{\int_\mathit{\Omega} {fH\left( \phi \right){\rm{d}}x} }}{{\int_\mathit{\Omega} {H\left( \phi \right){\rm{d}}x} }}, {u_2} = \frac{{\int_\mathit{\Omega} {f\left[ {1 - H\left( \phi \right)} \right]{\rm{d}}x} }}{{\int_\mathit{\Omega} {\left[ {1 - H\left( \phi \right)} \right]{\rm{d}}x} }} , $

(4)

u₁、u₂固定, 优化ϕ:

$ \begin{array}{l}\begin{array}{c}{\rm{Min}}\\ \varphi \end{array}\{E\left(\varphi \right)={\displaystyle {\int }_{\mathit{\Omega} }\left[{\alpha }_{1}{\left({u}_{1}-f\right)}^{2}\text{ }-{\alpha }_{2}{\left({u}_{2}-f\right)}^{2}\right]H\left(\varphi \right){\rm{d}}x}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\gamma {\displaystyle {\int }_{\mathit{\Omega} }\left|\nabla H\left(\varphi \right)\right|{\rm{d}}x}\}.\end{array} $

(5)

令R(u₁, u₂) = α₁(u₁ – f)² – α₂(u₂ – f)², 对式(5)求变分, 可得对应的Euler-Lagrange方程, 得到梯度下降流:

$ \frac{\partial \varphi }{\partial t}=\text{ }\left[\gamma {\rm{div}}\left(\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\right)-R\left({u}_{1}, {u}_{2}\right)\right]\delta \left(\varphi \right) . $

(6)

式(6)可利用有限差分方法得到近似解。

2.3 分裂Bregman投影算法(Split Bregman Projection Method, SBPM)

C-V模型在曲线演化过程中为保持符号距离函数特征, 需要定义一约束条件(式2), 解决办法通常需要重新初始化, 导致计算效率低下。GoldStein等人提出采用SBPM算法^[23]提高计算效率, 即引入一个辅助变量$\overrightarrow w $, 通过下述计算过程代替原变分水平集的求解。

对式(1), 令R(u₁, u₂) = α₁(u₁ – f)² – α₂(u₂ – f)², 并引入一辅助分裂变量$\overrightarrow w $, 令$\overrightarrow w = \nabla \phi $, 因$\int_\mathit{\Omega} {\left| {\nabla H\left( \phi \right)} \right|{\rm{d}}x} $与$\int_\mathit{\Omega} {\left| {\nabla \phi } \right|\delta \left( \phi \right){\rm{d}}x} $等价, 则式(1)与式(2)可改写为:

$ \begin{gathered} \begin{array}{*{20}{c}} {\rm{Min}} \\ {\phi , \vec w} \end{array}\left\{ {E\left( {\phi , \vec w} \right) = \int_\mathit{\Omega} {RH\left( \phi \right){\rm{d}}x} + \gamma \int_\mathit{\Omega} {\left| {\vec w} \right|\delta \left( \phi \right){\rm{d}}x} + } \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\frac{\theta }{2}\int {{{\left( {\vec w - \nabla \phi - \vec b} \right)}^2}{\rm{d}}x} } \right\}, \hfill \\ \end{gathered} $

(7)

约束条件:

$ \vec w = 1. $

(8)

式中, $ {\overrightarrow{b}}^{k+1}={\overrightarrow{b}}^{k}+ \nabla {\varphi }^{k}- {\overrightarrow{w}}^{k} $为Bregman迭代参数, θ为正的惩罚参数。

对(7)式可采用交替优化法分别得到$\phi $和$\overrightarrow w $:

$ \left\{ \begin{array}{l} R\left({u}_{1}, {u}_{2}\right)\delta \left(\varphi \right)+\gamma \left|{\overrightarrow{w}}^{k}\right|\frac{\partial \delta \left(\varphi \right)}{\partial \varphi }-\theta \nabla \cdot \left(\nabla \varphi +{\overrightarrow{b}}^{k+1}-{\overrightarrow{w}}^{k}\right)=0\; \;\;\;\;\;\;\;\;\; \text{在} \mathit{\Omega} \text{内部}\\ \left(\nabla \varphi +{\overrightarrow{b}}^{k+1}-{\overrightarrow{w}}^{k}\right)\cdot \overrightarrow{n}=0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \text{在}\partial \mathit{\Omega} \text{之上}\end{array} \right. , $

(9)

$ \left\{ \begin{array}{c}\gamma \frac{{\overrightarrow{w}}^{k}}{\left|{\overrightarrow{w}}^{k}\right|}\delta \left({\varphi }^{k+1}\right)+\theta \left({\overrightarrow{w}}^{k}-\nabla {\varphi }^{k+1}-{\overrightarrow{b}}^{k+1}\right)=0 \\ \text{约束条件}:\left|\overrightarrow{w}\right|=1\end{array} \right. . $

(10)

式(9)可采用半隐式高斯赛德尔(semi-implicit Gauss Seidel iteration)迭代方法进行求解, 得到式(11); 为提高计算效率, 式(10)上面的方程可采用解析广义软阈值公式进行求解, 得到式(12), 下面的约束公式可采用投影公式进行计算, 得到式(13):

$ {\varnothing}_{i, j}=\text{ }\frac{-\left[R\left({u}_{1}, {u}_{2}\right)-\theta \nabla \cdot {\overrightarrow{b}}^{k+1}+\text{ }\theta \nabla \cdot {\overrightarrow{w}}^{k+1}-\theta \alpha \right]}{4\theta } , $

(11)

$ \begin{array}{l} {\overrightarrow w ^{k + 1}} = {\rm{Max}}\left[ {\left| {\nabla {\varnothing ^{k + 1}} + {b^{k + 1}}} \right| - g\frac{\gamma }{\theta }\delta \left( {{\phi ^{k + 1}}} \right){\text{, }}\;0} \right] \hfill \\ \;\;\;\;\;\;\;\;\;\; \cdot \frac{{\nabla {\varnothing ^{k + 1}} + {b^{k + 1}}}}{{\left| {\nabla {\varnothing ^{k + 1}} + {b^{k + 1}}} \right|}}, \hfill \\ \end{array} $

(12)

$ {\overrightarrow{w}}^{k+1}=\text{ }\frac{{\overrightarrow{w}}^{k+1}}{{\left|{\overrightarrow{w}}^{k+1}\right|}_{\epsilon }} . $

(13)

式(11)中, α = ϕ_{i+1, j} +ϕ_{i–1, j} +ϕ_{i, j+1}+ϕ_{i, j–1}; 式(13)中ε为正的小量, 以防止分母为0溢出。

3 仿真实验与分析 3.1 仿真实验步骤

基于分裂Bregman算法, 快速、精确提取海冰信息的业务化应用实验, 主要步骤是: 首先选择能覆盖辽东湾、噪声少、清晰度高的卫星遥感图片(以下简称“卫片”)进行图像增强、不规则图像裁剪等预处理工作; 其次基于Matlab软件平台, 将上述变分及SBPM算法程序化, 针对不同的卫片调整初始参数, 当冰水分界线分割完好则停止程序计算; 最后对分割好的海冰结果进行量化, 得到最终的海冰分布总面积。按照计算机程序设计语言绘制流程图的原则, 绘制的工作流程图如图 3所示。

根据卫片的不同, 可以选择多种预处理方法。本文采用了图 3所示图像镶嵌等3种方法。为加快海冰信息提取的速度, 预处理卫片时, 将实验区域与陆地、岛屿、云等无关信息分离是一个关键步骤, 经过实验, 本文利用ENVI软件的感兴趣区(region of interest, ROI)工具, 通过裁剪不规则图像, 实现了实验区域与无关信息的分离。基于Matlab语言程序计算得到海冰外缘线后, 根据具体卫片的图像分辨率, 将提取的海冰像素点个数与单位像元的面积相乘即可完成所提取海冰总面积的量化工作(图 4)。

3.2 仿真实验结果与分析

所有的仿真实验工作在一台惠普Z620工作站上开展, 基于Matlab 2014a、ENVI4.5及ArcGIS10.2等软件平台完成。实验的数据情况见表 2。基于变分及SBPM算法的程序的几个主要初始参数设计为: α = 5, γ = 5, θ = 3 000, 迭代次数设计为15。

图 4是2018年2月4日MODIS遥感影像图。图 5是上述影像基于变分方法提取结果。为更清楚地查看提取结果, 在图 5上任意选取了4个区域1、2、3、4, 分别对应图 6中的(a)、(b)、(c)及(d)。图 7是2018年2月4日MODIS影像基于人工方法提取结果。

图 8是2018年2月6日“高分四号”卫星影像图。图 9是上述影像基于变分方法的提取结果。为更清楚地查看提取结果, 在图 9上任意选取了4个区域1、2、3、4, 分别对应图 10中的(a)、(b)、(c)及(d)。图 11是2018年2月6日“高分四号”影像基于人工方法的提取结果。

图 12是2018年1月29日HJ-1B遥感影像图。图 13是上述影像基于变分方法的提取结果。为更清楚地查看提取结果, 在图 13上任意选取了4个区域1、2、3、4, 分别对应图 14中的(a)、(b)、(c)及(d)。图 15是2018年1月29日HJ-1B影像基于人工方法的提取结果。

图 16是2018年1月12日HJ-1B遥感影像图。图 17是上述影像基于变分方法的提取结果。为更清楚地查看提取结果, 在图 17上任意选取了4个区域1、2、3、4, 分别对应图 18中的(a)、(b)、(c)及(d)。图 19是2018年1月12日HJ-1B影像基于人工方法的提取结果。

图 20是2018年2月15日“高分四号”卫星影像。图 21是上述影像基于变分方法提取结果。为更清楚地查看提取结果, 在图 21上任意选取了4个区域1、2、3、4, 分别对应图 22中的(a)、(b)、(c)及(d)。图 23是2018年2月15日“高分四号”影像基于人工方法的提取结果。

上面的图 5、6、9、10、13、14、17、18、21和22共10张图是基于变分及SBPM快速算法提取海冰信息的结果, 绿线即提取的海冰外缘线, 黑色及天蓝色背景(用PhotoShop软件处理过)是海水、陆地等非海冰区。在五景影像上分别选择了三处有代表性的区域进行重点、放大显示, 一是靠近陆地、有冰有水有陆地的区域(见图 5、9、13、17和21的1区); 二是靠近陆地有初生冰的区域(见图 5、9、13、17和21的3区); 三是远离陆地冰水分明的区域(见图 5、9、13、17和21的2区和4区)。从提取结果图片上看, 1区提取的效果较差, 原因是近辽东湾底海水结冰不明显, 冰水灰度接近, 变分法没有进行较好地识别; 3区的初生冰区域提取过细, 导致部分初生冰没有被识别; 2区和4区提取的效果最好, 海冰外缘线及冰间水道提取的精度都非常高。图 7、11、15、19及图 23五张图是基于人工解译方法提取的海冰信息结果, 辽东湾内的灰白色区域即为提取的海冰外缘线及分布面积。人工解译方法因为加入了专家知识, 在靠近陆地、冰水灰度不明显及初生冰较多的地方, 识别得较好, 但在远离陆地、冰水明显分离的冰间水道区域提取的不够精确。表 3是两种方法的数值比较结果, 其中的“提取时间”是指纯粹的海冰外缘线提取及海冰分布面积计算的时间。这个比较结果和图片结果基本吻合, 因人工解译方法对冰间水道及外缘线提取不精确, 造成数值结果大于变分方法; 从提取时间上看, 基于变分及SBPM方法最快仅需要15.568 9 s, 最慢也仅需要41.652 3 s, 而人工解译方法最快需要2 min 53 s, 最慢需要15 min 14 s, 在提取时间上变分方法有很大的优势。另外, 基于变分方法提取海冰外缘线结果稳定, 不受人为因素影响, 不会出现解译专家不同解译结果迥异的情况。

4 结论

本文提出基于变分水平集方法及快速分裂Bregman投影方法的海冰信息提取业务化工作解决方案。仿真实验结果表明, 基于变分水平集方法及快速分裂Bregman投影方法能够应用到业务化海冰监测工作中; 无论在提取海冰外缘线的精度上还是在提取海冰信息的时间上, 变分及SBPM方法都有很大的优势。