関数rayParameters_1Dscalar.html マニュアル

(The documentation of function rayParameters_1Dscalar.html)
Last Update: 2021/12/6


◆機能・用途(Purpose)

与えられた1次元速度構造のもとで2つの地点を結ぶ波線パラメータを推定する (波線パラメータは複数になる場合がある)。
Investigate ray parameters connecting two points using a given 1D velocity structure (there may be multiple ray parameters).

◆形式(Format)

#include <ray/calculate1Dscalar.h>
inline int rayParameters_1Dscalar
(const double x1,const double z1,const double x2,const double z2,
 const struct _1DvelocityStructure structure, double ∗∗pList)


◆引数(Arguments)

x1 計算に用いる一方の地点の水平位置。
The horizontal location of one point used for the calculation.
z1 \(x_1\)に対応する地点の深さ。
The depth of the point corresponding to \(x_1\).
x2 計算に用いる他方の地点の水平位置。
The horizontal location of the other point used for the calculation.
z2 \(x_2\)に対応する地点の深さ。
The depth of the point corresponding to \(x_2\).
structure 1次元速度構造。
A 1D velocity structure.
pList 推定した波線パラメータのリストの代入先。 宣言しただけのdouble ∗型変数を&を付けて与える。 関数内でこの変数が配列となり(動的メモリが確保され)、 その\(n\)番目の配列要素には \((x_1,z_1)\)と\((x_2,z_2)\)を結ぶ \(n\)番目の波線パラメータが代入される。
A memory into which the estimated ray parameter list will be inserted; prepare a double ∗-type variable that is only declared, and give it to this function with ‘&’. In this function, an array will be prepared for this function (by allocating a dynamic memory), and the \(n\)th ray parameter connecting \((x_1,z_1)\) and \((x_2,z_2)\) will be inserted into the \(n\)th component of the array.


◆戻り値(Return value)

\((x_1,z_1)\)と\((x_2,z_2)\)を結ぶ波線パラメータの個数。 波線パラメータが1つも見つからない場合はエラーとなる。
The number of ray parameters connecting \((x_1,z_1)\) and \((x_2,z_2)\). If there is no ray parameter that connects the two points, the program finishes as an error.

◆使用例(Example)

struct _1DvelocityStructure structure;
double ∗p;
int Np=rayParameters_1Dscalar(1000.0,100.0,3000.0,200.0,structure,&p);


◆計算方法(Computational algorithm)

簡単のため\(z_1<z_2\)の場合を例に説明する。 \(z_2\)における速度を\(v_2\)として、 \(p_{max}=1/v_2\)が波線パラメータの取りうる最大値である。 すなわち\(0\leq p\leq p_{max}\)の範囲で波線パラメータを探索すれば良い。
  1. \(z_1\)から\(z_2\)の深さに向かって深くなる一方の波線
  2. \(z_2\)よりも深部まで一旦潜ってから戻ってくる波線
で推定方法が大きく異なるので両者を分けて説明する。
For simplicity, we consider a situation where \(z_1<z_2\). The maximum possible ray parameter is given by \(p_{max}=1/v_2\), where \(v_2\) is the velocity at \(z_2\). Therefore the ray parameter can be searched within a range \(0\leq p\leq p_{max}\). In the following, we describe the algorithms for
  1. the ray that increases its depth monotonically from \(z_1\) to \(z_2\), and
  2. the ray that goes deeper than \(z_2\) and returns back.

1. \(z_1\)から\(z_2\)の深さに向かって深くなる一方の波線 (The ray that increases its depth monochromatically from \(z_1\) to \(z_2\))

\(z_2\)の深さに到達するまでの水平距離が 波線パラメータの単調増加関数であることが 用いている計算式の(15)式より分かる。 したがって、適当に波線パラメータを与えて計算した水平距離が \(|x_1-x_2|\)よりも小さければ波線パラメータを大きく、 \(|x_1-x_2|\)よりも大きければ波線パラメータを小さくする、 ということを繰り返せば波線パラメータを推定できる。
The horizontal distance until the ray reaches the depth of \(z_2\) is a monotonic increasing function of the ray parameter according to eq. (15) of the formulas. Thus, if the horizontal distance calculated with a trial ray parameter is less than (greater than ) \(|x_1-x_2|\), then the true ray parameter is greater than (less than) the trial value. Thus ray parameter can be investigated by iteratively correcting the trial value.
実際にはこの計算は次の要領で行う。 まず\(p=p_{max}/2\)を初期値として水平距離を計算する。 その距離が\(|x_1-x_2|\)よりも小さければ\(p\)を\(p_{max}/4\)だけ大きく、 \(|x_1-x_2|\)よりも大きければ\(p\)を\(p_{max}/4\)だけ小さくして もう一度水平距離を計算する。 その結果を再び\(|x_1-x_2|\)と比較し、今度は\(p_{max}/8\)だけ動かす。 以下\(p_{max}/16\)、\(p_{max}/32\)、というように 動かす幅を1/2倍しながら計算していき、 計算した水平距離と実際の\(|x_1-x_2|\)との差が \(|x_1-x_2|\)の\(10^{-10}\)倍よりも小さくなったら 解が収束したものとして計算を終了する。 また、50回繰り返しても解が収束しない場合は 適切な波線パラメータが無いものとして探索を終了する。
The investigation of the ray parameter for this case is conducted in the following way. First, \(p=p_{max}/2\) is used is the trial value to compute the horizontal distance. If the calculated distance is less than (greater than) \(|x_1-x_2|\), then the trial value is increased (decreased) by \(p_{max}/4\), and the horizontal distance is computed again. In the next trial, the \(p\) value is increased/decreased by \(p_{max}/8\) based on the comparison of the true and calculated distances. In the same way, the width of the correction of \(p\) becomes half in the evely step, as \(p_{max}/16\), \(p_{max}/32\), …, and if the difference between the calculated and real horizontal distances is less than \(10^{-10}\) times \(|x_1-x_2|\), then the trial ray parameter at that moment is regarded as the true ray parameter. If this condition is not met after more than 50 cycles of the trials, then the search finishes without giving the estimate of the ray parameter.

2. \(z_2\)よりも深部まで一旦潜ってから戻ってくる波線 (The ray that goes deeper than \(z_2\) and returns back)

この場合、速度構造によって波線パラメータと水平距離の関係は複雑に変化しうる。 解の個数すらも明白ではない。 そこで、\(0\leq p\leq p_{max}\)の範囲をまず10000等分し、 それぞれの波線パラメータ\(p_1, p_2, \cdots, p_{10000}\) について水平距離を計算する。 隣り合う2つの波線パラメータ(例えば\(p_1\)と\(p_2\))の間で 水平距離の理論値と\(|x_1-x_2|\)との大小関係が異なる場合に その間のどこかに解があるものと仮定し、 かつこの間の区間(\(p_1\leq p\leq p_2\))では波線パラメータと距離の関係が 単調増加または単調減少になっているものと仮定して、 1の場合と同じ考え方で波線パラメータの推定を行う。
For this case, the relation between the ray parameter and horizontal distance can be complicated depending on the velocity structure. Even the number of solutions is not obvious. In the program, the horizontal distances are computed for 10000 homogeneously sampled ray parameters \(p_1, p_2, \cdots, p_{10000}\) in the range \(0\leq p\leq p_{max}\). If the calculated horizontal distances between two adjacent ray parameters (e.g., \(p_1\) and \(p_2\)) are in the opposite sides of the true distance, then a true ray parameter is assumed to be between them, and also a monotonic relation between the ray parameter and the distance is assumed within this range (i.e., \(p_1\leq p\leq p_2\)). Then the ray parameter can be investigated in the same way as the case 1.

◆補足(Additional remarks)



◆検証(Validation)

この関数によって波線パラメータを正しく推定できることは以下の方法で検証した。 まず波線パラメータを与えて\((x,z)=(0,0)\)を出発する波線の位置を計算する。 ここで\(x\)軸は水平面内での波線の伝搬方向、\(z\)軸は鉛直下向きに取るものとする。 次にそれらの位置を与えてこの関数により波線パラメータを求める。 波線パラメータは複数得られる場合もあるが、 その中に与えた波線パラメータは必ず含まれるはずであるので、 それが含まれているか否かによってこの関数の動作を検証する。
Whether this function gives correct estimates of the ray parameters was examined by the following way. First, the locations on the ray starting from \((x,z)=(0,0)\) for a given parameter are calculated, where \(x\) is taken horizontally to the ray propagation direction, and \(z\) is taken downward. Then, using these locations, the ray parameters are investigated by this function. Although there may be multiple ray parameters, the given ray parameter must be included in the list of the ray parameters. This function is thus checked by evaluating whether the given ray parameter is in the list of the ray parameters estimated.
この検証では下記の速度構造を使用する。 なお深さ15000 m以深での速度勾配は 深さ13000-15000 mの速度勾配と同じにしている。 速度勾配が複雑に変化する構造であり、 この構造のもとで解析解とこの関数を用いた結果が偶然一致するとは考えにくい。
The velocity structure shown below is used for this check. The velocity slope below the depth of 15000 m is assumed to be equal to that in 13000-15000 m. Owing to the complicated velocity slope changes, it is unlikely that the travel times computed by this function are “by chance” consistent with the analytical solution.
深さ [m]
Depth [m] 		速度 [m/s]
Velocity [m/s]
02000
10002200
30003200
50004000
60005000
80005100
110006900
130007000
150007500


使用する波線パラメータ (The ray parameter used)

この関数の検証では位置\((x,z)=(0,0)\)から単調に深くなる波線だけでなく、 一旦深部まで潜って戻ってくる波線についても検証を行う必要がある。 そこで深さ14000 mにおいて水平になる波線を考える。 深さ14000 mにおける速度は7250 m/sであるので、 この波線の波線パラメータは \(p=\sin 90^{\circ}/7250=1/7250\) [s/m] である。深さ0 mでの射出角は\(\sin^{-1}(2000/7250)\approx 16^{\circ}\)となる。
In the validation of this function, not only a ray that increases its depth monochromatically from \((x,z)=(0,0)\) but also a ray that goes deep and returns back need to be checked. To realize it, let as consider a ray that becomes horizontal at a depth of 14000 m. The velocity at that depth is 7250 m/s, giving a ray parameter \(p=\sin 90^{\circ}/7250=1/7250\) [s/m]. The take-off angle at a depth of 0 m is \(\sin^{-1}(2000/7250)\approx 16^{\circ}\).

波線の位置の解析解 (Analytical solutions for the locations of the ray)

波線パラメータが決まったので、波線が各深さに到達するまでの水平距離を 計算式の(15)式によって計算できる (計算式自体は関数ray_horizontal_distance_1Dscalarの検証において チェック済みである)。 ここでは深さを1000 m毎に区切って求める (差分計算ではなく積分結果の解析解を用いるので刻みは粗くて構わない)。
As the ray parameter was determined, the horizontal distance of the ray that reaches each depth can be calculated by eq. (15) of the formula; the formula was checked in the validation of function ray_horizontal_distance_1Dscalar. Here, the horizontal distances are calculated for a depth interval of 1000 m; as the analytical solution of an integrated form is used (not a differential equation form), the interval needs not be fine.
まず波線が深さ1000 mに達するまでの水平距離を求める。 これは(15)式に\(z_s=0\) [m], \(z_d=1000\) [m], \(v_s=2000\) [m/s], \(v_d=2200\) [m/s], \(p=1/7250\) [s/m] を代入することで\(X_{0m\rightarrow 1000m}=302.662843\) [m]と計算できる。
First, let us calculate the horizontal distance when the ray reaches a depth of 1000 m. This is calculated by inserting \(z_s=0\) [m], \(z_d=1000\) [m], \(v_s=2000\) [m/s], \(v_d=2200\) [m/s], and \(p=1/7250\) [s/m] into eq. (15), resulting in \(X_{0m\rightarrow 1000m}=302.662843\) [m].
次に波線が深さ1000 mに到達してから2000 mに到達するまでの水平距離を求める。 これは(15)式に\(z_s=1000\) [m], \(z_d=2000\) [m], \(v_s=2200\) [m/s], \(v_d=2700\) [m/s], \(p=1/7250\) [s/m] を代入することで\(X_{1000m\rightarrow 2000m}=359.326273\) [m]と計算できる。 したがって深さ2000 mにおける波線の水平位置は \(X_{0m\rightarrow 1000m}+X_{1000m\rightarrow 2000m}=661.989116\) [m] となる。
Next, let us consider the horizontal distance of the ray as it propagates from depths of 1000 m to 2000 m. This is calculated by inserting \(z_s=1000\) [m], \(z_d=2000\) [m], \(v_s=2200\) [m/s], \(v_d=2700\) [m/s], and \(p=1/7250\) [s/m] into eq. (15), resulting in \(X_{1000m\rightarrow 2000m}=359.326273\) [m]. Therefore the horizontal location of the ray at a depth of 2000 m is calculated as \(X_{0m\rightarrow 1000m}+X_{1000m\rightarrow 2000m}=661.989116\) [m].
以下同様にして深さ14000 mまでの波線の水平位置を計算できる。 深さ14000 mで波線が水平になり、そこからは浅い方に戻っていくので 今度は深さを減らしながら同様の計算を繰り返していけば 波線が地表に戻るまでの水平位置のリストが得られる。
In the same way, horizontal distances of the ray until it reaches a depth of 14000 m can be calculated. At the depth of 14000 m, the ray becomes horizontal, and then it returns back to shallower depths. For this part, the calculation needs to be done for decreasing depths. In this way, a list of hozirontal locations of the ray until it returns back to the ground surface is obtained.
このようにして求めた波線上の位置のリストを下表に示す。 これらはymaeda_opentoolsの関数は使用せず、 (15)式を用いてawkで計算したものである。 なお小数点以下6桁での計算結果を示すが、 実際の計算は丸め誤差回避のためより多い桁数で計算している。
The locations on the ray calculated in this way are shown in the table below. These values were obtained by awk with eq. (15), without using the functions of ymaeda_opentools. Although 6 digits are shown for the dicimal part, more digits are used in the calculation to avoid rounding errors.
深さ[m]
Depth [m] 		1000 mの深さ変化に伴う水平距離の増分[m]
An increment of the horizontal distance corresponding to a depth change of 1000 m [m] 		水平距離(累積)[m]
Horizontal distance (cumulative) [m]
1000 \(X_{0m\rightarrow 1000m}=302.662843\) 302.662843
2000 \(X_{1000m\rightarrow 2000m}=359.326273\) 661.989115
3000 \(X_{2000m\rightarrow 3000m}=445.819390\) 1107.808505
4000 \(X_{3000m\rightarrow 4000m}=531.307022\) 1639.115527
5000 \(X_{4000m\rightarrow 5000m}=615.896032\) 2255.011559
6000 \(X_{5000m\rightarrow 6000m}=796.693311\) 3051.704870
7000 \(X_{6000m\rightarrow 7000m}=961.545571\) 4013.250441
8000 \(X_{7000m\rightarrow 8000m}=980.218465\) 4993.468906
9000 \(X_{8000m\rightarrow 9000m}=1121.129049\) 6114.597955
10000 \(X_{9000m\rightarrow 10000m}=1487.346214\) 7601.944169
11000 \(X_{10000m\rightarrow 11000m}=2270.675553\) 9872.619722
12000 \(X_{11000m\rightarrow 12000m}=3228.891287\) 13101.511010
13000 \(X_{12000m\rightarrow 13000m}=3530.362705\) 16631.873715
14000 \(X_{13000m\rightarrow 14000m}=7549.834435\) 24181.708150
13000 \(X_{14000m\rightarrow 13000m}=7549.834435\) 31731.542585
12000 \(X_{13000m\rightarrow 12000m}=3530.362705\) 35261.905290
11000 \(X_{12000m\rightarrow 11000m}=3228.891287\) 38490.796577
10000 \(X_{11000m\rightarrow 10000m}=2270.675553\) 40761.472130
9000 \(X_{10000m\rightarrow 9000m}=1487.346214\) 42248.818345
8000 \(X_{9000m\rightarrow 8000m}=1121.129049\) 43369.947394
7000 \(X_{8000m\rightarrow 7000m}=980.218465\) 44350.165859
6000 \(X_{7000m\rightarrow 6000m}=961.545571\) 45311.711430
5000 \(X_{6000m\rightarrow 5000m}=796.693311\) 46108.404741
4000 \(X_{5000m\rightarrow 4000m}=615.896032\) 46724.300773
3000 \(X_{4000m\rightarrow 3000m}=531.307022\) 47255.607794
2000 \(X_{3000m\rightarrow 2000m}=445.819390\) 47701.427184
1000 \(X_{2000m\rightarrow 1000m}=359.326273\) 48060.753457
0 \(X_{1000m\rightarrow 0m}=302.662843\) 48363.416300


波線パラメータの推定 (Estimating the ray parameters)

次に、\((x,z)=(0,0)\)と上で得た波線上の位置を結ぶ波線パラメータを この関数を用いて推定する。 例えば位置\((x,z)=(0,0)\)と\((x,z)=(302.662843,1000)\)を結ぶ波線パラメータは 以下のコードによって推定できる。
Next, let us investigate the ray parameters that connect \((x,z)=(0,0)\) and the locations on the ray obtained above using this function. For example, the ray parameters that connect the locations \((x,z)=(0,0)\) and \((x,z)=(302.662843,1000)\) can be estimated by the code below.
int Np;
double ∗pList;
struct _1DvelocityStructure structure;
Np=rayParameters_1Dscalar (0.0,0.0,302.662843,1000.0,structure,&pList);

ここでは根幹部分のみを抜き出して記載した。 実際には関数呼び出し前に変数strctureに上記の速度構造を設定する必要があり、 また推定した波線パラメータを書き出さなければ意味がない。
Here, only the most essential part of the code is shown. Indeed, the velocity structure given above must be set for the variable structure, and the estimated ray parameters must be shown.
推定された波線パラメータは以下の通りである。 いずれも出発地点の位置は\((x,z)=(0,0)\)として 到達地点のみを動かして推定した。 なお出発地点と到達地点が完全に同じ深さになる\((x,z)=(48363.416300,0)\)および 波線が完全に水平になる\((x,z)=(24181.708150,14000)\)は除外した。 直感的に理解しやすいように波線パラメータそのものではなくその逆数を示した。 これはその波線がちょうど水平になる深さでの速度を表している。
The estimated ray parameters are shown below, where the origin was fixed at \((x,z)=(0,0)\) and the various destination locations shown in the table were used. However, the destination locations \((x,z)=(48363.416300,0)\), which has exactly the same depth as the origin, and \((x,z)=(24181.708150,14000)\), for which the ray is completely horizontal, were not used. To easily understand the values, the inverse of the estimated ray parameters are shown, which represent the velocity at the depth where the ray becomes horizontal.
波線の到達地点の\((x,z)\)座標[m]
The \((x,z)\) coordinate of the destination location of the ray [m] 		推定された波線パラメータの逆数[m/s]
The inverse of the estimated ray parameters [m/s]
(302.662843, 1000.0) 7249.999998
(661.989115, 2000.0) 7250.000004
(1107.808505, 3000.0) 7250.000001
(1639.115527, 4000.0) 7250.000000
(2255.011559, 5000.0) 7249.999999
(3051.704870, 6000.0) 7250.000000
(4013.250441, 7000.0) 7250.000000
(4993.468906, 8000.0) 7250.000000
(6114.597955, 9000.0) 7250.000000
(7601.944169, 10000.0) 7250.000000
(9872.619722, 11000.0) 7250.000000
(13101.511010, 12000.0) 7250.000000
(16631.873715, 13000.0) 7250.000000
(31731.542585, 13000.0) 7002.045022, 7250.000000, 7008.049413
(35261.905290, 12000.0) 7250.000000, 7090.912958, 6956.175570
(38490.796577, 11000.0) 7250.000000, 7178.141855, 6931.476229
(40761.472130, 10000.0) 7249.999998, 7221.445871, 6927.875722
(42248.818345, 9000.0) 7249.999996, 7237.401867, 6927.220636
(43369.947394, 8000.0) 7250.000033, 7246.376812, 6926.907547
(44350.165859, 7000.0) 7253.454800, 7250.000031, 6926.678607
(45311.711430, 6000.0) 7260.554564, 7249.999997, 6926.459368, 5100.045011, 5096.119175
(46108.404741, 5000.0) 7265.798440, 7250.000001, 6926.306465, 5100.355535, 5089.309499
(46724.300773, 4000.0) 7269.268220, 7249.999998, 6926.210081, 5100.552120, 5086.835515
(47255.607794, 3000.0) 7272.101704, 7249.999996, 6926.133433, 5100.723526, 5085.064553
(47701.427184, 2000.0) 7274.358576, 7249.999998, 6926.073693, 5100.863952, 5083.788325
(48060.753457, 1000.0) 7276.091641, 7250.000003, 6926.028608, 5100.972775, 5082.881320

与えたパラメータの逆数が7250 [m/s]であったので、 推定された波線パラメータの逆数の中に7250 [m/s]が含まれていれば良いが、 最大で0.000033 [m/s]のずれで全観測点に7250 [m/s]が含まれている。 したがって波線パラメータの推定は小数点以下4桁(有効数字8桁)の精度で 正しく行われたと判断できる。
As the inverse of the ray parameter given was 7250 [m/s], the estimated ray parameters should consist of 7250 [m/s], and indeed this value was included for all the destination locations with a maximum error of 0.000033 [m/s]. Therefore it is considered that the ray parameters were investigated correctly with an accuracy of 4 digits in the decimal part (8 significant figures).
なお深部に潜ってから戻ってくる波線パラメータの推定においては \(0\leq p\leq p_{max}\)を10000等分した値を初期値としている (上記「計算方法」参照)。 これを1000等分に減らした場合、 到達地点(43369.947393, 8000.0)および(44350.165858, 7000.0) において7250 [m/s]に対応する波線パラメータが抜け落ちてしまった。 今回の例題では2000等分にすれば抜け落ちは防げたが、 速度構造と位置の組み合わせによっては抜け落ちが生じる恐れがあると考えて 余裕を見て10000等分を採用することにした。
In the estimation of ray parameters for rays that go deep and then return back, a total of 10000 initial values are homogeneously sampled in the range \(0\leq p\leq p_{max}\) as was explained in the “Computational algorithm” section above. When this number was reduced to 1000, the ray parameter corresponding to 7250 [m/s] was dropped for the destination points (43369.947393, 8000.0) and (44350.165858, 7000.0). In the current example, it was not dropped by using 2000 initial values. However, the dropping may occur depending on the velocity structure and locations. Therefore, the 10000 initial values are used to make the dropping less frequent.
また「計算方法」で述べたように、 波線パラメータを使って計算した水平距離と実際の水平距離とのずれが 水平距離の\(10^{-10}\)倍よりも小さくなったら 正しい波線パラメータが得られたと判断する仕組みにしている。 これを\(10^{-8}\)倍にした場合、波線パラメータの逆数が 7250 [m/s]から最大で0.001682 [m/s]だけずれた。 このことから\(10^{-8}\)倍を閾値とした場合には有効数字6桁の精度になってしまい、 しばしば用いられる%e表記において最終桁までを信用できる精度 (有効数字7桁)が得られない。 そのため\(10^{-10}\)倍を閾値とすることにした。
Also, as was explained in the “Computation algorithm” section, a trial ray parameter is regarded to be a correct ray parameter if the difference between the true and calculated horizontal distances is less than \(10^{-10}\) times the horizontal distance. When this threshold was lowered to \(10^{-8}\) times, the inverse of the estimated ray parameter had the maximum error of 0.001682 [m/s] from 7250 [m/s]. This means that the accuracy is only 6 significant figures when the threshold of \(10^{-8}\) times is used, which is not accurate enough as a “%e” representation (7 significant figures) is frequently used. Therefore the threshold of \(10^{-10}\) times is adopted.