ПРИМЕНЕНИЕ АСИМПТОТИЧЕСКОГО РЕШЕНИЯ УРАВНЕНИЯ ГЕЛЬМГОЛЬЦА В ПОЛНОВОЛНОВОМ ОБРАЩЕНИИ Текст научной статьи по специальности «Математика»

УДК 528.71

DOI 10.33764/2618-981X-2022-2-2-269-274

Применение асимптотического решения уравнения Гельмгольца

в полноволновом обращении

К. Г. Гадыльшин1, Д. А. Неклюдов1, M. И. Протасов1 *

1 Институт нетфегазовой геологии и геофизики СО РАН, г. Новосибирск,

Российская Федерация

* e-mail: protasovmi@ipgg.sbras.ru

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

Ключевые слова: полноволновое обращение, асимптотическое решение, уравнение Гельм-гольца

Application of asymptotic solution of Helmholtz equation in Full Waveform Inversion

K. G. Gadylshin1, D. A. Neklyudov1, M. I. Protasov1 *

1 Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk, Russian Federation

* e-mail: protasovmi@ipgg.sbras.ru

Abstract. This study considers the Full Waveform Inversion (FWI) method based on the asymptotic solution of the Helmholtz equation. We provide frequency-dependent ray tracing to get the wavefield used to compute the FWI gradient and calculate the modelled data. With a comparable quality of the inverse problem solution as applied to the standard finite-difference approach, the speed of calculations in the asymptotic method is an order of magnitude higher. A series of numerical experiments demonstrate the approach's effectiveness in reconstructing the macro velocity structure of complex media for low frequencies.

Keywords: full waveform inversion, asymptotic solution, Helmholtz equation

Введение

Требуемые вычислительные ресурсы для моделирования сейсмических данных для типичной системы наблюдения можно оценить примерно в 108 ядро часов (при расчете на центральных процессорах). Эти оценки верны для изотропной упругой среды. При усложнении модели значительно возрастет потребность в вычислительных ресурсах. Поэтому разработка методов численного моделирования волновых сейсмических полей и ускорение алгоритмов такого моделирования является чрезвычайно актуальной задачей.

Асимптотические методы работают намного быстрее, чем методы, основанные на конечных разностях или конечных элементах, которые используются в большин-

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

Метод

Пусть функция и(х, z; м) (полное волновое поле в частотно-временной области) удовлетворяет уравнению Гельмгольца:

(А + о)2т)и = —/(m)S(x — xs)S(z — zs), (1)

где f(<u) - спектральная характеристика сигнала в источнике, ш - временная частота, (xs, Zs) - координата источника, т = с~2(х, z) - квадрат медленности, и с(х, z) - скорость распространения волны. Введем оператор, вычисляющий волновое поле от одного точечного источника для фиксированной временной частоты в точках, соответствующих расположению приемников. Таким образом получается следующий оператор прямого моделирования:

F: M^D, (2)

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

dobs = F(mtrue), (3)

здесь dobs - наблюдённые данные, mtrue - истинная скоростная модель.

При полноволновом обращении применяется нелинейный метод наименьших квадратов к уравнению (3). Общая формулировка полноволнового обращения состоит в том, чтобы найти точку минимума функционала невязки, характеризующего среднеквадратичное отклонение наблюдённых данных от рассчитанных для текущей скоростной модели [3]:

m* = arg(min\\F(m)-d\\l). (4)

\тем /

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

тк+1 =mk + ßkSk, S0 = V0, (5)

с _ (У^Уг-Уц-!),, ,,,

где тк - модель на к-ой итерации. Градиент рассчитывается следующим образом:

Ук= Яе{Х5Г.ы 2С(х5. г5; х. г; ш; mk)G(x. г; хг. 2Г; ш; тк) 6с15 г} , (7)

где Яе - вещественная часть, С(х5. х, х; м; тк) - функция Грина, рассчитанная в модели тк для источника, расположенного в точке (х5. С(х, х; хг. ; м; тк) - функция Грина, рассчитанная в модели тк для точки приемника (хг. и точки источника (х. г), - невязка данных, зарегистрированная в точке приемника (хг. и сгенерированная в точке источника (х5. Таким образом, расчет градиента сводится к расчету функций Грина для всех положений источников и приемников, а значит, большой интерес представляет ускорение расчетов именно этой части. В этом подходе используется асимптотический метод, основанный на частотно-зависимом трассировании лучей (Ьошах, 1994). Функции Грина рассчитываются на основе времени прохождения, амплитуды и геометрического расхождения в результате частотно-зависимого трассирования лучей.

Численные эксперименты

Для тестирования предложенного подхода асимптотического полноволнового обращения используется модель Магшош1, которая является эталоном для тестирования алгоритмов обращения. Целевые размеры исследуемого участка составляют 9200 м по ширине и 3000 м по глубине (первые 500 м скоростной модели - слой воды, который считается известным). Расчетная сетка состоит из 921 точки по горизонтали и 351 точки по вертикали (шаг 10 м). Система наблюдения состоит из 91 источника и 459 сейсмоприемников, расположенных на поверхности воды с шагом 100 м и 20 м соответственно. В качестве начального приближения используется сглаженное по Гауссу распределение истинной скорости (рис. 1). Наши численные эксперименты направлены на определение границ применимости асимптотического метода Б^. Мы проводим численное сравнение градиентов на основе асимптотического решения и конечно-разностного решения в исходной модели. Затем мы приводим эксперименты, сравнивающие результаты инверсии на основе асимптотических и конечно-разностных решений.

Прежде чем проводить численные эксперименты для полноволнового обращения, мы рассчитываем и сравниваем градиенты. Градиенты вычисляются в исходной модели (Рисунок 1) для классической постановки и для асимптотического подхода. Полученные градиенты для частоты 3 Гц представлены на Рисунке 2. Как видно, градиенты визуально совпадают, что объясняется тем, что исходная модель гладкая, волновая картина в ней достаточно простая, и, следовательно, хорошо работает асимптотическое решение, основанное на лучевом приближении.

Рис. 1. Истинное распределение скоростей в среде (модель Marmousi) и начальная модель для полноволнового обращения

Рис. 2. Градиент классического полноволнового обращения для 3 Гц и градиент асимптотический обращения для 3 Гц, рассчитанные в стартовой модели

Далее с помощью следующего численного эксперимента мы демонстрируем возможность применения асимптотического подхода к полноволновому обращению для уточнения макроскоростной модели на частоте 3 Гц. Результаты обращения представлены на Рисунке 3. Как видно, модель, восстановленная с использованием асимптотического решения уравнения Гельмгольца, в целом аналогична результату классической инверсии. На основании этого сравнения можно сделать вывод, что даже для такой сложной модели, как МагшоиБ^ уточнение макроскоростной модели на низких частотах может быть достигнуто за счет использования асимптотического приближения.

Заключение

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

Рис. 3. Результат классического полноволнового обращения для частоты 3 Гц (вверху) и разработанного асимптотического метода обращения (внизу)

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

Данное исследование выполнено при финансовой поддержке РФФИ и ЧНФ в рамках научного проекта №20-55-26003. Acknowledgments: The reported study was funded by RFBR and GACR, project number 20-55-26003.

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

1. Lomax A. The wavelength-smoothing method for approximating broad-band wave propagation through complicated velocity structures: Geophys. J. Int. - 1994. - Vol.117. - P. 313-334.

2. Pratt R.G. Seismic waveform inversion in the frequency domain, part 1: Theory and verification in a physical scale model: Geophysics. - 1999. - Vol. 64. - P. 888-901.

3. Tarantola A. Inverse problem theory: Methods for data fitting and model parameter estimation. Elsevier Science Publ. Co., Inc. - 1987.

4. Gadylshin K. G., Cheverda V. A. Reconstruction of a depth velocity model by full-waveform inversion: Doklady Earth Sciences. - 2017. - Vol.476. - P. 1233-1237.

© К. Г. Гадыльшин, Д. А. Неклюдов, М. И. Протасов, 2022