Spatiotemporal Geomagnetic Intensity Mapping Using Gaussian Process for Underwater Vehicle Localization

Article information

J. Ocean Eng. Technol. 2026;40(2):227-237

Publication date (electronic) : 2026 April 1

doi : https://doi.org/10.26748/KSOE.2025.073

Dongwook Lee¹, Hyunkeun Cho², Jinwhan Kim³

¹Graduate Student, Robotics Program, Korea Advanced Institute of Science and Technology, Daejeon, Republic of Korea

²Graduate Student, Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Daejeon, Republic of Korea

³Professor, Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, Daejeon, Republic of Korea

Corresponding author Jinwhan Kim: +82-42-350-1519, jinwhan@kaist.ac.kr

Received 2025 December 2; Revised 2026 January 23; Accepted 2026 February 26.

Abstract

1. Introduction

Underwater exploration is a central pillar of ocean science and marine engineering, supporting numerous downstream applications including subsea cable and pipeline inspection, offshore renewableenergy surveys, marine ecosystem monitoring, and the prediction of seafloor and sub-seafloor processes. Although autonomous underwater vehicles (AUVs) have become increasingly capable, achieving robust, long-duration navigation without access to the global navigation satellite system (GNSS) remains a major research challenge. Industrial systems typically rely on underwater acoustic positioning technologies ―long-baseline (LBL), short-baseline (SBL), and ultra-short-baseline (USBL)―which can provide centimeter-level accuracy. Nevertheless, these systems require a pre-deployed infrastructure or external acoustic beacons, limiting scalability and constraining deployment in remote or dynamic environments. This gap motivates the development of fully sensor-based navigation frameworks.

The geomagnetic field, generated mainly by the dynamo process within the Earth’s outer core, permeates the planet and serves as a passive, naturally available signal. Although acoustic sonar has traditionally been the primary modality for underwater navigation, early work, such as Tyren (1987), first suggested that geomagnetic-field measurements may replace or complement terrain-based maps for underwater vehicle navigation. Unlike acoustic systems, magnetic sensing is entirely passive, enabling covert operation without disturbing marine life or revealing the presence of a vehicle in defense applications.

In addition to the global geomagnetic field, ferromagnetic structures within the Earth’s crust and anthropogenic installations generate local magnetic anomalies, which are spatially distinct signatures that can be exploited for mapping and localization. This observation has motivated extensive research on magnetic-field–based localization in indoor and enclosed environments. Early studies have shown that magnetic distortions from steel structures can serve as stable fingerprints for indoor localization. Siitsakulchai et al. (2000) used a two-axis electronic compass in a corridor environment to construct a magnetic map from heading discontinuities without interpolation. Haverinen and Kemppainen (2009) expanded this approach by constructing a building-scale magnetic map using three-axis magnetometer data and Monte Carlo localization. They reported that magnetic vectors and magnitudes remain stable over week-long intervals. Vallivaara et al. (2011) then introduced one of the first magnetic-field simultaneous localization and mapping (SLAM) systems, using a Rao–Blackwellized particle filter with a Gaussian process (GP) map with a radial basis function (RBF) kernel. Their method successfully captured the smooth spatial variations but assumed a static magnetic environment, neglecting temporal changes.

Compared to indoor magnetic SLAM, magnetic-field mapping in underwater or near-water environments has received far less attention. Jung et al. (2018) reported autonomous maritime magnetic mapping using an autonomous surface vehicle (ASV) equipped with a rigid tri-axial magnetometer. Their framework included comprehensive calibration, including soft-iron and hard-iron correction via 3D Kabsch alignment, and used GP regression with a spatial RBF kernel to construct a smooth geomagnetic intensity map. Nevertheless, the model remained purely spatial and did not incorporate temporal variability.

Although magnetic anomalies are widely regarded as stable spatial signatures, most existing work treats the geomagnetic field as static or quasi-static. Despite this, geophysical studies have shown that the geomagnetic field undergoes measurable temporal variations across multiple time scales, including diurnal oscillations, sub-daily disturbances, and storm-time fluctuations. Wardinski and Holme (2006) developed a time-dependent global geomagnetic model for 1980–2000 using cubic B-splines to represent Gauss coefficients, capturing long-term secular variation and geomagnetic jerks but not local spatial structure. Sutcliffe (2000) used feed-forward neural networks to model the diurnal variations in Southern Africa using inputs such as local time, day of year, latitude, and sunspot number, successfully reproducing daily and seasonal behavior but without modeling spatial interactions. More recently, the effective range of the station mapping (ERSM) model (Srinivasan & Nielsen, 2025) extended the temporal correction range of reference stations using regression methods, including linear models, k-nearest neighbors, and neural networks. This study focused exclusively on the temporal corrections without addressing spatial mapping.

Despite these advances, the temporal dynamics have not been integrated into local maritime magnetic mapping. No prior work has jointly modeled the spatial structures of the underwater geomagnetic intensity with diurnal or short-term temporal variability within a unified probabilistic framework. This paper addresses this gap by introducing a GP-based approach that models the spatial dependence across neighboring marine regions and the temporal variation induced by solar activity―whether cyclic or transient―directly through kernel design.

2. Maritime Geomagnetic Field

Earth’s geomagnetic field extends continuously through underwater, terrestrial, and atmospheric environments, but the characteristics of geomagnetic anomalies, which are defined as deviations from the global geomagnetic model, differ significantly across these environments. In indoor or terrestrial settings, anomalies are dominated by strong, spatially irregular signatures from ferromagnetic structures, electrical infrastructure, and heterogeneous near-surface geology. These shallow and human-made sources produce magnetic maps with sharp gradients and rapid spatial fluctuations over meter-scale distances.

By contrast, maritime magnetic anomalies arise primarily from deep-sea geological structures, with few ferromagnetic disturbances in the water column or on the sea surface. As a result, maritime geomagnetic intensity maps exhibit a smooth, low-frequency spatial structure with far more coherence than terrestrial magnetic maps. Fig. 1 supports this claim by showing that the urban magnetic intensity (right figure, overall) exhibits greater variance than the maritime magnetic intensity (right figure, within the blue boundary).

Fig. 1

Location and magnetic intensity map across the study area. The magnetic field was measured at an altitude of 1000 m during the aeromagnetic surveys. On the right figure, Red lines = country borders, blue polygon = Lake Kinneret, Il = Israel, Sy = Syria, Jo = Jordan (Schattner et al., 2019).

Many studies have relied on surface vessels or aerial platforms to map the maritime geomagnetic intensity because the measured geomagnetic field can be expressed as a function of distance using a mathematical transformation known as continuation. Nevertheless, the downward continuation of potential fields is fundamentally ill-posed, whereas magnetic anomalies can be continued reliably upward. This means that, in practice, maps collected at the sea surface or at higher altitudes lack the high-frequency components that exist near the seafloor. Therefore, underwater vehicles must conduct close-range magnetic surveys to capture the spatial details that would otherwise remain unresolved.

Maritime magnetic fields are also strongly influenced by natural temporal variations. Diurnal variations driven by solar radiation continuously modulate the geomagnetic field, while short-term disturbances arise from solar-wind activity and geomagnetic storms. An AUV may traverse large regions over extended periods. Therefore, explicitly modeling these temporal variations is essential for maintaining spatial consistency and achieving accurate regional geomagnetic mapping. Fig. 2 shows this behavior using real geomagnetic data from the Cheongyang observatory in the Republic of Korea, where the three-axis magnetic field exhibits clear diurnal patterns superimposed with short-term fluctuations.

Fig. 2

Geomagnetic Field Measurements at the Cheongyang INTERMAGNET Observatory. The plot displays geomagnetic intensity values along the x, y, and z axes, represented in red, green, and blue colors, respectively, with the normed intensity shown in black. The data was recorded from August 1^st to 7^th, 2025, at a sampling rate of one sample per minute.

3. Stochastic Variational Gaussian Process

GP provides a flexible Bayesian framework for regression by defining a distribution over functions with closed-form expressions for the marginal likelihood and predictive posterior. This enables automatic hyperparameter selection by maximizing the log marginal likelihood and yields well-calibrated uncertainty estimates. Nevertheless, the exact GP inference scales cubically with the number of training points, yielding a computational complexity of O(N³) (big O notation), making it impractical for large datasets.

This study addressed this limitation using sparse variational approximations with a smaller set of M inducing variables (where M ≪ N) to summarize the latent function, with N representing the total number of data points. Classical sparse GP regression (SGPR) analytically integrates out the optimal variational distribution and provides a deterministic evidence lower bound (ELBO), offering strong approximation quality but requiring full-batch computations.

Stochastic variational GP (SVGP) extends this framework to large-scale learning by replacing the analytic integration with an explicit variational distribution q(u) = N(m,S ) over the inducing variables u, where m and S are the mean vector and covariance matrix, respectively. This form allows the ELBO to factorize across data points, enabling unbiased stochastic gradient estimates using minibatches. As a result, the computational cost per update reduces to O(M³+ MB ), where B is the minibatch size, making SVGP suitable for datasets far too large for full-batch methods.

The flexibility of SVGP comes at the cost of a more complex optimization problem. The variational parameters (m,S), inducing locations, and kernel hyperparameters must be jointly optimized, and the stochastic nature of the gradients slows convergence relative to deterministic SGPR. Empirical studies have shown that the SVGP performance depends sensitively on the training hyperparameters, such as the minibatch size, learning rate schedules, and optimizer momentum, making it challenging to define a universally reliable, fully automatic training procedure.

Nevertheless, SVGP is a widely used approximation method because of its scalability and ability to accommodate non-Gaussian likelihoods. Its development represents an important milestone in making GP-based inference feasible for modern large-scale machine learning scenarios where deterministic variational approaches are computationally prohibitive.

4. Gaussian Kernel Design for Spatiotemporal Geomagnetic Intensity Mapping

The Gaussian kernel design for geomagnetic intensity mapping requires representing the spatial and temporal characteristics of the measured magnetic field. The temporally varying geomagnetic field contains three primary components: a spatially varying regional anomaly produced by seabed geology, a diurnal temporal variation generated by ionospheric currents and solar radiation, and short-term temporal fluctuations caused by solar–wind disturbances and geomagnetic storms. The GP covariance function is constructed as an additive kernel with distinct spatial and temporal terms because the measured magnetic intensity is the sum of these effects. The measured magnetic field B at a certain location x on a specific time t is expressed as B (x, t) = B_spatial(x) + B_diurnal(t) + B _short(t), which corresponds to the decomposition shown in Fig. 3.

Fig. 3

Composition of the measured magnetic field intensity. The plot illustrates the breakdown of the geomagnetic intensity, with the spatial intensity (in blue) representing the baseline geomagnetic field, the diurnal variation (in orange) capturing the daily cyclical changes, and the short-term variation (in black) reflecting abrupt fluctuations.

5. Experiment

5.1 Geomagnetic Field Map Emulation

The performance of the designed GP kernel was evaluated by emulating the test environment in MATLAB, which generated synthetic maritime magnetic field data resembling real-world conditions. The first step is to generate two random terrains using a diamond-square terrain generator, where each terrain is produced with a user-specified height mean, variation, and roughness parameters. One of the generated terrains is assigned as the upper bound, and the other as the lower bound, forming a set of three-dimensional prisms over the x, y grid. After the prism geometry is defined, the user specifies the physical density and the magnetization vector components in th e x, y, and z directions for each patch within the region. The final magnetic intensity map is derived by integrating the dipole field over irregular prisms using the assigned density and magnetization values for every patch. This procedure produces a spatially continuous magnetic map with a controllable geological structure and magnetic properties; the overall map-generation process is illustrated in Fig. 4. The AUV was assumed to operate at a fixed depth of 50m below the water surface, setting it as the z-location of observation points and mesh, because this study did not include the scope of downward or upward continuation.

Fig. 4

Magnetic field map generating emulation. Two terrains serve as upper and lower target boundaries, where the density and magnetization of the material are assigned within selected X-Y patches. Dipole integration at a certain observation height outputs the final magnetic field map (Arroyo et al., 2015).

5.2 AUV Mission Area and Trajectory Plan

In the emulator, a 10 × 10 km grid map composed of 100 × 100 cells was generated, where each cell represents a 100 × 100 m grid resolution. This grid resolution is consistent with common practices in marine magnetic surveys, where grid sizes of ~100 m are often used for regional-scale magnetic mapping (Saltus, 2023). Furthermore, regions within a ~100 km radius exhibit almost identical temporal variations in geomagnetic data, allowing for the superposition of a consistent temporal signal across the map area.

The AUV trajectory was defined using a left-to-right lawnmower pattern followed by a top-to-bottom pattern. This design ensured that the vehicle revisits identical spatial locations at different times, enabling spatiotemporal interpolation with sufficient coverage. Fig. 5 presents the generated geomagnetic intensity map and the planned path of the vehicle, based on the lawnmower pattern, for the sensing emulation.

Fig. 5

Generated geomagnetic intensity map and the planned lawnmower path. (a) shows the generated geomagnetic intensity map of the 10 × 10 km area, while (b) shows the planned vehicle path, following the lawnmower pattern, with red and green markers indicating the start and end points of the trajectory.

The mission area was scaled to various sizes ranging from 1 km² to 100 km² to show the scalability of the spatiotemporal SVGP across different mission time scales and map sizes. The test mission areas selected were 1 km², 9 km², 25 km², 49 km², and 100 km², and these map sizes are shown in Fig. 6. The maps were generated by cropping the raw map into the target size (rather than scaling) to preserve the spatial gradient resolution of the magnetic intensity map.

Fig. 6

Different scales of mission areas. Given a simulated magnetic intensity map, regions of varying sizes are cropped to generate maps corresponding to different mission time scales. The diagram shows mission areas of sizes 1 km², 4 km², 9 km², and 100 km², reflecting the scale of the data used in the validation process.

As an exception, the effects of temporal variations during a shorter mission period were tested by scaling down the 1 km² map to 0.25 km², effectively halving the mission period. The path shown in Fig. 5(b) was scaled to fit within each of these mission areas, maintaining the same pattern. Each map was assigned a unique map ID for evaluation, as shown in Table 1, along with the mission details.

Table 1

Map identification (ID) and mission details for different map sizes

5.3 Magnetometer Reading Emulation

A three-step process was followed to emulate magnetometer readings: (1) bilinear interpolation from a simulated grid map to obtain raw intensity values, (2) adding sensor noise based on specified sensor characteristics, and (3) incorporating temporal variations based on real geomagnetic data from the INTERMAGNET Data Viewer.

Magnetometer readings at each vehicle position were first obtained through bilinear interpolation from the underlying intensity grid. The SeaSpy2 maritime magnetometer model, which has a sensitivity of 0.01 nT and an absolute accuracy of 0.1 nT, was used to reflect real-world sensor characteristics. After obtaining the raw intensity values, the sensor noise (modeled as Gaussian noise with a standard deviation of 0.01 nT) and a bias of 0.1 nT were added to simulate the sensor’s behavior more accurately.

One week of 1 Hz geomagnetic data was retrieved from the INTERMAGNET Data Viewer to incorporate realistic temporal variation. Superimposing a single temporal signal across the entire simulated area was feasible because regions within a ~100 km radius exhibit similar temporal variations. A random mission time interval was then selected and added to the synthetic readings, mimicking the temporal behavior during the mission. Fig. 7 presents the pipeline of extracting temporal variation to be applied to the magnetic intensity measurements, while Eq. (5) describes the final mathematical notation of the emulated magnetometer reading (m: magnetic intensity):

Fig. 7

Pipeline explaining the extraction of the geomagnetic intensity temporal variation. The data retrieved from the INTERMAGNET Data Viewer was normalized and randomly windowed to produce the final temporal variation, which is then added to the simulated spatial geomagnetic readings.

(5) mmeasured(x, y, t)=mraw(x, y)+mnoise+mbias+mtemporal(t)

where

• m_raw(x, y) = B ilinear Interpolation(map, x, y),

• m_noise ~ N(0, 0.01), m_bias = 0.1,

• m_temporal (t) = INTERMAGNET(t)

A magnetic intensity measurement was taken at a 1Hz rate while the AUV traveled at 6 knots (≈ 3.1 m/s). Fig. 8 presents an example of the magnetic intensity measurements along the travel path for Map F before and after the addition of temporal variations.

Fig. 8

Example of emulated magnetometer readings before (a) and after (b) the addition of temporal variation. This figure presents an example of a map 100 km² in size (Map F).

5.4 Overall System Architecture

In this study, it was assumed that the vehicle follows the planned trajectory without any drift, as depicted in Fig. 9(a). This assumption simplifies the modeling process, allowing the study to focus on the core functionality of the geomagnetic intensity mapping. The trajectory and position of the vehicle are predefined, and the mapping module operates based on this idealized trajectory. In practical applications, however, the mapping module must be integrated with the inertial navigation system (INS) of the vehicle, which provides real-time estimates of the position, velocity, and orientation of the vehicle. This integration allows the system to account for real-world conditions, such as sensor noise, drift, and dynamic environmental factors, as shown in Fig. 9(b). The INS ensures that the state information of the vehicle is updated continuously during the mission, enabling more accurate and robust geomagnetic intensity mapping in real-time.

Fig. 9

System architecture of the geomagnetic intensity mapping module integrated with the navigation system of the vehicle. (a) Simulated pipeline where the vehicle follows a predefined trajectory without drift. (b) Real-world implementation pipeline where the inertial navigation system provides continuous updates on the position and motion of the vehicle.

6. Results

On each test map (of different sizes), a map regression test was conducted using two types of kernel designs: spatial kernel (using the RBF function only), and the spatiotemporal kernel (proposed, adding up three distinct functions). Training was done with 1,000 inducing points per map, with 300 epoch iterations. The GPU used for the training was a NVIDIA GeForce RTX 3090. Table 2 presents a quantitative comparison of using two kernel designs by calculating the RMSE on each map compared to the ground truth (GT) map, and Table 3 compares the computation time taken for both kernel designs on different map sizes.

Table 2

RMSE between the GT map and the estimated maps (nT)

Table 3

Computational cost (s)

After reviewing the RMSE values presented in Table 2, the map accuracy improved, starting from Map B (1 km² size, 4.33 h), where the accuracy consistently fell within the 2–8 nT range. For smaller map sizes (corresponding to shorter mission durations), however, the RBF function alone produced more accurate map estimates. Accordingly, for missions shorter than approximately 4 h, the temporal variation regression introduced by the spatiotemporal kernel is not beneficial and may even degrade the map accuracy. Regarding the computation time, a simpler RBF function results in faster training times because of its reduced complexity (Table 3). For well-working map sizes (Maps B to F), the proposed spatiotemporal kernel took approximately 1–4 min longer to calculate than the RBF kernel alone (5.6%–9.8% longer), while providing better map results.

These results clearly highlight the trade-off between the mission duration and the complexity of the kernel design. For shorter missions (i.e., smaller maps), the spatiotemporal added complexity of the kernel does not provide significant benefits and leads to longer computation times. By contrast, the RBF kernel offers faster training times and more accurate results in scenarios with relatively brief mission durations (less than 4 h).

For longer missions (i.e., larger maps), the proposed spatiotemporal kernel yields more accurate map estimates, justifying the slightly increased computation time. This trade-off between mission duration (which determines the map size and complexity) and the modeling complexity (which influences the computational cost) is crucial for understanding the practical applicability of the proposed method. The spatiotemporal added complexity of the kernel becomes more beneficial as the map size and mission duration increase, where the improved accuracy outweighs the additional computational cost. Fig. 10 qualitatively shows the estimated maps obtained with two kernel designs along with the GT map.

Fig. 10

Comparison of the GT map with the RBF-only kernel and the proposed kernel designs.

The estimated values were analyzed by plotting the error and variance of the GP estimate on the training points as a function of time, as shown in Fig. 11. In the first row of the plot, which shows the error between the GT and the RBF-only kernel estimate, the magnitude of the error increased significantly as the map size and mission time increased when using the RBF-only kernel. In addition, the estimated variance also increased. In the second row, however, where the proposed kernel was used for spatiotemporal regression, the error remained within a 20 nT bound, and the GP variance was controlled. Therefore, this kernel is stable in handling spatiotemporal regression problems.

Fig. 11

Error between the GT and the estimated magnetic intensity values at specific (x, y, t) locations, flattened along the time axis. The first row displays the error plot between GT and RBF-only kernel-based estimation, while the second row shows the error plot between the GT and the proposed kernel-based estimation.

Further analysis of the efficacy of this kernel can be conducted by estimating the temporal variation itself. Inputting a specific time t into the temporal kernel k_temporal of the trained model allows estimates of how well the temporal component of the geomagnetic intensity is captured. Fig. 12 shows the estimated temporal variation and the error magnitude for each map ID (A–F), using the k_temporal kernel of the proposed model.

Fig. 12

Estimated temporal variation and error magnitude for the proposed kernel-based estimation, showing the GT and estimated intensity at specific (x, y, t) locations, flattened along the time axis. The first row shows the estimated temporal variation (orange), while the second row shows the corresponding error magnitude (red).

7. Conclusions

This paper presents a GP-based regression framework for reconstructing maritime geomagnetic intensity maps, effectively capturing the spatial anomalies and temporal variations in geomagnetic signals. The proposed method integrates an additive kernel that combines an RBF spatial term, a periodic diurnal term, and a Matérn-1/2 short-term term, allowing accurate modeling of the key physical characteristics of geomagnetic fields.

A comprehensive simulation pipeline, which includes prism-based synthetic geomagnetic map generation, temporal variation injection using real-world INTERMAGNET data, and AUV trajectory and sensing emulation, was developed to evaluate the effectiveness of the proposed framework. The proposed spatiotemporal kernel outperformed the traditional RBF-only GP regression, particularly when measurements were obtained at different time intervals. Incorporating temporal variations mitigates the inconsistencies at revisited locations and ensures that regression errors remain within 10–20 nT. The final geomagnetic maps exhibit smooth, geophysically realistic structures with well-calibrated uncertainty estimates, maintaining an RMSE of 2–8 nT.

These results underscore the critical importance of considering temporal variation in real-time geomagnetic mapping, particularly for AUV navigation and localization. The proposed method offers a powerful tool for long-duration missions, where a precise reconstruction of the spatiotemporal geomagnetic field is essential. Future work will explore extending the model to handle larger, more complex mission areas and integrating it into real-time AUV navigation systems.

Article information Continued