(57)【要約】  動的システムにおける追従誤差を最小化するための方法であって、動的システムの所望の軌道データを取得することと、動的システムの所望の軌道および実際の軌道のうちの少なくとも一方に対する一組の制約を取得することと、既知の原Ｂ−スプライン基底関数を有するが未知のＢ−スプライン係数を有する一組の一様または非一様な有理Ｂ−スプラインを取得することと、順方向フィルタ処理をＢ−スプライン基底関数に適用することならびに原Ｂ−スプライン基底関数およびフィルタ処理されたＢ−スプライン基底関数を利用してＢ−スプラインの最適な係数を選択することを含む、所望の軌道データに軌道最適化プロセスを適用することと、所望の軌道および実際の軌道のうちの少なくとも一方に対する一組の制約を満たしながら、得られる実際の軌道が所望の軌道に実質的に等しくなるように、軌道最適化プロセスに応答して動的システムに対して最適な運動コマンド信号を出力することと、を含む、方法。
2018522356.tif 000002
【特許請求の範囲】
【請求項１】
動的システムにおける追従誤差を最小化するための方法であって、
前記動的システムの所望の軌道データを取得することと、
前記動的システムの前記所望の軌道および実際の軌道のうちの少なくとも一方に対する一組の制約を取得することと、
既知の原Ｂ−スプライン基底関数および少なくとも１つの未知のＢ−スプライン係数を有する、一組の一様または非一様な有理Ｂ−スプラインを取得することと、
順方向フィルタ処理をＢ−スプライン基底関数に適用することならびに前記原Ｂ−スプライン基底関数および前記フィルタ処理されたＢ−スプライン基底関数を利用して前記Ｂ−スプラインの最適な係数を選択することを含む、前記所望の軌道データに軌道最適化プロセスを適用することと、
前記所望の軌道および前記実際の軌道のうちの前記少なくとも一方に対する前記一組の制約を満たしながら、得られる実際の軌道が前記所望の軌道に実質的に等しくなるように、前記軌道最適化プロセスに応答して前記動的システムに対して最適な運動コマンド信号を出力することと、を含む、方法。
【請求項２】
軌道最適化プロセスを適用する前記ステップが、前記所望の軌道データの部分集合を使用して遂行され、それにより限定的先読み機能になる、請求項１に記載の方法。
【請求項３】
前記動的システムが、ＣＮＣ送り駆動装置であって、その所望の軌道および実際の軌道が、前記ＣＮＣ送り駆動装置の所望の位置および実際の位置を表す、ＣＮＣ送り駆動装置、を表し、前記制約が、前記ＣＮＣ送り駆動装置の前記所望の軌道および前記実際の軌道のうちの前記少なくとも一方に対する運動力学的な制約を表す、請求項１に記載の方法。
【請求項４】
前記運動力学的な制約が、位置制限、速度制限、加速度制限、および躍度制限のうちの少なくとも１つを含む、請求項３に記載の方法。
【請求項５】
前記軌道最適化プロセスが、前記制約を条件として前記所望の軌道に追従する際の前記誤差の２ノルムを最小化する、請求項１に記載の方法。
【請求項６】
前記出力される最適な運動コマンドが、ＣＮＣ送り駆動装置に渡される変更された位置コマンドである、請求項１に記載の方法。
【請求項７】
前記ＣＮＣ送り駆動装置が、不要な振動モードを有する、請求項６に記載の方法。
【請求項８】
ノットベクトル、制御点の数、および重みのうちの少なくとも１つを含む他のＢ−スプラインパラメータの最適化をさらに含む、請求項１に記載の方法。
【請求項９】
２ノルム以外の他の指標およびノルムに基づく最適化をさらに含む、請求項１に記載の方法。
【請求項１０】
前記動的システムおよび前記Ｂ−スプライン基底関数に対して制約を課すことをさらに含む、請求項１に記載の方法。
【請求項１１】
複数の動的システム、所望の軌道、およびＢ−スプラインを組み合わせることをさらに含む、請求項１に記載の方法。
【発明の詳細な説明】
【技術分野】
【０００１】
関連出願の関連出願の相互参照
本出願は、２０１５年７月２８日に出願された米国仮出願第６２／１９７，６８２号の利益を主張する。上記出願の開示の全体が、参照により本明細書に組み込まれる。
【０００２】
政府所有権
本発明は、米国国立科学財団により授与されたＣＭＭＩ１３５０２０２の下で政府の支援をうけてなされた。政府は、本発明に対して特定の権利を有する。
【０００３】
本開示は、運動誤差に関し、より具体的には、サーボ誘発性運動誤差を補償するための、一様または非一様な有理Ｂ−スプライン曲線から導出されたフィルタ処理された基底関数の使用に関する。
【背景技術】
【０００４】
本セクションは、必ずしも先行技術ではない、本開示に関係する背景情報を提供する。本セクションは、本開示の全般的な概要を提供するものであり、本開示の完全な範囲または本開示の特徴の全ての網羅的な開示ではない。
【発明の概要】
【発明が解決しようとする課題】
【０００５】
本教示は、フィルタ処理されたベーシススプライン（Ｂ−スプライン）を使用することにより、動的システム（非最小位相（ＮＭＰ）を有するものを含む）ゼロにおける追従誤差を最小化するための手法を提供する。
【課題を解決するための手段】
【０００６】
追従コントローラの出力は、未知の係数を有するＢ−スプラインの一次結合として表される。基底関数は、システムの動特性を使用して順方向にフィルタ処理され、基底関数の係数が、所与の軌道に追従する際の誤差を最小化するために選択される。数多くの他の利益のなかでも、Ｂ−スプラインの使用は、（ｉ）システムの入力および出力に対してスプライン係数の一次関数として制約が課されること、（ｉｉ）非ゼロの初期条件が課されること、（ｉｉｉ）より高階数の導関数の追従、ならびに（ｉｖ）フィルタ処理プロセスにおける限定的先読み能力の使用、を可能にする。
【０００７】
さらなる応用の分野は、本明細書で提供する説明から明らかになるであろう。本概要中の説明および具体的な例は、例示目的のみを意図し、本開示の範囲を限定することを意図しない。
本明細書で説明する図面は、選択した実施形態の例示目的のみのためのものであり、本開示の範囲を限定することを意図しない。
【図面の簡単な説明】
【０００８】
【図１】不要な振動モードを含むＣＮＣ送り軸の動特性を示す。
【図２】提案の手法の流れ図である。
【図３】実験的な構成、すなわち、可撓性の取付品を有する２軸ステージを示す。
【図４】実験的な構成のＸ軸およびＹ軸に関して、測定された周波数応答関数およびモデル化された周波数応答関数を示す。
【図５】実験のために使用される所望の経路を示す。
【図６】非補償の方法および提案の方法に関して、Ｘコマンド信号およびＹコマンド信号を示す。
【図７】非補償の軌道および提案の軌道に関して、取付品の測定された加速度信号の高速フーリエ変換（ＦＦＴ）を示す。
【図８】非補償の方法および提案の方法に関して、Ｘ軸およびＹ軸の位置追従誤差を示す。
【図９】非補償の軌道および提案の軌道に関して、輪郭誤差を示す。
【図１０】限定的先読み能力を分析するためのシミュレーションに使用される所望の経路を示す。
【図１１】提案の方法に関して、限定的先読みがある場合の追従誤差およびない場合の追従誤差を示す。
【発明を実施するための形態】
【０００９】
対応する参照数字は、図面のいくつかの図の全体にわたって、対応する部分を表す。
【００１０】
これより、添付図面を参照しながら、例示的な実施形態をより十全に説明する。例示的な実施形態は、本開示が徹底的なものになるように、かつ範囲を当業者に十全に伝達するように、提供される。本開示の実施形態の徹底的な理解を提供するために、具体的なコンポーネント、デバイス、および方法の例などの多数の具体的な詳細が述べられている。具体的な詳細は用いられなくてもよいこと、例示的な実施形態は数多くの異なる形態で実施され得ること、および例示的な実施形態は本開示の範囲を限定すると解釈されるべきではないことは、当業者に明らかであろう。一部の例示的な実施形態では、よく知られたプロセス、よく知られたデバイス構造、およびよく知られた技術は、詳細には説明されていない。
【００１１】
本明細書で使用する術語は、特定の例示的な実施形態のみを説明する目的のためのものであり、限定的であることを意図しない。本明細書で使用する場合、単数形「ａ」、「ａｎ」、および「ｔｈｅ」は、文脈上明らかに必要とされない限り、複数形も含むことを意図し得る。用語「含む（ｃｏｍｐｒｉｓｅｓ）」、「含む（ｃｏｍｐｒｉｓｉｎｇ）」、「含む（ｉｎｃｌｕｄｉｎｇ）」、および「有する（ｈａｖｉｎｇ）」は、包括的であり、したがって、言明した特徴、整数、ステップ、動作、要素、および／またはコンポーネントの存在を明示するが、１つ以上の他の特徴、整数、ステップ、動作、要素、コンポーネント、および／またはそれらの群の存在または追加を排除しない。本明細書に記載の方法ステップ、プロセス、および動作は、遂行の順序として明確に特定されない限り、論じた特定の順序または例示した特定の順序でそれらが遂行されることを必ず必要とすると解釈されるべきではない。追加のステップまたは代替的なステップが用いられてもよいことも、理解されるべきである。
【００１２】
入力または出力の制約を含む動的システムの追従
【００１３】
追従制御の目標は、システムの出力が所望の軌道をできるだけ接近して追尾するように強いることである。図１に示す、伝達関数Ｇを有する線形時不変な（ＬＴＩ）動的システムを考慮されたい。追従制御の目標は、χおよび／もしくはχｄ（またはそれらの時間微分）に対して課される制約を満たしながら、実際の軌道χが必ず所望の軌道χｄをできるだけ接近して追従するようにすることである。
【００１４】
さらに、Ｇが次式により与えられる不要な振動モードを有するＣＮＣマシンの送り軸の動特性を表す、具体的な例を考慮されたい。
【数１】
2018522356.tif 000003
【００１５】
式中、Ｍは、モードの数を表し、αｉおよびβｉは、残留物定数であり、ζｉは、減衰比であり、ωｎ，ｉは、ｒａｄ／ｓでのｉ番目のモードの固有振動数である。χｄの代わりに送り軸に適用されるとき、マシン軸の運動力学的な限界を満たしながらＧにより引き起こされる追従誤差を最小化するように、追従制御を使用してχｄｍ＊、つまり所望の位置軌道χｄの変更され最適化されたバージョンを生成することが望まれる。
【００１６】
図２は、フィルタ処理されたＢ−スプラインを使用してχｄｍ＊を生成するための、本教示による手法のブロック図を示す。χｄ＝［χｄ（０）χｄ（１）．．．χｄ（Ｅ）］Ｔを、ＣＮＣにより補間される工具経路のｘ成分のＥ＋１の離散時間ステップ（すなわち、所望の軌道）を表すものとしよう。我々は、ＣＮＣが、χｄのＥ＋１時間ステップが予め既知であるような先読み能力を有すると仮定する。さらに、我々は、変更されたが最適化されていない運動コマンドχｄｍ＝［χｄｍ（０）χｄｍ（１）．．．χｄｍ（Ｅ）］Ｔが、次式のようになるように、Ｂ−スプラインを使用してＣＮＣにおいてパラメータ化されると仮定する。
【数２】
2018522356.tif 000004
【００１７】
ｎ＋１≦Ｅ＋１の制御点ｐ０、ｐ１、．．．、ｐｎおよびノットベクトル［ｇ０ｇ１．．．ｇｍ＋ｎ＋１］Ｔを有する次数ｍのＢ−スプラインの場合、Ｎｊ，ｍ（ξ）は、次式により与えられる実数値の基底関数である。
【数３】
2018522356.tif 000005
【００１８】
式中、ｊ＝０、１、．．．、ｎであり、ξ∈[０、１]は、正規化時間を表すスプラインパラメータであり、正規化時間は、式（２）では、Ｅ＋１の一様間隔の点ξ０、ξ１、．．．、ξＥに離散化されている。ノットベクトルは、次式のようになるように、一様になるように選択することができる。
【数４】
2018522356.tif 000006
【００１９】
式（２）で与えられるχｄｍがχｄの代わりに運動コマンドとして図１のＣＮＣ送り軸に適用される場合、制御点ベクトルｐは、ｐは時間の関数ではないため、Ｇによる影響を受けないことに注意されたい。しかしながら、時間の関数である基底関数行列Ｎは、Ｇによりフィルタ処理される。したがって、得られる位置χ＝［χ（０）χ（１）．．．χ（Ε）］Ｔは、次式により近似することができる。
【数５】
2018522356.tif 000007
【数６】
2018522356.tif 000008
【００２０】
Ｂ−スプラインの１つの有用な数学的特性は、Ｂ−スプライン曲線の導関数もまた、原曲線の制御点を用いて線形に表すことができるＢ−スプライン曲線であることである。ゆえに、χｄｍのｒ番目の導関数は、次式のように書くことができる。
【数７】
2018522356.tif 000009
【００２１】
式中、Ｎｒは、位置のｒ番目の導関数の基底関数行列であり、その要素は、Ｎの要素から容易に計算することができる。したがって、本教示の言明した目標を達成するために、次式により与えられる目的関数を最小化するように計算された最適な制御点ｐ＊を使用して、χｄｍ＊＝Ｎｐ＊を生成することができる。
【数８】
2018522356.tif 000010
【００２２】
式中、Ｌｒは、χｄｍのｒ番目の導関数に関連付けられた運動力学的な限界（例えば、その速度限界、加速度限界、および躍度限界）を表し、｜｜・｜｜∞は、ベクトルの無限ノルム（すなわち、要素の最大絶対値）を意味する。
【００２３】
式（８）中のＮｒを
2018522356.tif 000011
（すなわち、Ｎｒのフィルタ処理されたバージョン）と交換することにより、制約を出力χおよびその導関数に対して同様に適用することができることに留意されたい。
【００２４】
運動力学的な制約が考慮されない場合、式（８）の解は、次式の従来の最小二乗フィッティング問題に単純化されることに注意されたい。
【数９】
2018522356.tif 000012
【００２５】
しかしながら、制約が含まれた状態では、式（８）は、制約が実行可能な限り保証された解を有する不等式制約を含む最小二乗問題になる。行列
2018522356.tif 000013
が最大階数である場合（すなわち、フィルタ処理された基底関数が線形独立である場合）、解は、唯一の解である。本方法は、業界標準のＢ−スプライン技法および基本的な最小二乗技法を使用して、軌道のパラメータ化、制約の適用、および最適化を同時に達成し、それを使用しやすくする。χｄおよびχｄｍ＊の両方が同じ長さ（すなわち、Ｅ＋１）を有することにも注意されたい。これは、提案の方法が運動スピードを犠牲にしないことを意味する。本方法を単一の軸（すなわち、χ）に関して論じてきたが、本方法は、各軸を別々に扱うことにより、ＣＮＣマシンの複数の軸に容易に適用可能である。したがって、本方法は、非常に体系的かつスケーラブルである。
【００２６】
非ゼロの初期条件の取り扱い
【００２７】
追従制御では、追従されるべき軌道は、軌道自体およびそのより高階数の導関数と関連付けられた非ゼロの初期条件を有し得る。これは、式（７）では、制御点とより高階数の導関数との間の線形関係を通じて実現される。例えば、初期時間ξ０（すなわち、ξ＝ξ０であるとき）におけるχｄｍの第１の導関数は、次式として表される。
【数１０】
2018522356.tif 000014
【００２８】
式（７）のより高階数の導関数は、類似した様式により得ることができる。ξ＝ξ０における第１の導関数の値は、最初の２つの制御点に依存することに留意されたい。第２の導関数の値は、最初の３つの制御点に依存することになり、同じように続く。その結果、所与の軌道のｒ番目の導関数に対する初期条件は、軌道に関連付けられた最初のｒ＋１個の制御点（ｒ＜ｍ）を制約する。ゆえに、χｄの初期条件（これは、ｘｄｍの初期条件と同じである）は、次式として表すことができる。
【数１１】
2018522356.tif 000015
【００２９】
式中、ｚは、初期条件のベクトルであり、Ｑは、式（７）および（１０）、ならびにより高階数の導関数のために導出された類似した式から得られた係数を含む行列である。それゆえに、式（９）および（１１）を組み合わせて、次式を得ることができる。
【数１２】
2018522356.tif 000016
【００３０】
式中、
2018522356.tif 000017
は、ラグランジュ乗数のベクトルである。式（１２）の解は、（式（９）に類似した様式で）、式（１１）により与えられた条件を満たしながら、最小の追従誤差をもたらす。
【００３１】
より高階数の導関数の追従向上
【００３２】
特定の追従用途、例えばロボットおよびＣＮＣマシンにおける位置追従では、所望の軌道、例えば速度および加速度のより高階数の導関数の綿密な追従を維持することも、目的となり得る。Ｂ−スプラインのフレームワークは、この観点において非常に有用である。
【００３３】
Ｇを使用してＮ１、Ｎ２、．．．Ｎｒをフィルタ処理すると、それぞれ
2018522356.tif 000018
となる。その結果、Ｘｄおよびその最初のｒ個の導関数の正確な追従は、次式により与えられる構成方程式を解くことにより達成することができる。
【数１３】
2018522356.tif 000019
【００３４】
式中、α１、α２、．．．αｒは、Ｘｄの各導関数に与えられたユーザ定義された重みであり、位置追従に関する各導関数の重要性を示す。
式（１４）の解は、式（９）に類似した様式で最小二乗を使用して得ることができる。
【００３５】
限定的先読み
【００３６】
提案の方法をここまで説明する上で、軌道ｘｄｍのＥ＋１個の全てのポイントが既知であること、および最適な係数ｐ＊がｘｄｍのＥ＋１個の全てのポイントを使用して解かれていることを、仮定してきた。基本的に、完全（全体）先読み能力を仮定してきた。
【００３７】
しかしながら、実務上、ｘｄｍは、その全体が知られているわけではない場合がある。そうではなくて、ｘｄｍは、順次に組み合わさって全体としてｘｄｍを形成するより小さいバッチ（部分集合）で入手可能な場合がある。ｘｄｍの全てが入手可能なときでさえ、ｐ＊を生成するためにｘｄｍの全てを一度に使用することは、（例えば、リアルタイム実装にとって）計算上法外に高価な場合がある。
【００３８】
このような場合、提案の方法の適用は、ｘｄｍのより小さい部分集合を使用して順次（再帰的に）遂行することができ、したがって、完全先読み能力に対立するものとしての限定的先読みを提供した。
【００３９】
これをするための１つの方法は、上述の行列Ｎおよび
2018522356.tif 000020
を生成し、次いでそれらを、ｘｄｍの入手可能な任意の部分集合に一致するように、正しいサイズのより小さい部分行列に分割することである。次いで、最適な係数ベクトルｐ＊の対応する部分集合を、より小さい部分行列を使用するために解く。このプロセスは、ｘｄｍの後続の部分集合に関して（それらが入手可能になるにつれて）順次繰り返される。
【００４０】
ｘｄｍの各バッチからの解の間の連続性は、ｄｍの連続する部分集合の部分行列および係数の間に多少の重なり合いを作成することにより、保証することができる。
【００４１】
実験例
【００４２】
本セクションでは、提案の技法を非補償の場合と比較するために、実験結果を提示する。図３に示す２軸（Ｘ−Ｙ）リニアモータ駆動ステージ（Ａｅｒｏｔｅｃｈ  ＡＬＳ  ２５０１０）を、実験のために使用する。ステージは、速度および加速度のフィードフォワードにより拡張された、従来型のＰ／ＰＩフィードバックコントローラを使用して制御される。コントローラを、１０ｋＨｚのサンプリング周波数を有するｄＳＰＡＣＥ  １１０３リアルタイム制御ボード上に実装する。その両方の軸の速度限界および加速度限界は、それぞれＬ１＝１００ｍｍ／ｓおよびＬ２＝８０００ｍｍ／ｓ２である。図３に示すように、ステージは、棒の上に装着されたブロックからなる取付品を備えている。ブロックは、ｘ位置およびｙ位置が、固有の構造的可撓性にもかかわらず、それらのそれぞれの所望の軌道ｘｄおよびｙｄを正確に追従すると予想される装置（例えば、工具、被加工物、または測定装置）を表すと仮定する。取付品の加速度は、図３に示す２つの一方向加速度計（ＰＣＢ  Ｐｉｅｚｏｔｒｏｎｉｃｓ  ３９３Ｂ０５）を使用して測定される。
【００４３】
図４は、掃引正弦加速コマンドをステージに適用し、取付品の対応する加速度（すなわち、および）を加速度計を使用して測定することにより生成された、ステージの各軸の動特性Ｇの周波数応答関数（ＦＲＦ）を示す。各軸の動特性は、図４では非常に明白である、２つの優勢モードを呈する。
【００４４】
各軸に関して、より低いモードは、ステージからのモードであり、より高いモードは、可撓性の取付品からのモードである（装着棒が四角形の断面を有するため、取付品の固有振動数は、ｘ方向とｙ方向とで異なる）。しかしながら、図４では明白ではないのが、各ＦＲＦ中の２つのはるかに小さい優勢モードである。言い換えれば、各軸の動特性には４つの不要なモードが存在し（Ｍ＝４）、それらのパラメータを、有理分数多項式法を使用して特定し、１にまとめる。図４は、Ｘ軸に関して、モデル化されたＦＲＦと実際のＦＲＦとの間の優れた適合を示している。Ｙ軸の２つの優勢モードは相対的により減衰しており、共鳴振動数がより近いため、Ｙ軸のモデルは、Ｘ軸のモデルほど正確ではなく、それらを曲線当嵌めするのはより困難である。しかしながら、Ｙ軸の当嵌め誤差は、実務上発生する可能性が高いいくつかのモデリング誤差の存在下で提案の方法の性能をテストする機会を提供する。
【表１】
2018522356.tif 000021
【００４５】
図５は、実験の所望の経路として使用される蝶形曲線を示す。図６は、所望の経路の位置軌道ｘｄおよびｙｄ（すなわち、非補償の運動コマンド）、ならびにその速度プロファイルおよび加速度プロファイルを時間の関数として示す。所望の軌道の持続時間は１ｓ（すなわち、１０ｋＨｚのサンプリング周波数に基づき、Ｅ＝１０４）であることに注意されたい。非補償の軌道を、提案の方法から得られた対応する変更された軌道と比較する。提案の方法を、ｎ＋１＝５１個の制御点を有するｍ＝５の階数のＢ−スプラインで実行する。図６から、提案の方法により生成された変更された軌道が、課される運動力学的な制約を自動的に満たし、いかなる遅延も有さない（すなわち、変更された軌道は、所望の軌道と厳密に同じ長さだけ存続する）ことに注意されたい。
【００４６】
非補償の方法および提案の方法からの位置軌道を、運動コマンドとしてステージのＸ軸およびＹ軸に適用する。図７は、（加速度計を使用して測定された）加速度信号の高速フーリエ変換（ＦＦＴ）を示す。それぞれ非補償の方法および提案の方法に関して、Ｘ軸の加速エネルギースペクトル密度は、６．３８×１０６ｍｍ２／ｓ４および１．９２×１０６ｍｍ２／ｓ４であり、Ｙ軸の加速エネルギースペクトル密度は、７．３２×１０６ｍｍ２／ｓ４および１．８１×１０６ｍｍ２／ｓ４である。これは、提案の方法が、非補償の方法よりも生成する振動が少ないことを意味する。このことを証明するために、３つの全ての場合の追従誤差、ｘｄ−ｘ、およびｙｄ−ｙを、観測器を使用して測定した加速度信号から導出した位置信号（ｘおよびｙ）に基づいて、図８にプロットする。それぞれ非補償の方法および提案の方法に関して、Ｘ軸のＲＭＳ追従誤差は、０．３０１４ｍｍおよび０．０８５２ｍｍであり、Ｙ軸のＲＭＳ追従誤差は、１．０１６５ｍｍおよび０．１１９８ｍｍである。図９は、３つの方法の輪郭誤差を比較するものである。それらのＲＭＳ値は、非補償の方法および提案の方法に関して、それぞれ０．５５０８ｍｍおよび０．０８８０ｍｍであり、追従誤差を最小化することを介して間接的に輪郭誤差を低減する提案の方法の能力を示している。
【００４７】
限定的先読みシミュレーションに関する例
【００４８】
図１１は、限定的先読みありで、および限定的先読みなしで蝶形状軌道（図１０に示す）を追従する位置決め性能を示し、統計を表２で提供する。追従性能は、同一であるが、限定的先読みの実装により計算時間の９６％が節約されている。
【表２】
2018522356.tif 000022
【図１】
2018522356.tif 000023
【図２】
2018522356.tif 000024
【図３】
2018522356.tif 000025
【図４】
2018522356.tif 000026
【図５】
2018522356.tif 000027
【図６】
2018522356.tif 000028
【図７】
2018522356.tif 000029
【図８】
2018522356.tif 000030
【図９】
2018522356.tif 000031
【図１０】
2018522356.tif 000032
【図１１】
2018522356.tif 000033
【国際調査報告】
2018522356.tif 000034
2018522356.tif 000035