Алгоритм решения прямой кинематической задачи сейсмики в трехмерных неоднородных изотропных средах Текст научной статьи по специальности «Математика»

Математические заметки СВФУ Январь—март, 2020. Том 27, № 1

УДК 519.6

АЛГОРИТМ РЕШЕНИЯ ПРЯМОЙ КИНЕМАТИЧЕСКОЙ ЗАДАЧИ СЕЙСМИКИ В ТРЕХМЕРНЫХ НЕОДНОРОДНЫХ ИЗОТРОПНЫХ СРЕДАХ А. А. Галактионова, А. С. Белоносов

Аннотация. Решение прямой кинематической задачи сейсмики — важный этап обработки данных полевых измерений в сейсмологии и сейсморазведке, например, при решении обратной кинематической задачи сейсмики методом нелинейной лучевой сейсмотомографии [1], в методе пространственно-временной миграции Кирхгофа и др. К настоящему моменту существует несколько подходов к решению данной задачи. К наиболее популярным из них относятся: 1) методы, основанные на решении двухточечной краевой задачи; 2) метод волновых фронтов; 3) разностные методы решения уравнения эйконала и др. [2,3]. Предлагаемый в работе метод относится к методам лучевого трассирования [4] и преодолевает некоторые ограничения, присущие вышеперечисленным методам.

Основой подхода, предложенного в данной работе, является итерационный метод Ньютона — Канторовича. Алгоритм трехмерной пристрелки прост в реализации и допускает «почти» линейное распараллеливание.

Вопросы сходимости метода Ньютона в применении к поставленной задаче в данной работе не рассматриваются. Также не рассматриваются вопросы применимости метода при наличии каустик, зон тени, дифракции.

Б01: 10.255877SVFU.2020.12.61.004

Ключевые слова: лучевой метод, лучевое трассирование, ЗО пристрелка.

1. Постановка задачи

Решается задача об определении времени прихода волны в трехмерной неоднородной изотропной среде, поэтому тип волны не имеет значения.

Считается заданным источник волны: это может быть точечный источник или заданное начальное положение фронта волны, а также скоростное строение среды. Мы рассматриваем несколько способов задания среды: если скорость в среде постоянна или линейна, скорость задается явной формулой, в этом случае и для времени вдоль луча определены явные формулы — этот случай в данной статье мы не описываем; если скорость задана гладкой функцией у(х, у, г) € С2, приведен алгоритм для вычисления лучей и производных; если в этой гладкой среде определены границы отражения или преломления, используется модификация предыдущего алгоритма; также приводится способ решения задачи, если скорость задана в узлах пространственной сетки.

© 2020 Галактионова А. А., Белоносов А. С.

Рис. 1. Постановка задачи.

Также задана 3Б сетка с узлами в точках ж^у^г^ определяющая расстановку приемников. В такой постановке может рассматриваться расстановка датчиков на поверхности, размещение их во внутренних точках среды (с погружением датчиков в скважины) или комбинация этих подходов. Требуется определить времена прихода волн в узлы этой пространственной сетки. Постановка задачи схематически показана на рис. 1.

2. Дифференциальное уравнение луча

Луч — кривая, которая в каждой точке ортогональна фронту волны. Для реализации лучевого метода [5] необходимо уметь рассчитывать луч в заданной среде. Последующий расчет времен будет производиться с помощью алгоритма лучевой пристрелки, описанного далее, для реализации которого потребуется знание производных функции луча и направляющего вектора в каждой точке среды. Дальнейшие рассуждения будут производиться с этим расчетом.

Приведем вывод дифференциальных уравнений луча. В качестве отправ-

ного возьмем известное соотношение

„ п

V* =

V

где ¿(ж, у, г) — скалярное поле времени в среде (время «пробега» волны от источника до точки г с координатами (ж, у, г)); п — единичный касательный вектор к лучу в этой точке; v(x,y, г) — скорость в точке г. Заметим, что

с1 /п

¿в V V

V

V

VI)

V2

Отсюда

V • п = п (Vv • п) — Vv (п • п).

Учитывая, что п — это единичный касательный вектор к лучу г, приходим к системе дифференциальных уравнений луча

• п) — Уи(п • п)].

1

V

3. Сглаживание скоростной модели с помощью низко-частотной фильтрации

В общем случае скоростная среда задана в узлах равномерной прямоугольной сетки точек Vijk. Чтобы построить гладкую функцию v(x,y,z), которая достаточно хорошо приближает сеточную, нужно использовать какой-либо алгоритм гладкой аппроксимации. Для удаления высокочастотных колебаний, вызванных случайными ошибками в данных, спектр базисной функции должен быстро уменьшаться в высокочастотной области. Мы используем оператор, который описывает работу цифро-аналогового преобразователя (ЦАП), хорошо известного в теории обработки цифровых сигналов [6]:

те

х

i= — oo

СО = х(и)Б(1 - и),

где Б — импульсная характеристика.

В качестве Б часто используют так называемые оконные Бшс-фильтры во временной области:

где — частота срезки, ^(х) — оконная функция, С — нормировочный коэф-

N ^

фициент такой, что / Б(х) = 1, таким образом Б(0) = 1. — N

х

Наиболее часто используемым является окно Блэкмана:

ш(ж) = 0.42 + 0.5 С08(2пх) + 0.08 сов(4пж), —0.5 < ж < 0.5; ш(ж) = 0, ж ф [—0.5, 0.5].

Но окно Блэкмана имеет только одну непрерывную производную, поэтому нам требуется другая функция окна. Для этих целей подойдет свертка окна Блэкма-на с самим собой. Она имеет четыре непрерывных производных, что достаточно для данного алгоритма.

Идея данного оператора легко обобщается на многомерный случай:

^ ^ ^

v(ж,У,z) = £ £ VijkВх(ж — жг)В2(у — у,-)Вз(г — гй),

Вх, В2, Вз могут быть различными.

Задача вычисления скорости и ее производных разбивается на два этапа. На первом этапе решается задача сглаживания и вычисления значений полученной функции и ее производных в узлах исходной сетки. После этого вычисление скорости и ее производных в произвольной точке производится путем трикуби-ческой интерполяции по следующей формуле:

ззз

/ (ж,у,г) = Х^о жу0 гк

г=0 0=0 к=0

где коэффициенты ао находятся с помощью значений в узлах сетки.

4. Система уравнений

В произвольных средах луч будем вычислять численно, а значит одного лишь решения системы (1) недостаточно для реализации алгоритма лучевой пристрелки. Поэтому к системе дифференциальных уравнений луча добавляются уравнения для производных векторов г и п по 0, Для начала перейдем в уравнении (1) от параметризации по длине дуги к параметризации по времени. Для этого воспользуемся простым соотношением:

¿в ¿в V'

где / — произвольная функция, из которого очевидно, что правая часть уравнений (1) умножится на V. Далее точка будет обозначать производную по времени. С учетом вышесказанного система теперь выглядит так:

( г = V • п,

• т ) (2)

[ п = n(Vv • п) — Vv.

Добавим уравнения для производных по 0, Знание этих производных в каждой точке луча позволяет воспользоваться методом Ньютона, что гарантирует квадратичную сходимость алгоритма пристрелки. С другой стороны, оно

позволяет рассчитывать лучевое расхождение вдоль луча, а следовательно, и амплитуды:

Г в = П0 • V + п • (V« • Те),

Гф = п^ • V + п • (V« • Тф), п = пе • (V« • п) + п •

11 ф = п^ • (V« • п) + п •

V« • п ) + (V« • пе)

£

Л9

-^-Уг; • п ) + (V« • пф)

А9 '

(3)

где

— V« = [(Уг^-гД (Уг;у-гф), (Уг;г-гф)], = [(У^-ге), (\7уу-гв), (Уг^-ге)].

а^ «У

Система уравнений (2), (3) позволяет находить значения г, п, тф, пф, Те, пе в произвольной точке среды. Для того чтобы решить систему, требуется задать граничные условия.

5. Граничные условия

В начальный момент времени положение фронта волны определяется поверхностью Ф(0, <^). Вектор г определяет точку на поверхности, а п — направление выхода луча. Зная параметризацию Ф(0, <^), несложно найти выражения для вектора г и его производных. Вектор п определяется следующим образом:

_ [г6Х1у] [гв хг^]| •

Для нахождения производных вектора п воспользуемся известными соотношениями [7] через коэффициенты первой и второй квадратичной формы:

пе

-ЕС + МЕ

ЕС - Е2 _ ЖЕ - МС * ~ ЕС-Е2

• Тф +

ЕЕ-МЕ ЕС — Е2

• те,

-ЖЕ + МЕ + ЕС-Е2 'Ге'

где Е, Е, С — коэффициенты первой квадратичной формы, Е, М, N — второй.

В случае точечного источника вектор т определяет координаты источника, из чего ясно, что производные тф, те равны нулю. Вектор п определяет направление луча в источнике, например, с помощью сферических координат. Зная выражение для вектора п, можно найти производные.

Полученная система решается методом Рунге — Кутты 5-го порядка с контролем точности на шаге [8].

6. Отражение/преломление

Пусть волна, распространяющаяся в среде со скоростью «(ж, у, г), падает на границу раздела сред, определяемую в некоторой окрестности точки М уравнением С(ж,у,г) = 0, а «*(ж, у, г) — скорость распространения выбранной волны в среде после отражения-преломления (рис. 3).

Пусть п — единичный направляющий вектор луча падающей волны в точке М, п — единичный направляющий вектор луча отраженной-преломленной волны, р — единичный вектор нормали к границе О в точке М, р =

\VG\-

Перед продолжением расчета необходимо переопределить направляющий вектор луча, а также производные этого вектора и вектора г по 9,

Обозначим через £(9,^) точку пересечения луча с поверхностью О, а через т(9, — время от источника до этой точки.

Значит, при Ь = т(9, будет выполнено

Рис. 3. Новая система координат.

G(x(T(в, р), в, р), у(т(в, р), в, р), z(т(в, р), р)) = 0. Дифференцируя соотношение (4) по в и по р, находим т' и т^:

(4)

р •г

Гд =--;

p ■ П p ■ rt

Найдем выражения для вектора направления луча и его производных в новой среде.

Согласно закону отражения-преломления векторы n, п и p лежат в одной плоскости, поэтому вектор п можно представить в виде линейной комбинации векторов p и p х (n х p):

п = Ci ■ p + C2 ■ p х (n х p) (5)

для некоторых Ci, C2. Найдем их значения. Умножив скалярно обе части равенства (5) на вектор p х (n х p) (ортогональный к p), используя свойства скалярного произведения, а также учитывая, что направления векторов n х p и п х p одинаковы, имеем

П ■ (p х (n х p)) = (п х p) ■ (n х p) = |п х p||n х p|

= C2 ■ |p х (n х p)|2 = C2 ■ |n х p|2, откуда, используя закон Снеллиуса, находим выражение для C2 в виде ^ \r] х р\ sin¡3 v*(M) 1 |n х p| sin a v(M) ш

Найдем Ci .

(6)

С\ = г) ■ р = 5 ■ sgn(n • р)\г) ■ р\ = 6 ■ sgn(n • р)л/1 — \г) х р\2, где 6 — условное число, определяемое по правилу

-1, в случае отражения, 1, в случае преломления.

Используя равенства (6), для С получаем следующее выражение:

С = 6 • • р)\/1 -

|п х р|2

£

ш'

где

£ = 6 • sgп(n • VС) 4 / ш2 —

УС

п

= 6 • sgп(n • VС)J |п •

УС |УС|

|2 + ш2 — 1

• ^ л/(п • ус)2 + (с2 - 1)(ус • ус)

|VС|

• (7)

Окончательно 1( VС (

V

с

+ п

V |VС| у

Вычислим производные вектора п по У,

1 / ус ш ^С|2

£— п^С + п . (8)

Пф

где

Ч • / 1

---V -'П-т +-

шш

+

" ус "

_|ус|2_ V

£ —

и.

с

V |VС|2

ус |УС|

ус v

+ п • тф + пф

е~ф = ?[(п • ус) (п • ус); + ш • ^|ус|2 + - 1)([ус]ф, ус)],

(п • VС)ф = ((п • т' + пф) • VG) + (п • ^С]ф),

(9)

(10)

/ = (V« — ш • V«*

■ т' + т'

у т

*1

с

[ус]ф • |ус|2 - 2ус • ([ус]ф • ус) |УС|4

[VС]ф

^С]фК • т' + ж ф] + ^С];[у' • т' + уф] + ^С]; ^ • т' + гф], ^С]Х[ж£ • т' + жф] + ^С];[у • т' + уф] + ^С];И • т' + гф], К • т' + жф] + ^С];[у' • т' + уф] + ^С];[г^ • т' + гф]

1 Г

= ( V«* - -V«, -

ш

•т ' + т '

V ^ V

(11) (12)

Обозначим через р вектор т после границы и найдем его производные на границе. Поскольку это координаты точки, очевидно, на границе т = р. Продифференцируем это равенство по

Р т • т' + т ф = р

тф + Р

рф = т' (Г — Р)+ т ф (13)

Векторы Г и р определяются с помощью уравнения (2):

т = « • п, р = « • п.

Выражения для производных р и п по У определяются с точностью до замены ^ на У в уравнениях (9), (13) и сопутствующих выражениях (10)-(12).

2

2

п

£ф —

п

ш

«

V

Г

7. Алгоритм ЗО пристрелки

Алгоритм двумерной пристрелки описан в работе [9]. Не будем на нем подробно останавливаться, отметим ключевые моменты и выделим основные отличия метода 3Б пристрелки [10] от него.

В алгоритме двумерной пристрелки предполагается, что сетка приемников задана на гладкой поверхности О(ж, у, г). Это накладывает существенные ограничения на задачу, поскольку реальные отражающие горизонты могут быть негладкими. Для предложенного алгоритма 3Б пристрелки не требуется задания функции О, необходимы лишь координаты приемников, ввиду чего это ограничение снимается.

Точку пространства в 3Б случае будем рассматривать как функцию трех независимых переменных г(^, 9, Ь), где (^,9) — параметры, характеризующие начальное положение точки на фронте волны (или углы выхода в случае точечного источника), Ь — время вдоль луча. В этом состоит существенное отличие от метода двумерной пристрелки, где вводилась функция точки пересечения луча с поверхностью О, а для расчета времени вдоль луча требовались дополнительные вычислительные ресурсы.

Рис. 4. Постановка задачи для (а) 20 пристрелки и (Ь) ЗО пристрелки.

Задача схематически показана на рис. 4 и может быть сформулирована следующим образом: известны значения параметров 9о, Ьо, при которых луч приходит в точку Мо = г(^о, 9о, ¿о)- Необходимо определить значения параметров, при которых луч придет в близкую точку Мь

В окрестности точки Мо вектор-функцию г(^, 9,Ь) можно разложить по формуле Тейлора, взяв всего один член:

дг дг д г

гв, I) « г^о, в0, ¿о) + + + — А1.

Считая точку г(<^, 9, Ь) достаточно близкой к точке г(<^о, 9о, ¿о), найдем приращения (Д<^, А9, ДЬ), решив линейную систему из трех уравнений с тремя неизвестными, и рассчитаем луч для этих приращений.

Получим точку Г! = т(^о + Д^, Уо + ДУ, ¿о + Д£). Повторяя эту операцию несколько раз, получим последовательность точек т^, стремящуюся к точке М!.

Алгоритм метода 3Б пристрелки представляет собой известный итерационный метод Ньютона — Канторовича, что при достаточно близком начальном приближении обеспечивает квадратичную сходимость алгоритма.

8. Результаты расчетов

Для начала продемонстрируем, как распространяется фронт волны в средах с границами, рассчитанный по описанному в п. 3.3 алгоритму. Была проведена серия экспериментов по расчету фронтов для нескольких различных моделей, отличающихся положением точечного источника, скоростями в среде и поверхностями, на которых происходит отражение или преломление. Здесь приводим только два наглядных примера.

Для наглядности среды были выбраны кусочно-однородными, конкретные значения скоростей указаны в подписях к рисункам. Положения фронта фиксировались в различные моменты времени.

Красной точкой на рисунках отмечен точечный источник, порождающий волну, синим цветом — фронт волны. На изображениях видно, как при прохождении границ фронт волны меняет свою форму (рис. 5, 6). Отметим, что приведенный алгоритм применим и для случая неплоских границ (см. рис. 5).

Далее представим результаты численных экспериментов для метода 3Б пристрелки.

В табл. 1 приведены результаты расчета точности пристрелки из одной точки в другую. Точностью здесь мы называем расстояние между точкой, к которой мы «пристреливаемся», и точкой, которая была получена с помощью алгоритма в зависимости от номера итерации (первый столбец) и исходного расстояния |МоМ!| (первая строка). Несложно заметить, что чем меньше исходное расстояние между точками, тем быстрее алгоритм сходится. Потому следует подобрать такой порядок обхода сетки, при котором каждый следующий узел будет близок к предыдущему (например, «змейкой»). Также результаты, представленные в таблице, численно подтверждают квадратичную сходимость алгоритма.

Ввиду быстрой сходимости алгоритма представляется возможным ограничиться вычислением всего 2-3 лучей для одного приемника. Поскольку нас интересует лишь время вдоль луча, после последней итерации достаточно рассчитать только приращения аргументов, а не сам луч.

Расчет времен вдоль луча также производился для нескольких моделей, отличающихся начальным положением фронта, скоростями в среде и поверхностью, на которой задана сетка приемников.

В качестве начального положения фронта использовался один из следующих вариантов:

• единичная сфера с центром в точке (1,1, —0.5), далее — сфера (рис. 7-9);

Ь = 2.1с Ь = 2.7с

Рис. 5. Преломление. Граница — внешняя сторона сферы, точечный источник, скорость до преломления 1.5 м/с, после преломления — 0.8 м/с.

Ь = 0.3 с

/

/ /

Ь = 1.2с

Ь = 2.1с

Ь = 3.9с

Рис. 6. Преломление. 2 плоские границы, точечный источник, скорость до первой границы 1.9 м/с, до второй — 1.5 м/с, далее — 0.9 м/с.

Таблица 1. Зависимость точности (расстояния между полученной точкой и точкой, к которой пристреливаемся) от номера итерации и исходного расстояния |М"оМх|

Шаг \|МоМ1| 1.19 .92 3.62 4.29

1 1.02е-001 5.84е-001 1.61е+000 2.75е+000

2 8.29е-004 2.99е-002 1.90е-001 5.25е-001

3 5.04е-008 7.30е-005 3.07е-003 3.07е-002

4 9.85е-014 4.10е-010 7.81е-007 7.13е-005

5 9.68е-014 7.57е-014 3.60е-014 3.94е-010

• точечный источник с координатами (1,1, -0.5), далее — точечный источник (рис. 10, 12);

(о +

• плоскость Ф(0, = < 2.25 + 0.50 — 3^, далее — плоскость (рис. 11).

I —0

В качестве скоростной модели среды использовался один из следующих вариантов:

• скорость, заданная аналитической формулой г>(х, у, г) = 1.5 + 0.1х + 0.1у + г, далее — линейная скорость (см. рис. 7-9, 12);

• скорость, заданная аналитической формулой г>(х, у, г) = 1.5 + 0.1х + 0.1у + г — 0.5соэ-^, далее — нелинейная скорость (см. рис. 10, 11).

В качестве сетки использовался один из следующих вариантов:

• сетка на поверхности С(ж, у, г) = 1 — + зш(47ту) — г, далее — гладкая поверхность (см. рис. 7, 8, 10, 11);

(г = —1, х < 5.7, у е Д, х = 5.7, г е ( — 1,1), у е Д, , да-г = 1, х > 5.7, у е Д,

лее — ступенька (см. рис. 9);

• сетка на поверхности С(х, у, г) = |3х — [3х] — 0.51, у е Д ([•] — целая часть числа), далее — «стиральная доска» (см. рис. 12).

Во всех случаях сетка задавалась размером [100 х 100] точек.

Цветом на изображенных поверхностях обозначено время прихода луча в данную точку: от синего к красному — от меньшего к большему. Синими точками на поверхности Ф отмечены начальные точки лучей, попадающих в приемники. На всех рисунках тонкой черной линией от источника до одного из приемников показана траектория одного из лучей.

Отметим, что проверка предложенного в работе алгоритма была произведена для сеток не только на гладких поверхностях, но и на поверхностях, претерпевающих разрывы производных (ступенька, «стиральная доска»). К таким

Рис. 7. Линейная скорость, начальное положение фронта — внешняя сторона сферы, гладкая поверхность (а) вид сверху и (Ь) вид под углом.

Рис. 8. Линейная скорость, начальное положение фронта — внутренняя сторона сферы, гладкая поверхность (а) вид сверху и (Ь) вид под углом

Рис. 9. Линейная скорость, начальное положение фронта — внешняя сторона сферы, ступенька (а) вид сверху и (Ь) вид под углом.

сеткам невозможно применить алгоритм двумерной пристрелки, а результат применения описанного алгоритма для них представлен на рис. 9, 12.

Заключение

Разработан алгоритм 3Б лучевой пристрелки как для случая точечного источника, так и в более общей постановке, когда начальное положение фронта волны задано в виде гладкой поверхности. По сравнению с методом двумерной пристрелки алгоритм проще в реализации, а также может быть использован в случае, когда «приемники» располагаются на негладкой поверхности.

(а) (Ь)

Рис. 10. Нелинейная скорость, начальное положение фронта — точечный источник, гладкая поверхность(а) вид сверху и (Ь) вид под углом.

(о) (б)

Рис. 11. Нелинейная скорость, начальное положение фронта — плоскость, гладкая поверхность (а) вид сверху и (Ь) вид под углом.

(а) (Ь)

Рис. 12. Линейная скорость, точечный источник, «стиральная доска» (а) вид сверху и (Ь) вид под углом.

В отличие от первого метода для реализации 3Б пристрелки требуются лишь координаты точек.

Метод может быть также применим при наличии в среде границ, на которых происходит отражение или преломление. Алгоритм расчета лучей, распространяющихся в средах с границами, выведен и представлен в работе.

Создан программный комплекс для расчета времен прихода волн в заданную 3Б сетку точек. Приводятся результаты численных экспериментов.

ЛИТЕРАТУРА

1. Юрченко М. А., Белоносова А. В., Белоносов А. С. Об одном алгоритме решения обратной кинематической задачи сейсмики // Сиб. электрон. мат. изв. 2013. № 10. C. 74—86.

2. Gjoystdal H., Iversen E., Laurain R. et al. Review of ray theory applications in modelling and imaging of seismic data // Stud. Geoph. Geod. 2002. V. 46. P. 113-164.

3. Cerveny V. Seismic Ray Theory. New York: Camb. Univ. Press, 2001.

4. Rawlinson N., Hauser J., Sambridge M. Seismic ray tracing and wavefront tracking in laterally heterogeneous media // Adv. Geophys. 2007. V. 49. P. 203-267.

5. Алексеев А. С., Гельчинский Б. Я. О лучевом методе вычисления полей волн в случае неоднородных сред с криволинейными границами раздела // Вопросы динамической теории распространения сейсмических волн. 1959. № 3. C. 107-159.

6. Rabiner L. R., Gold B. Theory and application of digital signal processing. Ch. 5. P. 300-302. Prentice Hall, 1975.

7. Погорелов А. В. Дифференциальная геометрия. М.: Наука, 1974.

8. Lapidus L. Numerical solution of ordinary differential equations. New York: Acad. Press, 1971.

9. Цецохо В. А., Белоносова А. В., Белоносов А. С. Расчетные формулы линейного геометрического расхождения при лучевом трассировании в трехмерной блоково-неоднородной градиентной среде // Сиб. журн. вычисл. математики. 2009. T. 12, № 3. C. 325-339.

10. Galaktionova A., Belonosov A.. Computation of seismic wave field kinematics in a three-dimensional heterogeneous isotropic medium // Application of mathematics in technical and natural sciences: AIP Conf. Proc. 2015. V. 1895. P. 120004.

Поступила в редакцию 3 ноября 2019 г. После доработки 6 февраля 2020 г. Принята к публикации 17 февраля 2020 г.

Галактионова Анастасия Андреевна, Белоносов Андрей Сергеевич

Институт вычислительной математики и математической геофизики СО РАН,

пр. Лаврентьева, 6, Новосибирск 630090;

Новосибирский государственный университет,

ул. Пирогова, 1, Новосибирск 630090;

a.galaktionova@g.nsu.ru, white@sscc.ru

Математические заметки СВФУ Январь—март, 2020. Том 27, № 1

UDC 519.6

THE ALGORITHM FOR SOLVING THE DIRECT KINEMATIC PROBLEM OF SEISMICS IN THREE-DIMENSIONAL HETEROGENEOUS ISOTROPIC MEDIA A. A. Galaktionova and A. S. Belonosov

Abstract: The solution of the direct kinematic problem of seismics is an important stage of the seismic data processing in seismology and seismic explorations, e.g., while solving the inverse kinematic problem of seismics by the method of nonlinear seismic tomography, in the method of the space-time Kirchhoff's migration, etc. There are several approaches now to the solution of this problem: 1) methods based on the solution of two-point boundary value problems; 2) the wavefront construction method; 3) finite-difference methods for solving eikonal differential equations, etc. [2—3]. The method proposed in this paper relates to ray tracing methods and overcomes some of the limitations inherent in the above methods.

The basis of the approach is the iterative Newton's method. The algorithm of three-dimensional shooting is simple to implement and allows "almost" linear speedup from parallelization. Questions of convergence of the Newton method applied to this problem are not considered in this paper. In addition, we do not consider the applicability of the method in the presence of caustics, shadow zones, and diffraction.

DOI: 10.25587/SVFU.2020.12.61.004

Keywords: ray-based method, ray tracing, 3D shooting.

REFERENCES

1. Yurchenko M. A., Belonosov A. S., and Belonosova A. V., "On an algorithm to solve an inverse kinematic problem of seismics," Sib. Electron. Math. Rep., 10, 74—86 (2013).

2. Gjoystdal H., Iversen E., Laurain R. et al., "Review of ray theory applications in modelling and imaging of seismic data," Stud. Geoph. Geod., 46, 113—164 (2002).

3. Cerveny V., Seismic Ray Theory, Camb. Univ. Press, New York (2001).

4. Rawlinson N., Hauser J., and Sambridge M., "Seismic ray tracing and wavefront tracking in laterally heterogeneous media," Adv. Geophys., 49, 203—273 (2007).

5. Alekseev A. S. and Gelchinsky B. Ya., "Ray-tracing method for the calculation of wave fields in inhomogeneous media with curved interfaces [in Russian]," in: Problems in Dynamic Theory of Seismic Wave Propagation, No. 3, 107-159 (1959).

6. Rabiner L. R. and Gold B., Theory and Application of Digital Signal Processing, ch. 5, pp. 300-302, Prentice Hall (1975).

7. Pogorelov A. V., Differential Geometry [in Russian], Nauka, Moscow (1974).

8. Lapidus L., Numerical Solution of Ordinary Differential Equations, Acad. Press, New York (1971).

9. Tsetsokho V. A., Belonosova A. V., and Belonosov A. S., "Calculation formulas of linear geometrical spreading at ray tracing in a 3D block-inhomogeneous gradient medium [in Russian]," Sib. Zh. Vychisl. Mat., 12, No. 3, 325-339 (2009).

© 2020 A. A. Galaktionova, A. S. Belonosov

10. Galaktionova A. and Belonosov A., "Computation of seismic wave field kinematics in a three-dimensional heterogeneous isotropic medium," in: Application of Mathematics in Technical and Natural Sciences, AIP Conf. Proc., 120004, 1895 (2017).

Submitted November 3, 2019 Revised February 6, 2020 Accepted February 17, 2020

Anastasiia A. Galaktionova and Andrey S. Belonosov

Institute of Computational Mathematics and Mathematical Geophysics,

6 Acad. Lavrentiev Avenue, Novosibirsk 630090, Russia;

Novosibirsk State University,

1 Pirogov Street, Novosibirsk 630090, Russia

a.galaktionova@g.nsu.ru, white@sscc.ru