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

Полный текст

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

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

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

Описание физического эксперимента

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

Оптический элемент помещен в вакуум-камеру, поскольку в воздухе очень сильно действие коэффициента теплообмена. Он имеет большую величину порядка ${10}^{-2} ~ 5\times {10}^{-3}$ и является нестабильным. Поэтому есть необходимость уменьшить влияние коэффициента теплообмена. В вакуум-камере ситуация стабильная, влияние изменения окружающей среды минимально и погрешность меньше. Вакуум-камера поддерживает стабильный вакуум и параметр теплообмена мало меняется (значение  становится близкой к константе). В дальнейшем он предполагается постоянной величиной.

Лазерный луч направляют через окно вакуум-рна оптический элемент. Оптический элемент частично поглощает лазерное излучение, вследствие чего нагревается. В силу теплопроводности происходит распространение тепла от зоны засветки до периферии. На периферии ( $r=R$) внутри оправы установлены датчики температуры, обладающие высокими теплоизолирующими свойствами. С помощью датчиков температуры проводят измерения теплового поля в процессе эксперимента.

В момент включения лазера температура окружающей среды характеризуется величиной $u_{с}$и считается постоянной. Луч входит в оптический элемент под углом к его граням, близким к 90 градусам. Оси пучка и оптического элемента совмещаются.

Оптический элемент представляет собой диск с параллельными гранями кругового сечения радиуса $R$.  из однородного материала толщиной $l$. Также он характеризуется коэффициентом поглощения $\beta$.

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

Математическая постановка прямой и обратной задачи

Математическая модель теплового процесса, протекающего в диске, описывается следующей системой [5]:

 

$\left\{\begin{array}{c}с\frac{\partial u}{\partial t}=k_T\left({\text{Δ}}_{r}u\right)-\frac{2\alpha }{l}u+I\left(r\right)\beta , \\ u{|}_{t=0}=\mathrm{0,} \\ {\left.\frac{\partial u}{\partial r}\right|}_{r=R}=\mathrm{0,} \end{array}\right. \begin{array}{c}0<r\le R, 0\le t\le T;\\ 0\le r\le R; \\ 0\le t\le T; \end{array}$

 

где $u=u\left(r,t\right)$ – разница между температурой оптического элемента и температурой окружающей среды,

$\alpha$ – коэффициент теплообмена,

$k_T$ – коэффициент теплопроводности материала диска,

$l$– толщина оптического элемента,

$c=c_y\rho$ – объемная теплоёмкость,

$I\left(r\right)$ – интенсивность лазерного излучения.

 

Рис. 1. Вид установки

 

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

$I\left(r\right)=\frac{P}{\pi a^2}{e}^{-{\left(\frac{r}{a}\right)}^2},$                                         

где  $a$ – эффективный радиус лазерного пучка,

$P$ – мощность излучения.

Решение задачи запишется в виде:

$u\left(r,t\right)=\underset{n=0}{\overset{+\infty }{}}\frac{\beta P{J}_0\left(\frac{{\mu }_{n}r}{R}\right)e^{-\frac{{\mu }_n^2{a}^2}{4R^2}}}{\pi R^2{J}_0^2\left({\mu }_n\right)\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)k_T}\left(1-e^{\frac{-\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)k_T}{c}t}\right).$

Контроль осуществляется с помощью оценки остатка ряда. Выяснено, что для обеспечения достаточной точности необходимо суммировать 11 элементов ряда [5].

Суть решаемой обратной задачи теплопроводности состоит в вычислении значений параметров уравнения теплопроводности (коэффициента поглощения $\beta$ и коэффициента теплоотвода $\alpha$ с граней оптического элемента во внешнюю среду) по результатам измерений температуры на периферии оптического элемента при $r=R$.

Размещение тепловых датчиков на периферии объясняется необходимостью избежать прямого воздействия на них лазерного излучения. Практический интерес представляет собой оценивание коэффициента поглощения, поскольку готовые образцы, выполненные из таких материалов, как NaCl, KCl, BaF₂ и др. имеют существенный разброс в значениях коэффициента поглощения (десятки процентов).

В качестве критерия оптимальности оценок параметров использовался следующий среднеквадратический критерий качества [6]:

$G\left(\alpha ,\beta \right)=\frac{1}{K}\underset{k=1}{\overset{K}{}}{\left(u\left(R, t_k, \alpha ,\beta \right)-{\vartheta }_k\right)}^2\to \underset{\alpha ,\beta }{\text{min}},$

где $u\left(R, t_k, \alpha ,\beta \right)$ – решение задачи теплопроводности в моменты времени $t_k$ при некоторых значения $\alpha ,\beta$.

${\vartheta }_k$ – результаты измерения температуры периферии элемента в моменты времени $t_k$, $k=\overline{\mathrm{1,}K}$.

Предполагается, что после включения лазера в течение времени $T$ производится регистрация температуры с помощью тепловых датчиков, размещенных на поверхности $r=R$ оптической пластины. При вычислении значения критерия качества $G\left(\alpha ,\beta \right)$ однократно решается прямая задача теплопроводности при значениях параметров $\alpha$ и $\beta$ уравнения, выбираемых в соответствии с алгоритмом оптимизации. 

Для определения оптимальных оценок параметров $\alpha$ и $\beta$ используется итерационный процесс, на каждом шаге которого рассчитывается значение критерия $G\left(\alpha ,\beta \right)$ и, следовательно, предварительно решается краевая задача теплопроводности при выбранных специальным образом значениях $\alpha$ и $\beta$. Для завершения итерационного процесса используется специальное условие остановки.

Исследование функции критерия качества

Для расчета $\alpha$ и $\beta$ уравнения теплопроводности проведем исследование функции цели $G\left(\alpha ,\beta \right)$. Здесь при проведении вычислительных экспериментов по исследованию критерия оптимальности в качестве «измерений» ${\vartheta }_k$ брались результаты расчета теплового поля при значениях параметров $\alpha =\mathrm{0,005}$ и $\beta =\mathrm{0,008}$ с отсутствием шума (рис. 2).

Исследование функции цели $G\left(\alpha ,\beta \right)$ проводилось методом сечений плоскостями $\alpha =const,$ для конкретного $\alpha$ из набора значений вычисляется функция $G\left(\alpha ,\beta \right)$ для различных $\beta$. Данное исследование показало, что функция цели имеет форму криволинейного оврага [7]. Далее была осуществлена минимизация функции $G\left(\alpha ,\beta \right)$.

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

 

Рис. 2. Графики сечений функции цели плоскостями

 

Рис. 3. График линии дна оврага функции

 

Рис. 4. Оптиматизация методом покоординатного спуска

 

Метод покоординатного спуска

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

Для работы метода необходимо задать начальную точку $\left({\alpha }_0,{\beta }_0\right)$. Поскольку в нашей задаче целевая функция зависит от $\alpha$ и $\beta$, тогда зафиксируем одну из переменных. Решая одномерную задачу оптимизации для этой функции, мы от начальной точки $\left({\alpha }_0,{\beta }_0\right)$ перейдем к точке $\left({\alpha }_1,{\beta }_0\right)$, в которой функция  принимает наименьшее значение по координате $\beta$ при фиксированной $\alpha$. В этом состоит первый шаг процесса оптимизации, состоящий в спуске по одной координате.

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

Рассмотрим работу покоординатного метода для исходной задачи из начальной точки $\alpha =\mathrm{0,01}$ и $\beta =\mathrm{0,01}$ (рис. 4).

Метод наискорейшего спуска

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

Шаг 1. Задаются начальное приближение для $\alpha , \beta$ и точность $\epsilon$.

Шаг 2. Находим значение градиента по формулам: $\frac{\partial G}{\partial \alpha }=\frac{\partial G}{\partial u}\times \frac{\partial u}{\partial \alpha }$  и , где $\frac{\partial G}{\partial \beta }=\frac{\partial G}{\partial u}\times \frac{\partial u}{\partial \beta }$

 

$\frac{\partial u}{\partial \beta }=\underset{n=0}{\overset{+\infty }{}}\frac{P{J}_0\left(\frac{{\mu }_{n}r}{R}\right)e^{-\frac{{\mu }_n^2{a}^2}{4R^2}}}{\pi R^2{J}_0^2\left({\mu }_n\right)\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)k_T}\left(1-e^{\frac{-\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)k_T}{c}t}\right);$

$\frac{\partial u}{\partial \alpha }=\underset{n=0}{\overset{+\infty }{}}\left[\frac{2\beta P{J}_0\left(\frac{{\mu }_{n}r}{R}\right)t{e}^{-\frac{{\mu }_n^2{a}^2}{4R^2}}}{\pi R^2{J}_0^2\left({\mu }_n\right)\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)k_{T}cl}{e}^{\frac{-\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)k_T}{c}t}\right.-$

$-\left.\frac{2\beta P{J}_0\left(\frac{{\mu }_{n}r}{R}\right)e^{-\frac{{\mu }_n^2{a}^2}{4R^2}}}{\pi R^2{J}_0^2\left({\mu }_n\right){\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)}^2{k_T}^{2}l}\left(1-e^{\frac{-\left(\frac{{\mu }_n^2}{R^2}+\frac{2\alpha }{k_{T}l}\right)k_T}{c}t}\right)\right].$                   

Шаг 3. Пересчитываем $\alpha , \beta$ по формулам:

${\alpha }^{\left[k+1\right]}={\alpha }^{\left[k\right]}-{\lambda }^{\left[k\right]}\frac{\partial G}{\partial \alpha },$

${\beta }^{\left[k+1\right]}={\beta }^{\left[k\right]}-{\lambda }^{\left[k\right]}\frac{\partial G}{\partial \beta },$ где ${\lambda }^{\left[k\right]}$ увеличивается в $p$ раз, пока значение функции цели уменьшается.

Шаг 4.

Если выполняется $\frac{{\alpha }^{\left[k+1\right]}-{\alpha }^{\left[k\right]}}{{\alpha }^{\left[k+1\right]}+{\alpha }^{*}}\le \epsilon и \frac{{\beta }^{\left[k+1\right]}-{\beta }^{\left[k\right]}}{{\beta }^{\left[k+1\right]}+{\beta }^{*}}\le \epsilon$, где в качестве ${\alpha }^{*} и {\beta }^{*}$ берется малое значение ${10}^{-7}$, то возвращаем текущие значения ${\alpha }^{\left[k+1\right]}, {\beta }^{\left[k+1\right]}$.

Иначе возврат к шагу 3.

Продемонстрируем работу метода наискорейшего спуска (рис. 5).

 

Модифицированный метод золотого сечения

Алгоритм поиска оптимальных значений параметров уравнения теплопроводности для нашей задачи является адаптацией одномерного метода золотого сечения к задаче двумерной оптимизации. Оптимальные значения оцениваемых коэффициентов ищутся как по одной, так и по другой переменной, однако здесь не идёт речь о классической покоординатной оптимизации. Оптимизация осуществлена вдоль линии оврага. Сначала по правилу золотого сечения выбираются точки на оси a без вычисления значений критерия качества. Затем для каждой точки на оси a по одномерному методу золотого сечения осуществляется спуск на дно оврага функции $G\left(\alpha ,\beta \right)$. При этом оптимизация проводится по параметру с использованием опять же одномерного метода золотого сечения. В полученных таким образом точках на дне оврага вычисляются значения критерия $G\left(\alpha ,\beta \right).$ В результате сравнения этих значений производится сужение интервалов неопределённости вдоль линии дна оврага и, следовательно, одновременно происходит сужение интервалов неопределенности как по параметру a, так и по параметру $\beta$ [7]. Продемонстрируем работу метода золотого сечения (рис. 6).

 

Рис. 5. Оптиматизация методом наискорейшего спуска

 

Рис. 6. Оптиматизация модифицированным методом золотого сечения

 

Рис. 7. Оптиматизация методом Гаусса

 

Рис. 8. Оптиматизация методом Гаусса с корректировкой шага, основанной на аппроксимации параболой

 

Рис. 9. Оптиматизация методом Гаусса с корректировкой шага, основанной на методе наискорейшего спуска

 

Метод Гаусса

Метод Гаусса относится к семейству градиентных методов, как и метод наискорейшего спуска. Итерационный метод Гаусса будет иметь следующий вид (рис. 7):

 

Однако неудачно выбранное начальное приближение может значительно увеличить число итераций или даже вызвать зацикливание. Можно заметить, что метод Гаусса стреляет с большим шагом, из-за чего иногда случаются «перелеты» и растет время работы программы. Этого можно избежать с помощью добавления корректировки для шага «выстрела».

Первая корректировка основана на аппроксимации параболой. Её суть заключается в том, чтобы уменьшать или увеличивать шаг в 2 раза и держать в памяти компьютера три оптимальные точки функции $G\left(\overline{x}\right)$. Из графика критерия качества заметно, что функция цели похожа на параболу, экстремум которой будет являться оптимальным значением. Для аппроксимации воспользуемся уравнением $G\left(\overline{x}\right)=a{\overline{x}}^2+b\overline{x}+c$. Подствляя полученные программой точки, решаем систему из трех уравнений и находим коэффициенты $a, b, c$, после чего находим экстремум полученной параболы и сравниваем значение функции цели в точке минимума $G\left(\overline{x}\right),$ с предложенной методом Гаусса, и выбираем наименьшую. Продемонстрируем работу метода Гаусса с применением данной корректировки шага из точки $\alpha =\mathrm{0,01,} b=\mathrm{0,02}$ (рис. 8).

В качестве эксперимента была также реализована вторая корректировка шага, которая применялась ранее в методе наискорейшего спуска. Продемонстрируем работу метода с такой корректировкой шага из точки $\alpha =\mathrm{0,02,} b=\mathrm{0,01}$ (рис. 9).

Обе корректировки не дали большого выигрыша во времени, но позволили методу Гаусса быстрее добраться до линии оврага.

Сравнительный анализ методов

В результат проведения вычислительных экспериментов на «точных» измерениях все методы позволили достичь оптимума с высокой точностью, хотя за разное время. Стоит отметить, что различие между полученными значениями $\alpha$, $\beta$ и ${\alpha }_{true}$, ${\beta }_{true}$  не превосходили 0,0005% при использовании модифицированного метода золотого сечения и 0,0001% при использовании метода Гаусса. Поскольку при проведении натурных экспериментов всегда присутствуют ошибки измерений, важно было проверить работу реализованных методов оптимизации на зашумленных измерениях. Для имитации зашумленных измерений использовалась аддитивная добавка Гауссовского шума с нулевым средним и среднеквадратическим отклонением, равным ${\sigma }^2$. В данной работе было рассмотрено $\sigma$ в диапазоне $\left[\mathrm{0,1};\mathrm{0,3}\right]$, что соответствует параметрам современных динамических датчиков температуры. Формула критерия оптимальности оценок параметров с учетом шума $v$:

$G\left(x,\xi \right)={\left[\xi -u\left(x\right)\right]}^T\left[\xi -u\left(x\right)\right]\to \underset{x}{\text{min}},$           

где $x=\left[\begin{array}{c}\alpha \\ \beta \end{array}\right]$,$\xi = u_{true}\left(x\right)+v$ , $u\left(x\right)=$$\left[\begin{array}{c}u\left(R,t_1,\alpha ,\beta \right)\\ u\left(R,t_2,\alpha ,\beta \right)\\ \dots \\ u\left(R,t_M,\alpha ,\beta \right)\end{array}\right].$

                                                        

В качестве критериев сравнения работы методов использовалось значение погрешности полученных оптимальных значений и среднее время работы программы (табл. 1 и 2). При этом каждый метод неоднократно запускался из разных начальных приближений.

  

Таблица 1

Сравнение средней погрешности разных методов оптимизации на зашумленных измерениях

Метод	Погрешность по $\alpha$  (${\epsilon }_{\alpha }$ ,%)			Погрешность по $\beta$ , ( ${\epsilon }_{\beta }$,%)
	$\sigma =\mathrm{0,1}$	$\sigma =\mathrm{0,2}$	$\sigma =\mathrm{0,3}$	$\sigma =\mathrm{0,1}$	$\sigma =\mathrm{0,2}$	$\sigma =\mathrm{0,3}$
Метод золотого сечения	0,73	2,04	3,25	0,46	0,99	1,87
Метод Гаусса	0,59	1,74	2,91	0,38	0,85	1,55

Таблица 2

Сравнение среднего времени работы разных методов оптимизации на зашумленных измерениях

Используемый метод	$\sigma =\mathrm{0,1}$	$\sigma =\mathrm{0,2}$	$\sigma =\mathrm{0,3}$
Метод покоординатного спуска	450 секунд	500 секунд	600 секунд
Метод наискорейшего спуска	250 секунд	300 секунд	400 секунд
Метод золотого сечения	115 секунд	115 секунд	115 секунд
Метод Гаусса	10 секунд	15 секунд	18 секунд
Метод Гаусса с корректировкой шага	3 секунды	5 секунд	8 секунд

 

На основе полученных результатов можно сделать вывод о том, что с ростом $\sigma$ возрастает погрешность полученных оптимальных значений.  Также стоит отметить, что градиентный метод Гаусса (без корректировки шага) может проявлять свойство зацикливания, не обеспечивая быстрого достижения оптимума, поэтому для решения задачи оптимизации лучше использовать метод Гаусса с корректировкой шага.

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

Можно заметить, что с ростом $\sigma$ возрастает среднее время работы всех градиентных методов. Самым стабильным по времени работы стал модифицированный метод золотого сечения, поскольку он зависит от длины рассматриваемого интервала значений.

Из всех методов самым быстрым и точным методом, работающем на зашумленных измерениях оказался метод Гаусса с корректировкой шага.

Об авторах

Алёна Константиновна Шляхова

Самарский университет

Автор, ответственный за переписку.
Email: gastingz@gmail.com

магистрант факультета информатики

Россия, 443086, Россия, г. Самара, Московское шоссе, 34

Александр Александрович Дегтярев

Самарский университет

Email: aadegt@gmail.com

доцент кафедры технической кибернетики

Россия, 443086, Россия, г. Самара, Московское шоссе, 34

Дополнительные файлы