JP4143481B2

JP4143481B2 - 移動発音体の運動分析方法及び装置

Info

Publication number: JP4143481B2
Application number: JP2003179880A
Authority: JP
Inventors: 裕輔中島; 義人井出; 一水谷口; 未央郡山
Original assignee: Hitachi Ltd
Current assignee: Hitachi Ltd
Priority date: 2003-06-24
Filing date: 2003-06-24
Publication date: 2008-09-03
Anticipated expiration: 2023-06-24
Also published as: JP2005017031A

Description

【０００１】
【発明の属する技術分野】
本発明は、移動発音体の運動分析方法に係り、特に、船舶等の何らかの音を発しながら移動する運動体からの音波を聴音器を用いて検出し、その運動パラメータ及び固有周波数の分析を行う移動発音体の運動分析方法に関する。
【０００２】
【従来の技術】
船舶等の移動発音体の運動分析方法に関する従来技術として、例えば、特許文献１等に記載された技術が知られている。この従来技術は、実際に聴音器で観測された音波の方位と、仮定した運動パターンで移動発音体が移動した場合に聴音器で観測されると思われる音波の方位との差から運動パターンの近似解を求めることにより、移動発音体の運動パラメータを解析するというものである。
【０００３】
【特許文献１】
特開平４−３７０７７９号公報
【０００４】
【発明が解決しようとする課題】
従来技術による移動発音体の運動分析方法は、周波数情報が使用されておらず、移動発音体の固有周波数を求めることができず、また、方位情報だけを使用しているため、最適な運動パラメータを一意的に求めることができない場合があるという問題点を有している。また、従来技術による運動分析方法は、運動パラメータをオペレータが推定しているが、推定すべきパラメータが、距離と速力、または、距離と針路の２つであり、これらをオペレータが推定し、パラメータを少しづつ変更しながら、最適解を得るには多くの時間を要するという問題点を有している。さらに、前述の従来技術は、各パラメータの全てのパターンに対する評価を行うことがオペレータでは実質的に不可能なため、方位残差量が十分小であるとオペレータが主観的に判断したパラメータの組み合わせを最適解としていたが、方位残差量がより小なパラメータの組み合わせが存在する可能性が残されており、求めた最適解が真の最適解とは言えないという問題点を有している。
【０００５】
本発明の目的は、前述した従来技術の問題点を解決し、移動発音体の運動パラメータを自動的に変更し、各パラメータの取り得る全てのパターンにおける方位残差量及び周波数残差量を数値的に比較することにより、運動パラメータ及び固有周波数の最適解を高速かつ正確に求めることができるようにした移動発音体の運動分析方法及び運動分析装置を提供することにある。
【０００６】
【課題を解決するための手段】
本発明によれば前記目的は、聴音器を用いて離散的に観測される、移動発音体の観測方位情報と観測周波数情報とから、前記移動発音体の運動パラメータと固有周波数とを分析する移動発音体の運動分析方法において、前記観測された移動発音体の前記観測方位情報と前記観測周波数情報とを入力された時間間隔で平滑化するステップと、前記平滑後の方位情報に基づく分析の基準データとなる方位情報の指定を受け付けるステップと、入力された仮定速力を入力された仮定速力の変更幅で順次変更するステップと、前記順次変更された仮定速力と基準データとから距離の変更範囲を計算するステップと、前回の仮定速力の分析結果から最適距離を予測し、距離変更範囲を修正するステップと、入力された仮定距離を入力された仮定距離の変更幅で順次変更するステップと、前記変更された仮定速力、前記変更された仮定距離、前記基準データから順次仮定運動パターンを計算するステップと、算出した前記仮定運動パターンから前記移動発音体の所定の時刻の方位情報を順次算出し、算出されたそれぞれの方位情報と同時刻の平滑後の前記方位情報との方位差分及び方位差分の二乗和である方位残差量を求めるステップと、前記各仮定速力における周波数残差量を最小とする最適固有周波数を計算するステップと、前記計算された固有周波数及び前記仮定運動パターンから前記移動発音体の所定の時刻の周波数情報を順次算出し、算出されたそれぞれの周波数情報と同時刻の平滑後の前記周波数情報との周波数差分及び周波数差分の二乗和である周波数残差量を求めるステップと、前記方位残差量と前記周波数残差量との和から最適運動パラメータを求めるステップとを有することにより達成される。
【０００８】
【発明の実施の形態】
以下、本発明による移動発音体の運動分析方法及び運動分析支援装置の実施形態を図面により詳細に説明する。
【０００９】
図１は本発明の一実施形態による移動発音体の運動分析支援装置の構成例を示すブロック図である。ここで説明する装置において、運動分析及び固有周波数分析の対象となる移動発音体を移動発音体６００とする。
【００１０】
本発明の実施形態による運動分析支援装置は、仮定運動パターン及び仮定固有周波数を評価し最適解を求める最適解計算装置１００と、求めた最適解及び方位と周波数との履歴情報を表示する表示装置２００と、オペレータからの入力を受け付けるオペレータ入力装置３００とを備えて構成される。そして、図１に示す実施形態は、前述の運動分析支援装置と、移動発音体６００の方位情報と周波数情報とを音圧として取り込む聴音器５００と、聴音器５００で取り込んだ、方位及び周波数に関する音圧を電気信号に変換する信号処理装置４００とを備えている。聴音器５００で収集され、信号処理装置４００を介して受け取った移動発音体６００に関する方位データ及び周波数データは、最適解計算装置１００内のデータベースに格納される。
【００１１】
最適解計算装置１００は、装置全体の制御を行い、本発明の実施形態における各種処理を実現するためのプログラムを実行する中央処理装置と、各種データ及びプログラムを格納する記憶部と、入出力のインターフェイスとを備えた一般的な構成を有する情報処理装置上に構築することが可能である。
【００１２】
表示装置２００は、分析結果表示部２１０と履歴情報表示部２２０とを備えたＣＲＴ、液晶ディスプレイ等の表示装置であり、後述する最適解計算装置１００の表示制御部１９５の制御に基づき、各表示部に最適解計算装置１００での処理の結果を表示する。また、履歴情報表示部２２０は、方位の履歴情報を表示する方位履歴情報表示部２２０ａと周波数の履歴情報を表示する周波数履歴情報表示部２２０ｂとを備えている。
【００１３】
オペレータ入力装置３００は、数値の入力を受け付けるテンキー３１０と、分析にあたって必要な各種設定を受け付けるスイッチ３２０と、表示装置２００に表示されたものに対する選択を受け付けるマウス等のポインティングディバイス３３０とを備える。ポインティングディバイス３３０は、ボールタブ等であってもよい。
【００１４】
最適解計算装置１００は、聴音器５００で収集され、信号処理装置４００を介して入手した観測データを、オペレータ入力装置３００で受け付けたオペレータによる設定値を用いて処理を行い、処理結果を表示装置２００に表示させるものである。
【００１５】
最適解計算装置１００の記憶部には、観測方位データベース７１０と、観測周波数データベース７２０と、分析用方位データベース７３０と、分析用周波数データベース７４０と、分析結果データベース７５０とが格納される。
【００１６】
観測用方位データベース７１０及び観測周波数データベース７２０は、聴音器５００より得られ、信号処理装置４００で信号化された、移動発音体６００の方位に関する離散時間的な観測データ及び周波数に関する離散時間的な観測データが格納されるデータベースである。
【００１７】
具体的には、観測方位データベース７１０には、観測時刻ｔ_i 毎の移動発音体６００の方位Ｂy_i及び聴音器５００（聴音器５００が設置されている観測船）の位置、速力、進行方向が格納され、観測周波数データベース７２０には、観測時刻ｔ_i 毎の移動発音体６００からの周波数及び聴音器（聴音器５００が設置されている観測船）の位置、速力、進行方向が格納される。
【００１８】
また、分析用方位データベース７３０及び分析用周波数データベース７４０は、それぞれ、後述する平滑方位計算処理部１１０ａ及び平滑周波数計算処理部１１０ｂにおける処理結果が格納されるデータベースである。分析結果データベース７５０は、各的速における最適距離、及びそれに対応する的針、固有周波数、方位残差量、周波数残差量が格納される。
【００１９】
最適解計算装置１００は、これらのデータベースに格納されているデータに基づき、移動発音体の運動パラメータ及び固有周波数の分析を行う以下に説明する機能をプログラムの実行によって実現する。プログラムの実行によって実現される機能は、平滑処理部１１０と、基準データ指定部１２０と、速力変更幅指定部１２５と、仮定速力変更部１３０と、距離変更幅指定部１３５と、距離範囲計算部１４０と、仮定距離変更部１４５と、仮定運動パターン計算部１５０と、方位残差分析部１６０と、仮定固有周波数計算部１７０と、周波数残差分析部１８０と、最適解判定部１９０と、表示制御部１９５とである。
【００２０】
以下、前述の各機能構成部のそれそれが実行する処理について説明する。
【００２１】
１．平滑処理部１１０
平滑処理部１１０は、ばらつきのある観測データから誤差を除去するため、観測データに平滑処理を施すもので、観測方位データベース７１０及び観測周波数データベース７２０に２Ｔ_h 秒内に蓄積されたデータから、オペレータが指定した一次平滑時隔Ｔ_h 秒毎の分析用データを算出する。算出された分析用データは、分析用方位データベース７３０及び分析用周波数データベース７４０に格納される。
【００２２】
一次平滑時隔Ｔ_h は、オペレータがテンキー３１０を用いて入力した指定データを受け付けることによって取得する。方位情報に対する平滑処理は平滑方位計算処理部１１０ａで、周波数情報に対する平滑処理は平滑周波数計算処理部１１０ｂで、それぞれ独立して実行される。以下に各部の処理を説明する。
【００２３】
１−ａ．平滑方位計算処理部１１０ａ
ここでの処理は、観測方位データベース７１０に２Ｔ_h 秒内に蓄積されたデータから以下に説明する線形回帰により平滑方位、平滑方位変化率及び平滑方位分散を求めるものである。
【００２４】
図２は平滑処理部で実施する一次平滑の原理を説明する図である。図２において、ｙ軸方向は時刻、ｘ軸方向は方位を示す。Ｔ_h は、前述のように、オペレータによって指定された一次平滑時隔であり、ｙ軸のｔ−３Ｔ_h 、ｔ−２Ｔ_h 、ｔ−Ｔ_h 、ｔは一次平滑時刻を示す。また、ｄ_jk（ｊ：平滑区間番号、ｋ：平滑区間内データ番号）は観測データＢy_iをプロットしたものである。
【００２５】
平滑は、Ｔ_h 秒間を１平滑区間とし、連続した２つの平滑区間のデータを用いて、１平滑区間ずつ行う。例えば、図２において、ｔ−３Ｔ_h からｔ−２Ｔ_h の区間を平滑区間１、ｔ−２Ｔ_h からｔ−Ｔ_h からの区間を平滑区間２、ｔ−２Ｔ_h からｔの区間を平滑区間３とすると、平滑区間１の平滑処理は、平滑区間１及び平滑区間２のデータを用いて行われ、平滑区間２の平滑処理は、平滑区間２及び平滑区間３のデータを用いて行われる。
【００２６】
ここで、平滑区間１における平滑処理について説明する。はじめに、平滑区間１及び平滑区間２の２つの平滑区間の観測データ（２Ｔ_h 秒内に蓄積された観測データ）を用いて平滑直線１₁₂を求める。この直線上における平滑区間１内の最新観測データｄ₁₆の観測時刻に対する平滑直線１₁₂上の方位を平滑方位Ｄk₁とする。同様に、平滑区間２及び平滑区間３の観測データを用いて、平滑直線１₂₃を求め、平滑方位Ｄk₂を計算する。
【００２７】
以下に平滑方位を計算する式を示す。
【００２８】
まず、一次平滑係数を以下に数１として示す式（１）〜式（５）により算出する。
【００２９】
【数１】

【００３０】
式（１）〜式（５）において、ｔ_i は観測時刻、ｂy_iは時刻ｔ_i における聴音器観測方位、ｔb は平滑基準時刻、Ｂb は平滑基準方位、Ａb は方位差成分の総和、Ａt は時間差成分の総和、Ａbtは方位差成分と時間差成分との積の総和、Ａt²は時間差成分の２乗総和、Ａb²は方位差成分の２乗総和、Ｎは被平滑区間及び次の平滑区間内の全観測データ数である。但し、観測方位は真北を０°とする絶対方位とする。また、ｂy_iは、観測方位データベース７１０内に格納されているデータであり、ｔb 及びＢb は、観測方位データベース７１０内の最古のデータの時刻及び方位である。
【００３１】
前述の式（１）〜式（５）により算出された一次平滑係数を用い、平滑方位変化率ＢＤＢ、平滑基準時刻ｔb における平滑方位ＳＢＹ、平滑方位Ｂy 及び平滑方位分散σＢy を、数２として示す式（６）〜式（９）により算出する。
【００３２】
【数２】

【００３３】
前述により求められた、平滑方位Ｂy 及び平滑方位分散σＢy は、平滑区間内最新観測時刻と共に、分析用方位データベース７３０に格納される。このとき、観測方位データベース７１０から取得した、その平滑区間内最新観測時刻における聴音器情報（位置、速力、進行方向）に関するデータも、同時に分析用方位データベース７３０に格納される。以後、ここで算出された平滑方位Ｂy を平滑方位と呼び、それに対応する平滑区間内最新観測時刻を方位平滑時刻と呼ぶ。
【００３４】
１−ｂ．平滑周波数計算処理部１１０ｂ
次に、観測周波数データベースに格納されているデータから、分析用周波数データベースに格納するデータを算出する平滑周波数計算処理部１１０ｂについて説明する。平滑周波数計算処理部１１０ｂは、観測周波数データベース７２０に２Ｔ_h 秒内に蓄積されたデータから、前述の線形回帰により平滑周波数、平滑周波数変化率及び平滑周波数分散を求めるものである。平滑処理の原理は、図２を参照して説明した方位の平滑の場合と同様である。
【００３５】
以下に平滑周波数を計算する式を示す。
【００３６】
まず、一次平滑係数を、数３として示す式（１０）〜式（１４）により算出する。
【００３７】
【数３】

【００３８】
式（１０）〜式（１４）において、ｔ_j は観測時刻、ｆ_j は時刻ｔ_j における聴音器観測周波数、ｔb は平滑基準時刻、Ｆb は平滑基準周波数、Ａf は周波数差成分の総和、Ａt は時間差成分の総和、Ａftは周波数差成分と時間差成分との積の総和、Ａt²は時間差成分の２乗総和、Ａf²は周波数差成分の２乗総和、Ｎは被平滑区間及び次の平滑区間内の全観測データ数である。また、ｆ_j は観測周波数データベース９２０に格納されているデータであり、ｔb 及びＦb は、観測周波数データベース９２０の最古のデータの時刻及び周波数である。
【００３９】
前述の式（１０）〜式（１４）で算出された一次平滑係数を用い、平滑周波数変化率ＳＤＦ、平滑基準時刻における平滑周波数Ｆ及び平滑周波数分散σＦを、以下に数４として示す式（１５）〜式（１８）により算出する。
【００４０】
【数４】

【００４１】
前述により求められた平滑周波数Ｆ及び平滑周波数分散σＦは、平滑区間内最新観測時刻と共に、分析用周波数データベース７４０に格納される。このとき、観測周波数データベースから取得した、その平滑区間内最新観測時刻における聴音器情報（位置、速力、進行方向）に関するデータも、同時に分析用周波数データベース７４０に格納される。以後、ここで算出された平滑周波数Ｆを、平滑周波数と呼び、それに対応する平滑区間内最新観測時刻を周波数平滑時刻と呼ぶ。
【００４２】
２．基準データ指定部１２０
基準データ指定部１２０は、オペレータがポインティングディバイス３３０を用いて選択した基準データ１及び基準データ２を受け付け、これらのデータを後述する距離範囲計算部１４０、仮定運動パターン計算部１５０、履歴情報表示部２２０に受け渡す。前述の基準データ１、２は、先に述べた方位履歴表示部２２０ａに表示されているデータの中からオペレータが分析の基準となるデータとして選択するものである。但し、基準データ１は、基準データ２より過去のデータとする。また、初期値として、最古の平滑方位データを基準データ１、最新の平滑方位データを基準データ２と設定する。後述するように、方位履歴表示部２２０ａに表示されているデータのうち、基準データとして指定を受け付けたデータには、破線円を重畳表示させる。
【００４３】
３．速力変更幅指定部１２５
速力変更幅指定部１２５は、オペレータがテンキー３１０を用いて入力した速力変更幅ΔＭt 及び移動発音体の予想最高速力Ｍt_maxを受け付け、これらを後述する仮定速力変更部１３０に受け渡す。
【００４４】
４．仮定速力変更部１３０
仮定速力変更部１３０は、運動パラメータの１つである速力を速力変更幅指定部１２５より受け付けた速力変更幅ΔＭt の間隔で増加させ、それを仮定速力Ｍt として、距離範囲計算部１４０及び仮定運動パターン計算部１５０に受け渡すものである。仮定速力Ｍt の変更範囲は、０ktから移動発音体の予想最高速力Ｍt_maxまでとする。
【００４５】
５．距離変更幅指定部１３５
距離変更幅指定部１３５は、オペレータがテンキー３１０を用いて入力した距離変更幅ΔＲを受け付け、これを後述する仮定距離変更部１４５に受け渡す。
【００４６】
６．距離範囲計算部１４０
距離範囲計算部１４０は、仮定速力Ｍt と基準データ１、２における平滑方位と方位平滑時刻とから、仮定速力Ｍt に対する距離の変更可能範囲を幾何学的条件により計算する。計算により求めた距離の変更可能範囲は、距離最小値Ｒ_min 及び距離最大値Ｒ_max として仮定距離変更部１４５に受け渡される。幾何学的条件による計算の方法は、仮定運動パターン計算部１５０の説明の部分で行うこととする。
【００４７】
７．仮定距離変更部１４５
仮定距離変更部１４５は、運動パラメータの１つである距離を距離変更幅指定部１３５より受け付けた速力変更幅ΔＲの間隔で増加させ、それを仮定距離Ｒとして仮定運動パターン計算部１５０に受け渡すものである。仮定距離Ｒの変更範囲は距離最小値Ｒ_min から距離最大値Ｒ_max までとする。
【００４８】
８．仮定運動パターン計算部１５０
仮定運動パターン計算部１５０は、運動パラメータと方位に関する基準データとを用い、仮定運動パターンを計算するものである。
【００４９】
仮定運動パターンの計算に用いる運動パラメータは、仮定速力変更部１３０から受け渡された仮定速力Ｍt 、仮定距離変更部１４５から受け渡された仮定距離Ｒであり、基準データは、基準データ指定部１２０においてオペレータが指定したものを受け付けた基準データ１及び基準データ２である。これらを用いて、仮定運動パターン計算部１５０は、基準データ１における平滑方位方向に移動発音体６００があると仮定し、移動発音体の針路Ｃt を算出する。
【００５０】
図３は仮定運動パターンを求める方法の原理を説明する図である。図３において、基準データ１における聴音器５００の位置を原点とし、基準データ２における聴音器５００の位置を（Ｘo ，Ｙo ）とする。基準データ１における平滑方位をｙ軸方向とし、原点からｙ軸上の距離Ｒの点に移動発音体６００を仮定する。その点を中心とし半径をＭt・Ｔとする円をＣとする。ここで、Ｒは仮定距離、Ｍt は仮定速力であり、Ｔは基準データ１から基準データ２までの時間間隔である。（Ｘo ，Ｙo ）点からｙ軸とＥb なる角度をなす直線をＬとし、円Ｃとの２つの交点をＬa （Ｘa ，Ｙa ）及びＬb （Ｘb ，Ｙb ）とする。ここで、図３に示すＥb は、基準データ１と基準データ２との平滑方位の差であり、直線Ｌの方向は、基準データ２の平滑方位の方向である。仮定運動パターン計算のアルゴリズムは、（Ｘo ，Ｙo ），Ｒ，Ｍt ，Ｔ及びＥb が与えられたとき、これらから円Ｃと直線Ｌとの交点Ｌa （Ｘa ，Ｙa ）またはＬb （Ｘb ，Ｙb ）を求め、さらに、それらの交点Ｌa またはＬb と円Ｃの中心とを結ぶ直線と、基準データ１における平滑方位とがなす角度Ｕa またはＵb を求めるものである。
【００５１】
前述において、基準データ１における聴音器５００の位置を原点とした基準データ２における聴音器の相対位置（Ｘo ，Ｙo ）は、分析用方位データベース７３０から取得した、基準データ１における平滑方位とそれぞれの基準データの聴音器位置から算出する。また、基準データ１から基準データ２までの時間間隔Ｔは、それぞれの基準データの方位平滑時刻から算出する。仮定距離Ｒ及び仮定速力Ｍt は、それぞれ基準データ１の時刻における移動発音体６００と聴音器５００とからの距離及び移動発音体６００の速力の仮定値であり、前述したように仮定速力変更部１３０及び仮定距離変更部１４５より与えられたものを受け付ける。これらを用いて求めた角度Ｕa またはＵb と基準データ１における平滑方位とから、以下のようにして移動発音体の針路を求めることができる。
【００５２】
図３において、円の方程式（１９）と直線の方程式（２０）とは、数５として示す通りであり、式（２０）を式（１９）に代入してについて整理すると、数５として示す式（２１）に示すようになる。
【００５３】
【数５】

【００５４】
前述の式（２１）をｘについて解くと、円Ｃと直線Ｌとの２つの交点（Ｘa ，Ｙa ）及び（Ｘb ，Ｙb ）は、数６として示す式（２２）〜式（２５）のように表される。
【００５５】
【数６】

【００５６】
また、基準データ２における聴音器とそれぞれの解との距離をＲa 及びＲb とすると、Ｒa 及びＲb は、数７として示す式（２６）、式（２７）のように表される。
【００５７】
【数７】

【００５８】
ここで、Ｒa ＞Ｒb ならば（Ｘa ，Ｙa ）が遠パターンの解、（Ｘb ，Ｙb ）が近パターンの解となる。逆に、Ｒa ＜Ｒb ならば、（Ｘa ，Ｙa ）が近パターンの解、（Ｘb ，Ｙb ）が遠パターンの解となる。（Ｘa ，Ｙa ）及び（Ｘb ，Ｙb ）をそれぞれ遠パターンまたは近パターンの場合の基準データ２における移動発音体６００の位置とする。円Ｃの中心と（Ｘa ，Ｙa ）及び（Ｘb ，Ｙb ）とを結ぶ直線と基準データ１における平滑方位とのなす角度をそれぞれＵa 及びＵb とすると、Ｕa 及びＵb は、数８に示す式（２８）、式（２９）のように表すことができる。また、基準データ１における平滑方位Ｂysを用い、移動発音体の針路Ｃta及びＣtbを求めると、移動発音体の針路Ｃta及びＣtbは、数８に示す式（３０）、式（３１）のように求められる。
【００５９】
【数８】

【００６０】
前述の２つの仮定運動パラメータの組合せ（Ｒ，Ｍt ，Ｃta）及び（Ｒ，Ｍt ，Ｃtb）に関して、それぞれ方位残差分析処理を行う。方位残差分析処理については、次の方位残差分析部１６０の説明の部分で説明する。
【００６１】
前述をふまえた上で、以下に距離範囲計算の説明を行う。前述の説明は、直線Ｌと円Ｃとが交点を持つ場合であったが、受け渡された仮定速力Ｍt 及び仮定距離Ｒによっては直線Ｌと円Ｃとが交わらない場合もある。図３において、Ｍt・Ｔ及びＥb を一定に保ったまま、距離Ｒを増加させていくと、いずれ直線Ｌと円Ｃとは交わらなくなる。
【００６２】
図４は速力における距離変更範囲の最大値について説明する図であり、前述した直線Ｌと円Ｃとが交わらなくなる臨界点、すなわち、図４に示すような直線Ｌと円Ｃとが接する点では式（２１）の解が重根となる。重根の場合、式（２１）に関する判別式は０となるので、判別式及びその判別式をＲに関して解いたものは、数９として示す式（３２）、式（３３）のようになる。
【００６３】
【数９】

【００６４】
式（３３）に示すＲが仮定速力Ｍt に対する距離の変更可能範囲の最大値、すなわち、距離最大値Ｒ_max である。但し、算出したＲ_max が聴音器の最大感知距離を越える場合には、聴音器の最大感知距離をＲ_max とする。
【００６５】
図５は速力における距離変更範囲の最小値について説明する図である。この図５に示すように、基準データ１と基準データ２との平滑方位線が交差する場合、直線Ｌと円Ｃとが接する点での判別式をＲに関して解いたものは、数１０として示す式（３４）のようになる
【数１０】

【００６６】
式（３４）に示すＲが仮定速力Ｍt に対する距離の変更可能範囲の最小値、すなわち、距離最小値Ｒ_min である。ただし、Ｒ_min が負となる場合には０をＲ_min とする。
【００６７】
９．方位残差分析部１６０
方位残差分析部１６０は、各方位平滑時刻において、前述の仮定運動パターン計算部１５０で求めた運動パターンで移動発音体６００が移動した場合に、聴音器５００で観測されると思われる方位（分析方位）と、実際に聴音器５００で観測された方位（平滑方位）との差を方位残差として求める。
【００６８】
基準データ１における聴音器５００に対する移動発音体６００の相対位置を（ｘ_s ，ｙ_s ）とし、基準データ１と各方位平滑時刻との時間間隔をＴ_i 、仮定運動パターン計算部１５０で求められた針路をＣ_t とすると、各方位平滑時刻における移動発音体６００の基準データ１に対する相対位置Ｐ_i （ｘ_i ，ｙ_i ）は、数１１として示す式（３５）、式（３６）のようになる。
【００６９】
【数１１】

【００７０】
図６は方位残差算出の原理を説明する図である。各観測時刻における聴音器５００からみた移動発音体６００の分析方位は、数１２に示す式（３７）の通りとなる。
【００７１】
【数１２】

【００７２】
式（３７）においてＸo_i及びＹo_iは、各観測時刻における聴音器５００の相対位置、すなわち、基準データ１における聴音器５００の位置を基準とした場合の東西（経度）成分及び南北（緯度）成分である。
【００７３】
各方位平滑時刻における実際の平滑方位をＢy_iとすると、方位残差jＢ_iは、数１３に示す式（３８）のようになる。
【００７４】
【数１３】

【００７５】
前述のようにして求められた各方位平滑時刻における方位残差と各方位平滑時刻とにおける平滑方位分散σＢy_iより、方位残差量δb を求める。この方位残差量δb は、数１４に示す式（３９）のようになる。
【００７６】
【数１４】

【００７７】
この方位残差量δb は、仮定した運動パターンが移動発音体６００の実際の運動に近いものかどうかの評価指標とされる。すなわち、方位残差量δb が微小値であるほど、仮定した運動パターンは移動発音体６００の実際の運動に近いものであると評価することができる。
【００７８】
方位残差分析は、遠パターンと近パターンとのそれぞれの針路Ｃt に対して行い、方位残差量δb の小さい方が仮定速力Ｍt における仮定距離の解析値として選択される。選択した方位残差量δb と仮定速力Ｍt における前回仮定距離までの最小方位残差量δb^〜とを比較し、今回の方位残差量δb が最小方位残差量δb^〜より小であるなら、仮定距離Ｒ、針路Ｃt 、今回の方位残差量δb を仮定速力Ｍt における最適解として、分析結果データベースに格納する。すなわち、今回の方位残差量δb が仮定速力Ｍt における最小方位残差量δb^〜となり、仮定距離Ｒ及び針路Ｃt が、それぞれ最適距離Ｒ^〜及び最適針路Ｃt^〜となる。
【００７９】
なお、前述において、「^〜」の記述は、^〜の直前の文字の上に^〜が付与されたものであることを示しており、以後も同様である。また、「^ 」の記述についても同様である。
【００８０】
前述した処理をＲ_min からＲ_max からまで仮定距離ＲをΔＲだけ増加しながら繰り返し行うことにより、仮定速力Ｍt における最適解（最小方位残差量δb´、最適距離Ｒ´、最適針路Ｃt´）を求めることができる。
【００８１】
また、前述の処理を、０ktからＭt_maxまで仮定速力Ｍt をΔＭt だけ増加しながら繰り返し行うことにより、各仮定速力Ｍt_kにおける最適解（最小方位残差量δb^〜、最適距離Ｒ^〜 _k、最適針路Ｃt^〜 _k ）を求めることができる。ここで、（ｋ＝１，２，３，‥‥‥，ｍ）であり、ｍは、Ｍt_maxをΔＭt で割った商の整数部である。
【００８２】
１０．最適距離予測部１６５
移動発音体６００の仮定速力Ｍt_kが聴音器５００の移動速力より大きい場合、仮定速力Ｍt_k，Ｍt_k+1と最適距離Ｒ^〜 _k，Ｒ^〜 _k+1との間に、数１５に示す式（４０）のような比例関係が成り立つ場合がある。
【００８３】
【数１５】

【００８４】
そこで、最適距離予測部１６５は、速力変更幅ΔＭt が１kt以上で、異なる３つの仮定速力Ｍt_kとＭt_k+1、Ｍt_k+1とＭt_k+2の間で式（４０）に示す比例関係が±５％以内で成り立つ場合に、Ｍt_k+3にも比例関係が成り立つと仮定し、Ｍt_k+3における予測最適距離Ｒ^_k+3を、数１６に示す式（４１）により求める。
【００８５】
【数１６】

【００８６】
最適距離Ｒ^〜 _k+3は、予測最適距離Ｒ^_k+3の付近に存在すると考えられるので、距離最大値Ｒ_maxk+3、距離最小値Ｒ_mink+3を、数１７に示す式（４２）、式（４３）により求める。
【００８７】
【数１７】

【００８８】
式（４２）、式（４３）におけるΔＲ^_k+3は、式（４４）に示すような予測変更幅であり、この予測変更幅ΔＲ^_k+3が距離変更幅ΔＲの５倍未満の場合、予測変更幅ΔＲ^_k+3を距離変更幅ΔＲの５倍とする。このようにして計算された距離最大値Ｒ_maxk+3及び距離最小値Ｒ_mink+3は、前述した仮定距離変更部１４５に受け渡される。
【００８９】
これにより距離Ｒの変更範囲が狭められることになり、よって各仮定速力において距離Ｒに対する反復計算回数を減少させることができ、最適解を得るまでの処理時間を短縮することができる。
【００９０】
１１．仮定固有周波数計算部１７０
仮定速力Ｍt_kにおける最適距離Ｒ^〜 _k、最適針路Ｃt^〜 _k が与えられた場合、周波数残差量δf を最小とする最適固有周波数Ｆ^〜 _kを計算する。最適固有周波数の計算の方法については、周波数残差分析部１８０で説明する。
【００９１】
１２．周波数残差分析部１８０
周波数残差分析部１８０は、各周波数平滑時刻において、前述した仮定固有周波数計算部１７０で求めた最適固有周波数を発っしつつ、仮定運動パターン計算部１５０で求めた運動パターンで移動発音体６００が移動した場合に、聴音器５００で観測されると思われる周波数（分析周波数）と、実際に聴音器５００で観測された周波数（平滑周波数）との差を周波数残差として求める。
【００９２】
図７は仮定運動パターンにおける聴音器と移動発音体との相対運動を説明する図である。
【００９３】
図７に示すように、移動発音体６００と聴音器５００とは相対運動をしているので、観測される周波数は、移動発音体の固有周波数からドップラー効果の分だけずれて観測される。
【００９４】
仮定固有周波数計算部１７０から与えられた最適固有周波数をＦ^〜とすると各周波数平滑時刻における分析周波数Ｆ´_jは、ドップラー効果の公式による数１８に示す式（４６）により求めることができる。
【００９５】
【数１８】

【００９６】
数１８に示す式（４６）において、Ｖは水中音速、Ｍo_jは聴音器５００の速力、Ｃo_jは聴音器５００の針路方向、Ｂｙ´_jは分析方位、Ｍt は移動発音体６００の仮定速力、Ｃt^〜は仮定速力Ｍt における移動発音体６００の最適針路である。
【００９７】
各周波数平滑時刻における実際の平滑周波数をＦ_j とすると、分析周波数Ｆ´_jを用い、周波数残差jＦ_jは、数１９に示す式（４７）により求められる。
【００９８】
【数１９】

【００９９】
前述のようにして求められた各周波数平滑時刻における周波数残差jＦ_jと各周波数平滑時刻における平滑周波数分散σＦ_j より、周波数残差量δf は、数２０に示す式（４８）により求められる。
【０１００】
【数２０】

【０１０１】
前述までで説明したようにして求められた各仮定的速Ｍt_k における周波数残差量δf_kは、最適固有周波数Ｆ^〜 _kと共に分析結果データベースに格納される。
【０１０２】
ここで、最適固有周波数Ｆ^〜の計算の方法を説明しておく。仮定速力Ｍt における最適距離Ｒ^〜、最適針路Ｃt^〜が与えられた場合、周波数残差量δf を最小とする仮定固有周波数Ｆ^〜は、前述の式（４８）より、数２１に示す式（４９）の方程式を満たすＦ^〜である。
【０１０３】
【数２１】

【０１０４】
この式（４９）を解くと、Ｆ^〜は、数２２に示す式（５０）となる。
【０１０５】
【数２２】

【０１０６】
このようにして求められた最適固有周波数Ｆ^〜を用いることにより、仮定速力Ｍt 、最適距離Ｒ^〜、最適針路Ｃt^〜における最小周波数残差量δf^〜を求めることができる。
【０１０７】
１３．最適解判定部１９０
最適解判定部１９０は、以下に説明するように演算を行う。まず、仮定速力Ｍt_jにおける最小方位残差量δb_jと最小周波数残差量δf^〜 _j の和を総残差量δs^〜 _j として、数２３に示す式（５１）により総残差量δs^〜 _j を求める。
【０１０８】
【数２３】

【０１０９】
求められた総残差量δs^〜 _j は、各仮定速力Ｍt_jにおける最適解の評価指標とされる。各仮定速力における総残差量を比較し、総残差量が最も小さい仮定速力における、最適距離、最適針路、最適固有周波数及びこの仮定速力を移動発音体の仮定運動パラメータの最適解とする。以後、この最適解を総最適解と呼び、総最適解における仮定速力を最適速力と呼ぶ。但し、周波数情報が聴音器で得られなかった場合には、各仮定速力における最小方位残差量を総残差量とし、総最適解を得るものとする。
【０１１０】
１４．表示制御部１９５
表示制御部１９５は、最適解計算装置１００の各部で行われた処理の結果を表示装置２００に表示させる。その際、以下に、１２−１として説明する分析結果表示制御、１２−２として説明する履歴情報表示制御の２種類の制御を行う。以下、それぞれの制御について説明する。
【０１１１】
１２−１．分析結果表示制御
前述の最適解判定部１９０で求めた総最適解の仮定速力Ｍt 、最適距離Ｒ^〜、最適針路Ｃt^〜、最適固有周波数Ｆ^〜、最小方位残差量δb^〜、最小周波数残差量δf^〜、総残差量δs^〜を数値情報として分析結果表示部２１０に表示する。また、総最適解の最適距離Ｒ^〜を用いて計算される基準データ１における移動発音体６００の相対位置を（ｘ_s ，ｙ_s ）とし、基準データ１と現在時刻の時間間隔をΔＴｎとして、総最適解の最適速力Ｍt^〜と最適針路Ｃt^〜とを用いて、現在時刻における移動発音体６００の相対位置Ｐ_n(ｘ_n ，ｙ_n ）は、数２４に示す式（５２）、式（５３）により求められる。
【０１１２】
【数２４】

【０１１３】
現在時刻における聴音器５００の相対位置をＯ_n(ｘo_n，ｙo_n）とし、Ｏ_n からみたＰ_n の方位を現在方位ＢＹ_n 、Ｏ_n とＰ_n との距離を現在距離ＲＲ_n とすると、現在方位ＢＹ_n 及び現在距離ＲＲ_n は、数２５に示す式（５４）、式（５５）により求められる。
【０１１４】
【数２５】

【０１１５】
前述のようにして求められた移動発音体６００の現在方位ＢＹ_n 及び現在距離ＲＲ_n も分析結果表示部２１０に併せて表示される。
【０１１６】
１２−２．履歴情報表示制御
履歴情報表示制御は、平滑処理部１１０により求められて、分析用方位データベース７３０および分析用周波数データベース７４０に格納されたデータを、表示装置２００の履歴情報表示部２２０に、時刻に対するプロットとして表示する制御を行う。
【０１１７】
１２−２−ａ．方位履歴情報表示処理２２０ａ
方位履歴情報表示処理２２０ａは、分析用方位データベース７３０に格納されている平滑方位Ｂy 及び平滑方位分散σＢy を、履歴情報表示部２２０の方位履歴情報表示部２２０ａに表示する処理を行う。
【０１１８】
図８は平滑方位Ｂy 及び聴音器針路方向を時間に対してプロットすることにより表示させた様子を示す図である。
【０１１９】
図８に示す表示例は、一方の軸を方位、もう一方の軸を時刻として各方位平滑時刻における平滑方位Ｂy をプロットして表示した場合の例である。各平滑方位Ｂy に対する平滑方位分散σＢy は、その平滑方位Ｂy を中心とし、幅２σＢy の誤差棒を付けて表示される。同時に、各方位平滑時刻における聴音器５００の針路方向も合わせて表示される。基準データ指定部１２０で受け付けた基準データ１及び基準データ２にあたる平滑方位Ｂy のプロットには破線円が重畳表示される。以下、このプロットを時刻方位履歴表示と呼ぶ。この時刻方位履歴表示の表示中心方位Ｂc 及び表示方位スケールＢs 、時間表示スケールＴs は、任意に設定可能である。この設定は、オペレータがオペレータ入力装置３００のテンキー３１０及びスイッチ３２０を用いて入力することにより行われる。なお、時間表示スケールＴs は、後に説明する時刻周波数履歴情報表示と連動して切り替えられる。
【０１２０】
１２−２−ｂ．周波数履歴情報表示処理２２０ｂ
周波数履歴情報表示処理２２０ｂは、分析用周波数データベース７４０に格納されている平滑周波数Ｆ及び平滑周波数分散σＦを、履歴情報表示部２２０の周波数履歴情報表示部２２０ｂに表示する処理を行う。
【０１２１】
図９は平滑周波数Ｆを時刻に対してプロットして表示させた様子を示す図である。
【０１２２】
図９に示す表示例は、一方の軸を周波数、もう一方の軸を時刻として各周波数平滑時刻における平滑周波数Ｆをプロットして表示した例である。各平滑周波数Ｆに対する平滑周波数分散σＦは、平滑周波数Ｆを中心とし、幅２σＦの誤差棒を付けて表示される。以下、このプロットを時刻周波数履歴表示と呼ぶ。この時刻周波数履歴表示の表示中心周波数Ｆc 及び表示周波数スケールＦs 、時間表示スケールＴs は任意に設定可能である。Ｆc 、Ｆs 、Ｔs の設定は、オペレータがオペレータ入力装置３００のテンキー３１０及びスイッチ３２０を用いて入力することにより行われる。また、総最適解が得られている場合には、最適固有周波数Ｆ^〜も併せて表示される。
【０１２３】
図１０は最適解計算装置１００の処理動作を説明するフローチャートであり、次に、これについて説明する。
【０１２４】
（１）まず、平滑処理部１１０は、オペレータによる平滑時隔の入力を受け付け、観測方位データベース７１０及び観測周波数データベース７２０に蓄積されているデータに対して、その平滑時隔を用いた平滑処理を行い、その処理結果を分析用方位データベース７３０及び分析用周波数データベース７４０に格納する（ステップＳ００１）。
【０１２５】
（２）次に、表示制御部１９５は、分析用方位データベース７３０及び分析用周波数データベース７４０に格納されたデータを時刻方位履歴表示及び時刻周波数履歴表示として、表示装置２００の履歴情報表示部２２０に表示する（ステップＳ００２）。
【０１２６】
（３）ここで、オペレータは、ステップＳ００２の処理で表示される時刻方位履歴表示を見て、基準データを選択することができる。基準データ指定部１２０は、オペレータ入力装置３００を介して、オペレータからの基準データに関する指定を受け付け、そのデータを表示制御部１９５、距離範囲計算部１４０、仮定運動パターン計算部１５０に受け渡す（ステップＳ００３）。
【０１２７】
（４）次に、オペレータは、仮定速力の速力変更幅と移動発音体の予想最高速力とを指定する。速力変更幅指定部１２５は、オペレータ入力装置３００を介して前述のオペレータの速力変更幅と移動発音体の予想最高速力とに関する指定を受け付け、それらの情報を仮定速力変更部１３０に受け渡す（ステップＳ００４）。
【０１２８】
（５）仮定速力の場合と同様に、オペレータは、仮定距離の距離変更幅を指定する。距離変更幅指定部１３５は、オペレータ入力装置３００を介して前述のオペレータの距離変更幅に関する指定を受け付け、その情報を仮定距離変更部１４５に受け渡す（ステップＳ００５）。
【０１２９】
（６）仮定速力変更部１３０は、仮定速力に速力変更幅指定部１２５より受け付けた速力変更幅を順次加算し、距離範囲計算部１４０及び仮定運動パターン計算部１５０に加算した結果の変更後の仮定速力順次を受け渡す。なお、この処理は、仮定速力０ktから始め、仮定速力が移動発音体の予想最高速力を上回るまで繰り返し、上回った時点で当処理を終了する（ステップＳ００６）。
【０１３０】
（７）距離範囲計算部１４０は、ステップＳ００６の処理で変更された仮定速力を順次受け取り、仮定速力が予想最高速力Ｍt_maxより小さいか否かを判定し、仮定速力が予想最高速力を越えた場合、後述するステップＳ０１４の処理に移行する（ステップＳ００６’）。
【０１３１】
（８）ステップＳ００６’の判定で、仮定速力が予想最高速力Ｍt_maxより小さかった場合、距離範囲計算部１４０は、基準データ指定部１２０から受け渡された基準データ情報と仮定速力変更部１３０から受け渡された仮定速力とから、その仮定速力に対する距離最小値及び距離最大値を計算し、計算結果を仮定距離変更部１４５に受け渡す（ステップＳ００７）。
【０１３２】
（９）ここで、最適距離予測部１６５は、分析結果データベース７５０を参照し、前回３回分の仮定速力と最適距離との間に比例関係が成り立っていた場合に、今回の仮定速力における最適距離を予測し、距離最小値及び距離最大値を予測した最適距離の近辺とし、これらを仮定距離変更部１４５に受け渡す（ステップＳ００８）。
【０１３３】
（10）仮定距離変更部１４５は、仮定距離に距離変更幅指定部１３５より受け付けた距離変更幅を順次加算し、加算結果を順次仮定運動パターン計算部１５０に受け渡す。なお、この処理は、距離最小値から始め、仮定距離が距離最大値を上回るまで繰り返し、上回った時点で当処理を終了する（ステップＳ００９）。
【０１３４】
（11）仮定運動パターン計算部１５０は、ステップＳ００9の処理で変更された仮定距離を順次受け取り、仮定距離が予想最高距離Ｒ_max より小さいか否かを判定し、仮定速力が予想最高速力を越えた場合、後述するステップＳ０１２の処理に移行する（ステップＳ００９’）。
【０１３５】
（12）ステップＳ００９’の判定で、仮定距離が予想最高距離Ｒ_max より小さかった場合、仮定運動パターン計算部１５０は、前述したステップＳ００３、Ｓ００６、Ｓ００９の処理で受け付けた、基準データ、仮定速力、仮定距離を用い、仮定運動パターンを計算し、その結果を方位残差分析部１６０に受け渡す（ステップＳ０１０）。
【０１３６】
（13）方位残差分析部１６０は、ステップＳ０１０の処理で計算された仮定運動パターンを基に、方位残差及び方位残差量を求める。この処理は、ステップＳ０００９からの処理を繰り返して行い、その仮定速力における前回までの最小方位残差量と比較することにより、最適解及び最適解での方位残差量を求める処理であり、求められた方位残差量を分析結果として、分析結果データベース７５０に受け渡す（ステップＳ０１１）。
【０１３７】
（14）仮定固有周波数計算部１７０は、ステップＳ００６で順次変更された各仮定的速における最適解から周波数残差量を最小とする最適固有周波数を計算し、周波数残差分析部１８０に受け渡す（ステップＳ０１２）。
【０１３８】
（15）周波数残差分析部１８０は、各仮定的速における最適解及びステップＳ０１２の処理で計算された最適固有周波数を基に、周波数残差及び周波数残差量を求める。このようにして求められた周波数残差量及び最適固有周波数は、分析結果データベース７５０に受け渡される（ステップＳ０１３）。
【０１３９】
（16）前述までの処理で全ての仮定速力における最適解及び最適解における方位残差量、周波数残差量を求められるので、その後、最適解判定部１９０は、各仮定速力における方位残差量と周波数残差量との和である総残差量を評価指標として、総最適解を求める。このようにして求められた総最適解（最適速力、最適距離、最適針路、最適固有周波数）及び総最適解における方位残差量、周波数残差量、総残差量は、表示制御部１９５を介して分析結果表示部２１０に受け渡される（ステップＳ０１４）。
【０１４０】
（17）分析結果表示部２１０は、ステップＳ０１４の処理で受け付けた総最適解の分析結果を数値情報として表示し、最適解計算装置１００での処理を終了する（ステップＳ０１５）。
【０１４１】
前述した本発明の実施形態における処理は、処理プログラムとして構成することができ、この処理プログラムは、ＨＤ、ＤＡＴ、ＦＤ、ＭＯ、ＤＶＤ−ＲＯＭ、ＣＤ−ＲＯＭ等の記録媒体に格納して提供することができる。
【０１４２】
前述した本発明の実施形態によれば、運動パラメータの一部を自動的に変更し、幾何学条件とドップラー効果とから他の運動パラメータ、固有周波数、方位残差量、周波数残差量を電子計算機を用いて計算することにより、無限に存在する仮定運動パラメータの組み合わせから最適解を高速に、かつ、一意的に求めることができる。
【０１４３】
【発明の効果】
以上説明したように本発明によれば、移動発音体の運動パラメータ及び固有周波数の最適解を高速かつ正確に求めることができる。
【図面の簡単な説明】
【図１】本発明の一実施形態による移動発音体の運動分析支援装置の構成例を示すブロック図である。
【図２】平滑処理部で実施する一次平滑の原理を説明する図である。
【図３】仮定運動パターンを求める方法の原理を説明する図である。
【図４】速力における距離変更範囲の最大値について説明する図である。
【図５】速力における距離変更範囲の最小値について説明する図である。
【図６】方位残差算出の原理を説明する図である。
【図７】仮定運動パターンにおける聴音器と移動発音体との相対運動を説明する図である。
【図８】平滑方位Ｂy 及び聴音器針路方向を時間に対してプロットすることにより表示させた様子を示す図である。
【図９】平滑周波数Ｆを時刻に対してプロットして表示させた様子を示す図である。
【図１０】最適解計算装置１００の処理動作を説明するフローチャートである。
【符号の説明】
１００最適解計算装置
１１０平滑処理部
１２０基準データ指定部
１２５速力変更幅指定部
１３０仮定速力変更部
１４０距離範囲計算部
１４５仮定距離変更部
１５０仮定運動パターン計算部
１６０方位残差分析部
１７０仮定固有周波数計算部
１８０周波数残差分析部
１９０最適解判定部
１９５表示制御部
２００表示装置
２１０分析結果表示部
２２０履歴情報表示部
３００オペレータ入力装置
３１０テンキー
３２０スイッチ
３３０ポインティングディバイス
４００信号処理装置
５００聴音器
６００移動発音体
７１０観測方位データベース
７２０観測周波数データベース
７３０分析用方位データベース
７４０分析用周波数データベース
７５０分析結果データベース

Claims

聴音器を用いて離散的に観測される、移動発音体の観測方位情報と観測周波数情報とから、前記移動発音体の運動パラメータと固有周波数とを分析する移動発音体の運動分析方法において、前記観測された移動発音体の前記観測方位情報と前記観測周波数情報とを入力された時間間隔で平滑化するステップと、前記平滑後の方位情報に基づく分析の基準データとなる方位情報の指定を受け付けるステップと、入力された仮定速力を入力された仮定速力の変更幅で順次変更するステップと、前記順次変更された仮定速力と基準データとから距離の変更範囲を計算するステップと、前回の仮定速力の分析結果から最適距離を予測し、距離変更範囲を修正するステップと、入力された仮定距離を入力された仮定距離の変更幅で順次変更するステップと、前記変更された仮定速力、前記変更された仮定距離、前記基準データから順次仮定運動パターンを計算するステップと、算出した前記仮定運動パターンから前記移動発音体の所定の時刻の方位情報を順次算出し、算出されたそれぞれの方位情報と同時刻の平滑後の前記方位情報との方位差分及び方位差分の二乗和である方位残差量を求めるステップと、前記各仮定速力における周波数残差量を最小とする最適固有周波数を計算するステップと、前記計算された固有周波数及び前記仮定運動パターンから前記移動発音体の所定の時刻の周波数情報を順次算出し、算出されたそれぞれの周波数情報と同時刻の平滑後の前記周波数情報との周波数差分及び周波数差分の二乗和である周波数残差量を求めるステップと、前記方位残差量と前記周波数残差量との和から最適運動パラメータを求めるステップと、を有することを特徴とする移動発音体の運動分析方法。
聴音器を用いて離散的に観測される、移動発音体の観測方位情報と観測周波数情報とから、前記移動発音体の運動パラメータと固有周波数とを分析する移動発音体の運動分析装置において、オペレータの入力を受け付ける入力受付手段と、前記離散的に観測された複数の観測方位情報及び観測周波数情報を記録する第１の記憶手段と、前記入力受付手段において受け付けた変更幅で前記運動パラメータを自動的に変更する手段と、前記運動パラメータを用い、前記移動発音体の運動分析を行う運動分析手段と、前記運動分析手段において分析した結果を表示する表示手段とを備え、前記運動分析手段は、前記複数の観測方位情報及び前記複数の観測周波数情報を平滑化する平滑処理手段と、前記平滑処理手段で平滑化した情報と自動的に変更される前記運動パラメータとを用いて、前記移動発音体の運動パターンを計算する運動パターン計算手段と、前記運動パターン計算手段で計算された運動パターンから算出した前記移動体の方位情報と、前記平滑手段により平滑化された前記方位情報との各々の差分及び差分の二乗和である方位残差量を求める方位残差計算手段と、方位残差量を最小とする運動パターンにおいて周波数の残差量を最小とする固有周波数を計算する手段と、その固有周波数と前記運動パターンとから算出した前記移動体の周波数情報と、前記平滑手段により平滑化された前記周波数情報とのそれぞれの差分及び差分の二乗和である周波数残差量を求める周波数残差計算手段と、前記方位残差量と前記周波数残差量との和を評価指標として最適な運動パラメータ及び固有周波数を求める最適パラメータ判定手段とを備えて構成されることを特徴とする移動発音体の運動分析装置。
前記最適パラメータ判定手段により求められた最適な運動パラメータ及び固有周波数を分析結果として前記表示手段に表示させる手段をさらに備えることを特徴とする請求項２記載の移動発音体の運動分析装置。