JP2015232513A - たわみ推定装置、及びプログラム - Google Patents

Abstract

【課題】柱状構造物のたわみ量を精度よく推定する。【解決手段】柱状構造物の範囲内の各部分領域に対応する３次元点群から選択された複数の３次元点から、たわみがないとした各部分領域の形状を規定する複数の部分領域モデルの複数の第１のモデルパラメータを推定し（Ｓ３＿２）、複数の第１のモデルパラメータから、全体領域の３次元点群に対し、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータを推定し、各部分領域モデルが各部分領域に存在する尤度値及び大域的モデルと部分領域モデルとの整合性の評価値から、第２のモデルパラメータの評価値を複数回計算し（Ｓ３＿３〜Ｓ３＿５）、評価値を最適にする第２のモデルパラメータから、柱状構造物のたわみ量を算出する。よって、全体領域の点群形状を考慮しつつ、部分領域に最適な形状モデルを推定でき、柱状構造物のたわみ量を精度よく推定できる。【選択図】図７

Description

本発明は、たわみ推定装置、及びプログラムに係り、特に、柱状構造物のたわみ量を推定するたわみ推定装置、及びプログラムに関する。

近年、トンネルや橋などのインフラ構造物の劣化・老朽化が問題となり、自治体において維持管理やリスク評価が課題となっている。それと同時に、インフラ設備の数の多さが問題となり、構造物の安全性の評価の効率化やコスト削減を目的とした技術開発が求められている。特に、電柱や信号などのコンクリートポールのような柱状構造物の点検業務は、その数の多さから全数点検作業の稼働コストが大きい。

構造物のひずみ等の形状変形に関する情報（以下、たわみと呼ぶ）は倒壊につながる重要な安全指標の一つであり、従来はひずみゲージとよばれるセンサを構造物表面にとりつけることで計測を行ってきた。近年では、コストの問題やセンサ自体の耐久性を考慮して、カメラやレーザー等の光学的なセンサの開発が進み非接触な計測が可能となった。これにより、ひずみゲージ等のセンサ設置のための仮設足場などが不要となるため、短時間で複数の位置の情報を計測することが可能となった（非特許文献１）。また、柱状の構造物については、非特許文献２のような専用の製品活用により、より効率的な点検が可能となった。これにより、ある基準以上のたわみ量を有するものを不安全と判定し、その柱状構造物について保守点検を優先して行うことができる。

ただし、上記の非特許文献１においては、測定場所に基準点を設定して計測するときに、人手の作業が依然として必要である。非特許文献２のような専用装置を用いても、計測箇所の指定や柱状物体の中心軸位置を決定するためには人手による指示作業が必要であり、膨大な数の柱状構造物の点検作業における稼働コストは大きい。

一方、モバイルマッピングシステム（ＭＭＳ）と呼ばれるカメラやレーザースキャナを搭載した車が、街中を走行することで道路周辺の構造物表面の形状を計測できるシステムが開発され利用されつつある。このシステムは、ＧＰＳ（全地球測位システム）や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）、柱状構造物の総数をＮＩと表示する。柱状構造物パラメータには、柱状構造物の下端部の中心位置の３次元座標p_{I_bottom}、中心軸の方向

、構造物の長さＬ_I、半径Ｒ_I等の柱状構造物に関する形状情報がある。

ステップＳ８において、繰り返し処理部１１７は、Ｉに１を代入する。

ステップＳ２において、取得部１１３は、柱状構造物Ｉの中心軸から半径Ｒcut（＝Ｒ_I＋ΔＲ）以内に含まれる３次元点群を取得する。この半径ΔＲは演算処理用パラメータ記憶部１１１に登録されているパラメータである。以降では、このようなパラメータを「実験的に決めるパラメータである」と記載する。本実施形態において、ΔＲ＝０．０５(ｍ)とした。

ここで、取得した３次元点群の３次元点の数をＮ_m、各３次元点を区別する番号をｍ(∈１,２,..,Ｎ_m)で表示する。注目点ｍの位置座標をベクトルで表すと、ベクトルｐ_mは次式によりＸ軸、Ｙ軸、Ｚ軸の成分で構成される。ただし、右上の添え字の記号「Ｔ」は転置を意味とする。

ただし、Ｘm、Ｙm、ＺmはＸ，Ｙ，Ｚ軸の成分である。

柱状構造物Ｉの３次元点群は、下記の式を満たす３次元点として求められる。

ステップＳ３において、幾何モデル当てはめ部１１４は、ステップ２で取得した柱状構造物Ｉの３次元点群について、柱状構造物Ｉをモデル化したたわみ付き円筒モデルの当てはめ処理を行い、たわみ付き円筒モデルのモデルパラメータを推定する。ここで、たわみ付き円筒のモデルパラメータには、柱状構造物の下端部の中心位置ｑ_{_bottom}、柱状構造物の上端部の中心位置ｑ_{_top}、凸関数で表現した柱状構造物の中心軸、柱状構造物の下端部における接線方向

、柱状構造物の下端部での半径の大きさＲ、柱状構造物の長さＬがある。中心軸を表現する凸関数および詳細な処理内容は後述する。

ステップＳ４において、たわみ量算出部１１５は、たわみ量を算出する。ステップ３で求められた柱状構造物の上端部の中心位置ｑ_top、柱状構造物の下端部の中心位置ｑ_bottomにおける傾き方向

について、たわみ量ｄ_hは、上端部の中心位置ｑ_topから傾き方向

までの最短距離としては次式により求められる。

もしくは、単位長さあたりの変形量として、次式により求めてもよい。柱状構造物Ｉの長さをＬとすると次式で求められる。

ステップＳ５で、繰り返し処理部１１７は、Ｉが柱状構造物の総数ＮＩ以上か否かを判断することにより、全ての柱状構造物について処理が終了したか否かを判断する。全ての柱状構造物について処理が終了していない場合には、ステップＳ７で、繰り返し処理部１１７は、Ｉ←Ｉ＋１とする。処理はステップＳ２へ戻る。なお、記号「←」は代入を意味するとし、数式「Ｉ←Ｉ＋１」はＩの数を１つ増やしてＩに代入するという意味である。全ての柱状構造物について処理が終了した場合には、ステップＳ６で、出力処理部１１６は、たわみ推定処理部１０５で推定したたわみ量を柱状構造物パラメータ記憶部１０３（柱状構造物パラメータ記憶部１１２）へ出力する。これにより、柱状構造物パラメータ記憶部１１２に記憶される。

以下、ステップＳ３について詳細な説明をする。本ステップＳ３では、幾何モデル当てはめ部１１４は、たわみ付き円筒モデルのモデルパラメータの推定処理を行う。以降、たわみ付き円筒モデルとは、円筒モデルの中心軸を凸関数で表現した３次元幾何モデルとする。

本実施形態においては、詳細には後述するが、柱状構造物をモデル化したたわみ付き円筒モデルのモデルパラメータを、次式の部分領域モデル評価項と、大域的モデルＧと部分領域モデルの整合性評価項と、から構成される評価関数Ｅ（Ｘ_s,Ｇ）の最適化問題として、推定する。

ただし、Ｍ_s（ｓ∈１,２,３,...Ｎs）は、たわみがないと仮定した場合の後述する部分領域ｓの形状を規定する部分領域モデルを示し、部分領域モデルＭ_sは、第１のモデルパラメータで表わされる。第１のモデルパラメータには、部分領域モデルＭ_sの中心位置、部分領域モデルＭ_sの半径、部分領域モデルＭ_sの中心軸方向がある。

Ｇは、柱状構造物の中心軸を凸関数で表現した３次元モデルであり、柱状構造物周辺の全点群に対する大域的なモデルである。大域的モデルＧは具体的には、各部分領域モデルＭ_sを凸関数で近似することで生成される凸形状の曲線モデルであり、第２のモデルパラメータで表わされる。

関数Ｄ（Ｍ_s）は、部分領域モデルＭ_sが部分領域sに存在する尤度であり、部分領域sの第１のモデルパラメータにより形状が表される柱状構造物の表面に、部分領域sに含まれる３次元点群が存在する度合いに基づいて評価する第１の値である。

関数Ｈ（Ｍ_s,Ｇ）は大域的モデルＧの第２のモデルパラメータと部分領域モデルＭ_sの第１のモデルパラメータとの整合性を評価する第２の値である。ｋは２つの項目の重みを変える調整用係数である。ｋの値が大きいほど、滑らかに連結されるような部分領域モデルが選択されやすい。逆に、ｋ値が小さいほど、各部分領域モデルは各部分領域の形状に最適なパラメータとなりやすい。

ここで、整合性とは、部分領域モデルＭ_Sの第１のモデルパラメータと、大域的モデルＧの第２のモデルパラメータとの類似度を意味する。本実施形態において、実験的なパラメータkは、部分領域モデル評価項と大域的モデル整合性評価項の比率を同じにするためｋ＝Ｎ_s(部分領域の総数)とした。

各部分領域については、たわみの無い円筒形状のたわみの無し円筒モデルのパラメータとして十分近似できると考え、各部分領域について、部分領域モデルとして、たわみ無し円筒モデルを採用して、第１のモデルパラメータを推定する。図８に示すように、部分領域ごとに、候補となるたわみ無し円筒モデル群を推定する。各部分領域で推定した全てのたわみ無し円筒モデルの中心位置を凸関数で近似したときに上記評価関数を最大にするような組み合わせを求めれば、上記の評価関数を最適にする結果として、たわみ付き円筒モデルのモデルパラメータが求められる。

しかしながら、全ての組み合わせを処理することは現実的な時間では処理できない。例えば、図８の部分領域の数がＮ_ｓ＝１０個、各部分領域での候補モデル数を１００個とすると、部分領域モデルの全ての組み合わせは１００の１０乗となるため、計算コストが膨大となってしまうために求められない。

そこで、本発明では、部分領域モデルと大域的モデルのそれぞれを逐次的に推定する処理を繰り返し行うことで、評価関数Ｅの最適解を求める（図７）。具体的には、部分領域モデル評価項を最適化するときは、大域的モデルＧを固定として、その固定した大域モデルに対して最も良い部分領域モデルＭ_sをそれぞれの領域で求める。次に、大域的モデルＧを算出する際には、上記で求めた各部分領域モデルＭ_sを固定としたときの最大尤度となる凸関数モデルを求める。この処理を収束するまで繰り返す。

図７のステップＳ３＿１において、幾何モデル当てはめ部１１４は、柱状構造物Ｉのパラメータとその柱状構造物Ｉから半径Ｒ_cut以内の周辺点群Ｐ_ｍを入力する。

ステップ３＿１においては更に、幾何モデル当てはめ部１１４は、複数の部分領域を次のように設定する。

まず、図８に示すような重複する部分領域を設定する。３次元点群Ｐ_ｍに含まれる３次元点をｐ_ｍ（ｍ∈１,２,３,...Ｎ_ｍ）、３次元点群内の３次元点の総数をＮ_ｍ、柱状構造物Ｉの中心軸方向を

軸、中心軸に垂直な２つの軸を

、下端部の中心位置ｐ_{I_bottom}を原点としたローカル座標系を設定すると、図８に示すように３次元点群を中心軸方向

に複数の領域へと分割する。ローカル座標原点ｐ_{I_bottom}から、中心軸方向へ長さΔｓづつ区切った区間を設定すると、部分領域ＲＳ_s(ｓ∈１,２,３,...Ｎs)に含まれる点は以下の式を満たす。

ただし、Δｓは実験的に決めるパラメータであり、本実施形態ではΔｓ＝０．５(ｍ)とする。

ここで数学的な記号の意味を記載すると、記号「→」は方向、つまり長さが１のベクトルであること（２ノルムが１であること）を明示するための記号とし、記号

はベクトルの２ノルムの大きさとする。記号ベクトル間の記号「・」はベクトルの内積演算、記号「×」はベクトルの外積、スカラー量とベクトル間の記号「・」は掛け算を意味する。

上式の物理的な意味は、柱状構造物Ｉの周辺点について、中心軸方向

に射影した点の位置が（ｓ・Δｓ）から±Δｓの範囲に入る点であることを意味する。

幾何モデル当てはめ部１１４は、後述する処理（Ｓ３＿２、Ｓ３＿３）を閾値ＴＨ_cnt回繰り返すが、ステップＳ３＿６で、繰り返し処理の回数を識別するｃｎｔに１を代入する。繰り返し処理回数の閾値ＴＨ_cntは実験的に決めるパラメータであり、本実施形態ではＴＨ_cnt＝１００００とした。

ステップＳ３＿２において、幾何モデル当てはめ部１１４は、部分領域モデル推定処理を実行し、評価関数Ｅのモデル評価項である第１の値Ｄ（Ｍ_s)を算出する。

具体的には、幾何モデル当てはめ部１１４は、各部分領域ｓにおいてＲＡＮＳＡＣ処理（図９）によりたわみ無し円筒当てはめ処理を行い、最もスコア値の高いたわみ無し円筒モデル（部分領域モデル）の第１のモデルパラメータのスコア値との比率（０．０以上１．０以下）として評価関数Ｅのモデル評価項の第１の値Ｄ（Ｍ_s)を算出する。上位Ｎｒａｎ個について、部分領域sにおけるたわみ無し円筒モデルの第１のモデルパラメータの候補として登録を行う。

図９を参照して、部分領域モデルパラメータ推定処理(ＲＡＮＳＡＣによる円筒当てはめ処理）について説明する。ＲＡＮＳＡＣによるたわみ無し円筒当てはめ処理（ステップＳ３＿２）は、各部分領域について実行する。

ステップＡ＿１において部分領域モデル推定部１２０は、部分領域ｓ_kの３次元点群（３次元点の総数はＮ_ｓ）を入力する。各３次元点ｐ_Ｓはｎ（１,２,・・・Ｎ_ｓ）で識別される。

部分領域モデル推定部１２０は、後述するステップＡ＿３からステップＡ＿６の処理をＫ_loop回繰り返すが、ステップＡ＿２において、繰り返し処理の回数を識別するｋに１を代入する。なお、回数Ｋ_loopは実験的に決めるパラメータであり、本実施形態ではＫ_loop＝１０００とした。

ステップＡ＿３において部分領域モデル推定部１２０は、部分領域ｓ_kにおける３次元点群からランダムに２個の３次元点を選択（サンプリング）し、２個の３次元点間の方向を処理ｋ回目における中心軸の方向

とする。具体的には２個の３次元点の座標の差分をとり、ノルムを１にするように長さを正規化すればよい。

の記号「ｋ」はステップＡ＿２からステップＡ＿７の繰り返し処理における処理のｋ(∈1,2,..,Ｋ_loop)番目の回数を意味する。

ステップＡ＿４において、部分領域モデル推定部１２０は、部分領域ｓ_kにおける３次元点群から更にランダムにＮｒ個を選択し、中心軸の方向

に垂直でかつ部分領域ｓ_kの重心を通る平面へ射影した座標位置を求める。平面の基底ベクトルは、

に直交して、かつ重心を通る２つのベクトルであればよい。例えば、グローバル座標のＸ軸とＹ軸について、

方向に直交するように、グラム・シュミットの正規直交化法を用いて変換してもよい。この平面上での座標位置は、平面の基底の軸方向

と注目点jの位置ベクトルｐ_jとの内積で求められる。ランダムに選択した点の位置座標をｐ_j、ランダムに選択したＮｒ個の点の重心位置をｑ_g ^kとすると、平面での座標［ｘ’_j、ｙ’_j］は次式で求められる。

ここで、Ｎｒは実験的に決めるパラメータであり、本実施形態ではＮｒ＝５とした。

ステップＡ_５において、部分領域モデル推定部１２０は、部分領域ｓ_kの第１のモデルパラメータを求める。具体的には、選択したＮｒ個を上記重心を通る平面上に射影した点群について、最小二乗法により円の当てはめを行い、半径r_e ^kと平面上での中心位置の座標

を部分領域ＲＳ_Ｋの第１のモデルパラメータの候補として求める。具体的には、以下の誤差関数Ｅ_straightを最小化するように求めればよい。

ステップＡ＿６において、部分領域モデル推定部１２０は、モデル当てはめスコア値を算出する。これは、ステップＡ＿５で求めたたわみ無し円筒モデルのモデルパラメータの評価を行うことに相当する。評価方法は幾つか存在するが、例えば円筒表面上に存在する点群の数や円筒表面上に存在する点の範囲（面積）等が用いられる。

本実施形態では、点の数をモデルの評価基準となるスコア値として求めた例を示す。部分領域ｓ_kに含まれる点ｐ_sについて、下記のスコア値Score_subregionにより求められる。

ここで、ＴＨ_estは実験的に決めるパラメータであり、本実施形態ではＴＨ_est=０．０２とした。上記式のＦはスカラーの値であることを意味する。

ステップＡ＿７において、部分領域モデル推定部１２０は、ｋがＫ_loop以上か否かを判断することにより、繰り返し処理がＫ_loop回以上終了したか判定を行う。Ｋ_loop回以上の処理が終了していない場合には、ステップＡ＿９で、部分領域モデル推定部１２０は、ｋに１加算した値をｋに代入する（ｋ←ｋ＋１）。処理は、ステップＡ＿３に戻る。

繰り返し処理がＫ_loop回以上終了したと判定された場合には、処理はステップＡ＿８に進む。

ステップＡ＿８において、部分領域モデル推定部１２０は、最も高いスコア値からＮｒａｎ個のたわみ無し円筒モデルの第１のモデルパラメータの候補と当該Ｎｒａｎ個のスコア値とを出力する。このとき出力するパラメータはたわみ無し円筒モデルの中心位置、半径、中心軸の方向である。中心位置については、平面上の座標系(２次元ベクトル)からグローバル座標系へ変換をした位置座標ｑ_s(３次元ベクトル)として大域的モデル推定部１２１に出力する。部分領域sにおいて、上位e個目のパラメータについて、次式のように登録する。

ステップＡ＿８では更に、上位Ｎｒａｎ個の第１のモデルパラメータの候補について、評価関数Ｅのモデル評価項である第１の値Ｄ（Ｍ_s）を算出して、大域的モデル推定部１２１に出力する。上位一番目のスコア値Score_subregionの値をScore_subregion_maxとすると、上位Ｎｒａｎ個のうちの予め定めた上位e番目の部分領域モデル（たわみ無し円筒モデル）Ｍ_s ^eのスコア値をScore_subregion(e)とすると、部分領域モデル（たわみ無し円筒モデル）Ｍ_s ^eにおけるモデル評価項である第１の値Ｄ（Ｍ_s ^e）は次式で求められる。なお、第１の値は、部分領域モデルＭ_s ^eが、部分領域ｓ_kの３次元点群を精度よく表現しているかを評価する値である。

そして、大域的モデルＧについて、部分領域sにおける中心位置をｑ_G、その中心位置での接線方向を

、部分領域モデルＭ_s ^eの中心位置をｐ_s ^e、中心軸方向を

と表記すると、整合性評価項である第２の値は次式で求められる。

ただし、αは重みを調整する係数であり実験的に決めるパラメータであり、本実施形態においては、例えばα＝０．５とした。ここで、部分領域ｓにおける中心位置をｑ_Gは、大域的モデルの中心軸上の点において、各部分領域ｓの部分領域モデルの中心位置からの最短位置とする。

このとき繰り返し処理cnt回目の部分領域ｓにおける部分領域モデルＭ_sの第１のモデルパラメータは次式で求められる。

また、繰り返し処理が初期の場合、大域的モデルＧを考慮しない次式で求める。

ステップＡ＿８の処理が終了すると、処理は、図７のステップＳ３＿３に進む。

ステップＳ３＿３において、大域的モデル推定部１２１は、大域的モデル推定処理を実行する。本ステップＳ３＿３では、大域的モデル推定部１２１は、柱状構造物Ｉの周辺の点群を用いずに、ステップＳ３−２で求めた部分領域モデルの中心位置を通過するＮ次の凸関数をＲＡＮＳＡＣ処理により求める。本実施形態では２次の凸関数（２次曲線）を例に説明する。

凸関数当てはめの詳細な説明の前に、図１０に中心軸推定処理の概要を説明する。まず、大域的モデル推定部１２１は、図１０（Ａ）に示すように、ステップＳ３−２で推定したたわみ無し円筒モデルのモデルパラメータのうちの中心位置、中心軸の方向（u１、ｕ２）をランダムに２つの部分領域分選び、柱状構造物Ｉの形状を近似するたわみ付き円筒モデルの中心軸(凸関数)が存在する平面を推定する。当該平面は、ｕ軸、ｖ軸により定まる。次に、大域的モデル推定部１２１は、図１０（Ｂ）に示すように、この平面上に、全ての部分領域のたわみ無し円筒モデルのモデルパラメータのうちからランダムに数個選択したモデルパラメータにおける中心位置の点を射影し、平面上（ローカル座標系）での座標を求める。大域的モデル推定部１２１は、この射影した点群を用いて、凸関数のパラメータを算出する。凸関数のパラメータから、たわみ開始位置を算出する。

以下、ステップＳ３＿３を実現する方法として、ＲＡＮＳＡＣによる２次関数のパラメータ推定処理を説明する（図１１参照）。

ステップＢ＿１において、大域的モデル推定部１２１は、ステップＳ３＿２で出力した部分領域のたわみ無し円筒モデルのモデルパラメータを入力する。

大域的モデル推定部１２１は、後述するステップＢ＿３からステップＢ＿６の処理をＴ_loop回繰り返すが、ステップＢ＿２において、繰り返し処理の回数を識別するｋに１を代入する。なお、Ｔ_loopは実験的に決めるパラメータであり、本実施形態ではＴ_loop＝１０００とした。

ステップＢ_３において、大域的モデル推定部１２１は、２次関数が存在する平面を推定する。まず、大域的モデル推定部１２１は、図１０（Ａ）に示すように、柱状構造物Ｉの下端位置での接線方向ｕ軸を推定する。部分領域で推定したたわみ無し円筒モデルの第１のモデルパラメータを区別する番号をc(c∈1,2,3,..,Nc)、第１のモデルパラメータcの中心軸方向をＵ_ｃ、部分領域において推定した第１のモデルパラメータの数をＮｃとすると、全ての第１のモデルパラメータからランダムに選んだ２つのモデルパラメータ（c1とc2と記載）を用いて、以下の式で柱状構造物Ｉの下端位置における接線方向Ｕ_kを求める。

ただし、f_randは０〜１の値をとる乱数とする。上記式は、ランダムに選んだ２つの部分領域のたわみ無し円筒モデルのモデルパラメータc1とc2の中心軸方向（u_C1、u_Ｃ２）を用いて接線方向Ｕ_kを計算することを意味する。

また、後処理のために平面の基底ベクトルであるＵ_kの第三成分が正となるように向きを調整する。Ｕ_kの第三成分、つまりグローバル座標系でのＺ軸の値が負の値であるとき、Ｕ_kの向きを変更する。

たわみ方向ｖ_k軸は、Ｕ_c1からＵ_kの方向成分を、グラム・シュミットの直交化法により０にして、ノルムが１となるように正規化することで求める。

平面の中心位置は、選択した２点（ｃ１とｃ２）の円筒パラメータの中心位置（ｑ_c1,q_c2）の重心とする。

ここで、添え字の「k」は、繰り返し処理がｋ回目のときの推定結果であることを意味する。

ステップＢ＿４において、大域的モデル推定部１２１は、図１０（Ｂ）に示すように、中心軸の存在する平面をローカル座標系とし、全ての部分領域モデルの第１のモデルパラメータからランダムに選択したＮｊ個の部分領域モデルの第１のモデルパラメータにおける中心位置を、上記中心軸が存在する平面へ射影して、その平面で基底される空間（たわみ円筒ローカル座標系）での座標を求める。

以下、このたわみ円筒ローカル座標系の軸はｕ軸、ｖ軸で表現し、その座標値（スカラー値）は、それぞれｕとｖで表現する。また位置ベクトルの右上の添え字として記号「’」(ダッシュ)があるものは、ローカル座標に射影した点の位置ベクトルもしくはローカル座標系での座標値であることを明示するためにつけている。

部分領域で推定したたわみ無し円筒モデルのモデルパラメータについて、ランダムに選択したモデルパラメータを区別する番号をj、その部分領域モデルパラメータの中心位置をｐ_jとすると、大域的モデルのローカル座標系での座標

への変換式は次式で求められる。

ステップＢ＿５において、大域的モデル推定部１２１は、図１０（Ｂ）に示すように、中心軸を表現する凸関数のパラメータを求める。凸関数は２次関数でも３次関数やそれ以上の次元の数でもよいが、本実施形態では２次関数とする。経験的に、安全・不安全の評価を行う柱状構造物は、形状の変形量が微小であるため、ノイズに影響されにくい２次か３次などの低い次元の関数の方が好ましい。本実施形態では２次関数を例として示す。２次数の係数をa,b,cとすると、下記の誤差関数Ｆ_errorを最小化するように係数として求める。

ここで、誤差関数Ｆ_errorについて、最小二乗法を用いて係数a,b,cを求めてもよいし、非線形の数値解析手法を用いて求めてもよい。

本ステップＢ＿５では更に、大域的モデル推定部１２１は、たわみ円筒モデルの下端部（図１０（Ｃ）参照）での中心位置ｑ_bottomと上端部での中心位置ｑ_topを求める。部分領域モデル推定処理Ｓ３−２で求めた中心位置を、上記推定した中心軸へ射影することで求められる。

高さ方向において最も低い位置にある部分領域ｓ＝０の位置におけるたわみ無し円筒モデル

の中心位置を

、高さ方向において最も高い位置にある部分領域ｓ＝Ｎｓの位置でのたわみ無し円筒モデル

の中心位置を

とすると、次式で下端部の中心位置ｑ_bottomと上端部の中心位置ｑ_topが求められる。

以上により求めたパラメータを、繰り返し処理ｋ回目の大域的モデルパラメータＧ^kとする。
ステップＢ_６において、大域的モデル推定部１２１は、評価関数Ｅの値を算出する。まず、ステップＳ３−２（図７参照）で推定した各部分領域のたわみ無し円筒モデルの最も高いスコア値に対応するモデルパラメータにおける中心位置が、上記凸関数として推定した柱状構造物Ｉの中心軸上に存在するかどうか判定する。

中心軸上に存在するかどうかは、注目点から中心軸への最短距離により判定する。3次元空間中の曲線モデルと３次元点との最短距離を算出する方法は幾つかあるが、本実施形態では、大域的モデル推定部１２１は、図１２に示すように、ある３次元上の点ｐ_mからステップＢ_５で設定した凸関数の基底軸（Ｕ軸、Ｖ軸）へ射影した点ｐ’_mと凸関数の中心軸との最短位置ｑ_hとの距離から三平方の定理により算出する。射影点

は、たわみ円筒ローカル座標系の座標として次式により求められる。

まず、大域的モデル推定部１２１は、点ｐ_mをＵＶ平面へ射影した射影点ｐ’_mとステップＢ_５で求めた中心軸について、最短距離ｄｉｓｔ1を求める。

ここで、幾何的な性質から、点からある曲線までの最短距離の位置まで引いた線分と、最短距離の位置での接線は直交するという関係性がある。よって、最短位置

での接線方向

と射影点ｐ’_mから最短位置ｑ_hまでの方向

の内積は０となり、次式が成立する。

未知数（変数）はｕ’_hのみであるから、ｕ’_hの方程式として解くことができる。この方程式はニュートン法や修正ニュートン法などの数値解析手法で解くことが可能である。

ただし、この方程式は高次の方程式のため、一般的に複数の解が求められる。そのため、求められる複数の解について、射影点ｐ’_mから最も近い位置が最短位置ｑ_hの位置とする。

最短距離ｄｉｓｔ1は次式で求められる。

次に、大域的モデル推定部１２１は、点ｐ_mと射影した射影点ｐ’_mとの距離ｄｉｓｔ２を求める。

点ｐ_mと垂線の足の位置ｑ_hとの距離ｄｉｓｔは、幾何的な関係から次式で求められる。

ここで、上記では幾何的な関係から最短距離を求めたが、この方法以外にも例えば、ラグランジュの未定乗数法を用いて解いてもよい。制約条件は、最短位置の点が２次関数上に存在することとし、射影点ｐ’_mと凸関数曲線上の点において距離が最小となる解（極値）を求めればよい。

上記最短距離の算出式ｄｉｓｔ（ｐ_m）用いて、部分領域ｑ_sと大域的モデルの中心軸との最短位置ｑ_hとの距離が求められるため、整合性項である第２の値は次式により求められる。

ただし、

は、最短位置ｑ_hにおける中心軸の接線方向

であるから、

よって、整合性項の値は次式で求められる。

そして、大域的モデル推定部１２１は、評価関数Ｅの評価値を、次の通り求める。なお、右辺の第１項の値は、部分領域モデル推定部１２０から入力されている。

ステップＢ_７において、大域的モデル推定部１２１は、ｋがＴ_loop以上か否かを判断することにより、繰り返し処理の回数がＴ_loop回以上か否かを判断する。繰り返し処理の回数がＴ_loop回以上の場合には、処理はステップＢ＿８に進む。繰り返し処理の回数がＴ_loop回以上でない場合には、ステップＢ＿９で、大域的モデル推定部１２１は、ｋ←ｋ＋１とする。処理はステップＢ_３へ戻る。

ステップＢ_８において、大域的モデル推定部１２１は、評価関数Ｅが最も大きいときのたわみ付き円筒モデルのモデルパラメータを、繰り返し処理ｃｎｔ回目のたわみ付き円筒モデルのモデルパラメータＧ^cntとして、たわみ量算出部１１５に出力する。たわみ付き円筒モデルのモデルパラメータは、たわみ中心軸の存在するローカル座標系(Ｕ軸、Ｖ軸、Ｗ軸)、凸関数の係数（Ｎ次関数のパラメータ）、たわみ付き円筒モデルの下端部および上端部の中心位置、下端部位置での凸関数（大域的モデルのモデルパラメータ）の接線方向（ローカル座標系のＵ軸）である。

ステップＢ＿８の処理が終了すると、処理は、図７のステップＳ３＿４に進む。ステップＳ３＿４について、幾何モデル当てはめ部１１４は、ｃｎｔが閾値ＴＨ_cnt以上か否かを判断することにより、繰り返し処理の回数ｃｎｔが閾値ＴＨ_cnt以上か否かを判断する。繰り返し処理の回数ｃｎｔが閾値ＴＨ_cnt以上の場合に、処理はステップＳ３＿７に進む。繰り返し処理の回数ｃｎｔが閾値ＴＨ_cnt以上でない場合には、ステップＳ３＿５で、幾何モデル当てはめ部１１４は、ｃｎｔ←ｃｎｔ+１とし、処理はステップＳ３＿２に戻る。

ステップＳ３＿７で、幾何モデル当てはめ部１１４は、評価関数Ｅを最大にするときのたわみ付き円筒モデルのモデルパラメータを出力する。ここで、たわみ付き円筒モデルのモデルパラメータとは、大域的モデルの第２のモデルパラメータであり、たわみ中心軸の存在するローカル座標系(Ｕ軸、Ｖ軸、Ｗ軸)および凸関数の係数および、下端部および上端部の中心位置、下端部位置での凸関数の接線方向（ローカル座標系のＵ軸）である。

Ｗ軸はＵ軸およびＶ軸に対して垂直な軸とし、グラムシュミットの直交化法として唯一の３次元ベクトルとして求められる。

（構造物安全性判定部１０６の動作説明）
最後に、構造物安全性判定部１０６について説明する。構造物安全性判定部１０６は、たわみ推定処理部から出力された柱状構造物パラメータを用いて、各柱状構造物について安全か不安全な状態であるかの判定を行う。

各柱状構造物において、柱状構造物パラメータであるたわみ量ｄ_hにより以下の方程式を満たすとき不安全な状態であると判定する。

ここで、閾値ＴＨ_dhは実験的に決めるパラメータであり、本実施形態では、ＴＨ_dh＝０．２（ｍ）とした。

構造物安全性判定部１０６は、不安全と判定された柱状構造物の座標情報、及びたわみ量を含む警告信号を出力部１０７に出力する。

（出力部１０７の動作説明）
出力部１０７は、警告信号を構造物安全性判定部１０６から取得した場合には、不安全と判定された柱状構造物の座標情報、たわみ量をリスト化したものを、表示したり印刷したりする。

［第２の実施形態］
次に、本発明の第２の実施形態を説明する。第２の実施形態のたわみ推定装置１００の構成は第１の実施形態と同様であるので、その説明を省略する。
次に、第２の実施形態の作用を説明する。第２の実施形態の作用は、第１の実施形態の作用と同様の部分があるので、主として異なる部分のみ説明する。

柱状構造物について屈曲（ある高さで急激な曲がりが生じている状態）が存在する場合、２次関数や３次関数などの滑らかな凸関数では中心軸は表現精度が低下する可能性がある。

そのような状態のときには、凸関数という制限のみを加えたモデル評価関数Ｅを用いることで、中心軸の推定精度を向上することが可能である。本ステップ３において、次式の部分領域モデル評価項と大域的モデルとの整合性項から構成される損失関数Ｅの最小化問題とみなす。

ただし、Ｘ_s=(ｓ∈１,２,・・・Ｎ_ｓ）は部分領域ｓで推定した部分領域モデルであり、ＧはＸ_sにより規定される大域的モデルである。関数Ｄ（Ｘ_s)は、部分領域モデルＸ_sにより表される部分領域の形状の表面に各部分領域における点群が存在する度合いに基づいて評価する第１の値であり、関数Ｈ（Ｘ_s,Ｘ_v）は各部分領域モデルを連結して形成される全体の形状が凸関数で近似される形状であるかどうかの整合性を評価する第３の値である。

ここで、凸関数で近似される形状としての整合性とは隣接する部分領域モデルＸ_sとＸ_vのモデルパラメータの半径の差分および中心軸の最短距離がなるべく小さいことと、各部分領域モデルＸ_sの中心軸方向の変化方向のベクトルが同一平面上に存在し、かつ一定の方向であることを良いと評価する項目とする。

上記の評価関数を最適にする結果は、部分領域モデル候補数の全数組み合わせにより理論上は求められるが、しかしながら実施形態同様に、組み合わせ爆発により求めることが困難である。

そこで、本発明では、組み合わせの最適化を効率的に解く方法として信頼度伝搬法（Belief Propagation）を用いることとする。(非特許文献４)このアルゴリズムにより、最適解は保障されないが、効率的に解くことが可能である。

整合性項Ｈ（Ｘ_s,Ｘ_v）である第３の値は、次式で算出される。

Ｓｍｏｏｔｈ（Ｘ_s,Ｘ_v）項は、隣接する部分領域モデルの中心軸方向が平行に近いほど大きな値を出力し、また部分領域モデルの中心軸が隣接する部分領域モデルの中心位置に近いほど大きな値を出力する項目である。

一方、Convec（Ｘ_s,Ｘ_v）項は、凸形状にするための制約項であり、注目する部分領域sと隣接する部分領域vのモデルそれぞれにおいて、隣接するモデル中心位置との相対的な位置の移動方向が同じであるかどうか判定する項目である。

［第３の実施形態］
次に、本発明の第３の実施形態を説明する。第３の実施形態のたわみ推定装置１００の構成は第１の実施形態と同様であるので、その説明を省略する。
次に、第３の実施形態の作用を説明する。第３の実施形態の作用は、第１の実施形態の作用と同様の部分があるので、主として異なる部分のみ説明する。

計測範囲が限られている被写体の点群を用いた場合、部分領域モデルで推定した中心位置には誤差が含まれやすい。そこで、図１４に示すように、ステップＳ３＿３＿２の大域的モデル推定処理（図７のステップＳ３＿３の処理に対応）の前に、大域的モデル推定部１２１は、ステップＳ３＿３＿１で柱状構造物のどの方向に計測点群が多く存在しているか調べ、その方向（代表分布方向）と部分領域モデルの第１のモデルパラメータにより表される柱状構造物の表面との交点（凸部交点）の位置を算出する。そして、大域的モデル推定部１２１は、算出した各部分領域の凸部交点に基づいて大域的モデルの第２のモデルパラメータを推定する。ただし、柱状構造物はたわんでいるため、表面に垂直な方向になるように変換した代表分布方向と部分領域モデルとの表面との交点（凸部交点）を用いる。

中心位置と異なり、計測点群が存在する表面の推定精度は高いため、計測範囲が限られた被写体であっても、凸部交点を用いることで中心軸の推定精度低下を抑制できる。

以下、図１３を参照しながら、説明する。

柱状構造物Ｉの周辺点群について、ＸＹ平面上に射影する。射影後の点群について、角度方向で分割した方向別の点群の頻度分布を作成し、最も頻度が高い方向を代表分布方向とする。例えば、柱状構造物Iの中心軸Ｚ_Iを中心回転軸として、Ｘ軸を基準方向（０度）としてＮｄｅｇ度刻みの方向別の頻度分布を作成すればよい。仮に、Ｘ軸からｉ・Ｎｄｅｇ度方向の頻度が最も多かったときは、ｉ・Ｎｄｅｇ度の方向を代表分布方向Ｖec_Iとすればよい。

ただし、ＰＩは円周率を表すとする。

次に、各部分領域モデルごとに、ＸＹに平行な平面（例えば、部分領域の中心位置を通る平面）で断面を作成する。中心軸から離れた位置に、代表分布方向とＸＹ平面に垂直な平面（計測面）を設定し、部分領域モデルの断面形状において最も計測面に近い位置を凸部交点の位置として求める。

ここで、部分領域モデルの中心軸ＵをＸＹ平面に射影した方向をＱｘ、ＸＹ平面上でＱｘに垂直な方向をＱｙ、部分領域モデルの半径をｒ、部分領域モデルの中心軸とグローバル座標の成す角をθとすると、楕円断面の長軸をＱｘ、短軸をＱｙとした次式の楕円方程式により断面形状が求められる。

最後に、計測面と上記断面形状が最も近い位置として、凸部交点を求める。代表分布方向Ｖec_Iに垂直な直線(Ｑｘ,Ｑｙ平面と計測面の交線)Ｇとの最短距離として求められる。具体的には、代表分布方向Ｖec_Iをローカル座標系Ｑｘ,Ｑｙ平面に射影した成分を（Ｖec_Qx、Ｖec_QY）とすると、直線Ｇは次式で求められる。

ただし、直線Ｇが楕円と直線Ｇが２点で交わらないために、Ｃには大きな値を入れればよい。本実施形態では、

とした。

凸部交点は、楕円状での接線方向と直線Ｇの傾き方向が平行になるときである。上記楕円方程式をＱｘについて微分して算出した傾きと直線Ｇの傾きが等しくなる。

上記３つの式の連立方程式により、断面のローカル座標系での位置として

として求められる。ただし、２組の解が求められるため、本実施形態では直線Ｇと近い方の組を正しい解とする。

３次元の座標として部分領域ｓの凸部交点Ｔｓは、次式で求められる。

部分領域ｓについて、Ｎｒａｎ個の第１のモデルパラメータの候補が登録されているが、それぞれについて独立して上記の計算により求めればよい。

ステップＳ３＿３＿２の大域的モデル推定処理において、大域的モデル推定部１２１は、各部分領域の中心位置に代えて、ステップＳ３＿３＿１で求めた各部分領域の凸部交点を用いて大域的モデルの第２のモデルパラメータを推定する。その他の処理は第１の実施形態と同様である。

［第４の実施形態］
次に、本発明の第４の実施形態を説明する。第４の実施形態のたわみ推定装置１００の構成は第１の実施形態と同様の部分であるので、同様の部分には同一の符号を付してその説明を省略し、異なる部分について説明する。図１６に示すように、幾何モデル当てはめ部１１４は、インライア推定部１１９を更に備えている。

次に、第４の実施形態の作用を説明する。第２の実施形態の作用は、第１の実施形態の作用と同様の部分があるので、主として異なる部分のみ説明する。

第４の実施形態では、図１８に示すように、第１の実施形態におけるステップ３＿１に対応するステップ３＿１＿１の次に、ステップＳ３＿１＿２で、インライア推定部１１９は、インライア推定処理を実行する。

詳細には後述するが、インライア推定処理により、柱状構造物Ｉの周辺点群の中から、アウトライアとなる点群を除去することができるため、後処理である大域モデル推定精度を向上させることができる。

なお、第４の実施形態の変形例として、図１７に示すように、幾何モデル当てはめ部１１４は、インライア推定処理（ステップＳ＿１＿２（図１８参照））を実行するインライア推定部１１９のみを備えることもできる。被計測表面積が上端部から下端部でも広く計測されているような状況では、例えば、柱状構造物を固定式のレーザーを用いて複数の方向から計測した結果を入力点群とした場合には、被計測表面積の不均一性が課題でなくなるため、本ステップＳ＿１＿２で推定したたわみ円筒モデルのパラメータベクトルを出力し、このパラメータベクトルを用いることで、たわみ量算出部１１５によりたわみ量を推定することも可能である。

次に、インライア推定処理（ステップＳ＿１＿２（図１８参照））を説明する。
インライア推定部１１９は、たわみ円筒モデルを３次元点群に当てはめることで、焼きなまし法（ＳＡ）を用いてたわみ付き円筒モデルのパラメータベクトルを推定する。推定するたわみ円筒モデルのパラメータベクトルについて説明する（図１９）。グローバル座標系(Ｘ軸、Ｙ軸、Ｚ軸)に対して、柱状構造物の下端部の位置をローカル座標系の原点位置Ｔとする座標系(Ｕ軸、Ｖ軸、Ｗ軸)を設定する。Ｕ軸は柱状構造物の傾き方向とし、Ｖ軸はたわみ円筒モデルの中心軸が存在する平面を規定する軸とする。Ｗ軸はＵＶ平面に垂直な方向の軸とする。

たわみ円筒モデルの中心軸は凸形状の曲線モデルとし、本実施形態ではローカル座標系の原点を通る２次関数Ｕ＝ａ・Ｖ²とする。柱状構造物の下端部での半径、すなわちローカル座標系の原点での半径の大きさをＲ、テーパの大きさをα（固定値）とする。本ステップにおいては、以上のパラメータをまとめたパラメータベクトル

について、最適化する。図２０にインライア推定処理のフローチャートを示す。

ステップＣ＿１において、インライア推定部１１９は、柱状構造物Iのパラメータ（下端部の中心位置、中心軸の方向、長さ、半径、凸関数の計数）および柱状構造物Ｉの周辺点群ｐ_mを入力する。ここで、個々の周辺点群を区別する番号をｍ、番号ｍの点の３次元位置をｐ_mと表す。

ただし、ｘ_m、ｙ_m、ｚ_mはグローバル座標系でのＸ軸、Ｙ軸、Ｚ軸の成分とする。

ステップＣ＿２において、インライア推定部１１９は、パラメータベクトルξを初期化し、そのときのスコアを次のように算出する。柱状構造物Ｉのパラメータを用いて、パラメータベクトル

の値を決定する。具体的には、柱状構造物Ｉの下端部の中心位置をローカル座標系の３次元位置Ｔに代入し、柱状構造物Ｉの傾き方向をＵ軸、柱状構造物Ｉの半径を下端部での半径Ｒ_Iとする。残りの軸の方向および凸関数の係数ａについては乱数により決定してよい。ただし、Ｕ、Ｗ、Ｖはノルムが１．０であり、かつ互いに直交する方向になるようにしておく。直交化する方法は、グラム・シュミットの正規直交化法などを用いればよい。本発明においては、事前知識により係数ａの値の大きさには制限０＜ａ＜０．０１としておく。これは、下端部から１０ｍの高さにおいて、上端部のたわみ量が１．０(ｍ)以下になるという物理的意味である。

次に、初期化したパラメータでのたわみ円筒モデルの当てはまりスコア値Score_bentを算出する。これは、たわみ円筒モデル表面上の点の数、もしくは点群の存在する面積の大きさと、各点の局所形状とたわみ円筒モデルの局所形状の類似度を用いて算出する。
本実施形態では、点の数および局所形状情報に基づくスコア値Score_bentの算出法を用いる。

ただし、FuncBent(ｐ_m)は点ｐ_mにおけるたわみ円筒モデル表面との距離（図１２参照）を意味し、関数absは絶対値を出力する関数とする。LocalGeometry(ｐ_m)は、ｐ_m点とその近傍点群の局所形状情報、例えば法線や曲率半径、最小曲率方向、最大曲率方向と、点ｐ_mの位置におけるたわみ円筒モデルの局所形状情報が類似度合いを出力する関数とする。また、ＴＨ_bentとＴＨ_localは実験的に決めるパラメータであり、本実施形態では、ＴＨ_bent＝０．０２（ｍ）、ＴＨ_local＝０．８とした。
ここで、局所形状の類似度合いの算出方法は幾つかあるが、本実施形態では次式で求める。

上記類似度の算出式の第一項は、点ｐ_mの位置での法線方向

と点ｐ_mから最も近いたわみ円筒モデル表面位置、つまりFuncBent(ｐ_m)＝０となる位置ｐ'_mでの、たわみ円筒モデル表面の法線方向

の内積が大きいほど大きな類似度を出力し、第二項は点ｐ_mの位置での曲率半径cur_pmと位置ｐ'_mにおけるたわみ円筒モデルの半径の大きさcur_modelの差分が小さい程大きな類似度を出力し、第三項は点ｐ_mの位置での最小曲率方向μ_mmと位置ｐ'_mにおけるたわみ円筒モデルの最小曲率方向、つまり２次曲線の中心軸の接線方向μ_modelの差分が小さい程大きな類似度を出力する。

ｗ_L1、ｗ_L2、ｗ_L3は局所形状を重視する重み係数であり、実験的に決めるパラメータである。本稿ではｗ_L1、＝０．５、ｗ_L2、＝ｗ_L3＝０．２５とした。点ｐ_mと近傍点群の局所形状情報は、非特許文献５のような手法で求まる。たわみ円筒モデルの局所形状は、図１２で示すように、曲率半径はｒ_h、最小曲率方向は位置ｑ_hでの中心軸の接線方向（２次関数をｕ軸に対して微分した方向）、法線方向は位置ｑ_hを基準とした点ｐ_mへの相対ベクトル（ｐ_m−ｑ_h）の方向として求まる（非特許文献５）

上記のスコア値Score_bentについて、最大スコアScore_{bent max}とし、以降の処理ではこの初期パラメータベクトルを
として記憶する。

インライア推定部１１９は、以下の処理（ステップＣ＿３〜Ｃ＿５）を閾値ＴＨ_k回繰り返すが、ステップＣ＿８で、繰り返し処理回数を識別するｋに１を代入する。なお、閾値ＴＨ_kは実験的に決めるパラメータであり、本実施形態においては、閾値ＴＨ_k＝１００００とした。

ステップＣ＿３において、インライア推定部１１９は、パラメータベクトル

の中から、パラメータをランダムに一つ選び、選択したパラメータの値を微小に変化させる。各パラメータの変動の仕方は任意でよいが、本実施形態においては以下のように、微小変動の値は、各パラメータごとに設定した閾値以内の値をランダムに発生させることとする。以降、繰り返し処理の回数をｋと表示する。

（１）Ｕ軸を選択した場合、U軸まわりにV軸とW軸を微小回転量

で回転したローカル座標系を求める。

とすると、ロドリゲスの公式により回転行列

は次式で求められる。

回転座標

を用いて、新たなローカル座標系Ｕ’、Ｖ’、Ｗ’の軸は次式で求められる。

ただし、微小回転量Δθの値は制限値±ＭＡＸ_Δθの中でランダムに決定を行う。本実施形態においては微小回転量Δθの制限値ＭＡＸ_Δθは０．１(度)とした。

（２）Ｖ軸を選択した場合、Ｕ軸と同様にしてＶ軸まわりに微小回転したローカル座標系を求める。

（３）Ｗ軸を選択した場合、Ｕ軸と同様にしてＷ軸まわりに微小回転したローカル座標系を求める。

（４）ローカル座標系の原点位置

を選択した場合、さらにランダムにＸ,Ｙ,Ｚ軸方向の一つを選択する。仮に、X軸方向を選択した場合には、次式により微小値

を加える。

新たな原点位置Ｔ^kは次式で求められる。

ただし、微小値ΔＴの値は制限値±ＭＡＸ_Tの中でランダムに決定を行う。本実施形態においては微小値ΔＴの制限値ＭＡＸ_Tは０．０１(m)とした。

同様にして、Ｙ軸、Ｚ軸を選択したときも、各軸の成分に微小値eを加えればよい。

（５）凸関数の係数パラメータａを選択した場合には、微小値Δａ（−ＭＡＸ_Δａ＜Δａ＜ＭＡＸ_Δａ）を加えればよい。ただし、微小値Δaの制限値ＭＡＸ_Δａは０．００１とした。ｋ回目の係数パラメータａは次式で求められる。

（６）凸関数の係数パラメータＲを選択した場合には、微小値ΔＲ（−ＭＡＸ_ΔR＜ΔＲ＜ＭＡＸ_ΔR）を加えればよい。ただし、微小値ΔRの制限値ＭＡＸ_ΔRは０．０１（ｍ）とした。ｋ回目の半径パラメータＲは次式で求められる。

ステップＣ＿４において、インライア推定部１１９は、ｋ回目のパラメータでのたわみ円筒モデルのスコア値を算出する。算出方法はステップＣ＿２と同様にして求められる。

ステップＣ＿５において、ｋ回目のパラメータのスコア値によって、更新の判定を行う。

ｋ回目のスコア値Score_bent ^kが、ｋ-１回目のスコア値Score_bent ^k-1よりも大きいときには、ｋ回目のスコア値Score_bent ^kにおけるパラメータは、ｋ-１回目のスコア値Score_bent ^k-1におけるパラメータよりも柱状構造物をより精度よく表現しているので、次のように更新を行う。

更に、ｋ−１回までの最高スコア値よりも大きいときには、最大スコア値Score_{bent_max}および最大スコア値のパラメータ

も更新を行う。

ｋ回目のスコア値Score_bent ^kが、ｋ−１回目のスコア値Score_bent ^k-1以下のときには、次式の確率的な判定により更新を行う。

ただし、random(0,1.0)は、０から１．０の乱数を発生させる関数とし、Temperatureは繰り返し処理が進む程小さい値をとる関数である。例えば、

で定義される。これは、繰り返し処理数kに比例する関数である。本実施形態では、Ｄａ＝０．８、Temperature⁰＝１００とした。上記更新判定式が満たされるときは、ｋ回目のスコア値Score_bent ^kおよびパラメータベクトルの更新を行う。

ただし、最大スコア値と最大スコアのパラメータベクトルの更新は行わない。

ステップＣ＿６において、インライア推定部１１９は、ｋがＴＨ＿ｋ以上か否かを判断することにより、繰り返し処理回数が閾値ＴＨ＿ｋ回数以上になったか否かを判断する。繰り返し処理回数が閾値ＴＨ＿ｋ回数以上の場合には、処理はステップＣ＿９に進む。繰り返し処理回数が閾値ＴＨ＿ｋ回数以上でない場合には、ステップＣ＿７において、インライア推定部１１９は、ｋ←ｋ＋１とする。処理はステップＣ＿３に戻る。

ステップＣ＿９において、インライア推定部１１９は、繰り返し処理を終了し、このとき、k回までの処理の中で最大スコア値のときのパラメータである最大スコアパラメータ

のときのスコアを加算したときの点群をインライアとみなして出力する。

［実施形態の効果］
本実施形態により、3次元点群を用いて柱状構造物のたわみ推定を自動化でき、柱状構造物の保守点検の稼働コストを減らす効果がある。また、レーザースキャナで計測したときに、被計測範囲が狭いため計測ノイズの影響が大きい状況でも、中心軸のたわみ量の誤推定結果を出力することを抑制できる。特に、オクルージョンの影響により計測表面積が不均一な状況や、看板やケーブルなどの付属品が接している状況でも高精度にたわみ量を推定することが可能な技術である。

［変形例］
なお、たわみ推定装置１００の各処理を実行するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、当該記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することにより、たわみ推定装置１００に係る上述した種々の処理を行ってもよい。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものであってもよい。また、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、フラッシュメモリ等の書き込み可能な不揮発性メモリ、ＣＤ−ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。

さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムが送信された場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリ（例えばＤＲＡＭ（Dynamic Random Access Memory））のように、一定時間プログラムを保持しているものも含むものとする。また、上記プログラムは、このプログラムを記憶装置等に格納したコンピュータシステムから、伝送媒体を介して、あるいは、伝送媒体中の伝送波により他のコンピュータシステムに伝送されてもよい。ここで、プログラムを伝送する「伝送媒体」は、インターネット等のネットワーク（通信網）や電話回線等の通信回線（通信線）のように情報を伝送する機能を有する媒体のことをいう。また、上記プログラムは、前述した機能の一部を実現するためのものであっても良い。さらに、前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるもの、いわゆる差分ファイル（差分プログラム）であっても良い。

以上、本発明者によってなされた発明を、前記実施形態に基づき具体的に説明したが、本発明は、前記実施形態に限定されるものではなく、その要旨を逸脱しない範囲において種々変更可能であることは勿論である。

１００ たわみ推定装置
１０１ 被写体計測部
１０２ 入力部
１０７ 出力部
１０３ 記憶部
１０４ 柱状構造物検出部
１０５ たわみ推定処理部
１０６ 構造物安全性判定部
１１０ ３次元点群記憶部
１１１ 演算処理用パラメータ記憶部
１１２ 柱状構造物パラメータ記憶部
１１２Ａ 入力処理部
１１３ 取得部
１１４ 幾何モデル当てはめ部
１１５ たわみ量算出部
１１６ 出力処理部
１１７ 繰り返し処理部
１１９ インライア推定部
１２０ 部分領域モデル推定部
１２１ 大域的モデル推定部

Claims (6)

1. 柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定し、
   前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定すること、及び
   前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記第１のモデルパラメータと、推定された第２のモデルパラメータとの整合性を評価する前記第２の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算すること
   を複数回行う幾何モデル当てはめ部と、
   前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、
   を含むたわみ推定装置。
2. 柱状構造物の位置に対応する範囲内に、各々前記柱状構造物の中心軸方向の隣と一部重複するように複数の部分領域を設定し、
   前記複数の部分領域の各々に対し、前記範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群のうち、前記部分領域に対応する３次元点群から選択された複数の３次元点に基づいて、前記部分領域にたわみがないと仮定した場合の前記部分領域の形状を規定する部分領域モデルの第１のモデルパラメータの候補を繰り返し推定し、繰り返し推定された前記第１のモデルパラメータの候補の各々についての、前記第１のモデルパラメータの候補により形状が表される柱状構造物の表面に前記３次元点群が存在する度合いに基づいて評価する第１の値と、前記第１のモデルパラメータの候補と、前記柱状構造物の全体領域の３次元点群に対して、柱状構造物の中心軸を凸関数で表現した大域的モデルの第２のモデルパラメータとの整合性を評価する第２の値とに基づいて、前記第１のモデルパラメータの候補から、前記部分領域の第１のモデルパラメータを推定すること、及び
   前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータに基づいて、前記第２のモデルパラメータを推定し、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータの各々についての、前記第１の値と、前記複数の部分領域の各々に対して推定された前記第１のモデルパラメータにより表わされる柱状構造物の全体形状が、凸関数で近似される形状であるかどうかの整合性を評価する第３の値とに基づいて、前記第２のモデルパラメータを評価する評価値を計算すること
   を複数回行う幾何モデル当てはめ部と、
   前記計算された複数の評価値の値に応じて選択された第２のモデルパラメータに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、
   を含むたわみ推定装置。
3. 前記幾何モデル当てはめ部は、
   前記複数の部分領域の各々に対し、前記部分領域に対応する３次元点群を、前記柱状構造物の高さ方向に垂直な平面に投影して得た複数の点の内、最も多くの３次元点が投影された点に、前記柱状構造物の中心軸から向かう代表分布方向を求め、
   前記複数の部分領域の各々について、前記部分領域の中心軸上の点から前記代表分布方向に向かう直線と、前記部分領域に対して推定された前記第１のモデルパラメータにより表される柱状構造物の表面との交点を求め、
   前記複数の部分領域の各々について求められた交点に基づいて、前記第２のモデルパラメータを推定する
   請求項１又は２記載のたわみ推定装置。
4. 前記幾何モデル当てはめ部は、
   前記柱状構造物の中心軸を凸関数で表現した３次元モデルを規定するパラメータベクトルのうちから選択されたパラメータを複数回変化させる毎に、前記パラメータを変化させた前記パラメータベクトルにより形状が表される柱状構造物の表面に前記範囲内の３次元点群で、かつ局所形状が前記柱状構造物と類似している３次元点群が存在する度合いを示すスコア値を求め、
   前記範囲内の３次元点群から、前記求められた最大のスコア値を求めるときに前記柱状構造物の表面に存在すると判断された３次元点群を抽出し、
   前記抽出した３次元点群に基づいて、前記複数の部分領域の各々に対して前記第１のモデルパラメータを推定する
   請求項１〜３の何れか１項に記載のたわみ推定装置。
5. 柱状構造物の位置に対応する範囲内に位置し、物体の表面上の位置を計測した計測結果である複数の位置を表す３次元点群に基づいて、前記柱状構造物の中心軸を凸関数で表現した３次元モデルを規定するパラメータベクトルのうちから選択されたパラメータを複数回変化させる毎に、前記パラメータを変化させた前記パラメータベクトルにより形状が表される柱状構造物の表面に前記範囲内の３次元点群で、かつ局所形状が前記柱状構造物と類似している３次元点群が存在する度合いを示すスコア値を求め、
   前記求められた最大のスコア値に対応する前記パラメータベクトルを出力する
   インライア推定部と、
   前記出力されたパラメータベクトルに基づいて、前記柱状構造物のたわみ量を算出するたわみ量算出部と、
   を含むたわみ推定装置。
6. 請求項１〜請求項５の何れか１項記載のたわみ推定装置の各部としてコンピュータを機能させるためのプログラム。