JP6029124B1 - 音源位置推定装置、方法、及びプログラム - Google Patents

Abstract

【課題】音源の位置を精度よく推定する音源位置推定装置、方法及びプログラムを提供すること。【解決手段】音源位置推定装置３０は、音源から受信した複数の到来波の尤度を算出する尤度算出部６２と、尤度算出部６２により算出された各到来波の各位置の尤度分布のシグマポイントを算定するシグマポイント算定部７０と、シグマポイント算定部７０により算定された各シグマポイントを用いて、組み合わせ可能な２つの到来波の分布間距離をそれぞれ求めた上でこれらの総和を算出すると共に、該総和が最小となる時刻を特定する収束時刻特定部８１と、収束時刻特定部８１により特定された時刻に基づいて音源の位置を推定する音源位置推定部８３と、を有する。【選択図】図３

Description

本発明は、音源から受信した音波をもとに音源位置を推定する音源位置推定装置、方法、及びプログラムに関する。

水中での音源−受波器間における音波伝搬では、音源から直接到来する直接波及び海面や海底で反射して到来する各マルチパス波の相互間において到来時刻にずれ（到来時間差）が生じる。従来から、音源の位置を推定する場合には、到来時間差を反転させた上で各マルチパス波の逆伝搬経路を推定し、推定した逆伝搬経路上の位置に各マルチパス波が存在する尤度を求めた上で、尤度の分布に基づいて音源位置を推定する、という手法が採られている。

例えば、特許文献１に開示された音源位置推定装置は、一連の到来波の到来順序とは逆の順序のパルス波列を、各到来波の到来ふ仰角と同じ角度で受波器群から送出した場合の伝搬経路を推定し、送出後の各時刻に各パルス波が各位置に存在する尤度を推定し、すべてのパルス波についての尤度の合計を求め、複数のパルス波についての尤度分布の間の距離（分布間距離）が最小となる時刻（収束時刻）において、尤度の合計が最大となる位置を音源位置と推定するものである。特許文献１の技術において、分布間距離の算出に用いるマハラノビス距離は、２つのマルチパス波の尤度分布の中心間距離を、各々の尤度分布の共分散行列によって規格化した平面で測るものである。

特開２０１３−１７０８１７号公報

しかしながら、特許文献１の音源位置推定装置は、マハラノビス距離の考え方を用いて２つのマルチパス波の分布間距離を求めるという手法を採っているため、分布間距離の算出精度が各マルチパス波の尤度分布の広がりの状態に左右され、収束時刻及び音源位置を精度よく推定できないという課題がある。この課題に鑑みて、従来から、音源の位置を精度よく推定する音源位置推定装置、方法、及びプログラムの提供が望まれていた。

本発明に係る音源位置推定装置は、音源から受信した複数のマルチパス波に基づいて、各マルチパス波の到来時間及び到来ふ仰角の情報を算出する受波情報算出部と、各マルチパス波の到来時間を反転する到来時間反転部と、到来時間反転部により反転された各到来時間と受波情報算出部により算出された到来ふ仰角の情報とに基づいて、各マルチパス波の各位置の尤度分布を算出する尤度算出部と、尤度算出部により算出された各マルチパス波の各位置の尤度分布について複数のシグマポイントを算定するシグマポイント算定部と、各マルチパス波の内の２つのマルチパス波を複数組み合わせた上で、シグマポイント算定部により算定された各シグマポイントを用いて各組み合わせの尤度分布間の距離を、一方のマルチパス波におけるシグマポイントと他方のマルチパス波におけるシグマポイントとの複数の組み合わせについて計算する分布間距離計算部と、分布間距離計算部により計算された各尤度分布間の距離の総和を算出すると共に、該総和が最小となる時刻を特定する収束時刻特定部と、収束時刻特定部により特定された時刻に基づいて音源の位置を推定する音源位置推定部と、を有するものである。

また、本発明に係る音源位置推定方法は、音源から受信した複数のマルチパス波に基づいて、各マルチパス波の到来時間及び到来ふ仰角の情報を算出し、各マルチパス波の到来時間を反転し、反転された各到来時間と算出された到来ふ仰角の情報とに基づいて、各マルチパス波の各位置の尤度分布を算出し、算出された各マルチパス波の各位置の尤度分布について複数のシグマポイントを算定し、各マルチパス波の内の２つのマルチパス波を複数組み合わせた上で、算定された各シグマポイントを用いて各組み合わせの尤度分布間の距離を、一方のマルチパス波におけるシグマポイントと他方のマルチパス波におけるシグマポイントとの複数の組み合わせについて計算し、計算された各尤度分布間の距離の総和を算出すると共に、該総和が最小となる時刻を特定し、特定された時刻に基づいて音源の位置を推定する。

さらに、本発明に係る音源位置推定プログラムは、音源から受信した複数のマルチパス波に基づいて、各マルチパス波の到来時間及び到来ふ仰角の情報を算出する受波情報算出手段、各マルチパス波の到来時間を反転する到来時間反転手段、到来時間反転手段により反転された各到来時間と受波情報算出手段により算出された到来ふ仰角の情報とに基づいて、各マルチパス波の各位置の尤度分布を算出する尤度算出手段、尤度算出手段により算出された各マルチパス波の各位置の尤度分布について複数のシグマポイントを算定するシグマポイント算定手段、各マルチパス波の内の２つのマルチパス波を複数組み合わせた上で、シグマポイント算定手段により算定された各シグマポイントを用いて各組み合わせの尤度分布間の距離を、一方のマルチパス波におけるシグマポイントと他方のマルチパス波におけるシグマポイントとの複数の組み合わせについて計算する分布間距離計算手段、分布間距離計算手段により計算された各尤度分布間の距離の総和を算出すると共に、該総和が最小となる時刻を特定する収束時刻特定手段、及び収束時刻特定手段により特定された時刻に基づいて音源の位置を推定する音源位置推定手段、としてコンピュータを機能させる。

本発明では、尤度算出部の出力側にシグマポイント算定部を設け、シグマポイント算定部で算定した各シグマポイントに基づいて収束時刻を特定するという構成を採ったため、収束時刻に係る誤判定を低減させることができ、音源の位置を精度よく推定することができる。

本発明の実施形態における音源位置推定装置の概略図である。 図１の受波器群の構成を例示する概略図である。 図１の音源位置推定装置の構成を示すブロック図である。 図３の到来時間反転部による各マルチパス波の到来時間差の順序反転を示す説明図である。 図３の尤度算出部の処理に係る音線の経路方向とこれに垂直な方向を表す模式図である。 図３の収束時刻特定部の具体的構成を示すブロック図である。 図３の収束時刻特定部による判定に係る２つのシグマポイント間の距離の加重平均に基づく分布間距離のイメージ図である。 図３等の音源位置推定装置の動作を示すフローチャートである。 図３のシグマポイント算定部にて算定したシグマポイントの間の距離が海底反射時に連続性を保つ状況を示す模式図である。 図３の収束時刻特定部によって合算した分布間距離の総和が、逆伝搬時間の経過に伴って変化する様子を示すグラフである。 本発明の実施形態において、誤検出を含む場合のシグマポイント間距離による分布間距離の総和の変動を示すグラフである。 マハラノビス距離による分布間距離のイメージを示す説明図である。 関連技術において、音波が海底反射する時にマハラノビス距離が不連続となる状況を示す説明図である。 関連技術において、マハラノビス距離による分布間距離の総和が、逆伝搬時間の経過に伴って変化する様子を示すグラフである。 関連技術において、誤検出を含む場合のマハラノビス距離による分布間距離の総和の変動を示すグラフである。

（実施形態の構成）
本発明における音源位置推定装置の一実施形態を図１〜図１１に基づいて説明する。図１に示すように、音源位置推定システム１０は、音源９０から放射されたパルス状音波である直接波及びマルチパス波（到来波）を受信すると共に、到来波の到来時間差に基づいて音源の位置を推定するものであって、二以上の受波器から成り各々において到来波を受信する受波器群２０と、これに接続された音源位置推定装置３０と、を有している。

図１では便宜上、海面から海底までの距離（水深）Ｚｗが一定（海底がフラット）である場合を想定し、海面から受波器群２０、音源９０までの垂直距離（深度）をそれぞれＺｒ、Ｚｓとする。音源９０と音源位置推定装置３０との間の水平距離はＲｗとする。

また、図１では、音源９０から到来する直接波及びマルチパス波（到来波）の伝搬経路（音線）の内の５つ（１〜５）を例示し、直接波１を実線、マルチパス波２〜５を破線でそれぞれ示す。実際には、水温・塩分濃度等の不均一性により、音速が一定とはならないため、伝搬経路は直線的にはならないが、ここでは便宜的に直線で示している。

一般に、音源と受波器群との間に他の音源等がない場合の音波伝搬は、自由空間の波動方程式で表され、自由空間の波動方程式の解は、相反性を有するＧｒｅｅｎ関数で表される。この相反性が成り立つことにより、音源から出て受波器群に到達した直接波及びマルチパス波の順序を逆にして（時間反転して）受波器群から逆方向に伝搬（逆伝搬）させると、受波器群から出た仮想音波は音源の位置に理論上収束する。

この音波の逆伝搬による音源位置への収束は、波動方程式を近似した音線計算においても成立する。すなわち、受信時の到来時間・到来角を観測し、到来波の順序を逆にして、到来角と同じ角度（同じふ仰角及び方角）で逆向きに受波器群からパルス波を送出したと仮定する。そして、送出された各パルス波の各時点における到達位置が音線計算により推定された場合、各パルス波の到達位置は時間とともに変化するが、すべてのパルス波は同時に音源位置に到達（収束）することになる。本実施形態における音源位置推定装置３０も、相反性に基づく到来波の時間反転を利用した逆伝搬処理によって音源位置を推定するように構成されている。

音源位置推定装置においては、マルチパス波の到来ふ仰角を判別するために、受波器群を二以上の受波器によって構成する。本実施形態における受波器群２０は、図２に例示するように、一定の間隔Ｄｈを隔てて配置された第１受波器２０Ａと第２受波器２０Ｂとを有する。第１受波器２０Ａ及び第２受波器２０Ｂは、到来波（直接波及びマルチパス波）のふ仰角を求めるために、同じ鉛直線上の異なる位置に設けられている。もっとも、受波器群２０は、三以上の受波器を有していてもよい。

図３に示すように、音源位置推定装置３０は、受波器群２０から入力される各マルチパス波の時間波形の情報等を解析して各マルチパス波の受波解析情報を算出する受波情報算出部４１と、受波情報算出部４１から入力されたマルチパス波の到来時間差を時間軸で反転する到来時間反転部（到来時間差反転部）４２と、入力条件に従って音線計算等を実施する伝搬計算処理部６１と、伝搬計算処理部６１が読み出して使用する各種の情報が格納された情報格納部５０と、を有している。

また、音源位置推定装置３０は、伝搬計算処理部６１から入力された各種計算の結果を用いて、マルチパス波ごとの尤度を算出する尤度算出部６２と、この尤度算出部６２から入力された各尤度の分布のシグマポイントを算定するシグマポイント算定部７０と、このシグマポイント算定部７０において算定した各シグマポイントを用いて、組み合わせ可能な２つの到来波の分布間距離をそれぞれ求めた上でこれらの総和を算出すると共に、該総和が最小となる時刻（収束時刻）を特定（判定）する収束時刻特定部８１と、を有している。

さらに、音源位置推定装置３０は、尤度算出部６２から取得した各尤度を統合（合算）する統合尤度計算部８２と、この統合尤度計算部８２から上記統合の結果である統合尤度を取得すると共に、収束時刻特定部８１において特定した時刻における全尤度の分布に基づいて音源９０の位置を推定する音源位置推定部８３と、を有している。

なお、図３のような音源位置推定装置３０の構成は、補助記憶装置に読み込まれた音源位置推定プログラムをマイコンもしくはコンピュータ（たとえばパーソナルコンピュータ等）上で実行することにより実現される。音源位置推定プログラムは、ＣＤ−ＲＯＭ等の情報記憶媒体に記憶され、もしくはインターネット等のネットワークを介して配布され、コンピュータにインストールされることになる。

受波情報算出部４１は、第１受波器２０Ａ及び第２受波器２０Ｂから取得した時間波形の情報に基づいて、マルチパス波の到来時間差時系列（到来時間差）、マルチパス波の到来ふ仰角（到来ふ仰角）、到来時間差誤差標準偏差、及び到来ふ仰角誤差標準偏差を、上記受波解析情報として算出するものである。

受波情報算出部４１は、各マルチパス波の入力を複数回繰り返して各回の到来時間差、到来ふ仰角の平均値をそれぞれ到来時間差、到来ふ仰角として到来時間反転部４２に出力する。このようにすれば、到来波の検出漏れや到来波の誤検出による測定誤差が生じた場合にも、必要十分な回数分のデータをもとに、信頼性の高い到来時間差及び到来ふ仰角を算出することができる。

到来時間差誤差標準偏差は、各マルチパス波について観測された到来時間差のばらつきの程度を表す指標であり、到来ふ仰角誤差標準偏差は、各マルチパス波について観測された到来ふ仰角のばらつきの程度を表す指標である。

ところで、上記の通り、第１受波器２０Ａと第２受波器２０Ｂとは、間隔Ｄｈを隔てて配置されているため、同じ到来波の受信時刻が異なることとなる。受波情報算出部４１は、この受信時刻の差異の他、第１受波器２０Ａと第２受波器２０Ｂとの間隔Ｄｈ、及び受波器群２０の位置における音速に基づいて到来ふ仰角を算出する。また、到来時間差の算出に際して受波情報算出部４１は、２つの受波器２０Ａ、２０Ｂにおける同じマルチパス波の受信時刻の中間値を受波器群２０における受信時刻とする。

到来時間差及び到来ふ仰角の算出には、例えば 、２つの受波器間に到来した時系列信号の相関処理といった一般的な手法を用いることができる。また、到来時間差又は到来ふ仰角を観測する方法は種々存するが、所望の情報を観測できるものであれば、何れの方法も適用可能である。例えば、２つの受波器間の情報を用いる方法の他、受波器アレイを構成して整相処理を行い、整相処理の結果から所望の情報を得るという手法を適用するようにしてもよい。

受波情報算出部４１は、上記のように算出した一連の到来波についての到来時間差及び到来ふ仰角の情報（受波解析情報）を到来時間反転部４２に出力する。到来時間反転部（到来時間差反転器）４２は、受波情報算出部４１から出力された到来ふ仰角、到来時間差誤差標準偏差、到来ふ仰角誤差標準偏差はそのままに、到来時間差を反転して伝搬計算処理部６１に出力する。

ここで、受波器群２０にｎ個のマルチパス波が到来したことを想定し（ｎは２以上の自然数）、これら各々の到来時間差（第１番目の到来波の受信時刻と第ｉ番目の到来波の受信時刻の差）をτ_ｉ（ｉ＝１、２、…、ｎ）とすると、各マルチパス波の到来時間差の順序を反転したτ’_ｉ（反転時間差）は、下記式（１）によって求めることができる（τ_ｉ＋１＞τ_ｉ）。

図４に例示するように、式（１）を用いた処理を経て、ｎ番目のマルチパス波の到来時間差τ_ｎは基準の０（反転時間差τ’_ｎ）に、（ｎ−１）番目のマルチパス波の到来時間差τ_ｎ−１は反転時間差τ’_ｎ−１に、（ｎ−２）番目のマルチパス波の到来時間差τ_ｎ−２は反転時間差τ’_ｎ−２になる、というように変換される。

すなわち、到来時間反転部４２は、上記式（１）に係る演算により、各マルチパス波の到来時間差を、順序を反転した上で内部メモリ等（図示せず）に保存すると共に、各マルチパス波の反転時間差を、到来ふ仰角、到来時間差誤差標準偏差、及び到来ふ仰角誤差標準偏差と共に伝搬計算処理部６１へ向けて出力する。

伝搬計算処理部６１は、第１〜第ｎの伝搬計算器６１−１〜６１−ｎを有し（ｎは２以上の自然数）、入力条件に従って音線計算を行うものである。すなわち、各マルチパス波に対応する伝搬計算器６１−１〜６１−ｎが、マルチパス波ごとに、到来ふ仰角の情報を用いて音線計算（伝搬経路の推定）を実行する。

より具体的には、深度方向をｚ、受波器群２０から音源９０に向かう方向をｒとした場合に、伝搬計算処理部６１は、媒介変数ｓを用いた下記式（２ａ）〜（２ｄ）により音線計算（伝搬経路の推定）を実行する。

各式（２ａ）〜（２ｄ）において、関数ζと関数ξは、送出ふ仰角θに基づいてそれぞれζ＝ｓｉｎθ／Ｃ、ξ＝ｃｏｓθ／Ｃにより求められる。ここで、Ｃは、受波器群の位置における音速である。

このように、伝搬計算処理部６１は、各マルチパス波毎に、各τ'_ｉについて、受波器群２０から音源９０に向かう音線計算（伝搬経路の推定）を式（２ａ）〜（２ｄ）によって実現する。また、伝搬計算処理部６１は、音線計算等を行うに際して、情報格納部５０から音速Ｃ等の情報を取得する。

情報格納部５０には、音速プロファイル５０Ａと、海底地形を示す情報である海底地形情報５０Ｂと、各受波器２０Ａ、２０Ｂの位置に係る情報である受波器位置情報５０Ｃとが予め格納されている。

本実施形態では、実測処理部５５が、温度分布及び塩分濃度を実測し、この実測値を用いて、経験式により音速プロファイル５０Ａを算出して格納するという好適な構成を採った。他に、音速を実測した上でこれを音速プロファイル５０Ａとして格納する機能構成を、実測処理部５５として設けてもよい。また、本実施形態では、実測処理部５５を音源位置推定装置３０の内部に設けたが、これと同様に機能する機器等を該装置の外部に設けるようにしてもよい。

上記のように、現場で実測した音速又は温度分布、塩分濃度から経験式により算出した音速プロファイルを音速プロファイル５０Ａとして格納することが望ましいが、かかる好適な構成を採らない場合（各実測値が得られない場合）でも、一般に市販されているデータベースの値や、通常よく用いられる１５００ｍ／ｓ等の一定の値を音速プロファイル５０Ａとして用いることができる。

海底地形情報５０Ｂとしては、観測地もしくは一般に市販されているデータベースの値などを用い、受波器位置情報５０Ｃとしては、設置した各受波器の深度を予め入力しておく。

また、伝搬計算処理部６１は、伝搬経路上の逆伝搬時間ｔ_i（s）を下記式（３）により計算するものである。

式（３）は、媒介変数ｓに対する積分であるため、式（３）を用いることにより、音線の経路に沿った形で音速の逆数を積分して逆伝搬時間を求めることができる。また、各式（２ａ）〜（２ｄ）を用いることで、音波伝搬の開始点（受波器群２０）からの経路長を求めることができる。

すなわち、伝搬計算処理部６１は、上記式（２）及び（３）を用いた計算により、予め設定された時間間隔の各マルチパス波の到来位置（伝搬経路）、受波器群２０からの経路長を求め、求めた伝搬経路及び経路長を、各受波器２０Ａ、２０Ｂで観測した各到来波の到来時間差誤差標準偏差、及び到来ふ仰角誤差標準偏差と共に尤度算出部６２へと出力する。

ところで、逆伝搬時間ｔにおいて、第ｉ番目のマルチパス波が位置（ｒ_ｉ（ｔ）、ｚ_ｉ（ｔ））に存する場合の尤度は、下記式（４）で定義される。

ここで、ｐ（Ｐａｔｈ_ｉ）は、第ｉ番目のマルチパス波が生じる（第ｉ番目のマルチパス波が受波器群に到達する）事前確率である。通常は、受波器群側から音源の指向性などの情報は得られないと考えられるため、ラプラスの不十分理由の原理により各到来波についての事前確率は等しいとすべきである。また、条件付き確率は、具体的に確率密度関数が明確となる情報が得られれば、それに従って構成すべきである。確率密度関数が明確になる情報がない場合は、図５に示すように、音線の経路方向をｓ方向、音線経路に垂直な方向をη方向として、（ｒ_ｉ（ｔ）、ｚ_ｉ（ｔ））を原点とする局所座標の距離を（ｓ_ｉ（ｔ）、η_ｉ（ｔ））とする。

本実施形態において尤度算出部６２は、第ｉ番のマルチパス波の経路長をｌ_ｉ（ｔ）としたときに、下記式（５）を用いて、各マルチパス波の尤度Ｌ_ｉ（ｔ）（ｉ＝１〜ｎ）を算出する。ここで、σ_ｔは到来時間差誤差標準偏差、σ_θは到来ふ仰角誤差標準偏差である。

すなわち、式（５）を用いた処理が、第１〜第ｎの尤度計算器６２−１〜６２−ｎの各々において個別に行われ、各尤度計算器６２−１〜６２−ｎは、それぞれにおいて求めた尤度Ｌ_ｉ（ｔ）（ｉ＝１〜ｎ）をシグマポイント算定部７０及び統合尤度計算部８２に送信する。

シグマポイント算定部７０は、第１〜第ｎの尤度計算器６２−１〜６２−ｎの各々に個別に対応する（音源からの各到来波に個別に対応する）第１〜第ｎのシグマポイント計算器７０−１〜７０−ｎを有し、これら各シグマポイント計算器７０−１〜７０−ｎは、対応する第１〜第ｎの尤度計算器６２−１〜６２−ｎから取得した尤度Ｌ_１（ｔ）〜Ｌ_ｎ（ｔ）に基づいて各シグマポイントを計算する。

より具体的に、シグマポイント算定部７０は、第ｉ番目のマルチパス波の尤度を表現する分布の各シグマポイント（ｘ_ｉ ^（ｋ））を下記式（６）によって算定する。

このように、シグマポイント算定部７０は、各シグマポイント計算器７０−１〜７０−ｎによって、マルチパス波の尤度分布を表す平均と共分散行列から、各マルチパス波の尤度分布のシグマポイントを計算するように構成されている。本実施形態では、２次元平面での分布を想定しているため「Ｌ＝２」とし、スケーリング定数λは、後述する重み係数「Ｗ_ｉ ^（ｋｉ）×Ｗ_ｊ ^（ｋｊ）」との兼ね合いから「λ＝１」とする。

収束時刻特定部８１は、第１〜第ｎの尤度計算器６２−１〜６２−ｎにおいて算出された尤度Ｌ_ｉ（ｔ）（ｉ＝１〜ｎ）の分布間距離を求めると共に、該距離が最小となる時刻（収束時刻）を特定（判定）するものである。尤度Ｌ_ｉ（ｔ）（ｉ＝１〜ｎ）の分布間距離を求めるに際して収束時刻特定部８１は、各マルチパス波の内の２つのマルチパス波を複数組み合わせた上で、シグマポイント算定部７０により算定された各シグマポイントを用いて該各組み合わせの尤度分布間の距離を計算する。

より具体的に、収束時刻特定部８１は、図６に示すように、第１〜第ｎのシグマポイント計算器７０−１〜７０−ｎの内から２つを選ぶ全組み合わせ（各シグマポイント計算器の全組み合わせ）に相当する数の分布間距離計算器８１−（１、２）〜８１−（ｎ−１、ｎ）を有しており、これら各分布間距離計算器８１−（１、２）〜８１−（ｎ−１、ｎ）が、各々個別に分布間距離（尤度分布間の距離）を求めるものである。

ここで、各シグマポイント計算器の全組み合わせに相当する数は、図６に示す通り、「１+２+…+ｎ−１＝（ｎ−１）・ｎ／２」であり、すなわち、本実施形態では、収束時刻特定部８１内に、各分布間距離計算器８１−（１、２）〜８１−（ｎ−１、ｎ）が「（ｎ−１）・ｎ／２」個設けられている。各分布間距離計算器８１−（１、２）〜８１−（ｎ−１、ｎ）は、各々に応じた分布間距離 ｄ（ｉ、ｊ）を、下記式（７）及び式（８）に基づいて計算する。この分布間距離ｄ（ｉ、ｊ）は、ｉ番目の送出波とｊ番目の送出波（ｉ≠ｊ）の尤度分布の間の距離である。

この式（８）において、Ｗ_ｉ ^（ｋｉ）×Ｗ_ｊ ^（ｋｊ）は重み係数であり、^ｔ（ｘ_ｉ ^（ｋｉ）−ｘ_ｊ ^（ｋｊ））は（ｘ_ｉ ^（ｋｉ）−ｘ_ｊ ^（ｋｊ））の転置行列である。もっとも、上記式（８）では、Ｗ_ｉ ^（ｋｉ）についてのみを示したが、Ｗ_ｊ ^（ｋｊ）についてもこれと同様である。

ここで、上記式（６）〜（８）によって求められる各シグマポイント間の距離（シグマポイント間距離）の加重平均に基づく分布間距離のイメージを図７に示す。ここでは、第ｉ番目のマルチパス波における各シグマポイントｘ_ｉ ^（０）〜ｘ_ｉ ^（４）と、第ｊ番目のマルチパス波における各シグマポイントｘ_ｊ ^（０）〜ｘ_ｊ ^（４）と、を簡単のため５つずつ示すと共に、これらを相互に組み合わせた場合における上記式（７）の重み係数の一部を例示する。例えば、重み係数Ｗ０４は「（１/３）×（１/６）」、重み係数Ｗ００は「（１/３）×（１/3）」、重み係数Ｗ３２は「（１/６）×（１/６）」となっている。

このように、本実施形態では、尤度分布の中心に位置するシグマポイント（ｘ_ｉ ^（０）、ｘ_ｊ ^（０））に大きな重み（１／３）を付与し、中心以外のシグマポイントには相対的に小さな重み（１／６）を付与すると共に、収束時刻特定部８１が、尤度分布における中心の重みと中心以外の重みとの積である重み係数を用いて、収束時刻の判定に影響する分布間距離を求めるように構成したため、音源位置の推定精度の向上を図ることができる。図７では、尤度分布の中心以外に位置する各シグマポイントに共通の重みを付与する構成を示したが、中心からの距離に応じて各シグマポイントの重みを変化させるようにしてもよい。

各分布間距離計算器８１−（１、２）〜８１−（ｎ−１、ｎ）において個別に算出された分布間距離は、各々から収束時刻計算器８１Ａに出力される。収束時刻計算器８１Ａは、各分布間距離計算器８１−（１、２）〜８１−（ｎ−１、ｎ）の各々から入力した分布間距離の総和ｄ_total（t）を、下記式（９）によって求めるものである。この分布間距離の総和ｄ_total（t）は、逆伝搬時間の関数である。

さらに、収束時刻計算器８１Ａは、分布間距離の総和ｄ_total（t）を用いて、逆伝搬時間軸上の収束時刻ｔ_{coｎvergeｎce}を、下記式（１０）によって求める。

そして、収束時刻特定部８１は、収束時刻ｔ_{coｎvergeｎce}を出力端子８１Ｂを介して音源位置推定部８３に出力する。

統合尤度計算部８２は、第１〜第ｎの尤度計算器６２−１〜６２−ｎの各々で計算した各マルチパス波の各逆伝搬時間の尤度Ｌ_１（ｔ）〜Ｌ_ｎ（ｔ）を入力すると共に、各尤度Ｌ_１（ｔ）〜Ｌ_ｎ（ｔ）を同じ逆伝搬時間のもの同士で統合（合算）して時系列情報を作成する。すなわち、統合尤度計算部８２は、上記統合の結果である統合尤度（尤度関数：Ｌｓ（ｔ））を、各逆伝搬時刻について下記式（１１）を用いて算出するものである。

音源位置推定部８３は、収束時刻特定部８１から収束時刻ｔ_{coｎvergeｎce}を入力すると共に、統合尤度計算部８２から統合尤度Ｌｓ（ｔ）を受け取り、収束時刻と統合尤度に基づいて、収束時刻の尤度関数の中で（ｒ、ｚ）に関して最も尤度が大きくなる場所を音源位置と推定し、その結果を出力端子８４を介して外部に出力するものである。

（動作説明）
続いて、本実施形態における音源位置推定装置３０の動作を、図８に示すフローチャートに基づいて説明する。まず、受波器群２０が第１受波器２０Ａ及び第２受波器２０Ｂによって、音源９０から放射された到来波を受信すると共に、該到来波の時間波形の情報等を受波情報算出部４１に出力する（図８：ステップＳ１０１）。次いで、受波情報算出部４１が、第１受波器２０Ａ及び第２受波器２０Ｂの各々から入力された時間波形の情報等を解析して受波解析情報（到来時間差、到来ふ仰角、到来時間差誤差標準偏差、及び到来ふ仰角誤差標準偏差）を算出する（図８：ステップＳ１０２）。

次に、到来時間反転部４２が、受波情報算出部４１から入力された到来時間差を上記式（１）により反転し、反転後の反転時間差を他の到来波に係る情報と共に伝搬計算処理部６１へと出力する（図８：ステップＳ１０３）。伝搬計算処理部６１が、到来時間反転部４２から入力された各情報及び情報格納部５０から取得した各情報に基づいて、上記式（２）及び式（３）を用いた音線計算等を実行し、予め設定された時間間隔における各マルチパス波の到来位置、受波器群２０からの経路長等の情報を尤度算出部６２に出力する（図８：ステップＳ１０４）。

次いで、尤度算出部６２が、上記式（５）の処理により各マルチパス波の尤度Ｌ_ｉ（ｔ）（ｉ＝１〜ｎ）を算出すると共に、尤度Ｌ_ｉ（ｔ）をシグマポイント算定部７０及び統合尤度計算部８２に出力する（図８：ステップＳ１０５）。次に、シグマポイント算定部７０が、上記式（６）の処理により第ｉ番目のマルチパス波の各尤度分布における各シグマポイントを算定すると共に、算定した各シグマポイントを収束時刻特定部８１に出力する（図８：ステップＳ１０６）。

各シグマポイントが入力された収束時刻特定部８１は、各分布間距離計算器８１−（１、２）〜８１−（ｎ−１、ｎ）によって上記式（７）及び（８）の処理を行うことで、組み合わせ可能な２つのマルチパス波の分布間距離を計算すると共に、収束時刻計算器８１Ａにおいて上記式（９）に係る処理を行うことで、分布間距離の総和を求める（図８：ステップＳ１０７）。

そして、収束時刻特定部８１は、上記式（１０）を用いた処理により、分布間距離の総和が最小となる時刻（収束時刻）を特定する（図８：ステップＳ１０８）。一方、統合尤度計算部８２は、上記式（１１）を用いて、尤度算出部６２から入力された尤度Ｌ_ｉ（ｔ）（ｉ＝１〜ｎ）を同じ逆伝搬時間のもの同士で統合（合算）した統合尤度を算出する（図８：ステップＳ１０９）。

音源位置推定部８３は、収束時刻特定部８１と統合尤度計算部８２とから入力された収束時刻と統合尤度とをもとに音源９０の位置を推定し、推定した音源位置を出力端子８４から外部に出力する。より具体的に音源位置推定部８３は、収束時刻特定部８１から入力された収束時刻における尤度の合計が最大となる位置を音源９０の位置と推定し、推定結果を外部に出力する（図８：ステップＳ１１０）。

（実施形態の効果）
本実施形態における音源位置推定装置３０では、各マルチパス波の尤度分布を表す情報をもとに、それぞれの尤度分布のシグマポイントを算定し、算定した各シグマポイント間の距離に基づいて分布間距離を計算するという構成を採ったため、高精度な収束時刻の判定処理を実現することができる。そして、音源位置推定部８３は、収束時刻での全ての尤度分布の中で、尤度が最大となる位置を音源位置と推定するように構成されているため、音源位置推定装置３０によれば、音源の位置を極めて精度よく推定することが可能となる。

ここで、本実施形態における上記構成によって、シグマポイント間距離が連続性を保つ状況を、図９を参照して説明する。図９では、海底反射が起こる直前の逆伝搬時間ｔ_ｋ及び海底反射が起こった直後の逆伝搬時間ｔ_ｋ＋１を、それぞれ破線で囲んだ領域で示すと共に、逆伝搬時間ｔ_ｋ、ｔ_ｋ＋１における２つの尤度分布であるｄｉｓ１とｄｉｓ２とをそれぞれ上下に示す。また、各尤度分布ｄｉｓ１、ｄｉｓ２は、黒丸で示した尤度の中心と、楕円状の一点鎖線で囲った誤差分布とから成る。さらに、各尤度分布ｄｉｓ１、ｄｉｓ２に対応する各マルチパス波の逆伝搬経路を示す音線ＳＲ１、音線ＳＲ２を例示する。

図９のように、尤度分布ｄｉｓ２の海底反射は起こっているが、図１０に示す通り、海底反射の前後において分布間距離の総和の連続性が担保されている。後述する従来技術では、海底反射の影響により、分布間距離の総和の不連続が存在していたが（図１４参照）、本実施形態における状況を示す図１０のグラフでは、不連続点がなくなっており、連続性が改善されていることがわかる。

特に、収束時刻特定部８１は、各シグマポイント間の距離をもとに算出する分布間距離を、すべてのシグマポイントの組み合わせに対してユークリッド距離を求めた上で加重平均をとって算出するため、尤度分布の形状を表現した状態で、分布間距離を連続的に変化させることができる。すなわち、図９では、尤度分布ｄｉｓ２が連続的に海底反射する様子が、各シグマポイントの位置によって表現されており、尤度分布ｄｉｓ１の各シグマポイントと尤度分布ｄｉｓ２の各シグマポイント間の距離の加重平均を計算することで、分布間距離も連続的に変化することとなる（後述する第２の問題点の解決）。

加えて、伝搬計算処理部６１による伝搬経路の推定をもとに尤度算出部６２が算出した尤度分布に基づいて、シグマポイント算定部７０が、各マルチパス波の尤度分布のシグマポイントを算定し、収束時刻特定部８１が、各シグマポイントを用いて各尤度分布間の距離を計算するという構成を採ったことから、共分散行列によって規格化するという処理が不要となったため、誤差分布が広がることによって距離が近づいている旨を誤判定するといった問題は生じない。

すなわち、図１０に示す逆伝搬時間領域Ｔａにおける変動の様子からも解るように、収束時刻特定部８１は、尤度分布が近づけば距離が近い旨を、尤度分布が遠のけば距離が遠い旨を高精度で判定することができる（後述する第１の問題点の解決）。また、逆伝搬した全てのマルチパス波が音源位置に到達する理想的な時刻（理想的な収束時刻）は逆伝搬時間ｔ_０であるため、図１０に示すグラフから、収束時刻特定部８１が、収束時刻の判定を精度よく実行できることを確認できる。

続いて、誤検出が含まれる場合における、シグマポイント間距離に基づく分布間距離の総和と逆伝搬時間との関係を図１１に示す。図１１の場合にも、誤検出による分布間距離の寄与が含まれているため、分布間距離の総和は一定以上の値を下回ることはないが、後述する従来技術の場合とは異なり（図１５参照）、不連続な値の変動が全域で存在しないことが確認できる。

また、後述する従来技術では、誤検出の影響による不連続な変動が存在し、理想的な収束時刻である逆伝搬時間ｔ_０において、分布間距離の総和が最小値にならないという不都合が生じる（図１５参照）。しかし、本実施形態では、図１１に示す通り、誤検出が発生したような場合にも不連続な変動がなくなるため、妥当な逆伝搬時間ｔ_０に分布間距離の総和が最小値となることから、収束時刻特定部８１は、理想的な逆伝搬時間ｔ_０を収束時刻として求めることができる。

ここで、音源位置推定装置３０によって得られる効果を更に詳細に説明するための比較例を、図１２〜図１５を参照して説明する。２つの尤度分布の中心と誤差分布のイメージを例示する図１２は、マハラノビス距離の考え方による分布間距離の判定処理に係る説明図である。図１２には、到来時間差の平均値及び到来ふ仰角の平均値に基づいて求めた尤度（尤度分布）の中心Ｃ１〜Ｃ８と、中心Ｃ１〜Ｃ８に対して標準偏差の値だけずれた位置にある誤差分布Ｅ１〜Ｅ８（楕円状の一点鎖線）とを例示する。特許文献１における技術では、各分布間距離の遠近を予め決められた距離（判定基準距離）との比較により判定するに際して誤差分布Ｅ１〜Ｅ８が考慮される。

図１２（ａ）の場合は、２つの尤度分布の中心Ｃ１と中心Ｃ２とは近くにあり、誤差分布Ｅ１と誤差分布Ｅ２は広がっていない。このため、マハラノビス距離の考え方によれば、２つの尤度分布は近いと判定される。図１２（ｂ）の場合は、２つの尤度分布の中心Ｃ３と中心Ｃ４とは遠くにあり、誤差分布Ｅ３と誤差分布Ｅ４は広がっていないため、２つの尤度分布は遠いと判定される。

図１２（ｃ）の場合は、２つの尤度分布の中心Ｃ５と中心Ｃ６とは遠くにあるが、誤差分布Ｅ５と誤差分布Ｅ６は各中心を結ぶ方向（互いの分布の中心方向）に沿って広がっている。このため、マハラノビス距離の考え方によれば、中心Ｃ５と中心Ｃ６との間の距離にかかわらず、２つの分布は近いと判定される。

図１２（ｄ）の場合は、２つの尤度分布の中心Ｃ７と中心Ｃ８とは遠くにあり、その誤差分布が広がっている。しかし、この場合における誤差分布Ｅ７と誤差分布Ｅ８の広がりは、各中心を結ぶ方向とは異なる方向に沿っているため、マハラノビス距離の考え方によれば、２つの分布は遠いと判定される。

図１２（ｂ）〜（ｄ）では、中心間距離が相互に等しい状況を想定しているが、マハラノビス距離により分布間距離を算出すると、尤度分布の中心間距離と誤差分布間距離との差が大きくなり得るため、特に図１２（ｃ）の場合のように、分布間距離の遠近を正しく判定できないという不都合が起こる。

なお、図１２では、誤差分布が、尤度分布の中心を結ぶ直線にほぼ沿う方向又はこれに垂直な方向に広がっている状況を例示したが、実際には、誤差分布の広がりの方向はこれらに限定されないため、尤度分布の中心間距離と誤差分布間距離との間に一定の関係が成立しないことから、尤度分布の遠近の判定結果は、誤差分布の広がりに大きく左右される非常に不安定なものになる。

したがって、分布間距離の算出にマハラノビス距離の考え方を採用した場合には、２つの尤度分布の中心がどこにあるのかという点の他に、共分散行列によって表現される誤差分布の広がりや向きの影響が非常に大きくなる。すなわち、マハラノビス距離の持つ特性には、マルチパス波の収束時刻の判定において、特に下記３つの問題がある。

まず第１の問題点は、図１２（ｃ）に示すように、誤差分布が広がることによって、２つの尤度分布の中心が離れているにもかかわらず、各分布間の距離が近いと判定してしまうことである。すなわち、２つの尤度分布の中心は離れているが、該各中心を結ぶ方向において、その中心間の距離と同程度に誤差分布が広がっている。このため、２つの尤度分布を表す共分散行列で規格化してしまうと、２つの尤度分布は近いと判定されてしまうため、２つの確率変数間の距離の期待値を、誤差分布の形状を考慮した上で表現して収束判定に用いることができない。

第２の問題点は、海面反射や海底反射によって、マハラノビス距離が不連続に変化するという不都合である。特許文献１に開示された手法によれば、逆伝搬しているマルチパス波が海面もしくは海底で反射した際に、マハラノビス距離に関数としての連続性が生じる。図１３は、２つのマルチパス波が逆伝搬している様子を示す模式図である。

図１３では、海底反射が起こる直前の逆伝搬時間ｔ_ｋ及び海底反射が起こった直後の逆伝搬時間ｔ_ｋ＋１を左右に破線で示すと共に、各逆伝搬時間ｔ_ｋ、ｔ_ｋ＋１における２つの尤度分布であるｄｉｓ０１とｄｉｓ０２とをそれぞれ上下に示している。また、各尤度分布ｄｉｓ０１、ｄｉｓ０２は、黒丸で示した尤度の中心と、楕円状の一点鎖線で囲った誤差分布と、から成る。

さらに、図１３では、各尤度分布ｄｉｓ０１、ｄｉｓ０２に対応する各マルチパス波の逆伝搬経路を示す音線ＳＲ０１、音線ＳＲ０２を例示している。尤度分布ｄｉｓ０２が海底反射した際に、マハラノビス距離の規格化に使用する２つの共分散行列のうち、尤度分布ｄｉｓ０２の共分散行列が不連続に変化するため、図１４に示すように、上記図１０に示す本実施形態の場合とは異なり、マハラノビス距離も不連続となる。

すなわち、尤度分布ｄｉｓ０２が海底反射した際に、マハラノビス距離の規格化に使用する尤度分布ｄｉｓ０２の共分散行列の項のうち、該分布の傾きを表す非対称項が、海底反射によって逆伝搬時間ｔ_ｋと逆伝搬時間ｔ_ｋ＋１とで不連続となる。このため、規格化に共分散行列を用いるマハラノビス距離も逆伝搬時間ｔ_ｋと逆伝搬時間ｔ_ｋ＋１で不連続となる。

ここで、逆伝搬時間に対するマハラノビス距離による分布間距離を示す図１４を参照して、上記２つの問題が起こっている状況を説明する。図１４には、逆伝搬時間ｔ_ｋ、ｔ_ｋ＋１で、図１３における尤度分布ｄｉｓ０２に該当するマルチパス波が海底反射しているグラフを示す。

このグラフには、逆伝搬時間ｔ_ｋ、ｔ_ｋ＋１（例えば０．０１８ｓｅｃ）の前後である逆伝搬時間ｔ_ｋから逆伝搬時間ｔ_ｋ＋１に亘って、分布間距離の総和が急激に増加していることが示されており、海底反射時にマハラノビス距離が不連続に変化している様子がわかる。この点、尤度分布ｄｉｓ２の海底反射が連続的である本実施形態とは異なる（図１０参照）。なお、理想的な収束時刻は逆伝搬時間ｔ_０である。

また、図１４では、海底反射後における逆伝搬時間領域Ｔｂのマハラノビス距離は、徐々に小さくなるという傾向を示している。このことから、図１２（ｃ）のように、互いの尤度分布の中心方向に沿って誤差分布が広がることによって、分布間距離が実際の距離よりも近いと判定してしまう事態が起こっていることを確認できる。すなわち、誤差分布が広がることによって、尤度分布の中心間距離の如何によらず、分布間距離が近づいているとの判定が行われていることがわかる。

ところで、マルチパス波の検出に誤検出が含まれていない場合には、検出精度の問題を除けば、ある逆伝搬時間に分布の中心が重なることが期待されるため、妥当な収束時刻においては、基本的にすべての逆伝搬するマルチパス波の分布間距離の総和が０となるはずである。しかしながら、誤検出が含まれている場合には、正当な検出によるマルチパス波と誤検出によるマルチパス波との分布間距離が、適切な値との隔たりを持つことになる。

すなわち、特許文献１に開示されたマハラノビス距離に係る手法によって分布間距離を求めた場合には、分布間の中心が離れている場合であっても、正当な検出によるマルチパス波と誤検出によるマルチパス波との分布間距離を近いと誤判定するといった問題や、海面や海底での反射によって分布間距離が不連続に変化するといった問題が起こる。こうした問題は、分布間距離の総和が最小となる逆伝搬時間（収束時刻）を求める際に大きな悪影響を及ぼす。

第３の問題点は、マルチパス波の検出に誤検出が含まれていた場合に、上記第１及び第２の問題発生が影響して、分布間距離の妥当な最小値を推定できない頻度が増えることである。ここで、マルチパス波の検出に誤検出が含まれている状況を示す図１５を参照して第３の問題点を説明する。

図１５では、全域に亘って誤検出が含まれている状況を想定しているため、各分布間距離の総和は、一定以上の０ではない分布間距離が不連続に変動しながら加算されたものである。図１５に示すように、理想的な収束時刻である逆伝搬時間ｔ_０での分布間距離の総和は、局所的には小さな値となっているが、全域的には最小値とはなっていない。

すなわち、誤検出の影響によって、逆伝搬時間０.０１５ｓｅｃ以降における分布間距離の総和の方が、理想的な収束時刻における分布間距離の総和よりも小さな値となっており、しかも伝搬時間０.０３ｓｅｃを過ぎた逆伝搬時間ｔ_ｂの時点において分布間距離の総和が最小であると判定される不都合が生じる。収束時刻の誤判定は、音源位置の推定結果の誤りに直結してしまうため、収束時刻の判定精度に影響する分布間距離を精度よく算出することは極めて重要である。

本実施形態における音源位置推定装置３０は、各マルチパス波の尤度分布から各シグマポイントを算定するシグマポイント算定部７０と、シグマポイント算定部７０において算定された各シグマポイント間の距離に基づいて各マルチパス波の分布間距離を計算する収束時刻特定部８１を有するため、収束時刻の判定を精度よく行うことができる（上記第１の問題点の解決）。また、収束時刻特定部８１内の各分布間距離計算器は、各シグマポイント間の距離の加重平均を尤度分布間の距離として求めるため、尤度分布間の距離が連続的に変化させることができる（上記第２の問題点の解決）。併せて、第１及び第２の問題点に関連する上記第３の問題点も改善することができる。

すなわち、尤度算出部６２の出力側にシグマポイント算定部７０を設け、シグマポイント算定部７０で算定した各シグマポイントに基づいて収束時刻特定部８１が収束時刻を判定し特定する、という構成を採った音源位置推定装置３０によれば、受波器群２０から入力した音波の情報をより有意に解析することで、収束時刻に係る誤判定を低減させることが可能となり、かつ該分布間距離が海面や海底での反射によって不連続に変動するといった状況を回避することができる。

なお、上述した実施形態は、音源位置推定装置、その方法における好適な具体例であり、技術的に好ましい種々の限定を付している場合もあるが、本発明の技術的範囲は、特に本発明を限定する記載がない限り、これらの態様に限定されるものではない。

例えば、収束時刻特定部８１が、シグマポイント算定部７０によって算定された各シグマポイント間の距離の平均を尤度分布間の距離として求めるようにしてもよい。また、収束時刻特定部８１の構成内容は、図６に示すものに限らず、すなわち、収束時刻特定部８１内に任意の数（１以上）の分布間距離計算器を設け、任意の数の分布間距離計算器が、第１〜第ｎのシグマポイント計算器から入力した各シグマポイントを相互に組み合わせて各尤度分布の間の距離を算出するようにしてもよい。

さらに、本実施形態では、図７等を参照して、各シグマポイント計算器７０−１〜７０−ｎが各尤度分布におけるシグマポイントを５つ算出するという内容を例示した。しかし、上記式（６）にも示す通り、各尤度分布におけるシグマポイントの個数は任意に設定することができるため、シグマポイント算定部７０が、それぞれの尤度分布において任意の個数（２以上）のシグマポイントを算出するようにしてもよい。

１ 直接波、２〜５ マルチパス波、ＳＲ１、ＳＲ２、ＳＲ０１、ＳＲ０２ 音線、１０ 音源位置推定システム、２０ 受波器群、２０Ａ 第１受波器、２０Ｂ 第２受波器、３０ 音源位置推定装置、４１ 受波情報算出部、４２ 到来時間反転部（到来時間差反転部）、５０ 情報格納部、５０Ａ 音速プロファイル、５０Ｂ 海底地形情報、５０Ｃ 受波器位置情報、５５ 実測処理部、６１ 伝搬計算処理部、６１−１〜６１−ｎ 第１〜第ｎの伝搬計算器、６２ 尤度算出部、６２−１〜６２−ｎ 第１〜第ｎの尤度計算器、７０ シグマポイント算定部、７０−１〜７０−ｎ 第１〜第ｎのシグマポイント計算器、８１ 収束時刻特定部、８１−（１、２）〜８１−（ｎ−１、ｎ） 分布間距離計算器、８１Ａ 収束時刻計算器、８１Ｂ 出力端子、８２ 統合尤度計算部、８３ 音源位置推定部、８４ 出力端子、９０ 音源、Ｃ１〜Ｃ８ 中心、Ｅ１〜Ｅ８ 誤差分布、尤度Ｌ_１（ｔ）〜Ｌｎ（ｔ）：尤度Ｌｉ（ｔ）（ｉ＝１〜ｎ） 尤度、Ｌｓ（ｔ） 統合尤度、Ｔａ、Ｔｂ 逆伝搬時間領域、Ｗ_ｉ ^（ｋｉ）×Ｗ_ｊ ^（ｋｊ）、Ｗ００、Ｗ０４、Ｗ３２ 重み係数、ｄ（ｉ、ｊ） 分布間距離、ｄｉｓ１、ｄｉｓ２、ｄｉｓ０１、ｄｉｓ０２ 尤度分布、ｄ_total（t） 総和、ｚ 深度方向、ｒ 受波器群から音源に向かう方向、ｓ 媒介変数、ｔ、ｔ_０、ｔ_ｂ、ｔ_ｉ、ｔ_ｋ、ｔ_ｋ＋１ 逆伝搬時間、ｔ_{coｎvergeｎce} 収束時間、ｘ_ｉ ^（０）〜ｘ_ｉ ^（４）、ｘ_ｊ ^（０）〜ｘ_ｊ ^（４） シグマポイント、θ 送出ふ仰角、Ｃ 受波器群の位置における音速、ζ ｓｉｎθ／Ｃ、ξ ｃｏｓθ／Ｃ、λ スケーリング定数、τ'_ｎ、τ'_ｎ−１、τ'_ｎ−２ 反転時間差、τ_ｎ、τ_ｎ−１、τ_ｎ−２ 到来時間差、Ｚｗ 水深、Ｚｒ 海面から受波器群までの垂直距離、Ｚｓ 海面から音源までの垂直距離、Ｒｗ 音源と音源位置推定装置との間の水平距離、Ｄｈ 間隔。

Claims (8)

1. 音源から受信した複数のマルチパス波に基づいて、各マルチパス波の到来時間及び到来ふ仰角の情報を算出する受波情報算出部と、
   前記各マルチパス波の到来時間を反転する到来時間反転部と、
   前記到来時間反転部により反転された各到来時間と前記受波情報算出部により算出された到来ふ仰角の情報とに基づいて、各マルチパス波の各位置の尤度分布を算出する尤度算出部と、
   前記尤度算出部により算出された各マルチパス波の各位置の尤度分布について複数のシグマポイントを算定するシグマポイント算定部と、
   前記シグマポイント算定部により算定された各シグマポイントを用いて、前記各マルチパス波の内の２つのマルチパス波を複数組み合わせて該各組み合わせの尤度分布間の距離を、一方のマルチパス波における前記シグマポイントと他方のマルチパス波における前記シグマポイントとの複数の組み合わせについて計算する分布間距離計算部と、
   前記分布間距離計算部により計算された各尤度分布間の距離の総和を算出し、該総和が最小となる時刻を特定する収束時刻特定部と、
   前記収束時刻特定部により特定された時刻に基づいて前記音源の位置を推定する音源位置推定部と、を有することを特徴とする音源位置推定装置。
2. 前記分布間距離計算部は、前記各シグマポイント間の距離の加重平均を前記尤度分布間の距離として求めるものであることを特徴とする請求項１に記載の音源位置推定装置。
3. 前記分布間距離計算部は、前記各シグマポイントの全組み合わせに相当する数の分布間距離計算器を有し、該各分布間距離計算器各々によって前記各尤度分布間の距離を求めるものであることを特徴とする請求項１又は２に記載の音源位置推定装置。
4. 前記シグマポイント算定部は、前記各位置の尤度分布と該各位置の尤度分布の共分散行列とに基づいて前記各シグマポイントを算定するものであることを特徴とする請求項１〜３の何れか１項に記載の音源位置推定装置。
5. 前記尤度算出部は、前記各マルチパス波の各位置の尤度分布をそれぞれ個別に計算する複数の尤度計算器を有し、
   前記シグマポイント算定部は、前記各尤度計算器に個別に対応する複数のシグマポイント計算器を有することを特徴とする請求項１〜４の何れか１項に記載の音源位置推定装置。
6. 前記音源位置推定部は、前記収束時刻特定部において特定された時刻における尤度の合計が最大となる位置を前記音源の位置と推定するものであることを特徴とする請求項１〜５の何れか１項に記載の音源位置推定装置。
7. 音源から受信した複数のマルチパス波に基づいて、各マルチパス波の到来時間及び到来ふ仰角の情報を算出し、
   前記各マルチパス波の到来時間を反転し、
   前記反転された各到来時間と前記算出された到来ふ仰角の情報とに基づいて、各マルチパス波の各位置の尤度分布を算出し、
   前記算出された各マルチパス波の各位置の尤度分布について複数のシグマポイントを算定し、
   前記各マルチパス波の内の２つのマルチパス波を複数組み合わせた上で、前記算定された各シグマポイントを用いて前記各組み合わせの尤度分布間の距離を、一方のマルチパス波における前記シグマポイントと他方のマルチパス波における前記シグマポイントとの複数の組み合わせについて計算し、
   前記計算された各尤度分布間の距離の総和を算出すると共に、該総和が最小となる時刻を特定し、
   前記特定された時刻に基づいて前記音源の位置を推定することを特徴とする音源位置推定方法。
8. 音源から受信した複数のマルチパス波に基づいて、各マルチパス波の到来時間及び到来ふ仰角の情報を算出する受波情報算出手段、
   前記各マルチパス波の到来時間を反転する到来時間反転手段、
   前記到来時間反転手段により反転された各到来時間と前記受波情報算出手段により算出された到来ふ仰角の情報とに基づいて、各マルチパス波の各位置の尤度分布を算出する尤度算出手段、
   前記尤度算出手段により算出された各マルチパス波の各位置の尤度分布について複数のシグマポイントを算定するシグマポイント算定手段、
   前記各マルチパス波の内の２つのマルチパス波を複数組み合わせた上で、前記シグマポイント算定手段により算定された各シグマポイントを用いて前記各組み合わせの尤度分布間の距離を、一方のマルチパス波における前記シグマポイントと他方のマルチパス波における前記シグマポイントとの複数の組み合わせについて計算する分布間距離計算手段、
   前記分布間距離計算手段により計算された各尤度分布間の距離の総和を算出すると共に、該総和が最小となる時刻を特定する収束時刻特定手段、
   及び前記収束時刻特定手段により特定された時刻に基づいて前記音源の位置を推定する音源位置推定手段、としてコンピュータを機能させることを特徴とする音源位置推定プログラム。

Images