図面 (/)

技術 異常診断装置および異常診断方法

出願人 株式会社明電舎
発明者 林孝則外山達斎
出願日 2016年4月28日 (4年8ヶ月経過) 出願番号 2016-091695
公開日 2017年11月2日 (3年1ヶ月経過) 公開番号 2017-198620
状態 特許登録済
技術分野 機械部品、その他の構造物または装置の試験 機械的振動・音波の測定
主要キーワード 参照データ数 診断フェーズ 統計値計算 合計演算 記号法 計算時刻 解析位置 散布状況
関連する未来課題
重要な関連分野

この項目の情報は公開日時点(2017年11月2日)のものです。
また、この項目は機械的に抽出しているため、正しく解析できていない場合があります

図面 (8)

課題

サンプルの多様性を確保し、長期間にわたって信頼性の高い異常診断を実施することができるようにする。

解決手段

予め異常診断対象から測定した時系列データを第1の定Q変換装置101によって定Q変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータ変換行列正規化係数)を作成する診断用パラメータ作成手段100と、新たに異常診断対象から測定した時系列データを第2の定Q変換装置201によって定Q変換して多変量サンプルを作成し、前記診断用パラメータ作成手段100により作成された診断用パラメータを用いて異常を診断する診断手段200と、を備えた。

概要

背景

異常診断対象機器の異常を診断する方法は、例えば特許文献1に記載のものが提案されている。

この特許文献1の方法では、時系列波形データを、短い時間毎に分割して、フーリエ変換して周波数成分の多変量データサンプルを多数作り、機器が正常と考えられる期間に予め計測した多数の多変量サンプルを主成分分析して主成分得点を求めるための固有ベクトルローテーション行列ともいう)を用意する。

診断の際には、計測データを前記同様にフーリエ変換して周波数成分の多変量データ・サンプルを作り、これを用意しておいたローテーション行列で変換して主成分得点を求め、それを正常時の主成分得点と比較することで異常診断を行っている。要するに診断毎に独立に主成分分析するのではなく、予め主成分分析しておいた結果を利用して、診断時には同じ基準で変換することにより比較を容易にするのが要点である。

尚、本発明で利用する定Q変換の成分演算方法は、例えば非特許文献1〜4に記載されている。

概要

サンプルの多様性を確保し、長期間にわたって信頼性の高い異常診断を実施することができるようにする。予め異常診断対象から測定した時系列データを第1の定Q変換装置101によって定Q変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータ変換行列正規化係数)を作成する診断用パラメータ作成手段100と、新たに異常診断対象から測定した時系列データを第2の定Q変換装置201によって定Q変換して多変量サンプルを作成し、前記診断用パラメータ作成手段100により作成された診断用パラメータを用いて異常を診断する診断手段200と、を備えた。

目的

本発明は上記課題を解決するものであり、その目的は、サンプルの多様性を確保し、長期間にわたって信頼性の高い異常診断を実施することができる異常診断装置および方法を提供する

効果

実績

技術文献被引用数
0件
牽制数
0件

この技術が所属する分野

ライセンス契約や譲渡などの可能性がある特許掲載中! 開放特許随時追加・更新中 詳しくはこちら

請求項1

予め異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段と、新たに異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段と、を備えたことを特徴とする異常診断装置。

請求項2

前記診断用パラメータ作成手段は、時系列データを定Q変換して多変量サンプルを作成する第1の定Q変換装置と、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部と、前記主成分得点を基にサンプルの指標となる統計量を得る第1の統計値計算部と、前記統計量を正規化して異常度とする第1の正規化部と、を備え、前記診断手段は、時系列データを定Q変換して多変量サンプルを作成する第2の定Q変換装置と、前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得る主成分計算部と、前記主成分得点を基に統計量を得る第2の統計値計算部と、前記統計量を、前記第1の正規化部で使用した正規化係数によって正規化して異常度を計算する第2の正規化部と、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする請求項1に記載の異常診断装置。

請求項3

前記診断用パラメータ作成手段および診断手段の第1、第2の定Q変換装置は、ともに、設定した周波数成分KM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Q変換における周波数軸係数との積和計算により求めた定Q変換成分の第1の演算量と、前記KMを超える高周波帯域において、時系列データと定Q変換における時間軸係数との積和計算により求めた定Q変換成分の第2の演算量とを含む合計演算量を算出し、前記KMを全周波数成分Kに各々設定して各KM毎の前記合計演算量を求め、前記各KM毎の合計演算量を評価し、該演算量が最小となるKMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、定Q変換における周波数成分Kが前記決定された周波数成分閾値KM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第1の計算法を選択し、前記周波数成分Kが前記決定された周波数成分閾値KMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第2の計算法を選択する計算法選択手段と、時系列データに対して、前記計算法選択手段により選択された第1の計算法、第2の計算法を実施して定Q変換の成分を求める定Q変換手段と、を備えたことを特徴とする請求項2に記載の異常診断装置。

請求項4

前記各定Q変換装置は、前記定Q変換における時間軸係数を、によって算出する時間軸係数算出部と、前記定Q変換における周波数軸係数を、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから(1/N)Sn,kとして算出する周波数軸係数算出部と、を備え、前記定Q変換手段は、前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第1の計算法を実施する低周波帯域計算部と、前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、を備えていることを特徴とする請求項3に記載の異常診断装置。

請求項5

前記各定Q変換装置は、を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)毎に算出して時間軸係数Tn,k,mを求め、該Tn,k,mを前記定Q変換における時間軸係数とする時間軸係数算出部と、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから得た(1/N)Sn,kを、前記解析位置pm毎に算出して周波数軸係数(1/N)Sn,k,mを求め、該(1/N)Sn,k,mを前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、前記定Q変換手段は、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で算出された周波数軸係数(1/N)Sn,k,mと前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施する低周波帯域計算部と、前記解析位置pm毎に、前記時間軸係数算出部で算出された時間軸係数Tn,k,mと時系列データxnとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、を備えていることを特徴とする請求項3に記載の異常診断装置。

請求項6

前記各定Q変換装置は、を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)のうち1つの解析位置p1において算出して時間軸係数Tn,k,1を求め、該Tn,k,1を前記定Q変換における時間軸係数とする時間軸係数算出部と、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換であるから得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求め、該(1/N)Sn,k,1を前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、前記定Q変換手段は、時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で求められた周波数軸係数(1/N)Sn,k,1と、前記解析位置pm毎の周波数軸係数(1/N)Sn,k,mは前記1つの解析位置p1における周波数軸係数(1/N)Sn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗と、前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施する低周波帯域計算部と、前記解析位置pm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Tn,k,1と時系列データxn-p1+pmとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、を備えていることを特徴とする請求項3に記載の異常診断装置。

請求項7

診断用パラメータ作成手段が、予め異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成ステップと、診断手段が、新たに異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断ステップと、を備えたことを特徴とする異常診断方法。

請求項8

前記診断用パラメータ作成ステップは、第1の定Q変換装置が、時系列データを定Q変換して多変量サンプルを作成するステップと、主成分分析部が、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算するステップと、第1の統計値計算部が、前記主成分得点を基にサンプルの指標となる統計量を得るステップと、第1の正規化部が、前記統計量を正規化して異常度とするステップと、を備え、前記診断ステップは、第2の定Q変換装置が、時系列データを定Q変換して多変量サンプルを作成するステップと、主成分計算部が、前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得るステップと、第2の統計値計算部が、前記主成分得点を基に統計量を得るステップと、第2の正規化部が、前記統計量を、前記第1の正規化部で使用した正規化係数によって正規化して異常度を計算するステップと、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする請求項7に記載の異常診断方法。

技術分野

0001

本発明は、異常診断対象機器の動作時の振動・音響などの波形データ(時系列データ)を周波数解析して異常診断を行う異常診断装置および方法に関する。

背景技術

0002

異常診断対象機器の異常を診断する方法は、例えば特許文献1に記載のものが提案されている。

0003

この特許文献1の方法では、時系列の波形データを、短い時間毎に分割して、フーリエ変換して周波数成分の多変量データサンプルを多数作り、機器が正常と考えられる期間に予め計測した多数の多変量サンプルを主成分分析して主成分得点を求めるための固有ベクトルローテーション行列ともいう)を用意する。

0004

診断の際には、計測データを前記同様にフーリエ変換して周波数成分の多変量データ・サンプルを作り、これを用意しておいたローテーション行列で変換して主成分得点を求め、それを正常時の主成分得点と比較することで異常診断を行っている。要するに診断毎に独立に主成分分析するのではなく、予め主成分分析しておいた結果を利用して、診断時には同じ基準で変換することにより比較を容易にするのが要点である。

0005

尚、本発明で利用する定Q変換の成分演算方法は、例えば非特許文献1〜4に記載されている。

0006

特許第3382240号公報

先行技術

0007

Judith C.Brown:Calculation of a constant Q spectral transform,J.Acoust.Soc.Am.89(1):425−434,1991
Judith C.Brown and Miller S.Puckette:An efficient algorithm for the calculation of a constant Q transform,J.Acoust.Soc.Am.92(5):2698−2701,1992
yukara_13:[Python]Constant−Q変換(対数周波数スペクトログラム)、音楽プログラミングの超入門(仮) in Hatena Blog,2013−12−01,2013、インターネット<URL:http://yukara−13.hatenablog.com/entry/2013/12/01/222742>.[2016−02−22アクセス]
yukara_13:[Python]高速なConstant−Q変換(withFFT)、音楽プログラミングの超入門(仮) in Hatena Blog,2014−01−05,2014、インターネット<URL:http://yukara−13.hatenablog.com/entry/2014/01/05/062414>.[2016−02−22 アクセス]

発明が解決しようとする課題

0008

フーリエ変換を基にした多変量サンプルで予め作成したローテーション行列によって長期にわたって診断しようとすると、事前解析には多くのデータを必要とする。例えば回転機の診断においては、少なくとも低い周波数領域では1Hz程度の周波数解像度が必要とされる。

0009

このためには各サンプルは少なくとも1秒程度の測定が必要であり、時間を重複してサンプル数を増やしても本質的には類似サンプルとなって多様性不足するため、それを基に診断すると、長期的には異常と単なる時間経過による変化の区別がつけられなくなる。

0010

以下に、その事例を述べる。この事例は、ある回転機の振動を長期間(約4ヶ月)にわたって測定して分析したものである。1時間に一回、約100kHzのサンプリング周波数で5秒間の測定を行い、これから時間を重複させながら1秒強の分割データを1測定から100サンプル作成して、各サンプルをフーリエ変換して1/6オクターブバンドの約80系列の多変量サンプルにした。この1週間分(ただし稼動している時刻のデータに絞ったため都合8千サンプルほど)を主成分分析してローテーション行列を計算し、これを基にホテリングT2/Q統計量ベースとした異常度を図6のように算出している。

0011

対象回転機は特に問題のないものであったが、ホテリングT2統計量をベースとしたIndicator1、Q統計量をベースとしたIndicator2とも、図6の「線形」の傾きから分かるように、時間の経過とともに増加していることが伺える。この増加の割合は、学習期間(7日)の10倍(70日)もあれば平均値が異常になるほどであり、長期的な診断にはとても使えない。

0012

この問題は異常度の定義に依存するものではなく、全期間データを主成分分析した図7の第一主成分×第二主成分の散布状況をみても、白色の丸印で示す学習期間のデータが、黒丸印で示す全データの範囲をカバーできていないことは明白である。

0013

尚、図7において、横軸PC1は第一主成分、縦軸PC2は第二主成分を各々示している。

0014

この問題は、時系列データから多変量データを作成するのにフーリエ変換を使っているためにデータの多様性が確保できないためと考えられる。

0015

本発明は上記課題を解決するものであり、その目的は、サンプルの多様性を確保し、長期間にわたって信頼性の高い異常診断を実施することができる異常診断装置および方法を提供することにある。

課題を解決するための手段

0016

上記課題を解決するための請求項1に記載の異常診断装置は、予め異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段と、
新たに異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段と、を備えたことを特徴とする。

0017

また、請求項2に記載の異常診断装置は、請求項1において、前記診断用パラメータ作成手段は、
時系列データを定Q変換して多変量サンプルを作成する第1の定Q変換装置と、
前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部と、
前記主成分得点を基にサンプルの指標となる統計量を得る第1の統計値計算部と、
前記統計量を正規化して異常度とする第1の正規化部と、を備え、
前記診断手段は、
時系列データを定Q変換して多変量サンプルを作成する第2の定Q変換装置と、
前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得る主成分計算部と、
前記主成分得点を基に統計量を得る第2の統計値計算部と、
前記統計量を、前記第1の正規化部で使用した正規化係数によって正規化して異常度を計算する第2の正規化部と、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする。

0018

また、請求項7に記載の異常診断方法は、診断用パラメータ作成手段が、予め異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成ステップと、
診断手段が、新たに異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断ステップと、を備えたことを特徴とする。

0019

また、請求項8に記載の異常診断方法は、請求項7において、前記診断用パラメータ作成ステップは、
第1の定Q変換装置が、時系列データを定Q変換して多変量サンプルを作成するステップと、
主成分分析部が、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算するステップと、
第1の統計値計算部が、前記主成分得点を基にサンプルの指標となる統計量を得るステップと、
第1の正規化部が、前記統計量を正規化して異常度とするステップと、を備え、
前記診断ステップは、
第2の定Q変換装置が、時系列データを定Q変換して多変量サンプルを作成するステップと、
主成分計算部が、前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得るステップと、
第2の統計値計算部が、前記主成分得点を基に統計量を得るステップと、
第2の正規化部が、前記統計量を、前記第1の正規化部で使用した正規化係数によって正規化して異常度を計算するステップと、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする。

0020

上記構成によれば、時系列データから周波数成分の多変量データを作成する際に、フーリエ変換でなく定Q変換を使ったことにより、診断用パラメータ作成手段で使用するデータが比較的少なくてもサンプルの多様性を確保できるようになり、長期間にわたって同じ基準で異常診断が可能になる。これによって異常診断の信頼性が向上する。

0021

また、請求項3に記載の異常診断装置は、請求項2において、前記診断用パラメータ作成手段および診断手段の第1、第2の定Q変換装置は、ともに、
設定した周波数成分KM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Q変換における周波数軸係数との積和計算により求めた定Q変換成分の第1の演算量と、前記KMを超える高周波帯域において、時系列データと定Q変換における時間軸係数との積和計算により求めた定Q変換成分の第2の演算量とを含む合計演算量を算出し、
前記KMを全周波数成分kに各々設定して各KM毎の前記合計演算量を求め、
前記各KM毎の合計演算量を評価し、該演算量が最小となるKMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Q変換における周波数成分kが前記決定された周波数成分閾値KM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第1の計算法を選択し、前記周波数成分kが前記決定された周波数成分閾値KMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第2の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第1の計算法、第2の計算法を実施して定Q変換の成分を求める定Q変換手段と、
を備えたことを特徴とする。

0022

また、請求項4に記載の異常診断装置は、請求項3において、前記各定Q変換装置は、
前記定Q変換における時間軸係数を、

0023

0024

によって算出する時間軸係数算出部と、
前記定Q変換における周波数軸係数を、時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0025

0026

から(1/N)Sn,kとして算出する周波数軸係数算出部と、を備え、
前記定Q変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第1の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、
を備えていることを特徴とする。

0027

上記構成によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Q変換の成分を演算することができ、演算の高速化が実現される。

0028

また、請求項5に記載の異常診断装置は、請求項3において、前記各定Q変換装置は、

0029

0030

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)毎に算出して時間軸係数Tn,k,mを求め、該Tn,k,mを前記定Q変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0031

0032

から得た(1/N)Sn,kを、前記解析位置pm毎に算出して周波数軸係数(1/N)Sn,k,mを求め、該(1/N)Sn,k,mを前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Q変換手段は、
時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、前記解析位置pm毎に、前記周波数軸係数算出部で算出された周波数軸係数(1/N)Sn,k,mと前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施する低周波帯域計算部と、
前記解析位置pm毎に、前記時間軸係数算出部で算出された時間軸係数Tn,k,mと時系列データxnとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、
を備えていることを特徴とする。

0033

上記構成によれば、複数の時刻で各々定Q変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、xn全体に対して一括して離散フーリエ変換したフーリエ係数Xnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。

0034

また、請求項6に記載の異常診断装置は、請求項3において、前記各定Q変換装置は、

0035

0036

を、予め設定した時間方向の複数の解析位置pm(m=1,…,M)のうち1つの解析位置p1において算出して時間軸係数Tn,k,1を求め、該Tn,k,1を前記定Q変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Tn,kをn=0,…,N−1の時系列データとみたときの離散フーリエ変換である

0037

0038

から得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求め、該(1/N)Sn,k,1を前記定Q変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Q変換手段は、
時系列データxn全体に対して離散フーリエ変換を行ってフーリエ係数Xnを求め、 前記解析位置pm毎に、
前記周波数軸係数算出部で求められた周波数軸係数(1/N)Sn,k,1と、
前記解析位置pm毎の周波数軸係数(1/N)Sn,k,mは前記1つの解析位置p1における周波数軸係数(1/N)Sn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗

0039

0040

と、
前記フーリエ係数Xnとを積和計算して前記第1の計算法を実施する低周波帯域計算部と、
前記解析位置pm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Tn,k,1と時系列データxn-p1+pmとを積和計算して前記第2の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Q変換の成分を求める合成部と、
を備えていることを特徴とする。

0041

上記構成によれば、複数の時刻で各々定Q変換を行う際に、1つの解析位置p1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎(複数の解析位置毎)に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。

発明の効果

0042

(1)請求項1〜8に記載の発明によれば、時系列データから周波数成分の多変量データを作成する際に、フーリエ変換でなく定Q変換を使ったことにより、診断用パラメータ作成手段で使用するデータが比較的少なくてもサンプルの多様性を確保できるようになり、長期間にわたって同じ基準で異常診断が可能になる。これによって異常診断の信頼性が向上する。
(2)請求項3、4に記載の発明によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Q変換の成分を演算することができ、演算の高速化が実現される。
(3)請求項5に記載の発明によれば、複数の時刻で各々定Q変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、xn全体に対して一括して離散フーリエ変換したフーリエ係数Xnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。
(4)請求項6に記載の発明によれば、複数の時刻で各々定Q変換を行う際に、1つの解析位置p1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎(複数の解析位置毎)に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。

図面の簡単な説明

0043

本発明の実施形態例による異常診断装置の構成を示すブロック図。
本発明の実施形態例による全期間データの第一主成分×第二主成分散布図
本発明の実施形態例による異常度の期間変化を示す説明図。
本発明の実施形態例による異常診断例を示す説明図。
本発明の実施形態例における定Q変換装置の一例を示す構成図。
従来手法をベースとした異常度の期間変化を示す説明図。
従来のフーリエ変換による全期間データの第一主成分×第二主成分散布図。

0044

以下、図面を参照しながら本発明の実施の形態を説明するが、本発明は下記の実施形態例に限定されるものではない。本発明では、時系列データ(波形データ)から周波数成分の多変量サンプルを生成する際に、従来のフーリエ変換の代わりに定Q変換を使う。

0045

定Q変換(Constant Q transform)は、フーリエ変換のように全ての周波数帯域で同じ期間のデータを解析するのではなく、全ての周波数帯域で同じ周期数になるように、周波数毎に参照するデータ数を変えて解析する。この方法では高周波帯域ほど短い期間のデータで解析するため、低周波帯域での周波数解像度を高く(例えば1Hz程度)するために長い期間(1秒以上)の測定を必要とする際に、期間を重複させ多数のサンプルをとっても高周波帯域では期間が重複しないため、多変量サンプルとしては多様性が確保される。

0046

具体的な診断手法の構成を図1に示す。図1において、本実施例1の手法は、予め測定して集めた時系列データから診断のためのパラメータを計算する準備フェーズと、新たに測定された時系列データを診断する診断フェーズからなる。

0047

準備フェーズは、予め異常診断対象から測定して集めた時系列データを定Q変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段100として構成している。

0048

診断フェーズは、新たに異常診断対象から測定した時系列データを定Q変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段200として構成している。

0049

診断用パラメータ作成手段100は、予め異常診断対象から測定して集めた時系列データを定Q変換して多変量サンプルを作成する第1の定Q変換装置101と、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部102と、前記主成分得点を基にサンプルの指標となる統計量を得る第1の統計値計算部103と、前記統計量を正規化して異常度とする第1の正規化部104と、を備えている。

0050

前記診断手段200は、新たに測定された時系列データを定Q変換して多変量サンプルを作成する第2の定Q変換装置201と、前記作成された多変量サンプルから、前記主成分分析部102で得られた変換行列を使って主成分得点を得る主成分計算部202と、前記主成分得点を基に統計量を得る第2の統計値計算部203と、前記統計量を、前記第1の正規化部104で使用した正規化係数によって正規化して異常度を計算する第2の正規化部204と、を備え、前記異常度に基づいて異常診断対象の異常を診断するものである。

0051

図1の異常診断装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばROM、RAM、CPU、入力装置出力装置通信インターフェースハードディスク記録媒体およびその駆動装置を備えている。

0052

このハードウェアリソースとソフトウェアリソース(OS、アプリケーションなど)との協働の結果、異常診断装置は、図1のように第1の定Q変換装置101、主成分分析部102、第1の統計値計算部103、第1の正規化部104、第2の定Q変換装置201、主成分計算部202、第2の統計値計算部203、第2の正規化部204を実装する。

0053

前記多変量サンプル、主成分得点、統計量、異常度、変換行列、正規化係数などのデータはハードディスクあるいはRAMなどの保存手段・記憶手段に格納されるものとする。

0054

次に、図1の装置の各部の動作を説明する。準備フェーズ(診断用パラメータ作成手段100)では、予め集めた時系列データを、第1の定Q変換装置101において順次定Q変換して多数の多変量サンプルを作成する。定Q変換装置101は、例えば非特許文献1〜4に開示されている方式を実行するものでもよく、また後述する実施例2〜4の高速演算が可能な定Q変換装置を採用しても良い。尚、実施例2〜4に示す方法で高速に求める場合には、時系列データのデータ長はそろっていることが望ましい。

0055

例えば、102.4kHzサンプリングで5秒間のデータが100件など。定Q変換の高周波領域での参照データ数の少なさにより、1回の計測から回転周期との位相電源周期との位相がそれぞれ異なる100〜400程度の多様なサンプルを得ることができる。

0056

次に、主成分分析部102において、多変量サンプルを主成分分析して変換行列を得るとともに、各サンプルの主成分得点を計算する。そして統計値計算部103が、主成分得点を基にホテリングT2/Q統計量のようなサンプルの指標となる統計量を得る。

0057

最後に、算出した統計量が準備フェーズの時系列データで平均0になるように、第1の正規化部104において統計量を正規化して異常度とする。この際に、正規化に使用した正規化係数を保存して診断フェーズ(診断手段200)で使う。尚、異常度の取り得る数値範囲を限定するために、異常度を統計量の対数値にすることが望ましい。

0058

診断フェーズ(診断手段200)では、新たに測定した時系列データを第2の定Q変換装置201において定Q変換して多変量サンプルを作成する。

0059

定Q変換装置201は、定Q変換装置101と同様に、例えば非特許文献1〜4に開示されている方式を実行するものでもよく、また後述する実施例2〜4の高速演算が可能な定Q変換装置を採用しても良い。

0060

次に、多変量サンプルに準備フェーズの主成分分析部102で作成した変換行列を使って、主成分計算部202で主成分得点を得る。そしてこれを基に、第2の統計値計算部203が準備フェーズと同様の計算で統計量を得る。最後に準備フェーズの第1の正規化部104で作成した正規化係数を使って、第2の正規化部204が異常度を計算し、その値により機器が正常か異常かを診断する。

0061

異常診断に使う異常値の閾値を決めるのは一般的に困難であるが、異常値を統計値対数で取る場合には、簡易の診断としては閾値をlog3、log5のような値としても良い。これは振動実効値での異常判断基準が通常時の3倍、5倍などとするのと同様である。

0062

本発明の効果を示すため、図7の第一主成分×第二主成分散布図と同じデータについて、定Q変換による多変量化データを主成分分析した場合の第一主成分×第二主成分散布図を図2に示す。図2によれば、図7の場合と異なり、白色の丸印で示す学習期間のデータが、黒丸印で示す全データの分布範囲を概ねカバーできていることが分かる。

0063

また、図6に示したものと同様の異常度を提案手法により計算した結果を図3に示す。図3によれば、図6の場合と異なり、時間がたっても異常度が上昇傾向にないことが分かる。線形回帰による傾きはほとんどゼロであり、正の傾きを取るIndicator1でも、異常度の平均が閾値に達するのに10年ほど掛かる。

0064

次に、提案手法による異常診断例を図4に示す。図4は、図6図3とは別の回転機の計算例である。こちらも同様に1時間に一回、約100kHzのサンプリング周波数で5秒間の測定を行ったが、測定開始から25日目に故障により緊急停止した。25日目は稼動直後から異常度が非常に大きくなっているが、良く見ると、13日目にIndicator1が注意ラインを越えていることが分かる。

0065

尚、この診断結果はフーリエ変換による手法でもほぼ同様の結果が得られる。

0066

本実施例2は、図1の第1、第2の定Q変換装置101、201を、より精度の高い定Q変換を、少ない演算量で高速に行うことができる図5の定Q変換装置で構成したものである。

0067

まず定Q変換の成分について説明する。

0068

非特許文献1によると、定Q変換は次で定式化される(表記は少し変更している)。

0069

0070

0以外の値は各kに対して、n=(N−Nk)/2,…,(N+Nk)/2−1の範囲でのみ取り、この範囲の窓関数の形状として具体的にはハミング窓ハニング窓、矩形窓など、多くの種類が提案されている
Nk:窓幅,Nk=Q(fs/fk)となるように取る
N≧N1≧…≧Nk…≧NK≧2Qとなる
Q:窓の周期数,定Q変換ではこれを全ての周波数で共通にする
非特許文献1、2ではQ=(21/B−1)-1程度に取るようにしている
Bはビン数で、1オクターブを分割する数のこと
fs:サンプリング周波数
fk:第k成分の周波数,fk=(21/B)k-1fmin
fmin=f1は定Q変換で解析する最小の周波数で、解析する範囲から決める
fkは解析する最大の周波数で、fk≦fs/2となるようにする。

0071

定Q変換は特に低周波帯域の成分で窓幅Nkが非常に大きくなることにより演算量が多い。

0072

これに対して非特許文献2では、高速フーリエ変換(FFT;Fast Fourier transform)を利用した演算の高速化方法が提案されている。

0073

この方法ではパーセバル公式により、定Q変換を次の数式(2)のように時間軸での計算から周波数軸での計算に置き換える。

0074

0075

尚、明細書中の以下の文章では、前記数式(1)、(2)の各左辺の定Q変換の成分を、「Xkcq」と表記することもある。

0076

ここでSn,kは(小さなkに対して特に)ほとんどの要素が0に近いので、一定値以下の項を0として疎行列表現すると(実質非ゼロ項の積和計算だけになるため)非常に高速に計算できるようになる。

0077

さらに(1/N)Sn,kは、計算するデータ長と解析する周波数範囲が決まっていれば、入力データ列xnに寄らず事前に計算できるためこれを予め用意しておけば、入力データの離散フーリエ変換はFFTで実行するとして、定Q変換を入力データに対する1回のFFTと、FFT結果ベクトルと予め用意した疎行列との積演算により計算できることになる。

0078

0079

離散でないフーリエ変換でも同等の公式が成り立ち、これらはパーセバルの定理プランシュレルの定理など様々な呼び方がされる。

0080

非特許文献2による高速化のアイデアは周波数軸での係数Sn,kが疎行列とみなせることに依拠しているが、実際にはSn,kが疎行列というのは粗い評価である。確かにその値は一部を除いて小さいが、精度を高めると必ずしも疎にならない。特に高周波帯域でこの傾向が顕著である。実際、非特許文献2、4で係数Sn,kが十分小さいとしている閾値はそれぞれ0.15、0.0054であり、通常の評価には問題ないが有効数字7桁以上の精度を得るには十分とは言えない。このため、要求精度によってはこの高速化手法がかえって遅くなる可能性がある。

0081

実際に周波数軸での計算でどの程度の演算量を必要とするかを調べるため、非特許文献4に準じた設定により時間軸での計算と周波数軸での計算でそれぞれ計算対象となる非ゼロの項数を評価する。

0082

設定パラメータは、本明細書の記号法で表記してB=24、Q=28、fs=16000、fmin=60、K=160とし、定Q変換をそれぞれ矩形窓、ハミング窓、ハニング窓で計算する場合において、時間軸での計算に必要な非ゼロ項数と、周波数軸での計算で所定の精度の係数まで計算する際に必要な項数を表1に示す。

0083

0084

なお、非特許文献1、2に合わせるとB=24のときQ=34となるが、非特許文献3、4では追加の係数fratioを導入してQの値を変えているので、その設定値Q=28を使っている。また解析する最大周波数6000HzからK=160を導いた。fsの値は非特許文献3,4に明記されていなかったので、最大周波数6000Hzまで解析できるサンプリング周波数12000Hz以上のWAV音声データとしてfs=16000Hzを選択した。

0085

表1中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の項数を基準としている。

0086

表1から、窓形状により程度が違うが、精度が高いほど(係数をゼロと評価する係数Sn,kの大きさの閾値が小さいほど)周波数軸での計算に必要な項数が増えることが分かる。特に矩形窓やハミング窓では精度が高くなると時間軸での計算より周波数軸での計算の方が積和すべき項数が大きくなり、周波数軸での計算で高速にならなくなる。

0087

尚、周波数軸での計算にはFFT計算と疎行列の非ゼロ係数抽出のオーバーヘッドもあるので、必要項数に差がなければ時間軸計算が有利である。

0088

より精度の高い定Q変換を高速に行うには工夫が必要になる。

0089

そこで本実施例2では、定Q変換を高精度でも演算量を少なく計算するために、定Q変換の成分Xkcqを数式(2)の時間軸での計算と周波数軸での計算のどちらの方法で計算するか、予め周波数成分k毎に演算量の小さい方を選択しておくことで実現する。

0090

定Q変換の成分Xkcqの計算はおおよそ低周波帯域では周波数軸での計算の演算量が小さく、高周波帯域では時間軸での計算の演算量が小さいから、kの境界値KMを決めて、k≦KMでは周波数軸での計算、k>KMでは時間軸での計算、というようにすれば良い。

0091

図5において、10は予め定Q変換を行うためのパラメータを決定して定Q変換計算の準備をする準備部であり、20は個別の時系列データに対する定Q変換を計算する個別計算部である。

0092

準備部10は、
設定した周波数成分KM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Q変換における周波数軸係数との積和計算により求めた定Q変換成分の第1の演算量と、前記KMを超える高周波帯域において、時系列データと定Q変換における時間軸係数との積和計算により求めた定Q変換成分の第2の演算量とを含む合計演算量を算出し、前記KMを全周波数成分kに各々設定して各KM毎の前記合計演算量を求め、前記各KM毎の合計演算量を評価し、該演算量が最小となるKMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段の機能と、
定Q変換における周波数成分kが前記決定された周波数成分閾値KM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第1の計算法を選択し、前記周波数成分kが前記決定された周波数成分閾値KMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第2の計算法を選択する計算法選択手段の機能とを具備している。

0093

前記計算法境界決定手段および計算法選択手段の各機能を実行する準備部10は、具体的には、定Q変換のパラメータを決定するパラメータ決定部11、時間軸係数Tn,kを計算する時間軸係数算出部12、周波数軸係数(1/N)Sn,kを計算する周波数軸係数算出部13、時間軸での計算と周波数軸での計算との境界のk値KMを決定する計算法境界決定部14からなる。

0094

個別計算部20は、時系列データに対して、前記計算法選択手段により選択された第1の計算法、第2の計算法を実施して定Q変換の成分を求める定Q変換手段の機能を具備している。

0095

前記定Q変換手段の機能を実行する個別計算部20は、具体的には、時系列データの長さを整える前処理部21、高周波帯域の定Q変換を計算する高周波帯域計算部22、低周波帯域の定Q変換を計算する低周波帯域計算部23、高周波帯域での結果と低周波帯域での結果をまとめる定Q変換合成部24からなる。

0096

図5の定Q変換の成分演算装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばROM、RAM、CPU、入力装置、出力装置、通信インターフェース、ハードディスク、記録媒体およびその駆動装置を備えている。

0097

このハードウェアリソースとソフトウェアリソース(OS、アプリケーションなど)との協働の結果、定Q変換の成分演算装置は、図5に示すように、パラメータ決定部11、時間軸係数算出部12、周波数軸係数算出部13、計算法境界決定部14、前処理部21、高周波帯域計算部22、低周波帯域計算部23、定Q変換合成部24を実装する。

0098

また各部での処理結果は、ハードディスク或いはRAMなどの保存手段・記憶手段に格納され、随時利用されるものとする。

0099

まず、準備部10のパラメータ決定部11では、予め定Q変換で解析する時系列データのサンプリング周波数fs、解析する周波数範囲(fmin,fmax)と1オクターブの周波数を区切るビン数B、解析対象データの周期数Q、窓関数を決定する。このときfmax≦fs/2が必要で、K=floor(B*log2(fmax/fmin))+1となる。これを基に計算する入力データの長さをN=N1=Q(fs/fmin)とする。

0100

また周波数軸係数をゼロとみなすべき係数の大きさの閾値Thもここで決定しておく。

0101

時間軸係数算出部12では、パラメータ決定部11で決定したパラメータを基に、時間軸係数Tn,kを次の定義式に従って計算する。

0102

0103

周波数軸係数算出部13では、パラメータ決定部11で決定したパラメータを基に、周波数軸係数(1/N)Sn,kを、

0104

0105

の定義式に従って計算する。この際、|Sn,k|≦Thとなる項は周波数軸係数(1/N)Sn,k=0とする。

0106

計算法境界決定部14では、周波数軸での計算を行う最大のk値KM(周波数成分閾値)を求める。KMは0,…,Kの何れかの値を持ち、KM=0は周波数軸での計算を行わずに非特許文献1、2と同様に全て時間軸での計算で求める場合を意味する。

0107

周波数成分閾値KMの決定方法としては、KM=0,…,Kの各場合について、定Q変換の成分Xkcqを求めるのに必要な演算量CX(KM)を評価してそれが最小となるようにKMを決定する。

0108

定Q変換の成分Xkcqの演算量CX(KM)(合計演算量)は、長さNのFFTの演算量CFFT(N)、非ゼロの時間軸係数毎の演算量CT、非ゼロの時間軸係数の数NT(k)、非ゼロの周波数軸係数毎の演算量CF、kで非ゼロとみなす周波数軸係数の数NF(k)を使っておおよそ

0109

0110

と表せる。

0111

ここでCFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、CFFT(N)=(1/4)Nlog2(N)、CT=1、CF=1としても良い。

0112

上記演算量CX(KM)を示す数式(3)において、右辺第1項は、KMが非ゼロのときCFFTとし、そうでないときは0とすることを表し、右辺第2項は、周波数成分kが周波数成分閾値KM以下の低周波帯域における周波数軸係数毎の演算量を表し、右辺第3項は、周波数成分kが周波数成分閾値KMを超える高周波帯域における時間軸係数毎の演算量を表している。

0113

このためKMが非ゼロのときは、数式(3)の右辺第1項のFFTの演算量および右辺第2項の演算量(周波数軸係数による演算量)からなる第1の演算量と、数式(3)の右辺第3項の時間軸係数による第2の演算量との合計が演算量CX(KM)となる。

0114

またKMがゼロのときは、数式(3)の右辺第1項および右辺第2項がともに0となって右辺第3項の時間軸係数による第2の演算量のみが演算量CX(KM)となる。

0115

この演算量CX(KM)は、非ゼロの時間軸係数の数NT(k)がkに関して指数的に減少し、非ゼロの周波数軸係数の数NF(k)がkに関しておおよそ指数的に増加するため、KM=1,…,Kの範囲ではおおよそ下に凸の特性カーブとなる。ただし、KM=0ではCFFTがなくなるためCX(0)が小さくなる。結果として演算量CX(KM)はKM=0あるいはKM=1,…,Kでの最小点の何れかで最小値を取ることになる。

0116

個別計算部20では、入力された時系列データの定Q変換を次のようにして求める。

0117

前処理部21では、入力された時系列データの長さがNより短ければ前後に0を補完して、Nより長ければ前後の値を除去して長さがNの時系列データxnになるようにする。このとき入力された時系列データの中央部分が時系列データxnの中央部分に一致するようにする。

0118

次に定Q変換の成分Xkcqを、前記計算法境界決定部14で決定した周波数成分閾値KMにより、k≦KMについては低周波帯域計算部23で計算し、k>KMについては高周波帯域計算部22で計算する。

0119

高周波帯域計算部22は、前記時間軸係数算出部12で算出された時間軸係数Tn,kと時系列データxnを積和計算して定Q変換の成分Xkcqを求める。各kについて非ゼロの時間軸係数Tn,kは高々Nk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

0120

低周波帯域計算部23は、まず時系列データxnをFFTしてフーリエ係数Xnを求める。次に前記周波数軸係数算出部13で算出された周波数軸係数(1/N)Sn,kとフーリエ係数Xnを積和計算して定Q変換の成分Xkcqを求める。各kについて非ゼロの周波数軸係数(1/N)Sn,kは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

0121

尚、周波数軸係数は多くの場合、両端(周波数ゼロ近辺とサンプリング周波数近辺)部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。

0122

最後に定Q変換合成部24で、高周波帯域計算部22と低周波帯域計算部23とで計算した定Q変換の成分Xkcqをまとめて、定Q変換を完了する。

0123

上記のように本実施例2によれば、従来の高速化技術の問題点である、高い精度の定Q変換を計算しようとすると演算量が膨れ上り通常の時間軸での計算よりもむしろ遅くなるという問題が改善される。

0124

その効果を示すため、定Q変換の計算方法(時間軸計算、周波数軸計算、提案手法)毎に、必要な演算量を表2のように試算した。試算は、時系列データのパラメータを表1と同じにして、演算量パラメータをCFFT(N)=(1/4)Nlog2(N)、CT=1、CF=1として計算した定Q変換の成分Xkcqの演算量CX(KM)を表2に示す。

0125

0126

表2中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の演算量を基準としている。

0127

表2では、FFT演算の演算量が原因で高速化の効果が小さめに出ているが、ハニング窓では閾値1.0e−5でもおおよそ1/5(22.5%)の演算量になっている。またFFT演算の演算量が重く周波数軸での計算で効果が出ない場合でも、自動的に時間軸での計算になるので反って遅くなるという問題は発生しない。

0128

実施例2は、ある一時刻における定Q変換の各成分値を求める際には効率的であるが、実際の定Q変換では、時系列データの周波数成分の時間変化を捉えるために少しずつ時刻をずらしながら繰り返し定Q成分を計算する。このときずらす時刻は、最短ではサンプリング周波数で、時系列データの情報を余すことなく計算するには最長でも最高周波数での窓幅にしたい。もちろん、時系列データの特徴サンプル抽出などでは、これよりも長い間隔で時刻をずらして定Q変換する場合もあるが、いずれにせよ時系列データ全体に対しては多くの繰り返し計算をすることになり、まだまだ負荷が高い。

0129

多数の時刻で同時に定Q変換を計算する際には更なる効率化が必要である。

0130

そこで本実施例3では、数式(2)の周波数軸での計算で定Q変換を行う場合に、定Q変換を計算する時刻ごとに計算に使用する時系列データの範囲でDFT(あるいはFFT)する代わりに、複数時刻での定Q変換に必要な範囲全体で一括してDFT(あるいはFFT)を行うことで高速化を図った。

0131

数式(1)で定義する定Q変換はデータ列xn,n=0,…,N−1の中央xN/2での成分Xkcqを計算するように定式化しているが、窓関数の非ゼロ範囲を時刻方向に平行移動することで他の時刻での成分を計算するように調整できる。すなわちTn,kとxnとの積和のインデックスをずらすことで実現できる。

0132

数式(2)による周波数軸での計算に関しては、初めに定Q変換する時刻全体で必要な全データに対して一括してDFT(あるいはFFT)を行い、予め各時刻用に用意した(1/N)Sn,kと順次積和することで実現する。

0133

本実施例3による定Q変換の成分演算装置の具体的な計算の構成は、図5と同様であるが、各部で実行される機能が以下に示すように異なる。

0134

まず、準備部10のパラメータ決定部11では、予め定Q変換で解析する時系列データのサンプリング周波数fs、解析する周波数範囲(fmin,fmax)と1オクターブの周波数を区切るビン数B、解析対象データの周期数Q、窓関数を決定する。このときfmax≦fs/2が必要で、K=floor(B*log2(fmax/fmin))+1となる。

0135

さらに入力データxnの全長Nと、解析位置pm,m=1,…,M(時間方向の解析位置;時刻)を決定しておく。このとき、各位置pmでの定Q変換はxpmを中心とする長さN1=Q(fs/fmin)のデータ列をもとに計算し、(N1/2)≦p1≦…≦pM≦N−(N1/2)であるものとする。すなわち、全ての解析位置pm,m=1,…,Mで定Q変換の計算に必要なデータがx0,…,xN-1に含まれるものとする。尚、入力データの全点で定Q変換を計算したいような場合は予め入力データを前後にゼロ拡張して上記条件を満たすように変換しておけばよい。

0136

また、周波数軸係数をゼロとみなすべき係数の大きさの閾値Thもここで決定しておく。

0137

時間軸係数算出部12ではパラメータ決定部11で決定したパラメータを基に、時間軸係数Tn,kを各解析位置pmについて次の定義式に従って計算する。

0138

0139

このときpmに対するTn,kをTn,k,mと書くと、非ゼロ項に関してTn,k,m=Tn-pm+p1,k,1である。

0140

周波数軸係数算出部13ではパラメータ決定部11で決定したパラメータを基に、各解析位置pmについて周波数軸係数(1/N)Sn,k,mを定義式に従って計算する。この際、|Sn,k,m|≦Thとなる項は周波数軸係数(1/N)Sn,k,m=0とする。

0141

計算法境界決定部14では、周波数軸での計算を行う最大のk値KM(周波数成分閾値)を求める。KMは0,…,Kの何れかの値を持ち、KM=0は周波数軸での計算を行わずに非特許文献1、2と同様に全て時間軸での計算で求める場合を意味する。

0142

周波数成分閾値KMの決定方法としては、KM=0,…,Kの各場合について、定Q変換の成分Xkcqを求めるのに必要な演算量CX(KM)を評価してそれが最小となるようにKMを決定する。尚、この条件は複数のpmに寄らないので、p1(1つの解析位置)について求めておけば十分である。

0143

定Q変換の成分Xkcqの演算量CX(KM)(合計演算量)は、長さNのFFTの演算量CFFT(N)、非ゼロの時間軸係数毎の演算量CT、非ゼロの時間軸係数の数NT(k)、非ゼロの周波数軸係数毎の演算量CF、kで非ゼロとみなす周波数軸係数の数NF(k)を使っておおよそ、前記実施例2で述べた数式(3)

0144

0145

と表せる。

0146

ここでCFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、CFFT(N)=(1/4)Nlog2(N)、CT=1、CF=1としても良い。

0147

この演算量CX(KM)は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数(周波数成分k;対数周波数)の関係が、非ゼロの時間軸係数の数NT(k)がkに関して指数的に減少し、非ゼロの周波数軸係数の数NF(k)がkに関しておおよそ指数的に増加するため、KM=1,…,Kの範囲ではおおよそ下に凸の特性カーブとなる。ただし、KM=0ではCFFTがなくなるためCX(0)が小さくなる。結果として演算量CX(KM)はKM=0あるいはKM=1,…,Kでの最小点の何れかで最小値を取ることになる。

0148

個別計算部20では入力された時系列の定Q変換を次のようにして求める。

0149

前処理部21では、入力された時系列データの長さがNより短ければ前後にゼロを補完して、Nより長ければ前後の値を除去して長さがNの時系列データxnになるようにする。ゼロを補完する場合には、入力された時系列データの中央部分が時系列データxnの中央部分に一致するようにする。データを除去する場合は、解析したい範囲に寄るが、単純に後ろのデータを除去するのでも良い。

0150

次に定Q変換の成分Xkcqを、計算法境界決定部14で決定したKMにより、k≦KMについては低周波帯域計算部23で、k>KMについては高周波帯域計算部22で計算する。

0151

0152

尚、前記定Q変換の成分は以下、「Xk,mcq」と表記することもある。

0153

各kについて、前記時間軸係数算出部12で算出された非ゼロの時間軸係数Tn,k,mは高々Nk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

0154

低周波帯域計算部23は、まず時系列データxnを全体に対して一度FFTしてフーリエ係数Xnを求める。次にpm毎に、前記周波数軸係数算出部13で算出された周波数軸係数(1/N)Sn,k,mとフーリエ係数Xnを積和計算して定Q変換の成分Xk,mcqを求める。

0155

各kについて非ゼロの周波数軸係数(1/N)Sn,k,mは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。

0156

尚、周波数軸係数は多くの場合、両端(周波数ゼロ近辺とサンプリング周波数近辺)部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。

0157

最後に定Q変換合成部24で、高周波帯域計算部22と低周波帯域計算部23とで計算した定Q変換の成分Xk,mcqをまとめて、定Q変換を完了する。

0158

前記実施例3では、時間軸係数Tn,k,m、周波数軸係数(1/N)Sn,k,mは解析位置pm毎に用意するので、事前に用意して保管しておくデータ量が大きくなる。本実施例4では、係数データのメモリ使用量を大幅に減らすように以下のような対策を講じた。

0159

0160

尚、数式(4)中の

0161

0162

は、複素指数関数の時刻および周波数毎に異なるべき乗を表している。

0163

本実施例4による定Q変換の成分演算装置の具体的な構成は図5と同様であり、以下に図5と異なる部分のみを説明する。

0164

時間軸係数算出部12では、パラメータ決定部11で決定したパラメータを基に、時間軸係数Tn,k,1を定義式に従って求める。すなわち、

0165

0166

を、予め設定した解析位置pmのうち1つの解析位置(例えば1番目の解析位置)p1において算出し、時間軸係数Tn,k,1を求める。

0167

周波数軸係数算出部13では、パラメータ決定部11で決定したパラメータを基に、周波数軸係数(1/N)Sn,k,1を定義式に従って計算する。

0168

すなわち、

0169

0170

から得た(1/N)Sn,kを、前記1つの解析位置p1において算出して周波数軸係数(1/N)Sn,k,1を求める。

0171

この際、|Sn,k,1|≦Thとなる項は周波数軸係数(1/N)Sn,k,1=0とする。

0172

計算法境界決定部14では、周波数軸計算の際に複素べき乗の積和が追加されるために非ゼロの周波数軸係数毎の演算量CFを実施例3の場合より大きく、例えば5(複素べき乗の演算に+3、複素べき乗の積和に+1)とする。

0173

高周波帯域計算部22は、複数の解析位置pm毎に範囲をずらして、前記時間軸係数算出部12で算出された時間軸係数Tn,k,1と時系列データxn-p1+pmを積和計算して定Q変換の成分Xk,mcqを求める。

0174

低周波帯域計算部23は、まず時系列データxnを全体に対して一度FFTしてフーリエ係数Xnを求める。

0175

0176

前記実施例3、4の手法によれば、定Q変換を高速に実行するにあたって、先行の高速化技術(実施例2)が時間軸での計算でも周波数軸での計算でもそれらを併用する計算方法でも、一時刻毎に定Q変換を計算するのに対して、複数時刻での定Q変換を同時に計算する際に、DFT(あるいはFFT)を一時刻の計算範囲毎ではなく全体に対して一度だけ実施することで、全体としての計算が高速化される。

実施例

0177

実施例4では、時間軸と周波数軸との係数データを計算時刻毎に持つ代わりに、それら複素数のべき乗計算に置き換えて(数式(4))、演算量は少し増える代わりに係数データのメモリ使用量を大幅に減らすことができる。

0178

10…準備部
11…パラメータ決定部
12…時間軸係数算出部
13…周波数軸係数算出部
14…計算法境界決定部
20…個別計算部
21…前処理部
22…高周波帯域計算部
23…低周波帯域計算部
24…定Q変換合成部
100…診断用パラメータ作成手段
101…第1の定Q変換装置
102…主成分分析部
103…第1の統計値計算部
104…第1の正規化部
200…診断手段
201…第2の定Q変換装置
202…主成分計算部
203…第2の統計値計算部
204…第2の正規化部

ページトップへ

この技術を出願した法人

この技術を発明した人物

ページトップへ

関連する挑戦したい社会課題

関連する公募課題

ページトップへ

技術視点だけで見ていませんか?

この技術の活用可能性がある分野

分野別動向を把握したい方- 事業化視点で見る -

ページトップへ

おススメ サービス

おススメ astavisionコンテンツ

新着 最近 公開された関連が強い技術

  • 三菱電機エンジニアリング株式会社の「 負荷試験装置」が 公開されました。( 2020/10/29)

    【課題】供試体の入力軸および出力軸にトルク制御を適用する際に、回転速度の制御が発散してしまうことを抑制する。【解決手段】負荷試験装置は、入力軸用モータ(SM1)と第1トルクメータ(TM1)と第1エンコ... 詳細

  • ミサワホーム株式会社の「 計測装置の取付構造」が 公開されました。( 2020/10/29)

    【課題】人や物が計測装置に接触しても正確な計測ができるようにする。【解決手段】振動を検出する装置本体11と、装置本体11を覆って装置本体11を周囲から遮蔽するケース12と、を備えた計測装置10の取付構... 詳細

  • 信藤薫の「 被対象物の状態モニタリングシステム及び状態モニタリング方法」が 公開されました。( 2020/10/29)

    【課題】初期コスト、メンテナンスコスト及び人的コスト等のコストを削減し、かつ河川や傾斜地といった国土表層、橋梁、隧道、建物等の被対象物の状態把握の精度を飛躍的に向上することができる被対象物の状態モニタ... 詳細

この 技術と関連性が強い人物

関連性が強い人物一覧

この 技術と関連する社会課題

関連する挑戦したい社会課題一覧

この 技術と関連する公募課題

関連する公募課題一覧

astavision 新着記事

サイト情報について

本サービスは、国が公開している情報(公開特許公報、特許整理標準化データ等)を元に構成されています。出典元のデータには一部間違いやノイズがあり、情報の正確さについては保証致しかねます。また一時的に、各データの収録範囲や更新周期によって、一部の情報が正しく表示されないことがございます。当サイトの情報を元にした諸問題、不利益等について当方は何ら責任を負いかねることを予めご承知おきのほど宜しくお願い申し上げます。

主たる情報の出典

特許情報…特許整理標準化データ(XML編)、公開特許公報、特許公報、審決公報、Patent Map Guidance System データ