Горный массив как нелинейная динамическая система. Математическое моделирование эволюции напряженно-деформированного состояния горного массива в окрестностях выработки

УДК 539.3, 539.4

Горный массив как нелинейная динамическая система. Математическое моделирование эволюции напряженно-деформированного состояния горного массива в окрестностях выработки

П.В. Макаров1'2, М.О. Еремин12

1 Национальный исследовательский Томский государственный университет, Томск, 634050, Россия 2 Институт физики прочности и материаловедения СО РАН, Томск, 634055, Россия

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

Ключевые слова: математическое моделирование, горный массив с выработками, нелинейные динамические системы, накопление повреждений, режим с обострением, самоорганизованная критичность

Rock mass as a nonlinear dynamic system. Mathematical modeling of stress-strain state evolution in the rock mass around a mine opening

P.V. Makarov12 and M.O. Eremin1-2

1 National Research Tomsk State University, Tomsk, 634050, Russia 2 Institute of Strength Physics and Materials Science SB RAS, Tomsk, 634055, Russia

The paper briefly reviews the fundamental (general) evolution properties of nonlinear dynamic systems. The stress-strain state evolution in a rock mass with mine openings has been numerically modeled, including the catastrophic stage of roof failure. The results of modeling the catastrophic failure of rock mass elements are analyzed in the framework of the theory of nonlinear dynamic systems. Solutions of solid mechanics equations are shown to exhibit all characteristic features of nonlinear dynamic system evolution, such as dynamic chaos, self-organized criticality, and catastrophic superfast stress-strain state evolution on the final stage of failure. The calculated seismic events comply with the Gutenberg-Richter law. The cut-off effect has been obtained in numerical computation (downward bending of the recurrence curve in the region of large-scale failure events). Prior to catastrophic failure, the diagrams of the probability density function of stress fluctuations deviate from the average trend, the slope of the recurrence curve of calculated seismic events becomes more gentle, seismic quiescence regions form in the central zones of the roof, and more active deformation begins at the periphery of the opening. These factors point to the increasing probability of a catastrophic event and can be considered as catastrophic failure precursors.

Keywords: mathematical modeling, rock mass with opening, nonlinear dynamic systems, damage accumulation, blow-up mode, self-organized criticality

© Макаров П.В., Еремин М.О., 2016

1. Введение. Горный массив с выработками — типичная нелинейная динамическая система. Возможен ли прогноз катастрофических разрушений?

1.1. Базовые свойства и имитационные модели нелинейных динамических систем

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

В качестве базовой модели эволюции геосреды обычно используется модель «кучи песка» (sand pile, SP), которую мы обсудим коротко ниже. В основе прогноза опасных динамических проявлений горного давления, включая горные удары, лежат данные, полученные при ведении горных работ и геофизическом мониторинге. Ведутся активные исследования по разработке прогностических теорий, основанных в том числе на фундаментальных законах эволюции нелинейных динамических систем [3-6]. Идея масштабной инвариантности и самоподобия разрушения на разных масштабах позволяет привлекать для анализа катастрофического разрушения горного массива с выработками обширные данные по землетрясениям, а также данные лабораторных экспериментов по разрушению малых образцов. Все больше появляется работ, посвященных разработке фундаментальной теории сейсмического процесса, независимо от масштаба событий. В большинстве таких работ в основу разрабатываемой теории положены хорошо изученные базовые уравнения нелинейной динамики. Однако многие детали формирования очага катастрофического разрушения, так же как и сейсмического процесса, недостаточно изучены. Не существует общепризнанной физической модели процессов разрушения. По этим причинам широкое распространение получили имитационные модели, которые демонстрируют некоторые свойства реальных динамических систем, например режимы с обострением, степенные распределения

63

и другие важные свойства эволюции динамических систем [2]. Наиболее популярными являются модели тепловых структур, в основе которых лежит нелинейное уравнение теплопроводности [7-9], а также различные модификации модели «кучи песка», первый вариант которой был предложен в работе [10].

Интерес к модели тепловых структур как базовой модели теории нелинейных динамических систем продиктован следующими фундаментальными свойствами решений нелинейного уравнения теплопроводности: 1) возможность развития самоорганизации — возникновения тепловых структур; 2) наличие режимов с обострением, т.е. катастроф. Режимы с обострением впервые были изучены при исследовании свойств решений нелинейного уравнения теплопроводности [9].

С нашей точки зрения, привлечение нелинейного уравнения теплопроводности для имитации и объяснения процессов автокаталитического катастрофического разрушения слабо аргументировано. Это уравнение имеет следующий вид, допускающий анализ общих свойств решений:

Т =

— < х <7, Т(х, 0) = Т0(х),

где к (Т) = к0Тс — нелинейная функция температуропроводности; Q(T) = q0Tв — нелинейный источник тепла, причем Р = ст+1 и Р>1.

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

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

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

Заметим также, что уравнение вида

¿и 7 ¿и

-= Q

& о ^и)

также приводит к режимам с обострением.

Неравенство в (2) означает конечность времени обострения и является необходимым и достаточным условием существования режима с обострением [7].

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

64

тах эволюции напряженно-деформированного состояния и разрушения в геосредах для задания скорости накопления в среде повреждений [11, 12].

Более аргументированным подходом для прогноза сильных землетрясений являются различные модификации модели «кучи песка». Разработанные на их основе количественные алгоритмы, построенные на основе различных сценариев землетрясений, например акти-визационного и антиактивизационного, адаптируются к прогнозу сильных землетрясений, а также горных ударов. Например, широко известен алгоритм М8 [1, 13]. Вариации моделей «куча песка» обладают свойством самоорганизованной критичности, в результате их эволюции получаются степенные распределения количества лавин N по их размерам (энергиям) Е [2, 13-15]:

N

где показатель степени Ь при соответствующем выборе параметров близок к единице, что отвечает фундаментальному закону повторяемости Гутенберга-Рихтера. Степенные распределения (3) универсальны, являются характерным признаком эволюции множества динамических систем, демонстрирующих сложное поведение, и выполняются практически без исключений при описании различного рода катастроф и стихийных бедствий [16-18]. Как увидим далее, степенным законам распределения соответствуют все многомасштабные разрушения элементов горного массива с выработками. Расчетные сейсмические события, полученные в результате численного решения системы нелинейных уравнений механики деформируемого твердого тела, также приводят к типичному графику повторяемости сейсмических событий (закону Гутенберга-Рихтера), который отражает степенное распределение сейсмических событий разных масштабов. Этот факт означает, во-первых, что эволюция напряженно-деформированного состояния в горном массиве с выработками соответствует типичным сценариям эволюции нелинейных динамических систем, а во-вторых, убедительно доказывает, что решения уравнений механики деформируемого твердого тела демонстрируют характерные для динамических систем сценарии, отвечающие степенным распределениям (в данном случае масштабов разрушения). Физическое содержание степенных законов распределения, как известно, заключается во взаимной зависимости происходящих в системе событий, что невозможно без информационного обмена между элементами динамической системы. Именно информационный обмен между различными частями системы обеспечивает ее самосогласованный когерентный отклик на воздействия.

Самоорганизованная критичность нелинейных динамических систем рассматривается как универсальный механизм катастроф [7, 10, 16], в такой системе катастрофы неизбежны и являются внутренним свойст-

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

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

Самоподобие, масштабная инвариантность — фундаментальные свойства эволюции динамических систем, обладающих свойством самоорганизованной критичности. Эти свойства в той или иной мере демонстрирует базовая имитационная модель «куча песка», а также различные ее модификации. Однако существует мнение, что системы, обладающие свойством самоорганизованной критичности, непредсказуемы [7, 16, 19]. Этот вопрос имеет принципиальное значение. Впервые он был поднят в дискуссии в журнале Nature [20] в связи с моделью ««куча песка» Бака-Тонга-Визенфельда [10] и развитой на ее основе теорией самоорганизованной критичности. Фактически дискуссия велась (и периодически возобновляется) по поводу возможности прогноза крупных событий в модели Бака-Тонга-Визен-фельда и ее модификациях. Эта дискуссия была перенесена на возможность прогноза катастрофических событий, прежде всего землетрясений, в реальных природных и других системах, для которых характерны степенные законы распределения и которые, следовательно, обладают свойством масштабной инвариантности и самоорганизованной критичности.

В настоящее время в различных вариантах моделей «куча песка», по мнению ряда авторов, предсказуемость существенно повысилась, так как модифицированные методики основываются на предвестниках, которые адаптированы к разработанной модельной динамике. В этих моделях используются как активизационные, так и антиактивизационные сценарии [2, 19, 21].

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

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

Рис. 1. Степенная статистика катастроф и бедствий по данным [23]; ранжировка техногенных катастроф по числу погибших (2047 крупнейших событий) (а); ранжировка стихийных бедствий по числу раненых (1084 события) (б). х(г) — размерность, г — ранг, упорядоченный по убыванию х

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

1.2. Проблема ««сМ-off»

При анализе эмпирических данных и изучении различных имитационных моделей самоорганизующихся систем исследователи встретились с так называемой проблемой «cut-off» или обрезанием энергии, выражающейся в загибе графика повторяемости вниз в области сильнейших сейсмических событий [2, 7, 22]. Этот эффект означает уменьшение вероятности крупных катастрофических событий. Особенности поведения хвоста зависимости распределения сейсмических событий в области редких сильнейших землетрясений имеют первостепенное значение для их прогноза [22]. Однако эта проблема не решается на основе наблюдений из-за крайней ограниченности эмпирических данных. По этой причине численные эксперименты, выполненные на основе различных моделей, приобретают особое значение.

При теоретическом анализе различных имитационных моделей эффект cut-off может быть связан с тем, что при ограниченных размерах расчетных областей система достигает масштабного предела. Отметим, что степенные распределения для многих реальных нелинейных процессов также демонстрируют эффект cutoff. На рис. 1 приведены распределения по статистике катастроф и бедствий [23], демонстрирующие этот эффект. К подобным зависимостям приводит, например, статистика заражения компьютеров вирусами, распространение эпидемий и т.д. [7].

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

лей. Для геосред эффект cut-off может быть связан с особенностями реологии и разрушения геосреды, а также с ограничениями, накладываемыми прочностными характеристиками и структурной организацией геосреды, на предельный масштаб возможный катастрофы. Вероятно, что определяющую роль играют эффекты дальнодействия концентраторов напряжений. Это явление также зависит от линейного размера L изучаемой системы: cut-off ~ L"22, однако такая зависимость больше имеет модельный характер и, скорее всего, ее нельзя распространить на системы с очень большими L [24, 25]. В численных расчетах процесса разрушения, представленных в настоящей работе и выполненных на основе уравнений механики деформируемого твердого тела, график повторяемости расчетных сейсмических событий также демонстрирует этот эффект.

1.3. Фликкер-шум

Фундаментальным свойством всех нелинейных динамических систем (реальных и искусственных) является генерация в них фликкер-шума 1/f. Спектр мощности фликкер-шума на низких частотах отвечает степенной зависимости

S

где в ~ 1.

Физическая природа 1/f шума считается невыясненной, несмотря на обширные исследования этого явления и огромное число публикаций по этой теме. Например, в электрических цепях 1/f шум связывают с наличием примесей. Интенсивный широкополосный шум наблюдается в областях критического неравновесного фазового перехода. Степенной характер спектра мощности означает отсутствие в системе характерных частот, а значит, и характерных времен. Он также показывает, что большая часть энергии отвечает медленным процессам, соответствующим флуктуациям очень большой интенсивности, т.е. крупномасштабным катастрофам. Также все события в системах, для которых спектр мощности имеет степенной характер, не отвечают каким-либо требованиям повторяемости и периодичнос-

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

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

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

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

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

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

1.4. Свойство целостности динамических систем

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

В классических моделях «куча песка», блоковых и других [10, 16, 17] полагается, что целостность системы достигается за счет самоорганизации в критическое состояние. Механизмом, обеспечивающим взаимодействия во всей системе в целом, является наличие двух противоположных процессов — естественная эволюция системы (например увеличение локальных наклонов кучи песка за счет добавления песчинок) и процесс отбраковки (лавины осыпания). Необходимо, чтобы скорость отбраковки была много больше скорости естественной эволюции. Фактического прямого информационного обмена в таких моделях нет. На наш взгляд, это существенный недостаток большинства имитационных базовых моделей динамических катастрофических событий.

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

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

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

В представленной численной модели введен параметр сжатия расчетного времени К, что позволяет легко пересчитать условное расчетное время на реальное время развития деформационного процесса:

*геа1 = К^са1е.

Из соображений устойчивости численных расчетов временные шаги очень малы, общее время расчетного деформационного процесса не превышает нескольких секунд или минут на современных ЭВМ при общем количестве временных шагов ~ п -106 (конечно, имеется в виду не затраченное машинное время, которое может составлять дни при 3D-моделировании). Параметр К выбирается из соображений обеспечения целостности системы — за каждый шаг нагружения (например приращения забойного пространства) сформировавшаяся волна напряжений должна пройти по расчетной области 1-3 раза, обеспечивая в среде информационный обмен и соответствующую подстройку напряженно-деформированного состояния во всей расчетной области так, как это происходит в реальности.

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

1.5. О геофизическом мониторинге в горных массивах с выработками

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

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

В горных массивах с выработками контроль геодинамической ситуации ведется на основе методов геодинамического мониторинга, важнейшими из которых являются методы регистрации микросейсмических событий и геоакустический метод [26].

В связи с отсутствием в действующих нормативных документах методик учета времени протекания геомеханических процессов и оценок нестационарного характера роста поврежденности в горных массивах (в частности в кровле) вблизи выработок, большое значение приобрели методы оперативного прогноза опасных проявлений горного давления. Например, методы сейсмической оценки, основанные на системах мониторинга, активно используются на шахтах России, Германии, Японии, Австралии, США, ЮАР, Китая [18, 27, 28]. Решения множества комплексных проблем по управлению геомеханическими процессами, возникающих в современных условиях больших нагрузок на шахтные забои, во многих случаях опираются также на методы физического моделирования геомеханической ситуации [26, 29], в том числе горных ударов [30], что крайне важно, но также ограничено отсутствием общей методологии, в частности количественной оценки предлагаемых критериев. В горной инженерии ведутся активные работы по учету поправок на нестационарность процесса нагружения горного массива с целью учета динамики геомеханических процессов, оценки угроз негативных динамических явлений и выработки мер и технологических решений по их предотвращению. Основные усилия сосредоточены на экспериментально-теоретических методах, методах физического моделирования и развитии методов современного геодинамического мониторинга.

Большинство систем сейсмического мониторинга хотя и обеспечивают непрерывность сейсмологических наблюдений за развитием макромасштабных и микросейсмических процессов на территориях горного отвода, позволяют отслеживать горнотектонические удары, посадки кровли и сопутствующие формы динамических разрушений [31-34], но нет обоснованных методик оценки близости горного массива к критическому состоянию на основе полученных данных. Так, система сейсмического мониторинга GITS, разработанная ВНИМИ и внедренная на ряде шахт Кузбасса, осуществляет контроль в пределах шахтного поля (5x5 км2) с энергетическим диапазоном событий от 100 Дж с интенсивностью потока событий до 100 событий в сутки. Выходной информацией системы сейсмического мониторинга GITS являются карты сейсмической активности с указанием удароопасных зон, однако критерии ударо-опасности необходимо устанавливать конкретно для каждого шахтного поля, а обоснованные методики подобной количественной оценки отсутствуют.

В работе [35] подробно излагаются современные методы геоакустического контроля удароопасных массивов горных пород, приведено до 20 прогностических признаков критического состояния массива по параметрам и характеру проявления акустической эмиссии. Как отмечает сам автор работы, эффективность геоакустического контроля во многом зависит от объективности интерпретации результатов измерений и обоснованности используемых критериев состояния горного массива.

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

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

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

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

Эти уравнения моделируют деформационные процессы, включая разрушение. Вместе с отрицательными и положительными обратными связями и определяющими уравнениями, описывающими скорости накопления в нагружаемой среде неупругих деформаций и/или повреждений, численные решения уравнений механики деформируемого твердого тела демонстрируют все характерные черты эволюции нелинейных динамических систем, проявляющих свойства целостности и самоорганизованной критичности, как это обсуждалось выше [6, 11, 12].

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

2. Математическая эволюционная модель нагружаемой среды

Если задача механического отклика прочной среды на нагружение сформулирована как эволюционная, то решения уравнений механики деформируемого твердого тела демонстрируют все характерные этапы эволюции нелинейных динамических систем, включая состояние самоорганизованной критичности среды и развитие разрушения на разных масштабах в режиме с обострением [6, 11, 12]. Фактически мы говорим об эволюции напряженно-деформированного состояния среды, по состоянию которой на каждом этапе нагружения необходимо понять, насколько среда близка к макроскопическому катастрофическому разрушению.

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

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

Р^ = Ро^о> Р1 = ст-,- + Р^ ,

2éj = V,j + , 2rój = vu - vj,t, P = -k (éT -ép), éT = et,

Si,+Sikró, - Sk,mik = 2 Jé¿ -1 éTSj -éj. I,

Jij ' Sik kj Skj ik *-Íj 33 ^ ^ij *-ij

f (aij) = -aP + J - Y, g (a) = J2 - ЛР(2Y + aP) + const, Y = Y0 (1 - D).

Здесь и далее р0, р — начальное и текущее значение плотности материала; ¥0, V — начальное и текущее значение объема материала; го1 — компоненты вектора скорости; Р — давление; ст- — компоненты тензора напряжений; Sij — компоненты девиатора тензора напряжений; Fi — компоненты вектора массовой силы; Ю- —

• т -

ротор вектора скорости; е- — компоненты тензора скорости деформации; 8- — символ Кронекера; X — пластический множитель; 11 — первый инвариант тензора напряжений; ^ = 1/2 SijSij — второй инвариант девиатора тензора напряжений; ер — компоненты тензора скорости неупругих деформаций; 6т — скорость объемной деформации; 6р — скорость объемной неупругой деформации; К — модуль объемного сжатия; ц — модуль сдвига; а — коэффициент внутреннего трения; Л — коэффициент дилатансии; g (ст-) — пластический потенциал; D — мера поврежденности; Н(х) — функция Хевисайда; ст1п( — интенсивность тензора напряжений; ст0, ст0 — начальные значения напряжений на упругой стадии, по достижении которых в материале начинается накопление повреждений в областях сжатия и растяжения соответственно; цст — коэффициент Лоде-На-даи; ст0* — параметр модели накопления повреждений; S1, S2, Sз — главные значения девиатора тензора напряжений; стс — обобщенное кулоновское напряжение; т — интенсивность тензора напряжений; Y — сцепление.

Для учета поворота элемента среды как целого используется коротационная производная по времени Яуманна

é Р = ij

(второе уравнение). Из основного соотношения теории пластичности ép = Xdg (ст.. )/da í. следует соотношение для компонент тензора скоростей неупругих деформаций [22]:

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

Процесс накопления повреждений в среде и соответствующая деградация механических параметров описаны через меру поврежденности D = D(t, , ay) (0 < D < 1), которая зависит от инварианта напряженного состояния ay и вида напряженного состояния, определяемого через коэффициент Лоде-Надаи ца [11, 12]:

D = f №q)

D =J 2гТТ^ N ^ TT^ NN t '

a* = ao* (1 + ) ' ^a = 2

aí [Hfoa)t* + (1 - H(^a))t*] So S-

- 1' ay = aint

-aP,

^

ст0, ст0 меняются в ходе деформирования по аналогичному закону, что и в последней формуле в (8) (деградируют по мере накопления повреждений). Причем ст0 много меньше ст0, таким образом повреждения в областях растяжения-сдвига (цст < 0) начинают накапливаться при существенно меньших напряжениях, чем при цст >0 в областях сжатия-сдвига. Скорости накопления повреждений для локальных областей, где цст < 0, также существенно выше (на 1 -2 порядка), чем в областях сжатия-сдвига (цст > 0). Этот процесс дополнительно регулируется параметром ст* в (10), формируя в квазихрупкой среде существенно меньшую прочность при сдвигах-растяжениях. Функция деградации в системе обес-

I I i i 1 i'

Вышележащие породы, 270 м

Основная кровля, 3-21 м

////////////// //////////////

////////////// ///////////

3

Непосредственная кровля, 7 м : Угольный пласт, 5 м

8>

Почва, 95 м

А Д А

ОО 00

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

Основная кровля

Рис. 3. Распределение локализованных неупругих деформаций в горном массиве в случае упруго-вязкопластического отклика

Основная кровля

Почва

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

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

3. Результаты моделирования эволюции напряженно-деформированного состояния горного массива с выработками

Для демонстрации основных черт эволюции напряженно-деформированного состояния геосреды рассмотрим наиболее простую структурную организацию гора

чм1 в

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

На рис. 3 показаны распределения в горном массиве локализованных неупругих деформаций в случае упру-го-вязкопластического отклика среды и фрактального разрушения в случае упругохрупкой реакции геосреды. В первом случае в геосреде сформировались блоки разных масштабов, разделенных полосами локализованной деформации. В этих полосах сдвига среда повреждена в различной степени (0 < В < 1) в зависимости от сложившейся локальной ситуации. В данном примере наглядно демонстрируется возможность моделирования формирования фрактальных структур при численном решении уравнений механики деформируемого твердого тела.

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

Рис. 4. Распределение кулоновских напряжений в горном массиве с выработкой. Мощность основной кровли 19 м, забой отошел от монтажной камеры на 55

Рис. 5. Запись тремора напряжений в горном массиве при движении очистной выработки на нескольких датчиках

а с = т/ (7 + аР), где т — интенсивность напряжений; Р — гидростатическое давление. При приближении ку-лоновских напряжений к 1, среда приближается к критическому состоянию.

На рис. 4 приведен пример распределения в кровле и почве кулоновских напряжений для двух моментов времени, когда забой отошел приблизительно на 55 м от монтажной камеры (повреждены пока только непосредственная кровля и почва, рис. 4, а) и на заключительном этапе перед генеральным обрушением ^ ~ ^ 150 м, рис. 4, б).

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

2 ^

И

-4 -2 0 2 1ё

М/_

Ма-СТоУ/и)

-4 -2 0

2 ^

-2 -

-3 -

к. А

/ / ч. \ ч

М тг- VI 1

у

-4 -2 0 2 1ё

-2 -

-3

-4 -4

Iе %

т * \ 1 :

V V

-2

0 2 1 ё

Рис. 6. Эволюция PDF-зaвисимости тремора напряжений при движении очистной выработки

Рис. 7. Эволюция АЧХ тремора напряжений при движении очистной выработки. А — амплитуда, F — частота

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

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

Рис. 8. График повторяемости расчетных сейсмических событий. N — число событий, Е — энергия

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

На начальных этапах деформирования PDF-зави-симости имеют характерный колоколообразный вид (рис. 6, а, б). Это свидетельствует о степенном распределении расчетных сейсмических событий и слабой зависимости событий друг от друга (наблюдаются события разных масштабов). Однако по мере формирования областей повреждений больших масштабов зависимости кардинально меняют свой вид, переходя от П-образ-ного распределения (рис. 6, в), демонстрирующего наличие в треморе пакета колебаний с близкими частотами, к ярко выраженному распределению с «тяжелыми хвостами» и с явным нарушением пространственно-временной симметрии (рис. 6, г, д) при формировании в системе критического состояния перед генеральным обрушением. После обрушения PDF-зависимость приобретает вид симметричного распределения, свидетельствующего о релаксации большого количества накопленной энергии и восстановлении пространственно-временной симметрии (рис. 6, е).

Эволюция амплитудно-частотной характеристики интересна с позиции изменения ее наклона. Расчеты показывают, что при формировании критического состояния в горном массиве, предваряющего генеральное обрушение, происходит уменьшение угла наклона амплитудно-частотной характеристики. Это означает, что в достаточно большом диапазоне масштабов разрушения происходит «выравнивание» событий по выделяемой энергии. Тот факт, что процесс разрушения активизировался на всех масштабных уровнях, свидетельствует об увеличении вероятности крупномасштабных разрушений (в данном случае генерального обрушения).

График повторяемости расчетных сейсмических событий приведен на рис. 8. Сейсмические события в гор-

Рис. 9. Пространственно-временное распределение расчетных сейсмических событий в горном массиве при движении очистной выработки

ном массиве при движении очистной выработки рассчитывались следующим образом. При наращивании выработанного пространства, элементы горного массива локально теряют устойчивость и переходят в неупругое состояние, как это наглядно демонстрирует рис. 4. Происходят диссипация накопленной упругой энергии и эмиссия сейсмических волн. В расчетах фиксировался уровень выделяемой энергии при акте неупругого деформирования и тем самым оценивался класс события по известной формуле: k = lg E, где E — энергия события, которая рассчитывается из соотношения E = 3G (de pj), где depl — приращение интенсивности пластических деформаций.

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

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

Рис. 10. Пространственное распределение расчетных сейсмических событий в горном массиве при движении очистной выработки

о 24 класс О 25 класс "Йг26 класс

0.15 0.10 0.05 0.00 -0.05 -0.10 -0.15

к

#*«"И " 1

о.ю-

0.050.00 -0.05-0.10-

9 10

—I-1-1-1-1-1-г-

Рис. 11. Анализ тремора напряжений методом вейвлет-анализа. Уровни разложения d1 (а), d5 (б), d8 (в), d9 (г) соответствуют различным частотным диапазонам исходного сигнала

ных исследований. Следует заметить, что для практически всех имитационных моделей такое поведение графика повторяемости (cut-off эффекта) является типичным [2, 3, 25].

Пространственная структура расчетных сейсмических событий, обусловленных разрушениями в кровле и почве, при движении очистной выработки на заключительном этапе (L ~ 150 м) приведена на рис. 10. Главные события переместились на периферию в область сочленения кровли с невыработанным пространством (отмечено на рис. 10 звездочками). В центре кровли ярко выражена область сейсмического затишья.

При вейвлет-анализе тремора напряжений на одном из датчиков обнаружено следующее. При разбиении исходного сигнала на 9 подуровней при помощи вейвлета Добеши 10, являющегося штатным средством программы MatLab, на «датчике», находящемся в непосредственной близости к области процесса разрушения, область сейсмического затишья наблюдается на всех частотах. Появление области сейсмического затишья предваряет катастрофическое разрушение элементов горного массива, соответствующее генеральному обрушению кровли. При этом на соответствующих уровнях разложения наблюдаются пики, значительно превышающие фоновую сейсмичность по амплитуде (рис. 11). После генерального обрушения также наблюдаются пики

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

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

4. Заключительные замечания

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

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

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

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

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

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

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

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

Работа выполнена при поддержке гранта Российского научного фонда № 14-17-00198.

Литература

1. Shebalin P., Keilis-Borok V.I., Gabrielov A., Zaliapin I., Turcotte D. Short-term earthquake prediction by reverse analysis of lithosphere dynamics // Tectonophysics. - 2006. - V. 413. - Р. 63-75.

2. Шаповал А.Б., Шнирман М.Г. Универсальность алгоритмического прогноза экстремальных событий временных рядов // Информационные технологии и вычислительные системы. - 2011. - №2 4.-С. 58-65.

3. Шаповал А.Б., Шнирман М.Г. Прогноз крупнейших событий в модели образования лавин с помощью предвестников землетрясений // Физика Земли. - 2009. - № 5. - С. 39-46.

4. Козырев А.А., СеменоваИ.Э., АветисянИ.М. Создание численной

геомеханической модели месторождения «Антей» как основа прогноза напряженно-деформированного состояния массива // Горный информационно-аналитический бюллетень. - 2014. -№ 4. - С. 33-40.

5. Опарин В.Н., Сашурин А.Д., Леонтьев А.В. и др. Деструкция земной коры и процессы самоорганизации в областях сильного техногенного воздействия / Отв. ред. Н.Н. Мельников. - Новосибирск: Изд-во СО РАН, 2012. - 632 с.

6. Макаров П.В. Математическая теория эволюции нагружаемых твер-

дых тел и сред // Физ. мезомех. - 2008. - Т. 11. - № 3. - C. 19-35.

7. Ахромеева Т.С., Курдюмов С.П., Малинецкий Г.Г., Самарский А.А.

Структуры и хаос в нелинейных средах. - М.: Физматлит, 2007. -488 с.

8. Курдюмов С.П., Куркина Е.С., Потапов А.Б., Самарский А.А. Архитектура многомерных тепловых структур // ДАН СССР. -1984. - Т. 274. - № 5. - С. 1071-1075.

9. Курдюмов С.П. Режимы с обострением. Эволюция идеи / Под ред.

Г.Г. Малинецкого. - М.: ФИЗМАТЛИТ, 2000. - 312 с.

10. BakP., Tang C., WiesenfeldK. Self-organized criticality: An explanation of 1/f noise // Phys. Rev. Lett. - 1987. - V. 59. - P. 381-384.

11. МакаровП.В., Еремин М.О. Модель разрушения хрупких и квазихрупких материалов и геосред // Физ. мезомех. - 2013. - Т. 16. -№ 1. - С. 5-26.

12. МакаровП.В., ЕреминМ.О., КостандовЮ.А. Определение времени предразрушения для образцов из габбро в модели накопления повреждений // Физ. мезомех. - 2013. - Т. 16. - № 5. - С. 35-40.

13. Keilis-Borok V.I., Kossobokov V.G. Preliminary activation of seismic flow: Algorithm M8 // Phys. Earth Planetary Interiors. - 1990. - V. 61.-Р. 73-83.

14. Keilis-Borok V.I., Rotwain I.M. Diagnosis of time of increased probability of strong earthquakes in different regions of the world: Algorithm CN // Phys. Earth Planetary Interiors. - 1990. - V 61. - No. 12. - Р. 57-72.

15. КособоковВ.Г. Прогноз землетрясений: основы, реализация, перспектива // Вычислительная сейсмология. - 2005. - Вып. 36. - С. 1175.

16. МалинецкийГ.Г., ПотаповА.Б. Современные проблемы нелинейной динамики. - М.: Эдиториал УРСС, 2000. - 336 с.

17. Подлазов А.В. Теория самоорганизованной критичности — наука о сложности // Будущее прикладной математики. Лекции для молодых исследователей. - М.: Эдиториал УРСС, 2005. - С. 404-426.

18. Капица С.П., Курдюмов С.П., Малинецкий Г.Г. Синергетика и прогнозы будущего. - М.: УРСС, 2001. - 288 с.

19. Bak P., Paczuski M. Complexity, contingency and criticality // Proc. Nat. Acad. Sci. USA. - 1995. - V. 92. - Р. 6689-6696.

20. Main J. Is the reliable prediction of individual earthquakes a realistic scientific goal? // Nature. - 22 February 1999. - http://www.nature. com/nature/debates/earthquake/equake_contents.html.

21. Keilis-Borok VI. Fundamentals of Earthquake Prediction: Four Paradigms // Nonlinear Dynamics of the Lithosphere and Earthquake Prediction / Ed. by V.I. Keilis-Borok, A.A. Soloviev. - Heidelberg: Springer-Verlag, 2003. - P. 1-36.

22. Родкин М.В., Писаренко В.Ф., Нго Тхи Лы, Рукавишникова Т. А. О возможных реализациях закона распределения редких сильных землетрясений // Геодинамика и текгонофизика. - 2014. - Т. 5. -№ 4. - С. 893-904.

23. EM-DAT: The OFDA/CRED international disaster database (по состоянию на 24.10.2003). Universite Catholique de Louvain, Brussels, Belgium. - http://www.cred.be/emdat.

24. Rundle J.B., Klein W. Scaling and critical phenomena in a cellular automaton slider-block model for earthquakes // J. Stat. Phys. - 1993. -V. 72. - No. 1-2. - P. 405-412.

25. Olami Z., Feder H.J.S., Christensen K. Self-organized criticality in a continuous, nonconservative cellular automaton modelling earthquakes // Phys. Rev. Lett. - 1992. - V. 68. - No. 8. - P. 1244-1248.

26. Зуев Б.Ю., Пальцев А.И. Научно-методические основы физического моделирования нелинейных геомеханических процессов при подземной разработке полезных ископаемых // Горный информационно-аналитический бюллетень. - 2010. - № 5. - С. 18-28.

27. ЗахаровВ.Н. и др. Опыт геодинамического районирования горных районов с использованием автоматизированного линеамент-ного анализа космических снимков // Матер. XXV Межд. научн. школы им. Христиановича. - 2015. - С. 57-64.

28. Zhou J., Yang X., Yang Zh, Li H., Zhou H. Micromechanics damage modeling of brittle rock failure processes under compression // Int. J. Comput. Meth. - 2013. - V. 10. - No. 6. - P. 1350034-1-1350034-18.

29. Зуев Б.Ю., Ромашкевич А.А., Ютяев Е.П., Логинов М.А. Исследование условий работы целиков и поддержания выработок при подготовке выемочных столбов спаренными выработками // Горный информационно-аналитический бюллетень (научно-технический журнал). - 2010. - № 4. - С. 1-11.

30. Беляков В.Г., Мирошниченко Н.А., Рубцова Е.В. Математические модели динамики блочной среды в базе знаний мониторинга горных ударов // Сибирский журнал индустриальной математики. -2001. - Т. 4. - № 1. - С. 3-13.

31. Alireza M, Juhlin Ch., Wijns Ch., Urosevic M., Valasti P., KoivistoE. 3D reflection seismic imaging for open-pit mine planning and deep exploration in the Kevitsa Ni-Cu-PGE deposit, northern Finland // Geophysics. - 2012. - V. 77(5). - P. WC95-WC108.

32. Arabasz W.J., Nava S.Y., McCarter M.K., Pankow K.L., Pechmann Y.C., Ake J., McGarr A. Coal-mining seismicity and ground-shaking hazard: A case study in the trail mountain area, Emery County // Utah Bull. Seismological Soc. Am. - 2005. - V. 95. - P. 18-30.

33. Lawrence M.G., Gochioco J.R., Ruev F., Jr. Coal geophysics expands with growing global demands for mine safety and productivity // The Leading Edge. - 2012. - V. 31. - P. 308-314.

34. Engelbrecht J., Inggs M.R., Makusha G. Detection and monitoring of surface subsidence associated with mining activities in the Witbank coalfields, South Africa, using differential radar interferometry // South African J. Geology. - 2011. - V. 114. - P. 77-94.

35. Рассказов И.Ю. Геоакустические предвестники горных ударов // Вестник инженерной школы Дальневосточного федерального университета. - 2011. - № 3-4(8-9). - С. 121-143.

36. Мухамедов В.А. Фрактальные характеристики высокочастотного сейсмического шума для задач прогноза землетрясений. - Ашхабад, 2001. - 48 с. / Препринт Туркменского гос. университета.

37. Wilkins M.L. Computer Simulation of Dynamic Phenomena. - BerlinHeidelberg: Springer-Verlag, 1999. - 246 р.

Поступила в редакцию 08.06.2016 г.

Сведения об авторах

Макаров Павел Васильевич, д.ф.-м.н., проф. ТГУ, зав. лаб. ИФПМ СО РАН, pvm@ispms.tsc.ru Еремин Михаил Олегович, к.ф.-м.н., мнс ТГУ, мнс ИФПМ СО РАН, eremin@ispms.tsc.ru