Двумерная вейвлет томография галактических полей Текст научной статьи по специальности «Физика»

УДК 621.391: 519.2

ДВУМЕРНАЯ ВЕЙВЛЕТ-ТОМОГРАФИЯ ГАЛАКТИЧЕСКИХ ПОЛЕЙ

Р.А. Степанов (Пермь)

Abstract

A special two-dimensional algorithm of wavelet-analysis is suggested, which reconstruct the large-scale structure of a field from measured integrals between the arbitrary point and the pole. The method is useful for the analysis of integrated quantities of the magnetic field of Galaxy (Faraday rotation, dispersion measures). The algorithm is tested by artificial examples. It is shown that the large-scale filtration is possible even under very inhomogeneous distribution of the data. The method does not use any information a priori and can be considered as an algorithm of wavelet-interpolation.

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

Примером задач такого типа является проблема реконструкции крупномасштабной составляющей магнитного поля нашей Галактики, где Земля представляет выделенную точку наблюдения. В первом приближении галактический диск можно представить как плоскость, поэтому сразу ограничимся рассмотрением только двумерных полей в галактической плоскости. Причем наличие центральной точки делает удобным применение полярных координат (полюс - Земля), которые и будут использоваться в дальнейшем. Измерение магнитного поля в межзвездной среде возможно только косвенным путем. Наиболее приемлемыми для анализа являются независимые измерения фарадеевских мер вращения RM и дисперсионных мер DM, определяемых выражениями [1]

L L

RM = Jn'tB • dl, DM = Jn'tdl, (1)

0 0

где B - магнитное поле, ne - плотность свободных электронов, L - расстояние до

источника поляризованного излучения, в направлении которого и ведется интегрирование. Существуют два типа источников: внегалактические (излучение проходит через все магнитное поле нашей Галактики и расстояние до источника не играет роли) и внутригалактические (как правило, это пульсары, при рассмотрении которых учет расстояния просто необходим). В данной работе рассматривается второй тип источников. Определить радиальную компоненту магнитного поля простым дифференцированием по верхнему пределу нельзя, поскольку измерения возможны далеко не в любой точке Галактики. Для простоты далее будет рассматриваться только DM. На данный момент преобладает мнение, что свободные электроны в межзвездном пространстве концентрируются в торе и примыкающим к ним спиральным ветвям (рис.1,а), а пульсары, доступные для измерений, концентрируются

10

+ + +

-10

* % < ІіУ»*. + + ++ +

5* ”Ш$Л * Ї +* + * * .

* -10

а б

Рис.1 Экспериментальные данные: (а) - предполагаемое поле пе (здесь и далее интенсивность черного цвета пропорциональна величине), (б) - распределение данных измерений.

в окрестности Солнца (рис.1,б) (данные на рис.1. представлены в галактической плоскости с Солнцем в начале координат и галактическим центром в точке (8,5;0)). В такой ситуации точного определения поля быть не может. Цель анализа заключается в выделении крупномасштабных структур поля из общей картины данных.

Вейвлеты на данный момент являются одним из самых удобных инструментов анализа многомерных полей путем разделения пространственных структур с различными размерами. Непрерывное двумерное вейвлет-преобразование в полярных координатах определяется выражением [2]

П то

(а г, ф) = Су а_1/2 11уо г,ф(г', ф,)/(г', ф?) г' ёг' ф ', (2)

уw(a, г, ф) = Су а а,г,ф

-п0

где /(г', ф') - анализируемая функция, сф - константа, уогф (г', ф') - семейство вейвлетов, которое получается масштабированием и сдвигом исходного вейвлета ¥(г) .

2 Г2

В дальнейшем будет использоваться вейвлет ¥(г) = (2 - г )ехр(- —) , называемый

«мексиканской шляпой».

Построение алгоритма. Предположим, что известны интегралы от / (г, ф)

Г

Я(г, Ф) = | /(г' , ф¥г'. (3)

о

тэ г, ч дя (г, ф)

В аналитической постановке задача восстановления решается легко: / (г, ф) =-----------.

дг

Однако проблема заключается в том, что реально измерения я (г, ф) заданы только в конечном наборе точек. В результате, необходимо интерполировать не только саму функцию, но и ее производную. Если теперь принять во внимание наличие ошибки измерения (а измерительных данных может оказаться недостаточно для статистики), то сразу сталкиваемся с проблемой устойчивости. Аналогичные трудности возникают в задачах восстановления двумерной функции по зашумленным проекционным данным [3]. Структура вейвлет-преобразования позволяет обойти эту проблему. Запишем выражение (2) с учетом (3),

w(a,r,Ф) = Су a"1'2 j/у(т\ф')r' dr’<kf . (4)

-п0

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

ТО

w(a,r,ф) = Су a“1/2(r ' a’r’^ 'ф ) + у(r' , Ф'))g(r' , Ф')

a,r,ф (r\ф') У«,г,ф (r ФО

(5)

ПТО

-1/2 f Г/дУ a,r,Ф (' , ф ) , У a,r,(f (' , ф Л , , , 7 , 7 ,

- СУ a JJ(--------------------------дТ.--------+-----------Г-----------------------------)g(r , ф ) r dr ф .

-п0

Л / f f\ ’ ’т ’ ’Т /^7\

У a,r,<p (Г , Ф ) = -----------------------------^----------------------------------3---------------. (7)

Первое слагаемое (5) обращается в ноль вследствие того, что, с одной стороны, вейвлет локализован в физическом пространстве, т.е. уa,r (r' , ф') ^ 0 при r' ^ то , с другой

стороны, по определению функции g(r, ф) выполняется условие g(0, ф) = 0. В конечном счете, имеем

п то

w(a, r, ф) = Су a ~m j J у a,r,ф (r ', ф ')g(r ', ф ') r 'dr Уф ' , (6)

-п 0

где у a,r,ф (r', ф' ) - новое семейство вейвлетов, которое связано с начальным

семейством соотношением

ду a,r,ф (r ' , ф ' ) у a,r,ф (r ' , ф 1 dr' r'

Таким образом, преобразуя семейство вейвлетов, удается избежать дифференцирования по функции g(r, ф). Необходимо отметить, что предлагаемое семейство вейвлетов (7) должно удовлетворять общим требованиям для вейвлет-функции [4]. Важными для дальнейшего рассмотрения являются спектральные свойства вейвлетов, а именно: нечувствительность вейвлет-разложения к среднему значению функции g(r, ф), а также менее востребованное аналогичное свойство по отношению к функции f (r, ф) . В математической форме они имеют следующий вид:

П ТО П ТО

j j у a,r,ф (r', ф')g(r', ф') r 'dr Уф' = 0, j j у a,r,ф (r', ф')g(r', ф') r'2dr Уф ' = 0. (8)

-п0 -п0

На практике значительные трудности возникают при вычислении интеграла (6). Данные измерений, которые требуется проанализировать, представляются в виде

дискретного набора { = g(rt, фг-)} . Более того, для астрофизических наблюдений они

имеют нерегулярное и неравномерное распределение (см. рис.1,б). Для численного нахождения вейвлет-коэффициентов интеграл (6) заменяется на сумму:

N

w(a,r, ф) = Суа_1/2 ^ у a,r,ф (r', ф' )g(r', ф' )g,dsi, (9)

i=1

где dsi - площадь, соответствующая i -й точке. Если данные распределены равномерно, то в первом приближении можно считать, что dsi = ds = const. Чем более неоднородно

распределение, тем ошибочнее результаты интегрирования. При этом важным остается сохранение на такой сетке спектральных свойств вейвлета (8). Для решения этой проблемы удобно воспользоваться идеей адаптивных вейвлетов [5,6]. Для реализации этой идеи представим исходный вейвлет Y(r) в виде Y(r) = h(r)Ф(г), где Ф(г) есть положительно определенная масштабная функция («оболочка»), а h(r) -

«заполняющая» функция. Адаптивный вейвлет Y(r) будем искать в виде

¥ (г) = (й(г) + С 0 + С1 г )Ф(г).

(10)

Константы С0, С! находятся для каждого масштаба а и положения вейвлета (т,ф), и

условия (8). Отметим еще один технический прием, позволяющий уменьшить ошибки, связанные с неоднородностью распределения данных. В сумму (9) вносят вклад только те точки, которые находятся в некоторой окрестности центра вейвлета. Можно оценить ds для каждого масштаба а и центра вейвлета (г, <р), используя масштабную функцию

ds =

/2 п а

N

ЕФ а,г,Ф (ГФ і )

і=1

(11)

Величина ds также может служить количественным критерием достаточности данных.

Тестирование алгоритма. Для иллюстрации алгоритма рассмотрим тестовые примеры, максимально приближенные к реальным данным. Мы будем использовать модельный сигнал (рис.2), воспроизводящий основные особенности предполагаемой структуры пе. Важным вопросом является выбор сетки, на которой задается сигнал. Из

общих соображений положим, что для задания некой двумерной структуры необходимо, чтобы расстояние между точками не превосходило половины её масштаба. Минимальный масштаб, которым мы интересуемся, ~1 килопарсека (кпс), следовательно, равномерная сетка с шагом 0,5 кпс вполне достаточна, тогда в области 20 х 20 кпс должно быть порядка 1600 точек. Реальность такова, что на данный момент число измерений в 4-5 раз меньше и их распределение далеко не равномерно. Поэтому рассмотрим три тестовых примера, переходя от «идеальной» сетки к реальной.

10

-10

р 10 Г^1И

5

0 \Мш

мн -5

л! 11Я -10

-10 -5

0

10

-10 -5

0

10

а б

Рис.2 Тестовый сигнал: (а) - поле пе , (б) - полученная из пе дисперсионная мера БЫ

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

а б

Рис.3 Вейвлет-коэффициенты для фиксированного масштаба: (а) -а=3 (выделение кольцевой структуры), (б) - а=1 (спиральная структура)

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

а б

Рис.4 Вейвлет-коэффициенты для фиксированного масштаба: (а) -а=3 (выделение кольцевой структуры), (б) - а=1 (спиральная структура)

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

С общих позиций проблема, исследуемая в статье, может рассматриваться как задача интерполяции. В этом плане предлагаемый алгоритм позволяет избежать процедуры прямого интерполирования как производной, так и самой исследуемой функции, а также использования какой-либо априорной информации. Невозможно что-либо сказать о структурах, размер которых сравним с размерами пробелов в данных.

а б

Рис.5 Вейвлет-коэффициенты для фиксированного масштаба: (а) -а=3 (выделение кольцевой структуры), (б) - а=1 (спиральная структура)

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

Библиографический список

1.Рузмайкин А.А., Соколов Д.Д., Шукуров А.М. Магнитные поля галактик. - М.: Наука, 1988. - 280 с.

2.Holschneider M. Wavelets. An analysis Tool. Oxford: Oxford University Press, 1995.

3.Патрикеев И. А. Фрик П.Г. Вейвлет-томография в условиях шума // Мат.моделир. сист. и процессов: Сб. научн. тр./ Перм. гос. техн. ун-т. Пермь, - 1997. - №5. - С. 8692

4.Фрик П.Г. Вейвлет-анализ и иерархические модели турбулентности: Препринт / ИМСС УоР РАН. Пермь, 1992

5.Галягин Д.К. Фрик П.Г. Адаптивные вейвлеты // Мат.моделир. сист. и процессов: Сб. научн. тр./ Перм. гос. техн. ун-т. Пермь, - 1997. - №2. - С. 20-28

6. Frick P., Grossmann A., Tchamichian Ph. Wavelet analysis of signals with gaps //Journal of Mathematical Physics, 1998. - Vol.39. - №.8. - P.4091-4107.

Получено 30.05.99.