JP2011527428A - データの周波数分析のための方法および装置 - Google Patents

Abstract

【課題】データの周波数分析のための方法および装置を実現する。
【解決手段】データの周波数を分析する方法は、第１のセンサから供給された信号を入力するステップ（310）と、少なくとも１つの第２のセンサから供給された信号を入力するステップであって、第２のセンサのそれぞれから供給された信号が第１のセンサから供給された信号と強く相関するように、第２のセンサのそれぞれが第１のセンサの近くに配置されるステップ（315）と、第１のセンサからの信号と第２のセンサのそれぞれからの信号との集合から作成された伝達関数またはモデルを、センサのそれぞれについて、推定するステップと、推定したモデルのそれぞれからシステムの構造的特性を抽出するステップ（320）と、を有する。
【選択図】図３

Description

本発明は、データの周波数分析のための方法および装置に関する。本発明は、特に、航空機の飛行領域拡大試験のデータ分析に適用される。

本発明は、特に、航空分野では、飛行の制御（例えば飛行中の構造体の振動モードの分析および制御）に適用され、自動車の分野では、車体の振動の調査研究および制御に適用され、電気力学（電気を発生させる機械の制御）に適用され、核の分野では、特に炉芯の振動の監視に適用され、機械学（可動部材の調査研究および制御）に適用され、地震学（石油探査で用いられる信号の調査研究）に適用され、動物学（動物が発する音の調査研究）に適用される。

本発明の目的は、試験中（航空機の場合には飛行中）に機体の特徴、特に、共鳴周波数およびスペクトル特性を推定（estimer）することである。言い換えると、機体に取り付けた複数のセンサから来る極めて大量の情報から、非常に迅速に、それどころかリアルタイムで、適切なサインを抽出しなければならない。

このようにして決定したシステムの特徴の集合により、設計者は、システムの構造を改善して、快適さを向上させたり、飛行領域を広げたり、燃料消費を向上させたりすることができる。

分析を所望する信号は、周波数および振幅が時間経過とともに変化する可能性のある１つまたは複数のサイン曲線形の信号にノイズが加わったもので構成されている。これらの周波数および振幅をリアルタイムで推定しなければならない。

空気力、慣性力および弾性力の間の相互作用を調査研究する空力弾性学の分野において、フラッタ現象は、航空機の構造体（翼、胴体、尾翼など）の振動が不安定になることであり、非常に危険である。なぜなら、振動が不安定になるとこの構造体の完全性に影響を及ぼし、構造体を損傷して破壊にまで至らせる可能性があるからである。航空機の構造体の性質が異なる２つ以上の運動の位相が適度にずれて結合すると、空気力学的な力がシステムにエネルギーを供給することがある。すると安定な現象は不安定な現象に変化し、エネルギーをこれ以上放散しない。これがモード結合である。最もよく引用される例は、撓みモードと捩じりモードとの間で結合が起こり、位相が１／４ずれて変位と同じ方向の揚力という空気力学的な力が生じて発散振動が起こることから生じる翼のフラッタである。

フラッタ現象を特徴づけるのに複数のパラメータが影響を及ぼす。例えば、構造体の重量、剛性、形状や、速度などの動作条件である。この現象に備えるため、航空機設計者は、この現象を調査研究し、この現象が存在する場合には、その出現閾値が動作最高速度（の１５％増し）の外に位置することを証明せねばならない。最初に航空機の構造体の風洞試験を実施した後、地上での振動試験を行なう。それから理論的調査研究によってフラッタ現象のない領域を規定し、規定した領域から、航空機の構造体を「励振（exciter）」しながら飛行領域の拡大全体が少しずつ実行されるであろう。

モーダルパラメータを特定する方法を利用し、ほぼリアルタイムで周波数および減衰値を抽出し、飛行領域におけるその変化を調べる。フラッタ試験から供給される時間のデータの分析は複雑である。データにはノイズが混ざっているため、信号処理（特にフィルタリングおよびサブサンプリング）によって成形しなければならない。現在では、空力弾性構造体のモーダルパラメータを「全体的」かつ自動的に抽出するのに複数のセンサが用いられている。

国際公開第０３／００５２２９号パンフレットが公知であり、センサから供給される信号の周波数分析システムが記載されている。しかし、この分析の精度は限定されている。

本発明はこれらの欠点を解決することを目的とする。

この効果のため、本発明は、第１の形態によれば、
第１のセンサから供給された信号を入力するステップと、
少なくとも１つの第２のセンサから供給された信号を入力するステップであって、第２のセンサのそれぞれから供給された信号が第１のセンサから供給された信号と強く相関するように、第２のセンサのそれぞれが第１のセンサの近くに配置されるステップと、
第１のセンサからの信号と第２のセンサのそれぞれからの信号との集合から作成された伝達関数またはモデルを、センサのそれぞれについて、推定するステップと、
推定したモデルのそれぞれからシステムの構造的特性を抽出するステップと、を有することを特徴とする、データの周波数分析を行う方法を対象とする。

１つのモデルにおいて、センサから供給される信号を、白色雑音による励振を受けたフィルタの出力と考える。構造的特性には、例えばスペクトル特性、周波数、振幅、初期の位相、減衰係数、モードが含まれる。

したがって、構造モードを表わすこのモデルは、線形であると見なされる。本発明を実現することによって、周波数／減衰係数の組をオンラインで追跡しつつ、リアルタイムで処理が実行される。本発明によってシステムの構造的特性がリアルタイムで処理されるため、リアルタイムで確実にシステム（例えば、航空機）の動作が期待通りになるようにすることができる。そのため、利用される分析法を改善し、時間の節約、したがってコストの節減という益々増大する制約に応じることができる。

特定の特徴によれば、抽出するステップにおいて、センサから供給された信号は、多項式と見なされる。これらの構成によって、多項式の係数の数が実現される信号のサンプル数よりもはるかに少なくなるため、信号の表示は、コンパクトになる。

特定の特徴によれば、抽出するステップにおいて、センサ間で、時間および空間に関してゆっくりと変化する係数を持つ線形回帰方程式を解く。

特定の特徴によれば、抽出するステップは、
時間、次数（ordre）およびセンサの空間に関して回帰式適応モデルを作成するステップと、
適応モデルを作成するステップの結果に応じて、それぞれの次数についてモードを推定するステップと、を有する。

特定の特徴によれば、センサから供給された信号を入力するステップのそれぞれは、適応モデルを作成するステップの前に、センサから供給された信号のノイズレベルをリアルタイムで低減するステップを有する。

特定の特徴によれば、モードを推定するステップは、適応モデルを作成するステップの結果に応じて、モデルのパラメータを抽出するステップを有する。

特定の特徴によれば、モデルのパラメータを抽出するステップは、次数がNで次元（dimension）がセンサの数に等しい多項式行列の逆行列を求めるステップを有する。

特定の特徴によれば、モードを推定するステップは、推定したモードの分散を小さくすることのできる冗長な情報の集合を構成するモデルのそれぞれに関するパラメータを供給するように適応される。

特定の特徴によれば、適応モデルを作成するステップにおいて、ＡＲＭＡタイプのモデルを作成する。

特定の特徴によれば、センサのそれぞれおよび考慮されるすべての次数について、各瞬間でＡＲＭＡタイプのモデルを作成する。

特定の特徴によれば、推定するステップは、主対角線上にセンサのそれぞれのパワースペクトル密度が表され、他の項にスペクトル間要素が表される対称スペクトル間行列である多項式行列の逆行列を求めるステップを有する。

特定の特徴によれば、適応モデルを作成するステップは、瞬間nについて、次数N＝0で初期化したモデルの次数N＝［1, 2, ..., Nmax］に関する
後ろ向きおよび前向き線形予測の誤差ベクトルを計算するステップと、
前向きおよび後ろ向き部分相関行列を計算するステップと、
後ろ向きおよび前向き線形予測の誤差共分散行列を計算するステップと、
前向き線形予測の誤差パワー行列を計算するステップと、
αはスカラであり、λは忘却因子であり、

である場合、利得ベクトル

を直接計算するステップと、
既知の

から反復によってベクトル

を計算するステップと、
k＝1〜Nについて、モデルを表わす行列A_kを計算ステップと、の再帰手続きを有する、時間に関する再帰手続きを有する。

特定の特徴によれば、本発明の対象の方法は、
同一のモデルからは、クラスごとにモードが１つだけ供給されることと、
推定（estimation）は、それぞれの推定の導出元とは独立にすべて同じ重みを有することと、の２つの制約のいずれか一方を実現し、モードを推定するステップから供給されたモードを分類するステップを有する。

特定の特徴によれば、信号は、航空機の構造体の加速度を表す。

第２の態様によれば、本発明は、情報システムにインストール可能であって、上で簡単に説明した方法を実現できる命令を有するコンピュータプログラムを対象とする。

このプログラムの特別な利点、目的、特徴は、本発明が対象とする、上で簡単に説明した方法の利点、目的、特徴と同様であるため、ここで改めて記載しない。

本発明の他の利点、目的、特徴は、添付の図面を参照した以下の説明で理解されるであろう。なおこの説明は例示を目的としたものであり、本発明がこの説明に限定されることはない。

本発明が対象とする方法を実現できる装置を備えた航空機の概略図である。 図１に示した装置の２つのセンサから供給される信号を表わす図である。 本発明が対象とする方法の第１の実施態様で実行されるステップを表すフローチャートである。 本発明が対象とする方法の第２の実施態様で実行されるステップを表すフローチャートである。 図４に示したステップのうちの１つを実行しているときのフィルタの構成を表す図である。 ノイズ低減システムの一実施態様で順番に実行される機能の概略を表す図である。 センサから供給され、図４に示した実施態様のアルゴリズムの入力となる各瞬間におけるサンプルを表す概略図である。 図４に示した第２の実施態様で実現される再帰手続きの概略を表す図である。 「ダイナミック雲」タイプの非監視分類法におけるクラスの展開を表す概略図である。 図４に示した第２の実施態様で実現される確認ウインドウを示す図である。

図１は、翼１２０よりも前方にあって互いに近接した２つのセンサ１１０および１１５と、翼１２０よりも後方にあって互いに近接した２つのセンサ１２５および１３０とが取り付けられた航空機１０５を示す。

説明のため、図１には、互いに近接した２組のセンサだけを示してある。しかし、本発明を実現する一実施形態では、３組以上が利用されることに留意されたい。

「互いに近接した」という用語は、ここでは、互いに強く相関した信号を受け取るセンサを意味する。

互いに近接したセンサは、時間的にずれていて減衰状態が異なるが、ほぼ線形の伝達関数に従ったほぼ同じ周波数を受け取る。ここで問題となるセンサは、例えば加速度計である。

図２から、１組のセンサの第１のセンサから供給された信号２０５は、ノイズ２１０と、２つのピーク２１５および２２０とを有すること、そして同じ組の第２のセンサから供給された信号２５５は、ノイズ２６０と、２つのピーク２６５および２７０とを有することが分かる。ピーク２６５は、時間経過とともに減衰してずれていくピーク２１５に対応する。ピーク２２０は、時間経過とともに減衰してずれていくピーク２７０に対応する。

容易に分かるように、互いに近接した少なくとも２つのセンサの集合（ここでは組）を実現することで、本発明により航空機の構造的特性を分析することができる。構造的特性には、例えば、スペクトル特性、周波数、振幅、元の位相、減衰、モードが含まれる。

図３から、第１の実施態様において、本発明が対象とする周波数分析法は、最初に、複数のセンサからなる各グループを、振動する傾向がある機械システムの構造体上に配置するステップ３０５を有することが分かる。センサの各グループは、少なくとも１つの「第２の」センサが「第１の」センサの近傍に配置される。

この機械システムの動作中に、１つのグループのセンサの第１のセンサから供給された信号が入力されるステップ３１０と、同じグループのセンサ群の少なくとも１つの第２のセンサから供給された信号が入力されるステップ３１５とを実行する。センサから供給された信号が入力される各ステップは、そのセンサから供給された信号のノイズレベルを低減させるステップを有する。このノイズの低減は、公知のようにしてセンサごとに実行するか、または、１つのセンサから供給された１つの信号を各成分として含むベクトルに関して実行することができる。このノイズ除去機能は、ウェーブレット基底上の分解によって確実に行われることが好ましい（ステファンマラーのアルゴリズム）。ステップ３１０および３１５では、例えば、16Hz未満の固有周波数を探索する構造体について、センサから供給された、例えば、図６の256Hzなどの周波数がはるかに大きい信号をサンプリングする。ウェーブレットを用いると、簡単かつ迅速な処理が可能になることに留意されたい。

ステップ３２０では、第１のセンサおよび第２の各センサから供給された信号を処理することによって伝達関数を抽出する。そのとき、構造モードを表すモデルは線形であると見なし、入力すなわち励振を考慮しないで周波数／減衰の組をオンラインで追跡する。

ステップ３２０では、センサから供給された信号を多項式と見なし、センサから供給された信号から情報を抽出する。このとき、異なるセンサ群の異なるセンサから供給された信号間の位相および上記情報の値の修正を考慮する。

ステップ３２０では、十分に安定した時間区間でモデルを推定できるようにするため、時間経過とともにゆっくりと変化する係数を持つ線形回帰方程式を解くステップ３２５を実行する。

ステップ３２５は、
−（各瞬間における）時間と（考慮される各次数に関する）次数（ordre）と（各センサに関する）センサ空間とに関して、ＡＲＭＡタイプの回帰式適応モデルを作成するステップ３３０と、
−適応モデルを作成するステップの結果に応じて、各次数に関するモードを推定するステップ３５０と、を有する。

ステップ３３０では、パラメータタイプのモデルとして自己回帰移動平均（autoregresive a moyenne ajustee；ＡＲＭＡ）タイプのモデルを作成する。センサから供給されたサンプルは、集められて１つのベクトルになる。このベクトルの成分の数は、考慮されるセンサの数である（図７を参照）。例えば、センサが４つ存在する場合は、考慮されるベクトルは、４次元（dimension）である。より一般的に、以下の説明では、「p」を、互いに近接したセンサからなる１つのグループに含まれるセンサの個数とする。

ステップ３３０では、「時間に関する再帰」モデルと呼ばれるモデルを作成する。なぜならこのステップでは、得られた最新の推定を用いてモデルのパラメータをアップデートするからである。したがって、この実施態様では、連続した２つの瞬間における関係は、相関しておりコヒーレントであるため、本発明によって両者間が関係付けられる。さまざまな実施態様において、時間に関する再帰手続きに利用できるのは、直前の２つの瞬間である。

モデル化の最適次数Nは必ずしも既知でなくてもよい。好ましくは、この最適次数Nを決定しないで、リアルタイムで結果を得る。対称的には、十分に大きい次数Nmaxまで再帰手続きを実行する。すると十分に多数のモデルの集合が得られるので、その中で探索している情報（構造モード）が表される。

図３に示した実施態様では、ステップ３３０は、主対角線上に各センサのパワースペクトル密度が表され、対角線以外の項にスペクトル間要素が表された対称なスペクトル間行列を処理するステップを有するステップ３３５を有する。

ステップ３３５は、瞬間nについて、次数N＝0で初期化してモデルの次数N＝［1, 2, ..., Nmax］に関して以下のステップの再帰手続きを行う、時間に関する再帰手続きを行うステップ３４０を有する。
−後ろ向きおよび前向き線形予測の誤差ベクトルを計算するステップ。
−前向きおよび後ろ向き部分相関行列を計算するステップ。
−後ろ向きおよび前向き線形予測の誤差共分散行列を計算するステップ。
−前向き線形予測の誤差パワー行列を計算するステップ。
−αはスカラであり、λは忘却係数であり、

である場合、次元がセンサの数である利得ベクトル

を直接計算するステップ。
−既知の

から反復によってベクトル

を計算するステップ。
−k＝1〜Nについてモデルを表わす行列A_kを計算するステップ。

ステップ３３０が終わったときに得られるモデルの集合から選別して構造モードを抽出する。モデルは、多項式で示されると考えられる（「自己回帰移動平均」の略号である「ＡＲＭＡ」と呼ばれるモデル）。構造モードを表わすパラメータ値を抽出するため、各瞬間および各次数に関して共分散行列E_aおよびE_bを最小化する。そのために、行列E_aの逆行列とベクトルε_aとの積および行列E_bの逆行列とベクトルε_bとの積を、それぞれ表すベクトルθ_aおよびθ_bを計算する。

これらの計算をリアルタイムで実行するため、1〜Nmaxの値すべてについてベクトル

を直接計算する。次に、

を知得したので、反復によって

を計算する。

1〜Nmaxのすべての次数について瞬間nで推定するモデルの係数行列である行列A_k（p次元の正方行列）の数は、Nmax×（Nmax−1）／2に等しく、N＝15であれば105である。

構造モードを分類するための原理は、所定の次数における１つのパラメータが適切なものであれば、そのパラメータは、より大きな次数で再発見されるであろうと考えることである。パラメータの分類は、監視されることなく、異なるモデルから供給された類似のオブジェクトを探索することからなる。後述するように、この分類のために軌跡を構築する。

ステップ３５０では、各モデルについて、推定したモードの分散を小さくすることのできる冗長な情報の集合を構成するパラメータを供給する。ステップ３５０は、適応モデルを作成するステップの結果に応じて、モデルのパラメータを抽出するステップ３５５を有する。ステップ３５５は、次数Nで次元がセンサの数に等しい多項式行列の逆行列を求めるステップ３６０を有する。

ステップ３６０は、コレスキ分解を行なうステップ３６５を有する。

本発明が対象とする方法の第２の特定の実施態様は、図４に記載した飛行領域拡大試験のデータを分析するためのリアルタイム手続きに従う。この実施態様により、構造体への励振を考慮せずに、次の情報が出現する前に各情報を処理することができる。

この方法の特徴は、全体として、冒頭に述べた安全性およびコスト削減の制約に応じており、この方法の特徴により、より高性能のモード分析を提供し飛行領域拡大の手続きを改良することができる。

分析する信号は、航空機の主構造体上で測定された加速度測定値である。

図４には、p個の加速度タイプの信号に関するリアルタイム分析手続きの各操作が記載されており、その各操作について以下の段落で詳細に説明する。

この分析法の第１のステップ４０５は、互いに近接したセンサからなる同じグループのp個のセンサから供給されたp個の信号の「ノイズ除去」を、例えば、Ｓ．マラーによって提案されたピラミッドアルゴリズムによって、ウェーブレットの正規直交基底を用いて実行することからなる。このアルゴリズムは、視覚および画像圧縮に関する1983年のBurtおよびAdelsonの研究に由来する。このアルゴリズムは実行するのが非常に簡単であり、計算量は処理するサンプルの数に比例する。

処理が非線形であることに加え、フィルタバンクは原因を示さないために実行が複雑になる。フィルタは、ウェーブレット変換を実行し、図５に示すように構成される。

マラーのアルゴリズムは常にリアルタイムであり、一般に、各瞬間にp個のサンプルを同時に処理する。考慮される原理は、信号のオーバサンプリングにより分析および合成による「ノイズ除去」を可能にする。

そこで、ステップ４０５では、最初に、信号をサブバンドに分解するステップ４０６を実行する。この分析すなわち分解の処理中、連続する信号のサンプリングによって形成される信号列は、最初の段階では、この信号の離散化に関する一定のスケールへの近似値であるとみなされる（サンプリングは、実際には最も細かいスケーリングに対応する）。

この近似値のスケールは、j＝0に対応することとする。したがって、部分空間V₀に属するデータから出発する。このとき、一連のサンプル｛x₀, x₁, ..., x_k, ...｝＝〈f | φ_0,k〉、k∈Zは、分析したいデータの集合を構成する。近似の部分空間同士の関係は、

である。したがって２つの部分空間V₁およびW₁上で離散信号すなわちサンプル列を分解し、分解能2^-1のサンプル列すなわち〈f | ψ_1,k〉および〈f | φ_1,k〉を得れば十分である。以下の関係によって、２つの連続した分解能に関する再帰的手続きが行われる。

これらの式は、それぞれ、畳み込みの積g_-n＊〈f | φ_j-1,n〉およびh_-n＊〈f | φ_j-1,n〉に続いて２でデシメーションしたものを表わす。したがって、スケールjにおける近似の係数〈f | φ_j,k〉および詳細の係数〈f | ψ_j,k〉は、スケールj-1で得られた係数から、フィルタ

による単なるフィルタリング操作に続いてデシメーションを行なうことによって計算される。

次に、分解係数の閾値処理というステップ４０７を実行する。このステップでは、第１のサブバンドの係数しか保持しない。

次に、合成、すなわち再構築のステップ４０８を実行する。このステップでは、近似の部分空間上への射影に関する知識から、信号が以下のように再構築される。

すなわち、

である。

再構築は、前の操作の双対操作である。再構築は、分解によって得られた近似係数および詳細係数における補間の後にデジタルフィルタリングを実行することによって得られる。

したがって、「ノイズ除去」の一般的構造は、図６に示したものである。考慮される最大周波数はサンプリング周波数fの半分であることに留意されたい。求める信号は、未処理の信号中で最も小さい周波数に位置する。信号を0〜f／16、f／16〜f／8、f／8〜f／4、f／4〜f／2の周波数帯に分解した後、ウェーブレットの係数をゼロにして閾値処理を行なう。次に合成を実行し、ノイズが除去された探索する信号を含む信号を得る。

次にステップ４１０では、ＡＲＭＡ（「自己回帰移動平均」）タイプの適応モデル化を実行する。信号のモデルを作成するこのステップは、パラメータタイプで実行される。このタイプによって、調査研究する信号のスペクトル分析と、時間、次数、センサの空間に関する再帰手続きとが可能になる。

ステップ４１０は、前向き線形予測ベクトルおよび後ろ向き線形予測ベクトルを決定するステップ４１１を有する。ベクトルの場合に拡張すると、予測ベクトルを決定するごとに、図７に示した表現を用いて、サンプルの最新のＮ個のベクトルの線形結合で

を表わす。この図７には、ｐ個のセンサの空間上にベクトル形式の前向き線形予測が示されている。

したがって、時間空間内で信号を直接モデル化しなければならない。x _nがサンプル列を表わすとすると、モデルは以下のようになる。

ベクトルx _nは、p個のセンサから来たp個の現サンプルを成分として持つことに留意されたい。この表現では、A_kは、センサの数に対応するp次元の行列である。なおNは、モデル化の次数である。

スカラ法を一般化することにより、誤差およびその予測に関する新しい表現式が与えられる。

すなわち、

である。

外部から見ると、e _nは、サンプルのベクトル列x _nによって励振されたフィルタFIRの出力である。線形性があるため手続きを逆転させることができる。したがってx _nはe _nによって励振されたフィルタの出力として現われる。このフィルタは、多項式行列Φ(z)の逆行列を求めることによって得られる。このフィルタは、安定でありIIRタイプ（「無限パルス応答」の略号）である。

後ろ向き線形予測の決定に関しては、x _n-Nの推定

は、ベクトル｛x _n-N+k｝_k=1,...,Nの線形結合として表現される。その特性は、前向き線形予測に関して得られた特性と同様である。

次にステップ４１２では、共分散行列の計算と、前向きおよび後ろ向き部分相関行列の計算と、前向きおよび後ろ向き予測の誤差パワー行列の計算とを実行する。

次に、以前に規定した前向きおよび後ろ向き線形予測の概念から構築された演繹法ならびに帰納法による前向きおよび後ろ向き線形予測の誤差共分散行列を最小化する。帰納法による前向き線形予測の誤差共分散行列は、以下のように規定される。

ただし、0<<λ≦1は、忘却因子または適応因子と呼ばれる。

ステップ４１５では、Nのすべての値についてN次モデルの係数行列が得られる。そのためステップ４１６では、共分散行列

が最小であるとき、すなわち、

であるとき、最適なベクトル

が得られる。

ステップ４１７では、演繹法による後ろ向き線形予測の誤差共分散行列の推定に基づいた同様の推論により、以下の式が得られる。

「ＥＳＡ」というアルゴリズムでは、マルチセンサの拡張の枠内で、以下に規定する２つのベクトルの計算を実行する。

アルゴリズム「ＥＳＡ」によって直接計算されるベクトル

は、１つの行列と１つのベクトルの積である。上記１つの行列は、ベクトル

のダイアディック積を瞬間n+1まで重みを付けて足し合わせたものであり、上記１つのベクトルは、瞬間n+1におけるこの同じベクトル、すなわち

である。

ベクトル

は、

から反復によって得られる。

次に、ステップ４２０では、逆の問題を解く。すなわち、モデルのパラメータを抽出する。マルチ-センサの場合に拡張したモデルに関する上記の式（１）および（２）から、以下のように表現される。

すなわち、

である。

ただし、

であり、これは、p×p次元で次数N＝［1, Nmax］のzに関する多項式行列である。なお、pは、分析するセンサの数である。

e(z)は、未知の誤差ベクトルだが、その共分散行列は推定可能である。したがって、その行列を明らかにするため、
e(z) ^t' e(z)＝Φ(z)x(z) ^t x(z)^tΦ(z)
を計算する。

すなわち、スペクトル間行列は、
x(z) ^t x(z)^t＝Φ(z)^-1 e(z) ^t' e(z)Φ(z) である。

この行列は対称であり、主対角線上に各センサのパワー密度スペクトルが現われる。対角線以外の項は、スペクトル間要素である。モデルのパラメータ抽出は、それぞれの次数について各瞬間にN次かつp次元の多項式行列Φ(z)の逆行列を求めることからなる。

そのため、行列e(z) ^t' e(z)＝L^t.Lのコレスキ分解を行なうと、
x(z) ^t x(z)＝Φ(z)^-1L(z)^tL(z)^tΦ(z)^-1＝(Φ(z)^-1L(z))^t(Φ(z)^-1L(z))
= (Λ(z)^-1U(z)L(z))^t(Λ(z)^-1U(z)L(z)) となる。

ただし、

多項式Fは、さまざまな伝達関数の分子を表わし、分母（重ね合わせの定理によれば実際には単一の分母）は、固有値である。

と置くと、

となる。

分子の係数は、（符号なしの）スペクトル間行列のN+1番目の列のN+1個の要素の平方根である。

適応モデル化の処理は、図８のように表わすことができる。このマルチセンサのモデル化処理に関する図では、３つの再帰手続き８０５、８１０および８１５が入れ子になっていることが分かる。

次にステップ４２５では、モードを分類する。各瞬間に次元｛1：p：N×p｝のN個のモデルを推定する（ただし、pは、センサの数である）。各モデルから抽出されたパラメータ（周波数および減衰の値）は、冗長な情報の集合であり、推定の分散を小さくすることができる。ここで、一定の基準という意味で以下の制約のもとで、モデルから供給されたデータを分類する。
−同一のモデルからは、クラスごとにモードが１つだけ供給される。
−推定は、それぞれの推定の導出元とは独立にすべて同じ重みを有する。

演繹法で（特に密度を）知得できないため、「ダイナミック雲（nuees dynamiques）」タイプの非監視分類法が好適である。この方法は、ラベルのないデータを集合させるための自然（naturel）（暗黙（implicite））クラスを見つけることからなる。

図９に示す方法は、以下のアルゴリズムによる要求の集合を十分に満たす。
−ステップ４２６では、最初にどのようにK個のクラスに分割するかを選択する。
−ステップ４２７では、各データについて最適のクラスをシステム的に探す。（例えば、ユークリッド距離またはカルバック-ライブラ距離を用いて）データから重心までの距離を計算し、クラスの中心が最も近くにあるクラスにデータを割り当てる。
−ステップ４２８では、「発信」クラスおよび「受信」クラスの重心をアップデートする。

次に、収束するまでステップ４２７および４２８に戻る。

ステップ４３０では、モードの軌跡を構築する。

問題は、個数が時間とともに変化するモードの周波数に対応するターゲットの集合の軌跡をリアルタイムで追跡することである。アルゴリズムの構造は、追跡する経路ごとに１つのカルマンフィルタの周囲に構成される。「測定値」の集合がモデルによって供給されるとき、それらの測定値を存在する経路と相関させる。その目的は、ここでは、受け取った測定値の中から、測定値を予測する元になったターゲットからのものである可能性のある測定値を選択することである。しばしば利用される原理は、なされた予測の周囲に「ゲーティング」と一般に呼ばれるウインドウを規定することからなる。

処理の一般的なアーキテクチャは、以下の原理に基づいている。
−既知の軌跡から瞬間nが分かっているため、瞬間n+1における状態を予測する。
−測定値−ターゲットの関連付けは、測定値を、既知の軌跡から予測される値と比較することからなる。この処理により、すでに存在する経路が維持されるだけでなく、新しい経路を開始するとともに、場合によっては観察空間から出たターゲットに対応する経路を除去することもできるはずである。経路追跡の性能は、これらの関数の品質に依存する。
−状態、カルマンフィルタの利得、共分散行列をアップデートしてフィルタリングする。

この手続きにおいて、カルマンフィルタにより、予測が基本的な役割を果たす複数のターゲットの追跡が可能である。このフィルタは、各ターゲットについて、分散の最小値という意味で状態に関するフィルタリングされた推定を供給し、状態を予測し、「ゲーティング」を計算することができる。

解は、予測方程式とフィルタリング方程式という２つの体系の集合によって構成される。すなわち、以下の表に示したようになる。

確認ウインドウにより、各ターゲットについて、ターゲットに属する可能性のある測定値を選択することができる。その原理は、予測される測定値の周囲に１つの領域（観測空間内における体積）を規定することである。この領域のサイズは、予測される測定値の統計的性質（ここではガウス型）によって規定される。一般に、この体積の「次元（dimension）」は、正しく選択されなければならない。実際、測定値の選別と、ターゲットを導出元とする測定値がこの体積の範囲を定める表面の内側にある確率とは、この体積に依存する。

図１０は、確認ウインドウを示す図である。

測定値−ターゲットを関連付ける技術は、ターゲットを追跡する手続きの中心的な部分である。多くの技術が存在しており、その中のいくつかは、経路の出現および消失を管理しない。したがって、この管理を実現するには追加メカニズムを考えなければならない。簡単な１つの方法は、以下の規則を採用することである。
−規則１：既存のいかなる経路とも関係しないあらゆる測定は、新しい１つの経路の開始と見なされる。
−規則２：１つの経路は、連続した少なくともNd回の測定と関連付けられた場合、確認（検出）される。
−規則３：１つの経路は、連続した少なくともNl回の測定と関連付けられない場合、消失したと見なされる。

このハンガリ法によって、以下の特別な解法を用いて最小コストを求めることにより、推定値−測定値を割り当てる問題を解決することができる。

m個の点に割り当てるm個の資源があり、Cが割り当てコストの行列であるとする。任意の１つの割り当ては、（x, y, ...u）を用いて（1, x）、（2, y）、..（k, t）、...（m, u）と表記されるm個の組によって規定される。すなわち｛1, 2, ..., m｝の順列である。特定の割り当ては、消費またはコストの合計：
F(x, y, ..t, .., u)＝C1,x＋C2,y＋...＋Ck,t＋ ... ＋Cm,u
に対応する。

したがって、問題は、Fが最小になるように（x, y, ..t, .., u）を決定することである。

本発明が対象とする方法を実現するため、本発明の好ましい一実施態様では、汎用コンピュータを用意する。このコンピュータには、このコンピュータにインストール可能で上に詳細に説明したステップおよびアルゴリズムを実現する命令を有するコンピュータプログラムが付随している。

Claims (15)

1. 第１のセンサから供給された信号を入力するステップと、
   少なくとも１つの第２のセンサから供給された信号を入力するステップであって、前記第２のセンサのそれぞれから供給された信号が前記第１のセンサから供給された信号と強く相関するように前記第２のセンサのそれぞれが前記第１のセンサの近くに配置されるステップと、
   前記第１のセンサからの信号と前記第２のセンサのそれぞれからの信号との集合から作成された伝達関数またはモデルを、前記センサのそれぞれについて、推定するステップと、
   前記推定したモデルのそれぞれからシステムの構造的特性を抽出するステップと、を有することを特徴とする、データの周波数分析を行う方法。
2. 前記抽出するステップにおいて、前記センサから供給された前記信号は、多項式と見なされることを特徴とする、請求項１に記載の方法。
3. 前記抽出するステップにおいて、前記センサ間で、時間および空間に関してゆっくりと変化する係数を持つ線形回帰方程式を解くことを特徴とする、請求項１または２に記載の方法。
4. 前記抽出するステップは、
   時間、次数およびセンサの空間に関して回帰式適応モデルを作成するステップと、
   前記適応モデルを作成するステップの結果に応じて、それぞれの次数についてモードを推定するステップと、を有することを特徴とする、請求項１〜３のいずれか１項に記載の方法。
5. 前記センサから供給された信号を入力するステップのそれぞれは、前記適応モデルを作成するステップの前に、前記センサから供給された信号のノイズレベルをリアルタイムで低減するステップを有することを特徴とする、請求項４に記載の方法。
6. 前記モードを推定するステップは、前記適応モデルを作成するステップの結果に応じて、モデルのパラメータを抽出するステップを有することを特徴とする、請求項４または５に記載の方法。
7. 前記モデルのパラメータを抽出するステップは、次数がNで次元がセンサの数に等しい多項式行列の逆行列を求めるステップを有することを特徴とする、請求項６に記載の方法。
8. 前記モードを推定するステップは、前記推定したモードの分散を小さくすることのできる冗長な情報の集合を構成する前記モデルのそれぞれに関するパラメータを供給するように適応されることを特徴とする、請求項４〜７のいずれか１項に記載の方法。
9. 前記適応モデルを作成するステップにおいて、パラメータタイプのモデルを作成することを特徴とする、請求項４〜８のいずれか１項に記載の方法。
10. 前記適応モデルを作成するステップにおいて、ＡＲＭＡタイプのモデルを作成することを特徴とする、請求項４〜９のいずれか１項に記載の方法。
11. 前記センサのそれぞれおよび考慮されるすべての次数について、各瞬間で前記ＡＲＭＡタイプのモデルを作成することを特徴とする、請求項１０に記載の方法。
12. 主対角線上に前記センサのそれぞれのパワースペクトル密度が表され、他の項にスペクトル間要素が表される対称スペクトル間行列である多項式行列の逆行列を求めるステップを有することを特徴とする、請求項４〜１１のいずれか１項に記載の方法。
13. 前記適応モデルを作成するステップは、瞬間nについて、前記次数N＝0で初期化した前記モデルの次数N＝［1, 2, ..., Nmax］に関する
    後ろ向きおよび前向き線形予測の誤差ベクトルを計算するステップと、
    前向きおよび後ろ向き部分相関行列を計算するステップと、
    後ろ向きおよび前向き線形予測の誤差共分散行列を計算するステップと、
    前向き線形予測の誤差パワー行列を計算するステップと、
    αはスカラであり、λは忘却因子であり、
    

    

    である場合、利得ベクトル
    

    

    を直接計算するステップと、
    

    

    を知得してから反復によってベクトル
    

    

    を計算するステップと、
    k＝1〜Nについて、前記モデルを表わす行列A_kを計算ステップと、の再帰手続きを有する、時間に関する再帰手続きを有することを特徴とする、請求項４〜１２のいずれか１項に記載の方法。
14. 同一のモデルからは、クラスごとにモードが１つだけ供給されることと、
    前記推定は、それぞれの推定の導出元とは独立にすべて同じ重みを有することと、の２つの制約のいずれか一方を実現し、前記モードを推定するステップから供給されたモードを分類するステップを有することを特徴とする、請求項１〜１３のいずれか１項に記載の方法。
15. 情報システムにインストール可能であって、請求項１〜１４のいずれか１項に記載の方法を実現できる命令を有するコンピュータプログラム。

Images