CN103106339A - 同步航空影像辅助的机载激光点云误差改正方法 - Google Patents

Abstract

一种同步航空影像辅助的机载激光点云误差改正方法，结合LiDAR点云和POS数据辅助进行同步航空影像高精度连接点自动匹配；然后结合少量地面控制点进行航空影像POS辅助光束法区域网平差，解算POS数据的***误差；最后基于LiDAR严格定位模型，结合每个LiDAR点的GPS时间、LiDAR点物方坐标、POS导航数据文件以及POS***误差，重新解算准确的LiDAR点云三维坐标。本发明解决了多尺度、多视角影像高精度自动匹配问题；采用少量地面控制点进行同步航空影像POS辅助光束法区域网平差，准确解算POS数据***误差；采用机载LiDAR定位严格模型和求解出的POS***误差改善LiDAR点云坐标精度。

Description

同步航空影像辅助的机载激光点云误差改正方法

技术领域

本发明属于测绘科学与技术领域，涉及一种同步航空影像辅助的机载激光点云误差改正方法，主要应用于高精度DEM/DSM生产、城市三维建模、真正射影像生产等领域。

背景技术

机载激光探测与测距（Airborne Light Detecting and Ranging,机载LiDAR）是一项集成激光扫描仪和定位定姿***（Positioning and Orientation System,POS）的三维地形数据快速获取技术。机载LiDAR设备以一定的频率向地面发射激光束并接收地物反射激光信号，根据发射信号与反射信号的时间差计算扫描仪与地物的距离，再结合POS设备获取的平台当前的姿态计算地物的三维坐标。这一过程称为直接对地定位，获取的具有三维坐标的密集地面点称为LiDAR点云，在一条航带飞行中获取的点云称为LiDAR点云条带。在LiDAR点云数据获取过程中，常存在多种误差源，如激光扫描仪测距误差、扫描镜控制误差、激光扫描仪与POS设备时间同步误差、激光扫描仪与POS设备的偏心以及安置角误差、POS观测值***误差等。其中POS观测值***误差为主要误差源，以上误差会直接影响所获取的LiDAR点云坐标精度。在LiDAR数据处理的实际生产中发现，点云的高程误差一般在0.3m至1.0m之间，特别是在POS设备精度不高以及GPS信号不理想时，点云坐标误差尤其不能忽视，常导致点云条带之间无法良好吻合、生成的DEM精度无法满足需求。

目前在实际生产和学术研究中解决LiDAR点云坐标精度问题时，主要有以下几种方法：①人工选定某一点云条带作为基准，在该条带内平坦处外业测量若干个高程控制点，将基准条带的高程值加减一个常量使其与控制点基本吻合，然后通过将其余点云条带的高程值加减一个常量使其与基准条带在高程方向上尽量吻合；②选定基准条带，在点云条带间重叠区内寻找同名的人字形屋顶面、边缘线特征、角点特征等作为控制数据采用三维空间相似变换改正其余条带变形误差，或是选择地面面片作为虚拟连接点采用二次多项式模型进行点云条带区域网平差；③基于LiDAR三角测量定位模型，在点云条带间手动或自动匹配同名平面片，以此作为控制数据采用最小二乘方法解算激光扫描仪与POS设备之间的安置角误差，在此基础上改正由安置角误差导致的点云坐标误差。①与②采用的点云改正模型并不符合点云条带误差规律，仅能对点云条带误差进行粗略改正。方法③虽基于LiDAR定位模型，但在改正点云时仅改正激光扫描仪与POS设备安置角误差，而LiDAR设备中激光扫描仪与POS设备为高度集成，二者间的安置误差对定位精度的影响一般可忽略。在实际数据获取中，LiDAR点云的主要误差源应为POS***的定位定姿误差，表现为航带漂移及线性漂移误差。由于目前的LiDAR***均可配备中等幅面数码相机，可在采集点云的同时获取高分辨率彩色航空影像，可结合少量地面控制点，将POS数据作为带权观测值引入航空影像光束法区域网平差中，解求POS数据的***误差，进而结合LiDAR严格定位模型改正点云坐标误差。

发明内容

本发明的目的是提供一种同步航空影像辅助的机载激光点云误差改正方法。

本发明的技术方案为一种同步航空影像辅助的机载激光点云误差改正方法，包括以下步骤，

步骤1，结合POS数据及LiDAR点云辅助进行同步航空影像连接点自动匹配；

步骤2，采用若干地面控制点，引入POS数据作为带权观测值，进行同步航空影像POS辅助光束法区域网平差，解求POS数据的***误差；所述地面控制点至少包括测区四角各一个控制点和测区中心的一个控制点；

步骤3，对于LiDAR点云数据，根据LiDAR点获取时的GPS时间、POS导航数据文件以及LiDAR点的物方坐标，恢复LiDAR点与LiDAR***的位置关系，结合在步骤2同步航空影像POS辅助光束法区域网平差中检校所得POS数据的***误差，对该时刻LiDAR***的位置与姿态进行改正，并通过重新解算LiDAR点坐标进行LiDAR点云坐标误差改正，实现LiDAR点云精度优化。

而且，步骤1实现方式为，由LiDAR点云数据内插生成数字表面模型，在各航空影像上采用Harris算子提取特征点，对于每一特征点根据POS数据和数字表面模型预测在各重叠影像上的同名点位，并采用仿射变换纠正影像间尺度和角度差异，最后依次采用相关系数匹配和最小二乘匹配获取特征点的准确同名点。

而且，步骤2所述同步航空影像POS辅助光束法区域网平差的公式如下，

V_X=A_Xq+Bx+CS_X-l,P

V_G=A_Gq+C_GS_G-l_G,P_G     (4)

V_I=A_Iq+C_IS_I-l_I,P_I

其中，V_X,V_G,V_I分别为像点坐标、GPS摄站坐标及IMU姿态角观测的改正数向量；q和x分别为影像外方位元素及加密点地面坐标改正数向量；s_X,s_G,s_I分别为像点***误差、GPS及IMU观测值***误差改正数向量；A_X,A_G,A_I,B,C,C_G,C_I分别为各未知数向量的系数矩阵；l,l_G,l_I分别为误差方程常数项向量；P,P_G,P_I分别为像点坐标、GPS摄站坐标及IMU姿态角的权阵。

而且，步骤3中LiDAR点云坐标误差改正，实现方式为对LiDAR点云中每一LiDAR点执行以下步骤，

①对于每一个LiDAR点，根据获取LiDAR点时的GPS时间在POS导航数据文件中内插出该时刻POS设备的位置与姿态角，获得t时刻内插获得的GPS天线中心坐标[X_G(t),Y_G(t),Z_G(t)]^T及t时刻内插获得的IMU角度观测值ω_I(t),φ_I(t),κ_I(t)组成的旋转矩阵R_I(t)；

②结合该LiDAR点的物方坐标[Y_P(t),Y_P(t),Z_P(t)]^T，采用下式计算当前时刻该LiDAR点在POS设备坐标系下的坐标[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T；

$\left[\begin{array}{c}{X}_{P\left(t\right)}\\ Y_{P\left(t\right)}\\ Z_{P\left(t\right)}\end{array}\right]=\left[\begin{array}{c}{X}_{G\left(t\right)}\\ Y_{G\left(t\right)}\\ Z_{G\left(t\right)}\end{array}\right]+R_{I\left(t\right)}\left[\begin{array}{c}{X}_{\mathrm{LPOS}\left(t\right)}\\ Y_{\mathrm{LPOS}\left(t\right)}\\ Z_{\mathrm{LPOS}\left(t\right)}\end{array}\right]---\left(6\right)$

其中，[Y_P(t),Y_P(t),Z_P(t)]^T为t时刻获取的LiDAR点的物方空间坐标；[X_G(t),Y_G(t),Z_G(t)]^T为t时刻内插获得的GPS天线中心坐标；[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T为当前时刻LiDAR点在POS设备坐标系下的坐标，R_I(t)为t时刻内插获得的IMU角度观测值组成的旋转矩阵；

③根据在同步航空影像区域网平差中解算获得的POS数据线元素和角元素***误差值，采用下式对POS数据线元素和角元素进行误差改正，

$\left[\begin{array}{c}{X}_{G\left(t\right)}^{\′}\\ Y_{G\left(t\right)}^{\′}\\ Z_{G\left(t\right)}^{\′}\end{array}\right]=\left[\begin{array}{c}{X}_{G\left(t\right)}\\ Y_{G\left(t\right)}\\ Z_{G\left(t\right)}\end{array}\right]-\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]-(t-t_0)\left[\begin{array}{c}{b}_X\\ b_Y\\ b_z\end{array}\right]---\left(8\right)$

其中，[X_G(t),Y_G(t),Z_G(t)]^T为改正前的POS数据的线元素，ω_I(t),φ_I(t),κ_I(t)为改正前的POS数据的角元素，(a_X，a_Y，a_Z),(b_X，b_Y，b_Z)为GPS观测值随航带的常数及漂移误差，t为当前影像的摄影时刻，t₀为参考时刻，

为IMU观测值，

为IMU观测值随航带的常数及漂移误差，X′_G(t),Y’_G(t),Z′_G(t)为改正后的POS数据的线元素，为改正后的POS数据的角元素；

④采用改正后的POS线元素[X′_G(t),Y’_G(t),Z’_G(t)]^T及角元素

组成的旋转矩阵R’_I(t)，结合②中计算出的[X_LPOS(t),Y_LPPOS(t),Z_LPOS(t)]^T采用公式（6）重新解算该LiDAR点的物方坐标[Y’_P(t),Y’_P(t),Z’_P(t)]^T。

本发明的优点是采用LiDAR点云和POS数据辅助航空影像连接点自动匹配，无需人工进行航带排列，解决了多尺度、多视角、高分辨率影像匹配问题；将POS观测值作为带权观测值引入航空影像光束法区域网平差中，在联合平差中仅需少量地面控制点即可准确解求POS***误差、影像外方位元素及加密点地面坐标；基于LiDAR对地定位严格模型，深入分析LiDAR点坐标误差产生机理，结合同步影像区域网平差中解算的POS数据***误差值，实现LiDAR点坐标误差改正。与现有方法相比，本方法理论上最严密、最符合LiDAR点云坐标误差规律。该方法实现了机载激光点云精度优化，为机载LiDAR在高精度DEM/DSM生产、正射影像生产及数字城市三维建模等方面的应用提供了有力保障。

附图说明

图1是本发明实施例的流程图。

图2是本发明实施例的POS数据及LiDAR点云辅助同步航空影像连接点自动匹配示意图。

具体实施方式

具体实施时，本发明技术方案可采用计算机软件技术实现自动运行。以下结合附图和实施例详细说明本发明技术方案。

参见图1，本发明实施例提供基于同步航空影像区域网平差的机载激光点云精度优化方法。对同步航空影像进行POS辅助光束法区域网平差，解算POS数据***误差，然后基于LiDAR对地定位严格模型逐个改正LiDAR点的坐标误差，具体实现方法包含以下步骤：

步骤1，进行POS数据及LiDAR点云辅助同步航空影像连接点自动匹配。

实施例由LiDAR点云数据内插生成数字表面模型（DSM），在各航空影像上采用Harris算子提取特征点，对于每一特征点根据POS数据和DSM预测其在各重叠影像上的同名点位，并采用仿射变换纠正影像间尺度和角度差异获取匹配搜索区，最后依次采用相关系数匹配和最小二乘匹配获取特征点的准确同名点。Harris算子为现有技术。

本方法采用POS数据和由LiDAR点云生成的数字表面模型（DSM）辅助进行航空影像连接点自动匹配，能够自动对同名点位进行预测，并基于仿射变换对搜索区进行重采样以消除影像间分辨率、角度差异，可解决多分辨率、多视角的交叉航带影像匹配问题。如图2所示，S，S₀，S_i为三幅影像，实施例的自动匹配步骤如下：

首先，将LiDAR点云数据按照一定的采样间隔（一般取点云平均间距）生成数字表面模型DSM；

其次，对每张航空影像，提取Harris特征点；

然后，逐一将航空影像当作参考影像与其重叠影像进行匹配，方法如下：

①对于某一张参考影像（如图2中影像S），首先在影像上打上3×3的格网，依次将每个格网点依据共线条件方程投影到DSM上获取格网点的概略地面坐标，并将概略地面坐标依据共线条件方程反投影到其他各航空影像上，通过判断反投影像点是否在有效像幅内获得与该参考影像重叠的所有影像。

②对于参考影像S上的某一特征点p，依据共线条件方程将其投影到DSM上，获得该点对应的概略物方坐标P(X_P,Y_P,Z_P)，再依次将点P依据共线条件方程反投影到各重叠影像上得到p点的同名点预测点位p_i(i=0,1…m,i为各重叠影像号，图2中m=1)。

③根据设置的搜索区大小，在参考影像S上以p点为中心取搜索区大小的影像块（如图2所示a,b,c,d四点为影像块的角点），将影像块四个角点投影到高程平面Z_P上，获取四点的物方坐标，并分别反投影到各重叠影像上获得像点（a_i,b_i,c_i,d_i）（i=0,1…m），根据（a,b,c,d）坐标与相应的（a_i,b_i,c_i,d_i）坐标计算参考影像窗口与重叠影像i上窗口的仿射变换参数，在重叠影像i上内插获得纠正的搜索窗口。然后在参考影像与纠正后的搜索窗口之间采用相关系数匹配和最小二乘匹配获得特征点p的同名点。

最后，整合所有影像的匹配结果，获得整个测区的影像加密点。

步骤2，结合少量地面控制点的POS辅助同步航空影像光束法区域网平差：结合少量地面控制点（至少测区四角、测区中心各一个，共五个控制点），引入POS数据作为带权观测值，进行公式（4）所示同步航空影像POS辅助光束法区域网平差，准确解求POS数据***误差。

实施例将POS数据作为带权观测值，采用考虑航带漂移、线性漂移的一次多项式拟合POS数据的线元素观测值及角元素观测值***误差，同时考虑到摄影中心与POS***中心的偏心分量误差、相机视准轴与IMU视准轴的偏置角误差以及影像残余的镜头畸变误差，结合少量地面控制点在区域网平差中一并解求POS***误差值与其它***误差值。

为便于实施参考起见，提供航空影像POS辅助光束法区域网平差详细说明如下：

以摄影测量中共线方程为模型，以影像连接点像方坐标、POS线元素和角元素为观测值进行联合区域网平差。实施例采用

转角***，任意一个像点与对应的影像外方位元素及地面点物方空间坐标间满足共线方程：

$x-x_0=-f\frac{a_1(X-X_S)+b_1(Y-Y_S)+c_1(Z-Z_S)}{a_3(X-X_S)+b_3(Y-Y_S)+c_3(Z-Z_S)}$      （1）

$y-y_0=-f\frac{a_2(X-X_S)+b_2(Y-Y_S)+c_2(Z-Z_S)}{a_3(X-X_S)+b_3(Y-Y_S)+c_3(Z-Z_S)}$

其中，x,y为像点的像平面坐标；x₀,y₀,f为影像内方位元素；X,Y,Z为地面点的空间坐标；X_S,Y_S,Z_S为外方位线元素；a_i,b_i,c_i(i=1,2,3)为外方位角元素

组成的旋转矩阵的各元素。

航空数码摄影的POS数据在后处理时一般均会对摄影中心与POS***中心的偏心分量及相机视准轴与IMU视准轴的偏置角进行改正，不过大量实践表明改正后仍然存在少量GPS残余偏心分量和IMU视准轴误差。实际数据处理时，可以认为POS提供的外方位元素除存在一定的残余偏心分量误差和残余视准轴偏置角误差外，其余***误差主要由GPS和IMU随时间的***漂移引起，已有研究表明在不超过15min飞行时间的航带内，GPS漂移误差随飞行时间t呈线性变化，IMU姿态角测量值漂移误差也具有类似规律。考虑GPS天线偏心残余分量、IMU视准轴残余误差的影响，并利用一次多项式表达GPS/IMU观测值随航带的常数误差和漂移误差，即可建立光束法区域网平差时POS观测值的***误差改正模型。

GPS观测值(X_G,Y_G,Z_G)的误差模型如下：

$\left[\begin{array}{c}{X}_G\\ Y_G\\ Z_G\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]+R\left[\begin{array}{c}u\\ v\\ w\end{array}\right]+\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]+(t-t_0)\left[\begin{array}{c}{b}_X\\ b_Y\\ b_z\end{array}\right]---\left(2\right)$

其中(u,v,w)为GPS天线分量，R为外方位角元素矩阵，(a_X，a_Y，a_Z),(b_X，b_Y，b_Z)为GPS观测值随航带的***漂移误差，t为当前影像的摄影时刻，t₀为参考时刻（一般取航带中心影像的摄影时刻）。

角元素的计算涉及IMU观测值

构成的旋转矩阵R_IMU、相机视准轴与IMU视准轴的偏置角

构成的旋转矩阵R_{MIS_C}及外方位角元素

构成的旋转矩阵R，三个旋转矩阵均为

转角***，关系为R_IMU=R_{MIS_C}R。考虑航带常数误差

漂移误差

IMU观测值的误差模型如下：

其中t与t₀与公式（2）中意义相同；a_i,b_i,c_i(i=1,2,3)为矩阵R_{MIS_C}R的各元素。

根据以上函数模型，区域网平差的未知数包括影像外方位元素、同名点地面坐标、摄影中心与POS***中心的偏心分量误差、相机视准轴与IMU视准轴的偏置角误差、GPS和IMU观测值随航带的常数误差及漂移误差等。如果考虑同名像点也存在一定的***误差，则还需添加相应的像点***误差改正模型，如12参数的Ebner模型及44参数的Gruen模型。区域网平差时联合公式（2）、（3）与公式（1），建立联合平差模型，对其线性化后得到如下的POS辅助光束法区域网平差误差方程式：

V_X=A_Xq+Bx+CS_X-l,P

V_G=A_Gq+C_GS_G-l_G,P_G      (4)

V_I=A_Iq+C_IS_I-l_I,P_I

其中V_X,V_G,V_I分别为像点坐标、GPS摄站坐标及IMU姿态角观测的改正数向量；q和x分别为影像外方位元素及加密点地面坐标改正数向量；s_X,s_G,s_I分别为像点***误差、GPS及IMU观测值***误差改正数向量；A_X,A_G,A_I,B,C,C_G,C_I分别为各未知数向量的系数矩阵；l,l_G,l_I分别为误差方程常数项向量；P,P_G,P_I分别为像点坐标、GPS摄站坐标及IMU姿态角的权阵。

设测区共有m条航线，N幅影像，n个加密点，则不考虑像点***误差时区域网平差的未知数总数为：3+9×m+3+9×m＋6×N+3×n。根据三类观测值的测量精度分别赋予相应权值，即可用最小二乘法解算各未知数的平差值。

步骤3，LiDAR点云坐标误差改正：对于LiDAR点云数据中的LiDAR点逐个进行坐标改正。根据LiDAR点获取时的GPS时间、POS导航数据文件以及LiDAR点的物方坐标，恢复LiDAR点与POS***的位置关系，根据在同步航空影像区域网平差中检校出的POS***误差值，对该时刻POS***的位置与姿态进行改正，并采用新的位置与姿态参数根据LiDAR对地定位严格模型重新解算LiDAR点的坐标，实现LiDAR点云精度优化的目的。

本方法涉及的基于同步航空影像POS辅助光束法区域网平差的机载激光点云精度优化方法原理如下。基于如下LiDAR对地定位严格模型进行点云误差改正：

$\left[\begin{array}{c}{X}_{P\left(t\right)}\\ Y_{P\left(t\right)}\\ Z_{P\left(t\right)}\end{array}\right]=\left[\begin{array}{c}{X}_{G\left(t\right)}\\ Y_{G\left(t\right)}\\ Z_{G\left(t\right)}\end{array}\right]+R_{I\left(t\right)}\left[\begin{array}{c}{X}_L\\ Y_L\\ Z_L\end{array}\right]+R_{I\left(t\right)}{R}_{\mathrm{MIS}\_L}{R}_{R\left(t\right)}\left[\begin{array}{c}0\\ 0\\ d\left(t\right)\end{array}\right]---\left(5\right)$

其中，[Y_P(t),Y_P(t),Z_P(t)]^T为t时刻获取的LiDAR点的物方空间坐标；[X_G(t),Y_G(t),Z_G(t)]^T为t时刻内插获得的GPS天线中心（即POS***中心）坐标；R_I(t)为t时刻内插获得的IMU角度观测值组成的旋转矩阵；[X_L,Y_L,Z_L]^T为激光扫描仪中心在POS坐标系下的偏心分量；R_{MIS_L}为激光扫描仪视准轴与IMU视准轴的偏置角组成的旋转矩阵；R_R(t)为t时刻激光束在扫描仪坐标系下的旋转角组成的旋转矩阵；d(t)为t时刻扫描仪测量的地面LiDAR点与扫描仪的距离。

由于LiDAR设备中激光扫描仪和POS设备为紧密集成，且会定期对二者的相对位置关系进行检校，故在公式（5）中，激光扫描仪与POS***的偏心分量[X_L,Y_L,Z_L]^T及偏置角组成的旋转矩阵R_{MIS_L}均可被认为不存在误差，此外激光扫描仪本身的测距和测角误差也可忽略，即对于一个LiDAR点可认为R_R(t)[0,0,d(t)]^T亦不存在误差。基于以上分析，公式（5）可简化成如下形式：

$\left[\begin{array}{c}{X}_{P\left(t\right)}\\ Y_{P\left(t\right)}\\ Z_{P\left(t\right)}\end{array}\right]=\left[\begin{array}{c}{X}_{G\left(t\right)}\\ Y_{G\left(t\right)}\\ Z_{G\left(t\right)}\end{array}\right]+R_{I\left(t\right)}\left[\begin{array}{c}{X}_{\mathrm{LPOS}\left(t\right)}\\ Y_{\mathrm{LPOS}\left(t\right)}\\ Z_{\mathrm{LPOS}\left(t\right)}\end{array}\right]---\left(6\right)$

其中，[Y_P(t),Y_P(t),Z_P(t)]^T为t时刻获取的LiDAR点的物方空间坐标；[X_G(t),Y_G(t),Z_G(t)]^T为t时刻内插获得的GPS天线中心坐标；[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T为当前时刻LiDAR点在POS设备坐标系下的坐标，R_I(t)为t时刻内插获得的IMU角度观测值组成的旋转矩阵；

当前时刻该LiDAR点在POS设备坐标系下的坐标[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T如公式（7）所示，在点云坐标改正中对于某一特定点为常量矩阵，故本方法在进行点云改正时无需获知激光扫描仪与POS***的偏心分量与偏置角。

$\left[\begin{array}{c}{X}_{\mathrm{LPOS}\left(t\right)}\\ Y_{\mathrm{LPOS}\left(t\right)}\\ Z_{\mathrm{LPOS}\left(t\right)}\end{array}\right]=\left(\begin{array}{c}\left[\begin{array}{c}{X}_L\\ Y_L\\ Z_L\end{array}\right]+R_{\mathrm{MIS}\_L}{R}_{R\left(t\right)}\left[\begin{array}{c}0\\ 0\\ d\left(t\right)\end{array}\right]\end{array}\right)---\left(7\right)$

基于以上分析设计了LiDAR点云坐标误差改正方法，步骤如下：

①对于每一个LiDAR点，根据其获取时的GPS时间在POS导航数据文件中内插出该时刻POS设备的位置与姿态角，获得t时刻内插获得的GPS天线中心坐标[X_G(t),Y_G(t),Z_G(t)]^T及t时刻内插获得的IMU角度观测值ω_I(t),φ_I(t),κ_I(t)组成的旋转矩阵R_I(t)；

②结合该LiDAR点的物方坐标[Y_P(t),Y_P(t),Z_P(t)]^T，采用公式（6）计算当前时刻该LiDAR点在POS设备坐标系下的坐标[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T；

③根据在同步航空影像区域网平差中解算获得的POS数据线元素和角元素***误差值，采用公式（8）和（9）对POS数据线元素和角元素进行误差改正。

$\left[\begin{array}{c}{X}_{G\left(t\right)}^{\′}\\ Y_{G\left(t\right)}^{\′}\\ Z_{G\left(t\right)}^{\′}\end{array}\right]=\left[\begin{array}{c}{X}_{G\left(t\right)}\\ Y_{G\left(t\right)}\\ Z_{G\left(t\right)}\end{array}\right]-\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]-(t-t_0)\left[\begin{array}{c}{b}_X\\ b_Y\\ b_z\end{array}\right]---\left(8\right)$

其中，[X_G(t),Y_G(t),Z_G(t)]^T为改正前的POS数据的线元素，ω_I(t),φ_I(t),κ_I(t)为改正前的POS数据的角元素，(a_X，a_Y,a_Z),(b_X，b_Y，b_Z)为GPS观测值随航带的常数及漂移误差，t为当前影像的摄影时刻，t₀为参考时刻，

为IMU观测值,

为IMU观测值随航带的常数及漂移误差，X′_G(t),Y’_G(t),Z′_G(t)为改正后的POS数据的线元素，

为改正后的POS数据的角元素。

④采用改正后的POS线元素X′_G(t),Y’_G(t),Z′_G(t)及角元素

组成的旋转矩阵R’_I(t)，结合②中计算出的[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T采用公式（6）重新解算该LiDAR点的物方坐标[Y’_P(t),Y’_P(t),Z′_P(t)]^T；

对每一LiDAR点均进行①-④步操作，直至所有点坐标改正完毕。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (3)

1.一种同步航空影像辅助的机载激光点云误差改正方法，其特征在于：包括以下步骤，

步骤1，结合POS数据及LiDAR点云辅助进行同步航空影像连接点自动匹配；

步骤2，采用若干地面控制点，引入POS数据作为带权观测值，进行同步航空影像POS辅助光束法区域网平差，解求POS数据的***误差；所述地面控制点至少包括测区四角各一个控制点和测区中心的一个控制点；

步骤3，对于LiDAR点云数据，根据LiDAR点获取时的GPS时间、POS导航数据文件以及LiDAR点的物方坐标，恢复LiDAR点与LiDAR***的位置关系，结合在步骤2同步航空影像POS辅助光束法区域网平差中检校所得POS数据的***误差，对该时刻LiDAR***的位置与姿态进行改正，并通过重新解算LiDAR点坐标进行LiDAR点云坐标误差改正，实现LiDAR点云精度优化。

2.根据权利要求1所述同步航空影像辅助的机载激光点云误差改正方法，其特征在于：步骤1实现方式为，由LiDAR点云数据内插生成数字表面模型，在各航空影像上采用Harris算子提取特征点，对于每一特征点根据POS数据和数字表面模型预测在各重叠影像上的同名点位，并采用仿射变换纠正影像间尺度和角度差异，最后依次采用相关系数匹配和最小二乘匹配获取特征点的准确同名点。

3.根据权利要求2所述同步航空影像辅助的机载激光点云误差改正方法，其特征在于：步骤2所述同步航空影像POS辅助光束法区域网平差的公式如下，

V_X=A_Xq+Bx+CS_X-l,P

V_G=A_Gq+C_GS_G-l_G,P_G    (4)

V_I=A_Iq+C_IS_I-l_I,P_I

其中，V_X,V_G,V_I分别为像点坐标、GPS摄站坐标及IMU姿态角观测的改正数向量；q和x分别为影像外方位元素及加密点地面坐标改正数向量；s_X,s_G,s_I分别为像点***误差、GPS及IMU观测值***误差改正数向量；A_X,A_G,A_I,B,C,C_G,C_I分别为各未知数向量的系数矩阵；l,l_G,l_I分别为误差方程常数项向量；P,P_G,P_I分别为像点坐标、GPS摄站坐标及IMU姿态角的权阵。4.根据权利要求3所述同步航空影像辅助的机载激光点云误差改正方法，其特征在于：步骤3中LiDAR点云坐标误差改正，实现方式为对LiDAR点云中每一LiDAR点执行以下步骤，

①对于每一个LiDAR点，根据获取LiDAR点时的GPS时间在POS导航数据文件中内插出该时刻POS设备的位置与姿态角，获得t时刻内插获得的GPS天线中心坐标[X_G(t),Y_G(t),Z_G(t)]^T及t时刻内插获得的IMU角度观测值ω_I(t),φ_I(t),κ_I(t)组成的旋转矩阵R_I(t)；

②结合该LiDAR点的物方坐标[Y_P(t),Y_P(t),Z_P(t)]^T，采用下式计算当前时刻该LiDAR点在POS设备坐标系下的坐标[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T；

$\left[\begin{array}{c}{X}_{P\left(t\right)}\\ Y_{P\left(t\right)}\\ Z_{P\left(t\right)}\end{array}\right]=\left[\begin{array}{c}{X}_{G\left(t\right)}\\ Y_{G\left(t\right)}\\ Z_{G\left(t\right)}\end{array}\right]+R_{I\left(t\right)}\left[\begin{array}{c}{X}_{\mathrm{LPOS}\left(t\right)}\\ Y_{\mathrm{LPOS}\left(t\right)}\\ Z_{\mathrm{LPOS}\left(t\right)}\end{array}\right]---\left(6\right)$

其中，[Y_P(t),Y_P(t),Z_P(t)]^T为t时刻获取的LiDAR点的物方空间坐标；[X_G(t),Y_G(t),Z_G(t)]^T为t时刻内插获得的GPS天线中心坐标；[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T为当前时刻LiDAR点在POS设备坐标系下的坐标，R_I(t)为t时刻内插获得的IMU角度观测值组成的旋转矩阵；

③根据在同步航空影像区域网平差中解算获得的POS数据线元素和角元素***误差值，采用下式对POS数据线元素和角元素进行误差改正，

$\left[\begin{array}{c}{X}_{G\left(t\right)}^{\′}\\ Y_{G\left(t\right)}^{\′}\\ Z_{G\left(t\right)}^{\′}\end{array}\right]=\left[\begin{array}{c}{X}_{G\left(t\right)}\\ Y_{G\left(t\right)}\\ Z_{G\left(t\right)}\end{array}\right]-\left[\begin{array}{c}{a}_X\\ a_Y\\ a_Z\end{array}\right]-(t-t_0)\left[\begin{array}{c}{b}_X\\ b_Y\\ b_z\end{array}\right]---\left(8\right)$

其中，[X_G(t),Y_G(t),Z_G(t)]^T为改正前的POS数据的线元素，ω_I(t),φ_I(t),κ_I(t)为改正前的POS数据的角元素，(a_X，a_Ya_Z),(b_X，b_Y，b_Z)为GPS观测值随航带的常数及漂移误差，t为当前影像的摄影时刻，t₀为参考时刻，

为IMU观测值，

为IMU观测值随航带的常数及漂移误差，X′_G(t),Y’_G(t),Z′_G(t)为改正后的POS数据的线元素，

为改正后的POS数据的角元素；

④采用改正后的POS线元素[X′_G(t),X’_G(t),Z’_G(t)]^T及角元素

组成的旋转矩阵R’_I(t)，结合②中计算出的[X_LPOS(t),Y_LPOS(t),Z_LPOS(t)]^T采用公式（6）重新解算该LiDAR点的物方坐标[Y’_P(t),Y’_P(t),Z′_P(t)]^T。

Images