CN113780085A - 一种近海单光子去噪分类方法 - Google Patents

Abstract

本发明涉及一种近海单光子去噪分类方法，包括初步提取海面点云光子信息，分段并拟合海平面，计算分段高程位置的点密度并建立密度曲线，结合密度曲线初步判定出海平面光子，进而划分出海平面上下的区域，并分别通过计算椭圆距离以及邻近光子距离和的方式完成岛礁光子和水下光子的去噪及分类。与现有技术相比，本发明通过先确定海平面，初步分出上下两部分进而进行去噪的方式，实现了对近海单光子的多级去噪以及分类，具有分类结果更精确，光子去噪效果好等优点。

Description

一种近海单光子去噪分类方法

技术领域

本发明涉及海洋测绘领域，尤其是涉及一种近海单光子去噪分类方法。

背景技术

海洋测绘是陆地测绘的延伸，海洋的位置测定与水深测定是构建海洋地理信息系统的关键，也是海洋最基本的信息元之一。随着科技的不断发展，声波测深仪器以及激光多波束系统开始广泛应用于海洋测深，其主要搭载在测深船上，水深测量精度有了大幅度提高，但受限于近浅海的水深，测量速度受到了很大的限制，测量船不能广泛活动在近海、浅海水域。

ICESat-2搭载的ATLAS系统是世界第一台多波束单光子激光测高系统，星载激光测高技术具有探测方向性好，测高精度高等特点，能够广泛运用在地形探测、生物量探测、极地冰盖高程变化检测等多个研究方向。在水深探测方面，根据其已公布的数据，ICESat-2数据能测得水下近30m数据，这为测深船以及无人机不能到达的岛礁、暗礁区域的浅海测量提供了解决方案。但是除了从地面返回的光子初始海平面光子之外，大量的噪声光子(太阳光背景噪声、仪器背景噪声等)也会被探测器探测到，这就使得原始的光子点云存在大量的噪声，在检测时误差较大，无法精确提取并分类。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种近海单光子去噪分类方法，解决了现有技术提取光子时噪声较大，无法精确分类的问题。

本发明的目的可以通过以下技术方案来实现：

一种近海单光子去噪分类方法，包括：

S1、获取点云光子，对所述点云光子的沿轨距离坐标和高程坐标进行尺度变换，使用OPTICS聚类方法得到噪声阈值，保留可达距离小于噪声阈值的点云光子为初始海平面光子；

S2、在沿轨方向上以一号长度值对所有初始海平面光子分段，将每一段的中点坐标与高程坐标作为拟合点；

S3、计算每一分段的高程坐标的点密度，以二号长度值在沿轨方向上对初始海平面光子重新进行分段，以连续的三段为一个探测单元，根据采样间隔阈值绘制每个探测单元的密度曲线；

S4、根据得到的密度曲线，对每个探测单元的中间段作如下判断：

若密度曲线在中间段中仅存在一个峰值，且该峰值与拟合海平面高度差值在第一高度阈值之内，则将该峰值所在处设定为海平面信号峰；

若密度曲线在中间段中存在两个峰值，选取与拟合海平面高度差值最小且在第一高度阈值之内的峰值所在处为海平面信号峰，将另一个峰值所在处判定为水下信号峰；

若密度曲线在中间段中不存在峰值或者不存在与海平面高度差值在第一高度阈值之内的峰值，则处理下一个探测单元；

S5、以海平面信号峰高程为原始高程，向上检测每个初始海平面光子的密度，当密度低于密度阈值时，停止检测，将该初始海平面光子所在的高程位置定义为上边界；

以海平面信号峰高程为原始高程，向下检测每个初始海平面光子的密度，当密度低于密度阈值时，停止检测，将该初始海平面光子所在的高程位置定义为下边界，同时，当所述下边界的高度低于海平面信号峰与水下信号峰之间的初始海平面光子的密度最低点的高程时，将所述密度最低点的高程位置定义为下边界；

将位于所述上边界高度以上的光子设为初始岛礁光子，位于下边界高度以下的光子设为初始水下光子；

S6、以每一个初始岛礁光子为中心绘制椭圆，计算椭圆内其它初始岛礁光子到中心的初始岛礁光子的椭圆距离，对椭圆距离进行加权计算，得到所有初始岛礁光子的光子密度特征值；

利用最大类间方差法对所述初始岛礁光子的光子密度特征值求取第一特征阈值并初步分类，光子密度特征值小于第一特征阈值的初始岛礁光子为第一噪声光子，光子密度特征值大于第一特征阈值的初始岛礁光子为岛礁光子；

计算出所有第一噪声光子的光子密度特征值均值和方差，求取第二特征阈值，将光子密度特征值大于第二特征阈值的第一噪声光子重新定义为岛礁光子，光子密度特征值小于第二阈值的第一噪声光子不进行处理；

S7、以多个固定长宽的方格对所有初始水下光子进行分割，并统计每一个方格内初始水下光子的个数，其中定义每个方格周围的八个方格为邻域范围，此时将该方格定义为中心方格，将所述邻域范围内光子数目最小的多个非零方格定义为噪声方格；

计算每一个噪声方格中每一个初始水下光子的若干个最邻近光子的距离和作为噪声特征值，并计算噪声方格中所有初始水下光子的噪声特征值的均值作为每一个噪声方格的局部噪声特征值，根据所述局部噪声特征值、信噪比，计算得到每一个噪声方格对应的中心方格的局部信号阈值，将在中心方格中噪声特征值大于局部信号阈值的初始水下光子设为水下光子，噪声特征值小于局部信号阈值的初始水下光子舍去；

对所有水下光子在探测区域内以一号距离阈值进行局部拟合，删除偏离拟合直线超过二号距离阈值的水下光子，直至遍历完探测区域。

进一步地，步骤S2中的拟合点可进行修正，具体步骤包括：

根据拟合点得到拟合海平面和拟合直线；

计算所有拟合点的标准差，以及计算拟合点和拟合直线的拟合残差，删除所有拟合点中拟合残差大于多倍标准差的拟合点；

根据高程修正公式对拟合海平面进行修正，删除所有拟合点中高于修正后拟合海平面第二高度阈值的拟合点，得到修正后拟合点；

对修正后拟合点重复以上步骤直至保留的所有拟合点的拟合残差均小于设定倍数的标准差，即确定为最终拟合点。

进一步地，步骤S6中所述椭圆距离dist(p,q,θ)计算表达式为：

式中，p表示中心的初始岛礁光子，q表示其它任意的一个初始岛礁光子，θ表示坡度，ΔX_θ代表在坡度为θ时p、q的水平距离，Δh_θ代表在坡度为θ时p、q的竖直距离，X_p、X_q分别代表p、q的沿轨坐标，h_p、h_q分别代表p、q的高程坐标，a与b分别表示椭圆的长轴长度和短轴长度。

进一步地，所述坡度通过变步长搜索的方法计算得到，具体包括：

在一号角度范围内，以第一角度阈值为定步长，初步确定坡度θ₀；

根据θ₀的值，在二号角度范围内以第二角度阈值为定步长，最终确定坡度。

进一步地，步骤S7中所述局部信号阈值threshold_local的计算表达式如下：

式中，m表示信噪比，

表示局部噪声特征值。

进一步地，步骤S7中所述信噪比的具体设定规则如下：

水深数值小于a1，信噪比设置m＞b1，

水深数值在a1～a2内，信噪比设置b2≤m≤b1，

水深数值大于a2，信噪比设置m＜b2。

进一步地，步骤S3中所述密度ρ_光子的计算表达式为：

ρ_光子＝n_total/(x_width·h_range)

式中，x_width表示分段包含的沿轨距离，h_range表示初始海平面光子分布区间高度，n_total表示在范围内初始海平面光子的总数。

进一步地，步骤S3中所述密度曲线的采样间隔小于h_range/2。

进一步地，步骤S5中所述密度阈值ρ_threshold计算表达式为：

ρ_threshold＝ρ_sea-0.5σ_sea

式中，ρ_sea表示密度的均值，σ_sea表示密度的标准差。

进一步地，步骤S6中所述第二特征阈值threshold计算表达式为：

threshold＝μ_noise+3·σ_noise

式中，μ_noise表示第一噪声光子的密度特征值均值，σ_noise代表第一噪声光子的密度特征值方差。

与现有技术相比，本发明具有以下优点：

1、本发明通过提取出点云光子信息，拟合海平面后求取密度曲线，根据密度曲线提取出海平面信号峰，利用海平面高度初步划分上下区域，并在岛礁光子和水下光子的去噪中分别计算了椭圆距离和最邻近光子距离和，通过不同区域光子的特点采用不同的计算方法求取光子信息，实现光子的分类，保证了数据的精确性，同时也对数据进行了多级修正，使得提取出的光子噪声较小；

2、在计算过程中所使用的参数如水下信噪比或坡度等都使用了特定算法进一步限定了这些参数的范围，在多级去噪的算法流程的基础上使用了更精确的参数，更能反映实际的情况，使得提取的单光子信息更为准确。

附图说明

图1为本发明的流程示意图。

图2为本发明处理初始海平面光子时的探测单元分段示意图。

图3为本发明处理初始水下光子时的方格划分示意图。

图4为本发明处理海面与水下信号分离情况时的效果示意图。

图5为本发明处理岛礁、海面与水下信号混合情况时的效果示意图。

图6为本发明处理海面与水下信号相切情况时的效果示意图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

本实施例提供了一种近海单光子多级去噪分类方法，流程如图1所示，具体包括：

使用ICEsat-2卫星获取光子点云数据，并对点云光子的沿轨距离坐标和高程坐标进行尺度变换，变换公式如下：

式中，x代表点云的沿轨距离坐标；y代表点云的高程坐标；a代表高程方向的转换因子，对于提取海平面信号，a的值为[80,120]；x′为转换后的沿轨距离坐标；y′为转换后的高程坐标。

使用OPTICS聚类方法得到噪声阈值，先计算每个点的核心距离(CD)和可达距离(RD)作为每个点的光子特征值，计算表达式如下：

式中，ε表示每一个点的搜索半径，因OPTICS算法对ε不敏感，通常该值可取为∞；MinPts表示在搜索范围内的最小点数；ε′表示从点o到达其最邻近第k个点的距离，N_ε(o)为在点o规定范围ε内的光子数；RD_ε,MinPts(q,o)代表了点o，q之间的可达距离。再使用最大类间方差法(OSTU)选取噪声阈值t，计算表达式如下：

式中，信号光子数与光子总数的比值记为ω₀，初始海平面光子的平均特征值为μ₀；噪声光子数与光子总数的比值为ω₁，噪声光子的平均特征值为μ₁。当

最大时，与之对应的t即为阈值。

保留可达距离(RD)小于噪声阈值的点云光子为初始海平面光子，其余为噪声光子。

在沿轨方向上以20米为一段对所有初始海平面光子分段，将每一段的中点坐标与高程坐标作为拟合点，对海平面进行拟合，计算表达式如下：

Y＝[y₁,y₂,…,y_k,…y_n]^T

式中，x_k代表第k个分段的沿轨距离的中心坐标，y_k代表第k个分段的统计到的初始海平面光子最大值，对X和Y构建线性方程，则有：

Y＝Xb

式中，b代表拟合系数，对海平面的拟合即是对b的求解，则b的求解公式如下：

b＝(X^TX)^-1X^TY

对于得到的拟合结果，有时海平面情况较复杂，需要对拟合结果进行修正，具体修正步骤如下：

根据拟合点得到拟合海平面和拟合直线；

计算所有拟合点的标准差，以及计算拟合点和拟合直线的拟合残差，计算表达式如下：

Δy_k＝y_k-y′_k

式中，y_k代表第k段的最大高程，y′_k代表第k段的拟合高程，Δy_k代表每一段的拟合残差。

删除所有拟合点中拟合残差大于3倍标准差的拟合点；

根据高程修正公式对拟合海平面进行修正，修正表达式如下：

y′＝y-a/2

式中，h代表了拟合的海平面高程，a代表了海面光子信号的平均展宽，k′代表修正后的拟合海平面。根据经验，a取0.2，可根据实际情况进行调节。

删除所有拟合点中高于修正后拟合海平面1米的拟合点，得到修正后拟合点；

对修正后拟合点重复以上步骤直至保留的所有拟合点的拟合残差均小于3倍的标准差，即确定为最终拟合点。

计算每一分段的高程坐标的点密度，计算表达式为：

ρ_光子＝n_total/(x_width·h_range)

式中，x_width代表分段包含的沿轨距离，h_range代表计算临近区间的范围，n_total代表在范围内光子的总数。h_range的选取与光子信号的分布有密切的关系，对于海洋表面信号，由于其展宽小，h_range不超过0.2m。

以40～60m的长度在沿轨方向上对初始海平面光子重新进行分段，以连续的三段为一个探测单元，如图2所示，根据小于h_range/2的采样间隔绘制每个探测单元的密度曲线。根据得到的密度曲线，对每个探测单元的中间段作如下判断：

若密度曲线在中间段中仅存在一个峰值，且该峰值与拟合海平面高度差值在1m之内，则将该峰值所在处设定为海平面信号峰；

若密度曲线在中间段中存在两个峰值，选取与拟合海平面高度差值最小且在第1m之内的峰值所在处为海平面信号峰，将另一个峰值所在处设为水下信号峰；

若密度曲线在中间段中不存在峰值、超过两个峰值或者不存在与海平面高度差值在1m之内的峰值，则处理下一个探测单元；

从信号峰的峰值开始，向上和向下进行搜寻，当光子密度低于ρ_threshold(表达式如下)：

ρ_threshold＝ρ_sea-0.5σ_sea

式中，ρ_sea表示密度的均值，σ_sea表示密度的标准差。

搜寻结束，将上下最后搜寻到的光子高程位置分别设定为海平面上边界和下边界，同时，当下边界的高度低于海平面信号峰与水下信号峰之间的初始海平面光子的密度最低点的高程时，将密度最低点的高程位置定义为下边界。

将位于上边界高度以上的光子设为初始岛礁光子，位于下边界高度以下的光子设为初始水下光子。

以每一个初始岛礁光子为中心绘制椭圆，计算椭圆内其它初始岛礁光子到中心的初始岛礁光子的椭圆距离，

式中，p表示中心的初始岛礁光子，q表示其它任意的一个初始岛礁光子，θ表示坡度，ΔX_θ代表在坡度为θ时p、q的水平距离，Δh_θ代表在坡度为θ时p、q的竖直距离，X_p、X_q分别代表p、q的沿轨坐标，h_p、h_q分别代表p、q的高程坐标，a与b分别表示椭圆的长轴长度和短轴长度。

对椭圆距离进行加权计算，得到所有初始岛礁光子的光子密度特征值；

利用最大类间方差法对初始岛礁光子密度特征值求取第一阈值并初步分类，光子密度特征值小于第一特征阈值的初始岛礁光子为第一噪声光子，光子密度特征值大于第一特征阈值的初始岛礁光子为岛礁光子；

计算出所有第一噪声光子的光子密度特征值均值和方差，求取第二特征阈值，计算表达式为：

threshold＝μ_noise+3·σ_noise

式中，μ_noise表示第一噪声光子的密度特征值均值，σ_noise代表第一噪声光子的密度特征值方差。将光子密度特征值大于第二特征阈值的第一噪声光子重新定义为岛礁光子，光子密度特征值小于第二阈值的第一噪声光子不进行处理。

以多个固定长宽的方格对所有初始水下光子进行分割，并统计每一个方格内初始水下光子的个数，其中定义每个方格周围的八个方格为邻域范围，此时将该方格定义为中心方格，将邻域范围内光子数目最小的3个非零方格定义为噪声方格,如图3所示。

计算每一个噪声方格中每一个初始水下光子的若干个最邻近光子的距离和作为噪声特征值，计算表达式如下：

式中，x′_p、x′_q分别代表p、q两点转换后的沿轨距离坐标，y′_p、y′_q分别代表p、q两点转换后的高程坐标，并计算噪声方格中所有初始水下光子的噪声特征值的均值作为每一个噪声方格的局部噪声特征值，计算表达式如下：

式中，D₁代表了中心方格区域，D₂代表噪声方格区域，n代表落在噪声方格区域内光子的个数。根据局部噪声特征值、信噪比，计算得到每一个噪声方格对应的中心方格的局部信号阈值，计算表达式如下：

式中，m表示该水深范围内的设定的信噪比。通常情况下，水深范围将会被划定为三个范围，分别为水深I区、水深II区、和水深III区，水深I区的水深范围通常<5m，信噪比设置通常m＞10，水深II区的水深大致在5～15m内，信噪比通常设置为3≤m≤10，III区指>15m水深的区域，信噪比设置m＜3。根据实际情况，水深区域的划分与信噪比的选择可以进行调整。

将在中心方格中噪声特征值大于局部信号阈值的初始水下光子设为水下光子，噪声特征值小于局部信号阈值的初始水下光子舍去；

对所有水下光子在探测区域内以沿轨距离为Δdist的范围内进行局部拟合，删除偏离拟合直线3m的水下光子，将拟合区域偏移Δdist/2，直至遍历完探测区域。

图4～图6展示了不同情况下本实施例对于光子信号的分类效果，分别为海面与水下信号分离情况、岛礁海面与水下信号混合情况以及海面与水下信号相切的情况，由图可知，本实施例在不同情况下均能实现对不同近海单光子的区分，且分类效果良好。

以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术人员无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。

Claims (10)

1.一种近海单光子去噪分类方法，其特征在于，包括：

S1、获取点云光子，对所述点云光子的沿轨距离坐标和高程坐标进行尺度变换，使用OPTICS聚类方法得到噪声阈值，保留可达距离小于噪声阈值的点云光子为初始海平面光子；

S2、在沿轨方向上以一号长度值对所有初始海平面光子分段，将每一段的中点坐标与高程坐标作为拟合点；

S3、计算每一分段的高程坐标的点密度，以二号长度值在沿轨方向上对初始海平面光子重新进行分段，以连续的三段为一个探测单元，根据采样间隔阈值绘制每个探测单元的密度曲线；

S4、根据得到的密度曲线，对每个探测单元的中间段作如下判断：

若密度曲线在中间段中仅存在一个峰值，且该峰值与拟合海平面高度差值在第一高度阈值之内，则将该峰值所在处设定为海平面信号峰；

若密度曲线在中间段中存在两个峰值，选取与拟合海平面高度差值最小且在第一高度阈值之内的峰值所在处为海平面信号峰，将另一个峰值所在处判定为水下信号峰；

若密度曲线在中间段中不存在峰值或者不存在与海平面高度差值在第一高度阈值之内的峰值，则处理下一个探测单元；

S5、以海平面信号峰高程为原始高程，向上检测每个初始海平面光子的密度，当密度低于密度阈值时，停止检测，将该初始海平面光子所在的高程位置定义为上边界；

以海平面信号峰高程为原始高程，向下检测每个初始海平面光子的密度，当密度低于密度阈值时，停止检测，将该初始海平面光子所在的高程位置定义为下边界，同时，当所述下边界的高度低于海平面信号峰与水下信号峰之间的初始海平面光子的密度最低点的高程时，将所述密度最低点的高程位置定义为下边界；

将位于所述上边界高度以上的光子设为初始岛礁光子，位于下边界高度以下的光子设为初始水下光子；

S6、以每一个初始岛礁光子为中心绘制椭圆，计算椭圆内其它初始岛礁光子到中心的初始岛礁光子的椭圆距离，对椭圆距离进行加权计算，得到所有初始岛礁光子的光子密度特征值；

利用最大类间方差法对所述初始岛礁光子的光子密度特征值求取第一特征阈值并初步分类，光子密度特征值小于第一特征阈值的初始岛礁光子为第一噪声光子，光子密度特征值大于第一特征阈值的初始岛礁光子为岛礁光子；

计算出所有第一噪声光子的光子密度特征值均值和方差，求取第二特征阈值，将光子密度特征值大于第二特征阈值的第一噪声光子重新定义为岛礁光子，光子密度特征值小于第二阈值的第一噪声光子不进行处理；

S7、以多个固定长宽的方格对所有初始水下光子进行分割，并统计每一个方格内初始水下光子的个数，其中定义每个方格周围的八个方格为邻域范围，此时将该方格定义为中心方格，将所述邻域范围内光子数目最小的多个非零方格定义为噪声方格；

计算每一个噪声方格中每一个初始水下光子的若干个最邻近光子的距离和作为噪声特征值，并计算噪声方格中所有初始水下光子的噪声特征值的均值作为每一个噪声方格的局部噪声特征值，根据所述局部噪声特征值、信噪比，计算得到每一个噪声方格对应的中心方格的局部信号阈值，将在中心方格中噪声特征值大于局部信号阈值的初始水下光子设为水下光子，噪声特征值小于局部信号阈值的初始水下光子舍去；

对所有水下光子在探测区域内以一号距离阈值进行局部拟合，删除偏离拟合直线超过二号距离阈值的水下光子，直至遍历完探测区域。

2.根据权利要求1所述的一种近海单光子去噪分类方法，其特征在于，步骤S2中的拟合点可进行修正，具体步骤包括：

根据拟合点得到拟合海平面和拟合直线；

计算所有拟合点的标准差，以及计算拟合点和拟合直线的拟合残差，删除所有拟合点中拟合残差大于多倍标准差的拟合点；

根据高程修正公式对拟合海平面进行修正，删除所有拟合点中高于修正后拟合海平面第二高度阈值的拟合点，得到修正后拟合点；

对修正后拟合点重复以上步骤直至保留的所有拟合点的拟合残差均小于设定倍数的标准差，即确定为最终拟合点。

3.根据权利要求1所述的一种近海单光子去噪分类方法，其特征在于，步骤S6中所述椭圆距离dist(p，q，θ)计算表达式为：

式中，p表示中心的初始岛礁光子，q表示其它任意的一个初始岛礁光子，θ表示坡度，ΔX_θ代表在坡度为θ时p、q的水平距离，Δh_θ代表在坡度为θ时p、q的竖直距离，X_p、X_q分别代表p、q的沿轨坐标，h_p、h_q分别代表p、q的高程坐标，a与b分别表示椭圆的长轴长度和短轴长度。

4.根据权利要求3所述的一种近海单光子去噪分类方法，其特征在于，所述坡度通过变步长搜索的方法计算得到，具体包括：

在一号角度范围内，以第一角度阈值为定步长，初步确定坡度θ₀；

根据θ₀的值，在二号角度范围内以第二角度阈值为定步长，最终确定坡度。

5.根据权利要求1所述的一种近海单光子去噪分类方法，其特征在于，步骤S7中所述局部信号阈值threshold_local的计算表达式如下：

式中，m表示信噪比，

表示局部噪声特征值。

6.根据权利要求1所述的一种近海单光子去噪分类方法，其特征在于，步骤S7中所述信噪比的具体设定规则如下：

水深数值小于a1，信噪比设置m＞b1，

水深数值在a1～a2内，信噪比设置b2≤m≤b1，

水深数值大于a2，信噪比设置m＜b2。

7.根据权利要求1所述的一种近海单光子去噪分类方法，其特征在于，步骤S3中所述密度ρ_光子的计算表达式为：

ρ_光子＝n_total/(x_width·h_range)

式中，x_width表示分段包含的沿轨距离，h_range表示初始海平面光子分布区间高度，n_total表示在范围内初始海平面光子的总数。

8.根据权利要求7所述的一种近海单光子去噪分类方法，其特征在于，步骤S3中所述密度曲线的采样间隔小于h_range/2。

9.根据权利要求1所述的一种近海单光子去噪分类方法，其特征在于，步骤S5中所述密度阈值ρ_threshold计算表达式为：

ρ_threshold＝ρ_sea-0.5σ_sea

式中，ρ_sea表示密度的均值，σ_sea表示密度的标准差。

10.根据权利要求1所述的一种近海单光子去噪分类方法，其特征在于，步骤S6中所述第二特征阈值threshold计算表达式为：

threshold＝μ_noise+3·σ_noise

式中，μ_noise表示第一噪声光子的密度特征值均值，σ_noise代表第一噪声光子的密度特征值方差。

Images