РЕШЕНИЕ ОБРАТНОЙ ДИНАМИЧЕСКОЙ ЗАДАЧИ СЕЙCМИКИ ДЛЯ СКАЛЯРНОГО ВОЛНОВОГО УРАВНЕНИЯ НА ОСНОВЕ ГЛУБОКОЙ СВЁРТОЧНОЙ НЕЙРОННОЙ СЕТИ Текст научной статьи по специальности «Компьютерные и информационные науки»

УДК 550. 3

DOI 10.33764/2618-981X-2022-2-2-292-304

Решение обратной динамической задачи сейcмики для скалярного волнового уравнения на основе глубокой свёрточной нейронной сети

Д. С. Братчиков1, К. Г. Гадыльшин1 * 1 Институт нефтегазовой геологии и геофизики им. А. А. Трофимука СО РАН, г. Новосибирск, Российская Федерация

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

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

Ключевые слова: сейсмический мониторинг, обращение полного волнового поля, глубокое обучение, свёрточные нейронные сети

Solution of an inverse dynamical seismic problem on the base of a deep convolutional neural network

D. S. Bratchikiov1 *, K. G. Gadylshin1

1 Trofimuk Institute of Petroleum Geology and Geophysics SB RAS, Novosibirsk,

Russian Federation

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

Abstract. The paper presents an original approach to solving the inverse dynamical seismic problem for the scalar wave equation in seismic monitoring. Instead of using the expensive full waveform inversion method, it is proposed to approximate the inverse problem operator, which translates the difference in the data into the models' perturbation, using a deep convolutional neural network of the Unet type architecture. The operability of this approach was demonstrated on a simple test example.

Keywords: seismic monitoring, full waveform inversion, deep learning, convolutional neural networks

Введение

Для оптимальной разработки месторождений углеводородов необходим быстрый и эффективный сейсмический мониторинг, обеспечивающий получение информации о совместном движении нефти и вытесняющего её агента. За последнее время одним из самых достоверных и высокоразрешённых способов построения глубинной скоростной модели Земли стал метод обращения полного волнового поля (FWI, от англ. Full Waveform Inversion). Однако существующие численные методы обращения полного волнового поля имеют высокую стои-

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

В данной работе для устранения этих недостатков предлагается оригинальный подход к решению обратной динамической задачи сейсмики для скалярного волнового уравнения с привлечением элементов машинного обучения. Метод основан на обучении глубокой свёрточной нейронной сети (CNN, от англ. Convolutional Neural Network) архитектуры типа Unet [1] для аппроксимации нелинейного оператора обратной задачи, задающего связь между рассчитанными волновыми сейсмограммами и скоростной моделью геологической среды. Для обучения CNN используется синтетическая обучающая выборка двумерных данных: изменение скоростных моделей и соответствующее им изменение сейсмических данных, вычисленных на основе конечно-разностных схем для волнового уравнения. Такой подход позволяет избежать моделирование оператора прямой задачи. Основным его преимуществом является экономия вычислительных ресурсов на дорогостоящих процессах обработки сейсмических данных на каждой стадии вытеснения углеводородов. В предложенном способе такие затраты производятся на самых первых этапах для создания репрезентативной обучающей выборки и обучения нейронной сети. На всех последующих этапах обработка заменяется использованием обученной нейронной сети, которая по изменчивости сейсмических данных на разных этапах разработки месторождения будет предсказывать изменения петрофизической модели резервуара. Серия проведенных численных экспериментов позволяет убедиться в успешности такого подхода.

Методы и материалы

1. Прямая задача моделирования волнового поля

1.1 Постановка задачи

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

1 кг /,ч

V2 W = А*,zU + f, (1)

где f - источник, U - волновое поле, V - скорость распространения сейсмических волн.

В качестве источника возьмем вейвлет Рикера:

/(х, г)ё, ^) (1-2 V - (2)

где ё(х, ) - дельта-функция Дирака, X,- координаты источника, г0 - сдвиг,

V - частота.

1.2 Конечно-разностная аппроксимация

Рассмотрим ограниченную область П = {(х, г) 6 М2:0 <х<Л, 0 <г< В]. Дискретизируем задачу на равномерной прямоугольной сетке. Для численного решения поставленной задачи воспользуемся явной конечно-разностной схемой. Аппроксимируем вторую частную производную по времени со вторым порядком точности, а вторые частные производные по пространственным координатам с четвертым порядком точности:

.и-1

а 2u _ и-1 - iuh+^+о (д, 2),

a 2t

_ Ui+2, j

At2 + 16u

i+1, j

-30un, +16un, .

_w_i 1, J

u

i-2, j

12 Ax2

+ O (Ax4),

a 2u__

a2 x

a 2u _ +2 + iKj+1- 3°<j+-1- Kj-2 +O (Az 4

(3)

a2 z

12 Az2

Подставляем формулы в уравнение (1) и явно выражаем решение на п +1 шаге. Получаем численную формулу для решения прямой задачи:

î

un++1 _ V, j

1

ut,, ,.+—u" ,.-—un

+v

t, j

1

4

3 '

4

__u _____

v 12 i+2,j 3 i+1,j 3 i-1,j 12 i-2,j j

-—u:;..,„ +—un

4 3'

4

1

Va ^

vAx j

1

+- u", i -—u"

+2

1 --vt

4 t,J

12 i,J+2 3 i,J+1 3 i,J-1 12 i,J-2

u;: j - un-1+j j At

Va ^

vAz j

+

+

5 ,2

r At2 At2 Ax2 +Az2

(4)

1.3 Граничные поглощающие условия

Для численного моделирования сейсмических волновых полей приходится рассматривать усеченную модель Земли. Чтобы устранить граничные отражения, вызванные этим усечением, применяются поглощающие граничные условия. Для данной задачи мы используем effective absorbing layer (EAL) [2], который является split perfectly matched layer (SPML) c упрощенным демпфирующим членом. Преимуществами этого граничного условия является его точность, экономичность и легкая реализуемость для высокопроизводительных процессов. Так же применение этого граничного условия для FWI приводит к сокращению вычислительных затрат и затрат на память.

Реализуется вышеописанное поглощающее граничное условие путем добавления в волновое уравнение (1) демпфирующего члена:

^+2D |U+D2u = V2 (и+f),

(5)

где Д х, 2) = ^

'1Л

3v

/ /\2

21

1_

- демпфирующая функция, Ь - ширина поглощаю-

щего слоя, / (х, г) - расстояние от точки в поглощающем слое до внутренней границы области, V (х, г) - значение скоростной модели в соответствующей точке в

поглощающем слое, О - отражающий коэффициент БРМЬ, который равен 0.001 в нашей задаче.

2. Обратная задача обращения полного волнового поля

Прямую задачу в уравнении (1) можно представить в виде нелинейного преобразования:

й = Т [ т ],

(6)

1

где Т:М^D - оператор прямой задачи из пространства моделей, т=— ЕМ -

медленность, в пространство данных, й=и(хГ ,t)eD - сейсмограммы, X = (X, 2Г) - координаты приемников.

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

распространения волнового поля в акустической среде, таких как обращение полного волнового поля, она решается с помощью минимизации функционала:

mmmJ[т] [т] = 2\\йоЫ — Т[т]У2,

(7)

Т 2

где \\й \\2 (7) &, йГ(7) = и£ (хГ,7), £ - индекс источника.

Г

Сначала происходит моделирование прямого оператора Т с помощью уравнения (5). Затем восстанавливаем скоростную модель среды по наблюдаемым данным путем минимизации целевого функционала.

Минимизация функционала (7) осуществляется с помощью градиентного спуска:

^ т(к) 1. (8)

т(к+1) = т(к) — о

SJ

5т

Обозначим ^ [т ]. Тогда из [3] для нахождения градиента известно, что: 5т

5£ 5ш

[ш] = ^*(т[ш] - а).

(9)

Таким образом для решения задачи (7) нужно найти сопряженный оператор производной Фреше прямого оператора. В [3] так же приведен вывод следующего утверждения:

д 2и0

()(х) = -£ |о д(х,х(х,х)ах

V а.

(10)

где и0, - волновое поле, соответствующее источнику , и стартовой модели ш0 :

д2 л

Ш дХ^-Д

ио, =/>

(11)

а сопряженное поле д, соответствующее источнику 5, находится из системы уравнений, которая решается обратно во времени:

т,

о

л- -л

V дХ У

Ч\<=Т дх

Ч , (x, Х) = ^ехМ (x, Х),

= о,

(12)

X=Т

где правая часть аех1з(х, X) = ^ аг,, (X)£(х - хг), Т - время записи наблюдаемых

г

сейсмограмм 0.,, (().

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

3. Свёрточная нейронная сеть

3.1. Постановка задачи

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

тт т-1 Ь (шхгие, т ~1[ааЬ* ])

где Ь - функция потерь, используемая при обучении.

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

Мы сконцентрируемся на обращении полного волнового поля на основе методов глубокого обучения. Будем аппроксимировать оператор Т-1 с помощью свёрточной нейронной сети. Согласно теореме в [4] любую непрерывную функцию нескольких переменных можно аппроксимировать сверточной нейронной сетью с любой наперёд заданной точностью. Поэтому если выборка будет достаточно репрезентативной, то отображение сейсмических данных в скоростную модель может быть выучено с довольно высокой точностью. В силу того, что мы предполагаем дифференцируемость оператора прямой задачи, то непрерывность оператора обратной задачи отсюда следует автоматически. После обучения мы применяем восстановленный оператор Т-1 на новых сейсмических данных Лпец! для прогнозирования скоростной модели среды тряеЛ.

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

Стоит заметить, что в нашей работе мы делаем переход к возмущениям в моделях и в данных. Считаем, что в момент времени ¿1 построена модель резервуара, то есть считаем Щ, Л нам известными. В какой-то момент времени ¿2 скоростная модель резервуара Щ2 и соответствующие им данные Л немного изменились. Вместо поиска Щ2 по данным Л>в постановке (6) будем искать возмущение моделей Ат = Щ2 — по невязке в данных АЛ = — Л. Тогда задача обучения (13) заменится следующей задачей обучения:

Ш1П Т—11 (А т, Т—1[А Л ]) (14)

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

3.2. Архитектура сверточной нейронной сети

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

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

Рис. 1. Архитектура свёрточной нейронной сети

Каждый свёрточный блок состоит из оператора свертки Conv, пакетной нормализации BN [5] и функции активации ReLU [6]. Размер свёрточного ядра 3 х 3, шаг скольжения ядра свертки равен 1, для сохранения размеров карты признаков после каждой свертки используется padding равный 1. Для понижения размерности карты признаков в кодировщике используется оператор агрегации max pooling. Для увеличения размерности карты признаков в декодировщике используется оператор max unpopling. В конце нейронной сети используется свертка с ядром размера 1х1 для уменьшения числа каналов и функция активации sigmoid для приведения значений выходной карты признаков к отрезку [0,1]. Ниже привден явный вид всех упомянутых функций:

x7+1

ReLU (BN | Conv (x1 ))),

Conv(x)(/, j) =xxx km,n,c ' X

( s —1)х/ + m,( s —1)х j+n ,c m n c

f Л

BN^(X/, j ,c ) = r

x/,j ,c Mb

£

V Ув + * у

+ A (15)

ReLU(x) = max(0, x), 1

sigmoid( x)

Здесь входные данные на свёрточный блок X и обучаемое ядро оператора свертки К - трехмерные тензоры с первыми двумя размерностями /, у и т, п отвечающими за пространственную координату, а третьей С за количество каналов, £ - шаг с которым ядро свертки пробегает по входным данным, у и р - два

тренировочных параметра, и - среднее значение и дисперсия, вычисленные для каждой карты признаков X по пакетам данных, Б - пренебрежимо малая константа, введённая для численной устойчивости метода.

3.3. Обучение сверточной нейронной сети

Таким образом с помощью свёрточной нейронной сети, архитектура которой описана выше, мы ищем оператора обратной задачи Т 1 в виде отображения:

F 1 ~ R = fn (•■• f2 ((1 ( wi) w2 , wn )

=

/i (d, wi), i = 1 (16)

i = n.

f (, w), i >1

где fi - отвечает за один свёрточный блок, d - входные сейсмические данные, W - веса сети, n - количество свёрточных блоков. В нашем случае количество блоков равно 14 и определяется из соображений вычислительной производительности и качества восстановления оператора.

Для обучения свёрточной нейронной сети необходимо выбрать функцию потерь L и оптимизатор, с помощью которого минимизируется функция потерь в (14). Для нашей задачи восстановления отображения мы выбрали наиболее популярную функцию среднюю квадратичную ошибку (17) и один из лучших оптимизаторов за последнее время с адаптивной скоростью обучения Adam [7].

1 N , Ч2

L = NХ(( - ^) , (17)

N i=1

где mfue - истинная модель, mfred = R(d°bs) - предсказанная модель, N - размер тренировочного набора данных.

Обучение происходит обновлением весов в (16) с помощью оптимизатора каждую эпоху обучения.

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

Для оценки точности качества обучения была использована метрика средняя абсолютная ошибка (MAPE от англ. mean absolute percentage error):

. ж-1 \т^ие - тРгеа МАРЕ = 1--^

N1

=о тах б.

И т^е\).

3.4. Обучающая выборка

Одной из важных задач при обучении свёрточной нейронной сети является создание репрезентативной обучающей выборки. Мы сгенерировали набор данных, состоящий из 120 обучающих примеров <изменение сейсмических данных, изменение скоростных моделей>, некоторые из которых изображены на (рис. 2).

В качестве базовой модели была выбрана однородная скоростная модель с резервуаром углеводородов в форме окружности. Предполагается, что после закачки вытесняющего агента, форма наблюдаемого резервуара становится эллипсом. Затем берется разность наблюдаемой и базовой модели, чтобы получить изменение моделей (рис. 3). Скоростные модели состоят из 300 х 300 точек с пространственным шагом дискретизации 10 м вдоль вертикальной и горизонтальной оси. Значения скоростных моделей варьируются от 3000 м/с до 3150 м/с. К описанным выше скоростным моделям для увеличения тренировочного набора применяются следующие трансформации: сдвиг центра окружности базовой модели, изменение радиуса окружности базовый модели, сдвиг центра эллипса наблюдаемой скоростной модели относительно центра окружности, поворот эллипса на 90°, изменение соотношения полуосей эллипсов. В конце изменение скоростных моделей нормируются к значениям отрезка [0,1] и превращаются в изображения размера 128 х128 пикселей для выхода нейронной сети. Таким образом получаем 120 изображений изменений скоростных моделей.

та

4-/

та тз

4-1

о. с

ф

тз о

Е

4-1

о.

4-/

о

=-

^^^ ;

с 0 и /■ч

о

Рис. 2. Обучающая выборка. Первая строка - данные, подающиеся на вход нейронной сети, вторая строка - модели, подающиеся на выход

нейронной сети

Рис. 3. Первая колонка - наблюдаемая скоростная модель, вторая колонка -базовая скоростная модель, третья колонка - изменение скоростной модели.

Белым цветом изображен резервуар углеводородов После генерации скоростных моделей считаются соответствующие им сейсмограммы с помощью конечно-разностных формул (4) и поглощающих граничных условий. Геометрия наблюдения состоит из 3 источников, находящихся на поверхности в левом углу, по центру и в правом углу, и 150 приемников для каждого источника, раскинутых по всей поверхности с шагом 20 м. Сейсмические данные моделируются с использованием вейвлета Рикера (2) с частотой 15 Гц для соответствующей базовой и наблюдаемой скоростной модели. Время записи 3 сек с временным шагом дискретизации 0.001 сек. Затем берется разность соответствующих сейсмических данных, чтобы получить изменение сейсмических данных (рис. 4). Мы так же нормируем значения изменений данных к отрезку [0,1] и обрезаем сейсмограммы по времени так, чтоб убрать области, где информация об изменение данных отсутствует и где содержится незначительная информация и шумы. В итоге на вход нейронной сети подаётся цветное изображение формата RGB, где каждому каналу цвета соответствует один источник. Таким образом получаем 120 изображений изменений сейсмических данных.

Рис. 4. Первая колонка - сейсмограммы для наблюдаемой модели, вторая колонка - сейсмограммы для базовой модели, третья колонка - изменений

сейсмических данных

Получившийся набор данных мы делим в соотношении 10:1:1 соответственно на тренировочную часть, на которой мы будем обучать нейронную сеть, валидационную часть, на которой мы будем смотреть, чтобы сеть не переобуча-

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

Результаты

Тренировка свёрточной нейронной сети делится на два этапа: тренировочный - здесь мы минимизируем функцию потерь, валидационный - следим за переобучением нейронной сети. Свёрточная нейронная сеть обучалась 150 эпох с начальной скоростью обучения 0.0005 и размером пакета данных 8. Процесс обучения изображен на (рис.5,6). На графике (рис. 5) изображено изменение функции потерь на тренировочной и валидационной выборке в процессе обучения. Видно, как с увеличением числа эпох функция потерь падает на обеих частях данных. К концу обучения для тренировочной части ошибка достигает зна-

—4 —2

чения 3 -10 , а на валидационной 10 . Примерно к 100 эпохе процесс обучения

не даёт уменьшения ошибки обобщения, поэтому для оценки качества обучения

на тестовом наборе зафиксирована модель с весами на 100 эпохе. На графике

(рис. 6) изображено изменение точности модели относительно метрики (18).

К концу обучения точность на валидационной выборке достигает 97%.

Epoch

Рис. 5. Функция потерь

Рис. 6. Метрика MAPE

После обучения нейронная сеть применялась к тестовому набору данных. Результаты представлены на (рис.7). Как видно, в целом, свёрточная нейронная сеть справляется с восстановлением скоростной модели по сейсмическим данным.

Заключение

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

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

го

-м ^——

лз

"О

-м

о; ——___

-м

Рис. 7. Предсказание нейронной сети. Первая, четвертая строка - входные данные; вторая, пятая - истинные модели; третья, шестая - предсказанные модели. Над предсказанными моделями указано значение метрики МАРЕ

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

Работа выполнена при поддержке Российского научного фона грант №22-21-00738.

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

1. Ronneberger O., Fisher P., Brox T. U-net: Convolutional networks for biomedical image segmentation [Электронный ресурс] - CoRR. - 2015. - режим доступа: https://arxiv.org/pdf/1505.04597.pdf, свободный. - (дата обращения: 25.05.2022).

2. Yao G., da Silva N.V., Wu D. An effective absorbing layer for the boundary condition in acoustic seismic wave simulation // Journal of Geophysics and Engineering. - 2018. - Vol. 15. - №2. - P.495-511.

3. Demanet L. Waves and imaging. Class notes - 18.367 [Электронный ресурс]. - 2021. -режим доступа: https://math.mit.edu/icg/resources/notes367.pdf, свободный. - (дата обращения: 25.05.2022).

4. Zhou D.-X. Universality of deep convolutional neural networks // Applied and Computational Harmonic Analysis. - 2019. - Vol.48. - №2. - P.787-794.

5. Ioffe S., Szegedy C. Batch normalization: Accelerating deep network training by reducing internal covariate shift // Proceedings of the 32nd International Conference on International Conference on Machine Learning - 2015.- Vol.37. - P.448-456.

6. Maas A. L., Hannun A. Y., Ng A. Y. Rectifier nonlinearities improve neural network acoustic models // Proceedings of the 30th International Conference on Machine Learning. - 2013. -Vol.28.

7. Kingma D. P., Ba L.J. Adam: A method for stochastic optimization // 3rd International Conference for Learning Representations. - San Diego, 2014. - 13 p.

© Д. С. Братчиков, К. Г. Гадыльшин, 2022