Модифицированный алгоритм расчета поля векторов смещений для оценки деформации Текст научной статьи по специальности «Математика»

УДК 004.932

МОДИФИЦИРОВАННЫЙ АЛГОРИТМ РАСЧЕТА ПОЛЯ ВЕКТОРОВ СМЕЩЕНИЙ

ДЛЯ ОЦЕНКИ ДЕФОРМАЦИИ

А.И. Солодушкин, В.В. Кибиткин, В.С. Плешанов

Институт физики прочности и материаловедения СО РАН, г. Томск E-mail: vvk@ispms.tsc.ru

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

Ключевые слова:

Векторное поле, метод корреляции, цифровые изображения, алгоритм шаблонного поиска, передискретизация. Key words:

Vector field, correlation technique, digital images, template search algorithm, resampling.

Введение

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

Метод корреляции цифровых изображений (digital image correlation) является одним из перспективных подходов для изучения процессов деформации и разрушения материалов. Он основан на анализе цифровых изображений поверхности. Данный метод впервые был предложен группой исследователей университета в Южной Каролине в 1982 г. [2].

На базе этого метода в 90-х гг. прошлого века в Институте физики прочности и материаловедения СО РАН, г. Томск, была разработана оптико-телевизионная измерительная система (ОТИС) «TOMSC». Методика измерений заключается в следующем. Изображение изучаемого участка поверхности объекта воспринимается оптическим микроскопом и с помощью цифровой камеры записывается на жесткий диск компьютера. Данное изображение будем называть эталонным. Далее образец подвергается воздействию механических нагрузок и снова с этого же участка поверхности фиксируется второе - текущее изображение.

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

S = Х f (Ii, I2),

Q

где I1 и I2 - массивы значений яркости изображений; /(/1,/2) - мера близости; Q - область сканирования.

Затем по известным значениям проекций ux, uy векторного поля ux(x,y)=ux(x,y)e+u'(x,y)ey (где ex, еу - единичные орты) рассчитываются локальные

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

Базовый алгоритм идентификации вектора смещения

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

Для расчета поля векторов смещений в ОТИС используется блочный алгоритм. Работа алгоритма осуществляется следующим образом:

1. Первое изображение разбивается на некоторое число эталонных участков (ЭУ) размером мхм, второе изображение делится на области поиска (ОП) размером ЯхЯ (К>ш).

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

3. По координатам глобального экстремума функционала находится искомый вектор (рис. 1).

Рис. 1. Схема, поясняющая работу блочного алгоритма идентификации

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

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

* =И|'и,У -1

2, ¡, IV

/=1 1-1

-ТТ('и, 1 '2,1,У)2> -

$1

1-1 1-1

где

$1 -ХХ(7и 1 - ки 1)(12,1,1 -12,1,1), 1-1 1-1

$2 - ТТ('и 1 -!и 1 )2' $3 - ТТ(к 1, 1 -12, 1, 1 )2, 1 -1 1 -1 1-1 1

/; и /2 - средние арифметические значения яркости изображений в пределах области размером тхт; г,] - текущие индексы по строкам и столбцам.

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

Оптимизация поиска по времени

эталонного участка изображения

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

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

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

кая вероятность найти локальный экстремум вместо глобального.

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

• 1 этап: передискретизация (понижение разрешения изображений) эталонного участка и соответствующей области поиска;

• 2 этап: нахождение приближенных значений координат глобального экстремума (харр,уарр) на основе базового алгоритма сканирования для ЭУ и ОП, полученных после первого этапа;

• 3 этап: поиск точных значений координат экстремума в пределах новой области (области С) исходным, т. е. не прошедшим процедуру передискретизации, эталонным участком с помощью метода шаблонного поиска. Область С принадлежит исходной ОП с координатами центра (харр,уарр) и размером йхй (т<й<К). Данную область С будем называть в дальнейшем областью шаблонного поиска (ОШП). Алгоритм передискретизации работает следующим образом. Сначала производится разбиение ЭУиОПна субблоки размером кхк (к - целое). Параметр к будем называть коэффициентом понижения разрешения (коэффициент передискретизации). В пределах каждого субблока вычисляется среднее арифметическое значение яркости пикселей. Далее субблок заменяется на пиксель со значением полученного среднего арифметического. Таким образом, после передискретизации и округления до целого размеры полученного эталонного участка и области поиска равны (т/к)х(т/к) и (Л/к)х(Л/к), соответственно.

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

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

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

Выбор параметров расчета в зависимости

от текстуры изображений

Покажем влияние текстуры изображений на выбор значения коэффициента передискретизации и размера ОШП путём численного моделирования методом Монте-Карло. Восьмиразрядные изображения формировали искусственно следующим образом. Путем розыгрыша с помощью генератора случайных чисел получали два массива размером тхт (первое изображение) и ЯхЯ (второе изображение), где т<Я. Гистограмму яркости задавали нормальным законом распределения с математическим ожиданием, равным половине предельного значения яркости, и дисперсией а,2.

Для первого изображения к значению яркости каждого пикселя добавляли ошибку, отражающую погрешность съёмки, и полученный участок изображения «вставляли» во второе с некоторым сдвигом. Таким образом, сдвиг изображений, равный вектору смещения, был заведомо известен. Полученное изображение имитировало область поиска. Эталонный участок изображения формировали из первого изображения путем добавления ошибки. Ошибку 1е разыгрывали случайным образом при заданном размахе аддитивной погрешности с нормальным законом распределения. Здесь /e=Gauss(л=0,а■e=7), где л - математическое ожидание, ае - среднеквадратическое отклонение, характеризующее данную погрешность.

Сравнивая рассчитанное отклонение амплитуды вектора от известного, можно оценить влияние текстуры а на значение коэффициента передискретизации (рис. 3). Видно, что эти зависимости имеют пороговый характер. Если гистограмма яркости ограничена, то коэффициент передискретизации не может быть более 2. Для более качественных типов текстур этот коэффициент в принципе может находиться на уровне £=10-15.

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

Рис. 3. Зависимости ошибки нахождения вектора смещения Аи от коэффициента передискретизации к при варьировании дисперсии яркости а2. Н=301, т=Ю1

Оценим влияние текстуры изображения на выбор размера ОШП (рис. 4).

Рис. 4. Зависимости количества локальных экстремумов п от размера в при варьировании дисперсии яркости а2

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

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

Оценка выигрыша по времени

Оценим выигрыш по времени для модифицированного алгоритма по сравнению с базовым. Будем полагать, что время расчета одного текущего значения функционала приблизительно постоянно и равно Лt (такт). Тогда время сканирования по базовому алгоритму будет равно:

tba = m2(R - m +1)2 At.

(1)

m = tda + t- =-

-At + 9m log2(d -1)At.

Следовательно, выигрыш по времени составляет: (R - m +1)2 к4

V

600-

450-

300-

150

04

___+

/ < -d = 115

J и t/ = 113

J • -</=111 -

J ▲ -J = 109

* -d= 107

12

15

Первый сомножитель соответствует количеству тактов при некотором фиксированном положении ЭУ в пределах ОП. Второй сомножитель равен числу позиций эталонного участка в пределах области поиска при сканировании.

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

^ «(т / к )2(( Я / к) - (т / к) + 1)2 Лt.

Время, затраченное на сканирование по алгоритму шаблонного поиска, рассчитывается по формуле tш^9mг[1ogг(d-1)]Лt. Первый сомножитель определяет количество положений ЭУ в пределах ОШП для текущей итерации, а второй - количество тактов для текущего положения ЭУ. Третий сомножитель отражает количество всех итераций. Таким образом, общее время сканирования по модифицированному алгоритму можно найти по формуле:

т2(Я -т + к)2

V = ^ = -

(R - m + к)2 + 9k4log2(d -1)

Рис. 5. Зависимость выигрыша по времени модифицированного алгоритма по сравнению с базовым от коэффициента передискретизации при варьировании d

На рис. 5 показана зависимость выигрыша по времени V от коэффициента передискретизации при разных значениях размера области шаблонного поиска. Значения коэффициента передискретизации и размер ОШП определяли из двух предыдущих графиков (рис. 3, 4) следующим образом. При определенном характере текстуры выбирали предельные значения коэффициента передискретизации, когда величина отклонения Лu еще равна нулю. Аналогично находили предельный размер d, при котором функционал еще имеет один экстремум.

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

Выводы

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

Приведен способ выбора коэффициента передискретизации и размера области шаблонного поиска путем численного моделирования. Аналитически показан выигрыш по времени модифицированного алгоритма по сравнению с базовым.

СПИСОК ЛИТЕРАТУРЫ

1. Экспериментальная механика / под ред. А. Кобаяси. - М.: Мир, 1990. - Т. 1. - 615 с. - Т. 2. - 551 с.

2. Peters W.H., Ranson WF. Digital imaging techniques in experimental stress analysis // Optics and Lasers in Engineering. - 1982. -V. 21. - №3. - P. 427-431.

3. Li E.B., Tieu A.K., Yuen W.Y.D. Application of digital image correlation technique to dynamic measurement of the velocity field in the deformation zone in cold rolling // Optics and Lasers in Engineering. - 2003. - V. 39. - № 4. - P. 479-488.

4. Измаилов А.Ф., Солодов М.В. Численные методы оптимизации. - М.: Физматлит, 2005. - 304 с.

Поступила 15.12.2010 г.

4