A Method for Predicting Radiated Acoustic Field in Shallow Sea Based on Wave Superposition and Ray

The traditional free-space and half-space analysis method ignore the reflection of the upper and lower boundaries of shallow sea and are not suitable for analyzing shallow sea problems especially at high frequency. Hence, a method combining ray theory and wave superposition theory is proposed in this paper to predict the high-frequency radiated acoustic field in shallow water. The proposed method takes into account the effect of channel boundaries on the acoustic field and has good adaptability to complex channel environments and accuracy of the calculated acoustic field.


Introduction
The accurate prediction of acoustic radiation of underwater vehicles is an important engineering application. The accurate simulation of radiated noise is also important for acoustic field characteristics analysis and acoustic measurement research. The oceans surrounding China mostly consist of shallow water and have a typical shallow water waveguide environment. In the shallow sea, the traditional free-space analysis method is limited due to the acoustic reflection phenomenon of the upper and lower boundary surfaces, and the acoustic radiation prediction in the shallow sea becomes very complicated. The acoustic propagation model used in sea includes ray theory [1], normal model theory [2,3], wavenumber integration method [4], and parabolic equation method [5]. The advantage of the ray theory is that the physical meaning is intuitive and has good adaptability to high frequency problems. The normal wave theory is suitable for low frequency problems and provides accurate calculation in the far field. The wave number integration method is an integral transformation method of a horizontal layered medium and provides fast calculation speed. The parabolic equation method performs better in dealing with the problem of sound propagation in the environment related to distance. This kind of research method that uses a point source, ignores the vibration distribution of the structure surface and the directivity of the acoustic radiation. Hence, it is not appropriate to regard the elastic structure as a point source with the same source strength and directivity. The Finite Element Method (FEM) is only applicable to the near field region. However, for the far acoustic field, and with the high frequency condition, the number of grid elements in the FEM is too large to calculate, and the FEM is inapplicable in this scenario. The wave superposition method [6,7] has excellent structural adaptability. Most of the current research focuses on free space and half-free space, but there are few researchers who concentrate on the sea waveguide [8][9][10][11].
For the acoustic radiation problem of the elastic structure in the shallow sea, some combined models are developed under the condition that it is difficult to meet the calculation requirements by a single method. Shang proposed a method that combines the wave superposition method and normal mode Green's function method, which effectively predicted the acoustic radiation of large elastic structures under shallow sea [12]. This method uses the normal mode Green's function to calculate the far field pressure, while the normal mode green function has a large order at high frequencies, which has a greater impact on the computational efficiency. Hence, this method is not appropriate for calculation speed at high frequencies. Qian proposed a method that combines the finite element method with the parabolic equation method to obtain cylindrical field data in the near field active region, and a complete far field solution in the far-field passive region, using the parabolic equation method [13]. The parabolic equation method is a recursive method, and far field calculation leads to error accumulation. In addition, it is difficult to directly obtain the data of a near field cylindrical field.
In this paper, a new method is proposed for the prediction of an acoustic field radiated from the structure excited at a frequency. Based on the combination of wave superposition and ray theory, the structure source is equivalent to a set of virtual sources that are reconstructed by matching the surface normal velocity of the structure source. Finally, the complete acoustic field is obtained by superimposing sound radiation fields from all virtual sources using ray theory.

Theoretical Model
The proposed theoretical model consists of the combination of wave superposition theory and ray theory, corresponding to the acoustic inverse problem and the acoustic positive problem respectively.

Wave Superposition Theory
The acoustic field radiated by the elastic structure in the fluid environment meet the Sommerfeld radiation condition in the far field, and the fluid-solid coupling boundary satisfies the condition of the normal vibration velocity. Additional boundary conditions for the Helmholtz equation can be expressed as: where, v n is the normal velocity of the fluid medium, and v n + is the normal vibration velocity of the structural surface in the fluid medium. Equation (3) shows the Sommerfeld far field radiation condition. The equation can be solved by the Helmholtz integral formula. For any shape of radiator, the radiated pressure at the field is the integral of the internal volume source of the structure as shown as: where, P and Q are coordinates of the field point and virtual source in the global coordinate system, respectively. The ρ 0 shows medium density in field, ω is the angular frequency, q(Q) is strength of the virtual source, and G(P, Q) is the Green's function under the corresponding boundary conditions. A schematic diagram is shown in Figure 1. Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 13 For the typical channel environment (i.e., the sea surface is the Dirichlet boundary and the sea floor is the Neumann boundary), the normal mode Green's function in cylindrical coordinates can be shown as [14]: where, s z shows the source depth, D is depth of the sea, zn k and rn k represent the vertical and horizontal components of the wave number k , respectively, and   zn k z  is Eigen function with different forms under different boundary conditions. The parameter r is the horizontal distance, and N is the normal order. Introducing the linear Euler formula to transform Equation (4) into an integral representation of normal vibration velocity on the surface can be expressed as: where, Equation (6) relates the surface normal velocity to the Green's function. In theory, there is no limit for the placement of virtual sources and they can be located anywhere within the structure. Related studies have shown that the virtual source surface and the surface of the structure are in the same shape, which can improve the inversion accuracy of the configured equivalent source. It can be assumed that the virtual source is distributed over a virtual spherical shell surface  of thickness   inside the structure; hence, Equation (6) can be written as: where P is the node coordinate on the surface of the structure. Dividing the virtual spherical shell into I parts, each surface of the structure is shown by i  . If each spherical shell crown i  is small enough, the summation coefficient   can be regarded as constant, and the normal vibration velocity at different points s P on the structural surface can be approximately written as: where, i Q represents the coordinate of the ith virtual source, and   i q Q represents the strength of the ith source on the virtual surface. Hence, Equation (8) can be written in the matrix as: The structural surface normal velocity matrix U in Equation (9) can be expressed by the source strength matrix Q and transfer matrix D , where D can be expressed as: For the typical channel environment (i.e., the sea surface is the Dirichlet boundary and the sea floor is the Neumann boundary), the normal mode Green's function in cylindrical coordinates can be shown as [14]: where, z s shows the source depth, D is depth of the sea, k zn and k rn represent the vertical and horizontal components of the wave number k, respectively, and Ψ(k zn z) is Eigen function with different forms under different boundary conditions. The parameter r is the horizontal distance, and N is the normal order. Introducing the linear Euler formula to transform Equation (4) into an integral representation of normal vibration velocity on the surface can be expressed as: where, Equation (6) relates the surface normal velocity to the Green's function. In theory, there is no limit for the placement of virtual sources and they can be located anywhere within the structure. Related studies have shown that the virtual source surface and the surface of the structure are in the same shape, which can improve the inversion accuracy of the configured equivalent source. It can be assumed that the virtual source is distributed over a virtual spherical shell surface σ of thickness δ τ inside the structure; hence, Equation (6) can be written as: where P is the node coordinate on the surface of the structure. Dividing the virtual spherical shell into I parts, each surface of the structure is shown by σ i . If each spherical shell crown σ i is small enough, the summation coefficient δ τ can be regarded as constant, and the normal vibration velocity at different points P s on the structural surface can be approximately written as: where, Q i represents the coordinate of the ith virtual source, and q(Q i ) represents the strength of the ith source on the virtual surface. Hence, Equation (8) can be written in the matrix as: The structural surface normal velocity matrix U in Equation (9) can be expressed by the source strength matrix Q and transfer matrix D, where D can be expressed as: It should be noted that in practical applications, the matrix D may be an ill-conditioned matrix, and slight disturbance of vibration velocity will cause a dramatic change in virtual source strength. According to the acoustic inverse problem observed in a literature study, the Tikhonov regularization method based on L curve can effectively suppress the disturbance and effectively deal with the ill-conditioned matrix before solving the inverse problem [15]. For the same discrete processing method in Equation (4), the pressure at point P can be expressed as the superposition of I discrete point source radiated acoustic fields as shown as:

Ray Theory
The ray solution exists in the form p(x, y, z, t) = A(x, y, z)e j[ωt−k 0 ϕ(x,y,z)] . Where k 0 = ω/c 0 represents the wave number at the initial position, c 0 is the sound speed of the initial position, n(x, y, z) = c 0 /c(x, y, z) is the index of refraction, and k = k 0 n (x, y, z) is the wave number at position (x, y, z). The A(x, y, z) represents the pressure amplitude and k 0 ϕ(x, y, z) shows the pressure phase. Substituting the ray solution into the wave equation, the following two equations can be obtained as [16]: When ∇ 2 A/A is much smaller than k 2 , then ∇ · A 2 ∇ϕ = 0 In ray theory, Equations (14) and (15) are called as eikonal equation and intensity equation respectively, and are the two basic equations of ray theory [17]. In ray theory, the propagation of waves are usually regarded as the propagation of an infinite number of acoustic rays. Each ray is perpendicular to the equiphase surface, and the propagation path in the sea is affected due to the medium and boundary. The sound ray carries energy, and the energy in each sound ray does not lose to the outside of the sound ray. These sound rays contribute to the acoustic field in ray theory.
Eikonal equation combined with Snell's law, to obtain the trajectory of the ray, satisfies the following relationship [18] as shown below: After giving the initial glancing angle α 0 and the sound speed function c(z), the solution of the above equation can achieve ray tracking and acquire ray trajectories. The eigenray are determined in the distribution trend of the rays, which are special rays that directly supply the power at a certain point in the acoustic field. The energy of each sound ray can be derived from the intensity equation as: where W represents the radiated power in unit solid angle. The energy of these rays reflects due to the boundary and distributes in the acoustic field in a steady state. After determining the position and strength of virtual sources by wave superposition method, the acoustic field is calculated by ray theory. Then the prediction of the structured radiation acoustic field is realized.

Accuracy Analysis of Combined Wave Superposition Method
The combined wave superposition method (combined WSM) is an approximation method, and the accuracy of the approximation is needed to be discussed. It approximates not only the normal mode order of the Green's function in sea, but also the process of inverting the virtual source strength. The ray theory is also an approximate acoustic propagation model in the acoustic positive problem. In the following sections, we will analyze the error caused by the local approximation in combined WSM, and some characteristics of the whole method in numerical calculation.

Accuracy Analysis of Green's Function in Channel
The basis of the wave superposition method is the Green's function. The accuracy of the Green's function, especially in the near field, directly determines the accuracy of the transfer matrix. The Green's function is expressed as the sum of the normal modes. The high-order normal mode has little contribution to the far field; however, it has a great influence on the near field. Therefore, the normal order determines the near field accuracy of the Green's function. When the order is too small, the result of the Green's function is inaccurate in the near field, but gradually approaches the exact value as the distance increases.
In order to illustrate the near field accuracy of the Green's function, a case is selected for analysis. In the shallow sea of 200 m depth, the sea surface is the Dirichlet boundary and the sea floor is the Neumann boundary. A point source with a frequency of 20 Hz is located at a depth of 25 m. According to the boundary conditions, the following constraints exist in Equation (5) as shown as: The eigen function Ψ that satisfies the constraints is Ψ(x) = sin(x). The specific form of the Green function can be obtained according to the constraint conditions, which is the analytical solution under boundary condition as shown as: where, z s shows the source depth, D indicates the depth of the sea, and k zn and k rn represent the vertical and horizontal components of the wave number k, respectively. The parameter r is the horizontal distance, and n is the normal order. The model is built using finite element software COMSOL and compared with the Green's function method when the mesh density is sufficiently dense, and the results are shown in Figure 2.
It can be seen that under far field conditions, the calculation result of the Green's function is nearly consistent with the finite element method. In fact, the situation is different under near field conditions, which is what needs to be discussed next.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 13 It can be seen that under far field conditions, the calculation result of the Green's function is nearly consistent with the finite element method. In fact, the situation is different under near field conditions, which is what needs to be discussed next. We assume the error threshold 1 dB, that is, the accuracy of the near field calculation is characterized by the horizontal distance when the absolute error between the transmission loss of the Green's function method and the finite element method is less than 1 dB. The further the distance is, the more accurate are the results of the near field of the Green's function we obtain. Figure 3 below shows the near field accuracy of the Green's function after changing the normal mode order of the Green's function at different frequencies or different depths, while the sea environment and point source parameters are consistent with the above case. The overall trend is that if the order is larger, the results of the near field of the Green's function are more accurate. In order to ensure that the wave superposition method can effectively reverse the acoustic field information, it needs to increase the order and to ensure that the horizontal distances from the equivalent sources to the positions of vibration velocity point are not too close. As the frequency increases, the energy of the acoustic field concentrates on the low-order normal modes, and loss due to the lack of energy of the high-order is also reduced. However, in order to get the accuracy of the near-field, it is necessary to keep the order large enough. In terms of greater depth, the requirement of the normal order is higher in this method. Therefore, theoretically the accuracy of the wave superposition method can be guaranteed with a small depth.  We assume the error threshold 1 dB, that is, the accuracy of the near field calculation is characterized by the horizontal distance when the absolute error between the transmission loss of the Green's function method and the finite element method is less than 1 dB. The further the distance is, the more accurate are the results of the near field of the Green's function we obtain. Figure 3 below shows the near field accuracy of the Green's function after changing the normal mode order of the Green's function at different frequencies or different depths, while the sea environment and point source parameters are consistent with the above case.

Applicability of Rays Superimposed
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 13 It can be seen that under far field conditions, the calculation result of the Green's function is nearly consistent with the finite element method. In fact, the situation is different under near field conditions, which is what needs to be discussed next. We assume the error threshold 1 dB, that is, the accuracy of the near field calculation is characterized by the horizontal distance when the absolute error between the transmission loss of the Green's function method and the finite element method is less than 1 dB. The further the distance is, the more accurate are the results of the near field of the Green's function we obtain. Figure 3 below shows the near field accuracy of the Green's function after changing the normal mode order of the Green's function at different frequencies or different depths, while the sea environment and point source parameters are consistent with the above case. The overall trend is that if the order is larger, the results of the near field of the Green's function are more accurate. In order to ensure that the wave superposition method can effectively reverse the acoustic field information, it needs to increase the order and to ensure that the horizontal distances from the equivalent sources to the positions of vibration velocity point are not too close. As the frequency increases, the energy of the acoustic field concentrates on the low-order normal modes, and loss due to the lack of energy of the high-order is also reduced. However, in order to get the accuracy of the near-field, it is necessary to keep the order large enough. In terms of greater depth, the requirement of the normal order is higher in this method. Therefore, theoretically the accuracy of the wave superposition method can be guaranteed with a small depth.  The overall trend is that if the order is larger, the results of the near field of the Green's function are more accurate. In order to ensure that the wave superposition method can effectively reverse the acoustic field information, it needs to increase the order and to ensure that the horizontal distances from the equivalent sources to the positions of vibration velocity point are not too close. As the frequency increases, the energy of the acoustic field concentrates on the low-order normal modes, and loss due to the lack of energy of the high-order is also reduced. However, in order to get the accuracy of the near-field, it is necessary to keep the order large enough. In terms of greater depth, the requirement of the normal order is higher in this method. Therefore, theoretically the accuracy of the wave superposition method can be guaranteed with a small depth.

Applicability of Rays Superimposed
For the Green's function requiring high-order superposition, the calculation time of each virtual source propagating the acoustic field increases rapidly with the increase of the calculated range and Appl. Sci. 2020, 10, 917 7 of 13 order, making the calculation efficiency low. The large number of virtual sources inversion by the wave superposition method makes the superposition with the Green's function not the most ideal choice. The ray theory performs ray tracing in the entire acoustic field and extracts the eigenray, which is less computationally expensive. Since the ray theory is derived under the high frequency approximation, it is necessary to verify that the ray solution can meet the accuracy requirements of the acoustic field superposition at different frequencies and depths.
The channel is a typical shallow sea waveguide (i.e., the sea surface boundary is the Dirichlet boundary; the sea bottom boundary is the Neumann boundary), and the six point sources are located at depths of 25 to 30 m at intervals of 1 m. The depth of water and frequency are changed, and the acoustic field is superimposed by the ray theory and the Green's function method, respectively. The field comparison chart after superposition is shown in Figure 4 as: For the Green's function requiring high-order superposition, the calculation time of each virtual source propagating the acoustic field increases rapidly with the increase of the calculated range and order, making the calculation efficiency low. The large number of virtual sources inversion by the wave superposition method makes the superposition with the Green's function not the most ideal choice. The ray theory performs ray tracing in the entire acoustic field and extracts the eigenray, which is less computationally expensive. Since the ray theory is derived under the high frequency approximation, it is necessary to verify that the ray solution can meet the accuracy requirements of the acoustic field superposition at different frequencies and depths.
The channel is a typical shallow sea waveguide (i.e., the sea surface boundary is the Dirichlet boundary; the sea bottom boundary is the Neumann boundary), and the six point sources are located at depths of 25 to 30 m at intervals of 1 m. The depth of water and frequency are changed, and the acoustic field is superimposed by the ray theory and the Green's function method, respectively. The field comparison chart after superposition is shown in Figure 4 as: The ray solution will have a slight oscillation near the true value at low frequency, and the superposition of multiple point sources will amplify the amplitude of the oscillation. A large number of virtual sources can cause an increase in the amplitude of the oscillation. The amplitude of the oscillations is different at different frequencies and depths. The error caused by oscillation is characterized by mean square error (MSE).   The ray solution will have a slight oscillation near the true value at low frequency, and the superposition of multiple point sources will amplify the amplitude of the oscillation. A large number of virtual sources can cause an increase in the amplitude of the oscillation. The amplitude of the oscillations is different at different frequencies and depths. The error caused by oscillation is characterized by mean square error (MSE).
where M represents the number of measurement points, and TL R and TL G are the transmission losses obtained by ray theory and the Green's function methods respectively. The table below gives the MSE for certain frequencies and depths.
Appl. Sci. 2020, 10, 917 8 of 13 In Table 1, in the case where the frequency is too low or the depth is too small, there is a large error in ray theory. At high frequency, the solution of the ray theory is very accurate. The field of virtual source can be directly superimposed, and the result is reliable. As a whole, the oscillation of the ray theory has a limited influence on the calculation results. The performance at high frequency is superior, and the application range of the wave superposition method can be expanded.

Calculation of Acoustic Radiation Field by Combined Method
The key to the combined method is the solution of the acoustic inverse problem. The results of the virtual source inversion have a significant impact on the calculation of the acoustic field [19]. Research shows that the accuracy of wave superposition inversion has a great relationship with virtual source configuration. These researches are based on free space or half space Green's function, but the Green's function in sea is a cylinder functions, which has its own particularity. In order to illustrate the calculation accuracy of the method proposed in this paper and analyze the influence of different virtual source configuration methods, an excited spherical shell is selected to calculate its far field pressure using the proposed method and compared with the finite element model in COMSOL software. The reason why this axisymmetric structure was selected is that the calculation distance of the three-dimensional field by COMSOL is limited by the performance of the computer. The selected spherical shell has a radius of 2 m and a thickness of 0.1 m. The Material elastic Young's modulus is E = 2 × 10 11 N/m 2 and Poisson's ratio is σ = 0.3; density is ρ s = 7850 kg/m 3 . The outside medium is water, in which the density is ρ 1 = 1024 kg/m 3 , and the speed of sound is c = 1500 m/s. The inside of the spherical shell is a vacuum. The center of the spherical shell is located at a water depth of 25 m. A normal force of 1 N with a frequency of 100 Hz is applied at the bottom of the spherical shell. The depth of the sea is 200 m, and the sea surface and seabed boundary conditions satisfy typical channel conditions. The schematic diagram of the COMSOL model is shown in Figure 5. The sound pressure level (SPL) of the radiated field of this model is calculated using COMSOL. As the spherical shell is circumferentially symmetrical, the vibration velocity of the two-dimensional structure surface calculated by COMSOL can be scanned every 36 degrees in the circumferential direction to obtain the three-dimensional structure surface vibration velocity distribution. We also calculate SPL of the radiated field using combined WSM based on 3D vibration velocity distribution data. Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 13 (a) (b) The virtual sources are distributed on a surface that has the same shape as the surface of the structure, which is an appropriate virtual source configuration. The virtual source surface is generally obtained by shrinking the surface of the structure inwardly. Here the radius ratio is chosen to be 0.5. Except for the upper and lower vertices, the virtual source surface is divided into 10 and 6 intervals in the circumferential and vertical directions respectively, that is, a total of 72 virtual source points (10 7  virtual sources) are created, as shown in Figure 6. The sound pressure level is calculated at different depths as shown in Figure 7: The virtual sources are distributed on a surface that has the same shape as the surface of the structure, which is an appropriate virtual source configuration. The virtual source surface is generally obtained by shrinking the surface of the structure inwardly. Here the radius ratio is chosen to be 0.5. Except for the upper and lower vertices, the virtual source surface is divided into 10 and 6 intervals in the circumferential and vertical directions respectively, that is, a total of 72 virtual source points (10 × 7 virtual sources) are created, as shown in Figure 6. The virtual sources are distributed on a surface that has the same shape as the surface of the structure, which is an appropriate virtual source configuration. The virtual source surface is generally obtained by shrinking the surface of the structure inwardly. Here the radius ratio is chosen to be 0.5. Except for the upper and lower vertices, the virtual source surface is divided into 10 and 6 intervals in the circumferential and vertical directions respectively, that is, a total of 72 virtual source points (10 7 × virtual sources) are created, as shown in Figure 6. The sound pressure level is calculated at different depths as shown in Figure 7:   Theoretically, when the number of virtual sources is large, the dispersion degree of the volume source will be higher, and the inversion of the sound field will be more accurate. In fact, for the where, K is the number of measurement points. The p r and p f are the sound pressures generated by combined WSM and FEM, at the ith measurement point, respectively. The 14 × 11 and 10 × 7 virtual sources (excluding the upper and lower vertices) are selected respectively, to form a virtual source surface, as shown in Figure 8, and analyzed the effect of the radius ratio λ on the calculation result at the receive depth of 25 m, as shown in Figure 9.  Theoretically, when the number of virtual sources is large, the dispersion degree of the volume source will be higher, and the inversion of the sound field will be more accurate. In fact, for the Theoretically, when the number of virtual sources is large, the dispersion degree of the volume source will be higher, and the inversion of the sound field will be more accurate. In fact, for the situation analyzed in this section, after the number of virtual sources meets the basic requirements of the calculation, it does have some impact on the accuracy of the calculation when the radius ratio is in the range of 0.7 to 1. The impact of this factor is far greater than the impact of the number of virtual sources. Figure 9 shows that in the case where the number of virtual sources is constant, the radius ratio is too small or too large, which leads to an increase in error. When the radius ratio approaches 1, because the Green's function is a cylindrical function, the horizontal distance between the virtual source and the surface node of the structure becomes too small, and the calculated D matrix singularity increases. When the radius ratio is too small, the difference in surface vibration velocity is not sufficiently reflected on the virtual source. In a word, the appropriate radius ratio plays a critical role in the accuracy of the calculation. For the spherical shell, the optimum ratio of the radius ratio is in the range of 0.5 to 0.9. When the shape or other parameters of the structure change, the configuration of the virtual sources should be adjusted.
The 10 × 7 virtual sources (excluding the upper and lower vertices) and a radius ratio of 0.5 are selected, which is same with the model in Figure 7 to get a two-dimensional color figure of SPL in 100 Hz using combined WSM and FEM respectively as shown in Figure 10.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 13 situation analyzed in this section, after the number of virtual sources meets the basic requirements of the calculation, it does have some impact on the accuracy of the calculation when the radius ratio is in the range of 0.7 to 1. The impact of this factor is far greater than the impact of the number of virtual sources. Figure 9 shows that in the case where the number of virtual sources is constant, the radius ratio is too small or too large, which leads to an increase in error. When the radius ratio approaches 1, because the Green's function is a cylindrical function, the horizontal distance between the virtual source and the surface node of the structure becomes too small, and the calculated D matrix singularity increases. When the radius ratio is too small, the difference in surface vibration velocity is not sufficiently reflected on the virtual source. In a word, the appropriate radius ratio plays a critical role in the accuracy of the calculation. For the spherical shell, the optimum ratio of the radius ratio is in the range of 0.5 to 0.9. When the shape or other parameters of the structure change, the configuration of the virtual sources should be adjusted. The 10 7  virtual sources (excluding the upper and lower vertices) and a radius ratio of 0.5 are selected, which is same with the model in Figure 7 to get a two-dimensional color figure of SPL in 100 Hz using combined WSM and FEM respectively as shown in Figure 10.
The field interference fringes in Figure 10a are clear, and the results are similar to the finite element results in Figure 10b. Same as in Figure 10b, Figure 10a can also describe the phenomenon caused by the interaction of radiated sound with the sea surface and the sea bottom during the propagation of the sound, that is, the sound rays reflect back and forth between the sea surface and the sea bottom. However, the difference is that compared with Figure 10b, the peak and trough amplitudes of Figure 10a are not clear enough, which is also reflected in Figure 7. It should be noted that in order for combined WSM to accurately predict the radiated field of the structure, the environment of the sea needs to meet the applicable range of ray theory. The following relationships can be used as a guide for high-frequency approximation of ray theory [1] as shown as: where f is the frequency, c is sound speed in water, and D is the depth of the water. This relationship can also be used as guidance for the applicable conditions of combined WSM. Although the guidance given by this relationship is not necessarily accurate, it can indeed reflect the potential of combined WSM for the calculation of high-frequency problems.

Conclusions
In this paper, a new method is proposed for the prediction of radiated acoustic field in shallow sea. Ray theory is a high frequency approximation theory with high accuracy and high computational efficiency. On the basis of the wave superposition method, the application range of the wave superposition method is expanded by introducing the ray theory. Since the Green's function in sea is The field interference fringes in Figure 10a are clear, and the results are similar to the finite element results in Figure 10b. Same as in Figure 10b, Figure 10a can also describe the phenomenon caused by the interaction of radiated sound with the sea surface and the sea bottom during the propagation of the sound, that is, the sound rays reflect back and forth between the sea surface and the sea bottom. However, the difference is that compared with Figure 10b, the peak and trough amplitudes of Figure 10a are not clear enough, which is also reflected in Figure 7.
It should be noted that in order for combined WSM to accurately predict the radiated field of the structure, the environment of the sea needs to meet the applicable range of ray theory. The following relationships can be used as a guide for high-frequency approximation of ray theory [1] as shown as: where f is the frequency, c is sound speed in water, and D is the depth of the water. This relationship can also be used as guidance for the applicable conditions of combined WSM. Although the guidance given by this relationship is not necessarily accurate, it can indeed reflect the potential of combined WSM for the calculation of high-frequency problems.

Conclusions
In this paper, a new method is proposed for the prediction of radiated acoustic field in shallow sea. Ray theory is a high frequency approximation theory with high accuracy and high computational efficiency. On the basis of the wave superposition method, the application range of the wave superposition method is expanded by introducing the ray theory. Since the Green's function in sea is a cylindrical function, a high-order superposition is required to ensure accuracy, while the distance between the virtual source surface and the surface of the structure cannot be too close. In addition, the structural surface vibration velocity nodes and the virtual sources should avoid overlapping in the vertical direction as much as possible, which should be paid attention to, while constructing the virtual source surface. After comparison with examples, it is verified that the proposed method is accurate. It can be used to analyze the acoustic radiation of complex structures under water and have potential for the calculation of high-frequency problems.

Conflicts of Interest:
The authors declare no conflict of interest.