(57)【要約】【課題】時系列データのみを用い、時系列に関わるその他の情報を一切使用することなく、二つの時系列の相互相関にかかわる知見を得る。【解決手段】本発明では、解析対象の２系列の時系列データを、夫々解析して、その一般化三角多項式パラメータを求める２系統の単一系列解析手段２ａ，２ｂと、夫々の単一系列解析手段により求められた２系列の時系列データの一般化三角多項式パラメータ１６ａ，１６ｂから、２系列の時系列データの間の相互相関を定量化する諸量を求める相互相関解析手段３とから構成している。本発明を利用して構成される予測支援システムは、多数の時系列中から一対の時系列を見いだすことにより、一方の時系列の現在までの挙動から他方の時系列の将来挙動を予測することが期待される。【選択図】図２２
2010108402.tif 000002
【特許請求の範囲】
【請求項１】
解析対象の２系列の時系列データを、夫々解析して、その一般化三角多項式パラメータを求める２系統の単一系列解析手段と、夫々の単一系列解析手段により求められた２系列の時系列データの一般化三角多項式パラメータから、２系列の時系列データの間の相互相関を定量化する諸量を求める相互相関解析手段とから構成し、単一系列解析手段は、解析対象の時系列データを解析するシーケンスを決定する解析シーケンス決定手段と、解析シーケンスに従い、最大エントロピー法及び非線形最小自乗法によりデータを解析して一般化三角多項式のパラメータと残差を求めるデータ解析手段とから構成し、相互相関解析手段は、一般化三角多項式パラメータからパワースペクトル密度を求めるパワースペクトル密度計算手段と、クロススペクトル（コスペクトル、クオドスペクトル）とコヒーレンスを含む相互相関を定量化する諸量を求める相互相関諸量計算手段と、相互相関諸量計算手段により求められたコヒーレンスを理論値と比較して、それらの差が所定値以下の場合に、その差をなくすようにクロススペクトルを補正する整合性評価補正手段と、整合性評価補正手段において評価した差が所定値以上の場合に、単一系列解析手段に再解析の指令を出力し、所定値以下の場合には、相互相関諸量計算手段により計算された諸量と、パワースペクトル密度計算手段により計算されたパワースペクトル密度及び整合性評価補正手段により補正されたクロススペクトルを出力する分岐出力手段を構成したことを特徴とする２系列の時系列データの解析装置。
【請求項２】
パワースペクトル密度計算手段は一般化三角多項式パラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データのスペクトルを求める構成とすることを特徴とする請求項１に記載の２系列の時系列データの解析装置。
【請求項３】
δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用することを特徴とする請求項２に記載の２系列の時系列データの解析装置。
【請求項４】
解析シーケンス決定手段は、解析対象の時系列データを、予め設定された最大データ点数を越えないようにセグメントに分割する分割処理過程と、各セグメントに対してデータ構造の均一性を評価する評価処理過程と、この評価に基づき、不均一のセグメントデータを均一の複数のセグメントに分割する分割処理過程と、分割された均一のセグメントの夫々を基本セグメントとして、それらに対して、記憶手段に記憶領域を設定する記憶領域設定処理過程を備えていることを特徴とする請求項１に記載の２系列の時系列データの解析装置。
【請求項５】
セグメントは、解析対象の時系列データの全期間データに対応する最上位レイヤーのセグメントから、順次２分割されて構成される、下位のレイヤーの夫々に含まれる２の累乗のセグメントと、最下位のレイヤーにおいて、データ構造が変化するセグメントを、データ構造の変化している個所において分割して順次形成されるサブレイヤーの各セグメントとから構成されることを特徴とする請求項４に記載の２系列の時系列データの解析装置。
【請求項６】
データ構造の均一性は、下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移により評価することを特徴とする請求項４に記載の２系列の時系列データの解析装置。
【数１】
2010108402.tif 000025
但し、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【請求項７】
解析シーケンス決定手段には、不等間隔の時系列データを等間隔化処理する等間隔化処理過程を設けたことを特徴とする請求項４に記載の２系列の時系列データの解析装置。
【請求項８】
データ解析手段は、解析対象の時系列データから取り出したセグメントデータに対して下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移を求める計算過程と、求めた推移から最大エントロピー法によるパワースペクトル密度の計算に用いるラグ値を決定するラグ値決定過程と、決定されたラグ値により上記パワースペクトル密度を求めるパワースペクトル密度計算過程と、求められたパワースペクトル密度から主要ピークの周波数と、そのパワーを抽出する抽出過程と、抽出された周波数とパワーとから非線形最小自乗法により当てはめ曲線のパラメータを求めるパラメータ導出過程と、求められた当てはめ曲線とセグメントデータとを比較して一致度を評価し、所定の一致度以上の場合には上記当てはめ曲線のパラメータをセグメントに対応する記憶手段に記憶して次の過程に移行し、所定の一致度に達しない場合には上記当てはめ曲線のパラメータを破棄して上記ラグ値を決定する過程に移行する評価過程とを備えており、ラグ値決定過程は、設定される解析の精度に応じて、より高精度の解析において、より大きいラグ値を選択するように構成し、抽出過程は、設定される解析の精度に応じて、より高精度の解析において、より多くのピークを主要ピークとして抽出するように構成すると共に、評価過程は、設定される解析の精度に応じて、より高精度の解析において、より高い一致度となるように構成し、相互相関解析手段からの再解析指令により、高精度の解析が設定されるように構成したことを特徴とする請求項４に記載の２系列の時系列データの解析装置。
【数２】
2010108402.tif 000026
但し、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【請求項９】
データ解析手段におけるセグメントデータの解析は、最上位のレイヤーから最下位のレイヤーに向かって順次行うものとし、この際、上位のレイヤーのセグメントのセグメントデータは、セグメント当たりのデータ点数が予め設定された最大データ点数を越えないように、解析対象の時系列データのデータ点数をバンチング処理により低減してセグメントデータを構成する処理を行う構成とすると共に、各レイヤーの各セグメントに対応して一時記憶領域を構成し、上位側のセグメントに対して求められた多項式パラメータを下位の各レイヤーの各セグメントの一時記憶領域に記憶すると共に、上記多項式パラメータから構成される当てはめ曲線と、原時系列データの対応部分との差演算により求められる残差時系列データを次の下位レイヤーのセグメントデータとして設定する処理を行うことを特徴とする請求項８に記載の２系列の時系列データの解析装置。
【請求項１０】
データ解析手段により求めた解析結果を評価し、条件に適合する場合に記憶手段に記憶する解析結果評価手段と、上記パラメータから得られるパワースペクトル密度と最大エントロピー法を用いて求める残差のパワースペクトル密度を合成して時系列データの再構成スペクトルを求める再構成スペクトル計算手段と、この計算手段を利用して、上記記憶手段から取り出した上記パラメータと残差により、解析対象の時系列データに対応する一般化三角多項式を再構成する解析結果構成手段とを備えたことを特徴とする請求項１に記載の２系列の時系列データの解析装置。
【請求項１１】
解析結果構成手段は、全ての記憶領域に記憶されている基本セグメントのパラメータと残差とから全データ領域を一括して再構成スペクトルを求めて記憶する一括解析結果構成手段と、各記憶領域に記憶されている基本セグメントのパラメータと残差とから、基本セグメント毎に再構成スペクトルを求めて記憶するセグメント解析結果構成手段とから構成すると共に、基本セグメント毎に頂位位相の周波数依存性を解析するための頂位位相推移計算手段を構成し、この頂位位相推移計算手段により計算した各基本セグメント毎の頂位位相の周波数依存をセグメント解析結果構成手段による再構成スペクトルと共に、基本セグメントの代表時刻の属性として記憶することを特徴とする請求項１０に記載の２系列の時系列データの解析装置。
【請求項１２】
再構成スペクトル計算手段は、得られたパラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データの再構成スペクトルを求める構成とすることを特徴とする請求項１１に記載の２系列の時系列データの解析装置。
【請求項１３】
δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用することを特徴とする請求項１２に記載の２系列の時系列データの解析装置。
【請求項１４】
解析対象の２系列の時系列データを、夫々解析して、その一般化三角多項式パラメータを求める単一系列解析手順と、単一系列解析手順により求められた２系列の時系列データの一般化三角多項式パラメータから、それらの時系列データの間の相互相関を定量化する諸量を求める相互相関解析手順とから構成し、単一系列解析手順は、解析対象の時系列データを解析するシーケンスを決定する解析シーケンス決定手順と、解析シーケンスに従い、最大エントロピー法及び非線形最小自乗法によりデータを解析して一般化三角多項式のパラメータと残差を求めるデータ解析手順とから構成し、相互相関解析手順は、一般化三角多項式パラメータからパワースペクトル密度を求めるパワースペクトル密度計算手順と、クロススペクトル（コスペクトル、クオドスペクトル）とコヒーレンスを含む相互相関を定量化する諸量を求める相互相関諸量計算手順と、相互相関諸量計算手順により求められたコヒーレンスを理論値と比較して、それらの差が所定値以下の場合に、その差をなくすようにクロススペクトルを補正する整合性評価補正手順と、整合性評価補正手順において評価した差が所定値以上の場合に、単一系列解析手順に再解析の指令を出力し、所定値以下の場合には、相互相関諸量計算手順により計算された諸量と、パワースペクトル密度計算手順により計算されたパワースペクトル密度及び整合性評価補正手順により補正されたクロススペクトルを出力する分岐出力手順を構成したことを特徴とする２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１５】
パワースペクトル密度計算手順は、一般化三角多項式パラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データのスペクトルを求める構成とする２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１６】
δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用することを特徴とする請求項１５に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１７】
解析シーケンス決定手順は、解析対象の時系列データを、予め設定された最大データ点数を越えないようにセグメントに分割する分割処理過程と、各セグメントに対してデータ構造の均一性を評価する評価処理過程と、この評価に基づき、不均一のセグメントデータを均一の複数のセグメントに分割する分割処理過程と、分割された均一のセグメントの夫々を基本セグメントとして、それらに対して、記憶手段に記憶領域を設定する記憶領域設定処理過程を備えていることを特徴とする請求項１４に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１８】
セグメントは、解析対象の時系列データの全期間データに対応する最上位レイヤーのセグメントから、順次２分割されて構成される、下位のレイヤーの夫々に含まれる２の累乗のセグメントと、最下位のレイヤーにおいて、データ構造が変化するセグメントを、データ構造の変化している個所において分割して順次形成されるサブレイヤーの各セグメントとから構成されることを特徴とする請求項１７に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項１９】
データ構造の均一性は、下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移により評価することを特徴とする請求項１７に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【数３】
2010108402.tif 000027
但し、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【請求項２０】
解析シーケンス決定手順には、不等間隔の時系列データを等間隔化処理する等間隔化処理過程を設けたことを特徴とする請求項１７に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項２１】
データ解析手順は、解析対象の時系列データから取り出したセグメントデータに対して下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移を求める計算過程と、求めた推移から最大エントロピー法によるパワースペクトル密度の計算に用いるラグ値を決定するラグ値決定過程と、決定されたラグ値により上記パワースペクトル密度を求めるパワースペクトル密度計算過程と、求められたパワースペクトル密度から主要ピークの周波数と、そのパワーを抽出する抽出過程と、抽出された周波数とパワーとから非線形最小自乗法により当てはめ曲線のパラメータを求めるパラメータ導出過程と、求められた当てはめ曲線とセグメントデータとを比較して一致度を評価し、所定の一致度以上の場合には上記当てはめ曲線のパラメータをセグメントに対応する記憶手段に記憶して次の過程に移行し、所定の一致度に達しない場合には上記当てはめ曲線のパラメータを破棄して上記ラグ値を決定する過程に移行する評価過程とを備えており、ラグ値決定過程は、設定される解析の精度に応じて、より高精度の解析において、より大きいラグ値を選択するように構成し、抽出過程は、設定される解析の精度に応じて、より高精度の解析において、より多くのピークを主要ピークとして抽出するように構成すると共に、評価過程は、設定される解析の精度に応じて、より高精度の解析において、より高い一致度となるように構成し、相互相関解析手順からの再解析指令により、高精度の解析が設定されるように構成したことを特徴とする請求項１７に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【数４】
2010108402.tif 000028
但し、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【請求項２２】
データ解析手順におけるセグメントデータの解析は、最上位のレイヤーから最下位のレイヤーに向かって順次行うものとし、この際、上位のレイヤーのセグメントのセグメントデータは、セグメント当たりのデータ点数が予め設定された最大データ点数を越えないように、解析対象の時系列データのデータ点数をバンチング処理により低減してセグメントデータを構成する処理を行う構成とすると共に、各レイヤーの各セグメントに対応して一時記憶領域を構成し、上位側のセグメントに対して求められた多項式パラメータを下位の各レイヤーの各セグメントの一時記憶領域に記憶すると共に、上記多項式パラメータから構成される当てはめ曲線と、原時系列データの対応部分との差演算により求められる残差時系列データを次の下位レイヤーのセグメントデータとして設定する処理を行うことを特徴とする請求項２１に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項２３】
データ解析手順により求めた解析結果を評価し、条件に適合する場合に記憶手段に記憶する解析結果評価手順と、上記パラメータから得られるパワースペクトル密度と最大エントロピー法を用いて求める残差のパワースペクトル密度を合成して時系列データの再構成スペクトルを求める再構成スペクトル計算手順と、この計算手順を利用して、上記記憶手段から取り出した上記パラメータと残差により、解析対象の時系列データに対応する一般化三角多項式を再構成する解析結果構成手順とを備えたことを特徴とする請求項１４に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項２４】
解析結果構成手順は、全ての記憶領域に記憶されている基本セグメントのパラメータと残差とから全データ領域を一括して再構成スペクトルを求めて記憶する一括解析結果構成手順と、各記憶領域に記憶されている基本セグメントのパラメータと残差とから、基本セグメント毎に再構成スペクトルを求めて記憶するセグメント解析結果構成手順とから構成すると共に、基本セグメント毎に頂位位相の周波数依存性を解析するための頂位位相推移計算手順を構成し、この頂位位相推移計算手順により計算した各基本セグメント毎の頂位位相の周波数依存をセグメント解析結果構成手順による再構成スペクトルと共に、基本セグメントの代表時刻の属性として記憶することを特徴とする請求項２３に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項２５】
再構成スペクトル計算手順は、得られたパラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データの再構成スペクトルを求める構成とすることを特徴とする請求項２４に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【請求項２６】
δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用することを特徴とする請求項２５に記載の２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、２系列の時系列データの解析装置及び２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体に関するものである。
【背景技術】
【０００２】
本発明者は、先に、時系列データを高精度に解析可能な解析装置として、例えば特許文献１、特許文献２に示されるように、最大エントロピー法（MEM）を利用した時系列データの解析手法及びその解析手法を適用したコンピュータプログラムを実行する解析装置を提案している。尚、最大エントロピー法を利用した時系列データの解析手法は、例えば非特許文献１、非特許文献2にも記載されている。
【０００３】
上述した解析手法は、解析対象の時系列データを一般化三角多項式により表現して解析を行うもので、概ね、時系列データに対して最大エントロピー法を適用して、高精度なパワースペクトル密度（PSD）を求め、その顕著なピークから、一般化三角多項式を解く初期値としての、項数と、各三角項の周期値を求める過程と、これらの初期値により、線形化された非線形最小自乗法を用いて、各三角項のパラメータを求めると共に、残差時系列データを求め、これらを加算して時系列データを再構成する過程とから構成されている。
【０００４】
この解析手法及びそれを適用した解析装置では、何らの数値モデルも想定せず、時系列データのみを用いて高精度の解析を行え、しかもその解析は、通常のパーソナルコンピュータ程度の演算能力でも十分に実用的であるという特徴を有している。
【特許文献１】特開平５−２１６１９５号公報
【特許文献２】特開平６−１４９８６３号公報
【非特許文献１】三宅浩次監修  高橋延昭・神山昭男・大友詔雄編、「生物リズムの構造  −MemCalcによる生物時系列データの解析−」、第１版、富士書院、1992年1月15日、p.19-39
【非特許文献２】常磐野和男・大友詔雄・田中幸雄著、「最大エントロピー法による時系列解析  …ＭｅｍＣａｌｃの理論と実際…」、第１版、北海道大学図書刊行会、2002年6月25日、p.79-99
【発明の開示】
【発明が解決しようとする課題】
【０００５】
しかしながら以上の文献に記載されている解析装置は、一系列の時系列データを解析するのに用いられるものであり、本発明者は、これらの解析装置に用いられている解析手法の上述した特徴を更に発揮させるべく、時系列データのみを用い、時系列に関わるその他の情報を一切使用することなく、二つの時系列の相互相関にかかわる知見を得ることのできる解析装置を提供することを課題とした。
【０００６】
ところで上述した解析装置で一系列の時系列データの解析を行う場合、演算速度、演算精度、記憶容量等の点でコンピュータの能力に限界があるため、長大な時系列データを一括して処理するのが困難であり、この場合には、長大な時系列データを分割して処理する必要がある。一方、時系列データには、観測期間内のデータ全体に渡ってデータ構造が均一に保たれるものの他、観測期間内でデータ構造が時々刻々と変化するものがあり、後者の場合には、観測期間内のデータ全体を一括して処理するのは困難である。
【０００７】
そこで本発明では、解析対象の時系列データを、そのデータ点数やデータ構造に応じて適切にセグメント化して処理を行うことにより、解析する時系列データによらず、常に適切に高精度な解析を行えるようにすることも課題とした。
【０００８】
尚、本発明が対象とする時系列データとは、実在の時系列データ、即ち、ある限られた期間に観測されたある点数の有限長離散時系列のデータであり、理論上の存在である無限長のデータや連続データは対象としない。また本発明が対象とする時系列データは、上述したような通常の時系列データの他、例えば地層断面の距離に対する輝度変化等の、空間的に連続するデータも対象とするものである。
【課題を解決するための手段】
【０００９】
上述した課題を解決するために、本発明では、解析対象の２系列の時系列データを、夫々解析して、その一般化三角多項式パラメータを求める２系統の単一系列解析手段と、夫々の単一系列解析手段により求められた２系列の時系列データの一般化三角多項式パラメータから、２系列の時系列データの間の相互相関を定量化する諸量を求める相互相関解析手段とから構成し、単一系列解析手段は、解析対象の時系列データを解析するシーケンスを決定する解析シーケンス決定手段と、解析シーケンスに従い、最大エントロピー法及び非線形最小自乗法によりデータを解析して一般化三角多項式のパラメータと残差を求めるデータ解析手段とから構成し、相互相関解析手段は、一般化三角多項式パラメータからパワースペクトル密度を求めるパワースペクトル密度計算手段と、クロススペクトル（コスペクトル、クオドスペクトル）とコヒーレンスを含む相互相関を定量化する諸量を求める相互相関諸量計算手段と、相互相関諸量計算手段により求められたコヒーレンスを理論値と比較して、それらの差が所定値以下の場合に、その差をなくすようにクロススペクトルを補正する整合性評価補正手段と、整合性評価補正手段において評価した差が所定値以上の場合に、単一系列解析手段に再解析の指令を出力し、所定値以下の場合には、相互相関諸量計算手段により計算された諸量と、パワースペクトル密度計算手段により計算されたパワースペクトル密度及び整合性評価補正手段により補正されたクロススペクトルを出力する分岐出力手段を構成した２系列の時系列データの解析装置を提案する。
【００１０】
また本発明では、上記の解析装置において、パワースペクトル密度計算手段は一般化三角多項式パラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データのスペクトルを求める構成とする２系列の時系列データの解析装置を提案する。
【００１１】
また本発明では、上記の解析装置において、δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用する２系列の時系列データの解析装置を提案する。
【００１２】
また本発明では、上記の解析装置において、解析シーケンス決定手段は、解析対象の時系列データを、予め設定された最大データ点数を越えないようにセグメントに分割する分割処理過程と、各セグメントに対してデータ構造の均一性を評価する評価処理過程と、この評価に基づき、不均一のセグメントデータを均一の複数のセグメントに分割する分割処理過程と、分割された均一のセグメントの夫々を基本セグメントとして、それらに対して、記憶手段に記憶領域を設定する記憶領域設定処理過程を備えている２系列の時系列データの解析装置を提案する。
【００１３】
また本発明では、上記の解析装置において、セグメントは、解析対象の時系列データの全期間データに対応する最上位レイヤーのセグメントから、順次２分割されて構成される、下位のレイヤーの夫々に含まれる２の累乗のセグメントと、最下位のレイヤーにおいて、データ構造が変化するセグメントを、データ構造の変化している個所において分割して順次形成されるサブレイヤーの各セグメントとから構成される２系列の時系列データの解析装置を提案する。
【００１４】
また本発明では、上記の解析装置において、データ構造の均一性は、下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移により評価する２系列の時系列データの解析装置を提案する。
【数１】
2010108402.tif 000003
但し、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【００１５】
また本発明では、上記の解析装置において、解析シーケンス決定手段には、不等間隔の時系列データを等間隔化処理する等間隔化処理過程を設けた２系列の時系列データの解析装置を提案する。
【００１６】
また本発明では、上記の解析装置において、データ解析手段は、解析対象の時系列データから取り出したセグメントデータに対して下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移を求める計算過程と、求めた推移から最大エントロピー法によるパワースペクトル密度の計算に用いるラグ値を決定するラグ値決定過程と、決定されたラグ値により上記パワースペクトル密度を求めるパワースペクトル密度計算過程と、求められたパワースペクトル密度から主要ピークの周波数と、そのパワーを抽出する抽出過程と、抽出された周波数とパワーとから非線形最小自乗法により当てはめ曲線のパラメータを求めるパラメータ導出過程と、求められた当てはめ曲線とセグメントデータとを比較して一致度を評価し、所定の一致度以上の場合には上記当てはめ曲線のパラメータをセグメントに対応する記憶手段に記憶して次の過程に移行し、所定の一致度に達しない場合には上記当てはめ曲線のパラメータを破棄して上記ラグ値を決定する過程に移行する評価過程とを備えており、ラグ値決定過程は、設定される解析の精度に応じて、より高精度の解析において、より大きいラグ値を選択するように構成し、抽出過程は、設定される解析の精度に応じて、より高精度の解析において、より多くのピークを主要ピークとして抽出するように構成すると共に、評価過程は、設定される解析の精度に応じて、より高精度の解析において、より高い一致度となるように構成し、相互相関解析手段からの再解析指令により、高精度の解析が設定されるように構成した２系列の時系列データの解析装置を提案する。
【００１７】
また本発明では、上記の解析装置において、データ解析手段におけるセグメントデータの解析は、最上位のレイヤーから最下位のレイヤーに向かって順次行うものとし、この際、上位のレイヤーのセグメントのセグメントデータは、セグメント当たりのデータ点数が予め設定された最大データ点数を越えないように、解析対象の時系列データのデータ点数をバンチング処理により低減してセグメントデータを構成する処理を行う構成とすると共に、各レイヤーの各セグメントに対応して一時記憶領域を構成し、上位側のセグメントに対して求められた多項式パラメータを下位の各レイヤーの各セグメントの一時記憶領域に記憶すると共に、上記多項式パラメータから構成される当てはめ曲線と、原時系列データの対応部分との差演算により求められる残差時系列データを次の下位レイヤーのセグメントデータとして設定する処理を行う２系列の時系列データの解析装置を提案する。
【００１８】
また本発明では、上記の解析装置において、データ解析手段により求めた解析結果を評価し、条件に適合する場合に記憶手段に記憶する解析結果評価手段と、上記パラメータから得られるパワースペクトル密度と最大エントロピー法を用いて求める残差のパワースペクトル密度を合成して時系列データの再構成スペクトルを求める再構成スペクトル計算手段と、この計算手段を利用して、上記記憶手段から取り出した上記パラメータと残差により、解析対象の時系列データに対応する一般化三角多項式を再構成する解析結果構成手段とを備えた２系列の時系列データの解析装置を提案する。
【００１９】
また本発明では、上記の解析装置において、解析結果構成手段は、全ての記憶領域に記憶されている基本セグメントのパラメータと残差とから全データ領域を一括して再構成スペクトルを求めて記憶する一括解析結果構成手段と、各記憶領域に記憶されている基本セグメントのパラメータと残差とから、基本セグメント毎に再構成スペクトルを求めて記憶するセグメント解析結果構成手段とから構成すると共に、基本セグメント毎に頂位位相の周波数依存性を解析するための頂位位相推移計算手段を構成し、この頂位位相推移計算手段により計算した各基本セグメント毎の頂位位相の周波数依存をセグメント解析結果構成手段による再構成スペクトルと共に、基本セグメントの代表時刻の属性として記憶する２系列の時系列データの解析装置を提案する。
【００２０】
また本発明では、上記の解析装置において、再構成スペクトル計算手段は、得られたパラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データの再構成スペクトルを求める構成とする２系列の時系列データの解析装置を提案する。
【００２１】
また本発明では、上記の構成において、δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用する２系列の時系列データの解析装置を提案する。
【００２２】
また本発明では、解析対象の２系列の時系列データを、夫々解析して、その一般化三角多項式パラメータを求める単一系列解析手順と、単一系列解析手順により求められた２系列の時系列データの一般化三角多項式パラメータから、それらの時系列データの間の相互相関を定量化する諸量を求める相互相関解析手順とから構成し、単一系列解析手順は、解析対象の時系列データを解析するシーケンスを決定する解析シーケンス決定手順と、解析シーケンスに従い、最大エントロピー法及び非線形最小自乗法によりデータを解析して一般化三角多項式のパラメータと残差を求めるデータ解析手順とから構成し、相互相関解析手順は、一般化三角多項式パラメータからパワースペクトル密度を求めるパワースペクトル密度計算手順と、クロススペクトル（コスペクトル、クオドスペクトル）とコヒーレンスを含む相互相関を定量化する諸量を求める相互相関諸量計算手順と、相互相関諸量計算手順により求められたコヒーレンスを理論値と比較して、それらの差が所定値以下の場合に、その差をなくすようにクロススペクトルを補正する整合性評価補正手順と、整合性評価補正手順において評価した差が所定値以上の場合に、単一系列解析手順に再解析の指令を出力し、所定値以下の場合には、相互相関諸量計算手順により計算された諸量と、パワースペクトル密度計算手順により計算されたパワースペクトル密度及び整合性評価補正手順により補正されたクロススペクトルを出力する分岐出力手順を構成した２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００２３】
また本発明では、上記の記録媒体において、パワースペクトル密度計算手順は、一般化三角多項式パラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データのスペクトルを求める構成とする２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００２４】
また本発明では、上記の記録媒体において、δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用する２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００２５】
また本発明では、上記の記録媒体において、解析シーケンス決定手順は、解析対象の時系列データを、予め設定された最大データ点数を越えないようにセグメントに分割する分割処理過程と、各セグメントに対してデータ構造の均一性を評価する評価処理過程と、この評価に基づき、不均一のセグメントデータを均一の複数のセグメントに分割する分割処理過程と、分割された均一のセグメントの夫々を基本セグメントとして、それらに対して、記憶手段に記憶領域を設定する記憶領域設定処理過程を備えている２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００２６】
また本発明では、上記の記録媒体において、セグメントは、解析対象の時系列データの全期間データに対応する最上位レイヤーのセグメントから、順次２分割されて構成される、下位のレイヤーの夫々に含まれる２の累乗のセグメントと、最下位のレイヤーにおいて、データ構造が変化するセグメントを、データ構造の変化している個所において分割して順次形成されるサブレイヤーの各セグメントとから構成される２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００２７】
また本発明では、上記の記録媒体において、データ構造の均一性は、下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移により評価する２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【数２】
2010108402.tif 000004
但し、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【００２８】
また本発明では、上記の記録媒体において、解析シーケンス決定手順には、不等間隔の時系列データを等間隔化処理する等間隔化処理過程を設けた２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００２９】
また本発明では、上記の記録媒体において、データ解析手順は、解析対象の時系列データから取り出したセグメントデータに対して下式による最大ラグ値までの最大エントロピー法のスペクトル密度の計算におけるバーグの係数Pmの、ラグ値に対する推移を求める計算過程と、求めた推移から最大エントロピー法によるパワースペクトル密度の計算に用いるラグ値を決定するラグ値決定過程と、決定されたラグ値により上記パワースペクトル密度を求めるパワースペクトル密度計算過程と、求められたパワースペクトル密度から主要ピークの周波数と、そのパワーを抽出する抽出過程と、抽出された周波数とパワーとから非線形最小自乗法により当てはめ曲線のパラメータを求めるパラメータ導出過程と、求められた当てはめ曲線とセグメントデータとを比較して一致度を評価し、所定の一致度以上の場合には上記当てはめ曲線のパラメータをセグメントに対応する記憶手段に記憶して次の過程に移行し、所定の一致度に達しない場合には上記当てはめ曲線のパラメータを破棄して上記ラグ値を決定する過程に移行する評価過程とを備えており、ラグ値決定過程は、設定される解析の精度に応じて、より高精度の解析において、より大きいラグ値を選択するように構成し、抽出過程は、設定される解析の精度に応じて、より高精度の解析において、より多くのピークを主要ピークとして抽出するように構成すると共に、評価過程は、設定される解析の精度に応じて、より高精度の解析において、より高い一致度となるように構成し、相互相関解析手順からの再解析指令により、高精度の解析が設定されるように構成した２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【数３】
2010108402.tif 000005
但し、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【００３０】
また本発明では、上記の記録媒体において、データ解析手順におけるセグメントデータの解析は、最上位のレイヤーから最下位のレイヤーに向かって順次行うものとし、この際、上位のレイヤーのセグメントのセグメントデータは、セグメント当たりのデータ点数が予め設定された最大データ点数を越えないように、解析対象の時系列データのデータ点数をバンチング処理により低減してセグメントデータを構成する処理を行う構成とすると共に、各レイヤーの各セグメントに対応して一時記憶領域を構成し、上位側のセグメントに対して求められた多項式パラメータを下位の各レイヤーの各セグメントの一時記憶領域に記憶すると共に、上記多項式パラメータから構成される当てはめ曲線と、原時系列データの対応部分との差演算により求められる残差時系列データを次の下位レイヤーのセグメントデータとして設定する処理を行う２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００３１】
また本発明では、上記の記録媒体において、データ解析手順により求めた解析結果を評価し、条件に適合する場合に記憶手段に記憶する解析結果評価手順と、上記パラメータから得られるパワースペクトル密度と最大エントロピー法を用いて求める残差のパワースペクトル密度を合成して時系列データの再構成スペクトルを求める再構成スペクトル計算手順と、この計算手順を利用して、上記記憶手段から取り出した上記パラメータと残差により、解析対象の時系列データに対応する一般化三角多項式を再構成する解析結果構成手順とを備えた２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００３２】
また本発明では、上記の記録媒体において、解析結果構成手順は、全ての記憶領域に記憶されている基本セグメントのパラメータと残差とから全データ領域を一括して再構成スペクトルを求めて記憶する一括解析結果構成手順と、各記憶領域に記憶されている基本セグメントのパラメータと残差とから、基本セグメント毎に再構成スペクトルを求めて記憶するセグメント解析結果構成手順とから構成すると共に、基本セグメント毎に頂位位相の周波数依存性を解析するための頂位位相推移計算手順を構成し、この頂位位相推移計算手順により計算した各基本セグメント毎の頂位位相の周波数依存をセグメント解析結果構成手順による再構成スペクトルと共に、基本セグメントの代表時刻の属性として記憶する２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００３３】
また本発明では、上記の記録媒体において、再構成スペクトル計算手順は、得られたパラメータからパワースペクトル密度を求める計算において、一般化三角多項式の夫々の項の夫々につき、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成し、その重畳として時系列データの再構成スペクトルを求める構成とする２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【００３４】
また本発明では、上記の記録媒体において、δ函数の原型として、余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルを利用する２系列の時系列データの解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提案する。
【発明の効果】
【００３５】
請求項１、１４の発明によれば、２系列の時系列データの夫々を単一系列解析手段において、まず最大エントロピー法を適用して高精度なパワースペクトル密度を求めると共に、その顕著なピークから、一般化三角多項式を解く初期値としての、項数と、各三角項の周期値を求め、次いで、これらの初期値により非線形最小自乗法を用いて、各三角項のパラメータを求めると共に、残差時系列データを求め、これらを加算して時系列データを再構成するという解析を行うことにより、何らの数値モデルも想定せず、夫々の時系列データのみを用いて高精度の解析を行え、これらの高精度の解析によって得られた２系列の時系列データの一般化三角多項式パラメータにより、相互相関解析手段において、２系列の時系列データの間の相互相関を定量化する諸量を求めることができる。
【００３６】
特に上述した発明によれば、相互相関解析手段においては、一般化三角多項式パラメータにより求められたコヒーレンスを理論値と比較して整合性を評価し、所定以上の整合性が得られない場合には、単一系列解析手段に対して再解析指令を発して、より高精度の解析を行わせるという、フィードバックと繰り返し動作を有するため、２系列の時系列データの間の相互相関を定量化する諸量を非常に高精度に求めることができる。
【００３７】
請求項２、３、１２、１３、１５、１６、２５及び２６の発明によれば、理論形式では有限離散値を扱うコンピュータで構成することができない線スペクトル列となる時系列データのスペクトルを、ピーク周波数から離れるに従って急減する連続関数、例えば余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルをδ函数の原型として利用して再構成することにより、コンピュータを利用して、理論と整合する時系列データの解析を行うことができる。また、これにより、２つの時系列データのクロススペクトルの絶対値と等しくなる、２つの時系列データの積を構成することも可能となる。
【００３８】
請求項４、１７の発明によれば、解析シーケンス決定手段又は手順により、解析対象の時系列データを、その全体のデータ点数に応じて分割して解析装置の演算能力に適合させると共に、データ構造が変化している部分においてセグメントに分割することにより、セグメント中にデータ構造が変化している部分を無くすことにより、これらの均一なセグメントにつき、データ解析手段又は手順により高精度に解析を行うことができる。
【００３９】
請求項５、１８の発明によれば、セグメントは、解析シーケンス決定手段又は手順により、解析対象の時系列データの全期間データに対応する、最上位レイヤーのセグメントから、順次２分割されて構成される、下位のレイヤーの夫々に含まれる２の累乗のセグメントと、最下位のレイヤーにおいて、データ構造が変化するセグメントを、データ構造の変化している個所において分割して順次形成されるサブレイヤーの各セグメントとから構成するため、解析装置の演算能力と時系列データのデータ構造の両方に対応して、合理的に時系列データのセグメント化を行うことができ、自動化も容易である。
【００４０】
請求項６、１９の発明によれば、データ構造の均一性を容易に評価することができる。
【００４１】
請求項７、２０の発明によれば、原時系列データが不等間隔の場合にも、等間隔化することにより、最大エントロピー法のスペクトル密度の計算を行うことができる。一方、非線形最小自乗法により当てはめ曲線のパラメータを求めるパラメータの導出においては、必要に応じて、等間隔化した時系列データか、不等間隔の原時系列データをそのまま使用することができる。
【００４２】
請求項８、２１の発明によれば、データ解析手段又は手順において、求められた当てはめ曲線とセグメントデータとを比較して一致度を評価し、所定の一致度以上の場合には上記当てはめ曲線のパラメータをセグメントに対応する記憶手段に記憶して次の過程に移行し、所定の一致度に達しない場合には上記当てはめ曲線のパラメータを破棄して上記ラグ値を決定する過程に移行する評価過程とを備えており、パラメータの導出と、その評価を繰り返すことにより、単一系列解析手段自体において、より高精度の当てはめ曲線を求めることができる。
【００４３】
更にこれらの発明によれば、ラグ値決定過程は、設定される解析の精度に応じて、より高精度の解析において、より大きいラグ値を選択するように構成され、また抽出過程は、設定される解析の精度に応じて、より高精度の解析において、より多くのピークを主要ピークとして抽出するように構成され、更に評価過程は、設定される解析の精度に応じて、より高精度の解析において、より高い一致度となるように構成されているので、相互相関解析手段からの再解析指令毎に、より高精度の解析が行われる。勿論、解析の精度は、相互相関解析手段からの再解析指令のみでなく、操作者が設定可能とすることもできる。
【００４４】
請求項９、２２の発明によれば、解析対象の時系列データの全体のデータ点数を一括して処理するのには、解析装置の演算能力が足りない場合であっても、バンチング処理によりデータ点数を低減することにより、最上位のレイヤーのセグメントデータ、即ち、解析対象の時系列データの全体の一括した解析が可能であり、従って時系列データの全体に渡る長周期の解析が可能であると共に、上位のレイヤーの解析により得られた当てはめ曲線と、原時系列データの対応部分との差演算により求められる残差時系列データを次の下位レイヤーのセグメントデータとして設定して、順次解析を行うので、最下位又はサブレイヤーのセグメントデータにより、バンチング処理によるデータ点数の低減がない、そのままの時系列データに対しての解析を行うことができる。
【００４５】
請求項１０、２３の発明によれば、解析結果評価手段又は手順において、得られた当てはめ曲線が、解析対象の時系列データの対応範囲と、所定の一致度に達しない場合にはデータを破棄して、セグメントデータを次のレイヤーのセグメントに移行して、処理に供されるので、上位のレイヤーにおいて一致度の高い当てはめ曲線が求まらない場合でも、不都合が生じず、下位のレイヤーにおいての処理により、一致度の高い当てはめ曲線を得ることができる。
【００４６】
請求項１１、２４の発明によれば、単一系列解析手段において、解析結果構成手段又は手順により、解析対象の時系列データを、一括した解析結果と、各セグメント毎の解析結果として出力することができる。
【発明を実施するための最良の形態】
【００４７】
次に本発明に係る２系列の時系列データの解析装置の最良の形態を、添付図面を参照して説明する。
まず、図１に示すように本発明に係る統合構造解析手段１は、解析対象の２系列の時系列データを、夫々解析して、その一般化三角多項式パラメータを求める２系統の単一系列解析手段２ａ，２ｂと、夫々の単一系列解析手段２ａ，２ｂにより求められた２系列の時系列データの一般化三角多項式パラメータから、２系列の時系列データの間の相互相関を定量化する諸量を求める相互相関解析手段３とから構成されている。夫々の単一系列解析手段２ａ，２ｂにより求められた２系列の時系列データの一般化三角多項式パラメータは、図中実線矢印に示すように相互相関解析手段３に入力され、一方、図中破線矢印で示すように相互相関解析手段３から夫々の単一系列解析手段２ａ，２ｂに必要に応じて再解析指令が発せられる。尚、図１において、符号４ａ，４ｂは２系列の夫々の時系列データ、符号５ａ，５ｂ，５ｃは本発明の解析装置から出力される諸量を示すものである。
【００４８】
単一系列解析手段２ａ，２ｂは、上述した解析手法を適用するため、解析対象の時系列データ（以降、必要に応じて原時系列データと称する）を一般化三角多項式により表現して解析を行うように構成されており、概ね、原時系列データに対して最大エントロピー法を適用して、高精度なパワースペクトル密度（PSD）を求め、その顕著なピークから、一般化三角多項式を解く初期値としての、項数と、各三角項の周期値を求める過程と、これらの初期値により、非線形最小自乗法を用いて、各三角項のパラメータを求めると共に、残差時系列データを求め、これらを加算して時系列データを再構成する過程を実行するように構成されている。
【００４９】
即ち、本発明に係る単一系列解析手段２ａ，２ｂにおいては、原時系列データは、次式の一般化三角多項式で表現される。
【数４】
2010108402.tif 000006
ここで、x(t)は原時系列データ、即ち、有限長離散時系列のデータを示し、Lは項数、ajは振幅、Tjは周期、φjは頂位位相、εx(t)は雑音成分を示しており、上述した一般化三角多項式のパラメータは、これらの項数(モード数)、各項の振幅aj、周期Tj及び頂位位相φjから成る。尚、頂位位相とは、そのモードの値が極大を示す時刻で、そのモードの周期ごとに現れる時刻のうちの一つを示している。
【００５０】
図２は単一系列解析手段２ａ，２ｂの実施の形態を模式的に表した全体構成の説明図であり、図中矢印により処理の流れを模式的に表している。
図２に示すように、この実施の形態に係る単一系列解析手段２ａ，２ｂは、原時系列データを解析するシーケンスを決定する解析シーケンス決定手段６と、この解析シーケンス決定手段６により決定された解析シーケンスに従い、最大エントロピー法及び非線形最小自乗法によりデータを解析して一般化三角多項式のパラメータと残差を求めるデータ解析手段７と、データ解析手段７により求めた解析結果を評価し、条件に適合する場合に、予め設定された記憶手段に記憶する解析結果評価手段８と、これらの解析シーケンス決定手段６と、データ解析手段７と、解析結果評価手段８を統合的に制御して、必要な繰り返し演算等を行わせるための制御手段９と、最大エントロピー法を用いて一般化三角多項式の再構成スペクトルを求める再構成スペクトル計算手段１０と、頂位位相の推移を計算する頂位位相推移計算手段１１と、出力するデータの準備等に係る計算を行う、その他計算手段１２と、これらの各計算手段１０〜１２を用いて、出力するデータ解析手段７の解析結果を構成する解析結果構成手段１３ａ，１３ｂとから構成されている。
【００５１】
後述するように、解析シーケンス決定手段６は、解析シーケンスの一つとして、原時系列データの全期間データのデータ点数やデータ構造等に応じて、それを複数レイヤーに構成した、データ構造の均一な複数のセグメントに分割すると共に、分割された均一なセグメントの夫々を基本セグメントとして、それらに対応して記憶手段に記憶領域を設定する過程を有し、この分割されたセグメントに取り込まれた原時系列データについてデータ解析手段７において解析されるものである。そしてデータ解析手段７において各セグメントの原時系列データについて得られた解析結果は、対象とする時系列データ全体を一括する一括解析結果と、セグメント解析結果として構成されるものであり、このために上記解析結果構成手段１３ａ，１３ｂは、前者に対応する一括解析結果構成手段１３ａと、後者に対応するセグメント解析結果構成手段１３ｂとから構成されている。
【００５２】
図３はこの実施の形態の単一系列解析手段２ａ，２ｂから出力される解析結果を模式的に示すものであり、原時系列データの全期間に渡ってデータ構造が保たれる一括の解析結果１４と、原時系列データの期間内でデータ構造が時々刻々と変動し、原時系列データの全期間を網羅する複数のセグメント毎に得られる解析結果１５とから構成される。従って、前者の解析結果１４は時間依存がなく、後者の解析結果１５は時間依存を有するものである。
【００５３】
解析結果１４において、一般化三角多項式パラメータａは、上述した一般化三角多項式の項数、振幅、周期及び頂位位相であり、これらのパラメータａから、周期的に変動する曲線、即ち、多項式パラメータによる最適あてはめ曲線ｂと、そのスペクトル、即ち、多項式パラメータによる再構成スペクトルｃが構成され、これに残差と残差のスペクトルｄによる残差成分を加算することにより、原時系列データとそのスペクトルが再構成される。
【００５４】
一方、解析結果１５においては、一般化三角多項式パラメータｅと、多項式パラメータによる再構成スペクトルｆと、多項式パラメータによる頂位位相ｇとが夫々各セグメント毎に構成され、従って再構成スペクトルｆと頂位位相ｇの時間変動が構成される。
【００５５】
次に、単一系列解析手段２ａ，２ｂの、上述した各構成要素の動作を説明する。まず図４、図５は解析シーケンス決定手段６における処理の流れを示す流れ図であり、図６は解析シーケンス決定手段６により決定されたレイヤーとセグメントの構成例を模式的に示すものである。
【００５６】
まず、単一系列解析手段２ａ，２ｂにステップＳ１において原時系列データが入力されると、解析シーケンス決定手段６は、まずステップＳ２において、一つのセグメントが配置されるレイヤーを設定する。この一つのセグメントが配置されるレイヤーは必ず設定されるものであるので、その設定時点はステップＳ２の時点でなくとも良い。
【００５７】
次にステップＳ３では、入力された原時系列データの全期間のデータ点数が、予め設定されている最大データ点数を越えているか、否かを判定する。
この最大データ点数は、後述する基本セグメントに取り込んで解析が行われる最大データ点数であるが、一般に、各セグメントに取り込んで同時に解析するデータ点数が小さすぎるとデータの構造の抽出が困難となり、逆に大きすぎると過剰なデータ構造が抽出されてしまうため、最適なデータ点数とするのが高精度の解析結果を得るためには好ましい。一方、解析に最適なデータ点数は、取り込んだデータの周期性が強いか否かによっても異なり、また単一系列解析手段２ａ，２ｂを構成する演算手段の演算速度や演算精度等によってもデータ点数が制限される。
【００５８】
このようなことから、最大データ点数としては、それまでの多数の時系列データの解析により蓄積された知見等から得られた値が設定される。尚、この最大データ点数は、後述する各評価手段による評価の結果をフィードバックして変更可能とし、変更された最大データ点数に基づいてステップＳ３及びそれ以降の処理を行うようにすることもできる。また最大データ点数は、操作者により設定可能とすることもできる。
【００５９】
ステップＳ３において、入力された時系列データの全体のデータ点数が、予め設定されている最大データ点数を越えていないと判定された場合には、一つのセグメントが配置された一つのレイヤーが設定された状態で、結合子Ａを経て図５のステップＳ６に移行する。
【００６０】
一方、ステップＳ３において、入力された時系列データの全体のデータ点数が、予め設定されている最大データ点数を越えていると判定された場合には、ステップＳ４に移行し、そのステップＳ４において、それまでに設定された上位のレイヤーの下位にレイヤーを追加し、その下位のレイヤーには、上位のレイヤーのセグメント数の２倍のセグメントを配置して、ステップＳ５に移行する。
【００６１】
ステップＳ５では、入力された原時系列データの全期間データが、ステップＳ４において配置された全セグメントに等分されて取り込まれた場合の各セグメント毎のデータ点数が上記最大データ点数を越えているか、否かを判定する。
【００６２】
ステップＳ５において、各セグメントのデータ点数が、上記最大データ点数を越えていると判定された場合には、ステップＳ４に移行して、更に下位のレイヤーが追加され、一方、各セグメントのデータ点数が、上記最大データ点数を越えていないと判定された場合には、結合子Ａを経て図５のステップＳ６に移行する。
【００６３】
以上の各ステップにおける処理を経て、原時系列データを、予め設定された最大データ点数を越えないようにセグメントに取り込むための、レイヤーとセグメントの構成が得られ、構成された単数又は順次配列された複数のセグメントが、原時系列データの全期間を網羅する。従って以上のステップの処理は、上述した分割処理過程に相当する。
【００６４】
次いでステップＳ６では、最下位のレイヤーのセグメントに原時系列データの対応部分が取り込まれ、セグメントデータが構成される。
【００６５】
次いでステップＳ７では、セグメントデータに対して、次式の最大エントロピー法により最小ラグ値から最大ラグ値までバーグの係数Pmを求め、ラグ値に対するPmの推移をステップＳ８において評価する。
【数５】
2010108402.tif 000007
ここで、Bm(f)はラグ値mにおけるスペクトル密度、Δtはサンプリング間隔、Pmとγm(k)はバーグの係数である。
【００６６】
後に例示するように、係数Pmは、ラグ値が増大するに従って単調減少する量であり、ホワイトノイズ的な時系列データでは、この減少は極めてゆるやかであるのに対して、三角関数波に一定量の雑音を重畳した時系列データや、その他の通常観測されるように周期性を有する成分を含む時系列データでは、係数Pmの減少はやや大きくなる。また、周期性を有する成分の他に雑音を含まないか、またはごくわずかに含む三角関数波として記述されるような極めて周期性の高い時系列データでは、係数Pmの減少は極めて急激におこり、最小ラグ値に対する係数Pmの値と、最大ラグ値（セグメントデータのデータ点数−１）に対する係数Pmの比は7〜10桁以上に達する。
【００６７】
そして、セグメントデータの全期間に渡ってデータ構造が変化せず、均一な時系列データでは、係数Pmは最小のラグ値から増大する際に、最初急激に減少した後に、プラトー領域に至って最大ラグ値まで緩やかに漸減するのに対して、セグメントデータの期間内においてデータ構造が変化する不均一な時系列データでは、プラトー領域において、係数Pmの急激な変化が生じる。
【００６８】
そこでステップＳ８においては、セグメントデータの評価として、ａ：最小ラグ値に対する係数Pmの値と、最大ラグ値に対する係数Pmの比を算出し、その値を設定値と比較して、セグメントデータの周期性の強さを評価することと、ｂ：プラトー領域における急激な値の変化を検出して、データ構造の変化を検出することが行われる。
【００６９】
ステップＳ８においてデータ構造が均一であると判定された場合には、ステップＳ９に移行して、データ構造を評価したセグメントを基本セグメントとして設定し、次いでステップＳ１０に移行して、全てのセグメントに対してセグメントデータのデータ構造の評価が完了したか否かを判定し、完了していない場合にはステップＳ６に移行して、次のセグメントに対しての上述した処理が行われる。
【００７０】
一方、ステップＳ８においてデータ構造が不均一であると判定された場合には、ステップＳ１１に移行し、ステップＳ１１においては、ステップＳ８において評価されたセグメントのレイヤーの下位にサブレイヤーを構成し、このサブレイヤーに、データ構造の変化している個所において２分割したセグメントを配置する。
【００７１】
次いでステップＳ６に移行して、ステップＳ１１において構成されたサブレイヤーに配置された２つのセグメントに対して、上述と同様な処理が行われる。即ち、ステップＳ８において、サブレイヤーのセグメントのデータ構造が均一であると判定された場合には、このセグメントもステップＳ１０において基本セグメントとして設定され、一方、データ構造が不均一であると判定された場合には、サブレイヤーの下位に更にサブレイヤーが構成され、このサブレイヤーに、データ構造の変化している個所において２分割したセグメントが配置され、更にステップＳ６に移行して上述と同様な処理が行われる。
【００７２】
こうしてステップＳ１０において、サブレイヤーのセグメントを含む全てのセグメントに対してセグメントデータのデータ構造の評価が完了したか否かを判定し、完了と判定された場合には、ステップＳ１２に移行し、ステップＳ９において基本セグメントとして設定された全てのセグメントに対応して、記憶手段に記憶領域を設定し、初期化がされて、解析シーケンスが完成する。
【００７３】
図６は、ある解析対象の時系列データに対して、解析シーケンス決定手段６が決定したレイヤーとセグメントの構成を模式的に示すものであり、レイヤーは、時系列データの全期間データに対応する一つのセグメントが配置されるレイヤー０の下位に、レイヤーＮ（Ｎ：自然数）まで、Ｎ＋１層のレイヤーが構成されており、また、レイヤーＮのセグメント３に対応してサブレイヤー１が構成されている。こうしてレイヤーＮの、セグメント３を除く全てのセグメントが基本セグメントとして設定され、またレイヤーＮのセグメント３が分割された、サブレイヤー１の２つのセグメントが基本セグメントとして設定されている。そしてこれらの基本セグメントの全てに対応して記憶手段に記憶領域が設定されている。
【００７４】
データ解析手段７は、解析シーケンス決定手段６により上述したように決定されたレイヤーとセグメントの構成に基づき、セグメントに取り込まれた原時系列データ、即ち、セグメントデータの解析を行う。
【００７５】
データ解析手段７によるセグメントデータの解析は、最上位のレイヤー０から最下位のレイヤーＮ及びサブレイヤー１に向かって順次同様な処理の流れで行われるが、この際、上位のレイヤーのセグメントのセグメントデータは、上述したようにセグメント当たりのデータ点数が予め設定された最大データ点数を越えないように、原時系列データのデータ点数をバンチング処理により低減してセグメントデータを構成して処理が行われる。例えばレイヤー０のセグメントには、バンチング処理により、原時系列データの全期間のデータ点数を１／２Ｎに低減したデータをセグメントデータとして取り込み、レイヤー１のセグメントには、原時系列データの全期間のデータ点数を１／２Ｎ−１に低減したデータをセグメントとして取り込む。
【００７６】
更に、上位のレイヤーにおける各セグメントの解析結果中、上位側のセグメントに対して求められたパラメータは、その上位のセグメントに対応する全ての基本セグメントの記憶領域に同様に順次記憶されて行く。一方、上記パラメータから構成される当てはめ曲線と、原時系列データの対応部分との差演算により求められる残差時系列データは次の下位レイヤーのセグメントデータとして設定して解析が行われる。
【００７７】
そこで、データ解析手段７及び解析結果評価手段８における各セグメントの解析処理の流れを図７、図８の流れ図につき説明する。
まずステップＳ１０１では、セグメントに原時系列データの対応部分が取り込まれて、セグメントデータが構成される。この際、レイヤー０〜Ｎ−１のセグメントの場合には、上述したとおり、バンチング処理によるデータ点数の低減が行われる。
【００７８】
次いでステップＳ１０２では、ステップＳ１０１において取り込まれたセグメントデータに対して、上述のステップＳ８と同様に最大エントロピー法により最小ラグ値から最大ラグ値までバーグの係数Pmを求め、夫々のラグ値に対するPmの推移を算出する。
【００７９】
次いでステップＳ１０３では、上述のステップＳ８と同様に、ラグ値に対するバーグの係数Pmの推移を評価し、それに基づいて最大エントロピー法によるパワースペクトル密度の計算に用いるＭ個のラグ値を決定する。このラグ値は単数又は複数で良い。
【００８０】
このラグ値の決定は、例えば、ａ：周期性の強いデータ、即ち最小ラグ値に対する係数Pmの値と最大ラグ値に対する係数Pmの比が大きい場合には、ラグ値を50％以内とし、データ点数が多い場合には、より小さなラグ値とする、ｂ：周期性の弱い場合には、ラグ値は50%以上で、時系列の構造を効率的に取り出すために75%までとし、データ点数が少ない場合には、より大きなラグ値とする、等のラグ値決定のための条件や、その他、数多くの時系列データにおける解析から得られた知見を知識ベースとして記憶しておき、上述したバーグの係数Pmの評価結果から知識ベースを用いて、単数又は複数のラグ値を自動的に決定する構成とすることができる。この他、単一系列解析手段２ａ，２ｂのインタラクティブな操作により人がラグ値の決定を行うようにすることができる。更に本発明では、このラグ値決定過程においては、設定される解析の精度に応じて、より高精度の解析において、より大きいラグ値を選択するように知識ベース等により構成されており、高精度の解析は、上述した操作者の操作に加えて、相互相関解析手段３からの再解析指令により設定されるように構成している。
【００８１】
なお、最大エントロピー法において時系列が雑音を含まず線スペクトル成分のみからなる特殊な場合を除き、ラグ値をデータ点数−１点まで任意に設定し得て、その際になんら制約の存在しないことは、前掲の非特許文献２において明かにされた重要な知見である。
【００８２】
次いでステップＳ１０３で決定されたＭ個のラグ値の夫々につき、ステップＳ１０４において、上述した文献等に開示される最大エントロピー法によるパワースペクトル密度（MEM-PSD）の計算を行う。
【００８３】
次いでステップＳ１０５では、ステップＳ１０４において算出したＭ個のパワースペクトル密度の夫々から、適数の主要なピークと、そのパワーを抽出する。この主要なピーク及びそのパワーの抽出動作においては、一般に、最初に抽出された主要はピークに近接しているピークを始めとして、より多くのピークを主要ピークとして抽出すると、原時系列データの再現性がより向上する。このため、この抽出過程においては、設定される解析の精度に応じて、より高精度の解析において、より多くのピークを主要ピークとして抽出するように知識ベース等により構成されており、高精度の解析は、操作者の操作に加えて、相互相関解析手段３からの再解析指令により設定されるように構成している。
【００８４】
次いでステップＳ１０６では、ステップＳ１０５で抽出した主要ピークの周波数の逆数を周期の初期値として、非線形最小自乗法により、複数の当てはめ曲線とそのパラメータが算出される。
【００８５】
次いでステップＳ１０７は、ステップＳ１０６で求められた当てはめ曲線のパラメータを、初期値と比較して一致度を評価する。この一致度の評価は、例えば、当てはめ曲線の各周期のモードのパワー（振幅の自乗の１／２）を、上記主要ピークのパワーで除した値により評価する。
【００８６】
このようにしてステップＳ１０７において、最も一致度の高い当てはめ曲線を最適当てはめ曲線として選択する。
【００８７】
次いでステップＳ１０８では、ステップＳ１０７において選択された最適当てはめ曲線のパラメータの一致度が、予め設定された一致度以上であるか否かを判定し、設定された一致度以上でない場合には、ステップＳ１１１においてラグ値の再設定指示を行ってステップＳ１０３に移行し、また予め設定された一致度以上である場合には、ステップＳ１０９に移行して、ステップＳ１０７において選択された最適当てはめ曲線を、解析結果評価手段８により評価する。
【００８８】
本発明では、この最適当てはめ曲線のパラメータの一致度の判定において、設定される解析の精度に応じて、より高精度の解析において、より高い一致度となるように構成しており、高精度の解析は、相互相関解析手段３からの再解析指令に加えて、操作者の操作によっても設定されるように構成している。
【００８９】
図７においてステップＳ１０９として示す解析結果評価手段８による評価は、図８の流れ図に示す流れで行われる。
まずステップＳ１１２では、ステップＳ１０７において選出されたパラメータにより、最適当てはめ曲線を再構成し、そのデータ期間に対応する原時系列データとの差の演算を行い、残差曲線を求める。また当てはめ曲線の再現性を評価する諸量、ここでは残差の標準偏差と残差パワーの時間依存が求められる。
【００９０】
次いでステップＳ１１３においては、ステップＳ１１２において求めた当てはめ曲線の再現性を評価する諸量に基づいて、最適当てはめ曲線と、原時系列データの対応部分が十分に一致しているか否かを判定する。
【００９１】
ステップＳ１１３において十分に一致していると判定された場合には、ステップＳ１１４に移行し、求められたパラメータは、そのセグメントに対応する全ての基本セグメントの記憶領域に同様に順次記憶されて行く。一方、上記パラメータから構成される当てはめ曲線と、原時系列データの対応部分との差演算により求められる残差時系列データは次の下位のレイヤーのセグメントデータとして設定された後、ステップＳ１０９の処理が終了する。
【００９２】
一方、ステップＳ１１３において十分に一致していないと判定された場合には、ステップＳ１１５に移行し、セグメント解析結果は破棄し、そのセグメントデータを、そのまま次のレイヤーのセグメントに設定した後、ステップＳ１０９の処理が終了する。
【００９３】
このような処理を行うことにより、一括した解析処理では、十分に一致する当てはめ曲線が得られないセグメントデータは、順次、下位の２つのセグメントに分割されて設定され、次は分割された夫々のセグメントデータとして解析処理されるため、一括した解析処理では、データ構造の変動により十分に一致する当てはめ曲線が得られないセグメントデータであっても、最終的にはサブレイヤーの基本セグメントにおいて、十分に一致する当てはめ曲線を求めることができる。
【００９４】
ステップＳ１０９の処理が終了すると、次にステップＳ１１０に移行して、全てのセグメントの解析が終了したか否かを判定し、終了していない場合には、ステップＳ１０１に移行して、次のセグメントの解析を行い、全てのセグメントの解析処理が終了している場合には、セグメント解析が終了する。
【００９５】
次に、上述したとおり、一括解析結果構成手段１３ａは、再構成スペクトル計算手段１０とその他計算手段１２を利用して図３に示すような解析結果１０を構成し、またセグメント解析結果構成手段１３ｂは、再構成スペクトル計算手段１０と、頂位位相推移計算手段１１と、その他計算手段１２を利用して図３に示すような解析結果１１を構成するものであり、そこでこれらの計算手段６、７の動作を次に説明する。
【００９６】
まず図９は一括解析結果構成手段１３ａにより利用されている場合の再構成スペクトル計算手段１０の動作を説明する模式的流れ図であり、まずステップＳ２０１では、全ての記憶領域１〜記憶領域２Ｎに保存されている全てのパラメータ（振幅aj及び周期Tj）と残差時系列データεx(t)を取り出す。尚、取り出されたパラメータの個数は図中に示すように便宜的にＬ個としている。
【００９７】
次いで、ステップＳ２０２において、時系列データのエネルギーを保存するため、Ｌ個のパラメータ、即ち、夫々Ｌ個の振幅ajの規格化を行う。この規格化は、解析シーケンス決定手段６よりサブレイヤーが構成されておらず、従って全ての基本セグメントのデータ期間が等しい長さの場合には、全てのモードの振幅ajの自乗を基本セグメントの個数２Ｎで除して、その平方根を新たな振幅ajとすることで規格化を行うことができる。一方、サブレイヤーが構成されており、従ってデータ期間の長さが異なる基本セグメントが含まれる場合には、全ての基本セグメントについて、対応する記憶域に記された全モードの振幅ajを自乗すると共に、この値に、モードの属する基本セグメントの長さと、全期間の長さの比を乗じ、得られた積の平方根を新たな振幅ajとすることで規格化を行うことができる。
【００９８】
このようにして振幅が規格化されたＬ個のパラメータ（振幅aj及び周期Tj）の夫々により、ステップＳ２０３において、中心周波数1/Tjにピークをもつ孤立ピーク列Pj、即ち図に示される孤立ピーク列（ P1 … Pj … PL ）を構成する。
【００９９】
一方、残差時系列データεx(t)については、ステップＳ２０４において、最大エントロピー法のスペクトル密度（MEM-PSD）の計算を行い、得られたスペクトル密度と、ステップＳ２０３で得られた孤立ピーク列をステップＳ２０５において合成することにより、再構成スペクトルが計算される。なお、Ｓ２０４でラグ０とおくと残差時系列のパワーが全周波数帯に均一に分布するスペクトルとなる。
【０１００】
この際、本発明では、理論形式では有限離散値を扱うコンピュータで構成することができない線スペクトル列となる時系列データのスペクトルを、ピーク周波数から離れるに従って急減する連続関数、例えば余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルをδ函数の原型として利用して再構成することにより、コンピュータを利用して、理論と整合する時系列データの解析を行うことができる。また、これにより、２つの時系列データのクロススペクトルの絶対値と等しくなる、２つの時系列データの積を構成することも可能となる。
【０１０１】
次に図１０はセグメント解析結果構成手段１３ｂにより利用されている場合の再構成スペクトル計算手段１０と頂位位相推移計算手段１１の動作を説明する模式的流れ図である。
【０１０２】
まずステップＳ３０１では、基本セグメントに対応する記憶領域１〜２Ｎのいずれかが選択され、それに記憶されている、対応する基本セグメントの代表時刻、パラメータ、残差データが取り出される。最初は最も早い代表時刻の記憶領域１の記憶内容が取り出され、次いで後述する繰り返し処理により、順次、時間が経過した基本セグメントに対応する記憶領域の内容が取り出されて、全ての記憶領域の内容が取り出される。
【０１０３】
取り出された記憶領域のデータ中、パラメータ（振幅aj及び周期Tj）と残差時系列データεx(t)は再構成スペクトル計算手段１０に入力され、またパラメータ（周期Tj及び頂位位相φj）は頂位位相推移計算手段１１に入力されて処理に供される。
【０１０４】
再構成スペクトル計算手段１０に入力されたパラメータ（振幅aj及び周期Tj）は、一括解析結果構成手段１３ａにより利用されている場合の上述したステップＳ２０２の振幅の規格化の処理は行われず、ステップＳ２０３と同様に、ステップＳ３０２において、中心周波数1/Tjにピークをもつ孤立ピーク列Pj、即ち図に示される孤立ピーク列（ P1 … PL ）を構成する。
【０１０５】
一方、残差時系列データεx(t)については、ステップＳ３０３において、最大エントロピー法のスペクトル密度（MEM-PSD）の計算を行い、得られたスペクトル密度と、ステップＳ３０２で得られた孤立ピーク列をステップＳ３０４において合成することにより、再構成スペクトルが計算される。なお、Ｓ３０３でラグ０とおくと残差時系列のパワーが全周波数帯に均一に分布するスペクトルとなる。
【０１０６】
一方、頂位位相推移計算手段１１に入力されたパラメータ（周期Tj及び頂位位相φj）は、ステップＳ３０５において解析されて、周波数と頂位位相との対応関係から、基本セグメントの代表時刻に対する頂位位相の周波数依存が得られる。
【０１０７】
これらの得られた解析結果は、ステップＳ３０６において所定の記録領域に、基本セグメントの代表時刻の属性として記録される。
【０１０８】
次いでステップＳ３０７では、全ての基本セグメントに対応する記憶領域のデータが解析されたか否かが判定され、解析されていない記憶領域が残っている場合には、結合子Ｂを経てステップＳ３０１に移行して、次の記憶領域の解析が行われる。一方、ステップＳ３０７において、全ての基本セグメントに対応する記憶領域のデータが解析されている場合には、終了処理が行われる。
【実施例１】
【０１０９】
次に以上に説明した本発明の解析装置を構成する単一系列解析手段の動作を、具体例につき説明する。
解析する時系列データは、図１１に示すもので、データは48時間201点の時系列データである。図から分かるように、24時間モードは、時刻24時においてその位相を反転させている。
【０１１０】
まず図１１に示される時系列データが解析装置１に入力されると、解析シーケンス決定手段２は上述したように図４、図５に示される処理過程を経て解析シーケンスを決定する。
【０１１１】
この解析装置１では、演算処理の精度上、例えば最大データ点数が300点に設定されているとすると、入力された原時系列データのデータ点数は201点であり、原時系列データのデータ点数が最大データ点数以下であるため、図４のステップＳ３を経て、レイヤーは一層のみ設定され、このレイヤーに一つのセグメントが配置される。
【０１１２】
次いでステップＳ６において、原系列データの全期間データがセグメントに取り込まれる。次いでステップＳ７において、最大ラグ値(200点)までバーグの係数Pmの推移が計算される。図１２が係数Pmの推移を示すものである。
【０１１３】
次いでステップＳ８において係数Pmの挙動が解析される。まず、最初と最後の係数Pmの比が計算され、その結果、それが10桁を超えるために「周期性の強いデータである」と判断される。
【０１１４】
次いで解析シーケンス決定手段２は、係数Pmの推移を解析し、通常、係数Pmは最初急激に減少したのちに最大ラグ値まで漸減するのに対して、入力された原時系列データでは、ラグ値3〜100(対応するデータ点数は4〜101点)の領域で殆ど係数Pmが減少せず、直後に急激に減少して、その後最大ラグ値ちかくまで再びプラトー領域を形成することから、101点を超える付近でデータ構造が変化していることを検出する。
【０１１５】
従って解析シーケンス決定手段２は、ステップＳ１１の処理に移行し、サブレイヤーを構成し、データ構造の変化している個所において２分割したセグメントを構成する。
【０１１６】
次いでステップＳ６に移行して、ステップＳ１１において分割されたセグメントの夫々につき、順次ステップＳ８により、データ構造の均一性を評価する。
【０１１７】
尚、分割されたサブレイヤーのセグメントの境界にデータが位置しているため、このデータは隣接するセグメントに重複して取り込まれ、結果としてサブレイヤーの各セグメントのデータ点数は、いずれも101点となる。
【０１１８】
まず時間が早い方のセグメントデータについてのステップＳ８の評価において、係数Pmの推移を計算すると図１３に示すように推移し、係数Pmの値はラグ値0から7までで14桁ほど減衰し、上述したようなプラトー領域での急激な減衰領域の存在を検出しない。従って解析シーケンス決定手段２は、ステップＳ９において、このセグメントを基本セグメントとして設定する。ステップＳ１０からステップＳ６を経て、時間が遅い方のセグメントデータについても同様な評価がなされ、このセグメントについてもステップＳ９において基本セグメントとして設定する。
【０１１９】
次いで、ステップＳ１０において評価が完了と判定され、ステップＳ１２に移行して、ステップＳ９において基本セグメントとして設定された全てのセグメントに対応して、記憶手段に記憶領域が設定され、初期化がされて、解析シーケンスが完成する。尚、記憶領域の初期化において、それに対応する基本セグメントがステップＳ８において「極めて周期性が強い」と判断された場合には、記憶領域にその情報が書き込まれる。
【０１２０】
図１４は解析シーケンス決定手段２が決定したレイヤーとセグメントの構成を模式的に示すものであり、レイヤーは、時系列データの全期間に対応する一つのセグメント１が配置されるレイヤー０の下位に、サブレイヤー１が構成され、このサブレイヤー１の２つのセグメントが基本セグメント１、２として設定されている。そしてこれらの基本セグメント１、２に対応して記憶手段に記憶領域１、２が設定されている。
【０１２１】
次いでデータ解析手段３は、図７の流れ図に示す処理過程により各セグメントの解析を行う。この実施例においては、まず最初にレイヤー０のセグメントについてセグメント解析が実行されるが、レイヤー０のセグメント１が含むサブレイヤーの２つの基本セグメント１、２はともに極めて周期性の強いデータであると判別されていることから、このレイヤー０のセグメントの解析は省略することができるが、以降の説明では、レイヤー０のセグメント１についても解析を行うものとしている。
【０１２２】
まずデータ解析手段３は、ステップＳ１０１においてセグメント１に原時系列データの対応するデータを取り込んでセグメントデータとする。この実施例においては、このセグメント１のデータは、原時系列データそのもので、２０１点のデータ点数を有している。
【０１２３】
ここでステップＳ１０２において、最大エントロピー法により最小ラグ値から最大ラグ値までバーグの係数Pmを求め、夫々のラグ値に対するPmの推移を算出する。この実施例においては、同様の計算を解析シーケンス決定手段２において行って、図１２に示されるような結果を得ているので、本発明の単一系列解析手段２ａ，２ｂにおいて、実行された計算を適宜期間保持するように構成することにより、再計算を行う必要がなくなる。
【０１２４】
次いでステップＳ１０３において、図１２に示される推移から特徴的な挙動を示すラグ値を複数選択する。図１２において最大の特徴はラグ値2〜101のプラトー領域の存在で、周期性の強いデータなので、データ解析手段３は、自体に構成している知識ベース等を利用した選択手段により、小さめのラグ値を選択することとし、例えばラグ値2・50・101の3つのラグ値を選択し、ステップＳ１０４において、これらのラグ値に対応する３つのMEM-PSDを計算する。その結果は、図１５に示すとおりであり、左から順にラグ値2・50・101での結果を示す。
【０１２５】
次いでステップＳ１０５では、主要ピークとしてラグ値2では周期の初期値24.66時間(パワー49.80)、ラグ値50では28.06時間(同49.76)、ラグ値101では2モードで34.48時間(同29.81)と18.63時間(同17.31)を抽出し、これを次の非線形最小自乗法による計算の初期値とする。
【０１２６】
次いでステップＳ１０６において、上述した3組の初期値に対して非線形最小自乗法により、夫々の場合の最適当てはめ曲線と曲線パラメータの探索が行われる。この場合、非線形最小自乗法解析においては、初期値の周期の近傍に周期値を探索する領域をあらかじめ設定し、初期値が真の値の近傍にあること、およびその近傍で残差自乗和が唯一の極小値を持つことを期待して計算が行われる。
【０１２７】
即ち、ステップＳ１０６においては、探索過程を監視し、探索中に近傍領域から逸脱する場合や、複数の極小値をもつ場合には、初期値そのものが不適切であると判定して、その結果を使用しない。そこで、ラグ値2と50の初期値は不適切であると判定され、唯一、ラグ値101の初期値に対してのみ計算結果を使用するものとして、最適あてはめ曲線とその曲線パラメータが算出される。図１６（ａ）は算出された最適あてはめ曲線であり、また（ｂ）は原時系列データとの残差曲線を示すものである。
【０１２８】
次いでステップＳ１０７においては、当てはめ曲線のパラメータを、初期値と比較して一致度を評価して、ステップＳ１０６で算出した最適あてはめ曲線の妥当性が判定される。即ち、図１７は、上記初期値のピーク周波数の逆数とパワーを、最終的に得られた曲線パラメータの周期・振幅とを比較したものである。尚、この場合、図中に示す一致度とは、上述したとおり、モードのパワー(振幅の自乗の1/2)をピークパワーで除した値としている。
【０１２９】
図１７から分かるように、一致度は36時間モードでは0.967、18時間モードでは0.803となる。18時間モードの一致度は余り高くないが、ステップＳ１０６において算出された唯一の最適あてはめ曲線とそのパラメータであるので、データ解析手段３は、その結果を、ステップＳ１０９、即ち、解析結果評価手段４による評価の処理過程に出力する。
【０１３０】
まずステップＳ１１２においては、ステップＳ１０６〜Ｓ１０８を経て得られた最適あてはめ曲線が再構成され、そのデータ期間に対応する原時系列データとの差の演算を行い、残差曲線を求めると共に、当てはめ曲線の再現性を評価する諸量、ここでは残差の標準偏差と残差パワーの時間依存が求められる。
【０１３１】
図１８は原時系列データ（点線）と最適あてはめ曲線（実線）とを比較したもの、また図１９は残差曲線を、２０個の区分領域毎に平均したものである（白丸）。残差の標準偏差は0.836と大きく、また図１８から、一見して当てはめが良好でないことが分かる。そして図１９から、残差はゼロの回りに一様に分布しておらず、時間依存性が認められる。
【０１３２】
このような残差の標準偏差の値の大小や、残差の時間依存性等から、データ解析手段３は、ステップＳ１１３において、最適当てはめ曲線と、原時系列データの対応部分が十分に一致していないと判定し、ステップＳ１１５に移行する。即ち、それまでの解析対象であったレイヤー０のセグメント１のセグメントデータは、そのまま次のレイヤー１のセグメント１、２の解析対象データとして繰り越して設定される。
【０１３３】
尚、上述したとおり、原時系列データのデータ点数をバンチング処理により低減してセグメントデータを構成する場合には、下位のレイヤーのセグメントに設定するデータは、必ずしも上位のレイヤーのセグメントデータと同一でなくとも良い。
【０１３４】
次いでデータ解析手段３は、ステップＳ１０９からステップＳ１１０を経て、ステップＳ１０１に移行し、次にサブレイヤー１の基本セグメント１、２の解析が順次行われる。
【０１３５】
まず基本セグメント１については、ステップＳ１０１においてセグメントデータとして原時系列データの1〜101点までのデータが取り込まれる。
【０１３６】
次いでステップＳ１０２における計算により、図１３に示すような係数Pmの推移が得られる。尚、この場合には、数値演算精度のためにラグ値7までの推移が求められる。
【０１３７】
次いでステップＳ１０３では、ラグ値として3と7の二つが決定され、ステップＳ１０４において決定されたラグ値でMEM-PSDが計算される。その結果は、図２０に示すとおりであり、左から順にラグ値３とラグ値7の結果である。
【０１３８】
次いでステップＳ１０５では、ラグ値2では初期周期値24.00時間(パワー49.66)を、ラグ値7では初期周期値24.00時間(同50.14)が抽出され、これらの二組の初期値を用いて、ステップＳ１０６において非線形最小自乗法により二本の最適あてはめ曲線と曲線パラメータが算出される。この場合には、いずれも周期24.00時間、振幅10.00、頂位位相6.00時が得られる。
【０１３９】
次いでステップＳ１０７では、上述した複数の曲線パラメータの組が比較されるのであるが、この実施例では、いずれのパラメータも一致し、且つ当てはめも良好であるので、その最適当てはめ曲線と、そのパラメータがステップＳ１０８を経てステップＳ１０９に入力される。
【０１４０】
ステップＳ１１２において得られる残差の標準偏差は、10−10と十分に小さく、且つ、残差のパワー推移に時間依存がないため、ステップＳ１１３において、今度は、最適当てはめ曲線と、原時系列データの対応部分が十分に一致していると判定し、ステップＳ１１４に移行する。
【０１４１】
ステップＳ１１４においては、解析結果のパラメータを記憶領域１に保存する。尚、この実施例と異なり、下位に更にサブレイヤーがある場合には、残差が計算されて、それを下位のサブレイヤーに設定するべく保存がされる。
【０１４２】
次いでステップＳ１１０を経て、ステップＳ１０１に移行して、次の基本セグメント２についても同様な解析処理が行われ、こうして、記憶領域１、２に、原時系列データを一般化三角多項式で記述し、且つ、以降の解析に必要な諸量の全てが蓄積される。
【０１４３】
そして記憶領域１、２に蓄積されている原時系列データの解析に必要な諸量は、上述したとおり、一括解析結果構成手段１３ａまたはセグメント解析結果構成手段１３ｂにより処理されて、図３に示すような解析結果１４と１５を構成することができる。
【０１４４】
一例として、一括解析結果構成手段１３ａにより利用されている場合の再構成スペクトル計算手段６による再構成スペクトルの算出の流れを説明する。
まず図９において、ステップＳ２０１では、記憶領域１と２からパラメータ（一般化三角多項式のパラメータ）を取り出される。この際、取り出されたパラメータは２つで、いずれも周期24時・振幅10・頂位位相6時である。
【０１４５】
次いでステップＳ２０２で振幅が規格化される。この規格化は、２つのパラメータそれぞれについて、振幅の自乗＝100を基本セグメント数(＝2)で割って、その平方根(＝√50)をあらたな振幅とすることで行われる。
【０１４６】
次いでステップＳ２０３における孤立ピークのPSDの計算により、同じピーク周波数(＝1/24時間)をもつ二つのピークスペクトルが得られる。この孤立ピークの計算は、ピーク周波数の逆数の周期をもつ三角関数のMEM-PSDを求めておき、その面積を各モードのパワーに一致させる規格化により実現する。図２１は余弦関数のMEM-PSDの例である。
【０１４７】
このように本発明では、孤立ピーク列は、δ函数の原型を用いて全周波数帯に裾の広がる夫々一つのスペクトルピークを構成することにより、パラメータから得られる各周波数に対応する孤立ピーク群を構成するもので、こうして理論形式では有限離散値を扱うコンピュータで構成することができない線スペクトル列となる時系列データのスペクトルを、ピーク周波数から離れるに従って急減する連続関数、例えば余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルをδ函数の原型として利用して再構成することにより、コンピュータを利用して、理論と整合する原時系列データの再構成スペクトルを求めることができる。尚、この手法は、後述するように、相互相関解析手段３における計算においても同様に用いられる。
【０１４８】
一方、各記憶領域１、２から取り出された基本セグメント１、２の残差成分から全測定長にわたる残差曲線を求め、次いでステップＳ２０４において、残差のMEM-PSDが計算される。次いで、ステップＳ２０５において、孤立ピーク群のスペクトルと残差のMEM-PSDとが重畳されて、最終的に原時系列データの一般化三角多項式表現と整合する再構成スペクトルが得られる。なお、Ｓ２０４でラグ０とおくと残差時系列のパワーが全周波数帯に均一に分布するスペクトルとなる。
【０１４９】
相互相関解析手段３の実施の形態を説明すると、図２２に示す相互相関解析手段３は、一般化三角多項式パラメータ１６ａ，１６ｂからパワースペクトル密度を求めるパワースペクトル密度計算手段１８ａ，１８ｂと、クロススペクトル（コスペクトル、クオドスペクトル）とコヒーレンスを含む相互相関を定量化する諸量２０を求める相互相関諸量計算手段１７と、相互相関諸量計算手段１７により求められたコヒーレンスを理論値と比較して、それらの差が所定値以下の場合に、その差をなくすようにクロススペクトルを補正する整合性評価補正手段２１と、整合性評価補正手段２１において評価した差が所定値以上の場合に、単一系列解析手段２ａ，２ｂに再解析の指令２３ａ，２３ｂを出力し、所定値以下の場合には、相互相関諸量計算手段１７により計算された諸量２０と、パワースペクトル密度計算手段１８ａ，１８ｂにより計算されたパワースペクトル密度１９ａ，１９ｂ及び整合性評価補正手段により補正されたクロススペクトルを出力する分岐出力手段２２を構成したものである。
【０１５０】
まず、相互相関解析手段３に入力される２系列の時系列データは、上述と同様に夫々、次式の一般化三角多項式で表現される。なお、クロススペクトルに関わる次式以下の展開は後述の非特許文献3で詳述される。
【数６】
2010108402.tif 000008
【０１５１】
上式のx(t)とy(t)は、上述したとおり有限長の時系列データを記述するものであるが、仮に、その記述を±∞に延長したとき、それら延長時系列データのパワースペクトル密度Ｓxx(f)とＳyy(f)は次式で表される。
【数７】
2010108402.tif 000009
上式において、δはディラックのデルタ函数、εxx(f)は雑音成分のスペクトルである。ここでΣ…の項はδ函数の引数が０となる周波数においてのみ、有効な値（∞）をもつ。
【０１５２】
即ち、パワースペクトル密度計算手段１８ａ，１８ｂにおいては、上式に基づき、夫々の時系列データのパワースペクトル密度Ｓxx(f)とＳyy(f)が求まる。この際には、上述したように、線スペクトル列となる時系列データのスペクトルを、ピーク周波数から離れるに従って急減する連続関数、例えば余弦函数に微弱な雑音を加えた時系列をラグ２にて最大エントロピー法のスペクトル密度の計算を行って得られたスペクトルをδ函数の原型として利用して再構成することにより、有限離散値を扱うコンピュータを利用して計算を行うことができる。このδ函数の原型は全周波数帯で正の値をもつ。
【０１５３】
すなわち上式はδ函数の原型を用いて次式となる。
【数８】
2010108402.tif 000010
ここに
【数９】
2010108402.tif 000011
と
【数１０】
2010108402.tif 000012
はそれぞれ周波数
【数１１】
2010108402.tif 000013
と
【数１２】
2010108402.tif 000014
に極大値をもつδ函数の原型である。ここでΣ…の項はD(…)函数の引数がゼロとなる周波数およびその近傍で有効な値（＜∞）をもつ。
【０１５４】
一方、相互相関を定量化する諸量のうち、そのクロススペクトルやコヒーレンス等は次式で表される。
【数１３】
2010108402.tif 000015
上式において、Ｋxy(f)はコスペクトル、Ｑxy(f)はクオドスペクトルである。また、
【数１４】
2010108402.tif 000016
と
【数１５】
2010108402.tif 000017
は時系列x(t)とy(t)で、一般化三角多項式成分と直交する純非決定論的成分のコスペクトルとクロススペクトルである。本発明による時系列のあてはめは常に良好であるので、これらは無視することができる。従って、自乗コヒーレンス（又は単にコヒーレンス）coh2(f)は、
【数１６】
2010108402.tif 000018
となる。フェイズθxy(f)と時間遅れτxy(f)は次式で与えられる。
【数１７】
2010108402.tif 000019
【０１５５】
同様に回転スペクトルについても、例えば反時計回りの回転スペクトルＳ＋(f)と回転係数ＣR(f)は次式で表される。
【数１８】
2010108402.tif 000020
その他の回転スペクトルの諸量もＫxy(f)とＱxy(f)の函数であることから、上述したクロススペクトルの諸量と同様に、x(t)とy(t)を記述するパラメータにより表される。
【０１５６】
以上の相互相関を定量化する諸量２０は、相互相関諸量計算手段１７において計算される。
【０１５７】
本発明では、相互相関解析手段３の整合性評価補正手段２１において、計算結果と理論上期待される値との一致を評価する。具体的には、まず、パワースペクトル密度計算手段１８ａ，１８ｂにおいて求められたパワースペクトル密度Ｓxx(f)，Ｓyy(f)と、相互相関諸量計算手段１７において求められたクロススペクトルＫxy(f)，Ｑxy(f)とから求められるコヒーレンスcoh2(f)を理論値１と比較し、それらの差が所定値以下の場合に整合性が高いと判定する。
【０１５８】
なお、理論上コヒーレンスは全周波数帯で１となる。この重要な知見は、２００８年１２月２５日刊行予定の下記非特許文献３において証明が記載される。すなわち、コヒーレンスは０〜１の間で連続的に変動するという従来の理解が誤りであったことが明らかとされる。
【非特許文献３】常磐野和男・大友詔雄・田中幸雄著、「最大エントロピー法による時系列解析  …ＭｅｍＣａｌｃの理論と実際…」、第２版、北海道大学図書刊行会、2008年12月25日（刊行予定）
【０１５９】
判定の結果、整合性が高くないと判定された場合には、それらの差を判定し、その差が微細と判定される場合は、整合性評価補正手段２１において、その差をなくすようにクロススペクトルを補正する。具体的には、コヒーレンスcoh2(f)，フェイズθxy(f)，時間遅れτxy(f)、パワースペクトル密度Sxx(f)，Syy(f)を固定して、物理量の間の定義式を満足するようにクロススペクトルKxy(f)，Qxy(f)を補正する。
【０１６０】
このようにして判定の結果、整合性が高いと判定された場合や、差が微細で補正がされた場合には、分岐出力手段２２により、相互相関諸量計算手段１７により計算された諸量と、パワースペクトル密度計算手段１８ａ，１８ｂにより計算されたパワースペクトル密度Sxx(f)，Syy(f)及び整合性評価補正手段２１により補正されたコスペクトルKxy(f)、クオドスペクトルQxy(f)が出力される。
【０１６１】
一方、整合性が高くないと判定された場合には、以上の計算結果は破棄されると共に、分岐出力手段２２を介して、夫々の単一系列解析手段２ａ，２ｂに上述した再解析指令を発する。
【０１６２】
こうして本発明では、相互相関解析手段３において、一般化三角多項式パラメータにより求められたコヒーレンスを理論値と比較して整合性を評価し、所定以上の整合性が得られない場合には、単一系列解析手段２ａ，２ｂに対して再解析指令を発して、より高精度の解析を行わせるという、フィードバックと繰り返し動作を有するため、２系列の時系列データの間の相互相関を定量化する諸量を非常に高精度に求めて出力することができる。
【０１６３】
以上説明したように、本発明においては、２系列の時系列データの夫々を単一系列解析手段２ａ，２ｂにおいて、まず最大エントロピー法を適用して高精度なパワースペクトル密度を求めると共に、その顕著なピークから、一般化三角多項式を解く初期値としての、項数と、各三角項の周期値を求め、次いで、これらの初期値により非線形最小自乗法を用いて、各三角項のパラメータを求めると共に、残差時系列データを求め、これらを加算して時系列データを再構成するという解析を行うことにより、何らの数値モデルも想定せず、夫々の時系列データのみを用いて高精度の解析を行え、これらの高精度の解析によって得られた２系列の時系列データの一般化三角多項式パラメータにより、相互相関解析手段３において、２系列の時系列データの間の相互相関を定量化する諸量を求めることができる。
【０１６４】
また本発明においては、２系列の時系列データの間の相互相関を定量化する諸量を求める他、予測(および解析)のための数値モデルを全く要さず、２系列の時系列データのみを用いて予測支援機能を発揮させることができる。
【０１６５】
即ち、図２３は本発明の装置を利用して構成した予測支援機能を模式的に示すもので、実線の矢印はデータの流れを示し、破線の矢印は制御（フィードバック）の流れを示すものである。
【０１６６】
なお、予測支援機能の目的は、多数の時系列データを解析してその中に、ある時系列Ａの将来挙動が別の時系列Ｂの現在までの挙動により予想し得る一対の時系列ＡとＢを見いだし、現在のＢによりＡの将来を予測することである。
【０１６７】
まず、予測支援機能の動作の概要を説明すると、図中符号２４ａ，２４ｂ，２４ｃ，２４ｄ，…は複数の系列１，２，３，４，…の時系列データと、本発明の装置による解析結果の諸量を示すもので、予測対象系列の参照系列の決定処理手段２５においては、多数の時系列とその解析結果を検査して共通の周期をもつ有用なモード(指標モード)をともに含む一対の時系列(予測対象系列と参照系列)を抽出し、次いで算出処理手段３２において、基底変動や指標モードの解析区間と予測区間における挙動など、予測支援に直接用いる諸量を求めて、それらの諸量３３を出力する。
【０１６８】
上記動作を更に詳細に説明すると、まず処理手段２６においては、同じ周期の三角函数波を含む一対の時系列データを探索する(それらを以下指標モードとする)。指標モードとしては、(a)長周期であること、(b)指標モードの振幅がその他のモードのそれに対して十分な大きさのあること、(c)頂位位相が接近または一方に半周期加えた位相が接近していること(これは必須ではない)、等が求められる。こうして位相が先行する指標を参照指標モード２８、後行する指標を予測対象指標モード２７として求める。
【０１６９】
そして本発明では、参照指標モードの現在の挙動から予測対象指標モードの将来の挙動を推測する。このため、それぞれの時系列について指標モードが時系列全体の挙動とその他の各モードの挙動に対して重要な立ち位置を占めていることが肝要である。処理手段２９においては、処理手段２６において見いだした指標モードをこの視点で評価し、図中右側の破線矢印に示すように、不適切ならば同じ時系列の組の別の指標モードを探索するか、別の時系列の組の探索に移る。
【０１７０】
更に処理手段２９においては、指標モードよりも長周期のモードの一群と短周期のモードの一群について、指標モードとの兼ね合いが検討される。(a)長周期モードの一群については、極長周期の少数個のモードで時系列全体の基本的挙動(基底変動)を記述できることが望ましい。なぜならば、そのような長周期モードは解析区間の外側、予測区間においても大きく変動しないことが期待されるからである。また、(b)短周期モードの一群については、時系列の挙動への寄与が極力少ないことが望ましい。例えば一週間のリズムなどの外部要因によることが明らかなモードを除いては、それらの寄与が大きい場合、予測区間において予測対象指標モードが期待通りに推移しても、全モードの和・実際の観測値の推移との間に大きな乖離が予想されるからである。
【０１７１】
以上の処理により、予測支援に用いる一対の時系列と解析結果３０、３１として得ることができ、次いで処理手段３２においては、こうして得られた２つの時系列から、実際に予測支援に用いるさまざまな値、即ち上述した諸量３３を得るための処理が行われる。この処理においても、結果が不満足な場合にはその一対の時系列は予測支援に適さないとして、図中左側の破線矢印で示すように処理手段２５に指令がフィードバックされる。
【０１７２】
処理手段３２においては、（１）に示すように、それぞれの時系列について、その三角函数波のうちで基底変動を構成するモードが抽出され、それに指標モードが加えられる。また（２）に示すように、少数個の長周期モードで観測期間全体にわたる時系列の基本的挙動を再現できるかが調べられ、できない場合には、時系列(の対)は破棄され、処理手段２５における処理に移行する。
【０１７３】
また、処理手段３２においては、（３）に示すように、二つの指標モードがそれぞれの時系列の解析区間全体において安定しているかが調べられる。即ち、解析区間の前半・真ん中・後半といった具合にデータを複数のセグメントに分割し、その各セグメントにおいて、与えられた周期値をもつ一般化三角多項式でデータに当てはめを行い、いずれのセグメントでも少なくとも指標モードが安定して存在するか、または変化がゆるやかであることを確認する。そして、加えて、（４）に示すよう、２つの指標モードがそれぞれの時系列の解析区間全体において安定した寄与を示すことを確認する。これは例えば、指標モードの一周期ごとに時系列の平均値を求め、その推移がゆるやかであるかなどを確認しておこなう。更に、その他、時系列に応じて(外部からの強制周期の有無、利用可能な付加情報の有無など)、二つの系列を予測対象系列・参照系列として扱って良いかが検討される。
【０１７４】
次に本発明に係る２系列の時系列データの解析装置による解析例を説明する。
【実施例２】
【０１７５】
まず図２４は、共に6つの正弦波に雑音を重畳した次式で示される２つの時系列を解析した例である。
【数１９】
2010108402.tif 000021
図中(a)の時系列データは、
【数２０】
2010108402.tif 000022
に(3,2)、(2,3)、(2,5)、(1,7)、(1,9)、(2,13)を、(b)の時系列データは(3,2)、(2,2.5)、(2,5)、(1,7)、(1,10)、(2,13)を夫々パラメータとして与え、ε(t)としてパワー10%の白色ノイズライクな雑音を重畳して得たもので、サンプリング間隔は0.1とし、時刻0〜99.9までの1000点のデータである。
そして図中(c)はクロススペクトルの絶対値、(d)はコスペクトル、(e)は頂位時刻の差(時間遅れ)、および(f)はクオドスペクトルの解析結果である。なお、時間遅れは余弦函数の頂位時刻について計算されている。
上述したように与えたパラメータから、共通の周期である2、5、7、13、すなわち周波数1/2、1/5、1/7、1/13で相互相関を持つことが期待されるが、図の(c)〜(f)ではまさにこの4つの周波数近傍に相互相関が認められる。
図中の(g)に、解析して得たそれぞれの最適あてはめ曲線のパラメータと時間遅れを示しているが、パラメータは時系列を発生する際に設定した値とよく一致していることが分かる。また時間遅れについても、設定値から期待される値によく一致している。大きな雑音の重畳にもかかわらず理論値と一致する結果が得られたことは、MEM-PSDのピークとピークパワーを初期値とする非線形最小自乗法を解くという、本発明に係る解析の手法が正しく機能しいてる証左である。
【実施例３】
【０１７６】
次に図２５は、図中(a)に示す振幅変調する時系列データと、(b)に示す正弦波時系列データの２つの時系列を解析した例であり、夫々の時系列データは次式により示される。
【数２１】
2010108402.tif 000023
(a)、(b)は、ｆ0=0.5、ｆｍ=0.25、サンプリング間隔0.1とし、時刻0〜29.9までの300点のデータである。
上式より、x(t)は周波数0.5に主たるピークを、周波数0.25と0.75に、振幅変調に由来するピークを持つ。また、y(t)は周波数0.5にピークを持つ。そして相互相関は周波数0.5においてのみ期待される。
図２４と同様に(c)はクロススペクトルの絶対値、(d)はコスペクトル、(e)は頂位時刻の差(時間遅れ)、および(f)はクオドスペクトルの解析結果である。
この解析で得たクロススペクトルの絶対値(＝二つのスペクトルの積の平方根)では、振幅変調に由来する周波数0.25と0.75の位置を含め、3つのピークが認められる。これは解析において、孤立ピークを全周波数帯にわたり構成した結果である(全周波数帯に対して構成しなければ振幅変調に由来するピークはつぶれてしまう)。
また、その他の3つの図、(d)のコスペクトル、(f)のクオドスペクトル及び(e)の頂位時刻の差(＝時間遅れ)では、いずれも周波数0.5においてのみ、相互相関が認められる。
図中(g)に示すように、解析で得られたそれぞれの最適あてはめ曲線のパラメータと時間遅れの値は、いずれも設定したパラメータ値をよく再現していることが分かる。
【実施例４】
【０１７７】
次に図２６は、次式で示される図中(a)のうなりの時系列x(t)と、周期的方形波y(t)の２つの時系列を解析した例である。
【数２２】
2010108402.tif 000024
ここにf1=1/50、f2=1/51である。また、y(t)は周期50、頂位時刻12.5で値1または-1からなる周期的方形波である。サンプリング間隔を1とし、時刻0〜1299まで1300点のデータを発生して解析を行った。
計算式より、周波数1/50においてのみ、相互相関が期待される。解析結果を示す(c)のクロススペクトルの絶対値はこの位置にピークを持つと共に、周期的方形波の高周波モードの位置にピークを持つ。ここで、注目すべきは、うなりのスペクトルは周波数1/50と1/51の近接ピークをもつが(非特許文献2、pp.198-199、または解析例3の図下段の数表参照)、クロススペクトルの絶対値の計算結果には理論式から期待される周波数1/50のピークのみが反映されていることである。本解析システムの周波数分解能の高さの証左である。
図中(g)に示すように、解析で得られたそれぞれの最適あてはめ曲線のパラメータと時間遅れの値は、いずれも設定したパラメータ値をよく再現していることが分かる。ここに時間遅れは+1.135となり、理論値0とは、周期に対する誤差約2%で一致している。
【産業上の利用可能性】
【０１７８】
以上に説明したように、本発明においては、２系列の時系列データの夫々を単一系列解析手段において、最大エントロピー法と非線形最小自乗法を用いて、各三角項のパラメータを求めると共に、残差時系列データを求め、これらを加算して時系列データを再構成するという解析を行うことにより、何らの数値モデルも想定せず、夫々の時系列データのみを用いて高精度の解析を行え、これらの高精度の解析によって得られた２系列の時系列データの一般化三角多項式パラメータにより、相互相関解析手段において、２系列の時系列データの間の相互相関を定量化する諸量を求めることができる
【０１７９】
本発明に係る解析装置においては、実在の時系列データに対して、時間軸上での一般化三角多項式表現と周波数軸上での再構成スペクトル表現を機械的に計算する。このとき、時系列データ・一般化三角多項式表現・再構成スペクトル表現それぞれの間には、コンピュータ上に構築された数値演算処理システムの能力の制限の元で、理論上の整合性が完全に保証される。すなわち、任意の実データに対して時間軸上および周波数軸上のもっとも整合的で透徹した理解が機械的に得られる。
【０１８０】
また、本発明に係る解析装置では、何らの数値モデルも想定せず、時系列データのみを用いて高精度の解析を行え、しかもその解析は、通常のパーソナルコンピュータ程度の演算能力でも十分に実用的であるという効果があり、特に、解析対象の時系列データを、そのデータ点数やデータ構造に応じて適切にセグメント化して処理を行うことにより、解析する時系列データによらず、常に適切に高精度な解析を行えるという格別なる効果がある。
【０１８１】
そして、本発明は、実在の時系列データ、即ち、ある限られた期間に観測されたある点数の有限長離散時系列のデータに加えて、例えば地層断面の距離に対する輝度変化等の、空間的に連続するデータの解析にも利用できるものである。
【０１８２】
本発明を利用して構成される予測支援システムは、多数の時系列中から一対の時系列を見いだすことにより、一方の時系列の現在までの挙動から他方の時系列の将来挙動を予測することが期待される。
【図面の簡単な説明】
【０１８３】
【図１】本発明に係る解析装置の実施の形態を模式的に表した全体構成の説明図である。
【図２】本発明に係る解析装置を構成する単一系列解析手段の実施の形態を模式的に表した全体構成の説明図である。
【図３】本発明に係る単一系列解析手段から出力される解析結果を模式的に示すものである。
【図４】解析シーケンス決定手段における処理の流れを示す流れ図である。
【図５】解析シーケンス決定手段における処理の流れを示す流れ図である。
【図６】解析シーケンス決定手段により決定されたレイヤーとセグメントの構成例を模式的に示すものである。
【図７】データ解析手段における各セグメントの解析処理の流れ図である。
【図８】解析結果評価手段における各セグメントの解析処理の流れ図である。
【図９】一括解析結果構成手段により利用されている場合の再構成スペクトル計算手段の動作を説明する模式的流れ図である。
【図１０】セグメント解析結果構成手段により利用されている場合の再構成スペクトル計算手段と頂位位相推移計算手段の動作を説明する模式的流れ図である。
【図１１】解析する時系列データの一例を示す説明図である。
【図１２】バーグの係数Pmの推移を示す説明図である。
【図１３】サブレイヤーの時間が早いほうのセグメントデータについてのバーグの係数Pmの推移を示す説明図である。
【図１４】決定されたレイヤーとセグメントの構成を模式的に示す説明図である。
【図１５】３つのラグ値に対応して計算された３つのMEM-PSDの一部を示す説明図である。
【図１６】（ａ）は算出された最適あてはめ曲線、（ｂ）は原時系列データとの残差曲線を示す説明図である。
【図１７】初期値のピーク周波数の逆数とパワーを、最終的に得られた曲線パラメータの周期・振幅とを比較した説明図である。
【図１８】原時系列データ（点線）と最適あてはめ曲線（実線）とを比較した説明図である。
【図１９】残差曲線を、２０個の区分領域毎に平均して示す説明図である。
【図２０】２つのラグ値に対応して計算された２つのMEM-PSDを示す説明図である。
【図２１】δ函数の原型として利用する、余弦関数のMEM-PSDを示す説明図である。
【図２２】本発明に係る解析装置を構成する相互相関解析手段の実施の形態を模式的に表した全体構成の説明図である。
【図２３】本発明の解析装置を利用して構成した予測支援機能を模式的に説明する説明図である。
【図２４】本発明の解析装置により解析した例をしめす説明図である。
【図２５】本発明の解析装置により解析した他の例をしめす説明図である。
【図２６】本発明の解析装置により解析した更に他の例をしめす説明図である。
【符号の説明】
【０１８４】
１                  統合構造解析手段
２ａ，２ｂ          単一系列解析手段
３                  相互相関解析手段
４ａ，４ｂ          時系列データ
５ａ，５ｂ，５ｃ    諸量
６                  解析シーケンス決定手段
７                  データ解析手段
８                  解析結果評価手段
９                  制御手段
１０                再構成スペクトル計算手段
１１                頂位位相推移計算手段
１２                その他計算手段
１３ａ              一括解析結果構成手段
１３ｂ              セグメント解析結果構成手段
１４                一括解析結果
１５                セグメント解析結果
１６ａ，１６ｂ      一般化三角多項式パラメータ
１７                相互相関諸量計算手段
１８ａ，１８ｂ      パワースペクトル密度計算手段
１９ａ，１９ｂ      パワースペクトル密度
２０                諸量
２１                整合性評価補正手段
２２                分岐出力手段
２３ａ，２３ｂ      再解析の指令
２４ａ，２４ｂ，２４ｃ，２４ｄ    諸量
２５                決定処理手段
２６                処理手段
２７                予測対象指標モード
２８                参照指標モード
２９                処理手段
３０，３１          解析結果
３２                処理手段
３３                諸量
【図１】
2010108402.tif 000029
【図２】
2010108402.tif 000030
【図３】
2010108402.tif 000031
【図４】
2010108402.tif 000032
【図５】
2010108402.tif 000033
【図６】
2010108402.tif 000034
【図７】
2010108402.tif 000035
【図８】
2010108402.tif 000036
【図９】
2010108402.tif 000037
【図１０】
2010108402.tif 000038
【図１１】
2010108402.tif 000039
【図１２】
2010108402.tif 000040
【図１３】
2010108402.tif 000041
【図１４】
2010108402.tif 000042
【図１５】
2010108402.tif 000043
【図１６】
2010108402.tif 000044
【図１７】
2010108402.tif 000045
【図１８】
2010108402.tif 000046
【図１９】
2010108402.tif 000047
【図２０】
2010108402.tif 000048
【図２１】
2010108402.tif 000049
【図２２】
2010108402.tif 000050
【図２３】
2010108402.tif 000051
【図２４】
2010108402.tif 000052
【図２５】
2010108402.tif 000053
【図２６】
2010108402.tif 000054