CN110031856B - 一种机载LiDAR测深数据的漫衰减系数提取方法 - Google Patents

Abstract

本发明公开了一种机载LiDAR测深数据的漫衰减系数提取方法，属于海洋测绘技术领域，其首先获取机载LiDAR测深原始波形数据，对其进行去噪；然后利用分层异构模型对去噪后的波形数据进行拟合，分解得到水体散射回波信号；再通过水体散射回波信号分别计算得到双段指数所对应的初始漫衰减系数；最后根据双段指数的采样时间确定相应权重，对初始漫衰减系数取加权平均，得到该激光测深点最终的漫衰减系数值。本发明通过这种方法，实现了机载LiDAR测深数据每一束激光的漫衰减系数提取，有效解决了目前漫衰减系数提取精度低、分辨率低的问题。

Description

一种机载LiDAR测深数据的漫衰减系数提取方法

技术领域

本发明属于海洋测绘技术领域，具体涉及一种机载LiDAR测深数据的漫衰减系数提取方法。

背景技术

漫衰减系数是水体的固有光学属性，反映了激光在水体中受散射、吸收所引起的能量衰减变化，是分析水质的重要参数之一。漫衰减系数主要取决于水体中的浮游生物含量(叶绿素浓度)、悬浮泥沙含量(混浊度)、营养盐含量(黄色物质、溶解有机物质、盐度指标)以及其他污染物等因素，研究漫衰减系数可对水体富营养化、水体混浊度骤增以及水体污染等环境变化实现实时、有效监测，在海洋科学研究、海洋工程应用、海洋渔业研究、海岸带环境监测等方面发挥重要的作用。

由于水体的流动性，导致了水中物质的非均匀分布，增加了水体中漫衰减系数的获取难度。目前国内外针对漫衰减系数的提取主要通过实地测量、被动遥感以及主动遥感3种方式进行。实地测量利用赛奇盘、透射计等工具测得漫衰减系数，这些方式需以测量船为载体，效率较低且覆盖范围小，难以获取大面积水域的漫衰减系数空间分布；通过卫星遥感影像数据可间接反演漫衰减系数，这种被动遥感方式可对大范围海域的漫衰减系数进行求取，但是无法在阴雨天气、夜间实施测量，且由于卫星设备本身的限制，测量精度难以保证，成果分辨率较低；运用全波形机载LiDAR(Light Detection And Ranging，激光雷达)测深这种主动遥感技术可实现漫衰减系数的快速反演，现有的反演方法受噪声影响严重，且需要水底底部回波强度和深度信息。

因此，有必要新提出一种机载LiDAR测深数据的漫衰减系数提取方法，以提高水体漫衰减系数的提取精度。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种机载LiDAR测深数据的漫衰减系数提取方法，设计合理，解决了在漫衰减系数提取过程中精度低、分辨率低的问题，实现了大面积水域漫衰减系数空间分布的精确获取。

为了实现上述目的，本发明采用如下技术方案：

一种机载LiDAR测深数据的漫衰减系数提取方法，包括以下步骤：

步骤1：获取机载LiDAR测深原始数据，读取测深数据中的全波形数据，对机载LiDAR测深原始波形数据进行去噪，得到去噪后的波形数据；

步骤2：通过分层异构模型的机载LiDAR波形拟合算法，利用步骤1中去噪后的波形数据，分解得到水体散射回波信号；具体包括如下步骤：

步骤2.1：采用极大值检测法对去噪后的波形数据进行峰值探测，确定包括水面和水底反射回波的初始峰值及其位置在内的初始参数；极大值检测法检测的极大值点位置如公式(1)所示：

k＝find(diff(sign(diff(w)))<0+1) (1)；

其中，w为去噪后的波形数据；diff表示差分运算；sign表示符号函数；find指查找满足条件的数据编号；

步骤2.2：针对步骤2.1中峰值探测所确定的初始参数，根据水面反射回波和水底反射回波的相应特性，采用高斯函数对水面反射回波和水底反射回波进行拟合，利用Levenberg-Marquardt算法进行参数优化，分解得到水面和水底反射回波信号；高斯函数如公式(2)所示：

其中，i＝s表示水面反射回波的拟合结果；i＝b表示水底反射回波的拟合结果；A、μ和σ分别为水面或水底高斯拟合函数的波峰值、波峰位置和半幅波宽；

步骤2.3：使用去噪后的波形数据剔除步骤2.2中的水面和水底反射回波信号，获得初始水体散射回波信号；

步骤2.4：根据步骤2.3中确定的初始水体散射回波信号，综合光在水体中的衰减特性和拟合效率问题，采用双指数函数对初始水体散射回波信号进行拟合，利用Levenberg-Marquardt算法进行参数优化，得到水体散射回波信号；双指数函数如公式(3)所示：

其中，a_x、b_x、c_x、d_x分别表示双指数函数4个顶点A、B、C、D的位置；b_y、c_y、d_y分别为B、C、D的振幅值；

步骤3：根据步骤2中获得的水体散射回波信号，分别计算得到双段指数对应的初始漫衰减系数；具体包括如下步骤：

步骤3.1：利用水体散射回波信号反演水体的漫衰减系数，确定机载LiDAR测深系统在深度为D_i处的水体散射回波强度方程；水体散射回波强度方程如公式(4)所示：

其中，D_i表示水层深度；P_c(D_i)表示深度D_i时的水体散射回波强度；P_e表示激光发射强度；

表示大气双程损失因子；A_R表示激光接收器的视场面积；η_e和η_R分别表示激光发射器的光学效率和激光接收器的光学效率；F表示接收视场角损失因子；L_S是通过表面传输的损失，即表面反照率；β(Φ)表示体散射函数；θ表示激光传入水体时的入射角，即天底角；θ_w表示射入水体后的折射角；n_w表示水体折射率；H表示飞机航高；K_d表示532nm波段处的漫衰减系数；c表示激光在真空的传播速度；t_i表示水层深度D_i处的时刻；t_s表示水面反射回波的峰值位置；

步骤3.2：根据步骤3.1中所述水体散射回波强度方程，将公式(5)代入公式(4)，并简化公式(4)，得到简化后的水体散射回波强度方程如公式(6)所示：

其中，

步骤3.3：将步骤2.4中所得双指数函数与步骤3.2中简化后的水体散射回波强度方程描述为指数函数的形式，如公式(7)所示：

y(t)＝m·e^nt (7)；

其中，y(t)为指数函数，m、n分别为指数函数的系数；

步骤3.4：根据步骤2.4中所得双指数函数，当b_x≤t<c_x时，f_c(t)＝P_c(t)，由公式(7)可知，公式(3)与公式(6)中的指数函数系数是对应相等的，推导得到公式(8)：

步骤3.5：对步骤3.4中公式(8)做方程变换，得到该段指数所对应的初始漫衰减系数K_d1，同理，计算得到另一段指数所对应的初始漫衰减系数K_d2，如公式(9)和公式(10)所示：

步骤4：根据步骤2中获得的水体散射回波信号，通过双段指数的采样时间确定相应权重，根据公式(13)，对步骤3中双段指数所对应的初始漫衰减系数K_d1与K_d2取加权平均，得到K_d1与K_d2的加权平均值，将此加权平均值作为该激光测深点最终的漫衰减系数值：

双段指数的采样时间如公式(11)和公式(12)所示：

Δt₁＝c_x-b_x (11)；

Δt₂＝d_x-c_x (12)；

初始漫衰减系数K_d1与K_d2的加权平均值如公式(13)所示：

本发明所带来的有益技术效果：

本发明方法，与现有技术相比，考虑到激光在水体中的衰减特性以及计算效率，通过分层异构模型的机载LiDAR波形拟合算法分解得到水体散射回波信号，并用其与水体散射回波强度方程构建相应模型，计算得到每个激光测深点的漫衰减系数，解决了在漫衰减系数提取过程中精度低、分辨率低的问题，实现了水体漫衰减系数的精确反演。

附图说明

图1为本发明一种机载LiDAR测深数据的漫衰减系数提取方法的流程图。

图2为机载LiDAR测深波形示意图。

图3为本发明中机载LiDAR测深波形去噪前后对比图。

图4为本发明中双指数函数拟合水体散射回波示意图。

图5为本发明中机载LiDAR测深波形拟合结果示意图。

图6为本发明中所提取的漫衰减系数空间分布图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

本发明提供了一种机载LiDAR测深数据的漫衰减系数提取方法，其流程如图1所示。机载LiDAR测深系统的回波波形包括水面反射回波、水体散射回波、水底反射回波以及噪声四部分，如图2所示。漫衰减系数提取方法具体包括以下步骤：

步骤1：获取机载LiDAR测深原始数据，读取测深数据中的全波形数据，对机载LiDAR测深原始波形数据进行去噪，得到去噪后的波形数据。

具体的，采用机载LiDAR测深原始波形数据。机载LiDAR测深技术具有高精度、高分辨率、灵活机动、快速高效的优势，其激光接收器能够以数字化的形式记录每个激光脉冲的全部回波信号。

具体实施时，由于机载LiDAR测深系统一般采用扫描方式测量，通过扫描镜的局部运动，实现测深点的条带式展宽，因而获得海量测深回波信号。水体中复杂的环境因素导致激光脉冲在水中传播时会产生大量噪声，其存在会导致拟合函数初始参数估计不准等问题，容易获得含有粗差的拟合波形，影响拟合精度，因此，需要对机载LiDAR测深回波信号进行相应的波形去噪，最终得到去噪后的测深波形数据，如图3所示。

步骤2：通过分层异构模型的机载LiDAR波形拟合算法，利用步骤1中去噪后的波形数据，分解得到水体散射回波信号。

具体的，波形拟合是机载LiDAR测深数据处理的重要环节，是漫衰减系数提取的基础。根据机载LiDAR测深回波信号各部分的波形特征，选择合适的函数构建分层异构模型(水面—高斯函数、水体—双指数函数及水底—高斯函数)进行拟合，从而将叠加在水面反射回波和水底反射回波中的水体散射回波精确地提取出来，得到水体散射回波所对应的函数。

进一步的实施例中，步骤2具体包括如下步骤：

步骤2.1：采用极大值检测法对去噪后的波形数据进行峰值探测，确定包括水面和水底反射回波的初始峰值及其位置在内的初始参数；极大值检测法检测的极大值点位置如公式(1)所示：

k＝find(diff(sign(diff(w)))<0+1) (1)；

其中，w为去噪后的波形数据；diff表示差分运算；sign表示符号函数；find指查找满足条件的数据编号；

步骤2.2：针对步骤2.1中峰值探测所确定的初始参数，根据水面反射回波和水底反射回波的相应特性，采用高斯函数对水面反射回波和水底反射回波进行拟合，利用Levenberg-Marquardt算法进行参数优化，分解得到水面和水底反射回波信号；高斯函数如公式(2)所示：

其中，i＝s表示水面反射回波的拟合结果；i＝b表示水底反射回波的拟合结果；A、μ和σ分别为水面或水底高斯拟合函数的波峰值、波峰位置和半幅波宽；

步骤2.3：使用去噪后的波形数据剔除步骤2.2中的水面和水底反射回波信号，获得初始水体散射回波信号；

步骤2.4：根据步骤2.3中确定的初始水体散射回波信号，综合光在水体中的衰减特性和拟合效率问题，采用双指数函数对初始水体散射回波信号进行拟合，利用Levenberg-Marquardt算法进行参数优化，得到水体散射回波信号；双指数函数如公式(3)所示：

其中，a_x、b_x、c_x、d_x分别表示双指数函数4个顶点A、B、C、D的位置；b_y、c_y、d_y分别为B、C、D的振幅值。

具体实施时，考虑到机载LiDAR测深系统发射和接收的激光脉冲信号近似服从高斯分布，因此，根据水面和水底反射回波相应特性分析，本发明采用高斯函数拟合得到水面、水底反射回波信号。从去噪后的机载LiDAR测深波形数据中剔除水面、水底反射回波信号，即可得到初始水体散射回波。由于水体浑浊度在垂直剖面上非均匀分布导致激光脉冲的衰减程度不同，综合激光在水体中衰减特性和拟合效率的问题，如图4所示，本发明将水体沿垂直剖面分为两层，利用双指数函数ABCD拟合得到水体散射回波信号。遵循非线性最小二乘方法，采用Levenberg-Marquardt算法对波形拟合的13个初始参数(水面-高斯函数3个、水底-高斯函数3个、水体-双指数函数7个)进行迭代优化，直到拟合偏差达到最小为止，获得各参数的精确值。将优化后的水体散射回波参数代入公式(3)，得到最优的水体散射回波信号。如图5所示，为运用分层异构模型的机载LiDAR波形拟合算法所得到的波形拟合结果。

步骤3：根据步骤2中获得的水体散射回波信号，分别计算得到双段指数对应的初始漫衰减系数。

具体的，根据水体散射回波强度方程分析可知，激光脉冲回波强度随水深深度和漫衰减系数增加呈指数衰减。因此，可以通过水体散射回波信号所对应的双指数函数反演水体的漫衰减系数。

进一步的实施例中，步骤3具体包括如下步骤：

步骤3.1：利用水体散射回波信号反演水体的漫衰减系数，确定机载LiDAR测深系统在深度为D_i处的水体散射回波强度方程；水体散射回波强度方程如公式(4)所示：

其中，D_i表示水层深度；P_c(D_i)表示深度D_i时的水体散射回波强度；P_e表示激光发射强度；

表示大气双程损失因子；A_R表示激光接收器的视场面积；η_e和η_R分别表示激光发射器的光学效率和激光接收器的光学效率；F表示接收视场角损失因子；L_S是通过表面传输的损失，即表面反照率；β(Φ)表示体散射函数；θ表示激光传入水体时的入射角，即天底角；θ_w表示射入水体后的折射角；n_w表示水体折射率；H表示飞机航高；K_d表示532nm波段处的漫衰减系数；c表示激光在真空的传播速度；t_i表示水层深度D_i处的时刻；t_s表示水面反射回波的峰值位置；

步骤3.2：根据步骤3.1中所述水体散射回波强度方程，将公式(5)代入公式(4)，并简化公式(4)，得到简化后的水体散射回波强度方程如公式(6)所示：

其中，

步骤3.3：将步骤2.4中所得双指数函数与步骤3.2中简化后的水体散射回波强度方程描述为指数函数的形式，如公式(7)所示：

y(t)＝m·e^nt (7)；

其中，y(t)为指数函数，m、n分别为指数函数的系数；

步骤3.4：根据步骤2.4中所得双指数函数，当b_x≤t<c_x时，f_c(t)＝P_c(t)，由公式(7)可知，公式(3)与公式(6)中的指数函数系数是对应相等的，推导得到公式(8)：

步骤3.5：对步骤3.4中公式(8)做方程变换，得到该段指数所对应的初始漫衰减系数K_d1，同理，计算得到另一段指数所对应的初始漫衰减系数K_d2，如公式(9)和公式(10)所示：

具体实施时，对机载LiDAR水体散射回波强度方程进行简化，用公式(6)来表示。由于近岸水深通常小于30m，与飞机航高H(通常为400～500m)相比，参数D_i对W的影响可以忽略不计，因此，假设外业测量过程中飞机航高保持稳定，激光发射器的天底角以及接收视场角保持恒定，且所有损耗因素均得到较好的控制，W可认为是一个常数。但是W涉及未知参数太多，不能直接确定。通过分析简化后的水体散射回波强度方程可知，其与公式(3)相应的指数函数系数n是对应相等的。利用上述关系，通过方程变换即可分别求得BC、CD段的初始漫衰减系数。

步骤4：根据步骤2中获得的水体散射回波信号，通过双段指数的采样时间确定相应权重，根据公式(13)，对步骤3中双段指数所对应的初始漫衰减系数K_d1与K_d2取加权平均，得到K_d1与K_d2的加权平均值，将此加权平均值作为该激光测深点最终的漫衰减系数值：

双段指数的采样时间如公式(11)和公式(12)所示：

Δt₁＝c_x-b_x (11)；

Δt₂＝d_x-c_x (12)；

初始漫衰减系数K_d1与K_d2的加权平均值如公式(13)所示：

具体的，计算激光测深点漫衰减系数的加权平均值。将水体散射回波信号的双段指数所对应的初始漫衰减系数取加权平均，即可得到该激光测深点最终的漫衰减系数值。

具体实施时，根据水体散射回波信号所对应的双指数函数，以BC段、CD段的采样时间差(Δt₁与Δt₂)作为权重，利用公式(13)，即可得到该激光测深点最终的漫衰减系数值。如图6所示，为应用本发明所提取的漫衰减系数空间分布图。

综上所述，本发明提供了一种机载LiDAR测深数据的漫衰减系数提取方法，考虑到激光在水体中的衰减特性以及计算效率的问题，方法包括：获取机载LiDAR测深原始波形数据，对其进行去噪；利用分层异构模型的机载LiDAR波形拟合算法拟合去噪后的波形数据，分解得到水体散射回波信号；通过水体散射回波信号分别计算得到双段指数所对应的初始漫衰减系数；根据双段指数的采样时间确定相应权重，对初始漫衰减系数取加权平均，得到该激光测深点最终的漫衰减系数值；本发明实现了机载LiDAR测深数据每一束激光脉冲的漫衰减系数提取，有效解决了目前提取漫衰减系数时精度低、分辨率低的问题。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (1)

1.一种机载LiDAR测深数据的漫衰减系数提取方法，其特征在于：包括以下步骤：

步骤1：获取机载LiDAR测深原始数据，读取测深数据中的全波形数据，对机载LiDAR测深原始波形数据进行去噪，得到去噪后的波形数据；

步骤2：通过分层异构模型的机载LiDAR波形拟合算法，利用步骤1中去噪后的波形数据，分解得到水体散射回波信号；具体包括如下步骤：

步骤2.1：采用极大值检测法对去噪后的波形数据进行峰值探测，确定包括水面和水底反射回波的初始峰值及其位置在内的初始参数；极大值检测法检测的极大值点位置如公式(1)所示：

k＝find(diff(sign(diff(w)))<0+1) (1)；

其中，w为去噪后的波形数据；diff表示差分运算；sign表示符号函数；find指查找满足条件的数据编号；

步骤2.2：针对步骤2.1中峰值探测所确定的初始参数，根据水面反射回波和水底反射回波的相应特性，采用高斯函数对水面反射回波和水底反射回波进行拟合，利用Levenberg-Marquardt算法进行参数优化，分解得到水面和水底反射回波信号；高斯函数如公式(2)所示：

其中，i＝s表示水面反射回波的拟合结果；i＝b表示水底反射回波的拟合结果；A、μ和σ分别为水面或水底高斯拟合函数的波峰值、波峰位置和半幅波宽；

步骤2.3：使用去噪后的波形数据剔除步骤2.2中的水面和水底反射回波信号，获得初始水体散射回波信号；

步骤2.4：根据步骤2.3中确定的初始水体散射回波信号，综合光在水体中的衰减特性和拟合效率问题，采用双指数函数对初始水体散射回波信号进行拟合，利用Levenberg-Marquardt算法进行参数优化，得到水体散射回波信号；双指数函数如公式(3)所示：

其中，a_x、b_x、c_x、d_x分别表示双指数函数4个顶点A、B、C、D的位置；b_y、c_y、d_y分别为B、C、D的振幅值；

步骤3：根据步骤2中获得的水体散射回波信号，分别计算得到双段指数对应的初始漫衰减系数；具体包括如下步骤：

步骤3.1：利用水体散射回波信号反演水体的漫衰减系数，确定机载LiDAR测深系统在深度为D_i处的水体散射回波强度方程；水体散射回波强度方程如公式(4)所示：

其中，D_i表示水层深度；P_c(D_i)表示深度D_i时的水体散射回波强度；P_e表示激光发射强度；

表示大气双程损失因子；A_R表示激光接收器的视场面积；η_e和η_R分别表示激光发射器的光学效率和激光接收器的光学效率；F表示接收视场角损失因子；L_S是通过表面传输的损失，即表面反照率；β(Φ)表示体散射函数；θ表示激光传入水体时的入射角，即天底角；θ_w表示射入水体后的折射角；n_w表示水体折射率；H表示飞机航高；K_d表示532nm波段处的漫衰减系数；c表示激光在真空的传播速度；t_i表示水层深度D_i处的时刻；t_s表示水面反射回波的峰值位置；

步骤3.2：根据步骤3.1中所述水体散射回波强度方程，将公式(5)代入公式(4)，并简化公式(4)，得到简化后的水体散射回波强度方程如公式(6)所示：

其中，

步骤3.3：将步骤2.4中所得双指数函数与步骤3.2中简化后的水体散射回波强度方程描述为指数函数的形式，如公式(7)所示：

y(t)＝m·e^nt (7)；

其中，y(t)为指数函数，m、n分别为指数函数的系数；

步骤3.4：根据步骤2.4中所得双指数函数，当b_x≤t<c_x时，f_c(t)＝P_c(t)，由公式(7)可知，公式(3)与公式(6)中的指数函数系数是对应相等的，推导得到公式(8)：

步骤3.5：对步骤3.4中公式(8)做方程变换，得到该段指数所对应的初始漫衰减系数K_d1，同理，计算得到另一段指数所对应的初始漫衰减系数K_d2，如公式(9)和公式(10)所示：

步骤4：根据步骤2中获得的水体散射回波信号，通过双段指数的采样时间确定相应权重，根据公式(13)，对步骤3中双段指数所对应的初始漫衰减系数K_d1与K_d2取加权平均，得到K_d1与K_d2的加权平均值，将此加权平均值作为该激光测深点最终的漫衰减系数值：

双段指数的采样时间如公式(11)和公式(12)所示：

Δt₁＝c_x-b_x (11)；

Δt₂＝d_x-c_x (12)；

初始漫衰减系数K_d1与K_d2的加权平均值如公式(13)所示：

Images