CN104280725A - 全波形激光雷达数据波形分解方法 - Google Patents

Abstract

本发明提供了一种全波形激光雷达数据波形分解方法。该方法利用AIC进行组分数估计，采用LM进行遗失组分探测，可快速获取后向散射回波波形中的组分数信息，算法简单，探测效率高，与一般的波形分解方法相比，兼顾模型组分数的初始化准确度和波形分解的完整性。

Description

全波形激光雷达数据波形分解方法

技术领域

本发明涉及遥感对地观测技术领域，尤其涉及一种全波形激光雷达数据波形分解方法。

背景技术

作为一种新型的遥感对地观测技术，机载激光雷达可直接获取地面目标三维信息，在数字城市、应急减灾、环境监测等各方面应用广泛，是传统摄影测量技术的有效补充。相对于传统激光雷达系统获取的地物相关信息量有限，全波形激光雷达系统以非常小的采样间隔记录激光回波的全部信息并以数字化存储，其本质上是激光束传播方向上地物的各次回波信息的叠加，能提供丰富的垂直结构信息和地物的反射信息，经过处理可以得到能够反映地物固有特性的潜在特征。

由于目标的散射特性和同激光脉冲的作用机理等因素的不同，全波形激光雷达记录的回波波形中回波信号混杂叠加在一起，不能直接获取目标的有价值信息，需要进行全波形数据处理，涵盖预处理、波形建模、波形分解、组分特征提取等。其中，精确获取各组分信息的关键在波形分解。目前常用的方法有基于期望最大化(Expectation Maximum，EM)法、反卷积法和B样条法等，其中EM法较为成熟，该算法需要进行模型初始化，模型初始化参数包括组分数和各组分的振幅、波形宽度等参数。以下给出应用EM法进行全波形激光雷达数据波形分解的具体实例。

参考文件1(Persson A，Soderman U，Topel J，et al.Visualization AndAnalysis Of Full-Waveform Airborne Laser Scanner Data.ISPRS WG III/3，III/4，V/3Workshop″Laser Scanning2005″，2005[C]，Enschede，TheNetherlands：Vosselman G，Brenner C，2005：103-108)提供了一种全波形激光雷达数据波形分解方法。此方法采用赤池信息量准则(AkaikeInformation Criterion，AIC)对组分数进行估计，进一步采用广义高斯函数为核函数，基于非线性最小二乘(Nonlinear least squares，NLS)方法对后向散射回波波形的模型参数进行初始化，完成波形模型参数初始化后采用EM算法完成波形分解。然而，上述方法仅通过EM算法中似然方程的收敛来判断波形分解是否完成，当后向散射回波波形较复杂或去噪未彻底时，可能会存在遗失组分的情况

为了改善遗失组分的情况，参考文件2(Chauve A，Vega C，Durrieu S，etal.Advanced Full-Waveform Lidar Data Echo Detection：Assessing QualityOf Derived Terrain And Tree Height Models In An Alpine Coniferous Forest[J].International Journal Of Remote Sensing，2009，30(19)：5211-5228)提出了在模型参数初始化时采用NLS方法，波形分解过程基于EM算法并结合列文伯格-马夸尔特法(Levenberg-Marquardt，简称LM)思想，可对激光脉冲后向散射回波波形进行增强探测的迭代算法，能有效地探测到回波中遗失的组分。

虽然参考文件2中方法能够在一定程度上探测到回波中遗失组分，但是，由于模型初始化时采用拟合残差最小为最优条件，因此引入的误差在波形分解前就带来对后续遗失组分探测的影响，从而导致组分可能不能被完整探测，影响波形分解的精度。

发明内容

(一)要解决的技术问题

鉴于上述技术问题，本发明提供了一种兼顾模型组分数的初始化准确度和波形分解的完整性的全波形激光雷达数据波形分解方法。

(二)技术方案

根据本发明的一个方面，提供了一种全波形激光雷达数据波形分解方法。该全波形激光雷达数据波形分解方法包括：

步骤A，对后向散射回波波形数据的组分数N进行初始化估计，得到组分数N的初始值令其中，该后向散射回波波形数据包括采样点数S、各采样点的采样时间x_i和强度值y_i，其中，x_i＝t₀+iΔt，i＝1、2、......、S，t₀为采样开始时间，Δt为采样间隔；

步骤B，基于后向散射回波波形数据以及组分数N，把S个采样点的强度值y_i作为f(x)代入如下广义高斯混合模型，确定广义高斯混合模型中各组分的初始参数，进行广义高斯混合模型的初始化：

$f\left(x\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{f}_j\left(x\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_j\exp (-\frac{{|x-{\mathrm{\μ}}_j|}^{{\mathrm{\α}}_j^2}}{{2\mathrm{\σ}}_j^2})$

其中，f(x)为采用广义高斯混合模型拟合的后向散射回波波形数据；f_j(x)为第j个组分的波形数据，其中x为采样点的采样时间，组分j的初始参数包括：振幅A_j，组分位置μ_j，波形宽度σ_j(1σ)，波形形状因子α_j，1≤j≤N。

步骤C，利用初始化的广义高斯混合模型，基于EM算法自迭代对广义高斯混合模型中的各个参数进行更新，完成后向散射回波波形数据的混合高斯分解，得到分解后的后向散射回波波形数据；以及

步骤D，判断分解后的后向散射回波波形数据是否存在遗失组分，如果不存在遗失组分，则波形分解完成，如果存在遗失组分，估计新的组分数N^*，用组分数N^*替换原有广义高斯混合模型中的N，重新执行步骤B～D。

(三)有益效果

从上述技术方案可以看出，本发明全波形激光雷达数据波形分解方法，与一般的波形分解方法相比，兼顾模型组分数的初始化准确度和波形分解的完整性，具有以下有益效果：

(1)利用AIC进行组分数估计，可快速获取后向散射回波波形中的组分数信息，算法简单，探测效率高；

(2)通过将估计后的组分数和组分的振幅、位置、波形宽度等模型参数导入到LM算法分析波形分解结果，可进一步探测到遗失组分，最大化的保证波形分解的完整性。

附图说明

图1为本发明实施例全波形激光雷达数据波形分解方法的流程图；

图2为分别采用本发明方法与参考文献1方法对测试数据进行波形分解后得到的波形数据f(x)的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。需要说明的是，在附图或说明书描述中，相似或相同的部分都使用相同的图号。附图中未绘示或描述的实现方式，为所属技术领域中普通技术人员所知的形式。另外，虽然本文可提供包含特定值的参数的示范，但应了解，参数无需确切等于相应的值，而是可在可接受的误差容限或设计约束内近似于相应的值。

本发明针对全波形激光雷达数据波形分解时存在的遗失组分问题，提出联合AIC和LM的全波形激光雷达数据波形分解方法，自迭代判断遗失组分数目及位置，从而可以有效提高波形分解精度。

在本发明的一个示例性实施例中，提供了一种全波形激光雷达数据波形分解方法。如图1所示，本实施例全波形激光雷达数据波形分解方法包括：

步骤A，对后向散射回波波形数据的组分数N进行初始化估计，得到组分数N的初始值令其中，该后向散射回波波形数据包括采样点数S、各采样点的采样时间x_i和强度值y_i，其中，x_i＝t₀+iΔt，i＝1、2、......、S，t₀为采样开始时间，Δt为采样间隔；

本步骤中，后向散射回波波形数据事先应当进行去噪平滑预处理。本步骤中，采用AIC方法对后向散射回波波形数据的组分数N进行初始化估计，得到组分数N的初始值进一步包括：

子步骤A1，对S个采样点强度值组成的采样强度序列采用一阶导数寻找极大值，令N₀＝M，其中M为极大值的个数；

子步骤A2，基于AIC方法估计后向散射回波波形数据组分数N进行初始化估计，得到组分数的初始值具体来讲，根据以下公式计算组分数初始值

$\hat{N}=\arg \min \log \left(v_N\left(N\right)\right)+2N/S---\left(1\right)$

$v_N\left(N\right)=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{(y_i-{\hat{y}}_i)}^2---\left(2\right)$

${\hat{y}}_i=\underset{j=1}{\overset{N_o}{\mathrm{\Σ}}}8*\exp (-\frac{{|x_i-{\mathrm{\μ}}_j|}^2}{50})---\left(3\right)$

其中，arg min log(v_N(N))为使v_N(N)取最小值时的组分数N，μ_j为采样强度序列中第j个极大值对应的采样时间。

步骤B，基于后向散射回波波形数据以及组分数N，把S个采样点的强度值y_i作为f(x)代入如下广义高斯混合模型，结合非线性最小二乘法NLS确定广义高斯混合模型中各组分的初始参数，进行广义高斯混合模型的初始化：

$f\left(x\right)={\mathrm{\Σ}}_{j=1}^N{f}_j\left(x\right)={\mathrm{\Σ}}_{j=1}^N{A}_j\exp (-\frac{{|x-{\mathrm{\μ}}_j|}^{{\mathrm{\α}}_j^2}}{{2\mathrm{\σ}}_j^2})---\left(4\right)$

其中，f(x)为采用广义高斯混合模型拟合的后向散射回波波形数据；f_j(x)为第j个组分的波形数据，其中x为采样点的采样时间，组分j的初始参数包括：振幅A_j，组分位置μ_j，波形宽度σ_j(1σ)，波形形状因子α_j，1≤j≤N。

步骤C，利用初始化的广义高斯混合模型，基于EM算法自迭代对广义高斯混合模型中的各个参数进行更新，完成后向散射回波波形数据的混合高斯分解，得到分解后的后向散射回波波形数据；

该步骤C进一步包括：

子步骤C1，利用广义高斯混合模型的初始参数，计算初始似然方程：

$L^{\left(0\right)}=\frac{1}{S}{\mathrm{\Σ}}_{i=1}^S\log \left({{\mathrm{\Σ}}_{j=1}^{N}f}_j\left(x_i\right|A_j^{\left(0\right)},{\mathrm{\μ}}_j^{\left(0\right)},{\mathrm{\σ}}_j^{\left(0\right)},{\mathrm{\α}}_j^{\left(0\right)})\right)---\left(5\right)$

其中，为第0次迭代时广义高斯混合模型中第j个组分的参数A_j，μ_j，σ_j，α_j，即初始模型参数，上标(0)表示迭代次数；为模型参数是的第j个组分的波形数据在采样时间x_i处的值。

子步骤C2，计算第i个采样点属于第j个组分的概率

${\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)}=\frac{f_j\left(x_i\right|A_j^{\left(m\right)},{\mathrm{\μ}}_j^{\left(m\right)},{\mathrm{\σ}}_j^{\left(m\right)},{\mathrm{\α}}_j^{\left(m\right)})}{{\mathrm{\Σ}}_{j=1}^N{f}_j\left(x_i\right|A_j^{\left(m\right)},{\mathrm{\μ}}_j^{\left(m\right)},{\mathrm{\σ}}_j^{\left(m\right)},{\mathrm{\α}}_j^{\left(m\right)})},i=1,....,S,m=\mathrm{0,1}...---\left(6\right)$

其中，为第m次迭代时广义高斯混合模型中第j个组分的参数，m初值为0，为模型参数是的第j个组分的波形数据在采样的时间x_i处理值。

子步骤C3，基于第m次迭代时的参数，计算获得第m+1次迭代时广义高斯混合模型的参数：

$A_j^{(m+1)}=\frac{n_j^{\left(m\right)}}{N}---\left(7\right)$

${\mathrm{\μ}}_j^{(m+1)}=\frac{1}{n_j^{\left(m\right)}}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)}{y}_i---\left(8\right)$

${\mathrm{\σ}}_j^{(m+1)}=\frac{1}{n_j^{\left(m\right)}}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)}{(y_i-{\mathrm{\μ}}_j^{(m+1)})}^2---\left(9\right)$

其中， $n_j^{\left(m\right)}={\mathrm{\Σ}}_{i=1}^S{\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)}.$

子步骤C4，基于第m+1次迭代时广义高斯混合模型的参数，计算第m+1次似然方程

$L^{(m+1)}=\frac{1}{S}{\mathrm{\Σ}}_{i=1}^S\log \left({\mathrm{\Σ}}_{j=1}^N{f}_j\left(x_i\right|A_j^{(m+1)},{\mathrm{\μ}}_j^{(m+1)},{\mathrm{\σ}}_j^{(m+1)},{\mathrm{\α}}_j^{\left(\text{m+1}\right)})\right)---\left(10\right)$

子步骤C5，判断是否满足终止迭代条件，设置阈值ε∈(0，1)，当|L^(m+1)-L^(m)|≤ε时，迭代终止，波形数据的混合高斯分解完毕，得到分解后各组分的波形数据和模型参数，进入步骤D；当|L^(m+1)-L^(m)|＞ε时，m＝m+1，重新进行子步骤C2，直到判断迭代终止。

步骤D，采用LM算法自迭代判断步骤C分解后的后向散射回波波形数据是否存在遗失组分，如果不存在遗失组分，则波形分解完成，如果存在遗失组分，估计新的组分数N^*，用组分数N^*替换原有广义高斯混合模型中的N，重新执行步骤B～D；

该基于LM算法自迭代判断遗失组分的步骤D进一步包括：

子步骤D1，计算非负常数 $\mathrm{\ϵ}\left(\mathrm{\θ}\right)=\frac{1}{S-\mathrm{dim\θ}}{\mathrm{\Σ}}_{i=1}^S{[y_i-f\left(x_i\right|\mathrm{\θ})]}^2;$

其中，θ＝(A₁，μ₁，σ₁，α₁，...，A_j，μ_j，σ_j，α_j，...，A_N，μ_N，σ_N，α_N)，j＝1，...，N为步骤C分解得到的广义高斯混合模型参数组成的参数向量，f(x_i|θ)是参数为θ时后向散射回波波形数据f(x)在采样时间x_i处的值，dimθ为参数向量θ的维数。

子步骤D2，设置阈值ε′∈(0，1)，当ε(θ)＞ε′时，判断存在遗失组分，估计新的组分数 $N*=\arg \min \left\{{\mathrm{\Σ}}_{i=1}^S{[y_i-f\left(x_i\right|\mathrm{\θ})]}^2\right\},$ 用组分数N^*替换原有广义高斯混合模型中的N，重新执行步骤B～D；当ε(θ)≤ε′时，判断不存在遗失组分，波形分解完成，得到N个组分波形数据及广义高斯混合模型参数组成的参数向量，包括各组分的振幅、位置、波形宽度和波形形状因子。

如图2所示为分别采用本发明方法与参考文献1方法对测试数据(后向散射回波波形数据采样点)进行波形分解后得到的波形数据f(x)的示意图。其中，采用本实施例方法得到的波形数据f(x)用实线表示，采用参考文献1方法得到的波形数据f(x)用虚线表示。从图2中可以看出，相对于参考文献1方法，采用本发明方法得到的波形数据f(x)更贴近于原始数据，有利于组分的完整探测，可确切反映原始数据的真实情况，包括组分数、各组分的振幅、组分位置等。

至此，已经结合附图对本实施例进行了详细描述。依据以上描述，本领域技术人员应当对本发明全波形激光雷达数据波形分解方法有了清楚的认识。

此外，上述对各步骤的定义并不仅限于实施方式中提到的各种具体方式，本领域的普通技术人员可对其进行简单地熟知地替换，例如：

(1)步骤A中AIC方法可以用LM方法来代替进行组分数初始化；

(2)步骤B中非线性最小二乘法可以用极大似然估计来代替进行广义高斯模型的初始化；

(3)步骤D中LM算法可以用惩罚最小匹配距离代替判断是否存在遗失组分。

综上所述，本发明提供一种顾及遗失组分的全波形激光雷达数据波形分解方法，与一般的波形分解方法相比，兼顾广义高斯混合模型中组分数的初始化准确度和波形分解的完整性，具有方法简单，探测效率高的优势，可以应用于特征提取、点云精细分类、DEM生成等领域。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种全波形激光雷达数据波形分解方法，其特征在于，包括：

步骤A，对后向散射回波波形数据的组分数N进行初始化估计，得到组分数N的初始值令其中，该后向散射回波波形数据包括采样点数S、各采样点的采样时间x_i和强度值y_i，其中，x_i＝t₀+iΔt，i＝1、2、......、S，t₀为采样开始时间，Δt为采样间隔；

步骤B，基于后向散射回波波形数据以及组分数N，把S个采样点的强度值y_i作为f(x)代入如下广义高斯混合模型，确定广义高斯混合模型中各组分的初始参数，进行广义高斯混合模型的初始化：

$f\left(x\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{f}_j\left(x\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_j\exp (-\frac{{|x-{\mathrm{\μ}}_j|}^{{\mathrm{\α}}_j^2}}{{2\mathrm{\σ}}_j^2})$

其中，f(x)为采用广义高斯混合模型拟合的后向散射回波波形数据；f_j(x)为第j个组分的波形数据，其中x为采样点的采样时间，组分j的初始参数包括：振幅A_j，组分位置μ_j，波形宽度σ_j(1σ)，波形形状因子α_j，1≤j≤N。

步骤C，利用初始化的广义高斯混合模型，基于EM算法自迭代对广义高斯混合模型中的各个参数进行更新，完成后向散射回波波形数据的混合高斯分解，得到分解后的后向散射回波波形数据；以及

步骤D，判断分解后的后向散射回波波形数据是否存在遗失组分，如果不存在遗失组分，则波形分解完成，如果存在遗失组分，估计新的组分数N^*，用组分数N^*替换原有广义高斯混合模型中的N，重新执行步骤B～D。

2.根据权利要求1所述的全波形激光雷达数据波形分解方法，其特征在于，所述步骤A中，采用AIC方法或LM方法对后向散射回波波形数据的组分数N进行初始化估计。

3.根据权利要求2所述的全波形激光雷达数据波形分解方法，其特征在于，采用AIC方法对后向散射回波波形数据的组分数N进行初始化估计，得到组分数N的初始值进一步包括：

子步骤A1，对S个采样点强度值组成的采样强度序列采用一阶导数寻找极大值，其中M为极大值的个数，令N₀＝M；

子步骤A2，基于AIC方法估计后向散射回波波形数据组分数N进行初始化估计，得到组分数的初始值具体来讲，根据以下公式计算组分数初始值

$\hat{N}=\arg \min \log \left(v_N\left(N\right)\right)+2N/S$

$v_N\left(N\right)=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{(y_i-{\hat{y}}_i)}^2$

${\hat{y}}_i=\underset{j=1}{\overset{N_o}{\mathrm{\Σ}}}8*\exp (-\frac{{|x_i-{\mathrm{\μ}}_j|}^2}{50})$

其中，arg min log(v_N(N))为使v_N(N)取最小值时的组分数N，μ_j为采样强度序列中第j个极大值对应的采样时间。

4.根据权利要求1所述的全波形激光雷达数据波形分解方法，其特征在于，所述步骤B中采用非线性最小二乘法或极大似然估计来进行广义高斯混合模型的初始化。

5.根据权利要求1所述的全波形激光雷达数据波形分解方法，其特征在于，所述步骤C进一步包括：

子步骤C1，利用广义高斯混合模型的初始参数，计算初始似然方程：

$L^{\left(0\right)}=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\log \left(\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{f}_j\left(x_i\right|A_j^{\left(0\right)},{\mathrm{\μ}}_j^{\left(0\right)},{\mathrm{\σ}}_j^{\left(0\right)},{\mathrm{\α}}_j^{\left(0\right)})\right)$

其中，为第0次迭代时广义高斯混合模型中第j个组分的参数A_j，μ_j，σ_j，α_j，即初始模型参数，上标(0)表示迭代次数；为模型参数是的第j个组分的波形数据在采样时间x_i处的值；

子步骤C2，计算第i个采样点属于第j个组分的概率

${\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)}=\frac{f_j\left(x_i\right|A_j^{\left(m\right)},{\mathrm{\μ}}_j^{\left(m\right)},{\mathrm{\σ}}_j^{\left(m\right)},{\mathrm{\α}}_j^{\left(m\right)})}{{\mathrm{\Σ}}_{j=1}^N{f}_j\left(x_i\right|A_j^{\left(m\right)},{\mathrm{\μ}}_j^{\left(m\right)},{\mathrm{\σ}}_j^{\left(m\right)},{\mathrm{\α}}_j^{\left(m\right)})},i=1,....,S,m=\mathrm{0,1}...$

其中，为第m次迭代时广义高斯混合模型中第j个组分的参数，m初值为0，为模型参数是的第j个组分的波形数据在采样时间x_i处的值；

子步骤C3，基于第m次迭代时的参数，计算获得第m+1次迭代时广义高斯混合模型的参数：

$A_j^{(m+1)}=\frac{n_j^{\left(m\right)}}{N}$

${\mathrm{\μ}}_j^{(m+1)}=\frac{1}{n_j^{\left(m\right)}}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)}{y}_i$

${\mathrm{\σ}}_j^{(m+1)}=\frac{1}{n_j^{\left(m\right)}}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)}{(y_i-{\mathrm{\μ}}_j^{(m+1)})}^2$

其中， $n_j^{\left(m\right)}={\mathrm{\Σ}}_{i=1}^S{\mathrm{\γ}}_{\mathrm{ij}}^{\left(m\right)};$

子步骤C4，基于第m+1次迭代时广义高斯混合模型的参数，计算第m+1次似然方程

$L^{(m+1)}=\frac{1}{S}\underset{i=1}{\overset{S}{\mathrm{\Σ}}}\log \left(\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{f}_j\left(x_i\right|A_j^{(m+1)},{\mathrm{\μ}}_j^{(m+1)},{\mathrm{\σ}}_j^{(m+1)},{\mathrm{\α}}_j^{\left(\text{m+1}\right)})\right)$

子步骤C5，判断是否满足终止迭代条件，设置阈值ε∈(0，1)，当|L(m+1)-L(m)|≤ε时，迭代终止，波形数据的混合高斯分解完毕，得到分解后各组分的波形数据和模型参数，进入步骤D；当|L^(m+1)-L^(m)|＞ε时，m＝m+1，重新进行子步骤C2，直到判断迭代终止。

6.根据权利要求1所述的全波形激光雷达数据波形分解方法，其特征在于，所述步骤D中采用LM算法或惩罚最小匹配距离法判断分解后的后向散射回波波形数据是否存在遗失组分。

7.根据权利要求6所述的全波形激光雷达数据波形分解方法，其特征在于，采用LM算法判断分解后的后向散射回波波形数据是否存在遗失组分具体包括：

子步骤D1，计算非负常数 $\mathrm{\ϵ}\left(\mathrm{\θ}\right)=\frac{1}{S-\mathrm{dim\θ}}{\mathrm{\Σ}}_{i=1}^S{[y_i-f\left(x_i\right|\mathrm{\θ})]}^2;$

其中，θ＝(A₁，μ₁，σ₁，α₁，...，A_j，μ_j，σ_j，α_j，...，A_N，μ_N，σ_N，α_N)，j＝1，...，N为步骤C分解得到的广义高斯混合模型参数组成的参数向量，f(x_i|θ)是参数为θ时后向散射回波波形数据f(x)在采样时间x_i处的值，dimθ为参数向量θ的维数；

子步骤D2，设置阈值ε′∈(0，1)，当ε(θ)＞ε′时，判断存在遗失组分，估计新的组分数 $N*=\arg \min \left\{{\mathrm{\Σ}}_{i=1}^S{[y_i-f\left(x_i\right|\mathrm{\θ})]}^2\right\},$ 用组分数N^*替换原有广义高斯混合模型中的N，重新执行步骤B～D；当ε(θ)≤ε′时，判断不存在遗失组分，波形分解完成，得到N个组分波形数据及广义高斯混合模型参数组成的参数向量，包括各组分的振幅、位置、波形宽度和波形形状因子。