"受波位置において、共通の鉛直線上に整列して配置された第1及び第2の受波器と、前記鉛直線からずれた位置に配置された第3の受波器とにおける到来波の受波時刻に基づいて音源の位置を推定する音源推定装置において、複数の仮想音源位置の各々について、該仮想音源位置に仮想音源が存在すると仮定したときの、前記受波位置における到来波のうちの直接波と海面1回反射波との到来時間差の理論値、並びに前記受波位置における直接波と海底1回反射波との到来時間差の理論値を格納する到来時間差理論値格納手段と、前記第1及び第2の受波器において観測した到来波のデータに基づいて、目標の音源からの直接波と海面1回反射波及び海底1回反射波との到来時間差の観測値を算出する到来時間差観測値算出手段と、前記海面1回反射波と前記直接波との到来時間差と、前記海面1回反射波について前記到来時間差の理論値との不一致度を第1の不一致度として求め、前記海底1回反射波と前記直接波との到来時間差と、前記海底1回反射波について前記到来時間差の理論値との不一致度を第2の不一致度として求め、前記第1の不一致度と前記第2の不一致度との和に基づいてコスト関数を求めるコスト関数算出手段と、前記コスト関数算出手段が算出した前記コストが最小となる前記仮想音源位置を前記目標の音源の位置と推定する音源位置推定手段とを備え、前記コスト関数算出手段は、前記到来波のうちの、前記直接波、前記海面1回反射波及び前記海底1回反射波の各々の、前記第1乃至第3の受波器における到来時刻観測値の前後関係に基づいて、前記仮想音源位置のうちで音源位置ではあり得ないと判定された位置以外の位置についての前記コスト関数に基づいて前記音源位置の推定を行うことを特徴とする音源位置推定装置。受波位置において、共通の鉛直線上に整列して配置された第1及び第2の受波器と、前記鉛直線からずれた位置に配置された第3の受波器とにおける到来波の受波時刻に基づいて音源の位置を推定する音源推定装置において、複数の仮想音源位置の各々について、該仮想音源位置に仮想音源が存在すると仮定したときの、前記受波位置における最初の到来波と2番目以降の到来波のうちの到来ふ仰角が上向きのものの各々との到来時間差の理論値、並びに前記受波位置における最初の到来波と2番目以降の到来波のうちの到来ふ仰角が下向きのものの各々との到来時間差の理論値を格納する到来時間差理論値格納手段と、前記第1及び第2の受波器において観測した到来波のデータに基づいて、目標の音源からの最初の到来波と2番目以降の到来波との到来時間差の観測値を算出する到来時間差観測値算出手段と、前記2番目以降の到来波の各々の到来ふ仰角が上向きか下向きかを判定し、前記2番目以降の到来波のうち、到来ふ仰角が上向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が上向きの到来波についての前記到来時間差の理論値との不一致度を第1の不一致度として求め、前記2番目以降の到来波のうち、到来ふ仰角が下向きの到来波の各々と、前記最初の到来波との到来時間差と、前記到来ふ仰角が下向きの到来波についての前記到来時間差の理論値との不一致度を第2の不一致度として求め、前記第1の不一致度と前記第2の不一致度との和に基づいてコスト関数を求めるコスト関数算出手段と、前記コスト関数算出手段が算出した前記コストが最小となる前記仮想音源位置を前記目標の音源の位置と推定する音源位置推定手段とを備え、前記コスト関数算出手段は、前記最初の到来波を直接波と判別し、前記ふ仰角上向きの到来波のうちの前記直接波の後に最初に到来したものを海面1回反射波と判別し、前記ふ仰角上向きの到来波のうちの前記直接波の後に最初に到来したものを海底1回反射波と判別し、前記到来波のうちの、前記直接波、前記海面1回反射波及び前記海底1回反射波の各々の、前記第1乃至第3の受波器における到来時刻の前後関係に基づいて、前記仮想音源位置のうちで音源位置ではあり得ないと判定された位置以外の位置についての前記コスト関数に基づいて前記音源位置の推定を行うことを特徴とする音源推定装置。前記コスト関数算出手段は、前記到来波のうちの、直接波、海面1回反射波及び海底1回反射波の各々の、前記第1乃至第3の受波器における到来時刻の前後関係に基づいて、前記仮想音源位置のうちで音源位置ではあり得ないと判定された位置以外の位置についての前記第1及び第2の不一致度を求めることを特徴とする請求項1又は2に記載の音源推定装置。前記コスト関数算出手段は、前記到来波のうちの、直接波、海面1回反射波及び海底1回反射波の各々の、前記第1乃至第3の受波器における到来時刻の前後関係に基づいて、前記仮想音源位置のうちで音源位置ではあり得ないと判定された位置について求められた前記和を除いたものを前記コスト関数とすることを特徴とする請求項1又は2に記載の音源位置推定装置。前記第1の受波器と第2の受波器を結ぶ直線と、前記第1の受波器と第3の受波器を結ぶ直線と、前記第2の受波器と第3の受波器を結ぶ直線とが互いに平行ではないように前記第1、第2及び第3の受波器が配置されていることを特徴とする請求項1又は2に記載の音源位置推定装置。受波器が4つ以上設けられ、そのうちの第1の受波器と第2の受波器とが鉛直線方向に整列されており、前記4つ以上の受波器のうちの任意の2つの受波器を結ぶ直線と、他の任意の2つの受波器を結ぶ直線とが互いに平行でないように前記受波器が配置されていることを特徴とする請求項1又は2に記載の音源位置推定装置。各仮想音源位置からの直接波、海面1回反射波、及び海底1回反射波が、前記第1乃至第3の受波器の各々への到来時刻の理論値を示す情報があらかじめ算出され、到来時刻テーブルとして、記憶されており、前記到来時刻の観測値の前後関係に基づく判定の際に上記到来時刻テーブルが参照されることを特徴とする請求項1又は2に記載の音源位置推定装置。前記コスト関数算出手段が、同じ到来波の、前記第1及び第2の受波器での受波時刻の差に基づいて、当該到来波の到来ふ仰角が上向きか下向きかの判定を行うことを特徴とする請求項1又は2に記載の音源推定装置。前記コスト関数算出手段は、前記海面1回反射波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第1の不一致度として計算し、前記海底1回反射波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第2の不一致度として計算することを特徴とする請求項1に記載の音源位置推定装置。前記コスト関数算出手段は、前記到来ふ仰角が上向きの到来波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第1の不一致度として計算し、前記到来ふ仰角が下向きの到来波についての前記到来時間差観測値と前記到来時間差理論値との差の二乗の総和を前記第2の不一致度として計算することを特徴とする請求項2に記載の音源位置推定装置。しかし、上記の従来の装置を用いた場合、音源と受波器を結ぶ線上のコスト関数の値が、音源の位置の前後のかなり広い範囲にわたり小さな値となるため、音源位置を正確に推定することができないという問題があった。即ち、マルチパス波到来時間差測定器102で求められる到来時間差は観測誤差を含むため、コスト関数の値が小さい場所で大きく推定値がずれる可能性が大きく、このため観測誤差が比較的大きいという問題があった。一例として、図3に示すように音速が一定で水深Zbが一定である場合を考える。図3において、Zsは音源の深度、Zrは受波点の深度、rは受波点WRから音源WSまでの水平距離を示す。到来時間差観測値に誤差がないとした場合のコスト関数を求めた結果を図4に示す。横軸は受波点WRからの水平距離r、縦軸は深度下方向ほど深度が大zを表す。図中の曲線はコスト関数の等高線であり、このコスト関数が最小となる位置を音源位置と推定する。到来時間差観測値に誤差がない場合は、真値を推定することができるが、図中のコスト関数は、受波器と音源とを結ぶ線の方向矢印DMCの方向に沿って、コスト関数の変化が小さいために、到来時間差の観測値に誤差が含まれる場合は、この方向に沿って、真の音源位置から推定結果がずれる可能性がある。経験的に知られている観測誤差を与えて、図3に示した条件でシミュレーションを30回試行したところ、距離の誤差標準偏差は、距離真値の約15と、かなり大きい値となった。このようにコスト関数の特質から音源の位置推定誤差が大きいという問題があった。そこで、本出願人は、到来ふ仰角が上向きのマルチパス波と到来ふ仰角が下向きのマルチパス波について別個に到来時間差の理論値を記憶しておき、受信したマルチパス波の各々について到来ふ仰角が上向きか下向きかの判定を行い、観測された到来時間差と理論値との不一致度を上向きのふ仰角を有するマルチパス波と下向きのふ仰角を有するマルチパス波の各々について別個に算出して、その和をコスト関数として、音源位置の推定を行うことを提案した特願2010-012186号。しかしながら、この方法でも、取得できるマルチパス波の数に少ない場合、特に海面1回反射波及び海底1回反射波のみの場合には、推定誤差を十分に小さくすることができないという問題があった。本発明の音源推定装置は、受波位置において、共通の鉛直線上に整列して配置された第1及び第2の受波器と、前記鉛直線からずれた位置に配置された第3の受波器とにおける到来波の受波時刻に基づいて音源の位置を推定する音源推定装置において、複数の仮想音源位置の各々について、該仮想音源位置に仮想音源が存在すると仮定したときの、前記受波位置における到来波のうちの直接波と海面1回反射波との到来時間差の理論値、並びに前記受波位置における直接波と海底1回反射波との到来時間差の理論値を格納する到来時間差理論値格納手段と、前記第1及び第2の受波器において観測した到来波のデータに基づいて、目標の音源からの直接波と海面1回反射波及び海底1回反射波との到来時間差の観測値を算出する到来時間差観測値算出手段と、前記海面1回反射波と前記直接波との到来時間差と、前記海面1回反射波について前記到来時間差の理論値との不一致度を第1の不一致度として求め、前記海底1回反射波と前記直接波との到来時間差と、前記海底1回反射波について前記到来時間差の理論値との不一致度を第2の不一致度として求め、前記第1の不一致度と前記第2の不一致度との和に基づいてコスト関数を求めるコスト関数算出手段と、前記コスト関数算出手段が算出した前記コストが最小となる前記仮想音源位置を前記目標の音源の位置と推定する音源位置推定手段とを備え、前記コスト関数算出手段は、前記到来波のうちの、前記直接波、前記海面1回反射波及び前記海底1回反射波の各々の、前記第1乃至第3の受波器における到来時刻観測値の前後関係に基づいて、前記仮想音源位置のうちで音源位置ではあり得ないと判定された位置以外の位置についての前記コスト関数に基づいて前記音源位置の推定を行うことを特徴とする。本発明によれば、到来時間差観測値に誤差があり、また取得できるマルチパス波の数が少ない場合にも、音源位置の推定誤差を小さくすることができる。実施の形態1.本発明の実施の形態1の音源位置推定装置を図5に示す。図示の音源位置推定装置は、音源から放射されたパルス状音波の直接波とマルチパス波を受信し、それらの到来時間差に基づいて音源の位置を推定するものであり、音源の推定位置は深度及び受波器からの水平距離で表される。なお、音源の位置する方位上方から見たときの方角を示す情報は別途得られるものとする。実施の形態1の音源位置推定装置は、第1、第2及び第3の受波器201、202及び301を備えている。第1の受波器201及び第2の受波器202は、後述のように、到来波のふ仰角が上向きか下向きかの情報を得るためにも用いられるものであり、図6に示すように、同じ鉛直線VL上の異なる位置に設けられている。このように配置されることにより、第1の受波器201及び第2の受波器202のいずれに、到来波信号がより早く到来したかに基づいて、到来波のふ仰角が上向きか下向きかの判断を行うことができる。第1、第2及び第3の受波器201、202及び301は、同じパルス状音波がいずれの受波器に先に到達したかの判定に用いられるので、図6に示すように、間隔DWRa、DWRb、DWRcを開けて配置されている。間隔DWRaは受信位置において想定される音速の最大値と第1及び第2の受波器201及び202で測定可能な最小時間差の積以上の値に設定される。間隔DWRbは受信位置において想定される音速の最大値と第1及び第3の受波器201及び301で測定可能な最小時間差の積以上の値に設定される。間隔DWRcは受信位置において想定される音速の最大値と第2及び第3の受波器202及び301で測定可能な最小時間差の積以上の値に設定される。第1及び第2の受波器201及び202はそれぞれ上側及び下側に配置されるので、上側の受波器、下側の受波器と呼ばれることもある。第3の受波器301は、第1の受波器201よりも下方で、第2の受波器202よりも上方に配置され、かつ第1の受波器201と第2の受波器202を結ぶ鉛直線VLからずれた位置、例えば音源の方向DWSにずれた位置に配置されている。上記のように、音源の方位は別途検出されるので、検出された方向に位置するように、第3の受波器301を、鉛直線VLを中心として回転させることで、第3の受波器301を音源の方向DWSに位置させることができる。また、音源の方向を検出するために、鉛直線を中心軸として円環或いは円筒面に沿って、或いは多数の方向に配置された多数の受波器を用いる場合には、上記多数の受波器のうちの、検出された音源の方向に位置する受波器を第3の受波器301として用いることとしても良い。なお、第3の受波器301は正確に音源の方向にずれた位置になくても良く、ずれが音源の方向への正又は負の成分を有すれば良い。即ち、音源の方向に対して直角な方向を除き、他の如何なる方向へずれていても良い。以下では、説明を簡略にするため、音源の方向にずれているものとする。また、図6に示す例では、第3の受波器301は、第1の受波器201と第2の受波器の中点を通り、鉛直線VLに対して水平な面BL1内にあるものとする。さらにまた、第1の受波器201と第2の受波器202を結ぶ直線VLと、第1の受波器201と第3の受波器301を結ぶ直線Lacと、第2の受波器202と第3の受波器301を結ぶ直線Lbcは互いに平行ではないという条件が満たされている。その意義は後に図15~図19を参照してなされる説明から明らかとなろう。第1、第2及び第3の受波器201、202及び301の出力は、入力端子203、204及び302を介してそれぞれ第1、第2及び第3の到来時刻決定部205、206及び303に入力される。第1、第2及び第3の到来時刻決定部205、206及び303は、それぞれ第1、第2及び第3の受波器201、202及び301の出力に基づいて到来時刻を決定乃至検出する。第1、第2及び第3の受波器201、202及び301は、それらで受信される音波をサンプリングすることで得られるサンプル値を表すデータを出力するが、音波をサンプリングすることで得られるサンプル値の絶対値を時間軸上に並べたものの包絡線は、例えば、図7に示す如くとなる。なお、サンプル値の絶対値の代りにサンプル値の自乗値の包絡線を用いる場合もある。第1、第2及び第3の到来時刻決定部205、206及び303は、それぞれ第1、第2及び第3の受波器201、202及び301からのデータの列に基づき到来波の受信時刻到来時刻を判断する。到来時刻の判断には、従来と同様に、レプリカ相関による方法、立ち上がり検出による方法などを用いることができる。これらの方法は、例えば上記特許文献1においてマルチパス波到来時間差測定器の動作の一部として説明されている。図8a、b及びcに、第1、第2及び第3の到来時刻決定部205、206及び303で決定乃至検出された第1、第2及び第3の受波器201、202及び301におけるパルス状波の到来時刻検出タイミングを示す。図8aは、第1の受波器201における検出タイミングを、早い順にTA1、TA2、TA3、...で示し、図8bは、第2の受波器202における検出タイミングを、早い順にTB1、TB2、TB3、...で示し、図8cは、第3の受波器301における検出タイミングを、早い順にTC1、TC2、TC3、...で示す。ノイズを信号と誤って検出したり、検出漏れが発生しなければ、最初の検出タイミングが直接波の検出タイミングであり、2番目以降の検出タイミングがマルチパス波の検出タイミングである。第1、第2及び第3の到来時刻決定部205、206及び303で検出された検出タイミングを示す情報TA、TB及びTCは、コスト関数算出部210に供給される。第1及び第2の到来時刻決定部205及び206で検出された検出タイミングを示す情報TA及びTBはまた、第1及び第2のマルチパス波到来時間差算出部207及び208に供給される。第1のマルチパス波到来時間差算出部207は、第1の到来時刻決定部205の出力に基づいて、第1の受波器201で受信したマルチパス波の各々についての到来時間差1番目の到来波との到来時間差を算出する。即ち、1番目の到来波PA1と2番目の到来波PA21番目のマルチパス波の到来時間差DU2、1番目の到来波PA1と3番目の到来波PA3の到来時間差DU3、一般化すれば、1番目の到来波PA1とn番目の到来波PAnn=2、3、...との到来時間差DUnを算出する。このようにして算出された到来時間差DUnを示すデータの列GDU=DU2、DU3、...はコスト関数算出部210に供給される。同様に、第2のマルチパス波到来時間差算出部208は、第2の到来時刻決定部206の出力に基づいて、第2の受波器202で受信したマルチパス波の各々についての到来時間差1番目の到来波との到来時間差を算出する。即ち、1番目の到来波PB1と2番目の到来波PB21番目のマルチパス波の到来時間差DL2、1番目の到来波PB1と3番目の到来波PB3の到来時間差DL3、一般化すれば、1番目の到来波PB1とn番目の到来波PBnn=2、3、...との到来時間差DLnを算出する。このようにして算出された到来時間差DLnを示すデータの列GDL=DL2、DL3、...はコスト関数算出部210に供給される。第1の到来時刻決定部205と第1のマルチパス波到来時間差算出部207の組合せは、図1のマルチパス波到来時間差測定器102と同様の機能を有し、第2の到来時刻決定部206と第2のマルチパス波到来時間差算出部208の組合せも、図1のマルチパス波到来時間差測定器102と同様の機能を有する。上記のように、第1のマルチパス波到来時間差算出部207と第2のマルチパス波到来時間差算出部208とは、受波位置において観測した到来波の到来時刻検出タイミングTA及びTBに基づいて、推定対象目標の音源からの最初の到来波と2番目以降の到来波の各々との到来時間差の観測値を算出する到来時間差観測値算出手段を構成し、そのうち、第1のマルチパス波到来時間差算出部207が、第1の受波器201への最初の到来波と2番目以降の到来波の各々との到来時間差の観測値を算出し、第2のマルチパス波到来時間差算出部208が、第2の受波器202への最初の到来波と2番目以降の到来波の各々との到来時間差の観測値を算出する。到来時間差テーブル作成部209は、入力端子213を介して供給される仮想音源の設定位置データWSPと、入力端子214を介して供給される音速プロファイルデータSVPと、入力端子を介して供給される受波器位置データWRPに基づいて、到来時間差テーブルを作成し、記憶部209Mに保存する。到来時間差テーブルを作成する際の仮想音源位置は例えば、図9に示すように設定される。図9に示す例では、探索領域を仮想的にメッシュ格子状に区切り、各交点を仮想音源位置p,qp,q=1,1~P,Qとする。なお、p、qはそれぞれ距離r、深度zを表す値そのものではなく、距離、深度に対応するインデックスであり、自然数1~P、1~Qで表される。メッシュ幅仮想音源位置の間隔は、所望の音源位置の推定精度に応じて設定される。メッシュ幅を狭くするほど高い精度で測定を行なうことができるが、データ量、演算量が多くなり、リアルタイムでの処理が困難となる。従って、高精度が要求される場合には、メッシュ幅を狭くし、高速性が重視される場合には、メッシュ幅を広くする。仮想音源位置WSからの伝搬経路は例えば図10に示すように想定される。図10には、簡単のため、5つの伝搬経路WP1~WP5のみが示されている。伝搬経路音線:音のエネルギーの流れが直線的でないのは、水温、塩分濃度などの不均一性のため、音速が一定ではないことによる。一度も反射をせずに受波点WRに達した波は直接波と呼ばれ、海面SAや海底SBで1回以上反射して受波器に達した波はマルチパス波と呼ばれる。マルチパス波のうち、海面SAで1回反射しただけのものは海面1回反射波と呼ばれ、海底で1回反射しただけのものは海底1回反射波と呼ばれる。このような想定に基づいて、仮想音源位置の各々からの到来時間差を推測し、推測値理論値を時間差理論値として求め、記憶部209Mに記憶しておく。到来時間差テーブル作成部209は、仮想音源位置の各々でパルス状音波が発せられたと仮定した場合に、上側の受波器201における第1の到来波直接波の推定到来時刻と、上側の受波器201に上向きのふ仰角でi番目i=1、2...に到来すると推定されるマルチパス波の到来時刻との差到来時間差の推測値乃至理論値EUii=1、2、...の組乃至列GEUEU1、EU2、...EUIを、すべての仮想音源について求め、記憶部209Mに格納し、同様に、下側の受波器202における第1の到来波直接波の推定到来時刻と、下側の受波器202に下向きのふ仰角でj番目j=1、2...に到来すると推定されるマルチパス波の到来時刻との差到来時間差の推測値乃至理論値ELjj=1、2、...の組乃至列GEL=EL1、EL2、...ELJを、すべての仮想音源について求め、記憶部209Mに格納する。このように、記憶部209Mは、到来時間差理論値格納手段としての役割を果たす。なお、Iは、経験的に受波器201で受信可能と推定される、ふ仰角が上向きのマルチパス波の最大数以上の値に、Jは、経験的に受波器202で受信可能と推定される、ふ仰角が下向きのマルチパス波の最大数以上の値に、それぞれ定められる。全ての仮想音源位置p,qp,q=1,1~P,Qの各々についての、上側の受波器201での到来時間差の列GEU、即ちEU1p,q、EU2p,q、...p,q=1,1~P,Qが得られ、これらがテーブル形式のデータとされ、上向きふ仰角到来時間差理論値として到来時間差テーブル作成部209の記憶部209M内に記憶されている。同様に、全ての仮想音源位置p,qの各々についての、下側の受波器202での到来時間差の列GEL、即ちEL1p,q、EL2p,q、...p,q=1,1~P,Qが得られ、これらがテーブル形式のデータとされ、下向きふ仰角到来時間差理論値として到来時間差テーブル作成部209の記憶部209M内に記憶されている。コスト関数算出部210は、第1及び第2の到来時刻決定部205、206の出力到来時刻を示す信号TA及びTBを入力とするとともに、第1及び第2のマルチパス波到来時間差算出部207及び208で得られた到来時間差データ観測値DU及びDLと、到来時間差テーブル作成部210で作成された到来時間差テーブルの値理論値EU及びELを入力とし、これらに基づいてコスト関数を求め、求めたコスト関数を音源位置推定部211に供給する。音源位置推定部211では、コスト関数算出部210で作成したコスト関数が最小となる仮想音源の距離及び深度を音源位置の推定値とし、該推定値を出力端子212から出力する。以下、到来時間差テーブル作成部209の一構成例につき、図11を参照してより詳細に説明する。入力端子213からは、仮想音源の距離rを表すインデックスPと深度zを表すインデックスQの設定値が仮想音源設定位置データWSPとして入力され、入力端子214からは、音速プロファイルデータSVPが入力され、入力端子215からは、受波器位置データWRPが入力される。これらのデータに基づき音線計算処理部703では音線計算を行う。音線計算の結果として、各到来波の到来時間の推測値理論値と到来ふ仰角の推測値理論値が求められる。ふ仰角判定部704では、音線計算処理部703における計算の結果に基づいて、各仮想音源位置からのマルチパス波の各々について到来ふ仰角が上向きか下向きかの判定を行い、判定結果を到来時間差算出部705に供給する。到来時間差算出部705は、音線計算処理部703における計算の結果に基づき、さらにふ仰角判定部704における判定結果に基づき、上向きと判定されたマルチパス波については、上側の受波器201における到来時間差直接波と当該マルチパス波の到来時間の差を算出し、算出した到来時間差推測値を、記憶部209M内のふ仰上向き到来時間差テーブル706に理論値として書き込み、下向きと判断されたマルチパス波については、下側の受波器202における到来時間差直接波と当該マルチパス波の到来時間の差を算出し、算出した到来時間差推測値を、記憶部209M内のふ仰下向き到来時間差テーブル707に理論値として書き込む。このような処理を行う結果、ふ仰上向き到来時間差テーブル706には、各仮想音源位置p,qから発し、上側の受波器201に上側から上向きのふ仰角で到来するマルチパス波の各々についての到来時間差を表すデータの組が、それぞれの仮想音源毎に格納され、ふ仰下向き到来時間差テーブル707には、各仮想音源位置p,qから発し、下側の受波器202に下側から下向きのふ仰角で到来するマルチパス波の各々についての到来時間差を表すデータの組が、それぞれの仮想音源毎に格納される。これらのテーブル706、707に格納されたデータは、後述のように、端子708、709を介してコスト関数算出部210に供給される。次に、コスト関数算出部210の一構成例につき、図12を参照してより詳細に説明する。コスト関数算出部210は、ふ仰角判定部805と、データ抽出部806及び807と、到来波判別部810と、領域フィルタ作成部305と、領域フィルタ適用部306及び307と、ふ仰上向き理論値-観測値差算出部808と、ふ仰下向き理論値-観測値差算出部809と、ふ仰上向き自乗和算出部812と、ふ仰下向き自乗和算出部813と、和算出部814とを有する。入力端子801には、第1の到来時刻決定部205の出力TAが供給され、入力端子802には、第2の到来時刻決定部206の出力TBが供給され、入力端子811には、第3の到来時刻決定部303の出力TCが供給され、入力端子803には、第1のマルチパス波到来時間差算出部207の出力DUが供給され、入力端子804には、第2のマルチパス波到来時間差算出部208の出力DLが供給される。ふ仰角判定部805は、第1及び第2の到来時刻決定部205、206の出力TA及びTBを、それぞれ入力端子801及び802を介して受け、それぞれ受波器201及び202で受信した到来波の各々についてふ仰角が上向きであるか下向きであるかを判定し、判定結果YULを出力する。図8a及びbに示す例では、1番目、2番目、4番目の到来波については、上側の受波器201のほうがより早いタイミングで受信しており、3番目、5番目の到来波については、下側の受波器のほうがより早いタイミングで受信している。これらのことから、1番目、2番目、4番目の到来波は、上方から到来したこと、即ちふ仰角が上向きであること、3番目、5番目の到来波は下方から到来したこと、即ちふ仰角が下向きであることが分かる。データ抽出部806は、ふ仰角判定部805における判定結果YULに基づいて、第1のマルチパス波到来時間差算出部207から出力される到来時間差を表すデータDUのうち、ふ仰角判定部805でふ仰角が上向きと判定されたマルチパス波についての到来時間差を示すデータFUのみを選択乃至抽出して、ふ仰上向き理論値-観測値差算出部808に供給する。データ抽出部807は、ふ仰角判定部805における判定結果YULに基づいて、第2のマルチパス波到来時間差算出部208から出力された到来時間差を表すデータDLのうち、ふ仰角判定部805でふ仰角が下向きと判定されたマルチパス波についての到来時間差を示すデータFLのみを選択乃至抽出して、ふ仰下向き理論値-観測値差算出部809に供給する。ふ仰角判定部805で上向きと判定されたマルチパス波についての、上側の受波器201における到来時間差をFUaa=1、2、...で表し、ふ仰角判定部805で下向きと判定されたマルチパス波についての、下側の受波器202における到来時間差をFLbb=1、2、...で表す。それぞれのマルチパスの伝搬経路が図10に示す通り即ち、推定通りであれば、1番目、2番目、4番目の到来波直接波と第1、第3のマルチパス波が上向きのふ仰角を有し、3番目、5番目の到来波第2、第4のマルチパス波が下向きのふ仰角を有することになり、したがって、FU1=DU2、FU2=DU4、FL1=DL3、FL2=DL5となる。到来時間差FU1、FU2、...がふ仰上向き理論値-観測値差算出部808に供給され、到来時間差FL1、FL2、...がふ仰下向き理論値-観測値差算出部809に供給される。なお、同じ到来波の受波器201及び202における受信のタイミングの差は極めて小さいことから、例えば検出タイミングの差が所定値以下であれば、同じ到来波が受信されたと判断することとしても良い。また、第1の到来時刻決定部205で到来時刻と決定された時点及びその前後にわたるサンプル値の列と、第2の到来時刻決定部206で到来時刻と決定された時点及びその前後にわたるサンプル値の列の相関が高い場合に、これらの到来時刻が同じ到来波の検出タイミングであると判断することとしても良い。到来波判別部810は、第1、第2及び第3の到来時刻決定部205、206及び303による到来時刻を示す情報を接続端子801、802及び811を介して受け、さらに仰角判定部805によるふ仰角の判定の結果YULを受け、第1、第2及び第3の到来時刻決定部205、206及び303で到来時刻が検出された到来波が直接波であるか、海面1回反射波であるか海底1回反射波であるかの判別を行い、判別結果DSAをフィルタ作成部305に供給する。領域フィルタ作成部305は、第1、第2及び第3の到来時刻決定部205、206及び303から出力された第1、第2及び第3の受波器201、202及び301における到来時刻検出タイミング示す情報TA、TB及びTC、並びに入力端子213を介して供給される仮想音源設定位置データWSP、入力端子214を介して供給される音速プロファイルデータSVP、及び入力端子215を介して供給される受波器位置データWRP、並びに入力端子327を介して供給される到来波判別部810による判別結果DSAに基づいて領域フィルタ情報RFLを作成し、領域フィルタ適用部306は、作成された領域フィルタ情報RFLに応じて、ふ仰上向き到来時間差テーブル706から入力端子708を介して供給される理論値EUに対して領域フィルタリングを掛け、領域フィルタ適用部307は、作成された領域フィルタ情報RFLに応じて、ふ仰下向き到来時間差テーブル707から入力端子709を介して供給される理論値ELに対して領域フィルタリングを掛ける。ここで言う領域フィルタリングは、音源の推定位置を制限することを意味する。即ち、領域フィルタイングを掛けることで、受波器201、202及び301における到来波の検出タイミングの前後関係を基に、音源の位置である可能性がない領域以下不可能領域或いは排除領域ということがあるを音源推定結果から排除し、或いは不可能領域内の仮想音源位置が音源推定結果とならないようにしている。例えば仮に水中の音速が一定であるとすれば、直接波については、第1及び第2の受波器201及び202における到来時刻の前後関係に基づいて、音源が、第1の受波器201と第2の受波器202の中点を通る水平面よりも上に位置するか下に位置するかが分かる。そこで、例えば、第1の受波器201における到来時刻観測値が第2の受波器202における到来時刻観測値よりも先であれば、水平面よりも下の領域を不可能領域として排除することができ、第2の受波器202における到来時刻が第1の受波器201における到来時刻よりも先であれば、水平面よりも上の領域を不可能領域として排除することができる。同様に、第1及び第3の受波器201及び301における到来時刻の前後関係に基づく不可能領域、並びに第2及び第3の受波器202及び301における到来時刻の前後関係に基づく不可能領域を求めることができる。水中の音速が一定でないことを考慮した場合には、上記の領域の境界の決定には、仮想音源位置の各々から、それぞれの受波器までの伝播時間を、仮想音源の位置、音速プロファイル及び受波器の位置に基づいて算出し、算出結果と、観測されたそれぞれの受波器での到来時刻の前後関係とに基づく判断をする必要がある。また、到来時刻の検出誤差や算出誤差を考慮して、音源位置である可能性がない領域として排除する領域を若干狭めるのが望ましい。音源位置である可能性のない領域以外は、音源位置である可能性のある領域であり、以下可能領域或いは通過領域と言うことがある。各仮想音源位置の各々が、可能領域内にあるか不可能領域内にあるかを示す情報は、仮想音源位置の各々p,qについての2値情報Fp,q=true/falseで構成されたデータテーブルの形を取る。領域フィルタ適用部306及び307によるフィルタリングは、上記の不可能領域内の仮想音源位置が、推定結果となることを防ぐための処理絞込み処理であり、コスト関数の算出に先立って、不可能領域内の仮想音源位置を排除する排除された仮想音源位置についてはコスト関数の算出をしないようにするか、コスト関数の算出後に不可能領域内の仮想音源位置についてのコスト関数を排除することで、不可能領域内の仮想音源位置が音源位置の推定結果とはならないようにするものであり、領域フィルタ作成部305で作成された、領域フィルタ情報可能領域、不可能領域の別を示す情報RFLに基づき絞込みが行われる。領域フィルタ作成部305については、後にさらに詳しく説明する。ふ仰上向き理論値-観測値差算出部808は、第1のマルチパス波到来時間差算出部207から入力端子803を介して供給される到来時間差の観測値DUのうち、データ抽出部806で抽出されたふ仰上向きのマルチパス波についての到来時間差の観測値の列GFUの各要素FU1、FU2、...と、ふ仰上向き到来時間差テーブル706から入力端子708を介して供給されるふ仰上向き到来時間差理論値EU1p,q、EU2p,q、...p,q=1,1~P,Qのうちのフィルタ適用部306で領域フィルタリングされたもの、即ち、EU1p,q、EU2p,q、...ただし、p,qは、1,1~P,QのうちのFp,q=trueを満たすものの各々の各要素とを受け、両者の対応する要素間の差を算出する。即ち、ふ仰上向き到来時間差理論値のうちの領域フィルタ適用部306を通過したもののi番目の要素をEUip,q、ふ仰上向きのマルチパス波についての到来時間差の観測値のi番目の要素をFUiとすると、次の式で表される計算を行う。△FUip,q=EUip,q-FUi...3ふ仰上向き自乗和算出部812では、式3の演算結果を自乗したものの、すべての、ふ仰角上向きのマルチパス波についての総和下記の式4で表されるを求める。式4で、I’は、I’≦Iであり、ふ仰角が上向きのマルチパス波がI番目まで検出された場合には、I’=Iであり、I番目まで検出されなかった場合には、検出されたふ仰角が上向きのマルチパス波の数を表す。同様に、ふ仰下向き理論値-観測値差算出部809は、第2のマルチパス波到来時間差算出部208から入力端子802を介して供給される到来時間差の観測値のうち、データ抽出部807で抽出されたふ仰下向きのマルチパス波についての到来時間差の観測値の列GFLの各要素FL1、FL2、...と、ふ仰下向き到来時間差テーブル707から入力端子709を介してさらに領域フィルタ適用部306を介して供給されるふ仰下向き到来時間差理論値EL1p,q、EL2p,q、...p,q=1,1~P,Qのうちのフィルタ適用部307で領域フィルタリングされたもの、即ち、EL1p,q、EL2p,q、...ただし、p,qは、1,1~P,QのうちのFp,q=trueを満たすものの各々の各要素とを受け、両者の対応する要素間の差を算出する。即ち、ふ仰下向き到来時間差理論値のうちの領域フィルタ適用部307を通過したもののi番目の要素をELip,q、ふ仰下向きのマルチパス波についての到来時間差の観測値のj番目の要素をELjとすると、次の式で表される計算を行う。△FLjp,q=ELjp,q-FLj...5ふ仰下向き自乗和算出部813では、式5の演算結果を自乗したものの、すべての、ふ仰角下向きのマルチパス波についての総和下記の式6で表されるを求める。式6で、J’は、J’≦Jであり、ふ仰角が下向きのマルチパス波がJ番目まで検出された場合には、J’=Jであり、J番目まで検出されなかった場合には、検出されたふ仰角が下向きのマルチパス波の数を表す。和算出部814では、式4及び式6の演算結果を用いて、次の演算を行い、演算結果を接続端子815から出力する。Kp,q=KUp,qKLp,q...7式7のKp,qが本実施の形態1で用いるコスト関数である。音源位置推定部211では式7で与えられるコスト関数が最小となる座標p,qを求め、求められた座標p,qが推定対象目標の音源位置WSを表すものであると推定する。式4、式6で求められる自乗和はともに、観測値と理論値テーブル値との一致度が大きいほど小さくなる。従って、式4、式6の自乗和の和で与えられる式7のコスト関数の値が小さいほど上記の一致度が大きいことを意味する。従って式7のコスト関数が最小となる座標p,qを求めることは、一致度が最も大きい不一致度が最も小さい座標p,qを求めることを意味する。以上のようにしてコスト関数を求め、該コスト関数が最小となる座標p,qを求めることで、音源位置の推定を正確に行うことができる。次に、領域フィルタ作成部305について説明する。図13は、領域フィルタ作成部305の構成例を示すブロック図である。図示の領域フィルタ作成部305は、直接波領域フィルタ作成部311と、海面反射波領域フィルタ作成部312と、海底反射波領域フィルタ作成部313と、領域フィルタ結合部314と、到来時刻テーブル作成部315とを有する。入力端子321からは、第1の到来時刻決定部205で検出された到来時刻示す情報TAが、直接波領域フィルタ作成部311及び海面反射波領域フィルタ作成部312に入力される。入力端子322からは、第2の到来時刻決定部206で検出された到来時刻示す情報TBが、直接波領域フィルタ作成部311及び海底反射波領域フィルタ作成部313に入力される。入力端子323からは、第3の到来時刻決定部303で検出された到来時刻示す情報TCが、直接波領域フィルタ作成部311、海面反射波領域フィルタ作成部312及び海底反射波領域フィルタ作成部313に入力される。入力端子213、214及び215からは、仮想音源設定位置データWSP、音速プロファイルデータSVP並びに第1、第2及び第3の受波器201、202及び301の位置を示す情報受波器位置データWRPが入力され、到来時刻テーブル作成部315に入力される。到来時刻テーブル作成部315は、記憶部315Mを有し、以下のように到来時刻データを作成し、テーブルの形で記憶部315Mに保持させる。到来時刻テーブル作成部315は、受波器位置データWRP、仮想音源設定位置データWSP、及び音速プロファイルSVPに基づいて、仮想音源位置p,qを発した音による直接波が、それぞれ第1、第2及び第3の受波器201、202及び301に到来する時刻相対時刻、例えば仮想音源で音が発生してからの時間、即ち伝播時間SAdp,q、SBdp,q及びSCdp,q、仮想音源位置p,qを発した音による海面1回反射波が、それぞれ第1及び第3の受波器201及び301に到来する時刻相対時刻、例えば仮想音源で音が発生してからの時間、即ち伝播時間SAap,q及びSCap,q、並びに仮想音源位置p,qを発した音による海底1回反射波が、それぞれ第2及び第3の受波器202及び301に到来する時刻相対時刻、例えば仮想音源で音が発生してからの時間、即ち伝播時間SBbp,q及びSCbp,qを、全ての仮想音源位置11~P,Qについて計算し、計算結果を、記憶部315Mにテーブルの形で記憶させておく。記憶部315Mに保持されている到来時刻テーブルのデータは直接波領域フィルタ作成部311、海面反射波領域フィルタ作成部312及び海底反射波領域フィルタ作成部313により読み出されて、それぞれの領域フィルタ情報の作成に利用される。直接波領域フィルタ作成部311は、第1、第2及び第3の到来時刻決定部205、206及び303から入力端子321、322及び323を介して供給される、第1、第2及び第3の受波器201、202及び301における到来波の検出タイミングを示す情報TA、TB及びTCを受け、さらに到来波判別部810による判別結果DSAを、入力端子327を介して受け、さらに到来時刻テーブル作成部315から到来時刻テーブルのデータSAdpq、SBdpq及びSCdpqを受け、到来検出タイミングを示す情報TA、TB及びTCのうち、判別結果DSAに基づいて直接波の到来検出タイミングTAd、TBd及びTCdを示す情報を識別し或いは抽出し、第1、第2及び第3の受波器201、201及び301における、直接波の到来時刻の前後関係に基づいて直接波領域フィルタ情報RFdを作成する。海面反射波領域フィルタ作成部312は、入力端子321及び323から供給される、第1及び第3の受波器201及び301における、到来波の検出タイミングを示す情報TA及びTC、到来時刻テーブル作成部315からの到来時刻テーブルのデータSAapq及びSCapq、並びに到来波判別部810による判別結果DSAに基づいて、到来検出タイミングを示す情報TA及びTCのうち、判別結果DSAに基づいて海面1回反射波の到来検出タイミングを示す情報TAa及びTCaを識別或いは抽出し、第1及び第3の受波器201及び301における、海面1回反射波の到来時刻の前後関係に基づいて海面反射波領域フィルタ情報RFaを作成する。海底反射波領域フィルタ作成部313は、入力端子322及び323から供給される、第2及び第3の受波器202及び301における、到来波の検出タイミングを示す情報TB及びTC、到来時刻テーブル作成部315からの到来時刻テーブルのデータSBbpq及びSCbpq、並びに到来波判別部810による判別結果DSAに基づいて、到来検出タイミングを示す情報TB及びTCのうち、判別結果DSAに基づいて海底1回反射波の到来検出タイミングを示す情報TBb及びTCbを識別或いは抽出し、第2及び第3の受波器202及び301における海底1回反射波について海底反射波領域フィルタ情報RFbを作成する。領域フィルタ結合部314は、直接波領域フィルタ作成部311で作成された直接波領域フィルタ情報RFd、海面反射波領域フィルタ作成部312で作成された海面反射波領域フィルタ情報RFa及び海底反射波領域フィルタ作成部313で作成された海底反射波領域フィルタ情報RFbを結合して領域フィルタ情報RFLを作成し、出力端子328から出力する。図14は、直接波領域フィルタ作成部311の構成例を示す。図示の直接波領域フィルタ作成部311は、第1の直接波領域フィルタ作成部331、第2の直接波領域フィルタ作成部332、第3の直接波領域フィルタ作成部333、及び直接波領域フィルタ結合部334を有する。第1の直接波領域フィルタ作成部331は、入力端子321及び322を介して供給される、第1及び第2の受波器201及び202における到来波の検出タイミングを示す情報TA及びTB、到来時刻テーブル作成部315からの到来時刻テーブルのデータSAdpq及びSBdpq、並びに到来波判別部810による判別結果DSAを入力とし、第1及び第2の到来時刻決定部205及び206で検出された直接波の到来時刻TAd及びTBdに基づいて第1の直接波領域フィルタ情報RFd1を作成する。まず、第1の直接波領域フィルタ作成部331は、第1の受波器201における直接波の到来時刻TAdと第2の受波器202における直接波の到来時刻TBdの前後関係を判定する。この判定は、到来時刻TAdと到来時刻TBdを比較することにより、即ち、TAd<TBdか否かによって行われる。そして判定結果に基づいて、以下の式8で表される演算を行って、仮想音源位置の各々pqが、可能領域内に位置するか否かを示す情報Fd1pqを生成する。ここで、¬はブール代数における論理否定演算を表し、Fd1p,qは以下の式で定義されるブール値を表す。SAdp,q、SBdp,qは到来時刻テーブル作成部315の記憶部315Mに記憶されている伝播時間の理論値である。上記の情報Fd1pqの集合座標pqのすべての値1,1~P,Qについての集合が第1の直接波領域フィルタ情報RFd1を構成する。仮に、水中における音速が一定である場合には、上記の式8によって定義される領域の境界は図15に線BL1で示されるごとくとなる。図15で線BL1は、第1の受波器201と第2の受波器202の中点鉛直線VL上にあり、第1の受波器201と第2の受波器202の中間に位置する点を通る水平面鉛直線VLに垂直な面を示す。直接波に関し、図8ことを意味する。従って式24のコスト関数が最小となる座標p,qを求めることは、一致度が最も大きい不一致度が最も小さい座標p,qを求めることを意味する。以上のように、領域フィルタ適用部306及び307で絞り込まれた仮想音源位置についてコスト関数を求め、上記コスト関数が最小となる座標p,qを求めることで、音源位置の推定を正確に行うことができる。以下、本発明の効果を説明する。まず、十分な数のマルチパス波が得られる場合の効果を説明する。従来例のコスト関数は、図4中の矢印DMCで示した方向に沿う領域でコスト関数の差が小さく、従って、到来時間差観測値に誤差がある場合に、音源位置の推定誤差が大きいという問題があったが、本実施の形態1では、上記のように、マルチパス波の到来ふ仰角の上向き、下向きのそれぞれについて別個の到来時間差のテーブルを用い、別個の不一致度を表す値CU、CLを求め、これらの和をコスト関数として、音源位置の推定を行なうこととしており、これにより推定誤差を低減することができる。一例として、従来例と同様に、図3に示すように音速が一定で水深の変化が無い場合を考える。従来例と同条件で到来時間差観測値に誤差がないとした場合のコスト関数を、図20に示す。従来例と同様に、横軸は受波点Zrからの水平距離r、縦軸は深度下方向ほど深度が大zを表す。図中の線はコスト関数の等高線である。図示のように、コスト関数の等高線は音源位置真値を中心に楕円状になり、受波器に近い場所はコストが大きくなるように変化している。従って、コスト関数の形状から定性的に改善していることが判る。また、従来例と同様に、経験的に知られている観測誤差を与えて、図3に示した条件でシミュレーションを30回試行したところ、距離の誤差標準偏差は距離真値の約5と従来例と比較して約1/3となり、音源の位置推定誤差に低減に寄与することが確認された。以上のように、図5の構成を用いることで、取得できるパルスの数が十分に多い場合には、到来時間差観測値に誤差がある場合にも、音源位置の推定誤差を小さくすることができる。しかしながら、取得できるパルスの数が少ない場合には、推定誤差を十分に小さくすることができないという問題があった。即ち、先願で開示された構成の場合にも、マルチパス波として海面1回反射波及び海底1回反射波のみが取得できた場合には、図21に示すように、コスト関数を表す等高線が一方向に延びた形状となり、コスト関数は図の矢印DMCの方向に緩やかな変化をするため、音源位置の推定誤差が大きいという問題があった。これに対して、上記のように、マルチパス波の到来ふ仰角の上向き、下向きのそれぞれについて別個の到来時間差のテーブルの理論値のうちの領域フィルタ適用部306及び307により絞り込まれた値を用い、ふ仰角上向き、下向きの各々について別個の不一致度を表す値CU、CLを求め、これらの和をコスト関数として、音源位置の推定を行なうことで、図22に示すように、推定誤差を低減することができる。なお、上記の実施の形態1では、ふ仰角の判定のために2つの受波器における到来時刻の差を用いているが、2つの受波器の代わりに、例えば、鉛直受波器アレイを用いることも可能である。なお、上記の例ではコスト関数の算出に当たり、音源の位置ではあり得ないと判定された仮想音源位置を除外して、除外された領域以外の仮想音源位置についてのみ、コスト関数を計算しているが、このようにする代わりに、一旦すべての仮想音源位置についてコスト関数を計算し、その後で、音源の位置ではあり得ない領域を除外したコスト関数に変換することとしても良い。例えば、コスト関数算出部210を図23に示すように構成し、和算出部814の出力側に領域フィルタ適用部816を挿入し、和算出部814で算出されたコスト関数Kp,qを以下のように変換する。Kfp,qは変換後のコスト関数である。音源位置推定部211は、この変換後のコスト関数が最小となる位置仮想音源位置のいずれかを音源位置と推定する。上記の実施の形態では、領域フィルタリングのために3つの受波器を用いたが、4つ以上の受波器を用いても良い。但し、2つの受波器を鉛直方向に整列させ、3つ目以降の受波器は、任意の2つの受波器を結ぶ直線が互いに平行にならないようにするのが望ましい。上記の例では、直接波のほかにマルチパス波として、海面1回反射波及び海底1回反射波のみが取得できた場合の動作について詳しく説明したが、マルチパス波が3つ以上取得できた場合にも上記と同様に領域フィルタリングによる絞込みを併用して音源位置推定を行うことができる。実施の形態2.上記の実施の形態1の方法では、第1及び第2の受波器201及び202で受信した音波が、信号のみから成る場合には問題ないが、ノイズが誤って信号と検出誤検出されたり、検出されるべき信号が検出されない場合には、問題がある。誤検出の場合には、誤った組合せで受波器間での到来時間差を計算し、誤ったふ仰角の上下判定の結果に基づいてコスト算出を行う可能性があり、コスト算出に誤りが生じ、音源位置推定を誤る可能性がある。検出されるべき信号が検出されない場合には、ふ仰角の判定のために用いるべき一対の信号受信タイミングが揃わず、そのため、到来時間差テーブルとの比較不一致度の算出が適切に行なえないという問題がある。領域フィルタ作成部305による領域フィルタ情報の作成に関しても同様の問題があり、信頼性の高い到来時刻データのみを用いて領域フィルタ情報の作成を行うのが望ましい。本発明の実施の形態2における音源位置推定装置の全体的構成は、図5に示される通りであるが、コスト関数算出部210の構成が、図24のごとくである点で異なる。図24に示されるコスト関数算出部210は、図12に示されるコスト関数算出部210と概して同じであるが、図24に示される部材に加えて、時間差判定部825、及びデータ抽出部858を備え、さらにデータ抽出部806及び807の代わりにデータ抽出部826、827を備えている点で異なる。図25は、時間差判定部825の構成例を示す。図示の時間差判定部825は、到来時間差最小組合せ算出部853と、受波器間最大到来時間差算出部855と、信頼性評価部854とを有する。到来時間差最小組合せ算出部853では、入力端子801及び802を介して供給される第1及び第2の到来時刻決定部205、206の出力に基づき、第1及び第2の受波器201及び202における到来波の到来時刻タイミングのうちで、差が最も小さい組合せ乃至対を求める。その動作を、図26a及びbを用いて説明する。第1の到来時刻決定部205で検出された、上側の受波器201における到来時刻タイミングをtAaa=1,2,...、第2の到来時刻決定部206で検出された下側の受波器202における到来時刻タイミングをtBbb=1,2,...とし、2つの受波器における到来時刻の内、差が最小となる組合せを信頼性の高い組合せ乃至対として求める。例えば、時刻tAaの各々について時刻tBbのうちで最も差の小さいものを信頼性の高い組合せ乃至対として求め、時刻tBbb=1、2、...のうちで、時刻tAaa=1、2、...のいずれについても差が最小と判断されなかったものは、信頼性の低いものとして除外される。また、ある時刻tBbが複数の時刻tAaについて差が最小と判断された場合には、当該時刻tBbとの差がより小さい時刻tAaが選択され、他の時刻tAaは信頼性の低いものとして除外される。図示の例では、時刻tB3は、時刻tA1~tA6のいずれについても、差が最小とはならないため、除外される。一方、時刻tA4と時刻tA5のいずれも、時刻tB5との差が最小であるが、時刻tA5よりも時刻tA4の方が時刻tB5との差が小さいので、時刻tA4が信頼性の高いものとして残され時刻tB5と対を成すものとして選択され、時刻tA5は信頼性の低いものとして除外される。以上の結果、時刻tA1と時刻tB1、時刻tA2と時刻tB2、時刻tA3と時刻tB4、時刻tA4と時刻tB5、時刻tA6と時刻tB6が互いに信頼性の高い対tG1、tG2、tG3、tG4、tG5を成すものとして従って同じ到来波の検出タイミングとして選択され、時刻tA5と時刻tB3はノイズによるもの信頼性の低いものと見なされ除外される。受波器間最大到来時間差算出部855では、入力端子215から入力した受波器位置データWRPに含まれる受波器間距離DWRa、DWRb、DWRcと、入力端子214から入力した音速プロファイルSVPの、受波器201、202及び301の近傍における音速の値Cから、tmab=DWRa/C...26aにより同じ伝搬経路を辿った音波の、第1の受波器201における到来のタイミングと第2の受波器202における到来のタイミングの時間差の理論的最大値を計算する。第1及び第2の受波器201及び202は鉛直方向に整列しているので、到来ふ仰角の絶対値が大きいほど上記の時間差は大きくなる。信頼性評価部854では、組み合わせ算出部853で決めた到来のタイミング時刻の対の時間差観測値dtabが受波器間最大到来時間差tmabを超えているか否か、即ち、dtab>tmabであるか否かを判定する。この判定結果は、第1及び第2の受波器201及び202における到来時刻の対が、同じ到来波の検出タイミングである蓋然性が高い信頼性が高いか否かを示し、dtab≦tmabであれば、上記蓋然性が高いことを意味し、従って、当該到来時刻のデータがふ仰角の判定、到来時間差の比較などに利用するのに適していることを意味し、逆にdtab>tmabであれば、上記蓋然性が低く、当該到来時刻のデータがふ仰角の判定、到来時間差の比較などに利用するのに適さないことを意味する。信頼性評価部954はさらに、第3の受波器301における到来波の検出タイミングtCの各々について、信頼性が高いと判断された第1及び第2の到来波の検出タイミングtA、tBの対に十分に近いかどうかの判断をし、十分に近いと判断されたもののみを信頼性が高いと判断して、到来波判別810及び領域フィルタ作成部305へ供給させる機能を有する。検出タイミングtCの各々について、信頼性が高いと判断された第1及び第2の到来波の検出タイミングtA、tBの対に近いか否かの判断は、例えば、対を構成する到来時刻タイミングtA、tBの各々との差の絶対値の和を、或いは対を構成する到来時刻タイミングtA、tBの各々との差の自乗の和を、対との時刻の差と定義し、この定義された時刻との差が例えば所定値と比べ、あるいは他のものと比較して小さいかどうかによって行っても良い。例えば、図26a~cに示す例において、時刻tAa及び時刻tBbの対tGgg=1、2、...の各々これらの対tGgを構成する時刻tA、tBは、時刻TA、TBとして抽出されるものであるについて、時刻tCcc=1、2、...のうちの最も差の小さいものを信頼性が高いものとして残し、時刻tCcのうちで、時刻tAa及び時刻tBbの対tGgのいずれについても差が最小と判断されなかったものは、信頼性が低いものとして除外される。図示の例では、時刻tC1は、時刻tA1及びtB1の対tG1に近いので信頼性が高いと判断されて、対tG1と組合せられ、時刻tC2は、時刻tA2及びtB2の対tG2に近いので信頼性が高いと判断されて、対tG2と組合せられ、時刻tC3は、時刻tA3及びtB4の対tG3に近いので信頼性が高いと判断されて、対tG3と組合せられるが、時刻tC4は、いずれの対とも近くないので、信頼性が低いと判断される。時刻tC5は、時刻tA4及びtB5の対tG4に近いので信頼性が高いと判断されて、対tG4と組合せられ、時刻tC6は、時刻tA6及びtB6の対tG5に近いので信頼性が高いと判断されて、対tG5と組合せられる。さらに、このようにして信頼性が高いと判断された時刻tCの各々について、組合せにされた受波器201、202における検出時刻との差が受波器201と受波器301、又は受波器202と受波器301について、時間差の理論的最大値以下であるかどうかの確認を行う確認された場合にそれらのデータ検出タイミングは信頼性が高いと判断することとしても良い。第1の受波器201における到来のタイミングと第3の受波器301における到来のタイミングの時間差の理論的最大値tmac、及び第2の受波器202における到来のタイミングと第3の受波器301における到来のタイミングの時間差の理論的最大値tmbcは、上記の式26aと同様に、それぞれtmac=DWRb/C...26btmbc=DWRc/C...26cにより求められる。信頼性評価部854は、時刻tCについての上記の判定結果を、時刻tA、tBについての上記の判定結果と組み合わせて、時刻tA、tB、tCの組が信頼性の高いものであるかどうかを総合的に判定し、その判定結果YVNを出力する。判定結果YVNは、データ抽出部858、806及び807に供給される。データ抽出部858では、入力端子801、802及び811を介して供給される第1、第2及び第3の到来時刻決定部205、206及び303の出力tA、tB、tCのうち、判定結果YVNにより信頼性が高いと認められたものを時刻TA、TB、TCとして抽出し、到来波判別部810、領域フィルタ作成部305、及びふ仰角判定部805に供給する。ふ仰角判定部805では、データ抽出部858から供給された時刻対を表すデータTA、TBについて、ふ仰角が上向きか下向きかの判定を行い、判定結果YULを出力する。判定結果YULは、データ抽出部826及び827、及び到来波判定部810に供給される。図24のデータ抽出部826は、ふ仰角判定部805による判定結果ふ仰角が上向きか下向きかの判定結果YUL及び時間差判定部825による判定結果各検出タイミングtA、tB、tCの信頼性についての判定結果YVNに応じて、第1のマルチパス波到来時間差算出部207の出力DU入力端子803を介して供給されるのうちの、データ抽出部858で抽出された到来時刻に対応し、かつふ仰角判定部805でふ仰角が上向きと判定された到来時間差データ値FUのみを抽出して、ふ仰上向き理論値-観測値差算出部808に供給する。同様に、データ抽出部827は、ふ仰角判定部805による判定結果ふ仰角が上向きか下向きかの判定結果YUL、及び時間差判定部825による判定結果各検出タイミングtA、tB、tCの信頼性についての判定結果YVNに応じて、第2のマルチパス波到来時間差算出部208の出力DL入力端子804を介して供給されるのうちの、データ抽出部858で抽出された到来時刻に対応し、かつふ仰角判定部805でふ仰角が下向きと判定された到来時間差データ値FLのみを抽出して、ふ仰下向き理論値-観測値差算出部809に供給する。図24のふ仰上向き理論値-観測値差算出部808及びふ仰下向き理論値-観測値差算出部809では、データ抽出部826及び827で抽出された到来時間差の観測値FU及びFLを用いて、差分の算出及びそれに基づくコスト関数の算出を行なう。コスト関数算出部210の他の部分の動作は実施の形態1について説明したのと同様である。実施の形態1では雑音の影響の無いほぼ理想的な状況にのみ到来音波のふ仰角の上下を判別したコスト関数を利用することができたが、本実施の形態2で説明したように到来音波のふ仰角の判定を行う信頼性の高いデータのみを用いてふ仰角の判定、到来時間差の比較不一致度の検出を行うことにより、SN比がさほど高くなく到来時間差に誤検出や未検出がある場合にも、到来音波のふ仰上下を判別したコスト関数を適用することが可能となる。このため、SN比がさほど高くなく到来時間差に誤検出や未検出がある場合にもコスト関数を円錐状に保つことができ、音源の位置推定精度を保つことができる。また、領域フィルタ情報の作成にも、信頼性の高いデータのみを用いるので、信頼性の高い領域フィルタ情報を作成することができる。尚、実施の形態2では、ふ仰角の計測に2つの受波器の時間差を用いているが、3つ以上の受波器を用いて、ふ仰角の上下の判定を行うこととしても良い。また、音速プロファイルを測定できなかった場合については、式26a、26b、26cの音速プロファイルの受波器201、202及び301の近傍における音速の値Cの代わりに代表的な音速1500m/sを用いて、tmab=DWRa/1500...27atmac=DWRb/1500...27btmbc=DWRc/1500...27cを利用して到来時間差の理論的最大値を求めることも可能である。"