Некоторые проблемы анализа радиоастрономических данных

Лекция для студентов (полная gzipped PostScript версия (950k) здесь)

Оглавление

• Введение
• Формула эксперимента
• Построение алгоритмов
  • Параметры распределений
    • Распределения
  • Компьютерное моделирование распределений
  • Фурье-анализ
  • Теорема отсчетов
  • Робастность
  • О восстановлении изображений
    • Метод максимальной энтропии
    • CLEAN
• Обзоры
  • Классификация обзоров
  • Статистические характеристики обзоров
  • Полнота
  • Имена объектов
    • Рекомендации МАС
• Проблемы архивизации
  • FITS
  • Расширения FITS
• Системы обработки данных
  • Linux и астрономическое программное обеспечение
  • Обработка данных на РАТАН-600
• Информация из Internet
• Схема получения астрофизического результата. Обзор проблем
• Литература
• Приложение

Введение

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

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

Формула эксперимента

Чтобы быть полностью корректным необходимо учесть нормировку на интеграл под аппаратной функцией (который равен 1 из условия термодинамического равновесия между антенной и средой в случае, когда яркостная температура постоянна в окружающем антенну пространстве). Кроме того, на выходе системы "антенна-радиометр" мы наблюдаем случайные колебания мощности N(x,y), вызванные нестабильностью работы приемника, антенны и излучением Земли. Тогда полностью наблюдаемое распределения яркости, мы можем представить в таком виде:

(2)

Упрощая запись, получим для описания процесса наблюдения такую формулу (назовем ее формулой эксперимента):

I = B * A + N         (3)

Здесь I - результат наблюдений, B - реальное распределение яркости объекта, A - аппаратная функция: прибор + глаз, N - многокомпонентный шум N = N_sky * A + N_sys, где N_sky - шум, связанный с прохождением сигнала сквозь галактический газ и земную атмосферу (сглаживается аппаратной функцией), N_sys - шум системы телескоп+приемник+глаз. Символ '*' (звездочка) обозначает операцию свертки. Процедура свертки изображения с аппаратной функцией является процедурой сглаживания изображения, т.е. подавления вариаций потока на масштабах, меньших, чем собственный интервал (грубо говоря, размер) аппаратной функции.

Функции в подынтегральном выражении имеют аргументы противоположных знаков, или, говоря по-другому, при свертке функции ``наползают'' друг на друга, причем в области их пересечения происходит интегрирование. Значение этого интеграла есть текущая точка свертки на данном шаге сглаживания. Сказанное демонстрируется рис.1. Нормировка на интеграл под аппаратной функцией позволяет нам сохранить реальную амплитуду сигнала. На рис.2 приведены 2 примера, демонстрирующие преобразования функций при свертке.

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

Что еще сказать о свертке ? Это коммутативная операция, т.е. интеграл не изменится, если функции поменять местами:

B * A = A * B,         (4)

это аддитивная операция, т.е. свертка суммы двух функций с функцией A, дает свертку каждой из функций с A:

(B1 + B2) * A = B1 * A + B2 * A.        (5)

Свертка с -функцией не меняет саму функцию (см. также рис. 2):

(6)

Рис.2. Схема, демонстрирующая свертку функций различных типов: прямоугольника и гауссианы (1) и прямоугольника и -функции (2).

Пожалуй, о свертке пока хватит.

Таким образом, можно понять, как выглядит объект (т.е. найти B), решив уравнение (3). Но проблема состоит в том, что, кроме свертки, в уравнении присутствует еще и шум N. Именно из-за него возникает проблема с решением (или, точнее говоря, невозможностью точного решения) системы (3). Даже без наличия шума эта задача является некорректной, т.е. характеризуется тем, что малым изменениям входных данных могут соответствовать сколь угодно большие изменения решения. Вообще говоря, практически все задачи наблюдательной астрофизики - задачи некорректные, т.к. по наблюдательным проявлениям объекта требуют установить его характеристики, недоступные прямым измерениям. Но, несмотря на то, что природа столь изощренна в своих запретах и ее, действительно, обмануть нельзя, все-таки используя некоторые подходы можно приблизиться к пониманию того, что же мы видим.

О том, как нам узнать все (или почти все) о функции B, мы поговорим позднее. А сейчас проведем обзор наших средств, способных нам помочь в решении задачи.

Начнем с первого.

Построение алгоритмов

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

Но прежде, чем говорить о возможности решения некорректных задач, для понимания процедур обработки необходимо упомянуть о нескольких разделах математики, которые пригодятся в дальнейшем при объяснении шагов обработки. Это основы статистики и Фурье-анализа.

Параметры распределений

Объединение (логическая сумма) E₁UE₂U... (или E₁+E₂+...) конечной или бесконечной последовательности событий E₁, E₂,..., есть событие, состоящее в осуществлении хотя бы одного из событий E₁, E₂, ...

Совмещение (логическое произведениe) E₁, E₂ E₁ E₂ (или E₁ E₂) двух E₁ и E₂, есть событие, состоящее в осуществлении и E₁ и E₂.

События E₁ и E₂ называются несовместными, если E₁UE₂=0, - логическое произведение есть невозможное событие (т.е. не осуществляется ни одно из событий класса E, к которому принадлежат события E₁ и E₂).

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

Вероятностью события (в статистическом смысле) является предел, к которому стремится относительная частота появления события при неограниченном увеличении числа испытаний (наблюдений) (Агекян, 1974).

Вероятностью P(E) события E называется определенная на классе событий E однозначная действительная функция, удовлетворяющая трем условиям (аксиомам вероятностей):
1) для любого события E из E;
2) P(I) = 1 для достоверного события I;
3) P(E₁UE₁U...) = P(E₁) + P(E₂)+... для любой (конечной или бесконечной) последовательности попарно несовместных событий E₁, E₂.... (Корн & Корн, 1984).

Условная вероятность P(E|E₁) события E при условии осуществления события E₁ определяется аксиомой P(EU E₁) = P(E₁) P(E|E₁) (правило умножения вероятностей).

Теорема Байеса дает выражение для условных вероятностей одного из событий E_k полной группы E₁, ..., E_n (т.е. группы несовместимых событий, одно из которых должно с достоверностью произойти) при условии, что произошло событие A:

(7)

или (см. детали у Агекяна, 1974)

(8)

Эта формула позволяет вычислить ``апостериорные'' вероятности P(E<sub>k</sub>|A) через ``априорные'' вероятности P(E_k|A) ``гипотез'' E_k. Здесь априорная вероятность характеризует предполагаемую возможность осуществления различных значений E_k до того, как проведен эксперимент. Эта предполагаемая возможность основывается на известных фактах, таких, как, например, результаты предыдущих экспериментов. Апостериорная вероятность P(E_k|A) характеризует возможность осуществления различных значений E_k, после того как к априорному знанию добавлено знание, извлеченное из экспериментальных данных A.

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

Апостериорная вероятность (небо | данные) ~ Априорная вероятность (небо) x Априорная~ вероятность (данные | небо)      (9)

Фактически последняя формула - это вероятностное представление уравнения (3), которое также используется для его решения.

Астрономы, говоря о чувствительности аппаратуры и параметрах объектов, используют две основные статистические характеристики выборки: среднее и дисперсию.

Выборка (x₁, x₂, ..., x_n) значений x рассматривается как результат n независимых измерений. Каждую выборку рассматривают как выборку из теоретически бесконечной генеральной совокупности.

Тогда основные параметры выборки - среднее и дисперсия определяются следующим образом. Среднее значение случайной величины x определяется как

(10)

где p(x_i) - распределение вероятности для дискретной последовательности, - плотность вероятности для непрерывной зависимости.

Дисперсия (и связанное с ней стандартное отклонение ) определяется как

(11)

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

Распределения

1. Нормальное (или гауссово) распределение.
   Случайная величина x распределена нормально со средним , и стандартным отклонением , если распределение вероятности подчиняется следующему закону:

   (12)

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

2. -функция.

   (13)

   Эта функция обладает свойством . Параметры распределения, которое описывается этой функцией, такие:

   .

   Это - важное распределение, введенное Дираком, характеризующее, например, положение точечного источника на карте.

3. Распределение Пуассона.

   (14)

   . Очень ``астрономическое'' распределение, описывающее статистику объектов на небе либо число регистрируемых фотонов в приемнике.

4. Распределение Коши.

   (15)

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

5. Равномерное (прямоугольное) распределение.

   (16)

   .

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

Компьютерное моделирование распределений

Псевдоравномерное распределение генерируется, например, функциями языка ``C'' rand() и drand48().

Чтобы получить гауссово распределение . мы должны воспользоваться следующим преобразованием:

(17)

Распределение Коши получается простым преобразованием: arctan(rand()).

Фурье-анализ

Прямое преобразование (в двумерном случае) функции I(x,y) это двумерный интеграл от функций комплексного переменного:

(18)

Обратное преобразование Фурье от V(u,v) выглядит похожим образом:

(19)

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

V = F(I),       (18a)

а обратное как

I = F^-1(V).        (19a)

I мы будем называть изображением, а F(I) - его Фурье-образом.

Преобразование Фурье чрезвычайно важная операция для описания и понимания процессов, происходящих в природе. Впервые оно, по-видимому, было использовано астрономом, который об этом даже и не догадывался. Птолемей, живший в 90-160-х годах нашей эры, предложил систему движения планет, обращающихся вокруг Земли по эпициклам. Математически это соответствовало разложению эллиптического движения по круговым функциям, т.е. в ряд Фурье.

Преобразование Фурье позволяет по-другому посмотреть на окружающий мир, где все дискретные предметы ``превращаются'' в набор пространственных волн. Астрономы для описания угловых масштабов источника при наблюдении на радиотелескопах тоже пользуются терминологией Фурье-преобразования. Область наблюдения описывается на языке комплексной (u,v)-плоскости, получаемой из действительной координатной плоскости преобразованием (18). Причем, качество получаемого изображения зависит от качества заполнения (u,v)-плоскости, т.е. наличия соответствующих пространственных частот при синтезе изображения (для более подробного изучения проблемы рекомендую монографию ``Радиотелескопы и радиометры'' Есепкиной и др., 1973).

``В быту'' астрономы используют модуль преобразования Фурье, называемый спектром мощности, и фазу, вычисляемую как -arctg(v/u). На рис. 3 приведены примеры спектров некоторых функций.

Рис.3. Примеры спектров мощности некоторых функций.

Для расчетов обычно используется так называемое быстрое преобразование Фурье (БПФ), основанное на малом количестве операций для массивов, чья размерность равна 2ⁿ. Существует несколько алгоритмов вычисления БПФ, среди которых один из наиболее распространенных - алгоритм Кули-Тьюки (текст программы на ФОРТРАНе можно найти в монографии Гончарского и др., 1985).

Теорема о свертке

(20)

На практике это дает значительное уменьшение количества операций при вычислении свертки, с одной стороны, и показывает возможность решения уравнения (3) при отсутствии N, с другой стороны.

Теорема отсчетов

Справедливо и обратное утверждение:
если функция существует только в промежутке (-x₀/2,   x₀/2), т.е. ограничена в пространстве, то для полного определения ее спектра достаточно измерять компоненты спектра, отстоящие друг от на интервал частот .

Приведем пример использования теоремы в радиоастрономии (Есепкина и др., 1973):
"Пространственно-частотная характеристика антенн ограничена по частоте (F_гр=D/, D - максимальный размер антенны), следовательно, измеряемые распределения яркости I_a(x) должны иметь ограниченный спектр. Несущие информацию отсчеты в функции I_a(x) расположены на интервалах, не меньших, чем . Этот интервал определяет разрешающую способность радиотелескопа; он получил название собственного интервала антенны."

Ниже мы еще вернемся к этой теореме.

Рис.4. Различные способы анализа данных (рисунок воспроизведен из книги Хампеля и др. (1986)).

Рис.5. Какого соответствия мы хотим ? (рисунок воспроизведен из книги Хампеля и др. (1986)).

Робастность

Из многих методов, именно робастные оценки, т.е. устойчивые к элементам влияния, приобрели особую значимость при анализе астрономических изображений, в том числе и в системе редукции данных на РАТАН-600.

История робастности в обработке данных началась еще в середине 30-х годов XX века, когда пришлось иметь дело с выпадающими данными и разрабатывать статистическую теорию устойчивых оценок. Революционные изменения произошли с широким распространением компьютеров в 70-е и 80-е годы, когда были реализованы на практике многие теоретические методы (Хьюбер, 1981; Хампель и др., 1986) На рис. 4 показано качественное сравнение методов анализа и устойчивости получаемого результата. На рис. 5 приведен пример робастного подхода к регрессионным оценкам.

Наиболее простые из робастных методов построены на ранговых алгоритмах. Рассмотрим возможности оценок параметров (моментов) распределения при гауссовом приближении статистики шумового сигнала (Шергин и др., 1997) для робастного подхода.

• Оценка среднего. Медиана.
  

  Рис.6. Пример работы простого рангового алгоритма оценки среднего - вычисление медианы путем сортировки данных: слева - до сортировки (найдите среднее), справа - после сортировки (вот оно !).

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

• Метод Ходжеса-Леманна
  Однако, не всегда выборка настолько велика, что медианная оценка работает устойчиво. Для усиления устойчивости часто используют размножение элементов выборки, а новый набор используют для оценки среднего. Один из таких методов - Ходжеса-Леманна (Хьюбер, 1981) - основан на разможении выборки путем вычисления всевозможных полусумм элементов из исходной выборки:

  (21)

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

• Оценка дисперсии (если мы предполагаем нормальный шум).
  Для оценки дисперсии гауссовых шумов используем момент второго порядка: медианное абсолютное отклонение, которое пропорционально дисперсии:

  (22)

  Абсолютное медианное отклонение вычисля