音源位置推定装置
【課題】到来時間差観測値に誤差がある場合にも、音源位置の推定誤差を小さくする。
【解決手段】到来ふ仰角が上向きのマルチパス波と到来ふ仰角が下向きのマルチパス波について別個に到来時間差のテーブル値を記憶しておき（２０９Ｍ）、受信したマルチパス波の各々について到来ふ仰角が上向きか下向きかの判定を行い、観測された到来時間差とテーブル値との不一致度を上向きのふ仰角を有するマルチパス波と下向きのふ仰角を有するマルチパス波の各々について別個に算出して、その和をコスト関数として（２１０）、音源位置の推定を行う（２１１）。
【解決手段】到来ふ仰角が上向きのマルチパス波と到来ふ仰角が下向きのマルチパス波について別個に到来時間差のテーブル値を記憶しておき（２０９Ｍ）、受信したマルチパス波の各々について到来ふ仰角が上向きか下向きかの判定を行い、観測された到来時間差とテーブル値との不一致度を上向きのふ仰角を有するマルチパス波と下向きのふ仰角を有するマルチパス波の各々について別個に算出して、その和をコスト関数として（２１０）、音源位置の推定を行う（２１１）。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は音源位置推定装置に関し、特に水中の音源から放射されたパルス状の音波の直接波とマルチパス波を受信し、これらの到来時間差から音源の位置を推定する装置に関する。
【背景技術】
【０００２】
水中の音波伝搬では、音源から直接受波器に到来する音波（直接波）の他に、海面や海底で反射して到来するマルチパス波があり、そのようなマルチパス波は、直接波と異なる時刻に到来する。
水中において音源からのパルス状の音波を受信して音源の位置を推定する技術として、下記の特許文献１に示されるものがある。
【０００３】
図１に、特許文献１に示される、従来の音源位置推定装置を示す。図示の音源位置推定装置においては、音響センサー（受波器）で受信した音波を表す信号が入力端子１０１を介してマルチパス波到来時間差測定器１０２に入力される。マルチパス波到来時間差測定器１０２では、直接波とマルチパス波の到来時間差を検出し、検出した到来時間差観測値を示すデータを出力する。
【０００４】
一方、到来時間差テーブル作成器１０３により、水中の音速プロファイルに基づき、各仮想音源からの音波の音線計算を行うことにより音波の到来時間を推定し、推定結果を到来時間差テーブル値（理論値）として、記憶しておく。
【０００５】
コスト関数算出器１０４は、マルチパス波到来時間差測定器１０２で求めたマルチパス波到来時間差の観測値と到来時間差テーブル作成器１０３で求めた到来時間差テーブル値より、それらの不一致度を反映したコスト関数の値を算出する。
【０００６】
コスト関数算出器１０４の詳細ブロック構成を図２に示す。
入力端子１１１にはマルチパス波到来時間差測定器１０２で求めたマルチパス波到来時間差観測値が入力され、入力端子１１２には到来時間差テーブル作成器１０３からのテーブル値が入力される。テーブル値−観測値算出器１１３は、到来時間差観測値と到来時間差テーブル値との差を算出する。
【０００７】
ここで、ｐを受波位置（受波点）から音源位置までの距離に対応するインデックス、ｑを音源位置の深度に対応するインデックスとし、マルチパス波到来時間差測定器１０２で得られるｎ番目の到来波についての到来時間差をＤｎ、到来時間差テーブル作成器１０３で求めた、仮想音源位置（ｐ,ｑ）の仮想音源からの、ｎ番目の到来時間差テーブル値をＥｎ（ｐ，ｑ）とした場合、テーブル値−観測値差算出器２０３では、テーブル値と観測値の差を次の式（１）で求める。
△Ｄｎ（ｐ，ｑ）＝ Ｅｎ（ｐ，ｑ）−Ｄｎ …（１）
【０００８】
自乗和算出器１１４は、式（１）の演算結果を自乗したものの、すべてのマルチパス波についての総和（下記の式（２）で表される）を、コスト関数として求め、求めたコスト関数を出力端子１１５から出力する。
【０００９】
【数１】
【００１０】
図１の音源位置推定器１０５では式（２）で与えられるコスト関数が最小となる（ｐ，ｑ）を求め、求められた（ｐ，ｑ）が推定対象（目標）の音源位置を表すものであると推定し、推定結果を、出力端子１０６を介して出力する。
【先行技術文献】
【特許文献】
【００１１】
【特許文献１】特開２００６−１９４６２７号公報
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかし、上記の従来の装置を用いた場合、音源と受波器を結ぶ線上のコスト関数の値が、音源の位置の前後のかなり広い範囲にわたり小さな値となるため、音源位置を正確に推定することができないという問題があった。即ち、マルチパス波到来時間差測定器１０２で求められる到来時間差は観測誤差を含むため、コスト関数の値が小さい場所で大きく推定値がずれる可能性が大きく、このため観測誤差が比較的大きいという問題があった。
【００１３】
一例として、図３に示すように音速が一定で水深Ｚｂが一定である場合を考える。図３において、Ｚｓは音源の深度、Ｚｒは受波点の深度、ｒは受波点ＷＲから音源ＷＳまでの水平方向距離を示す。
到来時間差観測値に誤差がないとした場合のコスト関数を求めた結果を図４に示す。横軸は受波点Ｚｒからの水平方向距離ｒ、縦軸は深度（下方向ほど深度が大）ｚを表す。図中の曲線はコスト関数の等高線であり、このコスト関数が最小となる位置を推定する。
【００１４】
到来時間差観測値に誤差がない場合は、真値を推定することができるが、図中のコスト関数は、受波器と音源とを結ぶ線の方向（矢印ＤＭＣの方向）に沿って、コスト関数の変化が小さいために、到来時間差の観測値に誤差が含まれる場合は、この方向に沿って、真の音源位置から推定結果がずれる可能性がある。経験的に知られている観測誤差を与えて、図３に示した条件でシミュレーションを３０回試行したところ、距離の誤差標準偏差は、距離真値の約１５％と、かなり大きい値となった。このようにコスト関数の特質から音源の位置推定誤差が大きいという問題があった。
【課題を解決するための手段】
【００１５】
本発明の音源推定装置は、
受波位置における到来波に基づいて音源の位置を推定する音源推定装置において、
複数の仮想音源位置の各々について、該仮想音源位置に仮想音源が存在すると仮定したときの、前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が上向きのものの各々との到来時間差の理論値、並びに前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が下向きのものの各々との到来時間差の理論値を格納する到来時間差理論値格納手段と、
前記受波位置において観測した到来波のデータに基づいて、目標の音源からの最初の到来波と２番目以降の到来波との到来時間差の観測値を算出する到来時間差観測値算出手段と、
前記２番目以降の到来波の各々の到来ふ仰角が上向きか下向きかを判定し、前記２番目以降の到来波のうち、到来ふ仰角が上向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が上向きの到来波についての前記到来時間差の理論値との不一致度を第１の不一致度として求め、前記２番目以降の到来波のうち、到来ふ仰角が下向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が下向きの到来波についての前記到来時間差の理論値との不一致度を第２の不一致度として求め、前記第１の不一致度と前記第２の不一致度との和をコストとして求めるコスト関数算出手段と、
前記コスト関数算出手段が算出した前記コストが最小となる前記仮想音源位置を前記目標の音源の位置と推定する音源位置推定手段と
を備えたことを特徴とする。
【発明の効果】
【００１６】
本発明によれば、到来時間差観測値に誤差がある場合にも、音源位置の推定誤差を小さくすることができる。
【図面の簡単な説明】
【００１７】
【図１】従来の音源位置推定装置を示すブロック図である。
【図２】図１のコスト関数算出器１０４の詳細を示すブロック図である。
【図３】音速が一定で水深が一定である場合の音波の伝搬経路を示す図である。
【図４】図１の装置で求めたコスト関数の等高線の一例を示す図である。
【図５】本発明の実施の形態１による音源位置推定装置を示すブロック図である。
【図６】実施の形態１における受波器の配置を示す図である。
【図７】受波器におけるサンプル値の絶対値の包絡線の一例を示すグラフである。
【図８】（ａ）及び（ｂ）は、一対の受波器における受信のタイミングの一例を示す図である。
【図９】到来時間差テーブルの作成に当たり使用される仮想音源位置の一例を示す図である。
【図１０】仮想音源位置ＷＳの音源からの伝搬経路の一例を示す図である。
【図１１】図５の到来時間差テーブル作成器２０９の一例を示す機能ブロック図である。
【図１２】図５のコスト関数算出器２１０の一例を示す機能ブロック図である。
【図１３】図５の装置で求めたコスト関数の等高線の一例を示す図である。
【図１４】本発明の実施の形態２による音源位置推定装置を示すブロック図である。
【図１５】図１４のコスト関数算出器２２０の一例を示す機能ブロック図である。
【図１６】図１５のふ仰角判定部８２５の一例を示す機能ブロック図である。
【図１７】（ａ）及び（ｂ）は、一対の受波器における受信のタイミングの一例を示す図である。
【発明を実施するための形態】
【００１８】
実施の形態１．
本発明の実施の形態１の音源位置推定装置を図５に示す。
図示の音源位置推定装置は、音源から放射されたパルス状音波の直接波とマルチパス波を受信し、それらの到来時間差に基づいて音源の位置を推定するものであり、音源の推定位置は深度及び受波器からの水平方向距離で表される。なお、音源の位置する方位を示す情報は別途得られるものとする。
【００１９】
実施の形態１の音源位置推定装置は、第１及び第２の受波器２０１及び２０２を備えている。
第１の受波器２０１及び第２の受波器２０２は、後述のように、到来波のふ仰角が上向きか下向きかの情報を得るために設けられたものであり、同じ鉛直線上の異なる位置に設けられている。このように配置されることにより、第１の受波器２０１及び第２の受波器２０２のいずれに、到来波（信号）がより早く到来したかに基づいて、到来波のふ仰角が上向きか下向きかの判断を行うことができる。
【００２０】
第１及び第２の受波器２０１及び２０２は、ふ仰角が上向きか下向きの判別が行えれば良いので、図６に示すように、間隔ＤＷＲを開けて配置されている。間隔ＤＷＲは受信位置において想定される音速の最大値と受波器２０１、２０２で測定可能な最小時間差の積以上の値に設定される。
第１及び第２の受波器２０１及び２０２はそれぞれ上側及び下側に配置されるので、「上側の受波器」、「下側の受波器」と呼ばれることもある。
【００２１】
上側及び下側の受波器２０１及び２０２の出力は、入力端子２０３及び２０４を介してそれぞれ第１及び第２の到来時刻決定部２０５及び２０６に入力される。
第１及び第２の到来時刻決定部２０５及び２０６は、それぞれ上側及び下側の受波器２０１及び２０１の出力に基づいて到来時刻を決定乃至検出する。
【００２２】
第１及び第２の受波器２０１及び２０２は、それらで受信される音波をサンプリングすることで得られるサンプル値を表すデータを出力するが、音波をサンプリングすることで得られるサンプル値の絶対値を時間軸上に並べたものの包絡線は、例えば、図７に示す如くとなる。なお、サンプル値の絶対値の代りにサンプル値の自乗値の包絡線を用いる場合もある。
【００２３】
第１及び第２の到来時刻決定部２０５及び２０６では、それぞれ上側及び下側の受波器２０１及び２０２からのデータの列に基づき到来波の受信時刻（到来時刻）を判断する。到来時刻の判断には、従来と同様に、レプリカ相関による方法、立ち上がり検出による方法などを用いることができる。これらの方法は、例えば上記特許文献１においてマルチパス波到来時間差測定器の動作の一部として説明されている。
【００２４】
図８（ａ）及び（ｂ）に、第１及び第２の到来時刻決定部２０５及び２０６で決定乃至検出された上側の受波器２０１及び下側の受波器２０２におけるパルス状波の到来時刻（検出タイミング）を示す。図８（ａ）は、上側の受波器２０１における検出タイミングを、早い順にＴＡ１、ＴＡ２、…で示し、図８（ｂ）は、下側の受波器２０２における検出タイミングを、早い順にＴＢ１、ＴＢ２、…で示す。
ノイズを信号と誤って検出したり、検出漏れが発生しなければ、最初の検出タイミングが直接波の検出タイミングであり、２番目以降の検出タイミングがマルチパス波の検出タイミングである。
【００２５】
第１のマルチパス波到来時間差算出部２０７は、第１の到来時刻決定部２０５の出力に基づいて、上側の受波器２０１で受信したマルチパス波の各々についての到来時間差（１番目の到来波との到来時間差）を算出する。即ち、１番目の到来波ＰＡ１と２番目の到来波ＰＡ２（１番目のマルチパス波）の到来時間差ＤＵ２、１番目の到来波ＰＡ１と３番目の到来波ＰＡ３の到来時間差ＤＵ３、一般化すれば、１番目の到来波ＰＡ１とｎ番目の到来波ＰＡｎ（ｎ＝２、３、…）との到来時間差ＤＵｎを算出する。このようにして算出された到来時間差ＤＵｎを示すデータの列[ＤＵ２、ＤＵ３、…］はコスト関数算出器２１０に供給される。
【００２６】
同様に、第２のマルチパス波到来時間差算出部２０８は、第２の到来時刻決定部２０６の出力に基づいて、下側の受波器２０２で受信したマルチパス波の各々についての到来時間差（１番目の到来波との到来時間差）を算出する。即ち、１番目の到来波ＰＢ１と２番目の到来波ＰＢ２（１番目のマルチパス波）の到来時間差ＤＬ２、１番目の到来波ＰＢ１と３番目の到来波ＰＢ３の到来時間差ＤＬ３、一般化すれば、１番目の到来波ＰＢ１とｎ番目の到来波ＰＢｎ（ｎ＝２、３、…）との到来時間差ＤＬｎを算出する。このようにして算出された到来時間差ＤＬｎを示すデータの列[ＤＬ２、ＤＬ３、…］はコスト関数算出器２１０に供給される。
【００２７】
第１の到来時刻決定部２０５と第１のマルチパス波到来時間差算出部２０７の組合せは、図１のマルチパス波到来時間差測定器１０２と同様の機能を有し、第２の到来時刻決定部２０６と第２のマルチパス波到来時間差算出部２０８の組合せも、図１のマルチパス波到来時間差測定器１０２と同様の機能を有する。
【００２８】
上記のように、第１のマルチパス波到来時間差算出部２０７と第２のマルチパス波到来時間差算出部２０８とは、受波位置において観測した到来波の到来時刻（検出タイミング）に基づいて、推定対象（目標）の音源からの最初の到来波と２番目以降の到来波の各々との到来時間差の観測値を算出する到来時間差観測値算出手段を構成し、そのうち、第１のマルチパス波到来時間差算出部２０７が、第１の受波器２０１への最初の到来波と２番目以降の到来波の各々との到来時間差の観測値を算出し、第２のマルチパス波到来時間差算出部２０８が、第２の受波器２０２への最初の到来波と２番目以降の到来波の各々との到来時間差の観測値を算出する。
【００２９】
到来時間差テーブル作成器２０９は、入力端子２１３を介して供給される仮想音源の設定位置データと入力端子２１４を介して供給される音速プロファイルデータに基づいて、到来時間差テーブルを作成し、記憶部２０９Ｍに保存する。
【００３０】
到来時間差テーブルを作成する際の仮想音源位置は例えば、図９に示すように設定される。
図９に示す例では、探索領域を仮想的にメッシュ（格子）状に区切り、各交点を仮想音源位置（ｐ，ｑ）（ｐ＝１，２，…Ｐ、ｑ＝１，２，…Ｑ）とする。なお、ｐ、ｑはそれぞれ距離ｒ、深度ｚを表す値そのものではなく、距離、深度に対応するインデックスであり、自然数１〜Ｐ、１〜Ｑで表される。
【００３１】
メッシュ幅（仮想音源位置の間隔）は、所望の音源位置の推定精度に応じて設定される。メッシュ幅を狭くするほど高い精度で測定を行なうことができるが、データ量、演算量が多くなり、リアルタイムでの処理が困難となる。従って、高精度が要求される場合には、メッシュ幅を狭くし、高速性が重視される場合には、メッシュ幅を広くする。
【００３２】
仮想音源位置ＷＳからの伝搬経路は例えば図１０に示すように想定される。図１０には、簡単のため、５つの伝搬経路ＷＰ１〜ＷＰ５のみが示されている。伝搬経路（音線：音のエネルギーの流れ）が直線的でないのは、水温、塩分濃度などの不均一性のため、音速が一定ではないことによる。
一度も反射をせずに受波点ＷＲに達した波は直接波と呼ばれ、海面や海底で１回以上反射して受波器に達した波はマルチパス波と呼ばれる。
【００３３】
このような想定に基づいて、仮想音源位置の各々からの到来時間差を推測し、推測値（理論値）を時間差テーブル値として求め、記憶部２０９Ｍに記憶しておく。到来時間差テーブル作成器２０９は、仮想音源位置の各々でパルス状音波が発せられたと仮定した場合に、上側の受波器２０１における第１の到来波（直接波）の推定到来時刻と、上側の受波器２０１に上向きのふ仰角でｉ番目（ｉ＝１、２…）に到来すると推定されるマルチパス波の到来時刻との差（到来時間差）の推測値乃至理論値ＥＵｉ（ｉ＝１、２、…）の組乃至列[ＥＵ１、ＥＵ２、…ＥＵＩ]を、すべての仮想音源について求め、記憶部２０９Ｍに格納し、同様に、下側の受波器２０２における第１の到来波（直接波）の推定到来時刻と、下側の受波器２０２に下向きのふ仰角でｊ番目（ｊ＝１、２…）に到来すると推定されるマルチパス波の到来時刻との差（到来時間差）の推測値乃至理論値ＥＬｊ（ｊ＝１、２、…）の組乃至列[ＥＬ１、ＥＬ２、…ＥＬＪ]を、すべての仮想音源について求め、記憶部２０９Ｍに格納する。このように、記憶部２０９Ｍは、到来時間差理論値格納手段としての役割を果たす。
【００３４】
なお、Ｉは、経験的に受波器２０１で受信可能と推定される、ふ仰角が上向きのマルチパス波の最大数以上の値に、Ｊは、経験的に受波器２０２で受信可能と推定される、ふ仰角が下向きのマルチパス波の最大数以上の値に、それぞれ定められる。
【００３５】
全ての仮想音源位置（ｐ，ｑ）（ｐ＝１、２、…Ｐ、ｑ＝１、２、…Ｑ）の各々についての、上側の受波器２０１での到来時間差の列
ＥＵ１（ｐ，ｑ）、ＥＵ２（ｐ，ｑ）、…（ｐ＝１、２、…Ｐ、ｑ＝１、２、…Ｑ）
が得られ、これらがテーブル形式のデータとされ、上向きふ仰角到来時間差テーブル値として到来時間差テーブル作成器２０９の記憶部２０９Ｍ内に記憶されている。
同様に、全ての仮想音源位置（ｐ，ｑ）の各々についての、下側の受波器２０２での到来時間差の列
ＥＬ１（ｐ，ｑ）、ＥＬ２（ｐ，ｑ）、…（ｐ＝１、２、…Ｐ、ｑ＝１、２、…Ｑ）
が得られ、これらがテーブル形式のデータとされ、下向きふ仰角到来時間差テーブル値として到来時間差テーブル作成器２０９の記憶部２０９Ｍ内に記憶されている。
【００３６】
コスト関数算出器２１０は、第１及び第２の到来時刻決定部２０５、２０６の出力（到来時刻を示す信号）を入力とするとともに、第１及び第２のマルチパス波到来時間差算出部２０７及び２０８で得られた到来時間差データ（観測値）と、到来時間差テーブル作成器２１０で作成された到来時間差テーブルの値（理論値）を入力とし、これらに基づいてコスト関数を求め、求めたコスト関数を音源位置推定器２１１に入力する。
【００３７】
音源位置推定器２１１では、コスト間数算出器２１０で作成したコスト関数が最小となる仮想音源の距離及び深度を音源位置の推定値とし、該推定値を出力端子２１２から出力する。
【００３８】
以下、到来時間差テーブル作成器２０９の動作を、図１１を参照してより詳細に説明する。
入力端子２１３からは、仮想音源の距離ｒを表すインデックスＰと深度ｚを表すインデックスＱの設定値が入力される。
入力端子２１４からは、音速プロファイルが入力される。
これらのデータに基づき音線計算処理部７０３では音線計算を行う。音線計算の結果として、各到来波の到来時間の推測値（理論値）と到来ふ仰角の推測値（理論値）が求められる。
【００３９】
ふ仰角判定部７０４では、音線計算処理部７０３における計算の結果に基づいて、各仮想音源位置からのマルチパス波の各々について到来ふ仰角が上向きか下向きかの判定を行い、判定結果を到来時間差算出部７０５に供給する。
【００４０】
到来時間差算出部７０５は、音線計算処理部７０３における計算の結果に基づき、さらにふ仰角判定部７０４における判定結果に基づき、上向きと判定されたマルチパス波については、上側の受波器２０１における到来時間差（直接波と当該マルチパス波の到来時間の差）を算出し、算出した到来時間差（推測値）を、記憶部２０９Ｍ内のふ仰上向き到来時間差テーブル７０６にテーブル値として書き込み、下向きと判断されたマルチパス波については、下側の受波器２０２における到来時間差（直接波と当該マルチパス波の到来時間の差）を算出し、算出した到来時間差（推測値）を、記憶部２０９Ｍ内のふ仰下向き到来時間差テーブル７０７にテーブル値として書き込む。
【００４１】
このような処理を行う結果、ふ仰上向き到来時間差テーブル７０６には、各仮想音源位置（ｐ，ｑ）から発し、上側の受波器２０１に上側から（上向きのふ仰角で）到来するマルチパス波の各々についての到来時間差を表すデータの組が、それぞれの仮想音源毎に格納され、ふ仰下向き到来時間差テーブル７０７には、各仮想音源位置（ｐ，ｑ）から発し、下側の受波器２０２に下側から（下向きのふ仰角で）到来するマルチパス波の各々についての到来時間差を表すデータの組が、それぞれの仮想音源毎に格納される。
これらのテーブル７０６、７０７に格納されたデータは、後述のように、端子７０８、７０９を介してコスト関数算出器２１０に供給される。
【００４２】
次に、コスト関数算出器２１０の動作を、図１２を参照してより詳細に説明する。
入力端子８０１には、第１の到来時刻決定部２０５の出力が供給され、入力端子８０２には、第２の到来時刻決定部２０６の出力が供給され、入力端子８０３には、第１のマルチパス波到来時間差算出部２０７の出力が供給され、入力端子８０４には、第２のマルチパス波到来時間差算出部２０８の出力が供給される。
【００４３】
ふ仰角判定部８０５は、第１及び第２の到来時刻決定部２０５、２０６の出力を、入力端子８０１、８０２を介して受け、受波器２０１、２０２で受信した到来波の各々についてふ仰角が上向きであるか下向きであるかを判定し、判定結果ＹＵＬを出力する。
【００４４】
図８（ａ）及び（ｂ）に示す例では、１番目、２番目、４番目の到来波については、上側の受波器２０１のほうがより早いタイミングで受信しており、３番目、５番目の到来波については、下側の受波器のほうがより早いタイミングで受信している。これらのことから、１番目、２番目、４番目の到来波は、上方から到来したこと、即ちふ仰角が上向きであること、３番目、５番目の到来波は下方から到来したこと、即ちふ仰角が下向きであることが分かる。
【００４５】
データ抽出部８０６は、ふ仰角判定部８０５における判定結果ＹＵＬに基づいて、第１のマルチパス波到来時間差算出部２０７から出力される到来時間差を表すデータのうち、ふ仰角判定部８０５でふ仰角が上向きと判定されたマルチパス波についての到来時間差のみを選択乃至抽出して、ふ仰上向きテーブル値一観測値差算出器８０８に供給する。
【００４６】
データ抽出部８０７は、ふ仰角判定部８０５における判定結果ＹＵＬに基づいて、第２のマルチパス波到来時間差算出部２０８から出力された到来時間差を表すデータのうち、ふ仰角判定部８０５でふ仰角が下向きと判定されたマルチパス波についての到来時間差のみを選択乃至抽出して、ふ仰下向きテーブル値一観測値差算出器８０９に供給する。
【００４７】
ふ仰角判定部８０５で上向きと判定されたマルチパス波についての、上側の受波器２０１における到来時間差を
ＦＵａ、（ａ＝１、２、…）で表し、
ふ仰角判定部８０５で下向きと判定されたマルチパス波についての、下側の受波器２０２における到来時間差を
ＦＬｂ、（ｂ＝１、２、…）で表す。
それぞれのマルチパスの伝搬経路が図１０に示す通り（即ち、推定通り）であれば、１番目、２番目、４番目の到来波（直接波と第１、第３のマルチパス波）が上向きのふ仰角を有し、３番目、５番目の到来波（第２、第４のマルチパス波）が下向きのふ仰角を有することになり、したがって、
ＦＵ１＝ＤＵ２、
ＦＵ２＝ＤＵ４、
ＦＬ１＝ＤＬ３、
ＦＬ２＝ＤＬ５
となる。
到来時間差ＦＵ１、ＦＵ２、…がふ仰上向きテーブル値一観測値差算出器８０８に供給され、到来時間差ＦＬ１、ＦＬ２、…がふ仰下向きテーブル値一観測値差算出器８０９に供給される。
【００４８】
なお、同じ到来波の受波器２０１、２０２における受信のタイミングの差は極めて小さいことから、例えば検出タイミングの差が所定値以下であれば、同じ到来波が受信されたと判断することとしても良い。また、第１の到来時刻決定部２０５で到来時刻と決定された時点及びその前後にわたるサンプル値の列と、第２の到来時刻決定部２０６で到来時刻と決定された時点及びその前後にわたるサンプル値の列の相関が高い場合に、これらの到来時刻が同じ到来波の検出タイミングであると判断することとしても良い。
【００４９】
ふ仰上向きテーブル値一観測値差算出器８０８は、第１のマルチパス波到来時間差算出部２０７から入力端子８０３を介して供給される到来時間差の観測値のうち、データ抽出部８０６で抽出されたふ仰上向きのマルチパス波についての到来時間差の観測値の列の各要素ＦＵ１、ＦＵ２、…と、ふ仰上向き到来時間差テーブル７０６から入力端子７０８を介して供給されるふ仰上向き到来時間差テーブル値ＥＵ１（ｐ，ｑ）、ＥＵ（ｐ，ｑ）２、…（ｐ＝１〜Ｐ，ｑ＝１〜Ｑ）の各々の各要素とを受け、両者の対応する要素間の差を算出する。即ち、ふ仰上向き到来時間差テーブル値のｉ番目の要素をＥＵｉ（ｐ，ｑ）、ふ仰上向きのマルチパス波についての到来時間差の観測値のｉ番目の要素をＦＵｉとすると、次の式で表される計算を行う。
△ＦＵｉ（ｐ，ｑ）＝ ＥＵｉ（ｐ，ｑ）−ＦＵｉ …（３）
【００５０】
ふ仰上向き自乗和算出器８１２では、式（３）の演算結果を自乗したものの、すべての、ふ仰角上向きのマルチパス波についての総和（下記の式（４）で表される）を求める。
【数２】
【００５１】
式（４）で、Ｉ’は、Ｉ’≦Ｉであり、ふ仰角が上向きのマルチパス波がＩ番目まで検出された場合には、Ｉ’＝Ｉであり、Ｉ番目まで検出されなかった場合には、検出されたふ仰角が上向きのマルチパス波の数を表す。
【００５２】
同様に、ふ仰下向きテーブル値一観測値差算出器８０９は、第２のマルチパス波到来時間差算出部２０８から入力端子８０２を介して供給される到来時間差の観測値のうち、データ抽出部８０７で抽出されたふ仰下向きのマルチパス波についての到来時間差の観測値の列の各要素ＦＬ１、ＦＬ２、…と、ふ仰下向き到来時間差テーブル７０７から入力端子７０９を介して供給されるふ仰下向き到来時間差テーブル値ＥＬ１（ｐ，ｑ）、ＥＬ２（ｐ，ｑ）、…（ｐ＝１〜Ｐ、ｑ＝１〜Ｑ）の各々の各要素とを受け、両者の対応する要素間の差を算出する。即ち、ふ仰下向き到来時間差テーブル値のｉ番目の要素をＥＬｉ（ｐ，ｑ）、ふ仰下向きのマルチパス波についての到来時間差の観測値のｉ番目の要素をＦＬｉとすると、次の式で表される計算を行う。
△ＦＬｊ（ｐ，ｑ）＝ ＥＬｊ（ｐ，ｑ）−ＦＬｊ …（５）
【００５３】
ふ仰下向き自乗和算出器８１３では、式（５）の演算結果を自乗したものの、すべての、ふ仰角下向きのマルチパス波についての総和（下記の式（６）で表される）を求める。
【数３】
【００５４】
式（６）で、Ｊ’は、Ｊ’≦Ｊであり、ふ仰角が下向きのマルチパス波がＪ番目まで検出された場合には、Ｊ’＝Ｊであり、Ｊ番目まで検出されなかった場合には、検出されたふ仰角が下向きのマルチパス波の数を表す。
【００５５】
和算出器８１４では、式（４）及び式（６）の演算結果を用いて、次の演算を行い、演算結果を接続端子８１５から出力する。
Ｊ（ｐ，ｑ）＝ＪＵ（ｐ，ｑ）＋ＪＬ（ｐ，ｑ） …（７）
式（７）のＪ（ｐ，ｑ）が本実施の形態１で用いるコスト関数である。
【００５６】
音源位置推定器２１１では式（７）で与えられるコスト関数が最小となる（ｐ，ｑ）を求め、求められた（ｐ，ｑ）が推定対象（目標）の音源位置ＷＳを表すものであると推定する。
【００５７】
式（４）、式（６）で求められる自乗和はともに、観測値と理論値（テーブル値）との一致度が大きいほど小さくなる。従って、式（４）、式（６）の自乗和の和で与えられる式（７）のコスト関数の値が小さいほど上記の一致度が大きい(不一致度が小さい)ことを意味する。従って式（７）のコスト関数が最小となる（ｐ，ｑ）を求めることは、一致度が最も大きい（不一致度が最も小さい）（ｐ，ｑ）を求めることを意味する。
【００５８】
以上のようにしてコスト関数を求め、該コスト関数が最小となる（ｐ，ｑ）を求めることで、音源位置の推定を正確に行うことができる。即ち、従来例のコスト関数は、図４中の矢印ＤＭＣで示した方向に沿う領域でコスト関数の差が小さく、従って、到来時間差観測値に誤差がある場合に、音源位置の推定誤差が大きいという問題があったが、本実施の形態１では、上記のように、マルチパス波の到来ふ仰角の上向き、下向きのそれぞれについて別個の到来時間差のテーブルを用い、別個の不一致度（を表す値）ＪＵ、ＪＬを求め、これらの和をコスト関数として、音源位置の推定を行なうこととしており、これにより推定誤差を低減することができる。
【００５９】
一例として、従来例と同様に、図３に示すように音速が一定で水深の変化が無い場合を考える。従来例と同条件で到来時間差観測値に誤差がないとした場合のコスト関数を、図１３に示す。従来例と同様に、横軸は受波点Ｚｒからの水平方向距離ｒ、縦軸は深度（下方向ほど深度が大）ｚを表す。図中の線はコスト関数の等高線である。図示のように、コスト関数の等高線は音源位置真値を中心に楕円状になり、受波器に近い場所はコストが大きくなるように変化している。従って、コスト関数の形状から定性的に改善していることが判る。
【００６０】
また、従来例と同様に、経験的に知られている観測誤差を与えて、図３に示した条件でシミュレーションを３０回試行したところ、距離の誤差標準偏差は距離真値の約５％と従来例と比較して約１／３となり、音源の位置推定誤差に低減に寄与することが確認された。
【００６１】
なお、上記の実施の形態１では、ふ仰角の判定のために２つの受波器における到来時刻の差を用いているが、２つの受波器の代わりに、例えば、鉛直受波器アレイを用いることも可能である。
【００６２】
実施の形態２．
上記の実施の形態１の方法では、各受波器で受信した音波が、信号のみから成る場合には問題ないが、ノイズが誤って信号と検出（誤検出）されたり、検出されるべき信号が検出されない場合には、問題がある。誤検出の場合には、誤った組合せで受波器間での到来時間差を計算し、誤ったふ仰角の上下判定の結果に基づいてコスト算出を行う可能性があり、コスト算出に誤りが生じ、音源位置推定を誤る可能性がある。検出されるべき信号が検出されない場合には、ふ仰角の判定のために用いるべき一対の信号（受信タイミング）が揃わず、そのため、到来時間差テーブルとの比較（不一致度の算出）が適切に行なえないという問題がある。
【００６３】
本発明の実施の形態２における音源位置推定装置を図１４に示す。図１４の音源位置推定装置は概して図５に示す音源位置推定装置と同じであるが、図５のコスト関数算出器２１０の代わりに、コスト関数算出器２２０が設けられている。
【００６４】
図１４のコスト関数算出器２２０は、例えば図１５に示すように構成されている。図１４及び図１５に示されるコスト関数算出器２２０は、図５及び図１２に示されるコスト関数算出器２１０と概して同じであるが、入力端子２１５を介して供給される受波器位置データ（受波器間距離ＤＷＲを表す）を受けるとともに、入力端子２１４を介して供給される音速プロファイルデータを受ける。また、図１５に示すように、ふ仰角判定部８０５の代わりに、ふ仰角判定部８２５を備えている。ふ仰角判定部８２５は、ふ仰角判定部８０５と同様の機能を有するほか、以下のような付加的な機能をも有する。
【００６５】
図１６は、ふ仰角判定部８２５の詳細を示す。図示のふ仰角判定部８２５は、到来時間差最小組合せ算出部８５３と、受波器間最大到来時間差算出部８５５と、観測値・最大値時間差比較部８５４と、データ抽出部８５８と、ふ仰角上下判定部８５９とを有する。
【００６６】
到来時間差最小組合せ算出８５３では、入力端子８０１、８０２を介して供給される第１及び第２の到来時刻決定部２０５、２０６の出力に基づき、２つの受波器２０１、２０２における到来波の到来時刻（タイミング）のうちで、差が最も小さい組合せを求める。その動作を、図１７（ａ）及び（ｂ）を用いて説明する。第１の到来時刻決定部２０５で検出された、上側の受波器２０１における到来時刻（タイミング）をｔＡｃ（ｃ＝１，２，…）、第２の到来時刻決定部２０６で検出された下側の受波器２０２における到来時刻（タイミング）をｔＢｄ（ｄ＝１，２，…）とし、２つの受波器における到来時刻の内、差が最小となる組合せを求める。
【００６７】
例えば、時刻ｔＡｃの各々について時刻ｔＢｄのうちで最も差の小さいものを求め、時刻ｔＢｄ（ｄ＝１、２、…）のうちで、時刻ｔＡｃ（ｃ＝１、２、…）のいずれについても差が最小と判断されなかったものは、除外される。
図示の例では、時刻ｔＢ３は、時刻ｔＡ１〜ｔＡ６のいずれについても、差が最小とはならないため、除外される。
【００６８】
一方、時刻ｔＡ４と時刻ｔＡ５のいずれも、時刻ｔＢ５との差が最小であるが、時刻ｔＡ５よりも時刻ｔＡ４の方が時刻ｔＢ５との差が小さいので、時刻ｔＡ４が残され（時刻ｔＢ５と対を成すものとして選択され）、時刻ｔＡ５は除外される。
【００６９】
以上の結果、時刻ｔＡ１と時刻ｔＢ１、時刻ｔＡ２と時刻ｔＢ２、時刻ｔＡ３と時刻ｔＢ４、時刻ｔＡ４と時刻ｔＢ５、時刻ｔＡ６と時刻ｔＢ６が互いに対を成すものとして（従って同じ到来波の検出タイミングとして）選択され、時刻ｔＡ５と時刻ｔＡ３はノイズによるものと見なされ除外される。
【００７０】
受波器間最大到来時間差算出部８５５では、入力端子２１５から入力した受波器間距離ＤＷＲと、入力端子２１４から入力した音速プロファイルの、受波器２０１、２０２の近傍における音速の値Ｃから、
ｔｍａｘ＝ＤＷＲ／Ｃ …（８）
により、同じ伝搬経路を辿った音波の、受波器２０１における到来のタイミングと、受波器２０２における到来のタイミングの時間差の理論的最大値ｔｍａｘを計算する。受波器２０１、２０２は鉛直方向に整列しているので、到来ふ仰角の絶対値が大きいほど上記の時間差は大きくなる。
【００７１】
観測値・最大値時間差比較部８５４では、組み合わせ算出部８５３で決めた到来のタイミング（時刻）の対の時間差（観測値）ｔｏｂｓが受波器間最大到来時間差ｔｍａｘを超えているか否か、即ち、
ｔｏｂｓ＞ｔｍａｘ
であるか否かを判定し、判定結果ＹＶＮを出力する。判定結果ＹＶＮは、データ抽出部８５８に供給されるほか、出力端子８６２を介して図１５のデータ抽出部８０６、８０７に供給される。判定結果ＹＶＮは、第１及び第２の受波器２０１、２０２における到来時刻の対が、同じ到来波の検出タイミングである蓋然性が高い（信頼性が高い）か否かを示し、ｔｏｂｓ≦ｔｍａｘであれば、上記蓋然性が高いことを意味し、従って、当該到来時刻のデータがふ仰角の判定、到来時間差の比較などに利用するのに適していることを意味し、逆にｔｏｂｓ＞ｔｍａｘであれば、上記蓋然性が低く、当該到来時刻のデータがふ仰角の判定、到来時間差の比較などに利用するのに適さないことを意味する。
【００７２】
データ抽出部８５８では、観測値・最大値時間差比較部８５４による判定結果ＹＶＮに基づいて、入力端子８０１、８０２を介して供給される第１及び第２の到来時刻決定部２０５、２０６の出力を、ふ仰角の判定に使用するかしないかを決定し、ふ仰角の判定に使用すると決定した時刻を表すデータ（到来時刻決定部２０５、２０６の出力）を抽出して、ふ仰角上下判定部８５９に供給する。即ち、各時刻対（互いに対をなすと判断された検出時刻）に関し、
ｔｏｂｓ＞ｔｍａｘ
の場合には、当該時刻対を、ふ仰角の判定に使用しないこととし、ふ仰角上下判定部８５９に供給しない。
【００７３】
一方、各時刻対に関し、
ｔｏｂｓ≦ｔｍａｘ
の場合には、当該時刻対を表すデータを抽出して、ふ仰角上下判定部８５９に供給する。
【００７４】
ふ仰角上下判定部８５９では、データ抽出部８５８から供給された時刻対を表すデータについて、ふ仰角が上向きか下向きかの判定を行い、判定結果ＹＵＬを、接続端子８６１を介して出力する。
【００７５】
データ抽出部８０６は、ふ仰角判定部８２５における判定結果（ふ仰角が上向きか下向きかの判定結果ＹＵＬ、及び各検出タイミングｔＡｃ，ｔＢｄの信頼性についての判定結果ＹＶＮ）に応じて、第１のマルチパス波到来時間差算出部２０７の出力（入力端子８０３を介して供給される）のうちの、データ抽出部８５８で抽出された到来時刻に対応し、かつふ仰角上下判定部８５９でふ仰角が上向きと判定された到来時間差データ値のみを抽出して、ふ仰上向きテーブル値一観測値差算出器８０８に供給する。
【００７６】
データ抽出部８０７は、ふ仰角判定部８０５における判定結果（ふ仰角が上向きか下向きかの判定結果ＹＵＬ、及び各検出タイミングｔＡｃ，ｔＢｄの信頼性についての判定結果ＹＶＮ）に応じて、第２のマルチパス波到来時間差算出部２０８の出力（入力端子８０４を介して供給される）のうちの、データ抽出部８５８で抽出された到来時刻に対応し、かつふ仰角上下判定部８５９でふ仰角が下向きと判定された到来時間差データ値のみを抽出して、ふ仰下向きテーブル値一観測値差算出器８０９に供給する。
【００７７】
図１５のふ仰上向きテーブル値一観測値差算出器８０８及びふ仰下向きテーブル値一観測値差算出器８０９では、データ抽出部８０６及び８０７で抽出された到来時間差の観測値を用いて、差分の算出及びそれに基づくコスト関数の算出を行なう。コスト関数算出器２２０の他の部分の動作は実施の形態１について説明したのと同様である。
【００７８】
実施の形態１では雑音の影響の無いほぼ理想的な状況にのみ到来音波のふ仰角の上下を判別したコスト関数を利用することができたが、本実施の形態２で説明したように到来音波のふ仰角の判定を行う（信頼性の高いデータのみを用いてふ仰角の判定、到来時間差の比較（不一致度の検出）を行う）ことにより、ＳＮ比がさほど高く無く到来時間差に誤検出や未検出がある場合にも、到来音波のふ仰上下を判別したコスト関数を適用することが可能となる。
【００７９】
このため、ＳＮ比がさほど高く無く到来時間差に誤検出や未検出がある場合にもコスト関数を円錐状に保つことができ、音源の位置推定精度を保つことができる。
【００８０】
尚、実施の形態２では、ふ仰角の計測に２つの受波器の時間差を用いているが、３つ以上の受波器を用いて、ふ仰角の上下の判定を行うこととしても良い。
【００８１】
また、音速プロファイルを測定できなかった場合については、式（８）の音速プロファイルの受波器２０１，２０２の近傍における音速の値Ｃの代わりに代表的な音速１５００ｍ／ｓを用いて、
ｔｍａｘ＝ＤＷＲ／１５００ …（９）
を利用することも可能である。
【符号の説明】
【００８２】
２０１、２０２ 受波器、 ２０５、２０６ 到来時刻決定部、 ２０７、２０８ マルチパス波到来時間差算出部、 ２０９ 到来時間差テーブル作成器、 ２１０ コスト関数算出器、 ２１１ 音源位置推定器。
【技術分野】
【０００１】
この発明は音源位置推定装置に関し、特に水中の音源から放射されたパルス状の音波の直接波とマルチパス波を受信し、これらの到来時間差から音源の位置を推定する装置に関する。
【背景技術】
【０００２】
水中の音波伝搬では、音源から直接受波器に到来する音波（直接波）の他に、海面や海底で反射して到来するマルチパス波があり、そのようなマルチパス波は、直接波と異なる時刻に到来する。
水中において音源からのパルス状の音波を受信して音源の位置を推定する技術として、下記の特許文献１に示されるものがある。
【０００３】
図１に、特許文献１に示される、従来の音源位置推定装置を示す。図示の音源位置推定装置においては、音響センサー（受波器）で受信した音波を表す信号が入力端子１０１を介してマルチパス波到来時間差測定器１０２に入力される。マルチパス波到来時間差測定器１０２では、直接波とマルチパス波の到来時間差を検出し、検出した到来時間差観測値を示すデータを出力する。
【０００４】
一方、到来時間差テーブル作成器１０３により、水中の音速プロファイルに基づき、各仮想音源からの音波の音線計算を行うことにより音波の到来時間を推定し、推定結果を到来時間差テーブル値（理論値）として、記憶しておく。
【０００５】
コスト関数算出器１０４は、マルチパス波到来時間差測定器１０２で求めたマルチパス波到来時間差の観測値と到来時間差テーブル作成器１０３で求めた到来時間差テーブル値より、それらの不一致度を反映したコスト関数の値を算出する。
【０００６】
コスト関数算出器１０４の詳細ブロック構成を図２に示す。
入力端子１１１にはマルチパス波到来時間差測定器１０２で求めたマルチパス波到来時間差観測値が入力され、入力端子１１２には到来時間差テーブル作成器１０３からのテーブル値が入力される。テーブル値−観測値算出器１１３は、到来時間差観測値と到来時間差テーブル値との差を算出する。
【０００７】
ここで、ｐを受波位置（受波点）から音源位置までの距離に対応するインデックス、ｑを音源位置の深度に対応するインデックスとし、マルチパス波到来時間差測定器１０２で得られるｎ番目の到来波についての到来時間差をＤｎ、到来時間差テーブル作成器１０３で求めた、仮想音源位置（ｐ,ｑ）の仮想音源からの、ｎ番目の到来時間差テーブル値をＥｎ（ｐ，ｑ）とした場合、テーブル値−観測値差算出器２０３では、テーブル値と観測値の差を次の式（１）で求める。
△Ｄｎ（ｐ，ｑ）＝ Ｅｎ（ｐ，ｑ）−Ｄｎ …（１）
【０００８】
自乗和算出器１１４は、式（１）の演算結果を自乗したものの、すべてのマルチパス波についての総和（下記の式（２）で表される）を、コスト関数として求め、求めたコスト関数を出力端子１１５から出力する。
【０００９】
【数１】
【００１０】
図１の音源位置推定器１０５では式（２）で与えられるコスト関数が最小となる（ｐ，ｑ）を求め、求められた（ｐ，ｑ）が推定対象（目標）の音源位置を表すものであると推定し、推定結果を、出力端子１０６を介して出力する。
【先行技術文献】
【特許文献】
【００１１】
【特許文献１】特開２００６−１９４６２７号公報
【発明の概要】
【発明が解決しようとする課題】
【００１２】
しかし、上記の従来の装置を用いた場合、音源と受波器を結ぶ線上のコスト関数の値が、音源の位置の前後のかなり広い範囲にわたり小さな値となるため、音源位置を正確に推定することができないという問題があった。即ち、マルチパス波到来時間差測定器１０２で求められる到来時間差は観測誤差を含むため、コスト関数の値が小さい場所で大きく推定値がずれる可能性が大きく、このため観測誤差が比較的大きいという問題があった。
【００１３】
一例として、図３に示すように音速が一定で水深Ｚｂが一定である場合を考える。図３において、Ｚｓは音源の深度、Ｚｒは受波点の深度、ｒは受波点ＷＲから音源ＷＳまでの水平方向距離を示す。
到来時間差観測値に誤差がないとした場合のコスト関数を求めた結果を図４に示す。横軸は受波点Ｚｒからの水平方向距離ｒ、縦軸は深度（下方向ほど深度が大）ｚを表す。図中の曲線はコスト関数の等高線であり、このコスト関数が最小となる位置を推定する。
【００１４】
到来時間差観測値に誤差がない場合は、真値を推定することができるが、図中のコスト関数は、受波器と音源とを結ぶ線の方向（矢印ＤＭＣの方向）に沿って、コスト関数の変化が小さいために、到来時間差の観測値に誤差が含まれる場合は、この方向に沿って、真の音源位置から推定結果がずれる可能性がある。経験的に知られている観測誤差を与えて、図３に示した条件でシミュレーションを３０回試行したところ、距離の誤差標準偏差は、距離真値の約１５％と、かなり大きい値となった。このようにコスト関数の特質から音源の位置推定誤差が大きいという問題があった。
【課題を解決するための手段】
【００１５】
本発明の音源推定装置は、
受波位置における到来波に基づいて音源の位置を推定する音源推定装置において、
複数の仮想音源位置の各々について、該仮想音源位置に仮想音源が存在すると仮定したときの、前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が上向きのものの各々との到来時間差の理論値、並びに前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が下向きのものの各々との到来時間差の理論値を格納する到来時間差理論値格納手段と、
前記受波位置において観測した到来波のデータに基づいて、目標の音源からの最初の到来波と２番目以降の到来波との到来時間差の観測値を算出する到来時間差観測値算出手段と、
前記２番目以降の到来波の各々の到来ふ仰角が上向きか下向きかを判定し、前記２番目以降の到来波のうち、到来ふ仰角が上向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が上向きの到来波についての前記到来時間差の理論値との不一致度を第１の不一致度として求め、前記２番目以降の到来波のうち、到来ふ仰角が下向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が下向きの到来波についての前記到来時間差の理論値との不一致度を第２の不一致度として求め、前記第１の不一致度と前記第２の不一致度との和をコストとして求めるコスト関数算出手段と、
前記コスト関数算出手段が算出した前記コストが最小となる前記仮想音源位置を前記目標の音源の位置と推定する音源位置推定手段と
を備えたことを特徴とする。
【発明の効果】
【００１６】
本発明によれば、到来時間差観測値に誤差がある場合にも、音源位置の推定誤差を小さくすることができる。
【図面の簡単な説明】
【００１７】
【図１】従来の音源位置推定装置を示すブロック図である。
【図２】図１のコスト関数算出器１０４の詳細を示すブロック図である。
【図３】音速が一定で水深が一定である場合の音波の伝搬経路を示す図である。
【図４】図１の装置で求めたコスト関数の等高線の一例を示す図である。
【図５】本発明の実施の形態１による音源位置推定装置を示すブロック図である。
【図６】実施の形態１における受波器の配置を示す図である。
【図７】受波器におけるサンプル値の絶対値の包絡線の一例を示すグラフである。
【図８】（ａ）及び（ｂ）は、一対の受波器における受信のタイミングの一例を示す図である。
【図９】到来時間差テーブルの作成に当たり使用される仮想音源位置の一例を示す図である。
【図１０】仮想音源位置ＷＳの音源からの伝搬経路の一例を示す図である。
【図１１】図５の到来時間差テーブル作成器２０９の一例を示す機能ブロック図である。
【図１２】図５のコスト関数算出器２１０の一例を示す機能ブロック図である。
【図１３】図５の装置で求めたコスト関数の等高線の一例を示す図である。
【図１４】本発明の実施の形態２による音源位置推定装置を示すブロック図である。
【図１５】図１４のコスト関数算出器２２０の一例を示す機能ブロック図である。
【図１６】図１５のふ仰角判定部８２５の一例を示す機能ブロック図である。
【図１７】（ａ）及び（ｂ）は、一対の受波器における受信のタイミングの一例を示す図である。
【発明を実施するための形態】
【００１８】
実施の形態１．
本発明の実施の形態１の音源位置推定装置を図５に示す。
図示の音源位置推定装置は、音源から放射されたパルス状音波の直接波とマルチパス波を受信し、それらの到来時間差に基づいて音源の位置を推定するものであり、音源の推定位置は深度及び受波器からの水平方向距離で表される。なお、音源の位置する方位を示す情報は別途得られるものとする。
【００１９】
実施の形態１の音源位置推定装置は、第１及び第２の受波器２０１及び２０２を備えている。
第１の受波器２０１及び第２の受波器２０２は、後述のように、到来波のふ仰角が上向きか下向きかの情報を得るために設けられたものであり、同じ鉛直線上の異なる位置に設けられている。このように配置されることにより、第１の受波器２０１及び第２の受波器２０２のいずれに、到来波（信号）がより早く到来したかに基づいて、到来波のふ仰角が上向きか下向きかの判断を行うことができる。
【００２０】
第１及び第２の受波器２０１及び２０２は、ふ仰角が上向きか下向きの判別が行えれば良いので、図６に示すように、間隔ＤＷＲを開けて配置されている。間隔ＤＷＲは受信位置において想定される音速の最大値と受波器２０１、２０２で測定可能な最小時間差の積以上の値に設定される。
第１及び第２の受波器２０１及び２０２はそれぞれ上側及び下側に配置されるので、「上側の受波器」、「下側の受波器」と呼ばれることもある。
【００２１】
上側及び下側の受波器２０１及び２０２の出力は、入力端子２０３及び２０４を介してそれぞれ第１及び第２の到来時刻決定部２０５及び２０６に入力される。
第１及び第２の到来時刻決定部２０５及び２０６は、それぞれ上側及び下側の受波器２０１及び２０１の出力に基づいて到来時刻を決定乃至検出する。
【００２２】
第１及び第２の受波器２０１及び２０２は、それらで受信される音波をサンプリングすることで得られるサンプル値を表すデータを出力するが、音波をサンプリングすることで得られるサンプル値の絶対値を時間軸上に並べたものの包絡線は、例えば、図７に示す如くとなる。なお、サンプル値の絶対値の代りにサンプル値の自乗値の包絡線を用いる場合もある。
【００２３】
第１及び第２の到来時刻決定部２０５及び２０６では、それぞれ上側及び下側の受波器２０１及び２０２からのデータの列に基づき到来波の受信時刻（到来時刻）を判断する。到来時刻の判断には、従来と同様に、レプリカ相関による方法、立ち上がり検出による方法などを用いることができる。これらの方法は、例えば上記特許文献１においてマルチパス波到来時間差測定器の動作の一部として説明されている。
【００２４】
図８（ａ）及び（ｂ）に、第１及び第２の到来時刻決定部２０５及び２０６で決定乃至検出された上側の受波器２０１及び下側の受波器２０２におけるパルス状波の到来時刻（検出タイミング）を示す。図８（ａ）は、上側の受波器２０１における検出タイミングを、早い順にＴＡ１、ＴＡ２、…で示し、図８（ｂ）は、下側の受波器２０２における検出タイミングを、早い順にＴＢ１、ＴＢ２、…で示す。
ノイズを信号と誤って検出したり、検出漏れが発生しなければ、最初の検出タイミングが直接波の検出タイミングであり、２番目以降の検出タイミングがマルチパス波の検出タイミングである。
【００２５】
第１のマルチパス波到来時間差算出部２０７は、第１の到来時刻決定部２０５の出力に基づいて、上側の受波器２０１で受信したマルチパス波の各々についての到来時間差（１番目の到来波との到来時間差）を算出する。即ち、１番目の到来波ＰＡ１と２番目の到来波ＰＡ２（１番目のマルチパス波）の到来時間差ＤＵ２、１番目の到来波ＰＡ１と３番目の到来波ＰＡ３の到来時間差ＤＵ３、一般化すれば、１番目の到来波ＰＡ１とｎ番目の到来波ＰＡｎ（ｎ＝２、３、…）との到来時間差ＤＵｎを算出する。このようにして算出された到来時間差ＤＵｎを示すデータの列[ＤＵ２、ＤＵ３、…］はコスト関数算出器２１０に供給される。
【００２６】
同様に、第２のマルチパス波到来時間差算出部２０８は、第２の到来時刻決定部２０６の出力に基づいて、下側の受波器２０２で受信したマルチパス波の各々についての到来時間差（１番目の到来波との到来時間差）を算出する。即ち、１番目の到来波ＰＢ１と２番目の到来波ＰＢ２（１番目のマルチパス波）の到来時間差ＤＬ２、１番目の到来波ＰＢ１と３番目の到来波ＰＢ３の到来時間差ＤＬ３、一般化すれば、１番目の到来波ＰＢ１とｎ番目の到来波ＰＢｎ（ｎ＝２、３、…）との到来時間差ＤＬｎを算出する。このようにして算出された到来時間差ＤＬｎを示すデータの列[ＤＬ２、ＤＬ３、…］はコスト関数算出器２１０に供給される。
【００２７】
第１の到来時刻決定部２０５と第１のマルチパス波到来時間差算出部２０７の組合せは、図１のマルチパス波到来時間差測定器１０２と同様の機能を有し、第２の到来時刻決定部２０６と第２のマルチパス波到来時間差算出部２０８の組合せも、図１のマルチパス波到来時間差測定器１０２と同様の機能を有する。
【００２８】
上記のように、第１のマルチパス波到来時間差算出部２０７と第２のマルチパス波到来時間差算出部２０８とは、受波位置において観測した到来波の到来時刻（検出タイミング）に基づいて、推定対象（目標）の音源からの最初の到来波と２番目以降の到来波の各々との到来時間差の観測値を算出する到来時間差観測値算出手段を構成し、そのうち、第１のマルチパス波到来時間差算出部２０７が、第１の受波器２０１への最初の到来波と２番目以降の到来波の各々との到来時間差の観測値を算出し、第２のマルチパス波到来時間差算出部２０８が、第２の受波器２０２への最初の到来波と２番目以降の到来波の各々との到来時間差の観測値を算出する。
【００２９】
到来時間差テーブル作成器２０９は、入力端子２１３を介して供給される仮想音源の設定位置データと入力端子２１４を介して供給される音速プロファイルデータに基づいて、到来時間差テーブルを作成し、記憶部２０９Ｍに保存する。
【００３０】
到来時間差テーブルを作成する際の仮想音源位置は例えば、図９に示すように設定される。
図９に示す例では、探索領域を仮想的にメッシュ（格子）状に区切り、各交点を仮想音源位置（ｐ，ｑ）（ｐ＝１，２，…Ｐ、ｑ＝１，２，…Ｑ）とする。なお、ｐ、ｑはそれぞれ距離ｒ、深度ｚを表す値そのものではなく、距離、深度に対応するインデックスであり、自然数１〜Ｐ、１〜Ｑで表される。
【００３１】
メッシュ幅（仮想音源位置の間隔）は、所望の音源位置の推定精度に応じて設定される。メッシュ幅を狭くするほど高い精度で測定を行なうことができるが、データ量、演算量が多くなり、リアルタイムでの処理が困難となる。従って、高精度が要求される場合には、メッシュ幅を狭くし、高速性が重視される場合には、メッシュ幅を広くする。
【００３２】
仮想音源位置ＷＳからの伝搬経路は例えば図１０に示すように想定される。図１０には、簡単のため、５つの伝搬経路ＷＰ１〜ＷＰ５のみが示されている。伝搬経路（音線：音のエネルギーの流れ）が直線的でないのは、水温、塩分濃度などの不均一性のため、音速が一定ではないことによる。
一度も反射をせずに受波点ＷＲに達した波は直接波と呼ばれ、海面や海底で１回以上反射して受波器に達した波はマルチパス波と呼ばれる。
【００３３】
このような想定に基づいて、仮想音源位置の各々からの到来時間差を推測し、推測値（理論値）を時間差テーブル値として求め、記憶部２０９Ｍに記憶しておく。到来時間差テーブル作成器２０９は、仮想音源位置の各々でパルス状音波が発せられたと仮定した場合に、上側の受波器２０１における第１の到来波（直接波）の推定到来時刻と、上側の受波器２０１に上向きのふ仰角でｉ番目（ｉ＝１、２…）に到来すると推定されるマルチパス波の到来時刻との差（到来時間差）の推測値乃至理論値ＥＵｉ（ｉ＝１、２、…）の組乃至列[ＥＵ１、ＥＵ２、…ＥＵＩ]を、すべての仮想音源について求め、記憶部２０９Ｍに格納し、同様に、下側の受波器２０２における第１の到来波（直接波）の推定到来時刻と、下側の受波器２０２に下向きのふ仰角でｊ番目（ｊ＝１、２…）に到来すると推定されるマルチパス波の到来時刻との差（到来時間差）の推測値乃至理論値ＥＬｊ（ｊ＝１、２、…）の組乃至列[ＥＬ１、ＥＬ２、…ＥＬＪ]を、すべての仮想音源について求め、記憶部２０９Ｍに格納する。このように、記憶部２０９Ｍは、到来時間差理論値格納手段としての役割を果たす。
【００３４】
なお、Ｉは、経験的に受波器２０１で受信可能と推定される、ふ仰角が上向きのマルチパス波の最大数以上の値に、Ｊは、経験的に受波器２０２で受信可能と推定される、ふ仰角が下向きのマルチパス波の最大数以上の値に、それぞれ定められる。
【００３５】
全ての仮想音源位置（ｐ，ｑ）（ｐ＝１、２、…Ｐ、ｑ＝１、２、…Ｑ）の各々についての、上側の受波器２０１での到来時間差の列
ＥＵ１（ｐ，ｑ）、ＥＵ２（ｐ，ｑ）、…（ｐ＝１、２、…Ｐ、ｑ＝１、２、…Ｑ）
が得られ、これらがテーブル形式のデータとされ、上向きふ仰角到来時間差テーブル値として到来時間差テーブル作成器２０９の記憶部２０９Ｍ内に記憶されている。
同様に、全ての仮想音源位置（ｐ，ｑ）の各々についての、下側の受波器２０２での到来時間差の列
ＥＬ１（ｐ，ｑ）、ＥＬ２（ｐ，ｑ）、…（ｐ＝１、２、…Ｐ、ｑ＝１、２、…Ｑ）
が得られ、これらがテーブル形式のデータとされ、下向きふ仰角到来時間差テーブル値として到来時間差テーブル作成器２０９の記憶部２０９Ｍ内に記憶されている。
【００３６】
コスト関数算出器２１０は、第１及び第２の到来時刻決定部２０５、２０６の出力（到来時刻を示す信号）を入力とするとともに、第１及び第２のマルチパス波到来時間差算出部２０７及び２０８で得られた到来時間差データ（観測値）と、到来時間差テーブル作成器２１０で作成された到来時間差テーブルの値（理論値）を入力とし、これらに基づいてコスト関数を求め、求めたコスト関数を音源位置推定器２１１に入力する。
【００３７】
音源位置推定器２１１では、コスト間数算出器２１０で作成したコスト関数が最小となる仮想音源の距離及び深度を音源位置の推定値とし、該推定値を出力端子２１２から出力する。
【００３８】
以下、到来時間差テーブル作成器２０９の動作を、図１１を参照してより詳細に説明する。
入力端子２１３からは、仮想音源の距離ｒを表すインデックスＰと深度ｚを表すインデックスＱの設定値が入力される。
入力端子２１４からは、音速プロファイルが入力される。
これらのデータに基づき音線計算処理部７０３では音線計算を行う。音線計算の結果として、各到来波の到来時間の推測値（理論値）と到来ふ仰角の推測値（理論値）が求められる。
【００３９】
ふ仰角判定部７０４では、音線計算処理部７０３における計算の結果に基づいて、各仮想音源位置からのマルチパス波の各々について到来ふ仰角が上向きか下向きかの判定を行い、判定結果を到来時間差算出部７０５に供給する。
【００４０】
到来時間差算出部７０５は、音線計算処理部７０３における計算の結果に基づき、さらにふ仰角判定部７０４における判定結果に基づき、上向きと判定されたマルチパス波については、上側の受波器２０１における到来時間差（直接波と当該マルチパス波の到来時間の差）を算出し、算出した到来時間差（推測値）を、記憶部２０９Ｍ内のふ仰上向き到来時間差テーブル７０６にテーブル値として書き込み、下向きと判断されたマルチパス波については、下側の受波器２０２における到来時間差（直接波と当該マルチパス波の到来時間の差）を算出し、算出した到来時間差（推測値）を、記憶部２０９Ｍ内のふ仰下向き到来時間差テーブル７０７にテーブル値として書き込む。
【００４１】
このような処理を行う結果、ふ仰上向き到来時間差テーブル７０６には、各仮想音源位置（ｐ，ｑ）から発し、上側の受波器２０１に上側から（上向きのふ仰角で）到来するマルチパス波の各々についての到来時間差を表すデータの組が、それぞれの仮想音源毎に格納され、ふ仰下向き到来時間差テーブル７０７には、各仮想音源位置（ｐ，ｑ）から発し、下側の受波器２０２に下側から（下向きのふ仰角で）到来するマルチパス波の各々についての到来時間差を表すデータの組が、それぞれの仮想音源毎に格納される。
これらのテーブル７０６、７０７に格納されたデータは、後述のように、端子７０８、７０９を介してコスト関数算出器２１０に供給される。
【００４２】
次に、コスト関数算出器２１０の動作を、図１２を参照してより詳細に説明する。
入力端子８０１には、第１の到来時刻決定部２０５の出力が供給され、入力端子８０２には、第２の到来時刻決定部２０６の出力が供給され、入力端子８０３には、第１のマルチパス波到来時間差算出部２０７の出力が供給され、入力端子８０４には、第２のマルチパス波到来時間差算出部２０８の出力が供給される。
【００４３】
ふ仰角判定部８０５は、第１及び第２の到来時刻決定部２０５、２０６の出力を、入力端子８０１、８０２を介して受け、受波器２０１、２０２で受信した到来波の各々についてふ仰角が上向きであるか下向きであるかを判定し、判定結果ＹＵＬを出力する。
【００４４】
図８（ａ）及び（ｂ）に示す例では、１番目、２番目、４番目の到来波については、上側の受波器２０１のほうがより早いタイミングで受信しており、３番目、５番目の到来波については、下側の受波器のほうがより早いタイミングで受信している。これらのことから、１番目、２番目、４番目の到来波は、上方から到来したこと、即ちふ仰角が上向きであること、３番目、５番目の到来波は下方から到来したこと、即ちふ仰角が下向きであることが分かる。
【００４５】
データ抽出部８０６は、ふ仰角判定部８０５における判定結果ＹＵＬに基づいて、第１のマルチパス波到来時間差算出部２０７から出力される到来時間差を表すデータのうち、ふ仰角判定部８０５でふ仰角が上向きと判定されたマルチパス波についての到来時間差のみを選択乃至抽出して、ふ仰上向きテーブル値一観測値差算出器８０８に供給する。
【００４６】
データ抽出部８０７は、ふ仰角判定部８０５における判定結果ＹＵＬに基づいて、第２のマルチパス波到来時間差算出部２０８から出力された到来時間差を表すデータのうち、ふ仰角判定部８０５でふ仰角が下向きと判定されたマルチパス波についての到来時間差のみを選択乃至抽出して、ふ仰下向きテーブル値一観測値差算出器８０９に供給する。
【００４７】
ふ仰角判定部８０５で上向きと判定されたマルチパス波についての、上側の受波器２０１における到来時間差を
ＦＵａ、（ａ＝１、２、…）で表し、
ふ仰角判定部８０５で下向きと判定されたマルチパス波についての、下側の受波器２０２における到来時間差を
ＦＬｂ、（ｂ＝１、２、…）で表す。
それぞれのマルチパスの伝搬経路が図１０に示す通り（即ち、推定通り）であれば、１番目、２番目、４番目の到来波（直接波と第１、第３のマルチパス波）が上向きのふ仰角を有し、３番目、５番目の到来波（第２、第４のマルチパス波）が下向きのふ仰角を有することになり、したがって、
ＦＵ１＝ＤＵ２、
ＦＵ２＝ＤＵ４、
ＦＬ１＝ＤＬ３、
ＦＬ２＝ＤＬ５
となる。
到来時間差ＦＵ１、ＦＵ２、…がふ仰上向きテーブル値一観測値差算出器８０８に供給され、到来時間差ＦＬ１、ＦＬ２、…がふ仰下向きテーブル値一観測値差算出器８０９に供給される。
【００４８】
なお、同じ到来波の受波器２０１、２０２における受信のタイミングの差は極めて小さいことから、例えば検出タイミングの差が所定値以下であれば、同じ到来波が受信されたと判断することとしても良い。また、第１の到来時刻決定部２０５で到来時刻と決定された時点及びその前後にわたるサンプル値の列と、第２の到来時刻決定部２０６で到来時刻と決定された時点及びその前後にわたるサンプル値の列の相関が高い場合に、これらの到来時刻が同じ到来波の検出タイミングであると判断することとしても良い。
【００４９】
ふ仰上向きテーブル値一観測値差算出器８０８は、第１のマルチパス波到来時間差算出部２０７から入力端子８０３を介して供給される到来時間差の観測値のうち、データ抽出部８０６で抽出されたふ仰上向きのマルチパス波についての到来時間差の観測値の列の各要素ＦＵ１、ＦＵ２、…と、ふ仰上向き到来時間差テーブル７０６から入力端子７０８を介して供給されるふ仰上向き到来時間差テーブル値ＥＵ１（ｐ，ｑ）、ＥＵ（ｐ，ｑ）２、…（ｐ＝１〜Ｐ，ｑ＝１〜Ｑ）の各々の各要素とを受け、両者の対応する要素間の差を算出する。即ち、ふ仰上向き到来時間差テーブル値のｉ番目の要素をＥＵｉ（ｐ，ｑ）、ふ仰上向きのマルチパス波についての到来時間差の観測値のｉ番目の要素をＦＵｉとすると、次の式で表される計算を行う。
△ＦＵｉ（ｐ，ｑ）＝ ＥＵｉ（ｐ，ｑ）−ＦＵｉ …（３）
【００５０】
ふ仰上向き自乗和算出器８１２では、式（３）の演算結果を自乗したものの、すべての、ふ仰角上向きのマルチパス波についての総和（下記の式（４）で表される）を求める。
【数２】
【００５１】
式（４）で、Ｉ’は、Ｉ’≦Ｉであり、ふ仰角が上向きのマルチパス波がＩ番目まで検出された場合には、Ｉ’＝Ｉであり、Ｉ番目まで検出されなかった場合には、検出されたふ仰角が上向きのマルチパス波の数を表す。
【００５２】
同様に、ふ仰下向きテーブル値一観測値差算出器８０９は、第２のマルチパス波到来時間差算出部２０８から入力端子８０２を介して供給される到来時間差の観測値のうち、データ抽出部８０７で抽出されたふ仰下向きのマルチパス波についての到来時間差の観測値の列の各要素ＦＬ１、ＦＬ２、…と、ふ仰下向き到来時間差テーブル７０７から入力端子７０９を介して供給されるふ仰下向き到来時間差テーブル値ＥＬ１（ｐ，ｑ）、ＥＬ２（ｐ，ｑ）、…（ｐ＝１〜Ｐ、ｑ＝１〜Ｑ）の各々の各要素とを受け、両者の対応する要素間の差を算出する。即ち、ふ仰下向き到来時間差テーブル値のｉ番目の要素をＥＬｉ（ｐ，ｑ）、ふ仰下向きのマルチパス波についての到来時間差の観測値のｉ番目の要素をＦＬｉとすると、次の式で表される計算を行う。
△ＦＬｊ（ｐ，ｑ）＝ ＥＬｊ（ｐ，ｑ）−ＦＬｊ …（５）
【００５３】
ふ仰下向き自乗和算出器８１３では、式（５）の演算結果を自乗したものの、すべての、ふ仰角下向きのマルチパス波についての総和（下記の式（６）で表される）を求める。
【数３】
【００５４】
式（６）で、Ｊ’は、Ｊ’≦Ｊであり、ふ仰角が下向きのマルチパス波がＪ番目まで検出された場合には、Ｊ’＝Ｊであり、Ｊ番目まで検出されなかった場合には、検出されたふ仰角が下向きのマルチパス波の数を表す。
【００５５】
和算出器８１４では、式（４）及び式（６）の演算結果を用いて、次の演算を行い、演算結果を接続端子８１５から出力する。
Ｊ（ｐ，ｑ）＝ＪＵ（ｐ，ｑ）＋ＪＬ（ｐ，ｑ） …（７）
式（７）のＪ（ｐ，ｑ）が本実施の形態１で用いるコスト関数である。
【００５６】
音源位置推定器２１１では式（７）で与えられるコスト関数が最小となる（ｐ，ｑ）を求め、求められた（ｐ，ｑ）が推定対象（目標）の音源位置ＷＳを表すものであると推定する。
【００５７】
式（４）、式（６）で求められる自乗和はともに、観測値と理論値（テーブル値）との一致度が大きいほど小さくなる。従って、式（４）、式（６）の自乗和の和で与えられる式（７）のコスト関数の値が小さいほど上記の一致度が大きい(不一致度が小さい)ことを意味する。従って式（７）のコスト関数が最小となる（ｐ，ｑ）を求めることは、一致度が最も大きい（不一致度が最も小さい）（ｐ，ｑ）を求めることを意味する。
【００５８】
以上のようにしてコスト関数を求め、該コスト関数が最小となる（ｐ，ｑ）を求めることで、音源位置の推定を正確に行うことができる。即ち、従来例のコスト関数は、図４中の矢印ＤＭＣで示した方向に沿う領域でコスト関数の差が小さく、従って、到来時間差観測値に誤差がある場合に、音源位置の推定誤差が大きいという問題があったが、本実施の形態１では、上記のように、マルチパス波の到来ふ仰角の上向き、下向きのそれぞれについて別個の到来時間差のテーブルを用い、別個の不一致度（を表す値）ＪＵ、ＪＬを求め、これらの和をコスト関数として、音源位置の推定を行なうこととしており、これにより推定誤差を低減することができる。
【００５９】
一例として、従来例と同様に、図３に示すように音速が一定で水深の変化が無い場合を考える。従来例と同条件で到来時間差観測値に誤差がないとした場合のコスト関数を、図１３に示す。従来例と同様に、横軸は受波点Ｚｒからの水平方向距離ｒ、縦軸は深度（下方向ほど深度が大）ｚを表す。図中の線はコスト関数の等高線である。図示のように、コスト関数の等高線は音源位置真値を中心に楕円状になり、受波器に近い場所はコストが大きくなるように変化している。従って、コスト関数の形状から定性的に改善していることが判る。
【００６０】
また、従来例と同様に、経験的に知られている観測誤差を与えて、図３に示した条件でシミュレーションを３０回試行したところ、距離の誤差標準偏差は距離真値の約５％と従来例と比較して約１／３となり、音源の位置推定誤差に低減に寄与することが確認された。
【００６１】
なお、上記の実施の形態１では、ふ仰角の判定のために２つの受波器における到来時刻の差を用いているが、２つの受波器の代わりに、例えば、鉛直受波器アレイを用いることも可能である。
【００６２】
実施の形態２．
上記の実施の形態１の方法では、各受波器で受信した音波が、信号のみから成る場合には問題ないが、ノイズが誤って信号と検出（誤検出）されたり、検出されるべき信号が検出されない場合には、問題がある。誤検出の場合には、誤った組合せで受波器間での到来時間差を計算し、誤ったふ仰角の上下判定の結果に基づいてコスト算出を行う可能性があり、コスト算出に誤りが生じ、音源位置推定を誤る可能性がある。検出されるべき信号が検出されない場合には、ふ仰角の判定のために用いるべき一対の信号（受信タイミング）が揃わず、そのため、到来時間差テーブルとの比較（不一致度の算出）が適切に行なえないという問題がある。
【００６３】
本発明の実施の形態２における音源位置推定装置を図１４に示す。図１４の音源位置推定装置は概して図５に示す音源位置推定装置と同じであるが、図５のコスト関数算出器２１０の代わりに、コスト関数算出器２２０が設けられている。
【００６４】
図１４のコスト関数算出器２２０は、例えば図１５に示すように構成されている。図１４及び図１５に示されるコスト関数算出器２２０は、図５及び図１２に示されるコスト関数算出器２１０と概して同じであるが、入力端子２１５を介して供給される受波器位置データ（受波器間距離ＤＷＲを表す）を受けるとともに、入力端子２１４を介して供給される音速プロファイルデータを受ける。また、図１５に示すように、ふ仰角判定部８０５の代わりに、ふ仰角判定部８２５を備えている。ふ仰角判定部８２５は、ふ仰角判定部８０５と同様の機能を有するほか、以下のような付加的な機能をも有する。
【００６５】
図１６は、ふ仰角判定部８２５の詳細を示す。図示のふ仰角判定部８２５は、到来時間差最小組合せ算出部８５３と、受波器間最大到来時間差算出部８５５と、観測値・最大値時間差比較部８５４と、データ抽出部８５８と、ふ仰角上下判定部８５９とを有する。
【００６６】
到来時間差最小組合せ算出８５３では、入力端子８０１、８０２を介して供給される第１及び第２の到来時刻決定部２０５、２０６の出力に基づき、２つの受波器２０１、２０２における到来波の到来時刻（タイミング）のうちで、差が最も小さい組合せを求める。その動作を、図１７（ａ）及び（ｂ）を用いて説明する。第１の到来時刻決定部２０５で検出された、上側の受波器２０１における到来時刻（タイミング）をｔＡｃ（ｃ＝１，２，…）、第２の到来時刻決定部２０６で検出された下側の受波器２０２における到来時刻（タイミング）をｔＢｄ（ｄ＝１，２，…）とし、２つの受波器における到来時刻の内、差が最小となる組合せを求める。
【００６７】
例えば、時刻ｔＡｃの各々について時刻ｔＢｄのうちで最も差の小さいものを求め、時刻ｔＢｄ（ｄ＝１、２、…）のうちで、時刻ｔＡｃ（ｃ＝１、２、…）のいずれについても差が最小と判断されなかったものは、除外される。
図示の例では、時刻ｔＢ３は、時刻ｔＡ１〜ｔＡ６のいずれについても、差が最小とはならないため、除外される。
【００６８】
一方、時刻ｔＡ４と時刻ｔＡ５のいずれも、時刻ｔＢ５との差が最小であるが、時刻ｔＡ５よりも時刻ｔＡ４の方が時刻ｔＢ５との差が小さいので、時刻ｔＡ４が残され（時刻ｔＢ５と対を成すものとして選択され）、時刻ｔＡ５は除外される。
【００６９】
以上の結果、時刻ｔＡ１と時刻ｔＢ１、時刻ｔＡ２と時刻ｔＢ２、時刻ｔＡ３と時刻ｔＢ４、時刻ｔＡ４と時刻ｔＢ５、時刻ｔＡ６と時刻ｔＢ６が互いに対を成すものとして（従って同じ到来波の検出タイミングとして）選択され、時刻ｔＡ５と時刻ｔＡ３はノイズによるものと見なされ除外される。
【００７０】
受波器間最大到来時間差算出部８５５では、入力端子２１５から入力した受波器間距離ＤＷＲと、入力端子２１４から入力した音速プロファイルの、受波器２０１、２０２の近傍における音速の値Ｃから、
ｔｍａｘ＝ＤＷＲ／Ｃ …（８）
により、同じ伝搬経路を辿った音波の、受波器２０１における到来のタイミングと、受波器２０２における到来のタイミングの時間差の理論的最大値ｔｍａｘを計算する。受波器２０１、２０２は鉛直方向に整列しているので、到来ふ仰角の絶対値が大きいほど上記の時間差は大きくなる。
【００７１】
観測値・最大値時間差比較部８５４では、組み合わせ算出部８５３で決めた到来のタイミング（時刻）の対の時間差（観測値）ｔｏｂｓが受波器間最大到来時間差ｔｍａｘを超えているか否か、即ち、
ｔｏｂｓ＞ｔｍａｘ
であるか否かを判定し、判定結果ＹＶＮを出力する。判定結果ＹＶＮは、データ抽出部８５８に供給されるほか、出力端子８６２を介して図１５のデータ抽出部８０６、８０７に供給される。判定結果ＹＶＮは、第１及び第２の受波器２０１、２０２における到来時刻の対が、同じ到来波の検出タイミングである蓋然性が高い（信頼性が高い）か否かを示し、ｔｏｂｓ≦ｔｍａｘであれば、上記蓋然性が高いことを意味し、従って、当該到来時刻のデータがふ仰角の判定、到来時間差の比較などに利用するのに適していることを意味し、逆にｔｏｂｓ＞ｔｍａｘであれば、上記蓋然性が低く、当該到来時刻のデータがふ仰角の判定、到来時間差の比較などに利用するのに適さないことを意味する。
【００７２】
データ抽出部８５８では、観測値・最大値時間差比較部８５４による判定結果ＹＶＮに基づいて、入力端子８０１、８０２を介して供給される第１及び第２の到来時刻決定部２０５、２０６の出力を、ふ仰角の判定に使用するかしないかを決定し、ふ仰角の判定に使用すると決定した時刻を表すデータ（到来時刻決定部２０５、２０６の出力）を抽出して、ふ仰角上下判定部８５９に供給する。即ち、各時刻対（互いに対をなすと判断された検出時刻）に関し、
ｔｏｂｓ＞ｔｍａｘ
の場合には、当該時刻対を、ふ仰角の判定に使用しないこととし、ふ仰角上下判定部８５９に供給しない。
【００７３】
一方、各時刻対に関し、
ｔｏｂｓ≦ｔｍａｘ
の場合には、当該時刻対を表すデータを抽出して、ふ仰角上下判定部８５９に供給する。
【００７４】
ふ仰角上下判定部８５９では、データ抽出部８５８から供給された時刻対を表すデータについて、ふ仰角が上向きか下向きかの判定を行い、判定結果ＹＵＬを、接続端子８６１を介して出力する。
【００７５】
データ抽出部８０６は、ふ仰角判定部８２５における判定結果（ふ仰角が上向きか下向きかの判定結果ＹＵＬ、及び各検出タイミングｔＡｃ，ｔＢｄの信頼性についての判定結果ＹＶＮ）に応じて、第１のマルチパス波到来時間差算出部２０７の出力（入力端子８０３を介して供給される）のうちの、データ抽出部８５８で抽出された到来時刻に対応し、かつふ仰角上下判定部８５９でふ仰角が上向きと判定された到来時間差データ値のみを抽出して、ふ仰上向きテーブル値一観測値差算出器８０８に供給する。
【００７６】
データ抽出部８０７は、ふ仰角判定部８０５における判定結果（ふ仰角が上向きか下向きかの判定結果ＹＵＬ、及び各検出タイミングｔＡｃ，ｔＢｄの信頼性についての判定結果ＹＶＮ）に応じて、第２のマルチパス波到来時間差算出部２０８の出力（入力端子８０４を介して供給される）のうちの、データ抽出部８５８で抽出された到来時刻に対応し、かつふ仰角上下判定部８５９でふ仰角が下向きと判定された到来時間差データ値のみを抽出して、ふ仰下向きテーブル値一観測値差算出器８０９に供給する。
【００７７】
図１５のふ仰上向きテーブル値一観測値差算出器８０８及びふ仰下向きテーブル値一観測値差算出器８０９では、データ抽出部８０６及び８０７で抽出された到来時間差の観測値を用いて、差分の算出及びそれに基づくコスト関数の算出を行なう。コスト関数算出器２２０の他の部分の動作は実施の形態１について説明したのと同様である。
【００７８】
実施の形態１では雑音の影響の無いほぼ理想的な状況にのみ到来音波のふ仰角の上下を判別したコスト関数を利用することができたが、本実施の形態２で説明したように到来音波のふ仰角の判定を行う（信頼性の高いデータのみを用いてふ仰角の判定、到来時間差の比較（不一致度の検出）を行う）ことにより、ＳＮ比がさほど高く無く到来時間差に誤検出や未検出がある場合にも、到来音波のふ仰上下を判別したコスト関数を適用することが可能となる。
【００７９】
このため、ＳＮ比がさほど高く無く到来時間差に誤検出や未検出がある場合にもコスト関数を円錐状に保つことができ、音源の位置推定精度を保つことができる。
【００８０】
尚、実施の形態２では、ふ仰角の計測に２つの受波器の時間差を用いているが、３つ以上の受波器を用いて、ふ仰角の上下の判定を行うこととしても良い。
【００８１】
また、音速プロファイルを測定できなかった場合については、式（８）の音速プロファイルの受波器２０１，２０２の近傍における音速の値Ｃの代わりに代表的な音速１５００ｍ／ｓを用いて、
ｔｍａｘ＝ＤＷＲ／１５００ …（９）
を利用することも可能である。
【符号の説明】
【００８２】
２０１、２０２ 受波器、 ２０５、２０６ 到来時刻決定部、 ２０７、２０８ マルチパス波到来時間差算出部、 ２０９ 到来時間差テーブル作成器、 ２１０ コスト関数算出器、 ２１１ 音源位置推定器。
【特許請求の範囲】
【請求項１】
受波位置における到来波に基づいて音源の位置を推定する音源推定装置において、
複数の仮想音源位置の各々について、該仮想音源位置に仮想音源が存在すると仮定したときの、前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が上向きのものの各々との到来時間差の理論値、並びに前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が下向きのものの各々との到来時間差の理論値を格納する到来時間差理論値格納手段と、
前記受波位置において観測した到来波のデータに基づいて、目標の音源からの最初の到来波と２番目以降の到来波との到来時間差の観測値を算出する到来時間差観測値算出手段と、
前記２番目以降の到来波の各々の到来ふ仰角が上向きか下向きかを判定し、前記２番目以降の到来波のうち、到来ふ仰角が上向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が上向きの到来波についての前記到来時間差の理論値との不一致度を第１の不一致度として求め、前記２番目以降の到来波のうち、到来ふ仰角が下向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が下向きの到来波についての前記到来時間差の理論値との不一致度を第２の不一致度として求め、前記第１の不一致度と前記第２の不一致度との和をコストとして求めるコスト関数算出手段と、
前記コスト関数算出手段が算出した前記コストが最小となる前記仮想音源位置を前記目標の音源の位置と推定する音源位置推定手段と
を備えた音源推定装置。
【請求項２】
前記受波位置に複数個の受波器が設けられ、
前記コスト関数算出手段が、前記複数個の受波器の出力に基づいて前記到来ふ仰角が上向きか下向きかの判定を行う
ことを特徴とする請求項１に記載の音源推定装置。
【請求項３】
前記受波位置の、深度方向の異なる位置に第１及び第２の受波器が設けられ、
前記コスト関数算出手段が、同じ到来波の、前記第１及び第２の受波器での受波時刻の差に基づいて、当該到来波の到来ふ仰角が上向きか下向きかの判定を行う
ことを特徴とする請求項１に記載の音源推定装置。
【請求項４】
前記第１の受波器で検出された到来波の各々について、その検出のタイミングと、
前記第２の受波器において検出された複数の到来波のうち、その検出のタイミングと前記第１の受波器における検出のタイミングとの差が最小であり、かつ該差が所定値以下であるものを、前記第１の受波器で検出された到来波と同じ到来波であると判定し、
該判定により同じと判定された到来波の検出のタイミングに基づいて前記到来ふ仰角が上向きか下向きかの判定を行う
ことを特徴とする請求項３に記載の音源位置推定装置。
【請求項５】
前記所定値が、２つの受波器間における前記検出のタイミングの差が、当該受波器間の距離及び前記受波位置における音速に基づいて計算される、前記検出のタイミングの差の理論的最大値以下であるか否かを判断するための値であることを特徴とする請求項４に記載の音源位置推定装置。
【請求項６】
前記同じと判定された到来波についての到来時間差の観測値のみを前記到来時間差の理論値との不一致度の算出に用いることを特徴とする請求項４に記載の音源位置推定装置。
【請求項７】
前記複数の受波器が鉛直受波器アレイを構成するものであることを特徴とする請求項２に記載の音源位置推定装置。
【請求項８】
前記到来時間差理論値格納手段に格納されている到来時間差の理論値が、受波位置の深度、仮想音源位置の距離、深度、及び音速プロファイルに基づいて算出されたものであることを特徴とする請求項１に記載の音源位置推定装置。
【請求項９】
前記コスト関数算出手段は、
前記到来ふ仰角が上向きの到来波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第１の不一致度として計算し、
前記到来ふ仰角が下向きの到来波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第２の不一致度として計算する
ことを特徴とする請求項１に記載の音源位置推定装置。
【請求項１】
受波位置における到来波に基づいて音源の位置を推定する音源推定装置において、
複数の仮想音源位置の各々について、該仮想音源位置に仮想音源が存在すると仮定したときの、前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が上向きのものの各々との到来時間差の理論値、並びに前記受波位置における最初の到来波と２番目以降の到来波のうちの到来ふ仰角が下向きのものの各々との到来時間差の理論値を格納する到来時間差理論値格納手段と、
前記受波位置において観測した到来波のデータに基づいて、目標の音源からの最初の到来波と２番目以降の到来波との到来時間差の観測値を算出する到来時間差観測値算出手段と、
前記２番目以降の到来波の各々の到来ふ仰角が上向きか下向きかを判定し、前記２番目以降の到来波のうち、到来ふ仰角が上向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が上向きの到来波についての前記到来時間差の理論値との不一致度を第１の不一致度として求め、前記２番目以降の到来波のうち、到来ふ仰角が下向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が下向きの到来波についての前記到来時間差の理論値との不一致度を第２の不一致度として求め、前記第１の不一致度と前記第２の不一致度との和をコストとして求めるコスト関数算出手段と、
前記コスト関数算出手段が算出した前記コストが最小となる前記仮想音源位置を前記目標の音源の位置と推定する音源位置推定手段と
を備えた音源推定装置。
【請求項２】
前記受波位置に複数個の受波器が設けられ、
前記コスト関数算出手段が、前記複数個の受波器の出力に基づいて前記到来ふ仰角が上向きか下向きかの判定を行う
ことを特徴とする請求項１に記載の音源推定装置。
【請求項３】
前記受波位置の、深度方向の異なる位置に第１及び第２の受波器が設けられ、
前記コスト関数算出手段が、同じ到来波の、前記第１及び第２の受波器での受波時刻の差に基づいて、当該到来波の到来ふ仰角が上向きか下向きかの判定を行う
ことを特徴とする請求項１に記載の音源推定装置。
【請求項４】
前記第１の受波器で検出された到来波の各々について、その検出のタイミングと、
前記第２の受波器において検出された複数の到来波のうち、その検出のタイミングと前記第１の受波器における検出のタイミングとの差が最小であり、かつ該差が所定値以下であるものを、前記第１の受波器で検出された到来波と同じ到来波であると判定し、
該判定により同じと判定された到来波の検出のタイミングに基づいて前記到来ふ仰角が上向きか下向きかの判定を行う
ことを特徴とする請求項３に記載の音源位置推定装置。
【請求項５】
前記所定値が、２つの受波器間における前記検出のタイミングの差が、当該受波器間の距離及び前記受波位置における音速に基づいて計算される、前記検出のタイミングの差の理論的最大値以下であるか否かを判断するための値であることを特徴とする請求項４に記載の音源位置推定装置。
【請求項６】
前記同じと判定された到来波についての到来時間差の観測値のみを前記到来時間差の理論値との不一致度の算出に用いることを特徴とする請求項４に記載の音源位置推定装置。
【請求項７】
前記複数の受波器が鉛直受波器アレイを構成するものであることを特徴とする請求項２に記載の音源位置推定装置。
【請求項８】
前記到来時間差理論値格納手段に格納されている到来時間差の理論値が、受波位置の深度、仮想音源位置の距離、深度、及び音速プロファイルに基づいて算出されたものであることを特徴とする請求項１に記載の音源位置推定装置。
【請求項９】
前記コスト関数算出手段は、
前記到来ふ仰角が上向きの到来波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第１の不一致度として計算し、
前記到来ふ仰角が下向きの到来波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第２の不一致度として計算する
ことを特徴とする請求項１に記載の音源位置推定装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【公開番号】特開２０１１−１４９８６４（Ｐ２０１１−１４９８６４Ａ）
【公開日】平成２３年８月４日（２０１１．８．４）
【国際特許分類】
【出願番号】特願２０１０−１２１８６（Ｐ２０１０−１２１８６）
【出願日】平成２２年１月２２日（２０１０．１．２２）
【出願人】（００００００２９５）沖電気工業株式会社 (6,645)
【Ｆターム（参考）】
【公開日】平成２３年８月４日（２０１１．８．４）
【国際特許分類】
【出願日】平成２２年１月２２日（２０１０．１．２２）
【出願人】（００００００２９５）沖電気工業株式会社 (6,645)
【Ｆターム（参考）】
[ Back to top ]