Проблема восстановления фазы световой волны является одной из классических задач оптики. Обсуждению тех или иных аспектов этой проблемы посвящено большое количество работ /1,2/. В классической постановке задачи фазовый профиль требуется найти по единственному распределению интенсивности. Несмотря на существующие дисперсионные соотношения, связывающие интенсивность и фазу поля (преобразование Гильберта) непосредственно воспользоваться этой связью в практических задачах оказывается довольно сложно /2,3/. По этой причине широкое распространение получил альтернативный подход к восстановлению фазы, приводящий к итерационным методам её восстановления. Рассмотрим этот подход подробнее.

Пусть на вход некоторой оптической системы поступает монохроматическая световая волна с комплексной амплитудой A(,z):

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

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

Одним из первых практических алгоритмов восстановления фазы был алгоритм Гершберга-Сакстона, предложенный в 1972 году /5/. Этот алгоритм является итерационным, т.е. фаза является пределом последовательности оценок фазовых функций , определяемых в процессе итераций по некоторому правилу. Процедура определения фазы в алгоритме Гершберга-Сакстона была предложена, исходя из эвристических соображений. Позднее было показано /6/, что при малых отклонениях начальной фазовой функции итерационного процесса от точного решения , алгоритм эквивалентен градиентной процедуре минимизации функционала ошибки

где - комплексная амплитуда поля в выходной плоскости z=z₀ на n-ой итерации, соответствующая фазовой функции . Более того, в процессе итераций значение функционала J не возрастает /5/. В целом же вопрос о математическом обосновании алгоритма Гершберга-Сакстона, его эффективности остается открытым. Помимо этого, чисто теоретического вопроса, белее важным являются практические потребности как-то модифицировать алгоритм. Дело в том, что алгоритм Гершберга-Саксона хорошо сходится лишь при достаточно малых отклонениях начальной фазы от искомой функции . Кроме того, метод оказался чувствительным к различного рода шумам, неизбежным при измерении интенсивности. В результате наличие сравнительно малых ошибок в определении интенсивности приводит к расходимости алгоритма. Другой, практически важный вопрос связан с самой постановкой задачи восстановления фазы. В практических задачах часто имеется возможность зарегистрировать не два, а, вообще говоря, несколько распределений интенсивности. Алгоритм Гершберга-Сакстона не дает рецепта, как воспользоваться этой дополнительной информацией.

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

Преобразование светового поля в системе будем описывать некоторым линейным оператором G: , где - комплексная амплитуда поля в выходной плоскости прибора. Будем считать, что при распространении волны через оптическую систему ее энергия W сохраняется, т.е. W=(A,A)=, где (f,g) означает скалярное произведение функций: (интегрирование ведется в бесконечных пределах). Легко показать, что из сохранения энергии W следует унитарность оператора: (сопряженный оператор равен обратному). Введенные ограничения на характер преобразования поля в оптической системе характерны для большинства оптических приборов, не содержащих внутренних диафрагм или других поглощающих энергию элементов. В математическом плане задача восстановления фазы сводится к решению нелинейного уравнения:

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

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

где h_n - длина шага метода, - градиент функционала на n-ой итерации. Начальная оценка выбирается произвольным образом. Найдем выражение для градиента функционала J. Для этого рассмотрим приращение ΔJ, вызванное малым изменением фазы . С точностью до членов второго порядка малости относительно приращения поля , получим

тогда вместо (6) запишем . Используя приведенные выше свойства оператора G, найдем:

В силу (1) , поэтому для градиента функционала получаем окончательное выражение:

Одним из наиболее эффективных градиентных методов является метод условного градиента /4/. В нем очередная оценка фазы на n-ой итерации определяется из условия минимума приращения ΔJ (8), который, как нетрудно показать, достигается при v_n+1=φ_n, где φ_n=argΨ⁽ⁿ⁾ - фаза функции Ψ на n-ой итерации.

Рассмотрим некоторые примеры использования предложенного подхода. Пусть оптическая система представляет собой тонкую линзу с фокусным расстоянием F. Пусть также плоскость регистрации расположена на некотором расстоянии z<F от апертуры линзы. Оператор G описывает в этом случае преобразование Френеля:

Рассмотренная выше процедура минимизации функционала J (4) методом условного градиента оказывается в этом случае эквивалентной алгоритму Гершберга-Сакстона. Тем самым рассмотренный подход к фазовой проблеме позволяет дать строгое обоснование этого алгоритма, выявляет его градиентный характер.

В практических задачах часто имеется возможность зарегистрировать сразу несколько распределений интенсивности I_k() (набор изображений), например, в плоскостях, расположенных на расстояниях z=z_k, k=1,...,N_I от входной плоскости системы (линзы). Стандартные алгоритмы восстановления фазы, как отмеча-

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

представляющий собой естественное обобщение функционала (4) и учитывающий отклонение вспомогательного распределения интенсивности от зарегистрированного в каждой из плоскостей регистрации. Не сложно показать, что градиент функционала также определяется выражением (10), если при этом функция имеет вид:

(оператор G_k описывает распространение поля до плоскости регистрации z=z_k).

Приведем некоторые примеры численного исследования рассмотренных алгоритмов. Рассмотрим линзу с радиусом апертуры b и фокусным расстоянием F. Фазовая функция во входной плоскости задавалась в виде

, М=9, (14)

где - полиномы Цернике, которые нормированы так, что . При этом среднеквадратичное отклонение σ фазы от плоской составляет, как не сложно показать, . Явный вид нормированных полиномов Цернике приведен, например, в работе /10/. Аберрационные коэффициенты характеризуют диапазон изменения фазовой функции . Амплитуду волны во входной плоскости будем считать постоянной: , . Все коэффи-

циенты полагались одинаковыми и равными a_m=a, m=1,...,M. Рассматривались три случая: 1 - малый диапазон изменения фазы (а=0,05 λ - слабые аберрации, максимальное по апертуре изменение фазы Δ=maxU-minU=1,06λ), 2 - аберрации средней силы (a=0,15 λ, Δ=3,2λ), 3 - сильные аберрации (a=0,25λ, Δ=5,3λ). Распределение интенсивности в плоскостях z=z_k рассчитывалось по формуле (11).

Рис.1. Динамика изменения критерия J_Σ при его минимизации методом условного градиента в итерационном процессе восстановления фазы по набору изображений, зарегистрированных на расстояниях {L_k}:

а = 0,15 λ по набору изображений. Рисунок отражает динамику изменения критерия J_Σ при его минимизации методом условного градиента. Расстояние до плоскости регистрации z_k удобно характеризовать безразмерными величинами, L_k, связанными с z_k - соотношениями:

где L_g=0,5kb² - дифракционная длина для заданного размера апертуры линзы b. При этом фокусу линзы z=F соответствует значение L=∞. Из рисунка видно, что использование всего двух изображений (L={0,1, ∞}, кривая 3) позволяет существенно ускорить сходимость процедуры по сравнению с восстановлением фазы по единственному распределению интенсивности (L=∞, кривая 1 и L = 0,1, кривая 2). Добавление третьей плоскости регистрации уже слабо улучшает эффективность восстановления.

Обработка набора из двух изображений (L={0,1 ,∞}) позволяет существенно расширить и диапазон восстановления. Теперь удается восстанавливать фазу, обладающую сравнительно сильными аберрациями (a=0,25 λ). При этом за 75 итераций критерий J_Σ уменьшается почти на 5 порядков (точность восстановления, оцениваемая как максимальное по апертуре отличие оценки фазы от ее точного значения составляет 0,0045λ). Отметим, что в традиционных схемах восстановления по одной плоскости регистрации (типа алгоритма Гершберга-Сакстона) удается восстанавливать фазу, обладающую значительно меньшими аберрациями (a=0,1λ для плоскости регистрации, лежащей в фокусе линзы L=∞, и a=0,15 λ для плоскости с L=0,1).

Вместе с тем, эффективность восстановления фазы и по двум распределениям интенсивности также может оказаться невысокой, если плоскости регистрации выбраны неудачно (рис.1, кривая 4). Это обстоятельство поясняет рис.2, где приведены распределения интенсивности для входного поля с сильными аберрациями (а=0,25 λ) волнового фронта, зарегистрированные на различием расстоянии L_k от апертуры линзы: L₁=0,01,

L₂=0,1 и L₃=∞ (фокус линзы). Как видно из рисунка, распределение интенсивности наиболее чувствительно к фазе входного поля в фокальной плоскости (L=∞). При сильных аберрациях это распределение имеет сложный нелокальный характер и резко отличается от соответствующего распределения при отсутствии аберраций. Это, по-видимому, и приводит к отсутствию сходимости итерационной процедуры восстановления. Использование еще одной плоскости регистрации, лежащей в зоне слабой дифракции ( L=0,01) не изменяет ситуацию, поскольку распределение интенсивности в этой плоскости слабо зависит от фазы входного поля. Улучшение сходимости может наступить, если вторая плоскость регистрации будет расположена в зоне умеренного проявления дифракционных эффектов (L=0,1), что и наблюдалось в численных примерах (см.рис.1, кривая 3). По мере приближения оценки фазы к ее точному значению , основной вклад в восстановление фазы начинает вносить фокальная плоскость регистрации, которая и обеспечивает необходимую точность.

Использование нескольких плоскостей регистрации позволяет также значительно повысить устойчивость процедуры восстановления по отношению к аддитивному шуму в зарегистрированной интенсивности. Это связано, по-видимому, с некоторым эффективным усреднением шума, которое происходит при обработке набора изображений. В таблицу сведены результаты восстановления волнового фронта при наличии шума в каждом из зарегистрированных распределений интенсивности , k=1,...,N₁, где - равномерно-распределенная на интервале [0,1] случайная величина: , β - отношение сигнал/шум, - неискаженное изображение.

На практике большой интерес представляет задача определения некоторого конечного числа M аберрационных коэффициентов волнового фронта a_m, m=1,...,N₁ непосредственно по распределению интенсивности (или набору распределений ) без предварительного восстановления самой фазы . В таких задачах использование градиентных методов

слабых аберраций волнового фронта (а=0,1 λ)

*) - приведены результаты для слабых аберраций волнового фронта (a=0.05λ).

дает особо значительные преимущества. В численных расчетах подлежащая анализу фаза выбиралась аналогично (14). Для определения аберрационных коэффициентов a_m использовался градиентный алгоритм:

где - оценка аберрационных коэффициентов на n-ом шаге итерационного процесса, - проекции градиента функционала J_Σ относительно функций , которые, как не сложно показать, анализируя изменение ΔJ_Σ, вызванное малым изменением Δа_m, определяются выражением:

На рис.3 приведена динамика изменения критерия J_Σ в процессе итераций для аберраций различной силы. Определение коэффициентов

Рис.3. Динамика сходимости итерационного процесса определения аберрационных коэффициентов волнового фронта по набору изображений L={0,1, ∞):

1 – a=0,1λ, 2 – a=0,20λ, 3 – a=0,25λ.

а_m производилось по набору из двух изображений (L={0.1, ∞}). Заметим, что теперь изменение критерия J_Σ (особенно в начале итерационного процесса) происходит немонотонным образом, что не наблюдалось при восстановлении всей фазы (см.рис.1).

Следует отметить также, что предлагаемый метод анализа волнового фронта имеет примерно в M/2 раз меньше вычислительные затраты на один итерационный шаг и в 5 раз больший диапазон восстановления, чем итерационный алгоритм, предложенный в /10,11/.

Все приведенные в работе результаты по восстановлению фазы и определению аберрационных коэффициентов волнового фронта получены при одной и той же начальной оценке фазы . Коэффициенты полагались равными 0,01λ (максимальный размах такой фазы по апертуре равен 0,06λ, т.е. она слабо отличается от плоского волнового фронта).

1. Аблеков В.К., Колядин С.А., Фролов А.В. Высокоразрешающие оптические системы. М., Машиностроение, 1985, с.176.