1. Введение

Гильберт-диагностика фазовых структур является результатом плодотворного синтеза методов, развиваемых в оптике и радиотехнике [1, 2]. Это интегральная операция, выполняющая перераспределение энергии сигнала в заданной полосе пространственных частот зондирующего поля, возмущённого исследуемой средой. При этом минимизируются энергетические потери оптического сигнала. Преобразование Гильберта в частотном пространстве сводится к определённому виду фазовой трансформации фурье-спектра сигнала с сохранением энергии в широкой полосе частот.

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

В настоящей работе для оптимизации итерационного алгоритма реконструкции фазовой функции при гильберт-диагностике предлагается применение метода Гаусса-Ньютона [5, 6]. Это – итерационный численный способом решения задачи наименьших квадратов, модификация метода Ньютона для поиска минимума целевой функции, в отличие от которого не требуется определения вторых производных, что упрощает и сокращает количество вычислений [7]. В [8–10] представлены различные модификации метода Гаусса-Ньютона, увеличивающие скорость сходимости и снижающие влияние некорректности в постановке задачи, в частности алгоритм Левенберга-Марквардта [11]. Метод Гаусса-Ньютона эффективен в решении многих задач оптимизации, он прост в реализации и присутствует в большинстве программных пакетов по прикладной математике.

2. Обратная задача реконструкции фазовой функции по гильбертограммам

Обратимся к рис. 1, где упрощённо показано выделенное сечение исследуемого фазового объекта с радиальным распределением показателя преломления , ось z совпадает с направлением зондирующего светового пучка.

 

Рис. 1. Схема диагностики фазового объекта (осесимметричное приближение): – радиальное распределение показателя преломления.

 

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

где – длина волны, – показатель преломления воздуха, и – точки входа и выходы луча из объекта.

Математической моделью гильберт-визуализации фазового сдвига является нелинейный интегральный оператор I рода:

(1)

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

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

В случае гладких полей искомая функция фазового возмущения в сечении в интервале моделируется полиномом Бернштейна n-го порядка (кривой Безье) [12]:

(2)

где – компоненты векторов опорных вершин, – базисные функции кривой Безье, называемые полиномами Бернштейна:

– степень полинома, – порядковый номер опорной вершины. Вне интервала функция = 0.

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

(3)

Поскольку , где , уравнение (3) можно представить, как

(4)

Задача оптимизации заключается в определении тех значений параметров (, …, ) и (, …, ), при которых достигается минимум целевой функции:

(5)

где – гильбертограмма, регистрируемая в эксперименте (эталонные данные).

3. Метод Гаусса-Ньютона

В общем виде задача оптимизации состоит в нахождении экстремума (минимума) целевой функции:

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

Матрица Якоби определяется по формуле [13, 14]:

Метод Гаусса-Ньютона заключается в выполнении последовательных приближений :

(6)

где – номер итерации, – коэффициент, служащий для регулирования «шага» оптимизации [6], – транспонированная матрица.

Обозначим

, .

Тогда в случае обработки гильбертограмм компоненты целевой функции (5) можно записать в виде

(7)

Для определения матрицы Якоби найдём значения частных производных по компонентам вектора от функции (4). Вначале вычислим значения производных по параметрам :

Тогда

Поскольку

и

получаем:

где

Таким образом,

или

где

	(8)

Преобразуем выражение (8) интегрированием по частям:

Поскольку

получаем:

В результате приходим к равенству:

Аналогично можно показать, что

Теперь найдём значения производных по параметрам :

Тогда

Поскольку

получаем:

или

Аналогично

В результате матрица Якоби будет иметь следующий вид:

(9)

С учётом (7) и (9) алгоритм Гаусса-Ньютона (6) будет определяться в следующем виде:

(10)

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

4. Применение алгоритма на примере тестовой функции

В случае осесимметричных объектов фазовую функцию достаточно определить только на интервале .

Смоделируем пример, когда диагностируется осесимметричный объект радиусом сечения = 30 мм. Зададим в пределах 0 ≤ x ≤ параметрически кривой Безье третей степени:

при значение 0, и продолжим её чётным образом в область
≤ x ≤ 0.

В проводимых экспериментах [3, 4] при регистрации оптического сигнала размер пикселя в изображении соответствовал значению 1/30 мм. Поэтому в интервале ] количество отсчетов = 1801.

На рис. 2 представлен график функции и её гильбертограммы при значениях опорных точек:

{ = 0, = 50},   { = 15, = 50},  { = 15, = 0}, { = 30, = 0}.

 

Рис. 2. Фазовая функция, представленная кривой Безье (с указанием опорных вершин), и рассчитанная от неё гильбертограмма. Гильбертограмма представлена в безразмерном виде с указанием точек локальных экстремумов.

 

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

 

Рис. 3. Эталонная и начальная (оптимизируемая) гильбертограммы и соответствующие им фазовые функции.

 

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

 

Таблица 1. Производные по и координатам

Координаты	Координаты

 

Матрица Якоби будет иметь следующий вид:

В результате применения метода Гаусса-Ньютона (10) через 30 итераций среднеквадратическая погрешность между и достигла величины:

при этом максимальное отклонение составило:

Полученная в результате оптимизации фазовая функция представлена на рис. 4.

 

Рис. 4. Результат оптимизации.

 

Максимальное отклонение от составило:

В рассмотренном примере начальное приближение было выбрано таким образом, чтобы количество гильберт-полос и совпадало. Попробуем сделать оценку «предельного» отклонения начального приближения , при котором алгоритм Гаусса-Ньютона начинает работать некорректно (рис. 5 и 6).

 

a

б

Рис. 5. (а) – Эталонная и начальная (оптимизируемая) гильбертограммы и соответствующие им фазовые функции (второе приближение); (б) – результат оптимизации.

 

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

Теперь зададим так, чтобы количество гильберт-полос было больше на одну в тех же переделах и (рис. 6. а).

 

a

б

Рис. 6. (а) – Эталонная и начальная (оптимизируемая) гильбертограммы и соответствующие им фазовые функции (третье приближение); (б) – результат оптимизации.

 

В результате через 41 итерацию было получено (рис. 6. б):

Таким образом, было установлено, что для сходимости алгоритма Гаусса-Ньютона необходимо, чтобы число гильберт-полос и было равным. На практике, конечно же, нужно стараться, чтобы гильберт-полосы экспериментальной и оптимизируемой гильбертограмм были насколько возможно «приближены» друг к другу.

5. Заключение

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

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

 

Работа первого автора выполнена в рамках государственного задания ИМ СО РАН № FWNF-2022-0009, остальных авторов – в рамках государственного задания ИТ СО РАН № 121031800217-8.

Список литературы

1. Дубнищев Ю. Н., Арбузов В. А. Методы гильберт-оптики в измерительных технологиях. Новосибирск: Изд-во НГТУ, 2007, 316 с.

2. Сороко Л. М. Гильберт-оптика. М.: Наука, 1981, 160 с.

3. Arbuzov V. A., Arbuzov E. V., Dubnishchev Yu. N., Lukashov V. V., Zolotukhina O. S. Method of polychromatic Hilbert diagnostics of phase and temperature perturbations of axisymmetric flames // CEUR Workshop Proceedings, Vol. 3027, 2021, pp. 369–378 (doi: 10.20948/graphicon-2021-3027-369-378).

4. Arbuzov V. A., Arbuzov E. V., Dubnishchev Yu. N., Lukashov V. V., Zolotukhina O. S., Tupikin A. V. Hilbert-optic diagnostics of hydrogen-oxygen inverse diffusion flame // Energies, Vol. 15 (24), 2022, P. 9566 (doi: 10.3390/en15249566).

5. Floudas C. A., Pardalos P. M. Encyclopedia of optimization. New York: Springer, 2008, 4622 p.

6. Алгоритм Левенберга-Марквардта для нелинейного метода наименьших квадратов и его реализация на Python, 2016. URL: http://habr.com/ru/articles/308626/.

7. Ловецкий К. П., Севастьянов Л. А. и др. Математический синтез оптических наноструктур. М.: РУДН, 2008, 145 с.

8. Юдин Н. Е. Модифицированный метод Гаусса-Ньютона для решения гладкой системы нелинейных уравнений // Компьютерные исследования и моделирование, Т. 13 (4), 2021, С. 697–723.

9. Loke M. H., Dahlin T. A comparison of the Gauss-Newton and quasi-Newton methods in resistivity imaging inversion // Journal of Applied Geophysics, Vol. 49 (3), 2002, pp. 149–162 (doi: 10.1016/S0926-9851(01)00106-9).

10. Bergou E. H., Diouane Y., Kungurtsev V. Convergence and complexity analysis of a Levenberg-Marquardt algorithm for inverse problems // Journal of Optimization Theory and Applications, Vol. 185, 2020, pp. 927–944 (doi: 10.1007/s10957-020-01666-1).

11. Васин В. В. Метод Левенберга-Марквардта для аппроксимации решений нерегулярных операторных уравнений // Автоматика и телемеханика, № 3, 2012, С. 28–38.

12. Роджерс Д., Адамс Дж. Математические основы машинной графики. М.: Мир, 2001, 604 с.

13. Зорич В. А. Математический анализ. Часть I. Изд. 6-е, доп. М.: МЦНМО, 2012, 702 с.

14. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Бином. Лаборатория знаний, 2003, 632 с.

---

Gauss-Newton Method in the Problem of Optimizing the Axisymmetric Phase Function Calculation Based on the Hilbert Diagnostic Data

Authors: E.V.  Arbuzov^1,А,В, V.A.  Arbuzov^2,A, Yu.N.  Dubnishchev^3,A, O.S.  Zolotukhina^4,A

^A Kutateladze Institute of Thermophysics of the Siberian Branch of the Russian Academy of Sciences

^B Sobolev Institute of Mathematics of the Siberian Branch of the Russian Academy of Sciences

¹ ORCID: 0000-0001-9488-8650, arbuzov@math.nsc.ru

² ORCID: 0000-0003-2404-326X , arbuzov@itp.nsc.ru

³ ORCID: 0000-0001-7874-039X, dubnistchev@itp.nsc.ru

⁴ ORCID: 0000-0003-3486-4459 , melexina-olga17@yandex.ru

 

Abstract

A method for reconstructing phase disturbances of a probing light field using the iterative Gauss-Newton algorithm is discussed as part of the Hilbert diagnostics development of gaseous, condensed and reacting media. In this case, the need to determine second derivatives is eliminated, which simplifies the calculations. The method consists of selecting a phase profile, which is specified by a Bezier curve, and hilbertogram calculating. The coincidence of the reference and reconstructed hilbertograms serves as a criterion for the results reliability. The Jacobian matrix for the nonlinear integral operator of Hilbert visualization is obtained. The algorithm is analyzed using a test function. The method development is associated with the algorithm application to the processing of experimental results, including the reconstruction of complex structures in which the phase function is described by several Bezier polynomials.

 

Keywords: Hilbert optics, phase function, optimization, Gauss-Newton method.

 

References

1. Arbuzov V. A., Dubnishchev Yu. N. Metody gil'bert-optiki v izmeritel'nykh tekhnologiyakh [Hilbert-optics methods in measurement technologies], NSTU University Publ., Novosibirsk, 2007, 316 p. (in Russian)

2. Soroko L. M. Gil'bert-optika [Hilbert optics], Nauka Publ., Moscow, 1981, 160 p. (in Russian)

3. Arbuzov V. A., Arbuzov E. V., Dubnishchev Yu. N., Lukashov V. V., Zolotukhina O. S. Method of polychromatic Hilbert diagnostics of phase and temperature perturbations of axisymmetric flames // CEUR Workshop Proceedings, Vol. 3027, 2021, pp. 369–378 (doi: 10.20948/graphicon-2021-3027-369-378).

4. Arbuzov V. A., Arbuzov E. V., Dubnishchev Yu. N., Lukashov V. V., Zolotukhina O. S., Tupikin A. V. Hilbert-optic diagnostics of hydrogen-oxygen inverse diffusion flame // Energies, Vol. 15 (24), 2022, P. 9566 (doi: 10.3390/en15249566).

5. Floudas C. A., Pardalos P. M. Encyclopedia of optimization. New York: Springer, 2008, 4622 p.

6. Levenberg-Marquardt algorithm for non-linear least squares method and its implementation in Python, 2016. URL: http://habr.com/ru/articles/308626/.

7. Lovetsky K. P., Sevastyanov L. A., Bikeev O. N., Paukshto M. V. Matematicheskiy sintez opticheskikh nanostruktur [Mathematical synthesis of optical nanostructures], RUDN University Publ., Moscow, 145 p. (in Russian)

8. Vasin V. V. The Levenberg-Marquardt method for approximation of solutions of irregular operator equations // Automation and remote control, Vol. 73, 2012, pp. 440–449 (doi: 10.1134/S0005117912030034).

9. Yudin N. E. Modified Gauss–Newton method for solving a smooth system of nonlinear equations // Computer research and modeling, Vol. 13, pp. 2021, pp. 697–723 (doi: 10.20537/2076-7633-2021-13-4-697-723).

10. Loke M. H., Dahlin T. A comparison of the Gauss-Newton and quasi-Newton methods in resistivity imaging inversion // Journal of Applied Geophysics, Vol. 49 (3), 2002, pp. 149–162 (doi: 10.1016/S0926-9851(01)00106-9).

11. Bergou E. H., Diouane Y., Kungurtsev V. Convergence and complexity analysis of a Levenberg-Marquardt algorithm for inverse problems // Journal of Optimization Theory and Applications, Vol. 185, 2020, pp. 927–944 (doi: 10.1007/s10957-020-01666-1).

12. Rogers D. F., Adams J. A. Matematical elements for computer graphics (Second Edition), McGraw-Hill, New York, 1990.

13. Zorich V. A. Matematicheskiy analiz [Mathematical analysis], part 1, 6st. ed., ICNMO Publ., Moscow, 2012, 702 p. (in Russian)

14. Bakhvalov N. S., Zhidkov N. P., Kobelkov G. M. Chislennyye metody [Numerical methods], Binom Publ., Moscow, 2003, 632 p. (in Russian)