JPH11118661A

JPH11118661A - 振動特性解析装置

Info

Publication number: JPH11118661A
Application number: JP9287161A
Authority: JP
Inventors: Mitsuo Iwahara; 光男岩原; Akio Nagamatsu; 昭男長松
Original assignee: Isuzu Motors Ltd
Current assignee: Isuzu Motors Ltd
Priority date: 1997-10-20
Filing date: 1997-10-20
Publication date: 1999-04-30

Abstract

(57)【要約】【課題】周波数応答関数データと該周波数応答関数の理
論値との重み付き残差二乗和を最小とする最小二乗法に
よりモード特性パラメータを計算することにより該被試
験物のモード特性を同定する演算装置が、該周波数応答
関数の選択された固有モードに対応する虚数部の極性に
＋１又は−１をそれぞれ乗じて正極性にした周波数応答
関数を重ね合わせることにより該固有モードに対応する
周波数応答関数をそれぞれ求め、該周波数応答関数の数
を該固有モードの数に縮小する振動特性解析装置におい
て、周波数応答関数の縮小に対応できる重み関数を理論
的に導出する。【解決手段】多数の該周波数応答関数を同時に考慮した
重み関数を求めて該残差二乗和に適用する。

Description

【発明の詳細な説明】

【発明の属する技術分野】本発明は振動特性解析装置に
関し、特にエンジンブロック等の被試験物を加振するこ
とにより応答を得てその応答関数から被試験物のモード
特性を同定する振動特性解析装置に関するものである。

【０００１】実際の機械構造物の動特性を解析するため
には、該構造物を加振試験することにより得られた振動
応答特性から該構造物の動特性を同定する必要がある
が、この場合、振動現象を直接表現するパラメータであ
り、現象を直接理解し易く、収束性が良好なモード特性
（固有振動数、モード減衰比、固有モード形状などのモ
ーダルパラメータ）を同定する手法が現在では主流とな
っている。

【０００２】

【従来の技術】図１は従来より知られている振動特性解
析装置の構成を概略的に示したもので、まず被試験物１
とこの被試験物１を加振するためのインパルスハンマー
２とを用意する。

【０００３】そして、ハンマー２には力検出用ロードセ
ル３を取り付け、また該被試験物１の任意の場所に該ハ
ンマー２の加振による振動応答の３次元（ｘ，ｙ，ｚ）
方向加速度を検出するための加速度センサ４を取り付け
る。

【０００４】そして、加速度センサ４の出力信号とロー
ドセル３の出力信号とを演算装置５に与えてＦＦＴ処理
を行い周波数応答関数（以下、ＦＲＦ（コンプライアン
ス）と略称することがある）データを求めモード特性を
計算して出力する構成を有している。

【０００５】この場合、加振実験は通常、ハンマー２に
よる被試験物１の加振場所を固定して行い、加速度セン
サ４は逐次移動させて複数の応答関数を演算装置５に与
えるようにする。なお、被試験物１は柔らかいバネ（図
示せず）等により固定されている。

【０００６】このような振動特性解析装置の演算装置
（実験モード解析装置）５のアルゴリズムとしては今ま
で多くの手法が提案されているが、この中で周波数領域
における『偏分反復法』は精度の良さで現在広く使用さ
れている。

【０００７】しかし、上記の改良反復法には次の２つの
欠点がある。特に応答関数が複数の場合、求めるパラメータ｛γ｝
も応答関数の増加に応じて増加するため、パラメータ
｛γ｝を修正する修正ベクトル｛Δγ｝を求める式の中
の逆行列の計算時間が急増してしまう。変数が増加するため反復計算が収束し難く発散し易く
なり、途中で計算不能になることがある。

【０００８】そこで本発明者は特開平８−５４５０号公
報において、発散し易く計算時間がかかるというアルゴ
リズムの本質に関わる問題点の改良を実現することを目
的として、演算装置５が、実測データとモード特性パラ
メータによる理論値との残差が最小になるまで該パラメ
ータを繰り返し変えて行く偏分反復法を線形項と非線形
項とに分けて計算すると共に該線形項及び非線形項が変
化しても該残差が発散しない条件を設けることにより該
モード特性を同定することを提案している。

【０００９】この改良型偏分反復法について以下に簡潔
に説明する。求めるモード特性を｛γ｝で表す。周波数
点ω＝ω_iにおける周波数応答関数ＦＲＦの理論式をｆ
_i、対応する実験データをｙ_iで表す。同時に参照する
ＦＲＦの数をｍ、周波数点の数をＮで表すと、ＦＲＦは
複素数なので全データ数は２Ｎｍとなる。

【００１０】ｙ_iに混入している誤差の分布が正規分布
とすると、最も確からしいモード特性は，次式の重み付
き残差二乗和Ｓを最小にする｛γ｝である。

【００１１】

【数１】

【００１２】ここで、重み係数Ｗ_iは、ｙ_iの確率密度
分布を表す式における母分散の逆数（１／σ²）に比例
することが知られている。

【００１３】全変数｛γ｝を線形項｛α｝（元数ｕ）と
非線形項｛β｝（元数ｖ）に分ける。｛α｝を消去する
考え方を示す。全変数空間のうちで｛β｝のみを変数と
して、次式が成立する領域を探索する。

【００１４】

【数２】

【００１５】式（２）が成立する｛α｝の極小領域内に
おいて、ＦＲＦの理論式ｆ_i(｛α｝,｛β｝)の変数
｛α｝を｛β｝で表したものをＨ_iとすれば、次式のよ
うになる。

【００１６】

【数３】

【００１７】ここで、式（１）より次式が得られる。

【００１８】

【数４】

【００１９】式（３）の微分係数が決まれば、次式の線
形近似による最小二乗法で修正ベクトル｛Δβ｝を求め
ることができる。

【００２０】

【数５】［Ｃ］^T［Ｗ］［Ｃ］｛Δβ｝＝［Ｃ］^T［Ｗ］｛ｙ_i−ｆ_i｝・・・式（５）［Ｃ］：∂Ｈ_i／∂β_kをｉｋ要素とする行列［Ｗ］：Ｗ_iを対角成分とする対角行列

【００２１】式（５）における行列［Ｃ］はヤコビアン
行列と呼ばれている。このヤコビアン行列［Ｃ］の大き
さはｖ×２Ｎｍであり、ＦＲＦの数ｍと共に増加する。

【００２２】ここで、近年、多チャンネルデータ処理装
置や多軸速度センサの発達に伴いこの演算装置４は、短
時間で多数の周波数応答関数（ＦＲＦ）を実験で得るこ
とが可能となった。

【００２３】例えば図１では、３軸の加速度センサ４か
らのＸ，Ｙ，Ｚの３チャンネルの情報及び加振入力を検
出する力検出用ロードセル３の１チャンネルの情報から
３個のＦＲＦを同時に得ることができる。測定時間等を
短縮するために加速度センサ４を複数並列に設けると、
その増加分だけチャンネルが増える。

【００２４】このＦＲＦは、図１６に示すように、縦軸
を周波数、横軸を各測定ポイント（点１，点２・・・）
の軸データＸ，Ｙ，Ｚを順次並べた２Ｎ行ｍ列の行列と
して扱うことができ、例えば１データに４バイトのメモ
リ容量が割り当てられる。

【００２５】しかし、参照される周波数応答関数の数ｍ
の増加に伴い、これらの情報は増加するので、計算機の
限られた計算速度及びメモリ容量や従来の演算手法で
は、迅速な結果を得ることが困難となりつつある。

【００２６】この問題を解消するためには、同時に処理
すべきデータはＦＲＦ行列として与えられるので、この
行列を縮小すればよい。

【００２７】このＦＲＦ行列を縮小する方法としては、
従来より特異値分解法が提案されている。

【００２８】この特異値分解法において、実験で得られ
たＦＲＦ行列［Ｐ］（２Ｎ行ｍ列）を特異値分解を利用
して縮小するときの手順を図１７（１）〜（４）により
説明する。

【００２９】行列［Ｂ］（ｍ行ｍ列）を次式で計算す
る（同図（１））。

【００３０】

【数６】［Ｂ］＝［Ｐ］^T［Ｐ］・・・式（６）［Ｐ］^T：行列［Ｐ］の転置行列

【００３１】行列［Ｂ］の固有値λと固有ベクトル
｛ｘ｝を求める（同図（２））。固有値λと固有ベク
トル｛ｘ｝の組を、固有値λの大きい順に並べる。そし
て、固有値λが他と比較して小さいものを省略する。行
列［Ｂ］の固有値はｍ個求まるが、ここで、モード数ｎ
個の固有値と固有ベクトル｛ｘ_i｝の組に縮小する（同
図（３））。

【００３２】ｎ個の固有ベクトル｛ｘ_i｝から次式で
ｎ個に縮小されたＦＲＦベクトル｛Ｐ_i｝を計算する
（同図（４））。

【００３３】

【数７】｛Ｐ_i｝＝［Ｐ］｛ｘ_i｝ｉ＝１〜ｎ・・・式（７）

【００３４】従って、この方法は式（６），（７）と
の固有値、固有ベクトルの計算に時間がかかる。また式
（６），（７）を計算するために行列［Ｐ］（２Ｎ行ｍ
列）と行列［Ｂ］（ｍ行ｍ列）を記憶する必要がある。

【００３５】また、固有値を演算し、その固有値及び固
有ベクトルの大きい順に並べ、その中で係数の小さいも
のを切り捨ててモード数ｎ個のデータに圧縮するので、
圧縮情報は物理的な意味が不明であった。

【００３６】そこで、本発明者は特願平８−６８２５１
号公報において、被試験物を加振するためのインパルス
ハンマーに取り付けられた力検出用ロードセルと、該被
試験物の任意の場所に取り付けられて該被試験物の該イ
ンパルスハンマーに起因する加振応答を測定する加速度
センサと、該ロードセルの出力信号及び該センサの出力
信号とを受けて周波数応答関数データを求め、該周波数
応答関数データと該周波数応答関数の理論値の重み付き
残差二乗和を最小とする最小二乗法によりモード特性パ
ラメータを計算することにより該被試験物のモード特性
を同定する演算装置と、を備えた振動特性解析装置にお
いて、特異値分解法に代え、正確に情報を圧縮すること
を目的として、該演算装置が、該周波数応答関数の選択
された固有モードに対応する虚数部の極性に＋１又は−
１をそれぞれ乗じて正極性にした周波数応答関数を重ね
合わせることにより該固有モードに対応する周波数応答
関数をそれぞれ求め、該周波数応答関数の数を該固有モ
ードの数に縮小することを提案した。

【００３７】以下に、この方法の物理的な意味を説明す
る。まず、各点のＦＲＦを列ベクトルとして並べたＦＲ
Ｆ行列［Ｐ］を次式で縮小する。

【００３８】

【数８】［Ｑ］＝［Ｐ］［Ｋ］・・・式（８）ただし、［Ｐ］；実験で得られたＦＲＦ行列２Ｎ行ｍ
列［Ｑ］；縮小されたＦＲＦ行列２Ｎ行ｎ列［Ｋ］；縮小のための行列ｍ行ｎ列ｎ；採用固有モード数（５〜１０）（ｍ＞ｎ）

【００３９】行列［Ｋ］の作り方は、まず行列［Ｐ］
（２Ｎ行ｍ列）のうち各固有モードの共振峰に対応する
虚数部分だけを取り出して行方向に並べた行列を作成す
る。次にこの行列の各成分の正負符号は同じで大きさだ
けを１に変える（図２（１）参照）。

【００４０】このようにして作成した行列［Ｋ］^Tを転
置して行列［Ｋ］を作成する（同図（２））。

【００４１】それから、式（８）によりＦＲＦ行列
［Ｐ］と行列［Ｋ］に基づいて縮小されたＦＲＦ行列
［Ｑ］を作成する（同図（３））。

【００４２】式（８）の計算は、ある固有モードの共振
峰付近の虚数部分の符号が同一になるように全ＦＲＦを
重ね合わせることになる。

【００４３】この意味は、ＦＲＦ行列を測定する時の加
振点を応答点（測定点）として、ある固有モードが最も
励起される様に各点の力の方向を調整して、同時に衝撃
加振する時のＦＲＦを測定することに対応する。従っ
て、固有モードの数に等しい数のＦＲＦに、元のＦＲＦ
行列を縮小できる。

【００４４】実際に、あるモードが最も励起される様に
各点を正しく加振することは大変困難である。しかし、
各点毎の打撃試験等でＦＲＦを多数収得しておけば、式
（８）を用いて同時衝撃加振に相当する測定結果を容易
に得ることができる。

【００４５】式（８）の処理は各ＦＲＦに−１または＋
１を掛けて加算することである。加算係数の絶対値が１
なので、加算後のＦＲＦにおいては不規則雑音は相殺さ
れて減少する。

【００４６】いかにＦＲＦの数が多くても、式（８）に
より、抽出する固有モードの数と同一の本数にＦＲＦを
縮小できる。そのための計算時間と記憶容量は特異値分
解、または主成分分析による方法よりはるかに少ない。
また全ＦＲＦを加算しているため、不規則雑音は相殺さ
れ、データのばらつきが低減される。

【００４７】

【発明が解決しようとする課題】このように本発明者
は、多数の周波数応答関数（ＦＲＦ）を同時に用いてモ
ード特性を同定する多点参照を効率よく行うために、入
力データの元の情報を保持したままでＦＲＦを縮小する
改良された方法を提案したが、統計学的に最も確からし
いモード特性を同定する最尤推定法を実行するために
は、上記の式（１）に示す如く、曲線適合における重み
関数として分散の逆数を用いる必要がある。

【００４８】ＦＲＦを縮小しない場合には、元のＦＲＦ
をそのまま用いることにより、各ＦＲＦ毎に分散を求め
ればよかったが、ＦＲＦを縮小した入力データを用いる
場合には、この方法によって分散を求めることはできな
いから、重み関数を求めることが出来ない。

【００４９】そこで、上記の方法においては、経験的に
求めた値（例えば、ＦＲＦの振幅の逆数の平方を用いた
値）を重み関数としていたが、理論的な根拠に欠けるた
め、その正確性を保証することができなかった。

【００５０】したがって、本発明は、被試験物を加振す
るためのインパルスハンマーに取り付けられた力検出用
ロードセルと、該被試験物の任意の場所に取り付けられ
て該被試験物の該インパルスハンマーに起因する加振応
答を測定する加速度センサと、該ロードセルの出力信号
及び該センサの出力信号とを受けて周波数応答関数デー
タを求め、該周波数応答関数データと該周波数応答関数
の理論値との重み付き残差二乗和を最小とする最小二乗
法によりモード特性パラメータを計算することにより該
被試験物のモード特性を同定する演算装置とを備え、該
演算装置が、該周波数応答関数の選択された固有モード
に対応する虚数部の極性に＋１又は−１をそれぞれ乗じ
て正極性にした周波数応答関数を重ね合わせることによ
り該固有モードに対応する周波数応答関数をそれぞれ求
め、該周波数応答関数の数を該固有モードの数に縮小す
る振動特性解析装置において、周波数応答関数の縮小に
対応できる重み関数を理論的に導出することを目的とす
る。

【００５１】

【課題を解決するための手段】上記の目的を達成するた
め、本発明では、重み関数行列を統計理論に従って導出
する。最尤推定法では、重み関数として応答の母分散の
逆数を用いる必要がある。有限回の実験では母分散を正
確に求めることはできないので、近似的に標本分散を用
いる。ｑ回の反復実験によって求めたＦＲＦを平均する
ことによって得られる平均化ＦＲＦの実部及び虚部は、
測定点をｊ、角振動数をω_iとすると、次式のように表
される。

【００５２】

【数９】

【数１０】

【００５３】入力に誤差が混入しない場合には、コヒー
レンス関数は次式で与えられる。

【００５４】

【数１１】

【００５５】また標本分散の実部と虚部は次式で表され
る。

【００５６】

【数１２】

【数１３】

【００５７】理想的には、コヒーレンス関数と平均化Ｆ
ＲＦから分散の実部と虚部を別々に導出できればよい
が、上記の情報量だけでは困難である。そこで、便宜的
に実部と虚部の重みを同一のものとすると、上記の式
（９）〜（１３）より次式が導かれる。

【００５８】

【数１４】

【００５９】したがって、式（１４）の逆数として次の
重み関数が求められる。

【００６０】

【数１５】

【００６１】次に多点参照の場合を述べる。実験で得ら
れるＦＲＦデータには多数の誤差要因があると考えられ
る。それらが各々同程度の大きさで、かつ互いに独立に
寄与しているとすると、誤差は正規分布に従うことが知
られている。また、統計理論より、ｍ本のＦＲＦの標本
和は正規分布に従う。

【００６２】上記の式（８）の計算では、各ＦＲＦの共
振峰に対応する点の虚部の符号が同一になるように加算
しているので、縮小化したＦＲＦデータの分散は式（１
４）より次式が得られる。

【００６３】

【数１６】

【００６４】よって、重み関数は次式のように導出でき
る。

【００６５】

【数１７】

【００６６】したがって本発明に係る振動特性解析装置
における演算装置は、さらに、多数の該周波数応答関数
を同時に考慮した重み関数を求めて該残差二乗和に適用
したものであり、この重み関数を用いれば、最小二乗法
は最も確からしいモード特性を推定する最尤推定法と同
一になる。

【００６７】

【発明の実施の形態】本発明に係る振動特性解析装置の
構成は上述した図１に示した従来例と同様のものを用い
ることができるが、従来例との差異は演算装置５におけ
る演算処理である。以下にこの演算処理について図３に
示したフローチャートに沿って説明する。

【００６８】＜パラメータの分割＞複数の応答関数を同
時に考慮する時、その応答関数の数をｍ、同定する固有
モードの数を上記と同様にｎとするとパラメータ｛γ｝
の数は２ｎ＋（２ｎ＋４）ｍである。そのパラメータ
｛γ｝のうち非線形項は２ｎ、線形項は（２ｎ＋４）ｍ
である。このパラメータ｛γ｝を次式のように線形項
｛α｝と非線形項｛β｝とに分割する。

【００６９】

【数１８】｛γ｝^T＝｛｛α｝^T,｛β｝^T｝｛α｝^T＝｛‥,Ｕ_r,Ｖ_r,‥,Ｃ,Ｄ,Ｅ,Ｆ,‥｝｛β｝^T＝｛‥,σ_r,ω_dr,‥｝・・・式（１８）

【００７０】ここで、線形項｛α｝の中のＵ_r，Ｖ_rは測
定点に対するｒ番目の留数の実部及び虚部を示し、１／
Ｃ，１／Ｄは注目していない固有モードの影響を剰余質
量として、１／Ｅ，１／Ｆは同様に剰余剛性として導入
したパラメータである。

【００７１】また、非線形項｛β｝の中のσ_rは各ピー
ク（ｎ個存在し得る山々）、つまり各モード特性におけ
る減衰率を示し、ω_drは各ピーク、つまり各モード特性
における減衰固有角振動数を示す。

【００７２】なお、ハンマー２で被試験物１を加振した
ときの演算装置５へ与えられるモード特性の応答関数
（理論値）は振幅と位相を持つので、次式のように実部
Ｇ_R（ω）と虚部Ｇ₁（ω）とから成る複素数で表される
理論モデルがよく知られている。

【００７３】

【数１９】

【００７４】この減衰率σ_rと固有振動数ω_drは応答点
（加速度センサ４の取付位置）が変わっても変化しな
い。つまり非線形項のパラメータ｛β｝の数２ｎは応答
関数の数ｍが増加しても山の位置（周波数）が変わらな
いため不変であることが特徴となっており、これにより
応答関数の数の影響を受けないことが分かる。

【００７５】そして、非線形項が求まれば線形項だけは
容易に求めることができる。

【００７６】そこで非線形項だけを先ず求めることを工
夫する。このため、まず式（１）における重み係数Ｗ_i
を式（１５）により計算するとともに非線形パラメータ
｛β｝としての減衰固有角振動数ω_drとモード減衰率σ
_rの初期値が設定される（図３のステップＳ１）。

【００７７】減衰の初期値はどの様な値でも収束するの
で、ＦＲＦの振幅を二乗和した関数のピークから減衰固
有角振動数の初期値だけを得てもよい。

【００７８】そして、実験データである周波数応答関数
ｙ_iを入力する（同ステップＳ２）。

【００７９】次に、以下の計算負荷を減少させるため
に、式（８）を用いて入力した周波数応答関数ＦＲＦを
縮小する（同ステップＳ３）。縮小率はｎ（採用モード
数）／ｍ（ＦＲＦの数）となる。

【００８０】非線形項パラメータ｛β｝が初期設定され
れば、次式により線形項パラメータ｛α｝を１回の計算
で求めることができる（同ステップＳ４）。

【００８１】

【数２０】｛α｝＝([Ｄ]^T[Ｗ][Ｄ])^-1[Ｄ]^T[Ｗ]｛ｙ_i｝・・・式（２０）ただし、Ｄ_ij＝∂ｆ_i／∂α_j （ｉ＝1〜2p,ｊ＝1〜(2
n+4)m)

【００８２】このようにして求めた線形項｛α｝と非線
形項｛β｝とを用いて式（１）により周波数応答関数Ｆ
ＲＦの理論値ｆ_iを計算する（同ステップＳ５）。な
お、最初の非線形項の値はステップＳ１で初期設定され
た値である。

【００８３】このようにして計算した応答関数の理論値
と実験値との残差二乗和Ｓを式（１）に従って計算する
（同ステップＳ６）。

【００８４】そして、この残差二乗和Ｓが前回求めた値
より小さくなっているか否かを判定する（同ステップＳ
７）。なお、この場合の最初の残差二乗和Ｓは適当な値
を設定しておけばよい。

【００８５】この結果、残差二乗和Ｓが前回の値より小
さくなっているときには残差二乗和Ｓはまだ最小値に達
していないことになるので、この残差二乗和Ｓを最小値
にするための次の非線形項パラメータ｛β｝を以下に述
べるとおり求める。

【００８６】ここで、上記の式（５）は上述の如くＦＲ
Ｆの数ｍと共に増加するので計算機の記憶容量と計算時
間の両方に対して大きな負荷となっており、このような
負荷を軽減するために、同式の計算方法を次のように改
良することができる。

【００８７】式（５）を解く方法としては、（１）正規
方程式と呼ばれるｖ元ｖ行の連立方程式として解く方法
と、（２）ヤコビアン行列［Ｃ］をハウスホールダ法等
により直交変換して解く方法とがこれまでに知られてお
り、一般には後者（２）の方が計算精度が良いと言われ
ている。

【００８８】実際に計算してみると、倍精度におけるこ
の場合の計算では両者の計算結果は同一であった。しか
も計算時間と記憶容量に関しては正規方程式を解く前者
（１）の方が有利であった。

【００８９】そこで本実施例では、正規方程式を用いて
計算方法を改良することに着目した。まず、式（５）の
正規方程式は次式のように書き直すことができる。

【００９０】

【数２１】［Ｇ］｛Δβ｝＝｛ｄ｝・・・式（２１）

【００９１】式（２１）における［Ｇ］，｛ｄ｝をヤコ
ビアン行列［Ｃ］から算出する際に、周波数応答関数の
理論式Ｈ_iの微分係数を、ニュートン法で使用する残差
二乗和Ｓの微分係数に変換する。

【００９２】式（５）に式（３）を代入すると、式（２
１）の右辺は次式のようになる。

【００９３】

【数２２】

【００９４】この式（２２）を展開して整理すると次式
のようになる。

【００９５】

【数２３】

【００９６】｛α｝は、繰り返し計算におけるある段階
において極小条件を満足するので、∂Ｓ／∂α_s＝０
（ｓ＝１〜ｕ）が成立する。従って、次式が得られる。

【００９７】

【数２４】

【００９８】次に、式（４）において２次微分項を省略
した量を次式で表すこととする。

【００９９】

【数２５】

【０１００】式（２１）の左辺の係数行列［Ｇ］の成分
Ｇ_jkは式（５）より次式のように求められる。

【０１０１】

【数２６】

【０１０２】式（３）を式（２６）に代入すれば、次式
が得られる。

【０１０３】

【数２７】

【０１０４】式（２７）の右辺＝Ａ＋Ｂ＋Ｃと表し、式
（１）と式（２１）を代入して整理すれば、Ａ，Ｂ，Ｃ
に関する次式がそれぞれ得られる。

【０１０５】

【数２８】

【０１０６】

【数２９】

【０１０７】

【数３０】

【０１０８】ここで、｛ｘ｝^T，｛ａ｝，｛ｂ｝を次式
のように定義する。

【０１０９】

【数３１】

【０１１０】式（３０）に式（３１）を代入すると次式
が成立する。

【０１１１】

【数３２】

【０１１２】従って、式（２８），（２９），（３２）
をＧ_ik＝Ａ＋Ｂ＋Ｃに代入すると次式が得られる。

【０１１３】

【数３３】

【０１１４】今までは、式（３）の右辺第２項の行列
［Ｆ］^-1の乗算回数が多く、計算時間に対する負担にな
っていた。しかも、この右辺第２項の計算をｖ×２Ｎｍ
回行っていた。

【０１１５】ところが、上記の式（３３）においてはＧ
_jkを計算するためにΣが除去されているので必要なこの
計算は３回でよく、［Ｇ］全体でも３ｖ²回となる。

【０１１６】非線形項は全ＦＲＦの共通項であり、その
数ｖは通常周波数点Ｎや参照ＦＲＦ数ｍと比べて小さい
ので、計算時間が短縮できる。また、今までは２Ｎｍ行
ｖ列の行列［Ｃ］をそのまま記憶していたが、ｖ行ｖ列
の［Ｇ］とｖ行の｛ｄ｝を記憶すればよくなり、必要な
記憶容量が減少する。

【０１１７】即ち、式（３３）を計算するためには、Ｓ
とＳ’のα_sとβ_kに関する微分係数を計算すればよいこ
とになる。

【０１１８】このように、式（２４）を用いてｄ_kを計
算し、式（４），（２５）を用いてＳとＳ’のαとβに
関する微分係数を計算する。そして、Ｇ_jkを計算する
（同ステップＳ８）。

【０１１９】以上の結果を式（２１）に代入して修正ベ
クトル｛Δβ｝を計算で求める（同ステップＳ９）。

【０１２０】＜側面制約の設定＞以上で非線形項の修正
ベクトル｛Δβ｝が求まったが、このまま使用すると発
散する場合があるので、更に側面制約を設定する。

【０１２１】発散防止の側面制約としては修正ベクトル
に縮小因子をかけてその方向は同じ大きさを縮小する方
法があるが、この場合、モードが異なれば非線形項相互
の影響は少なくなるので、非線形項各々に次式の制約を
設定する。

【０１２２】

【数３４】｜△ω_dr／ω_dr｜≦ Ｋ₁ ｜△σ_r／σ_r｜≦ Ｋ₂， σ_r＞０・・・式（３４）

【０１２３】｛Δβ｝がこの条件に合っていない場合は
式（２５）により修正する（同ステップＳ１１）。

【０１２４】このようにして求められた修正ベクトル
｛Δβ｝と前回の非線形項パラメータ｛β｝（最初は初
期値）とを加算して今回の新しい非線形項パラメータ
｛β^*｝を計算する（同ステップＳ１２）。

【０１２５】

【数３５】｛β^*｝＝｛β｝＋｛Δβ｝・・・式（３５）

【０１２６】このようにして求めた新しい非線形項パラ
メータ｛β^*｝を用いて上記のステップＳ４以降を繰り
返し実行する。この繰り返しループのステップＳ６で計
算した残差二乗和Ｓが前回の値より小さくなっていない
ときには、このＳが最小値に達したことになるので、
｛β｝の計算ループから出る（同ステップＳ７の“Ｎ
ｏ”）。

【０１２７】そして、再度縮小しないＦＲＦを入力して
（同ステップＳ１３）、このＦＲＦと既に決定された
｛β｝を式（１９）に代入して線形項｛α｝を計算する
（同ステップＳ１４）。

【０１２８】この過程は線形最小二乗法が適用できるた
め繰り返し計算は不要である。

【０１２９】そして、最適なモード特性を示す｛α｝と
｛β｝を出力する（同ステップＳ１５）。

【０１３０】〔モデルデータによる検証〕ここで、本発
明に係る振動特性解析装置をモデルデータにより以下に
検証する。

【０１３１】図４の一端を固定した１０自由度モデルに
おける各質点の質量と各質点間のバネ剛性の値を表１の
ように与える。比例粘性減衰を仮定し、減衰係数ｃ_i、
質量ｍ_i、バネ剛性ｋ_iの間に以下の関係があるとする。

【０１３２】

【表１】

【０１３３】

【数３６】ｃ_i＝αｍ_i＋βｋ_i （α＝0.15，β＝3.0×10^-6) ・・・式（３６）対象とする周波数域および周波数分解能Δｆは、0.1〜6
40.1(Hz),Δｆ＝0.1(Hz) とする。

【０１３４】この周波数領域には１次から６次までの固
有モードが存在する。加振点を質点１、応答点を質点１
から質点１０とする場合の１０個のＦＲＦをシミュレー
ションで求め、多点参照の入力データとする。

【０１３５】図５〜図１０は、質点１の自己ＦＲＦにお
ける同定結果を示している。まず、入力（加振力）には
雑音を入れず、応答に４％の正規分布雑音を混入させた
場合の同定を行ったときの入力データの位相と振幅とコ
ヒーレンス関数をそれぞれ図５（１）〜（３）に示す。
また、図６（１）及び（２）には、図５の入力データを
用いて従来より重み関数として一般的に用いられている
例えばＦＲＦの振幅の逆数の平方により同定した結果と
真値を位相とともに示す。図７（１）及び（２）には、
図５の入力データを用いて本発明による式（１５）の重
み関数によって同定した結果と真値を位相とともに示
す。

【０１３６】図６と図７を比較すると、入力データに雑
音を含まない時は両手法とも真値と同定結果の値とがほ
ぼ一致しており、真値を同定できることが分かる。

【０１３７】次に、加振力と応答のそれぞれに２％と４
％の正規分布に従った雑音を混入させて同定を行う。図
８（１）〜（３）には、同定に使用する入力データの位
相と振幅とコヒーレンス関数を示す。図９（１）及び
（２）には、図８の入力データを用い従来の重み関数に
より行った同定を行ったときの結果と真値を位相ととも
に示す。図１０（１）及び（２）には、図８の入力デー
タを用い本発明による重み関数で同定を行ったときの結
果と真値を位相とともに示す。

【０１３８】図９と図１０を比較すると、後者の本発明
による重み関数によって同定を行った場合には同定結果
と真値とがほぼ一致しているのに対し、前者の場合には
大きな差異が認められる。したがって、本発明の同定精
度の方が従来の重み関数を用いた方より明らかに優れて
いることが分かる。

【０１３９】このように、従来の重み関数による同定で
は、入力に雑音を含まない場合には同定精度は良いが、
入力に混入する雑音が大きくなるに従って同定精度は低
下する。それに対し、本発明による重み関数による同定
は全周波数応答関数より重み関数を求めているので偶然
誤差が相殺され、入力に雑音が混入していても同定精度
は良い。

【０１４０】〔簡略化ベアリングビームの同定結果〕図
１１に示す軟鋼製角棒を対象に振動試験を行う。周波数
応答関数の振幅の逆数の平方を用いる従来の重み関数
と、式（１５）に示した本発明による重み関数を用いて
同定を行い、両結果を比較する。

【０１４１】ただし、ここでは、３番の点をｙ方向に加
振し、応答は１番から１４番までを採用する。

【０１４２】図１２には、３番点の自己応答の入力デー
タの位相（同図（１））と、従来法及び本発明の振幅
（同図（２））を示す。図示の如く、従来法に比べて本
発明による同定結果の方が精度が良いことが分かる。

【０１４３】また、図１３〜図１５には、それぞれ、１
次固有モード（Ｚ方向の曲げ１次）と、３次固有モード
（Ｚ方向の曲げ３次）と、５次固有モード（Ｚ方向の曲
げ３次）のモード形状の同定結果が示されている。

【０１４４】図１３及び図１４では、両重み関数による
同定結果は共に真のモード形状を同定できているが、図
１５では従来の重み関数に比べ、本発明による重み関数
で同定した結果の方が３次の対称モード形（中心位置は
３９１mm）に近いことから真値に近いと考えられる。

【０１４５】これは本発明の手法が全ＦＲＦとコヒーレ
ンス関数より重み関数を求めているため、入力データに
含まれる雑音の影響が相殺され、全固有モードで精度の
良い同定結果が得られたと考えられる。

【０１４６】

【発明の効果】以上説明したように本発明に係る振動特
性解析装置によれば、周波数応答関数データと該周波数
応答関数の理論値との重み付き残差二乗和を最小とする
最小二乗法によりモード特性パラメータを計算すること
により該被試験物のモード特性を同定する演算装置が、
該周波数応答関数の選択された固有モードに対応する虚
数部の極性に＋１又は−１をそれぞれ乗じて正極性にし
た周波数応答関数を重ね合わせることにより該固有モー
ドに対応する周波数応答関数をそれぞれ求め、該周波数
応答関数の数を該固有モードの数に縮小するとともに、
多数の該周波数応答関数を同時に考慮した重み関数を求
めて該残差二乗和に適用するように構成したので、入力
に雑音が混入しても、従来の重み関数に比べ、同定精度
が向上する。

【０１４７】また、本発明による重み関数は全周波数応
答関数で共通であるため、多点参照において入力データ
が増加する場合に、計算機の記憶容量に対する負荷を小
さくすることができる。

【図面の簡単な説明】

【図１】本発明及び従来例に共通な振動特性解析装置を
示したブロック図である。

【図２】本発明に係る振動特性解析装置で用いられるＦ
ＲＦ（周波数応答関数）行列のデータ縮小を説明する図
である。

【図３】本発明に係る振動特性解析装置における演算装
置に格納され且つ実行される演算アルゴリズムを示した
フローチャート図である。

【図４】一端を固定した１０自由度モデルを示した図で
ある。

【図５】同定のために用いる入力ＦＲＦデータ（その
１：入力雑音を含まず応答に４％の正規分布雑音を混入
させたもの）の位相、振幅、及びコヒーレンス関数を示
したグラフ図である。

【図６】図５に示した入力データ（その１）に対して従
来の重み関数を用いて同定したときの位相と振幅を示し
たグラフ図である。

【図７】図５に示した入力データ（その１）に対して本
発明の重み関数を用いて同定したときの位相と振幅を示
したグラフ図である。

【図８】同定のために用いる入力ＦＲＦデータ（その
２：入力と応答にそれぞれ２％及び４％の正規分布雑音
を混入させたもの）の位相、振幅、及びコヒーレンス関
数を示したグラフ図である。

【図９】図８に示した入力データ（その２）に対して従
来の重み関数を用いて同定したときの位相と振幅を示し
たグラフ図である。

【図１０】図８に示した入力データ（その２）に対して
本発明の重み関数を用いて同定したときの位相と振幅を
示したグラフ図である。

【図１１】ベアリングビームを模擬した振動試験対象の
軟鋼製角棒を示した図である。

【図１２】図１１の角棒における３番点の自己応答の入
力データの位相及び振幅を示したグラフ図である。

【図１３】図１１の角棒を加振した時の１次固有モード
（Ｙ方向の曲げ１次）のモード形状の同定結果を示した
グラフ図である。

【図１４】図１１の角棒を加振した時の２次固有モード
（Ｙ方向の曲げ２次）のモード形状の同定結果を示した
グラフ図である。

【図１５】図１１の角棒を加振した時の３次固有モード
（Ｙ方向の曲げ３次）のモード形状の同定結果を示した
グラフ図である。

【図１６】ＦＲＦ行列の例を示した図である。

【図１７】従来の特異分解法によるデータ縮小例を示し
た図である。

【符号の説明】

１被試験物２インパルスハンマー３ロードセル４加速度センサ５演算装置図中、同一符号は同一または相当部分を示す。

Claims

【特許請求の範囲】

【請求項１】被試験物を加振するためのインパルスハン
マーに取り付けられた力検出用ロードセルと、該被試験
物の任意の場所に取り付けられて該被試験物の該インパ
ルスハンマーに起因する加振応答を測定する加速度セン
サと、該ロードセルの出力信号及び該センサの出力信号
とを受けて周波数応答関数データを求め、該周波数応答
関数データと該周波数応答関数の理論値との重み付き残
差二乗和を最小とする最小二乗法によりモード特性パラ
メータを計算することにより該被試験物のモード特性を
同定する演算装置とを備え、該演算装置が、該周波数応
答関数の選択された固有モードに対応する虚数部の極性
に＋１又は−１をそれぞれ乗じて正極性にした周波数応
答関数を重ね合わせることにより該固有モードに対応す
る周波数応答関数をそれぞれ求め、該周波数応答関数の
数を該固有モードの数に縮小する振動特性解析装置にお
いて、該演算装置が、さらに、多数の該周波数応答関数を同時
に考慮した重み関数を求めて該残差二乗和に適用するこ
とを特徴とした振動特性解析装置。