空間統計のモデリングをめぐる諸問題
尾形良彦
統計数理研究所 東京都港区南麻布４−６−７
ogata@ism.ac.jp
Modeling Issues in Spatial Statistics
Yosihiko Ogata
１．はじめに
空間統計は点、線、図形、画像などの様々な幾何学的対象の統計解析である。その様なデータから複雑性や生成メカニズムや不確定性を推定する。そのための統計的モデルの研究の歴史は浅いので、これからの発展性については洋々たるものが予感される。ただ、しばらくは、汎用性のあるモデルや一般的理論の大きな研究成果を期待するのは難しく、各分野の個別の問題に依拠したデータ解析や応用を極めるなどして方法論の深化が進むと考えられる。
演者は主に地震活動のメカニズムの解明や確率的予測の向上を目指した研究を進めているが、これまでの研究実践についてなるべく一般的で他分野でも関連しそうな事項を取り集めてみた。とくに本講演では専ら点配置に関係する統計解析を問題にする。ただし、種村正美氏との共同研究である点配置の相互作用の話題やマルコフ連鎖シミュレーション法の話題は割愛する。この問題を含み空間統計学としてここ20年前後から世界的に発展している内容についての邦書 [1] が最近出版されたのでご興味のある方は参照されたい。
空間統計における方法論として以下のような事を取り上げる。まずパラメタの推定値の誤差評価を可能にする最尤法によるモデリングについて述べたい。次に、空間統計のモデル化における簡素化の前提はデータの定常性 (stationary)、等方性 (isotropy) などであり、これでかなりモデルの自由度（パラメタ数）を絞り安定的かつ実用的なモデル推定が考えられてきた。しかし空間的データが豊富に揃うと非定常性(non-stationary)や異方性 (anisotropy)などの不均質性に関する議論が可能となる。このためには最尤法などの枠組を超えたベイズ理論が必要である。これがどの様に適用され、どのような結果が得られるのか、幾つかの具体例をもってお話したいと考える。
２．地震現象における逆ベキ統計法則をめぐる問題 ([2], [3], [4])
地震の大きさは通常マグニチュードで与えられているが、これは大雑把に言って、エネルギー（またはモーメント量）の対数スケールに対応する。大量の地震のマグニチュードに関する統計をとるときれいに指数分布に従う（Gutenberg-Richter則）。これはエネルギーやモーメントに関して逆ベキ分布であることを示している。この逆ベキ指数は地震活動を論ずる上で重要な役割を持っている。
地震破壊の開始点の位置は震源と呼ばれている。大きな地震が起きると多数の余震が近傍に続発する。これは本震の断層面付近の応力分布が急変し比較的小さな破壊が急速に起きたものと考えられ、それらの破壊開始点が余震の震源である。余震の震源分布は一様ランダムでなく集中クラスター型のパタンを示し、余震間の距離の分布としては逆ベキの減衰を示している。その指数からフラクタル次元が決まる。フラクタル次元は本震断層付近の破砕の複雑性を反映していると考えられる。余震発生数の時間的減衰も逆ベキ法則に従っている（改良大森公式）。
これらの時間・空間・規模のそれぞれのベキ指数は一つの破壊の緩和過程からの帰結なので何らかの相関があると期待できる。これを統計的に調べるにあたって逆ベキ指数の推定値の誤差を議論しなければならない。すなわち誤差評価をするための統計的枠組として最尤法が必要とされている。
３．点配置間の写像変換 ([5])
学部で微分幾何学を学んだ人は抽象的な意味での多様体の写像変換の概念をご存知と思うが、写像変換そのものを推定する問題に地震学でお目にかかった。長野県にある気象庁の松代精密地震観測室は半径5ｋｍ程度の円周上と中心付近にある７つの地震計のアレイによって世界中の地震の震央（緯度、経度）を推定している。したがってその位置の推定誤差は小さくない。しかし、驚く事には、その誤差は地震波の読み取り誤差によるものより、地球内部の地震波伝播速度の違いによる系統的な誤差（偏差）の方が圧倒的に大きいのである。従ってその空間的な偏りが推定できれば、松代で決めた位置をこれで逆変換すれば震源の補正が出来る。
ここでは変換写像をスプライン関数（曲面）で展開し、その係数を最小2乗法で推定したいのであるが、パラメタ（係数）の数が膨大であるので現実的な解が得られないという問題がある。そこでスプライン関数（曲面）の滑らかさの制約（ペナルティ）を入れてペナルティ付最小２乗法 [6] を考えるのであるが、制約の強さの加減をどのようにするかという問題が残る。ここでベイズ法の考え方とType II maximum likelihood の方法 [7] を紹介して、上記の震源補正への応用結果
[5] を述べたい。
４．点配置からアナログ画像を推定する
写真は見た目にはアナログ画像であるが、実際は大量の光粒子（photon）が乾板に衝突した跡の集合体であると考えられる。画像の明暗はそこに光粒子が頻繁に衝突したか否かによって出てくるものと考えられる。明るいところでは光粒子が大量にあるのでカメラの瞬間的なシャッターで写真が出来あがるが、暗黒の中の宇宙や暗室中の生体などに含まれた同位元素などの分布画像は長い時間をかけた粒子衝突の蓄積を必要とする。その間、被写体が動かなければ良いのであるが、そうでない場合も結構多い。
ここでは短時間の衝突による点配置から被写体の大まかな像を推定する方法を考えてみる。点配置（点過程）の強度関数という概念が紹介される。また地震の発生配置から地震活動度の空間分布も同じ概念のもとに推定される[8]。
前節のようにスプライン関数で強度関数（曲面すなわち画像）を推定したいところであるが、この問題に最小２乗法は使えない。このばあい代わりに点過程の対数尤度関数が考えられる。前節と同様にスプライン関数（曲面）のパラメタ（係数）の数が膨大であるので現実的な解を得るため滑らかさの制約（ペナルティ）を入れてペナルティ付対数尤度を考える。しかし、今回は非正規モデルなので Type II maximum likelihood 法の適用が近似的にならざるを得ない([8], [9])。
さらに、ここでは滑らかさの制約について様々な可能性（たとえば等方性・異方性）があり、またスプライン関数以外の別の（たとえばデロネ分割に基づく[10]）関数展開を考える可能性もある。このように、様々なオプションがある場合、赤池ベイズ情報量規準（ABIC, [11]）によって、これらのモデル同士の適合度を比較することが出来る。
５．階層的ベイズモデル
１節で紹介した地震のマグニチュード頻度の指数（b値）は場所や時間によって違う。これは地殻地質、断層群の複雑性や応力場変化度などを反映していると考えられている。b値の時間・空間的な違いが誤差を含めてベイズ法で客観的に推定できる([12], [13])。
以上は洩れなく検出されているような大きさの地震データに対して可能な推定である。ところが、例えば気象庁編集による地震の震源カタログ（1926~現在）の内容を見るとすぐ分かるのは、年代が進むに連れて収録されている地震の数が急速に多くなっていることである。これは小さい地震が検出されて地震の数が急速に増えているからである。また地域的にも地震の多いところとそうでないところがある。これも地震活動の高低のためとは限らない。地震計のネットワークが密で良く整備されている地域（たとえば内陸部）とそうでない所（たとえば沖合いの地域）では、地震の検出率が異なるからである。或る統計モデルを考えるとベイズ法に基づき各時期、各地域の地震検出率と真の地震活動（b値）が一挙に推定できる([14], [15])。ここで考えられているモデルは確率的に構造的な欠測のある時間・空間データの解析にも役に立つ[16]。
６．デロネ分割を使ったベイズモデルと残差解析 [10]
時空間ETAS点過程モデル[17]は、統計地震学における余震活動に関する研究成果に基づいて作成された地震活動予測モデルである。このモデルのパラメタ数は7個であるが、そのうちの重要な5つのパラメタを位置の関数であると考えて、夫々 2次元デロネ分割上の複三角面として展開し、推定をする（階層的時空間ETAS（HIST-ETAS）モデル）。推定すべき係数（パラメタ）の数は地震データ数の5倍以上になるので、三角面の傾きに関する制約のペナルティ付き対数尤度を考えて、ベイズ法によって5関数夫々の制約の強さを客観的に決めることになる。これを気象庁の震源データ (1926年以降、マグニチュード5以上)
にあてはめて、常時活動度、余震強度や余震減衰指数（p値）などのパラメタの空間変化を推定し、日本の地震活動の地域性に関する定量的な特徴付ができる。
ついで、このモデルで予測される地震活動度に対する実際の地震活動の相対比を見るため、時空間の3次元デロネ分割を使ったベイズ型「残差モデル」を構成する。気象庁震源データを再びあてはめて、地震活動の相対的な静穏化や活発化などの時空間領域を3次元画像によって示す。これによって大地震の前に相対的空白域が明瞭に見られる場合がある。
デロネ分割に基づいて関数を展開して推定する方法の特長は集中型・不均質のサンプリングに対応して効率よく関数の変化を表現できるところにある。密な所は詳しく、疎なところは大雑把にという具合である。特に空間が高次元になるにつれて局所的に細かい変化の表現はBスプラインなど他の展開法では困難と考えられる。
文献
尾形良彦 トップページ に戻る
English Page に戻る
統数研のページ
に戻る
Updated on