関数create_timefunc_pow34 マニュアル

(The documentation of function create_timefunc_pow34)
Last Update: 2021/12/8


◆機能・用途(Purpose)

時間関数pow3-4(下記(\ref{eq.definition})式)に基づく時系列データを作成する。 エラーチェック機能が不十分なのでこの関数を直接用いてはならない。 この関数は関数create_timefuncからの内部呼び出し用である。
Create a time series data based on a time function pow3-4 defined by eq. (\ref{eq.definition}). Do not directly call this function because error checks are insufficient. This function is for an internal call from function create_timefunc.


◆時間関数の定義 (Definition of the time function)

概要(Overview)
時刻\(t<0\)および\(t>\tau_p\)で0、最大値が1、 かつ全ての時刻で2階導関数までが連続となる関数のうち、 時刻\(t=\tau_p/2\)に関して非対称な最小次数の多項式である。
The polynominal with minimum possible degrees that satisfy the following conditions: zero for \(t<0\) and \(t>\tau_p\), unity at the maximum, second order derivative of the function is continuous at all the time, and asymmetric with respect to time \(t=\tau_p/2\).

定義(Definition)
\[\begin{equation} f(t)= \begin{cases} 0 & (t<0) \\ C\left(\frac{t}{\tau_p}\right)^3\left(1-\frac{t}{\tau_p}\right)^4 & (0\leq t \leq \tau_p) \\ 0 & (\tau_p<t) \end{cases} \label{eq.definition} \end{equation}\] \[\begin{equation} C=\left(\frac{7}{3}\right)^3\left(\frac{7}{4}\right)^4 \label{eq.amplitude} \end{equation}\]
以下の微分形と積分形の導出にあたっては \(0\leq t \leq \tau_p\)の範囲において(\ref{eq.definition})式が \[\begin{equation} f(t)=\frac{C}{\tau_p^7}t^3(t-\tau_p)^4 =\frac{C}{\tau_p^7}(t^7-4\tau_pt^6+6\tau_p^2t^5-4\tau_p^3t^4+\tau_p^4t^3) \label{eq.expansion} \end{equation}\] と展開できることを利用する。
In the derivations of differential and integral forms below, eq. (\ref{eq.expansion}) is used as an expanded form of eq. (\ref{eq.definition}) for the range \(0\leq t \leq \tau_p\).

1階導関数 (1st order derivative)
\[\begin{equation} \frac{df(t)}{dt}= \begin{cases} 0 & (t<0) \\ \frac{C}{\tau_p^7} (7t^6-24\tau_pt^5+30\tau_p^2t^4-16\tau_p^3t^3+3\tau_p^4t^2) & (0\leq t \leq \tau_p) \\ 0 & (\tau_p<t) \end{cases} \label{eq.derivative} \end{equation}\]

1階積分 (1st order integral)
\[\begin{eqnarray} F_1(t) &=& \int_{-\infty}^t f(t’)dt’ \nonumber \\ &=& \begin{cases} 0 & (t<0) \\ \frac{C}{\tau_p^7}\left(\frac{1}{8}t^8-\frac{4}{7}\tau_pt^7+\tau_p^2t^6 -\frac{4}{5}\tau_p^3t^5+\frac{1}{4}\tau_p^4t^4\right) & (0\leq t \leq \tau_p) \\ \frac{C}{280}\tau_p & (\tau_p<t) \end{cases} \label{eq.integral1} \end{eqnarray}\]
ここで積分定数は\(t<0\)で\(F_1(t)=0\)、 \(t=0\)および\(t=\tau_p\)で\(F_1(t)\)が連続となるように与えた。
Here, the integration constants are taken so that \(F_1(t)\) is zero for \(t<0\) and continuous at \(t=0\) and \(t=\tau_p\).

2階積分 (2nd order integral)
\[\begin{eqnarray} F_2(t) &=& \int_{-\infty}^t F_1(t’)dt’ \nonumber \\ &=& \begin{cases} 0 & (t<0) \\ \frac{C}{\tau_p^7} \left(\frac{1}{72}t^9-\frac{1}{14}\tau_pt^8+\frac{1}{7}\tau_p^2t^7 -\frac{2}{15}\tau_p^3t^6+\frac{1}{20}\tau_p^4t^5\right) & (0\leq t \leq \tau_p) \\ \frac{C}{504}\tau_p^2+\frac{C}{280}\tau_p(t-\tau_p) & (\tau_p<t) \end{cases} \label{eq.integral2} \end{eqnarray}\]
ここで積分定数は\(t<0\)で\(F_2(t)=0\)、 \(t=0\)および\(t=\tau_p\)で\(F_2(t)\)が連続となるように与えた。
Here, the integration constants are taken so that \(F_2(t)\) is zero for \(t<0\) and continuous at \(t=0\) and \(t=\tau_p\).

この関数の主な想定用途 (The main assumed purpose of this function)
この関数は波形インバージョンで用いるために考案したものである。 波形インバージョンでは観測波形のスペクトル\(U(\omega)\)が 震源におけるモーメントの時間関数のスペクトル\(M(\omega)\)と グリーン関数のスペクトル\(G(\omega)\)の積として \[\begin{equation} U(\omega)=M(\omega)G(\omega) \label{eq.inversion.obs} \end{equation}\] と書けることを利用する。 グリーン関数はインパルス応答であるが、 これを数値的に求めることは不可能であるので 適当な震源時間関数を与えて計算した理論波形で代用する。 計算に用いる震源時間関数のスペクトルを\(S(\omega)\)、 理論波形のスペクトルを\(G’(\omega)\)とすれば \[\begin{equation} G’(\omega)=G(\omega)S(\omega) \label{eq.Gdash} \end{equation}\] である。(\ref{eq.inversion.obs})式の両辺に\(S(\omega)\)を掛けて (\ref{eq.Gdash})を代入すると \[\begin{equation} U(\omega)S(\omega)=M(\omega)G’(\omega) \label{eq.inversion.obs.Gdash} \end{equation}\] という形が得られる。波形インバージョンでは
  1. 適当な\(S(\omega)\)を与えて\(G’(\omega)\)を計算し、
  2. \(U(\omega)\), \(S(\omega)\), \(G’(\omega)\)を用いて (\ref{eq.inversion.obs.Gdash})式を解くことで\(M(\omega)\)を求める
という手順で解析を行うことになる。 このときの\(S(\omega)\)に相応しい関数として考案したのがこのpow3-4である。
This function was developed for waveform inversion, which uses eq. (\ref{eq.inversion.obs}), where \(U(\omega)\), \(M(\omega)\), and \(G(\omega)\) are the Fourier spectra of the observed waveform, the time function of the moment at the source, and the Green function, respectively. The Green function is defined as an impulse response, which cannot be computed numerically. Instead, a theoretical waveform calculated with a given source time function is used. The Fourier spectrum of this theoretical waveform, denoted as \(G’(\omega)\), is related to the Green function as eq. (\ref{eq.Gdash}), where \(S(\omega)\) is the Fourier spectrum of the source time function used in the calculation of the theoretical waveform. Multiplying \(S(\omega)\) to the both sides of eq. (\ref{eq.inversion.obs}) and inserting eq. (\ref{eq.Gdash}) results in (\ref{eq.inversion.obs.Gdash}). The procedure of waveform inversion is composed of:
  1. computing \(G’(\omega)\) given \(S(\omega)\), and
  2. estimating \(M(\omega)\) using eq. (\ref{eq.inversion.obs.Gdash}) with given values of \(U(\omega)\), \(S(\omega)\), and \(G’(\omega)\).
The pow3-4 function was developed as a candidate for this \(S(\omega)\).
\(S(\omega)\)の条件として、まず\(G’(\omega)\)を安定に計算できる必要がある。 ここでは理論波形を速度・応力型の差分法を用いて計算することを想定している。 この場合、速度の遠地項が震源時間関数の2階導関数に比例する (Aki and Richards, 2002)ことから 2階導関数までが連続なものを震源時間関数として用いることが望ましい。 そのためpow3-4関数では2階導関数までが連続という条件を課している。
One requirement for \(S(\omega)\) is that \(G’(\omega)\) can be calculated stably. Here, a velocity-stress type finite difference method is assumed to be used for the calculation of the theoretical waveform. In this case, the second order derivative of the source time function should be continuous, because the far-field term of the velocity is proportional to the second order derivative of the source time function (Aki and Richards, 2002). This is the reason why the continuity to the second order derivative is demanded for function pow3-4.
\(S(\omega)\)のもう一つの条件として、スペクトルホールを持たないことが望ましい。 もしも\(S(\omega)\)がスペクトルホールを持つ場合、 \(G’(\omega)\)にもスペクトルホールが生じる(\ref{eq.Gdash}式)。 波形インバージョンは大まかに言えば周波数領域での割り算であり、 その分母に\(G’(\omega)\)が用いられる (\ref{eq.inversion.obs.Gdash}式)のであるから、 \(G’(\omega)\)がスペクトルホールを持つと インバージョンが不安定になると考えられる。
Another requirement for \(S(\omega)\) is that it should not have spectral holes. If \(S(\omega)\) has spectral holes, \(G’(\omega)\) also has spectral holes (eq. \ref{eq.Gdash}). Waveform inversion is roughly a division in frequency domain, where \(G’(\omega)\) is used as the denominator (eq. \ref{eq.inversion.obs.Gdash}), suggesting that spectral holes in \(G’(\omega)\) would result in instability in the inversion.
実はpow3-4関数の「\(t=\tau_p/2\)に関して非対称」という条件は スペクトルホールを持たないようにするためのものである。 対称な関数をフーリエ変換すると \(\exp(i\omega\tau_p)-\exp(-i\omega\tau_p)\)の組合せに由来する \(\sin\omega\tau_p\)が 振幅スペクトルに現れ、スペクトルホールが発生する。 非対称な関数であればこの組合せは現れないので スペクトルホールが現れないと期待される。 例えば非対称な関数 \[\begin{equation} f(t)= \begin{cases} 0 & (t<0, t>1) \\ t & (0 \leq t \leq 1) \end{cases} \label{eq.function_without_hole} \end{equation}\] はスペクトルホールを持たないが、\(t=1/2\)に関して対称な関数 \[\begin{equation} f(t)= \begin{cases} 0 & (t<0, t<1) \\ 2t & (0 \leq t \leq 1/2) \\ 2(1-t) & (1/2 < t \leq 1) \end{cases} \label{eq.function_with_hole} \end{equation}\] はスペクトルホールを持つ。 関数pow3-4においてはこのような理由から非対称であるという要件を設けた。
The requirement of pow3-4 function that it must be asymmetric with respect to \(t=\tau_p/2\) is related to avoiding spectral holes. The Fourier transformation of a symmetric function produces a combination of \(\exp(i\omega\tau_p)-\exp(-i\omega\tau_p)\), yielding \(\sin\omega\tau_p\) in the amplitude spectrum which causes the spectral holes. In case of an asymmetric function, this combination does not appear in the spectrum, so that the spectral holes are expected to be avoided. For example, an asymmetric function given by eq. (\ref{eq.function_without_hole}) does not have a spectral hole, while a symmetric function of eq. (\ref{eq.function_with_hole}) has spectral holes. The requirement of asymmetry for function pow3-4 is based on this reason.
この関数の作り方の概要は以下の通りである。 \[\begin{equation} f(t)= \begin{cases} 0 & (t<0, t>\tau_p) \\ At^m (\tau_p-t)^n & (0 \leq t \leq \tau_p) \end{cases} \label{eq.derive} \end{equation}\] (\(A\)は定数)の形の関数を考え、\(m \geq 3\), \(n \geq 3\)とすれば、 インパルス的な形かつ2階導関数までが連続となる。 但し\(m=n\)の場合には\(t=\tau_p/2\)に関して対称な形になって スペクトルホールが発生してしまう。 そこで\(m=3\), \(n=4\)とすることによってスペクトルホールの発生を回避し、 最大値が1になるように定数\(A\)を決めれば (\ref{eq.definition})式の関数形が得られる。
An outline of constructing this function is as follows. Consider a function of the form given by eq. (\ref{eq.derive}), where \(A\) is a constant. This function has an impulse-like shape and the second order derivative is continuous if \(m \geq 3\) and \(n \geq 3\). However, when \(m=n\), the function becomes symmetric with respect to \(t=\tau_p/2\) and thus has spectral holes. Avoiding the holes by using \(m=3\) and \(n=4\) and determining the constant \(A\) so that the maximum value of the function is 1, the function form of eq. (\ref{eq.definition}) is obtained.

◆形式(Format)

#include <sequence/timefunc.h>
inline struct sequence create_timefunc_pow34
(const double tp,const int size,const double t0,const double dt,
 const int Nintegral)


◆引数(Arguments)

tp (\ref{eq.definition})式における\(\tau_p\)の値。
The value of \(\tau_p\) in eq. (\ref{eq.definition}).
size 作成する時系列データのサンプル数。
The number of samples of the time series data to create.
t0 作成する時系列データの先頭時刻。
The beginning time of the time series data to create.
dt 作成する時系列データのサンプリング間隔。
The sampling interval of the time series data to create.
Nintegral 積分回数。以下の4つの値のいずれかを指定する。
The number of integrals, which must be one of the followings.
\(-1\) 時間関数の微分形(\ref{eq.derivative}式)を作成する。
Create the differential form (eq. \ref{eq.derivative}) of the time function.
0 時間関数そのもの(\ref{eq.definition}式)を作成する。
Create the time function itself (eq. \ref{eq.definition}).
1 時間関数の1階積分形(\ref{eq.integral1}式)を作成する。
Create the 1st order integral form (eq. \ref{eq.integral1}) of the time function.
2 時間関数の2階積分形(\ref{eq.integral2}式)を作成する。
Create the 2nd order integral form (eq. \ref{eq.integral2}) of the time function.


◆戻り値(Return value)

引数Nintegralの値に応じて (\ref{eq.definition}) (\ref{eq.derivative}) (\ref{eq.integral1}) (\ref{eq.integral2})の いずれかの時間関数を表す時系列データ。
A time series data which represents the time function of either of eqs. (\ref{eq.definition}), (\ref{eq.derivative}), (\ref{eq.integral1}), and (\ref{eq.integral2}) depending on argument Nintegral.

◆使用例(Example)

struct sequence timefunc=create_timefunc_pow34(5.0,10001,-2.0,0.01,1);

この例では(\ref{eq.integral1})式で\(\tau_p=5\) sとした時系列データが \(t\in [-2,98]\)の時間範囲において0.01 sのサンプリング間隔で作成される。
In this example, a time series data based on eq. (\ref{eq.integral1}) with \(\tau_p=5\) s is created in a time range \(t\in [-2,98]\) with 0.01 s sampling intervals.