Применение вейвлет-анализа для слоистого разреза Текст научной статьи по специальности «Физика»

УДК 550.834

Вестник СПбГУ. Сер. 4. 2012. Вып. 2

К. Ю. Санников, Е. Л. Лыскова, Г. В. Голикова

ПРИМЕНЕНИЕ ВЕЙВЛЕТ-АНАЛИЗА ДЛЯ СЛОИСТОГО РАЗРЕЗА

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

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

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

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

В качестве примера рассмотрим падение плоской продольной гармонической волны на верхнюю границу слоя под углом 81 к вертикали. Границы слоя считаем плоскими, идеально гладкими и параллельными друг другу. На рис. 1 видно, что, в точке Ж1

© К. Ю. Санников, Е. Л. Лыскова, Г. В. Голикова, 2012

^Ч9^ / N. Х1<Л \ // а

ХД ! кгф2! Х2 а2

^Д к

\ 1 / к

\9 '9 /

92 92 \у

Рис. 1. Падение плоской гармонической волны на тонкий слой

происходит разделение падающей волны на две когерентные волны. До точки Х2 падающая волна проходит путь si со скоростью ai, вторая волна проходит путь 2s2 со скоростью a2.

В точке Х2 падающая волна возбудит колебания

А\ cos wit---

ai

а отражённая от нижней границы слоя волна — колебание

( S2 \

Ао cos ш [t--- .

a2

Амплитуда результирующего колебания в точке Х2 определяется из выражения

A2 = Ai + A2 + 2Ai A2 cos ф,

а разность фаз колебаний

(2s2 sA ш / a2

ф = ш---= — 2s2--si

\a,2 ai J a2 \ ai

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

Учитывая, что S2 = h/ cos 82 и si = 2h tg 82 sin 6i, приходим к выражению для определения частоты максимума /max

2hm ( 1 a2 tg 62 sin 6 A 2hm , 2 x 2hm cos 62 2jt=- -----=-=- =-— (1 - sin" 02) =--,

a2 cos 62 ai a2 cos 62 a2

откуда

f — a'2 /max ~ 2ft. COS 02 '

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

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

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

сигнала на комплексные экспоненциальные функции различных частот

х(/)= J х(ь)ехр(—2л/ь)А,

—ж

где х(Ь) — сигнал; /(Ь) — частота; Ь — время. При этом ядро преобразования Фурье П = ехр(—%2п/Ь) предельно локализовано в частотной области, но абсолютно не локализовано во времени. Другими словами, в качестве базиса используется набор бесконечных во времени функций, осциллирующих со стационарными частотами и амплитудами. Как следствие, преобразование Фурье корректно передаёт информацию о периодичности сигналов при переходе из временной области в частотную. Эта особенность делает метод замечательным инструментом для исследования процессов, параметры которых — амплитуды и частоты — не меняются со временем.

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

Хя(/М)= ! х(Ь)д(Ь — Ь0)ехр(-г2п/Ь)вЬ,

где д(Ь — Ьо) — оконная функция, определяющая ширину окна.

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

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

А/ = д^-

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

Если сделать оконную функцию зависящей от частоты, то оконное преобразование Фурье перейдёт в новый класс преобразований — вейвлет-преобразований [3].

Фактически вейвлеты представляют собой короткие волновые пакеты той или иной формы, локализованные по оси времени, способные к сдвигу по ней и масштабированию (растяжению и сжатию).

Аналитически вейвлет-преобразование описывается выражением

IV(о, Ъ) = [ (Ч Л,

где x(t) — сигнал; t — время. Интегрирование ведут по всей числовой оси. Параметр a определяет ширину вейвлета и называется масштабом. Его аналогом в анализе Фурье является период гармонического колебания. Параметр b — положение вейвлета на оси времени и называется сдвигом. Функция у называется вейвлетом (базисным, материнским или анализирующим), при этом функция у, как правило, имеет компактный носитель либо её значения быстро стремятся к нулю.

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

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

Произведение Д/ • At = const справедливо для каждого из трёх преобразований. В частности, мы видим, что вейвлет-преобразование в области высоких частот обладает хорошим разрешением по времени, но плохим по частоте, в области низких частот — хорошим частотным разрешением, но плохим временным.

На рис. 3 представлен суммарный сигнал, составленный из трёх локализованных в различные моменты времени сигналов, с частотами / = 17 Гц, / = 27 Гц, /3 = 33 Гц, с добавлением некоррелированного шума. Представлены также спектр Фурье и окон-

fl-

fl

f

э

>-

t

f

f2

А б

/////

//////

/////

чип

/////

чип -*■

t

t

t

Рис. 2. Частотно-временная локализация преобразований с разными анализирующими функциями: о — преобразование Фурье; б — оконное преобразование Фурье; в — вейвлет-преобразование; г — анализирующие функции преобразования Фурье (вверху) и вейвлет-преобразования (внизу)

a

t

t

2

в

t

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

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

Ш(а, Ь) = \Ш(а, Ь)\ ехр (гф(а, Ь)) = Wm(a, Ь) ехр (гф(а, Ь)).

Для вейвлет-преобразования существует аналог равенства Парсеваля, в котором величина, стоящая справа, представляет собой энергию сигнала:

/ ./' / :'// / /|Ща,6)|2^,

— Ж —Ж —Ж

40

30

20

10

Рис. 3. Различные представления сигнала:

а — амплитудно-временное представление; б — три составляющие суммарного сигнала; в — скалограмма вейвлет-преобразования сигнала; г — двумерное представление спектра оконного преобразования Фурье; д — амплитудный спектр классического преобразования Фурье сигнала

соответственно величина Е(а, Ь) = \Ш(а, Ь)\2 представляет собой локальную плотность энергии.

Функцию Е(а^Ь]), описывающую распределение энергии по масштабам и сдвигам, называют скалограммой. Данная функция позволяет проследить изменение энергии сигнала с течением времени.

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

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

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

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

Особенностью данных трасс являлось наличие хорошо отражающих горизонтов на времени около 0,1, 0,2 и 1 с и слабо дифференцированная толща в диапазоне 1,5-3,3 с (рис. 4).

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

секунды

Рис. 4. Примеры сейсмических трасс исследуемого района, полученных методом общей глубинной точки на трёх пикетах

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

1 ( Л - Ь\ 1 (г - Ь

¥(*) = —¡= ехр г2я/с- ехр -- -

л/ая \ а I \ 2 \ а

где а — масштаб анализирующего вейвлета; Ь — временной сдвиг; /с — центральная частота вейвлета, связанная с масштабом а, шагом дискретизации с! и псевдочастотой Фурье / выражением / = /с/(а!).

Коэффициенты вейвлет-преобразования вычислялись по формуле

ц) = ¿ -¿ЕЕ-^-р ехр [Л

г=1 у ] = ! к=1 4 V

где N — длина временного ряда в отсчётах; г — наибольший масштаб анализирующего вейвлета.

На рис. 5 показано трёхмерное представление коэффициентов вейвлет-преобразования в собственных координатах а, Ь (масштаб, отсчёты). Видно, что максимумы интерференционных волн характеризуются четкими, хорошо различимыми пиками.

Рис. 5. Трёхмерное представление значений коэффициентов вейвлет-преобразования

10—| 8-

й

^ 6-

з -

§ 44 < -

240

10 0 0

200 300

ОГСЧ&ГЙ

Локальные максимумы ^млх(ал, Ьj) определялись по восьми ближайшим точкам

УУт(аиЪз-1) ^т ( Я; • bj) УУт(аиЪэ+1)

№т{аг-1, Ъз) Wm (а—1, Ь^+1)

с помощью формулы

^млх(аь Ь) = <

Шт(а^,Ь]), если'

№т(а±1 ,Ь±) < ~№т(аиЦ) > ~№т(а1Т1,Ьз±1) Wm(aiTl,Ьj) < Wm(ai,Ьi) > Wm(auЬj±l) в противном случае

2

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

Каждый локальный максимум характеризуется тремя параметрами или координатами на частотно-временной плоскости: временем вступления максимума, которое даёт координату верхней границы слоя (см. рис. 1); значением или амплитудой локального максимума, которые характеризует интенсивность отражённой от слоя волны; частотой локального максимума, являющейся функцией мощности слоя h и скорости волны в слое.

В результате, мы можем определить мощность слоя с точностью до некоторого множителя П:

a2 П

Г1 = - = -

2fmax COS 82 2fmax

Так как h = a2At, где At — время, которое требуется волне на прохождение расстояния h, то мощность слоя во временных координатах, с точностью до некоторого множителя Л, определится из выражения

1 1 Л

Ath =

2fmax COS 82 2fmax

Выполнив расчёт всех локальных максимумов по пикетам, получаем набор значений qkp(tkp, Athkp, Wkp) — координаты слоёв временного разреза, где tkp — собственно координаты слоёв во временном масштабе; Athkp — мощности слоёв; Wkp — интенсивности волн, отражённых от слоёв разреза. Как следует из приведённых формул, с глубиной мощности слоёв будут закономерно возрастать за счёт множителя, обратно пропорционального углу падения волны на слой. Вследствие этого мощности слоёв назовём кажущимися мощностями.

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

Как было сказано выше, по координатам локальных максимумов (tkp — время вступления, fkp — частоты) были рассчитаны значения кажущихся мощностей слоёв Athkp на всех пикетах и построен сейсмический профиль (рис. 7).

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

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

а

[-ч

Б

£

¡У

20 10

0 -10 -20

50 45

40

35

30

25

20

15

10

5

0,5 1 1,5 2 2,5 3 3,5 с

Рис. 6. Сейсмическая трасса (ОГТ) (а), скалограмма вейвлет-преобразования (б): крестиками отмечены локальные максимумы интерференционных волн

Рис. 7. Сейсмический профиль, построенный по методу локальных максимумов амплитуд интерференционных волн

0,5 1

3 3,5

100 150

Пикеты

Очистка профиля от подобного «шума» осуществлялась, исходя из следующих соображений. Сейсмический профиль, по сути, представляет собой визуализацию разреженной матрицы данных (к — количество строк или общее число временньгх отсчётов, p — число столбцов или общее количество пикетов), т. е. матрицы, состоящей из большого числа нулевых элементов, ненулевые элементы которой представляют локальные интенсивности или локальные энергии отражённых от слоёв волн. Из данной матрицы последовательно выбирались матрицы Mnm, при этом, если количество нулевых элементов было больше некоторого порогового значения, все элементы данной матрицы заменялись нулями. Очищенный подобным образом сейсмический профиль представлен на рис. 9.

0,5 1

■ 1,5

! 2 р

>2,5 3 3,5

100 Пикеты

Рис. 9. Очищенный от «шума» сейсмический профиль, построенный по методу локальных максимумов интерференционных

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

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

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

Литература

1. Силкин К. Ю., Дубянский А. И. Оценка возможности применения вейвлет-анализа к сейсмологическим данным // Вестн. ВГУ. Сер.: Геология. 2007. № 1. С. 138-140.

2. Филатова А. Е., Артемьев А. Е., Короновский А. А. и др. Успехи и перспективы применения вейвлетных преобразований для анализа нестационарных нелинейных данных в современной геофизике // Изд. вузов «ПНД». 2010. Т. 18, № 3. С. 3-23.

3. Витязев В. Вейвлет-анализ временных рядов. СПб.: Изд-во СПбГУ, 2001. 61 с.

4. АстафьеваН. Вейвлет-анализ: основы теории и примеры применения // Усп. физ. наук. 1996. Т. 166, № 11. С. 1145-1170.

5. Воробьёв В., Грибунин В. Теория и практика вейвлет-преобразования. СПб.: Изд-во ВУС, 1999. 204 с.

6. Добеши И. Десять лекций по вейвлетам / пер. с англ. Ижевск, 2001. 464 с.

7. ДреминИ., Иванов О., Нечитайло В. Вейвлеты и их использование // Усп. физ. наук. 2001. Т. 171, № 5. С. 465-561.

8. Чуи К. Введение в вейвлеты / пер. с англ. М.: Мир, 2001. 412 с.

Статья поступила в редакцию 27 декабря 2011 г.