机载激光扫描测距将姿态测量装置、差分GPS、激光测距仪等技术产品集成一体,可以直接获取地面三维信息[1]。将地形扫描激光波长由红外波段(1064nm)倍频至蓝绿波段(532nm),同时优化激光回波信号接收装置,可将机载激光扫描测距技术应用于浅水(0~50m,与水质有关)测量[2]。
激光测深系统的激光扫描装置通过扫描镜的局部运动,实现激光测深点的条带式测量。一般而言,激光测深扫描模式分为直线扫描、椭圆扫描和圆扫描,圆扫描相比直线扫描和椭圆扫描机械结构最为简单[3]。为了提高测量效率, 需要扫描装置高速运转, 而机械结构越复杂, 扫描运行就越不易平稳, 但为确保测量点精确定位, 必须保证扫描机构的平稳运行[4]。随着计算机技术的发展,数据处理的速度已经得到很大提升,因此采用机械结构最简单的圆扫描模式可以提高测量效率。
对激光测深技术而言,精确测量激光脚点在海底的三维坐标是其最基本也是最主要的目的。然而由于GPS、INS、LiDAR的高度集成,定位精度受制于硬件部分的精度(主要由生产厂商决定)和视准轴偏角等各种系统误差和偶然误差的影响,需要使用检校方法确定这些系统误差[5]。目前的LiDAR检校方法较为费时烦琐,主要是条带平差法。该方法利用条带间的重叠部分进行最小二乘平差[6-7],如基于连接点模型的LiDAR检校方法[8-9],随着特征提取方法的发展,基于路或建筑物表面等多个线面特征在相邻观测航线上提取特征并校正的方法也取得了一定的进展[10-12]。本文尝试使用参数加权最小二乘方法使点云拟合在一个平面上检校系统误差,这与文献[13—14]使用多个平面处理点云间位置关系的方法有所不同。
1 简单模式下的机载激光雷达测深系统定位方程图 1为简单模式下的机载激光测深雷达扫描方式及所在的激光扫描坐标系,在大气和海水中的传播距离分别为ρ1和ρ2,α为激光固定发射角,θsz为扫描角。图 2中激光光线方向矢量l1穿过水面折射后变为l2,θ1、θ2分别为激光方向矢量与海水表面法向量p的夹角。na、nw分别为大气和海水的折射率。海底激光脚点在激光扫描坐标系中的坐标可表示为
式中
根据文献[15]提出的折射定律矢量表达方程,得到静止海水表面的折射公式
则通过静止海面后激光光线的方向矢量为
激光测深系统将姿态测量系统、差分GPS接收机、激光扫描仪集成于一体。理想情况下,姿态测量系统坐标轴与激光扫描仪坐标轴平行,实际情况下却存在安置角常差,通常称为视准轴偏角。同样,差分GPS接收机天线中心与激光扫描中心也存在偏移量,其偏差称为偏心分量。设视准轴偏角对应的3个旋转角参数分别为航向偏角θix、滚动偏角θiz、俯仰偏角θiy,激光扫描中心至GPS/INS中心的偏心分量为ΔXil=xilyilzilT,则INS坐标系中激光脚点坐标为(详细的各个坐标系间关系的描述参考文献[16—17],坐标系间的数学转换关系可参考[18—20])
式中,R(·)为旋转角矩阵。
最后,根据GPS/INS测定的飞机在成图坐标系中的姿态变化,即航向角ψvx、俯仰角ψvy、滚动角ψvz,和GPS/INS中心在成图坐标系的位置ΔX=xIyIzIT,得到激光脚点在成图坐标系中的坐标为
因此,简单模式下圆扫描式机载激光测深系统的定位模型可表述为
点云可以看成是激光光线与平面交会产生的,海底激光点云经过了海表和海底两个平面,下面分别推导激光光线与海表面、海底面交会的数学表达式。
2.1 海表面激光点云模型激光光线和海表面的交会如图 3所示,激光光线于时间t=t0与海表面交会,[xLyLzL]为直线上任意一点,[aLbLcL]T为激光光线的方向,平面法向量为[aPbPcP]T,[xPyPzP]为平面上任意点,则激光光线和海表面的数学表达式分别为
解算得到
对于机载激光测深系统,假设某个时刻海表面法向量为Gn,飞机激光发射点的海表面垂点坐标为Npt,激光发射点坐标为Lpt,激光方向矢量为Lu,则激光光线与海表面交会的时间有如下表达
在球面坐标系下,海表面法向量可表示为
设飞机开始飞行时其正下方位于海表面的点为原点,则
由激光脚点定位模型得到激光发射点坐标
激光光线方向矢量为
根据以上表达式即可解算激光光线与海表面交会的位置。
2.2 水底激光点云模型激光穿透海面发生折射,通过折射公式(7)得到水中传输时的激光方向。由于激光与水面的交点已经确定,因此可以得到海水中激光光线直线方程,此时只需确定模型海底面的数学表达式即可得到海底激光脚点位置。解算步骤与海表面点云模型类似。
3 参数加权最小二乘平差模型对于加权最小二乘,观测值残差平方和可表达为
式中,r=[r1r2r3…rn]为观测值残差;Cl=diag(σ2l1, σ2l2, σ2l3, …, σ2ln);σ2li为第i个观测值对应的方差。
加入未知数加权改正数,得到
式中,δ为未知数改正值向量;Cx为未知数先验方差组成的对角矩阵。加入未知数权值,可以使最小二乘算法先调整权重高的参数,也使得最小二乘算法避免向无效数值收敛。
待估参数X与观测值L之间组成的函数关系式为
由泰勒级数展开得到一阶线性逼近
式中,g为模型初始值;r为观测值改正值;δ为未知数改正值;D为观测值对应的偏导矩阵;A为未知数对应的偏导矩阵。
由矩阵分块求解最终得到未知数δ和观测值残差r
在平面区域获取激光点云时,传感器的系统误差和随机误差使得本该共面的激光点云不共面,使用最小二乘模型将不共面的激光点云拟合到单个平面上,从而降低系统误差和随机误差的影响,进而纠正点云位置。
如图 4所示,n是平面法向量,xo是激光点云,xp为平面上的任意点,若所有激光点均位于平面上,则任意点xo与xp相减构成的向量与平面法向量n的向量积为0,得到
实际上,因为系统误差和随机误差造成在平面区域获取的激光点云不共面,使得式(26)不为0,需要使用检校算法纠正点云位置。
为了方便分析,这里将偏心分量设为0,设海平面是静止水平的,将海平面作为检校平面,本文定位模型进一步简化,即
对于平面法向量,使用球坐标系,n=[ρcosθsinϕρsinθsinϕ-ρcosϕ]T,设ρ=1,则n=[cosθsinϕsinθsinϕ-cosϕ]T,未知数由3个降为2个,则
式中,xp=[xpypzp]T;观测值向量是L=[ρ1θszψvxψvyψvzxIyIzI]T;未知数向量是X=[θixθiyθizφpθpxpypzp]T;其他参数参照前文。
根据参数加权最小二乘模型,下面需要分别求出模型初始值g,观测值对应的偏导矩阵D,未知数对应的偏导矩阵A。
模型初始值g为
观测值偏导矩阵D为
式中,D中非对角线部分每个元素0均为1×8零矩阵。
矩阵A为f对未知数X=[θixθiyθizpθpxpypzp]T的偏导矩阵。
观测值权矩阵为
式中,Cl中非对角线部分每个元素0均为1×8零矩阵。
未知数权矩阵为
根据式(29)—式(33)和式(24),即可解算未知数改正数。
假设存在n组观测值,则式(24)中(DClDT)-1为n×n矩阵的求逆,在本文后面模拟试验中,n最少等于1000,这对于计算机是巨大的工作量。由于D、Cl矩阵中均含有大量0元素,可以尝试将计算过程简化,由于Cl为对角线矩阵,得到
式中,DCl中非对角线部分每个元素0均为1×8零矩阵,Ti为对角线矩阵,不再详列。
式中,Si的计算式不再详列。
由于DClDT为对角线矩阵,则(DClDT)-1=diag(1/S1, 1/S2, …, 1/Sn)
使用上式对AT(DClDT)-1的计算过程优化,效率有质的提升。
5 点云模拟及检校模型数值计算下面使用文章中的定位模型、点云模拟模型仿真点云的生成过程(图 5—图 7)。以Optech公司的CZMIL机载激光测深系统[21-23]为例:航高H为400m,水文模式下激光脉冲发射频率Rprf为10kHz,棱镜上表面斜坡坡度为39.16°,脉冲固定发射角度θ为20°,扫描仪转速Rrps为27Hz(圈/s),飞行速度为140Kts(海里/小时)。空气折射率为1.0003,海水折射率为1.334,棱镜折射率为1.462,模拟点云500圈。本文主要对视准轴偏角进行探讨,假设不存在其他系统误差,分别使3个视准轴偏角设为1°,模拟在飞机平直飞行、航高从400m逐渐升至500m、飞机3个姿态角分别由0°至30°逐渐变化的点云分布。对于每个视准轴偏角,应有5对点云分布图,限于篇幅原因,只列出部分典型情况的点云分布图(所谓“典型”是指能够将本应共面的点云较大幅度的非共面分散,未列出的其他情况点云仍共面或非共面分散幅度较小),其中左图为立体图,右图为航向视角的YZ平面图,蓝色为存在视准轴偏角时获取的点云,红色为俯仰偏角未知时计算得到的点云。
基于以上点云的模拟得到以下检校思路:在存在视准轴俯仰偏角时,保持平直飞行,即可使本应共面的点云不共面,从而探测该偏角;在存在视准轴滚动偏角时,检校该偏角需要使航向姿态角发生变化;在存在视准轴航向偏角时,检校该偏角需要滚动姿态角发生变化。总的来说,在检校飞行时同时使航向角和滚动角发生动态变化即可基于本文检校思路解算3个视准轴偏角。
基于上述思路,采用文中激光点云定位方程和点云模拟模型,使飞行航向角和翻滚角同时由0°向10°均匀发生变化,为了防止最小二乘解算方程数据量过大,激光发射频率设为54Hz(扫描仪转速Rrps为27Hz,避免非整数近似操作),使3个视准轴偏角都设为5°,得到图 8所示的点云轨迹,共1000个点,同时记录1000对观测值。其中左图为立体图,右图为航向视角的YZ平面图,蓝色为存在视准轴偏角时获取的点云,红色为俯仰偏角未知时计算得到的点云。
参照国际上广泛应用的相关仪器系统指标,对各观测值和待求未知数标准差赋值。激光测距仪测距标准差设为0.01,扫描角标准差设为0.002,GPS/INS定位定姿系统的位置标准差和姿态标准差都采用广泛使用且精度较高的POS/AVTM 610后处理之后的系统参数(表 1)。视准轴偏角θix、θiy、θiz标准差使用文献[17—18]论文里的典型数值(表 2)。
根据上述获取的1000对观测值及对应的标准差(每对观测值标准差相同)和未知数标准差,采用文中参数加权最小二乘平差模型解算,结果如表 3所示。将频率分别增大到2倍、4倍、10倍、20倍、200倍,视准轴偏差标准差统计结果如表 4所示。将飞机姿态角航向角、滚动角变化幅度增大至2倍,视准轴偏差标准差统计结果如表 5所示。分别同方向和反方向飞行两条航带,航带间隔为1000m,视准轴偏差标准差统计结果如表 6所示。
频率 | 点数 | σθ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 |
姿态角变化 | 点数 | σθ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 |
航带 | 点数 | σθ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 分析与讨论
表 4中,随着激光发射频率的增加,点密度的提升能够降低随机误差的影响,逐渐提高视准轴偏角检校精度。由于CZMIL水文模式设计激光频率为10kHz,此时检校精度颇高,但计算时间也大幅度增加至约27min,因此在精度要求不是很高的情况下,可以重采样激光点云从而降低点云数以提升检校效率。
表 5中,在飞机航向姿态角变化至20°时,视准轴滚动偏角标准差有效降低,航向偏角和俯仰偏角的标准差几乎不变;同样在飞机滚动姿态角变化至20°时,视准轴航向偏角标准差有效降低,俯仰偏角和滚动偏角的标准差几乎不变,刚好验证了点云模拟时关于检校方式的结论。
表 6中,在飞机反方向飞行两次时,航向姿态角类似于突然改变了180°,视准轴滚动偏角检校精度有了质的提升,视准轴航向偏角和俯仰偏角类似于增加了点密度,检校精度也得到了提升。在飞机同方向飞行两次时,较单航带时点数量增多且分散范围更大,视准轴滚动偏角检校精度也得到了大幅提升,说明滚动偏角对于航向角的变化和点云的位置分散情况均相当敏感,视准轴航向偏角和俯仰偏角变化程度类似于反方向飞行两次的情况。
因此为了提高视准轴偏角检校的精度,需要在计算量允许的情况下保持较高的点密度,增加飞机航向姿态角和滚动姿态角变化幅度的同时飞行多个航带。
7 结束语本文在推导了简单模式下机载激光测深系统定位模型和模拟模型后,使用参数加权最小二乘方法使点云拟合在一个平面上检校机载激光测深系统误差,通过模拟试验初步验证了思路的正确性,并通过多项对比试验,对检校实施提出了有意义的建议。不足的地方有:只讨论了视准轴偏差存在的情况,对于其他可能的系统误差需要从定位模型开始继续深入研究,检校方法类似文中思路,只是未知数更多,对于各种系统误差之间的相关性也需要有针对性地深入研究。
[1] | 江月松. 机载GPS、姿态和激光扫描测距集成定位系统的精确定位方程、误差分析与精度评估[J]. 遥感学报 , 2001, 5 (4) : 241–247. JIANG Yuesong. A Rigorous Positioning Equation and It's Error Analysis and Precision Evaluation for Integrated Positioning System of Airborne GPS, INS and Laser Scanning Ranging[J]. Journal of Remote Sensing , 2001, 5 (4) : 241 –247. |
[2] | 翟国君, 王克平, 刘玉红. 机载激光测深技术[J]. 海洋测绘 , 2014, 34 (2) : 72–75. ZHAI Guojun, WANG Keping, LIU Yuhong. Technology of Airborne Laser Bathymetry[J]. Hydrographic Surveying and Charting , 2014, 34 (2) : 72 –75. |
[3] | 黄卫军, 李松. 基于表面回波的机载激光测深系统的最佳扫描方案[J]. 激光杂志 , 2001, 22 (6) : 54–56. HUANG Weijun, LI Song. The Perfect Scanning Configurations of Laser-based Airborne Hydrographic System[J]. Laser Journal , 2001, 22 (6) : 54 –56. |
[4] | 任来平, 赵俊生, 翟国君, 等. 机载激光测深海面扫描轨迹计算与分析[J]. 武汉大学学报(信息科学版) , 2002, 27 (2) : 138–142. REN Laiping, ZHAO Junsheng, ZHAI Guojun, et al. Scanning-track Computation and Analysis for Airborne Laser Depth Sounding[J]. Geomatics and Information Science of Wuhan University , 2002, 27 (2) : 138 –142. |
[5] | 王建军, 刘吉东. 影响机载激光扫描点云精度的测量误差因素分析及其影响大小排序[J]. 中国激光 , 2014, 41 (4) : 0414001. WANG Jianjun, LIU Jidong. Analysis and Sorting of Impacts of Measurement Errors on Positioning Accuracy of Laser Point Cloud Obtained from Airborne Laser Scanning[J]. Chinese Journal of Lasers , 2014, 41 (4) : 0414001 . DOI:10.3788/CJL |
[6] | 王丽英, 宋伟东. 机载LiDAR点云航带平差方法研究[J]. 武汉大学学报(信息科学版) , 2012, 37 (7) : 814–817. WANG Liying, SONG Weidong. Airborne LiDAR Point Cloud Strip Adjustment Method[J]. Geomatics and Information Science of Wuhan University , 2012, 37 (7) : 814 –817. |
[7] | 王致华, 张爱武, 王书民, 等. 基于重叠航带的机载激光雷达系统检校[J]. 中国激光 , 2014, 41 (2) : 0214003. WANG Zhihua, ZHANG Aiwu, WANG Shumin, et al. Airborne Radar Calibration System Based on the Overlap Strip[J]. Chinese Journal of Lasers , 2014, 41 (2) : 0214003 . DOI:10.3788/CJL |
[8] | 张靖, 江万寿. 基于虚拟连接点模型的机载LiDAR系统安置误差自检校[J]. 测绘学报 , 2011, 40 (6) : 762–769. ZHANG Jing, JIANG Wanshou. Self-calibration of LiDAR System Mounting Biases Using Virtual Tie Point Model[J]. Acta Geodaetica et Cartographica Sinica , 2011, 40 (6) : 762 –769. |
[9] | 张靖, 江万寿, 姜三. 基于虚拟点模型的机载LiDAR系统自动检校方法[J]. 测绘学报 , 2013, 42 (3) : 389–396. ZHANG Jing, JIANG Wanshou, JIANG San. Automated Airborne LiDAR System Calibration Using Virtual Tie Point Model[J]. Acta Geodaetica et Cartographica Sinica , 2013, 42 (3) : 389 –396. |
[10] | HABIB A F, BANG K I, SHIN S W, et al.LiDAR System Self-calibration Using Planar Patches from Photogrammetric Data[C]//Proceedings of the Fifth International Symposium on Mobile Mapping Technology.Padua:[s.n.], 2007. |
[11] | HABIB A F, AL-DURGHAM M, KERSTING A P, et al.Error Budget of LiDAR Systems and Quality Control of the Derived Point Cloud[C]//Proceedings of the 21st International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences.Beijing, China:ISPRS, 2008:203-209. |
[12] | HABIB A F, KERSTING A P, RUIFANG Z, et al.LiDAR Strip Adjustment Using Conjugate Linear Features in Overlapping Strips[C]//Proceedings of the 21st International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences.Beijing, China:ISPRS, 2008:385-390. |
[13] | VOSSELMAN G, DIJKMAN S.3D Building Model Reconstruction from Point Clouds and Ground Plans[C]//Proceedings of the International Archives of the Photogrammetry and Remote Sensing.Annapolis, MD:ISPRS, 2001:37-43. http://cn.bing.com/academic/profile?id=2141191968&encoded=0&v=paper_preview&mkt=zh-cn |
[14] | ALHARTHY A, BETHEL J.Detailed Building Reconstruction from Airborne Laser Data Using a Moving Surface Method[C]//Proceedings of the 17th International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences.Istanbul, Turkey:ISPRS, 2004:213-218. |
[15] | GLASSNER A S. An Introduction to Ray Tracing[M]. London: Academic Press, 1989 : 166 -190. |
[16] | FRIESS P.Toward a Rigorous Methodology for Airborne Laser Mapping[C]//Proceedings of the International Calibration and Orientation Workshop.Castelldefels, Spain:[s.n.], 2006:7-11. http://www.academia.edu/15379288/Rigorous_approach_to_bore-sight_self-calibration_in_airborne_laser_scanning |
[17] | SKALOUD J, SCHAER P.Towards Automated LiDAR Boresight Self-calibration[C]//Proceedings of the 5th International Symposium on Mobile Mapping Technology.Padova:MMS, 2007. http://www.isprs.org/proceedings/XXXVI/5-C55/ |
[18] | CHANG Guobin. On Least-squares Solution to 3D Similarity Transformation Problem Under Gauss-Helmert Model[J]. Journal of Geodesy , 2015, 89 (6) : 573 –576. DOI:10.1007/s00190-015-0799-z |
[19] | FANG Xing. Weighted Total Least-squares with Constraints:A Universal Formula for Geodetic Symmetrical Transformations[J]. Journal of Geodesy , 2015, 89 (5) : 459 –469. DOI:10.1007/s00190-015-0790-8 |
[20] | FANG Xing. A Total Least Squares Solution for Geodetic Datum Transformations[J]. Acta Geodaetica et Geophysica , 2014, 49 (2) : 189 –207. DOI:10.1007/s40328-014-0046-8 |
[21] | TUELL G, BARBOR K, WOZENCRAFT J.Overview of the Coastal Zone Mapping and Imaging LiDAR (CZMIL):A New Multisensor Airborne Mapping System for the U.S.Army Corps of Engineers[C]//Proceedings of SPIE 7695, Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVI.Orlando, Florida:SPIE, 2010:76950R. |
[22] | FUCHS E, TUELL G.Conceptual Design of the CZMIL Data Acquisition System (DAS):Integrating a New Bathymetric LiDAR with a Commercial Spectrometer and Metric Camera for Coastal Mapping Applications[C]//Proceedings of SPIE 7695, Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVI.Orlando, Florida:SPIE, 2010:76950U. |
[23] | PARK J Y, TUELL G.Conceptual Design of the CZMIL Data Processing System (DPS):Algorithms and Software for Fusing LiDAR, Hyperspectral Data, and Digital Images[C]//Proceedings of SPIE 7695, Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVI.Orlando, Florida:SPIE, 2010:7695010. |