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

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

В.В. Витязев

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

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

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

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

В54

Витязев В.В.

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

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

.


.

2


.

Содержание
Введение . . . . . . . . . . . . . . . . . . . . Периодограмма и спектр мощности . . . . . Спектральные окна . . . . . . . . . . . . . . Равномерная дискретность . . . . . . . Случайные пропуски . . . . . . . . . . Периодические пропуски . . . . . . . . Длинный пропуск . . . . . . . . . . . . LS-спектры . . . . . . . . . . . . . . . . . . . LS-спектр Барнинга . . . . . . . . . . . LS-спектр Ломба . . . . . . . . . . . . . Модификация Скаргла . . . . . . . . . LS-спектры и периодограмма Шустера Свойства периодограмм белого шума . . . . Периодограмма Ломба . . . . . . . . . Периодограмма Шустера . . . . . . . . Выделение сигнала из шума . . . . . . . . . Восстановление истинного спектра . . . . . Метод CLEAN и его модификации . . . . . Описание метода . . . . . . . . . . . . . Модифицированный алгоритм CLEAN Численные примеры . . . . . . . . . . Упражнения . . . . . . . . . . . . . . . . . . Литература . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 5 11 11 12 13 17 20 21 22 23 23 31 32 33 34 39 44 44 50 54 56 57

3


.

Введение
В настоящее время теория спектрально-корреляционного анализа временных рядов с исчерпывающей полнотой развита для равномерных рядов (Андерсон Т., 1976; Бендат Дж., Пирсол А., 1989; Блэкман и Тьюки, 1959; Отнес Р., Эноксон Л., 1982; Теребиж В.Ю., 1992 и др.). Однако получить измерения с постоянным шагом по времени удается весьма редко. Такая ситуация особенно характерна для наземной оптической астрономии, где ночные наблюдения прерываются светлой частью суток и из-за погодных условий могут возникать периодические или длинные пропуски наблюдений. Нерегулярность свойственна не только наблюдениям, выполняемым с поверхности Земли, она имеет место и для измерений, производящихся с борта ИСЗ: вхождение аппарата в радиационные пояса, экранирование наблюдаемого участка небесной сферы телом Земли и т. д. Ко всем таким явлениям добавляются чисто технические причины: отказы или сбои регистрирующей аппаратуры, длительные капитальные ремонты или реконструкции инструментов. Как было сказано выше, теория и эффективные методы получения спектральных оценок хорошо разработаны для равномерных рядов. Чтобы применить эти методы к неравномерным рядам, используются различные способы их редукции к регулярной сетке временных отсчетов: осреднение данных на интервалах постоянной длины, аппроксимация наблюдений различными математическими выражениями и т. д. Хотя употребляющиеся методы приведения наблюдений к постоянному шагу и могут быть достаточно хорошими при исследовании частных задач, они не решают проблему в общем случае. Настоящее учебное пособие предназначено для изучения теории и практических методов, которые позволили бы с достаточной полнотой проводить спектральный анализ временных рядов с произвольным распределением моментов наблюдений.

4

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

Теорема 1. Пусть временной ряд X (t) представляет собой цен-

трированный случайный процесс, заданный дискретным множеством реализаций
X (t) = {xp (t)}
M p=1

,

0 t T.

(1)

Будем считать, что мы производим измерения лишь в произвольных упорядоченных временных точках tk , k = 0, 1, . . . , N - 1. В этом случае между средней периодограммой (в дальнейшем угловые скобки ... символ математического ожидания)
D( ) = 1 N2
N -1 2

xp (tk )e
k=0

-i t

k

(2)

и спектральным окном
1 W ( ) = 2 N
N -1 2

e
k=0

-i t

k

(3)

существует соотношение
+

D ( ) =
-

g ( )W ( - ) d ,

(4)

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


Теорема 2. Пусть временной ряд представляет собой полигармоническую функцию
n

x(t) =
p=1

Ap cos(p t + p ),

0 t T < ,

(5)

где Ap , p , p соответственно амплитуда, частота и фаза p-й гармоники. Введем в рассмотрение корреляционную функцию и спектр мощности полигармонической функции (5):
k ( ) = lim 1 2
+

T

1 T

T

n

x(t)x(t + ) dt =
0 p=1

A2 p cos k , 2

(6)

g ( ) =

k ( )e
-

-i

d =

1 4

n

A2 [ ( + p ) + ( - p )], (7) p
p=1

где ( ) дельта-функция. В этом случае спектр мощности (7), периодограмма
1 D ( ) = 2 N
N -1 2

xk e
k=0

-i t

k

,

(8)

и спектральное окно
W ( ) = ( ) ( ),
(9)

где
( ) = 1 N

N -1

e
k=0

-i t

k

,

(10)

связаны между собой соотношением
+

D( ) =
-

g ( )W ( - ) d + D0 ( ),

(11)

6


в котором член D0 ( ) определяется формулой
4D 0 ( ) =
n p=1

A2 [( - p ) ( + p )e p + ( - p )( + p )e

i 2

p

+ ]+ +
(12)

-i2

p

+

n
p,q =1 p =q

Ap Aq [( - p ) ( + q )ei + ( - p )( + q )e

(p +q )

-i(p +q )

]+ +

+

n
p,q =1 p =q

Ap Aq [( - p ) ( - q )ei + ( + p ) ( + q )e

(p -q )

-i(p -q )

].

Сделаем теперь несколько замечаний к теоремам 1 и 2. 1. Множество функций (5) со случайными фазами p , распределенными равномерно внутри промежутка [0, 2 ], представляет собой стационарный случайный процесс. Осредняя (11) по множеству реализаций, мы вновь получим соотношение (4). 2. Сравнение (4) с выражением (11) показывает, что при переходе от множества реализаций к одной реализации структура периодограммы, вообще говоря, изменяется за счет появления слагаемого D0 ( ). 3. Влияние слагаемого D0 ( ) на периодограмму может быть существенным лишь только в некоторых специфических ситуациях, на которых мы остановимся позже. Пока же мы будем заниматься такими рядами, для которых членом D0 ( ) в выражении (11 ) можно пренебречь. В этом случае соотношения (4) и (11) позволяют достаточно просто объяснить особенности периодограмм временных рядов с пропусками. Действительно, считая, что в спектре мощности g ( ) 7


имеются концентрации мощности в бесконечно узких интервалах частот g ( ) = [ ( + 0 ) + ( - o )], (13) из уравнения (4) получаем

D( ) = [W ( + 0 ) + W ( - 0 )].

(14)

Теперь мы видим, каким образом истинный спектр мощности превращается в периодограмму. При этом важно подчеркнуть, что переход от (13) к (14) определяется только видом спектрального окна, и поэтому соотношения (13) и (14) справедливы для любого распределения отсчетов времени, в том числе и равномерного. В большинстве случаев нерегулярного распределения данных в спектральных окнах W ( ) можно обнаружить центральный пик на частоте = 0 и боковые пики на некоторых частотах +k , k = 1, 2 . . . ? Кроме этого, существует еще и структура боковых лепестков. Последние возникают по причине конечности интервала наблюдений, в то время как боковые пики определяются неравномерностью распределения данных. Как видно из (14), все эти детали будут присутствовать в периодограмме D( ). Действительно, центральный пик окна будет виден на частоте = 0 . Поскольку он соответствует реальной концентрации мощности, его следует назвать истинным пиком. В свою очередь боковые пики появятся на частотах | 0 + k |, где истинный спектр мощности не имеет никаких линий. По этой причине пики на частотах | 0 + k | называют духами, а периодограмму D( ), засоренную эффектами нерегулярности наблюдений, часто называется грязным спектром. На рис. 1 и 2 показаны результаты спектрального анализа двух рядов, составленных из двух косинусоид с периодами 1.3 и 0.7 года, на сетке из 120 равномерно распределенных точек (рис. 1) и для тех же точек с периодическими годичными пропусками (рис. 2). Эти рисунки демонстрируют специфическую особенность спектрального анализа неравномерных рядов появление ложных линий в периодограмме изучаемого процесса. Причины появления фиктивных пиков и способы их устранения из периодограмм рассматриваются в следующих параграфах.

8


Рис. 1: Спектральный анализ суммы двух косинусоид с равными ампли-

тудами и периодами P = 1.3 и P = 0.7 года; a регулярные наблюдения: 120 точек с шагом 0.1 года; b спектральное окно; c периодограмма ("чистый"спектр).

9


плитудами и периодами P1 = 1.3 и P2 = 0.7 года; a нерегулярные наблюдения с годичными пропусками; b спектральное окно; c периодограмма ("грязный"спектр).

Рис. 2: Спектральный анализ суммы двух косинусоид с равными ам-

10

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

w(t) =

1 N

N -1

(t - tk ).
k=0

(15)

В этом случае для спектрального окна имеем
+ 2

W ( ) =
-

w(t)e

-i t

dt .

(16)

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

W (0) = 1; W ( ) = W (- ); W ( ) 1.

(17)

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

Равномерная дискретность
В этом случае временные отсчеты образуют последовательность

tk = t k ,

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

(18)

где t = const. Подставляя (18) в (15) и (16), после несложных преобразований находим

W ( ) W0 ( , N , t) =
11

sin2 (N t/2) . N 2 sin2 ( t/2)

(19)


Рис. 3: Спектральное окно, соответствующее равномерно дискретному
распределению точек.

Важнейшим из свойств этой функции является ее периодичность с периодом 2c , где c = /t частота Найквиста. В свою очередь из периодичности спектрального окна (19) вытекают две важные особенности соответствующей периодограммы: ограниченность частотного диапазона и явление подмены частот (элайзинг). Перечислим еще два специфических свойства спектрального окна (19): 1) нули функции W0 расположены на частотах

k = 2 k /(N t),

k = +1, +2, . . . ;

(20)

2) ширина главного лепестка определяется величиной

= 21 = 4 /(N t).
Характерный вид спектрального окна (19) показан на рис. 3.

(21)

Случайные пропуски
Известно (Теребиж, 1992), что при случайных отклонениях моментов tk от равномерной сетки с шагом t спектральное окно, осредненное по множеству реализаций случайных отклонений, совпадает со спектральным окном, характерным для континуального заполнения наблюдениями интервала [t0 , tN -1 ]. Однако для каждого конкретного случая неравномерности (а именно с такими рядами мы встречаемся на практике) спектральное окно W ( ) имеет помимо 12


Рис. 4: Спектральное окно, соответствующее случайным отклонениям
моментов наблюдений от равномерной сетки.

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

Периодические пропуски
В астрономии, как правило, наблюдения производятся с периодическими или квазипериодическими пропусками. Приведем примеры нескольких типичных ситуаций: а) ночные наблюдения, которые длятся 2-6 часов подряд и при наличии ясной погоды разделяются соответственно перерывами в 2218 часов. В этом случае пропуски в наблюдениях имеют суточный период; б) фотометрические наблюдения переменных объектов, выполняемые обычно вдали от полнолуний, отчего в наблюдениях появляются пропуски с периодом, равным синодическому месяцу; в) позиционные и фотографические наблюдения планет, когда периодичность пропусков определяется взаимным расположением планеты, Земли и Солнца; 13

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

t постоянный шаг фиксации наблюдений; T = (n + p) t длительность цикла, в котором n первых точек заполнены наблюдениями, а p последующих точек образуют пропуск; m общее число циклов.
Рассмотрим теперь следующие временные окна:

w1 (t) =

1 n

n-1

[t - it],
i=0 m-1

(22)

w2 (t) =

1 m

[t - iT ].
i=0

(23)

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

w(t) =
-

w1 (t )w2 (t - t ) dt .

(24)

Вычислим теперь спектральное окно, соответствующее временному окну w(t). Подставляя (24) в (16) и используя теорему о преобразовании Фурье свертки двух функций, получаем

W ( ) = W0 ( , n, t)W0 ( , m, T ),
где

(25) (26)

W0 ( , N , ) =

sin2 ( N /2) . N 2 sin2 ( /2)

14

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

ми наблюдений. Штриховая линия график функции W0 ( , n, t) при n = 4, t = 1h ; тонкая линия график функции W0 ( , m, T ) при m = 5, T = 24h ; жирная линия результирующая кривая: последовательность наблюдений длиной в пять суток с ежесуточными наблюдениями, длящимися 4h .

Займемся исследованием функции W ( ). Мы видим, что она выражается через функцию (26), которая представляет собой стандартное окно для ряда из N точек, равномерно распределенных с шагом . Рассмотренные выше свойства функции W0 позволяют изучить сценарий развития спектрального окна W ( ) при накоплении наблюдений с пропусками. Характерные черты спектральных окон мы опишем в общих терминах, иллюстрации будут даны для типичной в астрономии ситуации ночные наблюдения внутри суток и перерывы наблюдений, исчисляемые сутками. a) Пусть проведено n равномерно распределенных наблюдений с шагом t. Полагая m = 1, из формулы (25) получаем

Wa ( ) = W0 ( , n, t).

(27)

Спектральные окна, соответствующие 4- и 6-часовым наблюдениям (n = 4, n = 6), выполненным с интервалом t = 1h внутри одних суток, показаны на рис. 5 и 6 штриховыми линиями. б) Теперь будем считать, что наблюдения, описанные в пункте a), 15


Рис. 6: Спектральное окно временного ряда с циклическими пропуска-

ми наблюдений. Штриховая линия график функции W0 ( , n, t) при n = 6, t = 1h ; тонкая линия график функции W0 ( , m, T ) при m = 5, T = 24h ; жирная линия результирующая кривая: последовательность наблюдений длиной в пять суток с ежесуточными наблюдениями, длящимися 6h .

повторены m раз. На основании формулы (25) имеем

Wb ( ) = W0 ( , n, t)W0 ( , m, T ).

(28)

Для исследования окна (28) заметим, что функция W0 ( , n, t) имеет главный пик шириной 4 /nt, а величина 2 /T равна периоду функции W0 ( , m, T ). В результате, в соответствии с формулой (28), получится сложный контур, в котором помимо главного пика на частоте = 0 будут видны убывающие по интенсивности пики на частотах

k = 2 k /T , ?

k = +1, +2, . . .

(29)

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

Pk = 2 /k = T /k , ?

k = 1, 2, . . .

(30)

Убывание интенсивностей пиков пропусков определяется огибающей Wa ( ). Таким образом, мы видим, что первый и наи16


более мощный пик пропусков лежит на частоте, соответствующей длительности цикла T . Всем остальным пикам соответствуют периоды T /2, T /3, ... На рис. 5, 6 жирными линиями показаны спектральные окна, соответствующие последовательности наблюдений длиной в пять суток (m = 5) с ежесуточными наблюдениями, длящимися по n = 4 и n = 6 часов. На этих рисунках мы видим пики пропусков, соответствующие периодам 24h , 12h , 8h и т. д., интенсивность которых зависит от числа наблюдений внутри суток. Рассмотренная картина образования спектральных окон временных рядов говорит о том, что положения всех пиков определяются цикличностью пропусков, а их интенсивность степенью заполненности циклов наблюдениями. Например, интенсивность первого пика определяется выражением

I1 = Wa (1 ) =

sin2 (n t/T ) . n2 sin2 ( t/T )

(31)

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

n1 , n2 число наблюдений в первом и втором сегментах, p число пропущенных наблюдений, образующих пропуск, M = n1 + n2 + p общее число точек, образующих весь ряд, N = n1 + n2 суммарное число наблюдений.
С учетом введенных обозначений для временного окна имеем n1 -1 M -1 1 (t - k t) + (t - k t) . (32) w(t) = N
k=0 k=n1 +p

17


Рис. 7: Спектральное окно двух восьмичасовых сеансов наблюдений, разделенных пропуском в 5 суток.

Подставляя (32) в (16), для соответствующего спектрального окна получаем

W ( ) =

n2 1 (n1 +n2 )

2

W 0 ( , n 1 , t ) +

n2 2 (n1 +n2 )2

W0 ( , n2 , t)+

W0 ( , n1 , t)W0 ( , n2 , t) cos(M + p) t/2. (33) Рассмотрим частный случай два одинаковых по длине сегмента, разделенных пропуском. Полагая в (33) n = n1 = n2 , находим 1 [1 + cos(n + p) t] W0 ( , n, t). (34) 2 Интересно отметить, что это же выражение можно получить из (28) при m = 2, поэтому можно сделать вывод, что структура спектральных окон при периодических пропусках и при существовании только одного пропуска фактически одинакова. В обоих случаях боковые пики располагаются на собственных частотах W ( ) = k = 2 k /(n + p)t, ? k = 1, 2, . . . ,
(35) где p длина пропуска. На рис. 7 показан пример спектрального окна, соответствующего длинному пропуску наблюдений. Заканчивая обзор спектральных окон для некоторых типичных в астрономии неравномерностей распределения наблюдений во времени, мы приводим на рис. 8 спектр скважности ряда наблюдений 18

2n1 n2 (n1 +n2 )2


Рис. 8: Спектр скважности мгновенных значений широты Пулкова
(ЗТЛ-180). Четырехчасовые наблюдения внутри суток, длина реализации 50.3 суток.

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

ћ период первого пика скважности соответствует длительности цикла "деньночь"; ћ периоды остальных пиков скважности совпадают с теоретическими значениями собственных периодов 12h , 8h , 6h ; ћ шумовой компонент в спектре скважности подтверждает нерегулярность выполнения ночных наблюдений.

19


LS-спектры
До сих пор, говоря об изучении спектров неравномерных рядов, мы основывались на том, что в качестве оценки спектра мощности выступает периодограмма Шустера. В работах Барнинга (1963), Ломба (1976), Скаргла (1982) для оценки спектра мощности неравномерных рядов предложены функции, отличные от традиционной периодограммы Шустера. В основе построения этих новых оценок лежит общая идея: аппроксимация временного ряда простой гармонической функцией. Рассмотрим подробнее эту задачу. Пусть временной ряд

xk = x(tk ),

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

задан на произвольном множестве временных отсчетов tk . В качестве модели этого ряда возьмем выражение
2

f (t) =
i=1

ai i (t),

(36)

где

1 (t) = cos t, 2 (t) = sin t.

(37)

Определим невязку аппроксимации

k = xk - f (tk )
и найдем коэффициенты модели (36) из условия

(38)



2

= min,

(39)

где использованы следующие обозначения:

(p, q ) =

1 N
2

N -1

p(tk )q (tk ),
k=0

(40) (41)

p

= (p, p).

Для определения коэффициентов a1 и a2 имеем систему нормальных уравнений

1 2 (2 , 1 )

(1 , 2 ) 2 2
20

a1 a2

=

(x, 1 ) (x, 2 )

,

(42)


решение которой имеет следующий вид:

a a
где

1 2

=

1

2 2 -(2 , 1 )
1 2

-(1 , 2 ) 1 2
2

(x, 1 ) (x, 2 )

,

(43) (44)

=



2

- [(1 , 2 )]2 .
2

Вычислим теперь величину

E

min

=

2

=x

2

-
i=1

ai (x, i )

(45)

и введем в рассмотрение функцию

P ( ) = x

2

- Emin ( ) 0,

(46)

где Emin ( ) имеет смысл минимума квадрата нормы невязки (45), а x 2 представляет собой дисперсию (полную мощность) ряда. Зафиксируем пробное значение частоты и вычислим коэффициенты a1 и a2 по формулам (43)(44). Очевидно, что значение функции P ( ) будет резко увеличиваться только в тех случаях, когда пробные частоты будут совпадать с частотами гармонических компонентов временного ряда xk , k = 0, 1, . . . , N - 1. По этой причине функцию P ( ) можно использовать в качестве оценки спектра мощности временных рядов. Поскольку эта оценка основана на методе наименьших квадратов, ее по предложению Ломба (1976) называют LS-спектром (LS Least Squares). Окончательная формула для вычисления LS-спектра следует из (46):

P ( ) =

1 2

2

ai (x, i ).
i=1

(47)

Множитель 1 введен сюда для удобства сравнения LS-спектра с пе2 риодограммой Шустера. Рассмотрим теперь конкретные представления LS-спектров.

LS-спектр Барнинга
Аналитическое выражение для этой функции, которую мы будем обозначать B ( ), следует из (47) при учете (43)-(44):

B ( ) =

1 2

1

2

(x, 2 )2 + 2 2 (x, 1 )2 - 2(1 , 2 )(x, 1 )(x, 2 ) , 1 2 2 2 - (1 , 2 )2 (48)
21


где



1

2

=

1 N

N -1

cos2 tk ,
k=0



2

2

=

1 N

N -1

sin2 tk ,
k=0

(49) (50)

(1 , 2 ) = (x, 1 ) = 1 N
N -1

1 N

N -1

sin tk cos tk ,
k=0

xk cos tk ,
k=0

(x, 2 ) =

1 N

N -1

xk sin tk .
k=0

(51)

LS-спектр Ломба
Для описания статистических свойств LS-спектра желательно иметь для него представление в виде суммы квадратов двух величин. Это можно сделать разными способами. Один из них был предложен Ломбом (1976). Суть этого метода заключается в том, что для каждой пробной частоты отсчеты tk смещают на величину ( ). Введем новые отсчеты времени

? tk = tk - ( ),

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

(52)

и наложим на величину (1 , 2 ) условие

(1 , 2 ) =

1 2

N -1

sin (tk - ( )) cos (tk - ( )) = 0.
k=0

(53)

Отсюда для величины ( ) находим
N -1

( ) =

1 arctg 2

sin 2 t

k

k=0 N -1 k=0

. cos 2 t
k

(54)

При сделанных предположениях для LS-спектра имеем следующее представление:

L( ) =
где

1 (x, 1 )2 ? (x, 2 )2 ? + , 2 2 1 ? 2 2 ?

(55) (56) (57)

? 1 (t) = cos( [tk - ( )]) cos( tk ), ? ? 2 (t) = sin( [tk - ( )]) sin( tk ). ?
22


Модификация Скаргла
Развивая идеи Ломба и желая сохранить информацию о фазе изучаемых колебаний, Скаргл (1982) ввел в рассмотрение специальный вид дискретного преобразования Фурье:

F T ( ) =

e

-i t



0

2

(x, 1 ) ? (x, 2 ) ? +i . 1 ? 2 ?

(58)

Легко понять, что периодограмма Ломба связана с преобразованием Скаргла следующим образом:

L( ) = |F T ( )|2 .

(59)

Сравнение рассмотренных нами LS-спектров показывает, что для любого временного ряда LS-спектры Барнинга (48) и Ломба (55) тождественны. Важно подчеркнуть, что эта тождественность имеет место для любых временных рядов, как равномерных, так и неравномерных.

LS-спектры и периодограмма Шустера
Во многих работах, где производилось сравнение классической периодограммы и LS-спектров, отмечается скорее их большое сходство, чем различие (Ломб, 1976; Карбонелл и др., 1992). Эти результаты, получаемые, как правило, для реальных наблюдений, представляют собой эмпирический факт и поэтому требуют объяснения. С этой целью мы проведем теоретический анализ рассмотренных выше спектральных оценок и сравним LS-спектры с периодограммой Шустера. Воспользуемся для представления LS-спектров формулой Барнинга (48), переписав ее в следующем виде:

B ( ) =
где

a2 ( )(x, 1 )2 + b2 ( )(x, 2 )2 - 2(x, 1 )(x, 2 )r , 1 - r2 a-2 ( ) = 2 b
-2 1 2 2

(60) (61) (62)

= 1 + Re (2 ), = 1 - Re (2 ),

( ) = 2 2

23


r=

(1 , 2 ) = 1 2

Im (2 ) [1 - Re 2 (2 )]

.

(63)

Мы определяем спектральное окно с помощью выражения

W ( ) = |( )|2 ,
где комплексное спектральное окно задается формулой

(64)

( ) =
Очевидно, что из условия

1 N

N -1

e
k=0

-i t

k

.

(65)

W (2 ) = 0
вытекает одновременное выполнение равенств

(66)

Re (2 ) = 0,

Im (2 ) = 0.

(67)

Кроме того, одновременное выполнение равенств (67) влечет за собой условие (66). Подставляя (66) и (67) в (60), имеем

B ( ) = (x, 1 )2 + (x, 2 )2 D( ).

(68)

Таким образом мы приходим к важному результату, который можно сформулировать следующим образом: если для конкретного распределения точек tk , k = 0, 1, . . . , N - 1, существуют такие частоты j , для которых выполняется условие W (2j ) = 0, (69)

то на этих частотах значения LS-спектров и периодограммы Шустера совпадают. Важным следствием этого результата является то обстоятельство, что степень сходства и различия LS-спектров и периодограммы Шустера определяется структурой спектрального окна. При этом следует ожидать, что существенное различие между LS-спектрами и периодограммой Шустера будет иметь место при таких распределениях точек tk , для которых спектральное окно W ( ) имеет четко выраженные боковые пики скважности. С другой стороны, эти отличия будут весьма малыми, если спектральное окно не имеет
24


Рис. 9: Локализация пиков функции W (2 ). Верхняя линия частоты,

на которых выполняется соотношение W (2 ) 0.01. Нижняя линия частоты, на которых W (2 ) > 0.01.

интенсивных пиков скважности. Рассмотрим подробнее этот вопрос для изученных выше типичных распределений временных точек.

Равномерное распределение точек
В этом случае для квадрата модуля спектрального окна (65), а также для его вещественной и мнимой части имеем

W ( ) W0 ( , N , t) = Re ( ) = Im ( ) =

sin2 (N t/2) , N 2 sin2 ( t/2)

(70) (71) (72)

sin(N t/2) cos ((N - 1) t/2) , N sin( t/2) sin(N t/2) sin ((N - 1) t/2) , N sin( t/2)

где t шаг выборки. Обычно значения периодограммы Шустера вычисляют на множестве фундаментальных частот

j =

2 j, N t

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

(73)

При нечетных N вычисления в формуле (73) ведутся до (N + 1)/2. Нетрудно показать, что эти частоты удовлетворяют соотношению (69). Таким образом, для равномерно распределенных временных рядов LS-спектры и периодограмма Шустера совпадают, если их значения вычислять на множестве фундаментальных частот.

25


Рис. 10: Периодограмма Шустера (a) и LS-спектр (b) на временной
сетке со случайно распределенными точками.

Случайные пропуски
Опыт показывает, что для случайного распределения наблюдений во времени спектральное окно W ( ) имеет центральный пик и не очень сильный шумовой компонент во всем спектральном диапазоне. Приведем в качестве примера рис. 9, где на верхней линии показаны частоты, для которых выполняется соотношение W (2 ) 0.01, а на нижней линии частоты, для которых это соотношение не выполняется. Как мы видим, последних точек не так много. По этой причине следует ожидать, что при случайном распределении наблюдений во времени вид LS-спектров не будет сильно отличаться от периодограммы Шустера. Это подтверждает рис. 10, на котором показаны периодограмма Шустера и LS-спектр, полученные по 50 точкам, случайным образом отклоняющимся от равномерно распределенных точек, для функции (80) с параметрами A1 = A2 = 1, 1 = 1.1, 2 = 3.3 Гц.

26


Рис. 11: Совпадение периодограммы Шустера (a) с LS-спектром (b) на
сетке с периодическими пропусками.

Периодические пропуски
Пусть на сетке с постоянным шагом t после каждых n наблюдений происходит p пропусков. Обозначим через m число групп n + p. Очевидно, что период пропусков составляет T = (n + p)t. Ранее было показано, что спектральное окно для такой временной сетки имеет вид W ( ) = W0 (n, , t)W0 (m, , T ), (74) где функция W0 , задаваемая выражением (70), есть временное окно, соответствующее равномерному распределению точек. Спектральное окно (74) имеет систему боковых пиков на собственных частотах

2 l , l = 1, 2, . . . (75) T Если в наших m группах, состоящих из n наблюдений и p пропусков, заполнить все пропуски наблюдениями, то для получившегося l = ?
27


Рис. 12: Отличие периодограммы Шустера (a) от LS-спектра (b) на
сетке с периодическими пропусками.

равномерного ряда, состоящего из m(n + p) точек, можно ввести систему фундаментальных частот

j =

2 j, m T

j = 1, . . . , m(n + p)/2.

(76)

Исключим отсюда частоты, кратные половине первой собственной частоты, т. е. образуем числа
j =

2 j, m T

(77)

где

j = 1, 2, . . . , m/2 - 1, m/2 + 1, . . . , m - 1, m + 1, . . . , m(n + p)/2 - 1.
Простой подстановкой теперь легко убедиться, что на этих частотах выполняется условие
W (2j ) = 0.

(78)

28


На этом основании можно утверждать, что при периодических пропусках наблюдений с частотой 1 = 2 /T , где T период про? пусков, периодограмма Шустера и LS-спектр будут одинаковы, если ни одна из частот гармонических компонентов исследуемого временного ряда не совпадает с частотами

j = ^

j, T

j = 1, 2, . . .

(79)

Например, для сетки с параметрами t = 0.1 с, n = 3, p = 7, m = 15, собственные частоты спектрального окна W ( ) образуют последовательность 1 = 1, 2 = 2, . . . , 5 = 5 Гц. На этой временной сетке ? ? ? были вычислены значения функции

x(t) = A1 cos 2 1 t + A2 cos 2 2 t.

(80)

На рис. 11 показаны периодограмма Шустера и LS-спектр для значений параметров A1 = A2 = 1, 1 = 1.1 Гц, 2 = 3.3 Гц. Мы видим, что периодограмма Шустера совпала с LS-спектром, поскольку значения 21 и 22 не совпадают ни с одним из значений j , j = 1, 2, . . . , 5. ? Противоположный случай показан на рис. 12, на котором вновь вычислены наши периодограммы для значений A1 = A2 = 1; 1 = 0.5 Гц, 2 = 3.3 Гц. Поскольку теперь имеет место равенство 21 = 1 = 1 Гц, периодограмма Шустера существенно отличается ? от соответствующего LS-спектра.

29


Рис. 13: Совпадение периодограммы Шустера (a) с LS-спектром (b) на
сетке длинный пропуск.

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

W ( ) = W0 ( , n, t)[1 + cos T ]/2,

(81)

где T = (n + p)t. Нетрудно показать, что условию W (2j ) = 0 удовлетворяют частоты
j =

T

j+

1 2

,

j = 0, 1, 2, . . . , n + p - 1.

(82)

Важно отметить, что собственные частоты спектрального окна

k = ?

2 k, T

k = 1, 2 . . . ,

(83)

не совпадают с частотами 2j .

30


Рис. 14: Отличие периодограммы Шустера (a) от LS-спектра (b) на
сетке типа длинный пропуск.

На рис. 13 показаны периодограмма D( ) и соответствующий LS-спектр функции (80) с параметрами A1 = A2 = 1, 1 = 1.1 Гц и 2 = 3.3 Гц, вычисленные на временной сетке с пропуском (n = 40, p = 40, t = 0.1 Гц). Обе периодограммы оказались одинаковыми, поскольку частоты 2 1 и 2 2 не совпадают ни с одной из частот j . Однако если взять 1 = 0.0625 Гц, 2 = 3.3 Гц, то наши периодограммы вновь окажутся различными (рис. 14), так как в этом примере удвоенная частота гармонического компонента выбрана равной частоте первого (и, следовательно, наибольшего) бокового пика.

Свойства периодограмм белого шума
Знание закона распределения отсчетов периодограммы белого шума является фундаментом для решения задачи о выделении сигнала из шума. Хорошо известно (Дженкинс и Ваттс, 1972), что для равномерно распределенных рядов отсчеты периодограммы распре31


делены по экспоненциальному закону. В настоящем параграфе этот вопрос будет пересмотрен для рядов с произвольно распределенными моментами времени как для LS-спектров, так и для периодограммы Шустера. Примем, что центрированный ряд {xk , tk }, k = 0, 1, 2, . . . , N - 1, представляет собой случайную величину, значения которой распределены по нормальному закону. Для такого ряда имеем

xp xq =

2 0 , p = q , 0, p = q ,

(84)

2 где 0 дисперсия шума (полная мощность). Приведем без доказательства основные положения, которые определяют законы распределения отсчетов периодограмм дискретного белого шума и степень их коррелированности. Доказательства этих положений содержатся в работах Ломба (1976), Скаргла (1982) и Витязева (1995).

Периодограмма Ломба
Значения нормированной периодограммы Ломба

L( ) 2 0 /N имеют экспоненциальное распределение l ( ) = p(z ) = e
-z

(85)

,

0 z < ,

= 0.

(86)

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

? ? xk = a cos 1 tk + b sin 1 tk ,

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

(87)

где a и b центрированные случайные величины. Пусть
( ) = (
k

L( ) , x2 )/N k

(88)

32


тогда коэффициент корреляции между частотами 1 и определяется формулой k (1 , ) = ( ). (89)

Периодограмма Шустера
Значения нормированной периодограммы Шустера

d( ) =
распределены по закону

D( ) 2 0 /N = 0,

(90)

p(z ) = e
где

-z

R(z , ), 0 z < ,

(91) (92) (93)

R (z , ) =

- 1 /2

exp[-z (1 - )/] I0 [-z (1 - )1/2 /], = 1 - W (2 ).

В уравнении (92) через I0 (x) обозначена модифицированная функция Бесселя первого рода. Для определения коррелированности отсчетов периодограммы Шустера может служить следующий результат (Витязев, 1995).

Пусть центрированный гауссовский белый шум задан на множестве произвольных временных отсчетов. В этом случае коэффициент корреляции между отсчетами периодограммы Шустера на частотах 1 и 2 имеет вид
k (1 , 2 ) = W (1 - 2 ) + W (1 + 2 ) 1 + W (21 ) 1 + W (22 ) .
(94)

Отметим, что если отсчеты периодограммы Шустера отнести к частотам j , j = 1, 2, ..., для которых выполняются соотношения

W (2j ) = 0, W (i + j ) = 0, i, j = 1, 2, ..., i = j,

(95) (96)

33


то эти отсчеты будут иметь экспоненциальное распределение и будут некоррелированными (независимыми). В анализе неравномерных рядов они играет ту же роль, какую фундаментальная система частот играет в теории и практике обработки равномерных рядов. Например, для рядов с периодическими пропусками такие частоты можно найти на любом интервале, не содержащем внутри себя точки 0.5 1 k , k = 1, 2, ..., где через k обозначены, как и ранее, соб? ? ственные частоты спектрального окна. В общем случае построение такой системы частот едва ли возможно, однако для практических целей вполне достаточно иметь квазифундаментальную систему частот (КФСЧ), для которой соотношения (95) и (96) выполняются приближенно. В отличие от регулярных рядов, для которых верхняя граница частотного интервала фиксируется найквистовой частотой c = /t, соответствующая граница для нерегулярных рядов не имеет однозначного определения. Часто ее назначают как max = / , где наименьшее или среднее расстояние между временными отсчетами. Внутри выбранного промежутка отсчеты периодограммы вычисляют на множестве частот, следующих друг за другом с постоянным шагом. Очевидно, что такая система частот может играть роль КФСЧ только в том случае, когда наибольший пик скважности соответствующего спектрального окна много меньше единицы. Такая ситуация осуществляется для временных отсчетов, случайным образом отклоняющихся от равномерной сетки, или для временных рядов с пропусками, если длина пропусков меньше половины длины ряда.

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

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

34

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

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

Равномерные временные ряды
Здесь считается, что временные точки задаются последовательностью tk = tk , k = 0, 1, . . . , N - 1, t = const. Ранее было показано, что на фундаментальной системе частот

k =

2 k , N t

k = 0, 1, . . . ,

N - 1, 2

(97)

периодограмма Шустера совпадает с LS-спектром. Добавим к этому, что суммы, разности и удвоенные фундаментальные частоты также являются фундаментальными частотами. Легко проверить, что при равномерном распределении точек выполняются равенства

W (k ) = 0,

k = 1, 2, . . . ,

N - 1, 2

(98)

35


поэтому в силу (94) отсчеты периодограмм на системе фундаментальных частот являются некоррелированными. Кроме того, отсчеты нормированной периодограммы

d(k ) =

D(k ) 2 0 /N N - 1. 2

(99)

распределены по экспоненциальному закону

p(z ) = e

-z

,

k = 1, 2, . . . ,

(100)

В ситуации H1 имеем: 1) вероятность того, что отсчет d(k ) при любом k не превосходит величины x, равна
x

e
0

-z

dz = 1 - e

-x

,

x 0;

(101)
N 2

2) вероятность того, что все независимые d(k ), k = 1, 2, . . . , не превосходят величины x, равна

- 1,

(1 - e

-x N/2-1

)

;

(102)

3) вероятность того, что хотя бы один отсчет d(k ) превзойдет величину x, определяется выражением

Q(x) = 1 - (1 - e

-x N/2-1

)

.

(103)

Мы получили функцию распределения наибольшего из отсчетов периодограммы случайного процесса в главном частотном интервале (распределение Уолкера, 1914). Решение уравнения

Q(X1 ) = q ,

(104)

где X1 порог обнаружения сигнала с вероятностью 1 - q , приводит к выражению X1 = - ln[1 - (1 - q )2/(N -2) ]. (105) Таким образом, те наибольшие отсчеты периодограммы, для которых выполняются неравенство

dmax > X1 ,
36

(106)


с вероятностью 1 - q могут считаться отсчетами, принадлежащими сигналу, а не шуму. В ситуации H2 легко получить вероятность того, что отсчет периодограммы на данной частоте превышает заданный уровень x:


Q(x) =
x

e

-z

dz = e

-x

.

(107)

Задавая вероятность редкого события q << 1 (уровень значимости), из уравнения q = e-X2 (108) получаем критический порог обнаружения сигнала в виде

X2 = - ln q .
Поэтому при выполнении неравенства

(109)

d(j ) > X

2

(110)

мы можем с вероятностью 1 - q подтвердить исходное предположение о том, что отсчет d(j ) принадлежит сигналу, а не шуму.

Неравномерные ряды
Как следует из (86), для отчетов нормированной периодограммы Ломба имеет место экспоненциальный закон распределения, однако, как было сказано выше, отсчеты периодограммы коррелированы, и, следовательно, зависимы. По этой причине переход от (101) к (102) в ситуации H1 неправомерен, что мешает получить строгое аналитическое выражение для функции распределения наибольшего отсчета. В равной степени все это касается и периодограммы Шустера. Вероятно, наиболее прямой путь получения функции Q(x) состоит в проведении численных экспериментов по построению периодограмм квазигауссовского белого шума с последующим выводом эмпирической функции распределения наибольших отсчетов. Сошлемся здесь на работу Хорна и Балиунаса (1986), в которой вместо соотношения (103) предлагается использовать эмпирическую формулу Q(x) = 1 - (1 - e-x )Ni , (111)

37

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

Ni = -6.362 + 1.193N + 0.00098N 2 .

(112)

Достаточно очевидно, что это соотношение нельзя считать универсальным, поэтому связь между Ni и N следует определять численными экспериментами для каждого неравномерного временного ряда. Что же касается ситуации H2 , то здесь для периодограммы Ломба критический порог обнаружения сигнала определяется формулой (109) при всех видах распределений точек tk . В случае периодограммы Шустера аналогом (107) выступает формула


Q(x) =
x

p(z ) dz ,

(113)

а величина X2 является корнем уравнения


q=
X2

e

-z

R(z , ) dz ,

(114)

где функция R(z , ) задается выражениями (92) и (93). Следует сказать, что если отсчеты периодограммы вычислять не на произвольной сетке частот, а на квазифундаментальной системе частот (см. с. 34), то ситуация резко меняется. Действительно, на этой системе частот различие между периодограммой Шустера и LS-спектрами весьма мало, отсчеты периодограммы становятся практически некоррелированными и распределенными по экспоненциальному закону. Это означает, что для отсчетов периодограммы Шустера на квазифундаментальной системе частот можно назначать порог обнаружения сигнала по формулам (105) и (109) в зависимости от того, пользуемся мы статистикой наибольшего отсчета (ситуация H1 ) или статистикой отсчета с заранее известной частотой (ситуация H2 ).

38

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

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

39

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

A2 A2 A2 [1+W (20 )]+ Re (20 ) cos(2)- Im (20 ) sin(2), 4 2 2 (115) в то время как для LS-спектра это выражение имеет другой вид: D(0 ) = B (0 ) = A2 [1 + Re (20 ) cos(2) - Im (20 ) sin(2)]. 4
(116)

Допустим, что частота нашей гармоники удовлетворяет условию

1 - W (20 ) << 1.

(117)

Ниже в таблице показаны значения D(0 ) и B (0 ) как функции фазы для синусоиды с параметрами A = 1, 0 = , вычисленной на сетке временных точек с пропусками (n = 3, p = 7, m = 20, t = 0.1 c). В этом случае W (20 ) = 0.762; Re (20 ) = 0.706; Im (20 ) = 0.513, так что условие (117) выполняется.

0 0 45 90 135 180

D(0 ) 0.793 0.184 0.087 0.697 0.793

B (0 ) 0.426 0.122 0.074 0.378 0.426

Из этой таблицы видно, что при выполнении неравенства (117) ни периодограмма Шустера, ни LS-спектр не могут применяться в качестве оценки спектра мощности, так как искомые интенсивности далеки от истинного значения 0.25 и показывают сильную зависимость от фазы. Для устранения этой зависимости в общем случае наложим на частоты p , p = 1, 2, ..., n, полигармонической функции следующие условия: 40


W (2p ) = 0, W (p + q ) = 0, p, q = 1, 2, ..., n, p = q .

(118) (119)

Из (12) видно, что при выполнении этих условий член D0 ( ) становится равным нулю, и структура периодограммы будет определяться первым слагаемым в уравнении (11), т. е. сверткой истинного спектра мощности и спектрального окна. Выясним теперь смысл соотношений (118) и (119). Очевидно, что при выполнении этих условий каждый гармонический компонент полигармонической функции отобразится в периодограмме в виде контура, форма которого определяется спектральным окном W ( ), причем боковые пики контуров не будут накладываться друг на друга и, следовательно, интенсивности как истинных, так и ложных пиков не будут искажаться (см. рис. 11 и 13). В этой ситуации возможно точное восстановление спектра (очищение периодограммы от ложных линий). По этой причине соотношения (118) и (119) следует назвать условиями восстановления истинного спектра полигармонической функции. Другими словами, определение частот и амплитуд полигармонической функции можно произвести точно, если удвоенные частоты полигармонической функции, а также их суммы и разности не совпадают с собственными частотами спектрального окна. Наоборот, нарушение хотя бы одного из условий (118) и (119) приведет к тому, что боковые пики контуров будут накладываться друг на друга (см. рис. 17 и 18). При этом, вообще говоря, произойдет искажение и истинных, и ложных пиков периодограммы. Если взаимодействие контуров линий происходит на слабых пиках скважности, то эти искажения будут не очень сильными и спектр может быть восстановлен лишь с некоторыми искажениями. Степень искажения интенсивностей пиков периодограммы зависит от фаз гармонических компонентов и от соотношения наблюдений и пропусков в общей длине временного ряда. Опыт показывает, что вполне приемлемые результаты могут быть получены даже при нарушении условий (118) и (119), если общая длина пропусков составляет не более половины ряда. Если спектральное окно имеет сильные пики скважности, то их наложения могут привести к различным патологиям: истинные пики могут стать менее интенсивными, чем ложные, некоторые пики на частотах p , p = 1, 2, ..., n, могут исчезнуть совсем, возможно 41


также появление лишних пиков и т. д. Вполне понятно, что восстановление истинного спектра в подобных ситуациях может оказаться невозможным. Следует сказать, что принципиальная невозможность восстановления амплитуд и частот полигармонической функции при нарушении условий (118) и (119) происходит не по причине несовершенства периодограммы, а из-за того, что амплитудно-частотный состав исследуемой функции может претерпеть изменения еще во временной области в результате пропусков. Дело в том, что термин "грязная"периодограмма относится к периодограмме функции x(t). По отношению же к произведению w(t)x(t) "грязная"периодограмма является ч и с т о й периодограммой, т. е. она показывает, что на самом деле произошло с нашей функцией после ее умножения на временное окно, зависящее от распределения пропусков наблюдений во времени. Приведем два примера. 1) Зададим P -периодическое спектральное окно в следующем виде: 1 2 j 1 [1 - cos j ] sin t. (120) w(t) = + 2 j =1 j P На временной сетке с постоянным шагом эта функция принимает значения 1 1 1 1 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 ... 2 2 2 2 Если исходная функция имеет представление

x(t) =

4 2 4 sin t + cos t, 3 P P

(121)

то для наблюденной функции имеем

y (t) = w(t) x(t) =

4 1 8 4 6 6 + ( - 2 ) cos t+ sin t... 3 2 2 9 P 5 P

(122)

Сравнивая (121) с (122), убеждаемся, что в наблюденной функции исчезла гармоника на частоте 2 /P. Очевидно, что эта потеря невосполнима, так как спектральному анализу подвергается функция (122), а не исходный процесс. 2) Теперь возьмем другое P -периодическое спектральное окно

1 w(t) = + 3

j =1

1 2 j 2 j 2 j 2 j (1 - cos ) sin t + sin cos t , (123) j 3 P 3 P
42


которое на дискретной сетке с постоянным шагом имеет вид

1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 ... 2 2 2 2
Задавая исходную функцию в виде

x(t) = sin

2 6 t - cos t, P P

(124)

для наблюденной функции получим

y (t) = w(t) x(t) = = 0.75 + 1.599 sin( + 1.385 sin(
6 t P 2 t P

+ 1 ) + 1.652 sin(
8 t P

4 t P

+ 2 )+

(125)

+ 3 ) + 0.964 sin(

+ 4 ) + . . .

Этот результат интересен тем, что в наблюденной функции появилась новая гармоника с частотой 4 /P, в то время как в исходной функции имеются только две гармоники с частотами 2 /P и 6 /P . Очевидно, что интерпретация "грязного"спектра, у которого наибольший пик является "духом", обречена на неудачу (если только у нас нет априорной информации о его ложности). В астрономии вращение Земли и ее обращение вокруг Солнца приводят к появлению суточных и годичных пропусков в наземных наблюдениях. Вряд ли можно считать, что периодичности, присущие различным процессам на Солнце, звездах, квазарах и т. п., могут быть как-то связаны с суточными и годичными периодами, свойственными Земле. По этой причине обнаружение нарушения условий (118) и (119) в наблюдениях указанных светил скорее всего свидетельствует о систематических воздействиях на результаты наблюдений со стороны, скажем, атмосферных явлений. Другое дело, если мы изучаем физические явления, присущие Земле, с самой Земли. Именно с этой ситуацией мы встречаемся в случае астрометрических наблюдений параметров вращения Земли, когда годичные и полугодичные колебания в неравномерности вращения Земли взаимодействуют с годичными пропусками наблюдений. Полученные здесь условия восстановления спектрального состава временных рядов, заданных на неравномерных сетках, будут использованы в качестве критериев правильности работы алгоритма CLEAN. 43


Метод CLEAN и его модификации
Описание метода
В настоящее время термином CLEAN обозначают два алгоритма обработки данных. Первый был предложен Хогбомом (1974) для чистки "грязных карт", получаемых при апертурном синтезе в радиоастрономии, второй является одномерной версией указанного метода и используется для получения "чистых"спектров в спектральном анализе временных рядов (Робертс и др., 1987). В дальнейшем мы будем использовать этот термин во втором смысле. Опишем вкратце основную идею метода CLEAN. Пусть имеется функция f = f (t), - < t < , значения которой могут измеряться на произвольном дискретном множестве точек tk , k = 0, 1, ћ ћ ћ , N - 1. Представим результаты наблюдений в виде d(t) = w(t)f (t), (126) где весовая функция w(t) может быть определена следующим образом: N -1 1 w(t) = (t - tk ). (127) N
k=0

Введем в рассмотрение следующие преобразования Фурье:
+

F ( ) =
- +

f (t)e

-i2 t

dt,

(128) (129) (130)

W ( ) =
- +

w(t)e d(t)e
-

-i2 t

dt,

D( ) =

-i2 t

dt.

Используя терминологию предыдущих параграфов, будем называть F ( ) спектром исходной функции f (t), W (t) спектральным окном, а D( ) "грязным спектром", однако в этом параграфе все эти функции являются комплекснозначными. По свойству преобразования Фурье произведения двух функций имеем
+

D( ) =
-

F ( )W ( - ) d ,

(131)

44


где

D( ) =

1 N

N -1

fk e-
k=0 N -1

i 2 t

k

,

(132) (133)


W ( ) =

1 N

e
k=0

-i2 t

k

.

При этом имеют место следующие соотношения (символом значена операция комплексного сопряжения):

обо(134) (135) (136)

D(- ) = D ( ), W (- ) = W ( ), W (0) = 1.

Допустим, что исходная функция f = f (t) является простейшей гармонической функцией

f (t) = A cos(2 t + ), ?
тогда ее спектр имеет вид

(137)

F ( ) = a ( - ) + a ( + ), ? ?
где комплексная амплитуда a определяется выражением

(138)

a=

A i e. 2

(139)

Подставляя (138) в (131) и пользуясь свойствами дельта-функции, находим D( ) = aW ( - ) + a W ( + ). ? ? (140) Нетрудно показать, что при известной частоте амплитуда A и фа? за исходного сигнала (137) восстанавливаются однозначно. Действительно, на основании (137) при учете (134), (135), (136) можно записать следующую систему уравнений:

D( ) = a + a W (2 ), ? ? D ( ) = aW (2 ) + a , ? ?
откуда следует, что

(141)

a=

D( ) - D ( )W (2 ) ? ? ? . 1 - |W (2 )|2 ?
45

(142)


Искомые амплитуда A и фаза теперь определяются из формул

A = 2|a|, = arctan

(143)

Im a , (144) Re a при этом четверть, в которой лежит фаза , находится по знакам Im a и Re a. Основная идея чистки спектра в методе CLEAN заключается в вычитании из "грязного"спектра D( ) вклада каждого гармонического компонента, определяемого формулой (140). Необходимая для этого величина приближенно оценивается по пикам пери? одограммы |D( )|2 , после чего значения a и a вычисляются по формуле (142). Для увеличения устойчивости алгоритма обычно из исходного спектра вычитают не весь вклад гармоники (140), а лишь его некоторую долю, определенную параметром g , 0 < g 1. Образующиеся при этом остатки, а также вклады от других гармоник убираются аналогичным образом в процессе итераций. Отметим, что введение в алгоритм параметра g не является принципиальным, поскольку производимую с его помощью компенсацию приближенного оценивания частоты можно заменить соответ? ствующим уменьшением частотного шага при вычислении функций D( ) и W ( ). В процессе последовательного съема гармоник образуется искусственный спектр сверхразрешения, состоящий из значений g D(j ), отнесенных к тем частотам j , на которых пери? ? одограмма |D( )|2 имела истинные пики. С помощью некоторого "чистого"спектрального окна B ( ), ширина которого определяется величиной порядка (tN -1 - t0 )-1 , искусственный спектр сверхразрешения преобразуется в "чистый"спектр с реальным разрешением. Интересной особенностью метода CLEAN является также то, что обратное преобразование Фурье чистого спектра позволяет "восстановить"значения функции f (t) в пропущенных точках. Авторы метода CLEAN продемонстрировали высокую эффективность своего метода на примерах обработки искусственных временных рядов без шума. Тем не менее, практика использования алгоритмов CLEAN в радиоастрономии (Томпсон, 1986) и в анализе временных рядов (Карбонелл и др., 1992) показала, что в отдельных случаях эти алгоритмы дают неверные результаты. По нашему мнению, это происходит по двум причинам:
46

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

1) Критерий применимости. Метод CLEAN основан на том, что наибольший отсчет в периодограмме |D( )|2 всегда соответствует истинному пику. Выше было показано, что при нарушении условий (118) и (119) наибольший отсчет периодограммы может соответствовать и ложному пику. В этом случае метод CLEAN сделает попытку вычесть из спектра несуществующий компонент, что приведет в конечном итоге не к чистке, а к дополнительному засорению спектра. Кроме того, при нарушении условий (118) и (119) может произойти потеря некоторых гармонических компонентов, в результате чего метод CLEAN их просто не найдет. Нетрудно показать, что нарушение условий (118) и (119) можно обнаружить при сравнении периодограммы |D( )|2 со спектральным окном |W ( )|2 . Если окажется, что все пики "грязной"периодограммы совпадают с пиками спектрального окна или симметрично смещены относительно них на одну и ту же величину, то у нас вообще нет гарантии, что спектральный состав исходной функции можно восстановить полностью, поэтому в подобных ситуациях к результатам чистки, в том числе и методом CLEAN, следует относиться критически. 2) Условие окончания итераций. Реальные наблюдения всегда зашумлены, поэтому в периодограмме |D( )|2 помимо ложных
47

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

Dq = D Xq ,

(145)

где Xq порог обнаружения сигнала в шумах с заданной вероятностью 1 - q ; D среднее значение периодограммы |D( )|2 . В ситуации H1 величину Xq следует назначить по результатам численных определений вероятности распределения наибольшего отсчета периодограммы белого шума. В ситуации H2 можно воспользоваться формулой (114). Разумеется, эти вычисления нужно проделывать для каждого конкретного расположения временных отсчетов. В отдельных случаях (см. с. 38) для вычисления порога можно использовать формулы (105) и (109).

Дополнительные возможности алгоритма CLEAN
Как было показано выше, работа метода CLEAN основана на получении искусственного спектра сверхразрешения, свертка которого со спектральным окном, не имеющим боковых пиков скважности, дает "чистый"спектр. Авторами метода было отмечено, что обращение "чистого"спектра позволяет восстановить значения функции f (t) в пропущенных точках. В дополнение к этому мы укажем еще на две задачи, которые можно решить методом CLEAN. а) Фильтрация данных. Как известно, целью фильтрации является преобразование исходного временного ряда путем изменения его спектрального состава. На практике используют низкочастотные, высокочастотные и полосовые фильтры, позволяющие оставить в спектрах исходных временных рядов спектральные линии в соответствующих участках спектра. Известно также, что при построении фильтров с помощью операции свертки на ко48

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

б) Коррелограмма. Получение коррелограммы неравномерного ряда это особая проблема, поскольку непосредственное вычисление по стандартным формулам требует предварительного осреднения данных на интервалах равной длины, а обращение по Фурье "грязного"спектра неминуемо даст "грязную"коррелограмму. Подобный метод обращения спектра Ломба предложен в работе Скаргла (1989), однако этот подход не является строгим, так как LS-спектр и коррелограмма не связаны между собой по Фурье. Однако если периодограмма очищена от "духов", то ее обращение по Фурье это корректный путь получения реальной оценки корреляционной функции неравномерного временного ряда. Очевидно, что все необходимые для этого условия реализуются в алгоритме CLEAN на той его стадии, когда получен "чистый"спектр с реальным разрешением, поэтому обращение "чистой"периодограммы позволяет очень просто получить искомую коррелограмму на равномерной сетке значений временного сдвига. Заметим, что такая коррелограмма не отягощена ошибками кругового эффекта, так как обычная мера борьбы с этим эффектом дополнение ряда нулями, приводящее к интерполяции в области частот, здесь выполняется автоматически за счет выбора соответствующего шага по частоте при вычислении всех спектральных оценок.

49


Модифицированный алгоритм CLEAN
С учетом вышесказанного модифицированный алгоритм CLEAN может быть представлен следующим образом. Исходные данные:

N число точек исходного временного ряда; f [k ] массив значений исходного временного ряда; t[k ] массив временных отсчетов k = 0, 1, . . . , N - 1.
Управляющие параметры: g доля вклада гармоники в "грязный"спектр (0 < g 1); Xq порог обнаружения сигнала в шумах с вероятностью 1 - q (см. с. 48). Исполнение алгоритма: 1. Найти
min

минимальное расстояние между точками tk .

2. Определить верхнюю границу диапазона частот

max =

1 . 2min

(146)

Замечание: может случиться, что найденное минимальное расстояние между отсчетами приведет к неразумному завышению верхней границы частотного диапазона. В этом случае величину max можно определить, исходя из среднего расстояния между точками t[k].
3. Вычислить параметр m, определяющий границы массивов спектральных оценок

m = max (t

N -1

- t0 ).

(147)

В работе Робертса и др. (1987) принято = 4, при желании это значение можно увеличить. 4. Вычислить массив значений "грязного"спектра

D [j ] =

1 N

N -1

fk e-
k=0

i 2 j t

k

, j = -m, . . . , 0, . . . , m.

(148)

50


5. Вычислить массив значений спектрального окна

1 W [j ] = N

N -1

e
k=0

-i2 j t

k

, j = -2m, . . . , 0, . . . , 2m.

(149)

В пп. 4 и 5 частоты j определены следующим образом:

j = (j /m)max .

(150)

6. Заполнить нулями элементы массива спектра сверхразрешения

C [j ] = 0, j = -m, . . . , 0, . . . , m.
7. Вычислить порог обнаружения сигнала в шумах

(151)

Dq = D Xq ,
где

(152)

D=

1 m+1

m

|D[j ]|2 .
j =0

(153)

8. Положить i = 1 (номер итерации). Изменяющийся в процессе итераций массив "грязного"спектра будем обозначать D(i) [j ], при этом D(0) [j ] D[j ]. 9. Определить в массиве |D(0) [j ]|2 наибольший отсчет D ветствующий ему номер элемента j = J > 0.
max

и соот-

10. Если Dmax < Dq , перейти к п. 16, в противном случае перейти к исполнению следующего шага. 11. Вычислить значение комплексной амплитуды

a=

D

(i-1)

[J ] - (Di-1 [J ]) W [2 J ] . 1 - |W [2 J ]2

(154)

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

D(i) [j ] = D

(i-1)

[j ] - g (aW [j - J ] + a W [j + J ]) ,
(155)

j = -m, . . . , 0, . . . , m.
51


13. Записать вклад снятой гармоники в массив C :

C [J ] = C [J ] + g a, C [-J ] = C [-J ] + g a .
(156)

14. Определить в периодограмме |D(i) [j ]|2 наибольший отсчет D и соответствующий ему номер элемента j = J > 0. 15. Перейти к п. 10, увеличив i на единицу.

max

16. Если i > 1, перейти к п. 17, в противном случае напечатать сообщение Исследуемый ряд с вероятностью 1-q не содержит гармонических компонентов и перейти к п. 22. 17. При необходимости выполнения фильтрации временного ряда соответствующим элементам массива C [i] присвоить нулевые значения. 18. Вычислить "чистое"спектральное окно, соответствующее равномерному распределению точек

1 B [j ] = N
где

N -1

e
k=0

-i2 j tk

, j = -2m, . . . , 0, . . . , 2m,

(157)

tN -1 - t0 k. (158) N -1 В оригинальном методе в качестве B используется функция, имеющая форму гауссианы. tk =
19. Вычислить "чистый"спектр
m

S [j ] =
k=-m

C [k ] B [j - k ], j = -m, . . . , 0, . . . , m.

(159)

20. Вычислить коррелограмму (смещенную оценку корреляционной функции) по формуле

(tk ) = ?

N m

m

|S [j ]|2 e
j =-m

i 2 j t

k

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

(160)

52


21. Вычислить временной ряд на равномерной сетке точек

f (tk ) =

N m

m

S [j ]e
j =-m

i 2 j t

k

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

(161)

22. Конец алгоритма.

Комментарии к исполнению алгоритма
Метод CLEAN предназначен для очищения "грязной"периодограммы от эффектов иррегулярности наблюдений, а также от действия шумов. Перечислим некоторые особенности процедуры чистки, которые надо иметь в виду при работе с программой метода CLEAN. 1. Программу метода CLEAN можно использовать для анализа равномерных рядов. В этом случае происходит очищение периодограммы от шумового компонента. Заметим, однако, что для равномерных рядов существуют более производительные программы, основанные на использовании процедур быстрого преобразования Фурье. 2. При обработке методом CLEAN рядов, у которых временные отсчеты случайно отклоняются от равномерной сетки и не образуют длинных пропусков, происходит очищение "грязной"периодограммы как от эффектов шума, так и от нерегулярности измерений. Для таких рядов, как правило, существует квазифундаментальная система частот и выполняются условия (118) и (119) восстановления истинного спектра. Поэтому никаких серьезных ограничений на применимость метода чистки теоретически не возникает. Надежность чистки таких рядов подтверждается многочисленными примерами. 3. Наиболее трудными для чистки являются ряды с периодическим пропусками или ряды с длинным пропуском. Дело в том, что для таких рядов из-за коррелированности отсчетов периодограммы дискретного белого шума существуют большие трудности для установления критерия выделения сигнала из шума. Кроме того, искомые частоты полигармонической функции могут не удовлетворять условиям (118) и (119). При слабой зашумленности ряда нарушения этих условий можно обнаружить по расположению линий "грязной"периодограммы: линии либо расположены на собственных частотах, либо сдвинуты относительно собственных частот на одну и ту же величину, либо 53

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

Численные примеры
Для иллюстрации возможностей алгоритма CLEAN мы применили его для анализа неравномерных рядов искусственного и естественного происхождения. На рис. 15 и 16 показаны результаты спектрального анализа временных рядов, состоящих из двух косинусоид и шума при отношении сигнал/шум, равном единице. В обоих примерах данные относятся к неравномерным сеткам временных отсчетов, однако условия восстановления истинного спектра (118) и (119) выполняются. В этом случае, как мы видим, метод CLEAN работает успешно, так как при переходе от "грязной"периодограммы к "чистой"ложные пики исчезают, а истинные пики восстанавливаются полностью. На этих же рисунках приводятся коррелограммы, а на рис. 16 мы показали, как можно использовать алгоритм CLEAN для низкочастотной фильтрации (при переходе от рис. 16,c к рис. 16,d спектральная линия на частоте 2.9 Гц была уничтожена). В противоположность этому рис. 17 и 18 показывают, что мо54


жет произойти при использовании чистки спектров в тех случаях, когда условия (118) и (119) нарушаются. Действительно, в первом случае (см. рис. 17) из двух заданных гармоник с частотами 1 и 2 Гц при чистке теряется гармоника с частотой 1 Гц, а во втором случае (см. рис. 18) из двух заданных гармоник с частотами 1.3 и 3.3 Гц при чистке возникает гармоника на частоте 2.3 Гц, а истинные линии исчезли. На рис. 19 мы показываем результаты анализа кривой блеска Плутона, полученной по наблюдениям на нормальном астрографе Пулковской обсерватории в 1954-1990 гг. Методика обработки наблюдательного материала и первый вариант его спектрального анализа описаны в работе Бронниковой и др. (1992). Здесь получена оценка векового падения блеска и показано, что фотометрическая кривая Плутона содержит компонент с периодом 7.8 лет. Эти результаты были получены на основе качественной интерпретации периодограммы без привлечения процедур чистки и без оценки статистической достоверности разделения сигнального и шумового компонентов. В нашем анализе с вероятностью 0.95 подтверждается существование периодического компонента с периодом 7.8 лет и показывается, что гармоники с периодами 2.71 и 1.59 года нельзя считать достоверными.

55


56


щью модифицированного метода CLEAN. a Исходный временной ряд; b модуль спектрального окна; c "грязный"спектр; d "чистый"спектр; e коррелограмма; f восстановленный на равномерной сетке временной ряд. Этот пример показывет успешную работу метода. Исходный временной ряд типа сигнал+шум задан на сетке с периодическими пропусками (n = 5, p = 5, m = 15) с шагом 0.1 c. Собственные частоты спектрального окна равны 1 = 1 и 1 = 3 Гц. Сигнал представлен суммой ? ? двух косинусоид с равными амплитудами, нулевыми фазами и частотами 1 = 1.2 и 2 = 3.1 Гц. Восстановление истинного спектра в данном случае гарантируется выполнением условий (118) и (119). Шум моделирован псевдослучайными числами, распределенными по нормальному закону. Отношение сигнал/шум равно единице. Разделение сигнала и шума производится с вероятностью 0.99. Фильтрация восстановленного ряда не применялась.

Рис. 15: Спектральный анализ неравномерного временного ряда с помо-

57


58


Рис. 16: Спектральный анализ неравномерного временного ряда с помо-

щью модифицированного метода CLEAN. a Исходный ряд; b модуль спектрального окна; c "грязный"спектр; d "чистый"спектр; e коррелограмма; f восстановленный на равномерной сетке временной ряд. Этот пример вновь демонстрирует усрешную работу метода CLEAN. Условия (118) и (119) выполняются. Исходный временной ряд типа сигнал+шум задан на сетке с длинным пропуском (n = 40, p = 40, n = 40) с шагом 0.1 c. Собственные частоты наиболее интенсивных пиков спектрального окна равны 1 = 0.125 и 2 = 0.250 Гц. Сигнал представлен ? ? суммой двух косинусоид с амплитудами, равными 1 и 0.7, нулевыми фазами и частотами 0.3 и 3.1 Гц. Шум моделирован псевдослучайными числами, распределенными по нормальному закону. Отношение сигнал/шум равно единице. Разделение сигнала и шума производится с вероятностью 0.99. Демонстрируется возможность узкополосной фильтрации: при образовании "чистого"спектра спектральная линия на частоте 3.1 Гц была уничтожена. Коррелограмма (e) и восстановленный временной ряд на равномерной сетке (f ) соответствуют отфильтрованному сигналу.

59


60


Рис. 17: Спектральный анализ неравномерного временного ряда с помо-

щью модифицированного метода CLEAN. a Исходный временной ряд; b модуль спектрального окна; c "грязный"спектр; d "чистый"спектр; e коррелограмма; f восстановленный на равномерной сетке временной ряд. В этом примере показана ситуация, когда чистка "грязного"спектра неравномерного ряда производится с частичным успехом. Исходный временной ряд задан на сетке с периодическими пропусками (n=3, p=7, m=15) с шагом 0.1 c. Собственные частоты наиболее интенсивных боковых пиков спектрального окна равны 1 = 1 и 2 = 2 Гц. Сигнал ? ? представлен суммой двух косинусоид с равными амплитудами, нулевыми фазами и частотами 1.0 и 2.0 Гц. Шум отсутствует. В этом случае частоты исследуемого ряда совпадают с собственными частотами спектрального окна, и, следовательно, условия восстановления спектра (118) и (119) нарушаются. Этот факт объясняет то, что в "чистом"спектре появилась лишь линия на частоте 2 Гц, а линия на частоте 1 Гц исчезла. Коррелограмма (e) и восстановленный временной ряд (f ) соответствуют "слишком чистому"спектру (d).

61


62


Рис. 18: Спектральный анализ неравномерного временного ряда с помо-

щью модифицированного метода CLEAN. a Исходный ряд; b модуль спектрального окна; c "грязный"спектр; d "чистый"спектр: e коррелограмма; f восстановленный на равномерной сетке временной ряд. Этот пример демонстрирует случай, когда результаты чистки методом CLEAN оказываются неверными. Исходный временной ряд задан на сетке с периодическими пропусками (n=3, p=7, m=15) с шагом 0.1 c. Собственные частоты наиболее интенсивных линий спектрального окна равны 1 и 2 Гц. Модельный ряд представлен суммой двух косинусоид с равными амплитудами, нулевыми фазами и частотами 1.3 и 3.3 Гц. Шум отсутствует. Условия восстановления истинного спектра (118) и (119) нарушены. Это приводит к тому, что в "грязном"спектре самая интенсивная линия на частоте 2.3 Гц оказывается ложной. Она же оказалась единственной линией и в "чистом"спектре (d). Коррелограмма (e) и восстановленный временной ряд (f ) соответствуют "чистому"спектру (d) с ложной гармоникой.

63


64


Рис. 19: Спектральный анализ кривой блеска Плутона модифицирован-

ным методом CLEAN. a Исходный временной ряд 100(m - 15), где m звездная величина. Фотометрические измерения даны с шагом 1 год на промежутке 19541990 гг.; b модуль спектрального окна; c "грязный"спектр; d "чистый"спектр; e коррелограмма; f компонент кривой блеска с периодом 7.8 лет.

65


Упражнения
1. Составьте программу анализа временных рядов методом CLEAN, пользуясь алгоритмом, приведенным на с. 5053. В качестве модельного ряда используйте функцию
xk = A1 cos(2 1 tk + 1 ) + A2 cos(2 2 tk + 2 ) + n nk ,

(162)

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

2. Задайте параметры модельного ряда (без шумового компонента) и

убедитесь в том, что при переходе к нерегулярным наблюдениям в периодограмме модельного ряда появляются "духи"(ложные линии). 2. Задайте параметры модельного ряда с шумом (отношение "сигнал к шуму" = 1) и убедитесь в том, что при переходе к зашумленным нерегулярным наблюдениям в периодограмме модельного ряда появляются "духи"двух типов: Духи шума и Духи иррегулярности. 3. Постройте модельный ряд (162) (без шумового компонента) с параметрами A1 = 1, 1 = 0.1 Гц, 1 = 0; A2 = 1, 2 = 0.2 Гц, 2 = 0 на сетке отсчетов с циклическими пропусками: n = 5, p = 5, m = 15, tk = t k , k = 0, 1, ..., N - 1, при N = 150, t = 1 c. Проведите чистку периодограммы. Почему в чистой периодограмме потерялась одна из заданных гармоник? 4. Получите периодограмму модельного ряда (162) (без шумового компонента) с параметрами A1 = 1, 1 = 0.01 Гц, 1 = 0; A2 = 1, 2 = 0.02 Гц, 2 = 0 на сетке временных отсчетов с длинным пропуском: n1 = 40, p = 60, n2 = 40, tk = t k, k = 0, 1, ..., N - 1, при N = 140, t = 1 c. Проведите чистку периодограммы. Почему результаты чистки оказываются неверными? 5. Модифицируйте программу метода CLEAN так, чтобы она давала возможность получения амплитуд и фаз гармоник "чистого спектра"(см. формулы (143) и (144)).

66


Литература
Андерсон Т., 1976. Статистический анализ временных рядов. М.: Мир. Барнинг, 1963. Barning F.J.M. The numerical analysis of the lightcurve of 12 Lacertae // Bull. Astr. Inst. Neth. Vol. 17, N 1. P. 22-28. Бендат Дж., Пирсол А., 1989. Прикладной анализ случайных данных. М.: Мир. Бронникова и др., 1992. Бронникова Н.М., Бершакова С.М., Прудникова Е.П. Исследование блеска Плутона // Астрофотография в исследовании Вселенной. СПб., С. 226-247. Блэкман и Тьюки, 1959. Blackman R.B., Tukey J.W. The Measurement of Power Spectra. NY.: Dower. Витязев, 1994. Vityazev V.V. The Time Interferometer: a new tool of spectral analysis of time series // Astron. and Astrophys. Trans. Vol. 5. P. 177-210. Витязев В.В., Прудникова Е.Я., 1994. Спектры скважности астрометрических рядов наблюдений // Вестн. С.-Петерб. ун-та. Сер. 1. Вып. 2, N 8. С. 78-86. Витязев, 1995. Vityazev V.V. Time Series Analysis of Unequally Spaced data: The statistical properties of the Schuster periodogram // Astron. and Astrophys. Trans. Vol. 11. P. 159-173. Диминг, 1975. Deeming T.J. Fourier analysis with unequally-spaced data // Astrophys. and Space Sci. Vol. 36. P. 137-158. Карбонел и др., 1992. Carbonell M., Oliver R., Ballester J.L. Power spectra of gapped time series: a comparison of several methods // Astron. and Astrophys. Vol. 264. P. 350-360. Ломб, 1976. Lomb N.R., Least-squares frequency analysis of unequally spaced data // Astrophys. Sp. Sci. Vol. 39. P. 447-462. Отнес и Эноксон, 1982. Прикладной анализ временных рядов. М.: Мир. Робертс и др., 1987. Roberts D.H., Lehar J., Dreher J.W. Time Series Analysis with CLEAN. I. Derivation of a Spectrum // Astrophys. J. N 4. P. 968-989. Скаргл, 1982. Scargle J.D. Studies in astronomical time series analysis. 2. Statistical aspects of spectral analysis of unevenly spaced data // Ap. J. Vol. 263. P. 835-853. Скаргл, 1989. Scargle J.D. Studies in astronomical time series analysis. 3. Fourier transforms. Autocorrelation function and cross-correlation functions of unevenly spaced data // Ap. J. Vol. 343. P. 874-887.

67


Теребиж В.Ю., 1992. Анализ временных рядов в астрофизике. М.: Наука. Томпсон и др., 1986. Tompson A.R., Moran J.M., Swenson Jr. G.W. Interferometry and Synthesis in Radio Astronomy // NY.: John Wiley & Sons. Уолкер, 1914. Walker G. // Calcutta Ind. Met. Mems. Vol. 21. N 9. Хогбом, 1974. Hogbom J.A. // Astron. and Astrophys. Suppl. Vol. 15. N 417. Хорн и Балиунас, 1986. Horne J.H., Baliunas S.L. A prescription for period analysis of unevenly spaced time series // Ap. J. Vol. 302. P.757-763.

68

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

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

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

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