Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://sed.sao.ru/~vo/rad_prob/rad_problems.html
Дата изменения: Tue Jun 5 19:18:10 2007 Дата индексирования: Mon Oct 1 20:40:29 2012 Кодировка: koi8-r Поисковые слова: m 15 |
Лекция для студентов (полная gzipped PostScript версия (950k) здесь)
Верходанов О.В.
Специальная Астрофизическая обсерватория, 1999 г.
Опыт последних поколений молодых людей, приходящих в астрономические институты, показывает, что первое, с чем им приходится сталкиваться, - это проблемы, связанные с техническим и программным обеспечением и обработкой полученных данных. От уровня реализации этих задач часто зависит успех молодого ученого. К счастью, многие проблемы автоматизации эксперимента на уровне программного обеспечения сейчас уже решены. И то, что занимало порой несколько недель (например, поиск необходимой статьи), теперь может быть решено за несколько минут.
Я освещу несколько моментов, являющихся, на мой взгляд, важными для понимания того, что же мы видим на небе и как получается конечный результат. Данная лекция не претендует на полноту информации, а только очерчивает круг проблем, с которыми приходится сталкиваться радио- (и не только) астрономам. В данной лекции я остановлюсь на проблемах обработки континуальных данных, т.е. полученных на широкополосных радиометрах, а когда буду говорить об изображении, то буду иметь в виду распределения яркости на небе (двумерное и одномерное).
Как и оптический телескоп, в радиоастрономии радиотелескоп играет роль линейного оператора, воздействующего на функцию, описывающую объект исследования. Каждая точка в распределении радиояркости по исследуемому объекту изображается в фокусе радиотелескопа в виде дифракционного пятна или пятна рассеяния (подход описан у Есепкиной и др., 1974). Используя только принцип суперпозиции этих пятен от различных точек объекта можно написать следующее фундаментальное уравнение антенного сглаживания:
(1) |
Где B(x,y) - распределение радиояркости по объекту, A(x',y') - диаграмма направленности радиотелескопа, I(x,y) - наблюдаемое распределение мощности в фокальной плоскости радиотелескопа, создаваемое источником излучения.
Чтобы быть полностью корректным необходимо учесть нормировку на интеграл под аппаратной функцией (который равен 1 из условия термодинамического равновесия между антенной и средой в случае, когда яркостная температура постоянна в окружающем антенну пространстве). Кроме того, на выходе системы "антенна-радиометр" мы наблюдаем случайные колебания мощности N(x,y), вызванные нестабильностью работы приемника, антенны и излучением Земли. Тогда полностью наблюдаемое распределения яркости, мы можем представить в таком виде:
(2) |
I = B * A + N (3)
Здесь I - результат наблюдений, B - реальное распределение яркости объекта, A - аппаратная функция: прибор + глаз, N - многокомпонентный шум N = Nsky * A + Nsys, где Nsky - шум, связанный с прохождением сигнала сквозь галактический газ и земную атмосферу (сглаживается аппаратной функцией), Nsys - шум системы телескоп+приемник+глаз. Символ '*' (звездочка) обозначает операцию свертки. Процедура свертки изображения с аппаратной функцией является процедурой сглаживания изображения, т.е. подавления вариаций потока на масштабах, меньших, чем собственный интервал (грубо говоря, размер) аппаратной функции.
Функции в подынтегральном выражении имеют аргументы противоположных знаков, или, говоря по-другому, при свертке функции ``наползают'' друг на друга, причем в области их пересечения происходит интегрирование. Значение этого интеграла есть текущая точка свертки на данном шаге сглаживания. Сказанное демонстрируется рис.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)) способствуют развитию методов в различных направлениях - науки от геофизики до медицины.
Но прежде, чем говорить о возможности решения некорректных задач, для понимания процедур обработки необходимо упомянуть о нескольких разделах математики, которые пригодятся в дальнейшем при объяснении шагов обработки. Это основы статистики и Фурье-анализа.
Объединение (логическая сумма) E1UE2U... (или E1+E2+...) конечной или бесконечной последовательности событий E1, E2,..., есть событие, состоящее в осуществлении хотя бы одного из событий E1, E2, ...
Совмещение (логическое произведениe) E1, E2 E1 E2 (или E1 E2) двух E1 и E2, есть событие, состоящее в осуществлении и E1 и E2.
События E1 и E2 называются несовместными, если E1UE2=0, - логическое произведение есть невозможное событие (т.е. не осуществляется ни одно из событий класса E, к которому принадлежат события E1 и E2).
Существует много определений вероятности. Упустив рассмотрение классического определения (как отношения числа благоприятных событий к числу равновозможных), остановимся на двух. Первое определение показывает смысл вероятности, а второе характеризует ее как функцию.
Вероятностью события (в статистическом смысле) является предел, к которому стремится относительная частота появления события при неограниченном увеличении числа испытаний (наблюдений) (Агекян, 1974).
Вероятностью P(E) события E
называется определенная на классе событий
E однозначная действительная функция, удовлетворяющая трем условиям
(аксиомам вероятностей):
1)
для любого события E из E;
2) P(I) = 1 для достоверного события I;
3) P(E1UE1U...) =
P(E1) + P(E2)+... для любой
(конечной или бесконечной) последовательности попарно несовместных событий
E1, E2.... (Корн & Корн, 1984).
Условная вероятность P(E|E1) события E при условии осуществления события E1 определяется аксиомой P(EU E1) = P(E1) P(E|E1) (правило умножения вероятностей).
Теорема Байеса дает выражение для условных вероятностей одного из событий Ek полной группы E1, ..., En (т.е. группы несовместимых событий, одно из которых должно с достоверностью произойти) при условии, что произошло событие A:
(7) |
или (см. детали у Агекяна, 1974)
(8) |
Эта формула позволяет вычислить ``апостериорные'' вероятности P(Ek|A) через ``априорные'' вероятности P(Ek|A) ``гипотез'' Ek. Здесь априорная вероятность характеризует предполагаемую возможность осуществления различных значений Ek до того, как проведен эксперимент. Эта предполагаемая возможность основывается на известных фактах, таких, как, например, результаты предыдущих экспериментов. Апостериорная вероятность P(Ek|A) характеризует возможность осуществления различных значений Ek, после того как к априорному знанию добавлено знание, извлеченное из экспериментальных данных A.
В нашем случае все знания о радионебе можно описать относительными
вероятностями для различных карт неба. Тогда теорема Байеса говорит
о том, как следует видоизменить наши знания или предварительные сведения
о небе, если будут получены некоторые дополнительные данные, и
может может быть записана так:
Апостериорная вероятность (небо | данные) ~ Априорная вероятность (небо) x Априорная~ вероятность (данные | небо) (9)
Фактически последняя формула - это вероятностное представление уравнения (3), которое также используется для его решения.
Астрономы, говоря о чувствительности аппаратуры и параметрах объектов, используют две основные статистические характеристики выборки: среднее и дисперсию.
Выборка (x1, x2, ..., xn) значений x рассматривается как результат n независимых измерений. Каждую выборку рассматривают как выборку из теоретически бесконечной генеральной совокупности.
Тогда основные параметры выборки - среднее и дисперсия определяются следующим образом. Среднее значение случайной величины x определяется как
(10) |
где p(xi) - распределение вероятности для дискретной последовательности, - плотность вероятности для непрерывной зависимости.
Дисперсия (и связанное с ней стандартное отклонение ) определяется как
(11) |
Вспомним, как ведут себя эти параметры для различных распределений случайной величины.
(12) |
Это, пожалуй, самое распространенное распределение в природе, которому подчиняются и распределение молекул по скоростям в газе, и шум радиометра и многие другие случайные события.
(13) |
Эта функция обладает свойством . Параметры распределения, которое описывается этой функцией, такие:
.
Это - важное распределение, введенное Дираком, характеризующее, например, положение точечного источника на карте.
(14) |
. Очень ``астрономическое'' распределение, описывающее статистику объектов на небе либо число регистрируемых фотонов в приемнике.
(15) |
Mx и Dx не существуют. Это "естественное" распределение в природе описывает случай, когда мы не можем предсказать, что же произойдет в следующий момент, например, пролет частицы или метеора через поле зрения телескопа или переговоры по мобильному телефону, дающие помехи в дециметровом диапазоне радиотелескопа.
(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. Примеры спектров мощности некоторых функций. |
Для расчетов обычно используется так называемое быстрое преобразование Фурье (БПФ), основанное на малом количестве операций для массивов, чья размерность равна 2n. Существует несколько алгоритмов вычисления БПФ, среди которых один из наиболее распространенных - алгоритм Кули-Тьюки (текст программы на ФОРТРАНе можно найти в монографии Гончарского и др., 1985).
(20) |
На практике это дает значительное уменьшение количества операций при вычислении свертки, с одной стороны, и показывает возможность решения уравнения (3) при отсутствии N, с другой стороны.
Справедливо и обратное утверждение:
если функция существует только в промежутке
(-x0/2, x0/2), т.е.
ограничена в пространстве, то для полного определения ее спектра достаточно
измерять компоненты спектра, отстоящие друг от на интервал частот
.
Приведем пример использования теоремы в радиоастрономии (Есепкина и др.,
1973):
"Пространственно-частотная характеристика антенн ограничена по частоте
(Fгр=D/,
D - максимальный размер антенны), следовательно,
измеряемые распределения яркости
Ia(x) должны иметь ограниченный спектр.
Несущие информацию отсчеты в функции
Ia(x) расположены на интервалах,
не меньших, чем
.
Этот интервал определяет
разрешающую способность радиотелескопа; он получил название
собственного интервала антенны."
Ниже мы еще вернемся к этой теореме.
Рис.4. Различные способы анализа данных (рисунок воспроизведен из книги Хампеля и др. (1986)). |
Рис.5. Какого соответствия мы хотим ? (рисунок воспроизведен из книги Хампеля и др. (1986)). |
Из многих методов, именно робастные оценки, т.е. устойчивые к элементам влияния, приобрели особую значимость при анализе астрономических изображений, в том числе и в системе редукции данных на РАТАН-600.
История робастности в обработке данных началась еще в середине 30-х годов XX века, когда пришлось иметь дело с выпадающими данными и разрабатывать статистическую теорию устойчивых оценок. Революционные изменения произошли с широким распространением компьютеров в 70-е и 80-е годы, когда были реализованы на практике многие теоретические методы (Хьюбер, 1981; Хампель и др., 1986) На рис. 4 показано качественное сравнение методов анализа и устойчивости получаемого результата. На рис. 5 приведен пример робастного подхода к регрессионным оценкам.
Наиболее простые из робастных методов построены на ранговых алгоритмах. Рассмотрим возможности оценок параметров (моментов) распределения при гауссовом приближении статистики шумового сигнала (Шергин и др., 1997) для робастного подхода.
Рис.6. Пример работы простого рангового алгоритма оценки среднего - вычисление медианы путем сортировки данных: слева - до сортировки (найдите среднее), справа - после сортировки (вот оно !). |
Пожалуй, самый распространенный из робастных методов - это простейший ранговый алгоритм с использованием медианы как оценки среднего (рис.6). Простота медианной оценки среднего привела к тому, что медианное осреднение и сглаживание включены во все стандартные системы обработки. Медиана замечательна еще и тем, что в качестве среднего используется реально наблюдаемая величина.
(21) |
Данный алгоритм оценки работает успешно и устойчиво уже для массива из пяти элементов.
(22) |
Абсолютное медианное отклонение вычисляется подобно медиане: ранжируются абсолютные значения отклонений от среднего . и берется элемент с номером n/2. Покажем некоторые стандартные способы использования робастных статистик.
Рис.7. Пример сглаживания записи источника с использованием оценки Ходжеса-Леманна в интервалах различной длины. |
Рис.8. Пример робастного осреднения примерно 20 записей, часть из которых показана на рисунке. По нулевому уровню проходит запись, полученная в результате осреднения. |