Next Article in Journal
Phase Extraction from Single Interferogram Including Closed-Fringe Using Deep Learning
Previous Article in Journal
GaN/Ga2O3 Core/Shell Nanowires Growth: Towards High Response Gas Sensors
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Use of a Weighted ICP Algorithm to Precisely Determine USV Movement Parameters

Institute of Navigation and Maritime Hydrography, Polish Naval Academy, 81-103 Gdynia, Poland
*
Authors to whom correspondence should be addressed.
Appl. Sci. 2019, 9(17), 3530; https://doi.org/10.3390/app9173530
Submission received: 1 August 2019 / Revised: 16 August 2019 / Accepted: 23 August 2019 / Published: 28 August 2019

Abstract

:
The purpose of this article is to present a study aimed at developing a method for the precise determination of unmanned surface vehicle (USV) movement parameters (heading (HDG), speed over ground (SOG) and rate of turn (ROT)) through appropriate processing. The technique employs a modified weighted ICP (Iterative Closest Point) algorithm and a 2D points layer arranged in the horizon plane obtained from measurements. This is performed with the help of Light Detection and Ranging (LIDAR). A new method of weighting is presented. It is based on a mean error in a given direction and the results of modified weighted ICP tests carried out on the basis of field measurement data. The first part of the paper characterizes LIDAR measuring errors and indicates the possibilities for their use in matching point clouds. The second part of the article deals with a method for determining the SOG and course over ground (COG), based on a modified weighted ICP algorithm. The main part of the paper reviews a test method aimed at evaluating the accuracy of determining the SOG and COG by the scan-matching method using a modified weighted ICP algorithm. The final part presents an analysis comparing the obtained SOG and COG results with reference results of GNSS RTK measurements and the resulting generalised conclusions.

1. Introduction

Heading (HDG), speed over ground (SOG) and rate of turn (ROT) are basic movement parameters which characterize the way a ship’s hull moves in relation to the Earth’s surface. They are crucial in the decision-making process, particularly when navigating in congested areas, including harbors and rivers. Currently, these parameters are most frequently determined on a ship using a gyrocompass and a log, i.e., navigational instruments whose measuring accuracy, according to the International Maritime Organization (IMO) requirements, is not required to be very high: “Errors in the indicated speed, when the ship is operating free from shallow water effect and from the effects of wind, current and tide, should not exceed 2% of the speed of the ship, or 0.2 knots, whichever is greater” [1], “The follow-up error for different rates of turn should be: less than ±0.5° at rates up to 10°/s; and less than ±1.5° between a rate of 10°/s and 20°/s” [2], although it determines the precision and the ability to control the ship’s movement which, consequently, may affect the navigational safety level. It is therefore reasonable to carry out scientific research aimed at searching for new measuring and data processing methods in order to minimize measuring errors in a ship’s movement parameters.
For many years, maritime navigation has placed increasing emphasis on the development of methods to determine position and navigational parameters that could offer an alternative to the GNSS system [3]. In the age of the development of unmanned surface vehicles (USVs) whose main tasks increasingly involve independent (autonomous) maneuvering in restricted areas, it is necessary to develop new methods to obtain precise (accurate) navigational movement parameters (in particular the HDG, SOG, course over ground (COG), and ROT) [4,5]. Basing the USV decision-making process exclusively on the parameters of movement and positions obtained from satellite systems is unjustified, particularly in port basins in which infrastructure facilities and other ships may suppress and reflect the system signals [6,7,8]. Another reason why one should not rely exclusively on such solutions is the new techniques of satellite signal jamming and spoofing [9].
A good way to avoid the satellite system inconveniences may be the development of locally operating systems whose actuator is located on a carrier ship, making their operation more difficult to interfere with. A system with a local range of operation needs no special, external infrastructure to generate correct data. Recently, the clear development of devices such as Light Detection and Ranging (LIDAR) [10] has been noted, with increased capacities in terms of the range, accuracy of generated data and the purchase price. Certain devices have a range above 500 m and high accuracy in angle and distance measurements, which enables more widespread use of such devices in sea-going ship navigation. It should also be noted that the measuring capacity of such sensors is not limited by the radar’s dead zone at close ranges. On the other hand, such a phenomenon is typical of pulse radars which are commonly used on almost every vessel [11].
This article presents a method for the precise determination of the SOG and COG parameters through processing, using a modified iterative closest point (ICP) algorithm [12,13,14,15] of a horizontal 2D layer of cloud points obtained from LIDAR measurements. The main process is based on the simultaneous localization and mapping (SLAM) method. In simple terms, this is a method allowing a map of an unknown environment to be built, which, after processing, is used for navigating in a particular area. This technique was isolated as a separated area of interest in 1996. The SLAM evolution was described in [16]. SLAM is rarely used in determining the position of surface vessels. Only a few applications of the method in the context of its use on an open sea surface are known [17]. On the other hand, SLAM has become very popular in determining the position of autonomous underwater vehicles (AUV). A review of the methods applied in determining the positions of AUVs can be found in [18]. Since SLAM is already a relatively well-studied method, many approaches to solving a particular position-determining problem have emerged. One of the most common solutions to the problem is the so-called scan-matching. This set of techniques involves matching the subsequent scans of the surrounding environment which are obtained while exploring a particular area [19]. There are numerous algorithms variants used to match scans obtained by mapping devices. These include the ICP algorithm [20,21], ICL [22] and algorithms based on matching histograms of the distance function [23]. However, given the different methods of matching point clouds, the ICP algorithm proposed for the first time in [24] outperforms other algorithms in its simplicity and the possibility for its continuous development.
Considering the use of the SLAM technique to determine the parameters of the vessel’s movement, it is worth looking at the 2011 article by Callmer et al. [25]. In this case, the authors use the object approach (landmark) based on the SIFT algorithm and the extended Kalman filter (EKF), and images from the navigation radar. The results of the work are promising; however, in this case there is a significant drift of coordinates of the position and parameters of the vessels movement. Another publication describing the implementation of the SLAM task at sea [17]. The approach presented here assumes that the map of the environment will be explored on a simultaneous basis by more than one vehicle/ship. This assumption allows obtaining a much more accurate model of the map of the environment (global map) by combining many sub-maps obtained by vehicles. This mapping is used to position each of them.
The approach presented in the article is a scan matching without odometry information. A similar solution, but based on another scan matching algorithm, was presented in [26]. It involves adjusting concurrent scans one after another without using dead reckoning data in the alignment algorithm. This article uses the ICP algorithm with the error metric in the point-to-point version [27]. Other versions of the metrics are: point-to-line [28], point-to-plane [29], and plane-to-plane [30]. The point-to-point metric was chosen because of the intuitive way of attaching importance to individual correspondences (based solely on the coordinate error of the position of a pair of points), without the need to consider the uncertainty of the position of the line or plane (for cases of three-dimensional point clouds).
In this article, the emphasis is put on the methods of weighting on the basis of estimated errors of the position coordinates of a pair of corresponding points. The problem of weighting of pairs in a navigational approach is presented in [31], and for other applications in [32,33]. This publication assumes that pairs with a higher distance value between points will be given lower weights compared to those at a smaller distance from one another. Weighting based on color was also accepted. Another approach presented in [21] is the weighting depending on the uncertainty of the position of points in the scan and the uncertainty of the position of the straight line (in the point-to-line variant). It was assumed that more weight will be put on pairs that have less influence on the error metric. A similar approach is presented in [34]. The authors carry out tests on three variants of the error measure: point, line and plane. The approach proposed in this article is similar to this; however, it assumes a different, more universal distribution of coordinates error. They can be easily implemented not only in the case of a lidar sensor, but also in the case of navigation radar.

2. Evaluation of the Accuracy of Determining the Coordinates of a Position Using LIDAR

The determination of horizontal coordinates using LIDAR is carried out based on the angle α and the distance r measured between the sensor and the target, i.e., the light signal reflection point (Figure 1).
The mean error M and the parameters of mean error ellipse (i.e., the lengths of semi-axes a and b ) of the determination of 2D horizontal coordinates ( x , y ) using LIDAR can be calculated by applying the law of mean error propagation [35,36]. When the functions of a single result of the following measurement are known:
α = arc   tg x y ,
r = x 2 + y 2 ,
equations of the mean errors of positional lines can be written as:
σ l ( α ) = σ α / ( α x ) 2 + ( α y ) 2 ,
σ l ( r ) = σ r / ( r x ) 2 + ( r y ) 2 ,
which enable the determination of mean error of position:
M = 1 sin θ σ l ( α ) 2 + σ l ( r ) 2 = ( σ α · r ) 2 + σ r 2 ,
a = σ α · r ,
b = σ r ,
where: θ —the angle of intersection of positional lines, which for LIDAR measurements is always equal to 90°,
σα—mean error of the angle measurement α,
σ r —mean error of the distance measurement r ,
r —distance between the sensor and the target,
a —the length of the long semi-axis of mean error ellipse (see Figure 2),
b —the length of the short semi-axis of mean error ellipse (see Figure 2).
It follows from Dependence (5) that the value of mean error M and the length of the long semi-axis a of the mean error ellipse are significantly affected by the distance r between the sensor and the target. The calculations assumed that the measurements were performed using the Scanse Sweep LIDAR [37] for which the manufacturer provides σ r = 0.004258 · r (rounded value). Thereby one takes into account the fact that the error value of the angle measurement does not only depend on the beam’s divergence (footprint of the beam increases as a function of divergence) and the angle of its incidence on the surface of the object, but also on the accuracy of measuring the angle change by the optical head encoder that can be affected by environmental conditions such as temperature and humidity, vibrations and imperfections in the performance of mechanical elements that cooperate with each other. This can be characteristic for entry level LIDAR devices. The spatial angle resolution can be derived from sampling frequency that is set up to 1000 Hz (when optical sensor is spinning with frequency 1 Hz the Lidar is able to sample up to 1000 point for a single spin, 5 Hz—200 points, 10 Hz—100 points). The range spatial resolution is around 0.1 cm.
Based on the known lengths of semi-axes a and b of the mean error ellipse, and based on the known direction of α measurement, one can simply generate a covariance matrix c o v of the ( x ,   y ) vector coordinates:
c o v ( x ,   y ) = [   a 2 · cos 2 α a + b 2 · sin 2 α a ( a 2 b 2 ) · sin α a · cos α a ( a 2 b 2 ) · sin α a · cos α a a 2 · sin 2 α a + b 2 · cos 2 α a ] ,
where α a is the angle between the semi-axis a of the mean error ellipse and the axis OY (for LIDAR measurements, equal to α   + 90 ° ).
Then, using this, to determine the mean square error σ β 2 of the determination of the position coordinates, but in the specified direction β:
σ β 2 = S β · c o v ( x ,   y ) · S β T = sin β · ( sin β ( r 2 · σ α 2 · sin 2 α a + σ r 2 · cos α a ) cos α a · cos β · sin α a ( σ r 2 r 2 · σ α 2 ) ) + + sin β · ( sin β · ( r 2 · σ α 2 · cos 2 α a + σ r 2 · sin 2 α a ) cos α a · sin α a · sin β · ( σ r 2 r 2 · σ α 2 ) ) ,
where S β = [ cos β sin β ] .

2.1. The Determination of USV Movement Parameters Using the Scan-Matching Method

The taxonomy of the ICP algorithm is carried out in the following successive stages of data processing [21]:
  • Selection (selection of points which are suitable for the alignment);
  • Matching (matching points from model cloud to reference points cloud);
  • Weighting (giving weight to corresponding points);
  • Rejecting (deleting some points based on the robust criterion function);
  • Assigning an error metric (choosing point-to-point, point-to-line, point-to-plane or plane-to-plane error metric);
  • Minimizing the metric between selected points.
Selection. Let us assume that a USV determines the change in its position through the translation vector T ( k + 1 ) and the rotation matrix R ( k + 1 ) determined based on two maps of the surrounding environment M ( k ) and M ( k + 1 ) , constructed at the moment ( k ) and the moment ( k + 1 ) following it (e.g., 2D horizontal layer of cloud points for LIDAR measurements). Let ( k ) and ( k + 1 )   represent two sets { ( r 1 ( k ) , α 1 ( k ) ) ,   , ( r n ( k ) , α n ( k ) ) } and { ( r 1 ( k + 1 ) , α 1 ( k + 1 ) ) ,   , ( r n ( k + 1 ) , α n ( k + 1 ) ) } of measurements of the distance r and direction α, carried out in relation to n field obstacles (see Figure 1). On their basis, set of points M ( k ) = { p 1 ( k ) ,   , p n ( k ) } and M ( k + 1 ) = { p 1 ( k + 1 ) ,   , p n ( k + 1 ) } will be determined, where the coordinates of each point in the formula are calculated via the following dependence:
p i = [ x i y i ] , x i = r i · cos ( α i ) , y i = r i · sin ( α i ) ,
where:
i = 1 , , n
They will be used to determine the translation vector T ( k + 1 ) and the rotation matrix R ( k + 1 ) using the singular value decomposition (SVD) method [38], which will be carried out in combination with a change in the position (in subsequent iteration steps) of the set of points from M ( k ) to M ( k + 1 ) using the modified weighted ICP method. In reality, there is very little probability that both sets of points M ( k ) and M ( k + 1 ) will contain the same number of points. Equalising the number of points is required for the family of ICP algorithms with point-to-point error metric and therefore, in the calculations they were carried out by reducing the number of points. From a more numerous set, the points were randomly removed so that their number would be equal to the number of the points from a smaller set.
Matching. Literature shows that there are many ways to match points in pairs. One of them is the k-d tree method presented in [39,40]. Also considered are methods based on the heuristic search of pairs to reduce the complexity of calculations [41]. In this article, correspondence between points is created based on the method presented in [42]. It is based on the Delaunay triangulation. The method generates a two-dimensional grid of triangles based on the reference set M ( k ) . Next, using the nearest neighbor method, the nearest apex of triangles to the points from the set M ( k + 1 ) are found (points from M ( k ) ) in such a way that for each point p i ( k ) the nearest point p i d ( i ) ( k + 1 )   is   assigned (see Figure 3).
Weighting of the pairs. The modification of the ICP method will involve the application of the weighting factor w 1 , , n ( k , k + 1 ) for calculations for each pair of points (nearest neighbors):
w i ( k , k + 1 ) = 1 G i ( k ,   k + 1 ) ,
where G i ( k ,   k + 1 ) is the sum of mean errors of the determination of the coordinates of points p i ( k ) and p i d ( i ) ( k + 1 )   in the specified direction β (see Figure 2):
G i ( k , k + 1 ) = σ β i 2 ( k ) + σ β i d ( i ) 2 ( k + 1 ) ,
where σ β i ( k ) is the mean error of the determination of the coordinates of the point p i ( k ) in the direction of β (from p i ( k ) to p i d ( i ) ( k + 1 ) ) , σ β i d ( i ) ( k + 1 ) is the mean error of the determination of the coordinates of the point p i ( k + 1 ) in the direction of β (from p i d ( i ) ( k + 1 ) to p i ( k ) ) , or G i ( k ,   k + 1 ) is the mean error of the determination of the coordinates of the points p i ( k ) and p i d ( i ) ( k + 1 ) (see Figure 1):
G i ( k , k + 1 ) = M i 2 ( k ) + M i d ( i ) 2 ( k + 1 ) ,
where M i ( k ) is the mean error of the determination of the coordinates of the point p i ( k ) , M i d ( i ) ( k + 1 ) is the mean error of the determination of the coordinates of the point p i d ( i ) ( k + 1 ) .
Weighting step gives specific calculated value for each pair. Then in further processing heavier pairs involve more the translation and rotation parameters (Figure 4).
Rejection. Due to the relatively high noise level of many scans, particularly in the sectors illustrating long distance measurements, the Huber function [43] which attenuated the outlying measurement results (affected by gross errors) was also applied in the calculations. This function was arbitrarily selected from among many functions known and commonly used for this purpose [21,31,44,45,46,47]. The robust function rejects or significantly reduces weights of the specific pairs due to robustness criterion (see Figure 5). The robust function is based on k H u factor presented in Algorithm 1. Rejecting or weight minimization is performed based on norm of two paired points. In this article the fusion of the robust function and proposed error weighting is used to derivate the final weight.
Assigning an error metric. The final function minimizing the value of the matching error E , used to accurately determine the translation vector T ( k + 1 ) and the rotation matrix R ( k + 1 ) , will take the following form:
E ( R ,   T ) =   i = 1 n w i ( k , k + 1 ) · R p i ( k ) + T p i d ( i ) ( k + 1 ) 2 .
Provided that points p i ( k ) and p i d ( i ) ( k + 1 ) are located the closest to one another, i.e., are the nearest neighbors. This method is based on the point-to-point measure of error.
Minimizing the metric between selected points. The method of minimizing the error between the corresponding points was carried out on the basis of the SVD (Singular Value Decomposition). The method works equally for the weighted and classical ICP algorithm. In this article a method identical to [47] was used. It is required to determine two weighted average coordinates from the cloud m ( k ) and m ( k + 1 ) , dependent on point clouds ( k ) , ( k + 1 ) and assigned weights in the previous stage:
m ( k ) = i = 1 n w i · p i ( k ) i = 1 n w i ,
m ( k + 1 ) = i = 1 n w i d ( i ) · p i d ( i ) ( k + 1 ) i = 1 n w i .
Subsequently, we can move on to determine the weighted covariance matrix, which will be used to calculate the innovation of rotation matrix and translation (in a given iteration of the algorithm)
C = i = 1 n ( w i · p i ( k ) · p id ( i ) T ( k + 1 ) ) i = 1 n w i m ( k ) · m T ( k + 1 ) .
Using the SVD decomposition method on the C matrix
C = U Σ V T .
where matrices U, Σ, V are characteristic matrices for the SVD method. By means of decomposition carried out on designated matrices, we obtain the innovations in given iteration— i t —of the R i t rotation matrix and the T i t translation:
R i t = V U T ,
T i t = m ( k ) R i t · m ( k + 1 ) .
The obtained values in the R i t e r and T i t e r matrices are partial values developed in a given iteration of the algorithm. Based on this, the set m(k + 1) is translated at last iteration (see Figure 6) and as a result the translation and rotation is obtained:
R ( k + 1 ) = R 1 · · R i t n ,
T ( k + 1 ) = ( R i t + 1 · T i t + T i t + 1 ) +     + ( R itn · T itn 1 + T itn ) .
where R i t = 1 and T i t = 1 are initialized values.
The translation vector T ( k + 1 ) = [ Δ x Δ y ] thus obtained will be used to determine
SOG ( k + 1 ) = Δ x 2 + Δ y 2 Δ t ,
where Δ t is the distance in time between the moment ( k ) of performing the measurements of ( k ) , and the moment ( k + 1 ) of performing subsequent measurements of ( k + 1 ) .
On the other hand, the final form of the rotation matrix R ( k + 1 ) = [ cos Δ θ sin Δ θ sin Δ θ cos Δ θ ] will be used to determine:
ROT ( k + 1 ) = arc   sin ( Δ θ ) Δ t ,
and
HDG ( k + 1 ) = HDG 0 + arc   sin ( Δ θ ) ,
where HDG 0 is the initial HDG value determined using other methods (devices) prior to the commencement of measurement.

2.2. Weighted Matching Point Clouds Using the Mean Direction Error

Figure 7 shows that the mean error value changes significantly as the function of angle β , and the gradient of these changes will increase while the eccentricity of the ellipse is approaching 1 (the ellipse will be becoming flat, which occurs where measurements are carried out using LIDAR). Owing to these particular properties, this error can usefully determine the value of translation and rotation established while minimalizing the error metric between the points from subsequent measurements (e.g., point clouds from LIDAR measurements). Figure 8 presents, in the graphical form, how the positions matched in pairs of points based on the corrections (e.g., those added to the components of the rotation matrix and translation vector) change; the value of the rotation and translation was made depending on the values of the sums of mean errors (Equation (9)) calculated along lines connecting the matched points from ( k ) and ( k + 1 ) cloud. In order to facilitate the interpretation of the phenomenon, it was assumed that the mean error ellipses were similar in size.
As can be seen in Figure 7, greater weight is introduced for the points whose mean errors in the specified direction (point-to-point direction) is lower–pairs marked as 2 and 3 in Figure 7. The heavier weight influences the rotation and translation more and that is why pair 2 and 3 are better aligned comparing with 1 and 4. The use of such a weighting coefficient allows taking into account the spatial distribution of errors according to the selected direction. For comparison, Figure 8 presents, in a graphical form, how the positions matched in pairs of points based on the corrections whose values depended on the values of the sums of mean errors (Equation (5)) calculated for the matched points. As can be seen each pair (1, 2, 3, 4) influences the rotation and translation in similar manner.
HDG, SOG, and ROT were determined by processing LIDAR scans archived in data packets in accordance with the methodology presented in Section 2.1. The pseudo-code of the program used for the calculations using the authors original weighting factor w ( k , k + 1 ) and Huber’s robust function coefficient k H u = c · Median p 1 , , n ( k ) p i d ( 1 , n ) ( k + 1 ) finally took the following form:
Algorithm 1 The weighted ICP
Applsci 09 03530 i001

3. The Study and the Analysis of the Obtained Results

The main aim of the study was to evaluate the accuracy of the determination of the HDG, SOG and ROT using a modified ICP algorithm, carried out as a result of their comparison with high-accuracy reference measurements. The study involved:
  • measurements carried out using LIDAR and a GNSS RTK receiver, and the synchronous recording of their results,
  • the determination of the HDG, SOG and ROT, based on LIDAR measurement results, by the scan-matching method using both classic and modified versions of the ICP algorithm,
  • an analysis comparing the HDG, SOG, and ROT with reference results of GNSS RTK measurements, based on the accuracy criterion.
The yacht port basin in Gdynia (Poland) was selected for the measurement area. Its boundaries include high, concrete wharves and a breakwater. Most of the time, good hydro-meteorological conditions prevail in the area (Figure 9). To carry out the measurements, a small USV (with a length of 1.62 m, width of 0.40 m, and a draught of 0.11 m) equipped with Scanse Sweep Lidar, a GNSS RTK receiver Leica Viva CS 15 [48], and an on-board computer connected to them with a RS-422/232C conversion cable [49], were used to record the measurement results (Figure 10). The measurements were carried out on a USV sailing with a speed of approximately 1.2 m s   ( 2.4   kts ) over a trajectory of approximately 220   m , as presented in Figure 10. All data was collected via a mobile computer wirelessly connected to AUSV (autonomous unmanned surface vehicle, see Figure 11.). While the USV was sailing, the on-board computer synchronously recorded (when generating a full scan using LIDAR) the measurement results making up the so-called data sets containing: a LIDAR scan (in the form of a measurement set = { ( r 1 , α 1 ) ,   , ( r n , α n ) } ), NMEA messages with the position coordinates of SOG and COG (course over ground) from the GNSS RTK receiver [50]. 865 data sets were thus collected by recording each subsequent sets every 0.2   s on average, i.e., after a change in the USV position by approx. 0.25   m .

Analysis of the Accuracy of Determining the COG and SOG Using a Modified ICP Algorithm

Table 1 presents statistical parameters (i.e., the mean value, standard deviation and the maximum value) which characterize the computational accuracy of each parameter. The following acronyms were used to describe ICP’s variants: K_ICP for algorithm without modification, H_ICP for robust version based on Huber’s function, HM_ICP for robust version based on Huber’s function with mean error weighting and HD_ICP for robust version based on Huber’s function with directional error weighting. The lower index ICP was used for values computed by ICP algorithm’s variants and lower index R was introduced to point the values surveyed with GNSS RTK receiver.
A comparison of the parameter values (obtained based on 865 collected data packets) listed in Table 1 shows that the HD_ICP algorithm was the best in terms of the accuracy of determining the USVs movement parameters; it is followed by HM_ICP, H_ICP and K_ICP. Figure 12 shows a graph of SOG_ICP-SOG_R differences.
Figure 12 shows that the accuracy of determining the SOG using the K_ICP method significantly decreases for scans in the intervals of 220 ,   365 and 720 ,   825 . These scans differ significantly from each other in rotation angles and include a greater number of erroneous measurement results than the others. This is undoubtedly due to the turns made by the USV at that time. The accuracy of determining the SOG using the other methods is, at the same time, considerably higher. These methods probably achieve that through the attenuation of measurements with gross errors using Huber’s function. The HD_ICP method proved to be the best in determining the SOG in a significant part of the graph course. This was also confirmed by the lowest values of the mean error, equal to 0.007 m s ( 0.014   kts ) and of the standard deviation, equal to 0.026   m s   ( 0.05   kts ) presented in Table 1.
The histogram representing the frequency of occurrence of SOG ICP SOG R differences, presented in Figure 13, shows even more clearly that in terms of the accuracy of determining the SOG, the HD_ICP method is better than the others. Most of the values of SOG ICP SOG R differences, determined using the HD_ICP, HM_ICP, and H_ICP methods, fall within the interval of −0.05~0.05.
Figure 14 shows that the accuracy of determination of the Δ COG using the K_ICP method significantly decreases, as in the case of the SOG, for scans in the intervals of 220 ,   365 and 720 ,   825 . The HD_ICP method also proved to be the best in determining the Δ COG in this case. This was also confirmed by the lowest values of the mean error, equal to 0.03 ° , and of the standard deviation, equal to 0.11 ° , presented in Table 1. Figure 15 shows a histogram representing the frequency of occurrence of Δ COG ICP Δ COG R increment differences.
The histogram representing the frequency of the occurrence of differences Δ COG ICP Δ COG R shown in Figure 15 enables the conclusion that the HD_ICP method is also the best in terms of the accuracy of determining the Δ COG . Its advantage over the others can be clearly seen in the interval of differences Δ COG ICP Δ COG R 0.05 ˚ ,   0.05 ˚ . Given that the Δ ROT results are identical with the Δ COG results, their presentation and analysis have been omitted. Differences in the algorithms’ accuracy are shown at Figure 16 (presenting results of H_ICP and HD_ICP). Differences between HD_ICP and HM_ICP versions are barely visible, so they were omitted.

4. Conclusions

As the study attempted to prove, an ICP algorithm that matches point clouds from a laser scanner can effectively generate USV movement parameters. The proposed weighting factor based on the mean error in the specified direction of the determination of the coordinates of the positions of scan points enables more realistic (corresponding to the actual movement) matching of subsequent maps of the surrounding environment. This helps to determine the USV movement parameters with a higher accuracy. This was confirmed by the study results—inter alia the values of statistical indices in the form of mean SOG ICP SOG R = 0.007   m s   ( 0.014   kts ) and standard deviation SOG ICP SOG R = 0.0025   m s   ( 0.05   kts ) and Δ COG ICP Δ COG R = 0.03 ° and standard deviation Δ COG ICP Δ COG R = 0.11 ° . However, it should be noted that the authors’ original weighting factor should be used in combination with a robust criterion function to reduce measurements with gross errors (e.g., based on Huber’s function), as only such a solution can be fully applicable. The indicated mean error values lead us to a generalized statement that the developed method allows measuring SOG and COG with the accuracy required by the IMO. It should be noted, however, that the test was carried out in confined waters (harbors, channels, berths, marinas), beyond which the method would not be effective due to the low LIDAR range. All data used in publication is available at [51].

Author Contributions

Ł.M., methodology, investigation, writing-original draft preparation, writing-review and editing, K.N., supervision, formal analysis, validation and funding acquisition.

Funding

The article has been funded by a statutory research on marine optical navigation systems conducted in Polish Naval Academy.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. IMO. Resolution A.824(19) Adoption of Amendments to Performance Standards for Devices to Measure and Indicate Speed and Distance; IMO: London, UK, 2000. [Google Scholar]
  2. IMO. Resolution MSC.116(73), Performance Standards for Marine Transmitting Heading Devices (thds); IMO: London, UK, 2000. [Google Scholar]
  3. Naus, K.; Waz, M. Precision in Determining Ship Position using the Method of Comparing an Omnidirectional Map to a Visual Shoreline Image. J. Navig. 2016, 69, 391–413. [Google Scholar] [CrossRef]
  4. Polvara, R.; Sharma, S.; Wan, J.; Manning, A.; Sutton, R. Obstacle Avoidance Approaches for Autonomous Navigation of Unmanned Surface Vehicles. J. Navig. 2018, 71, 241–256. [Google Scholar] [CrossRef]
  5. Allen, C.H. The Seabots are Coming Here: Should they be Treated as ‘Vessels’? J. Navig. 2012, 65, 749–752. [Google Scholar] [CrossRef]
  6. Nowak, A. The Proposal to “Snapshot” Raim Method for Gnss Vessel Receivers Working in Poor Space Segment Geometry. Pol. Marit. Res. 2015, 22, 3–8. [Google Scholar] [CrossRef]
  7. Xie, P.; Petovello, M.G. Measuring GNSS Multipath Distributions in Urban Canyon Environments. IEEE Trans. Instrum. Meas. 2015, 64, 366–377. [Google Scholar]
  8. Wang, Y.; Chen, X.; Liu, P. Statistical Multipath Model Based on Experimental GNSS Data in Static Urban Canyon Environment. Sensors 2018, 18, 1149. [Google Scholar] [CrossRef]
  9. Shafiee, E.; Mosavi, M.R.; Moazedi, M. Detection of Spoofing Attack using Machine Learning based on Multi-Layer Neural Network in Single-Frequency GPS Receivers. J. Navig. 2018, 71, 169–188. [Google Scholar] [CrossRef]
  10. Guo, R.; Sun, F.; Yuan, J. ICP based on Polar Point Matching with application to Graph-SLAM. In Proceedings of the 2009 International Conference on Mechatronics and Automation, Changchun, China, 9–12 August 2009; pp. 1122–1127. [Google Scholar]
  11. Navarro, W.; Vélez, J.C.; Orfila, A. Estimation of sea state parameters using X-band marine radar technology in coastal areas. In Proceedings of the Sixth International Conference on Remote Sensing and Geoinformation of the Environment (RSCy2018), Paphos, Cyprus, 26–29 March 2018. [Google Scholar]
  12. Mendes, E.; Koch, P.; Lacroix, S. ICP-based pose-graph SLAM. In Proceedings of the 2016 IEEE International Symposium on Safety, Security, and Rescue Robotics (SSRR), Lausanne, Switzerland, 23–27 October 2016; pp. 195–200. [Google Scholar]
  13. Yang, J.; Li, H.; Jia, Y. Go-ICP: Solving 3D Registration Efficiently and Globally Optimally. In Proceedings of the 2013 IEEE International Conference on Computer Vision, Sydney, Australia, 1–8 December 2013; pp. 1457–1464. [Google Scholar]
  14. Zhang, Z. Iterative point matching for registration of free-form curves and surfaces. Int. J. Comput. Vis. 1994, 13, 119–152. [Google Scholar] [CrossRef]
  15. Men, H.; Gebre, B.; Pochiraju, K. Color point cloud registration with 4D ICP algorithm. In Proceedings of the 2011 IEEE International Conference on Robotics and Automation, Shanghai, China, 9–13 May 2011; pp. 1511–1516. [Google Scholar]
  16. Durrant-Whyte, H.; Bailey, T. Simultaneous localization and mapping: Part I. IEEE Robot. Autom. Mag. 2006, 13, 99–110. [Google Scholar] [CrossRef]
  17. Moratuwage, M.D.P.; Wijesoma, W.S.; Kalyan, B.; Dong, J.F.; Namal Senarathne, P.G.C.; Hover, F.S.; Patrikalakis, N.M. Collaborative multi-vehicle localization and mapping in marine environments. In Proceedings of the OCEANS’10 IEEE SYDNEY, Sydney, Australia, 24–27 May 2010; pp. 1–6. [Google Scholar]
  18. Ribas, D.; Ridao, P.; Neira, J. Underwater SLAM for Structured Environments Using an Imaging Sonar; Springer Tracts in Advanced Robotics Springer Berlin Heidelberg: Berlin/Heidelberg, Germany, 2010; ISBN 978-3-642-14039-6. [Google Scholar]
  19. Nieto, J.; Bailey, T.; Nebot, E. Recursive scan-matching SLAM. Robot. Auton. Syst. 2007, 55, 39–49. [Google Scholar] [CrossRef]
  20. Pomerleau, F.; Colas, F.; Siegwart, R.; Magnenat, S. Comparing ICP variants on real-world data sets: Open-source library and experimental protocol. Auton. Robot. 2013, 34, 133–148. [Google Scholar] [CrossRef]
  21. Rusinkiewicz, S.; Levoy, M. Efficient variants of the ICP algorithm. In Proceedings of the Third International Conference on 3-D Digital Imaging and Modeling, Quebec City, QC, Canada, 28 May–1 June 2001; pp. 145–152. [Google Scholar]
  22. Olson, E. Robust and Efficient Robotic Mapping. Ph.D. Thesis, Department of Electrical Engineering and Computer Science, MIT, Cambridge, MA, USA, 2008. [Google Scholar]
  23. Konecny, J.; Prauzek, M.; Kromer, P.; Musilek, P. Novel Point-to-Point Scan Matching Algorithm Based on Cross-Correlation. Mob. Inf. Syst. 2016, 2016, 1–11. [Google Scholar] [CrossRef] [Green Version]
  24. Besl, P.J.; McKay, N.D. A method for registration of 3-D shapes. IEEE Trans. Pattern Anal. Mach. Intell. 1992, 14, 239–256. [Google Scholar] [CrossRef]
  25. Callmer, J.; Törnqvist, D.; Gustafsson, F.; Svensson, H.; Carlbom, P. Radar SLAM using visual features. EURASIP J. Adv. Signal Process. 2011, 2011, 71–82. [Google Scholar] [CrossRef]
  26. Amigoni, F.; Gasparini, S.; Gimi, M. Scan matching without odometry information. In Proceedings of the First International Conference on Informatics in Control, Automation and Robotics, SciTePress—Science and Technology Publications, Setúbal, Portugal, 25–28 August 2004; pp. 349–352. [Google Scholar]
  27. He, Y.; Liang, B.; Yang, J.; Li, S.; He, J. An Iterative Closest Points Algorithm for Registration of 3D Laser Scanner Point Clouds with Geometric Features. Sensors 2017, 17, 1862. [Google Scholar] [CrossRef] [PubMed]
  28. Censi, A. An ICP variant using a point-to-line metric. In Proceedings of the 2008 IEEE International Conference on Robotics and Automation, Pasadena, CA, USA, 19–23 May 2008; pp. 19–25. [Google Scholar]
  29. Segal, A.; Haehnel, D.; Thrun, S. Generalized-ICP. In Proceedings of the Robotics: Science and Systems V; Robotics: Science and Systems Foundation, Boston, MA, USA, 28 June–01 July 2009. [Google Scholar]
  30. Wang, Y.; Xiong, R.; Li, Q. Em-based point to plane icp for 3d simultaneous localization and mapping. Int. J. Robot. Autom. 2013, 28. [Google Scholar] [CrossRef]
  31. Godin, G.; Rioux, M.; Baribeau, R. Three-Dimensional Registration Using Range and Intensity Information; The International Society for Optical Engineering: Boston, MA, USA, 1994. [Google Scholar]
  32. Marinov, A.; Zlateva, N.; Dimov, D.; Marinov, D. Weighted ICP Algorithm for Alignment of Stars from Scanned Astronomical Photographic Plates. Serdica J. Comput. Sofia Bulg. 2012, 6, 101–110. [Google Scholar]
  33. Amamra, A.; Aouf, N.; Stuart, D.; Richardson, M. A recursive robust filtering approach for 3D registration. SignalImage Video Process. 2016, 10, 835–842. [Google Scholar] [CrossRef]
  34. Wang, R.; Geng, Z. WA-ICP algorithm for tackling ambiguous correspondence. In Proceedings of the 2015 3rd IAPR Asian Conference on Pattern Recognition (ACPR), Kuala Lumpur, Malaysia, 3–6 November 2015; pp. 76–80. [Google Scholar]
  35. Naus, K.; Wąż, M. Ripping Radar Image with the Sea Map Using the Alignment Method 2009–2019; Logistyka, Institute of Logistics and Warehousing: Poznan, Poland, 2010. [Google Scholar]
  36. Naus, K.; Wąż, M. Accuracy in fixing ship’s positions by camera survey of bearings. Geod. Cartogr. 2011, 60, 61–73. [Google Scholar] [CrossRef]
  37. Scansce. User’s Manual and Technical Specifications. Available online: www.robotshop.com/ media/files/pdf2/sweep_user_manual.pdf (accessed on 1 February 2019).
  38. Marden, S.; Guivant, J. Improving the Performance of ICP for Real-Time Applications using an Approximate Nearest Neighbour Search. In Proceedings of the Australasian Conference on Robotics and Automation, Wellington, New Zealand, 3–5 December 2012; pp. 3–5. [Google Scholar]
  39. Bentley, J.L. Multidimensional binary search trees used for associative searching. Commun. ACM 1975, 18, 509–517. [Google Scholar] [CrossRef]
  40. Greenspan, M.; Yurick, M. Approximate K-D tree search for efficient ICP. In Proceedings of the Fourth International Conference on 3-D Digital Imaging and Modeling, Banff, AB, Canada, 6–10 October 2003; pp. 442–448. [Google Scholar]
  41. Jost, T.; Hugli, H. A multi-resolution ICP with heuristic closest point search for fast and robust 3D registration of range images. In Proceedings of the Fourth International Conference on 3-D Digital Imaging and Modeling, Banff, AB, Canada, 6–13 October 2003; pp. 427–433. [Google Scholar]
  42. Eggert, D.; Dalyot, S. Octree-based simd strategy for icp registration and alignment of 3d point clouds. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci. 2012, I–3, 105–110. [Google Scholar] [CrossRef]
  43. Huber, P.J. Robust Estimation of a Location Parameter. In Breakthroughs in Statistics: Methodology and Distribution; Kotz, S., Johnson, N.L., Eds.; Springer New York: New York, NY, USA, 1992; pp. 492–518. ISBN 978-1-4612-4380-9. [Google Scholar]
  44. Zhao, Q.; Han, X.H.; Chen, Y.W. A robust registration method using Huber ICP and low rank and sparse decomposition. In Proceedings of the 2015 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA), Hong Kong, China, 16–19 December 2015; pp. 744–752. [Google Scholar]
  45. Zhllin, L. Robust surface matching for automated detection of local deformations using least-median-of-squares estimator. Photogramm. Eng. Remote Sens. 2001, 67, 1283–1292. [Google Scholar]
  46. Masuda, T.; Sakaue, K.; Yokoya, N. Registration and integration of multiple range images for 3-D model construction. In Proceedings of the 13th International Conference on Pattern Recognition, Vienna, Austria, 25–29 August 1996; pp. 879–883. [Google Scholar]
  47. Bergström, P.; Edlund, O. Robust registration of point sets using iteratively reweighted least squares. Comput. Optim. Appl. 2014, 58, 543–561. [Google Scholar] [CrossRef]
  48. Leica Geosystems. Leica GS10/GS15 User Manual. Available online: http://www.surveyequipment.com/ PDFs/Leica_Viva_GS10_GS15_ User_Manual.pdf (accessed on 1 February 2019).
  49. IEC. Maritime Navigation and Radiocommunication Equipment and Systems—Digital Interfaces—Part 3: Multiple Talker and Multiple Listeners. High Speed Network Bus; IEC: Geneve, Switzerland, 2014. [Google Scholar]
  50. NMEA. National Marine Electronics Association. Available online: https://www.nmea.org /content/nmea_standards/v411.asp (accessed on 5 February 2019).
  51. Data Repository. Available online: https://github.com/lukisp2/ICP_DATA_MARINA (accessed on 16 August 2019).
Figure 1. The principle of measuring 2D horizontal coordinates ( x ,   y ) using Light Detection and Ranging (LIDAR).
Figure 1. The principle of measuring 2D horizontal coordinates ( x ,   y ) using Light Detection and Ranging (LIDAR).
Applsci 09 03530 g001
Figure 2. Mean error of the determination of position coordinates in the specified direction. (calculated using the covariance matrix).
Figure 2. Mean error of the determination of position coordinates in the specified direction. (calculated using the covariance matrix).
Applsci 09 03530 g002
Figure 3. Left side—2-point clouds before Dalaunay triangulation; right side—Dalaunay triangulation and nearest neighbor match (dashed lines).
Figure 3. Left side—2-point clouds before Dalaunay triangulation; right side—Dalaunay triangulation and nearest neighbor match (dashed lines).
Applsci 09 03530 g003
Figure 4. Weighting due to proposed error factor. Bolted, dashed line equals more weight for specific pair.
Figure 4. Weighting due to proposed error factor. Bolted, dashed line equals more weight for specific pair.
Applsci 09 03530 g004
Figure 5. Rejecting pairs due to robust criterion.
Figure 5. Rejecting pairs due to robust criterion.
Applsci 09 03530 g005
Figure 6. Left side—2-point clouds; right side—2-point clouds aligned.
Figure 6. Left side—2-point clouds; right side—2-point clouds aligned.
Applsci 09 03530 g006
Figure 7. A graphic presentation of the corrected position of the points matched in pairs using the sums of mean errors calculated along the lines connecting the matched points.
Figure 7. A graphic presentation of the corrected position of the points matched in pairs using the sums of mean errors calculated along the lines connecting the matched points.
Applsci 09 03530 g007
Figure 8. A graphic presentation the corrected position using the sums of mean errors CEP (circular error probable) calculated for the matched points.
Figure 8. A graphic presentation the corrected position using the sums of mean errors CEP (circular error probable) calculated for the matched points.
Applsci 09 03530 g008
Figure 9. The trajectory covered by the USV at the time of measurements.
Figure 9. The trajectory covered by the USV at the time of measurements.
Applsci 09 03530 g009
Figure 10. The distribution of measuring instruments on the unmanned surface vehicle (USV).
Figure 10. The distribution of measuring instruments on the unmanned surface vehicle (USV).
Applsci 09 03530 g010
Figure 11. (a) USV during a collection data phase; (b) wireless control system.
Figure 11. (a) USV during a collection data phase; (b) wireless control system.
Applsci 09 03530 g011
Figure 12. A graph of SOG ICP SOG R differences.
Figure 12. A graph of SOG ICP SOG R differences.
Applsci 09 03530 g012
Figure 13. Histogram representing the frequency of occurrence of SOG ICP SOG R differences.
Figure 13. Histogram representing the frequency of occurrence of SOG ICP SOG R differences.
Applsci 09 03530 g013
Figure 14. A graph of Δ COG ICP Δ COG R increment differences.
Figure 14. A graph of Δ COG ICP Δ COG R increment differences.
Applsci 09 03530 g014
Figure 15. A histogram representing the frequency of occurrence of Δ COG ICP Δ COG R increments.
Figure 15. A histogram representing the frequency of occurrence of Δ COG ICP Δ COG R increments.
Applsci 09 03530 g015
Figure 16. (a) Grid map computed by ICP with proposed directional weighting—HD_ICP; (b) grid map computed standard H_ICP.
Figure 16. (a) Grid map computed by ICP with proposed directional weighting—HD_ICP; (b) grid map computed standard H_ICP.
Applsci 09 03530 g016
Table 1. Summary of comparative statistical parameters.
Table 1. Summary of comparative statistical parameters.
MethodMean ValueStandard DeviationMaximum Value
S O G I C P S O G R ( m / s ) Δ C O G I C P Δ C O G R ( ° ) S O G I C P S O G R ( m / s ) Δ C O G I C P Δ C O G R ( ° ) S O G I C P S O G R ( m / s ) Δ C O G I C P Δ C O G R ( ° )
K_ICP−0.0310.080.1130.57−0.524−3.5
H_ICP−0.0130.050.0310.14−0.3190.72
HM_ICP−0.0120.030.0310.13−0.185−0.48
HD_ICP−0.0070.030.0260.11−0.180.42

Share and Cite

MDPI and ACS Style

Naus, K.; Marchel, Ł. Use of a Weighted ICP Algorithm to Precisely Determine USV Movement Parameters. Appl. Sci. 2019, 9, 3530. https://doi.org/10.3390/app9173530

AMA Style

Naus K, Marchel Ł. Use of a Weighted ICP Algorithm to Precisely Determine USV Movement Parameters. Applied Sciences. 2019; 9(17):3530. https://doi.org/10.3390/app9173530

Chicago/Turabian Style

Naus, Krzysztof, and Łukasz Marchel. 2019. "Use of a Weighted ICP Algorithm to Precisely Determine USV Movement Parameters" Applied Sciences 9, no. 17: 3530. https://doi.org/10.3390/app9173530

APA Style

Naus, K., & Marchel, Ł. (2019). Use of a Weighted ICP Algorithm to Precisely Determine USV Movement Parameters. Applied Sciences, 9(17), 3530. https://doi.org/10.3390/app9173530

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop