Underwater Sonar Fusion for Enhanced 3D Reconstruction using Background Data and Voxel Initializer Mapping for AUV

2.1 Object Detection and Voxel Map Initialization

Creating a local voxel map across an entire environment is highly inefficient in terms of data storage and computational load. As observed in OctoMap [10], sparse voxel allocation significantly improves memory efficiency and computational performance by generating voxels only at required locations. Based on this principle, this study adopted an active object detection approach to implement algorithms using PS sensor data, generating local voxel maps only when objects are detected, instead of maintaining complete maps throughout the scanning process.

The PS rotated about the vertical axis and recorded the sonar signal intensity in the direction of the sensor faces. Although a full 360-degree rotation is generally possible, this study scanned only within a limited angular range to solve the problem of excessively long scan cycles during full rotation.

The object-detection algorithm using PS data is shown in Fig. 4. The data format in Fig. 4 (a) is from a single PS sequence, with the echo intensity recorded at all points within a distance of 0 m to the maximum range. The highest intensity point (marked with a green square) was identified as the selected point for that sequence. These selected points were then transformed from the local coordinate system of the sensor to global 3D coordinates by combining the position and orientation data of the AUV with the relative distance and angle measurements from the PS sensor. As shown in Fig. 4 (b), based on a repeated increase of the Z-coordinates of the selected points from each sequence, the system determines whether an object has been detected. In this process, when the PS detects an object, the system initializes a local voxel map with center at the location of the detected object. The voxel map covers a 10 × 10 × 10 m³ region with 100 × 100 × 100 resolution, providing a 10 cm voxel resolution for precise 3D mapping.

Fig. 4.

Object detection algorithm using PS data; (a) Point selection method in single sequence; (b) Object recognition pattern across multiple sequences.

2.2 FLS Data-based Voxel Map Update

Once the voxel map is initialized through object detection, the system begins to update the voxel map using the FLS data. FLS sensors provide 2D image data, which are fundamentally acoustic data in polar coordinates converted into Cartesian coordinate images, as illustrated in Fig. 5. Here, the vertical pixel information corresponds to the distance from the sensor, and horizontal pixel information includes angular information corresponding to the individual beams.

Fig. 5.

FLS image coordinate conversion.

This system does not convert all pixel data but selectively processes only pixels with high probability in object detection. Specifically, only pixels with intensity values exceeding the threshold of 150 (in the 0–255 intensity range) were considered as object detection data, and only these selected pixel indices (i, j) were processed for further computation. By selecting only high-intensity points, we accurately reconstructed the objects and seabed.

However, the pixels lack precise vertical information because of the inherent ambiguity of FLS. In Fig. 2 (b), the dark-gray shaded area represents this uncertainty zone, where the exact vertical position cannot be determined from a single measurement. To address this, we first converted the selected pixel indices (i, j) from Fig. 5 to physical measurements (r_i, θ_i) using Eqs. (1) and (2), before including the geometric parameters of tilt angle φ_t1, beamwidth φ_v, and FoV θ_b. Coordinate conversion formulae were designed to address this vertical ambiguity.

Eqs. (1) and (2) convert image indices into physical measurements. Eq. (1) converts the range information from FLS image coordinates to actual distance measurements. Here r_r is the distance from the sensor to the detected point. r_min and r_max are the minimum and maximum measurement distances of the sensor, respectively; h is the total number of vertical pixels; and r_i is the range index of the image.

Eq. (2) converts the horizontal angular information from image coordinates into beam angles; θ_F represents the horizontal FoV; θ_i is the angular index of the image; and ω is the total number of beams (horizontal pixels).

These converted measurements are then used for the 3D coordinate transformation, as illustrated in Fig. 2 (b). The vertical ambiguity is addressed by generating sample points across the vertical FoV as in Eq. (3), where k represents the k-th possible vertical position; φ_min and φ_max define the vertical scan angle range; and N_v is the number of sampling points. This equation distributes possible object positions across the entire range of vertical ambiguity.

Eq. (4) transforms each 2D image coordinate (i, j) with multiple vertical k samples into the final 3D coordinates. The parameter k generates sample points distributed across the vertical ambiguity range, whereas φ_t1 represents the FLS mounting tilt angle, which is a fixed sensor installation parameter.

Based on the reconstructed 3D points, the occupancy grid map is updated using the log-odds method. For each detected 3D point P = (x, y, z, intensity), we calculated the weight only for the voxel containing this 3D point.

In Eq. (5), used to calculate ΔL_o, P_o is an adjustable parameter that determines the strength of the occupied space updates. To assign a stronger weight to occupied regions, P_o is set closer to 1.0, whereas values closer to 0.5, provide more conservative updates.

For the FLS data processing, we updated only the voxels containing the detected 3D points in occupied spaces. Eq. (6) is the update formula where ΔL_o is the update weight; α_F adjusts the influence between PS and FLS updates; and I_F is the intensity value normalized to between 0 and 1.

The calculated L(x, y, z) values are limited to the upper and lower bounds expressed in Eq. (8) for numerical stability, where the clamp function restricts the values to a specified range, as in Eq. (7). In Eq. (8), this clamp function constrains the log-odds values to the range [-5.0, 5.0], to prevent numerical overflow and maintain computational stability during iterative updates.

Finally, the L(x, y, z) values are converted back to probabilities through a sigmoid function using Eq. (9), to calculate the final occupancy probability.

2.3 PS Data-based Voxel Map Update

As explained in Section 2.1, strong intensity signals in the PS sensor data are not necessarily definitive evidence that they originate from objects. Existing sonar data processing methods include using the maximum echo [11], selecting multiple peaks [12], and using all values above certain thresholds. While these approaches work well for general object detection, a more conservative approach is required for voxel map updates, considering the irreversible nature of the carving operation.

The critical difference lies in the purpose of data utilization. In Section 2.1, we focused on object detection, which tolerates false positives, whereas voxel carving permanently removes potential object regions. The most dangerous scenario occurs when an object is located nearby, but a distant seafloor produces stronger reflection signals. In intensity-based selection methods such as maximum echo, the system selects the distant seafloor point and incorrectly carves away the space with actual nearby objects.

Therefore, this study adopted a conservative closest distance-based selection strategy instead of conventional intensity-based methods. This approach ensures that actual objects are not carved away, prioritizing preservation over aggressive removal.

To implement closest distance-based selection on PS sensor data, we adopted a group-based selection method as illustrated in Fig. 6. The algorithm first applies a predefined threshold τ = 0.7 to select only high-reliability signals above this threshold (orange dashed line). These selected points are then sorted by distance and grouped, classifying consecutive points as belonging to the same group. Among all the formed groups, the group with the shortest distance measurements is selected. Within this closest group, the final detection point was determined by selecting the point with the highest signal intensity (marked with a black square). Through this process, a selected point was obtained from a single PS sequence.

Fig. 6.

Closest point selection algorithm from multiple echoes in PS data

The PS-based 3D coordinate conversion requires that the unique characteristics of the sensor be considered. As shown in Fig. 2 (a), the PS exhibits horizontal ambiguity (light gray fan-shaped area) while providing accurate vertical information, which is complementary to the FLS behavior. The geometric parameters including φ_t2, beamwidth θ_h, and profiling range φ_p define this sensing geometry. To handle the horizontal uncertainty, we distributed multiple candidate points across the ambiguous angular range at the selected distance r_s:

Eq. (10) describes the horizontal angle sampling used to handle the inherent horizontal ambiguity of PS; θ_min and θ_max are the horizontal FoV angles of the PS sensor; N_h is the number of sampling points distributed across this angular range; and k is the sampling index.

Eq. (11) transforms the PS data into 3D coordinates. Here, r_s is the chosen distance from the group-based selection method; φ_t2 is the PS sensor mounting tilt angle; φ_p is the PS profiling angle for the current measurement; and q_h(k) represents the k-th horizontal angle sample. Parameter k generates multiple candidate positions across the horizontal uncertainty range.

To identify free-space voxels between the sensor and detected points, we applied the 3D Bresenham algorithm [13], which extends the traditional 2D ray-tracing method to three-dimensional space. The Bresenham algorithm efficiently determined all the voxels along the straight-line path from the sensor position to each detected PS point. Because the acoustic signal successfully traversed this region without obstruction, these intermediate voxels represent the confirmed free space.

To each voxel identified by the Bresenham algorithm along the ray path , a weight update was applied as follows:

The ΔL_f is calculated according to Eq. (12), where P_f is an adjustable parameter that determines the strength of free space updates. To ensure conservative carving operations, P_f is typically set to less than 0.5, creating negative log-odds increments to mark voxels as free space.

Using Eq. (13), negative weight updates are applied with L_new(x, y, z) and L_old(x, y, z) representing the updated and previous voxel occupancy weights, respectively. ΔL_f is a negative value; α_P controls the carving intensity; and t represents the distance from sensor to point normalized to the 0–1 range to implement distance-based weighting. The exponential function precludes strong negative weights near the detection point due to potential noise while generating stronger negative weights at farther distances, where we can be more confident that no objects exist.

Through this PS weight update process, ambiguous regions predicted by FLS as potential object locations were systematically carved away owing to the PS's accurate vertical information, to complete the sensor fusion approach.

2.4 Voxel Map Completion and Object Segmentation

After completing the iterative voxel map updates from the FLS and PS sensors, the system obtains a comprehensive 3D occupancy probability map of the target region. This final stage realizes the core concept of “local mapping before object extraction.” The process consists of two main operations: probability thresholding and object detection.

For probability thresholding, the system converted a probabilistic occupancy map into a binary classification. Voxels with occupancy probability values above 0.7 are classified as occupied space, whereas those below this threshold are considered free space. This threshold value was determined by fine-tuning the occupancy weight parameters to achieve optimal reconstruction performance. The resulting binary map provides a clear representation of the reconstructed 3D scene.

Object detection applies a moving average-based approach to the reconstructed voxel map, following the same principle as that of the object detection algorithm described in Section 2.1. The system computes a moving average of the Z-coordinates from the occupied voxels to establish a robust seafloor baseline resistant to noise and small objects. When consecutive voxel heights deviated significantly from the moving-average baseline, the system identified these regions as objects. This moving-average filtering effectively suppresses noise and minor seafloor variations, while maintaining sensitivity to genuine objects of interest.

Through this four-stage process, the proposed method successfully demonstrates the “background-included mapping followed by segmentation” approach. Compared to existing FLS+PS methods that rely on ROI-based object detection before reconstruction, the proposed approach is expected to achieve superior reconstruction performance. Although this approach is effective in the current simulation environment, real-world marine applications would require additional preprocessing stages to handle increased noise levels and filter out smaller objects that are not of interest for specific mission objectives. These improvements will be considered in future studies.

3.1 Simulation Description

To validate the proposed method, we employed a ray tracing-based underwater acoustic simulator specifically designed for sonar simulations. The core operating principle of the simulator is based on ray tracing techniques [14], in which acoustic beams are represented as collections of rays. Object surfaces are modeled as triangular polygons; and the collisions between rays and polygons are detected. This simulator generates sonar data that closely resemble actual sonar images, thereby enabling a quantitative evaluation of algorithm performance. For FLS modeling, we utilized DIDSON specifications, whereas for PS modeling, 881A specifications were employed. The sensor specifications and experimental parameters used in our validation experiments are listed in Tables 1 and 2, respectively. Additional technical findings were reported in [15] and [16].

Table 1.

Specifications of forward Looking Sonar (DIDSON) used in experiment

Table 2.

Specifications of profiling sonar(881A) used in experiment

The mathematical implementation of the ray tracing information for the DIDSON and 881A specifications involves two key calculations: geometric intersection computation and acoustic intensity modeling. These key calculations enable accurate simulation of underwater acoustic propagation and reflection characteristics.

Additionally, to reproduce realistic sonar image characteristics, we applied speckle noise models that reflect the inherent acoustic properties of underwater environments [17,18]. The final intensity with added noise was calculated using a complex synthesis of coherent and incoherent echoes, as follows [18]:

Eq. (16) represents the final intensity calculation with speckle noise, where I_r is the measured intensity, which incorporates realistic noise characteristics. Through this comprehensive physical modeling approach, the simulator generated highly realistic sonar data that closely matched the actual underwater sonar characteristics, as shown in Fig. 7, which compares the real DIDSON data with our simulation results.

Fig. 7.

Real vs simulated FLS data comparison; (a) Actual and (b) Simulated DIDSON images.

3.2 Simulation Setup

In the simulation environment, the target object was placed on a seafloor with slightly irregular surface characteristics. The test object was a sloped box with 90-degree and 60-degree angled surfaces, as shown in Figs. 8 (d) and 9 (d), respectively, specifically designed for validating false slope suppression and reconstruction accuracy. This object geometry represents a common challenge in underwater sonar reconstruction, where the FLS vertical ambiguity typically causes false slope artifacts that degrade reconstruction quality.

Fig. 8.

Reconstruction results for 90-degree slope box; (a) FLS-based positive weights, (b) PS-based space-carving, (c) Fusion weight update, and (d) Scan direction and object geometry.

Fig. 9.

Reconstruction results for 60-degree slope box; (a) FLS-based positive weights, (b) PS-based space-carving, (c) Fusion weight update, and (d) Scan direction and object geometry.

The AUV was configured with FLS and PS sensors positioned at the same location. Scanning operations were conducted by approaching the sloped box from two opposite directions (0° and 180°) at a constant altitude of 3.12 m above the seabed. This dual-angle scanning approach enables comprehensive data collection from opposite directions, providing the multi-perspective information required for robust 3D reconstruction.

The simulation executed 1600 sequences with a 5 mm advancement per sequence. FLS data were recorded every 4th sequence (400 frames per scan), whereas PS data were acquired in every sequence (1600 measurements per scan).

3.3 Results

To validate the proposed method, we compared the reconstruction results of the FLS-only approach and the proposed FLS+PS fusion method using a slope box modeling file as the ground truth reference.

The reconstruction results for a slope box at 90-degree orientation are presented in Fig. 8. The positive weights estimated by the FLS (DIDSON) in Fig. 8 (a), demonstrate the vertical ambiguity inherent in the FLS measurements, with potential object locations distributed across the entire vertical FoV. Fig. 8 (b) displays the carving process on PS data, removing uncertain regions through space-carving operations, which effectively eliminates the false-positive areas marked in red. Fig. 8 (c) presents the final fusion reconstruction result, which successfully combines both sensor outputs, showing clear object boundaries and an accurate geometric representation. This result highlights the fundamental advantages of the occupancy-grid method. Instead of treating each voxel independently, our approach interprets the spatial relationship between voxels using the entire sensor-to-object path to carve away the empty space. This allows a single measurement to simultaneously inform on the state of multiple voxels, thereby maximizing information usage. Accumulating weighted evidence from multiple measurements makes the method robust against noise and false positives, which typically challenge point-based approaches.

The reconstruction results for a 60-degree slope orientation are displayed in Fig. 9. Although this configuration generates fewer false slope artifacts than the 90-degree case, inaccuracies still occur in the FLS-only reconstruction. The results demonstrate that the precise vertical information in PS allows the remaining false slope regions to be removed through accurate space-carving operations.

An important observation in both reconstruction results is the presence of shadow areas, indicated by the darker yellow regions in the object geometry diagrams (Figs. 8 (d) and 9 (d)). These shadow areas represent regions where objects physically exist but remain undetected during voxel reconstruction. This phenomenon occurs because of the directional nature of sonar propagation, in which acoustic beams cannot penetrate solid objects to gather information regarding the occluded areas behind them. This acoustic shadowing is an inherent limitation of sonar-based reconstruction systems and does not indicate failure of the proposed method. When these naturally occluded shadow regions are excluded from the evaluation, the reconstruction results accurately capture the visible object surfaces and boundaries.

Through a direct comparison between the ground truth and reconstructed voxel data, the results clearly demonstrate that the proposed FLS+PS fusion method successfully removes significant amounts of uncertain regions compared with the FLS-only reconstruction. The subsequent results show that our method addresses critical limitations observed in prior work. Unlike previous approaches that generated severe false slope artifacts in 90-degree box reconstructions and struggled to accurately estimate angles for inclined objects such as 60-degree slopes, the proposed method successfully eliminated these false slope problems and achieved an angle estimation accuracy similar to the ground truth across different object orientations.