JP6718098B2 - 位置推定装置及び方法 - Google Patents

Description

本発明は、海中の音源から発せられる音を受動的に受信して音源の位置を推定する位置推定装置及び方法に関する。

一般的に、水中では電磁波の減衰が大きいことから水中物体の探知には音波が利用されている。そして、音波による水中物体を探知する方式としては、アクティブ方式とパッシブ方式が知られている。

アクティブ方式は、反響定位を用い、音源から音を送波して物体から反射するエコー音を検出し、物体までの距離と方位を求めるものであるが、自ら音を出すために自身の存在を暴露してしまうことになる。一方、パッシブ方式は、水中物体の放射音を検出するものであり、自ら音を出さないため、秘密裏に水中物体を探知することができる。

しかしながら、パッシブ方式では、物体の方位を求めるのは容易であるが、物体までの距離を求めるためには目標運動解析やマッチドフィールド処理などによる複雑な計算が必要であった。目標運動解析では連続的な放射音の観測により物体までの距離を推定することができるが、一時的な過渡音では物体までの距離を求めることができない。また、マッチドフィールド処理は海洋環境を事前に把握しておく必要があるが、時間空間的に変化する海洋環境の複雑な特性をデータベース化するのは困難である。

一時的な過渡音で物体までの距離を求める方法として、音波の伝搬時間差を利用する方法が研究されている。このような方法として、下記特許文献１、２に開示される総当たり法と称される方法がある。

下記特許文献１、２の方法は、音源から直接水中を伝搬して到達する音波（直接波）と境界面で反射して到達する音波（反射波）との伝搬時間差を用いて目標位置を求めており、想定する距離や深度の範囲をデカルト座標や極座標において海洋空間全体を格子化し、格子点上の仮目標からの音波の伝搬時間差を総当たり的に計算して実目標からの音波の伝搬時間差（実測値）との差が最小となる距離や深度を目標位置として推定している。従って、総当たり法は、目標位置の推定範囲（深度及び距離）と計算刻み幅で決まる格子点数に対して反復計算を行う必要があるため、確実な手法である一方、推定範囲が広くなるほど計算時間を要するという問題がある。

そこで、上記問題を解決する他の手法として、例えば下記特許文献３に開示される双曲線法がある。この手法は、音源から受波器に直接到達した直接波、海水面で１回反射して受波器に到達した第１反射波、海底面で１回反射して受波器に到達した第２反射波を海水面と海底面との間に設置された受波器で受信し、音源と海水面に対する音源の鏡像とを焦点とする第１双曲線の式と、音源と海底面に対する音源の鏡像とを焦点とする第２双曲線の式から音源の座標である第１双曲線及び第２双曲線の共通する焦点座標を求めている。

特開２００６−１９４６２７号公報 特開２０１１−１４９８６４号公報 特開２０１５−１５２３７６号公報

しかしながら、特許文献３に開示される方法は、一様な環境を仮定しているため、音速分布に深度依存性があったり海底地形に距離依存性があったりするような海洋環境に適用した場合、音源位置が正確に計算できないという問題があった。そのため、音速分布の深度依存性や海底地形の距離依存性がある環境において、短時間で音源位置を推定する新たな手法の確立が求められている。

そこで、本発明は上記問題点に鑑みてなされたものであり、音速分布の深度依存性や海底地形の距離依存性の有無に拘わらず短時間で音源の位置推定が可能な位置推定装置及び方法を提供することを目的とする。

上記した目的を達成するため、本発明の第１の態様は、音速が深度方向に変化し、且つ海底が平坦な既知の海洋中に存在する音源の座標を推定する位置推定装置において、
前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波、海面反射波及び海底反射波の各音線パラメータを求める処理と、
摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求める処理と、
前記摂動量から推定される前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定する処理と、
を行う演算処理部を備えたことを特徴とする、位置推定装置である。

本発明の第２の態様は、音速が深度方向に変化し、且つ海底が起伏する既知の海洋中に存在する音源の座標を推定する位置推定装置において、
前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波と海面反射波の音線パラメータを求める処理と、
前記音源の仮座標と、前記受波器から海底反射点までの仮想水平距離から海底反射波の音線パラメータを求める処理と、
摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求める処理と、
前記摂動量から推定される前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定する処理と、
を行う演算処理部を備えたことを特徴とする、位置推定装置である。

本発明の第３の態様は、本発明の第１又は第２の態様に係る位置推定装置において、前記音線パラメータは、前記受波器の座標を始点とし、前記音源の仮座標を用いて非線形方程式の数値解法により求めることを特徴とする、位置推定装置である。

本発明の第４の態様は、本発明の第１〜第３の何れかの態様に係る位置推定装置において、前記直接波の音線パラメータが実数解を持たないときは、前記直接波の音線が転回点を経由していると仮定して転回点深度に基づき前記直接波の音線パラメータを求めて前記音源の座標を推定することを特徴とする、位置推定装置である。

本発明の第５の態様は、本発明の第２〜第４の何れかの態様に係る位置推定装置において、既知の海底地形座標と、前記受波器から海底反射点までの仮想水平距離とに基づき前記海底反射波の海底反射点座標と前記海底反射波の音線パラメータを求めることを特徴とする、位置推定装置である。

本発明の第６の態様は、音速が深度方向に変化し、且つ海底が平坦な既知の海洋中に存在する音源の座標を推定する位置推定方法において、
前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波、海面反射波及び海底反射波の各音線パラメータを求めるステップと、
摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求めるステップと、
前記摂動量から推定される前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定するステップと、
を含むことを特徴とする、位置推定方法である。

本発明の第７の態様は、音速が深度方向に変化し、且つ海底が起伏する既知の海洋中に存在する音源の座標を推定する位置推定方法において、
前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波と海面反射波の音線パラメータを求めるステップと、
前記音源の仮座標と、前記受波器から海底反射点までの仮想水平距離から海底反射波の音線パラメータを求めるステップと、
摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求めるステップと、
前記摂動量から推定される前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定するステップと、
を含むことを特徴とする、位置推定方法である。

本発明の第８の態様は、本発明の第６又は第７の態様に係る位置推定方法において、前記音線パラメータは、前記受波器の座標を始点とし、前記音源の仮座標を用いて非線形方程式の数値解法により求めることを特徴とする、位置推定方法である。

本発明の第９の態様は、本発明の第６〜第８の何れかの態様に係る位置推定方法において、前記直接波の音線パラメータが実数解を持たないときは、前記直接波の音線が転回点を経由していると仮定して転回点深度に基づき前記直接波の音線パラメータを求めることを特徴とする、位置推定方法である。

本発明の第１０の態様は、本発明の第６〜第８の何れかの態様に係る位置推定方法において、既知の海底地形座標と、前記受波器から海底反射点までの仮想水平距離とに基づき前記海底反射波の海底反射点座標と前記海底反射波の音線パラメータを求めることを特徴とする、位置推定方法である。

本発明によれば、音速分布の深度依存性や海底地形の距離依存性がある環境において、船舶や大型海洋生物海中移動体が発するパルス波を観測して位置を推定しながら行動を把握するエコーロケーションなどに適用すれば、従来の総当たり法のような膨大な演算処理を必要とせず短時間で音源位置を高精度に推定することができる。

本発明の位置推定装置における位置推定方法で用いる固有音線の一例を示す概要図である。 音速の深度分布と伝搬音の音線との関係を示す模式図である。 本発明の位置推定装置の電気的構成を示す概略機能ブロック図である。 第１実施形態の位置推定装置の位置推定処理の流れを示すフローチャートである。 転回点の一例を示した図である。 起伏のある海底面への音波の入射俯仰角と反射俯仰角の関係を示した模式図である。 音源Ｓから発せられたパルス波の音線と、音線のシャドーゾーンの例を示した模式図である。

以下、本発明を実施するための形態について、添付した図面を参照しながら詳細に説明する。また、この実施の形態によりこの発明が限定されるものではなく、この形態に基づいて当業者などによりなされる実施可能な他の形態、実施例及び運用技術などは全て本発明の範疇に含まれるものとする。

［本発明の概要］
図１は、本発明の位置推定装置１で用いる固有音線の一例を示す概要図である。 図１において、高さ方向（深さ方向）をｚ軸で表し、水平方向の距離をｒ軸で表す。以下では、図１に示すように、受波器Ｒが座標（ｒ₀，ｚ₀）にあるときに、音源Ｓの座標（ｒ，ｚ）を求める例である。

通常、音源Ｓから音波（パルス波）が放射されると、受波器Ｒには、受波器Ｒと音源Ｓとを結ぶ直線上を進んだ直接波ｄｗが到達し、次に、海水面Ｗ１で１回反射した海面反射波ｓｗと海底面Ｗ２で１回反射した海底反射波ｂｗが受波器Ｒに到達し、この後、海水面Ｗ１及び海底面Ｗ２での多重反射波が受波器Ｒに到達する。

本発明は、音速分布の深度依存性や海底地形の距離依存性がある海洋環境において、深度方向のみを格子化して海中の音源Ｓから音速分布の深度依存性や海底地形の距離依存性の有無に拘わらず短時間で音源Ｓの位置を推定するため、海中の未知の音源Ｓからのパルス波である直接波ｄｗ、海面反射波ｓｗ及び海底反射波ｂｗを受波器Ｒで受信し、各波の音線パラメータ用いて音源Ｓの位置（座標）を推定する。

［第１実施形態］
次に、本発明の第１実施形態について説明する。
第１実施形態では、音速は深度方向に変化し、海洋環境として海水面Ｗ１及び海底面Ｗ２は凹凸のない平面であり平行とする。また、海水面と海底面との間は、音波を伝える媒質（海水）で満たされているものとする。

図２は、音速の深度分布と伝搬音の音線との関係を示す図であり、海面位置を原点Ｏとしたときに高さ方向（深さ方向）をｚで表し、水平方向について、左図の図２（ａ）は音速ｃ（ｚ）を表し、右図の図２（ｂ）は距離ｒを表している。

一の音線に沿って、水平を基準とした俯仰角θ₀で深度ｚ_０から俯仰角θ（ｚ）となる深度ｚの点まで伝搬する音については、Ｓｎｅｌｌの法則により下記数１が成り立ち、ｐは、一の音線経路上で値が不変となることから、この値ｐを音線パラメータとし、このｐを使った表示を媒介変数表示とする。

また、このときの音線経路長の水平成分は距離ｒであり、媒介変数表示で下記数２として表すことができる。なお、この媒介変数表示による右辺の汎関数をＲ（ｚ₀，ｚ，ｐ）と表す。

ここで、音源Ｓの座標を推定するにあたり、水深ｈを一定、受波器Ｒ（受波点）の座標（ｒ₀，ｚ₀）を既知、音源Ｓの座標（ｒ，ｚ）を未知と仮定して音源Ｓと受波点の２点を通る固有音線を計算すると、相反定理から受波器Ｒから音源Ｓへ至る音線経路と、音源Ｓから受波器Ｒへ至る音線経路は同一になる。そのため、第１実施形態では受波点である受波器Ｒを始点として音源Ｓの位置座標を推定する処理を行う。

次に、第１実施形態に係る位置推定装置１のシステム構成について説明する。
図３に示すように、位置推定装置１は、音源Ｓから放射されたパルス波（直接波ｄｗ、海面反射波ｓｗ、海底反射波ｂｗ）を受信する受波器Ｒと、受波器Ｒが受信したパルス波による位置推定処理に必要な各種情報を入力するデータ入力部１０と、受波器Ｒが受信した各パルス波から得られる情報とデータ入力部１０で入力された各種情報とに基づき音源Ｓの推定される座標を演算する演算処理部２０とで概略構成される。

＜データ入力部＞
データ入力部１０は、音源Ｓの位置を推定するために必要な初期値や実測値などの各種データを入力する。このデータとしては、例えば音速の深度依存性を示す音速分布ｃ（ｚ）、受波器Ｒの座標（ｒ₀，ｚ₀ ）、受波器Ｒが受信した直接波ｄｗと海面反射波ｓｗとの伝搬時間差実測値（τ₁₂ ^(ob)）及び直接波ｄｗと海底反射波ｂｗとの伝搬時間差実測値（τ₁₃ ^(ob)）、例えば測定対象となる海洋における海底地形を所定の大きさのメッシュ状に分割してそれぞれのメッシュ位置の緯度・経度と各メッシュ位置の水深データとが関連付けされた海洋地形座標を含む海洋環境データと、音源Ｓの仮座標などである。
なお、仮座標とは、音源Ｓの推定座標を計算する際に設定される座標であり、推定処理が初回の場合は予め実験などで求められる音源Ｓの大凡の位置を示す初期値座標（ｒ⁽⁰⁾，ｚ⁽⁰⁾）が設定される。

＜演算処理部＞
演算処理部２０は、音線パラメータ計算部２１と、伝搬時間差計算部２２と、摂動計算部２３と、座標計算部２４と、判定出力部２５とを備えている。なお、音源Ｓから受波器Ｒに届く直接波ｄｗ、海面反射波ｓｗ、海底反射波ｂｗは、各パルス波の到達順序及び受波角度に基づいてパルス波の種類を同定する。

音線パラメータ計算部２１は、音源Ｓから放射されたパルス波（直接波ｄｗ、海面反射波ｓｗ、海底反射波ｂｗ）の音線パラメータｐ（ｐ₁，ｐ₂ ，ｐ₃ ）をそれぞれ計算する。ここで、「ｐ₁ 」は直接波ｄｗの音線パラメータ、「ｐ₂ 」は海面反射波ｓｗの音線パラメータ、「ｐ₃ 」は海底反射波ｂｗの音線パラメータを表す。

音線パラメータ計算部２１は、受波器Ｒを始点として直接波ｄｗが区間［ｚ₀ ，ｚ］内にある場合に、仮座標として音源Ｓの初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾ ）と上記数２から直接波ｄｗのｒ座標は下記数３として与えられる。

この汎関数は、非線形であるが、直接波ｄｗの音線パラメータｐ₁ の単調増加関数である。従って、ｒ，ｚの値に仮座標として音源Ｓの初期値座標を与えることで上記数３はｐ₁の関数となり、ニュートン法に代表される非線形方程式の数値解法により直接波ｄｗの音線パラメータｐ₁ が得られる。

また、直接波ｄｗと同様に、音源Ｓからの海面反射波ｓｗのｒ座標は下記数４、海底反射波ｂｗのｒ座標は下記数５として与えられる。

これらの汎関数も非線形であるが、音線パラメータｐ₂ ，ｐ₃ の単調増加関数であり、ｒ，ｚの値に仮座標として音源Ｓの初期値座標を与えることで上記数４、数５はｐ₂ ，ｐ₃ の関数となり、ニュートン法に代表される非線形方程式の数値解法により海面反射波ｓｗと海底反射波ｂｗの音線パラメータｐ₂ ，ｐ₃ がそれぞれ得られる。

なお、音線パラメータ計算部２１で計算される音線パラメータのうち、位置推定処理を開始した直後は、仮座標として予め設定した初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾ ）と上記数３〜５を適宜用いて各波の音線パラメータの初期値（ｐ₁ ⁽⁰⁾，ｐ₂ ⁽⁰⁾，ｐ₃ ⁽⁰⁾）を計算し、座標計算部２４で音源Ｓの推定座標が計算されると、この計算された音源Ｓの推定座標を仮座標として更新して使用する。つまり、推定処理が初回（１回目）の場合は予め与えられた初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾）を使用し、２回目以降は、座標計算部２４で計算した音源Ｓの推定座標（ｒ^(n-1) ，ｚ^(n-1) ）を使用する。ここで、ｎは推定処理の回数である。また、海面反射波ｓｗの音線パラメータｐ₂ が実数解を持たないときは、海面反射波ｓｗが存在しないものとして本装置の処理である音源Ｓの座標推定処理を終了する。

伝搬時間差計算部２２は、直接波ｄｗ、海面反射波ｓｗ及び海底反射波ｂｗの伝搬時間を計算するとともに、直接波ｄｗと海面反射波ｓｗとの伝搬時間差計算値τ₁₂と、直接波ｄｗと海底反射波ｂｗの伝搬時間差計算値τ₁₃を計算する。

音線に沿った伝搬経路ｓにおける深度ｚ₀ からｚまでの伝搬時間をｔとし、これを媒介変数表示として下記数６として表され、右辺の汎関数をＴ（ｚ₀ ，ｚ，ｐ）と表す。

そして、直接波ｄｗの伝搬時間Ｔ₁ 、海面反射波ｓｗの伝搬時間Ｔ₂ 、海底反射波ｂｗの伝搬時間Ｔ₃ はそれぞれ下記数７〜９として得られ、これら伝搬時間の差として直接波ｄｗと海面反射波ｓｗの伝搬時間差τ₁₂（＝Ｔ₂ −Ｔ₁ ）、直接波ｄｗと海底反射波ｂｗの伝搬時間差τ₁₃（＝Ｔ₃ −Ｔ₁ ）をそれぞれ伝搬時間差計算値として得る。

なお、推定処理が初回の場合は、仮座標である音源Ｓの初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾ ）から得られた各波の音線パラメータの初期値（ｐ₁ ⁽⁰⁾，ｐ₂ ⁽⁰⁾，ｐ₃ ⁽⁰⁾）に基づいて各波の伝搬時間差計算値の初期値(τ₁₂ ⁽⁰⁾,τ₁₃ ⁽⁰⁾)を計算する。また、２回目以降は、得られた推定座標（ｒ^(n-1) ，ｚ^(n-1) ）を仮座標として各波の伝搬時間差計算値 （τ₁₂ ^(n-1)，τ₁₃ ^(n-1)）を計算する。

摂動計算部２３は、摂動法により音源Ｓの推定座標の摂動量を計算する。上記５つの変数（ｒ，ｚ，ｐ₁ ，ｐ₂ ，ｐ₃ ）に対して５つの方程式があるため、原理的には方程式の解を求めることは可能であるが、一連の方程式は非線形であるため解析的に解を求めるのは困難となる。また、複雑な非線形方程式系の解を求めるため、Simulated annealing 法やGenetic Algorithm 法などの確率論的解法もあるが、これらの方法では一般に計算時間を要する上に解の収束判定の難しさがある。
このため、本発明では、摂動法により非線形方程式の解を求めている。

変数であるｒ，ｚ，ｐ₁ ，ｐ₂ ，ｐ₃ を微小変化させたときの変分量をそれぞれδr，δz，δp₁ ，δp₂ ，δp₃とし、テーラー展開の１次項までの摂動を考え上記数３〜５及び数７〜９に適用する。変分量δr，δz，δp₁ ，δp₂ ，δp₃のうち独立変数は２つであり、δp₁ ，δp₂ ，δp₃はδr，δzで表すことができるため、下記数１０に示す代数方程式が得られる。

ここで、数１０の「δr」は下記数１１に示すように変分量δrとδzを成分とする列ベクトルであり、Ｇは２行２列の行列で汎関数Ｒ及びＴのｚ，ｐ_i （ｉ＝１，２，３）に関する偏微分係数から構成される。また、数１０の「δτ」は、観測値である伝搬時間差実測値と、伝搬時間差計算部２２で計算された計算値である伝搬時間差計算値との差とする。例えば推定処理が初回の場合、下記数１２に示すような観測値である伝搬時間差実測値（τ₁₂ ^(ob)，τ₁₃ ^(ob)）と計算値である伝搬時間差計算値の初期値(τ₁₂ ⁽⁰⁾,τ₁₃ ⁽⁰⁾)との差となる。

また、列ベクトルは１次独立であるので、Ｇの逆行列が存在し、下記数１３を得る。

つまり、摂動計算部２３は、伝搬時間差計算部２２で計算した伝搬時間差計算値と、データ入力部１０で入力された直接波ｄｗと海面反射波ｓｗとの伝搬時間差実測値（τ₁₂ ^(ob)）及び直接波ｄｗと海底反射波ｂｗとの伝搬時間差実測値（τ₁₃ ^(ob)）との差に基づき、汎関数Ｒ，Ｔの偏微分係数から上記数１３の行列を計算して逆行列を求め、上記数１３から音源Ｓの推定座標の摂動量を計算する。

座標計算部２４は、摂動計算部２３で計算された摂動量と、下記数１４を用いて音源Ｓの推定座標（ｒ⁽ⁿ⁾ ，ｚ⁽ⁿ⁾ ）を求める。座標計算部２４で得られた座標は、現在の音源Ｓの位置として推定された推定座標となる。

判定出力部２５は、座標計算部２４で計算された音源Ｓの推定座標が正確な座標であるかを収束判定する。収束判定処理の内容としては、ｎ回目の推定座標を（ｒ⁽ⁿ⁾ ，ｚ⁽ⁿ⁾）とし、ｎ−１回目の値との誤差を計算し、この誤差値が実験などで予め得られた判定用閾値を下回るか否かを判断する。例えば、下記数１５に示すように差の絶対値の（Ｌ2 ノルム）をとって距離を正規化して誤差値を計算し、この誤差値が判定用閾値を下回るか否かを判断する。

判定出力部２５は、求めた誤差値が判定用閾値を超えていると判断すると、再度音源Ｓの推定座標の計算処理を行う指示を出力し、誤差値が判定用閾値を下回ると判断したときは、座標計算部２４で計算した音源Ｓの推定座標を外部の出力先に出力する。ここで、判定出力部２５が行う指示とは、具体的に音線パラメータ計算部２１に対し、座標計算部２４で得られた推定座標（ｒ⁽ⁿ⁾ ，ｚ⁽ⁿ⁾ ）を仮座標として再度計算を行わせることである。

未知の音源Ｓの座標を推定する場合、仮座標を与えた状態で音線パラメータ計算部２１、伝搬時間差計算部２２、摂動計算部２３及び座標計算部２４による計算を行っても得られた推定座標の値の信頼性は低いため、上記数１５の計算を複数回（例えば１０回程度）繰り返し行って誤差値を収束させることで推定値が真値に近付き信頼性が高まる。

なお、判定出力部２５による判定基準として判定用閾値を用いて判断しているが、推定値として信頼性が得られる程度の計算回数（例えば１０回程度）を満足したか否か、或いはこの計算のうち数１５の最小値を与える座標で推定座標の出力の可否を判断してもよい。この計算回数は、得られた推定座標の信頼性が高くなる回数であれば特に限定されない。この計算回数も判定用閾値と同様に実験などによって得られる。

次に、上述した第１実施形態の位置推定装置１の処理動作（音源Ｓの位置推定方法）について図４を参照しながら説明する。
なお、以下の処理は、位置推定装置１によって位置推定対象となる音源Ｓからの最初のパルス波が観測されたときの処理内容である。

図４に示すように、まずデータ入力部１０から推定処理に必要な各種データの入力を行う（ＳＴ１）。ここでは、音速分布ｃ（ｚ）、受波器Ｒの座標（ｒ₀，ｚ₀ ）と、受波器Ｒが受信した直接波ｄｗと海面反射波ｓｗとの伝搬時間差実測値（τ₁₂ ^(ob)）及び直接波ｄｗと海底反射波ｂｗとの伝搬時間差実測値（τ₁₃ ^(ob)）、海洋環境データ、仮座標である音源Ｓの初期値座標（ｒ⁽⁰⁾，ｚ⁽⁰⁾）を入力する。

次に、上記数２〜５の何れかの式を用いて音線パラメータ、すなわち直接波ｄｗの音線パラメータｐ₁、海面反射波ｓｗの音線パラメータｐ₂及び海底反射波ｂｗの音線パラメータｐ₃ を計算する（ＳＴ２）。
推定処理が初回であるため、音源Ｓの位置は未知なため、データ入力部１０で入力した初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾ ）と上記数３〜５を適宜用いて各波の音線パラメータの初期値（ｐ₁ ⁽⁰⁾，ｐ₂ ⁽⁰⁾，ｐ₃ ⁽⁰⁾）を求める。

次に、上記数６〜９の何れかの式を用いて伝搬時間を計算する（ＳＴ３）。
ここでは、伝搬時間差計算部２２が、直接波ｄｗの伝搬時間Ｔ₁，海面反射波ｓｗの伝搬時間Ｔ₂及び海底反射はｂｗの伝搬時間Ｔ₃を計算する。

次に、ＳＴ３で得た各波の伝搬時間（Ｔ₁，Ｔ₂，Ｔ₃）から直接波ｄｗと海面反射波ｓｗとの伝搬時間差である伝搬時間差計算値τ₁₂と、直接波ｄｗと海底反射波ｂｗとの伝搬時間差である伝搬時間差計算値τ₁₃を計算する（ＳＴ４）。
推定処理が初回であるため、音線パラメータの初期値（ｐ₁ ⁽⁰⁾，ｐ₂ ⁽⁰⁾，ｐ₃ ⁽⁰⁾）に基づいて伝搬時間差計算値の初期値τ₁₂ ⁽⁰⁾ ，τ₁₃ ⁽⁰⁾を計算する。

次に、データ入力部１０で入力された伝搬時間差実測値（τ₁₂ ^(ob)）及び伝搬時間差実測値（τ₁₃ ^(ob)）とＳＴ４で計算された伝搬時間差計算値τ₁₂ ⁽⁰⁾ ，τ₁₃ ⁽⁰⁾との差をδτとして置き換えた上記数１０、数１３を用いて音源Ｓの推定座標の摂動量を計算する（ＳＴ５）。

次に、ＳＴ５によって計算された摂動量に基づき上記数１４を用いて音源Ｓの推定座標を計算する（ＳＴ６）。

次に、ｎ回目の推定座標を（ｒ⁽ⁿ⁾ ，ｚ⁽ⁿ⁾）とし、ｎ−１回目の値との誤差を上記数１５に示すように誤差値を計算し、この誤差値が予め得られた判定用閾値を下回るか否かを判断する（ＳＴ７）。
推定処置が初回であるため、ＳＴ６で計算した推定座標（ｒ⁽¹⁾ ，ｚ⁽¹⁾）とデータ入力部１０で入力された音源Ｓの初期値座標（ｒ⁽⁰⁾，ｚ⁽⁰⁾）から誤差値を計算する。

ＳＴ７において、誤差値が判定用閾値を下回る場合は（ＳＴ７−Ｙｅｓ）、計算した推定座標（ｒ⁽¹⁾ ，ｚ⁽¹⁾）を最終的な音源Ｓの座標として外部に出力して処理を終了する。
一方、誤差値が判定用閾値を超えている場合は（ＳＴ７−Ｎｏ）、再度ＳＴ２に戻って２回目の推定処理を実施する。このとき、１回目で計算した推定座標（ｒ⁽¹⁾ ，ｚ⁽¹⁾）を音源Ｓの仮座標として更新して一連の計算処理を行う。また、ＳＴ２〜ＳＴ７の処理ループは、得られた誤差値が判定用閾値を下回るまで仮座標を更新しながら所定回数繰り返される。

なお、ＳＴ２において、直接波ｄｗの音線パラメータｐ₁や海面反射波ｓｗの音線パラメータｐ₂が実数解を持たない場合は、以降の処理を行わずに本処理を終了する。

［第２実施形態］
次に、本発明の第２実施形態について説明する。
なお、以下の説明では、上述した第１実施形態の同一構成については同符号を付してその構成や処理内容についての説明を省略し、異なる構成や処理機能を中心に説明する。

第２実施形態の位置推定装置１は、音線パラメータ計算部２１による計算処理において、直接波ｄｗの音線パラメータｐ₁ が実数解を持たない場合に、区間［ｚ₀ ，ｚ］外の転回点を経由して到達する音線経路が存在する可能性を考慮した構成である。第２実施形態の音線パラメータ計算部２１における直接波ｄｗの音線パラメータｐ₁ の計算処理及び伝搬時間差計算部２２における直接波ｄｗの伝搬時間計算処理が第１実施形態の処理とは異なる。

図５は、転回点の一例を示した模式図である。図示のように、転回点は音線の向きが下から上或いは上から下に変化する音線の深度の極大／極小点をｚ_tpで表し、転回点ではｐｃ（ｚ）＝１を満足する。図５において、（ａ）は転回点の極小点を示し、（ｂ）は極大点を示している。そして、図示のように転回点深度をｚ_tpとすると、音源Ｓのｒ座標は下記数１６で与えられる。

第２実施形態の音線パラメータ計算部２１は、音線パラメータｐ₁ の値が実数解を持たないときに直接波ｄｗが転回点を経由したものと仮定し、上記数３の代わりに上記数１６を用いて直接波ｄｗの音線パラメータｐ_1ｔｐを求める。計算処理を行うにあたり、推定処理が初回であれば初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾）を仮座標とし、２回目以降であれば座標計算部２４で求めた推定座標（ｒ^(n-1) ，ｚ^(n-1) ）を仮座標として用いて計算する。

また、第２実施形態の伝搬時間差計算部２２は、音線パラメータｐ₁ の値が実数解を持たないときに直接波ｄｗが転回点を経由したものと仮定し、上記数７の代わりに下記数１７を用いて直接波ｄｗの伝搬時間Ｔ₁ を求める。

なお、第２実施形態における音線パラメータ計算部２１は、直接波ｄｗが実数解であるか否かの判断を行いながら上記処理を２、３回程度行い、直接波ｄｗの音線パラメータが実数解を持たない場合は、直接波ｄｗが存在しないものとして本装置の処理である音源Ｓの座標推定処理を終了する。

第２実施形態の位置推定装置１の処理動作としては、上述したように音線パラメータ計算部２１の処理において、直接波ｄｗが実数解を持たないときに転回点を経由していると仮定しそれに対応する処理を行っている。
つまり、第２実施形態の処理動作は、使用する数式が多少異なるが基本的に上述した第１実施形態の処理動作と同様であり、ＳＴ２における直接波ｄｗの音線パラメータ処理が異なる。

第２実施形態における音線パラメータ計算処理（ＳＴ２）は、直接波ｄｗの音線パラメータｐ₁が実数解を持たないとき、上記数１６を用いて音線パラメータｐ_1ｔｐを計算する。このとき、推定処理が初回であれば初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾）を仮座標とし、２回目以降であれば座標計算部２４で求めた推定座標（ｒ^(n-1) ，ｚ^(n-1) ）を仮座標として用いて計算する。

そして、転回点経由と仮定された直接波ｄｗの音線パラメータｐ_1ｔｐが実数となれば次のステップ（ＳＴ４）へと進む。また、計算した結果が実数解を持たない場合は、再度直接波ｄｗの音線パラメータｐ_1ｔｐを計算し直す。この処理を２、３回程度繰り返し行った結果、直接波ｄｗの音線パラメータｐ_1ｔｐが実数解を持たないときは、直接波ｄｗが存在しないと判断して以降の処理を終了する。

［第３実施形態］
次に、本発明の第３実施形態について説明する。
なお、以下の説明では、上述した第１，第２実施形態の同一構成については同符号を付してその構成や処理内容についての説明を省略し、異なる構成や処理機能を中心に説明する。

海底に起伏がある場合に海底が平坦なときと海底反射波ｂｗの音線経路が異なる。そのため、第３実施形態では、海底に起伏があり、海底地形が既知で関数ｈ（ｒ）で求められる場合に、海底反射波ｂｗの音線パラメータｐ₃ と伝搬時間Ｔ₃を計算する処理が第１、第２実施形態と異なる。

図６は、起伏のある海底面への音波の入射俯仰角と反射俯仰角の関係を示した模式図である。受波器Ｒを基準とした海底地形は既知であり、関数ｈ（ｒ）で与えられるものとする。このとき、水平面に対して傾斜角α＝ｄｈ／ｄｒの海底面ｈ（ｒ）に俯仰角θ_inで入射する音波の反射俯仰角θ_out は下記数１８で与えられる。

従って、反射の前後で音線パラメータは変化するため、反射後の音線パラメータを上記数１と上記数１８から下記数１９で与えられる。ここで、受波点から海底反射点までのｒ座標は下記数２０で表し、海底反射点における音速を下記数２１で示す。

また、受波点、海底反射点及び音源Ｓを通る海底反射波ｂｗの固有音線について、音源Ｓまでの距離は下記数２２となるため、海底反射波ｂｗの伝搬時間Ｔ₃ は、下記数２３から求められる。

このように、第３実施形態では、海底に起伏がある場合の音線パラメータｐ₃ を決定するため、音源Ｓの仮座標（例えば、初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾））及び受波点から海底反射点までの仮想水平距離（予め実験などで得られる仮想的な距離）に基づき上記数２０を計算することで音線パラメータｐ₃の初期値を求める。また、上記数２２を解くことで、受波点から海底反射点までの水平距離の値を更新し、上記数２０及び数２２の計算を繰り返し行って、誤差が予め実験などによって得られた閾値範囲内に入るか、又は既定の計算回数（１０回程度）に達したときに上記水平距離が決定する。このため、海底反射波ｂｗの音線パラメータｐ₃のみならず、下記数２４に示す海底反射点の座標も求めることができる。

なお、海底反射点は、音線パラメータｐ₃で表される音線と海底面との交点を計算することでも求めることができる。

よって、第３実施形態の音線パラメータ計算部２１は、データ入力部１０による海洋環境データから海底地形が平坦ではなく変化していると判断した場合、上記数５の代わりに上記数２２を用い、推定処理が初回であれば初期値座標（ｒ⁽⁰⁾ ，ｚ⁽⁰⁾）を仮座標とし、２回目以降であれば座標計算部２４で求めた推定座標（ｒ^(n-1) ，ｚ^(n-1) ）を仮座標として音線パラメータｐ₃を計算する。

また、第３実施形態の伝搬時間差計算部２２は、海底地形が平坦ではなく変化していると判断した場合に、上記数９の代わりに上記数２３を用いて海底反射波ｂｗの伝搬時間Ｔ₃ を求める。

なお、第３実施形態の音線パラメータ計算部２１は、音線パラメータｐ₁ の値が実数解を持たないのときは、上述した第２実施形態の処理も実行する。すなわち、音線パラメータｐ₁ の値が実数解を持たないときは、直接波ｄｗが転回点を経由したものと推定し、上記数３の代わりに上記数１６を用いて直接波ｄｗの音線パラメータｐ₁を求める。

また、第３実施形態の伝搬時間差計算部２２は、音線パラメータｐ₁の値が実数解を持たないときは、上述した第２実施形態の処理も実行する。すなわち、転回点を経由した直接波ｄｗの音線パラメータｐ₁ を用い、上記数７の代わりに上記数１７を用いて直接波ｄｗの伝搬時間Ｔ₁を求める。

第３実施形態の位置推定装置１の処理動作としては、海底に起伏がある場合に音源Ｓから海底反射点を経由して受波器Ｒに届く海底反射波ｂｗの音線パラメータｐ₃と伝搬時間Ｔ₃を得るための対応する処理を行っている。

そのため、第３実施形態の処理動作は、使用する数式が多少異なるが基本的に上述した第１実施形態の処理動作と同様であり、ＳＴ２の処理で海底反射波ｂｗの音線パラメータを取得するにあたり、海底地形を考慮した音線パラメータｐ₃の計算処理の他、海底反射点の座標を計算する処理が追加される。

［第４実施形態］
次に、本発明の第４実施形態について説明する。
なお、以下の説明では、上述した第１〜第３実施形態の同一構成については同符号を付してその構成や処理内容についての説明を省略し、異なる構成や処理機能を中心に説明する。

音源Ｓの深度によっては、直接波ｄｗ，海面反射波ｓｗ，海底反射波ｂｗが受信器Ｒに到達しない領域が存在する。図７は、音源Ｓから発せられたパルス波の音線と、音線のシャドーゾーンの例を示した模式図である。
図示のように音源Ｓから発するパルス波は、海中を上下に屈折しながら伝搬するが、パルス波は海中領域全てを屈折しながら伝搬するわけではなく、所々音波の届かない領域（シャドーゾーン）も存在する。つまり、受波器Ｒがシャドーゾーン内にあると、第１〜第３実施形態の実施に必要なパルス波が観測できないため伝搬時間差による音源Ｓの座標を推定することができない。

そこで、第４実施形態では、音源Ｓの深度或いは受波器Ｒから音源Ｓまでの距離が大凡既知の場合、直接波ｄｗや海面反射波ｓｗ，海底反射波ｂｗのうち何れか１つの伝搬俯仰角θを観測することで音源Ｓのｒ座標或いはｚ座標を推定する。上述した第１〜第３実施形態では、音源Ｓの座標（距離，深度）が共に未知の場合の座標推定処理であるが、第４実施形態では音源Ｓの深度或いは受波器Ｒから音源Ｓまでの距離が既知の場合の座標推定処理である。

あるパルス波の伝搬俯仰角の観測値をθとし上記数１から音線パラメータｐは下記数２５で与えられる。ここで、伝搬俯仰角の観測値θは、受波器Ｒを海中に複数並設し、それぞれの受波器Ｒで受けたパルス波の受信時間（位相のずれ）を補正して得られる。

また、第４実施形態では、音源Ｓからのパルス波が直接波ｄｗ，海面反射波ｓｗ，海底反射波ｂｗのどのパルス波であるかを同定するため、受波器Ｒの上方から受信したパルス波を直接波ｄｗ又は海面反射波ｓｗであると同定し、受波器Ｒの下方から受信したパルス波を海底反射波ｂｗと同定する。よって、第４実施形態では、音源Ｓからのパルス波がどの種類の波であるかを同定する必要があるため、本形態では海底反射波ｂｗと同定できた場合に処理が可能となる。

音源Ｓの深度ｚが大凡既知である場合、第４実施形態の音線パラメータ計算部２１では、上記処理によって音線の種類を特定し、海底反射波ｂｗと同定できた場合は、海底反射波ｂｗのｒ座標を求める上記数５に上記数２５を代入して音源Ｓのｒ座標を求める。前記ｒ座標が求まれば自ずと音源Ｓの座標を得ることができる。つまり、受波器Ｒの座標における海底反射波ｂｗの伝搬俯仰角θの実測値と、音速の値とに基づき海底反射波ｂｗの音線パラメータｐ₃を求め、深度座標が既知の音源Ｓまでの距離を決定することで音源Ｓの座標を推定することができる。

音源Ｓまでの距離ｒが大凡既知である場合、第４実施形態の音線パラメータ計算部２１では、上記処理によって音線の種類を特定し、海底反射波ｂｗと同定できた場合は、海底反射波ｂｗのｒ座標を与える上記数５に上記数２５を代入し、ｚ座標について数値的に解くことにより音源Ｓのｚ座標を求める。前記ｚ座標が求まれば自ずと音源Ｓの座標を得ることができる。つまり、受波器Ｒの座標における海底反射波ｂｗの伝搬俯仰角θの実測値と、音速の値とに基づき海底反射波ｂｗの音線パラメータｐ₃を求め、距離座標が既知の音源Ｓの深度を決定することで音源Ｓの座標を推定することができる。

以上説明したように、本発明に係る位置推定装置１は、受波器Ｒを始点として音源Ｓが発するパルス波の直接波ｄｗ、海面反射波ｓｗ及び海底反射波ｂｗを受信して各波の音線パラメータから伝搬時間差を計算し、この伝搬時間差計算と受信した各波における伝搬時間差実測値との差の摂動量を求め、この摂動量に基づき得られた推定座標の誤差値を収束させながら音源Ｓの座標を推定している。

このため、従来のような総当たり法のように海洋空間をデカルト座標や極座標において海洋空間全体を格子化し、微分方程式による数値解法を用いて計算する必要がないため、短時間で音源Ｓの座標を推定することができる。

１…位置推定装置
１０…データ入力部
２０…演算処理部
２１…音線パラメータ計算部
２２…伝搬時間差計算部
２３…摂動計算部
２４…座標計算部
２５…判定出力部
Ｒ…受波器
Ｓ…音源
ｄｗ…直接波
ｓｗ…海面反射波
ｂｗ…海底反射波
W１…海水面
W２…海底面

Claims (10)

1. 音速が深度方向に変化し、且つ海底が平坦な既知の海洋中に存在する音源の座標を推定する位置推定装置において、
   前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波、海面反射波及び海底反射波の各音線パラメータを求める処理と、
   摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求める処理と、
   前記摂動量から推定される前記音源の推定座標を求める処理と、
   今回求めた前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定する処理と、
   を行う演算処理部を備えたことを特徴とする位置推定装置。
2. 音速が深度方向に変化し、且つ海底が起伏する既知の海洋中に存在する音源の座標を推定する位置推定装置において、
   前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波と海面反射波の音線パラメータを求める処理と、
   前記音源の仮座標と、前記受波器から海底反射点までの仮想水平距離から海底反射波の音線パラメータを求める処理と、
   摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求める処理と、
   前記摂動量から推定される前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定する処理と、
   を行う演算処理部を備えたことを特徴とする位置推定装置。
3. 前記音線パラメータは、前記受波器の座標を始点とし、前記音源の仮座標を用いて非線形方程式の数値解法により求めることを特徴とする請求項１又は２記載の位置推定装置。
4. 前記直接波の音線パラメータが実数解を持たないときは、前記直接波の音線が転回点を経由していると仮定して転回点深度に基づき前記直接波の音線パラメータを求めて前記音源の座標を推定することを特徴とする請求項１〜３の何れかに記載の位置推定装置。
5. 既知の海底地形座標と、前記受波器から海底反射点までの仮想水平距離とに基づき前記海底反射波の海底反射点座標と前記海底反射波の音線パラメータを求めることを特徴とする請求項２〜４の何れかに記載の位置推定装置。
6. 音速が深度方向に変化し、且つ海底が平坦な既知の海洋中に存在する音源の座標を推定する位置推定方法において、
   前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波、海面反射波及び海底反射波の各音線パラメータを求めるステップと、
   摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求めるステップと、
   前記摂動量から推定される前記音源の推定座標を求めるステップと、
   今回求めた前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定するステップと、
   を含むことを特徴とする位置推定方法。
7. 音速が深度方向に変化し、且つ海底が起伏する既知の海洋中に存在する音源の座標を推定する位置推定方法において、
   前記音源の仮座標と、前記音源が発するパルス波を受信する受波器の座標の２点を通る前記パルス波の固有音線である直接波と海面反射波の音線パラメータを求めるステップと、
   前記音源の仮座標と、前記受波器から海底反射点までの仮想水平距離から海底反射波の音線パラメータを求めるステップと、
   摂動法により、前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差実測値と、各波の音線パラメータから得られる前記直接波と前記海面反射波及び前記直接波と前記海底反射波の伝搬時間差計算値との差の摂動量を求めるステップと、
   前記摂動量から推定される前記音源の推定座標と前回求めた推定座標との誤差値が予め設定した判定用閾値を下回ったときに、今回求めた推定座標を前記音源の座標として決定するステップと、
   を含むことを特徴とする位置推定方法。
8. 前記音線パラメータは、前記受波器の座標を始点とし、前記音源の仮座標を用いて非線形方程式の数値解法により求めることを特徴とする請求項６又は７記載の位置推定方法。
9. 前記直接波の音線パラメータが実数解を持たないときは、前記直接波の音線が転回点を経由していると仮定して転回点深度に基づき前記直接波の音線パラメータを求めることを特徴とする請求項６〜８の何れかに記載の位置推定方法。
10. 既知の海底地形座標と、前記受波器から海底反射点までの仮想水平距離とに基づき前記海底反射波の海底反射点座標と前記海底反射波の音線パラメータを求めることを特徴とする請求項６〜８の何れかに記載の位置推定方法。

Images