PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Study on Simulation of Steady-state Sound Field Matrix in Complex
Undersea Environment and Detection sensitivity
Jun SUO 1; Yangyang ZHANG 2; Simi TANG 3
Naval Research Academy, China
ABSTRACT
It's crucial to simulate the acoustic detection effect timely and accurately between any two positions in thecomplex undersea environment in a real-time simulation platform. A numerical analysis method of sound ray
propagating in complex undersea environment based on the eikonal equation, non-coplanar sound ray andsound velocity gradient vector is presented in this paper. The total undersea concerned information undercertain accuracy requirement can be obtained in the algorithm. The multidimensional sound pressure array isestablished to simulate the whole sound space. The calculated amount for the real-time simulation of theacoustic detection effect is decreased from O (mnu) by the traditional eigenray solution to O (1) (m is thenumber of the non-layered structures, n is the grid quantity of the local layered structure, u is the number ofthe iterations in Newton method or Runger Kutta method). The real-time capability is highly improved by
the method mentioned above. By comparison, the results of the simulation and the marine measurement are
in good agreement. The concept of detection sensitivity is proposed to quantitatively describe the relationshipamong marine environment, detectors and sound sources. The concept can give good expression of therelationship among the three factors
Keywords: simulation algorithm, non-coplanar, sound ray, sound velocity gradient vector, steady-state soundfield matrix, detection sensitivity

1. INTRODUCTION
The real-time simulation for the detection in seawater becomes more widely applied in different
fields (1-3). It is the key technology to obtain the interaction among the units in time and accurately
under a complex ocean environment in the simulation platform. Underwater acoustic detection
simulation plays an important role in different physical fields. It is crucial to calculate the underwater
acoustic attenuation rapidly among positions in the ocean model on the platform.
The basic method to calculate the acoustic propagation under water is by wave equation. All othermethods are derived from the wave equation method and make simplification and approximation insome aspects. The wave equation method has a high precision compared with other methods, whilemuch more calculation time is paid. The slow computation speed limits the wave equation to be used
in a real-time simulation platform (4).
The ray-based method provides higher calculation speed. The traditional "shooting method"acquires huge computational cost to get the eigenray. More error accumulated in the long range showsthe "shooting method" is not an ideal method. The ray-span algorithm is accepted with higher speed
to calculate the eigenray (5). While it is just suitable for the layered ocean, but not for the complex
undersea environment. In most of the researches, the sound ray propagates in a certain plane with afixed bearing angle. Such situation is far from the real complex ocean environment (6, 7).
A numerical analysis method of sound ray propagation with non-coplanar sound velocity gradientvector in complex undersea environment is presented in this paper. A steady-state sound field isestablished by the ray-tracing method without the conception of eigenray. The concerned undersea
information to meet certain accuracy requirement and the sound relationship between any two points
can be obtained by the algorithm. Compared between the test and the calculation, the results coincidewith each other well. The validity of the algorithm is approved.
The concept of detection sensitivity is proposed to quantitatively describe the relationship1

mrsuo@163.com
1109402164@qq.com
3
oldcandy@163.com
23221

among marine environment, detectors and sound sources. The concept can give good expressionof the relationship among the three factors by comparing the results of the real-time simulation for
operation units in marine environment and the practical measurement.

2. NUMERICAL ANALYSIS METHOD OF SOUND RAY PROPAGATION IN
COMPLEX UNDERSEA ENVIRONMENT BASED ON NON-COPLANAR SOUNDRAY AND SOUND VELOCITY GRADIENT VECTOR
The medium is non-layered structure, and the seafloor is not flat in a complex undersea
environment. The non-layered structure leads to the sound velocity gradient vector non-coplanar withthe vertical axis. The fact of non-coplanar should be considered in the derivation process of the eikonal
equation.
Non-coplanar sound ray and sound velocity gradient vector are shown in Figure 1.

)&g

T \M &s

H
F

Figure 1 – Sketch of the non-coplanar sound ray and sound velocity gradient vector)&&s is the sound ray vector. H is the plane that the sound ray lies in. g is the sound velocity gradientvector which is non-coplanar with the sound ray vector.
The eikonal equation deduced from Helmholtz equation is
2

W =

1
c 2 ( z, r )(1)z is the depth. r is the horizontal distance. c(z, r) is the sound velocity at (z, r). W is the travel time.
W is perpendicular to the wave front. Then the sound track equation can be defined as adifferential equation,dr=c( s )Wds

(2)

s is the arc length of the sound ray
Since the sound velocity gradient vector is not in the plane of the depth and the horizontal distance,
the sound track equation in the plane of the depth and the horizontal distance isd § 1 dr ·1c ( z , r )¨¸=  2c ( z, r )ds © c( z , r ) ds ¹

(3)

The quantity of c( z, r ) is different from that of the sound velocity gradient here because the soundvelocity gradient vector is non-coplanar with the sound ray vector. Replaced by the sound velocity
gradient, Equation (3) in the sound ray plane can be expressed as follows)&
d § 1 dr ·1
g cos T cos\ tan M  1¨¸=  2ds © c( z , r ) ds ¹c ( z, r )(4))&
g is the sound velocity gradient vector. T is the angle between the sound velocity gradient vector

and the sound ray propagation plane. M is the angle between the sound ray vector and the verticalaxis. \ is the angle between the vertical axis and the projection of the sound velocity gradient on theplane of the sound ray.
The travel time of the sound ray track is deduced from the eikonal equation.

3222

W ( s ) W (0)  ³

s

0

1
ds '
c( s ')(5)

The sound ray track and the sound travel time can be deduced by solving Equation (4), (5) throughthe numerical analysis method.
Thus, the sound ray track and the sound travel time at any transmission angle can be deduced by
the method mentioned above.

3. SOUND RAY TRACK IN STEADY-STATE SOUND FIELD
3.1 Sound field matrix
A receiving ball should be established in the traditional sound ray track method to acquire certain
point's information by the eigenrays. Simulation error is influenced by the radius of the receiving ball.
The technology of sound ray track in steady-state sound field simulates the sound field distribution
without the receiving ball.
One sound ray represents a part of the total power of the sound source
P(T , M )P
) (T , M )N(6)P is the total power of the sound source. N is the number of the sound rays. T , M are the glancing
angle and the azimuthal angle of the sound rays, respectively. )(T , M ) is the directivity index of thesound rays.
A set of 3D energy array can be established for a determinate sound space. The energy distribution
can be imitated by a 3D array at certain accuracy requirement. The sound pressure level normalizedsound source's influence on other units can be calculated in the matrix. The transfer functions and thetravel time among different units are stored in the energy matrix. The factors independent of
frequencies are stored, such as the seafloor reflection angle, sediments, etc. The further information(reflection attenuation, absorption attenuation, etc.) is calculated later by considering the soundfrequency. In real-time simulation, many data can be read from the matrix directly with little
computation.
The key point of the algorithm above is to rapidly calculate the sound intensity attenuation on acertain sound track.
3.2 Calculation of travel time
Literature (8) gives the eikonal Equation (1)'s expression approximated at high frequency incomplex medium as below,
2

2

2

1
§ wt · § wt · § wt ·¨ ¸ +¨ ¸ +¨ ¸ = 2© wx ¹ © wy ¹ © wz ¹ c ( x, y, z )(7)t is the wavefront travel time. c(x, y, z) is the sound velocity at (x, y, z). Discretizing the 3D
complex ocean medium and using the second order difference equation
wt 3ti , j ,k  4ti 1, j ,k  ti  2, j ,k
=
wx
2'x

(8)

Equation (7) can be approximately expressed as
3ti , j ,k  4ti 1, j ,k  ti  2, j ,k2

+ 3ti , j ,k  4ti , j 1,k  ti , j  2,k2+ 3ti , j ,k  4ti , j , k 1  ti , j , k  2

2§ 2'x ·¨¨¸¸
© ci , j ,k ¹

2

(9)i, j, k are the discretized grid numbers along x, y, z axis, respectively. The side length of discretized
cubic grid is 'x . If ti 1, j ,k , ti  2, j ,k , ti , j 1,k , ti , j  2,k , ti , j ,k 1 , ti , j ,k  2 are known, the unknown wavefront travel
time ti , j ,k can be obtained from Equation (9).
The wavefront travel time at the adjacent nodes can be calculated by the discretized eikonal
equation as Equation (9) from the sound source through GMM (Group Marching Method) (9).
The connection line between the adjacent nodes in the frontwave narrowband S tends to beperpendicular to the normal of the frontwave if the difference between the frontwave travel time of
the adjacent nodes is smaller than below,

3223

G ts

'x
3vSˈmax(10)vSˈmax is the maximum sound velocity of the nodes in the frontwave narrowband. The secondarysound sources in the frontwave narrowband areR^ i, j , k  S : t i, j , k d t

S ,min

 G ts `(11)

tS ,min is the minimum frontwave travel time in S.

Under the condition to keep the causality of the travel time, several secondary sound sources are
determined to extend the frontwave. The travel time of all the nodes can be deduced more efficiently.
The sound attenuation between any points in the ocean model can be calculated by the travel time.
3.3 Calculation of sound velocity
The sound velocity of each node is the basic information to calculate the travel time. The data ofsound velocity is usually scarce in the ocean test region, which is hard to describe each node's
information. It is necessary to calculate each node's sound velocity according to the measured data tofulfil the deduction in Section 3.2.
s x, y, z is a point in the ocean model. Four points named s0 x0 , y0 , z0 , s1 x1 , y1 , z1 , s2 x2 , y2 , z2 , s3 x3 , y3 , z3&&&
are selected which are the nearest points to s. If c s3 -s0 is coplanar with a s1 -s0 䚸 b s2 -s0 , anotherpoint is reselected as s3 to make the three vectors non-coplanar.
p, q, k are the undetermined coefficients. According to the relationship of perpendicular vectors,
certain p, q, k can be solved to satisfy& & )& &­ a u b  d  kc
° & & )& &°® c u a  d  qb
° & & )&&
° b u c  d  pa
¯
)&
d

0

(12)

0
0

s  s0 .

So& & )&
buc d& & & ,qbuc ap& & )&cua d& & & ,kcu a b& & )&
aub d& & &aub c(13)

The relationship of the coefficient k and other elements is shown as below.
s3

& &
aub

c

)& &
d  kc

&

s

)&

kc&

d

s2
&
b

s0
&
a

s1

FigXUH 2 - Relationship of sound velocity data points, undetermined coefficient and point to besolvedv0 䚸 v1 䚸 v2 䚸 v3 are the sound velocities at s0 䚸 s1 䚸 s2 䚸 s3 , respectively. The vector of the sound)&velocity difference vs between s and s0 䚸 s1 䚸 s2 䚸 s3 is
)&
vs

&
&
&
pa
qb
kc
& v1  v0  & v2  v0  & v3  v0abcThen, the sound velocity at s is

3224

(14)

vs

)& )&
vs  dv0  )&d(15)

Through the above algorithm, the sound velocity at any point in the ocean can be calculated fromthe given data.

4. REAL-TIME SIMULATION IN COMPLEX UNDERSEA ENVIRONMENT
10000m×10000m×4000m ocean model is established based on the 3D sound velocity profile tested
in the seawater. The seafloor is half-space and liquid model. The basic parameters of the sediments
are shown as below.
Table 1 - Basic parameters of sediments
Type
Coarse sand
Silt
Slime
Density/g·cm-3
2.034
1.806
1.469
Sound velocity/m·s -1
1836
1668
1546
-1
Attenuation coefficient/dB·Ȝ
0.479
0.692
0.095
Figure 3 shows the sound velocity profiles at 2000m, 4000m, 6000m, 8000m (x, y directions),
respectively. Different sound velocity profiles and complex ocean structure are shown in the figure
because of the different properties in different water bodies. The simulated sound ray tracks are shownin Figure 4.

Figure 3 - Sound velocity profiles at 2000m,
4000m, 6000m, 8000m (x, y directions)Figure 4 - Simulated sound ray tracks
between two points in the ocean model

The sound ray tracks are not kept in the initial direction and shuttle in different planes.
The travel time, the sound ray intensities and the reflection angles of the eigenray between any two
points can be calculated by the platform. We can get the sound level at certain point if we consider
the particular frequency, shown in Table 2.
Table 2 - Acoustic parameters of points calculated from (2000, 2000, 200) m
No.
x/m
y/m
z/m
Travel time/sSound level attenuation/dB
1
3000
3000
200
1.016
63.5
2
4000
4000
300
2.075
69.8
3
5000
5000
400
3.199
73.2
4
6000
6000
200
4.072
75.9
5
7000
7000
300
5.185
77.8
6
8000
8000
400
6.336
79.6

5. CALCULATED AMOUNTS OF DIFFERENT ALGORITHMS
The calculated amounts cannot be compared accurately for the different approaches of the different
algorithms. While the calculated amounts can be compared qualitatively from each algorithm's
principle.
The traditional eigenray algorithm uses Newton method or Runger Kutta method to solve the eigen
equation. It only can be used in each local layered structure. The calculated amount of this method isrelated to the number of the non-layered structures, which can be expressed as O (mnu). m is thenumber of the non-layered structures. n is the grid quantity of the local layered structure. u is thenumber of the iterations in Newton method or Runger Kutta method.
In the real-time simulation, the information of certain points is extracted directly from the energy3225

matrix by the sound ray track technology in steady-state sound field. The sound level attenuation canbe deduced by the information of the sound frequency, without other information, such as the numberof the non-layered structures, the grid quantity of the local layered structure and the process of solving
the eigen equation, etc. So the calculated amount of this method is deemed to be O (1).
The two series of the different methods' calculation time are shown below. The CPU mainfrequency is 2GHz. The operation system is Windows 7. The sound level attenuation from (2000,
2000, 200) m to (6000, 6000, 400) m is investigated in the 10000m×10000m×4000m 3D ocean model.
The sound velocity profiles are listed in Fig.3. Different densities of the discrete grids are used.

Figure 5 - Relationship of CPU time and quantity of sound velocity zone
Figure 5 shows that the CPU time of the steady-state sound field matrix method is much shorter
than that of the eigenray method.

6. DETECTION SENSITIVITY
The strength of the acoustic signal detected by the acoustic detection device has a strong correlation
with the location of the detection device itself. It can be seen from Equation (1), (2), (3), (4) that therange of the acoustic attenuation is strongly influenced by the sound ray direction. The change of thesound attenuation along the sound ray is mainly caused by the change of the structure of the soundray along the propagation direction, while the change of the sound attenuation perpendicular to thesound ray is mainly caused by the change of the structure of the sound ray paralleled to the sound ray.
According to the analysis above, the definition of the detection sensitivity between the acoustic
attenuation and the sound speed profile at a typical position is proposed as a specific value, or aderivative, which is related to the normalized dependent variable, such as signal strengths/detection
distance, etc., and the normalized independent variable, such as sound velocity profiles.
“Normalized” means changing the amount of the variable in proportion within the value range.
Under the condition of different sound velocity profile, the properties of the detection sensitivity
are as follows:(1)Under the condition of 0 gradient of sound velocity profile, sensitivity is 0 in both horizontal
and vertical direction.
(2)Under the condition of positive gradient of sound velocity profile, the sound ray derived fromthe sound source become an upward spiral. With the increase of gradient, the curvature of the wave
crest envelope line becomes higher, so as the absolute value of the detection sensitivity.
(3)Under the condition of negative gradient of sound velocity profile, the sound ray derived fromthe sound source become a downward spiral. With the increase of gradient, the curvature of the wave
crest envelope line becomes higher, so as the absolute value of the detection sensitivity.

7. COMPARISON OF SIMULATION AND TEST
7.1 Analysis of practical measurement results
The test was carried out in South China Sea in July, 2016. The average depth of the test area isabout 3170m. 64 effective hydrophones were deployed from 30m to 1000m in the vertical directionof the ocean. The testing ship transmitted 10s linear frequency modulated signal (center frequency 1
kHz, bandwidth 100Hz). The depth profile of the test area is shown in Figure 6. The depth fluctuatedwith the distance.

3226

Figure 6 - Depth profile of test area

Figure 7 - Sound attenuations got by test andsimulation
The depth of the deep sound channel is about 1050m. The sound velocity profile above 1000m
obtained from the test. The profile below 1000m is from a database. The sediments thickness in the
double-layered liquid seafloor is 120m. The other parameters are listed in Table 1.
The sound attenuations got by the test and the simulation from the sound source to different
positions are shown in Figure 7.
Figure 7 shows that the trend of the simulation is consistent with that of the test. The relative
difference between them gets bigger with the crosswise distance growth because the testing distanceis within the middle range of the convergence region. The relative difference grows up and the soundintensity gets weaker when the distance close to the middle range of the convergence region. The
average error of the sound attenuation is within 3dB by the sound ray track technology in steady-statesound field.
Overall, the sound velocity profiles of the test area in this paper are relatively stable and show
obvious regularity. More fluctuated sound velocity profiles may cause higher uncertainty and bigger
error of the sound attenuation in the algorithm.
7.2 Analysis of detection sensitivity results
According to the test results of acoustic attenuation, combining with the information of thedifferent parts of the sound velocity profile, the relationship of variable quantities of the normalized
velocity gradients and those of the detection sensitivities is shown as below.

Figure 8 - Relationship of variable quantities of normalized velocity gradients and normalized
detection sensitivities
The relationship mentioned above shows that the normalized detection sensitivity is higher whenthe variable quantity of the normalized sound velocity gradient is closer to 0 (negative gradient areahere). The normalized detection sensitivity tends to the minimum when the variable quantity of thenormalized sound velocity gradient is closer to 0.64 (sound velocity gradient close to 0). The degreeof the acoustic attenuation change gets lower. When the velocity gradient is positive, the value of thedetection sensitivity will grow higher.

3227

8. CONCLUSIONS
A numerical analysis method of sound ray propagation in complex undersea environment based on
non-coplanar sound ray and sound velocity gradient vector is presented in this paper. The travel timefrom the sound source to any concerned position can be calculated through GMM, derived from thesound velocity calculated by the sound velocity gradient vector algorithm proposed in this paper. The
sound level attenuation can be further deduced by the sound travel time. The results calculated fromthe algorithm well fit the data of the test. The algorithm can be well applied in the real-time simulationof the detection in seawater. The quantitative relationship between physical quantities is welldescribed by normalized detection sensitivity according to the results of acoustic attenuation test.

REFERENCES
1. Catipovic J A. Performance limitations in underwater acoustic telemetry. IEEE J. Oceanic Eng. 1990;
15(3): 56-65.
2. Estes L E, Fain G, Carvalho D. A shallow water channel characterization for underwater acoustic
telemetry. Proc: Conf Oceans; 1993.
3. Kelf M A. Hardware-in-the-loop simulation for undersea vehicle applications. The International Society
for Optical Engineering. 2001; 4366: 1- 12.
4. Yang D, et al. An improved nearly analytical discrete method: an efficient tool to simulate the seismic
response of 2-D porous structures. J. Geophys. Eng. 2007; 4: 40-52.
5. Sun Z, Ma Y, Tu Q. Fast search for eigenrays in stratified ocean. Applied Acoustics 1997;4: 7-14.
6. Etter P C. Underwater acoustic modeling and simulation. London: Spon Press; 2003.
7. Han J, Huang J, Cao H. Ocean acoustic channel simulator HJRAY and its application on underwatercommunication. Journal of System Simulation 2007; 19(1): 35-37.
8. Huang Y, Zhang J. A three-dimensional sound ray tracing algorithm based on wavefront traveltimes
interpolation. Acta Acustica 2008; 33(1): 21-27.
9. Kim S. 3-D Eikonal Solvers: First-arrival traveltimes. Geophysics 2002; 67(4): 1225-1231.

3228