新潟山形県沖地震の揺れの解析 その１

新潟山形県沖地震の揺れの解析　その１

 

 

　先日の中国四川省の直下型地震で、Ｍ6.0，深さ「10 km」はやばい！なんてコメントしたら、18日に来ちゃいましたね。それもＭ6.8，深さ「10 km」が・・・。それでお見舞いというわけではないのですが、地震波を解析してみました。

 

１．震源情報

 

　気象庁発表によれば、発生時刻は2019年6月18日22時22分、震央は北緯38.6°，東経139.5°、Ｍ6.8、深さ10 kmです。ここでは防災科学研究所の強震観測網K-Netのデータで、震源に最も近い温海の観測ステーションの記録を使います。震央距離は8 kmです。ところで新潟は・・・震央距離約87 kmの計測震度4となっておりました。

 

 

２．地震加速度波形

 

　温海観測ステーションのEW（東西），NS（南北），UD（上下）の地震加速度波形を図-1にあげます。サンプリング周波数は100 Hz（0.01秒）、観測時間は148秒です。

　東西方向の最大加速度は約600 gal、南北方向は600 gal超え、上下方向でも200 galです。1 gal＝1 cm/s²ですから、水平方向には概ね0.6 g/cos45°、上下方向には0.2 gとなります。つまり水平方向に自分の体重の約8割がたの力で引っ張られ、同時に子供一人を肩に載せたくらいの衝撃力を、0.01秒間に食らった事になります。インタビューで誰かが言ってましたが、とても立ってはいられませんね。

 

３．基線補正

 

　加速度記録（cm/s²）があれば、1回積分して速度（cm/s）、2回積分して変位（cm）です。数値積分にはふつう台形公式が使われます。台形公式は振動と周期関数のたぐい対しては精度が良好な事、加速度記録はデジタルデータなので折れ線近似以上の情報は得られないから、というのがその理由です。しかしそのまま積分を実行すると、積分結果は見事に発散するのが普通です。そこで基線補正といわれるデータ補正を行います。ただし補正のやり方には諸説あって、ここでは一番単純な方法を採用します。

 

　まずNS加速度をそのまま積分した結果が、図-2です。

 

 

　積分結果は、速度の最大が375 cm/s，変位の最大が27551 cm。最終速度は時速13.5 km/h，移動距離は275.5 mになります。明らかに嘘ですよね？(^^;)。

　グラフを見直すと、速度はほとんど直線で等加速度運動状態です。対応して変位はほとんど放物線です。地震で等加速度運動はあり得ないでしょう。だって図-1の地震の主要動は概ね16～30秒の範囲で、それ以降は明らかに終息に向かってます。にも関わらず速度と変位がずっと増え続けるというのは、納得できません。宇宙空間じゃあるまいし・・・。

　加速度測定値に一定値の誤差があると考えられます。これを零点シフトと言って、現実の計測では必ず起こります。零点シフトの値を推定します。

　零点シフトは地震波がない時の測定値の状態を調べれば良いはずです。図-1のNS加速度を見ると、計測開始の最初の16秒間ほどは地震波が来てないので、そこの平均値をとり時系列全体から引けば良いと考えられます。しかし地震波の正確な始まりの時刻はどこでしょう？。もちろん厳密にはわかりません。また16秒くらいでやっても問題なさそうにも思えます。しかし数値積分はちょっとしたゴミのために図-2のような結果を招きます。出来る限りの事をするために、ここでは以下のようにします。

　K-Netの加速度計仕様によれば、分解能は0.001 galです。0.001 gal以内の変動は、＋0.001 galとも0 galとも－0.001 galとも記録される可能性があります。よって有意な変動は0.002 galよりも大きなものと考えられます。図-3は0.01秒ごとの加速度変動を計算し、0.002 galよりも大きければフラグ1を立て、以下であればフラグ0を立てて時系列順に並べたものです。

 

 

　有意変動は計測時間中に万遍なく分布していますが、非有意変動がほぼ連続的に密集分布している時間帯が、計測開始からしばらく続きます。その範囲と図-1のNS波形を見比べると、この部分は地震波がまだ来てない時間帯と判断できます。そこでその時間帯の最後の時刻t₀＝15.94秒を波形開始時刻と決めます。0～t₀にあるデータの平均値は、a_jを加速度データとすると台形公式で、

 

　　　（１）

 

で求められます。a₁＝a(0)，a_N0＝a(t₀)です。Ａ₀の値を加速度時系列全体から引き、積分をやり直した結果が図-4です。これを［基線補正-1］とします。

 

 

　地震波が来ていないと考えたt＝t₀までは確かに妥当です。速度，変位ともほとんど0をたもってます。しかし地震波が始まると、また等加速度状態に見えます。という事は波形開始後は、それ以前と別の零点シフトになったのでしょうか？。零点シフトが零点ドリフト（漂流）に移行する事も、じつはごく普通に起こります。

　加速度計は電気回路ですが、入力が小さい時と大きい時で回路の特性に微妙な変化が生じるのが現実です。それでも最終速度は15 cm/s，移動距離は9 m程度と大幅に改善されてます。人の歩行速度は1 m/sほどなので、［基線補正-1］は確かに効果があったと考えられます。ちなみにＡ₀＝－2.435 gal～－0.0025 gのゴミですが、速度と変位には大きな影響があります。この程度の加速度が本当に効いてるとすれば、無感地震程度で日本列島は移動しまくりです。

 

　図-4の速度波形を参考に波形開始後の零点シフトも開始前と同様に考えると、測定後半に速度の傾きとしてそれが明確に現れてるように見えます。図-4からは100秒以降で速度は十分直線的かな？と見えるのですが、NS波形の時系列全体を示すと、図-5です。

 

 

 

　100秒以降にも小さな有意振動が見えるんですよね。そこで非有意変動が時系列の終端でほぼ連続的に分布するようになるまで、有意の基準を上げてやります。

 

　図-6で時系列終端にある非有意データの密集時間帯の最初の時刻t₁を、いちおう地震波形の終わりとみなし、それ以降の平均Ａ₁（計算法はＡ₀と同じ）を、加速度時系列のt₀以降のデータから引きます。Ａ₀ではなくＡ₁を引きます。これを［基線補正-2］とします。

 

 

　再計算結果には、まだ少し等加速度状態が見えます。つまり0≦t₀とt₀≦t≦t₁とt₁≦tでは、本当に零点シフトが違うという事です。本当に零点ドリフトが起きてます。この中で、いちおう確実と思えるのは0≦t₀の零点シフトだけです。それでも終端速度は2 cm/s，移動距離は2 mになりました。Ａ₁＝－2.558 galで、Ａ₀＝－2.435 galとほとんど変わりません。それでもこれだけの違いがあります。

 

　こういう事態になると、もっと強い仮定を持ち込むべきだという意見も出てきます。今までやってきた事は、全てデータ処理に先立つ事前情報（仮定）に基づいています。地震が終息に向かえば速度と変位が増加し続けるわけないとか、計測器は時間的にほぼ一定な零点シフト誤差を必ず持つとかは、みなそうです。根本的な話として、加速度のデジタルデータの測定値を折れ線でつなぐのもそうです。これは、自然は少なくとも連続であろうという仮定です。より強い仮定の代表は、次の二つです。

 

　　1) 地震の終端で速度は0になる。つまり地震が終われば地盤は動かない。

　　2) 地震の終端で変位は0になる。つまり地震が終われば地盤は原位置に戻る。

 

　一見問題なさそうな（というか当たり前の）仮定に思えますが、地震の終端を計測データの終わりにとるしかないという問題があります。計測データの終わりで地震が止んだ、なんていう保証はないんですよ。図-5をみると、なんか150秒以降も地震は続いてるようにも見えますが、1)を適用してみると図-7が得られます。この方法は、時刻t₀以降のデータ平均をt₀以降のデータから引いた［基線補正-2］と同じです

 

　図-7の結果については「嘘っくせぇ～！」です(^^;)。何が気に食わないっていって、まず速度が（わずかですが）等加速度状態で0に収束するところです。自然現象がそんなに綺麗に挙動するわけねぇ～だろ～！。対応して変位は、時系列の終端でちょうど頂点に達する放物線になります。

　「嘘だこれは！」「人工的過ぎる！」。ただし[基線補正-1)]で納得できる結果が得られるケースも、けっこうあるのは事実です。

　さらに仮定2)まで適用すると、もっと人工的な結果になる事が多いです。仮定2)には明確な反証があります。例えば東北地方太平洋沖地震（いわゆる東日本大震災3.11）では、地盤が最大5.85 m動いた（動いたまま）と国土地理院がGPS観測に基づいて、はっきり認めています。

 

 

 

　そろそろデータ処理の結果としてどの辺りをめざすのか、決めるべき時です。まず計測データの終端で地震が止んだとは限りません。従って時系列終端で速度と変位が0になるとは限らない。しかし収束傾向は見えて欲しいし、終端速度と変位はあまり大きな値にもなって欲しくない。こんなところでしょう。というか、その辺りが精一杯では？(^^;)。

 

　いずれにしろ、図-6の等加速度状態を除く算段を考える必要があります。この等加速度状態は、t₀≦t≦t₁とt₁≦tでの零点シフトのわずかな差から、t₀≦tで平均的に起きた等加速度状態だと考えられます。それを言ったら、t₀≦t≦t₁とt₁≦tの零点シフトというものだってそれぞれの区間の平均値でしょう。正確な零点シフトの時間依存性（厳密な零点ドリフト）なんてわかりっこありません。

 

 

　でもいま実用的に問題になってるのは、t₀≦t全体の零点ドリフトの平均値としての零点シフトです。そこで気づくのが、図-4のt₀以降の変位が、ある2次関数に載りそうだという事です。変位の近似2次関数の2次の係数の2倍は、その等加速度状態の加速度です。図-4の変位のt₀以降を取り出し、近似2次関数を計算してみると、以下になります。

 

　　　（２）

　相関係数はR²＝0.9999972。ほぼ確実な近似です。この異常な相関性の良さが、零点ドリフトなどの系統誤差の存在を示唆します。

 

　(2)のtの2乗の項の係数の2倍が、零点ドリフトも考慮したt₀≦t全体における平均的な零点シフトＡ₁'＝－2×0.0562446＝－0.1124892 galだと思われます。この値を先の［基線補正-1］の加速度時系列のt₀以降のデータから引き、［基線補正-2'］とします。

 

 

　予想通り、ほぼ等速運動状態になりました。ちなみにＡ₁＝－2.558 gal、Ａ₀＋Ａ₁'＝－2.548 gal。両者の差は、わずかに0.01 gal＝0.1 mm/s²。たったこれだけのゴミで、図-6と図-9の劇的な違いが生まれます。

 

　さて図-9の等速運動状態の物理的解釈です。これはもちろん、(2)のtの1次の項に対応した等速度です。

もしt₀≦tで今度は直線近似を行えば、再び異常な相関性を示すでしょう。これを系統誤差と考えると速度の誤差なので、加速度波形開始からの加速度誤差の累積です。という事は時系列の後半はほとんど等速度で間違いないと思われるので、誤差の累積は波形開始の最初の部分で起きているはずです。特に疑われれるのは加速度波形の主要動部分です。主要動から後の時間では、等速度誤差に影響する誤差の累積はほとんどないと考えられます。

　時系列の終端t₂から逆向きに変位データをtまでたどり、t～t₂間の回帰直線を計算し変位データとの相関をとり、その変化をtのグラフとして表したのが図-10です。予想通り、異常な相関係数の高さです。最低値は0.9986になりました。終端t₂付近の相関係数の低下部は、変位変動がデータ数の少なさを上回った部分で、あまり意味がありません。十分データが増えると直線性が回復し0.9998～0.9999の間で増減します。データの始まり方の増加は、波形の主要動の影響が大きく単なる偶然でしょう。ここで注目するのは、直線性が回復したとみなせる、図-10のt＝t₃の時点です。さっきの話が妥当であれば、この時点の傾き（速度）v₀はかなり信頼できる等速度誤差の値とみなせます。t₃＝119.58秒，v＝0.780254 cm/s。

　カタツムリよりも遅いかもしれません(^^)。でもt₃＝119.58秒のタイミングは、図-6の結果であるt₁＝120.11秒とほとんど同じです。

 

 

 

 

　図-11は、さっきとは逆にデータの始端からt₀～t間の回帰直線を計算し、その傾きとv₀との差をtのグラフとして表したものです。ここで注目するのは、データ始端側で相関係数が最大になる時と、回帰直線の傾きがv₀に最も近くなるタイミングt₄です。さっきの話が妥当であれば、その時点までの加速度誤差の累積で等速度誤差が生じたと考える事は可能です。図-11では理想的に両者がほぼ一致してますが、相関最大時の方を優先します。t₄＝30.33秒となります。この結果と図-1のNS波形を確証バイアス付きの目で比較してやると、t₀＝15.94秒～t₄＝30.33秒の範囲、概ねt＝16～30秒の範囲が地震の主要動に見えてきたりします(^^;)。

　t₄～t₂の回帰直線の傾きはv＝0.780254 cm/sになります。［基線補正-3］として、

 

　　　（３）

をt₀≦t≦t₄の範囲で［基線補正-2'］の結果から引き、再計算します。図-12に結果を示します。図-12の中段のグラフです。

　終端速度と終端変位は、－0.111 cm/sと2.404 cm。この結果の評価ですが、震央距離にして8 km、深さにして僅か10 kmという至近距離でＭ6.0なんていう大規模地震が起こったんですから、本当に新潟県と山形県の沿岸部は2 cmくらい「ずれて止まった」のかも知れません。しかしデータ終端で地震が止まったという保証はないので、まだ変位中かも知れません。そういう思いで変位グラフの末尾を見ると、少しですが減少傾向（収束傾向）が見てとります。また終端速度と終端変位の値は、決して大きなものではありません。この状況は求めていたものです。自分には、かなり自然な結果に見えます。

　同様な方法で基線補正を行ったEWとUDの結果も、図-12に示します。