地震計の特性の補正
このページでは、データ提供のフォーマットの選択肢としてsacフォーマットが用意されている防災科学技術研究所のF-netデータを例に説明します。F-netデータは、地震計で記録された生データが公開されていて、地震計の特性はユーザーが補正するようになっています。地震計の特性は特性関数を記述するためのpoleとzero（いずれも複素数）という形で情報が与えられています。
F-netデータは、広帯域の地震計（速度計）で記録されたデータで、広い周波数対に渡って地震計の特性はフラットゆえ、特性の補正はしなくていいや、と僕は考えたりしましたが（僕の解析では長い方でもせいぜい周期30秒くらいが対象ですので）、それは浅はかでした。というのは、公開されている生データは地動速度という物理量ではなく、（おそらく）電圧の値そのものです。それならば感度（定数）をかけて地動速度に直せばいいかと思うと、感度は与えれておらず、与えられているのはより完全な情報である地震計の特性（を記述するpoleとzero）なのです。poleとzeroから特性のフラット部分の高さ（＝感度）を求める手間をかけるならば、それよりもいっそまじめに特性を補正した方がいいや、と考えてしまいます。というのは、SACにはtransferという地震計の特性を補正する（特性をdeconvolutionする）コマンドがちゃんと用意されているからです。
というわけで、以下ではSACのtransferコマンドを使って地震計の特性を補正する手順の例を紹介します。ただし、僕はSACを使った経験があまりないので、以下は誤りを含む可能性もあります。また、当然ながら解析の対象やデータによっては以下の手順では問題がある可能性もあります。したがって解析は自己責任で行ってください。
transferコマンドその他を使った地震計の特性の補正の手順の例
★以下の手順例でのファイル名の説明
・補正前のデータ（sacフォーマット）のファイル名：TSA_BLE
・poleとzeroのファイル名：TSA_BLE.zp ← F-netの場合波形データとともに提供されます。実に便利でありがたいです。
・補正後のデータ（sacフォーマット）を保存するファイル名：tsa.ew
★poleとzeroのファイルの中身は次のような感じです。
CONSTANT 2.936013e+10
ZEROS 2
0.000000E+00 0.000000E+00
0.000000E+00 0.000000E+00
POLES 3
-0.03614522 -0.03890150
-0.03614522 0.03890150
-7006.220 0.0
※注：F-netのデータとともに配布されてきたpoleとzeroのファイルは
CONSTANT 2.936013e+10
ZEROS 2
0.000000E+00 0.000000E+00
0.000000E+00 0.000000E+00
POLES 3
-3.614522E-02 -3.890150E-02
-3.614522E-02 3.890150E-02
-7.006220E+03 0.000000E+00
のように、POLESの数値が指数表示になっていました。しかし、これでは僕の環境（CPU: Intel Xeon, OS: Vine Linux 2.5CR, SAC 10.6f）ではtransferコマンド実行時にエラーが出て補正できませんでした。これを上のように浮動小数点表示にするとうまくいきました。また、CONSTANTのところを指数表示でなく浮動小数点表示にするとwarningが出て（桁が大きいゆえでしょうか？）、補正はするものの振幅と極性のおかしい結果になりました。このあたりの追求はしていませんので（この現象が普遍的なものか否かも含め）、詳細は不明です。情報をお持ちの方はお教えください。
（03/09/19追記）この辺りなかなか微妙なようで、ZEROSの値も小数表示にしないといけない場合に遭遇しました。これを機に（？）、awk を使って zero と pole の値の指数表示を小数表示に変える C shell script（zp_fmt.csh）を作ってみました。小数点以下の桁数も微妙なようで、なかなか難しいです。（03/09/19追記ここまで）
（ここより実際の手順：赤字がsacのコマンドです）
SAC> qdp off
SAC> bd x
補正前のデータを読み込みます。
SAC> r TSA_BLE
波形をウィンドウに表示します（注：これ以降はウインドウ表示の部分は省略）。
SAC> p
（拡大図）
補正を施す前にゼロ線からのずれ（オフセット）を取り除いておきます。波形データの平均値をさっ引くrmeanというコマンドを使いましょう。
SAC> rmean
（拡大図）
※平均値をさっ引いたのでは、うまくオフセットをゼロにできない場合も当然あります。その場合は、別の方法を考えねばなりません。泥臭い方法ですが、 一例を末尾に記しておきます。
transferコマンドで地震計の特性を補正します。
SAC> transfer from polezero s TSA_BLE.zp
※from：特性を取り除く（＝deconvolutionを行う）という意味。
polezero：地震計の特性がpoleとzeroを記したファイルで与えられることを示す。
s：このあとにpoleとzeroを記したファイル名を記す。
（拡大図）
これで補正が終わりました。補正を施すと多少のゼロ線からのずれや低周波（長周期）の信号が乗ったりします。まずrmeanでオフセットを取り除いておきましょう。
SAC> rmean
（拡大図）
次にこの例では、100秒以上の長周期（0.01 Hzより低周波）は解析の対象外だという状況にあるとして、ハイパスフィルターでそれより低周波（長周期）の信号をカットしてしまいます。
SAC> hp co 0.01
（拡大図）
※注：このような低周波（長周期）の信号をcutする処理を行うかどうかは、当然ながら状況によります。低周波（長周期）側の信号が信頼できるデータかどうかは地震計の種類（F-net観測点には、高感度の広帯域速度計と低感度の広帯域速度計＝強震計の2種類があって、前者の方がより長周期までカバーできる）や環境、それに特性を取り除く前のデータの前処理等によります。更に、関心のある周波数帯域は解析の対象によって当然ながら異なります。この辺りは御自分で判断してください。ここに示す手順は単なる一例だということにご注意ください。
F-netデータの場合、速度計なので振幅の次元は速度で振幅の単位はm/sです。僕は単位をcm/sにしたいので、mulコマンドで振幅の値を100倍することにします。
SAC> mul 100
（拡大図）
これですべての補正作業は終わりました。メモリーにあるデータをtsa.ewという名前のファイルに保存します。
SAC> w tsa.ew
●オフセットを取り除く別の方法（泥臭い！）
SAC> qdp off
SAC> bd x
P波が到着する前の部分（ここでは0秒から15秒の部分）を切り出します。この部分の振幅値がオフセットの値となるわけです。cutコマンドで切り出す範囲（0秒から15秒）を指定しておいてから波形データを読み込ませます。
SAC> cut 0 15
SAC> r TSA_BLE
SAC> p
（拡大図）
これでみると、この場合オフセットの値は -1666.5くらいであることが目視で（笑）わかります。これをデータからさっ引いてやることにします。
まず、0秒から15秒の切り出しをoffにしておいて、あらためて波形全体を読み込んでおき、
SAC> cut off
SAC> r TSA_BLE
subコマンドで -1666.5をさっ引いてやります。
SAC> sub -1665.5
これでオフセットを取り除くことができました。
＜謝辞＞
transferコマンドの使い方は防災科学技術研究所の久保篤規さんに教えていただきました。このページに示した例では、防災科学技術研究所のF-netの波形データを使わせていただいております。記して感謝致します。
「SAC10.6fの使い方」へ戻る
最終更新日：03/09/19