関数calculate_traveltime_multi_layers マニュアル

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


◆機能・用途(Purpose)

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

◆形式(Format)

#include <ray/calculateLayered.h>
inline double calculate_traveltime_multi_layers
(const double x1,const double z1,const double x2,const double z2,
 const int Nlayers,const double ∗v,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\).
Nlayers 層の数。
The number of layers.
v 各層の速度を浅い側から順に並べた配列。 配列要素数はNlayersであって 0<v[0]<v[1]<……v[Nlayers-1] でなければならない。
An array composed of the velocities in individual layers, ordered from shallow to deep. The number of array components must be Nlayers, and the values must satisfy 0<v[0]<v[1]<……v[Nlayers-1].
zb 層境界の標高を浅い側から順に並べた配列。 要素数はNlayers−1であって zb[0]>zb[1]>…>zb[Nlayers-2] でなければならない。
An array composed of the elevations of the layer boundaries, ordered from shallow to deep. The number of array components must be Nlayers−1, and the values must satisfy zb[0]>zb[1]>…>zb[Nlayers-2].


◆戻り値(Return value)

\((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)\).

◆使用例(Example)

double v[]={2000.0,3000.0,5000.0,6000.0};
double zb[]={1500.0,-1000.0,-5000.0};
double T =calculate_traveltime_multi_layers(0.0,-500.0,2000.0,1700.0,4,v,zb);


◆計算方法(Computation method)

2地点\((x_1,z_1)\), \((x_2,z_2)\)が同じ標高(\(z_1=z_2\))の場合、 2地点間の距離\(|x_1-x_2|\)を2地点が属する層の速度で割れば 走時が得られる。
In cases where the two points \((x_1,z_1)\) and \((x_2,z_2)\) are at the save elevation (i.e., \(z_1=z_2\)), the traveltime is obtained by dividing the distance between the two points \(|x_1-x_2|\) by the velocity of the layer which the two points belong to.
2地点の標高が等しくない(\(z_1\neq z_2\))場合、 2地点のうち標高が高い方を\((x_h,z_h)\)、低い方を\((x_l,z_l)\)とおく。 \((x_h,z_h)\)が属する層の速度を\(v_h\)、 \((x_l,z_l)\)が属する層の速度を\(v_l\)とおく。 2地点を結ぶ波線が\((x_h,z_h)\)において鉛直となす角を\(\theta_h\)、 \((x_l,z_l)\)において鉛直となす角を\(\theta_l\)とする。 これらの間には\(\sin\theta_h/v_h=\sin\theta_l/v_l\)の関係があり、 \(\theta_h\)と\(\theta_l\)のうちの片方が決まればもう片方も決まる。 この関数では\(\theta_l\)を独立なパラメータとして探索する。 2つの標高\(z_h\)と\(z_l\)を結ぶ波線の水平距離\(r\)は 仮定する\(\theta_l\)の単調増加関数であるので、 適当に\(\theta_l\)を仮定して計算してみた\(r\)が 真の水平距離\(|x_h-x_s|\)よりも大きければ 仮定した\(\theta_l\)が大きすぎたことになるので\(\theta_l\)を減らし、 逆に適当に\(\theta_l\)を仮定して計算してみた\(r\)が 真の水平距離\(|x_h-x_s|\)よりも小さければ 仮定した\(\theta_l\)が小さすぎたことになるので\(\theta_l\)を増やせば良い。 これを繰り返すことで2地点を結ぶ\(\theta_l\)を得ることができる。 この関数では\(0^{\circ}\leq\theta_l<90^{\circ}\)を探索範囲とし、 上の考え方のもとで2分法により\(\theta_l\)を求める。 ひとたび\(\theta_l\)が求まればスネルの法則により走時を計算できる。
In cases where the elevations of the two points are not same (i.e., \(z_1\neq z_2\)), let \((x_h,z_h)\) and \((x_l,z_l)\) be the higher and lower ones, respectively, of the two points. Let \(v_h\) and \(v_l\) be the velocities of the layers which \((x_h,z_h)\) and \((x_l,z_l)\), respectively, belong two. Let \(\theta_h\) and \(\theta_l\) be the angles (from the vertical) of the ray connecting the two points at \((x_h,z_h)\) and \((x_l,z_l)\), respectively. They are related by \(\sin\theta_h/v_h=\sin\theta_l=v_l\), meaning that \(\theta_h\) and \(\theta_l\) are calculated from one another. In this function, \(\theta_l\) is searched as an independent parameter. Because the horizontal distance \(r\) along the ray between the two elevations \(z_h\) and \(z_l\) monochromatically increases with the assumed \(\theta_l\), the optimal \(\theta_l\) can be investigated iteratively by reducing (increasing) it if \(r\) calculated from an assumed \(\theta_l\) was greater (less) than the actual horizontal distance \(|x_h-x_s|\). In this function, \(\theta_l\) is investigated by the bisection method with a search range \(0^{\circ}\leq\theta_l<90^{\circ}\). Once \(\theta_l\) was obtained, the travel time is calculated with the Snell's law.
以上は直達波の走時であり、実際にはhead waveも考慮しなければならない。 \(z_l\)よりも下の各層境界を通るhead waveの走時を計算し、 直達波とそれらのhead waveの走時の中の最小値を採用する。
The travel times calculated above are for the direct wave. Indeed, there also may be a head wave. Thus the travel time for the head wave that passes on each layer interface below \(z_l\) is computed, and the smallest travel time between the direct and head waves is finally adopted.

◆検証(Validation)

この関数が正しい答えを出すことのチェックのため、以下の構造を使用した。 十分に複雑な構造のため、偶然にして正解が得られることはまず無いと思われる。
To check that this function returns a correct answer, the structure below is used. This structure is complex enough to avoid the correct answer to be obtained by chance.
層
Layer 		標高範囲(m)
Elevation range (m) 		速度(m/s)
Velocity (m/s)
1 \(\geq\)20002500
2 500-20003500
3 0-5004000
4 \(-\)400-05500
5 \(\leq -\)4006000

この構造のもとで を考える。 各層におけるこれらの波線と鉛直とのなす角は下表のようになる。
Let us consider the rays The angles of these rays (from the vertical) in individual layers are given by the table below.
層
Layer 		点\(P_0\)からの波線の角度(°)
Angle of the ray from \(P_0\) (°) 		点\(Q_0\)からの波線の角度(°)
Angle of the ray from \(Q_0\) (°)
1 5.0000 38.508
2 7.0086 60.653
3 8.0160 85.000
4 11.054 -
5 12.074 -

また、各層の上端を通る屈折波(head wave)の波線と鉛直とのなす角は下表のようになる。
Angles of the ray (from the vertical) of the refracted (head) wave that passes through the top of each layer are given by the table below.
層
Layer 		第2層上端を通る屈折波の波線の角度(°)
Angle of the ray that passes through the top of the 2nd layer (°) 		第3層上端を通る屈折波の波線の角度(°)
Angle of the ray that passes through the top of the 3rd layer (°) 		第4層上端を通る屈折波の波線の角度(°)
Angle of the ray that passes through the top of the 4th layer (°) 		第5層上端を通る屈折波の波線の角度(°)
Angle of the ray that passes through the top of the 5th layer (°)
1 45.585 38.682 27.036 24.624
2 90.000 61.045 39.521 35.685
3 - 90.000 46.658 41.810
4 - - 90.000 66.444
5 - - - 90.000

これらの角度を用いて様々な標高における 点\(P_0\)から射出された波線上の点の座標と\(P_0\)からの走時の正解を計算すると 下表のようになる。 なお直達波の波線がどの屈折波(head wave)よりも高角度であるので 屈折波は存在しない。
The table below shows the coordinates and correct travel times at various elevations on the ray from \(P_0\), calculated from the angles given above. As the angle of the ray of the direct wave is steeper than those of all the refracted (head) waves, there is no refracted wave.
点
Point (m) 		直達波の走時(s)
Travel time of the direct wave (s)
\(P_1(126.2,2200.0)\) 0.1205
\(P_2(143.7,2000.0)\) 0.2008
\(P_3(266.7,1000.0)\) 0.4886
\(P_4(328.1,500.0)\) 0.6326
\(P_5(370.4,200.0)\) 0.7083
\(P_6(398.6,0.0)\) 0.7588
\(P_7(437.6,-200.0)\) 0.7958
\(P_8(476.7,-400.0)\) 0.8329
\(P_9(519.5,-600.0)\) 0.8670

一方、様々な標高における点\(Q_0\)から射出された波線上の点の座標と \(Q_0\)からの走時の正解を計算すると下表のようになる。 太字が走時の最小値を表す。
The table below shows the coordinates and correct travel times at various elevations on the ray from \(Q_0\). Here, bold represents the minimum value of the travel times.
点(m)
Point (m) 		直達波の走時(s)
Travel time of the direct wave (s) 		第4層上端を通る屈折波の走時(s)
Travel time of the refracted wave that passes through the top of the 4th layer (s) 		第5層上端を通る屈折波の走時(s)
Travel time of the refracted wave that passes through the top of the 5th layer (s)
\(Q_1(1343.0,300.0)\) 0.2868 0.2936 -
\(Q_2(2486.0,400.0)\) 0.5737 0.5186 -
\(Q_3(3057.5,450.0)\) 0.7171 0.6311 0.6555
\(Q_4(3629.0,500.0)\) 0.8605 0.7436 0.7601
\(Q_5(5407.6,1500.0)\) 1.4435 1.2873 1.2886
\(Q_6(6119.0,1900.0)\) 1.6767 1.5048 1.5000
\(Q_7(6296.9,2000.0)\) 1.7350 1.5592 1.5528
\(Q_8(6694.7,2500.0)\) 1.9906 1.8097 1.8009

これらの「走時の正解」と、この関数を用いて計算した走時との比較を下表に示す。 関数を用いた計算では出発地点と到達地点の位置および速度構造のみを与えており 射出角等は未知の状態で計算を行っている。 下表より正しい値が得られたことが分かる。
The table below indicates comparisons of these correct travel times and the travel times calculated with this function. Note that only the origin and destination locations as well as the velocity structure are given to the function, with the radiation angles unknown. Nevertheless, correct travel times are obtained, as the table below shows.
出発地点(m)
Origin (m) 		到達地点(m)
Destination (m) 		走時の正解(s)
Correct travel time (s) 		関数を用いて計算した走時(s)
Travel time calculated with this function (s)
\(P_0(100,2500)\) \(P_1(126.2,2200.0)\) 0.1205 0.1205
\(P_0(100,2500)\) \(P_2(143.7,2000.0)\) 0.2008 0.2008
\(P_0(100,2500)\) \(P_3(266.7,1000.0)\) 0.4886 0.4886
\(P_0(100,2500)\) \(P_4(328.1,500.0)\) 0.6326 0.6326
\(P_0(100,2500)\) \(P_5(370.4,200.0)\) 0.7083 0.7083
\(P_0(100,2500)\) \(P_6(398.6,0.0)\) 0.7588 0.7588
\(P_0(100,2500)\) \(P_7(437.6,-200.0)\) 0.7958 0.7958
\(P_0(100,2500)\) \(P_8(476.7,-400.0)\) 0.8329 0.8329
\(P_0(100,2500)\) \(P_9(519.5,-600.0)\) 0.8670 0.8670
\(Q_0(200,200)\) \(Q_1(1343.0,300.0)\) 0.2868 0.2868
\(Q_0(200,200)\) \(Q_2(2486.0,400.0)\) 0.5186 0.5186
\(Q_0(200,200)\) \(Q_3(3057.5,450.0)\) 0.6311 0.6311
\(Q_0(200,200)\) \(Q_4(3629.0,500.0)\) 0.7436 0.7436
\(Q_0(200,200)\) \(Q_5(5407.6,1500.0)\) 1.2873 1.2873
\(Q_0(200,200)\) \(Q_6(6119.0,1900.0)\) 1.5000 1.5000
\(Q_0(200,200)\) \(Q_7(6296.9,2000.0)\) 1.5528 1.5528
\(Q_0(200,200)\) \(Q_8(6694.7,2500.0)\) 1.8009 1.8009