Transcript 基礎数学 IV

2007.05.08-09
基礎数学 IV
0. 揺らぎとその解析例
１．周波数解析
２．データ処理における問題
３．高次スペクトル解析
参考文献：日野幹雄 スペクトル解析
杉山高一 多変量データ解析入門
榊原進 ウエーブレット
Fluctuation & Self-Organization
• Key words
2
J. Mather
G. Smoot
2006 Nobel Prize
宇宙からのマイクロ波の背景光の揺らぎ分布 => 宇宙の構造形成
3
NASA 2006
In ten sity (a.u .)
Visible/IR emission+ BB radiation
6 0 00
(a)
M OS
Plasma light emission
H
4 0 00
Spectrum analysis
M oI
2 0 00
H
•Line emission
BB
H
0
20 0
4 00
600
8 00
•Continuum emission
In ten sity (a .u .)
3 00
2 50
•BB emission
2 00
(b )
IR S p e ctro m e te r
1 50
1 00
In ten sity (a .u .)
10 00
1 20 0
1 40 0
1 6 00
Q 3+M oI
48 0
Q1
M oI
44 0
40 0
Q2
(c)
ACTON
Q4
Q5
36 0
60 0
6 01
60 2
6 03
6 04
6 05
W a velen gth (n m )
6 06
4
Raju, submitted to NF 2007
マイクロ波によるＢＢ連続放射の観測
密度揺らぎ計測
NASA
プラズマからの揺らぎ計測シス
5
テム（Wataya, Idei )
ｐｌａｓｍａ 分布の瞬時の変化
Ｘｒａｙカメラ
画像解析
#72402
2000
I
sx
( a .u .)
3000
sx128
sx131
1000
385
386
t(ms)
387
388
6
中心の突出した分布から
中空の分布へ突然変化
トーラスプラズマのエネルギー配列が突然変化
秩序から秩序への飛び
7
t=427
429.52
429.84
430.96
#72464
Z (mm)
Z (mm)
2D reconstructed emissivity at the internal disruption
m=2
DR (mm)
8
H. Zushi PPCF 1997
秩序の飛びと構造の爆発的成長
２ｍｓ
3
2
m=2 amplitude
p'(r/a=0.1)
2
1
0
-1
1
１ｍｓ
-2
t（ｍｓ）
385.0
構造の成長
386.0
387.0
385.0
386.0
387.0
9
time(ms)
time(ms)
秩序形成の過程
log(m=2,r=rs)
s
構
造
振
幅
lo g ( m = 2 ,r = r )
1
crash
0.5
before
crash
0
after crash
-0.5
3
5
p'(r~r )
s
中
空
秩
序
7
中
心
集
積
10
構造の間欠振動の画像解析例
m＝２成分の間欠的発達
m=2 amplitude
ｒ
time
中空分布
r
2.0
2.0
After crash
1.5
1.5
1.0
1.0
0.5
0.5
0.0
0.0
388
389
Beam pulse end
404
405 ms
11
lo g (
構
造
0
m =2
)
秩序と構造
-1
4
6
中
空
P'(q=2)
秩
序
8
中
心
集
積
10
何
が
中
心
部
へ
の
集
中
を
駆
動
す
る
か
12
Fluctuations
X  X  X
任意の物理量Xの時間発展を考える。
ここでＸは平均値Xの周りをXの変動の広がりをもって揺らぎながら変動する。
Define a quantity which represents the average range
of fluctuation!
X
X
？
13
Root Mean square (r.m.s.)
  X   
2
 X  X  
2
 X 
2

X  X 
2
 X  2XX  X
2
2
 X
2
 2XX  X
 X
2
  X   
2
X
 X
2
 X
2
2
2
14
揺らぎの解析方法
• フーリエ解析 （周期振動、相関、スペクトル）
• 多変量解析 （多変数の相関、主成分分析）
• ウエーブレット解析（バースト的な信号の解析）
15
フーリエ周波数解析
•
•
•
•
自己相関関数 auto correlation
パワースペクトル power spectrum
相互相関
cross correlation
クロススペクトル cross spectrum
基礎となるのはどのような周期関数も三角関数の級数和
で表すことができる というフーリエ理論である。
F (t ) 

a n cos( n  t )  b n sin( n  t )
n
16
時間-空間から周波数-空間へ
（実際には線スペクトルへの分解ではなく、連続的に変化する
周波数スペクトルへの変換）

x(t ) 

 X ( f )e
i 2  ft
df 

 X ( )e
i t
d


X (f ) 
 x ( t )e
 i 2  ft
dt ,

X ( ) 
1
2

 x ( t )e
 i t
dt ,

•特定の時間幅に存在する非周期関数を周波数の連続スペクトル
（単一の周波数成分が多数重なり合っている）で表現する。
•X(f)dfは周波数f~f+dfの波の振幅 全体の分布を周波数スペクトル
17
相関 correlation
解析の一つの目的は
１）変動が意味のある周期振動か？
２）変動が他の因子によって駆動され
たものか？
３）変動が他の因子に影響を与えてい
るものか？
４）変動が他の因子と相関があるもの
か？ などを調べることである。
r
E ( xy )
2
2
E ( x )E ( y )
rを相関係数、Eをアンサンブル平均
（母集団平均）という。
18
自己相関 self correlation
変動の周期性を調べるには、信号x(t)と少し時間をずらした
同じ信号x(t+tの相関を調べればよい。
C (t )  E [ x ( t ) x ( t  t )]
この場合だと時間平均のために
時間積分して、時間の間隔でわる。
自己相関係数
R (t )  C (t ) / C ( 0 )
19
周期関数の相関関数は
周期関数
ノイズ成分が大きくなると
周期関数の周期性はみる
ことができにくい
ノイズ成分の性質によって
相関時間が変化する。
無相関ノイズは関数となる
20
揺らぎの寿命と自己相関
相関関数の対称性（偶関数）
C ( t )  C (t )
不規則振動の性質の指標
C (t )  0 , t  
R e

t
tc
ここで、tcを相関時間といい、不規則振動の時間スケールの
目安である。たとえば揺らぎの寿命に相当する。
21
Power spectrum
S ( ) 
lim
T
 2

E
X ( ) X * ( ) 
 T

Eはensemble average
X()は複素フーリエ成分
Ｔは観測時間幅
Ｓは実数の偶関数
Wiener-Khintchine formula

C (t ) 


S ( )e
i t
d
S ( ) 
1
2

 C (t )e
 i t
dt

22
例１）1次Markov 過程
f ( t  D t )   f ( t )  N random ( t )
C (t )  E f ( t ) f ( t  t ) ,
C (t  D t )  E f ( t ) f ( t  t  D t ) 
 E f ( t )   f ( t  t )  N random ( t  t ) 
  E f ( t ) f ( t  t )   E f ( t ) N random ( t  t ) 
  C (t )
C (t  D t )  C (t )  C ' (t ) D t
  1
C ' (t )
C (t )

C ' (t )
C (t )
1 
Dt
C (t )  C ( 0 ) e
Dt,
  ,
 t
23
１次マルコフ過程のpower spectrum
S ( ) 


1


1
2
 C ( 0 )e
 t
e
 i t
dt


 C ( 0 )e
 t
cos t d t
0
C(0)


 
2
2
24
25
26
“1/f noise spectrum”は何を意味しているか
27
J. Appl. Phys. 1973
S (f )  f
1
E [ x ( t ) x ( t  t )]  t
0
28
揺らぎの中で大きなパワーを持つ
成分を調べる。
0
4
rfstray1_FFT_504184_3
? rfstray1
rfstray1_FFT_5041843
rfstray1_FFT_504188_3
rfstray1_FFT_5041903
rfstray1_FFT_504190_3
rfstray1_FFT_5041883
rfstray1_FFT_504192_3
8
frequency[Hz]
12x10
RFstray(a.u.)
#504192:0-0.05(0.0064)
3
RFstray(a.u.)
100
80
60
40
20
0
RFstray(a.u.)
power(rfstray)
RFstray(a.u.)
マイクロ波の漏洩強度計測
3.0
#504192
?
RFstray
2.0
1.0
0.0
3.00
10
20
30
40 50x10
#504190
?
-3
RFstray
time(s)
2.0
1.0
0.0
3.00
10
20
30
#504188
40 50x10
?
-3
RFstray
time(s)
2.0
プ
ラ
ズ
マ
密
度
の
増
加
1.0
0.0
3.00
10
20
#504184
30
time(s)
40 50x10
?
-3
RFstray
2.0
1.0
0.0
0
20
40
time(s)
-3
60x10
29
剛体球の弾性散乱の計算機シミュレーション （例）
（希薄な場合はBoltzmann 分布)
速
度 C(t)
相
関
関
数
低密度
拡
散
係
数
衝突時間で規格化
した時間
高密度
Ｔ
 t
2
~
~
~
C ( t )  v ( 0 ) v ( t )  v  exp   
 t 
~
~
d v (t )
v (t )

（１次でTaylor expansion)
dt
t
Alder PRL 1967
30
C(t)-CB(t)
高密度では長時間記憶が保持
衝突時間で規格化した時間
31
複数の信号間の相関
Cross correlation
C xy (t )  E [ x ( t ) y ( t  t )],
R xy (t )  E [ x ( t ) y ( t  t )] /
 C xy (t ) /

E ( x(t )x(t  t ) E ( y (t )y (t  t )

C x ( 0 )C y ( 0 )
32
相互相関の特性
一般にlag time=0に対して
非対称
33
Cross spectrum
 2

S xy (  )  E 
X * ( ) Y ( ) 
 T

S xy (  ) 
1
2

C xy (t ) 


C
xy
(t )e
 i t
dt

S xy (  ) e
i t
d

X,Y信号のうち相関の強い周波数帯の強度分布が解析可能。
34
Cross phase
 2

S xy (  )  E 
X * ( ) Y ( ) 
 T


2
T

E X ( ) Y ( ) e
 i xy (  )

ここで、xyは信号xおよび信号yの位相角の差を表す。
xy(＝ x(- y(
Ｘ、Ｙ揺らぎの信号がお互いに位相がずれていることを意味する。
35
密度と電位の揺らぎ
密度の空間分布をBoltzmann 分布とすると、電位の揺らぎによる
~
~
密度分布は
U (r )
e ( r )
n e  n 0e

T
 n 0e

T
となり、電位揺らぎが運動エネルギーに比べて十分小さい近似の下
では密度揺らぎは
~
e
~
ne  n 0
T
と書ける。
一方、電位変動による磁場に垂直方向の流体のドリフト速度は


~
 ~ 
E  B t
V drift 

ir
2
Bt
Bt
36
揺らぎに駆動された拡散束
~  v  ~v   n v  n ~
~v  n
~ ~v  n ~v  n
~ ~v
  n 0  n
v

n
0
0 0
0
0
0
~ ~v 
E    E n
平衡状態では流れがない
~
i  nv
~
~
~
~
nv  n v e
 n v cos(  nv )
もし、密度揺らぎと径方向速度揺らぎが90度位相のずれがあれば、
いくら揺らぎ振幅が大きくても、正味の拡散束は0となる。
すべての周波数帯の積分から、正方向の拡散束が得られれば、
“観測位置においては”、揺らぎによって粒子が“磁場のかご”から
漏れだしていると結論できる。
１）磁場のかごのすべての位置からの正味の損失を意味するもの
ではない。
37
揺らぎのcoherence
 ( ) 
2
 ( ) 
S xy (  )
2
S xx (  ) S yy (  )
,
S xy (  )
S xx (  ) S yy (  )
 2

S xy (  )  E 
X * ( ) Y ( ) 
 T


2
T

E X ( ) Y ( ) e
 i xy (  )
が１に近いということは二つの信号が極めて強い相関のある
揺らぎであることを意味している。
38

降雨と河川流出量の相互相関
一日程度の時間遅れ
流出量＝地表面流出＋
時間遅れのある滲透流成分
39
降雨と流出量のcoherence
f=0.25 day-1 でcoherenceに差がある。
１）即ち四日周期以下は表面流出と考えられる。そのピークはほぼ半日周期であり、
降雨後、６時間後に最大である。
40
-1
２）浸透流出は0.325,0.42d 程度の周期性を有す。
実空間と波数空間
１
二点の観測点が距離Lだけ離れていれば、
coherentな波動の揺動振幅とその周波数に
おける位相差から、
２
 U 
L 12 
 12
1

L 12 [ m ]f [ s ]
 12 [ rad ] / 2
揺動の伝播速度が求まる。
41
2-2b Ha の応答時間遅れ
（t=0:密度上昇時と定義）
Ha摂動のＦＦＴ解析
１）Ｐｏｌ 断面では
内側ー外側に変動振幅のピーク
基本周期は1.5 Hz
1.4
1.2
15.5
16.0
16.5
Poloidal arrayTime(s)
#86203 @1.7142[Hz]
ha1__FFT
ha2__FFT
ha3__FFT
ha4__FFT
ha5__FFT
ha6__FFT
ha7__FFT
IHa(a.u.)
power(Ha)
#86203
800 900 [mm]
0
2
4
6
frequency[Hz]
8
10
1.8
1.6
1.4
1.2
1.0
0.8
15.0
0.80
0.70
0.60
1.0
0.8
15.0
1000
800
600
400
200
0
'1D_r'
'2D_r'
'3D_r'
'5D_r'
'6D_r'
'7D_r'
'9D_r'
ne
Ne(1E19m-3)
Intensity (a.u.)
Toroidal array
0.50
17.0
0.80
0.70
'9A'
'9B'
'9C'
'9D'
'9E'
'9F'
'9G'
ne
15.5
16.0
16.5
0.60
0.50
17.0
Time(s)
閉じ込め緩和振動に伴うHa変動は同一トーラス位置で内外位相遅れを保持した
42
まま、トーラス全体のHaを変動させている。
揺らぎの伝播速度
D(degree)
80
60
40
20
0
-20
#85706-766
0
f~3-4 Hz
の決定例
85739
85737
85733
85735
85736
85738
85706
85760
85762
85766
1
2
3
toroidal distance[m]
4
-3
#85706-766
f~3-4 Hz
Dt(s)
50x10
40
30
20
10
0
0
1
2
3
toroidal distance[m]
85739
85737
85733
85735
85736
85738
85706
85760
85762
85766
4
43
空間モード数の決定（例）
４
５
円柱プラズマの周囲にポロイダル磁場を
計測するコイルを配置する。
３
２
６
E  
１
７
８
１１
９
１０
B
t
V  N
コイルの巻き線数N
磁束の時間変化
あるいは特定の周波数
に比例した電圧が発生
プラズマ中に揺らぎの電流密度が存在すると
~
~
  B  0 j
ある特定の周波数に対して、
im 
e
44
のポロイダル構造を持つ場合について考える。
低位のモード
１）m=0
膨張と収縮
2) m=1 cos
Horizontal shift
bent
3) m=2
cos2
４）m=3 cos3
楕円変形
△変形
45
モードの決定
１）coherentな周波数 ｆｃ を探す
２）基準のコイルに対して、cross phaseを fcの場合について
計算する。
11,12, 13………… 1n
３）コイルの設置角度を基準に、観測位相差を表示し、最小２乗
fitを行う。 m=infinite
位
相
差
m=2 m=1
1
コイル位置
n
m=0
46
モード数の最大値は？
• 設置可能なプローブ数N（たとえば均等配置）
と決定可能なモード数の関係は？
47
48
２．データ処理の問題点
•
•
•
•
•
０）digital 計測（サンプリング時間）に伴う誤差
１）有限観測時間と周波数解析
２）信号のdrift除去
３）アンプの周波数特性と周波数解析
４）FFT
49
1. 有限長のデータ列と分割平均
１） データ列の長さＴで決まるより低周波数(1/2T)は解析できない
2） データ列の長さＴを長さsの分割で解析する場合は1/2sより高
周波数は解析できない
50
Power (kJ/s)
Ultra Low Frequency Events
2.5
2
1.5
1
0.5
0
 Very low frequency semiregular oscillations are found
in signals on heat loads,
particle recycling, and
impurity influx and contents.
Antenna loss
ML
MoXIII(a.u.) H(a.u.)
50
 Frequency ~ 1-210-3 Hz
amplitude ~a few %– 100 %
40
30
20
1.0
0.8
0.6
0.4
0.2
0.0
 During the last ULF event, the
five hour discharge terminated.
0
5
10
time(s)
15
20x10
3
51
IH a.u.)
130
125
120
115
110
11.6
12.0
12.4
12.8x10
3
time(s)
IH a.u.)
124
120
116
112
11.60 11.65
IH a.u.)
130
11.70 11.75
11.8011.85x10
3
time(s)
trise=0.2s,tdecay=0.7s
125
120
115
110
11.748
11.750
11.752 11.754x10
time(s)
3
52
0. Aliasing in digital sampling
離散的なデータ点に基づくスペクトルには真のスペクトルと
２ｆＮ＝Dtの整数倍だけ異なった周波数が同一視される。
~
Pd ( f )  P ( f )  P ( 2 f N  f )  P ( 2 f N  f )  P ( 4 f N  f )  P ( 4 f N  f ) 
fN 
1
2Dt
Nyquist frequency
53
２．Drift 除去
低周波のdriftがある場合はあらかじめ
その成分を除去しておく
x ( t )  U ( t )  U drift ( t ),
U drift ( t )  a  kt  U 1 
U 2  U1
T
54
t
3. Amp の周波数応答
detector
Amp
Frequency response
X(t)
X(t)W(f)
ADC
fN=0.5 fadc
W(f)
-3dB
f
fN
55
4. FFT
1. N=2p 整数で且つ２のべき乗に選ぶ、たとえば、N=64, 128
…,足らなければ ０を足す。多ければ除く
２．時間を切り取るので滑らかなデータwindowを掛ける。
適切に処理することにより、負のスペクトル裾野の出現を押
さえる。（最後に強度<x2>が減ったことを補正）
３．スペクトルの平滑化
３－１）統計平均 （同様な条件下でのデータを平均する）
３－２）分割平均（１つのデータ長をｓ個の部分に分割し
各の区間のFFTスペクトルの平均をとる）
N=4096 s=64とすれば64個のスペクト
ルの平均となる）
３－３）周波数平均（スペクトルwindow Ｗ（ｆ）を用いてFFT
スペクトルを平滑化）
56
4. FFT II
4. Coherence :１個のデータのみからこれを計算すると、
(f)=1となる。従って、Ｘ（ｆ），Ｙ（ｆ）はアンサン
ブル平均（前述）か周波数windowによる平滑
化の処理を行った後、coherenceを計算する。
57
練習問題 Ｉ
• 電流密度と温度分布の逆相関の緩和振動が観測さ
れた。
• プラズマの熱輸送に電流密度分布が重要な役割を
果たしているという作業仮説を立てる。
• 電流分布が中心ピーク形状の時内部の熱輸送は大
きくなり、中心から熱が外部に運ばれる。そのため
中心にピークを持つ温度分布は緩和し、中心部が
平坦な分布となる。この温度分布は電流密度分布
に影響を及ぼし、電流密度分布を平坦化する。その
結果、再び熱輸送が減少するので、中心部の温度
は上昇し中心ピークの温度分布が形成される。
• こうして観測された電流密度と温度分布の緩和振動
が説明できる。
58
59
3
2
ne
dT e
dt
   n e  e  Te   Prf  Poh
j
j 

   V  B        
t
 

０次元の非線形モデル方程式
dTe
dt
dj
dt
  T Te    T Te j
  j j    j Te j
60
61