УМЕНЬШЕНИЕ ПОГРЕШНОСТИ ОПРЕДЕЛЕНИЯ ФАЗОВЫХ РАЗНОСТЕЙ ПРИ АНАЛИЗЕ ИНТЕРФЕРОГРАММ МЕТОДОМ ПОШАГОВОГОВОГО ФАЗОВОГО СДВИГА Текст научной статьи по специальности «Физика»

Уменьшение погрешности определения

фазовых разностей при анализе интерферограмм методом пошагового

фазового сдвига

В.И. Гужов, С.П. Ильиных, И.А. Сажин, Кузнецов Р.А., Хайдуков Д. С.

Новосибирский государственный технический университет (Россия)

Аннотация: В статье рассматривается алгоритм явного определения фазовых сдвигов из серии интерференционных картин, полученных методом пошагового фазового сдвига

Ключевые слова: Метод фазовых шагов, интерферограмма, уравнения расшифровки.

I. ВВЕДЕНИЕ

Наибольшее применение при построении интерференционных систем в последние годы получили методы получения и расшифровки интерферограмм на основе пошагового сдвига (пошаговая или фазо-сдвигающая интерферометрия, phase-sampling, phase-shifting interferometry) [1-12]. Метод пошагового фазового сдвига основан на регистрации нескольких интерференционных картин при изменении фазы опорной волны di на некоторые известные значения.

I, (x, y) = I0 {1 + V(x, y) cos [ j(x, y) + 5, ]} , (1)

где i=0,2, ... , m-1, m - число фазовых сдвигов.

Основное выражение для измерения фазовых разностей можно представить в виде:

m—1

Z k i1 i (x,y)sin6i

f(x,y) = arctanjm—j-

Z pJ i (x,y)cos6i = , (2)

где ki и pi . коэффициенты, определяемые из решения системы трансцендентных уравнений (1).

II. АНАЛИЗ ПОГРЕШНОСТЕЙ

Погрешность при измерении фазы зависит от погрешности при установке фазового сдвига и погрешности при измерении интенсивности [13,14]. Дифференцируя по 8t и получим:

Э j

Ж

Ik kk cos5k Z p.- Ii cos5, +

_1=0____

m —1 L

(Z k;I; (x,y)sin5I )2 +

1=0

m—1

+Ik Pks1n5k Z kiIi (x,y)s1n5I

m—1

Э j

+(Z Pi Ii (x,y)cos5i )2

1=0

m—1

kks1n5k Z р.- i, (x,y)cos5I —

(3)

1=0

ЭI

(Z (x, y) s1n5i )2 +

1=0

m—1

—Pk cos5k Z ki Iis1n5i

1=0

m—1

+(Z Pi I, (x,y)cos5I)

2

1=0

где к = 0,2, ... , т-1.

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

Af » Z

i=

т-1

Af » Z

Эф

эж

Эф

ЭЖ

(±Dj ) + Z

i=0

т-1

( ±Э j)+Z

Эф

j

Э h Эф

где Э^ =

ЭI, =

Э h

DI

(±A It) , (4) I,

j

I

Эj=j. Z

j i=0

j ff ЭЖ, j

m—1

+2:

1=0

j II

Э I, j

(±Э I, )

( ±Э IJ,

(5)

(±Эф) +

(6)

=0

i=0

a

В выражении для погрешности (4) коэффициенты

определения абсолютной

Dj Dli

9 11

имеют различный порядок. Из выражений (2) и (3) видно, что эти коэффициенты отличаются множителем 1к. Для оценки влияния погрешностей, вносимых неправильным определением интенсивности и погрешностей при установке фазового сдвига, перепишем (4) в виде:

III. ОПИСАНИЕ МЕТОДА И ПРАКТИЧЕСКАЯ РЕАЛИЗАЦИЯ

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

m-1

I

Df

i=0

Таким

дф

dd

(±Dj ) +1

дф д I

(±д Ii) , (7)

образом, оценка абсолютной погрешности фазы зависит от относительной

АI;

ошибки при измерении интенсивности —у— и

случайной флуктуации установки фазового

сдвига А 9;. Точность измерения фазы зависит от

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

Рис. 1 Две произвольные точки на интерферограмме, по которым определяется фазовый сдвиг

На Рис. 2 показаны профили интерференционных полос, получаемые при внесении фазовых сдвигов для двух точек интерферограмм А и В, соответственно.

Рис. 2(A). Профиль интерференционных полос в точке 1

Рис. 2(B). Профиль интерференционных полос в точке 2

i=0

необходимо определить коэффициенты

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

I,,1 = 10,1 [1 + +3,)]

I,,2 = 10,2 [1 + +3 )] , (8)

где ,=0,2, ... , т-1, т >= 5.

При изменении углов сдвига от 0 до 2п точка на комплексной плоскости описывает некоторую траекторию. Несложно показать, что данная траектория представляет собой центральную кривую второго порядка (произвольно ориентированный эллипс) с центром в точке (101,102).

100 150 200

Интенсивность в точке 1

Рис.3 Идеальная траектория, которую пробегают значения интенсивностей в двух точках

Рис.4 Определение параметров эллипса по 5 точкам. (й и й2 - главные оси эллипса, 101,1о,2 — координаты центра).

Для определения характеристик траектории

уравнения аппроксимирующей кривой а11х2 + а22у2 + 2а12ху+

+2 а13х + 2 а23у + а33 = 0 , (9)

где х и у - значения интенсивности в выбранных точках. Для определения параметров эллипса достаточно 5 точек (Рис.

4).

Коэффициенты а^ уравнения (9) можно получить путем вычисления определителя минора при соответствующем коэффициенте [15] (стр.64-76).

a11 a1 a a13 a 3

In I11 I 1 I 1 I11 I 1

In I1 I I 1 22 I1 I

I23 I13 I 3 I 3 I13 I 3

114 I14 I 4 I 4 I14 I 4

I is IisIis I". Iis I 25

a

33

(10)

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

Рис. 4 Траектория в комплексной плоскости реальных интерференционных сигналов

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

аппроксимации полиномом второго порядка.

Учитывая избыточность исходных данных, можно использовать для этой цели метод наименьших квадратов.

" N 2 х 2 у, 2 х,у, 2 х,2 2 у2" -1

а11 а12 а22 2 х, 2*2 2 х,у, 2х2 у, 2х3 2™2

2 у, 2 х,у, 2 у,2 ^х,2 у, 2 у,3

а13 2 х,у, 2х2 у, 2х,у2 2х,2 у,2 ^х,3 у, 2х,у3 е

а14 2х; 2 х3 2х2 у, 2х3 у, 2х4 ^х,2 у,2

2У? 2Х,У2 2у,3 2^? 22 2 хгуг 2 у,4

где е - единичный вектор с размерностью равной количеству коэффициентов уравнения (8).

Однако в результате применения данного алгоритма возникают существенные

вычислительные погрешности. Поэтому мы советуем применять следующий алгоритм для нахождения коэффициентов уравнения (9).

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

1) Формируем матрицу Б, состоящую из следующих элементов:

В =

х,

х.

х„

У х1 У

У, ХгУг

У х у

У п пУ п

х

У1 1

У,

х

Уп

1

(13)

где х, и у, — координаты ¡-ой точки на плоскости, т.е. значения интенсивностей в точке 1 и точке 2, ¡=1 ... п, п - общее количество точек.

2) Вычисляем матрицу S

Б=БТ В . (14)

3) Решается обобщенная задача на собственные значения

Бд=ХСд , (2-7) где матрица С имеет следующий вид:

С =

0 0 2 0 0 0

0 -1 0 0 00

2 0 0 0 00

0 0 0 0 00

0 0 0 0 00

0 0 0 0 0 0

шаге необходимо

(15)

На данном шаге необходимо найти собственные значения X и соответствующие им собственные векторы д, удовлетворяющие условию (14).

4) Из найденных собственных значений необходимо выбрать минимальное.

Собственный вектор, соответствующий минимальному собственному значению и будет

определять коэффициенты уравнения эллипса. Пусть д=[д1 д2 д3 д4 д5 дб ]Т — собственный вектор, соответствующий минимальному собственному значению. Тогда коэффициенты уравнения эллипса могут быть определены следующим образом:

С

ап = д1

а22= д2

азз= дб реализацией

а 12= дз /2 а13= д4 /2 а 23= д5 /2

первых двух пунктов

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

Вещественное число X и вектор д называются собственной парой матрицы А, если они удовлетворяют следующему условию:

Ад = Хд . (16)

В случае если вещественная матрица А размером N х N симметрична, у неё есть N собственных чисел (не обязательно различных) и N соответствующих им собственных векторов, образующих ортонормированный собственный базис.

Для решения уравнения (16) имеется много различных алгоритмов. Первым алгоритмом, решающим задачу собственных значений для симметричной матрицы размером N х N, был алгоритм Якоби, приводящий матрицу к диагональной форме ортогональными преобразованиями. По мере осуществления преобразований исходной матрицы, элементы за пределами главной диагонали уменьшались, а на главной диагонали - увеличивались. Естественным результатом этого процесса является матрица, у которой внедиагональные элементы равны нулю, а на диагонали находятся собственные значения. Есть и другие алгоритмы, например, семейство алгоритмов, основанных на рЬ/рЯ-итерации, применяемой к трехдиагональной матрице. Наиболее подробно проблема собственных значений рассматривается в [16].

Обобщенная задача на собственные значения, в отличие от задачи (16), состоит в нахождении такого числа X и вектора д, чтобы они удовлетворяли следующему условию:

Ад=ХВд , (17)

1

х

где B — матрица, относительно которой ищутся собственные значения матрицы A.

Для решения задачи (17) также имеются различные методы. Большая часть из этих методов требует, чтобы матрица, относительно которой ищется собственное значение, была положительно-определенной. В нашем случае мы ищем собственные значения матрицы S относительно матрицы C, а матрица C не является положительно-определенной. И еще, обычно используемые на практике алгоритмы поиска собственных значений носят итерационный характер. Поэтому при программной реализации данные алгоритмы будут давать приближенные собственные значения, что в нашем случае является нежелательным, так как это приведет к уменьшению точности рассматриваемого нами алгоритма расшифровки.

Поэтому задачу (14) будем решать следующим образом. Для нахождения собственных значений матрицы S относительно матрицы C необходимо найти корни следующего уравнения

| S - ÁC | = 0 . (18) Раскрываем определитель левой части уравнения (18), содержащий 720 членов (матрица размерности 6*6, количество элементов 6! = 720). В результате получаем уравнение третьей степени следующего вида

x3+ax2+bx+c=0. (19)

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

Q = (a2 - 3b)/9

R = (2a3 - 9ab + 27c)/54 Далее, если R2< Q3, то уравнение имеет три действительных корня, которые вычисляются по формулам

R

t = arceos—¡=

3VQ3

x = 2^Q ■ cos(t) - a /3 (20) x2 = 2^¡Q ■ cos(t + 2p) - a/3 (21)

x3 = 2JQ ■ cos(t - 2p) - a/3 (22)

В том случае, когда R2 > Q3, то действительных корней один (общий случай) или два (вырожденные случаи). Кроме действительного корня, имеется два комплексно-сопряженных. Для их нахождения вычисляются

A = sign( R)3¡\ R + B = Q / A при A Ф 0 или B = 0 при A = 0

Действительный корень:

xx = (A + B) - a/3

Комплексно-сопряженные корни:

A + B „../r A - B ■-a/3 ±/V 3-

x

2,3

(24)

2 2 В том случае, когда А=В , то комплексно-сопряженные корни вырождаются в действительный:

х2 =-А - а/3 (25) В нашем случае уравнение (17) будет иметь три действительных корня. Поэтому для нахождения корней достаточно

воспользоваться формулами (20) - (25).

При программной реализации алгоритма нахождения корней уравнения (17) возникли проблемы следующего характера. При вычислении коэффициентов а, Ь и свободного члена с кубического уравнения (17) возникла необходимость в более точных методах вычислений. Стандартные, встроенные в языки программирования типы данных с двойной точностью не справлялись: возникало переполнение в силу одновременного оперирования слишком большими и слишком малыми величинами. Для разрешения этой проблемы был реализован калькулятор неограниченной точности на основе рациональных чисел, при помощи которого были вычислены необходимые коэффициенты кубического уравнения.

На Рис. 5 представлен результат аппроксимации траектории интерференционных сигналов по предложенному выше алгоритму.

После определения коэффициентов уравнения эллипса необходимо рассчитать координаты его центра по следующим формулам:

x0 = I,

01

a13 a12

a14 a22

a11 a12

a12 a22

y0 = I

02

a11 a13

a12 a14

a11 a12

a12 a22

(26)

Углы фазового сдвига можно определить непосредственно, если преобразовать траекторию к виду круговой траектории.

100 150 200

Интенсивность в точке 1

Рис. 5. Сплошная линия - результат аппроксимации кривой второго порядка

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

а) привести центр эллипса к началу координат

х1 = х - х0, у1 = у - у0 ; (27)

б) развернуть эллипс параллельно одной из координатных осей.

(28)

На рис. 6 показаны исходная и скорректированная траектории интерференционных сигналов.

" x2" cos W sin W " x1

_ У 2_ - sin W cos W_ _ y1

Угол поворота определяется по формуле:

, <-л 2a tgW =

12

a,

a

(29).

*11 а22

в) растяжение эллипса до круга. Коэффициент растяжения у определяется из канонического уравнения эллипса, выраженного через его инварианты

10 х2 +1у2 +—

L

0,

(30)

где l уравнения

и

l

a11 -l a0

корни характеристического

Рис. 6 Исходная и скорректированная траектории интерференционных сигналов

Фазовые сдвиги определяются по координатам (круговой) траектории

a

a0

12 -l

0

с У 2 ¡

d = arctg —.

x2,.

(32)

21 22

из уравнения (2-5) видно, что отношение корней характеристического уравнения равно отношению диаметров эллипса

у.

l

Растяжение

эллипса

по

(31) координате

У 2

На следующих рисунках приведены результаты расчетов по модельным интерференционным картинам. Число фазовых сдвигов т = 32.

На Рис. 7 показана интерференционная картина с 256 уровнями интенсивности с двумя выбранными точками и траектория, построенная по 32 фазовым сдвигам. Фазовые сдвиги произвольные.

y выполняется следующим образом y2 =-.

g

Рис. 7. Интерференционная картина с двумя выбранными точками и траектория, построенная по 32 фазовым сдвигам

На Рис. 8 графики заданных и полученных фазовых сдвигов. Среднеквадратичное отклонение составило 0,00412 рад. Это вызвано дискретностью задания интенсивности.

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

0 г"|—I—I—I—|—I—I—I—I—|—I—I—I—I—|—I—I—I—I—|—I—I—I—I—|—I—I—I—I—|—I—I—I—ь

0 5 10 15 20 25 30 35

Отсчеты

Рис. 8. Графики заданных и рассчитанных фазовых сдвигов (32 сдвига)

IV. ВЫВОДЫ

Реализован эффективный метод

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

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

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

Достигнута точность расшифровки не менее 1% при исходной 20% погрешности установки фазовых сдвигов.

ЗАКЛЮЧЕНИЕ

Работа выполнена по заданию Министерства образования и науки по проекту «Исследование предельных точностей оптических методов измерения параметров движения и мехатронных методов управления движением и разработка новых робототехнических и электромеханических систем», Темплан, проект № 7.559.2011.

ЛИТЕРАТУРА

[1] Hariharan P., Oreb B.F., Brown N. Digital phase-measurement system for real-time holographic interferometry // Optics Communication.- Vol.41.- №6.-1982.- pp.393-398

[2] Wyant J.C., Creath K. Recent advances in interferometric optical testing // Laser Focus.- 1985. - pp. 118132.

[3] Wyant J.C. Interferometric optical metrology: basic system and principles // Laser Focus.- 1982.- pp.65-67.

[4] Creath K. Phase-shifting speckle interferometry // Applied Optics. 1985. V.24. P.3053-3058.

[5] J.E.Greivenkamp and J.H.Bruning, "Phase shifting interferometry," in Optical Shop Testing, Ed. byD.Malacara (Wiley, New York, 1992), Chapter 14, pp. 501-598.

[6] P. de Groot. Phase-shift calibration errors in interferometers with spherical Fizeau cavities // Applied Optics.-1994.-V.34.-No.16.-pp.2856-2863.

[7] P. de Groot. 101-frame algorithm for phase shifting interferometry. EUROPTO, 1997, Preprint 3098-33.

[8] J. Millerd, N. Brock, J. Hayes, et al., "Modern Approaches in Phase Measuring Metrology," Proc. SPIE. 5856, 14-22 (2004).

[9] P. Gao, B. Yao, N. Lindlein, et al., "Phase-Shift Extraction for Generalized Phase-Shifting Interferometry," Opt.Lett., 2009, 34 (22), 3553-3555.

[10] Гужов В.И., Ильиных С.П., Хайдуков Д.С., Вагизов А.Р. / Универсальный алгоритм расшифровки. // Научный вестник НГТУ. - 2010. - №4(41) - С. 51-58.

[11] Гужов В.И., Ильиных С.П., Хайдуков Д.С., Вагизов А.Р. / Устранение ошибок фазового сдвига в интерферометрии // Автометрия. - 2011. - Т. 47, №1.-С. 96101.

[12] V.I. Guzhov, S. P. Il'yinykh, D. S. Khaidukov and A. R. Vagizov Eliminating phase-shift errors in interferometry // Optoelectronics, Instrumentation and Data Processing.-2011., Vol.47, Nu.1.- pp. 76-80

[13] Гужов В.И., Солодкин Ю.Н. Анализ точности определения полной разности фаз в целочисленных интерферометрах // Автометрия.-1992.-№6.-С.24-30.

[14] Schmit J., Creath K. Extended averaging technique for derivation of error-compensating algorithms in phase-shifting interferometry. //Applied Optics.-1995.-V.34.-No.19.-pp.3610-3619.

[15] Г. Корн, Т. Корн Справочник по математике.-М., 1973, 832 с.

[16] Дж. Х. Уилкинсон. Алгебраическая проблема собственных значений.-М. Наука, Физматлит.-1970, 564 с.

Гужов Владимир

Иванович - декан факультета Автоматики и вычислительной техники Новосибирского Государственного Технического университета, профессор, доктор

технических наук. Он является автором 120 научных работ, в том числе является обладателем 4 патентов.

Область научных интересов: программные системы, высокоточные измерения.

E-mail: vig@nstu.edu.ru

Ильиных Петрович

кафедры тельная

Сергей

- доцент «Вычисли-техника» в

НГТУ, доцент, кандидат технических наук, автор более 100 научных статей, в том числе, 4 патента и 1 учебник высшей школы. Область научных Интересов и компетенций -лазерные измерительные системы, обработка изображений. E-mail: isp51@yandex.ru

Сажин Игорь

Александрович - доцент кафедры «Кафедра технической

теплофизики» в НГТУ, доцент, кандидат

технических наук.

Хайдуков Дмитрий Сергеевич - аспирант Новосибирского Государственного Технического университета.

Кузнецов Роман

Александрович -

аспирант Новоси-

бирского Государственного Технического университета.