Title

京都大学生存圏研究所(RISH)が保有する観測データの総合解析演習

Author(s)新堀, 淳樹Citation(2013)

Issue Date2013-02-27

URL

http://hdl.handle.net/2433/171279

Right/ This is not the published version. Please cite only thepublished version. この論文は出版社版でありません。引
用の際には出版社版をご確認ご利用ください。

TypePresentation

Textversion

author

Kyoto University

IUGONETデータ解析講習会@RISH
2013年2月27日(水)京都大学生存圏研究所(RISH)が保有する観測データの総合解析演習
担当：新堀淳樹(京大RISH)

1

0. 講習内容(1) 地球大気の鉛直構造とその観測測器(2) 京大生存研の保有する観測データについて(3) 熱圏・中間圏データの総合解析演習
基礎編：解析ソフト(UDAS)を用いたプロットの作成

応用編：スペクトル解析(コヒーレンス、S変換)
相互相関解析
他のデータセットとの統合解析

2

1. はじめに＜本講習会の目指すところ＞
1．IUGONETプロダクト(MDB、解析ソフトウェア)の実体験
太陽地球科学分野に属する研究者・学生に、この講習会を通じて実際に触れてもらうことによって、その便利さや凄さを実感してもらう。
2．利用者の獲得
IUGONETプロダクトの弱点を指摘してもらい、今後の開発過程の参考要素として取り組む。
3．成果物の増産
多くの研究者がIUGONETのプロダクトを利用することで、太陽地球系の多様な研究が促進されれば、我々開発者の成果だけでなく、各機関のデータ提供もとの成果にもなる
なお、この講習会で取り上げた題材で皆さんを共著に入れた論文が仕上がれば上出来である
3

2. 地球大気の鉛直構造とその観測測器

42. 地球大気の鉛直構造とその観測測器＜RISHのデータがカバーする領域(UDASで解析可能のもの)
＞
電離圏

熱圏

中間圏

赤道大気レーダー(FAI観測)流星レーダーMFレーダー

MUレーダー(電離圏)
アイオノゾンデ

MUレーダー(流星)

MUレーダー
(中間圏)

成層圏

対流圏

赤道大気レーダー(対流圏・成層圏)
境界層レーダーダーウィンゾンデ

赤道域(-10-10度)

MUレーダー(対流圏・成層圏)
下部対流圏レーダー
ウインドプロファイラ
境界層レーダー
自動気象観測装置

中緯度域(35度)

53. 京大生存研の保有する流星・MFレーダー＜インドネシアのMF、流星レーダー群＞
Meteor radar
Kototabang (2002-)

MF radarPontianak (1995-）

Regional network inIndonesia (1992- )流星レーダー(Jakarta, Koto Tabang, Biak)MFレーダー(Pontianak, Pameungpeuk)

Meteor radar
Jakarta (1992-1999)

MF radarPameungpeuk (2004-)

6

3. 京大生存研の保有する流星・MFレーダー＜ビアク、コトタバン流星レーダーを使用する意図＞
・ 2観測点ともにほぼ赤道上にあり、経度差は約36度
・ 観測システムは同一• 地球大気は東西方向に対してほ・ 高度70-110kmの風速を観測
ぼ一様であり、経度方向に離れた2点間で似た現象が観測されると期待される• しかし局地的な変動が加わることによる差異も生じうるので、とても興味深いデータ対である

流星レーダー
Koto Tabang(2002/11-)

南緯：0.203
東経：100.318★

★

流星レーダー
Biak (2011/5-)
南緯：1.175
東経：136.102

73. 京大生存研の保有する流星・MFレーダー＜流星レーダーについて＞

流星レーダーの諸元
ビアク：136.102度、-1.175度
コトタバン：100.320度、-0.204度
中心周波数: 33.3 MHz
37.7 MHz
送信出力: 12 kW
観測高度: 70km-110km (中間圏・下部熱圏風速)
5本のアンテナの干渉計で全天の流星飛跡の方向を計測し、その動きから中性風を導出

レーダーの原理図
8

4. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1]
ビアクとコトタバン流星レーダー観測から得られた熱帯域の熱圏下部の風速の高度時間プロットの作成
日時：2011/11/01-2012/02/29
観測高度範囲：70-110 km
時間分解能：1時間
1段目：ビアクにおける東西風2段目：ビアクにおける南北風
3段目：コトタバンにおける東西風4段目：コトタバンにおける南北風
下部熱圏(80-100 km)の南北成分において2日波が観測されている

9

4. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1](1) ビアクとコトタバン流星レーダーデータ解析-1●日時を指定して高度-時間プロットを作成する-1> timespan, '2011-11-1',152, /day(2011/10/01から152日分のデータの日時指定)> iug_load_meteor_rish, site = 'ktb', parameter = 'h2t60min00', length ='1_month'
ビアク：bik、コトタバン：ktb

bik:2011/5/31、ktb:2002/11/15以降パラメタは、4つ選択可で、‘h2t60min00’以外にも‘h2t60min00’、‘h2t60min30’、と‘h4t60min30’が選択可能
length：1日ファイル→ ‘1_day‘、 1ヶ月ファイル→'1_month'10

4. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1](1) ビアクとコトタバン流星レーダーデータ解析-1●日時を指定して高度-時間プロットを作成する-2> tplot_names （tplot変数名の確認）
風速2成分とその分散、流星痕のtplot変数がロードされていることを確認
なお、キーワード’/ verbose’をくわえて、 tplot_names,1,/verbose とす
れば、選択したtplot変数の中身が表示される。
> tplot, ‘tplot変数名’⇒指定したtplot変数のデータのプロットが出力される。
※複数のプロットをしたい場合は、以下のようにする。
> tplot, ['tplot1', 'tplot2',…]11

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1](2) ビアクとコトタバン流星レーダーデータ解析-2●簡単な解析の演習> zlim, ‘tplot変数名’, 最小値、最大値(カラーバー範囲変更)

例えば、東西風速のコンター範囲を-100から100m/sまでの範囲に変えたい場合、 tplot変数名のところに東西風速を表すtplot変数名を、最小、最大のところに、-100, 100をそれぞれいれればよい。
(ex. > zlim, 'iug_meteor_ktb_uwnd_h2t60min00‘, -100, 100)

> ylim, ‘tplot変数名’, 最小値、最大値(y軸の範囲変更)
例えば、東西風速の高度80－100kmの熱圏下部の部分に拡大したい場
合は、 tplot変数名のところに東西風速を表すtplot変数名を、最小、最大のところに、80, 100をそれぞれいれればよい。(他も同様)(ex. > ylim, 'iug_meteor_ktb_uwnd_h2t60min00‘, 80, 100)12

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1](2) ビアクとコトタバン流星レーダーデータ解析-2●簡単な解析の演習> tsmooth_in_time, ‘tplot変数名’, 時間(秒)例えば、東西風速の6時間の移動平均値を求めたい場合、tplot変数名のところに東西風速を表すtplot変数名を、時間のところに21600秒をいれる。
このまま実行すると、入力したtplot変数名の後に’_smoothed’という文字が付け加わるが、引数newnameを用いれば、お好みのtplot変数名がついたtplot変数を作成できる。
>tsmooth_in_time,'iug_meteor_ktb_uwnd_h2t60min00',21600,newname ='iug_meteor_ktb_uwnd_h2t60min00_ts21600‘

135. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1](3) ビアクとコトタバン流星レーダーデータ解析-3●簡単な解析の演習> tlimit, ‘tplot変数名’, 開始日時、終了日時(時刻範囲変更)(ex. > tlimit, ‘2011-12-01’, ‘2012-03-01’)
例えば、得られているプロットの図に対して、横軸の時刻範囲を2011年12月1日から2012年3月1日に拡大したい場合、このコマンドを用いる。
自分の好きな箇所を拡大したい場合は、tlimitだけを打つと、プロットにカーソルを持ってくると、十字の線が現れるので、各々のところで、右クリックをすれば、そこの時刻の拡大図が出来上がる。
zlim等でカラーバーの範囲をあらかじめ設定しておくと、時刻を拡大した時に、その範囲が自動で設定されなくなる。

14

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題
ビアクとコトタバンの観測点において、
1](3) コトタバン流星レーダーデータ解析-3●簡単な解析の演習下部熱圏(80-100 km)の南北成分に2日波が2012年1月の期間に出現している

東西風(ビアク)

南北風(ビアク)

東西風(コトタバン)
南北風(コトタバン)

15

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1](4) 解析やプロットに必要なtplot変数群
iug_meteor_bik_uwnd_h2t60min00
iug_meteor_bik_vwnd_h2t60min00

iug_meteor_ktb_uwnd_h2t60min00
iug_meteor_ktb_vwnd_h2t60min00※上から順に、ビアクとコトタバンで観測された東西、南北風が格納されたtplot変数群である

16

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [基礎課題1](5) 解析後のプロット等の保存●作成したプロットをpsやpngファイル等へ保存[psファイルへの保存]> popen, ‘ 保存するファイル名’ex. popen, ‘test’> tplot> pclose[pngファイルへの保存]> makepng, ‘保存するファイル名’ex. makepng, ‘test’●数値データの保存[tplot形式で保存]> tplot_save, ‘tplot変数名’,filename=‘**’
ex. tplot_save, ‘**’, filename= ‘**’[テキストファイルへの保存]
 tplot_ascii, ‘tplot変数名’ ex. tplot_ascii, ‘***’17

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](1) 風速のスペクトル解析●クロススペクトル、コヒーレンス解析
ビアクとコトタバンにおける中間圏・下部熱圏風の卓越周波数成分を調べるとともに、両者のスペクトルの相関を見るためにコヒーレンス解析を行う●S(Stockwell)変換
上記の解析で見られた卓越周波数成分の振幅変動の様子を捉えるとともに、時間とともに周波数が変化していないかを調べる
S変換は、ウェーブレット変換を一般化したものであり、時間的に非定常な現象のスペクトルを見るのに向いている18

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](2) クロススペクトル、コヒーレンス解析●データ補間、スペクトル計算の前準備> split_vec, ‘tplot変数名’
例えば、東西風速の各高度毎に分けたtplot変数を作成したい場合、tplot変数名のところに東西風速を表すtplot変数名をいれる。
このまま実行すると、入力したtplot変数名の後に’_**’という番号が付け加わったtplot変数名が生成される。
>split_vec,'iug_meteor_ktb_uwnd_h2t60min00‘※このパラメタの場合、高度70-110kmまでの範囲に2km毎にデータが並んでいるので、tplot変数名の最後の数字からどの高度かの見当がつく。
高度 = 70 + 番号*2

例えば、90 kmの場合は、’ iug_meteor_ktb_uwnd_h2t60min00_10’のtplot変数名のデータとなる19

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](2) クロススペクトル、コヒーレンス解析●データ補間> udata_interpolation, ‘tplot変数名1’, , ‘tplot変数名2’
ビアクとコトタバン上空の東西風速の各高度毎に分けたtplot変数に格納されたデータの補間を行う場合、tplot変数名のところに各観測点の東西風速を表すtplot変数名をいれる。
このまま実行すると、入力したtplot変数名の後に’_interpol’ が付け加わったtplot変数名が生成される。

>udata_interpolation, 'iug_meteor_bik_uwnd_h2t60min00‘,'iug_meteor_ktb_uwnd_h2t60min00‘
これを実行すれば、以下のtplot変数が生成される。
iug_meteor_bik_uwnd_h2t60min00_interpol
iug_meteor_ktb_uwnd_h2t60min00_interpol⇒南北風についても同じ作業を行う (vwndを含むtplot変数)20

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](2) クロススペクトル、コヒーレンス解析●データ補間Start Time : 2011-11-01/00:00:00
End Time : 2012-03-27/23:00:00
3600.0000 [sec] : Time Interval
3648 : number of iug_meteor_bik_uwnd_h2t60min00_10 data
3552 : number of iug_meteor_ktb_uwnd_h2t60min00_10 data

3552 : number of iug_meteor_bik_uwnd_h2t60min00_10_interpol data (interpolatediug_meteor_bik_uwnd_h2t60min00_10 data)
3552 : number of iug_meteor_ktb_uwnd_h2t60min00_10_interpol data (interpolatediug_meteor_ktb_uwnd_h2t60min00_10 data)
340 : interpolation number of iug_meteor_bik_uwnd_h2t60min00_10_interpol data
74 : interpolation number of iug_meteor_ktb_uwnd_h2t60min00_10_interpol data

上記のように、データ補間のコマンドを実行すると、コンソール上にデータ補間点数の情報が出力される→この情報を基にユーザーは、データの利用可否を判断できる

21

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](2) クロススペクトル、コヒーレンス解析●データ補間したものをプロット
東西風(ビアク)

東西風(コトタバン)
南北風(ビアク)

南北風(コトタバン)
22

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](2) クロススペクトル、コヒーレンス解析●コヒーレンス解析> uspec_coh, ‘tplot変数名1’, ‘tplot変数名2’, wd=1
データ補間で生成されたtplot変数に格納されたデータについてクロスス
ペクトル、コヒーレンス解析を行う場合、tplot変数名のところに各観測点の東西風速を表すtplot変数名をいれる。
このまま実行すると、ビアクとコトタバンの東西風速のパワースペクト
ル、コヒーレンス、及び位相のプロット図が現れる。

>uspec_coh, 'iug_meteor_bik_uwnd_h2t60min00_10_interpol‘,'iug_meteor_ktb_uwnd_h2t60min00_10_interpol‘,wd=1⇒南北風についても同じ作業を行う (vwndを含むtplot変数)※この解析を行う前に、必ずデータ補間の作業を行っておく。
wdはプロットwindowの番号であり、南北風では、wd=2とすると、別のwindow 2にプロットが描かれる

23

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](2) クロススペクトル、コヒーレンス解析●コヒーレンス解析

東西風90 km

ビアクの周波数スペク
トル

コトタバンの周波数スペク
トル

コヒーレン
ス

位相24

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](2) クロススペクトル、コヒーレンス解析●コヒーレンス解析

南北風90 km

ビアクの周波数スペク
トル

コトタバンの周波数スペク
トル

コヒーレン
ス

位相25

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析● S変換解析の利点

原油1バレルあたりの単価変動

非定常な変動現象を捉えるために地震分野(P-wave, S-wave解析)で開発された解析手法[Stockwell, 1996]
時間的に周波数と振幅が急激に変化する現象を捉えることが可能⇒地磁気分野：Pi 2脈動、磁気インパルス現象⇒ウェーブレット解析の欠点(ノイズに弱い、振幅の周波数依存)を克服した拡張版
26

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●実演習> ustrans_pwrspc, ‘tplot変数名’, /abs, /sampling
データ補間で生成されたtplot変数に格納されたデータについてS変換解析を行う場合、tplot変数名に各観測点の東西風速を表すtplot変数名をいれる。
引数/absと/samplingは必須である。ここで、 /absを/powerにすればS変換のパワースペクトルが得られる。

>ustrans_pwrspc,'iug_meteor_bik_uwnd_h2t60min00_10_interpol‘,/abs,/sampling
このまま実行すると、入力したtplot変数名の後に’_stpwrspcが付け加わったtplot変数名が生成される。
通常のプロットコマンド(tplot, ‘tplot変数’)を入れればS変換した時系列プロットが得られる⇒y, z軸の調整やy軸のキャプション名を変更する

27

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●実演習> ylim, ‘tplot変数名’, 0, 500
y軸の調整を行う。
そのままでは、高周波のスペクトルがつぶれて見えない
S変換で生成されるtplot変数群‘iug_meteor_bik_uwnd_h2t60min00_10_interpol_stpwrspc’‘iug_meteor_bik_vwnd_h2t60min00_10_interpol_stpwrspc’‘iug_meteor_ktb_uwnd_h2t60min00_10_interpol_stpwrspc’‘iug_meteor_ktb_vwnd_h2t60min00_10_interpol_stpwrspc’> options, ‘tplot変数名’, ytitle =‘MW-bik!Cuwind at 90 km!Cperiod [hour]’
y軸のキャプションを変更する。
> zlim, ‘tplot変数名’, 0, 50 ←カラーバーのスケール調整
28

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●実演習

周期：2-500 時間

東西風90 km

ビアク

コトタバ
ン

29

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●実演習

周期：2-500 時間

南北風90 km

ビアク

コトタバ
ン

30

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●実演習

周期：2-100 時間

東西風90 km

ビアク

コトタバ
ン

31

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●実演習

周期：2-100 時間

南北風90 km

ビアク

コトタバ
ン

32

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●進んだ使用法〇ustrans_pwrspcに許されている引数
選 helpexplains all the keywords and parameters
選 verbose
flags errors and size
必 samplingrate
if set returns array of frequency
選 maxfreq
maximum frequency performing the S transform選 minfreq
minimum frequency performing the S transform選 freqsmplingrate frequency interval
選必 power
returns the power spectrum
選必 absreturns the absolute value spectrum
選 Rremoveedge
removes the edge with a 5% taper, and takes out
least- squares fit parabola
33

5. 流星レーダーデータ解析演習<下部熱圏・中間圏の風速データ解析> [応用課題](3) S(Stockwell)変換解析●進んだ使用法〇ustrans_pwrspc,‘wdc_mag_kak_1min_xyz_x_avg’, /abs ,/sampling,/verbose[以下のようなメッセージがコンソール画面に出力]
Performing S transform:
Not complex data, finding analytic signal.
Minimum Frequency is 0.
Maximum Frequency is 372.
Frequency sampling rate is 1.
The number of frequency elements is 373.
Calculating Local Amplitude Spectra.
〇S変換はPCのメモリをかなり必要とするので、長期間のデータを解析するときには注意を要する
34

6. まとめ●今回の解析講習会では、赤道中間圏・下部熱圏領域の風速変動の統合解析演習を行った。ＴＤＡＳには実装されていない高度なスペクトル解析(コヒーレンス、S変換)を行ったが、これらのプログラムが将来的にUDASに実装される予定である。
●この解析講習会を通じて、各種の分野の異なるデータの統合解析の手法とUDASの使用法を学んだわけであるが、これ以外にもIUGONETプロジェクトでは、多くの超高層観測データ(MLT領域の風速場など)の検索、解析ソフトを公開しているので、皆さんは、ドシドシ使い倒してほしい。
●今後も解析ツールの修正や公開データの種類の増加に伴うアップデートがなされるので、IUGONETのホームページを時折、見て確認してほしい。
35