GNSS-Free Maritime Navigation using Radar and Digital Elevation Models

Modern maritime navigation is heavily dependent on satellite systems. Availability of an accurate position is critical for safe operations, but satellite-based navigation systems are vulnerable to interference, jamming, and spoofing. In this work, we propose a method for maritime navigation independent of GNSS, able to provide absolute positioning of the vessel based on marine radar scans. A measurement model is presented where a Digital Elevation Model is used to predict the output of a marine radar, given a hypothetical position. The model, as used by an on-line particle filter, is used to track the movements of a ship from real recorded data. This demonstrates the feasibility of this method for robust positioning, without the need of external positioning signals, in a maritime environment. The tracking only uses sensors commonly available on maritime vessels, and demonstrates its application using freely available elevation data.


I. INTRODUCTION
Positioning based on Global Navigation Satellite System (GNSS) has nothing less than revolutionized the localization of maritime vessels.However, its underlying technology is sensitive to disturbance and spoofing [1][2][3][4] and concerns are rising regarding the dependence on this single technology and its external transmitter-based approach.
It is thus of interest to find positioning solutions which do not rely on external technology, to increase the resilience against disturbance and adversarial actions.This would instead require a mapping of onboard sensor data to a global position, either as prior knowledge or by building a map on the go.The latter case is know as Simultaneous Localization And Mapping (SLAM) [5,6].SLAM is a family of algorithms which uses the sensor data as it goes to simultaneously build a map -and to position the sensor within the map built so far.Compared to using a Digital Elevation Model (DEM)-based model, SLAM could have the benefit of being less sensitive to the detailedness of underlying data while also adapting to changes more quickly.SLAM has been previously applied to radar data in a maritime setting [7].However, as there is no geographical reference to tie the SLAM map to the physical world, a SLAM algorithm would not be able to provide global positioning without additional localization sensors.Even with good initialization, SLAM is prone to drift unless combined with known global reference points.
Using traditional sea-charts, a similar application to the one presented in this article is presented in [8][9][10], with particle filter based maritime localization.Here, we instead propose the use of DEMs, for the following reasons.
Most importantly, elevation maps can generate more accurate predictions of the radar image.For instance, steep coastlines will give a much more distinct and larger radar return signal than a shallow beach.Moreover, hills and slopes further away from the coastline will also reflect the radar signal, which can explain radar echoes that cannot be predicted from a sea chart.Hence, a DEM model better predicts the multiple responses of each radar beam. Furthermore, • Sea charts generally only provide a binary map of sea vs. not sea.Hence, only one radar response can be modeled for each radar beam, representing the closest shoreline.• Coastlines are not static objects, they change over time.
Elevation maps generated by satellites can reflect these changes faster, and in a more accurate manner than traditional sea charts.• Elevation maps are global, while sea charts are most accurate in densely trafficked waters.In this paper we describe a system capable of localization in a known -but not necessarily previously visited -environment through the matching of DEMs with radar scans.The matching is performed using a particle filter utilizing a custom measurement model to model the expected radar response based on the DEM.This gives a method of localization which is resilient to interference, spoofing and jamming.
The system described in this paper is developed with, and tested against, recorded data from the ABB Marine and Ports sensor suite installed on a ferry in the Helsinki harbour.The data includes, but is not limited to, radar, electro-optical (EO) cameras, GNSS, inertial sensors and a digital compass.
The outline of the paper is as follows.In Section II, the background theory is presented along with an introduction of the data with which the implementation was made and tested.In Section III, the measurement model of the radar is derived, and a likelihood function is proposed based on the matching between modeled and measured radar responses.The implementation used for the evaluation of the proposed method is described in detail in Section IV.This is then used to demonstrate a real tracking scenario in Section V. Finally, Section VI concludes the paper.

II. BACKGROUND
In this section we review the background theory used in the work presented in this paper.

A. Filtering Theory
The basis for modern positioning is Bayesian state-space based filtering, in which the state x k at time k of the studied system is inferred from measurements z 1:k = {z 1 , . . ., z k } using the recursion The commonly known Kalman Filter (KF) [11] solves this recursion analytically for the linear Gaussian case.The Extended Kalman Filter (EKF) [12,13], which linearizes the system around the current estimate, has been successfully used to solve a more general class of problems.The Particle Filter (PF) [14] is a sequential Monte Carlo method to solve the Bayesian recursion.An overview of the method is presented here, see [15] for a more thorough description in the context of filtering.
The idea behind the PF is to use a number of random samples and weights (x i k and w i k ) to represent the sought distribution.This yields the approximation where N p is the number of samples, or particles, used in the approximation, all weights are positive w i k > 0, i w i k = 1, and δ x i k (x k ) is the Dirac delta function.The particle filter then updates the particles and weights over time to reflect the current information.In its basic form, also called the Bootstrap Filter, the filter boils down to the following steps all performed for i = 1, . . ., N p : • Initialization: Draw N p samples from the initial distribution, all with the same weight: • Prediction: Based on the expected motion, draw N p new particles to represent the state in the next time instance: • Measurement correction: Reweight the particles based on how well they fit the measurements: such that Np i=1 w i k = 1.• Resample: To avoid particle degeneration, i.e., all but one particle lose all their weight and importance, new particles are drawn from the old set with probability w i k , with replacement.There are two differences between the particle filter and the Kalman filter that makes the particle filter suitable for the studied application as outlined below.
Whereas the Kalman filter is limited to represent the state with a single mean and a covariance, the particle representation is more flexible allowing for several possible hypotheses to be perused simultaneously.Hence, resolving disambiguities in the estimated location can be delayed until more information is made available.

B. Marine Radar
Marine radar is today an indispensible tool for situational awareness on commercial ships, in particular in low-visibility environments.The marine radar is an instrument which measures electromagnetic (EM) reflectivity of its surroundings by active illumination with EM waves.Rotating around an axis, it yields samples at given angles, The raw radar returns have to undergo significant signal processing to make the relevant information accessible as imagery or individual detections [16].The result is the sum of all reflections from the illuminated areas as well as unwanted noise and interference phenomena speckle.
The returned signal strongly depends on the properties of the object causing the reflection, such as its reflectivity and its orientation in relation to the radar source.These properties are summarized in the Radar Cross-Section (RCS) [16], a measure of the object's reflectivity.The return signal power also is proportional to the inverse fourth power of the distance to the target, limiting the range of the radar [16].
Returns from an entire revolution yields a scan.In implementations, this is often stored as a matrix with rows and columns in accordance to angular and range resolution.However, when thresholded, this representation can be sparsified to yield individual measurements -reports -as 2D-points z c and their intensities, ι, In this paper, z c is defined in Cartesian coordinates centered around the ship radar.In each scan Z, N Z reports are detected.The radar's polar response is exemplified in Figure 1, with its Cartesian counterpart in Figure 2.

III. MEASUREMENT MODEL
This section outlines how the radar response is modeled from the DEM data, as needed to use the PF.For simplicity, only the position of the radar detections is considered, not the intensity.The model is separated in two parts; the first In the second part, further described below, the likelihood p(Z|x) = p(Z| Ẑ(x)) is derived.The former step generates potential radar responses, whereas the latter in essence defines a probability density over the similarity between two sets.In the paper, the orientation of the radar is assumed known, as measured by, e.g., the digital compass.If the orientation is not known, the state x could easily be extended to estimate the orientation.
The full pipeline of the resulting algorithm is illustrated in Figure 3 for a single particle position.The left column of the pipeline and the matching is then repeated for each particle in the particle filter.

A. Radar Response Model
The radar response model h (x) proposed in this paper bases its reports on the elevation surrounding the vessel, as given by a DEM.As radar measurements are taken in a radial fashion, we start by discussing how reports can be extracted for each individual angle, based on the elevation along a virtual beam radiating at an angle ψ from the radar, e ψ (r).Based on this elevation, we model the RCS at a given range.This singlebeam model is then applied in parallel to all angles reported by the radar.
The RCS in a given direction, ψ, and range, r, is approximated to be proportional to the area of the reflecting surface projected onto the plane perpendicular to the radar beam.This approximation makes it straightforward to compute σ ψ .Let d ψ (r) be the vector originating in the radar ending in the reflecting point on the surface, and e ⊥ ψ (r) the surface normal in the reflecting point.Then, as illustrated in Fig. 4, the modeled RCS follows from linear algebra as where •, • represents the scalar product and • the two norm.
To further model the radar's received response, the shadowing effect of the landscape is considered.To simplify the line of sight computations, all points along the direction of the ray are assumed visible if they are located higher than any preceding point otherwise not.Hence, the final RCS becomes This process is illustrated in Figure 5.
In the final step all points with sufficiently large RCS are identified as part of the report set.A threshold σ thres is applied to σψ (r) to obtain the set of modeled polar reports, where r max is the maximum range of the radar.This polar representation is then, through inverse polar transform, converted to Cartesian coordinates to give Ẑ(x), i.e., h (x).
Fig. 5: The model of the radar response in a single given direction is primarily a function of the incidence angle between the radar and the ground at range r.Here, σ thres = 0.

B. Point Cloud Matching
With the measured and modeled radar response available, Z and Ẑ(x), respectively, a measure of their correspondence can be created to provide the measurement likelihood needed in the PF measurement update (5).
Defining the difference between two point clouds is nontrivial.From the literature of point cloud registration, there are a number of error measures as summarized in, e.g., [17].One of the most basic measures is the mean (sometimes squared) distance between associated points.With unknown point-to-point association, a hypothetical association needs to be established to calculate the measure.An intuitive and, importantly, computationally efficient, method of association is the selection of the nearest neighbor of each point, from the other point set.In (12), NN z c k,j , Ẑ(x) is the nearest neighbor of z c k,j in Ẑ(x), by Euclidean norm.Based on the mean error, the following likelihood function is proposed: i.e., the mean distance of measured points to their nearest modeled point is considered Gaussian with zero mean and covariance σ 2 ptp .
Note that when comparing two point cloud sets, the association generally associates points from one set to points in the other, potentially leaving points in the second set without corresponding points.The question is hence whether to compare the measurements to the predicted points, vice versa or both.We found that due to the dominating cardinality of the modeled response, the one-sided matching of sensor data to the modelled points yielded a more stable tracking performance in this particular application -the map is in a sense more complete than the incoming radar scans.

C. Additional Models
More than modeling the radar response from the DEM, there are several other useful sources of information available on a ship, and more information to be gained from a DEM.
For example: radar returns not only come from natural sources -as modeled by the DEM -but also from other vessels in the vicinity.Commonly, nearby vessels report both position and size through the standardized Automatic Identification System (AIS) protocol.There are primarily two ways that ships can be included in the match between sensor data and the single-beam return model -either the returns can be exluded from the sensor data, or the returns can be included in the measurement model.For this, each ship could, for example, be modeled as an ellipsoid described by the reported AIS parameters, with an added approximate height.By adding this to the DEM prior to applying the radar response model, it would automatically be included in the beam model and detection matching.It should be noted, however, that AIS data is not fully reliable, and could potentially be subject to spoofing and jamming.
Furthermore; from the DEM, it is trivial to extract a virtual measurement of "being on water" contra being on land.Being a maritime vessel, it is natural to assume that the vessel is on water, i.e., not being on an elevation above water level as calibrated for tidal movements.Thus, given a tidal model, it is trivial to extract this data from a DEM, and to implement this as a particle filter measurement step, or as part of the prediction model.

IV. IMPLEMENTATION
In this chapter, we further describe the algorithm, and discuss practical aspects of the implementation that was made for the paper.The implentation of the particle filter and the likelihood function presented in Section III was made in the PYTHON programming language along with necessary supporting functions.For example -in the implementation, the tracking was performed with latitude-longitude state variables, but for each time update all particles were transformed into a local Cartesian north-east system, updated, then transformed back to latitude-longitude.

A. Radar Response Model
The DEM is given as a discrete map, expressed in this paper as the matrix D. The measurement model is most efficiently computed in the natural frame of the radar, using Fig. 6: Elevation maps provide further information compared to sea-charts -information which can be exploited to model the response from a radar sensor.polar coordinates.Hence, this data is transformed to polar coordinates to form the polar elevation matrix E, with one row per radar angle.The Cartesian map is illustrated in Figure 6.
The elevation for each angle ψ can then be defined as a lookup-function of the distance (range, r) from the radar, Here, ∆ ψ and ∆ r is the angular and range resolution respectively.
Pseudocode for the single-beam model from Section III-A is given in Algorithm 1.Note that the projection of ( 8) is performed with the trigonometric identity of tan (α Algorithm 1 Single beam model pseudocode Input: Matrix row i, RCS threshold, σ thres .The result from the algorithm applied to all angles yields the polar model of Figure 7, which is subsequently converted to Cartesian coordinates as in Figure 9a to match the format reported from the sensor.In the figures, the model is overlaid with the measured radar response to reveal that the model, based on the available elevation data, significantly overmodels the number of responses.This overmodeling, however, is counteracted in the choice of error measure in the point cloud matching between the sensed and modeled points.
In the testcases, the resolution of the DEM is 2-by-2 meters, limiting the accuracy of the positioning.Further error-sources include the rounding-effects of the polar transform of the discrete map.

B. Point Cloud Matching
The pseudocode implementation of the point cloud measure used in this paper is presented in Algorithm 2, where the measurement likelihood of (11) applied to the joint set of expected measurements.
The likelihood function of ( 11) is illustrated in Figure 8 for the example scan in Figure 9a over an m north/east area with a gridded 5 m resolution.It is clear that in this particular scan and DEM, a more or less nearby position provide an even better match for the received radar data.The corresponding match is illustrated in Figure 9b.A likely reason for this is the fact that the available elevation data does not include the houses along the shore which dominates the radar response within this particular harbor.This illustrates a weakness of the DEM approach compared to SLAM solutions which inherently is more robust to dicrepancies by not relying on external maps.

V. EXPERIMENTS
The algorithm was verified using data from a ferry running in the Helsinki harbor and archipelago.The ferry is equipped with Double Acting Technology, meaning its forward/aft direction can be fully reversed, leading to some additional considerations in algorithm implementations.Further, the ship has been retrofitted with a sensorsuite including radars, lidars, cameras and GNSS.These are connected through a backbone network to an onboard server for algorithm development.For the development presented in this article, multiple passes were recorded and post-processed using a real-time communication framework.This framework enables the same implementation that is developed offline to be used online without modification.Online tests are planned in the future.The marine radar in use performs pre-processing of the radar response.The radar data is sparsified by extracting individual points from a thresholded radar response, with a The particle filter was run with 100 particles initiated at first around a known GNSS position, with an added Gaussian noise of with a 10 m standard deviation in the north-and east direction.Further, a zero-velocity motion model was selected -not to optimize performance but to highlight the measurement model.In each timestep, process noise with covariance 7.2∆ t m2 was added to the particle filter, given the ∆ t time lapsed since the last measurement update.The RCS threshold was set at σ thres = 0.05, and σ 2 ptp = 5 m. Figure 10 shows an example run of the proposed algorithm.We considered 120 s of navigation, with the ferry cruising in south-southeast direction.In the figure, the reference GNSS track and the estimated position are shown on the map in red and yellow, respectively.The estimation error and the particles' standard deviation of spread are shown in the lower plots.The estimation error, post measurement updates, is kept below approximately 10 m to 15 m throughout the majority of the evaluation.
In this case, the tracking was performed on every 20th 3 Hz radar measurement, resulting in a runtime of approximately real time on a consumer computer.The runtime is dominated by the point cloud matching, indicating where future optimizations could be focused.

VI. CONCLUSION
In this work, we propose a method for maritime navigation independent of GNSS, able to provide absolute positioning of the vessel based on marine radar scans.We describe the method and demonstrate the feasibility to use of Digital Elevation Model (DEM) to model the reponse of a shipmounted marine radar.Tracking of the vessel is demonstrated using real data from Helsinki harbor.The quality of the tracking is, without surprise, observed to be dependent on the detailedness of the DEM.In particular, the exclusion of man-made structures in the DEM -in this case buildings along the harbor of Helsinki -sometimes skews the results by not providing a major source of detections present in the sensor data, instead comparing it to the more subtle coastline.Nevertheless, the model is capable of tracking the vessel in the harbor with a tracking error comparable to that of Global Navigation Satellite System (GNSS), maintaining an error in the harbor of less than about 10 m to 15 m throughout the majority of the testcase using only 100 particles.
At this point, the model does not factor in ground properties such as varying reflectivity, as only the position of expected returns are considered.
The implemented particle filter using the proposed likelihood was shown to be capable of real time tracking of the vessel.Real time tracking enables robust localization.This provides a redundant system to complement GNSS navigation.It can also be used to detect GNSS malfunctions or attacks, providing immunity to interference, jamming, and spoofing.The tracking performance could, with optimizations and increase in computing power, further be increased by including more of the radar scans at a higher rate or increasing the number of particles in the particle filter.

Fig. 1 :
Fig. 1: Radar reponse is sampled in polar coordinates centered at the radar.

Fig. 2 :
Fig.2: When the radar response is converted to Cartesian coordinates, the shape of the harbor becomes apparent to the human eye.

Fig. 3 :
Fig. 3: Measurement model pipeline.Figures and explanations thereof are found throughout the article.

Fig. 4 :
Fig. 4: Here, the relevant quantities of the radar reflection model are illustrated, with d ψ (r) being the vector originating in the radar ending in the reflecting point on the surface, and e ⊥ ψ (r) the surface normal in the reflecting point.

Fig. 7 :
Fig. 7: Measurement model with the corresponding polar radar response overlaid in red.

Algorithm 2 Fig. 8 :
Fig.8: The likelihood function from(11) evaluated around the reference GNSS position, in a 5 m resolution grid.In the plot, the GNSS reference position (white) as well as a local maxima (green) are marked.
Fig.9: In the likelihood plot in Figure8, the marked positions correspond to the matches illustrated above.The modeled response is plotted in black and the measured radar response in red.

Fig. 10 :
Fig. 10: Track and error of the particle filter tracker, as compared with the GNSS reference.The reference GNSS track and the estimated position are shown on the map in red and yellow, respectively.