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