WO2020151213A1

WO2020151213A1 - Air and ground combined intertidal zone integrated mapping method

Info

Publication number: WO2020151213A1
Application number: PCT/CN2019/098770
Authority: WO
Inventors: 卢秀山; 石波; 田茂义; 李国玉; 陈超; 孙海超; 李丁硕
Original assignee: 青岛秀山移动测量有限公司
Priority date: 2019-01-22
Filing date: 2019-08-01
Publication date: 2020-07-30
Also published as: CN109631863A

Abstract

An air and ground combined intertidal zone integrated mapping method. After original data is acquired and downloaded, point cloud spatiotemporal fusion is performed to complete transformation from an instrument coordinate system to an earth-centered earth-fixed coordinate system; a correspondence is established by obtaining point, line, plane and other geometric characteristics between single beam point cloud and laser scanner point cloud data and between laser scanner point cloud and airborne laser detection and measurement system point cloud data, so as to solve a rotation matrix and a translation amount for point cloud matching. An airborne laser detection and measurement system is integrated with a towing measurement system, and thus underwater topography and overwater topography data of an intertidal zone can be obtained, thereby solving the problem of intertidal zone topography measurement. The two systems can perfectly cover a measurement blind area, so as to achieve integrated and seamless intertidal zone measurement.

Description

一种空地结合的潮间带一体化测绘方法An integrated surveying and mapping method of intertidal zone combined with space and ground

技术领域Technical field

本发明属于海洋测绘技术领域，具体涉及以一种空地结合的方式进行潮间带区域水上和水下地形快速获取的测绘方法。The invention belongs to the technical field of marine surveying and mapping, and particularly relates to a surveying and mapping method for rapidly acquiring water and underwater topography in an intertidal zone area in an air-ground combination.

背景技术Background technique

潮间带属于海岸带的一部分，易受到潮汐的影响，测量条件比较困难。当前，海岸带测绘技术主要有人工实地测量、船载测量和航测遥感等方式。自然海岸带通常存在礁石、滩涂、湿地等特殊地理环境，以传统人工实地测量和船载测量方式难以高效地进行作业，甚至存在危及人身安全的作业风险，无法保障海岸带基础地理信息的快速更新。航测遥感技术手段具有范围广、频度高和实时的优点，被广泛应用于海洋测绘，但是由于受到卫星重复周期、海岸带云雨天气、分辨率低等影响，难以满足动态监测和水下地形的高精度测量需求，对潮水覆盖区域尤其是大面积浅滩和宽阔滩涂区域无能为力。The intertidal zone is a part of the coastal zone and is susceptible to the influence of tides, and the measurement conditions are more difficult. Currently, coastal surveying and mapping technologies mainly include manual field surveys, shipborne surveys, and aerial surveys and remote sensing. Natural coastal zones usually have special geographical environments such as reefs, tidal flats, wetlands, etc. It is difficult to perform operations efficiently with traditional manual field surveys and shipborne surveys, and there are even operational risks that endanger personal safety, and it is impossible to ensure the rapid update of basic geographic information in the coastal zone. . Aerial survey and remote sensing technology has the advantages of wide range, high frequency and real-time. It is widely used in marine surveying and mapping. However, it is difficult to meet the requirements of dynamic monitoring and underwater terrain due to satellite repetition cycles, coastal cloud and rain weather, and low resolution. The demand for high-precision measurement is powerless for areas covered by tidal water, especially large shoals and wide tidal flat areas.

船载多传感器水岸一体化综合测量技术的发展，为海岸带一体化测量提供了依据。但当海岸带地形以滩涂为主且潮差较小时，由于测量船和测深仪存在吃水问题使得基于船载激光测量无法对潮间带进行扫测，出现潮间带的测量“盲区”。The development of ship-borne multi-sensor integrated water-shore integrated measurement technology provides a basis for integrated coastal measurement. However, when the topography of the coastal zone is dominated by tidal flats and the tidal range is small, due to the draught problem of the survey ship and the echo sounder, it is impossible to scan the intertidal zone based on the shipborne laser measurement, and the measurement "blind zone" of the intertidal zone appears.

机载LiDAR技术具有高便捷、无控制、高精度、高密集、高效率、高分辨率、多回波等特点，作为精确、快速地获取地面三维数据的工具已得到广泛的认同，但由于潮间带地形特殊性，激光信号受地物遮挡、水面反射等影响，也会出现测量“盲区”。Airborne LiDAR technology has the characteristics of high convenience, no control, high precision, high density, high efficiency, high resolution, multiple echoes, etc. It has been widely recognized as a tool for accurately and quickly acquiring 3D ground data. Due to the particular terrain of the interzone, the laser signal is affected by the obstruction of ground objects and the reflection on the water surface, and measurement "blind areas" will also appear.

发明内容Summary of the invention

针对现有方法中存在的问题，本发明提出一种空地结合的潮间带一体化测绘方法，机载激光探测与测量***结合拖曳一体化测量***，解决了潮间带常见的测量问题，具有良好的效果。In view of the problems existing in the existing methods, the present invention proposes an integrated air-ground surveying and mapping method for intertidal zones. The airborne laser detection and measurement system is combined with the towing integrated measurement system, which solves the common measurement problems in the intertidal zone. good effect.

本发明采用以下的技术方案：The present invention adopts the following technical solutions:

一种空地结合的潮间带一体化测绘方法，包括测量装置，测量装置包括无人机和载体车，无人机上搭载激光探测与测量***，载体车上设置有激光扫描仪、单波束测深仪、惯性导航***和卫星定位***；An air-ground integrated intertidal zone surveying and mapping method includes a measuring device. The measuring device includes an unmanned aerial vehicle and a carrier vehicle. The unmanned aerial vehicle is equipped with a laser detection and measurement system. The carrier vehicle is equipped with a laser scanner and single beam sounding Instrument, inertial navigation system and satellite positioning system;

测绘方法包括以下步骤：The surveying method includes the following steps:

步骤1：原始数据采集下载后，检查各测量仪器观测数据的完整性、正确性，进行观测数据文件的格式转换，进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，通过点云数据预处理，形成标准LAS数据文件；Step 1: After the original data is collected and downloaded, check the completeness and correctness of the observation data of each measuring instrument, convert the format of the observation data file, perform high-precision time synchronization of each measuring instrument, and integrate the time and space of the data of each measuring instrument. Point cloud data preprocessing to form standard LAS data files;

步骤2：进行点云时空融合，完成从仪器坐标系到地心地固坐标系的转换；Step 2: Perform point cloud space-time fusion and complete the conversion from the instrument coordinate system to the geocentric ground-fixed coordinate system;

步骤3：通过获取单波束点云与激光扫描仪点云、激光扫描仪点云与机载LIDAR点云数据之间的点、线、面等几何特征来建立对应关系，求解旋转矩阵和平移量进行点云匹配。Step 3: Establish the correspondence relationship between the single-beam point cloud and the laser scanner point cloud, the laser scanner point cloud and the airborne LIDAR point cloud data, such as geometric features such as points, lines, and surfaces, and solve the rotation matrix and translation amount Perform point cloud matching.

优选地，所述步骤1包括以下子步骤：Preferably, the step 1 includes the following sub-steps:

步骤1.1：利用硬件同步控制器进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，硬件同步控制器的时间保持与卫星定位***的时间同步；硬件同步控制器控制各测量仪器的触发及监测其反馈信号，信号的传输都采用差分方式；Step 1.1: Use the hardware synchronization controller to synchronize the high-precision time of each measuring instrument, the time and space of the data of each measuring instrument are integrated, and the time of the hardware synchronization controller is synchronized with the time of the satellite positioning system; the hardware synchronization controller controls each measuring instrument Trigger and monitor its feedback signal, signal transmission adopts differential mode;

步骤1.2：利用IE进行基准站和流动站数据结算，输出POS数据；Step 1.2: Use IE to settle the data of base station and rover, and output POS data;

步骤1.3：利用VSursPROCESS软件进行点云数据预处理，形成标准LAS数据文件。Step 1.3: Use VsursPROCESS software to preprocess point cloud data to form a standard LAS data file.

优选地，所述步骤2包括以下子步骤：Preferably, the step 2 includes the following sub-steps:

步骤2.1：从激光扫描仪坐标系到载体坐标系转换；Step 2.1: Convert from the laser scanner coordinate system to the carrier coordinate system;

根据标定的公共点计算坐标系转换需要的的6参数：l _x、l _y、l _z、ω、

κ； Calculate the 6 parameters needed for coordinate system conversion according to the calibrated common point: l _x , l _y , l _z , ω,

κ;

l _x、l _y、l _z分别为坐标转换过程中点坐标沿载体坐标系x轴、y轴、z轴的平移量； l _x , l _y , and l _z are the translational amounts of the point coordinates along the x-axis, y-axis and z-axis of the carrier coordinate system during the coordinate conversion process;

ω、

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

κ is the rotation angle of the point coordinate around the x-axis, y-axis and z-axis of the carrier coordinate system during the coordinate conversion process;

设激光扫描仪坐标系下点坐标为

在载体坐标系下的坐标为

则有， Let the coordinates of the point in the laser scanner coordinate system be

The coordinates in the carrier coordinate system are

Then there is,

其中，

代表从激光扫描仪坐标系到载体坐标系转换的旋转矩阵，R _x(ω)、

R _z(κ)分别为绕载体坐标系x轴、y轴、z轴旋转的旋转矩阵； among them,

Represents the rotation matrix transformed from the laser scanner coordinate system to the carrier coordinate system, R _x (ω),

R _z (κ) is the rotation matrix rotating around the x-axis, y-axis, and z-axis of the carrier coordinate system;

步骤2.2：从载体坐标系到站心坐标系转换；Step 2.2: Convert from the carrier coordinate system to the station center coordinate system;

惯导坐标系前进方向为y轴、向右为x轴、向上为z轴；The forward direction of the inertial navigation coordinate system is the y-axis, the right is the x-axis, and the upward is the z-axis;

惯导记录姿态角：侧滚角Roll、俯仰角Pitch、偏航角Heading；Inertial navigation records the attitude angle: roll angle Roll, pitch angle Pitch, yaw angle Heading;

侧滚角Roll：惯导x轴与水平方向之间的夹角，载体右侧向下为正；Roll angle Roll: the angle between the x-axis of the inertial navigation and the horizontal direction, the right side of the carrier downwards is positive;

俯仰角Pitch：惯导y轴与水平方向之间的夹角，载体向上为正；Pitch angle: the angle between the y-axis of the inertial navigation and the horizontal direction, the carrier is positive upward;

偏航角Heading：惯导前进方向，即xy平面与正北方向之间的夹角，顺时针为正；Heading: the heading direction of the inertial navigation, that is, the angle between the xy plane and the true north direction, clockwise is positive;

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

载体坐标系下坐标向站心坐标系转换； Let Roll, Pitch, Heading be r, p, y respectively; let the coordinates of the point in the station center coordinate system be

Conversion of coordinates under the carrier coordinate system to the station center coordinate system;

①先绕z轴旋转y；① Rotate y around the z axis first;

②再绕x轴旋转p；② Then rotate p around the x axis;

③最后绕y轴旋转r；③Finally rotate r around the y axis;

则有，Then there is,

其中，

代表从载体坐标系到站心坐标系转换的旋转矩阵，R _x(p)、R _y(r)、R _z(y)分别为绕站心坐标系x轴、y轴、z轴旋转的旋转矩阵； among them,

Represents the rotation matrix converted from the carrier coordinate system to the station center coordinate system, R _x (p), R _y (r), R _z (y) are the rotations around the station center coordinate system x axis, y axis, and z axis respectively matrix;

步骤2.3：从站心坐标系到地心地固坐标系转换；Step 2.3: Transform from the station center coordinate system to the geocentric ground-fixed coordinate system;

站心坐标系原点在WGS84下的经纬度分别为L和B；The latitude and longitude of the origin of the station center coordinate system under WGS84 are L and B respectively;

设扫描点在地心地固坐标系下的坐标为

站心坐标系下的坐标转换为WGS84坐标系下的坐标： Suppose the coordinates of the scanning point in the geocentric and ground-fixed coordinate system are

The coordinates in the station center coordinate system are converted to the coordinates in the WGS84 coordinate system:

①先绕x轴旋转

① Rotate around the x axis first

②再绕z轴旋转

②Rotate around the z axis again

③最后将站心坐标系原点平移到WGS84坐标系原点；③ Finally, the origin of the station center coordinate system is translated to the origin of the WGS84 coordinate system;

则有，Then there is,

其中，

代表从站心坐标系到地心地固坐标系转换的旋转矩阵； among them,

Represents the rotation matrix transformed from the station-centered coordinate system to the earth-centered ground-fixed coordinate system;

其中，N代表卯酉圈曲率半径，h代表大地高，a、b代表椭球长、短半径；Among them, N represents the radius of curvature of the unitary circle, h represents the height of the earth, and a and b represent the length and short radius of the ellipsoid;

得Get

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

It is the spatial rectangular coordinates of the origin of the station center coordinate system in the geocentric ground-fixed coordinate system;

最后得，Finally,

优选地，所述步骤3包括：Preferably, the step 3 includes:

首先介绍一下所需要的理论基础：First introduce the required theoretical basis:

平面拟合原理Plane fitting principle

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

且a′ ²+b′ ²+c′ ²＝1，d≥0，(x y z)是平面上的任意点； And a′ ² +b′ ² +c′ ² =1, d≥0, (x y z) is any point on the plane;

其中，(a′ b′ c′)代表垂直于平面且远离坐标原点方向的方向余弦；Among them, (a′ b′ c′) represents the cosine of the direction perpendicular to the plane and away from the coordinate origin;

同样的，它也表示垂直于平面远离坐标原点的单位法向量的分量；Similarly, it also represents the component of the unit normal vector perpendicular to the plane away from the origin of the coordinate;

d表示了平面和坐标原点间的垂直距离；d represents the vertical distance between the plane and the origin of the coordinate;

最小二乘算法基本原理The basic principle of least squares algorithm

最小二乘法的思想就是估算一个最佳数学模型，求得最佳模型参数使估算模型点的值和实际观察值的差值的平方和最小；The idea of the least squares method is to estimate an optimal mathematical model, and obtain the optimal model parameters to minimize the sum of squares of the difference between the estimated model point value and the actual observation value;

在二维数据中，假设给出模型函数g(x _i)，g(x _i)与真实值(x _i,y _i)的误差可用三种形式表述： In two-dimensional data, assuming that the model function g(x _i ) is given, the error between g(x _i ) and the true value (x _i , y _i ) can be expressed in three forms:

①误差最大值，即

其中n为点总个数； ①The maximum error, namely

Where n is the total number of points;

②误差绝对值之和，即

②The sum of absolute error values, namely

③误差平方和，即

③The error sum of squares, namely

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

最小； The principle of the least squares algorithm is to find an optimal parameter model. For a given data (x _i , y _i ) (i=0, 1...n), the sum of squared errors

The smallest

最小二乘平面拟合算法Least squares plane fitting algorithm

基于最小二乘法的平面拟合算法是在三维点云数据的基础上进行的拟合算法，以公式(6)为三维模型，确定a′,b′,c′,d四个参数使三维模型的估算值和测量值之间的差值平方和E ²最小；由于E ²是非负数，因此它存在极小值，E ²对每个参数的偏导数为零，利用偏导数为零构建的方程组便可求解出a′,b′,c′,d，并使估算值和测量值之间的差值平方和E ²最小； The plane fitting algorithm based on the least squares method is a fitting algorithm based on three-dimensional point cloud data. Formula (6) is used as the three-dimensional model to determine the four parameters a′, b′, c′, and d to make the three-dimensional model The sum of squares of the difference between the estimated value and the measured value of E ^{2 is the} smallest; since E ² is a non-negative number, it has a minimum value. The partial derivative of E ² for each parameter is zero. An equation constructed with a partial derivative of zero The group can solve for a′, b′, c′, d, and minimize the sum of squared differences E ² between the estimated value and the measured value;

假设点云数据的点集为{p _i|p _i∈oxyz,i＝1,2,3...,N}，其中x _i和y _i设为无误差变化的自变量，z _i是包含误差的因变量，则z和x、y的函数相关性可以假设为函数表达式： Suppose the point set of the point cloud data is {p _i |p _i ∈ oxyz,i=1, 2, 3...,N}, where x _i and y _{i are} set as independent variables without error change, and z _i is the inclusion The dependent variable of the error, the functional correlation between z and x and y can be assumed as a functional expression:

z＝f(x,y；p,q,r)＝p+qx+ry (7)z=f(x,y;p,q,r)=p+qx+ry (7)

且其中的p,q,r可以通过最小二乘算法得到；And the p, q, r can be obtained by the least square algorithm;

根据公式(6)、(7)，令d＝c′p，a′＝-qc′，b′＝-rc′且

则公式(7)可以描述为平面方程； According to formulas (6) and (7), let d=c′p, a′=-qc′, b′=-rc′ and

Then formula (7) can be described as a plane equation;

根据上述定义f函数的公式，假如自变量的值为x _i和y _i，被测量值z应该为p+qx _i+ry _i，然而因为误差的存在测量值z为z _i，则第i个点的测量误差为： According to the above formula defining the f function, if the values of the independent variables are x _i and y _i , the measured value z should be p+qx _i +ry _i , but because of the error, the measured value z is z _i , then the i-th The measurement error of the point is:

e _i＝f(x _i, _i；p,q,r)-z _i (8) e _i ＝f(x _i , _i ; p,q,r)-z _i (8)

根据最小二乘的原理，p,q,r应该满足公式(8)的平方和最小，即E ²最小： According to the principle of least squares, p, q, r should satisfy the minimum sum of squares of formula (8), that is, E ^{2 is} minimum:

若使得E ²最小，应满足： If E ^{2 is} minimized, it should satisfy:

即：which is:

有：Have:

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

求出； According to formula (11), p, q, r can be solved. Similarly, the four parameters of a′, b′, c′, d in formula (6) can also be solved by d=c′p, a′=-qc ′, b′=-rc′,

Find out

求解旋转矩阵Solve the rotation matrix

在上述算法的计算过程中，提取3对以上对应的平面特征量后，需要求取其各自的平面法向量；In the calculation process of the above algorithm, after extracting more than 3 pairs of corresponding plane feature quantities, it is required to take their respective plane normal vectors;

假设求解的一对平面为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矩阵，其夹角为： Assuming that a pair of planes to be solved is a ₁ x+b ₁ y+z=d ₁ and a ₂ x+b ₂ y+z=d ₂ , the normal vector pair can be expressed as P=(a ₁ ,b ₁ , -1), Q=(a ₂ ,b ₂ ,-1), the rotation matrix R of the normal vector satisfies P=RQ and is solved by the Rodrigue rotation method; the above P and Q matrix, the included angle is:

P,Q矩阵的叉乘为：The cross product of P and Q matrices is:

若已知单位向量

利用法向量将其旋转θ角度后，根据罗德里格旋转法求解其三维旋转矩阵为： If the unit vector is known

After the normal vector is used to rotate it by the angle θ, the three-dimensional rotation matrix is solved according to the Rodrigue rotation method as:

即可求出

的三维变换矩阵C ₁； Can be found

The three-dimensional transformation matrix C ₁ ;

求解平移量Solve for translation

利用三维变换矩阵C ₁可以对点云数据进行旋转，相应的平面特征量对相互平行，在两片点云数据中寻找一对以上的特征点，或利用对应面的的质心坐标的差值，计算点云数据的平移量即： Using the three-dimensional transformation matrix C ₁ to rotate the point cloud data, the corresponding plane feature pairs are parallel to each other, find more than one pair of feature points in two pieces of point cloud data, or use the difference between the centroid coordinates of the corresponding planes, Calculate the translation amount of point cloud data:

本发明具有的有益效果是：The beneficial effects of the present invention are:

本发明提供的一种空地结合的潮间带一体化测绘方法，集成了机载激光探测与测量***和拖曳测量***，可以获取潮间带区域水下地形和水上地形数据，解决了潮间带地形测量问题。两者能完美覆盖测量“盲区”，实现潮间带区域一体化无缝测量。The invention provides an integrated air-ground intertidal zone surveying and mapping method, which integrates an airborne laser detection and measurement system and a towed measurement system, can obtain underwater topography and water topography data in the intertidal zone, and solves the problem of intertidal zone Topographic survey problems. The two can perfectly cover the measurement "blind zone", and realize the integrated seamless measurement of the intertidal zone.

附图说明Description of the drawings

图1为本发明提供的潮间带载体车外业采集流程图。Figure 1 is a flow chart of the field collection of intertidal carrier vehicles provided by the present invention.

图2为本发明提供的机载激光探测与测量***外业采集流程图。Fig. 2 is a flow chart of field collection of the airborne laser detection and measurement system provided by the present invention.

图3为本发明提供的单波束数据格式转换流程图。Fig. 3 is a flow chart of single-beam data format conversion provided by the present invention.

图4为本发明提供的数据解算流程图。Figure 4 is a flow chart of data solution provided by the present invention.

图5为本发明提供的单波束测深仪数据、激光扫描仪数据、机载激光探测与测量***数据空间整合流程图。Fig. 5 is a flow chart of spatial integration of single-beam sounder data, laser scanner data, and airborne laser detection and measurement system data provided by the present invention.

具体实施方式detailed description

下面结合附图和具体实施例对本发明的具体实施方式做进一步说明：The specific implementation of the present invention will be further described below in conjunction with the drawings and specific embodiments:

结合图1至图5，一种空地结合的潮间带一体化测绘方法，包括测量装置，测量装置包括无人机和载体车，无人机上搭载激光探测与测量***，载体车上设置有激光扫描仪、单波束测深仪、惯性导航***和卫星定位***。With reference to Figures 1 to 5, an integrated air-to-ground intertidal zone surveying and mapping method includes a measuring device. The measuring device includes a UAV and a carrier vehicle. The UAV is equipped with a laser detection and measurement system, and the carrier vehicle is equipped with a laser. Scanner, single beam sounder, inertial navigation system and satellite positioning system.

步骤1：原始数据采集下载后，检查各测量仪器观测数据的完整性、正确性，进行观测数据文件的格式转换，进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，通过点云数据预处理，形成标准LAS数据文件。Step 1: After the original data is collected and downloaded, check the completeness and correctness of the observation data of each measuring instrument, convert the format of the observation data file, perform high-precision time synchronization of each measuring instrument, and integrate the time and space of the data of each measuring instrument. Point cloud data preprocessing to form standard LAS data files.

各测量仪器观测数据包括激光扫描仪数据、激光探测与测量***数据、单波束测深仪数据、惯性导航***数据和卫星定位***数据等。The observation data of each measuring instrument includes laser scanner data, laser detection and measurement system data, single beam echo sounder data, inertial navigation system data and satellite positioning system data.

具体的，步骤1包括以下子步骤：Specifically, step 1 includes the following sub-steps:

步骤1.1：如图1和图2所示，分别为载体车和机载激光探测与测量***数据采集流程图。图3为单波束数据格式转换流程图。Step 1.1: As shown in Figure 1 and Figure 2, the data acquisition flow chart of the carrier vehicle and the airborne laser detection and measurement system are respectively. Figure 3 is a flow chart of single-beam data format conversion.

利用硬件同步控制器进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，它作用是为***中的激光扫描仪、单波束测深仪等仪器提供相应的同步信号，并记录相应的时间，为后期的数据处理提供一个时间同步基准，硬件同步控制器的时间保持与卫星定位***的时间同步；硬件同步控制器控制各测量仪器的触发及监测其反馈信号，考虑到环境因素复杂，为保证***的高可靠性，信号的传输都采用差分方式；The hardware synchronization controller is used to synchronize the high-precision time of each measuring instrument, and the time and space of the data of each measuring instrument are integrated. Its function is to provide the corresponding synchronization signal for the laser scanner, single beam sounder and other instruments in the system, and record it The corresponding time provides a time synchronization benchmark for later data processing. The time of the hardware synchronization controller is kept synchronized with the time of the satellite positioning system; the hardware synchronization controller controls the triggering of each measuring instrument and monitors its feedback signals, taking environmental factors into account Complex, in order to ensure the high reliability of the system, the signal transmission adopts the differential mode;

步骤1.2：如图4所示，利用IE进行基准站和流动站数据结算，输出POS数据；Step 1.2: As shown in Figure 4, use IE to settle the data of the base station and the rover, and output the POS data;

步骤2：进行点云时空融合，完成从仪器坐标系到地心地固坐标系的转换。Step 2: Perform point cloud space-time fusion to complete the conversion from the instrument coordinate system to the geocentric ground-fixed coordinate system.

步骤2包括以下子步骤：Step 2 includes the following sub-steps:

κ;

l _x、l _y、l _z为坐标转换过程中点坐标沿载体坐标系x轴、y轴、z轴的平移量； l _x , l _y , l _z are the translational amount of point coordinates along the x-axis, y-axis and z-axis of the carrier coordinate system during the coordinate conversion process;

ω、

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

设激光扫描仪坐标系下点坐标为

在载体坐标系下的坐标为

The coordinates in the carrier coordinate system are

Then there is,

其中，

R _z(κ)分别为绕载体坐标系x轴、y轴、z轴旋转的旋转矩阵。 among them,

R _z (κ) is the rotation matrix rotating around the x-axis, y-axis, and z-axis of the carrier coordinate system.

惯导记录姿态角：侧滚角Roll、俯仰角Pitch、偏航角Heading。The inertial navigation records the attitude angle: roll angle Roll, pitch angle Pitch, yaw angle Heading.

①先绕z轴旋转y；① Rotate y around the z axis first;

②再绕x轴旋转p；② Then rotate p around the x axis;

③最后绕y轴旋转r；③Finally rotate r around the y axis;

则有，Then there is,

其中，

代表从载体坐标系到站心坐标系转换的旋转矩阵，R _x(p)、R _y(r)、R _z(y)分别为绕站心坐标系x轴、y轴、z轴旋转的旋转矩阵。 among them,

Represents the rotation matrix converted from the carrier coordinate system to the station center coordinate system, R _x (p), R _y (r), R _z (y) are the rotations around the station center coordinate system x axis, y axis, and z axis respectively matrix.

设扫描点在地心地固坐标系下的坐标为

①先绕x轴旋转

① Rotate around the x axis first

②再绕z轴旋转

②Rotate around the z axis again

则有，Then there is,

代表从站心坐标系到地心地固坐标系转换的旋转矩阵。

Represents the rotation matrix transformed from the station-centered coordinate system to the earth-centered ground-fixed coordinate system.

其中，N代表卯酉圈曲率半径，h代表大地高，a、b代表椭球长、短半径。Among them, N represents the radius of curvature of the unitary circle, h represents the height of the earth, and a and b represent the long and short radii of the ellipsoid.

得Get

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

最后得，Finally,

步骤3包括：Step 3 includes:

平面拟合原理Plane fitting principle

已知平面方程：Known plane equation:

a′x+b′y+c′z＝d (6)a′x+b′y+c′z=d (6)

其中(a′ b′ c′)代表垂直于平面且远离坐标原点方向的方向余弦；Where (a′ b′ c′) represents the cosine of the direction perpendicular to the plane and away from the origin of the coordinates;

d表示了平面和坐标原点间的垂直距离。d represents the vertical distance between the plane and the origin of the coordinate.

最小二乘算法基本原理The basic principle of least squares algorithm

最小二乘法的思想就是估算一个最佳数学模型，求得最佳模型参数使估算模型点的值和实际观察值的差值的平方和最小；The idea of the least square method is to estimate an optimal mathematical model, and obtain the optimal model parameters to minimize the sum of squares of the difference between the estimated model point value and the actual observation value;

①误差最大值，即

其中n为点总个数； ①The maximum error, namely

Where n is the total number of points;

②误差绝对值之和，即

②The sum of absolute error values, namely

③误差平方和，即

③The error sum of squares, namely

The smallest

最小二乘平面拟合算法Least squares plane fitting algorithm

z＝f(x,y；p,q,r)＝p+qx+ry (7)z=f(x,y;p,q,r)=p+qx+ry (7)

根据公式(6)、(7)，令d＝c′p，a′＝-qc′，b′＝-rc′且

Then formula (7) can be described as a plane equation;

e _i＝f(x _i,y _i；p,q,r)-z _i (8) e _i ＝f(x _i ,y _i ; p,q,r)-z _i (8)

若使得E ²最小，应满足： If E ^{2 is} minimized, it should satisfy:

即：which is:

有：Have:

求出。 According to formula (11), p, q, r can be solved. Similarly, the four parameters of a′, b′, c′, d in formula (6) can also be solved by d=c′p, a′=-qc ′, b′=-rc′,

Find out.

求解旋转矩阵Solve the rotation matrix

P,Q矩阵的叉乘为：The cross product of P and Q matrices is:

若已知单位向量

即可求出

的三维变换矩阵C ₁； Can be found

The three-dimensional transformation matrix C ₁ ;

求解平移量Solve for translation

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the essential scope of the present invention shall also belong to the present invention. The scope of protection of the invention.

Claims

一种空地结合的潮间带一体化测绘方法，其特征在于，包括测量装置，测量装置包括无人机和载体车，无人机上搭载激光探测与测量***，载体车上设置有激光扫描仪、单波束测深仪、惯性导航***和卫星定位***；An air-ground integrated intertidal zone surveying and mapping method is characterized by comprising a measuring device, the measuring device includes an unmanned aerial vehicle and a carrier vehicle, the unmanned aerial vehicle is equipped with a laser detection and measurement system, and the carrier vehicle is equipped with a laser scanner, Single beam depth sounder, inertial navigation system and satellite positioning system;

测绘方法包括以下步骤：The surveying method includes the following steps:

步骤1：原始数据采集下载后，检查各测量仪器观测数据的完整性、正确性，进行观测数据文件的格式转换，进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，通过点云数据预处理，形成标准LAS数据文件；Step 1: After the original data is collected and downloaded, check the completeness and correctness of the observation data of each measuring instrument, convert the format of the observation data file, perform high-precision time synchronization of each measuring instrument, and integrate the time and space of the data of each measuring instrument. Point cloud data preprocessing to form standard LAS data files;

步骤2：进行点云时空融合，完成从仪器坐标系到地心地固坐标系的转换；Step 2: Perform point cloud space-time fusion and complete the conversion from the instrument coordinate system to the geocentric ground-fixed coordinate system;

步骤3：通过获取单波束点云与激光扫描仪点云、激光扫描仪点云与机载LIDAR点云数据之间的点、线、面等几何特征来建立对应关系，求解旋转矩阵和平移量进行点云匹配。Step 3: Establish the correspondence relationship between the single-beam point cloud and the laser scanner point cloud, the laser scanner point cloud and the airborne LIDAR point cloud data, such as geometric features such as points, lines, and surfaces, and solve the rotation matrix and translation amount Perform point cloud matching.
根据权利要求1所述的一种空地结合的潮间带一体化测绘方法，其特征在于，所述步骤1包括以下子步骤：The method for integrated surveying and mapping of the intertidal zone combined with air and ground according to claim 1, wherein said step 1 comprises the following sub-steps:

步骤1.1：利用硬件同步控制器进行各测量仪器高精度时间同步，各测量仪器数据的时、空融合，硬件同步控制器的时间保持与卫星定位***的时间同步；硬件同步控制器控制各测量仪器的触发及监测其反馈信号，信号的传输都采用差分方式；Step 1.1: Use the hardware synchronization controller to synchronize the high-precision time of each measuring instrument, the time and space of the data of each measuring instrument are integrated, and the time of the hardware synchronization controller is synchronized with the time of the satellite positioning system; the hardware synchronization controller controls each measuring instrument Trigger and monitor its feedback signal, signal transmission adopts differential mode;

步骤1.2：利用IE进行基准站和流动站数据结算，输出POS数据；Step 1.2: Use IE to settle the data of base station and rover, and output POS data;

步骤1.3：利用VSursPROCESS软件进行点云数据预处理，形成标准LAS数据文件。Step 1.3: Use VsursPROCESS software to preprocess point cloud data to form a standard LAS data file.
根据权利要求1所述的一种空地结合的潮间带一体化测绘方法，其特征在于，所述步骤2包括以下子步骤：The method for integrated surveying and mapping of intertidal zone combined with air and ground according to claim 1, wherein said step 2 includes the following sub-steps:

步骤2.1：从激光扫描仪坐标系到载体坐标系转换；Step 2.1: Convert from the laser scanner coordinate system to the carrier coordinate system;

根据标定的公共点计算坐标系转换需要的的6参数：l _x、l _y、l _z、ω、
κ； Calculate the 6 parameters needed for coordinate system conversion according to the calibrated common point: l _x , l _y , l _z , ω,
κ;

l _x、l _y、l _z分别为坐标转换过程中点坐标沿载体坐标系x轴、y轴、z轴的平移量； l _x , l _y , and l _z are the translational amounts of the point coordinates along the x-axis, y-axis and z-axis of the carrier coordinate system during the coordinate conversion process;

ω、
κ分别为坐标转换过程中点坐标绕载体坐标系x轴、y轴、z轴的旋转角度； ω,
κ is the rotation angle of the point coordinate around the x-axis, y-axis and z-axis of the carrier coordinate system during the coordinate conversion process;

设激光扫描仪坐标系下点坐标为
在载体坐标系下的坐标为
则有， Let the coordinates of the point in the laser scanner coordinate system be
The coordinates in the carrier coordinate system are
Then there is,

其中，
代表从激光扫描仪坐标系到载体坐标系转换的旋转矩阵，R _x(ω)、
R _z(κ)分别为绕载体坐标系x轴、y轴、z轴旋转的旋转矩阵； among them,
Represents the rotation matrix transformed from the laser scanner coordinate system to the carrier coordinate system, R _x (ω),
R _z (κ) is the rotation matrix rotating around the x-axis, y-axis, and z-axis of the carrier coordinate system;

步骤2.2：从载体坐标系到站心坐标系转换；Step 2.2: Convert from the carrier coordinate system to the station center coordinate system;

惯导坐标系前进方向为y轴、向右为x轴、向上为z轴；The forward direction of the inertial navigation coordinate system is the y-axis, the right is the x-axis, and the upward is the z-axis;

惯导记录姿态角：侧滚角Roll、俯仰角Pitch、偏航角Heading；Inertial navigation records the attitude angle: roll angle Roll, pitch angle Pitch, yaw angle Heading;

侧滚角Roll：惯导x轴与水平方向之间的夹角，载体右侧向下为正；Roll angle Roll: the angle between the x-axis of the inertial navigation and the horizontal direction, the right side of the carrier downwards is positive;

俯仰角Pitch：惯导y轴与水平方向之间的夹角，载体向上为正；Pitch angle: the angle between the y-axis of the inertial navigation and the horizontal direction, the carrier is positive upward;

偏航角Heading：惯导前进方向，即xy平面与正北方向之间的夹角，顺时针为正；Heading: the heading direction of the inertial navigation, that is, the angle between the xy plane and the true north direction, clockwise is positive;

设Roll、Pitch、Heading分别为r、p、y；设点在站心坐标系下的坐标为
载体坐标系下坐标向站心坐标系转换； Let Roll, Pitch, Heading be r, p, y respectively; let the coordinates of the point in the station center coordinate system be
Conversion of coordinates under the carrier coordinate system to the station center coordinate system;

①先绕z轴旋转y；① Rotate y around the z axis first;

②再绕x轴旋转p；② Then rotate p around the x axis;

③最后绕y轴旋转r；③Finally rotate r around the y axis;

则有，Then there is,

其中，
代表从载体坐标系到站心坐标系转换的旋转矩阵，R _x(p)、R _y(r)、R _z(y)分别为绕站心坐标系x轴、y轴、z轴旋转的旋转矩阵； among them,
Represents the rotation matrix converted from the carrier coordinate system to the station center coordinate system, R _x (p), R _y (r), R _z (y) are the rotations around the station center coordinate system x axis, y axis, and z axis respectively matrix;

步骤2.3：从站心坐标系到地心地固坐标系转换；Step 2.3: Transform from the station center coordinate system to the geocentric ground-fixed coordinate system;

站心坐标系原点在WGS84下的经纬度分别为L和B；The latitude and longitude of the origin of the station center coordinate system under WGS84 are L and B respectively;

设扫描点在地心地固坐标系下的坐标为
站心坐标系下的坐标转换为WGS84坐标系下的坐标： Suppose the coordinates of the scanning point in the geocentric and ground-fixed coordinate system are
The coordinates in the station center coordinate system are converted to the coordinates in the WGS84 coordinate system:

①先绕x轴旋转
① Rotate around the x axis first

②再绕z轴旋转
②Rotate around the z axis again

③最后将站心坐标系原点平移到WGS84坐标系原点；③ Finally, the origin of the station center coordinate system is translated to the origin of the WGS84 coordinate system;

则有，Then there is,

其中，
代表从站心坐标系到地心地固坐标系转换的旋转矩阵； among them,
Represents the rotation matrix transformed from the station-centered coordinate system to the earth-centered ground-fixed coordinate system;

其中，N代表卯酉圈曲率半径，h代表大地高，a、b代表椭球长、短半径；Among them, N represents the radius of curvature of the unitary circle, h represents the height of the earth, and a and b represent the length and short radius of the ellipsoid;

得Get

为站心坐标系原点在地心地固坐标系下的空间直角坐标；
It is the spatial rectangular coordinates of the origin of the station center coordinate system in the geocentric ground-fixed coordinate system;

最后得，Finally,
根据权利要求1所述的一种空地结合的潮间带一体化测绘方法，其特征在于，所述步骤3包括：The method for integrated surveying and mapping of intertidal zone combined with air and ground according to claim 1, wherein the step 3 comprises:

首先介绍一下所需要的理论基础：First introduce the required theoretical basis:

平面拟合原理Plane fitting principle

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

且a′ ²+b′ ²+c′ ²＝1，d≥0，(x y z)是平面上的任意点； And a′ ² +b′ ² +c′ ² =1, d≥0, (x y z) is any point on the plane;

其中，(a′ b′ c′)代表垂直于平面且远离坐标原点方向的方向余弦；Among them, (a′ b′ c′) represents the cosine of the direction perpendicular to the plane and away from the coordinate origin;

同样的，它也表示垂直于平面远离坐标原点的单位法向量的分量；Similarly, it also represents the component of the unit normal vector perpendicular to the plane away from the origin of the coordinate;

d表示了平面和坐标原点间的垂直距离；d represents the vertical distance between the plane and the origin of the coordinate;

最小二乘算法基本原理The basic principle of least squares algorithm

最小二乘法的思想就是估算一个最佳数学模型，求得最佳模型参数使估算模型点的值和实际观察值的差值的平方和最小；The idea of the least squares method is to estimate an optimal mathematical model, and obtain the optimal model parameters to minimize the sum of squares of the difference between the estimated model point value and the actual observation value;

在二维数据中，假设给出模型函数g(x _i)，g(x _i)与真实值(x _i,y _i)的误差可用三种形式表述： In two-dimensional data, assuming that the model function g(x _i ) is given, the error between g(x _i ) and the true value (x _i , y _i ) can be expressed in three forms:

①误差最大值，即
其中n为点总个数； ①The maximum error, namely
Where n is the total number of points;

②误差绝对值之和，即
②The sum of absolute error values, namely

③误差平方和，即
③The error sum of squares, namely

最小二乘算法的原理便是寻找一个最佳参数模型，对于给定的数据(x _i,y _i)(i＝0,1…n)，误差平方和
最小； The principle of the least squares algorithm is to find an optimal parameter model. For a given data (x _i , y _i ) (i=0, 1...n), the sum of squared errors
The smallest

最小二乘平面拟合算法Least squares plane fitting algorithm

基于最小二乘法的平面拟合算法是在三维点云数据的基础上进行的拟合算法，以公式(6)为三维模型，确定a′,b′,c′,d四个参数使三维模型的估算值和测量值之间的差值平方和E ²最小；由于E ²是非负数，因此它存在极小值，E ²对每个参数的偏导数为零，利用偏导数为零构建的方程组便可求解出a′,b′,c′,d，并使估算值和测量值之间的差值平方和E ²最小； The plane fitting algorithm based on the least squares method is a fitting algorithm based on three-dimensional point cloud data. Formula (6) is used as the three-dimensional model to determine the four parameters a′, b′, c′, and d to make the three-dimensional model The sum of squares of the difference between the estimated value and the measured value of E ^{2 is the} smallest; since E ² is a non-negative number, it has a minimum value. The partial derivative of E ² for each parameter is zero. An equation constructed with a partial derivative of zero The group can solve for a′, b′, c′, d, and minimize the sum of squared differences E ² between the estimated value and the measured value;

假设点云数据的点集为{p _i|p _i∈oxyz,i＝1,2,3...,N}，其中x _i和y _i设为无误差变化的自变量，z _i是包含误差的因变量，则z和x、y的函数相关性可以假设为函数表达式： Suppose the point set of the point cloud data is {p _i |p _i ∈ oxyz,i=1, 2, 3...,N}, where x _i and y _{i are} set as independent variables without error change, and z _i is the inclusion The dependent variable of the error, the functional correlation between z and x and y can be assumed as a functional expression:

z＝f(x,y；p,q,r)＝p+qx+ry      (7)z=f(x,y;p,q,r)=p+qx+ry (7)

且其中的p,q,r可以通过最小二乘算法得到；And the p, q, r can be obtained by the least square algorithm;

根据公式(6)、(7)，令d＝c′p，a′＝-qc′，b′＝-rc′且
则公式(7) 可以描述为平面方程； According to formulas (6) and (7), let d=c′p, a′=-qc′, b′=-rc′ and
Then formula (7) can be described as a plane equation;

根据上述定义f函数的公式，假如自变量的值为x _i和y _i，被测量值z应该为p+qx _i+ry _i，然而因为误差的存在测量值z为z _i，则第i个点的测量误差为： According to the above formula defining the f function, if the values of the independent variables are x _i and y _i , the measured value z should be p+qx _i +ry _i , but because of the error, the measured value z is z _i , then the i-th The measurement error of the point is:

e _i＝f(x _i,y _i；p,q,r)-z _i        (8) e _i ＝f(x _i ,y _i ; p,q,r)-z _i (8)

根据最小二乘的原理，p,q,r应该满足公式(8)的平方和最小，即E ²最小： According to the principle of least squares, p, q, r should satisfy the minimum sum of squares of formula (8), that is, E ^{2 is} minimum:

若使得E ²最小，应满足： If E ^{2 is} minimized, it should satisfy:

即：which is:

有：Have:

其中i＝1…n   (12)
Where i=1...n (12)

根据公式(11)可以求解出p,q,r，同样的，公式(6)中a′,b′,c′,d四个参数也就可通过d＝c′p，a′＝-qc′，b′＝-rc′，
求出； According to formula (11), p, q, r can be solved. Similarly, the four parameters of a′, b′, c′, d in formula (6) can also be solved by d=c′p, a′=-qc ′, b′=-rc′,
Find out

求解旋转矩阵Solve the rotation matrix

在上述算法的计算过程中，提取3对以上对应的平面特征量后，需要求取其各自的平面法向量；In the calculation process of the above algorithm, after extracting more than 3 pairs of corresponding plane feature quantities, it is required to take their respective plane normal vectors;

假设求解的一对平面为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矩阵，其夹角为： Assuming that a pair of planes to be solved is a ₁ x+b ₁ y+z=d ₁ and a ₂ x+b ₂ y+z=d ₂ , the normal vector pair can be expressed as P=(a ₁ ,b ₁ , -1), Q=(a ₂ ,b ₂ ,-1), the rotation matrix R of the normal vector satisfies P=RQ and is solved by the Rodrigue rotation method; the above P and Q matrix, the included angle is:

P,Q矩阵的叉乘为：The cross product of P and Q matrices is:

若已知单位向量
利用法向量将其旋转θ角度后，根据罗德里格旋转法求解其三维旋转矩阵为： If the unit vector is known
After the normal vector is used to rotate it by the angle θ, the three-dimensional rotation matrix is solved according to the Rodrigue rotation method as:

即可求出
的三维变换矩阵C ₁； Can be found
The three-dimensional transformation matrix C ₁ ;

求解平移量Solve for translation

利用三维变换矩阵C ₁可以对点云数据进行旋转，相应的平面特征量对相互平行，在两片点云数据中寻找一对以上的特征点，或利用对应面的的质心坐标的差值，计算点云数据的平移量即： Using the three-dimensional transformation matrix C ₁ to rotate the point cloud data, the corresponding plane feature pairs are parallel to each other, find more than one pair of feature points in two pieces of point cloud data, or use the difference between the centroid coordinates of the corresponding planes, Calculate the translation amount of point cloud data: