[发明专利]一种空地结合的潮间带一体化测绘方法

权利要求书

1.一种空地结合的潮间带一体化测绘方法，其特征在于，包括测量装置，测量装置包括无人机和载体车，无人机上搭载激光探测与测量系统，载体车上设置有激光扫描仪、单波束测深仪、惯性导航系统和卫星定位系统；

测绘方法包括以下步骤：

步骤1：原始数据采集下载后，检查各测量仪器观测数据的完整性、正确性，进行观测数据文件的格式转换，进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，通过点云数据预处理，形成标准LAS数据文件；

步骤2：进行点云时空融合，完成从仪器坐标系到地心地固坐标系的转换；

步骤3：通过获取单波束点云与激光扫描仪点云、激光扫描仪点云与机载LIDAR点云数据之间的点、线、面等几何特征来建立对应关系，求解旋转矩阵和平移量进行点云匹配。

2.根据权利要求1所述的一种空地结合的潮间带一体化测绘方法，其特征在于，所述步骤1包括以下子步骤：

步骤1.1：利用硬件同步控制器进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，硬件同步控制器的时间保持与卫星定位系统的时间同步；硬件同步控制器控制各测量仪器的触发及监测其反馈信号，信号的传输都采用差分方式；

步骤1.2：利用IE进行基准站和流动站数据结算，输出POS数据；

步骤1.3：利用VSursPROCESS软件进行点云数据预处理，形成标准LAS数据文件。

3.根据权利要求1所述的一种空地结合的潮间带一体化测绘方法，其特征在于，所述步骤2包括以下子步骤：

步骤2.1：从激光扫描仪坐标系到载体坐标系转换；

根据标定的公共点计算坐标系转换需要的的6参数：l_x、l_y、l_z、ω、κ；

l_x、l_y、l_z分别为坐标转换过程中点坐标沿载体坐标系x轴、y轴、z轴的平移量；

ω、κ分别为坐标转换过程中点坐标绕载体坐标系x轴、y轴、z轴的旋转角度；

设激光扫描仪坐标系下点坐标为在载体坐标系下的坐标为则有，

其中，代表从激光扫描仪坐标系到载体坐标系转换的旋转矩阵，R_x(ω)、R_z(κ)分别为绕载体坐标系x轴、y轴、z轴旋转的旋转矩阵；

步骤2.2：从载体坐标系到站心坐标系转换；

惯导坐标系前进方向为y轴、向右为x轴、向上为z轴；

惯导记录姿态角：侧滚角Roll、俯仰角Pitch、偏航角Heading；

侧滚角Roll：惯导x轴与水平方向之间的夹角，载体右侧向下为正；

俯仰角Pitch：惯导y轴与水平方向之间的夹角，载体向上为正；

偏航角Heading：惯导前进方向，即xy平面与正北方向之间的夹角，顺时针为正；

设Roll、Pitch、Heading分别为r、p、y；设点在站心坐标系下的坐标为载体坐标系下坐标向站心坐标系转换；

①先绕z轴旋转y；

②再绕x轴旋转p；

③最后绕y轴旋转r；

则有，

其中，代表从载体坐标系到站心坐标系转换的旋转矩阵，R_x(p)、R_y(r)、R_z(y)分别为绕站心坐标系x轴、y轴、z轴旋转的旋转矩阵；

步骤2.3：从站心坐标系到地心地固坐标系转换；

站心坐标系原点在WGS84下的经纬度分别为L和B；

设扫描点在地心地固坐标系下的坐标为站心坐标系下的坐标转换为WGS84坐标系下的坐标：

①先绕x轴旋转

②再绕z轴旋转

③最后将站心坐标系原点平移到WGS84坐标系原点；

则有，

其中，代表从站心坐标系到地心地固坐标系转换的旋转矩阵；

其中，N代表卯酉圈曲率半径，h代表大地高，a、b代表椭球长、短半径；

得

为站心坐标系原点在地心地固坐标系下的空间直角坐标；

最后得，

4.根据权利要求1所述的一种空地结合的潮间带一体化测绘方法，其特征在于，所述步骤3包括：

首先介绍一下所需要的理论基础：

平面拟合原理

已知平面方程：a′x+b′y+c′z＝d (6)

且a′²+b′²+c′²＝1，d≥0，(x y z)是平面上的任意点；

其中，(a′ b′ c′)代表垂直于平面且远离坐标原点方向的方向余弦；

同样的，它也表示垂直于平面远离坐标原点的单位法向量的分量；

d表示了平面和坐标原点间的垂直距离；

最小二乘算法基本原理

最小二乘法的思想就是估算一个最佳数学模型，求得最佳模型参数使估算模型点的值和实际观察值的差值的平方和最小；

在二维数据中，假设给出模型函数g(x_i)，g(x_i)与真实值(x_i，y_i)的误差可用三种形式表述：

①误差最大值，即其中n为点总个数；

②误差绝对值之和，即

③误差平方和，即

最小二乘算法的原理便是寻找一个最佳参数模型，对于给定的数据(x_i，y_i)(i＝0，1…n)，误差平方和最小；

最小二乘平面拟合算法

基于最小二乘法的平面拟合算法是在三维点云数据的基础上进行的拟合算法，以公式(6)为三维模型，确定a′，b′，c′，d四个参数使三维模型的估算值和测量值之间的差值平方和E²最小；由于E²是非负数，因此它存在极小值，E²对每个参数的偏导数为零，利用偏导数为零构建的方程组便可求解出a′，b′，c′，d，并使估算值和测量值之间的差值平方和E²最小；

假设点云数据的点集为{p_i|p_i∈oxyz，i＝1，2，3...，N}，其中x_i和y_i设为无误差变化的自变量，z_i是包含误差的因变量，则z和x、y的函数相关性可以假设为函数表达式：

z＝f(x，y；p，q，r)＝p+qx+ry (7)

且其中的p，q，r可以通过最小二乘算法得到；

根据公式(6)、(7)，令d＝c′p，a′＝-qc′，b′＝-rc′且则公式(7)可以描述为平面方程；

根据上述定义f函数的公式，假如自变量的值为x_i和y_i，被测量值z应该为p+qx_i+ry_i，然而因为误差的存在测量值z为z_i，则第i个点的测量误差为：

e_i＝f(x_i，y_i；p，q，r)-z_i (8)

根据最小二乘的原理，p，q，r应该满足公式(8)的平方和最小，即E²最小：

若使得E²最小，应满足：

即：

有：

根据公式(11)可以求解出p，q，r，同样的，公式(6)中a′，b′，c′，d四个参数也就可通过d＝c′p，a′＝-qc′，b′＝-rc′，求出；

求解旋转矩阵

在上述算法的计算过程中，提取3对以上对应的平面特征量后，需要求取其各自的平面法向量；

假设求解的一对平面为a₁x+b₁y+z＝d₁和a₂x+b₂y+z＝d₂，则其法向量对可表示为P＝(a₁，b₁，-1)、Q＝(a₂，b₂，-1)，法向量的旋转矩阵R满足P＝RQ，采用罗德里格旋转法求解；上述P，Q矩阵，其夹角为：

P，Q矩阵的叉乘为：

若已知单位向量利用法向量将其旋转θ角度后，根据罗德里格旋转法求解其三维旋转矩阵为：

即可求出的三维变换矩阵C₁；

求解平移量

利用三维变换矩阵C₁可以对点云数据进行旋转，相应的平面特征量对相互平行，在两片点云数据中寻找一对以上的特征点，或利用对应面的的质心坐标的差值，计算点云数据的平移量即：