Sound wave propagation calculation method, sound wave propagation calculation program and sound wave propagation calculation device

Abstract

To provide a sound wave propagation calculation method, a sound wave propagation calculation program and a sound wave propagation calculation device that can suppress excess or deficiency of a calculation angle to a propagation angle and efficiently calculate from a viewpoint of calculation accuracy and calculation amount.SOLUTION: A sound wave propagation calculation method calculates sound pressure by the PE (Parabolic Equation) method using the Split-Step Pade method. The sound wave propagation calculation method comprises: an optimum parameter calculating step which obtains optimum parameters including an approximation order and distance increment that differ depending on a calculated distance; an approximation coefficient calculating step which calculates an approximation coefficient using the optimum parameters; and a sound pressure calculating step which calculates the sound pressure using the approximation coefficient.SELECTED DRAWING: Figure 11

Description

[Technical Field]

[0001] The present invention relates to a sound wave propagation calculation method, a sound wave propagation calculation program, and a sound wave propagation calculation device for calculating sound wave propagation in underwater acoustics. [Background technology]

[0002] Conventionally, simulation of acoustic wave propagation characteristics in the ocean has been performed for various purposes, such as estimating marine environmental characteristics using sound or estimating target positions using sound emitted from targets. One method for simulating acoustic wave propagation characteristics is the "Parabolic Equation Method" (hereinafter referred to as the "PE Method") (see, for example, Non-Patent Document 1). The PE Method sequentially determines sound pressure with respect to distance based on the relationship obtained from the Helmholtz equation, which satisfies the sound pressure and wave number with respect to distance and depth.

[0003] Among the PE methods, the "Split-Step Pade method" is known as a method for calculating sound pressure using approximation (see, for example, Non-Patent Document 2). In the Split-Step Pade method, the accuracy of the calculated sound pressure and the amount of calculation required to calculate the sound pressure can be adjusted by setting the degree of approximation. [Prior art documents] [Non-patent literature]

[0004] [Non-Patent Document 1] FB Jensen et al., “6 Parabolic Equations,” in Computational Ocean Acoustics (Second Edition), 2011 [Non-patent document 2] Ishiwata et al., “Development of a parabolic sound propagation model for shallow water sound field analysis,” Journal of the Acoustical Society of Japan, 57(12), 2001 Summary of the Invention [Problem to be solved by the invention]

[0005] In the case of underwater sound propagation in a distance-independent environment where the sound speed profile of the sea and seafloor, the density and attenuation rate of the seafloor, and the seafloor topography do not change with distance, and where the sound speed is constant and does not depend on depth, the longer the propagation distance, the less the components with high depression angles contribute to the propagation, and the components with low depression angles become dominant. The conventional PE method uses this property to set the calculation angle from the propagation angle at the propagation distance for which the sound pressure is to be calculated.

[0006] However, when calculating sound pressure over a wide range from close to long distances in a distance-independent, constant sound speed environment, narrowing the calculation angle on the condition that calculations can be performed with sufficient accuracy at long distances poses the problem of worsening the calculation accuracy of sound wave propagation because components with high depression angles at close distances are not taken into account.On the other hand, if the calculation angle is set wide enough on the condition that calculations can be performed with sufficient accuracy at close distances, there is the problem of increased calculation load due to the calculation of components with high depression angles at long distances that do not actually need to be calculated.

[0007] Furthermore, in a distance-dependent environment where the sound speed profile or seabed topography changes with distance, the propagation angle narrows globally with distance, but the narrowing characteristic of the propagation angle varies with distance, making it difficult to estimate the propagation angle at the desired propagation distance.

[0008] The present invention has been made in light of the above-mentioned problems, and aims to provide a sound wave propagation calculation method, a sound wave propagation calculation program, and a sound wave propagation calculation device that can suppress excess or deficiency of the calculation angle relative to the propagation angle and perform calculations efficiently in terms of calculation accuracy and calculation amount. [Means for solving the problem]

[0009] The sound wave propagation calculation method according to the present invention is a sound wave propagation calculation method for calculating sound pressure by the PE (Parabolic Equation) method using the Split-Step Pade method, and includes an optimal parameter calculation step for determining optimal parameters including an approximation order and distance increment that differ depending on the calculation distance, an approximation coefficient calculation step for calculating approximation coefficients using the optimal parameters, and a sound pressure calculation step for calculating sound pressure using the approximation coefficients.

[0010] A sound wave propagation calculation program according to the present invention causes a processor of a sound wave propagation calculation device to execute the sound wave propagation calculation method.

[0011] Furthermore, the sound wave propagation calculation device according to the present invention is a sound wave propagation calculation device that calculates sound pressure by the PE (Parabolic Equation) method using the Split-Step Pade method, and includes an optimal parameter calculation unit that finds optimal parameters including an approximation order and distance increment that differ depending on the calculation distance, an approximation coefficient calculation unit that calculates approximation coefficients using the optimal parameters, and a sound pressure calculation unit that calculates sound pressure using the approximation coefficients. [Effects of the Invention]

[0012] According to the present invention, by calculating sound pressure using optimal parameters including an approximation order and distance increment that vary depending on the calculation distance, it is possible to suppress excess or deficiency of the calculation angle relative to the propagation angle, and to perform calculations efficiently in terms of calculation accuracy and calculation amount. [Brief explanation of the drawings]

[0013] [Figure 1] FIG. 10 is a diagram illustrating the relationship between a differential operator and a phase error. [Figure 2] FIG. 10 is a diagram showing the relationship between the approximation order, the distance step, and the phase error. [Figure 3] FIG. 1 is an image diagram of sound wave propagation calculation processing using the PE method in the prior art. [Figure 4] FIG. 1 is a schematic diagram illustrating the configuration of a sound wave propagation calculation device according to the prior art. [Figure 5]FIG. 10 is an image diagram of calculation of an initial sound field using the mirror image sound source method. [Figure 6] 1 is a flowchart of a sound wave propagation calculation process in the prior art. [Figure 7] FIG. 1 is a physical image diagram of calculation conditions in a sound wave propagation calculation process according to the prior art. [Figure 8] 10A and 10B are diagrams illustrating the difference in sound wave propagation distance depending on the depression angle. [Figure 9] FIG. 10 is a diagram showing the characteristics of the propagation angle with respect to the propagation distance. [Figure 10] 1 is a schematic configuration diagram of a sound wave propagation calculation device according to a first embodiment. [Figure 11] 4 is a flowchart of a sound wave propagation calculation process according to the first embodiment. [Figure 12] 4 is a flowchart showing the flow of an optimal parameter calculation process in the first embodiment. [Figure 13] FIG. 3 is a conceptual diagram of an angle table used in the optimum parameter calculation process in the first embodiment. [Figure 14] 1 is an image diagram of a sound wave propagation calculation process in a sound wave propagation calculation device according to the first embodiment. FIG. DETAILED DESCRIPTION OF THE INVENTION

[0014] First, we will explain the sound wave propagation calculation method using the general Split-Step Pade method. Basically, the wave model calculates sound wave propagation by solving the Helmholtz equation shown in Equation (1).

number

[0015] Here, ▽ 2 is a vector operator expressed by the following equation (2), where p represents sound pressure and k represents wave number.

number

[0016] In the PE (Parabolic Equation) method, the cylindrical diffusion component and the phase rotation component due to the central wave number are calculated using the Hankel function H0 (1) (k0r), and the other components are the envelope sound pressure ψ(r,z), and we assume a solution to the Helmholtz equation as shown in equation (3).

number

[0017] Here, r is the horizontal distance (hereafter, horizontal distance will be referred to as "distance"), z is the depth, and k0 is the central wave number. The central wave number k0 is the reference horizontal wave number, and is given by k0 = ω / c0. Also, c0 is the reference sound speed, ω (= 2πf) is the angular frequency, and f is the frequency. H0 (1) Since (k0r) can be easily calculated if the distance r is known, solving equation (3) mainly involves calculating the envelope sound pressure ψ(r,z). In the PE method, the envelope sound pressure ψ(r,z) is calculated by solving the differential equations (4) and (5).

number

number

[0018] Here, q represents the differential operator, Δr represents the distance increment, and n represents the refractive index. However, in the differential equation format of Equation (4), the differential operator q is contained within an exponential function, making it difficult to solve in this form. Therefore, in the PE method, the exponential function terms are approximated using Pade approximation shown in Equation (6). Pade approximation is an approximation method that expresses the differential operator q as a sum of rational functions as in Equation (6) by treating it as if it were a variable.

number

[0019] where a j,m ,b j,m(j=1,...m) represents the approximation coefficient, and m represents the order of approximation. Equation (7) is obtained by substituting equation (6) into equation (4).

number

[0020] In the PE method using the Split-Step Pade method, the envelope sound pressure ψ(r+Δr,z) at a distance Δr from the given envelope sound pressure ψ(r,z) is calculated by solving equation (7). Hereafter, the envelope sound pressure will be referred to as sound pressure ψ.

[0021] Next, we will explain the accuracy of the Pade approximation shown in equation (6). Figure 1 is a diagram showing the relationship between the differential operator and the phase error. Figure 1 shows the phase error for the differential operator q when the differential operator q is treated as if it were a variable. As shown in Figure 1, there is no error when q = 0, and the error increases as q moves away from 0. Next, we will explain the relationship between the differential operator q and the elevation angle θ. Assuming a plane wave solution such as equation (8) and applying the differential operator q, we obtain equation (9).

number

number

[0022] where n is the refractive index, k is the wave number, k0 is the central wave number, and k r is the horizontal wave number, k z is the vertical wave number, and θ is the elevation angle. From equation (9), we can see that the differential operator q has a relationship with the elevation angle θ as shown in equation (10).

number

[0023] From the above relational expression, we can see that the differential operator q depends on the elevation angle θ, with q approaching 0 as the propagation approaches the horizontal direction (0°), and q moving further away from 0 as the propagation approaches the vertical direction (±90° upwards is -, downwards is +). As shown in Figure 1, the closer q approaches 0, the smaller the approximation error, and the further q moves away from 0, the larger the approximation error, so the error is smaller for propagation at low elevation angles (near 0°) and increases as propagation at high elevation angles (near ±90°).

[0024] In other words, to minimize the error over a wide angular range, it is necessary to set the approximation order m and the distance increment Δr so that the approximation error is small even when q is large. Figure 2 shows the relationship between the approximation order and distance increment and the phase error. Figure 2(a) shows a comparison of the phase error when the distance increment Δr is fixed at Δr = 10λ / 3 and the approximation order m is changed from 1 to 15. Figure 2(b) shows a comparison of the phase error when the approximation order m is fixed at 4 and the distance increment Δr is changed from 16·(λ / 3), 15·(λ / 3), …, 1·(λ / 3). Here, λ (λ = c0 / f) is the wavelength. Figure 2 shows that to minimize the approximation error, it is necessary to either increase the approximation order m or decrease the distance increment Δr. However, both methods are expected to increase the amount of calculations.

[0025] Finally, the tolerance ε max The angle range that can be calculated within ε max When ≡0.002 is set, the angle that can be calculated within the allowable error for the combination of (m, Δr) = (4, 10λ / 3) is estimated to be 40° from Figure 2(a). For the combination of approximation order m and distance step Δr, the phase error for the elevation / depression angle θ is ε(m, Δr, θ). At this time, the maximum elevation / depression angle θ that satisfies equation (11) is called the "calculated angle θ". cal "

number

[0026] Next, the processing performed by the conventional sound wave propagation calculation device 500 will be explained using as an example the calculation of sound pressure in a distance-independent environment where the sound speed profile in the ocean and on the seabed, the density and attenuation rate of the seabed, and the seabed topography do not change with distance. In the following explanation, () indicates a function of the variable in (), and [] indicates an index reference of a column vector. For example, sound pressure ψ(r,z) indicates that it is a function of distance r and depth z, and ψ ri =[ψ(r i ,z1),…,ψ(r i ,z N )] T is the discretized depth sequence Z=[z1, ..., z N ] T The calculated distance r i The sound pressure vector at ψ ri [i z ] is the i z It shows the th element. ri and Z are column vectors. Hereafter, column vectors will be referred to as vectors.

[0027] First, we will explain the general processing flow of the PE method. Figure 3 is an illustration of the sound wave propagation calculation process using the PE method in the conventional technology. The PE method is a method for calculating the sound pressure ψ(r + Δr, z) at a distance Δr from a given sound pressure ψ(r, z). Δr is a parameter determined during initial setup. In the sound wave propagation calculation process of the conventional technology, a vertical sound field called the initial sound field is first created using the set conditions, and the vertical sound field at the first distance step is calculated using the created initial sound field. In Figure 3, the vertical sound field is the sound pressure for one row enclosed by the dashed rectangle. Furthermore, in the second and subsequent steps, the vertical sound field is calculated using the calculation results from the previous step. In this way, the PE method calculates the sound pressure sequentially while updating the distance at each step.

[0028] FIG. 4 is a schematic diagram of a sound wave propagation calculation device 500 according to the related art. As shown in FIG. 4, the sound wave propagation calculation device 500 according to the related art includes an initial sound field calculation unit 51, an approximation coefficient calculation unit 52, a sound pressure calculation unit 53, an update processing unit 54, a determination unit 55, and a storage unit 56. The initial sound field calculation unit 51, the approximation coefficient calculation unit 52, the sound pressure calculation unit 53, the update processing unit 54, and the determination unit 55 of the sound wave propagation calculation device 500 are functional units realized by dedicated equipment or dedicated processing circuits. Alternatively, the sound wave propagation calculation device 500 may include a processor such as a CPU, and the processor may execute a sound wave propagation calculation program stored in the storage unit 56 to realize each functional unit of the sound wave propagation calculation device 500. Alternatively, each functional unit of the sound wave propagation calculation device 500 may be realized by a combination of dedicated equipment or dedicated circuitry and software.

[0029] The storage unit 56 is composed of a volatile memory such as a RAM, a DRAM, or an SRAM, a non-volatile memory such as a ROM, a flash memory, or a hard disk drive, or a combination thereof. The storage unit 56 stores a sound wave propagation calculation program executed by the processor of the sound wave propagation calculation device 500, various data such as parameters required for executing the sound wave propagation calculation program, etc.

[0030] The initial sound field calculation unit 51 calculates the initial sound field distance r s , source depth z source , sound speed profile c profile =[c1,…,c Mw ,…,c MB ] T , the decay rate profile γ profile =[γ1,…,γ Mw ,…,γ MB ] T , depth series Z=[z1,…,z Mw ,…,z MB ] T The initial sound field ψ0 is calculated using the image source method described later. w is the index indicating the boundary position between the water and the seabed, M Bis the index indicating the bottom of the seafloor, and N is the index indicating the bottom of the calculation range of the sound field.

[0031] The image source method is a method for calculating sound pressure by superimposing waves emitted from a real sound source and a mirror image sound source. Figure 5 is an image of the calculation of the initial sound field using the image source method. The solid line in Figure 5 indicates the direct wave, and the dashed line indicates the reflected wave. Assuming that the sound wave travels in a straight line, the sound pressure created by the sound source is given by spherical diffusion (1 / l) x phase rotation (exp(iωl / c)). Here, l is the path length, with l1 being the path of the direct wave and l2 being the path of the reflected wave, ω is the angular frequency, and c is the average speed of sound from the sound source depth to the receiver depth. Considering the superposition of the real sound source and the mirror image sound source, the initial sound field ψ0 = [ψ(r s, Z[1]),…,ψ(r s, Z[N])] T i is an element of z The depth Z[i z ] distance r s The sound field ψ(r s, Z[i z ]) is the underwater region (Z[i z ]≦z bottom, i z =1,…,M w ) is calculated as in equation (12), where i z is the computation depth index.

number

[0032] The path l1 of the direct wave and the path l2 of the reflected wave are obtained from the following equations (13) and (14), respectively.

number

number

[0033] The average sound velocity c1 of the direct wave and the average sound velocity c2 of the reflected wave are calculated from the following equations (15) and (16), respectively.

number

number

[0034] where z s is the source depth, r s is the initial sound field calculation distance, and s is the sound source depth index. The initial sound field calculation unit 51 calculates the initial sound field distance for each element i in the underwater region of the depth series Z as described above. z =1,…,M w By applying equation (12) to the initial sound field ψ0, we can obtain the initial sound field ψ0. Also, when the receiving depth is below the seafloor (Z[i z ]>z bottom ,i z =M w +1,…,N) is multiplied by the attenuation that occurs in the sedimentary layer depending on the depth difference between the seafloor and the receiving depth, so the sound pressure is calculated using the following equation (17) instead of equation (12).

number

[0035] z in equation (17) diff is i z The depth Z[i z ] and seafloor depth z bottom and is shown in equation (18).

number

[0036] The initial sound field calculation unit 51 calculates the initial sound field ψ0 and updates the pre-update sound field ψ r0 and stores the calculated distance r in the storage unit 56 as the pre-update distance r0. s is stored in the storage unit 56. The sound field before updating ψ r0 is the sound field ψ at the distance r0 before updating r0 is.

[0037] The approximation coefficient calculation unit 52 calculates the approximation coefficient a using the approximation order m and the distance step Δr. j,m ,b j,m (j=1,...m) are calculated. The approximation coefficient a j,m ,b j,m The calculation method for (j=1, ...m) is explained below. Pade approximation is an approximation method that approximates a function with a rational function. When the numerator and denominator are approximated with a rational function whose denominator and numerator are m-th degree polynomials, the exponential function term in equation (4) is approximated as in equation (19).

number

[0038] Here, if we use the condition that the right side of equation (19) coincides with the 2m-th order Taylor expansion of the left side of equation (19), we obtain equation (20).

number

[0039] where coefficients c1,…,c 2m is the Taylor expansion up to order 2m at q=0 of the function to be approximated. The denominator (1+β1q+…+β m q m ) on both sides to obtain equation (21).

number

[0040] The simultaneous equations in equation (22) can be obtained from the condition that the coefficients of the 1st to 2mth order terms of q on both sides are equal.

number

[0041] By solving equation (22), the coefficients α1…α m ,β1…β m The approximation coefficient calculation unit 52 calculates the coefficients α1...α m ,β1…β mOnce this is determined, the right-hand side of equation (19) can be transformed into the right-hand side of equation (6) by partial fraction decomposition, and the approximation coefficient a j,m ,b j,m Calculate (j=1,...m).

[0042] The sound pressure calculation unit 53 calculates the pre-update sound field ψ r0 and the approximation coefficient a j,m ,b j,m (j=1,...m) and the environmental conditions (sound speed profile c profile , density profile ρ profile , the decay rate profile γ profile ) and the depth sequence Z, the sound field ψ at a distance of Δr r0+Δr Here, the distance step Δr and the approximation order m are calculated based on the allowable phase error condition in Equation (11). cal For example, when calculating a sound field in a distance-independent environment where the sound speed profile in the sea and seabed, the density and attenuation rate of the seabed, and the seabed topography do not change with distance, and furthermore, when calculating a sound field in a constant sound speed field where the sound speed is constant at depth, the sound ray path propagating at a high depression angle will repeatedly undergo multiple reflections with the seabed and the sea surface, and at sufficiently long distances, only components with low depression angles will remain. If the purpose is to calculate a sound field at a long distance, the calculation angle θ will be determined from the sound ray path that reaches that distance with the fewest reflections with the seabed. cal can be set.

[0043] Specifically, the sound pressure calculation unit 53 calculates the sound pressure based on the environmental conditions (sound speed profile c profile , density profile ρ profile , the decay rate profile γ profile ), approximation coefficient a j,m ,b j,m , sound field ψ before update r0 By solving the differential equation (7) using r0+Δr Calculate the vertical sound field ψ in equation (23). r0+Δr The differential equation for solving is shown below.

number

[0044] As shown in equation (24), the jth approximation coefficient a j,m ,b j,m The solution to φ j Then, equation (23) can be expressed as equation (25).

number

number

[0045] On both sides of equation (24) j,m q) and substituting equation (5) for q gives equation (26).

number

[0046] By calculating Equation (26) using the differencing described in Non-Patent Document 2, the j-th approximation coefficient a of the approximation order m is obtained. j,m and b j,m Solution φ for j In the conventional technology, the environmental conditions (sound speed profile c profile , density profile ρ profile , the decay rate profile γ profile ) are used, but they are not directly related to the present invention, so their explanation will be omitted. The sound pressure calculation unit 53 calculates all approximation coefficients a up to the approximation order m. j,m ,b j,m For (j=1,…m) φ j After calculating the depth Z[i z ] sound pressure ψ(r0+Δr,Z[i z ]) is calculated for all elements of the depth sequence Z. By calculating this for all elements of the depth sequence Z, the vertical sound field ψ r0+Δr =[ψ(r0+Δr,Z[1]),…,ψ(r0+Δr,Z[N])] T is obtained.

[0047] The update processing unit 54 updates the sound field before update ψ r0 And the pre-update distance r0 is updated as shown in equations (27) and (28).

number

number

[0048] The determination unit 55 determines the pre-update distance r0 and the calculated maximum distance R max Specifically, the determination unit 55 determines whether to end the sound wave propagation calculation process by determining whether the pre-update distance r0 is equal to or greater than the maximum distance R max In the following case (r0≦R max ) continues the sound pressure calculation in the sound pressure calculation unit 53, and the pre-update distance r0 is the maximum distance R max If it exceeds (r0>R max ) ends the sound wave propagation calculation process as it has satisfied the calculation distance range.

[0049] 6 is a flowchart of the sound wave propagation calculation process in the prior art. First, the initial sound field is calculated by the initial sound field calculation unit 51 (S101). The initial sound field calculation unit 51 calculates the initial sound field based on the initial sound field calculation distance r s , source depth z s , seabed depth z bottom , sound speed profile c profile , the decay rate profile γ profile , and the depth sequence Z is used to calculate the initial sound field ψ0, and the sound field ψ at the pre-update distance r0 is calculated. r0 and stores it in the storage unit 56 as

[0050] Next, the approximation coefficient calculation unit 52 calculates the approximation coefficients (S102). The approximation coefficient calculation unit 52 calculates the approximation coefficients a of the Pade approximation using the distance step Δr and the approximation order m set in advance. j,m ,b j,m Calculate (j=1,...m).

[0051] Next, the sound pressure calculation unit 53 calculates the sound pressure (S103). The sound pressure calculation unit 53 calculates the sound pressure before updating the sound field ψr0 and the approximation coefficient a j,m ,b j,m (j=1,...m) and the environmental conditions (sound speed profile c profile , density profile ρ profile , the decay rate profile γ profile ) and the depth sequence Z, the sound field ψ at a distance of Δr r0+Δr Calculate.

[0052] Next, the update processing unit 54 performs update processing (S104). The update processing unit 54 uses the sound field calculated by the sound pressure calculation unit 53 and the distance step Δr to calculate the pre-update sound field ψ stored in the storage unit 56. r0 , and the pre-update distance r0 is updated.

[0053] Thereafter, the determination unit 55 determines whether or not to end the sound wave propagation calculation process (S105). max If the sound wave propagation calculation process is not to be ended (S105: NO), the process returns to step S103, and the pre-update distance r0 is updated to the maximum distance R max The processes of steps S103 and S104 are repeated until the calculated value exceeds 1. On the other hand, if the sound wave propagation calculation process is to be ended (S105: YES), the sound wave propagation calculation process is ended.

[0054] FIG. 7 is a physical image diagram of the calculation conditions in the sound wave propagation calculation process of the prior art. In the above explanation, the sound speed profile c profile , the decay rate profile γ profile , density profile ρ profile The explanation is based on the assumption that each has a value corresponding to the depth series Z. However, each environmental condition (c profile , γ profile , ρ profile ) only has information at different depths, each calculation depth Z[i zWhen the values ​​of these environmental conditions corresponding to [ ] are required, the value at the relevant depth can be obtained by selecting the closest depth value from the given environmental condition vector each time, or by interpolating from the values ​​corresponding to the depths before and after that, so there is no loss of generality even if this assumption is made. Note that the creation of environmental conditions is an example of creating calculation conditions, and does not have to be processed within the device.

[0055] The above-described conventional sound propagation calculation device 500 has the following problems. Figure 8 is a diagram illustrating the difference in sound propagation distance depending on the depression angle. In Figure 8, the propagation of sound waves at low depression angles is shown by a dashed line, and the propagation of sound waves at high depression angles is shown by a solid line. As described above, in the case of underwater sound propagation in a distance-independent environment where the underwater and seafloor sound speed profile, seafloor density and attenuation rate, and seafloor topography do not change with distance, and in a constant sound speed environment where the speed of sound does not depend on depth, the components at high depression angles tend to contribute less to propagation as the propagation distance increases, and the components at low depression angles become dominant. This is because, as shown in Figure 8, almost all of the energy of a sound wave reflected at the sea surface returns into the water, whereas at the seafloor, some of the energy propagates into the seafloor and only a portion of the energy returns into the sea. Therefore, the more the seafloor reflections, the greater the attenuation of the sound wave energy. Therefore, at low depression angles, there is little reflection and the signal propagates over long distances, while at high depression angles, there is much reflection and the signal does not propagate over long distances. In other words, the range of elevation and depression angles that contains the dominant energy in propagation can be thought of as narrowing depending on the propagation distance (hereinafter, the range of elevation and depression angles that contains the dominant energy in propagation will be referred to as the "propagation angle").

[0056] FIG. 9 is a diagram showing the characteristics of the propagation angle versus the propagation distance. FIG. 9 is a diagram when the propagation angle is defined as an angle including the energy ratio (A%) that is desired to be secured in the propagation calculation. θa in FIG. 9 is the calculated angle when calculation is attempted with sufficient accuracy in a short distance, and θb is the calculated angle when calculation is attempted with sufficient accuracy in a long distance. Furthermore, when calculation is performed with θa as the calculation angle, region A1 is a region where the propagation angle is smaller than the calculation angle, and calculation is wasted. On the other hand, when calculation is performed with θb as the calculation angle, region B1 is a region where the propagation angle is larger than the calculation angle, and energy loss occurs. The conventional PE method calculates the calculation angle θ from the propagation angle at the propagation distance for which the sound field is desired according to the characteristics shown in FIG. cal and set the tolerance ε based on Eq. (11). max Taking this into consideration, the approximation order m and distance step Δr are set in advance.

[0057] However, even in the case of a distance-independent uniform sound velocity field, if you want to calculate a wide range of sound fields from close to long distances, you need to calculate the calculation angle θ cal If the angle is narrowed, the accuracy of the calculation of sound wave propagation will be deteriorated because the high depression angle component at close range is not taken into account. Therefore, the calculation angle θ cal If the angle is set sufficiently wide, there is a problem that the calculation load increases due to the calculation of components with high depression angles that do not need to be calculated at long distances.

[0058] Furthermore, in a distance-dependent environment where the sound speed profile or seabed topography changes with distance, the propagation angle narrows globally with distance, but the narrowing characteristic of the propagation angle varies with distance, making it difficult to estimate the propagation angle at the desired propagation distance.

[0059] Therefore, in the present invention, attention is paid to the phenomenon that the propagation angle narrows as the propagation distance increases globally, and by using a calculation angle that differs depending on the distance, the calculation angle θ cal This suppresses excess or deficiency in the calculation, and efficiently calculates sound wave propagation in terms of calculation accuracy and amount of calculation.

[0060] Embodiment 1 A sound wave propagation calculation device 100 according to the first embodiment will be described. FIG. 10 is a schematic diagram of the sound wave propagation calculation device 100 according to the first embodiment. As shown in FIG. 10, the sound wave propagation calculation device 100 according to the present embodiment includes an initial sound field calculation unit 1, an initialization unit 2, an optimal parameter calculation unit 3, an approximation coefficient calculation unit 4, a sound pressure calculation unit 5, an update processing unit 6, a determination unit 7, and a storage unit 8. The initial sound field calculation unit 1, the initialization unit 2, the optimal parameter calculation unit 3, the approximation coefficient calculation unit 4, the sound pressure calculation unit 5, the update processing unit 6, and the determination unit 7 of the sound wave propagation calculation device 100 are functional units realized by dedicated equipment or dedicated processing circuits. Alternatively, the sound wave propagation calculation device 100 may include a processor such as a CPU, and the processor may execute a sound wave propagation calculation program stored in the storage unit 8 to realize each functional unit of the sound wave propagation calculation device 100. Alternatively, each functional unit of the sound wave propagation calculation device 100 may be realized by a combination of dedicated equipment or dedicated circuitry and software.

[0061] The storage unit 8 is composed of a volatile memory such as a RAM, a DRAM, or an SRAM, a non-volatile memory such as a ROM, a flash memory, or a hard disk drive, or a combination thereof. The storage unit 8 stores a sound wave propagation calculation program executed by the processor of the sound wave propagation calculation device 100, various data such as parameters required for executing the sound wave propagation calculation program, etc.

[0062] The functions of the initial sound field calculation unit 1, the approximation coefficient calculation unit 4, the sound pressure calculation unit 5, and the update processing unit 6 are the same as the functions of the initial sound field calculation unit 51, the approximation coefficient calculation unit 52, the sound pressure calculation unit 53, and the update processing unit 54 in the prior art, so their explanation will be omitted.

[0063] The initialization unit 2 calculates the parameter update distance R end 0 and parameter update count n update are initialized as shown in equations (29) and (30), respectively. end 0 The initialization is set so that optimization is always performed for the first loop, and can be set to any value as long as this purpose is achieved.

number

number

[0064] The determination unit 7 performs the same end determination of the sound wave propagation calculation process and the optimization determination as in the prior art. In the optimization determination, the determination unit 7 determines the parameter update distance R end nupdate and the pre-update distance r0, it is determined whether to optimize the calculation angle, i.e., whether to change the distance step Δr and the approximation order m. Here, R end Superscript n update indicates the number of updates of the calculated angle. Specifically, the determination unit 7 determines whether the pre-update distance r0 is equal to or smaller than the parameter update distance R end nupdate If r0≦R end nupdate ), it is determined that optimization is not performed, and the pre-update distance r0 is the parameter update distance R end nupdate If greater than (r0>R end nupdate ) is determined to be optimized.

[0065] The optimal parameter calculation unit 3 calculates the pre-update sound field ψ r0 , and the distance before update r0 are used to calculate the calculated angle, and the distance step Δr nupdate and the approximation order m nupdate The optimal parameter calculation unit 3 calculates the optimal parameters including the parameter update distance R endnupdate The method for calculating the optimum parameters by the optimum parameter calculation unit 3 will be described in detail later.

[0066] 11 is a flowchart of the sound wave propagation calculation process in the first embodiment. First, the initial sound field calculation unit 1 calculates an initial sound field (S1). The initial sound field calculation unit 1 calculates an initial sound field based on an initial sound field calculation distance r s , source depth z s , seabed depth z bottom , sound speed profile c profile , the decay rate profile γ profile , and the depth sequence Z is used to calculate the initial sound field ψ0, and the pre-update sound field ψ at the pre-update distance r0 r0 and stores it in the storage unit 8 as

[0067] Subsequently, initialization is performed by the initialization unit 2 (S2). The initialization unit 2 calculates the parameter update distance R end nupdate and parameter update count n update Then, the determination unit 7 performs optimization determination (S3). Here, the determination unit 7 initializes the parameter update distance R end nupdate Compare the distance before update r0 and calculate the angle θ cal It is determined whether to optimize. Note that for the first loop (first step from the initial sound field), the optimization decision always determines that optimization should be performed.

[0068] If optimization is not to be performed (S3: NO), the process proceeds to step S6. On the other hand, if optimization is to be performed (S3: YES), the optimal parameter calculation unit 3 performs optimal parameter calculation processing (S4). In the optimal parameter calculation processing, the pre-update sound field ψ r0 and the pre-update distance r0 are input, and the optimal calculation angle, the distance step Δr and approximation order m to achieve that calculation angle, and the parameter update distance R, which is the distance until the next calculation angle optimization is required, are calculated. end nupdate is calculated.

[0069] FIG. 12 is a flowchart showing the flow of the optimum parameter calculation process in the first embodiment. The optimum parameter calculation process is performed by the optimum parameter calculation unit 3. First, a calculation angle for obtaining the optimum parameter is calculated (S41). The optimum parameter calculation unit 3 calculates the pre-update sound field ψ r0 , energy retention rate ε, and central wave number k0 are used to estimate the propagation angle, and the calculated angle θ cal Here, the calculation angle θ cal As one method for calculating ψ, we will explain how to calculate it from the wave number characteristics of the sound field. r0 The wavenumber characteristic Ψ p =[Ψ0,…,Ψ ζ / 2 Specifically, the pre-update sound field ψ r0 is converted into the wavenumber spectrum Ψ[s] by Fourier transform as shown in equation (31).

number

[0070] where N is the maximum depth index of the initial sound field ψ0, s is the vertical wave number index (s takes an integer value from -N / 2 to N / 2), and i z indicates the depth index. The vertical wavenumber index s can be associated with the elevation angle at which the sound waves propagate, with s<0 corresponding to an upward propagation path, s>0 corresponding to a downward propagation path, and s=0 corresponding to a horizontal propagation path. Since the PE method does not allow for different angles to be set for upward and downward propagation, the angle range is set as the absolute value of the angle. Therefore, the wavenumber spectrum Ψ[s] in equation (31) can be converted to the wavenumber characteristic Ψ in equation (32). p Convert to [s].

number

[0071] Next, the wavenumber characteristic Ψ p and the energy retention rate ε are used to calculate the wavenumber k cal Specifically, as shown in equation (33), the wavenumber characteristic Ψ in the total wavenumber space is calculated. pThe wave number characteristic Ψ for the total energy of p The minimum vertical wave number index s where the ratio of the total energy value from s=0 to a specific vertical wave number index s exceeds the energy retention rate ε min Ask for.

number

[0072] Next, the vertical wavenumber index s min Calculated wave number k corresponding to cal Calculate the wave number k cal can be obtained from the vertical wave number sequence K as shown in equation (34). The vertical wave number sequence K is expressed by equation (35).

number

number

[0073] where Δz is the difference between elements of the depth sequence Z. Finally, using equation (36), we calculate the wavenumber k cal Calculate the angle θ cal Convert to.

number

[0074] Here, Fourier transform is used to obtain the wave number characteristics from a sound field with equal depth intervals, but any technique can be used, even if the depth intervals are not equal, or even if Fourier transform is not used, as long as the minimum wave number that satisfies equation (33) can be obtained. For example, it is possible to consider the vertical sound field as sound waves input to an equally spaced vertical array, and obtain the calculated angle from the phased output of the vertical array. In this case, the standby direction of the standby phase processing that forms independent beams in each elevation direction is set to θ beam [w], w=0,1,…,w max Let θ beam [w] is the horizontal direction, which is 0°, and its absolute value |θ beamThe phased output B[w] can be calculated using equation (37), and the w[w] satisfies equation (38). min The desired calculated angle can be calculated from w min is the smallest phase orientation index at which the ratio of the sum of the energy from w=0 to a specific phase orientation w to the total energy in the phased output of all phase orientations exceeds the energy retention rate ε.

number

number

[0075] where w is the phase orientation index, w max indicates the phasing index equivalent to either 90° upward or downward, sorted in ascending order. In this way, if the method can obtain almost independent energy for each depression angle, the calculated angle can be calculated by using any phasing method. cal =θ beam [w min ]. The beamforming processor in equation (37) can be of any type as long as it is a beamforming processor that forms independent beams in each elevation and depression direction.

[0076] Next, the optimum parameter calculation unit 3 acquires the approximation order and the distance step (S42). cal and angle table θ TBL (m, Δr) and the approximation order m used to calculate the sound pressure nupdate and distance step Δr nupdate The angle table is a table that shows the combination of ε(m, Δr, θ cal )=ε max (i.e., the phase error is within the allowable range) TBL(m, Δr) Fig. 13 is a conceptual diagram of an angle table used in the optimum parameter calculation process in the first embodiment (however, in Fig. 13, Δr1>Δr2> . . . Δr9).

[0077] The approximation order m can be determined without unnecessary calculations due to excessive calculation angles and while ensuring a predetermined energy retention rate. nupdate and distance step Δr nupdate There are multiple combinations of these. In Fig. 13, the combinations in the range enclosed by the thick lines are selectable candidates, and areas A2 and B2 are combinations that cannot be selected. Area A2 is a combination where the calculation angle is excessive and unnecessary calculations are performed, and area B2 is a combination where the calculation accuracy is insufficient. The range enclosed by the thick lines solves the problems of the conventional technology in a distance-independent environment, and the approximation order m is calculated as shown in equation (39). nupdate and distance step Δr nupdate Any combination can be adopted from the combinations obtained.

number

[0078] However, combinations that do not satisfy the condition judgment (empty set condition) are removed from the combinations. From these multiple candidate combinations, update Approximation order and distance step (m nupdate ,Δr nupdate )=(m select ,Δr select ) is selected. For example, if reducing the processing load is important, the largest degree and largest distance increment among the selection candidates can be selected. Also, if this technology is applied to conditions where the environment changes gradually with distance, it is possible to ensure adaptability to gradual changes in the distance direction by selecting the smallest degree and smallest distance increment at the expense of a slight increase in the amount of calculation.

[0079] Next, the optimal parameter calculation unit 3 calculates the parameter update distance (S43).cal Optimize (i.e., the approximation order m nupdate and distance step Δr nupdate (changing the parameter update distance R end nupdate Calculate the parameter update distance R end nupdate Three methods for setting the distance interval R are explained below. int is set as a fixed value, and the parameter update distance R end nupdate This is a method of setting the initial sound field calculation distance r s , parameter update count n update Distance interval R according to int Since it is sufficient to add end nupdate is obtained by equation (40).

number

[0080] Although this method is simple, the distance interval R int If m is set too short, the frequency of calculating the optimal parameters increases, and even though the propagation angle is not narrow, the approximation order and distance step calculation are the same (m select ,Δr select ) combination, which may result in unnecessary calculation of optimal parameters and approximation coefficients. select ,Δr select ) is selected, some countermeasures can be taken, such as not performing calculations by the approximation coefficient calculation unit 4 and using the previous value. int If the angle is set too long, sound pressure will be calculated using the same calculation angle even though the propagation angle has narrowed, resulting in unnecessary calculations.

[0081] The second method is a countermeasure against the above-mentioned problem of the first method. The larger the change in the calculated angle, the smaller the parameter update distance R end nupdateThe smaller the change in the calculated angle, the smaller the parameter update distance R end nupdate Specifically, the previous calculated angle θ cal nupdate-1 and the calculated θ cal nupdate If equation (41) is satisfied (i.e., the angle change is large), R int If equation (42) is satisfied (i.e., the angle change is small), R int Increase the

number

number

[0082] where θ th is the threshold value for evaluating the magnitude of the calculated angle change. int After increasing or decreasing the parameter update distance R end nupdate By setting the calculated angle θ cal If the change in R is large, the parameter update distance can be set finely, and if the change is small, the parameter update distance can be set coarsely. int Regarding the method of increasing or decreasing , the method shown in equation (43) can be considered as an example, where α is an arbitrary constant.

number

[0083] The third method is to calculate the parameter update distance by predicting the distance at which the propagation angle becomes smaller than the calculated angle. More specifically, by predicting how the propagation angle narrows, the parameter update distance R is calculated finer for shorter distances and coarser for longer distances. end nupdate The propagation angle has the characteristic that it changes rapidly as the angle increases and gradually as the angle decreases. If it changes rapidly, the parameter update distance R is set at short intervals to accommodate the change. endnupdate In the case of slow changes, it is necessary to set the parameter update distance R end nupdate In addition, considering that the angle is larger for shorter distances and smaller for longer distances, the parameter update distance R can be set finer for shorter distances and coarser for longer distances. end nupdate There are various methods for predicting the propagation angle, but one example is a method for predicting the distance by assuming that the propagation angle exists on a curve inversely proportional to the calculated distance as shown below. First, the calculated angle θ cal , the minimum calculated angle θ min (minimum calculation angle that can be set), distance before update r0, maximum calculation distance R max Then, (x, y)=(r 0, θ cal ),(R max, θ min ) is obtained from equation (44).

number

[0084] α in equation (44) is obtained from equation (45), and K is obtained from equation (46).

number

number

[0085] Next, calculate the angle θ cal From the angle change amount θ del θ^ subtracted by cal (=θ cal -θ del ) and substitute it into equation (44), the parameter update distance R end nupdate The symbol "^" written next to θ is a hat indicating an estimated value, and is written above θ as shown in equation (47).

number

[0086] Finally, the minimum calculated angle θ min The method for setting the minimum calculated angle θ will be explained. min is set from the maximum critical angle obtained from the sound speed profile in the ocean, and is calculated from equation (48).

number

[0087] where c min is the sound speed profile c profile The minimum value in c max is the sound speed profile c profile As mentioned above, the parameter update distance R end nupdate By setting the parameter update distance, it is possible to set the parameter update distance at fine intervals in the short distance where the calculation angle is large, and at coarse intervals in the long distance where the calculation angle is small.

[0088] Next, the optimal parameter calculation unit 3 counts the number of updates (S44). Specifically, the optimal parameter calculation unit 3 counts the number of parameter updates n update is updated as shown in equation (49). After that, the optimum parameter calculation process is terminated, and the process proceeds to step S5 in FIG.

number

[0089] When the optimum parameters are calculated by the optimum parameter calculation unit 3, the approximation coefficient calculation unit 4 calculates the approximation coefficients a of the Pade approximation using the distance step Δr and the approximation order m calculated by the optimum parameter calculation unit 3 (S5). j,m ,b j,m Calculate (j=1,...m).

[0090] Thereafter, the sound pressure is calculated by the sound pressure calculation unit 5 (S6). The sound pressure calculation unit 5 calculates the sound pressure before updating ψ r0 and the approximation coefficient a j,m ,b j,m and the environmental conditions (sound speed profile c profile , density profile ρ profile , the decay rate profile γ profile ) and the depth sequence Z, the sound field ψ at a distance of Δr r0+Δr Calculate.

[0091] Next, the update processing unit 6 performs update processing (S7). The update processing unit 6 uses the sound field calculated by the sound pressure calculation unit 5 and the distance step Δr to calculate the sound field before update ψ r0 , and the pre-update distance r0 is updated. Thereafter, the determination unit 7 determines whether the process is to be completed (S8), and if so (S8: YES), the sound wave propagation calculation process is terminated, and if not (S8: NO), the process returns to step S3 and the subsequent processes are repeated.

[0092] FIG. 14 is an image diagram of the sound wave propagation calculation process in the sound wave propagation calculation device 100 according to the first embodiment. FIG. 14(a) shows an image of resetting the approximation order m and the distance increment Δr, and FIG. 14(b) shows an image of the calculation angle following the propagation angle. Note that the superscripts of the variables in FIG. 14 indicate the number of calculation angle updates. As described above, the sound wave propagation calculation device 100 according to this embodiment calculates the intermediate calculation distance (the distance R in FIG. 13) and the distance increment Δr. end 1 , R end 2 ), the approximation order m and distance step Δr (optimal parameters) are reset by reviewing the calculation angle. This allows the calculation angle to become smaller in stages, as shown in Figure 14(b), and the calculation angle can be optimized according to the distance. As a result, the vertical sound field can be calculated using a calculation angle that corresponds to the propagation angle, which narrows as the sound propagates further away, and it is possible to prevent the calculation angle from being too large or too small compared to the propagation angle.

[0093] As described above, according to the sound wave propagation calculation device 100 of this embodiment, a calculation angle corresponding to the propagation angle that narrows as the wave propagates to a long distance is calculated, and optimal parameters including the distance increment and the approximation order are determined to perform sound wave propagation calculations. In a distance-independent environment, sound wave propagation calculations can be performed efficiently at a high depression angle in the short distance to avoid losing energy, and at a low depression angle in the long distance to avoid wasting calculations in areas where no energy exists. In addition, in a distance-dependent environment, calculations are performed while sequentially resetting appropriate calculation angles according to the propagation conditions, so that sound wave propagation calculations can be performed at optimal calculation angles even if the propagation angle cannot be predicted in advance. From the above, by using the sound wave propagation calculation device 100 of this embodiment, the calculation angle θ for the propagation angle can be calculated in both a distance-independent environment and a distance-dependent environment. cal This makes it possible to suppress excess or deficiency in the calculation, and it becomes possible to efficiently perform sound wave propagation calculations in terms of calculation accuracy and calculation amount.

[0094] The above is a description of the embodiments of the present invention, but the present invention is not limited to the configurations of the above embodiments, and various modifications or combinations are possible within the scope of the technical concept. For example, in the first embodiment, when calculating the parameter update distance, the optimal parameter calculation unit 3 calculates the parameter update distance R by predicting the narrowing of the propagation angle. end nupdate In this case, the prediction method introduced is a linear prediction method based on the assumption that the propagation angle is inversely proportional to the distance, but other methods can also be used to predict the narrowing of the propagation angle.

[0095] For example, the calculated angle [θ cal 0 ,θ cal 1 ,…θ cal nupdate ] and its updated distance [r0,R end 1 ,…,R end nupdate] can be used to predict the narrowing of the propagation angle by polynomial fitting. Any method can be used as long as it can predict the narrowing of the propagation angle.

[0096] Furthermore, in the initial sound field calculation unit 1 of embodiment 1, an initial sound field generation method using the mirror sound source method of the prior art has been described as an example, but the initial sound field may be generated using any method as long as the initial sound field can be obtained.

[0097] Furthermore, in the first embodiment, the optimal parameter calculation unit 3 is explained as evaluating the magnitude of the angle change based on the difference between the calculated angles when calculating the parameter update distance. However, the ratio of the calculated angles (θ cal nupdate / θ cal nupdate-1 ) to evaluate the magnitude of the change. Any method can be used as long as it can measure the change in the calculated angle.

[0098] In the first embodiment, the optimal parameter calculation unit 3 calculates the parameter update distance by using the distance interval R int Although the method of multiplying or dividing by a constant ratio α has been described, the constant ratio α may be changed by multiplication or division, or the threshold θ th θ th1 and θ th2 Two types are provided, and the range of change is <θ th1 , θ th1 ≦Change range≦θ th2 , θ th2 <Change width, distance interval R int Increasing the distance between R int Maintain a distance R int A method of reducing the

[0099] Furthermore, in the first embodiment, the distance interval R int As a method of increasing or decreasing R, we have explained that it is multiplication or division by a fixed ratio. int =min(R int +ΔR i + ,R int max), R int =max(R int min ,R int -ΔR i - ) where R int max and R int min are the distance intervals R int Upper and lower limits of ΔR i + and ΔR i - are the ranges of addition and subtraction, respectively. Calculated angle θ cal nupdate Depending on the result of the change in the distance R int If you want to increase or decrease it, you can use any method of increase or decrease. [Explanation of symbols]

[0100] 1 initial sound field calculation unit, 2 initialization unit, 3 optimal parameter calculation unit, 4 approximation coefficient calculation unit, 5 sound pressure calculation unit, 6 update processing unit, 7 judgment unit, 8 memory unit, 51 initial sound field calculation unit, 52 approximation coefficient calculation unit, 53 sound pressure calculation unit, 54 update processing unit, 55 judgment unit, 56 memory unit, 100 sound wave propagation calculation device, 500 sound wave propagation calculation device.

Claims

1. A sound wave propagation calculation method for calculating sound pressure by a Parabolic Equation (PE) method using the Split-Step Pade method, an optimal parameter calculation step for calculating optimal parameters including an approximation order and distance step that vary depending on the calculation distance; an approximation coefficient calculation step of calculating approximation coefficients using the optimal parameters; a sound pressure calculation step of calculating the sound pressure using the approximation coefficient.

2. The optimal parameter calculation step includes:

2. The sound wave propagation calculation method according to claim 1, further comprising the step of calculating the optimum parameters according to a propagation angle that is a range of elevation and depression angles that includes dominant energy in the calculation distance.

3. The optimal parameter calculation step includes: Calculating a calculated angle from a vertical sound field including the sound pressure at a calculation distance one step before; The sound wave propagation calculation method according to claim 2 , further comprising the step of: obtaining an optimal combination of the approximation order and the distance step from the calculation angle.

4. The step of calculating the calculated angle includes: A step of Fourier transforming the vertical sound field to calculate wave number characteristics; 4. The sound wave propagation calculation method according to claim 3, further comprising the step of: calculating the calculation angle based on a ratio of the energy of the wave number characteristic in a specific wave number space to the energy of the wave number characteristic in the entire wave number space.

5. The step of calculating the calculated angle includes: calculating a phased output when the vertical sound field is used as an input for a vertical array; 4. The sound wave propagation calculation method according to claim 3, further comprising the step of: calculating the calculated angle based on a ratio of the energy of the phased output in a specific phased azimuth to the energy of the phased output in all phased azimuths.

6. The step of obtaining the combination includes: obtaining the combination using the calculated angle and an angle table; 4. The sound wave propagation calculation method according to claim 3, wherein the angle table contains preset combinations of the approximation orders and the distance increments that keep the phase error in the calculated angles within an allowable range.

7. The optimal parameter calculation step includes:

4. The sound wave propagation calculation method according to claim 3, further comprising a step of calculating a parameter update distance, which is a calculation distance for performing said optimum parameter calculation step.

8. The step of calculating the parameter update distance includes: calculating the parameter update distance by adding a distance interval to the previous parameter update distance; The sound wave propagation calculation method according to claim 7 , wherein the distance interval is calculated based on a preset fixed value.

9. The step of calculating the parameter update distance includes: calculating the parameter update distance by adding a distance interval to the previous parameter update distance; The sound wave propagation calculation method according to claim 7 , wherein the distance interval is a value that increases or decreases according to an angle change that is a difference between the current calculated angle and the previous calculated angle.

10. The step of calculating the parameter update distance includes: The sound wave propagation calculation method according to claim 7 , further comprising the step of calculating the parameter update distance by predicting a distance at which the propagation angle becomes smaller than the calculated angle.

11. The step of calculating the parameter update distance includes: The sound wave propagation calculation method according to claim 10, further comprising the step of predicting the distance by assuming that the propagation angle lies on a curve inversely proportional to the calculated distance.

12. A sound wave propagation calculation program that causes a processor of a sound wave propagation calculation device to execute the sound wave propagation calculation method according to any one of claims 1 to 11.

13. A sound wave propagation calculation device that calculates sound pressure by a PE (Parabolic Equation) method using the Split-Step Pade method, an optimal parameter calculation unit that calculates optimal parameters including an approximation order and distance increment that vary depending on the calculation distance; an approximation coefficient calculation unit that calculates approximation coefficients using the optimal parameters; a sound pressure calculation unit that calculates the sound pressure using the approximation coefficient.

Images