荷電粒子測定による統計的トモグラフィ再構成法
荷電粒子トモグラフィ・データから物体容積体散乱密度プロファイルを統計的に再構成する工程を含む荷電粒子検出用のシステム、装置、コンピュータ・プログラム製品及び方法であって、統計的多重散乱モデルに基づいて、荷電粒子散乱の確率分布を決定し、期待値最大化アルゴリズムを使用して物体容積体散乱密度の略最尤推定値を決定し、物体容積体の散乱密度を再構成する。注目の容積体を占める物体の存在と種類の両方又はそれらの一方を、最構成された容積体散乱密度プロファイルより特定することができる。荷電粒子トモグラフィ・データは、パッケージ、コンテナ、車両、又は貨物の走査用のミューオン追跡装置からの宇宙線ミューオントモグラフィ・データであってよい。このような方法は、コンピュータ上で実行されるコンピュータ・プログラムを使用して実施できる。
【発明の詳細な説明】
【技術分野】
【０００１】
本出願は、(1)「素粒子の検出と解析のためのシステム、方法及び装置、並びにフィールドへの展開」(SYSTEM, METHODS AND APPARATUS FOR PARTICLE DETECTION AND ANALYSIS AND DEPLOYMENT OF THE SAME)という名称の、２００６年１０月２７日出願の米国特許仮出願第６０/８５５、０６４号と、(2)「放射物ポータルモニターシステム及び方法」(RADIATION PORTAL MONITOR SYSTEM AND METHOD)という名称の、２００７年６月２９日出願の米国特許出願第１１/７７１、１６９号に付属する優先権を主張する。
【０００２】
上記の２つの出願の開示は、参照によりここに組み込んだものとする。
【０００３】
（米国連邦政府の権利に関する宣言）
本発明は、米国エネルギー省によって与えられた、契約番号DE-AC52-06NA25396に基づく政府支援でなされた。米国政府は、本発明に於いて一定の権利を保有する。
【０００４】
実施形態は、粒子の検出、解析、制御の分野に関し、特に、以下に限定されるわけではないが、複数の位置感知検出器を有する荷電粒子検出システムのデータを解析し、そして荷電粒子検出システムを通過する宇宙線ミューオンなどの荷電粒子の軌道を再構成する方法及び装置に関する。
【背景技術】
【０００５】
荷電粒子トモグラフィは、荷電粒子の散乱に基礎を置いている。荷電粒子トモグラフィの一形態は、宇宙線トモグラフィであり、宇宙線ミューオンの散乱を利用する。深宇宙からやって来る、安定な粒子、主にプロトンが、絶えず地球に降り注いでいる。これらの粒子は、上層大気圏の原子と相互作用して多数の短寿命のπオンを含む粒子シャワーを作り出し、それらπオンは崩壊し長寿命のミューオンを生み出す。ミューオンは、主にクーロン力により物質と相互作用し、核相互作用はせず、電子に比べると直ちに放射状に広がることはほとんどない。ミューオンは、電磁相互作用によって、ゆっくりとしかエネルギーを失わない。従って、ミューオンの多くは、透過度が高い荷電放射線として地球表面に到達する。海抜ゼロでのミューオン・フラックスは、約１ミューオン/ｃｍ2・ｍｉｎである。
【０００６】
ミューオンが物質を通過する時、原子構成要素の電荷とのクーロン散乱が、ミューオンの軌道を乱す。トータルの偏向は、いくつかの物質特性に依存するが、支配的なパラメータは、原子核の原子番号、Ｚ、と物質密度である。
【発明の概要】
【発明が解決しようとする課題】
【０００７】
注目する容積体を通過するミューオン又は他の荷電粒子からその容積体を再構成するための改良された方法及び装置を提供するニーズがある。
【課題を解決するための手段】
【０００８】
以下の発明の概要は、ミューオンのような荷電粒子などの粒子を検出し、荷電粒子トモグラフィ・データから物体容積体散乱密度プロファイルを統計的に再構成するための技術、装置、及びシステムに関するいくつかの技術的特徴の理解を容易にするために設けており、完全な説明を意図してはいない。本発明の様々な態様の完全な理解は、明細書全体、請求範囲、図面及び要約を全体として捉えることによって得ることができる。
【０００９】
本発明の上記の態様、他の目的及び利点は、ここに述べるように達成することができる。
【００１０】
一態様による、物体容積体を通過する荷電粒子を介して物体容積体を検出する検出システムが記載されている。このシステムは、物体容積体の第１の側に配置され、物体容積体に向かう入射荷電粒子の位置と角度を測定する、第１の位置感知検出器セットと、物体容積体の前記第１の側の反対側の第２の側に配置され、物体容積体から出て行く出射荷電粒子の位置と角度を測定する、第２の位置感知検出器セットと、前記第１の位置感知検出器セットからの測定信号のデータと前記第２の位置感知検出器セットからの測定信号のデータを受け取る信号処理ユニットとを含む。信号処理ユニットは、受信データを処理し、物体容積体内の容積体散乱密度分布を統計的に再構成する。
【００１１】
信号処理ユニットは、(ａ)物体容積体を通過する荷電粒子の散乱角と推定運動量に対応する荷電粒子トモグラフィ・データを取得し、(ｂ)統計的多重散乱モデルに基づいて、荷電粒子散乱密度の確率分布を用意し、(ｃ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて物体容積体散乱密度の略最尤推定値を決定し、(ｄ)前記略最尤推定値に基づいて再構成された物体容積体散乱密度を出力するように構成することができる。
【００１２】
別の態様によると、物体容積体から得られる荷電粒子トモグラフィ・データから物体容積体を検出する方法であって、(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、(ｂ)統計的多重散乱モデルに基づいて、荷電粒子散乱の確率分布を用意する工程と、(ｃ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）用いて、物体容積体散乱密度の略最尤推定値を決定する工程と、(ｄ)再構成された物体容積体散乱密度を出力する工程と、そして必要なら、(ｅ)再構成された物体容積体散乱密度に基づいて判定を行う工程と、を含む。
【００１３】
この方法によって、ユーザは、前記再構成された物体容積体散乱密度プロファイルより、注目する容積体を占める物体の存在と種類の両方又はいずれか一方を特定することが可能となる。いろいろな応用には、様々な自国保安の検査用途の宇宙線ミューオン・トモグラフィが含まれ、その用途では、車両又は貨物がミューオン追跡装置によって、走査され得る。
【００１４】
荷電粒子トモグラフィ・データは、宇宙線または他のなんらかの線源によって発生したミューオンなどの荷電粒子から収集されたトモグラフィ・データを含んでよい。
【００１５】
前記再構成された物体容積体散乱密度に基づいて判定を行う工程には、再構成された物体容積体散乱密度に基づいて、物体容積体内の目標物体の(1)存在と(2)種類の内のすくなくとも１つについて判定を行う工程を含めてよい。
【００１６】
期待値最大化アルゴリズム（ＭＬ/ＥＭ）で使用される、荷電粒子散乱の確率密度を用意する工程には、(ｇ)均一物体の所定の散乱密度に基づいて、荷電粒子の２次元確率分布を取得する工程と、(ｈ)前記２次元確率分布に基づいて、前記荷電粒子の３次元確率分布を取得する工程と、(ｉ)基底関数によって特徴付けられる非均一物体容積体による、多数の荷電粒子の散乱の確率分布を取得する工程と、(ｊ)その定義、及び荷電粒子の散乱と運動量の測定値に基づいて、多重散乱の確率分布を決定する工程と、を含めてよい。
【００１７】
均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元確率分布を取得する工程は、(ｋ)物質の散乱密度を、単位深さの物質による荷電粒子の散乱の２乗の期待平均値として決定する工程と、(ｌ)荷電粒子散乱分布をガウス分布モデルにより近似する工程と、(ｍ)荷電粒子線散乱及び変位を相関２次元ガウス分布により近似する工程と、を含んでよい。
【００１８】
２次元確率分布に基づいて、荷電粒子の３次元確率分布を得る工程は、座標を追加し３次元経路長を規定する工程と、３次元変位を計算する工程と、３次元共分散行列を規定する工程と、を含んでよい。
【００１９】
基底関数によって特徴付けられる非均一物体容積体による、多数の荷電粒子の散乱の確率分布を取得する工程は、離散的な散乱密度推定値を表す基底関数の３次元グリッドを設定する工程と、各個々のミューオンに対する散乱/変位の共分散行列を、宇宙線経路と散乱密度推定値の関数として決定する工程と、複数の荷電粒子に対する確率分布を個々の荷電粒子の確率の積として決定する工程と、を含んでよい。
【００２０】
前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の前記略最尤推定値を決定する工程は、各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、各荷電粒子と前記統計的散乱モデルの各基底関数との相互作用の幾何学を推定する工程と、荷電粒子と基底関数との各対に対して、重み行列：Ｗｉｊを決定する工程と、各基底関数の散乱密度を推測値で初期化する工程と、物体容積体内容に関して、近似の最尤解を反復して求める工程と、を含んでよい。ここで反復処理は、所定の反復回数で、又は解が所定の許容値未満に変わる時、停止する。
【００２１】
前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の略最尤推定値を決定する工程は、ｉ＝１からＭ(△θｘ、△θｙ、△ｘ、△ｙ、ｐｒ２)ｉまでの各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、各ミューオンとｊ＝１からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊを推定する工程と、荷電粒子とボクセルの各対に対して、
【数１】
【００２２】
として、重み行列：Ｗｉｊを計算する工程と、各ボクセルに於ける散乱密度λｊ,ｏｌｄの推測値を初期化する工程と、停止基準処理を使用してすべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗと置く工程と、を含む。
【００２３】
前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）は、平均値更新規則又は中央値更新規則を含んでよい。
【００２４】
さらなる別の態様によると、物体容積体から得られる荷電粒子トモグラフィ・データから物体容積体を検出する、コンピュータによって実行される方法であって、(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、(ｂ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）に使用される、統計的多重散乱モデルに基づく、荷電粒子散乱密度の確率分布を用意する工程と、(ｃ)前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の略最尤推定値を決定する工程と、(ｄ)再構成された物体容積体散乱密度を出力する工程と、を含む。判定は、前記再構成された物体容積体散乱密度に基づいて行うことができる。
【００２５】
さらなる別の態様によると、コンピュータ・プログラム製品は、コンピュータアクセス可能なデータ記憶媒体を含み、データ記憶媒体は、コンピュータによって実行される時荷電粒子トモグラフィ・データから物体容積体密度プロファイルを統計的に再構成する方法をコンピュータに実施させる命令を格納する。前記方法は、(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、(ｂ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）に使用される、統計的多重散乱モデルに基づく、荷電粒子散乱の確率分布を用意する工程と、(ｃ)前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体密度の略最尤推定値を決定する工程と、(ｄ)再構成された物体容積体散乱密度を出力する工程と、を含む。
【００２６】
添付図面於いては、同じ参照符号は、すべての図面を通じて、全く同じ又は機能的に同じ要素を参照している。添付図面は、組み込まれて明細書の一部を形成するが、それ以上に、本発明を図解し、本発明の詳細説明とともに、本発明の原理を説明する役割を果たす。
【図面の簡単な説明】
【００２７】
【図１】ミューオン・トモグラフィ概念の一例を図解している。
【図２】クーロン散乱を確定するのに使用される散乱と変位の２次元投影を図解している。
【図３】３次元散乱モデルを調整するのに使用される、散乱と変位の２次元投影のパラメータを図解している。
【図４】物質の多層レイヤによる散乱を図解している。
【図５】図３に示す投影の経路長計算に対する再近接点の使用を図解している。
【図６】ミューオン・トモグラフィの容積体散乱密度プロファイルを統計的に再構成する自動システムを図示している。
【図７】シミュレーションされた物体の斜視図である。
【図８】シミュレーションされた物体の上からみた図である。
【図９】ＰｏＣＡ点で接続される２個の直線を仮定して算定した経路長によるガウス分布散乱シミュレーションの再構成結果を図示している。
【図１０】非ガウス分布のテール部分の散乱によるシミュレーションデータの再構成結果を図示している。
【図１１】中央値法よる非ガウス分布のテール部分の散乱によるシミュレーションデータの再構成結果を図示している。
【図１２】シミュレーションされる乗用バンの主要物体を図解している。
【図１３】平均値法による乗用バンのミューオン暴露１分のシミュレーションの再構成結果を図示している。
【図１４】中央値法による乗用バンのシーンの再構成結果を図示している。
【図１５】一実施形態による、荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法の概略を説明しているフローチャートである。
【図１５Ａ】一実施形態による、多重・統計的散乱モデルを使用し、物体容積体を通過するミューオンの散乱の確率分布を推定する処理の一例の概略を説明しているフローチャートである。
【図１５Ｂ】一実施形態による、物体の所定散乱密度に基づいて、単一ミューオンに対する散乱の２次元期待確率分布を推定する処理の一例の概略を説明しているフローチャートである。
【図１５Ｃ】一実施形態による,統計的モデルの３次元への拡張する処理の一例の概略を説明しているフローチャートである。
【図１５Ｄ】一実施形態による,非均一物質による多数のミューオンの散乱及び変位の確率分布を決定する処理の一例の概略を説明しているフローチャートである。
【図１５Ｅ】一実施形態による,期待値最大化アルゴリズムを使用して物体容積体の推定密度プロファイルの尤度を最大化する処理の一例の概略を説明しているフローチャートである。
【発明を実施するための形態】
【００２８】
以下の非限定な例に於ける特定の値及び構成は、変更することができ、本発明の少なくとも一実施形態を説明するためにのみ出されており、発明の範囲を制限することを意図していない。
【００２９】
本願で説明している技術的特徴は、様々な粒子検出システムの構築に使用することができる。例えば、ミューオンを荷電粒子として検出する粒子検出システムは、検査対象の物体が配置される物体収容領域と、物体収容領域の第１の側に配置され、物体収容領域へ向かう入射ミューオンの位置及び角度を測定する、第１の位置感知検出器セットと、物体収容領域の第１の側の反対側の第２の側に配置され、物体収容領域から出ていく出射ミューオンの位置及び角度を測定する、第２の位置感知検出器セットと、第１の位置感知検出器セットからの入射ミューオンの測定信号のデータと第２の位置感知検出器セットからの出射ミューオンの測定信号のデータを受け取る、例えばマイクロプロセッサを備えてよい信号処理ユニットと、とを含むことができる。一例として、粒子検出器の第１及び第２のセットのそれぞれは、第１の方向の少なくとも３つの荷電粒子位置測定及び第１の方向と異なる第２の方向の少なくとも３つの荷電粒子位置測定が可能となるように配置されたドリフトチューブを含むように実現できる。信号処理ユニットは、ミューオンの測定された入射及び出射の位置及び角度に基づいて、物体収容領域内の物質におけるミューオン散乱によって生じるミューオンの散乱挙動を解析し、トモグラフィ・プロファイル又は物体収容領域内の散乱中心の空間分布を得るように構成される。得られたトモグラフィ・プロファイル又は散乱中心の空間分布は、物体収容領域における１つ以上の物体、例えば、核物質を含む大きな原子番号の物質や機器などの存否を明らかにするのに使用することができる。各位置感知ミューオン検出器は、ドリフトセル、例えば、ミューオンによってイオン化できるガスで充填されたドリフトチューブを含む様々な形状で実現できる。このようなシステムを使って、自然界の宇宙線ミューオンを、物体収容領域内の１つ以上の物体の検出用のミューオン線源として利用することができる。
【００３０】
一実施例においては、例示的実施形態による、荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法及びシステムは、物体を通過する宇宙線荷電粒子の散乱に基づいて、物体のイメージ又はモデルを再構成することが可能な解決法を提供する。
【００３１】
軌道は、よりありふれた物体（水、プラスチック、アルミニウム、スチールなど）を構成する物質より、特定の核物質（ＳＮＭ）及び良好なガンマ線遮蔽となる物質（鉛やタングステンなど）によって強く影響される。宇宙線荷電粒子、特に宇宙線ミューオンに関して、各ミューンは、それが透過した物体についての情報を運んでおり、多数のミューオンの散乱を測定することによって、これらの物体の特性を調べることができる。特に、より一般的な低Ｚや中間Ｚの物質の中の高Ｚの物体を検出することができる。
【００３２】
例示的実施形態による、荷電粒子トモグラフィの容積体密度プロファイルの統計的な再構成法の様々な技術的特徴を説明するために、先ず、図１にその例が図解してあるミューオントモグラフィ・コンセプトについて述べる。
【００３３】
位置感知検出器のセット10は、撮像対象の物体容積体の上方と下方に配置され、入射及び出射荷電粒子の飛跡12（矢印付き実線で示されている）の位置と角度を提供する。撮像対象の容積体上方に配置されている位置感知検出器10の２個以上のセットが、入射荷電粒子の飛跡の位置と角度を与える。これらの検出器は、２つの直交又は非直交座標に於ける荷電粒子の位置を測定する。位置感知検出器10の別のセットは出射荷電粒子の位置と角度を記録する。側面の検出器（図示せず）を使って、より水平に向いた荷電粒子飛跡を検出してもよい。各荷電粒子飛跡の散乱角は、符合する入射及び出射測定値より算出される。荷電粒子運動量は、検出器自身内で起こる僅かな散乱より算定される、又は位置感知検出器プレーンの２個のセット間に置かれた、既知の特性の散乱体のレイヤ内で起こる僅かな散乱より算定される。
【００３４】
位置感知荷電粒子検出器の一例は、作動ガスが充填されたドリフトチューブである。ドリフトチューブは、円筒状の管でよく、ミューオンなどの宇宙線荷電粒子の検出を可能とする、アルゴン・イソブタンなどの検出器ガスが充填される。管の外部表面をアースし、円筒状管の長さ方向に沿って延びている中央の陽極ワイヤに約２〜３ｋＶの正の高電圧が印加され、高電圧の静電場が生じる。荷電粒子がガス原子達と相互作用すると、管の弦(ｃｈｏｒｄ)によって一直線上の原子達から多くの電子が解放される。電子達の「ひも」は、静電場によって、正に帯電した陽極ワイヤへ向かってドリフトし、データ取得電子機器のＴＤＣＳ（時間デジタル変換器）で電子的に読み出される。検出器の各セットは、第１の方向の少なくとも３つの荷電粒子位置測定と、第１の方向とは異なり、第１の方向と直交してよい第２の方向の少なくとも３つの荷電粒子位置測定とを可能とするように配置されたドリフトチューブ達であってよい。
【００３５】
図１のシステムには、信号処理ユニット、例えばコンピュータが備えられ、物体容積体の上方の検出器による入射ミューオンの測定信号データと、物体容積体の下方の検出器による出射ミューオンの測定信号データとが受け取られる。この信号処理ユニットは、測定されたミューオンの入射及び出射の位置並びに角度に基づいて、容積体内の散乱によるミューオンの散乱挙動を解析し、トモグラフィ・プロファイル又は容積体内の散乱中心の空間分布を取得するように構成される。取得されたトモグラフィ・プロファイル又は容積体内の散乱中心の空間分布によって、容積体内の物体の存否を明らかにすることができる。ある実施例に於いては、追加のドリフトチューブ検出器を容積体の側面に配置して、システムの走査のために、パッケージ、乗物又は貨物コンテナを入れることができるボックス又は四面構造を形成することができる。このようにして、宇宙線ミューオンの多重散乱を使用して、通常の貨物のバックグラウンド内の高Ｚ物質(高原子番号)を選択的に検出することができる。利点として、この技法は、受動的であり、バックグラウンドを超える如何なる放射線量も加えず、高Ｚ高密度物質について選択的である。信号処理ユニットのトモグラフィ処理部は、検出器10と同じ場所にある現場のコンピュータに組み込んでもよい。あるいは、信号処理ユニットのトモグラフィ処理部は、プライベート・ネットワークのようなコンピュータ・ネットワーク、又はインターネットのような公共ネットワークに接続されているリモート・コンピュータに組み込んでもよい。
【００３６】
図１の例示的実施形態においては、荷電粒子は、宇宙線ミューオン又は他の宇宙線荷電粒子であり、位置感知検出器10は、荷電粒子を感知するための作動ガスが充填されたドリフトセルである。ドリフトセルは、例えば、中央の陽極ワイヤが管の長手方向に沿って走っているドリフトチューブで実現してよい。しかし、ミューオン以外の荷電粒子は、ドリフトセル以外の位置感知センサーを用いて検出してもよい。さらに荷電粒子は、宇宙線以外の線源によって発生してもよい。例えば、ミューオンは、加速器から低強度のビームとして発生することができる。
【００３７】
高密度物体を透過するミューオン（ブラックの飛跡）は、大気中を走るミューオン（グレーの飛跡）よりかなり強く散乱する。多数の飛跡測定より、物体形状と物質の電子密度を再構成することができる。容積体を通過するミューオンは、通過する物質に依存して散乱される。
【００３８】
図１のシステム用処理ユニットによる、検査中の容積体（例えば、パッケージ、コンテナ又は車両）に於ける宇宙線ミューオンの測定値の処理には、容積体を通るミューオンの軌道を再構成する工程と、容積体の各側方にある検出器からの信号を基に、入射ミューオンの運動量を測定する工程と、容積体の散乱密度の空間分布を決定する工程と、を含めてよい。これら及び他の処理結果を使用して、トモグラフィ・プロファイルを構築し、目標物体の検出などのような容積体の様々な特性を測定することができる。
【００３９】
例えば、ドリフトチューブ一式を有する検出器10を通過する荷電粒子の軌道の再構成には、(ａ)荷電粒子が衝突したドリフトセルの識別子と対応する衝突時刻を表す衝突信号を取得する工程と、(ｂ)検出器を通過する特定の荷電粒子の飛跡に対応していると同定された時間内のドリフトセル衝突をグループ化する工程と、(ｃ)前記特定の荷電粒子がドリフトセルに衝突する瞬間の時刻ゼロ値を最初に推定する工程と、(ｄ)時刻ゼロ値の推定値、ドリフト時間変換データ及び衝突時刻に基づいてドリフト半径達を決定する工程と、(ｅ)特定の時刻ゼロ値に対応するドリフト半径達に直線の飛跡をフィティングする工程と、(ｆ)特定の荷電粒子に対して行われたもっとも良い飛跡フィットに関連する時刻ゼロ値を検索して選択し、時刻ゼロと追跡パラメータに於ける誤差を算出する工程と、を含めてよい。時刻ゼロ・フィティングによるこのような飛跡の再構成法に於いては、高速な検出器（シンチレータ・パドル付光電子増倍管など）を、或いはミューオンの装置内の通過を数ナノ秒の単位で検出して時刻ゼロ値を与える何か他の高速検出器を使用することなく、荷電粒子検出器を通過する荷電粒子の再構成された直線軌道が得られる。
【００４０】
又、例えば、図１の検出器10からの信号に基づいて、入射ミューオン又は出射ミューオンの運動量を測定する処理には、(ａ)複数の位置感知検出器を配置し、通過する荷電粒子を散乱させる工程と、(ｂ)位置感知検出器内の荷電粒子の散乱を測定する工程と、(ｃ)位置測定値から少なくとも１つの荷電粒子の軌道を決定する工程と、(ｄ)前記少なくとも１つの軌道から荷電粒子のすくなくとも１つの運動量測定値を決定する工程と、を含めてよく、前記散乱の測定は少なくとも散乱荷電粒子の３つの位置測定値を取得する工程を含む。この技法を使用して、検出器内に追加のメタルプレートを使用することなく、位置感知検出器自身内での荷電粒子の散乱より決定される荷電粒子の軌道に基づいて、荷電粒子の運動量を決定することができる。
【００４１】
以下に、荷電粒子トモグラフィ・データからの物体容積体散乱密度プロファイルを統計的に再構成する例示的システム及び方法の詳細を示している。
【００４２】
一実施形態による、荷電粒子トモグラフィの物体容積体散乱密度プロファイルを統計的に再構成する自動化されたシステムの一例が、図６のブロック図に示してある。自動システム50は、荷電粒子トモグラフィ・データ54を受け取るように適合され構成されたコントローラ51を有している。荷電粒子トモグラフィ・データは、例えば、図１の荷電粒子検出器、或いは、容積体を通過する荷電粒子の追跡記録ができるように構成された位置感知検出器を有する他の荷電粒子検出器を使ったミューオンの測定値より決定されたミューオンのトモグラフィ・データでよい。結果として、ミューオン又は他の荷電粒子のトモグラフィ・データを使って、物体容積体を通過するミューオン又は他の荷電粒子の散乱角及び推定運動量を得る又は決定することができる。
【００４３】
自動システム50は、コントローラ上に格納された統計的再構成モジュール56を含む。再構成モジュール56は、ミューオン又は他の荷電粒子のトモグラフィの容積体散乱密度プロファイルの統計的再構成を司る。再構成モジュール56は、ソフトウエアモジュール又はハードウエアモジュールとして実現してよい。
【００４４】
図６の自動システム50の例示的実施形態に於いて、コントローラ51は、コンピュータ(ＰＣ)のような１つ以上の動作可能に連結された計算プロセッサユニット(ＣＰＵ)ベースのシステム、又はデジタル信号プロセッサベースのシステムのような他のマイクロプロセッサベースのシステムを使って、構成される。コントローラは、単一の通常のコンピュータでもよいが、しかし、実時間的に結果を得るために、コントローラは、典型的には、実時間で結果を得るのに必要な処理能力を提供するように、十分な数の並列処理コンピュータ（図示せず）という形態を含む。例えば、コントローラは、約２０個のＣＰＵを含んでよい。ミューオン検出器の走査容積体が大きくなるほど、そして所望の分解能が細かくなるほど、処理コンピュータの形態は大きくなる必要がある。
【００４５】
オペレーティングシステムは、コントローラ51上で動作し、そして、限定するものではないが、Ａｐｐｌｅ，Ｗｉｎｄｏｗｓ(登録商標)，Ｌｉｎｕｘ，又はＵｎｉｘ(登録商標)のオペレーティングシステム、又は将来開発される他のオペレーティングシステムを含め、市販又はオープン・ソースのオペレーティングシステムでよい。オペレーティングシステム、及びアプリケーション又はプログラムの命令は、ハードディスク・ドライブなどの記憶装置に格納される。また、自動システム50に於いては、飛跡の再構成モジュールは、コンピュータアクセス可能なデータ記憶媒体の形式のソフトウエアであって、データ記憶媒体は、コントローラによって実行される時例示的実施形態による荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法をコントローラに行わせしめる命令を格納している。モジュールは、図６に示されるように、現場のコントローラにインストールしてもよいし、コントローラに接続されたインターネット経由で遠隔地から起動してもよい。同業者なら、このようなモジュールの実現形態が複数あることは理解する。
【００４６】
自動システム50は、又、操作可能にコントローラ51に接続されたディスプレイ58を含む。ディスプレイ58はシステムによって再構成された物体密度プロファイルの画像又はデータを必要に応じて、ユーザに表示する。ユーザインタフェース（図示せず）が処理システムに操作可能に接続され、人のオペレータが必要に応じて、処理システムを操作できるようにしてもよい。
【００４７】
同業者であれば、図６の説明は、単に自動システム50の実施形態の一例を描いているにすぎず、実施形態はそれに限定されないことが理解できる。例えば、再構成モジュールの機能の一部又はすべては、マイクロプロセッサを使うことなく、アナログ又はデジタル回路などのハードウエアとして実現できる。
【００４８】
図１５を説明する。これは、一実施形態による、荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法の概要を説明するフローチャートを示している。方法100は、処理ステップ１０１に示されているように、物体容積体を通過する荷電粒子の飛跡位置、散乱角、推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得することで始まる。所定の荷電粒子トモグラフィ・データは、例えば、図１の検出器から取得できる。その後、処理ステップ１０２に示されているように、散乱密度の空間分布（以降に定義される）によって表される、物体容積体を通過する複数の荷電粒子の散乱の確率分布が、多重統計散乱モデルを基づいて、用意される。次に、処理ステップ１０３に示すように、期待値最大化アルゴリズムを使用して、物体容積体散乱密度プロファイルの最尤推定値が決定される。それから、再構成された容積体散乱密度プロファイルが判定用に出力される（処理ステップ１０４）。判定処理は、任意であり、処理ステップ１０５に示すように容積体を占有している物体の存在と種類の両方又はいずれか一方を特定する処理であってよい。判定処理は、物体容積体の再構成された密度プロファイルを表す画像の人による解釈と、さらなるアルゴリズムによる自動的な判定との両方又はそれらの一方を含んでよい。
【００４９】
本実施形態の方法と自動システムによって、多数の荷電粒子によって与えられるデータに基づいて、注目の容積体の離散的なトモグラフィ再構成が可能となる。反復的な期待値最大化(ＥＭ)アルゴリズムのインスタンスを使用して、物体の密度プロファイルの最尤推定値を求める。本実施形態の方法とシステムによって、ユーザは、再構成された容積体密度プロファイルから、注目の容積体を占める物体の存在と種類の両方又はそれらの一方を特定することが可能となる。いろいろな応用には、様々な自国保安の検査用途の宇宙線ミューオン・トモグラフィが含まれ、その用途では、車両又は貨物がミューオン追跡装置によって、走査される。得られるミューオンのトモグラフィ・データは、危険物を識別できる本例示的実施形態の方法と自動システムを用い、車両又は貨物の密度プロファイルを再構成し表示するのに使用できる。
【００５０】
最尤法は、医療画像再構成に於いて、特にＰＥＴ及びＳＰＥＴ再構成のために使用されているが、いくつかの重要な差異のために、それらの用途に開発された標準的な方法は使用できない。第１に、測定信号−散乱角−は確率的であり、平均値ゼロで、標準偏差は透過物質の特性によって規定される。第２に、宇宙線ミューオンは、確定した離散的な方向又は角度から来ることはなく、むしろ天頂の周りに幅広い、連続的な角度分布を持ち、その分布は殆ど水平にまで広がっている。最後に、ミューオンの軌道は、直線ではない。我々が、強い散乱をおこす物体のおおよその場所を見つけることができるのは、その湾曲である。ＥＭアルゴリズムは、柔軟性があり、計算効率が良い。そのアルゴリズムの複雑な幾何形状への適用について説明することができる。
【００５１】
これから、一実施形態に従って、処理ステップ１０２から１０４について説明する。データは、容積体を通過するミューオンを測定する図１の検出器から取得される宇宙線ミューオンのトモグラフィ・データである。
【００５２】
一実施形態による、多重統計散乱モデルを使用して物体容積体を通過するミューオンの散乱の推定確率分布を用意する処理（処理ステップ１０２）が、図１５Ａのフローチャートに概説されている。処理ステップ１１０から１１３で示されているように、本処理は、４つの主な構成部分を有している。最初に、均一物体の所定の散乱密度に基づいて単一ミューオンに対する２次元の確率分布が推定される（処理ステップ１１０）。そして２次元分布は３次元に拡張される（処理ステップ１１１）。次に、処理ステップ１１２に於いて、非均一物体容積体がボクセル基底関数を使用して表現され、ボクセル化された散乱密度を所与として、多数のミューオンの散乱の確率分布が表現される。最後に、確率分布表現は、ミューオン散乱及び運動量の測定値の有限精度に拡張される（処理ステップ１１３）。
【００５３】
処理ステップ１１０から１１３は、多重散乱統計モデルを使って実現される。そのモデルは、先ず均一物質の単一レイヤに於ける散乱に関して、そして非均一物質に於ける散乱に関して説明される。
【００５４】
物質を通過する宇宙線ミューオンは、図２に図解されているように多重クーロン散乱を受ける。図２は、多重クーロン散乱を記述するために使用される散乱と変位の２次元投影を図示している。この図及び他の図では、散乱の大きさは説明のために大幅に誇張されている。出射ミューオンは、入射ミューオンの方位と位置に対して測られた、散乱角と変位によって特徴付けることができる。典型的な散乱角は、数十ミリラジアン（１ミリラジアン≒０．０６度）であり、数度より大きい散乱角は非常にまれである。散乱角の中央９８%の分布は、平均値ゼロのガウス分布として近似してよい。
【数２】
【００５５】
実際の分布は、ガウス分布より重い或いは大きなテール部分をもっている。分布の幅は、物質特性の関数として近似的に表わすことができる。S.Eidelman et al.,"Review of particle physics" Phy.Lett.,vol.B592,p.1,2004にレビューされているように、多くの研究者が、散乱に対して、様々な物質の特性の関数として実験的に改良した式を提案している。その開示は、参照によりここに組み込んだものとする。特に単純な形は
【数３】
【００５６】
である。
【００５７】
ここで、ｐは単位がＭｅＶ/ｃの運動量、Ｈは物質の深度(ｄｅｐｔｈ)、Ｌｒａｄは物質の放射長、βｃは速度（ｃは光の速度）であり、β＝１の近似が使用される。放射長は、原子番号及び物質密度が大きくなるに従い減少する。公称のミューオンの運動量をｐ0とし、放射長Ｌｒａｄの物質の散乱密度を次式で定義する。
【数４】
【００５８】
物質の散乱密度、λ、は、その物質の単位深度を通過する公称運動量を持つミューオンの散乱角の２乗平均値を表す。いくつかの物質のミリラジアン2/センチメートル単位での散乱密度の値は、例えば、アルミニウムでは約３、鉄で１４、ウラニウムで７８である。それで、散乱密度λ、深度Ｈの物質を通過する運動量ｐのミューオンの散乱の分散は、
【数５】
【００５９】
である。ここで、
【数６】
【００６０】
と置くと、
【数７】
【００６１】
を得る。
【００６２】
変位△ｘは、散乱角△θと相関性がある。一緒に扱うと、散乱角と変位は、図１の経路の“屈曲(ｋｉｎｋ)”によって示唆されているように、大きな容積体中の局所散乱の寄与物の位置を示唆する情報を提供する。散乱角と変位の分布は、平均値がゼロで、
【数８】
【数９】
【００６３】
を持つ同時確率ガウス分布として特徴付けることができる。その共分散行列は次式で表すことができる。
【数１０】
【００６４】
ここで、
【数１１】
【００６５】
と置くと、
【数１２】
【００６６】
を得る。
【００６７】
上記を考慮して、一実施形態による、単一のミューオン散乱の２次元に於ける散乱確率分布の取得（処理ステップ１１０）は、図１５Ｂのフローチャートの概略のように説明できる。処理ステップ１５０に示されているように、物質の散乱密度は、式（３）によって、単位深度の物質によるｐ0＝３ＧｅＶ/ｃのミューオンの散乱の２乗の期待平均値として定義される。次に、処理ステップ１５１に示されているようにＲＭＳ散乱に対してガウス分布近似、式(１)、（５）,（６）が行われる。最後に、宇宙線の散乱及び変位の分布は、処理ステップ１５２に示されているように、式(10),(11)に要約されている平均値ゼロの相関２次元ガウス分布によって近似される。
【００６８】
３次元に於いては、散乱は、ｘ座標に直交するｙ座標を考慮することによって特徴付けられる。散乱角△θｘと△θｙ、そして変位△ｘと△ｙを参照されたい。ｘ平面とｙ平面への偏向は互いに独立で同一の分布に従う（Eidelman et alを参照）。前記展開は、入射ミューオンの方向に直交するように置かれた座標系に基づいている。３次元モデルにおいては、経路長を考慮し、変位測定値を入射ミューオン経路に直交する平面に投影する必要がある。図３は３次元散乱に対してモデルを適合させるのに使用されるパラメータを図示しているが、垂直方向に対して、投影角度θｘ,0で入射するミューオンが示されている。
【００６９】
この３次元散乱の理解の助けとなるように、図のページから外に向く直交するｙ座標における対応する投影角度θｙ,0を想像することは有用である。レイヤを通って投影（不散乱）点（ｘｐ、ｙｐ）までのミューオン経路の直線延長部（即ち、３次元経路長）は、
【数１３】
【００７０】
である。
【００７１】
出射ミューオンのｘ位置と角度を（ｘ１、θｘ,1）と定義し、そして
【数１４】
【００７２】
と置く。
【００７３】
測定されるｘ変位は、ｘｍ＝ｘ1−ｘｐとして計算されるであろう。しかし、この測定値を回転させて宇宙線経路に直交する平面に入れ、３次元経路長に対して調整する必要がある。変位は、次で定義される。
【数１５】
【００７４】
ここで、中央の２項は３次元経路長に相当し、最後の項は、測定値を適切な方向に投影するものである。
【００７５】
最後に、共分散重み行列を次のように再定義する。
【数１６】
【００７６】
そして散乱と変位に対して同様に進める必要があり、式(11)は、両方かつ座標散乱に対する共分散行列を定義している。散乱測定は、２つの直交、水平座標で独立に行われる。表記法を簡単にするために、ただ1つの座標について解析を展開する。２つの座標からの情報を結合するのは後で説明する。本モデルは、“小”散乱角と変位に対して有効であることに留意すべきである。大きな角度の散乱では、モデルの導出に於いて無視された２次の項が有意となり得る。
【００７７】
上記を考慮すると、一実施形態による図１５Ｃのフローチャートに於いて概説されているように、３次元に拡張された統計的モデルが得られる（処理ステップ１１１）。先ず、ｙ座標が追加され、３次元経路長が規定される（処理ステップ１６０、式(12）。次に、処理ステップ１６１に於いて、式(13),(14)によって３次元変位が計算される。最後に、式(15)によって、３次元共分散行列を表す（処理ステップ１６２）。
【００７８】
物質の非均一容積体に対しては、密度プロファイルは、再構成のために、係数{ν1、・・・νｊ、νＮ}と３次元の基底関数{φ1、・・・φｊ、ΦＮ}の線形結合として表現され、即ち；
【数１７】
【００７９】
である。
【００８０】
基底関数に対しては、多くの選択肢があるが、ここでは、矩形の３次元ボクセルに注意を払う。λｊは、第ｊ番目の基底関数の係数、即ち、第ｊ番のボクセルの散乱密度を示すのに使用される。図４を考えると、３つのレイヤ（又はボクセル）が示されており、宇宙線がスタックを通過し、観測情報△θと△ｘを届けている。第ｊ番のボクセルの“見えない”散乱及び変位がそれぞれ△θｊと△ｘｊとして示されている。又、図では、散乱の大きさは誇張されている。以下の式によって、観測されるデータを見えないデータに関係付けることができる。
【数１８】
【数１９】
【００８１】
ここで、２番目の式では、小角散乱という仮定に依拠しており、Ｔｊは、第ｊ番のボクセルの出口点から再構成容積体の出口点までの３次元宇宙線経路長として定義している。より一般的には、ボクセル集合χを通過する宇宙線に対して、
【数２０】
【数２１】
【００８２】
である。
【００８３】
最後に、先ず、第ｊ番のボクセルに対して、次式であることに留意することによって、総散乱/変位に対する共分散を表現することができる。
【数２２】
【００８４】
ここで、
【数２３】
【００８５】
であり、Ｌｉｊは、第ｊ番のボクセルを通過する第ｉ番の宇宙線の経路長であり、宇宙線の当らないボクセルに対しては、ゼロと規定される。式(19)から(22)を組み合わせると、次式を得る。
【数２４】
【００８６】
ここで、Ｎは全ボクセル数であり、各要素に対する単純ではあるが冗長な計算に基づいて、次の重み行列を定義している。
【数２５】
【００８７】
ボクセルを通る宇宙線の経路長を算定するために、未知のミューオン経路について仮定を行う。図５を参照すると、近似は、入射及び出射飛跡の再近接点(ＰｏＣＡ)の計算で始まる（ｘｃａ、ｙｃａ）。次に、入り口点とＰｏＣＡ、ＰｏＣＡと出口点を接続し、ボクセル経路長を算定する。
【００８８】
最後に、次のデータベクターを定義し、
【数２６】
【００８９】
ＤをＭ個のミューオンのすべての測定値を示すものとすると、密度プロファイルλが与えられたときのデータの尤度は次のように書ける。
【数２７】
【００９０】
上記の因子は、
【数２８】
【００９１】
である。
【００９２】
上記を考慮すると、図１５Ｄに概説されているように、一実施形態による非均一物質による多数のミューオンの散乱及び変位の確率分布を得ることができる（処理ステップ１１２）。先ず、ボクセル（又は他の基底関数）の３次元グリッドを定める（処理ステップ１７０）。次に処理ステップ１７１に於いて、各ミューオンに対して、散乱/変位の共分散行列が式(23),(24)によって計算される。最後に、所与の宇宙線経路とボクセル散乱密度に対して、すべてのミューオンに対する全体確率分布が式(25)から(27)によって計算される。これは処理ステップ１７２として示されている。
【００９３】
多重散乱統計モデルを説明したので、これから、実験効果に対するモデルの拡張について説明する（処理ステップ１１３）。現実のミューオン検出器は、有限な位置分解能を持っている。入射及び出射ミューオン飛跡は、多数の位置測定値に対する飛跡フィッティングにより得られた角度と位置によって特徴付けられる。従って、測定誤差が、ミューオントモグラフィ用データセットを構成する散乱角と変位の測定値に伝播する。特定の検出器の精度は、ＲＭＳ誤差ｅｐによって特徴付けられる。検出器の特定の配置に対して、誤差行列
【数２９】
【００９４】
は、誤差がどのように伝播するかに基づいて規定してよい。このような誤差は、反復的な再構成法では、扱いが比較的容易である。本特許の場合は、式(23)の共分散行列に追加することによって検出器誤差を考慮することができる。
【数３０】
【００９５】
このようにして、検出器誤差に起因するノイズは削減される。そうしないと、検出器誤差に起因するノイズが、再構成結果に現われる。検出器誤差に対するより正確なモデルは、運動量依存性を考慮しなければならない。追跡誤差源の1つは、検出器自身内の散乱であり、散乱は素粒子の運動量が増加すると、減少するからである。個々のミューオン運動量の推定値、
【００９６】
が利用できるなら、誤差行列
【００９７】
が各宇宙線に対して算定することができる。式(２)から明らかように、多重クーロン散乱分布の幅は、素粒子の運動量に依存する。異なるミューオンの運動量は、式(５)に因子ｐｒ2を導入することによって考慮されている。実際には、ミューオンの運動量は正確には知られていないが、個々のミューオンの運動量の推定値は、宇宙線ミューオンの既知のスペクトルなどの、既知の散乱体に於ける散乱の測定値より算定することができる。ここでは、各ミューオンに対して
【００９８】
という良い推定値があると仮定する。
【００９９】
物体容積体密度の最尤推定値は、期待値最大化アルゴリズムを使用して決定できる（方法100の処理ステップ１０３）。ＥＭアルゴリズムは、“不完全な”データの尤度を“全部揃った”データとして、即ち“見えない”データを加えた観測データとして表現することに依拠している。この適用では、観測データＤ＝｛Ｄｉ：１≦ｉ≦Ｍ｝は、測定された散乱である。見えないデータＨ＝｛Ｈｉｊ：１≦ｉ≦Ｍ＆１≦ｊ≦Ｎ｝は、第ｊ番のボクセルによる第ｉ番のミューオンの散乱角及び変位である。A.Dempster, N.Laird, D.Rubin,"Maximum likelihood from incomplete data via the EM algorithm", J.Roy.Statist.Soc.B, vol.39, pp.1-78,1977に、以下の補助関数によってアルゴリズムが説明されている。本文献の開示は参照によりここに組み込んだものとする。
【数３１】
【０１００】
この関数は、与えられたＨの条件付分布とパラメータベクトルλ(ｎ)に対して、パラメータベクトルλが与えられたときの観測データと未観測データの対数尤度の期待値である。アルゴリズムの各繰り返しは次の２ステップから構成される。
【０１０１】
Ｅステップ：見えないデータの条件付き分布Ｐ(Ｈ|Ｄ,λ(ｎ)）を推定又は特性化する。
【０１０２】
Ｍステップ：補助関数Ｑを最大化する。Ｑは、Ｅステップで特性化した分布に対する期待値である。
【０１０３】
本特許の場合、見えないデータは、観測データを一意に決めるので、より簡単な補助関数
【数３２】
【０１０４】
を使うことによって、ＱＤＬＲを使うことによって得られるのと同じ推定値列
【０１０５】
を得る。パラメータ推定値、λ(ｎ)よりアルゴリズムの1回の繰り返しにより新しい推定値、λ(ｎ+1)が以下により得られる。
【数３３】
【０１０６】
先ず、単一ボクセルによる単一ミューオンの散乱の確率分布は、単純に単一レイヤに対する統計的モデルから得られることに留意する。
【数３４】
【０１０７】
ここで、
【数３５】
【０１０８】
は、式(21)で定義されている。各ボクセルに於ける散乱の無条件確率分布は、他のボクセルの散乱とは独立であるので、見えないデータの全体集合の確率は、各要素の確率の積である。従って、対数尤度は、
【数３６】
【０１０９】
と書くことができる。ここでＣはλを含まない項を表している。条件付期待値を計算することによって、関数Ｑは、
【数３７】
【０１１０】
となり、その各加数は
【数３８】
【０１１１】
である。ここでＨｊは、Ｌｉｊ≠０である宇宙線の数（即ち、第ｊ番のボクセルに衝突する宇宙線の数）である。そしてＳｉｊ(ｎ)は次で定義される。
【数３９】
【０１１２】
式(36)のλｊについての導関数をゼロと置くと、補助関数を最大にする以下の繰り返し公式を得る（Ｍステップ）。
【数４０】
【０１１３】
Ｓｉｊの二次形式は、λ(ｎ+1)が正であることを保証する。条件付き期待値Ｓｉｊの計算が残っている。Ｘが確率変数Ｈｉｊ|Ｄｊを表すとすると、二次形式ＸＴＡ-1Ｘの期待値は、
【数４１】
【０１１４】
ここで、μｘとΣｘはＸの平均と共分散である。ＤｉはＨｉｊに一次従属しているので、これらは、同時確率ガウス分布である。またＤｉが与えられたときＨｉｊの条件付分布もガウス分布である。これは多変量分布理論からの結果である。この理論と、ＨｉｊとＤｉはそれぞれ平均値ゼロを持つという事実を使うと、以下を得る。
【数４２】
【数４３】
【０１１５】
ここで、観測されるデータの共分散、ΣＤｉ、は式(29)で与えられ、見えないデータ要素の共分散、ΣＨｉｊ、は式(21)を介して表される。観測されるー見えないデータの共分散、ΣＤｉＨｉｊを陽に書くよりも、簡単な（冗長だが）行列計算を行って、次式を示すことができる。
【数４４】
【０１１６】
式(39)から(42)からの結果を式(37)へ代入すると、次を得る。
【数４５】
【０１１７】
ここで、最後のステップで、Ｔｒ(ＡＢ)＝Ｔｒ(ＢＡ)を用いた。
【０１１８】
最後に、ｘとｙ座標の散乱データを組み入れるために、更新の式(38)に於いて、わかり易く次の平均を使用する。
【数４６】
【０１１９】
式(38)は、１つのボクセルに衝突する宇宙線に亙っての平均を表していることを留意する。この式の用法を今後、平均値法と呼ぶ。以下に、この更新式の別の形式が、局外のミューオンデータに起因するノイズの削減に有用であることを示す。アルゴリズムの中央値法は、改正した更新の次式によって定義される。
【数４７】
【０１２０】
上記を考慮して、期待値最大化アルゴリズムを使用する、物体容積体の推定密度プロファイルの尤度を最大化する処理（方法100の処理ステップ１０３）が、図１５Ｅに概説されている。処理ステップ１８０で示すように、ｉ＝1からＭ(△θｘ、△θｙ、△ｘ、△ｙ、ｐｒ2)ｉまでの各ミューオンに対して、散乱と運動量の測定値が集められる。各ミューオンとｊ＝1からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊが推定される（処理ステップ１８１）。ミューオンとボクセルの各対に対して、重み行列：Ｗｉｊが式(24)を使って計算される（処理ステップ１８２）。各ボクセルの散乱密度が、推測値λｊ,ｏｌｄで初期化される（処理ステップ１８３）。処理ステップ１８４に於いて、停止基準が以下のように示されている。各ミューオンに対して、式(29)を使用して
【０１２１】
が計算され、その逆行列が取られる。ミューオンとボクセルの各対に対して、条件付き平均値の項：Ｓｉｊが式(43)と(44)を使用して計算される。式(38)又は(45)を使用して、λｊ,ｎｅｗが計算され、すべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗと置かれる。
【０１２２】
さらに、方法100を説明するために、これから数値例について述べる。図１に示す構成と同じ構成がシミュレーションされた。最初の妥当性確認試験として、多重統計散乱モデルに厳密に適合するように設計された簡素なシミュレーションが使用された。単一の検出器プレーン（図に示してある３枚ではない）のサイズは２×２ｍ2にされた。上部と下部の検出器間の垂直距離は１．１ｍであった。これらの検出器はミューオンの位置と角度を完璧に記録した。運動量が５００−１０，０００ＭｅＶ/ｃに一様に分布するミューオンを持つ、単純化されたミューオンスペクトルが使用された。素粒子は、垂直から一様に広がる投影角度で上部の検出器プレーンの容積体に入射した。
【０１２３】
ミューオンの多重散乱と変位が処理ステップ１１０から１１３に従ってシミュレーションされた。図７と８に視覚化されているように、物体は、容積体の中央１．１×１．１×１．１ｍ3の部分内に置かれた。３個の１０×１０×１０ｍ3の立方体物質、散乱密度がそれぞれ７１．５、１４．２、２．８ｍｒａｄ2/ｃｍであるタングステン(Ｗ)、鉄(Ｆｅ)、アルミニウム(Ａｌ)がシミュレーションされた。シミュレーションは４００，０００個のミューオンを仮定し、これらは上部検出器スタックに入射し、約１０分の露出量に対応する。これらのミューオンの約１６０，０００個が下部検出器プレーンに衝突せず、再構成用には２４０、０００個が残った。再構成には、５×５×５ｍ3のボクセルサイズが使用され、上記で述べた平均値法が実装され、これは、各ミューオンに対する運動量の完全な知識を仮定している。シミュレーションは、空気が充填された容積体から開始し、アルゴリズムを１００回の繰り返しまで走らせた(ブロックの特徴の収束には十分である)。結果は図９に表示されている。３個の物体の各々の８個のボクセルに対して再構成された値の平均値は、Ｗ、Ｆｅ、Ａｌのブロックに対してそれぞれ７４．０、１４．７、２．７である。各立方体を構成する８個のボクセルの変動係数(ｒｍｓ/ｍｅａｎ)は１２．６％、１３．２％、１２．１％である。この結果は、シミュレーションとインバージョンモデルの一致を考えると、インバージョンアルゴリズムと実装を確証している。
【０１２４】
再構成結果は、物体シーンに一致しているように見える。３個の物体の各々の８個のボクセルに対して再構成された値の平均値は、Ｗ、Ｆｅ、Ａｌのブロックに対してそれぞれ７４．０、１４．７、２．７である。各立方体を構成する８個のボクセルの変動係数(ｒｍｓ/ｍｅａｎ)は１２．６％、１３．２％、１２．１％である。この結果は、シミュレーションとインバージョンモデルの一致を考えると、インバージョンアルゴリズムと実装を確証している。
【０１２５】
次に同じシーンをＧＥＡＮＴ４のモンテカルロパッケージを使用してシミュレーションを行った。ＧＥＡＮＴ４の詳細は、J.Allison, “Geant4 developments and applications”, IEEE Trans. Nucl. Sci., vol.53, no.1, pp.270-278, Feb. 2006の刊行物にある。この開示は参照によりここに組み込んだものとする。ＧＥＡＮＴ４は、多重散乱に対する、より完全で、正確でかつ妥当性が実証されたモデルを実装している。このモデルは、散乱分布の中央のガウス分布の部分の幅のより正確な計算、大きいテール部分の実装、物質を通過する時のミューオンのエネルギーロスのシミュレーションを含んでいる。ミューオンイベントジェネレータも使用した。これは、海面レベルの宇宙線ミューオンの角度及び運動量の分布を再現する。このシミュレーションでは、各ミューオンの完全な知識とともに、検出器は完璧であることが仮定された。宇宙線電子や飛跡２次粒子は含められなかった。結果は、図１０に表示されている。Ｗ,Ｆｅ,Ａｌブロックに対応するボクセル値の平均は、それぞれ６７４．４、６３．４、５．４である。
【０１２６】
ボクセル値はあまりにも大きく、中位及び低Ｚの領域のいくつかの誤判別は明らかである。すべてのボクセル値を約４で除して再構成値を正規化することによって、中位のＺのボクセルに対して正しい平均ボクセル値を得られるが、高及び低Ｚのボクセルに対しては正しい値を与えない、又はすべての誤判別の解消はない。この結果の原因は、ガウス分布モデルで旨く記述されていないように散乱する少数のミューオンである。散乱の投影角度分布の中央９８％は、ガウス分布としてよく近似されていると断言される。約２％のミューオンは、ここで説明した統計的モデルに比較して大きい角度で、即ち、ガウス分布の場合に見られるより非常に大きい角度で散乱する。散乱角度の２乗がフィッティングを決定するので、影響は劇的である。これらのテール部分にあるミューオン散乱は、非常に大きい散乱密度推定値を形成する。
【０１２７】
さらに、図１の機器内のミューオンの崩壊などのその他の過程や重大な検出器エラーは、非常に大きな散乱事象として誤って記録される可能性がある（とはいえ、本願に於けるシミュレーションには、これらの発生源は存在しなかった）。これは、容積体のどこでも発生し得るし、不合理に大きい散乱密度を持つ単独のボクセルを作り出す傾向がある。このような事象は、ＳＮＭについて誤った陽性表示を与えるので、除外する必要がある。
【０１２８】
ＥＭアルゴリズムが非ガウス分布データに対して耐性を持つようにするために、平均値更新規則の式(38)は、式(45)、即ち中央値法の使用で置き換えることができる。
【０１２９】
中央値法を使った結果が、図１１に示されている。Ｗ,Ｆｅ,Ａｌ領域に対するボクセル平均値は、それぞれ７９．２、１４．２、２．１であり、変動係数は２１．５％、２６．３％、２３．２％である。明らかに、中央値法を使うとインバージョンアルゴリズムの頑強性が向上する。
【０１３０】
図１２に、再構成される密度プロファイルのより実際的な例が示されており、乗用バンの詳細ＧＥＡＮＴ４シミュレーションを図解している。バン・ボディーが切り取られた主要構成要素のイラストが示されている。イラストの中央の赤いブロックは、１０×１０×１０ｃｍ3のタングステンの固体片、即ち高Ｚの危険物体の代用物を表す。この場合、そのシーンの４個の長い側面上に配置した模擬の検出器プレーンを使用し、より水平方向に向かうミューオンを利用した。宇宙線ミューオンに対する１分の暴露量がシミュレーションされ、平均値法及び中央値法に従い、５×５×５ｃｍ2のサイズのボクセルを使用して、データから再構成が実施された。図１３と１４は、それぞれ、平均値ＥＭ法と中央値ＥＭ法を使用して行った再構成結果を可視化したものを示している。このシーンの平均値法の再構成結果に於ける、非ガウス分布データの影響は、容易に判り、画像上に散らばったより黒いスポットとして表れている。中央値法の再構成結果に於いては、これらのアーチファクトは全く無く、バンのより密度の高い構成要素（エンジン、バッテリ、駆動系）が（低Ｚ）又は(中位のＺ)として表れており、他方、危険物体は、より黒く(高Ｚ)目立っている。中央値法の使用によって、非ガウス分布の散乱分布及び他の異常事象より生じる誤検出に対して堅牢な傾向の結果が得られる。
【０１３１】
ここに示した実施形態と例は、本発明とその実際の応用を最もよく説明し、それによって同業者が本発明を成し利用できるようにするために提示された。しかし同業者は、上記説明及び実施例は、説明と例のためにのみ提示されたことを認識しなければならない。
【０１３２】
本発明の他の変形及び変更（再構成の処理の調整など）は、同業者にとって明らかであり、このような変形と変更がカバーされることは、添付の請求範囲の意図である。
【０１３３】
記載した説明は、包括的であること、又は本発明の範囲を制限することを意図していない。上記教示を考慮すると、以下に続く請求範囲から逸脱することなく、多くの変更及び変形が可能である。本発明の使用には、様々な特徴を有する構成要素が含まれると考えられる。
【技術分野】
【０００１】
本出願は、(1)「素粒子の検出と解析のためのシステム、方法及び装置、並びにフィールドへの展開」(SYSTEM, METHODS AND APPARATUS FOR PARTICLE DETECTION AND ANALYSIS AND DEPLOYMENT OF THE SAME)という名称の、２００６年１０月２７日出願の米国特許仮出願第６０/８５５、０６４号と、(2)「放射物ポータルモニターシステム及び方法」(RADIATION PORTAL MONITOR SYSTEM AND METHOD)という名称の、２００７年６月２９日出願の米国特許出願第１１/７７１、１６９号に付属する優先権を主張する。
【０００２】
上記の２つの出願の開示は、参照によりここに組み込んだものとする。
【０００３】
（米国連邦政府の権利に関する宣言）
本発明は、米国エネルギー省によって与えられた、契約番号DE-AC52-06NA25396に基づく政府支援でなされた。米国政府は、本発明に於いて一定の権利を保有する。
【０００４】
実施形態は、粒子の検出、解析、制御の分野に関し、特に、以下に限定されるわけではないが、複数の位置感知検出器を有する荷電粒子検出システムのデータを解析し、そして荷電粒子検出システムを通過する宇宙線ミューオンなどの荷電粒子の軌道を再構成する方法及び装置に関する。
【背景技術】
【０００５】
荷電粒子トモグラフィは、荷電粒子の散乱に基礎を置いている。荷電粒子トモグラフィの一形態は、宇宙線トモグラフィであり、宇宙線ミューオンの散乱を利用する。深宇宙からやって来る、安定な粒子、主にプロトンが、絶えず地球に降り注いでいる。これらの粒子は、上層大気圏の原子と相互作用して多数の短寿命のπオンを含む粒子シャワーを作り出し、それらπオンは崩壊し長寿命のミューオンを生み出す。ミューオンは、主にクーロン力により物質と相互作用し、核相互作用はせず、電子に比べると直ちに放射状に広がることはほとんどない。ミューオンは、電磁相互作用によって、ゆっくりとしかエネルギーを失わない。従って、ミューオンの多くは、透過度が高い荷電放射線として地球表面に到達する。海抜ゼロでのミューオン・フラックスは、約１ミューオン/ｃｍ2・ｍｉｎである。
【０００６】
ミューオンが物質を通過する時、原子構成要素の電荷とのクーロン散乱が、ミューオンの軌道を乱す。トータルの偏向は、いくつかの物質特性に依存するが、支配的なパラメータは、原子核の原子番号、Ｚ、と物質密度である。
【発明の概要】
【発明が解決しようとする課題】
【０００７】
注目する容積体を通過するミューオン又は他の荷電粒子からその容積体を再構成するための改良された方法及び装置を提供するニーズがある。
【課題を解決するための手段】
【０００８】
以下の発明の概要は、ミューオンのような荷電粒子などの粒子を検出し、荷電粒子トモグラフィ・データから物体容積体散乱密度プロファイルを統計的に再構成するための技術、装置、及びシステムに関するいくつかの技術的特徴の理解を容易にするために設けており、完全な説明を意図してはいない。本発明の様々な態様の完全な理解は、明細書全体、請求範囲、図面及び要約を全体として捉えることによって得ることができる。
【０００９】
本発明の上記の態様、他の目的及び利点は、ここに述べるように達成することができる。
【００１０】
一態様による、物体容積体を通過する荷電粒子を介して物体容積体を検出する検出システムが記載されている。このシステムは、物体容積体の第１の側に配置され、物体容積体に向かう入射荷電粒子の位置と角度を測定する、第１の位置感知検出器セットと、物体容積体の前記第１の側の反対側の第２の側に配置され、物体容積体から出て行く出射荷電粒子の位置と角度を測定する、第２の位置感知検出器セットと、前記第１の位置感知検出器セットからの測定信号のデータと前記第２の位置感知検出器セットからの測定信号のデータを受け取る信号処理ユニットとを含む。信号処理ユニットは、受信データを処理し、物体容積体内の容積体散乱密度分布を統計的に再構成する。
【００１１】
信号処理ユニットは、(ａ)物体容積体を通過する荷電粒子の散乱角と推定運動量に対応する荷電粒子トモグラフィ・データを取得し、(ｂ)統計的多重散乱モデルに基づいて、荷電粒子散乱密度の確率分布を用意し、(ｃ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて物体容積体散乱密度の略最尤推定値を決定し、(ｄ)前記略最尤推定値に基づいて再構成された物体容積体散乱密度を出力するように構成することができる。
【００１２】
別の態様によると、物体容積体から得られる荷電粒子トモグラフィ・データから物体容積体を検出する方法であって、(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、(ｂ)統計的多重散乱モデルに基づいて、荷電粒子散乱の確率分布を用意する工程と、(ｃ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）用いて、物体容積体散乱密度の略最尤推定値を決定する工程と、(ｄ)再構成された物体容積体散乱密度を出力する工程と、そして必要なら、(ｅ)再構成された物体容積体散乱密度に基づいて判定を行う工程と、を含む。
【００１３】
この方法によって、ユーザは、前記再構成された物体容積体散乱密度プロファイルより、注目する容積体を占める物体の存在と種類の両方又はいずれか一方を特定することが可能となる。いろいろな応用には、様々な自国保安の検査用途の宇宙線ミューオン・トモグラフィが含まれ、その用途では、車両又は貨物がミューオン追跡装置によって、走査され得る。
【００１４】
荷電粒子トモグラフィ・データは、宇宙線または他のなんらかの線源によって発生したミューオンなどの荷電粒子から収集されたトモグラフィ・データを含んでよい。
【００１５】
前記再構成された物体容積体散乱密度に基づいて判定を行う工程には、再構成された物体容積体散乱密度に基づいて、物体容積体内の目標物体の(1)存在と(2)種類の内のすくなくとも１つについて判定を行う工程を含めてよい。
【００１６】
期待値最大化アルゴリズム（ＭＬ/ＥＭ）で使用される、荷電粒子散乱の確率密度を用意する工程には、(ｇ)均一物体の所定の散乱密度に基づいて、荷電粒子の２次元確率分布を取得する工程と、(ｈ)前記２次元確率分布に基づいて、前記荷電粒子の３次元確率分布を取得する工程と、(ｉ)基底関数によって特徴付けられる非均一物体容積体による、多数の荷電粒子の散乱の確率分布を取得する工程と、(ｊ)その定義、及び荷電粒子の散乱と運動量の測定値に基づいて、多重散乱の確率分布を決定する工程と、を含めてよい。
【００１７】
均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元確率分布を取得する工程は、(ｋ)物質の散乱密度を、単位深さの物質による荷電粒子の散乱の２乗の期待平均値として決定する工程と、(ｌ)荷電粒子散乱分布をガウス分布モデルにより近似する工程と、(ｍ)荷電粒子線散乱及び変位を相関２次元ガウス分布により近似する工程と、を含んでよい。
【００１８】
２次元確率分布に基づいて、荷電粒子の３次元確率分布を得る工程は、座標を追加し３次元経路長を規定する工程と、３次元変位を計算する工程と、３次元共分散行列を規定する工程と、を含んでよい。
【００１９】
基底関数によって特徴付けられる非均一物体容積体による、多数の荷電粒子の散乱の確率分布を取得する工程は、離散的な散乱密度推定値を表す基底関数の３次元グリッドを設定する工程と、各個々のミューオンに対する散乱/変位の共分散行列を、宇宙線経路と散乱密度推定値の関数として決定する工程と、複数の荷電粒子に対する確率分布を個々の荷電粒子の確率の積として決定する工程と、を含んでよい。
【００２０】
前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の前記略最尤推定値を決定する工程は、各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、各荷電粒子と前記統計的散乱モデルの各基底関数との相互作用の幾何学を推定する工程と、荷電粒子と基底関数との各対に対して、重み行列：Ｗｉｊを決定する工程と、各基底関数の散乱密度を推測値で初期化する工程と、物体容積体内容に関して、近似の最尤解を反復して求める工程と、を含んでよい。ここで反復処理は、所定の反復回数で、又は解が所定の許容値未満に変わる時、停止する。
【００２１】
前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の略最尤推定値を決定する工程は、ｉ＝１からＭ(△θｘ、△θｙ、△ｘ、△ｙ、ｐｒ２)ｉまでの各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、各ミューオンとｊ＝１からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊを推定する工程と、荷電粒子とボクセルの各対に対して、
【数１】
【００２２】
として、重み行列：Ｗｉｊを計算する工程と、各ボクセルに於ける散乱密度λｊ,ｏｌｄの推測値を初期化する工程と、停止基準処理を使用してすべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗと置く工程と、を含む。
【００２３】
前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）は、平均値更新規則又は中央値更新規則を含んでよい。
【００２４】
さらなる別の態様によると、物体容積体から得られる荷電粒子トモグラフィ・データから物体容積体を検出する、コンピュータによって実行される方法であって、(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、(ｂ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）に使用される、統計的多重散乱モデルに基づく、荷電粒子散乱密度の確率分布を用意する工程と、(ｃ)前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の略最尤推定値を決定する工程と、(ｄ)再構成された物体容積体散乱密度を出力する工程と、を含む。判定は、前記再構成された物体容積体散乱密度に基づいて行うことができる。
【００２５】
さらなる別の態様によると、コンピュータ・プログラム製品は、コンピュータアクセス可能なデータ記憶媒体を含み、データ記憶媒体は、コンピュータによって実行される時荷電粒子トモグラフィ・データから物体容積体密度プロファイルを統計的に再構成する方法をコンピュータに実施させる命令を格納する。前記方法は、(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、(ｂ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）に使用される、統計的多重散乱モデルに基づく、荷電粒子散乱の確率分布を用意する工程と、(ｃ)前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体密度の略最尤推定値を決定する工程と、(ｄ)再構成された物体容積体散乱密度を出力する工程と、を含む。
【００２６】
添付図面於いては、同じ参照符号は、すべての図面を通じて、全く同じ又は機能的に同じ要素を参照している。添付図面は、組み込まれて明細書の一部を形成するが、それ以上に、本発明を図解し、本発明の詳細説明とともに、本発明の原理を説明する役割を果たす。
【図面の簡単な説明】
【００２７】
【図１】ミューオン・トモグラフィ概念の一例を図解している。
【図２】クーロン散乱を確定するのに使用される散乱と変位の２次元投影を図解している。
【図３】３次元散乱モデルを調整するのに使用される、散乱と変位の２次元投影のパラメータを図解している。
【図４】物質の多層レイヤによる散乱を図解している。
【図５】図３に示す投影の経路長計算に対する再近接点の使用を図解している。
【図６】ミューオン・トモグラフィの容積体散乱密度プロファイルを統計的に再構成する自動システムを図示している。
【図７】シミュレーションされた物体の斜視図である。
【図８】シミュレーションされた物体の上からみた図である。
【図９】ＰｏＣＡ点で接続される２個の直線を仮定して算定した経路長によるガウス分布散乱シミュレーションの再構成結果を図示している。
【図１０】非ガウス分布のテール部分の散乱によるシミュレーションデータの再構成結果を図示している。
【図１１】中央値法よる非ガウス分布のテール部分の散乱によるシミュレーションデータの再構成結果を図示している。
【図１２】シミュレーションされる乗用バンの主要物体を図解している。
【図１３】平均値法による乗用バンのミューオン暴露１分のシミュレーションの再構成結果を図示している。
【図１４】中央値法による乗用バンのシーンの再構成結果を図示している。
【図１５】一実施形態による、荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法の概略を説明しているフローチャートである。
【図１５Ａ】一実施形態による、多重・統計的散乱モデルを使用し、物体容積体を通過するミューオンの散乱の確率分布を推定する処理の一例の概略を説明しているフローチャートである。
【図１５Ｂ】一実施形態による、物体の所定散乱密度に基づいて、単一ミューオンに対する散乱の２次元期待確率分布を推定する処理の一例の概略を説明しているフローチャートである。
【図１５Ｃ】一実施形態による,統計的モデルの３次元への拡張する処理の一例の概略を説明しているフローチャートである。
【図１５Ｄ】一実施形態による,非均一物質による多数のミューオンの散乱及び変位の確率分布を決定する処理の一例の概略を説明しているフローチャートである。
【図１５Ｅ】一実施形態による,期待値最大化アルゴリズムを使用して物体容積体の推定密度プロファイルの尤度を最大化する処理の一例の概略を説明しているフローチャートである。
【発明を実施するための形態】
【００２８】
以下の非限定な例に於ける特定の値及び構成は、変更することができ、本発明の少なくとも一実施形態を説明するためにのみ出されており、発明の範囲を制限することを意図していない。
【００２９】
本願で説明している技術的特徴は、様々な粒子検出システムの構築に使用することができる。例えば、ミューオンを荷電粒子として検出する粒子検出システムは、検査対象の物体が配置される物体収容領域と、物体収容領域の第１の側に配置され、物体収容領域へ向かう入射ミューオンの位置及び角度を測定する、第１の位置感知検出器セットと、物体収容領域の第１の側の反対側の第２の側に配置され、物体収容領域から出ていく出射ミューオンの位置及び角度を測定する、第２の位置感知検出器セットと、第１の位置感知検出器セットからの入射ミューオンの測定信号のデータと第２の位置感知検出器セットからの出射ミューオンの測定信号のデータを受け取る、例えばマイクロプロセッサを備えてよい信号処理ユニットと、とを含むことができる。一例として、粒子検出器の第１及び第２のセットのそれぞれは、第１の方向の少なくとも３つの荷電粒子位置測定及び第１の方向と異なる第２の方向の少なくとも３つの荷電粒子位置測定が可能となるように配置されたドリフトチューブを含むように実現できる。信号処理ユニットは、ミューオンの測定された入射及び出射の位置及び角度に基づいて、物体収容領域内の物質におけるミューオン散乱によって生じるミューオンの散乱挙動を解析し、トモグラフィ・プロファイル又は物体収容領域内の散乱中心の空間分布を得るように構成される。得られたトモグラフィ・プロファイル又は散乱中心の空間分布は、物体収容領域における１つ以上の物体、例えば、核物質を含む大きな原子番号の物質や機器などの存否を明らかにするのに使用することができる。各位置感知ミューオン検出器は、ドリフトセル、例えば、ミューオンによってイオン化できるガスで充填されたドリフトチューブを含む様々な形状で実現できる。このようなシステムを使って、自然界の宇宙線ミューオンを、物体収容領域内の１つ以上の物体の検出用のミューオン線源として利用することができる。
【００３０】
一実施例においては、例示的実施形態による、荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法及びシステムは、物体を通過する宇宙線荷電粒子の散乱に基づいて、物体のイメージ又はモデルを再構成することが可能な解決法を提供する。
【００３１】
軌道は、よりありふれた物体（水、プラスチック、アルミニウム、スチールなど）を構成する物質より、特定の核物質（ＳＮＭ）及び良好なガンマ線遮蔽となる物質（鉛やタングステンなど）によって強く影響される。宇宙線荷電粒子、特に宇宙線ミューオンに関して、各ミューンは、それが透過した物体についての情報を運んでおり、多数のミューオンの散乱を測定することによって、これらの物体の特性を調べることができる。特に、より一般的な低Ｚや中間Ｚの物質の中の高Ｚの物体を検出することができる。
【００３２】
例示的実施形態による、荷電粒子トモグラフィの容積体密度プロファイルの統計的な再構成法の様々な技術的特徴を説明するために、先ず、図１にその例が図解してあるミューオントモグラフィ・コンセプトについて述べる。
【００３３】
位置感知検出器のセット10は、撮像対象の物体容積体の上方と下方に配置され、入射及び出射荷電粒子の飛跡12（矢印付き実線で示されている）の位置と角度を提供する。撮像対象の容積体上方に配置されている位置感知検出器10の２個以上のセットが、入射荷電粒子の飛跡の位置と角度を与える。これらの検出器は、２つの直交又は非直交座標に於ける荷電粒子の位置を測定する。位置感知検出器10の別のセットは出射荷電粒子の位置と角度を記録する。側面の検出器（図示せず）を使って、より水平に向いた荷電粒子飛跡を検出してもよい。各荷電粒子飛跡の散乱角は、符合する入射及び出射測定値より算出される。荷電粒子運動量は、検出器自身内で起こる僅かな散乱より算定される、又は位置感知検出器プレーンの２個のセット間に置かれた、既知の特性の散乱体のレイヤ内で起こる僅かな散乱より算定される。
【００３４】
位置感知荷電粒子検出器の一例は、作動ガスが充填されたドリフトチューブである。ドリフトチューブは、円筒状の管でよく、ミューオンなどの宇宙線荷電粒子の検出を可能とする、アルゴン・イソブタンなどの検出器ガスが充填される。管の外部表面をアースし、円筒状管の長さ方向に沿って延びている中央の陽極ワイヤに約２〜３ｋＶの正の高電圧が印加され、高電圧の静電場が生じる。荷電粒子がガス原子達と相互作用すると、管の弦(ｃｈｏｒｄ)によって一直線上の原子達から多くの電子が解放される。電子達の「ひも」は、静電場によって、正に帯電した陽極ワイヤへ向かってドリフトし、データ取得電子機器のＴＤＣＳ（時間デジタル変換器）で電子的に読み出される。検出器の各セットは、第１の方向の少なくとも３つの荷電粒子位置測定と、第１の方向とは異なり、第１の方向と直交してよい第２の方向の少なくとも３つの荷電粒子位置測定とを可能とするように配置されたドリフトチューブ達であってよい。
【００３５】
図１のシステムには、信号処理ユニット、例えばコンピュータが備えられ、物体容積体の上方の検出器による入射ミューオンの測定信号データと、物体容積体の下方の検出器による出射ミューオンの測定信号データとが受け取られる。この信号処理ユニットは、測定されたミューオンの入射及び出射の位置並びに角度に基づいて、容積体内の散乱によるミューオンの散乱挙動を解析し、トモグラフィ・プロファイル又は容積体内の散乱中心の空間分布を取得するように構成される。取得されたトモグラフィ・プロファイル又は容積体内の散乱中心の空間分布によって、容積体内の物体の存否を明らかにすることができる。ある実施例に於いては、追加のドリフトチューブ検出器を容積体の側面に配置して、システムの走査のために、パッケージ、乗物又は貨物コンテナを入れることができるボックス又は四面構造を形成することができる。このようにして、宇宙線ミューオンの多重散乱を使用して、通常の貨物のバックグラウンド内の高Ｚ物質(高原子番号)を選択的に検出することができる。利点として、この技法は、受動的であり、バックグラウンドを超える如何なる放射線量も加えず、高Ｚ高密度物質について選択的である。信号処理ユニットのトモグラフィ処理部は、検出器10と同じ場所にある現場のコンピュータに組み込んでもよい。あるいは、信号処理ユニットのトモグラフィ処理部は、プライベート・ネットワークのようなコンピュータ・ネットワーク、又はインターネットのような公共ネットワークに接続されているリモート・コンピュータに組み込んでもよい。
【００３６】
図１の例示的実施形態においては、荷電粒子は、宇宙線ミューオン又は他の宇宙線荷電粒子であり、位置感知検出器10は、荷電粒子を感知するための作動ガスが充填されたドリフトセルである。ドリフトセルは、例えば、中央の陽極ワイヤが管の長手方向に沿って走っているドリフトチューブで実現してよい。しかし、ミューオン以外の荷電粒子は、ドリフトセル以外の位置感知センサーを用いて検出してもよい。さらに荷電粒子は、宇宙線以外の線源によって発生してもよい。例えば、ミューオンは、加速器から低強度のビームとして発生することができる。
【００３７】
高密度物体を透過するミューオン（ブラックの飛跡）は、大気中を走るミューオン（グレーの飛跡）よりかなり強く散乱する。多数の飛跡測定より、物体形状と物質の電子密度を再構成することができる。容積体を通過するミューオンは、通過する物質に依存して散乱される。
【００３８】
図１のシステム用処理ユニットによる、検査中の容積体（例えば、パッケージ、コンテナ又は車両）に於ける宇宙線ミューオンの測定値の処理には、容積体を通るミューオンの軌道を再構成する工程と、容積体の各側方にある検出器からの信号を基に、入射ミューオンの運動量を測定する工程と、容積体の散乱密度の空間分布を決定する工程と、を含めてよい。これら及び他の処理結果を使用して、トモグラフィ・プロファイルを構築し、目標物体の検出などのような容積体の様々な特性を測定することができる。
【００３９】
例えば、ドリフトチューブ一式を有する検出器10を通過する荷電粒子の軌道の再構成には、(ａ)荷電粒子が衝突したドリフトセルの識別子と対応する衝突時刻を表す衝突信号を取得する工程と、(ｂ)検出器を通過する特定の荷電粒子の飛跡に対応していると同定された時間内のドリフトセル衝突をグループ化する工程と、(ｃ)前記特定の荷電粒子がドリフトセルに衝突する瞬間の時刻ゼロ値を最初に推定する工程と、(ｄ)時刻ゼロ値の推定値、ドリフト時間変換データ及び衝突時刻に基づいてドリフト半径達を決定する工程と、(ｅ)特定の時刻ゼロ値に対応するドリフト半径達に直線の飛跡をフィティングする工程と、(ｆ)特定の荷電粒子に対して行われたもっとも良い飛跡フィットに関連する時刻ゼロ値を検索して選択し、時刻ゼロと追跡パラメータに於ける誤差を算出する工程と、を含めてよい。時刻ゼロ・フィティングによるこのような飛跡の再構成法に於いては、高速な検出器（シンチレータ・パドル付光電子増倍管など）を、或いはミューオンの装置内の通過を数ナノ秒の単位で検出して時刻ゼロ値を与える何か他の高速検出器を使用することなく、荷電粒子検出器を通過する荷電粒子の再構成された直線軌道が得られる。
【００４０】
又、例えば、図１の検出器10からの信号に基づいて、入射ミューオン又は出射ミューオンの運動量を測定する処理には、(ａ)複数の位置感知検出器を配置し、通過する荷電粒子を散乱させる工程と、(ｂ)位置感知検出器内の荷電粒子の散乱を測定する工程と、(ｃ)位置測定値から少なくとも１つの荷電粒子の軌道を決定する工程と、(ｄ)前記少なくとも１つの軌道から荷電粒子のすくなくとも１つの運動量測定値を決定する工程と、を含めてよく、前記散乱の測定は少なくとも散乱荷電粒子の３つの位置測定値を取得する工程を含む。この技法を使用して、検出器内に追加のメタルプレートを使用することなく、位置感知検出器自身内での荷電粒子の散乱より決定される荷電粒子の軌道に基づいて、荷電粒子の運動量を決定することができる。
【００４１】
以下に、荷電粒子トモグラフィ・データからの物体容積体散乱密度プロファイルを統計的に再構成する例示的システム及び方法の詳細を示している。
【００４２】
一実施形態による、荷電粒子トモグラフィの物体容積体散乱密度プロファイルを統計的に再構成する自動化されたシステムの一例が、図６のブロック図に示してある。自動システム50は、荷電粒子トモグラフィ・データ54を受け取るように適合され構成されたコントローラ51を有している。荷電粒子トモグラフィ・データは、例えば、図１の荷電粒子検出器、或いは、容積体を通過する荷電粒子の追跡記録ができるように構成された位置感知検出器を有する他の荷電粒子検出器を使ったミューオンの測定値より決定されたミューオンのトモグラフィ・データでよい。結果として、ミューオン又は他の荷電粒子のトモグラフィ・データを使って、物体容積体を通過するミューオン又は他の荷電粒子の散乱角及び推定運動量を得る又は決定することができる。
【００４３】
自動システム50は、コントローラ上に格納された統計的再構成モジュール56を含む。再構成モジュール56は、ミューオン又は他の荷電粒子のトモグラフィの容積体散乱密度プロファイルの統計的再構成を司る。再構成モジュール56は、ソフトウエアモジュール又はハードウエアモジュールとして実現してよい。
【００４４】
図６の自動システム50の例示的実施形態に於いて、コントローラ51は、コンピュータ(ＰＣ)のような１つ以上の動作可能に連結された計算プロセッサユニット(ＣＰＵ)ベースのシステム、又はデジタル信号プロセッサベースのシステムのような他のマイクロプロセッサベースのシステムを使って、構成される。コントローラは、単一の通常のコンピュータでもよいが、しかし、実時間的に結果を得るために、コントローラは、典型的には、実時間で結果を得るのに必要な処理能力を提供するように、十分な数の並列処理コンピュータ（図示せず）という形態を含む。例えば、コントローラは、約２０個のＣＰＵを含んでよい。ミューオン検出器の走査容積体が大きくなるほど、そして所望の分解能が細かくなるほど、処理コンピュータの形態は大きくなる必要がある。
【００４５】
オペレーティングシステムは、コントローラ51上で動作し、そして、限定するものではないが、Ａｐｐｌｅ，Ｗｉｎｄｏｗｓ(登録商標)，Ｌｉｎｕｘ，又はＵｎｉｘ(登録商標)のオペレーティングシステム、又は将来開発される他のオペレーティングシステムを含め、市販又はオープン・ソースのオペレーティングシステムでよい。オペレーティングシステム、及びアプリケーション又はプログラムの命令は、ハードディスク・ドライブなどの記憶装置に格納される。また、自動システム50に於いては、飛跡の再構成モジュールは、コンピュータアクセス可能なデータ記憶媒体の形式のソフトウエアであって、データ記憶媒体は、コントローラによって実行される時例示的実施形態による荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法をコントローラに行わせしめる命令を格納している。モジュールは、図６に示されるように、現場のコントローラにインストールしてもよいし、コントローラに接続されたインターネット経由で遠隔地から起動してもよい。同業者なら、このようなモジュールの実現形態が複数あることは理解する。
【００４６】
自動システム50は、又、操作可能にコントローラ51に接続されたディスプレイ58を含む。ディスプレイ58はシステムによって再構成された物体密度プロファイルの画像又はデータを必要に応じて、ユーザに表示する。ユーザインタフェース（図示せず）が処理システムに操作可能に接続され、人のオペレータが必要に応じて、処理システムを操作できるようにしてもよい。
【００４７】
同業者であれば、図６の説明は、単に自動システム50の実施形態の一例を描いているにすぎず、実施形態はそれに限定されないことが理解できる。例えば、再構成モジュールの機能の一部又はすべては、マイクロプロセッサを使うことなく、アナログ又はデジタル回路などのハードウエアとして実現できる。
【００４８】
図１５を説明する。これは、一実施形態による、荷電粒子トモグラフィの容積体散乱密度プロファイルを統計的に再構成する方法の概要を説明するフローチャートを示している。方法100は、処理ステップ１０１に示されているように、物体容積体を通過する荷電粒子の飛跡位置、散乱角、推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得することで始まる。所定の荷電粒子トモグラフィ・データは、例えば、図１の検出器から取得できる。その後、処理ステップ１０２に示されているように、散乱密度の空間分布（以降に定義される）によって表される、物体容積体を通過する複数の荷電粒子の散乱の確率分布が、多重統計散乱モデルを基づいて、用意される。次に、処理ステップ１０３に示すように、期待値最大化アルゴリズムを使用して、物体容積体散乱密度プロファイルの最尤推定値が決定される。それから、再構成された容積体散乱密度プロファイルが判定用に出力される（処理ステップ１０４）。判定処理は、任意であり、処理ステップ１０５に示すように容積体を占有している物体の存在と種類の両方又はいずれか一方を特定する処理であってよい。判定処理は、物体容積体の再構成された密度プロファイルを表す画像の人による解釈と、さらなるアルゴリズムによる自動的な判定との両方又はそれらの一方を含んでよい。
【００４９】
本実施形態の方法と自動システムによって、多数の荷電粒子によって与えられるデータに基づいて、注目の容積体の離散的なトモグラフィ再構成が可能となる。反復的な期待値最大化(ＥＭ)アルゴリズムのインスタンスを使用して、物体の密度プロファイルの最尤推定値を求める。本実施形態の方法とシステムによって、ユーザは、再構成された容積体密度プロファイルから、注目の容積体を占める物体の存在と種類の両方又はそれらの一方を特定することが可能となる。いろいろな応用には、様々な自国保安の検査用途の宇宙線ミューオン・トモグラフィが含まれ、その用途では、車両又は貨物がミューオン追跡装置によって、走査される。得られるミューオンのトモグラフィ・データは、危険物を識別できる本例示的実施形態の方法と自動システムを用い、車両又は貨物の密度プロファイルを再構成し表示するのに使用できる。
【００５０】
最尤法は、医療画像再構成に於いて、特にＰＥＴ及びＳＰＥＴ再構成のために使用されているが、いくつかの重要な差異のために、それらの用途に開発された標準的な方法は使用できない。第１に、測定信号−散乱角−は確率的であり、平均値ゼロで、標準偏差は透過物質の特性によって規定される。第２に、宇宙線ミューオンは、確定した離散的な方向又は角度から来ることはなく、むしろ天頂の周りに幅広い、連続的な角度分布を持ち、その分布は殆ど水平にまで広がっている。最後に、ミューオンの軌道は、直線ではない。我々が、強い散乱をおこす物体のおおよその場所を見つけることができるのは、その湾曲である。ＥＭアルゴリズムは、柔軟性があり、計算効率が良い。そのアルゴリズムの複雑な幾何形状への適用について説明することができる。
【００５１】
これから、一実施形態に従って、処理ステップ１０２から１０４について説明する。データは、容積体を通過するミューオンを測定する図１の検出器から取得される宇宙線ミューオンのトモグラフィ・データである。
【００５２】
一実施形態による、多重統計散乱モデルを使用して物体容積体を通過するミューオンの散乱の推定確率分布を用意する処理（処理ステップ１０２）が、図１５Ａのフローチャートに概説されている。処理ステップ１１０から１１３で示されているように、本処理は、４つの主な構成部分を有している。最初に、均一物体の所定の散乱密度に基づいて単一ミューオンに対する２次元の確率分布が推定される（処理ステップ１１０）。そして２次元分布は３次元に拡張される（処理ステップ１１１）。次に、処理ステップ１１２に於いて、非均一物体容積体がボクセル基底関数を使用して表現され、ボクセル化された散乱密度を所与として、多数のミューオンの散乱の確率分布が表現される。最後に、確率分布表現は、ミューオン散乱及び運動量の測定値の有限精度に拡張される（処理ステップ１１３）。
【００５３】
処理ステップ１１０から１１３は、多重散乱統計モデルを使って実現される。そのモデルは、先ず均一物質の単一レイヤに於ける散乱に関して、そして非均一物質に於ける散乱に関して説明される。
【００５４】
物質を通過する宇宙線ミューオンは、図２に図解されているように多重クーロン散乱を受ける。図２は、多重クーロン散乱を記述するために使用される散乱と変位の２次元投影を図示している。この図及び他の図では、散乱の大きさは説明のために大幅に誇張されている。出射ミューオンは、入射ミューオンの方位と位置に対して測られた、散乱角と変位によって特徴付けることができる。典型的な散乱角は、数十ミリラジアン（１ミリラジアン≒０．０６度）であり、数度より大きい散乱角は非常にまれである。散乱角の中央９８%の分布は、平均値ゼロのガウス分布として近似してよい。
【数２】
【００５５】
実際の分布は、ガウス分布より重い或いは大きなテール部分をもっている。分布の幅は、物質特性の関数として近似的に表わすことができる。S.Eidelman et al.,"Review of particle physics" Phy.Lett.,vol.B592,p.1,2004にレビューされているように、多くの研究者が、散乱に対して、様々な物質の特性の関数として実験的に改良した式を提案している。その開示は、参照によりここに組み込んだものとする。特に単純な形は
【数３】
【００５６】
である。
【００５７】
ここで、ｐは単位がＭｅＶ/ｃの運動量、Ｈは物質の深度(ｄｅｐｔｈ)、Ｌｒａｄは物質の放射長、βｃは速度（ｃは光の速度）であり、β＝１の近似が使用される。放射長は、原子番号及び物質密度が大きくなるに従い減少する。公称のミューオンの運動量をｐ0とし、放射長Ｌｒａｄの物質の散乱密度を次式で定義する。
【数４】
【００５８】
物質の散乱密度、λ、は、その物質の単位深度を通過する公称運動量を持つミューオンの散乱角の２乗平均値を表す。いくつかの物質のミリラジアン2/センチメートル単位での散乱密度の値は、例えば、アルミニウムでは約３、鉄で１４、ウラニウムで７８である。それで、散乱密度λ、深度Ｈの物質を通過する運動量ｐのミューオンの散乱の分散は、
【数５】
【００５９】
である。ここで、
【数６】
【００６０】
と置くと、
【数７】
【００６１】
を得る。
【００６２】
変位△ｘは、散乱角△θと相関性がある。一緒に扱うと、散乱角と変位は、図１の経路の“屈曲(ｋｉｎｋ)”によって示唆されているように、大きな容積体中の局所散乱の寄与物の位置を示唆する情報を提供する。散乱角と変位の分布は、平均値がゼロで、
【数８】
【数９】
【００６３】
を持つ同時確率ガウス分布として特徴付けることができる。その共分散行列は次式で表すことができる。
【数１０】
【００６４】
ここで、
【数１１】
【００６５】
と置くと、
【数１２】
【００６６】
を得る。
【００６７】
上記を考慮して、一実施形態による、単一のミューオン散乱の２次元に於ける散乱確率分布の取得（処理ステップ１１０）は、図１５Ｂのフローチャートの概略のように説明できる。処理ステップ１５０に示されているように、物質の散乱密度は、式（３）によって、単位深度の物質によるｐ0＝３ＧｅＶ/ｃのミューオンの散乱の２乗の期待平均値として定義される。次に、処理ステップ１５１に示されているようにＲＭＳ散乱に対してガウス分布近似、式(１)、（５）,（６）が行われる。最後に、宇宙線の散乱及び変位の分布は、処理ステップ１５２に示されているように、式(10),(11)に要約されている平均値ゼロの相関２次元ガウス分布によって近似される。
【００６８】
３次元に於いては、散乱は、ｘ座標に直交するｙ座標を考慮することによって特徴付けられる。散乱角△θｘと△θｙ、そして変位△ｘと△ｙを参照されたい。ｘ平面とｙ平面への偏向は互いに独立で同一の分布に従う（Eidelman et alを参照）。前記展開は、入射ミューオンの方向に直交するように置かれた座標系に基づいている。３次元モデルにおいては、経路長を考慮し、変位測定値を入射ミューオン経路に直交する平面に投影する必要がある。図３は３次元散乱に対してモデルを適合させるのに使用されるパラメータを図示しているが、垂直方向に対して、投影角度θｘ,0で入射するミューオンが示されている。
【００６９】
この３次元散乱の理解の助けとなるように、図のページから外に向く直交するｙ座標における対応する投影角度θｙ,0を想像することは有用である。レイヤを通って投影（不散乱）点（ｘｐ、ｙｐ）までのミューオン経路の直線延長部（即ち、３次元経路長）は、
【数１３】
【００７０】
である。
【００７１】
出射ミューオンのｘ位置と角度を（ｘ１、θｘ,1）と定義し、そして
【数１４】
【００７２】
と置く。
【００７３】
測定されるｘ変位は、ｘｍ＝ｘ1−ｘｐとして計算されるであろう。しかし、この測定値を回転させて宇宙線経路に直交する平面に入れ、３次元経路長に対して調整する必要がある。変位は、次で定義される。
【数１５】
【００７４】
ここで、中央の２項は３次元経路長に相当し、最後の項は、測定値を適切な方向に投影するものである。
【００７５】
最後に、共分散重み行列を次のように再定義する。
【数１６】
【００７６】
そして散乱と変位に対して同様に進める必要があり、式(11)は、両方かつ座標散乱に対する共分散行列を定義している。散乱測定は、２つの直交、水平座標で独立に行われる。表記法を簡単にするために、ただ1つの座標について解析を展開する。２つの座標からの情報を結合するのは後で説明する。本モデルは、“小”散乱角と変位に対して有効であることに留意すべきである。大きな角度の散乱では、モデルの導出に於いて無視された２次の項が有意となり得る。
【００７７】
上記を考慮すると、一実施形態による図１５Ｃのフローチャートに於いて概説されているように、３次元に拡張された統計的モデルが得られる（処理ステップ１１１）。先ず、ｙ座標が追加され、３次元経路長が規定される（処理ステップ１６０、式(12）。次に、処理ステップ１６１に於いて、式(13),(14)によって３次元変位が計算される。最後に、式(15)によって、３次元共分散行列を表す（処理ステップ１６２）。
【００７８】
物質の非均一容積体に対しては、密度プロファイルは、再構成のために、係数{ν1、・・・νｊ、νＮ}と３次元の基底関数{φ1、・・・φｊ、ΦＮ}の線形結合として表現され、即ち；
【数１７】
【００７９】
である。
【００８０】
基底関数に対しては、多くの選択肢があるが、ここでは、矩形の３次元ボクセルに注意を払う。λｊは、第ｊ番目の基底関数の係数、即ち、第ｊ番のボクセルの散乱密度を示すのに使用される。図４を考えると、３つのレイヤ（又はボクセル）が示されており、宇宙線がスタックを通過し、観測情報△θと△ｘを届けている。第ｊ番のボクセルの“見えない”散乱及び変位がそれぞれ△θｊと△ｘｊとして示されている。又、図では、散乱の大きさは誇張されている。以下の式によって、観測されるデータを見えないデータに関係付けることができる。
【数１８】
【数１９】
【００８１】
ここで、２番目の式では、小角散乱という仮定に依拠しており、Ｔｊは、第ｊ番のボクセルの出口点から再構成容積体の出口点までの３次元宇宙線経路長として定義している。より一般的には、ボクセル集合χを通過する宇宙線に対して、
【数２０】
【数２１】
【００８２】
である。
【００８３】
最後に、先ず、第ｊ番のボクセルに対して、次式であることに留意することによって、総散乱/変位に対する共分散を表現することができる。
【数２２】
【００８４】
ここで、
【数２３】
【００８５】
であり、Ｌｉｊは、第ｊ番のボクセルを通過する第ｉ番の宇宙線の経路長であり、宇宙線の当らないボクセルに対しては、ゼロと規定される。式(19)から(22)を組み合わせると、次式を得る。
【数２４】
【００８６】
ここで、Ｎは全ボクセル数であり、各要素に対する単純ではあるが冗長な計算に基づいて、次の重み行列を定義している。
【数２５】
【００８７】
ボクセルを通る宇宙線の経路長を算定するために、未知のミューオン経路について仮定を行う。図５を参照すると、近似は、入射及び出射飛跡の再近接点(ＰｏＣＡ)の計算で始まる（ｘｃａ、ｙｃａ）。次に、入り口点とＰｏＣＡ、ＰｏＣＡと出口点を接続し、ボクセル経路長を算定する。
【００８８】
最後に、次のデータベクターを定義し、
【数２６】
【００８９】
ＤをＭ個のミューオンのすべての測定値を示すものとすると、密度プロファイルλが与えられたときのデータの尤度は次のように書ける。
【数２７】
【００９０】
上記の因子は、
【数２８】
【００９１】
である。
【００９２】
上記を考慮すると、図１５Ｄに概説されているように、一実施形態による非均一物質による多数のミューオンの散乱及び変位の確率分布を得ることができる（処理ステップ１１２）。先ず、ボクセル（又は他の基底関数）の３次元グリッドを定める（処理ステップ１７０）。次に処理ステップ１７１に於いて、各ミューオンに対して、散乱/変位の共分散行列が式(23),(24)によって計算される。最後に、所与の宇宙線経路とボクセル散乱密度に対して、すべてのミューオンに対する全体確率分布が式(25)から(27)によって計算される。これは処理ステップ１７２として示されている。
【００９３】
多重散乱統計モデルを説明したので、これから、実験効果に対するモデルの拡張について説明する（処理ステップ１１３）。現実のミューオン検出器は、有限な位置分解能を持っている。入射及び出射ミューオン飛跡は、多数の位置測定値に対する飛跡フィッティングにより得られた角度と位置によって特徴付けられる。従って、測定誤差が、ミューオントモグラフィ用データセットを構成する散乱角と変位の測定値に伝播する。特定の検出器の精度は、ＲＭＳ誤差ｅｐによって特徴付けられる。検出器の特定の配置に対して、誤差行列
【数２９】
【００９４】
は、誤差がどのように伝播するかに基づいて規定してよい。このような誤差は、反復的な再構成法では、扱いが比較的容易である。本特許の場合は、式(23)の共分散行列に追加することによって検出器誤差を考慮することができる。
【数３０】
【００９５】
このようにして、検出器誤差に起因するノイズは削減される。そうしないと、検出器誤差に起因するノイズが、再構成結果に現われる。検出器誤差に対するより正確なモデルは、運動量依存性を考慮しなければならない。追跡誤差源の1つは、検出器自身内の散乱であり、散乱は素粒子の運動量が増加すると、減少するからである。個々のミューオン運動量の推定値、
【００９６】
が利用できるなら、誤差行列
【００９７】
が各宇宙線に対して算定することができる。式(２)から明らかように、多重クーロン散乱分布の幅は、素粒子の運動量に依存する。異なるミューオンの運動量は、式(５)に因子ｐｒ2を導入することによって考慮されている。実際には、ミューオンの運動量は正確には知られていないが、個々のミューオンの運動量の推定値は、宇宙線ミューオンの既知のスペクトルなどの、既知の散乱体に於ける散乱の測定値より算定することができる。ここでは、各ミューオンに対して
【００９８】
という良い推定値があると仮定する。
【００９９】
物体容積体密度の最尤推定値は、期待値最大化アルゴリズムを使用して決定できる（方法100の処理ステップ１０３）。ＥＭアルゴリズムは、“不完全な”データの尤度を“全部揃った”データとして、即ち“見えない”データを加えた観測データとして表現することに依拠している。この適用では、観測データＤ＝｛Ｄｉ：１≦ｉ≦Ｍ｝は、測定された散乱である。見えないデータＨ＝｛Ｈｉｊ：１≦ｉ≦Ｍ＆１≦ｊ≦Ｎ｝は、第ｊ番のボクセルによる第ｉ番のミューオンの散乱角及び変位である。A.Dempster, N.Laird, D.Rubin,"Maximum likelihood from incomplete data via the EM algorithm", J.Roy.Statist.Soc.B, vol.39, pp.1-78,1977に、以下の補助関数によってアルゴリズムが説明されている。本文献の開示は参照によりここに組み込んだものとする。
【数３１】
【０１００】
この関数は、与えられたＨの条件付分布とパラメータベクトルλ(ｎ)に対して、パラメータベクトルλが与えられたときの観測データと未観測データの対数尤度の期待値である。アルゴリズムの各繰り返しは次の２ステップから構成される。
【０１０１】
Ｅステップ：見えないデータの条件付き分布Ｐ(Ｈ|Ｄ,λ(ｎ)）を推定又は特性化する。
【０１０２】
Ｍステップ：補助関数Ｑを最大化する。Ｑは、Ｅステップで特性化した分布に対する期待値である。
【０１０３】
本特許の場合、見えないデータは、観測データを一意に決めるので、より簡単な補助関数
【数３２】
【０１０４】
を使うことによって、ＱＤＬＲを使うことによって得られるのと同じ推定値列
【０１０５】
を得る。パラメータ推定値、λ(ｎ)よりアルゴリズムの1回の繰り返しにより新しい推定値、λ(ｎ+1)が以下により得られる。
【数３３】
【０１０６】
先ず、単一ボクセルによる単一ミューオンの散乱の確率分布は、単純に単一レイヤに対する統計的モデルから得られることに留意する。
【数３４】
【０１０７】
ここで、
【数３５】
【０１０８】
は、式(21)で定義されている。各ボクセルに於ける散乱の無条件確率分布は、他のボクセルの散乱とは独立であるので、見えないデータの全体集合の確率は、各要素の確率の積である。従って、対数尤度は、
【数３６】
【０１０９】
と書くことができる。ここでＣはλを含まない項を表している。条件付期待値を計算することによって、関数Ｑは、
【数３７】
【０１１０】
となり、その各加数は
【数３８】
【０１１１】
である。ここでＨｊは、Ｌｉｊ≠０である宇宙線の数（即ち、第ｊ番のボクセルに衝突する宇宙線の数）である。そしてＳｉｊ(ｎ)は次で定義される。
【数３９】
【０１１２】
式(36)のλｊについての導関数をゼロと置くと、補助関数を最大にする以下の繰り返し公式を得る（Ｍステップ）。
【数４０】
【０１１３】
Ｓｉｊの二次形式は、λ(ｎ+1)が正であることを保証する。条件付き期待値Ｓｉｊの計算が残っている。Ｘが確率変数Ｈｉｊ|Ｄｊを表すとすると、二次形式ＸＴＡ-1Ｘの期待値は、
【数４１】
【０１１４】
ここで、μｘとΣｘはＸの平均と共分散である。ＤｉはＨｉｊに一次従属しているので、これらは、同時確率ガウス分布である。またＤｉが与えられたときＨｉｊの条件付分布もガウス分布である。これは多変量分布理論からの結果である。この理論と、ＨｉｊとＤｉはそれぞれ平均値ゼロを持つという事実を使うと、以下を得る。
【数４２】
【数４３】
【０１１５】
ここで、観測されるデータの共分散、ΣＤｉ、は式(29)で与えられ、見えないデータ要素の共分散、ΣＨｉｊ、は式(21)を介して表される。観測されるー見えないデータの共分散、ΣＤｉＨｉｊを陽に書くよりも、簡単な（冗長だが）行列計算を行って、次式を示すことができる。
【数４４】
【０１１６】
式(39)から(42)からの結果を式(37)へ代入すると、次を得る。
【数４５】
【０１１７】
ここで、最後のステップで、Ｔｒ(ＡＢ)＝Ｔｒ(ＢＡ)を用いた。
【０１１８】
最後に、ｘとｙ座標の散乱データを組み入れるために、更新の式(38)に於いて、わかり易く次の平均を使用する。
【数４６】
【０１１９】
式(38)は、１つのボクセルに衝突する宇宙線に亙っての平均を表していることを留意する。この式の用法を今後、平均値法と呼ぶ。以下に、この更新式の別の形式が、局外のミューオンデータに起因するノイズの削減に有用であることを示す。アルゴリズムの中央値法は、改正した更新の次式によって定義される。
【数４７】
【０１２０】
上記を考慮して、期待値最大化アルゴリズムを使用する、物体容積体の推定密度プロファイルの尤度を最大化する処理（方法100の処理ステップ１０３）が、図１５Ｅに概説されている。処理ステップ１８０で示すように、ｉ＝1からＭ(△θｘ、△θｙ、△ｘ、△ｙ、ｐｒ2)ｉまでの各ミューオンに対して、散乱と運動量の測定値が集められる。各ミューオンとｊ＝1からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊが推定される（処理ステップ１８１）。ミューオンとボクセルの各対に対して、重み行列：Ｗｉｊが式(24)を使って計算される（処理ステップ１８２）。各ボクセルの散乱密度が、推測値λｊ,ｏｌｄで初期化される（処理ステップ１８３）。処理ステップ１８４に於いて、停止基準が以下のように示されている。各ミューオンに対して、式(29)を使用して
【０１２１】
が計算され、その逆行列が取られる。ミューオンとボクセルの各対に対して、条件付き平均値の項：Ｓｉｊが式(43)と(44)を使用して計算される。式(38)又は(45)を使用して、λｊ,ｎｅｗが計算され、すべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗと置かれる。
【０１２２】
さらに、方法100を説明するために、これから数値例について述べる。図１に示す構成と同じ構成がシミュレーションされた。最初の妥当性確認試験として、多重統計散乱モデルに厳密に適合するように設計された簡素なシミュレーションが使用された。単一の検出器プレーン（図に示してある３枚ではない）のサイズは２×２ｍ2にされた。上部と下部の検出器間の垂直距離は１．１ｍであった。これらの検出器はミューオンの位置と角度を完璧に記録した。運動量が５００−１０，０００ＭｅＶ/ｃに一様に分布するミューオンを持つ、単純化されたミューオンスペクトルが使用された。素粒子は、垂直から一様に広がる投影角度で上部の検出器プレーンの容積体に入射した。
【０１２３】
ミューオンの多重散乱と変位が処理ステップ１１０から１１３に従ってシミュレーションされた。図７と８に視覚化されているように、物体は、容積体の中央１．１×１．１×１．１ｍ3の部分内に置かれた。３個の１０×１０×１０ｍ3の立方体物質、散乱密度がそれぞれ７１．５、１４．２、２．８ｍｒａｄ2/ｃｍであるタングステン(Ｗ)、鉄(Ｆｅ)、アルミニウム(Ａｌ)がシミュレーションされた。シミュレーションは４００，０００個のミューオンを仮定し、これらは上部検出器スタックに入射し、約１０分の露出量に対応する。これらのミューオンの約１６０，０００個が下部検出器プレーンに衝突せず、再構成用には２４０、０００個が残った。再構成には、５×５×５ｍ3のボクセルサイズが使用され、上記で述べた平均値法が実装され、これは、各ミューオンに対する運動量の完全な知識を仮定している。シミュレーションは、空気が充填された容積体から開始し、アルゴリズムを１００回の繰り返しまで走らせた(ブロックの特徴の収束には十分である)。結果は図９に表示されている。３個の物体の各々の８個のボクセルに対して再構成された値の平均値は、Ｗ、Ｆｅ、Ａｌのブロックに対してそれぞれ７４．０、１４．７、２．７である。各立方体を構成する８個のボクセルの変動係数(ｒｍｓ/ｍｅａｎ)は１２．６％、１３．２％、１２．１％である。この結果は、シミュレーションとインバージョンモデルの一致を考えると、インバージョンアルゴリズムと実装を確証している。
【０１２４】
再構成結果は、物体シーンに一致しているように見える。３個の物体の各々の８個のボクセルに対して再構成された値の平均値は、Ｗ、Ｆｅ、Ａｌのブロックに対してそれぞれ７４．０、１４．７、２．７である。各立方体を構成する８個のボクセルの変動係数(ｒｍｓ/ｍｅａｎ)は１２．６％、１３．２％、１２．１％である。この結果は、シミュレーションとインバージョンモデルの一致を考えると、インバージョンアルゴリズムと実装を確証している。
【０１２５】
次に同じシーンをＧＥＡＮＴ４のモンテカルロパッケージを使用してシミュレーションを行った。ＧＥＡＮＴ４の詳細は、J.Allison, “Geant4 developments and applications”, IEEE Trans. Nucl. Sci., vol.53, no.1, pp.270-278, Feb. 2006の刊行物にある。この開示は参照によりここに組み込んだものとする。ＧＥＡＮＴ４は、多重散乱に対する、より完全で、正確でかつ妥当性が実証されたモデルを実装している。このモデルは、散乱分布の中央のガウス分布の部分の幅のより正確な計算、大きいテール部分の実装、物質を通過する時のミューオンのエネルギーロスのシミュレーションを含んでいる。ミューオンイベントジェネレータも使用した。これは、海面レベルの宇宙線ミューオンの角度及び運動量の分布を再現する。このシミュレーションでは、各ミューオンの完全な知識とともに、検出器は完璧であることが仮定された。宇宙線電子や飛跡２次粒子は含められなかった。結果は、図１０に表示されている。Ｗ,Ｆｅ,Ａｌブロックに対応するボクセル値の平均は、それぞれ６７４．４、６３．４、５．４である。
【０１２６】
ボクセル値はあまりにも大きく、中位及び低Ｚの領域のいくつかの誤判別は明らかである。すべてのボクセル値を約４で除して再構成値を正規化することによって、中位のＺのボクセルに対して正しい平均ボクセル値を得られるが、高及び低Ｚのボクセルに対しては正しい値を与えない、又はすべての誤判別の解消はない。この結果の原因は、ガウス分布モデルで旨く記述されていないように散乱する少数のミューオンである。散乱の投影角度分布の中央９８％は、ガウス分布としてよく近似されていると断言される。約２％のミューオンは、ここで説明した統計的モデルに比較して大きい角度で、即ち、ガウス分布の場合に見られるより非常に大きい角度で散乱する。散乱角度の２乗がフィッティングを決定するので、影響は劇的である。これらのテール部分にあるミューオン散乱は、非常に大きい散乱密度推定値を形成する。
【０１２７】
さらに、図１の機器内のミューオンの崩壊などのその他の過程や重大な検出器エラーは、非常に大きな散乱事象として誤って記録される可能性がある（とはいえ、本願に於けるシミュレーションには、これらの発生源は存在しなかった）。これは、容積体のどこでも発生し得るし、不合理に大きい散乱密度を持つ単独のボクセルを作り出す傾向がある。このような事象は、ＳＮＭについて誤った陽性表示を与えるので、除外する必要がある。
【０１２８】
ＥＭアルゴリズムが非ガウス分布データに対して耐性を持つようにするために、平均値更新規則の式(38)は、式(45)、即ち中央値法の使用で置き換えることができる。
【０１２９】
中央値法を使った結果が、図１１に示されている。Ｗ,Ｆｅ,Ａｌ領域に対するボクセル平均値は、それぞれ７９．２、１４．２、２．１であり、変動係数は２１．５％、２６．３％、２３．２％である。明らかに、中央値法を使うとインバージョンアルゴリズムの頑強性が向上する。
【０１３０】
図１２に、再構成される密度プロファイルのより実際的な例が示されており、乗用バンの詳細ＧＥＡＮＴ４シミュレーションを図解している。バン・ボディーが切り取られた主要構成要素のイラストが示されている。イラストの中央の赤いブロックは、１０×１０×１０ｃｍ3のタングステンの固体片、即ち高Ｚの危険物体の代用物を表す。この場合、そのシーンの４個の長い側面上に配置した模擬の検出器プレーンを使用し、より水平方向に向かうミューオンを利用した。宇宙線ミューオンに対する１分の暴露量がシミュレーションされ、平均値法及び中央値法に従い、５×５×５ｃｍ2のサイズのボクセルを使用して、データから再構成が実施された。図１３と１４は、それぞれ、平均値ＥＭ法と中央値ＥＭ法を使用して行った再構成結果を可視化したものを示している。このシーンの平均値法の再構成結果に於ける、非ガウス分布データの影響は、容易に判り、画像上に散らばったより黒いスポットとして表れている。中央値法の再構成結果に於いては、これらのアーチファクトは全く無く、バンのより密度の高い構成要素（エンジン、バッテリ、駆動系）が（低Ｚ）又は(中位のＺ)として表れており、他方、危険物体は、より黒く(高Ｚ)目立っている。中央値法の使用によって、非ガウス分布の散乱分布及び他の異常事象より生じる誤検出に対して堅牢な傾向の結果が得られる。
【０１３１】
ここに示した実施形態と例は、本発明とその実際の応用を最もよく説明し、それによって同業者が本発明を成し利用できるようにするために提示された。しかし同業者は、上記説明及び実施例は、説明と例のためにのみ提示されたことを認識しなければならない。
【０１３２】
本発明の他の変形及び変更（再構成の処理の調整など）は、同業者にとって明らかであり、このような変形と変更がカバーされることは、添付の請求範囲の意図である。
【０１３３】
記載した説明は、包括的であること、又は本発明の範囲を制限することを意図していない。上記教示を考慮すると、以下に続く請求範囲から逸脱することなく、多くの変更及び変形が可能である。本発明の使用には、様々な特徴を有する構成要素が含まれると考えられる。
【特許請求の範囲】
【請求項１】
物体容積体から得られる荷電粒子トモグラフィ・データより物体容積体を検出する方法であって：
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、
(ｂ)期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される、統計的多重散乱モデルに基づく、荷電粒子散乱の確率分布を用意する工程と、
(ｃ)前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程と、
(ｄ)前記略最尤推定値を基に再構成された物体容積体散乱密度を出力する工程と、
を含む方法。
【請求項２】
前記再構成された物体容積体散乱密度に基づいて、物体容積体内の標的物体の(1)存在と(2)種類のうちの少なくとも１つに関して判定を行う工程を含む、請求項１に記載の方法。
【請求項３】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される荷電粒子散乱の確率分布を用意する工程は、
(ｇ)均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元の確率分布を得る工程と、
(ｈ)前記２次元の確率分布に基づいて、前記荷電粒子に対する３次元の確率分布を得る工程と、
(ｉ)基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程と、
(ｊ)その定義、及び前記荷電粒子の散乱及び運動量の測定値に基づいて、多重散乱に対する確率分布を得る工程と、
を含む、請求項１に記載の方法。
【請求項４】
前記均一物体の所定の散乱密度に基づいて荷電粒子に対する２次元の確率分布を得る工程は、
(ｋ)物質の散乱密度を、その物質の単位深度による荷電粒子の散乱の２乗の期待平均値として決定する工程と、
(ｌ)荷電粒子の散乱分布をガウス分布モデルに基づいて近似する工程と、
(ｍ)荷電粒子宇宙線の散乱及び変位を相関２次元ガウス分布に基づいて近似する工程と、を含む、請求項３に記載の方法。
【請求項５】
前記２次元の確率分布に基づいて、前記荷電粒子に対する３次元の確率分布を得る工程は、
座標を追加し、３次元経路長を規定する工程と、
３次元変位を計算する工程と、
３次元共分散行列を規定する工程と、
を含む、請求項３に記載の方法。
【請求項６】
前記基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程は、
離散的散乱密度推定値を表す基底関数の３次元グリッドを設ける工程と、
各ミューオンに対する前記散乱/変位の共分散行列を宇宙線経路と散乱密度の推定値の関数として決定する工程と、
複数の荷電粒子に対する確率分布を個々の荷電粒子に対する確率分布の積として決定する工程と、
を含む、請求項５に記載の方法。
【請求項７】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
各荷電粒子に対する散乱及び運動量の測定値を収集する工程と、
各荷電粒子と前記統計的散乱モデルの各基底関数との相互作用の幾何学を推定する工程と、
荷電粒子と基底関数の各対に対して、重み行列：Ｗｉｊを決定する工程と、
各基底関数の散乱密度を推測値で初期化する工程と、
物体容積体の内容に関して、略最尤解を反復的に求める工程と、
を含み、
前記反復処理は、所定の反復回数で、又は解が所定の許容値未満になったときに、停止される、請求項１に記載の方法。
【請求項８】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
ｉ＝１からＭ(△θｘ、△θｙ、△ｘ、△ｙ,ｐｒ２)ｉまでの各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、
各ミューオンとｊ＝１からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊを推定する工程と、
荷電粒子とボクセルの各対に対して
【数４８】
として、重み行列：Ｗｉｊを計算する工程と、
各ボクセルに於ける散乱密度λｊ,ｏｌｄの推測値を初期化する工程と、
停止基準処理を使用してすべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗ と置く工程と、
を含む、請求項１に記載の方法。
【請求項９】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムは、平均値更新規則又は中央値更新規則を含む、請求項１の方法。
【請求項１０】
前記荷電粒子トモグラフィ・データは、宇宙線ミューオントモグラフィ・データを含む、請求項１に記載の方法。
【請求項１１】
物体容積体から得られる荷電粒子トモグラフィ・データより物体容積体を検出するコンピュータによって実行される方法であって：
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、
(ｂ)期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される、統計的多重散乱モデルに基づく、荷電粒子散乱密度の確率分布を用意する工程と、
(ｃ)前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程と、
(ｄ)再構成された物体容積体散乱密度を出力する工程と、
を含む方法。
【請求項１２】
前記再構成された物体容積体散乱密度に基づいて、判定を行う工程をさらに含む、請求項１１に記載の方法。
【請求項１３】
期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される荷電粒子散乱密度の確率分布を用意する工程は、
(ｇ)均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元の確率分布を推定する工程と、
(ｈ)前記２次元の確率分布に基づいて、前記荷電粒子に対する３次元の確率分布を得る工程と、
(ｉ)基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程と、
(ｊ)その定義、及び前記荷電粒子の散乱及び運動量の測定値に基づいて、多重散乱に対する前記確率分布を決定する工程と、
を含む、請求項１１に記載の方法。
【請求項１４】
前記荷電粒子トモグラフィ・データは、宇宙線ミューオントモグラフィ・データを含む、請求項１３に記載の方法。
【請求項１５】
均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元の確率分布を推定する工程は、
物質の散乱密度を、公称運動量ｐ0＝３ＧｅＶ/ｃの宇宙線荷電粒子のその物質の単位深度による散乱の２乗の期待平均値として決定する工程と、
宇宙線荷電粒子散乱分布を、ガウス分布モデルを使用して近似する工程と、
宇宙線散乱及び変位を相関２次元ガウス分布によって近似する工程と、
を含み、
前記２次元確率分布に基づき、前記荷電粒子に対する３次元確率分布を得る工程は、
座標を追加して、３次元経路長を規定する工程と、
３次元変位を計算する工程と、
３次元共分散行列を規定する工程と、
を含み、
基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程は、
離散的散乱密度推定値を表す基底関数の３次元グリッドを設定する工程と、
個々の宇宙線荷電粒子に対する前記散乱/変位の共分散行列を宇宙線経路及び散乱密度推定値の関数として決定する工程と、
複数の宇宙線荷電粒子に対する確率分布を個々の荷電粒子の確率分布の積として決定する工程と、
を含む、請求項１４に記載の方法。
【請求項１６】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
各宇宙線荷電粒子に対する散乱及び運動量の測定値を収集する工程と、
各荷電粒子と前記統計的多重散乱モデルの各基底関数との相互作用の幾何学を推定する工程と、
荷電粒子と基底関数の各対に対して、重み行列：Ｗｉｊを決定する工程と、
各基底関数の散乱密度を推測値で初期化する工程と、
物体容積体の内容に関して、略最尤解を反復的に求める工程と、
を含み、
前記反復処理は、所定の反復回数で、又は解が所定の許容値未満になったときに、停止される、請求項１４に記載の方法。
【請求項１７】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
ｉ＝１からＭ(△θｘ、△θｙ、△ｘ、△ｙ,ｐｒ２)ｉまでの各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、
各荷電粒子とｊ＝１からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊを推定する工程と、
荷電粒子とボクセルの各対に対して
【数４９】
として、重み行列：Ｗｉｊを計算する工程と、
各ボクセルに於ける散乱密度λｊ,ｏｌｄの推測値を初期化する工程と、
停止基準処理を使用してすべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗと置く工程と、
を含む、請求項１０に記載の方法。
【請求項１８】
前記停止基準処理は、
各荷電粒子に対して、
【数５０】
を使用して、
を計算し、その逆行列を取る工程と、
荷電粒子とボクセルの各対に対して、
【数５１】
を使用して、条件付き期待値の項：Ｓｉｊを計算する工程と、
【数５２】
を使用して、ｘ座標及びｙ座標散乱データを更新規則に組み込む工程と、
を含み、
前記Ｓｉｊ(ｎ)の式の最後のステップに於いて、Ｔｒ(ＡＢ）＝Ｔｒ(ＢＡ)が適用される、請求項１７に記載の方法。
【請求項１９】
前記荷電粒子は、ミューオンを含む、請求項１８に記載の方法。
【請求項２０】
前記ＭＬ/ＥＭアルゴリズムは、
【数５３】
として定義される平均値更新規則又は、
【数５４】
として定義される中央値更新規則を含む、請求項１８に記載の方法。
【請求項２１】
コンピュータアクセス可能なデータ記憶媒体を含むコンピュータ・プログラム製品であって、そのデータ記憶媒体は、コンピュータによって実行される時荷電粒子トモグラフィ・データから物体容積体密度プロファイルを統計的に再構成する方法をコンピュータに実施させる命令を格納し、
前記方法は、
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、
(ｂ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）に使用される、統計的多重散乱モデルに基づく、荷電粒子散乱密度の確率分布を用意する工程と、
(ｃ)前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体密度の略最尤推定値を決定する工程と、
(ｄ)再構成された物体容積体散乱密度を出力する工程と、
を含む、コンピュータ・プログラム製品。
【請求項２２】
物体容積体を通過する荷電粒子によって物体容積体を検出する検出システムであって、
物体容積体の第１の側に配置され、物体容積体に向かう入射荷電粒子の位置と角度を測定する第１の位置感知検出器セットと、
物体容積体の前記第１の側の反対側の第２の側に配置され、物体容積体から出て行く出射荷電粒子の位置と角度を測定する第２の位置感知検出器セットと、
前記第１の位置感知検出器セットからの測定信号のデータと前記第２の位置感知検出器セットからの測定信号のデータを受け取る信号処理ユニットと、
を含み、
前記信号処理ユニットは、前記受信データを処理し、物体容積体内の容積体散乱密度分布を統計的に再構成する、前記検出システム。
【請求項２３】
前記信号処理ユニットは、
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する荷電粒子トモグラフィ・データを取得し、
(ｂ)統計的多重散乱モデルに基づいて、荷電粒子散乱密度の確率分布を用意し、
(ｃ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の略最尤推定値を決定し、
(ｄ)前記略最尤推定値に基づいて再構成された物体容積体散乱密度を出力するように構成される、請求項２２に記載のシステム。
【請求項２４】
素粒子検出器の前記第１及び第２のセットは各々、第１の方向の少なくとも３つの荷電粒子測定及び前記第１の方向とは異なる第２の方向の少なくとも３つの荷電粒子測定が可能となるように配置されたドリフトチューブを含む、請求項２２のシステム。
【請求項２５】
前記荷電粒子は、物体容積体に入射する自然の宇宙線ミューオンであって、そして前記信号処理ユニットは、物体容積体内の容積体散乱密度分布の統計的再構成結果に基づいて、物体容積体内に標的物体が存在するかどうかを示すように構成される、請求項２２に記載のシステム。
【請求項１】
物体容積体から得られる荷電粒子トモグラフィ・データより物体容積体を検出する方法であって：
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、
(ｂ)期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される、統計的多重散乱モデルに基づく、荷電粒子散乱の確率分布を用意する工程と、
(ｃ)前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程と、
(ｄ)前記略最尤推定値を基に再構成された物体容積体散乱密度を出力する工程と、
を含む方法。
【請求項２】
前記再構成された物体容積体散乱密度に基づいて、物体容積体内の標的物体の(1)存在と(2)種類のうちの少なくとも１つに関して判定を行う工程を含む、請求項１に記載の方法。
【請求項３】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される荷電粒子散乱の確率分布を用意する工程は、
(ｇ)均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元の確率分布を得る工程と、
(ｈ)前記２次元の確率分布に基づいて、前記荷電粒子に対する３次元の確率分布を得る工程と、
(ｉ)基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程と、
(ｊ)その定義、及び前記荷電粒子の散乱及び運動量の測定値に基づいて、多重散乱に対する確率分布を得る工程と、
を含む、請求項１に記載の方法。
【請求項４】
前記均一物体の所定の散乱密度に基づいて荷電粒子に対する２次元の確率分布を得る工程は、
(ｋ)物質の散乱密度を、その物質の単位深度による荷電粒子の散乱の２乗の期待平均値として決定する工程と、
(ｌ)荷電粒子の散乱分布をガウス分布モデルに基づいて近似する工程と、
(ｍ)荷電粒子宇宙線の散乱及び変位を相関２次元ガウス分布に基づいて近似する工程と、を含む、請求項３に記載の方法。
【請求項５】
前記２次元の確率分布に基づいて、前記荷電粒子に対する３次元の確率分布を得る工程は、
座標を追加し、３次元経路長を規定する工程と、
３次元変位を計算する工程と、
３次元共分散行列を規定する工程と、
を含む、請求項３に記載の方法。
【請求項６】
前記基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程は、
離散的散乱密度推定値を表す基底関数の３次元グリッドを設ける工程と、
各ミューオンに対する前記散乱/変位の共分散行列を宇宙線経路と散乱密度の推定値の関数として決定する工程と、
複数の荷電粒子に対する確率分布を個々の荷電粒子に対する確率分布の積として決定する工程と、
を含む、請求項５に記載の方法。
【請求項７】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
各荷電粒子に対する散乱及び運動量の測定値を収集する工程と、
各荷電粒子と前記統計的散乱モデルの各基底関数との相互作用の幾何学を推定する工程と、
荷電粒子と基底関数の各対に対して、重み行列：Ｗｉｊを決定する工程と、
各基底関数の散乱密度を推測値で初期化する工程と、
物体容積体の内容に関して、略最尤解を反復的に求める工程と、
を含み、
前記反復処理は、所定の反復回数で、又は解が所定の許容値未満になったときに、停止される、請求項１に記載の方法。
【請求項８】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
ｉ＝１からＭ(△θｘ、△θｙ、△ｘ、△ｙ,ｐｒ２)ｉまでの各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、
各ミューオンとｊ＝１からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊを推定する工程と、
荷電粒子とボクセルの各対に対して
【数４８】
として、重み行列：Ｗｉｊを計算する工程と、
各ボクセルに於ける散乱密度λｊ,ｏｌｄの推測値を初期化する工程と、
停止基準処理を使用してすべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗ と置く工程と、
を含む、請求項１に記載の方法。
【請求項９】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムは、平均値更新規則又は中央値更新規則を含む、請求項１の方法。
【請求項１０】
前記荷電粒子トモグラフィ・データは、宇宙線ミューオントモグラフィ・データを含む、請求項１に記載の方法。
【請求項１１】
物体容積体から得られる荷電粒子トモグラフィ・データより物体容積体を検出するコンピュータによって実行される方法であって：
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、
(ｂ)期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される、統計的多重散乱モデルに基づく、荷電粒子散乱密度の確率分布を用意する工程と、
(ｃ)前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程と、
(ｄ)再構成された物体容積体散乱密度を出力する工程と、
を含む方法。
【請求項１２】
前記再構成された物体容積体散乱密度に基づいて、判定を行う工程をさらに含む、請求項１１に記載の方法。
【請求項１３】
期待値最大化(ＭＬ/ＥＭ)アルゴリズムに於いて使用される荷電粒子散乱密度の確率分布を用意する工程は、
(ｇ)均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元の確率分布を推定する工程と、
(ｈ)前記２次元の確率分布に基づいて、前記荷電粒子に対する３次元の確率分布を得る工程と、
(ｉ)基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程と、
(ｊ)その定義、及び前記荷電粒子の散乱及び運動量の測定値に基づいて、多重散乱に対する前記確率分布を決定する工程と、
を含む、請求項１１に記載の方法。
【請求項１４】
前記荷電粒子トモグラフィ・データは、宇宙線ミューオントモグラフィ・データを含む、請求項１３に記載の方法。
【請求項１５】
均一物体の所定の散乱密度に基づいて、荷電粒子に対する２次元の確率分布を推定する工程は、
物質の散乱密度を、公称運動量ｐ0＝３ＧｅＶ/ｃの宇宙線荷電粒子のその物質の単位深度による散乱の２乗の期待平均値として決定する工程と、
宇宙線荷電粒子散乱分布を、ガウス分布モデルを使用して近似する工程と、
宇宙線散乱及び変位を相関２次元ガウス分布によって近似する工程と、
を含み、
前記２次元確率分布に基づき、前記荷電粒子に対する３次元確率分布を得る工程は、
座標を追加して、３次元経路長を規定する工程と、
３次元変位を計算する工程と、
３次元共分散行列を規定する工程と、
を含み、
基底関数によって特徴付けられる非均一物体容積体による多数の荷電粒子の散乱に対する確率分布を得る工程は、
離散的散乱密度推定値を表す基底関数の３次元グリッドを設定する工程と、
個々の宇宙線荷電粒子に対する前記散乱/変位の共分散行列を宇宙線経路及び散乱密度推定値の関数として決定する工程と、
複数の宇宙線荷電粒子に対する確率分布を個々の荷電粒子の確率分布の積として決定する工程と、
を含む、請求項１４に記載の方法。
【請求項１６】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
各宇宙線荷電粒子に対する散乱及び運動量の測定値を収集する工程と、
各荷電粒子と前記統計的多重散乱モデルの各基底関数との相互作用の幾何学を推定する工程と、
荷電粒子と基底関数の各対に対して、重み行列：Ｗｉｊを決定する工程と、
各基底関数の散乱密度を推測値で初期化する工程と、
物体容積体の内容に関して、略最尤解を反復的に求める工程と、
を含み、
前記反復処理は、所定の反復回数で、又は解が所定の許容値未満になったときに、停止される、請求項１４に記載の方法。
【請求項１７】
前記期待値最大化(ＭＬ/ＥＭ)アルゴリズムを使用して、物体容積体散乱密度の略最尤推定値を決定する工程は、
ｉ＝１からＭ(△θｘ、△θｙ、△ｘ、△ｙ,ｐｒ２)ｉまでの各荷電粒子に対して、散乱及び運動量の測定値を収集する工程と、
各荷電粒子とｊ＝１からＮまでの各ボクセルとの相互作用の幾何学：(Ｌ、Ｔ)ｉｊを推定する工程と、
荷電粒子とボクセルの各対に対して
【数４９】
として、重み行列：Ｗｉｊを計算する工程と、
各ボクセルに於ける散乱密度λｊ,ｏｌｄの推測値を初期化する工程と、
停止基準処理を使用してすべてのボクセルに対してλｊ,ｏｌｄ＝λｊ,ｎｅｗと置く工程と、
を含む、請求項１０に記載の方法。
【請求項１８】
前記停止基準処理は、
各荷電粒子に対して、
【数５０】
を使用して、
を計算し、その逆行列を取る工程と、
荷電粒子とボクセルの各対に対して、
【数５１】
を使用して、条件付き期待値の項：Ｓｉｊを計算する工程と、
【数５２】
を使用して、ｘ座標及びｙ座標散乱データを更新規則に組み込む工程と、
を含み、
前記Ｓｉｊ(ｎ)の式の最後のステップに於いて、Ｔｒ(ＡＢ）＝Ｔｒ(ＢＡ)が適用される、請求項１７に記載の方法。
【請求項１９】
前記荷電粒子は、ミューオンを含む、請求項１８に記載の方法。
【請求項２０】
前記ＭＬ/ＥＭアルゴリズムは、
【数５３】
として定義される平均値更新規則又は、
【数５４】
として定義される中央値更新規則を含む、請求項１８に記載の方法。
【請求項２１】
コンピュータアクセス可能なデータ記憶媒体を含むコンピュータ・プログラム製品であって、そのデータ記憶媒体は、コンピュータによって実行される時荷電粒子トモグラフィ・データから物体容積体密度プロファイルを統計的に再構成する方法をコンピュータに実施させる命令を格納し、
前記方法は、
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する所定の荷電粒子トモグラフィ・データを取得する工程と、
(ｂ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）に使用される、統計的多重散乱モデルに基づく、荷電粒子散乱密度の確率分布を用意する工程と、
(ｃ)前記期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体密度の略最尤推定値を決定する工程と、
(ｄ)再構成された物体容積体散乱密度を出力する工程と、
を含む、コンピュータ・プログラム製品。
【請求項２２】
物体容積体を通過する荷電粒子によって物体容積体を検出する検出システムであって、
物体容積体の第１の側に配置され、物体容積体に向かう入射荷電粒子の位置と角度を測定する第１の位置感知検出器セットと、
物体容積体の前記第１の側の反対側の第２の側に配置され、物体容積体から出て行く出射荷電粒子の位置と角度を測定する第２の位置感知検出器セットと、
前記第１の位置感知検出器セットからの測定信号のデータと前記第２の位置感知検出器セットからの測定信号のデータを受け取る信号処理ユニットと、
を含み、
前記信号処理ユニットは、前記受信データを処理し、物体容積体内の容積体散乱密度分布を統計的に再構成する、前記検出システム。
【請求項２３】
前記信号処理ユニットは、
(ａ)物体容積体を通過する荷電粒子の散乱角及び推定運動量に対応する荷電粒子トモグラフィ・データを取得し、
(ｂ)統計的多重散乱モデルに基づいて、荷電粒子散乱密度の確率分布を用意し、
(ｃ)期待値最大化アルゴリズム（ＭＬ/ＥＭ）を用いて、物体容積体散乱密度の略最尤推定値を決定し、
(ｄ)前記略最尤推定値に基づいて再構成された物体容積体散乱密度を出力するように構成される、請求項２２に記載のシステム。
【請求項２４】
素粒子検出器の前記第１及び第２のセットは各々、第１の方向の少なくとも３つの荷電粒子測定及び前記第１の方向とは異なる第２の方向の少なくとも３つの荷電粒子測定が可能となるように配置されたドリフトチューブを含む、請求項２２のシステム。
【請求項２５】
前記荷電粒子は、物体容積体に入射する自然の宇宙線ミューオンであって、そして前記信号処理ユニットは、物体容積体内の容積体散乱密度分布の統計的再構成結果に基づいて、物体容積体内に標的物体が存在するかどうかを示すように構成される、請求項２２に記載のシステム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１５Ａ】
【図１５Ｂ】
【図１５Ｃ】
【図１５Ｄ】
【図１５Ｅ】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１５Ａ】
【図１５Ｂ】
【図１５Ｃ】
【図１５Ｄ】
【図１５Ｅ】
【公表番号】特表２０１０−５０８５２２（Ｐ２０１０−５０８５２２Ａ）
【公表日】平成２２年３月１８日（２０１０．３．１８）
【国際特許分類】
【出願番号】特願２００９−５３４９１６（Ｐ２００９−５３４９１６）
【出願日】平成１９年１０月２６日（２００７．１０．２６）
【国際出願番号】ＰＣＴ／ＵＳ２００７／０８２７５３
【国際公開番号】ＷＯ２００８／１４０５６０
【国際公開日】平成２０年１１月２０日（２００８．１１．２０）
【出願人】（５０９１１７８５０）ロス アラモス ナショナル セキュリティー，エルエルシー (3)
【Ｆターム（参考）】
【公表日】平成２２年３月１８日（２０１０．３．１８）
【国際特許分類】
【出願日】平成１９年１０月２６日（２００７．１０．２６）
【国際出願番号】ＰＣＴ／ＵＳ２００７／０８２７５３
【国際公開番号】ＷＯ２００８／１４０５６０
【国際公開日】平成２０年１１月２０日（２００８．１１．２０）
【出願人】（５０９１１７８５０）ロス アラモス ナショナル セキュリティー，エルエルシー (3)
【Ｆターム（参考）】
[ Back to top ]