(57)【特許請求の範囲】
【請求項１】
非線形な物理モデルを用いてシミュレーションを行い、その出力と最適化対象の出力を比較し、その比較結果に基づき物理モデルのパラメータを繰り返し修正することにより、各パラメータを真値に漸近させていくパラメータの同定装置であって、
パラメータを変化させた時の出力評価関数値の変化量について、比較して強い相関の組み合わせがあるかを全てのパラメータ同士の出力評価関数値の変化量について判定し、
強く相関の組み合わせがある場合に物理モデルまたは最適化対象の出力が物理パラメータを真値に漸近させるための十分な情報を含んでいないと、判断することを特徴とするパラメータの同定装置。
【請求項２】
請求項１に記載のパラメータの同定装置において、
前記出力評価関数の変化量は、前記物理モデルの出力と最適化対象の出力の誤差についてのパラメータ変化の影響であることを特徴とするパラメータの同定装置。
【請求項３】
請求項２に記載のパラメータの同定装置において、
前記出力評価関数の変化量は、前記物理モデルの出力と最適化対象の出力の誤差の二乗和についての各パラメータ変化に対する勾配ベクトルであり、各パラメータの勾配ベクトル間の相関係数を求めることを特徴とするパラメータの同定装置。
【請求項４】
請求項１〜３のいずれか１つに記載のパラメータの同定装置において、
十分でないと判断された場合に、同定する物理パラメータを変更するか、パラメータの数を減らすか、同定に用いる最適化対象の出力データを変更することを特徴とするパラメータの同定装置。
【請求項５】
非線形な物理モデルを用いてシミュレーションを行い、その出力と最適化対象の出力を比較し、その比較結果に基づき物理モデルのパラメータを繰り返し修正することにより、各パラメータを真値に漸近させていくパラメータの同定処理をコンピュータに実行させるパラメータの同定プログラムであって、
各パラメータを変化させた時の出力評価関数値の変化量について、各パラメータの変化量の同士を比較して両者に強い相関の組み合わせがあるかを全てのパラメータ同士について判定する工程と、
強く相関の組み合わせがある場合に物理モデルまたは最適化対象の出力が物理パラメータを真値に漸近させるための十分な情報を含んでいないと、判断する工程と、
をコンピュータに実行させることを特徴とするパラメータの同定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非線形な物理モデルを用いてシミュレーションを行い、その出力と最適化対象の出力を比較し、その比較結果に基づき物理モデルのパラメータを繰り返し修正することにより、各パラメータを真値に漸近させていくパラメータの同定に関する。
【背景技術】
【０００２】
従来より、モデルを用いたシミュレーションなどが行われており、このためには実態に則したモデルが必要となる。このため、モデルの諸元を決定する必要があり、モデルの諸元を実験データから同定することが行われており、この場合は同定が精度良く行えたかを判定するプロセスが必要である。
【０００３】
統計の分野において、入出力の関係を近時的に一次関数で同定する時に、精度良く計数の同定ができない場合として、“多重共線性”の問題が知られている。これは、例えば式（１）の係数（ａ１，ａ２）を同定するために、時系列データ（ｙ，ｕ１，ｕ２）を用いる場合において、もし時系列データ（ｕ１，ｕ２）の間に強い相関があると、係数（ａ１，ａ２）は正しく同定できない、というものである。
ｙ＝ａ1・ｕ1＋ａ2・ｕ2
【０００４】
ここで、相関の強さは、相関係数を用いて調べることができ、相関係数は、１もしくは−１に近いほど相関が強いことを示す。なお、このような点については、非特許文献１に記載がある。
【０００５】
また、特許文献１には、対象設備や製品プロセスの異常状態を診断し、異常状態の改善のための方法が提案されている。ここでは、異常状態は、対象設備の状態（＝説明変数、式（１）のｕ１，ｕ２に対応する）と異常状態（＝目的変数、式（１）のｙに対応する）の関係を同定することにより診断するが、この場合不要説明変数を含めて同定すると同定精度が悪化する問題がある。そこで、特許文献１では、目的変数に対する説明変数の影響を要因分析（Ｆ値あるいはｔ値を計算）し、不要説明変数を破棄する。
【０００６】
【特許文献１】特開２０００−３１５１１１号公報
【非特許文献１】田中豊，脇本和昌(1985)：多変量統計解析法，現代数学社
【発明の開示】
【発明が解決しようとする課題】
【０００７】
ここで、非特許文献１は、線形のモデルを対象として、複数の入力時系列データ（ｕ１，ｕ２）が存在することを前提として、同定が行えるかを判定している。そこで、非線形のモデルについて、複数の時系列データが存在しない場合には適用できない。
【０００８】
また、特許文献１も非特許文献１と同様に、多入力の説明変数（ｕ１，ｕ２，…）から同定精度の良い説明変数を選択することを課題としており、このような前提が当てはまらないものには適用することができない。
【課題を解決するための手段】
【０００９】
本発明は、非線形な物理モデルを用いてシミュレーションを行い、その出力と最適化対象の出力を比較し、その比較結果に基づき物理モデルのパラメータを繰り返し修正することにより、各パラメータを真値に漸近させていくパラメータの同定装置であって、パラメータを変化させた時の出力評価関数値の変化量について、比較して強い相関の組み合わせがあるかを全てのパラメータ同士の出力評価関数値の変化量について判定し、強く相関の組み合わせがある場合に物理モデルまたは最適化対象の出力が物理パラメータを真値に漸近させるための十分な情報を含んでいないと、判断することを特徴とする。
【００１０】
また、前記出力評価関数の変化量は、前記物理モデルの出力と最適化対象の出力の誤差についてのパラメータ変化の影響であることが好適である。
【００１１】
また、前記出力評価関数の変化量は、前記物理モデルの出力と最適化対象の出力の誤差の二乗和についての各パラメータ変化に対する勾配ベクトルであり、各パラメータの勾配ベクトル間の相関係数を求めることが好適である。
【００１２】
また、十分でないと判断された場合に、同定する物理パラメータを変更するか、パラメータの数を減らすか、同定に用いる最適化対象の出力データを変更することが好適である。
【００１３】
さらに、本発明は、上記装置において実行されることで、上記処理をコンピュータの行わせるプログラムに関する。
【発明の効果】
【００１４】
このように、本発明によれば、出力の状態からパラメータが真値に収束可能か判断できる。また、その判定の際にどのパラメータの組かを特定することができる。
【発明を実施するための最良の形態】
【００１５】
以下、本発明の実施形態について、図面に基づいて説明する。
【００１６】
図２には、本実施形態において、同定の対象となる油圧系の概略構成を示す。エンジンには、自動変速機が接続されており、この自動変速機によりエンジン出力軸の回転が変更されてタイヤに伝達される。自動変速機は、油圧の供給を受け、変速動作を行う。例えば、車両の走行状態に応じて変速タイミングとなった場合には、ＥＣＵ（電子制御ユニット）がこれを検出し、変速指令に応じた電流を油圧制御回路に供給し、油圧制御回路自動変速機内のクラッチパックへの油圧を制御して変速動作が行われる。
【００１７】
ここで、この自動変速機の油圧系は、エンジンで駆動される油圧ポンプを有し、この油圧ポンプからのオイルが油圧ラインに供給される。この油圧ラインの油圧ポンプの近くには油圧調整弁が設けられており、これによって油圧ラインにおけるライン圧（元圧）が所定値に制御される。油圧ラインには、リニアソレノイド弁が設けられており、油圧ラインのライン圧のクラッチパックへの供給を制御している。なお、リニアソレノイド弁からクラッチパックへ至る流路には、オリフィスが設けられ、クラッチパックへの油圧の供給を制限している。また、クラッチパックは、油圧によってピストンをばね圧に抗して移動させることで、入力側のディスクと出力側のディスクの係合・解放を制御するものである。例えば、１速から２速への変速であれば、元々入力側のディスクと係合している１速の出力ディスクを解放し、２速の出力ディスクを係合させる。
【００１８】
このような制御対象となる油圧系を考えると、入力は油圧指令値であり、出力はクラッチパックへの供給油圧になる。
【００１９】
図３には、本実施形態の同定方法の概要を示してある。実機の油圧系出力とモデル出力との誤差を基に、油圧詳細モデル内の（ｉ）〜（ｉｖ）のパラメータを同定することを示している。この例では、パラメータ（ｉ）油圧系の時定数、（ｉｉ）クラッチパックの詰まり位置、（ｉｉｉ）クラッチ作動荷重、（ｉｖ）油圧不感帯幅である。油圧指令値は、油圧詳細モデルに供給される。この油圧詳細モデルは、設計条件モデルと、特性変動モデルからなっている。そして、油圧指令値を入力することによって、油圧詳細モデルからは、クラッチパックに供給される油圧が出力される。また、実機の油圧系にも油圧指令値を供給し実際に出力される油圧を計測する。この実機の出力が最適化対象の出力となる。そして、両者の誤差が算出され、これが同定（最適化）手法に供給される。同定手法では、油圧詳細モデルの中の特性変動モデルの特性パラメータを誤差に応じて変更し、パラメータ（ｉ）〜（ｉｖ）について新しいパラメータを特性変動モデルに設定する。これによって、次の油圧詳細モデルが設定されて、この出力を得るという動作が繰り返される。
【００２０】
また、図４Ａには、油圧系のブロック図が示されている。まず、変速指令がＥＣＵに供給される。ＥＣＵは、変速指令に応じてリニアソレノイド弁にこれが動作するための信号を供給する。なお、実際には、解放側の油圧ラインのリニアソレノイド弁に対する解放用の信号と、係合側の油圧ラインのリニアソレノイド弁に対する係合用の信号の両方である。これによって、オリフィス流路にリニアソレノイド弁の動作に伴うオイルがクラッチパックに供給され、クラッチ圧が決定される。なお、オリフィス流路の流量：ｑとし、クラッチパックにおけるクラッチピストンストローク：ｘｃ、クラッチ面積Ａｃ、クラッチ圧Ｐｃとする。
【００２１】
図４Ｂには、同定する油圧系の諸元について示してある。変速指令が入力されるとＥＣＵでは、この指令に対する遅れが発生する。この例では、遅れは１／（Ｔ・ｓ＋１）で表され、Ｔが時定数であり、これがパラメータ（ｉ）になる。また、リニアソレノイド弁には、動作開始までの不感帯があり、これがパラメータ（ｉｖ）となる。
【００２２】
オリフィスは、その流量ｑｃが、ｑｃ＝Ｃ0・Ａｃ√｛２（Ｐｌ−Ｐｃ）／ρ｝によって決定される。ここで、Ｐｌはオリフィスの上流側圧力、Ｐｃはクラッチパック側圧力、ρはオイルの密度、Ｃ0はオイルの粘度（開口比、レイノルズ数）に関する定数である。
【００２３】
そして、クラッチパックにおいては、そのクラッチピストンストロークｘｃ及びクラッチ圧ＰＣは、ｘｃ＝∫ｑｃｄｔ／Ａｃ，Ｐｃ＝Ｋｃ・ｘｃ／Ａｃで表される。Ｋｃは定数である。さらに、クラッチパックの動作は、パック詰まりによって大きく影響を受ける。そこで、パック詰まり位置をｘｃ＿ｍａｘとして、その位置でのクラッチ圧をＰｃ＿ｍａｘとする。従って、パック詰まり位置において、ｘｃ＝ｘｃ＿ｍａｘ，Ｐｃ＝Ｋｃ・ｘｃ＿ｍａｘ／Ａｃ＝ＰＣ＿ｍａｘとなる。
【００２４】
次に、本実施形態におけるパラメータの同定について、図１Ａ、図１Ｂのフローチャートに従って説明する。なお、説明では、実機の出力データの代わりに、パラメータ（ｉ）〜（ｉｖ）の値を適当に変動させた油圧詳細モデルの出力を用いる。
【００２５】
まず、ステップ１では、同定（最適化）するモデルパラメータの初期（ノミナル）値［ａ（１）  ｂ（１）  ｃ（１）  ｄ（１）］を決定する。このノミナル値は、通常考えられる正しいそうな値とすることが好ましい。なお、この例において、同定の対象とする油圧系諸元を図４Ｂに破線の丸印で示してある。
【００２６】
ステップ２、３では、モデルのシミュレーションを行い、図５のようにモデル出力と実機の出力との誤差を得る。本実施形態では、この誤差の二乗和を評価値としており、これを式（１）または式（２）により計算する。式（１）はアナログデータの連続系の場合であり、式（２）はサンプリングデータによる離散系の場合である。連続系の場合、時間ｔ＝０〜ｔｆまで、出力油圧Ｐｃの誤差の二乗を積分した値であり、離散系の場合、Ｎ個のサンプルの誤差の二乗和のＮ分の１となっている。
【００２７】
【数１】
0004851396.tif 000002
【００２８】
続いてステップ４において、式（１）または式（２）の誤差（評価値：誤差の二乗和）が閾値以内なら、目標とする出力が得られているため同定のステップは終了する。
【００２９】
ステップ４においてＮＯ、すなわち誤差（評価値）が閾値より大きい場合には、ステップ５に進み、モデルパラメータの初期値を適当にずらしたものを新たなパラメータに設定する（式（３））。
[a(1) b(1) c(1) ...]=[a(1)+δ1 b(1)+δ2 c(1)+δ3 ...]    （３）
【００３０】
ステップ６では、式（３）のパラメータにより、ステップ２と同様にモデルシミュレーションを行い、ステップ７で式（１）または式（２）により再度誤差（評価値）を計算して、ステップ８で誤差（評価値）が閾値以内か判定する。ステップ８の判定で閾値以内なら、同定のステップは終了する。一方、ステップ８の判定で、誤差（評価値）が閾値より大きい場合には、ステップ９へ進む。
【００３１】
ステップ９では、各パラメータの要因効果から勾配ベクトルを計算する。このために、まず図３、図４Ｂに示した４つのパラメータ（ｉ）〜（ｉｖ）に対して、表１に示すように２水準（１つのパラメータについて２つの値）のパラメータの組合せ、計１６通り用意し、各組合せごとに誤差の二乗和Ｊ１〜Ｊ１６を計算する。
【００３２】
【表１】
0004851396.tif 000003
【００３３】
求めた誤差から、表２に示すように、各パラメータの水準（２水準）ごとに、要因効果（水準ごとの平均値）を計算する。
【００３４】
【表２】
0004851396.tif 000004
【００３５】
すなわち、Ｊa（１）は、パラメータ（油圧系時定数）がａ（１）（水準１）である時の誤差の二乗和の平均値（評価値）、Ｊa（２）は、パラメータ（油圧系時定数）がａ（２）（水準２）である時の誤差の二乗和の平均値である。Ｊｂ（１），Ｊｂ（２），Ｊｃ（１），Ｊｃ（２），Ｊｄ（１），Ｊｄ（２）も同様に各パラメータについての水準１、２の際の評価値の平均値である。従って、この表２から各パラメータについて水準１，２の変化に対する効果がわかる。
【００３６】
このような表２の要因効果から、式（４）で定義するパラメータ全体の勾配ベクトル［Ｓａ（１）  Ｓｂ（１）  Ｓｃ（１）  Ｓｄ（１）］が求められる。すなわち、上記勾配ベクトルの１つ１つの要素が１つのパラメータの水準１，２に対する出力の勾配である。
【００３７】
【数２】
0004851396.tif 000005
【００３８】
ステップ１０では、勾配ベクトルＳ（ｉ）に基づき、パラメータの修正量を計算する。ここで、ｉは何回目の計算かを示しており、上記水準１，２は１回目、２回目の計算に用いたものを表しており、３回目の計算ではｉ＝３であり、水準１がｉ＝２、水準２がｉ＝３となる。
【００３９】
次のモデルパラメータは、スカラ値ｋを用いて式（５）により計算する。
【００４０】
【数３】
0004851396.tif 000006
【００４１】
このようにして、各パラメータについて、評価値に応じた修正量を加算して、次のパラメータが決定される。そして、次のパラメータが決定された場合には、ステップ６に戻って、誤差が閾値以下になるまで同じステップを繰り返す。
【００４２】
このようにして、誤差が所定値以下となるパラメータ（＝油圧諸元）を同定することができる。
【００４３】
「物理パラメータのチェック」
次に、このような油圧諸元の同定について、対象とする物理パラメータが適切かどうかを判定する処理について図１Ｂに基づいて説明する。
【００４４】
上述した図１Ａの処理により物理パラメータの繰り返し同定を実施する。まず、繰り返し回数ｎが予め設定した所定回数ｎ＿ｍｉｎ以上かを判定する（ステップ１１）。この判定結果においてＮＯであれば処理を終了する（適切かどうかを判定する処理はしない）。
【００４５】
そして、ステップ１１の判定においてＹｅｓとなった場合には、各勾配ベクトル間の相関係数γを演算する（ステップ１２）。
【００４６】
すなわち、各パラメータを変化させた時の出力評価関数値（例えば誤差の二乗和）の変化量として、式（４）の勾配ベクトル［Ｓａ（１）  Ｓｂ（１）  Ｓｃ（１）  Ｓｄ（１）］に着目し、各勾配ベクトル間の相関係数を図１Ｂ中ｍｐ式（６）のように求める。なお、式（６）においては、２つのパラメータ間の相関係数について示したが、図１Ｂに示すように、全てのパラメータ間の相関係数を演算算出する。
【００４７】
【数４】
0004851396.tif 000007
【００４８】
次に、全相関係数γi,jの値がγのしきい値（例えば、それ以上では強い相関と考えられるしきい値０．９）と比較する（ステップ１３）。この判定で、ＮＯであれば、問題となるパラメータは存在せず、処理を終了する。
【００４９】
一方、ステップ１３の判定でＹｅｓであれば、同定精度が悪いと判定する（ステップ１４）。そして、この場合には、（ｉ）同定パラメータの変更、（ｉｉ）同定パラメータの数の削減、（ｉｉｉ）同定に用いる実機出力データの変更、のいずれかを行う。なお、この（ｉ）〜（ｉｉｉ）の対処は、相関があるとされたパラメータを考慮して行う。例えば、（ｉ）（ｉｉ）であれば、相関が高かった２つのパラメータの一方を変更または削除する。また、（ｉｉｉ）の場合は、相関の高かった２つのパラメータの特性を考慮して、これらが独立して作用するように実機出力データを変更する。
【００５０】
「同定がうまくいく具体例」
まず、油圧系の諸元値が精度良く同定できる場合を示す。図６に同定前・後におけるモデルと実機の時間応答波形の比較を、図７に同定パラメータと誤差の二乗和の探索履歴を示す。モデルと実機はほぼ一致し、約１５回の探索の後、同定パラメータは真値に収束していることがわかる。
【００５１】
この時、各パラメータを変化させた時の出力評価関数値（例えば誤差の二乗和）の変化量として、式（４）の勾配ベクトル［Ｓａ（１）  Ｓｂ（１）  Ｓｃ（１）  Ｓｄ（１）］に着目し、各勾配ベクトル間の相関係数を上述の式（６）のように求める。
【００５２】
このようにして求めた各パラメータについての相関係数を表３に示す。
【００５３】
【表３】
0004851396.tif 000008
【００５４】
どの勾配ベクトル間の相関係数も１または−１に近くはない（相関係数の絶対値が約０．９より小さい値）ので、精度良く同定ができていると判定できる。なお、図７では、検証のために真値がわかっているモデルを対象に同定しているので真値に収束したか判定できるが、実際には当然わからない。
【００５５】
次に同定がうまくいかない場合の例を２つ示す。
【００５６】
「同定がうまくいかない場合１」
この例では、図３および図４Ｂにおいて、４つの同定パラメータのうち“油圧の不感帯の幅”を、“油圧応答のむだ時間”に変更する。
【００５７】
図８に同定前・後におけるモデルと実機の時間応答波形の比較を、図９に同定パラメータと誤差の二乗和の探索履歴を示す。図８からわかるように、モデルと実機の油圧出力は系はほぼ一致している。ところが、図９からわかるように、同定パラメータは真値に収束していない。
【００５８】
この時、式（４）の各勾配ベクトル間の相関係数を、式（６）により求めた結果を表４に示す。
【００５９】
【表４】
0004851396.tif 000009
【００６０】
勾配ベクトル（油圧系時定数）と（無駄時間）間の相関係数が０．９７と１に近いので、２つのパラメータは誤差の二乗和（評価関数）に対して同じ効果をもっていることがわかり、よって精度良く同定ができていないと判定できる。
【００６１】
「同定がうまくいかない場合２」
この例では、図３および図４Ｂにおいて、４つの同定パラメータに加え、パラメータ“オリフィスの流量係数”も同時に同定する場合を考える。
【００６２】
図１０に同定前・後におけるモデルと実機の時間応答波形の比較を、図１１に同定パラメータと誤差の二乗和の探索履歴を示す。この例においても、モデルと実機の出力波形はほぼ一致しているにもかかわらず、同定パラメータは真値に収束していない。
【００６３】
この時の、式（４）の各勾配ベクトル間の相関係数を、式（６）により求めた結果を表５に示す。
【００６４】
【表５】
0004851396.tif 000010
【００６５】
勾配ベクトル（パック詰り位置）と（オリフィスの流量係数）間の相関係数が−０．９９と−１に近い。このため、２つのパラメータは誤差の二乗和（評価関数）に対して同じ効果をもっていることがわかり、よって精度良く同定ができていないと判定できる。
【００６６】
「同定パラメータが真値に収束しないことの問題点」
最後に、同定パラメータが真値に収束していない場合に、どのような問題が生じるかについて説明する。
【００６７】
表６は、図１０、１１と同じ問題設定で、パラメータを同定した結果の一例である。
【００６８】
【表６】
0004851396.tif 000011
【００６９】
表６の同定パラメータを適用し、同定したのと同一の電流指令を用いてシミュレーションした結果を図１２に示す。同定したパラメータ値は真値とは異なっているにもかかわらず、油圧の時系列波形は、実機相当とモデルでほぼ一致している。
【００７０】
しかし、同定した時と別の電流指令を用いてシミュレーションした場合には、図１３に示すようにモデルの時系列波形が、実機相当の波形とずれてしまう。このように、同定パラメータが真値に収束しない場合には、パラメータの適切な設定が行われておらず、特定の条件のシミュレーションにしか対応できないという問題があることがわかる。
【００７１】
「本実施形態の作用効果」
このように、本実施形態では、同定精度悪化の原因である各パラメータ間の相関の強さを定量的に評価するため、各同定パラメータを変化させた時の出力の評価関数値（例えば誤差の二乗和）の変化量に着目した。そして、この変化量が各パラメータで同じ傾向を持つ（相関が大きい）場合に、同定が適切に行えないと判断する。これは、各パラメータで同じ傾向を持つ（相関が大きい）ほど（相関が強くなるほど）、両者を区別して同定ができなくなり、同定精度が悪化すると考えられるためである。
【００７２】
従って、本実施形態によれば、同定パラメータが真値に収束可能か判断できる。また、どの同定パラメータの組が、相関が大きいかに基づき、真値に収束しない原因かがわかる。すなわち，あるパラメータの組が評価関数値（例えば誤差の二乗和）に対して同じ効果を持つと同定パラメータは真値に収束できないが、この原因となるパラメータの組を特定できる。そこで、判定結果に応じて同定パラメータを変更したり、その数を減らすことで適切な同定作業を行うことが可能となる。
【図面の簡単な説明】
【００７３】
【図１Ａ】パラメータ同定の処理を示すフローチャートである。
【図１Ｂ】パラメータチェックの処理を示すフローチャートである。
【図２】油圧制御系の概要を示す図である。
【図３】パラメータの同定方法の概要を示す図である。
【図４Ａ】油圧系のブロック線図である。
【図４Ｂ】同定する油圧系の諸元を示す図である。
【図５】モデルのシミュレーション結果の例を示す図である。
【図６】同定結果１：時系列波形を示す図である。
【図７】同定結果２：パラメータ値と誤差の二乗和を示す図である。
【図８】同定結果１：時系列波形を示す図である。
【図９】同定結果２：パラメータ値と誤差の二乗和を示す図である。
【図１０】同定結果１：時系列波形を示す図である。
【図１１】同定結果２：パラメータ値と誤差の二乗和を示す図である。
【図１２】同定したときと同じ電流指令パターンを用いた場合の結果を示す図である。
【図１３】同定したときと等しい電流指令パターンを用いた場合の結果を示す図である。
【図１Ａ】
0004851396.tif 000012
【図１Ｂ】
0004851396.tif 000013
【図２】
0004851396.tif 000014
【図４Ａ】
0004851396.tif 000015
【図４Ｂ】
0004851396.tif 000016
【図５】
0004851396.tif 000017
【図６】
0004851396.tif 000018
【図７】
0004851396.tif 000019
【図８】
0004851396.tif 000020
【図９】
0004851396.tif 000021
【図１０】
0004851396.tif 000022
【図１１】
0004851396.tif 000023
【図１２】
0004851396.tif 000024
【図１３】
0004851396.tif 000025
【図３】
0004851396.tif 000026