АСТРОНОМИЧЕСКИЙ ЖУРНАЛ, 2019, том 96, № 1, с. 70-74

УДК 52-17

ВЛИЯНИЕ ЦИФРОВОГО ШУМА НА РЕЗУЛЬТАТЫ ИНТЕРПРЕТАЦИИ

ТРАНЗИТНОЙ КРИВОЙ БЛЕСКА

¹Московский государственный университет имени М. В. Ломоносова,

Государственный астрономический институт имени П. К. Штернберга, Москва, Россия

Поступила в редакцию 29.06.2018 г.; принята в печать 13.09.2018 г.

Часто используемые алгоритмы интерпретации кривых блеска, например, распространенный алгоритм

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

нефизическому вкладу в невязку (цифровому шуму). На примере интерпретации транзитной кривой

блеска двойной системы HD 209458 продемонстрирована необходимость учета цифрового шума.

Показано, что улучшение точности вычисления кривой блеска позволяет более надежно решать

вопрос об адекватности наблюдаемой кривой блеска используемой модели, поскольку исключается

нефизический вклад в невязку, связанный с ошибками вычислений. Интернет-ссылка на разработан-

ный алгоритм приведена в Заключении.

DOI: 10.1134/S0004629919020014

1. ВВЕДЕНИЕ

При этом во многих случаях при решении физи-

ческой задачи путем интерпретации кривой блеска

В настоящее время проводятся массовые вы-

приходится учитывать ошибки нефизической при-

сокоточные наблюдения транзитных кривых блес-

роды — погрешность вычислений, ошибки, связан-

ка, в связи с чем актуальна задача их надежной

ные с округлением (цифровой шум). Однако с

интерпретации [1, 2]. При этом представляется

учетом возможностей современных вычислитель-

возможным существенное увеличение возможно-

ных машин и точности современных наблюдений

стей интерпретации (анализа кривых блеска) по

блеска, содержащих десятки тысяч точек, необ-

сравнению с теми, что обеспечивают остальные

ходимы вычисления с лучшей точностью. Исполь-

алгоритмы. Прежде всего это связано с улучше-

зуя квадратурную формулу Гаусса для нахождения

нием точности вычисления кривой блеска за счет

численных значений интегралов, можно за при-

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

емлемое время получить результат с точностью,

Гаусса [3, 4].

соответствующей современной точности чисел с

плавающей точкой (около 18 значащих цифр).

Большинство используемых на данный момент

алгоритмов основано на алгоритме JKTEBOP, со-

Улучшение точности вычисления позволяет

зданном еще в начале 1980-х годов [5]. В алго-

прежде всего более качественно решать вопрос

ритме JKTEBOP вычисление интегралов в выра-

об адекватности наблюдаемой кривой блеска

жениях для блеска затменной двойной системы

используемой модели. Поскольку исключается

осуществляется непосредственно — путем прямо-

нефизический вклад в невязку, связанный с ошиб-

го суммирования по участкам, на которые разби-

ками вычислений, можно делать более надежные

вается область интегрирования. При этом ошибка

выводы о том, насколько значимо изменяется

в вычислении интеграла обратно пропорциональна

невязка в тех или иных случаях. Например,

можно рассмотреть, насколько значимо изменение

количеству производимых элементарных действий

минимального значения невязки при добавлении

(соответственно времени вычисления), а точность

параметров, по которым производится интерпрета-

результата существенно ограничена временем, ко-

ция, в том числе, как зависит минимальная невязка

торые можно отвести на вычисления. Такая точ-

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

ность соответствовала возможностям ЭВМ, кото-

рые существовали во время создания JKTEBOP, и

Следует отметить, что улучшение точности вы-

точности наблюдений того времени.

числения кривой блеска может существенно вли-

ять на значение минимальной невязки даже при

^*E-mail: marat@sai.msu.ru

отсутствии значимого изменения в значениях самих

^**E-mail: ngostev@mail.ru

искомых параметров. Это можно объяснить тем,

ВЛИЯНИЕ ЦИФРОВОГО ШУМА

1.000

0.995

0.990

0.985

0.980

165

170

175

180

185

190

195

Θ(°)

Рис. 1. Наблюдаемая (точки) и теоретическая (сплошная линия) кривые блеска двойной системы с экзопланетой

HD 209458 из работы [6]. Внизу показаны отклонения наблюдаемых значений блеска от теоретической кривой блеска,

рассчитанной в рамках модели с нелинейным (квадратичным) законом потемнения к краю.

что цифровой шум может учитываться в ошибках

1.13 × 10-4 до 2.47 × 10-4 (в долях внезатменной

искомых параметров, получаемых методом Монте-

интенсивности). Величины относительных ошибок

Карло. При таких обстоятельствах можно сделать

(в долях глубины затмения) лежат в пределах от

вывод, что алгоритм JKTEBOP в определенных

∼7 × 10-3 до ∼1.5 × 10-2. Интерпретация кривой

случаях позволяет получать на практике удовле-

блеска выполнена в линейном и квадратичном

творительные результаты, а именно, значения па-

законе потемнения звездного диска к краю.

раметров системы, оцениваемые методом Монте-

Карло (в предположении априорной адекватности

Кривая блеска двойной системы Kepler-15b,

также использованная нами для тестирования

модели наблюдательным данным). В то же время,

если речь идет об интерпретации без такого пред-

алгоритма, взята из каталога NASA Exoplanet

положения (о верности модели) — когда возможно

Archive. Использовались наблюдательные данные

пустое доверительное множество или когда реша-

короткой каденции (short cadence)

18-го сета

ется вопрос о статистически значимом изменении

наблюдений, выполненного с

сентября по

невязки, — требуется более точный расчет значе-

18 октября 2009 г.

ния невязки.

2. НАБЛЮДАТЕЛЬНЫЕ ДАННЫЕ

3. ОПИСАНИЕ МОДЕЛИ

Для демонстрации работы разработанного

Мы использовали модель двух сферических

нами алгоритма использовалась наблюдаемая

звезд на круговой орбите, предполагалось отсут-

кривая блеска двойной системы с экзопланетой

HD 209458 из работы [6]. Наблюдаемая кривая

ствие эффектов отражения и эллипсоидальности.

блеска, представленная в [6], получена на HST

Геометрия модели изображена на рис.

2, где

в апреле-мае 2000 г. Спектры получены с ис-

R_∗ — радиус звезды, R_o — радиус планеты, D —

пользованием спектрометра STIS со спектральной

расстояние между центрами дисков компонентов,

ρ, Ψ — полярный радиус и полярный угол точки

решеткой G750M. Наблюдения проводились в

на диске звезды соответственно. Начало координат

диапазоне 5813-6382

A с разрешением R =

находится в центре звезды, полярный угол отсчи-

= λ/Δλ = 5440

(подробнее см. [6]). Норми-

тывается против часовой стрелки.

рованная кривая блеска транзита экзопланеты

по диску звезды представлена на рис. 1. Кривая

Для круговой орбиты расстояние между центра-

блеска включает в себя 556 индивидуальных значе-

ми дисков звезд Δ зависит от фазы θ и угла наклона

ний блеска двойной системы. Среднеквадратичная

орбиты i как

√

ошибка индивидуального измерения σobsi для раз-

Δ(θ, i) = cos² i + sin² i sin² θ.

(1)

ных точек кривой блеска заключена в пределах от

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

АБУБЕКЕРОВ, ГОСТЕВ

R_*

R_o

dρ

Рис. 2. Модель затменной двойной системы. Проекция на картинную плоскость.

При расчете кривой блеска в качестве функции

Полный блеск звезды, который совпадает с пол-

распределения яркости по диску звезды использо-

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

вался линейный закон потемнения к краю диска с

∫

(

линейным коэффициентом потемнения к краю x

x⁾

(

√

)

L^full = 2π I(ρ)ρdρ = πR²I0

(4)

∗

I(ρ) = I₀

1-x+x

(2)

R2∗

в модели с линейным законом потемнения к краю и

и квадратичный закон потемнения к краю диска,

∫

(

y⁾

отличающийся от линейного дополнительным сла-

L^full = 2π I(ρ)ρdρ = πR²I0

(5)

∗

гаемым, содержащим квадратичный коэффициент

потемнения к краю y:

(

√

)

в модели с квадратичным законом потемнения к

краю.

I(ρ) = I₀

1-x

(3)

R2∗

Яркость в центре диска I₀ выбирается таким

(

√

)₂)

образом, чтобы полный блеск системы был равен

единице (условие нормировки).

-y

Падение блеска при затмении есть

R2∗

∫∫

Здесь ρ — полярное расстояние от центра диска

L^dec(Δ,R_∗,R_o,x,y) =

I(S)dS,

(6)

звезды, I₀ — яркость в центре диска, а R_∗ — ра-

S(Δ)

диус диска звезды. Яркость в центре диска ком-

понента 1 (звезды) далее будем обозначать как

где S(Δ) — область перекрытия дисков (завися-

I₀. Яркость I₀ в центре компонента 2 (планеты)

щая от расстояния между их центрами). Основной

и соответственно яркость в любой точке ее диска,

задачей при вычислении кривой блеска является

предполагается равной нулю. Компонент 2 (плане-

вычисление падения блеска для затмения.

та) в орбитальной фазе θ = π затмевает компонент

Модель двух сферических звезд для интерпре-

1 (звезду). Единицей измерения длины в наших мо-

тации транзитной кривой блеска, полученной в

делях является расстояние между центрами звезды

оптическом диапазоне, физически обоснована. В

и планеты a, a = 1, орбита считается круговой.

оптическом диапазоне звезда и планета обладают

“Третий свет” в модели отсутствует. Искомыми

резким краем. Эффекты, связанные с деформаци-

параметрами модели являются радиусы звезды R_∗

ей атмосферы экзопланеты звездным ветром, или

и планеты R_o, угол наклона орбиты i, коэффициент

эффект кометного хвоста, в оптическом диапазоне

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

не проявляются или пренебрежимо малы. Поэтому

закона потемнения к краю также и коэффициент

аппроксимация звезды и диска экзопланеты кругом

потемнения y.

для оптического диапазона удовлетворительна.

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

ВЛИЯНИЕ ЦИФРОВОГО ШУМА

1.000

0.998

0.996

0.994

0.992

0.990

176

178

180

182

184

Θ(°)

Рис. 3. Наблюдаемая (точки) и теоретические кривые блеска двойной системы с экзопланетой Kepler-15b. Сплошной

линией показана кривая блеска, рассчитанная в линейном законе потемнения к краю, штриховой линией — в квадратич-

ном.

4. ОПИСАНИЕ АЛГОРИТМА

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

краю использование нашего алгоритма дает невяз-

При интерпретации данных нами использовал-

ку χ² = 1.0134, что соответствует максимальному

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

уровню значимости 46%. Таким образом, исполь-

системы, подробно изложенный в работах [3, 4].

зование более точного алгоритма, подавляющего

В данном алгоритме интегрирование по области

цифровой шум, то есть без вклада, вызванного

перекрытия дисков выполняется путем аналитиче-

ошибками вычисления, позволяет сделать вывод

ского вычисления интегралов через эллиптические

о статистически значимом уменьшении невязки

функции (для вычисления которых с любой задан-

при введении в алгоритм квадратичного закона

ной точностью существуют эффективные алгорит-

потемнения к краю. Также в результате устране-

мы), либо путем применения формулы квадратур

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

Гаусса после определенных аналитических преоб-

возможность работать с непустым доверительным

разований подынтегрального выражения.

множеством на более высоком уровне значимо-

сти (более низком уровне доверия). Отметим, что

5. РАСЧЕТЫ

вывод о предпочтительности квадратичного закона

потемнения к краю для системы HD 209458 сделан

В работе Соузворза [5], в которой использу-

и в работе [7], однако там этот вывод в б ´ольшей

ется алгоритм интерпретации кривой блеска, ос-

степени основан на нефизической скоррелирован-

нованный на JKTEBOP, для системы HD 209458

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

значение невязки в линейном законе потемнения

нейном законе потемнения к краю. Уменьшение же

к краю составляет χ² = 1.1457, что соответствует

невязки при переходе от линейного к квадратично-

максимальному уровню значимости 1%. Отметим,

му закону потемнения к краю рассматривается как

что интерпретация в указанной работе ведется по

незначительное. Использование нашего алгорит-

следующим параметрам: радиусу звезды и планеты,

ма, позволяющего пренебречь ошибками округле-

наклонению орбиты, эксцентриситету и четырем

ния, показывает, что уменьшение невязки вызывает

коэффициентам потемнения к краю.

изменение максимального уровня значимости с 6%

При использовании же нашего алгоритма [3, 4]

до 46%, а это, очевидно, значительное изменение.

в предположении линейного закона потемнения к

Таким образом, устранение цифрового шума поз-

краю мы получаем значение невязки χ² = 1.103,

воляет выявить предпочтительность квадратичного

что соответствует максимальному уровню значи-

закона потемнения к краю для системы HD 209458

мости 6%. Хотя сами значения искомых пара-

по одному лишь уменьшению невязки.

метров согласуются в пределах соответствующих

ошибок, мы имеем заметное различие в минималь-

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

ных значениях невязок.

кривую блеска Kepler-15b (см. рис. 3), для которой

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019

АБУБЕКЕРОВ, ГОСТЕВ

не наблюдается значимого уменьшения невязки

Алгоритм предполагается использовать при

при переходе к квадратичному закону потемнения

массовой интерпретации транзитных кривых блес-

к краю, несмотря на видимые различия между

ка двойных с экзопланетами для получения эмпи-

соответствующими кривыми блеска.

рических значений коэффициентов потемнения к

краю звезд разных спектральных классов [2].

6. ОБСУЖДЕНИЕ

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

7. ЗАКЛЮЧЕНИЕ

ходимость применения алгоритмов, обеспечиваю-

щих высокоточное вычисление кривой блеска для

Разработанный авторами алгоритм интерпре-

ответственных суждений об адекватности моде-

тации транзитных кривых блеска общедоступен.

ли наблюдательным данным при интерпретации

Использование алгоритма позволит более каче-

транзитных кривых блеска. Разработанный нами

ственно анализировать богатый наблюдательный

алгоритм на основе квадратурных формул Гаусса

материал, поставляемый наземными и космически-

позволяет работать на точности, соответствующей

ми обсерваториями. С алгоритмом можно озна-

современной разрядности чисел с плавающей за-

комиться по адресу http://lnfm1.sai.msu.su/

пятой (19 значащих цифр). Такая точность вы-

~ngostev/algorithm.html.

числения кривой блеска позволяет более надежно

делать вывод об адекватности наблюдаемой кривой

блеска используемой модели.

Использование разработанного нами алгорит-

БЛАГОДАРНОСТИ

ма практически исключает нефизический вклад в

невязку, связанный с ошибками вычислений, и

Авторы выражают благодарность А.М. Черепа-

позволяет делать более надежные выводы о том,

щуку за плодотворное обсуждение работы.

насколько значимо изменяется невязка в тех или

иных случаях. Например, можно понять, насколько

значимо изменение минимального значения невяз-

СПИСОК ЛИТЕРАТУРЫ

ки при добавлении параметров, по которым про-

изводится интерпретация, в том числе, как зависит

1. М. К. Абубекеров, Н. Ю. Гостев, А. М. Черепащук,

минимальная невязка от используемых законов по-

Астрон. журн. 85(2), 121 (2008).

темнения к краю. Кроме того, интерес представляет

2. М. К. Абубекеров, Н. Ю. Гостев, А. М. Черепащук,

и сама по себе возможность работать на более

Астрон. журн. 87(12), 1199 (2010).

высоком уровне значимости (более низком уровне

3. M. K. Abubekerov and N. Yu. Gostev, Monthly Not.

доверия).

Roy. Astron. Soc. 432, 2216 (2013).

Следует отметить, что вопрос адекватности тео-

4. M. K. Abubekerov and N. Yu. Gostev, Monthly Not.

ретической кривой блеска наблюдательным дан-

Roy. Astron. Soc. 459, 2078 (2016).

ным тесно связан с обнаружением физических

5. J. Southworth, Monthly Not. Roy. Astron. Soc. 379,

эффектов, описываемых используемой моделью.

L11 (2007).

Высокая чувствительность критерия χ² к откло-

6. T. M. Brown, D. Charbonneau, R. L. Gilliland,

нениям наблюдательных данных от той или иной

R. W. Noyes, and A. Burrows, Astrophys. J. 552, 699

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

(2001).

физические эффекты, но и вместе с тем требует

7. J. Southworth, Monthly Not. Roy. Astron. Soc. 386,

достаточно точного вычисления статистики χ².

1644 (2008).

АСТРОНОМИЧЕСКИЙ ЖУРНАЛ том 96

№1

2019