УДК 535.42
ГРАДИЕНТНЫИ МЕТОД ОПТИМИЗАЦИИ В ЗАДАЧЕ СИНТЕЗА БИНАРНОЙ МИКРООПТИКИ
© 2001 В.В. Котляр1, Д.В. Нестеренко2
'Институт систем обработки изображений РАН, г. Самара 2Самарский государственный аэрокосмический университет
Рассмотрено применение гибридного метода граничных - конечных элементов для синтеза двумерных дифракционных оптических элементов (ДОЭ). Гибридный метод позволяет эффективно рассчитывать векторную дифракцию от элементов микрооптики. В статье обсуждается применение градиентного метода оптимизации к матричной записи гибридного метода. Рассмотрен синтез бинарных линз, имеющих размеры порядка длины волны падающего света.
Введение
Когда минимальный характерный размер дифракционных структур приближается к длине волны освещения, необходимо применять векторную модель для анализа дифракции от них. Также векторная модель требуется, когда представляющая интерес область находится вблизи дифракционного элемента и, тем более, внутри него. Аналитические векторные модели дифракции могут быть определены для некоторых периодических структур, для апериодических структур граничные условия на электромагнитные поля делают аналитическую модель невозможной. Следовательно, для моделирования дифракции от апериодических дифракционных оптических элементов (ДОЭ) должны использоваться численные технологии, основанные на решении уравнений Максвелла. В работах [1], [2] и [3] описаны подходы к решению задач ближнего поля в свободном пространстве и в волноводе соответственно. Показано, что при использовании рефракционных микролинз наблюдается смещение максимума интенсивности от геометрического фокуса. При построении дифракционных микролинз из рефракционных наблюдается сдвиг положения фокуса и снижение эффективности дифракционной микролинзы [4]. Невозможно использование рефракционных микролинз в задачах с фокусировкой вне оптической оси вследствие их симметричности. В связи с этим возникает необходимость применения оптимизационных методов для синтеза дифракционных оптических элементов,
способных определить оптимальный профиль для заданного распределения интенсивности, принимая во внимание ограничения на форму рельефа. В данной статье мы представляем градиентный алгоритм для синтеза бинарной микрооптики.
Метод расчета
Используя гибридный метод конечных элементов - граничных элементов [1], рассмотрим дифракцию света в свободном пространстве.
Для неоднородной диэлектрической двумерной среды для двух случаев поляризации можно записать уравнение Гельмгольца для монохроматической электромагнитной волны:
д_
дх
1 ду(х, у)
а
дх
+ -
ду
1 ду(х, у)
а
ду
+
+ к%ру( х, у) = 0,
(1)
где для ГЕ-поляризации у = Е2
а
т
то
2р
1о
- вол-
новое число в вакууме, 10 - длина волны падающего света, ю - циклическая частота излучения; е, е0, т, т0 - диэлектрические и магнитные проницаемости вакуума и среды;
Е, Н - вектора электрической и магнитной напряженности электрического и магнитного полей; для ГМ-поляризации: у = Н 2,
а = £(^), ь = т
£0 т0
Решение уравнения (1) равносильно минимизации функционала, рассчитываемого в ограниченной области О, включающей свою границу 5 [1 ]:
а
1
т( X У)
Уу(х, у)Уу(х, у) -
(2)
ЭР (у50) Эу
к
= о, (х, у) є а
(4)
где у = (у 1? ум); у к =у к + у к, у
к
Эп
производная поля по внеш-
ней нормали на границе 5.
Элементы матрицы А вычисляются по формулам:
ак,] = Я
Эх
Эх
+
+
ЭУ
эу
(6)
- ко2е(х, у)у 2 (х, у)\іхіу - 2 ЯЯ^усИ.
і Эп
Приведем уравнение (2) к системе линейных алгебраических уравнений, используя следующую расчетную схему. Покроем область а треугольной сеткой. После этого в каждом из треугольников построим линейный полином. Для однозначного определения полиномов достаточно задать их значения в вершинах треугольников. Тогда базисом будет система функций юк(х, у). Геометрически функция юк(х, у) представляется как треугольная пирамида с центром в точке
(х„ ук).
Представим поле у через данный базис как
- коекюкк'1 (x, у)ю] к(x, у)
іхіу,
к,] = 1, ..., N + И, где &к . - треугольная область, содержащая к иу узлы сетки.
Элементы матрицы В вычисляются по формулам:
Яю^т’ р^т-р ь = - Я т р а,
т,р і 7
т, р
(7)
т, р = 1, ..., М.
Получаем систему линейных уравнений, которую необходимо дополнить соотношением для поля и его производных на границе [6]:
N
у (х, у) = Ъукю<а (х, у), (3) Я
к=1
где к = 1, ..., N + И; N - число внутренних узлов сетки, И - число граничных узлов, приходим к системе линейных алгебраических уравнений:
Подставляя разложение (3) в функционал (2) и применяя разделение общего поля на сумму падающего и рассеянного полей, приравниваем нулю производные
Эу- ^ 5С Э02
------02 -у —2
Эп Эп
41 = - 1у*с,(х, у) є і ,(8)
где функция Грина для двумерных световых полей (цилиндрическая волна) равна функции Ханкеля второго рода нулевого порядка:
в( х х\ y, у О = - -4 н 02) ^к ^( х - х,)2 +(у - у,)2 ^ .(9)
Равенство (8) может быть представлено в матричном виде:
[С] Фе +[Б] V 5е = 0, (10)
где элементы матриц [С] и [Б]:
и приходим к системе линейных алгебраических уравнений [5]:
А(у а + у а) + в(уіс+уіп) = о, (5)
""т,р
X
у к - полное, неизвестное рассеянное и известное падающее световые поля; V = (у1, .,
дУк
- кЯю1 (р + [Рр+1 -Рр ]Х)
сх-
кІюР (-1 + [Рр-Рр-1]х)х (11)
Э^(Рт > Рр + [Рр+1 -Рр ]Х) Эп'
хЭв(Рт , рр-1 + [Рр -Рр-1]Х) сх + X , іх +
Эп' ь
+ у 8
т тр
ак, 1
ст,р = (рр + [Рр+і -Ррх)х
Х в(Рт , Рр + [Рр+і -Рр ХХ + +кІК (р-і + [Рр- Рр-іх )х хв(Рт, Рр-і + [рр -рр-1х¥х
(12)
т, р = 1, М,
где п’ - вектор внутренней нормали; рт = (хт, ут) е ^ - координаты т -ой точки; И - шаг сетки; параметр ут=фт/(2л), где фт - внутренний угол в точке (хт, ут) границы
После линеаризации уравнения (8) получаем полную систему линейных алгебраических уравнений вида:
[А в,в] [А 5,в] 0
[А в,5 ] [А 5,5 ] [в]
о [с] [в]
О
о
О
. БС ф в
, БС
, т
фв
, іп
(13)
[а в ,в] [а 5,в ] 0 ' , БС фв
т = [ в ,5 ] [ 5,5 ] [В] , фБС = , БС
_ 0 [с] [в]_ VБС
и
[А в,в ] [А 5,в ] 0 іп ф в
[Ав,5 ] [А5,5 ] [В] , ф іп = іп іп
ООО V
После нахождения поля уБС на границе 5 поле в любой точке окружающего пространства находится по формуле:
у(х, у) = у1П (х, у) -|
С1,
Эи Эи
(х, у) їО (і6)
В матричном виде поле в окружающем пространстве представляется в виде:
у = уіп + W уБС = уіп + WT-lU уіп. (17)
где матрица W размером Іх(^ + 2М) равна:
W = [о [е] [е]],
где элементы подматриц Е и Е:
Єр (х, у) = А ЦЮр (р + [Рр+і - Рр х )х Эв((х,у),Рр + [Рр+і -РрX)
х
Эп'
+ Кк(-і + [Рр - Рр-іХ)х (і8)
Эв((х,у),Рр-і + [Рр -Рр-іХ)
х
Эп'
где в - область О без границы Б, подматрица Ав в размерностью NN включает коэффициенты вклада внутренних узлов, подматрицы Ав5 и А5в размерностью NxM и MxN включают коэффициенты вклада внутренних узлов с внешними, подматрица А55 размерностью МхМ включает коэффициенты вклада внешних узлов, у в и у5 - вектора значений поля во внутренних и граничных узлах сетки.
Пред ставим систему уравнений (і3) в матричном виде:
туБС = и уіп, (і4)
у бс = т-іиуіп, (і5)
где
/р (х, у) = - Юр (рр + [р^+1 - Рр X )х
х в((X, у), Рр + [Рр+1 - Рр X -
-й/0ю^ррр + [рр -Рр-1#)х (19)
хв((х,у),рр-1 + [рр -Рр-1%№, р = 1, ..., М.
Метод оптимизации
Определение оптимального профиля состоит в расчете параметров (у1, ..ук) профиля бинарной микрооптики (рис.1). Это достигается разработкой алгоритма систематического поиска оптимального решения. В качестве метода оптимизации мы выбрали градиентный метод.
Для построения процедуры расчета введем функцию ошибки е (р), характеризующую отличие рассчитанных значений интенсивности I от требуемых значений I’ [7]:
е (р) = е (1(р), 1’),
где р = (ур ук) - вектор параметров бинар-
ного профиля; к - число параметров профи-
ля; (Іп) ММ, (I’) ММ - вектора рассчитанных и
требуемых значений интенсивности.
Градиентная процедура минимизации функции ошибки е (р) состоит в итерационной коррекции параметров профиля микрооптики по правилу:
Рп+і = Рп - ї -Vе СРпХ
где п - номер итерации, ї - шаг градиентного алгоритма,
V е (р) =
Эе(Р) Эе(Р)
ЭУl ,"', Эук
ные матрицы; у (р) - комплексный вектор. Из (і5) следует:
ЭфБС (Р) = т -1 эи(Р). п + Эт-‘(Р) иф п
Эр Эр1 'Эр1 ,
отсюда:
эт-і(Р) _ Эл-
Эф БС (Р) - т -і Эи(Р)
м эр
(иФіп) -і.(2О)
Из (і4) следует:
Э(т(Р)ф БС (Р)) = т Эф (Р) , Эт(Р)
Эр-
+-------- —~ Ф
Эр Эрі
Эи(Р)
Эр1
Фгп
>
отсюда ЭфБС (Р) Эр-
= т
-і
Эи(Р) ф іп Эт(Р) фБ,
Эр, ф Эр, ф
Подставляя (2і) в (2О), получим: эт-і(Р) „-і Эт(Р)
. (2і)
ФБС (иф іп ) -і. (22)
градиент функции ошибки.
Точно определить шаг t не представляется возможным, т.к. неизвестно значение минимума функции ошибки е (р). Поэтому после вычисления градиента функции ошибки V е (р) из исходной точки рп делается шаг заданной длины, находятся координаты новой точки р . Если значение е (ри+1) лучше, чем е (рп), то шаг считается удачным, и новая точка становится исходной, иначе - шаг неудачный, и его величина уменьшается. Если величина шага становится меньше определенного предела, то алгоритм останавливается, и последняя найденная точка - решение, иначе процедура уменьшения шага повторяется.
Вычисление градиента функции ошибки V е (р) может вестись двумя способами. Первый заключается в непосредственном
де(р)
нахождении частных производных ^х в
разностном виде.
Второй способ состоит в представлении е (р) через вектор электрического поля у (р), учитывая зависимость I = у у *.
Электрическое поле у определяется из системы (17), где Т(р), и(р), W - комплекс-
эр
Найдем выражение для производной поля в любой точке окружающего пространства по параметру Р из системы (і7), используя (22):
Эф (Р) = Э( \ут-і(Р)Ц(Р)) фіп =
^і эр
= W
Эт Чр) и + т-і Эи(Р) эр ЭЛ
Фгп
= WT
.-і Эт(Р) ф БС + Эи(Р) ф іп
_ эр эр
Учитывая (і3) можно записать:
(23)
Эф (Р) =
Фі
= WT -і
.-і ЭА11 (Р)
ЭА1,1 (Р) , іп ЭА1,/ (Р) , .
^------ф-----Э-----ф
Фі Эр
І % іп І БС І
1ф - ф }
(24)
Элементы вектора градиента , где
Фі
рі - компонента вектора Р, имеют вид: іО7
Эе(Р) =^Эе(1, I’) Э17- (Р)
Эрі
Э1,
Фі
Э/7.(р) = ( Эф*(р) + * Эф7.(р) (25)
= ф 3 ~Э + ф 3 Э . (25)
Ф/ Ф/ Ф/
Для определения V е (р) уравнение (24) подставляется в уравнение (25).
Отметим, что данный подход позволяет определять любые другие параметры, подлежащие оптимизации.
Проведен ряд экспериментов, показывающий эффективность использования градиентного метода оптимизации бинарных микролинз.
Конструкция дифракционной линзы
Рассмотрим цилиндрическую рефракционную микролинзу (рис.2).
Положим, что ее длина Ьх настолько больше ширины Ь, что длину линзы по сравнению с шириной можно считать бесконечной. Рассмотрим численный пример. Определим апертуру линзы как а = 8 мкм, радиус кривизны Я = 5 мкм, показатель преломления п = 2 для длины волны 10 = 1 мкм, что примерно соответствует хлориду серебра (пА§С1 = 2,02239 для 10 = 1 мкм).
Для данной цилиндрической линзы построим линзу с фазовой функцией, приведенной по модулю 2р, и дифракционную линзу с двумя степенями градации фазовой функции (рис.3).
3 МКМ
Рис. 2. Цилиндрическая рефракционная микролинза
1
0.5
1
0,5
Хмкм
У мкм
У мкм
У мкм
0 4
Рис. 3. Построение линзы с фазовой функцией, приведенной по модулю 2Р, и бинарной дифракционной линзы
Численное моделирование и оптимизация
Поставим вычислительный эксперимент, в котором плоская волна ТЕ-поляриза-ции падает на бинарные линзы в свободном пространстве. В качестве оптимизируемых параметров профиля брались координаты выступов.
На рис.4 показано распределение интенсивности дифракции плоской волны на построенной однофокусной бинарной линзе (рис.3) с дифракционной эффективностью 54,3 % и энергетической эффективностью 43,5 %. Ее профиль был использован в качестве начального для процедуры оптимизации и был оптимизирован до дифракционной эффективности 72,3 % и энергетической эффективности 64,0 %.
Рис. 4. Распределение интенсивности света при дифракции на бинарной линзе начального построения
-4-3-2-10 1 2 3 4
у, мкм
Рис. 6. Распределение интенсивности света в фокальной плоскости бинарных линз: с начальным профилем и оптимизированным
Рис. 5. Распределение интенсивности света при дифракции на оптимизированной бинарной линзе
Под энергетической эффективностью мы понимаем отношение энергии, попавшей в область фокуса, к энергии, падающей на микролинзу. Дифракционная эффективность - это отношение энергии, попавшей в область фокуса, к энергии, прошедшей микролинзу. Интенсивность - это квадрат модуля проекции электрического вектора на ось г.
Рассмотрим синтез бинарной линзы, дающей в плоскости, располагающейся на расстоянии 3,5 мкм, три максимума равной интенсивности шириной і мкм: один в центре и два на расстоянии 2,5 мкм от центра. Были использованы несколько начальных профилей линзы с тремя выступами, располагающимися симметрично относительно
оси симметрии, для нахождения по возможности всех локальных экстремумов функции ошибки. В результате работы процедуры оптимизации было найдено два вида профиля, дающих максимальные значения эффективности. На рис.7 показано распределение интенсивности дифракции плоской волны на оптимизированных трехфокусных бинарных линзах. Линза, показанная на рис.7а, дает дифракционную эффективность 86,2 % и энергетическую эффективность 79,3 %. Линза, показанная на рис.7б, дает дифракционную эффективность 85,8 % и энергетическую эффективность 84,1 %. Дифракционная эффективность в этом случае равна отношению энергии в трех заданных порядках к энергии, прошедшей микролинзу. Под энергией понимается сумма интенсивностей на выбранной линии вдоль оси у.
а) б)
Рис. 7. Распределение интенсивности света при дифракции на оптимизированных трехфокусных бинарных линзах
б
z, мкм
Рис. 8. Распределение интенсивности света при дифракции на оптимизированных трехфокусных
бинарных линзах в указанной плоскости
Таким образом, представленный алгоритм позволил найти профиль трехфокусной бинарной линзы (рис.7б) с высокой дифракционной эффективностью и глобальным максимумом энергетической эффективности, более простой в физической реализации, чем линзы с начальным профилем. Линза, показанная на рис.7а, с глобальным максимумом дифракционной эффективности дает форму поля дифракции наиболее соответствующую заданной.
Выводы
В работе рассмотрен градиентный метод расчета профиля диэлектрического микрообъекта с использованием совместного метода конечных - граничных элементов. Работоспособность метода продемонстрирована на примере расчета бинарных однофокусной и 7. трехфокусной микролинз.
СПИСОК ЛИТЕРАТУРЫ
1. Котляр В.В., Нестеренко Д.В. Анализ задачи дифракции света на микрооптике гибридным методом конечных элементов
- граничных элементов // Компьютерная оптика. 2000. Вып.20.
2. Головашкин Д.Л., Сойфер В.А. Анализ прохождения электромагнитного излучения через дифракционную линзу // Автометрия. 1999. №6.
3. Kotlyar V. V., Nesterenko D. V. Analyzing the light diffraction by binary micro-optics using a combination of boundary element method and finite element method / Laser physics and spectroscopy, Saratov fall meeting, Proceedings of SPIE. Vol. 4003. 2000.
4. Головашкин Д.Л., Котляр В.В., Нестеренко Д.В. Анализ дифракции света на микролинзах в свободном пространстве и волноводе // Компьютерная оптика. 2001. Вып.21.
5. Kotlyar V.V., Nesterenko D.V. Modeling the light diffraction by microoptics elements using the finite element method / Laser physics and spectroscopy, Saratov fall meeting, Proceedings of SPIE. Vol. 4002. 1999.
6. Бреббия К., Теллес Ж., Вроубел Л. Методы граничных элементов. М.: Мир, 1987.
Методы компьютерной оптики. Под ред. В.А. Сойфера. М.: Физматлит, 2000.
SYNTESIS OF BYNARY MICROOPTICS BY GRADIENT METHOD
OF OPTIMIZATION
© 2001 V.V. Kotlyar1, D.V. Nesterenko2
1 Image Processing Systems Institute of Russian Academy of Science, Samara
2 Samara State Aerospace University
We apply a gradient method as synthesis algorithm to the two-dimensional hybrid finite element-boundary element method as vector model of diffraction for the design of DOEs. The hybrid method is capable of modeling inhomogeneous DOEs in unbounded freespace in computationally efficient manner. In this paper we discuss the application of gradient method of optimization to matrix notation of hybrid method. In such application it is possible to analyze profiles with large number of features. This allows to overcome the limitations of time of calculation dependent on amount of modifications of DOE. We use the gradient method to design of binary-phase lenses that has subwavelength features. Although we have considered only binary-phase lenses, presented gradient method is capable of design inhomogeneous DOEs or coatings.