特開 2023-125197 音波伝搬計算方法及び音波伝搬計算装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2023125197

(43)【公開日】2023-09-07

(54)【発明の名称】音波伝搬計算方法及び音波伝搬計算装置

(51)【国際特許分類】

   G10K 15/00 20060101AFI20230831BHJP

【ＦＩ】

G10K15/00 L

【審査請求】未請求

【請求項の数】8

【出願形態】ＯＬ

(21)【出願番号】P 2022029172

(22)【出願日】2022-02-28

(71)【出願人】

【識別番号】000000295

【氏名又は名称】沖電気工業株式会社

(74)【代理人】

【識別番号】110001461

【氏名又は名称】弁理士法人きさ特許商標事務所

(72)【発明者】

【氏名】山▲崎▼ 慎太郎

(57)【要約】

【課題】音波の音波伝搬計算における計算時間の短縮を目的とする。
【解決手段】音波伝搬計算方法は、ＰＥ（Parabolic Equation）法を用いて、計算開始距離から所望計算距離までの間、広い距離刻みで音圧を計算する第１音圧計算を行うステップと、所望計算距離から計算終了距離までの間、第１音圧計算で計算された音圧を用いて、第１音圧計算よりも詳細に音圧を計算する第２音圧計算を行うステップと、を含む。
【選択図】図１

【特許請求の範囲】

【請求項1】

ＰＥ（Parabolic Equation）法を用いて、計算開始距離から所望計算距離までの間、広い距離刻みで音圧を計算する第１音圧計算を行うステップと、
前記所望計算距離から計算終了距離までの間、前記第１音圧計算で計算された音圧を用いて、前記第１音圧計算よりも詳細に音圧を計算する第２音圧計算を行うステップと、
を含む音波伝搬計算方法。

【請求項2】

前記第２音圧計算を行うステップは、
前記ＰＥ法を用いて前記第２音圧計算を行うものであり、
前記第１音圧計算と同じ角度範囲に相当し、前記第１音圧計算と異なる距離刻み及び近似次数を用いて、前記第２音圧計算を行う請求項１に記載の音波伝搬計算方法。

【請求項3】

前記第２音圧計算を行うステップは、
前記第１音圧計算で用いられる距離刻みよりも小さい距離刻みを用いて前記第２音圧計算を行う請求項２に記載の音波伝搬計算方法。

【請求項4】

前記第１音圧計算を行うステップは、
角度テーブルにおける最大の距離刻みと対応する次数とを用いて前記第１音圧計算を行うものであり、
前記第２音圧計算を行うステップは、
前記第１音圧計算の終了時における第１計算終了距離と、前記計算終了距離と、表示サンプル数とに基づき計算間隔を算出し、前記計算間隔に応じた前記距離刻みと対応する次数とを用いて前記第２音圧計算を行う請求項２又は３に記載の音波伝搬計算方法。

【請求項5】

前記第２音圧計算を行うステップは、
球面拡散と反射損失を用いて前記第２音圧計算を行うものであり、
前記第１音圧計算で求めた音圧を各深度での音源とし、複数の前記音源の音圧を距離方向に重ね合わせて前記音圧を計算する請求項１に記載の音波伝搬計算方法。

【請求項6】

ＰＥ（Parabolic Equation）法を用いて、計算開始距離から所望計算距離までの間、広い距離刻みで音圧を計算する第１音圧計算を行う第１音圧計算部と、
前記所望計算距離から計算終了距離までの間、前記第１音圧計算で計算された音圧を用いて、前記第１音圧計算よりも詳細に音圧を計算する第２音圧計算を行う第２音圧計算部と、
を備える音波伝搬計算装置。

【請求項7】

前記第２音圧計算部は、
前記ＰＥ法を用いて前記第２音圧計算を行うものであり、
前記第１音圧計算と同じ角度範囲に相当し、前記第１音圧計算と異なる距離刻み及び近似次数を用いて、前記第２音圧計算を行う請求項６に記載の音波伝搬計算装置。

【請求項8】

前記第２音圧計算部は、
球面拡散と反射損失を用いて前記第２音圧計算を行うものであり、
前記第１音圧計算で求めた音圧を各深度での音源とし、複数の前記音源の音圧を距離方向に重ね合わせて前記音圧を計算する請求項６に記載の音波伝搬計算装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、水中音響における音波伝搬を計算する音波伝搬計算方法及び音波伝搬計算装置に関するものである。

【背景技術】

【0002】

海中における音の利用を考える上で、海中の音の伝搬の様子を知ることは重要な要素となっている。海中の音の伝搬を模擬する技術に音波伝搬シミュレーションという技術がある。音波伝搬シミュレーションは音波伝搬理論を基に構成されたプログラム（伝搬モデル）を使って計算されており、現在、頻繁に用いられている伝搬モデルの一つに波動理論を使ったＰＥ（Parabolic Equation）モデルがある（非特許文献１）。

【先行技術文献】

【非特許文献】

【0003】

【非特許文献1】F. B. Jensen 他，“6 Parabolic Equations”， Computational Ocean Acoustics （Second Edition），2011

【発明の概要】

【発明が解決しようとする課題】

【0004】

音波伝搬シミュレーションを活用するにあたって、特定の距離以降の詳細なシミュレーション結果が重要になる場合がある。しかし、ＰＥ法は距離方向に逐次的に音圧を計算していく方式であるために、特定の距離以降だけをピンポイントで詳細に計算することができない。特定の距離以降を細かい距離刻みで計算する場合、すべての距離範囲を細かい距離刻みで計算する必要があり、計算時間が長くなるという問題がある。

【0005】

本発明は、上記のような課題を背景としたものであり、音波の音波伝搬計算における計算時間の短縮を目的とするものである。

【課題を解決するための手段】

【0006】

本発明に係る音波伝搬計算方法は、ＰＥ（Parabolic Equation）法を用いて、計算開始距離から所望計算距離までの間、広い距離刻みで音圧を計算する第１音圧計算を行うステップと、所望計算距離から計算終了距離までの間、第１音圧計算で計算された音圧を用いて、第１音圧計算よりも詳細に音圧を計算する第２音圧計算を行うステップと、を含む。
また、本発明に係る音波伝搬計算装置は、ＰＥ法を用いて、計算開始距離から所望計算距離までの間、広い距離刻みで音圧を計算する第１音圧計算を行う第１音圧計算部と、所望計算距離から計算終了距離までの間、第１音圧計算で計算された音圧を用いて、第１音圧計算よりも詳細に音圧を計算する第２音圧計算を行う第２音圧計算部と、を備える。

【発明の効果】

【0007】

本発明の音波伝搬計算方法及び音波伝搬計算装置によると、精度の異なる二段階の計算を行って音圧を計算することで、所望計算距離以降の詳細なシミュレーションを行いつつ、結果計算時間を短縮することができる。

【図面の簡単な説明】

【0008】

【図1】実施の形態１に係る音波伝搬計算装置における計算概念図である。

【図2】実施の形態１に係る音波伝搬計算装置の概略構成図である。

【図3】実施の形態１に係る近似次数・距離刻み算出部の処理フローである。

【図4】実施の形態１に係る音波伝搬計算装置における音波伝搬計算方法のフローである。

【図5】実施の形態１に係る音波伝搬計算装置における第１音圧計算及び第２音圧計算を説明する図である。

【図6】実施の形態２に係る音波伝搬計算装置の概略構成図である。

【図7】実施の形態２に係る音波伝搬計算装置の計算概念図である。

【図8】実施の形態２に係る音波伝搬計算装置における第１音圧計算及び第２音圧計算を説明する図である。

【図9】実施の形態２に係る簡易音圧算出部の処理フローである。

【図10】経路長の算出に用いられる各パラメータのイメージ図である。

【図11】（ａ）は、近似次数ｍ及び距離刻みΔｒの組み合わせを示すグラフであり、（ｂ）は、近似次数ｍ及び距離刻みΔｒの各組み合わせに対応する計算時間を示すグラフである。

【図12】従来技術に係る音波伝搬計算装置の概略構成図である。

【図13】近似次数ｍと距離刻みΔｒと角度範囲θの関係を示すパラメータテーブルである。

【図14】従来技術に係る近似次数・距離刻み算出部の処理フローである。

【図15】従来技術に係る音波伝搬計算装置における計算概念図である。

【図16】従来技術に係る音波伝搬計算装置における音圧計算を説明する図である。

【発明を実施するための形態】

【0009】

まず、従来のＰＥ（Parabolic Equation）法の高速化理論について説明する。ＰＥ法は、式（１）に示すヘルムホルツ方程式を解くことで音波伝搬を計算している。

【数1】

ここで、ｐは音圧、ｋは波数を表している。

【0010】

円筒拡散と中心波数による位相回転成分を、ハンケル関数を用いてＨ_０ ⁽¹⁾（ｋ_０ｒ）と表し、それ以外の成分をψ（ｒ，ｚ）とすることで、式（２）のようにヘルムホルツ方程式の解を仮定する。Ｈ_０ ⁽¹⁾（ｋ_０ｒ）は距離が分かれば簡単に求まるので、ψ（ｒ，ｚ）が計算できれば任意の距離と深度の音圧が計算できることになる。

【数2】

ここで、ｒは距離、ｚは深度、ｋ_０は中心波数を表している。

【0011】

式（１）に式（２）を代入し、式（３）で定義される微分作用素ｑを導入すると式（４）が得られる。

【数3】

【数4】

ここで、Δｒは距離刻み、ｎは屈折指標である。

【0012】

式（４）を解くにあたって、このままの形ではψ（ｒ＋Δｒ，ｚ）を求めるのは困難である。そこで、式（５）のようにＰａｄｅ近似で指数関数の項を近似することを考える。

【数5】

ここで、ａ_ｊ、ｂ_ｊは、近似係数を表しており、ｊは近似係数のインデックスを表している。

【0013】

この近似係数ａ_ｊ、ｂ_ｊは、距離刻みΔｒと近似次数ｍに依存しており、距離刻みΔｒと近似次数ｍが決まると、一意に求まることが分かっている。式（５）を式（４）に代入することで式（６）が得られる。

【数6】

【0014】

式（６）を解くことで、距離ｒからΔｒだけ進んだ点のψ（ｒ＋Δｒ，ｚ）を計算することができる。式（２）で説明しているように、ψ（ｒ＋Δｒ，ｚ）に対してＨ_０ ⁽¹⁾（ｋ_０ｒ＋Δｒ）を乗算すれば音圧ｐ（ｒ＋Δｒ，ｚ）を求めることができる。このようにして、ＰＥ法では距離方向に逐次的に音圧を計算している。ＰＥ法による計算の初期値（初期音場）の決定と、式（６）の差分化に関しては「石渡 他、“浅海域音場解析のための放物型音波伝搬モデルの開発”、日本音響学会誌、５７（１２）、２００１」に記載されている方法が考えられる。

【0015】

次に、近似による誤差について説明する。近似係数ａ_ｊ、ｂ_ｊが決定すると、式（５）の右辺の絶対値を式（７）のようなｑの関数と考えることができる。

【数7】

【0016】

また、ｑ＝－ｓｉｎ^２θという関係が成り立っていることが分かっているので、式（７）は式（８）のように角度θの関数として書き換えることができる。ここで、θは水平方向を０°とした時の平面波の進行方位を表している。

【数8】

【0017】

式（５）の左辺の絶対値は１なので、式（８）は１から離れるほど誤差が大きくなると言える。さらに、この誤差はθ（－９０°～９０°）の絶対値が大きくなるほど増大することが分かっている。また、この誤差は式（４）からも分かるようにΔｒ（１ステップ）だけ計算が進むごとに累積されていく。つまり、ある距離ｒまでの累積誤差ε（ｒ，θ）は式（９）のように表すことができる。

【数9】

【0018】

距離ｒまでの許容する誤差を設定することで、距離刻みΔｒ、近似次数ｍ、角度θの組み合わせが決まる。このときのθを、ＰＥ法で計算する角度範囲と称する。

【0019】

最後に、角度範囲θを固定した時の距離刻みΔｒ、近似次数ｍの組み合わせについて説明する。近似誤差ｕ（θ）は近似次数ｍが大きくなるほど、又は距離刻みΔｒが小さくなるほど誤差が小さくなる。反対に近似次数ｍが小さくなるほど、又は距離刻みΔｒが大きくなるほど誤差が大きくなることが分かっている。この関係を使うと距離刻みΔｒと近似次数ｍは互いの誤差の増加を打ち消し合うことができる。

【0020】

そのため、同じ値の累積誤差ε（ｒ，θ）でも複数の距離刻みΔｒ及び近似次数ｍの組み合わせが存在することになる。図１１（ａ）は、近似次数ｍ及び距離刻みΔｒの組み合わせを示すグラフであり、図１１（ｂ）は、近似次数ｍ及び距離刻みΔｒの各組み合わせに対応する計算時間を示すグラフである。図１１（ａ）及び（ｂ）は、ｆ＝１０ｋＨｚ、ｒ＝１０ｋｍ、θ＝３０°において、ε（ｒ，θ）＝－４ｄＢを誤差許容した場合の例である。ただし、距離刻みΔｒは波長比である。図１１から、近似次数ｍを下げて距離刻みΔｒを細かくするよりも、近似次数ｍを上げて距離刻みΔｒを粗くする方が同じ誤差の許容値で同じ角度範囲θを計算する場合でも、計算時間が短くなるということが分かる。

【0021】

次に、従来技術に係る音波伝搬計算装置２００について説明する。図１２は、従来技術に係る音波伝搬計算装置２００の概略構成図である。図１２に示すように、従来技術の音波伝搬計算装置２００は、初期条件記憶部１と、計算間隔算出部２と、近似次数・距離刻み算出部３と、近似係数算出部４と、音波伝搬計算部５と、制御部７と、を有している。音波伝搬計算装置２００の各部は、専用の機器又は専用の処理回路により実現される。もしくは、音波伝搬計算装置２００はＣＰＵなどのプロセッサとメモリとを備え、プロセッサがメモリに記憶されるプログラムを実行することにより音波伝搬計算装置２００の各部を実現してもよい。もしくは、音波伝搬計算装置２００の各部を、専用の機器又は専用の回路と、ソフトウェアとの組み合わせにより実現してもよい。

【0022】

初期条件記憶部１は、ＲＡＭ、ＲＯＭ又はフラッシュメモリなどのメモリであり、ユーザーから入力された初期条件を記憶する。初期条件は、環境条件、角度範囲θ、計算開始距離Ｒ_Ｓ、計算終了距離Ｒ_Ｅ、及び表示サンプル数Ｎである。環境条件は、例えば、地形及び音速などである。

【0023】

計算間隔算出部２は、制御部７を介して入力される計算開始距離Ｒ_Ｓ、計算終了距離Ｒ_Ｅ、及び表示サンプル数Ｎから計算間隔Δｒ_ｃを算出して出力する。計算間隔Δｒ_ｃは式（１０）のように計算できる。

【数10】

計算間隔算出部２から出力された計算間隔Δｒ_ｃは、制御部７を介して近似次数・距離刻み算出部３へ入力される。

【0024】

近似次数・距離刻み算出部３は、制御部７を介して入力される角度範囲θ及び計算間隔Δｒ_ｃから、音圧計算に使用する近似次数ｍと距離刻みΔｒとを算出して出力する。図１３は、近似次数ｍと距離刻みΔｒと角度範囲θの関係を示すパラメータテーブルである。近似次数・距離刻み算出部３は、図１３に示すパラメータテーブルを参照し、入力された角度範囲θ及び計算間隔Δｒ_ｃを用いて、近似次数ｍと距離刻みΔｒとを算出する。

【0025】

図１４は、従来技術に係る近似次数・距離刻み算出部３の処理フローである。図１４に示すように、近似次数・距離刻み算出部３は、まず、パラメータテーブルから角度範囲θに対応した行（角度テーブル）を抽出する（Ｓ１０１）。そして、近似次数・距離刻み算出部３は、抽出した角度テーブルから、計算間隔Δｒ_ｃ未満の要素を抽出し、新たな角度テーブルを生成する（Ｓ１０２）。近似次数・距離刻み算出部３は、生成した新たな角度テーブルから最大の距離刻みΔｒとそれに対応した近似次数ｍを取得する（Ｓ１０３）。そして、近似次数・距離刻み算出部３は、抽出した距離刻みΔｒと近似次数ｍを出力する（Ｓ１０４）。近似次数・距離刻み算出部３から出力された近似次数ｍと距離刻みΔｒは、制御部７を介して近似係数算出部４へ入力される。

【0026】

図１２に戻って、近似係数算出部４は、制御部７を介して入力される近似次数ｍと距離刻みΔｒから、近似係数ａ_ｊ、ｂ_ｊを算出して出力する。近似係数算出部４では、Ｐａｄｅ近似の元々の考え方である、近似した有理関数（分母と分子がｍ次の多項式）が近似元の関数の２ｍ次のテイラー展開に等しくなるという関係を用いて、距離刻みΔｒと近似次数ｍから近似係数ａ_ｊ、ｂ_ｊを算出する。近似係数算出部４から出力された近似係数ａ_ｊ、ｂ_ｊは、制御部７を介して音波伝搬計算部５へ入力される。

【0027】

音波伝搬計算部５は、制御部７を介して入力される近似係数ａ_ｊ、ｂ_ｊ、環境条件、計算開始距離Ｒ_Ｓ、計算終了距離Ｒ_Ｅから、計算距離系列Ｒｎｇ及び音圧ｐ（ｒ，ｚ）を算出して出力する。音波伝搬計算部５は、式（６）で示した微分方程式に対して差分化を用いることで、音圧ｐ（ｒ，ｚ）を距離方向に逐次的に計算する。計算距離系列Ｒｎｇは、計算を実施した距離の系列を表しており、初項を計算開始距離Ｒ_Ｓ、末項を計算終了距離Ｒ_Ｅ、交差を距離刻みΔｒとした時の等差数列で与えられる配列である。

【0028】

制御部７は、初期条件記憶部１、計算間隔算出部２、近似次数・距離刻み算出部３、近似係数算出部４、及び音波伝搬計算部５にそれぞれ接続されており、各部の入出力を制御する。

【0029】

図１５は、従来技術に係る音波伝搬計算装置２００における計算概念図である。図１５に示すように、従来技術に係る音波伝搬計算装置２００では、ＰＥ法を用いて、計算間隔Δｒ_ｃに応じた距離刻みΔｒと近似次数ｍで、全距離範囲（計算開始距離Ｒ_Ｓから計算終了距離Ｒ_Ｅまで）の音圧を計算する。

【0030】

図１６は、従来技術に係る音波伝搬計算装置２００における音圧計算を説明する図である。図１６においては、各部を仲介する制御部７は省略している。まず、ユーザーによって、初期条件が初期条件記憶部１に入力される。そして、初期条件記憶部１から計算開始距離Ｒ_Ｓ、計算終了距離Ｒ_Ｅ及び表示サンプル数Ｎが計算間隔算出部２に入力され、計算間隔算出部２から計算間隔Δｒ_ｃが出力される。そして、計算間隔算出部２で算出された計算間隔Δｒ_ｃと、初期条件記憶部１に記憶された角度範囲θとが近似次数・距離刻み算出部３に入力され、近似次数・距離刻み算出部３から距離刻みΔｒと近似次数ｍとが出力される。

【0031】

次に、近似次数・距離刻み算出部３から出力された距離刻みΔｒと近似次数ｍとが近似係数算出部４に入力され、近似係数算出部４から近似係数ａ_ｊ、ｂ_ｊが出力される。最後に、近似係数算出部４にて算出された近似係数ａ_ｊ、ｂ_ｊと、初期条件記憶部１に記憶された環境条件、計算開始距離Ｒ_Ｓ、及び計算終了距離Ｒ_Ｅとが音波伝搬計算部５に入力され、音波伝搬計算部５から計算距離系列Ｒｎｇと音圧ｐ（ｒ，ｚ）とが出力される。

【0032】

従来技術に係る音波伝搬計算装置２００においては、距離方向に逐次的に音圧を計算していくため、特定の距離以降だけをピンポイントで詳細に計算することができない。従来技術に係る音波伝搬計算装置２００において、特定の距離以降を細かい距離刻みで計算する場合、すべての距離範囲を細かい距離刻みで計算する必要があり、計算時間が長くなってしまう。このような問題を解決する本発明の実施の形態について、以下に説明する。

【0033】

実施の形態１．
図１は、実施の形態１に係る音波伝搬計算装置１００における計算概念図である。上記のように、従来方式では全距離範囲の計算を一回の計算（１種類の距離刻みΔｒ、近似次数ｍ）で実行していた。これに対し、本実施の形態の音波伝搬計算装置１００は、特定の距離までの音圧計算を行う第１音圧計算と、特定の距離以降の音圧計算を行う第２音圧計算の二段階の音波伝搬計算を行う構成を有している。

【0034】

図１に示すように、一段階目の第１音圧計算では、角度範囲θに対応した最大の距離刻みΔｒ^（１）、すなわち粗い距離刻みで、計算開始距離Ｒ_Ｓから所望計算距離Ｒ_Ｐの直前までが計算される。一方、二段階目の第２音圧計算では第１音圧計算の最後の距離である第１計算終了距離Ｒ_Ｅ ^（１）の音圧ｐ（Ｒ_Ｅ ^（１），z）を初期音場として用いて、計算間隔Δｒ_ｃに応じた細かい距離刻みΔｒ^（２）で所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでが計算される。

【0035】

図２は、実施の形態１に係る音波伝搬計算装置１００の概略構成図である。図２に示すように、本実施の形態の音波伝搬計算装置１００は、初期条件記憶部１０と、計算間隔算出部２０と、近似次数・距離刻み算出部３０と、近似係数算出部４０と、音波伝搬計算部５０と、計算結果記憶部６０と、制御部７０と、を有している。音波伝搬計算装置１００の各部は、専用の機器又は専用の処理回路により実現される。もしくは、音波伝搬計算装置１００はＣＰＵなどのプロセッサとメモリとを備え、プロセッサがメモリに記憶されるプログラムを実行することにより音波伝搬計算装置１００の各部を実現してもよい。もしくは、音波伝搬計算装置１００の各部を、専用の機器又は専用の回路と、ソフトウェアとの組み合わせにより実現してもよい。

【0036】

初期条件記憶部１０は、ＲＡＭ、ＲＯＭ又はフラッシュメモリなどのメモリであり、ユーザーから入力された初期条件を記憶する。初期条件は、環境条件、角度範囲θ、計算開始距離Ｒ_Ｓ、計算終了距離Ｒ_Ｅ、所望計算距離Ｒ_Ｐ、及び表示サンプル数Ｎである。所望計算距離Ｒ_Ｐは、第１音圧計算から第２音圧計算へと切替える距離である。また、初期条件記憶部１０は、さらに後述するフラグ、音波伝搬計算装置１００の処理に必要なプログラム、パラメータ及びデータなどを記憶している。

【0037】

計算間隔算出部２０は、制御部７０を介して入力される第１計算終了距離Ｒ_Ｅ ^（１）、計算終了距離Ｒ_Ｅ、及び表示サンプル数Ｎから計算間隔Δｒ_ｃを算出して出力する。第１計算終了距離Ｒ_Ｅ ^（１）は、第１音圧計算を行った距離であり、所望計算距離Ｒ_Ｐと略同じである。計算間隔算出部２０は、上記の式（１０）を用いて計算間隔Δｒ_ｃを算出する。ただし、式（１０）における計算開始距離Ｒ_Ｓに替えて第１計算終了距離Ｒ_Ｅ ^（１）が用いられる。計算間隔算出部２０から出力された計算間隔Δｒ_ｃは、制御部７０を介して近似次数・距離刻み算出部３０へ入力される。

【0038】

近似次数・距離刻み算出部３０は、制御部７０を介して入力される角度範囲θ、計算間隔Δｒ_ｃ、及びフラグに基づき、音圧計算に使用する近似次数ｍと距離刻みΔｒを算出して出力する。フラグは、第１音圧計算か第２音圧計算かを判断するフラグであり、フラグが「１」場合は第１音圧計算であり、「１」以外の場合は第２音圧計算であることを示す。

【0039】

図３は、実施の形態１に係る近似次数・距離刻み算出部３０の処理フローである。図３に示すように、近似次数・距離刻み算出部３０は、まずパラメータテーブルから角度範囲θに対応した行（角度テーブル）を抽出する（Ｓ１１）。そして、近似次数・距離刻み算出部３０は、フラグが１であるか否かを判断する（Ｓ１２）。フラグが１である場合（Ｓ１２：ＹＥＳ）、すなわち第１音圧計算の場合、近似次数・距離刻み算出部３０は、抽出した角度テーブルから最大の距離刻みΔｒとそれに対応した近似次数ｍを取得する（Ｓ１３）。

【0040】

一方、フラグが１でない場合（Ｓ１２：ＮＯ）、すなわち第２音圧計算の場合、近似次数・距離刻み算出部３０は、抽出した角度テーブルから、計算間隔Δｒ_ｃ未満の要素を抽出し、新たな角度テーブルを生成する（Ｓ１４）。そして、近似次数・距離刻み算出部３は、算出した新たな角度テーブルから最大の距離刻みΔｒとそれに対応した近似次数ｍを取得する（Ｓ１３）。最終的に近似次数・距離刻み算出部３は、距離刻みΔｒと近似次数ｍを出力する（Ｓ１５）。これにより、第１音圧計算においては、できるだけ広い距離刻みで計算が行われるのに対して、第２音圧計算においては、ユーザーが所望する計算間隔Δｒ_ｃに応じた、第１音圧計算よりも小さい距離刻みで計算を行うことができる。近似次数・距離刻み算出部３０から出力された近似次数ｍと距離刻みΔｒは、制御部７０を介して近似係数算出部４へ入力される。

【0041】

近似係数算出部４０は、制御部７０を介して入力される近似次数ｍと距離刻みΔｒから、近似係数ａ_ｊ、ｂ_ｊを算出して出力する。近似係数算出部４０では、Ｐａｄｅ近似の元々の考え方である、近似した有理関数（分母と分子がｍ次の多項式）が近似元の関数の２ｍ次のテイラー展開に等しくなるという関係を用いて、距離刻みΔｒと近似次数ｍから近似係数ａ_ｊ、ｂ_ｊを算出する。近似係数算出部４０から出力された近似係数ａ_ｊ、ｂ_ｊは、制御部７０を介して音波伝搬計算部５０へ入力される。

【0042】

音波伝搬計算部５０は、第１音圧計算においては、制御部７０を介して入力される近似係数ａ_ｊ、ｂ_ｊ、環境条件、計算開始距離Ｒ_Ｓ、所望計算距離Ｒ_Ｐから、計算距離系列Ｒｎｇ及び音圧ｐ（ｒ，ｚ）を算出して出力する。第１音圧計算において音波伝搬計算部５０から出力された音圧ｐ（ｒ，ｚ）及び計算距離系列Ｒｎｇは、制御部７０を介して計算結果記憶部６０に記憶される。

【0043】

また、音波伝搬計算部５０は、第２音圧計算においては、制御部７０を介して入力される近似係数ａ_ｊ、ｂ_ｊ、初期音場（第１音圧計算で最後に計算した音圧ｐ（Ｒ_Ｅ ^（１），ｚ））、環境条件、計算終了距離Ｒ_Ｅから、音圧ｐ（ｒ，ｚ）及び計算距離系列Ｒｎｇを算出して出力する。音波伝搬計算部５０は、式（６）で示した微分方程式に対して差分化を用いることで、音圧ｐ（ｒ，ｚ）を距離方向に逐次的に計算する。

【0044】

計算結果記憶部６０は、ＲＡＭ、ＲＯＭ又はフラッシュメモリなどのメモリであり、制御部７０を介して入力された第１計算終了距離Ｒ_Ｅ ^（１）、音圧ｐ（ｒ，ｚ）及び計算距離系列Ｒｎｇを保存する。第２音圧計算時には、制御部７０によって、初期音場（第１音圧計算で最後に計算した音圧ｐ（Ｒ_Ｅ ^（１），ｚ））と第１計算終了距離Ｒ_Ｅ ^（１）とが抽出される。なお、計算結果記憶部６０は、初期条件記憶部１０と別々のメモリにもうけられてもよいし、同じメモリであってもよい。

【0045】

制御部７０は、初期条件記憶部１０、計算間隔算出部２０、近似次数・距離刻み算出部３０、近似係数算出部４０、音波伝搬計算部５０と、及び計算結果記憶部６０にそれぞれ接続されており、各部の入出力を制御する。

【0046】

図４は、実施の形態１に係る音波伝搬計算装置１００における音波伝搬計算方法のフローである。図４に示すように、まず、角度テーブルの最大の距離刻みΔｒとそれに対応した近似次数ｍを用いた第１音圧計算が行われる（Ｓ１）。第１音圧計算では、音波伝搬計算部５０によって、０＜Ｒ_Ｐ－Ｎ^（１）Δｒ^（１）≦Δｒ^（１）を満たすようなＮ^（１）が算出される（Ｓ１０１）。Ｎ^（１）Δｒ^（１）は、第１計算終了距離Ｒ_Ｅ ^（１）である。そして、音波伝搬計算部５０は、距離ステップを更新しながらＮ^（１）回、音圧ｐ（ｒ，z）を計算する（Ｓ１０２）。これにより、所望計算距離Ｒ_Ｐまで第１音圧計算が行われる。

【0047】

続いて、ユーザーが所望する計算間隔Δｒ_ｃに応じた距離刻みΔｒとそれに対応した近似次数ｍを用いた第２音圧計算が行われる（Ｓ２）。そして、第２音圧計算が終了すると、音圧ｐ（ｒ，z）と計算距離系列Ｒｎｇが出力される（Ｓ３）。

【0048】

図５は、実施の形態１に係る音波伝搬計算装置１００における第１音圧計算及び第２音圧計算を説明する図である。図５において、破線で囲まれた各部が第１音圧計算を行う第１音圧計算部であり、一点鎖線で囲まれた各部が第２音圧計算を行う第２音圧計算部である。図５に示すように、まず、ユーザーによって初期条件が入力され、初期条件記憶部１０に記憶される。初期条件は、環境条件、角度範囲θ、計算開始距離Ｒ_Ｓ、計算終了距離Ｒ_Ｅ、所望計算距離Ｒ_Ｐ、及び表示サンプル数Ｎである。初期条件は、ユーザーが音波伝搬計算装置１００の操作部を操作することで入力されてもよいしユーザーから送信された初期条件を音波伝搬計算装置１００の通信部で受信することで入力されてもよい。

【0049】

そして、第１音圧計算において、初期条件記憶部１０から角度範囲θとフラグとが近似次数・距離刻み算出部３０に入力される。そして、近似次数・距離刻み算出部３０から角度範囲θに応じた最大の距離刻みΔｒ^（１）と、それに対応する近似次数ｍ^（１）とが出力される。

【0050】

続いて、近似次数・距離刻み算出部３０から出力された距離刻みΔｒ^（１）と近似次数ｍ^（１）が近似係数算出部４０に入力され、近似係数算出部４０から第１音圧計算で用いるための近似係数ａ_ｊ ^（１）、ｂ_ｊ ^（１）が出力される。

【0051】

続いて、近似係数算出部４０から出力された近似係数ａ_ｊ ^（１）、ｂ_ｊ ^（１）と、初期条件記憶部１０に記憶された環境条件、計算開始距離Ｒ_Ｓ、及び所望計算距離Ｒ_Ｐとが、音波伝搬計算部５０に入力される。そして、音波伝搬計算部５０においてＮ^（１）ステップまでの計算距離系列Ｒｎｇ^（１）と音圧ｐ（ｒ，z）が計算される。

【0052】

音波伝搬計算部５０において計算された計算距離系列Ｒｎｇ^（１）及び音圧ｐ（ｒ，z）と、第１計算終了距離Ｒ_Ｅ ^（１）とは、計算結果記憶部６０に記憶される。この音圧（ｒ，z）は、第２音圧計算で用いられる初期音圧ｐ（Ｒ_Ｅ ^（１），ｚ）である。その後、制御部７０によって、初期条件記憶部１０に記憶されるフラグが「１」以外（例えば「０」）に変更され、第１音圧計算が終了する。

【0053】

第２音圧計算においては、計算結果記憶部６０に記憶された第１計算終了距離Ｒ_Ｅ ^（１）と、初期条件記憶部１０に記憶された計算終了距離Ｒ_Ｅ及び表示サンプル数Ｎとが計算間隔算出部２０へ入力され、計算間隔算出部２０から計算間隔Δｒｃが出力される。

【0054】

続いて、計算間隔算出部２０から出力された計算間隔Δｒｃと、初期条件記憶部１０に記憶された角度範囲θ及びフラグとが、近似次数・距離刻み算出部３０へ入力される。そして、近似次数・距離刻み算出部３０から計算間隔Δｒｃに応じた距離刻みΔｒ^（２）と近似次数ｍ^（２）とが出力される。

【0055】

続いて、近似次数・距離刻み算出部３０から出力された距離刻みΔｒ^（２）と近似次数ｍ^（２）が近似係数算出部４０へ入力され、近似係数算出部４０から第２音圧計算で用いるための近似係数ａ_ｊ ^（２）、ｂ_ｊ ^（２）が出力される。

【0056】

続いて、近似係数算出部４０から出力された近似係数ａ_ｊ ^（２）、ｂ_ｊ ^（２）と、初期条件記憶部１０に記憶された環境条件及び計算終了距離Ｒ_Ｅと、計算結果記憶部６０に記憶された初期音圧ｐ（Ｒ_Ｅ ^（１），ｚ）及び第１計算終了距離Ｒ_Ｅ ^（１）とが音波伝搬計算部５０に入力される。そして、音波伝搬計算部５０によって所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでの音圧が詳細に計算され、音圧ｐ（ｒ，ｚ）と計算距離系列Ｒｎｇ^（２）とがユーザーに出力される。

【0057】

以上のように、本実施の形態の音波伝搬計算装置１００によれば、所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでを細かい距離刻みで計算し、それ以前の距離を粗い距離刻みで計算することになるので、従来方式のように全距離範囲を細かく計算する必要がなくなる。同じ角度範囲θの計算をする場合でも、距離刻みが粗いほど、計算時間が短くなるため、細かい距離刻みで計算する距離範囲が少ない本実施の形態の構成は、特定の距離以降の詳細なシミュレーションを行いつつ、計算時間を短縮することができる。

【0058】

実施の形態２．
図６は、実施の形態２に係る音波伝搬計算装置１００Ａの概略構成図である。実施の形態２では、実施の形態１の計算間隔算出部２０に替えて簡易音圧算出部８０を備える点において実施の形態１と相違する。音波伝搬計算装置１００Ａにおけるその他の構成は、実施の形態１と同じである。

【0059】

簡易音圧算出部８０は、初期音源、環境条件、計算終了距離Ｒ_Ｅ、表示サンプル数Ｎ、第１計算終了距離Ｒ_Ｅ ^（１）を用いて、任意の距離の音圧ｐ（ｒ，ｚ）と計算距離系列Ｒｎｇ^（２）を算出する。初期音源は、第１音圧計算で最後に計算した距離の計算深度毎の音圧ｐ（Ｒ_Ｅ ^（１），ｚ^ｉ）を、それぞれ各深度に対応した個々の音源と考えた場合の音源である。

【0060】

本実施の形態の音波伝搬計算装置１００Ａにおける全体の処理の流れについて説明する。本実施の形態の音波伝搬計算装置１００Ａにおいても、実施の形態１と同様に二段階の音波伝搬計算を実施するものであり、音波伝搬計算装置１００Ａにおける音波伝搬計算フローは、図４と同じである。

【0061】

図７は、実施の形態２に係る音波伝搬計算装置１００Ａの計算概念図である。本実施の形態では、第２音圧計算の計算方法が実施の形態１と相違する。本実施の形態の第２音圧計算では、簡易音圧算出部８０によって、等音速且つ、距離非依存環境を仮定し、球面拡散と反射損失のみで音圧ｐ（ｒ，ｚ）を計算する。

【0062】

図８は、実施の形態２に係る音波伝搬計算装置１００Ａにおける第１音圧計算及び第２音圧計算を説明する図である。図８において、破線で囲まれた各部が第１音圧計算を行う第１音圧計算部であり、一点鎖線で囲まれた簡易音圧算出部８０が第２音圧計算を行う第２音圧計算部である。実施の形態２における第１音圧計算は、実施の形態１と同じであるため、説明を省略する。

【0063】

本実施の形態の第２音圧計算では、初期条件記憶部１０に記憶された計算終了距離Ｒ_Ｅ、表示サンプル数Ｎ、及び環境条件と、計算結果記憶部６０に記憶された第１計算終了距離Ｒ_Ｅ ^（１）及び初期音圧ｐ（Ｒ_Ｅ ^（１），ｚ）とが、簡易音圧算出部８０へ入力される。そして、簡易音圧算出部８０は、所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでの音圧を計算し、音圧ｐ（ｒ，ｚ）と計算距離系列Ｒｎｇ^（２）を出力する。

【0064】

簡易音圧算出部８０は、球面拡散と反射損失を使って音圧を算出する。放射角θに対する音圧ｐ_θを算出するための式（１１）は下記のとおりである。

【数11】

ここで、θは放射角、ｐ_０は音源の音圧、ｌは経路長、ｃは音速、bottom_lossθは放射角θに対応した反射損失であり、海底反射していない場合はbottom_loss＝１となる。

【0065】

次に簡易音圧算出部８０の動作について説明する。簡易音圧算出部８０は、第１音圧計算で求めた鉛直音場を各深度での音源とし、複数の音源からの重ね合わせによって音圧を計算する。より詳しくは、簡易音圧算出部８０は、音源から受波点までの音圧を経路毎（直接経路、海面反射経路、海底反射経路）に計算する。そして、簡易音圧算出部８０は、経路毎の音圧を初期音源の数だけ計算し、最後に重ね合わせることで、受波点での音圧を求める。簡易音圧算出部８０は、深度方向に対して受波点での音圧の計算を行うことで、ある距離における音圧を計算する。そして、簡易音圧算出部８０は、ある距離での音圧を距離方向に計算していくことで、所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでの音圧を求めることができる。

【0066】

図９は、実施の形態２に係る簡易音圧算出部８０の処理フローである。簡易音圧算出部８０は、まず、所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでの計算に用いる計算距離系列Ｒｎｇ^（２）を算出する（Ｓ２１）。第１計算終了距離Ｒ_Ｅ ^（１）、所望計算距離Ｒ_Ｐ、計算終了距離Ｒ_Ｅ、表示サンプル数Ｎ、計算サンプル番号ｎとすると、計算距離系列Ｒｎｇ^（２）のｎ番目の要素ｒ_ｎは下記の式（１２）のように算出することができる。

【数12】

【0067】

式（１２）におけるｒ_０は、下記の式（１３）で求められる。

【数13】

【0068】

続いて、簡易音圧算出部８０は、音源から受波点までの直接経路、海面反射経路及び海底反射経路の各経路長を算出する（Ｓ２２）。音源深度をｚ_ｉ（ｉは音源深度のインデックス）、受波深度をｚ_ｒ、水深をｚ_ｈ、計算距離系列Ｒｎｇ^（２）のｎ番目の要素をｒ_ｎ、直接経路の放射角をθ_ＤＰ、海面反射経路の放射角をθ_ＳＲ、海底反射経路の放射角をθ_ＢＲとする。図１０は、経路長の算出に用いられる各パラメータのイメージ図である。経路長ｌは音源から受波点に到達する固有音線の放射角θを用いることで、式（１４）のように求めることができる。

【数14】

【0069】

以下に経路毎の放射角と経路長の算出方法を示す。
直接経路（Direct Pass）の経路長ｌ_ＤＰｉは、下記の式（１５）で算出される。

【数15】

【0070】

ここで、θ_ＤＰｉは直接経路の放射角であり、下記の式（１６）で求められる。

【数16】

【0071】

海面反射経路（Surface Reflection）の経路長ｌ_ＳＲｉは、下記の式（１７）で算出される。

【数17】

【0072】

ここで、θ_ＳＲｉは海面反射経路の放射角であり、下記の式（１８）で求められる。

【数18】

【0073】

海底反射経路（Bottom Reflection）の経路長ｌ_ＢＲｉは、下記の式（１９）で算出される。

【数19】

【0074】

ここで、θ_ＢＲｉは海底反射経路の放射角であり、下記の式（２０）で求められる。

【数20】

【0075】

続いて、簡易音圧算出部８０は、求めた受波点までの経路毎の経路長を使って、球面拡散と反射損失から、経路毎の音圧ｐ_ｉを算出する（Ｓ２３）。音源の音圧は第１音圧計算の最後に計算した距離の計算深度毎の音圧で与えられるので、ｐ（Ｒ_Ｅ ^（１），ｚ_ｉ）と表せる。受波点での音圧をｐ_ｉ（ｉは音源深度のインデックス）とすると、下記の式（２１）のように計算することができる。

【数21】

ここで、ｃは音速、bottom_lossはθ_ＢＲｉに対応した海底での反射損失である。簡易音圧算出部８０は、ステップＳ２２及びＳ２３の処理を初期音源の数だけループさせる。

【0076】

続いて、簡易音圧算出部８０は、経路毎の音圧ｐ_ｉを音源深度毎に加算することで、初期音源によって作られる受波点（ｒ，ｚ）での音圧ｐ（ｒ，ｚ）を算出する（Ｓ２４）。音圧ｐ（ｒ，ｚ）は、下記の式（２２）のように計算できる。

【数22】

ここで、ｎｓは音源の総数である。

【0077】

そして、簡易音圧算出部８０は、ステップＳ２４の処理を受波深度がｍとなるまでループさせる。そして、これを距離方向に計算していくことで、所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでの音圧ｐ（ｒ，ｚ）が求められる。

【0078】

以上のように、実施の形態２に係る音波伝搬計算装置１００Ａによれば、第２音圧計算を簡易的な音圧計算に置き換えることで、所望計算距離Ｒ_Ｐから計算終了距離Ｒ_Ｅまでの音圧を解析的に計算することができる。音圧を解析的に計算することが可能となることで、ピンポイントの座標の音圧計算が可能になる。また、実施の形態１の音波伝搬計算装置１００においてピンポイントの座標の音圧を計算する場合、計算サンプル間の補間が必要になるが、実施の形態２の音波伝搬計算装置１００Ａによれば、補間なしでピンポイントの座標の音圧を求めることができる。

【0079】

以上が本発明の実施の形態の説明であるが、本発明は、上記の実施の形態の構成に限定されるものではなく、その技術的思想の範囲内で様々な変形又は組み合わせが可能である。例えば、実施の形態２では、第２音圧計算を行う際に、第１音圧計算で計算した鉛直音場（一段階目の最後に計算した距離の計算深度毎の音圧ｐ（Ｒ_Ｅ ^（１），ｚ_ｉ））を深度方向に並んだ音源と捉えることで、球面拡散と反射損失を計算し音圧を求めた。この「鉛直音場を深度方向に並んだ複数の音源と捉える」という考え方を用いれば、簡易音圧算出部８０に替えて、音線モデル又はノーマルモードモデルなど別の伝搬モデルを使って、第２音圧計算を行うこともできる。

【符号の説明】

【0080】

１ 初期条件記憶部、２ 計算間隔算出部、３ 近似次数・距離刻み算出部、４ 近似係数算出部、５ 音波伝搬計算部、７ 制御部、１０ 初期条件記憶部、２０ 計算間隔算出部、３０ 近似次数・距離刻み算出部、４０ 近似係数算出部、５０ 音波伝搬計算部、６０ 計算結果記憶部、７０ 制御部、８０ 簡易音圧算出部、１００、１００Ａ、２００ 音波伝搬計算装置。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】