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

РЕАЛИЗАЦИЯ БЫСТРОГО АЛГОРИТМА ПРЕОБРАЗОВАНИЯ КИРХГОФА НА ПРИМЕРЕ БЕССЕЛЕВЫХ ПУЧКОВ

С.А. Балалаев, С.Н. Хонина Самарский государственный аэрокосмический университет, Институт систем обработки изображения РАН

Аннотация

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

Введение

Решение задачи дифракции при непараксиальном распространении света в свободном пространстве остается актуальной проблемой. Основная сложность связана с большой вычислительной трудоемкостью реализации интегральных преобразований в оптике, вытекающей из быстроосциллирую-щего характера ядра этих преобразований. Известны алгоритмы, использующие быстрое преобразование Фурье (БПФ), однако при этом теряется информация о реальном масштабе получаемых распределений. К тому же, в некоторых случаях БПФ применить невозможно, например, при реализации интеграла Кирхгофа.

В данной работе рассмотрены свойства симметрии ядра интеграла Кирхгофа и разработан алгоритм вычисления, значительно сокращающий время расчета по сравнению с реализацией численного интегрирования «в лоб». Эффективность разработанного алгоритма исследована на примере бесселевых пучков.

представляет собой несколько более точную форму интеграла Кирхгофа.

Для программной реализации (2) «в лоб» можно воспользоваться численным интегрированием по методу трапеций:

I = ([ ^ =

X

1 N N ,

4 Sd ^^ [ Рп_1,т-1 + т + К ,т_1 + ^п, т ]

(4)

где ^^ - элементарная площадь области между отчетами плоскости интегрирования.

Чтобы использовать (4) для реализации алгоритмов вычисления оператора (2), исходное и(£,п) и выходное и(х,у,2) поле разбивают по отчетам [2], реальные координаты которых вычисляются по сетке:

= Б/Ы(п _N/2); IX, = Б'/М( _N/2); П = О^(т _ N/2), [у = О'/N( _ N/2),

(5)

1. Описание алгоритма

Метод Кирхгофа решения дифракционных задач состоит в использовании интегральной теоремы, согласно которой значение функции и(х,у,2) на выходе, являющейся решением скалярного уравнения Гельмгольца:

д и(х, у, 2) + д и (х, у, 2) +

дх2

д2и(х, у, 2) к 2

ду2

(1)

д22

+ к2 и (х, у, 2) = 0

можно выразить через значение функции и (£, п, 0) на входе [1]:

и(х, у, 2) = _ П Я и(£,п) £ {¡к _ , (2)

где

г =,1 (х_#)2 + (у-п)2 +:

(3)

к = 2 я IX - волновое число, а X - длина волны светового пучка. Выражение (2) получатся через разложение светового поля по сферическим волнам и

где N х N - размер входной и выходной картины по отчетам; О х О - реальный размер входной и О' х О' - выходной картин.

Воспользуемся формулой численного интегрирования (4) и разбиением входных и выходных массивов сеткой (5) для выражений (2) и (3):

е'кг { 1

^ = и (£ , Пт ) — I ¡к _

(6)

где г = ;

2 О2

и (х , у , 2) =---- х

' 1 2п 4N

N N

[Рп_1т-1 + Рп_1т + Рп,т_1 + ^п,т }

Система (6) является простейшей реализацией непараксиального оператора распространения света на расстояние 2 от источника. Основным недостатком является большая ресурсозатратность алгоритма. Вызвано это тем, что подынтегральная функция е'кг из-за большого значения волнового числа к является быстроосциллирующей и приходиться выбирать очень мелкий шаг интегрирования.

В работе [3] была произведена оценка шага интегрирования, переписывая ее для рассматриваемого случая, получаем:

К <

1я.

2 + Г2 2

Яв N.

(7)

где ЯВ=В/2+В'/2, N - количество разбиений одного периода осцилляций Тт,п (для алгоритма (6) предполагалось N=10, см. рис. 1).

Выделим ядро интегрального преобразования (2):

О.. =-

¡к -1

(8)

и будем использовать для него шаг разбиения не крупнее, чем (7).

Рис. 1. Разбиение быстроосциллирующей функции Ог на области интегрирования, исходя из минимального периода биений Ттп

Необходимое число отсчетов не трудно вывести из (7):

М >

1

Яп

N.

N

Я

2 + г2 2

В

(9)

Выражение (9) показывает, какое число разбиений функции Ог является достаточным в области одного отчета входной функции для ее правильного «прописывания». На рис. 2 представлены дополнительные разбиения функции Ог в области отсчета входной функции.

О о

о о

О о о

о ® о

о о о

Рис. 2. Разбиение матрицы О на дополнительные отсчеты:

(a) соответствует отсутствию разбиения (М=1),

(b) двойному разбиению (М=2), (с) тройному (М=3)

Используя разбиение на дополнительные отчеты (см. рис. 2) для (8) имеем:

Г = 4(Х + х'к-4„ )2 + (У; + У\-Пт )2 + Г2

(10)

У =

В'

NM В NM

(к - (М +1)/2); ((- (М +1)/2).

(11)

г2;

Теперь не трудно переписать (6) используя (8), (10), (11):

1 М М Р&г ( 1 N

Ои,т = Мг ЕЕ VI ¡к - - I, где

М2 £! £1 г2 ^ Г )

г = 7 (х, + хк)2 +(у. + у',-Пт)2 +

Рпт = и (Мп ,Пт )Оп,т ; (12)

и( ) г В2

и (х , у , г) =---- X

' 1 2п 4N

N N

ХЕЕ [ Рп-1, т-1 + Еп-1, т + Еп,т-1 + Еп, т ].

При расчете выходного изображения легко заметить, что функция ОГ периодически повторяется, поэтому можно сократить время выполнения алгоритма (12) при помощи табуляции Опт. Действительно, из рис. 3 видно, что матрица табуляции значений ОГ имеет размерность 2Ы-1, что значительно сокращает время вычисления. Матрица О полностью покрывает все значения матрицы Е, необходимые для подсчета любого и(х,.,у,,г). На рис. 3 показано вычисление центрального значения и(0,0,г), когда суммируются произведения соответствующих элементов заштрихованной области. Для вычисления произвольных значений и(х,уц,г) соответственно матрица Е перемещается по всей области матрицы О. Таким образом, вместо N получаем 4.N операций.

где

Рис. 3. Вычисление выходных значений при помощи перемножения соответствующих значений подынтегральной функции Е на значения табулированной функции О

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

пользуется алгоритм радиальной развертки. В результате получаем оптимизацию вычисления в 8 раз.

различие с помощью среднеквадратического отклонения:

Рис. 4 Радиальная развертка матрицы О

2. Исследование алгоритма Рассмотрим теперь результаты работы алгоритма полученного на основе выражения (12), но реализованного различными способами - через пересчет всех значений (обычный) и с табуляцией. В качестве входной функции и(М,п) возьмем одномодовый бесселевый пучок:

и(г,Ф) = Ст ^т(аг)ехр(,тф), (13)

сформированный ДОЭ размером 5x5 мм, с параметрами моды т=1, а=7, Ст,а=1 (амплитудно-фазовая картина входного распределения показана на рис. 5).

(а) (б)

Рис. 5. Амплитуда (а) и фаза (б) моды Бесселя с параметрами: т=1, а=7, Ст,а=1 во входной плоскости

Длину волны монохроматического лазерного пучка освещающего ДОЭ примем равной 2=633 нм. В первую очередь нас будут интересовать результаты в непараксиальной области - на расстояниях, не превышающих г=70 мм от плоскости ДОЭ. Выберем количество отчетов входных/выходных изображений N=100. Исходя из (9) примем нижний предел числа внутренних разбиений экспоненциальной функции М=23. Для заданного г получим две выходные функции и0(х,у,г) и и(х,у,г), представленные на рис.6. Одна соответствует алгоритму с пересчетом всех значений (верхняя строка), другая алгоритму с табуляцией (нижняя строка). Визуально они ни чем не отличаются, однако мы можем оценить их

8 =

1

{{((0( х, у, г) - |и (х, у, г )|)) х а у

и0(х, у, г)|2а х а у

(14)

где и0(х,у,г) - комплексная амплитуда поля, вычисленная обычным алгоритмом численного интегрирования; и(х,у,г) - комплексная амплитуда поля, полученная оптимизированным методом. Для г=70 мм по формуле (14) получаем 8=10-11, при этом на одной и той же машине оптимизированный без развертки алгоритм рассчитывается в 845 раз быстрее обычного, а с разверткой в 1440 раз быстрее. Причем погрешность между результатами, полученными двумя последними алгоритмами равна машинному нулю.

Далее произведем расчет погрешности вычислений по формуле:

8м =

Ж2М (х, у, г) - рм (х, у, г)|)) х а у

и2М (х, у, г)|2а х а у

(15)

где иМ(х,у,г) - комплексная амплитуда поля, полученная с числом дополнительного разбиения М, и2М(х,у,г) - комплексная амплитуда поля с удвоенным числом разбиений. Для г=70 мм значение (15) получилось равным ¿М=0,05%.

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

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

Из таблицы 1 видно, что погрешность отлична от нуля там, где число дополнительных разбиений М> 1, т.е. операция развертывания матрицы О вносит некоторое (несущественное) отличие в результаты. Таким образом, алгоритм с табуляцией позволяет вычислять интеграл Кирхгофа с той же точностью что и при использовании обычного алгоритма (если не использовать радиальную развертку) с ускорением не менее, чем в 4 раза. При использовании же дополнительных разбиений точность вычислений повышается.

Е

Е

Е

Е

амплитуда

фаза

амплитуда

фаза

амплитуда

фаза

т

Ш

2 = 1 мм

Ж

т

2 = 10 мм

2 = 70 мм

Рис. 6. Комплексное распределение моды Бесселя с параметрами: т=1, а=7, Сща=1 на различных расстояниях от экрана при использовании обычного алгоритма (верхняя строка) и с табуляцией (нижняя строка)

Таблица 1. Сравнение параксиальных и непараксиальных операторов

Расстояние от ДОЭ 2, мм Число отсчетов входной и выходной функции, N Число разбиений экспоненциальной функции, М Время выполнения обычного алгоритма (час:мин:сек) Время выполнения оптимизированного алгоритма (час:мин:сек) СКО между амплитудами обычного и оптимизированного алгоритма 5, % Ускорение в К раз

1000 100 2 00:08:13 00:00:38 10"14 13

70 100 23 16:11:19 00:01:03 10"11 845

70 100 23 16:11:19 00:00:40* 10-11 1440

10 100 112 308:53:46 00:01:50 10"8 6836

1 100 152 - 00:02:59 - -

1000 100 2 00:08:15 00:00:38 10-14 13

700 100 3 00:18:19 00:00:38 10-14 29

500 100 4 00:30:44 00:00:40 10-14 46

1000 130 2 00:23:18 00:01:48 10-14 13

2000 130 1 00:07:13 00:01:47 0 4

3000 130 1 00:07:26 00:01:51 0 4

4000 130 1 00:07:12 00:01:45 0 4

* - помечаются то время, которое получается при использовании радиальной развертки

а

д

литу

л п

((((о))) в а о

530 щ ш щ м

2 = 0 мм 2 = 1000 мм 2 = 2000 мм 2 = 3000 мм 2 = 4000 мм

Рис. 7. Комплексное распределение моды Бесселя с параметрами: т=2, а=7, Ст,а=1 на различных расстояниях от экрана

Заключение

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

При реализации всех предложенных способов оптимизации алгоритма удалось получить значительное ускорение его работы, благодаря чему для отдельно взятой задачи результат можно получить не за неделю работы программы, а за считанные минуты (см. табл. 1) притом, что по-

грешность будет близка к погрешности машинных вычислений.

Благодарности

Работа выполнена при поддержке российско-американской программы «Фундаментальные исследования и высшее образование» (CRDF Project SA-014-02), гранта РФФИ 05-01-96505, а также гранта INTAS 04-77-7198.

Литература

1. Виноградова М.Б., Руденко О.В., Сухоруков А.П. Теория волн. - М.: Наука, 1979. - 383 с.

2. Бахвалов Н.С. Численные методы (анализ, алгебра, обыкновенные дифференциальные уравнения). - М.: Наука, 1975. - 630 с.

3. Дроздов М.А. и Хонина С.Н. Исследование границ применимости параксиального приближения для описания распространения лазерного света в свободном пространстве // Естествознание, экономика, менеджмент, 2003. № 4. С. 47-55.