JPH0129429B2 - - Google Patents

Description

Detailed Description of the Invention (Technical Field) The present invention relates to a sound wave propagation simulation method for analyzing the propagation characteristics of sound waves in water. This method relates to a simulation method for determining the received wave output when a sound wave is received by a wave transmitter.In particular, when the sound source or receiver has directivity, it is used to accurately represent the effect of the directivity. This provides a simulation method that allows for

(Background Art) The performance of marine acoustic equipment such as sonar varies greatly depending on the marine environment at the site of its operation. Therefore,
For effective operation of equipment, it is important to understand the characteristics of sound wave propagation in the ocean, and various sound wave propagation models have been developed.

Normal mode theory is a method for providing exact solutions to wave equations. According to this, sound waves radiated from a sound source are expressed as a superposition of cylindrical traveling waves (called normal mode waves) that form standing waves in the vertical direction.

This situation will be explained using the cylindrical coordinate system shown in FIG. A sound field in a stratified medium caused by a point source of unit intensity that vibrates harmonically at frequency f is expressed by the following Helmholtz equation.

▽ ² φ(r→)+K ² (r→)φ(r→)=−4πδ(r
→−r→ _S )
...(1) Here, ▽ ² is the Laplace operator, r→ is the coordinate of the receiving point, r→ _S = (O, O, Z _s ) is the coordinate of the sound source, φ(r→)
exp(-i2πft) is sound pressure, K(r→)=2πf/C(r)
is the wave number of the sound wave in the medium, and δ is the Dirac delta function. Hereafter, for simplicity, φ, which is the complex amplitude of sound pressure,
Let us simply call (r→) the sound pressure.

Here, assuming that K(r→) is a function of depth only, K(z _R ), equation (1) becomes symmetrical with respect to the Z axis of the solution.
Therefore, from now on, we will consider it on the plane of θ=0,
The notation of θ is omitted. The solution to equation (1) is as a superposition of sound waves spreading in all directions, φ (r, z) = 2∫ ^∞ ₀ G (z, z _S ; K) J _O (Kr) KdK
...It is expressed as (2). Here, Jo is a zero-order Betzel function. Also, G is Green's function, and the following Z
Equation (3) for , and the boundary conditions (sound pressure is zero at sea surface and infinite depth, sound pressure and vertical particle velocity are continuous at the layer interface) are satisfied. Further, K is a complex number on a complex plane.

^d2G / ^dz2 +[ ^K2 (z) ^−K2 ]G=−δ(z− _zs )
...(3) The normal mode solution is obtained by applying the residue theorem to equation (2), and is expressed by the cumulative sum of normal mode waves and some branch line integrals as follows.

φ (r, z) = πi 〓 〓 ⁿ H _O ⁽¹⁾ (K _o r) U _o (z) U _o (z _S ) + ∫ _BL G (z, z _S ; K
)H _O ⁽¹⁾ (Kr)KdK(4) Here, n is the mode number, and H _O ⁽¹⁾ represents the 0th order Wankel function of the first kind. Further, K _o is an eigenvalue of the n-th mode, and is a singular point on the complex plane of the Green's function G. U _o (z) is the normal mode eigenfunction in the normalized depth direction, and the homogeneous expression of equation (3) d ² U _o (z)/dz ² + [K ² (z) − K ² _o ] U _o (z)=0
(5) and the boundary conditions are satisfied. Also, the following orthogonality condition is satisfied ∫ ^∞ ₀ [ρ(z)/ρ(Zs)U _o (z)U _n (z)dz=δ
_on
...(6) Here, ρ(z) is the density of the medium, and δ _on represents the Kronecka delta function.

The second term on the right side of equation (4) represents the integral of the complex number K on the complex plane along the branch line of G. A branch line occurs when G becomes a multivalued function with respect to K and has a branch point in a half plane where the real part of K is positive. Therefore, we divide the ocean, including the ocean floor, into arbitrary layers, and express the square of the refractive index of the sound wave within each layer as a linear function with respect to depth (approximate the sound speed profile with a polygonal line).
Then, G becomes a single-valued function with respect to K, the integral term in equation (4) disappears, and a complete sound field is expressed only by the discrete normal mode waves of equation (7) below.

P(r, Z _R , Z _S ; f) iπ _∞ 〓 ⁿ⁼¹ U _o (Z _S ) U _o (Z _R ) H _O ⁽¹⁾ (K _o r) ……(7) Sound source and receiver For the case where is directional,
Conventionally, this effect was introduced in the form of the following equation.

P(r, Z _R , Z _S ; f)=iπ _∞ 〓 〓 ⁿ⁼¹ D _S (θ _oS ) D _R (θ _oR ) U _o (Z _S ) U _o (Z _R ) H _O ⁽¹⁾ ( K _o r
)...(8) Here, D _S (θ) and D _R (θ) are the vertical directivity coefficients of the sound source and receiver, respectively, and θ is the angle of depression with the downward direction being positive with respect to the horizontal direction. . θ _oS and θ _oR represent depression angles at the sound source and receiving point depths of an equivalent sound ray for the n-th mode normal mode wave (a sound ray having the same turning point depth as the n-th mode normal mode wave). That is, if the wave number in the medium is expressed as k(Z) as a function of depth Z, θ _oS = cos ^-1 [K _o /k (Z _S )], θ _oR = cos ^-1 [K _o /k ( Z _R ) ] ... is given by (9).

In other words, what Equation (8) means is that when each normal mode wave is radiated from the sound source, the weight of the sound source directivity coefficient D _S (θ _oS ) corresponding to the depression angle θ _oS at the sound source depth of the equivalent sound ray The receiving directivity D _R (θ _oR ) corresponding to the depression angle θ _oR at the reception depth of the equivalent sound ray when the wave is radiated with the same angle and received at the receiving point.
This means that it is weighted and output.
This method is based on the assumption that each normal mode wave is radiated in a direction along the equivalent acoustic line and arrives at the receiving point from the direction along the equivalent acoustic line.

However, this method has the following problems.

(a) As is clear from equation (7), each normal mode wave has an amplitude distribution of U _o (z) in the vertical direction and propagates in the horizontal direction (θ = 0), as shown in Figure 3. It does not propagate in the direction along the equivalent sound line as shown in . The propagation of energy in the direction of the equivalent sound ray is caused by mutual interference of several normal mode waves. Therefore, the above assumption does not hold.

(b) In order for a sound source or receiver to have directivity, it must have a finite size, and directivity is defined when a thing of that finite size is placed in a homogeneous medium. However, since the medium in the actual ocean is inhomogeneous, the above-mentioned directivity cannot be strictly achieved and there will be errors. Conventional methods cannot evaluate this error.

(c) There are two equivalent sound rays of the normal mode wave radiated from the sound source: one radiated from the sound ray in the +θ _oS direction and one radiated in the −θ _{oS direction} .
Therefore, in order to express the effect of sound source directivity using equation (8), D _S (θ) must be symmetrical with respect to θ = 0. The same can be said for received wave directivity D _R (θ). . In the conventional method, therefore,
In cases where the main axis of the directivity pattern has an angle of depression, the effect of directivity that is asymmetric with respect to θ=0 cannot be evaluated.

(Problems to be solved by the invention) In order to eliminate these drawbacks, the present invention forms a vertically serially arrayed sound source and a vertically arrayed receiver having vertical directivity equivalent to the directional sound source and the directional receiver, respectively. Then, the sound field due to the directional sound source is calculated as a superposition of the sound field due to each element of the vertically arrayed sound source, and the received wave output of the directional receiver is calculated as a superposition of the sound field due to each element of the vertically arrayed sound source. This is calculated as the sum of the received wave outputs, and will be explained in detail below.

(Structure and operation of the invention) First, the concept of the invention will be explained. directionality of sound source
Assume that D _S (θ) is defined in a homogeneous infinite medium with the speed of sound C _O. Consider a Cartesian coordinate system x-ξ with the sound source position as the origin and the ξ axis pointing vertically downward. If the weight distribution of a vertical linear sound source with vertical directivity equivalent to the directional sound source is W _S (ξ), then D _S (θ)
The following relationship holds true between and W _S (ξ).

D _S (θ)=∫ ^∞ _-∞ W _S (ξ)exp(i2πf/C _O ξsinθ)d
ξ ...(10) Here, if we replace η=f/C _O sinθ ...(11), equation (10) becomes the form of inverse Fourier transform, and (ξ) is given by the following Fourier transform. .

W _S (ξ) =∫ ^∞ _-∞ D[θ(n)]exp(−i2πηξ)daη,…
…(12) θ(η)=sin ^-1 (C _O η/f) …(13) Approximating equation (12) by discrete Fourier transform, we get
A weight distribution of a vertically arranged sound source having vertical directivity equivalent to that of the directional sound source is obtained, and is expressed by the following equation.

W _Sl = 1/N _S exp (−i2πη _SO ξ _S1 ) _NS-1 〓 〓 ^m=0 E _Sn exp (−i2π1m/N _S ), l=0, 1, 2, ..., N
_S −1……(14) Here, E _Sn =D _Sn exp(−i2πmξ _SO Δη _S ), D _Sn =D _S [θ(η _Sn )]……(15) ξ _Sl =ξ _SO +lΔξ _S , η _Sn = η _SO +mΔη _S , Δξ _S・Δη _S = 1/N _S ...(16) where ξ _Sl and n _Sn are the sample values of ξ and n, respectively,
Δξ _S and Δη _S represent sample intervals of ξ and η, and N _S represents the number of sample points.

Since the depth of the center of the sound source underwater is Z _S , the depth Z _S1 of each element of the vertically arranged sound source is Z _Sl =Z _S +ξ _Sl =Z _S ξ _SO +lΔξ _S (17).

Regarding the equivalent linear array receiver for the directional receiver, the position and weight (receiving sensitivity) of each receiving element can be determined in a similar manner. Z _Rp for each,
W _RP , p=0, 1, ..., N _R -1,
Sample data of receiving directivity corresponding to E _Sn and D _Sn in equation (15) are expressed as E _Rq , D _Rq , q=0, 1,...N _R
-1.

Next, the sound field and receiver output when the equivalent linear array sound source and the equivalent linear array receiver are placed in the actual ocean will be explained. The sound pressure at the position (r, 0, Z _Rp ) caused by a non-directional point sound source with a weight of 1 placed at the position (0, 0, Z _Sl ) is expressed by the following equation.

Ψ(r, Z _Rp Z _Sl ;f)=iπ _∞ 〓 ⁿ⁼⁰ U _o (Z _Sl ) U _o (Z _Rp ) H _O ⁽¹⁾ π(k _o r) ……(18) Therefore, the above By equivalent array sound source (r, 0,
The sound pressure at the position Z _Rp ) can be obtained by multiplying equation (18) by the weight W _Sl of the sound source element and superimposing the result with respect to l, as follows.

P(r, Z _Rp , Z _S :f)=iπ _∞ 〓 ⁿ⁼⁰ U _o (Z _Rp ) H _O ⁽¹⁾ (k _o r) _NS-1 〓 ^l=0 W _Sl U _o (Z _Sl ) ...(19) Using vector representation, W _S = (W _SO , W _S1 , ...W _S , N _S −1) ^T ...(20) U _oS = [U _o (Z _SO ), U _o (Z _S1 ), ..., U _o (Z _S , N _S −1)] ^T ... (21), then equation (17) can be rewritten as follows.

P(r, Z _Rp , Z _S :f)=iπ _∞ 〓 ⁿ⁼⁰ (U _oS , W ^* _S ) U _o (Z _Rp ) H _p ⁽¹⁾ (k _o r) ……(22) (* is a sign that means a conjugate complex number, the same applies hereinafter) The received wave output of the equivalent linear array receiver is (22)
The equation is multiplied by the weight W _Rp of the receiving element, and the result is superimposed for all receiving elements, resulting in the following equation.

V(r)=iπ _∞ 〓 ⁿ⁼⁰ (U _oS , W ^* _s ) H _p ⁽¹⁾ (k _o r) _NR-1 〓 ^p=0 W _Rp U _o (Z _rp ) ……(23) Vector Using the expression W _R = (W _RO , W _R1 ,...W _R , N _R −1) ^T
…(24) U _oR = [U _o (Z _RO ), U _o (Z _R1 ), …, U _o (Z _R , N _R −1)] ^T … (25) If expressed as (23) The formula can be rewritten as follows.

V(r)=iπ _∞ 〓 ⁿ⁼⁰ (U _oS W ^* _S ) (U _oR W ^* _R ) H _O ⁽¹⁾ (k _o r)=iπ _∞ 〓 ⁿ⁼⁰ G _o H _O ⁽¹⁾ ( k _o r ) ... (26) Here, G _o = (U _oS W ^* _S ) (U _{o R} W ^* _R ) ... (27). Preset distance r _j , j=0, 1,...
...Determining to find the value of V(r) for K-1, V=[V(r _O ), V(r ₁ ), ..., V(r _K-1 )] ^T
...(28) H _o = iπ[H _O ⁽¹⁾ (k _o r _O ), H _O ⁽¹⁾ (k _o r ₁ ),...
, H _O ⁽¹⁾ (k _o r _K-1 )] ^T ... (29), then equation (26) can be rewritten as follows.

V = _∞ 〓 ⁿ⁼⁰ G _o H _o = _∞ 〓 ⁿ⁼¹ V _o ... (30) For simplicity, we will also use vector representation for the following quantities.

D _S = (D _SO , D _S1 , ... D _S , N _S-1 ) ^T ... (31) D _R = (D _RO , D _R1 , ..., D _R , N _R-1 ) ^T ... (32) Z _SA = (Z _SO , ..., Z _Sl Z _S , N _S-1 ) ^T ... (33) Z _RA = (Z _RO , Z _R1 , ..., Z _R , N _R-1 ) ^T ...(34) Next, examples of the present invention will be described. 1st
The figure is an explanatory diagram of an embodiment of the present invention. In the figure, 1 is an input terminal for sound source information, 2 is an input terminal for receiver information, 3 and 4 are input terminals for the sound source/receiver information, and output sample values of directivity and equivalent linear array element positions. interpolator, 5, 6, 12-15,
19 is a register, 7 and 8 input directional sample values, and output the weight of each element of the equivalent linear array.
DFT circuit, 9 is an input terminal for environmental information, 10 is sound source depth, 10 is sound source depth, receiving depth, environmental information,
A normal mode analysis device inputs a normal mode number and calculates the values at the source, reception depth, and eigenvalue of the eigenfunction of the corresponding normal mode wave; 11 inputs the normal mode eigenvalue and calculates the vector of equation (29); A distance function calculation device for calculating H, 16 and 17 are vector multipliers, 18 is a scalar multiplier, 20 is a scalar vector multiplier, 21 is a counter, 22 is an adder, and 23 is an output terminal.

The operating principle of this device will be explained. First, the counter 21 is set to an initial value n=1, and the accumulator 22 is set to an initial value V=O. Sound source information input from input terminal 1 is sent to interpolator 3. The interpolator 3 creates sound source directivity sample data D _S with a spatial frequency interval Δη _S from the sound source information,
The element spacing Δξ and element depth vector Z _SA of the equivalent linear array sound source are calculated using equations (16) and (17).
The sound source directivity sample data D _S is transferred to the DFT circuit 7 along with Δη _S and Δξ _S , and the element depth vector is
Z _SA is transferred to register 5.

The DFT circuit 5 calculates the weight vector W S of the equivalent linear array sound source from the data D _S , Δη _S , Δξ _S transferred from the interpolator 3 according to equation (14), and calculates the weight vector W _S of the equivalent linear array sound source according to equation (14 ₎ .
is transferred to register 12.

On the other hand, receiver information input from the input terminal 2 is sent to the interpolator 4. The interpolator 4 generates receiver directivity sample data D _R with a spatial frequency interval Δη _R from the receiver information, and further calculates the element spacing Δξ _R and element depth Z _RA of the equivalent linear array receiver. . D _R is transferred to DFT rotation 8 along with Δη _R and Δξ _R , and Z _RA is sent to register 6.

Environmental information such as the sound speed profile, wind speed, seabed characteristics, etc. is inputted from the input terminal 9 and sent to the normal mode analysis device 10.

The normal mode analysis device 10 first reads the normal mode number n from the counter 21. The counter 21 has its initial value set to 1, and the accumulator 2
The value is increased by 1 by the control signal from 2. Next, a normal mode eigenvalue k _o is determined from the normal mode number n and the environment information, and k _o is transferred to the distance function calculation device 11 . Also, read the depth Z _SA of the arrayed sound source element from the register 5,
The value of the eigenfunction at each element depth is calculated to form a sound source depth eigenfunction vector U _oS [Equation (21)] and sent to the register 14. Furthermore, the depth Z _RA of the array receiver element is read from register 6, the value of the eigenfunction at each element depth is calculated to form the receiving depth eigenfunction vector U _oR [formula (25)], and
Send to 5.

The distance function calculation device 11 calculates a preset distance r _j , J=0, 1, from the normal mode eigenvalue k _o sent from the normal mode analysis device 10.
..., distance function value iπH _O ⁽¹⁾ (k _o r _j ) at K−1
is calculated, and a distance function vector H _o [formula (29)] is formed and stored in the register 19.

The vector multiplier 16 reads the weight vector W _S of the arrayed sound source from the register 12 and the sound source depth eigenfunction vector U _oS from the register 14, and calculates the sound source function value Sn=(U _oS , W ^* _S ). and transfers it to the scalar multiplier 18.

A vector multiplier 17 receives the weight vector W _R of the array receiver from the register 13 and receives the reception depth eigenfunction vector U _oR from the register 15.
is read, the receiver function value R _o = (U _oR , W ^* _R ) is calculated and transferred to the scalar multiplier 18 .

In the scalar multiplier 18, the vector multiplier 16
The product G _o [formula (27)] of the sound source function value S _o sent from the source function value S o and the receiver function value R _o is determined and transferred to the scalar vector multiplier 20 .

The scalar vector multiplier 20 uses the G _o sent from the scalar multiplier 18 and the register 19
The product v _o with the distance function vector H _o stored in is calculated and sent to the accumulator 22 . The accumulator 22 accumulates the v _o sent from the scalar vector multiplier 20 to V and sends a control signal to the counter 21 .

Upon receiving the control signal from the accumulator 22, the counter 21 increments the mode number n by 1 and sends it to the normal mode analysis device 10.

In this way, the output V for the normal mode number n is accumulated one after another. The normal mode analysis device 10 stops execution when the normal mode number n or the imaginary part I _n (k _o ) of the normal mode eigenvalue exceeds each preset reference value.

(Effects of the Invention) As explained above, in the embodiment of the present invention, an arrayed sound source/arrayed receiver having vertical directivity equivalent to a directional sound source/receiver is provided under the defining conditions of the directivity ( Since it is produced under a homogeneous infinite medium,
Disturbances in directivity due to inhomogeneity of the medium or discrepancy in sound speed with the defining conditions can be correctly evaluated, and the concept of an equivalent sound line for normal mode waves is not required, so it can be used in cases where the wave is asymmetric with respect to the horizontal axis. It has the advantage that the effect can be accurately evaluated for any directivity including .

In addition, since the received wave output is decomposed into normal mode waves and a form of variable separation is used for each mode as shown in equation (26), it is similar to the conventional sound wave propagation simulation method using a normal mode model. The device can be realized with a simple configuration.

The present invention has the advantage of being able to correctly evaluate the directivity effect of a directional sound source/receiver even under environments different from the defining conditions of the directivity, so it is suitable for use with audio equipment whose usage environment conditions vary widely. It is useful for on-site performance evaluation of sonar equipment used in the ocean, etc.

[Brief explanation of drawings]

FIG. 1 is an explanatory diagram of an embodiment of the present invention, FIG. 2 is an explanatory diagram of the present invention using cylindrical coordinates, and FIG. 3 is a diagram showing the relationship between directivity and a linear sound source. 1...Sound source information input terminal, 2...Receiver information input terminal, 3, 4...Interpolator, 5, 6, 12-1
5, 19...Register, 7, 8...DFT circuit,
9...Environmental information input terminal, 10...Normal mode analysis device, 11...Distance function calculation device, 16,
17... Vector multiplier, 18... Scalar multiplier, 20... Scalar vector multiplier, 21...
... Counter, 22 ... Accumulator, 23 ... Output terminal.

Claims (1)

[Claims] 1 In a sound wave propagation simulation method that simulates the sound wave propagation phenomenon in a heterogeneous medium,
Equivalent linear array calculation means for calculating the position and weight of each element of an equivalent vertical linear array sound source/receiver having the vertical directivity pattern for a given frequency and vertical directivity pattern of the directional sound source/receiver. Calculated by the frequency and the normal mode eigenvalue of the specified mode number under the given environmental conditions,
The eigenfunction value at the depth of each element of the array sound source/receiver is calculated by normal mode analysis means, and the normal mode distance function at a plurality of set distances is calculated as a distance function for the normal mode eigenvalue. The depth function calculation means calculates the cumulative sum of the product of the weight of each element of the array sound source and the eigenfunction value, and the cumulative sum of the product of the weight of each element of the array receiver and the eigenfunction value. and the product of the two cumulative sums (referred to as the depth function)
is calculated, the product of the depth function and the distance function is calculated using the output calculation means, the results are accumulated for all normal modes, and the output of the directional receiver placed in the sound field due to the directional sound source is calculated. A sound wave propagation simulation method characterized by calculating.