Оптимизационная методика выбора частот для получения RGB-представления результатов спектральной декомпозиции Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

Оптимизационная методика выбора частот

для получения ЯОБ-представления результатов спектральной декомпозиции

Ф.В.Краснов, А.В.Буторин

Аннотация— Авторы продолжают исследование применения методов машинного обучения к задачам геофизики. Фокусом данной работы стала процедура выбора частот для проведения КСБ-смешивания. Предыдущие работы авторов в этом направлении использовали эвристический подход к выбору частот для проведения КСБ-смешивания. Ключевой проблемой в рамках RGB-визуализации является выбор частот для смешивания. На сегодняшний день отсутствует обоснованный подход, позволяющий определить оптимальное сочетание частотных характеристик для получения наиболее информативного результата. При этом, оценка указанной информативности носит субъективный характер, что также затрудняет выбор методики для решения обозначенной проблемы.

В рамках статьи рассмотрено несколько подходов для выбора частотных характеристик с целью получения оптимального RGB-представления результатов спектральной декомпозиции. Получение оптимального выбора основано на различных предпосылках: субъективной оценке (данный подход получил название «эвристический»), распределении толщин пласта (данный подход получил название «геологический») и суммарной амплитудно-частотной характеристики используемых вейвлетов (данный подход получил название «оптимизационный»).

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

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

Ключевые слова— сейсмическое поле, оптимизационная задача, КСБ-смешивание, геофизика.

Статья получена 22 августа 2018.

Ф.В.Краснов, к.т.н., эксперт, ООО «Газпромнефть НТЦ», 190000 г. Санкт-Петербург, набережная реки Мойки д.75-79,. krasnov.fV@gazprom-neft.ru, orcid.org/0000-0002-9881-7371, РИНЦ 8650-1127

А.В.Буторин, к.г-м.н., ООО «Газпромнефть НТЦ», 190000 г. Санкт-Петербург, набережная реки Мойки д.75-79, butorin.av@gazpromneft-ntc.ru

I. Введение

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

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

Одной из наиболее распространенных методик интерпретации результатов непрерывного вейвлет-преобразования (НВП) является алгоритм цветового комбинирования RGB [8]. На вход алгоритма подается три массива, описывающих различные амплитудно-частотные характеристики. В рамках алгоритма каждому массиву присваивается свой цветокод: красного, зеленого или синего цвета. При этом отсутствие амплитуды гармоники характеризуется черным цветом, а ее максимальное значение -наибольшей насыщенностью. Далее в рамках алгоритма производится объединение цветовых каналов, таким образом, что выходной массив в каждой точке характеризуется тремя значениями амплитуды, каждой из которых соответствует свой цветовой канал. Цвет полученного на выходе дискрета определяется в рамках трехмерного цветового куба, который описывает все цвета путем комбинации красного, зеленого и синего цветового канала.

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

решения обозначенной проблемы.

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

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

II. Эвристическая методика

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

В общем случае отличия между разложением по вейвлетам будут обуславливаться разной длительностью импульса для разных частот. В случае вейвлета Риккера длительность импульса определяется следующим соотношением: dt(f) « 2/fd, где fd - доминантная частота вейвлета. В этом случае, очевидно, что наибольшие отличия между вейвлетами будут наблюдаться для низких значений доминантной частоты, с увеличением доминантной частоты отличия в длительности импульса будут снижаться, что обуславливает снижение различий между получаемыми спектральными характеристиками.

0,500

о

« 0,400 £

! 0,300 s

£ 0,200 и

0

1 0,100 О)

0,000

20 40 60 80 100 Доминантная частота, Гц

120

цветовой дифференциации объектов. Однако подобный подход является эвристическим и не обосновывается количественно, что является негативной стороной данной методики.

III. Геологическая методика

Для обоснования методики использована синтетическая модель волнового поля от контрастного выклинивающегося пласта, подробное описание модели дано в публикации [1].

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

0,015

d

> 0,01

га

Cl > 0,005

т

ли 0

S

С

_

0 50 100 150

Временная мощность клина, мс

10

20

30

40

50

Рис. 1. Зависимость длительности импульса от значения доминантной частоты.

Анализ графика на рис. 1 позволяется сделать несколько практических выводов:

• использование малого шага между доминантными значениями вейвлетов (менее 10 Гц) при разложении приводит к слабым отличиям в результатах декомпозиции;

• высокочастотные вейвлеты (/а > 50 Гц) приводят к получению близких по значениям результатов спектральной декомпозиции;

• низкочастотные спектральные характеристики (менее 10 Гц) характеризуются значительной длительностью импульса, что приводит к снижению детальности получаемого результата.

Исходя из приведенных выводов, были предложены оптимальные соотношения частот для смешивания 1525-35 Гц и 20-30-40 Гц. Данные сочетания в большинстве случаев позволяют получить приемлемое решение с точки зрения разрешающей способности и

Рис. 2. Графики распределения амплитуд гармоник в срединной части клина. Цветом показаны гармоники разных частот.

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

Таблица. 1.

Гц 10 20 30 40 50

Т, мс 0.1 0.05 0.033 0.025 0.02

Максимальное значение 0.035 0.022 0.018 0.016 0.014

Первый минимум 0.07 0.045 0.036 0.032 0.03

Побочный 0.097 0.059 0.049 0.043 0.041

Второй минимум 0.14 0.09 0.072 0.064 0.06

Побочный максимум 0.165 0.105 0.085 0.077 0.072

Начало флуктуаций 0.23 0.15 0.123 0.11 0.102

0

0,25

* 0'2 H

и

0

1 0,15 о

к (б I I

ш

0,1

£ 0,05

y = 1.6x + 0.07

0,05 0,1

Период гармоники, с

0,15

0,3 £ 0,25

H

с

но 0,2 ï о

s 0,15 к

а н

x 0 1 е

ме

1 0,05

y = 1.6 £ >x + 0.1

А/

y = 0.4x + —y-=o:24X" 0.031 + 0.014

0,05 0,1

Период гармоники, с

0,15

0

0

0

0

Рис. 3. Зависимость временной мощности, при которой наблюдается экстремум графика амплитуд, от периода анализируемой гармоники (для волнового поля с доминантной частотой 30 Гц).

Для установления закономерности возникающего «остатка» в уравнениях было выполнено повторное моделирование синтетического волнового поля с использованием сигнала Риккера с доминантной частотой 20 Гц. Далее все проделанные вычисления были повторно выполнены для нового сейсмического куба. В итоге была составлена таблица, отражающая закономерности флуктуаций в зависимости от временной мощности клина при доминантном значении сигнала 20 Гц.

Таблица. 2

Рис. 4. Зависимость временной мощности, при которой наблюдается экстремум графика амплитуд, от периода анализируемой гармоники (для волнового поля с доминантной частотой 20 Гц).

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

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

• = 1.6 * т + 0.1, для поля доминантной частоты 20 Гц

• = 1.6 * т + 0.07, для поля доминантной частоты 30 Гц

В анализируемых уравнениях значение (1.6 * т) , можно заменить величиной длительности сигнала Риккера (в соответствии с уравнением вейвлета Риккера). Нулевой остаток можно выразить через доминантное значение периода сейсмического сигнала. Таким образом, временная мощность, при которой начинают проявляться эффекты интерференции на результатах НВП могут быть примерно оценены по следующей формуле:

=2 * (Ьф + та) , где ^ - полу длительность сканирующего вейвлета Риккера (эквивалентно анализируемой частоте), а - доминантный период анализируемого волнового поля.

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

Hz 10 20 30 40 50

T 0.1 0.05 0.03 0.025 0.02

Максимальное значение 0.039 0.026 0.022 0.021 0.02

Первый минимум 0.072 0.051 0.045 0.041 0.04

Начало флуктуаций 0.254 0.18 0.15 0.137 0.128

Также целесообразным является определение временной мощности пласта, обуславливающей максимальное значение спектральной аномалии в результатах спектральной декомпозиции по алгоритму НВП. Для этого аналогичным образом были проанализированы уравнения для поля доминантной частоты 20 и 30 Гц:

• <ЛТтах = 0.25 * т + 0.014, для поля доминантной

частоты 20 Гц dTmax = 0.25 *т +0.09, частоты 30 Гц

для поля доминантной

(2 * ан1 * ) - т

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

IV. Оптимизационная методика

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

(f, в)= Y^iCi* Sf. (f),

(1)

где Sc (f) - это частотный спектр сейсмического куба, X = {inline, crosslines} Е И - индексы географических координат сейсмической трассы, sx (f) - частотный спектр сейсмической трассы с координатами X, а dim X - это размерность площади поверхности, составляющая ~ 106 ячейки 100 на 100 метров.

ячеек при размере

По аналогии с предыдущими выкладками, выразим нулевой остаток через значение доминантного периода и подставим в выражение:

dTmax = 0.25 * (Т + Td) , где т - период анализируемой гармоники (эквивалентно анализируемой частоте), а т^ - доминантный период анализируемого волнового поля.

Практическое использование полученных зависимостей подразумевает оценку характерных мощностей анализируемого пласта (dHt ) с использованием интервальной скорости (Vint), для которых определяются значения доминантной частоты вейвлета (/¿d ):

fd = 0.25 * 1

Теперь мы можем сформулировать задачу поиска наиболее точного соответствия (в) и Бс(/), как задачу поиска в для функционала (3) при котором ошибка соответствия Е минимальна:

mineE(Sc(f),Sw (f,в) ),

(3)

где Е - это функция соответствия. При этом на уравнение (3) должны быть наложены следующие ограничения (4):

frgb £ If 'minimum , fmaximum\ £ Ш

crgb £ [е, 1] £ Ш

(4)

Другими словами, частоты вейвлетов должны выбираться из определенного диапазона, который задается спектром Бс(/). Например, многие сейсмические кубы искусственно отфильтрованы в диапазоне частот выше 70 Гц. Так как, мы ищем не нулевые сгдЬ, вторым ограничением в (4) является

принадлежность crgb диапазону от положительного значения е ~ 0.01 до единицы.

некого

Поясним оптимизационную схематическом рисунке (рис. 5)

задачу

(3)

на

- _ .,.,(7.:. г--0 8

/, = 2(№,Сд = 0.2

Л = 15Яг,ц = 0.2

■■ SM.#)

где Sw (f, в) - это суммарный частотный спектр вейвлетов, s^ (f) - это частотный спектр вейвлета с частотой fi, i = {r, д, b} , а в = [crgb, frgb} Е Ш -параметры. Аргумент f показывает, что спектры являются функциями от частоты.

Волновое поле с точки зрения данных представлено сейсмическим кубом, содержащим сейсмические трассы для каждой точки координат с индексами inline, xline. Каждая сейсмическая трасса может быть преобразована с помощью дискретного преобразования Фурье в частотный спектр. Частотным спектром всего сейсмического куба будем считать усредненный частотный спектр всех содержащихся в нем сейсмических трасс (2).

Частота, Гц

Рис. 5. Иллюстрация постановки оптимизационной задачи.

На рис. 5 отображены спектры трех вейвлетов с частотами 55-20-15Гц с коэффициентами 0.8,0.2,0.2 соответственно. При их суммировании получается суммарный частотный спектр (/, в). Авторы рассмотрели следующие варианты для функции

соответствия (E).

• MSE (Mean Squared Error)

• MAE (Mean Absolute Error) На Рис. 6 приведены результаты решения оптимизационной задачи для модельного спектра (рис. 5) с использование различных функции соответствия (E).

>. о.иооь

Е 0.0004

--МАЕ

■■■■ MSE

-

100 150

частота. Гц

14

>, 12

гÎ if ю s

О 10

МАЕ

10 15 20 Итерации, шт.

25

30

Полученные оптимальные частоты составили 12-17.342.6 Гц. Частотные спектры вейвлетов и сейсмического куба представлены на рис. 8.

g 0.3

S >

^ 0.6

£ 0.2

12,0Hz 17,3Hz 42.6Hz

sлт

Sr(f)

0 50 100 150 200 250

Частота. Гц

Рис. 8. Частотные спектры вейвлетов и сейсмического куба после оптимизации.

Количественной оценкой для сравнения изображений является метрика, основанная энтропии Шэнона (5).

H = - ZkPk log2 (pfc ),

(5)

Рис. 6. Зависимость невязки от частоты для различных функций соответствия.

Поведение невязки при выборе MAE представляется более точным. Хотя точность оптимизации, достигаемая, как MAE, так и MSE не превышает 0.1%. Процедура оптимизации была выполнена с помощью алгоритма Sequential Least SQuares Programming (SLSQP), использующего метод Han-Powell quasiNewton с обновлением B-matrix [6] по алгоримту L-BFGS [7]. Выбор алгоритма SLSQP был сделан на основании исследования [5]. В качестве целевой функции была выбрана функция соответствия на основе MAE, значение которой после сходимости составило 6.59 (см. рис. 7).

16

где K - это количество оттенков серого, а р^ -вероятность оттенка с номером k. Для трех цветов (RGB) H будет суммироваться.

Н = HR + HG + Нв

Сравнение энтропии H для разных подходов к выбору частот смешивания представлено в таблице 3.

Таблица. 3

Подход Энтропия (H)

Эвристический 6.02 ± 0.3

Геологический 5.85 ± 0.1

Оптимизацио нный 5.61 ± 0.1

Рис. 7. Зависимость целевой функции от итераций.

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

V. Заключение

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

Для анализа выполнено построение RGB-карт, при этом выбор частот выполнен согласно описанным ранее методикам (рис. 9). В рамках эвристического подхода выбраны частоты 15, 25 и 35 Гц. Для геологической

оценки частот определена гистограмма значении мощности коллектора и использованы три характерных значения (20, 26 и 32 м), соответствующие доминантные частоты составили 22.4, 32.4 и 58.3 Гц. При этом необходимо отметить, что выбор характерных значений также носит субъективный характер, что является недостатком подхода. В рамках оптимизационного метода оценка частот выполнена в автоматическом режиме, согласно описанию, полученные оптимальные частоты: 12, 17.3, 42.6 Гц.

. 10 км .

Рис. 9. КОБ-карты по целевому интервалу: а) эвристический подход; б) геологический подход; в) оптимизационный подход.

Анализ полученных результатов во всех случаях позволяет закартировать перспективный геологический объект, связанный с эволюцией русла палеореки. В подобных условиях происходит формирование аккреционного комплекса, сложенного песчаными отложениями речной системы. На рис. 10 показано сопоставление ЯвВ-образа, полученного из сейсмических данных (рис. 10.1) с результатами моделирования эволюции речной системы, опубликованными в статье [12] (рис. 10.2). Таким образом, приведенное сравнение показывает геологическую интерпретацию получаемого ЯвБ-образа, полученного из волнового поля и описывающего эволюцию речной системы, которая развивалась на изучаемой территории в меловой период Мезозойской эры.

Преимуществом оптимизационного подхода является автоматизация процесса и минимизация субъективного вклада.

Библиография

1. Буторин А.В. Изучение спектральных характеристик волнового поля на примере модельных данных по результатам вейвлет-преобразования // Геофизика. 2016. №4. с 61-67

2. Витязев, В.В. Вейвлет-анализ временных рядов/ В.В. Витязев. -Санкт-Петербург: СПбГУ, 2001. -58 с.

3. Добеши, И. Десять лекций по вейвлетам. / И. Добеши. - Ижевск: РХД, 2001. — 464 с.

4. Яковлев, А.Н. Введение в вейвлет-преобразование/ А.Н. Яковлев. - Новосибирск: НГТУ, 2003. -104 с.

5. Lyu Z., Xu Z., Martins J. Benchmarking optimization algorithms for wing aerodynamic design optimization //Proceedings of the 8th International Conference on Computational Fluid Dynamics, Chengdu, Sichuan, China. - 2014. - Т. 11.

6. Powell M. J. D. A fast algorithm for nonlinearly constrained optimization calculations //Numerical analysis. - Springer, Berlin, Heidelberg, 1978. - С. 144-157.

7. Liu D. C., Nocedal J. On the limited memory BFGS method for large scale optimization //Mathematical programming. - 1989. - Т. 45. -№. 1-3. - С. 503-528.

8. Butorin A. V. et al. Spectral Inversion Methods and its Application for Wave Field Analysis (Russian) //SPE Russian Petroleum Technology Conference. - Society of Petroleum Engineers, 2017.

9. Буторин А. В., Краснов Ф. В. Сравнительный анализ методов спектральной инверсии волнового поля на примере модельных трасс //Геофизика. - 2016. - №. 4. - С. 68-76.

10. Буторин А. В., Краснов Ф. В. Возможности использования результатов спектральной инверсии при интерпретации сейсмических данных //Геофизика. - 2017. - №. 4. - С. 2-7.

11. Larkin K. G. Reflections on shannon information: In search of a natural information-entropy for images //arXiv preprint arXiv:1609.01117. - 2016.

12. Durkin P.R., Hubbard S.M., Holbrook J., Boyd R. Evolution of fluvial meander-belt deposits and implications for the completeness of the stratigraphic record // GSA Bulletin (2017) 130 (5-6): 721-739. https://doi.org/10.1130/B31699.1

Рис. 10. Сопоставление КОБ-изображения палеоречной системы (1) и моделирования эволюции современной реки, где

цветом показан относительный возраст русла (2) [12].

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

Optimization Methodology for the Selection of Frequencies to Produce an RGB Representation of the Results of Spectral Decomposition

Fedor Krasnov, Alexander Butorin

Abstract — The authors continue to study the application of machine learning methods to Geophysics problems. The focus of this work was the procedure for selecting frequencies for RGB-mixing. Previous work of the authors in this direction used a heuristic approach to the selection of frequencies for RGB-mixing. A key issue in RGB rendering is the choice of frequencies to mix. To date, there is no reasonable approach to determine the optimal combination of frequency characteristics to obtain the most informative result. At the same time, the evaluation of this information content is subjective, which also makes it difficult to choose a method for solving this problem.

The article considers several approaches to the selection of frequency characteristics to obtain optimal RGB-representation of the results of spectral decomposition. Achieving the optimal choice is based on various prerequisites: subjective assessment (this approach is called "heuristic"), the distribution of reservoir thicknesses (this approach is called "geological") and the total amplitude-frequency characteristics of the wavelets used (this approach is called "optimization").

In conclusion, the results of the practical application of the described techniques are presented, as well as the advantages of the choice of frequency characteristics by solving the optimization problem. At the same time, the authors note the absence of quantitative criteria for the evaluation of different approaches, since the use of any of the techniques allows to solve the main problem of RGB-analysis (i.e., to identify geological objects in the wave field).

Thus, the authors examined two well-known approaches to the choice of frequencies for carrying out an RGB mixing and proposed a new method based on the solution of the optimization problem, allowing you to get away from the subjective choices of the frequencies involved in the formation of a false color representation of the results of spectral decomposition.

Keywords— seismic field, optimization problem, RGB blending, Geophysics.