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

УДК 004.021:550.34

А. С. Сердюков, А. А. Дучков, А. С. Смелов

Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН пр. Академика Коптюга, 3, Новосибирск, 630090, Россия

Новосибирский государственный университет ул. Пирогова, 1, Новосибирск, 630090, Россия

[email protected] [email protected]. т, s.art.s@mail.т

АЛГОРИТМ СЕЙСМИЧЕСКОЙ МИГРАЦИИ В ОБРАТНОМ ВРЕМЕНИ НА ОСНОВЕ ЭФФЕКТИВНОГО ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ПЕРВЫХ ВСТУПЛЕНИЙ СЕЙСМИЧЕСКИХ ВОЛН *

Предлагается новый численный алгоритм моделирования распространения сейсмических волновых полей. Конечно-разностное решение уравнений упругости происходит только в полосе, следующей за фронтом первых вступлений, которую далее будем называть «окном». Положение «окна», т. е. фронта первых вступлений, рассчитывается путем численного решения уравнения эйконала. Предлагаемый подход «оконного» моделирования позволяет существенно ускорить вычисления (так как конечно-разностные вычисления проводятся в бегущем «окне», а не во всей расчетной области) и уменьшить объем оперативной памяти, необходимый для хранения поля. На тестах показана возможность ускорения до 40 раз и экономии памяти до 50 раз. Данный подход позволяет моделировать только волновые формы первых вступлений, но этого оказывается достаточно для многих процедур сейсмической миграции. В частности, в статье рассматривается распространенная процедура миграции в обращенном времени, реализация которой требует больших объемов вычислений. Показано, что кроме ускорения миграции, данный подход позволяет избежать появления «артефактов» - ложных изображений несуществующих границ.

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

Введение

Численное моделирование распространения упругих волн в неоднородных средах широко применяется в разных процедурах обработки и миграции сейсмических данных. При этом для решения уравнений упругости, т. е. моделирования волновых полей, широко применяются конечные разностные методы. В силу гиперболичности (t-гиперболичности) уравнений упругости для их численного решения можно использовать явные схемы, т. е. на каждом текущем шаге по времени решение рассчитывается на основе известных значений, рассчитанных на предыдущих шагах по времени. Наиболее широко распространенным подходом является схема второго порядка на сдвинутых сетках [1]. Для подавления отражений от границ расчетной области на них применяется искусственное поглощение, например, так называемые идеально согласованные слои (perfectly matched layers, или PML) [2]. Для реалистичных

* Работа выполнена при финансовой поддержке РФФИ (грант № 14-05-00862-а).

Сердюков А. С., Дучков А. А, Смелов А. С. Алгоритм сейсмической миграции в обратном времени на основе эффективного численного моделирования первых вступлений сейсмических волн // Вестн. Новосиб. гос. ун-та. Серия: Информационные технологии. 2016. Т. 14, № 4. С. 68-79.

ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2016. Том 14, № 4 © А. С. Сердюков, А. А. Дучков, А. С. Смелов, 2016

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

В работах [3, 4] была высказана идея, что для моделирования волновых форм первых вступлений сейсмических волн можно в каждый момент времени проводить конечно-разностные вычисления не во всей расчетной области, а только в полосе, следующей за фронтом первых вступлений, которую далее будем называть «окном». Для предварительного расчета фронтов, т.е. времен первых вступлений, используется численное решение уравнения эйконала, которое требует гораздо меньше вычислительных ресурсов. Существуют разные численные алгоритмы решения уравнения эйконала, включая Fast Marching Method [5] и Fast Sweeping Method [6]; последний метод допускает эффективную параллельную реализацию [7]. В статье мы предлагаем новый вычислительный алгоритм для реализации идеи моделирования волнового поля в «окне». Основной особенностью алгоритма является пересортировка используемых массивов в порядке возрастания времен прихода первых вступлений, что позволяет повысить локальность доступа к памяти при расчете в «окне» и существенно сократить время вычислений. Другим преимуществом предлагаемого подхода является существенное сокращение объемов памяти, необходимой для хранения волнового поля в окрестности первых вступлений с целью дальнейшего использования.

В качестве приложения в статье рассматривается миграция в обращенном времени (reverse time migration, или RTM) [8, 9] - метод, применяемый для трансформации сейсмических данных в изображение отражающих границ. Суть метода сводится к перемножению двух волновых полей - «прямого» поля от известного источника и «обращенного» поля, которое получающегося путем продолжения данных, зарегистрированных в приемниках, обратном времени. «Прямое» поле при этом необходимо хранить в памяти для дальнейшего построения изображения среды, но для него полезными являются только волновые формы первых вступлений. Таким образом, расчет и хранение «прямого» поля идеально подходят для реализации с помощью метода расчета в «окне».

Метод моделирования волновых форм

первых вступлений сейсмических волн

Для простоты изложения будем рассматривать скалярное волновое уравнение:

—- ип — Ли= 5 (х- xs)f(t), (1)

где и(х, £) - решение волнового уравнения с нулевыми начальными данными: и = 0, = 0 при £ = 0 (задача Коши). Правая часть соответствует точечному источнику, расположенному в точке х3, с(х)-скорость распространения сейсмических волн.

Идея моделирования в «окне» заключается в том, что расчет волнового поля выполняется на каждом шаге по времени не во всей области, а только полосе, следующей за фронтом первых вступлений. Для этого на первом этапе алгоритма происходит расчет времен пробега во всей расчетной области на основе конечно-разностного решения уравнения эйконала [5, 6]:

^ - тцту (2)

при условии инициации схемы в источнике х3:

t(Xs) = 0.

(3)

I_I_I_I_|_|

О 200 400 600 ВСЮ

Расстояние, м

200

Н 4500

1000

3000

1000 -

1200

2500

1200 "

2000

200 400 000 ВСЮ О

Расстояние, м

200 400 600 Расстояние, м

мУс

5000 О

400

- 4000

I

300

3500

юооЬ

1200 г

Рис. 1. Скоростная модель и фронты времен первых вступлений (а); «окна» для вычисления волнового поля в окрестности первых вступлений (б, в)

Примеры расчета фронтов волн первых вступлений (изолинии времен пробега т(х)) в контрастной модели показаны на рис. 1, а. Синтетическая скоростная модель представляет собой среду с высокоскоростным включением (с = 5000 м/с). Скоростная модель показана цветом, расходящиеся от источника фронты - черными линиями.

На втором этапе происходит численное решение волнового уравнения (1) в «окне» - полосе, следующей за фронтом первых вступлений, т. е. для момента времени ^ расчет ведется в полосе: 0< Ь — т(х) <ц. Соответствующие «окна» показаны для разных моментов времени на рис. 1, б-в (вычисления проводятся только внутри чёрной полосы). Видно, что в высокоскоростном слое ширина «окна» увеличивается, так как соответствует длительности сигнала /(£). Для примеров в статье в источнике использовался импульс Рикера с доминирующей частотой 30 Гц.

Заметим, что рассчитывать поле перед движущимся окном не имеет смысла, так как сигнал не мог туда попасть в силу принципа причинности. Волны в области поле окна нет смысла рассчитывать, так как они не смогут ее догнать и начать вносить вклад в волновое поле в движущемся «окне».

Численную реализацию конечно-разностных вычислений в «окне» рассмотрим для простоты изложения в двумерном случае. Пусть скоростная модель задана на регулярной прямоугольной сетке с узлами (¿,у). Применяя стандартные алгоритмы численного решения уравнения эйконала, получим значение времени пробега в каждой точке расчетной сетки т(£,у). Затем методом конечных-разностей решаем волновое уравнение на каждом шаге по времени только для точек сетки, удовлетворяющих условию:

пЛЬ < т(£,у) < пЛЬ + ц, (5)

где п, ЛЬ - номер и размер шага по времени, а константа д определяет ширину вычислительной полосы (определяется длительностью сигнала в источнике).

На каждом шаге по времени конечно-разностной схемы вычислительное окно может быть определено путем проверки условия (5). Именно такой подход был предложен в [10]. Хотя он позволяет существенно снизить объем вычислений конечно-разностной схемы, но все равно на каждом временном шаге для каждой точки (¿,у) проводится проверка условия (5). Кроме того, пропуск вычислений конечно-разностной схемы в некоторых точках приводит к тому, что доступ к памяти становится менее систематическим, что не позволяет системе организовать оптимальное кэширование доступа к памяти.

Для более эффективной реализации алгоритма расчета в «окне» мы предлагаем использовать альтернативное линейное индексирование расчетной сетки со взаимно-однозначным соответствием

(Ху) ^ к (6)

между пространственной декартовой индексацией (¿,у) и линейным индексом к, соответствующим возрастанию времени пробега, т.е. для каждого к:

т(¿(к + 1), ](к +1)) > т{Цк), ;(Ю). (7)

Все массивы, участвующие в конечно-разностных вычислениях, физически пересортировываются по возрастанию номера к. На каждом шаге по времени п вычислительное «окно» определяется последовательным диапазоном изменения индекса к:

кьедт — к < кепа, (8)

где начальные и конечные значения кЬед1п и кепа определяются на предварительном этапе сразу после численного решения уравнения эйконала с использованием списка (7) и условия (5).

72

А. С. Сердюков и др.

Решение волнового уравнения в текущий момент времени (в активном «окне») хранится в памяти последовательно (по возрастанию индекса к), и, таким образом, может быть в том же порядке последовательно записано в каждый момент времени в одномерный массив. Таким образом, при вычислении волнового поля в «окне» достигается локальность доступа к памяти, что значительно ускоряет вычисления. Решение уравнения (1) в «окне» рассчитывается методом конечных разностей с использованием явной схемы на сдвинутых сетках [1].

Было проведено тестирование алгоритма моделирования волнового поля в «окне» для анализа его эффективности по сравнению со стандартным методом конечных элементов при проведении вычислений во всех точках расчетной области. Результаты тестирования по ускорению вычислений показаны в табл. 1, а по экономии памяти (при сохранении волнового поля на каждом шаге по времени - в табл. 2. Рассматривалась однородная среда со скоростью 2000 м/с. Шаг пространственной сетки брался равным 2,7 м. В источнике использовался импульс Рикера с доминирующей частотой 30 Гц. В каждом из расчетов источник располагался в центре области, а время моделирования бралось равным 89 % от максимального времени пробега прямой волны для рассматриваемой модели.

Из результатов, приведенных в табл. 1, видно, что ускорение вычислений и экономия памяти для хранения поля увеличивается с ростом размеров модели. Для сетки 4000 X 4000 оконные вычисления оказываются в 40 раз быстрее, а памяти требуется почти в 50 раз меньше.

Таблица 1

Ускорение при расчёте поля в «окне»

Размер сетки Время вычислений, с Ускорение,

во всей области в полосе раз

500 х 500 2,2 0,15 15

1000 х 1000 15,9 0,65 24

2000 X 2000 117 3,36 35

4000 X 4000 877 21,7 40

Таблица 2

Экономии памяти при хранении поля в «окне»

Размер сетки Память Экономия

для полного поля для поля в полосе памяти, раз

500 х 500 467 Мб 73 Мб 6

1000 х 1000 3,7 Гб 0,3 Гб 12

2000 X 2000 30 Гб 1,22 Гб 25

4000 X 4000 240 Гб 4,9 Гб 49

Сейсмическая миграция в обращенном времени на основе моделирования в «окне»

Задача миграции состоит в том, чтобы трансформировать сейсмические данные в изображение отражающих границ, используя известную макроскоростную модель с0(х) (которая оценивается на предшествующих этапах обработки сейсмических данных). Процедура миграции в обращенном времени (ЯТМ, см. [8, 9]) применяется отдельно к каждой сейсмограмме заданного положения источника х3. Для нее вычисляется частичное изображение д(х, х3~) путем интегрирования по времени произведения двух волновых полей:

т

д(х,х3) = /0 щ(х, ^)иаг(х,

(10)

где щ (х, £) - «прямое» поле от точечного источника, расположенного в точке х3 с импульсом /(£), т.е. решение волнового уравнения (1) (для скоростной модели с0) с нулевыми начальными данными (задача Коши); изг(х, £) - «обращенное» поле, порожденное данными в точках приема и продолженное в обратном времени: ^ = Т — т. е. решение задачи Коши:

с нулевыми начальными данными и = 0 и и^ = 0 при £ = 0. Суммирование в правой части (11) ведется по всем приемникам. Функция в правой части ^г(х5, £) представляет собой данные, записанные приемниками в точках хг; Т - максимальное время наблюдений.

Для получения полного изображения необходимо сложить все частичные изображения (для разных положений источников):

где С(х) - полное изображение отражающих границ, которое затем используется при геологической интерпретации.

Заметим, что выбор макроскоростной модели с0(х) зависит от степени изученности района. Как правило, макроскоростная модель берется гладкой, что позволяет получать точные сейсмические изображения в большинстве геологических моделей. Исключение составляет случай солянокупольной тектоники, когда в скоростном разрезе присутствуют контрастные высокоскоростные слои, границы которых можно определить хорошо, но сильно падает точность восстановления целевых отражающих границ под ними.

На рис. 2, а представлена синтетическая скоростная модель среды и геометрия наблюдений: 17 источников типа центра расширения показаны треугольниками, 400 приемников были равномерно распределены вдоль зеленой линии. Контрастный высокоскоростной слой (имитирует слой соли) показан желтым цветом. Для этой модели были сделаны расчеты волновых полей методом конечных разностей для разных положений источников и для линии приемников (зеленая линия на рис. 2, а). В качестве зондирующего сигнала использовался импульс Рикера с центральной частотой 30 Гц. Одна из сейсмограмм показана на рис. 2, б. Эта же сейсмограмма после подавления кратных волн представлена на рис. 2, в. На ней четко видны три однократных отражения от трех границ в модели рис. 2, а. Заметим, что подавление кратных волн - этап предварительной обработки, необходимый перед применением процедуры миграции.

В предположении, что границы контрастного слоя известны можно рассмотреть два варианта задания макроскоростной модели с0(х), как показано на рис. 3. На рис. 3, а показана модель с контрастным слоем, на рис. 3, б - со сглаженным слоем. Именно второй вариант сглаженной модели обычно используется на практике, так как использование контрастной миграционной модели порождает артефакты (ложные изображения границ при миграции). К сейсмограммам после подавления кратных волн была применена процедура миграции в обращенном времени (10)-(12). На рис. 4, а показан результат стандартного подхода к миграции в обращенном времени, т. е. использование сглаженной скоростной модели (см. рис. 3, б) и стандартный расчет волновых полей («прямого» и «обращенного») во всех точках расчетной сетки. На рисунке четко видно изображение трех границ и отсутствие артефактов (ложных границ). На рис. 4, б показан результат стандартного алгоритма миграции, расчет волновых полей («прямого» и «обращенного») во всех точках расчетной сетки, к контрастной скоростной модели (см. рис. 3, а). На рисунке видно, что использование контрастной

1

(12)

500 1000

Расстояние, м

50 103 15» Я» 2!0 300 350 400 Нг*лер канат

Я) 1И 150 2Е0 300 ЙО «

Нпмяр канала

Рис. 2. Синтетическая модель: а - скоростная модель среды и геометрия наблюдений (треугольники - источники, зеленая линия - приемники); б - сейсмограмма для одного источника; в - та же сейсмограмма после подавления кратных волн

V = 2000 м/с р

р = 1600 кг/м'

V - 2000 м/с р

р — 1600 кг/м3

(с 1000 х

¿5

К

V = 4500 м/с

Р

р - 3600 кг/м

со 1000

з:

X

щ

с;

V = 4500 м/с

Р

р = 3600 кг/м

V = 2000 м/с Р

р ~ 1600 кг/м'

V = 2000 м/с Р

р - 1600 кг/м'

500 1000

Расстояние, м

500 1000

Расстояние,м

Рис. 3. Типы макроскоростной миграционной модели: а - модель с контрастным слоем; б - со сглаженным слоем

скоростной модели дает более четкое изображение верхней границы, но порождает большое количество артефактов, особенно в верхней части изображения. Наконец, на рис. 4, в показан результат нового алгоритма миграции в применении к контрастной скоростной модели (см. рис. 3, а). При этом расчет «прямого» воля от источника проводился в «окне». На рисунке видно, что четкое изображение верхней границы сохраняется, а артефакты в верхней части изображения пропадают.

Сравним сейсмические изображения на рисунке 4а и 4в более детально. На рисунке 5 они приведены в увеличенном масштабе и на них наложено истинное положение отражающих границ - пунктирные линии. Из рисунка 5а видно, что использование сглаженной модели приводит к "размазыванию" изображения и неточному положению границ (сдвигу по глубине на 20-30 м). Из рисунка 5б видно, что при использовании миграции с расчетом поля в «окне» изображение с лучшим разрешением и точным положением (по глубине). Заметим, что кроме улучшения качества изображения, миграция с моделированием «прямого» поля в «окне» оказывается в 2 раза быстрее и требует гораздо меньше памяти для хранения. Ускорение только в 2 раза связано с тем, что в «окне» считается только «прямое» поле (от источника), «обращенное» поле (от приемников) рассчитывается стандартным образом. Наличие в изображении на рисунке 5б остаточных артефактов также связано со стандартным подходом к расчету «обращенного» поля от приемников, так что при его продолжении в контрастной скоростной модели возникают кратные волны. Для устранения этих артефактов необходимо в будущем реализовать продолжение «обращенного» поля от приемников в варианте расчета в «окне».

200 400 600 ВСЮ 1000 1200 0 200 400 600 BDO 1000 1200 О 200 400 600 ЭОО 1000 1200

Расстояние, м Расстояние, м Расстояние, м

Рис. 4. Результат миграции в обращенном времени: а - сглаженная модель (см. рис. 3, б) и расчет поля во всей области; б - контрастная модель (см. рис. 3, а) и расчет поля во всей области;

в - контрастная модель (см. рис. 3, а) и расчет поля в «окне»

О 200 400 600 ВОТ 1000 1200 0 200 400 600 В00 1000 1200

Расстояние, м Расстояние, м

Рис. 5. Результат миграции в обращенном времени: а - сглаженная модель (см. рис. 3, б) и расчет поля во всей области; б - контрастная модель (см. рис. 3, а) и расчет поля в «окне»; пунктирные линии - истинное положение отражающих границ

Заключение

В работе предложен новый метод конечно-разностного моделирования и сохранения сейсмического волнового поля в окрестностях фронтов первых вступлений. На каждом шаге алгоритма конечно-разностные вычисления проводятся в «окне», т. е. в узкой полосе, которая следует за фронтом первых вступлений. В результате такого подхода удается существенно ускорить (до 40 раз) прямое моделирование волновых форм первых вступлений и существенно сократить (до 50 раз) объем памяти, необходимый для хранения волнового поля при реализации алгоритма миграции в обращенном времени.

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

Список литературы

1. Virieux J. Wave propagation in heterogeneous media: Velocity-stress finited fference method // Geophysics. 1986. Vol. 51. No. 4. P. 889-901.

2. Kristek J., Moczo P., Galis M. A brief summary of some PML formulations and discretizations for the velocity-stress equation of seismic motion // Studia Geophysica et Geodaetica. 2009. Vol. 53. No. 4. P. 459-474.

78

А. С. Сердюков и др.

3. Boore D. M. Finite difference methods for seismic wave propagation in heterogeneous materials // Methods of Computational Physics, vol. 2 / B. Alder, S. Fernbach, M. Rotenberg. (Eds). New York: Academic Press. 1972. P. 1-37.

4. Vidale J. Finite-difference calculation of travel times // Bulletin of the Seismological Society of America. 1988. Vol. 78. No. 6. P. 2062-2076

5. Sethian J. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge Uni. Press, 1999. 400 p.

6. Zhao H. A fast sweeping method for eikonal equations // Mathematics of computation. 2005. Vol. 74, № 250. Р. 603-627.

7. Никитин А. А., Сердюков А. С., Дучков А. А. Параллельный алгоритм решения уравнения эйконала для трехмерных задач сейсморазведки // Вестн. НГУ Серия: Информационные технологии. 2015. Т. 13, № 3. С. 19-28.

8. Baysal E., Kosloff D. D., Sherwood J. W. Reverse time migration // Geophysics. 1983. Vol. 48. No. 11. P. 1514-1524.

9. McMechan G. A. Migration by extrapolation of time-dependent boundary values // Geophysical Prospecting. 1983. Vol. 31. No. 3. P. 413-420.

10. Kvasnicka M., Zahradnfandëk J. A combined kinematic-dynamic method for fast computations of the first-arriving waveforms // Seismic Waves in Complex 3-D Structures, Dep. Geophys., Charles Uni., Prague, 1996. Rep. No. 4. P. 251-286.

Материал поступил в редколлегию 20.09.2016

A. S. Serdyukov, A. A. Duchkov, A. S. Smelov

A. A. Trofimuk Institute of Petroleum Geology and Geophysics SB RAS 3 Academician Koptyug, Novosibirsk, 630090, Russian Federation

Novosibirsk State University 1 Pirogov Str., Novosibirsk, Russian Federation

[email protected] DuchkovAA@ipgg. sbras. ru, s. art.s@mail. ru

SEISMIC REVERSE-TIME MIGRATION USING EFFICIENT NUMERICAL ALGORITHM FOR MODELING FIRST-ARRIVAL WAVEFORMS

We suggest a new algorithm for numerical modeling of seismic wavefields. Finite-difference solution of the wave equation occurs only in a strip (further called «window») following the front of first arrivals. The «window» position is determined from the numerical eikonal solver by computing first-arrival fronts. The proposed approach of the windowed modeling can significantly speed up calculation (since the finite-difference calculations are carried out in a moving «window» instead of computing in the whole domain) and reduce the amount of memory required for storing the field. The tests showed the possibility of the calculation accelerating up to 40 times and memory usage reduction up to 50 times. This approach allows us to model only the waveforms of the first arrivals, but this is enough for many seismic migration procedures. In particular, the article discusses the case of the reverse-time migration because it requires large amounts of computations. Tests show that in addition to accelerating the migration this approach allows us to avoid the appearance of artifacts, i.e. false images of non-existent seismic boundaries.

Keywords: numerical algorithms, finite-difference scheme, wave equation, eikonal equation, reverse-time seismic migration, parallel computing.

References

1. Virieux J. Wave propagation in heterogeneous media: Velocity-stress finited fference method // Geophysics. 1986. Vol. 51, No. 4. P. 889-901.

2. Kristek J., Moczo P., Galis M. A brief summary of some PML formulations and discretizations for the velocity-stress equation of seismic motion // Studia Geophysica et Geodaetica. 2009. Vol. 53, No. 4. P. 459-474.

3. Boore, D. M. Finite difference methods for seismic wave propagation in heterogeneous materials // Methods of Computational Physics, vol.2, / B. Alder, S. Fernbach, M. Rotenberg. (Eds). New York: Academic Press. 1972. P. 1-37.

4. Vidale, J. Finite-difference calculation of travel times // Bulletin of the Seismological Society of America. 1988. Vol. 78, No. 6. P. 2062-2076

5. Sethian J. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge Univ Pr. 1999. 400 p.

6. Zhao H. A fast sweeping method for eikonal equations // Mathematics of computation. 2005. Vol. 74, № 250. P. 603-627.

7. Nikitin A. A., Serdyukov A. S., Duchkov A. A. Parallel algorithm of three-dimensional eikonal solver for seismic applications // Vestnik NSU Series: Information Technologies. - 2015. -Volume 13, Issue No 3. - P. 19-28. - ISSN 1818-7900. (in Russian).

8. Baysal E., Kosloff D.D., Sherwood J.W. Reverse time migration // Geophysics. 1983. Vol. 48, No. 11. P. 1514-1524.

9. McMechan G. A. Migration by extrapolation of time-dependent boundary values // Geophysical Prospecting. 1983. Vol. 31, No. 3. P. 413-420.

10. Kvasnicka M., Zahradnrandek J. A combined kinematic-dynamic method for fast computations of the first-arriving waveforms // Seismic Waves in Complex 3-D Structures, Dep. Geophys., Charles Univ., Prague, 1996. Rep. No. 4. P. 251-286.