CN114488168A - 基于最大正向偏差的卫星激光测距全波形高斯拟合方法 - Google Patents

Abstract

本发明涉及一种基于最大正向偏差的卫星激光测距全波形高斯拟合方法，属于卫星激光测距领域。该方法的基本策略如下：初始状态拟合波形仅包含位于回波振幅最大值处的一个高斯分量，然后通过多次循环迭代，渐次地在拟合波形与实际回波的最大正向偏差处增加高斯分量，逐步缩小拟合波形与实际回波之间的偏差，直至满足拟合结束条件。本发明使用的波形拟合方法即适用于常规回波，也适用于复杂多峰、多个单脉冲合成类似单波的叠加波、面积较小的地物目标形成弱回波等情形，适用性更广。

Description

基于最大正向偏差的卫星激光测距全波形高斯拟合方法

所属领域：

本发明涉及一种基于最大正向偏差的卫星激光测距全波形高斯拟合方法，属于卫星激光测距领域。

背景技术：

激光测高是一种能够精确获取地面三维信息的主动遥感手段，卫星激光测高具备主动获取全球地表及目标三维信息的能力，不受白天黑夜的限制，能为快速获取包括境外地区在内的三维控制点以及立体测图提供服务，同时在极地冰盖测量、植被高度及生物量估测、云高测量、海面高度测量以及全球气候监测等方面都可以发挥重要作用。星载激光测高系统主要由激光测距仪、精确轨道测量系统和精密姿态测量系统组成，激光测距仪可以实现系统与测量目标之间的高精度距离测量，是星载激光测高系统的核心所在。

由于测距仪激光光束具有一定的发散角，并且大气对光束具有散射作用，因此激光光束在地面上照亮的足印是一个有一定面积的光斑，光斑直径约几十米甚至上百米。全波形激光测距仪能够记录足印范围内不同目标面回波脉冲的完整波形信息，回波脉冲是对足印内各点反射信号按时间先后顺序的记录，可以认为全波形回波脉冲是在一个激光足印中，不同距离的多个目标面散射后产生的单回波脉冲的叠加。通过分解全波形回波，能够获取对应不同目标面的单脉冲，进而由分析这些单脉冲的参数得到相应目标面的距离与其他特征信息，反映激光足印中地物的垂直分布特征。全波形激光测距仪可用于地表高程测量、植被特征参数反演、地表生物量估计、地表分类、广义高程点库建设等。

星载激光测距全波形回波中存在大量的由多脉冲合成的波形数据，全波形数据的分解是激光测距数据应用的关键步骤。目前常见的波形分解主要有三种方法：反卷积法、小波分析法和高斯分解法。由于反卷积运算本身存在病态问题且实现复杂，故其实际应用并不广泛；小波变换对于信号检测异常敏感，在叠加波检测中普遍会出现过检测现象，使用小波分析法对激光测距全波形回波进行波形分解的研究也相对较少；高斯分解法将激光测距仪接收的经地表反射后的回波波形看作一个或多个高斯脉冲的叠加对波形进行拟合，由于其模型相对简单且具有较好的拟合效果，分解所得高斯分量的个数和参数可反映足印范围内的多种特征，是目前波形分解中应用最为广泛的方法。然而，星载激光测距全波形信号信噪比低，且受地表坡度和复杂度的影响，极易出现复杂多峰、多个单脉冲合成类似单波的叠加波、面积较小的地物目标形成弱回波等情形；现有的高斯分解方法通常是根据回波相邻采样点二阶导数异号判断拐点，再根据拐点检测高斯分量，在不具有典型高斯曲线特征的回波处，非常容易漏掉某些地物目标对应的高斯分量或检测到错误的高斯分量，得到精度较低的拟合结果，导致提取的地物参数与实际情况并不相符。

发明内容：

针对以上问题，本发明提供一种基于最大正向偏差的卫星激光测距全波形高斯拟合方法，渐次地在拟合波形与实际回波的最大正向偏差处增加高斯分量，能在很大程度上减少高斯分量缺失和误检现象的发生。

本发明的技术方案是：提供一种基于最大正向偏差的星载大光斑激光测距全波形数据高斯拟合方法，在对复杂多峰回波进行波形拟合时往往容易出现高斯分量缺失的问题，若规定实际回波振幅值大于拟合回波振幅值时为正向振幅偏差，则理论上缺失的高斯分量均值应与最大正向偏差所在位置接近，即在最大正向偏差处新增合适的高斯分量可以最大限度地使拟合波形逼近实际回波。基于以上思想，本发明提出如下的全波形回波拟合策略：初始状态拟合波形仅包含位于回波振幅最大值处的一个高斯分量，然后通过多次循环迭代，渐次地在拟合波形与实际回波的最大正向偏差处增加高斯分量，逐步缩小拟合波形与实际回波之间的偏差，直至满足拟合结束条件。

星载激光测距全波形回波脉冲可以认为是在一个激光足印范围内不同距离的多个目标面散射后产生的单回波脉冲的叠加。激光发射波形一般近似为基模高斯形式，因此卫星接收的经地表反射后的回波波形可近似看作一个或多个高斯脉冲的叠加。Brenner提出使用多个高斯函数组合的方法对回波波形进行拟合，公式表达为

(其中：w(t)表示t时刻回波振幅；为估计的背景噪声均值；A_m、t_m、分别为第m个高斯函数的振幅值、均值和标准差)，本发明采用Brenner提出的多个高斯函数叠加公式来表达全波形回波脉冲。

本发明的技术方案包括以下步骤：

(1)对输入的原始回波进行预处理，包括：背景噪声估计得到背景噪声均值和标准差，使用一维高斯滤波器进行平滑滤波；

(2)回波有效信号检测，根据步骤(1)中估计得到的背景噪声均值和标准差，确定回波有效信号振幅阈值最下值Th_sig，分别从原始回波两端开始搜索第一个振幅值大于Th_sig的点作为有效信号起始位置sig和结束位置sig，将sig和sig之间的波形作为回波有效信号，后续处理若无说明，均是对sig和sig之间的波形进行操作；

(3)背景噪声重新估计，从平滑回波中取位于sig左侧的若干个点，计算这些点的振幅均值μ_no和标准差σ_no，作为高斯拟合波形的噪声估计参数；

(4)初始高斯分量检测及参数初值估计，搜索平滑回波振幅最大值点P_max，记点P_max对应采样点的索引序号为I_max，在采样点I_max处新增高斯分量，并估计该高斯分量的初始参数，此时拟合波形中有且仅有一个高斯分量；

(5)拟合波形参数调整优化，使用非线性最小二乘方法对拟合波形做参数调整，若调整后拟合残差变小，使用调整后参数进入后续处理，否则使用调整前参数；

(6)检测当前拟合波形中是否包含无效高斯分量，若有则剔除；

(7)计算拟合残差，判断是否满足拟合结束条件，结束拟合条件为：拟合次数大于等于某一设定值Iter_max或拟合残差变化量与拟合残差的比值小于0.01，对于GLAS数据Iter_max取值为12，若满足则结束拟合.方法结束；若不满足则进入下一步；

(8)计算拟合波形与平滑回波之间的振幅偏差；

(9)在最大正向偏差处，新增高斯分量并估计其初始参数；

(10)使用非线性最小二乘方法对拟合波形做参数调整，重新进入步骤(5)。

本发明的有益效果如下：

(1)使用拐点检测高斯分量对噪声非常敏感，波形的轻微抖动就极易造成高斯分量的遗漏和误检，本发明使用拟合波形与实际回波的最大正向偏差来确定新增高斯分量的位置，能很好地避免这一问题的发生；

(2)本发明首先在回波振幅最大值处增加一个高斯分量作为初始参数，可以避免激光足印范围内主要目标面对应高斯分量的遗漏；

(3)本发明通过多次循环迭代，渐次地在拟合波形与实际回波的最大正向偏差处增加高斯分量，逐步缩小拟合波形与实际回波之间的偏差，设置合理的拟合结束条件，可以最大限度地保留足印范围内各目标面对应的高斯分量，使提取的地物参数与实际情况尽量相符；

本发明使用的波形拟合方法即适用于常规回波，也适用于复杂多峰、多个单脉冲合成类似单波的叠加波、面积较小的地物目标形成弱回波等情形，适用性更广。

附图说明：

图1是实施例中基于最大正向偏差的波形拟合流程；

图2是实施例中标识为341405549-22的记录原始波形与平滑波形；

图3是实施例中标识为341405549-22的记录NSIDC拟合波形与本发明算法拟合波形；

图4是实施例中标识为341405549-22的记录本发明拟合波形与分解波形；

图5是实施例中标识为341405589-2的记录原始波形与平滑波形；

图6是实施例中标识为341405589-2的记录NSIDC拟合波形与本发明算法拟合波形；

图7是实施例中标识为341405589-2的记录本发明拟合波形与分解波形。

具体实施方式：

美国于2003年发射的ICESat对地观测卫星上，搭载了具备全波形记录功能的GLAS激光测高仪，用于两极冰盖监测、全球森林生物量估算、陆地高程测量等。GLAS在国内外一直备受关注，有不少学者投入了大量精力开展该数据的处理、分析与应用研究工作。GLAS数据是非常典型的星载大光斑激光测距全波形数据，因此，本发明以该数据为例进行波形拟合处理。

1)原始回波预处理，包括背景噪声估计和使用一维高斯滤波器进行平滑滤波

(1)背景噪声估计

回波信号的采集过程受到大气、云和系统自身等因素的影响，接收到的回波信号中常常混有背景噪声。回波波形前半部分对应的是脉冲还未到达地面的时间段，可以视为只包含背景噪声而不包含有效信号，故可取回波前半部分一定数量的点，计算其均值和标准差来估计背景噪声。

背景噪声估计的具体实施方式如下：

a)取回波前n(一般n＝100)个点作为背景噪声统计的初始点集P₀；

b)统计P₀中n个点的均值μ₀与标准差σ₀；

c)定义回波有效信号振幅阈值Th_{sig_min0}＝μ₀+N_sig*σ₀(N_sig为常数系数，一般取4.5)，由回波左侧开始搜索第一个振幅值大于Th_{sig_min0}的点作为有效信号起始位置sig_beg0；

d)由(sig_beg0-50)处开始向左侧取n个点作为统计背景噪声的最终点集P；

e)统计P中n个点的均值和标准差，即为背景噪声的均值和标准差，分别记为μ_{noise_ob}和σ_{noise_ob}。

(2)平滑滤波

激光回波的噪声可以看作非零均值高斯白噪声，即直流信号和零均值白噪声的叠加。正确估计噪声均值并作归零处理，使用宽度适宜的高斯滤波器即可去除高斯白噪声。高斯滤波是一种低通滤波器，一维高斯滤波的原理是将每个所述有效原始回波波形和高斯函数进行卷积处理，抑制高频的背景噪声，实现对所述有效原始回波波形的平滑，从而得到高斯滤波后的回波波形。

一维高斯滤波器的宽度选择为星载激光测距仪发射脉冲的脉冲宽度(也即脉冲的半高全宽，FWHM)，一维高斯滤波器窗口大小W＝2×ceil(3σ)(其中，ceil为函数，表示取大于或者等于x的最小整数)，一维高斯滤波器的高斯卷积核计算公式为

(i＝1，2，…，(其中，c为高斯模板中心位置，计算公式为

floor为函数，表示取小于或者等于x的最小整数)。

将高斯卷积核中的值进行归一化处理，计算公式如下：

确定了一维高斯滤波器的宽度、窗口大小以及高斯卷积核后，对原始回波波形中的采样点进行一维高斯滤波的计算公式为：

上式中，为滤波前的原始回波波形的第i个采样点的振幅值，Y为滤波后的回波波形第i个采样点的振幅值，N为原始回波波形中采样点的总个数，L＝floor。

2)回波有效信号检测

激光测距仪记录的全波形回波中，只有一部分是地面目标的后向散射回波形成的有效信号，还有一大部分回波记录的噪声信号，波形拟合时应该只对有效信号区域的回波进行，因此需要检测出有效信号的开始和结束位置。

定义回波有效信号振幅阈值Th_{sig_min}＝μ_{noise_ob}+N_sig*σ_nois(为常数系数，一般取4.5)，从原始回波左侧开始搜索第一个振幅值大于Th_sig的点作为有效信号起始位置sig，从原始回波右侧开始搜索第一个振幅值大于Th_sig的点作为有效信号结束位置sig，后续处理若无说明，均是对sig和sig之间的波形进行操作；

3)背景噪声重新估计

由于波形拟合是对平滑滤波后的回波进行的，如果拟合时使用原始回波进行背景噪声估计得到的均值μ_nois和方差σ_nois会引入一定的误差，因此，需要对平滑滤波后的回波重新进行背景噪声估计。对平滑滤波后回波，由索引(sig_beg-处开始向左侧取n个点(n一般等于100)，计算这些点的振幅均值μ_no和标准差σ_no，作为高斯拟合时的背景噪声估计值。

4)初始高斯分量检测及参数初值估计

本发明首先在平滑回波振幅最大值处增加一个高斯分量作为初始的拟合回波，检测初始高斯分量及参数初值估计的具体实施方式如下：

(1)搜索平滑回波振幅最大值点P，记点P对应采样点的索引序号为I，振幅值为A；

(2)初始高斯分量的振幅值估计值A₁＝A，均值估计值t₁＝I；

(3)由处开始分别向左右两侧搜索振幅值依次递减的边界采样点，其索引序号分别记为Y和I_r；

(4)定义常数系数、、、，分别表示与高斯分量均值相距、、、处的采样点振幅值与均值处采样点振幅值的比值，则有各系数值分别为r_a0.135335、r_a2＝0.606、r_a3＝0.882、r_a4＝0.969；

(5)由索引I处开始向左搜索第一个振幅值小于(r_a1·的点，直至索引号小于Y，若搜索成功则记该点的索引序号为，左侧高斯分量标准差估值为

继续下一步处理，若失败，则依次搜索第一个振幅值小于(r_a2·、(r_a3·、(r_a4·的点直至成功，若对应点的索引序号分别记为、、，对应的左侧高斯分量标准差估值计算公式分别为

若全部失败，则计算失败；

(6)由索引I处开始向右搜索第一个振幅值小于(r_a1·的点，直至索引号大于I_r，若搜索成功则记该点的索引序号为，右侧高斯分量标准差估值为

继续下一步处理，若失败，则依次搜索第一个振幅值小于(r_a2·、(r_a3·、(r_a4·的点直至成功，若对应点的索引序号分别记为、、，对应的右侧高斯分量标准差估值计算公式分别为

若全部失败，则计算失败；

(7)若和均计算成功，则高斯分量的标准差估计值取二者中的较小值，若只有一个计算成功，则取计算成功的值，若二者均计算失败，表示波形拟合失败，返回不再进行后续处理。

5)拟合波形参数调整优化

由Gauss-Newton(高斯-牛顿)方法演进而来的LM(Levenberg-Marquardt，列文伯格-马夸尔特)算法是非线性优化的“标准”方法，被广泛用于激光雷达回波数据拟合优化。通过引入阻尼系数，将Gauss-Newton法与最速下降法结合起来，在很小时，算法趋于Gauss-Newton法，拥有Gauss-Newton法的快速收敛性，而当在很大时，算法趋于最速下降法，拥有最速下降法的强适应性，使LM法可以更加快速有效的完成拟合优化。本发明使用LM方法来进行拟合参数的调整，具体实施方式如下：

(1)记平滑回波为

拟合波形为

N为参与拟合的波形采样点数目；

(2)记待拟合目标参数为

m为高斯分量个数；

(3)构建误差函数

(i＝1，2，…，构建残差函数

(4)设置阻尼系数的初始值(一般取0.01)，设置阻尼系数的变化系数的初始值(一般取10)，设置最大迭代次数Iter_max(处理GLAS数据时设为12)，目标参数的初始值为各高斯分量的初始参数估计值；

(5)对目标参数的第i次估值，计算误差函数的雅克比矩阵J，计算公式如下：

其中，

为m个高斯函数第i次的估值；

(6)计算误差函数f(的值以及残差F(；

(7)计算第i次拟合的目标参数增量，计算公式如下：

Δ_i＝(J^TJ+λI)^-1J^T[f₁，f₂，…，f_N]^T

(8)计算更新后的目标参数

(9)对更新后的目标参数

删除或合并不满足条件的分量得到

(10)计算目标参数更新后的残差

(11)若

则更新待拟合目标参数，并按下式更新阻尼系数；

(12)若

则不更新待拟合目标参数，并按下式更新阻尼系数。

λ′＝λβ

6)无效高斯分量剔除

在使用LM方法进行拟合参数的调整过程中，可能会引入无效的高斯分量，需要在参数调整后进行剔除，遍历当前拟合波形中的高斯分量，判断是否包含无效高斯分量并剔除的具体实施方式如下：

(1)若高斯分量振幅值A_i＜Th_sig，则当前高斯分量为无效高斯分量，将其提出并返回，否则继续判断；

(2)若高斯分量均值小于最左侧高斯分量的均值或大于最右侧高斯分量的均值，则当前高斯分量为无效高斯分量，将其提出并返回，否则继续判断；

(3)若高斯分量标准差σ₁＜Th_σ或σ₁＞Th_σ(对于GLAS数据，取Th_{σ_min}＝，Th_{σ_max}＝min((Sig_end-Sig_beg)/4，2)，则当前高斯分量为无效高斯分量，将其提出并返回，否则继续判断；

(4)计算当前高斯分量与左右相邻高斯分量均值距离，分别记为Δt和Δt，若Δt_left＜Th_{peak_inte}(对于GLAS数据，Th_{peak_interval}＝1)，则将当前高斯分量与左侧高斯分量合并，若Δ_tright＜Th_{peak_inte}，则将当前高斯分量与右侧高斯分量合并，否则结束判断并返回，高斯分量合并方法如下(以合并当前高斯分量与左侧相邻高斯分量为例)：

a)计算左侧高斯分量权重，W_left＝S_left/(S_left+S(其中，S和分别表示左侧高斯分量和当前高斯分量对应高斯曲线与横轴合围部分面积)；

b)计算当前高斯分量权重，W_cur＝1.0-W；

c)计算合并后高斯分量振幅值，A_merge＝max(A_left，A；

d)计算合并后高斯分量均值，t_merge＝W_left*t_left+W_cur*；

e)计算合并后高斯分量标准差，σ_merge＝W_left*σ_left+W_cur*。

7)计算拟合残差并判断是否结束拟合

使用当前拟合参数计算拟合残差

判断是否满足拟合结束条件(结束拟合的条件为拟合次数大于等于Iter_max或拟合残差变化量与拟合残差的比值小于0.01，对于GLAS数据Iter_max取值为12)，若满足则结束拟合，若不满足则进入下一步处理。8)计算拟合波形与平滑回波之间的振幅偏差

使用误差函数f计算当前拟合回波与平滑回波之间的振幅偏差，计算公式为

9)在最大正向偏差处新增高斯分量并估计其初始参数

本发明通过多次循环迭代，渐次地在拟合波形与实际回波的最大正向偏差处增加高斯分量，逐步缩小拟合波形与实际回波之间的偏差，直至满足拟合结束条件，以此来最大限度地保留足印范围内各目标面对应的高斯分量，使提取的地物参数与实际情况尽量相符。

在最大正向偏差处新增高斯分量并估计其初始参数的具体实施方式如下：

(1)记原有拟合目标参数对应的拟合残差为I；

(2)搜索当前拟合回波与平滑回波之间的最大正向偏差，记其索引序号为I_dev；

(3)向拟合目标参数中新增一个高斯分量，新增高斯分量的振幅值估计值等于I_dev处采样点的振幅值，均值估计值t_j＝I_dev，标准差估计值参照“4)初始高斯分量检测及参数初值估计”中的步骤(3)～(7)确定。

10)使用非线性最小二乘方法对拟合波形做参数调整

(1)将中所有高斯分量参数重置为初始估计参数，计算此时的拟合残差记为R_new，使用阻尼系数初始值按“5)拟合波形参数调整优化”中的步骤(1)～(8)进行一次拟合目标参数调整，计算参数调整后的拟合残差记为R_new；

(2)若R_{new_cur}＜R_new，则使用调整后的拟合目标参数，令阻尼系数λ＝λ，并转至“5)拟合波形参数调整优化”中的步骤(1)继续处理，否则转至本部分的步骤(4)；

(3)若R_{new_prev}＜I，则使用调整后的拟合目标参数，令阻尼系数λ＝，并转至“5)拟合波形参数调整优化”中的步骤(1)继续处理，否则转至本部分的步骤(4)；

(4)使用新增高斯分量前的原有拟合目标参数，令λ＝λ，转至“5)拟合波形参数调整优化”中的步骤(1)继续处理。

11)输出拟合结果

循环迭代直至满足拟合结束条件(拟合次数大于等于Iter_max或拟合残差变化量与拟合残差的比值小于0.01，对于GLAS数据Iter_max取值为12)，若拟合成功(拟合残差变化量与拟合残差的比值小于0.01)，输出最终的拟合结果及拟合残差，否则(拟合次数大于等于Iter_max)返回拟合失败消息。

为了验证本发明所提出的基于最大正向偏差的星载大光斑激光测距全波形数据高斯拟合方法，以GLAS官方公布的处理结果为参考，对多个单脉冲合成类似单波的叠加波、复杂多峰、弱回波三种较难处理的回波波形，选择对应记录与本发明方法的拟合结果进行对比，对比结果如图2～图7所示，其中NSIDC(NationalSnowandIceData Center)拟合波形代表GLAS公布的官方方法处理结果。

图2～图4是标识为341405549-22的记录处理结果，代表多个单脉冲合成类似单波的叠加波的情形。NSIDC拟合结果中只包含一个回波，拟合残差为0.022050；本发明方法拟合结果为四个回波，拟合残差为0.005658。通过图3中的波形曲线对比可发现，NSIDC拟合波形因为漏掉了强度较弱的几个回波，导致与原始波形有较大的偏差。

图5～图7是标识为341405589-2的记录处理结果，代表了复杂多峰和弱回波的情形，NSIDC拟合结果中只包含一个回波，拟合残差为0.043140；本发明方法拟合结果为三个回波，拟合残差为0.012308。通过对比可发现，NSIDC只检测出了强度最大的高斯分量，中间和最右侧的高斯分量应该是由于拐点检测失误而丢失，在拟合过程中导致强度最大的高斯分量均值位置向右侧偏移了约2ns，且由于丢失高斯分量，使得拟合结果不能真实地反映激光足印中地物的垂直分布特征，提取的地物参数也与实际情况并不相符。

以上结果表明，本发明的方法能很好地避免高斯分量遗漏和误检的问题，即适用于常规回波，也适用于复杂多峰、多个单脉冲合成类似单波的叠加波、面积较小的地物目标形成弱回波等情形，适用性广泛。

Claims (1)

1.基于最大正向偏差的卫星激光测距全波形高斯拟合方法，其特征在于，包括以下步骤：

(1)对输入的原始回波进行预处理，包括：背景噪声估计得到背景噪声均值和标准差，使用一维高斯滤波器进行平滑滤波；

(2)回波有效信号检测，根据步骤(1)中估计得到的背景噪声均值和标准差，确定回波有效信号振幅阈值最下值Th_{sig_}，分别从原始回波两端开始搜索第一个振幅值大于Th_{sig_}的点作为有效信号起始位置sig和结束位置sig，将sig和sig之间的波形作为回波有效信号；

(3)背景噪声重新估计，从平滑回波中取位于sig左侧的若干个点，计算这些点的振幅均值μ_no和标准差σ_noi，作为高斯拟合波形的噪声估计参数；

(4)初始高斯分量检测及参数初值估计，搜索平滑回波振幅最大值点P_max，记点P_max对应采样点的索引序号为I_max，在采样点I_max处新增高斯分量，并估计该高斯分量的初始参数，此时拟合波形中有且仅有一个高斯分量；

(5)拟合波形参数调整优化，使用非线性最小二乘方法对拟合波形做参数调整，若调整后拟合残差变小，使用调整后参数进入后续处理，否则使用调整前参数；

(6)检测当前拟合波形中是否包含无效高斯分量，若有则剔除；

(7)计算拟合残差，判断是否满足拟合结束条件，结束拟合条件为：拟合次数大于等于某一设定值Iter_max或拟合残差变化量与拟合残差的比值小于0.01，对于GLAS数据Iter_max取值为12，若满足则结束拟合，方法结束；若不满足则进入下一步；

(8)计算拟合波形与平滑回波之间的振幅偏差；

(9)在最大正向偏差处，新增高斯分量并估计其初始参数；

(10)使用非线性最小二乘方法对拟合波形做参数调整，重新进入步骤(5)。

Images