论文推荐|申二华：圆扫描式机载激光测深系统检校模型及仿真分析

    为了提升圆扫描式机载激光测深系统的定位精度，提出了一种检校思路：在平面区域获取激光点云时，系统误差和随机误差使得本应共面的激光点云不再共面，通过将激光点云拟合到单个平面上达到纠正点云位置的目的。首先推导了简单模式下圆扫描式机载激光测深系统定位模型，并从直线与平面交会的数学原理出发模拟激光光线与海面的交会过程，根据折射原理解算激光光线在水中的方向矢量，联合激光光线在水中的直线方程和海底面数学方程解算激光脚点的位置。然后，引入未知数先验方差，推导了参数加权最小二乘平差模型，为后面检校模型的解算奠定基础。最后，推导了用于检校的平差数学模型和详细的计算过程，对检校过程进行了模拟计算和分析讨论，得出了一些有益于检校过程的结论。

机载激光扫描测距将姿态测量装置、差分GPS、激光测距仪等技术产品集成一体，可以直接获取地面三维信息^[]。将地形扫描激光波长由红外波段(1064nm)倍频至蓝绿波段(532nm)，同时优化激光回波信号接收装置，可将机载激光扫描测距技术应用于浅水(0~50m，与水质有关)测量^[]。

激光测深系统的激光扫描装置通过扫描镜的局部运动，实现激光测深点的条带式测量。一般而言，激光测深扫描模式分为直线扫描、椭圆扫描和圆扫描，圆扫描相比直线扫描和椭圆扫描机械结构最为简单^[]。为了提高测量效率, 需要扫描装置高速运转, 而机械结构越复杂, 扫描运行就越不易平稳, 但为确保测量点精确定位, 必须保证扫描机构的平稳运行^[]。随着计算机技术的发展，数据处理的速度已经得到很大提升，因此采用机械结构最简单的圆扫描模式可以提高测量效率。

对激光测深技术而言，精确测量激光脚点在海底的三维坐标是其最基本也是最主要的目的。然而由于GPS、INS、LiDAR的高度集成，定位精度受制于硬件部分的精度(主要由生产厂商决定)和视准轴偏角等各种系统误差和偶然误差的影响，需要使用检校方法确定这些系统误差^[]。目前的LiDAR检校方法较为费时烦琐，主要是条带平差法。该方法利用条带间的重叠部分进行最小二乘平差^[-]，如基于连接点模型的LiDAR检校方法^[-]，随着特征提取方法的发展，基于路或建筑物表面等多个线面特征在相邻观测航线上提取特征并校正的方法也取得了一定的进展^[-]。本文尝试使用参数加权最小二乘方法使点云拟合在一个平面上检校系统误差，这与文献[—]使用多个平面处理点云间位置关系的方法有所不同。

1 简单模式下的机载激光雷达测深系统定位方程

为简单模式下的机载激光测深雷达扫描方式及所在的激光扫描坐标系，在大气和海水中的传播距离分别为ρ₁和ρ₂，α为激光固定发射角，θ_sz为扫描角。中激光光线方向矢量l₁穿过水面折射后变为l₂，θ₁、θ₂分别为激光方向矢量与海水表面法向量p的夹角。n_a、n_w分别为大气和海水的折射率。海底激光脚点在激光扫描坐标系中的坐标可表示为

图 1 简单模式下的机载激光测深雷达扫描方式Fig. 1 The airborne bathymetry LiDAR scanning pattern in the simple mode

图 2 激光光线矢量穿过海平面的方向变化Fig. 2 The laser vector change through the sea surface

 (1)

式中

 (2)

 (3)

根据文献[]提出的折射定律矢量表达方程，得到静止海水表面的折射公式

 (4)

 (5)

 (6)

则通过静止海面后激光光线的方向矢量为

 (7)

激光测深系统将姿态测量系统、差分GPS接收机、激光扫描仪集成于一体。理想情况下，姿态测量系统坐标轴与激光扫描仪坐标轴平行，实际情况下却存在安置角常差，通常称为视准轴偏角。同样，差分GPS接收机天线中心与激光扫描中心也存在偏移量，其偏差称为偏心分量。设视准轴偏角对应的3个旋转角参数分别为航向偏角θ_ix、滚动偏角θ_iz、俯仰偏角θ_iy，激光扫描中心至GPS/INS中心的偏心分量为ΔX_il=x_ily_ilz_il^T，则INS坐标系中激光脚点坐标为(详细的各个坐标系间关系的描述参考文献[—]，坐标系间的数学转换关系可参考[18—20])

 (8)

式中，R(·)为旋转角矩阵。

最后，根据GPS/INS测定的飞机在成图坐标系中的姿态变化，即航向角ψ_vx、俯仰角ψ_vy、滚动角ψ_vz，和GPS/INS中心在成图坐标系的位置ΔX=x_Iy_Iz_I^T，得到激光脚点在成图坐标系中的坐标为

 (9)

因此，简单模式下圆扫描式机载激光测深系统的定位模型可表述为

 (10)

2 点云模拟模型

点云可以看成是激光光线与平面交会产生的，海底激光点云经过了海表和海底两个平面，下面分别推导激光光线与海表面、海底面交会的数学表达式。

2.1 海表面激光点云模型

激光光线和海表面的交会如所示，激光光线于时间t=t₀与海表面交会，[x_Ly_Lz_L]为直线上任意一点，[a_Lb_Lc_L]^T为激光光线的方向，平面法向量为[a_Pb_Pc_P]^T，[x_Py_Pz_P]为平面上任意点，则激光光线和海表面的数学表达式分别为

图 3 激光光线和平面交会图Fig. 3 The intersection of laser line and the plane

 (11)

 (12)

解算得到

 (13)

对于机载激光测深系统，假设某个时刻海表面法向量为G_n，飞机激光发射点的海表面垂点坐标为N_pt，激光发射点坐标为L_pt，激光方向矢量为L_u，则激光光线与海表面交会的时间有如下表达

 (14)

在球面坐标系下，海表面法向量可表示为

 (15)

设飞机开始飞行时其正下方位于海表面的点为原点，则

 (16)

由激光脚点定位模型得到激光发射点坐标

 (17)

激光光线方向矢量为

 (18)

根据以上表达式即可解算激光光线与海表面交会的位置。

2.2 水底激光点云模型

激光穿透海面发生折射，通过折射公式(7)得到水中传输时的激光方向。由于激光与水面的交点已经确定，因此可以得到海水中激光光线直线方程，此时只需确定模型海底面的数学表达式即可得到海底激光脚点位置。解算步骤与海表面点云模型类似。

3 参数加权最小二乘平差模型

对于加权最小二乘，观测值残差平方和可表达为

 (19)

式中，r=[r₁r₂r₃…r_n]为观测值残差；C_l=diag(σ²_l1, σ²_l2, σ²_l3, …, σ²_ln)；σ²_li为第i个观测值对应的方差。

加入未知数加权改正数，得到

 (20)

式中，δ为未知数改正值向量；C_x为未知数先验方差组成的对角矩阵。加入未知数权值，可以使最小二乘算法先调整权重高的参数，也使得最小二乘算法避免向无效数值收敛。

待估参数X与观测值L之间组成的函数关系式为

 (21)

由泰勒级数展开得到一阶线性逼近

 (22)

 (23)

式中，g为模型初始值；r为观测值改正值；δ为未知数改正值；D为观测值对应的偏导矩阵；A为未知数对应的偏导矩阵。

由矩阵分块求解最终得到未知数δ和观测值残差r

 (24)

 (25)

4 检校模型及计算过程

在平面区域获取激光点云时，传感器的系统误差和随机误差使得本该共面的激光点云不共面，使用最小二乘模型将不共面的激光点云拟合到单个平面上，从而降低系统误差和随机误差的影响，进而纠正点云位置。

如所示，n是平面法向量，x_o是激光点云，x_p为平面上的任意点，若所有激光点均位于平面上，则任意点x_o与x_p相减构成的向量与平面法向量n的向量积为0，得到

图 4 激光点云和拟合平面示意图Fig. 4 The laser point clouds and fitting plane

 (26)

实际上，因为系统误差和随机误差造成在平面区域获取的激光点云不共面，使得式(26)不为0，需要使用检校算法纠正点云位置。

为了方便分析，这里将偏心分量设为0，设海平面是静止水平的，将海平面作为检校平面，本文定位模型进一步简化，即

 (27)

对于平面法向量，使用球坐标系，n=[ρcosθsinϕρsinθsinϕ-ρcosϕ]^T，设ρ=1，则n=[cosθsinϕsinθsinϕ-cosϕ]^T，未知数由3个降为2个，则

 (28)

式中，x_p=[x_py_pz_p]^T；观测值向量是L=[ρ₁θ_szψ_vxψ_vyψ_vzx_Iy_Iz_I]^T；未知数向量是X=[θ_ixθ_iyθ_izφ_pθ_px_py_pz_p]^T；其他参数参照前文。

根据参数加权最小二乘模型，下面需要分别求出模型初始值g，观测值对应的偏导矩阵D，未知数对应的偏导矩阵A。

模型初始值g为

 (29)

观测值偏导矩阵D为

 (30)

式中，D中非对角线部分每个元素0均为1×8零矩阵。

矩阵A为f对未知数X=[θ_ixθ_iyθ_izpθ_px_py_pz_p]^T的偏导矩阵。

观测值权矩阵为

 (31)

式中，C_l中非对角线部分每个元素0均为1×8零矩阵。

 (32)

未知数权矩阵为

 (33)

根据式(29)—式(33)和式(24)，即可解算未知数改正数。

假设存在n组观测值，则式(24)中(DC_lD^T)^-1为n×n矩阵的求逆，在本文后面模拟试验中，n最少等于1000，这对于计算机是巨大的工作量。由于D、C_l矩阵中均含有大量0元素，可以尝试将计算过程简化，由于C_l为对角线矩阵，得到

 (34)

式中，DC_l中非对角线部分每个元素0均为1×8零矩阵，T_i为对角线矩阵，不再详列。

 (35)

式中，S_i的计算式不再详列。

由于DC_lD^T为对角线矩阵，则(DC_lD^T)^-1=diag(1/S₁, 1/S₂, …, 1/S_n)

使用上式对A^T(DC_lD^T)^-1的计算过程优化，效率有质的提升。

5 点云模拟及检校模型数值计算

下面使用文章中的定位模型、点云模拟模型仿真点云的生成过程(—)。以Optech公司的CZMIL机载激光测深系统^[-]为例：航高H为400m，水文模式下激光脉冲发射频率R_prf为10kHz，棱镜上表面斜坡坡度为39.16°，脉冲固定发射角度θ为20°，扫描仪转速R_rps为27Hz(圈/s)，飞行速度为140Kts(海里/小时)。空气折射率为1.0003，海水折射率为1.334，棱镜折射率为1.462，模拟点云500圈。本文主要对视准轴偏角进行探讨，假设不存在其他系统误差，分别使3个视准轴偏角设为1°，模拟在飞机平直飞行、航高从400m逐渐升至500m、飞机3个姿态角分别由0°至30°逐渐变化的点云分布。对于每个视准轴偏角，应有5对点云分布图，限于篇幅原因，只列出部分典型情况的点云分布图(所谓“典型”是指能够将本应共面的点云较大幅度的非共面分散，未列出的其他情况点云仍共面或非共面分散幅度较小)，其中左图为立体图，右图为航向视角的YZ平面图，蓝色为存在视准轴偏角时获取的点云，红色为俯仰偏角未知时计算得到的点云。

图 5 视准轴俯仰偏角为1°时平直飞行、航高从400m逐渐升至500m的点云分布图Fig. 5 The point clouds map with the flying height level when the boresight misalignment pitching angle is 1°

图 6 视准轴滚动偏角为1°时飞机姿态的航向角由0°至30°的点云分布图Fig. 6 The point clouds map with the vessel heading 0° up to 30° when the boresight misalignment rolling angle is 1°

图 7 视准轴航向偏角为1°时飞机姿态的滚动角由0°至30°的点云分布图Fig. 7 The point clouds map with the vessel rolling 0° up to 30° when the boresight misalignment heading angle is 1°

基于以上点云的模拟得到以下检校思路：在存在视准轴俯仰偏角时，保持平直飞行，即可使本应共面的点云不共面，从而探测该偏角；在存在视准轴滚动偏角时，检校该偏角需要使航向姿态角发生变化；在存在视准轴航向偏角时，检校该偏角需要滚动姿态角发生变化。总的来说，在检校飞行时同时使航向角和滚动角发生动态变化即可基于本文检校思路解算3个视准轴偏角。

基于上述思路，采用文中激光点云定位方程和点云模拟模型，使飞行航向角和翻滚角同时由0°向10°均匀发生变化，为了防止最小二乘解算方程数据量过大，激光发射频率设为54Hz(扫描仪转速R_rps为27Hz，避免非整数近似操作)，使3个视准轴偏角都设为5°，得到所示的点云轨迹，共1000个点，同时记录1000对观测值。其中左图为立体图，右图为航向视角的YZ平面图，蓝色为存在视准轴偏角时获取的点云，红色为俯仰偏角未知时计算得到的点云。

图 8 3个视准轴偏角均为5°时飞机姿态的滚动角和航向角由0°至10°均匀变化的点云分布图Fig. 8 The point clouds map with the vessel rolling and vessel heading 0 degree up to 10 degree when the three boresight misalignment angle are all 5 degree

参照国际上广泛应用的相关仪器系统指标，对各观测值和待求未知数标准差赋值。激光测距仪测距标准差设为0.01，扫描角标准差设为0.002，GPS/INS定位定姿系统的位置标准差和姿态标准差都采用广泛使用且精度较高的POS/AV^TM 610后处理之后的系统参数()。视准轴偏角θ_ix、θ_iy、θ_iz标准差使用文献[—]论文里的典型数值()。

表 1 观测值标准差Tab. 1 The observed value STDev

σρ1	σθsz	σψvx	σψvy	σψvz	σxI	σyI	σzI
0.01	0.002	0.0028	0.0083	0.0083	0.05	0.05	0.075

表 2 未知数标准差Tab. 2 The unknown parameter STDev

σθix	σθiy	σθiz	σϕp	σθp	σxp	σyp	σzp
0.05	0.05	0.05	1	1	0.1	0.1	0.1

根据上述获取的1000对观测值及对应的标准差(每对观测值标准差相同)和未知数标准差，采用文中参数加权最小二乘平差模型解算，结果如所示。将频率分别增大到2倍、4倍、10倍、20倍、200倍，视准轴偏差标准差统计结果如所示。将飞机姿态角航向角、滚动角变化幅度增大至2倍，视准轴偏差标准差统计结果如所示。分别同方向和反方向飞行两条航带，航带间隔为1000m，视准轴偏差标准差统计结果如所示。

表 3 未知数计算结果Tab. 3 the computing results of the unknown parameters

θix	θiy	θiz	ϕp	θp	xp	yp	zp
4.995	4.998	4.996	-4.67e-10	-0.0252	-2.79e-08	-2.44e-10	2.14e-07

表 4 不同激光发射频率下的视准轴偏角标准差Tab. 4 The three boresight misalignment STDev computed with different laser emission frequency

频率	点数	σθix	σθiy	σθiz
54Hz	1000	0.042	0.008	0.027
108Hz	2000	0.029	0.006	0.019
216Hz	4000	0.022	0.004	0.014
540Hz	10000	0.013	0.003	0.008
1080Hz	20000	0.007	0.001	0.005
10800Hz	200000	0.003	0.0005	0.002

表 5 滚动姿态角和航向姿态角变化加倍时的视准轴偏角标准差Tab. 5 the three boresight misalignment STDev computed when the vessel rolling and heading angle changes are doubled

姿态角变化	点数	σθix	σθiy	σθiz
ψvx(0°-10°)ψvz(0°-10°)	1000	0.042	0.008	0.027
ψvx(0°-20°)ψvz(0°-10°)	1000	0.041	0.007	0.014
ψvx(0°-10°)ψvz(0°-20°)	1000	0.022	0.008	0.026
ψvx(0°-20°)ψvz(0°-20°)	1000	0.021	0.008	0.013

表 6 双航带时的视准轴偏角标准差Tab. 6 The three boresight misalignment STDev with the two flight strip

航带	点数	σθix	σθiy	σθiz
单航带	1000	0.042	0.008	0.027
反方向两航带	2000	0.024	0.005	0.004
同方向两航带	2000	0.026	0.006	0.011

6 分析与讨论

中，随着激光发射频率的增加，点密度的提升能够降低随机误差的影响，逐渐提高视准轴偏角检校精度。由于CZMIL水文模式设计激光频率为10kHz，此时检校精度颇高，但计算时间也大幅度增加至约27min，因此在精度要求不是很高的情况下，可以重采样激光点云从而降低点云数以提升检校效率。

中，在飞机航向姿态角变化至20°时，视准轴滚动偏角标准差有效降低，航向偏角和俯仰偏角的标准差几乎不变；同样在飞机滚动姿态角变化至20°时，视准轴航向偏角标准差有效降低，俯仰偏角和滚动偏角的标准差几乎不变，刚好验证了点云模拟时关于检校方式的结论。

中，在飞机反方向飞行两次时，航向姿态角类似于突然改变了180°，视准轴滚动偏角检校精度有了质的提升，视准轴航向偏角和俯仰偏角类似于增加了点密度，检校精度也得到了提升。在飞机同方向飞行两次时，较单航带时点数量增多且分散范围更大，视准轴滚动偏角检校精度也得到了大幅提升，说明滚动偏角对于航向角的变化和点云的位置分散情况均相当敏感，视准轴航向偏角和俯仰偏角变化程度类似于反方向飞行两次的情况。

因此为了提高视准轴偏角检校的精度，需要在计算量允许的情况下保持较高的点密度，增加飞机航向姿态角和滚动姿态角变化幅度的同时飞行多个航带。

7 结束语

本文在推导了简单模式下机载激光测深系统定位模型和模拟模型后，使用参数加权最小二乘方法使点云拟合在一个平面上检校机载激光测深系统误差，通过模拟试验初步验证了思路的正确性，并通过多项对比试验，对检校实施提出了有意义的建议。不足的地方有：只讨论了视准轴偏差存在的情况，对于其他可能的系统误差需要从定位模型开始继续深入研究，检校方法类似文中思路，只是未知数更多，对于各种系统误差之间的相关性也需要有针对性地深入研究。

全文内容请点击下方“阅读原文”，建议开启微信优化阅读。

申二华, 张永生, 李凯     

信息工程大学地理空间信息学院，河南 郑州 450052

收稿日期：2015-10-26； 修回日期：2016-06-20

第一作者简介：申二华(1988-)，男，博士生，研究方向为机载激光雷达测深。

E-mail:

更多精彩内容：