第7章「地殻活動シミュレーション手法」研究計画

 

 

１．はじめに

 

　地震発生の予測には地殻活動を監視しつつ，その時点での地殻内部の物理プロセスと応力の現在値を明らかにし，地震発生に至る地殻活動を評価する必要がある．これまでの地震予知研究においては主として観測に主眼がおかれ，地震の前兆現象を捕らえることを最大の目的としてきたが，最近の理論や数値シミュレーションの発展により，充分なデータがあれば地震の発生に至るまでのプロセスに対しある程度のシナリオを描けるところまで研究が進んできたと言えるであろう．もちろん，現在提唱されている物理モデルはまだ仮説の段階であるから，実データに基づいて検証しなくてはならない．実際の地球にこのようなモデルを直接当てはめるのは，地殻の構造の不均質性により，簡単にはいかない．このような不均質性を考慮しつつ，地震発生の物理モデルをもとに地殻活動データをモデルに同化して逐次検証とモデルの改良を進めていけばある程度定量的に地殻活動の予測が可能ではないだろうか．

このような考え方に基づき，新たに策定された建議「地震予知のための新たな観測研究計画の推進について」においては，地殻活動シミュレーション手法の開発を計画の基本方針の一つとしてとりあげるに至った．そこで，これを受けて，大学の研究者はいくつかの関連する研究課題を提案し，平成11年度より研究に着手した．研究開始から日が浅いことから,平成11年度の段階においては東京大学地震研究所からの課題しか提案されなかった．

 

２．平成１１年度の研究成果

 

（課題名：地殻応力・歪変化シミュレーション手法に関する研究（東京大学地震研究所［課題番号：0112］）

 

「地震予知のための新たな観測研究計画」においては，平成11年からはじまる５カ年計画の目標として“中長期・広域の地殻活動を対象とする「広域モデル」と地震発生準備過程の最終段階及び地震発生・波動伝播過程を対象とする短期・特定域の「特定域モデル」”の２種類の課題を提起している．本課題ではこれらのうち「広域モデル」をとりあげるものである．

「新たな計画」ではこの「広域モデル」においては，「主としてGPS全国観測データを用い，地殻変動シミュレーションモデルを構築する．さらに，このモデルにより数年後の予測を行い，その予測誤差の評価を行うとともに，観測，実験へのフィードバックを目指す」ことが目標とされている．

以上の全体計画をふまえ，本研究課題においては，主として地理院のGPS全国観測網に基づく変位速度ベクトルを基礎データとして，地殻構造やプレート運動に関するデータなどを取り込みつつ，地殻活動を監視する手法を開発することを主たる目標とした．具体的には，１）まず，基礎的な作業として，有限要素法(FEM)などによる日本列島３次元模型を作成し，日本列島の変位速度場がどの程度まで周囲のプレート運動の相互作用によって規定されているのか，を明らかにすること，またそれに基づき，２）日本列島の変位速度場から地殻内部の応力場を推定する手法を開発し，地震活動などの資料を取り込みつつ日本列島の応力場を推定する手法を開発すること，を目標に掲げた．なお，５年間のうちに，可能ならば，この手法をさらに「予測」にまで拡張してある程度の地殻応力変化予測とその検証を試みる，ところにまで踏み込む予定である．

平成11年度においては，次の２つの項目について研究を実施した．

１） プレート運動に基づく日本列島の変位場の再構築

２） GPSデータに基づく地殻応力の推定

以下，章を追って各項目毎の成果を述べる．

 

２−１．プレート運動に基づく日本列島の変位場の再構築

 

最近のGPS観測研究により，日本列島の周囲のプレートの相対運動が精密に明らかにされるようになった（例えば，Kato et al., 1998; 小竹他，1998）．小竹他（1998）はフィリピン海プレート内部で実施されたGPS観測に基づいてフィリピン海プレートのユーラシア大陸安定地塊に対するEuler極を精密に推定した．また，Heki et al.(1999)はアムールプレートのEuler極を同様に推定している．

このような極位置に基づいて，日本列島の東西からのプレートの収束過程がかなり精確にわかってきた．一方，日本列島内部の変位場は地理院のGPS全国観測網によって明らかにされている．もちろん，日本列島内部の変位場はプレート運動に起因するものだけでなく，局所的な力源や構造不均質によって複雑な変位場を形成していると考えられるが，大局的には日本列島の変位場は周囲から収束するプレートの相互作用によって規定されていることは明らかである．そこで，前述したプレート運動によって，日本列島の変位場がどこまで再現できるか，３次元有限要素法によって試行的に日本列島の変位場モデルを構築することを試みた．

このため，まず必要なFEM解析プログラムを導入すると共に，いくつかの簡単なモデルからかなり複雑なモデルまで段階を追って3次元有限要素模型を構築した．図1にこれまでに作成した最も詳細な模型の俯瞰図を示す．この模型は要素数・節点数が約8万で，８面体と６面体の複合体から構成されており，海洋プレートとその下（海側）は表現されていない．プレート上面の位置についてはKosuga et al. (1996), Ishida (1992), Iidaka et al. (1991), 山崎・大井田（1985）などを，また，モホ面位置はZhao et al. (1992)を参照している．現在，沈み込むプレートの面上に変位を与えて，列島の変位場をどこまで再現できるか，計算を実施中である．

なお，この数値モデルと比較すべき観測データとしては地理院の年平均変位速度ベクトル場が用いられるが，観測点はメッシュ上に存在しているわけではない．そこで，地理院の変位速度データにEl-Fiky et al.（1997）の最小二乗予測法を適用し，ノイズを取り除くと共に，格子点上の変位を算出して，計算値と直接比較できるようにした．

 

２−２．GPSデータに基づく地殻応力の推定

 

国土地理院のＧＰＳ全国連続観測網によって常時計測される大量かつ高精度の変位データに対し，より高度な解析をすることが望まれている．本研究項目では「地殻活動モニタ」を解析の目標と考えている．「地殻活動モニタ」とは，地殻応力と応力−ひずみ関係（構成則）を推定し，列島各地域でひずみ・応力等が刻々と変化していく様子をモニタすることをさす（図２参照）．

本研究における「地殻活動モニタ」の根幹は，応力推定の新しい逆解析理論である．これはひずみから応力を逆解析する全く新しい理論であり，元来は，複合材料等の材料サンプルの局所的な構成則を推定するために構築された．理論の一般性は高く，長さの尺度は全く異なるものの，衛星データから計測される変位を用いることで，列島の地域毎の応力推定に適用することが可能である．

以上を背景に，本研究では，衛星計測データを利用した応力分布と構成則を推定する数値計算手法を開発した．この手法では列島全域を対象とし，日々計測される地殻変動データを入力し，応力分布等を出力することで，「地殻活動モニタ」を実現することを試みる．

具体的な研究開発題目は，１）応力分布推定のコードの開発，２）構成則推定のコードの開発，の二つである．前者に関しては逆解析理論の精緻化と列島応力推定コードの開発，後者に関しては列島構成則推定コードの開発が研究項目である．

列島の応力分布推定は，実際に広域応力を測定することが不可能であるため，理論の厳密性が推定の正しさの鍵を握る．それと同時に，誤差を含む衛星データを数値計算に用いることにも考慮しなければならない．これを受けて，応力の代わりに応力ポテンシャル（Airy応力関数）を逆解析するよう，理論の精緻化を行った．応力ポテンシャルはその勾配が釣り合い状態にある応力分布を作るものであり，応力ポテンシャルに対する境界値問題を解くことで，応力が逆解析される．この理論をコード化した結果，より高精度・高速な応力分布推定が可能となった．実際，応力を直接推定する場合には応力テンソルの各成分を推定し，しかも変位の３階微分が入力として必要であったものが，応力ポテンシャルの境界値問題を解く場合には解析対象が１つの関数であり，しかも変位の１階微分のみが入力として必要である．

構成則推定では，図４に示すように，計測されたひずみと逆解析された応力の時系列データから列島各地点の構成則を推定するものである．なお，応力逆解析には，構成則パラメータを一つ設定しなければならない．開発された解析コードは，繰り返し計算によってパラメータを決定する．これは，１）時系列データから推定された構成則を使ってパラメータの更新量を計算，２）更新されたパラメータを使って応力と構成則を再推定，という二つの手順からなり，適当な範囲にパラメータが収束するまでこの作業を繰り返す．

開発されたコードは，数値計算を行う他に，データファイルのファイル処理が必要である．これは衛星より毎日相当のデータが送られるためである．このため，数値計算用にはFORTRAN，ファイル処理用にはc/c++を用いた．応力推定と構成則推定の二つをあわせたコードの構成を図５に示す．

 

３．まとめと今後の展開

 

　本年度では，地殻活動予測のためのシミュレーションの基礎的な課題を提起し，プレート運動と日本列島内部の変形，及び，地表変位に基づく地殻内部応力の推定を実施した．いずれもまだ試行段階であり，まだ成果が出ていない部分もある．

応力の推定に関しては，ここで実施したのは地理院のデータに基づく数年の応力変化だけなので，地殻の応力の絶対値がわかるわけではない．また，これから直ちに地殻活動に関して何らかの予測ができるわけではない，など残された課題は数多い．また，問題を２次元にしており，地殻内部を見ているわけではない，など計算上の問題がある．さらに，理論上の仮定が，実際の地殻の振る舞いとして正しいのかどうか，など検証しなくてはならない部分が残されている．今後さらに計算を進めあるいは解析手法を高度化することなどにより，５カ年の内に少なくともある程度実用的な地殻活動モニタプログラムを構築する必要がある．

もう少し大きな枠組みとして，地殻活動予測シミュレーションの研究部会においては，より広範な研究者の結集のもとに，地震発生に至る地殻の諸活動についてのシミュレーション研究を取り入れることが考えられており，一層の研究の展開が期待できよう．

　考慮すべき問題いくつかある．このような数値シミュレーションでは，結果の検証と観測へのフィードバックが必要になる．現在のところ，検証の方法についての模索はあるものの，確実な手法が見いだせているとは思えない．今後，地震サイクルのシミュレーションなども取り入れられる可能性があり，長時間を視野に入れたシミュレーションに対しどのような検証手法があるのか，またその科学的意義は何か，など，十分検討しておく必要があろう．一つのやり方としてアナログ実験的検証手段を取り込む方法も考えておいてよいだろう．また，他の観測計画とより密接に連携してどのような計算結果が観測計画へのフィードバックに有効か，についても検討すべきであろう．シミュレーション分野の研究はまだはじまったばかりなので，他の研究項目における観測や実験と相互に刺激しあいながら予知研究を進展させる必要があろう．

 

 

文献

 

El-Fiky, G. S., T. Kato, and Y. Fujii, Distribution of the vertical crustal movement rates in the Tohoku district, Japan, predicted by least squares collocation, J. Geodesy, 71, 432-442, 1997.

Heki, K., S. Miyazaki, H. Takahashi, M. Kasahara, F. Kimata, S. Miura, N. Vasilenko, A. Ivashcenko, and K.-D. An, The Amurian plate motion and current plate kinematics in eastern Asia, J. Geophys. Res., 104, 29147-29155, 1999.

Iidaka, T., M. Mizoue, I. Nakamura, T. Tsukuda, K. Sakai, M. Kobayashi, T. Haneda, and S. Hashimoto, The upper boundary of the Philippine Sea platebeneath the western Kanto region estimated from S-P and P-S converted waves,Bull. Earthq. Res. Inst., 66, 333-349, 1991.

Ishida, M., Geometry and relative motion of the Philippine Sea plate and Pacific plate beneath the Kanto-Tokai district, Japan, J. Geophys. Res., 97, 489-513, 1992.

Kato, T., Y. Kotake, S. Nakao, J. Beavan, K. Hirahara, M. Okada, M. Hoshiba, O. Kamigaichi, R. B. Feir, P. H. Park, M. D. Gerasimenko, and M. Kasahara, Initial results from WING, the continuous GPS network in the western Pacific area, Geophys. Res. Lett., 25, 369-372, 1998.

Kosuga, M. T. Sato, A. Hasegawa, T. Matsuzawa, S. Suzuki, and Y. Motoya,Spatial distribution of intermediate-depth earthquakes with horizontal or vertical model planes beneath northeastern Japan, PEPI, 93, 63-89, 1996.

小竹美子・加藤照之・宮崎真一・仙石新，GPS観測に基づくフィリピン海プレートの相対運動と西南日本のテクトニクス，地震II，51, 171-180, 1998.

山崎文人・大井田徹，中部地方におけるフィリピン海プレート沈み込みの形状，地震II, 38, 193-201, 1985.

Zhao, D., S. Horiuchi, and A. Hasegawa, Seismic velocity structure of the crust beneath the Japan Islands, Tectonophys., 212, 289-301, 1992.

＜図の説明＞

図１　日本列島の３次元有限要素モデルの俯瞰図．

図２　地殻モニタの例：1999/8/12〜8/14の間に計測されたひずみ(左)と推定された応力(右)．

図３　応力分布推定：統計処理によるＧＰＳデータからひずみを計算し，ついで応力を逆解析する．

図４　構成則推定：列島各地域の計測ひずみと推定応力の時系列データから構成則を推定する．

図５　解析コードの構成：青枠内がコードの作業内容，赤枠内が入力と出力．