音響診断のメカニズムと特徴,今後の展望 | コンディションモニタリングBOX | ジュンツウネット21

機械が発生する騒音(音響情報)は,容易にしかも機械に非接触で測定することが出来,診断システムを構築する際に有用である。従って設備の広域のモニタリングに適している。
ここでは,音響診断としてパターン認識及び最適化手法による音源同定法についてそれらの特長と解析結果の例,及び展望などについて述べる。

広島大学 大学院 中川 紀壽  2007/2

はじめに

近年の生産設備の自動化に伴い,機械の異常診断の自動化が進められているが,異常の判別には加速度等の振動情報が用いられることが多い。しかし,振動情報を用いて機械のモニタリングを行うには,異常信号を比較的敏感に検出出来るという利点があるが,大型設備の診断に際しては,センサを各機械に設置する必要があり,多大な投資が必要となる。

それに対して機械が発生する騒音(音響情報)は,容易にしかも機械に非接触で測定することが出来,診断システムを構築する際に有用である。従って設備の広域のモニタリングに適している。著者は*1,*2転がり軸受が発生する音の特徴周波数の時間的変化を用いた診断法に関する研究を行ってきている。

しかしながら,音響情報の測定は振動情報と比較して,S/N比が低く,正確な診断を行うためには,信号処理等の技術が必要となる。その一例として,Loukisら*3はパターン認識の手法を用いてガスタービンの異常診断を行っている。ここでは,音響診断としてパターン認識及び最適化手法による音源同定法についてそれらの特長と解析結果の例,及び展望などについて述べる。

1. パターン認識による音響診断*4

1.1 対象例

ここで示す事例の対象は図1に示すように,スピンドル,交流モータ及び軸受荷重負荷機構からなる。音響情報を測定するための供試軸受はスピンドル軸に取り付けられ,軸回転数は1,800rpmである。また,供試軸受にスラスト及びラジアル方向に軸受荷重を負荷するための機構があり,負荷荷重はロードセルを用いて測定する。また,実験は無響室内で行われている。

Experimental apparatus
図1 Experimental apparatus

供試軸受は呼び番号6010の単列深溝玉軸受で,傷のない正常な軸受及び外輪軌道面,内輪軌道面または転動体表面にそれぞれ1ヵ所の傷を有する軸受である。

1.2 音響診断

音響情報のパターン認識を用いて軸受の異常診断を行う際,パワースペクトルにおいて軸回転周波数及び軸受特有の周波数の整数倍成分に着目し,以下に示す縮小パターン及び参照パターンを用いて診断を行う。

i)縮小パターンの作成
 まず,軸受騒音の周波数分析を行った結果に次のような処理を施し,着目周波数に対して異常によるスペクトルの変化分を求める。

 (1)

 (2)

ただし,j=1,2,3,…, Nである。また,sp(f)及びsph(f)はそれぞれ異常軸受及び正常軸受のパワースペクトルである。fkは,縮小パターンの作成時に着目する周波数で,次節 ii)で述べる。ここで簡単のため,p( j )=p( j×fk )とおくと,

P=(p(1), p(2),…, p(N))  (3)

が得られる。式(2)のh(f)はパワースペクトルから周波数fkの整数倍の成分のみを抽出するフィルタ関数である。また,式(3)のベクトルPを縮小パターン(Reduced Pattern)と呼ぶ。

ii)着目周波数fkの選択
 縮小パターン作成時に着目する周波数fkには,fr,fo,fi 及びfbを考える。ただし,frは軸回転周波数であり,転動体の公転周波数fo,保持器に対する内輪の相対回転周波数fi,及び転動体の自転周波数fbはそれぞれ以下の式で表せる。

 (4)

ここで,Dはピッチ円直径,dは転動体直径,aは接触角である。外輪,内輪あるいは転動体に1ヵ所の傷がある場合に発生する振動及び音の周波数は,zを転動体の個数とすると,zfo,zfi 及び2fbであるが,ここでは,式(4)の値とその整数倍の成分に着目する。

iii)パターンによる認識の方法
 式(3)で得られる縮小パターンを各傷(内輪傷,外輪傷,転動体傷)及び様々な実験条件に対しあらかじめ求めておき,各傷ごとに式(5)を用いてまとめ,実際の診断に用いるパターンを作成しておく。このパターンを参照パターン(Reference Pattern)と呼び,RPで表す。実際に診断を行う際には,診断する軸受騒音の縮小パターンと各傷ごとの参照パターンを比較する。

ここで,RPi は傷i (3種類)の参照パターン,Pih は傷i を有する場合の実験条件hから得られる縮小パターン,及びNi は参照パターンを作成するために行った各傷ごとの実験条件数である。

パターン認識を用いた診断の際には,次に示す2つの判定値を用いる。

(1)ユークリッド幾何学的距離dE

 (6)

式(6)中のRPij は,傷i の参照パターンRPi の,及びPj (=p( j ))は診断を行う縮小パターンのそれぞれj 番目の要素である(j=1,2,3, … , N,Nは使用する周波数成分の個数,式(2)参照)。

(2)相関係数dC

 (7)

ここで,

である。

これら2つの判定値を用いて診断を行い,式(8)を満たすようなdERPiP)及びdCRPi P)を持つ傷i を診断判定結果とする。

 (8-1)

 (8-2)

ここでfは傷の種類を示す(f=1,2,3)。すなわち,式(8-1)ではユークリッド幾何学的距離の小さい事象が,また式(8-2)では相関係数の大きい事象が診断結果として採択されることになる。

iv)統計的パターンによる認識の方法
 この方法ではユークリッド幾何学的距離と相関係数の2つの確率変数を要素に持つベクトルX の同時確率密度関数を用いて診断を行う。

X={dE, dC}  (9)

ここでは,Xが多次元正規分布に従うものと考え,同時確率密度関数を次式で定義する。

 (10)

ここで,| | は行列式を,( )-1は逆行列を,及び( )T は転置行列を示す。n は式(9)の要素数である。また,

である。Xih は傷i を有する場合の条件hのデータから学習段階であらかじめ計算されるベクトルである。Ni は各傷ごとの実験条件数である。共分散行列Ci もまた学習段階でそれぞれの異常i に対し,次式から算出される。

 (11)

診断は,未知のXと学習段階で求めた,Ci を用いて式(10)を計算し,pr( X | i )が最大となるものを,原因i による異常として判定することで行う。

1.3 音響診断結果

まず,各傷の参照パターンと,それら参照パターンを計算する際に用いる縮小パターンの各周波数における標準偏差を図2に示す(fk =fr の場合)。図2(a)に示すように,参照パターンが傷ごとに異なることを利用して診断を行う。また,図2(b)に示す標準偏差は,周波数が1,000Hzを超えるあたりより値が大きくなる部分の存在することがわかる。ここでは,音響パワースペクトルの1,110Hzまでの値を用いて,縮小パターン及び参照パターンを作成し,診断に用いる。なお,音響パワースペクトルの更に高周波数の成分を用いて,実際に診断を試みたが,結果として診断精度が低くなった。

Reference pattern and its deviation/(a)Reference pattern
(a)Reference pattern
Reference pattern and its deviation/(b)Standard deviation
(b)Standard deviation
図2 Reference pattern and its deviation

また,診断時に用いた着目周波数を表1に,また参照パターンを求める際に行った実験の条件(Test count)を表2に示す。

表1 Selected frequency in diagnosis
 
Outer race fault
Inner race fault
Ball fault
Case-1
fr
fr
fr
Case-2
fo
fi
fb
表2 Test count and test condition
Bearing
Thrust N
40
80
120
160
200
Normal,O-1,I-1,B-1
1
2
3
4
5
O-2,I-2,B-2
6
7
8
9
10
O-3,I-3,,B-3
11
12
13
14
15

1.4 パターン認識による音響診断結果

 ユークリッド幾何学的距離と相関係数を用いて外輪傷軸受の診断を行った結果をそれぞれ図3(a)及び(b)に示す。着目周波数は軸回転周波数frである(Case-1,表1参照)。図の横軸のTest countは実験条件(傷の大きさと軸受に負荷するスラスト荷重)を変化させて行った実験番号である(表2参照)。

Outer race fault bearing (Case-1)/(a)Euclidean distance
(a)Euclidean distance
Outer race fault bearing (Case-1)/(b)Correlation cofficient
(b)Correlation cofficient
図3 Outer race fault bearing (Case-1)

図3(a)での診断結果は,0に一番近い値をとる場合に,その傷を有することになる。また,同図(b)では1に一番近い値をとる場合に,その傷を有することになる。

図3の記号のうち,中実の記号(●:外輪傷)が診断結果として正しいものであり,中空の記号(□:内輪傷,△:転動体傷)が診断結果として表れると誤診断となる。図(a)では15個のTest countのうち2ヵ所,図(b)では1ヵ所が誤診断であった。

図4に転動体傷を有する軸受の診断結果を示す。この場合にも着目周波数はfrである。記号のうち,▲が正解であり,○と□が診断結果として選択された場合には,誤診断となる。この場合には,ユークリッド幾何学的距離を使用した場合には,すべての場合に正解が得られた。しかし相関係数を使用した場合には,Test count 11が外輪傷と判断されており,誤診断となる。

Ball fault bearing (Case-1)/(a)Euclidean distance
(a)Euclidean distance
Ball fault bearing (Case-1)/(b)Correlation cofficient
(b)Correlation cofficient
図4 Ball fault bearing (Case-1)

なお,軸受に特徴的な周波数成分を用いた診断も行った(表1のCase-2の場合)が,軸回転周波数成分を用いた場合(表1のCase-1,図3及び4)と比較すると正解率は低くなった。

1.5 統計的パターン認識による診断結果

統計的パターン認識を用いた診断結果を図5及び図6に示す。図5は,着目周波数に軸回転数frを用いた場合である(Case-1)。また,図6は着目周波数に軸受特有の周波数であるfo,fi 及びfbを用いた場合である(Case-2)。両図とも(a)は外輪傷,(b)は内輪傷,そして(c)は転動体傷の診断結果を示す。pr( X | i )の値が一番大きな傷i が診断結果となる。両図とも良好な診断結果が得られていることがわかり,統計的手法を用いることが有効であることがわかる。

Statistical pattern (Case-1)/(a)Outer race fault bearing
(a)Outer race fault bearing
Statistical pattern (Case-1)/(b)Inner race fault bearing
(b)Inner race fault bearing
Statistical pattern (Case-1)/(c)Ball fault bearing
(c)Ball fault bearing
図5 Statistical pattern (Case-1)
Statistical pattern (Case-2)/(a)Outer race fault bearing
(a)Outer race fault bearing
Statistical pattern (Case-2)/(b)Inner race fault bearing
(b)Inner race fault bearing
Statistical pattern (Case-2)/(c)Ball fault bearing
(c)Ball fault bearing
図6 Statistical pattern (Case-2)

2. 観測点選択手法による音響診断(音源同定)

機械システムにおいて,通常とは異なる音を発する機械を発見することは,異常診断の第一歩であり,そのため,音源同定法は,機械の音響異常診断に必要なものである。そこで,ここでは,複数個の音源同定問題に関して,音源個数,音源位置ならびに振動速度を併せて同定する方法について述べる。

2.1 最適化手法(従来法及び観測点選択手法)

 ここでは音源同定問題を,評価関数を用いた最適化問題として扱う。評価関数は,数値シミュレーションにて仮想音源のパラメータを用いて境界要素法で求める,観測点位置での音圧と,観測した音圧との残差を用いて定義する。

従来の研究で用いられている評価関数Wは,

 (12)

である。ここで,Mは観測点数,pci とpmi はそれぞれ,観測点i での仮想音源から計算した音圧と,観測した音圧である。しかし式(12)で定義した評価関数を用いて,音源同定を行うと,最適問題を解く際に局所解問題により精度の悪いことが多い。

そこで,本研究では以下の評価関数と,同定計算用の観測点選択手法アルゴリズムを新たに用いる。観測点選択のための新しい評価関数としては,各観測点i の計算音圧と観測音圧の残差自乗和wi とする。同定計算の手順を次に示す。

(1) 式(12)の評価関数を用いて,計算が収束するまで同定を行う。

(2) 式(13)の評価関数と判断基準値を比較し,その重要度に応じて観測点を選択し,選択された観測点の音圧情報のみより再度評価関数を式(12)を用いて計算し直し,同定計算が収束するまで同定を行う。

(3) 再度すべての観測点情報を用いた計算手順(1)を行い,(2)と(1)の計算結果が収束するまで,(2)と(1)の計算を繰り返す。

手順(2)で行う,観測点選択時に用いる観測点iに対する評価関数を次式に示す。

wi = ( pfi - pmi )2  (13)

ここでpfi は,計算法としてはpci と同じであるが,計算手順(1)が収束した際の最終の値と定義する。観測点の選択に際しては,次式で定義する幾何平均値Gmを判断基準として用い,これよりもまだ残差wi の大きい観測点を重要度が高いと判断し,選択することにする。

 (14)

この幾何平均を用いることにより,各wi にオーダ的な差があり,算術平均では最大値を持つ観測点しか選択されない場合を,避けることが出来る。

また,上記の評価関数を最小化する際には共役勾配法を用いる。

2.2 数値実験による同定結果

以下に数値実験によって得られた同定結果を示す。

同定を行う音場と18点の観測点位置は,図7表3に示す通りである。各音源の大きさは,高さ0.1m,幅と奥行が0.2mの直方体であり,各音源は,上面のみが周波数100kHzで,振動速度振幅が0.1m/sでピストン振動し,音を放射している場合を考える。

Sound field model
図7 Sound field model
表3 Locations of observation points

No.
x,y,z m
1
0.3,0.3,0.5
2
0.3,0.8,0.5
3
0.3,1.3,0.5
4
0.8,0.3,0.5
5
0.8,0.8,0.5
6
0.8,1.3,0.5
7
1.3,0.3,0.5
8
1.3,0.8,0.5
9
1.3,1.3,0.5
10
0.3,0.3,1.0
11
0.3,0.8,1.0
12
0.3,1.3,1.0
13
0.8,0.3,1.0
14
0.8,0.8,1.0
15
0.8,1.3,1.0
16
1.3,0.3,1.0
17
1.3,0.8,1.0
18
1.3,1.3,1.0

同定結果としての音源位置を図8に示す。同図(a)は従来からの評価関数を用いた方法で,(b)は観測点選択手法の結果である。斜線を引いた四角が音源位置の同定結果で,格子を引いた四角が真の音源位置である。音源位置の初期値は,(x0 , y0)=(0.8 , 0.8)(m)である。

Identification of sound source/(a)Conventional method
(a)Conventional method
Identification of sound source/(b)Microphone selection method
(b)Microphone selection method
図8 Identification of sound source

同定過程の位置を見ると(図8),従来法では正確に同定できなかった結果の後,観測点選択手法を導入することで最終的に正確な位置に同定出来たことがわかる。

音源個数,音源の振動速度についても精度よく求めることが出来る(表4)。

表4 Identification of 4 sound sources
Sound source
Exact
Identified
x,y m
v m/s
x,y m
v m/s
A
1.300,1.200
0.095
1.300,1.202
0.095
B
0.300,0.400
0.090
0.299,0.400
0.090
C
0.500,1.300
0.100
0.500,1.300
0.100
D
1.200,0.500
0.105
1.200,0.501
0.106

3. 観測点重み付け法による音響診断(音源同定)

2.においては,評価関数で使用する測定値としてシミュレーション値を用いたが,ここでは実際に測定した値を用いることとし,更に,複数個の音源がそれぞれ異なった周波数成分を有する音を同時に発生している場合について,各音源の位置や音響スペクトルを同時に同定する。これは複数個の機械の同時リモート診断をも念頭においたものとなる。また,最適化のための評価関数としては観測点重み付け法を用い,観測点の重要度を重みとして評価関数の再構築を行う。

3.1 観測点重み付け法と評価関数

観測点重み付け法では,次に示す2通りの評価関数を交互に使用する。

 (15)

 (16)

ここで,Wは従来法の評価関数である。また,Mは観測点数,pci とpmi はそれぞれ仮想音源から計算によって求めた観測点i での音圧と,観測した音圧である。また式(16)中のwi は,

wi = (pfi - pmi2  (17)

で表す。式(17)中の記号は後述する。

観測点重み付け法の計算アルゴリズムは,以下の手順(1)~(3)から成っている。

(1)まず,式(15)で定義する従来法の評価関数Wを用いて,計算が収束するまで同定を行う。この段階での結果は,従来法を用いた同定結果と等しい。

(2)式(17)を用いて,その重要度に応じて観測点に重み付けを行った評価関数W2により,計算が収束するまで同定を行う。

(3)再度すべての観測点情報を等しい重み付けで用いた計算手順(1)を行い,(2)と(1)の計算結果が収束するまで,(2)と(1)の計算を繰り返す。

3.2 同定結果

実験装置を図9に示す。スピーカには,あらかじめコンピュータを用いて作成し,DATに録音した音を再生し,アンプにより増幅して使用する。なお,録音した音は,純音ではなく,音源ごとに異なる周波数成分を持つ音である。

Experimental apparatus
図9 Experimental apparatus

2.で述べた場合とほぼ同じ実験条件で,4つの音源が存在する場合において,音源の位置と振動速度の同時同定を行う。周波数成分は音源によって異なるが,100Hzから1,000Hzまでの100Hzきざみでいくつかのピークを有する音が,4音源同時に存在する場合である。同定結果を表5に示す。すべての音源が精度良く求められている様子がわかる。しかし,振動速度振幅の同定結果は,比較的精度が低い。上記の同定結果について,次にすべての周波数成分に関して,振動速度振幅を同定した結果を図10及び図11に示す。また図11には図10で同定した値を用いて計算した,音源より0.4m離れた場所での各音源の音圧を示す。特に,音圧での表示の際には,各周波数とも誤差は2,3dBであり,実用上問題の無い範囲の誤差である。

表5 Identification result
Sound source
Exact
Weight method
x,y m
v m/s
x,y m
v m/s
A
0.900,1.300
0.0275
0.928,1.317
0.0217
B
0.300,0.400
0.0237
0.243,0.388
0.0207
C
0.300,1.300
0.0255
0.343,1.405
0.0259
D
1.300,0.300
0.0243
1.282,0.215
0.0221
Velocity identification/(a)Sound source A
(a)Sound source A
Velocity identification/(b)Sound source C
(b)Sound source C
図10 Velocity identification
SPL identification/(a)Sound source A
(a)Sound source A
SPL identification/(b)Sound source C
(b)Sound source C
図11 SPL identification

従って,ここで示した観測点重み付け法は,複数個の音源が存在する音場において,音源位置の同時同定だけではなく,振動速度振幅及び音圧のスペクトル同定に際しても有用であることがわかった。この結果は,機械のリモート診断の際に有用な技術になると思われる。

4. 複数個音源のパワースペクトルの同時同定結果

4音源が同時に存在する場合で,100Hzから1000Hzの間に,10Hzきざみでピークが存在する複雑な場合の各音源のパワースペクトル同時同定結果を図12に示す。真値と同定値の間には数dBの差はあるが,周波数成分の有無や分布形状の同定は精度の高いことがわかる。

Sound spectrum identification/(a)Sound source A
(a)Sound source A
Sound spectrum identification/(b)Sound source B
(b)Sound source B
Sound spectrum identification/(c)Sound source C
(c)Sound source C
Sound spectrum identification/(d)Sound source D
(d)Sound source D
図12 Sound spectrum identification

おわりに

転がり軸受が発生する音を利用して異常診断を行う際に,音の周波数分析結果の着目周波数とその整数倍の成分を用い,縮小パターン及び参照パターンを作成し,診断軸受の縮小パターンと比較することで診断を行い,パターン認識を用いた診断法が有効であることを示した。また,統計的パターン認識の手法を用いると,着目周波数に無関係に精度良く診断が行えることがわかった。

音源同定においては,数値実験に対して,提案した観測点選択手法を用いて,音源個数,位置および振動速度の同時同定を行い,精度よく同定できることを示した。

評価関数と計算アルゴリズムに観測点重み付け法を使用すると,音源位置の同時同定に加え,各音源の振動速度振幅の周波数分布,あるいは音圧の周波数分布を求めることが出来る。

なお,この分野の研究にウェーブレット変換やカオスを取り入れた新しい研究がある。ただし,転がり軸受に生じる振動データを用いたものであるが,ウェーブレット変換を基にした解析により,傷の有無の判別が出来る*7。また,カオス解析において,適切な時間遅れを選択することで,軸受の運転状況をアトラクタとして可視化が出来,更に最大リアプノフ指数が異常診断に有効*8となる,などである。これらの研究が進めば,音響データを用いたリアルタイムでの診断も考えられることになると思われる。

本稿の執筆にあたりご協力を頂いた,広島大学大学院工学研究科の関口泰久助教授に感謝致します。

 

<参考文献>
*1 Nakagawa N. and Mahrenholtz O. , Diagnosis of Machines by Using Sound Information, Proc. 8th World Congress on the Theory of Machines and Mechanisms, Vol.6,( 1991),99-102.
*2 Nakagawa, N. , and Sekiguchi, Y. : Diagnosis of Ball Bearing by use of Sound Information, Proceedings of 11th International Acoustic Emission Symposium, (1992), 361-367
*3 Loukis, E. , Mathioudakis, K。, and Papailiou, K. D.: A Procedure for Automated Gas Turbine Blade Fault Identification Based on Spectral Pattern Analysis, Trans. ASME, J. Eng. Gas Turbines Power,114, 2( 1992), 201.
*4 関口泰久,中川紀壽,パターン認識による転がり軸受振動の音響異常診断,設計工学,(2004),第39巻,第4号,182-188.
*5 関口泰久,中川紀壽,増田健志,最適化手法を用いた複数個音源の同定(第1報 数値シミュレーション,観測点の重要度を考慮した同定法),設計工学,(2004),第39巻,第5号,274-280.
*6 関口泰久,小西琢哉,里信純,中川紀壽,最適化手法を用いた複数個音源の同定(第3報 音源のパワースペクトルの同時同定),設計工学,(2006),第41巻,第1号,42-47.
  関口泰久,中川紀壽,吉田博一,猿渡直行,転がり軸受振動のカオス解析と異常診断,設計工学,(2004),第39巻,第2号,100-107.
*7 関口泰久,板倉有吾,里信純,中川紀壽,実信号マザーウェーブレットによる転がり軸受振動の相関解析と異常診断, 日本設備管理学会誌,(2006),Vol.18,No.2,60-65.
*8 関口泰久,中川紀壽,吉田博一,猿渡直行,転がり軸受振動のカオス解析と異常診断,設計工学,(2004),第39巻,第2号,100-107.