(57)【要約】【課題】計測対象物のあらゆる位置においてひずみおよび変位を計測することができる変位計測方法および変位計測装置を提供する。【解決手段】変位計測方法は、環構造体１０の周囲に配置された光ファイバ２０ａおよび２０ｂから、環構造体１０のひずみに基づく散乱光のブリルアンゲインスペクトルを観測するスペクトル観測ステップと、ひずみの観測モデルを構築し、観測したブリルアンゲインスペクトルと観測モデルを用いてひずみを計測するひずみ計測ステップとを含み、観測モデルにおいて、ひずみをフーリエ級数で与える。【選択図】図１
2017058292.tif 000002
【特許請求の範囲】
【請求項１】
環構造体の周囲に配置された光ファイバから、前記環構造体のひずみに基づく散乱光のブリルアンゲインスペクトルを観測するスペクトル観測ステップと、
前記ひずみの観測モデルを構築し、観測した前記ブリルアンゲインスペクトルと前記観測モデルを用いて前記ひずみを計測するひずみ計測ステップとを含み、
前記観測モデルにおいて、前記ひずみをフーリエ級数で与える
変位計測方法。
【請求項２】
前記観測モデルは、以下の（式Ａ）で表される
請求項１に記載の変位計測方法。
【数１】
2017058292.tif 000025
【請求項３】
前記ひずみ計測ステップを行った後、前記ひずみ計測ステップで計測された前記ひずみに基づいて前記環構造体の変位を算出する変位算出ステップを含む
請求項１または２に記載の変位計測方法。
【請求項４】
前記変位算出ステップにおいて、前記環構造体を微小直線梁の集合とみなして前記環構造体の変位を算出する
請求項３に記載の変位計測方法。
【請求項５】
環構造体のひずみを計測するひずみ計測ステップと、
前記ひずみ計測ステップで計測された前記ひずみに基づいて前記環構造体の変位を算出する変位算出ステップとを含み、
前記変位算出ステップにおいて、前記環構造体を微小直線梁の集合とみなして前記環構造体の変位を算出する
変位計測方法。
【請求項６】
請求項１〜５のいずれか１項に記載の変位計測方法をコンピュータに実行させるための
プログラム。
【請求項７】
環構造体の周囲に配置される光ファイバと、
前記光ファイバにおける、前記環構造体のひずみに基づく散乱光のブリルアンゲインスペクトルを観測するスペクトル観測部と、
前記ひずみの観測モデルを構築し、観測した前記ブリルアンゲインスペクトルと前記ひずみの観測モデルを用いて前記ひずみを計測するひずみ計測部とを備え、
前記観測モデルにおいて、前記ひずみはフーリエ級数で与えられる
変位計測装置。
【請求項８】
前記変位計測装置は、前記ひずみに基づいて前記環構造体の変位を算出する変位算出部を備える
請求項７に記載の変位計測装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、構造体のひずみおよび変位を計測する変位計測方法および変位計測装置に関する。
【背景技術】
【０００２】
従来、構造体のひずみ計測システムとして、光ファイバに生じるブリルアンゲインスペクトル（ＢＧＳ）の周波数が発生位置でのひずみに依存して周波数シフトする物理現象を利用した、分布型光ファイバセンシングシステムの一例が提案されている（例えば、特許文献１参照）。
【０００３】
また、例えば特許文献２には、特許文献１に記載された技術を用いて、円環状の断面を有する環構造体について、直径変化などの構造体の変位を計測する技術が開示されている（特許文献２参照）。特許文献２に記載の技術では、環構造体に対する力学的解析に基づいてＢＧＳの形状がモデル化されている。詳細には、ＢＧＳの形状が非ローレンツ関数で与えられている。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２００１−１６５８０８号公報
【特許文献２】特開２０１２−０４７６９９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかし、特許文献２に記載の技術では、上述した解析を行う際に外力の作用点の位置の情報が必要であるため、計測対象物である環構造体が例えばトンネル等であり、それらの情報を得ることができない場合、例えば環構造体の内部に適当な数の支持点で支持された別の環状の構造物をセンサとして設置する。そのため、得られるひずみや変位は、環状の構造物の支持点におけるもののみに限られていた。したがって、計測対象の環構造体の内部における環状の構造物の支持点についてはひずみや変位を直接得ることはできるが、支持点以外の部分のひずみや変位を直接得ることはできず、支持点のひずみや変位から間接的に支持点以外の部分のひずみや変位を得る程度で、支持点以外では、重大なひずみや変位が起こったとしてもそれらを直接検出することができなかった。また、従来の方法では、支持点のような特定の位置のみのひずみや変位の取得しかできなかったため、計測対象物の例えば周方向におけるひずみや変位をともに分布的に取得することができなかった。
【０００６】
上記課題に鑑み、本発明は、計測対象物のあらゆる位置においてひずみおよび変位を計測することができる変位計測方法および変位計測装置を提供することを目的とする。
【課題を解決するための手段】
【０００７】
上記の課題を解決するため、本発明にかかる変位計測方法は、環構造体の周囲に配置された光ファイバから、前記環構造体のひずみに基づく散乱光のブリルアンゲインスペクトルを観測するスペクトル観測ステップと、前記ひずみの観測モデルを構築し、観測した前記ブリルアンゲインスペクトルと前記観測モデルを用いて前記ひずみを計測するひずみ計測ステップとを含み、前記観測モデルにおいて、前記ひずみをフーリエ級数で与える。
【０００８】
この構成によれば、計測対象の環構造体におけるひずみを分布的に検出することができる。したがって、計測対象物である環構造体のあらゆる位置において、ひずみを計測することができる。
【０００９】
また、前記観測モデルは、後に詳述する（式５）で表されてもよい。
【００１０】
この構成によれば、観測モデルである式Ａは、パラメータとしてフーリエ級数で表されるひずみを含むので、周期的に繰り返し観測される、環構造体のひずみを精度よく計測することができる。また、計測対象物である環構造体のあらゆる位置において、ひずみを計測することができる。
【００１１】
また、前記ひずみ計測ステップを行った後、前記ひずみ計測ステップで計測された前記ひずみに基づいて前記環構造体の変位を算出する変位算出ステップを含んでもよい。
【００１２】
この構成によれば、計測したひずみに基づいて環構造体の変位を算出するので、計測対象の環構造体における変位を分布的に算出することができる。したがって、計測対象物である環構造体のあらゆる位置において、変位を算出することができる。
【００１３】
また、前記変位算出ステップにおいては、前記環構造体を微小直線梁の集合とみなして前記環構造体の変位を算出してもよい。
【００１４】
この構成によれば、計測対象の環構造体における変位を分布的に算出することができる。したがって、計測対象物である環構造体のあらゆる位置において、変位を算出することができる。
【００１５】
また、上記の課題を解決するため、本発明にかかる変位計測方法は、環構造体のひずみを計測するひずみ計測ステップと、前記ひずみ計測ステップで計測された前記ひずみに基づいて前記環構造体の変位を算出する変位算出ステップとを含み、前記変位算出ステップにおいて、前記環構造体を微小直線梁の集合とみなして前記環構造体の変位を算出する。
【００１６】
この構成によれば、環構造体の変位の計測精度を向上することができる。
【００１７】
また、上記の課題を解決するため、本発明にかかるプログラムは、上述した特徴を有する変位計測方法をコンピュータに実行させる。
【００１８】
また、上記の課題を解決するため、本発明にかかる変位計測装置は、環構造体の周囲に配置される光ファイバと、前記光ファイバにおける、前記環構造体のひずみに基づく散乱光のブリルアンゲインスペクトルを観測するスペクトル観測部と、前記ひずみの観測モデルを構築し、観測した前記ブリルアンゲインスペクトルと前記ひずみの観測モデルを用いて前記ひずみを計測するひずみ計測部とを備え、前記観測モデルにおいて、前記ひずみはフーリエ級数で与えられる。
【００１９】
この構成によれば、計測対象の環構造体におけるひずみを分布的に直接得ることができる。したがって、計測対象物のあらゆる位置において、ひずみを計測することができる。
【００２０】
また、前記変位計測装置は、前記ひずみに基づいて前記環構造体の変位を算出する変位算出部を備えてもよい。
【００２１】
この構成によれば、計測したひずみに基づいて環構造体の変位を算出するので、計測対象の環構造体における変位を分布的に算出することができる。したがって、計測対象物である環構造体のあらゆる位置において、変位を算出することができる。
【発明の効果】
【００２２】
本発明により、計測対象物のあらゆる位置においてひずみおよび変位を計測することができる変位計測方法および変位計測装置を提供することができる。
【図面の簡単な説明】
【００２３】
【図１】実施の形態にかかる変位計測装置の構成を示す概略図
【図２】実施の形態にかかる変位計測装置について、環構造体へ光ファイバを設置した一例を示す概略図
【図３】実施の形態にかかる変位計測装置について、環構造体へ光ファイバを設置した一例を示す概略図
【図４】実施の形態にかかる変位計測方法の手順を示すフローチャート
【図５】実施の形態にかかる変位計測方法における環構造体の観測モデルを、環構造体の長手方向に見た図
【図６】実施の形態にかかるひずみ計測の手順を示すフローチャート
【図７】実施の形態にかかるひずみ計測について、評価関数Ｅの計算を模式的に示した図
【図８】実施の形態にかかる変位計測方法において環構造体の変位を算出するための観測モデルの図
【図９】微小直線梁を用いた直線梁の観測モデルの図
【図１０】微小直線梁におけるたわみ角の微小変化量について説明するための図
【図１１】実施の形態にかかる変位計測方法における変形前の環構造体の観測モデルを、環構造体の長手方向に見た図
【図１２】微小直線梁を用いた環構造体の観測モデルの図
【図１３】実施の形態にかかる変位計測方法による環構造体の外周のひずみのシミュレーション結果を示す図
【図１４】実施の形態にかかる変位計測方法による環構造体の内周のひずみのシミュレーション結果を示す図
【図１５】環構造体の変形前後の状態を示す概略図
【図１６】実施の形態にかかる変位計測方法による環構造体の変位のシミュレーション結果を示す図
【発明を実施するための形態】
【００２４】
以下、図面を用いて、本発明にかかる実施の形態について説明する。なお、図面において、同一の符号が付された構成要素は、同一または同種の構成要素を示す。
【００２５】
また、以下で説明する実施の形態は、本発明の好ましい一具体例を示す。以下の実施の形態で示される数値、形状、材料、構成要素、構成要素の配置位置、接続形態、ステップおよびステップの順序等は、一例であり、本発明を限定する主旨ではない。また、以下の実施の形態における構成要素のうち、本発明の最上位概念を示す独立請求項に記載されていない構成要素については、より望ましい形態を構成する任意の構成要素として説明される。
【００２６】
（実施の形態）
本実施の形態にかかる変位計測方法および変位計測装置について説明する。本実施の形態にかかる変位計測方法は、変位計測装置に含まれる光ファイバをトンネルなどの環構造体の円環の中立軸以外の位置に配置し、光ファイバに発生するブリルアン散乱光を利用することにより、環構造体のひずみを計測する方法および計測されたひずみから変位を計測する方法からなっている。後者のひずみから変位を計測するにあたっては、前者の方法に限定されず、例えばひずみゲージなど別の計測方法によって得られたひずみを用いてもよい。
【００２７】
光ファイバに光が入射されると、光の散乱現象によって光ファイバ内のあらゆる位置でレーリー散乱光、ラマン散乱光、ブリルアン散乱光などの散乱光が発生する。ブリルアン散乱光とは、光が物質中で音波と相互作用して発生した光であり、その周波数は、入射光の周波数よりわずかにずれている。ブリルアン散乱光のスペクトルであるブリルアンゲインスペクトル（ＢＧＳ）は、光ファイバの長手方向に生じるひずみに比例して周波数シフトするという性質を持つ。そこで、ＢＧＳの周波数変化を利用して、環構造体のひずみを計測する。また、計測した環構造体のひずみから環構造体の変位を算出する。以下、ＢＧＳを利用した、本実施の形態にかかる変位計測方法および変位計測装置について説明する。
【００２８】
［１．変位計測装置の構成］
はじめに、本実施の形態にかかる変位計測装置の構成について、図１〜図３を用いて説明する。図１は、本実施の形態にかかる変位計測装置の構成を示す概略図である。図２および図３は、本実施の形態にかかる変位計測装置について、環構造体へ光ファイバを設置した一例を示す概略図である。
【００２９】
図１に示すように、本実施の形態にかかる変位計測装置１００は、環状の計測対象物（環構造体）１０の外周および内周にそれぞれ配置される光ファイバ２０ａおよび２０ｂと、スペクトル観測部３０と、ひずみ計測部４０と、変位算出部５０とを備えている。
【００３０】
環構造体１０は、例えばトンネルであり、図２および図３に示すように、円柱状の構造物の内部に長手方向に沿った空洞を有している。図３の（ａ）に示すように、環構造体１０を長手方向と直交する平面で切断したときの断面は、環状、例えば円環となるように形成されている。
【００３１】
また、図２および図３の（ｂ）に示すように、光ファイバ２０ａは、環構造体１０の外周に配置され、光ファイバ２０ｂは環構造体１０の内周に配置される。
【００３２】
光ファイバ２０ａおよび２０ｂは、環構造体１０のひずみを測定するためのセンシング用光ファイバである。光ファイバ２０ａは、光ファイバ２０ａおよび２０ｂが巻き付けられた環構造体１０の外周を含む断面の円環が、環構造体１０の長手方向の軸にほぼ直交するように固定される。光ファイバ２０ｂは、光ファイバ２０ａが巻き付けられた環構造体１０の外周と同一の位置の内周に固定される。光ファイバ２０ａ、２０ｂの設置は、外周、内周に限られるものではなく、円環の中立軸以外の内部であってもよい。ただし、設置場所にかかわらず、これらの光ファイバ２０ａおよび２０ｂには、固定位置のひずみと同じひずみが生じているものとする。
【００３３】
なお、光ファイバ２０ａおよび２０ｂは、環構造体１０のひずみに代えて変位を算出することとしてもよい。環構造体１０のひずみの算出方法については、後に詳述する。また、光ファイバ２０ａおよび２０ｂが環構造体１０に固定される位置は一箇所の断面に限らず、複数箇所であってもよい。この場合、環構造体１０の複数の箇所に固定された光ファイバ２０ａは、一本の光ファイバが一筆書きになるように接続された構成であってもよい。また、光ファイバ２０ｂについても同様である。
【００３４】
また、光ファイバ２０ａおよび２０ｂの一端は、スペクトル観測部３０に接続されている。
【００３５】
スペクトル観測部３０は、光ファイバ２０ａおよび２０ｂに所定の周波数の光を入射し、光ファイバ２０ａおよび２０ｂの内部で散乱したＢＧＳをそれぞれ観測する。スペクトル観測部３０は、観測したＢＧＳをひずみ計測部４０に出力する。
【００３６】
ひずみ計測部４０は、上述したスペクトル観測部３０で観測されたＢＧＳから環構造体１０で発生するひずみを計測するものである。
【００３７】
変位算出部５０は、ひずみ計測部４０で計測された環構造体１０のひずみに基づいて変位を算出するものである。ただし、前述したように、変位算出部５０に入力されるひずみは、スペクトル観測部３０とひずみ計測部４０より得られたひずみだけではなく、例えばひずみゲージなど別な方法によって得られたものであってもよい。算出方法については後に詳述する。
【００３８】
なお、スペクトル観測部３０、ひずみ計測部４０および変位算出部５０は、例えばコンピュータで構成されている。また、コンピュータは、以下に説明する変位計測方法を実行するためのプログラムを備え、当該プログラムにより変位計測方法を行ってもよい。
【００３９】
［２．変位計測方法］
次に、上述した変位計測装置１００を用いた変位計測方法について説明する。図４は、変位計測方法の手順を示すフローチャートである。
【００４０】
図４に示すように、変位計測方法においては、はじめに、スペクトル観測部３０により、環構造体１０に設置された光ファイバ２０ａおよび２０ｂ（図３の（ａ）参照）のＢＧＳを観測する（ステップＳ１１）。これにより、図３の（ａ）に示した環構造体１０の外周と内周におけるＢＧＳが得られる。なお、ステップＳ１１は、本発明におけるスペクトル観測ステップに相当する。
【００４１】
次に、ひずみ計測部４０において、予め構築されていた観測モデル（式）が、観測したＢＧＳに近づくようにあてはめを行う。これにより、環構造体１０の外周および内周のひずみを計測する（ステップＳ１２）。なお、ステップＳ１２は、本発明におけるひずみ計測ステップに相当する。ステップＳ１１およびステップＳ１２の代わりに、別な方法によって、環構造体の外周および内周のひずみを計測するステップを入れてもよい。
【００４２】
さらに、変位算出部５０において、計測された環構造体１０の外周および内周のひずみから、環構造体１０の各位置における変位を算出する（ステップＳ１３）。なお、ステップＳ１３は、本発明における変位算出ステップに相当する。
【００４３】
以下、各ステップについて詳細に説明する。
【００４４】
［２−１．環構造体のひずみ計測の原理］
はじめに、ひずみ計測の原理について説明する。図５は、本実施の形態にかかる変位計測方法における環構造体１０の観測モデルを、環構造体１０の長手方向に見た図である。
【００４５】
ひずみ計測では、上述したように、予めＢＧＳのモデル（モデルＢＧＳ）を構築しておき、ＢＧＳを観測して、予め構築されていたモデルＢＧＳが観測したＢＧＳに近づくようにあてはめを行うものである。また、ＢＧＳの観測は、環構造体１０の断面ごとに行うものとする。
【００４６】
まず、計測対象の環構造体１０のモデルＢＧＳについて説明する。ここでは、環構造体１０として、環構造体１０の長手方向と垂直な方向の断面が円環である環構造体を観測モデルとする。なお、環構造体１０の断面は円環に限らず、非線形の連続する構造であれば、どのようなものであってもよい。例えば、断面が楕円形の環構造体であってもよい。
【００４７】
図５に示すように、環構造体１０の長手方向（ζ方向）の軸に垂直な面において、環構造体１０の断面である円環の中心を原点Ｏ、原点Ｏを通る一方向にｘ軸、原点Ｏを通りｘ軸方向に垂直な他の方向にｙ軸をとる。さらに、ｘ軸の負方向から右回りの角度をφとする。
【００４８】
ここで、環構造体１０の周囲に設置された光ファイバ２０ａおよび２０ｂ上（図５には示されていない）での、第ｉ番目のＢＧＳ観測点に対応する角度をφｉ（ｉ＝１，２，…，Ｉ、ＩはＢＧＳの観測点の数）、φｉを中心としたＢＧＳの観測区間をΔφとし、ＢＧＳを観測する周波数については、第ｊ番目の観測周波数をνｊ（ｊ=１，２，…，Ｊ、Ｊは観測周波数の数）とする。また、観測角度φｉ、観測周波数νｊで観測されたＢＧＳの値をＧｏｂｓ（νｊ，φｉ）、モデルＢＧＳより求められたＢＧＳの値をＧｍｏｄ（νｊ，φｉ，ｑ）とする。そして、このモデルＢＧＳが観測されたＢＧＳに最も良く一致する条件から、求めるべきパラメータｑを得る。具体的には、全観測角度、全観測周波数について、観測ＢＧＳとモデルＢＧＳの差の二乗和を評価関数Ｅとして次の（式１）で与え、評価関数Ｅが最小となるときのパラメータｑを非線形最小二乗法を用いて求める。なお、パラメータｑは、フーリエ係数で表されるパラメータｐを含むパラメータである。なお、パラメータｐ、ｑはベクトルであり、以下に示す式中では太字で示している。パラメータｐおよびｑについては後に詳述する。
【００４９】
【数１】
2017058292.tif 000003
………………………………………………………………………………………………（式１）
【００５０】
環構造体１０の断面における円環の円周上の観測角度φにおけるひずみε（φ）はε（φ＋２π）に等しく、ε（φ）は周期２πをもつ周期関数となっている。そこで、環構造体１０に生じるひずみを、周期２πのフーリエ級数として次の（式２）で与える。
【００５１】
【数２】
2017058292.tif 000004
………………………………………………………………………………………………（式２）
【００５２】
ここで、（式２）のフーリエ係数をまとめて以下の（式３）のように表す。
【００５３】
【数３】
2017058292.tif 000005
………………………………………………………………………………………………（式３）
【００５４】
ここで、Ｎはフーリエ係数の最高次数である。
【００５５】
これらのパラメータを明示的に示す場合にはひずみをε（φ，ｐ）と書くことにする。ε（φ，ｐ）に対応する周波数は、次の（式４）で表される。
【００５６】
【数４】
2017058292.tif 000006
【００５７】
環構造体１０の断面の円環の半径を１に規格化すると、φｉ、νｊにおけるモデルＢＧＳは次の（式５）で表される。
【００５８】
【数５】
2017058292.tif 000007
【００５９】
なお、（式５）は、本発明における（式Ａ）に相当する。νＢ（０）はひずみ０に対応する周波数、Ｃはひずみと光周波数の関係を与える比例係数である。
【００６０】
（式５）において、ｐは求めるべきひずみを与えるパラメータであり、ｈ、ｗ、νＢ（０）、ＣはＢＧＳの基本特性についてのパラメータである。
【００６１】
上述したように、環構造体１０の断面ごとにＢＧＳが観測され、ひとつの断面に設置された光ファイバ２０ａおよび２０ｂの長さは光損失が十分小さく無視できる程度であると考えられることから、環構造体１０の断面ごとの光ファイバ２０ａおよび２０ｂではｈは一定であると仮定することができる。しかし、環構造体１０の断面間ではその距離が長いこと、および、環構造体１０の断面間に光ファイバ２０ａおよび２０ｂの接続点がある場合があるなどの理由により、断面間では光ファイバ２０ａおよび２０ｂに生じる散乱光のパワーが異なることが考えられる。そのため、ｈについても環構造体１０の断面ごとに観測されたＢＧＳから推定する必要がある。したがって、フーリエ係数で表されるパラメータｐと各断面の環構造体１０におけるｈの両方を求めることとする。以下では、両パラメータを合わせてｑ=［ｐ  ｈ］と表す。
【００６２】
［２−２．ひずみの計測］
次に、上述したひずみ計測方法により、具体的に環構造体１０のひずみの計測を行う。図６は、本実施の形態にかかるひずみ計測の手順を示すフローチャートである。図７は、本実施の形態にかかるひずみ計測について、評価関数Ｅの計算を模式的に示した図である。
【００６３】
はじめに、上述したようにＢＳＧの観測を行う。ＢＧＳの観測では、まず、図６に示すように、ＢＧＳの観測条件を設定する（ステップＳ２１）。観測条件として、観測点と、観測周波数とを定める。
【００６４】
ＢＧＳを利用したひずみ計測では、ひずみは環構造体１０上の位置の関数として得られる。そこで、環構造体１０の長手方向である、例えばζ方向において、環構造体１０の第ｋ番目の断面である円環の円周上のＢＧＳ観測点の位置を、環構造体１０に設置されている光ファイバ２０ａおよび２０ｂの半径を用いて、その位置に対応する角度φｉに変換する。
【００６５】
以下では、簡単化のためＢＧＳ観測点の位置が等間隔である場合について考える。また、どの断面についても同様であることから、断面の番号については省略する。観測開始角度をφ１、全観測点の数をＩ、観測点の間隔をδφで表すと、φｉは次の（式６）で表される。
【００６６】
【数６】
2017058292.tif 000008
………………………………………………………………………………………………（式６）
【００６７】
これにより、ＢＧＳの観測点の位置が角度φｉで定められる。
【００６８】
また、観測開始周波数をν１、全観測周波数の数をＪとし、周波数についても等間隔δνで観測するものとすると、観測周波数νｊは次の（式７）で表される。
【００６９】
【数７】
2017058292.tif 000009
………………………………………………………………………………………………（式７）
【００７０】
次に、上述した観測条件の下でＢＧＳを観測し、観測周波数と角度をパラメータとしたＧｏｂｓ（νｊ，φｉ）（ｉ=１，２，…，Ｉ）（ｊ=１，２，…，Ｊ）を得る（ステップＳ２２）。Ｇｏｂｓ（νｊ，φｉ）は、観測周波数νｊ、角度φｉにおいて観測されたＢＧＳの値である。
【００７１】
次に、すべてのφｉおよびνｊについてモデルＢＧＳの値Ｇｍｏｄ（νｊ，φｉ，ｑ）を算出する（ステップＳ２３）。ここで、Ｇｍｏｄ（νｊ，φｉ，ｑ）の値は、（式５）の数学的な積分、または、数値積分によって求める。また、評価関数Ｅを最小化するｑを効率的に得るために、例えば、滑降シンプレックス法を用いる。滑降シンプレックス法は、求めるべき解ｑ=［ｐ  ｈ］の候補群を、求めるべきパラメータの個数より１個以上多く用意し、各解候補に対応する評価値が小さくなるように解候補群を更新していくことで最適解を得るという繰り返し計算方法である。
【００７２】
次に、算出したモデルＢＧＳと観測されたＢＧＳを用いて（式１）の評価関数Ｅを算出する。図７は、上述した（式１）における評価関数Ｅの計算を模式的に示した図であり、（ａ）は観測されたＢＧＳであるＧｏｂｓ（νｊ，φｉ）、（ｂ）はモデルより求められたＢＧＳであるＧｍｏｄ（νｊ，φｉ，ｑ）である。図７に示すように、モデルより求められたＢＧＳの値Ｇｍｏｄ（νｊ，φｉ，ｑ）（図７の（ｂ））を、観測されたＢＧＳの値Ｇｏｂｓ（νｊ，φｉ）（図７の（ａ））に近づける。より具体的には、滑降シンプレックス法を用いてパラメータｑを更新しながら、モデルＢＧＳと観測されたＢＧＳとの差を成分としてもつ評価関数Ｅを最小化するためのパラメータｑを求める（ステップＳ２４）。そして、求められたパラメータｑのうちのひずみにかかわるｐを（式２）に代入することによって、環構造体１０の外周側のひずみεｏｕｔ、および、内周側のひずみεｉｎが得られる。
【００７３】
なお、ひずみの計測方法は、上述したＢＧＳを利用した方法に限らず、他の方法であってもよい。他の方法の場合でも、ひずみの観測モデルを構築し、構築した観測モデルを用いて環構造体のひずみを計測すればよい。例えば、ＢＧＳ観測点に相当する位置にひずみゲージなどを設置し、それを用いてひずみを計測し、それらのひずみをεｏｕｔおよびεｉｎとしてもよい。
【００７４】
［２−３．変位の算出］
次に、環構造体１０の変位の算出を行う手順について説明する。
【００７５】
図８は、本実施の形態にかかる変位計測方法において環構造体１０の変位を算出するための観測モデルの図であり、（ａ）は観測モデルを模式的に示した図、（ｂ）は変位前後の観測位置の関係を示す図である。
【００７６】
環構造体１０の解析は複雑であるため、環構造体１０を後に詳述する微小直線梁（微小片）の集合とみなして環構造体１０の変位を算出する。詳細には、図８の（ａ）に示すように、環構造体１０は微小片６０の集合体であるとみなす。そして、図８の（ｂ）に示すように、環構造体１０の変形前の点Ｑに配置されていた微小片６０が、環構造体１０の変形に伴い点Ｑ’に移動したものとして、環構造体１０の変位を算出する。
【００７７】
ここでは、簡単化のために、図８の（ｂ）に示すように、環構造体１０の円環の中立軸は円環の中心線に一致しており、その半径をｒとする。また、この中心線から外周と内周までの距離は等しく、その距離をhとする。図８の（ｂ）では、図を見やすくするために外周と内周の一部のみが示されている。これにより、内周の半径は（ｒ−η）、外周の半径は（ｒ＋η）と表される。また、点Ｑの位置は、環構造体１０の外周と内周との中間を通る円環とｘ軸の負方向との交点である基準点Pから時計回りに角度φｍ回転させた位置とする。このとき、点Ｑから点Ｑ’へのｘ軸方向の変位はΔｘｍ、ｙ軸方向の変位はΔｙｍと表す。なお、角度φｍは、環構造体１０における第ｍ番目の変位算出角度とし、ｍを変更することにより環構造体１０上の各位置の変位を得ることができる。
【００７８】
以下、上述した、円環のひずみから環構造体１０の変位を算出する方法について説明する。
【００７９】
はじめに、理解を容易にするために、微小片（微小直線梁）が直線状に配置された直線梁を観測モデルとした従来の変位算出方法について説明する。図９は、微小直線梁を用いた直線梁の観測モデルの図である。
【００８０】
図９に示す観測モデルは、複数の微小直線梁の集合体である直線梁６１である。直線梁６１において、左端から第１番目の微小直線梁を微小直線梁６０ａ、第２番目の微小直線梁を微小直線梁６０ｂ、第ｍ番目の微小直線梁を微小直線梁６０ｃと表している。図９は、直線梁６１の左端を固定し、右端に変位を与えたときの変形後の微小直線梁６０ａ、微小直線梁６０ｂおよび微小直線梁６０ｃの状態を示している。変形前の直線梁６１の長手方向にｕ軸、変形前の直線梁６１の長手方向と直交する方向にν軸をとり、微小直線梁６０ａ、微小直線梁６０ｂおよび微小直線梁６０ｃのそれぞれの長さをｄｕとしている。この長さは変形前後で変わらない。
【００８１】
ここで、微小直線梁６０ｃのたわみ角、すなわち、微小直線梁６０ｃとｕ軸とのなす角度をωｍ、第ｍ−１番目の微小直線梁に対する微小直線梁６０ｃのたわみ角の微小変化量をｄωｍ、微小直線梁６０ｃのν軸方向の全変位をνｍ、第ｍ−１番目の微小直線梁に対する微小直線梁６０ｃの変位をｄνｍとし、これらの関係を求める。
【００８２】
ｍ＝１すなわち第１番目の微小直線梁６０ａでは、図９に示すように、たわみ角ω１はその微小変化量ｄω１となっている。また、その変位ｄν１も全変位ν１に等しい。したがって、第１番目の微小直線梁６０ａのν軸方向の全変位ν１は、次の（式８）のように表される。
【００８３】
【数８】
2017058292.tif 000010
……………………………………………………………………………………………（式８）
【００８４】
また、微小直線梁６０ｂのたわみ角ω２は、微小直線梁６０ａと微小直線梁６０ｂに生じた微小変化量ｄω１とｄω２の和であり、全変位ν２もそれぞれの変位ｄν１とｄν２の和となっている。したがって、第２番目の微小直線梁６０ｂの変位ｄν２は、次の（式９）のように表される。
【００８５】
【数９】
2017058292.tif 000011
……………………………………………………………………………………………（式９）
【００８６】
したがって、第２番目の微小直線梁６０ｂのν軸方向の全変位ν２は、次の（式１０）のように表すことができる。
【００８７】
【数１０】
2017058292.tif 000012
…………………………………………………………………………………………（式１０）
【００８８】
以下、同様に考えて、第２番目を第ｍ番目に置き換えると、第ｍ番目の微小直線梁６０ｃのたわみ角ωｍ、ν軸方向の全変位νｍ、第ｍ−１番目の微小直線梁に対する微小直線梁６０ｃの変位ｄνｍ、ν軸方向の全変位をνｍは、
【数１１】
2017058292.tif 000013
……………………………………………………………………………………………（式１１）
【数１２】
2017058292.tif 000014
……………………………………………………………………………………………（式１２）
【数１３】
2017058292.tif 000015
……………………………………………………………………………………………（式１３）
【数１４】
2017058292.tif 000016
……………………………………………………………………………………………（式１４）
と表すことができる。
【００８９】
次に、上述した直線梁６１の観測モデルにおいて、計測されたひずみから微小片（微小直線梁）６０のたわみ角の微小変化量ｄωを求める。図１０は、微小片６０におけるたわみ角の微小変化量ｄωについて説明するための図である。
【００９０】
図１０では、Ａ、Ｂ、Ｃ、Ｄを頂点とする微小直線梁６０ｄを示している。微小直線梁６０ｄの厚み、すなわち、図１０に示す辺ＡＤおよび辺ＢＣの長さを２ηとする。また、外周側の辺ＡＤおよび内周側の辺ＢＣに生じているひずみをそれぞれεｏｕｔ、εｉｎとする。また、微小直線梁６０ｄは、頂点Ａ、Ｂ、Ｃ、Ｄがそれぞれ頂点Ａ’、Ｂ’、Ｃ、Ｄの位置となるように変形するものとする。微小直線梁６０ｄにおいて、たわみ角の微小変化量ｄωが生じたとすると、Ａ’Ｄ=（１＋εｏｕｔ）ｄｕ、Ｂ’Ｃ=（１＋εｉｎ）ｄｕとなるので、たわみ角の微小変化量ｄωは次の（式１５）で求められる。
【００９１】
【数１５】
2017058292.tif 000017
……………………………………………………………………………………………（式１５）
【００９２】
なお、（式１５）の左辺においてｄωに負の符号がついているのは、図１０では右回りの角度を正にとっているため、左回りのたわみ角の微小変化量ｄωは負で示されるためである。
【００９３】
（式１４）と（式１５）より、ν軸方向の全変位νｍは、次の（式１６）のように書き換えられる。
【００９４】
【数１６】
2017058292.tif 000018
……………………………………………………………………………………………（式１６）
【００９５】
次に、上述した直線梁６１の観測モデルにおける変位計測方法を環構造体１０に応用した、微小直線梁を用いた環構造体１０の観測モデルについて説明する。
【００９６】
図１１は、本実施の形態にかかる変位計測方法における変形前の環構造体１０の観測モデルを、環構造体１０の長手方向に見た図である。図１２は、微小直線梁を用いた環構造体１０の観測モデルの図である。
【００９７】
図１１に示すように、変形前の環構造体１０の断面の円環において、頂点をＡ、Ｂ、Ｃ、Ｄとする微小片７０を考える。この微小片７０に対応する微小中心角をΔθとする。また、図８の（ｂ）と同様、環構造体１０の円環の中立軸は円環の中心線に一致しており、その半径をｒとする。また、この中心線から外周と内周までの距離を、いずれもηとする。このとき、環構造体１０の外周における弧ＡＤと内周における弧ＢＣの長さは、それぞれ（ｒ＋η）Δθ、（ｒ−η）Δθとなる。
【００９８】
ここで、ｒ≫ηとすると、両者はいずれもｒΔθと近似することができるので、この微小片７０を微小直線梁とみなすことができる。そこで、図１２に示すように、環構造体１０を微小片７０の集合として考える。図１２に示す一点鎖線は、変形前の環構造体１０の外周と内周との中央を通る円環の中心線である。変位算出に用いる角度θを、φと同じようにｘ軸負方向から右回りにとる。
【００９９】
図１２において、環構造体１０は、第１番目の微小片７０ａ、第２番目の微小片７０ｂ、・・・、第ｍ番目の微小片（図示せず）を有している。環構造体１０がＭ個の微小片で構成されている場合、第Ｍ番目の微小片と第１番目の微小片７０ａとは接続している。
【０１００】
ここで、第ｍ番目の微小片のν軸とｘ軸とのなす角をθｍ、ｘ軸方向およびｙ軸方向の全変位をそれぞれΔｘｍおよびΔｙｍ、第ｍ−１番目の微小片に対する第ｍ番目の微小片の変位をｄｘｍ、ｄｙｍとし、これらの関係を求める。上述した直線梁６１の観測モデルと同様に、ｍ＝１すなわち第１番目の微小片７０ａからｍ＝２すなわち第２番目の微小片７０ｂについて順に求め、第２番目を第ｍ番目に置き換えると、第ｍ番目の微小片のｘ軸方向の全変位Δｘｍおよびｙ軸方向の全変位Δｙｍは、それぞれ以下の（式１７）、（式１８）のように表される。
【０１０１】
【数１７】
2017058292.tif 000019
……………………………………………………………………………………………（式１７）
【０１０２】
【数１８】
2017058292.tif 000020
……………………………………………………………………………………………（式１８）
【０１０３】
円環の場合にはｄｕ=ｒΔθであり、また上述した微小直線梁６０ｃで得られている結果、すなわち、（式１３）を上式に代入すると、ｘｍおよびｙｍは以下の（式１９）および（式２０）ように表される。
【０１０４】
【数１９】
2017058292.tif 000021
……………………………………………………………………………………………（式１９）
【０１０５】
【数２０】
2017058292.tif 000022
……………………………………………………………………………………………（式２０）
【０１０６】
なお、εｏｕｔ，ｍ、εｉｎ，ｍは、それぞれ第ｍ番目の微小片の外周および内周に生じているひずみである。
【０１０７】
上述したように、適当な角度の間隔で外周側の辺ＡＤおよび内周側の辺ＢＣにおけるひずみεｉｎ、εｏｕｔを求め、（式１９）および（式２０）に代入することによって、環構造体１０の任意の位置の変位を得ることができる。
【０１０８】
以上の式を用いることにより、ひずみεｉｎ、εｏｕｔから環構造体１０のｘ軸方向の変位Δｘｍおよびｙ軸方向の変位Δｙｍを得ることができる。
【０１０９】
［３．効果等］
以下、上述した変位計測方法の有効性を確かめるために行ったシミュレーションの結果について説明する。図１３は、実施の形態にかかる変位計測方法による環構造体の外周のひずみのシミュレーション結果を示す図である。図１４は、実施の形態にかかる変位計測方法による環構造体の内周のひずみのシミュレーション結果を示す図である。図１５は、環構造体の変形前後の状態を説明するための概略図である。図１６は、本実施の形態にかかる変位計測方法による環構造体の変位のシミュレーション結果を示す図である。
【０１１０】
はじめに、ひずみ計測のシミュレーション結果について説明する。ひずみ計測のシミュレーションでは、以下のようにして観測ＢＧＳであるＧｏｂｓ（νｊ，φｉ）を生成する。生成したＧｏｂｓ（νｊ，φｉ）によりひずみを求めるシミュレーションを行った。
【０１１１】
ここで、環構造体１０の半径をｒ、環構造体１０の中立軸から環構造体１０の外周までの距離をη、断面積をＡ、断面係数をκ、力学的パラメータであるヤング率をＥＹとし、後に説明する図１５の（ａ）に示すように、環構造体１０のｙ軸に一対の荷重が作用している場合を考えると、環構造体１０の円周上に生じるひずみε（φ）は以下の（式２１）で近似される。
【０１１２】
【数２１】
2017058292.tif 000023
……………………………………………………………………………………………（式２１）
【０１１３】
（式２１）に示したひずみε（φ）を用い、設定されたＢＧＳ観測角度φｉ、観測周波数νｊについて算出されたＢＧＳを観測ＢＧＳ  Ｇｏｂｓ（νｊ，φｉ）とした。シミュレーションでは、（式２１）におけるｒ、η、κＥＹＡの値を、９７１．４［ｍｍ］、２３．６［ｍｍ］、５３６１×１０９［Ｎ・ｍｍ３］とした。ＢＧＳの基本特性であるｈ、νＢ（０）、ｗ、Ｃの値は、５．１０×１０−９［Ｗ］、１０．８２［ＧＨｚ］、０．０３［ＧＨｚ］、５０．０［ＧＨｚ］とした。これらの値は、観測に用いる環構造体１０の一例であり、光ファイバの典型的な値である。なお、νＢ（０）はひずみ０に対応する周波数、Ｃはひずみと光周波数の関係を与える比例係数である。さらに、ＢＧＳの観測条件にかかわるφ１、Ｉ、ν１、δν、Ｊはそれぞれ、０［°］、１３０、１０．８２［ＧＨｚ］、０．００２［ＧＨｚ］、１６１である。また、Δφおよびδφは、環構造体１０の内外周で異なっており、外周に対するΔφおよびδφは、それぞれ２．９［°］および５７．６［°］、内周に対するΔφおよびδφは、それぞれ３．０［°］および６０．５［°］とした。
【０１１４】
ここでは、Ｐを７０［Ｎ］、フーリエ級数の最大次数Ｎを１０としてひずみε（φ）を求めた。求められた環構造体１０の外周と内周のひずみをそれぞれ図１３および図１４に示す。実線は上述した本実施の形態におけるひずみ計測方法によって計測されたひずみ、細線は（式２１）で算出された理論的に求められたひずみであり、破線はＢＧＳをローレンツ関数にあてはめる従来の方法で得られたひずみである。外周および内周のいずれのひずみも、本実施の形態におけるひずみ計測方法により求められたひずみのほうが、従来方法により求められたひずみよりも、理論的に求められたひずみに近いことがわかる。
【０１１５】
９０［°］において、理論的ひずみ、従来の方法で得られた外周のひずみ、および、本実施の形態におけるひずみ計測方法で得られた外周のひずみはそれぞれ９．９４×１０−５、６．２０×１０−５、９．０５×１０−５であり、理論的ひずみと従来の方法で得られた外周のひずみの誤差、および、理論的ひずみと本実施の形態におけるひずみ計測方法で得られた外周のひずみの誤差は、それぞれ−３７．４×１０−６、−８．９×１０−６であった。また、１８０［°］においては、理論的ひずみ、従来の方法で得られた外周のひずみ、および、本実施の形態におけるひずみ計測方法で得られた外周のひずみはそれぞれ−５．５５×１０−５、−５．００×１０−５、−５．６７×１０−５であり、理論的ひずみと従来の方法で得られた外周のひずみの誤差、および、理論的ひずみと本実施の形態におけるひずみ計測方法で得られた外周のひずみの誤差は、それぞれ５．５×１０−６、−１．２×１０−６であった。
【０１１６】
また、同様に、９０［°］において、理論的ひずみ、従来の方法で得られた外周のひずみ、および、本実施の形態におけるひずみ計測方法で得られた外周のひずみはそれぞれ−１０．２７×１０−５、−６．００×１０−５、−９．３５×１０−５であり、理論的ひずみと従来方法で得られた内周のひずみの誤差、および、理論的ひずみと本実施の形態におけるひずみ計測方法で得られた内周のひずみの誤差は、それぞれ４２．７×１０−６、９．２×１０−６であった。また、１８０［°］においては、理論的ひずみ、従来の方法で得られた外周のひずみ、および、本実施の形態におけるひずみ計測方法で得られた外周のひずみはそれぞれ５．９９×１０−５、５．２０×１０−５、６．１０×１０−５であり、理論的ひずみと従来方法で得られた内周のひずみの誤差、および、理論的ひずみと本実施の形態におけるひずみ計測方法で得られた内周のひずみの誤差は、それぞれ−７．９×１０−６、１．１×１０−６であった。
【０１１７】
以上より、本実施の形態におけるひずみ計測方法では、誤差が数分の一に低減されていることがわかった。このように、本実施の形態におけるひずみ計測方法では、計測精度が向上していることが確認できた。
【０１１８】
次に、変位計測のシミュレーション結果について説明する。
【０１１９】
上述した環構造体１０の変位計測方法の有効性を確認するために、環構造体１０の変位計測シミュレーションを行った。ここでは、上述した本実施の形態におけるひずみ計測方法と同じ条件の下、環構造体１０の内周のひずみを与えるフーリエ係数ｐｉｎと外周のひずみを与えるフーリエ係数ｐｏｕｔとが得られているものとする。ひずみの計算角度間隔Δθを０．１［°］として、得られているフーリエ係数ｐｉｎおよびｐｏｕｔを（式２）に代入して円環の内周と外周のひずみを算出し、変位を求めた。
【０１２０】
まず、図１５の（ｂ）に示すように、ｙ軸とｘ軸に関して対象に変位が発生することから、環構造体１０を長手方向の軸（中心Ｏ）で固定したときの任意の位置における変位を求める。ここで、破線、実線の円はそれぞれ変形前、変形後の円環の中心線、破線円から実線円に伸びている矢印は変位である（各位置から円の外側方向の変位を正としている）。
【０１２１】
求められた環構造体１０の変位を図１６に示す。図１６における実線は、観測されたＢＧＳから環構造体１０の内周のひずみおよび外周のひずみを求めて変位を計測する本実施の形態における変位計測方法によって求めた変位、破線は（式２１）から理論的ひずみを求め、（式１９）および（式２０）により求めた変位である。ここで、ｘ、ｙ軸方向の直径変化Ｄｘ、Ｄｙは力学的解析によって次の（式２２）で与えられる。
【０１２２】
【数２２】
2017058292.tif 000024
……………………………………………………………………………………………（式２２）
【０１２３】
９０、１８０°における求めた円環の変位は、それぞれＤｙ／２、Ｄｘ／２の値に対応する。まず、（式２２）より求められた９０［°］におけるＤｙ／２と理論的ひずみから求めた変位はそれぞれ９．４３６×１０−１［ｍｍ］、９．４６２×１０−１［ｍｍ］であり、その誤差は２．６×１０−３［ｍｍ］、１８０［°］におけるＤｘ／２と理論的ひずみから求めた変位はそれぞれ−８．６６５×１０−１［ｍｍ］、−８．６３５×１０−１［ｍｍ］であり、その誤差は３．０×１０−３［ｍｍ］であることから、誤差は元の変位の１％より小さくなっていることがわかった。したがって、円環の内外周のひずみから変位を計測する方法の有効性が確認できた。
【０１２４】
次に、９０［°］における本方法で得られた変位は９．４５３×１０−１［ｍｍ］であり、上述のＤｙ／２との誤差は１．７×１０−３［ｍｍ］、１８０［°］において本方法で得られた変位は−８．６４３×１０−１［ｍｍ］であり、Ｄｘ／２との誤差は２．２×１０−３［ｍｍ］であった。このように、本実施の形態における変位計測方法による変位の誤差は、理論的ひずみから求めた変位の誤差と同程度小さくなり、計測精度が向上していることがわかった。したがって、本実施の形態における環構造体１０の変位計測方法の有効性が確認できた。
【０１２５】
以上、本実施の形態にかかる変位計測装置によると、観測モデルにおいてひずみをフーリエ級数で与えており、計測したＢＧＳから、観測モデルを用いてひずみを計測するので、環構造体においてセンサである環状の構造物を支持している点に限らず、支持している点以外の部分においてもひずみを直接得ることができる。また、計測したひずみに基づいて環構造体の変位を算出するので、環構造体における変位を分布的に算出することができる。したがって、計測対象物である環構造体のあらゆる位置において、ひずみおよび変位を算出することができる。
【０１２６】
これにより、計測対象物のあらゆる位置においてひずみおよび変位を計測することができる。
【０１２７】
なお、上述した実施の形態では、スペクトル観測ステップの後、ひずみ計測を行い、計測で得られたひずみに基づいて変位を算出する構成であったが、ひずみ計測を行わず、スペクトル観測ステップの後ひずみ計測を行わず、観測したＢＧＳに基づいて変位を算出する構成であっってもよい。
【０１２８】
以上、本発明にかかる変位計測方法および変位計測装置について、実施の形態に基づいて説明したが、本発明は実施の形態に限定されるものではない。実施の形態に対して当業者が思いつく変形を施して得られる形態、および、複数の実施の形態における構成要素を任意に組み合わせて実現される別の形態も本発明に含まれる。
【０１２９】
例えば、上述した実施の形態では、環構造体の断面形状を円形としているが、環構造体の断面形状はこれに限らず非線形の連続する形状であればよい。例えば、環構造体の断面形状は楕円形であってもよい。
【０１３０】
また、上述した実施の形態では、計測対象物の複数箇所に光ケーブルが配置される構成としたが、光ケーブルを配置する位置は一箇所であってもよいし、複数箇所であってもよい。また、光ケーブルを環構造体の周囲の複数箇所に配置する場合には、これらの光ケーブルは１本の光ケーブルで一筆書きとなるように環構造体の周囲に配置されていてもよい。
【０１３１】
また、上述した実施の形態では、計測対象物の環構造体の外周および内周に光ケーブルが配置される構成としたが、計測対象物のひずみのみを計測する場合には、外周のみに光ケーブルが配置される構成としてもよい。
【０１３２】
また、上述した実施の形態では、スペクトル観測部３０、ひずみ計測部４０および変位算出部５０は、一台のコンピュータで一体に構成されていてもよいし、複数のコンピュータまたは機器でそれぞれ構成されていてもよい。
【０１３３】
また、上述した実施の形態では、評価関数Ｅを最小化するｑを得るための方法として滑降シンプレックス法を用いることとしたが、滑降シンプレックス法に限らず他の方法を用いてもよい。
【産業上の利用可能性】
【０１３４】
本発明にかかる変位計測方法および変位計測装置は、トンネル、パイプライン、ガスタンク、航空機など断面が環状の構造物のひずみおよび変位の計測に有用である。
【符号の説明】
【０１３５】
１０  環構造体
２０ａ、２０ｂ  光ファイバ
３０  スペクトル観測部
４０  ひずみ計測部
５０  変位算出部
６０、７０、７０ａ、７０ｂ  微小片（微小直線梁）
６０ａ、６０ｂ、６０ｃ  微小直線梁
６１  直線梁
１００  変位計測装置
【図１】
2017058292.tif 000026
【図２】
2017058292.tif 000027
【図３】
2017058292.tif 000028
【図４】
2017058292.tif 000029
【図５】
2017058292.tif 000030
【図６】
2017058292.tif 000031
【図７】
2017058292.tif 000032
【図８】
2017058292.tif 000033
【図９】
2017058292.tif 000034
【図１０】
2017058292.tif 000035
【図１１】
2017058292.tif 000036
【図１２】
2017058292.tif 000037
【図１３】
2017058292.tif 000038
【図１４】
2017058292.tif 000039
【図１５】
2017058292.tif 000040
【図１６】
2017058292.tif 000041