JP6988050B2 - データ推定方法 - Google Patents

Description

本発明は、データ推定方法に関する。

従来、地震による構造物の被害状況を地震発生後に推定することができる構造物の損傷推定システムが知られている（例えば、特許文献１を参照）。特許文献１に記載の技術では、特性の異なる複数の模擬地震波と当該複数の模擬地震波のそれぞれに対する構造物の動的応答とを関連付けて記憶させ、その複数の模擬地震波の中から構造物の近傍で観測された実地震波に近似する模擬地震波を検索し、検索された模擬地震波に対する構造物の動的応答を表示する。

また、物性が急激に変化する現象に等価線形解析法を適用する動的応答解析方法が知られている（例えば、特許文献２を参照）。特許文献２に記載の技術では、入力地震動の分割が行われ、液状化前の第１グループと、液状化後の第２グループとにグループ分けされる。第１グループは、液状化前の動的変形特性を用いた等価線形解析が行われる。第２グループは、液状化後の動的変形特性を用い、入力地震動が、液状化発生時刻ｔから最後までとされる等価線形解析が行われる。そして、第１および第２グループの等価線形解析の解が得られると、解析結果(応答波形)を重ね合わせて手順が終了する。

特開2004-069598号公報 特開2001-116651号公報

しかし、上記特許文献１及び上記特許文献２に記載の技術では、逐次非線形性を有する非線形入出力システムを解析対象とする場合に、観測データから入力データを適切に推定することはできない。

例えば、非線形入出力システムに対応する解析モデルが地盤モデルである場合、地震時の地盤の挙動は逐次非線形性を示し、地盤の物性値は時々刻々変化する。

上記特許文献２に記載の技術においては、液状化前と液状化後とで等価線形解析が行われるが、液状化の前後に関わらず地盤の物性値は時々刻々変化するため、地盤モデルへの入力データである入力地震動を適切に推定することができない。

本発明は上記事実を考慮して、逐次非線形性を有する非線形入出力システムを解析対象とする場合に、観測データから非線形入出力システムの入力データを推定することを目的とする。

上記目的を達成するために、本発明のデータ推定方法は、入力ベクトルが入力されたときに出力ベクトルを出力する非線形入出力システムにおいて、任意の前記出力ベクトルが観測されたときの前記入力ベクトルを推定するデータ推定方法であって、前記入力ベクトルの候補を表す入力候補ベクトルを、前記非線形入出力システムへ入力したときに出力される前記出力ベクトルを表す出力候補ベクトルを計算し、観測された任意の前記出力ベクトルを目標ベクトルとして、前記目標ベクトルと前記出力候補ベクトルとの間の差分を表す差分ベクトルを計算し、前記入力候補ベクトルに含まれる複数の要素の各々について、前記要素を変化させたベクトルである入力候補変化ベクトルを計算し、複数の前記入力候補変化ベクトルの各々について、前記入力候補変化ベクトルを前記非線形入出力システムへ入力したときの前記出力ベクトルを表す出力候補変化ベクトルを計算し、複数の前記出力候補変化ベクトルの各々について、前記出力候補変化ベクトルと前記出力候補ベクトルとの間の差分から求まる変化率を表す係数ベクトルを計算し、複数の係数ベクトルを列ベクトルとした行列を表す応答行列を生成し、前記応答行列の行列分解結果から得られる複数のモードから特定のモードを選択し、選択された前記モードに対応する部分行列を用いて、前記応答行列の一般化逆行列を生成し、前記一般化逆行列と前記差分ベクトルとに基づいて、前記入力候補ベクトルに対する補正量を表す補正ベクトルを計算し、前記入力候補ベクトルと前記補正ベクトルとに基づいて、新たな前記入力候補ベクトルを生成し、前記差分ベクトルに対応する値が予め定められた閾値未満となるまで、前記出力候補ベクトルの計算、前記差分ベクトルの計算、前記入力候補変化ベクトルの計算、前記出力候補変化ベクトルの計算、前記応答行列の生成、前記一般化逆行列の生成、及び前記新たな前記入力候補ベクトルの生成の各処理を繰り返し、前記差分ベクトルに対応する値が前記閾値未満となった場合に、前記入力候補ベクトルを、観測された任意の前記出力ベクトルに対応する前記入力ベクトルとして取得する、処理をコンピュータに実行させる。これにより、逐次非線形性を有する非線形入出力システムを解析対象とする場合に、観測データから非線形入出力システムの入力データを推定することができる。

本発明のデータ推定方法は、前記応答行列の行列分解結果から得られる複数のモードから特定のモードを選択する際に、前記モードの寄与度を表す値が予め設定されたモードに関する閾値以上である前記モードを選択し、選択された前記モードに対応する部分行列を用いて、前記応答行列の一般化逆行列を生成するようにしてもよい。これにより、観測データから非線形入出力システムの入力データを推定する際に、繰り返し計算の収束時間を短くすることができる。

本発明のデータ推定方法は、前記応答行列の行列分解結果から得られる複数のモードから特定のモードを選択する際に、前記モードの寄与度を表す値が大きい方から予め設定された個数の前記モードを選択し、選択された前記モードに対応する部分行列を用いて、前記応答行列の一般化逆行列を生成するようにしてもよい。これにより、観測データから非線形入出力システムの入力データを推定する際に、繰り返し計算の収束時間を短くすることができる。

本発明のデータ推定方法は、前記非線形入出力システムは、入力ベクトルが入力されたときに各々異なる観測点での出力ベクトルを出力するものであって、各々異なる観測点で観測された任意の複数の前記出力ベクトルを、複数の前記目標ベクトルとして、各々異なる観測点で任意の複数の前記出力ベクトルが観測されたときの前記入力ベクトルを推定するようにしてもよい。これにより、観測データから非線形入出力システムの入力データを推定する場合に、複数の異なる観測データから、非線形入出力システムの入力データを推定することができる。

また、本発明のデータ推定方法は、前記入力ベクトルは、複数種類のベクトルを含むようにしてもよい。これにより、観測データから非線形入出力システムの入力データを推定する場合に、複数の異なる種類の入力データを推定することができる。

本発明のデータ推定方法における前記非線形入出力システムは、ニューロンの数が予め設定されたニューラルネットワークであり、前記入力ベクトルは、前記ニューラルネットワークのパラメータであり、前記出力ベクトルは、前記ニューラルネットワークの出力であり、前記差分ベクトルは、学習用の入力データが前記ニューラルネットワークに入力された際の前記出力ベクトルと、前記学習用の入力データに対する正解の前記出力ベクトルを表す前記目標ベクトルとの間の差分であるようにしてもよい。これにより、非線形入出力システムの一例であるニューラルネットワークのパラメータを推定することができる。

また、本発明のデータ推定方法は、前記各処理を繰り返す際に、前記モードに関する閾値を減少させることにより、選択される前記モードの数を増加させるようにしてもよい。これにより、観測データから非線形入出力システムの入力データを推定する際に、繰り返し計算の収束時間を短くすることができる。

また、本発明のデータ推定方法は、前記各処理を繰り返す際に、前記予め設定された個数を増加させることにより、選択される前記モードの数を増加させるようにしてもよい。これにより、観測データから非線形入出力システムの入力データを推定する際に、繰り返し計算の収束時間を短くすることができる。

また、本発明のデータ推定方法は、前記入力候補ベクトルと前記補正ベクトルとに基づいて、新たな前記入力候補ベクトルを生成する際に、前記入力候補ベクトルと前記補正ベクトルとから求められる新たな前記入力候補ベクトルの各要素の値が、予め設定された範囲内であるか否かを判定し、新たな前記入力候補ベクトルの各要素に、前記予め設定された範囲内でない前記要素が存在する場合に、該要素に対応する、前記補正ベクトルのうちの要素の値を調整して前記補正ベクトルを修正し、前記入力候補ベクトルと修正された前記補正ベクトルとに基づいて、新たな前記入力候補ベクトルを生成するようにしてもよい。これにより、非線形入出力システムの入力ベクトルの各要素の値に制限がある場合に、非線形入出力システムの入力データを適切に推定することができる。

本発明によれば、逐次非線形性を有する非線形入出力システムを解析対象とする場合に、観測データから非線形入出力システムの入力データを推定することができる、という効果が得られる。

本実施形態に係るデータ推定装置の構成の一例を示すブロック図である。 建物と地盤とからなる非線形入出力システムＦを説明するための説明図である。 演算マトリクスの接線勾配を説明するための説明図である。 差分ベクトルを説明するための説明図である。 微小の摂動に応じて得られる出力候補変化ベクトルを説明するための説明図である。 モードの選択を説明するための説明図である。 第１実施形態に係るデータ推定処理ルーチンの一例を示す図である。 第２実施形態に係るデータ推定処理ルーチンの一例を示す図である。 シミュレーション実験の内容を説明するための説明図である。 シミュレーション実験の内容を説明するための説明図である。 シミュレーション実験の内容を説明するための説明図である。 シミュレーション実験の内容を説明するための説明図である。 シミュレーション実験の内容を説明するための説明図である。 シミュレーション実験の内容を説明するための説明図である。 シミュレーション実験の内容を説明するための説明図である。 シミュレーション実験の内容を説明するための説明図である。 本実施形態の変形例を説明するための説明図である。 本実施形態の変形例を説明するための説明図である。 本実施形態の変形例を説明するための説明図である。

以下、図面を参照して、本発明の実施形態について詳細に説明する。

＜本発明の第１実施形態に係るデータ推定装置の構成＞

図１に、本発明の第１実施形態に係るデータ推定装置１００の構成の一例を示す。データ推定装置１００は、機能的には、図１に示されるように、データ受付部１０、コンピュータ２０、及び出力部５０を含んだ構成で表すことができる。本実施形態のデータ推定装置１００は、入力ベクトルが入力されたときに出力ベクトルを出力する非線形入出力システムＦにおいて、任意の出力ベクトルが観測されたときの入力ベクトルを推定する。

非線形入出力システムＦの一例としては、非線形地盤モデルが挙げられる。例えば、非線形地盤モデルの任意の観測点において地震動が観測された場合を考える。この場合、観測点における地震動を表す出力ベクトルに基づいて、その出力ベクトルを再現するための入力地震動である入力ベクトルを推定する、という問題がある。このような入力ベクトルを推定する手法としては様々な手法が存在する。

例えば、地震発生後の建物の健全性評価において、対象建物の基礎構造等の評価を行う際には、地震動のシミュレーション解析が利用される。地震動のシミュレーション解析を用いる場合、入力地震動の波形を表す入力ベクトルを再現する必要がある。

図２に、建物と地盤とからなる非線形入出力システムＦを説明するための説明図を示す。図２に示されるように、非線形入出力システムＦは、建物本体（Building）、杭（Pile）、及び地盤（Soil）からなる。この場合、図２に示される観測点Ｓにおいて、地震動の観測データを表す出力波形（Output Wave）が得られる。このとき、得られた出力波形からインバージョンによって入力地震動を表す入力波形（Input Wave）が推定される。この場合、シミュレーション解析で用いるための入力地震動を表す入力波形は、非線形地盤モデル上の観測点Ｓにおいてその出力波形が再現されるように推定する必要がある。

入力地震動の設定方法としては、信号の周波数領域における方法が挙げられる。しかし、周波数領域による方法は解析中の地盤の物性値が一定である場合にのみ適用可能である。一般的には、地震時の地盤の挙動は逐次非線形性を示し、地盤の物性値は時々刻々変化する。このため、周波数領域による方法を汎用的に用いることは難しい。

そのため、地震動のシミュレーション解析に用いる入力地震動を設定する際には、時間領域による方法を採用する必要がある。時間領域で入力地震動を設定する方法に関する従来技術は、対象システムが線形系である場合が多い。非線形を考慮した事例としては、基盤入力地震動を未知数として入力地震動を逆算する手法が知られているが、入力波に多くの高振動成分が生じてしまうという課題がある。

このため、解析対象の挙動が逐次非線形性を示す場合に、観測データである出力ベクトルに基づいて入力データである入力ベクトルを適切に推定する汎用的な手法は知られていない。そこで、本実施形態のデータ推定装置１００は、逐次非線形性を有する非線形入出力システムを解析対象とする場合に、観測データである出力ベクトルから非線形入出力システムの入力データである入力ベクトルを適切に推定する。

出力ベクトル｛ａ_ｏ｝と非線形入出力システムＦと入力ベクトル｛ａ_ｉ｝との関係は、以下の式（１）によって表される。

（１）

上記図２に示されるインバージョン（Inversion）は、上記式（１）の逆演算、すなわち、出力ベクトル｛ａ_ｏ｝から入力ベクトル｛ａ_ｉ｝を計算することに相当する。この逆演算を実施する際には２つの課題が存在する。１つ目の課題は、上記式（１）の演算は順方向に単射変換であるが、逆方向には単射変換とはならない点である。２つ目の課題は、対象となるシステムＦが線形システムではなく、入力ベクトル｛ａ_ｉ｝に応じて出力ベクトル｛ａ_ｏ｝が複雑に変化する非線形システムである点である。

入力ベクトル｛ａ_ｉ｝がＮ次元であり、出力ベクトル｛ａ_ｏ｝がＮ次元である場合には、非線形入出力システムＦはＮｘＮの演算マトリクスとみなすことができる。そこで、本実施形態では、非線形入出力システムＦが表す演算マトリクスの接線勾配を用いて、出力ベクトル｛ａ_ｏ｝から入力ベクトル｛ａ_ｉ｝を推定する。

図３に、演算マトリクスの接線勾配を説明するための説明図を示す。図３に示されるように、入力ベクトル｛ａ_ｉ｝の変化分｛Δａ_ｉ｝（Incremental Input Wave Vector）と演算マトリクスの接線勾配を表す行列［Ｋ］とから、出力ベクトル｛ａ_ｏ｝の変化分｛Δａ_ｏ｝（Incremental Output Wave Vector）が演算される。本実施形態では、入力ベクトル｛ａ_ｉ｝に対する変化分｛Δａ_ｉ｝を与え、行列［Ｋ］を生成する。そのため、本実施形態では行列［Ｋ］を摂動インパルス応答マトリクスと称する。

また、本実施形態では、摂動インパルス応答マトリクス［Ｋ］を特異値分解し、特異値分解の結果から入力ベクトル｛ａ_ｉ｝の推定に有意なモードを選択し、推定された入力ベクトル｛ａ_ｉ｝の推定誤差を反復的に修正する。

以下、具体的に説明する。

データ受付部１０は、観測された任意の出力ベクトルを受け付ける。データ受付部１０は、例えばキーボード、マウス、又は外部装置からの入力を受け付ける入出力装置等によって実現される。出力ベクトルは、例えば、建物周辺の所定の各箇所に設置されたセンサによって得られる地震動の各時刻の観測データである。

コンピュータ２０は、ＣＰＵ（Central Processing Unit）、各処理ルーチンを実現するためのプログラム等を記憶したＲＯＭ（Read Only Memory）、データを一時的に記憶するＲＡＭ（Random Access Memory）、記憶手段としてのメモリ、ネットワークインタフェース等を含んで構成されている。コンピュータ２０は、図１に示されるように、機能的には、データ記憶部２１と、入力候補ベクトル計算部２２と、出力候補ベクトル計算部２４と、差分ベクトル計算部２６と、判定部２７と、入力候補変化ベクトル計算部２８と、出力候補変化ベクトル計算部３０と、応答行列計算部３２と、一般化逆行列計算部３４と、補正ベクトル計算部３６と、更新部３８と、結果取得部４０とを備えている。

データ記憶部２１には、データ受付部１０によって受け付けた出力ベクトルが格納される。なお、データ記憶部２１に格納された出力ベクトルは、目標ベクトルとして用いられる。

入力候補ベクトル計算部２２は、入力ベクトルの候補を表す入力候補ベクトル｛ａ_ｉ｝を設定する。例えば、入力候補ベクトル計算部２２は、以下の式（２）に示されるように、入力候補ベクトル｛ａ_ｉ｝に対して初期ベクトル｛ａ_ｄｕｍ｝を設定する。

（２）

観測データが地震動のデータである場合、初期ベクトル｛ａ_ｄｕｍ｝は、例えば、地盤の初期剛性を用いて振動数領域の方法によって求めた入力地震動を設定することができる。または、入力が完全に０の地震動であってもよい。なお、後述する繰り返し計算の繰り返し回数を減らすためには、精度の高い初期予測解を用いることが好ましい。

出力候補ベクトル計算部２４は、以下の式（３）に従って、入力候補ベクトル｛ａ_ｉ｝を非線形入出力システムＦへ入力したときに出力される出力ベクトルである出力候補ベクトル｛ａ_ｏ｝を計算する。

（３）

差分ベクトル計算部２６は、データ記憶部２１に格納された出力ベクトルを目標ベクトル｛ａ_ｔａｒ｝として設定する。そして、差分ベクトル計算部２６は、以下の式（４）に従って、目標ベクトル｛ａ_ｔａｒ｝と出力候補ベクトル計算部２４によって計算された出力候補ベクトル｛ａ_ｏ｝との間の差分を表す差分ベクトル｛ｒ｝を計算する。

（４）

図４に、差分ベクトルを説明するための説明図を示す。図４に示されるように、差分ベクトル｛ｒ｝（Residual wave vector）の各要素（各時刻の値）は、目標ベクトル｛ａ_ｔａｒ｝（Target wave vector）の各要素（各時刻の値）と出力候補ベクトル｛ａ_ｏ｝（Calculated wave vector By predicted input wave）の各要素（各時刻の値）との間の差分である。そのため、この差分ベクトルの各成分の値が十分に小さければ、推定された入力候補ベクトル｛ａ_ｉ｝は、十分な精度を確保できていることになる。

判定部２７は、差分ベクトル計算部２６によって計算された差分ベクトルに対応する値が予め定められた閾値ε未満であるか否かを判定する。本実施形態において、判定部２７は、差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜を抽出し、以下の式（５）に従って、最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満であるか否かを判定する。なお、閾値εは解析対象のモデルや用途、単位系に応じて適宜設定する（例えば、0.01［cm/s2］など）。

（５）

そして、最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満である場合には、後述する入力候補変化ベクトル計算部２８、出力候補変化ベクトル計算部３０、応答行列計算部３２、一般化逆行列計算部３４、補正ベクトル計算部３６、及び更新部３８による繰り返し計算が終了し、後述する結果取得部４０による処理が開始される。一方、最大値｜｛ｒ｝_ｍａｘ｜が閾値ε以上である場合には、後述する入力候補変化ベクトル計算部２８、出力候補変化ベクトル計算部３０、応答行列計算部３２、一般化逆行列計算部３４、補正ベクトル計算部３６、及び更新部３８による繰り返し計算が継続される。

本実施形態では、差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満であるか否かに応じて、繰り返し計算を終了させるか否かを判定する場合を例に説明するが、判定部２７はどのように判定してもよい。例えば、判定部２７は、差分ベクトル｛ｒ｝のノルムが所定の閾値未満である場合に、繰り返し計算の終了を判定するようにしてもよい。また、判定部２７は、差分ベクトル｛ｒ｝の各要素の平均が所定の閾値未満である場合に、繰り返し計算の終了を判定するようにしてもよい。

入力候補変化ベクトル計算部２８は、判定部２７によって差分ベクトルのうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε以上であると判定された場合に処理を実行する。判定部２７によって差分ベクトルのうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε以上であると判定された場合、入力候補変化ベクトル計算部２８は、入力候補ベクトル計算部２２により設定された入力候補ベクトル｛ａ_ｉ｝に含まれる複数の要素の各々について、当該要素を変化させたベクトルである入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）を計算する。

そして、出力候補変化ベクトル計算部３０は、入力候補変化ベクトル計算部２８によって計算された複数の入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）の各々について、入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）を非線形入出力システムＦへ入力したときの出力ベクトルである出力候補変化ベクトル｛ａ_ｏ ^ｊ｝を計算する。

具体的には、出力候補変化ベクトル計算部３０は、以下の式（６）及び式（７）に従って、入力候補ベクトル｛ａ_ｉ｝の各要素に微小の摂動｛Δａ_ｉ ^ｊ｝を加えた場合の応答解析を実施して、摂動｛Δａ_ｉ ^ｊ｝を加えた場合の出力候補変化ベクトル｛ａ_ｏ ^ｊ｝を得る。

（６）

（７）

ここで、微小の摂動量Δは、推定対象の入力ベクトルの推定精度に応じて設定される。例えば、入力ベクトルが入力地震動である場合、0.1［cm/s²］単位での精度が必要であれば、同程度の値を設定することが好ましい。なお、摂動量Δの値が大きい場合には、収斂の仮定で振動が生じ収斂解が得られなくなる可能性が高くなる。一方、摂動量Δの値が小さい場合には、収斂に多くの時間を要することになる。

図５に、微小の摂動｛Δａ_ｉ ^ｊ｝に応じて得られる出力候補変化ベクトル｛ａ_ｏ ^ｊ｝を説明するための説明図を示す。図５に示されるように、入力候補ベクトル｛ａ_ｉ｝のｊ番目の要素に摂動｛Δａ_ｉ ^ｊ｝を与えた応答解析を、全ての要素に対して実施することにより、Ｎ個の出力候補変化ベクトル｛ａ_ｏ ^ｊ｝が得られる。この出力候補変化ベクトル｛ａ_ｏ ^ｊ｝と上記式（３）において求めた出力候補ベクトル｛ａ_ｏ｝との間の差分は、入力候補ベクトル｛ａ_ｉ｝の各要素での摂動が出力に与える影響を示している。

応答行列計算部３２は、以下の式（８）に従って、出力候補変化ベクトル計算部３０によって計算された複数の出力候補変化ベクトル｛ａ_ｏ ^ｊ｝の各々について、当該出力候補変化ベクトル｛ａ_ｏ ^ｊ｝と当該出力候補ベクトル｛ａ_ｏ｝との間の差分から求まる変化率を表す影響係数ベクトル｛ｋ^ｊ｝を計算する。影響係数ベクトル｛ｋ^ｊ｝は、本発明の係数ベクトルの一例である。

（８）

上記図５に示されるように、影響係数ベクトル｛ｋ^ｊ｝の各要素は、入力候補ベクトル｛ａ_ｉ｝のｊ番目の要素に摂動｛Δａ_ｉ ^ｊ｝が与えられた場合の、出力候補変化ベクトル｛ａ_ｏ ^ｊ｝の各要素と出力候補ベクトル｛ａ_ｏ｝の各要素との間の差分から求まる変化率である。

そして、応答行列計算部３２は、複数の影響係数ベクトル｛ｋ^ｊ｝を列ベクトルとした行列を表す摂動インパルス応答マトリクス［Ｋ］を生成する。摂動インパルス応答マトリクス［Ｋ］は、本発明の応答行列の一例である。

具体的には、応答行列計算部３２は、計算された複数の影響係数ベクトル｛ｋ^ｊ｝を列ベクトルとして、以下の式（９）に示されるように、全ての影響係数ベクトル｛ｋ^ｊ｝を並べることにより、摂動インパルス応答マトリクス［Ｋ］（正方行列Ｎ×Ｎ）を生成する。

（９）

なお、摂動インパルス応答マトリクス［Ｋ］は、以下の式（１０）に示されるように、現時点で候補となっている入力候補ベクトル｛ａ_ｉ｝に対して変化分｛Δａ_ｉ｝を与えた場合の出力ベクトル｛ａ_ｏ｝の変化分｛Δａ_ｏ｝を算定するマトリクスとなる。

（１０）

したがって、摂動インパルス応答マトリクス［Ｋ］は、非線形入出力システムＦが表す演算マトリクスの接線勾配であり、上記図３に示されるように、入力ベクトルの微小変分から出力ベクトルの微小変分を算出するマトリクスである。

摂動インパルス応答マトリクス［Ｋ］の各列は、入力ベクトルのある要素に摂動を与えた場合の出力ベクトルのある要素の摂動となる。そのため、この摂動インパルス応答マトリクス［Ｋ］の上三角部分は０となり、場合によっては対角部分も０となる。本実施形態では、摂動インパルス応答マトリクス［Ｋ］の逆行列を用いて、差分を解消するような入力候補ベクトル｛ａ_ｉ｝に対する変化分｛Δａ_ｉ｝を逆算し、修正を行うことで徐々に収斂解を近似的に得る。

摂動インパルス応答マトリクス［Ｋ］を用いた場合であっても、「逆行列が一意に定まらない」及び「マトリクスが時々刻々変化する」という課題は残る。そこで、本実施形態では、前者については、通常の逆行列ではなく一般化逆行列を用い、かつその主要モードのみを用いて修正を行う。また、後者については、摂動インパルス応答マトリクス［Ｋ］を反復計算毎に求めることで解決する。本実施形態では、一般化逆行列を算出する際には、特異値分解を用いる。特異値分解は行列分解手法の一例である。

一般化逆行列計算部３４は、行列分解手法の一例である特異値分解を用いて、応答行列計算部３２によって生成された摂動インパルス応答マトリクス［Ｋ］を行列分解する。摂動インパルス応答マトリクス［Ｋ］の特異値分解が行われると、摂動インパルス応答マトリクス［Ｋ］は、以下の式（１１）に示されるように分解される。

（１１）

ここで、上記式（１１）における［Ｕ］はユニタリ行列であり、各列は出力の正規直交基底で構成される。また、上記式（１１）における［Ｖ］もユニタリ行列であり、その各列は入力の正規直交基底を示す。上記式（１１）における［Σ］は特異値を対角項にもつ特異値行列である。なお、この特異値は分解前の行列のランクの数だけ得られ、特異値の大きさは各モードの寄与度の大きさを表している。

次に、一般化逆行列計算部３４は、摂動インパルス応答マトリクス［Ｋ］の行列分解結果から得られる複数のモードから特定のモードを選択する。具体的には、一般化逆行列計算部３４は、摂動インパルス応答マトリクス［Ｋ］の行列分解結果から得られる複数のモードから特定のモードを選択する際に、モードの寄与度を表す特異値が予め設定されたモードに関する閾値ＴＨ以上であるモードを選択する。したがって、一般化逆行列計算部３４は、特異値が０又は比較的小さなモードは除外し、Ｎ’次のモードまでを考慮するものとして、各部分行列を縮退させる。なお、部分行列は縮退するが全体としての影響係数マトリクスの次元は変わらない。一般化逆行列計算部３４によるモードの選択によって得られる摂動インパルス応答マトリクス［Ｋ’］は、以下の式（１２）によって表される。

（１２）

図６に、モードの選択を説明するための説明図を示す。一般化逆行列計算部３４は、図６に示されるように、特異値行列［Σ］から、特異値がモードに関する閾値ＴＨ以上であるモードを選択して、新たな特異値行列［Σ’］を生成する。

このように、特異値の小さなモードを無視するのは、特異値が０又は小さなモードを考慮した状態で逆行列を算出すると、同定された入力ベクトルに特異値（すなわち、影響度）の小さいモードの成分が多く含まれてしまうためである。また、数値計算上の０に近い値での除算等が生じ、多くの数値誤差が入ることを避ける意味合いもある。

本実施形態のデータ推定手法は、繰り返しによって誤差を小さくしていく近似解法であるため、微小なモードは無視し、有意なモードを選択して修正を行う。なお、選別する特異値の境界値は、対象とする問題に応じて適宜設定する必要がある。なお、微小なモードによる修正だけでは、要求される許容誤差以内に収まらない可能性がある。そのような場合には、繰り返し回数に応じて徐々に閾値を下げ、より多くのモードを用いて修正を行うことが必要になる。

次に、一般化逆行列計算部３４は、選択されたモードに対応する部分行列を用いて、摂動インパルス応答マトリクス［Ｋ’］の一般化逆行列［Ｋ’］^＋を生成する。

具体的には、一般化逆行列計算部３４は、上記式（１２）から、各行列のユニタリ性及び対角性を利用して、摂動インパルス応答マトリクス［Ｋ’］の一般化逆行列［Ｋ’］^＋を以下の式（１３）によって求める。

（１３）

補正ベクトル計算部３６は、一般化逆行列計算部３４によって生成された一般化逆行列［Ｋ’］^＋と、差分ベクトル計算部２６によって計算された差分ベクトル｛ｒ｝とに基づいて、以下の式（１４）に従って、入力候補ベクトル｛ａ_ｉ｝に対する補正量を表す補正ベクトル｛Δａ_ｉ｝を計算する。

（１４）

更新部３８は、入力候補ベクトル計算部２２により設定された入力候補ベクトル｛ａ_ｉ｝と補正ベクトル計算部３６により計算された補正ベクトル｛Δａ_ｉ｝とに基づいて、新たな入力候補ベクトル｛ａ_ｉ｝を生成する。

具体的には、更新部３８は、以下の式（１５）に示されるように、補正ベクトル｛Δａ_ｉ｝を入力候補ベクトル計算部２２により設定された入力候補ベクトル｛ａ_ｉ｝に足しこみ、新たな入力候補ベクトル｛ａ_ｉ｝とする。

（１５）

そして、差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満となるまで、出力候補ベクトル｛ａ_ｏ｝の計算、差分ベクトル｛ｒ｝の計算、入力候補変化（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）ベクトルの計算、出力候補変化ベクトル｛ａ_ｏ ^ｊ｝の計算、摂動インパルス応答マトリクス［Ｋ］の生成、一般化逆行列［Ｋ’］^＋の生成、及び新たな入力候補ベクトル｛ａ_ｉ｝の生成の各処理が繰り返される。

結果取得部４０は、差分ベクトル計算部２６によって計算される差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満となった場合に、更新部３８によって生成された入力候補ベクトル｛ａ_ｉ｝を、データ記憶部２１に格納された目標ベクトルに対応する入力ベクトルとして取得する。

出力部５０は、結果取得部４０によって取得された入力ベクトルを結果として出力する。例えば、出力部５０は、ディスプレイによって実現される。

＜データ推定装置１００の作用＞

次に、データ推定装置１００の作用を説明する。データ推定装置１００のデータ受付部１０が、観測された任意の出力ベクトルの入力を受け付けると、データ記憶部２１へ格納する。そして、データ推定装置１００のコンピュータ２０は、処理実行の指示信号を受け付けると、図７Ａに示すデータ推定処理ルーチンを実行する。なお、以下の処理では、摂動インパルス応答マトリクス［Ｋ］に対する一般化逆行列［Ｋ］^＋を生成した後に、特異値行列［Σ］から特定のモードを選択する場合を例に説明する。

ステップＳ１００において、入力候補ベクトル計算部２２は、上記式（２）に従って、入力ベクトルの候補を表す入力候補ベクトル｛ａ_ｉ｝を設定する。

ステップＳ１０２において、出力候補ベクトル計算部２４は、上記式（３）に従って、上記ステップＳ１００で設定された入力候補ベクトル｛ａ_ｉ｝を非線形入出力システムＦへ入力したときに出力される出力ベクトルである出力候補ベクトル｛ａ_ｏ｝を計算する。

ステップＳ１０４において、差分ベクトル計算部２６は、データ記憶部２１に格納された出力ベクトルを目標ベクトル｛ａ_ｔａｒ｝として設定する。そして、差分ベクトル計算部２６は、上記式（４）に従って、目標ベクトル｛ａ_ｔａｒ｝と上記ステップＳ１０２で計算された出力候補ベクトル｛ａ_ｏ｝との間の差分を表す差分ベクトル｛ｒ｝を計算する。

ステップＳ１０６において、判定部２７は、上記ステップＳ１０４で計算された差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜を抽出し、上記式（５）に従って、最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満であるか否かを判定する。そして、最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満である場合には、データ推定処理ルーチンを終了する。一方、最大値｜｛ｒ｝_ｍａｘ｜が閾値ε以上である場合には、ステップＳ１０８へ進む。

ステップＳ１０８において、入力候補変化ベクトル計算部２８は、上記ステップＳ１００で設定された入力候補ベクトル｛ａ_ｉ｝に含まれる複数の要素の各々について、上記式（７）に従って、当該要素を変化させたベクトルである入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）を計算する。

そして、ステップＳ１０８において、出力候補変化ベクトル計算部３０は、入力候補変化ベクトル計算部２８によって計算された複数の入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）の各々について、上記式（６）に従って、入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）を非線形入出力システムＦへ入力したときの出力候補変化ベクトル｛ａ_ｏ ^ｊ｝を計算する。

ステップＳ１１０において、応答行列計算部３２は、上記式（８）に従って、上記ステップＳ１０８で計算された複数の出力候補変化ベクトル｛ａ_ｏ ^ｊ｝の各々について、当該出力候補変化ベクトル｛ａ_ｏ ^ｊ｝と上記ステップＳ１０２で計算された出力候補ベクトル｛ａ_ｏ｝との間の差分から求まる影響係数ベクトル｛ｋ^ｊ｝を計算する。

ステップＳ１１２において、応答行列計算部３２は、上記式（９）に示されるように、上記ステップＳ１１０で計算された複数の影響係数ベクトル｛ｋ^ｊ｝を列ベクトルとして、摂動インパルス応答マトリクス［Ｋ］を生成する。

ステップＳ１１４において、一般化逆行列計算部３４は、上記ステップＳ１１２で生成された摂動インパルス応答マトリクス［Ｋ］を特異値分解し、特異値分解結果に応じて摂動インパルス応答マトリクス［Ｋ］の一般化逆行列［Ｋ］^＋を生成する。

ステップＳ１１６において、一般化逆行列計算部３４は、一般化逆行列［Ｋ］^＋を構成する特異値行列［Σ］から、特異値がモードに関する閾値ＴＨ以上であるモードを選択して、新たな特異値行列［Σ’］を生成する。また、一般化逆行列計算部３４は、選択された特異値に応じて、新たな行列［Ｕ’］と新たな行列［Ｖ’］とを生成する。そして、一般化逆行列計算部３４は、新たな特異値行列［Σ’］と新たな行列［Ｕ’］と新たな行列［Ｖ’］とに基づいて、新たな一般化逆行列［Ｋ’］^＋を生成する。

ステップＳ１１８において、補正ベクトル計算部３６は、上記ステップＳ１１６で生成された一般化逆行列［Ｋ’］^＋と、上記ステップＳ１０４で計算された差分ベクトル｛ｒ｝とに基づいて、上記式（１４）に従って、入力候補ベクトル｛ａ_ｉ｝に対する補正ベクトル｛Δａ_ｉ｝を計算する。

ステップＳ１２０において、更新部３８は、上記ステップＳ１００で設定された入力候補ベクトル｛ａ_ｉ｝又は前回のステップＳ１２０で更新された入力候補ベクトル｛ａ_ｉ｝と、上記ステップＳ１１８で計算された補正ベクトル｛Δａ_ｉ｝とに基づいて、上記式（１５）に従って、新たな入力候補ベクトル｛ａ_ｉ｝を生成する。

データ推定処理ルーチンでは、差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満となるまで、出力候補ベクトル｛ａ_ｏ｝の計算、差分ベクトル｛ｒ｝の計算、入力候補変化（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）ベクトルの計算、出力候補変化ベクトル｛ａ_ｏ ^ｊ｝の計算、摂動インパルス応答マトリクス［Ｋ］の生成、一般化逆行列［Ｋ’］^＋の生成、及び新たな入力候補ベクトル｛ａ_ｉ｝の生成の各処理が繰り返される。

結果取得部４０は、差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満となった場合に、更新部３８によって生成された入力候補ベクトル｛ａ_ｉ｝を、データ記憶部２１に格納された目標ベクトルに対応する入力ベクトルとして取得する。

出力部５０は、結果取得部４０によって取得された入力ベクトルを結果として出力する。

以上詳細に説明したように、本実施形態のデータ推定装置１００は、入力候補ベクトル｛ａ_ｉ｝を非線形入出力システムＦへ入力したときに出力される出力候補ベクトル｛ａ_ｏ｝を計算し、観測された任意の出力ベクトルを目標ベクトル｛ａ_ｔａｒ｝として、目標ベクトル｛ａ_ｔａｒ｝と出力候補ベクトル｛ａ_ｏ｝との間の差分を表す差分ベクトル｛ｒ｝を計算する。そして、データ推定装置１００は、入力候補ベクトル｛ａ_ｉ｝に含まれる複数の要素の各々について、当該要素を変化させたベクトルである入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）を計算し、複数の入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）の各々について、入力候補変化ベクトル（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）を非線形入出力システムＦへ入力したときの出力候補変化ベクトル｛ａ_ｏ ^ｊ｝を計算する。そして、データ推定装置１００は、複数の出力候補変化ベクトル｛ａ_ｏ ^ｊ｝の各々について、影響係数ベクトル｛ｋ｝を計算し、摂動インパルス応答マトリクス［Ｋ］を生成する。そして、データ推定装置１００は、摂動インパルス応答マトリクス［Ｋ］の行列分解結果から得られる複数のモードから特定のモードを選択し、選択されたモードに対応する部分行列を用いて、摂動インパルス応答マトリクス［Ｋ］の一般化逆行列［Ｋ’］^＋を生成する。そして、データ推定装置１００は、一般化逆行列［Ｋ’］^＋と差分ベクトル｛ｒ｝とに基づいて、入力候補ベクトル｛ａ_ｉ｝に対する補正ベクトル｛Δａ_ｉ｝を計算する。そして、データ推定装置１００は、入力候補ベクトル｛ａ_ｉ｝と補正ベクトル｛Δａ_ｉ｝とに基づいて、新たな入力候補ベクトル｛ａ_ｉ｝を生成する。また、データ推定装置１００は、差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満となるまで、出力候補ベクトル｛ａ_ｏ｝の計算、差分ベクトル｛ｒ｝の計算、入力候補変化（｛ａ_ｉ｝+｛Δａ_ｉ ^ｊ｝）ベクトルの計算、出力候補変化ベクトル｛ａ_ｏ ^ｊ｝の計算、摂動インパルス応答マトリクス［Ｋ］の生成、一般化逆行列［Ｋ’］^＋の生成、及び新たな前記入力候補ベクトル｛ａ_ｉ｝の生成の各処理を繰り返す。そして、データ推定装置１００は、差分ベクトル｛ｒ｝の各要素のうちの最大値｜｛ｒ｝_ｍａｘ｜が閾値ε未満となった場合に、入力候補ベクトル｛ａ_ｉ｝を、観測された任意の出力ベクトルに対応する入力ベクトルとして取得する。これにより、逐次非線形性を有する非線形入出力システムを解析対象とする場合に、観測データから非線形入出力システムの入力データを推定することができる。

＜本発明の第２実施形態に係るデータ推定装置の構成＞

次に、本発明の第２実施形態について説明する。第２実施形態のデータ推定装置は、入力候補ベクトル｛ａ_ｉ｝に対する補正ベクトル｛Δａ_ｉ｝を計算する際に、入力候補ベクトルと補正ベクトルとから求められる新たな入力候補ベクトルの各要素の値が、予め設定された範囲内であるか否かを判定する。そして、第２実施形態のデータ推定装置は、新たな入力候補ベクトルの各要素に、予め設定された範囲内でない要素が存在する場合には、補正ベクトルのうちの要素の値を調整して補正ベクトルを修正する。なお、第２の実施形態に係るデータ推定装置の構成は、第１実施形態と同様の構成となるため、同一符号を付して説明を省略する。

第２実施形態の補正ベクトル計算部３６は、上記式（１５）に従って、上記式（１４）に従って計算された補正ベクトル｛Δａ_ｉ｝と、入力候補ベクトル計算部２２により設定された入力候補ベクトル｛ａ_ｉ｝とに基づき、新たな入力候補ベクトル｛ａ_ｉ｝を一旦計算する。

次に、補正ベクトル計算部３６は、入力候補ベクトル｛ａ_ｉ｝と補正ベクトル｛Δａ_ｉ｝とから一旦求められた新たな入力候補ベクトル｛ａ_ｉ｝の各要素の値が、予め設定された範囲内であるか否かを判定する。

例えば、本実施形態のデータ推定装置を用いて同定を行う場合、入力候補ベクトル｛ａ_ｉ｝の各要素に、負になってはいけないパラメータ（例えば、ヤング率又は耐力等）が含まれている場合、各要素の値として負の値が候補として探索されてしまい、計算の収斂が停止してしまう場合がある。

そこで、第２実施形態のデータ推定装置は、入力候補ベクトル｛ａ_ｉ｝の各要素が予め設定された範囲外とならないように、補正ベクトル｛Δａ_ｉ｝を修正する。これにより、非線形入出力システムの入力ベクトルの各要素の値に制限がある場合に、非線形入出力システムの入力データを適切に推定することができる。

なお、入力候補ベクトル｛ａ_ｉ｝の各要素のうち、制約がない要素に関しては予め定められた範囲を設定する必要はない。

そして、補正ベクトル計算部３６は、一旦計算された新たな入力候補ベクトル｛ａ_ｉ｝の各要素に、予め設定された範囲内でない要素が存在する場合に、以下の式（１６）に従って、補正ベクトル｛Δａ_ｉ｝のうちの各要素の値を調整して補正ベクトル｛Δａ_ｉ｝を修正する。

（１６）

第２実施形態の更新部３８は、補正ベクトル計算部３６による修正が行われた後、上記式（１５）に示されるように、補正ベクトル｛Δａ_ｉ｝を入力候補ベクトル計算部２２により設定された入力候補ベクトル｛ａ_ｉ｝に足しこみ、新たな入力候補ベクトル｛ａ_ｉ｝とする。

＜データ推定装置の作用＞

次に、第２実施形態のデータ推定装置の作用を説明する。第２実施形態のデータ推定装置のデータ受付部１０が、観測された任意の出力ベクトルの入力を受け付けると、データ記憶部２１へ格納する。そして、第２実施形態のデータ推定装置のコンピュータ２０は、処理実行の指示信号を受け付けると、図７Ｂに示すデータ推定処理ルーチンを実行する。

ステップＳ２１９において、第２実施形態の補正ベクトル計算部３６は、上記式（１５）に従って、ステップＳ１００で設定された入力候補ベクトル｛ａ_ｉ｝又は前回のステップＳ１２０で更新された入力候補ベクトル｛ａ_ｉ｝と、ステップＳ１１８で計算された補正ベクトル｛Δａ_ｉ｝とに基づいて、新たな入力候補ベクトル｛ａ_ｉ｝を一旦計算する。

そして、ステップＳ２１９において、補正ベクトル計算部３６は、一旦求めた新たな入力候補ベクトル｛ａ_ｉ｝の各要素の値が、予め設定された範囲内であるか否かを判定する。これにより、補正ベクトル｛ａ_ｉ｝による入力候補ベクトル｛ａ_ｉ｝の変動が許容可能であるか否かが判定される。そして、一旦求めた新たな入力候補ベクトル｛ａ_ｉ｝の各要素の値が、予め設定された範囲外である場合には、ステップＳ２２０へ進む。一方、一旦求めた新たな入力候補ベクトル｛ａ_ｉ｝の各要素の値が、予め設定された範囲内である場合にはステップＳ１２０へ進む。

ステップＳ２２０において、補正ベクトル計算部３６は、上記式（１６）に従って、補正ベクトル｛Δａ_ｉ｝のうちの各要素の値を調整して補正ベクトル｛Δａ_ｉ｝を修正する。

ステップＳ１２０において、更新部３８は、上記ステップＳ１００で設定された入力候補ベクトル｛ａ_ｉ｝又は前回のステップＳ１２０で更新された入力候補ベクトル｛ａ_ｉ｝と、上記ステップＳ１１８で計算された補正ベクトル｛Δａ_ｉ｝又は上記ステップＳ２２０で修正された補正ベクトル｛Δａ_ｉ｝とに基づいて、上記式（１５）に従って、新たな入力候補ベクトル｛ａ_ｉ｝を生成する。

以上詳細に説明したように、第２実施形態のデータ推定装置は、入力候補ベクトル｛ａ_ｉ｝に対する補正ベクトル｛Δａ_ｉ｝を計算する際に、入力候補ベクトルと補正ベクトルとから求められる新たな入力候補ベクトルの各要素の値が、予め設定された範囲内であるか否かを判定する。そして、第２実施形態のデータ推定装置は、新たな入力候補ベクトルの各要素に、予め設定された範囲内でない要素が存在する場合に、補正ベクトルのうちの要素の値を調整して補正ベクトルを修正する。これにより、非線形入出力システムの入力ベクトルの各要素の値に制限がある場合に、非線形入出力システムの入力データを適切に推定することができる。

＜シミュレーション実験＞

次に、例題を用いたシミュレーション実験について説明する。

図８に、対象の解析モデルを示す。解析モデルは深さ50mの柱状の自由地盤モデルとする。また、地盤は軟弱な砂地盤を想定し、Vs値は100[m/s]、ポアソン比は0.33、単位体積重量は1.8[t/m³]とする。地盤の非線形特性にはRamberg-Osgood（RO）モデル（例えば、参考文献（Jennnings, P. C., "Periodic Response of General Yielding Structure", J. Eng. Mech. Div., ASCE, EM2, pp.131-163, 1964）を参照）を採用する。ROモデルのパラメータは、参考文献（Architectural Institute of Japan: "Seismic Response Analysis and Design of Buildings Considering Dynamic Soil-Structure Interaction", Maruzen, 2006(in Japanese)）に記載された砂質土の値を参照し、せん断剛性低下率（G/G0）が50%となるときのせん断ひずみγ_0.5=0.10%、最大減衰定数hmax=0.21となるように設定する。

また、地盤下部は、粘性境界とし、Vs値は300m/s、単位体積重量は1.8t/m3に相当する粘性ダンパーを設定している。時間積分にはNewmark-β法（β=1/4）を用い、解析時間刻みは0.01秒とする。解析プログラムにはMuDIAN（例えば、参考文献（塩見忠彦ほか3名，「地盤の非線形挙動を考慮した構造解析システムの開発」，竹中技術研究所報告，No.54，pp.1-8，1998）を参照）を用いる。

シミュレーションによって、この解析モデルに対して、図９（ａ）（ｂ）に示すTaft1952EW地震動（最大加速度1.76m/s²、継続時間10秒まで）をα倍（α＝1,2,3）したものを入力地震動（入力ベクトル）とする応答解析を実施し、地表面出力（出力ベクトル）を算出する。解析時間については、最も大きな３倍入力の場合で約３秒であった。

図１０に、算定された地盤の最大応答加速度（Max Acc）、最大変位（Max Disp）、最大せん断ひずみ（Max Shear Strain）の分布図を示す。図１０より、地盤の最大せん断ひずみは1倍入力の場合でも0.1%以上に達しており非線形化が生じていることがわかる。また、3倍入力の場合においては、さらにひずみが増大し、-30m以深においては1%以上のせん断ひずみが生じていることから、強非線形領域に達していることが予測される。このことは、加速度分布や変位分布が、単純な係数倍の応答となっていないことからも判断できる。

図１１に、入力地震動と地表面応答の加速度時刻歴の伝達関数（Transfer function）（ParzenWindow0.1Hz）を示す。図１１より、入力地震動が大きくなるにつれて、非線形化の影響により全周期帯に渡って伝達関数が変化していることが確認できる。

これらの解析から得られた地表面の加速度時刻歴を再現するような入力地震動が提案手法によって算定できるか、また算定された入力地震動は、この解析に用いた地震動とどの程度一致するのかを確認する。

なお、上記式（５）で設定する許容誤差は0.01m/s²、上記式（８）で設定する摂動量についても同じ0.01m/s²とし、繰り返し回数は100回で打ち切りとした。また、上記式（１２）で設定するモードの数については、初期のイテレーションでは「特異値1以上」のモードを抽出することとし、以降、上記式（１４）の加速度修正ベクトルの成分の最大値が0.5m/s²未満になった時点で、特異値の閾値を半分にしていく。これにより徐々に高次のモードを含んだ修正を行うようにしている。

図１２に、ターゲットとする地表面出力から同定された入力地震動（Calculated）を実線で示し、シミュレーション解析で用いた入力地震動（Target）を点線で示す。また、図１３に、図１２に示された入力地震動によって得られる地表面出力の加速度時刻歴（Calculated）を実線で示し、シミュレーション解析で用いた入力地震動から得られる地表面出力（Target）を点線で示す。なお、これらの同定結果はターゲットとする地表面出力とのRMS（Root Mean Square）誤差が概ね1.0を下回る場合のケースを抽出している。

図１２及び図１３より、いずれの入力倍率においても、同定された入力地震動に若干の差はあるものの、地表面出力はよい精度で一致していることがわかる。このように地表面出力は概ね整合している一方で、同定された入力地震動の差が大きい要因は、前述したようにこの逆問題が単射問題でないためである。すなわち、ターゲットとする地表面出力を実現する入力地震動は複数存在し、一意には定まらないため、今回同定された入力地震動も解の一つではあるものの、前節で入力地震動として用いた地震動と完全には一致しなかったものと考えられる。ただし、いずれの入力倍率においても、初期の入力が小さく略線形とみなせる領域においては、入力地震動の波形はよい精度で一致している。これは、線形状態での入出力関係は１対１対応するため、精度良く算定できたものと考えられる。

以下のTable１に、RMS誤差を1.0未満とするために必要な収斂回数および計算時間を示す。また、図１４には収斂によるRMS誤差の減少の様子を示す。なお、Table1の数値は、図１４における「pick up point」における値である。非線形性の小さな1倍入力の場合は収斂が早いが、非線形化が進む2倍入力、3倍入力では多くの収斂回数、収斂時間がかかっている。非線形化が進んだ解析においては応答が複雑になり、誤差が収束しにくくなったものと考えられる。

前述したように、本実施形態における提案手法の特徴は、反復的な修正を有意なモードに限定して行うことにある。しかしながら、考え方として全てのモードを用いる方法も当然考えられる。

そこで、図１５を用いて、入力1倍のケースについて、考慮するモードの数を変化させた場合の収束性の違いを確認する。比較するケースは、前述した特異値の最小値を徐々に漸減させる場合（Gradient Decrease）に加えて、1.0（const），10^-6で固定（Min. Singular value）とした場合を追記している。前者は、常に低次のモードのみで修正を行った場合、後者はほぼすべてのモードを用いて修正を行った場合に相当する。

図１５より、全モードを用いた10^-6で固定した場合においては、収斂せず逆に発散してしまっている。一方、最小値を1.0で固定した低次モードのみでの修正でも、誤差の収斂は遅く収斂回数を重ねても誤差は小さくならない。したがって、初めは低次モードで修正し、誤差が小さくならなくなる場合には徐々に高次のモードを追加して修正していくという手法が有効といえる。

なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

上記実施形態のデータ推定方法は、出力ベクトル｛ａ_ｏ｝と非線形入出力システムＦと入力ベクトル｛ａ_ｉ｝とが、どのようなものであっても適用することができる。

例えば、図１６に示されるように、非線形入出力システムは、ニューロンの数が予め設定されたニューラルネットワークＦであり、入力ベクトル｛ａ_ｉ｝はニューラルネットワークのパラメータ｛Ｗ１，Ｗ２，・・・，Ｗ１００｝であり、出力ベクトル｛ａ_ｏ｝はニューラルネットワークの出力であるように設定することもできる。この場合、差分ベクトル｛ｒ｝は、学習用の入力データ（TestData）がニューラルネットワークＦに入力された際の出力ベクトルと、学習用の入力データ（TestData）に対する正解の出力ベクトルを表す目標ベクトル（Target Output｛Test1,Test2,Test3,…｝）との間の差分であるように設定される。

また、上記実施形態では、摂動インパルス応答マトリクス［Ｋ］の行列分解結果から得られる複数のモードから特定のモードを選択する際に、モードの寄与度を表す値がモードに関する閾値ｘ以上であるモードを選択し、選択されたモードに対応する部分行列を用いて、摂動インパルス応答マトリクス［Ｋ］の一般化逆行列［Ｋ’］^＋を生成する場合を例に説明したが、これに限定されるものではない。例えば、摂動インパルス応答マトリクス［Ｋ］の行列分解結果から得られる複数のモードから特定のモードを選択する際に、モードの寄与度を表す値が大きい方から予め設定された個数のモードを選択し、選択されたモードに対応する部分行列を用いて、摂動インパルス応答マトリクス［Ｋ］の一般化逆行列［Ｋ’］^＋を生成するようにしてもよい。

また、図１７に示されるように、非線形入出力システムＦは、入力ベクトルが入力されたときに各々異なる観測点での出力ベクトルを出力するものであるようにしてもよい。この場合、データ推定装置１００は、各々異なる観測点で観測された任意の複数の出力ベクトルを、複数の目標ベクトル（Target Output1，Target Output2）として、各々異なる観測点で任意の複数の出力ベクトルが観測されたときの入力ベクトル（Input）を推定する。この場合、入力ベクトル｛ａ_ｉ｝の変化分｛Δａ_ｉ｝と摂動インパルス応答マトリクス［Ｋ］と出力ベクトル｛ａ_ｏ｝の変化分｛Δａ_ｏ｝との関係は、図１７に示されるようになり、出力ベクトル｛ａ_ｏ｝の変化分｛Δａ_ｏ｝には、複数の異なる観測点での出力ベクトル｛ａ_ｏ｝の変化分｛Δａ_ｏ｝が含まれる。

また、入力ベクトルは、複数種類のベクトルを含むようにしてもよい。例えば、図１８に示されるように、入力地震動（Input）と地盤剛性K0とを推定したい場合には、入力地震動を表す加速度データAccと地盤剛性K0を入力ベクトルとして設定することができる。この場合、入力ベクトル｛ａ_ｉ｝の変化分｛Δａ_ｉ｝と摂動インパルス応答マトリクス［Ｋ］と出力ベクトル｛ａ_ｏ｝の変化分｛Δａ_ｏ｝との関係は、図１８に示されるようになり、入力ベクトル｛ａ_ｉ｝の変化分｛Δａ_ｉ｝には、地盤剛性K0に関する摂動が含まれる。

また、各処理を繰り返す際に、選択されるモードの数を増加させるようにしてもよい。この場合には、例えば、各処理が繰り返される毎に、モードに関する閾値ｘを減少させることにより、選択されるモードの数を増加させることができる。または、各処理が繰り返される毎に、選択されるモードの個数を増加させることにより、選択されるモードの数を増加させることができる。初めは低次モードで修正し、徐々に高次のモードを追加して修正していくという手法は、図１５に示されるように有効といえる。この場合には、図１５に示されるように、観測データから非線形入出力システムの入力データを推定する際に、繰り返し計算の収束時間を短くすることができる。

また、上記ではプログラムが記憶部（図示省略）に予め記憶（インストール）されている態様を説明したが、プログラムは、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ及びマイクロＳＤカード等の記録媒体の何れかに記録されている形態で提供することも可能である。

１０ データ受付部
２０ コンピュータ
２１ データ記憶部
２２ 入力候補ベクトル計算部
２４ 出力候補ベクトル計算部
２６ 差分ベクトル計算部
２７ 判定部
２８ 入力候補変化ベクトル計算部
３０ 出力候補変化ベクトル計算部
３２ 応答行列計算部
３４ 一般化逆行列計算部
３６ 補正ベクトル計算部
３８ 更新部
４０ 結果取得部
１００ データ推定装置

Claims (7)

1. 入力ベクトルが入力されたときに出力ベクトルを出力する非線形入出力システムにおいて、任意の前記出力ベクトルが観測されたときの前記入力ベクトルを推定するデータ推定方法であって、
   前記入力ベクトルの候補を表す入力候補ベクトルを、前記非線形入出力システムへ入力したときに出力される前記出力ベクトルを表す出力候補ベクトルを計算し、
   観測された任意の前記出力ベクトルを目標ベクトルとして、前記目標ベクトルと前記出力候補ベクトルとの間の差分を表す差分ベクトルを計算し、
   前記入力候補ベクトルに含まれる複数の要素の各々について、前記要素を変化させたベクトルである入力候補変化ベクトルを計算し、
   複数の前記入力候補変化ベクトルの各々について、前記入力候補変化ベクトルを前記非線形入出力システムへ入力したときの前記出力ベクトルを表す出力候補変化ベクトルを計算し、
   複数の前記出力候補変化ベクトルの各々について、前記出力候補変化ベクトルと前記出力候補ベクトルとの間の差分から求まる変化率を表す係数ベクトルを計算し、複数の係数ベクトルを列ベクトルとした行列を表す応答行列を生成し、
   前記応答行列の行列分解結果から得られる複数のモードから特定のモードを選択し、選択された前記モードに対応する部分行列を用いて、前記応答行列の一般化逆行列を生成し、
   前記一般化逆行列と前記差分ベクトルとに基づいて、前記入力候補ベクトルに対する補正量を表す補正ベクトルを計算し、
   前記入力候補ベクトルと前記補正ベクトルとに基づいて、新たな前記入力候補ベクトルを生成し、
   前記差分ベクトルに対応する値が予め定められた閾値未満となるまで、前記出力候補ベクトルの計算、前記差分ベクトルの計算、前記入力候補変化ベクトルの計算、前記出力候補変化ベクトルの計算、前記応答行列の生成、前記一般化逆行列の生成、及び前記新たな前記入力候補ベクトルの生成の各処理を繰り返し、
   前記差分ベクトルに対応する値が前記閾値未満となった場合に、前記入力候補ベクトルを、観測された任意の前記出力ベクトルに対応する前記入力ベクトルとして取得する、
   処理をコンピュータに実行させるデータ推定方法。
2. 前記応答行列の行列分解結果から得られる複数のモードから特定のモードを選択する際に、前記モードの寄与度を表す値が予め設定されたモードに関する閾値以上である前記モードを選択し、選択された前記モードに対応する部分行列を用いて、前記応答行列の一般化逆行列を生成する、
   請求項１に記載のデータ推定方法。
3. 前記応答行列の行列分解結果から得られる複数のモードから特定のモードを選択する際に、前記モードの寄与度を表す値が大きい方から予め設定された個数の前記モードを選択し、選択された前記モードに対応する部分行列を用いて、前記応答行列の一般化逆行列を生成する、
   請求項１に記載のデータ推定方法。
4. 前記非線形入出力システムは、入力ベクトルが入力されたときに各々異なる観測点での出力ベクトルを出力するものであって、
   各々異なる観測点で観測された任意の複数の前記出力ベクトルを、複数の前記目標ベクトルとして、各々異なる観測点で任意の複数の前記出力ベクトルが観測されたときの前記入力ベクトルを推定する、
   請求項１〜請求項３の何れか１項に記載のデータ推定方法。
5. 前記入力ベクトルは、複数種類のベクトルを含む、
   請求項１〜請求項３の何れか１項に記載のデータ推定方法。
6. 前記入力候補ベクトルと前記補正ベクトルとに基づいて、新たな前記入力候補ベクトルを生成する際に、
   前記入力候補ベクトルと前記補正ベクトルとから求められる新たな前記入力候補ベクトルの各要素の値が、予め設定された範囲内であるか否かを判定し、
   新たな前記入力候補ベクトルの各要素に、前記予め設定された範囲内でない前記要素が存在する場合に、該要素に対応する、前記補正ベクトルのうちの要素の値を調整して前記補正ベクトルを修正し、
   前記入力候補ベクトルと修正された前記補正ベクトルとに基づいて、新たな前記入力候補ベクトルを生成する、
   請求項１〜請求項５の何れか１項に記載のデータ推定方法。
7. 前記非線形入出力システムは、ニューロンの数が予め設定されたニューラルネットワークであり、
   前記入力ベクトルは、前記ニューラルネットワークのパラメータであり、
   前記出力ベクトルは、前記ニューラルネットワークの出力であり、
   前記差分ベクトルは、学習用の入力データが前記ニューラルネットワークに入力された際の前記出力ベクトルと、前記学習用の入力データに対する正解の前記出力ベクトルを表す前記目標ベクトルとの間の差分である、
   請求項１〜請求項３の何れか１項に記載のデータ推定方法。

Images