一种基于聚类和边缘检测的遥感影像海岸线提取方法

技术领域

本发明涉及海岸线提取技术领域，具体涉及一种基于聚类和边缘检测的遥感影像海岸线提取方法。

背景技术：

海岸作为一种较为特殊的地类，对海岸线的提取及研究，可为沿海地区的资源开发、规划利用、生产生活、生态保护等提供宝贵的数据支撑及决策依据。遥感影像具有直观性、范围广、准确度高的特征，是用于研究海岸线提取的优良选择。目前，基于遥感影像的海岸线提取研究主要有两大类方法：目视解译和计算机解译。其中目视解译具有解译精度高的有点，但其效率低下且受制于人的主观因素的影响明显，因此已被计算机解译所取代，成为非主流方法。计算机解译目前主要有基于边缘检测、机器学习、超像素等算法。

遥感影像海岸线提取工作主要包括原始影像降噪处理和海岸线提取方式这两部分。

由于遥感影像拍摄受到传感器、天气状况等影像，所以要先对原始遥感影像进行降噪处理。目前降噪处理的主流方法主要有两种，一种是利用专业的影像处理软件对遥感影像进行诸如剪裁、拼接、大气校正、辐射校正等操作，如《高分影像砂质海岸线精细提取及校正方法》一文中按照辐射校正-几何校正-正射校正-影像融合的技术路线对原始影像进行预处理。这种方法可以依靠可靠的影像分析工具并充分利用遥感影像的多光谱信息，来实现较高的遥感影像处理效率，全面提升遥感影像的价值。另一种是基于滤波器的方法，通过设计不同的低通滤波器来去除图像噪声，如《遥感技术支持下海岸线提取算法》一文中使用了中值滤波和低通滤波两种滤波方法来实现遥感影像的降噪预处理，它们都是对像素点进行处理。以上所述降噪方法主要存在两个问题：一方面，为了得到更好的影像处理效果，这两种方法都需要人工手动去设置一些参数，而参数的设定受到人的经验判断的影响较大，无法保证最佳的参数设定，这就会导致影像细节丢失、降噪不充分等问题。另一方面，现有的主流降噪方法都没有很好地解决在处理遥感影像这类较大数据类型时的时间开销问题，比如影像的大气校正和辐射校正处理时间长，导致往往要实现花费大量的时间在降噪处理上，效率不高。

海岸线提取方式这部分，目前基于遥感影像进行海岸线提取的主流方法有阈值分割方法和边缘检测算子算法，阈值分割是一种简单、有效的图像分割方法，具有易于理解、实现简便、效果显著的特点。有许多经典的阈值分割算法可用于海岸线提取，但该方法阈值的选择具有一定的难度，而且往往由于海岸区域遥感影像中陆地与水体缺乏足够的灰度对比，会导致岸线提取效果不佳，且一般的阈值分割方法计算复杂度过高，难以适应一些高分辨率遥感影像海岸线的提取。所以与边缘检测算子算法相比，阈值分割算法总体上效果欠佳。边缘检测算子算法主要有Sobel算子、Roberts算子和Canny算子方法，如《基于边缘检测的海岸线自动提取研究》一文中，通过将传统Canny算子边缘检测和面向对象分类方法结合起来的手段，避免了阈值设置过于主观的问题，提高了提取精度，但仍存在提取的海岸线存在伪边缘现象、对影响噪声较为敏感、取出的岸线不连续的问题。

发明内容

首先，在原始影像降噪处理方面，本发明要解决现有两种降噪方法需要人工手动去设置一些参数，而参数的设定受到人的经验判断的影响较大，无法保证最佳的参数设定的问题，此问题会导致影像细节丢失、降噪不充分等问题。

进一步地，本发明要解决现有的主流降噪方法在处理遥感影像这类较大数据类型时的时间开销较大的问题，比如影像的大气校正和辐射校正处理时间长，导致往往要实现花费大量的时间在降噪处理上，效率不高。

更进一步地，在利用边缘检测方法提取海岸线方面，本发明要解决目前传统边缘检测方法提取的海岸线存在伪边缘现象、对影响噪声较为敏感、提取出的岸线有不连续的问题。

为解决以上问题，本发明提供一种基于聚类和边缘检测的遥感影像海岸线提取方法，创新性地提出一种改进型的K均值聚类分割算法用于遥感影像降噪处理，对遥感影像进行K均值聚类，降低影像复杂度，在较小的时间开销基础上完成对原始影像降噪处理的工作；然后利用聚类结果，通过构建Canny边缘检测算法，包括高斯滤波器的构造、像素点梯度的计算、非极大值抑制梯度幅值和双阈值选取检测，实现遥感影像中海岸线的提取工作。

技术方案：为解决上述技术问题，本发明采用的技术方案为：

本发明的技术方案分为海陆聚类分割和海岸线提取两部分，海陆分割部分包括原始数据获取、改进K均值算法实现遥感影像海陆聚类分割，海岸线提取部分包括Canny边缘检测算法的搭建，主要包括滤波器的设计构造、像素点梯度计算方向的选取和双阈值的选取。

在遥感影像海岸线提取工作中，为了克服传统影像降噪处理时所遇到的一些问题，本发明创新性地提出利用聚类分割手段先行处理原始数据，达到影像降噪的目的，选取使用K均值聚类算法实现海陆聚类分割，但面对传统K均值算法时间开销较大的问题，对算法进行相应改进，采用分批量训练子集的方式，将数据每一次划分限制在较小的子集中进行，在保证聚类效果的同时极大减少了算法的时间开销。

首先，获取高分辨率海岸地区遥感影像，作为原始数据进行处理。

其次，构建聚类算法，包括簇数K值的选取、训练子集的选取、质心的更新和算法终止的条件选择。簇数K值的选取是基于K均值聚类算法实现聚类分割的基础，K值的大小也就是要区分的地类的数量，K值的选择基于影像质量和研究区域地物情况，用于海岸线提取工作上我们只需区分出水体与陆地即可，一般取值2或3。

训练子集的选取关系到算法时间开销和最终效果，算法每一次迭代的处理都是在选出的训练子集中进行的，对于一个大的样本，首先随机抽取出部分对象构成一个训练子集，得到一个初步的K簇的训练模型，除了刚才已经选取过的对象，再选取部分对象构成新的训练子集，以此类推。这样一来，与传统K均值算法相比，在保证所有对象都能选取到的同时，避免了传统算法的遍历比较，减少了对象之间的重复计算，减少了算法时间开销。

质心的更新是算法结果不断优化的过程，对于每一次选出的训练子集，利用之前迭代得出的包含K个簇的训练模型(即K个质心的模型)，计算子集中每个对象到每个质心的距离，将每个对象划归到距离它最近的质心所属的簇类中，只要可以得到一个全新的包含K个簇的训练模型，再计算新模型中的每个簇的质心，得到更新后的质心。质心的计算遵循以下原则：每个簇的质心为簇中所有对象在各个维度的均值。

算法终止条件的选取也就是确定算法迭代次数，理论上迭代次数越多，效果越好，但这显然不现实，本发明的算法终止条件为：当所有对象都被选取到时或更新质心稳定时终止算法。终止时得到的最终的簇模型就作为最终聚类结果。

最后利用聚类结果，构建Canny边缘检测算法提取海岸线，包括滤波器的设计、像素点梯度计算、对梯度幅值进行非极大值抑制和设置双阈值检测。滤波器设计是Canny算法中对结果影响较为明显的一步，其目的是为了对数据进行进一步降噪处理，传统Canny算法使用的高斯滤波器，在对影像进行卷积的同时，产生的滤波结果会改变和丢失影像部分边缘，导致岸线提取结果不连续。滤波器的设计的目的，一是要保存像素点边界，二是对噪声进一步处理。像素点梯度计算的目的是为了比较找出影像中灰度强度变化最强的位置(即岸线位置)，传统梯度计算时只考虑水平与竖直方向上的梯度，这样会造成检测出许多假边缘，丢失真实边缘。像素点梯度计算方法是：一是对像素点进行一阶求导或二阶求导，二是计算方向不局限于水平竖直方向，应增加45度或135度方向梯度计算，这种模板求得的像素点梯度更加精确。对梯度幅值进行非极大值抑制的目的是为了清晰化基于梯度提取的边缘，主要步骤是：将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较，若当前像素梯度强度最大，则被保留为边缘像素，否则将被抑制。设置双阈值检测是为了对像素点阈值化，进行边缘连接。阈值化的方法是:如果边缘像素的梯度值高于高阈值，则将其标记为强边缘像素；如果边缘像素的梯度值小于高阈值并且大于低阈值，则将其标记为弱边缘像素；如果边缘像素的梯度值小于低阈值，则会被抑制。传统的双阈值的设置要人工根据输入影像的内容来决定，自适应较差，本发明将自适应算法引入阈值设定中来，可大大提高算法自适应性。

第一方面，提供一种基于聚类和边缘检测的遥感影像海岸线提取方法，包括：

S1、获取海岸地区遥感影像数据；

S2、利用预构建的改进K均值聚类算法对所述遥感影像数据进行海陆聚类分割处理得到海陆聚类结果；

S3、基于所述海陆聚类结果，利用预构建的Canny边缘检测算法提取海岸线。

在一些实施例中，步骤S2包括：对于遥感影像数据样本X，X＝{X

S21、第一次，随机抽取C

S22、第二次，在第一次的基础上，从X中随机抽取C

进一步地，在一些实施例中，簇数K值取值2或3。

进一步地，在一些实施例中，每个簇的质心为簇中所有对象在各个维度的均值。

进一步地，在一些实施例中，预设迭代停止条件为所有对象都被选取到或更新质心稳定。

在一些实施例中，步骤S3包括：

S31、采用中值滤波器对所述海陆聚类结果进行滤波降噪；

S32、对滤波降噪后的海陆聚类结果进行多角度计算像素点梯度，其中所述多角度计算包括水平、竖直、45°和135°四个方向；

S33、基于计算得到的像素点梯度，对梯度幅值进行非极大值抑制，得到边缘点；

S34、通过Otsu自适应算法确定设置的高、低阈值，根据边缘点的梯度大小和高、低阈值，判断每个边缘点的强弱性，按照设定原则对强边缘像素和弱边缘像素连接提取得到海岸线。

进一步地，在一些实施例中，步骤S31、采用中值滤波器对所述海陆聚类结果进行滤波降噪，包括：

设置平滑窗口大小为3×3或5×5，对窗口内所有像素值进行排序，用窗口内所有像素值的中值替换窗口中心像素的像素值；

按照平滑规则，依次将所有像素点的像素值替换为对应窗口的中值，完成滤波降噪，其中所述平滑规则为从左往右、从上往下依次进行。

进一步地，在一些实施例中，步骤S32，包括：

设定像素点梯度计算窗口大小，为3×3；

基于设定的像素点梯度计算窗口大小，使用Sobel算子计算由水平和竖直方向上计算得来的像素点梯度大小及方向；

基于设定的像素点梯度计算窗口大小，使用Rorberts算子计算由对角方向计算得来的像素点梯度大小及方向；

将由水平和竖直方向上计算得来的和由对角方向计算得来的像素点梯度大小及方向进行求平均，作为该像素点的梯度大小及方向。

进一步地，在一些实施例中，步骤S33，包括：基于计算得到的像素点梯度，对于每个像素点，在正、负梯度方向上分别选出一个距离该像素点最近的像素点，将该像素点与最近的两个像素点进行比较，若该像素点梯度最大，则保留作为边缘点，否则抑制该像素点；

进一步地，在一些实施例中，步骤S34，包括：

通过Otsu自适应算法确定设置的高、低阈值；根据边缘像素点的梯度大小G(x,y)和高、低阈值，判断每个边缘点的强弱性，判断原则为：

将所有强边缘像素连接起来，断开的地方用弱边缘像素补上，使边缘连续，得出所提取的海岸线。

第二方面，本发明提供了一种基于聚类和边缘检测的遥感影像海岸线提取装置，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据第一方面所述方法的步骤。

第三方面，本发明提供了一种存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现第一方面所述方法的步骤。

与现有技术相比，本发明的有益效果在于：

(1)本发明在利用遥感影像提取海岸线的工作中，创新性地用聚类分割手段代替传统的影像分析工具和一些低通滤波器，对遥感影像原始数据完成了降噪处理的工作。聚类分割在有效地区分出了影像中的水体与陆地，大大降低了影像复杂度，达到降噪目的同时，还避免了传统降噪方法需要人工设置参数导致影响细节丢失、降噪不充分的问题。

(2)本发明在利用聚类分割进行降噪处理时，改进了传统K均值聚类算法，针对传统K均值算法在更新簇时采取的遍历性的比较导致时间开销较大的问题，本发明采用分批量训练子集的方式，在保证聚类效果基础上，利用不同小的数据子集迭代更新质心的方式，大大减少了聚类的时间开销，平均时间缩短率达到了40％以上。聚类效果也没有减弱。

(3)本发明通过构建Canny边缘检测算法完成海岸线提取，设计了全新的中值滤波器来代替Canny算法传统意义上的高斯滤波器，克服了传统高斯滤波器会大量改变和丢失影像边缘，造成岸线提取结果不连续的问题。提取出的岸线连通率平均达到90％以上。

(4)本发明中像素梯度计算方向的增加，相较于传统只计算水平与竖直方向的梯度，这种模板求得的像素点梯度更加精确；本发明中自适应算法的引入，更加精确了双阈值的设定，省去了人工设置的烦恼，在提高了确定边缘点像素的准确率的同时，极大提高的算法的自适应度，适应了更多的应用场景。

附图说明

图1为本发明实施例提供的一种海岸线提取方法流程图。

图2为本发明实施例中对遥感影像进行海陆聚类分割的流程图。

图3为本发明实施例中构建Canny算法进行海岸线提取工作的流程图。

图4为本发明实施例中非极大值抑制流程图。

图5为本发明实施例中双阈值检测流程图。

具体实施方式

为使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施方式进一步阐述本发明。

在本发明的描述中，若干的含义是一个以上，多个的含义是两个以上，大于、小于、超过等理解为不包括本数，以上、以下、以内等理解为包括本数。如果有描述到第一、第二只是用于区分技术特征为目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量或者隐含指明所指示的技术特征的先后关系。

本文中术语“和/或”，仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况。另外，本文中字符"/"，一般表示前后关联对象是一种“或”的关系。

本发明的描述中，参考术语“一个实施例”、“一些实施例”、“示意性实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

实施例1

一种基于聚类和边缘检测的遥感影像海岸线提取方法，包括：

S1、获取海岸地区遥感影像数据；

S2、利用预构建的改进K均值聚类算法对所述遥感影像数据进行海陆聚类分割处理得到海陆聚类结果；

S3、基于所述海陆聚类结果，利用预构建的Canny边缘检测算法提取海岸线。

在一些具体实施例中，一种基于聚类和边缘检测的遥感影像海岸线提取方法，包括：

一、海陆聚类分割阶段：

获取高分辨率海岸地区遥感影像，作为原始数据；

构建聚类算法，包括簇数K值的选取、训练子集的选取、质心的更新和算法终止的条件选择，处理得到海陆聚类结果，作为岸线提取的输入。

二、海岸线提取阶段：

构建Canny边缘检测算法提取海岸线，包括滤波器的设计、像素点梯度计算、对梯度幅值进行非极大值抑制和设置双阈值检测；

滤波器的设计遵循原则：一是对于边界的保存情况要尽可能好，二是对于噪声的处理要有效。

像素点梯度计算：一是对像素点进行一阶求导或二阶求导，二是计算方向不局限于水平竖直方向，增加45度或135度方向梯度计算，这种模板求得的像素点梯度更加精确。

对梯度幅值进行非极大值抑制：基于所属像素点梯度计算的结果，将当前像素的梯度强度与沿正负梯度方向上的两个像素进行比较，若当前像素梯度强度最大，则被保留为边缘像素，否则将被抑制。

设置双阈值检测：采用自适应算法来确定双阈值，进而确定边缘点，连接边缘，得到最终岸线提取结果。

在本实施例中，提供一种基于聚类和边缘检测的遥感影像海岸线提取方法，具体实施平台为笔记本电脑，配置为：M1处理器(8核CPU，7核GPU)、16G内存，系统为MacOS系统(版本为13.0)。海岸地区遥感影像原始数据选择区域为盐城市射阳县至大丰县沿海地区，这一地区岸线连续，人为活动较少，适合用于海岸线提取研究。遥感影像数据来源于美国地质调查局(USGS)的Landsat-8卫星拍摄的遥感影像数据，Landsat-8卫星上有两种传感器，分别为OLI陆地成像仪和TIRS热红外传感器。此遥感数据共有11个波段，其中波段1-7和9的空间分辨率为30米，10-11波段为100米，波段8为15米分辨率的全色波段，卫星覆盖周期为16天。其中OLI陆地成像仪包含了一个专门用于海岸带观测的蓝色波段。此遥感影像数据具体参数如表1所示。

表1遥感影像原始数据具体参数表

首先对遥感影像原始数据进行海陆聚类分割，包括K值的选取、训练子集的选取、质心的更新和算法终止的条件选择。聚类分割的整体流程如图2所示。

构建聚类算法时K值的选取取决于想要区分的地类的数量，用于海岸线提取工作时只需区分出水体与陆地这两大类地类即可，K值的取值为2或3。

本实施例所使用的遥感地区的影像数据地类并不复杂，因此构建聚类算法时K值选择2，聚类区分出水体与陆地。

本发明改进的K均值算法，在本实施例中由于要区分的地物类别不多，训练子集的选取大小规定为50，选取方式为：对于一个大的样本X，X＝{X

本实施例在构建聚类算法是质心的更新，目的是要让簇中的相似度更高，遵循的目标函数为式(1)所示，

其中r

u

每次训练子集的迭代，划分对象到对应簇中的依据为对象到每一个簇质心的距离，这个距离有欧氏距离或曼哈顿距离。由于遥感影像数据为二维图像，因此不需要使用计算繁琐的欧氏距离，本实施例采用曼哈顿距离来进行衡量，可以简化计算。曼哈顿距离用来标明两个点在标准坐标系上的绝对轴距总和，坐标(X

S＝|X

其中，X

对象划分完成后就可以对簇模型中每个簇的质心进行更新，质心是簇内所有对象在各个维度的均值，计算公式如式(4)所示，

其中，C

算法的终止条件包括：

1.训练子集已完成所有对象的选取，这时为了避免对象之间的重复比较计算，应终止算法，输出结果；

2.簇模型中各个质心趋于稳定时，此时过多的迭代已失去意义。

本实施例将聚类算法的终止条件设定为训练子集完成对所有对象的选取后，训练子集的选取本就是为了避免对象之间的重复比较，因此当所有对象都被选取过后，停止迭代，最后得到的簇模型作为最终聚类结果输出。

然后，利用聚类结果，构建Canny边缘检测算法提取海岸线，包括滤波器的设计、像素点梯度计算、对梯度幅值进行非极大值抑制和设置双阈值检测，整体流程如图3所示。

本实例采用设计中值滤波器的方法来实现滤波器，具体步骤包括：规定像素点平滑规则，规定为从左往右、从上往下依次进行；设置平滑窗口大小，设置为3×3大小；对窗口内所有像素值进行排序，用这些像素值的中值替换窗口中心像素的像素值；按照平滑规则，依次将所有像素点的像素值替换为对应窗口的中值，完成滤波降噪。

本实施例中像素点梯度的计算，计算方向为水平、竖直、45°和135°四个方向，多角度计算梯度值会避免出现许多假边缘，计算结果会更加准确。具体计算方式包括：

首先设定像素点梯度计算窗口大小为3×3；其次使用Sobel算子计算由水平和竖直方向上计算梯度大小及方向；再用Rorberts算子计算由对角方向计算梯度大小及方向；最后两次结果平均，作为该像素点的梯度大小及方向。以水平和竖直方向为例，梯度计算方式为式(5)所示，

G

其中，G

由此，可以得出由水平和竖直方向计算出的该像素点梯度大小G和方向θ，计算方式如式(6)所示，

同理由对角方向的Rorberts算子S

根据上一步计算出的像素点梯度大小及方向，对所有像素点进行非极大值抑制，初步得出边缘点，其主要步骤包括：对应于每个像素点，选择出要和它进行比较的像素点；梯度值比较决定是否抑制，如图4所示。

本实施例中，对于每个像素点，非极大值抑制选择的比较对象为：在它的正负梯度方向上，与它距离最近的两个像素点。此三者进行比较，若此像素点梯度强度最大，则保留作为边缘点，否则抑制该像素点。

然后，设置双阈值检测，进行边缘点连接，得出提取的岸线，如图5所示。

本实施例通过Otsu自适应算法找出应设置的双阈值，比较得出每个边缘像素的强弱性，比较方法如式(7)所示，

其中，G(x,y)为该像素点的梯度值。

最后，将所有强边缘像素点连接起来，断开的地方尽量用弱边缘像素点补上，使边缘尽可能连续，得出所提取的海岸线。

实施例2

第二方面，本实施例提供了一种基于聚类和边缘检测的遥感影像海岸线提取装置，包括处理器及存储介质；

所述存储介质用于存储指令；

所述处理器用于根据所述指令进行操作以执行根据实施例1所述方法的步骤。

实施例3

第三方面，本实施例提供了一种存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现实施例1所述方法的步骤。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

由技术常识可知，本发明可以通过其它的不脱离其精神实质或必要特征的实施方案来实现。因此，上述公开的实施方案，就各方面而言，都只是举例说明，并不是仅有的。所有在本发明范围内或在等同于本发明的范围内的改变均被本发明包含。