関数calculate_traveltime_2layers マニュアル

(The documentation of function calculate_traveltime_2layers)
Last Update: 2021/12/6


◆機能・用途(Purpose)

水平2層構造のもとで2地点を結ぶ走時を計算する。
Compute the travel time between two points in a medium of two horizontal layers.

◆形式(Format)

#include <ray/calculateLayered.h>
inline double calculate_traveltime_2layers
(const double x1,const double z1,const double x2,const double z2,
 const double v1,const double v2,const double zb)


◆引数(Arguments)

x1 走時を計算する一方の地点の水平位置(\(x_1\))。
The horizontal location (\(x_1\)) of one point for the travel time calculation.
z1 \(x_1\)に対応する地点の標高(\(z_1\))。
The elevation (\(z_1\)) of the point corresponding to \(x_1\).
x2 走時を計算する他方の地点の水平位置(\(x_2\))。
The horizontal location (\(x_2\)) of the other point for the travel time calculation.
z2 \(x_2\)に対応する地点の標高(\(z_2\))。
The elevation (\(z_2\)) of the point corresponding to \(x_2\).
v1 第1層の速度(\(v_1\))。
The velocity (\(v_1\)) of the 1st layer.
v2 第2層の速度(\(v_2\))。
The velocity (\(v_2\)) of the 2nd layer.
zb 層境界の標高(\(z_b\))。
The elevation (\(z_b\)) of the layer boundary.


◆戻り値(Return value)

標高\(z>z_b\)で速度\(v_1\)、\(z<z_b\)で速度\(v_2\) の2層から成る媒質を介して \((x_1,z_1)\)から\((x_2,z_2)\)まで波が伝搬するのに要する時間。
The travel time of the wave propagating from \((x_1,z_1)\) to \((x_2,z_2)\) through a 2-layer medium for which the velocity is \(v_1\) and \(v_2\) for elevation ranges \(z>z_b\) and \(z<z_b\), respectively.

◆使用例(Example)

double T =calculate_traveltime_2layers(0.0,0.0,2000.0,0.0,2000.0,5000.0,-500.0);


◆用いている計算式(Formula used)

1. \(z_1 \geq z_b\)かつ\(z_2 \geq z_b\)の場合 (In cases where \(z_1 \geq z_b\) and \(z_2 \geq z_b\))

直達波の走時は \[\begin{equation} T_d=\frac{\sqrt{(x_2-x_1)^2+(z_2-z_1)^2}}{v_1} \label{eq.Tdirect.shallow} \end{equation}\] と書ける。
The travel time of the direct wave is written by eq. (\ref{eq.Tdirect.shallow}).
第2層の上端を通る屈折波(head wave)の走時は \[\begin{equation} T_r=\frac{l_{1b}}{v_1}+\frac{l_h}{v_2}+\frac{l_{b2}}{v_1} \label{eq.Trefracted.original} \end{equation}\] と書ける。ここで\(l_{1b}\)は 屈折波が\((x_1,z_1)\)を出発してから層境界に到達するまでの波線に沿った距離で、 層境界への入射角を\(\theta\)として \[\begin{equation} l_{1b}=\frac{z_1-z_b}{\cos\theta} \label{eq.l1b} \end{equation}\] と書ける。\(l_h\)は屈折波が第2層の上端を通る区間の長さで \[\begin{equation} l_h=|x_2-x_1|-(z_1-z_b)\tan\theta-(z_2-z_b)\tan\theta \label{eq.lh} \end{equation}\] と書ける。\(l_{b2}\)は 屈折波が層境界を離れてから\((x_2,z_2)\)に到達するまでの波線に沿った距離で \[\begin{equation} l_{b2}=\frac{z_2-z_b}{\cos\theta} \label{eq.lb2} \end{equation}\] と書ける。(\ref{eq.l1b})-(\ref{eq.lb2})を (\ref{eq.Trefracted.original})に代入して整理すると \[\begin{eqnarray} T_r &=& \frac{z_1-z_b}{v_1\cos\theta} +\frac{|x_2-x_1|-(z_1-z_b)\tan\theta-(z_2-z_b)\tan\theta}{v_2} +\frac{z_2-z_b}{v_1\cos\theta} \nonumber \\ &=& \frac{|x_2-x_1|}{v_2} +\frac{z_1-z_b}{v_1\cos\theta}\left(1-\frac{v_1}{v_2}\sin\theta\right) +\frac{z_2-z_b}{v_1\cos\theta}\left(1-\frac{v_1}{v_2}\sin\theta\right) \nonumber \\ &=& \frac{|x_2-x_1|}{v_2} +\frac{(z_1-z_b)+(z_2-z_b)}{v_1\cos\theta} \left(1-\frac{v_1}{v_2}\sin\theta\right) \label{eq.Trefracted.arrange1} \end{eqnarray}\] となる。スネルの法則より \[\begin{equation} \sin\theta=\frac{v_1}{v_2} \label{eq.snell.refracted} \end{equation}\] であるのでこの関係を用いて \[\begin{eqnarray} T_r &=& \frac{|x_2-x_1|}{v_2} +\frac{(z_1-z_b)+(z_2-z_b)}{v_1\cos\theta}\left(1-\sin^2\theta\right) \nonumber \\ &=& \frac{|x_2-x_1|}{v_2} +\frac{(z_1-z_b)+(z_2-z_b)}{v_1}\cos\theta \nonumber \\ &=& \frac{|x_2-x_1|}{v_2} +\frac{(z_1-z_b)+(z_2-z_b)}{v_1}\sqrt{1-\left(\frac{v_1}{v_2}\right)^2} \nonumber \\ &=& \frac{|x_2-x_1|}{v_2} +[(z_1-z_b)+(z_2-z_b)]\sqrt{\frac{1}{v_1^2}-\frac{1}{v_2^2}} \label{eq.Trefracted.final} \end{eqnarray}\] を得る。
The refracted (head) wave passing through the head of the 2nd layer is given by eq. (\ref{eq.Trefracted.original}). Here, \(l_{1b}\) is the distance along the ray of the refracted wave from \((x_1,z_1)\) to the boundary, given by eq. (\ref{eq.l1b}), where \(\theta\) is the incident angle of the ray to the boundary; \(l_h\) is the section length of the ray on the head of the 2nd layer, given by (\ref{eq.lh}); and \(l_{b2}\) is the distance along the ray from deparing the boundary to reach \((x_2,z_2)\), given by (\ref{eq.lb2}). Inserting eqs. (\ref{eq.l1b})-(\ref{eq.lb2}) into (\ref{eq.Trefracted.original}) reduces to (\ref{eq.Trefracted.arrange1}). Using the Snell's law given by eq. (\ref{eq.snell.refracted}), we finally obtain (\ref{eq.Trefracted.final}).
\(T_d\)と\(T_r\)の両方を計算し、小さい方が走時として採択される。
Both \(T_d\) and \(T_r\) are computed, and the smaller one is adopted as the travel time.

2. \(z_1<z_b\)かつ\(z_2<z_b\)の場合 (In cases where \(z_1<z_b\) and \(z_2<z_b\))

このときは直達波のみとなり、 \begin{equation} T_d=\frac{\sqrt{(x_2-x_1)^2+(z_2-z_1)^2}}{v_2} \label{eq.Tdirect.deep} \end{equation} と計算できる。
In this case, only the direct wave exists, and the travel time is computed by eq. (\ref{eq.Tdirect.deep}).

3. \(z_1 \geq z_b\)かつ\(z_2<z_b\)の場合 (In cases where \(z_1 \geq z_b\) and \(z_2<z_b)\)

波が層1から層2に入る地点の座標を\((x_b,z_b)\)とおくと、走時は\(x_b\)の関数として \[\begin{equation} T(x_b) =\frac{\sqrt{(x_b-x_1)^2+(z_1-z_b)^2}}{v_1} +\frac{\sqrt{(x_2-x_b)^2+(z_b-z_2)^2}}{v_2} \label{eq.Txb} \end{equation}\] と書ける。この走時を最小にするように\(x_b\)を求めることを考える。 これには\(dT(x_b)/dx_b=0\)とすれば良い。 \[\begin{eqnarray} \frac{dT(x_b)}{dx_b} &=& \frac{2(x_b-x_1)}{2v_1\sqrt{(x_b-x_1)^2+(z_1-z_b)^2}} +\frac{-2(x_2-x_b)}{2v_2\sqrt{(x_2-x_b)^2+(z_b-z_2)^2}} \nonumber \\ &=& \frac{x_b-x_1}{v_1\sqrt{(x_b-x_1)^2+(z_1-z_b)^2}} -\frac{x_2-x_b}{v_2\sqrt{(x_2-x_b)^2+(z_b-z_2)^2}} \label{eq.Txb.derivative} \end{eqnarray}\] であるので、\(dT(x_b)/dx_b=0\)はスネルの法則そのものである。 (\ref{eq.Txb.derivative})を更にもう1回微分すると \[\begin{eqnarray} & & \frac{d^2T(x_b)}{dx_b^2} \nonumber \\ &=& \frac{1}{v_1}\frac{d}{dx_b} \left\{(x_b-x_1)\left[(x_b-x_1)^2+(z_1-z_b)^2\right]^{-1/2}\right\} \nonumber \\ & & -\frac{1}{v_2}\frac{d}{dx_b} \left\{(x_2-x_b)\left[(x_2-x_b)^2+(z_b-z_2)^2\right]^{-1/2}\right\} \nonumber \\ &=& \frac{1}{v_1} \left\{ \left[(x_b-x_1)^2+(z_1-z_b)^2\right]^{-1/2} \right. \nonumber \\ & & \left. -\frac{1}{2}(x_b-x_1)\left[(x_b-x_1)^2+(z_1-z_b)^2\right]^{-3/2} \left[2(x_b-x_1)\right] \right\} \nonumber \\ & & -\frac{1}{v_2} \left\{ -\left[(x_2-x_b)^2+(z_b-z_2)^2\right]^{-1/2} \right. \nonumber \\ & & \left. -\frac{1}{2}(x_2-x_b)\left[(x_2-x_b)^2+(z_b-z_2)^2\right]^{-3/2} \left[-2(x_2-x_b)\right] \right\} \nonumber \\ &=& \frac{1}{v_1} \left\{ \frac{1}{\sqrt{(x_b-x_1)^2+(z_1-z_b)^2}} -\frac{(x_b-x_1)^2}{\left[\sqrt{(x_b-x_1)^2+(z_1-z_b)^2}\right]^3} \right\} \nonumber \\ & & +\frac{1}{v_2} \left\{ \frac{1}{\sqrt{(x_2-x_b)^2+(z_b-z_2)^2}} -\frac{(x_2-x_b)^2}{\left[\sqrt{(x_2-x_b)^2+(z_b-z_2)^2}\right]^3} \right\} \nonumber \\ &=& \frac{1}{v_1} \frac{\left[(x_b-x_1)^2+(z_1-z_b)^2\right]-(x_b-x_1)^2} {\left[\sqrt{(x_b-x_1)^2+(z_1-z_b)^2}\right]^3} \nonumber \\ & & +\frac{1}{v_2} \frac{\left[(x_2-x_b)^2+(z_b-z_2)^2\right]-(x_2-x_b)^2} {\left[\sqrt{(x_2-x_b)^2+(z_b-z_2)^2}\right]^3} \nonumber \\ &=& \frac{1}{v_1} \frac{(z_1-z_b)^2}{\left[\sqrt{(x_b-x_1)^2+(z_1-z_b)^2}\right]^3} +\frac{1}{v_2} \frac{(z_b-z_2)^2}{\left[\sqrt{(x_2-x_b)^2+(z_b-z_2)^2}\right]^3} \label{eq.Txb.derivative2} \end{eqnarray}\] となる。これより\(d^2T(x_b)/dx_b^2\geq 0\)であることが分かる。 したがって\(dT(x_b)/dx_b\)は\(x_b\)の単調増加関数となるので \(dT(x_b)/dx_b=0\)の解は2分法によって簡単に求められる。
Let \((x_b,z_b)\) be the point where the ray intersects the layer boundary. The travel time is expressed as a function of \(x_b\) in the form of eq. (\ref{eq.Txb}). The optimal \(x_b\) is obtained by minimizing this travel time. Differentiating the travel time with respect to \(x_b\) yields eq. (\ref{eq.Txb.derivative}), and the optimal \(x_b\) value is obtained from \(dT(x_b)/dx_b=0\), which is exactly same as the Snell's law. An additional differentiation for eq. (\ref{eq.Txb.derivative}) gives (\ref{eq.Txb.derivative2}), indicating that \(d^2T(x_b)/dx_b^2\geq 0\). Thus \(dT(x_b)/dx_b\) is a monochromatic increasing function of \(x_b\). Then the solution of \(dT(x_b)/dx_b=0\) can be investigated easily by the bisection method.

4. \(z_1<z_b\)かつ\(z_2 \geq z_b\)の場合 (In cases where \(z_1<z_b\) and \(z_2 \geq z_b\))

\((x_1,z_1)\)と\((x_2,z_2)\)を入れ替えて前節と同様の計算を行えば良い。
The same calculation as the previous section is conducted, with replacing \((x_1,z_1)\) and \((x_2,z_2)\) with each other.

◆検証(Validation)

\(z_1\geq z_b\)かつ\(z_2\geq z_b\)の場合の計算コードについては 以下の条件でチェックした。
The computation code for the cases of \(z_1\geq z_b\) and \(z_2\geq z_b\) was examined by the following condition:
\(z_1=z_2=0\), \(x_2>x_1=0\)であるので (\ref{eq.Tdirect.shallow})(\ref{eq.Trefracted.final})式はそれぞれ \[\begin{equation} T_d=\frac{x_2}{v_1} \label{eq.Tdirect.shallow.examine1} \end{equation}\] \[\begin{equation} T_r=\frac{x_2}{v_2}-2z_b\sqrt{\frac{1}{v_1^2}-\frac{1}{v_2^2}} \label{eq.Trefracted.examine1} \end{equation}\] となる。ごく初歩的な2層構造中の走時曲線の問題であり、 \(x_2\)が小さい場所では傾き\(1/v_1=1/300\)、 \(x_2\)が大きい場所では傾き\(1/v_2=1/500\)となる。 また折れ曲がりの場所は\(T_d=T_r\)となる条件 \[\begin{equation} \frac{x_2}{v_1}=\frac{x_2}{v_2}-2z_b\sqrt{\frac{1}{v_1^2}-\frac{1}{v_2^2}} \label{eq.bend_condition.examine1} \end{equation}\] から \[\begin{equation} x_2=-2z_b \frac{\sqrt{\frac{1}{v_1^2}-\frac{1}{v_2^2}}}{\frac{1}{v_1}-\frac{1}{v_2}} =-2z_b\frac{\sqrt{v_2^2-v_1^2}}{v_2-v_1} =-2z_b\sqrt{\frac{v_2+v_1}{v_2-v_1}} =2400 \label{eq.x2.examine1} \end{equation}\] となる。実際に上記の条件でこの関数を用いて計算し、 これらの理論通りの結果になることを確認できた。
Since \(z_1=z_2=0\) and \(x_2>x_1=0\), eqs. (\ref{eq.Tdirect.shallow}) and (\ref{eq.Trefracted.final}) would be simplified as (\ref{eq.Tdirect.shallow.examine1}) and (\ref{eq.Trefracted.examine1}), respectively. This is a very primary problem of the travel time curve in a 2-layer medium. The slopes are \(1/v_1=1/300\) and \(1/v_2=1/500\) for small and large \(x_2\), respectively. The bending position is given by the condition of \(T_d=T_r\), or explicitly, eq. (\ref{eq.bend_condition.examine1}), the solution of which is eq. (\ref{eq.x2.examine1}). A compute program using this function with the parameters given above yielded the results consistent with the theory.
次に\(z_1\geq z_b\)かつ\(z_2<z_b\)の場合について 以下のようにコードのチェックを行った。 \((x_1,z_1)=(0,0)\)から\(\theta_1=30^{\circ}\)の射出角で波を放射する場合を考える。 \(v_1=1000\), \(v_2=1000\sqrt{3}\), \(z_b=-500\), \(z_2=-1500\)とし、 波線が\(z=z_2\)に到達する場所の\(x\)座標を\(x_2\)とする。 この\(x_2\)の値および\((x_1,z_1)\)から\((x_2,z_2)\)までの走時は 次のように計算できる。 波線パラメータは\(p=\sin\theta_1/v_1=1/2000\)である。 第2層における射出角を\(\theta_2\)とおくと \(\sin\theta_2=pv_2=\sqrt{3}/2\)であるので \(\theta_2=60^{\circ}\)であることが分かる。 これより \[\begin{eqnarray} x_2 &=& (z_1-z_b)\tan\theta_1+(z_b-z_2)\tan\theta_2 \nonumber \\ &=& 500\tan 30^{\circ}+1000\tan 60^{\circ} \nonumber \\ &=& \frac{500}{\sqrt{3}}+1000\sqrt{3} \nonumber \\ &=& \frac{500+3\times 1000}{\sqrt{3}} \nonumber \\ &=& \frac{3500}{\sqrt{3}} \label{eq.x2.examine2} \end{eqnarray}\] と計算できる。また、このときの走時は \[\begin{eqnarray} T &=& \frac{z_1-z_b}{v_1\cos\theta_1}+\frac{z_2-z_b}{v_2\cos\theta_2} \nonumber \\ &=& \frac{z_1-z_b}{\frac{\sin\theta_1}{p}\cos\theta_1} +\frac{z_2-z_b}{\frac{\sin\theta_2}{p}\cos\theta_2} \nonumber \\ &=& \frac{p(z_1-z_b)}{\sin\theta_1\cos\theta_1} +\frac{p(z_2-z_b)}{\sin\theta_2\cos\theta_2} \nonumber \\ &=& \frac{\frac{1}{2000}\times 500}{\sin 30^{\circ}\cos 30^{\circ}} +\frac{\frac{1}{2000}\times 1000}{\sin 60^{\circ}\cos 60^{\circ}} \nonumber \\ &=& \frac{\frac{1}{4}}{\frac{1}{2}\frac{\sqrt{3}}{2}} +\frac{\frac{2}{4}}{\frac{\sqrt{3}}{2}\frac{1}{2}} \nonumber \\ &=& \frac{\frac{3}{4}}{\frac{\sqrt{3}}{4}} \nonumber \\ &=& \sqrt{3} \label{eq.T.examine2} \end{eqnarray}\] と計算できる。 上の計算では射出角が既知として\(x_2\)を計算で求めたのであるが、 この結果を用いれば\(x_2=3500/\sqrt{3}\)のときに \(T=\sqrt{3}\)となるべきであることが分かる。 実際にプログラムに\(x_1=z_1=0\), \(z_b=-500\), \(z_2=-1500\), \(v_1=1000\), \(v_2=1000\sqrt{3}\), \(x_2=3500/\sqrt{3}\)を与えて この関数を用いて走時を計算し、\(\sqrt{3}\)となることを確認した。
The program code for the cases of \(z_1\geq z_b\) and \(z_2<z_b\) was examined by the following way. Let the wave be emitted from \((x_1,z_1)=(0,0)\) with a takeoff angle of \(\theta_1=30^{\circ}\). Let the other parameters be \(v_1=1000\), \(v_2=1000\sqrt{3}\), \(z_b=-500\), and \(z_2=-1500\). The value \(x_2\) is given as the \(x\)-coordinate where the ray approaches \(z=z_2\). This \(x_2\) value and the travel time from \((x_1,z_1)\) to \((x_2,z_2)\) are theoretically calculated as follows. The ray parameter is \(p=\sin\theta_1/v_1=1/2000\). Let \(\theta_2\) be the takeoff angle in the 2nd layer. This angle satisfies \(\sin\theta_2=pv_2=\sqrt{3}/2\), indicating that \(\theta_2=60^{\circ}\). We therefore have eq. (\ref{eq.x2.examine2}) for \(x_2\) and eq. (\ref{eq.T.examine2}) for the travel time. In this example, \(x_2\) was calculated for a given takeoff angle. However, the results indicate that the travel time should be \(T=\sqrt{3}\) when \(x_2=3500/\sqrt{3}\). A computer program using this function, with \(x_1=z_1=0\), \(z_b=-500\), \(z_2=-1500\), \(v_1=1000\), \(v_2=1000\sqrt{3}\), and \(x_2=3500/\sqrt{3}\), indeed yielded the travel time of \(\sqrt{3}\).