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

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ДИНАМИКИ

ПЕРВЫХ ВСТУПЛЕНИЙ ДЛЯ МЕТОДА ВОЛНОВОЙ ТОМОГРАФИИ

Александр Сергеевич Сердюков

Новосибирский государственный университет, 630090, Россия, г. Новосибирск, ул. Пирогова, 2, кандидат физико-математических наук, доцент кафедры высшей математики физического факультета, тел. +7913-767-95-13, e-mail: AleksanderSerdyukov@yandex. ru

Антон Альбертович Дучков

Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, 630090, Россия, г. Новосибирск, пр. Ак. Коптюга, 3, кандидат физико-математических наук, зав. лаб., e-mail: DuchkovAA@ipgg.sbras.ru

Александр Алексеевич Никитин

Новосибирский государственный университет, 630090, Россия, г. Новосибирск, ул. Пирогова, 2, магистрант/ассистент кафедры параллельных вычислений факультета информационных технологий

В работе предлагается алгоритм моделирования и сохранения волнового поля в каждый момент времени в окрестности первых вступлений. На каждом шаге конечно-разностной схемы вычисления проводятся только в полосе, сдвигающейся вместе с положением фронта, которое заранее вычисляется с помощью конечно-разностного решения уравнения эйконала. Возможность хранения волнового поля позволяет проводить построение ядер чувствительности (sensitivity kernels) которые широко используются в ряде методов обращения волновых форм и времен вступлений.

Ключевые слова: сейсмические волны, численное моделирование, уравнение эйконала, волновая томография.

NUMERICAL MODELLING OF FIRST ARRIVAL WAVEFORMS FOR WAVE EQUATION TRAVELTIME INVERSION

Alexandr S. Serdyukov

Novosibirsk State University, 630090, Russia, Novosibirsk, Pirogova str. 2, PHD, assistant lecturer, department of Physics, tel. +7913-767-95-13, e-mail: AleksanderSerdyukov@yandex.ru

Anton A. Duchkov

A. A. Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 630090, Russia, Novosibirsk, Prospect Akad. Koptyuga 3, PHD, head of the laboratory, e-mail:

DuchkovAA@ipgg. sbras.ru

Alexandr A. Nikitin

Novosibirsk State University, 630090, Russia, Novosibirsk, Pirogova str. 2, master student, assistant

The efficient hybrid kinematic-dynamic technique for finite-difference (FD) computation and storing the seismic wavefields is presented. At every time step of the FD scheme calculations are done only in the moving zone around the first-arrival wave front that is determined by the time window, consistent with the wave length. The considered approach allows one to store the wavefield at every time step, that could be useful for constructing the 'banana-doughnut' sensitivity kernels that are widely used in inversion methods.

Key words: seismic waves, numerical modelling, eikonal equation, wave-equation tomography.

Введение

Конечно-разностные схемы решения уравнений упругости широко используются для моделирования распространения сейсмических волн в сложных неоднородных средах. Одним из основных недостатков данных методов является большое время вычислений, поэтому ускорение расчетов является актуальной задачей. Как было показано в работах [1] и [2], если есть необходимость в нахождении волновых форм, которые соответствуют первым вступлениям, вычисления могут проводиться только в окрестности соответствующего фронта. Мы предлагаем новый способ реализации конечно-разностных вычислений в окрестности первых вступлений.

Основной особенностью способа является пересортировка используемых массивов в порядке возрастания времени вступления, который естественным образом следует из схемы fast marching [3]. Этот порядок позволяет существенно сократить время вычислений. Другим преимуществом предлагаемого подхода является существенное сокращение ресурсов памяти, необходимой для хранения волнового поля в окрестности первых вступлений с целью дальнейшего использования. В данной работе мы рассматриваем применение предлагаемого алгоритма моделирования динамики первых вступлений в рамках метода волновой томографии (wave-equation traveltime tomography), предложенного в работе [4].

Конечно-разностные вычисления в окрестности первых вступлений

Будем рассматривать двумерный случай. Пусть скоростная модель задана на регулярной прямоугольной сетке с узлами (i, j). На первом этапе моделирования мы решаем уравнение эйконала и получаем времена первый вступлений tf (i, j). Затем методом конечных-разностей решается волновое

уравнение (уравнения теории упругости), при этом вычисления на каждом временном шаге схемы проводятся только для точек сетки, удовлетворяющих условию:

nAt < tf (i, j) < nAt + T, (0.1)

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

На каждом шаге конечно-разностной схемы зона вычислений может быть определена путем проверки условия (0.1) во всех узлах пространственной сетки [2]. Для того чтобы избежать проверки условия (0.1) на каждом временном шаге мы предлагаем новую стратегию реализации «оконных» вычислений. Рассмотрим взаимно-однозначное соответствие:

(i,j) ^ к, (°.2)

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

Мы сортируем все массивы, участвующие в конечно-разностных вычислениях, по возрастанию индекса к . На каждом временном шаге n схемы «окно» представляет собой все точки от номера к = kbeglп(п) до к = kend(n)

, которые определяются используя список (0.3) и условие (0.1). Решение уравнений на текущем временном слое в активном «окне» хранится в последовательном порядке (по возрастанию индекса к) и таким, образом, может быть в том же порядке последовательно записано в каждый момент времени в конец одномерного массива. Заметим, что «геометрическая» позиция может быть легко восстановлена с помощью соответствия (0.2). Для конечно-разностного решения волновых уравнений мы применяем классическую схему второго порядка на сдвинутых сетках (схема Вирье) [5].

Волновая томография

Мы рассматриваем метод волновой томографии (wave-equation tomography) для уравнений акустики, впервые предложенный в [4]. Метод основан на минимизации функционала невязки:

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

где р(x, t; xs) - давление вычисленное конечно-разностной схемой для априорной скоростной модели c(х), а р’(х, t, *s) - давление для сопряженной задачи, вычисленное путем продолжения в обратном времени «псевдоневязок» (pseudoresidual, см. Lu and Schuster [4]), которые играют роль импульсов в источниках x .

Для нахождения градиента (0.5) необходимо вычислить «прямое» поле р (xr, t; xs) используя наш конечно-разностный алгоритм и сохранить поле в окрестностях первых вступлений. Затем необходимо найти временные «псевдоневязки» путем кросс-корреляции наблюдаемых и синтетических сейсмограмм. В завершении следует вычислить сопряженное поле р'(х, t, xs) путем продолжения «псевдоневязок» в «обратном времени» (reverse time) и одновременно найти «ядро чувствительности» (один интеграл для фиксированного источника s из правой части (0.5)). Наш метод позволяет

(0.4)

где x = (xs,zs) и xr = (x,z) - вектора координат источников и приемников; Ат

(0.5)

считывать записанное «прямое» волновое поле в «обратном» времени, используя индексы кЬет(п) и (п) на каждом слое и биекцию (0.2).

Примеры. Пример прямого конечно-разностного моделирования показан на рис. 1. Для простоты мы рассматриваем двумерное уравнение акустики. Синтетическая скоростная модель представляет собой градиентную среду с высокоскоростным включением (У=5000м/с). Скоростная модель, положение источника и изохроны времен первых вступлений показаны на рис. 1, а. Для инициализации использовался импульс Риккера с доминирующей частотой 30 Ш. Окрестности времен первых вступлений («окна») в разные моменты времени на рис. 1, б-в (вычисления проводятся только в белой полосе). Мгновенный снимок волнового поля в вычислительном окне в момент времени г = 0,58 с. показан на рис.1, г. На рис. 1, д показан мгновенный снимок полного волнового поля в тот же самый момент времени г = 0,58 с. Из сравнения рис. 1, г и рис. 1, д видно, что благодаря использованию предложенного в работе алгоритма посчиталась только прямая волна, прошедшая через высококонтрастное включение. На рис. 1, е показана относительная разность между полным волновым полем и полем, рассчитанным в окрестности первых вступлений, в вычислительном окне при г = 0,58 с., она оставляет тысячные доли процента. Для данного примера поле в окрестностях первых вступлений при помощи предлагаемого алгоритма вычисляется в 11 раз быстрее полного волнового поля.

Рис. 1. Пример конечно-разностных вычислений в окрестности первых вступлений: а) скоростная модель и изохроны, б), в) «окна» для вычислений, г) мгновенный снимок волнового поля в «окне» при г = 0,58 с., д) снимок полного волнового поля посчитанного

во всей области при X = 0,58 с., е) относительная разность между полным волновым полем и полем, рассчитанным в окрестности первых вступлений при X = 0,58 с

100 200 300 100 200 300 100 200 300

х,т х,т х,т

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

Пример применения метода волновой томографии на основе предлагаемого конечно-разностного алгоритма моделирования показан на рис.2. Синтетическая скоростная модель изображена слева, 40 источников и сорок приемников расположены на противоположных боковых сторонах области, в источнике использовался импульс Риккера с доминирующей частотой 60 Ш; в центральной части рис. 2 показаны «ядра чувствительности» для одного из источников и нескольких приемников; справа показан результат инверсии синтетических данных после применения одного шага градиентного спуска метода волновой томографии. Использование нашего метода «оконных» вычислений и сохранения поля позволяет уменьшить в 7 раз объем памяти, необходимый для хранения прямого поля. Мы получили 12 кратное ускорение для одной итерации волновой томографии по сравнению со стандартным подходом, в котором конечно-разностное моделирование проводится 3 раза во всей целевой области.

Заключение

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

Благодарности

Работа выполнена при частичной поддержке гранта президента № МК-2598.2014.5 и гранта компании BP для молодых ученых.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Bore, D.M. [1972] Finite-differencemethods for seismic wave propagation in heterogeneousmaterials. Methods in computational physics, 2, 1-36.

2. Crandall, M., Lions, P. and CENTER., W.U.M.M.R. [1983] Viscosity solutions of hamilton-jacobi equations..

3. Luo, Y. and Schuster, G.T. [1991]Wave-equation traveltime inversion. Geophysics, 56(5), 645-653.4.

4. Vidale, J. [1988] Finite-difference calculation of travel times. Bulletin of the Seismological Society of America, 78(6), 2062-2076.

5. Virieux, J. [1986] P-sv wave propagation in heterogeneous media: Velocity-stress finite-difference method. Geophysics, 51(4), 889-901.

© А. С. Сердюков, А. A. Дучков, А. А. Никитин, 2014