気象庁精密地震観測室技術報告 第28号 51～58頁 平成23年３月
GEONETデータによる松代近傍の歪場推定
小池 哲司
Strain rate estimation near Matsushiro from GEONET data
Tetsuji KOIKE
１．はじめに
精密地震観測室では石英管伸縮計と水管傾斜計による地殻観測を行っている．また，石英管伸縮計による時系列と比較するために，国土地理院の電子基準点 GPS連続データ(GEONET)の「日々の座標値(F3解)《の日値から歪量を計算している[石川・小久保・山本 (2006)]．図１に松代近傍のGEONET観測局の位置を示し，図２にGEONETの長野観測局のデータを示す．観測開始からしばらくは変動が大きい が2000年以降は収まっている．これは2003年頃のアンテナ改修やGEONETのルーチン解析の向上等により解の精度が良くなっていると考えられ る．2002年以降の精度は1cm程度に収まり， 2004年にトレンドの変化が見えるが，これは長野以外のGEONETデータにも見られ，中越地震によるものと考えられる．変化傾向としては，南北方向は 南への変位を示す観測点が多い．東西方向は松代より北側や西側の観測点は東方向への変位を示しているが，佐久や群馬の観測点は東西方向への変化は小さい．
本報告ではGEONETが展開した1996年4月以降のデータを使用した．2003年にアンテナ交換によるギャップが発生しているが，観測点ごとに補正 を行った．松代近傍でGEONET展開当初から運用している観測局は長野(950267),東部(950268),大町(940046)で，現在歪量はこ の３点で計算を行っているが，その後観測点の移設や新規設置が行われており，2002年に松代に近い更埴(020984)と戸隠(020983)が新たに 設置されている．そのため，観測点を固定する３点法とは違う解析方法が望ましい．そこで，多点GPSデータから任意の場所の歪量・変位を求める解析を調査 した．また，欠測処理に有効な状態空間モデルによる推定方法との比較も行う．
２．GPSデータ３点による歪場推定（３点法）
３点の測地データにより歪場の算出方法としてFrank(1966)による以下の方法がよく用いられている(ノーテーションは笠原・杉村(1978)による)．
ある点P(x,y)での水平変位を(u,v)とする．点Pでから(u,v)だけ離れた転P’(x’,y’)の変位を(u’,v’)とすると，両変位の関係は
となる．これを書き換えると，
となる．ここで，場の歪量の直線歪exx，eyyと剪断歪exyおよび回転量ωを
と定義すると
u’＊u=exx δx＋exy δy＋ωδy
v’＊v=eyy δx＋exy δx＊ωδx
と場の基本歪量のみで簡潔に表すことができる．
また，面積膨張Θと最大剪断歪Σも良く使われる量で
と定義されている．exx, eyy, exy, ωと違いスカラー量のため座標軸の向きに上変である．これらから歪の最大値・最小値 γmax, γmin やその方向θmax, θminは
γmax = Θ/2＋Σ , γmin = Θ/2＊Σ,
θmax = tan-1 ((γmax - θxx)/exy), θmin = tan-1((γmin - θxx)/exy)
であり，これらを求めその歪場の方位を求め，それを評価する研究は盛んに行われている．
exx, eyy, exyが与えられると任意の方向 n への直線歪enは n のx,y軸への方向余弦を(cosθ,sinθ)とすると，
en = cos^2 θ ・ exx ＋ sin^2 θ ・ eyy ＋ 2cosθsinθ ・ exy
と表すことができる．これによりGPS等により３方向の直線歪 e1, e2，e3と方向余弦を得れば
の連立方程式を解くことにより，exx, eyy, exyが求めることができる．
e1, e2, e3はGPSにより解析区間の原点での水平基線長 l を求め，ある時刻の水平基線長 l’ から直線歪量
を観測点ごとに算出して行う．ここで，水平基線長の基準値は1996年4月の値をとる．
GEONETのデータは地心座標と緯度・経度・標高データであるので，緯度・経度を精密地震観測室の歪観測装置の 設置位置を原点(北緯36°32′33.63″, 東経138°12′25.25″)とした平面直角座標に変換し水平面での変化量としている． 平面直角座標の変換公式は国土地理院のHP記載(参考文献のアドレスを参照)の式を参照した． ただし，赤道からの子午線長を求める計算式は河瀬(2009)の式を用いた．
水平面解析区間の最初の座標値を基準とし，その変化量からの基線長を求めその量から直線歪を求めれば，exx, eyy, exyの時系列変化を求めることができる．
この方法では３点で作る領域内の歪場を示しており，任意の領域での歪場の算出には各三角網での歪量もしくは，各GPS点のデータから任意の点の変位量を最適内挿の方法等により補間する方法によって行われてた．
GEONETは観測点の新設や廃止の可能性があり，当初は長野北部の長野・東部・大町の３点で解析を行ってきたが，2002年に松代に近い更埴・戸隠の ２点が新設されている．また，松代に最も近い観測点である長野は2005年～2006年までのデータが乱れており，その部分を補完するには，他の三角網の データによる解析結果を用いたいが，三角網が異なれば歪の長期のトレンドも大きく異なる．できるだけ３点より多い観測データから歪基本量を計算できれば， 欠測や観測点移設に対処しやすくなる．
３．多点GPSデータによる任意の位置での歪場推定（多点法）
多点測地データからの歪基本量の算出方法としてShen et al. (1996)による方法がある．ある点P0 を基準点とする．この点での測地は行われておらず，変位量は未知でΔu0, Δｖ0とする．P0からΔx, Δyだけ離れた点Pでは測地が行われており，その観測された変位量をΔu , Δvとすると，前節で述べたexx, eyy, exyとの関係式は
Δu = Δu0 + exx Δx + exy Δy + ωΔy
Δv = Δv0 + eyy Δy + exy Δx - ωΔx
となる．これを行列で示すと
となる．未知数はΔu0, Δv0, exx，eyy，exy, ωの6個なので，3点以上の観測点があれば上記の方程式は解ける．観測点をN個使用する場合は
簡略して書くと
y = H x
と書ける．これを最小二乗法で解く．
ただし，この方法では各観測点での変位が同じように原点での変位・歪基本量に反映されてしまうため，できるだけ基準点から遠い観測点の影響を小さくするために
Wi = exp{-(Δxi^2 + Δyi^2 ) / σ^2 }
の重みをつける．ここでσは観測点の影響基準を示す定数である．観測点がσより遠くなると重みは大きく減衰することになる．重み行列Wを
W = diag (w1,w2,─,wN)
の対角行列に表すと y = H x は
W y = WH x
と書ける．一般にWHは正方行列ではないので両辺に左からtHをかけ正規方程式にする．
tHWHは正方行列なので逆行列を求め左からかけると
と各観測点での変位量 y から歪基本量 x を求めることができる．
これで観測点の位置に影響されない任意の点の変位量と歪基本量を求めることができるが，いくつか問題がある． 第一にσの値に任意性があるということである．今回の解析ではσの大きさによって大きくトレンドが変化した． 前節での三角形の基線長変化から求める方法および，石英管歪計との比較からσ～22.5kmとしたが， σを数kmかえるとトレンドの正負が逆転するほどの大きな変化である． ３点での推定では，トレンドは1×10＾-6[strain]程度に収まるはずであるが，本手法ではそれ以上の大きな変化をしている． もう１つ推定される原点での変位量がごく小さい問題がある．各点での変位から，原点での変位のオーダーが1×10^-2 [m]と 推定されるが，本手法では1×10^-10 [m]とかなり小さい．
また，誤差も考慮した方程式を考えると，
y = H x + e
であり， e は観測値 y の誤差である．Eを観測値の誤差分散行列とすると，この重みつき最小二乗法の解は
となる．この方法で求めた解は誤差を考慮しない場合とは，変位は若干大きくなる傾向があが，解の挙動は大きく変わらなかった．
４．状態空間モデルによる歪場推定
ここでは，非定常時系列データからトレンドを推定する方法として，状態空間モデルによる解析方法を用いる[Harrison and Stevens 1976, Kitagawa and Gersch 1984, Harvey 1989, 北川 2004]．
まず状態空間モデルの構築にあたり，３節で述べた方程式を基本方程式とする．
y ′= H′ x ′+ e ′
これに前節で用いた重み関数Wを左からかけると
W y ′= WH′ x ′+ W e ′
となる．Wは３節で用いたものと同じσ～22.5kmを使用する．
この式を簡単にし，かつ時系列であることから時刻も考慮にいれて
y (t) = H x (t) + e
とする．この方程式を観測方程式といい， y (t)は観測変数で実際に観測される値であるが， x (t)を状態変数と呼ばれる内部変数で，推定を行う対象である．
次に， y (t)が１期前データ x (t-1)との間に
x (t) = F x (t-1) + v
の関係があるとし，これを状態方程式という．
２階確率差分方程式のモデルでは
x (t) -2 x (t-1) + x (t-2) = v , v i～N(0, σi^2)
と表せる．この場合局所的に x (t)は直線的に変化し，x (t)-2 x (t-1) + x (t-2)～0となることを表現したモデルで，時系列変化の長期的な変化傾向を捉えている．このモデルではFは
と書け，ここでの状態変数 x (t)は
を用いる．状態変数がtとt-1の分，多点法より倍になるので，観測方程式のHは，
とt-1に関しては0になるよう構成される．
また，e の誤差分散行列をW，v の誤差分散行列をRとする．誤差が時間変化をするモデルも構築できるが，本報告での観測方程式，状態方程式の誤差分散の値は変化しないとする．
次に状態変数 x (t)の真値からの誤差の共分散行列をP(t)とする．ここで時刻tiの観測値までわかっていたとする，つまり観測値 Y i = (y1，y2,─,yi)が与えられている場合である．そのときの時刻tでの状態変数の推定値を x (t|ti)とし，状態変数の共分散行列をP(t|ti)とする．
時刻t-1まで観測値が与えられ，状態変数は x (t-1)まで，状態変数の誤差共分散行列もP(t-1)まで推定されていたとする．その推定値を x (t-1|t-1)，P(t-1|t-1)とすると，次の時刻tでの1期先予測値 x (t|t-1)とP(t|t-1)は次のようにかける．
もし時刻tでの観測値が存在していれば，次の漸化式で観測値による修正を行いながらの時刻tでの状態変数と誤差共分散の推定を行うことができる．
x (t|t) = x (t|t-1) + K(t)(y(t)-H x (t|t-1))
P(t|t)=(I-K(t)H)P(t|t-1)
ここでK(t)はカルマンゲインという．以上の漸化式が与えられると，以下の方法で状態方程式とその誤差分散行列の推定が行える．
①初期値 x (1|1)，P(1|1)を用意し， x (1|1)は原点なので0にする．ただしP(1|1)，計算の安定性のために予測される値より大きめにしておく必要があり，単位行列を初期値とした．
②次に１期後の x(2|1)，P(2|1)を漸化式により推定する．③観測値からK(2)を求め， x(2|2)，P(2|2)と②の推定値を修正する．④以降 x (2|1)，P(2|1)， x (2|1)，P(2|1),‥と漸化式より逐次時系列を求めていく．
以上で時系列が求まるが，欠測になった場合は，使用できる観測方程式の組数が少ないか，近傍の観測点が使え無い場合は③のカルマンゲインによる修正を行わず，そのまま
x (t|t)～ x (t-1|t-1)，P(t|t)～P(t-1|t-1)
とした．次に全時間Nで状態変数とその誤差分散行列を推定し，各時刻の状態変数と誤差共分散値の推定値 x (t|t)， x (t|t-1)，P(t|t)，P(t|t-1)を保存している場合，次の漸化式により通常のカルマンフィルタよりも精度の良い推定が行うことができる．
この漸化式はカルマンスムージングと呼ばれる．この方法では x (N|N), P(N|N)からはじまり， x (N-1|N)，P(N|N-1)，‥と時刻を１ステップずつ戻っていくのが特徴である．誤差分散行列の値はカルマンスムージングを行った後の方が小さい．よって状態変数の推定値も行った後の方がよくなっている．
５．計算と比較
(1)松代石英管歪計での歪観測値
精密地震観測室の石英管歪計は松代の大坑道内に，南北および東西の100m石英管の一方を固定し，自由端の変位量を差動トランスにより測定している．大坑道内温度はほぼ14℃で一定で測定値における温度変化の影響は少ないが，降雨や潮汐の影響がある．測定精度は1×10^-9 [strain]程度である． 図３に石英管歪計での日平均値を示す．ここでは潮汐や降雨の補正は行っていない．
(2) ３点の観測値から求めた歪量
次に３点のGPS観測値から求める方法で基本歪量を計算した．ここでx軸を南北方向，y軸を東西方向とし，以降は基本歪量のexx，eyy，exy をeNS，eEW，eNEをとして表す．長野・東部・大町および山ノ内・東部・大町のデータによる計算値を，図４と５に示す．
解析区間は1996～2009年までで，松代の伸縮計と比較のために図のスケールは2×10^-6 [strain]に統一してある．長野・東部・大町のデータからの解析値と石英管伸縮計の１次トレンドが，南北・東西両成分の直線歪は倍近く違う．また図３と４の一次トレンドも異なる．また山ノ内・東部・大町では2009年ごろにeNEトレンドの変化が現れているが，これは山ノ内局所的な変化の影響である．
(3) 多点GPSデータによる算出方法での解析
３の多点法による計算を行うが，比較として長野・東部・大町の３点と長野・東部・大町・山ノ内・妙高高原・白馬・豊科・和田・軽井沢・嬬恋の10点の場 合の結果を図６,７に示す．この手法では歪量と共に変位量も同時に計算を行うのに特徴があるが，計算結果で求まった変位量は10^ -12 [m]のオーダー程度と長野(950267)の10^-1～10^-2 [m] の変位変化とかなり異なる値である．
計算に使用するσの値は松代の石英管歪計(図３)の値や３点法(図４)の計算値を参考にσ～22.5kmとした．しかし，長野・東部・大町での計算の場合 はσを変えると１次トレンドがかなり変化する．場合によってはトレンドの正負が逆転する場合があり，１次トレンドの推定としては安定とは言えない．ま た，eNEは負のトレンドとなっているが，eNEが正のトレンドになるようσを変えてみたが，eEWの傾きは極端な負の傾きになった．
10点で計算した場合は３点の場合と傾きの大きな変化は無く短周期の変動も小さくなっている．σの変化に対する１次トレンドの変化は小さくなったが，観測 点を変えても大きく変化する問題は残ったままである．2002年の新設点の更埴・戸隠を入れて計算を行ったが，新設点のデータを入れたところの前後で ギャップが生じてしまった．解析を行うには同じ解析区間のデータ同士でないとトレンドを見るのはこの手法でも難しい． トレンドは異なるが，データのノイズは10点では半分程度に減少しており，短期的トレンドが見えやすくなっている．
(4) 状態空間モデルによる多点GPSデータによる推定方法での解析
状態空間モデルでも(3)と同じ３点と10点の観測点を用いて解析を行った．図８，９にその解析結果を示す．多点法と大きく異なるのは変位の推定値が解 析点近くのGEONETの長野(950267)観測値と同程度であり，３点使用時と10点使用時でもその値はほとんど変わらない．３点使用時ではeNSの１次トレンドは小さくなっているが，eEWの１次トレンドがかなり大きくなっている． 10点使用時はさらに短周期変動が小さくなっているが，eEWの１次トレンドも大きく，短期的なトレンドの変化が見られなくなっている．
６．まとめ
多点法と状態空間モデルにより，松代の歪量の推定を行った．多点法では変位量が実測とかなり異なる結果になった．状態空間モデルによる推定方法では変位量は実測に近い値を求めることができたが，σや観測点によって１次トレンドが上安定性となった． またeEWのトレンドは松代の石英管伸縮計の値と比べてかなり大きくなった．観測点が多い程，極短周期の変動は小さいがトレンドが見えづらくなるため，観測点の選別が重要になってくる．
今報告では状態空間モデル，２階差分のトレンドモデルのみで行ったが，ＡＲモデルや季節調整モデルの使用，誤差を時間によらず一定としたことなどまだまだ改良の余地は残っている．
参考文献
Frank, F. C., Deduction of Earth strains from survey data, 1966, Bull. Seismol. Soc. Am., 56, 35-42.
Harrison, P. J. and Stevens, C. F., Bayesian Forecasting (with discussion), 1976, Journal of the Royal Statistical Society, Ser. B, 334, 1-41
Harvey, A.C., Forecasting Structural Time Series Models and the Kalman Filter, 1989, Cambridge University Press.
石川有三・小久保一哉・ 山本剛靖, 北信地方の地殻変動, 2006, 精密地震観測室技術報告, 26,131-136.
Kitagawa, G. and Gersch, W.,1984, A smoothness priors-state space modeling of time series with trend and seasonally, Journal of the American Statistical Association, 79, No. 386, 378-389
北川源四郎,時系列解析入門, 2005, 岩波書店.
笠原慶一・杉村新, 変動する地球, 1987, 岩波講座地球科学.
河瀬和重,緯度を与えて赤道からの子午線弧長を求める一般的な計算式,2009,国土地理院時報,119,p45-55.
国土地理院, 緯度・経度から平面直角座標x,yおよび子午線収差角を求める計算, http:// vldb.gsi.go.jp/sokuchi/surveycalc/algorithm.
Shen, Z-K., Jackson, D.D., and Ge, B.X., 1996, Crustal deformation
across and beyond the Los Angeles basin from geodetic measurements, J.
Geophys. Res., 101, 27957 -27980.