基于KPCA光谱特征约束的水边线提取算法

陈飞宇，阮 鲲，胡友彬，曹 磊

CHEN Feiyu1,RUAN Kun2,HU Youbin1,CAO Lei3

1.解放军理工大学 气象海洋学院，南京 211101

2.南京师范大学 地理科学学院，南京 211101

3.江苏华高软件技术有限公司 大气海洋部，南京 211101

1.College of Meteorology and Oceanography,PLAUniversity of Science and Technology,Nanjing 211101,China

2.School of Geographic Science,Nanjing Normal University,Nanjing 211101,China

3.Department of Meteorology and Oceanography,Hhigh Software Technology（Jiangsu）Co.,Ltd.,Nanjing 211101,China

1 引言

水边线是潮汐波动下起伏不平的海面和陆地的瞬时交接线[1]。水边线信息提取对于后续的海岸线提取、岸线演变、影像定位[2]等研究具有重要意义。近些年来，遥感技术以其观测范围广、信息量大、获取信息快、更新周期短等优点[3-4]，逐渐成为水边线提取的重要方法和手段。

水边线提取实质上是遥感图像边缘的提取，关键技术是实现海陆分离[5]。除目视解译外，国内外学者相继提出了多种方法，大致可分为边缘检测方法、遥感分类方法和图像分割方法。边缘检测方法[6-8]主要利用遥感图像的边缘特征，使用边缘检测算子检测水体和非水体的边缘作为水边线。例如韩震等[9]以Landsat被动反射红外数据、热红外数据和合成孔径雷达（ERS-2 SAR）主动微波数据为信息源，分别采用Laplace-Gauss、Roberts和Sobel三种算法提取悬浮泥沙含量较高的长江口九段沙淤泥质滩涂水边线信息。马小峰等[10]使用Laplace-Gauss、Roberts、Prewitt、Sobel和Canny等算子分别对图像进行边缘检测实验，发现Canny算子对图像的提取效果最好，可以作为卫星图像中水边线自动提取的基本算法。张旭凯等[11]结合最小噪声分离变换、改进的归一化差异水体指数、数学形态学方法和改进的Canny边缘检测方法，以秦皇岛市海岸为例，采用高分辨率SPOT4卫星图像进行水边线提取。遥感分类方法主要基于遥感图像的光谱特征，通过各种遥感分类方法将图像分为水体和非水体，通过提取分类后水体区域和非水体区域的边界作为水边线。例如尧雪娟等[12]根据波谱相似度的思想，以水体独特的反射波谱曲线作为中心判据，将遥感图像中目标象元分为水体和陆地，然后利用形态学方法，对水陆分离后的图像进行海岸线提取。张汉女等[13]将福建海岸划分为五种典型海岸类型，并针对不同海岸类型的ETM+影像进行影像增强处理，然后采用SVM算法并选用Sigmoid核函数作为最佳分类器，使影像海陆边界更加清晰，有效提高了后续水边线提取的精度，最后通过目视解译提取水边线。邢坤等[14]首先使用模糊聚类和边缘跟踪的方法实现初始水边线轮廓的提取，然后使用改进的Snake模型提取最终的水边线。图像分割方法是近几年使用较多的一种方法，该方法通过水体与非水体的分割实现水边线的提取。图像分割的方法有很多，常用的方法有阈值分割、区域生长（种子点增长）法等。张华国等[15]利用IKONOS的近红外波段使用直方图法和阈值法来提取海岸线。翟辉琴等[16]针对水域的纹理比较细腻、灰度比较均匀的特点，提出了高、低帽变换区域生长法进行影像分割，实现对遥感影像上海洋的提取。Joo-Hyung Ryu等[17]根据TM4和TM5两个波段对陆地和水体反射率的不同，提出选择合适的阈值，采用密度分割的方法对韩国Gomso湾淤泥质海岸进行海岸线提取。另外，除了上述三种方法，一些结合高层先验知识和低层图像特征的分割方法也被应用于水边线的提取，其中，基于水平集方法的活动轮廓模型受到广泛关注。基于曲线演化理论与偏微分方程的活动轮廓模型利用能量泛函有效结合图像低层数据和高层信息（如目标边缘轮廓的连续性、闭合性、形状先验信息等）来完成分割，常见的算法包括测地线活动轮廓（Geodesic Active Contour，GAC）模型和基于距离正则化的水平集方法（Distance Regularized Level Set Evolution，DRLSE）。例如郭海涛等[18]提出了一种四叉树和测地线活动轮廓模型相结合的海陆影像分割方法，综合利用两种方法的优点，重构边界停止函数，演化水平集方程，实现海旅影像分割，取得了不错的效果。Yu等[19]以多种遥感影像作为数据源，首先使用不可分离小波变换检测得到初始水边线，然后通过DRLSE算法进行修正，获取最终的水边线。

在上述方法中，利用边缘检测方法提取出的边缘像素通常是不连续的，难以生成完整的水边线[20]，需要后续的边缘线段连接或使用活动轮廓模型方法进行轮廓线提取[21-22]。基于光谱特征分类的方法往往需要对所得的分类结果进行一些分类后处理，一般包括聚类分析、过滤处理和边缘处理等，过程复杂，精度不高。图像分割方法中，阈值法对于目标灰度不均匀的图像难以用统一的或局部的门限进行分割；区域生长法抗噪性能较差，种子点附近个别像元灰度值的突然变化极易造成水边线的偏移。活动轮廓模型融入了高层先验知识，对噪声不敏感，检测精度更高，但在遥感图像上水陆边界模糊、地形复杂多变的情况下，该方法的提取效果有待提高。另外，现有的模型和算法多停留在处理分割3个波段以下的图像，即仅针对灰度图像或者3个波段组成的彩色图像进行水边线提取，未能充分利用其他波段的信息。

针对以上方法存在的问题，为充分利用遥感影像的多光谱特征，进一步提高水边线的提取精度，本文提出了一种基于KPCA光谱特征约束的水边线提取算法。首先，利用KPCA提取水体样本的光谱特征，采用最大似然法估计特征空间中水体光谱特征概率密度函数的特征参数，进而构建水体的光谱特征项。其次，以测地线活动轮廓（Geodesic Active Contour，GAC）模型[23]为基础，建立图像数据项。然后，结合光谱特征项和图像数据项建立水边线提取模型。最后，为了验证模型的有效性，选取三个不同的实验区域分别进行水边线的提取，并从完整性、正确率、提取质量、运行时间、迭代次数等方面对提出的算法与未加入光谱特征的GAC模型和DRLSE算法进行比较。

2 基于KPCA的光谱特征项

令M={ψ1,ψ2,…,ψn}为遥感影像上水体样本的光谱矩阵，ψi(1≤i≤n)为m×1的列向量，代表第i个样本的多光谱矢量，m为波段数，n为水体样本数。利用KPCA方法对水体训练样本降维，定义映射函数为φ，采用高斯径向基函数（RBF）为核函数，其定义如下：得到样本协方差矩阵C的特征值λ和相应的特征向量V，满足V∈F{0},Cv=λv。

设测试样本矩阵为为m×1的列向量，代表第 j个样本的光谱矢量，m为波段数，t为测试样本数。测试样本在F空间向量Vk上的投影为：

其中为系数，k的选取是根据指定的各主成分方差累计的比例（在本文中是90%）确定的。

水体训练样本降到k维，在低维子空间中数据分布更紧凑并且更容易估计[24]。采用文献[25]提出的思想，假设核矩阵K′近似服从正态分布K′∼N(μ,Σ) ，μ 的最大似然估计为=0，Σ的最大似然估计为=Dk，其中，D为特征值矩阵，Dk由D的前k行k列组成，对角线上的特征值为 λ1,λ2,…,λk，满足 λ1≥λ2≥…≥λk。则水体光谱特征概率密度函数为：

根据贝叶斯公式，水体光谱特征概率密度函数即可度量影像中象元属于水体的后验概率，将其进行归一化。

为将光谱信息引入到水边线提取模型中，定义光谱特征项，用能量泛函表示为：

其中，ν为常数系数项，Rγ表示演化曲线包围的区域。通过变分法，可得到式（5）的梯度下降流为：

其中γ(t)为演化曲线，N为曲线的单位内法向量。

3 图像数据项

水平集分割模型中的图像数据项大多是基于边缘或基于区域的[26]。测地线活动轮廓模型（Geodesic Active Contour，GAC）是一种基于边缘的几何活动轮廓模型，其主要思想是利用黎曼空间中的测地线概念，把寻找图像中的边界线问题转化为寻找一条加权弧长最小值问题。该模型对噪声和伪边界有较强的鲁棒性，检测的边界可达到亚像素精度，并且模型所得到的光滑连续曲线，可弥补目标轮廓上的噪声、间隙以及其他不规则形状，非常适合于水边线这类具有复杂特征地物的轮廓提取。

设 γ(q):[0 ,1]→R2为平面参数曲线，I:[0 ,a]×[0,b]→R+为待检测图像，定义能量泛函：

其中，g(⋅)为停止函数。

其中，Gσ是方差为σ的高斯核函数，因此，g为关于图像梯度强度的递减函数。

计算能量泛函式（7）的变分，曲线演化的梯度下降流为：

其中，κ是欧式空间曲率，N是单位内法向量。

4 基于KPCA光谱特征约束的水边线提取模型

4.1 模型的建立

基于以上分析，定义KPCA光谱特征约束的水边线提取模型为：

其中，α∈[0,1]为光谱特征项的权重系数，ESPEC是光谱特征项，EGAC是图像数据项。由此可以看出，该模型综合利用了图像的底层特征和光谱特征。

由式（6）、（9）可知，式（10）的梯度下降流为：

其中，ε=-αν。

将曲线γ(0)隐式地表示为曲面的零水平集：{γ(0)=(x,y)|φ(x,y)=0}，这样就将曲线的演化转化为曲面的演化，曲线演化过程独立于曲线的参数化，因此可以自动处理演化过程中的拓扑结构的变化。考虑到数值计算的稳定性，一般用有符号距离函数来表示水平集，有符号距离函数的定义如下：

其中，D((x,y),γ)表示点(x,y)到曲线γ的最短距离，并规定曲线内部点的有符号距离函数值为负，外部点的有符号距离函数值为正。公式（11）的水平集演化方程如下：

4.2 模型的数值解格式

其中，i、j为图像坐标，n为迭代次数，Δt为迭代步长。

曲率κ采用中心差分格式离散化：

其中、和分别为：

4.3 算法流程

本文的水边线提取过程为：

（1）确定研究区域，选取遥感图像，并进行图像预处理。

（2）选取水体样本，进行KPCA变换，使用最大似然法估计水体概率密度函数的特征参数μ、Σ，进而构建光谱特征项。

（3）对遥感图像进行高斯滤波，计算梯度和停止函数，构建图像数据项。

（4）结合光谱特征项和图像数据项构建水边线提取模型，并利用水平集方法进行求解，输出水边线。

（5）对得到的水边线进行精度评价。

水边线提取流程图如图1所示。

图1 水边线提取的算法流程图

5 实验设计与结果分析

5.1 实验设计

为验证本文提出的水边线提取算法的正确性和有效性，选取三个不同区域的Landsat TM数据分别进行水边线的提取，并将实验结果与GAC模型和DRLSE算法的提取结果进行对比。

实验1中研究区域为浙江温州大门岛，其海岸类型是基岩海岸。影像成像时间为2014年12月6日。实验提取子图像(202×403)海岸带的水边线，所有波段以3×3为模板、σ=1进行高斯滤波。本文模型通过选择近海和近岸的水体样本565个，与TM影像的1～7共7个波段组成565×7的水体样本光谱矢量，高斯径向基函数参数，估计出水体类光谱特征概率密度函数的参数参数设置为ε=-5 000，时间步长t=3；为了获得较好的水边线提取结果，经过大量的实验，GAC模型时间步长t=1；DRLSE算法的参数设置为μ=1.5，λ=3，t=4，其中μ、λ分别为DRLSE算法中平衡能量项的权重系数[27]，t为时间步长。实验结果如图2所示。

实验2中的研究区域为浙江省温州市的霓屿岛，岛屿面积10.43 km2，岸线总长30.18 km，影像成像时间为2015年2月15日。与实验1相同，首先对遥感图像进行高斯滤波。其次，选取近海和近岸的水体样本363个，与影像的1～7共7个波段组成363×7的水体样本光谱矢量，高斯径向基函数参数，估计出水体类光谱特征概率密度函数的参数=0，=，其余参数同实验1。通过大量的实验，GAC模型时间步长t=0.8；DRLSE算法的参数设置为μ=3，λ=3，t=1。实验结果如图3所示。

图3 实验2结果图

实验3中的研究区域为福建省莆田市南日岛，其轮廓边界曲折复杂，影像成像时间为2015年4月13日。首先对遥感图像进行高斯滤波。其次，选取近海和近岸的水体样本374个，与影像的1～7共7个波段组成374×7的水体样本光谱矢量，高斯径向基函数参数估计出水体类光谱特征概率密度函数的参数=0，

其余参数同实验1。通过大量的实验，GAC模型时间步长t=0.8；DRLSE算法的参数设置为 μ=3，λ=3，t=6。实验结果如图4所示。

图4 实验3结果图

5.2 结果分析

如图2所示，图（a）是研究区大门岛6（R），5（G），3（B）波段合成的伪彩色图和初始演化曲线。图（b）是GAC模型的演化结果，可以看出GAC模型基本提取出了水边线的位置，但距离真实的水边线存在一定的距离，并且在一些细小的凹陷部位（如图（b）中的放大区域）提取不准确。图（c）是DRLSE算法的提取结果，相较于GAC模型，该算法的提取结果更平滑，更贴合实际的水边线，但部分区域存在过分割的情况（如图（c）中的放大区域）。图（d）是本文模型的提取结果，可以看出相比于上述两种模型，本文提出的方法在细小凹陷部位提取效果更好。

如图3所示，图（a）是研究区霓屿岛653（RGB）波段合成的伪彩色图和初始演化曲线。图（b）是GAC模型的提取结果，可以看出GAC模型提取的水边线距离岸线较远，并且放大区域的提取误差较大。图（c）是DRLSE算法提取的水边线，整体平滑，但部分区域由于过平滑导致尖角部分存在过分割的情况。图（d）是本文模型的提取结果，可以看出提取的水边线更加贴合实际岸线，并且尖角或者凹陷部位的水边线更加准确。

如图4所示，图（a）是研究区南日岛653（RGB）波段伪彩色合成图和初始演化曲线，相比于大门岛和霓屿岛，南日岛的岸线更加曲折复杂。图（b）是GAC模型的提取结果，可以看出，提取结果基本准确，但从放大区域的提取结果可以看出，仍存在演化不到位、地物分割不彻底等问题。图（c）是DRLSE算法提取的水边线，提取结果相比于GAC模型更加准确，但由图中放大的局部区域可以看出，部分区域的水边线存在着贴合度不够、过分割等现象。图（d）是本文模型的提取结果，通过对比可以看出，水边线提取效果更好，细节部位演化更彻底，结果更准确。

为了定量评价三种算法提取的水边线的准确性，将三种算法提取的水边线分别与人工解译的水边线进行比较，采用完整性(CP)、正确率(CR)、提取质量(QL)三个指标对结果进行评价，其计算方法为[28]：

式中，TN表示匹配的人工解译水边线长度，FN表示未匹配的人工解译水边线长度，TP表示匹配的提取水边线长度，FP表示未匹配的提取水边线长度。

另外，为了进一步分析算法性能，本文从运行时间和迭代次数等方面对三个模型进行了分析和比较，对比结果如表1所示。

表1 模型的定量比较

由表1可知，本文提出的方法相比于其他两种算法，提取的水边线在完整性、正确率和提取质量方面均最高，表明了算法的有效性。在运行时间方面，GAC模型需要的迭代次数和运行时间更多，本文方法次之，DRLSE算法最少。这是因为，GAC模型在利用水平集演化来获得水边线时，需要定期给水平集函数重新初始化为距离函数，而DRLSE算法引入了惩罚项，避免了重新初始化，大大提高了水边线的提取速度[27]。本文提出的方法与GAC模型的演化方式相同，但演化速度会随着式（13）中的项的变化而变化，在曲线的演化过程中，较大的会加速水边线向着边界演化，因此需要的演化次数和运行时间相对于GAC模型较少。

6 结论

本文提出了基于KPCA光谱特征约束的水边线提取算法，从遥感图像中提取近岸和近海的水体样本，构建光谱特征项。同时，利用GAC模型构建图像数据项，进而得到水边线的提取模型。Landsat TM数据集上的三个提取实验结果验证了算法的有效性，水边线的提取相比于GAC模型和DRLSE算法，在保证一定运行速度的同时，准确率更高。

参考文献：

[1]韩震，郭永飞.基于小波多分辨率分析提取长江口淤泥质潮滩水边线[J].海洋科学，2011，35（7）：67-70.

[2]郭海涛，申家双，黄辰虎，等.海岸带潮汐模型支持下的光束法区域网空中三角测量[J].测绘科学技术学报，2012，29（1）：33-37.

[3]王志一，徐素宁，姜艳辉，等.苏北废黄河三角洲岸线变迁与海岸冲淤动态遥感监测[J].南水北调与水利科技，2013（1）：136-139.

[4]Lillesand T M，Kiefer R W，Chipman J W，et al.遥感与图像解译[M].北京：电子工业出版社，2016.

[5]周良芬，何建农.弱交互式水边线提取新算法[J].计算机工程与设计，2013，34（4）：1442-1445.

[6]王方超，张旻，宫丽美.改进的Roberts图像边缘检测算法[J].探测与控制学报，2016，38（2）：88-92.

[7]Niedermeier A，Lehner S，Sanden J V D.Monitoring big riverestuariesusing SAR images[C]//Geoscienceand Remote Sensing Symposium，IGARSS’01，2001：1756-1758.

[8]Lee J S，Jurkevich I.Coastline detection and tracing in SAr images[J].IEEE Transactions on Geoscience& Remote Sensing，1990，28（4）：662-668.

[9]韩震，恽才兴.伶仃洋大铲湾潮滩冲淤遥感反演研究[J].海洋学报，2003，25（2）：58-63.

[10]马小峰，赵冬至，邢小罡，等.海岸线卫星遥感提取方法研究[J].海洋环境科学，2007，26（2）：185-189.

[11]张旭凯，张霞，杨邦会，等.结合海岸类型和潮位校正的海岸线遥感提取[J].国土资源遥感，2013，25（4）：91-97.

[12]尧雪娟.基于波谱特征曲线法的遥感图像海岸线提取[J].广东广播电视大学学报，2009，18（6）：104-107.

[13]Zhang Hannv，Jiang Qigang，Xu jiang.Coastline extraction using support vector machine from remote sensing image[J].Journal of Multimedia，2013，8（2）：175-182.

[14]Xing Kun，Fu Yili，Wang Shuguo，et al.Extraction of coastline in high-resolution remote sensing images based on the active contour model[J].Journal of Harbin Institute of Technology，2011，18（4）：13-18.

[15]张华国，黄韦艮，周长宝.应用IKONOS卫星遥感图像监测南麂列岛土地覆盖状况[J].遥感技术与应用，2003，18（5）：306-312.

[16]翟辉琴.基于数学形态学的遥感影像面状目标提取研究[D].郑州：解放军信息工程大学，2005.

[17]Ryu J H，Won J S，Min K D.Water extraction from landsat TM data in a tidal flat：A case study in Gomso Bay，Korea[J].Remote Sensing of Environment，2002，83（3）：442-456.

[18]郭海涛，孙磊，申家双，等.一种四叉树和测地线活动轮廓模型相结合的海陆影像分割方法[J].测绘学报，2016，45（1）：65-72.

[19]Yu Shujian，Mou Yi，Xu Duanquan，et al.A new algorithm for shoreline extraction from satellite imagery with non-separable wavelet and level set method[J].International Journal of Machine Learning and Computing，2013，3（1）：158-163.

[20]Liu H，Jezek K C.A complete high-resolution coastline ofAntarctica extracted from orthorectified radarsat SAR imagery[J].Photogrammetric Engineering&Remote Sensing，2004，70（5）：605-616.

[21]李林茹，高双喜，曹淑服.基于小波变换和梯度矢量流Snake模型的ERS-1 SAR图像的海岸线探测[J].河北工业科技，2004，21（4）：24-26.

[22]范典，郭华东，岳焕印，等.基于二进小波变换的SAR图像湖岸线提取[J].遥感学报，2002，6（6）：511-516.

[23]Caselles V，Kimmel R，Sapiro G.Geodesic active contours[C]//International Conference on Computer Vision，1995：61-79.

[24]Fang W，Chan K L.Incorporating shape prior into geodesic active contours for detecting partially occluded object[J].Pattern Recognition，2007，40（8）：2163-2172.

[25]Ruan K，Zha Y，Zhou Z，et al.A modified GAC model for extracting waterline from remotely sensed imagery[J].International Journal of Remote Sensing，2016，37（17）：3961-3973.

[26]沈霁，李元祥，周则明.参数自适应的KPCA先验形状约束目标分割[J].中国图象图形学报，2013，18（7）：783-789.

[27]Li C，Xu C，Gui C，et al.Distance regularized level set evolution and its application to image segmentation[J].IEEE Transactions on Image Processing，2010，19（12）.

[28]Wiedemann C，Heipke C，Mayer H，et al.Empirical evaluation of automatically extracted road axes[J].Empirical Evaluation Techniques in Computer Vision，1998：172-187.