LXX
读完需要
速读仅需 3 分钟
前言：
DBSCAN(Density-Based Spatial Clustering of Applications with Noise)是一个比较有代表性的基于密度的聚类算法。与划分和层次聚类方法不同，它将簇定义为密度相连的点的最大集合，能够把具有足够高密度的区域划分为簇，并可在噪声的空间数据库中发现任意形状的聚类。
但是在应用于 ICESat-2 光子点去噪过程中时，由于不同地形条件下光子数据点的分布情况，密度大小不同，在进行参数选择时会有一定的复杂度；同时传统的 DBSCAN 算法中的圆形搜索区域在地形复杂的条件下存在一定的局限性。因此有学者将其圆形搜索区域改为椭圆搜索区域提升去噪性能，本文分享这一改进的复现过程。
希望各位同学点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
1
DBSCAN 算法的核心思想是以数据中的某个数据点为计算中心，计算该数据点邻域范围内的其他数据点数量，以此来判断该数据点是否为核心点、边界点或噪声点。其中，当该数据点的邻域范围（ε）内至少存在 MinPts 个数据点时定义为核心点；当某个非核心数据点被包括在任意核心点的邻域范围内时被定义为边界点；当某个非核心数据点不被包括在任意核心点的邻域范围内时被定义为噪声点。具体来说，该算法首先需要设置两个参数：邻域半径(ε)和邻域最小点数(MinPts)，其中邻域半径定义了一个数据点的圆形搜索区域的邻域大小，邻域最小点数定义了一个核心点必须拥有的邻近点数量。若数据集中一个数据点的邻域内有至少 MinPts 个数据点，则将其标记为核心点，否则标记为噪声点。接下来，算法将从核心点开始向其邻近的每个点进行遍历搜索，将其邻域内的数据点标记为边界点，并将这些点加入到该核心点所在的簇中。如果某个边界点的邻域内有另一个核心点，则将这两个簇合并。最后重复这个过程，直到所有的数据点都被访问过为止。下图显示了 DBSCAN 算法中点类别的确定方法。
其原理为系统在众多样本点中随机选中一个，围绕这个被选中的样本点画一个圆，规定这个圆的半径以及圆内最少包含的样本点，如果在指定半径内有足够多的样本点在内，那么这个圆圈的圆心就转移到这个内部样本点，继续去圈附近其它的样本点。其中 A 为核心点，B、C 为边界点，D 为离群点。根据其原理，需要规定两个参数，即 epsilon：半径；minPoints：包含（圈住）样本点个数。半径的设定一般通过距离来判别，也就是要找到突变点；样本点的个数一般都是偏小一些，然后进行多次尝试。
算法具体过程为：
１）选定 DBSCAN 的两个参数，邻域半径（eps）和核心对象个数（minPoints）
２）以同一个研究区中每个光子作为聚类对象，以在给定半径的 eps 内包含的光子数量不能少于某一阈值 minPoints 作为聚类基础，在此基础上对具有密度连接特性的光子对象进行聚类
３）遍历研究区数据集中的光子云数据，只要其邻近区域内光子数据的密度超过 minPoints 则能够继续聚类，这些光子对象即为信号光子，否则，标识这些光子为噪声光子。
可视化算法过程如下：
https://www.naftaliharris.com/blog/visualizing-dbscan-clustering
【详细实现代码见前期】
DBSCAN
LXXtime，公众号：遥感小屋基于DBSCAN算法的ICESat-2光子点云去噪方法
2
传统的DBSCAN算法邻域定义为给定光子间的欧几里得距离，搜索域定义为圆形。通过观察分析噪声光子与信号光子分布可知，太阳光的噪声光子分布于整个研究区的上空与地面下，而且沿轨方向的光子云密度与垂轨方向的光子云密度差异很大，传统的圆形搜索域并不适用于星载光子云数据的去噪处理。
在森林研究区的光子不仅会集中在地面附近，而且会分布于冠层区域，在森林植被覆盖度适中情况下，地面光子密度会高于冠层光子密度，覆盖度较大情况下，地面光子密度会低于冠层光子密度。由于森林植被类型不同，其反射表面轮廓比人工制造建物的结构要复杂的多，这为光子云精准去噪提升了难度。因此根据光子分布形式，后续学者采用基于椭圆搜索域的DBSCAN 改进算法，算法将光子云数据搜索域改为水平方式的椭圆形状。
2.1
搜索区域
改进的DBSCAN算法核心为通过修改圆形搜索区域为椭圆形，从而更加适用于ICESat-2光子去噪工作。公式如下：
式中，x为光子沿轨距离，单位为 m，h为光子点云高程，单位为m，a为椭圆搜索域的长半轴，b为椭圆搜索域的短半轴，根据前人研究，在设置参数时，采用的长短半轴之比为6:1。
2.2
参数设置
椭圆搜索区域的参数设置如下：
其中along_track_eps为椭圆搜索域的长半轴，height_eps为椭圆搜索域的短半轴；而min_pts则为领域最小点数。
参数设置上，一般长短半轴之比为6:1，而领域最小点数min_pts需要多次实验确定。
# 设置参数along_track_eps = 18 # 6的倍数height_eps = 6 # 可调整min_pts = 10 # 可调整
3
1.改进DBSCAN算法：
"""# -*- coding: utf-8 -*-@Time: 2024/1/6 10:13@Author: LXX@File: improveDBSCAN.py@IDE：PyCharm@Motto：ABC(Always Be Coding)"""class ImprovedDBSCAN():#定义的椭圆搜索区域def _eps_neighborhood(self, p, q, eps):in_ellipse = ((q[0] - p[0])**2 / eps[0]**2 + (q[1] - p[1])**2 / eps[1]**2) <= 1return in_ellipsedef _region_query(self, m, point_id):n_points = m.shape[1]seeds = []......return seeds#进行聚类，与传统DBSCAN相同def _expand_cluster(self, m, classifications, point_id, cluster_id):seeds = self._region_query(m, point_id)if len(seeds) < self.min_pts:classifications[point_id] = cluster_idfor seed_id in seeds:classifications[seed_id] = cluster_idwhile len(seeds) > 0:results = self._region_query(m, current_point)if len(results) >= min_pts:result_point = results[i]if classifications[result_point] == classifications[result_point] == NOISE:if classifications[result_point] == UNCLASSIFIED:seeds.append(result_point)classifications[result_point] = cluster_idseeds = seeds[1:]return True......
2.传统DBSCAN算法：
# 数据标准化scaler = StandardScaler()scaled_data = scaler.fit_transform(required_data)# 使用DBSCAN算法进行聚类start_time = time.time()dbscan = DBSCAN(eps=0.3, min_samples=5)dbscan.fit(scaled_data)end_time = time.time()print("DBSCAN Execution Time: {:.2f} seconds".format(end_time - start_time))# 获取聚类结果labels = dbscan.labels_
3.加载数据：
# 读取CSV数据文件root = tk.Tk()root.withdraw()file_path = filedialog.askopenfilename(filetypes=[('CSV Files', '*.csv')])# 读取CSV文件data = pd.read_csv(file_path, header=0)# 获取原文件路径、文件名和后缀file_dir, file_name = os.path.split(file_path)file_prefix, file_suffix = os.path.splitext(file_name)required_data = data[['Along-Track (m)', 'Height (m HAE)']]
4.设置参数：
# 设置参数along_track_eps = 18 # 6的倍数height_eps = 6 # 可调整min_pts = 10 # 可调整
5.运行改进DBSCAN和传统DBSCAN并显示结果：
# 绘制图形plt.figure(figsize=(12, 6))# 绘制ImprovedDBSCAN结果图plt.subplot(1, 2, 1)plt.scatter(required_data['Along-Track (m)'], required_data['Height (m HAE)'], c=classifications_binary, cmap='RdYlGn')plt.rcParams['font.sans-serif'] = ['STSong']plt.rcParams['axes.unicode_minus'] = Falseplt.xlabel('Along-Track/m')plt.ylabel('Height/m')plt.title('改进DBSCAN')# 绘制DBSCAN结果图plt.subplot(1, 2, 2)plt.scatter(required_data['Along-Track (m)'], required_data['Height (m HAE)'], c=colors)plt.rcParams['font.sans-serif'] = ['STSong']plt.rcParams['axes.unicode_minus'] = Falseplt.xlabel('Along-Track/m')plt.ylabel('Height/m')plt.title('传统DBSCAN')plt.tight_layout()plt.show()
椭圆形搜索区域示意图
4
结果如下，可以看到相比于传统DBSCAN算法具有一定的提高。但仍存在以下几点问题：
1.领域最小点数min_pts需要多次实验确定，何时是最佳阈值仍无法确定；
2.搜索时仍是对所有点进行遍历，在地形复杂区域仍会产生一定的误差；
3.未考虑坡度因素；
4.椭圆应根据搜索方向旋转角度，水平搜索具有局限性。
后续改进......
假期结束，大家工作顺利！
