ホーム > 特許ランキング > 一般財団法人電力中央研究所
知財求人 - 知財ポータルサイト「IP Force」
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-08-23
(45)【発行日】2023-08-31
(54)【発明の名称】粒子のファブリックの解析方法、解析装置、及び解析プログラム
(51)【国際特許分類】
G06T 7/62 20170101AFI20230824BHJP
G01N 23/046 20180101ALI20230824BHJP
G06T 7/60 20170101ALI20230824BHJP
【ＦＩ】
G06T7/62
G01N23/046
G06T7/60 150S
【請求項の数】
24
(21)【出願番号】P 2019211917
(22)【出願日】2019-11-25
(65)【公開番号】P2021086160
(43)【公開日】2021-06-03
【審査請求日】2022-10-18
(73)【特許権者】
【識別番号】000173809
【氏名又は名称】一般財団法人電力中央研究所
(74)【代理人】
【識別番号】100087468
【弁理士】
【氏名又は名称】村瀬 一美
(72)【発明者】
【氏名】野原 慎太郎
【審査官】淀川 滉也
(56)【参考文献】
【文献】特開昭５４－１１５１８８（ＪＰ，Ａ）
【文献】特開２０１３－２２５５０２（ＪＰ，Ａ）
【文献】M BLOOM et al.，Imaging Systems and Algorithms for the Numerical Characterization of Three-Dimensional Shapes of Granular Particles，IEEE Transactions on Instrumentation and Measurement，Vol. 59，No. 9，2010年09月，p.2365-2375，DOI: 10.1109/TIM.2009.2034579
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｔ ７／６２
Ｇ０１Ｎ ２３／０４６
Ｇ０６Ｔ ７／６０
(57)【特許請求の範囲】
【請求項1】
粒子の集合体を非破壊で撮影することにより複数の断面画像を取得し、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出し、算出した近似楕円体モデルに基づいて粒子の大きさ、形状、配向であるファブリックを解析することを特徴とする粒子のファブリックの解析方法。
【請求項2】
算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルの粒子についてはファブリックの解析対象から除外することを特徴とする請求項１に記載の粒子のファブリックの解析方法。
【請求項3】
前記近似楕円体モデルの精度は、撮影により取得した断面画像での断面積と、前記断面画像に相当する前記近似楕円体モデルの断面積との誤差に基づいて検証することを特徴とする請求項１又は２に記載の粒子のファブリックの解析方法。
【請求項4】
前記近似楕円体モデルの精度は、Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、
Ｅｉの中央値Ｅｍｅｄ及び標準偏差Ｅｓｔｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外することを特徴とする請求項３に記載の粒子のファブリックの解析方法。
【請求項5】
前記近似楕円体モデルは、楕円体の方程式で表され、前記複数の断面画像は、少なくとも７枚であることを特徴とする請求項１から４のうちのいずれか一つに記載の粒子のファブリックの解析方法。
【請求項6】
前記近似楕円体モデルは、非線形最小二乗法を利用して算出することを特徴とする請求項１から５のうちのいずれか一つに記載の粒子のファブリックの解析方法。
【請求項7】
前記近似楕円体モデルは、中心位置が前記複数の断面画像の重心位置に基づいて設定されることを特徴とする請求項１から６のうちのいずれか一つに記載の粒子のファブリックの解析方法。
【請求項8】
前記複数の断面画像は、ＣＴ画像であることを特徴とする請求項１から６のうちのいずれか一つに記載の粒子のファブリックの解析方法。
【請求項9】
粒子の集合体を非破壊で撮影することにより複数の断面画像を取得する手段と、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出する手段と、算出した近似楕円体モデルに基づいて粒子の大きさ、形状、配向であるファブリックを解析する手段とを有することを特徴とする粒子のファブリックの解析装置。
【請求項10】
算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルの粒子についてはファブリックの解析対象から除外する手段を有することを特徴とする請求項９に記載の粒子のファブリックの解析装置。
【請求項11】
前記近似楕円体モデルの精度は、撮影により取得した断面画像での断面積と、前記断面画像に相当する前記近似楕円体モデルの断面積との誤差に基づいて検証することを特徴とする請求項９又は１０に記載の粒子のファブリックの解析装置。
【請求項12】
前記近似楕円体モデルの精度は、Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、
Ｅｉの中央値Ｅｍｅｄ及び標準偏差Ｅｓｔｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外することを特徴とする請求項１１に記載の粒子のファブリックの解析装置。
【請求項13】
前記近似楕円体モデルは、楕円体の方程式で表され、前記複数の断面画像は、少なくとも７枚であることを特徴とする請求項９から１２のうちのいずれか一つに記載の粒子のファブリックの解析装置。
【請求項14】
前記近似楕円体モデルは、非線形最小二乗法を利用して算出することを特徴とする請求項９から１３のうちのいずれか一つに記載の粒子のファブリックの解析装置。
【請求項15】
前記近似楕円体モデルは、中心位置が前記複数の断面画像の重心位置に基づいて設定されることを特徴とする請求項９から１４のうちのいずれか一つに記載の粒子のファブリックの解析装置。
【請求項16】
前記複数の断面画像は、ＣＴ画像であることを特徴とする請求項９から１５のうちのいずれか一つに記載の粒子のファブリックの解析装置。
【請求項17】
粒子の集合体を非破壊で撮影することにより複数の断面画像を取得する処理と、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出する処理と、算出した近似楕円体モデルに基づいて粒子の大きさ、形状、配向であるファブリックを解析する処理とをコンピュータに行わせることを特徴とする粒子のファブリックの解析プログラム。
【請求項18】
算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルの粒子についてはファブリックの解析対象から除外する処理をコンピュータに行わせることを特徴とする請求項１７に記載の粒子のファブリックの解析プログラム。
【請求項19】
前記近似楕円体モデルの精度は、撮影により取得した断面画像での断面積と、前記断面画像に相当する前記近似楕円体モデルの断面積との誤差に基づいて検証することを特徴とする請求項１７又は１８に記載の粒子のファブリックの解析プログラム。
【請求項20】
前記近似楕円体モデルの精度は、Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、
Ｅｉの中央値Ｅｍｅｄ及び標準偏差Ｅｓｔｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外することを特徴とする請求項１９に記載の粒子のファブリックの解析プログラム。
【請求項21】
前記近似楕円体モデルは、楕円体の方程式で表され、前記複数の断面画像は、少なくとも７枚であることを特徴とする請求項１７から２０のうちのいずれか一つに記載の粒子のファブリックの解析プログラム。
【請求項22】
前記近似楕円体モデルは、非線形最小二乗法を利用して算出することを特徴とする請求項１７から２１のうちのいずれか一つに記載の粒子のファブリックの解析プログラム。
【請求項23】
前記近似楕円体モデルは、中心位置が前記複数の断面画像の重心位置に基づいて設定されることを特徴とする請求項１７から２２のうちのいずれか一つに記載の粒子のファブリックの解析プログラム。
【請求項24】
前記複数の断面画像は、ＣＴ画像であることを特徴とする請求項１７から２３のうちのいずれか一つに記載の粒子のファブリックの解析プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、砂や岩石などの集合体を構成する粒子や，砂や岩石の空隙もしくは鋳物等の金属材料に混入した気泡のファブリックの解析方法、解析装置、及び解析プログラムに関する。さらに詳述すると、本発明は、例えば集合体を非破壊で撮影した断面の画像に基づいて、抽出した粒子のファブリックを解析する技術に関する。
【0002】
本発明において、粒子のファブリックとは、砂や岩石などの集合体を構成する粒子や，砂や岩石の空隙もしくは鋳物等の金属材料に混入した気泡の大きさ、形状、配向を意味するものとする。
【背景技術】
【0003】
Ｘ線ＣＴは、Ｘ線源からＸ線を照射し、サンプル透過後のＸ線を検出器により測定する。この時、測定されるのは線吸収係数の２次元投影画像であるが、この２次元投影画像を様々な角度で取得し、コンピュータにおいて再構成計算を行うことにより３次元画像が得られ、物体内部の構造を理解することができる。Ｘ線ＣＴの最大の特徴は、サンプルを物理的に壊すことなく内部構造を可視化できる点にある。人体以外を撮影対象とする場合は、高出力のＸ線を照射できる産業用Ｘ線ＣＴや、焦点寸法が数μｍオーダーの高解像度のマイクロフォーカスＸ線ＣＴを活用し、工業製品の品質検査、文化財の内部解析等、その適用先は多岐に渡っている。
【0004】
地質学、地盤工学、岩盤工学、資源工学などの地盤や岩石を研究対象とする分野においてもＸ線ＣＴは活用されており、地質構造解析に関する研究、地盤や岩石の力学特性や水理特性解析に関する研究、石油生産技術の開発、ＣＣＳ事業におけるＣＯ２の長期安定性に関する研究等に活用が図られている。地質試料の観察や分析の多くは、試料を物理的に加工することが必要であるのに対し、Ｘ線ＣＴは非破壊で物体内部の３次元構造を可視化できる。つまり、ボーリング調査で得たサンプルのように、替えが効かない試料であっても、サンプルを壊すことなく、様々な角度の断面をコンピュータのディスプレイ上に表示できる。従って、従来行われている地質試料の観察や分析と組み合わせることにより、地質構造の解析を更に効率化、高精度化できる可能性があることから、ＣＴ画像の解析手法について更なる発展が期待される。
【0005】
ここで、地質試料内に含まれる粒子に着目し、ＣＴ画像内の粒子について、その大きさ、形状、配向に関するファブリックを３次元的に解析する手法について検討する。地質試料内の粒子に着目した研究として、地盤や岩盤の力学特性と粒子形状に着目した研究、地質構造と粒子の配向に着目した研究等が行われている。土質力学においては、土粒子の形状が地盤の力学特性に及ぼす影響が大きいことが知られており、土粒子形状と力学特性の関係について主に実験的な研究が行われている。堆積学においては、堆積運搬時の環境を推定するうえで堆積物粒子の配向は貴重な情報であり、堆積物粒子の配向解析が行われている。また、堆積岩の異方性を解析する上でも、粒子の配向は重要な意味を成す。また、断層や地すべりの解析を行う上でも、粒子の形状や配向は重要であると考えられる。断層岩における断層面付近や地すべりのすべり面付近には、運動時に破砕された岩片粒子が多く存在することが知られており、その粒子について形状や配向を解析することによって多くの情報を引き出すことができれば、断層や地すべりの解析を行う上でも役立つ可能性がある。
【0006】
このように、地質試料内に含まれる粒子をキーマーカとして抽出し、そのファブリックを解析することは非常に重要である。粒子のファブリックを直接的に解析するには、地質サンプルを解析したい断面で物理的に切断し、写真や顕微鏡を使って得た断面の画像を得た後、粒子の輪郭をなぞって抽出する方法が最も単純である。しかし、数多くの粒子を抽出するには多大な時間や労力が必要である。画像解析ソフトウェアを使うことで効率化も図れるが、３次元での解析を行うことは困難である。従って、３次元の内部構造を可視化できるＣＴ画像やＭＲＩ画像を解析する方法が最も効率的である。
【0007】
３次元で粒子のファブリックを解析する試みとしては、ＣＴ画像を対象とし、画像の解析手法に関する研究がいくつか行われている。例えば、Ｓｔａｒ ａｎａｌｙｓｅｓと呼ばれる手法を使って、粒子のファブリックを解析する手法が提案されている（非特許文献１）。Ｓｔａｒ ａｎａｌｙｓｅｓとは、ノギスを使って粒子を測定するように、解析対象とする粒子内の任意の点から各方向に線を引き、粒子境界までの線分の長さを直接測定し、粒子の短軸や長軸の向きを算定する方法である。また、他の方法として、ＳＰｒｉｎｇ－８のＸ線ＣＴを使って土粒子をＣＴ撮影し、個々の土粒子を抽出してメッシュモデルに変換した後、所定の解析方法に基づき作成された解析プログラムを用いて、土粒子の形状解析を行なう方法が開発されている（非特許文献２）。
【先行技術文献】
【特許文献】
【0008】
【文献】Ketcham, R.A. Three-dimensional grain fabric measurements using high-resolution X-ray computed tomography (2005) Journal of Structural Geology,
【文献】片桐 淳など、ＳＰｒｉｎｇ－８のＸ線マイクロＣＴを用いた３次元土粒子形状の定量的評価、土木学会論文集Ｃ（地圏工学）２０１４
【発明の概要】
【発明が解決しようとする課題】
【0009】
しかしながら、上述した非特許文献１に記載の方法では、個々の粒子についてファブリックを解析することができるが、ノギスを使って粒子の長さを直接測定するような方法であることから、粒子表面の凹凸の影響を受けやすく、精度が不十分であった。また、球や楕円体から大きく外れた粒子の解析は、精度を確保するために慎重に行わなければならず、処理が煩雑で長時間を要した。また、上述した非特許文献２に記載の方法では、ＣＴ画像をＳＴＬ形式のファイルに変換する必要があることから、変換に使用するアルゴリズム等の検討が別途必要であり、処理が煩雑で長時間を要した。
【0010】
上述した非特許文献2では，本来二つもしくはそれ以上の粒子であるにも関わらず一つの粒子として誤認識された粒子や，あるいは本来一つの粒子であるにも関わらず二つ以上の粒子として誤認識された粒子が混在しており，それらの粒子を目視で一つ一つ確認することにより，評価精度の向上を図っているが，その過程で多くの時間を要していた。
【0011】
そこで、本発明は、岩石や金属などの集合体を構成する粒子のファブリックを、非破壊で撮影した断面画像を利用して簡易、かつ、高精度に解析することができるファブリックの解析方法、解析装置、及び解析プログラムを提供することを目的とする。
【課題を解決するための手段】
【0012】
かかる目的を達成するため、本発明の粒子のファブリックの解析方法は、粒子の集合体を非破壊で撮影することにより複数の断面画像を取得し、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出し、算出した近似楕円体モデルに基づいて粒子の大きさ、形状、配向であるファブリックを解析するようにしている。
【0013】
また、本発明の粒子のファブリックの解析装置は、粒子の集合体を非破壊で撮影することにより複数の断面画像を取得する手段と、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出する手段と、算出した近似楕円体モデルに基づいて粒子の大きさ、形状、配向であるファブリックを解析する手段とを有するようにしている。
【0014】
また、本発明の粒子のファブリックの解析プログラムは、粒子の集合体を非破壊で撮影することにより複数の断面画像を取得する処理と、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出する処理と、算出した近似楕円体モデルに基づいて粒子の大きさ、形状、配向であるファブリックを解析する処理とをコンピュータに行わせるようにしている。
【0015】
したがって、これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムによると、砂や岩石などの集合体を構成する粒子や，砂や岩石の空隙もしくは鋳物等の金属材料に混入した気泡のファブリックを、非破壊で撮影した断面画像を利用して簡易、かつ、高精度に解析することができる。
【0016】
本発明の粒子のファブリックの解析方法は、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルの粒子についてはファブリックの解析対象から除外するようにしても良い。
【0017】
また、本発明の粒子のファブリックの解析装置は、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルの粒子についてはファブリックの解析対象から除外する手段を有するようにしても良い。
【0018】
また、本発明の粒子のファブリックの解析プログラムは、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルの粒子についてはファブリックの解析対象から除外する処理をコンピュータに行わせるようにしても良い。
【0019】
これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムの場合には、所定の精度を満たさない近似楕円体モデルの粒子についてはファブリックの解析対象から除外するので、解析をより高精度に実現することができる。
【0020】
本発明の粒子のファブリックの解析方法は、前記近似楕円体モデルの精度は、撮影により取得した断面画像での断面積と、前記断面画像に相当する前記近似楕円体モデルの断面積との誤差に基づいて検証するようにしても良い。
【0021】
また、本発明の粒子のファブリックの解析装置は、前記近似楕円体モデルの精度は、撮影により取得した断面画像での断面積と、前記断面画像に相当する前記近似楕円体モデルの断面積との誤差に基づいて検証するようにしても良い。
【0022】
また、本発明の粒子のファブリックの解析プログラムは、前記近似楕円体モデルの精度は、撮影により取得した断面画像での断面積と、前記断面画像に相当する前記近似楕円体モデルの断面積との誤差に基づいて検証するようにしても良い。
【0023】
これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムの場合には、断面画像と近似楕円体モデルとの比較により検証するので、解析をより高精度に実現することができる。
【0024】
本発明の粒子のファブリックの解析方法は、前記近似楕円体モデルの精度は、
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、中央値Ｅｍｅｄ及びＥｉの中央値Ｅｍｅｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外するようにしても良い。
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、中央値Ｅｍｅｄ及びＥｉの中央値Ｅｍｅｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外するようにしても良い。
【0025】
また、本発明の粒子のファブリックの解析装置は、前記近似楕円体モデルの精度は、
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、
中央値Ｅｍｅｄ及びＥｉの中央値Ｅｍｅｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外するようにしても良い。
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、
中央値Ｅｍｅｄ及びＥｉの中央値Ｅｍｅｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外するようにしても良い。
【0026】
また、本発明の粒子のファブリックの解析プログラムは、前記近似楕円体モデルの精度は、
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、中央値Ｅｍｅｄ及びＥｉの中央値Ｅｍｅｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外するようにしても良い。
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ＳＳＶ：撮影により取得した断面画像での断面積
Ｓｅｐ：断面画像に相当する近似楕円体モデルの断面積
Ｅｉ：断面積ＳＳＶと断面積Ｓｅｐとの誤差に関する値
ｉ：断面数
とした場合に、中央値Ｅｍｅｄ及びＥｉの中央値Ｅｍｅｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外するようにしても良い。
【0027】
これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムの場合には、誤差に関する値の標準偏差を利用しているので、解析をより高精度に実現することができる。
【0028】
本発明の粒子のファブリックの解析方法は、前記近似楕円体モデルは、楕円体の方程式で表され、前記複数の断面画像は、少なくとも７枚であるようにしても良い。
【0029】
また、本発明の粒子のファブリックの解析装置は、前記近似楕円体モデルは、楕円体の方程式で表され、前記複数の断面画像は、少なくとも７枚であるようにしても良い。
【0030】
また、本発明の粒子のファブリックの解析プログラムは、前記近似楕円体モデルは、楕円体の方程式で表され、前記複数の断面画像は、少なくとも７枚であるようにしても良い。
【0031】
これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムの場合には、楕円体の方程式では未知数が７つであることから、７枚の断面画像によって全ての係数を算出することができる。
【0032】
本発明の粒子のファブリックの解析方法は、前記近似楕円体モデルは、非線形最小二乗法を利用して算出するようにしても良い。
【0033】
また、本発明の粒子のファブリックの解析装置は、前記近似楕円体モデルは、非線形最小二乗法を利用して算出するようにしても良い。
【0034】
また、本発明の粒子のファブリックの解析プログラムは、前記近似楕円体モデルは、非線形最小二乗法を利用して算出するようにしても良い。
【0035】
これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムの場合には、非線形最小二乗法を利用しているので、より高精度に解析を行うことができる。
【0036】
本発明の粒子のファブリックの解析方法は、前記近似楕円体モデルは、中心位置が前記複数の断面画像の重心位置に基づいて設定されるようにしても良い。
【0037】
また、本発明の粒子のファブリックの解析装置は、前記近似楕円体モデルは、中心位置が前記複数の断面画像の重心位置に基づいて設定されるようにしても良い。
【0038】
また、本発明の粒子のファブリックの解析プログラムは、前記近似楕円体モデルは、中心位置が前記複数の断面画像の重心位置に基づいて設定されるようにしても良い。
【0039】
これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムの場合には、近似楕円体モデルの重心位置を高精度に設定することができる。
【0040】
本発明の粒子のファブリックの解析方法は、前記複数の断面画像は、ＣＴ画像であるようにしても良い。
【0041】
また、本発明の粒子のファブリックの解析装置は、前記複数の断面画像は、ＣＴ画像であるようにしても良い。
【0042】
また、本発明の粒子のファブリックの解析プログラムは、前記複数の断面画像は、ＣＴ画像であるようにしても良い。
【0043】
これらの粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムの場合には、ＣＴ画像を用いて、解析により適した断面画像を得ることができる。
【発明の効果】
【0044】
本発明の粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムによれば、岩石や金属などの集合体を構成する粒子のファブリックを、非破壊で撮影した断面画像を利用して簡易、かつ、高精度に解析することが可能になる。
【図面の簡単な説明】
【0045】
【図1】本発明の粒子のファブリックの解析方法の実施形態の一例を説明するフローチャートである。
【図2】実施形態の粒子のファブリックの解析方法を粒子のファブリックの解析プログラムを用いて実施する場合の当該プログラムによって実現される粒子のファブリックの解析装置の機能ブロック図である。
【図3】ブロックサンプルの撮影の態様の例を示す模式図である。
【図4】切断面を示す座標系を示す説明図である。
【図5】実施形態の粒子の識別方法の一例を説明するフローチャートである。
【図6】ブロックサンプルの断面の一例を撮影したＣＴ画像である。
【図7】（ａ）はブロックサンプルの断面のサブボリューム画像であり、（ｂ）はそれを二値化処理した画像である。
【図8】複数の断面を抽出する状態を示す模式図である。
【図9】実施形態の評価解析手順の前半を説明するフローチャートである。
【図10】実施形態の評価解析手順の後半を説明するフローチャートである。
【図11】実施例の第１グループ～第４グループの数値的に作成した粒子を示すＣＴ画像である。
【図12】実施例の第５グループ及び第６グループの数値的に作成した粒子を示すＣＴ画像である。
【図13】実施例の第７グループの数値的に作成した粒子を示すＣＴ画像である。
【図14】実施例１の粒子の切断面ごとのＥｉの分布を示すグラフである。
【図15】実施例２の粒子の切断面ごとのＥｉの分布を示すグラフである。
【図16】比較例１の粒子の切断面ごとのＥｉの分布を示すグラフである。
【図17】比較例２の粒子の切断面ごとのＥｉの分布を示すグラフである。
【図18】近似楕円体の長軸と単位球の交点を計算し、長軸の走向角と傾斜角をシュミットネット下半球に投影した図であり、（ａ）は実施例４及び実施例５、（ｂ）は実施例６及び実施例７、（ｃ）は比較例３及び比較例４である。
【図19】実施例８の粒子と近似楕円体モデルを３次元表示した図である。
【図20】実施例８の粒子の切断面ごとのＥｉの分布を示すグラフである。
【図21】比較例５の粒子と近似楕円体モデルを３次元表示した図である。
【図22】比較例５の粒子の切断面ごとのＥｉの分布を示すグラフである。
【図23】比較例５の粒子と近似楕円体モデルの中央部の切断面を示す図である。
【図24】砂粒子の粒径加積曲線を示すグラフである。
【図25】砂粒子を容れた容器を撮影したＣＴ画像であり、（ａ）は水平断面、（ｂ）は鉛直断面、（ｃ）は水平断面の粒子識別結果、（ｄ）は鉛直断面の粒子識別結果である。
【図26】グレイバリューのヒストグラムである。
【図27】評価解析が収束した粒子のＥｍｅｄに関するヒストグラムである。
【図28】粒子のＣＴ画像と近似楕円体モデルの画像であり、（ａ）はＰ５６８８、（ｂ）はＰ２８６３、（ｃ）はＰ７８１１、（ｄ）は２３２１、（ｅ）は５４３０である。
【図29】粒子のＣＴ画像と近似楕円体モデルの画像であり、（ａ）はＰ９０２８、（ｂ）はＰ５９１５、（ｃ）はＰ７３０４、（ｄ）は５９７６である。
【図30】砂粒子の形状評価結果であり、（ａ）は楕円体の軸半径、（ｂ）はアスペクト比、細長比、扁平度を示すグラフである。
【図31】砂粒子の形状評価結果であり、（ａ）はＫｒｕｍｂｅｉｎの球形度、（ｂ）は傾斜角を示すグラフである。
【図32】砂粒子の形状評価結果であり、Ｚｉｎｇｇ Ｄｉａｇｒａｍのヒストグラムコンター図である。
【発明を実施するための形態】
【0046】
以下、本発明の構成を図面に示す実施の形態の一例に基づいて詳細に説明する。
【0047】
【0048】
本実施形態では、集合体の一例として、アクリル容器に充填した砂のＣＴ撮影を行い、砂粒子のファブリックを解析した場合を例に挙げて説明する。ここでは、粒径が３００～６００μｍを主体とする砂を、内径１０ｍｍの円筒形状のアクリル容器に自然落下法により充填した。
【0049】
本実施形態の粒子のファブリックの解析方法は、図１に示すように、粒子の集合体から非破壊で複数の断面画像を取得し（Ｓ１）、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出し（Ｓ２）、算出した近似楕円体モデルに基づいて粒子のファブリックを解析する（Ｓ４）ようにしている。また、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルについては粒子のファブリックの解析対象から除外する（Ｓ３）ようにしている。
【0050】
上記の粒子のファブリックの解析方法は、本実施形態に係る粒子のファブリックの解析装置によって実施され得る。本実施形態の粒子のファブリックの解析装置は、例えば図２に示すよう、粒子の集合体から非破壊で複数の断面画像を取得する手段（１１ａ）と、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出する手段（１１ｂ）と、算出した近似楕円体モデルに基づいて粒子のファブリックを解析する手段（１１ｄ）とを有するようにしている。また、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルについては粒子のファブリックの解析対象から除外する手段（１１ｃ）を有するようにしている。
【0051】
上記の粒子のファブリックの解析方法や粒子のファブリックの解析装置は、粒子のファブリックの解析プログラムがコンピュータ上で実行されることによっても実施されたり実現されたりし得る。本実施形態に係る説明では、粒子のファブリックの解析プログラムがコンピュータ上で実行されることによって粒子のファブリックの解析方法が実施されると共に粒子のファブリックの解析装置が実現される場合を取り上げて説明する。
【0052】
本実施形態の粒子のファブリックの解析プログラム１７を実行するためのコンピュータ１０（本実施形態では、粒子のファブリックの解析装置１０でもある）の全体構成を図２に示す。
【0053】
このコンピュータ１０（以下、「粒子のファブリックの解析装置１０」と表記する）は制御部１１，記憶部１２，入力部１３，表示部１４，及びメモリ１５を備え、これらが相互にバス等の信号回線によって接続されている。
【0054】
制御部１１は、記憶部１２に記憶されている粒子のファブリックの解析プログラム１７に従って粒子のファブリックの解析装置１０全体の制御並びにファブリックの解析に係る演算を行うものであり、例えばＣＰＵ（即ち、中央演算処理装置）である。
【0055】
記憶部１２は、少なくともデータやプログラムを記憶可能な装置であり、例えばハードディスクやソリッドステートドライブ（SSD;Solid State Drive）である。
【0056】
入力部１３は、少なくとも作業者の命令や種々の情報を制御部１１に与えるためのインターフェイス（即ち、情報入力の仕組み）であり、例えばキーボードやマウスである。なお、例えばキーボードとマウスとの両方のように複数種類のインターフェイスを入力部１３として有するようにしても良い。
【0057】
表示部１４は、制御部１１の制御によって文字や図形或いは画像等の描画・表示を行うものであり、例えばディスプレイである。
【0058】
メモリ１５は、制御部１１が種々の制御や演算を実行する際の作業領域であるメモリ空間となるものであり、例えばＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ の略）である。
【0059】
また、粒子のファブリックの解析装置１０に、必要に応じ、当該解析装置１０との間でデータや制御指令等の信号の送受信（即ち、出入力）が可能であるように、バスや広域ネットワーク回線等の信号回線により、データサーバ２０が接続されるようにしても良い。また、コンピュータ１０は、必要に応じ、インターネットなどのネットワークを介してクラウドサーバ（図示していない）にアクセス可能であるようにしても良い。
【0060】
なお、粒子のファブリックの解析装置１０とＣＴ装置（図２には図示していない）とが接続され、ＣＴ装置からの出力が粒子のファブリックの解析装置１０に直接入力されるようにしても良い。本実施形態では、粒子のファブリックの解析装置１０と産業用のＣＴ装置３０とが接続されている（図３参照）。
【0061】
そして、粒子のファブリックの解析装置１０の制御部１１には、粒子のファブリックの解析プログラム１７が実行されることにより、粒子の集合体から非破壊で複数の断面画像を取得する取得部１１ａと、取得した複数の断面画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出する算出部１１ｂと、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルについては粒子のファブリックの解析対象から除外する検証部１１ｃと、算出した近似楕円体モデルに基づいて粒子のファブリックを解析する解析部１１ｄと、が構成される。
【0062】
〈集合体の撮影〉
粒子のファブリックの解析方法の実施に係る手順として、まず、ＣＴ装置３０によりブロックサンプル３１の撮影が行われる（Ｓ１）。この処理では、ブロックサンプル３１の全体のＣＴ画像が取得される。
粒子のファブリックの解析方法の実施に係る手順として、まず、ＣＴ装置３０によりブロックサンプル３１の撮影が行われる（Ｓ１）。この処理では、ブロックサンプル３１の全体のＣＴ画像が取得される。
【0063】
〈近似楕円体モデルの算出〉
次に、Ｓ１の処理によって撮影されて取得された画像データが用いられて、粒子の近似楕円体モデルの算出が行われる（Ｓ２）。ここでの近似楕円体モデルの算出の手順としては、まず、算出対象となる粒子を識別し、次にその粒子について近似楕円体モデルとなる楕円体の方程式を算出する。
次に、Ｓ１の処理によって撮影されて取得された画像データが用いられて、粒子の近似楕円体モデルの算出が行われる（Ｓ２）。ここでの近似楕円体モデルの算出の手順としては、まず、算出対象となる粒子を識別し、次にその粒子について近似楕円体モデルとなる楕円体の方程式を算出する。
【0064】
ここで、近似楕円体モデルについて説明する。原点を中心とする楕円体の一般式は、数式１で表される。ここに、Ａ、Ｂ、Ｃ、Ｄ、Ｆ、Ｇ、Ｈは係数である。つまり、３次元空間上で楕円体を定義するには、７つのパラメータを評価することが必要である。
（数１）
Ａｘ２＋Ｂｙ２＋Ｃｚ２＋２Ｆｙｚ＋２Ｇｚｘ＋２Ｈｘｙ＋Ｄ＝０
（数１）
Ａｘ２＋Ｂｙ２＋Ｃｚ２＋２Ｆｙｚ＋２Ｇｚｘ＋２Ｈｘｙ＋Ｄ＝０
【0065】
一方、図４に示す座標系で、θとφを法線ベクトル４１ｎの角度とした場合に、原点を通り法線ベクトル４１ｎがｎｐｘ、ｎｐｙ、ｎｐｚである平面の方程式は、数式２で表される。
（数２）
ｎｐｘｘ＋ｎｐｙｙ＋ｎｐｚｚ＝０
但し、ｎｐｘ＝ｓｉｎθ・ｃｏｓφ
ｎｐｙ＝ｓｉｎθ・ｓｉｎφ
ｎｐｚ＝ｃｏｓθ
（数２）
ｎｐｘｘ＋ｎｐｙｙ＋ｎｐｚｚ＝０
但し、ｎｐｘ＝ｓｉｎθ・ｃｏｓφ
ｎｐｙ＝ｓｉｎθ・ｓｉｎφ
ｎｐｚ＝ｃｏｓθ
【0066】
図４に示すように、数式１で表される楕円体４０を、数式２の切断面４１で切断した時、切断面４１の断面積Ｓは数式３で算出される。
但し、ｅ１＝（ＢＣ－Ｆ２）ｎｐｘ
２
ｅ２＝（ＡＣ－Ｇ２）ｎｐｙ ２
ｅ３＝（ＡＢ－Ｈ２）ｎｐｚ ２
ｅ４＝２（ＦＧ－ＣＨ）ｎｐｘｎｐｙ
ｅ５＝－２（ＢＧ－ＦＨ）ｎｐｘｎｐｚ
ｅ６＝－２（ＡＦ－ＧＨ）ｎｐｙｎｐｚ
【数3】
ｅ２＝（ＡＣ－Ｇ２）ｎｐｙ ２
ｅ３＝（ＡＢ－Ｈ２）ｎｐｚ ２
ｅ４＝２（ＦＧ－ＣＨ）ｎｐｘｎｐｙ
ｅ５＝－２（ＢＧ－ＦＨ）ｎｐｘｎｐｚ
ｅ６＝－２（ＡＦ－ＧＨ）ｎｐｙｎｐｚ
【0067】
本実施形態では、ＣＴ画像から抽出した粒子の断面積を複数計算し、数式３で示される近似楕円体の断面積と比較することで、数式１中の７つの係数を算出する。これにより、粒子の近似楕円体モデルの算出が行われる。
【0068】
次に、算出対象となる粒子を識別する方法としては、２通りの方法が考えられる。１つは、領域成長、収縮／膨張フィルタ処理、オープニング／クロージングフィルタ処理のような画像解析ツールを使い、粒子１つ１つを手動で識別する方法である。この方法は、粒子の識別状況をＰＣ画面上で随時確認しながら行うため、解析対象とする粒子を高精度に識別することができる。一方で、粒子の識別には多くの労力・時間が必要であるため、抽出できる粒子の個数には限界がある。また、識別する粒子の形状や場所の選定には、少なからず主観が入ってしまう。
【0069】
もう１つの方法は、グレイバリューの範囲を指定して等値面を定義して関心領域を作成し、粒子を識別する方法である。この方法は、非常に簡単であり、短時間で多くの粒子を識別することができる。また、客観的に粒子を識別できることも特徴である。但し、粒子の識別精度は前者の方法に比べて悪く、本来２つもしくはそれ以上の粒子であるにも関わらず１つの粒子として識別してしまう場合や、本来１つの粒子であるにも関わらず２つ以上の粒子として識別してしまう可能性がある。目視によってそれら識別精度が悪い粒子を見分けることも可能であるが、粒子数が多いと確認作業に時間がかかってしまうことから、何らかの方法で識別精度が悪い粒子を見分けることが望まれる。そこで、本実施形態では、グレイバリューの範囲を指定して等値面を作成し、短時間でなるべく多くの粒子を識別すると共に、検証によって識別精度の悪い粒子を排除するようにする。
【0070】
算出対象となる粒子の識別の具体的な手順について、以下に説明する。まず、図５に示すように、ＣＴ画像を撮影し画像を取得する（Ｓ１０）。取得した画像の一例を図６に示す。そして、特定の粒子に着目し（図６中の円で囲んだ部分）、着目した粒子を包含するような代表体積要素をサブボリューム画像として抽出する（Ｓ１１）。サブボリューム画像の一例を図７（ａ）に示す。更に、サブボリューム画像に対して二値化処理を行って、粒子のみを抽出する（Ｓ１２）。二値化処理を行った粒子の画像の一例を図７（ｂ）に示す。
【0071】
続けて、識別した粒子の断面積の計算を行う（Ｓ１３）。本実施形態では、θ、φともに０～１８０度の範囲の中で、１５度毎に断面を抽出した（図８参照）。また、θ＝０度の場合は、φ＝０～１８０度のいずれの場合も同じ面になるので、同じ面が１２枚含まれてしまう。そこで、そのうちの１１枚を含めないようにするために、１２枚×１２－１１枚＝１３３枚として、粒子１個あたり計１３３組の断面積を計算する。
【0072】
次に、識別した粒子について計算された断面積データに基づいて、近似楕円体モデルとなる楕円体の方程式を算出する手順について、以下に説明する。図９及び図１０に示すように、本実施形態では、ガウスニュートン法に基づく最適化計算を実施し、数式１の係数を解析する。即ち、近似楕円体モデルは、ガウスニュートン法に基づく非線形最小二乗法を利用して算出する。評価解析は、最小反復計算回数を１５回、最大反復計算回数を１００回とする。評価解析では、まず初めに、第１段階の初期値として、軸半径は粒子を球近似した場合の球半径を、回転軸ｎｅｘ、ｎｅｙ、ｎｅｚについてはＺ軸を、回転角度θｅについては０度を与える（Ｓ２０）。
【0073】
ここでの解析における目的関数は、数式４の通りである。但し、ｆｏｂｊは目的関数、ＳＳＶはＣＴ画像から抽出した粒子の断面積、Ｓｅｐは近似楕円体の断面積をそれぞれ示す。
【数4】
【0074】
そして、ガウスニュートン法に基づき、残差行列の計算及びヤコビアン行列の計算を実行する（Ｓ２１）。このとき、反復計算回数のカウンタｉｔｒを１アップする。反復計算を進めるごとに、パラメータを算出して更新していく（Ｓ２２）。反復計算の過程で、数式３の右辺分母の平方根内が負になると計算が進行しなくなってしまう。このため、その時点でのパラメータを利用して、数式３の右辺分母の平方根内のｅが正の数であるか否かを判断する（Ｓ２３）。ｅが正の数であると判断された場合は（Ｓ２３のＹＥＳ）、反復計算が１５回以上行われたか否かを判断する（Ｓ２４）。反復計算が１５回以上行われていないと判断した場合は（Ｓ２４のＮＯ）、最小反復計算回数に満たないので、再び残差行列の計算及びヤコビアン行列の反復計算を実行する（Ｓ２１）。
【0075】
反復計算が１５回以上行われていると判断した場合は（Ｓ２４のＹＥＳ）、数式５で表されるＥｉを計算する（Ｓ２５）。Ｅｉは、法線ベクトルがｎｐｘ，ｉ、ｎｐｙ，ｉ、ｎｐｚ，ｉである断面で計算された断面積と近似楕円体との誤差に関する指標である。
（数５）
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ｉ＝１～１３３
（数５）
Ｅｉ＝｜１．０－ＳＳＶ／Ｓｅｐ｜ｉ
但し、ｉ＝１～１３３
【0076】
ここで、Ｅｉについては、反復計算回数が１５回を超えた場合に、５回前までのＥｉの平均値Ｅｉｔ＝ｋ～ｋ－５と、現在の計算過程において計算されたＥｉの平均値Ｅｉｔ＝ｋとの差分ΔＥｉを算出する。そして、その差分ΔＥｉが１．０×１０－５以下であるか否かを判断し（Ｓ２６）、パラメータの収束の判定を行う。即ち、それ以上、反復計算をすることでパラメータが殆ど変わらない程度に収束しているか否かを判断する。その差分ΔＥｉが１．０×１０－５以下でないと判断した場合は（Ｓ２６のＮＯ）、パラメータは収束していないと判断して、再び残差行列の計算及びヤコビアン行列の反復計算を実行する（Ｓ２１）。その差分ΔＥｉが１．０×１０－５以下であると判断した場合は（Ｓ２６のＹＥＳ）、パラメータは収束しており反復計算を続けても殆ど変わらないと判断する。
【0077】
そして、反復計算が１００回以下であるか否かを判断する（Ｓ２７）。反復回数が１００回以下であると判断された場合は（Ｓ２７のＹＥＳ）、評価解析は最大反復計算回数１００回以下という条件を満たしているので、そのときのパラメータに決定し（Ｓ２８）、処理を終了する。尚、近似楕円体モデルは、中心位置が複数の断面画像の重心位置に基づいて設定される。
【0078】
一方、Ｓ２３でｅが正の数でないと判断した場合（Ｓ２３のＮＯ）、あるいは、Ｓ２７で反復回数が１００回以下でないと判断した場合は（Ｓ２７のＹＥＳ）、第１段階の初期値では最大反復計算回数１００回以下の条件でパラメータが収束しなかったということで、第１段階の初期値を第２段階の初期値に変更して、近似楕円体モデルとなる楕円体の方程式を再び算出する（Ｓ２）。第２段階の初期値として、軸半径はＸ、Ｙ、Ｚ軸方向に粒子を投影した場合の辺長さを、回転軸と回転角度は第１段階と同じ値を与える（Ｓ３０）。
【0079】
そして、第１段階と同様に、ガウスニュートン法に基づき、残差行列の計算及びヤコビアン行列の計算を実行する（Ｓ３１）。このとき、反復計算回数のカウンタｉｔｒを１アップする。反復計算を進めるごとに、パラメータを算出して更新していく（Ｓ３２）。その時点でのパラメータを利用して、数式３の右辺分母の平方根内のｅが正の数であるか否かを判断する（Ｓ３３）。ｅが正の数であると判断された場合は（Ｓ３３のＹＥＳ）、反復計算が１５回以上行われたか否かを判断する（Ｓ３４）。反復計算が１５回以上行われていないと判断した場合は（Ｓ３４のＮＯ）、最小反復計算回数に満たないので、再び残差行列の計算及びヤコビアン行列の反復計算を実行する（Ｓ３１）。
【0080】
反復計算が１５回以上行われていると判断した場合は（Ｓ３４のＹＥＳ）、数式５で表されるＥｉを計算する（Ｓ３５）。Ｅｉについては、第１段階と同様に、差分ΔＥｉを算出し、その差分ΔＥｉが１．０×１０－５以下であるか否かを判断し（Ｓ３６）、パラメータの収束の判定を行う。その差分ΔＥｉが１．０×１０－５以下でないと判断した場合は（Ｓ３６のＮＯ）、パラメータは収束していないと判断して、再び残差行列の計算及びヤコビアン行列の反復計算を実行する（Ｓ３１）。その差分ΔＥｉが１．０×１０－５以下であると判断した場合は（Ｓ３６のＹＥＳ）、パラメータは収束しており反復計算を続けても殆ど変わらないと判断する。
【0081】
そして、反復計算が１００回以下であるか否かを判断する（Ｓ３７）。反復回数が１００回以下であると判断された場合は（Ｓ３７のＹＥＳ）、評価解析は最大反復計算回数１００回以下という条件を満たしているので、そのときのパラメータに決定し（Ｓ３８）、処理を終了する。この場合も、近似楕円体モデルは、中心位置が複数の断面画像の重心位置に基づいて設定される。
【0082】
一方、Ｓ３３でｅが正の数でないと判断した場合は（Ｓ３３のＮＯ）、計算が進行しなくなるので、処理を終了する。また、Ｓ３７で反復回数が１００回以下でないと判断した場合は（Ｓ３７のＹＥＳ）、第２段階の初期値であっても最大反復計算回数１００回以下の条件でパラメータが収束しなかったということで、処理を終了する。即ち、本実施形態では、近似楕円体モデルの精度は、撮影により取得した断面画像での断面積ＳＳＶと、断面画像に相当する近似楕円体モデルの断面積Ｓｅｐとの誤差（数式５参照）に基づいて検証する。
【0083】
上述したように、本実施形態では、図９及び図１０に示すように、Ｓ１の処理によって撮影されて取得された画像データが用いられて、第１段階の初期値を利用して粒子の近似楕円体モデルの算出が行われ（Ｓ２、Ｓ２０～Ｓ２４）、第１段階で除外された粒子の中から、第２段階の初期値を利用して粒子の近似楕円体モデルの算出が行われる（Ｓ２、Ｓ３０～Ｓ３４）。即ち、粒子の条件によっては、２段階に処理が行われる場合があるが、これには限られない。例えば、１段階のみとしたり、更に条件を変えて３段階以上の処理を実行するようにしてもよい。
【0084】
〈精度の検証〉
次に、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルについては粒子のファブリックの解析対象から除外する（Ｓ３）。所定の精度については、本実施形態では、粒子の半径の大きさ、Ｅｉの中央値Ｅｍｅｄ，Ｅｉの標準偏差Ｅｓｔｄを採用している。
次に、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルについては粒子のファブリックの解析対象から除外する（Ｓ３）。所定の精度については、本実施形態では、粒子の半径の大きさ、Ｅｉの中央値Ｅｍｅｄ，Ｅｉの標準偏差Ｅｓｔｄを採用している。
【0085】
まず、粒子の半径については、半径が５ｖｏｘｅｌ以上である粒子を解析対象とするようにし、半径が５ｖｏｘｅｌ未満の粒子については除外する。これは、評価対象とする粒子を球近似した場合の半径が５ｖｏｘｅｌより小さい場合、部分体積効果の影響によって評価精度が低下する可能性があるためである。但し、撮影条件などが異なれば、必ずしも５ｖｏｘｅｌ以上であることには限られない。
【0086】
また、Ｅｉの中央値Ｅｍｅｄ及び標準偏差Ｅｓｔｄについては、中央値Ｅｍｅｄ及びＥｓｔｄの少なくとも一方が０．１以上である近似楕円体モデルの粒子についてはファブリックの解析対象から除外する。これは、中央値Ｅｍｅｄ及び標準偏差Ｅｓｔｄの少なくとも一方が０．１以上の場合は、複数の粒子を単一の粒子として識別してしまった場合、あるいは単一の粒子を複数の粒子として識別してしまった場合、評価解析の計算過程で解が得られない場合などが含まれる可能性が高いためである。但し、他の条件などが異なれば、解析対象から除外するのは中央値Ｅｍｅｄ及び標準偏差Ｅｓｔｄの少なくとも一方が０．１以上である場合には限られない。
【0087】
〈ファブリックの解析〉
その後は、算出した近似楕円体モデルに基づいて粒子のファブリックを解析する（Ｓ４）。ここでの解析手法は、適宜な手法を適用することができる。
その後は、算出した近似楕円体モデルに基づいて粒子のファブリックを解析する（Ｓ４）。ここでの解析手法は、適宜な手法を適用することができる。
【0088】
従来から知られる粒子のファブリックの解析方法では、粒子の大きさを直接測定して粒子の形状や配向を評価する方法であるため、粒子表面の微小な凹凸が評価結果に及ぼす影響が大きかった。一方で、本実施形態で示した手法は、粒子の断面積に基づき楕円体として近似することから、粒子表面の形状の影響が緩和されると考えられる。また、あらゆる方向の断面積のデータ群から粒子を楕円体近似するため、複雑な形状の粒子であっても、平均的な楕円体として近似できる可能性がある。
【0089】
上述したように、粒子の識別方法としては主に２通りの方法が考えられるが、本実施形態では、グレイバリューの範囲を指定して等値面を作成し、短時間でなるべく多くの粒子を識別する。但し、識別された粒子のデータの中には、識別精度の悪い粒子が混在する。そこで、以下の方法で粒子を分類することができると考えた。第１に、複数の粒子を単一の粒子として識別してしまった場合、あるいは単一の粒子を複数の粒子として識別してしまった場合、その粒子は近似対象となる楕円体から大きく逸脱した形状となると考えられる。そのような粒子は、計算の過程で数式３の右辺分母の平方根内が負となり、計算が収束しない場合がある。従って、そのような粒子については、識別精度が悪い粒子と判断し、評価から除外する。第２に、計算が収束した場合であっても、数式５により計算されるＥｉを検証することにより、識別精度が悪い粒子を区別できる。つまり、Ｅｉは、ＣＴ画像から抽出した粒子と近似楕円体の比を表しており、抽出された粒子の断面積と近似楕円体の断面積の差が小さいほどＥｉは０に近づき、その差が大きいほどＥｉは大きくなる。従って、評価解析を行った後、Ｅｉの検証を行うことで、識別精度が悪い粒子を区別できる。
【0090】
以上のように構成された粒子のファブリックの解析方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムによれば、複数のＣＴ画像に基づいて粒子を楕円体に近似させた近似楕円体モデルを算出し、算出した近似楕円体モデルに基づいてファブリックを解析するので、岩石や金属などの集合体を構成する粒子のファブリックを、非破壊で撮影した断面画像を利用して簡易、かつ、高精度に解析することができる。特に、本実施形態の観察方法，粒子のファブリックの解析装置，粒子のファブリックの解析プログラムでは、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルについては粒子のファブリックの解析対象から除外するので、粒子のファブリックをより高精度に解析することができる。
【0091】
なお、上述の実施形態は本発明を実施する際の好適な形態の一例ではあるものの本発明の実施の形態が上述のものに限定されるものではなく、本発明の要旨を逸脱しない範囲において本発明は種々変形実施可能である。
【0092】
例えば、上述の実施形態では、算出した近似楕円体モデルの精度を検証し、所定の精度を満たさない近似楕円体モデルについては粒子のファブリックの解析対象から除外する場合について説明したが、これには限られず、算出した近似楕円体モデルの精度は必ずしも検証しなくてもよい。
【0093】
また、上述の実施形態では、近似楕円体モデルは、非線形最小二乗法を利用して算出する場合について説明したが、これには限られず、他の手法を利用して算出するようにしてもよい。
【0094】
また、上述の実施形態では、複数の断面画像としてＣＴ画像を適用した場合について説明したが、これには限られず、ＭＲＩ画像など、非破壊で透過画像を得られる撮影方法であれば他の撮影方法であってもよい。
【0095】
また、上述の実施形態では粒子のファブリックの解析プログラム１７がコンピュータ１０に予めインストールされているようにしているが、粒子のファブリックの解析プログラム１７がコンピュータ読み取り可能な記録媒体に格納されて提供されるようにしても良い。
【実施例】
【0096】
【0097】
本実施例では、数値的に作成するＣＴ画像は、画素数を５０×５０×５０ｖｏｘｅｌ、画素サイズを１ｍｍ／ｐｉｘｅｌとし、画像中心に粒子を配置した。粒子に該当する部分の画素値は２０００、それ以外の画素値は０とした。第１グループ～第４グループは形状が楕円体である粒子、第５グループ及び第６グループは形状が矩形である粒子、第７グループは２つの楕円体粒子が結合した粒子を表す。ＣＴ画像の作成条件を表１に示す。なお、楕円体の長軸の軸半径をｒ１、中間軸の軸半径をｒ２、短軸の軸半径をｒ３とした。また、軸の方向は、走向角（ｔｒｅｎｄ、方位角）と傾斜角（ｐｌｕｎｇｅ、傾き）とで示した。走向角は傾いた直線を水平面に投影したときの方向（方位角）であり、傾斜角は水平面からの傾きである。そして、楕円体の長軸の走向角をｔｒ１、中間軸の走向角をｔｒ２、短軸の走向角をｔｒ３、楕円体の長軸の傾斜角をｐｌ１、中間軸の傾斜角をｐｌ２、短軸の傾斜角をｐｌ３とした。
【表1】
【0098】
第１グループ～第４グループは、図１１に示すように、楕円体の粒子である。第１グループは、楕円体の軸半径について検証するために実施し、中間軸と短軸を各６通りずつ、計３６通りについて評価を行った。第２グループ～第４グループは、楕円体の長軸の走向角と傾斜角について検証するために実施した。それぞれ短軸の軸半径ｒ３が異なっており、走向角が１２通り、傾斜角が５通りの計６０通り実施した。第５グループ～第６グループは、第１グループ～第４グループとは粒子の形状が異なっており、図１２に示すように、１辺の長さがそれぞれ３０、２０、１０ｖｏｘｅｌである矩形形状の粒子であり、各軸の長さは各辺の半分の長さとなる。第６グループは、第５グループの矩形粒子についてＺ軸を中心に反時計回りに４５度回転させた粒子である。第７グループは、図１３に示すように、２つの楕円体が結合して単一の粒子を構成する粒子であり、楕円体の形状はそれぞれ同じであるが、中間軸方向に楕円体が結合する粒子である。
【0099】
以下、各グループごとに、条件を異ならせて実施例と比較例とを設定し、解析結果を比較した。方法としては、抽出粒子と近似楕円体モデルとのＥｉを算出し、その中央値Ｅｍｅｄと標準偏差Ｅｓｔｄとを算出することで、適合性を評価した。
【0100】
〈第１グループ〉
第１グループの解析結果（全３６組）のうち、軸半径の評価誤差が最も大きい粒子（実施例１）と、数式５により計算されるＥｉの中央値Ｅｍｅｄと標準偏差Ｅｓｔｄが最も大きい粒子（実施例２）の評価結果を表２に示す。実施例１のＥｉの計算結果を図１４に示し、実施例２のＥｉの計算結果を図１５に示す。図１４及び図１５は、横軸が断面積を計算した切断面の番号（１３３組）であり、縦軸がＥｉを示す。
第１グループの解析結果（全３６組）のうち、軸半径の評価誤差が最も大きい粒子（実施例１）と、数式５により計算されるＥｉの中央値Ｅｍｅｄと標準偏差Ｅｓｔｄが最も大きい粒子（実施例２）の評価結果を表２に示す。実施例１のＥｉの計算結果を図１４に示し、実施例２のＥｉの計算結果を図１５に示す。図１４及び図１５は、横軸が断面積を計算した切断面の番号（１３３組）であり、縦軸がＥｉを示す。
【表2】
【0101】
表２に示した通り、実施例１では、軸半径の最大誤差は短軸の軸半径ｒ３の０．５ｖｏｘｅｌであり、精度良く評価することができた。また、実施例１では、Ｅｍｅｄは０．０１１程度、Ｅｓｔｄは０．０１６程度であり、こちらについても十分に小さい結果となった。一方、実施例２では、Ｅｍｅｄは最大で０．０２１、Ｅｓｔｄは最大で０．０２７であり、こちらについても十分に小さい結果が得られた。
【0102】
次に、長軸の軸半径ｒ１を１５ｖｏｘｅｌ、中間軸ｒ２の軸半径を１１ｖｏｘｅｌとし、短軸の軸半径ｒ３を５（実施例３）、３（比較例１）、１（比較例２）ｖｏｘｅｌとした場合の評価結果を、表３に示す。比較例１のＥｉの計算結果を図１６に示し、比較例２のＥｉの計算結果を図１７に示す。実施例３では、短軸の軸半径ｒ３が５ｖｏｘｅｌであり、評価された軸半径も設定値と概ね等しく、Ｅｍｅｄ、Ｅｓｔｄも十分に小さい結果が得られた。一方、比較例１，２のように、短軸の軸半径ｒ３が小さくなるに従い、軸半径ｒ３の評価誤差は大きくなり、Ｅｍｅｄ、Ｅｓｔｄも大きくなる傾向にある。これは、部分体積効果の影響と考えられる。つまり、ＣＴ画像では、対象物を画素の集合体として表現するため、軸半径が３ｖｏｘｅｌ以下の場合は、楕円体を表現する画素数が不足する。このことにより、断面積を精度良く計算するのが困難になり、評価誤差が大きくなったと考えられる。以上のことから、本実施例で示した評価手法では、軸半径が５ｖｏｘｅｌ以上の粒子を評価対象とすることが望ましい。
【表3】
【0103】
〈第２グループ～第４グループ〉
第２グループ～第４グループの評価結果のうち、Ｅｍｅｄ、Ｅｓｔｄが最も大きい結果を、それぞれ表４に示す。また、近似楕円体モデルの長軸と単位球の交点を計算し、長軸の走向角と傾斜角をシュミットネット下半球に投影した結果を、第２グループの実施例４及び実施例５は図１８（ａ）、第３グループの実施例６及び実施例７は図１８（ｂ）、第４グループの比較例３及び比較例４は図１８（ｃ）にそれぞれ示す。
第２グループ～第４グループの評価結果のうち、Ｅｍｅｄ、Ｅｓｔｄが最も大きい結果を、それぞれ表４に示す。また、近似楕円体モデルの長軸と単位球の交点を計算し、長軸の走向角と傾斜角をシュミットネット下半球に投影した結果を、第２グループの実施例４及び実施例５は図１８（ａ）、第３グループの実施例６及び実施例７は図１８（ｂ）、第４グループの比較例３及び比較例４は図１８（ｃ）にそれぞれ示す。
【表4】
【0104】
短軸の軸半径ｒ３が７ｖｏｘｅｌである第２グループにおいては、走向角の最大誤差が２．７度、傾斜角の最大誤差が２．３度となり、短軸の軸半径ｒ３が５ｖｏｘｅｌである第３グループにおいては、走向角の最大誤差が３．１度、傾斜角の最大誤差が３．８度となり、短軸の軸半径ｒ３が３ｖｏｘｅｌである第３グループにおいては、走向角の最大誤差が３．５度、傾斜角の最大誤差が４．６度となり、短軸が小さくなるほど評価誤差が大きくなる傾向にある。また、Ｅｍｅｄ、Ｅｓｔｄについても同様に、軸半径が小さくなるほど大きくなっている。このことから、第１グループと同様に、本実施例で示した評価手法は、軸半径が５ｖｏｘｅｌ以上の粒子を評価対象とすることが望ましい。
【0105】
〈第５グループ、第６グループ〉
第５グループ及び第６グループは、粒子の形状が楕円体ではなく矩形の粒子である。第５グループの実施例８と第６グループの実施例９の評価結果を表５に示す。また、実施例８と近似楕円体モデルを３次元表示した結果を図１９に、Ｅｉの計算結果を図２０に示す。
第５グループ及び第６グループは、粒子の形状が楕円体ではなく矩形の粒子である。第５グループの実施例８と第６グループの実施例９の評価結果を表５に示す。また、実施例８と近似楕円体モデルを３次元表示した結果を図１９に、Ｅｉの計算結果を図２０に示す。
【表5】
【0106】
実施例８は、矩形粒子の辺長さが、５、１０、１５ｖｏｘｅｌであったのに対し、解析の結果、軸半径は６．５、１２．８、１９．１ｖｏｘｅｌといずれも大きく評価された。また、Ｅｍｅｄは０．０８程度、Ｅｓｔｄは０．０５程度と算定され、第１グループ～第４グループに比べるとやや大きい結果となった。但し、粒子の形状指標として代表的なアスペクト比（ｒ３／ｒ１）、細長比（ｒ２／ｒ１）、扁平度（ｒ３／ｒ２）を比較すると、矩形粒子のアスペクト比は０．３３、細長比は０．６７、扁平度は０．５０であったのに対し、近似楕円体モデルのアスペクト比は０．３４、細長比は０．６７、扁平度は０．５１となり、ほぼ等しい値となった。
【0107】
また、走向角と傾斜角については、矩形粒子の各辺にほぼ平行な方向に楕円体の長軸、中間軸、短軸が近似されており、楕円体の姿勢が良好に評価されている。粒子の大きさを直接測定する方法の場合、実施例８のような矩形形状の粒子については、対角線の方向が最も長くなるため、長軸の向きも対角線の方向に評価される可能性がある。このことから、本実施例で示す評価手法では、粒子の形状が楕円体とは異なる場合であっても、走向角と傾斜角を適切に評価できることが確認された。また、実施例８の矩形粒子をＺ軸を中心に４５度回転した実施例９についても、走向角と傾斜角を適切に評価できたことが確認された。
【0108】
〈第７グループ〉
第７グループの比較例５は、長軸、中間軸、短軸がそれぞれ１５、１１、５ｖｏｘｅｌの楕円体２つが結合することによって、１つの粒子を構成しているＣＴ画像である。比較例５の評価結果を表６に示し、比較例５と近似楕円体モデルを３次元表示した結果を図２１に、Ｅｉの計算結果を図２２に示す。２つの楕円体粒子を１つの楕円体粒子として評価しており、中間軸として設定した軸が長軸として、長軸として設定した軸が中間軸として評価された。
第７グループの比較例５は、長軸、中間軸、短軸がそれぞれ１５、１１、５ｖｏｘｅｌの楕円体２つが結合することによって、１つの粒子を構成しているＣＴ画像である。比較例５の評価結果を表６に示し、比較例５と近似楕円体モデルを３次元表示した結果を図２１に、Ｅｉの計算結果を図２２に示す。２つの楕円体粒子を１つの楕円体粒子として評価しており、中間軸として設定した軸が長軸として、長軸として設定した軸が中間軸として評価された。
【表6】
【0109】
Ｅｍｅｄは０．０５５と算定され、短軸の軸半径ｒ３が５ｖｏｘｅｌ以上の第１グループ～第４グループと比較するとやや大きいが、粒子形状が矩形である第５グループや第６グループと比べると小さい値となった。但し、Ｅｓｔｄは０．１０７と第１グループ～第６グループと比較しても大きい値となった。これは、図２３に示すように、複数の粒子が結合して１つの粒子として識別された場合、例えば、粒子の中央のくびれの部分では粒子の断面積と近似楕円体モデルの断面積の差が大きくなる断面が生じることにより、Ｅｍｅｄはそれほど大きくならないが、Ｅｓｔｄが大きくなったと考えられる。つまり、複数の粒子を１つの粒子として、あるいは１つの粒子を複数の粒子として誤って識別した粒子については、Ｅｍｅｄだけでなく、Ｅｓｔｄを指標とすることにより、区別できると考えられる。一例としては、Ｅｓｔｄが０．１を超えた場合に、その粒子は不自然に抽出されたものとして識別し、排除するようにできる。
【0110】
〈砂粒子の形状評価〉
次に、砂粒子の形状評価を行った。粒径が３００～６００μｍを主体とする砂をアクリル容器に充填した試料のＣＴ画像を対象に、上述した評価方法を用い、砂粒子の形状解析を実施した。砂は、内径１０ｍｍの円筒形状のアクリル容器に自然落下法により充填した。図２４に、使用した砂の粒径分布を示す。ＣＴ撮影には、本願出願人（（一財）電力中央研究所）所有のマイクロフォーカスＸ線ＣＴ ＴＸＳ－ＣＴ ４５０／１６０を用いた。ＴＸＳ－ＣＴ ４５０／１６０には特性の異なる２つのＸ線源が搭載されているが、今回は、最大管電圧が１６０ｋＶの線源（最小焦点寸法４μｍ）を使用した。ＣＴ撮影では、管電圧を７５ｋＶ、管電流を１５０μＡとし、露光時間は５００ｍｓｅｃ、プロジェクションは３６００とし、ＳＩＤは１０００ｍｍ、ＳＯＤは７５ｍｍとして行った。この時、画素サイズは１５μｍ／ｐｉｘｅｌとした。また、再構成範囲は１５００×１５００×８００ｖｏｘｅｌである。図２５（ａ）及び（ｂ）に、水平断面及び鉛直断面のＣＴ画像の模式図を示す。
次に、砂粒子の形状評価を行った。粒径が３００～６００μｍを主体とする砂をアクリル容器に充填した試料のＣＴ画像を対象に、上述した評価方法を用い、砂粒子の形状解析を実施した。砂は、内径１０ｍｍの円筒形状のアクリル容器に自然落下法により充填した。図２４に、使用した砂の粒径分布を示す。ＣＴ撮影には、本願出願人（（一財）電力中央研究所）所有のマイクロフォーカスＸ線ＣＴ ＴＸＳ－ＣＴ ４５０／１６０を用いた。ＴＸＳ－ＣＴ ４５０／１６０には特性の異なる２つのＸ線源が搭載されているが、今回は、最大管電圧が１６０ｋＶの線源（最小焦点寸法４μｍ）を使用した。ＣＴ撮影では、管電圧を７５ｋＶ、管電流を１５０μＡとし、露光時間は５００ｍｓｅｃ、プロジェクションは３６００とし、ＳＩＤは１０００ｍｍ、ＳＯＤは７５ｍｍとして行った。この時、画素サイズは１５μｍ／ｐｉｘｅｌとした。また、再構成範囲は１５００×１５００×８００ｖｏｘｅｌである。図２５（ａ）及び（ｂ）に、水平断面及び鉛直断面のＣＴ画像の模式図を示す。
【0111】
〈粒子の識別〉
評価解析を行う前に、画像解析ソフトウェアを使って粒子の識別を行った。画像解析ソフトウェアは、ＶＧＳｔｕｄｉｏＭａｘ３．２を用い、オプションモジュールとして欠陥介在物解析モジュールを使用した。粒子の識別は、以下の手順で行った。
（１）画像のノイズ低減のため、ＣＴ画像に対しフィルタ処理を実施する。
（２）領域成長ツールを使って、個々の粒子について関心領域を１０個程度作成する。
（３）（２）で作成された関心領域について、グレイバリュー解析を行う。
（４）（３）の結果を基に、グレイバリューの範囲を設定し、画像全体で等値面を定義し、関心領域を作成する。
尚、今回は、最小グレイバリューを１０８００、最大グレイバリューを１２６８０とした。
評価解析を行う前に、画像解析ソフトウェアを使って粒子の識別を行った。画像解析ソフトウェアは、ＶＧＳｔｕｄｉｏＭａｘ３．２を用い、オプションモジュールとして欠陥介在物解析モジュールを使用した。粒子の識別は、以下の手順で行った。
（１）画像のノイズ低減のため、ＣＴ画像に対しフィルタ処理を実施する。
（２）領域成長ツールを使って、個々の粒子について関心領域を１０個程度作成する。
（３）（２）で作成された関心領域について、グレイバリュー解析を行う。
（４）（３）の結果を基に、グレイバリューの範囲を設定し、画像全体で等値面を定義し、関心領域を作成する。
尚、今回は、最小グレイバリューを１０８００、最大グレイバリューを１２６８０とした。
【0112】
手順（１）～（４）を実施した結果、画像全体で計１００３５個の粒子が識別された。尚、全ての粒子を識別するのに要した時間は１０分程度であった。図２５（ａ）及び（ｂ）と同じ断面について、識別された粒子に相当する画素を白色で塗りつぶした結果の模式図を図２５（ｃ）及び（ｄ）に示し、画像全体のグレイバリューを解析した結果を図２６に示す。図２５（ｃ）及び（ｄ）から分かるように、手順（４）で指定した範囲を外れるグレイバリューを有する画素は識別されていない。また、目視でも複数の粒子が結合している様子が散見される。これらについては、関心領域を新たに作成し、既に作成された関心領域と結合することや、収縮処理を行って結合を切り離すことも可能であるが、ここでは敢えて、それらの処理を行わないで評価解析を実施した。
【0113】
〈評価解析〉
グレイバリューの範囲を指定することによって識別された１００３５個の粒子のうち、Ｚ方向の層厚が６５０ｖｏｘｅｌ（０．０１５ｍｍ／ｐｉｘｅｌ×６５０ｖｏｘｅｌ＝９．７５ｍｍ）に位置する粒子を解析対象とした。尚、前章での検討において、軸半径が５ｖｏｘｅｌより小さい粒子については楕円体の近似精度が低下したことから、ここでは、球近似をした場合の直径が１０ｖｏｘｅｌ（０．０１５ｍｍ／ｐｉｘｅｌ×１０ｖｏｘｅｌ＝０．１５ｍｍ）以上である粒子、計４４００個について評価解析を実施することとした。
グレイバリューの範囲を指定することによって識別された１００３５個の粒子のうち、Ｚ方向の層厚が６５０ｖｏｘｅｌ（０．０１５ｍｍ／ｐｉｘｅｌ×６５０ｖｏｘｅｌ＝９．７５ｍｍ）に位置する粒子を解析対象とした。尚、前章での検討において、軸半径が５ｖｏｘｅｌより小さい粒子については楕円体の近似精度が低下したことから、ここでは、球近似をした場合の直径が１０ｖｏｘｅｌ（０．０１５ｍｍ／ｐｉｘｅｌ×１０ｖｏｘｅｌ＝０．１５ｍｍ）以上である粒子、計４４００個について評価解析を実施することとした。
【0114】
評価解析の結果、４４００個の粒子のうち、３個の粒子については、数式３の右辺分母の平方根内が負となることにより計算が進行せず、解を得ることができなかった。評価解析が収束し、解を得られた粒子４３９７個について、Ｅｍｅｄのヒストグラムを作成した結果を図２７に示す。解が得られた粒子４３９７個のうち、Ｅｍｅｄが０．０８以下の粒子は４０４１個（９１％）、Ｅｍｅｄが０．０８より大きく、０．１０以下の粒子は１７７個（４％）、Ｅｍｅｄが０．１０以上となった粒子は１７９個（４％）であった。また、Ｅｍｅｄが０．０８以下の粒子４０４１個のうちＥｓｔｄが０．１０以下となる粒子は４０１６個であり、Ｅｍｅｄが０．０８より大きく０．１０以下となる粒子１７７個のうちＥｓｔｄが０．１０以下となる粒子は１４１個、Ｅｍｅｄが０．１０より大きい粒子１７９個のうちＥｓｔｄが０．１０以下となる粒子は３３個となった。
【0115】
〈考察〉
解析対象とした４４００個の粒子のうち、９個の粒子についてＥｍｅｄ及びＥｓｔｄの計算結果を表７に示す。また、図２８（ａ）～（ｅ）及び図２９（ａ）～（ｄ）には、評価対象とした粒子のＣＴ画像と、評価解析で得た近似楕円体モデルのＣＴ画像を示す。
解析対象とした４４００個の粒子のうち、９個の粒子についてＥｍｅｄ及びＥｓｔｄの計算結果を表７に示す。また、図２８（ａ）～（ｅ）及び図２９（ａ）～（ｄ）には、評価対象とした粒子のＣＴ画像と、評価解析で得た近似楕円体モデルのＣＴ画像を示す。
【表7】
【0116】
Ｐ５６８８は、評価解析において数式３の右辺分母の平方根内が負になることにより、解を得ることができなかった粒子である。Ｐ５６８８は、図２８（ａ）に示すように、２つの粒子が結合し１つの粒子として識別されている粒子である。つまり、楕円体から大きく逸脱した形状であることにより、反復計算の過程で解が収束しなかったと考えられる。尚、ここでの評価解析では、Ｐ５６８８以外にも、２個の粒子が数式３の右辺分母の平方根内が負となることにより解を得ることができなかったが、いずれも複数の粒子が結合した粒子であった。
【0117】
Ｐ２８６３は、解が得られた粒子４３９７個のうち、Ｅｍｅｄが最も小さい粒子である。Ｐ２８６３は、図２８（ｂ）に示すように、単一の粒子が抽出されている。また、Ｅｍｅｄ及びＥｓｔｄが十分に小さく、楕円体への近似も良好にされている。一方で、Ｐ７８１１はＥｍｅｄ及びＥｓｔｄが０．２０３、０．２４３と大きい値となった。Ｐ７８１１は、図２８（ｃ）に示すように、複数の粒子が結合することで１つの粒子として識別された粒子であるが、この傾向はＥｍｅｄ及びＥｓｔｄにも現れている。
【0118】
Ｐ２３２１、Ｐ９０２８、Ｐ７３０４は、Ｅｍｅｄがそれぞれ０．０５９、０．０８５、０．１０３であり、Ｅｓｔｄはいずれも０．１以下となった粒子である。ＣＴ画像を確認した結果、いずれも単一の粒子が抽出されており、楕円体への近似も良好にされている。一方、Ｐ５４３０、Ｐ５９１５、Ｐ５９７６は、ＥｍｅｄがそれぞれＰ２３２１、 Ｐ９０２８、Ｐ７３０４とほぼ等しい粒子であるが、Ｅｓｔｄはいずれも０．１以上となった粒子である。ＣＴ画像を確認した結果、いずれも複数の粒子が結合することで１つの粒子として識別されていた。
【0119】
このことから、Ｅｍｅｄが一定程度小さい粒子であっても、Ｅｓｔｄを指標とすることにより、識別精度が悪い粒子を区別することができると考えられる。上述したように、Ｅｍｅｄが０．１０以下で計算が収束した粒子は４２１８個あったが、Ｅｓｔｄが０．１を上回る粒子は６１個存在しており、Ｅｓｔｄを指標としてこれらを取り除くことで、評価精度を向上することができた。
【0120】
評価解析によって解が得られた粒子のうち、Ｅｍｅｄ及びＥｓｔｄがいずれも０．１０より小さい粒子（４１５７個）を対象に、粒子の形状や配向を評価した結果を表８に示す。形状指標としては、近似楕円体モデルの長軸の軸半径ｒ１、中間軸の軸半径ｒ２、短軸の軸半径ｒ３、アスペクト比ｒ３／ｒ１、細長比ｒ２／ｒ１、扁平度ｒ３／ｒ２、Ｋｒｕｍｂｅｉｎ球形度Ｓｋ、長軸の傾斜角ｐｌ１、短軸の傾斜角ｐｌ３を個々の粒子について整理し、その統計解析結果（最小値、最大値、中央値、標準偏差）を示した。
【表8】
【0121】
また、図３０（ａ）及び（ｂ）、図３１（ａ）及び（ｂ）には、それぞれの形状指標についてヒストグラムを示すと共に、図３２には、扁平度と細長比の関係を整理したＺｉｎｇｇ Ｄｉａｇｒａｍのヒストグラムコンター図を示す。尚、図３１（ａ）のＫｒｕｍｂｅｉｎ球形度Ｓｋは、数式６により算出した。
【数6】
【0122】
図３０（ａ）に示すように、近似楕円体軸半径の中央値は、長軸が０．３０ｍｍ、中間軸が０．２０ｍｍ、短軸が０．１３ｍｍとなった。図２４に示したように、ふるい分けによって得られた砂の粒径加積曲線では、粒径が３００～６００μｍを主体としていたが、今回のＣＴ画像を使った評価においても、概ね同じ結果が得られた。図３０（ｂ）に示すように、アスペクト比ｒ３／ｒ１、細長比ｒ２／ｒ１、扁平度ｒ３／ｒ２の中央値は、それぞれ０．４４、０．６９、０．６７となった。非特許文献２における豊浦砂の形状解析結果では、アスペクト比、細長比、扁平度がそれぞれ０．５２，０．７３，０．７３と評価されており、細長比及び扁平度についてはほぼ同じ値となったが、本実施例で得られたアスペクト比についてはやや小さい結果となった。また、Ｋｒｕｍｂｅｉｎ球形度についても、アスペクト比と同様に、非特許文献２における豊浦砂の形状解析結果よりもやや小さい結果が得られた。
【0123】
扁平度ｒ３／ｒ２と細長比ｒ２／ｒ１の関係は、Ｚｉｎｇｇ Ｄｉａｇｒａｍと呼ばれる。これは、扁平度、細長比のいずれも０．６６を基準としてゾーニングを行い、粒子を球状（Ｓｐｈｅｒｅ）、板状（Ｄｉｓｋｓ）、棒状（Ｒｏｄｓ）、小判状（Ｂｌａｄｅｓ）に分類する方法である。今回の結果から扁平度、細長比をプロットし、図３２に示すように、ヒストグラムコンター図を作成すると、今回評価を行った粒子は、板状もしくは球状に分類される粒子が多かった。
【0124】
図３１（ｂ）に示すように、長軸、短軸の傾斜角の中央値は、それぞれ２３．０度、４０．２度となった。また、同図より、長軸の傾斜角は水平に近い角度にピークが集中しているのに対し、短軸の傾斜角は明瞭なピークが確認されなかった。今回撮影したサンプルは、砂を自然落下法によって充填していることから、長軸の傾斜角は、安定勾配である水平に近い角度で充填されたためと考えられる。
【0125】
上述したように、本実施例では、ＣＴ画像内の粒子に着目し、粒子の形状や配向を評価する方法について検討を行った。まず初めに、粒子の形状や配向を評価する方法として、抽出した粒子の断面積を複数組計算することで楕円体一般式の係数を求め、粒子を楕円体として近似する方法を構築した。次に、数値的に作成したＣＴ画像を用いてパラメータスタディーを実施し、評価を行う上での留意点について整理した。最後に、粒径分布が３００～６００μｍの砂を対象に撮影されたＣＴ画像について、砂粒子の形状や配向を評価した。その結果、以下のことが明らかになった。
【0126】
（１）評価対象とする粒子を球近似した場合の半径が５ｖｏｘｅｌより小さい場合、部分体積効果の影響によって評価精度が低下する可能性があることが明らかになった。このことから、本実施例で提案する手法を用いて粒子の形状や配向を評価する場合、球半径が５ｖｏｘｅｌ以上である粒子を解析対象とすることが好ましい。但し、条件が異なれば、必ずしも５ｖｏｘｅｌ以上であることには限られない。
【0127】
（２）評価対象とする粒子の形状が矩形である場合、近似楕円体モデルの軸半径は大きめに評価される。但し、アスペクト比、細長比、扁平率は評価対象とする粒子とほぼ等しい結果となる。また、長軸、中間軸、短軸の走向角や傾斜角は、矩形要素の各辺の方向にほぼ等しい結果が得られることが明らかになった。
【0128】
（３）複数の粒子を単一の粒子として識別してしまった場合、あるいは単一の粒子を複数の粒子として識別してしまった場合、評価解析の計算過程で解が得られない場合があることが明らかになった。また、解が得られた場合であっても、抽出した粒子の断面積と近似楕円体モデルの断面積の比Ｅｉ（数式５参照）を計算した時、識別精度が悪い粒子では両者の差が大きくなる断面が存在する。このことにより、Ｅｉの統計解析を行い、Ｅｉの中央値ＥｍｅｄやＥｉの標準偏差Ｅｓｔｄを計算することにより、識別精度が悪い粒子を区別できることが明らかになった。
【0129】
（４）複雑な形状である砂粒子であっても、本実施例で構築した方法によって、良好に楕円体として近似できることを確認した。また、楕円体への近似結果から、アスペクト比、細長比、扁平率、球形度、傾斜角を求め、本実施例で提案する方法によって得られた結果が、別手法で得られる結果と大差ない結果が得られることを確認した。
【0130】
以上の結果から、本発明によれば、岩石や金属などの集合体を構成する粒子のファブリックを、非破壊で撮影した断面画像を利用して簡易、かつ、高精度に解析できることが確認された。
【産業上の利用可能性】
【0131】
本発明に係る粒子のファブリックの解析方法，粒子のファブリックの解析装置，及び粒子のファブリックの解析プログラムは、岩石や金属などの集合体を構成する粒子のファブリックを、非破壊で撮影した断面画像を利用して簡易、かつ、高精度に解析することができるので、あくまで一例として挙げると、地盤や岩盤の力学特性と粒子形状に着目した研究、地質構造と粒子の配向に着目した研究、土質力学、堆積学などの分野で利用価値が高い。
【符号の説明】
【0132】
１０ コンピュータ／粒子のファブリックの解析装置
１１ 制御部
１１ａ 取得部
１１ｂ 算出部
１１ｃ 検証部
１１ｄ 解析部
１２ 記憶部
１３ 入力部
１４ 表示部
１５ メモリ
１７ 粒子のファブリックの解析プログラム
２０ データサーバ
１１ 制御部
１１ａ 取得部
１１ｂ 算出部
１１ｃ 検証部
１１ｄ 解析部
１２ 記憶部
１３ 入力部
１４ 表示部
１５ メモリ
１７ 粒子のファブリックの解析プログラム
２０ データサーバ
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】
【図25】
【図26】
【図27】
【図28】
【図29】
【図30】
【図31】
【図32】