Метод анализа эволюции картин разрушения нагруженного мезообъема для идентификации формирующихся фрагментов Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 13, № 1, 2008

Метод анализа эволюции картин разрушения нагруженного мезообъема для идентификации формирующихся фрагментов*

М.С. Кириченко, И.Ю. Смолин, С. В. Панин Институт физики прочности и материаловедения СО РАН, Томск, Россия

e-mail: smolin@ispms.tsc.ru

A method for analysis of failure patterns of solids under loading is developed to identify forming fragments and to study features of their generation. 2D patterns exemplifying generation of discontinuities (cracks) in coal mesovolumes that consist of several ingredients and pores are considered as the images for analysis. The patterns were obtained with the help of numerical modeling of deformation and fracture of selected coal mesovolume. The first stage of the method is aimed to mark out the closed boundaries as an extension of the existing cracks, which determine failure fragments (dust particles). At the second stage, it is necessary to calculate and plot distributions of failure fragments according to their sizes. The necessity of pattern recognition arises due to the complex structure of the analyzed patterns since the most of the fragments are not enclosed by cracks and hence it is practically impossible to plot the required distributions. It is shown that the suggested approach allowed to reveal failure fragments and to estimate the characteristics of the particle-size distribution in the range of 1-20 ^m. The method developed is adaptable for other similar problems.

Введение

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

* Работа выполнена частично в рамках проекта 7.11.1.6 фундаментальных исследований СО РАН на 2007-2009 гг., междисциплинарного интеграционного проекта СО РАН № 90, проекта Министерства образования и науки по поддержке ведущих научных школ России НШ-394.2006.1 "Школа академика В.Е. Панина", Фонда содействия отечественной науке.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2008.

Механическое поведение твердых прочных сред (например, металлических или керамических материалов, горных пород и массивов) — сложный многоуровневый процесс. Его изучают на разных пространственных масштабах как экспериментально, так и теоретически, при этом большое значение имеет компьютерное моделирование процессов деформации и разрушения. Практическую значимость подобных работ сложно переоценить, поскольку даже приблизительные расчеты позволяют получить информацию, способную внести изменения в существующую технологию или предложить новые подходы. Например, можно спрогнозировать поведение горного массива в области извлечения полезного ископаемого. Появляется возможность предсказать поведение такого рода объектов в сложных условиях нагружения, что зачастую не удается осуществить в процессе натурного эксперимента.

В последние годы ИФПМ СО РАН (Томск) и НЦ ВостНИИ (Кемерово) проводят совместные исследования процессов деформации и разрушения угля в процессе его добычи. Одна из задач исследований — оценка характерного размера мелких фракций угля, формирующихся в процессе его разрушения.

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

С целью выработки подходов к решению этой задачи были выполнены расчеты на примере стругового метода добычи для конкретного состава угля в пласте [3, 4]. При добыче струговой установкой разрушение осуществляется одновременно несколькими резцами одного размера, расположенными примерно на одинаковом расстоянии друг от друга. Моделирование процесса резания угля производилось в двухмерной постановке в двух взаимно перпендикулярных плоскостях по ходу резания. Анализ полученных на этом этапе результатов позволил сделать вывод о том, что для зоны основного пы-леобразования характерным видом нагружения является комбинация сжатия и сдвига. Для определения фракционного состава пылевых частиц размером 1-100 мкм проведены расчеты деформации и разрушения характерных мезообъемов угля в условиях одновременного действия сжатия и сдвига.

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

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

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

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

1. Описание метода анализа изображений

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

Серия

компьютерных изображений

1 Получение бинарного образа

2 Идентификация объектов бинарных образов

3 Построение обобщенного бинарного эквивалента ^

4 Выделение объектов обобщенного бинарного эквивалента и вычисление их геометрических признаков -

Рис. 1. Упрощенная схема алгоритма обработки изображений

— получение бинарного образа;

— идентификация объектов бинарных образов;

— построение обобщенного бинарного эквивалента (создание суммарной картины разрушения образцов);

— выделение объектов (фрагментов разрушения) обобщенного бинарного эквивалента и вычисление их геометрических признаков.

1.1. Получение бинарного образа

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

Процедура бинаризации палитровых изображений выполнялась здесь нерекурсивным методом, т. е. с фиксированным порогом, — так называемая пороговая бинаризация по заданным цветам. Бинарный образ и исходное изображение имеют одинаковый размер. Алгоритм пороговой бинаризации выглядит следующим образом:

Г 0, если = 1р,

[ 1, если = 1р,

где БШ [г,]) — пиксел бинарного образа, Б [г,]) — соответствующий пиксел исходного изображения, 1в(г,о) — значение яркости исходного изображения, 1Р — пороговое значение яркости [6].

1.2. Идентификация объектов бинарных образов

Задача идентификации заключается в распознавании и наблюдении на различных изображениях меток, относящихся к одним и тем же по своей физической сути объектам. Наиболее известные методы идентификации объектов: метод трасс, зонно-комбинатор-ный и корреляционный методы, метод "редкой сетки", квазикорреляционный метод. Использование метода трасс в нашем алгоритме обусловливается его эффективностью при заданных расчетных характеристиках моделирования (отсутствие помех, малые нелинейные искажения объектов, отсутствие взаимного сдвига между ними) [7].

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

Под объектом здесь и далее будем понимать совокупность пикселов каждого следующего бинарного образа серии Б^, которые имеют одинаковую интенсивность и являются связанными, £ = 1, 2,... , где ^ — исходное количество изображений серии.

Под идентификацией понимаем совмещение объектов, присутствующих на каждом следующем бинарном образе серии Б^ изображений относительно предыдущего и установление их взаимно-однозначного соответствия. Для реализации алгоритма идентификации существенно наблюдение объектов в каждом из образов [7].

БШ [г,])

Рис. 2. Картины разрушения в мезообъеме, состоящем только из витринита, полученные с шагом по времени 0.2 мкс, приращение по компоненте деформации сжатия Аеуу = -0.128 % и по компоненте деформации сдвига Аеху = 0.064 %; замкнутые области, "залитые" серым цветом — поры; расчетная сетка — 574 х 626 ячеек; размер квадратного сечения мезообъема 57.4 х 62.6 мкм

Для идентификации удалим все объекты первого бинарного образа серии из каждого следующего. Таким образом, получим картину изменения с некоторой дискретой. Удаление объектов осуществляется простым вычитанием первого бинарного образа серии из последующих. В результате получим последовательность (серию) , состоящую из Л = ^ — 1 бинарных образов, каждый из которых представляет собой результирующую картину появления новых объектов.

Задача идентификации состоит в определении "траектории" появления объектов на каждом следующем бинарном образе относительно предыдущего.

Математическая модель задачи идентификации выглядит следующим образом [8]. Под меткой будем понимать точку (пиксел) изображения или совокупность точек. Если метки бинарного образа (кадра) пронумеровать в некоторой последовательности, то любой метке этого образа, обусловленной тем или иным объектом, можно поставить в соответствие тот же номер, что и метке следующего кадра 5п+1 серии. В данном случае максимальное правдоподобие обеспечивается таким объединением меток соседних бинарных образов, которое сообщает минимум функционалу вида

) = Е ЦЪ'Я)' (!)

г

где Яг и Яц — соответственно радиус-векторы г-й и ]-й меток бинарных образов Бп и 5п+1; Ь(Яг, Я) = Ьц — функция стоимости объединения меток г и ] в пару (г, ]), равная квадрату расстояния между объединяемыми метками.

При этом имеем в виду, что

M,

min Ф(г, j) = min V L(R,,Rj), (2)

(ij)en (,j)en^—J

г=1

где Mn — количество меток бинарного образа Sn, П — множество пар (i, j), удовлетворяющее условию не более чем попарного соединения меток Mn, т. е. если

(i,j) G П, то (i,k), (l,j) £ П; i,l = 1, 2,... , Mn; j, k = 1, 2,...,M,J+i; k = j; l = i. (3)

Алгоритм поиска абсолютной минимали функционала (2) при условии (3) имеет вид

M,-1

/mi^f^^'j) = min LMn j + miM 1 V Lj + (i,j)en jeMn+1 (¿,j)enMn-1 ^

+ min< min Lm„k - min Lm„j, min V" Lik - min L, > , (4)

1 • em^-1 jeiMn+i v (1>k)enMn-i ^ (U)enMn-1 ^ '

где Мп+1 — множество номеров меток бинарного образа 5п+1; ПМп 1 — множество пар

(l,j), удовлетворяющее условию (3) при l = 1, 2,... , Mn-1; j = 1, 2,... , Mn+1; M^+l 1 — то же, что и при Mn+1, при исключении номеров меток бинарного образа Sn+1, входящих в множество ПМп-1; Пмп-1 — множество пар (l, k), аналогичное множеству пар П, при исключении меток, входящих в пару (Mn, arg min^M +i Lm,j).

Представив в виде равенства (4) величины

M, -1 M,-2 2

min у L,, min у L,,..., min > L,,

(i,j)enMn-1 ^ (¿,j)enMn-2 ^ (¿,j)en2 ^

получим

m,

min j) = у min L,-, +

(i,j)e^ ¿^jeMn+1 j

,=1

M,-1 ( г i ^

+ V^ min < min Li+1jk - min L,+1,j ; min V^ L№ - min V^ L, > , (5) t! \keMi^+1 jeiMn+1 (1,k)eni ^ (и)еШ ^

где Пг — множество пар (l,j), удовлетворяющее условию (3), при l = 1, 2,...,i; j = 1, 2,... , Mn+1; MMn+1 = MMn+1 \ Мг, MM — множество номеров меток бинарного образа Sn+1, входящих в множество Пг; Пг — множество пар (l, k), удовлетворяющее условию (3) при l = 1, 2,... , i; k i Mn+1 \ arg minjeiMn+1 L,+1;j.

Определенные таким образом траектории помещаем в пространство меток Ck, которое содержит направляющие векторы движения каждого следующего бинарного образа относительно предыдущего [9].

1.3. Построение обобщенного бинарного эквивалента

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

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

Задача построения обобщенного бинарного эквивалента в нашем случае сводится к достраиванию объектов бинарных образов до получения их полной замкнутости.

Всякая область О плоскости бинарного образа серии Б^ содержит внутренние точки и точки контура (граничные точки). Первые из них обладают тем свойством, что не только они сами, но и их некоторая окрестность целиком принадлежит области О. Точки контура не являются внутренними, но в сколь угодно малой окрестности таких точек находятся внутренние точки области О и точки, не принадлежащие области О, — внешние (фоновые) точки. Область О обладает свойством связности, состоящим в том, что любые ее точки соединяются линией, целиком находящейся внутри О.

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

Таким образом, задача построения обобщенного бинарного эквивалента сводится к формированию контуров на образе серии Б^ [10]. Для формирования контура нет необходимости решать задачу обнаружения всех его точек с последующей проверкой обнаруженных точек на связность. С учетом сильной пространственной корреляции контурных точек и непрерывности линии контура, после обнаружения одной точки или небольшой группы точек резко сужается область пространства Ск, где располагается последующая контурная точка. В результате этого операция обнаружения начальной контурной точки и последующего продолжения контура сменяется операцией прослеживания контура, под этим понимается непрерывный переход от текущей точки контура к последующей до замыкания линии контура. Чтобы замкнуть линии контура, для любой точки каждого следующего бинарного образа строим операцию выбора области в пространстве измерений Ск, в которой следует извлекать свойства из любого вектора измерений в данной области. Для определения направления продолжения линий и конечной точки ее сходимости использовали установление эквивалентности совокупностей среди областей [11].

Процедура состоит из нескольких операций:

— установление эквивалентности ковариационных матриц выбранных областей пространства Ск;

— построение одной нелинейной функции, охватывающей кусочно-нелинейные функции;

— отыскание нелинейных соотношений для нескольких больших областей (можно потерять меньше информации, чем при ином способе);

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

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

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

1.4. Выделение объектов обобщенного бинарного эквивалента и вычисление их геометрических признаков

Для выделения объектов на обобщенном бинарном эквиваленте серии осуществляем поиск связных областей пикселов объектов и создаем матрицу меток Ь^, каждый элемент которой равен номеру объекта, которому принадлежит соответствующий пиксел изображения серии . Размер матрицы номеров объектов Ь^ равен размеру изображения серии . Объекты нумеруются по порядку, начиная с 1. Элементы, имеющие значение 1, относятся к первому объекту, имеющие значение 2 — относятся ко второму объекту и т.д. Если элемент в матрице Ь^ равен 0, то это означает, что соответствующий пиксел исходного изображения относится к фону (в данном случае соответствует трещинам на исходных изображениях рис. 2), т. е. не является объектом [13]. Для обхода контура используется критерий четырехсвязанности.

С целью дальнейшей визуализации отмеченных областей применили преобразование полученной матрицы меток Ь^ в палитровое изображение в градациях серого. Объект изображения серии окрашивается в определенный цвет в соответствии с номером объекта в матрице меток Ь^. Количество используемых цветов на единицу превышает количество объектов матрицы меток, фон соответствует красному цвету, для лучшей визуализации [14].

Под геометрическими признаками понимаем вычисление площади и размера объектов изображения серии .

На основании полученной матрицы меток вычисляем площадь каждого отмеченного в ней объекта в виде суммы

п т

А = ££ Ьц, (6)

г=1 j=1

где Ьц — значение бинарного изображения в точке, находящейся в г-й строке и ^'-м столбце. Здесь мы полагали, что изображение представляет собой матрицу с т столбцами и п строками. В данном случае под площадью понимаем сумму всех пикселов объекта в единицах площади элемента изображения. В результате получаем двумерный массив, содержащий значения площадей и соответствующее им количество объектов.

Под размером объекта (частицы) будем понимать длину малой оси эллипса, имеющего те же моменты инерции, что и плоская фигура соответствующего объекта [15-17].

Если принять, что N — количество пикселов, относящихся к объекту, а П — все множество пикселов объекта р(х,у), то координаты центра масс объекта определятся по формулам

х = N И х' ус = N Ц у■ (7)

р(х,у)£и р(х,у)£и

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

U = N Е (x - xc)2, U = N Е (У - Ус)2, (8)

p(x,y)en р(ж,у)еп

Uxy = nN E (x - Xc)(y - Ус), с ^(Ux - Uy)2 + 4Ux2y. (9)

p(x,y)en

Тогда размер в пикселах:

Amin = 2^2^Ux + Uy - C. (10)

В результате получаем массив, содержащий значения размера, площади и соответствующее им количество объектов [18].

2. Результаты моделирования разрушения мезообъема угля

Ранее было установлено [19, 20], что данные по тектонической структуре углей [21] дают основание в качестве изучаемого объема взять мезообъем с характерным линейным размером 500-600 мкм. В таком мезообъеме можно выделить и явно учесть в качестве элементов структуры отдельные микролитотипы (витринит, семивитринит, фюзинит, липтинит, минеральные примеси), а также систему пор размером 1-10 мкм. Пример подобного мезообъема показан на рис. 3. Влияние пор и трещин меньшего масштаба учитывается неявно, если в качестве определяющих соотношений для микролитотипов угля использовать упругопластическую модель среды с ненулевыми значениями коэффициентов внутреннего трения и дилатансии, например, модель Николаевского [22].

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

Рис. 3. Структура расчетного мезообъема, имеющего прямоугольное сечение размером 574 х 626 мкм

с разворотами отдельностей [21], то в качестве условия разрушения был использован комбинированный критерий. Проверялось выполнение двух условий:

1) достижение предельного значения неупругих деформаций (вне зависимости от напряжения);

2) достижение предельного значения растягивающего давления.

Предельные деформации варьировались в диапазоне е* = 0.01-0.13 % (е* ~ 0.02/а*, где а* — предел прочности на растяжение соответствующего компонента, мегапаскали), т. е. более прочные компоненты проявляют более хрупкое поведение. Фактически это эквивалентно заданию предельного значения работы напряжений на сдвиговых пластических деформациях. Предельные значения отрицательных давлений ограничивались величиной —0.5а*. Такой подход позволяет моделировать разрушение при сложном напряженном состоянии гетерогенной среды, когда в одних частицах преобладают сдвиги, в других — сжатие, а в третьих — растяжение.

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

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

Результаты проведенных расчетов представлены на рис. 4 и 5. Изображенная на рис. 4 картина разрушения в лагранжевом представлении (на исходной недеформиро-ванной сетке) определяется в первую очередь расположением пор, во вторую — геометрическими особенностями расположения структурных составляющих, а также способом нагружения. Приведенный результат получен на грубой сетке 574 х 626 ячеек, размер одной ячейки — 1 мкм, что не позволяет достаточно надежно выявить фрагменты разрушения размером меньше 10 мкм. Тем не менее некоторые выводы сделать можно.

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

Значения механических характеристик различных ингредиентов угля

Литотип Плот- Объемный Модуль Предел Коэффициент Коэффи-

ность, модуль сдвига, текучести внутреннего циент

г/см3 упругости, ГПа при растя- трения дила-

ГПа жении, МПа при нулевом давлении тансии

Витринит 1.25 33.33 7.14 2.30 0.5 0.17

Семивитринит 1.3 76.67 7.93 2.80 0.5 0.17

Липтинит 1.7 21.67 4.64 15.00 0.5 0.17

Фюзинит 1.4 383.3 39.66 135.00 0.5 0.17

Минеральные 2.5 33.82 29.74 19.00 0.5 0.17

примеси

Рис. 4. Схема нагружении в условиях одновременного действия сжатия и сдвига (а) и картины разрушения в изучаемом мезообъеме при использовании расчетной сетки 574 х 626 ячеек (б) и 1148 х 1252 ячеек (в). Общая деформация сдвига — 0.385 %, общая деформация сжатия — 0.77 %

Рис. 5. Структура другого мезообъема, имеющего прямоугольное сечение размером 574 х 626 мкм, (а) и полученная картина разрушения (б) (обозначения микролитотипов и схема нагру-жения те же, что и на рис. 3)

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

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

Получаемые результаты зависят от размера расчетной ячейки, и это накладывает определенные ограничения на параметры расчетной сетки для подобного рода расчетов. Поэтому для сравнения были проведены расчеты для того же мезообъема на расчетной сетке 1148 х 1252 ячеек (размер одной ячейки — 0.5 мкм).

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

Для проверки достоверности расчеты проводились для разных мезообъемов, состоящих из тех же микрокомпонентов, что и прежде, но с разным их процентным содержанием и пространственным расположением. Качественно результаты расчетов при этом остаются такими же. Пример расчета для другого мезообъема показан на рис. 5.

Кроме того, были проведены контрольные расчеты для мезообъемов меньшего размера (57.4 х 62.6 мкм2), контрольный образец был только из одного литотипа, но с явным учетом пор. Сетка использовалась с параметрами 574 х 626, т. е. шаг сетки равнялся 0.1 мкм и на квадрат размером 1 х 1 мкм приходится 100 ячеек. Таким образом, в этом случае также явно обнаруживаются куски размерами 1-50 мкм. Серия картин разрушения, полученная в одном из таких расчетов, представлена на рис. 2.

3. Результаты применения разработанной методики к анализу картин разрушения

Исходными данными для построения интегральной картины разрушения послужила серия (рис. 2) из ^ = 7 изображений (картин разрушения). Для создания бинарных образов использовали алгоритм пороговой бинаризации по заданным цветам. В данном случае для выделения областей разрушений (трещин) провели пороговую бинаризацию только для этого типа объектов (рис. 6, а—ж), для выделения объектов типа "пора" использовали аналогичный алгоритм бинаризации (рис. 6, з). То есть пикселу бинарного изображения присваивалось значение 1, если пиксел исходного изображения соответствовал заданному цвету, всем остальным пикселам изображения присваивалось значение 0.

Далее выполнялась процедура замыкания контуров бинарных образов серии изображений (рис. 6, а—ж), соответствующих разным последовательным моментам времени (процесса деформации), и было построено обобщенное изображение серии.

Эта процедура сводится к нескольким этапам.

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

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

3. Далее проводится выделение контуров. Считается, что контур должен иметь толщину не более 1 пиксела. Лишние точки контура удаляются, если толщина контура не превышает объект, поскольку контур может разделять два и более объектов, точки удаляются одинаково, с учетом наличия соседей. Если размер контура превышает

д е ж з

Рис. 6. Результат пороговой бинаризации картин разрушения (см. рис. 2) с идентификацией трещин (а—ж) и пор (з)

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

Следует отметить, что привлечение аппарата распознавания образов для построения обобщенного бинарного эквивалента повлекло за собой решение задачи об уменьшении размерности вектора измерений. То есть приведение совокупности измерений (с достаточно большим количеством данных) к совокупности признаков (с меньшим объемом данных). Решение задачи уменьшения размерности вектора измерений было вызвано необходимостью сокращения вычислений для идентификации контуров объектов бинарного изображения. В результате решения вопроса об уменьшении размерности [8, 11] перешли к следующему классу признаков для каждой точки изображения, из которых наиболее наглядными являются: наличие соседей (соседних точек) и пор, удаление ближайших соседей от точки контура, неудовлетворяющей критерию четырехсвязности, наличие направляющих векторов, направление роста контура (или увеличения совокупности точек) при вычислении отличительных признаков между двумя соседними изображениями. Пример построенных обобщенных изображений показан на рис. 7, а.

4. На следующем этапе находятся объекты, соответствующие порам, — они исключаются из рассмотрения (рис. 7, б). Затем на основании четырехсвязанной окрестности осуществляется поиск частиц разрушения на обобщенном бинарном изображении серии. В результате создается матрица меток, где каждый элемент равен номеру частицы, которому принадлежит соответствующий пиксел изображения.

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

После определения всех частиц на основании полученной матрицы меток вычисляются площади каждой частицы. В данном случае под площадью понимаем сумму всех пикселов объекта в метрическом эквиваленте (исходя из соответствия 1 пиксел = 0.01 мкм2). Для удобства сравнения полученных распределений объектов по площадям с экспериментальными данными округляем элементы массива до значений, соответствующих площадям квадратных ячеек сита с размером, равным целому числу микрометров, т. е. форматируем массив со значениями площадей выделенных объектов по следующему критерию: все элементы проверяются на соотношение г2 < А < {г + 1)2 и подсчитывается количество объектов, удовлетворяющих данному критерию, где А — площадь объекта, г = 1, 2,... (мкм). В результате получаем двумерный массив, содержащий значения площадей и соответствующее им количество частиц.

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

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

а

в

Рис. 7. Обобщенное изображение серии картин разрушения (а) с удалением из рассмотрения объектов типа "пора" (б) для визуализации отмеченных областей разрушений (в)

о4

|201 й 15^

т

10 50

0 5 10 15 20 25 30 Площадь частиц, мкм2

35

^ 30 «

я 251

Г 1 120

I 15

§ Ю и о о и

т

и

§ 51

0

0 1 2 3 4 5 Размер частиц, мкм

а

Рис. 8. Распределение выявленных частиц разрушения по площадям (а) и размерам (б)

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

пЩ

Е^

100.

БЕ

Здесь п^ — количество частиц размером Ог, — площадь к-й частицы размером Ог, Бе — площадь всех фрагментов разрушения. Следует отметить также следующий факт. Для диапазона от 1 до 10 мкм в эксперименте (седиментационный анализ) определяется дискретный набор размеров частиц, кратный 1 мкм. Далее в диапазоне от 10 до 100 мкм — с шагом 10 мкм. В расчетах же получаются частицы, размер которых содержат дробные доли микрометра. Поэтому при построении распределений частиц по размерам все фрагменты размером больше 1 мкм и меньшим или равным 2 мкм считаются соответствующими размеру 2 мкм. И так далее для других размеров.

Заключение

Отметим следующие особенности в применении предложенной методики.

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

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

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

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

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

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

[1] Дуда Р., Харт П. Распознавание образов и анализ сцен. М.: Мир, 1976. 513 с.

[2] Ищук И.Г., Трубицына Д.А. Особенности пылевыделения при различных нагрузках на комбайновые очистные забои (на примере шахт Кузбасса) // Рудничная аэрология и безопасность горного производства: Научные сообщения ННЦ ГП - ИГД им. А.А. Ско-чинского. 2005. Вып. 330. С. 47-58.

[3] Макаров П.В., Трубицын А.А., Трубицына Н.В. и др. Экспериментальное и теоретическое исследование разрушения углей и расчет выхода пылевых частиц. II. Численное изучение разрушения угля на мезо- и макроуровнях // Физ. мезомех. 2004. Т. 7. Спец. выпуск, ч. 2. С. 249-252.

[4] Макаров П.В., Трубицын А.А., Трубицына Н.В. и др. Численное изучение разрушения угля на мезо- и макроуровнях // Уголь. 2005. № 2. С. 33-36.

[5] Хорн Б.К.П. Зрение роботов: пер. с англ. М.: Мир, 1989. 487 с.

[6] Кириченко М.С., Панин С.В. Разработка адаптивного алгоритма оценки информативности динамических признаков для обработки и анализа изображений // Вычисл. технологии. 2005. Т. 10, № 1. С. 58-70.

[7] Претт У. Цифровая обработка изображений. М.: Мир, 1982. Кн. 2. 480 с.

[8] Анисимов Б.В., Курганов В.Д., Злобин В.К. Распознавание и цифровая обработка изображений: учеб. пособие для вузов. М.: Высшая школа, 1983. 295 с.

[9] Ту Дж., Гонсалес Р. Принципы распознавания образов: пер. с англ. М.: Мир, 1978. 403 с.

[10] ПАвлидис Т. Алгоритмы машинной графики и обработки изображений. М.: Радио и связь, 1986. 400 с.

[11] Патрик Э.А. Основы теории распознавания образов. М.: Советское радио, 1980. 407 с.

[12] Введение в контурный анализ; приложения к обработке изображений и сигналов / Я.А. Фурман, А.В. Кревецкий, А.К. Передреев, А.А. Роженцов, Р.Г. Хафизов, И.Л. Его-шина, А.Н. Леухин; под ред. Я.А. Фурмана. М.: Физматлит, 2003. 592 с.

[13] Писаревский А.Н. и др. Системы технического зрения: (Принципиальные основы, аппаратное и математическое обеспечение). Л.: Машиностроение, 1988. 424 с.

[14] Haraliok R.M., Shapiro L.G. Computer and Robot Vision. Reading: Addison-Wesley, 1992. Vol. 1. 1302 p.

[15] Горшков А.Г., Трошин В.Н., Шалашилин В.И. Сопротивление материалов. М.: Физматлит, 2005. 544 с.

[16] САргсян А.Е. Сопротивление материалов, теории упругости и пластичности. Основы теории с примерами расчетов. М.: Высшая школа, 2000. 286 с.

[17] ЯковЕнко Г.Н. Краткий курс теоретической механики. М.: БИНОМ. Лаборатория знаний, 2005. 125 с.

[18] Методы компьютерной обработки изображений / Под ред. В.А. Сойфера. М.: Физмат-лит, 2001. 784 с.

[19] Макаров П.В., Смолин И.Ю., Черепанов О.И. и др. Упруговязкопластическая деформация и разрушение угля на мезоскопическом масштабном уровне // Физ. мезомех. 2002. Т. 5, № 3. С. 63-87.

[20] Адаптация методов мезомеханики к исследованию процессов деформации и разрушения угля / А.А. Трубицын, П.В. Макаров, О.И. Черепанов, С.П. Ворошилов, Н.В. Трубицына, И.Ю. Смолин, В.В. Соболев, Я.С. Ворошилов, В.В. Киселев, С. Грюнинг. Кемерово: Кузбасс-ЦОТ, 2002. 116 с.

[21] Саранчук В.И., Айруни А.Т., Ковалев К.Е. Надмолекулярная организация, структура и свойства угля. Киев: Наукова думка, 1988. 192 с.

[22] Makaroy P.V., Stefanoy Y.P., Smolin I.Yu., Cherepanoy O.I. Modeling of mechanical behavior of geomaterials on the mesoscale // Intern. J. for Multiscale Computational Engineering. 2005. Vol. 3, Issue 2. P. 135-148.

Поступила в редакцию 1 октября 2007 г.