(57)【特許請求の範囲】
【請求項１】
被試験物を加振するためのインパルスハンマーに取り付けられた力検出用ロードセルと、該被試験物の任意の場所に取り付けられて該被試験物の該インパルスハンマーの加振による応答を測定する加速度センサと、該ロードセルの出力信号及び該センサの出力信号とを受けて周波数応答関数データを求め、該周波数応答関数データと該周波数応答関数の理論式の残差二乗和を最小とする最小二乗法において線形項と非線形項に分けてモード特性パラメータを計算することにより該被試験物のモード特性を同定する演算装置と、を備えた振動特性解析装置において、
該演算装置が、該周波数応答関数の選択された固有モードに対応する虚数部の極性に＋１又は−１を乗じて正極性にした周波数応答関数を重ね合わせることにより該固有モードに対応する該周波数応答関数を求め、該周波数応答関数の数を該固有モードの数に縮小することを特徴とする振動特性解析装置。
【請求項２】
請求項１において、
該演算装置が、線形近似による最小二乗法で修正ベクトルを求める正規方程式における正規行列を周波数応答関数の理論値と実験値との重み付き残差二乗和及びその２次微分項を省いた残差二乗和の線形項及び非線形項に関する微分係数で求めたことを特徴とする振動特性解析装置。
【発明の詳細な説明】
【発明の属する技術分野】
本発明は振動特性解析装置に関し、特にエンジンブロック等の被試験物を加振することにより応答を得てその応答関数から被試験物のモード特性を同定する振動特性解析装置に関するものである。
【０００１】
実際の機械構造物の動特性を解析するためには、該構造物を加振試験することにより得られた振動応答特性から該構造物の動特性を同定する必要があるが、この場合、振動現象を直接表現するパラメータであり、現象を直接理解し易く、収束性が良好なモード特性（固有振動数、モード減衰比、固有モード形状などのモーダルパラメータ）を同定する手法が現在では主流となっている。
【０００２】
【従来の技術】
図１は従来より知られている振動特性解析装置の構成を概略的に示したもので、まず被試験物１とこの被試験物１を加振するためのインパルスハンマー２とを用意する。
【０００３】
そして、ハンマー２には力検出用ロードセル３を取り付け、また該被試験物１の任意の場所に該ハンマー２の加振による振動応答の３次元（ｘ，ｙ，ｚ）方向加速度を検出するための加速度センサ４を取り付ける。
【０００４】
そして、加速度センサ４の出力信号とロードセル３の出力信号とを演算装置５に与えてＦＦＴ処理を行い周波数応答関数データ（ＦＲＦ：コンプライアンス）を求めモード特性を計算して出力する構成を有している。
【０００５】
この場合、加振実験は通常、ハンマー２による被試験物１の加振場所を固定して行い、加速度センサ４は逐次移動させて複数の応答関数を演算装置５に与えるようにする。なお、被試験物１は柔らかいバネ（図示せず）等により固定されている。
【０００６】
このような振動特性解析装置の演算装置（実験モード解析装置）５のアルゴリズムとしては今まで多くの手法が提案されているが、この中で周波数領域における『偏分反復法』は精度の良さで現在広く使用されている。
【０００７】
しかし、上記の偏分反復法には次の２つの欠点がある。
【０００８】
▲１▼特に応答関数が複数の場合、求めるパラメータ｛γ｝も応答関数の増加に応じて増加するため、パラメータ｛γ｝を修正する修正ベクトル｛Δγ｝を求める式の中の逆行列の計算時間が急増してしまう。
【０００９】
▲２▼変数が増加するため反復計算が収束し難く発散し易くなり、途中で計算不能になることがある。
【００１０】
そこで本発明者は特開平８−５４５０号公報において、発散し易く計算時間がかかるというアルゴリズムの本質に係わる問題点の改良を実現することを目的として、演算装置５が、実測データとモード特性パラメータによる理論値との残差が最小になるまで該パラメータを繰り返し変えて行く偏分反復法を線形項と非線形項とに分けて計算すると共に該線形項及び非線形項が変化しても該残差が発散しない条件を設けることにより該モード特性を同定することを提案している。
【００１１】
この改良型偏分反復法について以下に簡潔に説明する。
【００１２】
求めるモード特性を｛γ｝で表す。周波数点ω＝ωｉ  における周波数応答関数ＦＲＦの理論式をｆｉ  、対応する実験データをｙｉ  で表す。同時に参照するＦＲＦの数をｍ、周波数点の数をＮで表すと、ＦＲＦは複素数なので全データ数は２Ｎｍとなる。
【００１３】
ｙｉ  に混入している誤差の分布が正規分布とすると、最も確からしいモード特性は，次式の重み付き残差二乗和Ｓを最小にする｛γ｝である。
【００１４】
【数１】
0003599882.tif 000002
【００１５】
ここで、重み係数Ｗｉは、ｙｉ  の母分散の逆数に比例する。
【００１６】
全変数｛γ｝を線形項｛α｝（元数ｕ）と非線形項｛β｝（元数ｖ）に分ける。｛α｝を消去する考え方を示す。全変数空間のうちで｛β｝のみを変数として、次式が成立する領域を探索する。
【００１７】
【数２】
0003599882.tif 000003
【００１８】
式（２）が成立する｛α｝の極小領域内において、ＦＲＦの理論式ｆｉ（｛α｝，｛β｝）の変数｛α｝を｛β｝で表したものをＨｉとすれば、次式のようになる。
【００１９】
【数３】
0003599882.tif 000004
【００２０】
ここで、式（１）より次式が得られる。
【００２１】
【数４】
0003599882.tif 000005
【００２２】
式（３）の微分係数が決まれば、次式の線形近似による最小二乗法で修正ベクトル｛Δβ｝を求めることができる。
【００２３】
【数５】
［Ｃ］Ｔ［Ｗ］［Ｃ］｛Δβ｝＝［Ｃ］Ｔ［Ｗ］｛ｙｉ−ｆｉ｝・・・式（５）
［Ｃ］：∂Ｈｉ／∂βｋをｉｋ要素とする行列
［Ｗ］：Ｗｉを対角成分とする対角行列
【００２４】
式（５）における行列［Ｃ］はヤコビアン行列と呼ばれている。このヤコビアン行列［Ｃ］の大きさはｖ×２Ｎｍであり、ＦＲＦの数ｍと共に増加する。
【００２５】
【発明が解決しようとする課題】
ここで、近年、多チャンネルデータ処理装置や多軸速度センサの発達に伴い、この演算装置４は、短時間で多数の周波数応答関数（ＦＲＦ）を実験で得ることが可能となった。
【００２６】
例えば図１では、３軸の加速度センサ４からのＸ，Ｙ，Ｚの３チャンネルの情報及び加振入力を検出する力検出用ロードセル３の１チャンネルの情報＝合計４チャンネルの情報が入力される。測定時間等を短縮するために加速度センサ４を複数並列に設けると、その増加分だけチャンネルが増えることになる。
【００２７】
このＦＲＦは、図２０に示すように、縦軸を周波数、横軸を各測定ポイント（点１，点２・・・）の軸データＸ，Ｙ，Ｚを順次並べた２Ｎ行ｍ列の行列として扱うことができ、例えば１データに４バイトのメモリ容量が割り当てられる。
【００２８】
しかし、参照される周波数応答関数の数ｍの増加に伴い、これらの情報は増加するので、計算機の限られた計算速度及びメモリ容量や従来の演算手法では、迅速な結果を得ることが困難となりつつある。
【００２９】
この問題を解消するためには、同時に処理すべきデータはＦＲＦ行列として与えられるので、この行列を縮小すればよい。
【００３０】
このＦＲＦ行列を縮小する方法としては、従来より特異値分解法が提案されている。
【００３１】
この特異値分解法において、実験で得られたＦＲＦ行列［Ｐ］（２Ｎ行ｍ列）を特異値分解を利用して縮小するときの手順を図２１（１）〜（４）により説明する。
【００３２】
▲１▼行列［Ｂ］（ｍ行ｍ列）を次式で計算する（同図（１））。
【００３３】
【数６】
［Ｂ］＝［Ｐ］Ｔ［Ｐ］                              ・・・式（６）
［Ｐ］Ｔ：行列［Ｐ］の転置行列
【００３４】
▲２▼行列［Ｂ］の固有値λと固有ベクトル｛ｘ｝を求める（同図（２））。
【００３５】
▲３▼固有値λと固有ベクトル｛ｘ｝の組を、固有値λの大きい順に並べる。そして、固有値λが他と比較して小さいものを省略する。行列［Ｂ］の固有値はｍ個求まるが、ここで、モード数ｎ個の固有値と固有ベクトル｛ｘｉ｝の組に縮小する（同図（３））。
【００３６】
▲４▼ｎ個の固有ベクトル｛ｘｉ｝から次式でｎ個に縮小されたＦＲＦベクトル｛Ｐｉ｝を計算する（同図（４））。
【００３７】
【数７】
｛Ｐｉ｝＝［Ｐ］｛ｘｉ｝    ｉ＝１〜ｎ              ・・・式（７）
【００３８】
従って、この方法は式（６），（７）と▲２▼の固有値、固有ベクトルの計算に時間がかかる。また式（６），（７）を計算するために行列［Ｐ］（２Ｎ行ｍ列）と行列［Ｂ］（ｍ行ｍ列）を記憶する必要がある。
【００３９】
また、固有値を演算し、その固有値及び固有ベクトルの大きい順に並べ、その中で係数の小さいものを切り捨ててモード数ｎ個のデータに圧縮するので、圧縮情報は物理的な意味が不明であり、完全な情報として保存されず、従って正確なモード特性が得られないという問題点があった。
【００４０】
そこで、本発明は、被試験物を加振するためのインパルスハンマーに取り付けられた力検出用ロードセルと、該被試験物の任意の場所に取り付けられて該被試験物の該インパルスハンマーに起因する加振応答を測定する加速度センサと、該ロードセルの出力信号及び該センサの出力信号とを受けて周波数応答関数データを求め、該周波数応答関数データと該周波数応答関数の理論式の残差二乗和を最小とする最小二乗法において線形項と非線形項に分けてモード特性パラメータを計算することにより該被試験物のモード特性を同定する演算装置と、を備えた振動特性解析装置において、特異値分解法に代え、正確に扱う情報を圧縮することを目的とする。
【００４１】
【課題を解決するための手段】
〔１〕本発明に関し、入力データに含まれる情報を保存したままでＦＲＦを縮小するための、簡明で物理的意味が明確な方法を以下に説明する。
【００４２】
各点のＦＲＦを列ベクトルとして並べたＦＲＦ行列［Ｐ］を次式で縮小する。
【００４３】
【数８】
［Ｑ］＝［Ｐ］［Ｋ］                                ・・・式（８）
ただし、［Ｐ］；実験で得られたＦＲＦ行列  ２Ｎ行ｍ列
［Ｑ］；縮小されたＦＲＦ行列      ２Ｎ行ｎ列
［Ｋ］；縮小のための行列            ｍ行ｎ列
ｎ  ；採用固有モード数（５〜１０）（ｍ＞ｎ）
【００４４】
行列［Ｋ］の作り方は、まず行列［Ｐ］（２Ｎ行ｍ列）のうち各固有モードの共振峰に対応する虚数部分だけを取り出して行方向に並べた行列を作成する。次にこの行列の各成分の正負符号は同じで大きさだけを１に変える（図２（１））。
【００４５】
このようにして作成した行列［Ｋ］Ｔを転置して行列［Ｋ］を作成する（同図（２））。
【００４６】
それから、式（８）によりＦＲＦ行列［Ｐ］と行列［Ｋ］に基づいて縮小されたＦＲＦ行列［Ｑ］を作成する（同図（３））。
【００４７】
式（８）の計算は、ある固有モードの共振峰付近の虚数部分の符号が同一になるように全ＦＲＦを重ね合わせることになる。
【００４８】
この意味は、ＦＲＦ行列を測定する時の加振点を応答点（測定点）として、ある固有モードが最も励起される様に各点の力の方向を調整して、同時に衝撃加振する時のＦＲＦを測定することに対応する。従って、固有モードの数に等しい数のＦＲＦに、元のＦＲＦ行列を縮小できる。
【００４９】
実際に、あるモードが最も励起される様に各点を正しく加振することは大変困難である。しかし、各点毎の打撃試験等でＦＲＦを多数収得しておけば、式（８）を用いて同時衝撃加振に相当する測定結果を容易に得ることができる。
【００５０】
式（８）の処理は各ＦＲＦに−１または＋１を掛けて加算することである。加算係数の絶対値が１なので、加算後のＦＲＦにおいては不規則雑音は相殺されて減少する。
【００５１】
いかにＦＲＦの数が多くても、式（８）により、抽出する固有モードの数と同一の本数にＦＲＦを縮小できる。そのための計算時間と記憶容量は特異値分解、または主成分分析による方法よりはるかに少ない。
【００５２】
上記の原理を図２〜図７を参照してより具体的に説明する。
【００５３】
ＦＲＦ行列は、図２（１）に示すように、上半分を実部データ、下半分を虚部データとして配列したものである。この状態は、圧縮前の状態で、ｍ列分の情報量である。
【００５４】
これを採用固有モード数ｎ列の情報量に圧縮する手法について説明する。
【００５５】
採用固有モード数とは、図３（２）に示すように、ユーザーが当該技術分野における被試験物の観測したい周波数領域を区切った場合に、その範囲内における固有モード数のことである。
【００５６】
同図において、例えば０Ｈｚ近傍から１０００Ｈｚの範囲を観測範囲とすれば、固有モード数は固有モード▲１▼，▲２▼の２つとなる。一般には、５〜１０個である。
【００５７】
この固有モードに着目してＦＲＦの数を縮小する手順を図４〜図７を参照しつつ説明する。
【００５８】
図４（１）〜（４）において、被試験物１（図１参照）の各測定ポイント（１）〜（４）の波形を単純加算すると、例えば図５に示すように、固有モード▲１▼の形状が所望する共振峰を示す形状とならない。
【００５９】
そこで、図６（２）に示すように、位相が＋のものはそのままの振幅値で、位相が−のものはそれに−１を乗算した振幅値を用いて正極性とした後、これらを加算すると同図（１）に示すように所望の波形が得られる。
【００６０】
ただ、この操作では、固有モード▲１▼だけは所望する結果を得られるが、次の固有モード▲２▼以降はどんな波形になるか分からない。
【００６１】
そこで、次は、固有モード▲２▼の部分において同様な操作を行う（図７参照）。
【００６２】
以下順次、採用固有モード数分だけこの操作を実行する。
【００６３】
これにより、各固有モードに対して正確な情報で、全データ数ｍのＦＲＦが採用固有モード数ｎのＦＲＦに縮小される（ｎ＜ｍ）。
【００６４】
〔２〕一方、式（５）は上述の如くＦＲＦの数ｍと共に増加するので計算機の記憶容量と計算時間の両方に対して大きな負荷となっており、このような負荷を軽減するために、同式の計算方法を次のように改良する。
【００６５】
式（５）を解く方法としては、（１）正規方程式と呼ばれるｖ元ｖ行の連立方程式として解く方法と、（２）ヤコビアン行列［Ｃ］をハウスホールダ法等により直交変換して解く方法とがこれまでに知られており、一般には後者（２）の方が計算精度が良いと言われている。
【００６６】
実際に計算してみた結果、倍精度におけるこの場合の計算では両者の計算結果は同一であった。しかも計算時間と記憶容量に関しては正規方程式を解く前者（１）の方が有利であった。
【００６７】
そこで本発明においては、正規方程式を用いて計算方法を改良することに着目した。
【００６８】
まず、式（５）の正規方程式は次式のように書き直すことができる。
【００６９】
【数９】
［Ｇ］｛Δβ｝＝｛ｄ｝                              ・・・式（９）
【００７０】
式（９）における［Ｇ］，｛ｄ｝をヤコビアン行列［Ｃ］から算出する際に、周波数応答関数の理論式Ｈｉの微分係数を、ニュートン法で使用する残差二乗和Ｓの微分係数に変換する。
【００７１】
式（５）に式（３）を代入すると、式（９）の右辺は次式のようになる。
【００７２】
【数１０】
0003599882.tif 000006
【００７３】
式（１０）を展開して整理すると次式のようになる。
【００７４】
【数１１】
0003599882.tif 000007
【００７５】
｛α｝は、繰り返し計算におけるある段階において極小条件を満足するので、∂Ｓ／∂αｓ＝０（ｓ＝１〜ｕ）が成立する。従って、次式が得られる。
【００７６】
【数１２】
0003599882.tif 000008
【００７７】
次に、式（４）において２次微分項を省略した量を次式で表すこととする。
【００７８】
【数１３】
0003599882.tif 000009
【００７９】
式（９）の左辺の係数行列［Ｇ］の成分Ｇｊｋは式（５）より次式のように求められる。
【００８０】
【数１４】
0003599882.tif 000010
【００８１】
式（３）を式（１４）に代入すれば、次式が得られる。
【００８２】
【数１５】
0003599882.tif 000011
【００８３】
式（１５）の右辺＝Ａ＋Ｂ＋Ｃと表し、式（１）と式（９）を代入して整理すれば、Ａ，Ｂ，Ｃに関する次式がそれぞれ得られる。
【００８４】
【数１６】
0003599882.tif 000012
【００８５】
【数１７】
0003599882.tif 000013
【００８６】
【数１８】
0003599882.tif 000014
【００８７】
ここで、｛ｘ｝Ｔ，｛ａ｝，｛ｂ｝を次式のように定義する。
【００８８】
【数１９】
0003599882.tif 000015
【００８９】
式（１８）に式（１９）を代入すると次式が成立する。
【００９０】
【数２０】
0003599882.tif 000016
【００９１】
従って、式（１６），（１７），（２０）をＧｊｋ＝Ａ＋Ｂ＋Ｃに代入すると次式が得られる。
【００９２】
【数２１】
0003599882.tif 000017
【００９３】
今までは、式（３）の右辺第２項の行列［Ｆ］−１の乗算回数が多く、計算時間に対する負担になっていた。しかも、この右辺第２項の計算をｖ×２Ｎｍ回行っていた。
【００９４】
ところが、式（２１）においてはＧｊｋを計算するためにΣが除去されているので必要なこの計算は３回でよく、［Ｇ］全体でも３ｖ２回となる。
【００９５】
非線形項は全ＦＲＦの共通項であり、その数ｖは通常周波数点Ｎや参照ＦＲＦ数ｍと比べて小さいので、計算時間が短縮できる。また、今までは２Ｎｍ行ｖ列の行列［Ｃ］をそのまま記憶していたが、ｖ行ｖ列の［Ｇ］とｖ行の｛ｄ｝を記憶すればよくなり、必要な記憶容量が減少する。
【００９６】
即ち、式（２１）を計算するためには、ＳとＳ’のαｓとβｋに関する微分係数を計算すればよいことになる。
【００９７】
【発明の実施の形態】
本発明に係る振動特性解析装置の構成は上述した図１に示した従来例と同様のものを用いることができるが、従来例との差異は演算装置５における演算処理である。以下にこの演算処理について図８に示したフローチャートに沿って説明する。
【００９８】
＜パラメータの分割＞
複数の応答関数を同時に考慮する時、その応答関数の数をｍ、同定する固有モードの数を上記と同様にｎとするとパラメータ｛γ｝の数は２ｎ＋（２ｎ＋４）ｍである。そのパラメータ｛γ｝のうち非線形項は２ｎ、線形項は（２ｎ＋４）ｍである。このパラメータ｛γ｝を次式のように線形項｛α｝と非線形項｛β｝とに分割する。
【００９９】
【数２２】
｛γ｝Ｔ＝｛｛α｝Ｔ，｛β｝Ｔ｝
｛α｝Ｔ＝｛‥，Ｕｒｊ，Ｖｒｊ，‥，Ｃｊ，Ｄｊ，Ｅｊ，Ｆｊ，‥｝（ｒ＝１〜ｎ，ｊ＝１〜ｍ）
｛β｝Ｔ＝｛‥，σｒ，ωｄｒ，‥｝（ｒ＝１〜ｎ）                          ・・・式（２２）
【０１００】
ここで、線形項｛α｝の中のＵｒｊ，Ｖｒｊは測定点ｊに対するｒ番目の留数の実部及び虚部を示し、１／Ｃｊ，１／Ｄｊは注目していない固有モードの影響を剰余質量として、１／Ｅｊ，１／Ｆｊは同様に剰余剛性として導入したパラメータである。
【０１０１】
また、非線形項｛β｝の中のσｒは各ピーク（ｎ個存在し得る図９（２）に示す山々）、つまり各モード特性における減衰率を示し、ωｄｒは各ピーク、つまり各モード特性における減衰固有角振動数を示す。
【０１０２】
なお、ハンマー２で被試験物１を加振したときの演算装置５へ与えられるモード特性の応答関数（理論値）は振幅と位相を持つので、次式のように実部ＧＲｊ（ω）と虚部ＧＩｊ（ω）とから成る複素数で表される理論モデルがよく知られている。
【０１０３】
【数２３】
0003599882.tif 000018
【０１０４】
この減衰率σｒと固有振動数ωｄｒは応答点（加速度センサ４の取付位置）が変わっても変化しない。つまり非線形項のパラメータ｛β｝の数２ｎは応答関数の数ｍが増加しても図９（２）に示す山の位置（周波数）が変わらないため不変であることが特徴となっており、これにより応答関数の数の影響を受けないことが分かる。
【０１０５】
そして、非線形項が求まれば線形項は線形最小２乗法により、繰り返し計算なしで容易に求めることができる。
【０１０６】
そこで非線形項だけを先ず求めることを工夫する。このため、まず式（１）における重み係数Ｗｉ及び非線形項パラメータ｛β｝としての減衰固有角振動数ωｄｒとモード減衰率σｒの初期値が設定される（図８のステップＳ１）。
【０１０７】
減衰の初期値はどの様な値でも収束するので、ＦＲＦの振幅を二乗和した関数のピーク（図９（２））から減衰固有角振動数の初期値だけを得てもよい。
【０１０８】
そして、実験データである周波数応答関数ｙｉを入力する（ステップＳ２）。
【０１０９】
次に、以下の計算負荷を減少させるために、式（８）を用いて入力した周波数応答関数ＦＲＦを縮小する（ステップＳ３）。縮小率はｎ（採用モード数）／ｍ（ＦＲＦの数）となる。
【０１１０】
非線形項パラメータ｛β｝が初期設定されれば、次式により線形項パラメータ｛α｝を１回の計算で求めることができる（ステップＳ４）。
【０１１１】
【数２４】
｛α｝＝（［Ｄ］Ｔ［Ｗ］［Ｄ］）−１［Ｄ］Ｔ［Ｗ］｛ｙｉ｝    ・・・式（２４）
ただし、Ｄｉｊ＝∂ｆｉ／∂αｊ    （ｉ＝１〜２ｐ，ｊ＝１〜（２ｎ＋４）ｍ）
【０１１２】
このようにして求めた線形項｛α｝と非線形項｛β｝とを用いて式（２３）により周波数応答関数ＦＲＦの理論値ｆｉを計算する（ステップＳ５）。なお、最初の非線形項の値はステップＳ１で初期設定された値である。
【０１１３】
このようにして計算した応答関数の理論値と実験値との残差二乗和Ｓを式（１）に従って計算する（ステップＳ６）。
【０１１４】
そして、この残差二乗和Ｓが前回求めた値より小さくなっているか否かを判定する（ステップＳ７）。なお、この場合の最初の残差二乗和Ｓは充分大きな値を設定しておけばよい。
【０１１５】
この結果、残差二乗和Ｓが前回の値より小さくなっているときには残差二乗和Ｓはまだ最小値に達していないことになるので、この残差二乗和Ｓを最小値にするための次の非線形項パラメータ｛β｝を以下のとおり求める。
【０１１６】
即ち、式（１２）を用いてｄｋを計算し、式（４），（１３）を用いてＳとＳ’のαとβに関する微分係数を計算する。そして、式（２１）を用いてＧｊｋを計算する（ステップＳ８）。
【０１１７】
以上の結果を式（９）に代入して修正ベクトル｛Δβ｝を計算で求める（ステップＳ９）。
【０１１８】
＜側面制約の設定＞
以上で非線形項の修正ベクトル｛Δβ｝が求まったが、このまま使用すると発散する場合があるので、更に側面制約を設定する。
【０１１９】
発散防止の側面制約としては修正ベクトルに縮小因子をかけてその方向は同じ大きさを縮小する方法がある。この方法では各パラメータに同一の係数をかけることになる。ここでは、非線形項各々に別々の次式の制約を設定する。
【０１２０】
【数２５】
｜△ωｄｒ／ωｄｒ｜≦  Ｋ１  
｜△σｒ  ／σｒ  ｜≦  Ｋ２  ，  σｒ  ＞０                            ・・・式（２５）
【０１２１】
｛Δβ｝がこの条件に合っていない場合は式（２５）により修正する（ステップＳ１１）。
【０１２２】
このようにして求められた修正ベクトル｛Δβ｝と前回の非線形項パラメータ｛β｝（最初は初期値）とを加算して今回の新しい非線形項パラメータ｛β＊｝を計算する（ステップＳ１２）。
【０１２３】
【数２６】
｛β＊｝＝｛β｝＋｛Δβ｝                          ・・・式（２６）
【０１２４】
このようにして求めた新しい非線形項パラメータ｛β＊｝を用いて上記のステップＳ４以降を繰り返し実行する。この繰り返しループのステップＳ６で計算した残差二乗和Ｓが前回の値より小さくなっていないときには、このＳが最小値に達したことになるので、｛β｝の計算ループから出る（ステップＳ７の“Ｎｏ”）。
【０１２５】
そして、再度縮小しないＦＲＦを入力して（ステップＳ１３）、このＦＲＦと既に決定された｛β｝を式（２４）に代入して線形項｛α｝を計算する（ステップＳ１４）。
【０１２６】
この過程は線形最小二乗法が適用できるため繰り返し計算は不要である。
【０１２７】
そして、最適なモード特性を示す｛α｝と｛β｝を出力する（ステップＳ１５）。
【０１２８】
多数の周波数応答関数を同時に考慮してモード特性を抽出する場合の計算機負荷は、同定に必要な計算時間と記憶容量である。実験現場において安価な電算機を使用して同定するためには、この両負荷を大きく軽減する必要がある。ここでは、実際の実験データを使用して、本発明での計算時間と記憶容量に関する優位性を検証する。
【０１２９】
図１でインパルスハンマー２で被試験物１を加振し、３軸加速度センサ４で応答を測定して得た実験データを使用してモード特性の同定を行う。
【０１３０】
モード特性の実験同定の入力データとして、まず１４０Ｈｚから７３６ＨｚまでのＦＲＦ行列を縮小して使用する。この周波数領域内には９個の固有モードがそれぞれ存在する。応答点は６３点×３方向である。
【０１３１】
図９（１），（２）は縮小前の点（１）のＸ方向応答のＦＲＦの位相及び振幅をそれぞれ示している。実線が実験ＦＲＦで細破線が同定したモード特性から構築したＦＲＦである。固有モード▲１▼が存在する１８０Ｈｚ付近には不規則雑音が混入している。
【０１３２】
図１０（１），（２）は全ＦＲＦを単純に重ね合わせて１つのＦＲＦに縮小した位相及び振幅の例を示している。ここでは不規則雑音は殆ど除去されているが、固有モード▲４▼と▲５▼に相当する共振峰が不明瞭である。
【０１３３】
図１１（１），（２）〜図１９（１），（２）は本発明のＦＲＦ縮小手段を用いて固有モード▲１▼〜▲９▼を強調するように縮小したＦＲＦを示している。ここでは図１１では固有モード▲１▼、図１２では固有モード▲２▼、図１３では固有モード▲３▼、図１４では固有モード▲４▼、図１５では固有モード▲５▼、図１６では固有モード▲６▼、図１７では固有モード▲７▼、図１８では固有モード▲８▼、図１９では固有モード▲９▼、というように狙った各モードが強調されると共に、不規則雑音のレベルが低下しており、本発明の縮小方法の有効性を示している。
【０１３４】
＜計算時間＞
各ＦＲＦ毎に、５７３個の周波数点のデータを用いて５個の固有モードのモード特性を同定する場合の計算時間の比較を、表１に示す。
【０１３５】
この表１より分かる通り、本発明における行列計算の効率化により、従来の方法の約１／３に計算時間が短縮できている。また入力データ（ＦＲＦ行列）を縮小したときには、ＦＲＦの数が固有モードの数（この場合には５）以上になると有効になることがわかる。
【０１３６】
【表１】
0003599882.tif 000019
【０１３７】
＜記憶容量の低減＞
従来例と本発明における同定計算時に使用する計算機の記憶容量を表２に示す。同定条件は、同時参照ＦＲＦ数３００、採用固有モード数５、各ＦＲＦ周波数点数５７３である。
【０１３８】
表２から、本発明に係る行列計算方法の改良により必要記憶容量が約１／７に減少している。入力データの縮小を併用すればさらに約１／８に減少し、従来例の１．７％になっていることが分かる。また，入力データ縮小の場合はＦＲＦの数が増加しても必要記憶容量は変化しない。
【０１３９】
即ち、外部記憶容量により制限されるまでのＦＲＦの量を、この６００Ｋｂｉｔｅで処理できる。
【０１４０】
【表２】
0003599882.tif 000020
【０１４１】
【発明の効果】
以上説明したように本発明に係る振動特性解析装置によれば、周波数応答関数の選択された固有モードに対応する虚数部の極性に＋１又は−１を乗じて正極性にした周波数応答関数を重ね合わせることにより固有モードに対応する該周波数応答関数を求め、該周波数応答関数の数を固有モードの数に縮小するように構成したので、データ縮小による精度の低下を招くことなく計算効率を著しく向上させることができる。
【０１４２】
また、線形近似による最小二乗法で修正ベクトルを求める正規方程式における正規行列を周波数応答関数の理論値と実験値との重み付き残差二乗和及びその２次微分項を省いた残差二乗和の線形項及び非線形項に関する微分係数で求めれば、より一層計算効率を上げることができる。
【０１４３】
【図面の簡単な説明】
【図１】本発明及び従来例に共通な振動特性解析装置を示したブロック図である。
【図２】本発明に係る振動特性解析装置で用いられるＦＲＦ（周波数応答関数）行列のデータ縮小を説明する図である。
【図３】振動測定データ（１点の１回の１方向成分のデータ）のＦＲＦ位相及び振幅を示すグラフ図である。
【図４】４つのＦＲＦの振幅を示すグラフ図である。
【図５】図４の各ＦＲＦ振幅を単純加算したときのグラフ図である。
【図６】本発明に係る振動特性解析装置により固有モード▲１▼に注目したＦＲＦデータ縮小を行ったときのＦＲＦ振幅の改良を示すためのグラフ図である。
【図７】本発明に係る振動特性解析装置により固有モード▲２▼に注目したＦＲＦデータ縮小を行ったときのＦＲＦ振幅の改良を示すためのグラフ図である。
【図８】本発明に係る振動特性解析装置における演算装置に格納され且つ実行される演算アルゴリズムを示したフローチャート図である。
【図９】未縮小振動測定データのＦＲＦ位相及び振動を示すグラフ図である。
【図１０】全振動測定データの各ＦＲＦデータを単純加算したときのグラフ図である。
【図１１】全振動測定データの各ＦＲＦデータを固有モード▲１▼に注目して縮小したときのグラフ図である。
【図１２】全振動測定データの各ＦＲＦデータを固有モード▲２▼に注目して縮小したときのグラフ図である。
【図１３】全振動測定データの各ＦＲＦデータを固有モード▲３▼に注目して縮小したときのグラフ図である。
【図１４】全振動測定データの各ＦＲＦデータを固有モード▲４▼に注目して縮小したときのグラフ図である。
【図１５】全振動測定データの各ＦＲＦデータを固有モード▲５▼に注目して縮小したときのグラフ図である。
【図１６】全振動測定データの各ＦＲＦデータを固有モード▲６▼に注目して縮小したときのグラフ図である。
【図１７】全振動測定データの各ＦＲＦデータを固有モード▲７▼に注目して縮小したときのグラフ図である。
【図１８】全振動測定データの各ＦＲＦデータを固有モード▲８▼に注目して縮小したときのグラフ図である。
【図１９】全振動測定データの各ＦＲＦデータを固有モード▲９▼に注目して縮小したときのグラフ図である。
【図２０】ＦＲＦ行列の例を示した図である。
【図２１】従来の特異分解法によるデータ縮小例を示した図である。
【符号の説明】
１  被試験物
２  インパルスハンマー
３  ロードセル
４  加速度センサ
５  演算装置
図中、同一符号は同一または相当部分を示す。
【図１】
0003599882.tif 000021
【図２】
0003599882.tif 000022
【図３】
0003599882.tif 000023
【図４】
0003599882.tif 000024
【図５】
0003599882.tif 000025
【図６】
0003599882.tif 000026
【図７】
0003599882.tif 000027
【図８】
0003599882.tif 000028
【図９】
0003599882.tif 000029
【図１０】
0003599882.tif 000030
【図１１】
0003599882.tif 000031
【図１２】
0003599882.tif 000032
【図１３】
0003599882.tif 000033
【図１４】
0003599882.tif 000034
【図１５】
0003599882.tif 000035
【図１６】
0003599882.tif 000036
【図１７】
0003599882.tif 000037
【図１８】
0003599882.tif 000038
【図１９】
0003599882.tif 000039
【図２０】
0003599882.tif 000040
【図２１】
0003599882.tif 000041