Нелинейная регуляризация обращения экспоненциального преобразования Радона Текст научной статьи по специальности «Математика»

7. Корнилов В. И. Пространственные пристенные турбулентные течения в угловых конфигурациях. Новосибирск: Наука, 2001.

8. Самарский А. А., Н и к о л ае в Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.

9. Воеводин В. В., Воеводин Вл. В. Параллельные вычисления. СПб.: БХВ-Петербург, 2002.

10. Березин С. Б., Пасконов В.М. Компонентная система визуализации результатов расчетов на многопроцессорных вычислительных системах // Тр. конф. "Высокопроизводительные вычисления и их приложения". Черноголовка, 2000. С. 202-203.

11. Johnston J. On the three-dimensional boundary layer generated by secondary flow //J. Basic Engin. Trans. ASME. 1960. Ser. D. 1.

12. Пасконов В.М. Модификация уравнений Навье-Стокса для расчета вязких течений разностными методами // Обратные задачи естествознания. М.: Изд-во МГУ, 1997. С. 189-197.

13. Р a s к о п о v V. М. Modified Navier-Stockes equation for finite-difference computation of viscous flow // Сотр. Math, and Modeling. 1997. 8. N 4. P. 400-407.

Поступила в редакцию 10.10.05

УДК 519.2

0. В. Шестаков

НЕЛИНЕЙНАЯ РЕГУЛЯРИЗАЦИЯ ОБРАЩЕНИЯ ЭКСПОНЕНЦИАЛЬНОГО ПРЕОБРАЗОВАНИЯ РАДОНА

(кафедра математической статистики факультета ВМиК)

1. Введение. Задачи однофотонной эмиссионной томографии приводят к проблеме обращения так называемого экспоненциального преобразования Радона в Л2 [1], которое отличается от классического преобразования Радона наличием экспоненциального множителя у интегрируемой функции. В настоящее время существует множество методов восстановления функции по ее экспоненциальному преобразованию Радона, среди которых методы фильтрации обратных проекций и итерационные методы.

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

Предлагаемый в данной работе метод был разработан и применен для классического преобразования Радона в работе [2] и исследован численно в работе [3].

2. Экспоненциальное преобразование Радона. Пусть V = 51 X Л, где 51 — единичный круг в Л2. Экспоненциальным преобразованием Радона для ¡1 6 И. и функции /(ж, у), имеющей компактный носитель (без ограничения общности будем считать этим носителем круг II с единичным радиусом и центром в начале координат), называется преобразование вида

оо

J е»{х>е±)${х)<1х= J е^ЦзО + М^М, (0, в) 6 V.

{х,в) = 8 -ОО

Заметим, что если ¡1 = 0, то мы имеем классическое преобразование Радона. Сопряженным оператором к оператору является оператор обратного проецирования, задаваемый преобразованием вида

Е*д(х) = J е"{х'в±)д(0, {х,в))йв, х Е В2.

Б1

Многие результаты теории классического преобразования Радона переносятся на экспоненциальное преобразование Радона. Пусть / — преобразование Фурье функции /. Проекционная теорема обобщается на экспоненциальное преобразование Радона следующим образом:

^/(0,и) = (27г)1^/(и0 + ^0±), (1)

т.е. одномерное преобразование Фурье функции в) по второй переменной задает двумерное

преобразование Фурье функции /(ж) на некоторой поверхности в С2. Формула обращения экспоненциального преобразования Радона выглядит следующим образом:

(2)

где 1~1 — обобщенный потенциал Рисса:

1 0, иначе.

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

Начнем с описания базиса в £2(К). Фиксируем I 6 (2) и образуем все двоичные сжатия и сдвиги функций (р и ф:

<р11к(х) = 21/2ср(21х - к), ф^к(х) = 2>!2ф{2>х -к), ], к е (г), ] > I. Если определить вейвлет-коэффициенты

(Ч,к = (<Р1,к, /) = J <Р1,к(х)/(х)с1х, = (ф^к, /) = J ф^к{х)${х)<1х,

то функцию / 6 Ь2 (И) можно представить в виде

к 3^1 к

Первая сумма в (3) представляет собой "грубую аппроксимацию" функции /, а вторая сумма представляет собой добавляемые к этой аппроксимации детали.

В данной работе мы будем иметь дело с вейвлет-базисом в £2(К2). Ортонормальный вейвлет-базис в Ь2 (И.2) может быть составлен из всех двоичных сжатий и сдвигов функций Ф(ж, у) = (р(х)(р[у), (ж, у) = <р(х)ф(у), Ф№(х,у) = ф(х)<р(у) и фИ(ж,у) = ф(х)ф(у), где (р и ф соответственно "отцовская" и "материнская" вейвлет-функции. В £2(К2) справедливо разложение, аналогичное (3):

к т к

где

Щ,к = (&1,к,1) = ! J ^1,к(х,у)1(х,у)йхйу, = = ! J ^\тк(х,у)1(х,у) ¿хйу.

В силу компактности носителя функции / в сумме (4) для каждого I и ] сумма по к содержит лишь конечное число членов.

4. Метод реконструкции. Следуя [2], мы предлагаем метод реконструкции, основанный на выражении коэффициентов разложения в сумме (4) через Для этого необходимо найти такую последовательность что

<ф/,*,/> = д„/>

и

<ФЙ]./> = <сй],^/>.

(5)

(6)

Из формулы (2) следует, что Далее

Таким образом Аналогично

СЙ] = (4 тгГ^Д^/ФЙ. (7)

Если функции <*р т ф дважды непрерывно дифференцируемы и имеют компактный носитель, то формулы (6) и (7) определяют корректные £2-объекты. Кроме того, последовательность {ид} =

= {б,к/ 116,кII , С^к/ & } образует "почти" ортонормальный базис в где £ Ь2(В?) —

множество всех функций из £2(К2), имеющих носителем круг С/, т.е. существуют константы с > О, С > 0 такие, что для любой последовательности а\ Е I2:

2

а

а-

Если функция / выражается формулой (4), содержащей любое конечное число членов, то ее можно восстановить по Я^/ (слабо обратить оператор Я^) по формуле

/ = + ЕЕ Е<СЙ] > • (8)

к т к

4. Некорректность задачи реконструкции. Исследуем нормы элементов последовательности {&,к, С^ }• По теореме о среднем сечении имеем

(2тг)1/2

—»>0 ,к)

2'(2тг)1/2

(2тт)112 -Цшв-^е1- ,к) „ / и!0 — ¡пО^Л 1

— ^-£-Р 9.' ф I --- I — -Р ■>'

2'

2'

Ф(г) йг = /

= 2« (Д_м/2гФ)(^/2

Таким образом,

где I — функция индикатора. Аналогично

1

4тг2-?

Далее

(9) (10)

1М|2= (4^у Н2КН(Н) (Д_,/2!ф)(с/2',0)

з1 о

йи =

(4тг)

ОС

I I Ы21(|д|(|2Ч|)|(Д^ Ф)^!,^)2^! (11)

С1'

2-?

(4тг)

е — ^ / 1(1^1(12^1) (Д_^ФМ)(Ш1,0)

¿и л

(12)

откуда в силу компактности носителя и гладкости функций Ф и ф[т1 и компактности носителя функции / следует, что норма элементов присутствующих в разложении (8), ограничена при фиксированном /, а норма , присутствующих в разложении (8), возрастает при росте ] со скоростью ~ 2•?/2, что и характеризует некорректность задачи обращения экспоненциального преобразования Радона, заключающуюся в том, что даже малые изменения в наблюдаемых данных приводят к значительным отклонениям получаемого результата от истинного решения.

5. Пороговая обработка вейвлет-коэффициентов. Метод пороговой обработки вейвлет-коэффициентов был предложен для регуляризации задачи обращения классического преобразования Радона в работе [2] и численно исследован в работе [3]. Он имеет множество преимуществ по сравнению с традиционными линейными методами регуляризации как с теоретической, так и с практической точек зрения. Метод заключается в применении к вейвлет-коэффициентам в разложении (4) оператора мягкой пороговой обработки к (ж) = sgn(ж)(|ж| — + для коэффициентов и

ы. к (ж) = (ж)(|ж| - для коэффициентов ,

/ = Е Ки + Е Е Е ч. «<ё]> > (13)

к т к

т.е. выбирается некоторый порог — такой, что если абсолютная величина коэффициента меньше этого порога, то он обнуляется, а иначе его абсолютная величина уменьшается на величину порога (мягкая пороговая обработка). Можно также оставлять коэффициенты, превосходящие порог по абсолютной величине, неизменными (жесткая пороговая обработка), однако в этом случае функции пороговой обработки Ь,^ к{х) и получаются разрывными, что приводит к дополнительным трудностям [4].

Чтобы пояснить идею, лежащую в основе этого метода, отвлечемся от задачи обращения экспоненциального преобразования Радона и рассмотрим только задачу подавления шума в наблюдениях, исследованную в работе [2]:

Уг = 1{г/п) + г = 0,... ,га - 1,

где — независимые случайные величины, распределенные по закону А^О,*?2). Вычислим вейвлет-преобразование данных, применим к коэффициентам оператор мягкой пороговой обработки /1(п (ж) и вычислим обратное вейвлет-преобразование. Порог £п в этом случае выбирается равным (га)]1'2*?

[3], потому что для независимых случайных величин го,.. имеющих стандартное нормальное распределение,

Р(тах|,гг|)[2к^(га)]1/2) О

при га —т- оо. Поскольку вейвлет-преобразование является ортогональным, аддитивные шумовые добавки в вейвлет-коэффициентах также независимы и имеют нормальное распределение А^О,*?2). Следовательно, [2к^(га)]1/2*? является верхней вероятностной границей уровня шума в пространстве эмпирических вейвлет-коэффициентов [3]. Для неоднородных сигналов и изображений вейвлет-разложение

(4) содержит относительно небольшое число больших по абсолютной величине информативных коэффициентов, а шум равномерно распределяется по всем коэффициентам [3]. Следовательно, пороговая обработка удаляет большинство коэффициентов, соответствующих шуму, оставляя большинство информативных коэффициентов.

Перейдем теперь к определению порогов ¿/^ и tjtk для разложения (8). На практике мы имеем конечный набор данных вида

Уг,ь = + г = 0,..., Мапд - 1, £ = 0, ...,ЛГро8 - 1, (14)

где — независимые случайные величины, распределенные по закону М(0,а2). Эмпирические коэффициенты в разложении (8) вычисляются по формулам

/V —1 /V —1

(&,к,У) = —--- ^ Е

ро81У* апд г=0 г=0

и

/V —1 /V —1

, 1 * апд 1 1 * ро э

ров1у апд г=0 г=0

Из этого определения следует, что коэффициенты распределены нормально с математическими ожиданиями

и дисперсиями

3,к

С

N N 1

14 апд ров <

гДе смапд,мроа = 4к/(МапдМров).

Поскольку носителем функции / является единичный круг, в разложении (4) (и (8)) присутствует только 22'+1 коэффициентов а/^ и 22-7'"1"1 коэффициентов с?^ для каждого уровня ] и каждого т.

"Почти" ортогональность {ид} обеспечивает "почти" независимость эмпирических коэффициентов, и, как показано в работе [2], байесовский минимаксный риск рассматриваемой модели асимптотически эквивалентен байесовскому минимаксному риску модели с независимыми эмпирическими коэффициентами.

Наша модель в терминах и выглядит следующим образом:

(Ь,к, Г) « Д„/> + г™ = (Ф,,к, /) + г\1 <С™, У) « <С$, Д„/> + = <фЦ, /) +

где

3,к

С

N N

1 у апд

Таким образом, уровень шума не является постоянным. В частности, дисперсия растет со скоростью ~ 23, что является прямым следствием некорректности задачи.

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

для 22'+1 коэффициентов Д,*/) и

С1' 1

3,к

С

1/2

/V /V

1 * апд р

для 22э+1 коэффициентов на каждом уровне ] для каждого т. Вид этих порогов анало-

гичен используемым в работе [3], но зависит не только от уровня но и от положения к. Чтобы сократить вычисления, для больших ] можно точный вид tjtk заменить для к = (кх,ку), близких к (0,0), выражением

Ц,к « [21оё(22^1)]1/22^2(7Ао,оС^1ЛГро8,

где

ДпФ[г

где

Л,о =

с1в

(До ФН)(ш,9)

(4тг)2

з1 о

классическое преобразование Радона от фМ

, и

с!и)

1/2

/V '

(4тг)

_ ,к) е 23 ¿в

(До фН)(ш,9)

с!и)

1/2

для к = (кх,ку), далеких от (0,0), поскольку при выполнении условий, накладываемых на выбор функций Ф[т],

оо

I м21(Н(|2Ч)

о о

при j —7- оо.

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

1. Федоров Г.А. Вычислительная эмиссионная томография. М., 1990.

2. Don oh о D.L. Nonlinear solution of linear inverse problems by wavelet-vaguelette decomposition // Appl. Comput. Harm. Anal. 1995. 2. P. 101-126.

3. Kolaczyk E.D. A wavelet shrinkage approach to tomographic image reconstruction //J. Am. Statist. Assoc. 1996. 91. P. 1079-1090.

4. Johnstone I.M., Silverman B. W. Wavelet threshold estimators for data with correlated noise // J. Roy. Statist. Soc. Ser. 1997. B59. P. 319-351.

Поступила в редакцию 31.05.05

(R.^i $W)(w,i)

duj

(йоФН)(ы,0)

duj

УДК 519.214.4 Д. И. Жданов

НОРМАЛЬНАЯ СХОДИМОСТЬ САМОНОРМИРОВАННОЙ СТАТИСТИКИ ПИРСОНА

(кафедра математической статистики факультета ВМиК)

1. Введение и основные обозначения. Пусть производится п независимых наблюдений случайной величины £ и Но — гипотеза, состоящая в том, что результаты наблюдений образуют случайную выборку из генеральной совокупности с функцией распределения G(x), которая полностью определена. Одним из наиболее известных методов проверки статистической гипотезы Но является критерий хи-квадрат. Его суть состоит в следующем. Множество Е всех возможных значений случайной величины £ разбивается на s непересекающихся подмножеств Е\, Е2, ■ ■ ■, Es, и для оценки расхождения между теоретическим (гипотетическим) и эмпирическим распределениями используется статистика s 2

Пирсона xl= Y, — — Pk) 1 гДе Рк = / dG, р^/п — относительная частота наблюдений, принад-

к= 1 Рк " Ек

лежащих множеству Ekl к = 1,..., s. Известно [1, с. 455], что случайный вектор v = (fi,..., vs) имеет полиномиальное распределение, другими словами,

P{Vl = mb...,z/s = ms} = ------р™1 ...р™-

mi\m,2]. ■ ■ .ms\

для любых неотрицательных целых чисел rrii,..., ms, сумма которых равна п. Заметим, что сумма pi,...,ps равна единице. Проверка гипотезы Но сопряжена с вычислением распределения статистики Хп- Вычисления становятся чрезвычайно громоздкими при больших п и s. Поэтому распределение статистики Хп заменяют нормальным распредением. Эта замена оправдана не всегда. При определенных условиях такая замена обоснована теоремой Туманяна [2]. Проверка гипотезы Но может быть