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

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

В.В. Витязев

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

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

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

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

В54

Витязев В.В.

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

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


.

2


.

Содержание
Основные определения
Периодограмма и коррелограмма . . . . . . Периодограмма гармонического колебания Периодограмма дискретного белого шума . Сглаживание периодограммы . . . . . . . . Взаимные характеристики двух временных Дискретное преобразование Фурье . . Быстрое преобразование Фурье . . . . Дополнение ряда нулями . . . . . . . . Вычисление коррелограмм с помощью Спектрально-Корреляционный ( СКАВРя) . . . . . . . Взаимный Анализ Временных Упражнения . . . . . . . . . . .... .... .... .... рядов .... .... .... ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 8 10 12 14

4

Вычислительные процедуры

... ... ... БПФ

19
19 21 24 25

Алгоритмы

Анализ Временных ............ Рядов (ВАВРя) . . ............

Рядов ..... ..... .....

26
26 32 40

Приложение

Вариации продолжительности суток . . . . . . . . . . . . . Движение полюса в теле Земли . . . . . . . . . . . . . . . .

42

43 45

Литература

47

3


Основные определения
Спектральный анализ временных рядов это один из методов обработки результатов экспериментов и, в частности, данных астрономических наблюдений. Теория спектрального анализа основана на равносильности представления функций во временной и частотной областях с помощью преобразований Фурье. Это обстоятельство позволяет построить эффективные методы получения характеристик исследуемого сигнала в частотной области в тех случаях, когда во временной области это сделать трудно. В настоящее время методы спектрального анализа временных рядов широко применяются в различных областях знаний: в физике, астрономии, биологии, технике, экономике, медицине и т. д. Сделаем краткий обзор основных положений, лежащих в основе теории и практики спектрального анализа временных рядов. Обычно под временным рядом понимают упорядоченную во времени совокупность измерений определенной количественной характеристики какого-либо процесса. По способу записи результатов измерений временные ряды бывают непрерывными, дискретными, равномерными и неравномерными. Общим признаком временных рядов, легко распознаваемым при их простом обзоре, является наличие в них крупномасштабных изменений (тренд, периодические колебания) и мелкомасштабных компонентов, в структуре которых, как правило, отсутствуют какие-либо правильные изменения (рис. 1). Основной задачей анализа временных рядов является интерпретация наблюдательных данных с помощью некоторых моделей. В качестве моделей выбирают функции, свойства которых хорошо известны. Отмеченное выше наличие во временных рядах правильных и неправильных изменений оправдывает широко распространенную модель сигнал + шум. В рамках этой модели сигнал интерпретируется с помощью детерминированных функций: периодических, почти периодических (полигармонических) и непериодических, а шум с помощью особого вида случайных функций так называемого белого шума. В тех случаях, когда четкая грань между сигналом и шумом размыта, в качестве моделей временных рядов используются случайные процессы. Теоретическим фундаментом спектрального анализа временных рядов является теорема ВинераХинчина, которая устанавливает 4


Рис. 1: Временной ряд вариаций периода вращения Земли с 1962 по
1996 г. по данным Международной службы вращения Земли (IERS).

связь между двумя характеристиками случайного процесса плотностью спектра мощности и корреляционной функцией. Сама теорема формулируется следующим образом: Если центрированный стационарный случайный процесс задан на промежутке 0 t T множеством реализаций

X (t) = {xj (t)}j

N =1

,

то его корреляционная функция
K (t, t ) = E [xj (t)xj (t )] = C ( ), =t-t,
(1)

и плотность спектра мощности (в дальнейшем спектр мощности) 2 T 1 S ( ) = E lim xj (t)e-it dt (2) T 2 T 0 взаимно обратимы по Фурье:
S ( ) = 1 2
+

C ( )e
- +

-i

d .

(3)

C ( ) =
-

S ( )ei d .

(4)

5


В формулах (1) и (2) символом E обозначена операция нахождения математического ожидания. Доказательство этой теоремы можно найти, например, в монографии Дженкинса и Ваттса (1972). Приведем два примера. N 1. Для случайной функции X (t) = {A cos( t + j )}j =1 , у которой фазы равномерно распределены в интервале [0, 2 ], корреляционная функция и спектр мощности задаются выражениями

C ( ) = S ( ) =

A2 cos( ), 2

(5) (6)

A2 [ ( + 0 ) + ( - 0 )], 4 где ( ) дельта-функция Дирака.
2. Для белого шума соответственно имеем

C ( ) = ( ), S ( ) = const.

(7) (8)

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


во времени, приходится применять более сложные методы их исследования. Предположим, что центрированный эргодический случайный процесс x = x(t) задан конечным набором вещественных значений xk = x(tk ) на равномерной сетке tk = t k , k = 0, 1, . . . , N - 1, где постоянная величина t является шагом выборки. В этом случае в качестве оценки спектра мощности обычно берут периодограмму Шустера, определенную следующим образом:

1 D ( ) = 2 N

N -1

2

xk e
k=0

-i2 t

k

.

(9)

Заметим, что при переходе от формулы (2) к формуле (9) мы заменили круговые частоты обычными частотами , связанными между собой очевидным соотношением = 2 . Вычисления по формуле (9) можно производить для любого значения частоты , однако особый интерес представляют отсчеты периодограммы, отнесенные к так называемой фундаментальной системе частот (ФСЧ):

2c j j= , j = 0, 1, ..., N/2, N N t где c частота Найквиста, определенная равенством j = c =

(10)

1 (11) 2t. В формуле (10) мы предположили, что N четное число. Для нечетных N конечное значение индекса j следует заменить числом (N + 1)/2. С учетом этого замечания в дальнейшем мы везде будем считать, что N четное число.
Для оценки корреляционной функции вводят коррелограмму
N -m-1

cm c(t m) =

1 N -m

xk xk
k=0

+m

,

(12)

m = 0, 1, . . . , N - 1.
Отметим, что периодограмма и коррелограмма связаны между собой соотношением

D( ) =

1 2Re N

N -1

cm e ?
m=0

-i2 mt

- c0 , ?

(13)

7


Рис. 2: Различие между спектром мощности и периодограммой.
спектр мощности косинусоиды; b периодограмма косинусоиды.

a

где смещенная оценка корреляционной функции cm определяется ? следующим образом:

cm = ?

1 N -m cm = N N

N -m-1

xk x
k=0

k+m

,

(14)

m = 0, 1, . . . , N - 1.

Периодограмма гармонического колебания
Перечислим сначала свойства периодограммы гармонической функции x(t) = A cos(2 0 t + ). (15) Спектр мощности этой функции представляет собой две бесконечно узкие спектральные линии на частотах = +0 : 8


A2 [ ( + 0 ) + ( - 0 )]. (16) 4 Можно показать, что в случае задания функции (15) на дискретном множестве точек с постоянным шагом t ее периодограмма будет иметь следующее представление: S ( ) = A2 [W ( + 0 ) + W ( - 0 )], 4 где спектральное окно W ( ) задается выражением D( ) W ( ) = sin2 N t . N 2 sin2 t
(17)

(18)

Сопоставление спектра мощности и периодограммы косинусоидальной функции показано на рис. 2. Отметим некоторые свойства спектрального окна. Во-первых, эта функция периодична с основным периодом 2c , где c частота Найквиста, задаваемая выражением (11). Из периодичности окна следует, что периодограмма дискретного ряда является периодической функцией с тем же основным периодом. По этой причине значения периодограммы имеет смысл вычислять только на основном периоде [-c , c ]. К тому же периодограмма вещественной функции обладает свойством четности, поэтому обычно значения периодограммы вычисляют только в интервале [0, c ]. Можно показать, что ширина главного лепестка спектрального окна (18) обратно пропорциональна длине ряда

2 . (19) N t Следует иметь в виду, что на равномерной сетке с постоянным шагом t имеет место соотношение = cos(2 0 t + ) = cos[2 (0 + 2pc ) t + ], p = +1, +2, ...,

в силу которого периодограмма в основном диапазоне частот будет иметь один и тот же пик на частоте 0 для всех гармонических колебаний с чаcтотами 0 + 2pc . В специальной литературе явление неразличимости указанных частот называют элайзингом (от английского слова alias кличка преступника).

9

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

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

xk = x(tk ),

tk = t k ,

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

(20)

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

1 Dj = 2 N

N -1

2

xk e
k=0

-i

2 N

kj

,

(21)

отсчеты которой соответствуют фундаментальной системе частот (10), и введем теперь в рассмотрение нормированную периодограмму dj = N Dj .

Теорема (Шустер). Отсчеты dj нормированной периодограммы представляют собой случайную величину с плотностью распределения
p(x) = e
-x/2

e-x ,

j = 1, 2, ..., N/2 - 1,
(22)

/ 2 x, j = 0; j = N/2, 0 < x < .

10


Рис. 3: Периодограмма шума. Эта теорема является основой построения критериев выделения сигнала из шума. Шумовой компонент, имеющийся в измерениях, проявляет себя в периодограмме в виде многочисленных пиков (рис. 3). При сильной зашумленности ряда возникают ситуации, когда пики детерминированного компонента очень трудно отличить от пиков шумовой составляющей. В таких случаях применяют статистические критерии разделения сигнального и шумового компонентов. В общем случае, когда центрированный дискретный белый шум 2 имеет не единичную, а произвольную дисперсию 0 , критерий распознавания сигнала в шуме используют следующим образом. Задают положительное число q << 1, которое определяет вероятность (уровень значимости) очень редкого события появления сильного пика в периодограмме белого шума. При этом имеют место две ситуации. 1. Частота периодического компонента заранее не известна. В этом случае вполне естественно предположить, что наибольший отсчет периодограммы Dmax соответствует искомому периодическому компоненту. Если выполняется неравенство

Dmax
где

2 0 X1 , N

(23)

11


X1 = - ln(1 - (1 - q )

2 N -2

),

(24)

то с вероятностью p = 1 - q принимается утверждение, что отсчет Dmax принадлежит сигналу, а не шуму. 2. Если частота искомого сигнала заранее известна и соответствует отсчету периодограммы Dl , то при выполнении неравенства

Dl
где

2 0 X2 , N

(25) (26)

X2 = - ln q ,

с вероятностью p = 1 - q полагают, что отсчет Dl тоже принадлежит сигналу, а не шуму. Следует отметить, что теорема об экспоненциальном распределении отсчетов периодограммы центрированного дискретного белого шума справедлива для того случая, когда дисперсия ряда известна точно. На практике ее приходится оценивать по формуле
2 0 =

1 N -1

N -1

(xk )2 .
k=0

(27)

В этом случае предельные отсчеты (24) следует заменить величинами, приведенными в табл. 1 (Теребиж, 1992). Следует, однако, заметить, что при N 100 табличные значения весьма мало отличаются от величин, получаемых по формуле (24).

Сглаживание периодограммы
Из теоремы Шустера следует, что периодограмма как оценка спектра мощности дискретного белого шума не удовлетворяет критерию эффективности. Этим объясняется изрезанность периодограммы, не уменьшающаяся с увеличением длины ряда. Тем не менее, осреднение периодограммы дискретного белого шума в узких диапазонах частот приводит к получению вполне разумной гладкой оценки. На этом свойстве основан метод сглаживания периодограммы, известный как метод Блэкмана-Тьюки (Блэкман и Тьюки, 1959). 12


Таблица 1: Отсчеты нормированной периодограммы центрированного дискретного белого шума в зависимости от N и q .

N
6 8 10 12 14 16 18 20 22 32 42 52 62 82 102 122

q = 0.05
1.950 2.613 3.072 3.419 3.697 3.928 4.125 4.297 4.450 5.019 5.408 5.701 5.935 6.295 6.567 6.758

q = 0.01
1.990 2.827 3.457 3.943 4.331 4.651 4.921 5.154 5.358 6.103 6.594 6.955 7.237 7.663 7.977 8.225

N
142 162 182 202 302 402 502 602 702 802 1002 1202 1402 1602 1802 2002

q = 0.05
6.967 7.122 7.258 7.378 7.832 8.147 8.389 8.584 8.748 8.889 9.123 9.313 9.473 9.612 9.733 9.842

q = 0.01
8.428 8.601 8.750 8.882 9.372 9.707 9.960 10.164 10.334 10.480 10.721 10.916 11.079 11.220 11.344 11.454

В этом методе на первом шаге вычисляются N N значений коррелограммы

cm = ?

1 N

N -m-1

xk x
k=0

k+m

,

m = 0, 1, . . . , N - 1.

(28)

Для получения сглаженной периодограммы коррелограмму умножают на некоторую функцию, называемую корреляционным окном. В практике спектрального анализа получила распространение весовая функция m Wm = (1 - 2a) + 2a cos , (29) N получившая название "окно Тьюки". Эта функция содержит два управляющих параметра 0 a 0.25 и N N , смысл которых будет объяснен ниже. Взвешенная коррелограмма вычисляется сле13


дующим образом:

cm = cm Wm , ~ ?

m = 0, 1, . . . , N - 1.

(30)

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

1 ~ Dj = 2Re N

N -1

cm e ~
m=0

-i

2 N

jm

- c0 . ~

(31)

Значения сглаженной периодограммы относятся к системе частот

j =

j 2c j= , N N t

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

(32)

которая, как это видно из сравнения с (10), уже не является фундаментальной системой частот. В соответствии с замечанием, сделанным на с. 7, в формуле (32) предельное значение индекса j назначается равным (N + 1)/2, если N число нечетное. Величина параметра N управляет степенью сглаживания периодограммы. Если исходный временной ряд представляет собой широкополосный шум, то рекомендуемое оптимальное сглаживание происходит при N = 0.1N . Для временного ряда, содержащего гармонические компоненты и шум, значение этого параметра выбирается из диапазона 0.1N N N . С помощью второго параметра окна Тьюки производят изменение контура спектральной линии гармонических компонентов исходного ряда с целью устранения боковых лепестков. При этом происходит расширение контура линии и уменьшение ее интенсивности. Рекомендуемые значения таковы: a = 0.23 окно Хэмминга, a = 0.25 так называемое окно "Хэннинг". Выбор значений параметров N и a производится из расчета удовлетворения двум противоречивым требованиям: получению статистически состоятельной периодограммы и достижению разумной разрешающей способности. Обычно этот компромисс достигается на основе физических соображений или априорной информации.

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

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

xk = x(tk ), yk = y (tk ),

tk = t k , tk = t k ,

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

(33) (34)

Обозначая символом операцию комплексного сопряжения, определим взаимную периодограмму рядов xk и yk с помощью соотношения 1 Dxy ( ) = 2 X ( )Y ( ), (35) N где
N -1

X ( ) =
k=0 N -1

xk e

-i2 t

k

,

(36)

Y ( ) =
k=0

yk e-

i 2 t

k

.

(37)

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

cxy [m] = ?

1 N

N -m-1

xk yk
k=0

+m

,

(38)

m = 0, 1, . . . , N - 1.
При этом имеет место свойство

cxy [-m] = cyx [m]. ? ?

(39)

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

Dxy ( ) =

1 N

N -1

N -1

cxy [m]e ?
m=0

-i2 mt

+
m=0

cyx [m]e ?

i2 mt

- c0 . ?

(40) Очевидно, что в том случае, когда наши два ряда xk и yk тождественны, выражения (35) и (38) переходят в (9) и (14), т. е. взаимные 15

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

Dxy ( ) = Pxy ( ) - iQxy ( ),
где

(41) (42) (43)

Pxy ( ) = Re X ( ) Re Y ( ) + Im X ( ) Im Y ( ), Qxy ( ) = Re X ( ) Im Y ( ) - Re Y ( ) Im X ( ).

Соотношение (42) определяется суммой произведений одноименных компонентов, т. е. произведениями действительной на действительную и мнимой на мнимую части функций X ( ) и Y ( ), в то время как соотношение (43) определяется суммой произведений разноименных компонентов, т. е. произведениями действительной на мнимую часть тех же функций. Поскольку на комплексной плоскости одноименные компоненты направлены вдоль одной оси, а разноименные компоненты направлены вдоль взамно перпендикулярных осей, функцию Pxy ( ) называют коспектром, а функцию Qxy ( ) квадратурным спектром. Для получения конкретной информации об общих свойствах исследуемых рядов вводятся в рассмотрение модуль взаимной периодограммы

Axy ( ) =
и фазовый спектр

2 Pxy ( ) + Q2 ( ) xy

(44)

xy ( ) = -arctg

Qxy ( ) . Pxy ( )

(45)

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

xk = A1 cos(2 1 tk + 1 ),
16

(46)


yk = A2 cos(2 2 tk + 2 ), tk = t k , k = 0, 1, ..., N - 1.

(47)

Нетрудно показать, что при N взаимная коррелограмма функций (46) и (47) равна нулю в том случае, когда 1 = 2 независимо от фаз этих гармоник. Иначе говоря, две гармонические функции с различными частотами некоррелированы. Предположим теперь, что выполняются соотношения 1 = 2 и 1 = 2 . В этом случае

cxy [m] = ?
где

A1 A2 cos(2 1 tm + ), 2 = 2 - 1 .

(48) (49)

На основании равенства (40) для значения взаимной периодограммы на частоте 1 = 2 получаем

A1 A2 i e, (50) 4 откуда для значений модуля взаимной периодограммы и фазового спектра на этой же частоте имеем Dxy (1 ) = Axy (1 ) = A1 A2 , 4
(51) (52)

xy (1 ) = .

Теперь предположим, что наши временные ряды (36) и (37) представляют собой выборки центрированных некоррелированных случайных величин, распределенных по нормальному закону c единичной дисперсией. Используя формулу (35), вычислим взаимную коррелограмму Dxy (j ) двух наших дискретных белых шумов на множестве фундаментальных частот (10), после чего найдем значения квадрата нормированного модуля взаимной периодограммы

axy (j ) = N 2 A2 (j ). xy

(53)

Теорема. Числа axy (j ) представляют собой значения случайной
величины с плотностью распределения 2K0 (2 x), j = 1, 2, ..., N/2 - 1, p(x) = K0 ( x)/ x, j = 0; j = N/2, 0 < x < ,
17 (54)


где K0 (z ) функция Бесселя мнимого аргумента:

K (z ) =

1 2

z 2

0



exp(-t - z 2 /4t) t

- -1

dt, = 0, 1, ...

(55)

Эта теорема подобно теореме Шустера (54) для обычной периодограммы позволяет построить статистические критерии обнару2 жения общих гармоник двух исследуемых рядов с дисперсиями x и 2 y . Новые критерии вполне аналогичны тем, которые используются при выделении сигнального компонента периодограммы. Действительно, зададим уровень значимости q << 1, который определяет вероятность очень редкого события появления сильного пика среди отсчетов axy (j ), j = 1, 2, ..., N/2 - 1. При этом имеют место две ситуации. 1. Частота, на которой два ряда имеют общий гармонический компонент, заранее не известна. В этом случае принимается, что число a = Sup{axy (j )} соответствует искомому общему периодическому компоненту. Если выполняется неравенство

a

22 x y Z1 , N2

(56)

где величина Z1 является корнем уравнения

1 - 1 - 2 Z1 K1 (2

Z1 )

N/2-1

= q,

(57)

то с вероятностью p = 1 - q принимается утверждение, что отсчет amax порождается общим гармоническим компонентом обоих рядов, а не генерируется шумами. 2. Если частота искомого общего сигнала заранее известна и среди множества значений {axy (j )} ей соответствует отсчет al , то при выполнении неравенства

al

22 x x Z2 , N2

(58)

где величина Z2 является корнем уравнения

2 Z2 K1 (2 Z2 ) = q ,

(59)

с вероятностью p = 1 - q полагают, что отсчет al тоже порождается общими сигналами, а не шумами. 18


Отметим, что порог обнаружения сигнала Z2 не зависит от длины ряда. Для наиболее часто применяющихся уровней значимости q = 0.001; q = 0.01 и q = 0.05 соответствующие значения порога Z2 равны 15.5; 8.2 и 4.0. В противоположность этому значения порога Z1 зависят как от принятого уровня значимости, так и от длины ряда. В табл. 2 приведены значения порога Z1 как функции числа точек ряда при q = 0.05 и q = 0.01. Подобно обычной периодограмме взаимную периодограмму можно сгладить. Для этого нужно воспользоваться формулой (40), подставив в нее взвешенные значения взаимной коррелограммы

cxy [m] = cxy [m] Wm , ~ ?
где значения коэффициентов W
m

m = 0, 1, . . . , N - 1,

(60)

определяются формулой (29).

Вычислительные процедуры
Дискретное преобразование Фурье
Для получения численных значений спектральных и корреляционных характеристик равномерных временных рядов пользуются Таблица 2: Значения Z1 в зависимости от N и q .

N
50 100 150 200 250 300 350 400 450 500

q = 0.05
13.8 16.7 18.6 19.9 21.0 21.8 22.6 23.3 23.9 24.4

q = 0.01
20.8 24.4 26.2 28.1 29.2 30.5 31.3 32.0 32.6 33.2

N
600 700 800 900 1000 1200 1400 1600 1800 2000

q = 0.05
25.4 26.2 26.9 27.5 28.2 29.2 30.1 30.8 31.5 32.1

q = 0.01
34.6 35.5 36.3 36.9 37.5 39.0 40.0 40.8 41.6 42.2

19


дискретными преобразованиями Фурье (ДФТ), которые мы введем следующим образом. Пусть имеются две последовательности комплексных величин
{xk }, k = 0, 1, 2, . . . , N - 1, {Xj }, j = 0, 1, 2, . . . , N - 1.
Назовем прямым дискретным преобразованием Фурье (DFT) следующее выражение:
N -1

Xj = DF Tj {xk }

N -1 k=0

=
k=0

xk e

-i

2 N

kj

,

(61)

j = 0, 1, 2, . . . , N - 1.
При этом обратное дискретное преобразование (DF T вводится посредством соотношения
-1

) Фурье

xk = DF T

-1 k

{Xj }N -1 = j =0

1 N

N -1

Xj e
j =0

i

2 N

kj

,

(62)

k = 0, 1, 2, . . . , N - 1.
Мы не будем останавливаться здесь на свойствах дискретных преобразований Фурье, поскольку они хорошо описаны в многочисленных руководствах (Дженкинс Г., Ваттс Д., 1972; Отнес и Эноксон, 1982 и др.). Отметим только то, что с помощью дискретного преобразования Фурье вычисление периодограмм на фундаментальной системе частот (10) можно производить по следующим формулам

Dj = Dxy [j ] =

1 DF Tj {xk } N2

N -1 2 k=0

,
N -1 k=0

(63)

1 - DF Tj {xk }N=01 DF Tj {yk } k N2 j = 0, 1, 2, . . . , N/2.

,

(64)

20


Быстрое преобразование Фурье
Для вычисления преобразований (61) и (62) требуется выполнить N 2 основных операций. Это обстоятельство делает вычисление дискретного преобразования Фурье очень трудоемкой операцией, и поэтому в прошлом старались строить алгоритмы анализа данных таким образом, чтобы по возможности избежать вычисления ДПФ. Положение существенно изменилось после изобретения специфических вычислительных процедур (быстрых преобразований Фурье, БПФ или F F T ), которые позволили увеличить скорость вычислений в десятки и сотни раз. Один из самых распространенных алгоритмов быстрого преобразования Фурье (Кули и Тьюки, 1965) предполагает, что число точек ряда N представимо следующим образом:

N = 2p ,

(65)

где p целое положительное число. Этот алгоритм требует для своей реализации всего 2N log2 N операций, что при больших N существенно меньше, чем N 2 . Подчеркнем, однако, что быстрое преобразование Фурье не является каким-либо новым типом преобразования, а представляет собой лишь один из вариантов вычисления дискретных преобразований Фурье (61) и (62). Ниже приведена процедура быстрого вычисления дискретного преобразования Фурье для рядов, число точек которых представимо формулой N = 2p , где p целое положительное число (Ефимов и др., 1980). Текст программы соответствует p = 9 и N = 512. При желании значения этих параметров можно изменить. Для этого нужно внести вполне понятные исправления в операторах программы, отмеченных знаком {*}.

21


unit t; interface {Type arr = array[1..512];} {*}. { Процедура Быстрого Преобразования Фурье a,b - массивы преобразуемой функции; aa,bb - массивы преобразованной функции; kind - вид преобразования: kind = 0 - прямое преобразование; kind = 1 - обратное преобразование; } implementation const n=9; n1=512; {*} var aaa,bbb: arr; var v,m,k,j,i,ni,jm:longint; var c,si,co,ass,bob:double; { Вычисляет 2m } var y,k:longint; begin y:=1; if m=0 then y:=1 else for k:=1 to m do y:=2*y; f:=y end; {f}

procedure t (kind:longint;var a,b,aa,bb:arr);

procedure t(kind:longint;var a,b,aa,bb:arr);

function f(m:longint):longint;

22


begin {Сохранение входных массивов} for i:=1 to 512 do {*} begin aaa[i]:=a[i]; bbb[i]:=b[i] end; m:=1; repeat k:=0; repeat v:=0; repeat j:=f(m)*k+v+1; i:=f(m-1)*k+v+1; c:=3.141593*v/f(m-1); co:=cos(c); if kind=0 then si:=sin(c) else si:=-sin(c); ni:=f(n-1)+i; jm:=f(m-1)+j; ass:=a[ni]; bob:=b[ni]; aa[j]:=a[i]+ass*co+bob*si; bb[j]:=b[i]-ass*si+bob*co; aa[jm]:=a[i]-ass*co-bob*si; bb[jm]:=b[i]+ass*si-bob*co; v:=v+1; until v-f(m-1)0; k:=k+1; until k-f(n-m)0; m:=m+1; for i:=1 to n1 do begin a[i]:=aa[i]; b[i]:=bb[i]; end; until m-n>0; if kind>0 then for i:=1 to n1 do begin aa[i]:=aa[i]/n1; bb[i]:=bb[i]/n1; end; {Восстановление входных массивов} for i:=1 to 512 do {*} begin a[i]:=aaa[i]; b[i]:=bbb[i] end; end; end.

23


Дополнение ряда нулями
Поскольку в общем случае трудно ожидать выполнения равенства (65), обычно для применения быстрого преобразования Фурье исходный ряд дополняют нулями таким образом, чтобы длина нового ряда стала равной N1 = 2p N . Часто длину нового ряда делают даже равной N2 = 2N1 . Добавление нулей является единственным методическим новшеством, вытекающем из специфики алгоритма быстрого преобразования Фурье. Другими словами, после добавления нулей наши исходные данные будут иметь следующую структуру:

xk =

x(tk ), k = 0, 1, 2, . . . , N - 1, 0, k = N , . . . , N1 -1, . . . , N2 - 1.

(66)

Вычислим значения периодограммы

1 | F F Tj {xk } N2 для дискретного набора частот Dj =
j =

N2 -1 2 k=0 |

(67)

2c j, j = 0, 1, . . . , N2 /2, N2 следующих друг за другом с шагом = 2c /N2 .

(68)

(69)

Если в формуле (67) использовать дискретное преобразование Фурье и вычислять его непосредственно (без добавления нулей), то отсчеты периодограммы будут отнесены к фундаментальной системе частот 2c N j = j, j = 0, 1, . . . , , (70) N 2 отстоящих друг от друга на шаг

=
Из сравнения (69) и (71) находим

2c . N

(71)

N = < 1. N2
24

(72)

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

Вычисление коррелограмм с помощью БПФ
Отметим также, что для вычисления коррелограмм можно использовать не только соотношения (14) и (38), но и процедуру обратного дискретного преобразования Фурье. Рабочие формулы выглядят следующим образом:
- cm = N Re F F Tm 1 {Dj } ? -1 m -1 m N2 -1 j =0

, , ,

(73) (74) (75)

cxy [m] = N Re F F T ? cyx [m] = N Re F F T ?

{Dxy [j ]}
{Dxy [j ]}

N2 - 1 j =0 N2 - 1 j =0

m = 0, 1, . . . , N - 1.
При выполнении этих операций следует иметь в виду, что с целью исключения так называемого кругового эффекта ("The wraparound problem", Отнес и Эноксон, 1982) дополнение ряда нулями до значения N2 должно быть сделано обязательно, даже если соответствующие преобразования Фурье выполняются без использования алгоритма быстрого преобразования Фурье.

25


Алгоритмы
Спектрально-Корреляционный Анализ Временных Рядов ( СКАВРя)
Исходные данные: 1) Равномерный временной ряд

xk = x(tk ),

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

составленный из действительных чисел. Для использования процедур быстрого преобразования Фурье элементам массива мнимых частей наших данных следует присвоить нулевые значения; 2) X1 или X2 критические значения для разделения шумового и детерминированного компонентов временного ряда, отвечающие уровню значимости q (см. с. 1013); 3) a, N параметры корреляционного окна Тьюки (см. с. 1214). Ниже мы опишем основные шаги алгоритма и проиллюстрируем их исполнение для модельного временного ряда

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

(76)

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

где отношение "сигнал к шуму". Заданные параметры: t = 1 c; N = 230; q = 0.01; X1 = 9.0; A1 = 1; 1 = 0.1 Гц; 1 = 0; = 0.50; = 0.1; = 0.05.
Алгоритм спектрального анализа состоит из следующих шагов. 26


1. Графическое представление исходного ряда во вре-

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

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

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

1 m= N

N -1

xk ,
k=0

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

После исключения линейного тренда график центрированного временного ряда (76) представлен на рис. 5.

27


Рис. 4: График модельного ряда (76).

Рис. 5: График центрированного модельного ряда (76).

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

28


4. Вычисление периодограммы Чтобы воспользоваться процедурой БПФ, исходный ряд нужно дополнить нулями таким образом, чтобы длина нового ряда была N1 = 2p N . Для получения коррелограммы длину нового ряда делают равной N2 = 2N1 . Для такого ряда с помощью процедуры быстрого преобразования Фурье (c. 21) получают
N2 - 1

Xj = F F Tj {x } k

N2 - 1 k=0

=
k=0

x e k

-i

2 N2

kj

,

j = 0, 1 . . . , N2 - 1,
(78)

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

Dj =

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

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

Отсчеты периодограммы соответствуют частотам

j = j,
где

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

=

При обработке модельного ряда (76) были выбраны значения N1 = 512 и N2 = 1024.
5. Оценивание дисперсии ряда:
2 0 =

1 N -1

N -1

(x )2 . k
k=0

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

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

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

Для модельного ряда эти результаты показаны на рис. 6.
29


Если не ставится задача получения коррелограммы и сглаженной периодограммы, то на этом выполнение алгоритма можно прекратить. 7. Вычисление коррелограммы методом БПФ:

cm = ?

1 2 - Re F F Tm 1 [| Xj | ] N

N2 -1 j =0

,

(79)

m = 0, 1, 2, . . . , N - 1.
8. Графическое представление коррелограммы

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

W

m

= (1 - 2a) + 2a cos

m . N

(80)

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

cm = cm Wm , ~ ?

m = 0, 1, . . . , N - 1.

(81)

Для модельного ряда (76) были выбраны значения N = 0.1N и a = 0.25.
10. Вычисление сглаженной периодограммы Массив значений коррелограммы c? дополняют нулями и с m новым массивом длиной N2 осуществляют вычисление сглаженной периодограммы по формуле

(82) ~j относятся к той же системе чаОтсчеты периодограммы D стот, на которой была определена периодограмма Dj .

1 ~ Dj = [2 Re F F Tj {cm } ~ N

N2 -1 m=0

- c0 ], ~

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

30


Рис. 7: График смещенной коррелограммы модельного ряда (76).

Рис. 8: Сглаженная периодограмма ряда (76). Параметры сглаживания: N = 0.5N ; a = 0.25.

Рис. 9: Сглаженная периодограмма ряда (76). Параметры сглаживания: N = 0.1N ; a = 0.25.

31


11. Представление сглаженной периодограммы в графи-

ческом виде

Для модельного ряда (76) графики сглаженных периодограмм показаны на рис. 8 и 9. Сравнение рис. 6 с рис. 8 и 9 показывает, что по мере уменьшения значения параметра N в сглаженной периодограмме представление шумового компонента становится все более гладким, в то время как контур спектральной линии на частоте 0.1 Гц становится все менее четким.
12. Конец алгоритма.

Взаимный Анализ Временных Рядов (ВАВРя)
Ниже приводится алгоритм взаимного спектрально-корреляционного анализа двух рядов, который является обобщением соответствующего алгоритма анализа одного временного ряда. Исходные данные:

ћ t временной шаг задания рядов; ћ равномерные временные ряды xk = x(tk ), yk = y (tk ), tk = t k , k = 0, 1, ..., N - 1, tk = t k , k = 0, 1, ..., N - 1.

В дальнейшем мы будем считать, что наши ряды составлены из действительных чисел. Для использования процедур быстрого преобразования Фурье элементам массивов мнимых частей наших данных следует присвоить нулевые значения;

ћ Z1 или Z2 критические значения для выделения полезного сигнала из шума, отвечающие уровню значимости q (с. 18); ћ a, N параметры корреляционного окна Тьюки (с. 1214).

Ниже мы опишем основные шаги алгоритма и проиллюстрируем их исполнение для модельных временных рядов
xk = A1 cos(2 1 tk + 1 ) + x k ,
(83)

32


2

yk =
i=1

Bi cos(2 i tk + i ) + y k ,

(84)

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

где

t постоянный шаг выборки; A1 , 1 , 1 , амплитуда, частота и фаза ряда xk ; Bi , i , i , i = 1, 2, амплитуды, частоты и фазы ряда yk ; k , k , k = 0, 1, 2, . . . , N - 1, случайные величины с нулевым математическим ожиданием и единичной дисперсией, распределенные по нормальному закону (шумовые компоненты рядов xk и yk ); x , y , среднеквадратические отклонения шумовых компонентов рядов xk и yk , задаваемые с помощью соотношений x = A2 1 , 2x y =
2 2 B1 + B2 , 2y

(85)

где x , y отношения "сигнал к шуму"соответственно для первого и второго рядов. Заданные параметры: t = 1 c; N = 200; q = 0.01; Z1 = 28.1; A1 = 1; 1 = 0.11 Гц; 1 = 0; B1 = 1; 1 = 0.11 Гц; 1 = /3; B2 = 1; 2 = 0.18 Гц; 2 = 0; x = 0.1; y = 0.1;
Алгоритм спектрального анализа состоит из следующих шагов. 1. Графическое представление исходных рядов во вре-

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

Графики модельных рядов (83) и (84) показаны на рис. 10.
2. Центрирование рядов:

x = xk - mx , k
yk = xk - my ,

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

33


Рис. 10: Модельные ряды: сплошная линия ряд (83), штриховая линия
ряд (84).

где

mx =

1 N

N -1

xk ;
k=0

my =

1 N

N -1

yk .
k=0

3. Графическое представление центрированных рядов

Поскольку модельные ряды (83) и (84) не содержат постоянных составляющих, графики центрированных рядов повторяют графики на рис. 10.
4. Дополнение рядов нулями Для того чтобы воспользоваться процедурой БПФ, исходные ряды дополняют нулями. Выберем числа N1 = 2p N и N2 = 2N1 . Дополним каждый из исходных рядов нулями с тем, чтобы в процессе дальнейшей обработки оба ряда соcтояли бы из N2 точек.

При обработке модельных рядов (83) и (84) были выбраны значения N1 = 512 и N2 = 1024.
5. Вычисление взаимной периодограммы Сначала вычисляются дискретные преобразования Фурье по формулам
N2 - 1

Xj = F F Tj {x } k

N2 - 1 k=0

=
k=0

x e k

-i

2 N2

kj

,

j = 0, 1 . . . , N2 - 1,
(86)

34


N2 -1 Yj = F F Tj {yk } N2 -1 k=0

=
k=0

yk e

-i

2 N2

kj

,

j = 0, 1 . . . , N2 - 1,
(87)

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

Dxy [j ] =

1 X Yj , N2 j

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

Отсчеты взаимной периодограммы соответствуют частотам

j = j,
где

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

=

6. Вычисление коспектра и квадратурного спектра:

Pxy [j ] = Re Dxy [j ];

Qxy [j ] = -Im Dxy [j ].

7. Вычисление квадрата амплитуды взаимной периодо-

граммы:

2 A2 [j ] = Pxy [j ] + Q2 [j ]. xy xy

8. Вычисление фазового спектра:

xy [j ] = -arctg
9. Оценивание дисперсий рядов:
2 x =

Qxy [j ] . Pxy [j ]

1 N -1

N -1

(x )2 , k
k=0

2 y =

1 N -1

N -1 (yk )2 . k=0

10. Графическое представление квадрата амплитуды вза-

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


Рис. 11: Квадрат амплитуды взаимной периодограммы рядов (83) и
(84). Штриховая линия 99-процентный порог разделения шумовой и детерминированной составляющих обоих рядов.

22 уровень x y Z1 /N 2 , считаются значимыми, т. е. принадлежащими общему детерминированному компоненту рядов. Вероятность такого утверждения равна 1 - q .

Для модельных рядов эти результаты показаны на рис. 11. Мы видим, что в число значимых пиков с вероятностью 0.99 попадает только один пик на частоте 0.11 Гц. Остальные пики порождаются шумовыми компонентами обоих рядов.
11. Графическое представление фазового спектра Можно напечатать значения фазовогo спектра для всех значений j , однако целесообразно получить фазовый спектр только для тех частот, на которых обнаружены общие значимые гармоники. Другими словами, на печать выводятся только те значения фазового спектра, для которых выполняются соотношения 22 x y A2 [j ] Z1 . xy N2

Для модельных рядов фазовый спектр показан на рис. 12, где представлен только один пик на частоте 0.11 Гц. Высота этого пика, равная 0.89, хорошо соответствует заданной разности фаз /3. Заметим, что при увеличении значения N2 эту оценку можно улучшить за счет увеличения формального разрешения.
36


Рис. 12: Фазовый спектр рядов (83) и (84). Если не ставится задача получения взаимной коррелограммы и сглаженной взаимной периодограммы, то на этом выполнение алгоритма можно прекратить. 12. Вычисление взаимных коррелограмм методом БПФ:

1 - Re F F Tm 1 {Xj Yj }N2 -1 , j =0 N 1 - cyx [m] = ? Re F F Tm 1 {Xj Yj }N2 -1 , j =0 N cxy [m] = ? m = 0, 1, 2, . . . , N - 1.
13. Графическое представление коррелограмм

(88) (89)

Графики взаимной коррелограммы cxy [m] и cyx [m] для первых ? ? двадцати точек показаны на рис. 13.
14. Взвешенная коррелограмма Для получения сглаженной периодограммы коррелограммы умножают на весовую функцию Тьюки

W

m

= (1 - 2a) + 2a cos

m . N

(90)

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

cxy [m] = cxy [m]Wm , ~ ?
37

m = 0, 1, . . . , N - 1,

(91)


Рис. 13: Взаимные корреллограммы рядов (83) и (84). Сплошная линия
cxy [m], штриховая cyx [m]. ? ?

Рис. 14: Квадрат амплитуды сглаженной взаимной периодограммы рядов (83) и (84). Параметры сглаживания: N = 0.1N ; a = 0.25.

Рис. 15: Сглаженный фазовый спектр рядов (83) и (84).

38


cyx [m] = cyx [m]Wm , ~ ?

m = 0, 1, . . . , N - 1.

(92)

Для взвешивания взаимных коррелограмм модельных рядов (83) и (84) были выбраны значения N = 0.1N и a = 0.25.
15. Вычисление сглаженной периодограммы Массив значений коррелограмм cxy [m] и cyx [m] дополняют ну~ ~ лями и с новыми массивами длиной N2 осуществляют вычисление сглаженной взаимной периодограммы по формуле

1 ~ ~ Dxy [j ] = F F Tj {cxy [m]}N2 -1 + F F Tj {cyx [m]} ~ m=0 N

N2 - 1 m=0

- c0 . ?

16. Вычисление сглаженных коспектра и квадратурного

спектра:

~ ~ Pxy [j ] = Re Dxy [j ]; ~ ~ Qxy [j ] = -Im Dxy [j ].

17. Вычисление квадрата амплитуды взаимной коррело-

граммы:

~ ~2 ~ A2 [j ] = Pxy [j ] + Q2 [j ]. xy xy

18. Вычисление фазового спектра:

~ Qxy [j ] ~ xy [j ] = -arctg . ~ Pxy [j ]
Отсчеты сглаженных характеристик взаимной периодограм~ мы Dxy [j ] относятся к той же системе частот, на которой была определена периодограмма Dxy [j ]. 19. Графическое представление сглаженных оценок

Для модельных рядов (83) и (84) график квадрата амплитуды сглаженной взаимной периодограммы показан на рис. 14, а фазовый спектр сглаженной взаимной периодограммы на рис. 15. Мы видим, что, несмотря на сильную зашумленность обоих рядов, в их составе очень уверенно обнаруживается общая гармоника с частотой 0.1 Гц.
20. Конец алгоритма.

39


Упражнения
Составьте программу анализа временных рядов (СКАВРя), пользуясь алгоритмом, приведенным на с. 2632. Предусмотрите в программе возможность анализа модельного временного ряда

xk = A1 cos(2 1 tk - 1 ) + A2 cos(2 2 tk - 2 ) + x k , tk = t k , k = 0, 1, 2, . . . , N - 1,

(93) (94)

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

x =

A2 + A2 1 2 . 2

(95)

Задаваемые параметры ряда: 1) характеристики гармоник Ai , i , i , i = 1, 2; 2) отношение "сигнал к шуму" . Используя эту программу, выполните следующие упражнения.

ћ Оцените аналитически длину ряда N t, необходимую для полного разрешения в периодограмме двух близких спектральных линий, соответствующих гармоническим компонентам с частотами 1 и 2 . Используя программу СКАВРя, постройте периодограммы рядов с меньшей и большей длинами для того, чтобы увидеть эффекты слияния линий, а также картины их частичного и полного разрешения. ћ Почему на рис. 18 контуры линий видны, а на рис. 19 нет? ћ Для заданных амплитуд ряда (93) оцените аналитически минимальное значение отношения "сигнал к шуму", при котором с заданной вероятностью еще удается выделить сигнал из
40


шума. С помощью численных экспериментов проверьте Ваш результат. Почему иногда результаты программы СКАВРя не совпадают с теоретическими?

ћ Какие изменения нужно ввести в программу СКАВРя, чтобы по оцифровке оси ординат периодограммы можно было бы сразу определять численные значения амплитуды гармонических колебаний? ћ Добавьте к ряду (76) постоянное слагаемое C > Ai и постройте с помощью программы СКАВРя периодограммы нецентрированного и центрированного рядов. Объясните различный вид периодограмм, написав аналитическое выражение для периодограммы ряда xk = C, k = 0, 1, ..., N - 1, где C = const. ћ Дан временной ряд xk = x(t k ), k = 0, 1, ..., 2N - 1, у которого первое и второе, третье и четвертое и т. д. значения попарно совпадают. Докажите аналитически, что периодограмма такого ряда будет симметрична относительно частоты = 0.5c = 1/4t. Проиллюстрируйте это свойство с помощью программы СКАВРя. ћ Как определить амплитуды и фазы статистически значимых гармонических компонентов ряда, пользуясь только теми операциями, которые выполняются в программе СКАВРя ? ћ С помощью программы СКАВРя продемонстрируйте эффект неразличимости частот ("элайзинг"). ћ Вероятность того, что наибольший отсчет периодограммы дискретного белого шума превысит заданный уровень X1 , определяется уравнением (24). Как с помощью программы СКАВРя можно проиллюстрировать этот результат? ћ Спектр мощности и коррелограмма суть неотрицательные функции по определению. Почему в методе БлэкманаТьюки (особенно хорошо это заметно при a = 0) значения периодограммы на некоторых частотах бывают отрицательными? ћ Для гармонического компонента с амплитудой A и частотой 0 точное значение максимального отсчета периодограммы
41


равно A2 /4. При использовании программы СКАВРя Вы можете заметить, что это соотношение выполняется не всегда. Как это объяснить?

ћ Каким образом по коррелограмме временного ряда (76) можно оценить параметр (отношение "сигнал к шуму")? ћ Докажите соотношение (13). ћ Докажите соотношение (40). ћ Докажите соотношение (48). ћ Докажите соотношение (50). ћ Чем различаются фазовые спектры взаимных периодограмм Dxy ( ) и Dyx ( )?

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

LOD(t) = 86400 + LOD(t),
где LOD(t) (Length Of Day) период вращения Земли, измеряемый в секундах SI, и координаты полюса x(t) и y (t), которые задают положение мгновенного полюса вращения Земли в системе отсчета, жестко связанной с Землей (ось X направлена по касательной к Гринвичскому меридиану, ось Y по касательной к меридиану 90 западной долготы). В настоящее время эти три ряда а также еще два ряда, задающих движение полюса в небесной системе отсчета (так называемые Polar Osets), определяются различными методами (РСДБ, лазерная локация ИСЗ и Луны, GPS) и обрабатываются в Международной Службе Вращения Земли (IERS). Публикации этих рядов производятся в годичных отчетах IERS. Полные ряды можно найти в Интернете на сайте IERS: http://hpiers.obspm.fr.

42


Рис. 16: График ряда вариаций периода вращения Земли (LOD), дополненного нулями.

Вариации продолжительности суток
Ниже приводятся результаты обработки ряда вариаций периода вращения Земли (LOD) на промежутке 19621996 гг. по данным Международной службы вращения Земли (IERS, файл eopc04). Шаг выборки t = 1 сут., число точек ряда N = 12462. Исходные значения представлены на рис. 1. Дополненный нулями ряд вариаций LOD (N2 = 32 Ч 1024) показан на рис. 16. Результаты вычисления периодограммы Шустера в различных диапазонах частот демонстрируются на рис. 1719. В низкочастотной части спектра имеется мощный компонент с периодом 22.7 года и амплитудой 0.6 мс. Ее происхождение связано с вековыми колебаниями LOD. На рис. 18 видны два пика на частотах 1 и 2 цикла в год с амплитудами 0.35 мс. Эти пики соответствуют изменениям периода вращения Земли под действием сезонных атмосферных перемещений. На рис. 19 еще две спектральные линии имеют частоты 13.37 и 26.74 цикла в год и амплитуды 0.18 и 0.37 мс. Периоды этих колебаний, происходящих под действием лунных приливов, равны соответственно 27.32 и 13.66 сут. На рис. 1719 порог X1 выделения сигнала из шума, соответствующий уровню значимости q = 0.01, показан штриховой линией. Мы видим, что максимальные значения имеющихся в периодограмме пиков намного превосходят это значение.

43


Рис. 17: Низкочастотный компонент в спектре LOD.

Рис. 18: Сезонные компоненты в спектре LOD. Штриховая линия 99процентный порог разделения шумовой и детерминированной составляющих ряда.

Рис. 19: Приливные компоненты в спектре LOD. Штриховая линия
99-процентный порог разделения шумовой и детерминированной составляющих ряда.

44


Рис. 20: X-координата полюса (1962-1996).

Рис. 21: Y-координата полюса (1962-1996).

Движение полюса в теле Земли
Теперь приведем результаты обработки рядов x- и y -координат полюса на промежутке времени 19621996 гг. по данным Международной службы вращения Земли (IERS, массив eopc04). Шаг выборки t = 1 сут., число точек ряда N = 12462. В нашем анализе ряды координат полюса дополнялись нулями (N2 = 32 Ч 1024). Исходные ряды представлены на рис. 20 и 21 сплошными линиями. На этих же графиках штриховыми линиями показаны линейные тренды, являющиеся следствием векового движения полюса. По нашей выборке вековое движение происходит со скоростью 0 .44 в столетие в направлении 65 к западу от Гринвича. Заметим, что по более длительным выборкам эти оценки получаются иными. Например, по рядам координат полюса на промежутке 1899.71992.0 гг., редуцированным к системе каталога HIPPARCOS, оценка ско45


Рис. 22: Квадрат модуля взаимной периодограммы координат полюса.

Для уровня значимости q = 0.01 и N = 12462 решение уравнения (57) дает Z1 = 55.5, поэтому порог обнаружения значимых пиков составляет величину 1 .1 Ч 10 -10 и в масштабе графика не виден.

Рис. 23: Фазовый спектр координат полюса. рости векового движения полюса составляет 0 .34 в столетие в направлении 75.8 к западу от Гринвича (Вондрак Я., Рон К., 1995). Квадрат амплитуды взаимной периодограммы координат полюса показан на рис. 22. Здесь мы видим два пика на частотах Чандлеровской (0.85 год-1 ) и годичной составляющей (1 год-1 ). Фазовый спектр координат полюса (рис. 23) показывает значительную разность фаз обоих компонентов.

46


Литература
Блэкман и Тьюки, 1959. Blackman R.B., Tukey J.W. The Measurement of Power Spectra. N.Y.: Dower. Вондрак Я., Рон К., 1995. Vondrak J., Ron C. Earth Rotation, Reference Systems in Geodynamics and Solar System // Proc. JOURNEES 1995, Warsaw. Дженкинс Г., Ваттс Д., 1972. Спектральный анализ и его приложения. Т. 1, 2. М.: Мир. Ефимов А.В., Золотарев Ю.Г., Терпигорева В.М., 1980. Математический анализ (специальные разделы). Т.1, 2. М.: Высшая школа. Кули и Тьюки, 1965. Cooley J.W., Tukey J.W. // Math. Comput. Vol. 19. P. 297-301. Отнес и Эноксон, 1982. Прикладной анализ временных рядов. М.: Мир. Теребиж В.Ю., 1992. Анализ временных рядов в астрофизике. М.: Наука.

47

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

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

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

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