基于ICESat-2/ATLAS数据的近海岸水深提取

摘要

水深是表征海洋浅水和海岸环境的重要地形要素，光子计数激光雷达可穿透一定深度的水体，为水深信息提取提供可靠的数据支持。以我国南海岛礁为例，利用目前唯一在轨的星载光子计数激光雷达-ICESat-2/ATLAS数据开展岛礁浅水水深提取及精度评价研究。首先根据置信度参数对原始光子数据进行粗去噪，基于点密度分布差异分离水面和水底光子；然后对水面光子采用区间估计方法精去噪，利用RANSAC算法拟合水面高程；通过改进滤波参数，基于改进OPTICS算法对水底光子进行两次聚类，实现水底光子的精去噪，进而通过折射校正和潮汐校正提取近岸水深；最后利用机载LiDAR测深数据进行验证。实验结果表明，与ATL03高置信度光子数据和AVEBM去噪结果相比，该精去噪算法具有更高的F值，分别提高了约5.87％和3.38％；水深提取结果与机载测深数据的R²为0.91，均方根误差RMSE为0.53m。

关键词

ICESat-2/ATLAS; 光子去噪; 激光测深; 区间估计; 改进OPTICS; 近岸水深

近海岸水深测量是海洋测绘的基础工作，传统的借助于船载平台的单波束或多波束测深仪操作简单、直接，但人力物力消耗大，而且在浅水区域容易搁浅，测量效率低，难以满足长时间序列的大面积动态变化监测需要。光学遥感数据源丰富、覆盖范围广，在水深反演中发挥了一定的作用，但受成像技术等多种综合因素影响，测量精度较低^［1-2］。机载测深LiDAR（Light Detection And Ranging）可以兼顾精度与效率，实现对岛礁周边的水深探测，但数据获取成本高，且常受空域管控等限制。星载LiDAR覆盖范围广、重复频率高、受外界因素影响小，在湖泊、冰川和海洋遥感研究中得到广泛应用^［3-4］。

ICESat-2（Ice， Cloud and land Elevation Satellite-2）是美国NASA于2018年9月发射的新一代星载LiDAR卫星，搭载的ATLAS（Advanced Topographic Laser Altimeter System）传感器采用了多波束、微脉冲、光子计数技术^［5］，具有更高的脉冲重复频率，可以获取小光斑、高密度的光子点云数据（光斑直径约17m、同轨光斑间隔仅0.7m）。同时，该系统还采用了532nm的单脉冲激光器，水体在该波段散射弱、衰减系数小，激光可穿透1倍的塞克板深度（即透明度），其最大测深能力可达38m^［6］。由于ATLAS系统发射脉冲是弱信号，受太阳背景、系统自身、大气散射、地表覆盖等因素的影响，接收的回波中包含大量噪声信号，特别是白天获取的数据更为严重。同时受水中悬浮物散射与吸收、风浪等影响，水体垂直方向不同区域的噪声分布也存在一定差异。Ma等^［7］利用机载模拟光子数据（Multiple Altimeter Beam Experimental Lidar， MABEL）信号光子和噪声光子的分布特征，基于改进的DBSCAN（Density-Based Spatial Clustering of Applications with Noise）方法提取地面和水下的表面轮廓；楼乐乐^［8］提出一种基于光子点云密度统计并结合差分回归高斯自适应的去噪算法，在浅水区域取得了较高的精度，但对复杂水底地形的适用性有待进一步探索；Chen等^［9］针对光子在水面和水柱的分布提出了一种自适应椭圆滤波测深方法（Adaptive Variable Ellipse Filtering Bathymetric Method， AVEBM），需要依据一定的规则确定滤波阈值，且深水处易出现提取错误，滤波结果仍受椭圆滤波器大小和分块大小的参数影响。

由于水面光子点分布密集而水下光子噪声多，且随水深变化，信号光子密度在沿轨方向和高程方向也发生变化，如何有效剔除背景噪声是ICESat-2/ATLAS数据获取近岸水深数据的关键。但目前光子点云去噪的研究多针对陆地类型，水域方面的研究相对较少。本研究利用南海岛礁浅海区域的ATLAS光子数据开展近海岸水深提取方法研究，重点解决复杂水底地形及水体较深位置处水深信息的提取难题，并利用机载LiDAR测深数据验证本文方法的有效性和精度。

1 研究方法

ICESat-2/ATLAS光子数据存在大量噪声，在水深信息提取应用时首先要进行去噪处理，然后将去噪后的光子数据分为水面光子和水底光子；基于两者的密度分布差异，分别采用区间估计和改进OPTICS算法进行水体信号光子提取，并结合RANSAC算法获取水面高程信息。由于光子在水体传输中会产生折射以及海面受潮力产生周期性波动，为了提高水深反演精度，需进行折射校正和潮汐校正。最后，对实验结果进行精度评价，包括利用人工标注水体信号光子，对本文精去噪结果的定性和定量评价，并与 ATL03内置去噪算法、AVEBM算法结果进行对比，以及利用机载测深LiDAR数据对提取的水深信息进行的精度评定。总体技术流程如图1所示。

图1  基于ICESat-2/ATLAS光子数据的水深提取流程

Fig.1  Bathymetry workflow based on ICESat-2/ATLAS data

1.1　光子粗去噪

通过两步完成光子去噪，即粗去噪和精去噪。粗去噪即去除原始光子数据中比较明显的噪声，利用光子点云文件中的置信度参数“signal_conf_ph”去除参数为0的光子，但去噪后仍然存在大量噪声光子。

1.2　光子精去噪

考虑到水面和水底数据的分布密度差异较大，大部分信号光子集中在水面，因此通过构建高程分布直方图，利用混合高斯分布模型对其进行两次高斯函数拟合，以两个高斯曲线的交点来确定水面光子的范围，实现水面和水底光子的分离，其中二者精去噪的方法有所不同。

1.2.1　水面光子精去噪

由于水面光子分布较为集中，可利用区间估计法去除水底光子和其他噪点。统计水面光子的高程分位数，将［0，0.02］、［0.98，1］对应的光子点作为噪点去除，从而得到水面光子数据，然后利用RANSAC（random sample consensus）直线拟合确定水面高程，计算水面光子的最低高程点并将其设为阈值，低于该高程点的光子认为是水底光子。

1.2.2　水底光子精去噪

由于水底光子分布密度不均，随水深增加，信号光子密度减小。通过改变滤波参数进行两次精去噪，最大程度提取水底信息。具体步骤如下：

（1）基于改进密度聚类算法（ OPTICS）对水底光子精去噪^［10］，采用试错法，最终选定椭圆搜索半径a和b分别为11和1，利用式（1）计算输入点的阈值参数MinP_ts^［8］：

$$\left\{\begin{array}{c}{S}_{\mathrm{1}}=\frac{\mathrm{\pi }ab{N}_{\mathrm{1}}}{hl}\\ S_{\mathrm{2}}=\frac{\mathrm{\pi }ab{N}_{\mathrm{2}}}{{\mathrm{h}}_{\mathrm{2}}l}\\ \mathrm{M}\mathrm{i}\mathrm{n}{P}_{\mathrm{t}\mathrm{s}}=\frac{\mathrm{2}{S}_{\mathrm{1}}-S_{\mathrm{2}}}{\mathrm{l}\mathrm{n}\left(\raisebox{1ex}{$\mathrm{2}{S}_{\mathrm{1}}$}\!\left/ \!\raisebox{-1ex}{$S_{\mathrm{2}}$}\right.\right)}\end{array}\right.$$

（1）

式中：S₁是椭圆搜索区域内噪声光子和信号光子的总数；N₁是给定样本的总光子数；h和l分别是样本光子点云中垂直距离和沿轨距离的最大差值；S₂是椭圆搜索区域内包括的噪声光子数。假设最低水深处向上5m的范围为噪声光子，N₂为该范围内的光子数，则h₂为5m。一般将沿轨方向10 000光子点定为样本点，实际应用中水下光子总数往往小于10 000，故通常是将所有光子数据参与计算。

（2）进一步区分信号和噪声光子。OPTICS算法引入了两个特征，即核心距离（Core Distance， CD）和可达距离（Reachability Distance， RD）。C_d表示核心距离，是对于输入的MinPts使该点p成为核心点的最小距离，而当不满足输入参数时$C_{\mathrm{d}}(p)$无意义。R_d表示可达距离，针对核心对象，若p为核心对象，R_d（p，q）表示样本点q与样本点p之间核心距离，即$\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}(p,q)$的最大值；当p不为核心对象时，$R_{\mathrm{d}}(p)$无意义。

$$C_{\mathrm{d}}\left(p\right)=\left\{\begin{array}{c}\mathrm{U}\mathrm{N}\mathrm{D}\mathrm{E}\mathrm{F}\mathrm{I}\mathrm{N}\mathrm{E}\mathrm{D},N_{\epsilon }\left(p\right)<\mathrm{M}\mathrm{i}\mathrm{n}{P}_{\mathrm{t}\mathrm{s}}\\ \mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\left(x,N_{\epsilon }\left(p\right)\right),N_{\epsilon }\left(p\right)\ge \mathrm{M}\mathrm{i}\mathrm{n}P\mathrm{t}\mathrm{s}\end{array}\right.$$

（2）

$$R_{\mathrm{d}}\left(q,p\right)=\left\{\begin{array}{c}\begin{array}{l}\mathrm{U}\mathrm{N}\mathrm{D}\mathrm{E}\mathrm{F}\mathrm{I}\mathrm{N}\mathrm{E}\mathrm{D},\\ N_{\epsilon }\left(p\right)<\mathrm{M}\mathrm{i}\mathrm{n}P\mathrm{t}\mathrm{s}\end{array}\\ \begin{array}{l}\mathrm{m}\mathrm{a}\mathrm{x}\left(C_{\mathrm{d}}\left(p\right),\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{t}\left(p,q\right)\right),\\ N_{\epsilon }\left(p\right)\ge \mathrm{M}\mathrm{i}\mathrm{n}P\mathrm{t}\mathrm{s}\end{array}\end{array}\right.$$

（3）

式中：ε是邻域半径；$N_{\epsilon }\left(p\right)$是以p为圆心；ε为半径的邻域内的点数。

R_d值越大表明该光子离其他光子越远，即为噪声光子，反之为信号光子。本文采用最大类间方差法（OTSU）确定信号和噪声光子，通过对光子数据进行分块并统计每个分块中光子数N。假设前t个光子为信号光子，剩下N-t个光子为噪声；当类间方差g最大时，将此时的R_d设为阈值R_dth；当光子R_d>R_dth时，此光子记为噪声光子，反之记为信号光子。

（3）利用改进的OPTICS算法进行第一次滤波，然后对滤波结果进行纵向分块，构建高程统计直方图。假设最低水深处向上5m的范围为噪声光子，并以此为最低高度确定二次滤波范围。以直方图峰值为中心，向两边搜索高程点数小于平均高程点数的分块，以一次滤波确定的最小邻域光子数为基础，扩大椭圆搜索半径进行二次滤波。

1.3　水深信息提取

ICESat-2/ATLAS数据产品（ATL03）仅考虑了激光在空气介质中的传播。而实际应用中，激光脉冲从空气进入水体会发生折射，导致光子位置偏移，可根据斯涅尔定律对水底光子点进行折射校正，获取光子点在沿轨方向和垂直方向的偏移量，然后基于水底信号光子计算各光子点位置的水深数据。折射校正时，需要根据激光单位指向矢量方位角将光子坐标的经纬度转换到其对应的位置。此外，由于海水受到潮力影响会产生周期性波动，且光子数据提取的水深值是获取时刻的瞬时水深，为了得到基于深度基准面的水深值，本文利用T_TIDE模型^［11］选取8分潮（$M_{\mathrm{2}}、S_{\mathrm{2}}、K_{\mathrm{2}}、N_{\mathrm{2}}、K_{\mathrm{1}}、O_{\mathrm{1}}、P_{\mathrm{1}}、Q_{\mathrm{1}}$）的潮汐调和常数对水深数据进行潮汐校正。

1.4　精度评价

采用定性分析和定量指标对精去噪结果进行评价。定性评价主要基于目视识别，即通过噪声和信号光子点的分布情况，人为判断是否存在错分现象；定量指标是利用人工标注的水体信号光子，用召回率（Recall， R）、准确率（Precision， P）和F值来评价。其中召回率表示被正确提取的有效信号光子个数占原有信号光子总数的比例，准确率是指被正确提取的有效信号光子个数与提取到的有效信号光子总数的比值，F值是召回率和正确率的调和平均值。

最后，为了验证光子数据提取的水深精度，将机载测深LiDAR获取的水深数据作为真值进行比较。实际中由于光子位置与机载点云位置无法完全重合，本文将距离光子位置最近且距离值小于2m的机载测深LiDAR点云水深值作为真值，通过绘制散点图直观地显示验证结果。

2 实验结果与分析

2.1　实验区与数据

选取2019—2021年南海永乐环礁的ICESat-2/ATLAS数据（共19个条带）进行实验，包括羚羊礁、甘泉礁、晋卿岛、甘泉岛、珊瑚岛等岛礁。该区域水质清澈含沙量低，退潮时会露出大面积的礁石，台风大潮时会被海水淹没。研究区的ATL03光子数据分布如图2所示，每一条ATL03数据包括6个条带，分别为GT1L、GT1R、GT2L、GT2R、GT3L和GT3R，并记录了每个光子的经纬度、高程、光子往返时间、激光器位置和姿态角等信息，并对大多数误差进行了纠正（如固体潮汐、大气延迟等），可获取光子高度的最佳估计^［12］。

图2  实验区ICESat-2/ATLAS光子数据分布情况

Fig.2  ICESat-2/ATLAS data distribution in the study areas

机载LiDAR测深数据由中国科学院上海光学精密机械研究所于2018年8月28日获取，覆盖范围为甘泉岛和羚羊礁区域，所用设备为机载双频雷达激光测深系统Mapper5 000，最大和最小测深分别为51m和0.25m，测深精度0.23m^［13］。由于机载数据与ATLAS数据的坐标系不同，在精度评价前对其进行坐标系统转换。

2.2　三种去噪方法结果分析

ICESat-2/ATLAS光子点云噪声点分布范围较广，通过粗去噪可以去除明显的噪声点（图3）。其中20190524GT3L数据水底地势平坦，20200419GT3R

图3  光子粗去噪结果

Fig. 3  Coarse noise photons removal results

数据水底地形复杂。图4分别显示了人工标注水体信号光子、ATL03内置精去噪（置信度参数“signal_conf_ph=4”光子数据）、AVEBM和本文方法精去噪结果。其中，图4 c~f中的圆圈表示ATL03去噪和AVEBM方法对光子错提和漏提位置。可以看出，针对平缓的水底地形，两种方法均可有效去除噪声点，但ATL03内置去噪算法难以提取较深区域的信号光子，且容易在水底地形起伏较大的区域出现信号光子缺失；AVEBM对较深处的噪声光子易产生错提取。

图4  光子精去噪结果比较

Fig. 4  Comparison of fine noise photons removal results

进一步进行定量分析（表1），本文方法提取的水体信号光子可获得更高的F值，平均为97.98%，而ATL03内置精去噪算法为92.55％，AVEBM为94.78％，分别提高了约5.87％和3.38％。

表1  三种方法精去噪统计指标结果

Tab. 1  Statistical indicators of fine noise photons removal results based on the three methods

ICESat-2/ATLAS	ATL03内置精去噪算法			AVEBM			本文方法
	R	P	F	R	P	F	R	P	F
20190222GT1L	0.970	0.982	0.976	0.932	0.999	0.965	0.972	0.999	0.986
20190222GT3L	1	0.972	0.986	0.921	0.996	0.957	0.989	0.994	0.991
20190421GT1L	0.943	0.993	0.968	0.915	0.986	0.949	0.973	0.987	0.980
20190421GT3L	0.902	0.995	0.946	0.881	0.975	0.926	0.974	0.997	0.985
20190524GT3L	0.962	0.994	0.976	0.946	0.997	0.971	0.992	0.993	0.992
20190721GT2L	0.999	0.955	0.977	0.897	0.987	0.940	0.974	0.991	0.982
20190721GT3L	0.997	0.954	0.975	0.909	0.982	0.944	0.969	0.987	0.976
20190819GT1L	0.997	0.897	0.944	0.909	0.953	0.931	0.944	0.973	0.958
20190819GT2L	0.956	0.947	0.951	0.925	0.988	0.956	0.944	0.993	0.968
20191024GT1R	0.787	0.987	0.876	0.852	0.990	0.916	0.946	0.987	0.966
20191118GT1R	0.998	0.921	0.958	0.881	0.975	0.926	0.934	0.990	0.961
20191118GT2R	0.999	0.913	0.954	0.917	0.963	0.940	0.948	0.990	0.969
20200119GT1R	0.968	0.994	0.981	0.932	0.995	0.963	0.989	0.994	0.992
20200119GT3R	0.891	0.987	0.937	0.856	0.990	0.918	0.980	0.986	0.983
20200419GT1R	0.592	0.995	0.743	0.926	0.992	0.958	0.987	0.982	0.985
20200419GT3R	0.502	0.998	0.668	0.971	0.981	0.976	0.986	0.993	0.990
20200719GT2L	0.781	0.999	0.876	0.936	0.998	0.966	0.971	0.997	0.984
20200820GT3L	0.998	0.966	0.982	0.954	0.988	0.971	0.990	0.994	0.992
20210519GT2R	0.838	0.996	0.910	0.884	0.998	0.937	0.972	0.983	0.977
平均值	0.900	0.971	0.925	0.913	0.986	0.948	0.970	0.990	0.980

2.3　水深信息精度评价

利用机载LiDAR测深数据对本文光子数据提取的水深信息进行精度评价。通过获取的水面高程和折射校正后水底信号光子高程值，计算水底信号光子对应位置的水深值。以20190524GT3L数据为例，对获取的水面信号光子和经过折射校正前后的水底信号光子采用B样条曲线对其进行拟合，得到连续的水面曲线和水底地形（图5）。

图5  20190524GT3L数据折射校正前后结果图

Fig.5  Signal photons before and after refraction correction of 20190524GT3L data

由于浅海区域的细沙、碎沙较多，且机载测深受底质物以及透明度的影响，分析中对较浅区域（水深值≤0.5m）的数据进行了剔除。本文以羚羊礁区域20190524GT3L和20200820GT3L数据提取的水深结果为例进行精度评价（图6）。结果显示，最大探测深度7.83m，最大误差约1.69m，最小约0.001m，差值较大，可能因机载数据和ICESat-2/ATLAS数据位置无法完全重合所导致；决定系数R²为0.91，均方根误差RMSE为0.53m。

图6  ATLAS测深结果与机载LiDAR测量验证

Fig. 6  Accuracy assessment of the bathymetric results obtained by ATLAS and airborne LiDAR

3 结论

ICESat-2/ATLAS光子计数数据应用的关键是去除原始数据中大量的噪声光子，本文以水体清澈、受杂质影响较小的南海为实验区，将改进OPTICS算法用于水底光子信息精去噪，降低了滤波参数的敏感性。通过对一次水底滤波结果分块、扩大搜索半径等，快速、高精度提取水体深处更密集的信号光子，提高了复杂水底地形区域的水体光子提取精度。由于光子受水体及其他悬浮物散射的影响，获取的水面数据存在一定的高程差异，本文在水面精去噪后重新计算了水底数据的范围，虽一定程度上降低对浅水区域水底信息错分，但受高斯拟合曲线的影响，在水面波动较大的区域仍存在部分水面光子错分为水底光子的现象；同时二次精去噪的滤波参数不能实现自动化设置，还需进一步提高方法的自动化程度。此外，考虑到大部分近海区域水质并非如南海岛礁水质好、透明度高，本文方法的适应性还有待进一步研究，这也是目前激光雷达测深应用的难点。

作者贡献声明

习晓环：学术指导，提出选题，设计论文框架，论文审阅与修改。

王子家：整理文献，方法实现，数据整理，数据处理及精度验证，论文撰写与修改。

王成：学术指导，论文审阅与修改。