使用Spherical线性内插方法解算移动LiDAR点云 | ![]() |
2. 中国石油天然气管道工程有限公司,河北 廊坊,065000
2. China Petroleum Pipeline Engineering Corporation, Langfang 065000, China
国家石油天然气管网集团有限公司的成立标志着深化油气体制改革迈出关键一步,公司主要负责全国油气干线管道、部分储气调峰设施的投资建设、干线管道互联互通及与社会管道联通。面对中国油气能源进口和输送的强大需求,针对管道设计和运营中地理信息数据采集的难点,掌握移动激光雷达(light detection and ranging,LiDAR)集成技术,有助于加快油气管道中线的断面测量及管道设计进度,为后期智慧管网建设提供高质量的站场和阀室三维地理信息数据。随着LiDAR技术的普及,硬件制造商将激光传感器设计得既小巧又实用,这使其应用范围更广泛。针对不同的应用需求,在集成LiDAR系统时,激光发射器的安装位置和角度越发灵活多样。为了适应LiDAR定制化的市场要求,工程研究人员有必要掌握惯性测量单元(inertial measurement unit,IMU)、GPS和激光发射器的硬件集成方案,并编写程序将这三者的数据联合解算成点云[1]。
本文使用Spherical线性内插方法解决IMU与激光扫描角度的万向节锁问题,对机载和背包Li- DAR进行实验,首先将定位测姿系统(position and orientation system,POS)与激光发射器刚性固连[2, 3],将集成设备安装在无人机和一台推车上,采集完实验数据后,编写代码将POS和激光数据进行联合解算,得到机载和背包点云数据[4]。
1 研究方法 1.1 激光器与POS之间的坐标关系实验集成的设备有POS设备和激光器。LiDAR在集成过程中涉及到成图坐标系统、IMU坐标系统、激光器坐标系统[5]。
激光打到的地面点I可以转换为地面坐标rIm:
$ \boldsymbol{r}_I^m=\boldsymbol{r}_b^m(t)+\boldsymbol{R}_b^m(t) \boldsymbol{r}_l^b+\boldsymbol{R}_b^m(t) \boldsymbol{R}_l^b \boldsymbol{r}_I^l(t) $ | (1) |
式中,t表示激光发射时间;rbm表示IMU原点的地面坐标;Rbm表示姿态信息;rlb、Rlb分别表示激光坐标系统与IMU坐标系统间的偏移量和旋转矩阵;rIl表示激光坐标系统中的点云坐标。
在激光器坐标系统里,原点定义在激光束发射点,Z轴沿着激光发射器的旋转轴,每一束激光发射时都有一个垂直角β、水平角α和激光发射点到地物之间的距离ρ。
$ \boldsymbol{r}_I^l(t)=\left(\begin{array}{l} x \\ y \\ z \end{array}\right)=\left(\begin{array}{c} \rho(t) \cos \beta(t) \cos \alpha(t) \\ \rho(t) \cos \beta(t) \sin \alpha(t) \\ \rho(t) \sin \beta(t) \end{array}\right) $ | (2) |
相对于IMU,激光发射器安装时会有一定角度,要引入虚拟激光器坐标系统L',虚拟激光器坐标系统和IMU坐标系统基本保持对齐,所以将式(1)修改为:
$ \boldsymbol{r}_I^m=\boldsymbol{r}_b^m(t)+\boldsymbol{R}_b^m(t) \boldsymbol{r}_l^b+\boldsymbol{R}_b^m(t) \boldsymbol{R}_{L^{\prime}}^b \boldsymbol{R}_L^{L^{\prime}} \boldsymbol{r}_I^l(t) $ | (3) |
式中,RLL'是激光器相对于IMU的安装旋转角度。
以上公式可从数学上解释本实验硬件集成系统解算点云的过程[8]。
万向节锁是LiDAR系统解算点云的一个重要问题。在三维中常用的欧拉角坐标定向系统用绕3个轴旋转的角度来表示物体的朝向(Rx、Ry、Rz),注意3个轴是针对物体坐标系的。例如,飞行器从初始位置先绕Xl轴转30°,然后绕Yl轴转90°,最后绕Zl轴转−40°,对应地,用欧拉角坐标来跟踪是(30°、90°、−40°)。若飞行器又绕Xl轴转1°,那么欧拉坐标应该变成(31°、90°、−40°),使用这个坐标,飞行器坐标位置无法对应。实际上,使用3个量来表示三维空间朝向的系统都会遇到这个问题,除非用4个量来描述方向,如四元数方法[9, 10]。
1.2 四元数方法一般3个向量的旋转矩阵可表示如下:
$ \boldsymbol{R}_\omega \boldsymbol{R}_{\varphi} \boldsymbol{R}_k=\boldsymbol{R}=\left[\begin{array}{lll} r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \end{array}\right] $ | (4) |
现在用4个向量q=q0+qxi+qyj+qzk来描述激光器发射激光的方向,四元数旋转方向量将3个量的旋转矩阵转变为4个量的旋转矩阵,就可以避免万向锁的问题。
1.3 Spherical线性内插假设在q1和q2两个激光发射瞬间,POS记录下了激光发射器的位置姿态信息,这样很容易计算q1和q2的激光发射位置和角度。在q1和q2之间还有很多个qi的激光发射位置和角度信息需要计算[11]。
Spherical线性内插描述qi如下:
$ \boldsymbol{q}_i=\frac{\sin \theta_2}{\sin \theta} \boldsymbol{q}_1+\frac{\sin \theta_1}{\sin \theta} \boldsymbol{q}_2 $ | (5) |
式中,θ是q1和q2的夹角;θ1是q1和qi的夹角;θ2是q2和qi的夹角。将qi的四元数反算为激光传感器的三维旋转矩阵:
$ \left\{\begin{array}{l} r_{11}=\boldsymbol{q}_x^2+\boldsymbol{q}_0^2-\boldsymbol{q}_z^2-\boldsymbol{q}_y^2 \\ r_{12}=2 \boldsymbol{q}_x \boldsymbol{q}_y-2 \boldsymbol{q}_0 \boldsymbol{q}_z \\ r_{13}=2 \boldsymbol{q}_x \boldsymbol{q}_z+2 \boldsymbol{q}_0 \boldsymbol{q}_y \\ r_{21}=2 \boldsymbol{q}_x \boldsymbol{q}_y+2 \boldsymbol{q}_0 \boldsymbol{q}_z \\ r_{22}=\boldsymbol{q}_y^2-\boldsymbol{q}_z^2+\boldsymbol{q}_0^2-\boldsymbol{q}_x^2 \\ r_{23}=2 \boldsymbol{q}_y \boldsymbol{q}_z-2 \boldsymbol{q}_0 \boldsymbol{q}_x \\ r_{31}=2 \boldsymbol{q}_x \boldsymbol{q}_z-2 \boldsymbol{q}_0 \boldsymbol{q}_y \\ r_{32}=2 \boldsymbol{q}_y \boldsymbol{q}_z+2 \boldsymbol{q}_0 \boldsymbol{q}_x \\ r_{33}=\boldsymbol{q}_z^2-\boldsymbol{q}_y^2-\boldsymbol{q}_x^2+\boldsymbol{q}_0^2 \end{array}\right. $ | (6) |
所以有:
$ \left\{\begin{array}{l} \varphi=\arcsin \left(r_{13}\right) \\ \omega=\arctan \left(\frac{-r_{23}}{\cos \varphi}, \frac{r_{33}}{\cos \varphi}\right) \times \frac{180}{\pi} \\ k=\arctan \left(\frac{-r_{12}}{\cos \varphi}, \frac{r_{11}}{\cos \varphi}\right) \times \frac{180}{\pi} \end{array}\right. $ | (7) |
式中,φ、ω、k分别是qi的3个方向角度。由此可以内插出qi的激光发射位置和角度信息,再计算激光打到地面的大地坐标。
2 实验与分析将集成的LiDAR系统安装在无人机和小推车上采集数据。本文使用的是旋翼无人机,采用前后往返飞行方式初始化IMU后,进入航线采集数据。激光传感器与POS之间的安装角度是90°,停车场的照明灯和车都清晰可见,点密度均匀。
图 1展示了使用简单线性内插方法与Spherical优化方法得到的点云数据,图 1(b)中车的轮廓更细腻,这表明使用Spherical方法得到的点云质量更好。
![]() |
图 1 使用两种方法得到的点云数据 Fig.1 Point Cloud Data Obtained by Two Methods |
将集成的LiDAR系统放置在小推车上,沿着路从房屋一侧推到另一侧,由图 2可知,抽稀后路边房子的点云和房屋两侧的点云密度均匀、结构完整,使用Spherical线性内插方法可以很好地为每一束发射激光插入位置姿态信息,并解算出点云大地坐标。集成的机载和背包LiDAR系统都属于移动LiDAR,可直接解算测区的整体点云信息,缩短工期。
![]() |
图 2 背包LiDAR系统的点云 Fig.2 Point Clouds of Backpack LiDAR System |
在石油天然气管道设计阶段需要准确测量管道中线的断面,如果管道必须经过山区林地,那就给测量工作提出了新要求。例如,已经开工的中俄管线在山区林地使用的航空数据采集方式无法有效穿透植被测量管道中线断面信息,采用移动背包LiDAR能在测量困难地区减少工作量。
油气管道每隔一段距离会有站场和阀室控制石油和天然气的压力等,在油气管网设计阶段要测绘站场和阀室选址地的1∶500地形图,中国智慧油气管网的建设需要在GIS系统里建立站场和阀室的真实三维模型,移动LiDAR是现阶段快速测绘小面积1∶500地形图和三维建模的有效手段。
3 结束语针对油气管网多样化的需求,激光器的安装必将灵活多变。激光器小型化的趋势使其使用范围不局限于测量领域。在激光器和POS联合解算点云的过程中,机载和背包LiDAR都无法回避为发射激光内插位置姿态信息的问题,Spherical线性内插方法可以很好地解决大视场角激光传感器与POS之间的定姿定位问题,确保点云解算程序流畅运行。此项技术的研究将可以为中国油气管道的站场、阀室和中线测量提供技术支撑。
[1] |
李奇, 马洪超. 基于激光雷达波形数据的点云生产[J]. 测绘学报, 2008, 37(3): 349-354. DOI:10.3321/j.issn:1001-1595.2008.03.014 |
[2] |
周伟, 李奇, 李畅. 利用激光扫描技术监测大型古建筑变形的研究[J]. 测绘通报, 2012(4): 52-54. |
[3] |
邬建伟, 马洪超, 李奇. 顾及语义的机载LiDAR点云格网化方法[J]. 测绘科学技术学报, 2008, 25(2): 87-89. |
[4] |
Li Q, Ma H C, Wu J W, et al. Filter Algorithm for Airborne LiDAR Data[C]. MIPPR 2007: Multispectral Image Processing, Bellingham, USA, 2007
|
[5] |
Zhou W, Chen F L, Guo H D, et al. UAV Laser Scanning Technology: A Potential Cost-Effective Tool for Micro-topography Detection over Wooded Areas for Archaeological Prospection[J]. International Journal of Digital Earth, 2020, 13(11): 1279-1301. DOI:10.1080/17538947.2019.1711209 |
[6] |
Chiang K W, Tsai M L, Naser E S, et al. New Calibration Method Using Low Cost MEM IMUs to Verify the Performance of UAV-Borne MMS Payloads[J]. Sensors, 2015, 15(3): 6560-6585. DOI:10.3390/s150306560 |
[7] |
Kersting A P, Habib A. A Comparative Analysis Between Rigorous and Approximate Approaches for LiDAR System Calibration[J]. Journal of the Korean Society of Surveying, Geodesy, Photogrammetry and Cartography, 2012, 30(6_2): 593-605. DOI:10.7848/ksgpc.2012.30.6-2.593 |
[8] |
Yan W Y, Shaker A, Habib A, et al. Improving Classification Accuracy of Airborne LiDAR Intensity Data by Geometric Calibration and Radiometric Correction[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 67: 35-44. DOI:10.1016/j.isprsjprs.2011.10.005 |
[9] |
Li Q, Zhou W, Li C. Use of Airborne LiDAR to Estimate Forest Stand Characteristics[J]. IOP Conference Series: Earth and Environmental Science, 2014, 17: 012252. |
[10] |
Cramer M, Stallmann D. System Calibration for Direct Georeferencing[C]. Commission Ⅲ Symposium"Photogrammetric Computer Vision", Graz, Austria, 2002
|
[11] |
Bäumker M, Heimes F J. New Calibration and Computing Method for Direct Georeferencing of Image and Scanner Data Using the Position and Angular Data of an Hybrid Inertial Navigation System[C]. Proceeding of OEEPE Workshop on Integrated Sensor Orientation, Hannover, Germany, 2002
|