明細書 状態判定方法と状態予測方法及び装置
一. 技術分野
本発明は設備診断、 医療診断などにおいて、 対象物の状態変化の有無を判定す るための状態判定装置及びオンライン状態監視 ·診断システムに関するものであ る。 二. 背景技術
( 1) 従来、 状態監視の対象物から測定した波形データの確率密度関数は、 正 常状態時に 「正規分布」 に従うと仮定して状態判定を行う ([ 1 ]特許公開 2000 一 171291、 [2] 豊田利夫、 二保知也: 「正常時の振動波形のみを用いた回転機 械の異常診断」、 日本設備管理学会誌、 Vol. 11, No. 1, 1999年、 p.4- 11. [3] 劉 信芳、 豊田利夫、 陳 鵬、 馮 芳、 二保知也: information Divergenceに よる回転機械の異常診断」、 精密工学会誌、 Vol.66、 No.1, 2000年、 ρ· I57- 1δ2· )
(2) 対象物の状態変化の有無を判定するために用いる特徴パラメータを正規 分布の確率変数と仮定して統計検定を行う ([4] 河部佳樹、 豊田利夫、 江口 透、福田和久: 「無次元兆候パラメータを用いた回転機械の劣化傾向管理(ΠΙ)」、 平成 7年度日本設備管理学会秋季研究発表大会論文集、 Vol. 2、 1995年、 p.32. ) 状態監視の対象物から測定した信号が脈動信号の場合は、 特徴パラメータや確 率密度分布関数などにより状態判定が困難である。 従来、 脈動信号の中に含まれ る特異成分 (あるいは、 異常成分) を検出する方法として、
(3) 時間一周波数解析手法 (ウェーブレット解析、 ウイダナ分布解析、 短時 間 FFTなど) を用いた局所異常の検出法がある ([ 5] 章 忠、 中堀智之、 川畑 洋昭: 「高速ウェーブレツト変換およびその脳波解析への応用」、 日本機械学会 論文集 (C編)、 Vol.65, No.632, 1999年、 p. 1915-1921.)
(4) 脈動信号のピーク値などに合わせた所定の間隔でサンプリングした値は 所定の値以下、 または、 前回サンプリングした値との差が所定の値以上の場合に 異常と判定する。 ([6] 特開平 6-317215)
( 5 ) 脈動信号のピーク値などに合わせた所定の間隔でサンプリング毎の代表 値のばらつき (標準偏差) が所定のしきい値以上の場合に異常と判定する。 ([ 7 ] 特開平 7- 293311) 三. 発明が解決しょうとする課題
しかしながら、 前記従来の諸手法には次のような問題があった。
第 1と第 2の方法では、 波形データと特徴パラメータの値が正規分布に従うと 仮定しているが、実際に正常状態でも必ずしも正規分布に従うとは限らないから、 波形データと特徴パラメータの値が正規分布に従うと仮定して対象物の状態変化 を統計検定により判定すると、 誤判定をもたらす。
第 3の方法では、 処理に必要な時間が長いため、 アルタイムの特異点の検出が 困難である。
第 4と第 5の方法では、 波形データを所定の間隔でサンプリングする場合、 回 転数の変化や負荷の変化などによるピーク位置が変化した時、 サンプリング値は 正常でも大きく変化するから、 誤判定をもたらす。 また、 サンプリング毎の代表 値のばらつき (標準偏差) を計算すると、 リアルタイムの特異検出が困難である。 四. 課題を解決するための手段
上記に述べたような問題点を解決するために、 本発明においては、 測定した波 形データを既知確率分布 (例えば、 正規分布) の波形データに変換した後、 ある いは、 波形データから算出された特徴パラメータを既知確率分布 (例えば、 正規 分布) の特徴パラメータに変換した後、 統計検定や可能性理論や情報理論などに より対象物の状態判定を行う。
また、 測定した対象物の信号が脈動信号である場合は、 特徴パラメータゃスぺ ク トルなどによる特異成分の抽出及び状態判定が困難であるので、 雑音を除去し た信号の包絡線波形データを求め、 正規化された包絡線波形データおよび周期パ ルスにより信号中の特異成分を検出して状態判定を行う。 または、 上記に雑音を 除去した脈動信号のピーク波形データを求め、 ピーク波形データの正規化処理を 行い、 正規化されたピーク波形データおよび周期パルス波形データにより脈動信 号中の特異成分を検出して状態判定を行う。
五. 発明の効果
対象物の状態監視のために測定された波形データや特徴パラメータは必ずしも 正規分布に従うと限らないため、 これらの波形データあるいは特徴パラメータを 正規分布に従うと仮定して対象物の状態変化の有無を統計検定により判定し、 或 いは、 状態予測を行うと、 大きな誤差をもたらす。 本発明においては、 測定した 波形データを既知確率分布 (例えば、 正規分布) の波形データに変換した後、 あ るいは、 波形データから算出された特徴パラメータを既知確率分布 (例えば、 正 規分布) の特徴パラメータに変換した後、 統計検定や可能性理論や情報理論など により対象物の状態判定を行う。 従って、 本発明の手法は状態判定、 或いは、 状 態予測の精度が従来正規分布に従うと仮定した場合より高い。
脈動信号の場合、 測定した脈動信号の波形データには一部のピークが欠落し、 または波形データの局所に微小な特異成分が存在する時、 従来の判定法 (特徴パ ラメータやスペク トルや確率密度関数など) による検出が困難である。 本発明に おいては、 雑音を除去した脈動信号の包絡線波形データを求め、 正規化された包 絡線波形データおよび周期パルスにより信号中の特異成分を検出して状態判定を 行う。 または、 上記に雑音を除去した脈動信号のピーク波形データを求め、 ピー ク波形データの正規化処理を行い、 正規化されたピーク波形データおよび周期パ ルス波形データにより脈動信号中の特異成分を検出して状態判定を行う。 本発明 の方法では、主な信号処理はハードウエアでも実現でき、数値処理装置(計算機) にかかる負担が小さいから、 リアルタイムに異常の検出が実現できる。 また、 正 常状態の波形データを基準としないため、 脈動周期が変化しても検出と判定の結 果に与える影響が少なレ、。 六. 発明を実施するための最良の形態
1 . 特徴パラメータについて
状態判定のために用いる特徴パラメータは時間領域と周波数領域の特徴パラメ ータがある。 周波数領域の特徴パラメータは (参考文献 1 ) で定義されているも のがある。 ([ 8 ] 陳 鵬, 豊田利夫:遺伝的プログラミングによる周波数領域 の特徴パラメータの自己再組織化, 日本機械論文集 (C編), Vol. 65 No. 633, p
p.1946-1953, 1998. ) ここで時間領域の特徴パラメータについて詳述する。
対象物の状態変化の有無を判定するために用いる時間領域の特徴パラメータは 次の通りである。
1) 無次元特徴パラメータ
測定した時系列波形データから、 フィルタを用いて、 低、 中、 高周波数領域の 波形データを抽出する。 抽出した波形データ X ( t ) を次式で正規化する。
X;
ここで、 χ' iは A/D変換後の X (t)の離散値であり、 と Sはそれぞれ x' iの平 均値と標準偏差である。
従来用いられる無次元特徴パラメータは式 (2) 〜式 (1 3) に示す ([9] Peng CHEN, Toshio TOYOTA, Yueton LIN, Feiyue WANG : FAILURE DIAGNOSIS OF MACHINERY BY SELF-REORGANIZATION OF SYMPTOM PARAMETERS IN TIME DOMAIN U SING GENETIC ALGORITHMS, International Journal of Intelligent Control an d System, Vol.3, No.4, pp.571—585, 1999. )
変動率: ρ^σ/μ^ (2) ここで、 μ"¾=》 - |λ-,.|/Ν 3 ) は絶対平均値で, Νはデータの総数である,
1 (4)
N-1 は標準偏差である。 -l (5) 2= (Ν-1)σ3
尖度: ' (6)
(Ν-Ι)σ4
4=μノ μ ¾ ( 7 ) こで、 μΡは波形の極大値 (ピーク値) の平均値である。
5=|μ, |/μρ ( 8 ) ここで、 | Iは波形の 1 0個の最大値の平均値である。
ここで、 σ Pは極大値の標準偏差値である。
Pl=^L
/0L ( 1 0) ここで、 Lと σ まそれぞれ極小値 (谷値) の平均値と標準偏差値である。
N
Exf
p9=― ( 1 2 )
式 (2) 〜式 (1 3) は従来の特徴パラメータであるが、 数値計算で容易に高 速計算するために、 「区間特徴パラメータ」 は式 (1 4) 〜式 (2 1 ) のように 新たに提案する。
'·=ι ( 1 4)
但し、
X i≥ k a kは任意に設定できるが、 例えば、 k =0.5 1 2。
klは
X iの平均値である。 tは任意に設定できるが、 例えば、 t =0.5 1 2 3 4
h\
N/,i
但し、 x
h — 1ι σ、 hは任意に設定できるが、 例えば、 k =0.5、 1、 2。 μ
hl は x
hの平均値である。 tは任意に設定できるが、 例えば、 t =0.5、 1、 2、 3、 4。
ho
15 hP ( 1 8 ) 但し、 h。は時系列波形が単位時間あたり 0 レベルをクロースする頻度、 hPは単 位時間あたり時系列波形のピークの数である。
Pl ( 1 9 )
P ( 2 0 ) 但し、 h„„は波形が単位時間あたり ησ レベルをクロースする頻度である。 ηは 任意に設定できるが、 例えば、 η =0.5、 1、 2。
-na ( 2 1 ) 18 ho 但し、 h-n。は波形が単位時間あたり- ησ レベルをクロースする頻度。 である。 ηは任意に設定できるが、 例えば、 η =0.5、 1、 2。
2 ) 有次元特徴パラメータ
有次元特徴パラメータを計算するとき、測定した波形データに対して、式(1 ) のような正規化をしない。
波形データの絶対平均値: ( 2 2)
Ν 波形データの実効値: —^ ( 2 3)
^ Γ
波形データの絶対値のピーク平均値: ζ^ 1"^ ( 2 4 ) ここで、 | xiしは波形データの絶対値のピーク値 (極大値)、 NPはピーク値の総 数である。 波形データの絶対値のピーク実効値: ( 2 5 )
Pd4= ~" なお、 上記の特徴パラメータ以外にも多くの特徴パラメータが定義できるが、 本方法を応用するとき、 まず上記の特徴パラメータで試し、 もし状態識別の効果 が良くなければ、 更に他の特徴パラメータを追加定義すればよい。
2 . 波形データと特徴パラメータを指定の確率分布への変換
測定した波形データを で、 波形データから算出した特徴パラメータを で表 す。 X*iあるいは
P*iを用いて、 統計理論により状態判定や状態予測を行う場合、 A あるいは ϊΛがどのような確率分布に従うか事前に知る必要がある。 しカゝし、 χ あ るいは P*iがどのような確率分布従うかは事前に知らない場合が多い。 そこで、 指 定の既知確率分布関数を∑とすると、 次式を用いて x*iあるいは ιΛを∑に従う確率 変数 Xiあるいは Piに変換することができる。
ここで、 Fxi (x*i)と FPi ( )はそれぞれ x*iと ιΛとの累積確率分布 (或いは、 累積度 数分布) であり、 ∑ は∑の逆関数である。 例えば、 ∑は正規分布、 ワイブル分 布、 指数分布、 ガンマ分布などである。 X あるいは を あるいは Piに変換した 後、 統計検定などを用いて状態判定や状態予測を行う。
なお、 元の波形データ x*iは図 1のように 4種類に分けられる。 つまり、 平均値 より大きいデータ Xi+、 平均値より小さいデータ -、 式 (1 ) による正規化後の絶 対値データ i xi |及び全体波形データ χ*ίΛであり、それぞれの累積確率分布(或いは、 累積度数分布) は FX (X i+)、 Fxi— (Xi -)、 ?^ ( | | )及ぴ1^ )で表すが、 以下、 特 に指定しなければ、 統一に Fxi (x*i)で表す。 また、 特徴パラメータ p*iは正規化後の
Xi+と Xi -、 I Xi I及び X*iAのいずれかで計算される。
ここで、 例として、 ∑を正規分布に指定した場合、 X*iあるいは p*iを正規分布 に従うように変換する方法について詳述する。
ここで、 は確率変数 tの平均値で、 σは標準偏差である。
( 1 ) 基準状態を参照する場合
測定対象の基準状態、 例えば、 1回目測定した時の状態、 を決めて、 この時の 波形データ X*i。あるいは特徴パラメータ 。の確率密度関数と累積確率分布関数を それぞれ fx。(x*i。)と Fx。(x*i。)あるいは fP。(p*i。)と FP。(p*i。)とする。 平均値 μが 0で標準 偏差 σが 1である標準正規分布の確率密度関数を φ (χ 、 標準正規分布の確率分 布関数を Φ 0 )とする。 なお、 離散データである X あるいは p*i。の確率密度関数 の替わりに 「頻度分布関数」 や「ヒスとグラム」 を用いてもよいが、 以下には「確 率密度関数」 を用いて説明する。
1 ) 正規分布に従う平均値に変換する方法
状態 kでの x*
ikと p*
ikを次式でそれぞれ正規分布に従う平均値
xik。と μ
Pik。に変換 さ
ここで、 Φ-1は Φの逆関数で、 σ χί0と σ Ρπ)はそれぞれ正規分布に変換された x*ik と p*ikとの標準偏差であり、 次式で求められる。
< ^ 1畴) (3 1) ^ (嘛
( 3 2) 平均値 μ xik。 と μ piko を用いて状態判定や状態予測を行う。
2) 直接変換法
状態 kでの x と p を次式で正規分布の確率変数に変換する。
φ
-1 (Fo (½ ))+ ^ (3 3)
(Pik ))+μ (34) ここで、 S
xkと S
Pkはそれぞれ x*
ikと p*
ikの標準偏差で、 と; u
Pkはそれぞれ x*
ikと p*
ik の平均値である。 x' あるレ、は ρ' を用レ、て状態判定や状態予測を行う。
(2) 基準状態を参照しない場合
まず、 状態 kでの波形データ x*ikと特徴パラメータ p との確率密度関数 (ある いは、 頻度分布) をそれぞれ d)と f (p*ik)とし、 確率分布関数 (あるいは、 累積度数分布) を Fxk(x*ik)と FPk(p*ik)とする。
1) 直接変数変換法
波形データ x*ikと 特徴パラメータ p*ikを次式で正規分布の確率変数に変換する。
二 -
1 (F^ (
Xik ))+μ^ (3 5)
ここで、 Sxkと SPkはそれぞれ x*ikと p*ikとの標準偏差で、 / xkと /zPkは x*ikと p*ikとの平 均値である。 x' ikあるいは ikを用いて状態判定ゃ状態予測を行う。
2) 正規分布に従う平均値に変換する方法
正規分布に従う平均値は次式で求められる。
^pik=Pik ^Pik^'l^pk(p*k ))
(38) ここで、、 とび pikは次式で求められる。 σ^7^ φ(φ"Ι(Ρ ¾ )))
σ = ^ φ(φ-1 )))
Jpk ik) (40) Aixiltと は正規分布に従うので、 /z あるいは AiPikを用いて状態判定や状態予 測を行う。
3) 間接変換法
M個の特徴パラメータ p*ikの最小値と最大値をそれぞれ (p*ik) ^と (p*ik)raaxとす
る。 (p*ik) ^から(p*ik) maxまで N個の等間隔区間に分割する。 各区間の代表値 (例え ば、 中央値) を p*iWとする。 ここで、 j=l〜N。 p*ikjを p*ikの替わりに式 (36) ま たは式 (38) に代入すると、 N個の p',ikまたは ' が得られる。 p"ikまたは μ ' を用いて状態判定や状態予測を行う。
同様に、 Μ個の波形データ x*ikの最小値と最大値をそれぞれ (x*ik)^と (x Laxと する。 (x*ik)minから( ) まで N個の等間隔区間に分割する。 各区間の代表値 (例 えば、 中央値) を x*ikjとする。 ここで、 j=l〜K。 x*ikjを X の替わりに式 (35) または式 (34) に代入すると、 K個の x' ' ikまたは '^が得られる。 χ'Άまた は を用いて状態判定や状態予測を行う。
4) 特徴パラメータ p*ikの区間平均値を求める方法
N個特徴パラメータ p*ikを求めた後、 M組に分割して第 j組の平均値は次のように 求める。
'·=ι (41) ここで、 Ν' は第 j組にある p*ikの数である。 は近似的に正規分布に従うので、 μ を用いて状態判定や状態予測を行う。
波形データ x*ikを式 (29)、 式 (33)、 式 (35)、 式 (37) で正規確率分 布の波形データに変換した ^ 、 X' ik。、 ' ik 、 /i xik , X" ik , ^は 「正規分布 の波形データ」 と呼ぶ。 また、 特徴パラメータ p*iを式 (30)、 式 (34)、 式 (36)、 式 (38) で正規確率変数に変換した μ piko 、 p' 。 、 p' ik 、 μ Plk 、 p' ' は 「正規分布の特徴パラメータ」 と呼ぶ c
3. 正規分布の特徴パラメータによる状態判別方法
ここで、 正規分布の特徴パラメータを用いて、 対象物の状態を判別する方法を 述べる。
(1) 統計理論による判別
1) 正規分布の特徴パラメータの平均値の検定
状態 kと状態 yで求めた正規分布の特徴パラメータをそれぞれ pikと Piyとする。 ここで、 i= l〜M、 Mは使用する正規分布の特徴パラメータの総数を表す。 と piyとの平均値をそれぞれ と iy、 pikと piyとの標準偏差をそれぞれ Sikと Siyとす
る。 一般に J個の Pjの平均値 ίと標準偏差 Sは次の式で計算する。
J
s2=
J-l (43) μ と iyが等しいか否かの検定は次のように行う (参考文献 3参照)
文 3 ) K. A. Brownlee. Statistics丄 Theory and Methodology in Scie nce and Engineering, Second Edition, The University Chicago, 1965
|μ¾- , |>- -/α/2(ΐ-ΐ)
V J (44) が成立すれば、 有意水準ひで 「 /i ikと / Z iyとが等しくない」 と判定する。 ここで、 ta/2(J- 1)は自由度 J-1の t分布の確率密度関数が下側確率 α/2に対するパーセント 点である。
2) 正規分布の特徴パラメータの分散の検定
Sikと Siyが同じか否かの検定は次のように行う (参考文献 4参照)。
(参考文献 4ノ K. A. Brownlee. Statistical Theory and Metnodology in Scie nce and Engineering, Second Edition, The University Chicago, 1965
- ->Fa/2(J-l,J-l) or >Fa/2(J-l,J-l)
^iy1 ¾ (45) が成立すれば、 有意水準 aで 「Sikと Siyとが等しくない」 と判定する。 ここで、 Fa/2(J-1, J-l)は自由度 J- 1の F分布の確率密度関数が下側確率 a /2に対するパー セント点である。
有意水準 aを変えたとき、 式 (44) または式 (45) を満足するか否かを確 認することにより、 状態 yが状態 kに対する状態変化の程度を決める。 状態変化 の程度を有意水準 aにより決定する例は表 1に示す。 なお、 設備診断の場合、 状 態 kを正常状態とすれば、 状態 yは正常状態か、 注意状態か、 危険状態かの判別 については、 表 1のように 「正常」 ( ,)、 「注意」 (ひ 2)、 及び 「危険」 (a 3) の ように設定し検定することができる。 つまり、 式 (44) または式 (45) はひ 'の時に成立しなければ 「正常」 と判定する。 また、 式 (44) または式 (45)
は α; 2の時に成立すれば 「注意」 と、 α 3の時に成立すれば 「危険」 と判定する。 なお、表 1中の αの数値範囲は例であり、設備の重要度などによって決められる。 数個の特徴パラメータを用いて状態判定を行う場合、 判定の結果は最も状態変 化が大きいと示した特徴パラメータの判定結果に従われる。 例えば、 3つの特徴 パラメータ ρ ρ2、 Ρ3 を用いて判定するとき、 Plの判定結果が 「注意」、 ρ2の判定 結果が 「正常」、 ρ3の判定結果が 「危険」 である場合、 最終の判定結果は 「危険」 とする。
表 1
3) 信頼区間による判定
基準時点で測定した波形データから求めた正規分布の特徴パラメータの平均値 を /ii。とし、 他の時点で測定した波形データから求めた正規分布の特徴パラメ一 タの平均値を
ikとすると、
β ί0 の信頼区間は次式で与えられる。
(46) ここで、 /
2(J- 1)は自由度 J- 1の t分布の確率密度関数が下側確率 ct/2に対する パーセント点である。 Si。は波形データから求めた
Pi。の標準偏差である。 は式 (46) に示す区間内にあれば、 1-ひの確率で A 。との差がないという。 /ii。の 99
°/0信頼区間は、 J>10のとき、 近似的に次のようになる。
従って、 は式 (4 7) の範囲を超えれば、 99%の確率で i。と違うと判定す る。 更に、 測定した波形データで求めた の信頼区間は、 次式で求められる。 μ,. ±3S' VJ (48) ここで、 Sikは波形データから求めた pikの標準偏差である。
表 1の αを式 (46) に代入すれば、 次のような信頼区間が得られる。 状態変化がない区間 (正常) :
(49)
状態変化が中である区間 (注意) :
(5 0) 状態変化が大きい区間 (危険) μ!'
0± 17}
σ3' (5 1) kがこれらの区間内にあるか否かによって状態の判定を行う。
数個の特徴パラメータを用いて状態判定を行う場合、 最終の判定結果は式 (4 9), (50), (5 1) に記述される内容と同じである。
(2) 可能性理論による判別
1) 可能性分布関数の作成
状態 kにおける波形データで正規分布の特徴パラメータ P iの値を算出した後、 P iの確率密度関数 f k (p から可能性分布関数 Pk (P i) を式 (5 2) で求 める。 可能性理論によれば、 p iがどのような確率分布に従っても、 その可能性 分布関数が求められる。 p iが正規分布に従う場合、 N段の可能性分布関数 p k (p は次のように求める (参考文献 5参照),
但し、 上式中、 p i
x = min{ p i}+xX (max{pj _min{ p J) N、 x= l〜N, 3
1は の標準偏差, はの pi平均値である。
(参考文献 5) L. Davis: HANDBOOK OF GENETIC ALGORITHMS, Van Nostrand Rei nhold, A Division of Wadsworth, Inc ( 1 990)
2) 可能性の求め方
図 2のように、 状態 kと状態 yで求めた正規分布の特徴パラメータ Piの可能性 分布関数を Pk (pO と Py (ρ とし、 状態 yで求めた正規分布の特徴パラメータ の値を P'iとすると、 「状態 yが状態 kと同じである」 という可能性 wは次のよう
に求められる。
a) p' iの平均値 p' im e a nと Pk (pi) とのマッチングによる wの決定、 b) Py (Pi) と Pk (Pi) とのマッチングによる wの決定。
なお、 Py (Pi) と Pk (Pi) とのマッチングによる wを求める式は次に示す。
+∞ 一 (5 5)
3) 状態変化の判別
状態 kで求めた正規分布の特徴パラメータ Piの可能性分布関数 p k (Pi) が得 られた後、 左右両側の 「状態変化が小さい」 の可能性分布関数 (pcl (Pi) と p c2 (ρΟ)、 および 「状態変化が大きい」 の可能性分布関数 (pdl (Pi) と pd2 (P
0) は図 2に示すように決定する。 境界値の
μ土 ,, μ士ゾ 'S!' (56) における i, jはユーザ入力により決定するが、 標準値として i=3、 j=6とする。 設備診断の場合'、 正常状態の可能性分布関数を P k (p 、 注意状態の可能性 分布関数を p cl (pi) と p c2 (Pi)、 危険状態の可能性分布関数を (pi) と p d2 (pi) とする。 実際の識別時に得られた 「正常」 と 「注意」 と 「危険」 との可 能性は図 2のように表示する。 また、 「危険」 と判定した場合、 警報を出すこと も可能である。
( 3 ) 情報理論による状態判別法
対象物の基準状態における正規分布の特徴パラメータの確率密度関数を f P。 (Pi) し、 基準状態以外における正規分布の特徴パラメータの確率密度関数を (Pi)と する。 なお、 基準状態以外の時点での状態を 「テスト状態」 と呼ぶ。 テスト状態 が基準状態と同じ状態であるか否かは次の ("Kullback- Leibler Information (力 ルバック情報量、 KI)」 と 「Information Divergence (情報ダイパージエンス、 ID) J で判定できる。 -ゾ》 '·
1 , 丄 σ Φθ2+0½ο-μ )2
_ 1 OipO^dlipO-^ipkf^^H^ipO-V-ipk)2
=T{ O ~ipi 一22 "2} ( 5 8 )
Kipと IDP による状態判別法は [1 0] に詳細に記述しているので、 ここで省略 する。 ([ 1 0] 劉 信芳、 豊田利夫、 陳 鵬、 馮 芳、 二保知也: 「Informatio n Divergenceによる回転機械の異常診断」、 精密工学会誌、 Vol.66、 No.1、 2000 年、 .157-162. )
( 4 ) 複数の特徴パラメータを統合することによる状態判別法
数個の特徴パラメータを統合して、 状態変化の有無や状態予測を行うこともで きる。 特徴パラメータの統合法は、 主成分分析法や KL展開法などがあるが、 特徴 パラメータの統合法により求めた新たな特徴パラメータは「統合特徴パラメータ」 とよぶ。 ここで、 主成分分析法の例を示す。 ([ 1 1 ] 大津, 栗田, 関田著:パ ターン認識, 朝倉書店、 1996.)
例えば、設備診断の場合、設備の正常状態下で求めた無次元特徴パラメータ(Pl, p2, '··,ρ∞) に対して、 m個の主成分は次のように表されます。
ζι =άηΡι +α ρ2 + -- + almpm
z2 = a2ipx + a22p2 +■■■ + a2mpm ( 5 o ) zm = am\P\ + amlPl ,"
各主成分 Zl〜Z» は 「統合特徴パラメータ」 とも呼ぶ。
相関行列は次のように求められます。
'11 r12 "' r\n
'21 Γ22 "" ·2η
R = (6 0) ml
rm2 ここで、 正常状態で収録した n組のデータを {p
lk, p
2k, ··· p
mk} k=l,2, ·'·,η とすると、
相関行列 Rの固有べク トル、 ぇ= (λ ,, λ ', λ m) , λ λ λ とすると、
固有値え i に対応する式 (5 9) の係数が求まり、 第 i番目の主成分は次のよ うに求められます。
= αηρλ + ai2p2 +■■■ + aimPm (6 3) 主成分を用いて次式のように状態判定を行う。
一つ 2
≤Z2(K, a)
(6 4) ここで、 Kは使用する主成分の数、 αは有意水準、 %2 (Κ, a) は自由度 Kの力 ィ 2乗分布の確率密度関数が上側確率 aに対するパーセント点である。 αの決定 法は表 1同じである。 なお、 Κ=3に対して α=0· 05のとき、 2 (3,0.05) =7.815 である。
なお、 式 (6 4) に示す判定方法以外に、 特徴パラメータの統合法により求め た統合特徴パラメータ、 例えば、 式 (6 3) に示す主成分、 を正規分布の確率変 数に変換した後、 統計検定や可能性理論により状態判定を行うこともできる。 ここで、 実例を示す。 図 3はある回転機械の正常状態時 (図 3 (a)) と回転軸 ミスァライメント状態時(図 3 (b)) に測定した振動の加速度波形データである。 図 4 (a)及び図 5 (a)は式 (2) 〜式 (1 3) 及び式 (1 8) 〜式 (2 1) に示す 1 4個の無次元特徴パラメータを 6 0回求め、 正規分布変換を行う前に、 K=3の ときに求めた式 (6 4) 右側の値 (6 0個) を示す。 図 4 (b)及び図 5 (b)は式 (3 1) に基づいて正規分布変換を行った後、 K=3のときに求めた式 (64) 右側の値
(9個) を示す。 これらの図により、 正規分布変換後の結果は正規分布変換前の 結果より良好であることがわかる。 なお、 「良好である」 とは、 式 (6 4) 左側 の値は正常状態の時に小さく、 異常状態の時に大きいということである。
4. 正規分布の波形データによる状態判別方法
(1 ) 特徴パラメータによる状態判別法
対象物から計測された波形データ x*i (例えば、 図 6 (a)、 (c)) を式 (3 3) に
より正規分布の波形データ x'ik。 (例えば、 図 6 (b)、 (d)) に変換した後、 正規分 布の波形データ ik。を用いて前述の正規分布の特徴パラメータによる状態判別方 法で状態判定を行うことができる。 なお、 この方法は正規分布の波形データ β xiko と X' ik。に適用できる。
(2) 情報理論による状態判別法
基準状態の正規分布の波形データの確率密度関数を fx。(x とし、 基準状態以外 の時点での正規分布の波形データの確率密度関数を fxk(Xi)とする。 なお、 基準状 態以外の時点での状態を 「テスト状態」 と呼ぶ。 テスト状態が基準状態と同じ状 態であるか否かは次の「Kullback - Leibler Information (力ルバック情報量、 KI) J と 「Information Divergence (情報ダイパージエンス、 ID)」 で判定できる。
„ KTI J c°° t
f^— -dxi
"
∞ fxk{xi) t= ) l ¾C¾rPhi)
2 、o o
一 、 ixi 2 " (66)
KLと IDxによる状態判別法はに詳細に記述しているので、 ここで省略する。 ([1 2] 劉 信芳、 豊田利夫、 陳 鵬、 馮 芳、 二保知也: 「Information Divergen ceによる回転機械の異常診断」、 精密工学会誌、 Vol.66、 No.1、 2000年、 p.157- 162. )
ここで、 実例を示す。 図 6 (a)はある回転機械の正常状態時に測定した振動加 速度の波形データである。 図 6 (c) は同じ回転機械のアンバランス時に測定し た波形データである。
平均値以上の波形データ X 、 平均値以下の波形データ Xi -、 及び絶対値データ |x i| を正規分布の波形データに変換した後、 それぞれ x'i+、 x'i -、 及び il'で表す。
、 X' i -、 及び |xi|,の平均値と標準偏差は次の通りである。
x' ί+:正常時の平均値 =0.69、 アンバランス時の平均値 =1.96
正常時の標準偏差 = 1.03、 寺の標準偏差 =2.37
:正常時の平均値 = -0.37、 寺の平均値 = -0.66
正常時の標準偏差 =0.47、 I寺の標準偏差 =0.69
|Xiに :正常時の平均値 =1.45、 アンバランス時の平均値 = 3·29
正常時の標準偏差 =1.52、 アンバランス時の標準偏差 =3.53
(参考文献 8) によると、 過検出率ひ =0.15、 見逃率 =0.15とすると、 ΚΙ〉1. 21或いは ID〉2.43ならば、 85%ηの確率で 「テス ト状態は基準状態と違う」 と判定 できる。 上記の x'i+、 x,i -、 及び |xi|'の KIと IDは次の通りである。
x' i+: KI=0.519, ID=2.24
x, i-: KI=0.53, ID=0.40
|xi|' : KI=0.57, ID=2.67
|xi|,の ID〉2.43から、 85%nの確率で 「図 6 (a)の状態は図 6 (c)の状態と違う」 と判定できる。
また、 式 (44) と式 (45) に示す平均値と分散の検定により状態判定を行 つてもよレヽ。
更に、 上記の x'i+、 x'i -、 及び |xi|,を用いて、 式 (2) 〜式 (25) で特徴パ ラメータを求めて、 統計検定や可能性理論及び特徴パラメータの統合法により状 態判定を行うこともできる。
5. 状態予測
正規分布の特徴パラメータ Pi が求められた後、 従来の状態予測手法を用いて、 測定対象の状態を予測することが行える ([ 1 3] 石川真澄、 武藤博道:予測手 法、 計測と制御、 1982.3. [14] Ogawa, M. : Time series analysis and sto chastic prediction, Bull. Math. Stat. , 8, 8—72, 1958. [1 5] B.ピカニォ ル著、 小村賢二、 柴山宫恵子訳:予測のための統計学、 晃洋書房、 1987.) 図 7は状態予測の方法を示す。 各測定時点(Xl〜x で測定した特徴パラメータ 値、 或いは、 主成分値を正規確率分布に従うように変換した後、 回帰解析により 予測曲線とその信頼区間を求め、 「寿命限界」 との交点で 「最短寿命」、 「平均寿 命 J 及び 「最長寿命」 を求める。
図 8はある回転機械が正常状態からアンバランス状態へ変化していく間に 8回 測定した波形例である。 測定 1は正常状態であり、 測定 8は最も程度の重いアン バランス状態である。 図 9は式 (2) 〜式 (1 3) で求めた、 各測定における特 徴パラメータの平均値を示す。 この例でわかるように、 特徴パラメータは異常状 態の程度が重くなつていくにつれて値が単調増加の特徴パラメータ (例えば、 p4,
p5) があれば、 値が単調減少の特徴パラメータ (例えば、 ρ2, Ρ8, Ρ9, Ρπ») もある。 また、 異常状態の程度が重くなつても値がほぼ不変な特徴パラメータ (例えば、 p., ρτ) もある。 従って、 値が単調増加 (或いは、 単調減少) の特徴パラメータを 選んで、 寿命予測を行うべきである。 また、 主成分により寿命予測を行う時、 値 が単調増加の特徴パラメータのみ (或いは、 単調減少の特徴パラメータのみ) を 選んで、 主成分を求めて寿命予測を行うべきである。
上記に述べた波形データと特徴パラメータを正規分布の確率変数に変換する手 法、 状態判定手法及び状態予測手法を実現するための計測と処理の流れは図 1 6 (a)に示す。 また、 図 1 6 (a)を実現するための波形データ計測及び状態判定装置 の回路は図 1 7に示す。
6 . 脈動信号から特異成分の検出法及び状態判定法
ここで、 ガスエンジン失火診断の例を用いて説明する。
図 1 0 (a)、 図 1 1 (a)、 図 1 2 (a)、 図 1 3 (a)はガスエンジン失火 (図中のピ ーク値の異常個所) が発生したときのシリンダ圧力波形データである。 なお、 図 12 (a) , 図 13 (a)は多数の異常ピークが連続して発生する例を示す。 このような脈 動信号に発生する局所異常 (特異成分) に対して、 オンライン特異点の検出及び 状態判定の流れは図 1 5に示す。 以下、 その手順を説明する。
1 ) . 診断ための準備 (図 1 5 (a) )
( 1 ) 診断対象の脈動信号と回転パルス信号を同時に計測する。 なお、 回転パ ルス信号は周期パルス信号とも言い、 脈動信号の各ピーク値のタイミングを決め るときに用いられる。
( 2 ) 雑音を除去するために, ローパスフィルタリングを行う。 ローパスフィ ルタの力ットオフ周波数 fLは次のように決定する。
ここで、 nは軸の回転数 (rPm)、 zは 1回転毎のピーク数 (ピーク / 1回転)、 fOは余裕周波数 (> n/60、 フィルタリング後のノイズ除去効果を観察することに より決定する) である。
( 3 ) ローパスフィルタリング後の包絡線波形データ、 または、 ピーク波形デ ータ、 または、 移動平均波形データを求める。
(4) 上記の包絡線波形データ、 または、 ピーク波形データ、 または、 移動平 均波形データの正規化を次式のように行う。
' (り- μ' (り
χ(り
S (68) ここで、 x(t)は正規化後の波形データ、χ' (t)は元の波形データ、 μ, (t)は χ' (t) の平均値、 sは x' (t)の標準偏差である。 なお、 ;U ' (t) と sを求めるとき、 正常状 態の波形データが望ましいが、 正常状態の波形データでなくてもよい。 また、 μ ' (t)と sを求めるとき、 運転条件 (回転数や負荷) をできるだけ判定のときと同 じょうに設定する。 但し、 運転条件が多少変化しても判定結果に大きく影響しな い。
(5) 回転パルス信号のピークと脈動信号のピークとの相対位置関係を調べ て置き、異常を識別するための閾値を決定する。閾値の決定は次のように考える。
I x(t) I >ka (69) のとき、 異常が発生したと判定する。 ここで、 k (=2〜4) は診断対象物によって 決定する。
I x(t) I 力 より大きい間、 回転パルス信号に対応する脈動信号のピークの 所に異常があると判定する。
2). オンライン診断 (図 1 5 (b))
( 1 ) 判定対象物の脈動信号と回転パルス信号を同時に計測する。
(2) ローパスフィ^/タリングを行う。 口一パスフィ^^タのカッ トオフ周波数 fLは式 (66) のように決定する。
(3) ローパスフィルタリング後の包絡線波形データ、 または、 ピーク波形デ ータを求める。
(4) 上記の包絡線波形データ、 または、 ピーク波形データの正規化は式 (6 8) により求める。
(5) 正規化された包絡線波形データ、 または、 ピーク波形データの絶対値は 閾値 (ka) より大きいか否かを監視する。 大きくなければ、 異常がないと判定 する。 大きければ、 異常が発生したと判定し、 | x(t) I 力 ¾びより大きい間、 前 もって測って置いた回転パルス信号と脈動信号のピークとの関係により、 異常の 箇所を判定する。
上記の手順に基づく信号処理の例を図 1 0、 図 1 1、 図 1 2、 図 1 3に示す。 図 1 0 (a) (b;)、 図 1 1 (a) (b)、 図 1 2 (a) (b)、 図 1 3 (a) (b)ガスエンジン失火
(図中の異常箇所) が発生したときのシリンダ圧力の時系列波形データとスぺク トルである。
図 1 0 (c) (d)、 図 1 1 (c) (d)、 図 1 2 (c) (d)、 図 1 3 (c) (d)はローパスフィル タにより雑音を除去した時系列波形データとスぺク トルである。
図 1 0 (e)、図 1 1 (e)、図 1 2 (e)、図 1 3 (e)は回転パルス波形データである。 図 1 0 (f)、 図 1 2 (f)は包絡線波形データである。
図 l l (f)、 図 1 3 (f)はピーク波形データである。
図 1 0, 1 1, 1 2 , 1 3の(f)で分かるように、 異常 (ピークの欠落、 また は、 小さくなる) が発生した場合、 その絶対値 1 x(t) Iが 2 σより大きくなるこ とが判明できる。
また、 絶対値 I x(t) Iが 2 σより大きくなる箇所と回転パルス波形データとの 対応関係により異常箇所を判明できる。
上記に包絡線波形データとピーク波形データの例を示したが、 図 1 4に示すよ うに、 移動平均波形データ ( ) を求めたら、 x(t)の替わりに u(t)を用いて上記 の手順と同様に脈動信号の特異点検出ができる。 なお、 カットオフ周波数 (f。) は式 (6 7) により求めたら、 移動平均の点数 (M) が求められる。
上記に述べた脈動信号の特異成分の検出と判定方法を実現するための計測と処 理の流れは図 1 6 (b)に示す。 また図 1 6 (b)を実現するための脈動信号計測及び 状態判定装置の回路は図 1 8に示す。
7. 状態判定装置、 または、 オンライン状態判定システムの処理流れ
状態判定装置装置、 または、 オンライン状態判定システムの処理流れは図 1 6 に示す。 図 1 6 (a)に示すように、 測定波形データからノイズを除去し、 特徴パ ラメータを求め、 その特徴パラメータを正規分布の確率変数に変換した後、 統計 検定や可能性理論や特徴パラメータの統合により状態変化の有無判定を行う。 こ れらの処理は計算機または専用装置で実現できる。 また、 図 1 6 (b)に示すよう に、脈動信号からの特異成分検出及び状態判定は、 ローパスフィルタ、包絡線(ま たは、 ピーク値) はハードウェアで実現できる。 数値演算装置(または、 計算機) は正規化した x(t)、 I x(t) > o Iの判定、 特異箇所の検出と判定および判定結
果の表示を行えばよいから、 リアルタイムに行える。 七. 図面の簡単な説明
図 1は、 測定した波形データを 4種類、 つまり、 平均値より大きいデータ Xi+、 平 均値より小さいデータ X; -、 式 (1 ) による正規化後の絶対値データ 及び全体 波形データ X*Mに分ける例を示すグラフである。
図 2は可能性分布関数の例を示すグラフである。
図 3は正常状態と回転軸ミスァライメント状態の振動波形データの例を示すグラ フである。
図 4は正常状態の振動波形データに対する主成分分析結果を示すグラフである。 図 5は回転軸ミスァライメント状態の振動波形データに対する主成分分析結果を 示すグラフである。
図 6は正常状態とァンバランス状態の振動波形データの例を示すグラフである。 図 7は状態予測方法を示すグラフである。
図 8は正常状態からァンバランス状態へ変化していく時の振動波形データの例を 示すグラフである。
図 9は各測定の波形の特徴パラメータ値を示すグラフである。
図 1 0は包絡線波形データによる 1個のピークが欠落した場合の信号処理例を示 すグラフである。
図 1 1はピーク波形データによる 1個のピークが欠落した場合の信号処理例を示 すグラフである。
図 1 2は包絡線波形データによる多数のピーク異常の場合の信号処理例を示すグ ラフである。
図 1 3はピーク波形データによる多数のピーク異常の場合の信号処理例を示すグ ラフである。
図 1 4は移動平均波形データの求め方を示すグラフである。
図 1 5は脈動信号の特異成分検出及び異常診断処理の流れを示すフローチヤ一ト である。
図 1 6は状態判定装置、 または、 オンライン状態判定システムの処理流れを示す フローチヤ一トである。
図 1 7は信号計測及び状態判定装置の回路の一例を示す回路図であり、 図中の符 号は次の通り ある。 .
1 センサ、 2 アンプ、 3 フィルタ、 4 信号処理 '計算装置、 5 表示出 力装置、 6 データ用 RAM、 7 AD変換器、 8 DCポート、 9 S C I (シ リアル . コミュニケーション . インタフェース)、 1 0 CPU、 1 1 フラッ シュ ROM、 1 2 外部コンピュータ.
図 1 8は脈動信号計測及び状態判定装置の回路の一例を示す回路図であり、 図中 の符号は次の通りである。
1 センサ、 2 アンプ、 3 フィルタ、 4 包絡線処理部、 5 周期パルスセ ンサ、 6 表示出力装置、 7 AD変換器、 8 DCポート、 9 SC I (シリ アル · コミュニケーション ' インタフェース)、 1 0 CPU, 1 1 フラッシ ュ ROM、 1 2 外部コンピュータ.、 1 3 データ用 RAM、 14 信号処理 ' 計算装置.