激光测高卫星波形饱和识别与测高误差改正初步研究

郭金权，李国元，裴亮，么嘉棋，聂胜

1.辽宁工程技术大学 测绘与地理科学学院,阜新 123000;2.自然资源部国土卫星遥感应用中心,北京 100048;3.山东科技大学 测绘科学与工程学院,青岛 266590;4.中国科学院空天信息创新研究院 数字地球重点实验室,北京 100094

卫星激光测高作为一种新型的地理空间信息获取手段，因其可以快速获取大范围的地表高精度三维信息而被广泛关注（李国元等，2017）。NASA于2003年发射的冰，云和陆地高程卫星/地球科学激光测高系统ICESat/GLAS（Ice，Cloud and land Elevation Satellite/Geo-science Laser Altimeter System）在其任务期获取了大量数据，取得了很好的科学应用效果（王成等，2015）。之后NASA 于2018年发射了ICESat/ATLAS，继续依靠激光测高相关技术对极地冰盖、森林树高、海冰覆盖等方面进行科学研究。中国已经发射的高分七号卫星也装备具有全波形记录和采样功能的激光测高仪（唐新明和李国元，2017；
李国元，2018）。

针对激光测高全波形数据处理国内外学者进行了广泛的研究，但大部分研究聚焦于全波形数据的分解、波形噪声滤波和全波形数据在极地、森林等方面的应用（Lefsky 等，2007；
李松 等，2007；
邢艳秋和王立海，2009；
李增元等，2015；
李国元 等，2017；
汪垚 等，2018；
李洪鹏 等，2019），对于全波形数据存在的波形饱和现象研究较少，而回波波形中不可避免地存在饱和现象，对饱和波形进行识别和误差改正非常必要。Fricker 等（2005）通过经验公式对回波波形能量进行校正，再改正激光飞行时间，改正后高程偏差从9.6 cm 减小到1.9 cm，标准差从4.9 cm 减小到3.2 cm，但是该方法只对增益值为13 的波形数据才有效，局限性较大。Kwok 等（2006）将回波能量分级，按照经验公式分别计算不同回波能量下的激光飞行时间偏差，进而改正高程，该方法在波形饱和程度不高时效果较好，对于严重的饱和波形其改正效果很差。Sun 等（2017）对回波能量、增益值、和激光飞行时间偏差进行分析建模，改正后0—40 fJ脉冲能量的最大误差为2 cm，40—140 fJ 脉冲能量的最大误差为7 cm，但其计算过程较为繁琐，每一个增益值就需计算相关的12 个参数，且仅作为GLAS 数据经验公式，对于其他卫星测高数据的改正效果尚未可知。朱世贤等（2018）提出窗宽自适应形心修正算法，计算自适应窗宽内形心位置改正激光飞行时间和距离，但其方法适用于严谨的高斯回波波形，对于卫星激光测高仪的高斯回波波形的处理效果不太理想。针对全波形数据饱和现象，本文提出了一种饱和识别方法及结合高斯拟合与波形形心的饱和改正方法，并采用ICESat/GLAS数据对算法有效性进行验证。

2.1 波形饱和原因

波形饱和是指返回脉冲的峰值功率超出了接收机的线性动态范围，导致波形扭曲。发射脉冲高能量、地表高反射率、接收机高增益值等均可导致波形饱和。

（1）发射脉冲高能量和地表高反射率。对于接收机接收的回波脉冲能量（累计的光子数）可由下列公式表示：

式中，G为增益值；
e为电子荷；
Pt为激光发射功率；
AR为接受孔径面积；
R为激光器与目标间距离；
ρ为地表反射率；
θ为地表反射面法向量与望远镜视场方向夹角；
Ta为单程大气影响；
η为光子探测效率（门华涛等，2019；
蓝晓萍等，2015）。由公式推导出回波脉冲能量与地表反射率成正相关关系，当发射脉冲能量一定时，地表高反射率是造成回波高能量的主要原因。地表高反射率多发生在晴朗天空条件下的被冰雪覆盖的极地和湖泊水面。例如，GLAS 记录的最强回波脉冲能量发生在玻利维亚的乌尤尼盐湖，干涸湖床的浅水层发生镜面反射，返回脉冲能量大约为1000 fJ，达到探测器损伤水平的十分之一。

（2）高增益值。GLAS 接收机通过接收到的原始回波能量对增益值进行线性调整以输出增益后的能量值，当卫星飞行路径上经过反射率相差较大的区域时，如从陆地到水面，接收器需要一定的时间来通过降低增益以控制电子饱和，GLAS 最大调整持续时间为0.25 s，最长距离为1.5 km，接收机调整增益的时间取决于卫星指向角，大于标称值的指向角会降低饱和概率，但会增加高程误差（Baghdadi等，2011）。

2.2 波形饱和影响

波形饱和使得波形失真，发生“平顶”现象并伴随着波形展宽。饱和导致的波形失真使得通过波形计算脉冲到达地物地表的时间产生较大误差，进而影响高程测量精度。不进行饱和校正的GLAS 高程与参考高程存在0.1—0.9 m 的误差，极端情况下误差可达1.5 m（Baghdadi 等，2011）。并且波形饱和广泛发生于GLAS 回波数据中，尤其是极地冰川地区，如图1所示。

图1 格陵兰冰盖区域不同时间段饱和率Fig.1 Comparison of saturation rates of an orbit in the Greenland ice sheet during each mission period

GLAS14 数据产品中提供了饱和度改正标签，标签取值范围为0—5，其中0 为未饱和；
1 为虽然发生饱和但是影响可忽略不计；
2 为发生饱和，可以改正；
3 为发生饱和，改正量无法计算；
4 和5为饱和改正模型不适用，数据不能被改正。GLAS为不同的增益值规定了不同的饱和电压阈值，如表1所示。但GLAS提供的标志不一定完全准确。

表1 GLAS提供的增益值与饱和电压阈值Table 1 Gain value and saturation voltage threshold provided by GLAS

峰度的概念来自于对零均值随机变量分布特征的描述，峰度用来描述随机变量概率密度函数的尖峭程度。峰度公式如下：

式中，x为统计变量；
f为变量拟合的高斯函数；
β为峰度值；
V4为四阶中心动差；
σ4为标准差的四次方（单聪淼等，2018）。以四阶中心动差V4为测定依据，将V4除以其标准差的四次方σ4，以消除单位量纲的影响，便于不同次数分布曲线的峰度比较，从而得到以无量纲数表示的相对数，即为峰度的测定值β。在实际应用中，为保证标准高斯分布的峰度值为0，通常将计算的峰度值做减3处理，本文计算波形峰度值时均进行减3处理，以统一进行评价。

峰度是和正态分布相比较而言的统计量，描述的是数据分布的陡峭程度。当β＞0时，表示分布曲线呈尖顶峰度，为尖顶曲线，说明变量值的次数较为密集地分布在众数的周围，β值越大于0，分布曲线的顶端越尖峭。当β＜0时，表示分布曲线呈平顶峰度，为平顶曲线。当β值接近于-1.2 时，分布曲线呈水平矩形分布形态，当β小于-1.2 时，则代表发生了“U”型分布。

结合峰度知识，以及对饱和波形的观察，在GLAS 波形饱和识别的基础上，利用峰度系数对未识别的波形进一步判断，以检测遗漏饱和波形：

（1）首先判断波形中是否有大于或等于饱和阈值的采样点数据，若有则直接视为饱和波形，若没有则进行下一步。

（2）判断回波数据中是否存在电压值大于最低饱和阈值0.525 V 的采样点数据，若存在则进行第3步，若不存在则视为非饱和波形。

（3）对有效波形内的采样点求取峰度系数，当峰度系数小于-1.2时，则认为其发生了饱和。

电压值0.525 V 为GLAS 提供的最小饱和电压，利用该值可避免因波形整体偏低而导致的峰度系数偏小而导致的误识别情况；
以峰度系数-1.2 作为阈值是因为均匀分布函数更加接近饱和波形，并且饱和波形在饱和区域内的采样点数值会出现上下起伏情况，形成局部的“U”型分布，导致峰度系数更一步减小。

GLAS 根据卫星发射前的地面调试，确定了激光脉冲能量Er与波形饱和造成的时间偏差Δt之间的经验公式，如下所示：

式中，时间单位为ns，能量单位为十万亿分之一焦耳fJ。在ICESat发射后的各项研究表明，该公式对于饱和程度较小的点有较好的改正效果，对于饱和程度较大的波形效果很差（Kwok等，2006）。

在发现上式公式效果较差的情况后，NASA 团队于2017年又推出了新的经验公式对波形饱和进行改正，公式如下：

式中，Er是由高斯拟合计算的返回脉冲能量；
G是VGA 增益；
并且Eth(G)是给定增益值下的饱和开始时的阈值能量；
常数a0=0.25 和b0=0.0625 控制Eth(G)附近的δttof的变化率；
并且由对测量的拟合确定。参数α(G)控制行程时间偏差相对于饱和后的接收能量的斜率。函数Eth(G)和α(G)近似为

式中，c1—c10是根据测试数据确定的常数，不同的增益值需要重新计算（Sun等，2017）。

本文以波形为依据，对原始波形进行高斯拟合，以高斯拟合波形与原始波形的交点作为界限，计算范围内高斯拟合形心位置，与原始高斯拟合波形形心位置做差计算时间偏差，进而计算饱和波形高程校正值。具体流程如下：

（1）依据上述饱和识别方法判断回波波形是否饱和，若识别为饱和波形即进行饱和改正；

（2）采用高斯拟合方法对饱和波形进行拟合，并计算其初始形心位置。对原始波形的采样点数据进行曲线拟合，控制拟合精度大于0.98，以保证后续提取交点的准确性；

（3）计算高斯拟合与原始波形的交点，选取两交点对高斯波形进一步处理。交点选取满足：1）两交点分列与高斯拟合波形中轴线两侧；
2）两交点不在原始波形受饱和影响范围内；
3）为更好地体现波形拟合偏差，两交点应在有效波形内距离高斯拟合中轴线最远。

（4）连接选取的两个交点形成线段，计算该线段与高斯拟合波形围成的封闭区域的形心位置，即图2中的阴影区域。

图2 饱和波形高斯拟合差异Fig.2 Difference between gauss fitting and saturated waveform

（5）高斯拟合形心位置与封闭区域形心位置做差即为波形饱和引起的测高时间偏差，利用该时间偏差再乘以二分之一光速即为饱和改正值。

对于饱和波形，依据波形的特性，即其波形存在“平峰”现象，保证了其存在交点，其因饱和产生的波形展宽现象保证了其存在有限个交点。对于未发生饱和的波形，因其保有较好的高斯特性，其波形在高斯拟合时会出现大量的交点，甚至重合的曲线范围，对于饱和波形而言因为高斯波形性质被饱和现象破坏，导致其在有效波形范围内只能产生有限个交点。

如图2 所示，其中红色虚线为高斯拟合波形，蓝色实线为原始波形，黑色点为两波形交点，箭头指向为本次改正所选取的交点，图2中阴影区域即为高斯拟合波形在交点范围内的区域，求取阴影范围内的高斯拟合形心。再计算与初始高斯拟合形心位置差异即为时间偏差，进而计算改正高程。

误差的正负号由形心位置决定，若交点范围区域的形心位置相较于高斯拟合形心位置前移，则证明波形饱和使得通过形心计算的高程偏低，误差改正值为正；
若范围内的形心位置相较于高斯拟合形心位置后移，则证明波形饱和使得通过形心计算的高程偏低，误差改正值为负。

5.1 研究区域介绍

选取青海湖、纳木错湖、色林错湖3 个区域GLAS 数据进行实验分析。青海湖位于青海省西北部的青海湖盆地内，既是中国最大的内陆湖泊，也是中国最大的咸水湖，总面积约4583 km2。纳木错位于西藏自治区中部，是西藏第一大湖泊，也是中国第二大的咸水湖，总面积约1920 km2。色林错地处西藏自治区申扎、班戈和尼玛三县交界处，位于岗底斯山北麓，申扎县以北，总面积约239 km2。如图2 所示，为3 个湖泊区域概况，其中红色点为受饱和影响激光点，黄色为正常激光点。

选择湖泊区域进行实验主要有以下几点原因：

（1）水体因存在耀斑现象，易发生波形饱和情况，湖泊范围内可获得多个饱和波形数据。

（2）GLAS 光斑面积较大，回波波形受到光斑内复杂地物的影响，而水体回波波形一般为单峰波形，可以更好的对因饱和产生的误差进行分析，避免其他因素影响。

（3）现实中因风力等因素水体高程值存在差异，但静态非流动的湖面在同一时间段内水面高程在小范围内基本保持一致。

图3 研究区域概览Fig.3 Overview of the research region

5.2 实验数据介绍

在选取实验用GLAS数据时应满足以下原则：

（1）因为实验区域皆为高原湖泊，为避免结冰期与风期对高程数据的影响，应选用收集时间为4月—10月的GLAS数据。

（2）湖泊与陆地相接处易发生饱和情况，但是湖岸边复杂的地形情况导致无法获得高程的近似真值，所以选取的饱和数据应在湖泊深水区域内。

（3）每组内饱和点数不应超过总点数的三分之一，各未饱和数据高程值波动程度小，以保证通过未饱和数据计算的高程真值的准确性。

依据以上原则，本文实验共选取4组数据，分为A、B、C、D 共4 组。组内数据按A1、A2、A3至A12 排序，每组数据总长度约为2 km。A 组为2004年8月青海湖GLAS 数据，A 组数据内包含3 个饱和点，分别为A1、A2、A5。B 组为2008年8月青海湖GLAS 数据，B 组数据内包含3 个饱和点，分别为B7、B8、B9。C 组为2007年4月纳木错GLAS 数据，C 组数据内包含3 个饱和点分别为C4、C5、C6。D 组为2005年6月色林错GLAS 数据，D 组数据内包含3 个饱和点分别为D9、D10、D11，其中D9为本文算法探测到的GLAS未注明为饱和的激光点。

6.1 饱和识别结果与分析

对4组数据进行饱和识别实验，按前文方法先判断回波波形是否有采样点电压值大于GLAS 饱和阈值的情况，有则视为饱和波形，然后判断最大采样点电压值是否大于0.525 V，小于则视为未饱和波形，大于则计算其有效波形峰度系数，若峰度系数值小于-1.2 则视为饱和。经过试验发现，本文方法对A、B、C 3 组数据进行饱和识别结果与GLAS 提供结果一致，D 组中D9 点GLAS 未认为发生饱和，但通过本文方法识别出来，如图4所示。

图4 波形饱和识别及误差改正实验流程图Fig.4 Flow chart of waveform saturation identification and error correction experiment

D9 点波形最大电压值大于0.525 V，峰度系数为-1.62小于阈值-1.2，为本文方法发现的GLAS未将其标志为饱和的数据点，下面我们利用饱和波形相关特性对其准确性进行评估。首先D9 点波形最大电压值为1.36 V，明显高于同组内其他未饱和激光点波形最大电压值，且与GLAS 认为饱和的D10、D11 点最大电压值十分接近；
以波形形状来说，D9点出现了明显的“平峰”，且“平峰”范围较大，符合饱和波形形状特点。

6.2 饱和改正实验及精度分析

使用GLAS14 提供的地面高程值作为原始高程值，以每组内未饱和数据的平均值作为水面高程近似真值，对于饱和数据分别将GLAS 数据提供的饱和校正量与本文算法计算的饱和校正量加入到原始高程值中，计算饱和改正后高程值，再与未饱和数据计算的近似水面高程值进行比较，以评价精度。

改正前后数据均值、饱和波形改正前后高程对比、改正后误差对比如下表3、4所示。

表3 改正前后数据均值比较Table 3 Comparison of data mean before and after correction /m

在共12 个饱和点数据中，D10 数据点在改正后与近似真值差距大于0.1 m，B9 点改正后偏离近似真值，下面对这两点进行分析：

D10 点本文改正量为0.041 m，改正后与近似真值偏差为0.124 m，虽然相较于GLAS 改正后的偏差值0.153 m 较小，但是相对于近似真值来说误差较大。但GLAS 提供的D10 点饱和改正量仅为0.012 m，即GLAS认为D10点由饱和引起的测量误差值不大，所以我们推测D10 点与近似真值之间的误差中不仅仅包含由饱和引起的偏差，可能还包含其他原因造成的偏差，导致仅仅利用饱和改正无法获得理想结果。

B9 点本文给出的改正量为0.027 m，GLAS 给出的改正量为0.167 m。相较于GLAS 给出的改正数值，本文给出的值要小很多，虽然改正后B9高程值更加偏离近似真值，但是差值为0.052 m，仍在可接受的范围内。我们推测B9 点的高程可能因为水面风浪等原因，造成局部水面的升高，以致于高程值高于利用均值计算的近似水面真值。

总体来看，本文算法饱和改正最大值为0.136 m，饱和改正最小值为-0.008 m，改正后的误差的绝对值均值为0.045 m，相较于GLAS 改正量的误差均值0.119 m 有了很大提高，改正后精度可达0.05 m

6.3 特殊波形分析

D9 点是本文识别算法发现的GLAS 遗漏的饱和波形，下面对其进行分析

首先，D9点的波形形状满足饱和波形的特征，从图5（a）中可以看出，D9 点波形出现明显的“平峰”，并且波形相较于未发生饱和的D12点明显展宽，且波形幅值与GLAS 探测到的发生饱和的D10、D11 点十分接近，高于未饱和波形，所以认为D9点发生饱和，并且被本文识别算法正确识别。

图5 D组数据饱和波形识别结果Fig.5 Identification results of saturated waveform in group D

表2 波形数据峰度系数Table 2 Kurtosis coefficient of waveform data

利用本文算法进行饱和改正，得出了-0.008 m的改正量，即波形饱和导致高斯拟合的形心位置前移，使得测量高程值加大。

这与前人的研究相反，以前的研究都认为波形饱和会使得高斯拟合位置后移而使得测量高程偏低，但是前人的研究都是以GLAS 提供的饱和标志作为依据，仅仅对GLAS 探测到的饱和波形进行分析，而D9 点GLAS 并不认为其发生波形饱和，而是通过本文的饱和识别方法筛选出来的，所以与前人的研究有差异也是可能存在的，并且可信的。

表4 饱和波形改正后比较Table 4 Elevation comparison after correction of saturated waveform /m

图6 改正前后高程数据对比Fig.6 Elevation data comparison before and after correction.

作为水面上的激光点，其高程应该与相邻激光点在小范围内保持一致，由于D10、D11 点发生饱和，其原始高程不可信，将D9 点与相邻的D7点、D8 点高程相比较，以评价饱和改正量是否可信。

从表6 中可以看出，D9 点高程不仅高于以组为单位计算的近似真值，并且高于相邻的D7 点、D8 点高程，所以本文认为GLAS 给出的原始高程值要高于D9 点真实高程值，本文改正方法给出的-0.008 m 的饱和改正值也是可信的，虽然改正量偏小，但是改正后高程值确是更加接近于真实值。根据D9 点的情况认为波形饱和是在大多数情况下导致测量的高程值偏低，但不是绝对的，存在少量高程值偏高的情况。

表6 D7点，D8点，D9点原始高程值Table 6 D7 point，D8 point，D9 point original elevation value

图7 D9点高斯拟合差异Fig.7 Difference of gaussian fitting at D9 point

针对激光测高仪回波饱和现象，本文结合饱和阈值与波形峰度系数对饱和波形进行了有效识别；
并利用原始波形数据形状特征信息，计算形心位置差异，进而确定饱和改正值，以3个湖泊区域的4 组12 个GLAS 数据进行实验，证明了算法的有效性。

相较于前人的研究，本文提出的饱和识别算法更加准确，发现了GLAS 利用固定阈值识别所遗漏的饱和波形；
本文提出的饱和改正算法相较于Sun 等（2017）提出的方法更加简单、易于实现，相较于朱世贤等提出的方法改正后精度更高。且仅依靠波形数据形状特征进行计算，不依靠经验公式，所以更具有普遍适用性。与前人研究不同的是，本文实验中还发现波形饱和对高程测量影响不总是偏低，也会造成偏高的情况，并通过改正后的高程与真实高程的比较证明了该种情况确实存在。

中国已经发射高分七号并计划发射陆地生态系统碳监测等搭载全波形激光测高仪的国产卫星，通过对饱和波形数据的有效改正可以大大增加可用激光点数目及精度。本文所提方法对国产激光测高仪回波数据饱和处理能提供一定参考，但主要适用于湖泊、裸地、极地冰盖等单波形饱和数据的有效地识别和改正，对于森林等复杂区域多波形叠加的饱和数据如何识别及改正仍待研究。

猜你喜欢 高斯波形高程 红树林宜林地滩面高程及潮水退干时间时长的测量方法绿色科技(2022年16期)2022-09-15基于SVD的雷达嵌入式通信波形设计方法航空学报(2022年7期)2022-09-05正面碰撞车身加速度对乘员腿部损伤的影响汽车零部件(2022年7期)2022-08-038848.86m珠峰新高程当代陕西(2020年23期)2021-01-07通用6T系列变速器离合器鼓失效的解决方案汽车与驾驶维修（维修版）(2020年2期)2020-01-07数学王子高斯小天使·二年级语数英综合(2019年4期)2019-10-06基于二次曲面函数的高程拟合研究消费导刊(2017年8期)2018-01-18基于MATLAB的GPS高程拟合程序设计建筑建材装饰(2017年6期)2017-07-04从自卑到自信 瑞恩·高斯林电影故事(2015年16期)2015-07-14任意波形编辑软件中手动任意绘制功能的设计现代电子技术(2009年8期)2009-06-25

推荐访问:波形 误差 饱和