(57)【要約】【課題】物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の実数部及び虚数部の一方のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換する。【解決手段】地盤の動的剛性をインパルス応答へ変換する際に、動的剛性データに実数部のデータが含まれていない等の場合は、虚数部のデータのみを用いてインパルス応答を演算し(110,112)、当該演算で算出されない同時成分ｋ０,ｍ０を別途の演算等によって決定する(114〜124)。また、動的剛性データに虚数部のデータが含まれていない等の場合は、実数部のデータのみを用いてインパルス応答を演算し(126,128)、当該演算で算出されない同時成分ｃ０を別途の演算等によって決定する(130〜140)。そして、算出されたインパルス応答を用いて建物の時刻歴応答解析を行う(142)。【選択図】図８
2008151643.tif 000002
【特許請求の範囲】
【請求項１】
物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するにあたり、
物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔＪ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として
【数１】
2008151643.tif 000014
上記(１)式を用いると共に、前記物体の動的剛性の虚数成分ＳＩ(ω)を規定する数式として
【数２】
2008151643.tif 000015
上記(２)式を用い（但し、Ｎ１＝Ｎ／２,Ｎ２＝Ｎ／２−１、又は、Ｎ１＝Ｎ２＝(Ｎ−１)／２）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を
【数３】
2008151643.tif 000016
上記(３)式及び(４)式へ代入して演算することで、前記インパルス応答を求めることを特徴とするインパルス応答演算方法。
【請求項２】
前記物体の動的剛性Ｓ(ω)を規定する
【数４】
2008151643.tif 000017
上記(５)式に、請求項１記載のインパルス応答演算方法によって求めた前記インパルス応答の各成分を代入すると共に、(５)式における同時成分ｋ０,ｍ０に０を設定し、前記振動がＮ３種の周波数のときの物体の動的剛性の実数成分ＳＲ(ωｉ)を各々演算し、
演算した実数成分ＳＲ(ωｉ)と、当該実数成分ＳＲ(ωｉ)と同一周波数での動的剛性の実数成分を表すＮ３個の実数部データＤＲ(ωｉ)を
【数５】
2008151643.tif 000018
上記(６),(７)式に代入して同時成分ｋ０,ｍ０を演算することを特徴とする請求項１記載のインパルス応答演算方法。
【請求項３】
物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するにあたり、
物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として
【数６】
2008151643.tif 000019
上記(１)式を用いると共に、前記物体の動的剛性の実数成分ＳＲ(ω)を規定する数式として
【数７】
2008151643.tif 000020
上記(８)式を用い（但し、Ｎ１＝Ｎ２＝(Ｎ／２)−１、又は、Ｎ１＝(Ｎ−１)／２,Ｎ２＝((Ｎ−１)／２)−１）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を
【数８】
2008151643.tif 000021
上記(９)式及び(10)式へ代入して演算することで、前記インパルス応答を求めることを特徴とするインパルス応答演算方法。
【請求項４】
前記物体の動的剛性Ｓ(ω)を規定する
【数９】
2008151643.tif 000022
上記(５)式に、請求項３記載のインパルス応答演算方法によって求めた前記インパルス応答の各成分を代入すると共に、(５)式における同時成分ｃ０に０を設定し、前記振動がＮ４種の周波数のときの物体の動的剛性の虚数成分ＳＩ(ωｉ)を各々演算し、
演算した虚数成分ＳＩ(ωｉ)と、当該虚数成分ＳＩ(ωｉ)と同一周波数での動的剛性の虚数成分を表すＮ４個の虚数部データＤＩ(ωｉ)を
【数１０】
2008151643.tif 000023
上記(11),(12)式に代入して同時成分ｃ０を演算することを特徴とする請求項３記載のインパルス応答演算方法。
【請求項５】
物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置であって、
物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として
【数１１】
2008151643.tif 000024
上記(１)式を用いると共に、前記物体の動的剛性の虚数成分ＳＩ(ω)を規定する数式として
【数１２】
2008151643.tif 000025
上記(２)式を用い（但し、Ｎ１＝Ｎ／２,Ｎ２＝Ｎ／２−１、又は、Ｎ１＝Ｎ２＝(Ｎ−１)／２）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を
【数１３】
2008151643.tif 000026
上記(３)式及び(４)式へ代入して演算することで、前記インパルス応答を求める第１演算手段を備えたことを特徴とするインパルス応答演算装置。
【請求項６】
物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置であって、
物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として
【数１４】
2008151643.tif 000027
上記(１)式を用いると共に、前記物体の動的剛性の実数成分ＳＲ(ω)を規定する数式として
【数１５】
2008151643.tif 000028
上記(８)式を用い（但し、Ｎ１＝Ｎ２＝(Ｎ／２)−１、又は、Ｎ１＝(Ｎ−１)／２,Ｎ２＝((Ｎ−１)／２)−１）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を
【数１６】
2008151643.tif 000029
上記(９)式及び(10)式へ代入して演算することで、前記インパルス応答を求める第２演算手段を備えたことを特徴とするインパルス応答演算装置。
【請求項７】
コンピュータを、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置として機能させるインパルス応答演算プログラムであって、
前記コンピュータを、
物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として
【数１７】
2008151643.tif 000030
上記(１)式を用いると共に、前記物体の動的剛性の虚数成分ＳＩ(ω)を規定する数式として
【数１８】
2008151643.tif 000031
上記(２)式を用い（但し、Ｎ１＝Ｎ／２,Ｎ２＝Ｎ／２−１、又は、Ｎ１＝Ｎ２＝(Ｎ−１)／２）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を
【数１９】
2008151643.tif 000032
上記(３)式及び(４)式へ代入して演算することで、前記インパルス応答を求める第１演算手段として機能させることを特徴とするインパルス応答演算プログラム。
【請求項８】
コンピュータを、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置として機能させるインパルス応答演算プログラムであって、
前記コンピュータを、
物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として
【数２０】
2008151643.tif 000033
上記(１)式を用いると共に、前記物体の動的剛性の実数成分ＳＲ(ω)を規定する数式として
【数２１】
2008151643.tif 000034
上記(８)式を用い（但し、Ｎ１＝Ｎ２＝(Ｎ／２)−１、又は、Ｎ１＝(Ｎ−１)／２,Ｎ２＝((Ｎ−１)／２)−１）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を
【数２２】
2008151643.tif 000035
上記(９)式及び(10)式へ代入して演算することで、前記インパルス応答を求める第２演算手段として機能させることを特徴とするインパルス応答演算プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明はインパルス応答演算方法、装置及びプログラムに係り、特に、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算方法、該インパルス応答演算方法を適用可能なインパルス応答演算装置、コンピュータをインパルス応答演算装置として機能させるためのインパルス応答演算プログラムに関する。
【背景技術】
【０００２】
地震時に震源から地盤を伝播した地震動（地震波）は、表層地盤を介して建物の基礎に入射し、その一部は基礎によって反射され残りは建物上部へ伝達されて建物の振動を引き起こすが、地震動によって振動された建物は新たな震源となって地盤へ波動を放出する。このように、地震時には建物と地盤が互いに影響し合って振動するので、地震時の建物の挙動や耐震安全性等を地震応答解析によって解析・評価する際には、地震動に対する地盤の挙動（反力の特性）も考慮する必要がある。地震動と地盤の挙動との関係は地盤の動的剛性（地盤インピーダンスともいう）で表すことができる。地震応答解析は、周波数領域で応答解析を行う周波数応答解析と、時間領域で応答解析を行う時刻歴応答解析とに大別されるが、地盤の動的剛性は、振動の周波数の変化に応じて実数部及び虚数部の値が変化する周波数領域の複素関数であるため、地震応答解析では地盤の動的剛性をそのまま利用可能な周波数応答解析も多用されている。
【０００３】
但し、大きなエネルギーが構造物に入力される大地震時等には、そのエネルギーの一部が、構造物を構成する各部材の内部に亀裂を生じさせたり各部材を部分的に塑性化させる等によって消費されると共に、この亀裂発生や部分的な塑性化等に伴い各部材の強度が低下することが繰り返されるというプロセスを経るため、構造物の挙動は非線形性を有している。このため、地震時の構造物の挙動等を高精度に予測解析するためには、時刻歴応答解析により、地震時の各時点での構造物の状態（過去にどのような力が加わり、その力によってどのような状態になっているか）を考慮して各時点での構造物の挙動を解析する必要がある。時刻歴応答解析を行うにあたっては、周波数領域の複素関数である地盤の動的剛性を、時間領域で表されるインパルス応答へ変換して用いる必要がある。
【０００４】
このため、本願発明者は、物体を振動させる外力と物体の挙動との関係を時間領域で表すインパルス応答を規定する数式として、物体の変位に依存し同時成分と時間遅れ成分から成る剛性項と、物体の速度に依存し同時成分と時間遅れ成分から成る減衰項と、物体の加速度に依存し少なくとも同時成分を含んで成る質量項を含む数式を用い、前記振動が各周波数のときの前記関係を周波数領域で表す動的剛性の値に基づいてインパルス応答を求める演算方法を提案している（特許文献１も参照）。
【０００５】
また、特許文献１に記載の演算方法は、変換対象の動的剛性に比較的大きな履歴減衰成分が含まれている場合に、算出されるインパルス応答の精度が低下するという欠点を有しており、本願発明者は、これを解決するために、特許文献１に記載のインパルス応答演算方法を適用してインパルス応答を求めた後に、求めた時刻ｔ＝ｔ0〜ｔnの期間のインパルス応答を表すデータのうち、時刻ｔ０に相当する同時成分のデータと、時刻ｔ１〜ｔｎ'の期間に相当するｎ'個の時間遅れ成分（但しｎ'＜ｎ、ｔｎ'＜Δｔ・ｎ'）のデータを用いて物体の動的剛性を再現する演算を行い、再現した動的剛性と元の動的剛性の偏差に基づいて、インパルス応答を表すデータのうち剛性項の同時成分のデータ及び質量項の同時成分のデータを修正する演算方法も提案している（特許文献２も参照）。
【０００６】
一方、非特許文献１には、地盤の動的剛性等の複素関数のうちの虚数部のデータのみから逆フーリエ変換によって近似的にインパルス応答を算出する方法が提案されている。この方法では、得られたインパルス応答から元の複素関数を再現すると、元の虚数部が再現されるのと同時に、変換に用いられなかった未知の実数部も同時に再現されるので、虚数部のデータから当該虚数部に対応しかつ時間領域で因果律を満たす実数部のデータを得ることができる。また、非特許文献１に記載の方法によれば、同様に、実数部のみのデータから、当該実数部に対応しかつ時間領域で因果律を満たす虚数部のデータを得ることも可能である。
【特許文献１】特許第３８７８６２６号公報
【特許文献２】特開２００６−０７８２９０号公報
【非特許文献１】林  康裕，勝倉  裕，「周波数領域の地盤剛性や地盤柔性を用いた地盤・建屋相互作用系の時刻歴応答解析」，日本建築学会構造系論文報告集，１９９０年７月，第４１３号，ｐ．６５−７４
【発明の開示】
【発明が解決しようとする課題】
【０００７】
図６(Ａ)にも示すように、特許文献１及び特許文献２に記載の演算方法では、変換対象の動的剛性の実数部のデータ及び虚数部のデータが既知であれば、これらのデータからインパルス応答を直接求めることができる。しかし、特許文献１及び特許文献２に記載の演算方法では、例えば物体の動的剛性が所望の特性を示すように新たなモデルを創出してインパルス応答を求める等の場合にも、動的剛性の実数部及び虚数部のデータを各々必要とするので、新たなモデルを創出する等の作業が煩雑であるという問題がある。
【０００８】
すなわち、物体の動的剛性をインパルス応答へ変換して時刻歴応答解析に用いるためには、動的剛性が時間領域で因果律を満たしている必要がある。新たなモデルの創出にあたって上記の条件を満足するように動的剛性を定めることは、動的剛性の実数部及び虚数部のうちの一方を、所望の特性を示すように定めた後に、動的剛性の実数部と虚数部が周波数領域でヒルベルト(Hilbert)変換対を形成するように、動的剛性の実数部及び虚数部のうちの既知の一方から未知の他方を演算によって求めることで実現できる。
【０００９】
ヒルベルト変換の演算式を以下に示す。なお、以下の演算式において、ＺＲ(ω),ＺＩ(ω)は時間領域で因果律を満たす複素関数の実数部及び虚数部、ＰはCauchyの主値を表す。
【００１０】
【数１】
2008151643.tif 000003
【００１１】
上記の演算式からも明らかなように、ヒルトベルト変換は積分区間が−∞から∞に亘っているので、この演算式そのままでは数値演算は実質的に不可能である。このため、ヒルベルト変換の数値演算は、積分区間を変化させながら演算を繰り返して演算結果を比較し、積分区間の変化に対する演算結果の変化が閾値以下となる積分区間を判定し、当該積分区間での演算結果を最終的な演算結果とする処理を、各振動数について各々行うことによって実現される。このため、新たなモデルの創出にあたって時間領域で因果律を満たすように動的剛性を定めるための演算が非常に煩雑であるという問題がある。また、上記のヒルベルト変換は複素関数が正則成分のみから成る状態（特異成分を含まない状態）で成立する（ヒルベルト変換対を形成する）が、物体の動的剛性等の多くの複素関数には正則成分と特異成分が混在しているので、多くの場合、ヒルベルト変換の前に複素関数から特異成分を分離・除去する必要があり、この特異成分の分離・除去も非常に煩雑である。
【００１２】
一方、非特許文献１に記載の方法は、図６(Ｂ)にも示すように、実数部のデータ及び虚数部のデータの一方が既知であれば適用可能であり、煩雑なヒルベルト変換も不要になるという利点を有している。しかし、非特許文献１に記載の方法における逆フーリエ変換についても、正則成分と特異成分が混在している複素関数を対象とする場合には、ヒルベルト変換と同様に特異成分を分離・除去する必要がある。また、非特許文献１に「ω→０における漸近項」や「ω→∞において発散する項」と記載されているように、非特許文献１に記載の方法では変換方法等の条件によって特異成分の現れ方が異なり、特異成分の現れ方に応じて特異成分の適切な分離・除去方法が相違すると共に、特異成分の分離・除去の良否が変換精度に影響する。そして、特異成分の現れ方に応じた特異成分の適切な分離・除去方法の選択には経験を要するので、実際には、個々の分離・除去方法を順に試しその条件での最適な分離・除去方法を選択等の試行錯誤的な演算が必要となり、演算が非常に煩雑であるという問題があった。
【００１３】
本発明は上記事実を考慮して成されたもので、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の実数部及び虚数部の一方のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換できるインパルス応答演算方法、インパルス応答演算装置及びインパルス応答演算プログラムを得ることが目的である。
【課題を解決するための手段】
【００１４】
前述のように、時間領域で因果律を満たす複素関数の実数部及び虚数部は周波数領域でヒルベルト変換対を形成し、複素関数の実数部及び虚数部のうちの虚数部のみが既知の場合にも、ヒルベルト変換により、当該虚数部から当該虚数部に対応しかつ時間領域で因果律を満たす実数部を一意に算出できるので、虚数部には実数部の成分も含まれていると見なすことができる。同様に、複素関数の実数部及び虚数部のうちの実数部のみが既知の場合にも、ヒルベルト変換により、当該実数部から当該実数部に対応しかつ時間領域で因果律を満たす虚数部を一意に算出できるので、実数部にも虚数部の成分が含まれていると見なすことができる。また、時間領域で因果律を満たさない（ヒルベルト変換対を形成しない）複素関数についても、実数部及び虚数部のうちの一方が既知であれば、未知の他方も近似によっておおよそ求まるので、時間領域で因果律を満たす複素関数と同様に、虚数部には実数部の成分が、実数部には虚数部の成分が含まれていると見なすことができる。
【００１５】
本願発明者は上記事実に事実に着目し、特許文献１に記載のインパルス応答演算方法を基礎として、物体の動的剛性のうちの虚数成分及び実数成分の一方のみを用いて時間領域変換（インパルス応答の演算）を行うことに想到し、当該演算の演算式を導出した。より詳しくは、特許文献１に記載のインパルス応答演算方法では、物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性Ｓ(ω)を次の(１)式で、前記関係を時間領域で表すインパルス応答Ｆ(t)を次の(２)式で各々規定しており、
【００１６】
【数２】
2008151643.tif 000004
【００１７】
このうちの(１)式の動的剛性Ｓ(ω)を実数成分ＳＲ(ω)と虚数成分ＳＩ(ω)に分けて表すと次の(３)式のようになる。
【００１８】
【数３】
2008151643.tif 000005
【００１９】
なお、上記の(４),(５)式は、(１)式で表される関係を動的剛性のＮ個のデータ(Ｄ(ω１)〜Ｄ(ωＮ))に対してマトリクス表示したものであり、この２Ｎ×２Ｎの係数マトリクスを有する連立方程式を解くことで、未知数であるインパルス応答の値（ｋ０〜ｋＮ−１,ｃ０〜ｃＮ−２,ｍ０）を求めることができる。
【００２０】
本願発明者は、上記の(３)式を基礎として、物体の動的剛性の虚数成分ＳＩ(ω)のみを用いて時間領域変換（インパルス応答の演算）を行う数式として以下の(６)式を定め、
【００２１】
【数４】
2008151643.tif 000006
【００２２】
(６)式で表される関係を、動的剛性のＮ個の虚数部のデータ(ＤＩ(ω１)〜ＤＩ(ωＮ))のデータに対してマトリクス表示した、以下の(７),(８)式を導出した。
【００２３】
【数５】
2008151643.tif 000007
【００２４】
また本願発明者は、先の (３)式を基礎として、物体の動的剛性の実数成分ＳＲ(ω)のみを用いて時間領域変換（インパルス応答の演算）を行う数式として以下の(９)式を定め、
【００２５】
【数６】
2008151643.tif 000008
【００２６】
(９)式で表される関係を、動的剛性のＮ個の実数部のデータ(ＤＲ(ω１)〜ＤＲ(ωＮ))のデータに対してマトリクス表示した、以下の(10),(11)式を導出した。
【００２７】
【数７】
2008151643.tif 000009
【００２８】
そして本願発明者は、上記の(６)〜(８)式の動的剛性の虚数成分のみを用いたインパルス応答の演算（以下、虚部変換法と称する）及び上記の(９)〜(11)式の動的剛性の実数成分のみを用いたインパルス応答の演算（以下、実部変換法と称する）の有効性を確認するために、以下のような解析検討を行った。
【００２９】
この解析検討では、解析モデルとして、図１(Ａ)に示すように、地震波の伝播速度Vs2＝800(ｍ/ｓ)、ポアソン比ν＝0.45、密度ρ＝2.0(ｔ/ｍ３)、減衰率ｈ２＝５％の地層と、地震波の伝播速度Vs1＝400(ｍ/ｓ)、ポアソン比ν＝0.45、密度ρ＝2.0(ｔ/ｍ３)、減衰率ｈ１＝０％又は20％、層厚Ｈ＝20(ｍ)の地層から成る地盤の上に50×50(ｍ)の正方形の基礎が設けられているモデルを用いた。この解析モデルにおいて、減衰率ｈ１＝０％及びｈ１＝20％の各条件毎に薄層要素法によって求めた地盤の動的剛性の一例を図１(Ｂ)に示す。なお図１(Ｂ)において、"Real"は動的剛性の実部を、"Imag"は動的剛性の虚部を各々意味している。
【００３０】
本願発明者は、図１(Ｂ)に示す減衰率ｈ1＝０％での地盤の動的剛性を用い、０〜20(Hz)の周波数範囲内における地盤の動的剛性の複素データを所定個抽出し、抽出した複素データを特許文献１に記載のインパルス応答演算方法（旧変換法）によりΔｔ＝0.05秒の条件で時間領域へ変換することでインパルス応答を求めると共に、抽出した複素データを前述の虚部変換法によりΔｔ＝0.05秒の条件で時間領域へ変換することでインパルス応答を求め、更に、抽出した複素データを前述の実部変換法によりΔｔ＝0.05秒の条件で時間領域へ変換することでインパルス応答を求めた。なお、旧変換法では複素データの実数部データ及び虚数部データが演算に各々用いられるのに対し、虚部変換法及び実部変換法では実数部データ及び虚数部データの一方のみが演算に用いられるので、精度を揃えるために、虚部変換法及び実部変換法で演算に用いたデータ点の数を旧変換法で演算に用いたデータ点の数の約２倍とした（虚部変換法及び実部変換法でのデータ点数＝41、旧変換法でのデータ点数＝21）。
【００３１】
上記演算の結果を、剛性項(物体の変位に依存するインパルス応答)、減衰項(物体の速度に依存するインパルス応答)及び質量項(物体の加速度に依存するインパルス応答)に分けて図２(Ａ)〜(Ｃ)に示す。なお、図２では旧変換法は"Both"、虚部変換法は"Imag"、実部変換法は"Real"と表記している。また、虚部変換法では剛性項の同時成分ｋ０と質量項の同時成分ｍ０が算出されず((６)〜(８)式参照)不定となり、実部変換法では減衰項の同時成分ｃ０が算出されず((９)〜(11)式参照)不定となるため、図２(Ａ)〜(Ｃ)ではこれらを"0"として示している。図２(Ａ)〜(Ｃ)を参照しても明らかなように、減衰率ｈ1＝０％の条件では、旧変換法、虚部変換法及び実部変換法によって算出されるインパルス応答はおおよそ一致している。
【００３２】
また本願発明者は、減衰率ｈ1＝０％の条件で各変換法によって得られたインパルス応答を周波数領域へ各々変換して動的剛性を再現する演算を行った。この演算によって再現された動的剛性を各変換法毎に図３(Ａ)に示す。なお、上記の演算にあたり、虚部変換法で算出されない剛性項の同時成分ｋ０と質量項の同時成分ｍ０、実部変換法で算出されない減衰項の同時成分ｃ０については、旧変換法で算出された値を用いた。図３(Ａ)からも明らかなように、各変換法によって算出されたインパルス応答から再現される動的剛性は元の複素データ（データ点）におよそ一致している。従って、減衰率ｈ1＝０％の条件では、虚部変換法を適用し虚数部データのみを用いてインパルス応答を算出した場合にも良好な精度で実数部データが得られると共に、実部変換法を適用し実数部データのみを用いてインパルス応答を算出した場合にも良好な精度で虚数部データが得られることが確認された。
【００３３】
また本願発明者は、図１(Ｂ)に示す減衰率ｈ1＝20％での地盤の動的剛性についても、０〜20(Hz)の周波数範囲内における地盤の動的剛性の複素データを所定個抽出し、抽出した複素データを旧変換法、虚部変換法及び実部変換法によりΔｔ＝0.05秒の条件で時間領域へ変換することでインパルス応答を各々求めた。この減衰率ｈ1＝20％の条件においても、精度を揃えるために、虚部変換法及び実部変換法で演算に用いたデータ点の数を旧変換法で演算に用いたデータ点の数の約２倍とした（虚部変換法及び実部変換法でのデータ点数＝41、旧変換法でのデータ点数＝21）。
【００３４】
上記演算の結果を、剛性項、減衰項及び質量項(物体の加速度に依存するインパルス応答)に分けて図２(Ｄ)〜(Ｆ)に示す。なお、図２(Ｄ)〜(Ｆ)においても、虚部変換法で算出されない剛性項の同時成分ｋ０と質量項の同時成分ｍ０、実部変換法で算出されない減衰項の同時成分ｃ０を各々"0"として示している。図２(Ｄ)〜(Ｆ)を参照しても明らかなように、減衰率ｈ1＝20％の条件で虚部変換法及び実部変換法によって得られるインパルス応答は、旧変換法によって得られるインパルス応答に対し、剛性項の0.05〜0.10秒の範囲の時間遅れ成分及び質量項の同時成分の偏差が比較的大きくなっているものの、旧変換法、虚部変換法及び実部変換法によって算出されるインパルス応答は全体としてはおおよそ近似している。
【００３５】
更に本願発明者は、減衰率ｈ1＝20％の条件で各変換法によって得られたインパルス応答についても周波数領域へ各々変換して動的剛性を再現する演算を行った。この演算によって再現された動的剛性を各変換法毎に図３(Ｂ)に示す。なお、この演算においても、虚部変換法や実部変換法で算出されない同時成分については旧変換法で算出された値を用いた。図３(Ｂ)からも明らかなように、各変換法によって算出されたインパルス応答から再現される動的剛性は全体としては元の複素データ（データ点）におよそ近似しているが、低振動数域では偏差が若干大きくなっている。これは、減衰率ｈ1＝20％の条件で虚部変換法及び実部変換法によって得られるインパルス応答が、旧変換法によって得られるインパルス応答に対し、剛性項の0.05〜0.10秒の範囲の時間遅れ成分の偏差が比較的大きいことが影響しているものと推察される。
【００３６】
図３(Ｂ)のうちの０〜５(Hz)の振動数域を図３(Ｃ)に拡大して示す。図３(Ｃ)より明らかなように、虚部変換法及び実部変換法では、虚数部及び実数部のうちインパルス応答の演算にデータを用いている側（虚部変換法であれば虚数部、実部変換法であれば実数部）については元の複素データとの偏差が比較的小さい一方、その逆の側については元の複素データとの偏差が比較的大きいことがわかる。旧変換法では実数部、虚数部共にそれらの中間程度の精度となっている。
【００３７】
物体の動的剛性が所望の特性を示す新たなモデルを創出する等の場合に従来用いられていたヒルベルト変換は、時間領域で因果律を満たす複素関数（因果的な動的剛性）についてのみ成り立つ関係であり、時間領域で因果律を満たさない複素関数（非因果的な動的剛性、例えば減衰率ｈ＞0％での動的剛性）については本来は成り立たず、動的剛性の実数部及び虚数部のうちの既知の一方からヒルベルト変換に基づいて未知の他方を演算によって求めるにあたり、特異成分を分離・除去するという煩雑な処理を必要とし、またヒルベルト変換に基づく演算自体も煩雑である。また、非特許文献１の技術についても特異成分を分離・除去するという煩雑な処理が必要となる。
【００３８】
これに対し、上述した虚部変換法や実部変換法は、図２(Ｄ)〜(Ｆ)及び図３(Ｂ),(Ｃ)に示す結果からも明らかなように、物体の動的剛性が非因果的であったとしても、動的剛性の虚数部及び実数部のうちの既知の一方のみから、実用上十分な精度でインパルス応答を算出することができる。また、虚部変換法の演算式((６)〜(８)式参照)には特異成分に対応する減衰項の同時成分ｃ０が未知数として含まれており、実部変換法の演算式((９)〜(11)式参照)には特異成分に対応する剛性項の同時成分ｋ０及び質量項の同時成分ｍ０が未知数として含まれているので、図６(Ｃ)にも示すように、特異成分を分離・除去するという煩雑な処理を行う必要もない（なお、図６(Ｃ)は一例として虚部変換法を示しているが、実部変換法も同様である）。更に、図６(Ｃ)に示すように、虚部変換法又は実部変換法によって求めたインパルス応答を周波数領域へ変換して動的剛性を再現することで、動的剛性の虚数部及び実数部のうちの未知の他方を求めることも可能である。
【００３９】
上記に基づき請求項１記載の発明に係るインパルス応答演算方法は、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するにあたり、物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として前出の(２)式を用いると共に、前記物体の動的剛性の虚数成分ＳＩ(ω)を規定する数式として前出の(６)式を用い（但し、Ｎ１＝Ｎ／２,Ｎ２＝Ｎ／２−１、又は、Ｎ１＝Ｎ２＝(Ｎ−１)／２）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を前出の(７)式及び(８)式へ代入して演算することで、前記インパルス応答を求めることを特徴としている。
【００４０】
請求項１記載の発明は前述の虚部変換法に対応しており、振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を(７)式及び(８)式へ代入して演算することでインパルス応答を求めているので、前述のように、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の虚数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換することができる。
【００４１】
ところで、前述のように虚部変換法では剛性項の同時成分ｋ０と質量項の同時成分ｍ０が算出されず不定となる。このため、インパルス応答演算対象の動的剛性が同時成分ｋ０,ｍ０が元々０又は０に近い値を示す特性である場合には、請求項１記載の発明によって算出したインパルス応答のみを用いて（同時成分ｋ０,ｍ０を０とおいて）時刻歴応答解析等を行うことも可能であるが、インパルス応答演算対象の動的剛性が上記のような特性でない場合は、同時成分ｋ０,ｍ０を別途算出する必要が生ずる。
【００４２】
上記を考慮すると、請求項１記載の発明において、例えば請求項２に記載したように、物体の動的剛性Ｓ(ω)を規定する前出の(１)式に、請求項１記載のインパルス応答演算方法によって求めたインパルス応答の各成分を代入すると共に、(１)式における同時成分ｋ０,ｍ０に０を設定し、振動がＮ３種の周波数のときの物体の動的剛性の実数成分ＳＲ(ωｉ)を各々演算し、演算した実数成分ＳＲ(ωｉ)と、当該実数成分ＳＲ(ωｉ)と同一周波数での動的剛性の実数成分を表すＮ３個の実数部データＤＲ(ωｉ)を
【００４３】
【数８】
2008151643.tif 000010
【００４４】
上記(12),(13)式に代入して同時成分ｋ０,ｍ０を演算することが好ましい。
【００４５】
請求項２記載の発明は、同時成分ｋ０,ｍ０に０を設定して演算したＮ３個の動的剛性の実数成分ＳＲ(ωｉ)と、当該実数成分ＳＲ(ωｉ)と同一周波数での動的剛性の実数成分を表すＮ３個の実数部データＤＲ(ωｉ)の偏差が最小となるように、最小二乗法によって同時成分ｋ０,ｍ０を演算するものである。請求項２記載の発明によれば、Ｎ３≧２であればよく、ごく小数の実数部データＤＲ(ωｉ)から同時成分ｋ０,ｍ０を演算することができる。本願発明者は、虚部変換法によるインパルス応答の演算（図２）で不定となっている同時成分ｋ０,ｍ０を、請求項２記載の発明を適用して演算し、虚部変換法によって算出されたインパルス応答及び請求項２記載の発明を適用して算出された同時成分ｋ０,ｍ０を用いて動的剛性を再現し、再現された動的剛性の実数部を元の動的剛性の実数部と比較することを、Ｎ３＝2,5、減衰率ｈ1＝0,20％の各条件について各々行う解析検討を行った。この解析検討の結果を図４に示す（なお、図４ではＮ３を"Ｎ"と表記している）。
【００４６】
図４(Ａ)に示すように、減衰率ｈ1＝0％の条件ではＮ３の値に拘わらず、再現された動的剛性の実数部は元の動的剛性の実数部におおよそ一致している。また、図４(Ｂ)に示すように、減衰率ｈ1＝20％の場合は、Ｎ３＝２の条件で、元の動的剛性の実数部に対して再現された動的剛性の実数部に若干の偏差が生じているが、Ｎ３＝５の条件では再現された動的剛性の実数部が元の動的剛性の実数部におおよそ一致している。従って、請求項２記載の発明によれば、請求項１記載の発明（虚部変換法）を適用してインパルス応答を演算した場合に不定となる同時成分ｋ０,ｍ０を、比較的少数の実数部データから良好な精度で算出することができる。
【００４７】
請求項３記載の発明に係るインパルス応答演算方法は、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するにあたり、物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として前出の(２)式を用いると共に、前記物体の動的剛性の実数成分ＳＲ(ω)を規定する数式として前出の(９)式を用い（但し、Ｎ１＝Ｎ２＝(Ｎ／２)−１、又は、Ｎ１＝(Ｎ−１)／２,Ｎ２＝((Ｎ−１)／２)−１）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を前出の(10)式及び(11)式へ代入して演算することで、前記インパルス応答を求めることを特徴としている。
【００４８】
請求項３記載の発明は前述の実部変換法に対応しており、振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を(10)式及び(11)式へ代入して演算することでインパルス応答を求めているので、前述のように、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の実数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換することができる。
【００４９】
また、前述のように実部変換法では減衰項の同時成分ｃ０が算出されず不定となる。このため、インパルス応答演算対象の動的剛性が同時成分ｃ０が元々０又は０に近い値を示す特性である場合には、請求項３記載の発明によって算出したインパルス応答のみを用いて（同時成分ｃ０を０とおいて）時刻歴応答解析等を行うことも可能であるが、インパルス応答演算対象の動的剛性が上記のような特性でない場合は、同時成分ｃ０を別途算出する必要が生ずる。
【００５０】
上記を考慮すると、請求項３記載の発明において、例えば請求項４に記載したように、
物体の動的剛性Ｓ(ω)を規定する前出の(１)式に、請求項３記載のインパルス応答演算方法によって求めたインパルス応答の各成分を代入すると共に、(１)式における同時成分ｃ０に０を設定し、振動がＮ４種の周波数のときの物体の動的剛性の虚数成分ＳＩ(ωｉ)を各々演算し、演算した虚数成分ＳＩ(ωｉ)と、当該虚数成分ＳＩ(ωｉ)と同一周波数での動的剛性の虚数成分を表すＮ４個の虚数部データＤＩ(ωｉ)を
【００５１】
【数９】
2008151643.tif 000011
【００５２】
上記(14),(15)式に代入して同時成分ｃ０を演算することが好ましい。
【００５３】
請求項４記載の発明は、同時成分ｃ０に０を設定して演算したＮ４個の動的剛性の虚数成分ＳＩ(ωｉ)と、当該虚数成分ＳＩ(ωｉ)と同一周波数での動的剛性の実数成分を表すＮ４個の虚数部データＤＩ(ωｉ)の偏差が最小となるように、最小二乗法によって同時成分ｃ０を演算するものである。請求項４記載の発明によれば、Ｎ４≧１であればよく、ごく小数の虚数部データＤＩ(ωｉ)から同時成分ｃ０を演算することができる。本願発明者は、実部変換法によるインパルス応答の演算（図２）で不定となっている同時成分ｃ０を、請求項４記載の発明を適用して演算し、実部変換法によって算出されたインパルス応答及び請求項４記載の発明を適用して算出された同時成分ｃ０を用いて動的剛性を再現し、再現された動的剛性の虚数部を元の動的剛性の虚数部と比較することを、Ｎ４＝2,5、減衰率ｈ1＝0,20％の各条件について各々行う解析検討を行った。この解析検討の結果を図５に示す（なお、図５ではＮ４を"Ｎ"と表記している）。
【００５４】
図５(Ａ)に示すように、減衰率ｈ1＝20％の条件ではＮ４の値に拘わらず、再現された動的剛性の虚数部は元の動的剛性の虚数部におおよそ一致している。また、図５(Ａ)に示すように、減衰率ｈ1＝0％の場合は、Ｎ４＝２の条件で、元の動的剛性の虚数部に対して再現された動的剛性の虚数部に若干の偏差が生じているが、Ｎ４＝５の条件では再現された動的剛性の虚数部が元の動的剛性の虚数部におおよそ一致している。従って、請求項４記載の発明によれば、請求項３記載の発明（実部変換法）を適用してインパルス応答を演算した場合に不定となる同時成分ｃ０を、比較的少数の虚数部データから良好な精度で算出することができる。
【００５５】
請求項５記載の発明に係るインパルス応答演算装置は、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置であって、物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として前出の(２)式を用いると共に、前記物体の動的剛性の虚数成分ＳＩ(ω)を規定する数式として前出の(６)式を用い（但し、Ｎ１＝Ｎ／２,Ｎ２＝Ｎ／２−１、又は、Ｎ１＝Ｎ２＝(Ｎ−１)／２）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を前出の(７)式及び(８)式へ代入して演算することで、前記インパルス応答を求める第１演算手段を備えたことを特徴としているので、請求項１記載の発明と同様に、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の虚数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換することができる。
【００５６】
請求項６記載の発明に係るインパルス応答演算装置は、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置であって、物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として前出の(２)式を用いると共に、前記物体の動的剛性の実数成分ＳＲ(ω)を規定する数式として上記(９)式を用い（但し、Ｎ１＝Ｎ２＝(Ｎ／２)−１、又は、Ｎ１＝(Ｎ−１)／２,Ｎ２＝((Ｎ−１)／２)−１）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を前出の(10)式及び(11)式へ代入して演算することで、前記インパルス応答を求める第２演算手段を備えたことを特徴としているので、請求項３記載の発明と同様に、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の実数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換することができる。
【００５７】
請求項７記載の発明に係るインパルス応答演算プログラムは、コンピュータを、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置として機能させるインパルス応答演算プログラムであって、前記コンピュータを、物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として前出の(２)式を用いると共に、前記物体の動的剛性の虚数成分ＳＩ(ω)を規定する数式として前出の(６)式を用い（但し、Ｎ１＝Ｎ／２,Ｎ２＝Ｎ／２−１、又は、Ｎ１＝Ｎ２＝(Ｎ−１)／２）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を前出の(７)式及び(８)式へ代入して演算することで、前記インパルス応答を求める第１演算手段として機能させることを特徴としている。
【００５８】
請求項７記載の発明に係るインパルス応答演算プログラムは、コンピュータを上記の第１演算手段として機能させるためのプログラムであるので、コンピュータが請求項７記載の発明に係るインパルス応答演算プログラムを実行することで、コンピュータが請求項５に記載のインパルス応答演算装置として機能することになり、請求項１及び請求項５記載の発明と同様に、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の虚数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換することができる。
【００５９】
請求項８記載の発明に係るインパルス応答演算プログラムは、コンピュータを、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、前記関係を時間領域で表すインパルス応答へ変換するインパルス応答演算装置として機能させるインパルス応答演算プログラムであって、前記コンピュータを、物体の変位に依存するインパルス応答の同時成分をｋ０、物体の速度に依存するインパルス応答の同時成分をｃ０、物体の加速度に依存するインパルス応答の同時成分をｍ０、物体の変位に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｋｊ、物体の速度に依存するインパルス応答のΔｔ刻みの時間遅れ成分をｃｊ（但しｊは自然数でｋｊ,ｃｊは時刻ｔｊ＝Δｔ・ｊにおける成分ｋ,ｃの値）、時間領域での物体の変位をu(t)、速度をu'(t)、加速度をu"(t)としたときに、前記インパルス応答Ｆ(t)を規定する数式として前出の(２)式を用いると共に、前記物体の動的剛性の実数成分ＳＲ(ω)を規定する数式として前出の(９)式を用い（但し、Ｎ１＝Ｎ２＝(Ｎ／２)−１、又は、Ｎ１＝(Ｎ−１)／２,Ｎ２＝((Ｎ−１)／２)−１）、物体の動的剛性のデータとして、前記振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を前出の(10)式及び(11)式へ代入して演算することで、前記インパルス応答を求める第２演算手段として機能させることを特徴としている。
【００６０】
請求項８記載の発明に係るインパルス応答演算プログラムは、コンピュータを上記の第２演算手段として機能させるためのプログラムであるので、コンピュータが請求項８記載の発明に係るインパルス応答演算プログラムを実行することで、コンピュータが請求項６に記載のインパルス応答演算装置として機能することになり、請求項３及び請求項６記載の発明と同様に、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の実数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換することができる。
【発明の効果】
【００６１】
以上説明したように本発明は、物体を振動させる外力と物体の挙動との関係を時間領域で表すインパルス応答Ｆ(t)を規定する数式として(２)式を用いると共に、前記関係を周波数領域で表す物体の動的剛性の虚数成分ＳＩ(ω)を規定する数式として(６)式を用い、物体の動的剛性のデータとして、振動がＮ種の周波数のときの動的剛性の虚数成分を表すＮ個の虚数部データＤＩ(ω１),…, ＤＩ(ωＮ)を(７)式及び(８)式へ代入して演算することでインパルス応答を求めるようにしたので、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の虚数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換できる、という優れた効果を有する。
【００６２】
また本発明は、物体を振動させる外力と物体の挙動との関係を時間領域で表すインパルス応答Ｆ(t)を規定する数式として(２)式を用いると共に、前記関係を周波数領域で表す物体の動的剛性の実数成分ＳＲ(ω)を規定する数式として(９)式を用い、物体の動的剛性のデータとして、振動がＮ種の周波数のときの動的剛性の実数成分を表すＮ個の実数部データＤＲ(ω１),…, ＤＲ(ωＮ)を(10)式及び(11)式へ代入して演算することでインパルス応答を求めるようにしたので、物体を振動させる外力と物体の挙動との関係を周波数領域で表す動的剛性を、当該動的剛性の実数部のデータから、前記関係を時間領域で表すインパルス応答へ容易かつ精度良く変換できる、という優れた効果を有する。
【発明を実施するための最良の形態】
【００６３】
以下、図面を参照して本発明の実施形態の一例を詳細に説明する。図８には本発明を適用可能なパーソナル・コンピュータ（ＰＣ)１０が示されている。ＰＣ１０は、ＣＰＵ１０Ａ、ＲＯＭ１０Ｂ、ＲＡＭ１０Ｃ及び入出力ポート１０Ｄが、データバス、制御バス、アドレスバス等から成るバス１０Ｅを介して互いに接続されて構成されている。また入出力ポート１０Ｄには、各種の入出力機器として、ＣＲＴ又はＬＣＤから成るディスプレイ１２、キーボード１４、マウス１６、プリンタ１８、ハードディスクドライブ（ＨＤＤ）２０、ＣＤ−ＲＯＭ２２からの情報の読み出しを行うＣＤ−ＲＯＭドライブ２４が各々接続されている。
【００６４】
ＰＣ１０のＨＤＤ２０には、後述する地震応答解析処理を行うための地震応答解析プログラムがインストールされている。この地震応答解析プログラムは、請求項７及び請求項８に記載の発明に係るインパルス応答演算プログラムを含んで構成されている。地震応答解析プログラムをＰＣ１０にインストール（移入）するには幾つかの方法があるが、例えば地震応答解析プログラムをセットアッププログラムと共にＣＤ−ＲＯＭ２２に記録しておき、該ＣＤ−ＲＯＭ２２をＣＤ−ＲＯＭドライブ２４にセットし、ＣＰＵ１０Ａに対して前記セットアッププログラムの実行を指示すれば、ＣＤ−ＲＯＭ２２から地震応答解析プログラムが順に読み出され、読み出された地震応答解析プログラムがＨＤＤ２０に順に書き込まれることで、インパルス応答演算プログラムを含む地震応答解析プログラムのインストールが行われる。
【００６５】
ＰＣ１０は請求項７及び請求項８に記載のコンピュータに対応しており、ＣＰＵ１０Ａが、インパルス応答演算プログラムを含む地震応答解析プログラムを実行することで、請求項５及び請求項６に記載の発明に係るインパルス応答演算装置として機能する。なお、請求項７及び請求項８に記載のコンピュータはＰＣ１０に限られるものではなく、例えばワークステーションであってもよいし、汎用の大型コンピュータであってもよい。
【００６６】
次に本実施形態の作用を説明する。ＰＣ１０のＨＤＤ２０には動的剛性ＤＢ（データベース）が記憶されており、この動的剛性ＤＢには、地盤を振動させる外力（地震動）と地盤の挙動との関係を周波数領域で表す動的剛性のデータが複数種の地盤について各々登録されている。
【００６７】
動的剛性データは、例えば以下のようにして得ることができる。すなわち、動的剛性の演算対象の地盤が、例えば、後述する地震応答解析処理を行う建物（解析対象の建物）の建設を予定している地点、或いは、解析対象の建物の建設地の候補に挙がっている地点における地盤（演算対象の地盤）である場合には、まず前記地点でボーリングを行い、このボーリングによって得られたサンプルに対して所定の試験を行うことで、演算対象の地盤を図１(Ａ)に示したような地盤モデルとしてモデル化し、演算対象の地盤（地盤モデル）のデータとして、地盤の地層構成、各地層の層厚Ｈ、地震波の伝播速度Vs、ポアソン比ν、密度ρ、減衰率ｈを取得する。そして、取得したデータに解析対象の建物の基礎の形状やサイズ等を加えたデータを動的剛性の演算条件とし、当該演算条件に基づき、例えば薄層要素法等の演算方法を適用して演算対象の地盤の動的剛性を演算する。これにより、例として図１(Ｂ),(Ｃ)に示すように、地盤を振動させる外力（地震動）と地盤の挙動との関係を周波数領域で表す動的剛性のデータを得ることができる。
【００６８】
なお、上記の演算方法において、動的剛性の演算条件に含まれる地盤のデータは、サンプルに対して所定の試験を行うことで得られた地盤のデータをそのまま用いることに限られるものではなく、例えば比較的強い地震が起こった後の余震に対する解析対象の建物の挙動を解析したい等の場合には、比較的強い地震により演算対象の地盤の特性が変化することを想定し、所定の試験によって得られたデータに対し、特性変化に相当する値の変更を加えたデータを演算条件データとして用いてもよい。また、動的剛性のデータは演算によって求めることに限られるものではなく、実験を行って求めることも可能である。
【００６９】
また、動的剛性ＤＢに登録される動的剛性データは、現実の地盤（をモデル化した地盤モデル）の動的剛性を表すデータに限られるものではない。地震応答解析に対しては、例えば動的剛性が特定の特性を示す（より詳しくは、動的剛性の虚数部又は実数部が特定の特性を示す）仮想的な地盤上に建物を建設した場合の当該建物の地震時の挙動を確認したい、等のニーズが存在している。このため動的剛性ＤＢには、上記のようなニーズを満たす目的で創出された仮想的な地盤モデルの動的剛性を表す動的剛性データも登録される。この種の動的剛性データは、従来、動的剛性の虚数部及び実数部の一方が特定の特性を示すように虚数部及び実数部の一方の特性が設定された後に、対応するインパルス応答を演算して地震応答解析（時刻歴応答解析）を行うために、動的剛性の虚数部及び実数部がヒルベルと変換対を形成するように（時間領域で因果律を満たすように）、動的剛性の虚数部及び実数部の他方の特性が演算によって設定されていた。しかし本実施形態では、本発明に係る虚部変換法又は実部変換法を適用してインパルス応答が演算されることを前提に、動的剛性の虚数部及び実数部の一方の特性のみが設定され、当該特性を表す虚数部又は実数部のデータのみが動的剛性データとして動的剛性ＤＢに登録されている。
【００７０】
次に、解析対象の建物の地震応答解析の実行を所望しているオペレータによってキーボード１４又はマウス１６が操作されることで、地震応答解析プログラムの実行が指示されると共に、動的剛性ＤＢに登録される動的剛性データのうち地震応答解析に用いる動的剛性データ（処理対象の動的剛性データ）が指定され、これに伴ってＰＣ１０のＣＰＵ１０Ａで実行される地震応答解析処理について、図８のフローチャートを参照して説明する。
【００７１】
地震応答解析処理では、まずステップ１００において、オペレータによって指定された処理対象の動的剛性データを動的剛性ＤＢから読み出し、読み出した動的剛性データをメモリ又はＨＤＤ２０の所定の領域に記憶させる。本実施形態に係る地震応答解析処理では、インパルス応答の演算方法として旧変換法、虚部変換法及び実部変換法の３種類の演算方法が用意されており（ＨＤＤ２０には、動的剛性をインパルス応答へ変換するための各演算方法に対応する連立方程式が各々記憶されている）、ステップ１０２では読み出した処理対象の動的剛性データに基づいて、処理対象の動的剛性データに適したインパルス応答の演算方法を判定し、ステップ１０２の判定の結果に応じて次のステップ１０４で分岐する。
【００７２】
なお、ステップ１０２の判定は、具体的には、例えば読み出した処理対象の動的剛性データに虚数部及び実数部のデータが各々含まれており、かつ虚数部のデータの情報量と実数部のデータの情報量に大きな偏りも無い場合は、処理対象の動的剛性データに適したインパルス応答の演算方法として旧変換法が選択され、読み出した処理対象の動的剛性データが虚数部のデータのみから構成されている場合や、処理対象の動的剛性データに虚数部及び実数部のデータが各々含まれているものの、実数部のデータの情報量が虚数部のデータの情報量よりも大幅に少ない場合には、処理対象の動的剛性データに適したインパルス応答の演算方法として虚部変換法が選択され、処理対象の動的剛性データが実数部のデータのみから構成されている場合や、処理対象の動的剛性データに虚数部及び実数部のデータが各々含まれているものの、虚数部のデータの情報量が実数部のデータの情報量よりも大幅に少ない場合には、処理対象の動的剛性データに適したインパルス応答の演算方法として実部変換法が選択されるように行うことができる。
【００７３】
ステップ１０２において、処理対象の動的剛性データに適したインパルス応答の演算方法として旧変換法が選択された場合には、ステップ１０４からステップ１０６へ移行し、処理対象の動的剛性データから、予め設定された演算対象の周波数範囲内のＮ種の周波数（Ｎ種の角振動数ω１,…,ωＮ）における動的剛性の値を表すＮ個の複素データＤ(ω１),…,Ｄ(ωＮ)を各々抽出し、抽出した複素データをメモリ又はＨＤＤ２０に記憶させる。なお、演算対象の周波数範囲としては例えば０〜20(Hz)の範囲が好適であるが、他の範囲であってもよい。また、複素データの抽出を行うＮ種の周波数は、例えば演算対象の周波数範囲の上限に相当する周波数（例えば演算対象の周波数範囲が０〜20(Hz)であれば、上限周波数である20(Hz)）を含むように設定することができる。
【００７４】
また、地盤の動的剛性のデータから抽出した複素データはインパルス応答の演算に用いられ、この演算により時刻ｔ＝Δｔ・ｊ（ｊ=0,1,2,…,n)の各時刻における地盤のインパルス応答を表すインパルス応答データが得られるが、得られるインパルス応答データの個数は演算に用いる複素データの個数に応じて定まり（すなわちｊの最大値ｎは複素データの個数−１(＝Ｎ−１))、得られるインパルス応答データによって表される地盤のインパルス応答の時刻範囲も演算に用いる複素データの個数に応じて定まる（例えば複素データの個数が21個、Δｔ＝0.05秒とすると、ｔn＝Δｔ・ｊmax＝0.05×20＝１秒となり、時刻ｔ＝０〜１秒の時刻範囲の地盤のインパルス応答を表す21個のインパルス応答データが得られる）ことになるので、地盤の動的剛性から抽出する複素データの個数（複素データの抽出を行う周波数の種類数）は、地盤のインパルス応答を算出すべき時刻範囲の長さも勘案して予め定めておくことができる。
【００７５】
次のステップ１０８では、動的剛性をインパルス応答へ変換するための旧変換法の連立方程式（前出の(４),(５)式）をＨＤＤ２０から読み出し、読み出した連立方程式に、ステップ１０６で抽出したＮ個の複素データＤ(ω１),…,Ｄ(ωＮ)を代入し、この連立方程式の解を求めることで、地盤のインパルス応答を表すインパルス応答データを、旧変換法によりΔｔ刻みで演算する。この演算により、地盤のインパルス応答を表すインパルス応答データとして、剛性項（地盤の変位に依存するインパルス応答）の同時成分ｋ０、減衰項（地盤の速度に依存するインパルス応答）の同時成分ｃ０、質量項（地盤の加速度に依存するインパルス応答）の同時成分ｍ０のデータが得られると共に、剛性項の時間遅れ成分ｋｊのデータがΔｔ刻みでｎ個（ｎ＝Ｎ−１）得られ、減衰項の時間遅れ成分ｃｊのデータがΔｔ刻みでｎ−１個得られることになる。そして、得られたインパルス応答データをメモリ又はＨＤＤ２０に一旦記憶するとステップ１４２へ移行する。
【００７６】
一方、ステップ１０２において、処理対象の動的剛性データに適したインパルス応答の演算方法として虚部変換法が選択された場合には、ステップ１０４からステップ１１０へ移行し、ステップ１１０以降で虚部変換法を適用してインパルス応答を演算する。すなわち、ステップ１１０では処理対象の動的剛性データから、演算対象の周波数範囲内のＮ種の周波数（Ｎ種の角振動数ω１,…,ωＮ）における動的剛性の虚数成分の値を表すＮ個の虚数部データＤＩ(ω１),…,ＤＩ(ωＮ)を各々抽出し、抽出した虚数部データをメモリ又はＨＤＤ２０に記憶させる。なお、ステップ１１０における虚数部データの抽出において、演算対象の周波数範囲については旧変換法（ステップ１０６における複素データの抽出）と同一でよいが、虚部変換法を適用したインパルス応答の演算において、旧変換法を適用したインパルス応答の演算と同程度の演算精度を確保するために、例えばΔｔの値を旧変換法におけるΔｔの値の1/2にすることで、抽出する虚数部データの個数Ｎを旧変換法における個数Ｎの２倍程度の値とすることが望ましい。
【００７７】
次のステップ１１２では、動的剛性をインパルス応答へ変換するための虚部変換法の連立方程式（前出の(７),(８)式）をＨＤＤ２０から読み出し、読み出した連立方程式に、ステップ１１０で抽出したＮ個の虚数部データＤＩ(ω１),…,ＤＩ(ωＮ)を代入し、この連立方程式の解を求めることで、地盤のインパルス応答を表すインパルス応答データを、虚部変換法によりΔｔ刻みで演算する。この演算により、地盤のインパルス応答を表すインパルス応答データとして、減衰項の同時成分ｃ０のデータが得られると共に、剛性項の時間遅れ成分ｋｊのデータがΔｔ刻みでｎ個（ｎ＝Ｎ−１）得られ、減衰項の時間遅れ成分ｃｊのデータがΔｔ刻みでｎ−１個得られる。そして、得られたインパルス応答データはメモリ又はＨＤＤ２０に一旦記憶される。なお、上記のステップ１１０，１１２は請求項５，７に記載の第１演算手段に対応している。
【００７８】
ところで、虚部変換法の連立方程式（前出の(７),(８)式）では剛性項の同時成分ｋ０及び質量項の同時成分ｍ０が未知数として設定されていないため、ステップ１１２のインパルス応答の演算でも剛性項の同時成分ｋ０及び質量項の同時成分ｍ０は算出されず、これらの値が不定となる。このため、次のステップ１１４〜１２４では剛性項の同時成分ｋ０及び質量項の同時成分ｍ０を設定する処理を行う。
【００７９】
すなわち、まずステップ１１４では処理対象の動的剛性データに実数部のデータが含まれているか否か判定する。判定が肯定された場合はステップ１１６へ移行し、処理対象の動的剛性データから、演算対象の周波数範囲内のＮ３種の周波数（Ｎ３種の角振動数ω１,…,ωＮ３）における動的剛性の実数成分の値を表すＮ３個の実数部データＤＲ(ω１),…,ＤＲ(ωＮ３)を抽出し、抽出したＮ３個の実数部データをメモリ又はＨＤＤ２０に記憶させた後にステップ１２２へ移行する。なお、Ｎ３は２以上の値であればよいが、より大きな値（例えば５程度の値）とした方が後述する同時成分ｋ０,ｍ０の演算精度を確保するために望ましく、処理対象の動的剛性データに含まれる実数部のデータの情報量等に応じてＮ３の値を切り替えるようにしてもよい。
【００８０】
また、処理対象の動的剛性データに実数部のデータが含まれていない場合はステップ１１４の判定が否定されてステップ１１８へ移行する。この場合、処理対象の動的剛性データから同時成分ｋ０,ｍ０を演算することは困難であるので、ステップ１１８では、ディスプレイ１２に所定のメッセージを表示させる等によってオペレータに対処を問い合わせ、オペレータによってキーボード１４又はマウス１６が操作されることで対処が指示されると、指示された対処を認識して処理を分岐する。本実施形態において、虚部変換法を適用してインパルス応答を演算した処理対象の動的剛性データに実数部のデータが含まれていない場合の対処としては「実数部データの追加入力」と「同時成分ｋ０,ｍ０を０に設定」があり、上記対処として「実数部データの追加入力」がオペレータによって指示された場合はステップ１２０へ移行し、オペレータによって引き続きキーボード１４が操作されることで追加入力されたＮ３個の実数部データＤＲ(ω１),…,ＤＲ(ωＮ３)をメモリ又はＨＤＤ２０に記憶させる。
【００８１】
次のステップ１２２では、ステップ１１２の演算によって得られたインパルス応答データと、ステップ１１６又はステップ１２０によって得られたＮ３個の実数部データに基づいて同時成分ｋ０,ｍ０を演算する。すなわち、まずステップ１１２の演算によって得られたインパルス応答データ（減衰項の同時成分ｃ０及び時間遅れ成分ｃｊ、剛性項の時間遅れ成分ｋｊ）を前出の(１)式に代入すると共に、(１)式における同時成分ｋ０,ｍ０に０を設定し、Ｎ３個の実数部データＤＲ(ω１),…,ＤＲ(ωＮ３)と同一のＮ３種の周波数における物体の動的剛性の実数成分ＳＲ(ω１),…,ＳＲ(ωＮ３)を各々演算する。そして、上記演算によって得られたＮ３個の実数成分ＳＲ(ω１),…,ＳＲ(ωＮ３)と、ステップ１１６又はステップ１２０によって得られたＮ３個の実数部データＤＲ(ω１),…,ＤＲ(ωＮ３)を前出の (12),(13)式に代入して同時成分ｋ０,ｍ０を演算し、算出された同時成分ｋ０,ｍ０をメモリ又はＨＤＤ２０に記憶させた後にステップ１４２へ移行する。なお、上述したステップ１１６、１２０、１２２は請求項２記載の発明に対応している。
【００８２】
また、虚部変換法を適用してインパルス応答を演算した処理対象の動的剛性データに実数部のデータが含まれていない場合の対処として「同時成分ｋ０,ｍ０を０に設定」がオペレータによって指示された場合はステップ１１８からステップ１２４へ移行し、同時成分ｋ０,ｍ０として各々０を設定し、０を設定した同時成分ｋ０,ｍ０をメモリ又はＨＤＤ２０に記憶させた後にステップ１４２へ移行する。以上の処理により、虚部変換法を適用したインパルス応答の演算において、時刻歴応答解析に必要な全てのインパルス応答（剛性項の同時成分ｋ０及び時間遅れ成分ｋｊ、減衰項の同時成分ｃ０及び時間遅れ成分ｃｊ、質量項の同時成分ｍ０）が得られることになる。
【００８３】
また、ステップ１０２において、処理対象の動的剛性データに適したインパルス応答の演算方法として実部変換法が選択された場合には、ステップ１０４からステップ１２６へ移行し、ステップ１２６以降で実部変換法を適用してインパルス応答を演算する。すなわち、ステップ１２６では処理対象の動的剛性データから、演算対象の周波数範囲内のＮ種の周波数（Ｎ種の角振動数ω１,…,ωＮ）における動的剛性の実数成分の値を表すＮ個の実数部データＤＲ(ω１),…,ＤＲ(ωＮ)を各々抽出し、抽出した実数部データをメモリ又はＨＤＤ２０に記憶させる。なお、ステップ１２６における実数部データの抽出においても、演算対象の周波数範囲については旧変換法（ステップ１０６における複素データの抽出）と同一でよいが、実部変換法を適用したインパルス応答の演算において、旧変換法を適用したインパルス応答の演算と同程度の演算精度を確保するために、例えばΔｔの値を旧変換法におけるΔｔの値の1/2にすることで、抽出する実数部データの個数Ｎを旧変換法における個数Ｎの２倍程度の値とすることが望ましい。
【００８４】
次のステップ１２８では、動的剛性をインパルス応答へ変換するための実部変換法の連立方程式（前出の(10),(11)式）をＨＤＤ２０から読み出し、読み出した連立方程式に、ステップ１２６で抽出したＮ個の実数部データＤＲ(ω１),…,ＤＲ(ωＮ)を代入し、この連立方程式の解を求めることで、地盤のインパルス応答を表すインパルス応答データを、実部変換法によりΔｔ刻みで演算する。この演算により、地盤のインパルス応答を表すインパルス応答データとして、剛性項の同時成分ｋ０のデータ及び質量項の同時成分ｍ０のデータが得られると共に、剛性項の時間遅れ成分ｋｊのデータがΔｔ刻みでｎ個（ｎ＝Ｎ−１）得られ、減衰項の時間遅れ成分ｃｊのデータがΔｔ刻みでｎ−１個得られる。そして、得られたインパルス応答データはメモリ又はＨＤＤ２０に一旦記憶される。なお、上記のステップ１２６，１２８は請求項６，８に記載の第２演算手段に対応している。
【００８５】
ところで、実部変換法の連立方程式（前出の(10),(11)式）では減衰項の同時成分ｃ０が未知数として設定されていないため、ステップ１２８のインパルス応答の演算でも減衰項の同時成分ｃ０は算出されず、減衰項の同時成分ｃ０の値が不定となる。このため、次のステップ１３０〜１４０では減衰項の同時成分ｃ０を設定する処理を行う。
【００８６】
すなわち、まずステップ１３０では処理対象の動的剛性データに虚数部のデータが含まれているか否か判定する。判定が肯定された場合はステップ１３２へ移行し、処理対象の動的剛性データから、演算対象の周波数範囲内のＮ４種の周波数（Ｎ４種の角振動数ω１,…,ωＮ４）における動的剛性の虚数成分の値を表すＮ４個の虚数部データＤＩ(ω１),…,ＤＩ(ωＮ４)を抽出し、抽出したＮ４個の虚数部データをメモリ又はＨＤＤ２０に記憶させた後にステップ１３８へ移行する。なお、Ｎ４は１以上の値であればよいが、より大きな値（例えば５程度の値）とした方が後述する同時成分ｃ０の演算精度を確保するために望ましく、処理対象の動的剛性データに含まれる虚数部のデータの情報量等に応じてＮ４の値を切り替えるようにしてもよい。
【００８７】
また、処理対象の動的剛性データに虚数部のデータが含まれていない場合はステップ１３０の判定が否定されてステップ１３４へ移行する。この場合、処理対象の動的剛性データから同時成分ｃ０を演算することは困難であるので、ステップ１３４では、ディスプレイ１２に所定のメッセージを表示させる等によってオペレータに対処を問い合わせ、オペレータによってキーボード１４又はマウス１６が操作されることで対処が指示されると、指示された対処を認識して処理を分岐する。本実施形態において、実部変換法を適用してインパルス応答を演算した処理対象の動的剛性データに虚数部のデータが含まれていない場合の対処としては「虚数部データの追加入力」と「同時成分ｃ０を０に設定」があり、上記対処として「虚数部データの追加入力」がオペレータによって指示された場合はステップ１３６へ移行し、オペレータによって引き続きキーボード１４が操作されることで追加入力されたＮ４個の虚数部データＤＩ(ω１),…,ＤＩ(ωＮ４)をメモリ又はＨＤＤ２０に記憶させる。
【００８８】
次のステップ１３８では、ステップ１２８の演算によって得られたインパルス応答データと、ステップ１３２又はステップ１３６によって得られたＮ４個の虚数部データに基づいて同時成分ｃ０を演算する。すなわち、まずステップ１２８の演算によって得られたインパルス応答データ（剛性項の同時成分ｋ０及び時間遅れ成分ｋｊ、減衰項の時間遅れ成分ｃｊ、質量項の同時成分ｍ０）を前出の(１)式に代入すると共に、(１)式における同時成分ｃ０に０を設定し、Ｎ４個の虚数部データＤＩ(ω１),…,ＤＩ(ωＮ４)と同一のＮ４種の周波数における物体の動的剛性の虚数成分ＳＩ(ω１),…,ＳＩ(ωＮ４)を各々演算する。そして、上記演算によって得られたＮ４個の虚数成分ＳＩ(ω１),…,ＳＩ(ωＮ４)と、ステップ１３２又はステップ１３６によって得られたＮ４個の虚数部データＤＩ(ω１),…,ＤＩ(ωＮ４)を前出の(14),(15)式に代入して同時成分ｃ０を演算し、算出された同時成分ｃ０をメモリ又はＨＤＤ２０に記憶させた後にステップ１４２へ移行する。なお、上述したステップ１３２、１３６、１３８は請求項４記載の発明に対応している。
【００８９】
また、実部変換法を適用してインパルス応答を演算した処理対象の動的剛性データに虚数部のデータが含まれていない場合の対処として「同時成分ｃ０を０に設定」がオペレータによって指示された場合はステップ１３４からステップ１４０へ移行し、同時成分ｃ０として０を設定し、０を設定した同時成分ｃ０をメモリ又はＨＤＤ２０に記憶させた後にステップ１４２へ移行する。以上の処理により、実部変換法を適用したインパルス応答の演算においても、時刻歴応答解析に必要な全てのインパルス応答（剛性項の同時成分ｋ０及び時間遅れ成分ｋｊ、減衰項の同時成分ｃ０及び時間遅れ成分ｃｊ、質量項の同時成分ｋ０,ｍ０）が得られることになる。
【００９０】
そしてステップ１４２では、メモリ又はＨＤＤ２０に記憶されているインパルス応答データ（剛性項の同時成分ｋ０及び時間遅れ成分ｋｊ、減衰項の同時成分ｃ０及び時間遅れ成分ｃｊ、質量項の同時成分ｋ０,ｍ０）読み出し、読み出したインパルス応答データを用いて解析対象の建物の時刻歴地震応答解析を行う。この時刻歴地震応答解析に用いるインパルス応答データは、処理対象の動的剛性データに実数部のデータが含まれていなかったり、実数部のデータの情報量が虚数部のデータの情報量よりも大幅に少ない場合は虚部変換法が適用され、処理対象の動的剛性データに虚数部のデータが含まれていなかったり、虚数部のデータの情報量が実数部のデータの情報量よりも大幅に少ない場合は実部変換法が適用されて演算されたデータであるので、上記のような場合や、処理対象の動的剛性データが表す地盤の動的剛性が、比較的大きな履歴減衰成分を含んでいる場合にも、地震動と地盤の挙動との関係を時間領域で精度良く表すデータであり、解析対象の建物の時刻歴地震応答解析を精度良く行うことができる。
【００９１】
なお、上記では地盤の動的剛性から本発明を適用して地盤のインパルス応答を求め、求めたインパルス応答を解析対象の建物の時刻歴地震応答解析に用いる例を説明したが、これに限定されるものではなく、例えば建物の制振装置として用いられている粘性ダンパー等の他の物体の動的剛性をインパルス応答へ変換する際に本発明を適用することも可能であることは言うまでもない。
【００９２】
また、上記では求めたインパルス応答を時刻歴応答解析に用いることを前提として、物体の動的剛性から物体のインパルス応答を演算する場合に本発明を適用した態様を説明したが、これに限定されるものではなく、例えば虚数部及び実数部のうちの一方が既知で、他方が未知又は他方に関しての情報が不足している状態の動的剛性が存在しており、この動的剛性の虚数部及び実数部のうちの前記他方の情報を取得したい等の場合に本発明を適用する（具体的には、本発明に係る虚部変換法又は実部変換法を適用してインパルス応答を演算し、当該演算で算出されない同時成分(ｋ０,ｍ０又はｃ０)に値を設定した後に、この同時成分を含むインパルス応答を周波数領域へ変換して動的剛性を再現することで、前記動的剛性の虚数部及び実数部のうちの未知の他方を求める(図６(Ｃ)も参照))ことも可能である。
【実施例】
【００９３】
次に、本発明の効果を確認するために本願発明者が実施した解析検討の結果について説明する。減衰定数ｈが振動数非依存特性を示す物体の減衰を正確に表す履歴減衰モデル、すなわち履歴減衰を有する動的剛性Ｓ(ω)は、次の(16)式で表される。なお、(16)式においてＫ０は剛性、ｈは減衰定数、ｉは虚数単位である。
Ｓ(ω)＝Ｋ０(１＋２ｈ・ｉ)  …(16)
ここで、(16)式における虚数単位ｉに代えて(18)式で表される単位虚数関数Ｚ(ω)を用いれば、先の(16)式は(17)式で表される。
Ｓ(ω)＝Ｋ０(１＋２ｈ・Ｚ(ω))  …(17)
Ｚ(ω)＝ＺＲ(ω＋ＺＩ(ω)・ｉ  …(18)
但し、ＺＲ(ω)は単位虚数関数Ｚ(ω)の実数部、ＺＩ(ω)は単位虚数関数Ｚ(ω)の虚数部であり、それぞれ以下の(19)式で表される。
【００９４】
【数１０】
2008151643.tif 000012
【００９５】
但し、上記(18)式で表される単位虚数関数Ｚ(ω)は、実数部ＺＲ(ω)と虚数部ＺＩ(ω)がヒルベルト変換対を形成しないために時間領域で因果律を満たさず、時間領域への厳密な変換は不可能である。
【００９６】
このため本願発明者は、地震応答解析等の時刻歴応答解析における解析対象の振動数範囲が限られていることに着目し、『中村  尚弘，「因果的履歴減衰の実用的評価法−複素減衰の簡便な時間領域近似表現の試み−」，日本建築学会構造系論文報告集，２００５年１０月，第５９６号，ｐ．３３−３９』（以下、非特許文献２と称する）において、図９に示す正則成分Ｚ'ＩＲ(ω)と特異成分Ｚ'ＩＳ(ω)の和で表され、階段状に変化すると共に一定の振動数範囲(検討振動数範囲)内における値が単位虚数関数Ｚ(ω)と等しい関数を虚部Ｚ'Ｉ(ω)とおき、この虚部Ｚ'Ｉ(ω)の正則成分Ｚ'ＩＲ(ω)からヒルベルト変換を用いて因果的な実部Ｚ'Ｒ(ω)を算定することで、
Ｚ'(ω)＝Ｚ'Ｒ(ω)＋Ｚ'Ｉ(ω)
なる因果的単位虚数関数Ｚ'(ω)を提案している（図１０も参照）。
【００９７】
本願発明者が実施した解析検討では、インパルス応答演算対象の動的剛性として上記の因果的単位虚数関数Ｚ'(ω)を用い、因果的単位虚数関数Ｚ'(ω)に対し、複素データ（実数部及び虚数部）のサンプリングを行って時間領域へ変換する変換法（旧変換法）と、虚数部データのみのサンプリングを行って時間領域へ変換する変換法（本発明に係る虚部変換法）を各々適用してインパルス応答値を各々求めた。また、旧変換法については、サンプリングを行う複素データの数（データ点の数）として「９点」(旧変換法１と表記)及び「１９点」(旧変換法２と表記)を各々適用して演算を行った。各変換法における演算条件の詳細を次の表１に示し、各変換法によって算出されたインパルス応答の剛性項を図１１(Ａ)に、減衰項を図１１(Ｂ)に示す。
【００９８】
【表１】
2008151643.tif 000013
【００９９】
図１１(Ａ),(Ｂ)に示すように、剛性項の同時成分ｋ０は各変換法共におよそ０になっており、減衰項については各変換法共に同時成分ｃ０のみで構成されている。従って、因果的単位虚数関数Ｚ'(ω)のインパルス応答は、実質的に、剛性項の時間遅れ成分と減衰項の同時成分、すなわち現在の速度と過去の変位に依存する量（成分）で構成されていると見なすことができる。そして図１１(Ａ),(Ｂ)からも明らかなように、虚部変換法によって算出されたインパルス応答は、旧変換法１、２によって算出されたインパルス応答に非常に近い値を示している。
【０１００】
また、図１１(Ｃ),(Ｄ)には、各変換法によって算出されたインパルス応答を用い、前出の(１)式に基づいて再現した複素関数が、元の因果的単位虚数関数Ｚ'(ω)のデータ点（旧変換法２におけるデータ点）と共に示されている。図１１(Ｃ),(Ｄ)に示すように、虚部変換法によって算出されたインパルス応答から再現された複素関数は、旧変換法１によって算出されたインパルス応答から再現された複素関数と、旧変換法２によって算出されたインパルス応答から再現された複素関数の間に相当する値を示しており、データ点に対する偏差もおおよそ小さい。
【０１０１】
なお、虚部変換法によって算出されたインパルス応答から再現された複素関数のうちの虚数部(図では"Imag"と表記)は、図１１(Ｄ)に示すように、0.5(Hz)付近でデータ点に対する偏差が比較的大きくなっているが、これは虚数部データのみを用いてインパルス応答を演算していることの影響と考えられ、偏差自体も顕著な大きさではない。このように、本願発明者が行った解析検討により、本発明に係る虚部変換法（や実部変換法）を適用した場合にも、虚数部データ及び実数部データを各々用いる旧変換法と同等の演算精度でインパルス応答を算出できることが確認された。
【図面の簡単な説明】
【０１０２】
【図１】(Ａ)は本願発明者による解析検討に用いたモデルを示すイメージ図、(Ｂ)は(Ａ)のモデルから演算によって求めた地盤の動的剛性の一例を示す線図である。
【図２】本発明に係る虚部変換法及び実部変換法によるインパルス応答演算結果の一例を示す線図である。
【図３】虚部変換法及び実部変換法を適用して演算したインパルス応答から再現した動的剛性の一例を示す線図である。
【図４】虚部変換法を適用して演算したインパルス応答及び別途演算した同時成分ｋ０,ｍ０を用いて再現した動的剛性の一例を示す線図である。
【図５】実部変換法を適用して演算したインパルス応答及び別途演算した同時成分ｃ０を用いて再現した動的剛性の一例を示す線図である。
【図６】(Ａ)は特許文献１,２に記載の方法、(Ｂ)は非特許文献１に記載の方法、(Ｃ)は本発明に係る虚部変換法の概略を各々示すフローチャートである。
【図７】本実施形態に係るＰＣの概略構成を示すブロック図である。
【図８】地震応答解析処理の内容を示すフローチャートである。
【図９】因果的単位虚数関数の虚数部全体の特性を示す線図である。
【図１０】検討振動数範囲内での因果的単位虚数関数の特性を示す線図である。
【図１１】本願発明者が実施した解析検討において、(Ａ),(Ｂ)は算出されたインパルス応答、(Ｃ),(Ｄ)は再現された因果的単位虚数関数を各々示す線図である。
【符号の説明】
【０１０３】
１０  ＰＣ
１２  ディスプレイ
１４  キーボード
１６  マウス
２０  ＨＤＤ
【図２】
2008151643.tif 000036
【図３】
2008151643.tif 000037
【図４】
2008151643.tif 000038
【図５】
2008151643.tif 000039
【図６】
2008151643.tif 000040
【図７】
2008151643.tif 000041
【図８】
2008151643.tif 000042
【図９】
2008151643.tif 000043
【図１０】
2008151643.tif 000044
【図１１】
2008151643.tif 000045
【図１】
2008151643.tif 000046