数理解析研究所講究録
第 892 巻 1995 年 26-3726

シリコン融液流動による融液温度変動の時系列解析

住友金属工業 (株)未来技術研究所

宮野尚哉(Takaya Miyano)

稲見修一(Shu-ichi lnami)

新谷昭.

(Akira Shintani)

1. 序論

Czochralski 成長法は、外部より電熱ヒーターで加熱され、 かつ、機械的に定速転している f 甘渦内の原料融液表面に種結晶を置き、 この結晶を増塙とは逆方向にさせつつ上方に引き上げることによって単結晶材料を製造する単結晶成長技術である。原料融液は、熱対流と強制対流のために非常に複雑な流動状態にある。単結晶の品質は、融液の流動に強い影響を受ける。結晶成長時の融液流動状態を解析するために、 Navier-Stokes 方程式の数値計算が行われているが [1, 2]、現実の流動状態を定量的に再現することは難しいようである。筆者らは、適切な初期条件、境界条件の下で基礎方程式を解くのではなく、観測されたデータから流動のダイナミックスを関数近似によって再構成し、未来における状態を予測するというアプローチにおいて融液流動を解析することを試みた [3, 4]。系のダイナミックスが非線形で、 カオス的振るまいを示す場合には、 このような手法は特に有効である [5-7]。本報告は、定速i 甘渦内のシリコン融液流動状態を融液温度変動として観測し、 時系列データの非線形予測により流動の複雑さを解析した結果について述べるものである。

27

2. 実験
実験方法の概要を図 1 に示す。

温度変動を白金一白金ロジウム熱電対で測定する。 i 甘渦直径は 406cm、坦渦深さは356 $cm$、融液重量は 65 kg、増渦を 5cm 間隔に並列し、融液表面から約 2 cm の深さで、 増蝸中心、 中心から$5$、$10$ 、 $15$cm 離れた位置に挿入する。熱電対は坦渦ではなく実験室床に固定されているので、 デ$-5$ は回転座標系から見たものではない。融液温度は、 1 $\sec$ 周期で 600 $\sec$ 間測定されるので、温度時系列は 600 点のデータ点から成る。

熱電対

図1

実験装置の概要

28

各熱電対で測定された融液温度変動の 0–200 $\sec$ の部分を図 2 に示す。縦軸は熱起電力であるが、熱電対の出力を増幅する各アンプのオフセット電圧が互いに異るので、 図 2a $-d$ の縦軸は共通の基準値に対する値ではない。 しかしながら、 本研究では変動の乱雑さを解析することが目的であるので、 これは解析結果には影響しない。

$>E$$\frac{o}{arrow 0}Q)$$\geq q)$$\tilde{\frac{arrow 0\in 0}{o}}$$\sim^{\frac{Q)}{=\vdash Q)\frac{\circ E}{\Phi}}}$Time /SeC

Time $/\sec$Time $/\sec$Time $/\sec$$>E$$\veearrow\circ\circ\Phi$$\Phi\geq$$\underline{\tilde{o\in 0}}$$\tilde{\frac{\Phi\circ}{\Sigma\vdash\frac{\in}{\Phi}\circ\Phi}}$

図2(熱電対位置)観測された融液温度a: 増渦中心, $b:5cm$ , $c:10cm$ , $d:15cm$

29

温度変動時系列の自己相関関数を図 3 に示す。 f 甘塙中心では、 規則的変動は観測されないが、 中心から離れると規則的な変動が認められる。変動の周期は・ $b-d$ いずれについても約 28 $\sec$ であり、 増渦Time ag $/sec$

Time lag $/sec$

Time lag $/sec$

Time lag $/\sec$

図3(熱電対位置)

融液温度変動の自己相関関数$a$: 増禍中心, $b:5$ cm, $c:10$ cm, $d:15$ cm

303. 非線形予測モデル$N$個のデータ点からなる離散時系列$\{x(t)\}_{I\cdot 1}N$について考える。 $x(t)$ の変動が因果的であるならば、 即ち、 変動が測定ノイズのような制御不可能な要因に基づくものでな

いならば、 $t+\tau$ の未来における値 $x(t+\tau)$ は、 過去の変動パターンに関する連続で滑らかな関数$F$y(t但し、$\tau$$+$によって近似ことができるであろう。
(1)r) $=F(t, t- l, t- 2, \cdots , t-(D- 1))+\epsilon_{r}$は予測時間間隔、数である。 予測誤差$\epsilon_{r}$相関係数 $\gamma=\gamma(D, \tau)$$D$は埋め込み次元 $[8|$、$\epsilon_{r}$は予測誤差を表わすランダム変の統計的振るまいを、 予測値 $y(t+\tau)$ と実測値 x(t r) との間の、$+$

または、 root-mean-squared error を実測値の標準偏差で規格化して得られる normalized root-mean-squared error $E=E(D, \tau)$ で表現する。 $E=0$
ならば、 予測は完壁であるが、 $E=1$ ならば、 予測は、 時刻$\tau$の平均値であると予測するのと同程度の意義しかもたない。

における値が元の時系列$\gamma$または$E$の挙動を解析すると、 時系列のダイナミックスに関する情報が得られる [5-7]。例えば、 白色ノイズに汚染された正弦波は一見複雑に変動しているように見えるが、 長期予測が可能で$\gamma,$$E$は$\tau$に依存せず一定値をとる。一方、 時系列が少数自由度のカオスならば、短期予測のみ可能で、$\gamma$は$\tau$の増加に伴って 1 に近い値から急速に に減少するであろう$0$[5]。この場合、 系の自由度は、$\tau$を固定した時にる。 カオス時系列の最大リアプノフ数ができる $[7|$$\lambda>0$

は、$E$$\gamma$が最大となるようなの$\tau$$D$に相当すに対する挙動から推定すること。

(2)$E(\tau)\propto E(1)\exp(\lambda\tau)$

パワースペクトルが周波数 に対して伊 ($f$$\alpha$はパワースペクトル指数) で特徴づけら31れる自己相似ランダム時系列 (有色ノイズまたは非整数ブラウン運動とも呼ばれる)の場合には、 カオス時系列であるかのように$\gamma$が$\tau$の増加に伴って減少するが、 (2) 式は成立しない [9, 10]。これは以下のようにして示すことができる [10]。時系列が次の性質をもつならば、 自己相似ランダムである。
$\tau^{-H}<x(t+\tau)-x(t)>M=<x(t+1)-x(t)>M$
ここで、

$H$

はスケーリング係数で $\alpha=2H+1$、(3)$<>M$ は統計モーメントであり、 (3) 式は統計的性質が等しいことを表わす。 もし、 予測モデル$F$が、 自己相似ランダム時系列の自己相似性をよく再現しているならば、$\tau^{-2H}<\{[y(t+\tau)-x(t)]-[x(t+\tau)-x(t)]\}^{2}>$$=\tau^{-2H}<\{y(t+\tau)-x(t+\tau)\}^{2}>$(4)$=<\{y(t+1)-x(t+1)\}^{2}>$従って、$E(\tau)=E(1)\tau^{H}$(5)

即ち、 $\log T^{-\log[E(\tau)/E(1)]}$ プロットが線形相関をもつならば、 時系列は自己相似ランダムであって、 その傾きはスケーリング係数に相当する。一方、 $\tau-\log[E(\tau)/$$E(1)]$ プロットが線形相関をもつならば、 時系列はカオスであり、 その傾きが最大リアプノブ数である。 (5) 式は、正弦波に自己相似ランダムを重畳した時系列についても成立する [10]。

近似関数$F$を構成するには様々な手法があるが、近似関数がもつべき性質は、与えられたデータの因果性をよく再現できるように a priori に決定しなければならない [11,12]。例えば、 ニューラルネットワークによって$F$を構成することができる [3, 4, 11$-$

3214]。しかしながら、 ニューラルネットワークは、 その構造を最適化するのに比較的長時間を要する。 ここでは計算時間を節約するために Sugihara と May によって開発されたシンプレックス投影法 [5] (Nonlinear regression of a simplex projection method)
を用いて$F$を構成する。 この非線形

LIBRARY, PREDICT の 2 グループに分ける。各グループについて空間のベクトルを支

$=$

$D$

次元ユークリッド$(x(t), x(t- 1), x(t- 2), \cdots, x(t- (D- l )))$ とし、xt に対するターゲ$*$

ットを $x(t+\tau)$ と決める。予測は、 LIBRARY グループのデータを参照して PREDICTグループについて行なわれる。 x p$\tilde$$y_{\mu\tau}$$\in$

PREDICT に対して$\tau$

時間 ス$\overline{\tau}$

ップ未来の予測値

は、$y_{\mu\tau}=Z_{k-1^{D+1}}x(i(k)+\tau)\exp(-d_{ip})/Z_{k\cdot 1}^{D+1}\exp(-d_{ip})$(6)$d_{ip}=|\dot{\tilde{x}}_{\langle k)}$(7)

一$\vec{x}_{p}|$と与えることにする。 ここで、頂点にあって、$\tilde{X}_{\beta}$$\overline{x}_{tk)}^{*}\in$

LIBRARY は荒$\in$

PREDICT を取り囲む多面体のに (7) 式で与えられる距離の意味で最近接の $D+1$ 個のベクトルである $(k=1, \cdots, D+1)$ 。非線形予測により時系列のダイナミックスを解析する手法の顕著な利点は、$G$

rassg rg r-P cacci a のアルゴリズム [15] のようなフラクタル次元を$e$$e$$re$計算する手法に比べて少数のデータ点数から信頼性の高い診断を下すことができることである。 このような長所は、近似関数の滑らかさを a priori に決めたことによってもたらされたと考えられる。

4.

予測結果および考察

シンプレックス投影法を予測モデルに用いて、 図 2a–$d$の温度時系列を予測す33る。 時系列の前半 300 点を LIBRARY グループとし、 後半 300 点を PREDICT グループとする。予測値と実測値の相関係数を埋め込み次元の関数として表わした結果を図 4 に示す。相関係数は、埋め込み次元の選択に余り影響されない。融液の自由度は実際上無限大とみなしてよいので、 これは予想された結果である。熱電対の位置の相違による差はほとんど認められない。

$\underline{\sim\subset\frac{\Phi}{o}}$$\vee\check{\Phi\circ O}$$\underline{\subset 0}\tilde{\varpi}$$\overline{o\Phi 0=}$

Embedding dimension

図4

温度時系列の非線形予測: 埋め込み次元依存性(熱電対位置).

$\bullet$: 増渦中心, $X:5$ cm,$\triangle:10$cm,$\square :15$cm

図 5 は予測値と実測値との相関係数を予測時間間隔に対してプロットした結果を表わす。予測方法は、 図 4 と同様であるが、 埋め込み次元は $D=10$ と固定されている。堀
渦中心付近では、 相関係数は時間$\overline{\tau}$ップの増加に伴って減少し、 10 $\sec$ 後には完全に予測不可能になる。一方、 増渦中心から 10cm 以上離れると、 5 $\sec$ 後までは相関係数は減少するが、 それ以後ほぼ一定値をとり、 ある程度長期予測可能となる。温度時

34

系列がカオスであるかどうか判定するためには、 予測誤差を $E(\tau)$ で表現し、 予測時間間隔$\tau$に対して semi $-\log$ または log–log プロットをとればよい。

$\underline{\vee\subset\frac{\Phi}{o}}$$\vee\check{\Phi oO}$$\frac{\underline{\subset 0}\tilde{\varpi}}{\Phi,=}$$Oo$Time step $/\sec$
図5

温度時系列の非線形予測: 予測時間間隔依存性(熱電対位置)$\bullet$: 増蝸中心, $X:5$ cm,$\triangle:10$

cm,$\coprod;15$cm

図 6 と図 7 は、 それぞれ、 予測時間間隔と予測誤差の semi– $\log$ および log–log プロットである。埋め込み次元は、 $D=10$ である。 プロ、ットの線形相関は、 semi $-\log$ プロットよりも、 $\log-\log$ プロットのほうが良い。即ち、 融液温度は乱雑に変動しているが、 カオスではなく自己相似ランダムである。最小自乗法によって得られた $\log-\log$

プロットの傾き、即ち、 スケーリング係数を表 1 に示す。

35$-$$\check{t1J}-$.

$-$$\}$$\check{\mu\lrcorner}$$\mathring{\underline{o}}’$$\tau/\sec$

図6

温度時系列の非線形予測: 予測間隔依存性 (semi $-|og$ プロット)(熱電対位置)$\bullet$: 柑渦中心,$\cross:5$cm,$\triangle:10$

cm,$\square :15$

cm

$-$

量
$-$
$b$

$\check{\coprod J}$$\underline{\overline{\mathring{o}’}}$$\log r$

図7

温度時系列の非線形予測: 予測時間間隔依存性 (lOg-lOg プロット)(熱電対位置)$\bullet$: 増塙中心, $X:5$ cm,$\triangle:10$

cm,$\square :15$

cm

36

表1

温度変動の自己相似性におけるスケーリング係数の推定値

熱電対位置

スケーリング係数 Hd2

ま甘渦中心0.3160.999$5cm$0.493

0.999

10cmO.4370.991

15cmO.4330.994

表 1 において、$7^{d^{2}}$はプロットの線形相関、 即ち、 図 7 のデータ点の直線への fitting の程度を表わす correlation of determination である。増塙中心付近で、$H$が小さいことは、 f 甘渦の中心に向かうほど温度変動がより乱雑になること反映している。 自己相関関数を考慮すると、 融液の流動は、 増蝸中心付近では自己相似ランダムであり、 中心から離れると、 自己相似ランダムに波動を重畳したような状態にあると推測される。

融液の流動は、 どの部分も時間と温度変動を適切にスケーリングしない限り、 平均値や分散が時間に依存して変動するという意味で、 非定常性を有している。 このような状態が Navier $-$ Stokes の方程式から導出され得るのかどうか筆者らにはわからない。 自己相似ランダムは、 地球気温変動についても認められるようである [16]。多数自由度系である流体の一つの特徴を表わしているのかも知れない。

参考文献[1] J. R. Ristorcelli and J. L. Lumley, . $Cr\gamma sta/Growth,$$d$$\underline{129},249$(1993).

37[2] K. Kakimoto, M. Watanabe, M. Eguchi, and T. Hibiya, J. $Cr\gamma sta/Growth,$$\underline{126},435$(1993).
[3] T. Miyano and A. Shlntani, Appl.

$Ph\gamma s$. Lelt.,$\underline{63}$, 3574, (1993).

[4] T. Miyano, H. Morita, A. Shintani, K. Kanda, and M. Hourai, submitted to $f$.
$Ph\gamma s$$A\rho p/$.

.

[5] G. Sugihara and R. M. May, Nature, 344, 734 (1990).

[6] J. D. Farmer and J. J. Sidorowich,[7] M. Casdagli, Physica $D,$$\underline{35},335$$Ph\gamma s$. Rev. Len., 59, 845 (1987).

(1989).

[S] F. Takens, Lect. $Notes/nMathemat/cs,$[9] A. A. Tsonis and J. B. EIsner, Nature,$\underline{898}$, pp.366-381 (1981).

$\underline{358},217$(1992).

[10] T Miyano, submitted to Nature.
[11] T. Poggio and F. Girosi, A./. Memo No.1140 (Artificial lntelligence Laboratory,

Massachusetts lnstitute of Technology, 1989).
[12] T. Poggio and F. Girosi, Proceedings of lEEE, 78, 1481 (1990).

[13] D. E. Rumelhart, J. L. McCIelland, and the PDP Research Group, $Para//e/$Oistnbuted Processing, pp.31 S-362 (MIT Press, Cambridge, 1986).
[14] T. Miyano and F. Girosi, A./. Memo No.1447 (Artificial lntelligence Laboratory,

MaSSaChuSetts lnstitute of Technology, 1994).
[15] P. Grassberger and 1. Procaccla, $Ph\gamma sicaD,$[16] T. Miyano and F. Girosi, unpublished.

$\underline{9},189$(1983).