Программное обеспечение для оптимизации сетки БВР в условиях сложной анизотропии трещиноватости

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

Для решения задачи учёта сложной анизотропии трещиноватости (3 системы трещин и более) в ОАО «ВИОГЕМ» разработаны математические модели распределения энергии взрыва в трещиноватой среде [1], методы визуализации трещиноватости, модели сеток скважин, методы и алгоритмы измерения кусковатости по снимкам лазерного сканера и методы сопоставления распределений размеров кусков и величин энергии для построения зависимости между ними. Эти математические модели и методы реализованы в виде автоматизированной системы научных исследований (АСНИ) «Облако 3D» и программного комплекса, состоящего из следующих модулей:

1. визуализация трещиноватости;
2. построение гистограмм распределений;
3. сопоставление распределений (инструмент построения
   функциональной зависимости по распределениям величин);
4. анизотропная оптимизация БВР в условиях трещиноватости;
5. оконтуривание фотоснимков развала;
6. анализ контуров кусков;
7. изотропная оптимизация БВР.

Интенсивность трещиноватости L в направлении единичного вектора определяется [2] как:

\(L(\vec{e})=\sum_{i=1}^{n} |np_{e} \cdot \vec{w_{i}}| = \sum_{i=1}^{n} |\vec{w_{i}} \cdot \vec{e}|\),

где n – число трещин, – \(\vec{w}\) вектор i-той системы трещин (перпендикулярный трещинам в системе, длина вектора равна числу трещин на единицу длины [2]).

Среднее расстояние до ближайшей трещины в направлении \(\vec{e}\) определяется как:

\(d(\vec{e})=1/L(\vec{e})\)

Для двумерной и трёхмерной визуализации интенсивности трещиноватости \(L(\vec{e})\) и диаграммы среднего расстояния \(d(\vec{e})\) до ближайшей трещины разработан модуль «Анализ трещиноватости» (рис. 1).

а)
б)

Рисунок 1 – Снимки экрана модуля анализа трещиноватости.
a) 2D б) 3D

Оптимизация сетки скважин осуществляется в модуле «Анизотропная оптимизация БВР» (рис. 2).

Рисунок 2 – Снимки экрана модуля анизотропной оптимизации БВР

При выполнении оптимизации буровзрывных работ максимизируется выходная доля породы в кусках размером в заданном диапазоне (dmin; dmax). Условие максимальной доли выхода кусков требуемых размеров в объёме V буровзрывного блока можно записать как:

\(\int_{v} C_{d}(d_{cp}(dV))dV \rightarrow min\) (1)

где \(dcp(dV)\) – размер среднего выходного куска в объёме \(dV\) блока; \(Cd(x)\) –штрафная функция. График примера такой функции приведён на рис. 3а.

a)

б)

Рисунок 3 – Графики штрафных функций от размера куска; б) от номинальной энергии

Если кусок имеет требуемый размер (dmin ≤ x ≤ dmax), штрафная функция принимает меньшие значения. Условие (1) можно выразить через энергию, переданную объёму куска dV:

\(\int_{v} C_{d}(E(dV))dV \rightarrow min\)

где E(dV) – энергия, переданная объёму dV, CE(E(dV)) – штрафная функция (аналогичная Cd(x), см. рис. 3б), принимающая свои наименьшие значения в диапазоне (Emin,Emax).

Штрафная функция может строиться на основе двух сигмоид:

\(CE(E)=\frac{1}{1+e^{-(E-E_{min})/k_{min}}}+\frac{1}{1+e^{-(E-E_{max})/k_{max}}}\) (2)

kmin < 0; kmax > 0,где E – значение номинальной энергии (Дж/м3); Emin – нижний порог допустимой номинальной энергии; Emax – верхний порог допустимой номинальной энергии; kmin и kmax – параметры растяжения сигмоид, определяющие ширины переходов от 0 к 1. Для определения порогов Emin и Emax предлагается сопоставлять [1] замеренные распределения кусков по размерам и рассчитанные распределения участков блока по уровням номинальной энергии, для чего разработан модуль «Сопоставление распределений» (рис. 4).

Рисунок 4 – Снимок экрана модуля построения функциональной зависимости величин по их распределениям

Согласно разработанным математическим моделям распределения и затухания энергии в массиве, энергия взрывов от скважин в точке P записывается как:

\(L(\vec{e})=\sum_{i=1}^{n}|np_{e}\cdot \vec{w_{i}}|=\sum_{i=1}^{n}|\vec{w_{i}}\cdot \vec{e}|\)

\(f(\vec{r_{ij}})=|\vec{r_{ij}}|L(\vec{r_{ij}})\)

где E0ij – начальная энергия взрыва в ij-той скважине, – \(\vec{r_{ij}}\)вектор направления от точки P до ij-той скважины, re – расстояние, на котором в данной среде энергия взрыва убывает в e раз, \(f(\vec{r_{ij}})\) – число трещин, пройденных волной от скважины до точки P, – \(L(\vec{r_{ij}})\) интенсивность трещиноватости в направлении \(\vec{r_{ij}}\); k – коэффициент ослабевания энергии волны при прохождении трещины (соотношение энергии после прохождения трещины к энергии до
прохождения, 0<k<1). Фактор слоистости породы можно учесть, сделав расстояние re функцией направления на скважину.

При моделировании двумерных сеток буровзрывных скважин вводится понятие «порождающая матрица регулярной сетки скважин»[1]. Координаты (хs;ys)ij ij-той скважины предлагается вычислять умножением индексов скважины ij на матрицу 2х2:

\(\binom{x_{s}}{y_{s}}=M_{c}\binom{i}{j}=\binom{m_{11}m_{12}}{m_{12}m_{22}}\binom{i}{j}=\binom{m_{11^{i}} \underline{m_{11^{j}}}}{m_{21^{i}}{m_{22^{j}}}}\)

где MC – вещественная матрица размером 2×2; i, j – целочисленные индексы скважины. Для ввода порождающей матрицы сетки скважин используется разработанный предметно-ориентированный язык [3] (англ. domain-specific language, DSL), включающий операции поворота, масштабирования, ввода общепринятых сеток скважин (шахматной, прямоугольной и т.д.).

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

\(M_{,}^{c}=\binom{cos\beta -sin\beta}{sin\beta\ cos\beta}\binom{a c}{0 b}\) (3)

где β – угол поворота, градусы; a – расстояние между скважинами в ряду (РМС), м; b – расстояние между рядами (РМР), м; с – (абсолютный) сдвиг соседнего ряда, м. Значения параметров a, b, c, β наглядны, и их можно рассчитать для любой порождающей матрицы M.

Рисунок 5 – Шахматная сетка скважин с поворотом и сдвигом

Для удобства расчётов вводятся параметры сетки скважин \(d=\sqrt{b/a}\)\(q=^{\sqrt{ab}},h=c/a.g\) (м) равен расстоянию между рядами и между скважинами при квадратной сетке, дающей равный удельный расход с шахматной сеткой (a,b);
d описывает соотношение сторон шахматной сетки; h равен соотношению параметров c и a.

Параметр g вычисляется из требуемого удельного расхода и вместимости скважин и постоянен при минимизации целевой функции, параметры d,h,β могут меняться. Целевая функция записывается так:

\(F(d,h,\beta)=\int_{0}^{q/d}\int_{0}^{qd}C_{E}\left ( \sum_{ij}E_{0ij} exp \left [ -\frac{\vec{|r_{ij}}|}{r_{e}} \right ] k^{f\ast (\vec{r}_{ij})}\right )dydx\rightarrow min\)

\(\vec{r}_{ij} = \left \{ x_{ij} – x; y_{ij}-y\right \}\)

\(\binom{x_{s}}{y_{s}}=M_{cмод}\binom{i}{j}\)

\(M _{c }=q\begin{pmatrix}
1/d & h/d \\
0 & d
\end{pmatrix}\)

где {xsij;ysij} – координаты ij-той скважины; CE– штрафная функция (2) номинальной энергии; Мс мод.– порождающая матрица сетки на основе матрицы (3), но без поворота; i – номер ряда (i=-m,-m+1,…,0,…,m-1,m), j – номер скважины в ряду (k=-n,-n+1,…0,…,n-1,n), a,b,c – параметры сетки скважин; – число трещин в направлении вектора , повёрнутого на угол (-β). Размеры сетки m,n берутся минимально такими, чтобы влияние крайних скважин на прямоугольник (0,0)-(a,b) было незначительным.

Необходимые для точного измерения кусковатости алгоритмы автоматического выделения и измерения кусков реализованы в АСНИ «Облако 3D» [4] (рис.6), предназначенной для обработки облаков точек.

Рисунок 6 – Снимок экрана АСНИ «Облако 3D»

Для измерения кусковатости в развале предлагается делать снимок развала лазерным 3D сканером (облако точек Q={(xi, yi, zi)}, i=1,2,…,n, где xi, yi, zi – координаты точек в декартовой системе координат, n – число точек в облаке) b определять размеры кусков на основе математических моделей кусков.

Из облака точек Q вручную или автоматически выделяются точки, принадлежащие куску, – эти точки образуют «модель облака точек» Lpc={={(xi, yi, zi)}, i=1,2,…n, , где (xi, yi, zi) – координаты i-той точки, принадлежащей поверхности куска, в трёхмерном пространстве; n – число точек в модели. Для измерения размеров наиболее пригодны сферическая модель куска Lq={d}, где d – линейный размер куска, и эллиптическая модель Le={a; b; c}, где a, b, c – линейные размеры куска вдоль трёх осей.

Размер куска в эллиптической модели может рассчитываться несколькими способами: как максимальный размер эллипсоида (или максимальный размер стороны параллелепипеда): d=max(a,b,c); как максимальная диагональ параллелепипеда \(d=\sqrt{a^2+b^2+c^2}\) ; как среднее арифметическое размеров вдоль осей:\( d=a+b+c/3\); как среднее геометрическое размеров вдоль осей (размер шара или куба, эквивалентного по объёму данному эллипсоиду или параллелепипеду): \(d=^{\sqrt[3]{abc}}\) .

В эллиптической модели можно определить меру близости куска к шарообразной, плоской или игольчатой форме, например в виде числа m в диапазоне [0;1] (1 – кусок близок к выбранной форме, 0 – не близок). Для этого упорядочим размеры a,b,c по неубыванию и назовём их d1, d2, d3 ({d1,d2,d3}={a,b,c}, d1 ≤ d2 ≤ d3) . Мера близости к форме шара: ш=d1/d3 ; к плоской форме: п=1-d1/d2 ; к игольчатой форме: и=1-d2/d3. Лещадность куска можно определить, вычислив соотношение l максимальных и минимальных его размеров: l=d3/d1. Если l > 3, кусок считается лещадным, и распределение, имеющее высокую долю таких кусков, будет лешадным.

Для аппроксимации облака точек куска L сферой минимизируется сумма квадратов расстояний от точек облака L до =

= \(\sum_{\forall PeL}\left ( r-\sqrt{(x_{i}-c_{x})^2+(y_{i}-c_{y})^2+(z_{i}-c_{z})^2} \right )^{2}\rightarrow min\) (5)

где C=(cx;cy;cz)– координаты центра сферы; P=(xi,yi,zi)– точки из облака точек L куска; r – радиус сферы. Для минимизации применяется метод Левенберга-Марквардта [5], использующий частные производные функции (5) по параметрам cx,cy,cz,r. Параметр d сферической модели (размер куска) определяется по формуле:

\(d=1^{|r|}\).

Для аппроксимации облака точек куска L эллипсоидом минимизируется сумма квадратов расстояний от точек облака L до поверхности эллипсоида в направлении его центра:

\(F(C,M)=\sum_{\forall PeL}(\frac{|P-C|}{\sqrt{(P-C)^{T}M(P-C)}} – |P-C|)^{2}\rightarrow min\) (6)

где C=(cx;cy;cz) – координаты центра эллипсоида; P=(xi,yi,zi) – точки из облака точек L куска;

\(M=\begin{pmatrix}
m_{11} & m_{12} & m_{13}\\
m_{21} & m_{22} & m_{23}\\
m_{31} & m_{32} & m_{33}
\end{pmatrix}\) – матрица квадратичной формы эллипсоида ;

(PE-C)TM(PE-C)=1; PE=(x;y;z) – точка на поверхности эллипсоида. Для минимизации применяется метод Левенберга-Марквардта [5], использующий частные производные функции (6) по параметрам cx,cy,cz,m11,m22,m33,m12,m23,m13. Размеры осей эллипсоида a,b,c (или d1≤d2≤d3), a>0,b>0,c>0 рассчитываются из собственных чисел \(\lambda 1\geq \lambda 2\geq \lambda 3\) матрицы M квадратичной формы:

\(a=d_{1}\frac{2}{\sqrt{\lambda _{1}}};a=d_{2}\frac{2}{\sqrt{\lambda _{2}}};a=d_{3}\frac{2}{\sqrt{\lambda _{3}}};\)

Разработанный алгоритм выделения кусков на снимках лазерного 3Dсканера основан на многократном разномасштабном применении нейронной сети обнаружения границ кусков (обозначим её «EdgeNet», см. рис 7в) к участкам масштабированных карт глубин, выделении кусков заливкой замкнутых областей в результатах EdgeNet, фильтрации недостоверно выделенных кусков и уточнении границ кусков с помощью нейронной сети маски куска (обозначим её «MaskNet», см.рис. 7г).

Этот подход базируется на идее архитектуры нейронной сети Mask R-CNN [6], используемой для семантической сегментации и содержащей сеть Faster R-CNN [7] для определения расположения кусков и свёрточную сеть-автокодировщик [8,9] для определения маски куска, но сеть Faster R-CNN из-за неуспешного обучения на данной задаче заменена сетью определения границ кусков EdgeNet и алгоритмом комбинирования результатов EdgeNet.

Разработанные свёрточные сети EdgeNet и MaskNet имеют архитектуру автокодировщиков, 8 слоёв, ядро свёртки 5х5, шаг свёртки 2х2, сигмоидальную функцию активации на последнем слое и ReLU на промежуточных, принимают на вход нормализованные в диапазон [0;1] карты глубин.

Входные и выходные размеры сети MaskNet 93х93х1, EdgeNet – 253х253х1 на входе и 253х253х3 на выходе (3 канала «ребро/кусок/не кусок»). Сеть EdgeNet достоверно выделяет границы внутри «активной области» размером примерно 160х160px в центре, поэтому алгоритм перебирает различные варианты области её применения и выделяет куски, одинаково определённые на разных областях (рис.7д).

Рисунок 7 – Исходные данные и результаты алгоритма выделения кусков
а) точки снимка в естественных цветах, чёрный фон; б) карта глубин участка (чёрный – ближе, белый – дальше, сплошная рамка – граница области рис.7в, внутренняя пунктирная рамка – граница области рис.7г); в) результат работы EdgeNet (красный канал – «ребро», зелёный канал – «кусок», синий канал – «не кусок», рамкой обозначена граница активной области); г) результат работы MaskNet (синий – 0.0 / не кусок, красный – 1.0 / кусок); д) участок с выделенными кусками (ортопроекция на экран, точки
кусков отмечены жёлтым цветом, синяя рамка указывает расположение рис.7а,б)

Разработанный алгоритм позволяет корректно выделить порядка 80% кусков. На основании изложенных методов предлагается следующая методика оптимизации сетки буровзрывных скважин:

• Определение коэффициентов затухания энергии re, k для данного типа пород (может не выполняться, если эти коэффициенты уже известны);
• Определение параметров штрафной функции CE (пороговых значений номинальной энергии Emin, Emax и параметров ширины порогов kmin, kmax).
• Определение требуемого удельного расхода ВВ и расчёт количества ВВ в скважинах, энергий E0i, выделенных в скважинах и параметров равномерной сетки a,b (и g);
• Определение экстремумов интенсивности трещиноватости для определения формы зоны разрушения скважины;
• Подбор различных наборов параметров d, h, β с целью минимизации функции (4) с использования численного способа оптимизации, построение и визуальная оценка сетки скважин, выбор приемлемого набора параметров;
• Расчёт новых рекомендуемых значений a, b, c (РМС, РМР, сдвиг ряда соответственно) по параметрам q, d, h для сетки скважин (3) (см. рис. 5).

Литература:

1. Коновалов А. В. Учёт анизотропии трещиноватости пород в математических моделях расстановки буровзрывных скважин /
   А. В. Коновалов, Г. М. Редькин // Научные ведомости Белгородского государственного университета. Экономика. Информатика. Том 45, №3, сентябрь 2018 – С.523-536.
2. Редькин, Г. М. Нестационарное анизотропное математическое моделирование неоднородностей систем минерального сырья /
   Г. М. Редькин. – М.: Издательство Ассоциации строительных вузов, 2007. – 500 с.
3. Предметно-ориентированный язык / Википедия, свободная энциклопедия // Электронный ресурс URL: https://ru.wikipedia.org/wiki/Предметно-ориентированный_язык
4. Свидетельство о государственной регистрации программы для ЭВМ № 2016610687 «Автоматизированная система научных исследований Облако 3D (АСНИ Облако 3D)». Коновалов А. В., Игнатенко И. М., Агарков И.Б. Заявка № 2015662094, дата поступления 30 ноября 2015 г. Зарегистрировано в Реестре программ для ЭВМ 18 января 2016 г.
5. Marquardt D. An Algorithm for Least-Squares Estimation of Nonlinear Parameters / SIAM Journal on Applied Mathematics, 1963, 11 (2): p.431–441.
6. He K. Mask R-CNN / K. He, G. Gkioxari, P. Dollar, R. Girshick // International Conference on Computer Vision (ICCV), October 22, 2017. URL:
   https://arxiv.org/abs/1703.06870v3
7. Ren Sh. Faster R-CNN: Towards Real-Time Object Detection with Region Proposal Networks / Sh. Ren, K. He, R. Girshick, J. Sun // URL:
   https://arxiv.org/abs/1506.01497v3
8. Николенко С. Глубокое обучение. / С. Николенко, А. Кадурин, Е. Архангельская // СПб.: Питер, 2018. – 480 с.: ил.
9. Goodfellow I. Deep Learning / Ian Goodfellow, Yoshua Bengio, Aaron Courville // MIT Press, 2016. – 800 p.