JPS63282679A - Sonic wave propagation simulation system - Google Patents
Sonic wave propagation simulation systemInfo
- Publication number
- JPS63282679A JPS63282679A JP11772787A JP11772787A JPS63282679A JP S63282679 A JPS63282679 A JP S63282679A JP 11772787 A JP11772787 A JP 11772787A JP 11772787 A JP11772787 A JP 11772787A JP S63282679 A JPS63282679 A JP S63282679A
- Authority
- JP
- Japan
- Prior art keywords
- depth
- sound
- function
- calculation means
- normal mode
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000004088 simulation Methods 0.000 title claims abstract description 15
- 238000004364 calculation method Methods 0.000 claims abstract description 74
- 238000000034 method Methods 0.000 claims abstract description 20
- 230000006866 deterioration Effects 0.000 abstract description 2
- 238000009825 accumulation Methods 0.000 description 9
- 230000008859 change Effects 0.000 description 5
- 238000007796 conventional method Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 5
- 238000012545 processing Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 238000013459 approach Methods 0.000 description 3
- 238000010606 normalization Methods 0.000 description 3
- 230000005856 abnormality Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000004913 activation Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 210000003127 knee Anatomy 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000009304 pastoral farming Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000005855 radiation Effects 0.000 description 1
- 239000013535 sea water Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Landscapes
- Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
Abstract
Description
【発明の詳細な説明】
〔産業上の利用分野〕
本発明は、水中の音波伝搬特性を解析する音波伝搬シミ
ュレーション方式に関し、特にノーマルモード論理に基
づいた音波伝搬シミュレーション方式に関するものであ
る。DETAILED DESCRIPTION OF THE INVENTION [Field of Industrial Application] The present invention relates to a sound wave propagation simulation method for analyzing underwater sound wave propagation characteristics, and particularly to a sound wave propagation simulation method based on normal mode logic.
ソーナー等の海洋音響機器の性能は、その運用現場の海
洋環境により大きく変化する。従って、機器の効果的運
用のために、海洋中の音波伝搬特性を把握することが重
要であり、種々の音波伝搬モデルが開発されている。The performance of marine acoustic equipment such as sonar varies greatly depending on the marine environment in which it is operated. 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.
ノーマルモード・モデルもその−っでアル。角周波数ω
の単位強度の無指向性点音源による音場は、媒質特性が
水平方向に変化しない場合、ノーマルモード波の累積和
として次式のように表わされる。The normal mode model is also the same. angular frequency ω
When the medium characteristics do not change in the horizontal direction, the sound field caused by a non-directional point sound source with unit intensity is expressed as the cumulative sum of normal mode waves as shown in the following equation.
ψ(r、z 職; z s・ ω )=πiΣA、、
U、(z 5)Ua(z m)Ha(1)n−/
(k、r) ・・・・(1)ここ
で、座標系として円筒座標系<r、φ、z)を用い、音
源直上の海面を原点として、Z軸を深さ方向にとってい
る。音源位置は(0,O,zs)、受波器位置は(r、
0.2++)とした。nはノーマルモード波のモード番
号、k、は第nモードの固有値、Ufi(z)は第nモ
ードの固有関数、A、は正規化定数であり、H,(1)
は0次の第1種ハンゲル関数である。ψ(r, z position; z s・ω)=πiΣA,,
U, (z 5) Ua (z m) Ha (1) n-/ (k, r) ... (1) Here, using the cylindrical coordinate system < r, φ, z) as the coordinate system, the sound source The origin is the sea surface directly above, and the Z axis is taken in the depth direction. The sound source position is (0, O, zs), and the receiver position is (r,
0.2++). n is the mode number of the normal mode wave, k is the eigenvalue of the nth mode, Ufi(z) is the eigenfunction of the nth mode, A is a normalization constant, and H, (1)
is a Hungel function of the first kind of order 0.
U。(z)は次式
%式%
及び以下に述べる境界条件を満足する。k(z)は角周
波数ωの音波に対する媒質の波数である。境界条件は次
のように設定きれる。U. (z) satisfies the following formula % formula % and the boundary conditions described below. k(z) is the wave number of the medium for a sound wave of angular frequency ω. The boundary conditions can be set as follows.
(1)海面で音圧Oとなる。(1) The sound pressure becomes O at sea level.
(2)0<z<ooの全領域でU、(z)及びdU、/
dzはそれぞれ連続である。(2) U, (z) and dU, / in the entire range 0<z<oo
Each of dz is continuous.
(3ン2→■でSommerfeldの放射条件を満足
する。但し、海底特性が海底表面z=zllにおける反
射係数の形で与えられている場合は、U。/(dU、、
/dz)z=zllが設定された値をとる。(3n2→■ satisfies Sommerfeld's radiation condition. However, if the seafloor characteristics are given in the form of the reflection coefficient at the seafloor surface z=zll, then U./(dU, ,
/dz) z=zll takes the set value.
波数の鉛直分布k(Z)C以下、これを波数プロファイ
ルと呼ぶ)は簡単な関数で近似されるが、通常はO<z
<ωの全範囲を1つの関数で近似することが困難であり
、媒質をいくつかの層に分割して、各層内をそれぞれ別
の関数で表わすという方法がとられている。これに伴い
、前記固有関数U。(Z)も各層でそれぞれ関数を用い
て表わされる。The vertical distribution of wavenumbers k(Z)C (hereinafter referred to as wavenumber profile) is approximated by a simple function, but usually O<z
It is difficult to approximate the entire range of <ω with one function, so a method is used in which the medium is divided into several layers and each layer is represented by a different function. Accordingly, the eigenfunction U. (Z) is also expressed using a function in each layer.
しかし、層数が多い場合は、ノーマルモード固有(a
k 、を厳密に求めるには多大な計算時間とダイナミッ
ク・レンジが必要であり、そのためしばしばWKB近似
による固有値及び固有関数の近似計算が行なわれる。However, when the number of layers is large, normal mode-specific (a
Exactly determining k requires a large amount of calculation time and dynamic range, and therefore approximate calculations of eigenvalues and eigenfunctions are often performed using WKB approximation.
固有値k。を実部と虚部とに分け
に、、=ξo+iα1 ・・・・(
3)とおくと、実部ξ。は次の方程式を解くことに近似
的に求められる。Eigenvalue k. Divide into real part and imaginary part, =ξo+iα1...(
3), the real part ξ. can be approximately found by solving the following equation.
εS(ξ。)+ε、(ξ。)−2πn ・・・べ
4)ここで、z4(ξ、)はξ、に対応した等価音線の
最上端(海面又は回転点)の深度であり、za(ξ。)
は上記等価音線の最下端(海底又は回転点)の深度であ
る。また、ε、(ξハ、ε、(ξ。)はそれぞれ、上記
等価音線の最上端及び最下端における位相変化量を表わ
す。ε、(ξ。)は通常、等価音線が海面反射する場合
は一π、海面反射せず回転点(apex)を形成する場
合は−π/2と近似される。また、ε8(ξ。)は等価
音線が海底反射する転点(nadir)を形成する場合
は、−π/2とする。εS (ξ.) + ε, (ξ.) - 2πn ... Be4) Here, z4 (ξ,) is the depth of the top end (sea surface or rotation point) of the equivalent sound ray corresponding to ξ, za (ξ.)
is the depth of the lowest end (seafloor or rotation point) of the above equivalent sound ray. In addition, ε, (ξc, ε, (ξ.) represent the amount of phase change at the top and bottom ends of the above equivalent sound ray, respectively. ε, (ξ.) usually represents the amount of phase change when the equivalent sound ray is reflected from the sea surface. It is approximated as -π/2 when the sound ray does not reflect on the sea surface and forms a rotation point (apex). Also, ε8 (ξ.) forms a turning point (nadir) where the equivalent sound ray reflects on the sea floor. If so, set it to -π/2.
固有値の虚部α。は次式で得られる等価音線のサイクル
距離り。Imaginary part α of the eigenvalue. is the cycle distance of the equivalent sound ray obtained by the following formula.
・・・・く5)
を用いて、
α。=−P。IR8(ξ。〉・R3(ξ、)l/D、・
・・(6)
で近似できる。ここでR5(ξ。〉は海面反射係数であ
る(高周波で水中における吸収減衰量が大きい場合には
さらにその量を付加しなければならない)。...5) Using α. =-P. IR8(ξ.〉・R3(ξ,)l/D,・
...(6) can be approximated. Here, R5(ξ.) is the sea surface reflection coefficient (if the amount of absorption and attenuation in water is large at high frequencies, that amount must be added).
固有関数U。(z)はWKB近似は次式の形で与えられ
る。Eigenfunction U. The WKB approximation for (z) is given in the following form.
tJ、(z)−ω””Ck”(z )−ξ% :)−
1/4(B+exp[j 1)(Z、、Z、ξ。〉〕+
B*exp[ニーjp(Z、、z、ξa))・・・・(
7)
ここでp(z、、z、ξ、)は次式で定義される。tJ, (z)-ω""Ck"(z)-ξ%:)-
1/4(B+exp[j 1)(Z,,Z,ξ.〉]+
B*exp[knee jp(Z,, z, ξa))...(
7) Here, p(z,,z,ξ,) is defined by the following equation.
・・・・く8) Bt、Btは境界条件により定まる定数である。...ku8) Bt and Bt are constants determined by boundary conditions.
等価音線が海面反射するようなモード波の場合、2.=
0となり、Bl= 1/2 j 、B!−−12jとお
くことができて、U、(z )はU、(z)= ω
l/ 2 (k 1 (z )−ξ 、z 〕 −
”’sin[1)(0,z、ξ、)〕
・・・・(9)と簡単になる。また、等側音線が海面
に達せずaexpを形成するようなモードの場合には、
B1−exp(−j π/ 4 )/ 2、B 、=
−exp(jπ/4)/2とおくことができU、(z)
は、
Un(Z )= ω”2[k ”(Z )−ξa”)−
”’!B1n[p(z、、z、ξ。)+π/4〕・・・
・(10)となる。In the case of a mode wave in which the equivalent sound ray is reflected from the sea surface, 2. =
0, Bl= 1/2 j, B! −−12j, and U, (z) is U, (z) = ω
l/ 2 (k 1 (z)−ξ, z] −
”'sin [1) (0, z, ξ,)]
...(9) becomes simple. In addition, in the case of a mode in which the isolateral sound ray does not reach the sea surface and forms an aexp,
B1-exp(-jπ/4)/2,B,=
−exp(jπ/4)/2 can be set as U, (z)
is Un(Z)=ω”2[k”(Z)−ξa”)−
”'!B1n[p(z,,z,ξ.)+π/4]...
・It becomes (10).
上述のWKB近似法はノーマルモード解を厳密に求める
場合に比べて計算が非常に容易であり、また送受深度が
等側音線の上下端から十分離れた場合には良好な結果が
得られる。また、等側音線が海面や海底で反射している
場合には、送受波深度がその反射点に近づいても精度は
保たれるという特徴を有している。しかしながら等側音
線の上下端のいずれかが回転点となり、しかも送波深度
又は受波深度の少なくとも一方がその回転点付近にある
場合には、上記(9)、(10)式の〔k2(z)−ξ
g )−1/4の項が非常に大きな値になり、回転点
に近付くとともに無限大に発散してしまいノーマルモー
ド波の計算ができなくなってしまう。特に低周波では、
上述の振幅の異常が生じる深度範囲は広くなる。従って
、広い領域にわたって音圧分布に誤差が生じてしまうと
いう問題があった。The above-mentioned WKB approximation method is much easier to calculate than strictly determining the normal mode solution, and good results can be obtained when the transmitting and receiving depths are sufficiently far from the upper and lower ends of the isolateral sound ray. Furthermore, when the isolateral sound ray is reflected from the sea surface or the seabed, accuracy is maintained even if the depth of transmission and reception approaches the reflection point. However, if either the upper or lower end of the isolateral sound ray becomes a rotation point, and at least one of the transmission depth or reception depth is near the rotation point, then in equations (9) and (10) above, [k2 (z)−ξ
g) -1/4 becomes a very large value and diverges to infinity as it approaches the rotation point, making it impossible to calculate the normal mode wave. Especially at low frequencies,
The depth range in which the above-mentioned amplitude abnormality occurs becomes wider. Therefore, there is a problem in that an error occurs in the sound pressure distribution over a wide area.
本発明は上述の点に鑑みてなされたもので、上記回転点
近傍で計算不能が生じたり、近似誤差が拡大したりする
という問題点を除去し、簡単な処理によって高い精度が
実現できる実用的な音波伝搬シミュレーション方式を提
供することにある。The present invention has been made in view of the above-mentioned points, and is a practical method that eliminates problems such as inability to calculate near the rotation point and expansion of approximation errors, and achieves high accuracy through simple processing. The purpose of this invention is to provide a sound wave propagation simulation method.
上記問題点を解決するため本発明は、音波伝搬シミュレ
ーション方式を、音源深度及び受波点深度におけるノー
マルモード波の固有関数の振幅値をノーマルモード波の
回転点深度から音源深度及び受波点深度までの位相積分
値の低次関数を用いて表わす固有関数計算手段を具備し
、該固有関数計算手段によって得られる音源深度及び受
波点深度におけるノーマルモード関数値を用いて音場計
算を行なうようにした。In order to solve the above problems, the present invention uses a sound wave propagation simulation method to calculate the amplitude value of the eigenfunction of the normal mode wave at the sound source depth and the receiving point depth from the rotation point depth of the normal mode wave to the sound source depth and the receiving point depth. The sound field is calculated using the normal mode function values at the sound source depth and the receiving point depth obtained by the eigenfunction calculation means. I made it.
上記の如く音波伝搬シミュレーション方式を構成するこ
とにより、固有関数を回転点付近でも発散したり不安定
になったりすることがない安定した精度の良い関数で近
似しているので、従来の方式のように回転点近傍で精度
が劣化し、値が発散してしまい計算不能になるという問
題点を除(ことができる。By configuring the sound wave propagation simulation method as described above, the eigenfunction is approximated by a stable and accurate function that does not diverge or become unstable even near the rotation point, so it is different from the conventional method. This eliminates the problem that accuracy deteriorates near the rotation point, causing values to diverge and making calculations impossible.
また、固有関数の近似を非常に簡単な処理により実現し
ているので、従来の方式に比較して処理時間の増加は殆
どなくソーナーの運用現場での性能評価等、実時間性を
要求されるシステムに利用することが可能となる。In addition, since the approximation of the eigenfunction is achieved through very simple processing, there is almost no increase in processing time compared to conventional methods, and real-time performance is required, such as performance evaluation in the field of sonar operation. It becomes possible to use it in the system.
以下、本発明の詳細な説明する。 The present invention will be explained in detail below.
まず、始めに本発明の詳細な説明する。上記(9)、(
10)式による近似は、位相については回転点の近傍で
も十分な精度が保たれるで、振幅項についてのみ精度を
改善することを考える。First, the present invention will be explained in detail. (9) above, (
The approximation using equation 10) maintains sufficient accuracy for the phase even in the vicinity of the rotation point, so consider improving the accuracy only for the amplitude term.
音源又は受波点が上側の回転点の近傍にある場合を例に
とり説明する。An example will be explained in which the sound source or wave receiving point is near the upper rotation point.
回転点深度z=z、の近傍で、音速C(z)=ω/k(
k)が
C−”(z)=C,−”〔1−2g、(z−z、、)/
C,)・・・(11)
の形で近似できるものとする。ここで、C0=ω/ξ。Near the rotation point depth z=z, the sound speed C(z)=ω/k(
k) is C-”(z)=C,-”[1-2g, (zz-z,,)/
C, )...(11) It is assumed that it can be approximated in the following form. Here, C0=ω/ξ.
は基準音速である。is the reference speed of sound.
次の置きかえ
a3=−2g、ξn”/ Ca、77 = a(z
zu)・・・・(12)
を行なうと、(2)式は次のAi ryの方程式に簡略
化される。Next replacement a3 = -2g, ξn”/Ca, 77 = a(z
zu)...(12) When the equation (2) is performed, the equation (2) is simplified to the following Air ry equation.
21”Ulつη2−ηU=0 ・・・・(
13)このとき、(7)式に対応する(12)式の解は
次式で与えられる。21"Ul η2-ηU=0...(
13) At this time, the solution of equation (12) corresponding to equation (7) is given by the following equation.
U、(z)=(πω/a)””(βIcBI(η)〕+
βtcBI(7)−jAI(7)〕)
・・・・(14)
但し、AI(77)、 Bl(77)はAi ry関数
である。U, (z)=(πω/a)””(βIcBI(η))]+
βtcBI(7)-jAI(7)]) (14) However, AI(77) and Bl(77) are Airy functions.
上記(14)式は海面における境界条件を満たす必要が
あるが、海面は回転点から十分離れているとすると、β
1=−βz= 1 / 2 Jとおくことができ、U、
、(z )= (π(L) / a )””A +(7
7)−(15)となる。Equation (14) above needs to satisfy the boundary condition at the sea level, but assuming that the sea level is sufficiently far from the rotation point, β
1=-βz=1/2 J, and U,
, (z)=(π(L)/a)””A+(7
7)-(15).
ここでa、ηを別の形で表わすことを考える。Let us now consider expressing a and η in another form.
(11)式を(8)式に代入すると、
C,) ””d z = 2/ 3 (−(2ξ、”g
、)/C、) ””< z −z 、)””
・=(16)従って、(12)、(16)式により、次
式が得られる。Substituting equation (11) into equation (8), C,) ""d z = 2/3 (-(2ξ, "g
,)/C,) ””< z −z ,)””
.=(16) Therefore, from equations (12) and (16), the following equation is obtained.
η=−(3/2p(z。、z、ξ。))””・・・・(
17)また、(11)、(12)式から、次式を得る。η=-(3/2p(z., z, ξ.))""...(
17) Also, from equations (11) and (12), the following equation is obtained.
a−(k”(z)−ξ * ) l / ! / (n
) l / t・・・・(18)
(18)式を用いると(15)式は次のように書き改め
られる。a-(k”(z)-ξ*) l/!/(n
) l/t (18) Using equation (18), equation (15) can be rewritten as follows.
U、(z)−(yrω)””(−77/(k”(z)−
ξ、′)) ”’A、(η) ・・・・(
19)(19)式には(9)、(10)式と同様(k2
(2)−ξ。り−1/4の項が含まれているが、この項
の回転点付近での増大は(−η)1′4の項により相殺
され、従来の方式の問題点は基本的に解消される。また
、ηを(17)式の形で与えたことにより、波数分布が
(11)式からはずれた部分でも用いることができる。U, (z)-(yrω)””(-77/(k”(z)-
ξ, ′)) ”'A, (η) ・・・(
19) In equation (19), as in equations (9) and (10), (k2
(2) −ξ. Although a term of -1/4 is included, the increase of this term near the rotation point is canceled out by the term of (-η)1'4, and the problems of the conventional method are basically solved. . Further, by giving η in the form of equation (17), it can be used even in a portion where the wave number distribution deviates from equation (11).
しかしく19)式中の(−η/(k”(z)−ξ 1〕
)l/4(7)部分はξ。=k(z)とすると不定にな
ってしまうため、別の形で表現することが必要である。However, (−η/(k”(z)−ξ 1) in formula 19)
) l/4 (7) The part is ξ. =k(z), it becomes undefined, so it is necessary to express it in another form.
これは(12)式を(15)に代入すると容易に得られ
る。This can be easily obtained by substituting equation (12) into (15).
U、、(Z)=”(WC,)””Iω/ 2 g −1
”’A +(η)・・・・(20)
従って、ηの値によって(19)、(20)式を使い分
けるようにすれば、固有関数U。(z)を精度よく求め
ることができ、それを(1)式に代入れすることによっ
て、音場を高精度でシミュレートできる。U,, (Z)=”(WC,)””Iω/ 2 g −1
``'A + (η)...(20) Therefore, by using equations (19) and (20) depending on the value of η, the eigenfunction U.(z) can be calculated with high accuracy, By substituting it into equation (1), the sound field can be simulated with high accuracy.
ここで、(19)式の計算を更に簡略化することを考え
る。(19)式と(10)式との比をとり、次の量を定
義する。Here, let us consider further simplifying the calculation of equation (19). The following quantity is defined by taking the ratio between equation (19) and equation (10).
q(η)−π1″(−η)’ ′4A + (77)/
5in(φ十π/4〉 ・・・・
(21)ここでφはz6から2までの位相積分値φ−p
(z、、z、ξ。)=2/3(−η)3″・・・・(2
2)
である、qをηの関数としてプロットすると第2図のよ
うになる。q(η)−π1″(−η)′′4A + (77)/
5in (φ1π/4〉 ・・・・
(21) Here φ is the phase integral value φ−p from z6 to 2
(z,,z,ξ.)=2/3(-η)3″...(2
2) When q is plotted as a function of η, the result is shown in Figure 2.
第2図かられかるように、q(η)はηの非常に滑らか
な関数になっており、ηが負の大きな値になると1に漸
近していく。従っ工、q(η〉を簡単な関数形で近似し
ておけば、(10)式のWKB近似を次のように補正す
ることにより、簡単な計算で高精度の固有関数を得るこ
とができる。As can be seen from Figure 2, q(η) is a very smooth function of η, and approaches 1 as η becomes a large negative value. Therefore, if q(η〉) is approximated by a simple functional form, a highly accurate eigenfunction can be obtained by simple calculation by correcting the WKB approximation of equation (10) as follows. .
U、(Z)=(L)””(k”(z)−ξm” 〕−”
’5IJI (φ(η)+π/4)q(η) ・・
・・(23)一方、(20)式は音速プロファイルが次
式C”’1(z)−C+−”(1−2g +(z −z
H>/Cr )l z、≦2≦z l+t * j
w 1 +−−−−r N・・・・(24)
の形で表わされている場合にはさらに簡単になる。回転
点がj層内(z+≦2≦21*+)にある場合を考える
、(24)式を2で微分し整理すると、g(z)−〇C
(z)/’2)Z”g+c3<z)/C+”となり、z
”z、を代入すると、次式を得る。U, (Z)=(L)""(k"(z)-ξm"]-"
'5IJI (φ(η)+π/4)q(η)...
...(23) On the other hand, in equation (20), the sound velocity profile is as follows
H>/Cr)l z, ≦2≦z l+t * j
It becomes even simpler if it is expressed in the form w 1 +−−−−r N (24). Considering the case where the rotation point is in the j layer (z+≦2≦21*+), differentiating equation (24) by 2 and rearranging it, we get g(z)−〇C
(z)/'2)Z"g+c3<z)/C+", and z
By substituting “z,” we obtain the following equation.
g、=g(z、)−gバc、/c、)” ・−−−
(25)(24)式を(20)式に代入すると2が2.
≦ZSZ、やtの範囲内にある場合には、U、(z)−
(πG、)”” I ω/ 2 g + 1 ”’AI
(77)・・・・(26)
となることがわかる。g,=g(z,)-gbac,/c,)” ・---
(25) Substituting equation (24) into equation (20), 2 becomes 2.
If it is within the range of ≦ZSZ, or t, then U, (z)−
(πG,)"" I ω/ 2 g + 1 "'AI
It can be seen that (77)...(26).
固有関数U。(z)を(26)式の形で扱わねばならな
いのはη−0の近傍の非常に狭い範囲である。この範囲
ではA、(η)は極めて滑らかな関数となり、1次式等
の低次式で十分近似することができる。Eigenfunction U. It is in a very narrow range around η-0 that (z) must be treated in the form of equation (26). In this range, A and (η) become extremely smooth functions, and can be sufficiently approximated by a low-order equation such as a linear equation.
以上から、深度2が回転点深度zuに十分近い場合には
(26)式を用い、ある程度以上離れている場合には(
23)式の補正式を用いることにすれば、従来の方法で
問題となった回転点近傍での精度劣化や計算異常(無限
大発散)の現象が解消され、簡単な計算で精度のよい固
有関数値を得ることができる。また(26)式はη〉0
の領域にも適用できる。From the above, if depth 2 is sufficiently close to the rotation point depth zu, use equation (26), and if it is a certain distance or more away, use (
23) By using the correction formula shown in Equation 23), the problem of accuracy deterioration near the rotation point and calculation abnormalities (infinite divergence), which were problems with the conventional method, can be solved, and accurate eigenvalues can be obtained with simple calculations. You can get the function value. Also, equation (26) is η〉0
It can also be applied to the following areas.
音源又は受波器深度が下側の回転点(nadirと呼ぶ
)深度z4の近くにある場合にも、同様の方法で精度の
よい固有関数が得られる。この場合、φ、ηは次式で与
えられる。If the sound source or receiver depth is near the lower rotation point (referred to as nadir) depth z4, a highly accurate eigenfunction can be obtained in a similar manner. In this case, φ and η are given by the following equations.
φ−p(Z、Za、ξ、) −p (z 、、 z a
、ξ、)−p(z、、z、ξ。) ・・
・・(27)η−−(3/2φ)8′3
・・・・(28)本発明の音波伝搬シミュレーション方
式は上記原理に基づいてなされたもので、次に音波伝搬
シミュレーション方式を実施する例を図面に基づいて説
明する。φ−p(Z, Za, ξ,) −p(z ,, za
,ξ,)−p(z,,z,ξ.)...
...(27)η--(3/2φ)8'3
(28) The sound wave propagation simulation method of the present invention is based on the above principle. Next, an example of implementing the sound wave propagation simulation method will be described with reference to the drawings.
第1図は、本発明に係る音波伝搬シミュレーション方式
を実施する装置の構成を示すブロック図である。図示す
るように該装置は、入力端子1、初期設定手段2、固有
値実部計算手段3、位相積分計算手段4、深度関数生成
手段5、固有関数計算手段6、固有値虚部計算手段7、
距離関数生成手段8、音場累積手段9及び出力端子10
で構成されている。FIG. 1 is a block diagram showing the configuration of an apparatus that implements a sound wave propagation simulation method according to the present invention. As shown in the figure, the device includes an input terminal 1, initial setting means 2, eigenvalue real part calculation means 3, phase integral calculation means 4, depth function generation means 5, eigenfunction calculation means 6, eigenvalue imaginary part calculation means 7,
Distance function generation means 8, sound field accumulation means 9 and output terminal 10
It is made up of.
初期設定手段2は、前記入力端子1に入力される環境条
件及び周波数等の入力データに基づき、音速プロファイ
ル、音波の減衰係数、海面・海底反射係数等を設定する
ものである。The initial setting means 2 sets the sound speed profile, sound wave attenuation coefficient, sea surface/bottom reflection coefficient, etc. based on input data such as environmental conditions and frequencies inputted to the input terminal 1.
固有値実部計算手段3は、前記初期設定手段2から送ら
れる音速プロファイル、海底反射係数、周波数等に基づ
きノーマルモード固有値の実部を計算するものである。The eigenvalue real part calculating means 3 calculates the real part of the normal mode eigenvalue based on the sound speed profile, seabed reflection coefficient, frequency, etc. sent from the initial setting means 2.
位相積分計算手段4は、前記初期設定手段2から送られ
てくる音速プロファイル、周波数及び前記固有値実部計
算手段3又は固有関数計算手段6から送られてくる固有
値実部及び積分の上下限深度の値から、位相積分値を算
出するものである。The phase integral calculation means 4 calculates the sound velocity profile and frequency sent from the initial setting means 2, the eigenvalue real part sent from the eigenvalue real part calculation means 3 or the eigenfunction calculation means 6, and the upper and lower limit depths of integration. The phase integral value is calculated from the value.
深度関数生成手段5は、ノーマルモード波の正規化定数
A。を計算し、更に固有関数計算手段6を用いて音源及
び受波器深度における固有関数値を計算して深度関数2
゜(Z S r Z * )を生成するものである。The depth function generating means 5 generates a normalization constant A for normal mode waves. Then, the eigenfunction calculation means 6 is used to calculate the eigenfunction values at the sound source and receiver depths to obtain the depth function 2.
゜(Z S r Z *).
固有関数計算手段6は、前記固有値実部計算手段3から
送られてくる固有値実部及び回転点深度の値、及び深度
関数生成手段5から与えられる音源又は受波器深度の値
から与えられた深度における固有関数値U、(z)を計
算するものである。The eigenfunction calculation means 6 is given from the eigenvalue real part and rotation point depth values sent from the eigenvalue real part calculation means 3 and the sound source or receiver depth value given from the depth function generation means 5. The eigenfunction value U, (z) at depth is calculated.
固有値虚部計算手段7は、前記初期設定手段2から送ら
れる音速プロファイル、音波の減衰係数、海面・海底反
射係数、前記固有値実部計算手段3から送られてくるノ
ーマルモード固有値の実部の値からノーマルモード固有
値の虚部を算出するものである。The eigenvalue imaginary part calculation means 7 calculates the sound speed profile, the acoustic wave attenuation coefficient, the sea surface/bottom reflection coefficient sent from the initial setting means 2, and the value of the real part of the normal mode eigenvalue sent from the eigenvalue real part calculation means 3. The imaginary part of the normal mode eigenvalue is calculated from
距離関数生成手段8は、前記固有値実部計算手段3及び
固有値虚部計算手段7から送られる固有値虚部ノ値から
ノーマルモード波の距離関数の各距離サンプルにおける
値を計算するものである。The distance function generating means 8 calculates the value of the distance function of the normal mode wave at each distance sample from the eigenvalue imaginary part values sent from the eigenvalue real part calculating means 3 and the eigenvalue imaginary part calculating means 7.
音場累積手段9は、前記深度関数生成手段5から送られ
る深度関数値を前記距離関数生成手段8から送られてく
る距離関数のサンプル値に乗じてノーマルモード波の音
波を算出し、すべてのノーマルモード波の寄与を各サン
プル点について累積していくものである。The sound field accumulation means 9 calculates the sound waves of normal mode waves by multiplying the depth function value sent from the depth function generation means 5 by the sample value of the distance function sent from the distance function generation means 8. The contribution of normal mode waves is accumulated for each sample point.
出力端子10は、前記音場累積手段9で得られた音場計
算結果を出力する出力端子である。The output terminal 10 is an output terminal for outputting the sound field calculation result obtained by the sound field accumulation means 9.
次に上記構成の装置の動作原理を説明する。初期設定手
段2は、入力端子1から環境条件等の計算条件を受ける
と音速プロファイルに(24)式の形の関数近似を施し
、(26)式右辺の7、=(πc、)”川ω/2g+l
”“、 j=1 。Next, the operating principle of the apparatus having the above configuration will be explained. When the initial setting means 2 receives calculation conditions such as environmental conditions from the input terminal 1, it applies a function approximation in the form of equation (24) to the sound speed profile, and calculates 7, = (πc,)" on the right side of equation (26). /2g+l
"", j=1.
2、 ・・・・、N ・・・・(
29)の計算を行なう。また、入力端子1から送られた
周波数fの音波の海面・海底反射係数のグレージング角
特性、海水中の音波の減衰係数、音波、受波器深度にお
ける波数値の2乗ks” = k”(z s)及びに*
”=k”(z*)を算出する。2,...,N...(
29). In addition, the grazing angle characteristics of the sea surface/bottom reflection coefficient of the sound wave of frequency f sent from input terminal 1, the attenuation coefficient of the sound wave in seawater, the square of the wave number at the receiver depth, ks" = k" ( z s) and *
"=k" (z*) is calculated.
初期設定手段2は上記計算を終えると、音速プロファイ
ルを固有値実部計算手段3、位相積分計箕手段4、深度
関数生成手段5、固有関数計算手段6及び固有値虚部計
算手段7に以下のようにデータを送る。即ち、
音源、受波点に関する情報を深度関数生成手段5に、
海底反射における位相変化特性を固有値実部計算手段3
及び深度関数生成手段5に、
海面反射係数、音波の減衰係数を固有値虚部計算手段7
に、
(29)式のパラメータ71.J=1.2.・・・、N
を固有関数計算手段6に、
そして、角周波数ω=2πrを固有値実部計算手段3、
位相積分計算手段4、深度関数生成手段5及び距離関数
生成手段8に、
また、ω1′2を固有関数固有関数計算手段6に送る。After completing the above calculation, the initial setting means 2 outputs the sound velocity profile to the eigenvalue real part calculation means 3, the phase integrator measurement means 4, the depth function generation means 5, the eigenfunction calculation means 6, and the eigenvalue imaginary part calculation means 7 as follows. Send data to. That is, the information regarding the sound source and the receiving point is sent to the depth function generation means 5, and the phase change characteristics in the seabed reflection is sent to the eigenvalue real part calculation means 3.
and the depth function generation means 5, and the eigenvalue imaginary part calculation means 7 to calculate the sea surface reflection coefficient and the attenuation coefficient of the sound wave.
Then, parameter 71 of equation (29). J=1.2. ..., N
to the eigenfunction calculation means 6, and the angular frequency ω=2πr to the eigenvalue real part calculation means 3,
It sends ω1'2 to the phase integral calculation means 4, the depth function generation means 5, and the distance function generation means 8, and also to the eigenfunction eigenfunction calculation means 6.
固有値実部計算手段3では、初期設定手段2から送られ
た音速プロファイル情報、海底反射における位相変化特
性及び角周波数ωを用いて、各ノーマル番号nに対応し
た固有値の実部ξ4を順に求めていく。これは上記(4
)を満足するξ。The eigenvalue real part calculating means 3 sequentially calculates the real part ξ4 of the eigenvalue corresponding to each normal number n using the sound speed profile information sent from the initial setting means 2, the phase change characteristics in the seabed reflection, and the angular frequency ω. go. This is the above (4)
).
を繰り返し法により求めることにより得られる。can be obtained by finding it by an iterative method.
(4)式中の積分値
一ξ。Z ) l / 2 dz ・・・
・(30)は、位相積分計算手段4によって実行される
。(4) Integral value 1ξ in formula. Z ) l / 2 dz...
- (30) is executed by the phase integral calculation means 4.
固有値実部計算手段3は固有値実部ξ。を求めめると、
これを深度関数生成手段5、固有値虚部計算手段7及び
距離関数生成手段8に、またξ、′を固有関数計算手段
6に送る。また、このとき回転点深度2゜+Zdの値を
深度関数生成手段5、固有関数計算手段6及び固有値虚
部計算手段7に送る。The eigenvalue real part calculation means 3 calculates the eigenvalue real part ξ. When you ask for
This is sent to the depth function generation means 5, the eigenvalue imaginary part calculation means 7, and the distance function generation means 8, and ξ,' is sent to the eigenfunction calculation means 6. Further, at this time, the value of the rotation point depth 2°+Zd is sent to the depth function generation means 5, the eigenfunction calculation means 6, and the eigenvalue imaginary part calculation means 7.
深度関数生成手段5は、固有値実部計算手段3から送ら
れる固有値実部のξ。及び初期設定手段2から送られる
音速プロファイル、角周波数ω及び海底反射特性から深
度関数の正規化定数A。の計算を行なう。定数A、は例
えば次式の積分を実行することにより算出できる。The depth function generation means 5 receives the eigenvalue real part ξ sent from the eigenvalue real part calculation means 3. and a normalization constant A of the depth function from the sound speed profile, angular frequency ω, and seabed reflection characteristics sent from the initial setting means 2. Perform the calculation. The constant A can be calculated, for example, by performing the integration of the following equation.
U、(z)は(9)式又は(10)式のWKB近似形を
用いても問題は生じない。また1、l1llは次式によ
り与えられる。No problem will occur even if the WKB approximate form of equation (9) or equation (10) is used for U and (z). Further, 1 and l1ll are given by the following equation.
・・・・(32) ここで、Rは海底面における音波の反射係数である。...(32) Here, R is the reflection coefficient of the sound wave on the seafloor surface.
深度関数生成手段5は、更に固有関数計算手段6に音源
深度2.及び音源深度2.における波数の2乗値ks!
を送って音源深度における固有関数U。<z=)を、ま
た受波器深度z、l及び受波器深度2゜における波数の
2乗値に*2を送って受波器深度における固有関数U、
(Zi)を、それぞれ固有関数計算手段6で計算させ、
深度関数2゜(Z5. Zi)=π1AjJ。(z、)
U。(2え〉を算出して、その結果を音場累積手段9に
送る。The depth function generation means 5 further provides the eigenfunction calculation means 6 with the sound source depth 2. and sound source depth 2. The square value of the wave number at ks!
and the eigenfunction U at the source depth. <z=), and send *2 to the square value of the wave number at the receiver depths z, l and receiver depth 2° to obtain the eigenfunction U at the receiver depth,
(Zi) are respectively calculated by the eigenfunction calculation means 6,
Depth function 2° (Z5. Zi) = π1AjJ. (z,)
U. (2e) and sends the result to the sound field accumulation means 9.
固有関数計算手段6は、深度関数生成手段5から深度2
を、また固有値実部計算手段3からξ。、z、、Z、を
受けると、zuから2までの位相積分値φ、及び2から
2.までの位相積分i4をそれぞれ位相積分計算手段4
で計算させ、φ、とφ4の値を用いて、(23)式又は
(26)式により固有関数U。(z)を計算する。The eigenfunction calculation means 6 receives the depth 2 from the depth function generation means 5.
and ξ from the eigenvalue real part calculation means 3. ,z, ,Z, the phase integral value φ from zu to 2, and from 2 to 2. The phase integral calculation means 4 calculates the phase integral i4 up to
Using the values of φ and φ4, the eigenfunction U is calculated by equation (23) or equation (26). Calculate (z).
固有関数計算手段6は、固有関数U、(z)の計算を終
えると、その結果を深度関数生成手段5に送る。When the eigenfunction calculation means 6 finishes calculating the eigenfunctions U and (z), it sends the results to the depth function generation means 5.
一方、固有値虚部計算手段7は、固有関数計算手段3か
ら固有値実部ξ。及び回転点深度2゜。On the other hand, the eigenvalue imaginary part calculation means 7 receives the eigenvalue real part ξ from the eigenvalue calculation means 3. and rotation point depth 2°.
2、を送られ、初期設定手段2から送られた音速プロフ
ァイル等を用いて(5)式にノーマルモードの等側音線
のサイクル距離り、を計算する。さらに海面、海底反射
が反射する場合はそれぞれ反射率IR5(ξ。)1及び
IR8(ξ。)1を求めて(反射せず屈折する場合は1
とする)、(6)式により固有値虚部α、を計算し、そ
の結果を距離関数生成手段8に送る。2, and using the sound velocity profile etc. sent from the initial setting means 2, the cycle distance of the isolateral sound ray in the normal mode is calculated using equation (5). Furthermore, if the sea surface or seabed is reflected, calculate the reflectance IR5 (ξ.) 1 and IR8 (ξ.) 1, respectively (if it is not reflected but refracted, 1
), the imaginary part α of the eigenvalue is calculated using equation (6), and the result is sent to the distance function generating means 8.
距離関数生成手段8では、固有値実部計算手段3から固
有値実部ξ。を、また固有値虚部計算手段7から固有値
虚部α、を受けると(3)式により複数固有値k。を求
め、所望の距離r、、、 m”1.2.・・・・9Mに
対する距離関数ベクトルpfiT=(ρ。1.ρ。2.
・・・・、ρ。M)(ρ。、=H,(1)(H。The distance function generation means 8 receives the eigenvalue real part ξ from the eigenvalue real part calculation means 3. , and receiving the imaginary part α of the eigenvalue from the imaginary eigenvalue calculation means 7, the imaginary part α of the eigenvalue is obtained from the equation (3). Find the distance function vector pfiT=(ρ.1.ρ.2.
..., ρ. M)(ρ.,=H,(1)(H.
(1) (k、 r 、 )を生成する。k、r、>>
1の場合には、H0(1)(k、r 、、)は次の漸近
形で近似することができる。(1) Generate (k, r, ). k, r, >>
1, H0(1)(k, r , ,) can be approximated by the following asymptotic form.
H、c′)(k、r 、、): (2/ yr k、r
、)exp (i(k、r、−π/4)〕
距離関数生成手段8は、距離関数値の組みρ。、。H, c′) (k, r , , ): (2/ yr k, r
, ) exp (i(k, r, -π/4))] The distance function generating means 8 generates a set of distance function values ρ.
m=1.2.・・・・9Mを生成すると、これを音場累
積手段9に送る。m=1.2. ...9M is generated and sent to the sound field accumulation means 9.
音場累積手段9では、深度関数生成手段5から深度関数
(zs、zえ)をまた距離関数生成手段8から距離関数
ベクトルρを送られると両者の積Z。The sound field accumulation means 9 receives the depth functions (zs, ze) from the depth function generation means 5 and the distance function vector ρ from the distance function generation means 8, and calculates the product Z of both.
p。であられされるモード番号nのノーマルモード音場
を計算し、それ以前のノーマルモード音場の和に累加し
ていく。p. The normal mode sound field of mode number n is calculated and added to the sum of the previous normal mode sound fields.
こうして全てのノーマルモード音場についての累加を終
えると、音場累積手段9はその結果1<”zs、z次)
t)(z、、z、)=Σza(zs、z*)ρ ・・
・・(33)れ
を出力端子10から出力する。After completing the accumulation for all normal mode sound fields, the sound field accumulation means 9 obtains the result 1<"zs, z-order) t) (z,, z,) = Σza (zs, z*) ρ...
...(33) Output it from the output terminal 10.
次に、固有関数計算手段6の実施例として、(21)式
のq(η)及び(26)式のAi(η)を次
式のように近似した場合の実施例について説明する。Next, as an example of the eigenfunction calculation means 6, an example will be described in which q(η) of equation (21) and Ai(η) of equation (26) are approximated as shown in the following equation.
q(η)=α、+β、η 、η、くη≦ηに+1 。q(η) = α, +β, η, η, +1 for η≦η.
k=0.・・・・、に、−1・・・・(34)AI(η
)−αに+βにη、ηk〈ηくηに+1+に寓にゆ、・
・・・、k ・・・・(35〉即ち、η=η
3゜を境界としてη≦η、。の範囲では(23)式を、
η〉η、。の範囲では(26)式を用いることにし、そ
れぞれq(η)及びA、(η)をα、+β、ηの形で折
れ線近似する。ただしηが負の十分大きな値になるとq
(77):1となるのでα。−1,β=Oとし、η1に
対応して
φ、、、、−2/3(−η1〉8′2 ・・・
・(36)を定めておく。また、ηが正の十分大きな値
をとる場合には、A、(η)キ0となるので、α、=β
、=Oとする。この場合−は虚数となる。k=0. ......, -1... (34) AI(η
) −α to +β to η, ηk〈η to η to +1+ to Fableniyu,・
..., k ... (35〉, that is, η=η
η≦η, with 3° as the boundary. In the range of , formula (23) becomes
η〉η,. In the range of , Equation (26) is used, and q(η) and A, (η) are approximated by polygonal lines in the form of α, +β, η, respectively. However, when η becomes a sufficiently large negative value, q
(77):1, so α. -1, β=O, corresponding to η1, φ, , , -2/3(-η1〉8'2...
- Define (36). Also, if η takes a sufficiently large positive value, A, (η) becomes 0, so α, = β
,=O. In this case, - becomes an imaginary number.
第3図は固有関数計算手段6の具体例を示すブロック図
である。同図において、各番号は下記部分を示す。FIG. 3 is a block diagram showing a specific example of the eigenfunction calculation means 6. In the figure, each number indicates the following parts.
20.22.28は、それぞれ前記初期設定手段2から
送られてくる海底の境界番号j++=N+1、(29)
式の71 + j” 1 + 2 、・・・・、N及び
角周波数ωの1/2乗を記憶するレジスタ、12.44
.14は、それぞれ前記固有値実部計算手段3から送ら
れてくる固有値の2乗ξ。2゜回転点深度Zu+Zd及
び回転点深度の層番号j1.j4を記憶するレジスタ、
13.15は、それぞれ前記位相積分計算手段4から送
られてくるZ。からZまでの位相積分値φ、及び2から
z6までの位相積分値i、を記憶するレジスタ、
11.45は、それぞれ前記深度関数生成手段5から送
られてくるkSi又はkIIffi及び2.又は2゜を
記憶するレジスタ、
16.18.40は、加算器、
31.32,37.43は、乗算器、
17は、正弦関数計算における位相補正値を格納するレ
ジスタ、
19は、海面の層番号(j=0と定義)を格納するレジ
スタ、
21は、位相の基準値φ1&、を格納するレジスタ、
23は、−1/4乗算器、
24.27は、切換スイッチ、
25.26は、それぞれ海面と上側回転点の層番号及び
海底と下側回転点の層番号の比較に基つき前記切換スイ
ッチ24.27の制御信号を出力する比較器、
29は、正弦関数発生器、
30は、入力きれた2つを比較し、その最小値及び最小
値を、取る入力番号を出力する比較器、33は、−(3
/2(・)]””演算器、34は、q(η)の折れ線近
似の境界値を格納するレジスタ、
35は、前記−[3/2(・))”演算器33の出力と
前記折れ線近似の境界値を格納するレジスタ34の内容
とを比較し、用いるべき近似関数の番号を決定する近似
関数判定手段、
36は、q(η)とAI(η)との切換えの境界値η、
。の番号に0を格納するレジスタ、39は前記レジスタ
36の内容と前記近似関数判定手段35の出力を比較し
、U、、(z)の計算を(23)式により行なうか(2
6)式により行なうかを決定する比較器、
38.41は、それぞれq(η)、AI(η)の近似関
数の係数値β、及びα、を格納するレジスタ、42は、
前記比較器39の制御に基づき前記乗算器32の出力と
レジスタ22の出力との切換えを行なう切換スイッチ、
46は固有関数U、(z)の出力端子、次に上記構成の
固有関数計算手段の動作原理を説明する。20.22.28 are the seabed boundary numbers j++=N+1, (29) sent from the initial setting means 2, respectively.
71 + j" 1 + 2, ..., N of the formula and a register for storing the 1/2 power of the angular frequency ω, 12.44
.. 14 are the squares ξ of the eigenvalues sent from the eigenvalue real part calculation means 3, respectively. 2° rotation point depth Zu+Zd and layer number j1 of the rotation point depth. Registers 13 and 15 that store j4 are Z sent from the phase integral calculation means 4, respectively. A register 11.45 stores the phase integral value φ from 2 to Z and the phase integral value i from 2 to z6, respectively. 16.18.40 is an adder; 31.32, 37.43 are multipliers; 17 is a register that stores a phase correction value in sine function calculation; 21 is a register that stores the phase reference value φ1&, 23 is a -1/4 multiplier, 24.27 is a changeover switch, 25.26 is a register that stores the layer number (defined as j = 0), , a comparator that outputs a control signal for the changeover switches 24 and 27 based on a comparison between the layer numbers of the sea surface and the upper rotation point, and the layer numbers of the seabed and the lower rotation point, respectively; 29 is a sine function generator; 30 is a sine function generator; , the comparator 33 outputs the input number which compares the two inputs and takes the minimum value and the minimum value, -(3
34 is a register that stores the boundary value of the polygonal line approximation of q(η); 35 is a register that stores the output of the −[3/2(・))” operator 33 and the above-mentioned Approximation function determining means for determining the number of the approximation function to be used by comparing the contents of a register 34 that stores the boundary value of the polygonal line approximation; 36 is the boundary value η for switching between q(η) and AI(η); ,
. A register 39 stores 0 in the number of , and compares the contents of the register 36 with the output of the approximate function determining means 35, and calculates U, , (z) by equation (23) or (2).
6) A comparator that determines whether to perform the process according to the formula; 38.41 is a register that stores coefficient values β and α of the approximation functions of q(η) and AI(η), respectively; 42 is a register that stores
a changeover switch for switching between the output of the multiplier 32 and the output of the register 22 based on the control of the comparator 39; 46 is an output terminal of the eigenfunction U, (z); The operating principle will be explained.
まず、初期設定手段2から送られてくる海底の層番号j
s及び角周波数ωの平方根ω1″がそれぞれレジスタ2
0及び28に格納される。また、固有値実部計算手段3
から送られた固有値の2乗ξI9回転点深度Z、、Za
及び回転点深度の層番号j、、Lがそれぞれレジスタ1
2,44.14に格納される。First, the seabed layer number j sent from the initial setting means 2
s and the square root ω1″ of the angular frequency ω are respectively stored in register 2.
0 and 28. In addition, the eigenvalue real part calculation means 3
The square of the eigenvalue ξI9 rotation point depth Z, , Za
and the layer numbers j, , L of the rotation point depth are respectively registered in register 1.
2,44.14.
次に、深度関数生成手段5から深度2及び深度における
波数の2乗k !(z )を送られると、それぞれレジ
スタ45及び11に格納し、レジスタ44の2.とレジ
スタ45の2とを位相積分計算手段4に送って位相積分
値φ。を受け、レジスタ13に格納する。Next, from the depth function generating means 5, the square of the wave number at depth 2 and the depth k! When (z) is sent, it is stored in registers 45 and 11, respectively, and 2. and 2 of the register 45 are sent to the phase integral calculation means 4 to obtain the phase integral value φ. is received and stored in the register 13.
次に比較器25は、レジスタ14の3段レジスタ19の
内容0とを比較し、j、=0の場合にはレジスタ14か
ら0を、またju≠0の場合にはレジスタ14からπ/
4を読出して加算器18に送る。きらに比較器25は、
j。=0の場合にはレジスタ21の内容φrmaxが、
またju≠0の場合にはレジスタ13の内容φ。が比較
器30に送られるように、切換スイッチ24を制御し、
該切換スイッチ24の出力を−tとする。Next, the comparator 25 compares the contents of the register 14 with the contents 0 of the three-stage register 19, and reads 0 from the register 14 if j = 0, and reads π// from the register 14 if ju≠0.
4 is read out and sent to the adder 18. Kirani comparator 25 is
j. = 0, the content φrmax of the register 21 is
If ju≠0, the contents of the register 13 are φ. is sent to the comparator 30,
The output of the changeover switch 24 is assumed to be -t.
加算器18は、レジスタ13の内容φ、にレジスタ17
から送られた値を加算し、正弦関数発生器29におくる
。該正弦関数発生器29は、加算器18から送られてき
た位相値の正弦(−啜め、その結果を乗算器32に送る
。The adder 18 adds the contents φ of the register 13 to the register 17
The values sent from are added and sent to the sine function generator 29. The sine function generator 29 subtracts the sine value of the phase value sent from the adder 18 and sends the result to the multiplier 32 .
次に、レジスタ45の2としレジスタ44の26とが位
相積分計算手段4に送ら′れ、該位相積分計算手段4で
位相積分値φ、が計算されてレジスタ15に格納される
と、比較器26はレジスタ14の内容jdとレジスタ2
0の内容j、とを比較し、Ja−Jllの場合には、レ
ジスタ21の内容φ0.工を、またJ、+≠jsの場合
には、レジスタ15の内容φ6を比較器30に送るよう
に切換スイッチ27を制御し、切換スイッチ27の出力
をφ4′とする。比較器30では、切換スイッチ24及
び27からそれぞれ位相データを受は取ると、両者を比
較し、次の判断に基づいて一方を選択する。■φu’l
d′が共に実数の場合には値の小さい方を選択する。Next, 2 of the register 45 and 26 of the register 44 are sent to the phase integral calculation means 4, and when the phase integral value φ is calculated by the phase integral calculation means 4 and stored in the register 15, the comparator 26 is the contents jd of register 14 and register 2
0 contents j, and in the case of Ja-Jll, the contents φ0. If J, +≠js, the changeover switch 27 is controlled to send the contents φ6 of the register 15 to the comparator 30, and the output of the changeover switch 27 is set to φ4'. When the comparator 30 receives the phase data from the changeover switches 24 and 27, it compares them and selects one based on the next judgment. ■φu'l
If both d' are real numbers, the one with the smaller value is selected.
■φ、°、φ4゛の一方が実数、他方が虚数の場合には
虚数の方を選択する。■If one of φ, °, and φ4′ is a real number and the other is an imaginary number, select the imaginary number.
■φ8′、φ4′がともに虚数の場合には絶対値の大き
い方を選択する。(2) If both φ8' and φ4' are imaginary numbers, select the one with the larger absolute value.
比較器30はφu’+’bd′の一方を選択するとその
値を−(3/2(・) ) ””演算器33に送ると共
に、選んだφ゛に対応した回転点の層番号の値をレジス
タ14から読出して、その番号に対応した7、の値をレ
ジスタ群22から読出し、切換スイッチ42に送る。When the comparator 30 selects one of φu'+'bd', it sends that value to the -(3/2(・))" calculation unit 33, and also outputs the value of the layer number of the rotation point corresponding to the selected φ". is read from the register 14, and the value 7 corresponding to that number is read from the register group 22 and sent to the changeover switch 42.
−(3/2(・))””演算器33は比較器30から送
られた位相値に対応に対して、次式7式%(37)
によりηの計算を行なう。()””は3つの解を持つが
、実数解をとることにする。−〔3/2(・) ) ”
”演算器33はηを算出すると、ηの値を近似関数判定
手段35及び乗算器37に送る。-(3/2(·))"" The computing unit 33 calculates η according to the following equation 7, %(37), corresponding to the phase value sent from the comparator 30. ()"" has three solutions, but we will take the real number solution. −[3/2(・) )”
After calculating η, the calculator 33 sends the value of η to the approximate function determining means 35 and the multiplier 37.
近似関数判定手段35は−(3/2(・) :] ”’
演算器33かも送られた刀の値とレジスタ34に格納さ
れた関数近似の境界値と比較し、η、くη≦η、1
・・・・(38)となるkを選び出し
て、その番号を比較器39に送るとともにkに対応した
近似関数の係数β、をレジスタ38から選び出して乗算
器37に、そして係数α、をレジスタ41から選び出し
て加算器40に送る。The approximate function determining means 35 is -(3/2(・):] ”'
The arithmetic unit 33 also compares the sent sword value with the boundary value of the function approximation stored in the register 34, and calculates η, η≦η, 1
...(38), and send that number to the comparator 39, select the coefficient β of the approximation function corresponding to k from the register 38 and send it to the multiplier 37, and send the coefficient α to the register. 41 and sends it to the adder 40.
乗算器37は前記−(3/2(・) ) ”’演算器3
3から送られたηにレジスタ38から送られたβ、を乗
し、その結果βに77を加算器40に送る。The multiplier 37 is the above-mentioned −(3/2(・)) ”' arithmetic unit 3
.eta. sent from register 38 is multiplied by .beta. sent from register 38, and 77 is sent to the adder 40 as a result.
加算器40は乗算器37からβ、ηの値を送られると、
これにレジスタ41から送られたαを加え、その結果を
乗算器43に送る。When the adder 40 receives the values of β and η from the multiplier 37,
α sent from the register 41 is added to this, and the result is sent to the multiplier 43.
比較器39は、近似関数判定手段35から送られた近似
関数の番号にと、レジスタ36に格納されたq(η)、
AI(77>との境界の番号に0とを比較し、k <
k oの場合は加算器16に起動信号を送ると共に、乗
算器32の出力が乗算器43に行くように、切換スイッ
チ42を制御する。またに≧に0の場合にはレジスタ2
2の出力が乗算器43に行くように切換スイッチ42を
制御する。The comparator 39 calculates the approximate function number sent from the approximate function determining means 35, q(η) stored in the register 36,
Compare the number of the boundary with AI(77> with 0 and find k <
In the case of k o, a start signal is sent to the adder 16, and the changeover switch 42 is controlled so that the output of the multiplier 32 goes to the multiplier 43. Also, if ≧ is 0, register 2
The changeover switch 42 is controlled so that the output of the output signal 2 goes to the multiplier 43.
加算器16は比較器39から起動信号を受けると、レジ
スタ11に格納された波数の2乗値に!(2)から、レ
ジスタ12に格納きれた固有値の2乗ξ。2を減じ、そ
の結果k”(z)−ξ。′を一1/4乗算器23に送る
。−1/4乗算器23は加算器16からk”(z)−ξ
、、!の値を受けると、これを−1/4乗し、その結果
を乗算器31に送る。When the adder 16 receives the activation signal from the comparator 39, it becomes the square value of the wave number stored in the register 11! From (2), the square of the eigenvalue stored in the register 12 is ξ. 2 and sends the result k"(z)-ξ.' to the 1/4 multiplier 23. The -1/4 multiplier 23 subtracts k"(z)-ξ.
,,! When it receives the value, it raises it to the -1/4 power and sends the result to the multiplier 31.
乗算器31は一1/4乗算器23から(kl(2)−ξ
、2 ) 1/4の値を受けると、これにレジスタ28
に格納されたω1″の値を乗じ、その結果を乗算器32
に送る。乗算器32は乗算器31から送られたωl /
l (kl (z )−ξ% 〕−1/4の値に、正
弦関数発生器29から送られた正弦関数値を乗し、その
結果を切換スイッチ42に送る。The multiplier 31 receives (kl(2)-ξ
, 2) When the value of 1/4 is received, register 28 is added to it.
is multiplied by the value of ω1″ stored in
send to Multiplier 32 receives ωl/
The value of l(kl(z)-ξ%)-1/4 is multiplied by the sine function value sent from the sine function generator 29, and the result is sent to the changeover switch 42.
切換スイッチ42は比較器39の制御に基づき、乗算器
32から送られる値、又はレジスタ22から送られる値
のいずれかを通過させて乗算器43に送る。乗算器43
は切換スイッチ42から送られる値に加算器40から送
られる近似関数値を乗じ、その結果である固有関数値U
、(z)を出力端子46から出力する。The changeover switch 42 passes either the value sent from the multiplier 32 or the value sent from the register 22 and sends it to the multiplier 43 under the control of the comparator 39 . Multiplier 43
multiplies the value sent from the changeover switch 42 by the approximate function value sent from the adder 40, and the resulting eigenfunction value U
, (z) are output from the output terminal 46.
以上説明したように本発明に係る音波伝搬シミュレーシ
ョン方式によれば、固有関数を回転点付近でも発散した
り不安定になったりすることがない安定した精度の良い
関数で近似しているので、従来の方式のように回転点近
傍で精度が劣化し、値が発散してしまい計算不能になる
という問題点を除くことができ、特に低周波において計
算精度を大幅に向上きせることができるという優れた効
果が得られる。As explained above, according to the sound wave propagation simulation method according to the present invention, the eigenfunction is approximated by a stable and accurate function that does not diverge or become unstable even near the rotation point. This method has the advantage of being able to significantly improve calculation accuracy, especially at low frequencies, by eliminating the problem that accuracy deteriorates near the rotation point and the values diverge, making calculations impossible. Effects can be obtained.
また、本発明に係る音波伝搬シミュレーション方式は固
有関数の近似を非常に簡単な処理により実現しているの
で、従来の方式に比較して処理時間の増加は殆どなくソ
ーナーの運用現場での性能評価等、実時間性を要求され
るシステムに利用することかできるという優れた効果も
有する。In addition, since the sound wave propagation simulation method according to the present invention realizes the approximation of the eigenfunction through extremely simple processing, there is almost no increase in processing time compared to conventional methods, and performance evaluation in the field of sonar operation is possible. It also has the excellent effect of being able to be used in systems that require real-time performance.
第1図は本発明に係る音波伝搬シミュレーション方式を
実施する装置の構成を示すブロック図、第2図は関数q
(η)の振舞を示す説明図、第3第1図の固有関数計算
手段の具体的構成例を示す図である。
図中、1・・・・入力端子、2・・・・初期設定手段、
3・・・・固有値実部計算手段、4・・・・位相積分計
算手段、5・・・・深度関数生成手段、6・・・・固有
関数計算手段、7・・・・固有値虚部計算手段、8・・
・・距離関数生成手段、9・・・・音場累積手段、10
,46・・・・出力端子、11 、12 、13 、1
5 、17.19,20,21,22,28,34.3
B、38,39,44,45.・・・・レジスタ、16
.18.20・・・・加算器、23・・・・−1/4乗
算器、24,27.42・・・・切換スイッチ、25.
26,30.39・・・・比較器、29・・・・正弦関
数発生器、31,32.43・・・・乗算器、33・・
・・−(3/2(・))”3演算器、35・・・・近似
関数判定手段。FIG. 1 is a block diagram showing the configuration of a device that implements the sound wave propagation simulation method according to the present invention, and FIG.
FIG. 3 is an explanatory diagram showing the behavior of (η); FIG. 3 is a diagram showing a specific configuration example of the eigenfunction calculation means in FIG. 1; In the figure, 1: input terminal, 2: initial setting means,
3...Eigenvalue real part calculation means, 4...Phase integral calculation means, 5...Depth function generation means, 6...Eigenfunction calculation means, 7...Eigenvalue imaginary part calculation Means, 8...
...Distance function generation means, 9...Sound field accumulation means, 10
,46...Output terminal, 11, 12, 13, 1
5, 17.19, 20, 21, 22, 28, 34.3
B, 38, 39, 44, 45. ...Register, 16
.. 18.20... Adder, 23...-1/4 multiplier, 24, 27.42... Changeover switch, 25.
26, 30.39... Comparator, 29... Sine function generator, 31, 32.43... Multiplier, 33...
...-(3/2(-))"3 arithmetic units, 35...Approximate function determination means.
Claims (1)
マルモード論理に基づく音波伝搬シミュレーション方式
において、 音源深度及び受波点深度におけるノーマルモード波の固
有関数の振幅値をノーマルモード波の回転点深度から音
源深度及び受波点深度までの位相積分値の低次関数を用
いて表わす固有関数計算手段を具備し、該固有関数計算
手段によって得られる音源深度及び受波点深度における
ノーマルモード関数値を用いて音場計算を行なうことを
特徴とする音波伝搬シミュレーション方式。[Claims] In a sound wave propagation simulation method based on normal mode logic that simulates a sound wave propagation phenomenon in a heterogeneous medium, the amplitude value of the eigenfunction of a normal mode wave at a sound source depth and a receiving point depth is defined as a normal mode. Eigenfunction calculation means is provided, which is expressed using a low-order function of the phase integral value from the rotation point depth of the wave to the sound source depth and the wave reception point depth, and A sound wave propagation simulation method characterized by performing sound field calculations using normal mode function values.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP11772787A JPH0746137B2 (en) | 1987-05-14 | 1987-05-14 | Sound wave propagation simulation method |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP11772787A JPH0746137B2 (en) | 1987-05-14 | 1987-05-14 | Sound wave propagation simulation method |
Publications (2)
Publication Number | Publication Date |
---|---|
JPS63282679A true JPS63282679A (en) | 1988-11-18 |
JPH0746137B2 JPH0746137B2 (en) | 1995-05-17 |
Family
ID=14718791
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP11772787A Expired - Lifetime JPH0746137B2 (en) | 1987-05-14 | 1987-05-14 | Sound wave propagation simulation method |
Country Status (1)
Country | Link |
---|---|
JP (1) | JPH0746137B2 (en) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114355287A (en) * | 2022-01-04 | 2022-04-15 | 湖南大学 | Ultra-short baseline underwater acoustic ranging method and system |
CN115014703A (en) * | 2022-05-26 | 2022-09-06 | 中国人民解放军海军工程大学 | Method for forecasting seawater surface wave excited by midpoint sound source in actual ocean waveguide |
CN117249894A (en) * | 2023-11-16 | 2023-12-19 | 自然资源部第一海洋研究所 | Diagnosis method for transmission thickness of underwater far-field sound propagation on seabed |
-
1987
- 1987-05-14 JP JP11772787A patent/JPH0746137B2/en not_active Expired - Lifetime
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114355287A (en) * | 2022-01-04 | 2022-04-15 | 湖南大学 | Ultra-short baseline underwater acoustic ranging method and system |
CN114355287B (en) * | 2022-01-04 | 2023-08-15 | 湖南大学 | Ultra-short baseline underwater sound distance measurement method and system |
CN115014703A (en) * | 2022-05-26 | 2022-09-06 | 中国人民解放军海军工程大学 | Method for forecasting seawater surface wave excited by midpoint sound source in actual ocean waveguide |
CN115014703B (en) * | 2022-05-26 | 2024-06-07 | 中国人民解放军海军工程大学 | Sea surface wave forecasting method for excitation of point sound source in actual ocean waveguide |
CN117249894A (en) * | 2023-11-16 | 2023-12-19 | 自然资源部第一海洋研究所 | Diagnosis method for transmission thickness of underwater far-field sound propagation on seabed |
CN117249894B (en) * | 2023-11-16 | 2024-04-05 | 自然资源部第一海洋研究所 | Diagnosis method for transmission thickness of underwater far-field sound propagation on seabed |
Also Published As
Publication number | Publication date |
---|---|
JPH0746137B2 (en) | 1995-05-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US5947903A (en) | Method to estimate planar flow from Doppler velocity distribution in an observation plane | |
NO153021B (en) | DEVICE FOR MEASUREMENT OF WATER FLOW | |
Rudolph et al. | Doppler Velocity Log theory and preliminary considerations for design and construction | |
CN114779170A (en) | Shallow sea near-field sound source positioning method | |
CN110703205B (en) | Ultra-short baseline positioning method based on self-adaptive unscented Kalman filtering | |
JPS63282679A (en) | Sonic wave propagation simulation system | |
CN107907591B (en) | Ultrasonic detection system and method for component concentration of multi-component solid-liquid two-phase mixture | |
Baxley et al. | Comparison of beam tracing algorithms | |
Skobelt’syn et al. | Finding, by means of a scattered sound, the geometric parameters of a finite elastic cylinder located near the half-space border | |
JPH10170638A (en) | Underwater sound simulator | |
Reilly et al. | Computing acoustic transmission loss using 3D Gaussian ray bundles in geodetic coordinates | |
JP7243531B2 (en) | Distance estimation device, distance estimation method, and distance estimation program | |
Longbottom et al. | Ultrasonic multiple-sensor solid level measurements | |
US7920821B2 (en) | All ship speeds correlation SONAR simulator | |
JP3112158B2 (en) | Sound source level measuring method and sound source level measuring device | |
JPH052193B2 (en) | ||
Varadan et al. | Ultrasonic wave scattering by a subsurface flaw in joined fluid-solid half spaces | |
Štramberský et al. | Model of the Far-Field Acoustic Localisation | |
Gray et al. | Blade rate radiation from cavitation volume modulation on merchant vessels | |
Adegeest | Linear and nonlinear hydrodynamic forces acting on a segmented heaving Wigley model | |
Fox | Computation of linear ultrasound fields using 2D Fourier-Bessel series | |
Cooper et al. | Session Summary of “Parabolic Equation Methods Across Acoustics” | |
Loeser | Determination of Maneuvering Properties in Shallow Water by Impulse Response Techniques | |
Nahraoui et al. | Modeling the Cut-off Frequency of Acoustic Signal with an Adaptative Neuro-Fuzzy Inference System (ANFIS) | |
Bergem et al. | Comparison of a time domain spectral collocation method and the time-evolution backscatter model BORIS |