Содержание к диссертации
Введение
2 Поиск периодичностей в наблюдательных данных 13
2.1 Введение и определения 13
2.2 Вероятность ложной тревоги 16
2.3 Численная проверка методом Монте-Карло 20
2.4 Поиск экзопланет 25
3 Астрофизическое и инструментальное дрожание лучевой скорости 29
3.1 Традиционный подход 29
3.2 Переход к методу максимального правдоподобия 32
3.3 Смещение ММП-оценки дрожания лучевой скорости 34
3.4 Влияние негауссовых ошибок 36
3.5 Статистическое сравнение орбитальных моделей 39
3.6 Периодограммы отношения правдоподобия 43
3.7 Численная максимизация функции правдоподобия 44
3.8 Периодические систематические ошибки 46
4 Эффекты нелинейности моделей 53
4.1 Введение 53
4.2 Смещение оценок орбитальных параметров 53
4.3 Оценка надежности нелинейных оценок 56
5 Планетная система у звезды HD37124 59
5.1 Введение 59
5.2 Доступные данные 60
5.3 О взаимных возмущениях в системе 62
5.4 Оценка орбитальных параметров «в лоб» 64
5.5 Учет требования динамической устойчивости 69
5.6 Значение апсидальных коротаций 74
5.7 Проверка существования других планет 81
5.8 Выводы 82
6 Оптимальное планирование наблюдений 84
6.1 Предпосылки 84
6.2 Критерии оптимальности 85
6.3 Уточнение оценок 89
6.4 Дискриминация моделей 92
6.5 Планирование нескольких наблюдений 93
6.6 Область применимости 94
6.7 Практическая эффективность 96
Заключение 103
Приложения 105
- Вероятность ложной тревоги
- Переход к методу максимального правдоподобия
- Смещение оценок орбитальных параметров
- Учет требования динамической устойчивости
Введение к работе
Актуальность исследования 4
Цели работы 5
Научная и практическая ценность результатов 7
Новизна и достоверность результатов 8
Результаты, выносимые на защиту 9
Апробация работы 9
Структура диссертации 9
Публикации по результатам работы 11
Некоторые математические обозначения 11
2 Поиск периодичностей в наблюдательных данных 13
Введение и определения 13
Вероятность ложной тревоги 16
Численная проверка методом Монте-Карло 20
Поиск экзопланет 25
3 Астрофизическое и инструментальное дрожание лучевой ско
рости 29
Традиционный подход 29
Переход к методу максимального правдоподобия 32
Смещение ММП-оценки дрожания лучевой скорости 34
Влияние негауссовых ошибок 36
Статистическое сравнение орбитальных моделей 39
Периодограммы отношения правдоподобия 43
Численная максимизация функции правдоподобия 44
Периодические систематические ошибки 46
4 Эффекты нелинейности моделей 53
Введение 53
Смещение оценок орбитальных параметров 53
Оценка надежности нелинейных оценок 56
5 Планетная система у звезды HD37124 59
Введение 59
Доступные данные 60
О взаимных возмущениях в системе 62
Оценка орбитальных параметров «в лоб» 64
Учет требования динамической устойчивости 69
Значение апсидальных коротаций 74
Проверка существования других планет 81
Выводы 82
6 Оптимальное планирование наблюдений 84
Предпосылки 84
Критерии оптимальности 85
Уточнение оценок 89
Дискриминация моделей 92
Планирование нескольких наблюдений 93
Область применимости 94
Практическая эффективность 96
7 Заключение 103
Приложения 105
А Метод Раиса и ЛМНК-периодограммы 105
АЛ Идея метода 105
А.2 Применение к ЛМНК-периодограммам 107
8 Асимптотические свойства ММП-оценок при нестандартных
условиях 117
С Нахождение коротационных орбитальных конфигураций 121
D Уравнения в вариациях (уравнения чувствительности) 125
Е Некоторые матричные тождества 127
Литература 129
Список иллюстраций 137
Список таблиц 138
1. Введение
1.1. Актуальность исследования
Актуальность данной работы во многом определяется чрезвычайной важностью исследований внесолнечных планегных систем. Вопрос о существовании планет у других звезд солнечного типа волновал умы человечества на протяжении по крайней мере четырех столетий. Открытие первой такой внесолнечной планеты (экзопланеты) в 1995 году [75] поставило перед астрономическим научным сообществом множество новых вопросов и породило новую область научных исследований. Эти вопросы имеют как самостоятельное научное значение, так и входят составной частью в другие научные задачи. Одной из этих задач является, например, проблема распространенности жизни во Вселенной.
Со времени открытия первой внесолнечной планеты, обращающейся вокруг звезды солнечного типа, прошло 14 лет. В течение этого промежутка времени число известных экзопланет непрерывно росло и уже превысило 300. К настоящему времени уже известно около 30 планетных систем, содержащих две или более планеты ([32, 8, 9,15,1], см. также профессиональный интернет-ресурс The Extrasolar Planets Encyclopaedia, ). До сих пор большая часть экзопланет открывается методом лучевых скоростей. Этот метод поиска экзопланет оказался весьма эффективным на практике. Сейчас он дает новую планету в среднем каждую неделю. Суть данного метода состоит в том, что планета обнаруживается по колебанию лучевой скорости своей звезды, наблюдаемому с Земли. Лучевая скорость звезды, вокруг которой обращается невидимый (планетный) спутник, дается следующим выраженем (см., например, [48, 15, 16]):
v = К [cos (д + v) + е cos д] + с. (1.1)
Здесь К — амплитуда кривой лучевой скорости, д — аргумент перицентра орбиты планеты, е — эксцентриситет, v — истинная аномалия, с — лучевая скорость барицентра системы. Истинная аномалия зависит от времени, эксцентриситета е, средней долготы Л в заданный фиксированный момент времени, а также от орбитального периода Р.
В выбранной параметризации величины К, Р, д, е, Л представляют собой первичные параметры, определяемые в ходе подгонки модельной кривой лучевой скорости (1.1) к имеющимся наблюдениям. Существуют также вторичные параметры — величины, определяемые как функции от первичных параметров. Ко вторичным параметрам относятся, в частности, большая полуось а ее орбиты и минимальная масса планеты msinz, где г — угол наклона плоскости орбиты к картинной плоскости. Приближенно их можно вычислить как
msini ~ К (~У3 = МКР^М?3, а с {^^-\" = ЛР2"МІ'\ (1.2)
Здесь G — гравитационная постоянная, М* — масса звезды, а К = К л/1 — е2 — модифицированная амплитуда. Если период Р измеряется в сутках, М* — в массах Солнца М, а К — в метрах в секунду, то постоянные множители Л4 « 4.919 10-3[М3ирМ^2/3м-1сут~^3сек] и Л « 1.957 10-2[Мд1/3а.е.сут-2/3]. Приближения (1.2) действительны при т <С М*. Как известно, информацию об орбитальном наклоне г (а значит и об истинной массе планеты т) нельзя получить из наблюдений одних лучевых скоростей, если только взаимные возмущения планет в системе не удается регистрировать напрямую в кривой лучевой скорости.
Сложившаяся ситуация «доминирования» метода лучевых скоростей делает улучшение уже имеющихся и разработку новых, более эффективных, алгоритмов статистического анализа и планирования наблюдений лучевых скоростей звезд в программах поиска экзопланет весьма актуальной задачей.
1.2. Цели работы
Несмотря на неоспоримые успехи программ поиска внесолнечных планет, использующих метод лучевых скоростей, до сих пор большинство открытых экзопланет официально именуются всего лишь «кандидатами» в экзопланеты. Эта осторожность связана с косвенной природой данного метода. Например, на практике иногда оказывется, что наблюдаемая переменность лучевой скорости звезды вызвана не наличием у этой звезды планеты, а какими-либо явлениями активности в звездной атмосфере, например, пятнами (см. [67]). Но параметры планет часто определяются ненадежно просто с точки зрения теории статистической обработки наблюдений. К сожалению, в погоне за быстрыми и массовыми открытиями, наблюдатели часто забывают уделять должное внимание контролю качества не только применяемых методов обработки и планирования наблюдений, но даже и результатов самих этих наблюдений. В конечном счете может оказаться, что рассматриваемые программы наблюдений в существенной части работают «вхолостую». Главная цель данной работы — представить подробный анализ качества измерений лучевых скоростей, получаемых в современных программах поиска экзопланет, алгоритмов обработки и планирования наблюдений, используемых в этих программах, а также предложить ряд практически значимых усовершенствований указанных алгоритмов. Далее мы опишем поставленные задачи более подробно.
Выделение периодичностей из зашумленных данных. Первый этап анализа измерений лучевых скоростей звезд в программах поиска экзопланет — поиск в них периодических колебаний, потенциально отражающих наличие у звезды невидимых спутников (в том числе и планет). На сегодня уже разработаны весьма эффективмные методы такого анализа. В частности, в этой связи стоит упомянуть периодограмму Ломба-Скаргла [74, 84] и различные ее обобщения. Эти периодограммы позволяют легко и эффективно находить реальные периодичности в наблюдательных данных, однако не могут (сами по себе) оценивать статистическую значимость выявляемых сигналов (то есть, отличать ре-
альные выявляемые сигналы от шумовых эффектов). Используемые обычно на практике методы оценки статистической значимости пиков периодограмм обладают существенными недостатками (например, требуют слишком большого машинного времени или обладают недостаточной математической строгостью и обоснованностью и, следовательно, дают ненадежные результаты). Таким образом, необходима какая-то новая методика оценки значимости пиков периодограмм, которая позволяла бы получать требуемые оценки быстро, с достаточной точностью и с какими-либо теоретическими гарантиями надежности.
Дрожание лучевых скоростей звезд. При анализе временных рядов лучевых скоростей звезд в программах поиска экзопланет необходимо всегда иметь ввиду, что оценки инструментальной погрешности измерений не описывают полную дисперсию ошибок. Эта полная дисперсия содержит значительный добавочный компонент, который обычно называют «дрожание» («jitter»). Дополнительное кажущееся дрожание лучевой скорости изменяет соотношение статистических весов наблюдений и в конечном счете влияет на получаемые значения оценок орбитальных параметров и масс экзопланет.
В своей астрофизической части это дрожание вызывается различными явлениями активности в звездной атмосфере (и самой звезде), которые приводят к кажущейся нестабильности измеряемой лучевой скорости звезды. Наблюдаемая величина дрожания может существенно зависеть от характеристик спектрографа, условий и методики проведения наблюдений. Например, достаточно длительная экспозиция спектра (~ 20 мин) позволяет осреднить колебания видимой лучевой скорости, вызванные нерадиальными звездными осцилляциями [90, 79]. Однако, как будет показано в диссертации, свой вклад в дрожание вносят и чисто инструментальные эффекты. Например, оценки «внутренних» погрешностей измерения лучевой скорости, даваемые самими наблюдателями, совсем не обязаны быть точными. Если при расчетах не учитывался какой-то дополнительный шумовой эффект, то получаемые значения погрешностей измерений могут быть систематически заниженными. Наличие в измерительных ошибках неисключенных систематических компонент также увеличивает эффективную величину дрожания. Нельзя исключать и возможности, что даваемые погрешности измерений в действительности завышены.
Указанные причины приводят к тому, что эффективные значения дрожаний лучевых скоростей звезд сильно разнятся для различных инструментов (даже для одной и той же звезды). Это означает, что при совместном анализе данных, полученных на разных обсерваториях, мы должны назначать разным массивам данных существенно различные статистические веса с использований достаточно точных оценок эффективной величины дрожания (раздельно для каждого массива).
Обычно величина дрожания оценивается априорно на основе некоторых эмпирических корреляций с различными астрофизическими параметрами звезд, такими как, например, скорость вращения и уровень активности, измеряемый по Н и К линиям кальция [82, 83, 100]. Однако таким образом мы можем принять во внимание только астрофизическую часть дрожания, а как обсуждено выше, инструментальная компонента также значительна. При этом в литературе отсутствуют работы, посвященные исследованию
систематических погрешностей наблюдательных данных, получаемых в современных до-пплеровских программах поиска планет. Таким образом, при анализе таких временных рядов лучевых скоростей мы стоим перед проблемой корректного учета того факта, что дисперсии измерений содержат некоторый дополнительный компонент, величина которого различна для разных звезд и разных обсерваторий и плохо известна a priori. Этот факт необходимо корректно учитывать при нахождении оценок масс и орбитальных параметров планет.
Систематические инструментальные ошибки. На практике иногда оказывется, что наблюдаемая переменность лучевой скорости звезды вызвана не наличием у этой звезды планеты, а какими-либо явлениями активности в звездной атмосфере, например пятнами (см. например [67]). Однако в литературе фактически отсутствуют какие-либо достаточно подробные исследования инструментальных систематических ошибок высокоточных измерений лучевых скоростей, получаемых в программах поиска экзопланет. Вследствие таких ошибок получаемые оценки орбитальных параметров экзопланет могут оказаться сильно искаженными и ненадежными. В связи с этим возникает естественная необходимость оценки статистического качества таких наблюдений. Было бы полезно проанализировать возможность наличия систематических ошибок в опубликованных измерениях лучевых скоростей звезд с планетами, а также оценить величину и характер этих систематических ошибок.
Планирование наблюдений. При поиске внесолнечных планет методом лучевых скоростей практически не применялись какие-либо методы оптимального планирования наблюдений. Лишь в совсем недавней работе [53] опубликовано некоторое теоретическое исследование по этой теме. При этом на практике орбитальные конфигурации многих внесолнечных планетных систем определены довольно плохо, несмотря на наличие для этих систем уже многолетних рядов наблюдений. В математической статистике теория оптимального планирования экспериментов разработана уже весьма глубоко и подробно. Таким образом, естественно возникает стремление применить общие методы оптимального планирования экспериментов для планирования наблюдениіі лучевых скоростей в программах поиска экзопланет. Разработка таких алгоритмов и оценка их практической эффективности также составляет одну из целей данной работы.
1.3. Научная и практическая ценность результатов
Разработанные в диссертации методы и алгоритмы статистического анализа и планирования наблюдений имеют весьма важное значение для наблюдательных программ поиска экзопланет. Разработанные методы могут существенно повысить эффективность этих наблюдений. В результате можно повысить точность и надежность получаемых оценок планетных параметров и частоту открытий планет.
Построенные методики обработки наблюдений могут использоваться в других областях астрономии и других наук. В частности, новая методика оценки выявления (поиска)
периодичностей в наблюдательных (экспериментальных) данных, описанная в главе 2, имеет весьма универсальную область применения. Алгоритм оценки параметров модели с учетом неизвестного «дрожания» измеряемой величины (глава 3) можно использовать и в других схожих научных задачах, связанных с обработкой наблюдательных (экспериментальных) данных. Потенциально широкую область применения имеют и предложенные в диссертации алгоритмы оптимального планирования наблюдений (глава б).
Самостоятельную ценность имеют результаты анализа наблюдательных данных трехпланетной системы HD37124 (глава 5). Эти результаты имеют большое значение для небесной механики экзопланетных систем, а также для теории образования и миграции экзопланет.
1.4. Новизна и достоверность результатов
Сам объект исследования — внесолнечные планеты — является новым, так как массовые открытия экзопланет начались лишь около десятилетия назад. Отметим следующие наиболее важные новые результаты, полученные в данной работе:
Впервые получены практически эффективные (быстро вычисляемые ц имеющие хорошую точность) оценки статистической значимости периодичностей, выявляемых в зашумленных наблюдательных данных при помощи периодограммы Ломба-Скаргла и ее обобщений.
Впервые построена строгая методика учета эффективной величины «дрожания» лучевых скоростей звезд, основанная на методе максимального правдоподобия.
Впервые указано на то, что в высокоточных измерениях лучевых скоростей, получаемых в программах поиска экзопланет, присутствуют существенные систематические ошибки.
В ходе анализа наблюдательных данных для планетной системы HD37124 обнаружено качественно новое семейство допустимых орбитальных планетных конфигураций, соответствующее резонансу 2/1 между двумя внешними планетами.
Фактически впервые были разработаны строгие методы оптимального планирования наблюдений лучевых скоростей звезд для программ поиска внесолнечных планет.1
Достоверность полученных результатов обеспечивается использованием адекватных математических моделей и статистических методов, их проверкой различными численных тестами и сравнением с результатами других авторов в сопоставимых случаях.
гЛишь в 2008 г. по данной теме была опубликована работа зарубежных исследователей [53]. Эта работа была выполнена независимо и опубликована почти одновременно с соответствующей статьей соискателя. Впрочем, в указанной работе задача оптимального планирования наблюдений лучевых скоростей звезд решается существенно иными методами, так что аналогичные результаты соискателя не теряют своего
значения.
1.5. Результаты, выносимые на защиту
Новая эффективная методика оценки статистической значимости периодичностей, выявляемых при помощи периодограммы Ломба-Скаргла и ее обобщений.
Алгоритм учета кажущегося дрожания лучевых скоростей звезд, основанный на методе максимального правдоподобия и позволяющий получать более точные и надежные оценки параметров планет.
Обнаружение систематических ошибок годичного периода в измерениях лучевых скоростей, получаемых в современных программах поиска экзопланет.
Применение разработанных алгоритмов обработки наблюдений к опубликованным рядам измерений лучевой скорости звезды HD37124, вокруг которой обращаются (как минимум) три планеты-гиганта. Обоснование того, что две внешние планеты этой системы могут двигаться по орбитам с большими эксцентриситетами, не нарушая динамической устойчивости системы.
Алгоритмы оптимального планирования наблюдений лучевых скоростей в программах поиска экзопланет, наиболее эффективные для систем, содержащих две планеты или более, особенно при наличии резонансов.
1.6. Апробация работы
Результаты, полученные в ходе данного исследования, докладывались на семинарах Кафедры небесной механики СПбГУ, Института прикладной астрономии РАН, Главной (Пулковской) астрономической обсерватории РАН, ГАИШ МГУ, а также на нескольких всероссийских и международных научных конференциях: на 36-й, 37-й и 38-й международных студенческих научных конференциях «Физика космоса» (г. Екатеринбург, 2007, 2008 и 2009 г.); на всероссийской молодежной научной конференции «Физика и прогресс» (г. Санкт-Петербург, 2007 г.); на 249-м симпозиуме Международного Астрономического Союза «Exoplanets: Detection, Formation, and Dynamics» (г. Сучжоу, КНР, 2007 г.); на международной научной конференции «Extrasolar planets in multi-body systems: Theory and observations» (г. Торунь, Польша, 2008 г.).
1.7. Структура диссертации
Диссертация состоит из введения, пяти глав основного содержания, заключения, пяти приложений, списка литературы (102 ссылки), списка иллюстраций (24 рисунка) и списка таблиц (4 таблицы). В приложениях обсуждаются и решаются некоторые специальные математические задачи, возникшие в основном содержании. Общий объем диссертации составляет 138 страниц.
Во Введении (глава 1) описываются поставленные научные задачи и цели работы, обосновывается их актуальность, описываются основные научные результаты работы. Приводится описание математических обозначений, используемых в диссертации.
Вероятность ложной тревоги
Описанные периодограммы имели бы малую ценность без метода оценки статистической значимости их пиков. Поскольку наблюдения отягощены случайными ошибками, мы всегда находимся под угрозой совершения ошибок двух типов: ошибочное отвержение основной гипотезы Н (т.е. ошибочный вывод о существовании сигнала, ложноположи-тельный результат, «ложная тревога») и ошибочное неотвержение 7Ї (т.е. невыявление реальноіі периодичности, ложноотрицательный результат). Как правило, ложные тревоги считаются более опасными, поэтому критерием статистической значимости обычно выступает «вероятность ложной тревоги», то есть вероятность отвергнуть гипотезу Н, когда в действительности она верна. Мы будем обозначать эту вероятность как FAP (False Alarm Probability). Задавшись некоторым малым критическим значением FAP (скажем, 1%), мы можем заявить что найденная периодичность статистически значима (когда FAP FAP ) или незначима (когда FAP FAP ). Чтобы определить эту вероятность, нам необходимо найти функцию распределение той величины, которая является основой для используемого критерия (при этом считая гипотезу Н верной). В случае, когда частота сигнала заранее известна, критерием значимости выступает значение периодограммы z(f) на данной частоте /. В случае, когда частота заранее не известна, таким критерием выступает максимальный отсчет периодограммы в заданном частотном диапазоне. Обычно периодограммы строятся в диапазоне частот от нуля до некоторого предельного значения fmax. Поэтому нас будет интересовать распределение величины maxo / /mox z(f). Функцию распределения одиночных отсчетов какой-либо из введенных периодограмм мы будем обозначать как PSmg\e(z), а распределение указанного максимального отсчета — как Pmax(z, fmax). При этом соответствующие вероятности ложной тревоги выражаются как FAPSing]e( ) = 1 — Psing\o(z) " Г Л1 maxV- i /max/ = і мпах і /max/ Одиночные значения введенных выше ЛМНК-периодограммы и ее модификаций представляют собой (с точностью до некоторых простых коэффициентов) классические статистики из задачи линейной регрессии, так что их функции распределения хорошо известны [10, 7.1]. Мы можем использовать тот факт, что распределение одиночных отсчетов 2z(f) есть х2 распределение с d степенями свободы, распределение одиночных значений 2z2/d есть F-распределение с d и NJC степенями свободы, а одиночных значений 2ZI/N-H — і?(бета)-распределение с теми же числами степеней свободы.
Явное выражение вероятности ложной тревоги для одиночных отсчетов базовой ЛМНК-периодограммы выглядит как где Тх(п) есть (нормализованная) неполная гамма-функция. В частности, в случае периодограммы Ломба-Скаргла (d = 2 параметра синусоидального сигнала) мы получаем широко известное экспоненциальное выражение FAPSing]e = e z. Для модифицированной ЛМНК-периодограммы zi(f) аналогичное выражение выглядит как где IX(P,Q) есть (нормализованная) неполная бета-функция. Для периодограмм z2 и z3 соответствующие выражения можно построить по (2.12), просто выражая z\ через , или z3 с использованием связей (2.10). Заметим, что для случая d = 2 значения периодограммы , оказываются распределены в точности так же, как и значения периодограммы z (т.е., экспоненционально). Для распределения максимального отсчета даже простейшей периодограммы Ломба-Скаргла не существует точного выражения. Довольно часто на практике применяется формула «повторных испытаний» которая была бы верной, если бы максимальный отсчет периодограммы вычислялся среди конечного дискретного набора из Л статистически независимых отсчетов. К сожалению, такая система «независимых частот» существует только в случае равномерного распределения наблюдений во времени. Для неравномерных временных рядов такую систему построить вряд ли возможно. Неравномерное распределение моментов наблюдений вызывает явление элайзинга (aliasing), которое проявляется, в частности, в коррелированности различных (даже довольно удаленных) отсчетов периодограммы [4]. Например, в работе [70] автору не удалось найти в неравномерном случае более двух независимых отсчетов периодограммы. Кроме того, даже в случае равномерного временного ряда, «независимые частоты» слишком разрежены, так что существенные детали периодограммы проваливаются в промежутки между ними. Поэтому периодограммы строятся, как правило, на значительно более плотной сетке, которая практически эквивалентна непрерывному отрезку частот.
В этом случае понятие «числа независимых частот» теряет физический смысл, а формула (2.13) становится лишь эвристической аппроксимацией [55]. При этом мы уже не имеем сколько-нибудь приемлемой аналитической оценки для параметра Nint\ и, следуя авторам работы [68], вынуждены оценивать этот параметр методом Монте-Карло, используя (2.13) как априорную модель для получаемой эмпирической функции распределения максимального отсчета периодограммы. При этом у нас нет никаких теоретических гарантий того, что формула (2.13) вообще верна хотя бы асимптотически в каких-то предельных условиях. В статье соискателя [21] предложен другой метод оценки искомого распределения, основанный на результатах теории экстремальных значений случайных процессов. Основная идея этого подхода — рассматривать периодограмму наблюдений, отягощенных случайных ошибками, как случайный процесс (см. детали в приложении А в настоящей диссертации). Для оценки вероятности ложной тревоги, связанной с максимальным пиком
Переход к методу максимального правдоподобия
В работах соискателя [25, 24] разработан другой алгоритм одновременного оценивания дрожания сг и параметров 9, который позволяет решить описанные выше проблемы. Этот алгоритм основан на методе максимального правдоподобия. В случае, когда ошибки лучевой скорости нормально распределены и статистически независимы, логарифм функции правдоподобия может быть записан где Oi даются квадратично-аддитивной моделью с неизвестным параметром р = сг . В более общем случае мы можем работать со смешанными данными, полученными на различных обсерваториях. Пусть на j-й обсерватории получено N3- наблюдений лучевой скорости данной звезды. Всего имеется s таких однородных подмассивов. При этом N = YlSj=iNj. Соответствующие моменты наблюдений, сами лучевые скорости, их предполагаемые и полные инструментальные погрешности индексируются как tji, VJU crmeaSjJj, aji (j = 1, 2,..., s; і — 1, 2,..., Nj). Для каждого из s подмассивов мы должны ввести свой параметр дрожания pj, которые все вместе составляют вектор р. В этом случае функция правдоподобия зависит от полного вектора из d+s параметров = {0,р} и записывается как в котором постоянные члены COJ различны для разных массивов данных. Точка максимума функции правдоподобия () (или функции 1п(), что эквивалентно), дает требуемую оценку вектора параметров . Эти оценки мы будем называть ММП-оценки (оценки, полученные методом максимального правдоподобия).
Заметим, что если бы мы приняли мультипликативную модель для CTJ (с единственным однородным массивом данных, то есть при s — 1), то задача максимизации функции правдоподобия свелась бы к задаче минимизации функции х {&) = ку?{в,к) по в. При этом в = argminx2(#) есть МНК-оценка, а ММП-оцснка для к дается явным выражением /с = їгт/N. Таким образом, в каждом конкретном случае мы можем получать одновременную оценку для обычных параметров кривой лучевой скорости и для дрожания лучевой скорости, не привязываясь к каким-либо эмпирическим моделям для сг . При этом оценки параметров в и р автоматически учитывают информацию друг о друге. Похожий подход, основанный на методе максимального правдоподобия, применялся в работе [14] для оценки дисперсий скоростей в ансамблях звезд. Итак, мы можем строить совместные ММП-оценки для d + s параметров в и р. Однако любая оценка имеет малую ценность без сопутствующей оценки статистической погрешности. В данном случае нам нужно оценить ковариационную матрицу полного вектора . Ковариационная матрица ММП-оценок обычно определяется на основе информационной матрицы Фишера где математические ожидания должны вычисляться в предположении, что данное значение в истинно. Обратная матрица, F-1, представляет собой асимптотическое (при N — со, то есть для больших выборок) приближение ковариационной матрицы [11, 6.4]. В ходе вычисления мы обнаруживаем, что полная матрица F представляется в диагонально-блочной форме, причем первый блок в (левом-верхнем углу F) равен имеет размер d х d и совпадает с матрицей Фишера для МНК-оценок в. Остальные ненулевые элементы F находятся только на диагонали и соответствуют параметрам дрожания. Это означает, что ММП-оценки параметров в и р асимптотически некоррелированы, равно как и отдельные оценки Pj.
При этом ковариационная матрица оценок обычных параметров кривой лучевой скорости асимптотически выражается как Var0 Q-1, то есть так же, как и в случае МНК-оценок, а дисперсии оценок дрожания выражаются (асимптотически) как ММП-оценки обладают многими полезными асимптотическими свойствами при N —» со. Они являются асимптотически несмещенными (математические ожидания стремятся к истинным значениям параметров). Их распределение асимптотически нормально, с указанными выше выражениями для соответствующей ковариационной матрицы. Они являются асимптотически эффективными, то есть их статистические погрешности в некотором смысле минимальны [11, глава б].1 Несложно показать, что дисперсия ММ-оценки параметра р3 равна 2 а / А N, г Таким образом, отношение дисперсий ММ- и ММП-оценок pj равно Вследствие неравенства Коши-Шварца это отношение не меньше единицы. Таким образом, метод максимального правдоподобия дает статистически более точную оценку дрожания, чем метод моментов. Метод максимального правдоподобия организован так, что получающееся значение функции X2(6 )/N всегда близко к единице. Это означает, что мы уже не можем использовать х2-функцию как меру качества подгонки наблюдаемых данных моделью ц. Вместо этого, мы должны построить другую меру согласия модели с наблюдениями, основанную на функции правдоподобия (). Чтобы эта мера была наглядной, она должна измеряться в тех же единицах, что и лучевая скорость (то есть, в м/с). Значит, такая величина должна быть пропорциональна C llN . Чтобы найти подходящий нормализующий множитель, предположим на секунду, что х2 = N в точности. Тогда -1 — егєеотег/2\/2Іт, где сгєеот — геометрическое среднее полных дисперсий jjj. Значит, уже в общем случае, в качестве меры качества модели кривой лучевой скорости удобно использовать функцию При этом оценка максимального правдоподобия доставляет минимум /(). Впрочем, в целях сравнения с работами других авторов, мы будем использовать также более традиционную меру согласия — среднеквадратичное уклонение невязок (r.m.s.).
Смещение оценок орбитальных параметров
Хорошо известно, что в линейном случае МНК-оценки являются несмещенными. В случае нелинейных моделей, а также при переходе к методу максимального правдоподобия, может возникать некоторое смещение. Также известно, что ММП-оценки обладают свойством асимптотической несмещенности, то есть их смещение стремится к нулю при N — со. Несмотря на это свойство, на практике смещение ММП-оценок может оказаться значительным, даже если число наблюдений на первый взгляд довольно велико. При этом с ростом числа наблюдений смещение все же убывает как const /N, но величина соответствующей константы оказывается довольно большой. Эта величина зависит от ряда факторов: числа свободных параметров, степени и характера нелинейности подгоночной модели. Какие эффекты может давать смещение оценок орбитальных параметров планетных систем? Один из легко предсказуемых эффектов — переоценка значений орбитальных эксцентриситетов. Рассмотрим пример планеты на строго круговой орбите. Со ответствующии эксцентриситет в точности равен нулю, но его оценка всегда будет строго положительна. При этом, очевидно, оценка эксцентриситета будет обладать некоторым положительным смещением. В случае, когда число наблюдений велико, амплитуда колебания лучевой скорости существенно больше погрешностей, а планета — единственная в системе, величина смещения оценки ее эксцентриситета вряд ли будет существенна. Однако если система содержит несколько планет, а число наблюдений не очень велико, то оценки орбитальных эксцентриситетов могут оказаться значительно смещенными. Ситуация становится еще хуже, если в системе есть резонансные планетные пары. Тогда параметры резонансных планет оказываются сильно коррелированы, их точность заметно падает, а статистическое смещение соответственно растет. Таким образом, возникает естественная потребность в редукции смещения оценок. Можно предложить несколько методов редукции смещения оценок: 1. Коррективный аналитический. Если заданы уравнения подгоночной модели, то можно вычислить смещение оценок напрямую (разумеется, с точностью до некоторого — обычно лишь первого — порядка малости). Это можно сделать, используя степенные разложения логарифма функции правдоподобия в окрестности истинных значений параметров [33]. В приложении В такая работа выполняется для модифицированной функции правдоподобия (3.15) с учетом возможной негауссовости ошибок измерений. Затем эту аналитическую оценку смещения можно вычесть из текущих оценок, то есть провести коррекцию смещения постфактум. К сожалению, аналитические выражения смещения сложно использовать на практике, так как они содержат производные второго порядка от уравнений подгоночной модели.
Выражения вторых производных даже простейшей кеплеровской модели (1.1) слишком сложны и громоздки, а получить их трудно. Впрочем, выражение для смещения оценки дрожания лучевой скорости (3.16) все же можно использовать на практике. 2. Превентивный ауьалитический. Можно попытаться внести коррекцию в саму функцию правдоподобия так, чтобы смещение оценок получаемых в результате максимизации новой функции, имело возможно меньший порядок малости. Такой подход можно назвать превентивным [49]. Мы использовали этот подход для превентивной редукции смещения оценки дрожания лучевой скорости. К сожалению, возможности данного подхода сильно ограничены тем, что на практике очень редко удается внести правильную коррекцию в функцию правдоподобия. Какого-либо общего алгоритма построения такой коррекции не существует. 3. Численный. На основе текущих оценок параметров можно провести моделирование Монте-Карло с тем, чтобы эмпирически оценить смещение полученной оценки и затем вычесть его. Главный недостаток этого метода очевиден — это непомерно большой объем вычислений. По-видимому, этот метод нет смысла применять в практике определения орбит экзопланет методом лучевых скоростей. 4. Метод Кенуя. Этот метод имеет много названий, но впервые он, вероятно, был предложен Кенуем [81]. В англоязычной литературе он чаще называется «jackknife» или «leave-one-out (loo)» [42], а на русский язык его обычно переводят как «метод складного ножа» [17]. Этот метод полуаналитический — он требует повторных решений задачи максимизации функции правдоподобия с некоторыми модельными временными рядами, но содержит элементы аналитики, которые позволили сильно сократить число таких испытаний по сравнению собычным моделированием Монте-Карло. Этот метод мы опишем ниже. Использовать метод Кенуя при оценке параметров экзопланет было предложено соискателем в статьях [24, 25]. Теоретическим элементом метода Кенуя является асимптотическое представление смещения в виде ряда
Идея состоит в том, чтобы эмпирически сравнить величину смещения оценок, полученных по массивам данных с разными N и затем по формуле (4.2) провести экстраполяцию оценок на несмещенный случай 1/N — 0. Простейшая форма алгоритма Кенуя (первого порядка) имеет следующий вид. 1. Вычислить основную оценку на основе полного временного ряда. Эта оценка смещена на величину bi/N. 2. Построить N новых временных рядов путем выкидывания ровно одного из N наблюдений. Таким образом, каждый такой временной ряд будет содержать (N — 1) точек. 3. Для каждого из этих сокращенных временных рядов построить собственную оценку параметров (г = 1,2,..., N), используя в точности тот же алгоритм, что и для получения основной оценки . Полученные оценки будут смещены на чуть большую величину bi/(N — 1). По отношению к основной оценке они будут в среднем смещены на величину b\(l/N — 1/(N — 1)) by/N2. 4. Составить сумму b\ = І=:1( І— ) Эта величина (уже сама по себе, без деления на N) представляет собой оценку смещения первого порядка, 0(1/N). Исправленная оценка — &i будет смещена лишь на величину порядка 0(l/N2). Достоинства алгоритма Кенуя состоят в том, что его практическая реализация проста и модельно независима. Кроме того, указанный алгоритм не требует дополнительной информации, такой как, например, степень нормальности измерительных ошибок. Его можно применять к оценкам любых параметров: как обычных параметров кривой лучевой скорости, так и параметров дрожания лучевой скорости. Недостатком алгоритма является все еще довольно большой объем необходимых вычислений, требующий решения серии нелинейных задач на максимизацию.
Учет требования динамической устойчивости
Сначала мы применим приблизительно такой же подход, как использованный в работах [59, 57]. Мы продолжим использовать описанный выше метод построения двумерных графиков частично минимизированной статистики I, но теперь дополним его проверкой динамической устойчивости получаемых орбитальных конфигураций. Для каждой точки на плоскости выбранных параметров мы не только находим соответствующие ММП-оценки оставшихся параметров, но и будем проводить численное интегрирование уравнений движения планет с тем, чтобы отсеять быстро разрушающиеся планетные конфигурации. Чтобы выполнять эту процедуру интегрирования мы должны знать истинные массы планет тп. Как можно видеть из выражений (1.2), эти массы зависят от массы звезды М . Следуя работе [51], мы принимаем значение М и О.78М0 с возможной погрешностью 10%. Кроме того, мы должны принять какое-то значение для наклона орбит системы. Этот наклон нельзя определить из измерений лучевых скоростей, по крайней мере пока взаимные планетные возмущения не начнут достаточно ясно проявляться в наблюдаемой кривой лучевой скорости. Поэтому лучшее, что мы можем сделать — это считать, что орбиты в системе лежат в одной плоскости и видны с ребра, то есть гп = 90. Если действительный угол наклона менее 90, то истинные массы планет будут больше, чем в случае гп — 90. Это увеличение масс вызовет некоторое сжатие областей устойчивости. Аналогичный эффект вызывается ненулевыми взаимными наклонами орбит планет (из-за известного явления сопряжения наклонов с эксцентриситетами). Таким образом, те области устойчивости, которые мы получим из предположения гп = 90, будут представлять наиболее «оптимистичный» из возможных вариантов.
Интегрирование уравнений движения трехпланетной системы выполнялось сим-плектическими интеграторами Йошиды 4, б и 8 порядков [102]. Как известно, симплекти-ческие интеграторы обладают свойством более качественного сохранения интеграла энергии и фазового портрета динамической системы на больших временах. Ошибка интеграла энергии почти не накапливается со временем. Поэтому эти интеграторы удобны нам здесь для проверки устойчивости планетной системы на относительно больших временых. При этом их слабостью является невозможность изменения длины шага непосредственно в ходе интегрирования. В частности, при обработке тесных сближений тел нельзя уменьшить этот шаг с тем, чтобы сохранить приемлемую точность. Критерием остановки интегрирования («распада» планетной системы) считалось осуществление одного из следующих условий: Столь большие изменения областей движения планет вряд ли совместимы с сохранением структуры планетной системы. Наибольшие трудности в получении надежной орбитальной конфигурации системы связаны со слишком большим эксцентриситетом внешней планеты. Поэтому рассмотрим сначала плоскость переменных (eccosgc,ecsingc). Соответствующие графики построены на рис. 5.5. На них хорошо видна сложная структура поверхности правдоподобия. На графиках, соответствующих семейству решениіі А, можно наііти по 2-3 локальных минимума /; все они соответствуют довольно большим значениям ес. На графиках для семейства В есть только один минимум, но также соответствующий большому ес. Ни один из этих локальных минимумов не соответствует устойчивой конфигурации. Устойчивые решения занимают область малых и умеренных ес. Важно отметить, что слабо эксцентрические решения соответствуют только конфигурациям типа А: с уменьшением ес решения типа В приближаются к резонансу 2/1 и плавно переходят в класс А. По рис. 5.5 можно легко найти устойчивые планетные конфигурации из класса А. Такие конфигурации могут давать значения I « 8.6 м/с (для модели I), 8.2 м/с (модель II) и 8.1 м/с (модель III).
Интересно, что область устойчивости в классе А не является строго симметричной и несколько коррелирует с одним из локальных минимумов /, расположенном в третьей четверти координатной плоскости. Что касается класса В, то соответствующая область устойчивости очень узка и находится у границы с классом А. Тот факт, что почти все решения класса В неустойчивы, можно объяснить тем, что отношение формальных оценок орбитальных периодов Pc/Pd 2.1 — 2.3 довольно мало и без посторонней помощи не попадает на какой-либо резонанс низкого порядка с достаточной точностью. Для того, чтобы найти область устойчивости в классе В, мы используем другую пару переменных, а именно пару (Рс,ес). Соответствующие графики частично минимизированной функции I построены на рис. 5.6. Можно заметить, что значения І в зоне малого ес и большого Рс довольно велики, так что соответствующие орбитальные конфигурации маловероятны. Это показано более наглядно на рис. 5.7. При построении этих зависимостей минимизация I осуществлялась так, что эксцентриситет ес был жестко зафиксирован на безопасном значении 0.15, а отношение орбитальных периодов фиксировалось на значениях, отображенных на оси абсцисс графиков. Так как эксцентриситет ес фиксирован, мы можем ожидать, что эффекты нелинейности теперь заметно подавлены, и что можно попробовать применить для величины Pc/Pd стандартную теорию построения асимптотических (N — оо) доверительных интервалов, основанную на статистике отношения правдоподобия и ее асимптотическом распредении. Как следует из рис. 5.7, при использовании модели I значения Рс/Ра 2.5 и Pc/Pd 1-87 находятся за пределами 99-процентного доверительного интервала. Впрочем соответствующий график на рис. 5.7 имеет два локальных минимума, что указывает на все еще существенную нелинейность.
