JPS63282679A - Sonic wave propagation simulation system - Google Patents

Abstract

PURPOSE:To realize a practical simulation system capable of realizing high accuracy, by performing the calculation of a sound field using the normal mode function values at the depth of a sound source and at the depth of a receiving point obtained by an inherent function calculation means. CONSTITUTION:An inherent function calculation means showing the amplitude values of the inherent functions of a normal mode wave at the depth of a sound source and at the depth of a receiving point using the low order functions of the position integrated values from the depth of the rotary point of the normal mode wave to the depth of the sound source and that of the receiving point is provided. Then, a sound field is calculated by using the normal mode function values at the depth of the sound source and at the depth of the receiving point obtained by said inherent function calculation means. By this method, since the inherent function is made approximate as a stable good accuracy function not diffused or becoming unstable in the vicinity of the rotary point, it is eliminated that a value is diffused by the deterioration of accuracy in the vicinity of the rotary point and calculation is rendered impossible.

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.

〔従来技術〕[Prior art]

ソーナー等の海洋音響機器の性能は、その運用現場の海
洋環境により大きく変化する。従って、機器の効果的運
用のために、海洋中の音波伝搬特性を把握することが重
要であり、種々の音波伝搬モデルが開発されている。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 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) 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) The sound pressure becomes O at sea level.

（２）０＜ｚ＜ｏｏの全領域でＵ、（ｚ）及びｄＵ、／
ｄｚはそれぞれ連続である。(2) U, (z) and dU, / in the entire range 0<z<oo
Each of dz is continuous.

（３ン２→■でＳｏｍｍｅｒｆｅｌｄの放射条件を満足
する。但し、海底特性が海底表面ｚ＝ｚｌｌにおける反
射係数の形で与えられている場合は、Ｕ。／（ｄＵ、、
／ｄｚ）ｚ＝ｚｌｌが設定された値をとる。(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.

波数の鉛直分布ｋ（Ｚ）Ｃ以下、これを波数プロファイ
ルと呼ぶ）は簡単な関数で近似されるが、通常はＯ＜ｚ
＜ωの全範囲を１つの関数で近似することが困難であり
、媒質をいくつかの層に分割して、各層内をそれぞれ別
の関数で表わすという方法がとられている。これに伴い
、前記固有関数Ｕ。（Ｚ）も各層でそれぞれ関数を用い
て表わされる。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.

しかし、層数が多い場合は、ノーマルモード固有（ａ　
ｋ　、を厳密に求めるには多大な計算時間とダイナミッ
ク・レンジが必要であり、そのためしばしばＷＫＢ近似
による固有値及び固有関数の近似計算が行なわれる。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.

固有値ｋ。を実部と虚部とに分け に、、＝ξｏ＋ｉα１　　　　　　　　　　・・・・（
３）とおくと、実部ξ。は次の方程式を解くことに近似
的に求められる。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 ... 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) 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).

固有関数Ｕ。（ｚ）はＷＫＢ近似は次式の形で与えられ
る。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[knee jp(Z,, z, ξa))...(
7) Here, p(z,,z,ξ,) is defined by the following equation.

・・・・く８） Ｂｔ、Ｂｔは境界条件により定まる定数である。...ku8) Bt and Bt are constants determined by boundary conditions.

等価音線が海面反射するようなモード波の場合、２．＝
０となり、Ｂｌ＝　１／２　ｊ　、Ｂ！−−１２ｊとお
くことができて、Ｕ、（ｚ　）はＵ、（ｚ）＝　　ω　
ｌ／　２　（ｋ　１　（ｚ　　）−ξ　、ｚ　　〕　−
”’ｓｉｎ［１）（０，ｚ、ξ、）〕　　　　　　　　
　・・・・（９）と簡単になる。また、等側音線が海面
に達せずａｅｘｐを形成するようなモードの場合には、
Ｂ１−ｅｘｐ（−ｊ　π／　４　）／　２、Ｂ　、＝　
−ｅｘｐ（ｊπ／４）／２とおくことができＵ、（ｚ）
は、 Ｕｎ（Ｚ　）＝　ω”２［ｋ　”（Ｚ　）−ξａ”）−
”’！Ｂ１ｎ［ｐ（ｚ、、ｚ、ξ。）＋π／４〕・・・
・（１０）となる。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).

〔発明が解決しようとする問題点〕[Problem that the invention seeks to solve]

上述のＷＫＢ近似法はノーマルモード解を厳密に求める
場合に比べて計算が非常に容易であり、また送受深度が
等側音線の上下端から十分離れた場合には良好な結果が
得られる。また、等側音線が海面や海底で反射している
場合には、送受波深度がその反射点に近づいても精度は
保たれるという特徴を有している。しかしながら等側音
線の上下端のいずれかが回転点となり、しかも送波深度
又は受波深度の少なくとも一方がその回転点付近にある
場合には、上記（９）、（１０）式の〔ｋ２（ｚ）−ξ
　ｇ　）−１／４の項が非常に大きな値になり、回転点
に近付くとともに無限大に発散してしまいノーマルモー
ド波の計算ができなくなってしまう。特に低周波では、
上述の振幅の異常が生じる深度範囲は広くなる。従って
、広い領域にわたって音圧分布に誤差が生じてしまうと
いう問題があった。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.

〔問題点を解決するための手段〕[Means for solving problems]

上記問題点を解決するため本発明は、音波伝搬シミュレ
ーション方式を、音源深度及び受波点深度におけるノー
マルモード波の固有関数の振幅値をノーマルモード波の
回転点深度から音源深度及び受波点深度までの位相積分
値の低次関数を用いて表わす固有関数計算手段を具備し
、該固有関数計算手段によって得られる音源深度及び受
波点深度におけるノーマルモード関数値を用いて音場計
算を行なうようにした。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.

〔作用〕[Effect]

上記の如く音波伝搬シミュレーション方式を構成するこ
とにより、固有関数を回転点付近でも発散したり不安定
になったりすることがない安定した精度の良い関数で近
似しているので、従来の方式のように回転点近傍で精度
が劣化し、値が発散してしまい計算不能になるという問
題点を除（ことができる。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.

〔実施例〕〔Example〕

以下、本発明の詳細な説明する。 The present invention will be explained in detail below.

まず、始めに本発明の詳細な説明する。上記（９）、（
１０）式による近似は、位相については回転点の近傍で
も十分な精度が保たれるで、振幅項についてのみ精度を
改善することを考える。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.

回転点深度ｚ＝ｚ、の近傍で、音速Ｃ（ｚ）＝ω／ｋ（
ｋ）が Ｃ−”（ｚ）＝Ｃ，−”〔１−２ｇ、（ｚ−ｚ、、）／
Ｃ，）・・・（１１） の形で近似できるものとする。ここで、Ｃ０＝ω／ξ。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.

次の置きかえ ａ３＝−２ｇ、ξｎ”／　Ｃａ、７７　＝　　ａ（ｚ　
　ｚｕ）・・・・（１２） を行なうと、（２）式は次のＡｉ　ｒｙの方程式に簡略
化される。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) 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) However, AI(77) and Bl(77) are Airy functions.

上記（１４）式は海面における境界条件を満たす必要が
あるが、海面は回転点から十分離れているとすると、β
１＝−βｚ＝　１　／　２　Ｊとおくことができ、Ｕ、
、（ｚ　）＝　（π（Ｌ）　／　ａ　）””Ａ　＋（７
７）−（１５）となる。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).

ここでａ、ηを別の形で表わすことを考える。Let us now consider expressing a and η in another form.

（１１）式を（８）式に代入すると、 Ｃ，）　””ｄ　ｚ　＝　２／　３　（−（２ξ、”ｇ
、）／Ｃ、）　””＜　ｚ　−ｚ　、）””　　　　　
・＝（１６）従って、（１２）、（１６）式により、次
式が得られる。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) Also, from equations (11) and (12), the following equation is obtained.

ａ−（ｋ”（ｚ）−ξ　＊　）　ｌ　／　！　／　（ｎ
　）　ｌ　／　ｔ・・・・（１８） （１８）式を用いると（１５）式は次のように書き改め
られる。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) 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).

しかしく１９）式中の（−η／（ｋ”（ｚ）−ξ　１〕
）ｌ／４（７）部分はξ。＝ｋ（ｚ）とすると不定にな
ってしまうため、別の形で表現することが必要である。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.

これは（１２）式を（１５）に代入すると容易に得られ
る。This can be easily obtained by substituting equation (12) into (15).

Ｕ、、（Ｚ）＝”（ＷＣ，）””Ｉω／　２　ｇ　−１
”’Ａ　＋（η）・・・・（２０） 従って、ηの値によって（１９）、（２０）式を使い分
けるようにすれば、固有関数Ｕ。（ｚ）を精度よく求め
ることができ、それを（１）式に代入れすることによっ
て、音場を高精度でシミュレートできる。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.

ここで、（１９）式の計算を更に簡略化することを考え
る。（１９）式と（１０）式との比をとり、次の量を定
義する。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 (φ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.

第２図かられかるように、ｑ（η）はηの非常に滑らか
な関数になっており、ηが負の大きな値になると１に漸
近していく。従っ工、ｑ（η〉を簡単な関数形で近似し
ておけば、（１０）式のＷＫＢ近似を次のように補正す
ることにより、簡単な計算で高精度の固有関数を得るこ
とができる。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) 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,)-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).

固有関数Ｕ。（ｚ）を（２６）式の形で扱わねばならな
いのはη−０の近傍の非常に狭い範囲である。この範囲
ではＡ、（η）は極めて滑らかな関数となり、１次式等
の低次式で十分近似することができる。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.

以上から、深度２が回転点深度ｚｕに十分近い場合には
（２６）式を用い、ある程度以上離れている場合には（
２３）式の補正式を用いることにすれば、従来の方法で
問題となった回転点近傍での精度劣化や計算異常（無限
大発散）の現象が解消され、簡単な計算で精度のよい固
有関数値を得ることができる。また（２６）式はη〉０
の領域にも適用できる。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.

音源又は受波器深度が下側の回転点（ｎａｄｉｒと呼ぶ
）深度ｚ４の近くにある場合にも、同様の方法で精度の
よい固有関数が得られる。この場合、φ、ηは次式で与
えられる。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 ,, 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.

第１図は、本発明に係る音波伝搬シミュレーション方式
を実施する装置の構成を示すブロック図である。図示す
るように該装置は、入力端子１、初期設定手段２、固有
値実部計算手段３、位相積分計算手段４、深度関数生成
手段５、固有関数計算手段６、固有値虚部計算手段７、
距離関数生成手段８、音場累積手段９及び出力端子１０
で構成されている。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.

初期設定手段２は、前記入力端子１に入力される環境条
件及び周波数等の入力データに基づき、音速プロファイ
ル、音波の減衰係数、海面・海底反射係数等を設定する
ものである。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.

固有値実部計算手段３は、前記初期設定手段２から送ら
れる音速プロファイル、海底反射係数、周波数等に基づ
きノーマルモード固有値の実部を計算するものである。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.

位相積分計算手段４は、前記初期設定手段２から送られ
てくる音速プロファイル、周波数及び前記固有値実部計
算手段３又は固有関数計算手段６から送られてくる固有
値実部及び積分の上下限深度の値から、位相積分値を算
出するものである。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.

深度関数生成手段５は、ノーマルモード波の正規化定数
Ａ。を計算し、更に固有関数計算手段６を用いて音源及
び受波器深度における固有関数値を計算して深度関数２
゜（Ｚ　Ｓ　ｒ　Ｚ　＊　）を生成するものである。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 *).

固有関数計算手段６は、前記固有値実部計算手段３から
送られてくる固有値実部及び回転点深度の値、及び深度
関数生成手段５から与えられる音源又は受波器深度の値
から与えられた深度における固有関数値Ｕ、（ｚ）を計
算するものである。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.

固有値虚部計算手段７は、前記初期設定手段２から送ら
れる音速プロファイル、音波の減衰係数、海面・海底反
射係数、前記固有値実部計算手段３から送られてくるノ
ーマルモード固有値の実部の値からノーマルモード固有
値の虚部を算出するものである。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

距離関数生成手段８は、前記固有値実部計算手段３及び
固有値虚部計算手段７から送られる固有値虚部ノ値から
ノーマルモード波の距離関数の各距離サンプルにおける
値を計算するものである。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.

音場累積手段９は、前記深度関数生成手段５から送られ
る深度関数値を前記距離関数生成手段８から送られてく
る距離関数のサンプル値に乗じてノーマルモード波の音
波を算出し、すべてのノーマルモード波の寄与を各サン
プル点について累積していくものである。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.

出力端子１０は、前記音場累積手段９で得られた音場計
算結果を出力する出力端子である。The output terminal 10 is an output terminal for outputting the sound field calculation result obtained by the sound field accumulation means 9.

次に上記構成の装置の動作原理を説明する。初期設定手
段２は、入力端子１から環境条件等の計算条件を受ける
と音速プロファイルに（２４）式の形の関数近似を施し
、（２６）式右辺の７、＝（πｃ、）”川ω／２ｇ＋ｌ
”“、　ｊ＝１　。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). 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.

初期設定手段２は上記計算を終えると、音速プロファイ
ルを固有値実部計算手段３、位相積分計箕手段４、深度
関数生成手段５、固有関数計算手段６及び固有値虚部計
算手段７に以下のようにデータを送る。即ち、 音源、受波点に関する情報を深度関数生成手段５に、 海底反射における位相変化特性を固有値実部計算手段３
及び深度関数生成手段５に、 海面反射係数、音波の減衰係数を固有値虚部計算手段７
に、 （２９）式のパラメータ７１．Ｊ＝１．２．・・・、Ｎ
を固有関数計算手段６に、 そして、角周波数ω＝２πｒを固有値実部計算手段３、
位相積分計算手段４、深度関数生成手段５及び距離関数
生成手段８に、 また、ω１′２を固有関数固有関数計算手段６に送る。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.

固有値実部計算手段３では、初期設定手段２から送られ
た音速プロファイル情報、海底反射における位相変化特
性及び角周波数ωを用いて、各ノーマル番号ｎに対応し
た固有値の実部ξ４を順に求めていく。これは上記（４
）を満足するξ。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) Integral value 1ξ in formula. Z ) l / 2 dz...
- (30) is executed by the phase integral calculation means 4.

固有値実部計算手段３は固有値実部ξ。を求めめると、
これを深度関数生成手段５、固有値虚部計算手段７及び
距離関数生成手段８に、またξ、′を固有関数計算手段
６に送る。また、このとき回転点深度２゜＋Ｚｄの値を
深度関数生成手段５、固有関数計算手段６及び固有値虚
部計算手段７に送る。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.

深度関数生成手段５は、固有値実部計算手段３から送ら
れる固有値実部のξ。及び初期設定手段２から送られる
音速プロファイル、角周波数ω及び海底反射特性から深
度関数の正規化定数Ａ。の計算を行なう。定数Ａ、は例
えば次式の積分を実行することにより算出できる。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.

Ｕ、（ｚ）は（９）式又は（１０）式のＷＫＢ近似形を
用いても問題は生じない。また１、ｌ１ｌｌは次式によ
り与えられる。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) Here, R is the reflection coefficient of the sound wave on the seafloor surface.

深度関数生成手段５は、更に固有関数計算手段６に音源
深度２．及び音源深度２．における波数の２乗値ｋｓ！
を送って音源深度における固有関数Ｕ。＜ｚ＝）を、ま
た受波器深度ｚ、ｌ及び受波器深度２゜における波数の
２乗値に＊２を送って受波器深度における固有関数Ｕ、
（Ｚｉ）を、それぞれ固有関数計算手段６で計算させ、
深度関数２゜（Ｚ５．　Ｚｉ）＝π１ＡｊＪ。（ｚ、）
Ｕ。（２え〉を算出して、その結果を音場累積手段９に
送る。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.

固有関数計算手段６は、深度関数生成手段５から深度２
を、また固有値実部計算手段３からξ。、ｚ、、Ｚ、を
受けると、ｚｕから２までの位相積分値φ、及び２から
２．までの位相積分ｉ４をそれぞれ位相積分計算手段４
で計算させ、φ、とφ４の値を用いて、（２３）式又は
（２６）式により固有関数Ｕ。（ｚ）を計算する。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).

固有関数計算手段６は、固有関数Ｕ、（ｚ）の計算を終
えると、その結果を深度関数生成手段５に送る。When the eigenfunction calculation means 6 finishes calculating the eigenfunctions U and (z), it sends the results to the depth function generation means 5.

一方、固有値虚部計算手段７は、固有関数計算手段３か
ら固有値実部ξ。及び回転点深度２゜。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, 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.

距離関数生成手段８では、固有値実部計算手段３から固
有値実部ξ。を、また固有値虚部計算手段７から固有値
虚部α、を受けると（３）式により複数固有値ｋ。を求
め、所望の距離ｒ、、、　ｍ”１．２．・・・・９Ｍに
対する距離関数ベクトルｐｆｉＴ＝（ρ。１．ρ。２．
・・・・、ρ。Ｍ）（ρ。、＝Ｈ，（１）（Ｈ。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) 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))] The distance function generating means 8 generates a set of distance function values ρ.

ｍ＝１．２．・・・・９Ｍを生成すると、これを音場累
積手段９に送る。m=1.2. ...9M is generated and sent to the sound field accumulation means 9.

音場累積手段９では、深度関数生成手段５から深度関数
（ｚｓ、ｚえ）をまた距離関数生成手段８から距離関数
ベクトルρを送られると両者の積Ｚ。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. The normal mode sound field of mode number n is calculated and added to the sum of the previous normal mode sound fields.

こうして全てのノーマルモード音場についての累加を終
えると、音場累積手段９はその結果１＜”ｚｓ、ｚ次） ｔ）（ｚ、、ｚ、）＝Σｚａ（ｚｓ、ｚ＊）ρ　　・・
・・（３３）れ を出力端子１０から出力する。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.

次に、固有関数計算手段６の実施例として、（２１）式
のｑ（η）及び（２６）式のＡｉ（η）を次 式のように近似した場合の実施例について説明する。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 for η≦η.

ｋ＝０．・・・・、に、−１・・・・（３４）ＡＩ（η
）−αに＋βにη、ηｋ〈ηくηに＋１＋に寓にゆ、・
・・・、ｋ　　　　　　・・・・（３５〉即ち、η＝η
３゜を境界としてη≦η、。の範囲では（２３）式を、
η〉η、。の範囲では（２６）式を用いることにし、そ
れぞれｑ（η）及びＡ、（η）をα、＋β、ηの形で折
れ線近似する。ただしηが負の十分大きな値になるとｑ
（７７）：１となるのでα。−１，β＝Ｏとし、η１に
対応して φ、、、、−２／３（−η１〉８′２　　　　　・・・
・（３６）を定めておく。また、ηが正の十分大きな値
をとる場合には、Ａ、（η）キ０となるので、α、＝β
、＝Ｏとする。この場合−は虚数となる。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.

第３図は固有関数計算手段６の具体例を示すブロック図
である。同図において、各番号は下記部分を示す。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 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.

まず、初期設定手段２から送られてくる海底の層番号ｊ
ｓ及び角周波数ωの平方根ω１″がそれぞれレジスタ２
０及び２８に格納される。また、固有値実部計算手段３
から送られた固有値の２乗ξＩ９回転点深度Ｚ、、Ｚａ
及び回転点深度の層番号ｊ、、Ｌがそれぞれレジスタ１
２，４４．１４に格納される。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.

次に、深度関数生成手段５から深度２及び深度における
波数の２乗ｋ　！（ｚ　）を送られると、それぞれレジ
スタ４５及び１１に格納し、レジスタ４４の２．とレジ
スタ４５の２とを位相積分計算手段４に送って位相積分
値φ。を受け、レジスタ１３に格納する。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.

次に比較器２５は、レジスタ１４の３段レジスタ１９の
内容０とを比較し、ｊ、＝０の場合にはレジスタ１４か
ら０を、またｊｕ≠０の場合にはレジスタ１４からπ／
４を読出して加算器１８に送る。きらに比較器２５は、
ｊ。＝０の場合にはレジスタ２１の内容φｒｍａｘが、
またｊｕ≠０の場合にはレジスタ１３の内容φ。が比較
器３０に送られるように、切換スイッチ２４を制御し、
該切換スイッチ２４の出力を−ｔとする。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.

加算器１８は、レジスタ１３の内容φ、にレジスタ１７
から送られた値を加算し、正弦関数発生器２９におくる
。該正弦関数発生器２９は、加算器１８から送られてき
た位相値の正弦（−啜め、その結果を乗算器３２に送る
。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 .

次に、レジスタ４５の２としレジスタ４４の２６とが位
相積分計算手段４に送ら′れ、該位相積分計算手段４で
位相積分値φ、が計算されてレジスタ１５に格納される
と、比較器２６はレジスタ１４の内容ｊｄとレジスタ２
０の内容ｊ、とを比較し、Ｊａ−Ｊｌｌの場合には、レ
ジスタ２１の内容φ０．工を、またＪ、＋≠ｊｓの場合
には、レジスタ１５の内容φ６を比較器３０に送るよう
に切換スイッチ２７を制御し、切換スイッチ２７の出力
をφ４′とする。比較器３０では、切換スイッチ２４及
び２７からそれぞれ位相データを受は取ると、両者を比
較し、次の判断に基づいて一方を選択する。■φｕ’ｌ
ｄ′が共に実数の場合には値の小さい方を選択する。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.

■φ、°、φ４゛の一方が実数、他方が虚数の場合には
虚数の方を選択する。■If one of φ, °, and φ4′ is a real number and the other is an imaginary number, select the imaginary number.

■φ８′、φ４′がともに虚数の場合には絶対値の大き
い方を選択する。(2) If both φ8' and φ4' are imaginary numbers, select the one with the larger absolute value.

比較器３０はφｕ’＋’ｂｄ′の一方を選択するとその
値を−（３／２（・）　）　””演算器３３に送ると共
に、選んだφ゛に対応した回転点の層番号の値をレジス
タ１４から読出して、その番号に対応した７、の値をレ
ジスタ群２２から読出し、切換スイッチ４２に送る。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(·))"" 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.

近似関数判定手段３５は−（３／２（・）　：］　”’
演算器３３かも送られた刀の値とレジスタ３４に格納さ
れた関数近似の境界値と比較し、η、くη≦η、１　　
　　　　　　　　・・・・（３８）となるｋを選び出し
て、その番号を比較器３９に送るとともにｋに対応した
近似関数の係数β、をレジスタ３８から選び出して乗算
器３７に、そして係数α、をレジスタ４１から選び出し
て加算器４０に送る。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.

乗算器３７は前記−（３／２（・）　）　”’演算器３
３から送られたηにレジスタ３８から送られたβ、を乗
し、その結果βに７７を加算器４０に送る。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.

加算器４０は乗算器３７からβ、ηの値を送られると、
これにレジスタ４１から送られたαを加え、その結果を
乗算器４３に送る。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.

比較器３９は、近似関数判定手段３５から送られた近似
関数の番号にと、レジスタ３６に格納されたｑ（η）、
ＡＩ（７７＞との境界の番号に０とを比較し、ｋ　＜　
ｋ　ｏの場合は加算器１６に起動信号を送ると共に、乗
算器３２の出力が乗算器４３に行くように、切換スイッ
チ４２を制御する。またに≧に０の場合にはレジスタ２
２の出力が乗算器４３に行くように切換スイッチ４２を
制御する。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.

加算器１６は比較器３９から起動信号を受けると、レジ
スタ１１に格納された波数の２乗値に！（２）から、レ
ジスタ１２に格納きれた固有値の２乗ξ。２を減じ、そ
の結果ｋ”（ｚ）−ξ。′を一１／４乗算器２３に送る
。−１／４乗算器２３は加算器１６からｋ”（ｚ）−ξ
、、！の値を受けると、これを−１／４乗し、その結果
を乗算器３１に送る。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.

乗算器３１は一１／４乗算器２３から（ｋｌ（２）−ξ
、２　）　１／４の値を受けると、これにレジスタ２８
に格納されたω１″の値を乗じ、その結果を乗算器３２
に送る。乗算器３２は乗算器３１から送られたωｌ　／
　ｌ　（ｋｌ　（ｚ　）−ξ％　〕−１／４の値に、正
弦関数発生器２９から送られた正弦関数値を乗し、その
結果を切換スイッチ４２に送る。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.

切換スイッチ４２は比較器３９の制御に基づき、乗算器
３２から送られる値、又はレジスタ２２から送られる値
のいずれかを通過させて乗算器４３に送る。乗算器４３
は切換スイッチ４２から送られる値に加算器４０から送
られる近似関数値を乗じ、その結果である固有関数値Ｕ
、（ｚ）を出力端子４６から出力する。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.

〔発明の効果〕〔Effect of the invention〕

以上説明したように本発明に係る音波伝搬シミュレーシ
ョン方式によれば、固有関数を回転点付近でも発散した
り不安定になったりすることがない安定した精度の良い
関数で近似しているので、従来の方式のように回転点近
傍で精度が劣化し、値が発散してしまい計算不能になる
という問題点を除くことができ、特に低周波において計
算精度を大幅に向上きせることができるという優れた効
果が得られる。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.

【図面の簡単な説明】[Brief explanation of drawings]

第１図は本発明に係る音波伝搬シミュレーション方式を
実施する装置の構成を示すブロック図、第２図は関数ｑ
（η）の振舞を示す説明図、第３第１図の固有関数計算
手段の具体的構成例を示す図である。 図中、１・・・・入力端子、２・・・・初期設定手段、
３・・・・固有値実部計算手段、４・・・・位相積分計
算手段、５・・・・深度関数生成手段、６・・・・固有
関数計算手段、７・・・・固有値虚部計算手段、８・・
・・距離関数生成手段、９・・・・音場累積手段、１０
，４６・・・・出力端子、１１　、１２　、１３　、１
５　、１７．１９，２０，２１，２２，２８，３４．３
Ｂ、３８，３９，４４，４５．・・・・レジスタ、１６
．１８．２０・・・・加算器、２３・・・・−１／４乗
算器、２４，２７．４２・・・・切換スイッチ、２５．
２６，３０．３９・・・・比較器、２９・・・・正弦関
数発生器、３１，３２．４３・・・・乗算器、３３・・
・・−（３／２（・））”３演算器、３５・・・・近似
関数判定手段。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.