Доклады Российской академии наук. Математика, информатика, процессы управления, 2023, T. 512, № 1, стр. 58-64

1. ВВЕДЕНИЕ

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

Здесь $G \in {{\mathbb{R}}^{n}}$ – область в n-мерном вещественном пространстве, $u$ – набор переменных состояния системы, $A$ – оператор, задающий динамику среды, зависящий от конечного числа неизвестных функций-коэффициентов $\{ {{c}_{j}}(\vec {r})\} $. Известно конечное количество функций – наблюдений за средой $~{{f}_{m}}\left( {\vec {r}} \right)$, которые связаны с вектором переменных $u\left( {{{t}_{m}}} \right)$ в различные моменты времени ${{t}_{m}}$ с помощью оператора проекции $P$. Такая постановка является некорректной с математической точки зрения, задача решается численно путем оптимизации функционала, сравнивающего данные измерений с результатами моделирования динамики среды. Решение такой задачи инверсии обычно требует наличия тех или иных априорных предположений о структуре среды.

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

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

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

Методы машинного обучения были ранее рассмотрены рядом авторских коллективов для волнового уравнения и его различных приложений [4–6], в том числе конкретно для сейсмической инверсии [7, 8]. Вклад данной работы заключается в том, что рассмотрено использование нейронной сети, обученной на простых данных, для получения начального приближения для достаточно сложной модели Мармузи.

2. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ И ЧИСЛЕННЫЙ МЕТОД

В данной работе предполагается, что распространение импульса в исследуемой сплошной среде подчиняется системе уравнений акустики (2), что является корректным приближением для сред с сильным затуханием сдвиговых волн [9].

Здесь $\rho \left( {\vec {r}} \right)$ – плотность среды, $v\left( {\vec {r},t} \right)$ – вектор скорости, $p\left( {\vec {r},t} \right)$ – поле давления, $c\left( {\vec {r}} \right)$ – скорость звука в среде, $\alpha \left( {\vec {r}} \right)$ – коэффициент затухания. В данной работе предполагается, что затухание отсутствует. Для параметризации среды в данной работе используется величина ее акустического импеданса $~m = 1{\text{/}}{{c}^{2}}\left( {\vec {r}} \right)$, обратно пропорциональная квадрату скорости звука в среде. Классическая постановка задачи сейсморазведки подразумевает распространение импульса в полуплоскости, поэтому на верхней границе расчетной области ставится условие свободной границы, а на всех остальных – поглощающее граничное условие.

Для решения прямой задачи могут быть использованы различные численные методы. Так, в работе [10] для аналогичной задачи для модели акустики предложена компактная схема с использованием сеточно-характеристического численного метода. В данной работе используется явная конечно-разностная схема на регулярной сетке, обеспечивающая четвертый порядок аппроксимации по пространству и второй по времени, предложенная в [11]. На верхней границе расчетной области используется условие свободной границы. Отсутствие отражений волн от прочих границ расчетной области обеспечивает поглощающее условие типа C-PML [12].

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

Программная реализация решения прямой задачи в данной работе использует библиотеки PyTorch [13] и TorchFWI [14, 15] и поддерживает расчет на видеокартах на базе технологии CUDA.

В качестве метода численной оптимизации в данной работе используется L-BFGS-B [16] – квазиньютоновский алгоритм, использующий вместо матрицы Гессе ее двухточечную аппроксимацию, что значительно повышает его привлекательность в применении к крупномасштабной оптимизации. Взятие градиента функционала ошибки по параметрам среды реализуется с помощью метода сопряженных переменных состояния [17], получить который можно методом множителей Лагранжа, рассматривая целевую задачу как задачу условной минимизации с ограничениями в виде системы уравнений в частных производных. Применение описанного метода к системе уравнений акустики дает вид (3) для градиента функционала ошибки J по параметрам системы $m~$ в точке $\vec {r}~$ расчетной области.

Здесь $T$ – итоговое время расчета, а сумма ведется по конечному числу $s$ источников возмущения. Величина ${{q}_{s}}$ определяется из сопряженной системы уравнений, которая, в случае использования среднеквадратичного отклонения в качестве функционала потерь, имеет вид (4).

Здесь ${{P}_{{s,r}}}$ – оператор проекции, характеризующий снятие измерений, a ${{d}_{{s,r}}}$ – наблюдаемая сейсмограмма.

Для предсказания распределений параметров в среде по входному сигналу используется полносверточная нейронная сеть архитектуры UNet. Используется классическая архитектура UNet [18], содержащая 30 миллионов параметров. Ранее сеть такой архитектуры уже была успешно использована для более простых геологических задач, в работе [19] с ее помощью решалась задача локализации трещиноватых зон по сейсмическому отклику от них. Следует отметить, что уменьшение размерности сети без значимого снижения качества предсказаний представляется реалистичной задачей, но на данном этапе она не рассматривалась.

3. РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

Для обучения и тестирования нейронной сети используется набор данных из 20 тысяч пар “распределение свойств в среде – сейсмограмма”, полученный методами численного моделирования.

Геометрический размер области каждого образца составляет 9192 м по горизонтали и 3192 м по вертикали. Для каждого образца сначала генерируется некоторое распределение свойств в среде. Генерация происходит процедурным образом, алгоритм генерации обеспечивает получение случайных слоистых структур. Каждое такое распределение (изображение) содержит 384 узла сетки (пиксела) по горизонтали и 134 по вертикали. Значение в каждом узле (пикселе) – значение скорости звука в данной точке пространства. Для полученного распределения скоростей звука в среде численно решается прямая задача и получается расчетная сейсмограмма.

При решении прямой задачи шаг регулярной расчетной сетки по пространству равен 5 м. Количество датчиков на верхней границе расчетной области, выполняющих запись сейсмограммы – 384 штуки. Сканирующий импульс задается в центре верхней границы расчетной области.

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

Рис. 1.

Примеры сгенерированных случайных слоистых гетерогенных сред (верхний ряд) и сейсмограмм для них (нижний ряд).

Полученный набор из 20 000 пар разделяется на обучающий, проверочный и тестовый наборы в соотношении: 16 000, 3000, 1000 образцов.

На полученных данных обучается сверточная нейронная сеть. Обучение нейронной сети велось с помощью алгоритма оптимизации Adam с постоянным шагом 0.001. В качестве критерия остановки оптимизатора выступала неизменность функции потерь на тестовой выборке в пределах 5% от ее текущего значения в течение 5 эпох обучения, после чего производился перезапуск алгоритма с уменьшенным в 10 раз шагом обучения. Для обучения использовались GPU NVIDIA Tesla K80 12 Gb. В конце обучения точность предсказания, выраженная как попиксельный MSE, составляет 0.099 для проверочного набора и 0.113 для тестового набора.

На рис. 2 приведены примеры предсказаний обученной сети, которые в дальнейшем использовались как начальные точки для работы алгоритма градиентной оптимизации. Справа показаны реальные конфигурации сред, а слева – предсказания их структуры, сделанные нейронной сетью на основании сейсмограммы. По осям отложены геометрические размеры в метрах, цветом показана скорость звука в данной точке области, выраженная в м/с.

Рис. 2.

Примеры предсказаний нейронной сети (слева) и реальных конфигураций сред (справа).

Чистое время предсказания для одного образца составляло 0.4 с. С учетом загрузки данных с диска в память и выполнения технических процедур полное время получения предсказания составляло 1.2 с для одного образца.

Для тестирования совместной работы нейронной сети и градиентной оптимизации была использована геологическая модель Мармузи [20], которая широко используется при отработке алгоритмов инверсии. На рис. 3 показано истинное распределение параметров среды для данной модели. По осям отложены геометрические размеры в метрах, цветом показана скорость звука в данной точки области, выраженная в м/с.

Рис. 3.

Распределение скоростей в среде для модели Мармузи.

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

Для данной сейсмограммы с помощью обученной ранее нейронной сети было получено начальное приближение структуры среды. Далее на базе этого начального приближения была выполнена  последовательность итераций методом L-BFGS-B. Отдельно следует отметить, что данные Мармузи никаким образом не использовались при обучении нейронной сети. То есть начальное приближение было получено нейронной сетью, которая обучалась исключительно на синтетических данных для случайных слоистых структур.

На рис. 4 приведен натуральный логарифм функционала потерь при проведении 50 итераций градиентной оптимизации алгоритма L-BFGS-B с различными видами начального приближения. Рассматривается начальное приближение, выданное нейронной сетью, и простое линейное распределение скорости по глубине. Следует отметить, что данный функционал потерь специфичен для метода градиентной оптимизации и никак не связан с функционалом потерь для самой нейронной сети, предоставившей начальное приближение.

Рис. 4.

История минимизации функционала при различных типах инициализации при решении задачи FWI для модели Мармузи с помощью алгоритма L-BFGS-B.

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

Из представленных на рис. 4 результатов численных  экспериментов видно, что скорость сходимости метода градиентной оптимизации L-BFGS-B при использовании предложенного метода получения начального приближения не изменяется по сравнению с его скоростью сходимости при использовании линейной интерполяции в качестве начального приближения. Однако абсолютное значение функционала потерь при использовании предложенного метода меньше на каждой итерации. Таким образом, для достижения одинакового значения функционала потерь при использовании предложенной инициализации потребуется меньшее количество итераций градиентного метода.

На рис. 5 показана траектория оптимизатора в пространстве параметров для первых 50 итераций. По осям отложены геометрические размеры в метрах, цветом показана скорость звука в данной точке области, выраженная в м/с. Результат нулевой итерации – предсказание нейронной сети. На качественном уровне видно, что сеть определяет принципиальную локализацию крупных структур. Обеспечить тонкое разрешение структуры слоев и разломов нейронная сеть оказывается не в состоянии, что достаточно ожидаемо с учетом ее обучения на исключительно синтетических геометриях. При этом предсказание нейронной сети оказывается удачным начальным приближением, что иллюстрируют последующие шаги алгоритма, который уверенно сходится к реальной структуре среды.

Рис. 5.

Траектория L-BFGS-B при решении задачи FWI для модели Мармузи при инициализации предсказанием нейронной сети.

4. ЗАКЛЮЧЕНИЕ

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

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

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