《测绘学报》
构建与学术的桥梁 拉近与权威的距离
机载激光测深波形分解中LM与EM参数优化方法比较
郭锴1,2, 刘焱雄1, 徐文学1
1. 自然资源部第一海洋研究所, 山东 青岛 266061;
2. 深圳大学海岸带地理环境监测自然资源部重点实验室, 广东 深圳 518060;
3. 河海大学地球科学与工程学院, 江苏 南京 210098;
4. 国网经济技术研究院有限公司, 北京 102209;
5. 国网山东省电力公司经济技术研究院, 山东 济南 250022
收稿日期：2018-05-23；修回日期：2019-04-14
基金项目：国家自然科学基金（41871381；41401537）；国家电网公司科技项目（JYYKJXM[2017]003）；中央级公益性科研院所基本科研业务费专项资金（2015P13）
第一作者简介：郭锴(1987-), 男, 博士, 研究方向为机载激光测深数据处理及应用、海陆一体化测量。E-mail:guokai4545@163.com
通信作者：徐文学, E-mail：xuwx@fio.org.cn
摘要：受作业环境与扫描条件的影响，机载激光测深全波形数据处理所得的参数分量初值精度通常较低，需采用优化方法在此基础上进一步调整分量参数以提高波形分解精度。本文代表性地选取了非线性阻尼最小二乘方法（Levenberg-Marquardt，LM）和期望最大化方法（expectation-maximization，EM）两种不同的参数优化方案，针对机载激光测深全波形数据的波形分解参数优化问题，结合实测和模拟波形数据对两种优化方法在相同初值条件下的波形模拟精度、扫描条件及其响应，以及水深估计偏差等方面的特征进行了详细的对比分析，着重探讨了两种参数优化方法所得结果的准确度与稳定性，并对其主要技术特征与效果差异进行了总结。
关键词：机载激光测深 水下地形 海洋测绘 非线性阻尼最小二乘方法 期望最大化方法
Comparison of LM and EM parameter optimization methods for airborne laser bathymetric full-waveform decomposition
GUO Kai1,2, LIU Yanxiong1, XU Wenxue1, HE Xiufeng3, ZHANG Ziyin4, LIU Haitao5
1. The First Institute of Oceanography, Ministry of Natural Resources, Qingdao 266061 China;
2. Key Laboratory for Geo-Environmental Monitoring of Coastal Zone of the Ministry of Natural Resources, Shenzhen University, Shenzhen 518060, China;
3. School of Earth Science and Engineering, Hohai University, Nanjing 210098, China;
4. State Grid Economic and Technological Research Institute Co., Ltd., Beijing 102209, China;
5. Economic and Technology Research Institute of State Grid Shandong Electric Power Company, Jinan 250022, China
Foundation support: The National Natural Science Foundation of China (Nos. 41871381; 41401537); The Technology Project of State Grid Corporation of China (No. JYYKJXM [2017] 003); The Basic Scientific Fund for National Public Research Institutes of China (No. 2015P13)
First author: GUO Kai (1987—), male, PhD, majors in data processing and application of airborne laser bathymetry, integrated surveying and mapping technology of ocean and land.E-mail:guokai4545@163.com.
Corresponding author: XU Wenxue, E-mail : xuwx@fio.org.cn.
Abstract: As operating environment and scanning conditions would affect the accuracy of the parameters of initial waveform decomposition, it is necessary to use some optimization method to adjust the initial parameters for making the results more precise. With the measured and simulated data, the strategies of Levenberg-Marquardt (LM) and expectation-maximization (EM) have been chosen typically in this paper. It is studied that the features of the waveform fitting precision, the scanning condition response and the depth of the deviation in the same initial parameters condition. The accuracy and stability of the data processing results has been discussed in detail, and the main technical features and the difference of the results of the two methods has been summarized in this paper.
Key words: airborne laser bathymetry underwater topographic survey hydrographic surveying Levenberg-Marquardt expectation-maximization
机载激光测深(airborne laser bathymetry, ALB)技术是一种准确、高效的水深探测技术，具有海陆一体化、全天时以及高分辨率等特性，在水下地形测量、岸线提取、岸滩地貌描绘以及环境动态监测等方面均有广泛应用[1-2]。回波信号的全波形记录是机载激光测深技术的基础[3-4]。扫描获得的全波形数据可通过回波能量在时序上的分布特征及其度量参数反映目标物的空间与物理特征[5-7]。目标物的位置及特征信息可通过波形分解获取。波形分解方法的实质是将原始信号理解为多个波形的叠加[8]，并根据激光发射信号和雷达方程对全波形数据的各部分进行有效分解，从而得到各回波分量准确的时间位置及相关参数。
目前，按照分解后波形分量的函数类型主要分为高斯分解与三角函数、四边形函数等非高斯分解[8-12]。文献[8, 13-14]认为系统的激光发射脉冲与目标的反射回波能量均具有类似高斯分布的相关特征；文献[15]在对Wa-LiD测深回波模拟工具进行研究时认为，水体回波波形可表示为各深度位置反射分量的叠加。在此基础上，文献[16-17]将其运用至相关分析试验中，并取得了良好的效果。因此，采用高斯函数相互叠加的方式可有效实现激光回波波形的模拟。采用波形分解的缺点在于建模过程复杂，其效果容易受到初值和估计方法的影响[5, 8]。因此，对波形进行首次分解时，回波信号噪声与不良的波形分解过程，往往使得其结果与真实波形数据在整体上产生明显差异。为减小波形分解结果与实际回波信号的不符值，提高波形分解精度，可采用对参数进行多次优化的方法直至其结果趋于稳定。
参数优化是在获得参数初值的基础上，按照波形拟合最优化标准对模型参数进行调整，使其拟合效果达到最佳的过程[8]。受传播路径的复杂性影响，在单一波形分解中，回波分量的参数优化是保证测量精度的关键，也是波形分解的重要环节。性能良好的优化方法对提高探测精度，准确反映目标的空间与物理特征具有重要意义[18]。目前，基于非线性最小二乘的分层高斯函数波形拟合[8]、高斯-牛顿法[7]、期望最大化(expectation-maximization, EM)算法[19]，以及改进EM波形分析技术[20]等方法均受到了广泛认可。
参数优化方法之间的主要区别在于对回波波形的处理思路不同。典型的参数优化方法，如非线性阻尼最小二乘法(Levenberg-Marquardt, LM)[8]，认为全波形数据直观上是一个“时间-强度”函数，通过分解该函数得到各反射分量的时空形态从而实现接收波形的分解[21-23]。另一种方法认为，全波形数据是某一时刻接收到来自目标物回波能量的随机函数。当概率分布的值越大时，说明在对应时间位置处回波为地物反射的可能性越大[19, 24]，如期望最大化法[25-26]以及可逆跳变的马尔可夫蒙特卡洛方法(reversible jump Markov chain Monte Carlo, RJMCMC)[18-19, 27-28]等。目前，波形分解研究及应用多集中在陆地机载激光测量数据处理中，对测深全波形数据的波形分解研究相对不足。为了解不同参数优化方法在处理机载激光测深全波形数据分解过程和拟合效果方面的异同，本文通过采用Wa-LiD(water LiDAR simulator)工具在排除环境影响的情况下对不同系统条件的全波形数据进行模拟[15-16, 29]，并在确定高斯分量初值的基础上，代表性地选取了LM方法和EM方法中的参数优化部分进行分别处理。通过比较两种方法在处理过程与分解效果的主要差异，分析不同类型参数优化方法的技术特点及其与扫描条件间的响应关系。此外，根据实测区域范围内波形数据的处理结果，探讨不同参数优化方法对水深测量结果的影响，并进一步考查参数优化方法在水深估算方面的准确度与适用性，为全波形数据处理方法的深入研究提供参考。
1 回波波形的数学模拟
按照机载激光的测深原理，激光器发射脉冲信号经水面和水底反射后被接收器接收，以水面与水底反射发生的时间差作为计算激光在水中传播路径长度的依据，并通过三角函数归算至垂直方向。采用全波形技术的扫描测深系统可较为完整地记录回波能量在时间上的分布特征，为估算目标的空间位置并同时获取更为丰富的反射信息，建立能量变化与反射特征间的函数关系，并采用数学模型对回波波形进行模拟具有重要意义。
根据研究[8]，发射激光和回波能量之间的关系可表示为
(1)
式中，i=1, 2, …,n；Pr(ti)为信号的接收能量；Ptj(ti)为激光发射信号的强度；n表示激光脉冲在大于距离分辨率的目标间传播反射过程中的系统的采样个数；k为反射发生的次数；Dr为激光发射器的孔径(直径)；Ri指激光器到目标的距离；ωi为发射电波波束宽度；Φj(ti)表示接收器脉冲函数；Aj(ti)表示微分后向散射截面在微小距离dR条件下的响应函数，该函数显示了散射体的散射特征；依据激光测距的基本原理，Ri与传播时间ti之间的关系应为
，其中vg为光在介质中的传播速度；*为卷积运算。
假设雷达散射体的散射特征可以用高斯函数进行描述，则单个高斯函数可以表示为如下形式
(2)
式中，
(3)
式中，
(4)
(5)
通过以上分析可知，激光的接收波形可看作k个高斯函数的叠加。如果可以准确估计各个高斯回波分量的参数Aj,μj,σj, 即可实现回波波形的模拟。因此，从波形数据处理的本质上讲，回波波形的模拟就是高斯分量的参数估计及其最优化过程。
2 波形分解预处理与初值估计
采用高斯函数对全波形测深数据进行分解之前，需对原始波形数据进行相应的预处理，包括波形有效部分的选取与噪声滤波以减少原始波形中的冗余数据并消除信号中随机噪声的干扰[30]。之后，针对处理后的波形数据估算高斯分量的初值参数，通常初值分量的拟合波形与真实的回波波形之间存在较为明显的差异，需要通过参数优化方法调整相关参数的具体数值，直至满足波形拟合的最优化条件。图 1为采用波形分解方法处理全波形数据的基本流程。
图 1 全波形数据分解基本流程Fig. 1 The basic process of full-waveform data decomposition
图选项
2.1 有效波形选取
ALB系统通过对所接收的回波信号按照固定频率进行采样来实现数据采集，实际作业过程中获得的全波形数据往往存在大量冗余。另外，有效波形以外部分所包含的噪声干扰对反射时间位置的准确判断也造成了不利影响。因此，在对全波形数据进行处理之前，为减少数据冗余，锁定回波波形的所在时间区间是提高算法效率，保证结果质量所进行的必要步骤。通常以某一阈值作为信号有效部分截取的限定依据[26]。虽然该方法操作简单，易于实现，但在选取效果上存在信号中断且无法顾及噪声的整体性影响等问题，容易产生回波有效部分的误判。
为了尽量保证有效波形的完整性，使有效部分的首尾准确截断在局部最小值处，本文从原始信号的首端和尾端分别向有效部分方向的极值进行滑动取值，形成极值序列。在不断更新极值序列的过程中，以3倍中误差为置信区间，即当新序列的标准差超过上一步序列的3倍中误差时，认为已达到反射信号的有效部分，同时停止更新序列，并认为前一个加入的极小值之后的信号是激光反射回波的有效部分[31]。
2.2 数据滤波
由于激光在发射、传播以及接收过程中受到来自传感器或大气散射等原因造成的不确定性影响，机载激光所接收的回波波形中往往存在大量噪声，不同程度的噪声使得波形对反射信息的有效表达产生不良影响，因此，在对波形进行分解之前需首先对信号中存在的噪声进行抑制。通常而言，一维数据的滤波方法均可用于全波形数据的去噪。本文根据文献[32-34]中保留信号矩的去噪方法，采用离散滤波系统对回波信号进行滤波，在削弱波形中随机噪声干扰的同时，尽量保留信号中回波信息的有效成分。相关参数的计算标准在于使滤波前原始信号矩和滤波后的信号矩相同。通常随着滤波阶数的增大，波形整体更为平缓，说明滤波有效抑制了回波信号中的高频信息，在波形滤波处理过程中应根据实测数据选择适当的滤波器阶数，以保证得到良好的抑噪效果。当滤波阶数为12阶时对本文中波形数据的去噪效果达到最佳，因此，本文波形试验均采用这一滤波阶数进行处理。
2.3 参数初值估计
波形分解中，参数的初值估计具有重要的作用，特别是对于参数优化而言，多数方法需要采用迭代计算。因此，采用较为准确的参数初值不仅可以加快系统的收敛过程，同时也可在一定程度上避免算法陷入局部最优解。常规高斯分解首先通过峰值探测估计波形中高斯分量的位置和数量，且在分解过程主要遵循强度顺序，即按照波形中波峰的强度大小依次进行分解。但这样的方法容易因反射目标间距离过近而引起峰值偏移现象[33-34]。因此本文采用以全波形数据接收时间顺序为基础的高斯半波长递进分解方法进行初值估计[31]与反射分量的选取(图 2(b))。该方法可有效避免由于回波分量叠加而造成的反射峰值位置偏移。
图 2 高斯分解的初值估计(数字代表分解顺序)Fig. 2 Initial value estimation of Gaussian decomposition (the number represent the decomposition order of the waveform components)
图选项
3 参数优化
一般情况下，通过初值估计直接得到高精度波形分量的可能性很小。通常需要在初值估计的基础上，按照一定的最优化标准对分量的参数进行多次调整，直至满足波形的拟合精度。参数优化过程是确保波形分解有效性以及激光特征提取精度的重要步骤。目前，普遍采用的参数优化方法主要包括两类：最小二乘方法和概率分布方法。前者将全波形数据直观地表述为一种时间域函数，其中，最为普遍的是非线性阻尼最小二乘法，后者将波形数据理解为一段时间内传感器所接收回波能量的概率密度，强度较大的时刻对应接收回波能量概率较大，具有代表性的是期望最大化方法。本文主要选取了这两种方法进行研究，对二者在计算思路和方法效果方面进行了对比分析，从而为机载激光测深波形分解方法研究方面提供参考。
3.1 EM方法
设在某时段内，时间采样样本为{ti}，其中i=1, 2, …,n。接收反射能量的概率分布可表示为包含k个隐含分布的混合函数，隐含的分布函数表示为Zj(ti)，j=1, 2, …,k。则在该时间范围内，反射所发生时间位置在第j次分布下的概率即样本来自于Zj(ti)分布的概率可表示为
(6)
设ti时接收反射能量概率分布为y(ti)
(7)
两边取对数得
(8)
这里假设P{Zj(ti)}服从高斯分布，其参数为Aj,μj,σj，按照极大似然估计(maximum likelihood estimation, MLE)[35]，迭代中各参数可表示如下
(9)
(10)
(11)
(12)
式中，fj∈N(μj,σj2)，根据以上公式进行循环迭代，直至收敛，至此得到最终优化后的参数结果。
3.2 LM方法
LM方法是传统高斯-牛顿法的改进，其基本思想是采用泰勒级数的前二阶作为函数的近似进行迭代计算。不同之处在于，LM方法认为海塞矩阵为一个对称矩阵，在求逆过程中容易出现矩阵奇异而造成方程无唯一解的情况。因此，可根据数据具体情况在海塞矩阵对角线上加一个阻尼系数，保证原函数存在唯一解。高斯函数为
(13)
式中，i=1, 2, …,n，j=1, 2, …,k，有关参数包括θj=[Ajμjσj]。对其分别求偏导数可得
(14)
(15)
(16)
由于采用高斯函数作为拟合函数，在进行迭代计算时，高斯函数所需的迭代初值为θ0：(A0,μ0,σ0)，迭代过程中Jacobian矩阵为
(17)
Levenberg-Marquart步
(18)
ΔθN可作为参数向量改正序列。实现以上迭代过程需要以下数据：[t,y(t)]为原始数据在二维空间中的位置，λ为阻尼系数，J为Jacobian矩阵，N为迭代次数，如N=0表示高斯函数参数的迭代初值，其中，包括振幅Aj(0)、波峰位置μj(0)，波形展宽(标准差)σj(0)。v1为放大系数，一般取(1 <v1< 10)。v2为收缩系数，一般取(0.1 <v2< 0.9)。LM算法的计算过程如下：
(1) 确定待估参数向量初值θ(0)，以及允许误差ε，可作为停止迭代的阈值条件。
(2) 在θ(N)处计算[y(t)－f(t|θ(N))，原始波形和拟合波形的总残差为
(19)
(3) 由式(18)迭代计算各个位置的[θ(N)－θ(N－1)]，若
S(θ(N))≤S(θ(N－1))，且|θ(N)－θ(N－1)| <ε则接受θ(N)，整个计算结束；
S(θ(N))≤S(θ(N－1))，但|θ(N)－θ(N－1)| >ε，则以μjNv2代替原参数μjN，更新θ(N)，转到第(2)步再次迭代；S(θ(N))>S(θ(N－1))，则取μjNv1代替原参数μjN，更新θ(N)，转到第(2)步重新迭代。
4 试验与分析4.1 单一波形处理效果比较
为比较不同方法在全波形数据分解过程中的参数优化效果，本文分别对实测与模拟两种波形数据进行了有效处理。前者主要从直观上反映LM方法与EM方法在拟合效果上的区别，后者通过模拟均匀介质条件下的全波形数据，对波形拟合的准确度以及所选方法对水深估计的影响进行量化研究。
4.1.1 实测波形数据分析
由于对波形分解的理解方式不同，EM和LM方法采用了完全不同的参数优化思路，即使采用相同的初值，二者的处理结果也应有所不同。本文以Aquarius系统在中国南海某海岛周边水域的实测全波形数据为对象进行处理分析。所选区域底质以珊瑚礁等生物沉积为主且较为稳定，观测当天天气晴朗，海况良好。对单一实测波形的相关分析中需首先排除参数初值对优化方法的影响，为此本文在参数初值的确定上均采用基于时间顺序的高斯半波长递进分解方法，以2.3节中所得结果作为波形分解中参数优化前的分量初值。
图 3显示了经EM和LM两种参数优化方法处理后波形的拟合效果，可以看出两种方法均对参数初值所对应的波形分量有所调整。二者所得反射分量的时间位置较为接近，说明其最终得到的水深估计值差异性较小，但从振幅强度角度看EM方法处理得到的反射分量较LM方法处理后的分量强度值略高。表 1为经过参数优化处理后所得各高斯分量的具体参数值。
图 3 两种波形参数优化方法处理效果Fig. 3 Results processed by the two optimization methods of the waveform parameters
图选项
表 1 两种参数优化方法处理后各分量参数(波形见图 3)Tab. 1 The parameters of each component are processed by two parameter optimization methods (waveform:Fig. 3)
| 分量 | 初值参数 | LM优化后参数 | EM优化后参数 | ||||||||
| A | μ | σ | A | μ | σ | A | μ | σ | |||
| 1 | 109.100 | 25.159 | 3.417 | 112.42 | 25.398 | 3.725 | 108.510 | 25.219 | 3.726 | ||
| 2 | 60.542 | 33.483 | 3.266 | 50.143 | 32.840 | 2.528 | 60.281 | 33.633 | 3.326 | ||
| 3 | 29.684 | 39.314 | 1.928 | 31.158 | 37.565 | 2.938 | 25.654 | 39.615 | 2.271 | ||
| 4 | 33.378 | 44.949 | 2.799 | 32.048 | 43.491 | 3.682 | 34.356 | 44.871 | 2.990 | ||
| 5 | 16.123 | 50.030 | 2.060 | 13.725 | 48.328 | 2.434 | 14.932 | 49.978 | 2.482 | ||
| 6 | 18.467 | 56.084 | 2.569 | 21.231 | 56.148 | 4.420 | 21.075 | 56.873 | 3.514 | ||
| 7 | 283.030 | 71.697 | 3.293 | 271.270 | 71.305 | 3.008 | 288.54 | 71.759 | 3.250 | ||
| 8 | 93.075 | 77.914 | 2.329 | 101.500 | 76.809 | 2.759 | 94.153 | 77.955 | 2.373 | ||
| 9 | 26.639 | 82.346 | 1.893 | 35.376 | 80.939 | 3.940 | 29.428 | 83.157 | 2.584 | ||
| 与原始波形关联度
LM优化中r0LM=0.978；EM优化中r0EM=0.974 | |||||||||||
表选项
对比单一波形处理效果可以看出，LM方法和EM方法对相同初值分量进行处理后产生了不同的优化效果。就优化后波形与原始波形的接近程度而言，LM参数优化后的波形与原始波形的关联度为0.978，略高于EM优化后的0.974。Hausdorff距离能够度量两个波形间的最大不匹配程度，依据文献[36-37]中双向Hausdorff距离，分别对两种优化方法与原波形之间的Hausdorff距离进行计算。LM优化后所得的参数包括θj：(Aj,μj,σj)，拟合波形与原始波形的Hausdorff距离为1.969，对应EM优化后的结果与原始波形的Hausdorff距离为3.900。以上结果说明，从形态角度看LM方法与EM方法均能较好地反映真实的波形特征，即LM方法的波形拟合准确程度优于EM方法，主要原因在于LM参数优化方法的迭代终止条件主要以原始波形的形态差异为标准，波形的拟合效果较好。
以下随机选取了350个实测波形数据进行高斯分解，分别采用LM和EM两种参数优化方法进行处理，并按照上述方法计算拟合后波形与原波形之间的Hausdorff距离和关联度。为便于比较，此处计算二者所得结果之差，以其差值衡量LM和EM参数优化方法
式中，dEM(i)为采用EM参数优化方法对第i个波形进行处理后所得结果与原波形之间的Hausdorff距离；dLM(i)为采用LM方法所得波形的Hausdorff距离；同理，rEM(i)和rLM(i)分别为EM方法和LM方法所得波形与原波形的关联度。图 4为Hausdorff距离之差和关联度差值的直方图。
图 4 LM与EM处理后波形拟合效果比较Fig. 4 Comparison of the fitting effect after LM and EM treatment
图选项
图 4显示了350个实测波形数据经两种参数优化方法处理后，在波形拟合精度方面所表现出的特征。总体上看，差值较小的部分主要集中在零值附近，偏离零点较远的数据相对较少。说明两种方法处理后其波形的拟合效果在总体上是相近的。参考(图 4(a))以及具体数值后可以看出，采用EM方法所获得的拟合结果与原波形的Hausdorff距离明显大于采用LM方法处理后的波形。该试验中dEM(i)>dLM(i)的波形数共249个，占试验数据总数的71.1%，说明采用LM方法所得结果在叠加后的形态更接近于原波形。此外，通过关联度分析也得到了类似结论(如图 4(b))，其中rEM(i) <rLM(i)的波形数为257个，占试验数据总数的73.4%。说明多数经LM参数优化后的波形与原波形的关联度大于采用EM方法后所得结果，这与Hausdorff距离分析结果是一致的。
4.1.2 模拟数据分析
机载激光测深波形数据模拟工具(Wa-LiD)可通过建立波形发射、传播、反射等过程中激光脉冲的能量变化模型，经过卷积运算模拟出特定系统与环境条件下的机载激光测深系统接收的波形数据[15-16]。该算法避免了实际测量中难以量化的不确定性影响，模拟了均匀介质条件下的波形数据，其结果符合机载激光测深的基本原理，并具有较好的模拟精度[15]。
以下通过采用Wa-LiD波形模拟工具[15]，调节系统参数与水体环境参数以模拟不同条件下的接收波形(图 5)。分别采用LM和EM两种参数优化方法对相同的模拟波形进行处理并分析结果，从而得出两种方法在不同探测条件下的量化效果。图 6表现了拟合波形与原始波形间拟合误差的均方根，均方根的计算方法如下
(22)
图 5 Wa-LiD模拟波形Fig. 5 The schematic diagram of Wa-LiD simulation effect
图选项
图 6 LM和EM参数优化方法的波形拟合误差均方差Fig. 6 The mean square error of waveform fitting error for LM and EM parameter optimization methods
图选项
式中，
通过对比两种参数优化方法在波形拟合方面的效果可以看出，两种拟合方式在不同水深、航高以及扫描角条件下，对单个波形的拟合效果存在一定差别，但该差别并不存在趋势性，且与不同的测量条件相比没有明显的相关性。说明在所选测量条件范围内，相关条件变化对参数优化效果没有明显的影响。另外，由于最小二乘方法主要以波形拟合后形态与原波形之间的差异为参考标准，因此LM方法与EM方法相比其对噪声影响下的波形质量响应较为敏感。从图 7可以看出，LM方法处理后拟合误差的均方差存在较为明显的跳变值，说明LM方法在结果的稳定度方面略逊于EM方法。
图 7 LM和EM参数优化方法的水深估计误差Fig. 7 The mean square error of depth estimation error for LM and EM parameter optimization methods
图选项
模拟数据的分析反映了两种参数优化方法不同的水深估计结果，基本呈随机状态，说明参数优化方法在水深值估算方面与以上3个测量条件之间没有明显的相关性。
图 7以及表 2中的统计结果显示，在传播介质均匀的理想状态下，LM参数优化方法与EM参数优化方法的水深估计精度相当，其估计误差较差均在1 cm范围之内，且呈现出一定的随机性。以上结果表明，所测水深在1~12 m，最大扫描角在-23~ 23°以及航高为200~300 m的范围内，采用全波形数据进行分解处理并获取水深估值时，尽管两种参数优化方法对处理后结果有一定影响，但二者的处理精度相当。其中LM参数优化方法对数据处理后水深估计值的不符值较EM方法处理后结果更具波动性，且表 2中水深估计误差的标准差LM处理方法略高于EM方法，说明EM参数优化方法处理后的结果更为稳定。
表 2 两种参数优化方法水深估计误差统计结果比较Tab. 2 Comparison of statistical results of depth estimation errors between two parameter optimization methods
| m | |||||||||||||||||||||||||||||
| 优化方法统计量 | LM | EM | |||||||||||||||||||||||||||
| 均值 | 标准差 | 均值 | 标准差 | ||||||||||||||||||||||||||
| 深度 | 0.020 | 0.009 | 0.010 | 0.006 | |||||||||||||||||||||||||
| 扫描角 | 0.018 | 0.012 | 0.012 | 0.007 | |||||||||||||||||||||||||
| 航高 | 0.021 | 0.011 | 0.019 | 0.009 | |||||||||||||||||||||||||
表选项
4.2 区域测深数据分析
深水区地形数据可结合单波束、多波束等声学设备获取，浅水区域由于水浅、礁多，一直是海洋测绘中的困难区域，采用机载激光测深技术可有效获取该区域的数据。但浅水区域回波数据存在波形叠加、目标反射回波难以分辨等问题是机载激光测深数据处理的难点。为了验证机载激光测深系统在浅水区域的性能和在中国海域的适应性，自然资源部第一海洋研究所与Optech公司合作，于2012年12月底采用Aquarius测深激光扫描仪对中国南海部分海岛及其周边水域进行了测深试验。Aquarius是Optech公司在Gemini的基础上于2011年推出的浅水型单波段机载激光测深系统，其设计最大测深15 m，采用直线扫描、扫描角宽约为±20[38]。本文采用的全波形数据来自于两个不同底质类型的试验区，其中试验区A位于海岛西北沙嘴(图 8)，该区域底质主要为砂质，水下地形随季节变化而改变，但整体较为平坦，水深最大处约19 m；试验区B位于南海某海岛北部，其底质主要为珊瑚礁与珊瑚礁碎块、碎屑，水深地形变化较为复杂。扫描后覆盖区域离岸最远可达1240 m，试验区内水体透明度约8 m，最大水深约20 m，飞行试验均选择天气晴朗，海况平稳的条件下进行。为削弱系统扫描角与海面波动间的相互作用对测深结果产生的影响，选择图 8和图 9中沿飞行方向，同时基本贯穿测区的地形剖面数据作为实测数据进行对比分析。按照全波形数据处理过程，分别采用LM和EM参数优化方法处理所获得的波形分量，通过与Aquarius随机软件处理后水深数据的结果比较，对两种方法在区域水深估计的影响效果进行分析与评价。
图 8 试验区域A及截取断面位置Fig. 8 The experimental areaAand the position of test profiles
图选项
图 9 试验区域B及截取断面位置Fig. 9 The experimental areaBand the position of test profile
图选项
首先，从两种方法处理后的水深结果互差可以直观看出(图 10)，除有少量差值明显较大的水深点外，大部分的差值均集中在±0.1 m的范围内，即总体上采用两种不同参数优化方法所造成的差异较小且十分有限。另外，如果将波形分解处理后得到的水深值与对应的IWD处理后的水深值进行直接比较，利用二者之差可直接反映采用特定优化方法处理后的结果在水深估计方面的响应
图 10 LM参数优化方法与EM参数优化方法水深估计较差Fig. 10 Depth estimation of LM and EM optimization methods
图选项
式中，Δ为水深估计的误差；|Δ|为水深估计的绝对误差；dIWD为采用成熟商业软件IWD处理后的水深值；dM为采用LM或EM参数优化方法处理后所得水深值。
通过选取以上两段剖面进行处理可以明显看出，图 10中试验区域A的水下地形较为平坦，其深度沿着远离岛屿的方向逐渐增加，剖面中并未出现明显的地形起伏。试验区域B的底质情况相比试验区域A的情况明显不同，其水底底质以珊瑚礁为主，地形变化复杂且不具有明显的规律[39]。
图 11-图 14为以上各个剖面在不同波形处理方法条件下所得结果与IWD的偏差。值得注意的是，其中试验区域A中的数据处理结果的稳定性较好，而采用波形分解方法所得数据结果与系统随机软件IWD处理的结果虽有一定差异，但具有良好的符合程度(图 11和图 12)。而试验区域B中的两条剖面的水深偏差值稳定性相对较差，且剖面中出现水深估计异常点的分布较为集中。图 13和图 14中水深异常出现明显增多的区段分别对应于图 9中的ab段和cd段，此部分为两剖面对应水深地形最为复杂的部分，底质主要为珊瑚和珊瑚礁。造成这一结果的主要原因在于复杂的底质环境致使对应区域波形数据的起伏形态以及展宽受地形起伏变化和多路径干扰严重，从而降低波形分解过程中参数初值的稳定性，继而对水深归算精度产生影响。
图 11 试验区A内剖面1的处理结果与IWD处理结果对比Fig. 11 Comparison of the processing result of the profile 1 in experimentAand IWD
图选项
图 12 试验区A内剖面2处理结果与IWD处理结果对比Fig. 12 Comparison of the processing result of the profile 2 in experimentAand IWD
图选项
图 13 试验区B内剖面1处理结果与IWD处理结果对比Fig. 13 Comparison of the processing result of the profile 1 in experimentBand IWD
图选项
图 14 试验区B内剖面2处理结果IWD处理结果对比Fig. 14 Comparison of the processing result of the profile 2 in experimentBand IWD
图选项
以上试验表明海底状况对测深结果的精度以及水深估计稳定性具有较为明显的影响，但在实际扫描测深过程之前无法提前选择或避免复杂的海底地形，故可通过后期点云数据的相关处理在一定程度上减少数据异常对结果造成的不利影响。
表 3中，E(Δ)和σΔ主要表现了深度估计误差的总体偏差以及偏差的幅度。通过表中数据可以发现，尽管LM方法的测深误差略小于EM方法，但两种方法的误差偏差总体上均在0.1 m范围之内，说明两种波形分解参数优化方法的结果精度相当。另外，LM处理后的数据偏差的标准差略高于EM，说明EM所得结果的数据稳定性较好，这与模拟试验结果是一致的。LM方法在计算结果的精度稳定性上略逊于EM方法，可能的原因在于所处理波形的质量对参数优化方面的干扰，进而影响到波形分解乃至水深估计的效果。地形的剧烈变化、混响以及水质状况不佳等因素是降低波形质量的重要因素，因此，针对测深环境的具体情况选择适合的参数优化方法，可以在一定程度上提高水深探测的数据结果质量。
表 3 水深估计偏差统计Tab. 3 Statistical results of water depth estimation deviation
| m | |||||||||||||||||||||||||||||
| 试验区 | 剖面 | 优化方法 | E(Δ) | σΔ2 | σΔ | ||||||||||||||||||||||||
| 试验区A | 剖面1 | EM | 0.044 3 | 0.063 152 | 0.245 7 | ||||||||||||||||||||||||
| LM | 0.031 3 | 0.060 368 | 0.251 3 | ||||||||||||||||||||||||||
| 剖面2 | EM | 0.059 8 | 0.006 691 | 0.081 8 | |||||||||||||||||||||||||
| LM | 0.048 8 | 0.010 63 | 0.103 1 | ||||||||||||||||||||||||||
| 试验区B | 剖面1 | EM | 0.091 3 | 0.071 663 | 0.267 7 | ||||||||||||||||||||||||
| LM | 0.061 1 | 0.082 025 | 0.286 4 | ||||||||||||||||||||||||||
| 剖面2 | EM | 0.124 7 | 0.061 405 | 0.247 8 | |||||||||||||||||||||||||
| LM | 0.096 4 | 0.066 358 | 0.257 6 | ||||||||||||||||||||||||||
表选项
从总体上看，采用两种参数优化方法的水深估计结果与较为成熟的水深数据处理系统IWD的结果偏差远小于国际海道测量组织(International Hydrographic Organization, IHO)对水深测量的精度标准[40]。因此，两种参数优化方法在水深估计方面的响应不会对区域水深估计效果造成明显影响，采用波形分解方法提取水深值，可符合相关应用实际对测量结果的要求。
5 结论
机载激光测深全波形数据分解中的参数优化方法是准确估计分量参数，提高数据处理精度的重要环节。为进一步了解波形分解过程中不同参数优化方法在处理步骤和拟合效果方面的主要区别，本文采用适合浅水区域测量的Aquarius测深系统数据进行试验分析，代表性地选取了LM和EM两种参数优化方法进行对比研究。通过相关试验分析对LM和EM方法在技术原理及应用效果方面的主要特征进行了详细探讨。试验中，波形分解所得反射分量经两种参数优化方法处理后在振幅强度、波宽以及反射位置等重要特征方面均有所不同，其主要原因在于两种方法所采用的优化标准不同。就波形拟合效果而言，模拟数据分析和实测数据试验均说明LM方法通常可得到较EM方法更接近原始波形的拟合结果，有利于进一步的波形分析，但受波形数据质量的影响较大，容易造成结果的不稳定。另外，从总体效果看，采用LM或EM参数优化方法处理后所得的水深差异值并没有在具体的扫描探测条件下表现出显著的规律性。经实测数据对比，两种方法对水深探测结果的影响程度相当，而采用EM方法所得水深估计结果较LM方法更为稳健。因此，当扫描区域内水下地形状况复杂或系统电流干扰较为明显时，接收器得到的回波波形容易受到随机噪声的影响，此时应尽量选择数据处理结果更加稳定的EM方法，而当激光探测环境较好且系统较为稳定时，为保证波形分解后取得较好的拟合结果与水深探测精度，可采用LM方法作为高斯波形分解后参数初值的优化方法。相关结论也同样适用于不同机载激光测深设备对浅水区域测量获取的回波波形数据处理。
【引文格式】郭锴, 刘焱雄, 徐文学, 等. 机载激光测深波形分解中LM与EM参数优化方法比较. 测绘学报，2020，49(1)：117-131. DOI: 10.11947/j.AGCS.2020.20180242
