JP5963161B2 - 内部状態解析方法およびプログラム並びに内部状態解析システム - Google Patents

Description

本発明は、内部状態解析方法およびプログラム並びに内部状態解析システムに関する。

従来、１２本×１２本のミュオンセンサーモジュールが直交して並べられてなる２層の検出器からなるホドスコープと、各センサーモジュールが接続されたミュオンリードアウトモジュールと、を有するミュオンの計測装置が提案されている（例えば、特許文献１参照）。この計測装置では、ミュオンリードアウトモジュールを、基盤と、基盤に実装されると共にミュオンセンサーモジュールに接続されてミュオンセンサからの信号を処理して角度分布ヒストグラムを生成してメモリに蓄積するミュオンリードアウト処理回路と、基盤に実装されてヒストグラムデータをメモリから呼び出して外部システムに出力するイーサネット（登録商標）インタフェースと、を有するモジュールとすることにより、巨大物体の内部構造を遠隔地からリアルタイムで可視化して解析することができるようにしている。

特開２０１０−１０１８９２号公報

一般に、ミュオン（ミュー粒子）などの高透過性を有する高透過性粒子を検出可能な検出部（上述のホドスコープなど）を山や産業プラントなどの測定対象の周辺に配置して検出部による検出結果を用いて測定対象の内部の状態を解析する場合、検出部の各部分のバラツキなどによる影響をより抑制して、測定対象の内部の状態をより適正に解析できるようにすることが望まれている。

本発明の内部状態解析方法およびプログラム並びに内部状態解析システムは、測定対象の内部の状態をより適正に解析できるようにすることを主目的とする。

本発明の内部状態解析方法およびプログラム並びに内部状態解析システムは、上述の主目的を達成するために以下の手段を採った。

本発明の内部状態解析方法は、
所定の高透過性を有する高透過性粒子の通過を検出可能な検出部を有する検出装置が測定対象の周辺に配置されたときの前記検出部による検出結果を用いて前記測定対象の内部の状態を解析する内部状態解析方法であって、
前記検出装置が前記測定対象の周辺に配置された第１状態での前記検出部による検出結果に応じた前記高透過性粒子のフラックスである第１状態検出フラックスと、前記検出部で前記高透過性粒子の通過を検出可能な該検出部から見た方向の範囲である検出方向範囲が前記第１状態から変更された第２状態での前記検出部による検出結果に応じた前記高透過性粒子のフラックスである第２状態検出フラックスと、の比を用いて前記測定対象の密度関連物理量を演算する、
ことを特徴とする。

この本発明の内部状態解析方法では、所定の高透過性を有する高透過性粒子の通過を検出可能な検出部を有する検出装置が測定対象の周辺に配置された第１状態での検出部による検出結果に応じた高透過性粒子のフラックスである第１状態検出フラックスと、検出部で高透過性粒子の通過を検出可能な検出部から見た方向の範囲である検出方向範囲が第１状態から変更された第２状態での検出部による検出結果に応じた高透過性粒子のフラックスである第２状態検出フラックスと、の比を用いて測定対象の密度関連物理量を演算する。即ち、第１状態検出フラックスと第２状態検出フラックスとの比を用いて測定対象の密度関連物理量を演算するのである。これにより、検出部の各部分のバラツキなどによる影響を抑制することができ、測定対象の内部の状態をより適正に解析することができる。ここで、「高透過性粒子」は、ミュオンである、ものとすることもできる。また、「密度関連物理量」は、密度（方位角方向や仰角方向の平均密度を含む）や、密度分布（上限密度，下限密度）である、ものとすることもできる。

こうした本発明の内部状態解析方法において、前記第２状態は、前記第１状態から前記検出部が方位角方向に所定角度だけ回転された状態であり、前記第１状態における前記検出方向範囲である第１方向範囲の各方向の前記第１状態検出フラックスと、前記第２状態における前記検出方向範囲である第２方向範囲の各方向の前記第２状態検出フラックスと、のうち互いに前記所定角度だけ異なる方向の前記第１状態検出フラックスと前記第２状態検出フラックスとの比である複数の検出比を用いて前記測定対象の密度関連物理量を演算する、ものとすることもできる。ここで、「方位角方向」は、前記第１状態を基準として前記検出部を中心とする水平角度の方向を意味する。

この複数の検出比を用いて測定対象の密度関連物理量を演算する態様の本発明の内部状態解析方法において、前記複数の検出比と、前記第１方向範囲と前記第２方向範囲とからなる全体方向範囲の各方向の前記測定対象の通過経路長に基づく該全体方向範囲の各方向の前記高透過性粒子のフラックスである経路長起因フラックスのうち、互いに前記所定角度だけ異なる方向の前記経路長起因フラックスの比である複数の解析比と、を用いて前記測定対象の密度関連物理量を演算する、ものとすることもできる。

この複数の検出比と複数の解析比とを用いて測定対象の密度関連物理量を演算する態様の本発明の内部状態解析法において、前記全体方向範囲の各仰角について、各方位角の前記測定対象の通過経路長と前記測定対象の方位角方向の平均密度とに基づく各方位角の前記経路長起因フラックスを用いて得られる前記複数の解析比の方位角方向についての総和と、前記複数の検出比の方位角方向についての総和と、の差分が第１許容値以下となるよう前記測定対象の方位角方向の平均密度を演算する、ものとすることもできる。こうすれば、全体方向範囲の各仰角について、測定対象の方位角方向の密度をより適正に演算することができる。ここで、「第１許容値」は、第１状態での検出部による検出結果と第２状態での検出部による検出結果とが共にガウス分布に従うと仮定したときの各検出比の統計誤差の二乗の和の平方根である、ものとすることもできる。

また、複数の検出比と複数の解析比とを用いて全体方向範囲の各仰角における測定対象の方位角方向の密度を演算する態様の本発明の内部状態解析法において、前記第１方向範囲の各方位角，各仰角ついて、前記演算した測定対象の方位角方向の平均密度と前記測定対象の方位角方向の密度勾配と前記測定対象の通過経路長とに基づく前記経路長起因フラックスを用いて得られる前記解析比である勾配起因解析比が前記検出比より大きく該検出比に第２許容値を加えた値より小さくなるよう前記所定角度だけ異なる方位角との前記測定対象の密度比の上限としての上限密度比を演算すると共に前記勾配起因解析比が前記検出比より小さく該検出比から前記第２許容値を減じた値より大きくなるよう前記所定角度だけ異なる方位角との前記測定対象の密度比の下限としての下限密度比を演算し、前記全体方向範囲の各方位角，各仰角について、前記演算した上限密度比を用いて該測定対象の密度分布の上限としての上限密度を演算すると共に前記演算した下限密度比を用いて該測定対象の密度分布の下限としての下限密度を演算する、ものとすることもできる。こうすれば、全体方向範囲の各方位角，各仰角の密度分布（上限密度，下限密度）をより適正に演算することができる。ここで、「第２許容値」は、第１状態での検出部による検出結果と第２状態での検出部による検出結果とが共にガウス分布に従うと仮定したときの検出比の統計誤差である、ものとすることもできる。

さらに、複数の検出比と複数の解析比とを用いて全体方向範囲の各仰角における測定対象の方位角方向の密度を演算する態様の本発明の内部状態解析法において、前記第１方向範囲の各方位角，各仰角ついて、前記演算した測定対象の方位角方向の平均密度と前記測定対象の方位角方向の密度勾配と前記測定対象の通過経路長とに基づく前記経路長起因フラックスを用いて得られる前記解析比である勾配起因解析比と前記検出比との差分が第３許容値以下となるよう前記測定対象の方位角方向の密度勾配を演算し、前記全体方向範囲の各方位角，各仰角について、前記演算した測定対象の方位角方向の密度勾配を用いて前記測定対象の密度を演算する、ものとすることもできる。こうすれば、全体方向範囲の各方位角，各仰角の密度をより適正に演算することができる。ここで、「第３許容値」は、第１状態での検出部による検出結果と第２状態での検出部による検出結果とが共にガウス分布に従うと仮定したときの検出比の統計誤差である、ものとすることもできる。

本発明の内部状態解析方法において、前記第２状態は、前記第１状態から前記検出部が仰角方向に所定角度だけ回転された状態であり、前記第１状態における前記検出方向範囲である第１方向範囲の各方向の前記第１状態検出フラックスと、前記第２状態における前記検出方向範囲である第２方向範囲の各方向の前記第２状態検出フラックスと、のうち互いに前記所定角度だけ異なる方向の前記第１状態検出フラックスと前記第２状態検出フラックスとの比である複数の検出比を用いて前記測定対象の密度関連物理量を演算する、ものとすることもできる。この場合、前記複数の検出比と、前記第１方向範囲と前記第２方向範囲とからなる全体方向範囲の各方向の前記測定対象の通過経路長に基づく該全体方向範囲の各方向の前記高透過性粒子のフラックスである経路長起因フラックスのうち、互いに前記所定角度だけ異なる方向の前記経路長起因フラックスの比である複数の解析比と、を用いて前記測定対象の密度関連物理量を演算する、ものとすることもできる。ここで、「仰角方向」は、前記第１状態を基準として前記検出部を中心とする仰角角度の方向を意味する。

また、本発明の内部状態解析方法において、前記第２状態は、前記第１状態から前記検出部が水平方向および／または鉛直方向に所定距離だけ平行移動された状態であり、前記第１状態における前記検出方向範囲である第１方向範囲の各方向の前記第１状態検出フラックスと、前記第２状態における前記検出方向範囲である第２方向範囲の各方向の前記第２状態検出フラックスと、のうち同一方向の前記第１状態検出フラックスと前記第２状態検出フラックスとの比である複数の検出比を用いて前記測定対象の密度関連物理量を演算する、ものとすることもできる。この場合、前記複数の検出比と、前記第１状態での前記検出部の位置における前記第１方向範囲の各方向の前記測定対象の通過経路長に基づく該第１方向範囲の各方向の前記高透過性粒子のフラックスと、前記第２状態での前記検出部の位置における前記第２方向範囲の各方向の前記測定対象の通過経路長に基づく該第２方向範囲の各方向の前記高透過性粒子のフラックスと、のうち同一方向の前記高透過性粒子のフラックスの比である複数の解析比と、を用いて前記測定対象の密度関連物理量を演算する、ものとすることもできる。ここで、「平行移動」は、第１状態か第２状態かに拘わらず検出方向範囲が同一となる移動を意味する。

本発明のプログラムは、
コンピュータを、所定の高透過性を有する高透過性粒子の通過を検出可能な検出部を有する検出装置が測定対象の周辺に配置されたときの前記検出部による検出結果を用いて前記測定対象の内部の状態を解析する装置として機能させるプログラムであって、
前記検出装置が前記測定対象の周辺に配置された第１状態での前記検出部による検出結果に応じた前記高透過性粒子のフラックスである第１状態検出フラックスと、前記検出部で前記高透過性粒子の通過を検出可能な該検出部から見た方向の範囲である検出方向範囲が前記第１状態から変更された第２状態での前記検出部による検出結果に応じた前記高透過性粒子のフラックスである第２状態検出フラックスと、の比を用いて前記測定対象の密度関連物理量を演算するモジュール、
を備えることを特徴とする。

この本発明のプログラムでは、所定の高透過性を有する高透過性粒子の通過を検出可能な検出部を有する検出装置が測定対象の周辺に配置された第１状態での検出部による検出結果に応じた高透過性粒子のフラックスである第１状態検出フラックスと、検出部で高透過性粒子の通過を検出可能な検出部から見た方向の範囲である検出方向範囲が第１状態から変更された第２状態での検出部による検出結果に応じた高透過性粒子のフラックスである第２状態検出フラックスと、の比を用いて測定対象の密度関連物理量を演算する。即ち、第１状態検出フラックスと第２状態検出フラックスとの比を用いて測定対象の密度関連物理量を演算するのである。これにより、検出部の各部分のバラツキなどによる影響を抑制することができ、測定対象の内部の状態をより適正に解析することができる。ここで、「高透過性粒子」は、ミュオンである、ものとすることもできる。また、「密度関連物理量」は、密度（方位角方向や仰角方向の平均密度を含む）や、密度分布（上限密度，下限密度）である、ものとすることもできる。

本発明の内部状態解析システムは、
所定の高透過性を有する高透過性粒子の通過を検出可能な検出部を有する検出装置と、前記検出装置が測定対象の周辺に配置されたときの前記検出部による検出結果を用いて前記測定対象の内部の状態を解析する解析手段と、を備える内部状態解析システムであって、
前記解析手段は、前記検出装置が前記測定対象の周辺に配置された第１状態での前記検出部による検出結果に応じた前記高透過性粒子のフラックスである第１状態検出フラックスと、前記検出部で前記高透過性粒子の通過を検出可能な該検出部から見た方向の範囲である検出方向範囲が前記第１状態から変更された第２状態での前記検出部による検出結果に応じた前記高透過性粒子のフラックスである第２状態検出フラックスと、の比を用いて前記測定対象の密度関連物理量を演算する手段である、
ことを特徴とする。

この本発明の内部状態解析システムでは、所定の高透過性を有する高透過性粒子の通過を検出可能な検出部を有する検出装置が測定対象の周辺に配置された第１状態での検出部による検出結果に応じた高透過性粒子のフラックスである第１状態検出フラックスと、検出部で高透過性粒子の通過を検出可能な検出部から見た方向の範囲である検出方向範囲が第１状態から変更された第２状態での検出部による検出結果に応じた高透過性粒子のフラックスである第２状態検出フラックスと、の比を用いて測定対象の密度関連物理量を演算する。即ち、第１状態検出フラックスと第２状態検出フラックスとの比を用いて測定対象の密度関連物理量を演算するのである。これにより、検出部の各部分のバラツキなどによる影響を抑制することができ、測定対象の内部の状態をより適正に解析することができる。ここで、「高透過性粒子」は、ミュオンである、ものとすることもできる。また、「密度関連物理量」は、密度（方位角方向や仰角方向の平均密度を含む）や、密度分布（上限密度，下限密度）である、ものとすることもできる。

本発明の一実施例としての内部状態解析システム２０の構成の概略を示す構成図である。 内部状態解析システム２０の検出装置２２の構成の概略を示す構成図である。 検出装置２２（ホドスコープ２４）や測定対象を天頂から見た様子を模式的に示す模式図である。 実施例の内部状態解析ルーチンの一例を示すフローチャートである。 平均密度計算処理の一例を示すフローチャートである。 第１密度長起因フラックス計算処理の一例を示すフローチャートである。 上下限密度比計算処理の一例を示すフローチャートである。 第２密度長起因フラックス計算処理の一例を示すフローチャートである。 障害物（水相当）の厚さとその障害物を通過するミュオンのフラックスとの関係の一例を示す説明図である。 ホドスコープ２４でミュオンの通過を検出可能な方向範囲（第１方向範囲Ｒｄｉや第２方向範囲Ｒｄｉ２）の各方位角φ，各仰角θのホドスコープ２４の検出効率ηの一例を示す説明図である。 ミュオンのエネルギとパラメータａ，ｂと通過経路長ＤＬとの関係の一例を示すテーブルである。 実施例の内部状態解析ルーチンの一例を示すフローチャートである。 密度勾配計算処理の一例を示すフローチャートである。 検出装置２２（ホドスコープ２４）や測定対象を天頂から見た様子を模式的に示す模式図である。

次に、本発明を実施するための形態を実施例を用いて説明する。

図１は、本発明の一実施例として、所定の高透過性を有する高透過性粒子としてのミュオン（ミュー粒子）の性質（障害物を通過する性質）を用いて測定対象（例えば、山や産業プラント，ダムなど）の内部の状態を解析する内部状態解析システム２０の構成の概略を示す構成図あり、図２は、内部状態解析システム２０の検出装置２２の構成の概略を示す構成図である。実施例の内部状態解析システム２０は、図１に示すように、測定対象の周辺に配置されてミュオンの通過を検出する検出装置２２と、検出装置２２による検出結果を用いてミュオンの数をカウントするカウント装置４０と、カウント装置４０によるカウント結果を用いて測定対象の内部の状態を解析する解析装置５０と、を備える。

検出装置２２は、図２に示すように、上面が方位角方向に回転可能な回転架台２３と、回転架台２３の上面に固定されてミュオンを検出するホドスコープ２４と、を備える。ホドスコープ２４は、略平板状の検出器２６，２８が若干の間隔を開けて平行となる（対向する）よう回転架台２３に固定されて構成されている。検出器２６，２８は、それぞれ、略四角柱のシンチレータ３１とシンチレータ３１の長手方向の端部に取り付けられた光電子増倍管３２とを有するセンサモジュール３０を複数用いて構成されており、１０個のセンサモジュール３０が配列されたセンサモジュール群２６ａ，２８ａと１０個のセンサモジュール３０が配列されたセンサモジュール群２６ｂ，２８ｂとが互いの配列方向が略直交するよう隣接配置されて全体として略平板状となるよう構成されている。各センサモジュール３０では、シンチレータ３１にミュオンが入射されると、シンチレータ３１内でフォトンが発生し、そのフォトンが光電子増倍管３２で電子に変換され更に増幅されて、カウント装置４０に比較的大きな振幅（波高）の信号を出力する。なお、実施例では、１ｍ×１０ｃｍ×２ｃｍなどのシンチレータ３１を用いると共に検出器２６と検出器２８との間隔を１ｍなどとするものとした。即ち、ホドスコープ２４は、測定対象に対して十分に小さいとみなすことができるものとした。

カウント装置４０は、検出器２６，２８からの出力（信号）に対してノイズを除去して出力するディスクリミネータ４２と、ディスクリミネータ４２からの出力（信号）のうち検出器２６，２８で同時にミュオンの通過を検出したとみなされる信号を選択して出力するコインシデンス検出部４４と、コインシデンス検出部４４からの出力（信号）を用いてホドスコープ２４から見た各方向（方位角φ，仰角θ）のミュオンの数を検出データＯｂｓ［φ，θ］としてカウントするカウント部４６と、を備える。ここで、方位角φは、ホドスコープ２４から見て北方向を基準（０°）として時計回りに正の角度を定めるものとし、仰角θは、ホドスコープ２４から見て水平方向を基準（０°）として上向き（天頂）側に正の角度を定めるものとした。

コインシデンス検出部４４は、ディスクリミネータ４２からの出力（信号）を用いて、センサモジュール群２６ａ，２６ｂのコインシデンスを取ったり（検出器２６におけるミュオンの通過位置α１を特定したり）、センサモジュール群２８ａ，２８ｂのコインシデンスを取ったり（検出器２８におけるミュオンの通過位置α２を特定したり）、センサモジュール群２６ａ，２６ｂ，２８ａ，２８ｂのコインシデンスを取ったり（通過位置α１，α２に対応する方位角φ，仰角θを特定したり）する。実施例では、検出器２６，２８の通過位置α１，α２をそれぞれ［ｘ１，ｙ１］，［ｘ２，ｙ２］（ただし、ｘ１，ｘ２については図２中左側から右側に向けて１〜１０とし、ｙ１，ｙ２については図２中下側から上側向けて１〜１０とする）として方位角φ，仰角θを特定するものとした。例えば、図２のミュオンの場合には、通過位置α１［８，９］，通過位置α２［４，４］に対応する方位角φ，仰角θを特定するものとした。

解析装置５０は、図１に示すように、汎用のコンピュータ５２にアプリケーションソフトウエアとしての内部状態解析プログラム６０がインストールされたものとして構成されている。コンピュータ５２は、図示しないＣＰＵやＲＯＭ，ＲＡＭ，グラフィックプロセッサ（ＧＰＵ），グラフィックメモリ（ＶＲＡＭ），システムバス，ハードディスクドライブ（ＨＤＤ）などを備え、ハードディスクドライブに内部状態解析プログラム６０などが記憶されている。内部状態解析プログラム６０は、データを入力する入力モジュール６２と、入力されたデータを用いて測定対象の内部の状態を解析する解析モジュール６４と、解析結果を出力する出力モジュール６６と、によって構成されている。なお、コンピュータ５０には、表示装置としてのディスプレイ７０や、入力装置としてのキーボード７２やマウス７４などが接続されている。

次に、こうして構成された実施例の内部状態解析システム２０の動作、特に、測定対象の内部状態を解析する処理について説明する。なお、実施例では、測定対象の内部の状態を解析する前に以下の準備処理を実行するものとした。

図３は、検出装置２２（ホドスコープ２４）や測定対象を天頂から見た様子を模式的に示す模式図である。図中、「Ｌ［φ，θ］」は、ホドスコープ２４から見た各方向（方位角φ，仰角θ）の測定対象の通過経路の長さ（以下、通過経路長という）を示す。この通過経路長Ｌ［φ，θ］は、山などを測定対象とする場合には地図などを用いて求めることができ、産業プラントやダムなどを測定対象とする場合には設計図などを用いて求めることができる。また、図中、実線のホドスコープ２４は、ホドスコープ２４を方位角方向に回転させる前（第１状態）を示し、点線のホドスコープ２４は、ホドスコープ２４を方位角方向に回転させた後（第２状態）を示す。

準備処理としては、まず、検出装置２２（ホドスコープ２４）を測定対象の周辺（例えば、長手方向に数ｋｍ程度の山などを測定対象とする場合には測定対象の中心から測定対象の短手方向に数百ｍ程度離れた位置など）に配置して（第１状態として）、第１方位角範囲Ｒａｚ１（φ：φ１ｍｉｎ〜φ１ｍａｘ）（図３参照），仰角範囲Ｒｅｌ（θ：θｍｉｎ〜θｍａｘ）（図示せず）からなる第１方向範囲Ｒｄｉ１のミュオンを検出できるようにする。ここで、下限方位角φ１ｍｉｎは、第１状態において検出器２６，２８の通過位置α１，α２のｘ座標がｘ１＝１，ｘ２＝１０となるミュオンの進行方向に相当する方位角φであり、上限方位角φ１ｍａｘは、第１状態において検出器２６，２８の通過位置α１，α２のｘ座標がｘ１＝１０，ｘ２＝１となるミュオンの進行方向に相当する方位角φである。また、下限仰角θｍｉｎは、第１状態および後述の第２状態において検出器２６，２８の通過位置α１，α２のｙ座標がｙ１＝ｙ２となるミュオンの進行方向に相当する仰角θであり、上限仰角θｍａｘは、第１状態および後述の第２状態において検出器２６，２８の通過位置α１，α２のｙ座標がｙ１＝１０，ｙ２＝１となるミュオンの進行方向に相当する仰角θである。

そして、その第１状態で、ホドスコープ２４によって期間Ｔ１（例えば、数週間〜数ヶ月程度）に亘ってミュオンを検出すると共に、ホドスコープ２４による検出結果を用いてカウント装置４０によって第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θのミュオンの数を検出データＯｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ）としてカウント（検出）する。

その後、ホドスコープ２４を第１状態から方位角方向に所定角度Δφ（例えば、数度〜十数度程度）だけ回転させて（第２状態として）、第２方位角範囲Ｒａｚ２（φ＋Δφ：φ２ｍｉｎ（＝φ１ｍｉｎ＋Δφ）〜φ２ｍａｘ（＝φ１ｍａｘ＋Δφ））（図３参照），仰角範囲Ｒｅｌ（θ：θｍｉｎ〜θｍａｘ）（図示せず）からなる第２方向範囲Ｒｄｉ２のミュオンを検出できるようにする。ここで、下限方位角φ２ｍｉｎは、第２状態において検出器２６，２８の通過位置α１，α２のｘ座標がｘ１＝１，ｘ２＝１０となるミュオンの進行方向に相当する方位角φであり、上限方位角φ２ｍａｘは、第２状態において検出器２６，２８の通過位置α１，α２のｘ座標がｘ１＝１０，ｘ２＝１となるミュオンの進行方向に相当する方位角φである。また、上述したように、下限仰角θｍｉｎ，上限仰角θｍａｘは、第１状態と第２状態とで同一である。

そして、その第２状態で、ホドスコープ２４によって期間Ｔ２（例えば、数週間〜数ヶ月程度）に亘ってミュオンを検出すると共に、ホドスコープ２４による検出結果を用いてカウント装置４０によって第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θのミュオンの数を検出データＯｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ２ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）としてカウント（検出）する。

ここで、第２状態は、ホドスコープ２４を第１状態から方位角方向に所定角度Δφだけ回転させた状態であるから、第１状態における方向［φ，θ］と第２状態における方向［φ＋Δφ，θ］とは、ホドスコープ２４の検出器２６，２８における同一のセンサモジュール３０の組み合わせ又はそれと同一方向のセンサモジュール３０の組み合わせでミュオンを検出したものとして考えることができる。具体的には、同一のセンサモジュール３０の組み合わせとは、例えば、第１状態における通過位置α１［８，９］，α２［４，４］の組み合わせと第２の状態における通過位置α１［８，９］，α２［４，４］の組み合わせとをいい、同一方向のセンサモジュール３０の組み合わせとは、例えば、第１状態における通過位置α１［８，９］，α２［４，４］の組み合わせと、第２状態における通過位置α１［６，９］，α２［２，４］の組み合わせや通過位置α１［１０，９］，α２［６，４］の組み合わせ，通過位置α１［８，７］，α２［４，２］の組み合わせ，通過位置α１［１０，７］，α２［６，２］の組み合わせなどをいう。

次に、測定対象の内部の状態を解析する処理について説明する。図４は、測定対象の内部の状態の解析に用いられる内部状態解析ルーチンの一例を示すフローチャートである。このルーチンは、ユーザによって内部状態解析プログラム６０の実行が指示されたときに、ＣＰＵにより、ＲＡＭの所定アドレスに書き込まれて実行される。以下の説明では、第１方位角範囲Ｒａｚ１（φ１ｍｉｎ〜φ１ｍａｘ）と第２方位角範囲Ｒａｚ２（φ２ｍｉｎ〜φ２ｍａｘ）とからなる方位角φの範囲を全体方位角範囲Ｒａｚｔｏ（φ１ｍｉｎ〜φ２ｍａｘ）（図３参照）といい、第１方向範囲Ｒｄｉ１（φ１ｍｉｎ〜φ１ｍａｘ，θｍｉｎ〜θｍａｘ）と第２方向範囲Ｒｄｉ２（φ２ｍｉｎ〜φ２ｍａｘ，θｍｉｎ〜θｍａｘ）とからなる方向の範囲を全体方向範囲Ｒｄｉｔｏ（φ１ｍｉｎ〜φ２ｍａｘ，θｍｉｎ〜θｍａｘ）という。なお、仰角範囲Ｒｅｌ（θｍｉｎ〜θｍａｘ）は、第１状態と第２状態とで同一である。

内部状態解析プログラム６０が実行されると、ＣＰＵは、まず、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ、θの範囲については以下同様）や、第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ２ｍｉｎ〜φ２ｍａｘ），全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの測定対象の通過経路長Ｌ［φ，θ］（φ：φ１ｍｉｎ〜φ２ｍａｘ）などのデータを入力する処理を実行する（ステップＳ１００）。ここで、検出データＯｂｓ１［φ，θ］や検出データＯｂｓ２［φ＋Δφ，θ］は、インターネットなどの外部ネットワークを介してカウント装置４０から受信して用いたり、カウント装置４０から図示しない記憶媒体（例えば、ハードディスクやＣＤ，ＤＶＤなど）に記憶させた後にその記憶媒体を解析装置５０に接続して読み出して用いたりすることができる。また、通過経路長Ｌ［φ，θ］は、地図データや設計図のデータなどを読み出して用いたりすることができる。

こうしてデータを入力すると、次式（１）に示すように、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］をそれぞれ第１状態の期間Ｔ１で規格化する（除する）ことにより、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θのミュオンのフラックス（以下、検出フラックスという）Ｆｏｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）を計算すると共に、式（２）に示すように、第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］をそれぞれ第２状態の期間Ｔ２で規格化する（除する）ことにより、第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出フラックスＦｏｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ１ｍｉｎ〜φ２ｍａｘ）を計算する（ステップＳ１１０）。

続いて、式（３）に示すように、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出フラックスＦｏｂｓ１［φ，θ］をそれぞれ方位角方向に所定角度Δφだけ異なる方向の検出フラックスＦｏｂｓ２［φ＋Δφ，θ］で除することにより、各仰角θにおける各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）を計算する（ステップＳ１２０）。この処理は、例えば、第１方向範囲Ｒｄｉ１の下限方位角φ１ｍｉｎ，下限仰角θｍｉｎの検出フラックスＦｏｂｓ１［φ１ｍｉｎ，θｍｉｎ］を、第２方向範囲Ｒｄｉ２の下限方位角φ２ｍｉｎ，下限仰角θｍｉｎの検出フラックスＦｏｂｓ２［φ２ｍｉｎ，θｍｉｎ］で除することにより、下限仰角θｍｉｎにおける方位角φ１ｍｉｎ，φ１ｍｉｎ＋Δφ間の検出比Ｒｏｂｓ［φ１ｍｉｎ，Δφ，θｍｉｎ］を計算する処理である。

そして、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］と第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］とが共にガウス分布に従うと仮定して、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］と第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］と各仰角θにおける各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ］とを用いて、次式（４）により、各仰角θにおける各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ］についての統計誤差ＥＲｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）を計算する（ステップＳ１３０）。

次に、仰角範囲Ｒｅｌのうち解析の対象とする仰角θとしての対象仰角θ＊に下限仰角θｍｉｎを設定し（ステップＳ１４０）、設定した対象仰角θ＊における測定対象の方位角方向（全体方位角範囲Ｒａｚｔｏ）の平均密度ρｍ［θ＊］を図５に例示する平均密度計算処理により計算する（ステップＳ１５０）。以下、図４の内部状態解析ルーチンの説明を一旦中断し、図５の平均密度計算処理について説明する。

図５の平均密度計算処理では、まず、対象仰角θ＊における測定対象の方位角方向（全体方位角範囲Ｒａｚｔｏ）の平均密度ρｍ［θ＊］に初期値ρｍ０（例えば、０．０ｇ／ｃｃなど）を設定し（ステップＳ３００）、次式（５）に示すように、全体方位角範囲Ｒａｚｔｏの各方位角φの測定対象の通過経路長Ｌ［φ，θ＊］（φ：φ１ｍｉｎ〜φ２ｍａｘ）のそれぞれに方位角方向の平均密度ρｍ［θ＊］を乗じたものを、全体方位角範囲Ｒａｚｔｏの各方位角φの測定対象の第１密度長ＤＬ１［φ，θ＊］（φ：φ１ｍｉｎ〜φ２ｍａｘ）として計算する（ステップＳ３１０）。そして、計算した全体方位角範囲Ｒａｚｔｏの各方位角φの測定対象の第１密度長ＤＬ１［φ，θ＊］を用いて、全体方位角範囲Ｒａｚｔｏの各方位角φのミュオンのフラックス（以下、第１密度長起因フラックスという）Ｎ１［φ，θ＊］（φ：φ１ｍｉｎ〜φ２ｍａｘ）を図６に例示する第１密度長起因フラックス計算処理により計算する（ステップＳ３２０）。以下、図５の平均密度計算処理の説明を一旦中断し、図６の第１密度長起因フラックス計算処理について説明する。

図６の第１密度長起因フラックス計算処理では、まず、全体方位角範囲Ｒａｚｔｏの各方位角φの測定対象の第１密度長ＤＬ１［φ，θ＊］（φ：φ１ｍｉｎ〜φ２ｍａｘ）を用いて、ミュオンのカットオフエネルギ（以下、第１カットオフエネルギという）Ｅｃｕｔ１［φ，θ＊］の計算に用いるパラメータａ１［ＤＬ１［φ，θ＊］］，ｂ１［ＤＬ１［φ，θ＊］］（φ：φ１ｍｉｎ〜φ２ｍａｘ）を設定する（ステップＳ４００）。パラメータａ１［ＤＬ１［φ，θ＊］］，ｂ１［ＤＬ１［φ，θ＊］］の設定は、実施例では、第１密度長ＤＬ１［φ，θ＊］とパラメータａ１［ＤＬ１［φ，θ＊］］，ｂ１［ＤＬ１［φ，θ＊］］との予め定められた関係に第１密度長ＤＬ１［φ，θ＊］を適用してパラメータａ１［ＤＬ１［φ，θ＊］］，ｂ１［ＤＬ１［φ，θ＊］］を設定するものとした。例えば、第１密度長ＤＬ１［φ，θ＊］が４９２０ｇ／ｃｍ２のときには、パラメータａ１［ＤＬ１［φ，θ＊］］，ｂ１［ＤＬ１［φ，θ＊］］にはそれぞれ２．１７０ＭｅＶ・ｃｍ２／ｇ，０．０１９ＭｅＶ・ｃｍ２／ｇを設定し、第１密度長ＤＬ１［φ，θ＊］が５５３.４ｇ／ｃｍ２のときには、パラメータａ１［ＤＬ１［φ，θ＊］］，ｂ１［ＤＬ１［φ，θ＊］］にはそれぞれ１．８０８ＭｅＶ・ｃｍ２／ｇ，０．０００ＭｅＶ・ｃｍ２／ｇを設定することができる。

続いて、設定したパラメータａ１［ＤＬ１［φ，θ＊］］，ｂ１［ＤＬ１［φ，θ＊］］を用いて、次式（６）により、全体方位角範囲Ｒａｚｔｏの各方位角φの第１カットオフエネルギＥｃｕｔ１［φ，θ＊］（φ：φ１ｍｉｎ〜φ２ｍａｘ）を計算する（ステップＳ４１０）。

そして、計算した全体方位角範囲Ｒａｚｔｏの各方位角φの第１カットオフエネルギＥｃｕｔ１［φ，θ＊］を用いて、次式（７）により、全体方位角範囲Ｒａｚｔｏの各方位角φの第１密度長起因フラックスＮ１［φ，θ＊］（φ：φ１ｍｉｎ〜φ２ｍａｘ）を計算して（ステップＳ４２０）、本ルーチンを終了する。ここで、式（７）中、「Ｅ」はミュオンのエネルギであり、「θｚ＊」は対象仰角θ＊に対応する天頂角（π／２−θ＊）である。また、式（７）中、「ｄＮ［Ｅ，θｚ＊］／ｄＥ」は式（８）により表わすことができる。式（８）中、「γ」は値２．７を用いることができ、「Ａ」は式（９）により表わすことができ、「Ｗ」は式（１０）により表わすことができる。式（８）および式（１０）中、「ΔＥ」は式（１１）により表わすことができる。式（１０）および式（１１）中、「Ｌ」は式（１２）により表わすことができる。

以上、図６の第１密度長起因フラックス計算処理について説明した。図５の平均密度計算処理の説明に戻る。ステップＳ３２０で全体方位角範囲Ｒａｚｔｏの各方位角φの第１密度長起因フラックスＮ１［φ，θ＊］（φ：φ１ｍｉｎ〜φ２ｍａｘ）を計算すると、次式（１３）に示すように、第１方位角範囲Ｒａｚ１の各方位角φの第１密度長起因フラックスＮ１［φ，θ＊］（φ：φ１ｍｉｎ〜φ１ｍａｘ）をそれぞれ対応する第１密度長起因フラックスＮ１［φ＋Δφ，θ＊］で除することにより、全体方位角範囲Ｒａｚｔｏの各方位角φ，φ＋Δφ間の第１解析比Ｒｅｓｔ１［φ，Δφ，θ＊］（φ：φ１ｍｉｎ〜φ１ｍａｘ）を計算する（ステップＳ３３０）。この処理は、例えば、方位角φ１ｍｉｎの第１密度長起因フラックスＮ１［φ１ｍｉｎ，θ＊］を方位角φ２ｍｉｎ（＝φ１ｍｉｎ＋Δφ）の第１密度長起因フラックスＮ１［φ２ｍｉｎ，θ＊］で除することにより、方位角φ１ｍｉｎ，φ１ｍｉｎ＋Δφ間の第１解析比Ｒｅｓｔ１［φ１ｍｉｎ，Δφ，θ＊］を計算する処理である。

続いて、次式（１４）に示すように、全体方位角範囲Ｒａｚｔｏの各方位角φ，φ＋Δφ間の第１解析比Ｒｅｓｔ１［φ，Δφ，θ＊］（φ：φ１ｍｉｎ〜φ１ｍａｘ）の方位角方向についての総和（方位角φについての総和）を解析比和Ｓｅｓｔ［θ＊］として計算し（ステップＳ３４０）、式（１５）に示すように、ステップＳ１２０で計算した全体方位角範囲Ｒａｚｔｏの各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ＊］（φ：φ１ｍｉｎ〜φ１ｍａｘ）の方位角方向についての総和（方位角φについての総和）を検出比和Ｓｏｂｓ［θ＊］として計算する（ステップＳ３５０）。

そして、ステップＳ１３０で計算した全体方位角範囲Ｒａｚｔｏの各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ］についての統計誤差ＥＲｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）を用いて、次式（１６）により、検出比和Ｓｏｂｓ［θ＊］についての統計誤差ＥＳｏｂｓ［θ＊］を計算する（ステップＳ３６０）。

次に、検出比和Ｓｏｂｓ［θ＊］と解析比和Ｓｅｓｔ［θ＊］との差分が検出比和Ｓｏｂｓ［θ＊］についての統計誤差ＥＳｏｂｓ［θ＊］未満であるか否かを判定する（ステップＳ３７０）。この処理は、現在の平均密度ρｍ［θ＊］が適正であるか否かを判定する処理である。

検出比和Ｓｏｂｓ［θ＊］と解析比和Ｓｅｓｔ［θ＊］との差分が検出比和Ｓｏｂｓ［θ＊］についての統計誤差ＥＳｏｂｓ［θ＊］以上であると判定されたときには、現在の平均密度ρｍ［θ＊］は適正でないと判断し、現在の平均密度ρｍ［θ＊］に所定値ｄρｍ０（例えば、０．０１ｇ／ｃｃなど）を加えたものを新たな平均密度ρｍ［θ＊］に設定して（ステップＳ３７０）、ステップＳ３１０に戻る。

こうしてステップＳ３１０〜Ｓ３８０の処理を繰り返し実行して、ステップＳ３７０で検出比和Ｓｏｂｓ［θ＊］と解析比和Ｓｅｓｔ［θ＊］との差分が検出比和Ｓｏｂｓ［θ＊］についての統計誤差ＥＳｏｂｓ［θ＊］未満であると判定されたときには、現在の平均密度ρｍ［θ＊］は適正であると判断し、それを対象仰角θ＊における平均密度ρｍ［θ＊］として確定して（ステップＳ３９０）、本ルーチンを終了する。

以上、図５の平均密度計算処理について説明した。図４の内部状態解析ルーチンの説明に戻る。ステップＳ１５０で対象仰角θ＊における測定対象の方位角方向（全体方位角範囲Ｒａｚｔｏ）の平均密度ρｍ［θ＊］を計算すると、計算した平均密度ρｍ［θ＊］を全体方位角範囲Ｒａｚｔｏ（φ１ｍｉｎ〜φ２ｍａｘ）の中心方位角φｃの密度ρ［φｃ，θ＊］に設定する（ステップＳ１６０）。

続いて、第１方位角範囲Ｒａｚ１のうち解析の対象とする方位角φとしての対象方位角φ＊に下限方位角φ１ｍｉｎを設定し（ステップＳ１７０）、対象仰角θ＊における方位角φ＊，φ＊＋Δφ間の密度比の上下限としての上下限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］，ρｆｍｉｎ［φ＊，Δφ，θ＊］を図７に例示する上下限密度比計算処理により計算する（ステップＳ１８０）。以下、図４の内部状態解析ルーチンの説明を一旦中断し、図７の上下限密度比計算処理について説明する。

図７の上下限密度比計算処理では、まず、対象仰角θ＊における方位角φ＊，φ＊＋Δφの方位角方向の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］にある程度大きな正の初期値Δρ０（例えば、０．５０ｇ／ｃｃなど）を設定し（ステップＳ５００）、設定した密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］とステップＳ１６０で設定した中心方位角φｃの密度ρ［φｃ，θ＊］とを用いて、次式（１７），（１８）により、方位角φ＊，φ＊＋Δφの密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］を計算する（ステップＳ５１０）。

続いて、次式（１９），（２０）に示すように、方位角φ＊，φ＊＋Δφの通過経路長Ｌ［φ＊，θ＊］，Ｌ［φ＊＋Δφ，θ＊］に方位角φ＊，φ＊＋Δφの密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］を乗じたものを方位角φ＊，φ＊＋Δφの測定対象の第２密度長ＤＬ２［φ＊，θ＊］，ＤＬ２［φ＊＋Δφ，θ＊］として計算し（ステップＳ５２０）、計算した方位角φ＊，φ＊＋Δφの測定対象の第２密度長ＤＬ２［φ＊，θ＊］，ＤＬ２［φ＊＋Δφ，θ＊］を用いて、方位角φ＊，φ＊＋Δφのミュオンのフラックス（以下、第２密度長起因フラックスという）Ｎ２［φ＊，θ＊］，Ｎ２［φ＊＋Δφ，θ＊］を図８に例示する第２密度長起因フラックス計算処理により計算する（ステップＳ５３０）。以下、図７の上下限密度比計算処理の説明を一旦中断し、図８の第２密度長起因フラックス計算処理について説明する。

図８の第２密度長起因フラックス計算処理では、まず、方位角φ＊，φ＊＋Δφの測定対象の第２密度長ＤＬ２［φ＊，θ＊］，ＤＬ２［φ＊＋Δφ，θ＊］を用いてミュオンのカットオフエネルギ（以下、第２カットオフエネルギという）Ｅｃｕｔ２［φ＊，θ＊］，Ｅｃｕｔ２［φ＊＋Δφ，θ＊］の計算に用いるパラメータａ２［ＤＬ２［φ＊，θ＊］］，ｂ２［ＤＬ２［φ＊，θ＊］］，ａ２［ＤＬ２［φ＊＋Δφ，θ＊］］，ｂ２［ＤＬ２［φ＊＋Δφ，θ＊］］を設定する（ステップＳ７００）。パラメータａ２［ＤＬ２［φ＊，θ＊］］，ｂ２［ＤＬ２［φ＊，θ＊］］，ａ２［ＤＬ２［φ＊＋Δφ，θ＊］］，ｂ２［ＤＬ２［φ＊＋Δφ，θ＊］］の設定は、図６の第１密度長起因処理のステップＳ４００の処理と同様に行なうことができる。

続いて、設定したパラメータａ２［ＤＬ２［φ＊，θ＊］］，ｂ２［ＤＬ２［φ＊，θ＊］］，ａ２［ＤＬ２［φ＊＋Δφ，θ＊］］，ｂ２［ＤＬ２［φ＊＋Δφ，θ＊］］を用いて、次式（２１），（２２）により、方位角φ＊，φ＊＋Δφの第２カットオフエネルギＥｃｕｔ２［φ＊，θ＊］，Ｅｃｕｔ２［φ＊＋Δφ，θ＊］を計算する（ステップＳ７１０）。

そして、方位角φ＊，φ＊＋Δφの第２カットオフエネルギＥｃｕｔ２［φ＊，θ＊］，Ｅｃｕｔ２［φ＊＋Δφ，θ＊］を用いて、次式（２３），（２４）により、方位角φ＊，φ＊＋Δφの密度長起因フラックスＮ２［φ＊，θ＊］，Ｎ２［φ＊＋Δφ，θ＊］を計算して（ステップＳ７２０）、本ルーチンを終了する。ここで、式（２３），（２４）は、それぞれ、上述の式（７）の積分開始の値を「Ｅ＝Ｅｃｕｔ１［φ，θ＊］」から「Ｅ＝Ｅｃｕｔ２［φ＊，θ＊］」，「Ｅ＝Ｅｃｕｔ２［φ＊＋Δφ，θ＊］」に置き換えた点を除いて、式（７）と同一である。

以上、図８の第２密度長起因フラックス計算処理について説明した。図７の上下限密度比計算処理の説明に戻る。ステップＳ５３０で方位角φ＊，φ＊＋Δφの第２密度長起因フラックスＮ２［φ＊，θ＊］，Ｎ２［φ＊＋Δφ，θ＊］を計算すると、次式（２５）に示すように、第２密度長起因フラックスＮ２［φ＊，θ＊］を対応する第２密度長起因フラックスＮ２［φ＊＋Δφ，θ＊］で除することにより、方位角φ＊，φ＊＋Δφ間の第２解析比Ｒｅｓｔ２［φ＊，Δφ，θ＊］を計算する（ステップＳ５４０）。

続いて、方位角φ＊，φ＊＋Δφ間の上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］を設定済みであるか否かを判定し（ステップＳ５５０）、上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］を設定済みでないと判定されたときには、方位角φ＊，φ＊＋Δφ間の第２解析比Ｒｅｓｔ２［φ＊，Δφ，θ］が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］より大きく検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］と検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］についての統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］との和未満であるか否かを判定する（ステップＳ５６０）。この処理は、上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］の設定において現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］や密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］が適正であるか否かを判定する処理である。

第２解析比Ｒｅｓｔ２［φ＊，Δφ，θ］が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］以下であると判定されたときや、第２解析比Ｒｅｓｔ２［φ＊，Δφ，θ］が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］と統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］との和以上であると判定されたときには、上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］の設定において現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］や密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］は適正でないと判断し、現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］から所定値ｄΔρ（例えば、０．０１ｇ／ｃｃなど）を減じたものを新たな密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］に設定して（ステップＳ５７０）、ステップＳ５１０に戻る。

こうしてステップＳ５１０〜Ｓ５７０の処理を繰り返し実行して、ステップＳ５６０で第２解析比Ｒｅｓｔ２［φ＊，Δφ，θ］が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］より大きく検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］と統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］との和未満であると判定されたときには、上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］の設定において現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］や密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］は適正であると判断し、次式（２６）に示すように、密度ρ［φ＊，θ＊］を密度ρ［φ＊＋Δφ，θ＊］で除することにより、方位角φ＊，φ＊＋Δφ間の上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］を計算して（ステップＳ５８０）、ステップＳ５１０に戻る。

こうして方位角φ＊，φ＊＋Δφ間の上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］を設定した後には、ステップＳ５５０で上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］を設定済みであると判定され、方位角φ＊，φ＊＋Δφ間の解析比Ｒｅｓｔ２［φ＊，Δφ，θ］が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］から検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］についての統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］を減じた値より大きく検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］未満であるか否かを判定する（ステップＳ５９０）。この処理は、下限密度比ρｆｍｉｎ［φ＊，Δφ，θ＊］の設定において現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］や密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］が適正であるか否かを判定する処理である。

解析比Ｒｅｓｔ２［φ＊，Δφ，θ］が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］から統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］を減じた値以下または検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］以上であると判定されたときには、下限密度比ρｆｍｉｎ［φ＊，Δφ，θ＊］の設定において現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］や密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］は適正でないと判断し、現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］から所定値ｄΔρ（例えば、０．０１ｇ／ｃｃなど）を減じたものを新たな密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］に設定して（ステップＳ６００）、ステップＳ５１０に戻る。

こうしてステップＳ５１０〜Ｓ５５０，Ｓ５９０，Ｓ６００の処理を繰り返し実行して、ステップＳ５９０で解析比Ｒｅｓｔ２［φ＊，Δφ，θ］が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］から統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］を減じた値より大きく検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］未満であると判定されたときには、下限密度比ρｆｍｉｎ［φ＊，Δφ，θ＊］の設定において現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］や密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］は適正であると判断し、次式（２７）に示すように、密度ρ［φ＊，θ＊］を密度ρ［φ＊＋Δφ，θ＊］で除することにより方位角φ＊，φ＊＋Δφ間の下限密度比ρｆｍｉｎ［φ＊，Δφ，θ＊］を計算して（ステップＳ６１０）、本ルーチンを終了する。

以上、図７の上下限密度比計算処理について説明した。図４の内部状態解析ルーチンの説明に戻る。ステップＳ１８０で対象仰角θ＊における対象方位角φ＊と方位角φ＊＋Δ＊との上下限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］，ρｆｍｉｎ［φ＊，Δφ，θ＊］を計算すると、対象仰角φ＊が上限方位角φ１ｍａｘであるか否かを判定し（ステップＳ１９０）、対象仰角φ＊が上限方位角φ１ｍａｘでないと判定されたときには、現在の対象仰角φ＊に方位角方向の角度分解能としての角度ｄφを加えたものを新たな対象仰角φ＊に設定して（ステップＳ２００）、ステップＳ１８０に戻る。ここで、角度ｄφは、例えば、検出器２６，２８のｘ方向の座標ｘ１，ｘ２を用いて得られる「ｘ１−ｘ２」が値１だけ大きくなるのに相当する角度（値−９から値−８になるとき，値−４から値−３になるとき，値１から値２になるときなどのそれぞれに応じた角度）を用いることができる。

こうしてステップＳ１８０〜Ｓ２００の処理を繰り返し実行して、ステップＳ１９０で対象方位角φ＊が上限方位角φ１ｍａｘであると判定されたときには、ステップＳ１６０で設定した中心方位角φｃの密度ρ［φｃ，θ＊］（＝ρｍ［θ＊］）と次式（２８）〜（３１）とを用いて、全体方向範囲Ｒｄｉｔｏの対象仰角θ＊における各方位角φの密度分布の上下限としての上下限密度ρｍａｘ［φ，θ＊］，ρｍｉｎ［φ，θ＊］（φ１ｍｉｎ〜φ２ｍａｘ）を計算する（ステップＳ２１０）。具体的には、式（２８）中、「ρ［φ＋Δφ，θ＊］」，「ρｆｍａｘ［φ，Δφ，θ＊］」にそれぞれ中心方位角φｃの密度ρ［φｃ，θ＊］，上限密度比ρｆｍａｘ［φｃ，Δφ，θ＊］を代入して方位角φｃ−Δφの上限密度ρｍａｘ［φｃ−Δφ，θ＊］を計算すると共に同様に方位角φｃ−２Δφ〜下限方位角φ１ｍｉｎの上限密度ρｍａｘ［φｃ−２Δφ，θ＊］〜ρｍａｘ［φ１ｍｉｎ，θ＊］を計算し、式（２９）中、「ρ［φ，θ＊］」，「ρｆｍａｘ［φ，Δφ，θ＊］」にそれぞれ中心方位角φｃの密度ρ［φｃ，θ＊］，上限密度比ρｆｍａｘ［φｃ，Δφ，θ＊］を代入して方位角φｃ＋Δφの上限密度ρｍａｘ［φｃ＋Δφ，θ＊］を計算すると共に同様に方位角φｃ＋２Δφ〜上限方位角φ２ｍａｘの上限密度ρｍａｘ［φｃ＋２Δφ，θ＊］〜ρｍａｘ［φ２ｍａｘ，θ＊］を計算する。また、式（３０）中、「ρ［φ＋Δφ，θ＊］」，「ρｆｍｉｎ［φ，Δφ，θ＊］」にそれぞれ中心方位角φｃの密度ρ［φｃ，θ＊］，下限密度比ρｆｍｉｎ［φｃ，Δφ，θ＊］を代入して方位角φｃ−Δφの下限密度ρｍｉｎ［φｃ−Δφ，θ＊］を計算すると共に同様に方位角φｃ−２Δφ〜下限方位角φ１ｍｉｎの下限密度ρｍｉｎ［φｃ−２Δφ，θ＊］〜ρｍｉｎ［φ１ｍｉｎ，θ＊］を計算し、式（３１）中、「ρ［φ，θ＊］」，「ρｆｍｉｎ［φ，Δφ，θ＊］」にそれぞれ中心方位角φｃの密度ρ［φｃ，θ＊］，下限密度比ρｆｍｉｎ［φｃ，Δφ，θ＊］を代入して方位角φｃ＋Δφの下限密度ρｍｉｎ［φｃ＋Δφ，θ＊］を計算すると共に同様に方位角φｃ＋２Δφ〜上限方位角φ２ｍａｘの下限密度ρｍｉｎ［φｃ＋２Δφ，θ＊］〜ρｍｉｎ［φ２ｍａｘ，θ＊］を計算する。

次に、対象仰角θ＊が上限仰角θｍａｘであるか否かを判定し（ステップＳ２２０）、対象仰角θ＊が上限仰角θｍａｘでないと判定されたときには、現在の対象仰角θ＊に仰角方向の角度分解能としての角度ｄθを加えたものを新たな対象方位角φ＊に設定して（ステップＳ２３０）、ステップＳ１５０に戻る。ここで、角度ｄθは、例えば、検出器２６，２８のｙ方向の座標ｙ１，ｙ２を用いて得られる「ｙ１−ｙ２」が値１だけ大きくなるのに相当する角度（値１から値２になるとき，値４から値５になるとき，値８から値９になるときなどそれぞれに応じた角度）を用いることができる。こうしてステップＳ１５０〜Ｓ２３０の処理を繰り返し実行して、ステップ２２０で対象仰角θ＊が上限仰角θｍａｘであると判定されたときには、計算結果をディスプレイ７０に出力したりハードディスクドライブに記憶させたりして（ステップＳ２４０）、本ルーチンを終了する。計算結果の出力では、全体方向範囲Ｒｄｉｔｏの各仰角θにおける測定対象の方位角方向の平均密度ρｍ［θ］や、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの測定対象の上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］などを出力する。

ここで、第１状態（ホドスコープ２４を方位角方向に回転させる前）の第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）に応じた検出フラックスＦｏｂｓ１［φ，θ］と、第２状態（ホドスコープ２４を方位角方向に回転させた後）の第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ２ｍｉｎ〜φ２ｍａｘ）に応じた検出フラックスＦｏｂｓ２［φ，θ］とを用いて各仰角θにおける各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）を計算し、計算した複数の検出比Ｒｏｂｓ［φ，Δφ，θ］を用いて、全体方向範囲Ｒｄｉｔｏの各仰角θにおける測定対象の方位角方向の平均密度ρｍ［θ］を計算したり、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの測定対象の上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］を計算したりする理由について説明する。

図９は、障害物（水相当）の厚さとその障害物を通過するミュオンのフラックスとの関係の一例を示す説明図であり、図１０は、ホドスコープ２４でミュオンの通過を検出可能な方向範囲（第１方向範囲Ｒｄｉや第２方向範囲Ｒｄｉ２）の各方位角φ，各仰角θのホドスコープ２４の検出効率η［φ，θ］の一例を示す説明図であり、図１１は、ミュオンのエネルギとパラメータａ，ｂと通過経路長ＤＬとの関係の一例を示すテーブル（参考文献１参照）である。図９中、［ｋｍ・ｗ・ｅ］は水相当の物質の厚さを示し、ハッチングは、ミュオンのフラックスの不定性（バラツキ）を示す。図１１中、「T」はミュオンのエネルギ（運動エネルギ）を示し、「Ionization」（電離）はカットオフエネルギＥｃｕｔの計算に用いるパラメータａを示し、「Brems」（制動輻射）と「Pair prod」（直接電子−陽電子対生成）と「Photonucl」（光核反応）との和はカットオフエネルギＥｃｕｔの計算に用いるパラメータｂを示し、「CSDA range」は障害物（測定対象）の通過経路長ＤＬを示す。
参考文献１：「Table 281: Muons in Standard rock」、［On line］、［平成２４年１月２３日検索］、インターネット
＜http://pdg.lbl.gov/2011/AtomicNuclearProperties/MUON_ELOSS_TABLES/muonloss_281.pdf＞

一般に、ミュオンは、地球磁場による方位角依存性を有すると共に、中間エネルギ領域（図９では０．５ｋｍ・ｗ・ｅ〜３．０ｋｍ・ｗ・ｅ程度を通過するミュオンのエネルギ領域，図１１では２９７ＭｅＶ〜６９３ＧｅＶ程度の領域）より高エネルギ側の領域や低エネルギ側の領域でフラックスに不定性を有する（バラツキが生じやすい）ことが分かっている（図９参照）。また、ホドスコープ２４（検出器２６，２８）は、通常、センサモジュール３０の非一様性（製造バラツキ）に起因して各部分に検出効率ηのバラツキなどが生じる（図１０参照）。

高エネルギ領域や低エネルギ領域のミュオンのフラックスの不定性は、ミュオンのカットオフエネルギＥｃｕｔの決定精度に依存する。以下、具体的に説明する。まず、高エネルギ領域のミュオンのフラックスの不定性は、カットオフエネルギＥｃｕｔの計算に用いるパラメータｂの確度に依存する。パラメータａ，ｂと障害物の密度長ＤＬとその障害物（密度長ＤＬ）を通過（透過）するのに要する最低エネルギ（カットオフエネルギ）Ｅｃｕｔとの関係を関係を示す次式（３２）（参考文献２参照）の両辺を積分すると、式（３３）が得られる（参考文献３参照）。ここで、式（３３）中、「Ｅ０」は、初期ミュオンエネルギであり、「Ｅ」は残存ミュオンエネルギである。この残存ミュオンエネルギＥが０ＧｅＶのときに、初期ミュオンエネルギＥ０とカットオフエネルギＥｃｕｔとが等しくなる。また、式（３３）中、パラメータｂは、図１１（参考文献１参照）などでは、確定値ではなく、平均値で記載されている。パラメータｂは、ある確率分布に従ってミュオンがエネルギを失う指標となるが、図１１を用いてパラメータａ，ｂを設定する場合、所定エネルギ（Muon critical energyに相当するエネルギ（６９３ＧｅＶ））より高い領域ではカットオフエネルギＥｃｕｔの決定においてパラメータａに比してパラメータｂの影響が大きくなる。パラメータａは、確率的にエネルギを失うのではなく、関係式によってどの程度エネルギを失うのか分かっている値である。所定エネルギ（６９３ＧｅＶ）以下の中エネルギ領域（２９７ＭｅＶ〜６９３ＧｅＶ）では、カットオフエネルギＥｃｕｔの決定にパラメータａが大きく関係し、このパラメータａは、所定エネルギ以下の中エネルギ領域では、それほど大きく変化しないから、カットオフエネルギＥｃｕｔを精度よく求めることができる。一方、所定エネルギ（６９３ＧｅＶ）より高い高エネルギ領域では、カットオフエネルギＥｃｕｔの決定にパラメータｂの影響が大きくなるため、カットオフエネルギＥｃｕｔを精度よく求めようとするとモンテカルロシミュレーションを行なう必要があり、シミュレーションを行なう毎にカットオフエネルギＥｃｕｔの値にバラツキが生じる。したがって、高エネルギのミュオンではカットオフエネルギＥｃｕｔを求める際に不定性が生じるため、ミュオンのエネルギと通過できる障害物の厚さ（密度長ＤＬ）との関係にバラツキが生じる（図９参照）。ミュオンのエネルギスペクトルモデルを一つの関係式で確定的に表わせたとしてもカットオフエネルギＥｃｕｔに不定性が生じるため、ミュオンのフラックスのシミュレーションにおいても不定性が生じる。また、低エネルギ領域の不定性は、ミュオンのエネルギが低い場合には地球磁場の影響によって粒子の進行方向が曲がりやすい、という理由に基づく。図１１では、ミュオンのエネルギが２９７ＭｅＶ未満の領域でミュオンのエネルギが低いほど急激にパラメータａの値が大きくなっている。したがって、低エネルギのミュオンほど通過経路中の物質量の変化に大きく影響されてカットオフエネルギＥｃｕｔの値が変化しやすくなり、大気密度の変化に敏感に反応しやすくなる。一方、２９７ＭｅＶ以上ではパラメータａがそれほど大きく変化しないため、大気密度の変化に反応しにくくなり、カットオフエネルギＥｃｕｔが精度よく決定されやすくなる。以上を踏まえると、低エネルギ領域（例えば２９７ＭｅＶ未満の領域）では地球磁場や大気密度に敏感に反応してカットオフエネルギＥｃｕｔの決定精度が低くなり、高エネルギ領域（例えば６９３ＧｅＶより大きな領域）では確率的にエネルギを失う過程の影響が大きくなりカットオフエネルギの決定精度が低くなり、中エネルギ領域（例えば２９７ＭｅＶ〜６９３ＧｅＶの領域）では低エネルギ領域や高エネルギ領域に比してカットオフエネルギＥｃｕｔの決定精度が高くなる（不定性が小さくなる）と言える。
参考文献２：DONALD E. GROOM, NIKOLAI V. MOKHOV, SERGEI I. STRIGANOV, MUON STOPPING POWER AND RANGE TABLES 10 MeV・100 TeV, Atomic Data and Nuclear Data Tables, Vol. 76, No. 2, July 2001, p. 3
参考文献３：K. Nakamura et al. (Particle Data Group), J. Phys. G 37, 075021 (2010), p. 9,［On line］、［平成２４年１月２３日検索］、インターネット
＜http://pdg.lbl.gov/2011/reviews/rpp2011-rev-cosmic-rays.pdf＞

また、ホドスコープ２４の各部分の検出効率η［φ，θ］のバラツキは、検出器２６，２８におけるミュオンの通過位置α１［ｘ１，ｙ１］，α２［ｘ２，ｙ２］に応じた各方位角φ，各仰角θの検出効率η［φ，θ］が検出器２６の各座標［ｘ１，ｙ１］の検出効率η１[ｘ１，ｙ１]と検出器２８の各座標［ｘ２，ｙ２］の検出効率η２［ｘ２，ｙ２］とに応じて定まる、ためである。

これらを考慮すると、ホドスコープ２４を回転させずに、ホドスコープ２４によりミュオンの通過を検出してカウント装置４０によりカウントした各方位角φ，各仰角θの検出データＯｂｓ［φ，θ］を用いて測定対象の内部の状態を解析する従来例の場合、例えば、測定対象を通過したミュオンのフラックスと測定対象を通過せずに飛来するミュオンのフラックスとを比較して測定対象の内部の状態を解析する従来例の場合、上述の種々の影響を十分に抑制できずにこれらの影響が測定対象の内部の状態の解析結果にある程度反映されていたと考えられる。

一方、実施例では、複数の検出比Ｒｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）を用いて、全体方向範囲Ｒｄｉｔｏの各仰角θにおける測定対象の方位角方向の平均密度ρｍ［θ］を計算したり、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの測定対象の上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］を計算したりすることにより、上述の影響（ホドスコープ２４の回転前か回転後かに拘わらず生じる影響）を抑制して測定対象の内部の状態を解析することができる。

しかも、実施例では、各仰角θについて、各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）の方位角方向の総和としての検出比和Ｓｏｂｓ［θ］と、各方位角φ，φ＋Δφ間の第１解析比Ｒｅｓｔ１［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）の方位角方向についての総和としての解析比和Ｓｅｓｔ［θ］と、を用いて方位角方向の測定対象の平均密度ρｍ［θ］を計算したり、各方位角φ，φ＋Δφ間の検出比Ｒｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）と、各方位角φ，φ＋Δφ間の第２解析比Ｒｅｓｔ２［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ）と、を用いて各方位角φの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］を計算したりするから、ユーザの経験則に基づくパラメータの設定を必要とせず、測定対象の内部の状態の解析をより容易に行なう（解析のより自動化を図る）ことができる。

以上説明した実施例の内部状態解析システム２０によれば、第１状態（ホドスコープ２４の方位角方向の回転前）の第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ）に応じた検出フラックスＦｏｂｓ１［φ，θ］と、第２状態（ホドスコープ２４の方位角方向の回転後）の第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ２ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）に応じた検出フラックスＦｏｂｓ２［φ，θ］と、の比としての複数の検出比Ｒｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ）を用いて、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］（φ：φ１ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）を計算するから、ミュオンの性質（地球磁場による方位角依存性や、中間エネルギ領域より低エネルギ側の領域や高エネルギ側の領域のフラックスのバラツキなど）やホドスコープ２４の各部分のバラツキなどによる影響を抑制して上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］を計算することができ、測定対象の内部の状態をより適正に解析することができる。

実施例の内部状態解析システム２０では、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］（φ：φ１ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）を計算するものとしたが、全体方向範囲Ｒｄｉｔｏの各仰角θにおける測定対象の方位角方向（全体方位角範囲Ｒａｚｔｏ）の平均密度ρｍ［θ］を計算するものとしてもよい。この場合、図４の内部状態解析ルーチンのステップＳ１６０〜Ｓ２１０の処理を実行しないものとすればよい。この場合でも、ミュオンの性質（地球磁場による方位角依存性や、中間エネルギ領域より低エネルギ側の領域や高エネルギ側の領域のフラックスのバラツキなど）やホドスコープ２４の各部分のバラツキなどによる影響を抑制して各仰角θにおける測定対象の方位角方向の平均密度ρｍ［θ］を計算することができ、測定対象の内部の状態をより適正に解析することができる。

実施例の内部状態解析システム２０では、上下限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］，ρｆｍｉｎ［φ＊，Δφ，θ＊］の計算において、密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］をある程度大きな正の初期値Δρ０（例えば、０．５０ｇ／ｃｃなど）から徐々に小さくしながら上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］，下限密度比ρｆｍｉｎ［φ＊，Δφ，θ＊］の順に計算するものとしたが、密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］をある程度小さな負の初期値Δρ１（例えば、−０．５０ｇ／ｃｃなど）から徐々に大きくしながら下限密度比ρｆｍｉｎ［φ＊，Δφ，θ＊］，上限密度比ρｆｍａｘ［φ＊，Δφ，θ＊］の順に計算するものとしてもよい。

実施例の内部状態解析システム２０では、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］（φ：φ１ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）を計算するものとしたが、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの密度ρ［φ，θ］（φ：φ１ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）を計算するものとしてもよい。この場合の内部状態解析ルーチンの一例を図１２に示す。図１２の内部状態解析ルーチンは、図４の内部状態解析ルーチンのステップＳ１８０，Ｓ２１０の処理に代えて、ステップＳ１８０ｂ，Ｓ２１０ｂの処理を実行する点を除いて、図４の内部状態解析ルーチンと同一である。したがって、同一の処理については同一のステップ番号を付し、その詳細な説明は省略する。図１２の内部状態解析ルーチンでは、対象方位角φ＊に下限方位角φ１ｍｉｎを設定すると（ステップＳ１７０）、対象仰角φ＊，対象仰角θ＊における測定対象の方位角方向の密度勾配Δρ［φ＊，θ＊］を図１３に例示する密度勾配計算処理により計算する（ステップＳ１８０ｂ）。以下、図１２の内部状態解析ルーチンの説明を一旦中断し、図１３の密度勾配計算処理について説明する。

図１３の密度勾配計算処理では、まず、対象仰角θ＊における方位角φ＊，φ＊＋Δφの方位角方向の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］に初期値Δρ２（例えば、０．００ｇ／ｃｃなど）を設定する（ステップＳ８００）。

続いて、図７の上下限密度比計算処理のステップＳ５１０〜Ｓ５４０の処理と同様に、方位角φ＊，φ＊＋Δφの密度ρ［φ＊，θ＊］，ρ［φ＊＋Δφ，θ＊］を計算し（ステップＳ８１０）、方位角φ＊，φ＊＋Δφの測定対象の第２密度長ＤＬ２［φ＊，θ＊］，ＤＬ２［φ＊＋Δφ，θ＊］を計算し（ステップＳ８２０）、方位角φ＊，φ＊＋Δφの第２密度長起因フラックスＮ２［φ＊，θ＊］，Ｎ２［φ＊＋Δφ，θ＊］を計算し（ステップＳ８３０）、方位角φ＊，φ＊＋Δφ間の解析比Ｒｅｓｔ２［φ＊，Δφ，θ＊］を計算する（ステップＳ８４０）。

そして、方位角φ＊，φ＊＋Δφ間の検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］と解析比Ｒｅｓｔ２［φ＊，Δφ，θ］との差分が検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］についての統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］未満であるか否かを判定する（ステップＳ８６０）。この処理は、現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］が適正であるか否かを判定する処理である。

検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］と解析比Ｒｅｓｔ２［φ＊，Δφ，θ］との差分が統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］以上であると判定されたときには、現在の密度勾配Δρ［φ＊、θ＊］が適正でないと判断し、現在の密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］を変更して（ステップＳ８７０）、ステップＳ８１０に戻る。ここで、密度勾配Δρ［φ＊，θ＊］，Δρ［φ＊＋Δφ，θ＊］の変更は、例えば、０．００ｇ／ｃｃ，０．０１ｇ／ｃｃ，−０．０１ｇ・ｃｃ，０．０２ｇ／ｃｃ，−０．０２ｇ／ｃｃ・・・の順に変化するように変更することができる。

こうしてステップＳ８１０〜Ｓ８７０の処理を繰り返し実行して、ステップＳ８６０で検出比Ｒｏｂｓ［φ＊，Δφ，θ＊］と解析比Ｒｅｓｔ２［φ＊，Δφ，θ］との差分が統計誤差ＥＲｏｂｓ［φ＊，Δφ，θ＊］未満であると判定されたときには、現在の密度勾配Δρ［φ＊、θ＊］は適正であると判断し、それを対象方位角φ＊，対象仰角θ＊における密度勾配Δρ［φ＊，θ＊］として確定して（ステップＳ８８０）、本ルーチンを終了する。

以上、図１３の密度勾配計算処理について説明した。図１２の内部状態解析ルーチンの説明に戻る。ステップＳ１８０ｂで対象方位角φ＊，対象仰角θ＊における測定対象の方位角方向の密度勾配Δρ［φ＊、θ＊］を計算すると、対象仰角φ＊が上限方位角φ１ｍａｘであるか否かを判定し（ステップＳ１９０）、対象仰角φ＊が上限方位角φ１ｍａｘでないと判定されたときには、現在の対象仰角φ＊に角度ｄφを加えたものを新たな対象仰角φ＊に設定して（ステップＳ２００）、ステップＳ１８０ｂに戻る。

こうしてステップＳ１８０ｂ〜Ｓ２００の処理を繰り返し実行して、ステップＳ１９０で対象方位角φ＊が上限方位角φ１ｍａｘであると判定されたときには、ステップＳ１６０で設定した中心方位角φｃの密度ρ［φｃ，θ＊］（＝ρｍ［θ＊］）と次式（３４），（３５）を用いて、対象仰角θ＊における各方位角φ（φ１ｍｉｎ〜φ２ｍａｘ）の密度ρ［φ，θ］を計算する（ステップＳ２１０）。具体的には、式（３４）中、「ρ［φ＋Δφ，θ＊］」，「Δρ［φ，θ＊］」にそれぞれ中心方位角φｃの密度ρ［φｃ，θ＊］，密度勾配Δρ［φｃ，θ＊］を代入して方位角φｃ−Δφの密度ρ［φｃ−Δφ，θ＊］を計算すると共に同様に方位角φｃ−２Δφ〜下限方位角φ１ｍｉｎの密度ρ［φｃ−２Δφ，θ＊］〜ρ［φ１ｍｉｎ，θ＊］を計算し、式（３５）中、「ρ［φ，θ＊］」，「Δρ［φ，θ＊］」にそれぞれ中心方位角φｃの密度ρ［φｃ，θ＊］，密度勾配Δρ［φｃ，θ＊］を代入して方位角φｃ＋Δφの密度ρ［φｃ＋Δφ，θ＊］を計算すると共に同様に方位角φｃ＋２Δφ〜上限方位角φ２ｍａｘの密度ρ［φｃ＋２Δφ，θ＊］〜ρｍａｘ［φ２ｍａｘ，θ＊］を計算する。こうした処理により、ミュオンの性質（地球磁場による方位角依存性や、中間エネルギ領域より低エネルギ側の領域や高エネルギ側の領域のフラックスのバラツキなど）やホドスコープ２４の各部分のバラツキなどによる影響を抑制して密度ρ［φ，θ］を計算することができ、測定対象の内部の状態をより適正に解析することができる。

次に、対象仰角θ＊が上限仰角θｍａｘであるか否かを判定し（ステップＳ２２０）、対象仰角θ＊が上限仰角θｍａｘでないと判定されたときには、現在の対象仰角θ＊に仰角方向の角度分解能としての角度ｄθを加えたものを新たな対象方位角φ＊に設定して（ステップＳ２３０）、ステップＳ１５０に戻る。こうしてステップＳ１５０〜Ｓ２３０の処理を繰り返し実行して、ステップ２２０で対象仰角θ＊が上限仰角θｍａｘであると判定されたときには、計算結果をディスプレイ７０に出力したりハードディスクドライブに記憶させたりして（ステップＳ２４０）、本ルーチンを終了する。計算結果の出力では、全体方向範囲Ｒｄｉｔｏの各仰角θにおける測定対象の方位角方向の平均密度ρ［θ］や、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの測定対象の密度ρ［φ，θ］などを出力する。

この変形例によれば、第１状態（ホドスコープ２４の回転前）の第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ）に応じた検出フラックスＦｏｂｓ１［φ，θ］と、第２状態（ホドスコープ２４の回転後）の第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ２ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）に応じた検出フラックスＦｏｂｓ２［φ＋Δφ，θ］と、の比としての複数の検出比Ｒｏｂｓ［φ，Δφ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ）を用いて、全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの密度ρ［φ，θ］（φ：φ１ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）を計算するから、ミュオンの性質（地球磁場による方位角依存性や、中間エネルギ領域より低エネルギ側の領域や高エネルギ側の領域のフラックスのバラツキなど）やホドスコープ２４の各部分のバラツキなどによる影響を抑制して密度ρ［φ，θ］を計算することができ、測定対象の内部の状態をより適正に解析することができる。

実施例の内部状態解析システム２０では、第１状態の第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ）と、ホドスコープ２４を第１状態から方位角方向に所定角度Δφだけ回転させた第２状態の第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ２ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）と、を測定対象の内部の状態の解析に用いるものとしたが、これに代えて、第１状態の第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］（φ：φｍｉｎ〜φｍａｘ，θ：θ１ｍｉｎ〜θ１ｍａｘ）と、ホドスコープ２４を第１状態から仰角方向に所定角度Δθ（例えば、数度〜十数度程度）だけ回転させた第２状態の第２方向範囲Ｒｄｉ２の各方位角φ，各仰角θ＋Δθの検出データＯｂｓ２［φ，θ＋Δθ］（φ：φｍｉｎ〜φｍａｘ，θ＋Δθ：θ２ｍｉｎ（＝θ１ｍｉｎ＋Δθ）〜θ２ｍａｘ（＝θ１ｍａｘ＋Δθ））と、を測定対象の内部の状態の解析に用いるものとしてもよい。なお、この場合、第１状態か第２状態かに拘わらず方位角範囲（φｍｉｎ〜φｍａｘ）は同一となる。

この場合、検出データＯｂｓ１［φ，θ］を規格化して検出フラックスＦｏｂｓ１［φ，θ］を計算すると共に検出データＯｂｓ２［φ，θ＋Δθ］を規格化して検出フラックスＦｏｂｓ２［φ，θ＋Δθ］を計算し、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出フラックスＦｏｂｓ１［φ，θ］をそれぞれ仰角方向に所定角度Δθだけ異なる方向の検出フラックスＦｏｂｓ２［φ，θ＋Δθ］で除して複数の検出比Ｒｏｂｓ［φ，θ，Δθ］を計算し、計算した複数の検出比Ｒｏｂｓ［φ，θ，Δθ］を用いて測定対象の内部の状態（各方位角φにおける測定対象の仰角方向の平均密度ρ［φ］や、各方位角φ，各仰角θの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］，各方位角φ，各仰角θの密度ρ［φ，θ］など）を解析するものとしてもよい。

また、この場合、複数の検出比Ｒｏｂｓ［φ，θ，Δθ］を計算するのに加えて、第１方向範囲Ｒｄｉ１と第２方向範囲Ｒｄｉ２とからなる全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの測定対象の通過経路長Ｌ［φ，θ］（φ：φｍｉｎ〜φｍａｘ，θ：θ１ｍｉｎ〜θ２ｍａｘ）に基づいて実施例と同様の手法により全体方向範囲Ｒｄｉｔｏの各方位角φ，各仰角θの密度長起因フラックスＮ［φ，θ］を計算し、第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの密度長起因フラックスＮ［φ，θ］をそれぞれ所定角度Δθだけ異なる方向の密度長起因フラックスＮ［φ，θ＋Δθ］で除して複数の解析比Ｒｅｓｔ［φ，θ，Δθ］を計算し、計算した複数の検出比Ｒｏｂｓ［φ，θ，Δθ］と複数の解析比Ｒｅｓｔ［φ，θ，Δθ］とを用いて測定対象の内部の状態（各方位角φにおける測定対象の仰角方向の平均密度ρ［φ］や、各方位角φ，各仰角θの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］，各方位角φ，各仰角θの密度ρ［φ，θ］など）を解析するものとしてもよい。

実施例の内部状態解析システム２０では、第１状態の第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの検出データＯｂｓ１［φ，θ］（φ：φ１ｍｉｎ〜φ１ｍａｘ，θ：θｍｉｎ〜θｍａｘ）と、ホドスコープ２４を第１状態から方位角方向に所定角度Δφだけ回転させた第２状態の第２方向範囲Ｒｄｉ２の各方位角φ＋Δφ，各仰角θの検出データＯｂｓ２［φ＋Δφ，θ］（φ＋Δφ：φ２ｍｉｎ〜φ２ｍａｘ，θ：θｍｉｎ〜θｍａｘ）と、を測定対象の内部の状態の解析に用いるものとしたが、これに代えて、第１状態の方向範囲Ｒｄｉの各方位角φ，各仰角θの検出データＯｂｓ１［Ｐ，φ，θ］（Ｐ：第１状態でのホドスコープ２４の位置，φ：φｍｉｎ〜φｍａｘ，θ：θｍｉｎ〜θｍａｘ）と、ホドスコープ２４を第１状態から水平方向や鉛直方向に所定距離ΔＰ（例えば、３０ｍや５０ｍ，７０ｍなど）だけ平行移動させた第２状態の方向範囲Ｒｄｉの各方位角φ，各仰角θの検出データＯｂｓ２［Ｐ＋ΔＰ，φ，θ］（Ｐ＋ΔＰ：第２状態でのホドスコープ２４の位置，φ：φｍｉｎ〜φｍａｘ，θ：θｍｉｎ〜θｍａｘ）と、を測定対象の内部の状態の解析に用いるものとしてもよい。ここで、平行移動とは、第１状態か第２状態かに拘わらず方位角範囲Ｒａｚ（φｍｉｎ〜φｍａｘ）や仰角範囲Ｒｅｌ（θｍｉｎ〜θｍａｘ）が同一となる移動をいう。図１４は、検出装置２２（ホドスコープ２４）を方位角方向に平行移動させる場合における検出装置２２（ホドスコープ２４）や測定対象を天頂から見た様子を模式的に示す模式図である。図中、実線のホドスコープ２４や第１状態のホドスコープ２４を示し、点線のホドスコープ２４は第２状態のホドスコープ２４を示す。また、平行移動の距離として３０ｍや５０ｍ，７０ｍなどを用いるのは、以下の理由による。例えば、山の径が１ｋｍ程度で、方位角方向の角度分解能が１００ｍｒａｄ（約６°）の検出装置２２（ホドスコープ２４）を山の中心から５００ｍ程度離れた位置に配置する場合を考える。この場合、第１方位角範囲Ｒａｚ１（φ１ｍｉｎ〜φ１ｍａｘ）の中心方位角と中心方位角から１００ｍｒａｄだけ異なる方向の方位角とは、５００ｍ先で５０ｍ程度ズレることになる。したがって、このズレに相当する距離だけホドスコー２４を平行移動させることにより、実施例と同様に、測定対象の内部の状態を解析できると考えられる。

この場合、検出データＯｂｓ１［Ｐ，φ，θ］を規格化して検出フラックスＦｏｂｓ１［φ，θ］を計算すると共に検出データＯｂｓ２［Ｐ＋ΔＰ，φ，θ］を規格化して検出フラックスＦｏｂｓ２［Ｐ＋ΔＰ，φ，θ］を計算し、各方位角φ，各仰角θの検出フラックスＦｏｂｓ１［Ｐ，φ，θ］をそれぞれ同一方向の検出フラックスＦｏｂｓ２［Ｐ＋ΔＰ，φ，θ］で除して複数の検出比Ｒｏｂｓ［Ｐ，ΔＰ，φ，θ］を計算し、計算した複数の検出比Ｒｏｂｓ［φ，θ］を用いて測定対象の内部の状態（各方位角φ，各仰角θの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］，各方位角φ，各仰角θの密度ρ［φ，θ］など）を解析するものとしてもよい。

また、この場合、複数の検出比Ｒｏｂｓ［φ，θ］を計算するのに加えて、第１状態でのホドスコープ２４の位置Ｐにおける第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの測定対象の通過経路長Ｌ［Ｐ，φ，θ］（φ：φｍｉｎ〜φｍａｘ，θ：θｍｉｎ〜θｍａｘ）に基づいて実施例と同様の手法により第１状態でのホドスコープ２４の位置Ｐにおける第１方向範囲Ｒｄｉ１の各方位角φ，各仰角θの密度長起因フラックスＮ［Ｐ，φ，θ］を計算し、第２状態でのホドスコープ２４の位置Ｐ＋ΔＰにおける第２方向範囲Ｒｄｉ２の各方位角φ，各仰角θの測定対象の通過経路長Ｌ［Ｐ＋ΔＰ，φ，θ］（φ：φｍｉｎ〜φｍａｘ，θ：θｍｉｎ〜θｍａｘ）に基づいて実施例と同様の方法により第２状態でのホドスコープ２４の位置Ｐ＋ΔＰにおける第２方向範囲Ｒｄｉ２（＝Ｒｄｉ１）の各方位角φ，各仰角θの密度長起因フラックスＮ［Ｐ＋ΔＰ，φ，θ］を計算し、各方位角φ，各仰角θの密度長起因フラックスＮ［Ｐ，φ，θ］を同一方向の密度長起因フラックスＮ［Ｐ＋ΔＰ，φ，θ］で除して複数の解析比Ｒｅｓｔ［Ｐ，ΔＰ，φ，θ］を計算し、計算した複数の検出比Ｒｏｂｓ［φ，θ］と複数の解析比Ｒｅｓｔ［Ｐ，ΔＰ，φ，θ］とを用いて測定対象の内部の状態（各方位角φ，各仰角θの上下限密度ρｍａｘ［φ，θ］，ρｍｉｎ［φ，θ］や各方位角φ，各仰角θの密度ρ［φ，θ］など）を解析するものとしてもよい。

実施例の内部状態解析システム２０では、ホドスコープ２４を用いるものとしたが、ネオンガスが封入された複数のガラス管を電極（金属板）間に積み重ねたいわゆるホドスコープ・チェンバーを用いるものとしてもよい。このホドスコープチェンバーでは、ミュオンの通過の検出に応じて高電圧パルスが電極間に印加されると、複数のガラス管のうちミュオンが通過したガラス管内に放電が生じる。したがって、この放電が生じたガラス管の並びによってミュオンの進行方向（方位角φ，仰角θ）を検出することができる。

実施例の内部状態解析システム２０では、検出装置２２は、回転架台２３とホドスコープ２４とを備えるものとしたが、回転架台２３に代えて回転できない基盤を用いるものとしてもよい。

実施例の内部状態解析システム２０では、ホドスコープ２４は、検出器２６，２８を回転架台２３の上面に固定して構成するものとしたが、検出器２６，２８間の距離を調節できるように回転架台２３に取り付けて構成するものとしてもよい。

実施例の内部状態解析システム２０では、ホドスコープ２４は、２つの検出器２６，２８を備えるものとしたが、検出器の数は２つに限られず、３つ以上としてもよい。

実施例の内部状態解析システム２０では、ホドスコープ２３の検出器２６，２８のセンサモジュール群２６ａ，２６ｂ，２８ａ，２８ｂのそれぞれのセンサモジュール３０の数は、１０個に限定されるものではなく、８個や１２個，１４個などしてもよい。また、センサモジュール群２６ａ，２６ｂ，２８ａ，２８ｂのそれぞれのセンサモジュール３０の数は、同一に限定されるものではなく、互いに異なるものとしてもよい。

実施例の内部状態解析システム２０では、高透過性粒子としてのミュオン（ミュー粒子）の性質を用いて測定対象（例えば、山や産業プラント，ダムなど）の内部の状態を解析するものとして説明したが、ミュオンに限られず、所定の高透過性（例えば、１メートルの鉄を透過する透過性）を有する粒子（例えば、中性子やＸ線など）の性質を用いて測定対象の内部の状態を解析するものであればよい。

実施例では、内部状態解析システム２０として説明したが、内部状態解析方法の形態としてもよいし、内部状態解析方法を１以上のコンピュータに実現させるためのプログラム（内部状態解析プログラム６０）の形態としてもよいし、このプログラムを記憶する記憶媒体の形態としてもよい。

実施例の主要な要素と課題を解決するための手段の欄に記載した発明の主要な要素との対応関係について説明する。実施例では、回転架台２３とホドスコープ２４とを備える検出装置２２が「検出装置」に相当し、カウント装置４０と解析装置５０とが「解析手段」に相当する。

なお、実施例の主要な要素と課題を解決するための手段の欄に記載した発明の主要な要素との対応関係は、実施例が課題を解決するための手段の欄に記載した発明を実施するための形態を具体的に説明するための一例であることから、課題を解決するための手段の欄に記載した発明の要素を限定するものではない。即ち、課題を解決するための手段の欄に記載した発明についての解釈はその欄の記載に基づいて行なわれるべきものであり、実施例は課題を解決するための手段の欄に記載した発明の具体的な一例に過ぎないものである。

以上、本発明を実施するための形態について実施例を用いて説明したが、本発明はこうした実施例に何等限定されるものではなく、本発明の要旨を逸脱しない範囲内において、種々なる形態で実施し得ることは勿論である。

本発明は、内部状態解析システムの製造産業や、内部状態解析方法のステップをコンピュータに実現させるためのプログラムの製造産業などに利用可能である。

２０ 内部状態解析システム、２２ 検出装置、２３ 回転架台、２４ ホドスコープ、２６，２８ 検出器、２６ａ，２６ｂ，２８ａ，２８ｂ センサモジュール群、３０ センサモジュール、３１ シンチレータ、３２ 光電子増倍管、４０ リードアウトモジュール、４２ ディスクリミネータ、４４ コインシデンス検出部、４６ フラックス演算部、５０ 解析装置、５２ コンピュータ、６０ 内部状態解析プログラム、６２ 入力モジュール、６４ 解析モジュール、６６ 出力モジュール、７０ ディスプレイ、７２ キーボード、７４ マウス。

Claims (10)

1. 複数のシンチレータを有し且つ該シンチレータの長手方向とは異なる方向で該シンチレータを通過するミュオンを検出可能な検出器を複数備える検出部を有する検出装置が測定対象の周辺に配置されたときの前記検出部による検出結果を用いて前記測定対象の内部の状態を解析する内部状態解析方法であって、
   前記検出装置が前記測定対象の周辺に配置された第１状態での前記検出部による検出結果に応じた前記ミュオンのフラックスである第１状態検出フラックスと、前記検出部で前記ミュオンの通過を検出可能な該検出部から見た方向の範囲である検出方向範囲が前記第１状態から変更された第２状態での前記検出部による検出結果に応じた前記ミュオンのフラックスである第２状態検出フラックスと、の比を用いて前記測定対象の密度関連物理量を演算する、
   ことを特徴とする内部状態解析方法。
2. 請求項１記載の内部状態解析方法であって、
   前記第２状態は、前記第１状態から前記検出部が方位角方向に所定角度だけ回転された状態であり、
   前記第１状態における前記検出方向範囲である第１方向範囲の各方向の前記第１状態検出フラックスと、前記第２状態における前記検出方向範囲である第２方向範囲の各方向の前記第２状態検出フラックスと、のうち互いに前記所定角度だけ異なる方向の前記第１状態検出フラックスと前記第２状態検出フラックスとの比である複数の検出比を用いて前記測定対象の密度関連物理量を演算する、
   内部状態解析方法。
3. 請求項２記載の内部状態解析方法であって、
   前記複数の検出比と、前記第１方向範囲と前記第２方向範囲とからなる全体方向範囲の各方向の前記測定対象の通過経路長に基づく該全体方向範囲の各方向の前記ミュオンのフラックスである経路長起因フラックスのうち、互いに前記所定角度だけ異なる方向の前記経路長起因フラックスの比である複数の解析比と、を用いて前記測定対象の密度関連物理量を演算する、
   内部状態解析方法。
4. 請求項３記載の内部状態解析方法であって、
   前記全体方向範囲の各仰角について、各方位角の前記測定対象の通過経路長と前記測定対象の方位角方向の平均密度とに基づく各方位角の前記経路長起因フラックスを用いて得られる前記複数の解析比の方位角方向についての総和と、前記複数の検出比の方位角方向についての総和と、の差分が第１許容値以下となるよう前記測定対象の方位角方向の平均密度を演算する、
   内部状態解析方法。
5. 請求項４記載の内部状態解析方法であって、
   前記第１方向範囲の各方位角，各仰角ついて、前記演算した測定対象の方位角方向の平均密度と前記測定対象の方位角方向の密度勾配と前記測定対象の通過経路長とに基づく前記経路長起因フラックスを用いて得られる前記解析比である勾配起因解析比が前記検出比より大きく該検出比に第２許容値を加えた値より小さくなるよう前記所定角度だけ異なる方位角との前記測定対象の密度比の上限としての上限密度比を演算すると共に前記勾配起因解析比が前記検出比より小さく該検出比から前記第２許容値を減じた値より大きくなるよう前記所定角度だけ異なる方位角との前記測定対象の密度比の下限としての下限密度比を演算し、前記全体方向範囲の各方位角，各仰角について、前記演算した上限密度比を用いて該測定対象の密度分布の上限としての上限密度を演算すると共に前記演算した下限密度比を用いて該測定対象の密度分布の下限としての下限密度を演算する、
   内部状態解析方法。
6. 請求項４記載の内部状態解析方法であって、
   前記第１方向範囲の各方位角，各仰角ついて、前記演算した測定対象の方位角方向の平均密度と前記測定対象の方位角方向の密度勾配と前記測定対象の通過経路長とに基づく前記経路長起因フラックスを用いて得られる前記解析比である勾配起因解析比と前記検出比との差分が第３許容値以下となるよう前記測定対象の方位角方向の密度勾配を演算し、前記全体方向範囲の各方位角，各仰角について、前記演算した測定対象の方位角方向の密度勾配を用いて前記測定対象の密度を演算する、
   内部状態解析方法。
7. 請求項１記載の内部状態解析方法であって、
   前記第２状態は、前記第１状態から前記検出部が仰角方向に所定角度だけ回転された状態であり、
   前記第１状態における前記検出方向範囲である第１方向範囲の各方向の前記第１状態検出フラックスと、前記第２状態における前記検出方向範囲である第２方向範囲の各方向の前記第２状態検出フラックスと、のうち互いに前記所定角度だけ異なる方向の前記第１状態検出フラックスと前記第２状態検出フラックスとの比である複数の検出比を用いて前記測定対象の密度関連物理量を演算する、
   内部状態解析方法。
8. 請求項１記載の内部状態解析方法であって、
   前記第２状態は、前記第１状態から前記検出部が水平方向および／または鉛直方向に所定距離だけ平行移動された状態であり、
   前記第１状態における前記検出方向範囲である第１方向範囲の各方向の前記第１状態検出フラックスと、前記第２状態における前記検出方向範囲である第２方向範囲の各方向の前記第２状態検出フラックスと、のうち同一方向の前記第１状態検出フラックスと前記第２状態検出フラックスとの比である複数の検出比を用いて前記測定対象の密度関連物理量を演算する、
   内部状態解析方法。
9. コンピュータを、複数のシンチレータを有し且つ該シンチレータの長手方向とは異なる方向で該シンチレータを通過するミュオンを検出可能な検出器を複数備える検出部を有する検出装置が測定対象の周辺に配置されたときの前記検出部による検出結果を用いて前記測定対象の内部の状態を解析する装置として機能させるプログラムであって、
   前記検出装置が前記測定対象の周辺に配置された第１状態での前記検出部による検出結果に応じた前記ミュオンのフラックスである第１状態検出フラックスと、前記検出部で前記ミュオンの通過を検出可能な該検出部から見た方向の範囲である検出方向範囲が前記第１状態から変更された第２状態での前記検出部による検出結果に応じた前記ミュオンのフラックスである第２状態検出フラックスと、の比を用いて前記測定対象の密度関連物理量を演算するモジュール、
   を備えるプログラム。
10. 複数のシンチレータを有し且つ該シンチレータの長手方向とは異なる方向で該シンチレータを通過するミュオンを検出可能な検出器を複数備える検出部を有する検出装置と、前記検出装置が測定対象の周辺に配置されたときの前記検出部による検出結果を用いて前記測定対象の内部の状態を解析する解析手段と、を備える内部状態解析システムであって、
    前記解析手段は、前記検出装置が前記測定対象の周辺に配置された第１状態での前記検出部による検出結果に応じた前記ミュオンのフラックスである第１状態検出フラックスと、前記検出部で前記ミュオンの通過を検出可能な該検出部から見た方向の範囲である検出方向範囲が前記第１状態から変更された第２状態での前記検出部による検出結果に応じた前記ミュオンのフラックスである第２状態検出フラックスと、の比を用いて前記測定対象の密度関連物理量を演算する手段である、
    ことを特徴とする内部状態解析システム。