(57)【特許請求の範囲】
【請求項１】
複数の受波器から成る受波器群で受信した到来波に基づいて音源の位置を推定する音源推定装置において、
一連の到来波列における到来波間の到来時間差の検出値又は該到来時間差の複数回の測定結果の平均値と、各到来波の到来ふ仰角の検出値又は該到来ふ仰角の複数回の測定結果の平均値を求めるとともに、前記各到来波の到来時間差の標準偏差及び前記各到来波の到来ふ仰角の標準偏差を出力する到来波解析部と、
前記一連の到来波列における到来波の到来順序とは逆の順序で、かつ相前後するパルス波間の時間差を前記到来波の到来時間差又はその平均値と同じにしたパルス波列を構成する複数のパルス波を、それぞれ対応する到来波の到来ふ仰角又はその平均値と同じ送出ふ仰角で、前記受波器群から順次送出した場合の各パルス波の伝搬経路を音線計算により推定し、各時刻に上記推定した伝搬経路上及びその周囲の各位置に、前記各パルス波が存在する尤度を推定する逆伝搬推定部と、
前記パルス波列の先頭のパルス波の送出後の各時刻において、前記パルス波列のそれぞれのパルス波が前記各位置に存在する尤度の、すべてのパルス波についての合計を求める統合尤度計算部と、
前記逆伝搬推定部で、前記複数のパルス波について、それぞれ推定された各時刻の尤度の分布間の距離を求めて、該距離が最小となる時刻を収束時刻と推定する収束時刻判定部と、
前記収束時刻判定部で推定された収束時刻において、前記統合尤度計算部で算出された前記尤度の合計が最大となる位置を音源位置と推定する収束位置判定部と
を有する音源位置推定装置。
【請求項２】
前記逆伝搬推定部は、
前記推定された伝搬経路上の位置を尤度の中心とし、前記標準偏差に基づいて、前記尤度の中心及びその周囲の各位置に、各パルス波が存在する尤度を求め、各位置と尤度との関係を尤度分布として出力する
ことを特徴とする請求項１に記載の音源位置推定装置。
【請求項３】
前記収束時刻判定部が、
前記複数のパルス波から２つを取り出したすべての組合せの各々を構成する２つのパルス波についての尤度の分布間の距離を求め、該距離の、すべての組み合わせについての合計を求めて、前記複数のパルス波についての尤度の分布間の距離として出力する
ことを特徴とする請求項１又は２に記載の音源位置推定装置。
【請求項４】
前記逆伝搬推定部は、
前記一連の到来波列における第１乃至第ｎの到来波の到来順序とは逆の順序で、かつ相前後するパルス波間の時間差を前記到来波の到来時間差又はその平均値と同じにしたパルス波列を構成する第１乃至第ｎのパルス波の仮想的送出のタイミングを示す情報及び該パルス波の到来ふ仰角を示す情報を生成する順序反転部と、
前記順序反転部で生成された第１乃至第ｎのパルス波についての前記情報をそれぞれ入力とし、それぞれ第１乃至第ｎのパルス波について、それぞれ前記仮想的送出のタイミングで送出した場合の、それぞれのパルス波の伝搬経路の計算を行う第１乃至第ｎの伝搬計算部と、
前記第１乃至第ｎの伝搬計算部による前記第１乃至第ｎのパルス波についての伝搬経路の計算結果に基づいて各位置における尤度を計算する第１乃至第ｎの尤度計算部とを有し、
第ｉの伝搬計算部（ｉは１乃至ｎのいずれか）は、第ｉのパルス波を、対応する第ｉの到来波の到来時間差に対応した前記仮想的送出のタイミングで、前記第ｉの到来波の到来ふ仰角又はその平均値と同じ送出ふ仰角で、前記受波器群から、送出した場合の第ｉのパルス波の伝搬経路を音線計算により推定し、
第ｉの尤度計算部は、前記第ｉの伝搬計算部で推定された第ｉのパルス波の伝搬経路上の位置を尤度の中心とし、前記第ｉの到来波の到来時間差の標準偏差及び前記第ｉの到来波の到来ふ仰角の標準偏差に基づいて、前記第ｉのパルス波の前記尤度の中心及びその周囲の各位置に、前記第ｉのパルス波が存在する尤度を求め、各位置と尤度との関係を第ｉのパルス波についての尤度分布として出力する
ことを特徴とする請求項１に記載の音源位置推定装置。
【請求項５】
前記収束時刻判定部が、
前記第１乃至第ｎのパルス波から２つを取り出したすべての組合せについて各時刻における尤度の分布間の距離を求める複数の分布間距離計算部と、
前記複数の分布間距離計算部で求めた分布間距離の合計を求める距離加算部とを有し、
各分布間距離計算部は、対応する組合せを構成する第ｉのパルス波についての尤度の分布と、第ｊ（ｊは１からｎのいずれか、但しｊ≠ｉ）のパルス波についての尤度の分布の距離を求める
ことを特徴とする請求項４に記載の音源位置推定装置。
【請求項６】
前記尤度の分布が正規分布であると仮定して前記各位置の尤度を求める
ことを特徴とする請求項２又は４に記載の音源位置推定装置。
【請求項７】
前記到来波解析部が、
一連の到来波列における到来波間の到来時間差の検出値の複数回の測定結果に基づいて各到来波の到来時間差の標準偏差を求め、
一連の到来波列における各到来波のふ仰角の複数回の測定結果に基づいて各到来波の到来ふ仰角の標準偏差を求める
ことを特徴とする請求項１乃至６のいずれかに記載の音源位置推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は音源位置推定装置に関し、特に水中の音源から放射されたパルス状の音波の直接波とマルチパス波を受信し、これらの到来時間差から音源の位置を推定する装置に関する。
【背景技術】
【０００２】
水中の音波伝搬では、音源から直接受波器に到来する音波（直接波）の他に、海面や海底で反射して到来するマルチパス波があり、そのようなマルチパス波は、直接波と異なる時刻に到来する。水中において音源からのパルス状の音波（パルス波）を受信して音源の位置を推定する技術として、下記の特許文献１に示されるものがある。
【０００３】
特許文献１に示される従来の音源位置推定装置においては、音響センサー（受波器）で受信した音波（到来波）を表す信号に基づき直接波とマルチパス波の到来時間差を検出し、検出した到来時間差観測値として出力する。一方、水中の音速プロファイルに基づき、複数の仮想音源の各々からの音波の音線計算を行うことにより音波の到来時間を推定し、推定結果を到来時間差テーブル値（理論値）として、記憶しておく。そして、到来時間差の観測値とテーブル値より、それらの不一致度を反映したコスト関数の値を算出し、該コスト関数が最小となる仮想音源位置を、推定対象（目標）の音源位置と推定する。
【０００４】
上記の従来の装置を用いた場合、音源と受波器を結ぶ線上のコスト関数の値が、音源の位置の前後のかなり広い範囲にわたり小さな値となるため、音源位置を正確に推定することが困難であるという問題があった。即ち、到来時間差は観測誤差を含むため、コスト関数の値が小さい場所で大きく推定位置がずれる可能性が大きく、このため観測誤差が比較的大きいという問題があった。
【０００５】
そこで、到来時間差観測値に誤差がある場合にも、音源位置の推定誤差を小さくするため、到来ふ仰角が上向きのマルチパス波と到来ふ仰角が下向きのマルチパス波について別個に到来時間差のテーブル値を記憶しておき、受信したマルチパス波の各々について到来ふ仰角が上向きか下向きかの判定を行い、観測された到来時間差とテーブル値との不一致度を上向きのふ仰角を有するマルチパス波と下向きのふ仰角を有するマルチパス波の各々について別個に算出して、その和をコスト関数として、音源位置の推定を行うことが提案されている（特許文献２）。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開２００６−１９４６２７号公報
【特許文献２】特開２０１１−１４９８６４号公報
【非特許文献】
【０００７】
【非特許文献１】Ｆｉｎｎ  Ｂ．Ｊｅｎｓｅｎ，Ｗｉｌｌｉａｍ  Ａ．  Ｋｕｐｅｒｍａｎ，  Ｍｉｃｈａｅｌ  Ｂ．  Ｐｏｒｔｅｒａｎｄ  Ｈｅｎｒｉｋ  Ｓｃｈｍｉｄｔ，“Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｏｃｅａｎ  Ａｃｏｕｓｔｉｃ”，Ｃｈｐ．３，  Ａｌｐ  Ｐｒｅｓｓ
【０００８】
なお、非特許文献１については後に言及する。
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、特許文献２の方法でも、音源位置の推定を正確に行うことができない場合があった。その要因としては、以下の２つのことが挙げられる。
一つは、到来時間差の観測値とテーブルとの対応関係の誤りである。これは、特許文献２の方法を用いることで、特許文献１の方法よりも減少させることができるが、十分とはいえない。他の一つは、到来波についての、誤検出や検出漏れによる推定結果の誤りである。
【００１０】
本発明の目的は、到来波の誤検出や検出漏れがあっても、音源位置の正確な推定を可能にすることにある。
【課題を解決するための手段】
【００１１】
本発明の音源推定装置は、
複数の受波器から成る受波器群で受信した到来波に基づいて音源の位置を推定する音源推定装置において、
一連の到来波列における到来波間の到来時間差の検出値又は該到来時間差の複数回の測定結果の平均値と、各到来波の到来ふ仰角の検出値又は該到来ふ仰角の複数回の測定結果の平均値を求めるとともに、前記各到来波の到来時間差の標準偏差及び前記各到来波の到来ふ仰角の標準偏差を出力する到来波解析部と、
前記一連の到来波列における到来波の到来順序とは逆の順序で、かつ相前後するパルス波間の時間差を前記到来波の到来時間差又はその平均値と同じにしたパルス波列を構成する複数のパルス波を、それぞれ対応する到来波の到来ふ仰角又はその平均値と同じ送出ふ仰角で、前記受波器群から順次送出した場合の各パルス波の伝搬経路を音線計算により推定し、各時刻に上記推定した伝搬経路上及びその周囲の各位置に、前記各パルス波が存在する尤度を推定する逆伝搬推定部と、
前記パルス波列の先頭のパルス波の送出後の各時刻において、前記パルス波列のそれぞれのパルス波が前記各位置に存在する尤度の、すべてのパルス波についての合計を求める統合尤度計算部と、
前記逆伝搬推定部で、前記複数のパルス波について、それぞれ推定された各時刻の尤度の分布間の距離を求めて、該距離が最小となる時刻を収束時刻と推定する収束時刻判定部と、
前記収束時刻判定部で推定された収束時刻において、前記統合尤度計算部で算出された前記尤度の合計が最大となる位置を音源位置と推定する収束位置判定部と
を有することを特徴とする。
【発明の効果】
【００１２】
本発明によれば、誤検出や検出漏れを考慮に入れた尤度分布に基づいて音源位置の推定を行っているので、伝搬経路到来波の誤検出や検出漏れがあっても誤差を大きくすることなく音源位置を正確に推定することができる。
【図面の簡単な説明】
【００１３】
【図１】本発明の実施の形態１の音源位置推定装置を示すブロック図である。
【図２】音源位置ＷＳの音源からの伝搬経路の一例を示す図である。
【図３】実施の形態１における受波器の配置を示す図である。
【図４】受波器におけるサンプル値の絶対値の包絡線の一例を示すグラフである。
【図５】（ａ）及び（ｂ）は、一対の受波器における受信のタイミングの一例を示す図である。
【図６】図１の逆伝搬推定部の構成例を示すブロック図である。
【図７】（ａ）及び（ｂ）は、図１の逆伝搬推定部により検出された到来波列、及び図６の順序反転部により生成される順序を反転したパルス波列を示す図である。
【図８】音線計算のための局所座標を示す図である。
【図９】２つのパルス波についての到達位置の尤度分布の例を示す図である。
【図１０】図１の収束時刻判定部の構成例を示す図である。
【図１１】シミュレーションのために仮定した音源位置と受波器群の位置の例を示す図である。
【図１２】図１１の配置の場合に得られる、逆伝搬時間と分布間距離の総和の関係を示すグラフである。
【図１３】（ａ）〜（ｃ）は、パルス波の送出からの時間経過に対する、到達位置の尤度分布の変化を示す図である。
【図１４】図１４は、図１３（ｃ）の推定音源位置の部分を拡大して示す図である。
【発明を実施するための形態】
【００１４】
実施の形態１．
本発明の実施の形態１の音源位置推定装置を図１に示す。図示の音源位置推定装置は、音源から放射されたパルス状音波（パルス波）の直接波とマルチパス波を受信し、それらの到来時間差に基づいて音源の位置を推定するものであり、音源の推定位置は深度及び受波器からの水平方向距離で表される。なお、音源の位置する方角を示す情報は別途得られるものとする。
【００１５】
図示の音源位置推定装置は、第１及び第２の受波器２０１ａ、２０１ｂと、第１及び第２の到来時刻決定部２０２ａ及び２０２ｂと、到来波解析部２０３と、逆伝搬推定部２０５と、統合尤度計算部２０６と、収束時刻判定部２０７と、収束位置判定部２０８と、データベース２０９とを有する。
【００１６】
データベース２０９には、音速プロファイル２０９ａ、海底地形２０９ｂ、受波器位置２０９ｃを示す情報が保存されている。
音速プロファイル２０９ａを示す情報としては、「現場で実測した音速」、又は「温度分布、塩分濃度から経験式により算出した音速プロファイル」を用いることが望ましい。これらが得られない場合は、一般に市販されているデータベースの値を用いるか、又は音速は一定の値として処理する（例えば通常良く用いられている１５００ｍ／ｓ一定値を用いる）ことができる。
【００１７】
海底地形２０９ｂを示す情報としては、実測値もしくは一般に市販されているデータベースの値などを用いる。
受波器位置２０９ｃを示す情報としては、受波器２０１ａ、２０１ｂが設定されている深度、及び水平方向の位置を示す情報が入力される。
【００１８】
音源ＷＳからの伝搬経路は例えば図２に示すように想定される。図２には、簡単のため、５つの伝搬経路ＷＰ１〜ＷＰ５のみが示されている。伝搬経路（音線：音のエネルギーの流れ）が直線的でないのは、水温、塩分濃度などの不均一性のため、音速が一定ではないことによる。
一度も反射をせずに受波点ＷＲに達した波は直接波と呼ばれ、海面や海底で１回以上反射して受波器に達した波はマルチパス波と呼ばれる。
【００１９】
音源と受波器群との間に他の音源がない状態であると、音波伝搬は自由空間の波動方程式で表される。自由空間の波動方程式の解はＧｒｅｅｎ関数で表される。Ｇｒｅｅｎ関数の一つの大きな特徴として、相反性がある。本発明に関する特性としては、相反性が成り立つため、音源から出て受波器群に到達した直接波及びマルチパス波の順序を逆にして（時間反転して）受波器群から逆方向に伝搬させると、受波器群から出た音波は音源位置に収束する。
【００２０】
これは、波動方程式を近似した音線計算でも成り立ち、受信時の到来時間、到来角を観測し、到来波の順序を逆にして、到来角と同じ角度（同じふ仰角および方角）で逆向きに受波器群から音源方向にパルス波を送出したと仮定して、送出された各パルス波の各時点における到達位置を音線計算により推定すると、各パルス波の到達位置は時間とともに変化するが、すべてのパルス波が同時に音源位置に到達する（収束）するはずである。本発明では、この性質を利用して、音源位置を推定する。
【００２１】
第１の受波器２０１ａ及び第２の受波器２０１ｂは、到来波を受信するためのものであり、到来波のふ仰角を求めるため、同じ鉛直線上の異なる位置に設けられている。受波器２０１ａ及び２０１ｂにより受波器群２０１が構成されている。
【００２２】
第１及び第２の受波器２０１ａ及び２０１ｂは、例えば、図３に示すように、間隔ＤＷＲを開けて配置されている。間隔ＤＷＲは受信位置において想定される音速の最大値と受波器２０１ａ、２０１ｂで測定可能な最小時間差の積以上の値に設定される。
第１及び第２の受波器２０１ａ及び２０１ｂはそれぞれ上側及び下側に配置されるので、「上側の受波器」、「下側の受波器」と呼ばれることもある。
【００２３】
上側及び下側の受波器２０１ａ及び２０１ｂの出力は、それぞれ第１及び第２の到来時刻決定部２０２ａ及び２０２ｂに入力される。
【００２４】
第１及び第２の受波器２０１ａ及び２０１ｂは、それらで受信される音波をサンプリングすることで得られるサンプル値を表すデータを出力するが、音波をサンプリングすることで得られるサンプル値の絶対値を時間軸上に並べたものの包絡線は、例えば、図４に示す如くとなる。なお、サンプル値の絶対値の代りにサンプル値の自乗値の包絡線を用いる場合もある。
【００２５】
第１及び第２の到来時刻決定部２０２ａ及び２０２ｂでは、それぞれ上側及び下側の受波器２０１ａ及び２０１ｂからのデータの列に基づき到来波の受信時刻（到来時刻）を判断する。到来時刻の判断には、従来と同様に、レプリカ相関による方法、立ち上がり検出による方法などを用いることができる。これらの方法は、例えば上記特許文献１において到来時間差測定部の動作の一部として説明されている。
【００２６】
図５（ａ）及び（ｂ）に、第１及び第２の到来時刻決定部２０２ａ及び２０２ｂで決定乃至検出された上側の受波器２０１ａ及び下側の受波器２０１ｂにおけるパルス波の到来時刻（検出タイミング）を示す。図５（ａ）は、上側の受波器２０１ａにおける検出タイミングを、早い順にＴＡ１、ＴＡ２、…で示し、図５（ｂ）は、下側の受波器２０１ｂにおける検出タイミングを、早い順にＴＢ１、ＴＢ２、…で示す。
【００２７】
なお、到来波の到来時間差と到来ふ仰角を観測する方法は、所望の情報が観測できる方法であれば何れの方法も適用可能で、２つの受波器間での受信タイミングに関する情報を用いる他に、受波器アレイを構成して整相処理を行い、その結果から所望の情報を得ることも可能である。
【００２８】
第１及び第２の到来時刻決定部２０２ａ及び２０２ｂの出力は、到来波解析部２０３に入力される。
【００２９】
到来波解析部２０３は、到来時刻決定部２０２ａ及び２０２ｂの出力を受け、一連のパルス波を構成する複数のパルス波（直接波及びマルチパス波）の各々についての到来時間差及び到来ふ仰角を算出する。到来時間差は、一連のパルス波のうちの先頭のパルス波の受信（到来）から、各パルス波の各々の受信（到来）までの時間差を意味する。先頭のパルス波については、到来時間差はゼロであるとして扱えば、以降の処理を一律に行うことができる。到来時間差及び到来ふ仰角の算出に当たり、以下のような処理をする。
【００３０】
まず、上側の受波器２０１ａと下側の受波器２０１ｂとでは同じパルス波の受信時刻が異なるが、この受信時刻の差異、並びに受波器２０１ａ、２０１ｂ間の間隔、及び受波器群の配置されている位置における音速に基づいてふ仰角を算出する。また、到来時間差の算出に当たっては、２つの受波器２０１ａ及び２０１ｂにおける受信時刻の中間値を受波器群２０１における受信時刻とする。
【００３１】
一連のパルス波のみに基づいて（一連のパルス波の先頭のパルス波と後続のパルス波の受信に基づいて）、到来時間差や到来ふ仰角を求めることもできるが、そのようにした場合には、到来波の検出漏れ、到来波の誤検出による測定誤差があり得る。そこで、本実施の形態では、一連のパルス波の受信を複数回繰り返した結果に基づいて到来時間差の平均値を取り、平均値を「到来時間差」として用いる。到来ふ仰角についても同様に、複数回の受信結果の平均値を用いる。
これとともに、複数回の測定結果に基づき、各到来波について観測された到来時間差のばらつきの程度を表す指標としての到来時間差の標準偏差ｉσｔと、各到来波について観測された到来ふ仰角のばらつきの程度を表す指標として到来ふ仰角の標準偏差ｉσθを求める。
【００３２】
なお、複数回の測定結果を行う間に音源又は受波器群の位置が変化する可能性があるが、その変化は十分に小さいものとする。逆に言えば、複数回の測定を、音源や受波器群の位置の変化に対して十分に短い時間で行うものとする。
【００３３】
到来波解析部２０３で求められた一連の到来波についての到来時間差（平均値）及び到来ふ仰角（平均値）を示す情報は、逆伝搬推定部２０５に供給される。
【００３４】
逆伝搬推定部２０５は、到来波解析部２０３で求められた一連の到来波列における到来波の到来順序とは逆の順序で、かつ相前後するパルス波間の時間差を、到来波の到来時間差（平均値）と同じにしたパルス波列を構成する複数のパルス波を、それぞれ対応する到来波の到来ふ仰角（平均値）と同じ送出ふ仰角で、受波器群２０１から順次送出した場合の各パルス波の伝搬経路を音線計算により推定し、各時刻に、推定した伝搬経路上及びその周囲の各位置に、各パルス波が存在する尤度を推定する。具体的には、推定された伝搬経路上の位置を尤度の中心とし、標準偏差ｉσｔ、ｉσθに基づいて、尤度の中心及びその周囲の各位置に、各パルス波が存在する尤度を求め、各位置と尤度との関係を尤度分布として出力する。この場合、尤度の分布が正規分布であると仮定する。
【００３５】
上記の伝搬経路の推定に当たっては、データベース２０９に記憶されている音速プロファイル２０９ａ、海底地形２０９ｂ、及び受波器位置２０９ｃを示す情報をも利用する。
【００３６】
逆伝搬推定部２０５は、例えば、図６に示すように、順序反転部５１と、第１乃至第Ｎの伝搬計算部５２−１〜５２−Ｎと、第１乃至第Ｎの尤度計算部５３−１〜５３−Ｎと、第１乃至第Ｎの分散共分散計算部５４−１〜５４−Ｎとを有する。
【００３７】
順序反転部５１は、到来波解析部２０３から逆伝搬推定部２０５に入力された到来時間差及び到来ふ仰角についての情報に基づいて、一連のパルス波（到来波）の順序を逆にしたパルス波列を表す情報を生成する。
即ち、到来波解析部２０３から供給された情報が、第１、第２、…第Ｎのパルス波を順に受信した場合、順序を逆転し、第Ｎ、第Ｎ−１、…第１のパルス波の順に並べ替えたパルス波列を表す情報を生成する。この場合第ｉ（ｉは１〜Ｎ−１のいずれか）のパルス波と、第ｉ＋１のパルス波の間隔は、受信した到来波列における間隔と同じにする。後述のように、これらのパルス波についての情報は、該パルス波を仮想的に送出した場合の伝搬経路の推定に利用されるものである。
【００３８】
具体的には、到来波列の第ｉの到来波（ｉ番目に到来したパルス波）についての到来時間差（先頭の到来波の受信時刻と第ｉの到来波の受信時刻の差）をＦｉとし、第Ｎ（最後）の到来波についての到来時間差をＦＮとすると、当該第ｉの到来波に対応する第ｉのパルス波の時間差（先頭のパルス波との時間差）Ｇｉは、
Ｇｉ＝ＦＮ−Ｆｉ      …（１）
とする。この時間差Ｇｉは、後述伝搬経路計算のための仮想的送出のタイミングを決めるために利用されるものであり、送出時間差と呼ぶ。
このような処理をする結果、到来波列における第ｉの到来波に対応する、パルス波列における第ｉのパルス波は、パルス波列において、後ろから数えたときにｉ番目に位置するものとなる。
【００３９】
例えば、到来波列を構成する到来波の数が５であり、図７（ａ）に示されるように、到来波ＦＷ１〜ＦＷ５が図示のタイミングで受信された場合、順序を逆にし、隣り合うパルス波間の時間間隔は元のままとして、図７（ｂ）に示すようなパルス波ＧＷ５〜ＧＷ１から成るパルス波列を示す情報を生成する。
式（１）より、
Ｇ５＝Ｆ５−Ｆ５＝０
Ｇ４＝Ｆ５−Ｆ４
Ｇ３＝Ｆ５−Ｆ３
Ｇ２＝Ｆ５−Ｆ２
Ｇ１＝Ｆ５−Ｆ１＝Ｆ５
の関係がある。
【００４０】
順序反転部５１では上記のようなパルス波列を構成するパルス波の送出時間差についての情報を生成するとともに、各パルス波のふ仰角を表す情報を合わせて生成する。この場合、各到来波列の各到来波についての到来ふ仰角をそのままパルス波列の対応するパルス波についての送出ふ仰角として生成する。
【００４１】
順序反転部５１の出力（パルス波列の各パルス波について送出時間差及び送出ふ仰角についての情報）は、それぞれ第１乃至第Ｎの伝搬計算部５２−１乃至５２−Ｎに供給される。即ち、第ｉのパルス波についての情報は、第ｉの伝搬計算部５２−ｉに供給される。
【００４２】
第１乃至第Ｎの伝搬計算部５２−１〜５２−Ｎは、それぞれ順序反転部５１で生成された第１乃至第Ｎのパルス波についての情報と、データベース２０９に記憶されている音速プロファイル２０９ａ、海底地形２０９ｂ、及び受波器位置２０９ｃを示す情報を入力として、それぞれ第１乃至第Ｎのパルス波について、それぞれ送出時間差で決められる仮想的送出のタイミングで送出した場合の、それぞれのパルス波の伝搬経路の計算を行う。
【００４３】
第１乃至第Ｎの伝搬計算部５２−１〜５２−Ｎの各々、即ち、第ｉの伝搬計算部５２−ｉ（ｉは１乃至Ｎのいずれか）は、順序反転部５１から供給された第ｉのパルス波についての到来時間差及び到来ふ仰角についての情報と、該第ｉのパルス波（パルス状音波）を、当該パルス波についての送出時間差Ｇｉに対応するタイミングで、当該パルス波についての送出ふ仰角で音源の方角に（到来波と同じ方角に）仮想的に送出（発射）した場合の、該パルス波の伝搬経路（音線）を推定し、各時刻（最初に送出したパルス波の送出後の経過時間）における当該パルス波の到達位置を推定する。なお、上記の「最初に送出したパルス波の送出後の経過時間」を単に逆伝搬時間ということもある。
「送出時間差Ｇｉに対応するタイミング」は、第１乃至第Ｎの伝搬計算部５２−１〜５２−Ｎ間で互いに関係付けられる。即ち、第ｉの伝搬計算部５２−ｉにおいては第ｉのパルス波について設定される送出のタイミングは、第Ｎの伝搬計算部５２−Ｎで第Ｎのパルス波について設定される送出のタイミングよりも時間差Ｇｉだけ後である。
【００４４】
伝搬経路の推定（音線計算）の方法としては、例えば、上記の非特許文献１に記載された方法を用いることができる。具体的には、深度をｚ、受波器群２０１から音源に向かう方向をｒで表し、次の式（２ａ）〜（２ｄ）を用いて計算を行う。
【数１】
0005919869.tif 000002
【００４５】
式（２ａ）乃至（２ｄ）において、Ｃは音速である。
ζ、ξは、送出ふ仰角θに基づいて
ζ＝ｓｉｎθ／Ｃ（０）
ξ＝ｃｏｓθ／Ｃ（０）
により求められる。ここで、Ｃ（０）は、受波器群の位置における音速である。
ｓは媒介変数で伝搬経路に沿う伝搬経路の長さを表す。
【００４６】
伝搬経路上の各位置に各パルス波が到達する時間は、次の式（３）により計算する。
【数２】
0005919869.tif 000003
【００４７】
式（３）の積分項は、伝搬経路に沿う伝搬時間を表す。すなわち、この積分項は、媒介変数ｓについての積分であるので、伝搬経路に沿った形で音速の逆数を積分する形で求められる。また、式（２ａ）〜（２ｄ）からは音波伝搬の開始点（パルス波が送出された受波器群）からの経路長（伝搬経路に沿う長さ）も求めることができる。
【００４８】
第ｉの伝搬計算部５２−ｉは、上記の式のｒ、ｚの代わりに、対応する第ｉのパルス波についての伝搬経路上の位置を表すｒｉｈ、ｚｉｈを用いて、適切な時間毎の、即ち、先頭パルス波（第Ｎのパルス波）送出後の所定の時間が経過する毎の、各パルス波の到達位置（ｒｉｈ（ｔ），ｚｉｈ（ｔ））、及び受波器群２０１からの経路長ｓｉ（ｔ）を数値計算で求めて、対応する第ｉの尤度計算部５３−ｉ及び第ｉの分散共分散計算部５４−ｉに出力する。
【００４９】
一方、到来波解析部２０３からは、受波器群２０１で観測した第ｉのパルス波についての到来時間差の標準偏差ｉσｔ及び到来ふ仰角の標準偏差ｉσθを示す情報が尤度計算部５３−ｉ及び第ｉの分散共分散計算部５４−ｉに供給される。
【００５０】
上記のように、第ｉの伝搬計算部５２−ｉでは、第ｉのパルス波について、対応する第ｉの到来波の到来時間差（平均値）及び到来ふ仰角（平均値）に基づいて各時刻における到達位置（尤度の中心）が算出されるが、到来波解析部２０３での解析の結果には、パルス波の誤検出、検出漏れによるものがあり、受波器群２０１で受信される到来波の到来時間差及び到来ふ仰角にはバラツキがあり、これらを考慮に入れると、各時刻における到達位置は、伝搬計算部５２−ｉで推定された伝搬経路上の位置のみならず、その周囲の位置に存在する可能性もある。そこで、各時刻に各パルス波が存在する（到達している）尤度は、上記の推定された伝搬経路上の位置を中心としてその周囲にも分布しているとみなして、これらの位置にパルス波が存在する尤度の分布を求める。
【００５１】
具体的には、到来波解析部２０３で求めた到来時間差の標準偏差ｉσｔ及び到来ふ仰角の標準偏差ｉσθを用いて、到達位置の尤度を算出する。そのようにして得られる尤度は、対応する第ｉの到来波の到来時間差（平均値）及び到来ふ仰角（平均値）に基づいて音線計算により求めた各時刻における到達位置（尤度の中心）を中心とする分布である。以下では、該分布が上記のように、正規分布であるとして説明する。
【００５２】
時刻ｔ（先頭パルス波の送出時刻からの経過時間）に第ｉのパルス波の到達位置が、位置（ｒｉ（ｔ），ｚｉ（ｔ））である尤度は次の式（４）で定義される。
Ｌｉ（ｔ）＝ｐｉ（ｒｉ（ｔ），ｚｉ（ｔ）｜Ｐａｔｈｉ）ｐ（Ｐａｔｈｉ）    （４）
【００５３】
ここで、ｐ（Ｐａｔｈｉ）は、第ｉのパルス波に対応する第ｉの到来波が生じる（第ｉの到来波が受波器群に到達する）事前確率である。
通常は、受波器群側から音源の指向性などの情報は得られないと考えられるので、ラプラスの不十分理由の原理により各到来波についての事前確率は等しいとすべきである。
【００５４】
また、条件付き確率は、具体的に確率密度関数が明確となる情報が得られれば、それに従って構成すべきである。
上記の確率密度関数が明確になる情報がない場合は、図８に示すように音線経路方向をｓ、音線経路に垂直な方向をηとして、到来時間差（平均値）及び到来ふ仰角（平均値）に対応する送出時間差及び送出ふ仰角を用いて音線計算により求めた位置（ｒｉｈ（ｔ），ｚｉｈ（ｔ））を原点とする局所座標（ｕｉ（ｔ），ｖｉ（ｔ））で表される位置（位置（ｒｉ（ｔ），ｚｉ（ｔ）で表される位置）が第ｉの送出波の到達位置である尤度は下記の式（５）で表される。尤度計算部５３−ｉでは下記の式による演算を行うことで尤度Ｌｉ（ｔ）を求める。
【００５５】
【数３】
0005919869.tif 000004
【００５６】
式（５）において、
ｓｉ（ｔ）は、上記の位置（ｒｉｈ（ｔ），ｚｉｈ（ｔ））までの伝搬経路の長さである。
ｉσｔは第ｉのパルス波についての到来時間差標準偏差、ｉσθは第ｉのパルス波についての到来ふ仰角標準偏差（ラジアンで表したもの）である。
また、式（５）の導出に当たっては、これ以降の処理において、尤度の絶対値ではなく、相対的な大きさを問題とするので、式（４）のｐ（Ｐａｔｈｉ）＝１として処理している。
【００５７】
式（５）で表される尤度分布を図９に概念的に示す。
図９では、２つパルス波（第１の到来波、第２の到来波に対応する送出パルス波）が、時刻（先頭のパルス波の送出からの経過時間）ｔに、各々各位置に到達する尤度を等高線で示している。即ち、第１の到来波に対応する送出パルス波の、時刻ｔにおける到達位置が各位置（ｒ，ｚ）である尤度ｐ１（ｒ１（ｔ），ｚ１（ｔ）｜Ｐａｔｈ１）と、第２の到来波に対応する送出パルス波の、時刻ｔにおける到達位置が各位置（ｒ，ｚ）である尤度ｐ２（ｒ２（ｔ），ｚ２（ｔ）｜Ｐａｔｈ２）とが示されている。
【００５８】
上記のように、到来時間差の平均値及び到来ふ仰角の平均値に基づいて求めた位置（尤度の中心）に対して標準偏差に対応する値だけずれた位置が鎖線で示されている。鎖線で示すように、標準偏差に対応する位置は、楕円状である。
【００５９】
上記の処理が第１乃至第Ｎの尤度計算部５３−１〜５３−Ｎの各々で行われ、それぞれの尤度計算部５３−１〜５３−Ｎで求められた尤度Ｌｉ（ｔ）が、統合尤度計算部２０６及び収束時刻判定部２０７に出力される。
【００６０】
分散共分散計算部５４−ｉでは、下記の式（６Ａ）〜（６Ｃ）により分散及び共分散を求める。
【数４】
0005919869.tif 000005
【００６１】
式（６Ａ）〜（６Ｃ）で表される第ｉのパルス波についての分散及び共分散をまとめて符号ＶＣｉで表す。
上記の処理が第１乃至第Ｎの分散共分散計算部５４−１〜５４−Ｎの各々で行われ、それぞれの分散共分散計算部５４−１〜５４−Ｎで求められた分散及び共分散ＶＣｉが、収束時刻計算部２０７に出力される。
【００６２】
統合尤度計算部２０６は、第１乃至第Ｎの尤度計算部５３−１乃至５３−Ｎから出力される第１乃至第Ｎの尤度Ｌ１（ｔ）〜ＬＮ（ｔ）を入力として、尤度を統合し、統合結果を統合尤度Ｌｓ（ｔ）として、収束位置判定部２０８に出力する。
統合尤度計算部２０６では、各時刻ｔにおける第１乃至第Ｎの尤度Ｌ１（ｔ）〜ＬＮ（ｔ）を下記の式（７）により加算することで上記の統合を行う。
【００６３】
【数５】
0005919869.tif 000006
式（７）の計算は、各時刻ｔについてそれぞれ行われ、その結果尤度Ｌｓ（ｔ）の時系列情報
Ｌｓ（ｔ１）、Ｌｓ（ｔ２）、Ｌｓ（ｔ３）、…
が得られる。
【００６４】
収束時刻判定部２０７は、逆伝搬推定部２０５の尤度計算部５３−１〜５３−Ｎで推定された尤度Ｌ１（ｔ）〜ＬＮ（ｔ）の分布間の距離を求めて、該距離が最小となる時刻を推定する。
具体的には、収束時刻判定部２０７は、複数のパルス波から２つを取り出したすべての組合せの各々を構成する２つのパルス波についての尤度の分布間の距離を求め、該距離の、すべての組み合わせについての合計を求めて、複数のパルス波についての尤度の分布間の距離として出力する。
【００６５】
収束時刻判定部２０７は例えば、図１０に示されるように、複数個の分布間距離計算部７１−（１，２）〜７１（Ｎ−１，Ｎ）と、距離加算部７２とを有する。分布間距離計算部７１−（１，２）〜７１（Ｎ−１，Ｎ）の総数Ｕとパルス波の数Ｎとの間には、
Ｕ＝（Ｎ−１）・Ｎ／２
の関係がある。
【００６６】
各分布間距離計算部７１−（ｉ，ｊ）（ｉは、１から（Ｎ−１）までのいずれか、ｊは（ｉ＋１）からＮまでのいずれか、従ってｊ≠ｉ）は、第ｉの尤度計算部５３−ｉから出力される尤度Ｌｉ（ｔ）と、第ｊの尤度計算部５３−ｊから出力される尤度Ｌｊ（ｔ）と、第ｉの分散共分散計算部５４−ｉから出力される分散及び共分散ＶＣｉと、第ｊの分散共分散計算部５４−ｊから出力される分散及び共分散ＶＣｊとを受け、第ｉの送出波の尤度分布と第ｊの送出波の尤度分布の分布間距離ｄ（ｉ，ｊ）を求める。
【００６７】
分布間距離ｄ（ｉ，ｊ）は式（８）によって求められる
【数６】
0005919869.tif 000007
ここで、
【数７】
0005919869.tif 000008
は第ｉの送出波の尤度の中心位置、
【数８】
0005919869.tif 000009
は第ｊの送出波の尤度の中心位置、
Σｉは第ｉの送出波の尤度の分散共分散行列であり、
【数９】
0005919869.tif 000010
で表され、
Σｊは第ｊの送出波の尤度の分散共分散行列であり、
【数１０】
0005919869.tif 000011
で表される。
また、ｔ（μｉ−μｊ）は（μｉ−μｊ）の転置ベクトルである。
【００６８】
式（８）で求められる値は、一般にマハラノビス距離と呼ばれるものである。図９は、このマハラノビス距離を概念的に説明するものである。
【００６９】
図９には、特定の時刻ｔにおける、２つの尤度分布Ｌ１（ｔ）、Ｌ２（ｔ）間のマハラノビス距離が符号ｄ（１，２）で示されている。
【００７０】
分布間距離計算部７１−（ｉ，ｊ）（ｉ＝１〜（Ｎ−１）、ｊ＝（ｉ＋１）〜Ｎ）の各々は、それぞれ入力される尤度の分布間の距離（送出波分布間距離）ｄ（ｉ，ｊ）を計算し、計算結果を距離加算部７２に出力する。
【００７１】
距離加算部７２は、分布間距離計算部７１（１，２）〜７１−（Ｎ−１，Ｎ）で求められた分布間距離ｄ（１，２）〜ｄ（Ｎ−１，Ｎ）を入力として、収束時刻を計算し、計算結果を出力する。即ち、距離加算部７２は、分布間距離ｄ（１，２）〜ｄ（Ｎ−１，Ｎ）を入力とし、式（１１）によって分布間距離の総和を求める。
【数１１】
0005919869.tif 000012
【００７２】
分布間距離の総和Ｄｓ（ｔ）は逆伝搬時間（先頭の送出波の送出からの経過時間）ｔの関数である。この分布間距離の総和Ｄｓ（ｔ）を用い、収束時刻ｔｃを式（１２）によって求める。
ｔｃ = t such that Ｄｓ(t) → minimum        …（１２）
即ち、Ｄｓ(t)が最小となる時刻が収束時刻ｔｃである。
収束時刻判定部２０７は、求めた収束時刻ｔｃを表す情報を出力する。
【００７３】
収束位置判定部２０８は、統合尤度計算部２０６から尤度関数Ｌｓ（ｔ）を受け、収束時刻判定部２０７から収束時刻ｔｃを受け、時刻ｔｃにおける尤度関数Ｌｓ（ｔ）の中で（異なる（ｒ,ｚ）の値に対する尤度Ｌｓ（ｔ）の分布の中で）、最も尤度Ｌｓ（ｔ）の値が大きくなる（ｒ，ｚ）の値を求め、求めた（ｒ，ｚ）の値で表される位置を収束位置を表すものとして出力する。
【００７４】
図１１に示すように、受波器群ＷＲは距離ｒ＝０ｍ、深度ｚ＝７０ｍにあり、音源ＷＳは距離ｒ＝１５０ｍ、深度ｚ＝１０ｍにあり、海底は７５ｍでフラット、音速分布は一定で１５００ｍ／ｓであるとしてシミュレーションを行った場合の、逆伝搬時間ｔと分布間距離Ｄｉ（ｔ）の総和Ｄｓ（ｔ）の関係は図１２に示す如くである。
図１１の位置関係の場合、送出波の伝搬を開始してから１１３ｍｓｅｃの時点で、音線経路が収束したと判定されるのが正しい。
【００７５】
図１２に示されるように、分布間距離の総和は、すべての送出波の音線（伝搬経路）が互いに最も近い位置にある（最も収束した）時刻に最小となっており、分布間距離の総和に基づいて、最も収束した時刻を正しく判定できることが分かる。
【００７６】
収束位置判定部２０８は、上記のように、統合尤度計算部２０６から出力される統合尤度Ｌｓ（ｔ）及び収束時刻判定部２０７から出力される収束時刻ｔｃを入力として、収束位置を検出し、収束位置を表す情報を出力する。
【００７７】
図１３（ａ）〜（ｃ）に到来時間差の標準偏差と到来ふ仰角の標準偏差に相当する範囲を楕円で表して、該範囲が逆伝搬時間とともにどのように変化するかを示す。図１４は、図１３（ｃ）の、推定音源位置の部分を拡大して示す。図１３（ａ）〜（ｃ）及び図１４も、音源と受波器群の位置関係は上記の図１１に示すごとくであると仮定した場合のシミュレーションの結果を示す。
【００７８】
上記の説明では、逆伝搬推定部２０５に設けられた伝搬計算部５２−１〜５２−Ｎの数に等しい数の到来波から成る到来波列が得られた場合を想定しているが、到来波解析部２０３で得られた到来波列を構成する到来波の数ｎが、伝搬計算部５２−１〜５２−Ｎの数Ｎよりも少ない場合には、伝搬計算部５２−１〜５２−Ｎのうちのｎ個のみを用いて伝搬経路の計算を行い、尤度計算部５３−１〜５３−Ｎのうちのｎ個のみを用いて尤度の計算を行う。同様に、収束時刻判定部２０７では、ｎ個の尤度Ｌ１（ｔ）〜Ｌｎ（ｔ）に基づいて、分布間距離計算部７１−（ｉ,ｊ）〜７１−（Ｎ−１,Ｎ）のうちの一部である分布間距離計算部７１−（ｉ,ｊ）〜７１−（ｎ−１,ｎ）のみを用いて分布間距離ｄ（１，２）〜ｄ（ｎ−１，ｎ）の計算を行う。また、順序反転部５１では、第１乃至第ｎの到来波を反転して第１乃至第ｎのパルス波を生成する。この場合、上記の式（１）の代わりに、下記の式（１３）により送出のタイミングを求める。
Ｇｉ＝Ｆｎ−Ｆｉ      （１３）
ここで、Ｆｎは、第ｎの到来波の到来時間差である。
【００７９】
上記の実施の形態による方法を用いる場合、従来の装置で用いていた到来時間差−観測値のテーブルを用いない。従って、観測した到来波と到来時間差−観測値のテーブルを対応付ける必要がなく、従来で生じていた対応付けの誤りによる測距誤差が生じ得ない。
【００８０】
上記の実施の形態では音源からの音波を受信するために上下に配置された受波器から成る受波器群を用いているが、代わりに受信アレイを用いても良く、要するに、到来時間差と到来ふ仰角を観測する構成であれば良い。
【００８１】
また、１つの受波器群ではなく、複数の受波器群を設けることとしても良い。例えば複数の受波器群を、深度の異なる位置に配置してそれぞれの位置で情報収集を行えば、音源の推定精度を向上させることができる。
【００８２】
また、誤差分散を実測により求めることとしているが、同じ状況の観測値が少なく、標準偏差を算出するに足るサンプル数を確保できない場合は、受信時のＳＮ比を測定して、Ｃｒａｍｅｒ−Ｒａｏの不等式から得られる値を、標準偏差として用いることも可能である。この場合の上記の到来波解析部２０３として、上記のようにして標準偏差を求めるように構成したものを用いる。
【符号の説明】
【００８３】
５１  順序反転部、  ５２−１〜５２−Ｎ  伝搬計算部、  ５３−１〜５３−Ｎ  尤度計算部、  ７１−（１,２）〜７１−（ｎ−１,ｎ）  分布間距離計算部、  ７２  距離加算部、  ２０１ａ、２０１ｂ  受波器、  ２０２ａ、２０２ｂ  到来時刻決定部、  ２０３  到来波解析部、  ２０５  逆伝搬推定部、  ２０６  統合尤度計算部、  ２０７  収束時刻判定部、  ２０８  収束位置判定部。
【図１】
0005919869.tif 000013
【図２】
0005919869.tif 000014
【図３】
0005919869.tif 000015
【図４】
0005919869.tif 000016
【図５】
0005919869.tif 000017
【図６】
0005919869.tif 000018
【図７】
0005919869.tif 000019
【図８】
0005919869.tif 000020
【図９】
0005919869.tif 000021
【図１０】
0005919869.tif 000022
【図１１】
0005919869.tif 000023
【図１２】
0005919869.tif 000024
【図１３】
0005919869.tif 000025
【図１４】
0005919869.tif 000026