JP2009083528A - 線状被曳航体形状推定方法及び装置 - Google Patents

Abstract

【課題】 誤差の大きな観測値が突発的に生じても形状を推定できる、運動方程式モデルに求められる要求を緩和できる線状被曳航体形状推定方法を提供する。
【解決手段】 本発明は、曳航体によって水中若しくは水面上を曳航される柔軟性を有する線状の被曳航体の形状を推定する線状被曳航体形状推定方法に関する。そして、観測値取込手段が、被曳航体の形状推定に利用する観測値の取り込みに、一定時間の観測値と形状推定結果の誤差を最小にする、弱拘束ラグランジアンを用いた方法を適用していることを特徴とする。
【選択図】 図１

Description

本発明は、線状被曳航体形状推定方法及び装置に関し、例えば、目標体から発信された音響を、海中航走体又は海上航走体に取り付けられた海中を曳航されるケーブル状の受波器アレイで受信し、目標体の位置や速度などの状態量を推定する目標体運動解析方法や目標体運動解析装置に適用し得るものである。

目標体が海中を進行するような場合、複数の受波器（受波センサ（例えば音響センサ））でなる受波器アレイを船舶の底部などに設けて、目標体から発信された音響を捕捉して目標体の状態量を推定する方法のほか、線状被曳航体に設けられた受波器アレイを利用して目標体の状態量を推定する方法も有効である。

図３は、受波器アレイを装備した線状被曳航体の概略構成を示す説明図である。図３において、線状被曳航体１０は、船舶や潜水船が該当する曳航体１１に対し、曳航ケーブル１２を介して接続され、曳航体１１によって曳航される。線状被曳航体１０は、受波器アレイ部１３と複数（図３では２個の場合を示している）の方位・深度センサ部１４−１、１４−２とでなる。線状被曳航体１０は、曳航体１１からの音響が雑音として受波器アレイ部１３に影響を及ぼさないように、曳航体１１から十分な距離だけ離されている。

受波器アレイ部１３における各受波器はそれぞれ、音響を捕捉して、曳航ケーブル１２を介して曳航体１１に与えるものである。各方位・深度センサ部１４−１、１４−２はそれぞれ、方位センサと深度センサとを搭載している部分であり、受波器アレイ部１３の所定箇所の方位と深度を測定し、曳航ケーブル１２を介して曳航体１１に与えるものである。図３では、方位・深度センサ部（１４−１、１４−２）を、受波器アレイ部１３の両端に設けた場合を示しているが、設置箇所や設置数は、図３のものに限定されない。

線状被曳航体１０（特に、受波器アレイ部１３の部分）や曳航ケーブル１２は、一般的には、フレキシブルな材料で作成されている。そのため、曳航体１１が直進している場合には、受波器アレイ部１３は概ね直線になっているが、曳航体１１が変針した場合には受波器アレイ部１３は変形する。変形した場合には、受波器アレイ部１３のゲインが劣化する。すなわち、受波器アレイ部１３における各受波器の位置などが正確に分からないまま、受波器アレイ部１３が受信した信号を処理した場合には、アレイゲインが劣化したと同様な悪影響が出てしまう。

アレイゲインの劣化を防ぐために、変針時の受波器アレイ部１３（線状被曳航体１０）の形状の推定技術が開発されている（非特許文献１及び２参照）。

図４は、線状被曳航体１０（受波器アレイ部１３）の形状を推定し、受波器受信信号を修正する従来の被曳航体形状推定・受波信号修正装置の構成を示すブロック図である。図４の構成要素の多くは、曳航体１１に設けられている。また、複数存在する構成要素も１個のブロックで示している。

図４において、方位データ２１は、図３の方位・深度センサ部１４−ｉの方位センサによって計測された方位のデータである。方位センサとしては、地磁気を測定し方位を計測するものを適用できる。深度データ２２は、図３の方位・深度センサ部１４−ｉの深度センサが計測した深度のデータである。曳航体運動諸元データ２３は、曳航体１１の針路や速力等の各データを示しており、曳航体１１の図示しない駆動制御部などが出力したものである。潮流データ２４は、曳航体１１の図示しない潮流計測部で観測された海洋の潮流を示している。地磁気環境データ２５は、予め記憶されている検出方位の修正用のデータである。地磁気は、真北を指すのではなく方位方向及び俯仰方向に偏差を持っている。この偏差は、方位センサ２１（線状被曳航体１０）の位置により定まるので、偏差のテーブル（地磁気環境データ２５）を持っていれば、方位計の観測値とあわせて真の方位を求めることができる。

アレイ形状推定装置２６は、上述した各種データ２１〜２５を用いて受波器アレイ部１３の形状を推定するものである。この推定方法（従来方法）については後述する。

受波器位置データ２７は、アレイ形状推定装置２６による受波器アレイ部１３の形状の推定結果から、各受波器の位置を求めた結果である。線状の受波器アレイ部１３のどこに受波器が設けられているかの基準データを、当該被曳航体形状推定・受波信号修正装置２０は保持しており、この基準データと、受波器アレイ部１３の形状の推定結果とを照合して、各受波器の位置を得る。受波器位置ベクトル２８は、各受波器の位置を要素としたベクトルであり、例えば、テーブルデータとして格納される。受波器受信データ２９は、受波器アレイ部１３の各受波器が捕捉した受信データ（受信信号）である。整相処理部３０は、受波器位置ベクトルの情報を用いて、受波器受信データの整相処理を行うものである。例えば、受波器アレイ部１３の出力を処理する目標体運動解析装置が前提としている、受波器アレイ部１３の基本形状（直線状や円筒状）に従う受波器受信データになるように整相処理を行うものである。ビーム出力３１は、整相した結果のデータであり、次段の処理装置（目標体運動解析装置）に与えられる。

次に、アレイ形状推定装置２６による受波器アレイ部１３の形状の従来の推定方法を説明する。図５は、非特許文献１及び２に示されている、カルマンフィルタを用いた従来の形状推定方法を模式図として示したものである。図５において、ｔはいずれも時間のインデックスを表している。

観測値導出部４３では、方位・深度に対するセンサ情報の時系列４１と、図４の説明で述べた地磁気の補正根拠のデータである磁場環境情報４２とを合わせて正しい方位と深度との観測値を導出する。曳航体１１の運動諸元である速力と針路のデータ４４や、海洋環境情報で潮流（流向・流速）のデータ４５や、時間を１ステップ戻す演算子Ｚ^−１の操作部４６の出力などから、遷移行列４７が形成される。遷移行列４７は、具体的には柔軟な線状被曳航体１０（以下、場合によっては単に被曳航体と呼ぶ）の運動方程式である。例えば、非特許文献３に、柔軟な被曳航体が中性浮力の場合の運動方程式が示されている。中性浮力を仮定しているのは、運動が水平面上に限定されるという意味である。鉛直方向の運動は、柔軟な被曳航体に働く力に重力と浮力とが加わるだけでその他は同じであるので、ここでは、柔軟な被曳航体が中性浮力で水平方向に変位する例を用いて説明する。

図６は、従来の推定方式での被曳航体１０の変位と座標系を示している。

図６において、被曳航体１０の深度における平面を考える。また、深度が、海面から離れていて、海流が無視できる程度の深さであるとする。被曳航体１０における流速（相対的には曳航速度）Ｕの方向をｘ方向とし、同一深度における直交方向をｙ方向とし、ｙ方向に柔軟な被曳航体１０が変形することを考える。

図７は、柔軟な被曳航体１０の微小長さδｘ当りにかかる力を示している。Ｔ（ｘ）は張力であり、柔軟な被曳航体１０の長さ方向にかかる。Ｆ_Ａは、流体中の柔軟な被曳航体１０が運動（変形）することにより、流体を加速する力を表している。Ｆ_Ｎ及びＦ_Ｔはそれぞれ、柔軟な被曳航体１０の法線方向及び接線方向に働く粘性による力である。柔軟な被曳航体１０は、長さ方向に伸び縮みしないと近似（仮定）して、ｘ方向の力のバランスから、（１）式が得られる。また、ｙ方向について、運動方程式を（２）式で表すことができる。（２）式において、ｍは柔軟な被曳航体１０の単位長さ当りの質量であり、Ｂは曲げ剛性（ｂｅｎｄｉｎｇ ｓｔｉｆｆｎｅｓｓ）である。

（２）式において、さらに、流体加速力Ｆ_Ａは、非特許文献４によると、（３）式で表され、また、粘性による力Ｆ_Ｎ及びＦ_Ｔはそれぞれ、非特許文献５によると、（４）式、(５)式で表される。（３）式において、ρ_０は流体の密度、ａ_Ａは柔軟な被曳航体１０の半径である（被曳航体１０が円筒形状のものと仮定している）。（４）式及び(５)式において、Ｃ_Ｄはｆｏｒｍ ｄｒａｇ係数、Ｃ_Ｎ及びＣＴはそれぞれ法線及び接線方向のｆｒｉｃｔｉｏｎ ｄｒａｇ係数、ｉは柔軟な被曳航体１０の微小部分の傾きである。

ｉが小さいと、ｓｉｎｉについて（６）式の近似式が成り立ち、ｓｉｎ^２ｉが十分小さく無視できるとすると、（４）式及び(５)式はそれぞれ、（７）式及び(８)式に変形することができる。 (８)式を(１)式に適用して整理すると、張力Ｔ（ｘ）を、（９）式で表すことができ、また、以上の全体を整理すると（１０）式が得られる。(９)式において、ｌ_Ａは柔軟な被曳航体１０の長さである。

図５に戻り、観測行列４８は、モデルの推定値と観測値とを結びつける関係式から求められる。観測値方位に対して、モデルの推定値はｄｙ／ｄｘに相当する。時間軸及び空間軸に対して離散系で書くとすると、（１１）式で表される。ここで、１ｏｒ０はｙ_ｔの対応する位置に観測値がある場合は１、ない場合は０をとるものである。ｖ_ｔは観測誤差の標準偏差である。

カルマンゲイン４９は、一連のカルマンフィルタの計算中に求められるものである。カルマンフィルタの計算式（カルマンフィルタの取り扱うシステム）を、(１２)式に示す状態方程式及び（１３）式に示す観測方程式で定義する。(１２)式及び（１３）式において、Ｆ_ｔは状態遷移行列、ｗ_ｔはシステムの誤差、Ｇ_ｔは駆動行列、Ｈ_ｔは観測行列、ｖ_ｔは観測誤差である。

状態遷移行列Ｆ_ｔは、（１０）式を離散化したときに得られる行列である。各誤差に対しては、相互に無相関の不偏ガウシアン白色雑音を仮定する。このような雑音は、（１４）式で表現することができる。（１４）式において、Ｅ｛Ｘ｝はＸの期待値を表しており、左肩のｔは行列（又はベクトル）の転置を表している。状態量ｙ_ｔが、誤差と相関がないと仮定すると、（１５）式及び（１６）式が成り立つ。初期状態ｙ_０は、Ｎ（／ｙ_０，Σ_０）に従うガウス確率ベクトルであるとする（なお、／ｙ_０は平均値を表すものとする）。

このようなシステムの場合、カルマンゲイン、カルマンフィルタ及びカルマンフィルタによる推定値などは、以下の（１７）式〜（２３）式で表される。

以上のカルマンフィルタによっての被曳航体１０の形状推定値が得られ、被曳航体１０中の受波器位置２７が求められる。上記構成の従来方式では、方位、深度情報の形状推定への取り込みにカルマンフィルタを用いている。方位・深度センサの観測値は誤差を持っている。カルマンフィルタは、誤差の大きな観測値が突発的に入った場合には発散する可能性が大きい。海上交通の盛んな領域が曳航領域である場合には、意図しない船舶や潜水艦や漁業機具などから、突発的な音響が発せられることも多い。

これに対して、観測値の取り込み方法として、変分法若しくはアジョイント（Ａｄｊｏｉｎｔ）法（以下では、単に変分法と呼ぶ）と呼ばれる、一定時間の観測値と形状推定結果の誤差を最小にする方法を採用することを、本願発明者は考えてみた。以下、アレイ形状推定装置２６（図４参照）が、変分法を利用した推定方法を適用した場合について説明する。

図７は、アレイ形状推定装置２６が適用可能な、変分法を利用したアレイ形状推定方法の構成を示すブロック図である。

この推定方法でも、センサ情報（方位・深度の時系列）６１、磁場環境情報６２、曳航体船舶運動諸元の時系列６４、潮流の流向・流速に関する海洋環境情報６５が入力情報となる。センサ情報６１及び磁場環境情報６２から、磁場環境情報を反映して補正した方位情報と深度情報の時系列でなる観測値６３を得る。初期時刻形状の初期値６８も入力される。

ラグランジアン（Ｌａｇｒａｎｇｉａｎ）計算部６６は、ラグランジアン（ラグランジアン関数）を計算すると共に、ラグランジアンから運動方程式とアジョイント（Ａｄｊｏｉｎｔ）方程式を生成する。運動方程式時間前進計算部６７は、生成された運動方程式に、その初期値として適当な値である初期時刻形状の初期値６８を用いて時間発展を計算する。運動方程式の時間発展の時系列情報は全て保存しておく。アジョイント方程式時間後退計算部６９は、生成されたアジョイント方程式を用い、運動方程式時間前進計算部６７で得られた形状の時間発展時系列と観測値６３を用いて、時間後退を計算する。収束判定部７０は、ラグランジアンの勾配による収束判定を行い、収束していない場合は、初期時刻形状更新部７１で初期時刻の形状を計算して、運動方程式時間前進計算部６７から処理を繰り返す。収束している場合は、形状推定結果７２が得られ、さらに、受波器位置２７が計算される。

一連の計算を、カルマンフィルタを適用した場合と同様に、柔軟な被曳航体１０が中性浮力の場合で、水平方向の運動方程式について具体的に示す。

方位の推定値をξとし、ｉ番目の観測値をξ^ｉ _ｏｂｓとする。方位推定値ξと変位ｙとの間には、（２４）式に示す関係がある。

ラグランジアン計算部６６は、（２５）式に示すラグランジアンを計算する。積分範囲のΣは時間ｔと空間ｘに対する範囲である。

（２５）式の積分中の第１項は、観測値と推定値との差異の２乗和を、観測値の誤差分散σ^２で割ったものである。δ（ｘ）はδ関数で、ｘ_ｉは観測値のある位置（方位センサがある位置）を表している。第１項は、観測値と推定値との誤差を表す項であるので、観測種別にあわせて適宜変更する必要がある。第２項〜第４項のλ_１及びλ_２がかかっている項は、ラグランジュ（Ｌａｇｒａｎｇｅ）の未定乗数法の拘束条件に該当し、λ_１及びλ_２は、未定係数と呼ばれている。第２項及び第３項の括弧内は運動方程式である。ここで用いている運動方程式は、中性浮力の場合の例であり、柔軟な被曳航体１０の形状を精度良く表せる方程式であれば良い。また、第４項は、観測値の推定値に対応する方位ξと変位ｙとの関係を表す（２４）式が入る。誤差が最小となるような形状では、（２５）式の値は最小値になる。そこで、（２５）式の変分をとった、(２６)式の変分式を利用する。

（２６）式を、δ_ξ、δλ_１、δλ_２及びδ_ｙの時空間微分ではない形になるように部分積分を行う。得られた式に対して、δ_ξ、δλ_１、δλ_２及びδ_ｙの係数が０となるように、ξ、ｙ、λ_１、λ_２を求めることで形状推定結果を得る。
Ｄｏｕｇｌａｓ Ａ．Ｇｒａｙ ｅｔ ａｌ．，"Ｔｏｗｅｄ Ａｒｒａｙ Ｓｈａｐｅ Ｅｓｔｉｍａｔｉｏｎ Ｕｓｉｎｇ Ｋａｌｍａｎ Ｆｉｌｔｅｒｓ−Ｔｈｅｏｒｅｔｉｃａｌ Ｍｏｄｅｌｓ"，ＩＥＥＥ Ｏｃｅａｎ Ｅｎｇ．１８（１９９３）５４３ Ｊｏｈｎ Ｌ．Ｒｉｌｅｙ ａｎｄ Ｄｏｕｇｌａｓ Ａ．Ｇｒａｙ，"Ｔｏｗｅｄ Ａｒｒａｙ Ｓｈａｐｅ Ｅｓｔｉｍａｔｉｏｎ Ｕｓｉｎｇ Ｋａｌｍａｎ Ｆｉｌｔｅｒｓ−Ｅｘｐｅｒｉｍｅｎｔａｌ Ｉｎｖｅｓｔｉｇａｔｉｏｎｓ"，ＩＥＥＥ Ｏｃｅａｎ Ｅｎｇ．１８（１９９３）５７２ Ａ．Ｐ．ＤＯＷＬＩＮＧ，"Ｔｈｅ ｄｙｎａｍｉｃｓ ｏｆ ｔｏｗｅｄ ｆｌｅｘｉｂｌｅ ｃｙｌｉｎｄｅｒｓ Ｐａｒｔ １． Ｎｅｕｔｒａｌｌｙ ｂｕｏｙａｎｔ ｅｌｅｍｅｎｔｓ"，Ｊ．Ｆｌｕｉｄ Ｍｅｃｈ．１８７（１９８８）５０７ Ｍ．Ｊ．ＬＩＧＨＴＩＬＬ，"Ｎｏｔｅ ｏｎ ｔｈｅ ｓｗｉｍｍｉｎｇ ｏｆ ｓｌｅｎｄｅｒ ｆｉｓｈ．"，Ｊ．Ｆｌｕｉｄ Ｍｅｃｈ．９（１９６０）３０５ Ｇ．Ｉ．ＴＡＹＬＯＲ，"Ａｎａｌｙｓｉｓ ｏｆ ｔｈｅ ｓｗｉｍｍｉｎｇ ｏｆ ｌｏｎｇ ａｎｄ ｎａｒｒｏｗ ａｎｉｍａｌｓ．"，Ｐｒｏｃ．Ｒ．Ｓｏｃ．Ｌｏｎｄ．Ａ２１４（１９５２）１５８

上述したように、方位、深度情報の形状推定への取り込みにカルマンフィルタを用いた場合には、方位・深度センサからの観測値として、誤差の大きな観測値が突発的に入ったときには、カルマンフィルタが発散する可能性が大きいという問題があった。

そこで、本願発明者は、方位、深度情報の形状推定への取り込みに変分法を適用することを考えた。このようにした場合には、突発的に入った異常に誤差の大きい観測値に対しても安定して形状を推定できる。しかしながら、ラグランジュの未定乗数法では運動方程式が厳密に合わなければならず、モデルが不完全な場合には、観測値の情報が生かされず、正しい形状推定を行うことができない。

そのため、誤差の大きな観測値が突発的に生じても形状を推定できる、しかも、運動方程式モデルに求められる要求を緩和できる線状被曳航体形状推定方法及び装置が望まれている。

第１の本発明は、曳航体によって水中若しくは水面上を曳航される柔軟性を有する線状の被曳航体の形状を推定する線状被曳航体形状推定方法において、観測値取込手段が、上記被曳航体の形状推定に利用する観測値の取り込みに、一定時間の観測値と形状推定結果の誤差を最小にする、弱拘束ラグランジアンを用いた方法を適用していることを特徴とする。

第２の本発明は、曳航体によって水中若しくは水面上を曳航される柔軟性を有する線状の被曳航体の形状を推定する線状被曳航体形状推定装置において、上記被曳航体の形状推定に利用する観測値の取り込みに、一定時間の観測値と形状推定結果の誤差を最小にする、弱拘束ラグランジアンを用いた方法を適用している観測値取込手段を有することを特徴とする。

本発明の線状被曳航体形状推定方法及び装置によれば、誤差の大きな観測値が突発的に生じても形状を推定できる、しかも、運動方程式モデルに求められる要求を緩和できるようになる。

（Ａ）主たる実施形態
以下、本発明による線状被曳航体形状推定方法及び装置の一実施形態を、図面を参照しながら詳述する。

（Ａ−１）実施形態の構成
この実施形態の線状被曳航体形状推定装置が、形状を推定する対象（線状被曳航体１０）も、上述した図３に示すような曳航状態にあるものである。

図２は、実施形態の線状被曳航体形状推定装置（アレイ形状推定装置）を有する被曳航体形状推定・受波信号修正装置の構成を示すブロック図であり、従来に係る図４との同一、対応部分には同一符号を付して示している。

図２に示すように、実施形態のアレイ形状推定装置（線状被曳航体形状推定装置）２６Ａも、方位データ２１、深度データ２２、曳航体運動諸元データ２３、潮流データ２４及び地磁気環境データ２５を用いて受波器アレイ部１３の形状を推定するものである。

図１は、実施形態のアレイ形状推定装置（線状被曳航体形状推定装置）２６Ａの機能的構成を示すブロック図である。なお、図１は、実施形態のアレイ形状推定方法（線状被曳航体形状推定方法）における処理の流れを示すフローチャートとして見ることもできる。例えば、図示しないセンサや受波器などの測定構成を除けば、プログラム及びその実行構成（ＣＰＵ、メモリ等）によって線状被曳航体形状推定装置が実現されるが、機能的には、図１で表すことができる。また、別の見方をすれば、図１は、線状被曳航体形状推定プログラムの構成を示したものとして見ることもできる。

図１において、実施形態のアレイ形状推定装置２６Ａも、図８に示した構成とほぼ同様な構成を有するが、ラグランジアン計算部６６に代え、弱拘束ラグランジアン計算部６６Ａを適用している点が、図８に示したものと異なっている。

（Ａ−２）実施形態の動作
以下、実施形態のアレイ形状推定装置（線状被曳航体形状推定装置）２６Ａの動作を説明する。

方位・深度のセンサ情報の時系列６１、及び、利用地点に対する地磁気の強度・方位・俯仰のテーブル値を意味する磁場環境情報６２が入力され、これらを処理し、観測値６３として、磁場環境情報を反映して補正した方位情報と深度情報の時系列６３が得られ、弱拘束ラグランジアン計算部６６Ａに与えられる。また、柔軟な被曳航体１０を曳航する曳航体１１の速力及び針路などの曳航体運動諸元の時系列情報６４や、潮流の流向、流速のデータでなる海洋環境情報６５も、弱拘束ラグランジアン計算部６６Ａに与えられる。

弱拘束ラグランジアン計算部６６Ａは、弱拘束ラグランジアン（弱拘束ラグランジアン関数）を生成すると共に、弱拘束ラグランジアンから運動方程式とアジョイント方程式を生成する。弱拘束ラグランジアン計算部６６Ａは、生成した運動方程式を運動方程式時間前進部６７に与え、生成したアジョイント方程式をアジョイント方程式時間後退部６９に与える。

運動方程式時間前進部６７は、生成された運動方程式に、その初期値として、外部から入力された適当な値である初期時刻形状の初期値６８を用いて時間発展を計算する。運動方程式時間前進部６７は、運動方程式の時間発展の時系列情報を全て保存しておく。

アジョイント方程式時間後退部６９は、生成されたアジョイント方程式を用い、運動方程式時間前進部６７で得られた形状の時間発展時系列と観測値６３とを用いて、時間後退を計算する。

収束判定部７０は、弱拘束ラグランジアンの勾配による収束判定を行う。収束していない場合は、初期時刻形状更新部７１が初期時刻の形状を計算し、運動方程式時間前進計算部からの処理を繰り返させる。収束している場合は、形状推定結果７２を求め、受波器位置計算部２７が各受波器の位置を求める。

以下、実施形態のアレイ形状推定装置２６Ａにおける一連の計算処理を、上述の場合（図８）と同様に、柔軟な被曳航体１０が中性浮力の場合とし、水平方向の運動方程式について具体的に説明する。

方位の推定値をξとし、ｉ番目の観測値をξ^ｉ _ｏｂｓとする。方位推定値ξと変位ｙとの間には、（２７）式（（２４）式と同一である）に示す関係がある。

弱拘束ラグランジアン関数計算部６６Ａでは、まず、(２８)式に示す弱拘束ラグランジアン関数を計算する。積分範囲のΣは時間ｔと空間ｘに対する範囲である。

（２８）式の積分処理対象中の第１項は、観測値と推定値との差異の２乗和を観測値の誤差分散σ_１ ^２で割ったものである。δ（ｘ）は、δ関数でｘ_ｉは観測値のある位置（方位センサがある位置）を表している。第１項は、観測値と推定値との誤差を表す項であるので、観測種別にあわせて適宜変更する必要がある。第２項は、運動方程式を、運動方程式の誤差分散σ_２ ^２で割った形をとっている。より具体的に言えば、上述した変分法と呼ばれる方式で説明した運動方程式（左辺−右辺）の２乗を、運動方程式の誤差分散σ_２ ^２で割った形をとっている。第２項及び第３項の未定係数λ_１及びλ_２がかかっている項は、ラグランジュの未定乗数法の拘束条件に該当する。第２項の括弧内は、運動方程式である。ここで用いている運動方程式は、中性浮力の場合の例であり、柔軟な被曳航体１０の形状を精度良く表せる運動方程式であれば他の形式の運動方程式であっても良い。また、第３項は、観測値の推定値に対応する方位ξと変位ｙとの関係を表すように式に導入されたものである。

上述した変分法と呼ばれる方式と異なるのは、運動方程式が、運動方程式の２乗を誤差分散で割った形で導入されていることである。運動方程式の２乗を誤差分散で割った形にしたことにより、運動方程式にかかる項が、誤差についてあたかも正規化されたようになっており、運動方程式モデルが仮に誤差が入り込みやすくても計算を行うことができ、モデルに求められる要求を小さくできるという効果をもたらすことができる。

誤差が最小となるような形状では、（２８）式の値が最小値になる。そこで、（２８）式の変分をとった、(２９)式の変分式を利用する。

Σ_ｔ及びΣ_ｘをそれぞれ、時間ｔ及び空間ｘに対する積分範囲とする。また、Γ_ｔ及びΓ_ｘをそれぞれ、時間ｔ及び空間ｘに対する境界値とする。上述した変分法の場合と同様に、部分積分を行って、δ_ξ、δ_ｙ、δ_λにかかる項に整理し、δ_ξ、δ_ｙ、δ_λに対するΓ_ｔ及びΓ_ｘでの境界条件が０とすると、観測値を反映した上で最適な形状推定結果ｙの値が得られる。

ここで、各種の部分積分δＪ_１〜δＪ_１０をそれぞれ、（３０）式〜（３９）式のように表すこととする。

部分積分δＪ_２については、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４０）式が得られる。

部分積分δＪ_３についても同様に、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４１）式が得られる。

部分積分δＪ_４についても同様に、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４２）式が得られる。

部分積分δＪ_５についても同様に、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４３）式が得られる。

部分積分δＪ_６についても同様に、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４４）式が得られる。

部分積分δＪ_７についても同様に、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４５）式が得られる。

部分積分δＪ_８についても同様に、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４６）式が得られる。

部分積分δＪ_１０についても同様に、δｙに対する境界条件が０になることを考慮して、部分積分を行うと、（４７）式が得られる。

Ｊが最小値になるためには、δＪ＝０、すなわち、δ_ξ、δ_ｙ、δλ_１及びδλ_２にかかる各項が０となれば良いので、(４８)式〜（５１）式の方程式を満たす。なお、（４９）式は、長い式であるので６分割して、Ａ＋Ｂ＋Ｃ＋Ｄ＋Ｅ＋Ｆ＋Ｇ＝０と置換え、Ａ、Ｂ、Ｃ、Ｄ、Ｅ、Ｆ、Ｇの部分ごとに書き出している。

これらの式を、変分法（アジョイント法）と呼ばれる方式で解くには、（４８）式〜（５１）式の方程式を差分化し、（５０）式に関しては、適当な初期値を用いて時間前進方向にΣ区間で解き、ｙについては、その（５０）式の計算結果を用い、（４８）式、（４９）式、（５０）式をλに対して時間方向の逆方向に解けば良い。得られた結果を用いて（５０）式を解くことを繰り返し、解が収束したところで推定値を得る。

（Ａ−３）実施形態の効果
上記実施形態によれば、運動方程式の誤差が最小になるが、運動方程式を厳密には満たさなくとも良くなるため、運動方程式（モデル）が完全ではない場合でも適用することができ、そのような場合にも、線状被曳航体形状の良好な推定値を得ることができる。

（Ｂ）他の実施形態
上記実施形態では、中性浮力の柔軟な被曳航体１０の変位を表す式としてＰａｉｄｏｕｕｓｓｉｓ方程式（（１）式、（２）式参照）を適用している場合を説明したが、被曳航体１０の運動を表すことができる全ての方程式に対して、本発明を適用することができる。例えば、中性浮力・重力、水平方向の運動・鉛直方向の運動など、形状を表現できる方程式を用いれば、この方法を用いることにより、耐性のある柔軟な被曳航体１０の形状推定が可能である。

また、上記実施形態では、主に、方位の計測値に基づいて推定する方法を説明したが、被曳航体の形状を観測する装置であれば、利用する観測値の種類は問われないものである。例えば、深度計、ジャイロ及び張力計の計測値などを利用して形状を推定するようにしても良い。

さらに、上記実施形態では、線状の被曳航体に受波器アレイが設けられているものを示したが、線状の被曳航体の用途などは、これに限定されるものではない。線状の被曳航体は、水中を曳航されるものに限定されず、水面上を曳航されるものであっても良い。

実施形態の線状被曳航体形状推定装置の機能的構成を示すブロック図である。 実施形態に係る被曳航体形状推定・受波信号修正装置の構成を示すブロック図である。 受波器アレイを装備した線状被曳航体の概略構成を示す説明図である。 従来の被曳航体形状推定・受波信号修正装置の構成を示すブロック図である。 従来のカルマンフィルタを用いたアレイ形状推定方法を示すブロック図である。 アレイ形状推定方法が適用する座標系の説明図である。 柔軟な被曳航体の微小長さ当りにかかる力を示す説明図である。 変分法を用いたアレイ形状推定方法を示すブロック図である。

符号の説明

２６Ａ…線状被曳航体形状推定装置（アレイ形状推定装置）、６１…センサ情報（方位・深度の時系列）、６２…磁場環境情報、６３…観測値（方位情報と深度情報の時系列）６４…曳航体船舶運動諸元の時系列、６５…海洋環境情報（潮流の流向・流速）、６６Ａ…弱拘束ラグランジアン（Ｌａｇｒａｎｇｉａｎ）関数計算部、６７…運動方程式時間前進計算部、６８…初期時刻形状の初期値、６９…アジョイント方程式時間後退計算部、７０…収束判定部、７１…初期時刻形状更新部、７２…形状推定結果。

Claims (6)

1. 曳航体によって水中若しくは水面上を曳航される柔軟性を有する線状の被曳航体の形状を推定する線状被曳航体形状推定方法において、
   観測値取込手段が、上記被曳航体の形状推定に利用する観測値の取り込みに、一定時間の観測値と形状推定結果の誤差を最小にする、弱拘束ラグランジアンを用いた方法を適用していることを特徴とする線状被曳航体形状推定方法。
2. 上記弱拘束ラグランジアンは、運動方程式の（左辺−右辺）の２乗を、運動方程式の誤差分散で割った関数を用いていることを特徴とする請求項１に記載の線状被曳航体形状推定方法。
3. 上記被曳航体が受波器アレイを有するものであることを特徴とする請求項１又は２に記載の線状被曳航体形状推定方法。
4. 曳航体によって水中若しくは水面上を曳航される柔軟性を有する線状の被曳航体の形状を推定する線状被曳航体形状推定装置において、
   上記被曳航体の形状推定に利用する観測値の取り込みに、一定時間の観測値と形状推定結果の誤差を最小にする、弱拘束ラグランジアンを用いた方法を適用している観測値取込手段を有することを特徴とする線状被曳航体形状推定装置。
5. 上記弱拘束ラグランジアンは、運動方程式の（左辺−右辺）の２乗を、運動方程式の誤差分散で割った関数を用いていることを特徴とする請求項５に記載の線状被曳航体形状推定装置。
6. 上記被曳航体が受波器アレイを有するものであることを特徴とする請求項４又は５に記載の線状被曳航体形状推定装置。

Images