CN110673109B - 一种星载大光斑激光雷达全波形数据分解方法 - Google Patents

Abstract

本发明公开了一种星载大光斑激光雷达全波形数据分解方法，该方法包括以下步骤：步骤S1，对星载大光斑激光雷达全波形数据包括的每个原始回波波形数据进行预处理，得到滤波后的回波波形；步骤S2，检测所述滤波后的回波波形中的明显峰值点和与所述明显峰值点对应的拐点，并估计有效峰值点对应的高斯分量的参数的初始值；步骤S3，检测有效明显峰值点两侧的叠加高斯分量，并估计叠加高斯分量的参数初始值；步骤S4，对滤波后的回波波形中检测到的高斯分量的总个数进行更新并对所述高斯分量的参数进行最优化估计。本发明结合拐点匹配和迭代优化算法，改进高斯分量的初始检测策略，得到了较好的回波波形的高斯分量初始估计。

Description

一种星载大光斑激光雷达全波形数据分解方法

技术领域

本发明涉及星载激光测高数据处理技术领域，特别涉及一种星载大光斑激光雷达全波形数据分解方法。

背景技术

全波形激光雷达技术能够记录地物回波脉冲的完整波形信息，通过分析后向散射波形，可以得到包括激光脚点在内的地物的物理属性，反映地物的垂直分布特征。随着激光雷达技术的发展，星载全波形激光雷达技术成为一种有效的主动对地观测手段。1999年，NASA设计了一套机载大光斑激光测高仪系统LVIS(Laser Vegetation Imaging Sensor)，该系统能够数字化采样激光发射脉冲和回波脉冲，获取子树冠和冠层顶部地形数据。NASA在2003年发射了第一颗主要用于极地冰量测量的冰、云和陆地海拔测量卫星ICESat(Ice,Cloud,and land Elevation Satellite)，该卫星搭载了地球激光测高系统GLAS(Geoscience Laser Altimeter System)。GLAS能够记录被探测地物返回脉冲的全波形信息，可用于地表高程测量和地物高度提取等。

对于全波形激光雷达技术而言，全波形数据的处理是一个关键步骤，而全波形数据的处理方法主要有三种，即阈值法，反卷积法和波形分解法。相比较前两种方法，常采用基于高斯分解的波形分解方法对接收波形进行处理。

现有技术主要针对的是地面和机载激光雷达的全波形数据处理。对于星载的大光斑激光雷达全波形数据而言，由于其全波形信号信噪比低，且易受地表坡度和复杂度的影响。因此，星载激光雷达全波形数据的处理是一个关键问题。

对于星载的大光斑激光雷达全波形数据处理而言，目前主要以面向应用为主，在全波形数据的处理上，还未充分的考虑到不同的地表覆盖状况，且没有同时针对叠加波形、复杂多峰波形以及弱回波波形的分解难点的技术方案。

发明内容

针对上述问题，本发明提供了一种星载大光斑激光雷达全波形数据分解方法，解决不同地表覆盖类型条件下的星载激光雷达全波形数据的分解问题，特别是接收波形数据中的叠加波形、复杂多峰波形以及弱回波波形的检测和分解问题，以提高地物的波形特征参数的提取精度，进而提高卫星激光高程控制点的提取精度和地物参数反演的精度。由此可见，本发明的技术方案可充分地消除由于现有技术的限制和缺陷导致的一个或多个问题。

本发明提供了一种星载大光斑激光雷达全波形数据分解方法，其特征在于，所述方法包括以下步骤：

步骤S1，对星载大光斑激光雷达全波形数据包括的每个原始回波波形数据进行预处理，以判断每个所述原始回波波形是否为有效原始回波波形，并对所述有效原始回波波形进行一维高斯滤波处理，得到滤波后的回波波形；

步骤S2，检测所述滤波后的回波波形中的明显峰值点和与所述明显峰值点对应的拐点，基于检测到的明显峰值点和与所述明显峰值点对应的拐点，判断检测到的所述明显峰值点是否有效，并估计有效峰值点对应的高斯分量的参数的初始值；

步骤S3，检测所述有效明显峰值点两侧的叠加高斯分量，并估计所述叠加高斯分量的参数初始值；

步骤S4，对所述滤波后的回波波形中包含的高斯分量的个数进行更新并对所述高斯分量的参数进行最优化估计。

优选的，步骤S1具体包括以下子步骤：

步骤S1.1，估计所述星载大光斑激光雷达全波形数据中的每个所述原始回波波形的背景噪声水平；

步骤S1.2，基于所述背景噪声水平设置背景噪声阈值，并将所述原始回波波形的采样点的最大振幅值与所述背景噪声阈值进行比较，以判断所述原始回波波形是否为有效原始回波波形，并仅保留所述有效原始回波波形；

步骤S1.3，对所述有效原始回波波形进行一维高斯滤波处理，得到滤波后的回波波形。

步骤S1.4，计算滤波后的回波波形的采样点的最大振幅值，并将所述最大振幅值与该回波波形的背景噪声阈值进行比较，如果所述最大振幅值小于等于所述背景噪声阈值，则保留滤波前的原始有效回波波形，否则，保留滤波后的回波波形。

优选的，在步骤S1.1中，所述原始回波波形的背景噪声水平包括所述原始回波波形的背景噪声的平均值μ_n和背景噪声的标准差δ_n，其中：

背景噪声的平均值μ_n的计算公式如下式所示：

其中，N为选取的某个原始回波波形中采样点数的总和，n(优选的，n取20)为在所述原始回波波形的采样点中开头部分的和结尾部分分别选取的采样点数，V_i是第i个采样点的振幅。

背景噪声的标准差δ_n的计算公式如下式所示：

其中，N为选取的某个原始回波波形中采样点数的总和，n(优选的，n取20)为在所述原始回波波形的采样点中开头部分的和结尾部分分别选取的采样点数，V_i是第i个采样点的振幅。

优选的，在步骤S1.1中，根据以下公式设置背景噪声的阈值(th)：

th＝μ_n+m·δ_n

其中，th为背景噪声阈值，μ_n为背景噪声的平均值，δ_n为背景噪声的标准差，m为常数值(优选的，m取4.5)。

其中，如果原始回波波形中采样点的最大振幅值大于背景噪声阈值(th)，则该原始回波波形可被判定为有效原始回波波形，否则将该原始回波波形判定为噪声回波波形，并仅保留有效原始回波波形。

优选的，步骤S1.3具体包括以下子步骤：

步骤S1.3.1，确定一维高斯滤波器的宽度和一维高斯滤波器的窗口大小，其中，一维高斯滤波器的宽度(σ)选择为所述星载激光雷达的发射波形的半高宽(FWHM)；一维高斯滤波器窗口W大小如下式计算：

W＝1+2·ceil(3·σ)

其中，ceil(x)为函数，其功能是返回大于或者等于指定表达式的最小整数，W为一维高斯滤波器的窗口大小，σ为一维高斯滤波器的宽度；

步骤S1.3.2，计算一维高斯滤波器的高斯卷积核，具体计算如下式所示：

C＝floor(W/2)+1

其中，K(i)为一维高斯滤波器的高斯卷积核，W为一维高斯滤波器的窗口大小，C为高斯模板中心位置，σ为一维高斯滤波器的宽度，floor(x)为函数，其功能是取不大于x的最大整数。

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

其中，K(i)为一维高斯滤波器的高斯卷积核，W为一维高斯滤波器的窗口大小，C为高斯模板中心位置，σ为一维高斯滤波器的宽度。

步骤S1.3.3，基于确定的一维高斯滤波器的宽度和窗口大小以及高斯卷积核，对所述有效原始回波波形中的所有的采样点进行一维高斯滤波。

具体的，根据下述公式对有效原始回波波形中所有的采样点进行一维高斯滤波：

其中，y(i)为滤波前的有效原始回波波形的第i个采样点的振幅值，Y(i)为滤波后的回波波形的第i个采样点的振幅值，L＝floor(W/2)，W为一维高斯滤波器的窗口大小，K(L+j+1)为一维高斯滤波的高斯卷积核，N为有效原始回波波形中采样点的总个数。

优选的，步骤S2具体包括以下子步骤：

步骤S2.1，通过邻域窗口遍历所述滤波后的回波波形，来检测所述滤波后的回波波形中的明显峰值点；

步骤S2.2，计算所述滤波后的回波波形的采样点的二阶导数，并根据相邻的两个采样点的二阶导数是否异号的准则，检测所述明显峰值点左右两侧的拐点；

步骤S2.3，基于检测到的明显峰值点和与所述明显峰值点对应的拐点，判断所述明显峰值点是否有效；

步骤S2.4，估计确定的所述有效明显峰值点对应的高斯分量的参数的初始值。

优选的，步骤S2.3具体包括以下子步骤：

步骤S2.3.1，如果明显峰值点的左侧拐点和右侧拐点均不存在，则判定所述明显峰值点为无效峰值点；如果仅存在左侧拐点，则执行步骤S2.3.2-S2.3.4，如果仅存在右侧拐点，则执行步骤S2.3.5-S2.3.7，如果同时存在左侧拐点和右侧拐点，则执行步骤S2.3.8-S2.3.10；

步骤S2.3.2，计算所述明显峰值点的所有左侧拐点的时间坐标的平均值；

其中，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，m为明显峰值点的左侧拐点的总数，

为明显峰值点左侧第k个拐点的时间坐标值。

步骤S2.3.3，基于所述所有左侧拐点的时间坐标的平均值和所述明显峰值点的时间坐标值，计算所述左侧拐点到所述明显峰值点的时间距离；

d_L＝|μ^L-X(i)|

其中，d_L为左侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值。

步骤S2.3.4，判断所述左侧拐点到所述明显峰值点的时间距离是否大于等于星载激光雷达的发射波形的半高宽的一半，如果是，则确定所述明显峰值点为有效峰值点，否则为无效峰值点；

其中，d_L为左侧拐点到明显峰值点的时间距离，FWHM为所述星载激光雷达的发射波形的半高宽。

步骤S2.3.5，计算所述明显峰值点的所有右侧拐点的时间坐标的平均值；

其中，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值，n为明显峰值点的右侧拐点的总数，

为明显峰值点的右侧第k个拐点的时间坐标值。

步骤S2.3.6，基于所述所有右侧拐点的时间坐标的平均值和所述明显峰值点的时间坐标值，计算所述右侧拐点到所述明显峰值点的时间距离；

d_R＝|μ^R-X(i)|

其中，d_R为右侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值。

步骤S2.3.7，判断所述右侧拐点到所述明显峰值点的时间距离是否大于等于星载激光雷达的发射波形的半高宽的一半，如果是，则确定所述明显峰值点为有效峰值点，否则为无效峰值点；

其中，d_R为右侧拐点到明显峰值点的时间距离，FWHM为所述星载激光雷达的发射波形的半高宽。

步骤S2.3.8，分别计算所述明显峰值点的所有左侧拐点和所有右侧拐点的时间坐标的平均值，具体计算如下式：

其中，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值，m为明显峰值点的左侧拐点的总数，n为明显峰值点的右侧拐点的总数，

为明显峰值点左侧第k个拐点的时间坐标值，

为明显峰值点的右侧第k个拐点的时间坐标值。

步骤S2.3.9，分别计算左侧拐点和右侧拐点到明显峰值点的时间距离。具体的，根据下式计算左侧拐点和右侧拐点到明显峰值点的时间距离：

其中，d_L为左侧拐点到明显峰值点的时间距离，d_R为右侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值。

步骤S2.3.10，分别判断所述左侧拐点和右侧拐点到所述明显峰值点的时间距离是否大于等于激光雷达的发射波形的半高宽的一半，如果所述左侧拐点和右侧拐点中的至少一者到所述明显峰值点的时间距离大于等于激光雷达的发射波形的半高宽的一半，则确定所述明显峰值点为有效峰值点，否则为无效峰值点。

优选的，步骤S3具体包括以下子步骤：

步骤S3.1，设置区间[X(i)-D,X(i)+D]，其中，X(i)为所述有效明显峰值点对应的高斯分量对应的时间坐标值，D为估计的所述有效明显峰值点对应的高斯分量的标准差的初始值；

步骤S3.2，基于设置的所述区间，根据以下公式检测所述有效明显峰值点左侧的叠加高斯分量：

且y(j1)＜y(k1)

其中，j1、k1分别为在所述有效明显峰值点左侧单调区间的设置的所述区间之外检测到的两个相邻的拐点在X时间坐标数据中对应的位置，X(j1)和X(k1)分别为在所述有效明显峰值点左侧单调区间设置的所述区间之外检测到的两个相邻的拐点的时间坐标值，y(j1)和y(k1)分别为在所述有效明显峰值点左侧单调区域检测到的两个相邻的拐点在滤波前的有效原始回波波形中的振幅值，FWHM表示星载激光雷达的发射波形的半高宽；如果满足上式的条件，则在所述有效明显峰值点的左侧，在X时间坐标数据中对应序号分别为j1，k1的两个相邻的拐点区域检测到一个叠加的高斯分量；

步骤S3.3，基于设置的所述区间，根据以下公式检测所述有效明显峰值点右侧的叠加高斯分量：

且y(j2)＞y(k2) (26)

其中，j2、k2分别为在所述有效明显峰值点右侧单调区间的设置的所述区间之外检测到的两个相邻的拐点在X时间坐标数据中对应的位置，X(j2)和X(k2)分别为在所述有效明显峰值点右侧单调区间的设置的所述区间之外检测到的两个相邻的拐点的时间坐标值，y(j2)和y(j2)分别为在所述有效明显峰值点右侧单调区域检测到的两个相邻的拐点在滤波前的有效原始回波波形中的振幅值，FWHM表示发射波形的半高宽；如果满足上式的条件，则在所述有效明显峰值点的右侧，在X时间坐标数据中对应序号分别为j2，k2的两个相邻的拐点区域检测到一个叠加的高斯分量。

步骤S3.4，估计所述叠加高斯分量的参数的初始值；步骤S3.4具体包括以下子步骤：

步骤S3.4.1，估计所述有效明显峰值点左侧叠加高斯分量的参数初始值；其中，将该左侧叠加高斯分量的振幅A估计为所述有效明显峰值点左侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点的振幅值；将该叠加高斯分量的中心时间坐标μ值估计为X(k1)的值；将该叠加高斯分量的标准差σ估计为所述有效明显峰值点左侧的位于设定的所述区间之外的第一个相邻的两个拐点之间距离的一半；

步骤S3.4.2，估计所述有效明显峰值点右侧叠加高斯分量的参数初始值；其中，将该右侧叠加高斯分量的振幅A估计为所述明显峰值点右侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点的振幅值；将该叠加高斯分量的中心时间坐标μ值估计为X(j2)的值；将该叠加高斯分量的标准差σ估计为所述有效明显峰值点右侧的位于设定的所述区间之外的第一个相邻的两个拐点之间距离的一半；

步骤S3.5，计算所述滤波后的回波波形中包含的高斯分量个数。

优选的，步骤S4具体包括以下子步骤：

步骤S4.1，基于LM算法，对所述滤波后的回波波形中检测到的高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的参数的初始值进行迭代优化。其中，步骤S4.1具体包括以下子步骤：

步骤S4.1.1，建立拟合函数和在所述的有效原始回波波形(即实际的观测波形)中选择的观测值之间的误差函数和目标函数。

步骤S4.1.2，基于所述目标函数和误差函数，对所述滤波后的回波波形中检

测

到的所有高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的参数进行迭代优化。

步骤S4.2，结合拟合残差RMSE的阈值th_RMSE对所述滤波后的回波波形中检测到高斯分量的参数值进一步的迭代优化。其中，步骤S4.2具体包括以下子步骤：

步骤S4.2.1计算所述有效原始回波波形中的采样点的振幅值和LM算法迭代优化后得到的所述采样点的振幅值的拟合值之间的残差RMSE。

步骤S4.2.2，通过拟合残差RMSE的阈值找到漏掉的地物目标的高斯分量，更新高斯分量的个数并对找到的漏掉的高斯分量进行初始值估计。

步骤S4.2.3，对所述滤波后的回波波形中检测到的高斯分量的参数进一步的迭代优化，直到拟合残差RMSE小于等于拟合残差RMSE阈值th_RMSE时，停止迭代。

本发明还提供了一种计算机产品，其特征在于，通过执行计算机程序来执行权利要求1-9所述的方法。

本发明由于采用了上述的技术方案，使之与现有技术相比，具有以下的优点和积极效果:

(1)本发明结合拐点匹配和迭代优化算法，改进高斯分量的初始检测策略，包括检测明显峰值点对应的高斯分量和叠加高斯分量，从而得到较好的回波波形的高斯分量初始估计。在此基础上，改进高斯分量参数的最优化估计策略，在高斯分量的参数迭代优化过程中，设置拟合残差RMSE的阈值，不断迭代的去寻找漏掉的高斯分量(包括复杂多峰波形中的高斯分量)，并设置无效高斯分量剔除的约束条件，最终高斯分量的分解结果更加可靠、有效、准确。

(2)针对某些原始回波波形而言，在高斯滤波之前，原始回波波形采样点的最大振幅值大于其计算得到的背景噪声阈值，而在高斯滤波之后，该回波波形采样点的最大振幅值小于等于背景噪声阈值。这使得滤波后的回波波形无法用于后续的明显峰值点和拐点的检测。对于这种情形，本发明保留滤波前的原始有效回波波形并将其作为滤波后的回波波形，将其用于后续的明显峰值点和拐点的检测，从而提高了原始回波波形的有效利用率。

(3)本发明通过高斯滤波对原始回波波形进行去除噪声处理以降低背景噪声的影响，对得到的滤波后的回波波形进行明显峰值点和拐点的检测以及拐点的匹配，能够检测到易于区分的地物目标的高斯分量，主要对应于裸地、农田区、建筑区等地表覆盖区域，并通过明显峰值点两侧拐点的匹配来检测地物目标的叠加高斯分量，这种情况对应于在垂直方向上相隔较近的地物目标的区域，主要对应于垂直高度相近的建筑区、植被区等地表覆盖区域，通过以上的检测，完成对滤波后回波波形中检测到的高斯分量的参数初始值估计。之后，通过非线性的迭代优化算法进行高斯分量的参数迭代优化，在这个过程中，通过拟合残差RMSE的阈值寻找漏掉的地物目标的高斯分量，并通过附加的约束条件以剔除无效的高斯分量。基于以上过程，在地物结构复杂的区域，本发明同样能够完成对地物目标的高斯分量的检测和分解，如森林植被区等。所以在不同的地表地物覆盖条件下，如城市建筑区、农田区，森林植被区等，本发明提出的方法均具有一定的适用性和可靠性。

(4)本发明通过对原始回波波形进行高斯滤波等预处理，对得到的滤波后的回波波形进行明显峰值点和拐点的检测以及拐点的匹配，完成对滤波后回波波形中高斯分量的检测和相应的高斯分量的参数初始值估计，并通过非线性迭代优化算法进行高斯分量参数的最优化估计，在这个最优化的过程中，通过拟合残差RMSE的阈值寻找漏掉的地物目标的高斯分量，并结合附加的约束条件以剔除无效的高斯分量。通过以上手段，能够有效的解决原始回波波形数据分解中的难点问题，如含叠加高斯分量的原始回波波形、复杂多峰的原始回波波形以及含弱回波分量的原始回波波形的分解问题。

附图说明

图1为根据本发明实施例的一种星载大光斑激光雷达全波形数据分解方法的流程图。

图2(a)为根据本发明实施例的GLAS测高系统的发射波形。

图2(b)为根据本发明实施例的发射波形对应的原始回波波形，该波形在接收时经过了压缩处理。

图2(c)为根据本发明实施例的解压缩后的原始回波波形及背景噪声估计示意图。

图2(d)为根据本发明实施例的高斯滤波结果示意图。

图2(e)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果示意图。

图3(a)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果，其中示出了叠加的原始回波波形的分解结果。

图3(b)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果，其中示出了含有弱回波分量的原始回波波形的分解结果。

图3(c)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果，其中示出了复杂多峰的原始回波波形的分解结果。

图4(a)为根据本发明实施例的北京平原试验区的大量波形分解结果的R²统计直方图。

图4(b)为根据本发明实施例的黑龙江省大兴安岭试验区的大量波形分解结果的R²统计直方图。

具体实施方式

下面对本发明进行更全面的描述，其中说明本发明的示例性实施例。

本发明实施例提供了一种结合拐点匹配和迭代优化算法的星载大光斑激光雷达全波形数据分解方法，如图1所示，该方法主要包括以下步骤：

步骤S1，对星载大光斑激光雷达全波形数据包括的每个原始回波波形数据进行预处理，以判断每个所述原始回波波形是否为有效原始回波波形，并对所述有效原始回波波形进行一维高斯滤波处理，得到滤波后的回波波形。

星载大光斑激光雷达全波形数据可以从国外的官方网站下载获取，如GLAS数据可以从美国国家冰雪数据中心(NSDIC)网站下载获取，该网站提供的GLAS数据中的GLA01数据包含星载的全波形激光雷达数据，而GLA14产品则包含对全波形数据处理后的结果参数。星载大光斑激光雷达全波形数据包含一定数量待处理的原始回波波形，而原始回波波形在接收的过程中会受到背景噪声的影响，因此，需要对星载大光斑激光雷达全波形数据中的每个原始回波波形分别进行预处理操作。所述预处理操作具体为：判断每个所述原始回波波形是否为有效原始回波波形，并对所述有效原始回波波形进行一维高斯滤波处理，得到滤波后的回波波形。

具体来说，步骤S1具体包括以下子步骤：

步骤S1.1，估计所述星载大光斑激光雷达全波形数据中的每个所述原始回波波形的背景噪声水平。

原始回波波形由N个离散的采样点组成，而其中开头部分和结尾部分的采样点通常是背景噪声的随机起伏产生的，即属于背景噪声信号，为了进行背景噪声水平的估计，传统的方法一般选择原始回波波形中开头部分的采样点进行背景噪声估计，本发明改进之处在于，从原始回波波形的采样点中分别选取开头部分的和结尾部分的各n个采样点(n为自然数)。将所选取的开头部分的和结尾部分的各n个采样点合并后统一进行统计，这样的选择方式能够更充分、更合理的估计背景噪声水平。优选的，n取20，即共选择2n(40)个采样点，计算其均值(μ_n)和标准差(δ_n)，并将μ_n和δ_n分别作为背景噪声的平均值和标准差。

具体来说，所述原始回波波形的背景噪声水平包括所述原始回波波形的背景噪声的平均值μ_n和背景噪声的标准差δ_n，其中：

背景噪声的平均值μ_n的计算公式如下式所示：

其中，N为选取的某个原始回波波形中采样点数的总和，n为在所述原始回波波形的采样点中开头部分的和结尾部分分别选取的采样点数，V_i是第i个采样点的振幅。

背景噪声的标准差δ_n的计算公式如下式所示：

其中，N为选取的某个原始回波波形中采样点数的总和，n为在所述原始回波波形的采样点中开头部分的和结尾部分分别选取的采样点数，V_i是第i个采样点的振幅。

通过以上计算，得到背景噪声的平均值和标准差的良好估计，μ_n和δ_n这两个参数值可以体现背景噪声的水平。

步骤S1.2，基于所述背景噪声水平设置背景噪声阈值，并将所述原始回波波形的采样点的最大振幅值与所述背景噪声阈值进行比较，以判断所述原始回波波形是否为有效原始回波波形，并仅保留所述有效原始回波波形。

在本步骤中设置背景噪声阈值，用于将一个原始回波波形作为整体判断其为有效原始回波波形或噪声回波波形。

根据以下公式设置背景噪声的阈值(th)：

th＝μ_n+m·δ_n (3)

其中，th为背景噪声阈值，μ_n为背景噪声的平均值，δ_n为背景噪声的标准差，m为常数值(经验值，根据实验和实际需求设置)。为了降低背景噪声对回波波形的检测和后续处理的影响，优选的，m取4.5。

如果原始回波波形中采样点的最大振幅值大于背景噪声阈值(th)，则该原始回波波形可被判定为有效原始回波波形，否则将该原始回波波形判定为噪声回波波形，并仅保留有效原始回波波形。需要说明的是：由于每个原始回波波形都是独立的，因此每个回波波形都单独对应一个背景噪声阈值。

根据本发明的优选实施例，考虑到卫星激光载荷的软硬件的不同，对于每个原始回波波形而言，如果在其获取时基于减少数据的冗余考虑而回波波形已经过压缩处理，则在对原始回波波形数据进行预处理之前，首先需要对其进行解压缩处理，通过解压缩处理恢复为完整的回波信号。

例如，GLAS的回波波形数据该数据是经过压缩处理后得到的，因此需要先进行解压缩处理恢复得到完整的回波信号。具体来讲，根据以下公式对原始回波波形数据进行解压缩处理：

其中，p、q、L均为Npq型压缩常数，i为原始回波波形中的采样点的序号，N_gates为原始回波波形中采样点的总个数，t(i)为原始回波波形中第i个采样点在压缩前的序号。

在数据压缩前，GLAS的每个原始回波波形均包含N(系统给定值为1000)个采样点，根据Npq型压缩原理，对于1～N个采样点，连续多个采样点用其平均值代替，直到回波波形的采样点数达到N_gates(系统给定值为544)，如果1000个采样点不足以填充整个回波波形，则对回波波形的剩余采样点以0填充。

对于经过压缩处理的原始回波波形而言，从1～L-1个采样点，每个采样点的振幅值是通过对压缩前的原始回波波形的连续的p个采样点进行平均得到的，从L～N_gates个采样点，每个采样点的振幅值是通过对压缩前的原始回波波形的连续的q个采样点进行平均得到的。

设压缩前GLAS每个原始回波波形数据对应的时间坐标为X(i)，其中i＝1,2，…,N，而对应的振幅值为Y(i)，其中i＝1,2，…,N。设经过压缩处理得到的原始回波波形数据对应的时间坐标为X_Npq(i)，其中i＝1,2，…,N_gates，而对应的振幅值为Y_Npq(i)，其中i＝1,2，…,N_gates。具体来讲，根据上式(4)，经过解压缩过程恢复得到的原始回波波形数据为：

X[t(i)]＝X_Npq(i)，其中i＝1,2，…,N_gates。

Y[t(i)]＝Y_Npq(i)，其中i＝1,2，…,N_gates。

在以上赋值的基础上，对于压缩前原始回波波形的其他采样点的赋值，则有当i＝1,2，…,N_gates时，

如果i≤L-1，则有：

如果L≤i≤N_gates，则有：

其中p、q、L均为Npq型压缩常数，i为原始回波波形中的采样点的序号，N_gates(系统给定值为544)为原始回波波形采样点的总个数，t(i)为原始回波波形中第i个采样点在压缩前的序号，X[t(i)]为压缩前原始回波波形的第t(i)个采样点对应的时间坐标，Y[t(i)]为压缩前原始回波波形的第t(i)个采样点对应的振幅值。根据以上解压缩的过程，恢复得到压缩前的每个原始回波波形。

附图2(a)为根据本发明实施例的一个GLAS测高系统的发射波形，附图2(b)为根据本发明实施例的发射波形对应的原始回波波形，该波形在接收时经过了压缩处理。附图2(c)为根据本发明实施例的解压缩后的原始回波波形及背景噪声估计示意图，具体来说，附图2(c)为对附图2(b)的原始回波波形经过解压缩后得到的原始回波波形，对该波形进行背景噪声平均值和标准差的计算，并进行噪声阈值的估算，其中，黑色水平实线为估算的背景噪声阈值。

步骤S1.3，对所述有效原始回波波形进行一维高斯滤波处理，得到滤波后的回波波形。

由于背景噪声的影响，每个所述原始有效回波波形在振幅上会出现抖动，如果将所述有效原始回波波形直接用于后续步骤S2中的明显峰值点和拐点的检测时，会得到大量具有低振幅和低脉冲宽度的明显峰值点和拐点，这些明显峰值点和拐点基本都是由背景噪声引起的，从而影响所述有效的原始回波波形中包含的高斯分量的检测结果。因此，需要通过对其进行一维的高斯滤波处理，对所述有效回波波形进行平滑，从而降低背景噪声对所述有效回波波形的影响。

高斯滤波是一种低通滤波器，一维的高斯滤波的原理是将每个所述有效原始回波波形和高斯函数进行卷积处理，抑制高频的背景噪声，实现对所述有效原始回波波形的平滑，从而得到高斯滤波后的回波波形。

一维高斯滤波函数g(x,σ)如下式所示：

其中，σ为滤波器宽度，该参数可变，t为滤波器在横坐标轴上的自变量。

设f(t)代表待处理的信号，则一维高斯滤波后的信号s(t,σ)表示为：

s(t,σ)＝f(t)*g(t,σ) (6)

式中*为卷积运算符，σ为滤波器宽度，t为待处理信号的位置。

步骤S1.3具体包括以下子步骤：

步骤S1.3.1，确定一维高斯滤波器的宽度和一维高斯滤波器的窗口大小。

对于星载大光斑激光雷达而言，由于卫星轨道高，且受到大气和地表的影响，背景噪声对接收到的回波波形影响大，造成星载激光雷达的回波信号信噪比低，因此，需要选取更大的滤波器宽度以降低背景噪声的影响，为此，优选的，一维高斯滤波器的宽度(即，σ)选择为所述星载激光雷达的发射波形的半高宽(FWHM)。这样选择是避免了因一维高斯滤波器宽度小而去噪效果差，回波波形平滑程度不够，也避免了因为滤波器宽度过大而使得原始回波波形被过度平滑，损失回波波形的距离精度。例如，对于GLAS测高系统而言，其发射波形的半高宽(FWHM)数值为4ns，则一维高斯滤波器宽度选择为4ns。

对于高斯函数而言，一维高斯滤波器宽度σ确定之后，需要确定一维高斯滤波器的窗口大小。高斯函数是钟形曲线，钟型曲线在区间(μ_n-3σ,μ_n+3σ)范围内的面积占曲线下总面积的99.7％，其中μ_n为高斯函数的均值，因此，选择选择半径为3σ的一维高斯滤波器窗口，这样一维高斯滤波器窗口能够几乎包含高斯函数的有效取值范围。

一维高斯滤波器窗口W大小如下式计算：

W＝1+2·ceil(3·σ) (7)

其中，ceil(x)为函数，其功能是返回大于或者等于指定表达式的最小整数，W为一维高斯滤波器的窗口大小，σ为一维高斯滤波器的宽度，其中，优选的，一维高斯滤波器的宽度σ选择为发射波形的半高宽(FWHM)。

步骤S1.3.2，计算一维高斯滤波器的高斯卷积核。

通过上述步骤S1.3.1，确定了一维高斯滤波器宽度和滤波器窗口大小，本步骤中需要利用这两个参数，计算高斯模板中心位置C，通过(5)式中的滤波函数g(x,σ)，建立一维高斯滤波器的高斯卷积核K(i)，其中i＝1,2,…,W，高斯的卷积核是对连续高斯的离散近似。具体计算如下式所示：

C＝floor(W/2)+1 (8)

其中，K(i)为一维高斯滤波器的高斯卷积核，W为一维高斯滤波器的窗口大小，C为高斯模板中心位置，σ为一维高斯滤波器的宽度，floor(x)为函数，其功能是取不大于x的最大整数。

计算完成之后，将高斯的卷积核中的值进行归一化处理，这样做的目的在于需要保证高斯的卷积核中各个数值的和必须为1。

具体计算如下式所示：

其中，K(i)为一维高斯滤波器的高斯卷积核，W为一维高斯滤波器的窗口大小，C为高斯模板中心位置，σ为一维高斯滤波器的宽度。

步骤S1.3.3，基于确定的一维高斯滤波器的宽度和窗口大小以及高斯卷积核，对所述有效原始回波波形中的所有的采样点进行一维高斯滤波。

通过上述步骤S1.3.1和步骤S1.3.2，得到一维高斯滤波需要的高斯卷积核的各个值，这些值作为输入用于本步骤的一维高斯滤波中。

具体的，根据下述公式对有效原始回波波形中所有的采样点进行一维高斯滤波：

其中，y(i)为滤波前的有效原始回波波形的第i个采样点的振幅值，Y(i)为滤波后的回波波形的第i个采样点的振幅值，L＝floor(W/2)，W为一维高斯滤波器的窗口大小(floor(x)为函数，其功能是取不大于x的最大整数，L通过滤波器窗口大小的一半取整后得到的，因此L可视为滤波器窗口的半径)，K(L+j+1)为一维高斯滤波的高斯卷积核，N为有效原始回波波形中采样点的总个数。

附图2(d)显示了对有效原始回波波形进行一维高斯滤波处理后得到的滤波后的回波波形，可以看出实现了回波波形平滑的效果，降低了背景噪声的影响。

通过以上步骤得到了一维高斯滤波后的回波波形，将其用于步骤S2中的明显峰值点和拐点的检测，可减少检测到的具有低振幅和低脉冲宽度的明显峰值点和拐点的数量，同时总体上也减少了明显峰值点和拐点的检测个数，从而能够提高以下步骤S2的结果的准确性。

根本发明的优选实施例，在步骤S1.3之后，还包括以下步骤S1.4：

步骤S1.4，计算滤波后的回波波形的采样点的最大振幅值，并将所述最大振幅值与该回波波形的背景噪声阈值进行比较，如果所述最大振幅值小于等于所述背景噪声阈值，则保留滤波前的原始有效回波波形并将其作为滤波后的回波波形(即，不对该原始有效回波波形进行滤波处理，而将该滤波前的原始有效回波波形直接作为滤波后的回波波形)，否则，保留滤波后的回波波形。

需要注意的是，对于某些回波波形而言，在高斯滤波之前，原始回波波形采样点的最大振幅值大于噪声阈值，而在高斯滤波之后，该回波波形采样点的最大振幅值小于等于噪声阈值。由于每个原始回波波形都是独立的，因此所述的每个回波波形都单独对应一个背景噪声阈值。对于以上情况，利用高斯滤波后的回波波形无法进行后续步骤中的明显峰值点和拐点的检测。为此，具体分以下情况进行处理：

计算滤波后的回波波形数据Y(i)(其中i＝1,2,…,N)的采样点中的最大振幅值，通过max()函数计算得到，其中max()表示计算一组数据的最大值，即求取max(Y)。

(1)如果max(Y)≤th，则保留高斯滤波前的有效原始回波波形数据y(i)(其中i＝1,2,…,N)，用于下述步骤S2的明显峰值点检测和明显峰值点左右两侧拐点的检测。

(2)如果max(Y)＞th，则保留滤波后的回波波形数据Y(i)(其中i＝1,2,…,N)，用于下述步骤S2的明显峰值点检测和明显峰值点左右两侧的拐点检测。大部分情况下，所述的原始回波波形都属于这种情况。

通过执行步骤S1.4，可以进一步提高有效原始回波波形的利用率，这也是本发明的改进点之一。

步骤S2，检测所述滤波后的回波波形中的明显峰值点和与所述明显峰值点对应的拐点，基于检测到的明显峰值点和与所述明显峰值点对应的拐点，判断检测到的所述明显峰值点是否有效，并估计有效峰值点对应的高斯分量的参数的初始值。

由于滤波后的回波波形中比较容易检测的是明显的峰值点，它对应着一个潜在的地物目标的高斯分量。因此，本步骤中，首先检测滤波后的回波波形的明显峰值点，在此基础上，根据滤波后的回波波形中采样点的二阶导数在拐点两侧异号的原理，检测明显峰值点左右两侧对应的拐点。一般而言，一个明显峰值点匹配左右两个拐点，考虑到回波波形的复杂性，明显峰值点的两侧应至少应有一侧存在拐点。为了避免背景噪声的影响，只对其振幅值在噪声阈值以上的采样点进行明显峰值点和拐点的检测。

通过对滤波后的回波波形中的明显峰值点和其左右两侧的拐点的检测，能够反映有效原始回波波形中包含的高斯分量，因此通过滤波后的回波波形的检测结果可初步判断所述原始回波波形中的高斯分量的个数，利用滤波后的回波波形中检测到的明显峰值点和匹配的拐点，估计所述有效原始回波波形中包含的高斯分量的振幅、中心位置、标准差等参数的初始值。

步骤S2具体包括以下子步骤：

步骤S2.1，通过邻域窗口遍历所述滤波后的回波波形，来检测所述滤波后的回波波形中的明显峰值点。

步骤S2.1具体如下：对于滤波后的回波波形Y(i)(其中i＝1,2,…,N)，利用局部最大值法，选择2r+1个采样点作为邻域窗口大小(其中，r为所选邻域窗口的半径大小，优选的，r＝2)，以所述邻域窗口遍历滤波后的回波波形，若邻域窗口所覆盖的滤波后的回波波形的2r+1个采样点的振幅值均大于该回波波形的背景噪声阈值(需要说明的是：由于每个原始回波波形都是独立的，因此每个回波波形都单独对应一个背景噪声阈值)，滤波后的回波波形的第i个采样点处在邻域窗口的中心位置且该第i个采样点的振幅值Y(i)在2r+1个采样点的振幅值中最大，且在第i个采样点的左侧，从第i-1个采样点开始，依次向左，相应的采样点对应的振幅值逐渐减小；而在第i个采样点的右侧，第i+1个采样点的振幅值小于等于第i个采样点的振幅值，从第i+1个采样点开始，依次向右，相应的采样点对应的振幅值逐渐减小。若满足以上条件，则该第i个采样点为局部最大值点，即，将该第i个采样点确定为明显峰值点。

具体的，根据下式确定明显峰值点：

其中，Y(i)为滤波后的回波波形的第i个采样点的振幅值，r为邻域窗口的半径大小，th为背景噪声阈值，N为滤波后的回波波形的总采样点数。若滤波后回波波形的采样点i满足上述公式(12)，则将该采样点i确定为明显峰值点。

步骤S2.2，计算所述滤波后的回波波形的采样点的二阶导数，并根据相邻的两个采样点的二阶导数是否异号的准则，检测所述明显峰值点左右两侧的拐点。

在检测得到多个明显峰值点后，每个明显峰值点对应着一个潜在的地物目标的高斯分量，对于该高斯分量而言，一般而言，应存在一个明显峰值点和左、右两个拐点，即明显峰值点和左右拐点匹配一个高斯分量。考虑到回波波形的复杂性，明显峰值点的两侧应至少有一侧存在拐点。

本发明的改进之处在于：在每个明显峰值点左右两侧进行拐点检测，计算滤波后的回波波形的采样点的二阶导数，并通过相邻两个采样点的二阶导数是否异号的准则来检测滤波后的有效回波波形中的拐点。而一般的拐点检测方法为：根据拐点两侧曲线凹凸性相反的原理，对于待检测的采样点而言，在其左右两侧区域，选择一定大小的邻域窗口(如6个采样点)，将这些采样点分成两部分，每个部分有3个采样点，在这3个采样点中，将两端采样点的振幅值之和与中间的采样点的振幅值的两倍相减，之后将两个部分分别得到的计算结果值相乘，若该相乘的结果小于0，则找到一个拐点。相比较一般的拐点检测方法，本发明的拐点检测方法更加准确、有效。

具体的，根据以下公式，检测每个明显峰值点左右两侧的拐点：

Y″(i)·Y″(i+1)＜0且Y(i)＞th且Y(i+1)＞th (14)

其中，Y(i)为高斯滤波后的回波波形的第i个采样点的振幅值，Y″(i)为高斯滤波后的回波波形的第i个采样点的二阶导数值，th为背景噪声阈值；如果明显峰值点两侧的大于背景噪声阈值的相邻的两个采样点Y(i)和Y(i+1)的二阶导数Y″(i)和Y″(i+1)满足上式(14)，则Y(i)对应的采样点为拐点。

上式(13)是离散点二阶导数的五点数值微分公式，用于计算采样点的二阶导数值，其中的常数值均为离散点二阶导数的五点数值微分公式所规定的。

上式(14)是拐点的判别公式，用于检测每个明显峰值点左右两侧的拐点。在每个明显峰值点的左右两侧的单调区域内，分别进行拐点的检测，在检测的过程中，只考虑大于背景噪声阈值th的滤波后回波波形的采样点，而小于等于背景噪声阈值th的滤波后回波波形的采样点则不做考虑。需要说明的是：由于每个原始回波波形都是独立的，因此每个回波波形都单独对应一个背景噪声阈值。根据高斯分量的二阶导数在其拐点两侧异号的原理，利用公式(14)进行拐点的计算。

如果明显峰值点的左侧单调区域内的大于背景噪声阈值的相邻的两个采样点Y(i)和Y(i+1)的二阶导数Y″(i)和Y″(i+1)满足上式(14)，则Y(i)对应的采样点为拐点，其对应的时间坐标值为X(i)。同理，如果明显峰值点的右侧单调区域内的相邻的Y(i)和Y(i+1)对应的两个采样点的二阶导数Y″(i)和Y″(i+1)满足上式(14)，则Y(i)对应的采样点为拐点，其时间坐标值为X(i)。假设找到一个拐点其对应的振幅值为Y(i)，则其对应的时间坐标值为X(i)，通过这种方式找到明显峰值点两侧所有的拐点。

利用以上步骤完成所有明显峰值点和拐点的检测。对于每个高斯分量的明显峰值点而言，其两侧检测到的拐点可能会不止一个，假设明显峰值点左侧单调区域内检测到m个拐点，设其在横坐标的时间值分别为：X(i1)，X(i2)，…，X(im)(i1,i2,…,im为检测到的拐点在X数据中的序号)，可将其表示为：

而明显峰值点右侧单调区域内检测到n个拐点，其在横坐标的时间值分别为：X(j1)，X(j2)，…，X(jn)(j1,j2,…,jn为检测到的拐点在X数据中的序号)，可将其表示为：

步骤S2.3，基于检测到的明显峰值点和与所述明显峰值点对应的拐点，判断所述明显峰值点是否有效。

在上述步骤完成之后，检测到所有的明显峰值点和对应的拐点，作为本步骤的输入数据。由于一个明显峰值点并不一定对应地物目标的高斯分量，即明显峰值点有可能是无效的，为此，对明显峰值点的有效性进行判别。

本步骤具体包括以下子步骤：

步骤S2.3.1，如果明显峰值点的左侧拐点和右侧拐点均不存在，则判定所述明显峰值点为无效峰值点；如果仅存在左侧拐点，则执行步骤S2.3.2-S2.3.4，如果仅存在右侧拐点，则执行步骤S2.3.5-S2.3.7，如果同时存在左侧拐点和右侧拐点，则执行步骤S2.3.8-S2.3.10；

如果明显峰值点两侧的左侧拐点和右侧拐点均不存在，则判定该明显峰值点并不对应高斯分量，判别为无效的明显峰值点，将其舍弃。做出以上判断的原理是对于地物目标的高斯分量而言，应包含一个明显峰值点，而明显峰值点左右两侧则分别对应一个左拐点和一个右拐点。在实际情况下，受背景噪声的影响，明显峰值点的左侧和右侧可能会出现多个拐点，如果滤波后的回波波形中明显峰值点两侧的左侧拐点和右侧拐点均不存在，则该明显峰值点很有可能是背景噪声产生的，在这种情况下，应将其舍弃，从而避免引入无效的明显峰值点。

步骤S2.3.2，计算所述明显峰值点的所有左侧拐点的时间坐标的平均值；

其中，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，m为明显峰值点的左侧拐点的总数，

为明显峰值点左侧第k个拐点的时间坐标值。

步骤S2.3.3，基于所述所有左侧拐点的时间坐标的平均值和所述明显峰值点的时间坐标值，计算所述左侧拐点到所述明显峰值点的时间距离；

d_L＝|μ^L-X(i)| (16)

其中，d_L为左侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值。

步骤S2.3.4，判断所述左侧拐点到所述明显峰值点的时间距离是否大于等于星载激光雷达的发射波形的半高宽的一半，如果是，则确定所述明显峰值点为有效峰值点，否则为无效峰值点；

其中，d_L为左侧拐点到明显峰值点的时间距离，FWHM为所述星载激光雷达的发射波形的半高宽。

步骤S2.3.5，计算所述明显峰值点的所有右侧拐点的时间坐标的平均值；

其中，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值，n为明显峰值点的右侧拐点的总数，

为明显峰值点的右侧第k个拐点的时间坐标值。

步骤S2.3.6，基于所述所有右侧拐点的时间坐标的平均值和所述明显峰值点的时间坐标值，计算所述右侧拐点到所述明显峰值点的时间距离；

d_R＝|μ^R-X(i)| (19)

其中，d_R为右侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值。

步骤S2.3.7，判断所述右侧拐点到所述明显峰值点的时间距离是否大于等于星载激光雷达的发射波形的半高宽的一半，如果是，则确定所述明显峰值点为有效峰值点，否则为无效峰值点；

其中，d_R为右侧拐点到明显峰值点的时间距离，FWHM为所述星载激光雷达的发射波形的半高宽。

步骤S2.3.8，分别计算所述明显峰值点的所有左侧拐点和所有右侧拐点的时间坐标的平均值。

如果明显峰值点左侧单调区域内检测到的拐点个数大于等于1，则所述明显峰值点存在左拐点，否则，明显峰值点不存在左拐点。如果明显峰值点存在左拐点，则将明显峰值点左侧所有拐点的时间坐标的平均值作为该明显峰值点的左拐点在横坐标的时间坐标，记为μ^L。同理，如果明显峰值点右侧单调区域内检测到的拐点个数大于等于1，则所述明显峰值点存在右拐点，否则，明显峰值点不存在右拐点。如果明显峰值点存在右拐点，则将明显峰值点右侧所有拐点的时间坐标的平均值作为该明显峰值点的右拐点在横坐标的时间坐标，记为μ^R。具体计算如下式：

其中，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值，m为明显峰值点的左侧拐点的总数，n为明显峰值点的右侧拐点的总数，

为明显峰值点左侧第k个拐点的时间坐标值，

为明显峰值点的右侧第k个拐点的时间坐标值。

步骤S2.3.9，分别计算左侧拐点和右侧拐点到明显峰值点的时间距离。

具体的，根据下式计算左侧拐点和右侧拐点到明显峰值点的时间距离：

其中，d_L为左侧拐点到明显峰值点的时间距离，d_R为右侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值。

步骤S2.3.10，分别判断所述左侧拐点和右侧拐点到所述明显峰值点的时间距离是否大于等于激光雷达的发射波形的半高宽的一半，如果所述左侧拐点和右侧拐点中的至少一者到所述明显峰值点的时间距离大于等于激光雷达的发射波形的半高宽的一半，则确定所述明显峰值点为有效峰值点，否则为无效峰值点。

如果判别为有效明显峰值点，则该峰值点对应一个潜在的地物目标的高斯分量。根据本发明的优选实施例，需要注意的是，对于滤波后的回波波形，且满足max(Y)＞th，如果经过明显峰值点检测和明显峰值点左右两侧拐点检测等步骤，并没有寻找到相应的高斯分量，对于这种情况，本发明采取的做法如下：

计算滤波后的回波波形的采样点的最大振幅值，即max(Y)，并认为该最大振幅值的位置对应一个高斯分量，即，将该最大振幅值对应的采样点的振幅值、时间坐标以及

(标准差的估计值)值作为该高斯分量的初始参数(其中，FWHM为发射波形的半高宽)。这是由于理论上高斯分量的标准差σ与FWHM之间的数值关系为

并进入后续步骤的高斯分量的最优化估计中。

另外，需要说明的是，发射波形和回波波形之间的对应关系在于：回波波形是发射波形和地物目标之间响应后形成的，回波波形是多个地物目标对应的高斯分量的叠加之和，而回波波形中的每个高斯分量反映了对应地物目标的特性。一般而言，回波波形中每个高斯分量的半高宽大于等于发射波形的半高宽。

步骤S2.4，估计确定的所述有效明显峰值点对应的高斯分量的参数的初始值。

在经过上述步骤后，确定了有效的明显峰值点，将其作为本步骤的输入数据。对于每个确定的有效明显峰值点而言，均对应一个潜在的有效原始回波波形中的地物目标的高斯分量。假设Y(i)对应的采样点为明显峰值点，则对其对应的高斯分量进行参数初始值估计。

由于单个地物目标的原始回波波形近似于高斯分布，因此可以通过高斯函数对其进行数学表达，而卫星激光脉冲传输路径上可能会有多个地物，仪器接收到的回波波形则是多个地物目标回波的总和，为此，通过高斯混合模型来表达接收到的每个原始回波波形，如下式所示：

理论上，原始回波波形可看作多个高斯函数的叠加，可用高斯混合模型来表达，因此，将接收到的原始回波波形表达为k个高斯分量的叠加

其中，y为接收到的原始回波波形的近似波形，t为通过离散化的形式接收的原始回波波形的采样点对应的时间坐标数据，k为原始回波波形中高斯分量的个数，A_i、μ_i、σ_i分别为第i个高斯分量的振幅、中心位置(即时间坐标)、标准差，ε为原始回波波形的背景噪声，可用背景噪声的平均值μ_n直接估算。在建立高斯混合模型的基础上，利用本发明提出的全波形数据分解方法完成对波形的拟合和高斯分量的分解。

对于公式(23)中的高斯混合模型而言，其中的k、A_i、μ_i、σ_i(i＝1，2，…，k)是原始回波波形处理最终需要得到的结果参数，而这些参数值是未知的，由于高斯混合模型是非线性的模型，因此需要利用后续步骤中的非线性最小二乘法对这些参数进行不断地迭代和优化，最终得到这些参数的最优值。而在这些参数的迭代和优化之前，需要估计接收到的原始回波波形中包含的高斯分量的个数k的初值，以及每个高斯分量的参数A_i、μ_i、σ_i(i＝1，2，…，k)的初始值，将这些初始值作为初始输入参数。

基于上述的原理，对于所述有效明显峰值点对应的高斯分量而言，它是原始回波波形中包含的高斯分量中的一个，该高斯分量的数学表现形式为：

其中的高斯分量的参数A_i、μ_i、σ_i是未知的，因此需要首先估计这些参数的初始值。

具体如下：

将与有效明显峰值点对应的高斯分量的振幅A_i的初始值估计为所述有效明显峰值点的振幅值Y(i)；

将与有效明显峰值点对应的高斯分量的中心位置μ_i的初始值估计为滤波后的回波波形中所述有效明显峰值点对应的第i个采样点的时间坐标值；

将与有效明显峰值点对应的高斯分量的标准差σ_i的初始值D估计如下：

(1)如果

且右拐点不存在，则D＝d_L；

(2)如果

且左拐点不存在，则D＝d_R；

(3)如果

且

则D＝d_L；

(4)如果

且

则D＝d_R；

(5)如果

且

对此进一步判断，如果d_L≥d_R，则D＝d_R，否则D＝d_L。

上述的判断原理是：根据上述步骤明显峰值点的判别，得到有效的明显峰值点。而每个有效明显峰值点所处位置则对应了一个有效原始回波波形包含的地物目标的高斯分量，结合上述发射波形和回波波形之间的关系说明，则所述明显峰值点对应的高斯分量的左拐点或右拐点到明显峰值点之间的距离d_L和d_R在数值上应大于等于FWHM的一半，但考虑到背景噪声的影响，距离值d_L和d_R至少应有1个大于等于FWHM的一半。

通过以上步骤中对高斯滤波后的回波波形中的明显峰值点和拐点进行检测，能够检测到易于区分的地物目标的高斯分量，对应于裸地、农田区、建筑区等地表覆盖区域。之后完成对应的所述原始回波波形中包含的高斯分量的初始判断和相应的高斯分量的参数初始值估计。

步骤S3，检测所述有效明显峰值点两侧的叠加高斯分量，并估计所述叠加高斯分量的参数初始值。

通过上述步骤S2，检测出高斯滤波后的回波波形中包含的有效明显峰值点，并估计得到所有所述有效明显峰值点对应的高斯分量的振幅、中心位置、标准差等参数的初始值，将这些结果作为输入，用于本步骤的叠加高斯分量检测中。

在明显峰值点的两侧，存在可能的叠加高斯分量，需要进一步的检测回波波形中的叠加高斯分量。叠加高斯分量的产生是由于两个地物目标相隔距离较近，如小于1.5米时，在其对应的原始回波波形中，这两个地物目标对应的高斯分量会叠加在一起，可能会产生只存在一个明显峰值点的叠加回波波形，在这种情况下，则需要进行叠加高斯分量的检测。含叠加高斯分量的回波波形主要对应于垂直高度相近的建筑区、植被区等地表覆盖区域。

在本步骤中，对于所述有效明显峰值点对应的高斯分量而言，设置区间

[X(i)-D,X(i)+D] (24)

其中，X(i)为所述有效明显峰值点对应的高斯分量对应的时间坐标值，D为估计的所述有效明显峰值点对应的高斯分量的标准差的初始值。

设置上式中的区间是认为在该高斯分量的峰值点两侧，潜在的叠加高斯分量应分布在该区间之外，从而避免检测到假的叠加高斯分量。因此，在该高斯分量的峰值点两侧且在上式的区间之外，寻找叠加高斯分量。理论上认为有效的叠加高斯分量所在的位置应距离临近的高斯分量至少1个脉冲宽度(FWHM)的距离,而2·D可近似的等于FWHM。因此，对于明显峰值点左侧存在的叠加高斯分量而言，其存在的拐点的位置都应该是在X(i)-D的左侧位置，对于明显峰值点右侧存在的叠加高斯分量而言，其存在的拐点的位置都应该是在X(i)+D的右侧位置。

本发明提出改进的叠加高斯分量检测方案，通过检测所述有效明显峰值点两侧的相邻的两个拐点，并增加这两个相邻的拐点在原始回波波形中对应的振幅值的比较，判别是否存在叠加高斯分量。而所述有效明显峰值点两侧检测的叠加高斯分量属于对应的所述有效原始回波波形。

步骤S3具体包括以下子步骤：

步骤S3.1，设置区间[X(i)-D,X(i)+D]，其中，X(i)为所述有效明显峰值点对应的高斯分量对应的时间坐标值，D为估计的所述有效明显峰值点对应的高斯分量的标准差的初始值。

步骤S3.2，基于设置的所述区间，检测所述有效明显峰值点左侧的叠加高斯分量。

上述步骤S2中完成对高斯滤波后的回波波形中的所述有效明显峰值点的判别，并对所述有效明显峰值点对应的属于所述的有效原始回波波形中的高斯分量进行初始检测和参数初始值估计。之后，本步骤中，对于每个所述有效明显峰值点对应的高斯分量，在其所述有效明显峰值点的左侧单调区域，检测叠加高斯分量。滤波前的有效原始回波波形数据y(i)(其中i＝1,2,…,N)、滤波后的回波波形数据Y(i)(其中i＝1,2,…,N)和对应的时间坐标数据X(i)(其中i＝1,2,…,N)以及初步得到的有效明显峰值点对应的高斯分量和对应的参数作为输入数据。

在所述有效明显峰值点左侧单调区域的采样点中，而且在[X(i)-D,X(i)+D]的区间外，根据公式(13)和(14)相结合的拐点检测方法进行拐点的检测，如果检测到第一个相邻的两个拐点，记这两个相邻的拐点在X时间坐标数据中对应的序号分别为j1，k1，则其对应的时间坐标值分别为X(j1)和X(k1)。

为了提高叠加高斯分量的准确性，规定X(j1)和X(k1)之间的距离大于发射波形半高宽的一半，并且这两个相邻的拐点中的位于左边的拐点在原始回波波形中的振幅值y(j1)应小于位于右边的拐点在原始回波波形中对应的振幅值y(k1)。

具体来说，根据以下公式检测所述有效明显峰值点左侧的叠加高斯分量：

且y(j1)＜y(k1) (25)

其中，j1、k1分别为在所述有效明显峰值点左侧单调区间的设置的所述区间之外检测到的两个相邻的拐点在X时间坐标数据中对应的位置，X(j1)和X(k1)分别为在所述有效明显峰值点左侧单调区间设置的所述区间之外检测到的两个相邻的拐点的时间坐标值，y(j1)和y(k1)分别为在所述有效明显峰值点左侧单调区域检测到的两个相邻的拐点在滤波前的有效原始回波波形中的振幅值，FWHM表示星载激光雷达的发射波形的半高宽。

选择滤波前的有效原始回波波形中的振幅值加入判别是因为相比较滤波后的回波波形，前者的采样点振幅值没有因为滤波处理而损失，能够作为参考进行辅助判别，这样能够更为有效的判别所述有效明显峰值点左侧是否存在叠加高斯分量。如果满足上式的条件，则在所述有效明显峰值点的左侧，在X时间坐标数据中对应序号分别为j1，k1的两个相邻的拐点区域检测到一个叠加的高斯分量，将这两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点的时间坐标值作为该叠加的高斯分量的中心位置。

步骤S3.3，基于设置的所述区间，检测所述有效明显峰值点右侧的叠加高斯分量。

在上述步骤S2中的有效明显峰值点左侧叠加高斯分量检测完成之后，在有效明显峰值点右侧单调区域，检测叠加高斯分量，其检测方式与上一步骤S2.2.1类似。滤波前的有效原始回波波形数据y(i)(其中i＝1,2,…,N)、滤波后的回波波形数据Y(i)(其中i＝1,2,…,N)和对应的时间坐标数据X(i)(其中i＝1,2,…,N)以及初步得到的有效明显峰值点对应的高斯分量和对应的参数作为输入数据。

具体的，在所述有效明显峰值点右侧单调区域的采样点中，而且在[X(i)-D,X(i)+D]的区间外，根据公式(13)和(14)相结合的拐点检测方法进行拐点的检测，如果出现第一个相邻的两个拐点，记这两个相邻的拐点在X时间坐标数据中对应的位置分别为j2，k2，则其对应的时间坐标值分别为X(j2)和X(k2)。如果满足以下条件：

且y(j2)＞y(k2) (26)

其中，j2、k2分别为在所述有效明显峰值点右侧单调区间的设置的所述区间之外检测到的两个相邻的拐点在X时间坐标数据中对应的位置，X(j2)和X(k2)分别为在所述有效明显峰值点右侧单调区间的设置的所述区间之外检测到的两个相邻的拐点的时间坐标值，y(j2)和y(j2)分别为在所述有效明显峰值点右侧单调区域检测到的两个相邻的拐点在滤波前的有效原始回波波形中的振幅值，FWHM表示发射波形的半高宽。选择滤波前的有效原始回波波形中的振幅值加入判别是因为相比较滤波后的回波波形，前者的振幅值没有因为滤波处理而损失，能够辅助进行判别，这样能够更为有效的判别所述有效明显峰值点右侧是否存在叠加高斯分量。如果满足上式的条件，则在所述有效明显峰值点的右侧，在X时间坐标数据中对应序号分别为j2，k2的两个相邻的拐点区域寻找到一个叠加的高斯分量，将这两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点的时间坐标值作为该叠加的高斯分量的中心位置。

步骤S3.4，估计所述叠加高斯分量的参数的初始值。

在完成所有有效明显峰值点左右两侧叠加高斯分量的检测之后，对叠加高斯分量的参数进行初始值估计。

步骤S3.4具体包括以下子步骤：

步骤S3.4.1，估计所述有效明显峰值点左侧叠加高斯分量的参数初始值。

对于有效明显峰值点左侧叠加高斯分量而言，检测到第一个相邻的两个拐点，记这两个相邻拐点的时间坐标值为X(j1)和X(k1)，则有：

将该左侧叠加高斯分量的振幅A估计为所述有效明显峰值点左侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点的振幅值，即A＝Y(k1)，其中，k1为第一个相邻的两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点在X时间坐标数据中的序号，Y(k1)为第一个相邻的两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点的振幅值；

将该叠加高斯分量的中心时间坐标μ值估计为X(k1)的值，即μ＝X(k1)，其中，k1为所述有效明显峰值点左侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点在X时间坐标数据中的序号，X(k1)为所述第一个相邻的两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点所对应的时间坐标值；

将该叠加高斯分量的标准差σ估计为所述有效明显峰值点左侧的位于设定的所述区间之外的第一个相邻的两个拐点之间距离的一半，即σ＝|X(j1)-X(k1)|/2，其中，j1和k1分别为第一个相邻的两个拐点在X时间坐标数据中的序号，X(j1)和X(k1)分别为第一个相邻的两个拐点中位于左边位置和位于右边位置的拐点所对应的时间坐标值。

步骤S3.4.2，估计所述有效明显峰值点右侧叠加高斯分量的参数初始值。

对于有效明显峰值点右侧叠加高斯分量而言，检测到第一个相邻的两个拐点，记这两个相邻拐点的时间坐标值为X(j2)和X(k2)，则有：

将该右侧叠加高斯分量的振幅A估计为所述明显峰值点右侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点的振幅值，即A＝Y(j2)，其中j2为第一个相邻的两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点在X时间坐标数据中的序号，Y(j2)为第一个相邻的两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点的振幅值；

将该叠加高斯分量的中心时间坐标μ值估计为X(j2)的值，即μ＝X(j2)，其中j2为所述有效明显峰值点右侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点在X时间坐标数据中的序号，X(j2)为第一个相邻的两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点所对应的时间坐标值；

该叠加高斯分量的标准差σ估计为所述有效明显峰值点右侧的位于设定的所述区间之外的第一个相邻的两个拐点之间距离的一半，即σ＝|X(j2)-X(k2)|/2，其中j2和k2分别为第一个相邻的两个拐点在X时间坐标数据中的序号，X(j2)和X(k2)分别为第一个相邻的两个拐点中所对应的时间坐标值。

步骤S3.5，计算所述滤波后的回波波形中检测到的高斯分量的个数。

具体的，计算上述检测到的所述有效明显峰值点的个数和所述有效明显峰值点两侧叠加高斯分量的个数之和，从而得到所述滤波后的回波波形中包含的高斯分量个数。

表1显示了所选的原始回波波形的高斯分量参数如振幅、中心位置、标准差等的初始值估计结果，从表1中看出，经过上述步骤S2，对于所选择的一个原始回波波形案例而言，共检测出9个初始的高斯分量，高斯分量的参数如表1所示。

表1附图1中对应的原始回波波形包含的高斯分量的参数初始值估计结果

在以上步骤完成之后，需要通过以下步骤S4对步骤S2、S3得到的滤波后的回波波形中检测到的所有高斯分量的参数初始值进行最优化估计。对于叠加的回波波形、复杂多峰的回波波形以及弱回波波形而言，通过以上步骤还难以得到原始回波波形中本应包含的所有高斯分量，即存在漏掉地物目标对应的高斯分量的情况。因此，在步骤S4的滤波后的回波波形中检测到的高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的最优化估计中，需要通过波形拟合残差阈值来寻找步骤S2、S3中未能检测出的地物目标的高斯分量，即漏掉的高斯分量。

步骤S4，对所述滤波后的回波波形中检测到的高斯分量的总个数进行更新并对所述高斯分量的参数进行最优化估计。

在上述步骤S2、S3完成之后，得到所有有效明显峰值点对应的高斯分量以及叠加高斯分量，计算得到步骤S2、S3检测到的高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的总个数，并且估计得到所述的每个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数初始值，将以上结果以及有效原始回波波形数据y(i)(其中i＝1,2,…,N)和相应的时间坐标数据X(i)(其中i＝1,2,…,N)作为输入用于本步骤中。

步骤S2、S3中得到的检测结果代表了所述有效原始回波波形中的得到的初始结果，因此，根据步骤S2、S3中的结果，通过本步骤的对这些高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的总个数进行更新以及高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量以及漏掉的地物目标的高斯分量)的参数的最优化估计，并将拟合的回波波形去逼近所述有效原始回波波形，从而得到所述有效原始回波波形中包含的高斯分量个数的最终结果以及高斯分量的参数的最优化结果。

本发明的改进之处在于，在所述滤波后的回波波形中进行地物目标高斯分量的初始判断和叠加高斯分量的检测以及相应的高斯分量的参数初始值估计完成之后，将所有的高斯分量的参数初始值合为一个整体后作为输入数据，在以上所述的高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的参数的迭代优化过程中，设置波形拟合残差RMSE的阈值来不断的寻找步骤S2、S3中未能检测到而导致漏掉的属于所述有效原始回波波形包含的高斯分量，包括复杂多峰回波波形中的高斯分量，并其设置约束条件，剔除无效的高斯分量，直到满足回波波形分解的终止条件，最终完成所述滤波后的回波波形中检测到的高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量以及漏掉的地物目标的高斯分量)的最优化估计，从而得到所述有效原始回波波形中包含的高斯分量总个数的最终结果以及高斯分量的参数的最优化结果。步骤S4具体包括以下子步骤：

步骤S4.1，基于LM算法，对所述滤波后的回波波形中检测到的高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的参数的初始值进行迭代优化。

本步骤中需要用到LM参数迭代优化算法，LM算法是由Levenberg和Marquardt提出的，该方法是在高斯-牛顿法的基础上引入了阻尼，从而形成了阻尼高斯-牛顿法。LM算法首先建立拟合函数与观测值之间的误差函数f(x)，其中x为需要求解的目标参数，之后建立目标函数F(x)：

通过非线性的迭代优化，使得目标函数的结果值达到最小化，从而得到目标参数的最优结果。

根据LM算法原理，有：

(J^TJ+μI)h_lm＝-g (27)

h_lm＝-(J^TJ+μI)^-1g (28)

其中，J为误差函数f(x)的雅克比矩阵，g＝J^Tf(x)，μ为阻尼系数(μ＞0)，h_lm为待求目标参数x的更新量，I为与矩阵J^TJ的行列数相同的单位矩阵。

LM算法迭代的具体过程如下：

(1)设置μ的初始值，设置系数β，设置最大迭代次数m，估计目标参数x的初始值x₀，计算雅克比矩阵J和矩阵J^TJ。设置已执行的迭代次数为j，初始的j＝0。

(2)根据目标参数x的初始值x₀，表示为x_j，计算误差函数f(x_j)的值，计算目标函数F(x_j)的值e_j，以及计算g＝J^Tf的值。

(3)根据上式(27)和(28)，求得x_j的增量h_lm，设x_j+1＝x_j+h_lm，将x_j+1代入误差函数f(x_j+1)，计算目标函数F(x_j+1)的值e_j+1。

(4)如果e_j+1＜e_j，则更新目标参数x的值，即x_j+1＝x_j+h_lm，则设μ＝μ/β；如果e_j+1≥e_j，则不更新目标参数x的值，即x_j+1＝x_j，则设μ＝μ·β。

根据以上迭代过程，(1)-(4)的过程为迭代1次，则更新已执行的迭代次数为j＝j+1，当LM算法的迭代次数j达到最大迭代次数限制时，则整个迭代过程结束，最终得到待求目标参数x的最优结果x_m。

上述为LM算法的迭代优化原理，下面为具体的方法步骤。步骤S4.1具体包括以下子步骤：

步骤S4.1.1，建立拟合函数和在所述的有效原始回波波形(即实际的观测波形)中选择的观测值之间的误差函数和目标函数。

根据上述的高斯分量参数的迭代原理，由于原始回波波形是上述的高斯混合模型需要去逼近的实际观测波形，为了避免背景噪声对所述高斯分量的参数迭代的影响，将所述有效原始回波波形中振幅值大于背景噪声阈值的采样点的振幅值作为观测值，参与到以下对所述有效原始回波波形进行拟合的迭代优化过程，设有效原始回波波形中满足所述条件的观测值有n个，其在时间坐标数据X中对应的序号为p_i，其中i＝1,2,…,n。

建立拟合函数和所述观测值之间的误差函数，具体公式表示如下：

设需要求取的目标参数x＝(A₁,μ₁,σ₁,A₂,μ₂,σ₂,…,A_k,μ_k,σ_k)，则

其中，

为将所述有效原始回波波形中振幅值大于背景噪声阈值的采样点的时间坐标数据作为自变量的拟合函数，f_i(x)为误差函数，代表有效原始回波波形中第p_i个采样点的振幅值和拟合振幅值之间的误差，其中x为待求的目标参数，即x＝(A₁,μ₁,σ₁,A₂,μ₂,σ₂,…,A_k,μ_k,σ_k)；n为从所述有效原始回波波形中选择的采样点作为观测值的个数，y(p_i)为所述有效原始回波波形中第p_i个采样点的振幅值，p_i为从所述有效原始回波波形中选择的采样点作为观测值所对应的序号，其中i＝1,2,…,n，X(p_i)为所述有效原始回波波形中第p_i个采样点的时间坐标值，k为高斯分量的个数，A_j、μ_j、σ_j分别表示第j个高斯分量的振幅、中心位置、标准差，μ_n为背景噪声的平均值。

将所述有效原始回波波形中第p_i个采样点的时间坐标值X(p_i)(其中i＝1,2,…,n)，所述有效原始回波波形中第p_i个采样点的振幅值y(p_i)，以及背景噪声的平均值μ代入上式(30)中，形成含有待求目标参数的误差函数。

根据上述的拟合函数和所述误差函数，建立目标函数F(x)，具体公式如下：

其中，F(x)为目标函数，f_i(x)为误差函数，代表所述有效原始回波波形中所选择的参与拟合的第p_i个采样点的振幅值和拟合振幅值之间的误差，其中x＝(A₁,μ₁,σ₁,A₂,μ₂,σ₂,…,A_k,μ_k,σ_k)，n为从所述有效原始回波波形中选择的采样点作为观测值的个数。

上述公式(31)的含义是计算拟合回波波形和所述的有效原始回波波形的观测值之间的误差平方和的一半。而采用LM算法的目的就是使得目标函数F(x)的值最小化。

步骤S4.1.2，基于所述目标函数和误差函数，对所述滤波后的回波波形中检测到的所有高斯分量(包括有效峰值点对应的高斯分量和上述的叠加高斯分量)的参数进行迭代优化。

根据上述步骤S2获得的k个高斯分量和每个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数的初始值，结合上述步骤S4.1.1中建立的含有待求目标参数的拟合函数

误差函数f_i(x)和目标函数F(x)，以及所述有效原始回波波形中采样点的时间坐标数据X和上述步骤S4.1.1所选择的参与拟合过程的n个观测值(在时间坐标数据X中对应的序号为p_i，其中i＝1,2,…,n)。将这些结果作为输入用于本步骤中。

本步骤具体包括以下子步骤：

步骤S4.1.2.1，根据上述的LM算法原理和迭代过程，设置阻尼系数μ和系数β的初始值、最大迭代次数以及待求的k个高斯分量的参数的初始值。

根据上述的LM算法原理和迭代的具体过程，在采用LM算法之前，设置μ的初始值，一般的，μ取值为0.01，设置系数β的初始值，一般的，β取值为10，设置最大迭代次数m，一般的，m取30。设已执行的迭代次数为r，则在LM算法迭代过程开始之前，有r＝0。

步骤S4.1.2.2，计算所述的拟合函数和所述误差函数的雅克比矩阵。

根据LM算法的迭代优化原理和上述公式(29)和(30)，通过上述步骤S4.1.1建立了含有待求目标参数的误差函数f_i(x)(其中i＝1,2,…,n)，求算所述拟合函数和所述误差函数f_i(x)的雅克比矩阵J(x₀)，J(x₀)的计算方式如下：

其中，

为待求的k个高斯分量的参数x的初始估计值，k为高斯分量的初始个数，f₁(x)，f₂(x)，…，f_n(x)为误差函数，代表从所述有效原始回波波形(也称为观测波形)中所选择的参与拟合过程的n个采样点的振幅值和拟合振幅值之间的误差函数，其中i＝1,2,…,N，

其中，

为误差函数f₁(x)对目标参数中的振幅参数A₁求偏导数得到的值，其中x＝x₀，

为误差函数f₁(x)对目标参数中的中心位置(时间坐标)参数μ₁求偏导数得到的值，其中x＝x₀，

为误差函数f₁(x)对目标参数中的标准差参数σ₁求偏导数得到的值，其中x＝x₀，同理，其他的偏导数的含义类似。

上述的待求的k个高斯分量的参数x的初次估计值x₀，用上述步骤S2检测得到的所述的有效原始回波波形中包含的k个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数的初始值进行估计，即

其中k为步骤S2中得到的高斯分量的初始个数，

为振幅A_i的初次估计值，

为中心位置(时间坐标)μ_i的初次估计值，

为标准差σ_i的初次估计值。

将初始估计值x₀代入上述公式(32)中，计算所述的拟合函数和所述的有效原始回波波形中所选择的采样点的观测值之间的误差函数的雅克比矩阵J(x₀)。之后，计算J(x₀)^TJ(x₀)的值，其中上标T表示对雅克比矩阵J(x₀)求转置。

步骤S4.1.2.3计算步骤S4.1.1所建立的拟合函数和所述有效原始回波波形中所选择的采样点的观测值之间的误差函数的值和目标函数的值。

x的初次估计值

以及所述有效原始回波波形中采样点的时间坐标数据X和上述步骤S4.1.1所选择的参与拟合过程的n个观测值(在时间坐标数据X中对应的序号为p_i，其中i＝1,2,…,n)代入上述公式(30)中，计算拟合函数和所述有效原始回波波形中所选择的采样点的观测值之间的误差函数f_i(x₀)的值，其中i＝1,2,…,n。之后，将得到的所述误差函数f_i(x)的值代入上述公式(31)中，计算上述步骤S4.1.1中的目标函数F(x₀)的值，即

设

以上的计算公式可统一表达为：

其中F(x_r)表示待求的k个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数在LM算法的第r次迭代优化开始时的估计值

所对应的目标函数值，用e_r表示，即e_r＝F(x_r)。

步骤S4.1.2.4，计算待求的k个高斯分量的参数估计值的增量。

根据上式(28)的原理，求得待求k个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数估计值的增量h_lm，即h_lm＝(ΔA₁,Δμ₁,Δσ₁,ΔA₂,Δμ₂,Δσ₂,…,ΔA_k,Δμ_k,Δσ_k)，具体计算公式如下：

h_lm＝-(J(x_r)^TJ(x_r)+μI)^-1*J(x_r)^T*[f₁(x_r)，f₂(x_r)，…，f_n(x_r)]^T (34)

其中h_lm为待求k个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数估计值的增量，J(x_r)为在LM算法的第r次迭代优化开始时，待求的k个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数的估计值x_r所对应的雅克比矩阵，计算方式与初次估计值x₀所对应的雅克比矩阵相似，只是参数估计值不同。

f₁(x_r)，f₂(x_r)，…，f_n(x_r)为误差函数，代表所述的有效原始回波波形中所选择的参与拟合的采样点的振幅值和拟合振幅值之间的误差。设x_r+1＝x_r+h_lm，将x_r+1代入上述公式(30)中计算误差函数f_i(x_r+1)(其中i＝1,2,…,n)，之后，计算目标函数F(x_r+1)的值，即

设e_r+1＝F(x_r+1)。

步骤S4.1.2.5判定是否更新待求k个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数的估计值。

如果e_r+1＜e_r，则更新目标参数x的估计值，即x_r+1＝x_r+h_lm，更新后

根据LM算法的迭代原理，则有μ＝μ/β；如果e_r+1≥e_r，则不更新目标参数A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数的值，即x_r+1＝x_r，

根据LM算法的迭代原理，则有μ＝μ·β。

根据以上迭代过程，步骤S4.1.2.1-步骤S4.1.2.5的过程为LM算法迭代1次，则更新已执行的迭代次数为r＝r+1，当LM算法已执行的迭代次数r等于最大迭代次数m时，一般的，m取30，则整个迭代过程结束。根据LM算法的以上迭代过程，得到目标参数x的优化结果，即

作为优化后的振幅值

中心位置(时间坐标)值

和标准差值

即得到所述的k个高斯分量的A_i、μ_i、σ_i(其中i＝1,2,…,k)等参数的优化结果，并且得到所述的有效原始回波波形的近似波形

可表示为：

其中，μ_n为背景噪声的平均值，k为高斯分量个数，

和

分别为LM算法迭代优化后的第i个高斯分量的振幅、中心位置(即时间坐标)、标准差参数。

步骤S4.2，结合拟合残差RMSE的阈值th_RMSE对所述滤波后的回波波形中检测到的高斯分量的参数值进一步的迭代优化。

通过上述步骤S4.1，得到k个高斯分量的优化后的振幅值

中心位置(即时间坐标)值

标准差值

以及所述有效原始回波波形的拟合波形的拟合振幅值

其中i＝1,2,…,N，将这些结果作为输入用于本步骤中。由于上述步骤S2、S3中在进行所述滤波后的回波波形中检测到的高斯分量的初始判断和所检测到的叠加高斯分量的检测后，有可能存在地物目标的高斯分量因未被检测出来而导致漏掉的情况，因此，在步骤S4.1之后，还需要通过本步骤的设置拟合残差RMSE的阈值th_RMSE去判断是否存在漏掉地物目标的高斯分量的情况，如果存在漏掉地物目标的高斯分量，则需要更新所述滤波后的回波波形中检测到的高斯分量总的个数和参数值，对所述滤波后的回波波形中检测到的高斯分量的参数值进一步迭代优化。

步骤S4.2具体包括以下子步骤：

步骤S4.2.1计算所述有效原始回波波形中的采样点的振幅值和LM算法迭代优化后得到的所述采样点的振幅值的拟合值之间的残差RMSE。

在上述步骤S4.1的LM算法迭代优化过程之后，计算所述有效原始回波波形中的采样点的振幅值和LM算法迭代优化后得到的所述采样点的振幅值的拟合值之间的残差RMSE，计算如下式所示：

其中，n表示参与拟合的所述的有效原始回波波形数据中采样点的个数，只选择所述的有效原始回波波形中大于等于背景噪声阈值的采样点进行拟合近似，y(i)为所述有效原始回波波形中第i个采样点的振幅值，

为经过LM算法迭代优化后得到的第i个所述有效原始回波波形采样点的拟合振幅值。

通过上述步骤后，设置拟合残差阈值th_RMSE，用于确定步骤S4.1中漏掉的高斯分量。具体的，对于星载激光测高系统而言，激光回波信号的信噪比低且易受噪声影响，设置拟合残差阈值th_RMSE，阈值大小如下式所示：

th_RMSE＝4.5·δ_n (37)

其中，δ_n表示背景噪声的标准差。

步骤S4.2.2，通过拟合残差RMSE的阈值找到漏掉的地物目标的高斯分量，更新高斯分量的个数并对找到的漏掉的高斯分量进行初始值估计。

根据上述步骤S4.1中得到的拟合残差RMSE，以及设置的拟合残差RMSE阈值th_RMSE，确定是否漏掉有效的高斯分量，并进行高斯分量个数的更新，并重新进行高斯分量的优化。具体如下：

如果残差RMSE＞th_RMSE,则说明找到一个新的高斯分量，计算

的误差，找出误差的最大值在所述有效原始回波波形中的所在位置，记为ii，并在该位置添加一个新的高斯分量，将高斯分量的振幅、中心位置、标准差分别设为A_ii，μ_ii，σ_ii，对这3个参数进行赋初值：

A_ii＝y(ii)，μ_ii＝X(ii)，

其中A_ii为找到的新的高斯分量的振幅初值，y(ii)为找到的新的高斯分量在所述有效原始回波波形的采样点中的第ii位置对应的振幅值，μ_ii为找到的新的高斯分量的中心位置的初值，X(ii)为找到的新的高斯分量在所述有效原始回波波形的采样点中的第ii位置对应的时间坐标，σ_ii为找到的新的高斯分量的标准差初值。

对滤波后的回波波形中检测到的高斯分量的个数进行更新，即将所述的高斯分量的总个数增加1，即k＝k+1个。以上对找到的新的高斯分量的初值估计原理在于，新的高斯分量的振幅初值可用其在所述有效原始回波波形中的所在位置对应的振幅值进行估计，新的高斯分量的中心位置初值可用其在所述有效原始回波波形中的所在位置对应的时间坐标进行估计，而新的高斯分量的标准差初值则可用星载激光雷达的发射波形的半高宽除以

进行估计，这是由于理论上高斯分量的标准差σ与FWHM之间的数值关系为

步骤S4.2.3，对所述滤波后的回波波形中检测到的高斯分量的参数进一步的迭代优化，直到拟合残差RMSE小于等于拟合残差RMSE阈值th_RMSE时，停止迭代。

上述步骤S4.2.1和步骤S4.2.2确定了漏掉的高斯分量，并得到该高斯分量的参数初值，更新了滤波后的回波波形中检测到的高斯分量的个数，将更新后的滤波后的回波波形中检测到的高斯分量的个数和高斯分量的参数值作为输入，用于本步骤中。将该高斯分量加入到公式(23)的高斯混合模型中，再次进行上述步骤S4.1.2的LM算法的迭代优化过程。

通过上述步骤S4.2.1和步骤S4.2.2，得到再次迭代优化后的高斯分量个数和高斯分量的参数。之后进行如下判断：

(1)如果RMSE＞th_RMSE，则加入新的高斯分量，更新高斯分量总的个数和参数。之后，重复步骤S4.1中的子步骤S4.1.2的高斯分量的参数迭代过程，完成之后，继续执行上述步骤S4.2.1、步骤S4.2.2和步骤S4.2.3，在高斯分量参数的优化的过程中加入下述的附加的约束条件，用于剔除无效的高斯分量，并更新高斯分量总的个数和高斯分量的参数。直到RMSE≤th_RMSE时，最终所述的有效原始回波波形中包含的所有高斯分量的参数迭代优化的整个过程停止。

(2)如果RMSE≤th_RMSE，则加入下述的附加的约束条件，用于剔除无效的高斯分量，并更新高斯分量个数和高斯分量的参数。之后，重新进行步骤S4.1中的子步骤S4.1.2的高斯分量的参数迭代过程，直到所有的无效高斯分量均被剔除，整个的高斯分量的参数迭代过程终止。

附加的约束条件：

在上述步骤S4.2的高斯分量个数更新和高斯分量的参数的迭代优化过程中，设置附加的约束条件，剔除无效的高斯分量，条件如下：

(1)高斯分量之间的最小距离>最小距离阈值，单位ns；

(2)高斯分量的最小振幅>噪声阈值th；

(3)高斯分量的最小脉冲宽度>发射波形的半高宽(FWHM)；

(4)高斯分量个数n≤最大个数阈值。

设置上述条件的原因在于，设置条件(1)的最小距离阈值是因为受噪声影响，如果高斯分量之间距离设置太小，存在非真实地物的信号被误判为有效的高斯分量的情况，因此需要设置(1)的条件，提高高斯分量的分解准确性，优选的，最小距离阈值为10ns，对应1.5米的距离。设置条件(2)是为了避免背景噪声对高斯分量检测的影响。设置条件(3)是为了避免无效高斯分量的影响，设置条件(4)是考虑到卫星激光传输路径上地物垂直分布的特性以及回波波形信号本身特性的限制，从而提出的限制条件，通常最大高斯分量个数阈值取值为6。如果有高斯分量不满足任意一个条件，则判别为无效的高斯分量并进行剔除，之后重新进行回波波形的拟合和参数的优化。通过以上的对高斯分量的参数值进一步的迭代优化步骤，在地物结构复杂的区域，本发明同样能够完成对地物目标的高斯分量的检测和分解，如森林植被区等。通过以上步骤S1、S2、S3、S4，最终本发明完成对星载大光斑激光雷达全波形数据的分解，从所述的有效原始回波波形中，将其中包含的每个地物目标的高斯分量从中全部分解出来，从而提取得到每个所述的有效原始回波波形中的每个高斯分量的振幅A_i、中心位置μ_i、标准差σ_i(其中i＝1,2,…,k)等参数的最终优化结果和高斯分量个数的最终结果。

附图2(e)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果示意图，图中显示了经过步骤S3得到的最终的原始回波波形分解结果。经过本步骤的回波波形中高斯分量参数的最优化估计，最终分解得到6个高斯分量，获得高斯分量的个数，以及每个高斯分量的振幅、中心位置、标准差等参数。附图2(e)对应的原始回波波形的分解结果值如表2所示。

表2附图2对应的原始回波波形的分解结果

以GLAS数据为例，附图3(a)-(c)显示了3种较难处理的原始回波波形的分解结果典型样例，该图中显示了叠加回波波形、弱回波波形以及复杂多峰回波波形等3种较难处理的回波波形。附图3(a)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果典型样例之一，代表了叠加的原始回波波形的分解结果。附图3(b)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果典型样例之二，代表了含有弱回波分量的原始回波波形的分解结果。附图3(c)为根据本发明实施例的解压缩后的原始回波波形对应的波形分解结果典型样例之三，代表了复杂多峰的原始回波波形的分解结果。可以看出，本发明均有效的将这些回波波形进行了分解，得到了每个原始回波波形包含的高斯分量及其参数。

为了验证本发明所提出的星载大光斑激光雷达全波形数据分解方法，以GLAS的回波波形数据作为输入数据，分别选择北京平原地区和黑龙江省大兴安岭地区作为试验区域，进行了大量的全波形数据分解试验。表3和表4分别显示了北京平原地区和黑龙江大兴安岭地区的波形分解结果的高斯分量个数统计。其中NSIDC(National Snow and IceCenter)算法代表GLAS的官方波形数据分解方法。

从中可看出，在北京平原地区，本发明提出的波形分解方法和NSIDC方法相比，含有高斯分量个数为1的原始回波波形个数减少了444个，含有高斯分量个数为2、3、4、5、6的原始回波波形个数均增多。以上结果表明：本发明提出的波形分解方法在平原地区能够获得更多的地物波形分量。在黑龙江大兴安岭地区，含有高斯分量个数为1的原始回波波形个数减少了309个，含有高斯分量个数在4以上的原始回波波形个数均减少，含有高斯分量个数为2、3的原始回波回波波形个数增多。以上结果表明：本发明提出的波形分解方法在保证高斯分量分解的拟合精度的同时，在森林植被区域，能够以更少的高斯分量来代表激光光斑内地物的垂直分布情况。

表3北京平原地区的波形分解结果的高斯分量个数统计

表4黑龙江大兴安岭地区的波形分解结果的高斯分量个数统计

附图4(a)-(b)显示了试验对应的波形分解结果的R²统计直方图。附图4(a)为根据本发明实施例的大量波形分解结果的R²统计直方图，图中代表的为北京平原试验区。附图4(b)为根据本发明实施例的大量波形分解结果的R²统计直方图，图中代表的为黑龙江省大兴安岭试验区。

从中可以看出，波形拟合决定系数R²主要分布在0.95以上，表明了本发明的波形分解方法的拟合结果精度高。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims (9)

1.一种星载大光斑激光雷达全波形数据分解方法，其特征在于，所述方法包括以下步骤：

步骤S1，对星载大光斑激光雷达全波形数据包括的每个原始回波波形数据进行预处理，以判断每个所述原始回波波形是否为有效原始回波波形，并对所述有效原始回波波形进行一维高斯滤波处理，得到滤波后的回波波形；

步骤S2，检测所述滤波后的回波波形中的明显峰值点和与所述明显峰值点对应的拐点，基于检测到的明显峰值点和与所述明显峰值点对应的拐点，判断检测到的所述明显峰值点是否有效，并估计有效峰值点对应的高斯分量的参数的初始值；

步骤S3，检测所述有效明显峰值点两侧的叠加高斯分量，并估计所述叠加高斯分量的参数初始值；

步骤S4，对所述滤波后的回波波形中检测到的高斯分量的总个数进行更新并对所述高斯分量的参数进行最优化估计；其中：

步骤S1具体包括以下子步骤：

步骤S1.1，估计所述星载大光斑激光雷达全波形数据中的每个所述原始回波波形的背景噪声水平；

步骤S1.2，基于所述背景噪声水平设置背景噪声阈值，并将所述原始回波波形的采样点的最大振幅值与所述背景噪声阈值进行比较，以判断所述原始回波波形是否为有效原始回波波形，并仅保留所述有效原始回波波形；

步骤S1.3，对所述有效原始回波波形进行一维高斯滤波处理，得到滤波后的回波波形；

步骤S1.4，计算滤波后的回波波形的采样点的最大振幅值，并将所述最大振幅值与该回波波形的背景噪声阈值进行比较，如果所述最大振幅值小于等于所述背景噪声阈值，则保留滤波前的原始有效回波波形，否则，保留滤波后的回波波形。

2.根据权利要求1所述的方法，其特征在于，在步骤S1.1中，所述原始回波波形的背景噪声水平包括所述原始回波波形的背景噪声的平均值μ_n和背景噪声的标准差δ_n，其中：

背景噪声的平均值μ_n的计算公式如下式所示：

其中，N为选取的某个原始回波波形中采样点数的总和，n为在所述原始回波波形的采样点中开头部分的和结尾部分分别选取的采样点数，n取为20，V_i是第i个采样点的振幅；

背景噪声的标准差δ_n的计算公式如下式所示：

其中，N为选取的某个原始回波波形中采样点数的总和，n为在所述原始回波波形的采样点中开头部分的和结尾部分分别选取的采样点数，n取为20，V_i是第i个采样点的振幅。

3.根据权利要求1所述的方法，其特征在于，在步骤S1.1中，根据以下公式设置背景噪声的阈值th：

th＝μ_n+m·δ_n

其中，th为背景噪声阈值，μ_n为背景噪声的平均值，δ_n为背景噪声的标准差，m为常数值，m取为4.5；

其中，如果原始回波波形中采样点的最大振幅值大于背景噪声阈值th，则该原始回波波形可被判定为有效原始回波波形，否则将该原始回波波形判定为噪声回波波形，并仅保留有效原始回波波形。

4.根据权利要求1所述的方法，其特征在于，步骤S1.3具体包括以下子步骤：

步骤S1.3.1，确定一维高斯滤波器的宽度和一维高斯滤波器的窗口大小，其中，一维高斯滤波器的宽度σ选择为所述星载激光雷达的发射波形的半高宽FWHM；一维高斯滤波器窗口W大小如下式计算：

W＝1+2·ceil(3·σ)

其中，ceil(x)为函数，其功能是返回大于或者等于指定表达式的最小整数，W为一维高斯滤波器的窗口大小，σ为一维高斯滤波器的宽度；

步骤S1.3.2，计算一维高斯滤波器的高斯卷积核，具体计算如下式所示：

C＝floor(W/2)+1

其中，K(i)为一维高斯滤波器的高斯卷积核，W为一维高斯滤波器的窗口大小，C为高斯模板中心位置，σ为一维高斯滤波器的宽度，floor(x)为函数，其功能是取不大于x的最大整数；

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

其中，K(i)为一维高斯滤波器的高斯卷积核，W为一维高斯滤波器的窗口大小，C为高斯模板中心位置，σ为一维高斯滤波器的宽度；

步骤S1.3.3，基于确定的一维高斯滤波器的宽度和窗口大小以及高斯卷积核，对所述有效原始回波波形中的所有的采样点进行一维高斯滤波；

具体的，根据下述公式对有效原始回波波形中所有的采样点进行一维高斯滤波：

其中，y(i)为滤波前的有效原始回波波形的第i个采样点的振幅值，Y(i)为滤波后的回波波形的第i个采样点的振幅值，L＝floor(W/2)，W为一维高斯滤波器的窗口大小，K(L+j+1)为一维高斯滤波的高斯卷积核，N为有效原始回波波形中采样点的总个数。

5.根据权利要求1所述的方法，其特征在于，步骤S2具体包括以下子步骤：

步骤S2.1，通过邻域窗口遍历所述滤波后的回波波形，来检测所述滤波后的回波波形中的明显峰值点；

步骤S2.2，计算所述滤波后的回波波形的采样点的二阶导数，并根据相邻的两个采样点的二阶导数是否异号的准则，检测所述明显峰值点左右两侧的拐点；

步骤S2.3，基于检测到的明显峰值点和与所述明显峰值点对应的拐点，判断所述明显峰值点是否有效；

步骤S2.4，估计确定的所述有效明显峰值点对应的高斯分量的参数的初始值。

6.根据权利要求5所述的方法，其特征在于，步骤S2.3具体包括以下子步骤：

步骤S2.3.1，如果明显峰值点的左侧拐点和右侧拐点均不存在，则判定所述明显峰值点为无效峰值点；如果仅存在左侧拐点，则执行步骤S2.3.2-S2.3.4，如果仅存在右侧拐点，则执行步骤S2.3.5-S2.3.7，如果同时存在左侧拐点和右侧拐点，则执行步骤S2.3.8-S2.3.10；

步骤S2.3.2，计算所述明显峰值点的所有左侧拐点的时间坐标的平均值；

其中，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，m为明显峰值点的左侧拐点的总数，

为明显峰值点左侧第k个拐点的时间坐标值；

步骤S2.3.3，基于所述所有左侧拐点的时间坐标的平均值和所述明显峰值点的时间坐标值，计算所述左侧拐点到所述明显峰值点的时间距离；

d_L＝|μ^L-X(i)|

其中，d_L为左侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值；

步骤S2.3.4，判断所述左侧拐点到所述明显峰值点的时间距离是否大于等于星载激光雷达的发射波形的半高宽的一半，如果是，则确定所述明显峰值点为有效峰值点，否则为无效峰值点；

其中，d_L为左侧拐点到明显峰值点的时间距离，FWHM为所述星载激光雷达的发射波形的半高宽；

步骤S2.3.5，计算所述明显峰值点的所有右侧拐点的时间坐标的平均值；

其中，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值，n为明显峰值点的右侧拐点的总数，

为明显峰值点的右侧第k个拐点的时间坐标值；

步骤S2.3.6，基于所述所有右侧拐点的时间坐标的平均值和所述明显峰值点的时间坐标值，计算所述右侧拐点到所述明显峰值点的时间距离；

d_R＝|μ^R-X(i)|

其中，d_R为右侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值；

步骤S2.3.7，判断所述右侧拐点到所述明显峰值点的时间距离是否大于等于星载激光雷达的发射波形的半高宽的一半，如果是，则确定所述明显峰值点为有效峰值点，否则为无效峰值点；

其中，d_R为右侧拐点到明显峰值点的时间距离，FWHM为所述星载激光雷达的发射波形的半高宽；

步骤S2.3.8，分别计算所述明显峰值点的所有左侧拐点和所有右侧拐点的时间坐标的平均值，具体计算如下式：

其中，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值，m为明显峰值点的左侧拐点的总数，n为明显峰值点的右侧拐点的总数，

为明显峰值点左侧第k个拐点的时间坐标值，

为明显峰值点的右侧第k个拐点的时间坐标值；

步骤S2.3.9，分别计算左侧拐点和右侧拐点到明显峰值点的时间距离；具体的，根据下式计算左侧拐点和右侧拐点到明显峰值点的时间距离：

其中，d_L为左侧拐点到明显峰值点的时间距离，d_R为右侧拐点到明显峰值点的时间距离，X(i)为明显峰值点Y(i)对应的时间坐标，μ^L为明显峰值点的所有左侧拐点的时间坐标的平均值，μ^R为明显峰值点的所有右侧拐点的时间坐标的平均值；

步骤S2.3.10，分别判断所述左侧拐点和右侧拐点到所述明显峰值点的时间距离是否大于等于激光雷达的发射波形的半高宽的一半，如果所述左侧拐点和右侧拐点中的至少一者到所述明显峰值点的时间距离大于等于激光雷达的发射波形的半高宽的一半，则确定所述明显峰值点为有效峰值点，否则为无效峰值点。

7.根据权利要求1所述的方法，其特征在于，步骤S3具体包括以下子步骤：

步骤S3.1，设置区间[X(i)-D,X(i)+D]，其中，X(i)为所述有效明显峰值点对应的高斯分量对应的时间坐标值，D为估计的所述有效明显峰值点对应的高斯分量的标准差的初始值；

步骤S3.2，基于设置的所述区间，根据以下公式检测所述有效明显峰值点左侧的叠加高斯分量：

且y(j1)＜y(k1)

其中，j1、k1分别为在所述有效明显峰值点左侧单调区间的设置的所述区间之外检测到的两个相邻的拐点在X时间坐标数据中对应的位置，X(j1)和X(k1)分别为在所述有效明显峰值点左侧单调区间设置的所述区间之外检测到的两个相邻的拐点的时间坐标值，y(j1)和y(k1)分别为在所述有效明显峰值点左侧单调区域检测到的两个相邻的拐点在滤波前的有效原始回波波形中的振幅值，FWHM表示星载激光雷达的发射波形的半高宽；如果满足上式的条件，则在所述有效明显峰值点的左侧，在X时间坐标数据中对应序号分别为j1，k1的两个相邻的拐点区域检测到一个叠加的高斯分量；

步骤S3.3，基于设置的所述区间，根据以下公式检测所述有效明显峰值点右侧的叠加高斯分量：

且y(j2)＞y(k2)

其中，j2、k2分别为在所述有效明显峰值点右侧单调区间的设置的所述区间之外检测到的两个相邻的拐点在X时间坐标数据中对应的位置，X(j2)和X(k2)分别为在所述有效明显峰值点右侧单调区间的设置的所述区间之外检测到的两个相邻的拐点的时间坐标值，y(j2)和y(k2)分别为在所述有效明显峰值点右侧单调区域检测到的两个相邻的拐点在滤波前的有效原始回波波形中的振幅值，FWHM表示发射波形的半高宽；如果满足上式的条件，则在所述有效明显峰值点的右侧，在X时间坐标数据中对应序号分别为j2，k2的两个相邻的拐点区域检测到一个叠加的高斯分量；

步骤S3.4，估计所述叠加高斯分量的参数的初始值；步骤S3.4具体包括以下子步骤：

步骤S3.4.1，估计所述有效明显峰值点左侧叠加高斯分量的参数初始值；其中，将该左侧叠加高斯分量的振幅A估计为所述有效明显峰值点左侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于右边位置的更靠近所述有效明显峰值点的拐点的振幅值；将该叠加高斯分量的中心时间坐标μ值估计为X(k1)的值；将该叠加高斯分量的标准差σ估计为所述有效明显峰值点左侧的位于设定的所述区间之外的第一个相邻的两个拐点之间距离的一半；

步骤S3.4.2，估计所述有效明显峰值点右侧叠加高斯分量的参数初始值；其中，将该右侧叠加高斯分量的振幅A估计为所述明显峰值点右侧的位于设定的所述区间之外的第一个相邻的两个拐点中位于左边位置的更靠近所述有效明显峰值点的拐点的振幅值；将该叠加高斯分量的中心时间坐标μ值估计为X(j2)的值；将该叠加高斯分量的标准差σ估计为所述有效明显峰值点右侧的位于设定的所述区间之外的第一个相邻的两个拐点之间距离的一半；

步骤S3.5，计算所述滤波后的回波波形中检测到的高斯分量的个数。

8.根据权利要求1所述的方法，其特征在于，步骤S3具体包括以下子步骤：

步骤S4.1，基于LM算法，对所述滤波后的回波波形中检测到的高斯分量的参数的初始值进行迭代优化，所述高斯分量包括有效峰值点对应的高斯分量和上述的叠加高斯分量；其中，步骤S4.1具体包括以下子步骤：

步骤S4.1.1，建立拟合函数和在所述有效原始回波波形中选择的观测值之间的误差函数和目标函数；

步骤S4.1.2，基于所述目标函数和误差函数，对所述滤波后的回波波形中检测到的所有高斯分量的参数进行迭代优化，所述高斯分量包括有效峰值点对应的高斯分量和上述的叠加高斯分量；

步骤S4.2，结合拟合残差RMSE的阈值th_RMSE对所述滤波后的回波波形中包含的高斯分量的参数值进一步的迭代优化；其中，步骤S4.2具体包括以下子步骤：

步骤S4.2.1计算所述滤波后的回波波形中的采样点的振幅值和LM算法迭代优化后得到的所述采样点的振幅值的拟合值之间的残差RMSE；

步骤S4.2.2，通过拟合残差RMSE的阈值找到漏掉的地物目标的高斯分量，更新高斯分量的个数并对找到的漏掉的高斯分量进行初始值估计；

步骤S4.2.3，对所述滤波后的回波波形中检测到的高斯分量的参数进一步的迭代优化，直到拟合残差RMSE小于等于拟合残差RMSE阈值th_RMSE时，停止迭代。

9.一种计算机产品，其特征在于，通过执行计算机程序来执行权利要求1-8所述的方法。

Images