Full Text

Введение Современный уровень развития технологий автоматической обработки, анализа и интерпретации визуальной информации обусловливает их широкое применение во многих отраслях: микробиологии, медицине, анализе структуры (в том числе наноматериалов), автоматизированных системах для выполнения широкого ряда производственных операций и др. [1]. В рамках экспериментальной механики также активно развивается и применяется ряд подходов, основанных на обработке и анализе изображений, получаемых в процессе нагружения объектов исследований. К настоящему времени в научно-технической литературе получил широкое распространение термин «оптический поток» («optical flow»): кажущееся движение яркостной картины, наблюдаемое при перемещении объектов перед камерой либо камеры в неподвижной окружающей обстановке. Работы по данному направлению активно развиваются с начала 80-х гг. XX в. Среди основных методов определения оптического потока выделяют: дифференциальные алгоритмы [2-7], фазовые (частотные) алгоритмы [8-12]; блочные алгоритмы [13-16]. Первые дифференциальные алгоритмы были предложены B.K. P. Horn и B.G. Schunck [4], а также B.D. Lucas и T. Kanade [5, 6]. Развитие алгоритмов в первую очередь было основано на усложнении модели движения, которая учитывает афинные преобразования, поворот [17-19]. Проблемы фильтрации векторных полей были затронуты в работах [20, 21]. Проблему ошибок определения поля оптического потока часто решают с использованием пространственной регуляризации, которая впоследствии была расширена до пространственно-временной [22-26]. Алгоритмы вычисления оптического потока активно реализуются с использованием параллельных вычислений [27-32]. В то же время известен ряд исследований, посвященных непосредственно алгоритмической оптимизации вычислений при определении оптического потока [33-35]. Также в работах [36-38] было рассмотрено применение многомасштабного иерархического подхода для обработки изображений (пирамиды изображений) при вычислении оптического потока. Также на расчете оптического потока основан и получивший широкое распространение в экспериментальной механике для оценки деформации твердых тел метод корреляции цифровых изображений (DIC - Digital Image Correlation). В основе его действия лежит построение векторов перемещений и дальнейшей расчет компонент деформации с помощью процедуры численного дифференцирования. В настоящее время метод корреляции цифровых изображений является одним из наиболее распространенных подходов к изучению процессов деформации и разрушения структурно-неоднородных материалов. Существенный вклад в его развитие внесли M.A. Sutton [13, 39-43], F. Hild, S. Roux [44-47], B. Pan [48-52] и др. В последние годы аппаратно-программные комплексы, принцип действия которых основан на корреляции цифровых изображений, получили широкое коммерческое распространение, прежде всего продукция Correlated Solutions Inc. (VIC 3D), LaVision GmbH (StrainMaster) и др. Указанные системы, как правило, оборудуются штатными оптическими системами, калибровка которых проводится перед проведением эксперимента в зависимости от условий видео (фото) наблюдения. Исследования с помощью указанных методов проводятся и в Центре экспериментальной механики ПНИПУ [53]. В то же время компьютерный код, поставляемый в составе таких систем, является закрытым, что не позволяет эффективно решать ряд прикладных задач. Поэтому рядом коллективов проводятся исследования по развитию собственных программно-алгоритмических решений, а также и аппаратных средств для изучения процессов деформации и разрушения, включая расчет стандартных параметров механики разрушения на основе метода корреляции цифровых изображений [46, 54-62]. В данной работе описан ряд задач подобного класса, решенных в последнее время в лаборатории механики полимерных композиционных материалов ИФПМ СО РАН. 1. Развитие аппаратных средств метода корреляции цифровых изображений В методе корреляции цифровых изображений (DIC) алгоритмы определения движения на серии изображений в общем виде основаны на минимизации функционала: (1) где w - функция формы; S - область изображения, по которой минимизируется E; ED - мера подобия блоков изображений; Es - мера подобия векторов в оптическом потоке; α - регуляризационный коэффициент. Меры близости ED двух участков текущего изображения (ТИ) и эталонного изображения (ЭИ) могут определяться различными функциями, условием которых является наличие экстремума, обычно минимума. Наиболее простой мерой близости является сумма квадратов разности [63, 64]: (2) где , - значения яркости пикселов сравниваемых участков ТИ и ЭИ соответственно; X - координаты пикселов участка S; w - функция формы. Одной из наиболее сложных мер является нормированная сумма квадратов разности с нулевым средним [65, 66]: (3) где и - среднеарифметические значения участков ТИ и ЭИ соответственно. Существует два пути нахождения минимума функционала E. 1. Блочный метод. Область S начального изображения сравнивается в скользящем окне со вторым изображением. Часто данный подход используется для предварительного поиска больших перемещений с пиксельной точностью. Однако он редко применяется при вычислениях с субпиксельной точностью и со сложной функцией формы, так как приводит к неоправданно большим вычислительным затратам. 2. Применение оптимизации. Используется какой-либо численный оптимизационный алгоритм, например алгоритм Ньютона-Рафсона или Гаусса-Ньютона. Оптимизационный подход, напротив, используется в субпиксельной области, поскольку обеспечивает устойчивую сходимость в диапазоне [-1…+1] пикселей. Оптимизацию методом Ньютона-Рафсона [42, 67] можно записать в следующем виде: (4Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.) где k - номер итерации; p - параметрический вектор, содержащий коэффициенты функции формы. Система уравнений, описывающая перемещения в области S, основанная на аффинной функции формы, имеет следующий вид [13, 67]: (5) Был реализован лабораторный макет системы технического стереозрения (СТСЗ), который включает две цифровые фотокамеры Canon EOS 700D (разрешение КМОП (CMOS) матриц 5184×3456 пикселов при физическом размере сенсора 22,3×14,9 мм), закрепленные на штанге с возможностью регулировки расстояния и угла между камерами. На рис. 1 приведен лабораторный макет системы технического стереозрения с камерами на расстоянии 40 см друг от друга и углом поворота каждой камеры относительно центральной оптической оси системы равным 10°. Камеры были подключены к ЭВМ с целью их синхронизации и получения изображений в память ЭВМ в процессе работы системы. Система позволяет записывать пару изображений (от левой и правой камеры) в текущий момент времени и может регистрировать серию изображений до 5 кадров в секунду при уставленных параметрах фокуса и выдержки. Разработанная система технического стереозрения для задач оценки деформации включает следующие этапы функционирования (рис. 2): калибровка камер, включающая поиск узловых точек калибровочного шаблона на изображениях; ректификация изображений стереопар; вычисление карт диспаратности; восстановление трехмерных координат точек стереопары; вычисление трехмерных перемещений точек; вычисление деформаций в терминах экспериментальной механики. Рис. 1. Лабораторный макет системы технического стереозрения Fig. 1. The laboratory scale stereo machine vision system Поиск узлов калибровочного шаблона? Калибровка камер (K, K’, R, t...) {opencv} Изображения стереопары в момент времени t i Ректификация изображений {opencv} Поиск перемещений опорных точек в момент времени t i Восстановление трехмерных координат в момент времени t i {opencv} Калибровочные стереопары Вычисление деформации Вычисление трехмерных перемещений Определение опорных точек на левом изображении в момент времени t 0 Изображения стереопары в момент времени t i+1 Поиск перемещений опорных точек в момент времени t i+1 Восстановление трехмерных координат в момент времени t i+1 {opencv} Рис. 2. Алгоритм работы системы технического стереозрения Fig. 2. The operation algorithm of the stereo machine vision system Сопоставление методик оценки деформации при 2D- и 3D-расчете. Сопоставление результатов расчета деформации в режиме моно- и стереозрения проведено на модели изгиба тонкой пластины, опертой по краям. Такая модель выбрана в первую очередь по причине присутствия в ней внеплоскостной деформации. Была получена последовательность модельных пар изображений, отражающих внеплоскостные перемещения. Оценивали среднюю абсолютную ошибку вычисления продольного и поперечного тензора деформации как для случая плоско-деформированной поверхности, так и для случая внеплоскостных деформаций. Среднюю абсолютную ошибку оценивали по формуле (6) где - заданное значение тензора деформации; - рассчитанное значение тензора деформации. Ошибки поперечной и продольной компоненты тензора обозначались и соответственно. Эксперимент показал, что система с одной камерой не позволяет достоверно определять деформацию при наличии внеплоскостных перемещений (рис. 3). Ошибка вычисления компонент тензора деформации в случае системы с одной камерой на два порядка больше, чем при использовании системы стереозрения. Кроме того, в случае системы с одной камерой ошибка на порядок больше самой величины деформации. а б Рис. 3. Зависимость ошибки определения компонент тензора деформации Eεx (а) и Eεy (б) от номера стереопары для системы технического зрения с одной (2D) и двумя камерами (3D) Fig. 3. Error graph of determining the Eεx (a) and Eεy (b) strain tensor components vs the stereopair number for a vision system with one (2D) and two cameras (3D) Пример восстановленной поверхности представлен на рис. 4. Рис. 4. Поле глубины, отображенное в координатах левого изображения при изгибе поверхности опертой по краям Fig. 4. The depth field displayed in the coordinates of the left image under bending of a surface supported over the edges Пример поля компоненты тензора деформации в случае изгиба поверхности, опертой по краям, представлен на рис. 5. Рис. 5. Поле компоненты εxx тензора деформации, отображенное в координатах левого изображения при изгибе поверхности, опертой по краям Fig. 5. The field of the εxx strain tensor component displayed in the coordinates of the left image under bending of the surface supported over the edges 2. Инкрементальный алгоритм оценки перемещений Определение пространственных координат заданных точек объекта по последовательности стереопар подразумевает определение положения этих точек на изображениях. Для этого необходимо знать положения этих точек в начальный момент времени. Определение перемещения точек между начальным и текущим изображением может быть затруднено по причине значительной деформации объекта, изменения условий съемки и т.д. При изменении положения объекта в пространстве на проекциях происходят значительные изменения. Поскольку между соседними изображениями в серии такие изменения значительно меньше, чем между начальным и текущим изображениями, то помехоустойчивость и точность определения перемещений значительно возрастает. По этой причине был предложен инкрементальный алгоритм, расширенный для обработки последовательности стереопар. Используется сквозная последовательная нумерация изображений стереопар (рис. 6). Рис. 6. Нумерация изображений стереопар Fig. 6. The principle of a stereopair image numbering Схема инкрементального алгоритма приведена на рис. 7, в ней используются обозначения, описанные ниже: 1. I0, I1, I2 и т.д. изображения, как на рис. 6. 2. VF01P, VF01S - поля векторов перемещений, посчитанные с пиксельной и субпиксельной точностью между изображениями I0 и I1. 3. Для обобщения п. 1, 2 и нумерации изображений и векторных полей в зависимости от текущего номера стереопары введены обозначения с фигурными скобками вида I{a}, VF{a}{b}P, VF{a}{b}S, где {a}, {b} - номера изображений пары в общем виде, например I{2(i-1)}, VF{2(i-1)}{2i}P. 4. VFP, VFS - процедуры определения полей перемещений с пиксельной и субпиксельной точностью соответственно. а б Рис. 7. Схема инкрементального алгоритма определения полей векторов перемещений для последовательностей стереопар Fig. 7. The scheme of the incremental algorithm for constructing the displacement vectors fields for stereopairs in image series Стереопара, соответствующая начальному моменту времени, обозначена как 0-я стереопара, изображения I0 и I1 (рис. 7, а). На первом этапе строится поле векторов перемещений с пиксельной точностью, по классическому алгоритму, основанному на быстром преобразовании Фурье и обозначенному на схеме как VFP. Затем поле векторов уточняется до субпиксельной точности (блок VFS) с использованием дифференциального алгоритма. Поле векторов перемещений VF01S между кадрами стереопары является картой диспаратности, а также используется как начальное приближение для определения перемещений в паре I1, I3 с пиксельной точностью. Для первой стереопары входное векторное поле в левой части схемы VF00S соответствует нулевому полю и не используется. Блоки VFP используют предыдущее и текущее изображения для определения поля векторов перемещений. Перемещения определяются отдельно для левых изображений I{2(i-1)}, I{2i} и правых I{2(i-1)+1}, I{2i+1} (рис. 7, б). Также используется векторное поле, рассчитанное для предыдущей стереопары, перемещения которого используются для определения положения ядра корреляции. Перемещения с субпиксельной точностью (блоки VFS) определяются по изображению I0 (начальное левое изображение) и текущему изображению стереопары, левому I{2i}, либо правому I{2i+1}. Для начального приближения двух первых элементов р0 и р1 искомого вектора P (7) используется векторное поле с пиксельной точностью VF{2(i-1)}{2i}P либо VF{2(i-1)+1}{2i+1}P (см. рис. 7, б). Вычитая поле VF0{2i}S из VF0{2i+1}S, получаем векторное поле VF{2i}{2i+1}S, которое соответствует карте диспратности i-й стереопары изображений I{2i} и I{2i+1}. Была сгенерирована серия модельных изображений, отражающих изменение рельефа, при одновременном растяжении поверхности материала. Для изменения характера рельефа каждое изображение формировалось из пары изображений двух модельных серий для двухосного растяжения при заданном приращении деформации. Построенные поля векторов перемещений (рис. 8, а, б) показывают значительное преимущество инкрементального алгоритма над обычным подходом. Зависимость коэффициента корреляции модельного и построенного векторного поля от номера кадра серии изображений (рис. 8, в) показывает, что при изменении рельефа поверхности более чем на 50 % от исходного «классический» алгоритм не справляется с задачей сопоставления изображений. Более того, при изменении рельефа на 75 % и более сопоставить начальное и текущее изображения не представляется возможным. Инкрементальный подход позволяет сопоставлять изображения во всем диапазоне приращения деформаций. а б в Рис. 8. Векторные поля, полученные по «классическому» алгоритму (а) и инкрементальному алгоритму (б). Зависимость коэффициента корреляции kr векторного поля от номера кадра серии изображений (в): 1 - классический алгоритм; 2 - инкрементальный алгоритм Fig. 8. Vector fields calculated by the classical (a) and incremental algorithms (b). Dependence of the correlation coefficient kr of a vector field vs the frame number of an image series (c): 1 is the classical algorithm, 2 is the incremental algorithm 3. Развитие алгоритма выделения трещины В рамках механики усталостного разрушения одной из ключевых является задача оценки параметров распространения трещины (длины, скорости роста, раскрытия берегов и т.д.), поэтому актуальной технической задачей является автоматизация процесса их измерения, что обусловливает необходимость слежения за ее вершиной. С позиции оценки сопротивления материалов зарождению и росту усталостных трещин [68] значительный интерес представляют задачи, когда длина трещины не превышает десятков микрон. В этом случае применение макрорежима видеофиксации не позволяет достичь требуемого пространственного разрешения [69]. Для решения указанной проблемы используются оптические микроскопы, оснащенные фото- или видеокамерой. Однако при больших увеличениях оптической системы в процессе роста (усталостной) трещины ее вершина будет постоянно выходить из поля зрения, поэтому необходимо создавать аппаратно-программные комплексы, позволяющие (в режиме реального времени) управлять механическим сканером, перемещающим микроскоп и видеокамеру относительно образца, на поверхности которого контролируется рост трещины. Ключевыми работами, направленными на автоматическое детектирование трещин на оптических изображениях, можно считать [70, 71]. За прошедшее время сформировалось направление прикладных исследований, нацеленных на автоматическое выделение трещин на изображениях. Значительный практический интерес представляло детектирование и количественная характеризация трещин на поверхности дорожного полотна [71-73], в то время как другая заметная часть работ посвящена оценке состояния инженерных конструкций и выявлению трещин в них [74-76]. В данном разделе работы решалась задача разработки и тестирования алгоритма детектирования трещины на оптических изображениях поверхности нагруженного материала и определения координат ее вершины. Алгоритм выделения трещины. Алгоритм детектирования трещины на серии изображений и определения координат ее вершины содержит следующие основные этапы: 1) анализ последовательности изображений; 2) выделение объекта максимальной площади, соответствующего области наибольших неравномерностей на картине оптического потока; 3) анализ индивидуальных изображений; 4) построение скелетона трещины; 5) поиск вершины трещины. Анализ последовательности изображений. При анализе серии изображений (либо видеопоследовательности), отображающих распространение трещины, ключевую информацию содержат данные о перемещениях участков (фрагментов) на них. В этом случае характер текстуры на поверхности важен лишь для повышения достоверности (помехоустойчивости и точности) определения перемещений. Для дальнейшего анализа используются данные об оптическом потоке: , где и - изображения в текущий и предыдущий моменты времени. Таким образом, объектом анализа является оптический поток, визуализация которого в виде цветового кодирования показана на рис. 9, а [77]. При этом разрыв или резкое изменение перемещений в его пределах (на берегах раскрывающейся трещины) свидетельствует о наличии трещины. Для определения оптического потока используется ранее описанный в подразд. 1 корреляционный алгоритм. Однако с точки зрения детектирования трещины алгоритм также должен обеспечивать устойчивость при поиске перемещений на ее границах (берегах). Данное требование, предъявляемое к подобного рода алгоритмам, обосновано необходимостью точнее выделять область интереса (здесь и далее под этим термином подразумевается фрагмент изображения, содержащий трещину). Для этого могут быть использованы алгоритмы, основанные на расчете взвешенной меры подобия блоков изображений, поскольку они позволяют в полной мере «сохранять» границы между объектами [78]. а б в Рис. 9. Результаты анализа последовательности изображений на основе определения поля перемещений: изображение оптического потока (а); изображение дифференцированного оптического потока (б); алгоритмически определенное расположение трещины (в) Fig. 9. Results of the image series analysis based on the constructing the displacement fields: optical flow pattern (a); differentiated optical flow pattern (b); algorithmically determined crack location (c) Поскольку у «коротких» трещин (длиной до десятков микрон) величина раскрытия может составлять единицы пикселов, с целью дополнительного контрастирования и подчеркивания границ применяется процедура дифференцирования поля оптического потока (рис. 9, б), что реализуется использованием оператора Собеля для окна размером 3×3 [79]. Выделение объекта максимальной площади. Для определения области интереса полученное изображение дифференцированного оптического потока подвергается нормализации (приведению к диапазону [0, 255]) и пороговой бинаризации. В качестве порога бинаризации принимается значение яркости, равное единице, т.е. «отбрасываются» только самые малые перемещения. Далее к бинаризованному изображению применяется операция морфологического замыкания [79] с целью заполнения мелких промежутков внутри сплошных объектов. Процесс завершается поиском контуров, оставшихся на изображении объектов, и оценкой их формы. В качестве трещины принимается объект наибольшей площади. Данный метод применяется для первичного (грубого) детектирования фрагмента изображения, наиболее вероятно содержащего трещину (рис. 9, в). Анализ индивидуальных изображений. При обработке одиночных изображений (без учета остальных изображений в серии/видеопотоке) в основе анализа данных лежит информация о текстуре поверхности. В процессе слежения за трещиной решаются две ключевых задачи: ее поиск на изображении и последующее определение местоположения ее вершины (детектирование и локация). Для этого проводится контурный анализ, основанный на разбиении изображения на отдельные объекты и фон, т.е. процедура сегментации [79]. Использование контурного анализа позволяет выделить объекты, максимально подобные трещине. Априори считается, что трещина (как область, соответствующая отсутствующему материалу - несплошности) должна являться наиболее темным участком (объектом) на изображении. С целью снижения вычислительных затрат и повышения помехоустойчивости работы алгоритма детектирования трещин используется процедура бинаризации. Это позволяет заметно сократить количество выделяемых и анализируемых объектов. Для бинаризации применяется алгоритм, в котором для определения глобального порогового значения яркости анализируется соответствующая гистограмма. Последующей операцией в обработке полученного черно-белого изображения является морфологическое замыкание [79]. Его целью является удаление мелких (темных) фрагментов на изображении (формо- и размерозависимая фильтрация), соответствующих элементам текстуры на поверхности, не относящихся к трещине. Затем полученное черно-белое изображение подвергается контурному анализу с целью детектирования объекта (объектов), геометрические характеристики которого максимально соответствуют таковым для трещины. В качестве трещины принимается объект максимальной площади. Построение скелетона трещины. Для поиска и определения координат вершины несплошности используется построение скелетона [80] объекта типа «трещина». Для этого оценивается кратчайшее расстояние от точки внутри трещины до ее границ. Применение данной операции основано на том, что скелетом несплошности является кривая, объединяющая точки, являющиеся равноудаленными от ее границ. Поиск вершины трещины. Для определения координат вершины несплошности проводится оценка величины раскрытия вдоль оси скелетона. Решение данной задачи упрощается тем, что все смещения уже были рассчитаны в первой части работы алгоритма. Искомая точка принадлежит участку, характеризуемому наименьшими смещениями. Тестирование алгоритма. Тестирование алгоритма проводили с использованием сгенерированной серии модельных изображений. На первом этапе функционирование предложенного метода исследовали в условиях отсутствия возмущающих воздействий (шум, размытие, имитирующее расфокусировку, и др.). Серию модельных изображений формировали в несколько этапов в виде текстурированной поверхности, имитирующей напыленный спекл [13, 81], на которую программным способом накладывалась трещина (рис. 10, а): 1) создание модели поверхности материала (формирование изображения спекла); 2) задание параметров несплошности (длина, раскрытие, форма берегов, имитирующая ее ветвление); 3) изменение модели материала (пересчет изображения) с учетом изменения параметров трещины. а б Рис. 10. Модельное (а) и экспериментально полученное (б) изображения, характеризующие распространение усталостной трещины Fig. 10. The model (a) and experimentally obtained image (b) characterizing the fatigue crack propagation Для определения устойчивости работы алгоритма на втором этапе тестирования было сгенерировано 8 пар модельных изображений, отражающих 4 варианта траектории роста трещины и 2 типа текстуры. Каждая траектория получена случайным разбиением линии раскрытия трещины. Длина последней составляла половину ширины изображения, а величина раскрытия между двумя соседними изображениями в последовательности, снятыми в моменты времени t и t + 1, составляла 30 и 60 пикселов. Использованные варианты текстуры различались размером элементов напыленного спекла, когда максимальный диаметр «пятен» равнялся 20 и 30 пикселов. Сгенерированные модельные изображения подвергали размытию и зашумлению (для имитации искажений, характерных для лабораторных экспериментов). Зашумление изображений проводили путем наложения на каждый пиксель (значения яркости которого лежат в диапазоне [0; 255]) случайного значения, имеющего нормальное распределение (что подобно характеру шумов для реальных видеодатчиков) и различное стандартное отклонение . Для размытия использовали нормальное Гауссово размытие со среднеквадратичным отклонением . Дополнительно проведены эксперименты по детектированию трещины при обработке экспериментальных изображений [69], зарегистрированных при циклическом растяжении образца алюминиевого сплава Д16 толщиной 3 мм (рис. 10, б). Для оценки эффективности работы предложенного алгоритма были использованы два критерия: 1. Ошибка определения местоположения вершины трещины, равная расстоянию между ее координатами на модельном изображении (с их известным значением) и таковыми, определенными с помощью предложенного алгоритма: (8Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.) где и - модельные и найденные координаты соответственно. Заметим, что причиной данной ошибки в первую очередь может являться потеря информации в области вершины трещины, когда алгоритм в процессе бинаризации «отсекает» фрагмент, также относящийся к несплошности. 2. Ошибка определения площади трещины, выделенной с помощью предложенного алгоритма (см. выделение объекта максимальной площади), рассчитываемая из отношения найденной площади (Sf) к ее величине, заданной в модели: (9) О повышении точности функционирования алгоритма должно свидетельствовать достижение минимального значения каждого из этих двух показателей. Результаты тестирования и их обсуждение. На рис. 11 представлены графики зависимостей ошибок определения параметров трещины от степени зашумления (s) (рис. 11, а) и размытия (рис. 11, б) исходных изображений. Показано, что для незашумленных изображений (см. рис. 11; при стандартном отклонении шума (СО) и среднеквадратичном отклонении (СКО) фильтра размытия, равных нулю) средняя ошибка определения площади трещины составила 1,93 %, в то время как средняя ошибка определения координат ее вершины равняется 56 пикселей. При этом стандартное отклонение определения площади трещины и координат ее вершины составляет 0,21 % и 15 пикселей соответственно (согласно приведенным на графиках величинам разброса анализируемых параметров). Зашумление изображений. Результаты анализа всех 8 пар изображений, представленных на графиках на рис. 11, а, показывают, что использование предложенного алгоритма в условиях значительного зашумления обусловливает максимальную величину ошибки при детектировании (области) трещины, dES = 5,2 % (кривая 1). При определении координат ее вершины (кривая 2) возрастание ошибки происходит практически линейно, и максимум данного параметра D = 90,2 (с учетом разброса) также наблюдается при максимальном СО шума . По мнению авторов, в рамках задачи слежения за трещиной при позиционировании оптической системы относительно поверхности образца с учетом поставленных в начале раздела технических требований подобные величины являются вполне приемлемыми. а б Рис. 11. Графики зависимости ошибок определения площади трещины и координат ее вершины от уровня шума (а) и размытия (б); 1 - ошибка определения площади трещины (dES); 2 - ошибка определения местоположения вершины трещины (D) Fig. 11. Errors graphs on the crack area and crack tip coordinates determination vs the noise level (a) and blurring (b); The determination error of the crack tip location (D) - 2; the crack area (dES) - 1 Размытие изображений. Размытие изображений более критично влияет на величину ошибки (см. рис. 11, б), которая возрастает в значительно большей степени и достигает dES = 36 % при размытии фильтром с (кривая 1). При определении координат вершины трещины (кривая 2) возрастание ошибки происходит подобным образом; максимум данного параметра D = 280,2 также наблюдается при максимальном СКО фильтра размытия . В постановке задачи исследований было указано, что точность определения участка, содержащего вершину трещины, должна составлять порядка 20 % от горизонтального размера растра, т.е. при изображениях размером не менее 2000×1000 пикселей максимально допустимая величина ошибки определения координат составляет 400 пикселей. В случае всего исследованного диапазона зашумления (см. рис. 11, а) этот порог не оказывается превышенным, в то время как при размытии уже при ошибка при работе алгоритма превышает допустимое значение (см. рис. 11, б). Это обусловлено тем, что большая часть трещины имеет расстояние между берегами большее, чем апертура фильтра размытия. Данное обстоятельство, в свою очередь, накладывает ограничения на минимальное значение раскрытия трещины, которое может быть определено с помощью предложенного алгоритма. С другой стороны, введено ограничение на максимальную ошибку определения координат вершины трещины, не превышающую 15 % от их истинного положения. Результаты детектирования трещины в условиях зашумления полностью удовлетворяют данному условию, так как ошибка не превышает dES £ 5,5 %. В то же время, при имитации расфокусировки, предельное значение этой ошибки достигается уже при . Это означает, что алгоритм может быть эффективно использован только для случаев небольших расфокусировок оптической системы. Результат выделения трещины на экспериментальных изображениях представлен на рис. 12. Видно, что алгоритмически выделенный контур четко повторяет границы трещины, что свидетельствует в пользу возможности применения предложенного алгоритма для решения практических задач выделения трещин на реальных оптических изображениях. Рис. 12. Выделенная трещина на экспериментальном изображении Fig. 12. A crack detected on the experimental image 4. Использование DIC для изучения процессов деформирования и разрушения образцов с трещиной при статических и циклических нагрузках Несмотря на более чем столетнюю историю исследований, многие аспекты усталостного разрушения материалов до сих пор остаются не до конца изученными. В частности, пристальное внимание уделяется процессам, происходящим впереди вершины трещины на микроуровне, а именно формированию пластических монотонной и циклической зон, их эволюции в процессе распространения трещины, особенно в случаях неоднородного нагружения. Одним из наиболее распространенных методов изучения усталостного разрушения является фрактографический анализ изображений поверхностей изломов при помощи оптического и растрового электронного микроскопа. Однако этот метод имеет ряд существенных недостатков, прежде всего получение результатов только на уже разрушенном образце и, как следствие, невозможность мониторинга эволюции трещины и количественной оценки параметров ее роста во время циклического нагружения, а также необходимость дорогостоящего оборудования, в частности растрового электронного микроскопа, и высококвалифицированного оператора. При этом данный метод не позволяет напрямую количественно измерять многие параметры, используемые в механике разрушения. Другим известным способом исследования процессов роста трещины является применение датчика раскрытия трещины (в англоязычной литературе - crack opening displacement, COD), который также не лишен ряда недостатков. Например, разрешающая способность COD-датчика не позволяет получить оценку изменений, происходящих в материале при неоднородном нагружении, и для решения подобных задач данный подход не эффективен. Перспективным способом преодоления описанных ограничений является применение метода корреляции цифровых изображений, который способен измерять механический отклик материала впереди вершины трещины при приложении как единичных циклов перегрузки, так и при случайном спектре нагружения. Однако для успешного применения данного метода необходима дальнейшая разработка аппаратно-программных средств с целью а) повышения точности и стабильности получаемых результатов (количественных оценок); б) разработки методик подготовки поверхности и регистрации изображений поверхности для дальнейшего расчета полей деформаций; в) разработки методов расчета параметров, используемых в механике разрушения, на основе построения полей деформаций для проведения последующей оценки состояния материала и прогнозирования остаточного ресурса. Успешные исследования в области определения величины закрытия трещины при помощи корреляции цифровых изображений в процессе циклических испытаний (при наличии единичного цикла перегрузки, для образцов разной толщины, с целью исследования механизма распространения трещины) продемонстрированы в работах [54, 55, 82, 83]. Другой широко используемый подход основан на предположении о том, что на процесс роста трещины оказывает влияние зона пластичности и ее эволюция. Концепция зоны пластичности служит для объяснения эффекта замедления трещины вследствие приложения единичного цикла перегрузки во время циклического испытания с постоянной амплитудой. Аналитическая формулировка процесса роста трещины на микроуровне была разработана авторами [84]. Если при помощи традиционно используемых моделей закон роста трещины формулируется в зависимости от количества циклов нагружения da/dN, то описываемый подход основан на определении кинетики роста трещины в каждый конкретный момент времени в процессе нагружения, т.е. с использованием параметра da/dt. Основным предположением данной модели является то, что рост трещины зависит от взаимодействия монотонной и циклической пластических зон, размеры которых зависят от уровня закрытия трещины. Получение данных о распределении деформаций и напряжений вблизи вершины трещины при помощи оптических методов проводилось рядом авторов [56, 85-87]. Дальнейшим продолжением развития метода корреляции цифровых изображений является использование его не только для измерения величины закрытия или построения полей деформаций, но и для расчета параметров механики разрушения, таких как J-интеграл и коэффициент интенсивности напряжений K. J-интеграл, предложенный независимо Г.П. Черепановым и Дж. Райсом [88, 89], определяется как контурный интеграл, не зависящий от пути интегрирования, который характеризует интенсивность выделения энергии при распространении трещины. Обычно J-интеграл определяется по любому замкнутому контуру, включающему в себя вершину трещины. Преимуществом применения J-интеграла в расчетах является возможность его использования в случаях больших пластических деформаций материала (когда размер пластической зоны превышает 20 % от длины трещины). В работах [58-60, 90] авторы использовали метод DIC для расчета полей деформаций на поверхности материала, а затем импортировали это поле в конечно-элементную модель для последующего расчета J-интеграла. Определение КИН с использованием DIC может быть выполнено с использованием трех подходов: 1) на основании полей смещений, полученных методом DIC, совмещенных с результатами решения степенных рядов Вестергарда (Westergaard) и Вильямса (Williams) [57, 61, 91-93]; 2) на основании показаний программных датчиков деформаций в методе DIC, размещаемых на берегах раскрывающейся трещины [54, 62, 94]; 3) на основе расчетов J-интеграла [46, 95, 96]. 4.1. Методика расчета параметров механики усталостного разрушения при помощи метода корреляции цифровых изображений Механический отклик материала. Одним из параметров, позволяющих качественно охарактеризовать поведение материала вблизи вершины трещины, является петля механического гистерезиса, где по оси Oy откладывается величина нагрузки, определенная по датчику нагрузки ИМ, а по оси Ox - величина раскрытия, определяемая по программному экстензометру, размещенному на поле смешений (DIC). Диаграмма нагружения и разгрузки принимает вид петли вследствие нелинейности отклика материала (рис. 13). При этом увеличение доли пластической деформации ведет к увеличению «ширины» петли, тогда как ее наклонный характер вызван влиянием открытия/закрытия трещины. Описываемая петля механического гистерезиса состоит из трех участков: · начальный линейный участок, когда трещина закрыта; · участок нелинейности, вызванный постепенным открытием трещины на полуцикле нагрузки и ее постепенным закрытием на полуцикле разгрузки; · конечный линейный участок, когда трещина полностью открыта. Таким образом, анализируя внешний вид петли, можно давать качественную характеристику происходящим в материале процессам. Численной оценкой может выступать максимальное значение деформации (εmax), достигаемое в каждой петле, что косвенно характеризует раскрытие трещины через податливость образца. Оценка уровня открытия/закрытия трещины. Один из важных численных параметров, непосредственно влияющих на скорость распространения трещины, особенно на начальном этапе ее роста и в режиме Пэриса, является уровень открытия/закрытия трещины, который можно определить как эффективную асимметрию цикла: (10) где Pop - нагрузка, при которой трещина открывается; Pmax - максимальная нагрузка в цикле. В идеальном случае, когда закрытие трещины отсутствует, петля механического гистерезиса должна иметь симметричную форму. То есть если провести вертикальною или горизонтальную прямую, которая будет делить петлю на две равные по площади части, то эта прямая и будет осью симметрии петли. При закрытии площадь нижней части петли станет меньше площади верхней. На этом предположении был разработан метод оценки уровня закрытия усталостной трещины на основе расчета площади петли механического гистерезиса. Процедура расчета (см. рис. 13) уровня закрытия состоит: 1) из определения площади, очерченной петлей; 2) из определения такого уровня нагрузки Pmean, линия которой «разрезала» бы петлю на две равные по площади части - S1 и S2; 3) из расчета эффективной асимметрии цикла (Reff), которая находится из выражения Рис. 13. Схема расчета уровня открытия/закрытия трещины Fig. 13. The scheme illustrating calculation principle of crack opening/closure level Оценка параметров роста усталостной трещины после единичной перегрузки. Для проведения циклических испытаний использовался образец алюминиевого сплава Д16Т. Форма и основные размеры образца представлены на рис. 14. Рис. 14. Форма образца. Толщина 3 мм Fig. 14. Specimen’s shape and dimension. Thickness 3 mm Данный образец имеет боковой надпил, в вершине которого происходит рост усталостной трещины. Перед началом испытаний трещина была предварительно выращена на длину 1,7 мм. Область фотографирования составляла 2,5 мм, что составляет в физическом масштабе 1,04 пикселей/мкм. Схема нагружения представлена на рис. 15 и состоит из следующих этапов: 1) производится циклическое нагружение в базовом диапазоне нагрузок (0,25-2,5 кН) до тех пор, пока трещина не достигнет величины 2,5 мм (т.е. прирост должен составить 0,8 мм); 2) после проращивания трещины прикладывается цикл нагружения, в течение которого камера производит фотографирование поверхности образца с частотой 3 кадра/с, этот цикл - опорный; 3) затем прикладывается перегрузка с амплитудой на 50 % больше базовой, до 3,75 кН, во время которой также регистрируются изображения; 4) после приложения перегрузки происходит дальнейшее циклическое нагружение на базовом уровне, однако в отдельных циклах запускается регистрация фотографий, чтобы проследить изменение характера роста усталостной трещины. Номера циклов, в которых производится фотографирование, после перегрузки: 1, 2, 3, 10, 25, 100, 500, 4000, 10000, 15000, 24000, 28000 и 37000 циклов. Рис. 15. Схема нагружения для испытания с единичной перегрузкой. Параметры базового цикла Pmax=2,5 кН, Pmin=0,28 кН, асимметрия R=0,1. Цикл перегрузки POvL=3,75 кН. Длина трещины a=14,5 мм. Коэффициенты интенсивности Kmax= 12,86 МПа√м, Kmin= 1,29 МПа√м, KOvL=19,3 МПа√м Fig. 15. The loading sequence scheme for testing with a single overload. The parameters of the basic cycle Pmax = 2.5 kN, Pmin = 0.28 kN, asymmetry R = 0.1. The overload cycle POvL = 3.75 kN. The crack length is a = 14.5 mm. The stress intensity factors Kmax = 12.86 MPa√m, Kmin = 1.29 MPa√m, KOvL = 19.3 MPa√m Циклы, в которых регистрировались фотографии образца, имели частоту нагружения 0,05 Гц с одновременным фотографированием с частотой 3 кадра/с. Таким образом, записывали 60 кадров в течение цикла «нагрузки-разгрузки». Между циклами измерения продолжалось «базовое» циклирование, но с большей частотой - 5 Гц. Полученные серии испытания изображений далее обрабатывают при помощи метода корреляции цифровых изображений с целью вычисления полей деформаций и параметров механики разрушения. Результаты исследований и их обсуждение. Параметром, характеризующим изменение реакции материала на прикладываемую нагрузку, является максимальное значение деформации в вершине трещины - εmax, представленное на графике рис. 16. Видно, что εmax после перегрузки достаточно быстро снижается, пока трещина не достигнет длины в 20% rOvL, где стесненность материала вблизи вершины трещины максимальная - материал деформируется на минимальную величину при постоянных прикладываемых нагрузках. Затем следует достаточно быстрый рост значений εmax до уровня выше начального (до перегрузки). Результаты оценки уровня открытия/закрытия трещины представлены на рис. 17. Горизонтальной линией отмечено значение асимметрии при базовом уровне нагрузок R=0,1, т.е. значения меньше 0,1 означают, что трещина полностью открыта в течение всего цикла. До перегрузки закрытие трещины минимально и составляет всего 0,115 от Pmax. Перегрузка вызывает притупление вершины трещины, вследствие чего трещина остается полностью открытой во время всего цикла нагрузки-разгрузки. Однако, как уже отмечалось выше, открытие и закрытие трещины процесс протяженный во времени, поэтому точки, лежащие вблизи R = 0,1, могут оказывать влияние на скорость роста трещины, что и наблюдается в данном случае. Максимальный уровень закрытия достигается при ≈30% rOvL. При дальнейшем росте трещины закрытие снижается и остается на уровне немного выше базового (до перегрузки). Рис. 16. Зависимость максимальной деформации εmax от расстояния, на которое трещина проросла в пластическую зону, вызванную перегрузкой Fig. 16. Dependence of the maximum strain εmax vs the crack extension into the overloading induced plastic zone Рис. 17. Зависимость Reff от расстояния, на которое трещина проросла в пластическую зону, вызванную перегрузкой Fig. 17. Dependence of the Reff vs the crack extension into the overloading induced plastic zone 4.2. Развитие алгоритма измерения J-интеграла J-интеграл [88, 89] представляет собой независимый от пути (path-independent) линейный интеграл, который характеризует напряжения и деформации в вершине трещины, возникающие под действием упругих сил. Расчет J-интеграла происходит по заданному контуру интегрирования, который должен охватывать вершину трещины. Одним из преимуществ применения J-интеграла является возможность использования измеренных полей перемещений в областях, расположенных достаточно далеко от вершины трещины, в которых преобладают линейные упругие деформации и которые с достаточно высокой точностью могут быть определены при помощи DIC. Поэтому в расчетах могут быть использованы принципы Linear Elastic Fracture Mechanics (LEFM) и из уравнений теории упругости можно вычислить соответствующие напряжения и деформации, необходимые при оценке J-интеграла. В представленном исследовании описан алгоритм измерения J-интеграла с использованием метода корреляции цифровых изображений, и проведена оценка помехоустойчивости расчета в рамках применения метода корреляции цифровых изображений. Процесс измерения J-интеграла и выделение пластической зоны на базе вычисления оптического потока при механических испытаниях материалов можно разделить на четыре этапа: 1) расчет полей перемещений, основанный на сопоставлении участков изображений, оригинальный алгоритм построения которых развит и описан авторами в работах [81, 97]; 2) постобработка полей перемещений (сглаживание оптического потока); 3) расчет компонент деформации и напряжений; 4) вычисление J-интеграла и зоны пластичности в пределах построенного векторного поля. Сглаживание оптического потока. Наличие шумов на изображении, искажение поверхности, неоднородность освещения и пр. обусловливают появление ошибок при оценке перемещений, что, в свою очередь, приводит к значительной зашумленности полей деформации вследствие применения процедуры численного дифференцирования. Частично данная проблема может быть решена посткорректировкой поля векторов перемещений (пространственная фильтрация и сглаживание) [98, 99]. Расчет компонент деформации, напряжений и J-интеграла. Дифференцируя векторное поле, можно вычислить компоненты деформации [100]. Выражения для продольной εxx, поперечной εyy и сдвиговой εxy компонент тензора дисторсии имеют вид (Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.11) где , - компоненты вектора перемещений участков изображения по осям x, y. Вычисление производных проводили с привлечением процедуры численного дифференцирования. Безразностная формула для его реализации имеет вид [101] (Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.12) где k - номер отсчета функции; n - величина, определяющая уровень приближения; - величина, рассчитываемая по следующей формуле: (Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.13) где m - натуральное число и t - произвольное действительное число. При известных смещениях, деформациях и напряжениях J-интеграл может быть рассчитан по известной формуле [88, 89] (Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.14) где w - плотность энергии деформации, (15) Т - вектор нагрузок, ; Г - контур произвольной формы, охватывающий вершину трещины; n - вектор нормали к Г; σ, ε, u - напряжения, деформации и поля перемещений соответственно. В рамках линейно-упругой механики и в условиях плосконапряженного состояния тензор напряжений можно определить на основе закона Гука из тензора деформации: (Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.16) где D - матрица упругих модулей при плоском напряженном состоянии; E - модуль упругости; v - коэффициент Пуассона. Таким образом, когда известны поля тензора напряжений и деформаций, J-интеграл можно оценить из уравнения (14) следующим образом [102]: (Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.17) Здесь Г - контур, охватывающий вершину трещины, обход по которому осуществляется против часов стрелки; , - компоненты вектора нормали к Г; u, v - компоненты вектора перемещения. Для верификации исследуемого метода измерения J-интеграла использовался программный комплекс ANSYS для численного моделирования процесса статического растяжения образца с заданными геометрическими размерами (рис. 18). С помощью ANSYS были получены модельные поля перемещений (рис. 19), деформаций и напряжений, соответствующие величинам прикладываемой нагрузки, а также рассчитаны значения J-интеграла (таблица, рис. 20). Рис. 18. Расчетная модель образца (расчет проведен с помощью пакета ANSYS) Fig. 18. Calculated model of a specimen (the calculation was carried out using the ANSYS package) а б Рис. 19. Результаты моделирования в ANSYS: directional deformation (X axis (а), Y axis (б)) Fig. 19. ANSYS simulation results : directional deformation (X axis (a), Y axis (b)) Оценка ошибки вычисления J-интеграла. Для получения количественной оценки «точности» и помехоустойчивости расчета J-интеграла предложено использовать следующую расчетную величину: (Ошибка! Используйте вкладку "Главная" для применения к тексту, который должен здесь отображаться.18) где - величина J-интеграла, полученная с помощью пакета ANSYS; - рассчитанное значение J-интеграла. Из модельных перемещений по формуле (2) были рассчитаны значения J-интеграла (таблица, рис. 20). Видно, что с увеличением нагрузки возрастает отклонение измеряемой величины от модельных данных. Ошибка увеличивается от 0,06 до 3,95 %, при этом средняя ошибка для используемых величин нагрузки составляет 1,75 %, что говорит о применимости метода корреляции цифровых изображений для измерения J-интеграла. Значения J-интеграла и ошибки вычисления J-интеграла The values and calculation errors of the J-integral Величина Нагрузка, кН 2 3 4 5 6 7 8 9 10 11 12 13 , Н/м 463 1043 1868 2941 4276 5885 7789 10011 12554 15461 18765 23698 , Н/м 463 1041 1859 2919 4230 5801 7646 9783 12224 14994 18133 22761 Err, % 0,06 0,18 0,45 0,74 1,07 1,43 1,84 2,28 2,62 3,02 3,36 3,95 Рис. 20. График изменения величины J-интеграла при возрастании прикладываемой нагрузки: 1 - расчет при помощи пакета ANSYS; 2 - по предлагаемой методике Fig. 20. The J-integral variation with increasing applied load: 1 - calculation with the ANSYS; 2 - calculations based on the proposed method Заключение В рамках развития метода корреляции цифровых изображений разработана компьютерная система стереозрения и проведено тестирование на сериях стереопар, отражающих изменение положения объекта в пространстве, плоскостные и внеплоскостные деформации. Показано, что ошибка определения пространственных координат не превышает 0,75 ед. при угле между осями камер 20°. Предложенные алгоритмы позволяют определять плоскую и внеплоскостную деформацию с ошибкой, не превышающей 4·10-3. Ошибка вычисления компонент тензора деформации в случае системы с одной камерой на два порядка больше, чем при использовании системы стереозрения в случае внеплоскостных деформаций. Предложен модифицированный инкрементальный алгоритм определения перемещений на серии стереопар, который позволяет оценивать большие перемещения при последовательной обработке изображений. Представлен алгоритм вычисления деформации на поверхности пространственного объекта с использованием восстановленных профилей поверхности объекта. Предложен и протестирован алгоритм автоматического детектирования трещины на оптических изображениях и определения координат ее вершины. Результаты тестирования с использованием модельных изображений (включая размытие и зашумление) показали, что алгоритм позволяет определить координаты вершины трещины со средней ошибкой порядка 56 пикселов. Средняя ошибка определения площади трещины (включая выделение ее берегов) не превышает 1,93 %. Результаты тестирования алгоритма на экспериментальных изображениях показали возможности его эффективного применения для решения задач выделения трещин на реальных (оптических) изображениях. Проведена количественная характеризация процесса роста усталостной трещины в металлических сплавах с помощью методики, основанной на методе корреляции цифровых изображений (в терминах механики разрушения), включающей расчет: а) скорости роста усталостной трещины da/dN; б) максимального значения деформации (εmax) и в) эффективной асимметрии цикла (Reff). При испытании алюминиевого сплава Д16АТ показано, что в режиме Пэриса при единичной перегрузке минимальная скорость роста трещины достигается при длине трещины 30 % от размера пластической зоны, сформированной циклом перегрузки, и вызвана закрытием трещины, увеличивающим эффективную асимметрию цикла до Reff = 0,32. Развит алгоритм измерения J-интеграла с использованием метода корреляции цифровых изображений, заключающийся в измерении векторов перемещений, вычислении деформаций и напряжений с последующим расчетом значения J-интеграла. Его тестирование показало, что отклонение рассчитываемых значений J-интеграла от модельных составляет в среднем 1,75 % при минимальной величине отклонения 0,06 % и максимальной 3,95 %.

About the authors

P S Lyubutin

Institute of Strength Physics and Materials Science

S V Panin

Institute of Strength Physics and Materials Science; National Research Tomsk Polytechnic University

V V Titkov

Institute of Strength Physics and Materials Science

A V Eremin

Institute of Strength Physics and Materials Science

R Sunder

National Research Tomsk Polytechnic University

References

Statistics

Views

Abstract - 611

PDF (Russian) - 295