ray/calculate1Dscalar.hで用いている計算式

(Formula used in ray/calculate1Dscalar.h)

Last Update: 2021/12/6

1. 以下の説明で用いる記号 (Notation used in the description below)

z

軸を鉛直下向きに取り、深さ

z

における速度を

v (z)

、波線が

z

軸となす角を

θ (z)

とする。速度の定義点の深さを

z = z_{1}, \dots, z_{N}

、それらの深さでの速度を

v_{1}, \dots, v_{N}

とする。このとき

\begin{matrix} (1) & v_{n} = v (z_{n}), n = 1, \dots, N \end{matrix}

である。

z_{n} \leq z \leq z_{n + 1}

(

n = 1, \dots, N - 1

)における速度勾配を

\begin{matrix} (2) & s_{n} = \frac{v_{n + 1} - v_{n}}{z_{n + 1} - z_{n}} \end{matrix}

とする。また

z < z_{1}

における速度勾配を

s_{0} (> 0)

、

z > z_{N}

における速度勾配を

s_{N} (> 0)

とする。波線パラメータを

\begin{matrix} (3) & p = \frac{\sin θ (z)}{v (z)} \end{matrix}

とする。
We assume the

z

axis to lie downward, and

v (z)

and

θ (z)

be the velocity and the angle between the ray and the

z

axis, respectively. We assume the velocities

v_{1}, \dots, v_{N}

to be given at depths of

z = z_{1}, \dots, z_{N}

, suggesting eq. (

1

). Let

s_{n}

(

n = 1, \dots, N - 1

) be the velocity slope in

z_{n} \leq z \leq z_{n + 1}

given by eq. (

2

), and

s_{0} (> 0)

and

s_{N} (> 0)

be the velocity slopes in

z < z_{1}

and

z > z_{N}

, respectively. We denote the ray parameter by

p

, which satisfies eq. (

3

).

以上の記号とプログラムで用いている変数との関係は下表の通りである。
The table below represents the relation between the notation in this documentation and variables in the program code.

本マニュアルでの記号 Notation in this documentation	プログラム中の記号 Variables in the program code
$N$	structure.N
$z_{n}$ ( $n = 1, \dots, N$ )	structure.depths[n-1]
$v_{n}$ ( $n = 1, \dots, N$ )	structure.velocities[n-1]
$s_{0}$	structure.slope_shallow
$s_{N}$	structure.slope_deep
$p$	p

2. 波線の最大到達深度 (The maximum depth of the ray)

波線の最大到達深度は

\sin θ (z) = 1

となる

z

として定義される。 (

3

)式より、この条件は

\begin{matrix} (4) & p v (z) = 1 \end{matrix}

と等価である。

s_{0} > 0

s_{N} > 0

という条件を付けていることで

v (z)

は0から

\infty

までの全範囲の値を取るので、どのような波線パラメータについても (

4

)式を満たす

z

が最低1つは存在することが保証される。 (

4

)式の解が複数存在する場合はその中の一番小さい

z

を波線の最大到達深度と見なす。なぜならば、ひとたび角度が90°になったらそれ以上深くは潜れないからである。
The maximum depth of the ray is defined as the solution of

\sin θ (z) = 1

, which is equivalent to eq. (

4

) because of eq. (

3

). Since

s_{0} > 0

and

s_{N} > 0

v (z)

takes all values from 0 to

\infty

, which ensure that at least one solution of eq. (

4

) exists for an arbitrary ray parameter. When two or more solutions of eq. (

4

) exists, the smallest

z

among these solutions is adopted as the maximum depth because once the angle becomes 90°, the ray does not go deeper.

波線の最大到達深度は関数ray_max_depth_1Dscalarにおいて以下の3つのケースに場合分けして計算される。
The maximum depth of the ray is computed by function ray_max_depth_1Dscalar switching for the following three cases.

$p v_{1} > 1$ のとき、( $4$ )式の最小解は $z < z_{1}$ の範囲にある。この範囲の速度は $\begin{matrix} (5) & v (z) = v_{1} + s_{0} (z - z_{1}) \end{matrix}$ と表され、この関係を( $4$ )式に代入して $z$ について解くと $\begin{matrix} (6) & z = z_{1} + \frac{1}{s_{0}} (\frac{1}{p} - v_{1}) \end{matrix}$ が得られる。
When $p v_{1} > 1$ , the minimum solution of eq. ( $4$ ) lies in a region $z < z_{1}$ . The velocity in this region is represented by eq. ( $5$ ), and inserting this relation into eq. ( $4$ ) and solving for $z$ gives eq. ( $6$ ).
$p v_{n} > 1$ となる最初の $n$ が $m$ ( $2 \leq m \leq N$ )のとき、 ( $4$ )式の最小解は $z_{m - 1} < z < z_{m}$ の範囲にある。この範囲の速度は $\begin{matrix} (7) & v (z) = v_{m} + s_{m - 1} (z - z_{m}) \end{matrix}$ と表され、この関係を( $4$ )式に代入して $z$ について解くと $\begin{matrix} (8) & z = z_{m} + \frac{1}{s_{m}} (\frac{1}{p} - v_{m}) \end{matrix}$ が得られる。
When $p v_{n} > 1$ first occurs for $n = m$ , where $2 \leq m \leq N$ , the minimum solution of eq. ( $4$ ) lies in a region $z_{m - 1} < z < z_{m}$ . The velocity in this region is represented by eq. ( $7$ ), and inserting this relation into eq. ( $4$ ) and solving for $z$ gives eq. ( $8$ ).
いかなる $n$ についても $p v_{n} < 1$ のとき、 ( $4$ )式の最小解は $z > z_{N}$ の範囲にある。この範囲の速度は $\begin{matrix} (9) & v (z) = v_{N} + s_{N} (z - z_{N}) \end{matrix}$ と表され、この関係を( $4$ )式に代入して $z$ について解くと $\begin{matrix} (10) & z = z_{N} + \frac{1}{s_{N}} (\frac{1}{p} - v_{N}) \end{matrix}$ が得られる。
When $p v_{n} < 1$ for all $n$ , the minimum solution of eq. ( $4$ ) lies in a region $z > z_{N}$ . The velocity in this region is represented by eq. ( $7$ ), and inserting this relation into eq. ( $4$ ) and solving for $z$ gives eq. ( $8$ ).

3. 速度が深さの1次関数で表される2つの深さの間の波線に沿った水平距離・実距離・走時 (Horizontal distance, travel distance, and travel time along the ray from one depth to another between which the velocity is a linear function of depth)

2つの深さ

z = z_{s}

z = z_{d}

(

z_{s} < z_{d}

)があり、

z_{s} < z < z_{d}

の範囲では速度が深さの1次関数(速度勾配

s

)で与えられているとする。また、波線の最大到達深度は

z \geq z_{d}

の範囲にあるものとする。この場合に

z = z_{s}

を出発した波線が最初に

z = z_{d}

に到達するまでの波線に沿った水平距離および実距離と走時を求めることを考える。
Let

z = z_{s}

and

z = z_{d}

be two depths (

z_{s} < z_{d}

), between which the velocity is a linear function of depth with the slope

s

, We assume the maximum depth of the ray to be

z \geq z_{d}

. Then, let us calculate the horizontal and travel distances as well as the travel time along the ray as it travels from

z = z_{s}

and first reaches

z = z_{d}

3.1. 水平距離 (Horizontal distance)

一般に、深さ

z

において波線が微小距離

Δ z

だけ深さ方向に進むときの水平方向距離は

\begin{matrix} (11) & Δ x = Δ z \tan θ (z) = Δ z \frac{\sin θ (z)}{\sqrt{1 - \sin^{2} θ (z)}} \end{matrix}

と表される。この関係は(

3

)式を用いると

\begin{matrix} (12) & Δ x = Δ z \frac{p v (z)}{\sqrt{1 - p^{2} v (z)^{2}}} \end{matrix}

と変形できるので、これを

z = z_{s}

から

z = z_{d}

まで積分することにより、水平距離の計算式として

\begin{matrix} (13) & X_{z_{s} \to z_{d}} = \int_{z_{s}}^{z_{d}} \frac{p v (z)}{\sqrt{1 - p^{2} v (z)^{2}}} d z \end{matrix}

を得る。
In general, the horizontal distance along the ray as it moves a small vertical distance

Δ z

at a depth

z

is given by eq. (

11

). Inserting eq. (

3

) into (

11

) results in (

12

), and integrating this equation from

z = z_{s}

z = z_{d}

gives (

13

) as the formula for the horizontal distance.

この積分を実行するにあたり、

z_{s} < z < z_{d}

の範囲で

v (z)

が

z

の1次関数で書けることを利用する。

v (z_{s}) = v_{s}

v (z_{d}) = v_{d}

とおき、速度勾配が

d v / d z = s

であることを用いると

\begin{array}{rcl} X_{z_{s} \to z_{d}} & = & \int_{v_{s}}^{v_{d}} \frac{p v}{\sqrt{1 - p^{2} v^{2}}} \frac{d z}{d v} d v \\ = & - \frac{1}{s p} \int_{v_{s}}^{v_{d}} \frac{- p^{2} v}{\sqrt{1 - p^{2} v^{2}}} d v \\ = & - \frac{1}{s p} {[\sqrt{1 - p^{2} v^{2}}]}_{v_{s}}^{v_{d}} \\ = & - \frac{1}{s p} (\sqrt{1 - p^{2} v_{d}^{2}} - \sqrt{1 - p^{2} v_{s}^{2}}) \\ (14) & = & \frac{\sqrt{1 - p^{2} v_{s}^{2}} - \sqrt{1 - p^{2} v_{d}^{2}}}{s p} \end{array}

と積分を解析的に計算できる。
Now, remember that

v (z)

is a linear function of

z

in the

z_{s} < z < z_{d}

range. Using

v (z_{s}) = v_{s}

v (z_{d}) = v_{d}

, and the velocity slope

d v / d z = s

, the integral of eq. (

13

) can be analytically calculated as eq. (

14

).

(

14

)式は

v_{s} = v_{d}

のときに0/0の不定値になってしまう。これを避けるため、

s = (v_{d} - v_{s}) / (z_{d} - z_{s})

であることを用いて変形を行うと

\begin{array}{rcl} X_{z_{s} \to z_{d}} & = & \frac{\sqrt{1 - p^{2} v_{s}^{2}} - \sqrt{1 - p^{2} v_{d}^{2}}}{s p} \frac{\sqrt{1 - p^{2} v_{s}^{2}} + \sqrt{1 - p^{2} v_{d}^{2}}}{\sqrt{1 - p^{2} v_{s}^{2}} + \sqrt{1 - p^{2} v_{d}^{2}}} \\ = & \frac{(1 - p^{2} v_{s}^{2}) - (1 - p^{2} v_{d}^{2})}{s p (\sqrt{1 - p^{2} v_{s}^{2}} + \sqrt{1 - p^{2} v_{d}^{2}})} \\ = & \frac{p^{2} (v_{d}^{2} - v_{s}^{2})}{\frac{v_{d} - v_{s}}{z_{d} - z_{s}} p (\sqrt{1 - p^{2} v_{s}^{2}} + \sqrt{1 - p^{2} v_{d}^{2}})} \\ (15) & = & \frac{p (v_{d} + v_{s}) (z_{d} - z_{s})}{\sqrt{1 - p^{2} v_{s}^{2}} + \sqrt{1 - p^{2} v_{d}^{2}}} \end{array}

を得る。この式ならば

v_{s} = v_{d}

の場合も適用可能であり、 (

13

)式で

v (z) = v_{s} = v_{d}

(定数)とおいて直接得られる式と一致する。
Eq. (

14

) gives 0/0 when

v_{s} = v_{d}

. This problem is avoided by arranging the equation as (

15

), where

s = (v_{d} - v_{s}) / (z_{d} - z_{s})

was used. Eq. (

15

) is applicable to the case of

v_{s} = v_{d}

, and is consistent with the direct result from eq. (

13

) assuming

v (z) = v_{s} = v_{d}

(constant).

関数ray_horizontal_distance_1Dscalar_linearでは (

15

)式を用いて水平距離の計算を行う。
Function ray_horizontal_distance_1Dscalar_linear computes the horizontal distance using eq. (

15

3.2. 実距離 (Travel distance)

深さ

z

において波線が微小距離

Δ z

だけ深さ方向に進むとき、水平方向を考慮に入れた波線に沿った微小距離は

\begin{matrix} (16) & Δ l = \frac{Δ z}{\cos θ (z)} = \frac{Δ z}{\sqrt{1 - \sin^{2} θ (z)}} \end{matrix}

と表される。(

3

)式を代入すると

\begin{matrix} (17) & Δ l = \frac{Δ z}{\sqrt{1 - p^{2} v (z)^{2}}} \end{matrix}

となる。これを

z = z_{s}

から

z = z_{d}

まで積分すれば、

z = z_{s}

を出発した波線が最初に

z = z_{d}

に到達するまでの波線に沿った実距離は

\begin{matrix} (18) & L_{z_{s} \to z_{d}} = \int_{z_{s}}^{z_{d}} \frac{1}{\sqrt{1 - p^{2} v (z)^{2}}} d z \end{matrix}

と書ける。
The total distance along the ray as it moves a small vertical distance

Δ z

at a depth

z

is given by eq. (

16

). Inserting (

3

) into (

16

) results in (

17

). Integrating this equation from

z = z_{s}

z = z_{d}

yields (

18

) as the formula for the travel distance of the ray that propagates from the depth

z = z_{s}

z = z_{d}

.

この積分を実行するにあたり、

z_{s} < z < z_{d}

の範囲で

v (z)

が

z

の1次関数で書けることを利用する。

v (z_{s}) = v_{s}

v (z_{d}) = v_{d}

とおき、速度勾配が

d v / d z = s

であることを用いると

\begin{array}{rcl} L_{z_{s} \to z_{d}} & = & \int_{v_{s}}^{v_{d}} \frac{1}{\sqrt{1 - p^{2} v^{2}}} \frac{d z}{d v} d v \\ (19) & = & \frac{1}{s} \int_{v_{s}}^{v_{d}} \frac{1}{\sqrt{1 - p^{2} v^{2}}} d v \end{array}

と変形できる。更に

v_{s} = \sin θ_{s} / p

v_{d} = \sin θ_{d} / p

とおき、

v = \sin θ / p

の変数変換を施すと

\begin{array}{rcl} L_{z_{s} \to z_{d}} & = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} \frac{1}{\sqrt{1 - \sin^{2} θ}} \frac{d v}{d θ} d θ \\ = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} \frac{1}{\cos θ} \frac{\cos θ}{p} d θ \\ = & \frac{1}{s p} \int_{θ_{s}}^{θ_{d}} d θ \\ (20) & = & \frac{θ_{d} - θ_{s}}{s p} \end{array}

となり、積分を解析的に計算できる。
Now, remember that

v (z)

is a linear function of

z

in the

z_{s} < z < z_{d}

range. Using

v (z_{s}) = v_{s}

v (z_{d}) = v_{d}

, and the velocity slope

d v / d z = s

, the integral of eq. (

13

) can be arranged as (

19

). Further converting the integral variable from

v

θ

using a relation

v = \sin θ / p

, with defining

v_{s} = p \sin θ_{s}

and

v_{d} = p \sin θ_{d}

, the integral can be analytically calculated. The result is (

20

).

(

20

)式は

v_{s} = v_{d}

のときに0/0の不定値になってしまう。これを避けるため、

s = (v_{d} - v_{s}) / (z_{d} - z_{s})

であることを用いて変形を行うと

\begin{array}{rcl} L_{z_{s} \to z_{d}} & = & \frac{θ_{d} - θ_{s}}{\frac{v_{d} - v_{s}}{z_{d} - z_{s}} p} \\ = & \frac{(z_{d} - z_{s}) (θ_{d} - θ_{s})}{(v_{d} - v_{s}) p} \\ = & \frac{(z_{d} - z_{s}) (θ_{d} - θ_{s})}{v_{d} p - v_{s} p} \\ = & \frac{(z_{d} - z_{s}) (θ_{d} - θ_{s})}{\sin θ_{d} - \sin θ_{s}} \\ = & \frac{(z_{d} - z_{s}) (θ_{d} - θ_{s})}{2 \sin \frac{θ_{d} - θ_{s}}{2} \cos \frac{θ_{d} + θ_{s}}{2}} \\ (21) & = & \frac{z_{d} - z_{s}}{\cos \frac{θ_{d} + θ_{s}}{2}} \frac{\frac{θ_{d} - θ_{s}}{2}}{\sin \frac{θ_{d} - θ_{s}}{2}} \end{array}

となる。

v_{d} \to v_{s}

とした極限では

θ_{d} \to θ_{s}

となり、このとき

\begin{matrix} (22) & \frac{\frac{θ_{d} - θ_{s}}{2}}{\sin \frac{θ_{d} - θ_{s}}{2}} \to 1 \end{matrix}

より

L_{z_{s} \to z_{d}} \to (z_{d} - z_{s}) / \cos θ_{s}

となって波線を直線にした場合の式と一致する。
Eq. (

20

) gives 0/0 when

v_{s} = v_{d}

. This problem is avoided by arranging the equation as (

21

). At the limit of

v_{d} \to v_{s}

, asymptotic behaviours

θ_{d} \to θ_{s}

and (

22

) hold, thus giving

L_{z_{s} \to z_{d}} \to (z_{d} - z_{s}) / \cos θ_{s}

, which is consistent with the direct result assuming a straight ray.

関数ray_travel_distance_1Dscalar_linearでは (

21

)式を用いて波線に沿った実距離の計算を行う。
Function ray_travel_distance_1Dscalar_linear computes the travel distance along the ray using eq. (

21

3.3. 走時 (Travel time)

深さ

z

において波線が微小距離

Δ z

だけ深さ方向に進むとき、水平方向を考慮に入れた波線に沿った微小距離は(

17

)で表される。この微小距離を進むのにかかる時間は

\begin{matrix} (23) & Δ t = \frac{Δ l}{v (z)} = \frac{Δ z}{v (z) \sqrt{1 - p^{2} v (z)^{2}}} \end{matrix}

と書ける。これを

z = z_{s}

から

z = z_{d}

まで積分すれば、

z = z_{s}

を出発した波線が最初に

z = z_{d}

に到達するまでの所要時間は

\begin{matrix} (24) & T_{z_{s} \to z_{d}} = \int_{z_{s}}^{z_{d}} \frac{1}{v (z) \sqrt{1 - p^{2} v (z)^{2}}} d z \end{matrix}

と書ける。
The total distance along the ray as it moves a small vertical distance

Δ z

at a depth

z

is given by eq. (

17

). The time needed for the wave to propagate by this distance is given by eq. (

23

). Integrating this equation from

z = z_{s}

z = z_{d}

yields (

24

) as the formula for the travel time of the ray that propagates from the depth

z = z_{s}

z = z_{d}

.

この積分を実行するにあたり、

z_{s} < z < z_{d}

の範囲で

v (z)

が

z

の1次関数で書けることを利用する。

v (z_{s}) = v_{s}

v (z_{d}) = v_{d}

とおき、速度勾配が

d v / d z = s

であることを用いると

\begin{array}{rcl} T_{z_{s} \to z_{d}} & = & \int_{v_{s}}^{v_{d}} \frac{1}{v \sqrt{1 - p^{2} v^{2}}} \frac{d z}{d v} d v \\ (25) & = & \frac{1}{s} \int_{v_{s}}^{v_{d}} \frac{1}{v \sqrt{1 - p^{2} v^{2}}} d v \end{array}

と変形できる。更に

v_{s} = \sin θ_{s} / p

v_{d} = \sin θ_{d} / p

とおき、

v = \sin θ / p

の変数変換を施すと

\begin{array}{rcl} T_{z_{s} \to z_{d}} & = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} \frac{1}{\frac{\sin θ}{p} \sqrt{1 - \sin^{2} θ}} \frac{d v}{d θ} d θ \\ = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} \frac{p}{\sin θ \cos θ} \frac{\cos θ}{p} d θ \\ = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} \frac{d θ}{\sin θ} \\ = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} \frac{d θ}{2 \sin \frac{θ}{2} \cos \frac{θ}{2}} \\ = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} \frac{1}{2} (\frac{\cos \frac{θ}{2}}{\sin \frac{θ}{2}} + \frac{\sin \frac{θ}{2}}{\cos \frac{θ}{2}}) d θ \\ = & \frac{1}{s} \int_{θ_{s}}^{θ_{d}} (\frac{1}{2} \frac{\cos \frac{θ}{2}}{\sin \frac{θ}{2}} - \frac{1}{2} \frac{- \sin \frac{θ}{2}}{\cos \frac{θ}{2}}) d θ \\ = & \frac{1}{s} {[\ln (\sin \frac{θ}{2}) - \ln (\cos \frac{θ}{2})]}_{θ_{s}}^{θ_{d}} \\ = & \frac{1}{s} {[\ln \sqrt{\frac{1 - \cos θ}{2}} - \ln \sqrt{\frac{1 + \cos θ}{2}}]}_{θ_{s}}^{θ_{d}} \\ = & \frac{1}{s} {[\ln \sqrt{\frac{1 - \cos θ}{1 + \cos θ}}]}_{θ_{s}}^{θ_{d}} \\ = & \frac{1}{s} {[\ln \sqrt{\frac{(1 - \cos θ)^{2}}{(1 + \cos θ) (1 - \cos θ)}}]}_{θ_{s}}^{θ_{d}} \\ = & \frac{1}{s} {[\ln \sqrt{\frac{(1 - \cos θ)^{2}}{1 - \cos^{2} θ}}]}_{θ_{s}}^{θ_{d}} \\ = & \frac{1}{s} {[\ln \sqrt{\frac{(1 - \cos θ)^{2}}{\sin^{2} θ}}]}_{θ_{s}}^{θ_{d}} \\ = & \frac{1}{s} {[\ln (\frac{1 - \cos θ}{\sin θ})]}_{θ_{s}}^{θ_{d}} \\ = & \frac{1}{s} [\ln (\frac{1 - \cos θ_{d}}{\sin θ_{d}}) - \ln (\frac{1 - \cos θ_{d}}{\sin θ_{d}})] \\ = & \frac{1}{s} [\ln (\frac{1 - \sqrt{1 - p^{2} v_{d}^{2}}}{p v_{d}}) - \ln (\frac{1 - \sqrt{1 - p^{2} v_{s}^{2}}}{p v_{s}})] \\ (26) & = & \frac{1}{s} [\ln (\frac{1}{p v_{d}} - \sqrt{\frac{1}{p^{2} v_{d}^{2}} - 1}) - \ln (\frac{1}{p v_{s}} - \sqrt{\frac{1}{p^{2} v_{s}^{2}} - 1})] \end{array}

となり、積分を解析的に計算できる。
Now, remember that

v (z)

is a linear function of

z

in the

z_{s} < z < z_{d}

range. Using

v (z_{s}) = v_{s}

v (z_{d}) = v_{d}

, and the velocity slope

d v / d z = s

, the integral of eq. (

13

) can be arranged as (

25

). Further converting the integral variable from

v

θ

using a relation

v = \sin θ / p

, with defining

v_{s} = p \sin θ_{s}

and

v_{d} = p \sin θ_{d}

, the integral can be analytically calculated. The result is (

26

).

(

26

)式は

v_{d} = v_{s}

の場合に0/0の不定値になってしまう。また0/0が現れない形に変形するのも容易ではない。そこで、(

26

)式において

v_{d} \to v_{s}

とした極限での振る舞いを考える。そのために関数

\begin{matrix} (27) & f (v) \equiv \ln (\frac{1}{p v} - \sqrt{\frac{1}{p^{2} v^{2}} - 1}) \end{matrix}

を定義するのが便利である。この関数

f

の微分形は

\begin{array}{rcl} f^{'} (v) & = & \frac{1}{\frac{1}{p v} - \sqrt{\frac{1}{p^{2} v^{2}} - 1}} [- \frac{1}{p v^{2}} - \frac{1}{2 \sqrt{\frac{1}{p^{2} v^{2}} - 1}} (- \frac{2}{p^{2} v^{3}})] \\ = & \frac{1}{\frac{1}{p v} - \sqrt{\frac{1}{p^{2} v^{2}} - 1}} (- \frac{1}{p v^{2}} + \frac{1}{p^{2} v^{3} \sqrt{\frac{1}{p^{2} v^{2}} - 1}}) \\ = & \frac{1}{\frac{1}{p v} - \sqrt{\frac{1}{p^{2} v^{2}} - 1}} \frac{1}{p v^{2} \sqrt{\frac{1}{p^{2} v^{2}} - 1}} (- \sqrt{\frac{1}{p^{2} v^{2}} - 1} + \frac{1}{p v}) \\ = & \frac{1}{p v^{2} \sqrt{\frac{1}{p^{2} v^{2}} - 1}} \\ (28) & = & \frac{1}{v \sqrt{1 - p^{2} v^{2}}} \end{array}

\begin{array}{rcl} f^{″} (v) & = & (- \frac{1}{v^{2}}) \frac{1}{\sqrt{1 - p^{2} v^{2}}} + \frac{1}{v} [- \frac{1}{2} {(1 - p^{2} v^{2})}^{- 3 / 2}] (- 2 p^{2} v) \\ = & - \frac{1}{v^{2} \sqrt{1 - p^{2} v^{2}}} + \frac{p^{2}}{{(1 - p^{2} v^{2})}^{3 / 2}} \\ = & \frac{- (1 - p^{2} v^{2}) + p^{2} v^{2}}{v^{2} {(1 - p^{2} v^{2})}^{3 / 2}} \\ (29) & = & \frac{2 p^{2} v^{2} - 1}{v^{2} {(1 - p^{2} v^{2})}^{3 / 2}} \end{array}

\begin{array}{rcl} f^{(3)} (v) & = & \frac{4 p^{2} v \cdot v^{2} {(1 - p^{2} v^{2})}^{3 / 2} - (2 p^{2} v^{2} - 1)}{v^{4} {(1 - p^{2} v^{2})}^{3}} \\ [2 v {(1 - p^{2} v^{2})}^{3 / 2} + v^{2} \cdot \frac{3}{2} {(1 - p^{2} v^{2})}^{1 / 2} (- 2 p^{2} v)] \\ = & \frac{4 p^{2} v^{3} {(1 - p^{2} v^{2})}^{3 / 2} - (2 p^{2} v^{2} - 1)}{v^{4} {(1 - p^{2} v^{2})}^{3}} \\ [2 v {(1 - p^{2} v^{2})}^{3 / 2} - 3 p^{2} v^{3} {(1 - p^{2} v^{2})}^{1 / 2}] \\ = & \frac{1}{v^{4} {(1 - p^{2} v^{2})}^{3}} [4 p^{2} v^{3} {(1 - p^{2} v^{2})}^{3 / 2} - 2 v (2 p^{2} v^{2} - 1) {(1 - p^{2} v^{2})}^{3 / 2} \\ + 3 p^{2} v^{3} (2 p^{2} v^{2} - 1) {(1 - p^{2} v^{2})}^{1 / 2}] \\ = & \frac{4 p^{2} v^{2} (1 - p^{2} v^{2}) - 2 (2 p^{2} v^{2} - 1) (1 - p^{2} v^{2}) + 3 p^{2} v^{2} (2 p^{2} v^{2} - 1)}{v^{3} {(1 - p^{2} v^{2})}^{5 / 2}} \\ = & \frac{1}{v^{3} {(1 - p^{2} v^{2})}^{5 / 2}} [(4 p^{2} v^{2} - 4 p^{4} v^{4}) - (4 p^{2} v^{2} - 4 p^{4} v^{4} - 2 + 2 p^{2} v^{2}) \\ + (6 p^{4} v^{4} - 3 p^{2} v^{2})] \\ (30) & = & \frac{2 - 5 p^{2} v^{2} + 6 p^{4} v^{4}}{v^{3} {(1 - p^{2} v^{2})}^{5 / 2}} \end{array}

と書ける。関数

f

を用いると(

26

)式は

\begin{array}{rcl} T_{z_{s} \to z_{d}} & = & \frac{1}{s} [f (v_{d}) - f (v_{s})] \\ = & \frac{z_{d} - z_{s}}{v_{d} - v_{s}} [f (v_{d}) - f (v_{s})] \\ (31) & = & (z_{d} - z_{s}) \frac{f (v_{d}) - f (v_{s})}{v_{d} - v_{s}} \end{array}

と書けるので、

v_{d} \to v_{s}

の極限では

\begin{matrix} (32) & T_{z_{s} \to z_{d}} = (z_{d} - z_{s}) f^{'} (v_{s}) = \frac{z_{d} - z_{s}}{v_{s} \sqrt{1 - p^{2} v_{s}^{2}}} \end{matrix}

となる。これは(

24

)式において

v (z)

を

z

によらない定数と見なして積分した形になっている。ここでは速度が深さの線形関数で書ける場合を考えているのであるから

v_{d} \to v_{s}

の極限では速度は定数となる。したがって(

26

)式で

v_{d} \to v_{s}

とした極限 (

32

)が (

24

)式で

v (z)

を定数と見なした場合の式と一致するのは妥当な結果である。
Eq. (

26

) yields 0/0 when

v_{d} = v_{s}

, and arranging this equation to avoid 0/0 is difficult. Thus, let us consider an approximate form of eq. (

26

) at the limit of

v_{d} \to v_{s}

. A function

f (v)

defined by eq. (

27

) is convenient for this purpose. Differential forms of this function are given by eqs. (

28

), (

29

), and (

30

). Using this function, the travel time is given by (

31

), yielding (

32

) at the limit of

v_{d} \to v_{s}

. This equation is equivalent to the expression obtained by (

24

) assuming that

v (z)

is not dependent on

z

. Remember that our formulation here assumed a linear velocity change with depth. Thus the limit

v_{d} \to v_{s}

means a constant velocity. It is therefore reasonable that the limit of eq. (

26

) at

v_{d} \to v_{s}

, given by (

32

), is equivalent to the expression obtained by assuming a constant velocity in eq. (

24

).

ところでプログラミングでは一般に2つの実数の等号比較は不安定の原因であり、

v_{d} = v_{s}

のときとそうでないときで場合分けするのは良いコードとは言えない。ある幅を持たせて

v_{d}

と

v_{s}

が十分に近い場合とそうでない場合で場合分けする方がより安定な計算ができる。そのためには

v_{d}

が

v_{s}

と「完全には一致しないが十分近い」場合に成り立つ近似式が必要になる。そこで、

f (v_{d})

を

v_{s}

のまわりでテイラー展開して3次の項まで残すと

\begin{matrix} (33) & f (v_{d}) \sim f (v_{s}) + f^{'} (v_{s}) (v_{d} - v_{s}) + \frac{f^{″} (v_{s})}{2} (v_{d} - v_{s})^{2} + \frac{f^{(3)} (v_{s})}{6} (v_{d} - v_{s})^{3} \end{matrix}

と書ける。このとき(

31

)式より

\begin{array}{rcl} T_{z_{s} \to z_{d}} & \sim & (z_{d} - z_{s}) \frac{f^{'} (v_{s}) (v_{d} - v_{s}) + \frac{f^{″} (v_{s})}{2} (v_{d} - v_{s})^{2} + \frac{f^{(3)} (v_{s})}{6} (v_{d} - v_{s})^{3}}{v_{d} - v_{s}} \\ = & f^{'} (v_{s}) (z_{d} - z_{s}) + \frac{f^{″} (v_{s})}{2} (z_{d} - z_{s}) (v_{d} - v_{s}) \\ (34) & + \frac{f^{(3)} (v_{s})}{6} (z_{d} - z_{s}) (v_{d} - v_{s})^{2} \end{array}

となり、(

28

)-(

30

)を代入して

\begin{array}{rcl} T_{z_{s} \to z_{d}} & \sim & \frac{z_{d} - z_{s}}{v_{s} \sqrt{1 - p^{2} v_{s}^{2}}} + \frac{(2 p^{2} v_{s}^{2} - 1) (z_{d} - z_{s})}{2 v_{s}^{2} {(1 - p^{2} v_{s}^{2})}^{3 / 2}} (v_{d} - v_{s}) \\ (35) & + \frac{(2 - 5 p^{2} v_{s}^{2} + 6 p^{4} v_{s}^{4}) (z_{d} - z_{s})}{6 v_{s}^{3} {(1 - p^{2} v_{s}^{2})}^{5 / 2}} (v_{d} - v_{s})^{2} \end{array}

を得る。この最後の項を無視できる条件は

\begin{matrix} (36) & | \frac{(2 p^{2} v_{s}^{2} - 1) (z_{d} - z_{s})}{2 v_{s}^{2} {(1 - p^{2} v_{s}^{2})}^{3 / 2}} (v_{d} - v_{s}) | ≫ | \frac{(2 - 5 p^{2} v_{s}^{2} + 6 p^{4} v_{s}^{4}) (z_{d} - z_{s})}{6 v_{s}^{3} {(1 - p^{2} v_{s}^{2})}^{5 / 2}} (v_{d} - v_{s})^{2} | \end{matrix}

であり、共通因子を消去すると

\begin{matrix} (37) & | 2 p^{2} v_{s}^{2} - 1 | ≫ | \frac{2 - 5 p^{2} v_{s}^{2} + 6 p^{4} v_{s}^{4}}{3 v_{s} (1 - p^{2} v_{s}^{2})} (v_{d} - v_{s}) | \end{matrix}

整理して

\begin{array}{rcl} | \frac{v_{d} - v_{s}}{3 v_{s}} | & ≪ & | \frac{(1 - p^{2} v_{s}^{2}) (2 p^{2} v_{s}^{2} - 1)}{2 - 5 p^{2} v_{s}^{2} + 6 p^{4} v_{s}^{4}} | \\ = & | \frac{- 1 + 3 p^{2} v_{s}^{2} - 2 p^{4} v_{s}^{4}}{2 - 5 p^{2} v_{s}^{2} + 6 p^{4} v_{s}^{4}} | \\ (38) & = & | \frac{1 - 3 p^{2} v_{s}^{2} + 2 p^{4} v_{s}^{4}}{2 - 5 p^{2} v_{s}^{2} + 6 p^{4} v_{s}^{4}} | \end{array}

を得る。この条件が成り立つとき、(

35

)式の最後の項を無視した近似式

\begin{matrix} (39) & T_{z_{s} \to z_{d}} \sim \frac{z_{d} - z_{s}}{v_{s} \sqrt{1 - p^{2} v_{s}^{2}}} + \frac{(2 p^{2} v_{s}^{2} - 1) (z_{d} - z_{s})}{2 v_{s}^{2} {(1 - p^{2} v_{s}^{2})}^{3 / 2}} (v_{d} - v_{s}) \end{matrix}

を利用できる。
In the computer programs, two real numbers should not be compared by an equal sign because it would result in a numerical instability. Thus branching based on whether

v_{d} = v_{s}

or not is not a good code. Rather, the branching should be based on whether

v_{d}

and

v_{s}

are sufficiently close or not. To implement this branching, an approximated equation for the travel time that is valid when

v_{d}

and

v_{s}

are ``sufficiently close to each other but not exactly same'' is needed. We now derive this equation. A Taylor expansion of

f (v_{d})

with respect to

v_{s}

, remaining until the 3rd-order term, results in eq. (

33

). Using this relation, eq. (

31

) reduces to (

34

), and inserting eqs. (

28

)-(

30

) into it gives (

35

). The final term of this equation can be neglected when (

36

) holds. Removing the common factors from the both hand sides of (

36

) results in (

37

), which can further be arranged as (

38

). When this condition is met, it is valid to neglect the final term of (

35

), resulting in (

39

).

(

26

)式と(

39

)式の比較の例を図1に示す。ここでは(

26

)式に基づいて計算した走時を

T_{z_{s} \to z_{d}}^{e x a c t}

、 (

39

)式に基づいて計算した走時を

T_{z_{s} \to z_{d}}^{a p p r o x}

、 (

38

)式の左辺を右辺で割った値を

M_{a p p r o x}

とする。理論上は

M_{a p p r o x} ≪ 1

であれば近似が成り立ち、

M_{a p p r o x}

の値が小さいほど理論上の近似精度が良い。また

T_{z_{s} \to z_{d}}^{e x a c t}

と

T_{z_{s} \to z_{d}}^{a p p r o x}

が近いほど実際の近似精度が良い。これらを比較することで、実際に

M_{a p p r o x}

がどの程度小さければ

T_{z_{s} \to z_{d}}^{a p p r o x}

がどれくらい良い近似になるのかを評価できる。図1cは

M_{a p p r o x}

と近似式の誤差

| T_{z_{s} \to z_{d}}^{a p p r o x} - T_{z_{s} \to z_{d}}^{e x a c t} | / T_{z_{s} \to z_{d}}^{e x a c t}

との関係を示したものであり、

M_{a p p r o x} = 0.01

であれば近似式の誤差は

10^{- 4}

のオーダーになることが分かる。この誤差レベルは波線パラメータ

p

や層厚

z_{d} - z_{s}

にも依存すると思われるので、プログラムでは安全側に余裕を見て

M_{a p p r o x} < 0.0001

のときに近似式を用いるようにしている。
Fig. 1 is a comparison of the travel times computed by eqs. (

26

) and (

39

). Here,

T_{z_{s} \to z_{d}}^{e x a c t}

and

T_{z_{s} \to z_{d}}^{a p p r o x}

are the travel times based on eqs. (

26

) and (

39

), respectively, and

M_{a p p r o x}

is the ratio of left- to right-hand-side of eq. (

38

). The approximation is theoretically valid when

M_{a p p r o x} ≪ 1

; theoretically, the smaller the

M_{a p p r o x}

, the better the approximation. The real precision of the approximation is evaluated by comparing

T_{z_{s} \to z_{d}}^{e x a c t}

and

T_{z_{s} \to z_{d}}^{a p p r o x}

. In Fig. 1c, a relation between the theoretical goodness of the approximation (

M_{a p p r o x}

) and the real goodness of the approximation (

| T_{z_{s} \to z_{d}}^{a p p r o x} - T_{z_{s} \to z_{d}}^{e x a c t} | / T_{z_{s} \to z_{d}}^{e x a c t}

) is shown, indicating that a value of

M_{a p p r o x} = 0.01

results in a travel time error in order of

10^{- 4}

by the approximation. This error level may depend on the ray parameter

p

and the layer thickness

z_{d} - z_{s}

. To have a margin for the safety side, the program switches the exact and approximate equations at

M_{a p p r o x} = 0.0001

図1.

z_{s} = 0

z_{d} = 1000

p = 0.0001

s m

^{- 1}

v_{s} = 2000

m s

^{- 1}

のときの(

26

)式と(

39

)式の比較。

v_{d}

を1800 m s

^{- 1}

から2200 m s

^{- 1}

まで 0.1 m s

^{- 1}

刻みで動かして計算した。 (a)

v_{d}

と

T_{z_{s} \to z_{d}}

の関係。赤線は(

26

)式による厳密解

T_{z_{s} \to z_{d}}^{e x a c t}

、青線は(

39

)式による近似解

T_{z_{s} \to z_{d}}^{a p p r o x}

。 (b)

v_{d}

と理論上の近似精度

M_{a p p r o x}

との関係。この値が小さいほど理論上の近似精度が良い。 (c)

M_{a p p r o x}

と

| T_{z_{s} \to z_{d}}^{a p p r o x} - T_{z_{s} \to z_{d}}^{e x a c t} | / T_{z_{s} \to z_{d}}^{e x a c t}

との関係。
Fig. 1. A comparison of travel times computed by eqs. (

26

) and (

39

) for a case of

z_{s} = 0

z_{d} = 1000

p = 0.0001

s m

^{- 1}

v_{s} = 2000

m s

^{- 1}

, and moving

v_{d}

from 1800 m s

^{- 1}

to 2200 m s

^{- 1}

at an interval of 0.1 m s

^{- 1}

. (a) A relation between

v_{d}

and

T_{z_{s} \to z_{d}}

. The red and blue lines represent the exact solution

T_{z_{s} \to z_{d}}^{e x a c t}

from eq. (

26

) and an approximate form

T_{z_{s} \to z_{d}}^{a p p r o x}

based on eq. (

39

), respectively. (b) A relation between

v_{d}

and

M_{a p p r o x}

. The value of

M_{a p p r o x}

represents a theoretical goodness of the approximation; smaller

M_{a p p r o x}

values correspond to better approximations. (c) A relation between

M_{a p p r o x}

and

| T_{z_{s} \to z_{d}}^{a p p r o x} - T_{z_{s} \to z_{d}}^{e x a c t} | / T_{z_{s} \to z_{d}}^{e x a c t}

4. 速度定義点をまたぐ2つの深さの間の波線に沿った水平距離 (Horizontal distance along the ray between two depths over a velocity definition point)

前節では

z_{s} < z < z_{d}

において速度が深さの1次関数で表される場合を考えた。これは

z_{s} < z < z_{d}

の範囲に速度定義点が存在しない状況に対応する。

z_{s} < z < z_{d}

の範囲に速度定義点が含まれる場合、その地点で速度勾配が変わる。このようなケースでは

z_{s} < z < z_{d}

の範囲を速度定義点で小区間に分割し、それらの小区間ごとに(

15

)式を適用すれば良い。
In the previous section, the velocity was assumed to be a linear function in the

z_{s} < z < z_{d}

range, which corresponds to a situation where there is no velocity definition point in this range. If the velocity definition points are present in

z_{s} < z < z_{d}

, the velocity slope changes at these points. In this case, eq. (

15

) is applicable to subsections bounded by the velocity definition points.

いま、

z_{s}

よりも深部側の

z_{s}

に最も近い速度定義点を

z = z_{l}

、

z_{d}

よりも浅部側の

z_{d}

に最も近い速度定義点を

z = z_{m}

とする。

z_{s} < z < z_{d}

の間に少なくとも1つの速度定義点が含まれる状況を考えれば

l \geq m

となる。このとき、水平距離は

\begin{matrix} (40) & X_{z_{s} \to z_{d}} = X_{z_{s} \to z_{l}} + \sum_{n = l}^{m - 1} X_{z_{n} \to z_{n + 1}} + X_{z_{m} \to z_{d}} \end{matrix}

となる。この式の個々の項(

X_{z_{s} \to z_{l}}

など)は (

15

)式によって計算できる。
Let

z = z_{l}

and

z = z_{m}

be the nearest velocity definition points downward from

z_{s}

and upward from

z_{d}

, respectively. When at least one velocity definition point exists in

z_{s} < z < z_{d}

l \geq m

holds. In this case, the horizontal distance is calculated by eq. (

40

). Individual terms in this equation (

X_{z_{s} \to z_{l}}

, etc.) can be calculated with eq. (

15

).

関数ray_horizontal_distance_1Dscalarでは (

40

)式と(

15

)式を組み合わせて水平距離の計算を行う。
Function ray_horizontal_distance_1Dscalar computes the horizontal distance using eqs (

40

) and (

15

ray/calculate1Dscalar.hで用いている計算式

1. 以下の説明で用いる記号 (Notation used in the description below)

2. 波線の最大到達深度 (The maximum depth of the ray)

3. 速度が深さの1次関数で表される2つの深さの間の 波線に沿った水平距離・実距離・走時 (Horizontal distance, travel distance, and travel time along the ray from one depth to another between which the velocity is a linear function of depth)

3.1. 水平距離 (Horizontal distance)

3.2. 実距離 (Travel distance)

3.3. 走時 (Travel time)

4. 速度定義点をまたぐ2つの深さの間の波線に沿った水平距離 (Horizontal distance along the ray between two depths over a velocity definition point)

3. 速度が深さの1次関数で表される2つの深さの間の波線に沿った水平距離・実距離・走時 (Horizontal distance, travel distance, and travel time along the ray from one depth to another between which the velocity is a linear function of depth)