Сравнительный анализ методов спектральной инверсии волнового поля на примере модельных трасс

А.В. Буторин, Ф.В. Краснов, Научно-Технический Центр «Газпром нефти» (ООО «Газпромнефть НТЦ»)

Журнал «Геофизика»

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

Спектральная информация находит применение в различных областях интерпретации волнового поля для решения разнообразных задач: прогнозирования мощности коллектора [14], анализа тонких геологических объектов, таких как палеоканалы [12, 8, 9] и рифовые постройки [10], оценки затухания сейсмического сигнала [15], а также для возможного прогнозирования углеводородов по особенностям частотного состава [4,11]. Отдельно можно выделить метод СВАН-анализа, который стоит у стоков развития спектрального направления в отечественной сейсморазведке. Метод СВАН заключается в непрерывной развертке сейсмической трассы по частоте и времени, при этом разложение отвечает некоторому скользящему временному окну. Получаемые СВАН-диаграммы позволяют выполнять разделение различных типов разреза по их спектральному представлению [16]. Таким образом, широкое применение методов частотно-временного анализа делает актуальным вопрос изучения алгоритмов спектральной декомпозиции, используемых для его выполнения.

В общем случае методы частотно-временного анализа могут быть разделены на несколько классов.

Первый класс алгоритмов включает в себя технологии, основанные на использовании преобразования Гэбора, или оконного преобразования Фурье. Данные технологии используют для восстановления спектра записи преобразование Фурье, применяемое в локальном скользящем окне. [7]. Недостатком подобного анализа является значительная зависимость результатов от выбора окна анализа. Необходимо отметить, что общим недостатком данного класса методов является использование бесконечных гармонических функций для восстановления спектра сигнала, поэтому многие исследователи считают применение алгоритмов Фурье неоптимальным для изучения спектральных характеристик динамических процессов.

Второй класс алгоритмов основан на использовании вейвлет-анализа для изучения локального спектра записи. Развитие данного метода привело к появлению нового направления в области анализа сигналов — непрерывного и дискретного вейвлет-преобразования, его использование в рамках изучения сейсмических сигналов описано в работе [6]. Алгоритм вейвлет-преобразования предусматривает разложение сейсмической трассы по заданным вейвлетам определенной частоты, расчет выполняется путем подискретной корреляции участка трассы с вейвлетом заданной частоты. В рамках данной технологии используется большое количество вейвлетов, что позволяет получать различную информацию о составе сложного сигнала. Среди многообразия функций, которые могут выступать базисом разложения в алгоритме вейвлет-преобразования, можно выделить вейвлеты Риккера, Морле и Гаусса, которые традиционно используются для решения подобной задачи. Однако так же стоит отметить широкий класс дискретных вейвлетов — Daubechies, Coiflet, Symlet.

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

В рамках данной статьи рассмотрено применение методов спектральной инверсии для изучения частотного состава волнового поля, исследование выполнено с применением собственных разработок Компании «Газпромнефть НТЦ» на базе языка программирования Python.

Методы спектральной инверсии

В рамках постановки задачи спектральной инверсии используется сверточная модель, которая описывает сейсмическую трассу (s(t)) как результат свертки трассы коэффициентов отражения (r(t)) с некоторым вейвлетом (w(t)):

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

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

Представленная сверточная модель может быть выражена в матричной форме, где D — матрица вейвлетов (библиотека вейвлетов), m — матрица соответствующих вейвлет-зависимых коэффициентов отражения, n — аддитивный шум:

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

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

Под символом ||m||Lp понимается Lp норма вида , накладывающая определенное ограничение на результат решения. В этой связи обычно рассматривают несколько типов регуляризации [1]:

  • L0-регуляризация — задает количество весов отличных от нуля, в этом случае решение обладает свойством разряженности;
  • L1-регуляризация — задает суммарное значение весов в рамках решения, в этом случае решение аналогично характеризуется разряженностью;
  • L2-регуляризация — задает значение суммарной энергии весов, данный тип регуляризации был разработан А.Н. Тихоновым. В рамках регуляризации Тихонова строго нулевые веса в оптимальном решении практически невозможны.

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

  • Orthogonal matching pursuit (OMP)
  • Lasso using coordinate descent (Lasso)
  • Thresholding

Наиболее распространенным методом решения спектральной инверсии является алгоритм «согласованного преследования» (Matching pursuit), предложенный Mallat S. и Zhang Z. [13]. Суть алгоритма сводится к итеративному поиску элементов словаря, минимизирующих на каждом шаге ошибку аппроксимации. При использовании ортогонального базиса, данный метод представляет «ортогональное согласованное преследование» (Orthogonal matching pursuit) [1,2]. Данный метод решения задачи спектральной инверсии относится к L0-регуляризации и обеспечивает наиболее разреженное решение.

Среди методов L1-регуляризации можно выделить Lasso, как один из наиболее распространенных алгоритмов подобного класса. На выходе также возможно получение разреженного решения.

Также для сравнения был выбран алгоритм L2-регуляризации — Thresholding, решение, в рамках которого не характеризуется разреженностью.

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

Реализация алгоритмов на языке Python

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

В ходе исследования особое внимание было уделено структуре и содержанию словаря атомов. Определено, что наиболее точный результат аппроксимации достигается при генерации словаря по всему интервалу частот до частоты Найквиста (Nf). Определение интервала информативной части спектра проводилось на основании Фурье-образа трассы. Таким образом, размерность словаря составила [Nf, Resolution,Components], где Resolution — это число дискретов в трассе, а Components — необходимое число атомов для данной частоты. Опытным путем было установлено, что наилучшим соотношением между количеством атомов для одной частоты и числом дискретов в трассе, является 1:1. При использовании меньшего количества атомов ошибка кодирования трассы резко возрастает.

Каждый атом словаря характеризуется следующими параметрами: центр, размерность и ширина. Центры атомом словаря варьируются в пределах количества компонентов (Components). Размерность каждого атома равна количеству дискретов в трассе (Resolution). Вычисление масштаба атома (а) на основании частоты (f) производится по следующей формуле с учетом шага дискретизации сейсмической трассы (\ Δt):

Чтобы понять масштаб данных, используемых в словаре атомов, приведем следующий пример: если число дискретов в трассе 501, а информативная часть спектра варьирует от 1 до 100 Гц, то в словаре будет содержаться около 25 млн значений.

На основании словаря атомов создается инстанс кодировщика. При создании инстанса кодировщика необходимо выбрать один из алгоритмов спектральной инверсии и указать параметры, определяющие точность вычисления. Размерность получаемого в результате работы кодировщика массива коэффициентов составляет [Nf, Components] для одной трассы. Количество ненулевых коэффициентов характеризует плотность полученного массива и в зависимости от алгоритмов спектральной инверсии может быть, как задано (для L0 алгоритмов), так и определяться автоматически (для L1, L2 -алгоритмов).

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

Сравнение алгоритмов на примере синтетической трассы редких отражений

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

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

Созданная синтетическая трасса редких отражений подавалась на вход разработанного алгоритма спектральной инверсии для ее аппроксимации различными методами Lp-регуляризации. Сравнение методов с точки зрения точности аппроксимации показано на рис. 1. Как видно из иллюстрации наиболее точное восстановление обеспечивается алгоритмом OMP, алгоритм Thresholding, относящийся к методам L2-регуляризации, как отмечалось в теоретической части, не обеспечивает разреженного решения, в результате аппроксимация достигается использованием значительного количества вейвлетов.

Рис. 1. Аппроксимация исходной синтетической трассы (сверху) различными методами аппроксимации (OMP, Lasso, Thresholding)


Наличие аппроксимации библиотекой вейвлета позволяет оценить спектр трассы через распределение Винера-Вилля. Спектры, рассчитанные для каждого метода спектральной инверсии, изображены на рис. 2. Как видно из сравнения, OMP и Lasso характеризуются примерно схожим характером частотно-временного спектра, однако результат OMP лучше соотносится с исходной модельной трассой. Спектр, полученный по методу L2-регуляризации, характеризуется менее локализованными спектральными максимумами, при этом наблюдается значительное разрастание аномалий в область высоких частот.

Рис. 2. Спектры синтетичсекого сигнала, полученные по рассматриваемым алгоритмам


На сегодняшний день наиболее распространенными подходами спектрального анализа при интерпретации сейсмических данных остаются оконное преобразование Фурье и непрерывное вейвлет-преобразование. Исходя из этого, в рамках исследования было выполнено сравнение спектров, получаемых в рамках метода OMP, вейвлет-преобразования и оконного преобразования Фурье (рассмотрен вариант распределения Винера-Вилля). Результат сравнения методик, позволяет однозначно установить значительно более детальный и точный результат спектральной инверсии по сравнению со стандартными подходами декомпозиции (рис. 3).

Рис. 3. Сравнение спектров синтетической трассы, полученных по алгоритму спектральной инверсии (OMP), вейвлет-преобразованию (CWT) и оконному преобразованию Фурье (WVD)


Сравнение алгоритмов на примере интерференционной синтетической трассы

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

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

Рис. 4. Результат восстановления исходной синтетической трассы (слева) различными методами аппроксимации (справа)


Спектры суммарной трассы, полученные при помощи методов Lp-регуляризации, характеризуются схожей структурой в распределении энергии (рис. 5). Метод OMP обеспечивает наиболее локальное решение, а Thresholding, относящийся к алгоритмам L2-регуляризации, напротив характеризуется широкими аномалиями, при этом наблюдается их развитие в высокочастотную область спектра.

Рис. 5. Спектры сложной синтетической трассы, полученные с использованием различных алгоритмов (OMP, Lasso, Thresholding)


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

Оценка оптимальных параметров

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

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

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

Оценка оптимального значения L0-нормы основано на вычислении среднеквадратичного значения невязки между трассой и ее аппроксимацией. Как показывает результат тестирования, увеличение числа ненулевых коэффициентов приводит к экспоненциальному убыванию ошибки аппроксимации (рис. 6). Оптимальным значением L0-нормы является 10-20% от количества дискретов. При этом увеличение точности аппроксимации значительно сказывается на детальности результирующего спектра сигнала.

Рис. 6. Аппроксимация исходной трассы методом OMP при различном значении L0-нормы (слева), график изменения ошибки аппроксимации с увеличением значения нормы (справа)


Для алгоритма L1-регуляризации, к которым относится метод Lasso, ограничением выступает значение L1-нормы, отражающий суммарное значение ненулевых коэффициентов. В рамках исследования были протестированы значения L1-нормы от 0.05 до 0.95. Малые значения суммарного коэффициента отвечают наиболее полной аппроксимации с использованием большого числа вейвлетов, тогда как большие значения L1-нормы приводят к более разреженной аппроксимации, характеризующейся большей погрешностью. Оптимальным значением параметра регуляризации может быть принято 0.35-0.65, в данном случае обеспечивается как точность аппроксимации, так и разреженность решения.

Положения вейвлетов и их амплитуды в зависимости от значения L1-нормы изменяются, поэтому данное решение является неустойчивым, в отличие от OMP. При этом характер результирующего спектра незначительно меняется в зависимости от L1-нормы (рис. 7).

Рис. 7. Аппроксимация исходной трассы методом Lasso при различном значении L1-нормы (слева), график изменения ошибки аппроксимации с увеличением значения нормы (справа)


Для алгоритмов L2-регуляризации, к которым относится Thresholding, ограничивающим параметром выступает значение энергии коэффициентов. В рамках данного исследования был выполнен аналогичный перебор значений L2-нормы, с целью сравнения результатов при различных ограничениях решения. В данном случае минимальную ошибку обеспечивают значения нормы 0.75-0.95, при этом результат решения является нестабильным, то есть набор, и характеристики вейвлетов значительно меняются. Спектральные характеристики не значительно меняются в зависимости от значения нормы. Общей характеристикой получаемых спектров является их низкая локализация и завышение амплитуд в высокочастотной области спектра (рис 8).

Рис. 8. Аппроксимация исходной трассы методом Thresholding при различном значении L2-нормы (слева), график изменения ошибки аппроксимации с увеличением значения нормы (справа)


Влияние шумовой компоненты

Общим ограничением методов Lp-регуляризации является уровень шумовой компоненты. Исходя из этого, в рамках исследования было проведено тестирование влияния шума в данных на точность и устойчивость решения задачи спектральной инверсии. С этой целью, в трассу был добавлен случайный шум, амплитуда которого составляла 20% от максимальной амплитуды полезного сигнала. Для этого была рассчитана трасса, содержащая в себе только Гауссовский шум с амплитудой 20% от полезного сигнала. Далее полученная трасса применялась как аддитивная помеха путем суммирования трассы модельных сигналов с шумовой трассой. Полученный результат подавался на вход алгоритмов спектральной инверсии.

В результате применения описанных ранее методов спектральной инверсии, были получены частотно-временные спектры синтетической трассы с добавлением шумовой компоненты (рис. 9). Влияние шума на результаты L0- и L1-регуляризации минимально, спектры отличаются незначительно. В случае использования L2-регуляризации, влияние шума критически сказывается на результате аппроксимации, восстановление спектра сигнала в присутствии шума оказывается затруднительным.

Рис. 9. Частотно-временные спектры по различным алгоритмам без шумовой компоненты (сверху) и с добавлением 20% шума (снизу)


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

Заключение

В рамках выполненного исследования проведено сравнение методов спектральной декомпозиции на примере модельных сейсмических трасс при помощи разработанного алгоритма на языке программирования Python с использованием высокопроизводительных серверов на ОС Linux.

Методы спектральной декомпозиции могут быть разделены на три класса: первый класс включает в себя методы, основанные на использовании преобразования Фурье, ко второму классу относятся методы вейвлет-преобразования, третий класс объединяет методы спектральной инверсии, основанные на решении задачи Lp-регуляризации.

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

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

Список литературы

1) Граничин О.Н. «Рандомизация измерений и L1 оптимизация», СПбГУ, 2009

2) Луковенкова О.О. «Сравнение методов разреженной аппроксимации на примере сигналов геоакустичсекой эмиссии», Вестник КРАУНЦ. Физ.-мат. науки. 2014. № 2(9). C. 59-67. ISSN 2079-6641

3) Ягола А.Г. Обратные задачи и методы их решения. Приложения к геофизике / А.Г. Ягола, Ван Янфей, И.Э. Степанова, В.Н. Титаренко / -М.: БИНОМ. Лаборатория знаний, 2014. −216 с.

4) Castagna, J. P., S. Sun, and R. W. Siegfried, 2003, Instantaneous spectral analysis: Detection of low frequency shadows associated with hydrocarbons: The Leading Edge, 22, 120–127.

5) John P. Castagna, Shengjie Sun, Comparison of spectral decomposition methods, first break volume 24, March 2006

6) Chakraborty, A., and D. Okaya, 1995, Frequency-time decomposition of seismic data using wavelet based methods: Geophysics, 60, 1906-1916.

7) Cohen, L., 1995, Time frequency analysis: Prentice Hall PTR.

8) Giroldi, L. and F. Alegria. «Using spectral decomposition to identify and characterize glacial valleys and fluvial channels within the Carboniferous section in Bolivia.» The Leading Edge 24 (2005): 1152–1159.

9) Liu, J. and K. Marfurt. «Instantaneous spectral attributes to detect channels.» Geophysics 72 (2007): P23—P31.

10) Li, Y. and X. Zheng. «Spectral decomposition using Wigner-Ville distribution with applications to carbonate reservoir characterization.» The Leading Edge 27 (2008): 1050–1057.

11) Li, Y., X. Zheng, and Y. Zhang. «High-frequency anomalies in carbonate reservoir characterization using spectral decomposition.» Geophysics 76 (2011): V47—V57.

12) Marfurt, K. and R. Kirlin. «Narrow-band spectral analysis and thin-bed tuning.» Geophysics 66 (2001): 1274–1283.

13) Mallat, S., and Z. Zhang, 1992, Matching pursuit with time-frequency dictionaries: Technical Report 619, IEEE Transactions in Signal Processing, 41, 3397-3415.

14) Partyka, G., J. Gridley, and J. Lopez. «Interpretational applications of spectral decomposition in reservoir characterization.» The Leading Edge 18 (1999): 353–360.

15) Reine, C., M. van der Baan, and R. Clark. «The robustness of seismic attenuation measurements using fixed- and variable-window time-frequency transforms.» Geophysics 74 (2009): WA123—WA135.

16) Патент РФ № 2001123284/28, 21.08.2001. Копилевич Е.А., Давыдова Е.А., Славкин В.С., Мушин И.А., Шик Н.С. Закрытое акционерное общество «МиМГО» // Патент России № 2183335. 2003. Бюл. № 4.

Возврат к списку