(57)【要約】【課題】柱状構造物のたわみ量を精度よく推定する。【解決手段】柱状構造物の範囲内の各部分領域に対応する３次元点群から選択された複数の３次元点から、たわみがないとした各部分領域の形状を規定する複数の部分領域モデルの複数の第１のモデルパラメータを推定し（Ｓ３＿２）、複数の第１のモデルパラメータから、全体領域の３次元点群に対し、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータを推定し、各部分領域モデルが各部分領域に存在する尤度値及び大域的モデルと部分領域モデルとの整合性の評価値から、第２のモデルパラメータの評価値を複数回計算し（Ｓ３＿３〜Ｓ３＿５）、評価値を最適にする第２のモデルパラメータから、柱状構造物のたわみ量を算出する。よって、全体領域の点群形状を考慮しつつ、部分領域に最適な形状モデルを推定でき、柱状構造物のたわみ量を精度よく推定できる。【選択図】図７
2015232513.tif 000002
【特許請求の範囲】
【請求項１】
柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定し、
前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定すること、及び
前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記第１のモデルパラメータと、推定された第２のモデルパラメータとの整合性を評価する前記第２の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算すること
を複数回行う幾何モデル当てはめ部と、
前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、
を含むたわみ推定装置。
【請求項２】
柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定し、
前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定すること、及び
前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータにより表わされる柱状構造物の全体形状が、凸関数で近似される形状であるかどうかの整合性を評価する第３の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算すること
を複数回行う幾何モデル当てはめ部と、
前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、
を含むたわみ推定装置。
【請求項３】
前記幾何モデル当てはめ部は、
前記複数の部分領域の各々に対し、前記部分領域に対応する３次元点群を、前記柱状構造物の高さ方向に垂直な平面に投影して得た複数の点の内、最も多くの３次元点が投影された点に、前記柱状構造物の中心軸から向かう代表分布方向を求め、
前記複数の部分領域の各々について、前記部分領域の中心軸上の点から前記代表分布方向に向かう直線と、前記部分領域に対して推定された前記第１のモデルパラメータにより表される柱状構造物の表面との交点を求め、
前記複数の部分領域の各々について求められた交点に基づいて、前記第２のモデルパラメータを推定する
請求項１又は２記載のたわみ推定装置。
【請求項４】
前記幾何モデル当てはめ部は、
前記柱状構造物の中心軸を凸関数で表現した３次元モデルを規定するパラメータベクトルのうちから選択されたパラメータを複数回変化させる毎に、前記パラメータを変化させた前記パラメータベクトルにより形状が表される柱状構造物の表面に前記範囲内の３次元点群で、かつ局所形状が前記柱状構造物と類似している３次元点群が存在する度合いを示すスコア値を求め、
前記範囲内の３次元点群から、前記求められた最大のスコア値を求めるときに前記柱状構造物の表面に存在すると判断された３次元点群を抽出し、
前記抽出した３次元点群に基づいて、前記複数の部分領域の各々に対して前記第１のモデルパラメータを推定する
請求項１〜３の何れか１項に記載のたわみ推定装置。
【請求項５】
柱状構造物の位置に対応する範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群に基づいて、前記柱状構造物の中心軸を凸関数で表現した３次元モデルを規定するパラメータベクトルのうちから選択されたパラメータを複数回変化させる毎に、前記パラメータを変化させた前記パラメータベクトルにより形状が表される柱状構造物の表面に前記範囲内の３次元点群で、かつ局所形状が前記柱状構造物と類似している３次元点群が存在する度合いを示すスコア値を求め、
前記求められた最大のスコア値に対応する前記パラメータベクトルを出力する
インライア推定部と、
前記出力されたパラメータベクトルに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、
を含むたわみ推定装置。
【請求項６】
請求項１〜請求項５の何れか１項記載のたわみ推定装置の各部としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、たわみ推定装置、及びプログラムに係り、特に、柱状構造物のたわみ量を推定するたわみ推定装置、及びプログラムに関する。
【背景技術】
【０００２】
近年、トンネルや橋などのインフラ構造物の劣化・老朽化が問題となり、自治体において維持管理やリスク評価が課題となっている。それと同時に、インフラ設備の数の多さが問題となり、構造物の安全性の評価の効率化やコスト削減を目的とした技術開発が求められている。特に、電柱や信号などのコンクリートポールのような柱状構造物の点検業務は、その数の多さから全数点検作業の稼働コストが大きい。
【０００３】
構造物のひずみ等の形状変形に関する情報（以下、たわみと呼ぶ）は倒壊につながる重要な安全指標の一つであり、従来はひずみゲージとよばれるセンサを構造物表面にとりつけることで計測を行ってきた。近年では、コストの問題やセンサ自体の耐久性を考慮して、カメラやレーザー等の光学的なセンサの開発が進み非接触な計測が可能となった。これにより、ひずみゲージ等のセンサ設置のための仮設足場などが不要となるため、短時間で複数の位置の情報を計測することが可能となった（非特許文献１）。また、柱状の構造物については、非特許文献２のような専用の製品活用により、より効率的な点検が可能となった。これにより、ある基準以上のたわみ量を有するものを不安全と判定し、その柱状構造物について保守点検を優先して行うことができる。
【０００４】
ただし、上記の非特許文献１においては、測定場所に基準点を設定して計測するときに、人手の作業が依然として必要である。非特許文献２のような専用装置を用いても、計測箇所の指定や柱状物体の中心軸位置を決定するためには人手による指示作業が必要であり、膨大な数の柱状構造物の点検作業における稼働コストは大きい。
【０００５】
一方、モバイルマッピングシステム（ＭＭＳ）と呼ばれるカメラやレーザースキャナを搭載した車が、街中を走行することで道路周辺の構造物表面の形状を計測できるシステムが開発され利用されつつある。このシステムは、ＧＰＳ（全地球測位システム）やIMS(慣性計測装置)を用いて物体の表面を３次元の座標情報として記録できる。
【０００６】
また、近年では、非特許文献３のような技術により柱状構造物を検出する技術があり、この技術とＭＭＳを用いることで、屋外環境下での柱状構造物の位置を自動検出することは可能である。しかしながら、ＭＭＳで計測した点群からたわみ量は推定できないため、その柱状構造物が安全か不安全であるかの判定ができない。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】「光学的全視野計測法を用いた鋼・コンクリート部材の変形・ひずみ計測に関する基礎的研究」 : Nagasaki University (長崎大学), 博士(工学) (2012-09-20). 記述: 長崎大学学位論文
【非特許文献２】アイサンテクノロジ柱状物体計測ツール「バーム」、[online]、[平成２６年５月１９日検索]、インターネット＜ＵＲＬ：http://www.aisantec.co.jp/products-services/atstation/baumstation.html＞
【非特許文献３】新垣、島村、新井、谷口、「３次元点群からの局所形状と大域的な形状モデルを用いた柱状物体検出」、信学技報, vol. 112, no. 441, PRMU2012-131, pp. 7-12, 2013年2月.
【非特許文献４】林和則、「確率伝搬法とその応用」、数理解析研究所講究録 1616, 16-40, 2008-10
【非特許文献５】Yongwei Miao, Jieqing Feng and Qunsheng Peng, 「Curvature Estimation of Point-Sampled Surfaces and Its Applications」、Computational Science and Its Applications−ICCSA 2005, Lecture Notes in Computer Science, pp. 1023-1032.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
従来、柱状構造物のたわみ計測をするための稼働時間が、設備の数の膨大さから問題視されてきた。現地での計測作業の時間を削減するために、点群データをもとに電柱や信号等の柱状構造物のたわみ形状について、3次元ＣＡＤ装置などの専用編集ツールの利用が考えられる。しかし、こういったＣＡＤツールでは、目視による中心軸の位置の指定作業が必須となり、手動による位置指定の精度が低い場合にはたわみ量の推定精度が低くなる。そのため、ＣＡＤ等を用いた計測には熟練の操作者が必要である。また、画面上での操作は現場への移動時間を減らせるが、画面上での細かな作業のため作業者に大きな負担をしいる。
【０００９】
例えば、図１に示すように、柱状構造物について、様々な高さの断面を作成し、その断面における中心位置を手動により指定することで、柱状構造物の中心軸の形状を推定することも可能であるが、非計測領域が狭い場合には中心位置を精度よく求めることは難しい。
【００１０】
特に、ＭＭＳのような移動する車に搭載されたセンサーシステムで屋外の構造物を計測した場合、原理的に車が走行する道路側の被写体の表面形状は測定できるが、道路反対側の表面形状は計測できないため、計測領域が狭い（表面積が小さい）ことが多い。このため、図１の断面Ａ−Ａのように、計測範囲が限られるため点群の数は少なくなるため、計測ノイズや隣接する物体（図２（Ａ）も参照）の影響により、中心位置の推定精度は低くなる。
【００１１】
そのため、本来はたわみが生じていない真っ直ぐな柱状物体であっても、推定した中心位置（図２（Ｂ））を、図２（Ｃ）に示すように滑らかに連結した曲線モデルでは、例えばスプライン曲線で近似した場合、凹凸が生じてしまう。つまり、図２（Ｃ）のような、連結した形状からでは、柱状構造物の形状表現精度が低いために、精度よくたわみ判定をすることが困難である。以下では、たわみ判定とは、設備管理者が目視点検によって不安全と判定するほど大きな所定量の形状変形（屈曲、歪等）を有しているかどうか判定することとする。
【００１２】
本発明は、このような事情を考慮してなされたものであり、その目的は、被計測表面積が狭い場合でも、計測ノイズに頑健に柱状構造物のたわみ量を精度よく推定することができるたわみ推定装置、及びプログラムを提供することにある。
【課題を解決するための手段】
【００１３】
上記の目的を達成するために請求項１に記載のたわみ推定装置は、柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定し、前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定すること、及び前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記第１のモデルパラメータと、推定された第２のモデルパラメータとの整合性を評価する前記第２の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算することを複数回行う幾何モデル当てはめ部と、前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、を含む。
【００１４】
幾何モデル当てはめ部は、柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定する。
また、幾何モデル当てはめ部は、以下の第１及び第２を複数回行う。
第１に、幾何モデル当てはめ部は、前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定する。
第２に、幾何モデル当てはめ部は、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記第１のモデルパラメータと、推定された第２のモデルパラメータとの整合性を評価する前記第２の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算する。
【００１５】
たわみ量算出部は、前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出する。
【００１６】
このように本発明は、複数回上記処理を繰り返すことにより、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータの複数の評価値を計算し、複数回繰り返すことにより得られた複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、柱状構造物のたわみ量を算出するので、被計測表面積が狭い場合でも、計測ノイズに頑健に柱状構造物のたわみ量を精度よく推定することができる。
【００１７】
請求項２のたわみ推定装置は、柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定し、前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定すること、及び前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータにより表わされる柱状構造物の全体形状が、凸関数で近似される形状であるかどうかの整合性を評価する第３の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算することを複数回行う幾何モデル当てはめ部と、前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、を含む。
【００１８】
幾何モデル当てはめ部は、柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定する。
また、幾何モデル当てはめ部は、以下の第１及び第２を複数回行う。
第１に、幾何モデル当てはめ部は、前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定する。
第２に、幾何モデル当てはめ部は、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータにより表わされる柱状構造物の全体形状が、凸関数で近似される形状であるかどうかの整合性を評価する第３の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算する。
【００１９】
たわみ量算出部は、前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出する。
【００２０】
このように本発明は、複数回上記処理を繰り返すことにより、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータの複数の評価値を計算し、複数回繰り返すことにより得られた複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、柱状構造物のたわみ量を算出するので、被計測表面積が狭い場合でも、計測ノイズに頑健に柱状構造物のたわみ量を精度よく推定することができる。
【００２１】
請求項５に記載のたわみ推定装置は、柱状構造物の位置に対応する範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群に基づいて、前記柱状構造物の中心軸を凸関数で表現した３次元モデルを規定するパラメータベクトルのうちから選択されたパラメータを複数回変化させる毎に、前記パラメータを変化させた前記パラメータベクトルにより形状が表される柱状構造物の表面に前記範囲内の３次元点群で、かつ局所形状が前記柱状構造物と類似している３次元点群が存在する度合いを示すスコア値を求め、前記範囲内の３次元点群から、前記求められた最大のスコア値を求めるときに前記柱状構造物の表面に存在すると判断された３次元点群を抽出し、前記抽出した３次元点群に基づいて、前記柱状構造物の形状を規定する柱状構造物モデルパラメータを計算するインライア推定部と、前記計算された柱状構造物モデルパラメータに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、を含む。
【００２２】
インライア推定部は、柱状構造物の位置に対応する範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群に基づいて、前記柱状構造物の中心軸を凸関数で表現した３次元モデルを規定するパラメータベクトルのうちから選択されたパラメータを複数回変化させる毎に、前記パラメータを変化させた前記パラメータベクトルにより形状が表される柱状構造物の表面に前記範囲内の３次元点群で、かつ局所形状が前記柱状構造物と類似している３次元点群が存在する度合いを示すスコア値を求め、前記範囲内の３次元点群から、前記求められた最大のスコア値を求めるときに前記柱状構造物の表面に存在すると判断された３次元点群を抽出し、前記抽出した３次元点群に基づいて、前記柱状構造物の形状を規定する柱状構造物モデルパラメータを計算する。
【００２３】
たわみ量算出部は、前記計算された柱状構造物モデルパラメータに基づいて、前記柱状構造物のたわみ量を算出する。
【００２４】
このように本願発明は、パラメータを変化させた前記パラメータベクトルにより形状が表される柱状構造物の表面に前記範囲内の３次元点群で、かつ局所形状が前記柱状構造物と類似している３次元点群が存在する度合いを示すスコア値を求め、３次元点群から、最大のスコア値を求めるときに柱状構造物の表面に存在すると判断された３次元点群を抽出し、抽出した３次元点群に基づいて、柱状構造物の形状を規定する柱状構造物モデルパラメータを計算し、計算された柱状構造物モデルパラメータに基づいて、柱状構造物のたわみ量を算出するので、被計測表面積が狭い場合でも、計測ノイズに頑健に柱状構造物のたわみ量を精度よく推定することができる。
【００２５】
本発明に係るプログラムは、請求項１〜請求項５の何れか１項記載のたわみ推定装置としてコンピュータを機能させるためのプログラムである。
【発明の効果】
【００２６】
以上説明したように、本発明のたわみ推定装置、及びプログラムによれば、被計測表面積が狭い場合でも、計測ノイズに頑健に柱状構造物のたわみ量を精度よく推定することができる、という効果が得られる。
【図面の簡単な説明】
【００２７】
【図１】（Ａ）は、柱状構造物における３次元点群を、（Ｂ）は、高さＺ=Ｚｈでの円筒当てはめ結果の概念を示す図である。
【図２】（Ａ）は、柱状構造物に付属品が取り付けられている様子を、（Ｂ）は、各高さで独立して推定した中心軸の位置を示し、（Ｃ）は、（Ｂ）における中心軸を滑らかに連結した様子を示す図である。
【図３】（Ａ）及び（Ｂ）は、電柱のたわみの原因の概念を示し、（Ｃ）は、柱状構造物の中心軸を凸関数で表現した曲線モデルのパラメータ示す図である。
【図４】柱状構造物のたわみ量を示す図である。
【図５】第１の実施形態のたわみ推定装置１００のブロック図である。
【図６】たわみ推定処理を示すフローチャートである。
【図７】図６のステップＳ３における幾何モデル当てはめ処理のフローチャートである。
【図８】重複した部分領域を設定する概念を示す図である。
【図９】図７のステップＳ３＿２の部分領域モデルパラメータ推定処理(円筒モデル当てはめ）のフローチャートである。
【図１０】（Ａ）〜（Ｃ）は、各部分領域モデルの中心位置のプロットと凸関数を示す図である。
【図１１】図７のステップＳ３＿３の大域的モデルパラメータ推定処理のフローチャートである。
【図１２】３次元点ｐｍと凸関数の中心軸との距離の概念を示す図である。
【図１３】（Ａ）〜（Ｃ）は、代表分布方向及び凸部交点を求める方法を示す図である。
【図１４】第２の実施形態における図６のステップＳ３における幾何モデル当てはめ処理のフローチャートである。
【図１５】（Ａ）及び（Ｂ）は、高さ方向の位置の違いによるアウトライアの割合の違いの概念を示す図である。
【図１６】第４の実施形態におけるたわみ推定装置１００のブロック図である。
【図１７】第４の実施形態の変形例におけるたわみ推定装置１００のブロック図である。
【図１８】第４の実施形態における図６のステップＳ３における幾何モデル当てはめ処理のフローチャートである。
【図１９】たわみ円筒モデルパラメータの概念を示す図である。
【図２０】第４の実施形態における図１８のステップＳ３＿１＿２におけるインライア推定処理のフローチャートである。
【発明を実施するための形態】
【００２８】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
［発明の概要］
（第１の概要）
本発明では柱状構造物のたわみに関する事前知識を導入することにより、計測ノイズに頑健に形状情報を取り出す。以下の実施形態では柱状構造物を電柱として説明するが、それ以外の柱状構造物でもよい。
【００２９】
電柱のような柱状構造物がたわむ原因は、図３（Ａ）に示すように、電力線や通信線などのケーブルの張力の影響や変圧器などの付属品の荷重である。図３（Ｂ）に示すように、ケーブル等の張力および付属品の荷重等の合力Ｆにより生じた柱状構造物への曲げモーメントによりたわみが生じる。
【００３０】
このとき、柱状構造物の中心軸は直線と仮定するよりも曲線モデルで表現した方が精度よく表現することが可能である。本発明では、図３（Ｃ）に示すように個々の柱状構造物ごとに設定したローカル座標系（Ｕ−Ｖ座標）において、柱状構造物の中心軸を凸関数で表現した曲線モデルのパラメータを推定することで、たわみ量を求める。これにより、複数の凸形状が生じるような物理的に生じにくい形状、例えば複数の凹凸形状を有する形状を誤って推定することを避けることができる。
【００３１】
本発明では、部分的に区切った領域ごとに推定した中心軸上の位置ベクトル群に対して、曲線モデルを当てはめることでたわみ量を推定する。ここで、当てはめるとは、最適な形状モデルパラメータを探索することであり、評価関数の値が高くなるような繰り返し処理を行い、最も評価値が良いときのモデルパラメータを求めることを意味する。
【００３２】
（第２の概要）
しかしながら、レーザー等の計測機器での測定では、柱状構造物の下端部ほど計測密度が高くなり上端部ほど計測密度が小さくなりやすく、また看板やケーブルなどの隣接物がオクルージョンとなり、計測表面積も場所により異なる。単純にパラメータの数を増やすと、被計測表面積が大きい部分に過度に影響を受けて曲線モデルが当てはまる可能性がある。つまり、被計測表面積が不均一性により、曲線モデルパラメータの推定精度が低下するという新たな課題がある。
【００３３】
そこで、本発明では上記の被計測表面積が不均一性の課題について、柱状構造物の点群を、部分的に区切った領域（部分領域）ごとに設定した正規化したモデル当てはめスコアを用いることで精度が低下することを抑制する。具体的には、部分領域ごとに最適なモデルパラメータ候補を複数検出し、それらのモデルパラメータ（例えば、位置および中心軸）が凸関数の中心軸と整合性をとれるような最適な組み合わせを求めることで、中心軸の形状を推定する。
【００３４】
（第３の概要（第３の実施形態））
ＭＭＳで計測した場合，図１に示すようにある特定の方向から計測した表面の点群が多く取得されやすい。それが原因となり、部分領域モデルの中心位置の精度が低くなりやすい。特に、柱状構造物の長さが短い場合には、大域的なモデルの大きさも短くなるため、曲線モデルのパラメータ推定精度が低下する可能性がある。
【００３５】
そこで、本発明では、事前に計測方向を推定し、その方向上の最も突出している部分の位置の軌跡から中心軸を推定する。つまり、部分領域モデルの推定した中心位置ベクトル群の代わりに凸部交点の位置ベクトル群を用いて曲線モデルパラメータを推定する（図１３参照）。
【００３６】
ただし、凸部交点の位置は、部分領域モデル表面とＸＹ平面に平行な平面との断面形状（楕円）と計測面（代表分布方向と垂直な面）との最短距離として求められる。
【００３７】
（第４の概要（第４の実施形態））
第１の実施形態〜第３の実施形態により、被計測領域の表面積の不均一性にも頑健にたわみ量を推定することが可能である。しかしながら、図１５に示すように、柱状構造物の柱に沿うようにつけられた電力線や配管や看板等の付属品が多い位置においては、繰り返し処理の初期の段階での部分領域モデルのパラメータ推定精度が低下しやすい。初期の段階で推定精度が低い部分領域モデルの数が多いと、大域的なモデルの推定精度も低下する可能性がある。そのため、推定精度を低下させるアウトライアを、たわみ量を推定する事前に除去することが望ましい。
【００３８】
そこで、本発明において、インライア推定部１１９により柱状構造物の表面以外の点群を事前に除去することで、部分領域モデルの推定精度の低下を抑制する。事前処理として、入力点群全体を用いてたわみ円筒モデル当てはめを行い、求めた円筒モデルの表面上に存在する点群をインライアと判定する処理を追加する。
【００３９】
［実施形態の概説］
本発明は、電柱や信号などの柱状構造物をレーザースキャナで計測し、得られた３次元点群を用いてその柱状構造物のたわみ量を推定する技術である。
【００４０】
たわみ推定の前処理として、非特許文献３のような既存技術を用いて、柱状構造物を円筒モデルと見立てて検出し、位置と柱状構造物の中心軸の傾きを求める。次に、柱状構造物の中心軸周辺の点群を用いて、その構造物のたわみ量を推定する。
【００４１】
ここで、本実施形態における柱状構造物のたわみとは、図４に示すように、柱状物体の上端部の中心位置から下端部での中心位置ｑbottomにおける接線（柱状物体の傾きの方向）への垂線の足までの長さｄhのこととする。
【００４２】
以下の実施形態では、具体的な例としてレーザーレンジファインダにより取得した点群を用いた電柱のたわみ推定方法について説明するが、信号柱、標識柱などのそれ以外の柱状構造物のたわみ推定も可能な技術である。
【００４３】
実施形態での３次元情報とは、緯度、経度、海抜（高さ）情報でもよいし、ユーザーが設定した特定の位置を原点とした３次元ユークリッド座標系でも極座標系でもよい。以下の例では、ユーザーが設定した原点における３次元ユークリッド座標系（各方向をＸ，Ｙ，Ｚ座標とする）を想定する。各座標の単位はメートル(m)やセンチメートル(ｃｍ)、ミリメートル(ｍｍ)で表現するが、他の単位でもよい。３次元点とは、各点に上記の３次元座標に、その点群が撮影された時刻や、レーザーの反射強度や赤・青・緑などの色情報等が付与されている点である。３次元点に付与される情報に制限はないが、少なくとも位置情報（Ｘ，Ｙ，Ｚ座標）は付与されたものであり、３次元点群とはその３次元点が2点以上集まった集合である。また、本実施形態において、ユークリッド座標系のZ軸の方向は鉛直（高さ）方向を意味し、Ｘ軸とＹ軸の張る二次元平面は水平面を意味するとする。
【００４４】
[第１の実施の形態]
（たわみ推定装置全体の説明）
図５は、本発明の第１の実施形態による柱状構造物のたわみ推定装置１００の構成を示すブロック図である。たわみ推定装置１００は、ＣＰＵ（Central Processing Unit）、ＲＡＭ、２次記憶装置、プログラムを記憶するＲＯＭを備えたコンピュータで構成される。たわみ推定装置１００には、被写体計測部１０１、入力部１０２、及び出力部１０７が接続されている。
【００４５】
図５において、被写体計測部１０１は、レーザーレンジファインダや、赤外線センサ、または超音波センサなど、被写体とセンサとの距離を測定可能な装置であり、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群を出力する。例えば、レーザーレンジファインダをＧＰＳ（Global Positioning System）が搭載された車の上、もしくはＧＰＳの搭載された飛行機に搭載し、移動しながら計測することで、屋外の環境の人工物、例えばワイヤ、建物、木、道路や道路以外の地面など不特定多数の被写体の３次元位置を計測するシステムである。本実施形態では、被写体計測部１０１として、車上にＧＰＳとレーザーレンジファインダとが搭載されているＭＭＳ（Mobile Mapping System）を想定している。ただし、被写体計測部１０１は、ある特定の位置（交差点など）１箇所からの計測部であってもよい。
【００４６】
入力部１０２は、マウスやキーボードなどのユーザーインターフェースであり、たわみ推定処理装置で使用するパラメータを入力するものである。また、パラメータを記憶したＵＳＢメモリなどの外部記憶媒体でもよく、記憶部にパラメータを供給する。例えば、地図上の柱状構造物の位置情報などが登録されている。
【００４７】
たわみ推定装置１００は、上記２次記憶装置に設けられた記憶部１０３、柱状構造物検出部１０４、たわみ推定処理部１０５、及び構造物安全性判定部１０６を備えている。
【００４８】
記憶部１０３は、３次元点群記憶部１１０と演算処理用パラメータ記憶部１１１と柱状構造物パラメータ記憶部１１２とから構成される。３次元点群記憶部１１０は、被写体計測部１０１から取得した３次元点群を記憶し、柱状構造物検出部１０４とたわみ推定処理部１０５に供給する。演算処理用パラメータ記憶部１１１は、入力部１０２から取得したパラメータの値を記憶し、柱状構造物検出部１０４とたわみ推定処理部１０５に供給する。柱状構造物パラメータ記憶部１１２は、たわみ推定処理部１０５から取得したたわみ量を記憶し、構造物安全性判定部１０６へ供給する。これら３つの手段により構成される記憶部１０３は、例えば、ＨＤＤ（Hard Disk Drive）や、ＳＳＤ（Solid State Drive）等のハードウェアによる記憶装置である。
【００４９】
柱状構造物検出部１０４は、３次元点群記憶部１１０から点群を取得して、柱状構造物を検出し、柱状構造物の位置と傾きと長さ等のパラメータを推定し、柱状構造物パラメータ記憶部１１２へ供給する。本実施形態では、非特許文献３のような既存技術を用いることとする。もしくは、既に設備管理ＤＢなどに登録されている場合には、入力部１０２から情報を得られるため、柱状構造物検出部１０４の処理を省くことも可能である。
【００５０】
たわみ推定処理部１０５は、入力処理部１１２Ａ、取得部１１３、幾何モデル当てはめ部１１４、たわみ量算出部１１５、出力処理部１１６、及び繰り返し処理部１１７を備えている。たわみ推定処理部１０５は、記憶部１０３の３次元点群記憶部１１０から３次元点群を取得し、柱状構造物パラメータ記憶部１１２から柱状構造物パラメータを取得し、推定したたわみ量を記憶部１０３へ供給する。
【００５１】
入力処理部１１２Ａは、３次元点群、及び後述するパラメータを入力する。取得部１１３は、３次元点群と柱状構造物パラメータとに基づいて、柱状構造物Iの中心軸から半径Ｒcut（=ＲI+ΔR）以内に含まれる３次元点群を取得し、幾何モデル当てはめ部１１４へ供給する。
【００５２】
幾何モデル当てはめ部１１４は、部分領域モデル推定処理を実行する部分領域モデル推定部１２０と、大域的モデル推定処理を実行する大域的モデル推定部１２１とを備えている。大域的モデル推定処理により、柱状構造物のモデルパラメータが得られる。モデルパラメータはたわみ量算出部１１５へ供給される。ここで、モデルパラメータとは、柱状構造物の位置、傾き、テーパ、中心軸の曲がり具合など、その柱状構造物の形状に関する物理的な量のことを意味する。
【００５３】
たわみ量算出部１１５では、幾何モデル当てはめ部１１４で推定した柱状構造物Ｉのモデルパラメータから、たわみ量を算出する処理を行い、算出したたわみ量を記憶部１０３の柱状構造物パラメータ記憶部１１２へ供給する。
【００５４】
構造物安全性判定部１０６は、柱状構造物パラメータ記憶部１１２から柱状構造物のモデルパラメータを取得し、出力部１０７へ供給する。このとき、たわみ量から不安全と判定されたときには、不安全と判定された柱状構造物全ての座標および警告信号を出力部１０７へ供給する。
【００５５】
出力部１０７は、ＣＲＴ（Cathode Ray Tube）、ＬＣＤ（Liquid Crystal Display）、ＰＤＰ（Plasma Display Panel）等のディスプレイやプリンタなどであり、本装置により求めたたわみ量および警告表示の有無を画面表示もしくは印刷等により保守管理担当者へ示すことに用いる。ここで、警告信号を構造物安全性判定部１０６から取得した場合には、不安全と判定された柱状構造物の座標情報、たわみ量をリスト化したものを表示してもよい。
【００５６】
次に、本実施形態の作用を説明する。
（たわみ推定処理部１０５の動作説明）
まず、本実施形態のたわみ推定装置１００によるたわみ推定処理部１０５全体の処理を、図６を参照して説明する。
【００５７】
ステップＳ１において、入力処理部１１２Ａは、３次元点群記憶部１１０から３次元点群を、演算処理用パラメータ記憶部１１１からたわみ推定処理部の演算処理で使用する演算処理用パラメータを、柱状構造物パラメータ記憶部１１２から柱状構造物パラメータをそれぞれ入力する。
【００５８】
ここで、各柱状構造物を区別する番号をＩ（∈１,２,・・・ＮI）、柱状構造物の総数をＮＩと表示する。柱状構造物パラメータには、柱状構造物の下端部の中心位置の３次元座標pI_bottom、中心軸の方向
【００５９】
【数１】
2015232513.tif 000003
【００６０】
、構造物の長さＬI、半径ＲI等の柱状構造物に関する形状情報がある。
【００６１】
ステップＳ８において、繰り返し処理部１１７は、Ｉに１を代入する。
【００６２】
ステップＳ２において、取得部１１３は、柱状構造物Ｉの中心軸から半径Ｒcut（＝ＲI＋ΔＲ）以内に含まれる３次元点群を取得する。この半径ΔＲは演算処理用パラメータ記憶部１１１に登録されているパラメータである。以降では、このようなパラメータを「実験的に決めるパラメータである」と記載する。本実施形態において、ΔＲ＝０．０５(ｍ)とした。
【００６３】
ここで、取得した３次元点群の３次元点の数をＮm、各３次元点を区別する番号をｍ(∈１,２,..,Ｎm)で表示する。注目点ｍの位置座標をベクトルで表すと、ベクトルｐmは次式によりＸ軸、Ｙ軸、Ｚ軸の成分で構成される。ただし、右上の添え字の記号「Ｔ」は転置を意味とする。
【００６４】
【数２】
2015232513.tif 000004
【００６５】
ただし、Ｘm、Ｙm、ＺmはＸ，Ｙ，Ｚ軸の成分である。
【００６６】
柱状構造物Ｉの３次元点群は、下記の式を満たす３次元点として求められる。
【００６７】
【数３】
2015232513.tif 000005
【００６８】
ステップＳ３において、幾何モデル当てはめ部１１４は、ステップ２で取得した柱状構造物Ｉの３次元点群について、柱状構造物Ｉをモデル化したたわみ付き円筒モデルの当てはめ処理を行い、たわみ付き円筒モデルのモデルパラメータを推定する。ここで、たわみ付き円筒のモデルパラメータには、柱状構造物の下端部の中心位置ｑ_bottom、柱状構造物の上端部の中心位置ｑ_top、凸関数で表現した柱状構造物の中心軸、柱状構造物の下端部における接線方向
【００６９】
【数４】
2015232513.tif 000006
【００７０】
、柱状構造物の下端部での半径の大きさＲ、柱状構造物の長さＬがある。中心軸を表現する凸関数および詳細な処理内容は後述する。
【００７１】
ステップＳ４において、たわみ量算出部１１５は、たわみ量を算出する。ステップ３で求められた柱状構造物の上端部の中心位置ｑtop、柱状構造物の下端部の中心位置ｑbottomにおける傾き方向
【００７２】
【数５】
2015232513.tif 000007
【００７３】
について、たわみ量ｄhは、上端部の中心位置ｑtopから傾き方向
【００７４】
【数６】
2015232513.tif 000008
【００７５】
までの最短距離としては次式により求められる。
【００７６】
【数７】
2015232513.tif 000009
【００７７】
もしくは、単位長さあたりの変形量として、次式により求めてもよい。柱状構造物Ｉの長さをＬとすると次式で求められる。
【００７８】
【数８】
2015232513.tif 000010
【００７９】
ステップＳ５で、繰り返し処理部１１７は、Ｉが柱状構造物の総数ＮＩ以上か否かを判断することにより、全ての柱状構造物について処理が終了したか否かを判断する。全ての柱状構造物について処理が終了していない場合には、ステップＳ７で、繰り返し処理部１１７は、Ｉ←Ｉ＋１とする。処理はステップＳ２へ戻る。なお、記号「←」は代入を意味するとし、数式「Ｉ←Ｉ＋１」はＩの数を１つ増やしてＩに代入するという意味である。全ての柱状構造物について処理が終了した場合には、ステップＳ６で、出力処理部１１６は、たわみ推定処理部１０５で推定したたわみ量を柱状構造物パラメータ記憶部１０３（柱状構造物パラメータ記憶部１１２）へ出力する。これにより、柱状構造物パラメータ記憶部１１２に記憶される。
【００８０】
以下、ステップＳ３について詳細な説明をする。本ステップＳ３では、幾何モデル当てはめ部１１４は、たわみ付き円筒モデルのモデルパラメータの推定処理を行う。以降、たわみ付き円筒モデルとは、円筒モデルの中心軸を凸関数で表現した３次元幾何モデルとする。
【００８１】
本実施形態においては、詳細には後述するが、柱状構造物をモデル化したたわみ付き円筒モデルのモデルパラメータを、次式の部分領域モデル評価項と、大域的モデルＧと部分領域モデルの整合性評価項と、から構成される評価関数Ｅ（Ｘs,Ｇ）の最適化問題として、推定する。
【００８２】
【数９】
2015232513.tif 000011
【００８３】
ただし、Ｍs（ｓ∈１,２,３,...Ｎs）は、たわみがないと仮定した場合の後述する部分領域ｓの形状を規定する部分領域モデルを示し、部分領域モデルＭsは、第１のモデルパラメータで表わされる。第１のモデルパラメータには、部分領域モデルＭsの中心位置、部分領域モデルＭsの半径、部分領域モデルＭsの中心軸方向がある。
【００８４】
Ｇは、柱状構造物の中心軸を凸関数で表現した３次元モデルであり、柱状構造物周辺の全点群に対する大域的なモデルである。大域的モデルＧは具体的には、各部分領域モデルＭsを凸関数で近似することで生成される凸形状の曲線モデルであり、第２のモデルパラメータで表わされる。
【００８５】
関数Ｄ（Ｍs）は、部分領域モデルＭsが部分領域sに存在する尤度であり、部分領域sの第１のモデルパラメータにより形状が表される柱状構造物の表面に、部分領域sに含まれる３次元点群が存在する度合いに基づいて評価する第１の値である。
【００８６】
関数Ｈ（Ｍs,Ｇ）は大域的モデルＧの第２のモデルパラメータと部分領域モデルＭsの第１のモデルパラメータとの整合性を評価する第２の値である。ｋは２つの項目の重みを変える調整用係数である。ｋの値が大きいほど、滑らかに連結されるような部分領域モデルが選択されやすい。逆に、ｋ値が小さいほど、各部分領域モデルは各部分領域の形状に最適なパラメータとなりやすい。
【００８７】
ここで、整合性とは、部分領域モデルＭSの第１のモデルパラメータと、大域的モデルＧの第２のモデルパラメータとの類似度を意味する。本実施形態において、実験的なパラメータkは、部分領域モデル評価項と大域的モデル整合性評価項の比率を同じにするためｋ＝Ｎs(部分領域の総数)とした。
【００８８】
各部分領域については、たわみの無い円筒形状のたわみの無し円筒モデルのパラメータとして十分近似できると考え、各部分領域について、部分領域モデルとして、たわみ無し円筒モデルを採用して、第１のモデルパラメータを推定する。図８に示すように、部分領域ごとに、候補となるたわみ無し円筒モデル群を推定する。各部分領域で推定した全てのたわみ無し円筒モデルの中心位置を凸関数で近似したときに上記評価関数を最大にするような組み合わせを求めれば、上記の評価関数を最適にする結果として、たわみ付き円筒モデルのモデルパラメータが求められる。
【００８９】
しかしながら、全ての組み合わせを処理することは現実的な時間では処理できない。例えば、図８の部分領域の数がＮｓ＝１０個、各部分領域での候補モデル数を１００個とすると、部分領域モデルの全ての組み合わせは１００の１０乗となるため、計算コストが膨大となってしまうために求められない。
【００９０】
そこで、本発明では、部分領域モデルと大域的モデルのそれぞれを逐次的に推定する処理を繰り返し行うことで、評価関数Ｅの最適解を求める（図７）。具体的には、部分領域モデル評価項を最適化するときは、大域的モデルＧを固定として、その固定した大域モデルに対して最も良い部分領域モデルＭsをそれぞれの領域で求める。次に、大域的モデルＧを算出する際には、上記で求めた各部分領域モデルＭsを固定としたときの最大尤度となる凸関数モデルを求める。この処理を収束するまで繰り返す。
【００９１】
図７のステップＳ３＿１において、幾何モデル当てはめ部１１４は、柱状構造物Ｉのパラメータとその柱状構造物Ｉから半径Ｒcut以内の周辺点群Ｐｍを入力する。
【００９２】
ステップ３＿１においては更に、幾何モデル当てはめ部１１４は、複数の部分領域を次のように設定する。
【００９３】
まず、図８に示すような重複する部分領域を設定する。３次元点群Ｐｍに含まれる３次元点をｐｍ（ｍ∈１,２,３,...Ｎｍ）、３次元点群内の３次元点の総数をＮｍ、柱状構造物Ｉの中心軸方向を
【００９４】
【数１０】
2015232513.tif 000012
【００９５】
軸、中心軸に垂直な２つの軸を
【００９６】
【数１１】
2015232513.tif 000013
【００９７】
、下端部の中心位置ｐI_bottomを原点としたローカル座標系を設定すると、図８に示すように３次元点群を中心軸方向
【００９８】
【数１２】
2015232513.tif 000014
【００９９】
に複数の領域へと分割する。ローカル座標原点ｐI_bottomから、中心軸方向へ長さΔｓづつ区切った区間を設定すると、部分領域ＲＳs(ｓ∈１,２,３,...Ｎs)に含まれる点は以下の式を満たす。
【０１００】
【数１３】
2015232513.tif 000015
【０１０１】
ただし、Δｓは実験的に決めるパラメータであり、本実施形態ではΔｓ＝０．５(ｍ)とする。
【０１０２】
ここで数学的な記号の意味を記載すると、記号「→」は方向、つまり長さが１のベクトルであること（２ノルムが１であること）を明示するための記号とし、記号
【０１０３】
【数１４】
2015232513.tif 000016
【０１０４】
はベクトルの２ノルムの大きさとする。記号ベクトル間の記号「・」はベクトルの内積演算、記号「×」はベクトルの外積、スカラー量とベクトル間の記号「・」は掛け算を意味する。
【０１０５】
上式の物理的な意味は、柱状構造物Ｉの周辺点について、中心軸方向
【０１０６】
【数１５】
2015232513.tif 000017
【０１０７】
に射影した点の位置が（ｓ・Δｓ）から±Δｓの範囲に入る点であることを意味する。
【０１０８】
幾何モデル当てはめ部１１４は、後述する処理（Ｓ３＿２、Ｓ３＿３）を閾値ＴＨcnt回繰り返すが、ステップＳ３＿６で、繰り返し処理の回数を識別するｃｎｔに１を代入する。繰り返し処理回数の閾値ＴＨcntは実験的に決めるパラメータであり、本実施形態ではＴＨcnt＝１００００とした。
【０１０９】
ステップＳ３＿２において、幾何モデル当てはめ部１１４は、部分領域モデル推定処理を実行し、評価関数Ｅのモデル評価項である第１の値Ｄ（Ｍs)を算出する。
【０１１０】
具体的には、幾何モデル当てはめ部１１４は、各部分領域ｓにおいてＲＡＮＳＡＣ処理（図９）によりたわみ無し円筒当てはめ処理を行い、最もスコア値の高いたわみ無し円筒モデル（部分領域モデル）の第１のモデルパラメータのスコア値との比率（０．０以上１．０以下）として評価関数Ｅのモデル評価項の第１の値Ｄ（Ｍs)を算出する。上位Ｎｒａｎ個について、部分領域sにおけるたわみ無し円筒モデルの第１のモデルパラメータの候補として登録を行う。
【０１１１】
図９を参照して、部分領域モデルパラメータ推定処理(ＲＡＮＳＡＣによる円筒当てはめ処理）について説明する。ＲＡＮＳＡＣによるたわみ無し円筒当てはめ処理（ステップＳ３＿２）は、各部分領域について実行する。
【０１１２】
ステップＡ＿１において部分領域モデル推定部１２０は、部分領域ｓkの３次元点群（３次元点の総数はＮｓ）を入力する。各３次元点ｐＳはｎ（１,２,・・・Ｎｓ）で識別される。
【０１１３】
部分領域モデル推定部１２０は、後述するステップＡ＿３からステップＡ＿６の処理をＫloop回繰り返すが、ステップＡ＿２において、繰り返し処理の回数を識別するｋに１を代入する。なお、回数Ｋloopは実験的に決めるパラメータであり、本実施形態ではＫloop＝１０００とした。
【０１１４】
ステップＡ＿３において部分領域モデル推定部１２０は、部分領域ｓkにおける３次元点群からランダムに２個の３次元点を選択（サンプリング）し、２個の３次元点間の方向を処理ｋ回目における中心軸の方向
【０１１５】
【数１６】
2015232513.tif 000018
【０１１６】
とする。具体的には２個の３次元点の座標の差分をとり、ノルムを１にするように長さを正規化すればよい。
【０１１７】
【数１７】
2015232513.tif 000019
【０１１８】
の記号「ｋ」はステップＡ＿２からステップＡ＿７の繰り返し処理における処理のｋ(∈1,2,..,Ｋloop)番目の回数を意味する。
【０１１９】
ステップＡ＿４において、部分領域モデル推定部１２０は、部分領域ｓkにおける３次元点群から更にランダムにＮｒ個を選択し、中心軸の方向
【０１２０】
【数１８】
2015232513.tif 000020
【０１２１】
に垂直でかつ部分領域ｓkの重心を通る平面へ射影した座標位置を求める。平面の基底ベクトルは、
【０１２２】
【数１９】
2015232513.tif 000021
【０１２３】
に直交して、かつ重心を通る２つのベクトルであればよい。例えば、グローバル座標のＸ軸とＹ軸について、
【０１２４】
【数２０】
2015232513.tif 000022
【０１２５】
方向に直交するように、グラム・シュミットの正規直交化法を用いて変換してもよい。この平面上での座標位置は、平面の基底の軸方向
【０１２６】
【数２１】
2015232513.tif 000023
【０１２７】
と注目点jの位置ベクトルｐjとの内積で求められる。ランダムに選択した点の位置座標をｐj、ランダムに選択したＮｒ個の点の重心位置をｑgkとすると、平面での座標［ｘ’j、ｙ’j］は次式で求められる。
【０１２８】
【数２２】
2015232513.tif 000024
【０１２９】
ここで、Ｎｒは実験的に決めるパラメータであり、本実施形態ではＮｒ＝５とした。
【０１３０】
ステップＡ_５において、部分領域モデル推定部１２０は、部分領域ｓkの第１のモデルパラメータを求める。具体的には、選択したＮｒ個を上記重心を通る平面上に射影した点群について、最小二乗法により円の当てはめを行い、半径rekと平面上での中心位置の座標
【０１３１】
【数２３】
2015232513.tif 000025
【０１３２】
を部分領域ＲＳＫの第１のモデルパラメータの候補として求める。具体的には、以下の誤差関数Ｅstraightを最小化するように求めればよい。
【０１３３】
【数２４】
2015232513.tif 000026
【０１３４】
ステップＡ＿６において、部分領域モデル推定部１２０は、モデル当てはめスコア値を算出する。これは、ステップＡ＿５で求めたたわみ無し円筒モデルのモデルパラメータの評価を行うことに相当する。評価方法は幾つか存在するが、例えば円筒表面上に存在する点群の数や円筒表面上に存在する点の範囲（面積）等が用いられる。
【０１３５】
本実施形態では、点の数をモデルの評価基準となるスコア値として求めた例を示す。部分領域ｓkに含まれる点ｐsについて、下記のスコア値Score_subregionにより求められる。
【０１３６】
【数２５】
2015232513.tif 000027
【０１３７】
ここで、ＴＨestは実験的に決めるパラメータであり、本実施形態ではＴＨest=０．０２とした。上記式のＦはスカラーの値であることを意味する。
【０１３８】
ステップＡ＿７において、部分領域モデル推定部１２０は、ｋがＫloop以上か否かを判断することにより、繰り返し処理がＫloop回以上終了したか判定を行う。Ｋloop回以上の処理が終了していない場合には、ステップＡ＿９で、部分領域モデル推定部１２０は、ｋに１加算した値をｋに代入する（ｋ←ｋ＋１）。処理は、ステップＡ＿３に戻る。
【０１３９】
繰り返し処理がＫloop回以上終了したと判定された場合には、処理はステップＡ＿８に進む。
【０１４０】
ステップＡ＿８において、部分領域モデル推定部１２０は、最も高いスコア値からＮｒａｎ個のたわみ無し円筒モデルの第１のモデルパラメータの候補と当該Ｎｒａｎ個のスコア値とを出力する。このとき出力するパラメータはたわみ無し円筒モデルの中心位置、半径、中心軸の方向である。中心位置については、平面上の座標系(２次元ベクトル)からグローバル座標系へ変換をした位置座標ｑs(３次元ベクトル)として大域的モデル推定部１２１に出力する。部分領域sにおいて、上位e個目のパラメータについて、次式のように登録する。
【０１４１】
【数２６】
2015232513.tif 000028
【０１４２】
ステップＡ＿８では更に、上位Ｎｒａｎ個の第１のモデルパラメータの候補について、評価関数Ｅのモデル評価項である第１の値Ｄ（Ｍs）を算出して、大域的モデル推定部１２１に出力する。上位一番目のスコア値Score_subregionの値をScore_subregion_maxとすると、上位Ｎｒａｎ個のうちの予め定めた上位e番目の部分領域モデル（たわみ無し円筒モデル）Ｍseのスコア値をScore_subregion(e)とすると、部分領域モデル（たわみ無し円筒モデル）Ｍseにおけるモデル評価項である第１の値Ｄ（Ｍse）は次式で求められる。なお、第１の値は、部分領域モデルＭseが、部分領域ｓkの３次元点群を精度よく表現しているかを評価する値である。
【０１４３】
【数２７】
2015232513.tif 000029
【０１４４】
そして、大域的モデルＧについて、部分領域sにおける中心位置をｑG、その中心位置での接線方向を
【０１４５】
【数２８】
2015232513.tif 000030
【０１４６】
、部分領域モデルＭseの中心位置をｐse、中心軸方向を
【０１４７】
【数２９】
2015232513.tif 000031
【０１４８】
と表記すると、整合性評価項である第２の値は次式で求められる。
【０１４９】
【数３０】
2015232513.tif 000032
【０１５０】
ただし、αは重みを調整する係数であり実験的に決めるパラメータであり、本実施形態においては、例えばα＝０．５とした。ここで、部分領域ｓにおける中心位置をｑGは、大域的モデルの中心軸上の点において、各部分領域ｓの部分領域モデルの中心位置からの最短位置とする。
【０１５１】
このとき繰り返し処理cnt回目の部分領域ｓにおける部分領域モデルＭsの第１のモデルパラメータは次式で求められる。
【０１５２】
【数３１】
2015232513.tif 000033
【０１５３】
また、繰り返し処理が初期の場合、大域的モデルＧを考慮しない次式で求める。
【０１５４】
【数３２】
2015232513.tif 000034
【０１５５】
ステップＡ＿８の処理が終了すると、処理は、図７のステップＳ３＿３に進む。
【０１５６】
ステップＳ３＿３において、大域的モデル推定部１２１は、大域的モデル推定処理を実行する。本ステップＳ３＿３では、大域的モデル推定部１２１は、柱状構造物Ｉの周辺の点群を用いずに、ステップＳ３−２で求めた部分領域モデルの中心位置を通過するＮ次の凸関数をＲＡＮＳＡＣ処理により求める。本実施形態では２次の凸関数（２次曲線）を例に説明する。
【０１５７】
凸関数当てはめの詳細な説明の前に、図１０に中心軸推定処理の概要を説明する。まず、大域的モデル推定部１２１は、図１０（Ａ）に示すように、ステップＳ３−２で推定したたわみ無し円筒モデルのモデルパラメータのうちの中心位置、中心軸の方向（u１、ｕ２）をランダムに２つの部分領域分選び、柱状構造物Ｉの形状を近似するたわみ付き円筒モデルの中心軸(凸関数)が存在する平面を推定する。当該平面は、ｕ軸、ｖ軸により定まる。次に、大域的モデル推定部１２１は、図１０（Ｂ）に示すように、この平面上に、全ての部分領域のたわみ無し円筒モデルのモデルパラメータのうちからランダムに数個選択したモデルパラメータにおける中心位置の点を射影し、平面上（ローカル座標系）での座標を求める。大域的モデル推定部１２１は、この射影した点群を用いて、凸関数のパラメータを算出する。凸関数のパラメータから、たわみ開始位置を算出する。
【０１５８】
以下、ステップＳ３＿３を実現する方法として、ＲＡＮＳＡＣによる２次関数のパラメータ推定処理を説明する（図１１参照）。
【０１５９】
ステップＢ＿１において、大域的モデル推定部１２１は、ステップＳ３＿２で出力した部分領域のたわみ無し円筒モデルのモデルパラメータを入力する。
【０１６０】
大域的モデル推定部１２１は、後述するステップＢ＿３からステップＢ＿６の処理をＴloop回繰り返すが、ステップＢ＿２において、繰り返し処理の回数を識別するｋに１を代入する。なお、Ｔloopは実験的に決めるパラメータであり、本実施形態ではＴloop＝１０００とした。
【０１６１】
ステップＢ_３において、大域的モデル推定部１２１は、２次関数が存在する平面を推定する。まず、大域的モデル推定部１２１は、図１０（Ａ）に示すように、柱状構造物Ｉの下端位置での接線方向ｕ軸を推定する。部分領域で推定したたわみ無し円筒モデルの第１のモデルパラメータを区別する番号をc(c∈1,2,3,..,Nc)、第１のモデルパラメータcの中心軸方向をＵｃ、部分領域において推定した第１のモデルパラメータの数をＮｃとすると、全ての第１のモデルパラメータからランダムに選んだ２つのモデルパラメータ（c1とc2と記載）を用いて、以下の式で柱状構造物Ｉの下端位置における接線方向Ｕkを求める。
【０１６２】
【数３３】
2015232513.tif 000035
【０１６３】
ただし、frandは０〜１の値をとる乱数とする。上記式は、ランダムに選んだ２つの部分領域のたわみ無し円筒モデルのモデルパラメータc1とc2の中心軸方向（uC1、uＣ２）を用いて接線方向Ｕkを計算することを意味する。
【０１６４】
また、後処理のために平面の基底ベクトルであるＵkの第三成分が正となるように向きを調整する。Ｕkの第三成分、つまりグローバル座標系でのＺ軸の値が負の値であるとき、Ｕkの向きを変更する。
【０１６５】
【数３４】
2015232513.tif 000036
【０１６６】
たわみ方向ｖk軸は、Ｕc1からＵkの方向成分を、グラム・シュミットの直交化法により０にして、ノルムが１となるように正規化することで求める。
【０１６７】
【数３５】
2015232513.tif 000037
【０１６８】
平面の中心位置は、選択した２点（ｃ１とｃ２）の円筒パラメータの中心位置（ｑc1,qc2）の重心とする。
【０１６９】
【数３６】
2015232513.tif 000038
【０１７０】
ここで、添え字の「k」は、繰り返し処理がｋ回目のときの推定結果であることを意味する。
【０１７１】
ステップＢ＿４において、大域的モデル推定部１２１は、図１０（Ｂ）に示すように、中心軸の存在する平面をローカル座標系とし、全ての部分領域モデルの第１のモデルパラメータからランダムに選択したＮｊ個の部分領域モデルの第１のモデルパラメータにおける中心位置を、上記中心軸が存在する平面へ射影して、その平面で基底される空間（たわみ円筒ローカル座標系）での座標を求める。
【０１７２】
以下、このたわみ円筒ローカル座標系の軸はｕ軸、ｖ軸で表現し、その座標値（スカラー値）は、それぞれｕとｖで表現する。また位置ベクトルの右上の添え字として記号「’」(ダッシュ)があるものは、ローカル座標に射影した点の位置ベクトルもしくはローカル座標系での座標値であることを明示するためにつけている。
【０１７３】
部分領域で推定したたわみ無し円筒モデルのモデルパラメータについて、ランダムに選択したモデルパラメータを区別する番号をj、その部分領域モデルパラメータの中心位置をｐjとすると、大域的モデルのローカル座標系での座標
【０１７４】
【数３７】
2015232513.tif 000039
【０１７５】
への変換式は次式で求められる。
【０１７６】
【数３８】
2015232513.tif 000040
【０１７７】
ステップＢ＿５において、大域的モデル推定部１２１は、図１０（Ｂ）に示すように、中心軸を表現する凸関数のパラメータを求める。凸関数は２次関数でも３次関数やそれ以上の次元の数でもよいが、本実施形態では２次関数とする。経験的に、安全・不安全の評価を行う柱状構造物は、形状の変形量が微小であるため、ノイズに影響されにくい２次か３次などの低い次元の関数の方が好ましい。本実施形態では２次関数を例として示す。２次数の係数をa,b,cとすると、下記の誤差関数Ｆerrorを最小化するように係数として求める。
【０１７８】
【数３９】
2015232513.tif 000041
【０１７９】
ここで、誤差関数Ｆerrorについて、最小二乗法を用いて係数a,b,cを求めてもよいし、非線形の数値解析手法を用いて求めてもよい。
【０１８０】
本ステップＢ＿５では更に、大域的モデル推定部１２１は、たわみ円筒モデルの下端部（図１０（Ｃ）参照）での中心位置ｑbottomと上端部での中心位置ｑtopを求める。部分領域モデル推定処理Ｓ３−２で求めた中心位置を、上記推定した中心軸へ射影することで求められる。
【０１８１】
高さ方向において最も低い位置にある部分領域ｓ＝０の位置におけるたわみ無し円筒モデル
【０１８２】
【数４０】
2015232513.tif 000042
【０１８３】
の中心位置を
【０１８４】
【数４１】
2015232513.tif 000043
【０１８５】
、高さ方向において最も高い位置にある部分領域ｓ＝Ｎｓの位置でのたわみ無し円筒モデル
【０１８６】
【数４２】
2015232513.tif 000044
【０１８７】
の中心位置を
【０１８８】
【数４３】
2015232513.tif 000045
【０１８９】
とすると、次式で下端部の中心位置ｑbottomと上端部の中心位置ｑtopが求められる。
【０１９０】
【数４４】
2015232513.tif 000046
【０１９１】
以上により求めたパラメータを、繰り返し処理ｋ回目の大域的モデルパラメータＧkとする。
ステップＢ_６において、大域的モデル推定部１２１は、評価関数Ｅの値を算出する。まず、ステップＳ３−２（図７参照）で推定した各部分領域のたわみ無し円筒モデルの最も高いスコア値に対応するモデルパラメータにおける中心位置が、上記凸関数として推定した柱状構造物Ｉの中心軸上に存在するかどうか判定する。
【０１９２】
中心軸上に存在するかどうかは、注目点から中心軸への最短距離により判定する。3次元空間中の曲線モデルと３次元点との最短距離を算出する方法は幾つかあるが、本実施形態では、大域的モデル推定部１２１は、図１２に示すように、ある３次元上の点ｐmからステップＢ_５で設定した凸関数の基底軸（Ｕ軸、Ｖ軸）へ射影した点ｐ’mと凸関数の中心軸との最短位置ｑhとの距離から三平方の定理により算出する。射影点
【０１９３】
【数４５】
2015232513.tif 000047
【０１９４】
は、たわみ円筒ローカル座標系の座標として次式により求められる。
【０１９５】
【数４６】
2015232513.tif 000048
【０１９６】
まず、大域的モデル推定部１２１は、点ｐmをＵＶ平面へ射影した射影点ｐ’mとステップＢ_５で求めた中心軸について、最短距離ｄｉｓｔ1を求める。
【０１９７】
ここで、幾何的な性質から、点からある曲線までの最短距離の位置まで引いた線分と、最短距離の位置での接線は直交するという関係性がある。よって、最短位置
【０１９８】
【数４７】
2015232513.tif 000049
【０１９９】
での接線方向
【０２００】
【数４８】
2015232513.tif 000050
【０２０１】
と射影点ｐ’mから最短位置ｑhまでの方向
【０２０２】
【数４９】
2015232513.tif 000051
【０２０３】
の内積は０となり、次式が成立する。
【０２０４】
【数５０】
2015232513.tif 000052
【０２０５】
未知数（変数）はｕ’hのみであるから、ｕ’hの方程式として解くことができる。この方程式はニュートン法や修正ニュートン法などの数値解析手法で解くことが可能である。
【０２０６】
ただし、この方程式は高次の方程式のため、一般的に複数の解が求められる。そのため、求められる複数の解について、射影点ｐ’mから最も近い位置が最短位置ｑhの位置とする。
【０２０７】
最短距離ｄｉｓｔ1は次式で求められる。
【０２０８】
【数５１】
2015232513.tif 000053
【０２０９】
次に、大域的モデル推定部１２１は、点ｐmと射影した射影点ｐ’mとの距離ｄｉｓｔ２を求める。
【０２１０】
【数５２】
2015232513.tif 000054
【０２１１】
点ｐmと垂線の足の位置ｑhとの距離ｄｉｓｔは、幾何的な関係から次式で求められる。
【０２１２】
【数５３】
2015232513.tif 000055
【０２１３】
ここで、上記では幾何的な関係から最短距離を求めたが、この方法以外にも例えば、ラグランジュの未定乗数法を用いて解いてもよい。制約条件は、最短位置の点が２次関数上に存在することとし、射影点ｐ’mと凸関数曲線上の点において距離が最小となる解（極値）を求めればよい。
【０２１４】
上記最短距離の算出式ｄｉｓｔ（ｐm）用いて、部分領域ｑsと大域的モデルの中心軸との最短位置ｑhとの距離が求められるため、整合性項である第２の値は次式により求められる。
【０２１５】
【数５４】
2015232513.tif 000056
ただし、
【数５５】
2015232513.tif 000057
【０２１６】
は、最短位置ｑhにおける中心軸の接線方向
【０２１７】
【数５６】
2015232513.tif 000058
【０２１８】
であるから、
【０２１９】
【数５７】
2015232513.tif 000059
【０２２０】
よって、整合性項の値は次式で求められる。
【０２２１】
【数５８】
2015232513.tif 000060
【０２２２】
そして、大域的モデル推定部１２１は、評価関数Ｅの評価値を、次の通り求める。なお、右辺の第１項の値は、部分領域モデル推定部１２０から入力されている。
【０２２３】
【数５９】
2015232513.tif 000061
【０２２４】
ステップＢ_７において、大域的モデル推定部１２１は、ｋがＴloop以上か否かを判断することにより、繰り返し処理の回数がＴloop回以上か否かを判断する。繰り返し処理の回数がＴloop回以上の場合には、処理はステップＢ＿８に進む。繰り返し処理の回数がＴloop回以上でない場合には、ステップＢ＿９で、大域的モデル推定部１２１は、ｋ←ｋ＋１とする。処理はステップＢ_３へ戻る。
【０２２５】
ステップＢ_８において、大域的モデル推定部１２１は、評価関数Ｅが最も大きいときのたわみ付き円筒モデルのモデルパラメータを、繰り返し処理ｃｎｔ回目のたわみ付き円筒モデルのモデルパラメータＧcntとして、たわみ量算出部１１５に出力する。たわみ付き円筒モデルのモデルパラメータは、たわみ中心軸の存在するローカル座標系(Ｕ軸、Ｖ軸、Ｗ軸)、凸関数の係数（Ｎ次関数のパラメータ）、たわみ付き円筒モデルの下端部および上端部の中心位置、下端部位置での凸関数（大域的モデルのモデルパラメータ）の接線方向（ローカル座標系のＵ軸）である。
【０２２６】
ステップＢ＿８の処理が終了すると、処理は、図７のステップＳ３＿４に進む。ステップＳ３＿４について、幾何モデル当てはめ部１１４は、ｃｎｔが閾値ＴＨcnt以上か否かを判断することにより、繰り返し処理の回数ｃｎｔが閾値ＴＨcnt以上か否かを判断する。繰り返し処理の回数ｃｎｔが閾値ＴＨcnt以上の場合に、処理はステップＳ３＿７に進む。繰り返し処理の回数ｃｎｔが閾値ＴＨcnt以上でない場合には、ステップＳ３＿５で、幾何モデル当てはめ部１１４は、ｃｎｔ←ｃｎｔ+１とし、処理はステップＳ３＿２に戻る。
【０２２７】
ステップＳ３＿７で、幾何モデル当てはめ部１１４は、評価関数Ｅを最大にするときのたわみ付き円筒モデルのモデルパラメータを出力する。ここで、たわみ付き円筒モデルのモデルパラメータとは、大域的モデルの第２のモデルパラメータであり、たわみ中心軸の存在するローカル座標系(Ｕ軸、Ｖ軸、Ｗ軸)および凸関数の係数および、下端部および上端部の中心位置、下端部位置での凸関数の接線方向（ローカル座標系のＵ軸）である。
【０２２８】
Ｗ軸はＵ軸およびＶ軸に対して垂直な軸とし、グラムシュミットの直交化法として唯一の３次元ベクトルとして求められる。
【０２２９】
（構造物安全性判定部１０６の動作説明）
最後に、構造物安全性判定部１０６について説明する。構造物安全性判定部１０６は、たわみ推定処理部から出力された柱状構造物パラメータを用いて、各柱状構造物について安全か不安全な状態であるかの判定を行う。
【０２３０】
各柱状構造物において、柱状構造物パラメータであるたわみ量ｄhにより以下の方程式を満たすとき不安全な状態であると判定する。
【０２３１】
【数６０】
2015232513.tif 000062
【０２３２】
ここで、閾値ＴＨdhは実験的に決めるパラメータであり、本実施形態では、ＴＨdh＝０．２（ｍ）とした。
【０２３３】
構造物安全性判定部１０６は、不安全と判定された柱状構造物の座標情報、及びたわみ量を含む警告信号を出力部１０７に出力する。
【０２３４】
（出力部１０７の動作説明）
出力部１０７は、警告信号を構造物安全性判定部１０６から取得した場合には、不安全と判定された柱状構造物の座標情報、たわみ量をリスト化したものを、表示したり印刷したりする。
【０２３５】
［第２の実施形態］
次に、本発明の第２の実施形態を説明する。第２の実施形態のたわみ推定装置１００の構成は第１の実施形態と同様であるので、その説明を省略する。
次に、第２の実施形態の作用を説明する。第２の実施形態の作用は、第１の実施形態の作用と同様の部分があるので、主として異なる部分のみ説明する。
【０２３６】
柱状構造物について屈曲（ある高さで急激な曲がりが生じている状態）が存在する場合、２次関数や３次関数などの滑らかな凸関数では中心軸は表現精度が低下する可能性がある。
【０２３７】
そのような状態のときには、凸関数という制限のみを加えたモデル評価関数Ｅを用いることで、中心軸の推定精度を向上することが可能である。本ステップ３において、次式の部分領域モデル評価項と大域的モデルとの整合性項から構成される損失関数Ｅの最小化問題とみなす。
【０２３８】
【数６１】
2015232513.tif 000063
【０２３９】
ただし、Ｘs=(ｓ∈１,２,・・・Ｎｓ）は部分領域ｓで推定した部分領域モデルであり、ＧはＸsにより規定される大域的モデルである。関数Ｄ（Ｘs)は、部分領域モデルＸsにより表される部分領域の形状の表面に各部分領域における点群が存在する度合いに基づいて評価する第１の値であり、関数Ｈ（Ｘs,Ｘv）は各部分領域モデルを連結して形成される全体の形状が凸関数で近似される形状であるかどうかの整合性を評価する第３の値である。
【０２４０】
ここで、凸関数で近似される形状としての整合性とは隣接する部分領域モデルＸsとＸvのモデルパラメータの半径の差分および中心軸の最短距離がなるべく小さいことと、各部分領域モデルＸsの中心軸方向の変化方向のベクトルが同一平面上に存在し、かつ一定の方向であることを良いと評価する項目とする。
【０２４１】
上記の評価関数を最適にする結果は、部分領域モデル候補数の全数組み合わせにより理論上は求められるが、しかしながら実施形態同様に、組み合わせ爆発により求めることが困難である。
【０２４２】
そこで、本発明では、組み合わせの最適化を効率的に解く方法として信頼度伝搬法（Belief Propagation）を用いることとする。(非特許文献４)このアルゴリズムにより、最適解は保障されないが、効率的に解くことが可能である。
【０２４３】
整合性項Ｈ（Ｘs,Ｘv）である第３の値は、次式で算出される。
【０２４４】
【数６２】
2015232513.tif 000064
【０２４５】
Ｓｍｏｏｔｈ（Ｘs,Ｘv）項は、隣接する部分領域モデルの中心軸方向が平行に近いほど大きな値を出力し、また部分領域モデルの中心軸が隣接する部分領域モデルの中心位置に近いほど大きな値を出力する項目である。
【０２４６】
一方、Convec（Ｘs,Ｘv）項は、凸形状にするための制約項であり、注目する部分領域sと隣接する部分領域vのモデルそれぞれにおいて、隣接するモデル中心位置との相対的な位置の移動方向が同じであるかどうか判定する項目である。
【０２４７】
［第３の実施形態］
次に、本発明の第３の実施形態を説明する。第３の実施形態のたわみ推定装置１００の構成は第１の実施形態と同様であるので、その説明を省略する。
次に、第３の実施形態の作用を説明する。第３の実施形態の作用は、第１の実施形態の作用と同様の部分があるので、主として異なる部分のみ説明する。
【０２４８】
計測範囲が限られている被写体の点群を用いた場合、部分領域モデルで推定した中心位置には誤差が含まれやすい。そこで、図１４に示すように、ステップＳ３＿３＿２の大域的モデル推定処理（図７のステップＳ３＿３の処理に対応）の前に、大域的モデル推定部１２１は、ステップＳ３＿３＿１で柱状構造物のどの方向に計測点群が多く存在しているか調べ、その方向（代表分布方向）と部分領域モデルの第１のモデルパラメータにより表される柱状構造物の表面との交点（凸部交点）の位置を算出する。そして、大域的モデル推定部１２１は、算出した各部分領域の凸部交点に基づいて大域的モデルの第２のモデルパラメータを推定する。ただし、柱状構造物はたわんでいるため、表面に垂直な方向になるように変換した代表分布方向と部分領域モデルとの表面との交点（凸部交点）を用いる。
【０２４９】
中心位置と異なり、計測点群が存在する表面の推定精度は高いため、計測範囲が限られた被写体であっても、凸部交点を用いることで中心軸の推定精度低下を抑制できる。
【０２５０】
以下、図１３を参照しながら、説明する。
【０２５１】
柱状構造物Ｉの周辺点群について、ＸＹ平面上に射影する。射影後の点群について、角度方向で分割した方向別の点群の頻度分布を作成し、最も頻度が高い方向を代表分布方向とする。例えば、柱状構造物Iの中心軸ＺIを中心回転軸として、Ｘ軸を基準方向（０度）としてＮｄｅｇ度刻みの方向別の頻度分布を作成すればよい。仮に、Ｘ軸からｉ・Ｎｄｅｇ度方向の頻度が最も多かったときは、ｉ・Ｎｄｅｇ度の方向を代表分布方向ＶecIとすればよい。
【０２５２】
【数６３】
2015232513.tif 000065
【０２５３】
ただし、ＰＩは円周率を表すとする。
【０２５４】
次に、各部分領域モデルごとに、ＸＹに平行な平面（例えば、部分領域の中心位置を通る平面）で断面を作成する。中心軸から離れた位置に、代表分布方向とＸＹ平面に垂直な平面（計測面）を設定し、部分領域モデルの断面形状において最も計測面に近い位置を凸部交点の位置として求める。
【０２５５】
ここで、部分領域モデルの中心軸ＵをＸＹ平面に射影した方向をＱｘ、ＸＹ平面上でＱｘに垂直な方向をＱｙ、部分領域モデルの半径をｒ、部分領域モデルの中心軸とグローバル座標の成す角をθとすると、楕円断面の長軸をＱｘ、短軸をＱｙとした次式の楕円方程式により断面形状が求められる。
【０２５６】
【数６４】
2015232513.tif 000066
【０２５７】
最後に、計測面と上記断面形状が最も近い位置として、凸部交点を求める。代表分布方向ＶecIに垂直な直線(Ｑｘ,Ｑｙ平面と計測面の交線)Ｇとの最短距離として求められる。具体的には、代表分布方向ＶecIをローカル座標系Ｑｘ,Ｑｙ平面に射影した成分を（ＶecQx、ＶecQY）とすると、直線Ｇは次式で求められる。
【０２５８】
【数６５】
2015232513.tif 000067
【０２５９】
ただし、直線Ｇが楕円と直線Ｇが２点で交わらないために、Ｃには大きな値を入れればよい。本実施形態では、
【０２６０】
【数６６】
2015232513.tif 000068
【０２６１】
とした。
【０２６２】
凸部交点は、楕円状での接線方向と直線Ｇの傾き方向が平行になるときである。上記楕円方程式をＱｘについて微分して算出した傾きと直線Ｇの傾きが等しくなる。
【０２６３】
【数６７】
2015232513.tif 000069
【０２６４】
上記３つの式の連立方程式により、断面のローカル座標系での位置として
【０２６５】
【数６８】
2015232513.tif 000070
【０２６６】
として求められる。ただし、２組の解が求められるため、本実施形態では直線Ｇと近い方の組を正しい解とする。
【０２６７】
３次元の座標として部分領域ｓの凸部交点Ｔｓは、次式で求められる。
【０２６８】
【数６９】
2015232513.tif 000071
【０２６９】
部分領域ｓについて、Ｎｒａｎ個の第１のモデルパラメータの候補が登録されているが、それぞれについて独立して上記の計算により求めればよい。
【０２７０】
ステップＳ３＿３＿２の大域的モデル推定処理において、大域的モデル推定部１２１は、各部分領域の中心位置に代えて、ステップＳ３＿３＿１で求めた各部分領域の凸部交点を用いて大域的モデルの第２のモデルパラメータを推定する。その他の処理は第１の実施形態と同様である。
【０２７１】
［第４の実施形態］
次に、本発明の第４の実施形態を説明する。第４の実施形態のたわみ推定装置１００の構成は第１の実施形態と同様の部分であるので、同様の部分には同一の符号を付してその説明を省略し、異なる部分について説明する。図１６に示すように、幾何モデル当てはめ部１１４は、インライア推定部１１９を更に備えている。
【０２７２】
次に、第４の実施形態の作用を説明する。第２の実施形態の作用は、第１の実施形態の作用と同様の部分があるので、主として異なる部分のみ説明する。
【０２７３】
第４の実施形態では、図１８に示すように、第１の実施形態におけるステップ３＿１に対応するステップ３＿１＿１の次に、ステップＳ３＿１＿２で、インライア推定部１１９は、インライア推定処理を実行する。
【０２７４】
詳細には後述するが、インライア推定処理により、柱状構造物Ｉの周辺点群の中から、アウトライアとなる点群を除去することができるため、後処理である大域モデル推定精度を向上させることができる。
【０２７５】
なお、第４の実施形態の変形例として、図１７に示すように、幾何モデル当てはめ部１１４は、インライア推定処理（ステップＳ＿１＿２（図１８参照））を実行するインライア推定部１１９のみを備えることもできる。被計測表面積が上端部から下端部でも広く計測されているような状況では、例えば、柱状構造物を固定式のレーザーを用いて複数の方向から計測した結果を入力点群とした場合には、被計測表面積の不均一性が課題でなくなるため、本ステップＳ＿１＿２で推定したたわみ円筒モデルのパラメータベクトルを出力し、このパラメータベクトルを用いることで、たわみ量算出部１１５によりたわみ量を推定することも可能である。
【０２７６】
次に、インライア推定処理（ステップＳ＿１＿２（図１８参照））を説明する。
インライア推定部１１９は、たわみ円筒モデルを３次元点群に当てはめることで、焼きなまし法（ＳＡ）を用いてたわみ付き円筒モデルのパラメータベクトルを推定する。推定するたわみ円筒モデルのパラメータベクトルについて説明する（図１９）。グローバル座標系(Ｘ軸、Ｙ軸、Ｚ軸)に対して、柱状構造物の下端部の位置をローカル座標系の原点位置Ｔとする座標系(Ｕ軸、Ｖ軸、Ｗ軸)を設定する。Ｕ軸は柱状構造物の傾き方向とし、Ｖ軸はたわみ円筒モデルの中心軸が存在する平面を規定する軸とする。Ｗ軸はＵＶ平面に垂直な方向の軸とする。
【０２７７】
たわみ円筒モデルの中心軸は凸形状の曲線モデルとし、本実施形態ではローカル座標系の原点を通る２次関数Ｕ＝ａ・Ｖ2とする。柱状構造物の下端部での半径、すなわちローカル座標系の原点での半径の大きさをＲ、テーパの大きさをα（固定値）とする。本ステップにおいては、以上のパラメータをまとめたパラメータベクトル
【０２７８】
【数７０】
2015232513.tif 000072
【０２７９】
について、最適化する。図２０にインライア推定処理のフローチャートを示す。
【０２８０】
ステップＣ＿１において、インライア推定部１１９は、柱状構造物Iのパラメータ（下端部の中心位置、中心軸の方向、長さ、半径、凸関数の計数）および柱状構造物Ｉの周辺点群ｐmを入力する。ここで、個々の周辺点群を区別する番号をｍ、番号ｍの点の３次元位置をｐmと表す。
【０２８１】
【数７１】
2015232513.tif 000073
【０２８２】
ただし、ｘm、ｙm、ｚmはグローバル座標系でのＸ軸、Ｙ軸、Ｚ軸の成分とする。
【０２８３】
ステップＣ＿２において、インライア推定部１１９は、パラメータベクトルξを初期化し、そのときのスコアを次のように算出する。柱状構造物Ｉのパラメータを用いて、パラメータベクトル
【０２８４】
【数７２】
2015232513.tif 000074
【０２８５】
の値を決定する。具体的には、柱状構造物Ｉの下端部の中心位置をローカル座標系の３次元位置Ｔに代入し、柱状構造物Ｉの傾き方向をＵ軸、柱状構造物Ｉの半径を下端部での半径ＲIとする。残りの軸の方向および凸関数の係数ａについては乱数により決定してよい。ただし、Ｕ、Ｗ、Ｖはノルムが１．０であり、かつ互いに直交する方向になるようにしておく。直交化する方法は、グラム・シュミットの正規直交化法などを用いればよい。本発明においては、事前知識により係数ａの値の大きさには制限０＜ａ＜０．０１としておく。これは、下端部から１０ｍの高さにおいて、上端部のたわみ量が１．０(ｍ)以下になるという物理的意味である。
【０２８６】
次に、初期化したパラメータでのたわみ円筒モデルの当てはまりスコア値Scorebentを算出する。これは、たわみ円筒モデル表面上の点の数、もしくは点群の存在する面積の大きさと、各点の局所形状とたわみ円筒モデルの局所形状の類似度を用いて算出する。
本実施形態では、点の数および局所形状情報に基づくスコア値Scorebentの算出法を用いる。
【０２８７】
【数７３】
2015232513.tif 000075
【０２８８】
ただし、FuncBent(ｐm)は点ｐmにおけるたわみ円筒モデル表面との距離（図１２参照）を意味し、関数absは絶対値を出力する関数とする。LocalGeometry(ｐm)は、ｐm点とその近傍点群の局所形状情報、例えば法線や曲率半径、最小曲率方向、最大曲率方向と、点ｐmの位置におけるたわみ円筒モデルの局所形状情報が類似度合いを出力する関数とする。また、ＴＨbentとＴＨlocalは実験的に決めるパラメータであり、本実施形態では、ＴＨbent＝０．０２（ｍ）、ＴＨlocal＝０．８とした。
ここで、局所形状の類似度合いの算出方法は幾つかあるが、本実施形態では次式で求める。
【０２８９】
【数７４】
2015232513.tif 000076
【０２９０】
上記類似度の算出式の第一項は、点ｐmの位置での法線方向
【０２９１】
【数７５】
2015232513.tif 000077
と点ｐmから最も近いたわみ円筒モデル表面位置、つまりFuncBent(ｐm)＝０となる位置ｐ'mでの、たわみ円筒モデル表面の法線方向
【０２９２】
【数７６】
2015232513.tif 000078
の内積が大きいほど大きな類似度を出力し、第二項は点ｐmの位置での曲率半径curpmと位置ｐ'mにおけるたわみ円筒モデルの半径の大きさcurmodelの差分が小さい程大きな類似度を出力し、第三項は点ｐmの位置での最小曲率方向μmmと位置ｐ'mにおけるたわみ円筒モデルの最小曲率方向、つまり２次曲線の中心軸の接線方向μmodelの差分が小さい程大きな類似度を出力する。
【０２９３】
ｗL1、ｗL2、ｗL3は局所形状を重視する重み係数であり、実験的に決めるパラメータである。本稿ではｗL1、＝０．５、ｗL2、＝ｗL3＝０．２５とした。点ｐmと近傍点群の局所形状情報は、非特許文献５のような手法で求まる。たわみ円筒モデルの局所形状は、図１２で示すように、曲率半径はｒh、最小曲率方向は位置ｑhでの中心軸の接線方向（２次関数をｕ軸に対して微分した方向）、法線方向は位置ｑhを基準とした点ｐmへの相対ベクトル（ｐm−ｑh）の方向として求まる（非特許文献５）
【０２９４】
上記のスコア値Scorebentについて、最大スコアScorebent maxとし、以降の処理ではこの初期パラメータベクトルを
【数７７】
2015232513.tif 000079
として記憶する。
【０２９５】
インライア推定部１１９は、以下の処理（ステップＣ＿３〜Ｃ＿５）を閾値ＴＨ_k回繰り返すが、ステップＣ＿８で、繰り返し処理回数を識別するｋに１を代入する。なお、閾値ＴＨ_kは実験的に決めるパラメータであり、本実施形態においては、閾値ＴＨ_k＝１００００とした。
【０２９６】
ステップＣ＿３において、インライア推定部１１９は、パラメータベクトル
【０２９７】
【数７８】
2015232513.tif 000080
【０２９８】
の中から、パラメータをランダムに一つ選び、選択したパラメータの値を微小に変化させる。各パラメータの変動の仕方は任意でよいが、本実施形態においては以下のように、微小変動の値は、各パラメータごとに設定した閾値以内の値をランダムに発生させることとする。以降、繰り返し処理の回数をｋと表示する。
【０２９９】
（１）Ｕ軸を選択した場合、U軸まわりにV軸とW軸を微小回転量
【０３００】
【数７９】
2015232513.tif 000081
【０３０１】
で回転したローカル座標系を求める。
【０３０２】
【数８０】
2015232513.tif 000082
【０３０３】
とすると、ロドリゲスの公式により回転行列
【０３０４】
【数８１】
2015232513.tif 000083
【０３０５】
は次式で求められる。
【０３０６】
【数８２】
2015232513.tif 000084
【０３０７】
回転座標
【０３０８】
【数８３】
2015232513.tif 000085
【０３０９】
を用いて、新たなローカル座標系Ｕ’、Ｖ’、Ｗ’の軸は次式で求められる。
【０３１０】
【数８４】
2015232513.tif 000086
【０３１１】
ただし、微小回転量Δθの値は制限値±ＭＡＸΔθの中でランダムに決定を行う。本実施形態においては微小回転量Δθの制限値ＭＡＸΔθは０．１(度)とした。
【０３１２】
（２）Ｖ軸を選択した場合、Ｕ軸と同様にしてＶ軸まわりに微小回転したローカル座標系を求める。
【０３１３】
（３）Ｗ軸を選択した場合、Ｕ軸と同様にしてＷ軸まわりに微小回転したローカル座標系を求める。
【０３１４】
（４）ローカル座標系の原点位置
【０３１５】
【数８５】
2015232513.tif 000087
【０３１６】
を選択した場合、さらにランダムにＸ,Ｙ,Ｚ軸方向の一つを選択する。仮に、X軸方向を選択した場合には、次式により微小値
【０３１７】
【数８６】
2015232513.tif 000088
【０３１８】
を加える。
【０３１９】
新たな原点位置Ｔkは次式で求められる。
【０３２０】
【数８７】
2015232513.tif 000089
【０３２１】
ただし、微小値ΔＴの値は制限値±ＭＡＸTの中でランダムに決定を行う。本実施形態においては微小値ΔＴの制限値ＭＡＸTは０．０１(m)とした。
【０３２２】
同様にして、Ｙ軸、Ｚ軸を選択したときも、各軸の成分に微小値eを加えればよい。
【０３２３】
（５）凸関数の係数パラメータａを選択した場合には、微小値Δａ（−ＭＡＸΔａ＜Δａ＜ＭＡＸΔａ）を加えればよい。ただし、微小値Δaの制限値ＭＡＸΔａは０．００１とした。ｋ回目の係数パラメータａは次式で求められる。
【０３２４】
【数８８】
2015232513.tif 000090
【０３２５】
（６）凸関数の係数パラメータＲを選択した場合には、微小値ΔＲ（−ＭＡＸΔR＜ΔＲ＜ＭＡＸΔR）を加えればよい。ただし、微小値ΔRの制限値ＭＡＸΔRは０．０１（ｍ）とした。ｋ回目の半径パラメータＲは次式で求められる。
【０３２６】
【数８９】
2015232513.tif 000091
【０３２７】
ステップＣ＿４において、インライア推定部１１９は、ｋ回目のパラメータでのたわみ円筒モデルのスコア値を算出する。算出方法はステップＣ＿２と同様にして求められる。
【０３２８】
ステップＣ＿５において、ｋ回目のパラメータのスコア値によって、更新の判定を行う。
【０３２９】
ｋ回目のスコア値Scorebentkが、ｋ-１回目のスコア値Scorebentk-1よりも大きいときには、ｋ回目のスコア値Scorebentkにおけるパラメータは、ｋ-１回目のスコア値Scorebentk-1におけるパラメータよりも柱状構造物をより精度よく表現しているので、次のように更新を行う。
【０３３０】
【数９０】
2015232513.tif 000092
【０３３１】
更に、ｋ−１回までの最高スコア値よりも大きいときには、最大スコア値Scorebent_maxおよび最大スコア値のパラメータ
【０３３２】
【数９１】
2015232513.tif 000093
【０３３３】
も更新を行う。
【０３３４】
【数９２】
2015232513.tif 000094
【０３３５】
ｋ回目のスコア値Scorebentkが、ｋ−１回目のスコア値Scorebentk-1以下のときには、次式の確率的な判定により更新を行う。
【０３３６】
【数９３】
2015232513.tif 000095
【０３３７】
ただし、random(0,1.0)は、０から１．０の乱数を発生させる関数とし、Temperatureは繰り返し処理が進む程小さい値をとる関数である。例えば、
【０３３８】
【数９４】
2015232513.tif 000096
【０３３９】
で定義される。これは、繰り返し処理数kに比例する関数である。本実施形態では、Ｄａ＝０．８、Temperature0＝１００とした。上記更新判定式が満たされるときは、ｋ回目のスコア値Scorebentkおよびパラメータベクトルの更新を行う。
【０３４０】
【数９５】
2015232513.tif 000097
【０３４１】
ただし、最大スコア値と最大スコアのパラメータベクトルの更新は行わない。
【０３４２】
ステップＣ＿６において、インライア推定部１１９は、ｋがＴＨ＿ｋ以上か否かを判断することにより、繰り返し処理回数が閾値ＴＨ＿ｋ回数以上になったか否かを判断する。繰り返し処理回数が閾値ＴＨ＿ｋ回数以上の場合には、処理はステップＣ＿９に進む。繰り返し処理回数が閾値ＴＨ＿ｋ回数以上でない場合には、ステップＣ＿７において、インライア推定部１１９は、ｋ←ｋ＋１とする。処理はステップＣ＿３に戻る。
【０３４３】
ステップＣ＿９において、インライア推定部１１９は、繰り返し処理を終了し、このとき、k回までの処理の中で最大スコア値のときのパラメータである最大スコアパラメータ
【０３４４】
【数９６】
2015232513.tif 000098
【０３４５】
のときのスコアを加算したときの点群をインライアとみなして出力する。
【０３４６】
［実施形態の効果］
本実施形態により、3次元点群を用いて柱状構造物のたわみ推定を自動化でき、柱状構造物の保守点検の稼働コストを減らす効果がある。また、レーザースキャナで計測したときに、被計測範囲が狭いため計測ノイズの影響が大きい状況でも、中心軸のたわみ量の誤推定結果を出力することを抑制できる。特に、オクルージョンの影響により計測表面積が不均一な状況や、看板やケーブルなどの付属品が接している状況でも高精度にたわみ量を推定することが可能な技術である。
【０３４７】
［変形例］
なお、たわみ推定装置１００の各処理を実行するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、当該記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより、たわみ推定装置１００に係る上述した種々の処理を行ってもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものであってもよい。また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、フラッシュメモリ等の書き込み可能な不揮発性メモリ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。
【０３４８】
さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムが送信された場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリ（例えばＤＲＡＭ（Dynamic Random Access Memory））のように、一定時間プログラムを保持しているものも含むものとする。また、上記プログラムは、このプログラムを記憶装置等に格納したコンピュータシステムから、伝送媒体を介して、あるいは、伝送媒体中の伝送波により他のコンピュータシステムに伝送されてもよい。ここで、プログラムを伝送する「伝送媒体」は、インターネット等のネットワーク（通信網）や電話回線等の通信回線（通信線）のように情報を伝送する機能を有する媒体のことをいう。また、上記プログラムは、前述した機能の一部を実現するためのものであっても良い。さらに、前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるもの、いわゆる差分ファイル（差分プログラム）であっても良い。
【０３４９】
以上、本発明者によってなされた発明を、前記実施形態に基づき具体的に説明したが、本発明は、前記実施形態に限定されるものではなく、その要旨を逸脱しない範囲において種々変更可能であることは勿論である。
【符号の説明】
【０３５０】
１００    たわみ推定装置
１０１    被写体計測部
１０２    入力部
１０７    出力部
１０３    記憶部
１０４    柱状構造物検出部
１０５    たわみ推定処理部
１０６    構造物安全性判定部
１１０    ３次元点群記憶部
１１１    演算処理用パラメータ記憶部
１１２    柱状構造物パラメータ記憶部
１１２Ａ    入力処理部
１１３    取得部
１１４    幾何モデル当てはめ部
１１５    たわみ量算出部
１１６    出力処理部
１１７    繰り返し処理部
１１９    インライア推定部
１２０    部分領域モデル推定部
１２１    大域的モデル推定部
【図５】
2015232513.tif 000099
【図６】
2015232513.tif 000100
【図７】
2015232513.tif 000101
【図９】
2015232513.tif 000102
【図１１】
2015232513.tif 000103
【図１４】
2015232513.tif 000104
【図１６】
2015232513.tif 000105
【図１７】
2015232513.tif 000106
【図１８】
2015232513.tif 000107
【図２０】
2015232513.tif 000108
【図１】
2015232513.tif 000109
【図２】
2015232513.tif 000110
【図３】
2015232513.tif 000111
【図４】
2015232513.tif 000112
【図８】
2015232513.tif 000113
【図１０】
2015232513.tif 000114
【図１２】
2015232513.tif 000115
【図１３】
2015232513.tif 000116
【図１５】
2015232513.tif 000117
【図１９】
2015232513.tif 000118