Specification :(In Japanese)外力推定装置及び方法、並びにコンピュータプログラム

Country	(In Japanese)日本国特許庁(JP)
Gazette	(In Japanese)特許公報(B2)
Patent Number	P6648883
Publication number	P2017-000623A
Date of registration	Jan 20, 2020
Date of issue	Feb 14, 2020
Date of publication of application	Jan 5, 2017
Title of the invention, or title of the device	(In Japanese)外力推定装置及び方法、並びにコンピュータプログラム
IPC (International Patent Classification)	A61B  34/10        (2016.01) G06F  30/23        (2020.01)
FI (File Index)	A61B 34/10 G06F 17/50 612H
Number of claims or invention	5
Total pages	17
Application Number	P2015-120581
Date of filing	Jun 15, 2015
Exceptions to lack of novelty of invention	(In Japanese)特許法第30条第2項適用 生体医用画像研究会 第2回若手発表会 平成27年3月14日開催 〔刊行物等〕 平成27年3月14日掲載 p：／／biomedimg．umin．jp／proceedings－2ndConfMBI2014．pdf 〔刊行物等〕 第54回日本生体医工学大会 平成27年5月7－9日開催（5月8日発表） 〔刊行物等〕 平成27年5月（平成27年5月7日配布）生体医工学 第53巻 第54回日本生体医工学大会 プログラム・抄録集 第54回日本生体医工学会大会事務局 名古屋工業大学大学院 未来医療介護健康情報学研究所
Date of request for substantive examination	Jun 1, 2018
Patentee, or owner of utility model right	(In Japanese)【識別番号】504132272 【氏名又は名称】国立大学法人京都大学
Inventor, or creator of device	(In Japanese)【氏名】中尾 恵 【氏名】松田 哲也 【氏名】坂田 良平
Representative	(In Japanese)【識別番号】110000280、【氏名又は名称】特許業務法人サンクレスト国際特許事務所
Examiner	(In Japanese)【審査官】宮下 浩次
Document or reference	(In Japanese)田中 利幸，圧縮センシングの数理，電子情報通信学会 基礎・境界ソサイエティ，日本，一般社団法人 電子情報通信学会，2010年，Vol.4 No.1，p39-47
Field of search	Ａ６１Ｂ ３４／００ － ３４／３７ Ｇ０６Ｆ １７／５０
Scope of claims	(In Japanese)【請求項1】 スパース性を有する外力が付与される弾性体における視認可能な第1領域と視認不能な第2領域のうち、前記第1領域の第1変位の情報を取得する取得部と、 前記弾性体を表す所定のメッシュモデルにおいて、前記第1領域の第1変位と、前記第1領域の第1外力及び前記第2領域の第2外力との間に成立する有限要素法の関係式に、l1再構成を適用して最適化問題に変換することにより、前記第1外力及び前記第2外力の双方を推定する推定部と、を備えている、外力推定装置。 【請求項2】 前記推定部は、 推定された前記第2外力を用いて、前記第2領域の第2変位を推定する、請求項1に記載の外力推定装置。 【請求項3】 前記弾性体に付与された既知の外力と、前記推定部が推定した前記第1外力との差に基づいて、前記メッシュモデルの剛性情報を修正する修正部を更に備えている、請求項1又は2に記載の外力推定装置。 【請求項4】 コンピュータが、スパース性を有する外力が付与される弾性体における視認可能な第1領域と視認不能な第2領域のうち、前記第1領域の第1変位の情報を取得するステップと、 前記コンピュータが、前記弾性体を表す所定のメッシュモデルにおいて、前記第1領域の第1変位と、前記第1領域の第1外力及び前記第2領域の第2外力との間に成立する有限要素法の関係式に、l1再構成を適用して最適化問題に変換することにより、前記第1外力及び前記第2外力の双方を推定するステップとを含む、外力推定方法。 【請求項5】 スパース性を有する外力が付与される弾性体における視認可能な第1領域と視認不能な第2領域のうち、前記第1領域の第1変位の情報を取得する機能、及び、 前記弾性体を表す所定のメッシュモデルにおいて、前記第1領域の第1変位と、前記第1領域の第1外力及び前記第2領域の第2外力との間に成立する有限要素法の関係式に、l1再構成を適用して最適化問題に変換することにより、前記第1外力及び前記第2外力の双方を推定する機能を、コンピュータに実現させるコンピュータプログラム。
Detailed description of the invention	(In Japanese)【技術分野】 【0001】 本発明は、外力推定装置及び方法、並びにコンピュータプログラムに関する。 【背景技術】 【0002】 近年、患者への侵襲を少なくし術後のQOL向上等を図ることを目的として、体に形成した小さな穴から内視鏡カメラや術具等を挿入し、内視鏡カメラで術野を確認しながら術具により手術を行う内視鏡下手術が数多く実施されている。また、最近では、遠隔からの操縦によってロボットに術具を操作させるロボット手術も実施されてきている。 【0003】 しかし、内視鏡下手術やロボット手術においては、術野が狭いために手術対象となる臓器の全体像を視認することができない。また、ロボット手術においては、鉗子を用いて臓器に力を加えたとしてもその反力が術者へは提示されないので、術者は、狭い術野の中で臓器の変形を確認しながら臓器にかかる力を想像しなければならない。そのため、鉗子から臓器に過剰な力が加わって臓器や血管を傷つけてしまうおそれがあり、術具の操作には熟練した技能が求められる。 このような実情に鑑み、鉗子に加わる力を計測し、手術の技能分析に役立てる研究が行われている（例えば、非特許文献1参照）。 【先行技術文献】 【0004】 【非特許文献1】吉田健志、外6名、「腹腔鏡手術用鉗子先端に加わる作用力の計測と剥離操作における技能分析」、生体医工学、Vol．48（2010）、No．1、P25－32 【発明の概要】 【発明が解決しようとする課題】 【0005】 ロボット手術による術中に、非特許文献1に記載の技術を用いて鉗子に加わる力を計測することができれば、臓器に加わる力を客観的に把握でき、臓器等の損傷を未然に防ぐことができると考えられるが、術野の範囲外で鉗子を操作することは困難であるため、術野の範囲外において臓器にどのような力が加わっているのかを把握することはできない。また、術野の範囲内で鉗子により臓器を引っ張った場合、術野の範囲外で他の組織に繋がった部分にも力が加わると考えられるが、そのような力も当然に把握することができない。 【0006】 そこで、本発明は、弾性体の局所的な変位情報から、弾性体の全体に加わる外力を推定することができる外力推定装置及び方法、並びにコンピュータプログラムを提供することを目的とする。 【課題を解決するための手段】 【0007】 （1） 本発明に係る外力推定装置は、スパース性を有する外力が付与される弾性体における視認可能な第1領域と視認不能な第2領域のうち、前記第1領域の第1変位の情報を取得する取得部と、 前記弾性体を表す所定のメッシュモデルにおいて、前記第1領域の第1変位と、前記第1領域の第1外力及び前記第2領域の第2外力との間に成立する有限要素法の関係式に、l1再構成を適用して最適化問題に変換することにより、前記第1外力及び前記第2外力の双方を推定する推定部と、を備えたものである。 【0008】 上記の本発明は、物体の構造解析等に多用されている有限要素法を用いて、弾性体の変位と弾性体の剛性情報とから弾性体に付与される外力を推定する。本来、有限要素法において、弾性体の全ての要素に付与される外力を取得するには、全ての要素における変位の情報が必要である。しかし、内視鏡下手術やロボット手術等では、術野の範囲外から変位の情報を取得することは困難である。 【0009】 本出願の発明者は、手術対象となる臓器等の弾性体に鉗子等の小さな術具によって力を付与すると、術具が直接的に接する部分や、他の組織に繋がる部分において外力が付与され、外力が付与される領域の体積は、弾性体全体の体積よりも十分に小さく、弾性体に付与される外力にはスパース性があることに気づいた。そのため、弾性体の多くの要素に付与される外力をゼロと仮定でき、少ない変位の情報からでも非ゼロの外力を求めることが可能であると考えた。そして、このような考えを実現するため、少ない変位の情報を用いて有限要素法による関係式を構築するとともに、当該関係式にl1－再構成を適用し最適化問題に変換するという本発明を完成するに到った。このような発明によって、弾性体の一部の領域の変位の情報から、簡単な計算で弾性体全体の外力を推定することが可能となる。 【0010】 （2） 前記推定部は、推定された前記第2外力を用いて、前記第2領域の第2変位を推定することが好ましい。 上記のように推定部は、弾性体の一部の領域である第1領域の変位を用いて、弾性体全体に付与される第1外力及び第2外力を推定するが、その推定された第2外力をさらに用いることによって、前記関係式から弾性体の第2領域の第2変位をも推定することができる。したがって、例えば内視鏡カメラによって視認できない領域も含めて弾性体全体の変位を把握することができる。 【0011】 （3） 外力推定装置は、前記弾性体に付与された既知の外力と、前記推定部が推定した前記第1外力との差に基づいて、前記メッシュモデルの剛性情報を修正する修正部を更に備えていることが好ましい。 剛性情報は、例えば弾性体のヤング率やポアソン比等から算出することが可能であるが、弾性体が手術対象となる臓器等の場合には、その臓器から直接的にこれらのパラメータを取得するのは困難であるため、その臓器の標準的な剛性情報を採用することが考えられる。しかしながら、標準的な剛性情報と実際の剛性情報との相違によって、推定部によって推定された外力と実際に付与された外力とにずれが生じるおそれもある。本発明の外力推定装置は、推定部によって推定された外力と、既知の外力とを比較し、両者の差に基づいて剛性情報を修正することで、より正確に外力を推定することが可能となる。 【0012】 （4） 本発明の外力推定方法は、コンピュータが、スパース性を有する外力が付与される弾性体における視認可能な第1領域と視認不能な第2領域のうち、前記第1領域の第1変位の情報を取得するステップと、 前記コンピュータが、前記弾性体を表す所定のメッシュモデルにおいて、前記第1領域の第1変位と、前記第1領域の第1外力及び前記第2領域の第2外力との間に成立する有限要素法の関係式に、l1再構成を適用して最適化問題に変換することにより、前記第1外力及び前記第2外力の双方を推定するステップとを含むものである。 【0013】 （5） 本発明のコンピュータプログラムは、スパース性を有する外力が付与される弾性体における視認可能な第1領域と視認不能な第2領域のうち、前記第1領域の第1変位の情報を取得する機能、及び、前記弾性体を表す所定のメッシュモデルにおいて、前記第1領域の第1変位と、前記第1領域の第1外力及び前記第2領域の第2外力との間に成立する有限要素法の関係式に、l1再構成を適用して最適化問題に変換することにより、前記第1外力及び前記第2外力の双方を推定する機能を、コンピュータに実現させるものである。 【発明の効果】 【0014】 本発明によれば、弾性体の一部の領域における変位情報から、弾性体の全体に付与される外力を推定することができる。 【図面の簡単な説明】 【0015】 【図1】本発明の実施形態に係る外力推定装置の構成図である。 【図2】弾性体を説明する図である。 【図3】表示装置に表示する推定結果の一例を示す図である。 【図4】有限要素法で用いる対象物のメッシュモデルの一例を示す図である。 【図5】評価実験に用いる肝臓のモデルを示す図である。 【図6】評価実験に用いる肝臓のモデルに推定対象となる外力を付与した状態を示す図である。 【図7】評価実験において、1領域に外力を加えた場合の推定結果を示す図である。 【図8】評価実験において、2領域に外力を加えた場合の推定結果を示す図である。 【図9】評価実験において、1領域に外力を加えた場合の推定誤差を示すグラフである。 【図10】評価実験において、2領域に外力を加えた場合の推定誤差を示すグラフである。 【発明を実施するための形態】 【0016】 ＜1．実施形態＞ 以下、図面を参照して本発明の実施形態を詳細に説明する。 ［外力推定システムの全体構成］ 図1は、本発明の実施形態に係る外力推定システムの説明図である。 本実施形態の外力推定システム10は、弾性体からなる対象物の画像を撮影するカメラ（撮像装置）11と、カメラ11により撮影された画像の信号が入力され、所定の処理を行うことによって対象物に付与される外力を推定する外力推定装置12と、外力推定装置12の出力結果等を表示する表示装置13とを備えている。 【0017】 この外力推定システム10は、例えば、内視鏡下手術やロボット手術の支援のために用いることができる。 カメラ11は、CCDやCMOS等の撮像素子を備えたものであり、内視鏡に取り付けられ、患者の体腔内に挿入されることによって、対象物としての臓器等の術部を撮影する。カメラ11は、1つであってもよいし、複数であってもよい。 【0018】 外力推定装置12は、カメラ11によって撮影された術部の画像から対象物に付与される外力を推定するものである。外力の情報として、外力の位置、大きさ、向きが推定される。また、外力推定装置12は、カメラ11によって撮影された部分だけでなく撮影されていない部分に付与される外力をも推定する。例えば、図2に示すように、対象物Wの一部の領域A1は、カメラ11の画像から視認できる状態であるが、他の領域A2は、カメラ11によって撮影されず、視認できない状態にあるものとする。この場合、外力推定装置12は、視認領域A1に付与される外力F1だけでなく、非視認領域A2に付与される外力F2も推定する。推定された外力の情報は、表示装置13に表示される。 【0019】 なお、手術中、臓器には、鉗子等の小さな術具によって部分的に外力が付与されることが多く、また、術具が直接的に作用する部分以外は、他の組織と繋がっている部分だけに外力が付与されると考えられる。したがって、臓器に対して外力が付与される領域の体積は、臓器全体の体積よりも十分に小さいと考えられ、臓器に付与される外力には「スパース性」があるといえる。外力推定装置12は、このような外力のスパース性に着目して、臓器の一部の領域の画像から外力を推定するものとなっている。外力推定装置12の具体的な構成と、外力推定のための処理アルゴリズムについては後述する。 【0020】 表示装置13は、術者又は手術スタッフ（以下、「術者等」ともいう）の視界に配置されたモニター、又は、術者等が装着したヘッドマウントディスプレイ等とすることができる。そのため、術者は、表示装置13に表示された外力の情報をもとに臓器に付与される外力を客観的に把握することができ、臓器に対して過剰な力を付与しないように配慮して手術を行うことが可能になる。 【0021】 図3は、表示装置13への表示例を示す。対象物を模した図に対して外力が付与されている部分Rと外力の大きさを色彩の変化又は濃淡の変化等によって示したものである。対象物の図には、カメラを介して視認できる領域A1と視認できない領域A2とを識別できるように、2つの領域A1，A2が区画線Cによって区画されている。このような図を表示装置13に表示させることによって、術者等は外力が生じている部位と外力の大きさとを一目で把握することができる。 【0022】 ［外力推定装置の具体的構成］ 外力推定装置12は、例えばパーソナルコンピュータから構成され、CPU等の演算部と、ROM、RAM、ハードディスク等の記憶部と、各種入出力インタフェース等を有している。カメラ11によって撮影された画像の情報は、入出力インタフェースを介して外力推定装置12に入力される。そして、外力推定装置12は、記憶部にインストールされたプログラムを演算部が実行することによって、入力された情報を用いて外力を推定し、推定結果を表示装置13に表示させる。記憶部には、外力の推定に用いられる種々の情報も記憶されている。 【0023】 外力推定装置12は、その機能構成として、カメラ11から入力された画像情報を処理し、撮影された領域の変位についての情報を取得する取得部21と、取得部21において取得された変位についての情報を用いて対象物に付与される外力を推定する推定部22と、推定部22の処理に用いられる所定の定数を修正する修正部23とを有している。 【0024】 取得部21は、時系列で取得された画像において特徴点の位置を変化を追うことによって変位情報を取得する。この取得部21における具体的な処理は、従来公知の技術を用いることができる。例えば、手術中に取得したカメラ画像によって弾性体表面の形状変化を追う技術（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.）等を用いることができる。 また、取得部21は、例えば、術前に撮影されたCT（Computed Tomography）画像を用いて3次元の初期モデルを生成し、カメラによって撮影された臓器の画像に基づいて初期モデルの形状を順次更新していき、その更新の度に変位を取得する方法を用いてもよい。 【0025】 推定部22は、取得部21によって取得された変位の情報に基づいて、弾性体に付与される外力を推定する。以下、推定部22における具体的な処理アルゴリズムについて説明する。 ［推定部22における処理アルゴリズム］ 推定部22は、有限要素法を用いることによって外力を推定する。有限要素法では、弾性体を多数の要素に分割し、各要素の頂点における変位と外力との関係として、以下の式（1）を用いる。 【0026】 【数1】 【0027】 ここで、fは、外力ベクトル、Kは、剛性マトリクス（剛性情報）、uは、変位ベクトルである。剛性マトリクスKは、弾性体のヤング率、ポアソン比等から求められる定数である。 また、式（1）は、次の式（2）のように変形することができる。 【0028】 【数2】 【0029】 図4に、弾性体を複数の要素に分割したメッシュモデルの一例を示す。この弾性体は立方体形状であり、図4（a）は、外力が付与される前の状態、図4（b）は、外力が付与された状態を示す。メッシュモデルの頂点t1（黒い点）は固定点、頂点t2（濃いグレーの点）は自由点、頂点t3（グレーの点）は観測点、頂点t4（薄いグレーの点）は外力が加わる点を示す。頂点t4には、矢印S1で示すように推定対象となる外力が付与されている。観測点t3は、図2に示すように、外部から視認できる領域A1に配置される頂点に相当する。 【0030】 本実施形態では、既知条件として、図4のメッシュモデルの観測点t3における変位ベクトルuoと、剛性マトリクスKとを与え、メッシュモデルの全頂点に加わるスパース性のある外力ベクトルfを推定する。変位ベクトルuoは、外力推定装置12の取得部21によって取得される。剛性マトリクスKは、定数として予め与えられる。 【0031】 メッシュモデルの頂点を、視認できる頂点oと視認できない頂点iに分類し、マトリクスLをブロック分割すると、式（2）は、次の式（3）のように表すことができる。 【0032】 【数3】 【0033】 ただし、uo及びfoは、視認できる頂点o（観測点t3）における変位ベクトル及び外力ベクトルであり、ui及びfiは、視認できない頂点iにおける変位ベクトル及び外力ベクトルである。 そして、視認できる頂点oにおける変位ベクトルuoが既知量であるため、次の式（4）のように、当該変位ベクトルuoのみから外力ベクトルfo、fiを求めることを考える。 【0034】 【数4】 【0035】 式（4）は、既知量uoの成分数よりも、未知量fo，fiの成分数の方が多い多元連立一次方程式である。そのため、何らかの拘束条件なしに一意に解くことはできない不良設定問題となる。本実施形態では、式（4）に対してl1再構成を適用することによって、外力ベクトルfの推定を可能にする。なお、l1再構成の一般的な考え方については、後の＜3．l1再構成について＞の項目において述べることとする。 【0036】 まず、外力ベクトルfの1－ノルムを式（5）のように定義すると、上記の式（4）は、次の式（6）のような最適化問題を解くことになる。 【0037】 【数5】 【0038】 【数6】 【0039】 そして、上記最適化問題は、以下のような線形計画問題として記述することができる。 具体的には、式（6）は、fのi番目の要素fiについて、非負整数f＋i，f－iを用いると、式（7）のように記述することができる。 【0040】 【数7】 【0041】 そして、f＝f＋－f－より、次の式（8）を得ることができる。 【0042】 【数8】 【0043】 なお、式（7）及び式（8）を導出する詳細な手順は、後述の＜3．l1再構成について＞の項目において説明する。 【0044】 推定部22は、上記のような線形計画問題を公知の解法を用いて解き、外力ベクトルfを推定する。公知の解法として、例えばシンプレックス法又は内点法を採用することができる。推定された外力は、図3に示すように図式化され、表示装置13に表示される。 【0045】 また、推定部22は、推定された外力ベクトルfを用いることによって、視認できない頂点における変位ベクトルuiを上記の式（2）により求める。これにより、既知量である変位ベクトルuoに加え、未知量である変位ベクトルuiをも求めることができ、弾性体全体の外力だけでなく、変位の状態をも把握することができる。 【0046】 外力推定装置12の修正部23は、定数である剛性マトリクスKを適正なかたちに修正するものである。体内の臓器を手術するような場合、その臓器自体の剛性マトリクスKを直接的に求めることは困難であるため、外力推定装置12の推定部22は、臓器の標準的な剛性マトリクスを用いる。しかしながら、標準的な剛性マトリクスと、実際の臓器の剛性マトリクスとは厳密に一致しないので、推定される外力ベクトルfや変位ベクトルuiが実際の外力や変位から若干ずれる可能性がある。 【0047】 このような問題を解決するため、修正部23は、例えば、手術前に臓器に既知の外力を付与し、その外力による変位uoをカメラ11の画像から取得し、この変位uoを用いて上記の式（1）～（8）を用いて外力ベクトルfを推定する。そして、次の式（9）を用いることによって、既知の外力f’と、推定された外力ベクトルf（＝Ku）と差分を最小化するようなKの値を求める。これにより、実際の対象物により適した剛性マトリクスKを求めることができ、より正確な外力及び未知の変位の推定が可能となる。 【0048】 【数9】 【0049】 以上、説明したように、本実施形態の外力推定装置においては、弾性体にスパースな外力が付与される場合に、有限要素法の関係式にl1再構成を適用し、最適化問題（線形計画問題）に変換することで、弾性体の一部の領域における変位を用いて、弾性体全体の外力を推定する。 したがって、内視鏡下手術やロボット手術において、狭い術野内で臓器の一部しか視認することができず、その一部の変位しか取得できない場合であっても、視認できる臓器の一部に加え、視認できない臓器の他の部分についても外力を推定することができる。そのため、鉗子等の術具によって臓器に力を付与する場合に、視認できる部分だけでなく視認できない部分においても損傷が生じないように配慮して手術を行うことができる。 【0050】 ［変形例］ 本発明の外力推定装置は、上述したように、内視鏡下手術やロボット手術において手術対象となる臓器等の外力を推定するために好適に利用することができるが、これに限定されるものではなく、あらゆる物体に付与される外力を推定するために好適に利用することができる。 例えば、本発明は、木造やコンクリート造等の構造物や建築物に対して付与されるスパースな外力を推定するために用いることができる。この場合、視認できない部分の外力を推定することで、構造物等の破損、倒壊等を予測することが可能となる。また、地質等、表面しか観測できないような対象物において、内部に付与される外力を推定することができる。また、本発明は、生体の細胞等のような微小な物体にも適用することができる。 また、対象物の画像を撮影する撮像装置は、上述したようなカメラに限らず、対象物の一部を画像化できるものであれば特に限定されない。例えば、手術中であっても撮影が可能な、超音波検査装置、X線撮影装置、CT撮影装置、MRI（Magnetic Resonance Imaging）撮影装置等を採用することもできる。また、対象物が構造物や建築物の場合には、レンジセンサ（デプスセンサ）を採用することもできる。 もっとも、特許請求の範囲に記載の「弾性体」には、上記の木造やコンクリート造等の構造物や建築物は含まれない。 【0051】 ＜2．評価実験＞ 本出願の発明者は、上記において説明した処理アルゴリズムの有効性を確認するため、評価実験を行った。この評価実験では、対象物として、図5に示すような肝臓のメッシュモデルを用いた。図5において、領域r1内の頂点t1は固定点であり、その他の頂点t2は自由点である。固定点t1は、肝臓に下大動脈が接続されていることを想定して設定している。このメッシュモデルは、固定点t1を除いて表面に178個、内部に85個の頂点があり、合計263個の頂点を有する。図5には、3次元の座標系も併せて示す。 【0052】 第1の評価実験として、図6（a）に示すように、メッシュモデルの視認できる1つの領域に、鉗子で引っ張り力を与えた場合を想定して－z方向に外力を加えた。具体的には、領域r4内の9個の頂点t4に外力を加えた。その外力を矢印S1で示す。 また、第2の評価実験として、図6（b）に示すように、メッシュモデルの視認できる1つの領域に、鉗子で引っ張り力を与えた場合を想定して－z方向の外力を加え、それによって、視認できない他の領域にも外力が加わったと想定して、＋x方向に外力を加えた。これらの外力を矢印S1で示す。－z方向の外力は領域r4内の9個の頂点t4に加え、＋x方向の外力は領域r5内の7個の頂点t5に加えた。 【0053】 そして、第1，第2の評価実験の双方とも、視認できる領域において頂点の変位を取得し、その変位から対象物全体の外力を推定した。そして、推定された外力と、実際に与えた外力とを比較し、推定精度を評価した。 【0054】 図7は、第1の評価実験における推定結果の1例を示す図である。図7には、実際に加えた外力を黒い矢印S1で示し、推定された外力をグレーの矢印S2で示している。この推定結果は、固定点t1を除くメッシュモデルの263頂点のうち、z座標の小さい頂点から順に53点（領域r3内の頂点t3）を観測した例を示す。 【0055】 図7（b）は、図7（a）のA部を拡大したものである。図7（b）から明らかなように、推定された外力S2と、実際に加えた外力S1とは、概ね一致していることがわかる。 また、9個の頂点に実際に加えた外力の大きさは5Nであったのに対して、その頂点において推定された外力の大きさは、3．7～6．3Nであった。それ以外の頂点の外力は、これよりも10～102のオーダーで小さい値であった。よって、実際に外力を加えた頂点について、実際に近い大きさの外力を高精度に推定することができた。 【0056】 図8は、第2の評価実験における推定結果の1例を示す図である。図8には、実際に加えた外力を黒い矢印S1で示し、推定された外力をグレーの矢印S2で示している。この推定結果は、固定点t1を除くメッシュモデルの263頂点のうち、z座標の小さい頂点から順に79点（領域r3内の頂点t3）を観測した例を示す。 【0057】 図8（b）（c）は、それぞれ図8（a）のB部及びC部を拡大したものである。図8（b）（c）から明らかなように、推定された外力S2と、実際に加えた外力S1とは、概ね一致していることがわかる。 また、16個の頂点に実際に加えた外力の大きさは5Nであったのに対して、その頂点において推定された外力の大きさは、3．9～5．3Nであった。それ以外の頂点で推定された外力は、これよりも10～102のオーダーで小さい値であった。よって、実際に外力を加えた頂点について、実際に近い大きさの外力を高精度に推定することができた。 【0058】 図9（a）は、第1の評価実験でメッシュモデルの表面のみを観測した場合について、表面の頂点178点のうち、観測点をz座標の小さいものから順に9点～178点まで増加させた場合の推定誤差を示している。観測点を9点以上としたのは、外力を加えた点数以上とするためである。 同様に、図9（b）は、第1の評価実験でメッシュモデルの内部も含めて観測した場合について、表面及び内部の頂点263点のうち、観測点をz座標の小さいものから順に9点～234点まで増加させた場合の推定誤差を示している。 【0059】 ここで用いる推定誤差は、与えた外力と推定外力との差の二乗平均平方根（Root Mean Square；以下、「RMS」という）とし、次の式（10）により求めた。 【0060】 【数10】 【0061】 fxi，fyi，fziは、実際に加えた外力fの第i要素fiについてのxyz成分である。f’xi，f’yi，f’ziは、推定された外力f’の第i要素f’iについてのxyz成分である。 Nは、観測する頂点をメッシュモデルの表面のみとした場合、メッシュモデルの全頂点のうち、固定点と内部の点を除いた頂点の数であり、観測する頂点にメッシュモデルの内部を含めた場合、メッシュモデルの全頂点のうち、固定点を除いた頂点の数である。 【0062】 RMSは、加えた外力と推定された外力の要素の相違が生じている場合や、加えた外力fと推定された外力f’のxyz成分ごとの大きさのずれが生じている場合に、増加する。言い換えると、正しく推定されるほど、RMSの値は小さくなる。 【0063】 図9（a）（b）に示すように、観測する頂点数に関わらず、ほぼRMSは1以下の値となった。また、観測する頂点数が多くなるほど推定誤差RMSが小さくなることが分かる。 【0064】 図10（a）は、第2の評価実験でメッシュモデルの表面のみを観測した場合について、表面の頂点178点のうち、観測点をz座標の小さいものから順に16点～178点まで増加させた場合の推定誤差を示している。観測点を16点以上としたのは、外力を加えた点数以上とするためである。 【0065】 同様に、図10（b）は、第2の評価実験でメッシュモデルの内部も含めて観測した場合について、表面及び内部の頂点263点のうち、観測点をz座標の小さいものから順に16点～234点まで増加させた場合の推定誤差を示している。 【0066】 図10（a）（b）に示すように、観測する頂点数に関わらず、ほぼRMSは1以下の値となった。また、図10（a）に示すように、観測する頂点数が多くなるほど、概ね推定誤差RMSが小さくなることが分かる。特に、観測点が58点（33％）以上の場合、RMSは安定して小さい値をとっていることが分かる。 【0067】 ＜3．l1再構成について＞ 最後に、有限要素法による関係式から最適化問題に変換するために適用したl1再構成について説明する。 l1再構成は、従来より圧縮センシングの手法として適用されている。画像、音、声、動画などのデータは概してサイズが大きいため、保存の際に圧縮されることが多い。このようなデータ圧縮では、通常、大規模データがもつ冗長さを利用するが、この冗長さを見つけるためには、一旦、原データを保存する必要があり、大規模な記憶装置および多大な計測時間が必要となる。また、取得したデータの冗長な部分は、圧縮に際して失われるだけであるため、データを予め圧縮したかたちで取得し、それを原データと同質に復元することができれば、記憶容量やデータ取得の時間の削減に繋がり、より効率的である。 【0068】 圧縮センシングの基本的な問題設定を考えるにあたって、以下の式（11）のような多元連立一次方程式を考える。 【0069】 【数11】 【0070】 ただし、y1・・・yMは観測情報、x1・・・xNは推定対象となる原情報である。式（11）は、N＞Mであれば何らかの拘束なしには解を一意に得ることはできない。しかし、原情報ベクトルxが多くのゼロ成分をもつ「スパース」なベクトルであると仮定すると、実質的に非ゼロ成分数だけ観測値yがあればxを推定することができる。例えば、式（11）において、N＝7，M＝3とした場合、x1～x7のうち、x2＝x3＝x4＝x6＝x7＝0であったとすると、式（11）は以下の式（12）のように書き変えることができる。 【0071】 【数12】 【0072】 式（12）では，xの非ゼロ成分数2だけ観測値（例えばy1，y2）があればxを推定することができる。つまり、元の測定情報yは3成分あったが、xの非ゼロ成分の数、すなわち2成分にまで圧縮して観測すればよいことになる。 M×Nセンシング行列Aと、M次元測定値ベクトルyが与えられたとき、線形方程式 y＝Ax を満たすN次元ベクトルxを求めるような問題において、xのスパース性を仮定できるとき、復元に際して測定値数M（＜N）をどの程度まで小さくできるかという点が圧縮センシングの基本的な問題設定となる。 【0073】 xを求めるにあたり、まず以下のようなl0再構成を考える。ベクトルxの0－ノルムをxの非ゼロ要素の個数と定義すると、式（13）のように、線形制約 y＝Ax のもとで0－ノルムを最小化することによってxの推定値を得るという推定法が考えられ、これがl0再構成と呼ばれる。 【0074】 【数13】 【0075】 しかし、l0再構成は各要素を非ゼロか否か調べることになるため，一般に組合せ最適化問題となり、NP困難である。したがって、計算時間がNに関して指数関数的に増大するという欠点がある。 【0076】 そこで、次のように0－ノルムから1－ノルムへ問題を緩和して考え、ベクトルxの1－ノルムを式（14）のように定義すると、式（15）のような最適化問題とみなすことができる。これをl1再構成と呼ぶ。 【0077】 【数14】 【0078】 【数15】 【0079】 l1再構成は線形計画問題として定式化でき、シンプレックス法又は内点法等の公知の解法によって効率的に解くことができる。それらの解法を用いるため、1－ノルムの絶対値を外すには、xのi番目の要素xiは、非負整数x＋i、x－iを用いて式（16）のように表現することができる。 【0080】 【数16】 【0081】 ここで、xi≧0ならばx－i＝0、xi≦0ならばx＋i＝0と考えると、式（17）が成り立ち、式（16）及び式（17）は、それぞれ式（18）、式（19）のように表現できる。 【0082】 【数17】 【0083】 【数18】 【0084】 【数19】 【0085】 そのため、式（15）は、次の式（20）となる。 また、x＝x＋－x－（式（17））より、次の式（21）のような線形計画問題に帰着する。 【0086】 【数20】 【0087】 【数21】 【0088】 本発明では、このようなl1再構成の考え方を有限要素解析に応用することによって、不良設定問題を解くことを可能としている。 【符号の説明】 【0089】 10 ：外力推定システム 12 ：外力推定装置 21 ：取得部 22 ：推定部 23 ：修正部 A1 ：視認領域 A2 ：非視認領域 K ：剛性マトリクス（剛性情報） W ：対象物（弾性体） f ：外力ベクトル u ：変位ベクトル
Drawing	(In Japanese)【図1】 (In Japanese)【図2】 (In Japanese)【図3】 (In Japanese)【図4】 (In Japanese)【図5】 (In Japanese)【図6】 (In Japanese)【図7】 (In Japanese)【図8】 (In Japanese)【図9】 (In Japanese)【図10】