一种渐进剥离与整体拟合的子波形提取方法

谭建伟，程春泉，王志勇，徐志达

(1.山东科技大学 测绘科学与工程学院，山东 青岛 266590；2.中国测绘科学研究院，北京 100036)

0 引言

冰、云和陆地高程测量卫星(ice，cloud and land elevation satellite，ICESat)搭载地球科学激光测高系统(the geoscience laser altimeter system，GLAS)，是全球首个具备全波形记录和采样功能的大光斑激光雷达测高卫星[1]。其主要科学任务是监测极地冰盖变化情况，另外也用于准确估计全球地区的森林结构特征，如林区冠高、生物量和碳含量等[2]。美国航空航天局(National Aeronautics and Space Administration，NASA)于2018年发射的ICESat-2卫星进一步监测海冰厚度变化，并收集有关森林生长和地形高度的数据[3]。此外，作为我国高分辨率对地观测系统重大专项之一的高分七号已经发射升空。该星具备优于1 m分辨率的立体观测能力以及精度优于1 m的激光测高能力，主要用于1∶50 000、1∶10 000比例尺地图测图和更新，以及海冰和陆地植被信息等资源的调查工作。为更好了解并利用卫星测高数据产品进行极地冰盖变化监测、陆地资源调查以及提升自主处理星载激光雷达数据的能力，对星载全波形数据进行数据处理研究工作显得非常必要。

由于激光雷达具有穿透植被的特性，回波样本中记录的是多个目标反射波形相互叠加的结果[4]。为更好获取地物特征参数，需对回波波形进行分解，中外学者对波形分解的方法进行了大量而深入的探讨。Hofton等[5]提出经典的高斯分解算法，通过相邻奇偶拐点确定高斯分量，该方法简单有效，但受噪声影响较大，会产生许多“伪”拐点。王成等[6]采用二阶偏导数求取拐点来解算高斯分量的中心位置、振幅和半宽，并对最大振幅处的高斯分量位置和半宽予以调整，用来削弱由于饱和前向散射使波形曲线扭曲的影响，最后采用非线性最小二乘来拟合优化，其分解精度的高低依赖于高斯分量个数初始估计的准确度。赖旭东等[7]采用逐层剥离的策略不断分解出波形分量，直到最大峰值小于设定的阈值，该方法分解速度快，对弱波的探测能力强，但其利用1/2峰值位置确定半宽值的方法，导致在连续波形震荡处往往得不到正确结果。李国元等[8]提出一种基于波峰自动识别的全波形高斯分解算法，将峰值检测法和奇偶拐点法进行结合，大大提高了分解的速度和精度，但快速获取有效拐点依然受到限制。段乙好等[9]提出一种高斯拐点匹配法，利用平面曲线离散点、集拐点的快速查找算法检测拐点，通过相邻左、右拐点来确定高斯分量，该方法能够删减大量的“伪拐点”，对弱波的检测能力较强，但其分量峰值位置会产生一定的偏移。

针对目前全波形分解方法存在的“伪”拐点数量多、弱子波形提取困难以及拟合精度不高等问题，本文提出一种在总波形约束下的子波形渐进剥离分解方法。该方法能快速、准确探测分量参数初值，有效削弱次波对主波的叠加效应影响，对弱波也有一定的探测能力，同时在总波形的约束下进行所有提取子波形的整体拟合，拟合精度较高。

1 次波对主波的影响分析

复杂地物的反射信号通常是一个连续震荡的波形，子波形之间的相互叠加效应使真实地物回波信号产生展宽及峰值位置偏移等现象。子波对主波强度的瞬时影响绝对量和相对量与该时刻次波的强度及次波与主波的相对强度直接相关，离子波交点越近的时刻，主波波形的相对和绝对影响均越大。也可认为，波形拐点和峰值位置的偏移程度与主、次波之间的叠加面积大小成正比，次波与主波之间的重叠面积越大，对主波的拐点和峰值位置影响也就越大；反之，则越小。

如图1所示，当次波位于主波的左侧时，峰值位置左侧的拐点位置受影响程度较大，往左偏离了主波波形，而峰值位置右侧拐点偏移程度很小；同样当次波位于主波的右侧时，对峰值右侧拐点位置影响较大，因此根据整体波形的奇偶拐点测定各波形分量的实际峰值位置和半宽值误差较大，会产生偏移和波形展宽等现象。本文透过次波对主波的影响探讨子波形峰值位置偏移和波形展宽的原因，从提升波形参数的初始精度和减少“伪”拐点的错误数2个方面入手，利用受次波影响较小的主波峰值和近峰值拐点代替左、右拐点提取波形参数初值，并与峰值配合多种策略，减少“伪”拐点的数目，稳健提取各级子波形，实现高精度整体拟合。

图1 次波对主波的影响

2 子波形的分解与拟合

2.1 高斯波形及其预处理

一般认为，复杂地物的接收激光脉冲信号是若干个单一地物反射的高斯信号相互叠加的结果，接收器接收的返回激光脉冲可认为是发射脉冲与探测目标在时间域上的卷积[10]。经典的高斯函数模型和叠加情况如式(1)、式(2)所示。

(1)

(2)

式中：φk(xi)为第k个高斯分量函数；Ak、Tk和σk分别为第k个高斯分量的振幅、中心位置和半宽值；N为高斯分量个数，一般不大于6；noise为回波信号的背景噪声值。

波形复原是将原始回波记录样本中反向记录的数据通过波形倒置并进行电压值转换，将其以真实强度值正向输出来。本文采用5×5的高斯滤波模板从左往右逐步扫描，将邻域内的加权平均值作为模板中心点的幅值来平滑滤波整个回波信号，以提高波形的信噪比[11]，利用首尾各20帧波形数据的均值作为背景噪声均值并统计其标准差，以2、3倍标准差与噪声均值之和作为整个波形的振幅阈值[12]，这种方法对符合高斯分布的波形尤为适用。回波信号的去噪效果如图2所示，其中，横坐标序为等间距1 ns的时间序列；纵坐标为波形脉冲的振幅大小；绿线表示原始回波信号；蓝线表示去噪后波形；红线表示振幅阈值线。

图2 GLA01波形去噪图

2.2 子波形的渐进剥离

由于GLAS回波信号是以横坐标为等间距1 ns的离散点集构成，本文采用平面曲线离散点集拐点的快速查找方法[13]来探测拐点的位置及数量。在拐点处作一条正向切线，其左右凹凸曲线上的离散点集都分别位于切线的两侧，因此需要4个点就能判断拐点的位置。设4个点的坐标分别为(xk-2，yk-2)、(xk-1，yk-1)、(xk，yk)、(xk+1，yk+1)，当回波信号相邻4点坐标满足式(3)时，(xk，yk)则为拐点。表达式如式(3)所示。

(yk-2+yk-2yk-1)*(yk-1+yk+1-2yk)<0

(3)

通过以上算法可探测整个回波信号中的拐点。由于复杂地物背景噪声值的干扰会产生许多的“伪”拐点，需要予以删减。首先，将处于振幅阈值以下的拐点全部删除；其次，将检测出来的拐点位于波峰位置左侧称为左拐点，位于右侧的称为右拐点，通过判断波形分量峰值位置，两侧拐点跟其后一个波形点的单调性，将左、右拐点之间的异常点进行剔除，当存在多个左、右拐点时，取其平均值作为最终的左、右有效拐点值。由于噪声和波形叠加的影响，左右拐点经常不对称，通常取左、右拐点分别到波峰位置处的距离较小值作为半宽值。

本文通过计算回波局部最大值来确定高斯分量的峰值和中心位置[14]，采用高斯拐点匹配法求取左、右有效拐点，取左、右拐点到分量中心位置的较小值为半宽值，渐进剥离波形直到剩余波形低于振幅阈值线。具体步骤如下。

1)计算滤波后回波的局部最大值，即为第1个高斯分量的振幅值A1，其对应的横坐标值为高斯分量的中心位置T1。

2)根据高斯拐点匹配法，判断波峰两侧的左、右有效拐点Tm-1和Tm+1，半宽值σ取|Tm-1-T|和|Tm+1-T|的较小值，如式(4)所示。

σ=min(|Tm-1-T|，|Tm+1-T|)

(4)

3)根据3个基本高斯参数即可确定第1个高斯分量，并从波形中剥离出第1个高斯分量，将剩余波形重复步骤1)、步骤2)得到第2个高斯分量，直到剥离的剩余波形局部最大值小于设定的振幅阈值。

4)一般认为，一个复杂地物的回波信号最多分离出6个高斯分量，当探测的高斯分量多余6个时，根据脉冲和面积合并原则进行高斯分量合并[15]，一般将半宽小于发射脉宽的一半或者面积较小的分量合并到其邻近面积最大的高斯分量上，合并后高斯分量振幅取其中较大的高斯分量振幅值，中心位置和半宽值皆取二者合并前的平均值。

5)将剥离出的高斯分量进行整体拟合，求取拟合后的波形与原始波形之间的均方根误差，判断其是否满足拟合精度要求。

2.3 子波形的整体拟合

当激光足印内存在坡地或建筑物等具有复杂地物类型时，其回波波形通常较为复杂且伴有波形展宽现象，以及采用高斯滤波去噪会使各子波形振幅变小，因此需要对分解的子波形进行整体优化拟合，以使拟合波形最大程度逼近原始回波信号。最小二乘拟合法主要是通过最小化误差的平方和找寻数据的最佳函数匹配，通过判断原始数据和拟合结果的差值来调节迭代系数的大小，以达到最优的拟合结果。其在一定程度上具有弹性拟合的优势，被广泛应用于非线性方程的求解过程中。

本文采取最小二乘拟合算法对高斯分量参数进行整体拟合优化工作，使优化后的拟合波形尽可能逼近原始回波。首先将叠加式(2)展开如式(5)所示。

(5)

式中：An、Tn和σn分别代表第n个高斯分量的振幅值、峰值中心位置和半宽值。

对其进行全微分，表达式如式(6)所示。

f(x)=f(x0)+a1·dA1+b1·dT1+c1·dσ1+

a2·dA2+b2·dT2+c2·dσ2+…

(6)

误差方程如式(7)所示。

V=AX-L

(7)

A=[a1b1c1a2b2c2…]

(8)

式中：[anbncn]分别为回波函数f(x)对各子波的振幅A、中心位置T和半宽σ的一阶导数值。

L=[f(x)-f(x0)]

(9)

式中：f(x0)为子波形拟合后的回波信号。

求解误差方程如式(10)所示。

X=(ATA)-1ATL

(10)

在拟合前后差值L达到最小的情况下，通过解算式(10)可获得各高斯分量的最佳拟合参数。渐进剥离与整体拟合方法的流程如图3所示。

图3 渐进剥离与整体拟合流程图

3 实验结果与分析

3.1 实验数据

GLA01全球测高数据产品由美国冰雪数据中心(National Snow and Ice Data Center，NSIDC)发布，记录回波样本中包括激光足印的位置、脉冲强度及接收时间等信息。其1 s内的40个激光脚点共享一个地理坐标，地表每个激光足印的直径为72 m，相邻足印中心间距170 m[16]，从2017年8月份开始，数据产品由二进制格式改为HDF5格式进行分发。本文实验采用的是1A级GLA01数据产品，其主要属性如表1所示。

表1 实验数据主要属性

3.2 拐点提取实验

本文拐点提取实验以子号为20的波形数据为例，利用平面曲线离散点集拐点的快速查找算法对去噪后波形进行拐点提取，然后根据高斯拐点匹配法的相关原则，将低于振幅阈值线以及不符合左、右拐点条件的“伪”拐点予以删减，具体方法步骤参考上文。提取的整个波形拐点和经删减后获得的波形有效拐点分别如图4、图5所示，图中红色星号代表提取的拐点位置，蓝线代表去噪后波形，红线是根据背景噪声值设置的振幅阈值线。

图4 波形拐点图

图5 波形有效拐点图

3.3 子波形的渐进剥离实验

首先，进行波形峰值最大值检测提取第1个高斯分量，将波峰最大值及对应的波形中心位置作为分解的第1个高斯分量的振幅和中心位置，半宽则依据高斯拐点匹配法确定的左、右拐点来确定，分解的第1个高斯分量如图6(a)所示。接着，从完整去噪波形中剥离出第1个高斯分量，将剩余波形进行新一轮的波峰最大值检测和高斯拐点匹配，以提取第2个高斯分量，在图6(a)的基础上剥离的第2个高斯分量图如图6(b)所示。依此类推，渐进剥离出的高斯分量分别如图6(c)、图6(d)、图6(e)所示，直至剩余波形的最大振幅值低于设定的振幅阈值。图6(f)是渐进剥离分解后的剩余波形全部位于振幅阈值线以下的结果。图6(a)至图6(f)中的横坐标和纵坐标分别代表回波数据的相对时间和电压幅值，蓝线代表原始去噪后波形以及渐进剥离后的剩余波形，绿线代表分解出的高斯分量，红线代表振幅阈值线。

图6 子波形渐进剥离分解图

3.4 子波形的整体拟合实验

子波形的整体拟合实验选取主峰在中间、左侧以及右侧3个典型位置的复杂激光测高波形数据，分别命名为波形1、波形2和波形3，其索引子号分别为20、15和22。复杂原始回波通过波形分解可获得若干个子波形，将各子波形的3个高斯基本参数进行最小二乘拟合，获取精化后的高斯参数，以尽可能使拟合后波形最大限度地逼近原始回波信号，从而达到整体拟合的目的。图7为主峰在中间位置的波形1经本文方法分解后，获得的各个高斯分量及拟合波形与原始去噪后回波信号的对比情况。图8和图9分别为主峰位于左侧的波形2和位于右侧的波形3经渐进剥离分解并拟合后的结果。图7至图9中的蓝线代表原始去噪后波形，绿线代表经渐进剥离后获得的高斯分量，红线代表整体拟合后波形。

图7 波形1高斯分量拟合图

图8 波形2高斯分量拟合图

图9 波形3高斯分量拟合图

3.5 精度对比

为了验证本文方法在ICESat-GLAS回波信号分解的可靠性和适用性，通过高斯分量数目以及拟合波形与去噪回波间的均方根误差、相关系数和拟合优度等评判指标进行定量综合对比。均方根误差指的是实际观测值与拟合值之间的偏差，常用来分析数据的可靠性[17]，其计算如式(11)所示。

(11)

式中：N为波形的采样点个数；fest(x)为拟合波形；f(x)为原始去噪后波形。

相关系数反映的是拟合波形和原始去噪后回波之间的密切程度，取值位于[0，1]之间，越逼近1代表二者相关性越密切[18]，其计算如式(12)所示。

(12)

式中：Cov()和Var[]分别指f(x)和fest(x)的协方差和方差计算。

拟合优度也称为决定系数，其值是相关系数的平方，拟合优度值越接近1，代表拟合效果越好。

文中选取3种典型的激光测高数据波形1、波形2和波形3为研究对象进行测试，以经典的奇偶高斯拐点法、峰值检测法、波峰自动识别的高斯分解法作对比，以验证本文方法的可靠性和实用性。其中，经典的奇偶高斯拐点法由相邻的奇偶拐点来确定高斯参数；峰值检测法是通过查找波形中比左右相邻元素值都大的峰值点以确定其位置和幅值，再人为定义半宽取值区间为[3，6]；波峰自动识别法是结合上述2种方法的优势，分解速度提高且精度理论上更为可靠，各方法得到的结果统计如表2所示。同时以波形1为例，通过计算提取的高斯参数初值与精化后参数值间的误差，定量对比分析本文方法分解子波形参数初值的准确度，其中，ΔA、ΔT和Δσ分别代表分量振幅、中心位置和半宽精化前后的差值，具体结果如表3所示。

表2 不同波形和分解方法的结果统计

表3 高斯分量参数初值与精化值对比

3.6 实验分析

提取数目上，以波形1为例，奇偶高斯拐点法得到的高斯分量个数最多，达到8个，可见该方法受噪声影响较大，因噪声产生的“伪”拐点没有得到有效去除；峰值检测法和波峰自动识别法确定的高斯分量相对较少，均为4个，“伪”拐点问题得到了有效解决，但忽略了波形首尾的部分弱子波形；本文方法确定的高斯分量个数为5个，包括渐进剥离分解出4个较为明显的子波形和1个隐含的弱子波形，剩余波形的噪声属性明显，表明分解结果的合理性与可靠性。

提取精度上，与其他3种方法相比，均方根误差大小分别下降75%、66.85%和64.1%；相关系数分别提升11.16%、1.53%和0.81%；拟合优度分别提升23.54%、3.12%和1.64%。由此可见，通过本文方法得到的拟合波形和原始去噪后回波信号偏差较小，二者间的相关性也较为密切。

通过本文方法在典型位置波形1中提取的高斯参数初值与精化后值对比，各子波形电压幅值初值与精化值的最大差值为0.186 V，最小差值为0.010 V；峰值中心位置误差最大差值为1.514 ns，最小差值为0.016 ns；半宽误差最大差值为1.479 ns，最小差值为0.366 ns。可以看出，估计的高斯分量参数初值与精化值误差较小，表明本文方法能较为准确地测定高斯分量参数初值，同时经过参数精化和整体拟合后的波形也能近似原始回波信号，从而达到在总波形的约束下精确分解高斯分量的目的。

4 结束语

每个波形分量都对应地物的一次反射特征，波形分解的效果直接影响反演地物特征参数的精度。本文在分析经典波形分解的基本原理以及改进方法的基础上，利用渐进剥离波形与整体拟合的方法对ICESat-GLAS原始全波形信号进行波形分解，通过渐进剥离波形寻找局部最大值来确定分量峰值和位置，并结合高斯拐点匹配法左、右拐点原则来确定半宽值，最后对3种典型位置的波形利用4种分解方法实验并进行精度定量对比。实验结果表明，本文方法能准确测定分量参数初值和部分弱波，有效避免了传统分解算法中因噪声拐点产生的错误分量过多以及不合理给定半宽值等问题，同时对高斯分量参数初值进行精化和整体拟合，使分解的高斯分量受到总波形的约束，从而提高波形分解的精度。

此外，由于次波对主波的叠加效应影响,导致拟合波形主峰峰值变大,以及噪声去除不完全都会对高斯分量提取和拟合的精度产生一定影响。因此，接下来的研究工作将继续深入探讨次波对主波的影响因素，如采用最优的拟牛顿算法拟合子波形来降低主波峰值偏大的问题，以及结合更高效的波形去噪方法，如用改进的小波变换去噪方法,来降低因高斯滤波去噪引起的各子波振幅变小误差。