基于统计直方图两步法的星载单光子数据去噪

摘要

针对星载单光子激光有效回波信号混杂于噪声中难以区分的问题，提出一种基于统计直方图两步法的星载单光子数据去噪方法。先后采用沿轨小窗口直方图粗去噪与距离平方统计直方图精去噪，实现星载单光子回波数据中的噪声光子有效剔除。利用该方法对美国星载单光子激光雷达先进地形激光测高系统(ATLAS)的强与弱波束、白天和夜间、平地与山地3种典型情况下的回波光子数据进行实验，结合ATLAS官方去噪结果，基于混淆矩阵统计去噪精确度。实验结果表明，强波束数据去噪精确度为98.86%，弱波束数据去噪精确度为96.94%；夜间数据去噪精确度为99.02%，白天数据去噪精度为98.86%；山地数据去噪精确度为96.28%，平地数据去噪精确度为96.94%。说明本文方法适用于常见的以上3种典型情况下的星载单光子数据去噪。

关键词

星载单光子数据; 先进地形激光测高系统; 统计直方图; 去噪精确度; 混淆矩阵

卫星是我国空间资产的重要组成部分，在测绘导航等方面有巨大的应用价值^[1]。随着星载激光雷达硬件技术的不断提升，具有微脉冲、高重频、光子探测水平的新一代星载激光雷达已成功用于对地观测。2018年9月美国成功发射的第二代冰、云和陆地高程卫星(Ice, Cloud, and land Elevation Satellite-2，ICESat-2)搭载了全球首台星载单光子激光雷达先进地形激光测高系统(ATLAS)^[2]，其采用强弱3对波束同步对地测量^[3]，目前已获取大量的白天/夜间，平地/山地，强/弱波束等多类情况下的回波光子数据，并广泛用于冰盖、海冰、地形和树冠高以及内陆水等多个领域。

星载单光子激光雷达具有多波束、低耗能、高重频等明显优势，但受太阳背景噪声、大气、硬件自身暗计数等影响，其回波光子数据中存在大量的噪声，有效回波信号光子淹没其中。因此，星载单光子数据去噪成为单光子数据深加工和应用前的首要工作。研究学者针对ATLAS数据开展了多种去噪实验，文献[4-5]提出一种改进空间密度聚类的方法，对森林地区单光子数据开展了去噪实验；文献[6]则提出方向自适应算法，并对山地地形开展相关实验；马跃等先后提出基于离群因子的去噪方法^[7-8]，对水体开展有效光子信号提取。以上去噪方法取得了一定成果，但主要是针对单一地形或情景下的，目前尚未同时针对强/弱波束或白天/夜间或不同地形等情况进行实验。

本文针对单光子背景噪声与有效信号光子混杂难以区分的问题，提出一种基于统计直方图两步法的星载单光子数据去噪方法，对强与弱激光波束、白天与夜间、平地与山地3类典型情景的星载单光子回波数据开展实验。首先根据沿轨距离与高程，将原始数据分割成多个连续局部小窗口，基于小窗口进行频数直方图统计实现粗去噪；统计剩余所有光子的距离平方频数直方图，基于正态分布特性，实现星载单光子点云数据的精去噪；最后，结合ATLAS官方去噪结果，分别统计本文方法在强光束与弱光束、白天与夜间、平地与山地3类典型情景下的单光子数据去噪精确度。

1 单光子数据去噪方法原理

星载单光子数据具备三维地理信息，单波束光子数据沿飞行方向呈条带状点云分布，因此在实际光子去噪时，通过降维方式，将三维地理坐标光子数据转换为“沿轨距离-高程”的二维数据，进而开展单光子激光数据去噪。去噪方法流程如图1所示。

图1  星载单光子数据去噪流程图

Fig.1  Denoising flow chart of spaceborne single photon data

1.1 “沿轨时间”到“沿轨距离”转换

星载单光子激光记录了每个光子的返回时刻及光子时间，根据光子时间与速度可将光子数据从“沿轨时间-高程”转换为“沿轨距离-高程”，常用转换方法主要为以下2种：

基于卫星固定速度转换^[9]：

$$d=v_{\mathrm{0}}{t}_{\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}}$$

(1)

式中：$v_{\mathrm{0}}=\mathrm{7}\mathrm{k}\mathrm{m}/\mathrm{s}$为卫星固定速度常量；$t_{\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}}$为光子时间。

基于卫星实时速度分量转换^[10]：

$$d=\sqrt[]{{v_{\perp }}^{\mathrm{2}}+{v_{//}}^{\mathrm{2}}}{t}_{\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}}$$

(2)

式中：$v_{\perp }$为卫星北方向速度分量；$v_{//}$为卫星东方向速度分量。显然，第2种转换方法较第1种方法精确度更高，更可靠，故本文采用第2种方法进行数据转换。

1.2 基于沿轨小窗口直方图统计的粗去噪

地表返回的有效信号光子与噪声光子密度在水平方向和垂直方向上均显著不同，综合地形起伏等因素，本文提出基于局部小窗口直方图的粗去噪。即将目标数据按沿轨距离与高程2个维度细分为多个小窗口，统计每个小窗口内光子数据，实现粗去噪。

首先将目标数据按固定窗口$D_{\mathrm{a}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}$划分为$N$个沿轨距离窗口，随后对单个沿轨距离窗口按高程窗口$D_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\_\mathrm{d}\mathrm{i}\mathrm{r}}$划分为M个高程子窗口。高程窗口与沿轨窗口关系如下：

$$D_{\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\_\mathrm{d}\mathrm{i}\mathrm{r}}=D_{\mathrm{a}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}\times \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{2}\alpha$$

(3)

式中$\alpha$为沿轨窗口内地形粗略坡度。沿轨距离窗口由地形坡度确定，可随坡度的增大而缩小，通常设置为10~100 m。

随后，统计单个沿轨局部窗口内部每个高程子窗口内的光子数量，构建“高程—光子数”频数直方图。提取直方图中最大3个频数，默认直方图中频数最大的高程子窗口内光子数据为有效信号，并根据式(4)判断另外2个频数对应高程窗口光子数据是否为有效信号。

$$signal\_photon=\left\{\begin{array}{l}S\_N_{\mathrm{1}},N_{\mathrm{2}}<{\sigma }_{\mathrm{1}}{N}_{\mathrm{1}}\\ S\_N_{\mathrm{1}}+S\_N_{\mathrm{2}},N_{\mathrm{2}}{\sigma }_{\mathrm{1}}\ge N_{\mathrm{1}}\&N_{\mathrm{3}}<{\sigma }_{\mathrm{2}}{N}_{\mathrm{1}}\\ S\_N_{\mathrm{1}}+S\_N_{\mathrm{2}}+S\_N_{\mathrm{3}},N_{\mathrm{3}}\ge {\sigma }_{\mathrm{2}}{N}_{\mathrm{1}}\end{array}\right.$$

(4)

式中：$signal\_photon$为单个沿轨窗口内有效信号光子数据；N₁, N₂, N₃为该沿轨窗口内最大3个频数(N₁≥N₂≥N₃)；S_N₁, S_N₂ , S_N₃为3个频数对应的高程区间内光子数据；σ₁,σ₂为经验值。

1.3 基于距离平方统计直方图的精去噪

粗去噪后，光子数据中仅含少量噪声，剩余有效信号光子之间距离小于有效信号光子与噪声光子或噪声光子与噪声光子之间的距离。为显著区分噪声光子，本文基于式(5)计算每个光子的K邻域距离平方的均值作为衡量标准。

$$Dist\_av{e}_K(i)=\frac{\mathrm{1}}{K}\sum _{i\in N_K(i)}\left[(x_i-x_j{)}^{\mathrm{2}}+(h_i-h_j{)}^{\mathrm{2}}\right];i\in D$$

(5)

式中：$i$为粗去噪后剩余光子数据集D中的每个光子对象；$Dist\_av{e}_K(i)$为第$i$点与其K邻域内光子距离平方的均值；$x_i$,$x_j$分别为光子i, j沿轨距离；$h_i$,$h_j$为光子i, j高程。

将统计数据集D内所有光子的K邻域距离平方的均值，按其范围划分为dn个间隔，统计每个间隔内的光子数量，绘制距离平方频数直方图，该直方图近似呈单侧正态分布，故同样满足正态分布特性^[11]。以$\mathrm{1}\times \sigma$置信度作为信号与噪声的分界点，对应有效信号光子置信区间如式(6)所示：

$$Inter\_Dis{t}_{\mathrm{S}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{l}}=\left\{{\displaystyle \sum _{p=\mathrm{1}}^{dn}}Nu{m}_{p\in dn}(p)<\sigma \times Tota{l}_{\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}}\right\}$$

(6)

式中：$Inter\_Dis{t}_{\mathrm{S}\mathrm{i}\mathrm{g}\mathrm{n}\mathrm{a}\mathrm{l}}$为有效信号光子区间；${\displaystyle \sum _{p=\mathrm{1}}^{dn}}Nu{m}_{p\in dn}(p)$为距离平方频数直方图第$p$个的频数；$\sigma =\mathrm{0.682}\mathrm{6}$；$Tota{l}_{\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}}$为数据集$D$内所有光子数量。有效信号区间外所有光子数据均认为是噪声光子。

1.4 基于混淆矩阵的精确度验证

参考遥感图像分类精确度评价方法^[12]，本文采用混淆矩阵来统计和分析单光子去噪精确度，主要通过将分类后每个光子属性与真实光子属性比较，统计去噪结果，如式(7)所示。

$$AC=\frac{TP+TN}{TP+FN+FP+TN}$$

(7)

式中：AC为去噪精确度；TP为算法判定与实际均为有效信号光子的总数；TN为算法判定与实际均为噪声的光子总数；FN为真值为有效信号光子但算法判定其为噪声的光子总数；FP为分类为有效信号光子但真值为噪声的光子总数。

2 实验数据与实验结果

2.1 实验数据

为验证本文算法的可行性，实验选取ATLAS全球地理定位数据ATL03级产品，该数据包括去噪前的光子数据的经度、纬度和椭球高。光子数据的差异性主要是受地形特征和观测时间为白天或夜间的影响，因此，本文选取地形较为平坦的内蒙古苏尼特右旗以及地形较为陡峭的贺兰山脉作为实验区域。其中，在苏尼特右旗实验区，主要开展强/弱波束，以及白天/夜间光子实验，具体数据如表1所示；在贺兰山实验区，开展夜间数据实验，并和平坦地形对比。最后将本文算法实验结果结合NASA官网结果，评价算法的精确度。

表1  实验数据列表

Table1  The list of experimental data

	study area	time(Beijing)	data type	data name
experiment 1	Sunite Right Banner	11:22	strong beam	ATL03_20190929032220_00340506_003_01_gt3r
			weak beam	ATL03_20190929032220_00340506_003_01_gt3l
experiment 2	Sunite Right Banner	23:14	night	ATL03_20190926151450_13830402_003_01_gt1l
		11:22	daytime	ATL03_20190929032220_00340506_003_01_gt3r
experiment 3	Sunite Right Banner	11:22	flat ground	ATL03_20190929032220_00340506_003_01_gt3l
	Helan Mountain	21:51	mountain land	ATL03_20191105135133_06060502_003_01_gt1r

第1组实验所用数据相关信息如图2所示，该地区地势相对平坦，最大高程差约20 m。将ATLAS数据由“沿轨时间—高程”转换为“沿轨距离—高程”，转换后结果如图2(b)所示。根据经纬度截取沿轨长度14 km作为实验数据。

图2  苏尼特右旗光子数据

Fig.2  Photons data of Sunite Right Banner

2.2 实验与分析

1) 强/弱波束实验结果

ATLAS为适应不同测量目标的反射率，减少由单次回波光子数过多而导致的地表反射率反演失真等问题，采用3对6波束，每对波束包含能量比为1:4的强/弱波束。本文选取表1中的强/弱波束进行实验，首先对数据进行粗去噪处理。实验中沿轨距离窗口$D_{\mathrm{a}\mathrm{l}\mathrm{o}\mathrm{n}\mathrm{g}\_\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{c}\mathrm{k}}$设置为100 m，α=5°，式(4)中的${\sigma }_{\mathrm{1}}$,${\sigma }_{\mathrm{2}}$分别设置为0.9和0.85。结果如图3所示。其中蓝色的点表示被剔除的噪声光子，红色的点表示粗去噪识别的有效信号光子。从图中可以看出，本文提出的粗去噪方法可以剔除偏离有效信号光子较远的明显噪声，为精去噪提供数据基础。

图3  粗去噪结果

Fig.3  Rough denoising results

在粗去噪结果的基础上，对剩余光子进行局部距离平方直方图统计。经实验，当式（5）中K设置较大时，会造成过去噪，导致光子信息不连续；设置较小时，会导致噪声光子过度保留，因此设置K=5，即统计距离每个光子最近的5个光子的距离平方均值，实现精去噪。结果如图4所示，从图4（b）看出，尽管弱波束精去噪后结果可以清晰表达出地形起伏情况，但仍存在较为明显的尚未剔除的“毛刺”。相比之下，强波束精去噪结果较好。为直观对比2种情况去噪水平，图4(c)~4(d)为同一区域NASA给出的光子置信度和本文算法处理结果。不同的颜色代表不同置信水平，选取每个光子最大的置信水平代表其置信度，黑色点代表缓冲区，红色点代表置信水平为低的光子，蓝色点代表置信水平为中的光子，绿色点代表置信水平高的光子，青蓝色点代表本文算法识别的有效信号光子。可以看出，本文算法识别的光子基本包含中高置信水平的光子。

图4  苏尼特右旗白天数据去噪结果及与官方对比

Fig.4  Denoising results of daytime data in Sunite Right Banner and comparison with official results

为定量评价本文算法的去噪精确度，以及对比算法针对强/弱波束的不同效果，将本文算法与NASA官网结果同时判定为有效信号光子的数据作为真值，获得强/弱波束光子以及真值的混淆矩阵，如表2所示。通过式(7)统计可得(文中所有混淆矩阵TP、TN、FN、FP位置如表2所示)，强波束的去噪精确度为98.86%，弱波束的去噪精确度为96.94%。由此可见，强波束具有更好的去噪效果。

表2  强波束、弱波束光子以及真值的混淆矩阵

Table2  Strong beam, weak beam photon and confusion matrix of true value

data		signal	denoise	counts
strong beam	signal	3 960 (TP)	347 (FP)	4 307
	denoise	0 (FN)	26 226 (TN)	26 226
	counts	3 960	26 573	30 533
weak beam	signal	10 788 (TP)	532 (FP)	11 320
	denoise	665 (FN)	27 226 (TN)	27 891
	counts	11 453	27 758	39 211

2) 白天与夜间实验结果

光子噪声主要包括太阳背景与大气噪声，往往夜间数据信噪比较高。本文采用苏尼特右旗2019年9月26日的数据进行实验，其中ATL03_20190926151450_13830402_003_01(gt1l)为强波束数据，ATL03_20190926151450_13830402_003_01(gt1r)为弱波束数据。图5可以明显看出夜间数据噪声较少(图中蓝色的点表示噪声光子，红色的点表示本文算法识别出的有效信号光子)，因此直接采用精去噪的方法进行实验，并将其结果与官网提供结果进行对比。

图5  苏尼特右旗夜间数据去噪结果及与官方对比

Fig.5  Denoising results of nighttime data in Sunite Right Banner and comparison with official results

从图5(d)标注的位置可以看出，官方判定为置信水平高的光子疑似噪声光子，而本文算法能够较好地将其剔除。将其与第1组实验中的白天数据比较，夜间数据由于背景噪声率较低，去噪效果比白天好，且精去噪结果的连续性也好于白天的数据。为量化评价夜间和白天情况下算法去噪精确度，统计夜间强波束与真值的混淆矩阵，如表3所示。利用式(7)得到去噪精确度为99.02%，因此可以看出，在夜间噪声背景噪声率低的情况下，同白天强波束对比，其去噪精确度较高于白天的光子数据。

表3  夜间与真值的混淆矩阵

Table 3  Confusion matrix of nighttime strong beam and true value

data		signal	denoise	counts
strong beam	signal	8 903	36	8 939
	denoise	57	494	551
	counts	8 960	530	9 490

3) 平地与山地对比

上述实验均选取地形较为平坦的区域，为验证本文算法在山地的适应性，选取贺兰山区域光子数据进行实验并分析，其中ATL03_20191105135133_06060502_003_01_gt1l(gt1l)为弱波束数据，ATL03_20191105135133_06060502_003_01_gt1r(gt1r)为强波束数据。挑选地形起伏较大的区域，对该数据进行去噪处理实验中，设置K=10，结果如图6所示。无论是强波束还是弱波束，在1.326×10⁶ m附近，均出现噪声密集从而被判定为信号的情况。但在其他地方，即使山脊或山谷，本文算法也可以较好地识别信号光子。图6(c)~(d)为与NASA官网对比结果。针对各实验数据分析，当波束能量较弱时，会引入更多的背景噪声，而夜间数据由于包含较少太阳辐射噪声及大气散射噪声，能够采取较为简单的方法剔除噪声。采用与第一、第二组类似的方法获得山地弱波束与真值的混淆矩阵，如表4所示，经计算可知，陡峭地形的去噪精确度为96.28%，相较平地弱波束数据的去噪结果，去噪精确度较差。

图6  贺兰山实验结果

Fig.6  Experimental results of Helan Mountain

表4  山地弱波束与真值的混淆矩阵

Table4  Confusion matrix of mountain weak beam and true value

data		signal	denoise	counts
weak beam	signal	47 014	1 038	48 052
	denoise	886	2 807	3 693
	counts	47 900	3 845	51 745

3 结论

本文为解决星载单光子激光雷达有效信号光子与噪声难以区分的问题，提出了一种基于统计直方图两步法的星载单光子去噪方法。为验证算法的适应性，选取强波束与波束、白天与夜间、平地与山地3类对立典型情况下的回波光子数据进行实验。选取以上3类典型情况下的ATL03数据，采用直方图统计的方法先后进行沿轨小窗口粗去噪与整体距离平方精去噪实验，得到结论如下：

1) 相比于弱波束数据，强波束数据去噪精确度更高，白天情况下，强波束去噪精确度达到98.86%，弱波束去噪精确度达到96.94%；

2) 相对于白天情况，夜间数据去噪效果更好，针对强波束数据进行分析，白天数据去噪精确度为98.86%，夜间数据去噪精确度为99.02%；

3) 针对平地与山地地形数据，在弱波束情况下，平地地形数据去噪精确度高于山地数据，分别为96.94%和96.28%。

综上，本文方法在强、弱波束与白天、夜间以及平地、山地情况下均开展了去噪实验，在3类典型对立情况下整体去噪精确度均较高，可以看出本文方法适应性强，可用于星载单光子激光绝大多数情况下的回波数据噪声滤除。