(57)【要約】      （修正有）【課題】構造物の振動測定において、測定周波数領域外に存在する固有振動モードの影響を考慮した実験モード解析を実現する振動測定装置を提供する。【解決手段】加振力を示す加振力信号および加振力により生じた測定対象物の測定点での加速度を示す加速度信号に基づいて、所定の周波数帯域における測定対象物の伝達関数を算出する高速フーリエ変換部と、伝達関数に基づいて、測定対象物のモード・パラメータを算出する低次モード演算部と、モード・パラメータに基づいて、周波数帯域外の固有振動数に対応する固有振動モードベクトルを推定する高次モードベクトル推定部と、固有振動モードベクトルに基づいて算出した第１の高次モード加速度応答、および、測定対象物の加速度応答と推定した周波数帯域外の固有振動モードベクトルとを乗算して得た第２の高次モード加速度応答の差を、周波数帯域において最小とする演算を行う最適化演算部とを具備する。【選択図】図１
2017044673.tif 000002
【特許請求の範囲】
【請求項１】
測定対象物の加振点への加振力を示す加振力信号および前記加振力により生じた前記測定対象物の測定点での加速度を示す加速度信号に基づいて、所定の周波数帯域における前記測定対象物の伝達関数を算出する高速フーリエ変換部と、
前記伝達関数に基づいて、前記周波数帯域における前記測定対象物のモード・パラメータを算出する低次モード演算部と、
前記周波数帯域における前記モード・パラメータに基づいて、前記周波数帯域外の固有振動数に対応する固有振動モードベクトルを推定する高次モードベクトル推定部と、
前記固有振動モードベクトルに基づいて前記周波数帯域外における前記測定対象物のモード・パラメータを算出し、該算出したモード・パラメータを用いて表した第１の高次モード加速度応答、および、前記測定対象物の伝達関数における加速度応答と前記高次モードベクトル推定部が推定した前記周波数帯域外の固有振動数に対応する固有振動モードベクトルとを乗算して得た第２の高次モード加速度応答の差を、前記周波数帯域において最小とする演算を行う最適化演算部と
を具備する振動測定装置。
【請求項２】
前記高次モードベクトル推定部は、前記周波数帯域における前記モード・パラメータのうち、前記周波数帯域の固有振動数に対応する固有振動モードベクトルとの直交化式を用いて前記周波数帯域外の固有振動数に対応する固有振動モードベクトルを推定することを特徴とする請求項１記載の振動測定装置。
【請求項３】
前記最適化演算部は、前記第１の高次モード加速度応答および前記第２の高次モード加速度応答それぞれの実数部の最小化と虚数部の最小化とを個別に最小化することを特徴とする請求項１または２に記載の振動測定装置。
【請求項４】
前記最適化演算部は、前記最小とする演算に先立ち、前記周波数帯域外のモード質量を算出することを特徴とする請求項１ないし３のいずれか１項に記載の振動測定装置。
【請求項５】
前記周波数帯域および前記周波数帯域外における前記測定対象物のモード・パラメータに基づいて、前記周波数帯域および前記周波数帯域外における前記測定対象物の特性行列を算出する行列演算部をさらに備えたことを特徴とする請求項１ないし４のいずれか1項に記載の振動測定装置。
【請求項６】
前記行列演算部は、前記測定対象物のモード・パラメータのうち、剛性行列で連成が生じない成分を０として演算することを特徴とする請求項５に記載の振動測定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は構造物の振動測定装置に関するものである。
【背景技術】
【０００２】
構造物の固有振動特性を実験的に把握しておくことは、回転機械など加振力が作用する機械において、共振や強制振動などの不適合要因を抽出し、機器の低振動化、低騒音化を図る上で重要である。固有振動測定において実験モード解析は、固有振動特性を周波数領域に加えて振動の形態（モード）を一体で評価できる非常に有効な実験解析手法として普及している。
【０００３】
実験モード解析は、基本的にＦＦＴアナライザによる周波数領域伝達関数に基づいて評価している。したがって、評価できる固有振動モードの周波数範囲に制約が存在する。例えば、アナライザの測定限界周波数（最高周波数）により、測定周波数範囲の上限が設定されるから、ＦＦＴアナライザを用いて振動測定を行う場合には、最高周波数および周波数分解能を指定する必要がある。すなわち、必要な分解能を保持し、かつ注目される周波数に適切な余裕を持たせた測定条件で測定を実施しなくてはならない。
【０００４】
しかし、構造物の固有振動数は一般的に多数（原理的には無限個）存在しているから、測定に際して設定した測定周波数範囲よりも周波数が高い高周波領域にも固有振動数が複数個存在することになる。したがって、ＦＦＴアナライザで測定された振動データには、最高周波数より高い固有振動数の応答成分も含まれている。
【０００５】
実験モード解析では、ＦＦＴアナライザの測定周波数範囲に存在する固有振動数に対応するモード・パラメータを算出する。モード・パラメータが得られれば、対象構造物の振動応答をモード応答の重ね合わせにより求めて対象構造物の振動特性を評価することができる。一方、再解析された構造物の振動応答は測定周波数範囲の固有振動数に対応するモード・パラメータから構成されており、測定周波数範囲外の固有振動数のモード応答は考慮されていない。このため、測定周波数範囲の上限あるは下限においては、実測の振動応答と再解析による振動応答の間に数値の乖離が生じてしまう。そこで、実験モード解析では測定範囲外のモード応答の影響を上限では剰余剛性、下限では剰余質量を用いて補正することが一般に行われている。
【０００６】
ＦＦＴアナライザは、通常０Ｈｚから測定することができるので、測定周波数範囲の下限が問題になることは稀である。一方上限については、測定範囲に影響を与える高周波固有振動数は、本来複数存在することが一般的であるにもかかわらず、剰余剛性は定数に過ぎない。実測値と再解析値との間の誤差をある程度縮小することはできるものの、完全に補正することは困難である。特に、測定周波数上限に近接して複数個の高周波固有振動数が存在する場合には誤差が拡大してしまう。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開平７−１６７７３６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
このように、従来の振動測定装置では、測定周波数の領域外に存在する固有振動モードの影響を加味した実験モード解析を行うことが困難であるという問題がある。実施形態は、かかる課題を解決するためになされたもので、構造物の振動測定において、測定周波数領域外、とりわけ高周波領域に存在する固有振動モードの影響を考慮した実験モード解析を実現する振動測定装置を提供することを目的とする。
【課題を解決するための手段】
【０００９】
実施形態の振動測定装置は、測定対象物の加振点への加振力を示す加振力信号および前記加振力により生じた前記測定対象物の測定点での加速度を示す加速度信号に基づいて、所定の周波数帯域における前記測定対象物の伝達関数を算出する高速フーリエ変換部と、前記伝達関数に基づいて、前記周波数帯域における前記測定対象物のモード・パラメータを算出する低次モード演算部と、前記周波数帯域における前記モード・パラメータに基づいて、前記周波数帯域外の固有振動数に対応する固有振動モードベクトルを推定する高次モードベクトル推定部と、前記固有振動モードベクトルに基づいて前記周波数帯域外における前記測定対象物のモード・パラメータを算出し、該算出したモード・パラメータを用いて表した第１の高次モード加速度応答、および、前記測定対象物の伝達関数における加速度応答と前記高次モードベクトル推定部が推定した前記周波数帯域外の固有振動数に対応する固有振動モードベクトルとを乗算して得た第２の高次モード加速度応答の差を、前記周波数帯域において最小とする演算を行う最適化演算部とを具備する。
【図面の簡単な説明】
【００１０】
【図１】第１実施形態の振動測定装置の構成を示す概略図である。
【図２】第１実施形態の振動測定装置の測定対象である構造物の構成を示す概略図である
【図３】第１実施形態の振動測定装置の動作を示すフローチャートである。
【図４Ａ】実測された図２の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図４Ｂ】実測された図２の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図５Ａ】第１実施形態の振動測定装置により求められた図２に示す構造物の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図５Ｂ】第１実施形態の振動測定装置により求められた図２に示す構造物の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図６Ａ】従来技術による振動測定装置により求められた図２に示す構造物の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図６Ｂ】従来技術による振動測定装置により求められた図２に示す構造物の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図７Ａ】従来技術による振動測定装置により求められた図２に示す構造物の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図７Ｂ】従来技術による振動測定装置により求められた図２に示す構造物の系の自己伝達関数（アクセレランス）を示す周波数応答関数図である。
【図８】第４実施形態の振動測定装置におけるモード・パラメータ推定方法を説明する自己伝達関数（コンプライアンス）を示す周波数応答関数図である。
【発明を実施するための形態】
【００１１】
（測定原理）
測定対象の構造物Oの運動方程式は、数式１により表すことができる。
【数１】
2017044673.tif 000003
ここで、［Ｍ］、［Ｃ］、［Ｋ］は、それぞれ質量行列、減衰行列、剛性行列であり、｛ｘ｝は変位ベクトル、｛ｆ｝は外力ベクトルである。行列の大きさ（自由度）は、測定点数（応答数）Ｎｏに等しい。質量行列、減衰行列および剛性行列を、総称して特性行列と呼ぶ。
【００１２】
数式１に対応して、固有振動モードベクトル｛φｒ｝；ｒ＝１，２，…，Ｎｉは、数式２の固有値問題を解くことで得ることができ、Ｎｏ行Ｎｉ列の大きさを持つ数式３に示す行列［Φ］を求める。
【数２】
2017044673.tif 000004
【数３】
2017044673.tif 000005
ここで、Ｎｉは入力数であり、解析では採用モード数、実験では測定モード数に相当する。
【００１３】
固有振動モードベクトル｛φｒ｝は正規直交基底をなすので、減衰行列が比例減衰行列である場合、数式４が成り立つ。
【数４】
2017044673.tif 000006
ここで、記号“「」”は対角行列を意味し、「ｍ」、「ｃ」および「ｋ」の対角成分はそれぞれモード質量、モード減衰、モード剛性と呼ばれ、数式５が成り立つ。
【数５】
2017044673.tif 000007
数式５では、ｍｒ＝１となるように固有振動モードベクトル｛φｒ｝を正規化している。この場合、モード剛性ｋｒは固有値ωｒ２と一致する。角固有振動数ωｒは固有値の平方根であり、固有振動数は数式６により求めることができる。
【数６】
2017044673.tif 000008
また、減衰率ζｒは数式７により求めることができる。
【数７】
2017044673.tif 000009
固有振動モードベクトル、モード質量、モード減衰、モード剛性およびこれらより誘導される定数を総称してモード・パラメータと呼ぶ。
【００１４】
実験モード解析では、数式１の運動方程式をＦＦＴアナライザの測定結果に合わせて周波数応答関数の形式で取り扱う。すなわち、１点加振の周波数応答関数｛Ｈ｝は、加振力｛ｆ｝の加振点成分をｆ＝ｅｓｔ（ｓ＝ｊω）、他の成分をゼロとして求めた応答｛ｘ｝のラプラス変換である。この場合、多点加振の周波数応答関数は、応答数・入力数（Ｎｏ×Ｎｉ）の行列となり、数式１は数式８と表せる。
【数８】
2017044673.tif 000010
ここで、行列［Ｉｄ］の各列｛Ｉｄ｝は、加振点が１、その他は０のベクトルである。実験モード解析は、数式８に基づいて周波数応答関数｛Ｈ｝にカーブ・フィッティングという技法を用いてモード・パラメータを算出する解析手法であり、原理上は入力数Ｎｉ対までのモード・パラメータを特定することが可能である。
【００１５】
いま、実験モード解析によってＮｉ対のモード・パラメータが求められているとする。このとき、固有振動モードベクトル｛φｒ｝；ｒ＝１，２，…，Ｎｉは正規直交基底をなすから、これら既知の固有振動モードベクトルに直交化式を適用することによって、測定周波数領域より高い周波数領域にある未知なる固有振動モードベクトル｛φｒ｝；ｒ＝（Ｎｉ＋１），（Ｎｉ＋２），…，Ｎｏを予測することができる。すなわち、直交化式として代表的なグラム・シュミットの手法を用いると、任意の初期ベクトル｛ｙｒ｝に対して数式９の演算を繰り返す。
【数９】
2017044673.tif 000011
数式９によって得られる予測固有振動モードベクトル｛φｒ｝；ｒ＝（Ｎｉ＋１），（Ｎｉ＋２），…，Ｎｏを含めて固有振動モードベクトル｛φｒ｝；ｒ＝１，２，…，Ｎｏは正規直交基底をなす。
【００１６】
一方、モード質量ｍｒ、モード減衰ｃｒ、モード剛性ｋｒについては、ｒ＝（Ｎｉ＋１），（Ｎｉ＋２），…，Ｎｏについて数式１０に示す加振力｛ｆ｝に対するモード応答に対して、数式１１に示すモード質量ｍｒ、モード減衰ｃｒ、モード剛性ｋｒによるモード座標系の１自由度運動方程式と、数式１２に示すその応答を考えることになる。
【数１０】
2017044673.tif 000012
【数１１】
2017044673.tif 000013
【数１２】
2017044673.tif 000014
すなわち、数式１０と数式１２の差が最小となるようにモード質量ｍｒ、モード減衰ｃｒ、モード剛性ｋｒを決定する。
【００１７】
さらに、入力数が応答数に等しい場合、すなわち、Ｎｉ＝Ｎｏのとき、モード行列［Φ］はＮｏ行Ｎｏ列の正方行列をなすから、数式４より各特性行列は数式１３により求めることができる。
【数１３】
2017044673.tif 000015
【００１８】
以上の手順によって、モード質量ｍｒ、モード減衰ｃｒ、モード剛性ｋｒ、固有振動モードベクトル｛φｒ｝；ｒ＝（Ｎｉ＋１），（Ｎｉ＋２），…，Ｎｏが決定され、測定周波数領域より高い周波数領域にある固有振動モードの影響を考慮することが可能となる。数式８の表記に従えば、｛ｘ｝は周波数応答関数｛Ｈ｝に、加振力｛ｆ｝は｛Ｉｄ｝にそれぞれ相当する。
【００１９】
（第１実施形態）
以下、図１を参照して、第１の実施形態に係る振動測定装置１の構成を説明する。図１に示すように、実施形態の振動測定装置は、インパルスハンマー１０、インパルスハンマー１０に内蔵されたロードセル１５、加速度ピックアップ２１…２４、信号ケーブル３０、ＦＦＴアナライザ４０、および演算部５０を備えている。
【００２０】
インパルスハンマー１０は、振動測定装置１が測定する振動を発生させるハンマーである。インパルスハンマー１０は、自身の加振力を測定するロードセル１５を備えており、ハンマーにより加えた振動を測定することができる。
【００２１】
加速度ピックアップ２１…２４は、測定対象たる構造物Ｏに配設される振動センサであり、構造物Ｏにおける所定の位置における振動応答を測定する。信号ケーブル３０は、ロードセル１５や加速度ピックアップ２１…２４からの測定信号を振動測定装置１へ導くケーブルである。
【００２２】
ＦＦＴアナライザ４０は、信号ケーブル３０から受けた測定信号に基づき高速フーリエ変換（ＦＦＴ）を行って伝達関数を算出する。演算部５０は、ＦＦＴアナライザ４０が算出した伝達関数に基づき、モード・パラメータや質量行列、剛性行列、減衰行列などを演算するコンピュータである。
【００２３】
図１に示すように、演算部５０は、低次モード演算部５１、高次ベクトル演算部５２、最適化演算部５３および行列演算部５４を備えている。低次モード演算部５１は、ＦＦＴアナライザ４０の測定周波数範囲内におけるモード・パラメータを演算するプロセッサである。高次ベクトル演算部５２は、測定周波数範囲外における高次固有振動モードベクトルを演算するプロセッサである。最適化演算部５３は、全モードのパラメータを演算するプロセッサである。行列演算部５４は、質量行列、剛性行列および減衰行列を演算するプロセッサである。
【００２４】
続いて、図１ないし３を参照して、実施形態の振動測定装置の動作を説明する。まず、測定対象たる構造物Ｏについて、測定点および加振点を設定する（ステップ１００。以下「Ｓ１００」のように称する）。以下の説明において、構造物Ｏは図２に示すようにＸ軸方向に運動する４つの自由度ｘ１、ｘ２、ｘ３、ｘ４からなる４自由度系と仮定する。構造物Ｏは、固定端６０に自由度ｘ１の質量６１を剛性６２と減衰６３を介して接続し、質量６１に対して自由度ｘ２の質量６４と自由度ｘ３の質量６７を剛性６５、減衰６６および剛性６８、減衰６９を介して直列に接続するものとする。また、自由度ｘ４の質量７０を質量６１に対して剛性７１、減衰７２を介して自由度ｘ２と並列に接続するものとする。各自由度の質量をｍ１、ｍ２、ｍ３、ｍ４、剛性をｋ１、ｋ２、ｋ３、ｋ４、減衰をｃ１、ｃ２、ｃ３、ｃ４とおく。
【００２５】
自由度ｘ１すなわち質量６１に対してインパルスハンマー１０により打撃加振する。打撃加振に応じて自由度ｘ１〜ｘ４に発生する加速度を加速度ピックアップ２１…２４で計測する。ロードセル１５が検出した加振力信号および加速度ピックアップ２１…２４が検出した各自由度の加速度信号は、信号ケーブル３０を介してＦＦＴアナライザ４０に送られる。ＦＦＴアナライザ４０は、受け取った加振力信号および加速度信号を用いて伝達関数を算出する（Ｓ１１０）。算出された伝達関数は、演算部５０に送られる。
【００２６】
ここで、構造物Ｏの諸元を表１に示す値とする。
【表１】
2017044673.tif 000016
ただし、減衰行列［Ｃ］については、比例減衰を数式１４に示すものとして、α＝５．０、β＝５．０×１０−５とする。
【数１４】
2017044673.tif 000017
このとき、数式１の運動方程式は、数式１５にて表される。
【数１５】
2017044673.tif 000018
【００２７】
数式１５に対応する数式２の固有値問題は、数式１６に示す通りとなる。
【数１６】
2017044673.tif 000019
低次モード演算部５１は、数式１６を演算して固有値、固有振動モードベクトルなどモード・パラメータを算出する。算出結果例を表２に示す。
【表２】
2017044673.tif 000020
【００２８】
低次モード演算部５１による解析では、表２に示すようにすべてのモードのモード・パラメータが求められる。しかし実際の測定では、ＦＦＴアナライザ４０の測定周波数領域の制約が存在する。そこで、ＦＦＴアナライザ４０の測定周波数範囲を０〜５０Ｈｚとすると、低次モード演算部５１は、表２に示す５０Ｈｚ以下の２つの固有振動数について、計測された伝達関数に基づいてモード・パラメータを決定する。
【００２９】
実測される伝達関数の一例として、図４Ａおよび図４Ｂに自己アクセレランス（加振点における応答を加速度で評価した伝達関数、すなわち自由度ｘ１における伝達関数を加速度表記したもの）を示す。図４Ａおよび図４Ｂは実測されたものであるから、測定周波数範囲外にある固有振動モードの応答も含まれている。
【００３０】
アクセレランスは動数領域において各自由度の応答加速度を加振力で除したものであるから、数式１５において数式１７に示す単位調波加振力Ｆを入力した場合の応答加速度に等しい。
【数１７】
2017044673.tif 000021
【００３１】
図５Ａおよび図５Ｂに示すように、１次および２次固有振動数でアクセレランス線図はピーク値を示す。測定周波数範囲が０〜５０Ｈｚの伝達関数を用いて実験モード解析を実施すれば、測定周波数範囲に固有振動数が存在する１次、２次モードについてモード・パラメータが数式１８、数式１９のように求められる。
【数１８】
2017044673.tif 000022
【数１９】
2017044673.tif 000023
【００３２】
表２に示すように、測定周波数範囲の上限である５０Ｈｚ以上には３次、４次の２つの固有振動数が存在する。しかし、測定周波数範囲外であるため測定周波数範囲を０〜５０Ｈｚとして測定された伝達関数の測定結果から３次、４次の固有振動数については求めることはできない。すなわち、３次、４次固有振動数に対応するモード・パラメータは未知数である。
【００３３】
数式１８、数式１９で得られたモード・パラメータを用いて自己アクセレランスを再解析すると、図６Ａおよび図６Ｂに示すようになる。図６Ａおよび図６Ｂにおいて破線で示す実測値と実線で示す演算結果とを比較すると、１次および２次固有振動数におけるピーク近傍で両者はよく一致しているが、それ以外の帯域では振幅に大きなずれが見られる。位相についても１次と２次固有振動数の中間で差が生じている。この相違は、測定振動数範囲以上にある３次および４次固有振動数に対応するモード応答が考慮されていないことによるものである。測定周波数範囲より高い固有振動数の影響を考慮して精度向上を図るため、一般に剰余剛性が用いられる。剰余剛性は数式２０にて求めることができる。
【数２０】
2017044673.tif 000024
【００３４】
剰余剛性を考慮した自己アクセレランスは図７Ａおよび図７Ｂに示すようになる。剰余剛性によって１次および２次固有振動数の中間部における相違は小さくなるものの、２次固有振動数を越えた周波数帯域では大きな差が残ることがわかる。これは、剰余剛性が数式２０に示すように３次および４次固有振動数に対応するモード剛性の逆数を総和して算出したものであり、各固有振動数を個々に評価した算出したものではないからである。
【００３５】
そこで、実施形態の振動測定装置は、３次および４次固有振動数に対応する固有振動モードベクトルを推定する高次ベクトル演算部５２を備えて、それぞれのモード・パラメータを算出する（Ｓ１３０）。高次ベクトル演算部５２は、３次および４次固有振動数に対応する固有振動モードベクトルを推定するに当たり、低次モード演算部５１が算出した１次および２次固有振動数の固有振動モードベクトルとの直交化式を利用する。ここでは直交化式として代表的なグラム・シュミットの式を用いる。
【００３６】
３次固有振動モードベクトル｛φ３｝は、数式９において初期ベクトルに｛１  １  ０  １｝Ｔを選べば、数式２１のように求められる。
【数２１】
2017044673.tif 000025
｛φ３｝Ｔ｛φ３｝＝１となるように正規化を施せば、数式２２のように表される。
【数２２】
2017044673.tif 000026
【００３７】
４次固有振動モードベクトル｛φ４｝は、初期ベクトルに｛１  １  １  １｝Ｔを選べば、同様に数式２３のように表される。
【数２３】
2017044673.tif 000027
このように、実施形態の振動測定装置によれば、グラム・シュミットの直交化式を用いて既知の１次、２次固有振動モードベクトルに対して直交する未知の３次、４次固有振動モードベクトルを推定するので、ＦＦＴアナライザ４０の測定周波数範囲外の固有振動数を考慮した精度の高いモード・パラメータや特性行列を得ることができる。１次、２次固有振動モードベクトルと併せてモード行列を作れば、数式２４が得られる。
【数２４】
2017044673.tif 000028
【００３８】
さらに、高次ベクトル演算部５２は、３次、４次のモード質量、モード減衰、モード剛性について、応答（伝達関数）ベクトルに固有振動モードベクトルを乗じて得られるモード応答と曲線適合させる。伝達関数は数式１７で表わされる加振力に対する応答として求められるから、数式２５に示すように、測定で得られた加速度応答に３次固有振動モードベクトル｛φ３｝Ｔを乗じて３次モード加速度応答を求める。
【数２５】
2017044673.tif 000029
一方、３次のモード・パラメータを用いて３次モード加速度応答を表すと、数式１５、数式１７より数式２６のようになる。
【数２６】
2017044673.tif 000030
【００３９】
最適化演算部５３は、数式２５で得られる３次モード加速度応答と、数式２６で求められる３次モード加速度応答との差（目的関数）が０〜５０Ｈｚの範囲で最小となるように最適化を行い、数式２６に現れるモード・パラメータm３、ω３、ζ３を決定する（Ｓ１４０）。
【００４０】
ω３、ζ３はモード質量m３、モード減衰ｃ３、モード剛性ｋ３から誘導できるから、数式２６は、モード質量ｍ３、モード減衰ｃ３、モード剛性ｋ３で定義される。すなわち、この問題の未知数は、モード質量ｍ３、モード減衰ｃ３、モード剛性ｋ３の３つであるから、適当な３点（たとえば１０Ｈｚ、３０Ｈｚ、５０Ｈｚ）における値で方程式を立てて算出する代数的な方法がまず考えられる。
【００４１】
しかし、実際の問題においては伝達関数に測定値を用いるので、図５Ａおよび図５Ｂに示すように理想的な特性を示すとは限らないため、このような代数的な算法では誤差が混入する恐れがある。そこでモード・パラメータｍ３、ｃ３、ｋ３を設計変数とし、０〜５０Ｈｚの周波数範囲で目的関数（３次モード加速度応答の差の絶対値）の総和を最小化する最適化問題に帰着させる。すなわち、数式２７にて表されるように０〜５０Ｈｚの範囲で目的関数の最小化を図る。
【数２７】
2017044673.tif 000031
【００４２】
同様に、４次モードについてもモード加速度応答について実測と解析の差を最小化することで４次のモードパラメータモードパラメータm４、ω４、ζ４を求める。これによって求められたモード・パラメータ例を表３に示す。
【表３】
2017044673.tif 000032
【００４３】
ここで１次、２次モードの値は、実験モード解析によって求められているから、表２と同一であるので表３では省略した。表３のモード・パラメータを用いて自己アクセレランスを計算した結果は図５Ａおよび図５Ｂに示す。図５Ａおよび図５Ｂにおいて、破線は実測の応答曲線であり、図４Ａおよび図４Ｂのものと同一である。実線で示す計算結果は破線の実測結果とほぼ一致しており、推定された３次、４次モード・パラメータの精度は非常に良好であることがわかる。実施形態の振動測定装置によれば、ＦＦＴアナライザ４０の測定周波数範囲外の固有振動数について、ベクトル直交化式による高次固有振動モードベクトルを算出し、モード応答の残差を最小化するので、ＦＦＴアナライザ４０の測定周波数範囲外においても制度の高い測定を可能にする。
【００４４】
続いて、行列演算部５４は、特性行列を算出する（Ｓ１５０）。まず、行列演算部５４は、表３の値を用いて、数式１３の行列を算出する。算出結果例を数式２８ないし３２に示す。
【数２８】
2017044673.tif 000033
【数２９】
2017044673.tif 000034
【数３０】
2017044673.tif 000035
【数３１】
2017044673.tif 000036
【数３２】
2017044673.tif 000037
【００４５】
続いて、行列演算部５４は、数式２８ないし３２を数式１３に代入することで、数式３３ないし３５に示すような特性行列を算出する。
【数３３】
2017044673.tif 000038
【数３４】
2017044673.tif 000039
【数３５】
2017044673.tif 000040
【００４６】
数式３３ないし３５に示す特性行列を真の値である数式１５の行列と比較すると、行列の対称性が保たれ、かつ主たる項（数値の大きい成分）の値は非常によく一致していることから実施形態の測定手順の妥当性が確認できる。
【００４７】
なお、伝達関数の再解析だけであれば、数式２５のように加速度応答と推定されたモードに対応する固有振動モードベクトル｛φｒ｝の内積を重ね合わせれば十分である。すなわち、数式３６に示すモード応答に対して、数式３７に示すアクセレランス・ベクトルにより求めることができる。
【数３６】
2017044673.tif 000041
【数３７】
2017044673.tif 000042
【００４８】
（第２実施形態）
第１実施形態では、最適化演算部５３は、数式２７のようにモード応答の実数部と虚数部の二乗平均を目的関数としている。この場合、最適化の過程で不適切な局所解に収束すると、実数部または虚数部のみが最小化されて真の最適解から外れてしまう可能性がある。そこで、第２の実施形態に係る振動測定装置では、最適化演算部５３が、実数部および虚数部を同時に最小化する。
【００４９】
すなわち、この実施形態の最適化演算部５３は、設計変数ｍ３、ｃ３、ｋ３に対して、数式３８にて表される制約条件のもとに、数式３９にて表される目的関数を最小化する演算を行う。４次モード・パラメータについても同様である。
【数３８】
2017044673.tif 000043
【数３９】
2017044673.tif 000044
これにより実数部の差と虚数部の差を同時に最小化することができるので、より正確な曲線適合が可能となる。
【００５０】
（第３実施形態）
剛性行列や減衰行列は、構造物Ｏの各自由度間の接合状態などによって影響を受けるから、解析的に特定することが一般に困難である。一方、質量行列は、構成から直接的に値を求めることが比較的容易である。すなわち、曲線適合する前に質量行列を決定しておくことが可能である。質量行列［Ｍ］が既知であれば、推定したモードベクトルとの三重積の対角項からモード質量ｍ３、ｍ４を算出することができる。ここで、１次、２次モードのモード・パラメータは既知であるから、この実施形態の最適化演算部５３は、数式４０により３次、４次モード質量のみ求めている。
【数４０】
2017044673.tif 000045
これによりｍ３、ｍ４が決定されるので、最適化演算部５３が演算する数式２７の最適化問題における未知数を３個から２個に減らすことができる。したがって、この実施形態の振動測定装置によれば、計算時間の削減や、計算精度の向上を図ることができる。
【００５１】
なお、第２の実施形態における数式３８および数式３９の演算について適用することもできる。すなわち、第２の実施形態において、最適化演算部５３が、３次、４次モード質量のみ先行して算出し、次いで数式３８および数式３９の演算を行ってもよい。
【００５２】
（第４実施形態）
数式２６のモード加速度に替えてモード変位を考えると、数式４１にように書き表すことができる。
【数４１】
2017044673.tif 000046
ここで、Ω＝０とおくと、数式４２が得られる。
【数４２】
2017044673.tif 000047
これは、伝達関数としてコンプライアンスの周波数軸切片の値からモード剛性ｋ３が計算できることを意味している。
【００５３】
図８は、３次および４次モード応答を、コンプライアンスで表記したものである。図中のＰ３およびＰ４の値からモード剛性ｋ３、ｋ４を求めると、数式４３および数式４４のようになり、これは表２に示す値をほぼ等しくなることがわかる。
【数４３】
2017044673.tif 000048
【数４４】
2017044673.tif 000049
ただし、コンプライアンスを測定すると０Ｈｚ近傍でノイズの影響を受けやすいので、実際には低周波領域でノイズの少ない平滑な領域を補外して０Ｈｚの値を推定することになる。
【００５４】
これによりｋ３、ｋ４が決定されるので、最適化演算部５３が演算する数式２７の最適化問題における未知数を３個から２個に減らすことができる。したがって、計算時間の削減や計算精度の向上を図ることができる。
【００５５】
なお、第２の実施形態における数式３８および数式３９の演算について適用することもできる。すなわち、第２の実施形態において、最適化演算部５３が、３次、４次モード剛性のみ先行して算出し、次いで数式３８および数式３９の演算を行ってもよい。
【００５６】
（第５実施形態）
第１実施形態では特性行列として数式３３ないし３５に示すような値が得られた。ここで、数式３４の減衰行列と数式３５の剛性行列を詳細に見ると、真の値である数式１５の行列では０となるべき成分（１行３列、３行１列、２行４列、４行２列、３行４列、４行３列）の値について、微小ではあるものの誤差が生じている。直接接続されていない自由度の間（自由度ｘ１と自由度ｘ３、自由度ｘ２と自由度ｘ４、自由度ｘ３と自由度ｘ４）では相互に連成が生じないから、連成が生じない上記の成分は、構造物Ｏの接続構成から特定が可能である。そこで、数式３５の剛性行列で連成が生じない成分を０に設定する。すなわち、この実施形態では行列演算部５４は、数式４５を演算する。
【数４５】
2017044673.tif 000050
【００５７】
数式３３と数式４５を用いて数式２の固有値問題を解くと、モード行列として数式４６が得られる。
【数４６】
2017044673.tif 000051
【００５８】
グラム・シュミットの方法では、数式９に見るように初期ベクトルが必要であり、初期ベクトルが真のベクトルに近似しているほど精度が高まる。そこで、数式４５で得られた固有振動モードベクトルを初期値として直交化式である数式９に代入してモード行列の改良をはかると、減衰行列、剛性行列はそれぞれ数式４７および数式４８のように表される。
【数４７】
2017044673.tif 000052
【数４８】
2017044673.tif 000053
すなわち、非連成項の値が０に近付き精度の改善が見られた。必要に応じて行列演算部５４がこの操作を繰り返せば減衰行列［Ｃ］、剛性行列［Ｋ］を真の値に漸近させることができる。これにより計算精度のさらなる向上がはかられる。
【００５９】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。
【符号の説明】
【００６０】
１…振動測定装置、１０…インパルスハンマー、１５…ロードセル、２１〜２４…加速度ピックアップ、３０…信号ケーブル、４０…ＦＦＴアナライザ、５０…演算部、５１…低次モード演算部、５２…高次ベクトル演算部、５３…最適化演算部、５４…行列演算部。
【図１】
2017044673.tif 000054
【図２】
2017044673.tif 000055
【図３】
2017044673.tif 000056
【図４Ａ】
2017044673.tif 000057
【図４Ｂ】
2017044673.tif 000058
【図５Ａ】
2017044673.tif 000059
【図５Ｂ】
2017044673.tif 000060
【図６Ａ】
2017044673.tif 000061
【図６Ｂ】
2017044673.tif 000062
【図７Ａ】
2017044673.tif 000063
【図７Ｂ】
2017044673.tif 000064
【図８】
2017044673.tif 000065