(57)【特許請求の範囲】
【請求項１】
スパース性を有する外力が付与される弾性体における視認可能な第１領域と視認不能な第２領域のうち、前記第１領域の第１変位の情報を取得する取得部と、
前記弾性体を表す所定のメッシュモデルにおいて、前記第１領域の第１変位と、前記第１領域の第１外力及び前記第２領域の第２外力との間に成立する有限要素法の関係式に、ｌ１再構成を適用して最適化問題に変換することにより、前記第１外力及び前記第２外力の双方を推定する推定部と、を備えている、外力推定装置。
【請求項２】
前記推定部は、
推定された前記第２外力を用いて、前記第２領域の第２変位を推定する、請求項１に記載の外力推定装置。
【請求項３】
前記弾性体に付与された既知の外力と、前記推定部が推定した前記第１外力との差に基づいて、前記メッシュモデルの剛性情報を修正する修正部を更に備えている、請求項１又は２に記載の外力推定装置。
【請求項４】
コンピュータが、スパース性を有する外力が付与される弾性体における視認可能な第１領域と視認不能な第２領域のうち、前記第１領域の第１変位の情報を取得するステップと、
前記コンピュータが、前記弾性体を表す所定のメッシュモデルにおいて、前記第１領域の第１変位と、前記第１領域の第１外力及び前記第２領域の第２外力との間に成立する有限要素法の関係式に、ｌ１再構成を適用して最適化問題に変換することにより、前記第１外力及び前記第２外力の双方を推定するステップとを含む、外力推定方法。
【請求項５】
スパース性を有する外力が付与される弾性体における視認可能な第１領域と視認不能な第２領域のうち、前記第１領域の第１変位の情報を取得する機能、及び、
前記弾性体を表す所定のメッシュモデルにおいて、前記第１領域の第１変位と、前記第１領域の第１外力及び前記第２領域の第２外力との間に成立する有限要素法の関係式に、ｌ１再構成を適用して最適化問題に変換することにより、前記第１外力及び前記第２外力の双方を推定する機能を、コンピュータに実現させるコンピュータプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、外力推定装置及び方法、並びにコンピュータプログラムに関する。
【背景技術】
【０００２】
近年、患者への侵襲を少なくし術後のＱＯＬ向上等を図ることを目的として、体に形成した小さな穴から内視鏡カメラや術具等を挿入し、内視鏡カメラで術野を確認しながら術具により手術を行う内視鏡下手術が数多く実施されている。また、最近では、遠隔からの操縦によってロボットに術具を操作させるロボット手術も実施されてきている。
【０００３】
しかし、内視鏡下手術やロボット手術においては、術野が狭いために手術対象となる臓器の全体像を視認することができない。また、ロボット手術においては、鉗子を用いて臓器に力を加えたとしてもその反力が術者へは提示されないので、術者は、狭い術野の中で臓器の変形を確認しながら臓器にかかる力を想像しなければならない。そのため、鉗子から臓器に過剰な力が加わって臓器や血管を傷つけてしまうおそれがあり、術具の操作には熟練した技能が求められる。
このような実情に鑑み、鉗子に加わる力を計測し、手術の技能分析に役立てる研究が行われている（例えば、非特許文献１参照）。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】吉田健志、外６名、「腹腔鏡手術用鉗子先端に加わる作用力の計測と剥離操作における技能分析」、生体医工学、Ｖｏｌ．４８（２０１０）、Ｎｏ．１、Ｐ２５−３２
【発明の概要】
【発明が解決しようとする課題】
【０００５】
ロボット手術による術中に、非特許文献１に記載の技術を用いて鉗子に加わる力を計測することができれば、臓器に加わる力を客観的に把握でき、臓器等の損傷を未然に防ぐことができると考えられるが、術野の範囲外で鉗子を操作することは困難であるため、術野の範囲外において臓器にどのような力が加わっているのかを把握することはできない。また、術野の範囲内で鉗子により臓器を引っ張った場合、術野の範囲外で他の組織に繋がった部分にも力が加わると考えられるが、そのような力も当然に把握することができない。
【０００６】
そこで、本発明は、弾性体の局所的な変位情報から、弾性体の全体に加わる外力を推定することができる外力推定装置及び方法、並びにコンピュータプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
（１）  本発明に係る外力推定装置は、スパース性を有する外力が付与される弾性体における視認可能な第１領域と視認不能な第２領域のうち、前記第１領域の第１変位の情報を取得する取得部と、
前記弾性体を表す所定のメッシュモデルにおいて、前記第１領域の第１変位と、前記第１領域の第１外力及び前記第２領域の第２外力との間に成立する有限要素法の関係式に、ｌ１再構成を適用して最適化問題に変換することにより、前記第１外力及び前記第２外力の双方を推定する推定部と、を備えたものである。
【０００８】
上記の本発明は、物体の構造解析等に多用されている有限要素法を用いて、弾性体の変位と弾性体の剛性情報とから弾性体に付与される外力を推定する。本来、有限要素法において、弾性体の全ての要素に付与される外力を取得するには、全ての要素における変位の情報が必要である。しかし、内視鏡下手術やロボット手術等では、術野の範囲外から変位の情報を取得することは困難である。
【０００９】
本出願の発明者は、手術対象となる臓器等の弾性体に鉗子等の小さな術具によって力を付与すると、術具が直接的に接する部分や、他の組織に繋がる部分において外力が付与され、外力が付与される領域の体積は、弾性体全体の体積よりも十分に小さく、弾性体に付与される外力にはスパース性があることに気づいた。そのため、弾性体の多くの要素に付与される外力をゼロと仮定でき、少ない変位の情報からでも非ゼロの外力を求めることが可能であると考えた。そして、このような考えを実現するため、少ない変位の情報を用いて有限要素法による関係式を構築するとともに、当該関係式にｌ１−再構成を適用し最適化問題に変換するという本発明を完成するに到った。このような発明によって、弾性体の一部の領域の変位の情報から、簡単な計算で弾性体全体の外力を推定することが可能となる。
【００１０】
（２）  前記推定部は、推定された前記第２外力を用いて、前記第２領域の第２変位を推定することが好ましい。
上記のように推定部は、弾性体の一部の領域である第１領域の変位を用いて、弾性体全体に付与される第１外力及び第２外力を推定するが、その推定された第２外力をさらに用いることによって、前記関係式から弾性体の第２領域の第２変位をも推定することができる。したがって、例えば内視鏡カメラによって視認できない領域も含めて弾性体全体の変位を把握することができる。
【００１１】
（３）  外力推定装置は、前記弾性体に付与された既知の外力と、前記推定部が推定した前記第１外力との差に基づいて、前記メッシュモデルの剛性情報を修正する修正部を更に備えていることが好ましい。
剛性情報は、例えば弾性体のヤング率やポアソン比等から算出することが可能であるが、弾性体が手術対象となる臓器等の場合には、その臓器から直接的にこれらのパラメータを取得するのは困難であるため、その臓器の標準的な剛性情報を採用することが考えられる。しかしながら、標準的な剛性情報と実際の剛性情報との相違によって、推定部によって推定された外力と実際に付与された外力とにずれが生じるおそれもある。本発明の外力推定装置は、推定部によって推定された外力と、既知の外力とを比較し、両者の差に基づいて剛性情報を修正することで、より正確に外力を推定することが可能となる。
【００１２】
（４）  本発明の外力推定方法は、コンピュータが、スパース性を有する外力が付与される弾性体における視認可能な第１領域と視認不能な第２領域のうち、前記第１領域の第１変位の情報を取得するステップと、
前記コンピュータが、前記弾性体を表す所定のメッシュモデルにおいて、前記第１領域の第１変位と、前記第１領域の第１外力及び前記第２領域の第２外力との間に成立する有限要素法の関係式に、ｌ１再構成を適用して最適化問題に変換することにより、前記第１外力及び前記第２外力の双方を推定するステップとを含むものである。
【００１３】
（５）  本発明のコンピュータプログラムは、スパース性を有する外力が付与される弾性体における視認可能な第１領域と視認不能な第２領域のうち、前記第１領域の第１変位の情報を取得する機能、及び、前記弾性体を表す所定のメッシュモデルにおいて、前記第１領域の第１変位と、前記第１領域の第１外力及び前記第２領域の第２外力との間に成立する有限要素法の関係式に、ｌ１再構成を適用して最適化問題に変換することにより、前記第１外力及び前記第２外力の双方を推定する機能を、コンピュータに実現させるものである。
【発明の効果】
【００１４】
本発明によれば、弾性体の一部の領域における変位情報から、弾性体の全体に付与される外力を推定することができる。
【図面の簡単な説明】
【００１５】
【図１】本発明の実施形態に係る外力推定装置の構成図である。
【図２】弾性体を説明する図である。
【図３】表示装置に表示する推定結果の一例を示す図である。
【図４】有限要素法で用いる対象物のメッシュモデルの一例を示す図である。
【図５】評価実験に用いる肝臓のモデルを示す図である。
【図６】評価実験に用いる肝臓のモデルに推定対象となる外力を付与した状態を示す図である。
【図７】評価実験において、１領域に外力を加えた場合の推定結果を示す図である。
【図８】評価実験において、２領域に外力を加えた場合の推定結果を示す図である。
【図９】評価実験において、１領域に外力を加えた場合の推定誤差を示すグラフである。
【図１０】評価実験において、２領域に外力を加えた場合の推定誤差を示すグラフである。
【発明を実施するための形態】
【００１６】
＜１．実施形態＞
以下、図面を参照して本発明の実施形態を詳細に説明する。
［外力推定システムの全体構成］
図１は、本発明の実施形態に係る外力推定システムの説明図である。
本実施形態の外力推定システム１０は、弾性体からなる対象物の画像を撮影するカメラ（撮像装置）１１と、カメラ１１により撮影された画像の信号が入力され、所定の処理を行うことによって対象物に付与される外力を推定する外力推定装置１２と、外力推定装置１２の出力結果等を表示する表示装置１３とを備えている。
【００１７】
この外力推定システム１０は、例えば、内視鏡下手術やロボット手術の支援のために用いることができる。
カメラ１１は、ＣＣＤやＣＭＯＳ等の撮像素子を備えたものであり、内視鏡に取り付けられ、患者の体腔内に挿入されることによって、対象物としての臓器等の術部を撮影する。カメラ１１は、１つであってもよいし、複数であってもよい。
【００１８】
外力推定装置１２は、カメラ１１によって撮影された術部の画像から対象物に付与される外力を推定するものである。外力の情報として、外力の位置、大きさ、向きが推定される。また、外力推定装置１２は、カメラ１１によって撮影された部分だけでなく撮影されていない部分に付与される外力をも推定する。例えば、図２に示すように、対象物Ｗの一部の領域Ａ１は、カメラ１１の画像から視認できる状態であるが、他の領域Ａ２は、カメラ１１によって撮影されず、視認できない状態にあるものとする。この場合、外力推定装置１２は、視認領域Ａ１に付与される外力Ｆ１だけでなく、非視認領域Ａ２に付与される外力Ｆ２も推定する。推定された外力の情報は、表示装置１３に表示される。
【００１９】
なお、手術中、臓器には、鉗子等の小さな術具によって部分的に外力が付与されることが多く、また、術具が直接的に作用する部分以外は、他の組織と繋がっている部分だけに外力が付与されると考えられる。したがって、臓器に対して外力が付与される領域の体積は、臓器全体の体積よりも十分に小さいと考えられ、臓器に付与される外力には「スパース性」があるといえる。外力推定装置１２は、このような外力のスパース性に着目して、臓器の一部の領域の画像から外力を推定するものとなっている。外力推定装置１２の具体的な構成と、外力推定のための処理アルゴリズムについては後述する。
【００２０】
表示装置１３は、術者又は手術スタッフ（以下、「術者等」ともいう）の視界に配置されたモニター、又は、術者等が装着したヘッドマウントディスプレイ等とすることができる。そのため、術者は、表示装置１３に表示された外力の情報をもとに臓器に付与される外力を客観的に把握することができ、臓器に対して過剰な力を付与しないように配慮して手術を行うことが可能になる。
【００２１】
図３は、表示装置１３への表示例を示す。対象物を模した図に対して外力が付与されている部分Ｒと外力の大きさを色彩の変化又は濃淡の変化等によって示したものである。対象物の図には、カメラを介して視認できる領域Ａ１と視認できない領域Ａ２とを識別できるように、２つの領域Ａ１，Ａ２が区画線Ｃによって区画されている。このような図を表示装置１３に表示させることによって、術者等は外力が生じている部位と外力の大きさとを一目で把握することができる。
【００２２】
［外力推定装置の具体的構成］
外力推定装置１２は、例えばパーソナルコンピュータから構成され、ＣＰＵ等の演算部と、ＲＯＭ、ＲＡＭ、ハードディスク等の記憶部と、各種入出力インタフェース等を有している。カメラ１１によって撮影された画像の情報は、入出力インタフェースを介して外力推定装置１２に入力される。そして、外力推定装置１２は、記憶部にインストールされたプログラムを演算部が実行することによって、入力された情報を用いて外力を推定し、推定結果を表示装置１３に表示させる。記憶部には、外力の推定に用いられる種々の情報も記憶されている。
【００２３】
外力推定装置１２は、その機能構成として、カメラ１１から入力された画像情報を処理し、撮影された領域の変位についての情報を取得する取得部２１と、取得部２１において取得された変位についての情報を用いて対象物に付与される外力を推定する推定部２２と、推定部２２の処理に用いられる所定の定数を修正する修正部２３とを有している。
【００２４】
取得部２１は、時系列で取得された画像において特徴点の位置を変化を追うことによって変位情報を取得する。この取得部２１における具体的な処理は、従来公知の技術を用いることができる。例えば、手術中に取得したカメラ画像によって弾性体表面の形状変化を追う技術（Nazim Haouchine等、「Image-guided simulation of heterogeneous tissue deformation for augmented reality during hepatic surgery」、Mixed and Augmented Reality(ISMAR), 2013 IEEE International Symposium、(IEEE, 2013)、P199-208.）等を用いることができる。
また、取得部２１は、例えば、術前に撮影されたＣＴ（Ｃｏｍｐｕｔｅｄ  Ｔｏｍｏｇｒａｐｈｙ）画像を用いて３次元の初期モデルを生成し、カメラによって撮影された臓器の画像に基づいて初期モデルの形状を順次更新していき、その更新の度に変位を取得する方法を用いてもよい。
【００２５】
推定部２２は、取得部２１によって取得された変位の情報に基づいて、弾性体に付与される外力を推定する。以下、推定部２２における具体的な処理アルゴリズムについて説明する。
［推定部２２における処理アルゴリズム］
推定部２２は、有限要素法を用いることによって外力を推定する。有限要素法では、弾性体を多数の要素に分割し、各要素の頂点における変位と外力との関係として、以下の式（１）を用いる。
【００２６】
【数１】
0006648883.tif 000002
【００２７】
ここで、ｆは、外力ベクトル、Ｋは、剛性マトリクス（剛性情報）、ｕは、変位ベクトルである。剛性マトリクスＫは、弾性体のヤング率、ポアソン比等から求められる定数である。
また、式（１）は、次の式（２）のように変形することができる。
【００２８】
【数２】
0006648883.tif 000003
【００２９】
図４に、弾性体を複数の要素に分割したメッシュモデルの一例を示す。この弾性体は立方体形状であり、図４（ａ）は、外力が付与される前の状態、図４（ｂ）は、外力が付与された状態を示す。メッシュモデルの頂点ｔ１（黒い点）は固定点、頂点ｔ２（濃いグレーの点）は自由点、頂点ｔ３（グレーの点）は観測点、頂点ｔ４（薄いグレーの点）は外力が加わる点を示す。頂点ｔ４には、矢印Ｓ１で示すように推定対象となる外力が付与されている。観測点ｔ３は、図２に示すように、外部から視認できる領域Ａ１に配置される頂点に相当する。
【００３０】
本実施形態では、既知条件として、図４のメッシュモデルの観測点ｔ３における変位ベクトルｕｏと、剛性マトリクスＫとを与え、メッシュモデルの全頂点に加わるスパース性のある外力ベクトルｆを推定する。変位ベクトルｕｏは、外力推定装置１２の取得部２１によって取得される。剛性マトリクスＫは、定数として予め与えられる。
【００３１】
メッシュモデルの頂点を、視認できる頂点ｏと視認できない頂点ｉに分類し、マトリクスＬをブロック分割すると、式（２）は、次の式（３）のように表すことができる。
【００３２】
【数３】
0006648883.tif 000004
【００３３】
ただし、ｕｏ及びｆｏは、視認できる頂点ｏ（観測点ｔ３）における変位ベクトル及び外力ベクトルであり、ｕｉ及びｆｉは、視認できない頂点ｉにおける変位ベクトル及び外力ベクトルである。
そして、視認できる頂点ｏにおける変位ベクトルｕｏが既知量であるため、次の式（４）のように、当該変位ベクトルｕｏのみから外力ベクトルｆｏ、ｆｉを求めることを考える。
【００３４】
【数４】
0006648883.tif 000005
【００３５】
式（４）は、既知量ｕｏの成分数よりも、未知量ｆｏ，ｆｉの成分数の方が多い多元連立一次方程式である。そのため、何らかの拘束条件なしに一意に解くことはできない不良設定問題となる。本実施形態では、式（４）に対してｌ１再構成を適用することによって、外力ベクトルｆの推定を可能にする。なお、ｌ１再構成の一般的な考え方については、後の＜３．ｌ１再構成について＞の項目において述べることとする。
【００３６】
まず、外力ベクトルｆの１−ノルムを式（５）のように定義すると、上記の式（４）は、次の式（６）のような最適化問題を解くことになる。
【００３７】
【数５】
0006648883.tif 000006
【００３８】
【数６】
0006648883.tif 000007
【００３９】
そして、上記最適化問題は、以下のような線形計画問題として記述することができる。
具体的には、式（６）は、ｆのｉ番目の要素ｆｉについて、非負整数ｆ＋ｉ，ｆ−ｉを用いると、式（７）のように記述することができる。
【００４０】
【数７】
0006648883.tif 000008
【００４１】
そして、ｆ＝ｆ＋−ｆ−より、次の式（８）を得ることができる。
【００４２】
【数８】
0006648883.tif 000009
【００４３】
なお、式（７）及び式（８）を導出する詳細な手順は、後述の＜３．ｌ１再構成について＞の項目において説明する。
【００４４】
推定部２２は、上記のような線形計画問題を公知の解法を用いて解き、外力ベクトルｆを推定する。公知の解法として、例えばシンプレックス法又は内点法を採用することができる。推定された外力は、図３に示すように図式化され、表示装置１３に表示される。
【００４５】
また、推定部２２は、推定された外力ベクトルｆを用いることによって、視認できない頂点における変位ベクトルｕｉを上記の式（２）により求める。これにより、既知量である変位ベクトルｕｏに加え、未知量である変位ベクトルｕｉをも求めることができ、弾性体全体の外力だけでなく、変位の状態をも把握することができる。
【００４６】
外力推定装置１２の修正部２３は、定数である剛性マトリクスＫを適正なかたちに修正するものである。体内の臓器を手術するような場合、その臓器自体の剛性マトリクスＫを直接的に求めることは困難であるため、外力推定装置１２の推定部２２は、臓器の標準的な剛性マトリクスを用いる。しかしながら、標準的な剛性マトリクスと、実際の臓器の剛性マトリクスとは厳密に一致しないので、推定される外力ベクトルｆや変位ベクトルｕｉが実際の外力や変位から若干ずれる可能性がある。
【００４７】
このような問題を解決するため、修正部２３は、例えば、手術前に臓器に既知の外力を付与し、その外力による変位ｕｏをカメラ１１の画像から取得し、この変位ｕｏを用いて上記の式（１）〜（８）を用いて外力ベクトルｆを推定する。そして、次の式（９）を用いることによって、既知の外力ｆ’と、推定された外力ベクトルｆ（＝Ｋｕ）と差分を最小化するようなＫの値を求める。これにより、実際の対象物により適した剛性マトリクスＫを求めることができ、より正確な外力及び未知の変位の推定が可能となる。
【００４８】
【数９】
0006648883.tif 000010
【００４９】
以上、説明したように、本実施形態の外力推定装置においては、弾性体にスパースな外力が付与される場合に、有限要素法の関係式にｌ１再構成を適用し、最適化問題（線形計画問題）に変換することで、弾性体の一部の領域における変位を用いて、弾性体全体の外力を推定する。
したがって、内視鏡下手術やロボット手術において、狭い術野内で臓器の一部しか視認することができず、その一部の変位しか取得できない場合であっても、視認できる臓器の一部に加え、視認できない臓器の他の部分についても外力を推定することができる。そのため、鉗子等の術具によって臓器に力を付与する場合に、視認できる部分だけでなく視認できない部分においても損傷が生じないように配慮して手術を行うことができる。
【００５０】
［変形例］
本発明の外力推定装置は、上述したように、内視鏡下手術やロボット手術において手術対象となる臓器等の外力を推定するために好適に利用することができるが、これに限定されるものではなく、あらゆる物体に付与される外力を推定するために好適に利用することができる。
例えば、本発明は、木造やコンクリート造等の構造物や建築物に対して付与されるスパースな外力を推定するために用いることができる。この場合、視認できない部分の外力を推定することで、構造物等の破損、倒壊等を予測することが可能となる。また、地質等、表面しか観測できないような対象物において、内部に付与される外力を推定することができる。また、本発明は、生体の細胞等のような微小な物体にも適用することができる。
また、対象物の画像を撮影する撮像装置は、上述したようなカメラに限らず、対象物の一部を画像化できるものであれば特に限定されない。例えば、手術中であっても撮影が可能な、超音波検査装置、Ｘ線撮影装置、ＣＴ撮影装置、ＭＲＩ（Ｍａｇｎｅｔｉｃ  Ｒｅｓｏｎａｎｃｅ  Ｉｍａｇｉｎｇ）撮影装置等を採用することもできる。また、対象物が構造物や建築物の場合には、レンジセンサ（デプスセンサ）を採用することもできる。
もっとも、特許請求の範囲に記載の「弾性体」には、上記の木造やコンクリート造等の構造物や建築物は含まれない。
【００５１】
＜２．評価実験＞
本出願の発明者は、上記において説明した処理アルゴリズムの有効性を確認するため、評価実験を行った。この評価実験では、対象物として、図５に示すような肝臓のメッシュモデルを用いた。図５において、領域ｒ１内の頂点ｔ１は固定点であり、その他の頂点ｔ２は自由点である。固定点ｔ１は、肝臓に下大動脈が接続されていることを想定して設定している。このメッシュモデルは、固定点ｔ１を除いて表面に１７８個、内部に８５個の頂点があり、合計２６３個の頂点を有する。図５には、３次元の座標系も併せて示す。
【００５２】
第１の評価実験として、図６（ａ）に示すように、メッシュモデルの視認できる１つの領域に、鉗子で引っ張り力を与えた場合を想定して−ｚ方向に外力を加えた。具体的には、領域ｒ４内の９個の頂点ｔ４に外力を加えた。その外力を矢印Ｓ１で示す。
また、第２の評価実験として、図６（ｂ）に示すように、メッシュモデルの視認できる１つの領域に、鉗子で引っ張り力を与えた場合を想定して−ｚ方向の外力を加え、それによって、視認できない他の領域にも外力が加わったと想定して、＋ｘ方向に外力を加えた。これらの外力を矢印Ｓ１で示す。−ｚ方向の外力は領域ｒ４内の９個の頂点ｔ４に加え、＋ｘ方向の外力は領域ｒ５内の７個の頂点ｔ５に加えた。
【００５３】
そして、第１，第２の評価実験の双方とも、視認できる領域において頂点の変位を取得し、その変位から対象物全体の外力を推定した。そして、推定された外力と、実際に与えた外力とを比較し、推定精度を評価した。
【００５４】
図７は、第１の評価実験における推定結果の１例を示す図である。図７には、実際に加えた外力を黒い矢印Ｓ１で示し、推定された外力をグレーの矢印Ｓ２で示している。この推定結果は、固定点ｔ１を除くメッシュモデルの２６３頂点のうち、ｚ座標の小さい頂点から順に５３点（領域ｒ３内の頂点ｔ３）を観測した例を示す。
【００５５】
図７（ｂ）は、図７（ａ）のＡ部を拡大したものである。図７（ｂ）から明らかなように、推定された外力Ｓ２と、実際に加えた外力Ｓ１とは、概ね一致していることがわかる。
また、９個の頂点に実際に加えた外力の大きさは５Ｎであったのに対して、その頂点において推定された外力の大きさは、３．７〜６．３Ｎであった。それ以外の頂点の外力は、これよりも１０〜１０２のオーダーで小さい値であった。よって、実際に外力を加えた頂点について、実際に近い大きさの外力を高精度に推定することができた。
【００５６】
図８は、第２の評価実験における推定結果の１例を示す図である。図８には、実際に加えた外力を黒い矢印Ｓ１で示し、推定された外力をグレーの矢印Ｓ２で示している。この推定結果は、固定点ｔ１を除くメッシュモデルの２６３頂点のうち、ｚ座標の小さい頂点から順に７９点（領域ｒ３内の頂点ｔ３）を観測した例を示す。
【００５７】
図８（ｂ）（ｃ）は、それぞれ図８（ａ）のＢ部及びＣ部を拡大したものである。図８（ｂ）（ｃ）から明らかなように、推定された外力Ｓ２と、実際に加えた外力Ｓ１とは、概ね一致していることがわかる。
また、１６個の頂点に実際に加えた外力の大きさは５Ｎであったのに対して、その頂点において推定された外力の大きさは、３．９〜５．３Ｎであった。それ以外の頂点で推定された外力は、これよりも１０〜１０２のオーダーで小さい値であった。よって、実際に外力を加えた頂点について、実際に近い大きさの外力を高精度に推定することができた。
【００５８】
図９（ａ）は、第１の評価実験でメッシュモデルの表面のみを観測した場合について、表面の頂点１７８点のうち、観測点をｚ座標の小さいものから順に９点〜１７８点まで増加させた場合の推定誤差を示している。観測点を９点以上としたのは、外力を加えた点数以上とするためである。
同様に、図９（ｂ）は、第１の評価実験でメッシュモデルの内部も含めて観測した場合について、表面及び内部の頂点２６３点のうち、観測点をｚ座標の小さいものから順に９点〜２３４点まで増加させた場合の推定誤差を示している。
【００５９】
ここで用いる推定誤差は、与えた外力と推定外力との差の二乗平均平方根（Ｒｏｏｔ  Ｍｅａｎ  Ｓｑｕａｒｅ；以下、「ＲＭＳ」という）とし、次の式（１０）により求めた。
【００６０】
【数１０】
0006648883.tif 000011
【００６１】
ｆｘｉ，ｆｙｉ，ｆｚｉは、実際に加えた外力ｆの第ｉ要素ｆｉについてのｘｙｚ成分である。ｆ’ｘｉ，ｆ’ｙｉ，ｆ’ｚｉは、推定された外力ｆ’の第ｉ要素ｆ’ｉについてのｘｙｚ成分である。
Ｎは、観測する頂点をメッシュモデルの表面のみとした場合、メッシュモデルの全頂点のうち、固定点と内部の点を除いた頂点の数であり、観測する頂点にメッシュモデルの内部を含めた場合、メッシュモデルの全頂点のうち、固定点を除いた頂点の数である。
【００６２】
ＲＭＳは、加えた外力と推定された外力の要素の相違が生じている場合や、加えた外力ｆと推定された外力ｆ’のｘｙｚ成分ごとの大きさのずれが生じている場合に、増加する。言い換えると、正しく推定されるほど、ＲＭＳの値は小さくなる。
【００６３】
図９（ａ）（ｂ）に示すように、観測する頂点数に関わらず、ほぼＲＭＳは１以下の値となった。また、観測する頂点数が多くなるほど推定誤差ＲＭＳが小さくなることが分かる。
【００６４】
図１０（ａ）は、第２の評価実験でメッシュモデルの表面のみを観測した場合について、表面の頂点１７８点のうち、観測点をｚ座標の小さいものから順に１６点〜１７８点まで増加させた場合の推定誤差を示している。観測点を１６点以上としたのは、外力を加えた点数以上とするためである。
【００６５】
同様に、図１０（ｂ）は、第２の評価実験でメッシュモデルの内部も含めて観測した場合について、表面及び内部の頂点２６３点のうち、観測点をｚ座標の小さいものから順に１６点〜２３４点まで増加させた場合の推定誤差を示している。
【００６６】
図１０（ａ）（ｂ）に示すように、観測する頂点数に関わらず、ほぼＲＭＳは１以下の値となった。また、図１０（ａ）に示すように、観測する頂点数が多くなるほど、概ね推定誤差ＲＭＳが小さくなることが分かる。特に、観測点が５８点（３３％）以上の場合、ＲＭＳは安定して小さい値をとっていることが分かる。
【００６７】
＜３．ｌ１再構成について＞
最後に、有限要素法による関係式から最適化問題に変換するために適用したｌ１再構成について説明する。
ｌ１再構成は、従来より圧縮センシングの手法として適用されている。画像、音、声、動画などのデータは概してサイズが大きいため、保存の際に圧縮されることが多い。このようなデータ圧縮では、通常、大規模データがもつ冗長さを利用するが、この冗長さを見つけるためには、一旦、原データを保存する必要があり、大規模な記憶装置および多大な計測時間が必要となる。また、取得したデータの冗長な部分は、圧縮に際して失われるだけであるため、データを予め圧縮したかたちで取得し、それを原データと同質に復元することができれば、記憶容量やデータ取得の時間の削減に繋がり、より効率的である。
【００６８】
圧縮センシングの基本的な問題設定を考えるにあたって、以下の式（１１）のような多元連立一次方程式を考える。
【００６９】
【数１１】
0006648883.tif 000012
【００７０】
ただし、ｙ１・・・ｙＭは観測情報、ｘ１・・・ｘＮは推定対象となる原情報である。式（１１）は、Ｎ＞Ｍであれば何らかの拘束なしには解を一意に得ることはできない。しかし、原情報ベクトルｘが多くのゼロ成分をもつ「スパース」なベクトルであると仮定すると、実質的に非ゼロ成分数だけ観測値ｙがあればｘを推定することができる。例えば、式（１１）において、Ｎ＝７，Ｍ＝３とした場合、ｘ１〜ｘ７のうち、ｘ２＝ｘ３＝ｘ４＝ｘ６＝ｘ７＝０であったとすると、式（１１）は以下の式（１２）のように書き変えることができる。
【００７１】
【数１２】
0006648883.tif 000013
【００７２】
式（１２）では，ｘの非ゼロ成分数２だけ観測値（例えばｙ１，ｙ２）があればｘを推定することができる。つまり、元の測定情報ｙは３成分あったが、ｘの非ゼロ成分の数、すなわち２成分にまで圧縮して観測すればよいことになる。
Ｍ×Ｎセンシング行列Ａと、Ｍ次元測定値ベクトルｙが与えられたとき、線形方程式  ｙ＝Ａｘ  を満たすＮ次元ベクトルｘを求めるような問題において、ｘのスパース性を仮定できるとき、復元に際して測定値数Ｍ（＜Ｎ）をどの程度まで小さくできるかという点が圧縮センシングの基本的な問題設定となる。
【００７３】
ｘを求めるにあたり、まず以下のようなｌ０再構成を考える。ベクトルｘの０−ノルムをｘの非ゼロ要素の個数と定義すると、式（１３）のように、線形制約  ｙ＝Ａｘ  のもとで０−ノルムを最小化することによってｘの推定値を得るという推定法が考えられ、これがｌ０再構成と呼ばれる。
【００７４】
【数１３】
0006648883.tif 000014
【００７５】
しかし、ｌ０再構成は各要素を非ゼロか否か調べることになるため，一般に組合せ最適化問題となり、ＮＰ困難である。したがって、計算時間がＮに関して指数関数的に増大するという欠点がある。
【００７６】
そこで、次のように０−ノルムから１−ノルムへ問題を緩和して考え、ベクトルｘの１−ノルムを式（１４）のように定義すると、式（１５）のような最適化問題とみなすことができる。これをｌ１再構成と呼ぶ。
【００７７】
【数１４】
0006648883.tif 000015
【００７８】
【数１５】
0006648883.tif 000016
【００７９】
ｌ１再構成は線形計画問題として定式化でき、シンプレックス法又は内点法等の公知の解法によって効率的に解くことができる。それらの解法を用いるため、１−ノルムの絶対値を外すには、ｘのｉ番目の要素ｘｉは、非負整数ｘ＋ｉ、ｘ−ｉを用いて式（１６）のように表現することができる。
【００８０】
【数１６】
0006648883.tif 000017
【００８１】
ここで、ｘｉ≧０ならばｘ−ｉ＝０、ｘｉ≦０ならばｘ＋ｉ＝０と考えると、式（１７）が成り立ち、式（１６）及び式（１７）は、それぞれ式（１８）、式（１９）のように表現できる。
【００８２】
【数１７】
0006648883.tif 000018
【００８３】
【数１８】
0006648883.tif 000019
【００８４】
【数１９】
0006648883.tif 000020
【００８５】
そのため、式（１５）は、次の式（２０）となる。
また、ｘ＝ｘ＋−ｘ−（式（１７））より、次の式（２１）のような線形計画問題に帰着する。
【００８６】
【数２０】
0006648883.tif 000021
【００８７】
【数２１】
0006648883.tif 000022
【００８８】
本発明では、このようなｌ１再構成の考え方を有限要素解析に応用することによって、不良設定問題を解くことを可能としている。
【符号の説明】
【００８９】
１０      ：外力推定システム
１２      ：外力推定装置
２１      ：取得部
２２      ：推定部
２３      ：修正部
Ａ１      ：視認領域
Ａ２      ：非視認領域
Ｋ        ：剛性マトリクス（剛性情報）
Ｗ        ：対象物（弾性体）
ｆ        ：外力ベクトル
ｕ      ：変位ベクトル
【図１】
0006648883.tif 000023
【図２】
0006648883.tif 000024
【図３】
0006648883.tif 000025
【図４】
0006648883.tif 000026
【図５】
0006648883.tif 000027
【図６】
0006648883.tif 000028
【図７】
0006648883.tif 000029
【図８】
0006648883.tif 000030
【図９】
0006648883.tif 000031
【図１０】
0006648883.tif 000032