Анализ двумерных астрономических изображений методами Фурье и вейвлет анализа

Костенков Александр Евгеньевич

Анализ двумерных астрономических изображений методами Фурье и вейвлет анализа

В работе рассмотрено применение методов Фурье и вейвлет-анализа в рамках многоуровневых моделей для обработки двумерных изображений. Представлены алгоритмы фильтрации данных и их адаптация для работы с астрономическими изображениями. Изложено применение многоуровневого подхода для поиска точечных и протяженных объектов на снимках. В ходе исследования были найдены звездные скопления в галактике М33. Методами статистического анализа доказана их связь с LBV-звездами в галактике.

Математика

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36635f1be77c40d58e40

UUID: 76e0b7b0-6b40-4896-832e-0ab0ffcd05bd

Язык: Русский

Опубликовано: больше 7 лет назад

Просмотры: 59

Костенков Александр Евгеньевич

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 8708955 bytes

Поделиться работой

ПРАВИТЕЛЬСТВО РОССИЙСКОЙ ФЕДЕРАЦИИ ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧЕРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ «САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ» КАФЕДРА АСТРОФИЗИКИ Костенков Александр Евгеньевич Анализ двумерных астрономических изображений методами Фурье и вейвлет анализа ДИПЛОМНАЯ РАБОТА Допущен к защите. Заведующий кафедрой астрофизики: д.ф.-м.н., профессор Гаген-Торн В.А. Научный руководитель: д.ф.-м.н., профессор Холтыгин А.Ф. Научный консультант: д.ф.-м.н., профессор С.Н. Фабрика Рецензент: д.ф.-м.н., профессор М.А. Погодин Санкт-Петербург, 2016

SAINT-PETERSBURG STATE UNIVERSITY DEPARTMENT OF ASTROPHYSICS Kostenkov Alexander Analysis of Two-dimensional Astronomical Images by Methods of Fourier and Wavelet Analysis GRADUATION THESIS Admitted for defence. Head of Department: Doctor of Sciences, Professor V.A. Hagen-Thorn Scientific supervisor: Doctor of Sciences, Professor A.F. Kholtygin Scientific supervisor: Doctor of Sciences, Professor S.N. Fabrika Reviewer: Doctor of Sciences, Professor M.A. Pogodin Saint-Petersburg, 2016

Оглавление 1 Основные преобразования и соотношения 1.1 Пространственная фильтрация . . . . . . . . . . . . . . . . . 1.1.1 Сглаживающие пространственные фильтры . . . . . . 2 Частотная фильтрация 7 8 9 11 2.1 Дискретное преобразование Фурье . . . . . . . . . . . . . . . 11 2.2 Вейвлет преобразование . . . . . . . . . . . . . . . . . . . . . 15 2.2.1 Общие понятия вейвлет анализа . . . . . . . . . . . . 16 2.2.2 Кратномасштабный анализ . . . . . . . . . . . . . . . 16 2.2.3 Дискретное вейвлет преобразование . . . . . . . . . . 19 2.2.4 Двумерный КМА . . . . . . . . . . . . . . . . . . . . . 21 2.2.5 Альтернативные алгоритмы вейвлет преобразований 22 3 Trous-вейвлет преобразование 23 3.1 Trous-вейвлет алгоритм . . . . . . . . . . . . . . . . . . . . . 23 3.2 Одномерные многоуровневые преобразования . . . . . . . . . 27 3.3 Двумерные многоуровневые преобразования . . . . . . . . . 31 4 Многоуровневая обработка астрономических снимков 4.1 Модель сигнал-шум . . . . . . . . . . . . . . . . . . . . . . . 35 35 4.1.1 Soft and hard thresholding . . . . . . . . . . . . . . . . 36 4.1.2 3𝜎-clipping и винеровская фильтрация . . . . . . . . . 37 1

5 Многоуровневое представление данных 41 5.1 Описание многоуровневой модели . . . . . . . . . . . . . . . 41 5.2 Поиск объектов и структур . . . . . . . . . . . . . . . . . . . 43 5.3 Алгоритм поиска объектов . . . . . . . . . . . . . . . . . . . 44 5.4 SExtractor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 5.4.1 Поиск рассеянных структур и скоплений . . . . . . . 48 5.4.2 Проблемы выделения звезд и скоплений . . . . . . . . 49 6 Поиск звездных скоплений в М33 52 6.1 Luminous blue variable . . . . . . . . . . . . . . . . . . . . . . 52 6.2 Выбор исходных данных . . . . . . . . . . . . . . . . . . . . . 54 6.3 Определение параметров групп звезд . . . . . . . . . . . . . 55 6.4 Идентификация объектов . . . . . . . . . . . . . . . . . . . . 57 Приложение А 70 Приложение Б 72 Приложение В 74 Приложение Г 76 Приложение Д 77 2

Абстракт В дипломной работе представлены способы обработки двумерных астрономических изображений, основанные на различных методах матричной свертки, фильтрации, а так же двумерном вейвлет анализе. Представлен подход к декомпозиции изображений адаптированный для обработки астрономических снимков. Отладка алгоритмов производилась с использованием модельных решений. Первая часть диплома - описание математического аппарата и существующих способов, изложение теории многоуровневого представления в астрономии. Второй этап работы является демонстрацией использования различных методов анализа изображений на модельных решениях, сравнение результатов, полученных разными способами.Третья часть - применение вейвлет анализа в идентификации скоплений, адаптивный поиск объектов на снимках. Используя этот подход, были выделены скопления в галактике М33, которые могут быть связаны с LBV-звездами . Относительно скоплений было составлено распределение расстояний до звезд. Такой же подход можно применить и для других галактик Местной группы. 3

Введение Обработка и анализ данных становятся все более значимыми в астрономии. Это можно объяснить эволюцией приемников излучения - в наше время размеры ПЗС матриц могут достигать 16000 × 16000 пикселей и более. Следствием этого является огромное количество получаемой информации, которую невозможно обрабатывать вручную. Поэтому нужны методы и алгоритмы, позволяющие автоматически обрабатывать полученные данные. Большую часть данных получают из снимков, о них и пойдет речь в этой работе. Чаще всего изображения сильно зашумлены, объекты в низком разрешении на них сложно идентифицировать, присутствует эффект проекции на снимке, особенно в переполненных звездами полях, что осложняет классификацию. Для того, чтобы извлечь информацию о каждом объекте, применяют методы математической обработки снимков, начиная от простых преобразований, например, матричной свертки, до сложных декомпозиций изображений при помощи вейвлет преобразований. Это позволяет анализировать снимок на разных масштабах, составлять многоуровневое представление объекта - выделить самые мелкие и более крупные образования. Этот подход является очень сильным инструментом в анализе данных. В работе ставилась цель выделить характеристики скоплений, связанных с LBV, путем обработки снимков, идентифицировать аналогичные объекты на всем изображении. Решение этой задачи можно разделить на три этапа: 4

∙ Анализ и сбор данных - первичная редукция ∙ Обработка снимков - математическая обработка,например, очистка от шума, выделение исследуемых объектов, фильтрация изображения ∙ Выделение необходимых критериев для идентификации скоплений, классификация объектов на снимке В данный момент существует большое количество алгоритмов и моделей, которые используются в обработке снимков. Применяют как простые преобразования изображений с помощью фильтров низких и высоких частот (НЧ и ВЧ), так и обработка в частотной области при помощи Фурье и вейвлет анализа. Широко используется многоуровневое представление данных. Отличным примером реализации такого подхода является SExtractor (Source Extractor, см. [3]). Это программное обеспечение, которое выделяет и классифицирует астрономические объекты на снимках. В ней заложен интересный подход к анализу данных: поиск объектов и оценка фона осуществляется после свертки изображения с выбранным фильтров, вариации фильтров позволяют искать объекты разной формы и размеров (например, Point Spread Function (PSF) для поиска одиночных звезд, для протяженных объектов широкий НЧ фильтр). Иногда специфика задачи не предполагает точных критериев - хорошим примером являются поиск молодых скоплений, которые могут занимать довольно обширную площадь, а применение широких фильтров приводит к размыванию и захвату лишних звезд в переполненых звездами областях ("crowded fields"). В данной работе предлагается использовать многоуровневый подход в задаче - декомпозиция изображения на структуры разного уровня, поиск интересующих объектов с помощью набора фильтров на каждом уровне. Так можно найти как крупные, так и мелкие структуры. Проблема поиска скоплений в переполненных областях все еще остается открытой, нашу задачу упрощает тот факт, что производился поиск 5

молодых скоплений рядом с которыми находится эмиссионная туманность (область HII с выcоким излучением в линии 𝐻𝛼 , возраст < 108 -109 лет). 6

Глава 1 Основные преобразования и соотношения Информация, содержащаяся в снимках, может быть представлена разными способами. В данной работе будет рассматриваться пространственная интерпретация. Растровые изображения представляются пространственным распределением освещенности на плоскости. Ее можно описать как функцию двух переменных: 𝑓 = (𝑥, 𝑦) (1.1) где 𝑓 - освещенность клетки с координатами 𝑥, 𝑦. Таким образом изображение 𝐴 может быть представлено в виде двумерных массивов: ⎛ 𝑎 𝑎12 ⎜ 11 ⎜ ⎜ 𝑎21 𝑎22 𝐴=⎜ ⎜ .. .. ⎜ . . ⎝ 𝑎𝑚1 𝑎𝑚2 · · · 𝑎1𝑛 ⎞ ⎟ ⎟ · · · 𝑎2𝑛 ⎟ ⎟ .. ⎟ ... . ⎟ ⎠ · · · 𝑎𝑚𝑛 (1.2) где 𝑎𝑚𝑛 - освещенность пикселя с заданными кординатами. В дипломной работе рассматриваются полутоновые и бинарные изображения. Полутоновое изображение кодируется с помощью битовой сетки, 7

то есть каждый пиксель описывается разным количеством бит, что и определяет количество полутонов. Например, 2 бита - 4 полутона, 3 хранят 8 и т.д. Бинарное изображение - вырожденный случай полутонового (1 бит на 1 пиксель), пиксель передает лишь 2 полутона - черный или белый (в цифровом представлении матрица из 0 и 1). Далее под обработкой снимка будет подразумеваться преобразование матрицы изображения для какой-то конкретной задачи. Большинство таких методов делится на две большие категории: пространственная обработка (например, фильтрация по маске) и частотная обработка (преобразование Фурье и т.д.). Ниже будет рассмотрен математический аппарат тех способов, которые были использованы в работе (описание взято из [33], [31]). 1.1 Пространственная фильтрация Некоторые локальные преобразования оперируют одновременно как со значениями пикселей в окрестности, так и с соотвествующими им значениями некоторой матрицы, имеющий те же размеры, что и окрестность. Соседние клетки определяются стандартно: все комбинации(𝑥 ± 1, 𝑦 ± 1), где (𝑥, 𝑦) координаты выбранного пикселя в матрице. Схема пространственной фильтрации основана на простом перемещении маски фильтра от точки к точке изображения. Обработка изображения 𝑓 , имеющего размеры 𝑀 × 𝑁 , с помощью фильтра размерами 𝑚 × 𝑛 задается выражением общего вида: 𝑔(𝑥, 𝑦) = 𝑎 ∑︁ 𝑏 ∑︁ 𝑤(𝑠, 𝑡)𝑓 (𝑥 + 𝑠, 𝑦 + 𝑡) (1.3) 𝑠=−𝑎 𝑡=−𝑏 где, как следует из предыдущего абзаца, 𝑎 = (𝑚 − 1)/2 и 𝑏 = (𝑛 − 1)/2. При вычислении итоговой интенсивности всего изображения данная формула должны быть вычислена для всех сочетаний 𝑥 = 0, 1, 2 . . . 𝑀 − 1 и 8

𝑦 = 0, 1, 2 . . . 𝑁 −1. Это означает что все элементы изображения будут обработаны по заданной маске. Процедура линейной фильтрации, задаваемая этим уравнением в частотной области аналогична процедуре свертки в Фурье пространстве. Поэтому эту процедуру называет "свертку маски(ядра) с изображением". Важным вопрос при применении формулы 1.3 остается ситуация, когда центр маски находится вблизи границы изображения. Можно ограничить перемещение фильтра по границе изображения либо достраивать изображение с помощью нулей, а потом убрать их. В зависимости от выбора ядра можно сгладить изображение или,наоборот, выделить границы. 1.1.1 Сглаживающие пространственные фильтры Сглаживающие пространственные фильтры применяются для расфокусировки изображения и подавления шума. Расфокусировка может применяться как предварительный шаг для обработки изображений, например для удаления мелких деталей для обнаружения более крупных объектов или для устранения разрывов в линиях или деталях. Примером линейного низкочастного фильтра является размытие по Гауссу. Производится свертка матрицы изображения с ядром вида: −(𝑥2 +𝑦 2 ) 1 𝑔(𝑠, 𝑡) = · 𝑒 2𝜎2 2𝜋𝜎 2 (1.4) где (𝑥, 𝑦) - координаты в маске, 𝜎 - среднеквадратичное отклонение. Медианный фильтр - один из видов нелинейных цифровых фильтров, широко используемый в обработке сигналов изображений. Медианная фильтрация в меньшей степени сглаживает границы изображения,чем любая линейная свертка. Значения отсчётов внутри окна фильтра сортируются в порядке возрастания (убывания) и значение, находящееся в середине упорядоченного списка, поступает на выход фильтра. В случае четного чис9

ла отсчетов в окне выходное значение фильтра равно среднему значению двух отсчетов в середине упорядоченного списка. Окно перемещается вдоль фильтруемого сигнала и вычисления повторяются. Низкочастотные фильтры используются в астрономии для выделения более крупных структур или объектов на снимках (например, галактики или скопления). Они очень просты в использовании, но существуют более тонкие алгоритмы по выделению структур в частотной области (об этом в следующей главе). 10

Глава 2 Частотная фильтрация В данной главе представлена математическая интерпретация методов, основанных на обработке в частотной области - Фурье и вейвлет преобразованиях. Они представляют собой очень гибкие инструменты для фильтрации и анализа данных. Для визуализации примеров в этой главе использовались их реализации среды OCTAVE. Теория изложена в [32], [33]. 2.1 Дискретное преобразование Фурье Дискретное преобразование Фурье переводит конечную последовательность вещественных чисел в конечную последовательность коэффициентов Фурье. Пусть {𝑥𝑖 } , 𝑖 = 0, · · · , 𝑁 −1 - последовательность вещественных чисел - например, отчеты яркости пикселов по строке изображения. Эту последовательность можно представить в виде комбинации конечных сумм вида 𝑁/2 ∑︁ (︂ 2𝜋𝑛𝑖 𝑥 𝑖 = 𝑎0 + 𝑎𝑛 cos 𝑁 𝑛=1 )︂ 𝑁/2 ∑︁ (︂ 2𝜋𝑛𝑖 + 𝑏𝑛 sin 𝑁 𝑛=1 )︂ (2.1) где 𝑁 −1 1 ∑︁ 𝑎0 = 𝑥𝑖 , 𝑁 𝑖=0 𝑎𝑁/2 𝑁 −1 1 ∑︁ = 𝑥𝑖 (−1)𝑖 , 𝑁 𝑖=0 11 (︂ )︂ 𝑁 −1 2 ∑︁ 2𝜋𝑖𝑘 (2.2) 𝑎𝑘 = 𝑥𝑖 cos 𝑁 𝑖=0 𝑁

(︂ )︂ 𝑁 −1 2 ∑︁ 2𝜋𝑖𝑘 𝑏𝑘 = 𝑥𝑖 sin , 𝑁 𝑖=0 𝑁 1 ≤ 𝑘 < 𝑁/2. (2.3) Часто используют комплексную запись, основанную на формуле Эйлера: 𝑁 −1 −2𝜋𝑖𝑚𝑛 1 ∑︁ 𝐺𝑚 = 𝑥𝑛 𝑒 𝑁 𝑁 𝑛=1 (2.4) Дискретное преобразование Фурье для двумерного массива чисел размера 𝑀 × 𝑁 определяется следующим образом: 𝐺𝑢𝑤 𝑁 −1 𝑀 −1 𝑚𝑢 𝑛𝑤 1 ∑︁ ∑︁ = 𝑥𝑚𝑛 𝑒−2𝜋𝑖[ 𝑀 + 𝑁 ] 𝑁 𝑀 𝑛=1 𝑚=1 (2.5) Двумерное преобразование Фурье может быть преставлено в виде композиции двух одномерных преобразований: 𝐺𝑢𝑤 [︃ 𝑀 −1 ]︃ 𝑁 −1 −2𝜋𝑖𝑚𝑤 −2𝜋𝑖𝑛𝑢 1 ∑︁ 1 ∑︁ = 𝑥𝑚𝑛 𝑒 𝑀 𝑒 𝑁 𝑁 𝑛=1 𝑀 𝑚=0 (2.6) Здесь выражение в квадратных скобках есть одномерное преобразование строки матрицы данных. Таким образом, для получения двумерного преобразования Фурье нужно сначала вычислить одномерные преобразования строк, записать результаты в исходную матрицу и вычислить одномерные преобразования для столбцов полученной матрицы. Алгоритмы фильтрации в частотной области основываются на свертке. В двумерном случае преобразование свертки выглядит следующим образом: 𝐺 (𝑢, 𝑣) = 𝐻 (𝑢, 𝑣) 𝐹 (𝑢, 𝑣) (2.7) где 𝐺 - Фурье-образ результата свертки, 𝐻 - Фурье-образ фильтра, а 𝐹 Фурье-образ исходного изображения. 12

На Рис. 2.1 приведен пример обработки снимка в частотной области при помощи НЧ фильтра Гаусса (𝜎 = 3). На изображении представлены два скопления - NGC 869 и NGC 884 (более известные как 𝜒 и ℎ Персея, любительский снимок - SBIG ST-8300M, Astro-Tech AT90DT, f/6.7 603 мм, выдержка 285 минут, Dardenne Prairie, Missouri). Используя широкий фильтр можно убрать высокочастотные шумы для того, чтобы идентифицировать скопления. Можно выбрать маску на границе обработанного изображения и бинаризовать по порогу, чтобы определить точное положение. Чем уже профиль фильтра в частотной области (чем больше 𝜎), тем он шире в пространственной. а) б) Рис. 2.1: а) Исходное изображение б) Результат свертки с НЧ фильтром 13

С помощью аналогичных преобразований можно выделить ВЧ (например, фон, см. Рис. 2.2 ) Комбинируя разные маски 𝐻(𝑢, 𝑣) можно выделить почти любую часть изображения. Легко подобрать фильтр индивидуально, т.к. в частотном пространстве двумерная свертка заменяется поэлементным перемножением образов исходного изображения и соответствующего ядра. Рис. 2.2: Вклад ВЧ в исходное изображение С позиции точного представления произвольных сигналов преобразование Фурье имеет ряд ограничений и недостатков. Оно не учитывает колебания сигнала во времени.Локальные особенности сигнала (разрывы, ступеньки, пики и т.п.) дают едва заметные составляющие спектра, по которым обнаружить эти особенности, и тем более их место и характер, практически невозможно. Эти ограничения можно избежать, используя оконное преобразование Фурье, однако появляется ограничение по времени и частоте для всех точек, которое не может быть адаптированно к локальным свойствам сигнала. 14

Рис. 2.3: Мощность спектра исходного изображения, логарифмическая шкала Такие методы, как преобразование Фурье, дискретное косинус преобразование и т.д., позволяют восстанавить исходный сигнал по определенной частотной полосе, большая часть информации о нем содержится в нескольких коэффицентах(спектр мощности 2.3). Это существенный недостаток. Многоуровневые вейвлет-преобразования позволяют фильтровать изображения по всех масштабам, независимо друг от друга. Благодаря выявлению локальных особенностей сигнала, принципиально отсутствующему у преобразования Фурье, вейвлет преобразование нашло широкое применение для анализа тонкой структуры сигналов и изображений, для их сжатия и очистки от шума, что важно и полезно в радиотехнике, электронике, гидроакустике, геофизике, медицине и других областях науки и техники. 2.2 Вейвлет преобразование В этом разделе будет изложен математический аппарат, который в дальнейшем понадобится для описание многоуровневого представления (Multiscale Vision Model) данных. А так же основы вейвлет-анализа и основные формулы для дискретного преобразования. 15

2.2.1 Общие понятия вейвлет анализа Вейвлет преобразованием сигнала будем называть его представление в виде обощенного ряда по системе базисных функций вида: 1 𝜓𝑎𝑏 (𝑡) = √ 𝜓 𝑎 (︂ 𝑡−𝑏 𝑎 )︂ (2.8) сконструированных из материнского (исходного) вейвлета 𝜓(t) , обладающего определенными свойствами за счет операций сдвига во времени (𝑏) и √ изменения временного масштаба (𝑎). Множитель 1/ 𝑎 обеспечивает независимость нормы этих функций от масштабирующего числa. При непрерывном изменении параметров 𝑎 и 𝑏 для расчета вейвлетспектра необходимы большие вычислительные затраты. Множество функций 𝜓𝑎𝑏 (𝑡) избыточно. Необходима дискретизация этих параметров при сохранении возможности восстановления сигнала из его преобразования. Дискретизация, как правило, осуществляется через степень двойки: 𝜓𝑚𝑘 1 =√ 𝜓 𝑎 (︂ 𝑡−𝑏 𝑎 )︂ 1 = √ 𝜓(2−𝑚 𝑛 − 𝑘) 2𝑚 (2.9) где 𝑎 = 2𝑚 , 𝑏 = 𝑘 · 2𝑚 , 𝑚 и 𝑘 - целые числа. В этом случае плоскость 𝑎, 𝑏 превращается в соответствующую сетку 𝑚, 𝑘. Параметр 𝑚 называется параметром масштаба, 𝑛 - это сдвиг. 2.2.2 Кратномасштабный анализ При анализе изображения зачастую удобно представить его в виде совокупности последовательных приближений (Описание взято из [29]). Кратномасштабный анализ (КМА) обладает целым рядом полезных свойств, главным из которых является возможность выделения из исходного сигнала его деталей различных масштабов. Это позволяет адапитровать изображение к исходной задаче. 16

Кратномасштабным анализом называется описание пространства 𝐿2 (𝑅) через иерархически вложенные подпространства 𝑉𝑚 , которые не пересекаются (ортогональность по условию задания), и объединение которых дает в пределе всё 𝐿2 (𝑅), т.е. . . . ⊂ 𝑉2 ⊂ 𝑉1 ⊂ 𝑉0 ⊂ 𝑉−1 ⊂ 𝑉−2 ⊂ . . . , ⋂︁ ⋃︁ 𝑉𝑚 = {0} , 𝑚∈𝑍 = 𝐿2 (𝑅) (2.10) 𝑚∈𝑍 Далее, эти пространства имеют следующее свойство: для любой функции 𝑓 (𝑥) ∈ 𝑉𝑚 , ее сжатая версия будет пренадлежать 𝑉𝑚−1 , 𝑓 (𝑥) ∈ 𝑉𝑚 ←→ 𝑓 (2𝑥) ∈ 𝑉𝑚−1 (2.11) А так же существует такая функция 𝜑(𝑥) ∈ 𝑉0 , что ее сдвиги 𝜑0,𝑛 (𝑥) = 𝜑(𝑥 − 𝑛), 𝑛 ∈ 𝑍 образуют ортонормированный базис пространства 𝑉0 . Из последнего свойства вытекает, что функции 𝜑𝑚,𝑛 (𝑥) = 2−𝑚/2 , 𝜑(2−𝑚 𝑥− 𝑛) образуют ортонормированный базис пространства 𝑉𝑚 . Эти базисные функции называются масштабирующими. Таким образом, любая функция 𝑓 (𝑥) ∈ 𝐿2 (𝑅) может быть представлена множеством последовательных ее приближений 𝑓𝑚 (𝑥) ∈ 𝑉𝑚 . Другими словами, функция 𝑓 (𝑥) есть предел апроксимаций 𝑓𝑚 (𝑥) ∈ 𝑉𝑚 : 𝑓 (𝑥) = lim 𝑓𝑚 (𝑥) (2.12) 𝑥→∞ Из определения кратномасштабного анализа следует, что функция 𝑓𝑚 (𝑥) есть ортогональная проекция 𝑓 (𝑥) на подпространство 𝑉𝑚 , то есть 𝑓 (𝑥) = ∑︁ ⟨𝜑𝑚,𝑛 , 𝑓 (𝑥)⟩𝜑𝑚,𝑛 (𝑥) = 𝑛 ∑︁ 𝑛 17 𝑐𝑚,𝑛 𝜑𝑚,𝑛 (𝑥) (2.13)

Так как 𝜑(𝑥) = 𝜑0,0 (𝑥) ∈ 𝑉0 ⊂ 𝑉−1 , то ∑︁ √ ∑︁ ℎ𝑛 𝜑(2𝑥 − 𝑛) ℎ𝑛 𝜑−1,𝑛 (𝑥) = 2 𝜑0,0 (𝑥) = 2 (2.14) 𝑛 𝑛 где ℎ𝑛 - некоторая последовательность. Это равенство принято называть масштабирующим уравнением. Имеют место следующие соотношения: ∑︁ ℎ𝑛 = 1 (2.15) 𝑛 𝛿0,𝑘 = ⟨𝜑0,0 (𝑥)𝜑0,𝑘 (𝑥)⟩ = 2 ∑︁ ℎ𝑛 ℎ𝑛+2𝑘 (2.16) 𝑛 и в спектральной области |𝐻(𝜔)|2 +|𝐻(𝜔+𝜋)| = 1, что приводит к условиям 𝐻(0) = 1 и 𝐻(𝜋) = 0. Введем в рассмотрение еще один объект - пространства 𝑊𝑚 : 𝑉𝑚−1 = 𝑉𝑚 ⊕ 𝑊𝑚 , ⋂︁ 𝑊𝑚 = {0} , 𝑚∈𝑍 ⋃︁ 𝑊𝑚 = 𝐿2 (𝑅) (2.17) 𝑚∈𝑍 Пусть 𝜓(𝑥) = 𝜓0,0 (𝑥) - базисная функция 𝑊0 , тогда √ ∑︁ 𝜓0,0 (𝑥) = 2 𝑔𝑛 𝜑−1,𝑛 (𝑥) (2.18) 𝑛 для некоторой последовательности 𝑔𝑛 . Определим семейство вейвлет-функций 𝜓𝑚,𝑛 (𝑥) = 2−𝑚/2 𝜓(2−𝑚 𝑥 − 𝑛). ∑︀ Имеет место соотношение: 0 = ⟨𝜑0,0 , 𝜓0,0 ⟩ = 2 ℎ𝑛 𝑔𝑛 . Легко видеть, что 𝑛 выбор 𝑔𝑛 = (−1)𝑛 = ℎ−𝑛+2𝑡+1 корректен для всех целых 𝑡. Определения вейвлет-функций позволяют записать любую функцию ∞ ∑︀ 2 𝑓 (𝑥) ∈ 𝐿 (𝑅) в виде суммы ее проекций на 𝑊𝑗 , 𝑗 ∈ 𝑅: 𝑓 (𝑥) = 𝑒𝑗 (𝑥) 𝑗=−∞ Если осуществить анализ до некоторого масштаба 𝑚, то 𝑓𝑚 (𝑥) будет 18

представлена суммой ее грубой апроксмации 𝑓𝑚 (𝑥) ∈ 𝑉𝑚 и множества деталей 𝑒𝑗 ∈ 𝑊𝑗 : 𝑓 (𝑥) = 𝑓0 (𝑥) + 𝑚 ∑︁ 𝑒𝑗 (𝑥) = ∑︁ 𝑐𝑚,𝑛 𝜑𝑚,𝑛 (𝑥) + 𝑛 𝑗=−∞ 2.2.3 𝑚 ∑︁ 𝑑𝑗,𝑘 𝜓𝑗,𝑘 (𝑥) (2.19) 𝑗=−∞ Дискретное вейвлет преобразование Пусть имеется некоторый дискретный сигнал 𝑐𝑛 . Интерпретируем его как коэффиценты разложения некоторой функции 𝑓0 (𝑥) ∈ 𝑉0 по базису масштабирующих функций пространства 𝑉0 : 𝑓0 (𝑥) = ∑︁ 𝑐0 𝜑0,𝑛 (𝑥) (2.20) 𝑛 где 𝑐0,𝑛 = 𝑐𝑛 . Тогда мы можем вычислить апроксимации этой функции, принадлежащие пространствам 𝑉1 , 𝑉2 , . . .. Согласно идее кратномасштабного анализа функция 𝑓0 (𝑥) раскладывается на сумму: 𝑓0 (𝑥) = 𝑓1 (𝑥) + 𝑒1 (𝑥) = ∑︁ 𝑐1,𝑘 𝜑1,𝑘 (𝑥) + 𝑘 ∑︁ 𝑑1,𝑘 𝜓1,𝑘 (𝑥) (2.21) 𝑘 Таким образом мы получили две новые последовательности коэффицентов 𝑐1,𝑛 и 𝑑1,𝑛 . Этот процесс может быть продолжен разложением 𝑓1 (𝑥), 𝑓2 (𝑥), . . ., и функция 𝑓0 (𝑥) (а следовательно и исходный дискретный сигнал 𝑐𝑛 ) будет представлена совокупностью коэффицентов 𝑑𝑚,𝑛 , 𝑚 ∈ 𝑍 + , 𝑛 ∈ 𝑍. Так определяются дискретные ряды вейвлетов. Учитывая, что масштабирующие функции образуют базисы соотвествующих пространств, можно показать, что: 𝑐𝑗,𝑘 = √ ∑︁ √ ∑︁ 2 𝑐𝑗−1,𝑛 ℎ𝑛+2𝑘 и 𝑑𝑗,𝑘 = 2 𝑐𝑗−1,𝑛 𝑔𝑛+2𝑘 𝑛 𝑛 А обратный преобразование: 19 (2.22)

𝑐𝑗−1,𝑛 = ∑︁ 𝑐𝑗,𝑘 ℎ𝑛+2𝑘 + ∑︁ 𝑑𝑗,𝑘 𝑔𝑛+2𝑘 (2.23) 𝑘 𝑘 Последовательности ℎ𝑛 и 𝑔𝑛 называются фильтрами. Они связаны с базисными функциями 𝜓(𝑥) и 𝜑(𝑥): 𝜑(𝑥) = 2 ∑︁ ℎ𝑘 𝜑(2𝑥 − 𝑘) (2.24) 𝑘 𝜓(𝑥) = 2 ∑︁ (−1)𝑘 ℎ1−𝑘 𝜑(2𝑥 − 𝑘) = 2 ∑︁ 𝑔𝑘 𝜑(2𝑥 − 𝑘) (2.25) 𝑘 𝑘 𝑔𝑘 = (−1)𝑘 ℎ2𝑛−1−𝑘 (2.26) где 𝑘 = 0, 1, 2 . . . 2𝑛−1, 𝑛 - порядок вейвлета. Например, для вейвлетов Добеши (подробнее в [27]), решая уравнение 2.25, можно получить следующий фильтр (для 4 точек): (︃ ℎ= √ √ √ √ )︃ 1+ 3 3+ 3 3− 3 1− 3 √ , √ , √ , √ 4 2 4 2 4 2 4 2 (2.27) Для ВЧ фильтра можно получить соотвеncnвенно 𝑔0 = ℎ3 , 𝑔1 = −ℎ2 , 𝑔2 = ℎ1 , 𝑔3 = −ℎ0 . Так как фильтры пропускают только половину всех частотных компонентов сигнала, то не попавшие в полосу прозрачности составляющие могут быть удалены. Поэтому выполняется децимация ↓ 2, т.е. прореживание в два раза. На стадии реконструкции сигнала осуществляется операция интерполяции ↑ 2 , обратная децимации ↓ 2, путем увеличения в два раза числа составляющих добавлением нулевых компонентов вперемежку с имеющимися. При сложении сигналов, полученных на выходе фильтров, будем иметь сигнал, близкий к исходному, т.е. произойдет его реконструкция на 20

начальном уровне. 2.2.4 Двумерный КМА На практике двумерный КМА строят как тензорное произведение одномерных КМА. При таком подходе отцовский и материнский вейвлеты будут сформированы следующим образом: 𝜑(𝑥, 𝑦) = 𝜑(𝑥)𝜑(𝑦) (2.28) 𝜓𝐿𝐻 (𝑥, 𝑦) = 𝜑(𝑥)𝜓(𝑥) и 𝜓𝐻𝐿 (𝑥, 𝑦) = 𝜓(𝑥)𝜑(𝑦) (2.29) 𝜓𝐻𝐻 (𝑥, 𝑦) = 𝜓(𝑥)𝜓(𝑦) (2.30) где индексы 𝐻 и 𝐿 означают реализацию фильтров ВЧ и НЧ составляющих. Предположим, что мы имеем изображение размером 𝑁 × 𝑁 . первоначально каждая строчка фильтруется на НЧ и ВЧ половины. В результате получаются два изображения размером 𝑁 × 𝑁/2. Далее каждый столбец делится таким же образом. В итоге получаются четыре изображения размером 𝑁/2 × 𝑁/2. НЧ по горизонтали и вертикали (НЧНЧ1 ), ВЧ по горизонтали и вертикали (ВЧВЧ1 ), НЧ по горизонтали и ВЧ по вертикали (НЧВЧ1 ) и ВЧ по горизонтали и НЧ по вертикали (ВЧНЧ1 ). Первое из указанных выше изображений делится аналогичным образом на следующем шаге (уровне) преобразования. Это называется пирамидальным алгоритмом Маала. На Рис. 2.4 показан пример разложения изображения скоплений NGC 869 и NGC 884 на коэффиценты 1 уровня с помощью вейвлета Хаара (для удобства вейвлет коэффиценты нормированы по яркости, выполнено в среде OCTAVE, [35]). 21

a) b) c) d) Рис. 2.4: Коэффиценты пирамидального алгоритма Маала: a) Апроксимация исходного изображения 1-ого уровня (НЧНЧ1 ) b) НЧВЧ1 c) ВЧНЧ1 d) ВЧВЧ1 2.2.5 Альтернативные алгоритмы вейвлет преобразований С практической точки зрения существуют более адаптированные для астрономических задач методы, чем алгоритм Маала. Например, trous-вейвлет преобразование (название связано с особым способом свертки). Оно является многоуровневым: на выходе N изображений такого же размера как и исходное (N задано заранее), причем каждое из них отображает уровень исходного изображения - от высокочастотных составляющих до самых низкочастотных, суммируя все уровни получаем исходный снимок. Это позволяет фильтровать данные на всех уровнях, а затем просто складывать все уровни преобразования. Подробнее о реализации и тестировании будет рассказано в следующей главе. 22

Глава 3 Trous-вейвлет преобразование Для реализации алгоритмов, представленных в этой главе, была выбрана среда OCTAVE. Однако, использованы только стандартные алгебраические библиотеки, поэтому в случае необходимости методы можно легко перенести на другой язык (например, Python+Numpy). 3.1 Trous-вейвлет алгоритм Основной идеей trous-вейвлет преобразование является отсутствие децимации, то есть пропуская сигнал через НЧ и ВЧ фильтры, мы не прореживаем его, а заполняем "дырки"(отсюда и название с французского "algorithme a trous"). В работе был использован алгоритм, изложенный в [11]. На выходе у данного метода 𝐽 изображений одинакового размера (того же как исходное), где 𝐽 задано заранее. Каждое изображение с номером 𝑗 будет называться 𝑗 уровнем и обозначаться 𝑤𝑗 . Точку уровня 𝑗 с номером 𝑙 обозначим как 𝑤𝑗,𝑙 . Уровни упорядочены так, что первый отвечает за самую высокочастотную составляющую изображения, а последний за низкочастотную. Введем исходный сигнал {𝑐0,𝑙 } и масштабирующую функцию 23

𝜑(𝑥), которой соответствует НЧ фильтр. Вейвлет-функция удовлетворяет масштабирующему уравнению 2.14: 1 (︁ 𝑥 )︁ ∑︁ 𝑔(𝑘)𝜑(𝑥 − 𝑘) 𝜓 = 2 2 (3.1) 𝑘 Коэффиценты диискретного вейвлет преобразование вычисляются как: 𝑤𝑗+1,𝑙 = ∑︁ 𝑔(𝑘)𝑐𝑗,𝑙+2𝑗 𝑘 (3.2) 𝑘 В итоге каждый уровень вейвлет преобразования представляет собой разность сигналов: 𝑤𝑗+1 = 𝑐𝑗,𝑙 − 𝑐𝑗+1,𝑙 (3.3) Разность сигналов {𝑐0,𝑙 }−{𝑐1,𝑙 } содержит информацию о разности масштабов и этот ряд связан с вейвлет преобразованием, которое соответствует 𝜑(𝑥). Соответствующая вейвлет-функция определяется, как 1 (︁ 𝑥 )︁ 1 (︁ 𝑥 )︁ 𝜓 = 𝜑(𝑥) − 𝜑 2 2 2 2 (3.4) В итоге сигнал вычисляется, как 𝑐𝑗,𝑙 = ∑︁ ℎ(𝑘)𝑐𝑗−1,𝑙+2𝑗 𝑘 (3.5) 𝑘 Таким образом происходит свертка предыдущего уровня с фильтром, при этом расстояние между пикселями в свертке увеличивает как степень двойки. Коэффиценты {ℎ(𝑘)} определятся через масштабирующую функцию 𝜑(𝑥): 1 (︁ 𝑥 )︁ ∑︁ 𝜑 = ℎ(𝑙)𝜑(𝑥 − 𝑙) 2 2 𝑙 Полный алгоритм trous-вейвлет преобразования выглядит так: 24 (3.6)

1. j инициализируется нулем, а 𝑐0,𝑙 начальным сигналом. 2. Вычисляем свертку уровня 𝑗 с фильтром ℎ. Расстояние между пикселями в свертке равно 2𝑗 . 3. Вычисляем коэффиценты дискретного вейвлет-преобразоания, как разность между соседними уровнями: 𝑤𝑗+1,𝑙 = 𝑐𝑗,𝑙 − 𝑐𝑗+1,𝑙 (3.7) 4. Если 𝑗 не достигло заданного значения, увеличиваем 𝑗 на единицу и переходим к шагу 2.. 5. Множество 𝑊 = 𝑤1 , 𝑤2 , . . . , 𝑤𝐽 , 𝑐𝐽 является вейвлет преобразованием исходного сигнала. Можно легко осуществить обратное преобразование: 𝑐0,𝑙 = 𝑐𝐽,𝑙 + 𝐽 ∑︁ 𝑤𝑗,𝑙 (3.8) 𝑗 В работе в качестве масштабирующей функции использовался кубический 𝐵3 (𝑥)-сплайн (рис. 3.1): 𝜑(𝑥) = 𝐵3 (𝑥) = 1 (|𝑥 − 2|3 − 4|𝑥 − 1|3 + 6|𝑥|3 − 4|𝑥 + 1|3 + |𝑥 + 2|3 ) (3.9) 12 В этом случае коэффиценты матрицы свертки равны в одномерном слу- 25

б) а) Рис. 3.1: a) Масштабирующая функция 𝜑(𝑥) (кубический сплайн) б) Соотвествующая вейвлет функция 𝜓(𝑥) чае )︀ 1 3 1 1 , , , , 16 4 8 4 16 , а в двумерном: (︀ 1 ⎛ ⎞ 1 ⎜ 16 ⎟ ⎜1⎟ ⎜4⎟ ⎜ ⎟ ⎜ ⎟ (1/16, 1/4, 3/8, 1/4, 1/16) ⊗ ⎜ 38 ⎟ ⎜ ⎟ ⎜1⎟ ⎜4⎟ ⎝ ⎠ (3.10) 1 16 У данного алгоритма есть недостатки - отрицательные значения коэффицентов и отсутствие устойчивости. Среднее значение вейвлет коэффицентов - ноль. Все уровни засоряются шумами и при фильтрации можно отсечь значимые коэффиценты, что приводит при реконструкции к изменению яркости объектов. Поэтому для выполнения работы так же рассматривались многоуровневые преобразования, основанные не на вейвлет-анализе. Например, многоуровневое медианное преобразование, основанное на соответствующем нелинейном НЧ фильтре. Обозначим исходное изображение как 𝑓 , 𝑆 - заданное количество уровней разложения, тогда точный алгоритм выглядит так: 26

1. Обозначим 𝑐𝑗 = 𝑓 , где 𝑗 = 1. Пусть 𝑠 = 1. 2. Следующий уровень 𝑐𝑗+1 определяется как результат медианной фильтрации 𝑐𝑗 с окном 2𝑠 + 1. 3. Многоуровневые коэффиценты 𝑤𝑗+1 определяются разницей уровней: 𝑤𝑗+1 = 𝑐𝑗 − 𝑐𝑗+1 . 4. Перейдем на следующий уровень 𝑗 = 𝑗 + 1, увеличим окно в два раза 𝑠 = 2𝑠. Возвращаемся ко второму шагу, если 𝑗 < 𝑆. Исходный сигнал восстанавливается аналогично trous-вейвлет преобразованию: 𝑓 = 𝑐𝑝 + ∑︁ 𝑤𝑗 (3.11) 𝑗 Здесь коэффиценты 𝑤𝑗 не близки к нулю, поэтому яркость восстановленного изображения будет соотвествовать исходной. Представленные в данной главе методы декомпозиции являются многоуровневыми, то есть они переводят одно изображение в набор из 𝐽 +1 изображений, каждое из которых отвечает за свой масштаб исходного снимка. Рассмотрим наглядно примеры таких преобразований. 3.2 Одномерные многоуровневые преобразования Начнем с одномерного случая. Смоделируем нормальное распределение вида: −(𝑥−𝜇)2 1 𝑓 (𝑥, 𝜎, 𝜇) = √ 𝑒 2𝜎2 𝜎 2𝜋 (3.12) Возьмем 𝑁 = 10000 точек из промежутка [0, 10]. Пусть 𝜎 = 2, 𝜇 = 5. Добавим аддитивный белый гауссовский шум к получишемуся сигналу (Рис. 3.2). 27

1 1.2 0.9 1 0.8 0.8 0.7 0.6 0.6 0.5 0.4 0.4 0.3 0.2 0.2 0 0.1 0 -0.2 0 2 4 6 8 10 0 2 4 а) 6 8 10 б) Рис. 3.2: a) Исходное распределение б) Распределение с добавлением белого шума Результатом многоуровневого преобразования исходного распределения будут 6 графиков, каждый из которых отвечает за свой масштаб изображения (Рис. 3.3, 6 масштабов). 0.15 0.05 0.1 0.05 0 0 -0.05 -0.1 -0.15 -0.05 0 2 4 6 8 10 0 2 4 а) 6 8 10 6 8 10 6 8 10 б) 0.04 0.03 0.03 0.02 0.02 0.01 0.01 0 0 -0.01 -0.01 -0.02 -0.02 -0.03 -0.04 -0.03 0 2 4 6 8 10 0 2 4 в) г) 0.015 0.15 0.01 0.1 0.005 0.05 0 0 -0.005 -0.05 -0.01 -0.1 -0.015 -0.15 0 2 4 6 8 10 0 д) 2 4 е) Рис. 3.3: a) - д) - масштабы распределения, начиная от самых высоких частот, е) - сумма всех шумовых компонент На Рис. 3.4 избражен 6 уровень преобразования (самая НЧ компонента). 28

Сумма всех уровней в точности представляет собой исходное распределение. Интересно что, частоты уровней белого шума отличаются, это может быть связано с работой генератора случайных чисел при моделировании. 1 1.2 0.9 1 0.8 0.8 0.7 0.6 0.6 0.5 0.4 0.4 0.3 0.2 0.2 0 0.1 0 -0.2 0 2 4 6 8 10 0 2 4 а) 6 8 10 б) Рис. 3.4: a) 6 гармоника преобразования б) Сумма всех уровней Как видно из 3.3 и 3.4, используя trous-вейвлет алгоритм и подобные ему методы, можно рассматривать каждый масштаб отдельно в пространственных координатах. Это удобнее с практической точки зрения. Для того чтобы получить исходную гладкую кривую, можно сгладить НЧ масштаб, предварительно убрав все остальные гармоники. Рассмотрим пример экстракции линий в спектре ультраяркого ренгеновского источника (NGC 4395, ULX-1, Рис. 3.5) 3500 3000 2500 I(λ) 2000 1500 1000 500 0 -500 4000 4200 4400 4600 4800 5000 5200 λ Рис. 3.5: Спектр ULX-1, NGC 4395 29 5400 5600

Аналогично модельному решению, разделим спектр на 6 гармоник. Сразу отбросим самые высокочастотные и низкочастотные части (1 и 6). Пусть 𝑤𝑗,𝑙 - коэффиценты 𝐽 уровня (𝐽 = 2 . . . 5). Для каждого 𝐽 выберем порог,отвечающие пикам в спектрам и ниже которого обнулим все вейвлет-коэффиценты уровня. После этого сложим обратно все 𝐽 уровни преобразования (Рис. 3.6). Так как мы намеренно не суммировали самую ВЧ гармонику(Рис. 3.7) в итоге 𝐼(𝜆) изменилась примерно на 7 − 8%. Примеры, представленные в этом разделе были вычислены с помощью многоуровневого медианного преобразования. Существеный недостаток его с по сравнению с trous-вейвлет алгоритмом заключается во времени работы, как правило это 𝑂(𝑛2 log(𝑛)) в отличии от 𝑂(𝑛 log(𝑛)) у вейвлет преобразования. Если в одномерном случае для 𝑛 = 104 − 105 это не так критично, то для изображений, где число элементов больше 108 это существенная разница (об этом в следующем разделе). 3500 3000 2500 I(λ) 2000 1500 1000 500 0 4000 4200 4400 4600 4800 5000 5200 5400 λ Рис. 3.6: Отфильтрованный спектр ULX-1, NGC 4395 30 5600

300 200 I(λ) 100 0 -100 -200 -300 4000 4200 4400 4600 4800 5000 5200 5400 5600 λ Рис. 3.7: ВЧ компонента спектра 3.3 Двумерные многоуровневые преобразования Аналогично смоделируем гауссово распределение для двумерного случая и заполним поле белым шумом (Рис. 3.8, размер 500 × 500 пикселей). Для того чтобы выделить скопления можно воспользоваться бинаризацией по маске (выбрать на НЧ гармонике маску, далее бинаризовать по этому значению, все что ниже этого значения будет приравнено к 0). На рис. 3.8 представлено trous-вейвлет разложение (с использованием 𝐵3 -сплайна в качестве масштабирующей функции). Многоуровневое медианное преобразование показывает аналогичные результаты, средние значения коэффицентов в этих алгоритмах почти равны, это связано с отсутствием сильной ВЧ зашумленности при моделировании. 31

50 50 100 100 150 150 200 200 250 250 300 300 350 350 400 400 450 450 50 100 150 200 250 300 350 400 450 500 50 100 150 200 a) 250 300 350 400 450 500 300 350 400 450 500 б) 50 50 100 100 150 150 200 200 250 250 300 300 350 350 400 400 450 450 50 100 150 200 250 300 350 400 450 500 50 100 150 200 в) 250 г) 50 50 100 100 150 150 200 200 250 250 300 300 350 350 400 400 450 450 50 100 150 200 250 300 350 400 450 500 50 д) 100 150 200 250 300 350 400 450 500 е) Рис. 3.8: a) Исходная модель, б) - е) Гармоники, начиная с НЧ, масштабированные по яркости, в сумме дают первоначальное распределение Пусть 𝑊1 . . . 𝑊𝑘 - набор вейвлет коэффицентов каждого уровня, 𝑓 - исходное изображение. Сравним преобразования по нескольким критериям: 1. Средние значения вейвлет коэффицентов на каждом уровне 𝐽: < 𝑤𝑗,𝑘 >. Для trous-вейвлет приобразования среднее значение коэффицентов 0. Это может приводить к неустойчивости метода. Например, нарушению интенсивности объектов на итоговом изображении. 32

2. Разница между суммой вейвлет-шкал и интенсивностью исходного изоб∑︀ ражения: < 𝑊𝑘 − 𝑓 >. Это проверка первого критерия. 𝑘 3. Скорость работы в секундах. Важный параметр, особенно на больших 𝑁 (> 108 ). В качестве примеров снимков объектов будем использовать NGC 2997 (DFOSC, Danish 1.54m telescope, 27/02/98, размер 2036 × 2042 пикселей) и NGC 5194 (Kitt Peak National Observatory, 2.1m, 2001/03/28, 2013 × 2822). Исходные изображения на Рис. 3.9. а) б) Рис. 3.9: a) NGC 2997 б) NGC 5194 (M51, Галактика Водоворот) Выполним декомпозицию этих снимков на 𝐽 = 6 уровней. Все уровни разложения представлены в Приложении А для вейвлет-преобразования и в Приложении Б для медианного метода. Заметим, что в случае trous-вейвлет преобразования (Рис. А.2 и А.1) в местах ярких объектов на масштабах 𝐽 = 1 . . . 5 появляются "провалы"(нулевая интенсивность в окрестностях), которые заполняются при суммировании с НЧ гармоникой. Для изображения NGC 2997 метод является устойчивым: все < 𝑤𝑗,𝑘 > больше 0.2 для 𝑗 = 1 (чем больше уровень, тем больше среднее значение коэффицентов), суммирование гормоник да33

ет в точности исходное изображение: < ∑︀ 𝑊𝑘 − 𝑓 >= 0 для всех эле- 𝑘 ментов 𝑊𝑘 . В случае NGC 5194 среднее значение коэффицентов равно ∑︀ < 𝑤1,𝑘 >= 1.4 · 10−6 , однако < 𝑊𝑘 − 𝑓 >= 10−6 , то есть получаем прак𝑘 тически абсолютные точные интенсивности для всех элементов. В среде MATLAB вейвлет преобразование занимает 1.113 секунд для NGC 2997 и 1.141 секунд для NGC 5194 (core i3-2330M, 2.2GHz, 4 GB RAM). Скорость работы одно из главных преимуществ этого подхода. Многоуровневый медианный алгоритм (Рис. Б.2 и Б.1) устойчив, т.к. ко∑︀ эффиценты на всех уровнях больше нуля и для всех элементов: < 𝑊𝑘 − 𝑘 𝑓 >= 0. Преобразование выполняется за 163.294 секунды для NGC 2997 и 247.211 секунд для NGC 5194. Время выполнения можно сократить, используя быстрые реализации медианного фильтра (в работе использовался базовый алгоритм). Оба алгоритма показывают себя хорошо,однако у каждого есть свои недостатки: возможная неустойчивость и провалы на ярких объектах у вейвлет-преобразования и медленная работа медианного метода. Несмотря на это, они вполне применимы на практике. 34

Глава 4 Многоуровневая обработка астрономических снимков В этой главе рассматривается подход к фильтрации изображений, используя многоуровневые разложения, построенные с помощью trous-вейвлет преобразования или медианного алгоритма. 4.1 Модель сигнал-шум Пусть 𝐼(𝑥, 𝑦) — исходное изображение. Тогда его можно представить в виде структуры: 𝐼(𝑥, 𝑦) = 𝑛 ∑︁ 𝑂𝑖 (𝑥, 𝑦) + 𝐹 (𝑥, 𝑦) + 𝐵(𝑥, 𝑦) (4.1) 𝑖=1 где 𝑂1 . . . 𝑂𝑛 — набор из 𝑛 объектов снимка, 𝐹 (𝑥, 𝑦) — фон, 𝐵(𝑥, 𝑦) — шумовая компонента. С другой стороны изображение можно представить как сумму вейвлет шкал: 𝐼(𝑥, 𝑦) = 𝐽 ∑︁ 𝑗=1 35 𝑊𝑗 (4.2)

Здесь 𝑊𝑗 — коэффиценты 𝑗 масштаба, 𝐽 - количество уровней разложения. Все наборы коэффиценетов {𝑊1 . . . 𝑊𝐽 } определены в пространственных координатах. Как уже упоминалось ранее, большинство снимков сильно зашумленны, это приводит к переполнению всех уровней коэффицентами, которые не являются значимыми (т.е. не несут информацию об объектах изображения). Основная идея подхода к фильтрации снимка заключается в последовательной очистке каждого уровня от элементов, которые относятся к шумовой компоненте. Затем изображение восстанавливается при сложении всех масштабов. Многоуровневые разложения являются мощным инструментом в обрабтке астрономических изображений, они позволяют легко анализировать любые масштабы снимка. Ниже будут предложены несколько вариантов применения такого подхода. 4.1.1 Soft and hard thresholding 𝐻𝑎𝑟𝑑 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑𝑖𝑛𝑔 — алгоритм фильтрации изображения по пороговому значению. Все вейвлет коэффиценты ниже определенного значения 𝑇𝑗 приравниваются к 0 (считаются незначимыми): 𝑤 ̃︂ 𝑗,𝑘 = ⎧ ⎪ ⎨𝑤𝑗,𝑘 если |𝑤𝑗,𝑘 | ≥ 𝑇𝑗 ⎪ ⎩0 если 𝑤𝑗,𝑘 < 𝑇𝑗 (4.3) где 𝑤𝑗,𝑘 - вейвлет коэффицент масштаба 𝑗 в позиции 𝑘. 𝑆𝑜𝑓 𝑡 𝑡ℎ𝑟𝑒𝑠ℎ𝑜𝑙𝑑𝑖𝑛𝑔 преполагает замену каждого коэффицента на значение 𝑤 ̂︂ 𝑗,𝑘 : 𝑤 ̂︂ 𝑗,𝑘 = ⎧ ⎪ ⎨sign(𝑤𝑗,𝑘 )(𝑤𝑗,𝑘 − 𝑇𝑗 ) если |𝑤𝑗,𝑘 | ≥ 𝑇𝑗 ⎪ ⎩0 если 𝑤𝑗,𝑘 < 𝑇𝑗 36 (4.4)

В случае Гауссовского шума 𝑇𝑗 = 𝐾𝜎𝑗 , 𝑗 - уровень масштаба, 𝜎𝑗 — отклонение от шума на шкале 𝑗. Значение 𝐾 чаще всего принимается за 3. Рассмотрим снимок Юпитера (HST, 2007-04-27, 6730Å). Заполним изображение равномерным Гауссовским шумом. Разложим снимок на 6 гармоник при помощи медианного преобразования и применим фильтрацию по порогу с 4𝜎 к каждому масштабу. На Рис. 4.1 видно, что алгоритм показывает хорошие результаты: изображение немного потеряло в детализации, но границы сохранены (большое преимущество в сравнении с сглаживающими фильтрами типа Гаусса). а) б) б) Рис. 4.1: a) Оригинальный снимок б) Зашумленное изображение в) Изображение после применения пороговой обработки 4.1.2 3𝜎-clipping и винеровская фильтрация Одним из самых эффективных методов по экстракции фона изображения является 3𝜎-clipping, этот подход используется в программном обеспечении DAOPHOT (Stetson’s DAOPHOT program, [6]). Cнимок итеративно делится на части, пока элементы находящиеся в ±3𝜎 от среднего значения (Median) не меняются больше, чем на 20% за проход. Тогда можно сказать, что поле не является переполненным звездами ("crowded fields") и вычислить среднее значение фона (Mean). Отфильтрованное изображение представляется в виде: 37

Image = 2.5 × Median − 1.5 × Mean (4.5) На практике используют формулу вида: Image = 3 × Median − 2 × Mean (4.6) Аналогичный подход использован для выделения фона в SExtractor (Bertin E., S. Arnouts, [3]). Так же хорошие результаты показывает алгоритм адаптивной винеровской фильтрации для подавления аддитивного гауссова белого шума. Данный метод основан на статистических оценках фрагментов изображения в пределах скользящего окна размера 𝑚 × 𝑛 пикселей (подробнее в [31], [32]). Для всех положений скользящего окна с центральным пикселем и координатами (x,y) среднее значение яркости равно: 𝜇= 1 ∑︁ 𝐼(𝑥, 𝑦) 𝑚𝑛 (4.7) (𝑥,𝑦)∈𝜂 где 𝜂 — окрестность 𝑚 × 𝑛 вокруг каждого пикселя изображения 𝐼. Дисперсия вычисляется как: 1 ∑︁ 2 𝜎 = 𝐼 (𝑥, 𝑦) − 𝜇2 𝑚𝑛 2 (4.8) (𝑥,𝑦)∈𝜂 Данная формула применяется для всех положений окна. Отфильтрованное изображение записывается в виде: 𝜎2 − 𝜈 2 𝐹 (𝑥, 𝑦) = 𝜇 + (𝐼(𝑥, 𝑦) − 𝜇) 𝜎2 (4.9) Здесь 𝜈 2 — мощность белого шума. Можно оценить 𝜈 как среднее по всем 𝜎 на снимке. 38

Приведем пример оценки самых мелкомасштабных колебаний интенсивности фона снимка. Рассмотрим изображение M33 (Рис. 4.2, KPNO 4.0m telescope, 2000-10-04, Massey, фильтр V 5386.76Å, размер снимка 8278 × 8384) Рис. 4.2: M33, 5386.76Å, Massey, ds9 scale 92.5% В Приложении В приведены разложения снимка по 6 гармоникам с помощью изложенных выше алгоритмов. Рассмотрим самый высокочастотный масштаб: обнулим очевидные шумы (< 2 × Mean) и применим винеровскую фильтрацию с окном 5 × 5, оценивая мощность шума, как среднее отколение по участку. В итоге на Рис. 4.3 видны ВЧ гармоники и фон изображения. На изображении стало заметно больше деталей - наблюдается спиральная структура вблизи ядра, ярко выраженные рукава галактики. Интересная особенность - после фильтрации видно, что снимок разделен на 8 частей (в шапке FITS-файла указано: Number of CCDs = 8). 39

Рис. 4.3: M33, ВЧ гармоника, ds9 scale 92.5% 40

Глава 5 Многоуровневое представление данных Одно из основных применений математических алгоритмов попиксельной обработки в астрономии - это поиск и отождествление объектов на снимках с целью составления каталогов с основными параметрами (видимая звездная величина, размеры, координаты) и их дальнейшей классификации. Ниже будет представлен подход к идентификации объектов на изображении при помощи многоуровневого представления данных. 5.1 Описание многоуровневой модели Пусть 𝑓 — исходное изображение. Построим его бинарную модель из многоуровневого представления: если на уровне 𝑗 содержится информация в пикселе с координатами (𝑖, 𝑘), то его значение 𝑀𝑖,𝑘,𝑗 = 1, остальные элементы приравниваются к 0. Таким образом модель строится следующим образом: 1. Преобразование исходного изображения 𝑓 с помощью многоуровневого алгоритма в структуру вида {𝑤1 , . . . , 𝑤𝑗 }, где 𝑤𝑘 — набор коэффицентов для масштаба 𝑘. 41

2. Определение значимых элементов на каждом уровне. 3. Бинаризация всех уровней по определенному порогу значимости. В итоге многоуровневая модель данных (Muliscale Vision Model — MVM) представляет собой следующий набор структур: ∙ Наборы значимых вейвлет коэффицентов {𝑊1 , . . . , 𝑊𝑗 }, определяемых при помощью алгоритмов, приведенных в Главе 4. ∙ Структуры 𝑆𝑗,𝑘 , которые представляют собой наборы значимых вейвлет коэффицентов на одной шкале 𝑗: 𝑆𝑗,𝑘 = {𝑤𝑗 [𝑘1 , 𝑙1 ], 𝑤𝑗 [𝑘2 , 𝑙2 ], . . . , 𝑤𝑗 [𝑘𝑝 , 𝑙𝑝 ]} (5.1) где p — число значимых коэффицентов, входящих в структуру 𝑆𝑗,𝑘 , 𝑤𝑗,𝑥𝑖 ,𝑦𝑖 — вейвлет коэффицент шкалы 𝑗 в позиции (𝑥𝑖 , 𝑦𝑖 ). ∙ Объекты — наборы структур 𝑆𝑗,𝑘 : 𝑂𝑙 = {𝑆𝑗1 ,𝑘1 , . . . , 𝑆𝑗𝑛 ,𝑘𝑛 } (5.2) Так же определим оператор ℒ, отображающий принадлежность структуры данному объекту: ℒ(𝑆𝑖,𝑗 ) = 𝑙, если 𝑆𝑖,𝑘 ∈ 𝑂𝑙 и ℒ(𝑆𝑖,𝑗 ) = 0 в остальных случаях. ∙ Уровень масштаба объекта, который определяется уровнем вейвлет коэффицента с максимальным значением. ∙ Набор связей между структурами: набор правил по которым можно отнести несколько структур к одному объекту. ∙ Подобъекты, которые имеют свои и максимумы в вейвлет пространстве, но сами являются частью более крупного объекта. 42

5.2 Поиск объектов и структур После применения многоуровневого представления снимка мы получим несколько бинарных изображений. Как уже было оговорено ранее, если вейвлет коэффицент в позиции (𝑥, 𝑦) является значимым, то соотвествующее ему значение пикселя равно 1, в остальных случаях 0. Последовательности ненулевых пикселей 𝑖 образуют структуры 𝑆𝑖,𝑗 на масштабе 𝑗. Любой объект представляет собой набор иерархических структур. Правила, которые определяют какие структуры входят в объект определяются набором связей между структурами (см. Раздел 5.1). Если пиксель с максимальным значением структуры 𝑆𝑖,𝑗 входит в структуру 𝑆𝑖,𝑗+1 , тогда эти структуры связаны. Полные наборы связанных структур являются объектами. Данная система наглядно представлена на Рис. 5.1. Много объектов и связных структур появляется из-за зашумленности данных, поэтому важно применять эффективную фильтрацию к каждому уровню. На Рис. 5.2 изображена структура суммы с 1 по 5 уровня разложения снимка M33 из Приложения В и отдельно пятая гармоника. Можно видеть как объекты по отдельности входят в структуры. В вейвлет пространстве решения о разделении и слиянии структур, определяется по локальным максимум на различных масштабах. Cтрук𝑚 тура определяет объект на уровне 𝑗, если выполняется условие: 𝑤𝑗𝑚 > 𝑤𝑗−1 𝑚 и 𝑤𝑗𝑚 > 𝑤𝑗+1 , где: ∙ 𝑤𝑗𝑚 — масимальный вейвлет коэффицент структуры 𝑗 𝑚 ∙ 𝑤𝑗−1 = 0, если коэффицент на связан ни с какой структурой масштаба 𝑗 − 1 или максимальное значение 𝑆𝑗−1,𝑙 , где ℒ(𝑆𝑗−1,𝑘′ ) = ℒ(𝑆𝑗,𝑘 ). 𝑚 ∙ 𝑤𝑗+1 = 𝑚𝑎𝑥 {𝑤𝑗+1,𝑥1 ,𝑦1 , . . . , 𝑤𝑗+1,𝑥𝑛 ,𝑦𝑛 }, где 𝑤𝑗,𝑥,𝑦 — коэффиценты 𝑆𝑗,𝑘 (все 𝑤𝑗,𝑥,𝑦 ∈ 𝑆𝑗,𝑘 ). 43

Рис. 5.1: Схематическое представление связности структур. Пример взят из работы [14] б) а) Рис. 5.2: a) сумма пяти гармоник разложения М33 б) Пятый уровень преобразования 5.3 Алгоритм поиска объектов Метод поиска объектов на снимках можно разделить на несколько этапов: 44

1. Выполнить декомпозицию изображения, получив набор коэффицентов {𝑤1 , . . . , 𝑤𝑗 } для каждого уровня. 2. Определить отклонение шума на каждом масштабе 𝑗. 3. Определить пороговые значения для шума на каждом уровне. 4. Отдельно выполнить фильтрацию для каждого масштаба. 5. Определить межуровневые отношения по бинаризованным уровням. 6. Найти все максимальные значения коэффицентов для каждой структуры уровня. 7. Вычислить все связанные структуры и разделить их. 5.4 SExtractor SExtractor (Source-Extractor, Bertin E., S. Arnouts, [3]) — программное обеспечение для идентификации объектов на астрономических снимках и составления каталогов. Алгоритм экстракции объектов, использованный авторами, представляет собой некоторое подобие много уровневого подхода (схематично метод изображен на Рис. 5.3): ∙ Large Image — исходный снимок в FITS-формате. ∙ Background Algorithm — использован 3𝜎-clipping алгоритм для экстракции фона изображения (описан в разделе 4.1.2, использована формула 4.6) ∙ Thresholding включается в себя свертку с фильтром (SExtractor хорошо оптимизирован для свертки матриц с маленьким ядром) и опредление порогового значения для объектов (подробнее об этом ниже) 45

∙ Deblending — разделение структур на источники, используя граф связности и априорные представления об искомых объектах. Размеры объектов вычисляются по бинарной карте. ∙ Photometry — фотометрия выделенных объектов. Рис. 5.3: Алгоритм работы Source-Extractor, Bertin E., S. Arnouts, [3] Перед определением пороговых значений выполняется свертка изображения с одним из следующих фильтров: ∙ Набор сглаживающих фильтров Гаусса (3 × 3, 5 × 5 и 9 × 9). Применяются для идентификации объектов с 1.5 < FWHM < 5. ∙ Top-hat — морфологическая фильтрация изображения, вычитает разомкнутое изображение (морфологическое замыкание) из исходного. Используется для поиска слабых протяженных объектов (FWHM < 5). 46

∙ Mex-hat — преобразования при помощи набора вейвлет функций типа "мексиканская шляпа". Применяется в переполненых звездами полях. ∙ Small pyramidal function - быстрая свертка (по умолчанию). В зависимости от типа искомых объектов, может быть выбран узкий фильтр (например, для звезд это PSF), либо широкий (для небольших скоплений матрица Гаусса 9 × 9, для туманностей хорошие результаты показывает Top-hat). Таким образом, комбинацией разных фильтров можно выделить все структуры снимка, в этом заключается некое подобие многоуровневого подхода, примененного в SExtractor. Настройка экстракции производится через конфигурационный файл, который содержит в себе такие основные параметры, как: ∙ BACK_SIZE — размер элемента сетки при выделении фона. ∙ DETECT_THRESH — пороговое значение для определения значимости пикселей. ∙ DETECT_MINAREA и DETECT_MAXAREA — границы размеров определяемых объектов. ∙ FILTER_NAME и FILTER — выбор типа фильтра для свертки. ∙ DEBLEND_MINCONT и DEBLEND_NTHRESH — минимальная разница в яркости для раздеделения объекта на струкуры и максимальное количество разделений. Весь список параметров представлен в [7] и [3]. Особое внимание стоит уделить выбору порога и фильтра, так как эти характеристики очень сильно влияют на качество выборки. 47

5.4.1 Поиск рассеянных структур и скоплений В Приложении Г приведен пример работы SExtractor в поиске галактик в Hubble Ultra Deep Field 2009 (F105W). Изначально это программное обеспечение создавалось для выделения галактик и подобных структур на снимках, поэтому на данном поле не возникает проблем с поиском и используются стандартные параметры. Пример выделения рассеяных структур и небольших скоплений (𝐻𝛼 туманностей М33) приведен на Рис. 5.4. Часто SExtractor используют для поиска небольших скоплений на больших картах. Например, такая работа представлена в [8]. Выделение скоплений производилось по нескольким признакам: ∙ Параметр DETECT_THRESH для определение порога индентификациии. ∙ Размеры выделяемого объекта: FWHM, Area. ∙ Вытянутость — отношение малой полуоси к большой. ∙ Видимая звездная величина объекта. Рис. 5.4: SExtractor, пример выделения объектов и туманностей M33, 𝐻𝛼 , Massey, KPNO 48

5.4.2 Проблемы выделения звезд и скоплений Существует проблема выделения объектов в переполненых звездами полях, особенно если стоит задача выделения отдельных звезд или небольших групп. Причем понижание порога для алгоритма 4.6 или комбинирование фильтров (уменьшение размера сглаживающей матрицы Гаусса) не очень помогают в этой ситуации. Пример работы программы с параметром -OBJECTS для разных фильтров представлен на Рис. 5.5 (в качестве примера использована центральная область галактики M33, Massey, KPNO). Первые три изображения на Рис. 5.5 — результат поиска с применением малых сглаживающих фильтров. Как можно видеть, они показывают хорошие идентичные результаты, но отождествить отдельные звезды в таких заполенных полях проблематично. Авторы программы Bertin E. и Arnouts S рекомендуют в [3] использовать Mex-hat преобразование для скученных полей, если стоит цель отождествления малых объектов и структур (изображение г), Рис. 5.5). Метод позволяет идентифицировать некоторые звезды по отдельности, но большая часть информации потеряна. Уменьшение порога < 1𝜎 не улучшает результат. Попробуем использовать для этого многоуровневое вейвлет преобразование, представленное в разделе 3.3, и виенеровскую фильтрацию 4.6 и пороговую обработку 4.1.2 на разных масштабаъ исходного изображения (Рис. 5.7). Комбинацией параметров фильтров можно выделять или скрывать источники. Чтобы выделить группы звезд, можно применить Гауссовское сглаживание, описанное в 1.1.1. 49

a) б) г) в) Рис. 5.5: Вывод промемужточных изображений с помощью SExtractor (OBJECTS) с использованием разных фильтров: a) Матрица Гаусса 3x3 б) Top-hat в) Без использования фильтра г) Mex-hat преобразование Как видно на Рис. 5.6 trous-вейвлет преобразование с применением виенеровской фильтрации показывает хорошие результаты в выделении точечных источников, лучше чем предложенное в SExtractor Mex-hat преобразование с экстракцией фона при помощи 3𝜎-clipping алгоритма. Преимущество второго подхода заключается в высокой скорости работы, большей автоматизированности и универсальности. 50

Рис. 5.6: Центр снимка M33, trous-вейвлет разложение (Рис. 5.7), виенеровская фильтрация, пороговая обработка а) б) Рис. 5.7: a) M33, Massey, оригинал б) Отфильтрованный снимок 51

Глава 6 Поиск звездных скоплений в М33 В данной главе будет рассмотрен пример поиска скоплений в галактике М33. Основной задачей стояло выделение молодых групп звезд, которые могут быть связаны с яркими голубыми переменными (LBV — luminous blue variable). 6.1 Luminous blue variable Яркие голубые переменные (LBVs) — одна из наименее понятных стадий в эволюции массивных звезд (до 150 𝑀⊙ ). Трудности в понимании этих объектов определяются как малочисленностью известных LBV-звезд, их разнообразием, сильной спектральной и фотометрической переменностью [19, 20], так и неточным знанием их фундаментальных параметров. Очевидно, что массивная звезда может проявить себя как LBV после стадии О-звезды Главной Последовательности и перед стадией поздней WRзвезды азотной последовательности (WNL). Однако, более точное положение LBV-звезд в переходе OV → WNL пока неясно [15]. Весьма вероятно, что стадия LBV и соответствующие этой стадии неустойчивости вет52

ров определяется содержанием водорода в атмосфере звезды. Однако, как определение содержания водорода, так и точное определение светимости и температуры звезды задача нелегкая. Определение фундаментальных параметров массивных звезд в близких галактиках, например в богатой массивными звездами галактике М 33, может быть более уверенным чем в нашей Галактике, так как расстояния до галактик измерены с достаточной точностью. У LBV-звезд наблюдаются сильные изменения блеска и спектра при примерно постоянной болометрической светимости. Характерные времена такой глобальной переменности составляют месяцы–годы. При увеличении визуального блеска LBV-звезды температура фотосферы понижается до 9000-10000 К, размер фотосферы звезды увеличивается, это высокое/холодное состояние. Здесь и далее под фотосферой мы понимаем ”псевдофотосферу”, то есть то место в ветре звезды, где формируется наблюдаемое континуальное излучение. При понижении визуального блеска температура заметно возрастает до 35000 K или даже выше [16, 17, 19], размер звезды уменьшается, это низкое/горячее состояние. В горячем состоянии спектры LBV весьма похожи на спектры WNL-звезд, однако, истинные звезды типа WN7/8/9 отличаются от LBV в низком/горячем состоянии рядом спектральных особенностей [15, 21] и низким содержанием водорода. Так называемые переходные (in transition или "slash") звезды типа Ofpe/WN9, выделенные в отдельный класс звезд Walborn [22], похожие на звезды LBV в низком/горячем состоянии. Ofpe/WN9-звезды характеризуются смешанным спектром с эмиссионными линиями Of и WNL звезд. На примере двух галактических звезд типа LBV AG Car и He 3–519, которые в низком/горячем состоянии соответствовали типу Ofpe/WN9, Smith, Crowther & Prinja [21] предложили продление WNL-классификации до WN10 и WN11. Тем не менее, до сих пор остается открытым вопрос, отличаются ли ис53

тинные звезды LBV в низком/горячем состоянии от самых поздних звезд WN9/10/11. Другими словами, возможно ли, что все звезды самых поздних типов WNL являются звездами LBV, которые за время их наблюдений (всего несколько десятков лет) не проявили сильной переменности блеска. Такие звезды, которые можно подозревать в принадлежности к типу LBV, принято называть спящими (dormant) LBV. Заметим, что спящие LBVзвезды могут быть не только среди звезд типа WNL. Две классическме звезды LBV, 𝜂 Car и 𝑃 Cyg, сейчас имеют спектр OB- гипергигантов, переменность их блеска мала. Однако, несколько сот лет назад у них наблюдались гигантские вспышки (giant eruptions) [19]. Если бы не эти гигантские поярчания 𝜂 Car и 𝑃 Cyg и не знаменитая туманность вокруг звезды 𝜂 Car, сейчас мы не смогли бы заподозрить их в принадлежности к классу LBV. N-body simulation предлагают следующую модель: на стадии образования скопления, плотный центр выбрасывает массивную звезду(LBV), и она начинает удалятся от скопления. Для того, чтобы проверить эту теорию, предлагается идентифицировать все скопления, которые могут связаны с LBV и статистически оценить вероятность такого распределения. LBV — массивные и молодые звезды с возрастом < 107 лет. Поэтому в работе велся поиск групп звезд, рядом с которыми находятся области HII с высоким излучением в 𝐻𝛼 . Если звезды формировались в этих областях, то их возраст не больше 107 − 108 лет. 6.2 Выбор исходных данных В качестве исходных были использованы снимки M33, выполненные группой Massey (KPNO, 4.0m) в B, V, 𝐻𝛼 фильтрах. Изображение в диапазоне V уже было использовано в предыдущих разделах (например, Рис.4.2). Снимки в фильтрах B (4381Å) и 𝐻𝛼 (6563Å),рассмотренные в работе, представлены на Рис. 6.1. 54

а) б) Рис. 6.1: a) B (4381Å) б) 𝐻𝛼 (6563Å) M33, Massey, ds9 scale 92.5% 6.3 Определение параметров групп звезд Для работы были выбраны два списка звезда: подтвержденные LBV (8 звезд) и 29 кандидатов в LBV-звезды в М33(данные были взяты из [24, 25] В начале было необходимо вручную (визуально оценить размеры по снимкам 𝐵, 𝐵, 𝐻𝛼 в RGB цветах) замерить характерные параметры скоплений, которые могут связаны с LBV-звездами. Все характеристики представлены в Приложении Д(Табл.Д.1,Д.2 для кандидатов, Д.3,Д.4 для подтвержденных LBV). Гистограммы распределение параметров представлены на Рис.6.2 для кандидатов и Рис.6.3 для подтвержденных LBV. Эти диаграммы далее использовались для подбора критериев поиска объектов на снимках. 55

14 20 12 10 Nebula Clusters 15 10 8 6 4 5 2 0 0 0 5 10 15 20 Cluster size, arcseconds 25 30 0 10 15 20 25 Nebula size, arcseconds 30 35 б) 20 20 15 15 Clusters Clusters а) 5 10 5 10 5 0 0 0 0.2 0.4 2Deltax/(a+b) 0.6 0.8 0 в) 15 30 45 60 75 90 Distance to nearest LBV, arc seconds г) Рис. 6.2: a), б) — распределение размеров ((𝑎 + 𝑏)/2 из Табл. Д.3,Д.4), в) — гистограммма для 2Δ𝑥/(𝑎 + 𝑏), г) — распределение расстояний для ближайшей LBV 56

5 4 4 3 3 Nebula Clusters 5 2 2 1 1 0 0 0 5 10 15 20 Cluster size, arcseconds 25 30 0 5 30 35 б) 4 4 3 3 Clusters Clusters а) 10 15 20 25 Nebula size, arcseconds 2 1 2 1 0 0 0 0.2 0.4 2Delta x/(a+b) 0.6 0.8 0 15 30 45 60 75 90 Distance to nearest LBV, arc seconds в) г) Рис. 6.3: a), б) — распределение размеров ((𝑎 + 𝑏)/2 из Табл. Д.1,Д.2), в) — гистограммма для 2Δ𝑥/(𝑎 + 𝑏), г) — распределение расстояний для ближайшей LBV 6.4 Идентификация объектов Эта часть диплома является продолжением работы по выделению скоплений М33, начатой в прошлом году. В начале предполагалось из каталога звезд M33 (около 150000 звезд, Massey) выделить скопления,которые могут быть связаны с LBV-звездами. Модель заключалась в том, что каждая звезда представляла собой точку с заданными координатами (𝛼, 𝛿). Первым был реализован самый очевидный способ — метод звездных подсчетов. Он заключается в том, чтобы поделить карту на участки и по57

считать количество звезд в каждом их них, а затем вычислить среднеквадратическое отклонение этой величины от значений окружающих ее участков (например, в данном случае брали участок 3×3 с центром в определяющей клетке). Таким образом, скопления идентифицировались как участки с локальным повышением плотности звезд. Схематично это представлено на Рис.6.4. Выбирались секторы с отклонением > 5𝜎, центры которых наносились на карту вместе с LBV-кандидатами. В этой работе вся область была поделена по участкам размером 10′′ , это оптимально, чтобы в каждой клетке было не более одного скопления. Рис. 6.4: Метод звездных подсчетов, схематичное представление Однако этот алгоритм плохо работал в областях с высокой плотностью звезд, метод удалось улучшить при помощи отсечения гармоник по секторам в переполненных звездами зонах, однако одним из его главных недостатков являлось то, что он не давал никакой информации о скоплении кроме координат сектора, в котором оно находится. После этого был использован альтернативный подход. Применялась аналогичная прдыдущему методу модель: на карту наносились звезды, как точки с координатами (𝛼, 𝛿), были выбраны голубые звезды c 𝐵 −𝑉 < 0.5 и 𝑚 < 22𝑚 (Рис. 6.5). К полученной карте применялось вейвлет-преобразование 58

и подбирался оптимальный масштаб для выделения скоплений. Далее применялся hard-thresholding, маска выбиралась экспериментально по границе предполагаемых скоплений (Например, Рис.6.6). Для каждой обособленной группы пикселей можно определить центр масс (это будет центра скопления) и площадь, которую она занимает (некий критерий отбора). Подход показывает хорошие результаты в областях с высокой плотностью (Рис. 6.6), однако большинство предполагаемых скоплений выделить не удалось — если взять недостаточно крупный масштаб, то после пороговой обработки останутся все точечные источники. При увеличении масштаба сглаживания сильно теряется качество выборки, в таком случае невозможно выбрать скопления по формальным критериям. Рис. 6.5: Карта звезд по каталогу Massey 59

Рис. 6.6: Центральная часть M33, пороговая обработка Можно идентифицировать изображения по снимкам М33. Как видно из Рис.6.2,6.3, искомые объекты имеют сравнительно крупные размеры. Используем SExtractor для экстракции объектов по сглаженным и отфильтрованным изображениям (надо подобрать максимально широкий фильтр, не теряя информации о размерах объектов в центральных областях). Основные критерии поиска скоплений — размеры и близость к 𝐻𝛼 туманности. В начале были выделены все возможные объекты на снимке в B и 𝐻𝛼 фильтрах, а затем выбраны туманности по критериям из Рис.6.2,6.3(размеры от 8′′ до 23′′ ). После этого составлено распределение (Рис.6.7) 2Δ𝑥/(𝑎 + 𝑏) по каждому скоплению, где Δ𝑥 — расстояние от скопления до ближайшей туманности, 𝑎, 𝑏 — размеры туманности. Будем считать, что если этот параметр меньше 0.7, то скопление находится рядом с туманностью. 60

Рис. 6.7: Распределение 2Δ𝑥/(𝑎 + 𝑏) по всем выделенным объектам Далее были выделены предполагаемые скопления по размеру (𝑎, 𝑏 взяты из Рис.6.2,6.3, 9′′ < (𝑎 + 𝑏)/2 < 14′′ ). После получено распределение по расстояниям до ближашей LBV для каждого скопления (Рис.6.8). Всего было выделено 109 объектов. Рис. 6.8: Распределение расстояний до ближайших LBV по всем выделенным объектам Для того, чтобы расчитать вероятность такого распределения, используем метод моделирования Монте-Карло. На Рис.6.9 представлен снимок в фильтре 𝐻𝛼 c нанесенными скоплениями, LBV-звездами и кандидатами. 61

Рис. 6.9: M33, 𝐻𝛼 , Massey. × — смоделированные скопления, ∘ — LBV и кандидаты. Нанесем на карту 𝑁 = 1000 случайных точек и составим гистограмму относительно расстояний до ближайшей LBV для каждой точки. Попробуем аппроксимировать получившиеся графики разными распределениями (Рис.6.10). Как видно, распределение растет быстро до 100′′ , а затем медленно убывает, поэтому точно аппроксимировать его нормальным распределением не получается (Рис.6.10). На Рис.6.11 представлены диаграммы из Рис.6.2 с соответствующими им гамма и бета распределения. Из графиков видно, что они существенно отличаются от смоделированных методом Монте-Карло. Это можно строго доказать. Классическим примером критерия принадлежности выбранных данных к какому-либо закону распредления является критерий 𝜒2 (критерий согласия Пирсона). Для сравнения данных полученных практическим путем и 62

теоретического закона распредления используется формула 2 𝜒 =𝑁 𝑁 ∑︁ ((𝑝𝑡𝑖 − 𝑝𝑒𝑖 )2 /𝑝𝑡𝑖 ) (6.1) 𝑖=1 где 𝑝𝑡𝑖 и 𝑝𝑒𝑖 — вероятности попадания величин в выбранный интервал 𝑖. Главным недостатком критерия 𝜒2 является необходимость строгой привязки начальных данных к интервалам и объединении интервалов при малом числе данных. Поэтому в работе для доказательства неслучайного распределения скоплений использовался критерий Колмогорова-Смирнова (KS test). а) б) в) г) Рис. 6.10: a) — распределение для моделирования Монте-Карло 𝑁 = 1000 б) — аппроксимация нормальным распределением в) — гаммараспределение г) — бета-распределение 63

а) б) Рис. 6.11: a),б) — распределение для LBV и кандидатов из Рис.6.2, аппроксимированные гамма и бета распределениями Формула для вычисления критерия Колмогорова-Смирнова представляет собой: 𝐷𝑚,𝑛 = 𝑠𝑢𝑝|𝐹𝑚 (𝑥) − 𝐺𝑚 (𝑥)| (6.2) где 𝐷𝑚,𝑛 — расхождение между распределениями, 𝐹𝑚 (𝑥),𝐺𝑚 (𝑥) — распределение частоты встречаемости значений. Для KS-теста использовался язык статистической обработки данных R(см. [36]). Статистический анализ двух распределений проводился при помощи p-value (уровень значимости). Проведем численный эксперимент: 100 раз смоделируем Монте-Карло точки на карте (𝑁 = 1000) и составим распределение p-value. Для каждой итерации значение p-value не превышает 10−5 . Это значит, что с вероятностью > 99.99% эти распределения не соответствуют друг другу. На этом основании можно сделать вывод, что распределение выделенных объектов не случайно относительно LBV-звезд. Это c высокой долей вероятности означает, что они вылетели из молодых скоплений на стадии их формирования. 64

Заключение В работе были изложены основные методы обработки астрономических снимков при помощи пространственных и адаптивных фильтров,пороговой обработки, Фурье и вейвлет-анализа. В первой части диплома на примерах показано преимущество многоуровневых преобразований над стандартными методами фильтрации изображений (например, Фурье). Далее представлен алгоритм trous-вейвлет преобразования, алгоритм и его отладка на модельных решениях — одномерных, двумерных(гауссов шум) и реальных данных (спектр NGC 4395, снимки NGC 2997, NGC 5194). Во второй части работы изложены алгоритмы, используемые в современном астрономическом программном обеспечении, например 3𝜎-clipping или пороговая обработка(hard-thresholding), а так же их адаптация в рамках многоуровневой модели. Такой подход к фильтрации показал очень хорошие результаты, что показано в Гл.4. Далее в рамках многоуровневой модели был представлен алгоритм поиска объектов на снимках, а так же примеры использования ПО(SExtractor) с аналогичным подходом, который был улучшен в некоторых случаях с помощью представленных преобразований (например, при поиске объектов в переполненных звездами областях). В последней главе представленные инструменты были использованы для поиска молодых звездных образований в галактике М33. Статистически доказана связь LBV-звезд с найденными скоплениями. 65

Литература [1] Albert Bijaoui, A multiscale vision model adapted to the astronomical images, Signal Processing 46 (1995) 345-362 [2] Albert Bijaoui, Christophe Benoistand, Eric Slezak, Analysis of multiband astronomical images using multiscale tools, Casiopee laboratory, B.P. 4229, 06304 Nice Cedex 4, France [3] Bertin E., S. Arnouts SExtractor: Software for source extraction — Astronomy & Astrophysics Supplemen, 1996 [4] Bertin, E. and Arnouts, S.: 1996, Astronomy and Astrophysics, Supplement Series 117, 393 [5] Eric H. Neilsen, Jr., Notes on the Essentials of Astronomy Data, 2014-0616 [6] DAOPHOT: A COMPUTER PROGRAM FOR CROWDED-FIELD STELLAR PHOTOMETRY, PETER B. STETSON, Publications of the astronomical society of the pacific, March 1987 [7] Dr. Benne, W. Holwerda, Source extractor for dummies, Baltimore MD 21218 [8] Izaskun San Roman and Ata Sarajedini, Photometric Properties of the M33 Star Cluster System, Department of Astronomy, University of Florida, 211 Bryant Space Science Center, Gainesville, FL 32611-2055 66

[9] Mallat S. A Theory for Multiresolution Signal Decomposition: The Wavelet Representation // IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989. — Vol. 11. — P. 674–693 [10] Nathan Smith, and Ryan Tombleson, Luminous blue variables are antisocial: their isolation implies that they are kicked mass gainers in binary evolution, MNRAS 447, 598–617 (2015) [11] JL. Starck, F. Murtagh Handbook of Astronomical Data Analysis - 2002 [12] JL. Starck, F. Murtagh, Albert Bijaoui, Image processing and data analysis: The multiscale approach [13] JL. Starck, Astronomical image representation by the curvelet transform, A&A 398, 785–800 (2003) [14] Jean-Luc Starck, Fionn Murtagh, Mario Bertero, Starlet Transform in Astronomical Data Processing, Springer Science+Business Media LLC 2011 [15] N. Smith and P. S. Conti, Astrophys. J. 679, 1467 (2008). [16] T. Szeifert, LIACo 33, 459 (1996). [17] J. S. Clark, V. M. Larionov, and A. Arkharov, Astronom. and Astrophys. 435, 239 (2005). [18] Roberta M. Humphreys, Kerstin Weis, Kris Davidson and Michael S. Gordon, On the Social Traits of LBVs, Minnesota Institute for Astrophysics, University of Minnesota, Astronomical Institute, RuhrUniversitaet Bochum, Germany, Minnesota Institute for Astrophysics, University of Minnesota, 2016 [19] R. M. Humphreys and K. Davidson, Publ. Astronom. Soc. Pacific 106, 1025 (1994). 67

[20] A. M. van Genderen, Astronom. and Astrophys. 366, 508 (2001). [21] L. J. Smith, P. A. Crowther, and R. K. Prinja, Astronom. and Astrophys. 281, 833 (1994). [22] N. R. Walborn, Astrophys. J. 215, 53 (1977). [23] Shashikiran Ganesh, Physical Research Laboratory, Wavelets in the Astronomical Context, 2003 [24] Humphreys, Roberta M.; Weis, Kerstin; Davidson, Kris; Bomans, D. J.; Burggraf, Birgitta, Luminous and Variable Stars in M31 and M33. II. Luminous Blue Variables, Candidate LBVs, Fe II Emission Line Stars, and Other Supergiants, The Astrophysical Journal, 790: 48 (21pp), 2014 July [25] Humphreys, Roberta M.; Martin, John C.; Gordon, Michael S., A New Luminous Blue Variable in M31, PUBLICATIONS OF THE ASTRONOMICAL SOCIETY OF THE PACIFIC, 127: 347–350, 2015 April [26] А.Н. Яковлев, Введение в вейвлет-преобразования, Новосибирск: Изд-во НГТУ, 2003 [27] Добеши И., Десять лекций по вейвлетам, 2001, 464 стр. [28] И.М. Дремин, О.В. Иванов, В.А. Нечитайло, Вейвлеты и их использование, Успехи физических наук, май 2001, том 171, № 5 [29] М. А. Иванов, Применение вейвлет-преобразований в кодировании изображений, работа по гранту № 01-01-794 [30] Н.Ф. Астафьева, Вейвлет-анализ: основы теории и примеры применения, Успехи физических наук, ноябрь 1996, том 166, № 11 68

[31] Прэтт У Цифровая обработка изображений В 2 т М Мир, 198 [32] Писаревский А. Н. и др. Системы технического зрения (Принципиальные основы, аппаратное и математическое обеспечение) Л: Машиностроение, 1988. [33] Р. Гонсалес, Р. Вудс, Цифровая обработка изображений, ТЕХНОСФЕРА, Москва, 2005 [34] О.Н. Шолухова, С.Н. Фабрика, А.В. Жарова, А.Ф. Валеев, В.П. Горанский, Спетральная переменность LBV-звезды V 532 (ROMANO’S STAR), Аастрофизический Бюллетень, 2011, том 66, № 2, с. 135–157 [35] https://www.gnu.org/software/octave/ [36] https://www.r-project.org/ 69

Приложение А а) б) в) г) д) е) Рис. А.1: a) - е) - разложение снимка NGC 5194 (trous-вейвлет преобразование) 70

а) б) в) г) д) е) Рис. А.2: a) - е) - разложение снимка NGC 2997 (trous-вейвлет преобразование) 71

Приложение Б а) б) в) г) д) е) Рис. Б.1: a) - е) - разложение снимка NGC 5194 (медианный алгоритм) 72

а) б) в) г) д) е) Рис. Б.2: a) - е) - разложение снимка NGC 2997 (медианный алгоритм) 73

Приложение В а) б) в) г) д) е) Рис. В.1: a) - е) - разложение снимка M33 (trous-вейвлет преобразование) 74

а) б) в) г) д) е) Рис. В.2: a) - е) - разложение снимка M33 (медианный алгоритм) 75

Приложение Г а) б) в) г) д) е) Рис. Г.1: Hubble Ultra Deep Field 2009 (F105W): a) Original б) Apertures (выделение границ) в) OBJECTS г) -OBJECTS (снимок без объектов) д) SEGMENTATION (бинаризованные уровни) e) BACKGROUND (фон) 76

Приложение Д Таблица Д.1: Измеренные вручную размеры скоплений для подтвержденных LBV и кандидатов, 𝛼,𝛿 — координаты в формате h:m:s, 𝑎,𝑏 — малая и большая полуось в угловых секундах 𝛼 скопления 𝛿 скопления 𝑎 скопления 𝑏 скопления 1:32:35 30:30:26 8.26 11.6 1:32:35 30:30:26 8.26 11.6 1:33:26 30:24:07 11.86 15.0 1:33:41 30:22:37 12.0 15.8 1:33:59 30:23:20 12.0 19.0 1:34:16 30:28:26 12.0 14.0 1:34:18 30:31:37 13.0 17.5 1:34:42 30:31:40 5.5 13.0 1:35:07 30:41:46 12.5 25.0 1:35:00 30:37:18 5.0 8.5 1:32:37 30:40:07 15.0 16.0 1:32:45 30:38:58 24.0 36.0 1:32:45 30:38:58 24.0 36.0 1:33:11 30:48:57 12.8 30.0 77

1:33:28 30:31:31 12.0 15.0 1:32:59 30:33:33 10.5 17.5 1:34:22 30:33:13 11.5 15.5 1:34:20 30:33:57 8.0 14.0 1:34:15 30:33:41 12.0 15.0 1:34:18 30:34:10 9.5 10.0 1:34:10 30:34:38 9.5 15.0 1:34:16 30:36:42 14.0 16.0 1:34:18 30:38:36 7.0 12.5 1:34:32 30:47:10 25.0 29.0 1:33:40 30:45:55 12.0 16.0 1:33:26 30:39:50 12.0 17.0 1:33:33 30:41:34 24.0 25.0 1:33:40 30:41:38 10.0 11.0 1:33:35 30:36:28 9.0 12.0 1:33:55 30:45:26 8.0 14.0 1:34:06 30:41:44 12.0 14.0 1:33:51 30:40:56 7.0 11.0 1:33:33 30:33:44 7.5 14.0 1:33:50 30:41:30 9.0 14.0 1:33:47 30:38:39 8.0 10.0 1:33:53 30:38:53 4.8 7.5 1:33:52 30:39:15 5.0 6.0 78

Таблица Д.2: Измеренные вручную размеры туманностей для подтвержденных LBV и кандидатов, 𝛼,𝛿 — координаты в формате h:m:s, 𝑎,𝑏 — малая и большая полуось в угловых секундах, Δ𝑥 — расстояние между центрами скопления и туманности в угловых секундах 𝛼 скопления 𝛿 скопления 𝑎 скопления 𝑏 скопления Δ𝑥 1:32:35 30:30:28 12.9 17.3 30.0 1:32:35 30:30:28 12.9 17.3 10.0 1:33:26 30:24:08 7.7 9.51 46.0 1:33:41 30:22:42 9.5 15.6 4.0 1:33:59 30:23:20 15.3 21.4 80.0 1:34:16 30:28:22 17.0 22.0 15.0 1:34:18 30:31:36 14.0 15.0 34.0 1:34:42 30:31:40 9.0 13.0 36.0 1:35:07 30:41:51 15.0 28.8 35.0 1:35:00 30:37:20 6.0 10.8 21.0 1:32:38 30:39:58 49.0 63.0 1.0 1:32:45 30:38:59 38.0 45.0 5.0 1:32:45 30:38:59 38.0 45.0 64.0 1:33:11 30:48:54 40.0 45.0 64.0 1:33:28 30:31:32 14.0 23.0 64.0 1:33:00 30:33:38 21.0 24.0 2.0 1:34:22 30:33:17 24.0 24.5 34.0 1:34:19 30:33:58 7.0 14.0 82.0 1:34:15 30:33:45 16.0 21.0 7.0 1:34:18 30:34:11 15.0 18.0 1.0 1:34:10 30:34:34 30.0 37.0 1.0 79

1:34:16 30:36:41 14.0 22.0 4.5 1:34:18 30:38:36 9.0 14.0 4.0 1:34:32 30:47:10 65.0 95.0 30.0 1:33:40 30:45:58 11.5 16.5 19.0 1:33:26 30:39:06 14.5 22.5 9.5 1:33:33 30:41:34 28.0 32.0 16.0 1:33:40 30:41:37 14.0 15.0 3.0 1:33:35 30:36:28 9.0 12.0 30.0 1:33:55 30:45:26 16.0 27.0 8.0 1:34:06 30:41:45 15.0 20.0 4.0 1:33:51 30:40:56 11.0 16.0 7.0 1:33:33 30:33:46 16.0 8.0 1.0 1:33:50 30:41:30 14.0 20.0 3.0 1:33:47 30:38:41 14.0 21.0 35.0 1:33:53 30:38:52 16.0 20.0 4.0 1:33:52 30:39:15 7.0 7.5 7.0 Таблица Д.3: Измеренные вручную размеры скоплений для подтвержденных LBV, 𝛼,𝛿 — координаты в формате h:m:s, 𝑎,𝑏 — малая и большая полуось в угловых секундах 𝛼 скопления 𝛿 скопления 𝑎 скопления 𝑏 скопления 1:32:35 30:30:26 8.26 11.6 1:34:16 30:28:26 12.0 14.0 1:35:07 30:41:46 12.5 25.0 1:33:28 30:31:31 12.0 15.0 1:34:15 30:33:41 12.0 15.0 1:34:10 30:34:38 9.5 15.0 80

1:33:35 30:36:28 9.0 12.0 1:33:47 30:38:39 8.0 10.0 1:33:52 30:39:15 5.0 6.0 Таблица Д.4: Измеренные вручную размеры туманностей для подтвержденных LBV, 𝛼,𝛿 — координаты в формате h:m:s, 𝑎,𝑏 — малая и большая полуось в угловых секундах, Δ𝑥 — расстояние между центрами скопления и туманности в угловых секундах 𝛼 скопления 𝛿 скопления 𝑎 скопления 𝑏 скопления Δ𝑥 1:32:35 30:30:28 12.9 17.3 30.0 1:34:16 30:28:22 17.0 22.0 15.0 1:35:07 30:41:51 15.0 28.8 35.0 1:33:28 30:31:32 14.0 23.0 64.0 1:34:15 30:33:45 16.0 21.0 7.0 1:34:10 30:34:34 30.0 37.0 1.0 1:33:35 30:36:28 9.0 12.0 30.0 1:33:47 30:38:41 14.0 21.0 35.0 1:33:52 30:39:15 7.0 7.5 7.0 81

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв