Регуляризация обратной динамической задачи сейсмики посредством усечения SVD Текст научной статьи по специальности «Математика»

УДК 550.34.01

DOI: 10.18303/2618-981X-2018-3-30-35

РЕГУЛЯРИЗАЦИЯ ОБРАТНОЙ ДИНАМИЧЕСКОЙ ЗАДАЧИ СЕЙСМИКИ ПОСРЕДСТВОМ УСЕЧЕНИЯ SVD

Алексей Александрович Василенко

Новосибирский национальный исследовательский государственный университет, 630090, Россия, г. Новосибирск, ул. Пирогова, 2, студент, e-mail: alexvas1996@ya.ru

Владимир Альбертович Чеверда

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

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

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

REGULARIZATION OF THE INVERSE DYNAMIC SEISMIC PROBLEM BY TRUNCATED SVD

Alexey A. Vasilenko

Novosibirsk National Research State University, 2, Pirogova St., Novosibirsk, 630073, Russia, Student, e-mail: alexvas1996@ya.ru

Vladimir A. Cheverda

Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, 3, Prospect Аkademik Koptyug St., Novosibirsk, 630090, Russia, D. Sc., Professor, phone: (383)333-00-54, e-mail: CheverdaVA@ipgg.sbras.ru

The inverse dynamic problem of seismic consists in recovering the velocity model based on the solution of the wave equation known at some acquisition. In the process of solving this problem a system of linear algebraic equations arises. As a rule, its matrix is ill-conditioned and needs SVD-analysis to perform correct pseudoinversion. This paper presents and discusses the results of numerical experiments.

Key words: geophysics, computational mathematics, inverse problem, singular value decomposition.

Введение

Общепризнано, что одним из самых общих методов извлечения количественной информации об упругих параметрах среды из сейсмических данных яв-

ляется обращение полных волновых полей (синонимы: решение обратной динамической задачи сейсмики (ОДЗ), full waveform inversion (FWI) и сейсмическая инверсия). Как правило, ОДЗ формулируется как нелинейная задача минимизации целевого функционала, характеризующего (в некоторой норме) невязки сейсмограмм, зарегистрированных в ходе полевых работ, а также рассчитанных синтетических сейсмограмм [1, 2].

Однако прямое применение FWI к решению ОДЗ устойчиво восстанавливает только резко осциллирующую компоненту решения и неспособно восстановить гладкую скоростную модель. Здесь мы изучаем эту проблему на основе анализа сингулярного разложения гессиана оператора обратной задачи.

Постановка задачи и метод решения

Волновой процесс на отрезке [0, L] описывается следующей начально-краевой задачей:

дй dz

1 д2й д2й г, , ч с2(z) dt2 dz2

= 0,

z=0

^ дй 1 дйЛ

--1---

ydz сда dt j

= 0,й|,=0=0, ut |f=0=0,

z=L

где с(2) - скорость распространения волны, zs - координата источника, /(/) -функция, определяющая характер импульса в точке источника (в нашем случае возьмем импульс Рикера с доминирующей частотой 25 Гц), 8 (2 - ) - дельта-функция Дирака.

При переходе в частотную область будет справедливо уравнение Гельм-гольца:

d 2u

du dz

ю

dz 2 c2(z)

u = - F (ю)8( z - zs),

= 0,

z=0

du ю --h i— u

dz

c

X j

= 0.

z=L

Обратная динамическая задача сейсмики состоит в определении функции c(z) по заданному решению ф(ю) данной краевой задачи в некоторой точке zr - приемнике [1, 3].

Для ее решения используется метод обращения полного волнового поля (англ. Full Waveform Inversion или FWI). Этот метод заключается в минимизации разности наблюдаемых и смоделированных (синтетических) данных в среднеквадратичном смысле. Введем оператор B[c], который является реше-

нием задачи в точке 2Г при некоторой скорости с(2). Таким образом, нужно минимизировать следующий функционал:

Ф[с] = В [с] - ф1 2 = (В[с] - ф,Б[с] - ф),

где норма берется в гильбертовом пространстве .

Стандартным подходом здесь служит организация итерационного процесса отыскания точки минимума. В рамках нашей задачи применяются две достаточно известные реализации локальных градиентных методов: метод сопряженных градиентов и модифицированный метод Ньютона [4, 5]. Обе техники предполагают численное вычисление градиента и гессиана исходного нелинейного функционала.

Модифицированный метод Ньютона сводится к следующему итерационному процессу:

НМ ск+1 - ск) = -УФ[с^ ].

Каждое следующее приближение ищем в виде ск+1 = ск + 5ск. Пренебрегая слагаемыми второго порядка, для Ф[ ск+1] справедливо

Ф[ск + 5ск ] = Ф[ск ] + (УФ[ск ], Ъск) +1 (Н [ск ] (5с^, Ъск), Ф[ск + 5ск ] = Ф[со] + 2Яе (ОВ*[сЬ^В[со] - ф, 5с) + (ВВ^рВ^^ 5с), 5с).

Отсюда получаем выражения для градиента и гессиана функционала, соответственно:

УФ[ск ] = 2ЯеВВ*[ск ^В[ск ] -ф,

*

Н[ск] = 2БВ [ск]БВ[ск],

где БВ[с] - производная Фреше оператора В[с], которая, в свою очередь, является интегральным оператором.

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

Результаты численных экспериментов

*

Пусть матрица А соответствует оператору И[с0] и А = Ц£У - сингулярное разложение матрицы А, где и, У - матрицы, состоящие из левых (ип) и правых (уп) сингулярных векторов, соответственно, £ - матрица с неотрицательными элементами, у которой элементы, лежащие на главной диагонали -сингулярные числа, а все остальные элементы - нулевые. Предположим, что а! > <з2 >...> а N, т. е. сингулярные числа в матрице £ расположены по убыванию.

Числом обусловленности матрицы А называется величина

и(А) = 8ТО||_АХШ1 и атах(А) _ а!(А) Ц( ) ||1 А11 II XI ] атп(А) ам(А)'

Так как рассматриваемый оператор компактен, то число ц(А) достаточно велико. Характер убывания сингулярных чисел виден на рис. 1 [2, 6]. Введем величину

Ц, (А) = •

аг ( А)

Рис. 1. Сингулярные числа гессиана в логарифмическом масштабе

На практике номер г выбирают так, чтобы величина цг (А) имела порядок

10 ■ 104. При таком выборе первые г правых сингулярных векторов образуют базис устойчивого подпространства Хг в пространстве моделей X [3, 7].

Таким образом, для матрицы А можно построить г - псевдообратную матрицу ] такую, что если Ах = у, то

Х\

[г ]

=4] уV

п=1 ап

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

Рис. 2. Скоростные модели

а) б)

Рис. 3. Проекции скоростей на устойчивое подпространство Хг

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

Заключение

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

Данная работа выполнена при поддержке РНФ, проект 17-17-01128.

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

1. Gadylshin K., Cheverda V. Reconstruction of a depth velocity model by full waveform inversion // Doklady Akademii Nauk. - 2017. - Vol. 476, No. 6. - P. 693-697.

2. Silvestrov I., Tcheverda V. SVD analysis in application to full waveform inversion of mul-ticomponent seismic data // J. Phys. : Conf. Ser. 290 012014. - 2011.

3. Cheverda V., Clement F., Khaidukov V., Kostin V. Linearized inversion of data of multi-offset data for vertically inhomogeneous background // J. Inv. Ill-Posed Problems. - 1998. - Vol. 6, No. 5. - P. 453-484.

4. Васильев Ф. П. Численные методы решения экстремальных задач. - 2-е изд., переработ. и доп. - М. : Наука, 1988. - 552 с.

5. Канторович Л. В., Акилов Г. П. Функциональный анализ. - М. : Наука, 1984. - 752 с.

6. Chavent G., Papanicolaou G., Sacks P., Symes W. Inverse Problems in Wave Propagation. Springer Science & Business Media, 1997. - Р. 277-294.

7. Костин В. И., Чеверда В. А. R-решения уравнения первого рода с компактным оператором в гильбертовых пространствах: существование и устойчивость // Доклады РАН. -1997. - Т. 355. - № 3.

REFERENCES

1. Gadylshin K., Cheverda V. Reconstruction of a depth velocity model by full waveform inversion // Doklady Akademii Nauk. - 2017. - Vol. 476, No. 6. - P. 693-697.

2. Silvestrov I., Tcheverda V. SVD analysis in application to full waveform inversion of mul-ticomponent seismic data // J. Phys. : Conf. Ser. 290 012014. - 2011.

3. Cheverda V., Clement F., Khaidukov V., Kostin V. Linearized inversion of data of multi-offset data for vertically inhomogeneous background // J. Inv. Ill-Posed Problems. - 1998. - Vol. 6, No. 5. - P. 453-484.

4. Vasil'ev F. P. Chislennye metody resheniya ekstremal'nyh zadach. - 2-e izd., pererabot. i dop. - M. : Nauka, 1988. - 552 s.

5. Kantorovich L. V., Akilov G. P. Funkcional'nyj analiz. - M. : Nauka, 1984. - 752 s.

6. Shavent G., Papanicolaou G., Sacks P., Symes W. Inverse Problems in Wave Propagation. Springer Science & Business Media, 1997. - Р. 277-294.

7. Kostin V. I., Cheverda V. A. R-resheniya uravneniya pervogo roda s kompaktnym operatorom v gil'bertovyh prostranstvah: sushchestvovanie i ustojchivost' // Doklady RAN. -1997. - T. 355. - № 3.

© А. А. Василенко, В. А. Чеверда, 2018