Численное исследование центроидирования изображения оптического излучения в плоскости регистрации

DOI:

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

При измерении искажений излучения по оценкам, полученным датчиком Шэка – Гартмана в адаптивной оптической системе [8], и стабилизации излучения на входной апертуре в случае следящей адаптивной оптической системы [13] оценивание светового поля выполняется вычислением координат энергетического центра тяжести фокального пятна (центроида) в плоскости регистрации посредством алгоритма центроидирования (the centroiding algorithm). Смещения центроида в двух ортогональных направлениях пропорциональны наклонам волнового фронта в пределах субапертур в случае работы датчика [8], и в пределах входной апертуры системы при фокусировке оптического пучка [10].

К наиболее применимым на сегодняшний день алгоритмам центроидирования [14] относятся алгоритм взвешенных центров тяжести, дающий наилучший результат вне зависимости от отношения сигнал/шум для точечных источников, квадрантный алгоритм, а также корреляционные алгоритмы. Авторы работы [14] показали, что при низком уровне шума алгоритмы центроидирования могут работать со средним потоком всего десяти фотонов на субапертуру при шуме считывания в три электрона. При высоких уровнях потока фотонов доминирующие ошибки происходят из-за нелинейности отклика, усечений и искажений фокального пятна. Алгоритм взвешенных центров тяжести и корреляционный метод эквивалентны и дают в случае оптимальных параметров лучшие результаты, чем квадрантный алгоритм, эффективность которого зависит от параметров оптической части датчика [8] или системы [13]. Основным недостатком корреляционных алгоритмов является невысокая по сравнению с другими алгоритмами вычислительная скорость.

Показано [14], что варьирование потока фотонов, количества шума считывания и дискретизации пикселей, не являются достаточным условием для разработки рекомендаций по выбору наилучшего алгоритма центроидирования изображения излучения в фокальной плоскости, и наилучшего размера пикселя светочувствительной матрицы видеокамеры.

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

С помощью численной модели датчика волнового фронта Шэка–Гартмана, свободной от шума, обусловленного оптической частью датчика, легко проверить, что именно

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

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

1.    Численные модели

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

В численной модели датчика волнового фронта Шэка–Гартмана прохождение пучка через линзу растра описывается решением параболического уравнения в приближении квазиоптики относительно комплексной амплитуды светового поля, выраженной в пределах субапертуры амплитудным и фазовым распределениями Е(х,у) = Re(E(х,у) + + Ит(Е(х,у)) = А(х,у)ехр(1Ф(х,у')') [3]:

. дЕ

2ik—— = Л^Е + 2k2п(х, у)Е, дг где г - направление распространения; k = 2п/Л - волновое число; Л - длина волны; Л± = д22 + дд2 - оператор Лапласа. Нормирование продольной координаты г на дифракционную длину Z^ = ka0 представляющую пространственный масштаб в направлении распространения излучения, и поперечных координат на начальный радиус входного пучка а0 , заданный в масштабах расчетной сетки в пределах субапертуры, приводит параметры в уравнении (1) к безразмерному виду. Уравнение (1) решается методом расщепления по параметрам с использованием алгоритма быстрого преобразования Фурье (БПФ).

Компонента п (х,у) задает случайное поле флуктуаций показателя преломления, обусловленное турбулентной атмосферой. В инерционном интервале Z ₀ ^ р ^ L ₀ случайное поле флуктуаций показателя преломления с учетом выражения, связывающего структурную характеристику флуктуаций показателя преломления С ² с радиусом Фрида l           3

г ₀ для плоской волны выражением г ₀ = (0,423k²/ C 2 (l)dl) 5 , соответствует случайно- 0

му полю фазовых флуктуаций, заданному спектром фон Кармана:

F _s ⁽ X k,m ) = 0,489г о ^{5 / 3 (} х Р,м + Х 2 ) ^П/6 ^ехР ^{( -} Х р,м /Х^,

где Хк,м = Аха/(К2 + М2) ; К, М - индексы узлов расчетной сетки размерностью N х N; Ах = 2n/(NAp) - шаг дискретизации по пространственной частоте; А"^ -расстояния между разнесенными точками; р = (х, у); Хт = 5,92/10 и Хо = 2n/L0 - спектральные частоты, соответствующие внутреннему 10 и внешнему L0 масштабам турбулентных неоднородностей, соответственно. Фазовые флуктуации светового поля, обусловленные атмосферной турбулентностью, характеризуются структурной функцией фазы [3]:

П ф (? )=6, 88( 1?"| /r o ) 5/3 . (3)

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

Ф(х,у) = [F _s (х к,м )] 1 ^{/ 2} &xexp(t $ rand ), (4)

где ф r_an_d — равномерно распределенная фаза Фурье - образа в интервале [0,2 п ] . Параметр г ₀ в выражении (2) задает поперечный размер неоднородностей, в пределах которого сохраняется когерентность излучения, или устойчивые фазовые соотношения между пространственно разнесенными точками. Другими словами, г ₀ характеризует интенсивность флуктуаций показателя преломления или интенсивность турбулентных искажений фазового распределения Ф(х, у) . В данных исследованиях значения радиуса Фрида варьируются от 1 см (сильная турбулентность) до 10 см (слабая турбулентность).

Распределение интенсивности I(х,у) = Re 2 (E(х,у)) + 1т²(Е(х,у)) в плоскости регистрации светового поля в пределах субапертуры соответствует изображению на светочувствительной матрице видеокамеры, которое содержит фокальное пятно с характерным профилем интенсивности или плотности световой энергии в заданный момент времени.

Координаты центра тяжести распределения интенсивности I(х,у) , представляющего массив пикселей светочувствительной матрицы, каждый из которых характеризуется интенсивностью i _i,j , могут быть записаны в виде [13]:

Na     Na         Na     Na         Na          Na ь = Ex Eii,n. = EJij Eij,ii = Eiij,ij = Eiij,        (5)

i=1     i=1           j=1      j=1          j=1           i=1

где k = 1, 2, ...,N _iens ; N _a xN _a - размерность распределения в пределах к -ой субапертуры линзового растра размерностью N _lens = (N/N _a ) х (N/N _a ) в предположении плотной упаковки линз.

Алгоритм взвешенных центров тяжести (5) не зависит от параметров датчика, представляет координаты центроида в безразмерном виде как результат суммирования весов пикселей i^ , умноженных на расстояние между координатой центральной точки распределения и координатой взвешиваемого пикселя распределения.

Локальные наклоны измеряемого волнового фронта Ф(х,у) пропорциональны смещениям вычисленных координат центроидов k k , п . относительно координат опорного волнового фронта (например, плоского) k . , n . :

^дф( х,у) = k k -^k _в = дф(х,у) = п. ^- П. дх        F iens , . ду        F iens ,

где F _lens - фокусное расстояние линзы.

Искомое фазовое распределение W(х, у) представляется в виде разложения волно вой функции по полиномам Цернике Zl, ортогональным в круге единичного радиуса:

W(х,у) = £ c i Z l (х,у) , где c i - коэффициент, пропорциональный вкладу l - го поли- i=1

нома в волновую функцию; I = 1, 2,..., N b_asis . Реконструкция волнового фронта W (х,у)

осуществляется подгонкой локальных наклонов разложения по полиномам Цернике:

3 Ф(х,у) = _Х' дZ l (x,y)   3 Ф(х,у) = _Х' 3Z i (х,у)

дх ^ дх ^Ch ду ^ 3x ^Cl, 1=1                                   1=1

к оценкам волнового фронта в виде смещений центроидов (6).

Оценивание алгоритма взвешенных центров тяжести выполняется вычислением среднеквадратических отклонений (СКО) моделируемого по формуле (4) распределения Ф(х, у) и распределения W(x,y) , реконструированного по смещениям центроидов, измеренным в плоскости регистрации [9].

СКО разности фазовых распределений Ф(х, у) и W(х,у) , выраженных массивами значений Ф^- и Wi^ , характеризует остаточную ошибку реконструкции волнового фрон-

N                              2

та: Л = J ^2—1    [(Ф^- - Wi-) - (Ф^- - W^>] , где (Ф^- - W^) - среднее значение разности фазовых распределений Ф(х,у) и W(х,у).

Для любого набора случайных неоднородностей показателя преломления (4) зависимость остаточной ошибки реконструкции волнового фронта от интенсивности турбу- - ⁵

лентных искажений фазового распределения пропорциональна величине г ₀ ³ [9]. Другими словами, с увеличением значения радиуса Фрида интенсивность турбулентных искажений фазового распределения уменьшается, соответственно, уменьшается и остаточная ошибка реконструкции.

Для отношения D/r ₀ = 10 , характеризующего разрешение системы (телескопа) с радиусом входной апертуры D = 10 см и параметром Фрида r ₀ = 1 см, остаточная ошибка реконструкции волнового фронта по 35 – ти полиномам Цернике согласно с таблицей Нолла[16] соответствует величине Л = 0, 0135(D/r ₀ ) 5/³ = 0, 6286 радиан.

Рис. 1. Остаточная ошибка реконструкции волнового фронта для случайных фазовых экранов с номером N _scr в зависимости от наличия окна анализа: ряд 1 соответствует алгоритму взвешенных центров тяжести без окна анализа; ряд 2 – с окном анализа.

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

N _a N _a                              N _a N _a

^ = { £ E^ - ^2 + Cm - JГ)1, ]/ E E 4 , } ^{1 / 2} .           (8)

.=1 ,=1                                   .=1 ,=1

Область интегрирования [2 V _a ] x [2 V _a ] в алгоритме (5) меняется на область [ — Ф, Ф] х х [ — Ф, Ф] с центром в точке с координатами i _m ,j _m .

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

Если в плоскости регистрации отсутствует фокальное пятно, то усилив слабоосве-щенные пиксели посредством пороговой фильтрации, можно устранить шумовую составляющую в распределении интенсивности [17].

Далее, умножение структурированного распределения интенсивности и распределения функции Гаусса позволяет выявить фокальное пятно и оценить эффективный радиус в пределах анализируемой субапертуры, который соответствует формуле (8).

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

Введение весовой функции 1 0 (х,у) (кривая 2) позволяет сгладить боковые пики в распределении интенсивности I(х,у) . Положим, что весовая функция или распределение интенсивности I 0 (х,у) , является гауссовым. Результат умножения распределений I(х,у) и 1 ₀ (х,у) соответствует кривой 3. Очевидно, что боковые пики в распределении интенсивности сглажены. Удвоение произведения распределений I(х,у) и 1 ₀ (х,у) не изменяет координаты максимального пика, хотя при этом удваивает значение интенсивности, и делает пик более выраженным (кривая 4).

Если начальное распределение интенсивности является гауссовым, то алгоритм взвешенных центров тяжести, обеспечивает положение середины центроидов для распределений интенсивности, выраженных в пиксельном представлении как ^ и I ., .

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

Рис. 2. Профили, проходящие через середину регистрируемого структурированного распределения интенсивности I ( х,у ) (кривая 1); кривая 2 соответствует гауссову распределению в качестве весовой функции I 0 ( x,y) ; кривая 3 - произведению распределений I ( х,у ) и I 0 ( x,y) ; 4 - удвоенному произведению распределений I ( х,у ) и I 0 ( x,y) .

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

Na     Na      Na      Na             Na Na      Na      Na ak = ко £ г1г/£ 1г - £ Uj/ £ Ij, в. = ко £ jI./ £ I. - £ jty £ 70

i=1 г=1      г=1      г=1              .=1 .=1      j=1

Na             Na             Na где 4 = E kjIjj -I = E 1ц Ij;;Ij = £ I.IIkJ® = £ Ik jk; ко > 1 - коэффициент j=1              i=1              j=1

усиления.

Смена в алгоритме (9) области интегрирования [N_a] х [N_a] на область

[—N _a /6, N _a /6] х [—N _a /6, N _a /6] для распределения двумерной функции Гаусса I ₀ (x,y) = = 1/ у/2 п ехр[—((х — x _m ) 2 + (y — y _m ) 2 )/2] , где x _m и y _m соответствуют координатам центра оцениваемого распределения интенсивности в масштабах расчетной сетки, задает оптимальный размер окна анализа (N _a /3) х (N _a /3) , который составляет практически 99,7 от функции I ₀ (x,y) , представляющей плотность стандартного нормального распределения.

Результаты введения весовой гауссовой функции для оценки локальных наклонов в пределах окна анализа размером, определенным в соответствии с нормальным распределением Гаусса, представлены на рис. 3.

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

3.    Искусственные нейронные сети для центроидирования изображения оптического излучения в плоскости регистрации

Искусственные нейронные сети (ИНС) [4] [5] [11] являются перспективным математическим подходом для решения ряда задач адаптивной оптики. Алгоритм центрои-дирования на основе ИНС как и алгоритм взвешенных центров тяжести (5), его модификация (9), базируемые на измерениях волнового фронта датчиком Шэка–Гартмана, не

Рис. 3. Остаточная ошибка реконструкции волнового фронта в зависимости от набора случайных неоднородностей показателя преломления, где ряд 1 – с помощью алгоритма взвешенных центров тяжести; 2 представляет результат введения в (5) весовой гауссовой функции для оценки локальных наклонов.

нуждается в знании ни параметров атмосферы[2][12], ни параметров оптической части датчика [8]. Общим для всех алгоритмов является необходимость вычисления ошибок измерения, вносимых светочувствительной матрицей камеры в датчике Шэка – Гартмана.

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

Количество нейронов во входном слое будет 4096 и в выходном слое – одна координата. Количество нейронов в скрытом слое уточняется путем эксперимента на основе среднеквадратичной ошибки обучения. Архитектура ИНС для центроидирования изображения оптического излучения в плоскости, изображенная на рис. 4, предполагает ленточное представление распределений интенсивности в пределах субапертуры. Для первой координаты центроида нейроны во входном слое записаны столбцами распределения интенсивности, для второй координаты – строками того же распределения.

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

При определении координаты центроида ^^ первый слой вычисляет слагаемые сум-Na мы £1 = £ Hi , где горизонтальные пиксельные координаты i в каждом столбце рас-t=i пределения интенсивности исполняют роль весов: w^j = i; i,j = 1, 2,...,Na, а интенсив- ности It представляют входные данные, функция активации о = о(£1) равна единице.

алгоритм как гра-

• [2 ]                                                     ^N

Во втором слое w \j = 1 ; функция активации составляет £ ₂ = ^2 I _t .

t=i

Наиболее приемлемым вариантом обучения в данном случае является обратного распространения ошибки (ОРО) Румелхарта, который реализуется диентный метод поиска минимума функции ошибки с рассмотрением сигналов ошибки от выходов ИНС к ее входам, функция ошибки представляет собой сумму квадратов рассогласования (ошибки) желаемого выхода сети и реального.

Рис. 4. Архитектура ИНС для оценивания излучения в фокальной плоскости, где w^ ] - вес нейрона с номером i,j в п - ом слое.

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

Обучение ИНС для центроидирования изображения оптического излучения в плоскости регистрации выполняется по гартманограмме, полученной на этапе калибровки численной модели датчика Шэка–Гартмана, когда фокальные пятна в плоскости регистрации имеют нормальное гауссово распределение. Центры фокальных пятен совпадают с центрами двумерных массивов пикселей, представляющих распределения интенсивности. Для каждого из распределений значения координат калибровочного центроида в результате обучения уменьшаются на три порядка, стремясь к нулю.

Использование ИНС для центроидирования изображения оптического излучения в фокальной плоскости обнаруживает ту же зависимость остаточной ошибки реконструкции волнового фронта от интенсивности турбулентных искажений оптического излучения, что представлено на рис. 5. Применение классических нейронных сетей для цен-троидирования изображения оптического излучения в плоскости регистрации обладает недостатком: большое количество параметров приводит сеть к переобучаемости.

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

Для выборки из 50-ти зарегистрированных гартманограмм (64 субапертуры размером 4096 значений), полученных для искаженного турбулентностью волнового фронта, была проведена проверка эффективности ИНС с архитектурой, представленной на рис. 4.

Рис. 5. Остаточная ошибка реконструкции волнового фронта в зависимости от набора случайных неоднородностей показателя преломления: ряд 1 – с использованием ИНС; 2 – алгоритмом (5) с весовой гауссовой функцией; 3 – посредством алгоритма (5); 4 – с использованием ИНС к произведению регистрируемого распределения интенсивности и весовой функции Гаусса в пределах окна анализа.

Для оценки эффективности определялось число истинно положительных распознаваний значений N _tp и число ложно положительных распознаваний значений N f_p . Рассчитывались метрика, характеризующая долю положительных случаев: N _tp /(N _tp + N f_p ) . Общая точность нейронных сетей для выборки из 50-ти зарегистрированных гартмано-грамм превысила значение 0,87; с учетом окна анализа составила около 100%.

Следует также отметить, что использование ИНС в системах адаптианой оптики не ограничивается алгоритмами центроидирования изображения излучения. Способность нейронной сети к обобщению и обнаружению скрытых зависимостей между входными и выходными данными предполагает ее использование в задачах построения прогноза фазовых искажений излучения, обусловленных турбулентной атмосферой. В качестве объектов для модели прогнозирования [6] используются наборы значений координат центроидов, вычисляемых алгоритмом взвешенных центров тяжести по измерениям датчика Шэка–Гартмана в последовательные моменты времени. Таким образом, формируются временные ряды, которые используются для опережающего формирования фазовой поверхности лазерного излучения в системах адаптивной оптики [1], в том числе и на основе ИНС.

Заключение

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

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

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

Минимизировать фоновую засветку на стенде адаптивной оптической системы рекомендуется оптимальным расположением ее элементов и применением контура опережающей коррекции, что подтверждено полученным Патентом РФ на техническую реализацию следящей адаптивной оптической системы с контуром опережающей коррекции [1].