CC BY
117
АСТРОНОМИЯ
УДК 520.87, 519.25
О ПОИСКЕ ПЕРИОДИЧЕСКИХ КОМПОНЕНТ В НАБЛЮДАТЕЛЬНЫХ ДАННЫХ*
Р. В. Балуев
С.-Петербургский государственный университет, аспирант, roman@astro.spbu.ru
Введение
В ходе анализа астрономических наблюдений часто возникает задача выявления в имеющемся временном ряде периодичностей. Исходные наблюдательные данные могут иметь разную природу, а физические источники возможных периодичностей могут быть весьма различными. Например, нашей целью может быть поиск переменных звезд на основе фотометрических массивов данных или поиск внесолнечных планет на основе высокоточных измерений лучевых скоростей. Между тем, математические алгоритмы анализа этих данных и выявления в них периодических компонент остаются в целом одинаковыми. Данная статья представляет собой краткий обзор определенного семейства популярных алгоритмов поиска периодичностей, основанных на периодограмме Ломба—Скаргла [1, 2]. Значительное внимание будет уделено вопросу оценки статистической значимости находимых периодичностей. Поскольку любые измерения отягощены случайными погрешностями, всегда остается некоторая ненулевая вероятность того, что найденная периодичность могла быть вызвана случайными ошибками данных. Если эта вероятность достаточно мала (скажем, меньше 1%), то мы можем считать найденную периодичность статистически значимой, и нет — в противном случае. Задача практической оценки статистической значимости найденных периодичностей не является тривиальной и не была решена с достаточной точностью на протяжении трех десятилетий. В дополнение к работам [1, 2] здесь можно упомянуть, в частности, статьи [3-8]. Эта задача требует довольно серьезного математического аппарата теории вероятностей и математической статистики, включая теорию экстремальных значений случайных процессов. В настоящей статье описываются недавние практически важные результаты, полученные автором в этой области.
* Работа выполнена при финансовой поддержке РФФИ (грант №06-02-16795) и Совета по грантам Президента РФ для государственной поддержки молодых российских ученых и ведущих научных школ (грант №НШ-1323.2008.2).
© Р.В.Балуев, 2009
Периодограмма Ломба—Скаргла и ее обобщения
Пусть мы имеем временной ряд, состоящий из N однородных наблюдений. Каждое наблюдение характеризуется тройкой чисел (ti, xi, ai), где ti есть момент i-го наблюдения, xi —его количественный результат, а ai — стандартное отклонение его ошибки. На протяжении всей статьи мы будем предполагать, что ошибки измерений статистически независимы и нормально распределены с нулевым средним и указанной дисперсией.
Сформулируем нашу задачу в терминах теории проверки статистических гипотез и статистических критериев. Пусть базовая гипотеза H, которую мы собираемся проверить, состоит в том, что наблюдаемая физическая величина x постоянна и равна нулю. Альтернативная гипотеза K состоит в том, что эта физическая величина меняется синусоидальным образом в соответствии с моделью x = p(t,a,b) = a cos wt + b sin wt = Acos(wt + а). Можно заметить, что если бы круговая частота w = 2nf была заранее известна, то мы бы имели частный случай классической линейной регрессионной задачи математической статистики. Оптимальный критерий для проверки соответствующих гипотез хорошо известен и основан на линейном методе наименьших квадратов [9, глава 7]. Введем функцию
которая определяет, насколько хорошо наша синусоидальная модель описывает реальные измерения (за разъясниниями используемых в статье математических обозначений, таких как треугольные скобки (•}, см. Приложение). Определение (1) можно привести к более удобной классической форме периодограммы Ломба—Скаргла [1, 2]:
периодического сигнала заданной частоты. Построив функцию г(/) в некотором заранее заданном диапозоне частот, мы можем найти наибольший пик ее графика. Этот пик соответствует наиболее правдоподобной частоте искомого сигнала, а его высота характеризует статистическую значимость последнего. Большие значения тах г(/) заставляют нас отвергнуть гипотезу Н (нет сигнала) в пользу альтернативы К (есть сигнал).
Определение периодограммы Ломба—Скаргла (1)-(2) можно обобщать в разных направлениях. Прежде всего, модель сигнала не обязательно должна быть синусоидальной, а может даваться в общем виде как ^(£, в, /), где вектор в содержит ! неизвестных параметров сигнала. Моделируемый сигнал отсутствует при в = 0: ^1в=о = 0. Кроме того, проверяемая гипотеза Н также может содержать некоторую модель измеряемой физической величины ^н(£, вн), где вектор вн содержит заранее неизвестных параметров базовой модели. Таким образом, альтернативная гипотеза К описывает измеряемую величину моделью ^к(£, вк, /) = Мн(^, вн) + в, /), где вектор вк = {вн, в} содержит все + ! параметров, подлежащих определению (помимо частоты /).
Мы хотим проверить, согласуется ли с наблюдениями основная гипотеза Н : в = 0 (нет сигнала), или ее следует отвергнуть в пользу альтернативы К : в = 0 (есть периодический сигнал).
Как правило, модель представляет собой просто константу или полиномиальный тренд. Здесь мы не должны забывать, что на практике часто бывает недостаточно
(1)
Таким образом, чем больше г(/) на данной частоте /, тем более правдоподобно наличие
простого предварительного центрирования временного ряда (т. е., простого вычитания из него среднего значения или наилучшей модели полиномиального тренда). Как показано в работе [10], при неравномерном распределении наблюдений во времени, синусы и косинусы в модели сигнала оказываются неортогональными константе и степеням времени (в смысле скалярного произведения, задаваемого операцией (Ф1Ф2}). При этом, если мы используем простую предварительную центрировку временного ряда, чувствительность периодограммы может значительно снижаться [7]. Для получения более эффективной периодограммы необходимо ввести следующее обобщение определения (1):
г(/)= [хН - Хк(/)2] /2, Хн = тіп((ж - мн)2) , Хк(/) = т1п((х - Мк)2) . (3)
ви ви,в
Рассмотрим более подробно случай, когда модели и ц линейны по неизвестным параметрам: ^н(4, вн) = вн • (^) и ^(£, в, /) = в • ^(£, /), где векторы и ^ содер-
жат и ! фиксированных базисных функций времени, так что при фиксированной частоте / мы имеем обычную задачу линейной регрессии. В этом случае статистический критерий, основанный на статистике г(/), является равномерно наиболее мощным инвариантным критерием [9, глава 7]. Иными словами, чувствительность значений г(/) к слабым периодическим сигналам (следующим заданной модели) максимальна среди всех возможных критериев, инвариантных относительно перенумерации наблюдений. Данный вид линейных периодограмм, построенных по методу наименьших квадратов, рассматривался в довольно общем виде (правда, без учета возможности непустой базовой модели у«н) в работах [5, 6, 11]. Линейность моделей по параметрам позволяет вычислить минимумы в формулах (3) в явном виде, а также придать таким периодограммам элегантную трактовку в терминах проекций массива наблюдений на линейные подпространства, задаваемые гипотезами Н и К. Чтобы выписать явное выражение для г(/) в общем линейном случае, введем следующие информационные матрицы Фишера, связанные с параметрами вн и в:
Оее = (<? ® <?), Овпв = (<Рн ® <?), Оввн = ®впв, Овпвп = (Фн ® <Рн>- (4)
Кроме того, сформируем полную информационную матрицу Овквк:
Система общих формул (4)—(8) выглядит несколько громоздко, однако ее можно существенно упростить путем ортонормировки базиса [5, 6]. После такой ортонормировки мы получим, в частности, Оввн = Овнв = 0, О = Овв = I и д = дв.
Определение периодограммы (3) требует, чтобы дисперсии ошибок <г* были известны точно. Такая ситуация на практике встречается редко. В лучшем случае известны лишь
(5)
и определим следующие вектора, зависящие от наблюдений ж*:
Я ви = (xlPн), Я вк = (х^к} Я в = (ж^}-
(6)
Тогда легко проверить, что
*(/) = ЯТ0-1Я, (7)
где
О = Овв — Овв
Овив, Я = Яв — Овви ОвНви Яви •
(8)
веса наблюдений ж 1/^2, так что сами дисперсии известны с точностью до постоянного множителя. На практике из этого затруднения выходят, нормализуя периодограмму (3) на выборочную дисперсию невязок наилучшей (в смысле метода минимума X2) модели, соответствующей гипотезе Н или К. Мы введем следующие модификации базисной периодограммы (3):
хн-хНЛ т дг хн-хНЛ т х2н
' 2Д ’ = ^ 2ХНЛ •
Периодограммы ^х(/) и г2(/) представляют собой те самые нормализации г(/), а периодограмма гэ(/) связана с отношением правдоподобия. Все три периодограммы (9) эквивалентны, так как они представляют собой однозначные и монотонные функции друг друга:
— = 1 - е~2гз/Мк — = е2^ _ 1 (1 ~ 2гЛ (\ + —1] = 1 (ю)
Вероятность ложной тревоги
Описанные периодограммы имели бы малую ценность без метода оценки статистической значимости их пиков. Поскольку наблюдения отягощены случайными ошибками, мы всегда находимся под угрозой совершения ошибок двух типов: ошибочное отвержение основной гипотезы Н (т. е. ошибочный вывод о существовании сигнала, «ложная тревога») и ошибочное неотвержение Н (т.е. невыявление реальной периодичности). Как правило, ложные тревоги считаются более опасными, поэтому критерием статистической значимости обычно выступает «вероятность ложной тревоги», то есть вероятность отвергнуть гипотезу Н, когда в действительности она верна. Мы будем обозначать эту вероятность как РАР. Задавшись некоторым малым критическим значением РАР* (скажем, 1%), мы можем заявить что найденная периодичность статистически значима (когда РАР < РАР*) или незначима (когда РАР > РАР*).
Чтобы определить эту вероятность, нам необходимо найти функцию распределения той величины, которая является основой для используемого критерия (при этом считая гипотезу Н верной). В случае, когда частота сигнала заранее известна, критерием значимости выступает значение периодограммы ,г(/) на данной частоте /. В случае, когда частота заранее не известна, таким критерием выступает максимальный отсчет периодограммы в заданном частотном диапозоне. Обычно периодограммы строятся в диапозоне частот от нуля до некоторого предельного значения /тах. Поэтому нас будет интересовать распределение величины шахо^/^/тах г(/).
В случае заранее известной частоты сигнала мы можем использовать тот факт, что распределение одиночных отсчетов 2г(/) есть х2 распределение с й степенями свободы (см. [9], §7.1). В случае периодограммы Ломба—Скаргла (й =2 параметра синусоидального сигнала) мы получаем широко известное экспоненциальное выражение РАРН1Пв1е = в-г. Распределения значений модифицированных периодограмм (9) суть ^-распределение с й и Nк степенями свободы для 2^/^ и В-распределение с теми же числами степеней свободы для 2^1/№н. Функцию распределения периодограммы ^3 можно получить из функций распределения для ^1,2 с использованием соотношений (10). Для случая й =2 значения периодограммы ^э распределены в точности так же, как и значения периодограммы г (т. е. экспоненционально).
Для распределения максимального отсчета даже простейшей периодограммы Ломба—Скаргла не существует точного выражения. На практике довольно часто применяется формула «повторных испытаний»
РАРтах = 1 — (1 — РА-Рв1п§1е) ^1П(^ « NindРAPsing1e, (11)
которая была бы верной, если бы максимальный отсчет периодограммы вычислялся среди конечного дискретного набора из N1^ статистически независимых отсчетов. К сожалению, такая система «независимых частот» существует только в случае равномерного распределения наблюдений во времени. Для неравномерных временных рядов такую систему построить вряд ли возможно. Например, автору работы [4] не удалось найти в неравномерном случае более двух независимых отсчетов периодограммы. Кроме того, даже в случае равномерного временного ряда «независимые частоты» слишком разрежены, так что существенные детали периодограммы проваливаются в промежутки между ними. Поэтому на практике, как правило, периодограммы строятся на значительно более плотной сетке, которая практически эквивалентна непрерывному отрезку частот. В этом случае понятие «число независимых частот» теряет физический смысл, а формула (11) становится лишь эвристической аппроксимацией [8]. При этом мы уже не имеем сколько-нибудь приемлемой аналитической оценки для параметра N1^ и, следуя авторам работы [3], вынуждены оценивать этот параметр методом Монте-Карло, используя (11) как априорную модель для получаемой эмпирической функции распределения максимального отсчета периодограммы. При этом у нас нет никаких теоретических гарантий что формула (11) вообще верна хотя бы асимптотически, в каких-то предельных условиях.
В работе [12] предложен другой метод оценки искомого распределения, основанный на результатах теории экстремальных значений случайных процессов. Основная идея этого подхода — рассматривать периодограмму наблюдений, отягощенных случайными ошибками, как случайный процесс. При этом для оценки распределения его максимумов применяется так называемый метод Райса [13] (см. также ссылки в работе [12]). Для периодограмм Ломба—Скаргла, например, получена оценка
РАртах(>, /тах) \Уе~*у/г, ^ = /тах^еЯ, (12)
где Тед есть эффективная длина временного ряда (пропорциональная взвешенной дисперсии моментов наблюдений и обычно близкая к его фактической длине). Для линейных периодограмм общего вида имеем
г Ч (Й-1)/2
рАРтах(^ /шах) £ А(иах)~ {-) , (13)
где функция А(/тах), определяющая зависимость от ширины частотной полосы поиска, зависит от регрессионного базиса , от распределения моментов наблюдений и их весов. Точное выражение для А(/тах) довольно громоздко, и мы не будем его здесь приводить. Однако при этом стоит отметить, что с достаточной точностью А(/тах), как правило, приблизительно пропорциональна /тах. Формулы, аналогичные (12)—(13), но для модифицированных периодограмм 21,2,э(/), можно найти в [12].
Тот факт, что формулы (12)—(13) ограничивают вероятность ложной тревоги сверху, весьма важен. Как показано в [12] на примере периодограммы Ломба—Скаргла, эти формулы весьма точны для случая равномерного временного ряда, содержащего большое количество наблюдений. Неравномерное распределение моментов наблюдений
и/или малое их число могут вызывать некоторые отклонения эмпирических функций распределения от теоретических, но при этом теоретические выражения (12)—(13) всегда недооценивают статистическую значимость. Это означает, что возможные ошибки данных формул не должны благоприятствовать «ложным тревогам». При этом для больших г (меньших РАР) ошибки вычисления вероятности ложной тревоги существенно уменьшаются. Это тоже важно, так как мы наиболее заинтересованы в точном вычислении именно малых значений РАР, обычно в диапозоне 0.001-0.01. Кроме того, выражения (12)—(13)—замкнутые и аналитические, и не требуют никакого предварительного численного моделирования.
Общее выражение (13) может быть использовано для различных обобщений периодограммы Ломба—Скаргла. Например, в работе [11] рассматривались мультигармони-ческие периодограммы, для которых в качестве модели сигнала выступал отрезок ряда Фурье заданной степени п. Мультигармонические периодограммы более эффективны в тех случаях, когда предполагаемый периодический сигнал может иметь сильно несинусоидальную форму, как, например, при поиске переменных звезд. В этом случае базис ^ включает в себя п первых гармоник функций из базиса Фурье. Данная ситуация рассмотрена в [14], где показано, что для мультигармонической периодограммы г(/)
Аналогичные формулы получены и для модифицированных мультигармонических периодограмм £1,2,з(/). Заметим, что «1 = 1, «2 = 14/9 « 1.556, «з « 1.062, «5 « 0.136.
Результаты численного моделирования
В работах [12, 14] аналитические оценки (12), (14), а также их аналоги для соответствующих модифицированных периодограмм, проверены при помощи численного Монте-Карло моделирования при различных входных условиях. Полученные результаты можно свести к следующим позициям.
1. Указанные аналитические оценки являются весьма точными, если в заданном частотном диапозоне эффекты неравномерного распределения моментов наблюдений (например, элайзинг) малы.
2. Если в рассматриваемом диапозоне частот проявляются существенные эффекты неравномерности наблюдений (например, явление элайзинга в случае периодических пропусков наблюдений или «шумоподобная» утечка мощности из главного пика в случае малого числа наблюдений), то точность аналитических оценок вероятности ложной тревоги ухудшается. Но даже при сильно неравномерном распределении точек это ухудшение точности сравнительно невелико для практически значимых малых уровней РАР.
3. Указанные аналитические оценки в любых условиях действительно ограничивают вероятность ложной тревоги сверху.
Таким образом, численное моделирование показывает в целом высокую эффективность оценок (12), (14).
Заключение
В работе приведен краткий обзор методов поиска (выявления) периодических компонент в наблюдательных данных. Введен довольно общий класс периодограмм, построенных на основе метода наименьших квадратов. Значительное внимание уделено задаче оценки статистической значимости пиков периодограмм. Случаи периодограммы Ломба—Скаргла и мультигармонических периодограмм рассмотрены более подробно. Для этих частных случаев описаны последние теоретические результаты исследования их распределений, имеющие важнейшее значение в задаче оценки статистической значимости периодичностей, выявленных при помощи указанных периодограмм. В качестве конкретного астрономического приложения описанной здесь методологии можно упомянуть работу [15], посвященную анализу временного ряда лучевых скоростей звезды НБ37124, обладающей планетной системой из (не менее) трех планет-гигантов.
Приложение: Некоторые обозначения
Введем следующую осредняющую операцию:
Здесь ф(£) произвольная функция, а tj и аг представляет собой моменты и стандартные ошибки наблюдений. Аналогичные определения можно ввести для дискретных последовательностей фг, определенных для N узлов tj. Для простоты формул аргумент t и индекс i в основном тексте обычно опускаются. Заметим, что операцию (^>i(t^2(t)) можно также рассматривать как скалярное произведение функций в гильбертовом пространстве.
Пусть x и у —произвольные многомерные вектора. Тогда их диадное произведение x <g> у = xyT представляет собой матрицу, составленную из попарных произведений xiyj.
I есть единичная матрица.
Литература
1. Lomb N. R. Least-squares frequency analysis of unequally spaced data // Astrophys. & Space Sci. 1976. Vol. 39. P. 447-462.
2. Scargle J. D. Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data // Astrophys. J. 1982. Vol. 263. P. 835-853.
3. Horne J. H., Baliunas S. L. A prescription for period analysis of unevenly spaced time series // Astrophys. J. 1986. Vol. 302. P. 757-763.
4. Koen C. Significance testing of periodogram ordinates // Astrophys. J. 1990. Vol. 348. P. 700702.
5. Schwarzenberg-Czerny A. The distribution of empirical periodograms: Lomb-Scargle and PDM specta // Mon. Not. Roy. Astron. Soc. 1998. Vol. 301. P. 831-840.
6. Schwarzenberg-Czerny A. Period search in large datasets // Baltic Astron. 1998. Vol. 7. P. 43-
7. Cumming A., Marcy G. W., Butler R. P. The Lick planet search: detectability and mass thresholds // Astrophys. J. 1999. Vol. 526. P. 890-915.
8. Frescura F.A.M., Engelbrecht C.A., Frank B.S. Significance of periodogram peaks and a pulsation mode analysis of the Beta Cephei star V403Car // Mon. Not. Roy. Astron. Soc. 2008. Vol. 388. P. 1693-1707.
N
(15)
69.
9. Леман Э. Л. Проверка статистических гипотез. М.: Наука, 1979.
10. Ferraz-Mello S. Estimation of periods from unequally spaced observations // Astron. J. 1981. Vol. 86. P. 619-624.
11. Schwarzenberg-Czerny A. Fast and statistically optimal period search in uneven sampled observations // Astrophys. J. 1996. Vol. 460. P. L831-840.
12. Baluev R. V. Assessing the statistical significance of periodogram peaks // Mon. Not. Roy. Astron. Soc. 2008. Vol. 385. P. 1279-1285.
13. Azai's J.-M., Wschebor M. The distribution of the maximum of a Gaussian process: Rice method revisited / Sidoravicius V., ed. // In and Out of Equilibrium: Probability with a Physics Flavor. Vol. 51 of Prog. Prob. Birkhauser book. Boston, 2002. P. 321-348.
14. Baluev R. V. Detecting non-sinusoidal periodicities in observational data using multiharmonic periodograms // Mon. Not. Roy. Astron. Soc., in preparation.
15. Baluev R. V. Resonances of low orders in the planetary system of HD37124 // Celest. Mech. Dynam. Astron., accepted.
Статья поступила в редакцию 3 октября 2008 г.
