Обработка данных методом преобразования значений производных функций на сетке в коэффициенты Фурье Текст научной статьи по специальности «Математика»

УДК 519.6, 519.612, 519.651

Обработка данных методом преобразования значений производных функций на сетке в коэффициенты Фурье

Л. А. Севастьянов, М. Г. Кокотчикова, Д. С. Кулябов

Кафедра систем телекоммуникаций Российский университет дружбы народов ул. Миклухо-Маклая, д. 6, Москва, Россия, 117198

В работе рассмотрены две практические задачи восстановления функции на круге по возмущенным значениям функции и её частных производных на дискретной сетке. Предложены формулировки устойчивого восстановления функции методом регуляризо-ванного разложения в ряд Фурье по полиномам Цернике. Вычисление стабилизирующих функционалов реализованно в модуле для системы компьютерной алгебры Axiom.

Ключевые слова: восстановление функции по дискретному следу, метод регуляризации, аналитические вычисления.

1. Введение

При разработке математической модели процесса экранируемого вакуумного напыления возникла вспомогательная задача: восстановление коэффициентов Фурье функции эффективного распределения источника напыления по измеренным с погрешностью значениям функции на конечной сетке (не определенной заранее) [1,2].

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

Обе задачи являются неустойчивыми ко входным данным. Следовательно, к ним нужно применить метод регуляризации А.Н.Тихонова [5,6].

1.1. Устойчивое восстановление функции по её следу на

сетке

Рассмотрим задачу отыскания вектора коэффициентов Фурье в разложении по полиномам Цернике } на единичном круге функции f, значения которой известны лишь в точках некоторой сетки Тт = {¿1,..., Ьт} С .

Полиномы Цернике образуют ортонормированный базис: ^^) = в гильбертовом простанстве Ь2 (О), который задает отображение Р Ь2 (О) в гильбертово

\

пространство 12: ^ Е С]Фз(?) = {сз}Г".

\з=1 )

Измеряются значения f в точках }, т.е. оценивается абсолютная величина возмущения функции (tj)| = /(tj) — /(^) . Этому соответствует метрика рс:

Рс(1,1) = таХ ^ — ^ в пространстве С (О) непрерывных на круге Q функций. Решать надо обратную задачу относительно отображения А : 12 ^

} ^ Е Сз <рэ (д). 3 = 1

С (Q) :

Статья поступила в редакцию 25 октября 2008 г.

Работа выполнена при частичной поддержке грантов РФФИ №07-01-00738а, №08-01-00800а.

Данная задача эквивалентна задаче с оператором А : ^ С (О). В силу компактности вложения С С (О) [7], регуляризация последней задачи до-

стигается сужением А на компакт Р$ = | / е Ь (О) : ||/||^2 ^ $} •

1.2. Устойчивое восстановление функции по следу её

градиента на сетке

В задачах обработки гартманограмм результатом наблюдения является вектор приближенных значений частных производных в конечном наборе точек апертуры Тт [3, 4]. Рассматриваемые нами функции являются один раз непрерывно дифференцируемыми / е С1(О) С [(О). Полиномы Цернике образуют полный набор функций в Ш2(О), а отображение Р (Е а,]<£](&)) = {а^ на Ш2 индуцирует

отображение БР (д)) = \ Ет^аЛ С 12 ©I2.

I к к )

Рассматривается отображение Б : С 1(О) ^ С (О) ©С (О), где Б(/) = , ,

так как требуется оценка её частных производных ^ и ^ в метрике пространства непрерывных функций:

д/

дх

с (Я)

тах чеЯ

д/, , дх(д)

д

д

с(Я)

тах чеЯ

д/, , ¥(д)

Задача восстановления выглядит следующим образом

тах чеЯ

тах чеЯ

£

а> ^ (д) - дх («)

д

^(д) - дУ(д)

д

^ тт, аеЬ

^ тт.

аек

Таким образом, требуется для линейного оператора В : 12 ^ С (О) : ©С (О)

В (а)

решить обратную задачу

4 ду

с учётом того факта, что вектор а е и>2 обеспечивает принадлежность Еа^^(д) е С1(О). Метод регуляризации А.Н.Тихонова требует сужения отображения В на компакт [5, 6]. Воспользуемся компактностью вложения [7] Ш2(О) С С1(О) и зададим компактное множество следующим образом:

^ = {/ еС1 : ||/|и2з < 6} • 2. Необходимые сведения

к(1)

Дифференцирование функции / = Е а^^^ либо в полярных, либо в декар-

3 = 1

товых координатах понижает степень полиномиальной функции на единицу. По-к (1-1) к (1-1) этому =2^ Оз'^з и = с^^з, то же самое можно сказать и про и

. Здесь к(1) = (I + 1)(1 + 2)/2 [8].

3 = 1

3 = 1

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

дРз ^ х дРз ^ у

= ЮХкРк; = Т.т%Рк ■

дх

з<1

(1)

з<1

Матрицы {т^к} и |тук} являются треугольными и могут быть построены по формулам

тхк = I РкАхАУ; = I Рк.

Я

Я

(2)

Норма в пространстве №3(0!) задаётся соотношением

2

W°

/Н2+(

Я у у

д/, , дх^

+

д/ , , ду®

+

+И

д2/( ) дХ2{(1)

+ 2

д2

и(д)

+ 3

д х д д3

(а)

д х2 д

( )

+

+

+3

д2/( ) д3

+

д х д 2

( )

+

0 (^

ё д

и индуцирует для функции /(д) = ^2ауРу(д) норму в пространстве и>2 числовых последовательностей, вычисляемую с помощью треугольных матриц тх и ту дифференцирования полиномов Фурье. А именно, первые производные имеют вид

%(<!) = ЕЬХРз(^ %(1) = ЕЬзРз(Я)>

где

ЬХ = Е т°хкак, ЬУу = Е тЗкак

Вторые производные имеют вид

д2/ = ,,ххт

дх2 = ^ Су Рз'

д2

д2

Еху д и _ V""4 уу

дхду С3 Рз' = ^ С3 Рз'

где

схх = ^ тХкТыа1, сУуУ = Е ТукТ^аь сХу = Е тХкТка = Е Т%7ыа1 = с:

3,к

„Ух

~з ■

3,к у,к

Аналогично, третьи производные имеют вид

д х д3/

3,к

д х2 д д3/ д х д 2

д3£ д 3

Т^^аз, <х = ^тХкТыТгта™,

Е^ХУаз, <ху = ^тХкТыТш0-™

Е<уХаз > <ух = ЕТукТыТ^ш

Е^уа3, = у у у / ЛзкТЫТ1таш.

2

2

2

2

2

3

Опираясь на указанные соотношения, сформируем квадратичные формы следующим образом

П0О =[5го ],ьз = !,•••, к(1) = (1 + 1)(1 + 2)/2,

П1 =тх(тх )Т + ту (ту )т,

П2 =тхтх(тхт Х)Т + т уту (туту )Т + (тхту (тхтУ )Т + т утх(тутх)т), п3 =тхтхт х(т хтхтх)т + тут уту (тутут у )т + т хт хту (тхт хту )т+ + т хтутх(тхт утх )т + тутхт х(т утхтх)т + тхт уту (тхтуту )т+ + т Утхту (тУт хту )т + тутут х(т УтУтх)Т •

С их помощью нормы в пространствах числовых последовательностей {а^ } = а е 12 = ^2,^2,^2 и ю3 задаются формулами

(3)

N1= Еап%ак, N11} = Еа,п)как + N1, N1 ^2 = 2 (п2(а,(а) + ||а11

||а|| 1з = (П0а, а) + (П1а, а) + 2 (П2а, а) + 2 (П3а, а) . 2 2 6

3. Постановка задачи

|

3.1. Задача 1. Устойчивое восстановление функции по её

следу на сетке

Задача устойчивого восстановления вектора коэффициентов f формулируется следующим образом

N

а,1 = тах I Е (д) -

J яеЯ

^ тт •

(4)

Задача (4) эквивалентна задаче минимизации сглаживающего функционала

Мс

= N

+ а ||а11 ш2 ^ т]п

(5)

с параметром а, согласованным с уровнем погрешности ё задания /.

Для любого а > 0 существует единственный вектор аа, доставляющий минимум функционала (5).

Уравнение Эйлера задачи (5) имеет вид:

Е(8Ц + аро ¿уП«) а°? = сг.

(6)

где с\ = ((£{, /) — коэффициенты разложения f по ортонормированному в \ / ¿2 Ь2(О) базису {(}.

Дискретный аналог задачи (4) формулируется следующим образом

N[а,тт(/)] = тах !Еа((дк) - /(дк)

Чк €Т,

^ тт

аеР5

Дискретизация уравнения (6) приводит к уравнению

к(1) / 2 1 \

+ а Е а? = 2

к \ „а _ хт

г = 1,...,к(1).

(7)

Здесь с™ — дискретное преобразование функции f с сетки Тт в коэффициенты Фурье по полиномам Цернике.

3.2. Задача 2. Устойчивое восстановление функции по следу её градиента на сетке

Задача устойчивого восстановления вектора коэффициентов а$ по измеренным с точностью 5 частным производным С^ , ^функции / = ЕазРз формулируется следующим образом

а;

Здесь

т

д1д1

' дх1 ду

Е<

д£ д£

д х д

д£±_ д! д х д х

^ Ш1П .

(8)

+

с

Е'

дРз д/

д д

с

Задача (8) эквивалентна задаче минимизации сглаживающего функционала

М3

а;

д£ д£

д х д

= N1

а;

д_[_ д£ д х д

+ /з \\а\\;

,з —► Ш1П

(9)

с параметром @, согласованным с уровнем погрешности 5. Уравнение Эйлера задачи (9) имеет вид

е( % + ^ Е Е

з \ г=0 /к

Е( ^5: ^ Е

з \ г=0 /к

13ка3 =

у 3 _ ~У

7Зкак = ^ ,

(10)

где с? = Ц) = щ)ь

^ > — коэффициенты разложения 'у по ор-

тонормированному в базису {рз}.

Дискретизация уравнения (10) имеет вид

Щ) / 3 \ к(1 -1

Е % + ^Е^Гз) е

3 = 1 \ г=0 ' ) к=1

к(1) / 3 , \ к(1 -1

Е 6и + /1Е¿4- Е

3 = 1 \ г=0 ) к=1

у у

Гукак = Су,

г = 1,...,к(1), ъ = 1,...,к(1).

(11)

7зкак = с:

!

=

к(1) / з 1

=1

4. Алгоритм аналитического построения матриц стабилизирующего функционала

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

Испорченная оптическая поверхность содержит аберрации такого типа как, астигматизм, кома и др. Полиномы Цернике разного порядка описывают один из типов аберраций (например, полиномы с порядковым номером 5 и 6 — описывает аберрацию типа кома). Если же поверхность содержит несколько типов искажений, то её аналитическая функция представляет собой сумму различный полиномов Цернике с определенными коэффициентам.

При восстановлении функции, приближенно заданной на дискретной сетке или с приближенно заданным градиентом на сетке, возникает нестабильность, для

устранении которой необходимо построить матрицы стабилизирующего функционала По, fii, 0,2 и Пэ, как указано выше. Для построения этих матриц был разработан алгоритм, реализованный в системе компьютерной алгебры Axiom. Алгоритм описан ниже пошагово и с пояснениями.

Входные данные: N — максимальная степень полиномов Цернике, через

(N+1)(N+2)

которую определяется количество полиномов — -—1—2——.

1. Сначала вычисляется радиальная составляющая полиномов Цернике:

п—т

Rm( ) = _( — — _nn-2s

Rn \Р) — Z^ s| ( п+т — ^ п-т — p ,

где п, m £ {N U 0}, m < п, (п — m) — чётное число. Для вычисления коэффициентов Rnm(p) существует рекуррентная формула:

( п + m)R?(p) — 2npRZ-1(p) + (п — m)R?-2(p) = 0. (12)

Данная формула применима при любых п > 1 и т > 1. Использование рекуррентной формулы (12) позволяет существенно сократить время вычисления.

2. На следующем шаге вычисляются полиномы Цернике:

Z0 = Vu+1R°n(p), m = 0, (13)

ZHe4j = V п + 1R™(p)V2sinme, m = 0, (14)

= VaTlR™(p)V2cosme, m = 0, (15)

где п, m £ {N U 0}, m < п, (п — m) — чётное число. Полиномы задаются тремя формулами — для m = 0 (13), а для m > 0 полиномы задаются двумя формулами (14) и (15).

Обычно двумерные полиномы нумеруют двумя индексами. Однако при вычислениях удобнее пользоваться одноиндексной нумерацией. Соответствующее определение было введено в статье [8].

Полиномы Цернике при одноиндексной нумерации задаются следующими соотношениями:

{у/пj + 1 R™*(p)cosmjp, (—1) > 0, mj > 0, л/п-ТГН^з (p)sinmJp, (—1У < 0, m0 > 0, (16)

vO+TR^ (p), mj = 0,

При этом перевод мультииндексов п и m по индексу осуществляются:

■ Г, n (к + 1)( к + 2) I

щ = mini к > 0 Л-^-> j

2 (17)

mj = min{k > 0 : (—1)k = ( — 1)n?, к > 1 —j— п (п + ^ },

0 ^ т ^ п; (п — т) — чётное число. Для заданной максимальной степени полиномов N индекс сквозной нумерации, изменяется от нуля до (N + 1)(N + 2)/2 (число полиномов). 3. Далее необходимо разложить производные по ж и по у от базисных функций (1) в базисе полиномов Цернике.

В силу ортогональности системы полиномов Цернике коэффициенты разложения (1) выражаются через скалярные произведения (2). Скалярные произведения вычислялись следующим образом:

— аналитическое дифференцирование базисных функций на первом этапе;

— процедура scalar_product_circle1 нахождения скалярного произведения для двух заданных функций вычисляет скалярное произведение повторным интегралом, пределы интегрирования в круге радиуса 1: X = -1,..., 1, у = —V1 — х2,... 1 — х2.

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

4. После этого заполняются вспомогательные матрицы 7х, 7у разложения производных от полиномов Цернике по полиномам Цернике (2).

5. В итоге матрицы стабилизирующего функционала П0, П1, П2 и П3 через вспомогательные матрицы 7х и 7х задаются формулами (3).

Выходные данные: В качестве результата получаем матрицы стабилизирующего функционала П0, П1, П2, П3 (3), которые необходимы в задаче восстановления (11). В качестве результата представлена матрица П1 (рис. 1), полученная при использовании разработанного алгоритма для Axiom.

5. Заключение

Итак, решение задачи 1 реализуется решением системы линейных алгебраических уравнений (7) для вычисления вектора коэффициентов аа по заданному вектору с:

Аа о аа = с. (18)

В свою очередь решение задачи 2 реализуется решением системы (вдвое большей размерности) линейных алгебраических уравнений (11):

Г В33 о а3 = С1,

ио аз = (19)

В работе [9] рассмотрен вопрос об устойчивом дискретном преобразовании следа функции на сетке Т (/), измеренного с точностью ¿1, в вектор коэффициентов Фурье с:

С5 = Б- оТ(/).

Решаем один раз с а(51) систему (18) и получаем окончательный результат:

аа = ^1 оТ(¡), Л = А-(1) о Б-

с матрицей для серийной обработки экспериментальных данных Т(/) на сетке Тт (фиксированной), измеренных с фиксированной точностью ¿1.

В работе [9] рассмотрен вопрос об устойчивом дискретном преобразовании следа градиента функции на сетке БТ(/) (в два раза больше данных, чем число точек сетки), измеренного с точностью ё2, в вектор коэффициентов Фурье

(с1,<2)Т:

(С1,С2)1 = ББ- о БТ(/)

(количество строк в матрице ББ- в два раза больше числа точек сетки, а количество строк в матрице Б- равно числу точек сетки).

Решаем один раз с 0(52) систему (19) и получаем окончательный результат:

аз = ББ2 о БТ(/), ББ2 = ^ о БББ-

с матрицей ББ2 для серийной обработки экспериментальных БТ(/) данных на сетке Тт (фиксированной), измеренных с фиксированной точностью 82.

оооооооооооооооооооо

ооооооооооооооооооо^о

ооооооооо_<'ооооооооюоо

ОООООООО ооооооооюооо

оо> оооо ?

4 00

о ? ооооооооооооо

о > оооо

оооооооооооооо

оооооооооооооо^оооооо

ооооооооооооо^ооооооо

ооооо"\ооооооооооооооо

о о о о ^ о

ооооооооооооооо

000^000000010000000000

ооооооооо,^, ОООООООО _<'оо

ОООООООО^, ОООООООО < ООО

О О > ООО

О^ОООООООО £ оооо

ОООООООО < ооооо

о > оооо

ооооо^оооооо^оооооооо

о к л

4

ф

н

5 V 2

4 а а

5

к

®

с

ф

н

0

'В

1

ю

о ч

а о

оооо^оооооо^ооооооооо

ООО ,<,000000^0000000000

о о ^ о о о

о > о

ооооооо > оооо

о^оооо> оооооооо> ооооо

ООО

о к

5

ч о

с

«

ч ч

с

о

к

Рч

о о о о о о

оооооооооооо

с

Литература

1. Жидков Е. П., Севастьянов Л. А. Макропараметры эффективного распределения в математической модели экранируемого напыления // Математическое моделирование. — 1998. — Т. 10, № 9.

2. Sevastianov L. A., Zhidkov E. P. Analysis of Problems in Mathematical Model for Shadowed Sputtering // Comp. Phys. Comm. — 2000. — Vol. 130.

3. Войцехович В. В. Влияние атмосферной турбулентности на точность определения параметров волнового фронта // Препринт ИКИ АН СССР, Пр-862. — 1984.

4. Севастьянов Л. А. и др.. Програмно-математический комплекс для обработки киногартманограмм // Препринт ИКИ АН СССР, Пр-1209. — 1987.

5. Тихонов А. Н. О решении некорректно поставленных задач и методе регуляризации // ДАН СССР. — 1963. — Т. 151, № 3. — С. 501-504.

6. Тихонов А. Н. О регуляризации некорректно поставленных задач // ДАН СССР. — 1963. — Т. 153, № 1. — С. 49-52.

7. Соболев С. А. Избранные вопросы теории функциональных пространств и обо-щенных функций. — М.: Наука, 1989.

8. Noll R. J. Zernike Polynomials and Atmosheric Turbulence // JOSA. — 1976. — Vol. 66, No 2.

9. Севастьянов Л. А., Ловецкий К., Кокотчикова М. Г. Дискретное преобразование значений функций на сетке в коэффициенты полиномов Цернике // Вестник РУДН, Серия «Математика. Информатика. Физика». — 2007. — № 34.

UDC 519.6, 519.612, 519.651 Data Processing by Method of Transformation of Functions and its Derivatives Values on Grids into Fourier Coefficients

L. A. Sevastianov, M. G. Kokotchikova, D.S. Kulyabov

Telecommunication Systems Department Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., Moscow, Russia, 117198

In this paper two practical problems are considered: function reconstruction from its or its partial derivatives perturbed values on a discrete grid. Formulation of stable reconstruction of a function via method of regularized expansion in Fourier series of Zernike polynomials is proposed. Calculations of stabilizing functionals are realized in computer algebra system Axiom.