JPH052193B2

JPH052193B2 -

Info

Publication number: JPH052193B2
Application number: JP23645185A
Authority: JP
Inventors: Hisayasu Ootsubo; Shunji Ozaki
Original assignee: Oki Electric Industry Co Ltd
Current assignee: Oki Electric Industry Co Ltd
Priority date: 1985-10-24
Filing date: 1985-10-24
Publication date: 1993-01-11
Also published as: JPS6296878A

Description

【発明の詳細な説明】（産業上の利用分野）本発明は、水中の音波伝搬特性を解析する音波
伝搬シミユレーシヨン方法に関し、特に音線理論
に基づいたシミユレーシヨン方法に関するもので
ある。（従来の技術）従来のこの種の技術は、「成層媒体内の波動」
（L.M.Brekhovskikh著，「Waves in Latered
Media」第２版，Academic Press（米国）発行，
第192−199頁）、及び（日本音響学会講演論文集，
56〔10〕（昭和56年10月），１−６−15，尾崎俊二
著，「低周波カツトオフ効果の音線理論への一導
入法」，頁535）に記載されている。ソーナー等の海洋音響機器の性能は、その運用
現場の海洋環境により、大きく変化する。したが
つて、機器の効果的運用のために、海洋中の音波
伝搬特性を把握することは重要であり、種々の音
波伝搬モデルが開発されている。音線モデルもその１つである。角周波数ωの単
位強度の無指向性点音源による音場は、媒質特性
が水平方向に変化しない場合、多重経路の累積和
として次式のように表わされる。 Ψ（ｒ，z_R；z_S，ω）＝_∞ 〓ⁿ⁼⁰ Ψ⁽ⁿ⁾〓（ｒ，z_R；z_S，ω） (1) Ψ⁽ⁿ⁾〓（ｒ，z_R；z_S，ω）∫^∞ ₀（iωλ／2πr）^1/2F⁽〓⁾ ₁（z_S；λ，ω）F⁽〓⁾ ₂（z_R；λ，ω）γ⁽ⁿ⁾〓(λ
,ω)exp(iωλr)dλ(2) ここで、座標系として円筒座標系（ｒ，φ，
ｚ）を用い、音源直上の海面を原点として、ｚ軸
を深さ方向にとつている。音源位置は（ｏ，ｏ，
z_S），受波器位置は（ｒ，ｏ，z_R）とした。ｎは
音線が経路の最下点（海底または下側の転回点）
を通過する回数、νは音線のタイプ（音源から上
下いずれの方向に放射されるか、受波点に上下い
ずれの方向から到来するかにより定まる）を表わ
すパラメータであり、ωλは波数の水平成分を表
わす。またγ⁽ⁿ⁾〓（λ）は音線の最上点（海面または
上側の転回点）および最下点での境界条件に関係
する関数であり、F⁽〓⁾ ₁（z_S；λ，ω），F⁽〓⁾ ₂（z_R；
λ，
ω）は音場の深度変化に関係する関数である。以
下ではλを音線パラメータと呼ぶことにする。 F⁽〓⁾ ₁（z_R；λ，ω），F⁽〓⁾ ₂（z_S；λ，ω）はとも
に、
次の方程式を満足する。 ∂²F（ｚ；λ，ω）／∂z² ＋ω²〔c^-2(z)−λ²〕Ｆ（ｚ；λ，ω）＝０ (3) ここでｃ(z)は媒質中の音速分布である。 (3)式を解き、関数Ｆ（ｚ；λ，ω）を求めるの
はかなり複雑であり、音速プロフアイルの形によ
つては解析解が得られない場合がある。そのため
従来は、簡単な近似手法として、次のWKBの１
次近似の形がよく用いられている。 F_±（ｚ；λ，ω）＝〔c^-2(z)−λ²〕^-1/4 exp〔±iωQ（z_u，ｚ；λ）〕 (4) Ｑ（z_u，ｚ；λ）＝∫^z _zu〔C^-2（ζ）−λ²〕^1/2dζ(
5) ここでz_uは基準深度であり、通常は音線の最上
点深度がとられる。位相項の複号は−が上向き、
＋が下向きの波を表わす。 (2)式中のF⁽〓⁾ ₁，F⁽〓⁾ ₂はνの値によりF_-，F₊のい
ずれかの形をとる。たとえば音線のタイプνを第
３図にｎ＝１の例で示すように ν＝１：音源から下向きに出て、受波点に下から
到来する。 ν＝２：音源から下向きに出て、受波点に上から
到来する。 ν＝３：音源から上向きに出て、受波点に下から
到来する。 ν＝４：音源から上向きに出て、受波点に上から
到来する。というように定義すると、F⁽〓⁾ ₁F⁽〓⁾ ₂はν＝１〜４に
ついて、それぞれ次のように表わされる。 F⁽¹⁾ ₁（z_S；λ，ω），F⁽¹⁾ ₂（z_R；λ，ω）＝ｈ（z_S，z_R；λ，ω）exp 〔−iωQ₊（z_S，z_R；λ）〕 F⁽²⁾ ₁（z_S；λ，ω），F⁽²⁾ ₂（z_R；λ，ω）＝ｈ（z_S，z_R；λ，ω）exp 〔−iωQ_-（z_S，z_R；λ）〕 F⁽³⁾ ₁（z_S；λ，ω），F⁽³⁾ ₂（z_R；λ，ω）＝ｈ（z_S，z_R；λ，ω）exp 〔＋iωQ_-（z_S，z_R；λ）〕 F⁽⁴⁾ ₁（z_S；λ，ω），F⁽⁴⁾ ₂（z_R；λ，ω）＝ｈ（z_S，z_R；λ，ω）exp 〔＋iωQ₊（z_S，z_R；λ）〕 (6) ここで Q₊（z_S，z_R；λ）＝Ｑ（z_u，z_S；λ）＋Ｑ（z_u，z_R；λ） Q₊（z_S，z_R；λ）＝Ｑ（z_u，z_S；λ）＋Ｑ（z_u，z_R；λ） Q_-（z_S，z_R；λ）＝Ｑ（z_u，z_S；λ）−Ｑ（z_u，z_R；λ
）(7) である。また、ｈ（z_S，z_R；λ，ω）＝ｆ（z_S；λ，ω）ｆ（z_R；λ，ω） (8) であり、ｆ（ｚ；λ，ω）はF_±（ｚ；λ，ω）の
振幅項を表わす。 (6)式を(2)式に代入すると、 Ψ⁽ⁿ⁾〓（ｒ，z_R；z_S，ω）∫^∞ ₀（iωλ
／2πr）^1/2ｈ（z_S，z_R；λ，ω） γ⁽ⁿ⁾〓（λ，ω）×exp｛iω〔λr±iωQ_±（z_S，z_R；λ）〕｝dλ(9) が得られる。(9)式の積分は一般に解析的に解くこ
とはできず、通常はλをサウンド・チヤネルに対
した区間に分割して区間ごとに停留位相近似や数
値積分などを用いて計算される。その際、λの適
当な間隔のサンプル値に対応した音線の計算が必
要になる。（発明が解決しようとする問題点）しかしF⁽〓⁾ ₁，F⁽〓⁾ ₂を(4)式で近似することは、次
の
２つの問題点を持つている。１ (4)式の右辺第１項はλ→c^-1(z)とすると無限
大になり、計算できなくなる。またλ＝c^-1(z)
の近傍で近似精度が急速に悪くなる。２ F_±（ｚ；λ，ω）の振幅の周波数依存性が(4)
式の近似では失われてしまう。したがつて、周波数カツトオフ効果等も表わす
ことができない。従来の音線モデルでは、この問題を次のように
して逃がれていた。１ λ＝c^-1(z)の近傍の微小範囲のλについては
音線計算を行わない。２周波数カツトオフ効果の補正を個々の音線に
つき行う。しかし、この方法でもλ＝c^-1(z)の近傍での精
度の劣化は解消されず、さらに、λ＝c^-1(z)の近
傍の微小範囲を計算しないことで音場の一部が計
算されないという別の誤差を生むことになる。ま
た、カツトオフ効果の補正値は音速プロフアイル
ｃ(z)とω，λとの関係だけで決まり、音源や受波
点深度に対する依存性は表わせないという問題が
あつた。この発明は、以上述べた転回点近傍で計算不能
が生じたり、近似誤差が拡大したりする問題点、
および周波数カツトオフ効果等音場の周波数依存
性を十分に表わせないという問題点を除去し、簡
単な処理で高い精度が実現できる実用的な音波伝
搬シミユレーシヨン方法を提供することを目的と
する。（問題点を解決するための手段）前記目的を達成するための本発明の特徴は、不
均質な媒質中の音波伝搬現象をシミユレートする
音線理論に基づく音波伝搬シミユレーシヨン方法
において、深度関数計算手段により、音線パラメ
ータのサンプル値に対応した深度関数の振幅値
を、転回点から音源および受波点深度までの位相
積分値の低次関数の積として算出し、深度関数計
算手段によつて得られた周波数依存性を有する深
度関数の振幅値に従つて音場計算を行う音波伝搬
シミユレーシヨン方法にある。（作用）入力データ及び初期設定値に従つて発生する音
線パラメータに従つて、基本経路及び深度関数が
計算され、これらの出力を入力として音場計算手
段により音場のシミユレーシヨンが行なわれる。
深度関数計算手段は、音線パラメータのサンプル
値に対応した深度関数の振幅値を、転回点から音
源及び受波点深度までの位相積分値の低次関数の
積として算出し、この結果の周波数依存性を有す
る深度関数の振幅値に従つて音場計算手段が音場
計算を行なう。深度関数は転回点付近でも発散し
たり不定になることのない安定な精度のよい関数
で近似される。（発明の原理）まず、本発明の原理を説明する。(4)式による近
似は位相については転回点の近傍でも十分な精度
が保たれているので、振幅項について、精度を改
善することを考える。音源または受波点が上側の
転回点の近傍にある場合を例にとり説明する。転回点z_uの近傍で音速が c^-2(z)＝c^-2 ₀〔１−2g（ｚ−z_u）／c₀〕 (10) の形で近似できるものとする。ここでc₀＝λ^-1で
ある。このとき、次の置き換え a³＝−2gω²／c₀ ³，η＝−ａ（ｚ−z_u） (11) を行うと、(3)式は次のAiryの方程式に書き改め
られる。 ∂²F／∂η²−ηF＝０ (12) (4)式に対応する(12)式の解はAiry関数を用いて、
次の形で与えられる。Ｆ± Ｆ± （ｚ；λ，ω）＝（πω／ａ）^1/2〔Bi（η）±iAi（
η）〕（13）位相は(4)式の位相項で近似できるものとし、振
幅項のみを考えると、振幅は次式のように表わす
ことができる。ｆ（ｚ；λ，ω）＝（πω／ａ）^1/2〔B_i ²（
η）＋A² _i（η）〕^1/2（14）ここで、ａ，ηを別の形で表わすことを考え
る。(10)式を(5)式に代入すると、Ｑ（z_u，ｚ；λ）＝∫^z _zu〔−2g（ｚ−z_u）／c₀ ³〕^1/2dz＝２／３（−2g／c₀ ³）^1/2（ｚ−z_u）^3/2（15）
したがつて、(11)，（15）式より、次式を得る。 η＝−［２／３ωQ（z_u，ｚ；λ）］^2/3（16）また、(10)，(11)式からａ＝ω〔c^-2(z)−λ²〕^1/2／（−η）^1/2 （17）が得られる。（17）式を（14）式に代入すると、ｆ（ｚ；λ，ω）＝π^1/2｛−η／〔c^-2(z)
−λ²〕｝^1/4〔B_i ²（η）＋A² _i（η）〕^1/2（18）を得る。（18）式は、ηを介して角周波数ωの関数とな
つており、また〔c^-2(z)−λ²〕^-1/4の項の転回点付
近での増大は、（−η）^1/4の項により相殺されて、
従来の方法の問題点は基本的に解消される。ま
た、ηを（16）式の形で与えたことにより、音速
分布が(10)式からはずれた部分でも用いることがで
きる。しかし、（18）式の振幅項中の｛−η／
〔c^-2(z)−λ²〕｝^1/4の部分はλ＝c^-(z)とすると不定
に
なつてしまうため、別の形で表現することが必要
である。これは(11)式を（14）式に代入すると容易
に得られる。ｆ（ｚ；λ，ω）＝（π／λ）^1/2｜ω／2g
｜^1/6〔B_i ²（η）＋A² _i（η）〕^1/2（19）したがつて、ηの値によつて、（18），（19）式
を使い分けるようにすれば、関数Ｆ± （ｚ；λ）
を精度よく求めることができ、それを(2)式に代入
することができる。ここで（18）式の計算をさらに簡略化すること
を考える。次式の形におく。ｆ（ｚ；λ，ω）＝ｐ（η）〔c^-2(z)−λ²〕｝^-1/4（2
0）そうするとｐ（η）は次のように表わされる。ｐ（η）＝π^1/2（−η）^1/4〔B_i ²（η）＋A_i ²（η）
〕^1/2
（21）ｐをηでなく次の位相積分値 ζ＝ωQ（z_u，ｚ；λ）＝２／３（−η）^3/2（22）の関数としてプロツトすると第４図のようにな
る。第４図からわかるようにｐ（ζ）はζの非常
に滑らかな関数になつており、ζの値が大きくな
ると１に漸近していく。すなわち、（18）式は(4)
式の振幅項に漸近する。ｐをζでなくηの関数と
考えても、同様のことが言える。したがつて、ｐ
（ζ）またはｐ（η）を簡単な関数形で近似してお
けば、（20）式により高精度のｆ（ｚ；λ，ω）を
容易に計算することができる。一方、（19）式は音速プロフアイルが次式 e^-2(z)＝c^-2 _j〔１−2g_j（ｚ−z_j）／c_j，z_j
ｚz_j+1，Ｊ＝１，……，Ｎ（23）の形で表わされている場合には、さらに簡単にな
る。転回点がｊ層内（z_jz_uz_j+1）にある場合
を考える。（23）式を微分し整理すると、 ∂c(z)／∂z＝g_jc³(z)／c_j ³ （24）（24）式にｚ＝z_uを代入すると、次式を得る。ｇ＝∂c（z_u）／∂z＝g_j（c_jλ）^-3 （25）（25）式を（19）式に代入すると、ｚがｊ層内
にある場合には、ｆ（ｚ；λ，ω）＝（πc_j）^1/2｜ω／2g_j｜^1/6〔B
_i ²（η）＋A_i ²（η）〕^1/2（26）となることがわかる。（26）式中の〔B_i ²（η）＋A_i ²（η）〕^1/2の項はη
またはζの極めて滑らかな関数となり、１次式等
の低次式で十分近似できる。したがつて、（26）
式は各層ごとにζ（またはη）の低次式として近
似できることになる。したがつて、深度ｚが転回点深度z_uと同じ層内
にある場合には、振幅項として（26）式を用い、
異なつた層内にある場合には（18）式を用いるこ
とにすれば、簡単な計算で精度のよい深度関数値
を得ることができる。音源または受波点深度ｚが下側の転回点深度z_d
の近くにある場合にも、同様の方法で精度のよい
深度関数値が得られる。この場合、ζおよびηは
次のようになる。 ζ＝ωQ（ｚ，z_d；λ）＝ω〔Ｑ（z_u，z_d；λ）−Ｑ
（z_u，ｚ；λ）〕（27） η＝−（３／２ζ）^2/3＝−｛３／２ω〔Ｑ（z_u，z_d
；λ）−Ｑ（z_u，ｚ；λ）〕｝^2/3（28）（実施例）次に本発明の実施例について説明する。第１図
は、本発明の実施例の説明図である。図において
１は環境条件、角周波数等のデータの入力端子、
２は前記入力端子１からの入力データに基づき音
速プロフアイル、音波の減衰係数、海面・海底損
失等を設定する初期設定手段、３は前記初期設定
手段２から送られる音速プロフアイルおよび音
源／受波器深度情報に基づき、音場計算に用いる
音線パラメータλのサンプル値を発生させる音線
パラメータ発生手段、４は前記初期設定手段２か
ら送られる音速プロフアイル、音源／受波器深度
情報、および前記音線パラメータ発生手段３から
送られる音線パラメータλに基づき、音線の基本
経路に関する量を算出する基本経路計算手段、５
は前記初期設定手段２から送られる音源および受
波点深度情報および角周波数情報、および前記基
本経路計算手段４から送られるλ²，転回点の層番
号、およびＱ（z_u，z_S；λ），Ｑ（z_u，z_R；λ），Ｑ
（z_u，z_d；λ）の値から、深度関数の振幅を算出
する深度関数計算手段、６は前記初期設定手段２
から送られる角周波数ω、音波の減衰係数、海
面・海底損失、前記音線パラメータ発生手段３か
ら送られる音線パラメータλ、前記基本経路計算
手段４から送られる基本経路に関する種々の量、
および前記深度関数計算手段５から送られる深度
関数の振幅を用い、音場を計算する音場計算手
段、７は前記音場計算手段６で得られた音場計算
結果を出力する出力端子である。この装置の動作原理を説明する。初期設定手段
２は、入力端子１から環境条件等の計算条件が入
力されると、音速プロフアイルに関数近似を施
し、（26）式の右辺の a_j＝（πc_j）^1/2｜ω／2g_j｜^1/6，ｊ＝１，２，……，
Ｎ（29）の計算を行う。また、角周波数ωの音波の海面・
海面損失特性、減衰係数を算出する。さらに、音
源、受波点深度の音速c_S，c_Rおよび層番号j_S，j_R
を決定する。初期設定手段２は、これらの計算を
終えると、音速プロフアイルと音源／受波点深度
情報を音線パラメータ発生手段３、および基本経
路計算手段４に、音源、受波点の音速の逆２乗c_S
^-2，c_R ^-2および層番号j_S，j_Rを基本経路計算手段４
および深度関数計算手段５に、角周波数ωを深度
関数計算手段５および音場計算手段６に、そし
て、音波の減衰係数、海面・海底損失等を音場計
算手段６に送る。音線パラメータ発生手段３では、初期設定手段
２から送られた音速プロフアイル、および音源、
受波器深度に基づき、適当な個数の音線パラメー
タλのサンプル値λ_i，ｉ＝１，……，Ｍを決定
し、基本経路計算手段４および音場計算手段６に
送る。基本経路計算手段４では、初期設定手段２から
音速プロフアイルと音源、受波点深度情報を、そ
して音線パラメータ発生手段３から音線パラメー
タλ_i，ｉ＝１，……，Ｍを送られると、ｉ＝１か
らＭまで順に、次の処理を行う。１音線の最上端z_u、および最下端z_dが含まれる
層番号j_u，j_dを求める。２音線の３つの基本経路に関する種々の量の計
算を行う。３つの基本経路の例としては、ａ音線の最上点z_uから音源深度z_Sまで、ｂ音線の最上点z_uから受波点深度z_Rまで、ｃ音線の最上点z_uから最下点z_dまでの３つがある。この３つの基本経路のそれぞれ
について、次の３つの量、ａ水平距離ｒ（λ_i）ｂ水平距離ｒのλによる微分値 ∂r／∂λ〓_i ｃ (5)式のＱ（z_u，ｚ；λ_i）を求めておけば、音源深度と受波深度とを結
び、音源パラメータλ_iの任意の音線に関する量
は、これらの量の線形結合として表わすことが
できる。３基本経路に関する計算結果のすべての量を音
場計算手段６に送るとともに、j_u，j_d，λ² _iおよ
びＱについての計算結果を深度関数計算手段５
に転送する。深度関数計算手段５では初期設定手段２から送
られた音源受波点深度情報、角周波数情報、およ
び基本経路計算手段４から送られたj_u，j_d，λ² _i，
およびＱ（z_u，z_S；λ_i），Ｑ（z_u，z_R；λ_i），Ｑ（z_u，
z_d；λ_i）を用いて、音源および受波点深度におけ
るζの値を求め、（18）または（26）式および(8)
式を用いてｈ（z_S，z_R；λ_i，ω）を計算する。計算
を終えると、深度関数計算手段５はｈ（z_S，z_R；
λ_i，ω）の値を音場計算手段６に送る。音場計算手段６は、初期設定手段２から送られ
た角周波数、海面・海底反射特性、および減衰係
数、音線パラメータ発生手段３から送られた音線
パラメータのサンプル値λ_i，ｉ＝１，……，Ｍ，
基本経路計算手段４から送られた基本経路に関す
る前記の種々の量、および深度関数計算手段５か
ら送られた深度関数ｈ（z_S，z_R；λ_i，ω），ｉ＝１，
……，Ｍを用いて、(1)，(2)式により音場の計算を
行う。計算結果は、出力端子７から出力される。次に、前記深度関数計算手段５の実施例として
（21）式および（26）式をそれぞれ次式のように
近似した場合の実施例について説明する。ｐ（ζ）＝α_k＋β_kζ，ζ_kζζ_k+1，ｋ＝１，…
…，Ｋ（30）ｆ（ｚ；λ，ω）＝a_j（α_p＋β_pζ），z_jｚz_j+1，ｊ＝１，……，Ｎ（31）ただし、ζ_k+1＝∞とする。第２図は深度関数計算手段５の実施例の説明図
である。図において、８は前記初期設定手段２か
らの入力端子、９は前記基本経路計算手段４から
の入力端子、１０，１１，１２，１３はそれぞれ
角周波数ω、音源／受波点深度の層番号、音源／
受波点深度における音速の逆数の２乗、および
（29）式のa_j，ｊ＝１，……，Ｎを格納するレジ
スタ、１４，１５，１６，１７はそれぞれ、λ² _i
およびQ_S＝Ｑ（z_u，z_S；λ）およびQ_R＝Ｑ（z_u，
z_R；λ），Q_C＝Ｑ（z_u，z_d；λ）、およびj_u，j_dの値
を格納するレジスタ、１８，２１，２３，２９は
加算器、１９は比較器、２０，２８，３１，３３
は乗算器、２２は入力された値の−1/4乗を求め
るベキ乗器、２４はζの値等によつて振幅の計算
に用いる近似式を選択する近似関数判定器、２
５，２６，２７は、それぞれζ_k，ｋ＝１，……，
Ｋ、β_k，ｋ＝０，……，Ｋ、α_k，ｋ＝１，……，
Ｋの値を格納しておくレジスタ、３０はa_jおよび
〔c^-2(z)−λ²〕^1/4の値を格納しておくレジスタ、３
２は振幅関数値ｆ（z_S；λ，ω）およびｆ（z_R；
λ，ω）を格納するレジスタ、そして３４は結果
の出力端子である。この装置の動作原理を説明する。まず、レジス
タ群２５には、（30）式によるｐ（ζ）の曲線補間
の境界値を、またレジスタ群２６，２７には、そ
れぞれ（30）式の補間関数値の係数を格納してお
く。入力端子８から入力されたω，j_Sおよびj_R，
c_S ^-2およびc_R ^-2，a_jS，a_jRの情報はそれぞれレジス
タ（群）１０，１１，１２，１３に格納される。
入力端子９空、λ_i，Q_SおよびQ_R，Q_C，j_uおよびj_d
の情報を受けると、それぞれの値が、レジスタ
（群）１４，１５，１６，１７に格納される。レ
ジスタ１５内のQ_Sとレジスタ１６内のQ_Cとは加
算器１８に送られ、Q_C−Q_Sが作り出されて、比
較器１９に送られる。比較器１９は、加算器１８
から値が送られてくると、レジスタ１５の内容
Q_Sと比較し、小さい方の値を乗算器２０に送り
出すとともに、Q_SQ_C−Q_Sの場合はレジスタ１
７からj_uを、またQ_S＞Q_C−Q_Sの場合はレジスタ１
７からj_dを加算器２３に送る。一方、レジスタ１
１の内容j_Sにより、レジスタ群１３から、j_S番目
の内容a_jSが選び出されてレジスタ３０に格納さ
れる。また、加算器２１はレジスタ１２および１
４から、それぞれc_S ^-2およびλ_i ²を読み出してc_S ^-2
−λ_i ²を作り出し、ベキ乗器２２に送る。ベキ乗
器２２は加算器２１から送られた情報を−1/4乗
し、レジスタ３０に格納する。乗算器２０は比較器１９から送られたＱの最小
値にレジスタ１０の内容ωを乗じてζを作り、近
似関数判定器２４、および乗算器２８に送る。加
算器２３はレジスタ１１および１７から送られた
内容の差をとり、その結果を近似関数判定器２４
に送る。近似関数判定器２４は次の基準に基づい
てレジスタ群２６，２７およびレジスタ３０の内
容を選び出し、それぞれ乗算器２８、加算器２９
および乗算器３１に送る。１ ζζ_kの場合、レジスタ群２６の中から、β_k
を、またレジスタ群２７からα_kを選び出し、ま
た、レジスタ３０から、｛c_S ^-2−λ_i ²｝^-1/4を選び
出してそれぞれ乗算器２８、加算器２９および
乗算器３１に送る。２ ζ＜ζ_kかつ加算器２３から送られてくる内容
が０の場合には、レジスタ群２６からβ_pを、レ
ジスタ群２７からα_pを、そして、レジスタ３０
からa_jSを選び出して、それぞれ乗算器２８、
加算器２９、および乗算器３１に送る。３ ζ＜ζ_kでかつ加算器２３から送られてくる内
容が０以外の場合には、ζ_kζζ_k+1を満足す
るβ_kおよびα_kをそれぞれレジスタ群２６および
２７から選び、またレジスタ３０から｛c_S ^-2−
λ_i ²｝^-1/4を選び出して、それぞれ乗算器２８、
加算器２９、および乗算器３１に送る。乗算器２８は乗算器２０から送られた内容にレ
ジスタ群２６から送られた内容を乗じてβ_kζを作
り、加算器２９に送る。加算器２９では、乗算器
２８から送られた内容β_kζに、レジスタ２７から
送られたα_kを加え、（α_k＋β_kζ）を作つて、乗算
器３１に送る。乗算器３１は、加算器２９から送
られた内容に、レジスタ３０から送られた内容
〔a_j又は｛c_S ^-2−λ_i ²｝^-1/4〕を乗じ、f_S＝ｆ（z_S；λ
，
ω）を作り出してレジスタ３２に送る。レジスタ３２にf_Sが送られると、音源深度に関
する計算が終了し、次は受波深度についての計算
に移る。レジスタ１１，１２、および１５は内容
がそれぞれ１つずつシフトされ、音源に関する量
と受波点に関する量とが入れ替えられる。以下音
源について行つたと同様の処理手順によりf_R＝ｆ
（z_R；λ，ω）が計算されて、結果がレジスタ３
２に送られてくる。レジスタ３２は内容が１つシ
フトされ、f_Sとf_Rとが並ぶことになる。乗算器３
３はレジスタ３２にf_Sとf_Rの両方がそろうと、両
者の積をとり、その結果ｈ（z_S，z_R；λ，ω）を
出力端子３４から出力して音場計算手段６に送
る。（発明の効果）以上説明したように、本発明では深度関数を、
転回点付近でも発散したり不定になることのない
安定した精度のよい関数で近似しているので、従
来の方式の持つ、転回点近傍で精度が悪くなる、
周波数カツトオフ効果を表わすことができない等
の問題点を除くことができ適用周波数範囲を大幅
に拡大することができる。また、深度関数の近似を非常に簡単な処理によ
り実現しているので、従来の方式に比べて処理時
間の増加はほとんどなく、ソーナーの運用現場で
の性能評価等、実時間性を要求されるシステムに
も利用することができる。 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 simulation method based on acoustic ray theory. (Conventional technology) This type of conventional technology uses "wave motion in a stratified medium"
(Waves in Latered by LMBrekhovskikh)
Media” 2nd edition, published by Academic Press (USA),
pp. 192-199) and (Acoustical Society of Japan Proceedings,
56 [10] (October 1982), 1-6-15, Shunji Ozaki, "A method for introducing the low-frequency cut-off effect into sound ray theory," p. 535). 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 sound ray model is one of them. If the medium characteristics do not change in the horizontal direction, the sound field caused by a non-directional point source of unit intensity with an angular frequency ω can be expressed as the cumulative sum of multiple paths as shown in the following equation. Ψ (r, z _R ; z _S , ω) = _∞ 〓 ⁿ⁼⁰ Ψ ⁽ⁿ⁾ 〓 (r, z _R ; z _S , ω) (1) Ψ ⁽ⁿ⁾ 〓 (r, z _R ; z _S ,ω)∫ ^∞ ₀ (iωλ／2πr) ^1/2 F ⁽ 〓 ⁾ ₁ (z _S ; λ, ω) F ⁽ 〓 ⁾ ₂ (z _R ; λ, ω) γ ⁽ⁿ⁾ 〓(λ
,ω)exp(iωλr)dλ(2) Here, the coordinate system is a cylindrical coordinate system (r, φ,
z), the origin is the sea surface directly above the sound source, and the z-axis is taken in the depth direction. The sound source position is (o, o,
z _S ), and the receiver position was (r, o, z _R ). n is the lowest point on the path of the sound ray (turning point on the ocean floor or lower side)
, ν is a parameter representing the type of sound ray (determined by whether it is radiated upward or downward from the sound source, or whether it arrives at the receiving point from above or below), and ωλ is the horizontal wave number. Represents a component. In addition, γ ⁽ⁿ⁾ 〓(λ) is a function related to the boundary conditions at the top point (sea level or upper turning point) and bottom point of the sound ray, and F ⁽ 〓 ⁾ ₁ (z _S ; λ, ω ), F ⁽ 〓 ⁾ ₂ (z _R ;
λ,
ω) is a function related to changes in the depth of the sound field. In the following, λ will be referred to as a sound ray parameter. Both F ⁽ 〓 ⁾ ₁ (z _R ; λ, ω) and F ⁽ 〓 ⁾ ₂ (z _S ; λ, ω) are
The following equation is satisfied. ∂ ² F (z; λ, ω) / ∂z ² + ω ² [c ^-2 (z)−λ ² ] F (z; λ, ω) = 0 (3) Here, c(z) is the This is the sound velocity distribution. Solving equation (3) and finding the function F(z;λ,ω) is quite complicated, and depending on the shape of the sound speed profile, an analytical solution may not be obtained. Therefore, conventionally, as a simple approximation method, the following WKB 1
The form of second approximation is often used. F _± (z; λ, ω) = [c ^-2 (z)−λ ² ] ^-1/4 exp [±iωQ (z _u , z; λ)] (4) Q (z _u , z; λ) =∫ ^z _zu [C ^-2 (ζ)−λ ² ] ^1/2 dζ(
5) Here, z _u is the reference depth, and usually the depth of the highest point of the sound ray is taken. The sign of the phase term is - with upward direction,
+ represents a downward wave. F ⁽ 〓 ⁾ ₁ and F ⁽ 〓 ⁾ ₂ in equation (2) take the form of either F _- or F ₊ depending on the value of ν. For example, the sound ray type ν is shown in FIG. 3 as an example of n=1: ν=1: the sound ray exits downward from the sound source and arrives at the receiving point from below. ν=2: The sound comes out downward from the sound source and arrives at the receiving point from above. ν=3: The sound comes out upward from the sound source and arrives at the receiving point from below. ν=4: The sound comes out upward from the sound source and arrives at the receiving point from above. When defined as follows, F ⁽ 〓 ⁾ ₁ F ⁽ 〓 ⁾ ₂ can be expressed as follows for ν = 1 to 4, respectively. F ⁽¹⁾ ₁ (z _S ; λ, ω), F ⁽¹⁾ ₂ (z _R ; λ, ω) = h(z _S , z _R ; λ, ω) exp [−iωQ ₊ (z _S , z _R ; λ)] F ⁽²⁾ ₁ (z _S ; λ, ω), F ⁽²⁾ ₂ (z _R ; λ, ω) = h(z _S , z _R ; λ, ω) exp [−iωQ _- (z _S , z _R ; λ)] F ⁽³⁾ ₁ (z _S ; λ, ω), F ⁽³⁾ ₂ (z _R ; λ, ω) = h(z _S , z _R ; λ, ω) exp [+iωQ _- (z _S , z _R ; λ)] F ⁽⁴⁾ ₁ (z _S ; λ, ω), F ⁽⁴⁾ ₂ (z _R ; λ, ω) = h(z _S , z _R ; λ, ω) exp [+iωQ ₊ (z _S , z _R ; λ)] (6) Here, Q ₊ (z _S , z _R ; λ) = Q(z _u , z _S ; λ) + Q(z _u , _zR ;λ) Q ₊ (z _S , z _R ; λ) = Q (z _u , z _S ; λ) + Q (z _u , z _R ;λ) Q _- (z _S , z _R ; λ) = Q (z _u , z _S ; λ) − Q (z _u , z _R ; λ
)(7). Also, h (z _S , z _R ; λ, ω) = f (z _S ; λ, ω) f (z _R ; λ, ω) (8), and f (z; λ, ω) is F _± represents the amplitude term of (z; λ, ω). Substituting equation (6) into equation (2), Ψ ⁽ⁿ⁾ 〓(r, z _R ; z _S , ω)∫ ^∞ ₀ (iωλ
/2πr) ^1/2 h(z _S , z _R ; λ, ω) γ ⁽ⁿ⁾ 〓(λ, ω)×exp{iω[λr±iωQ _± (z _S , z _R ;λ)]}dλ(9) is obtained. The integral in equation (9) cannot generally be solved analytically, and is usually calculated by dividing λ into intervals for the sound channel and using stationary phase approximation, numerical integration, etc. for each interval. In this case, it is necessary to calculate sound rays corresponding to sample values at appropriate intervals of λ. (Problems to be Solved by the Invention) However, approximating F ⁽ 〓 ⁾ ₁ and F ⁽ 〓 ⁾ ₂ by equation (4) has the following two problems. 1 The first term on the right side of equation (4) becomes infinite when λ→c ^-1 (z) and cannot be calculated. Also, λ=c ^-1 (z)
The approximation accuracy deteriorates rapidly near . The frequency dependence of the amplitude of 2 F _± (z; λ, ω) is (4)
It is lost in the approximation of the equation. Therefore, frequency cutoff effects and the like cannot be expressed. In the conventional sound ray model, this problem was avoided in the following way. 1. Sound ray calculations are not performed for λ in a very small range near λ=c ^-1 (z). 2. Correct the frequency cutoff effect for each sound ray. However, even with this method, the deterioration of accuracy in the vicinity of λ = c ^-1 (z) is not resolved, and furthermore, because the minute range in the vicinity of λ = c ^-1 (z) is not calculated, part of the sound field is This will result in another error in that it will not be calculated. Furthermore, the correction value for the cut-off effect is determined only by the relationship between the sound velocity profile c(z) and ω and λ, and there is a problem in that the dependence on the sound source and the depth of the receiving point cannot be expressed. This invention solves the above-mentioned problem of inability to calculate near the turning point or increase of approximation error,
The present invention aims to provide a practical sound wave propagation simulation method that can achieve high accuracy with simple processing and eliminates problems such as frequency cut-off effects and other problems in which the frequency dependence of a sound field cannot be sufficiently represented. (Means for Solving the Problems) A feature of the present invention for achieving the above object is that, in a sound wave propagation simulation method based on sound ray theory that simulates a sound wave propagation phenomenon in a heterogeneous medium, a depth function calculation means is used. The amplitude value of the depth function corresponding to the sample value of the sound ray parameter is calculated as the product of lower-order functions of the phase integral values from the turning point to the depth of the sound source and the receiving point, and is obtained by the depth function calculation means. The present invention provides a sound wave propagation simulation method for calculating a sound field according to the amplitude value of a depth function having frequency dependence. (Operation) A basic path and a depth function are calculated according to sound ray parameters generated according to input data and initial setting values, and a sound field simulation is performed by the sound field calculation means using these outputs as input.
The depth function calculation means calculates the amplitude value of the depth function corresponding to the sample value of the sound ray parameter as a product of lower-order functions of the phase integral values from the turning point to the depth of the sound source and the receiving point, and calculates the frequency of this result. The sound field calculation means performs sound field calculation according to the amplitude value of the dependent depth function. The depth function is approximated by a stable and accurate function that does not diverge or become unstable even near the turning point. (Principle of the invention) First, the principle of the invention will be explained. Since the approximation by Equation (4) maintains sufficient accuracy for the phase even near the turning point, we will consider improving the accuracy for the amplitude term. An example will be explained in which the sound source or wave receiving point is located near the upper turning point. It is assumed that the sound speed near the turning point z _u can be approximated in the form c ^-2 (z) = c ^-2 ₀ [1-2g (z - z _u )/c ₀ ] (10). Here, c ₀ =λ ⁻¹ . At this time, by performing the following substitution a ³ =-2gω ² /c ₀ ³ , η=-a(z-z _u ) (11), equation (3) can be rewritten as the following Airy equation. ∂ ² F／∂η ² −ηF=0 (12) The solution to equation (12) corresponding to equation (4) is given by using the Airy function,
It is given in the form: F± F± (z;λ,ω)=(πω/a) ^1/2 [Bi(η)±iAi(
η)] (13) Assuming that the phase can be approximated by the phase term in equation (4), and considering only the amplitude term, the amplitude can be expressed as in the following equation. f(z;λ,ω)=(πω/a) ^1/2 [B _i ² (
η)+A ² _i (η)] ^1/2 (14) Now, consider expressing a and η in another form. Substituting equation (10) into equation (5), Q(z _u , z; λ)=∫ ^z _zu [−2g(z−z _u )/c ₀ ³ ] ^1/2 dz=2/3 (-2g/c ₀ ³ ) ^1/2 (z−z _u ) ^3/2 (15)
Therefore, from equations (11) and (15), we obtain the following equation. η=−[2/3ωQ(z _u , z; λ)] ^2/3 (16) Also, from equations (10) and (11), a=ω[c ^-2 (z)−λ ² ] ^1/2 /(−η) ^1/2 (17) is obtained. Substituting equation (17) into equation (14), f(z;λ,ω)=π ^1/2 {−η/[c ^-2 (z)
−λ ² ]} ^1/4 [B _i ² (η) + A ² _i (η)] ^1/2 (18) is obtained. Equation (18) is a function of the angular frequency ω via η, and the increase of the term [c ^-2 (z)−λ ² ] ^-1/4 near the turning point is (−η ) offset by the ^1/4 term,
The problems of conventional methods are basically eliminated. Furthermore, by giving η in the form of equation (16), it can be used even in areas where the sound velocity distribution deviates from equation (10). However, {−η/ in the amplitude term of equation (18)
[c ^-2 (z)−λ ² ]} ^{The 1/4} part becomes indeterminate if λ=c ^- (z), so it is necessary to express it in another form. This can be easily obtained by substituting equation (11) into equation (14). f (z; λ, ω) = (π/λ) ^1/2 | ω/2g
| ^1/6 [B _i ² (η) + A ² _i (η)] ^1/2 (19) Therefore, if you use equations (18) and (19) differently depending on the value of η, Function F± (z;λ)
can be calculated with high accuracy and substituted into equation (2). Let us now consider further simplifying the calculation of equation (18). Let it be in the form of the following equation. f(z;λ,ω)=p(η) [c ^-2 (z)−λ ² ]} ^-1/4 (2
0) Then p(η) can be expressed as follows. p(η)=π ^1/2 (−η) ^1/4 [B _i ² (η) + A _i ² (η)
] ^1/2
(21) If p is plotted not as a function of η but as a function of the next phase integral value ζ = ωQ (z _u , z; λ) = 2/3 (-η) ^3/2 (22), the result is as shown in Figure 4. . As can be seen from FIG. 4, p(ζ) is a very smooth function of ζ, and approaches 1 as the value of ζ increases. In other words, equation (18) becomes (4)
Asymptotes to the amplitude term in Eq. The same thing can be said if p is considered as a function of η instead of ζ. Therefore, p
If (ζ) or p(η) is approximated by a simple functional form, it is possible to easily calculate f(z;λ, ω) with high precision using equation (20). On the other hand, in equation (19), the sound velocity profile is expressed as follows: e ^-2 (z)=c ^-2 _j [1-2g _j (z-z _j )/c _j , z _j
It becomes even simpler if it is expressed in the form zz _j+1 , J=1,...,N(23). Consider the case where the turning point is within the j layer (z _j z _u z _j+1 ). Differentiating and rearranging equation (23), ∂c(z)/∂z=g _j c ³ (z)/c _j ³ (24) Substituting z=z _u into equation (24), we get the following equation. . g=∂c(z _u )/∂z=g _j (c _j λ) ^-3 (25) Substituting equation (25) into equation (19), if z is in layer j, f( z; λ, ω) = (πc _j ) ^1/2 | ω／2g _j | ^1/6 [B
It can be seen that _i ² (η) + A _i ² (η)] ^1/2 (26). The term [B _i ² (η) + A _i ² (η)] ^1/2 in equation (26) is η
Alternatively, it becomes an extremely smooth function of ζ, and can be sufficiently approximated by a low-order equation such as a linear equation. Therefore, (26)
The equation can be approximated as a lower-order equation of ζ (or η) for each layer. Therefore, if the depth z is in the same layer as the turning point depth z _u , use equation (26) as the amplitude term,
If they are in different layers, by using equation (18), a highly accurate depth function value can be obtained with simple calculations. Turning point depth z _d where the sound source or receiving point depth z is lower
Even when the depth function is close to , a highly accurate depth function value can be obtained using the same method. In this case, ζ and η are as follows. ζ=ωQ(z, _zd ;λ)=ω[Q( _zu , _zd ;λ)−Q
(z _u , z; λ)] (27) η=-(3/2ζ) ^2/3 =-{3/2ω[Q(z _u , z _d
;λ)−Q(z _u , z;λ)]} ^2/3 (28) (Example) Next, an example of the present invention will be described. FIG. 1 is an explanatory diagram of an embodiment of the present invention. In the figure, 1 is an input terminal for data such as environmental conditions and angular frequency;
2 is an initial setting means for setting the sound speed profile, sound wave attenuation coefficient, sea surface/bottom loss, etc. based on input data from the input terminal 1; 3 is a sound speed profile and sound source/receiver sent from the initial setting means 2; a sound ray parameter generating means for generating a sample value of a sound ray parameter λ used for sound field calculation based on the wave depth information; 4 is a sound velocity profile sent from the initial setting means 2; sound source/receiver depth information; and basic path calculation means 5 for calculating a quantity related to the basic path of the sound ray based on the sound ray parameter λ sent from the sound ray parameter generation means 3.
are the sound source and receiving point depth information and angular frequency information sent from the initial setting means 2, λ ² sent from the basic path calculation means 4, the layer number of the turning point, and Q(z _u , z _S ; λ ), Q (z _u , z _R ; λ), Q
Depth function calculation means for calculating the amplitude of the depth function from the values of (z _u , z _d ; λ); 6 is the initial setting means 2;
angular frequency ω sent from the sound ray attenuation coefficient, sea surface/bottom loss, sound ray parameter λ sent from the sound ray parameter generation means 3, various quantities related to the basic path sent from the basic path calculation means 4,
and a sound field calculation means for calculating a sound field using the amplitude of the depth function sent from the depth function calculation means 5, and 7 is an output terminal for outputting the sound field calculation result obtained by the sound field calculation means 6. . The operating principle of this device will be explained. When the calculation conditions such as environmental conditions are input from the input terminal 1, the initial setting means 2 applies a function approximation to the sound velocity profile, and calculates a _j = (πc _j ) ^1/2 | ω on the right side of equation (26). /2g _j | ^1/6 ,j=1,2,...,
Calculate N (29). Also, the sea surface of the sound wave with angular frequency ω
Calculate sea level loss characteristics and attenuation coefficient. Furthermore, the sound velocities c _S , c _R and the layer numbers j _S , j _R
Determine. After completing these calculations, the initial setting means 2 sends the sound velocity profile and the sound source/receiving point depth information to the sound ray parameter generating means 3 and the basic path calculating means 4, which are the inverse of the sound speed at the sound source and the receiving point. squared c _S
^-2 , c _R ^-2 and layer numbers j _S , j _R as basic route calculation means 4
The angular frequency ω is sent to the depth function calculation means 5 and the sound field calculation means 6, and the attenuation coefficient of the sound wave, sea surface/bottom loss, etc. are sent to the sound field calculation means 6. The sound ray parameter generation means 3 receives the sound velocity profile sent from the initial setting means 2, the sound source,
Based on the receiver depth, an appropriate number of sample values λ _i , i=1, . The basic path calculation means 4 receives the sound velocity profile, sound source, and receiving point depth information from the initial setting means 2, and the sound ray parameters λ _i , i=1, ..., M from the sound ray parameter generation means 3. Then, the following processing is performed in order from i=1 to M. 1 Find the layer numbers j _u and j _d that include the top end _{z u} and the bottom end z _d of the sound ray. 2. Perform calculations of various quantities regarding the three basic paths of sound rays. Examples of the three basic paths are: a From the highest point of the sound ray z _u to the sound source depth z _S , b From the highest point of the sound ray z _u to the receiver depth z _R , c From the highest point of the sound ray z _u There are three points up to the lowest point z _{and d} . For each of these three basic paths, the following three quantities are calculated: a Horizontal distance r(λ _i ) b Differential value of horizontal distance r with respect to λ ∂r/∂λ〓 _i c Q(z _u , z; λ _i ), the sound source depth and reception depth can be connected, and the quantity of the sound source parameter λ _i for any sound ray can be expressed as a linear combination of these quantities. 3. Send all the calculation results regarding the basic path to the sound field calculation means 6, and send the calculation results regarding j _u , j _d , λ ² _i and Q to the depth function calculation means 5.
Transfer to. The depth function calculating means 5 uses the sound source/receiving point depth information and angular frequency information sent from the initial setting means 2, and j _u , j _d , λ ² _i , sent from the basic path calculating means 4 .
and Q(z _u , z _S ; λ _i ), Q(z _u , z _R ; λ _i ), Q(z _u ，
z _d ; λ _i ), find the value of ζ at the depth of the sound source and receiving point, and use equations (18) or (26) and (8)
Calculate h(z _S , z _R ; λ _i , ω) using the formula. After completing the calculation, the depth function calculation means 5 calculates h(z _S , z _R ;
The values of λ _i , ω) are sent to the sound field calculation means 6. The sound field calculation means 6 calculates the angular frequency, sea surface/bottom reflection characteristics, and attenuation coefficient sent from the initial setting means 2, and the sample value λ _i of the sound ray parameter sent from the sound ray parameter generation means 3, i=1. ,...,M,
The aforementioned various quantities related to the basic route sent from the basic path calculation means 4 and the depth function h (z _S , z _R ; λ _i , ω), i=1, sent from the depth function calculation means 5
..., M is used to calculate the sound field using equations (1) and (2). The calculation result is output from the output terminal 7. Next, as an example of the depth function calculating means 5, an example in which equations (21) and (26) are approximated as shown in the following equations will be described. p(ζ)=α _k +β _k ζ, ζ _k ζζ _k+1 , k=1,...
..., K (30) f (z; λ, ω) = a _j (α _p + β _p ζ), z _j zz _j+1 , j=1, ..., N (31) However, ζ _k+1 = ∞. FIG. 2 is an explanatory diagram of an embodiment of the depth function calculation means 5. In the figure, 8 is an input terminal from the initial setting means 2, 9 is an input terminal from the basic path calculating means 4, and 10, 11, 12, and 13 are the angular frequency ω and the layer number of the sound source/receiver point depth, respectively. ,sound source/
Registers 14, 15, 16, and 17 that store the square of the reciprocal of the sound speed at the depth of the receiving point and a _j , j=1, ..., N in equation (29) are λ ² _i
and Q _S =Q(z _u , z _S ;λ) and Q _R =Q(z _u ,
z _R ; λ), Q _C =Q (z _u , z _d ; λ), and registers that store the values of j _u and j _d ; 18, 21, 23, and 29 are adders; 19 is a comparator; 20 ,28,31,33
2 is a multiplier; 22 is a power multiplier that calculates the −1/4th power of the input value; 24 is an approximate function determiner that selects an approximate formula to be used for amplitude calculation based on the value of ζ, etc.;
5, 26, 27 are respectively ζ _k , k=1,...,
K, β _k , k=0, ..., K, α _k , k=1, ...,
A register 30 stores the value of K, and a register 3 stores the value of a _j and [c ^-2 (z)−λ ² ] ^1/4 .
2 is the amplitude function value f(z _S ; λ, ω) and f(z _R ;
λ, ω), and 34 is a result output terminal. The operating principle of this device will be explained. First, the register group 25 stores the boundary value of the curve interpolation of p(ζ) according to the equation (30), and the register groups 26 and 27 respectively store the coefficients of the interpolation function value of the equation (30). . ω, j _S and j _R input from input terminal 8,
Information on c _S ^-2 and c _R ^-2 , a _jS , and a _jR is stored in registers (groups) 10, 11, 12, and 13, respectively.
Input terminal 9 empty, λ _i , Q _S and Q _R , Q _C , j _u and j _d
Upon receiving the information, the respective values are stored in registers (groups) 14, 15, 16, and 17. Q _S in register 15 and Q _C in register 16 are sent to adder 18 to produce Q _C -Q _S and sent to comparator 19. The comparator 19 is the adder 18
When the value is sent from , the contents of register 15
Q _S and sends the smaller value to the multiplier 20, and if Q _S Q _C −Q _S ,
7 to j _u , and if Q _S > Q _C − Q _S , register 1
7 to send j _d to the adder 23. On the other hand, register 1
Based on the first content j _S , the j _S -th content a _jS is selected from the register group 13 and stored in the register 30 . Further, the adder 21 is connected to the registers 12 and 1.
4, read c _S ^-2 and λ _i ² respectively and obtain c _S ^-2
−λ _i ² is generated and sent to the power multiplier 22. The power multiplier 22 raises the information sent from the adder 21 to the -1/4 power and stores it in the register 30. Multiplier 20 multiplies the minimum value of Q sent from comparator 19 by the content ω of register 10 to create ζ, and sends it to approximate function determiner 24 and multiplier 28 . The adder 23 takes the difference between the contents sent from the registers 11 and 17, and sends the result to the approximate function determiner 24.
send to The approximate function determiner 24 selects the contents of the register groups 26 and 27 and the register 30 based on the following criteria, and selects the contents of the register groups 26 and 27 and the register 30, and selects the contents of the register groups 26 and 27 and the register 30, respectively.
and sent to multiplier 31. 1 ζζ _k , from among the register group 26, β _k
and α _k are selected from the register group 27, and {c _S ⁻² −λ _i ² } ^−1/4 is selected from the register 30 and sent to the multiplier 28, adder 29, and multiplier 31, respectively. 2 If ζ<ζ _k and the content sent from the adder 23 is 0, β _p is sent from the register group 26, α _p is sent from the register group 27, and the register 30
Select a _jS from the multiplier 28,
It is sent to an adder 29 and a multiplier 31. 3 If ζ<ζ _k and the content sent from the adder 23 is other than 0, select β _k and α _k that satisfy ζ _k ζζ _k+1 from register groups 26 and 27, respectively, and register From 30 {c _S ^-2 −
λ _i ² } ^-1/4 are selected, and multiplier 28,
It is sent to an adder 29 and a multiplier 31. Multiplier 28 multiplies the content sent from multiplier 20 by the content sent from register group 26 to create β _k ζ, and sends it to adder 29 . The adder 29 adds α _k sent from the register 27 to the content β _k ζ sent from the multiplier 28 to create (α _k +β _k ζ) and sends it to the multiplier 31 . The multiplier 31 multiplies the content sent from the adder 29 by the content [a _j or {c _S ^-2 −λ _i ² } ^-1/4 ] sent from the register 30, and calculates f _S = f (z _S ;λ
，
ω) and sends it to the register 32. When f _S is sent to the register 32, the calculation regarding the sound source depth is completed, and the next step is to calculate the receiving depth. The contents of registers 11, 12, and 15 are each shifted by one, and the quantities related to the sound source and the quantities related to the wave receiving point are exchanged. Following the same processing procedure as that for the sound source, f _R = f
(z _R ; λ, ω) is calculated and the result is stored in register 3.
It will be sent to 2. The contents of the register 32 are shifted by one, so that f _S and f _R are aligned. Multiplier 3
3, when both f _S and f _R are present in the register 32, the product of the two is taken, and the result h (z _S , z _R ; λ, ω) is outputted from the output terminal 34 to the sound field calculation means 6. send. (Effect of the invention) As explained above, in the present invention, the depth function is
Since the approximation is performed using a stable and accurate function that does not diverge or become unstable even near the turning point, it is possible to avoid the problem of conventional methods, which have poor accuracy near the turning point.
Problems such as the inability to express frequency cut-off effects can be eliminated, and the applicable frequency range can be greatly expanded. In addition, since the approximation of the depth function 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 sonar operation sites. It can also be used in systems.

[Brief explanation of the drawing]

第１図は本発明の実施例を示すブロツク図、第
２図は深度関数計算手段の実施例の説明図、第３
図は音線の４つのタイプを示す説明図、第４図は
関数ｐ（ζ）の振舞を示す説明図である。１，８，９……入力端子、２……初期設定手
段、３……音線パラメータ発生手段、４……基本
経路計算手段、５……深度関数計算手段、６……
音場計算手段、７，３４……出力端子、１０，１
３，１４，１６，１７，２５，２６，２７，３０
……レジスタ（群）、１１，１２，１５，３２…
…シフト・レジスタ、１８，２１，２３，２９…
…加算器、１９……比較器、２０，２８，３１，
３３……乗算器、２２……ベキ乗器、２４……近
似関数判定器。 FIG. 1 is a block diagram showing an embodiment of the present invention, FIG. 2 is an explanatory diagram of an embodiment of the depth function calculation means, and FIG.
The figure is an explanatory diagram showing four types of sound rays, and FIG. 4 is an explanatory diagram showing the behavior of the function p(ζ). 1, 8, 9...Input terminal, 2...Initial setting means, 3...Sound ray parameter generation means, 4...Basic path calculation means, 5...Depth function calculation means, 6...
Sound field calculation means, 7, 34... Output terminal, 10, 1
3, 14, 16, 17, 25, 26, 27, 30
...Register (group), 11, 12, 15, 32...
...Shift register, 18, 21, 23, 29...
...Adder, 19...Comparator, 20, 28, 31,
33... Multiplier, 22... Power multiplier, 24... Approximate function determiner.

Claims

[Scope of Claims] 1. In a sound wave propagation simulation method based on sound ray theory that simulates the sound wave propagation phenomenon in a heterogeneous medium, the sound ray propagation from the turning point of the sound ray to the sound source depth corresponds to the sample value of the sound ray parameter. The phase integral value ζ _R of is calculated, and the amplitude function f(z _S ; λ, ω) of the sound source and the amplitude function f(z _R ;
λ, ω) are calculated as lower-order functions of the phase integral values ζ _S and ζ _R , respectively, and the amplitude value of the depth function is calculated as the product of the amplitude function of the sound source and the amplitude function of the receiving point, and A sound wave propagation sinulation method characterized by performing sound field calculations using amplitude values of functions. 2 The phase integral values ζ _S and ζ _R are respectively ζ _S = ω∫ ^zs _zu [C ^-2 (ξ) - λ ² ] ^1/2 dξ ζ _R = ω∫ ^zR _zu [C ^-2 (ξ) - λ ² ] ^1/2 dξ, and the amplitude function of the sound source f(z _S ; λ,
ω) and the amplitude function f(z _R ; λ, ω) at the receiving point are f(z, ζ; λ, ω) = P(ζ) [C ^-2 (z
) −λ ² ] ^-1/4 , calculated by a _j (z _u , ω) q (ζ) when ζ > ζ ₁ , and P (ζ), q (ζ) when 0≦ζ≦ζ ₁ . are each 1
Calculated by the following formula P (ζ) = α _k + β _k ζ ζ _k ≦ζ≦ζ _k+1 , k = 1, ..., K q (ζ) = α ₀ + β ₀ ζ 0≦ζ≦ζ ₁ , Here, C(z) is the depth characteristic of the speed of sound, z is the depth, λ is the sound ray parameter where the horizontal component of the wave number is expressed by ωλ,
z _u is the turning point depth where C(z _u ) is equal to λ, z _S is the sound source depth, z _R is the receiving point depth, ζ _k (k=1,...,
K) is the interval boundary value in the function approximation of the amplitude function of the sound source and receiving point, α _k + β _k (k = 1, ..., K) is the coefficient value of the approximation function, a _j (z _u , ω) is the rotation 2. The sound wave propagation simulation method according to claim 1, wherein the coefficient is determined from the sound velocity at a point, the sound velocity gradient, and the angular frequency.