Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/~rastor/Study/Vityazev_Wavelets.pdf Дата изменения: Mon Feb 2 16:32:51 2009 Дата индексирования: Mon Oct 1 23:03:08 2012 Кодировка: Windows-1251 Поисковые слова: п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п п р п р п р п р п р п р п р п р п р п р п р п р п р п р п р п

Санкт-Петербургский государственный университет

В.В. Витязев

ВЕЙВЛЕТ-АНАЛИЗ ВРЕМЕННЫХ РЯДОВ
Учебное пособие

Издательство С.-Петербургского университета
2001


ББК 22.6 В54 Рецензенты: д-р физ.-мат. наук, проф. В.А. Гаген-Торн (C.-Петерб. гос. ун-т), д-р физ.-мат. наук, проф. К.В. Холшевников (C.-Петерб. гос. ун-т)

Печатается по постановлению Редакционно-издательского совета С.-Петербургского государственного университета

В54

Витязев В.В.

Вейвлет-анализ временных рядов: Учеб. пособие. СПб.: Изд-во С.-Петерб. ун-та, 2001. 58 с.
Рассмотрены практические аспекты использования вейвлет-преобразования в анализе временных рядов. Показаны аналогии и различия между Фурье-преобразованием и вейвлет-преобразованием. Изучены свойства скалограмм (вейвлет-спектров) для дискретного белого шума и полигармонических функций, получены критерии выделения сигнала из шума. Приводятся скалограммы типичных временных рядов: синусоидального сигнала, дискретного белого шума и авторегрессионного процесса. Разработан алгоритм анализа временных рядов, основанный на применении вейвлета Морле. С помощью этого алгоритма проведен анализ среднегодовых чисел Вольфа и кривой блеска квазара 3С273. Даны примеры и упражнения. Пособие предназначено для студентов, аспирантов и научных сотрудников, занимающихся обработкой экспериментальных данных.

ББК 22.6
c В.В.Витязев, 2001 c Издательство С.-Петербургского университета, 2001

.


.

2


Wavelets had... no denition in the glossary. Thus I provide my own crude one: wavelet transforms are a lot like Fourier transforms, but you get more choices. Virginia Trimble (1997).

Введение
Настоящее учебное пособие посвящено изучению временных рядов с помощью вейвлет-преобразования, т. е. разложения одномерного сигнала по базису, сконструированному из солитоноподобных функций (вейвлетов), посредством их масштабных изменений и переносов вдоль оси времени. Считается, что принцип вейвлет-анализа был впервые изложен в работе Гроссмана и Морле (1984), и с тех пор вейвлет-анализ стал одним из самых популярных разделов математики и ее приложений. Помимо временных рядов вейвлеты нашли широкое применение в задачах фильтрации и чистки многомерных сигналов, в анализе изображений, в сжатии больших массивов информации и т. д. Теория и практика вейвлет-анализа описана в многочисленных книгах и статьях (Добечи 1992; Койфман, 1992; Фостер, 1996; Скаргл, 1997). Созданы мощные программные пакеты (например, MIDAS, MATHLAB), позволяющие производить различные операции, свойственные вейвлет-анализу. Основная идея вейвлет-преобразования отвечает специфике многих временных рядов, демонстрирующих эволюцию во времени своих основных характеристик среднего значения, дисперсии, периодов, амплитуд и фаз гармонических компонентов. Подавляющее число процессов, изучаемых в астрономии, обладают такими свойствами: блеск квазаров, солнечная активность, неравномерность вращения Земли вот далеко не полный перечень примеров. Из литературы о вейвлетах, опубликованной русском языке, следует отметить обстоятельный обзор Н.М.Астафьевой (1996). Тем не менее, сейчас ощущается острый дефицит литературы о вейвлетах учебного характера. Настоящее пособие, посвященное практической стороне использования вейвлет-преобразований и ориентированное на исследование астрономических временных рядов, призвано в какой-то мере устранить этот недостаток.

3


От Фурье-анализа к вейвлетам
Классический анализ Фурье основан на возможности исследования функций во временной (|t| < ) и частотной (| | < ) областях с помощью прямого и обратного преобразований Фурье


^ f ( ) =
-

f (t)e

-i2 t

dt,

(1)

f (t) =
-

^ f ( )e

i 2 t

d.

(2)

Например, для функции

f (t) = A cos(2 0 t)
преобразование Фурье имеет следующий вид:

(3)

^ f ( ) = A [ ( - 0 ) + ( + 0 )],

(4)

где ( ) дельта-функция. Этот пример демонстрирует замечательную способность преобразования Фурье фокусировать в точку "размазанную"по времени информацию о периодичности функции при переходе из временной области в частотную. Достигается это за счет того, что ядро преобразования Фурье, т. е. функция e-i 2 t , не локализовано во времени, но имеет предельную локализацию в частотной области. Это обстоятельство и делает преобразование Фурье прекрасным инструментом для изучения процессов, свойства которых не меняются со временем (в нашем примере это обеспечивается условиями A = const и 0 = const). Однако именно это обстоятельство делает преобразование Фурье плохим методом для исследования иррегулярных функций, т. е. функций, характеристики которых эволюционируют во времени. Например, преобразование Фурье не отличает сигнал, представляющий собой сумму двух синусоид, от сигнала, состоящего из тех же синусоид, но включающихся последовательно (рис.1). Для устранения этого недостатка нужно локализовать преобразование Фурье на промежутках конечной длины. Таким приемом 4


Рис. 1: Неоднозначность преобразования Фурье: a модельный ряд

сумма двух синусоид c частотами 1 = 0.062 и 2 = 0.105 Гц; b периодограмма суммы этих синусоид; с те же синусоиды, включающиеся последовательно; d периодограмма синусоид, включающихся последовательно; e квадрат модуля преобразования Гэбора как функция частоты и времени при = 10.

5

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


GT (, b, ) =
-

f (t)e

-

(t-b)2 2

e

-i2 t

dt.

(5)

Сравнивая это выражение с (1), мы видим, что введенное под знак интеграла временное окно e- 2 выделяет лишь небольшой отрезок исходного временного ряда с центром в произвольной точке b и, тем самым, позволяет получить эволюцию спектра во времени (рис.1). Здесь важно подчеркнуть, что окно преобразования Гэбора имеет постоянную ширину, которая определяется параметром . Эффективная ширина окна определяет длину интервала T , который дает главный вклад в значение интеграла в выражении (5). Длина T является мерой временного разрешения, в то время как ширина спектральной линии определяет меру частотного разрешения. Известно, что обе эти характеристики связаны между собой соотношением
(t-b)2

1 . (6) T Отсюда видно, что естественное для анализа нерегулярных сигналов стремление повысить временное разрешение всегда приводит к уменьшению разрешающей способности в области частот. К этому следует добавить, что при использовании преобразования Гэбора возникает проблема выбора ширины окна во временной области. Слишком широкое окно может обеспечить разумное представление низкочастотных компонентов ряда, но его ширина будет избыточной для гармоник с высокой частотой, поскольку все интересные нерегулярности в высокочастотной области спектра сгладятся. Наоборот, достаточно узкое окно даст возможность изучить вариации во времени высокочастотных компонентов, но оно не будет адекватным для низкочастотных гармоник. Если сделать оконную функцию зависящей от частоты так, чтобы для низких частот окно становилось шире, а для высоких уже,

6


то оконное преобразование Фурье переходит в новый класс преобразований, который и получил название вейвлет-преобразования. Термин "вейвлет"происходит от английского слова wavelet, буквальный перевод которого означает маленькая волна (сравните booklet, starlet). Хорошим аналогом в русском языке является слово всплеск. Переход от слова "волна"(wave) к слову "всплеск"(wavelet) отражает суть перехода от преобразования Фурье к вейвлет-преобразованию (рис. 2). И. Добечи (1992) начинает свою книгу "Ten lectures on wavelets"такими словами: "The wavelet transform is a tool that cuts up data or functions or operators into dierent frequency components, and then studies each component with a resolution matched to its scale". Основное здесь то, что вейвлет-преобразование не просто "режет"исследуемый объект на куски, а выделяет из него компоненты разных масштабов и что каждый компонент анализируется с той степенью детальности, которая соответствует его масштабу. С этим свойством вейвлет-анализа мы встречаемся не только при исследовании временных рядов. Например, рассматривая Исаакиевский собор издалека, мы не увидим украшений его стен, зато сможем оценить общую композицию, т. е. крупномасштабную структуру сооружения. Наоборот, приблизившись к нему вплотную, мы сможем разглядеть в деталях все украшения, но потеряем при этом ощущение формы самого собора.

Основы вейвлет-анализа
Дадим теперь формальные определения, лежащие в основе вейвлетанализа. Определение 1. Интегральным вейвлет-преобразованием функции f (t) L2 (R) называется выражение

1 W (a, b) = 1 |a|
где a, b R, a = 0.

/2 -

f (t)



t-b a

dt,

(7)

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

7


Рис. 2: Волны и всплески. Преобразование Фурье это корреляция меж-

ду исходным рядом a и волной b. Волна характеризуется значением частоты, поэтому преобразование Фурье является функцией одной переменной частоты. Вейвлет-преобразование это корреляция между исходным рядом a и всплеском c. Вейвлет (всплеск) характеризуется масштабом и локализацией на оси времени, поэтому вейвлетпреобразование зависит от двух аргументов масштаба вейвлета и его положения на оси времени.

8


вейвлетом). Заметим, что в формуле (7) символом обозначена процедура комплексного сопряжения. Параметр a определяет размер вейвлета и называется масштабом (scale). Его аналогом в Фурье-анализе является период (частота) гармонического колебания. Следует сказать, что понятие масштаба более широкое (хотя и менее наглядное), чем понятие периода. Связано это с тем, что в Фурье-преобразовании функциональный вид ядра преобразования зафиксирован раз и навсегда, в то время как вейвлет-преобразование одной и той же функции можно получить с помощью различных базисных вейвлетов (т. е. в разных системах масштабов). Параметр b задает временную локализацию вейвлета и называется сдвигом (shift). Этот параметр не имеет аналога в Фурьепреобразовании. Определение 5. Обратное интегральное вейвлет-преобразование задается выражением


f (t) = C

-1 - -

W (a, b)

t-b a

1 dadb , a2 a
1 /2

(8)

где C нормирующий коэффициент:


C =
-

^ | |2 | |-1 d < .

(9)

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

9


Свойства вейвлетов
Частотно-временная локализация
Как было сказано выше, частотно-временная локализация является характерной особенностью анализирующих вейвлетов. Это означа^ ет, что вейвлеты (t) и их преобразования Фурье ( ) существенно отличаются от нуля лишь на малых интервалах времени и частоты и очень мало отличаются от нуля (или просто равны нулю) вне этих интервалов. Количественной мерой локализации функции z (t) L2 (R) могут служить ее центр t и радиус t :

1 t= ||z || 1 ||z ||2

2 -

t |z (t)|2 dt,

(10)



2 = t

[t - t ]2 |z (t)|2 dt.
-

(11)

При этом эффективная ширина вейвлета принимается равной 2t .

Нулевые моменты
Чтобы обеспечить обратимость вейвлет-преобразования, вейвлет должен удовлетворять соотношению (9). Это обеспечивается условием

^ (0) = 0,
откуда следует


(12)

(t)dt = 0.
-

(13)

Для приложений бывает важно, чтобы не только нулевой момент, но и m старших моментов были равны нулю:


(t)tm dt = 0.
-

(14)

10

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

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

MHAT-вейвлет
Этот вейвлет получается двукратным дифференцированием функции Гаусса: 2 2 d2 (t) = 2 e-t /2 = (1 - t2 )e-t /2 . (15) dt Преобразование Фурье этого вейвлета имеет вид 2 ^ ( ) = 2 2 e- /2 . (16) Графики этих функций показаны на рис. 3. Название этого вейвлета MHAT (Mexican HAT) произошло от характерного вида его графика, напоминающего сомбреро. MHAT-вейвлеты хорошо локализованы как во временной, так и в частотной области. Центры и радиусы локализации в обеих

Рис. 3: Вейвлет MHAT (слева) и его преобразование Фурье (справа). 11


областях имеют следующие значения:

t = 0; = 1.51;

t = 1.08; = 0.49.

(17) (18)

Отметим также, что MHAT-вейвлет имеет нулевые значения нулевого и первого моментов.

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

(t) = e

-t2 /

2

[e

ik0 t

-e

2 -k0 2 /4

], ].

(19) (20)

^ ( ) = [e

-2 (k0 - )2 /4

- e-

2 (k0 + )2 /4

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

(t) = e

-t2 /

2

e

i 2 t

, .

(21) (22)

^ ( ) = e

-2 (2 - )2 /4

Рис. 4: Вещественная часть вейвлета Морле (слева) и его преобразование Фурье (справа).

12


Графики этих функций показаны на рис. 4. Центр и радиус локализации вейвлета Морле во временной области определяются значениями

t = 0;

t = /2.

(23)

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

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


(t) =
-

h(t ) f (t - t )dt ,

(24)

где h(t) весовая функция фильтра, удовлетворяющая условию


h(t)dt = 1.
-

(25)

Преобразуя (24) по Фурье, получаем

^ ^ ^ ( ) = h( ) f ( ).

(26)

^ В этом равенстве h( ) передаточная функция фильтра (24), которая является Фурье-образом весовой функции:


^ h( ) =
-

h(t)e-

i t

dt.

(27)

Например, для ступенчатой весовой функции 1 2T , |t| T , h(t) = 0, |t| > T , 13

(28)


передаточная функция имеет вид

sin(T ) ^ h( ) = . (29) T ^ В общем случае из (25) следует h(0) = 1. Таким образом, оператор (24) задает низкочастотную фильтрацию. Отметим, что ширина пропускания передаточной функции (29) определяется только длиной весовой функции фильтра, т. е. значением параметра T . Сравнивая преобразования (7) и (24), мы видим, что интегральное вейвлет-преобразование (7) можно также считать фильтром исходной функции f (t). Применяя к (7) преобразование Фурье по переменной b, получаем ^ ^ ^ W (a, ) = a (a ) f ( ), (30)
откуда видно, что передаточная функция вейвлет-преобразования имеет вид ^ H ( ) = a (a ). (31) В силу (12) из последнего равенства имеем

H (0) = 0.

(32)

К этому надо добавить, что для каждого значения масштаба a > 0 максимум функции H ( ) лежит на частоте

max =

C , a

(33)

где величина C зависит от взятого вейвлета (например, для MHAT вейвлета C = 2, для вейвлета Морле C = k0 = 2 ). Кроме того, эффективная ширина функции H ( ) убывает с ростом масштаба a. Другими словами, вейвлет-преобразование это полосовой фильтр с переменной шириной полосы пропускания (узкой для больших масштабов и широкой для малых масштабов). Сравнение передаточных функций (29) и (31) показано на рис. 5. Таким образом, для каждого значения масштабного коэффициента a > 0 величины W (a, b), определенные с помощью вейвлетпреобразования (7), представляют собой результат фильтрации исходной функции в диапазоне частот, центр которого определяется значением масштабного коэффициента, а размер свойствами принятого анализирующего вейвлета. 14


ция низкочастотного фильтра (29) при T = 1 (штриховая линия) и при T = 2 (сплошная линия). Справа: передаточная функция (31) полосового фильтра вейвлет-преобразования MHAT при a = 0.25 (сплошная линия) и при a = 1 (штриховая линия).

Рис. 5: Передаточные функции фильтров. Слева: передаточная функ-

Спектральные характеристики
^ Будем считать, что две функции f (t) и f ( ) связаны между собой соотношениями Фурье (1)-(2). В этом случае справедливо равенство Парсеваля


|f (t)|2 dt =
-

1 2



^ |f ( )|2 d .
-

(34)

Интеграл, стоящий в левой части этого равенства, называется энергией сигнала f (t). На этом основании величину

E ( ) =

1^ |f ( )| 2

2

(35)

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


|f (t)| dt = C
-

2

-1 - -

|W (a, b)|

2

da db . a2

(36)

Поэтому величину

E (a, b) = |W (a, b)|
15

2

(37)


уместно также назвать плотностью спектра энергии, однако, в отличие от (35), величина E (a, b) определяет спектральную характеристику не только для заданного масштаба, но и для параметра сдвига b. По этой причине ее называют локальным спектром энергии. В противоположность этому величину


Ew (a) =
-

|W (a, b)|2 db

(38)

называют глобальным спектром энергии. Эта характеристика показывает распределение энергии по масштабам и является аналогом плотности спектра энергии E ( ). Можно показать, что обе эти характеристики связаны между собой следующим соотношением:


Ew (a) = a
-

^ E ( )| (a )|2 d .

(39)

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

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

fk = f (tk ),

tk = t k ,

k = 0, 1, ..., N - 1.

(40)

16


Для оценки вейвлет-преобразования этой последовательности воспользуемся следующим выражением:

WA (a, b) =
где

1 n(a, b)

N -1

fk
k=0



tk - b a
2

,

(41)

N -1

n(a, b) =
k=0

e

-

1 B

tk -b a

,

(42)

причем B = 2 для MHAT-вейвлета и B = 2 для вейвлета Морле. При переходе от (7) к (41) мы убрали из знаменателя формулы (7) множитель a, заменив его выражением


e
-

-

(t-b)2 a2 B

dt = a B .

(43)

Этим самым была устранена зависимость амплитуд гармонических компонентов от параметра a, что обычно мешает правильно оценить их относительные интенсивности по графическому представлению вейвлет-спектров. Кроме того, функция n(a, b) как аппроксимация (43) позволяет "уравнять"для различных значений масштабного коэффициента a число отсчетов исходной функции, участвующее в вычислении. Следуя Фостеру (1996), будем называть оценку (41) амплитудной вейвлет-функцией. Эта функция вычисляется на дискретном множестве значений аргументов ai и bj , i = 0, ...Na - 1; j = 0, ...Nb - 1 (конкретные способы дискретизации параметров a и b обсуждаются подробно в п. 8 основного алгоритма (см. с. 27). Используя (41), введем оценку локального спектра энергии

S (ai , bj ) = |WA (ai , bj )|2 .

(44)

Эту функцию обычно называют скалограммой (scalogram), подчеркивая тем самым ее способность описывать распределение энергии по масштабам. Поскольку это распределение локализовано во времени с помощью параметра сдвига b, уместно называть (44) локальной скалограммой, однако такой термин не нашел широкого распространения. 17


Очевидно, что на основе скалограммы S (ai , bj ) можно ввести также и оценку глобального спектра энергии

G(ai ) =

1 N

S (ai , bj ),
j

(45)

где N число точек, по которому осуществляется осреднение (подробности изложены в п. 13 основного алгоритма (см. с. 31). По предложению Скаргла (1993, 1997) функцию (45) называют скейлограммой (scalеgram). В силу соотношения (39) скейлограмма в вейвлет-анализе является прямым аналогом сглаженной периодограммы в Фурье-анализе.

Вейвлет-преобразование синусоиды
Фурье-анализ и вейвлет-анализ временных рядов основаны на переводе имеющейся информации из временной области в частотную область или область, задаваемую параметрами a, b вейвлет-преобразования. При анализе полигармонических функций необходимо знать, какой вид будет иметь простая синусоида в этих областях. В Фурье-анализе показывается, что образом синусоиды в периодограмме Шустера является спектральная линия, размер и форма которой определяются длиной ряда и расположением точек внутри него. Изучим аналогичную ситуацию в вейвлет-анализе. Для этого рассмотрим синусоиду

f (t) = A sin(0 t) =

A [e 2i

i0 t

-e

-i0 t

].

(46)

Для Фурье-образа этой функции имеем ^ f ( ) = A [ ( - 0 ) - ( + 0 )] . i Применим к (30) обратное преобразование Фурье:

(47)

a W (a, b) = A 2





^ ^ (a ) f ( )e
-

i b

d .

(48)

Подстановка (47) в (48) дает a i0 b ^ W (a, b) = A e (a0 ) - e 2i 18

-i0 b

^ (-a0 ) .

(49)


^ В тех случаях, когда функция ( ) действительная и четная, выражение (49) упрощается: ^ W (a, b) = A a (a0 ) sin(0 b).
(50)

Формулы (49) и (50) соответствуют определению непрерывного вейвлетпреобразования (7). Численные оценки выполняются по формуле (41). Очевидно, что с учетом описанных на с. 17 обстоятельств вывода формулы (41) выражения (49) и (50) следует переписать следующим образом:

WA (a, b) =

A e 2i B

i0 b

^ (a0 ) - e

-i0 b

^ (-a0 ) ,

(51) (52)

A^ WA (a, b) = (a0 ) sin(0 b), B

где B = 2 для MHAT-вейвлета и B = 2 для вейвлета Морле. Рассмотрим теперь эти формулы для конкретных вейвлетов.

MHAT-вейвлет
Подставляя (16) в (52), получаем
2 WA (a, b) = A a2 0 e
2 - a 2 0 /2

sin 0 b,

(53) (54)

4 S (a, b) = A2 a4 0 e

- a2

2 0

sin2 0 b.

Мы видим, что при фиксированном значении частоты 0 нашей синусоиды ее MHAT-преобразование есть синусоида той же частоты с амплитудой, зависящей от текущего значения масштабного коэффициента. Легко найти максимумы функций (53) и (54) по шкале масштабов. Действительно, приравнивая к нулю результаты дифференцирования этих функций по параметру a, получаем

amax 0 =



2,

(55)

откуда для периода P0 = 2 /0 находим P0 = 2 amax 4.44a 19

max

.

(56)


Вейвлет Морле
Подставляя (22) в (52), получаем

WA (a, b) =

A e 2i

i0 b -

e

2 4

(2 -a0 )2

-e

-i0 b -

e

2 4

(2 +a0 )2

,

(57)

2 2 A2 - 2 (2-a0 )2 e2 + e- 2 (2+a0 ) . (58) 4 В отличие от предыдущего случая мы видим, что локальный спектр энергии синусоиды не зависит от сдвига b и определяется симметричным контуром, максимум которого связан с периодом синусоиды следующим образом:

S (a, b) =

P0 = amax ,
при этом интенсивность линии в точке максимума равна

(59)

S

max

(amax , b) =

A2 . 4

(60)

Здесь полезно подчеркнуть, что хотя в общем случае связь периода синусоиды и значения amax , а также связь амплитуды и максимума линии зависят от вида анализирующего вейвлета, для вейвлета Морле эти связи точно такие же, как и для обычного Фурьепреобразования. Это обстоятельство говорит о том, что вейвлет Морле предпочтительнее использовать в анализе временных рядов, чем MHAT-вейвлет. Примеры скалограмм синусоид приводятся на с. 37 (рис.11).

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

20


Будем считать, что временной ряд

xk = x(tk ),

tk = t k ,

k = 0, 1, ..., N - 1,

(61)

представляет собой выборку некоррелированных значений случайной величины, распределенной по нормальному закону, с нулевым 2 математическим ожиданием и дисперсией 0 (дискретный белый шум): 2 0 , p = q , xp xq = (62) 0, p = q . Получим закон распределения отсчетов скалограммы по масштабам для вейвлета Морле. Введя обозначение

tk - b , a представим скалограмму следующим образом: k = S (a, b) = P 2 (a, b) + Q2 (a, b),
где
N -1

(63)

(64)

P (a, b) =

1 n(a, b) 1 n(a, b)

xk e
k=0 N -1

-k /

2

cos 2 k ,

(65)

Q(a, b) =

xk e
k=0

-k /

2

sin 2 k .

(66)

Нетрудно показать, что при фиксированном значении параметра b функции P (a, b) и Q(a, b) представляют собой нормально распределенные случайные величины с нулевым математическим ожиданием и дисперсией 2 2 = 0 Z (a, b), (67) 2 где N -1 2 1 Z (a, b) = 2 e-2k / . (68) n (a, b)
k=0

Введем в рассмотрение нормированную скалограмму

s(a, b) =

S (a, b) . 2 0 Z (a, b)
21

(69)


Величина 2s(a, b) имеет хи-квадрат распределение с двумя степенями свободы, откуда следует, что

нормированная скалограмма s(a, b) представляет собой случайную величину по переменной a, закон распределения которой имеет вид p(x) = e-x , 0 < x < . (70)
Знание закона распределения нормированной скалограммы позволяет получить критерий выделения сигнала из шума. Действительно, вероятность того, что отсчеты нормированной скалограммы превзойдут заданную величину T , определяется формулой


P r{s(a, b) > T } =
T

p(x) dx = e

-T

.

(71)

Задавая теперь уровень значимости q << 1, определяющий вероятность редкого события превышения отсчета нормированной скалограммы заданной величины, для порога обнаружения сигнала имеем Tq = - ln q . (72) Другими словами, с вероятностью p = 1 - q можно констатировать, что те значения скалограммы S (a, b), которые удовлетворяют неравенству 2 (73) S (a, b) > 0 Z (a, b) Tq , генерируется не шумом, а сигналом. Отметим, что функция Z (a, b) практически не зависит от значений сдвига b. Характерный вид функции, стоящей в правой части формулы (73), показан на рис. 8 (см. с. 32) при N=100 и t = 1.

22


Алгоритм вейвлет-анализа временных рядов
Исходные данные: 1) равномерный временной ряд

xk = x(tk ),

tk = t k , k = 0, 1, ..., N - 1,

где t шаг выборки, N число точек ряда; 2) критический порог для разделения шумового и детерминированного компонентов временного ряда Tq , отвечающий уровню значимости q (см. соотношение ( 72)).

Ниже мы опишем основные шаги алгоритма и проиллюстрируем их исполнение для модельного временного ряда
xk = A1 cos(2 1 tk - 1 ) + n nk , tk = t k , k = 0, 1, 2, . . . , N - 1,
(74)

где

A1 , 1 , 1 , амплитуда, частота и фаза гармонического компонента; nk , k = 0, 1, 2, . . . , N - 1, значения случайной величины с нулевым математическим ожиданием и единичной дисперсией, распределенной по нормальному закону (шумовая составляющая ряда); n среднеквадратическое отклонение шумового компонента nk , задаваемое с помощью соотношения n = A2 1 , 2
(75)

где отношение "сигнал к шуму".

Заданные параметры ряда: t = 1c; N = 100; q = 0.01; A1 = 1; 1 = 0.1 Гц; 1 = 0;
23


= 0.50;
Алгоритм спектрального анализа состоит из следующих шагов: 1. Графическое представление исходного ряда во вре-

менной области

Обычно визуальное изучение графика исходного ряда позволяет обнаружить либо присутствие в данных постоянного слагаемого, либо низкочастотный компонент (тренд). Хотя вейвлеты обладают нулевыми моментами и, следовательно, игнорируют постоянное слагаемое, тем не менее по отношению к тренду они таким свойством обладают не всегда. Для изучения высокочастотных компонентов ряда обе эти составляющие полезно исключить из данных. 2. Исключение тренда и центрирование ряда Для исключения тренда необходимо задать его модель. Если природа тренда имеет теоретическое объяснение, то моделирование тренда производится на основе этой теории. Чаще всего мы не знаем природу тренда. В таких случаях в качестве формальной модели используют аппроксимацию тренда с помощью линейной комбинации каких-нибудь полиномов. При этом в состав такого выражения входит и свободный член. Параметры взятой модели тренда определяются с помощью метода наименьших квадратов, а затем значения тренда вычитаются из исходных данных. В простейшем случае такая операция сводится к исключению постоянного слагаемого (центрированию ряда). При этом среднее значение ряда определяется по формуле

m=

1 N

N -1

xk ,
k=0

а центрированный ряд получается из исходного следующим образом: x = xk - m, k = 0, 1, ..., N - 1. k 3. Графическое представление центрированного ряда

Для модельного ряда (74) это представление показано на рис. 6 (см. с 33).
24


4. Оценивание дисперсии ряда:
2 0 =

1 N -1

N -1

(x )2 . k
k=0

Перед выполнением вейвлет-анализа полезно провести предварительное изучение временного ряда с помощью Фурье-преобразования. Если это не требуется, можно сразу перейти к п. 7. 5. Вычисление периодограммы Здесь приводится упрощенный вариант Фурье-анализа нашего ряда без использования алгоритма быстрого преобразования Фурье. Сначала вычисляют дискретное преобразование Фурье
N -1

Xj =
k=0

x e k

-i

2 N

kj

,

j = 0 , 1 . . . , N - 1,

(76)

после чего вычисляется периодограмма

Dj =

1 [(Re Xj )2 + (Im Xj )2 ], N2

j = 0, 1 . . . , N/2

(при нечетном N вычисления ведут до (N + 1)/2). Отсчеты периодограммы соответствуют частотам

j = j,
где

j = 0, 1 . . . , N/2, 1 . N t

=

6. Графическое представление периодограммы и порога

обнаружения сигнала

Этот график позволяет отождествить значимые спектральные линии. Для этого следует задать уровень значимости q << 1 и определить соответствующий ему порог обнаружениям сигнала. Если нам не известны периоды, присутствующие в

25

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

Xq = - ln[1 -

N -2

(1 - q )2 ].

(77)

Если, наоборот, частота гармонического колебания нам известна, то

Xq = - ln q .

(78)

Все пики периодограммы, превышающие критический уро2 вень 0 Xq /N , считаются значимыми, т. е. принадлежат детерминированному компоненту ряда. Вероятность такого утверждения равна 1 - q .

Для модельного ряда эти результаты показаны на рис. 7 (см. с. 33).
7. Вычисление вейвлет-преобразования Дискретные значения амплитудной вейвлет-функции вычисляются по следующим формулам:

1 WA (ai , bj ) = n(ai , bj )
N -1

N -1

x k
k=0



tk - bj ai
2

,

(79)

n(ai , bj ) =
k=0

exp -

1 B

tk - bj ai

.

(80)

В этих формулах (t) принятый анализирующий вейвлет. Дальнейшее исполнение алгоритма основано на использовании вейвлета Морле с параметром 2 :

(t) = e

-t2 /

2

e

i 2 t

,

(81)

при этом в формуле (80) следует положить B = 2 .

26

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

da = 2t a,
где t радиус вейвлета, определяемый формулой (11).

(82)

Функция WA (a, b) определяет корреляцию между анализирующим вейвлетом, помещенным в точку b, и какой-либо деталью временного ряда размеров da с центром в точке b. Модуль этой функции принимает наибольшее значение в том случае, когда размер вейвлета совпадает с размером "текущей"детали временного ряда. В случае полигармонической функции естественной мерой масштаба ее деталей является период гармонических компонентов, мера же протяженности вейвлета da определяется значением масштабного коэффициента a. Для полигармонической функции, заданной на сетке с шагом t = const, диапазон периодов гармоник определяется величинами Pmin = 2t и Pmax = (N - 1) t. В соответствии с этим выберем наибольшее и наименьшее значения масштабного коэффициента из условия соответствия размеров вейвлета и предельных периодов гармонических компонентов:

2t amin = Pmin ,
откуда имеем

2t amax = Pmax ,

(83)

amin =

t , t

a

max

=

(N - 1) t . 2t

(84)

Для вейвлета Морле эти формулы переписываются следующим образом:

amin =

2t ,

a

max

=

(N - 1) t .

(85)

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


Введем теперь шаг дискретизации масштабов

a =

a

max

- amin , Na - 1

(86)

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

ai = a

min

+ a i,

i = 0, 1, ..., Na - 1.

(87)

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

bmin 0;

bmax (N - 1) t,

(88)

после чего дискретные значения сдвигов можно вычислить по формуле

bj = b

min

+ b j ,

j = 0, 1, ..., Nb - 1,

(89)

где Nb принятое число узлов сетки, следующих друг за другом с постоянном шагом

b =

bmax - bmin . Nb - 1

(90)

Легко понять, что при таком способе дискретизации параметра b вблизи границ bmin и bmax величины WA (ai , bj ) будут вычисляться с погрешностями, так как около границ невозможно использовать всю длину анализирующего вейвлета. Для исключения граничных эффектов следует вычислять вейвлетпреобразование только для значений сдвига, отстоящих от границ на величину, равную текущему радиусу вейвлета a t . При таком подходе вместо формулы (89) имеем

bj = b

min

+ b j ,

j = Ja + 1, ..., Nb - Ja - 1,

(91)

28


где Ja выраженный в единицах b радиус вейвлета, соответствующий текущему значению масштаба a: Ja =

ai . 2 b

(92)

В этой формуле квадратными скобками обозначена операция округления до ближайшего целого. Множество узлов дискретной сетки, определенное формулами (87) и (91), называется треугольником достоверности. Заметим, что очень часто краевые эффекты игнорируются и результаты вейвлет-анализа представляют просто в прямоугольной области узлов (87) и (89). 9. Вычисление скалограммы Значения скалограммы для всех принятых узлов сетки производится по формуле

S (ai , bj ) = |WA (ai , bj )|2 .
10. Визуализация скалограммы

(93)

Значения скалограммы S (ai , bj ) лежат на поверхности, для изображения которой используют два основных способа: 1) представление поверхности в трехмерном пространстве координат (a, b, S ). Графические пакеты MATHCAD, SURFER и др. позволяют построить такие графики и воспроизвести их при любой степени разворота и наклона системы координат (например, опция Surface пакета SURFER 32); 2) плоское изображение поверхности S (a, b) в координатах (a, b) в виде топографической карты. Масштаб третьего измерения указывается либо оцифровкой изолиний, либо раскраской областей между изолиниями. При этом градации цвета или последовательность различных цветов оцифровываются с помощью специальной цветовой палитры, являющейся приложением к графику вейвлета (например, опции Contour и Image пакета SURFER 32).

Представление скалограммы модельного ряда (74) в виде контурной карты показано на рис.9,a, в виде поверхности на рис.10,a (см. с. 34-35).
29


11. Вычисление скелетона Бывает так, что широкие контуры линий близких по частоте гармоник мешают проследить за эволюцией их частот во времени. Чтобы отсечь влияние контуров, можно выделить те точки скалограммы, в которых она имеет максимумы по переменным a и b:

Sij , если Si-1,j < Sij > Si+1,j или Si,j -1 < Sij > Si,j +1 , Sc (ai , bj ) = 0, в противном случае.

(94)

В этой формуле использовано обозначение Sij S (ai , bj ). Функцию (94) мы будем называть скелетоном. В случае синусоидального сигнала точки скелетона располагаются вдоль линий, идущих параллельно оси времени. Если в данных имеются гармонические или квазигармонические компоненты, то топографическая карта скелетона будет состоять из линий, ориентированных вдоль оси b. В случае шумового компонента линии скелетона вытягиваются в перпендикулярном направлении, т. е. параллельно оси a. Если в данных присутствуют и гармонические компоненты, и шум, то карта скелетона позволяет увидеть их раздельно.

Представление скелетона модельного ряда (74) в виде контурной карты показано на рис.9,b, а в виде поверхности на рис.10,b (см. с. 34-35).
12. Выделение сигнала из шума Критерий выделения сигнала из шума (73) позволяет построить массив значений скалограммы, которые с заданной вероятностью 1 - q принадлежат сигнальному, а не шумовому компоненту:

Sq (ai , bj ) =

2 S (ai , bj ), S (ai , bj ) 0 Z (ai , bj ) Tq , 0,
2 S (ai , bj ) < 0 Z (ai , bj ) Tq ,

(95)

30


где

Z (ai , bj ) =

1 n2 (ai , bj )

N -1

exp -
k=0

2 tk - bj . 2 ai

(96)

Иногда бывает целесообразно выделять сигнал не из самой скалограммы, а из ее скелетона. В этом случае в правую часть формулы (95) вместо величины S (ai , bj ) следует подставить Sc (ai , bj ) 1

Представление "сигнальных"линий скелетона модельного ряда (74) в виде поверхности показано на рис.9,c, а в виде контурной карты на рис.10,c (см. c. 34-35). Интересно, что проведенное нами восстановление сигнала из шума оказалось неполным: линия синусоиды с заданным периодом P = 10 c получилась слегка изогнутой. По-видимому, эта деформация произошла за счет взаимодействия нашей синусоиды с шумом.
13. Вычисление скейлограммы Эта функция вычисляется по формуле
Nb -Ja -1

1 Gi = Nb - 2J

a

S (ai , bj ),
j =Ja

i = 0, ...Na - 1,

(97)

где Ja задается формулой (92), если вычисления ведутся в треугольнике достоверности. В противном случае следует по ложить Ja = 0.

График скейлограммы модельного ряда (74) показан на рис. 8 (см. с. 33). Здесь же штриховой линией прочерчен порог об2 наружения сигнала 0 Z (a, b) Tq как функция параметра a при значении сдвига b, соответствующего середине ряда.
14. Конец алгоритма.
1 Следует отметить, что выделение сигнала из шума, описанное в этом пункте, можно применять лишь в тех случаях, когда сигнал засорен именно белым (на практике широкополосным) шумом. Такой подход более соответствует Фурье-анализу, чем вейвлетам, поскольку смысл вейвлет-анализа заключается скорее в изучении нерегулярностей, чем в их исключении.

31

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

Синусоида
Рис. 11 дает возможность сравнить две скалограммы синусоиды с периодом 20 с, заданной на равномерной сетке с шагом t = 1 с при N = 100. Мы видим, что использование MHAT-вейвлета дает специфическую скалограмму, состоящую из выпуклых поверхностей, следующих друг за другом с интервалом 10 с. Объясняется это тем, что с точки зрения MHAT-вейвлета (который сам имеет форму шляпы, пусть и мексиканской) характерной деталью на графике полигармонической функции является выпуклость или вогнутость синусоиды, т. е. полупериод колебания. В противоположность этому вейвлет Морле, имеющий форму гармонического колебания (хотя и затухающего), распознает во временном ряду именно колебания полного периода. Этим и объясняется характерная форма синусоиды в координатах (a, b, S ), похожая на равномерную горную гряду (см. рис. 11,e). Такое представление скалограммы синусоиды соответствует ее Фурье-образу, поэтому вейвлет Морле более предпочтителен для анализа гармонических функций, чем MHAT-вейвлет.

Разрешение двух синусоид
Следующие два примера посвящены вейвлет-анализу временного ряда, представляющего собой сумму двух синусоид с периодами P1 = 20 с и P2 = 35 с при t = 1 c и N = 100. На рис.12 показаны график исходной функции, ее Фурье-периодограмма, скалограммы и скейлограммы, полученные с помощью MHAT-вейвлета и вейвлета Морле. Легко заметить (см. рис. 12,c,d), что с помощью MHAT-преобразования не удается разрешить два компонента исходного ряда, несмотря на то, что их периоды существенно различны. Этот факт объясняется сильной локализацией MHAT-вейвлета во временной области, в результате чего разрешающая способность

32


Рис. 6: Модельный ряд (74).

Рис. 7: Периодограмма модельного ряда (74).

Штриховая линия 99-процентный порог обнаружения сигнала в шумах.

Рис. 8: Скейлограмма модельного ряда (74). Штриховой линией показан
99-процентный порог обнаружения сигнала в шумах.

33


Рис. 9: Результаты вейвлет-анализа модельного ряда (74), представленные в виде топографических карт. a скалограмма; b скелетон скалограммы; c сигнальный компонент скелетона, выделенный из шума с вероятностью 0.99. Вычисления произведены в треугольнике достоверности.

34


Рис. 10: Результаты вейвлет-анализа модельного ряда (74), представленные в виде поверхностей в трехмерном пространстве. a скалограмма; b скелетон скалограммы; c сигнальный компонент скелетона, выделенный из шума с вероятностью 0.99. Вычисления произведены во всей области определения параметров a и b.

35


36


Рис. 11: Вейвлет-анализ синусоиды: a исходный ряд; b Фурьепериодограмма; c скалограмма (MHAT-вейвлет); d скейлограмма (MHAT-вейвлет); e скалограмма (вейвлет Морле); f скейлограмма (вейвлет Морле).

в области масштабов становится очень низкой. С этой точки зрения результаты вейвлет-анализа нашей синусоиды с помощью вейвлета Морле оказываются более благоприятными. Теперь мы видим, что обе линии успешно разрешились как в скалограмме, так и в скейлограмме (см. рис. 12,e,f ). Объясняется этот факт просто: вейвлет Морле не так сильно локализован во времени (он всегда охватывает периодов), в результате чего при 2 = 2 его разрешающая способность в области масштабов больше, чем у MHAT-вейвлета. Тем не менее, могут возникнуть ситуации, когда две близкие по частоте гармоники не будут разрешаться и с помощью вейвлета Морле при стандартном значении = 2. В таких случаях значение параметра нужно увеличить. При этом мы, естественно, уменьшим временное разрешение вейвлет-преобразования за счет увеличения его разрешающей способности в области масштабов (периодов).

37


38


Рис. 12: Разрешение двух синусоид: a исходный ряд; b Фурьепериодограмма; c скалограмма (MHAT-вейвлет); d скалограмма (вейвлет Морле); e скейлограмма (MHAT-вейвлет); f скейлограмма (вейвлет Морле).

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

39


Рис. 13: Вейвлет-анализ синусоиды с переменным периодом: a
исходный ряд; b Фурье-периодограмма; c скалограмма (вейвлет Морле); d скейлограмма (вейвлет Морле).

40


Рис. 14: Вейвлет-анализ шумового сигнала: a дискретный белый шум;
b скалограмма (вейвлет Морле); c скелетон скалограммы (локализация максимумов вдоль шкалы масштабов; d скелетон скалограммы (локализация максимумов вдоль оси времени.)

41


Дискретный белый шум
На рис. 14 показан результат вейвлет-анализа последовательности 100 случайных чисел, полученных с помощью одной из стандартных программ генерирования значений нормально распределенной случайной величины, имеющей нулевое математическое ожидание и единичную дисперсию. На рис. 14,a представлен временной ряд, а на рис. 14,b картина скалограммы, имеющей топографию сложного горного массива. Рис. 14,c,d показывают линии скелетонов в направлениях осей времени (сдвига) и периодов (масштабов). Можно видеть, что в состав псевдослучайной последовательности входят как квазигармонические компоненты, эволюционирующие во времени (линии, идущие по направлению оси времени), так и чисто шумовые компоненты, т. е. линии, идущие перпендикулярно оси времени, на которых значения скалограммы остаются приблизительно одинаковыми в широком диапазоне масштабов.

Авторегрессионная модель
Для изучения многих процессов, промежуточных по своей природе между периодическими и чисто случайными, используются так называемые авторегрессионые модели AR(p), задаваемые в дискретном случае с помощью соотношения

xk = a1 x

k-1

+ a2 xk

-2

+ ... + ap xk

-p

+ k

k = p + 1, ... ,

(98)

где a1 , a2 ..., ap параметры модели, p порядок модели, k гауссовский белый шум с нулевым средним и единичной дисперсией. Например, модель AR(0) это просто шум, AR(1) марковский процесс, AR(2) псевдоколебания и т.д. Последовательности типа (98) характеризуются спектрами, в которых имеются концентрации мощности в узких интервалах частот. Следует отметить, что сами AR-модели не вскрывают внутреннюю структуру временного ряда. Очевидно, что такую возможность нам дает вейвлет-анализ. На рис. 15 и 16 показаны результаты Фурье-анализа и вейвлет-анализа авторегрессионных модельных рядов AR(2) и AR(3):

xk = xk xk = xk
-1

-1

- 0.4xk
-2

-2

+

k

k = p + 1, ...,
k

(99) (100)

- 0.5xk

+ 0.4x

k-3

+

k = p + 1, ...

42


Рис. 15: Анализ авторегрессионного процесса (99): a исходный ряд; b
периодограмма; с скелетон; d скейлограмма.

43


Рис. 16: Вейвлет-анализ авторегрессионного процесса (100): a исходный ряд; b скелетон.

На периодограмме Шустера ряда (99) мы видим концентрацию мощности в пределах от 0 до 0.10.2 Гц, однако в этом диапазоне линии плохо разрешаются. Это свидетельствует о существовании в структуре ряда квазипериодических компонентов. Скелетоны рядов (99) и (100) подтверждают это предположение, показывая одновременно, как эти псевдогармонические компоненты развиваются во времени. В скелетонах наших рядов имеются длинные и короткие синусоподобные линии с переменными интенсивностями и масштабами (периодами). Интересно отметить, что на картах скелетонов практически нет линий, направленных перпендикулярно оси времени. Основные линии скелетонов ориентированы вдоль оси времени, свидетельствуя о преобладании в AR-процессах скорее колебательных, чем стохастических компонентов.

44


Анализ чисел Вольфа
Для описания пятнообразовательной деятельности Солнца используют числа Вольфа W = k (10 g + f ), (101) где g число групп пятен, видимых на диске Солнца, f общее число пятен, как отдельных, так и принадлежащих группам, k коэффициент для приведения результатов наблюдений в единую систему. С физической точки зрения числа Вольфа не являются исчерпывающей характеристикой солнечной активности. Тем не менее, только они дают нам астрономическую информацию о вариациях солнечной активности в прошлом, поскольку обработка зарисовок диска Солнца позволила составить равномерный ряд среднегодовых чисел Вольфа, начиная с 1700 г., а отдельные значения этого ряда имеются вплоть до начала 17 столетия, когда солнечные пятна были открыты Галилеем. Ряд чисел Вольфа это знаменитый временной ряд. Во-первых, он играет ключевую роль в проблеме солнечно-земных связей, вовторых, именно он использовался Шустером (1906) и Юлом (1928) для построения фундаментальных понятий анализа временных рядов периодограммы и модели авторегрессии. Не претендуя на исчерпывающий анализ физических причин солнечной активности, используем ряд чисел Вольфа для демонстрации возможностей вейвлет-преобразования. На рис. 17,a показаны среднегодовые значения чисел Вольфа на промежутке от 1700 до 1999 г. Изучение этого ряда во временной области показывает много характерных особенностей. Во-первых, четко видна повторяемость максимумов и минимумов с характерным периодом приблизительно 11 лет (более тонкое изучение показывает, что расстояния между максимумами не остаются постоянными, а претерпевают небольшие изменения во времени). Во-вторых, легко усматривается различие максимальных значений ряда (так, наибольшее значение приходится на 1957 г.). Все это говорит о том, что характеристики ряда меняются во времени и что Фурье-анализ не является адекватным методом для его исследования. Действительно, в периодограмме Шустера (рис. 17,b) нашего ряда, освобожденного от линейного

45


Рис. 17: Анализ чисел Вольфа: a ряд среднегодовых чисел Вольфа 1700
1999; b периодограмма Шустера (cплошная линия), сглаженная периодограмма (штриховая линия), 99-процентный порог обнаружения сигнала в шумах (штрих-пунктирная линия ); c скейлограмма в диапазоне 5-120 лет; d скейлограмма в диапазоне 5-15 лет.

46


Рис. 18: Вейвлет-анализ ряда чисел Вольфа: a cкелетон в диапазоне
периодов (масштабов) 5-120 лет; b скалограмма в диапазоне периодов 5-15 лет; c скелетон в диапазоне периодов 5-15 лет.

47


тренда

tr(t) = (33.9 + 2.3) + (0.105 + 0.013)(t - 1700),

(102)

мы видим две значимые концентрации мощности, приходящиеся на частоты 0.00-0.02 и 0.08-0.10 цикла в год. В первой полосе есть пики, соответствующие периодам 204.5; 102.3 и 56.8 года. Во второй полосе имеются четыре максимума с периодами 12.04; 11.00; 10.55 и 10.03 года. В обеих полосах частот периодограмма сильно изрезана. Изрезанность периодограммы свидетельствует о значительной стохастичности, присущей процессу солнечной активности. В подобных случаях в спектральном анализе прибегают к сглаживанию периодограммы. Результаты сглаживания с помощью окна Тьюки при параметрах сглаживания N = 150 и a = 0.25 показаны штриховой линией на рис. 17,b. Мы видим, что и периодограмма Шустера, и ее сглаженная модификация говорят о том, что ряд солнечной активности состоит из трех квазипериодических компонентов с периодами приблизительно 100, 57 и 11 лет, и это все, что Фурье-анализ может дать при исследовании этого ряда. В противоположность этому вейвлет-анализ позволяет увидеть не только концентрации мощности на известных масштабах, но и проследить за их развитием во времени. На рис. 18,a представлена скалограмма, вычисленная с вейвлетом Морле в диапазоне масштабов от 5 до 120 лет. Здесь имеются три спектральные линии, соответствующие масштабам (периодам) 100, 54 и 11 годам, однако, в отличие от Фурье-спектра, мы видим, что и периоды, и амплитуды этих линий изменяются во времени. На рис.18,b показан только участок скалограммы в диапазоне масштабов от 5 до 15 лет. Здесь эти изменения видны особенно ясно. Характерной особенностью является резкое падение интенсивности пятнообразования с 1800 по 1830 г., сопровождающееся одновременным изменением периода. Еще более четко изменение периода во времени показано на рис. 18,c, где приводится не контуры спектральных линий, а только линии скелетона. Мы видим, что основной, 11-летний цикл Солнца, имеющий вид извилистой линии, идущей вдоль оси времени, пересекается шумовыми полосами, вытянутыми в перпендикулярном направлении. Этот факт говорит о том, что пятнообразовательная деятельность Солнца характеризуется не только периодическим механизмом с переменными периодом и амплитудой, но и аддитивными стохастическими компонентами типа белого шума. 48


Кривая блеска квазара 3C273
Объект 3С273 был идентифицирован как квазар в 1963 г. Это самый близкий и самый яркий из известных квазаров. На рис.19,a показана его кривая блеска в полосе V на интервале времени от 1887 до 1997 г. с шагом 100 сут. Первая часть (1887.57 1967.24) этой кривой была составлена Кункелем (1967) по фотографическим и фотоэлектрическим наблюдениям. Вторая часть (1968.09 1997.27) основана на базе данных (Терлер и др., 1999; http:// obswww.unige.ch/3c273/), с помощью которых блеск в полосе V был вычислен по формуле

V = 16.5 - 2.5 lg(F ),

(103)

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

tr(t) = (12.737 + 0.004) + (1.852 + 0.004)10-3 (t - 1900)

(104)

мы вычислили периодограмму Шустера (рис. 19,b), на которой имеется 7 пиков, превышающих 99-процентный порог обнаружения сигнала в шумах. Эти пики могут свидетельствовать о существовании в блеске квазара 3С273 гармонических компонентов с периодами 56.0, 33.0, 25.5, 17.0, 14.4, 12.5 и 11.2 года. Однако сильная изрезанность периодограммы свидетельствует, как правило, о нестабильности этих гармоник. Результаты сглаживания периодограммы методом Тьюки (N = 200; a = 0.25) показаны штриховой линией на рис. 19,b. Таким образом, анализ Фурье позволяет сделать заключение о том, что кривая блеска квазара 3С273 по своей природе представляет собой узкополосный шум в полосе частот от 0 до 0.1 цикла в год. На рис. 20,a показаны линии скелетона в диапазоне периодов от 1 до 65 лет. Основные ее детали сконцентрированы в области периодов от 1 до 35 лет, которую мы показываем отдельно в виде скалограммы (рис. 20,b) и в виде карты ее скелетона (рис. 20,с). Мы видим, что линии скелетона ориентированы преимущественно вдоль оси времени, что говорит в пользу существования в кри49


Рис. 19: Анализ блеска квазара 3С273: a кривая блеска квазара 3С273
(сплошная линия) и ее тренд (пунктрная линия); b периодограмма Шустера (сплошная линия), сглаженная периодограмма (штриховая линия), 99-процентный порог обнаружения сигнала в шумах (штрихпунктирная линия); c скейлограмма.

50


Рис. 20: Вейвлет-анализ кривой блеска квазара 3С273: a cкелетон в
диапазоне периодов (масштабов) 1-65 лет; b скалограмма в диапазоне периодов 1-35 лет; c скелетон в диапазоне периодов 1-35 лет;

51


вой блеска квазара 3C273 синусоподобных компонентов. Более того, значения всех периодов, полученных по периодограмме Шустера, подтверждаются скалограммой. Это обстоятельство является очень важным, поскольку в нашем случае можно сказать, что изрезанность периодограммы объясняется не тем, что она не является состоятельной оценкой спектральной плотности (это верно только для белого шума), а тем, что периодограмма кривой блеска квазара 3С273 на самом деле регистрирует квазипериодические сигналы, которые возникают и исчезают случайным образом. Этот факт говорит о том, что периодограмма Шустера и скалограмма дополняют друг друга: первая обнаруживает в данных гармонические или квазигармонические компоненты с высоким частотным, но с нулевым временным разрешением, а вторая позволяет локализовать гармоники во времени (с плохим разрешением по масштабу). В астрофизической литературе долго обсуждался вопрос о физической модели кривой блеска квазара 3С273. Одни авторы объясняли кривую блеска суперпозицией гармонических колебаний, другие стохастическими эффектами типа авторегрессионных моделей. Эти альтернативы (периодические колебания или стохастичесие вспышки) были основаны на специфике применявшихся методов преобразовании Фурье и авторегрессионного анализа. Мы видим, что специфика вейвлет-анализа в известной степени объединяет эти два подхода, позволяя объяснить кривую блеска квазара 3С273 с помощью механизма стохастических периодичностей, т. е. квазипериодических процессов, возникающих и исчезающих случайным образом. Интересно также отметить, что отсутствие на рис. 20,с линий, ориентированных вдоль оси периодов, говорит о том, что в кривой блеска квазара 3C273 отсутствует аддитивный шумовой компонент, т. е. процесс стохастических колебаний является основным механизмом, объясняющим характер изменения блеска этого квазара.

52


Упражнения
1. Составьте программу анализа временных рядов с помощью вейвлетпреобразования, пользуясь алгоритмом, приведенным на с. 23-31. 2. Модифицируйте программу так, чтобы при вычислении вейвлетпреобразования перейти от масштаба a к частоте = 2 /a. Используя эти программы, проведите анализ следующих модельных временных рядов.

ћ Скачок периода
Этот пример является своеобразной визитной карточкой вейвлет-преобразования. Графические результаты, которые Вы получите, выполняя это упражнение, можно встретить во многих руководствах по вейвлетам. Смоделируйте для четного N два следующих ряда:

xk = A1 cos(2 1 tk - 1 ) + A1 cos(2 2 tk - 2 ), A1 cos(2 1 tk - 1 ), k N/2, A1 cos(2 2 tk - 2 ), k > N/2, k = 0, 1, 2, . . . , N - 1.

(105)

yk =

(106) (107)

tk = t k ,

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

ћ Скачок амплитуды
Для четного значения N cмоделируйте два ряда

xk = A1 cos(2 1 tk - 1 ), A1 cos(2 1 tk - 1 ), k N/2, A2 cos(2 1 tk - 2 ), k > N/2, k = 0, 1, 2, . . . , N - 1.

(108)

yk =

(109) (110)

tk = t k ,

Объясните различие периодограмм и скалограмм двух рядов при A1 = A2 . 53


ћ Дрейф частоты
Получите модельный ряд

xk = A cos

2 tk Pk

,

k = 0, 1, 2, . . . , N - 1,

(111)

где A амплитуда, Pk переменный период:

Pk = P 1 +

1 sin 2

2 tk P0

,

k = 0, 1, 2, . . . , N - 1. (112)

Здесь P начальный период гармонического компонента; P0 период изменения величины P . Постройте скалограмму при следующих значениях параметров:

N = 100,

t = 1 c, A = 1, P = 15 c, P0 = 200 c.

Почему ширина линии гармонического компонента с переменным периодом меняется с изменением значения параметра b?

ћ Синусоиды в белом шуме
Получите модельный ряд

xk = A1 cos(2 1 tk - 1 ) + A2 cos(2 2 tk - 2 ) + n nk , (113) tk = t k ,
где

k = 0, 1, 2, . . . , N - 1,

(114)

t постоянный шаг выборки; Ai , i , i , i = 1, 2, амплитуды, частоты и фазы двух гармонических компонентов; nk , k = 0, 1, 2, . . . , N - 1, значения случайной величины с нулевым математическим ожиданием и единичной дисперсией, распределенной по нормальному закону (шумовой компонент ряда); n среднеквадратическое отклонение шумового компонента nk , которое можно вычислить через отношение "сигнал к шуму" с помощью соотношения
54


n =
Задаваемые параметры ряда:

A2 + A2 1 2 . 2

(115)

1) характеристики гармоник Ai , i , i , i = 1, 2; 2) отношение "сигнал к шуму" . Выполните следующие упражнения и ответьте на заданные вопросы. 1. Получите скалограммы ряда (113) при N = 100 и N = 1000. Почему разрешающая способность вейвлет-преобразования не зависит от длины ряда? 2. Как обеспечить разрешение в скалограмме двух близких по частоте гармоник (скажем, 1 = 10 Гц и 2 = 12 Гц при t = 1 с)? 3. Используя программу вейвлет-анализа, постройте скалограммы ряда, состоящего из двух близких по частоте гармоник, так, чтобы увидеть эффекты слияния линий, а также картины их частичного и полного разрешения.

ћ Линейный тренд
Получите модельный ряд (116) (117)

xk = a + b tk , tk = t k , k = 0, 1, 2, . . . , N - 1,

где t постоянный шаг выборки; a, b параметры тренда. Выполните следующие упражнения и ответьте на заданные вопросы. 1. Постройте скалограммы и скелетоны линейного тренда при a = 0 и b = 0. 2. Постройте скалограмму ряда (116) при b = 0. В каких случаях эффект постоянного слагаемого обнаруживается в скалограмме вопреки соотношению (13)?

55


ћ Вейвлет-анализ шума
Получите модельный ряд (118)

xk = k ,

k = 0, 1, 2, . . . , N - 1,

где k , k = 0, 1, 2, . . . , N - 1, выборка случайной величины с нулевым математическим ожиданием и единичной дисперсией, распределенной по нормальному закону; Получите скалограммы нескольких реализаций шума, пользуясь стандартными программами генерации случайной величины, распределенной по нормальному закону. Почему в картах скелетонов всегда присутствуют линии, ориентированные и вдоль оси времени, и вдоль оси масштабов?

ћ Псевдоколебания
Постройте следующие модельные ряды:

xk = 1.372 x xk = 1.246 x
где

k-1 k-1

- 0.693 xk

-2

+ 16.55 k , k = 2, 3, ...N - 1, (119) + 0.052
k

- 0.344 xk

-2

k = 2, 3, ...N - 1, (120)

x0 = 0, x1 = 0 начальные значения; k , k = 0, 1, 2, . . . , N - 1, выборка из случайной величины с нулевым математическим ожиданием и единичной дисперсией, распределенной по нормальному закону.
Получите несколько реализаций рядов (119) и (120), пользуясь стандартными программами генерации нормально распределенной случайной величины, и постройте их скелетоны. Ряды (119) и (120) это авторегрессионные модели ряда чисел Вольфа и кривой блеска квазара 3С273 (освобожденные от линейных трендов). Сравните результаты, полученные Вами для моделей (119) и (120), с характеристиками вейвлетанализа исходных рядов, которые приведены в двух предыдущих параграфах.

56


Литература
Гроссман и Морле, 1984. Grossman A., Morlet J. Decomposition of Hardy functions into square integrable wavelets of constant shape // SIAM J. Math. P.723-736. Гэбор,1946. Gabor D. // J. Inst. Elect. Eng. Vol. 93. P. 429. Добечи, 1992. Daubechies I. Ten lectures on wavelets. Society for industrial and applied mathematics. Philadelphia, Pennsylvania, 1992. Койфман, 1992. Wavelets and their applications /Ed. R.Coifman. Boston: Jones and Barlett Publ. Кункель, 1967. Kunkel W.E. A harmonic analysis of the Light Variations of 3C273. // Astron. J. Vol. 72. N10. P. 1341-1348. Скаргл и др., 1993. Scargle J.D. The Quasi-Periodic Oscillations and Very Low Frequency Noise of Scorpius X-1 as Transient Chaos: A Dripping Handrail? // Ap. J. Vol. 411. L91-L94. Скаргл, 1997. Scargle J.D. Wavelet and Other Multi-resolution Methods for Time Series Analysis. Statistical Challenges in Modern Astronomy II /Ed. G.J.Babu and E.D.Feigelson. P. 333-347. N.Y.: Springer-Verlag. Терлер, 1999. Turler et al. 30 years of multi-wavelength observations of 3C273 // A&AS. Vol. 134. P. 89-101. Тримбл В. 1997. Trimble V. Late-Night Thoughts of a Classical Astronomer. Statistical Challenges in Modern Astronomy II /Ed. G.J.Babu and E.D.Feigelson. P. 365-385. N.Y.: Springer-Verlag. Фостер, 1996. Foster G. Wavelets for period analysis of unevenly sampled time series // Astron. J. Vol. 112. N4. P. 1709-1729. Шустер, 1906. Schuster A. On the Periodicities of Sun Spots // Trans. R. Soc. London. Ser A. Vol. 206. P. 69-100. Юл, 1928. Yule J.U. On a Method of Investigation Periodicities in Disturbed Series with Special Reference to Wolfer's Sunspot Numbers. // Phylos. Trans. R. Soc. London. Ser. A. Vol. 226. P.267-298.

57


.

Содержание
Введение От Фурье-анализа к вейвлетам Основы вейвлет-анализа
Свойства вейвлетов . . . . . . . . . . . . Примеры вейвлетов . . . . . . . . . . . . Вейвлет-преобразование как фильтрация Спектральные характеристики . . . . . ... ... .. ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3 4
10 11 13 15

7

Практические аспекты вейвлет-анализа

Скалограмма и скейлограмма . . . . . . . . . . . . . . . . . Вейвлет-преобразование синусоиды . . . . . . . . . . . . . . Вейвлет-преобразование шума . . . . . . . . . . . . . . . . .

16

16 18 20

Алгоритм вейвлет-анализа временных рядов Вейвлет-анализ модельных функций
Синусоида . . . . . . . . . . . . . . Разрешение двух синусоид . . . . . Синусоида с переменным периодом Дискретный белый шум . . . . . . Авторегрессионная модель . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23 32
32 32 42 42 42

Анализ чисел Вольфа Кривая блеска квазара 3C273 Упражнения Литература

46 50 54 58


Учебное издание

Вениамин Владимирович Витязев
ВЕЙВЛЕТ-АНАЛИЗ ВРЕМЕННЫХ РЯДОВ Учебное пособие

Зав. редакцией Г.И. Чередниченко Обложка Е.А.Соловьевой

Лицензия ЛР N 040050 от 15.08.1996 Подписано к печати с оригинал-макета 03.05.2001. Печать офсетная. Ф-т 60 Ч 84/16. Усл. печ. л. 3,49. Уч.-изд. л. 3,44. Тираж 50 экз. Заказ N Редакция оперативной подготовки учебно-методических и научных изданий Издательства С.-Петербургского университета. 199034, С.-Петербург, Университетская наб., 7/9. ЦОП типографии Издательства СПбГУ. 199034, С.-Петербург, наб. Макарова, 6.


. Для заметок =======================================