Thermal diffusivity estimation by Scilab

Nov. 6, 2007

目的と概要

2深度で測定した地温の時系列データから熱拡散率：κ（=熱拡散係数、温度伝導度、thermal diffusivity）を推定する。推定した熱拡散率は、地温プロファイルのモデル化（任意の深さ、時間における地温の推定）、他の物性値の推定（特に土壌水分、気相率）などに利用できる。

推定法自体は新しい物でも何でもない（例えばHeusinkveld, et al., 2004. Agr. For. Meteorol.のAppendixと同じ）。プログラムを作ろうと思った動機はデータ処理の自動化。大量のデータを読み込んで熱拡散率（の経時変化）を出すところまで、一本のプログラムで終わらせるようにする。

手順

測定した2深度のうち、浅い方（z=z0）の地温データ（T_z0）を有限フーリエ級数の形で表す。フーリエ係数は高速フーリエ変換（FFT）を用いて決定する。
有限フーリエ級数で表したT_z0を熱伝導方程式に代入し、任意時間、深さにおける地温の解析解を求める。解析解の中には未知数であるκが含まれている。
深い方（z=z）の地温データ（T_z）に対して解析解をフィッティングし（非線形最小二乗法）、最も当てはまりがよくなるようにκを定める。

熱伝導方程式

${\frac{\partial T}{\partial t}=\kappa\frac{{\partial}^2T}{\partial z^2}}, \hspace{20} \normalsize{ 0<z<\infty \,\, -\infty<t<\infty } \hspace{100pt} (1)$

zは深さ方向に正、と定義。Tは地温、tは時間、κは熱拡散率、ここで必ずしもz=0を地表面とする必要はない。

境界条件と解析解の導出

上部と下部、二つの境界条件を与えると式(1)は解析的にとける。一つの境界条件は簡単で、 $T({\infty,\,t)=\bar{T} \hspace{20pt} (2)$ とする。ここで、 $\bar{T}$ は対象期間中の平均地温で、どの深さでも一定と仮定する。逆に、解析に用いるデータは平均地温がどの深さでも一定とみなせることが必要^※1。

2つ目の境界条件はz=0における地温データ( $T(0,\,t)$ )^※2。解析解を得ようと思ったら、境界条件はなんらかの式で表す必要があるが、実際に得られるデータは、非連続な離散データなので、それをどう連続した関数形で現すか？？が問題。

これにはよく知られている、離散フーリエ変換（Discrete Fourier Transform）というものが使える。地温をいくつもの波（三角関数）の重ね合わせで表現する方法である^※3。解析解はそれぞれの波に対する解の重ね合わせで表現される。この操作を調和解析（Harmonic Analysis）という。

$T(0,\,t)= \bar{T}+\sum _{n=1}^{M} \left{a_n \sin \left(n \omega t \right) + b_n \cos \left(n \omega t \right) \right} = \bar{T}+\sum _{n=1}^{M}\left{A_n\sin\left( n\omega t+\Phi_n \right) \right} \hspace{50pt} (3)$

ここで、 $\omega$ は角周波数（=角速度）で2πN^-1、 $A_n = \sqrt{a_n ^2 + b_n ^2}, \hspace{5} \Phi_n=\arctan \left( \frac{b_n}{a_n} \right)$ である。 $a_n,\,\,b_n$ はフーリエ係数と呼ばれ、これを計算するところでScilabに備わっている高速フーリエ変換関数"fft"が活躍する。これに別に書き記すことにして、まずは(3)式が求まったとして解析解を出してしまう。解析解はラプラス変換を使って求めることができる。ここでは結果のみ記す。(1)式を(2), (3)式のもとで解いた解析解は(4)式で与えられることが知られている。

$T(z,t)=\bar{T}+\sum _{n=1}^{M}\left{A_n \exp \left(-z\sqrt{\frac{n\omega}{2\kappa}} \right) \sin \left(n\omega t + \Phi_n - z\sqrt{\frac{n\omega}{2\kappa}} \right) \right} \hspace{100} (4)$

この解析解を別の深さで測定した地温データにフィッティングさせることにより、κを求めることができる。フィッティングについてはScilabのlsqrsolve関数を使うのが良い。

フーリエ係数（a_nb_n）の求め方-fft関数の利用

Scilabに組み込まれている高速フーリエ変換（Fast Fourier Transform）の関数"fft"を利用して離散フーリエ変換（Discrete Fourier Transform）^※4を行う。fftによって、一発でフーリエ係数ベクトル（ $C_n = \beta_n + i\,\alpha_n$ ）がもとまる！Scilabのfftの仕様だと式(3)との対応では、 $a_n = -\frac{2\alpha_n}{N},\hspace{10} b_n = \frac{2\beta_n}{N} \hspace{10} \small{(1\le n \le N-1)}$ である^※5。また $\bar{T}=\frac{\alpha_0}{N}$ である。ちなみにfftの結果得られたフーリエ係数からExcelの関数をつかって(3)式を表すと、

$T(0,t) = \frac{\beta _0}{N} + \sum _{n=1}^{M} \frac{2 \sqrt{\alpha_n ^2+ \beta_n ^2}}{N} \sin \left( \frac{2\pi tn}{N} + ATAN2 \left( -\alpha_n, \, \beta_n \right) \right)$

ここで0=<t=<N, Mは重ね合わせる波の数でN/2を超えない最大の整数まで計算できる。

※1 地温変化には日変化の他に年変化もあり、位相が異なるので、あまり長い期間を指定するとこの仮定に反することになる。これは要注意。

※2 地表面の温度測定は難しいので、現実的には数mm-数cm程度の深さの地温を使えばよい。

※3 地温のような周期関数は三角関数で現すことができる。とくに無限級数で現した形をフーリエ級数という。地温に対する適用は以下の文献を参照。Horton, R., P.J. Wierenga, and D.R. Nielsen, 1983. Evaluation of Methods for Determining the Apparent Thermal Diffusivity of Soil Near the Surface Soils. Science Society of American Journal. Vol. 47. 25-32.

※4 高速フーリエ変換は2のべき乗データ数を高速にフーリエ級数展開（フーリエ係数を求めることと同値）するものだが、Scilabのfftでは点数が2のべき乗の時は高速なアルゴリズムを、そうでないときは低速なアルゴリズムを使用する。

※5 FFTの正規化係数に決まりはないが、通常はDFTで1、Inverse DFTで1/nとしている場合が多い。Scilabの"fft"このやり方であり、正規化するためにはFFTの得られた振幅を2/N倍する必要がある。またScilabのfftでは、負の符号のiを使用して計算している。

目的と概要

手順

熱伝導方程式

境界条件と解析解の導出

フーリエ係数（anbn）の求め方-fft関数の利用

フーリエ係数（a_nb_n）の求め方-fft関数の利用