JP7485076B2 - 推定装置、推定方法及び推定プログラム - Google Patents
推定装置、推定方法及び推定プログラム Download PDFInfo
- Publication number
- JP7485076B2 JP7485076B2 JP2022560628A JP2022560628A JP7485076B2 JP 7485076 B2 JP7485076 B2 JP 7485076B2 JP 2022560628 A JP2022560628 A JP 2022560628A JP 2022560628 A JP2022560628 A JP 2022560628A JP 7485076 B2 JP7485076 B2 JP 7485076B2
- Authority
- JP
- Japan
- Prior art keywords
- value
- dimensional space
- layer
- error
- volume
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims description 44
- 239000000463 material Substances 0.000 claims description 125
- 230000006870 function Effects 0.000 claims description 96
- 239000002245 particle Substances 0.000 claims description 96
- 238000004364 calculation method Methods 0.000 claims description 85
- 238000013527 convolutional neural network Methods 0.000 claims description 58
- 238000009434 installation Methods 0.000 claims description 42
- 238000005259 measurement Methods 0.000 claims description 38
- 230000008569 process Effects 0.000 claims description 27
- 238000010586 diagram Methods 0.000 description 29
- 230000004907 flux Effects 0.000 description 24
- 101150099612 Esrrb gene Proteins 0.000 description 11
- 230000015654 memory Effects 0.000 description 10
- 238000003325 tomography Methods 0.000 description 9
- 238000010801 machine learning Methods 0.000 description 8
- 238000012545 processing Methods 0.000 description 6
- 230000004044 response Effects 0.000 description 6
- 238000004590 computer program Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 230000010365 information processing Effects 0.000 description 4
- 230000009466 transformation Effects 0.000 description 4
- 238000004891 communication Methods 0.000 description 3
- 101100521334 Mus musculus Prom1 gene Proteins 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- JJWKPURADFRFRB-UHFFFAOYSA-N carbonyl sulfide Chemical compound O=C=S JJWKPURADFRFRB-UHFFFAOYSA-N 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 210000002569 neuron Anatomy 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000013307 optical fiber Substances 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V5/00—Prospecting or detecting by the use of ionising radiation, e.g. of natural or induced radioactivity
Landscapes
- Physics & Mathematics (AREA)
- High Energy & Nuclear Physics (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Image Analysis (AREA)
Description
本発明は、推定装置、推定方法及び非一時的なコンピュータ可読媒体に関し、特に、探査対象物を含む三次元空間の密度情報の分布を推定するための推定装置、推定方法及び非一時的なコンピュータ可読媒体に関する。
地下構造物（廃坑など）の探査にミュオグラフィを用いる際に、探査対象の構造物の形状をトモグラフィ的に三次元に復元するニーズが高まっている。例えば、特許文献１には、ミュー粒子を利用した三次元地盤探査システムに関する技術が開示されている。
ここで、ミューオンセンサにより測定されたミューオン粒子に関する測定データのみからトモグラフィを行う場合、推定すべきパラメータ数に対してデータが不足している（ill-posed）問題がある。そのため、ill-posedなトモグラフィ問題を解くことができず、トモグラフィ結果の精度が不十分であるという課題があった。尚、特許文献１にかかる技術には、測定されたミューオン粒子の情報から三次元に復元するためのモデルの具体的な構成やモデルの具体的な修正方法について開示されておらず、上述した課題を解決することができない。
本開示は、このような問題点を解決するためになされたものであり、探査対象物を含む三次元空間の密度情報の分布を精度良く推定するための推定装置、推定方法及び非一時的なコンピュータ可読媒体を提供することを目的とする。
本開示の第１の態様にかかる推定装置は、
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付手段と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定手段と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成手段と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出手段と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習手段と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算手段と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定手段と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力手段と、
を備える。
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付手段と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定手段と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成手段と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出手段と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習手段と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算手段と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定手段と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力手段と、
を備える。
本開示の第２の態様にかかる推定方法は、
コンピュータが、
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付け、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定し、
前記畳み込み演算により前記最終層の値を取得し、
前記取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成し、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出し、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習し、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得し、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定し、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する。
コンピュータが、
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付け、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定し、
前記畳み込み演算により前記最終層の値を取得し、
前記取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成し、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出し、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習し、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得し、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定し、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する。
本開示の第３の態様にかかる推定プログラムが格納された非一時的なコンピュータ可読媒体は、
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付処理と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定処理と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成処理と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出処理と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習処理と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算処理と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定処理と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力処理と、
をコンピュータに実行させる。
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付処理と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定処理と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成処理と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出処理と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習処理と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算処理と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定処理と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力処理と、
をコンピュータに実行させる。
本開示により、探査対象物を含む三次元空間の密度情報の分布を精度良く推定するための推定装置、推定方法及び非一時的なコンピュータ可読媒体を提供することができる。
以下では、本開示の実施形態について、図面を参照しながら詳細に説明する。各図面において、同一又は対応する要素には同一の符号が付されており、説明の明確化のため、必要に応じて重複説明は省略される。
＜実施形態１＞
図１は、本実施形態１にかかる推定装置１の構成を示すブロック図である。推定装置１は、探査対象物を含む三次元空間を通過した荷電粒子の測定結果と、三次元空間の事前知識とを用いて三次元空間における密度情報の分布の推定値を推定するための情報処理装置である。前提として、荷電粒子の複数の測定器が三次元空間の地下等の異なる場所に設置されているものとする。そして、各測定器は、探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する。尚、荷電粒子とは、例えば、ミューオン粒子であるが、これに限定されない。三次元空間とは、地上及び地下を含んでよい。そのため、三次元空間は、空気層や複数の地層を含む異なる複数の物質層を含むものである。そして、三次元空間の事前知識とは、三次元空間に含まれる各物質層の構成、順序、種別等である。また、探査対象物とは、地下構造物であり、例えば、廃坑等である。
図１は、本実施形態１にかかる推定装置１の構成を示すブロック図である。推定装置１は、探査対象物を含む三次元空間を通過した荷電粒子の測定結果と、三次元空間の事前知識とを用いて三次元空間における密度情報の分布の推定値を推定するための情報処理装置である。前提として、荷電粒子の複数の測定器が三次元空間の地下等の異なる場所に設置されているものとする。そして、各測定器は、探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する。尚、荷電粒子とは、例えば、ミューオン粒子であるが、これに限定されない。三次元空間とは、地上及び地下を含んでよい。そのため、三次元空間は、空気層や複数の地層を含む異なる複数の物質層を含むものである。そして、三次元空間の事前知識とは、三次元空間に含まれる各物質層の構成、順序、種別等である。また、探査対象物とは、地下構造物であり、例えば、廃坑等である。
推定装置１は、受付部１１と、設定部１２と、再構成部１３と、誤差算出部１４と、学習部１５と、演算部１６と、推定部１７と、出力部１８とを備える。受付部１１は、複数の測定器のそれぞれによる測定値と、各測定器の設置位置と、三次元空間に含まれる複数の物質層の種別との入力を受け付ける。尚、測定値は、三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向との測定結果そのものであるか、測定結果から公知技術により加工された画像データであってもよい。
設定部１２は、第１のボリューム内の各ボクセルに荷電粒子の密度情報として任意の値を設定する。ここで、ボリュームとは、三次元空間を複数のボクセルの集合で表現した情報であり、例えば、三次元配列である。そして、ボクセルは、三次元空間を構成する単位空間（例えば、立方体）である。また、第１のボリュームは、所定の畳み込みニューラルネットワーク（Convolutional Neural Network（ＣＮＮ））における第１層のニューロン（素子）とする。ここで、当該畳み込みニューラルネットワークは、第１のボリュームに対して所定の畳み込み演算を行うモデルである。また、設定部１２は、当該畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、測定器と物質層の組合せに対応させて設定する。つまり、設定部１２は、最終層におけるボリュームの数を測定器と物質層の組合せの数として設定する。そして、各第２のボリュームは、測定器と物質層の組合せに対応することとなる。例えば、測定器数が３、物質層数が５である場合、最終層のボリューム数は１５となる。そして、例えば、１番目のボリュームは、測定器１と物質層１に対応し、２番目のボリュームは、測定器１と物質層２に対応する。また、６番目のボリュームは、測定器２と物質層１に対応し、１０番目のボリュームは、測定器２と物質層５に対応する。以降同様に、１５番目のボリュームは、測定器３と物質層５に対応することとなる。
再構成部１３は、畳み込み演算により取得された最終層の値のうち、測定器ごとの各物質層に対応する第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する設置位置とを用いて、測定器ごとに測定値の期待値を再構成する。
誤差算出部１４は、測定器ごとの再構成された期待値と測定値との二乗誤差を含む誤差関数値を算出する。つまり、誤差関数は、測定器ごとの再構成された期待値と測定値との二乗誤差の項を含む。
学習部１５は、誤差関数値が所定値以上である場合に、誤差関数値を最小化するように、畳み込みニューラルネットワークのパラメータを学習する。
演算部１６は、学習後のパラメータを用いて畳み込み演算により最終層の値を取得する。
推定部１７は、取得された最終層の値と設置位置を用いて三次元空間における密度情報の分布の推定値を推定する。
出力部１８は、誤差関数値が所定値未満である場合に、推定値を出力する。尚、出力部１８は、推定値を三次元のトモグラフィデータとして出力してもよい。
図２は、本実施形態１にかかる推定方法の流れを示すフローチャートである。まず、受付部１１は、各測定器による測定値と、各測定器の設置位置と、三次元空間に含まれる複数の物質層の種別との入力を受け付ける（Ｓ１１）。
次に、設定部１２は、ＣＮＮの第１のボリューム内の各ボクセルに荷電粒子の密度情報として任意の値を設定する。また、設定部１２は、ＣＮＮの最終層の各第２のボリュームを、測定器と物質層の組合せに対応させて設定する（Ｓ１２）。尚、設定部１２は、外部からの指示に応じて、ＣＮＮの階層数、各階層のボリューム数、ボクセルの初期値や畳み込み演算に用いるパラメータ等を設定してもよい。
その後、演算部１６は、ＣＮＮの畳み込み演算により最終層の値を取得する（Ｓ１３）。そして、再構成部１３は、取得された最終層の値のうち、測定器ごとの各物質層に対応する第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、測定器に対応する設置位置とを用いて、測定器ごとに測定値の期待値を再構成する（Ｓ１４）。続いて、誤差算出部１４は、測定器ごとの再構成された期待値と測定値との二乗誤差を含む誤差関数値を算出する（Ｓ１５）。そして、誤差算出部１４は、誤差関数値が所定値以上であるか否かを判定する（Ｓ１６）。
誤差関数値が所定値以上である場合（ステップＳ１６でＹＥＳ）、学習部１５は、誤差関数値を最小化するように、ＣＮＮのパラメータを学習する（Ｓ１７）。そして、演算部１６は、学習後のパラメータを用いてＣＮＮの畳み込み演算により最終層の値を取得する（Ｓ１３）。以降、上記同様に、推定装置１は、ステップＳ１４からＳ１６を実行する。
誤差関数値が所定値未満である場合（ステップＳ１６でＮＯ）、推定部１７は、直近のステップＳ１３で取得された最終層の値と設置位置を用いて三次元空間における密度情報の分布の推定値を推定する（Ｓ１８）。そして、出力部１８は、推定値を出力する（Ｓ１９）。
このように本実施形態では、探査対象物を含む三次元領域の密度情報の分布を精度良く推定することができる。特に、荷電粒子の測定器と設置位置に加えて、観測対象の三次元領域の物質層の種別の情報を用いることで、上述したill-posed問題を解決できる。例えば、荷電粒子の測定器と観測対象の三次元領域の物質層との組わせを、ＣＮＮの最終層の各ボリュームに対応付けておくことで、期待値を再構成する際に、物質の固有密度係数を適切に用いることができる。そのため、期待値の再構成（算出）精度も向上する。それ故、ＣＮＮの機械学習の効率も向上する。よって、ＣＮＮの機械学習を安定的に収束させることができる。すなわち、解の安定性、収束性が向上する。また、ＣＮＮの最終層の各ボリュームの値から三次元領域の密度情報を推定する際に、測定器の設置位置を加味することで、推定精度を向上できる。
尚、推定装置１は、図示しない構成としてプロセッサ、メモリ及び記憶装置を備えるものである。また、当該記憶装置には、本実施形態にかかる推定方法の処理が実装されたコンピュータプログラムが記憶されている。そして、当該プロセッサは、記憶装置からコンピュータプログラムを前記メモリへ読み込ませ、当該コンピュータプログラムを実行する。これにより、前記プロセッサは、受付部１１、設定部１２、再構成部１３、誤差算出部１４、学習部１５、演算部１６、推定部１７及び出力部１８の機能を実現する。
または、受付部１１、設定部１２、再構成部１３、誤差算出部１４、学習部１５、演算部１６、推定部１７及び出力部１８は、それぞれが専用のハードウェアで実現されていてもよい。また、各装置の各構成要素の一部又は全部は、汎用または専用の回路（circuitry）、プロセッサ等やこれらの組合せによって実現されもよい。これらは、単一のチップによって構成されてもよいし、バスを介して接続される複数のチップによって構成されてもよい。各装置の各構成要素の一部又は全部は、上述した回路等とプログラムとの組合せによって実現されてもよい。また、プロセッサとして、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）、ＦＰＧＡ（field-programmable gate array）、量子プロセッサ（量子コンピュータ制御チップ）等を用いることができる。
また、推定装置１の各構成要素の一部又は全部が複数の情報処理装置や回路等により実現される場合には、複数の情報処理装置や回路等は、集中配置されてもよいし、分散配置されてもよい。例えば、情報処理装置や回路等は、クライアントサーバシステム、クラウドコンピューティングシステム等、各々が通信ネットワークを介して接続される形態として実現されてもよい。また、推定装置１の機能がＳａａＳ（Software as a Service）形式で提供されてもよい。
＜実施形態２＞
本実施形態２は、上述した実施形態１の具体的な一実施例である。まず、本実施形態２にかかる受付部は、探査対象物の形状に関する形状情報の入力をさらに受け付ける。そして、誤差算出部は、推定値のうち探査対象物の形状に関する推定値と、入力された形状情報との誤差をさらに含めて誤差関数値を算出する。つまり、誤差関数は、推定値のうち探査対象物の形状に関する推定値と、入力された形状情報との誤差の項を含む。このように、探査対象物の形状に関する事前知識を効果的に用いることで、ill-posedなトモグラフィ問題を解くことができ、トモグラフィ結果の精度をさらに向上させることができる。
本実施形態２は、上述した実施形態１の具体的な一実施例である。まず、本実施形態２にかかる受付部は、探査対象物の形状に関する形状情報の入力をさらに受け付ける。そして、誤差算出部は、推定値のうち探査対象物の形状に関する推定値と、入力された形状情報との誤差をさらに含めて誤差関数値を算出する。つまり、誤差関数は、推定値のうち探査対象物の形状に関する推定値と、入力された形状情報との誤差の項を含む。このように、探査対象物の形状に関する事前知識を効果的に用いることで、ill-posedなトモグラフィ問題を解くことができ、トモグラフィ結果の精度をさらに向上させることができる。
さらに、受付部は、探査対象物の三次元空間内の想定位置の入力をさらに受け付けるとよい。その場合、誤差算出部は、探査対象物の形状に関する推定値に対して、想定位置とのずれを調整するように誤差関数値を算出する。これにより、機械学習の効率、解の安定性、収束性がさらに向上する。そして、推定精度がさらに向上する。
ここで、本実施形態２では、地下の空洞（廃坑）や埋設物等の探査対象物の探査に宇宙線ミューオン粒子の測定データを用いるケースを説明する。そして、探査対象物の形状や埋設されている地層の位置等は、ある程度判明しているものの、正確な形状や埋設位置が不明であるものとする。例えば、探査対象物が地下１０ｍ程度に存在することがわかっているが、正確な位置が不明とする。そこで、地下２０ｍ程度の位置に複数の測定器を埋設するといったことである。また、探査対象物を含む地下空間の地層情報（物質層の構成、位置、種別）等は、別途、公知の物理探査により取得可能であるものとする。尚、本実施形態２は、他のケースにも適用可能である。
図３を用いて、本実施形態２にかかるミューオン粒子の測定の概念を説明する。三次元空間２は、複数の地層２１から２６を含み、地表上の空気層や岩盤層等も含む。地層２１等や空気層は、異なる複数の物質層の一例である。また、図３の例では、三次元空間２の地層２４に探査対象物２０が存在することが判明しているものとする。但し、探査対象物２０が地層２４に存在することが事前に判明していなくてもよい。また、地層２６には、複数の測定器２１１から２１ｎ（ｎは２以上の自然数。）が埋設されている。測定器２１１等は、ミューオン粒子のセンサである。測定器２１１等は、三次元空間２を通過してセンサ面に到達したミューオン粒子の数と入射方向とを測定する。尚、測定器２１１等は、公知技術を用いて実現可能である。
軌跡３１から３４は、上空から三次元空間２を通過するミューオン粒子の軌跡の一例である。ミューオン粒子の軌跡３１は、地上の空気層、地層２１から２３を通過し、地層２４内で探査対象物２０を通過し、地層２５及び２６を通過して測定器２１ｎに到達したことを示す。ミューオン粒子の軌跡３３は、同様に、探査対象物２０を通過した上で、測定器２１１に到達したことを示す。一方、ミューオン粒子の軌跡３２や３４は、探査対象物２０を通過せずに、測定器２１１や２１ｎに到達したことを示す。
図４は、本実施形態２にかかるミューオン粒子のエネルギーとミューオンフラックス（粒子数）との関係を説明するための図である。横軸は、ミューオン粒子のエネルギーであり、縦軸は、ミューオンフラックスである。図４に示すように、エネルギーの最小値Ｅｍｉｎ以下のミューオン粒子は、物質との相互作用で失われる。そのため、観測（測定）されるミューオン粒子数は、最小値Ｅｍｉｎ以上である。
ここで、ミューオン粒子のエネルギー損失式は、経験的に以下の式（１）により定義される。
・・・式（１）
ここで、ａ（Ｅ）とｂ（Ｅ）とは物質の種類により定まる。そして、物質層の長さｌを通過できるミューオン粒子のエネルギーの最小値Ｅｍｉｎは、以下の式（２）により求められる。
・・・式（２）
ここで、ａ（Ｅ）とｂ（Ｅ）とは物質の種類により定まる。そして、物質層の長さｌを通過できるミューオン粒子のエネルギーの最小値Ｅｍｉｎは、以下の式（２）により求められる。
そして、測定器２１１等による測定結果は、公知技術を用いてミューオンフラックス画像（測定値の一例）に変換が可能である。図５は、本実施形態２にかかる測定値から生成されるミューオンフラックス画像の例を示す図である。横軸は、測定器のセンサ面における水平方向の周囲１８０度のタンジェント値であり、縦軸は、天頂方向からの方位角のタンジェント値である。また、各画素は、ミューオン粒子数の密度の大小を濃淡で示している。
ここで、測定器fのセンサ面上の位置（座標）ｊにおいて測定されるミューオン粒子数の期待値Ｆｅｘｐ（ｆ、ｊ）は、以下の式（３）により求められる。
・・・式（３）
ここで、ｄ（ｊ）は、センサ面上の位置（座標）ｊにおいて測定されるミューオン粒子の入射方向である。Ωは、入射方向（立体角、方位角）の変数である。関数Ｎ（Ｅ、ｄ（ｊ））は、測定器ｆのセンサ面上の位置（座標）ｊにおいて測定されるミューオン粒子のエネルギーＥ及び入射方向ｄ（ｊ）におけるミューオン粒子数の関数である。以下の式（４）は、ミューオン粒子のエネルギーＥμ及び方位角θにおける関数Ｎの微分強度（differential intensity）の定義である。
・・・式（４）
ここで、ｐμは、ミューオンの運動量であり光速に近い速度を持つ場合はエネルギーＥμと等価になる。また、ｇ（ｐμ、θ）は、以下の式（５）で定義される。
・・・式（５）
ここで、ｄ（ｊ）は、センサ面上の位置（座標）ｊにおいて測定されるミューオン粒子の入射方向である。Ωは、入射方向（立体角、方位角）の変数である。関数Ｎ（Ｅ、ｄ（ｊ））は、測定器ｆのセンサ面上の位置（座標）ｊにおいて測定されるミューオン粒子のエネルギーＥ及び入射方向ｄ（ｊ）におけるミューオン粒子数の関数である。以下の式（４）は、ミューオン粒子のエネルギーＥμ及び方位角θにおける関数Ｎの微分強度（differential intensity）の定義である。
ここで、ｐμは、ミューオンの運動量であり光速に近い速度を持つ場合はエネルギーＥμと等価になる。また、ｇ（ｐμ、θ）は、以下の式（５）で定義される。
図６は、本実施形態２にかかる推定装置４００の構成を示すブロック図である。推定装置４００は、上述した推定装置１の機能を包含した上で、本実施形態２に特有の機能を備える。尚、推定装置４００は、複数台のコンピュータに冗長化されてもよく、各機能ブロックが複数台のコンピュータに分散して実現されてもよい。
推定装置４００は、記憶部４１０と、メモリ４２０と、ＩＦ（InterFace）部４３０と、制御部４４０とを備える。記憶部４１０は、ハードディスク、フラッシュメモリ等の記憶装置の一例である。記憶部４１０は、推定プログラム４１１と、ＣＮＮ４１２と、物質種別４１３と、固有密度係数４１４と、最終層定義情報４１５を記憶する。推定プログラム４１１は、本実施形態２にかかる推定方法の処理が実装されたコンピュータプログラムである。
ＣＮＮ４１２は、畳み込みニューラルネットワークの一例である。ＣＮＮ４１２は、第１層４１２１、・・・最終層４１２Ａ（Ａは、３以上の自然数。）とパラメータ４１２０を含む。第１層４１２１は、ボリュームｖｏ１（第１のボリュームの一例）を含む。ボリュームｖｏ１は、三次元空間２を複数のボクセルｖｘ１１からｖｘ１ｋ（ｋは、２以上の自然数。）で表現した三次元配列である。各ボクセルには、ミューオン粒子の密度情報を示す値が設定される。最終層４１２Ａは、ボリュームｖｏＡ１１からｖｏＡ１Ｎｍ、・・・ｖｏＡＮｆ１からｖｏＡＮｆＮｍを含む。これらは、複数の第２のボリュームの一例である。ここで、Ｎｍは、物質層の総数であり、Ｎｆは、探査に用いられた測定器の総数、言い換えると、入力されるミューオンフラックス画像の総数である。ボリュームｖｏＡ１１は、測定器１の物質層１に対応する。ボリュームｖｏＡ１Ｎｍは、測定器１の物質層Ｎｍに対応する。ボリュームｖｏＡＮｆ１は、測定器Ｎｆの物質層１に対応する。ボリュームｖｏＡＮｆＮｍは、測定器Ｎｆの物質層Ｎｍに対応する。パラメータ４１２０は、ＣＮＮ４１２の各層の間の変換を行う際に用いられる重み付け係数及びバイアス値の集合である。
図７は、本実施形態２にかかるＣＮＮ４１２の構成、期待値の再構成、密度情報の推定の概念を説明するための図である。ＣＮＮ４１２は、第１層、第２層、・・・第Ａ－１層、第Ａ層（最終層）で構成される畳み込みニューラルネットワークである。第１層には、ボリュームｖｏ１が属する。尚、第１層には２以上のボリュームが属しても良い。第２層には、ボリュームｖｏ２１、ｖｏ２２及びｖｏ２３が属する。尚、第２層には、２又は４以上のボリュームが属しても良い。ボリュームｖｏ１の各ボクセルの値は、重み付け係数とバイアス値を用いて変換され、ボリュームｖｏ２１の対応する各ボクセルに設定される。同様に、ボリュームｖｏ１の各ボクセルの値は、上記とは異なる重み付け係数とバイアス値を用いて変換され、ボリュームｖｏ２２及びｖｏ２３のそれぞれの対応する各ボクセルに設定される。尚、重み付け係数とバイアス値は、ＣＮＮ４１２のパラメータ４１２０である。同様に、第２層の各ボリュームは、第３層の各ボリュームにそれぞれのパラメータを用いて変換され、第Ａ－１層まで変換される。ここでは、第２層から特定層までは、ボリュームの粒度が大きくなるように変換し、特定層から第Ａ－１層までは、ボリュームの粒度が小さくなるように変換する。言い換えると、第２層から特定層までは、ボクセルのスケールが大きくするように変換し、特定層から第Ａ－１層までは、ボクセルのスケールが小さくなるように変換する。また、第１層、第２層及び第Ａ－１層に属するボリュームのボクセル数は同じとするが、これに限定されない。第Ａ－１層の各ボリュームは、第Ａ層の各ボリュームへ変換される。例えば、ボリュームｖｏＡ－１１の各ボクセルの値は、異なるパラメータを用いて変換され、ボリュームｖｏＡ１１からｖｏＡ１Ｎｍ、・・・ｖｏＡＮｆ１からｖｏＡＮｆＮｍの各ボクセルに設定される。
図６に戻り説明を続ける。物質種別４１３は、物質層の物質の種別を示す情報である。固有密度係数４１４は、物質の固有密度係数である。物質種別４１３と固有密度係数４１４とは対応付けられている。
最終層定義情報４１５は、最終層４１２Ａに含まれる各ボリュームのボリュームＩＤ４１５１に対応付けられた測定器ＩＤ４１５２及び物質種別４１５３の組合せを定義した情報である。
メモリ４２０は、ＲＡＭ（Random Access Memory）等の揮発性記憶装置であり、制御部４４０の動作時に一時的に情報を保持するための記憶領域である。ＩＦ部４３０は、
推定装置４００の外部との入出力を行うインタフェースである。例えば、ＩＦ部４３０は、キーボード、マウス、タッチパネル等の入力デバイス（不図示）を介して、ユーザの操作を受け付け、受け付けた操作内容を制御部４４０へ出力する。また、ＩＦ部４３０は、制御部４４０からの指示に応じて、タッチパネル、表示装置、プリンタ等（不図示）へ出力を行う。
推定装置４００の外部との入出力を行うインタフェースである。例えば、ＩＦ部４３０は、キーボード、マウス、タッチパネル等の入力デバイス（不図示）を介して、ユーザの操作を受け付け、受け付けた操作内容を制御部４４０へ出力する。また、ＩＦ部４３０は、制御部４４０からの指示に応じて、タッチパネル、表示装置、プリンタ等（不図示）へ出力を行う。
制御部４４０は、推定装置４００の各構成を制御するプロセッサつまり制御装置である。制御部４４０は、記憶部４１０から推定プログラム４１１をメモリ４２０へ読み込ませ、推定プログラム４１１を実行する。これにより、制御部４４０は、受付部４４１、設定部４４２、再構成部４４３、誤差算出部４４４、学習部４４５、演算部４４６、推定部４４７及び出力部４４８の機能を実現する。
受付部４４１は、上述した受付部１１の一例である。受付部４４１は、複数のミューオンフラックス画像と、測定器２１１から２１ｎの設置位置と、三次元空間２に関する地層情報と、探査対象物２０の形状情報及び想定位置との入力を受け付ける。入力されるミューオンフラックス画像は、測定器２１１等のそれぞれの測定結果から変換された画像データである。つまり、各ミューオンフラックス画像は、測定器を一意に特定し、設置位置と一対一に対応する情報である。また、ミューオンフラックス画像は、上述した通り、測定器のセンサ面の各位置におけるミューオン粒子の測定数と、各ミューオン粒子の入射方向を含む情報である。設置位置は、三次元空間２における三次元の位置座標の集合である。つまり、設置位置は、対応する測定器のセンサ面上の位置ｊの三次元の位置座標の集合である。地層情報は、地表上の空気層や岩盤層、地下の岩盤層の構成（位置関係、層の厚さ、物質層の種別、含水量等）を含むものとする。形状情報は、探査対象物２０の想定される形状（縦、横、高さ等）を示す情報である。想定位置は、三次元空間２内で探査対象物２０が存在すると想定される三次元位置座標の集合である。
設定部４４２は、上述した設定部１２の一例である。設定部４４２は、上記入力に応じてＣＮＮ４１２の各種設定を行う。すなわち、設定部４４２は、ＣＮＮ４１２の第１層のボリュームｖｏ１の各ボクセルに初期値を設定する。尚、設定部４４２は、ユーザから入力された初期値を各ボクセルに設定してもよい。また、設定部４４２は、入力されたミューオンフラックス画像数と地層情報内の物質層の数とを乗算した数を算出し、ＣＮＮ４１２の最終層のボリューム数を、乗算した数として設定する。また、設定部４４２は、最終層４１２Ａに含まれる各ボリュームのボリュームＩＤ４１５１に、測定器ＩＤ４１５２及び物質種別４１５３の組合せを対応付けて最終層定義情報４１５として記憶部４１０に保存する。その他、設定部４４２は、ユーザの入力に応じて、ＣＮＮ４１２の階層数、各階層のボリューム数、ボクセルのサイズ等を設定しても良い。
また、設定部４４２は、上記入力に応じて後述する誤差関数の設定を行う。例えば、設定部４４２は、入力されたミューオンフラックス画像の総数（測定器数Ｎｆ）、物質層の総数Ｎｍ等を設定する。また、設定部４４２は、入力された地層情報に含まれる物質層の物質種別４１３に対応付けられた固有密度係数４１４を取得し、誤差関数に設定する。尚、設定部４４２は、入力された地質情報に含まれる含水量を加味して固有密度係数４１４を調整して誤差関数に設定するとよい。これにより、再構成の精度や密度情報の推定精度が向上する。また、設定部４４２は、入力された形状情報及び想定位置に基づく設定を誤差関数に行う。例えば、設定部４４２は、形状情報に対応する形状フィルタ（後述）を誤差関数に設定する。
演算部４４６は、上述した演算部１６の一例である。演算部４４６は、設定されたＣＮＮ４１２の畳み込み演算により最終層の値を取得する。また、演算部４４６は、学習部４４５による学習後のパラメータを用いてＣＮＮ４１２の畳み込み演算により最終層の値を取得する。
再構成部４４３は、上述した再構成部１３の一例である。再構成部４４３は、畳み込み演算により導出された最終層の各ボリュームのボクセル値を取得し、測定器単位に、ミューオンフラックス画像の期待値を再構成する。例えば、図７に示すように、再構成部４４３は、測定器１に対応するボリュームｖｏＡ１１からボリュームｖｏＡ１Ｎｍを対象に、同じ位置のボクセル値を用いて、ボクセルごとのミューオン粒子の密度情報を算出する。そして、再構成部４４３は、算出した密度情報を合成して再構成画像ｒ１を生成する。また、再構成部４４３は、測定器２に対応するボリュームｖｏＡ２１からボリュームｖｏＡ２Ｎｍを対象に、再構成画像ｒ２を生成する。以降同様に、再構成部４４３は、測定器Ｎｆに対応するボリュームｖｏＡＮｆ１からボリュームｖｏＡＮｆＮｍを対象に、再構成画像ｒＮｆを生成する。
推定部４４７は、上述した推定部１７の一例である。推定部４４７は、畳み込み演算により導出された最終層の各ボリュームのボクセル値を取得し、各ボクセル値と設置位置とを用いて三次元空間２における密度情報の分布の推定値を推定する。例えば、図７に示すように、推定部４４７は、ボリュームｖｏＡ１１からボリュームｖｏＡＮｆＮｍの各ボクセルの値と、入力された各測定器の設置位置とを用いて三次元復元密度情報ｒｄを推定値として推定する。
誤差算出部４４４は、上述した誤差算出部１４の一例である。誤差算出部４４４は、設定された誤差関数に、測定器ごとの再構成画像ｒ１からｒＮｆ及びミューオンフラックス画像ｏｂ１からｏｂＮｆを入力し、これらの二乗誤差を算出する。また、誤差算出部４４４は、設定された誤差関数に、推定値と入力された形状情報、想定位置等を入力し、密度情報の誤差を算出する。そして、誤差算出部４４４は、二乗誤差及び密度情報の誤差により誤差関数値を算出する。このとき、誤差算出部４４４は、推定値のうち探査対象物２０の形状に関する推定値と、入力された形状情報との誤差をさらに含めて誤差関数値を算出する。さらに、誤差算出部４４４は、探査対象物２０の形状に関する推定値に対して、想定位置とのずれを調整するように誤差関数値を算出する。また、誤差算出部４４４は、誤差関数値が所定値以上であるか否かを判定する。つまり、誤差算出部４４４は、機械学習の収束条件を満たすか否かを判定する。
学習部４４５は、上述した学習部１５の一例である。学習部４４５は、誤差関数値が所定値以上である場合、誤差関数値を最小化するように、ＣＮＮ４１２のパラメータ４１２０を学習する。特に、学習部４４５は、上記誤差関数の偏微分を用いた誤差逆伝搬法により第Ａ層（最終層４１２Ａ）から第Ａ－１層への逆伝搬を行い、第Ａ層と第Ａ－１層の間のパラメータを更新する。また、学習部４４５は、一般的な誤差逆伝搬法により第Ａ－１層から第１層４１２１への逆伝搬を行い、第Ａ－１層から第１層４１２１の間のパラメータを更新する。
出力部４４８は、上述した出力部１８の一例である。出力部４４８は、誤差関数値が所定値未満である場合、推定部４４７により直近に推定された推定値を出力する。
図８は、本実施形態２にかかる推定処理の流れを示すフローチャートである。まず、受付部４４１は、測定器ごとのミューオンフラックス画像と、各測定器の設置位置と、地層情報と、探査対象物２０の形状情報及び想定位置との入力を受け付ける（Ｓ２０１）。
次に、設定部４４２は、ＣＮＮ４１２の第１層４１２１のボリュームｖｏ１に初期値を設定し、最終層４１２Ａの各ボリュームを、測定器と物質層の組合せに対応させて設定する（Ｓ２０２）。このとき、上述したように、設定部４４２は、ユーザの入力に応じてＣＮＮ４１２の各種設定を行っても良い。また、設定部４４２は、ＣＮＮ４１２のパラメータ４１２０の任意の初期値を設定する。併せて、設定部４４２は、上述したように誤差関数を設定する。
その後、演算部４４６は、ＣＮＮ４１２の畳み込み演算を実行する（Ｓ２０３）。初回の畳み込み演算において、演算部４４６は、パラメータ４１２０の初期値を用いてＣＮＮ４１２の畳み込み演算を実行する。これにより、演算部４４６は、最終層４１２Ａの各ボリュームの各ボクセル値を取得する。
ステップＳ２０３の後、再構成部４４３は、測定器ごとの期待値の再構成（復元、算出）を行う（Ｓ２０４）。例えば、上述したように、再構成部４４３は、測定器１について再構成画像ｒ１を生成し、測定器２について再構成画像ｒ２を生成し、以下同様に、測定器Ｎｆについて再構成画像ｒＮｆを生成する。
再構成部４４３は、例えば、以下の式（６）により、測定器ｆに対応する各ボリューム内のｉ番目のボクセルｖｉにおける密度情報ρｆ（ｖｉ）を算出（再構成）する。
・・・式（６）
ここで、
は、物質種別ｍにおける固有密度係数である。αは任意の係数である。また、
は、測定器ｆ及び物質種別ｍに対応するボリューム内のｉ番目のボクセルの値である。そして、再構成部４４３は、以下の式（７）により、ボクセルｖｉ及び入射方向ｄ（ｊ）における最小エネルギーＥｍｉｎ（ｆ、ｄ（ｊ））を算出する。
・・・式（７）
ここで、ａｍ及びｂｍは、物質種別ｍにおけるエネルギー損失の係数である。また、
は、三次元空間２内のボクセルｖｉを入射方向ｄ（ｊ）で通過して測定器ｆに到達したミューオン粒子について、ボクセルｖｉ内を通過した軌跡の長さである。例えば、再構成部４４３は、測定器のセンサ面の位置座標と対象ボクセルの位置座標と入射方向から、画素ｊで測定されたミューオン粒子が対象ボクセルを通過したか否かを判定する。そして、通過したと判定した場合、再構成部４４３は、対象ボクセルに対するミューオン粒子の入射位置及び出射位置を特定し、入射位置及び出射位置の長さＬを算出する。但し、ボクセル内を通過したミューオン粒子の軌跡の長さの算出方法は、これに限定されない。
ここで、
ここで、ａｍ及びｂｍは、物質種別ｍにおけるエネルギー損失の係数である。また、
図９は、本実施形態２にかかる所定のボクセルｖｉを通過する複数のミューオン粒子の軌跡と、各ミューオン粒子を測定する複数の測定器との関係を説明するための図である。ここでは、図９の左下にある測定されたミューオンフラックス画像ｆ（測定器ｆに対応）内の画素ｊ（センサ面の位置）に入射方向ｄ（ｊ）で到達したミューオン粒子について、ボクセルｖｉを通過した軌跡の長さが
であることを示す。同様に、図９の右下にある測定されたミューオンフラックス画像ｆ’（測定器ｆ’に対応）内の画素ｊ’に入射方向ｄ（ｊ’）で到達したミューオン粒子について、ボクセルｖｉを通過した軌跡の長さが
であることを示す。
図１０は、本実施形態２にかかる所定のボクセルｖｉを通過する複数のミューオン粒子の軌跡と、各ミューオン粒子を測定する単一の測定器との関係を説明するための図である。ここでは、図１０の下にある測定されたミューオンフラックス画像ｆ内の画素集合
に様々な入射方向で到達した複数のミューオン粒子について、ボクセルｖｉを通過した軌跡を示す。各軌跡の長さは、上述したように算出可能である。
そして、再構成部４４３は、式（７）により算出された最小エネルギーＥｍｉｎ（ｆ、ｄ（ｊ））を式（３）に代入して、期待値Ｆｅｘｐ（ｆ、ｊ）を算出する。尚、期待値Ｆｅｘｐ（ｆ、ｊ）は、測定器fのセンサ面上の位置（座標）ｊにおいて測定されるミューオン粒子数の密度情報ともいえる。
また、ステップＳ２０３の後、推定部４４７は、三次元空間２の密度情報の分布の推定値を推定する（Ｓ２０５）。そのため、推定値には、探査対象物２０の形状に関する推定値も含まれる。ここで、推定部４４７は、以下の式（８）により、三次元空間２に対応する立体空間内のボクセルｖｉにおける密度情報ρ（ｖｉ）を算出（推定）する。
・・・式（８）
そして、推定部４４７は、全ボクセルの密度情報ρ（ｖｉ）をまとめて三次元空間２の密度情報の分布の推定値とする。
ステップＳ２０４及びＳ２０５の後、誤差算出部４４４は、誤差関数値を算出する（Ｓ２０６）。例えば、誤差算出部４４４は、以下の式（９）の誤差関数により誤差関数値Ｅｒｒを算出する。
Ｅｒｒ＝ｃ１Ｅｒｒ１＋ｃ２Ｅｒｒ２ ・・・式（９）
尚、ｃ１及びｃ２は、任意の係数である。
Ｅｒｒ＝ｃ１Ｅｒｒ１＋ｃ２Ｅｒｒ２ ・・・式（９）
尚、ｃ１及びｃ２は、任意の係数である。
ここで、誤差Ｅｒｒ１は、ミューオンフラックス画像の期待値と測定値との誤差を示す。言い換えると、誤差Ｅｒｒ１は、各測定器において測定されたミューオン粒子の画素ごとの密度情報（測定値）と期待値との二乗誤差を、全画素及び全測定器について合計した値である。例えば、誤差算出部４４４は、以下の式（１０）により誤差Ｅｒｒ１を算出できる。
・・・式（１０）
ここで、Ｎｐは、測定器ｆのセンサ面の画素番号の最大値である。期待値Ｆｅｘｐ（ｆ、ｊ）は、ステップＳ２０４により算出された値である。測定値Ｆｏｂｓ（ｆ、ｊ）は、測定器ｆのセンサ面の画素ｊにおいて測定されたミューオン粒子数、つまり密度情報の測定値である。
ここで、Ｎｐは、測定器ｆのセンサ面の画素番号の最大値である。期待値Ｆｅｘｐ（ｆ、ｊ）は、ステップＳ２０４により算出された値である。測定値Ｆｏｂｓ（ｆ、ｊ）は、測定器ｆのセンサ面の画素ｊにおいて測定されたミューオン粒子数、つまり密度情報の測定値である。
また、式（９）の誤差Ｅｒｒ２は、推定値のうち探査対象物の形状に関する推定値と、入力された形状情報との誤差を含む。例えば、誤差算出部４４４は、以下の式（１１）により誤差Ｅｒｒ２を算出できる。
・・・式（１１）
ここで、ｉはボクセル番号、Ｎｖはボクセル番号の最大値である。ｒは探査対象物２０に当てはめる形状フィルタの水平方向の回転角、Ｑｈはこの水平方向探索回転角の最大値である。Ｘｒ（ｖｉ）は、以下の式（１２）で定義される。
・・・式（１２）
ここで、ｒ’は水平方向探索回転角度ｒに直交するように定めれられた探査対象物２０に当てはめる形状フィルタの鉛直方向の回転角、Ｑｖはこの鉛直方向探索回転角の最大値である。ｒ⊥はｒ及びｒ’と直交する軸である。図１１は、本実施形態２にかかるＥｒｒ２を算出する際の探査対象物の軸の回転方向を説明するための図である。図１２は、本実施形態２にかかるＥｒｒ２を算出する際の形状・密度勾配フィルタの回転方向を説明するための図である。
ここで、ｉはボクセル番号、Ｎｖはボクセル番号の最大値である。ｒは探査対象物２０に当てはめる形状フィルタの水平方向の回転角、Ｑｈはこの水平方向探索回転角の最大値である。Ｘｒ（ｖｉ）は、以下の式（１２）で定義される。
ここで、ｒ’は水平方向探索回転角度ｒに直交するように定めれられた探査対象物２０に当てはめる形状フィルタの鉛直方向の回転角、Ｑｖはこの鉛直方向探索回転角の最大値である。ｒ⊥はｒ及びｒ’と直交する軸である。図１１は、本実施形態２にかかるＥｒｒ２を算出する際の探査対象物の軸の回転方向を説明するための図である。図１２は、本実施形態２にかかるＥｒｒ２を算出する際の形状・密度勾配フィルタの回転方向を説明するための図である。
そして、式（１２）の３つのＺは、以下の式（１３）、式（１４）及び式（１５）で定義される。
・・・式（１３）
・・・式（１４）
・・・式（１５）
ここで、ｗ１、ｗ２及びｗ３は、任意の係数である。ｊは、図１０に示した画素集合に含まれる画素の位置を示す変数である。関数Ｆは、探査対象物２０の形状の推定値の境界と、入力された形状データの境界とのずれを調整するための関数である。そして、式（１３）の
及び式（１４）の
は、まとめて
と表記し、以下の式（１６）で定義される。
・・・式（１６）
ここで、μｍは物質種別ｍにおける密度関数（密度モデル）である。尚、μｍは物質種別４１３に対応付けられて予め保存されているものとする。そして、Ｖは以下の式（１７）で定義される。
・・・式（１７）
ここで、Ｍは形状フィルタの要素番号（ボクセル数）の最大値（フィルタのサイズ）である。ｃｒは、回転角度ｒでの座標変換関数である。そして、Δｖは以下の式（１８）で定義される。
・・・式（１８）
ここで、Ｍ’は密度勾配（算出）フィルタの要素番号（ボクセル数）の最大値（フィルタのサイズ）である。そして、ｇｒ／ｒｒ‘（ｖｉ）は、以下の式（１９）で定義される。
・・・式（１９）
ここで、ｊは、密度勾配（算出）フィルタにおける要素番号である。Ｇｒ／ｒｒ‘（ｊ）は、要素番号ｊにおける密度勾配（算出）フィルタ関数である。尚、Ｍ、Ｍ’及びＧといった形状フィルタ、密度勾配（算出）フィルタに関する情報は、入力される形状情報に含まれる。また、ρＲは、密度の正規化係数である。図１３は、本実施形態２にかかるＥｒｒ２を算出する際の形状・密度勾配フィルタと座標変換を説明するための図である。なお形状フィルタと密度勾配フィルタの両中心は一致させた上でフィルタ処理を行う。
ここで、μｍは物質種別ｍにおける密度関数（密度モデル）である。尚、μｍは物質種別４１３に対応付けられて予め保存されているものとする。そして、Ｖは以下の式（１７）で定義される。
ここで、Ｍは形状フィルタの要素番号（ボクセル数）の最大値（フィルタのサイズ）である。ｃｒは、回転角度ｒでの座標変換関数である。そして、Δｖは以下の式（１８）で定義される。
ここで、Ｍ’は密度勾配（算出）フィルタの要素番号（ボクセル数）の最大値（フィルタのサイズ）である。そして、ｇｒ／ｒｒ‘（ｖｉ）は、以下の式（１９）で定義される。
ここで、ｊは、密度勾配（算出）フィルタにおける要素番号である。Ｇｒ／ｒｒ‘（ｊ）は、要素番号ｊにおける密度勾配（算出）フィルタ関数である。尚、Ｍ、Ｍ’及びＧといった形状フィルタ、密度勾配（算出）フィルタに関する情報は、入力される形状情報に含まれる。また、ρＲは、密度の正規化係数である。図１３は、本実施形態２にかかるＥｒｒ２を算出する際の形状・密度勾配フィルタと座標変換を説明するための図である。なお形状フィルタと密度勾配フィルタの両中心は一致させた上でフィルタ処理を行う。
図８に戻り説明を続ける。ステップＳ２０６の後、誤差算出部４４４は、算出した誤差関数値が所定値以上であるか否かを判定する（Ｓ２０７）。誤差関数値が所定値以上である場合（ステップＳ２０７でＹＥＳ）、学習部４４５は、ＣＮＮ４１２のパラメータを学習する（Ｓ２０８）。
まず、学習部４４５は、最終層である第Ａ層から第Ａ－１層へ誤差を逆伝搬して第Ａ層と第Ａ－１層の間のパラメータを更新する。具体的には、まず、誤差Ｅｒｒ１の誤差逆伝搬量の計算について説明する。学習部４４５は、以下の式（２０）により、誤差Ｅｒｒ１の偏微分値を算出する。
・・・式（２０）
ここで、最小エネルギーＥｍｉｎ（ｆ、ｄ（ｊ））の偏微分値は以下の式（２１）により算出できる。
・・・式（２１）
ここで、ｍ’は物質種別の変数である。ｆ’はミューオンフラックス画像（又は測定器）を識別するための変数である。
ここで、最小エネルギーＥｍｉｎ（ｆ、ｄ（ｊ））の偏微分値は以下の式（２１）により算出できる。
ここで、ｍ’は物質種別の変数である。ｆ’はミューオンフラックス画像（又は測定器）を識別するための変数である。
次に、誤差Ｅｒｒ２の誤差逆伝搬量の計算について説明する。学習部４４５は、以下の式（２２）により、誤差Ｅｒｒ２の偏微分値を算出する。
・・・式（２２）
ここで、右辺の第１項は、注目ボクセルが密度勾配フィルタ計算の際のフィルタ中心位置にある項である。右辺の第２項は、注目ボクセルが密度勾配フィルタ計算の際のフィルタ中心位置にない項である。
ここで、右辺の第１項は、注目ボクセルが密度勾配フィルタ計算の際のフィルタ中心位置にある項である。右辺の第２項は、注目ボクセルが密度勾配フィルタ計算の際のフィルタ中心位置にない項である。
そして、式（２２）の右辺の第１項は、以下の式（２３）で定義される。
・・・式（２３）
ここで、Ｚは、上記式（１３）から式（１５）に示したものである。また、Ｘは、上記式（１２）で示したものである。例えば物質０の中に物質１がある一定の形状で存在している場合のｒｒ’、ｒ、ｒ⊥のそれぞれにおけるＥｒｒ２の偏微分の項は、以下の式（２４）、式（２５）、式（２６）のそれぞれで定義される。
・・・式（２４）
・・・式（２５）
・・・式（２６）
ここで、μ０は、物質０の中での物質１の形状を表す形状フィルタである。μ１は、物質０のみが存在している領域を表す形状フィルタである。Ｆ’は、上述した関数Ｆの偏微分である。また、式（２４）から式（２６）のそれぞれに含まれる
は、以下の式（２７）で定義される。
・・・式（２７）
また、式（２４）及び式（２５）のそれぞれに含まれる
は、以下の式（２８）で定義される。
・・・式（２８）
ここで、Ｚは、上記式（１３）から式（１５）に示したものである。また、Ｘは、上記式（１２）で示したものである。例えば物質０の中に物質１がある一定の形状で存在している場合のｒｒ’、ｒ、ｒ⊥のそれぞれにおけるＥｒｒ２の偏微分の項は、以下の式（２４）、式（２５）、式（２６）のそれぞれで定義される。
ここで、μ０は、物質０の中での物質１の形状を表す形状フィルタである。μ１は、物質０のみが存在している領域を表す形状フィルタである。Ｆ’は、上述した関数Ｆの偏微分である。また、式（２４）から式（２６）のそれぞれに含まれる
また、式（２４）及び式（２５）のそれぞれに含まれる
そして、式（２２）の右辺の第２項は、上記式（２３）から式（２８）においてｊをフィルタ中心のボクセル番号と読み替えて、以下の式（２９）で定義される。
・・・式（２９）
ここで、ｃａｓｅ１は、注目ボクセルｉがボクセルｊを中心とした密度勾配フィルタ計算に含まれない場合である。また、ｃａｓｅ２は、注目ボクセルｉがボクセルｊを中心とした密度勾配フィルタ計算に含まれる場合である。尚、式（２９）のｃ－１ ｒ／ｒｒ‘（ｊ）は式（１７）等で用いた座標変換関数ｃｒの逆関数である。
ここで、ｃａｓｅ１は、注目ボクセルｉがボクセルｊを中心とした密度勾配フィルタ計算に含まれない場合である。また、ｃａｓｅ２は、注目ボクセルｉがボクセルｊを中心とした密度勾配フィルタ計算に含まれる場合である。尚、式（２９）のｃ－１ ｒ／ｒｒ‘（ｊ）は式（１７）等で用いた座標変換関数ｃｒの逆関数である。
図１４は、本実施形態２にかかる形状フィルタμ０と復元（推定）された密度とが一致する場合を説明するための図である。形状フィルタμ０は、探査対象物２０の形状のうち中央に空洞がある場合のボクセル列に対応する形状情報の一例である。図１４の上部に示すように、形状フィルタμ０は、ボクセル列の長さ（ボクセル数）がＭであり、中央付近のボクセルの密度ρ０が両端のボクセルの密度ρ１よりも低くなるように定義されている。図１４の中央部は、探査対象物２０の形状として推定された密度（推定密度）を示す。この例は形状推定したボクセルｖｉに形状フィルタの中心を対応させて当てはめたものであり、形状フィルタμ０と推定密度の位置関係に誤差がないことを示す（図１４の下部）。
図１５は、本実施形態２にかかる形状フィルタμ０と復元（推定）された密度とのずれの調整の概念を説明するための図である。この例は形状推定したボクセルｖｉ’に形状フィルタの中心を対応させて当てはめたものであり、推定された探査対象物２０の位置と形状フィルタμ０の位置とに誤差（ずれ）があることを示す。そして、学習部４４５は、誤差を縮めるために、推定された探査対象物２０の位置をずらして形状フィルタμ０の位置と一致させたことを示す（図１５の下部）。
図１６は、本実施形態２にかかる他の形状フィルタμ１の例を示す図である。形状フィルタμ１は、探査対象物２０の形状のうちボクセルの密度が一定値ρ１であるケースを示す。
図１７は、本実施形態２にかかる重み付け関数を用いた、形状フィルタと復元（推定）された密度との境界のずれの調整の概念を説明するための図である。図１７は、上述した式（１３）や式（１４）の関数Ｆによる作用を示したものである。つまり、関数Ｆは、探査対象物２０の形状の推定値における境界位置と入力された形状フィルタの境界位置の相違を一定程度許容するものである。
続いて、学習部４４５は、第Ａ－１層から第１層まで順番に、一般的な誤差逆伝搬法により第Ａ－１層から第１層へ誤差の逆伝搬を行い、第Ａ－１層から第１層の間のパラメータを更新する。
その後、演算部４４６は、ステップＳ２０８で学習後（更新後）のパラメータを用いてＣＮＮ４１２の畳み込み演算を実行する（Ｓ２０３）。そして、上記同様に、推定装置４００は、ステップＳ２０４からＳ２０７を実行する。
誤差関数値が所定値未満である場合（ステップＳ２０７でＮＯ）、出力部４４８は、直近のステップＳ２０５で推定された推定値（三次元復元密度情報ｒｄ）を出力する（Ｓ２０９）。
このように本実施形態では、誤差関数に探査対象物の形状情報と推定値との誤差を最小化するための項（第２項）が追加されている。そのため、探査対象物の形状に関する事前知識を効果的に用いて、ill-posedなトモグラフィ問題を解くことができる。よって、上述した実施形態１の効果に加えて、推定値に基づくトモグラフィ結果の精度をさらに向上させることができる。
さらに、探査対象物の想定位置を用いて、探査対象物の形状に関する推定値の位置を補正するように機械学習を行う。よって、機械学習の効率、解の安定性、収束性がさらに向上し、推定精度がさらに向上する。
尚、本実施形態２は、次のような処理を行っても良い。例えば、受付部は、三次元空間のうち地表面より下の物質層の指定をさらに受け付けるとよい。その場合、設定部は、第１のボリュームのうち指定された物質層に対応するボクセルの初期値を、指定された物質層以外の物質層に対応するボクセルの初期値よりも高くするように設定する。これにより、探査対象外である地表面上のボクセルによる誤差の影響を抑制し、機械学習の効率がさらに向上する。
また、受付部は、複数の物質層のうち探査対象物が存在する可能性のある物質層の種別の指定をさらに受け付けるとよい。その場合、学習部は、各ボクセルのうち指定されなかった物質層に対応するボクセルを学習の対象外としてパラメータの学習を行う。すなわち、学習部は、指定外の地層に該当するボクセルをｄｒｏｐｏｕｔして学習を行う。これにより、機械学習の収束性及び安定性と、推定精度がさらに向上する。
＜その他の実施形態＞
尚、上述の実施形態では、ハードウェアの構成として説明したが、これに限定されるものではない。本開示は、任意の処理を、ＣＰＵにコンピュータプログラムを実行させることにより実現することも可能である。
尚、上述の実施形態では、ハードウェアの構成として説明したが、これに限定されるものではない。本開示は、任意の処理を、ＣＰＵにコンピュータプログラムを実行させることにより実現することも可能である。
上述の例において、プログラムは、様々なタイプの非一時的なコンピュータ可読媒体（non-transitory computer readable medium）を用いて格納され、コンピュータに供給することができる。非一時的なコンピュータ可読媒体は、様々なタイプの実体のある記録媒体（tangible storage medium）を含む。非一時的なコンピュータ可読媒体の例は、磁気記録媒体（例えばフレキシブルディスク、磁気テープ、ハードディスクドライブ）、光磁気記録媒体（例えば光磁気ディスク）、ＣＤ－ＲＯＭ（Read Only Memory）、ＣＤ－Ｒ、ＣＤ－Ｒ／Ｗ、ＤＶＤ（Digital Versatile Disc）、半導体メモリ（例えば、マスクＲＯＭ、ＰＲＯＭ（Programmable ROM）、ＥＰＲＯＭ（Erasable PROM）、フラッシュＲＯＭ、ＲＡＭ（Random Access Memory））を含む。また、プログラムは、様々なタイプの一時的なコンピュータ可読媒体（transitory computer readable medium）によってコンピュータに供給されてもよい。一時的なコンピュータ可読媒体の例は、電気信号、光信号、及び電磁波を含む。一時的なコンピュータ可読媒体は、電線及び光ファイバ等の有線通信路、又は無線通信路を介して、プログラムをコンピュータに供給できる。
なお、本開示は上記実施形態に限られたものではなく、趣旨を逸脱しない範囲で適宜変更することが可能である。また、本開示は、それぞれの実施形態を適宜組み合わせて実施されてもよい。
上記の実施形態の一部又は全部は、以下の付記のようにも記載され得るが、以下には限られない。
（付記Ａ１）
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付手段と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定手段と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成手段と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出手段と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習手段と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算手段と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定手段と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力手段と、
を備える推定装置。
（付記Ａ２）
前記受付手段は、前記探査対象物の形状に関する形状情報の入力をさらに受け付け、
前記誤差算出手段は、前記推定値のうち前記探査対象物の形状に関する推定値と、前記入力された形状情報との誤差をさらに含めて前記誤差関数値を算出する
付記Ａ１に記載の推定装置。
（付記Ａ３）
前記受付手段は、前記探査対象物の前記三次元空間内の想定位置の入力をさらに受け付け、
前記誤差算出手段は、前記探査対象物の形状に関する推定値に対して、前記想定位置とのずれを調整するように前記誤差関数値を算出する
付記Ａ２に記載の推定装置。
（付記Ａ４）
前記受付手段は、前記三次元空間のうち地表面より下の物質層の指定をさらに受け付け、
前記設定手段は、前記第１のボリュームのうち前記指定された物質層に対応するボクセルの初期値を、前記指定された物質層以外の物質層に対応するボクセルの初期値よりも高くするように設定する
付記Ａ１乃至Ａ３のいずれか１項に記載の推定装置。
（付記Ａ５）
前記受付手段は、前記複数の物質層のうち前記探査対象物が存在する可能性のある物質層の種別の指定をさらに受け付け、
前記学習手段は、各ボクセルのうち前記指定されなかった物質層に対応するボクセルを学習の対象外として前記パラメータの学習を行う
付記Ａ１乃至Ａ４のいずれか１項に記載の推定装置。
（付記Ａ６）
前記荷電粒子は、ミューオン粒子である
付記Ａ１乃至Ａ５のいずれか１項に記載の推定装置。
（付記Ｂ１）
コンピュータが、
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付け、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定し、
前記畳み込み演算により前記最終層の値を取得し、
前記取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成し、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出し、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習し、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得し、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定し、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する、
推定方法。
（付記Ｃ１）
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付処理と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定処理と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成処理と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出処理と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習処理と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算処理と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定処理と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力処理と、
をコンピュータに実行させる推定プログラムが格納された非一時的なコンピュータ可読媒体。
（付記Ａ１）
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付手段と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定手段と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成手段と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出手段と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習手段と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算手段と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定手段と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力手段と、
を備える推定装置。
（付記Ａ２）
前記受付手段は、前記探査対象物の形状に関する形状情報の入力をさらに受け付け、
前記誤差算出手段は、前記推定値のうち前記探査対象物の形状に関する推定値と、前記入力された形状情報との誤差をさらに含めて前記誤差関数値を算出する
付記Ａ１に記載の推定装置。
（付記Ａ３）
前記受付手段は、前記探査対象物の前記三次元空間内の想定位置の入力をさらに受け付け、
前記誤差算出手段は、前記探査対象物の形状に関する推定値に対して、前記想定位置とのずれを調整するように前記誤差関数値を算出する
付記Ａ２に記載の推定装置。
（付記Ａ４）
前記受付手段は、前記三次元空間のうち地表面より下の物質層の指定をさらに受け付け、
前記設定手段は、前記第１のボリュームのうち前記指定された物質層に対応するボクセルの初期値を、前記指定された物質層以外の物質層に対応するボクセルの初期値よりも高くするように設定する
付記Ａ１乃至Ａ３のいずれか１項に記載の推定装置。
（付記Ａ５）
前記受付手段は、前記複数の物質層のうち前記探査対象物が存在する可能性のある物質層の種別の指定をさらに受け付け、
前記学習手段は、各ボクセルのうち前記指定されなかった物質層に対応するボクセルを学習の対象外として前記パラメータの学習を行う
付記Ａ１乃至Ａ４のいずれか１項に記載の推定装置。
（付記Ａ６）
前記荷電粒子は、ミューオン粒子である
付記Ａ１乃至Ａ５のいずれか１項に記載の推定装置。
（付記Ｂ１）
コンピュータが、
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付け、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定し、
前記畳み込み演算により前記最終層の値を取得し、
前記取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成し、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出し、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習し、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得し、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定し、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する、
推定方法。
（付記Ｃ１）
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付処理と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定処理と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成処理と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出処理と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習処理と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算処理と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の推定値を推定する推定処理と、
前記誤差関数値が所定値未満である場合に、前記推定値を出力する出力処理と、
をコンピュータに実行させる推定プログラムが格納された非一時的なコンピュータ可読媒体。
以上、実施形態（及び実施例）を参照して本願発明を説明したが、本願発明は上記実施形態（及び実施例）に限定されものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
１ 推定装置
１１ 受付部
１２ 設定部
１３ 再構成部
１４ 誤差算出部
１５ 学習部
１６ 演算部
１７ 推定部
１８ 出力部
２ 三次元空間
２０ 探査対象物
２１～２６ 地層
２１１ 測定器
２１ｎ 測定器
３１～３４ ミューオン粒子の軌跡
４００ 推定装置
４１０ 記憶部
４１１ 推定プログラム
４１２ ＣＮＮ
４１２０ パラメータ
４１２１ 第１層
４１２Ａ 最終層
４１３ 物質種別
４１４ 固有密度係数
４１５ 最終層定義情報
４１５１ ボリュームＩＤ
４１５２ 測定器ＩＤ
４１５３ 物質種別
４２０ メモリ
４３０ ＩＦ部
４４０ 制御部
４４１ 受付部
４４２ 設定部
４４３ 再構成部
４４４ 誤差算出部
４４５ 学習部
４４６ 演算部
４４７ 推定部
４４８ 出力部
ｖｏ１ ボリューム
ｖｘ１１ ボクセル
ｖｘ１ｋ ボクセル
ｖｏ２１ ボリューム
ｖｏ２２ ボリューム
ｖｏ２３ ボリューム
ｖｏＡ－１１ ボリューム
ｖｏＡ－１２ ボリューム
ｖｏＡ－１３ ボリューム
ｖｏＡ１１ ボリューム
ｖｏＡ１Ｎｍ ボリューム
ｖｏＡＮｆ１ ボリューム
ｖｏＡＮｆＮｍ ボリューム
ｒ１ 再構成画像
ｒ２ 再構成画像
ｒＮｆ 再構成画像
ｏｂ１ ミューオンフラックス画像
ｏｂ２ ミューオンフラックス画像
ｏｂＮｆ ミューオンフラックス画像
ｒｄ 三次元復元密度情報
１１ 受付部
１２ 設定部
１３ 再構成部
１４ 誤差算出部
１５ 学習部
１６ 演算部
１７ 推定部
１８ 出力部
２ 三次元空間
２０ 探査対象物
２１～２６ 地層
２１１ 測定器
２１ｎ 測定器
３１～３４ ミューオン粒子の軌跡
４００ 推定装置
４１０ 記憶部
４１１ 推定プログラム
４１２ ＣＮＮ
４１２０ パラメータ
４１２１ 第１層
４１２Ａ 最終層
４１３ 物質種別
４１４ 固有密度係数
４１５ 最終層定義情報
４１５１ ボリュームＩＤ
４１５２ 測定器ＩＤ
４１５３ 物質種別
４２０ メモリ
４３０ ＩＦ部
４４０ 制御部
４４１ 受付部
４４２ 設定部
４４３ 再構成部
４４４ 誤差算出部
４４５ 学習部
４４６ 演算部
４４７ 推定部
４４８ 出力部
ｖｏ１ ボリューム
ｖｘ１１ ボクセル
ｖｘ１ｋ ボクセル
ｖｏ２１ ボリューム
ｖｏ２２ ボリューム
ｖｏ２３ ボリューム
ｖｏＡ－１１ ボリューム
ｖｏＡ－１２ ボリューム
ｖｏＡ－１３ ボリューム
ｖｏＡ１１ ボリューム
ｖｏＡ１Ｎｍ ボリューム
ｖｏＡＮｆ１ ボリューム
ｖｏＡＮｆＮｍ ボリューム
ｒ１ 再構成画像
ｒ２ 再構成画像
ｒＮｆ 再構成画像
ｏｂ１ ミューオンフラックス画像
ｏｂ２ ミューオンフラックス画像
ｏｂＮｆ ミューオンフラックス画像
ｒｄ 三次元復元密度情報
Claims (8)
- 探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付手段と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定手段と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成手段と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出手段と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習手段と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算手段と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の第１の推定値を推定する推定手段と、
前記誤差関数値が所定値未満である場合に、前記第１の推定値を出力する出力手段と、
を備える推定装置。 - 前記受付手段は、前記探査対象物の形状に関する形状情報の入力をさらに受け付け、
前記推定手段は、前記三次元空間に対応する立体空間内の各ボクセルにおける前記密度情報を推定し、全ボクセルの前記密度情報を前記三次元空間にまとめて、前記三次元空間における前記密度情報の分布の第１の推定値として推定し、
前記誤差算出手段は、前記第１の推定値から推定される前記探査対象物の形状に関する第２の推定値と、前記入力された形状情報との誤差をさらに含めて前記誤差関数値を算出する
請求項１に記載の推定装置。 - 前記受付手段は、前記探査対象物の前記三次元空間内の想定位置の入力をさらに受け付け、
前記誤差算出手段は、前記探査対象物の形状に関する第２の推定値に対して、前記想定位置とのずれを調整するように前記誤差関数値を算出する
請求項２に記載の推定装置。 - 前記受付手段は、前記三次元空間のうち地表面より下の物質層の指定をさらに受け付け、
前記設定手段は、前記第１のボリュームのうち前記指定された物質層に対応するボクセルの初期値を、前記指定された物質層以外の物質層に対応するボクセルの初期値よりも高くするように設定する
請求項１乃至３のいずれか１項に記載の推定装置。 - 前記受付手段は、前記複数の物質層のうち前記探査対象物が存在する可能性のある物質層の種別の指定をさらに受け付け、
前記学習手段は、各ボクセルのうち前記指定されなかった物質層に対応するボクセルを学習の対象外として前記パラメータの学習を行う
請求項１乃至４のいずれか１項に記載の推定装置。 - 前記荷電粒子は、ミューオン粒子である
請求項１乃至５のいずれか１項に記載の推定装置。 - コンピュータが、
探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付け、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定し、
前記畳み込み演算により前記最終層の値を取得し、
前記取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成し、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出し、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習し、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得し、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の第１の推定値を推定し、
前記誤差関数値が所定値未満である場合に、前記第１の推定値を出力する、
推定方法。 - 探査対象物を含む三次元空間を通過した荷電粒子の数と当該荷電粒子の入射方向とを測定する複数の測定器による測定値と、各測定器の設置位置と、前記三次元空間に含まれる複数の物質層の種別との入力を受け付ける受付処理と、
前記三次元空間を複数のボクセルの集合で表現した第１のボリューム内の各ボクセルに前記荷電粒子の密度情報として任意の値を設定し、前記第１のボリュームに対して所定の畳み込み演算を行う畳み込みニューラルネットワークにおける最終層に含まれる複数の第２のボリュームのそれぞれを、前記測定器と物質層の組合せに対応させて設定する設定処理と、
前記畳み込み演算により取得された最終層の値のうち、前記測定器ごとの各物質層に対応する前記第２のボリュームの値と、各物質層の種別に対応する物質の固有密度係数と、当該測定器に対応する前記設置位置とを用いて、前記測定器ごとに前記測定値の期待値を再構成する再構成処理と、
前記測定器ごとの前記再構成された期待値と前記測定値との二乗誤差を含む誤差関数値を算出する誤差算出処理と、
前記誤差関数値が所定値以上である場合に、前記誤差関数値を最小化するように、前記畳み込みニューラルネットワークのパラメータを学習する学習処理と、
前記学習後のパラメータを用いて前記畳み込み演算により最終層の値を取得する演算処理と、
前記取得された最終層の値と前記設置位置を用いて前記三次元空間における前記密度情報の分布の第１の推定値を推定する推定処理と、
前記誤差関数値が所定値未満である場合に、前記第１の推定値を出力する出力処理と、
をコンピュータに実行させる推定プログラム。
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| PCT/JP2020/041762 WO2022097305A1 (ja) | 2020-11-09 | 2020-11-09 | 推定装置、推定方法及び非一時的なコンピュータ可読媒体 |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JPWO2022097305A1 JPWO2022097305A1 (ja) | 2022-05-12 |
| JPWO2022097305A5 JPWO2022097305A5 (ja) | 2023-07-14 |
| JP7485076B2 true JP7485076B2 (ja) | 2024-05-16 |
Family
ID=81457708
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2022560628A Active JP7485076B2 (ja) | 2020-11-09 | 2020-11-09 | 推定装置、推定方法及び推定プログラム |
Country Status (2)
| Country | Link |
|---|---|
| JP (1) | JP7485076B2 (ja) |
| WO (1) | WO2022097305A1 (ja) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP7372611B1 (ja) | 2022-09-12 | 2023-11-01 | 日本電気株式会社 | 構造体観測装置、モデル構築方法及びモデル構築プログラム |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2013002830A (ja) | 2011-06-13 | 2013-01-07 | Kawasaki Geological Engineering Co Ltd | ミュー粒子を利用した三次元地盤探査システム |
| JP2018189428A (ja) | 2017-04-28 | 2018-11-29 | 国立大学法人東北大学 | ３次元像生成方法、３次元像生成システムおよび３次元像生成装置 |
| JP2019522803A (ja) | 2016-05-11 | 2019-08-15 | サントル ナシオナル ドゥ ラ ルシェルシェ シアンティフィクＣｅｎｔｒｅ Ｎａｔｉｏｎａｌ Ｄｅ Ｌａ Ｒｅｃｈｅｒｃｈｅ Ｓｃｉｅｎｔｉｆｉｑｕｅ | 岩石体積または人工建造物の密度を判断する方法および装置 |
-
2020
- 2020-11-09 WO PCT/JP2020/041762 patent/WO2022097305A1/ja active Application Filing
- 2020-11-09 JP JP2022560628A patent/JP7485076B2/ja active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2013002830A (ja) | 2011-06-13 | 2013-01-07 | Kawasaki Geological Engineering Co Ltd | ミュー粒子を利用した三次元地盤探査システム |
| JP2019522803A (ja) | 2016-05-11 | 2019-08-15 | サントル ナシオナル ドゥ ラ ルシェルシェ シアンティフィクＣｅｎｔｒｅ Ｎａｔｉｏｎａｌ Ｄｅ Ｌａ Ｒｅｃｈｅｒｃｈｅ Ｓｃｉｅｎｔｉｆｉｑｕｅ | 岩石体積または人工建造物の密度を判断する方法および装置 |
| JP2018189428A (ja) | 2017-04-28 | 2018-11-29 | 国立大学法人東北大学 | ３次元像生成方法、３次元像生成システムおよび３次元像生成装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| WO2022097305A1 (ja) | 2022-05-12 |
| JPWO2022097305A1 (ja) | 2022-05-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US8983162B2 (en) | Method and apparatus for estimating monte-carlo simulation gamma-ray scattering in positron emission tomography using graphics processing unit | |
| Barnes et al. | Geological analysis of Martian rover‐derived digital outcrop models using the 3‐D visualization tool, Planetary Robotics 3‐D Viewer—Pro3D | |
| CN109416408B (zh) | 震中距估计装置、震中距估计方法以及计算机可读记录介质 | |
| CN108072892B (zh) | 一种自动化的地质构造约束层析反演方法 | |
| Ding et al. | New parameterized model for GPS water vapor tomography | |
| CN112363236A (zh) | 一种基于pde的重力场数据等效源延拓与数据类型转换方法 | |
| WO2016001697A1 (en) | Systems and methods for geologic surface reconstruction using implicit functions | |
| JP7485076B2 (ja) | 推定装置、推定方法及び推定プログラム | |
| BR112015000879B1 (pt) | Sistema e método para modelagem de velocidade de migração | |
| Hou et al. | 3D density inversion of gravity gradiometry data with a multilevel hybrid parallel algorithm | |
| Raymund et al. | Model‐assisted ionospheric tomography: A new algorithm | |
| CN106471392B (zh) | 图像重构处理方法 | |
| Balco | Simple computer code for estimating cosmic-ray shielding by oddly shaped objects | |
| CN117011476B (zh) | 大气切伦科夫望远镜阵列的位型布局方法、设备及介质 | |
| CN112346139B (zh) | 一种重力数据多层等效源延拓与数据转换方法 | |
| CN116127314B (zh) | 基于自适应多尺度深度学习网络预测地下密度的方法 | |
| Hou et al. | 3D inversion of vertical gravity gradient with multiple graphics processing units based on matrix compression | |
| CN117092702A (zh) | 孔-隧激发极化探水结构的施工方法及反演探水方法 | |
| WO2021115917A1 (en) | Machine learning-based scintillator response modelling for increased spatial resolution in nuclear imaging | |
| Lesaint et al. | GCC and FBCC for linear tomosynthesis | |
| WO2012021938A1 (en) | A method of analysing data obtained using a gravity gradiometer | |
| KR20110124685A (ko) | Ｇｐｕ를 이용한 양전자 방출 단층 촬영 영상에서의 감마선 산란 추정 방법 및 장치 | |
| CN113591030B (zh) | 基于多gpu的重力梯度数据灵敏度矩阵压缩及调用方法 | |
| Yilmaz | Artificial neural networks pruning approach for geodetic velocity field determination | |
| Yin et al. | A fast 3D gravity forward algorithm based on circular convolution |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20230428 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20230428 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20231226 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20240221 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20240402 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20240415 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 7485076 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |