[S15P-17] 速度計波形から震度をリアルタイム演算するフィルターの設計について
気象庁の計測震度は，3成分の加速度波形をフーリエ変換し，周波数領域で周波数に応じた重みを掛け，逆フーリエ変換で時間領域に戻した波形から推定される．この位相遅れのない波形は因果律を満たさないため，厳密にはリアルタイムで計測震度を推定することはできない．一方，リアルタイムで震度の相当値を用いるニーズはあり，例えば，緊急地震速報のPLUM法で用いられている．リアルタイムでの推定には，いくつかの提案があるもののの，功刀・他（2008, 2013）の方法が，本来の定義に極めて近い演算を行っている．そこでは，周波数に応じて掛ける重みとほぼ同じ周波数依存性を，因果律を満たす（ために，位相遅れが生じる）フィルタで表現することにより，本来の計測震度と極めて近い値を推定することが出来る．本来の計測震度との違いは，位相遅れのない波形と，（リアルタイム処理が可能なように）位相遅れがある波形という相違であり，もはや，この位相遅れが異なるのみと言える．
さて，現在，PLUM法では震度相当値をリアルタイムで推定して用いており，加速度計の出力波形を用いている．一方，リアルタイムで伝送されている波形データの数は，加速度計よりも速度計の出力波形の方が多い．緊急地震速報では，観測点密度を上げることが迅速化につながる．このため，速度計の出力波形から，震度相当値をリアルタイム演算することを考える．なお，このリアルタイム演算を用いることにより，速度計と加速度計が併設している観測点では，速度計から推定した波形と，加速度計からの波形とが（振切れ等が無ければ）一致するはずのものであり，観測点の品質管理をリアルタイムで行うことも期待できる．
速度計の出力波形から，震度相当値をリアルタイム演算するには，速度計の計器周波数特性の補正， ((s2+2h0ω0s+ω02)/s2．ここで，s=iωで, ω0とh0は速度計の固有角周波数と減衰定数）と，速度から加速度の変換（つまり，微分），s，の演算を行うフィルタとを合わせた (s2+2h0ω0 s +ω02)/s を考える（高周波側ではsに漸近する）．これを，功刀・他（2013）が提案するフィルタ（以下，K2012フィルタ，高周波でs-6に漸近する）に追加して施せば良さそうである．高周波側ではsとs-6とで都合s-5で漸近すると期待され発散は避けられそうである．しかし，実際に試みると，高周波で発散するようになる．その原因を考察する．
ｓ-ｚ変換として，功刀・他（2013）では，プリワーピング（ωc→2/ΔT·tan(ΔTωc/2)）を行わずに，s-1=ΔT/2·(1+z-1)/(1-z-1), s-2=ΔT2/12·(1+10z-1+z-2)/(1-z-1)2を用いている（ここで，ΔTはデジタイズのサンプリング間隔，z=exp(iωΔT), ωcはハイカットやローカットの遮断角周波数）.この操作で得たｚを用いるフィルタを，以下K2012Aフィルタと呼ぶ）．一方，プリワーピングを施し，s-1=ΔT/2·(1+z-1)/(1-z-1), s-2=ΔT2/4·(1+z-1) 2/(1-z-1)2を用いる方法も一般的に知られている（以下，K2012Bフィルタ）．K2012AとK2012Bフィルタの違いは，K2012A では，ナイキスト周波数(1/(2ΔT))に近づくにつれて（ω→π/ΔT），s-1→0, s-2→-ΔT2/6なのに対して，K2012B では s-1→0, s-2→0であり，高周波でのs-2に現れる．K2012フィルタでは，3つの2次のハイカットフィルタ，ωc2/(s2+2hcωc s +ωc2)，を用いている（遮断周波数ωc/2πは，12,20,30Hz）.これに，K2012Aフィルタを施すと，低周波側ω→0で1である一方，高周波側でナイキスト周波数(1/(2ΔT))に近づくにつれて，-ωc2ΔT2/(6-ωc2ΔT2)となり，0には近づかない．ΔTが0.01秒で遮断周波数が30Hzの場合には，この値の絶対値は約1.45であり，1より大きい（図）．つまり，高周波になるほど振幅が大きくなり，高周波が遮断されずハイカットは機能していない．一方，K2012Bフィルタを施すと，低周波側ω→0で1，高周波側でナイキスト周波数(1/(2ΔT))に近づくにつれて0となり，高周波が遮断される．
さて，2次のハイカットフィルター，ωc2/(s2+2hcωcs +ωc2)，と速度計の計器特性補正と微分特性を組合わせたもの (s2+2h0ω0s+ω02)/s を，合わせると，ωc2(s2+2h0ω0s+ω02)/(s2+2hcωcs +ωc2) /sとなり，ω≫ω0かつω≫ωcの高周波では，ωc2/sで近似され，高周波になるほど（sが大きくなるほど），0に近づくことが予想される．しかし，K2012Aフィルタで2次のハイカットフィルターを通すと，ハイカットは機能しておらず，これに，速度計の計器特性補正と微分特性を合わせたもの（(s2+2h0ω0 s +ω02)/s，高周波でsに漸近する）を施すと，発散する特性になる．一方，K2012Bフィルタでの2次のハイカットフィルターでは高周波でs-2に漸近し，さらに，ナイキスト周波数(1/(2ΔT))に近づくにつれて0に近づくので，(s2+2h0ω0+ω02)/sを施しても発散しない．
以上見てきたように，速度計の出力波形から震度相当値をリアルタイムで求めるフィルタとして，K2012フィルタに速度計の計器特性補正と微分特性を合わせたものを施した際に発散する原因は，ｓ-ｚ変換の際にs-2=ΔT2/12·(1+10z-1+z-2)/(1-z-1)2を用いたことに起因していると言える．この発散は，s-2=ΔT2/4·(1+z-1) 2/(1-z-1)2を用いることで回避することが出来る．よって，速度計の出力波形から，震度相当値をリアルタイム演算する際には，K2012Bフィルタ（プリワーピングを施し，s-1=ΔT/2·(1+z-1)/(1-z-1), s-2=ΔT2/4·(1+z-1) 2/(1-z-1)2）を用いるのが良いと思われる．
さて，現在，PLUM法では震度相当値をリアルタイムで推定して用いており，加速度計の出力波形を用いている．一方，リアルタイムで伝送されている波形データの数は，加速度計よりも速度計の出力波形の方が多い．緊急地震速報では，観測点密度を上げることが迅速化につながる．このため，速度計の出力波形から，震度相当値をリアルタイム演算することを考える．なお，このリアルタイム演算を用いることにより，速度計と加速度計が併設している観測点では，速度計から推定した波形と，加速度計からの波形とが（振切れ等が無ければ）一致するはずのものであり，観測点の品質管理をリアルタイムで行うことも期待できる．
速度計の出力波形から，震度相当値をリアルタイム演算するには，速度計の計器周波数特性の補正， ((s2+2h0ω0s+ω02)/s2．ここで，s=iωで, ω0とh0は速度計の固有角周波数と減衰定数）と，速度から加速度の変換（つまり，微分），s，の演算を行うフィルタとを合わせた (s2+2h0ω0 s +ω02)/s を考える（高周波側ではsに漸近する）．これを，功刀・他（2013）が提案するフィルタ（以下，K2012フィルタ，高周波でs-6に漸近する）に追加して施せば良さそうである．高周波側ではsとs-6とで都合s-5で漸近すると期待され発散は避けられそうである．しかし，実際に試みると，高周波で発散するようになる．その原因を考察する．
ｓ-ｚ変換として，功刀・他（2013）では，プリワーピング（ωc→2/ΔT·tan(ΔTωc/2)）を行わずに，s-1=ΔT/2·(1+z-1)/(1-z-1), s-2=ΔT2/12·(1+10z-1+z-2)/(1-z-1)2を用いている（ここで，ΔTはデジタイズのサンプリング間隔，z=exp(iωΔT), ωcはハイカットやローカットの遮断角周波数）.この操作で得たｚを用いるフィルタを，以下K2012Aフィルタと呼ぶ）．一方，プリワーピングを施し，s-1=ΔT/2·(1+z-1)/(1-z-1), s-2=ΔT2/4·(1+z-1) 2/(1-z-1)2を用いる方法も一般的に知られている（以下，K2012Bフィルタ）．K2012AとK2012Bフィルタの違いは，K2012A では，ナイキスト周波数(1/(2ΔT))に近づくにつれて（ω→π/ΔT），s-1→0, s-2→-ΔT2/6なのに対して，K2012B では s-1→0, s-2→0であり，高周波でのs-2に現れる．K2012フィルタでは，3つの2次のハイカットフィルタ，ωc2/(s2+2hcωc s +ωc2)，を用いている（遮断周波数ωc/2πは，12,20,30Hz）.これに，K2012Aフィルタを施すと，低周波側ω→0で1である一方，高周波側でナイキスト周波数(1/(2ΔT))に近づくにつれて，-ωc2ΔT2/(6-ωc2ΔT2)となり，0には近づかない．ΔTが0.01秒で遮断周波数が30Hzの場合には，この値の絶対値は約1.45であり，1より大きい（図）．つまり，高周波になるほど振幅が大きくなり，高周波が遮断されずハイカットは機能していない．一方，K2012Bフィルタを施すと，低周波側ω→0で1，高周波側でナイキスト周波数(1/(2ΔT))に近づくにつれて0となり，高周波が遮断される．
さて，2次のハイカットフィルター，ωc2/(s2+2hcωcs +ωc2)，と速度計の計器特性補正と微分特性を組合わせたもの (s2+2h0ω0s+ω02)/s を，合わせると，ωc2(s2+2h0ω0s+ω02)/(s2+2hcωcs +ωc2) /sとなり，ω≫ω0かつω≫ωcの高周波では，ωc2/sで近似され，高周波になるほど（sが大きくなるほど），0に近づくことが予想される．しかし，K2012Aフィルタで2次のハイカットフィルターを通すと，ハイカットは機能しておらず，これに，速度計の計器特性補正と微分特性を合わせたもの（(s2+2h0ω0 s +ω02)/s，高周波でsに漸近する）を施すと，発散する特性になる．一方，K2012Bフィルタでの2次のハイカットフィルターでは高周波でs-2に漸近し，さらに，ナイキスト周波数(1/(2ΔT))に近づくにつれて0に近づくので，(s2+2h0ω0+ω02)/sを施しても発散しない．
以上見てきたように，速度計の出力波形から震度相当値をリアルタイムで求めるフィルタとして，K2012フィルタに速度計の計器特性補正と微分特性を合わせたものを施した際に発散する原因は，ｓ-ｚ変換の際にs-2=ΔT2/12·(1+10z-1+z-2)/(1-z-1)2を用いたことに起因していると言える．この発散は，s-2=ΔT2/4·(1+z-1) 2/(1-z-1)2を用いることで回避することが出来る．よって，速度計の出力波形から，震度相当値をリアルタイム演算する際には，K2012Bフィルタ（プリワーピングを施し，s-1=ΔT/2·(1+z-1)/(1-z-1), s-2=ΔT2/4·(1+z-1) 2/(1-z-1)2）を用いるのが良いと思われる．