常宇佳,邱志伟,2*,牛 原,吴振宇, 潘春天
(1.江苏海洋大学 海洋技术与测绘学院,江苏 连云港 222000;2. 江苏省海洋资源开发研究院, 江苏 连云港 222000)
海岸线的检测对于沿海地区是非常重要的,对于地理地图、海岸侵蚀监测都有重大意义[1-5]。近年来,国内外专家学者对于SAR 影像海岸线检测提出了诸多方法,经过大量文献查阅,大致分为5 种。
1)区域型方法检测海岸线。该方法利用区域合并算法,通过连续自适应区域合并实现海岸线的分割。近些年,通过谢明鸿等对区域合并算法的不断优化,该算法已有效减少了计算量,并有较高的提取精度,但仍有缺陷,如海面上的高亮区域、地面与海面的相似区域,算法导致出现孤立的小区域,极大的影响了海岸线的分割效果。
2)阈值检测提取海岸线。采取合适的阈值实现海岸线的检测,只有当海岸线区域像素差异明显,并且相干斑噪声较小时,才能比较理想的提取出海岸线部分。该方法需要人为设置参数,并且对于阈值极为敏感,容错性不高。
3)马尔可夫随机场方法。马尔可夫随机场是建立在马尔可夫模型和贝叶斯理论基础之上的。自2012 年,Fabio 等提出基于马尔可夫随机场的海岸线检测算法,虽然该算法只能检测高分辨率单视点复合SAR 影像,但也可以提取出诸多海岸线的细节信息。之后Fabio 等又在该算法进行了改进,提出一种无监督,低报错率的海岸线检测算法。但是该方法需要大量的计算量,且提取精度不理想。
4)人工神经网络方法。该方法模拟人类视觉信息处理提高海岸线检测精度,有较好的学习能力和自动更新能力,但是在海面纹理数据差异较大时,该算法的处理能力相对较弱,且该算法耗时长,提取精度受样本数量和特征类型的影响。
5)主动轮廓模型法。该方法将海岸线区域与主动轮廓线相结合,给出基于海岸线检测区域的水平集法。该方法可得到连续的海岸分割线,但该算法的参数需要人工调试,且对于海面粗糙度不均匀的区域效果不是很理想。本文在众多理论基础上改进并提出一种适用于Sentinel-1 影像的海岸线自动提取方法,该方法快速、简便、提取精度较高,实用性较强。
连云港是全国为数不多拥有港口的城市,海岸线全长176.5 km,因此对连云港的海岸线进行监测分析有着非常重大的意义。本次实验区域为接近连岛附近的海岸线,全长约为32 km。Sentinel-1 影像分辨率为30 m,重访周期为12 d。
本文使用的海岸线提取流程如图1 所示,首先获取连云港SAR 多年的影像,然后进行图像的预处理,将其相对应的轨道数据与之结合,处理成清晰的影像[6],接着运用MATLAB 软件使用K-means 算法分割处理,然后进行形态学处理,最后运用canny 算子提取海岸线[7],最终根据多年海岸线,进行时空变化分析及精度验证。
图1 海岸线提取流程图
对于Sentinel-1 合成孔径雷达影像,首先需要对图像进行图像预处理,利用ENVI SARscape 进行简单的预处理,然后进行降噪处理,否则图像中的一些无关噪声点可能会影响海岸线的提取[8]。因此,本文在尽可能保证海岸线信息完整的基础上,消除如船舶、海浪等噪声点的干扰。
降噪结束后,对图像进行NDWI 水体增强处理,根据水体指数直方图,选取合适的阈值,最后得到水体指数影像,如图2 所示。
图2 水体指数影像
本文所采用的二值化是基于Kittler 改进的二值化方法,本文所采用的最小错误准则公式为:
对于给定的p(g)和p(g|i),存在一个灰度g,其灰度g满足:
式中,t为最小错误阈值。通过以下改进公式得到最优解,该公式可以更简单的计算出准则函数,大大节省了时间。
通过改进的二值化方法得到二值化图像3c,并与Otsu 二值化3a 和kittle 二值化3b 作比较,如图3 所示,明显可以看出本文的二值化方法效果更好一些。
图3 3 种二值化方法比较
在降噪结束后,对图像进行K-means 聚类分割处理。在众多的聚类算法中,K-means 聚类分割是最简单的聚类算法之一,它属于一种无监督的学习算法。Kmeans 算法的难点在于需要设置不同的k值来获得不同的聚类结果,k值的不确定性正是该算法的一个缺点。为了达到好的实验结果,需要进行多次尝试才能够选取最优的k值。在本次实验中,选定的K值为10。
形态学处理是基于图像形状对图像进行的一系列处理操作,基本操作步骤包括:腐蚀-膨胀-填充-开运算-闭运算。
本文采用的算法是先对图像填充小的孔洞,补充小的像素点,使海岸线更加的光滑,然后进行去噪操作。接下来进行对图像膨胀操作,填补图像的高亮区域,然后进行腐蚀操作,将海岸线分解成像素元素,把不属于海岸线的元素变为个体,然后再进行填充,使海岸线变得更加光滑、连续。
边缘检测是图像处理和计算机视觉中的基本问题,边缘检测的目的是标识数字图像中亮度变化明显的点[9],就是确定像素点的状态来判断像素点是否在边缘。
常见的边缘检测算子有Roberts 算子、Prewitt 算子、Sobel 算子、Canny 算子。本文经实验研究实例对比,Canny 边缘检测算子为最优检测算子,选择Canny边缘检测算子对Sentinel-1 图像进行边缘检测。
Canny 算子在处理边缘检测时,首先利用高斯滤波器平滑图像,然后利用一阶偏导有限差分计算梯度幅值和方向,对梯度幅值进行非极大值抑制,用双阈值检测和连接边缘[10-12]。对形态学处理之后的图像进行Canny 算子检测边缘,在此,本文对算法做了进一步的优化,将图像在matlab 中提取出的笛卡尔坐标直接在算法中进行转化,使其转换为WGS84 坐标。
本文实验首先从欧空局获取原始图像,并用Envi进行图像预处理,裁剪后的结果如图4 所示。
图4 裁剪后的实验区域图片
使用算法进行图像处理,结果如图5 所示。图5展示了实验区域海岸线提取的一个实例过程,从图中可以看出图像经过K-means 算法分割(a)、二值化降噪处理(b)、形态学处理(c)、边缘检测算子检测(d)等操作后,海岸线准确的被提取出来,最后与原始岸线做对比(e)。从图5 中可看出经形态学处理后岸线分割完整,提取的海岸线与原本岸线基本完全契合。
图5 实验区域海岸线提取过程
图6 为多年海岸线提取对比图,实验数据采用不同年份相同季度的SAR 影像。从实验结果来看,2016 年与2015 年岸线变化并不明显,而自2017-2020 年,岸线变化逐年明显,向内陆延伸明显。为了验证预测,随后在获取到提取的海岸线数据后,可以计算出所研究区域的海岸线包围的陆地面积,从而计算出陆地变化的趋势。提取出的海岸线离散点的坐标,把图像边界像素点也加进散点数据集,所有的点连接成一个任意多边形,然后利用多边形面积公式计算出陆地的像素面积。公式如下:
图6 多年海岸线提取对比图
式中,xi、yi为海岸线离散点的坐标,计算出像素面积后,则可通过比例关系,计算出陆地真实的面积。
经过图像处理并提取数据,计算出实验区域连云港海岸陆地面积变化时间在2015-2020 年,如图7 所示。表中表明,自2015 年,连云港内陆面积在逐年减少,可能是由于近年来全球气候变暖,海平面升高以及多年海水侵蚀所致。经实验数据表明,实验区域连云港沿岸陆地面积年平均缩减率为4.209 km2/a。由此例可以看出,本文提出的海岸线提取算法是可以基于哨兵一号数据对海岸线数据进行提取的。因此进行了分析和计算,极大的避免了传统野外测量效率低、周期长的弊端,证实了该方法的在海洋工程方面的可行性和有效性。
图7 2015-2020 年实验区域海岸陆地面积变化
以本文的精度验证进行实地测量,测量时间为2020-05-17,并下载2020-05-23 的Sentinel-1 数据。本文选取30 个特征点,其中最大像素偏差为5,最小像素偏差为1,像素平均偏差为2.9,如图8 所示。从精度对比结果分析得出结论,海岸线提取精度较高。
图8 像素偏差直方图
本文首先获取了连云港多年同一季度SAR 影像,并进行一系列的数据处理操作,然后应用Matlab 软件对多年影像进行海岸线自动提取分析。实验结果显示,利用该方法对Sentinel-1 影像提取边界线效果良好,从中可以提取出准确海陆边界线,对生态环境保护、农林牧渔业的发展、海洋资源的规划管理及航海安全等方面有着重要的影响,并为基于SAR 影像提取边界线进行时空变换分析做出了补充。