Enhanced Multi-Beam Echo Sounder Simulation through Distance-Aided and Height-Aided Sound Ray Marching Algorithms

: The study proposes two innovative algorithms in the field of multi-beam echo sounder (MBES) simulation: distance-aided sound ray marching (DASRM) and height-aided sound ray marching (HASRM). These algorithms aim to enhance the efficiency and accuracy of MBES simulations, particularly when dealing with long-distance propagation and real-time processing limitations. DASRM addresses issues related to simulation accuracy by efficiently utilizing the KD-tree for spatial indexing and intersection detection instead of the signed distance field (SDF). Building upon the further analysis of DASRM, HASRM is proposed, which improves the search strategy for ray intersections and utilizes a height field pyramid for sampling and retrieval, thereby reducing memory usage while enhancing indexing efficiency. The experimental results demonstrate that both algorithms significantly outperform traditional methods in terms of simulation time, with HASRM exhibiting particular advantages in parallel computing due to its data structure and improved strategies. Additionally, DASRM is well suited for applications requiring complex scene construction, while HASRM proves especially effective in simulating MBES with a focus on underwater terrain due to its effectiveness in handling large incident angles and long-distance propagation.


Introduction
The ocean covers approximately 71% of the Earth and significantly influences the climate, the economy, and transportation.Thus, it is crucial for humans to correctly and efficiently understand, study, and develop the ocean.However, due to the different transmission media between the ocean and land regions, signals such as light and electromagnetic waves used on land are rapidly attenuated in water.As a result, underwater tasks typically require the use of acoustic signals, which have a wide range of applications, including underwater communication, location, distance measurement, speed measurement, and imaging.
Acoustic signals have a wide range of applications in underwater localization, acoustic imaging, and topographic surveys [1].In addition, experiments are necessary for testing and development.However, underwater experiments are limited by financial burdens, time, and manpower costs; unexpected accidents sometimes hinder the progress of projects.The significance of simulation experiments for cost reduction, risk reduction, and development efficiency is self-evident, and more realistic simulations can be applied to monitor system states and reproduce the environment [2].
Simulation is an essential component of the development and improvement of underwater acoustic devices, and algorithms with good simulation capabilities can also be applied to digital twin technology and virtual reality technology, allowing for the effective reduction of experimental costs, improving development efficiency, monitoring the working state of equipment, and simulating and analyzing equipment data.

Related Works
The primary objective of current research in simulating underwater acoustic equipment is to improve both accuracy and efficiency.The first type of investigation focuses on the physical properties of underwater acoustic signals, such as energy attenuation, reverberation, noise, and refraction by media.By employing computer graphics techniques to simulate underwater acoustic equipment, this line of research aims to improve simulation accuracy and provide more precise experimental methods for corresponding equipment studies.The second type of inquiry concentrates on optimizing algorithm efficiency through computer graphics techniques like rasterization, ray casting, BVH (bounding volume hierarchy), AABB (Axis-Aligned Bounding Box) [3], and ray tracing.This optimization enables a faster simulation of the physical properties of underwater acoustic signals and enhances the overall efficiency of simulating underwater acoustic equipment.Consequently, it establishes a foundation for real-time and efficient simulation algorithm research.The main objective pursued in this article is primarily focused on enhancing simulation efficiency.
Bell and Linnett [4] proposed a ray-tracing-based side-scan sonar simulation algorithm, incorporating acoustic signal propagation.They accounted for the curved propagation path of sound waves influenced by the sound velocity profile and employed a height field ray-casting technique.By employing larger steps outside the bounding box of the terrain height field and conducting fine searching inside it, they effectively reduced the computational load associated with ray tracing due to complex propagation paths.Guériot et al. [5] proposed a tube tracing technique that reduces the number of rays required by utilizing only four rays to achieve tube tracing.Additionally, they separated the acoustic signal propagation and rendering processes, enabling the simulation of multiple sonar devices using a single set of acoustic signal propagation histories.Coiras and Groen [6] investigated the side-looking sonar imaging process under the assumption of constant sound velocity and signal propagation along straight lines.They successfully demonstrated its application through simulation and three-dimensional sonar image reconstruction.Gu et al. [7] used ray-casting to observe the beam emitted by the sonar as a group of straight rays distributed according to the beam shape.By calculating the intersection points with the object's triangles, they simulated the image sonar.Kwak et al. [8] improved Gu et al.'s [7] method and introduced sound attenuation effects to generate grayscale sonar images.Saz et al. [9] introduced an acoustic model based on the assumption that the sound velocity is constant, combining ray tracing and frequency domain methods to generate high-quality sonar images, but with poor real-time performance.Aykin and Negahdaripour [10] used a ray-casting-based method to analyze quadratic surfaces and provide a better ray distribution pattern to address the problem of undersampling and oversampling caused by uniformly distributed rays, but mainly based on the assumption of ray propagation along straight lines.DeMarco et al. [11] combined the Gazebo simulator and ROS to design the FLS simulator [12].They used ray casting to generate point clouds and then converted them into sonar images, but did not consider the environmental influence on signal propagation paths.Gwon et al. [13] developed the missing SSS module in UWSim, using a simplified Lambertian diffusion model with rayleigh noise and speckle noise, but did not consider the influence of sound velocity on propagation paths.Cerqueira et al. [14,15] combined rasterization and ray tracing to optimize simulated sonar reflection, reducing the ray tracing area required through rasterization.They also used BVH and AABB algorithms to accelerate rendering [16], and employed osgOcean for the construction of underwater scenes, thereby achieving a real-time application of the simulator.However, only sound wave propagation along straight lines was considered.Ding, Rui, and Shiguang Liu et al. [17] proposed a novel method for simulating sound propagation in underwater scenes by incorporating an enhanced ray tracing technique that accounts for the unique characteristics of the underwater environment, along with a threshold-based approach to compute impulse response in the high-frequency domain.This approach efficiently calculates underwater sound propagation while yielding simulation results that closely align with real-world values.The validity of this method was demonstrated through various experiments conducted in underwater scenes, marking its pioneering contribution as the first tailored sound propagation model specifically designed for virtual reality applications within an aquatic setting.
Comparing the underwater sonar simulation with underwater acoustic communication simulation, it can be observed that due to the necessity of considering the actual propagation path of underwater sound rays, softwares such as WOSS [18] and COMSOL [19] are often used for underwater acoustic channel simulation.Methods like the Bellhop model [20] or finite element analysis are employed to simulate the propagation of underwater acoustic signals.However, since these simulations primarily focus on accuracy rather than real-time performance, they do not meet the requirements for real-time sonar simulations.
Only a few existing methods for real-time sonar simulation take into account the curved propagation of sound waves influenced by the velocity of sound.Among them, references [4,5,17] have proposed simulation methods for acoustic signals along complex curved paths; however, they all have certain limitations.In reference [4], a height field bounding box was utilized to reduce the computational cost of ray tracing, which required representing the scene solely with a height field and limited the types and complexity of simulated objects.This approach is more suitable for scenarios involving exclusively terrain height fields.Additionally, if the bounding box of the height field data in the scene is large, rays still need to perform numerous small step movements and corresponding detections, resulting in significant computational costs for ray tracing.Reference [5] adopted a method that separates rendering from acoustic signal propagation by initially recording the propagation history of acoustic signals and subsequently conducting rendering and simulation based on sonar equipment characteristics.However, this method primarily suits scenarios where simulation conditions are infrequently altered and only changes in sonar device model definitions occur; it is unsuitable for situations requiring modifications in acoustic signal propagation within scenes.Reference [17] employed a threshold method and utilized three distinct sound speed models, namely constant gradient within layers, constant sound speed within layers, and fixed sound speed.The selection of acoustic models with varying levels of precision was based on the variation in the incidence angle of the sound wave.A higher-precision acoustic model was employed when significant changes in angle occurred, while a lower-precision model or straight-line propagation was used for minor angle variations.However, certain assumptions regarding mid-high frequency and shallow water areas were still necessary to enhance efficiency.Moreover, the lack of addressing the fundamental issue of intersection detection for curved rays is evident in the approximations made for different acoustic models.Therefore, in order to enhance the practicality of simulation algorithms and accurately simulate scenarios involving long propagation distances (whether the distance is long or not depends on the actual distribution of sound velocity profiles, typically referring to situations where the actual path of sound rays deviates significantly from a straight line), particularly for MBES devices, it is crucial to explore new simulation methods.

Contribution
The paper addresses the simulation artifacts and real-time issues associated with sound ray simulation for the MBES over long propagation distances.Initially, widely-used underwater sonar simulation frameworks are presented.Based on the framework and leveraging the distance-aided ray marching algorithm, a novel approach called distanceaided sound ray marching (DASRM) is introduced.This method utilizes the KD-tree instead of the signed distance field (SDF) to mitigate issues related to simulation accuracy and real-time performance.Following this, the paper explores unresolved challenges within DASRM and introduces a height-aided sound ray marching (HASRM) algorithm.HASRM employs a height field pyramid for terrain profile sampling and spatial indexing, effectively reducing memory usage and enhancing indexing efficiency.By focusing on calculating target position depth instead of propagation distance, this algorithm minimizes required iterations, thus improving simulation efficiency.Ultimately, experiments are presented in this paper comparing traditional intersection test algorithms based on BVH and AABB with DASRM and HASRM, providing a critical analysis of their respective strengths and weaknesses while suggesting potential areas for future improvement.

Basic Sonar Simulation Framework
Underwater sonar simulation commonly employs ray theory, which precisely models the acoustic beam as an array of rays with various incident angles.This approach facilitates the accurate simulation of the sonar beam's measurement process.Initially, following standard sonar measurement principles, the local coordinate system (n-frame) is designated as the geographic coordinate system, while the body coordinate system (b-frame) is configured with respect to the vehicle and oriented forward-right-down.Based on these coordinate system definitions, the transformation matrix C n b is calculated to facilitate converting the ray direction u b and ray coordinates w b from the b-frame to n-frame within the local reference.
As shown in Figure 1, the basic simulation framework begins by generating multiple rays in the b-frame to simulate the acoustic beam, based on the MBES's beam shape, opening angle, and configured sonar settings.An error is introduced to account for sensor calibration, represented by σ sensor , as well as the installation position of the MBES at p mbe and its associated installation error denoted by σ install .Both σ sensor and σ install primarily affect the beam pointing error, which is essentially similar to the attitude error.Therefore, these two types of errors are combined with the attitude error to calculate the coordinate transformation matrix In the n-frame, the sound ray tracing method uses the actual sound velocity profile (SVP) svp true to calculate the actual trajectory of a sound ray.Then, by using intersection detection techniques such as ray tracing or ray casting, the intersection points w inter between the ray and the object being tested, along with the propagation time t prop and echo intensity I echo , are calculated.
Measurement outcomes are obtained by performing calculations based on simulated propagation time t prop and echo intensity I echo , using the measurement principles of the simulated sonar equipment.
MBES commonly employs the constant gradient sound ray tracing algorithm to obtain measurement results, which assumes an underwater environment with multiple layers characterized by varying sound velocity and a constant gradient, as shown in Figure 2. Within these layers, rays propagate along arcs with constant curvature determined by Equation (2).
In this equation, θ represents the incident angle of the sound ray element ds, α denotes its glancing angle, and c signifies the sound velocity at location ds.Equation (3) expresses the curvature radius of the arc.
Within the i-th layer, the horizontal distance increment ∆h is given by Equation (4).
Finally, the propagation time in the i-th layer can be calculated using Equation (5).
According to Equations ( 4) and ( 5), the time of propagation t i and trajectory (∆h, z i ) of the sound wave in each layer can be calculated.By subtracting t i step by step from the measured propagation time until reaching the N-th layer, the total propagation time of will be exhausted.Finally, based on the remaining propagation time in the N-th layer, we can obtain the measurement.
The simulation of sound rays is reversed compared to the process of the aforementioned sound ray tracing algorithms, which require obtaining the propagation time of sound rays based on their propagation trajectory for simulating sonar measurements.
Existing sonar simulation algorithms, which are commonly developed based on 3D rendering engines, assume the linear propagation of rays and utilize BVH and AABB algorithms to accelerate intersection detection.To accurately simulate sound ray propagation, it is necessary to employ ray tracing algorithms that follow curved paths instead of straight ones.This requires approximating a curved path with short straight segments, resulting in low efficiency and high computational costs.Therefore, developing a novel algorithm for MBES simulation is necessary in order to enhance the real-time performance and accuracy of the simulation.

Distance-Aided Sound Ray Marching
Distance-aided ray marching algorithms are extensively used for rendering various phenomena such as liquids, deformations, mixtures, and volumetric clouds [21,22].The algorithm iteratively moves for each ray based on the signed distance to the nearest object surface provided by the SDF [23], repeating the process until an object intersection is encountered.The distance can also be interpreted as the radius of a sphere centered at the current position that is free of other objects, which is why the method is also referred to as sphere tracing.This technique ensures that rays move safely and guarantees that no potential intersections with other objects are missed, as shown in Figure 3.The propagation paths of sound rays are curved due to the influence of the SVP.Constructing an SDF based on straight-line propagation paths leads to inaccuracies.Additionally, MBES simulation scenarios often involve large spatial scales, which require significant storage and construction time.Therefore, constructing an SDF for MBES simulation is not cost-effective.As a result, we replace the SDF with a KD-tree.The KD-tree is a data structure commonly used for spatial indexing in high-dimensional space searches [24], enabling the efficient search for the nearest point G near and the distance D near , thereby reducing overall computational and construction time.

……
The distance D near derived from the KD-tree, as depicted in Figure 4, does not precisely represent the actual distance from the current position to the underwater terrain due to the spatial resolution of the data, which introduces an associated error e r .However, considering that curved paths also introduce some error and e r is a relatively minor positive value, it is acceptable to initially ignore e r .Subsequently, by establishing appropriate iterative stopping criteria, the influence of e r can be minimized.The KD-tree is constructed using the terrain data, as depicted in Figure 4. Within the KD-tree, we can obtain the nearest point G near and its distance D near from the current position w i .Then, the sound ray moves forward along its current direction u i .Due to the inherent properties of curved paths, D near 's move along the ray ensures that the subsequent position w i+1 remains within a spherical boundary centered at G near with a radius of D near .
When calculating w i+1 based on a given propagation distance D near , we can first calculate its possible maximum change ∆θ in value of θ (∆θ = D near /R i ).If |∆θ| is smaller than or equal to the difference |θ i+1 − θ w i | in the θ values corresponding to the boundary of the sound velocity layer, it implies that w i+1 lies within the current sound velocity layer.Finally, the corresponding position changes in coordinates (∆h, ∆z) can then be calculated using Equation (6).
If ∆θ exceeds |θ i+1 − θ w i |, it suggests that w i+1 is located outside the current layer.Subsequently, we set ∆θ as ∆θ = θ i+1 − θ w i and propagate it to the boundary of the current sound velocity layer according to Equation (6).We update the remaining propagation distance as D near = D near − R i ∆θ, and continue iterating until the remaining distance D near equals 0. Simultaneously, we compute the propagation time t prop using Equation (5).
When the sound ray approaches the terrain, it is necessary to calculate the intersection point.Since the KD-tree only provides an unsigned distance, which is insufficient to determine if w i is inside the terrain, directly applying the ray marching algorithm for intersection point computation is not feasible.Therefore, it is necessary to establish an appropriate distance threshold D limit .When D near is less than or equal to D limit , we can ignore the curved nature of the sound ray and search for surrounding n triangles centered on G near .Then, we apply the Möller-Trumbore ray-triangle intersection algorithm to calculate the potential intersection point w inter [25].
As n increases, accuracy improves, but computational performance decreases.Therefore, when choosing an appropriate value for n, it is important to balance the desired spatial resolution with the available computational resources.Generally, a value of n greater than or equal to 4 is recommended (where 4 represents the nearest neighbor triangle and its associated triangles).Based on our experience, 6 is a preferable choice.However, if the computational device's performance allows it, using a larger value of n can also be advantageous in finding intersection points more quickly.For larger values of n, considering calculating the AABB for the respective triangles may improve the computational performance.
The value of D limit should be set carefully to strike a balance between algorithm efficiency and accuracy.Setting an excessively large D limit may result in unnecessary intersection tests, which can reduce performance.Conversely, setting an overly small D limit could lead to missed intersections, thereby affecting the accuracy.
When the proximity is sufficient, the ray can be approximated as a straight line that intersects with the triangle at the point w inter .In this scenario, w i+1 is located on an arc with the center at w i , a radius of D i , and a chord length of l i .The D i denotes the distance from point w i to the nearest vertex G near of the triangle.The distance from w i+1 to G near is denoted as D i+1 .Since the ray has intersected with the triangle, the chord length l i resides within the plane of the triangle.
Given that the other vertices of the triangle are not the nearest ones, the distances between them to w i exceed D i , hence verifying that l i is no greater than the maximum edge length l tri of the triangle.Furthermore, since w i+1 lies on the arc, it follows that D i+1 is no greater than l i .Consequently, an upper bound for D i+1 is established: D i+1 ≤ l tri .To prevent the omission of intersection points, the distance threshold D limit should be equal to l tri .

Height-Aided Sound Ray Marching
The DASRM algorithm can address the simulation artifacts and real-time processing constraints that arise from long propagation distances; however, it still faces several challenges in practical scenarios due to the unique characteristics of underwater environments.
Firstly, in the simulation of MBES, there is a need to model and perform sonar simulations for large underwater terrains.Although using a KD-tree can decrease memory usage and speed up construction compared to SDF, constructing it still involves significant computational expenses.Moreover, the data structure of the KD-tree is not suitable for parallel processing, which makes GPU acceleration impractical.As a result, the DASRM still faces significant computational cost issues.
In the DASRM, during intersection testing, it is still necessary to employ the traditional ray-triangle intersection method.This requirement necessitates the need for storing the corresponding triangle meshes.Given that terrain scenes in simulations are often large, the required storage space is increased, which significantly diminishes the practicality of the simulation.Consequently, there is a need to refine the relevant data structures and iterative strategies to better accommodate MBES simulation requirements, thereby reducing both spatial and temporal complexity while enabling parallel processing to enhance simulation speed.
According to the constant gradient sound ray tracing theory [26], which assumes that sound rays are only influenced by changes in sound velocity along the z axis, horizontal variations in sound velocity are small and do not affect the trajectory of sound rays.Therefore, the trajectory of the sound ray in the horizontal direction only involves changes in horizontal coordinates, while maintaining its original direction.This means that the intersection point's horizontal coordinate must lie on the projection ray within the xoy plane.As a result, there is no need to search through the entire terrain dataset; instead, one can simply sample the terrain profile at the corresponding horizontal coordinate along the projection line in order to minimize retrieval overhead.Furthermore, when accelerated by GPU, it becomes possible to obtain results directly with a single search operation without requiring multiple iterations.This approach offers advantages over the KD-tree as it saves storage space and accelerates retrieval rate.
In the context of computational simulations, it is noted that sound rays with larger θ angles require a greater number of iterations to converge near the actual intersection point if they continue to move according to D near as depicted in Figure 4 (the sound ray with a larger incident angle took four iterations to reach the vicinity of the intersection, while the vertically downward sound ray only took two iterations), which leads to a deceleration of the simulation.
Given that the depth at the intersection point w inter is denoted as z inter , and assuming that the sound ray propagates nearly in a straight line along this segment, the distance D inter between w inter and w is z inter −z w cos θ .Equation ( 7) expresses the distance error ε w between D inter and D near .
Meanwhile, according to Equation ( 8), the terrain depth z is expressed as the sum of its mean z and random noise ν, which follows a normal distribution.
Underwater terrains are predominantly flat, and these data typically have low resolution.Therefore, when the sound ray is at a considerable distance from the seafloor, the nearest point G near is generally located just below or in close proximity to the w.In such cases, D near can be approximated as z near − z w .By substituting Equation (8) into Equation (7), ε w can be determined by Equation (9).
The error in Equation ( 9) is mainly influenced by three factors: • The first factor, 1 cos θ − 1, is dependent on θ.The larger the value of θ, the greater the value of ε w becomes.When θ is zero, the sound ray propagates vertically downward and G near and w inter coincide, resulting in a zero value for ε w ; • The second factor, represented by z − z w , indicates the distance from w to the seafloor.The greater this distance, the higher the value of ε w .When the sound ray approaches close to the seafloor, this error can be disregarded; • The third factor is random noise ν, mainly generated by the complex variations of the terrain, which primarily affects intersection testing.Algorithmic optimization is required to prevent the sound rays from missing the actual intersection positions.
The primary limitation for the second part of ε w lies in the distance between MBES and the seafloor.In order to minimize ε w and reduce the number of iterations, it is essential to concentrate on mitigating the influence of the first part on ε w .
According to Equation (10), it is evident that the coefficient preceding z near should be equivalent to that of z inter , which is 1 cos θ .As a result, the first factor of ε w will naturally decrease to zero, and the constant error in the second factor, z − z w , will also be eliminated.This approach ensures that only the stochastic noise error component remains within ε w .
When the propagation distance is set to z near −z w cos θ , this essentially utilizes z near as the target position depth z aim and calculates the corresponding horizontal coordinate at the depth of z aim , as given in Equation (11).If the maximum depth of the current sound velocity layer is less than z aim , then we propagate the sound ray to the boundary of this layer.We repeat this procedure until z w = z aim .
However, the sound ray may be affected by its curvature and the effects of ε w , which could result in premature intersections with the terrain and computational inaccuracies.Therefore, it is crucial to define a suitable target depth z aim instead of using z near directly, in order to minimize the influence of residual random noise in ε w and account for the curvature of the sound ray.
According to Snell's law, when the incident angle (θ) is 0 • , the sound ray enters the medium vertically without being affected by the velocity of sound and propagates downward in a vertical manner.As θ approaches 90 • , there is a tendency for the horizontal propagation of the sound ray.If there are further changes in sound velocity causing bending, gradual propagation towards decreasing depths occurs.Therefore, achieving the theoretical sweeping width corresponding to its opening angle is challenging in MBES due to both bending and attenuation effects on sound rays.As a result, it can be considered that incident angles of sound rays always fall within 0 • to 90 • .Consequently, the depth z increases monotonically with the horizontal distance h, which can be represented as the continuous function z w = F(h w ) for the trajectory of the sound ray.Similarly, the terrain profile, derived from sampling the terrain height field, can also be described by a continuous function d = G(h).
Both functions, F and G, are continuous.The intersection point between the sound ray trajectory and the terrain profile corresponds to the first zero value of the continuous function P(h) = F(h) − G(h), as indicated in Equation (12).
However, both the terrain profile and sound ray trajectory are complex curves, and cannot be precisely represented by functional expressions.Consequently, directly calculating the intersection points through zero values is not feasible.It is essential to conduct a detailed analysis of the intersection between the terrain profile and sound ray trajectory to narrow down the potential search range of intersection points.Then, the parametric approach can be utilized to approximate the functions and obtain the intersection points.
The terrain profile data are derived from sampling the height field of the terrain.During the sampling process, various interpolation algorithms are typically applied, such as linear interpolation, bilinear interpolation, or cubic interpolation.Each of these interpolation algorithms follows a specific set of mathematical rules.Generally, the height value at an interpolated point is a linear combination of the height values of adjacent sampling points.As a result, the terrain profile is usually continuous and differentiable.
Therefore, as shown in Figure 5, the intersection between the terrain profile and sound ray can be broadly classified into three scenarios:

•
G ′ (h) = 0, which is precisely located at the position of a local minimum point G locmin ; • G ′ (h) < 0, indicating a decreasing trend in terrain depth, which is opposite in sign to the derivative of sound propagation F ′ (h); • G ′ (h) > 0, indicating an upward trend in the depth of the terrain, while When G ′ (h) < 0, the depth of the terrain decreases as h increases.However, it is not possible for the underwater terrain depth to decrease indefinitely; therefore, there will inevitably be a local minimum G locmin .In this case, once we accurately identify G locmin , we can find w inter within a smaller range of horizontal coordinates.
When G ′ (h) > 0, the depth of the terrain increases with an increase in h.However, it is not possible for the depth of the terrain to infinitely increase.Therefore, the depth of the terrain will inevitably reach a maximum value and then decrease to G locmin .
Similarly, when G ′ (h) = 0, regardless of whether the intersection is located at a local maximum or minimum point, it must lie within the interior region of a G locmin .However, due to limitations in our search range, there may be situations where the terrain profile within this range changes monotonically without containing any extrema.Therefore, we should also consider values at the boundaries of our search range.As long as these values are smaller than those on their inner side, they can also qualify as local minima.
Therefore, the problem can be efficiently detected by dividing it into two phases.Firstly, find the nearest local minimum value G locmin of w inter , which helps narrow down the search range to (h w , h G min ).Subsequently, we conduct an accurate search to obtain the w inter .
Before establishing the search range (h w , h G min ), we can initially set the search range as (h w , n(G(h w ) − z w )), where n is determined by the actual performance of the MBES's sweeping width.This ensures that the subsequent position w i+1 remains within the search range.
According to Equation ( 13), the target depth z aim should be established as the shallowest depth G min within the local minima G locmin of the function G(h) within the current search range.It is important to note that any minima with depths less than or equal to z w are excluded from consideration.This avoids scenarios where the sound ray is unable to iteratively move towards the intersection point after its position has been adjusted according to the local minimum.These scenarios, as shown in Figure 6, include instances where h w equals, exceeds, or falls short of h G min .Figure 6.The sound ray movement strategy.Three sound rays initiate their search from (h w0 , z w0 ).After the first search, z aim for these sound rays is determined to be the local minimum G 1 locmin .Sound ray 1 then moves directly to the intersection point.For sound rays 2 and 3, since their horizontal positions satisfy , it is clear that G 1 locmin is not the target minimum for them, and they require further movement.After setting their target depth to a new local minimum G 2 locmin and executing the movement, it is observed that sound ray 2 intersects with the interval where h 2 w2 is less than the horizontal position of G 2 locmin .This successfully obtains the precise search range for the w 2 inter .
When h w equals h G min , the sound ray intersects with the terrain.If h w exceeds h G min , the local minimum G min is considered to be outside of the current search range and should be excluded from further consideration.Therefore, only scenarios where h w is less than h G min need to be analyzed.
In this scenario, since the movement of the sound ray is guided by the local minimum G min , it will not intersect with the terrain before reaching the intersection point.Therefore, it is confirmed that the intersection point h inter lies within the interval (h w , h G min ).However, G min may not be the closest local minimum to the actual intersection point.To refine the search range and find a closer minimum to w inter , we should move the sound ray towards deeper local minimum values within this range.
As Equation ( 13) illustrates, when calculating the movement distance, only local minimum values G locmin with depths exceeding the current sound ray depth z w should be taken into account because shallower minimum values have already been evaluated.When there are no more G locmin to consider, it indicates that the sound ray has completed its initial search for w inter .
Once the search range for the intersection point is established, to streamline subsequent computational processes, the maximum terrain depth G max is designated as the maximum depth z max , and the current sound ray depth z w is designated as the minimum depth z min .This depth interval (z min , z max ) is then used to represent the search range.If this range is sufficiently narrow, a parametric approach can be applied to compute the location of the intersection point.
In cases where the search range is larger, a threshold z lim for the search range can be established.When the search range exceeds this threshold, the dichotomy method can be used to further narrow down the range.The target depth z aim is set as the midpoint between the current z max and z min , (z max + z min )/2.After relocating the sound ray, adjustments are made to either z max or z min based on the difference between z w and the terrain depth G(h w ), as well as between h w and h G min , effectively reducing the search range by half in each iteration.
The threshold z lim is set to ensure the accuracy of the parametric approach in determining the search range.The accuracy of calculating intersection points using parametric approaches mainly depends on the precision of approximating the difference P(h) between terrain profile curve and sound ray curve.The simpler the P(h) curve, naturally, the higher the accuracy of parametric approximation.Generally speaking, since the search distance near intersection points is not too large, sound ray curvature has a relatively small impact on P(h).The complexity of P(h) is primarily influenced by variations in the terrain.
Therefore, setting thresholds is related to the terrain resolution and interpolation method.For example, in this paper, 10 m resolution terrain data were used for simulation with bilinear interpolation for terrain sampling.Therefore, the horizontal range threshold h lim for the search range would be 10 m.Bilinear interpolation uses four points for interpolation, and the distance between the nearest four terrain grid points around one terrain grid point is 10 m.Thus, when the horizontal range limit is set to 10 bilinear interpolation can ensure that the sampled values of the terrain profile within conform to a relatively simple curve.If other interpolation methods are used, h lim can be redefined based on the window size and calculation method of the interpolation method.Additionally, the depth range threshold z lim would be calculated as h lim /tan(θ).
As shown in Equation ( 14), when the size of the search range is relatively small, we can employ a parametric approach to handle the depth difference function P(h) between sound rays and the terrain.Specifically, we can approximate this function, compute relevant coefficients through sampling, and ultimately calculate the zero positions of P(h) to obtain intersection points.This approach takes advantage of the fact that the search range is already sufficiently small, ensuring both accuracy in approximation and efficiency in computation.
In Equation ( 14), the coefficients k 1 , k 2 , . . ., b have yet to be determined and need to be solved by sampling.Due to the gentle nature of underwater terrain and the low curvature characteristics of the sound ray, a first-or second-order Taylor expansion is sufficient for obtaining accurate estimates of intersection positions.On the other hand, higher-order expansions are more complex when solving for zero values and are not suitable for realtime computation.
By following a three-step strategy that begins with searching for the target location of the minimum value, narrowing down the search range using the dichotomy method, and concluding with precise intersection calculations through the parametric approach, we ensure both precision and real-time capabilities in the simulation.This approach significantly reduces storage requirements by utilizing only height field data, eliminating the need for additional KD-tree and mesh data storage.
However, when conducting terrain profile sampling based on the height field, a coarse-to-fine strategy is also necessary.This requirement arises from the need to synchronize the sampling of multiple terrain profiles when using parallelized computation techniques.To maintain this synchronization, a uniform number of sampling points is required.Consequently, for a height field with a set resolution, a large search range may result in undersampling and potentially omitting the minimum value points of interest, thus affecting the accuracy.
Therefore, the use of a pyramid data structure is essential.The terrain pyramid is constructed by repeatedly applying a 2 × 2 window and downsampling the minimum depth value of the terrain with a stride of 2 until a pyramid with n levels of progressively lower resolution terrain data is established.resolution terrain data are utilized for sampling in large search ranges, while higher resolution terrain data are preferred for relatively small search ranges.
As depicted in Figure 7, the dashed line represents the terrain profile achieved through minimum downsampling.While this technique reduces the terrain's resolution, making it appear flatter, all minimum values are maintained, preventing premature intersections.When the search range is narrow enough, employing high-resolution terrain data in conjunction with interpolation algorithms enhances the data's resolution and precision, ensuring accurate intersection calculations and reducing the need for pre-interpolation.locmin within the low-resolution terrain profile.Once the target depth is determined, the sound ray is directed to that depth.Then, utilizing the high-resolution data, it identifies the subsequent local minimum G 2 locmin within a constrained search range and relocates to w 2 .Since there are no additional local minima, the interval (h w2 , h G min ) is established to finalize the first phase.In the second phase, the depth range (z min , z max ) is defined by a horizontal range, and through the dichotomy method, the sound ray moves to positions w 3 , w 4 , and w 5 , thereby reducing the search range.Finally, by employing a parametric approach and using data samples from w 3 , w 4 , and w 5 , the difference function P(h) is computed to pinpoint the intersection point (w inter ).
The basic principle of the HASRM algorithm is illustrated in Figures 7 and 8, which combines all the mentioned intersection calculation strategies and terrain profile sampling methods.

Results
The real-time performance of the proposed algorithms was validated by employing a traditional ray-casting algorithm based on BVH and AABB as the comparative algorithm.The rays propagate according to constant gradient sound ray tracing and advance in fixed increments of 20 m at each step.If the AABB distance exceeds D limit , the ray will move to the next position without further AABB testing or ray-triangle intersection testing.Detection is considered successful when the distance between w inter and w i is smaller than D limit .
To ensure a fair assessment, the BVH, DASRM, and HASRM were implemented in Python 3.9 and executed serially on a CPU.The experiments were conducted on a platform with a 32 GB memory and an Intel(R) Core(TM) i9-13980HX processor.
The study employed three simulated underwater terrains with average depths of 100 m, 500 m, and 1000 respectively.The simulation utilized a beam configuration of four rays, with 128 beams oriented at an opening angle of 1 • × 1 • , while the MBES had an opening angle of 120 • .
The MBES simulation was conducted with 100 trials for each of these algorithms, and the average simulation times are presented in Table 1.The results presented in Table 1 indicate that DASRM and HASRM have significant advantages over the traditional method in terms of average simulation time.Furthermore, as the average depth increases, the advantage becomes more pronounced.This is because the traditional algorithm based on fixed distance movement requires more iterations to reach the intersection point with increasing depth.In contrast, DASRM or HASRM can calculate more reasonable propagation distances based on the actual terrain distribution, enabling faster arrival at the intersection point.Due to its superior movement strategy and faster data retrieval speed, HASRM has a greater advantage.
Furthermore, HASRM is particularly well suited for leveraging the parallel computing capabilities of GPUs, thereby significantly enhancing simulation speed.We employed the PyTorch 1.13 framework to effectively parallelize computations involved in both DASRM and HASRM (excluding the KD-tree).By employing an RTX 4070 laptop GPU, we validated the performance of DASRM and HASRM and compared these outcomes with those obtained from their CPU-based implementations.
By configuring the MBES parameters to include 128 beams with the MBES opening angle of 120 • and a beam opening angle of 1 • × 1 • , We conducted simulations on beams with varying ray quantities, performing 100 trials for each quantity and averaging the simulation times on a with an average depth of 1000 m.The obtained results are presented in Table 2.The simulation results presented in Table 2 demonstrate that, for both algorithms, the GPU versions outperform the CPU versions in terms of speed when the number of rays is higher.This improvement is due to the GPU versions' increased utilization of parallel processing capabilities as the number of rays grows.However, the performance enhancement achievable with GPUs is restricted in DASRM due to its reliance on the KD-tree for terrain data retrieval through the CPU before processing by the GPU, which sets it apart from HASRM.
Despite its relatively limited parallel processing capabilities, the CPU still maintains an advantage in executing complex computations.As a result, when the number of rays is low, the CPU versions outperform their GPU counterparts in terms of speed.
The proposed algorithms can simulate long-range detection of MBES thanks to improved algorithm strategies.Assuming 15 beams with an opening angle of 120 • and a flat ground surface, the simulation effect is shown in Figure 9.The red solid line represents the central path of the sound beam, light blue represents the edge of the beam, green dashed lines represent rays along each beam's initial straight direction at its starting point, black dots and solid lines represent actual intersections between sound beams and underwater terrain, while pink dots represent measured terrain obtained by a sound beam tracking algorithm based on SVP calculation.The SVP used in the simulation is shown on the left side.
The trajectory of the beam deviates from its initial straight path in response to variations in SVP, as shown in Figure 9. Beyond a depth of 300 m, the path starts bending upwards and then transitions into a downward curve at depths exceeding 600 m, ultimately resulting in straight-line propagation.The actual measurement value is generated by the intersection between the beam edge and the terrain, which is influenced by the conical shape and specific opening angle of the acoustic beam in MBES.Consequently, the depth calculated along the central ray tends to be slightly shallower than the actual depth, with a subtle upward bend observed at the edge of computed terrain.In MBES measurements, inaccuracies in the SVP can result in increased measurement errors.HASRM effectively simulates this phenomenon.Figure 10 illustrates calculations using SVPs with varying gradients: (a-c).In Figure 10a, the employed SVP accurately matches the simulated one, thereby attributing errors primarily to variations in beam opening angle that progressively increase with propagation distance.This error pattern is analogous to that observed in Figure 9.In Figure 10b, the utilization of a higher gradient in sound velocity and average sound velocity leads to increased curvature along the propagation path, resulting in computed terrain depths that exceed the actual depth with an upward curvature.Conversely, Figure 10c exhibits a distinctly contrasting scenario.We maintained consistent simulation conditions and conducted terrain measurements across randomly generated terrain, characterized by an average depth of 500 m and a spatial resolution of 10 m × 10 m.We utilized the SVP illustrated in Figure 10 for simulation purposes, and computed the measurement outcomes for each sound velocity profile presented in Figure 10a-c.The obtained results are depicted in Figure 11.The observed trend of measurement errors displayed in the figure aligns with that demonstrated in Figure 10.

Discussion
The experimental results demonstrate that the DASRM and HASRM algorithms introduced in this study outperform traditional methods in terms of real-time performance.They are also capable of accurately simulating the complex propagation paths of acoustic signals, enabling precise simulation of MBES measurement errors.The DASRM algorithm, an adaptive improvement based on sphere tracing, enhances the iterative efficiency of finding intersection points.In addition to the DASRM algorithm, the HASRM algorithm further focuses on underwater terrain measurement scenarios of MBES.It analyzes factors that affect the accuracy of sound ray movement in iteration and optimizes the sound ray simulation process by adopting a progressive intersection detection strategy-from searching for minimum value positions to narrowing down search ranges using dichotomy method, and finally, calculating intersections using parametric approach.This algorithm uses target depth to move sound rays, replacing the nearest neighbor distance in DASRM, thereby reducing the negative impact of large incident angles on simulation speed.
The adoption of the HASRM algorithm ensures a relatively consistent number of iterations for each sound ray, regardless of the incident angle or propagation distance, as shown in Figure 12.Instead, it is solely influenced by the actual distribution of terrain profiles.Consequently, intersection detection can be completed more quickly than with DASRM.
In the context of data structures used for scene construction, DASRM incorporates the KD-tree as a replacement for SDF and maintains the mesh data of the scene for intersection testing.Conversely, the HASRM algorithm utilizes height field pyramids for retrieving and sampling terrain data, effectively replacing the KD-tree and mesh data.This approach not only further reduces memory usage but also enhances retrieval speed, resulting in a significant improvement in parallel computation performance.However, since the HASRM algorithm relies on height field data for scene construction, it is unable to accommodate relatively complex objects due to the inherent limitations of height fields in representing complex geometries.Nevertheless, this restriction does not significantly impair the ef-fectiveness of the simulation algorithm given that MBES devices are primarily used for underwater terrain measuring.Therefore, when considering the application of the algorithms presented in this paper to other sonar systems or more complex scenarios, DASRM might be a more suitable choice.This is due to DASRM's implementation of data structures, which offer better scene description capabilities.To expand the scope of HASRM, one could explore the construction of multidimensional height fields across various dimensions such as x, y, and z, or consider adopting more sophisticated data structures to augment HASRM's capacity for simulating complex scenes.
Furthermore, there is potential for further refinement and advancement in the algorithms introduced in this paper.For instance, incorporating the spectral characteristics of acoustic signals could result in a more precise simulation of the underwater acoustic field, thereby improving the applicability of these algorithms.Additionally, in order to enhance the precision of the simulation, it is crucial to consider the influence of horizontal velocities on sound rays and revise HASRM's terrain profile sampling strategy, as the computation of sound ray paths currently assumes that rays are only affected by the vertical component of SVP.This may necessitate an independent calculation of the projected trajectory of sound rays within the xoy plane.

Figure 3 .
Figure 3. Distance-aided ray marching.The solid black lines in the figure represent objects in the scene, while the dashed black circles represent non-intersecting spherical regions corresponding to distances given by SDF.The red dots represent the positions after each step forward along the ray, and the black dots represent points on the objects closest to the red dots.The blue dots indicate the final intersection point location.

Figure 5 .
Figure 5.The situations of intersection points.

Figure 7 .
Figure7.HASRM.The sound ray originates at w 0 and searches for G 1 locmin within the low-resolution terrain profile.Once the target depth is determined, the sound ray is directed to that depth.Then, utilizing the high-resolution data, it identifies the subsequent local minimum G 2 locmin within a constrained search range and relocates to w 2 .Since there are no additional local minima, the interval (h w2 , h G min ) is established to finalize the first phase.In the second phase, the depth range (z min , z max ) is defined by a horizontal range, and through the dichotomy method, the sound ray moves to positions w 3 , w 4 , and w 5 , thereby reducing the search range.Finally, by employing a parametric approach and using data samples from w 3 , w 4 , and w 5 , the difference function P(h) is computed to pinpoint the intersection point (w inter ).

Figure 9 .
Figure 9. Simulation of sound ray trajectory by HASRM.Due to the change in SVP, the sound ray will undergo distortion.Utilizing HASRM can effectively simulate this bending of the sound ray.

Figure 10 .
Figure 10.Simulated measurement errors corresponding to different SVPs.(a) Simulated measurement error of zero-error SVP.(b) Simulated measurement error of greater SVP.(c) Simulated measurement error of smaller SVP.The utilization of HASRM and DASRM enables precise simulation of acoustic signal propagation time, facilitating an accurate sound ray tracing algorithm for pinpointing measurement positions.This approach offers a more realistic simulation of MBES measurement errors rather than imposing a fixed error term based on the system's measurement principles.

Figure 11 .
Figure 11.Comparison of simulated measurement errors for different SVPs in MBES.(a) Real terrain.(b) Simulated measurement error of zero-error SVP.(c) Simulated measurement error of greater SVP.(d) Simulated measurement error of smaller SVP.

Figure 12 .
Figure 12.The iteration number of sound rays with different θ in HASRM.

Table 1 .
Average simulation times of different algorithms.

Table 2 .
Average simulation time.