特許 7533277 音波伝搬特性算出方法および音波伝搬特性算出装置

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

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-08-05

(45)【発行日】2024-08-14

(54)【発明の名称】音波伝搬特性算出方法および音波伝搬特性算出装置

(51)【国際特許分類】

   G01S 7/526 20060101AFI20240806BHJP

【ＦＩ】

G01S7/526 Z

【請求項の数】 6

(21)【出願番号】P 2021029340

(22)【出願日】2021-02-26

(65)【公開番号】P2022130771

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

【審査請求日】2023-11-09

(73)【特許権者】

【識別番号】000000295

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

(74)【代理人】

【識別番号】110001461

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

(72)【発明者】

【氏名】森下 到

【審査官】佐藤 宙子

(56)【参考文献】

【文献】特開２０１６－０３１２３５（ＪＰ，Ａ）

【文献】特開２００３－１６１７７３（ＪＰ，Ａ）

【文献】特開２００５－０８３９３２（ＪＰ，Ａ）

【文献】米国特許第１０３７９２１８（ＵＳ，Ｂ１）

【文献】米国特許出願公開第２０１４／０３７２０４９（ＵＳ，Ａ１）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｓ １／７２－ １／８２

Ｇ０１Ｓ ３／８０－ ３／８６

Ｇ０１Ｓ ５／１８－ ５／３０

Ｇ０１Ｓ ７／５２－ ７／６４

Ｇ０１Ｓ １５／００－１５／９６

(57)【特許請求の範囲】

【請求項1】

近似式を用いて音圧を算出する音波伝搬特性算出方法であって、
予め決めた前記音圧の計算時間である許容時間で計算が可能な、近似の次数と、距離方向に前記音圧を逐次的に計算する際の距離間隔との複数の組み合わせに対して、それぞれ計算精度を評価し、
前記計算精度が最も高い前記次数および前記距離間隔の組み合わせを選択し、
選択された前記組み合わせにおける前記次数および前記距離間隔に基づき、前記近似式の係数である近似係数を算出し、
算出された前記近似係数に基づき、前記音圧を算出する
音波伝搬特性算出方法。

【請求項2】

前記近似式における真値と近似値との位相誤差に基づき、前記計算精度を評価する
請求項１に記載の音波伝搬特性算出方法。

【請求項3】

前記位相誤差が予め設定された閾値以下であるか否かに基づき、前記計算精度を評価する
請求項２に記載の音波伝搬特性算出方法。

【請求項4】

前記次数を１から順に１ずつ増加させながら、対応する前記距離間隔を算出する
請求項１～３のいずれか一項に記載の音波伝搬特性算出方法。

【請求項5】

前記複数の組み合わせにおいて、前記距離間隔は、前記次数に比例して増加する
請求項１～４のいずれか一項に記載の音波伝搬特性算出方法。

【請求項6】

近似式を用いて音圧を算出する音波伝搬特性算出装置であって、
近似の次数および予め決めた前記音圧の計算時間である許容時間を出力する制御装置と、
前記次数および前記許容時間に基づき、前記近似式の係数である近似係数と、距離方向に前記音圧を逐次的に計算する際の距離間隔とを出力するパラメータ算出装置と、
前記近似係数に基づき、前記音圧を算出する音圧算出装置と
を備え、
前記パラメータ算出装置は、
前記許容時間で計算が可能な、前記次数および前記距離間隔の複数の組み合わせに対して、それぞれ計算精度を評価し、
前記次数および前記距離間隔に基づき前記近似係数を算出し、
前記制御装置は、
前記計算精度が最も高い前記次数および前記距離間隔の組み合わせを選択し、
前記音圧算出装置は、
選択された前記組み合わせに対応する前記近似係数に基づき、前記音圧を算出する
音波伝搬特性算出装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、海洋中の音波伝搬特性を算出する音波伝搬特性算出方法および音波伝搬特性算出装置に関するものである。

【背景技術】

【0002】

従来、海洋中の音波伝搬特性の模擬は、音を使った海洋環境特性の推定やターゲットから放射される音を使ったターゲット位置の推定等、様々な目的のために行われている（例えば、非特許文献１参照）。音波伝搬特性を模擬する方法の一つとして、例えば、非特許文献２に記載されているように、「放物型方程式（Parabolic Equation）法（以下、「ＰＥ法」と適宜称する）」が知られている（例えば、非特許文献２参照）。ＰＥ法は、距離および深度に対する音圧および波数が満足するヘルムホルツ方程式から得られる関係に基づいて、距離について音圧を逐次的に求めるものである。

【0003】

ＰＥ法のうち、近似を用いて音圧を求める方法として、「Ｓｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法」が知られている（例えば、非特許文献３参照）。Ｓｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法では、近似の度合いを設定することにより、算出される音圧の精度と、音圧を算出する際の計算量を調整することができる。

【先行技術文献】

【非特許文献】

【0004】

【文献】T. F. Duda 他， “Computational Acoustics in Oceanography: The Research Roles of Sound Field Simulations,” Acoustics Today, 15(3), 2019

【文献】F. B. Jensen 他, “6 Parabolic Equations,” in Computational Ocean Acoustics (Second Edition), 2011

【文献】石渡 他， “浅海域音場解析のための放物型音波伝搬モデルの開発，” 日本音響学会誌, 57(12), 2001

【発明の概要】

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

【0005】

ところで、Ｓｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法では、近似度合いの設定と、算出される音圧の精度との関係を定量的に把握することができない。そのため、精度よく音圧を算出することを目的として必要以上に計算量を増やしたり、計算量を減らすことによって模擬する音波伝搬特性の俯仰角方向に対する角度範囲が狭くなり、必要な角度範囲を達成できなかったりするという課題があった。すなわち、従来のＳｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法では、音圧計算の際の精度と計算量とのトレードオフを図ることが困難である。

【0006】

そこで、音圧計算に対して許容できる計算時間を満足しながら、俯仰角方向に対する最大の角度範囲で音波伝搬特性を精度よく算出することができる音波伝搬特性算出方法および音波伝搬特性算出装置が望まれている。

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

【0007】

本発明に係る音波伝搬特性算出方法は、近似式を用いて音圧を算出する音波伝搬特性算出方法であって、予め決めた前記音圧の計算時間である許容時間で計算が可能な、近似の次数と、距離方向に前記音圧を逐次的に計算する際の距離間隔との複数の組み合わせに対して、それぞれ計算精度を評価し、前記計算精度が最も高い前記次数および前記距離間隔の組み合わせを選択し、選択された前記組み合わせにおける前記次数および前記距離間隔に基づき、前記近似式の係数である近似係数を算出し、算出された前記近似係数に基づき、前記音圧を算出するものである。

【0008】

また、本発明に係る音波伝搬特性算出装置は、近似式を用いて音圧を算出する音波伝搬特性算出装置であって、近似の次数および予め決めた前記音圧の計算時間である許容時間を出力する制御装置と、前記次数および前記許容時間に基づき、前記近似式の係数である近似係数と、距離方向に音圧を逐次的に計算する際の距離間隔とを出力するパラメータ算出装置と、前記近似係数に基づき、音圧を算出する音圧算出装置とを備え、前記パラメータ算出装置は、前記許容時間で計算が可能な、前記次数および前記距離間隔の複数の組み合わせに対して、それぞれ計算精度を評価し、前記次数および前記距離間隔に基づき前記近似係数を算出し、前記制御装置は、前記計算精度が最も高い前記次数および前記距離間隔の組み合わせを選択し、前記音圧算出装置は、選択された前記組み合わせに対応する前記近似係数に基づき、前記音圧を算出するものである。

【発明の効果】

【0009】

以上のように、本発明によれば、計算精度が最も高い次数および距離間隔の組み合わせが選択され、選択された組み合わせにおける次数および距離間隔に基づき近似係数が算出される。そして、算出された近似係数に基づき音圧が算出される。そのため、音圧計算に対して許容できる許容時間を満足しながら、俯仰角方向に対する最大の角度範囲で音波伝搬特性を精度よく算出することができる。

【図面の簡単な説明】

【0010】

【図1】実施の形態１に係る音波伝搬特性算出装置の構成の一例を示すブロック図である。

【図2】図１のパラメータ算出装置の構成の一例を示すブロック図である。

【図3】図２のパラメータ算出装置の構成の一例を示すハードウェア構成図である。

【図4】図２のパラメータ算出装置の構成の他の例を示すハードウェア構成図である。

【図5】実施の形態１に係る制御装置による処理の流れの一例を示すフローチャートである。

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

【0011】

以下、本発明の実施の形態について、図面を参照して説明する。本発明は、以下の実施の形態に限定されるものではなく、本発明の主旨を逸脱しない範囲で種々に変形することが可能である。また、本発明は、以下の各実施の形態に示す構成のうち、組合せ可能な構成のあらゆる組合せを含むものである。また、各図において、同一の符号を付したものは、同一のまたはこれに相当するものであり、これは明細書の全文において共通している。

【0012】

実施の形態１．
本実施の形態１に係る音波伝搬特性算出装置について説明する。本実施の形態１に係る音波伝搬特性算出装置は、ＰＥ法の特にＳｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法を用いて、俯仰角方向の角度範囲についての音波伝搬特性を算出するものである。

【0013】

［ＰＥ法およびＳｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法の概要］
まず、本実施の形態１に係る音波伝搬特性算出装置について説明する前に、音波伝搬特性を算出する際に用いられるＰＥ法およびＳｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法について説明する。ＰＥ法は、音源からの距離ｒおよび深度ｚに対する音圧ｐ（ｒ，ｚ）および波数ｋ（ｒ，ｚ）が満足する、式（１）に示すヘルムホルツ方程式から得られる式（２）の関係に基づき、音圧ｐを距離ｒについて逐次的に求めるものである。

【0014】

【数1】

【0015】

【数2】

【0016】

式（２）において、関数ψ（ｒ，ｚ）は、式（３）の関係を満足するものである。式（３）において、「Ｈ_０ ^（１）（ｋ_０ｒ）」は第１種の０次ハンケル関数を示す。「Δｒ」は距離方向に音圧を逐次的に計算する際の距離間隔を示す。「ｋ_０」は基準波数を示す。また、「ｑ」は、式（４）で表される微分作用素を示す。

【0017】

【数3】

【0018】

【数4】

【0019】

ここで、式（２）の右辺にある指数関数の項は、式（５）に示すように近似することができる。式（５）において、「ａ_ｊ，ｍ」および「ｂ_ｊ，ｍ」（ｊ＝１，２，・・・，ｍ）は近似関数の係数（以下、「近似係数」と称する）を示す。

【0020】

【数5】

【0021】

この近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍは、近似の次数ｍおよび距離間隔Δｒに依存して決定される。このような近似を用いて音圧を算出する方法は、特に「Ｓｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法」と称される。

【0022】

式（５）における「ｑ」は、本来、微分作用素を示すが、近似精度を評価するために変数とみなすと、変数ｑの値が「０」である場合、式（５）の近似式における両辺の値は共に「１」となるため、両辺の値の間に誤差は生じない。一方、「ｑ」が値「０」から離れるに従って、式（５）の両辺の値の差が大きくなるため、両辺の値の間の誤差が大きくなる。また、変数ｑは、伝搬する音波の俯仰角θとの間に、式（６）に示す関係がある（例えば、非特許文献２参照）。

【0023】

【数6】

【0024】

海洋中の音波は、水平方向に近い俯仰角から垂直方向に近い俯仰角の範囲に広がって伝搬するため、０°から±９０°までの角度範囲θを有している。そのため、式（６）において、俯仰角θの取り得る範囲は「－９０°≦θ≦９０°」となり、これにより、変数ｑの取り得る範囲は、「－１≦ｑ≦０」となる。したがって、俯仰角θが「０°」（ｑ＝０）に近いほど誤差が小さく、俯仰角θが「±９０°」（ｑ＝－１）に近いほど誤差が大きくなる。

【0025】

ここで、広い角度範囲で音圧計算の誤差を小さくするためには、式（５）の近似式を用いることによって生じる両辺の値の間の誤差を、変数ｑの広い範囲で小さくする必要がある。そのためには、近似の次数ｍを大きく、あるいは、距離間隔Δｒを小さくすればよい。

【0026】

しかしながら、音圧計算による音波伝搬特性の算出の際の計算量は、次数ｍに比例し、距離間隔Δｒに反比例するため、いずれの場合も計算量は増加してしまう。次数ｍを必要以上に大きくする、あるいは、距離間隔Δｒを必要以上に小さくすると、無駄に計算量が増加したり、次数ｍの不足または大きな距離間隔Δｒの設定によって俯仰角θに対する必要な角度範囲が達成できなかったりしてしまう。

【0027】

また、音波伝搬特性の計算は、際限なく時間をかけてよいわけではなく、それに費やすことが可能な計算時間の上限が存在することが多い。そのため、音波伝搬特性の計算に対して許容できる計算時間（以下、「許容時間」と適宜称する）で、最大の角度範囲を得られるように、次数ｍおよび距離間隔Δｒを適切に設定する必要がある。

【0028】

そこで、本実施の形態１に係る音波伝搬特性算出装置では、音波伝搬特性を算出する際の許容時間において最大の角度範囲を得るための次数ｍおよび距離間隔Δｒを適切かつ自動的に設定し、音波伝搬特性を精度よく算出することができるようにする。

【0029】

［音波伝搬特性算出装置１００の構成］
図１は、本実施の形態１に係る音波伝搬特性算出装置の構成の一例を示すブロック図である。図１に示すように、音波伝搬特性算出装置１００は、パラメータ算出装置１０、制御装置２０および音圧算出装置３０を含んで構成されている。パラメータ算出装置１０、制御装置２０および音圧算出装置３０のそれぞれは、ソフトウェアを実行することにより各種機能を実現するマイクロコンピュータ等の演算装置、もしくは各種機能に対応する回路デバイス等のハードウェアで構成されている。

【0030】

パラメータ算出装置１０は、制御装置２０から次数ｍ（ｍ＝１，２，・・・，Ｍ）および許容時間ｔ_ｍａｘを受け取り、受け取った次数ｍおよび許容時間ｔ_ｍａｘに基づき、近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ（ｊ＝１，２，・・・，ｍ）、角度範囲θならびに距離間隔Δｒをパラメータとして算出する。そして、パラメータ算出装置１０は、算出した近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ、角度範囲θならびに距離間隔Δｒを制御装置２０に供給する。ここで、「Ｍ」は制御装置２０から送られる次数ｍの最大値である。

【0031】

制御装置２０は、パラメータ算出装置１０に対して、パラメータ算出装置１０が各種のパラメータを算出する際に必要な次数ｍおよび許容時間ｔ_ｍａｘを送信する。また、制御装置２０は、次数ｍ毎にパラメータ算出装置１０から近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ、角度範囲θならびに距離間隔Δｒを受け取り、受け取った角度範囲θのうち最大のものに対応する近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを音圧算出装置３０に供給する。

【0032】

音圧算出装置３０は、制御装置２０から供給された近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍに基づき、音圧ｐを算出する。

【0033】

（パラメータ算出装置１０）
図２は、図１のパラメータ算出装置の構成の一例を示すブロック図である。図２に示すように、パラメータ算出装置１０は、距離間隔算出部１１、近似パラメータ算出部１２、真値算出部１３、位相誤差算出部１４および角度範囲算出部１５を備えている。

【0034】

距離間隔算出部１１は、制御装置２０から供給される次数ｍおよび許容時間ｔ_ｍａｘに基づき、距離間隔Δｒを算出する。距離間隔算出部１１は、算出した距離間隔Δｒを近似パラメータ算出部１２、真値算出部１３および制御装置２０に供給する。

【0035】

近似パラメータ算出部１２は、制御装置２０から供給される次数ｍと、距離間隔算出部１１から供給される距離間隔Δｒとに基づき、式（５）で表される近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを算出する。また、近似パラメータ算出部１２は、算出した近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍに基づき、式（５）の右辺に示される近似値ｆ＾を算出する。近似パラメータ算出部１２は、算出した近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを制御装置２０に供給し、算出した近似値ｆ＾を位相誤差算出部１４に供給する。

【0036】

真値算出部１３は、距離間隔算出部１１から供給される距離間隔Δｒに基づき、式（５）の左辺に示される真値ｆを算出する。真値算出部１３は、算出した真値ｆを位相誤差算出部１４に供給する。

【0037】

位相誤差算出部１４は、近似パラメータ算出部１２から供給される近似値ｆ＾と、真値算出部１３から供給される真値ｆとに基づき、位相誤差εを算出する。位相誤差算出部１４は、算出した位相誤差εを角度範囲算出部１５に供給する。

【0038】

角度範囲算出部１５は、位相誤差算出部１４から供給される位相誤差εに基づき、閾値ε_ｍａｘに対してε≦ε_ｍａｘとなる最小のｑを算出する。また、角度範囲算出部１５は、式（６）の関係によって角度範囲θを算出する。

【0039】

図３は、図２のパラメータ算出装置の構成の一例を示すハードウェア構成図である。パラメータ算出装置１０の各種機能がハードウェアで実行される場合、図２のパラメータ算出装置１０は、図３に示すように、処理回路４１で構成される。図２のパラメータ算出装置１０において、距離間隔算出部１１、近似パラメータ算出部１２、真値算出部１３、位相誤差算出部１４および角度範囲算出部１５の各機能は、処理回路４１により実現される。

【0040】

各機能がハードウェアで実行される場合、処理回路４１は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Application Specific Integrated Circuit）、ＦＰＧＡ（Field-Programmable Gate Array）、またはこれらを組み合わせたものが該当する。パラメータ算出装置１０は、距離間隔算出部１１、近似パラメータ算出部１２、真値算出部１３、位相誤差算出部１４および角度範囲算出部１５の各部の機能それぞれを処理回路４１で実現してもよいし、各部の機能を１つの処理回路４１で実現してもよい。

【0041】

図４は、図２のパラメータ算出装置の構成の他の例を示すハードウェア構成図である。パラメータ算出装置１０の各種機能がソフトウェアで実行される場合、図２のパラメータ算出装置１０は、図４に示すように、プロセッサ４２およびメモリ４３で構成される。パラメータ算出装置１０において、距離間隔算出部１１、近似パラメータ算出部１２、真値算出部１３、位相誤差算出部１４および角度範囲算出部１５の各機能は、プロセッサ４２およびメモリ４３により実現される。

【0042】

各機能がソフトウェアで実行される場合、パラメータ算出装置１０において、距離間隔算出部１１、近似パラメータ算出部１２、真値算出部１３、位相誤差算出部１４および角度範囲算出部１５の機能は、ソフトウェア、ファームウェア、またはソフトウェアとファームウェアとの組み合わせにより実現される。ソフトウェアおよびファームウェアは、プログラムとして記述され、メモリ４３に格納される。プロセッサ４２は、メモリ４３に記憶されたプログラムを読み出して実行することにより、各部の機能を実現する。

【0043】

メモリ４３として、例えば、ＲＡＭ（Random Access Memory）、ＲＯＭ（Read Only Memory）、フラッシュメモリ、ＥＰＲＯＭ（Erasable and Programmable ROM）およびＥＥＰＲＯＭ（Electrically Erasable and Programmable ROM）等の不揮発性または揮発性の半導体メモリ等が用いられる。また、メモリ４３として、例えば、磁気ディスク、フレキシブルディスク、光ディスク、ＣＤ（Compact Disc）、ＭＤ（Mini Disc）およびＤＶＤ（Digital Versatile Disc）等の着脱可能な記録媒体が用いられてもよい。

【0044】

このように、パラメータ算出装置１０は、ハードウェア、ソフトウェア、ファームウェアまたはこれらの組み合わせによって、上述した各機能を実現することができる。なお、制御装置２０および音圧算出装置３０のそれぞれについても、パラメータ算出装置１０と同様に、ハードウェア、ソフトウェア、ファームウェアまたはこれらの組み合わせによって、上述した各機能を実現することができる。

【0045】

［音波伝搬特性算出装置１００の動作］
次に、本実施の形態１に係る音波伝搬特性算出装置１００の動作について説明する。ここでは、パラメータ算出装置１０による近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ、角度範囲θならびに距離間隔Δｒの算出方法と、制御装置２０による最大角度範囲θ_ｍａｘ、ならびに近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍの決定方法と、音圧算出装置３０による音圧ｐの算出方法とについて説明する。

【0046】

（近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ、角度範囲θならびに距離間隔Δｒの算出）
まず、パラメータ算出装置１０が制御装置２０から次数ｍおよび許容時間ｔ_ｍａｘを受け取ると、距離間隔算出部１１は、次数ｍおよび許容時間ｔ_ｍａｘに基づき距離間隔Δｒを算出する。Ｓｐｌｉｔ－ｓｔｅｐ Ｐａｄｅ法による音圧計算の際の計算時間ｔは、距離方向のグリッド数Ｋ（＝ｒ_ｍａｘ／Δｒ）、深度方向のグリッド数Ｌ（＝ｚ_ｍａｘ／Δｚ）および次数ｍに比例する。ここで、「ｒ_ｍａｘ」は計算最大距離を示し、「ｚ_ｍａｘ」は計算最大深度を示す。「Δｚ」は計算の深度間隔を示す。したがって、予め決められた基準条件での計算時間ｔ^（０）に対して、計算時間ｔは、式（７）に基づき算出される。

【0047】

【数7】

【0048】

式（７）において、「Ｋ^（０）」は基準条件での距離方向のグリッド数を示し、「Ｋ^（０）＝ｒ^（０） _ｍａｘ／Δｒ^（０）」である。「Ｌ^（０）」は基準条件での深度方向のグリッド数を示し、「Ｌ^（０）＝ｚ^（０） _ｍａｘ／Δｚ^（０）」である。「ｍ^（０）」は基準条件での次数を示す。また、「ｒ^（０） _ｍａｘ」は基準条件での計算最大距離を示し、「ｒ^（０）」は基準条件での距離間隔を示す。「ｚ^（０） _ｍａｘ」は基準条件での計算最大深度を示し、「ｚ^（０）」は基準条件での深度間隔を示す。

【0049】

さらに、基準条件での計算時間ｔ^（０）は、これらの基準条件に基づいて音圧算出装置３０として用いられる計算機上で音圧ｐを算出した際の計算時間である。このことから、式（７）は、式（８）に示すように変形することができる。

【0050】

【数8】

【0051】

式（８）において、計算時間ｔを許容時間ｔ_ｍａｘとした場合、式（８）は、式（９）に示すように、許容時間ｔ_ｍａｘに対する距離間隔Δｒを算出するための関係式に変形することができる。

【0052】

【数9】

【0053】

このようにして、距離間隔算出部１１は、次数ｍおよび許容時間ｔ_ｍａｘに基づき距離間隔Δｒを算出する。式（９）においては、次数ｍが大きいほど、それに比例して距離間隔Δｒも大きくなり、計算時間が一定に保たれる関係となっている。なお、上述した基準条件である各種の値「ｒ^（０） _ｍａｘ」、「ｚ^（０） _ｍａｘ」、「Δｒ^（０）」、「Δｚ^（０）」および「ｍ^（０）」は、基準条件での計算時間ｔ^（０）を得るための暫定的な条件である。この場合、異なる基準条件を使ったとしても、それによって計測した計算時間をｔ^（０）とすれば、式（９）によって同じ距離間隔Δｒを得ることができる。

【0054】

次に、近似パラメータ算出部１２は、制御装置２０から受け取った次数ｍと、距離間隔算出部１１で算出された距離間隔Δｒとに基づき、式（５）で表される近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ（ｊ＝１，２，・・・，ｍ）を算出する。また、近似パラメータ算出部１２は、算出した近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを用いて、変数ｑ_ｎに対して、式（５）の右辺による近似値ｆ＾（ｑ_ｎ）を、式（１０）に基づき算出する。ここで、変数ｑ_ｎは、「ｑ_０＝－１」、「ｑ_１＝ｑ_０＋Δｑ」、・・・、「ｑ_Ｎ＝ｑ_０＋ＮΔｑ」となる値である。「Ｎ」は、予め決めた自然数を示し、「Δｑ＝１／Ｎ」である。

【0055】

【数10】

【0056】

また、真値算出部１３は、距離間隔算出部１１で算出された距離間隔Δｒを用いて、変数ｑ_ｎに対して、式（５）の左辺による真値ｆ（ｑ_ｎ）を、式（１１）に基づき算出する。

【0057】

【数11】

【0058】

次に、位相誤差算出部１４は、近似パラメータ算出部１２で算出された近似値ｆ＾（ｑ_ｎ）（ｎ＝０，１，・・・，Ｎ）と、真値算出部１３で算出された－１≦ｑ≦０に対する真値ｆ（ｑ_ｎ）（ｎ＝０，１，・・・，Ｎ）とに基づき、式（１２）を用いて位相誤差ε（ｑ_ｎ）を算出する。式（１２）において、「φ［・］」は、複素数ζがその振幅Ａと位相Ｐによってζ＝Ａｅ^ｉＰと表されるときに、φ［ζ］＝Ｐとなるような関数である。

【0059】

【数12】

【0060】

角度範囲算出部１５は、位相誤差算出部１４で算出された位相誤差ε（ｑ_ｎ）（ｎ＝０，１，・・・，Ｎ）と予め設定された閾値ε_ｍａｘとを比較し、「ε（ｑ_ｎ）≦ε_ｍａｘ」となるときの最小の変数ｑ_ｎを算出する。そして、角度範囲算出部１５は、算出した変数ｑ_ｎに基づき、式（６）の関係によって角度範囲θを算出する。すなわち、角度範囲θは、式（６）の関係から得られる式（１３）によって算出される。

【0061】

【数13】

【0062】

このようにして、パラメータ算出装置１０は、最大の角度範囲θと、最大の角度範囲θで音圧ｐを適切に算出するための近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍとを制御装置２０で決定するための近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ、角度範囲θならびに距離間隔Δｒを算出する。なお、式（５）の近似式に含まれる「ｑ」は、本来は微分作用素だが、式（５）の近似誤差を定量的に表すために、パラメータ算出装置１０では、「ｑ」を「－１≦ｑ≦０」の範囲の値をとる変数として扱っている。

【0063】

（最大角度範囲θ_ｍａｘおよび近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍの決定）
図５は、本実施の形態１に係る制御装置による処理の流れの一例を示すフローチャートである。まず、制御装置２０は、ステップＳ１において、角度範囲θの最大値である最大角度範囲θ_ｍａｘを初期化し、値を「０」とする。また、ステップＳ２において、制御装置２０は、近似の次数ｍを初期化し、値を「０」とする。

【0064】

ステップＳ３において、制御装置２０は、次数ｍの値を「１」から順に１ずつインクリメントする。そして、制御装置２０は、設定した次数ｍと、予め設定された許容時間ｔ_ｍａｘをパラメータ算出装置１０に供給する。その後、供給された次数ｍおよび許容時間ｔ_ｍａｘに基づき、パラメータ算出装置１０によって各種のパラメータが算出されると、制御装置２０は、パラメータ算出装置１０から近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ、角度範囲θならびに距離間隔Δｒを受け取る。

【0065】

次に、ステップＳ４において、制御装置２０は、パラメータ算出装置１０から受け取った距離間隔Δｒと、予め設定された距離間隔Δｒの最大値である最大距離間隔Δｒ_ｍａｘとを比較する。比較の結果、距離間隔Δｒが最大距離間隔Δｒ_ｍａｘ以下である場合（ステップＳ４：Ｙｅｓ）には、処理がステップＳ５に移行する。

【0066】

一方、距離間隔Δｒが最大距離間隔Δｒ_ｍａｘよりも大きい場合（ステップＳ４：Ｎｏ）、制御装置２０は、保持している近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを音圧算出装置３０に供給する。そして、一連の処理が終了する。

【0067】

ステップＳ５において、制御装置２０は、パラメータ算出装置１０から受け取った角度範囲θと、予め設定された最大角度範囲θ_ｍａｘとを比較する。比較の結果、角度範囲θが最大角度範囲θ_ｍａｘよりも大きい場合（ステップＳ５：Ｙｅｓ）、制御装置２０は、ステップＳ６において、現在の角度範囲θを最大角度範囲θ_ｍａｘとして更新する。そして、ステップＳ７において、制御装置２０は、パラメータ算出装置１０から受け取った近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを保持する。その後、処理がステップＳ３に戻る。

【0068】

一方、ステップＳ５において、角度範囲θが最大角度範囲θ_ｍａｘ以下である場合（ステップＳ５：Ｎｏ）には、処理がステップＳ３に戻る。

【0069】

ここで、上述した式（９）で示されるように、パラメータ算出装置１０から送られる距離間隔Δｒは、次数ｍに比例して大きくなる。そのため、次数ｍが１ずつ大きくなっていけば、いずれ「Δｒ＞Δｒ_ｍａｘ」となり、ステップＳ３～ステップＳ７の処理のループは、必ず有限の回数で終了する。つまり、１ずつ大きくなっていく次数ｍの最大値Ｍは、式（９）の左辺をΔｒからΔｒ_ｍａｘに変えて、次数ｍについて解いた式によって得られる。ただし、次数ｍは整数でなければならないため、最大値Ｍは、式（１４）に基づき決定される。式（１４）において、ｆｌｏｏｒ［・］は小数点以下を切り捨てる関数である。

【0070】

【数14】

【0071】

図５のステップＳ３～ステップＳ７によるループ中の処理では、まず、次数ｍに１を加算し、次数ｍを１だけ大きくする（ステップＳ３）。そして、制御装置２０は、次数ｍおよび許容時間ｔ_ｍａｘをパラメータ算出装置１０に供給する。

【0072】

次に、制御装置２０は、パラメータ算出装置１０から供給される近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍ、角度範囲θおよび距離間隔Δｒを受け取る。制御装置２０は、距離間隔Δｒと最大距離間隔Δｒ_ｍａｘとを比較し（ステップＳ４）、「Δｒ＞Δｒ_ｍａｘ」であれば（ステップＳ４：Ｎｏ）、ループを抜ける。

【0073】

一方、「Δｒ≦Δｒ_ｍａｘ」であれば（ステップＳ４：Ｙｅｓ）、制御装置２０は、パラメータ算出装置１０から供給された角度範囲θと最大角度範囲θ_ｍａｘとを比較する（ステップＳ５）。もし、「θ＞θ_ｍａｘ」であれば（ステップＳ５：Ｙｅｓ）、制御装置２０は、最大角度範囲θ_ｍａｘを「θ_ｍａｘ＝θ」とし、そのときの近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを保持する。そして、処理がループの先頭（ステップＳ３）に戻る。また、「θ≦θ_ｍａｘ」であれば（ステップＳ５：Ｎｏ）、そのままループの先頭に戻る。

【0074】

以上のようにして、距離間隔Δｒが「Δｒ＞Δｒ_ｍａｘ」となるまでステップＳ３～ステップＳ７のループが継続され、ループを抜けた時点で保持されている近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍが音圧算出装置３０に供給される。

【0075】

（音圧ｐの算出）
音圧算出装置３０は、上述したようにして制御装置２０で決定された近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを受け取ると、式（２）の指数項を式（５）で近似した式（１５）と、式（３）の関係とに基づき、音圧ｐ（ｒ，ｚ）を算出する。

【0076】

【数15】

【0077】

このように、音波伝搬特性算出装置１００では、許容時間ｔ_ｍａｘを満足する近似の次数ｍと、距離間隔Δｒとの組み合わせのうち、最大の角度範囲θ_ｍａｘを有する場合の組み合わせ、つまり精度の最もよくなる組み合わせが自動的に設定される。そして、設定された次数ｍと距離間隔Δｒとに基づき算出された近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍを用いて音圧ｐが算出される。そのため、音圧計算に対する許容時間を満足しながら、俯仰角方向に対する最大の角度範囲で音波伝搬特性を精度よく算出することができる。

【0078】

なお、上述の例では、変数ｑの値を「－１≦ｑ≦０」の範囲で等間隔となるように「ｑ_ｎ＝ｑ_０＋ｎΔｑ」（ｎ＝０，１，２，・・・，Ｎ）と設定したが、これはこの例に限られない。例えば、角度範囲θの値を「０°≦θ≦９０°」の範囲で「θ_０＝０」、「θ_１＝Δθ」、・・・、「θ_Ｎ＝ＮΔθ＝９０°」と設定し、式（６）に基づき、変数ｑ_ｎを「ｑ_ｎ＝－ｓｉｎ^２θ_ｎ」（ｎ＝０，１，２，・・・，Ｎ）として導出してもよい。この場合、「Δθ」は、「Δθ＝９０／Ｎ」から得られる値である。

【0079】

また、変数ｑ_ｎ（ｎ＝０，１，２，・・・，Ｎ）は、上述した方法以外の方法によって導出されてもよい。例えば、式（１２）に基づき得られる位相誤差ε（ｑ_ｎ）は、「－１≦ｑ≦０」で単調変化するため、単調変化する関数の根を見つける一般的な数値計算手法、例えば二分法を用いてもよい。この場合、式（１６）の根を求めることにより、「ε（ｑ）＝ε_ｍａｘ」となる変数ｑを得ることができる。そして、角度範囲θは、「θ＝ｓｉｎ^－１√－ｑ」によって算出することができる。

【0080】

【数16】

【0081】

以上のように、本実施の形態１に係る音波伝搬特性算出装置１００では、計算精度が最も高い次数ｍおよび距離間隔Δｒの組み合わせが選択され、選択された組み合わせにおける次数ｍおよび距離間隔Δｒに基づき近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍが算出される。そして、算出された近似係数ａ_ｊ，ｍおよびｂ_ｊ，ｍに基づき音圧ｐが算出される。そのため、音圧計算に対して許容できる許容時間ｔ_ｍａｘを満足しながら、俯仰角方向に対する最大の角度範囲θ_ｍａｘで音波伝搬特性を精度よく算出することができる。

【0082】

音波伝搬特性算出装置１００では、近似式における真値ｆ（ｑ_ｎ）と近似値ｆ＾（ｑ_ｎ）との位相誤差ε（ｑ_ｎ）に基づき、計算精度を評価する。このとき、位相誤差ε（ｑ_ｎ）が予め設定された閾値ε_ｍａｘ以下であるか否かに基づき、計算精度を評価してもよい。真値ｆ（ｑ_ｎ）と近似値ｆ＾（ｑ_ｎ）との位相誤差ε（ｑ_ｎ）が小さく、閾値ε_ｍａｘ以下であれば、誤差が小さく精度が高いと判断することができる。

【符号の説明】

【0083】

１０ パラメータ算出装置、１１ 距離間隔算出部、１２ 近似パラメータ算出部、１３ 真値算出部、１４ 位相誤差算出部、１５ 角度範囲算出部、２０ 制御装置、３０ 音圧算出装置、４１ 処理回路、４２ プロセッサ、４３ メモリ、１００ 音波伝搬特性算出装置。

【図1】

【図2】

【図3】

【図4】

【図5】