CC BY
52
удк 550 34 а.в. Азаров, А.С. Сердюков, А.В. Яблоков
МЕТОДИКА ОПРЕДЕЛЕНИЯ МЕХАНИЗМОВ ОЧАГОВ МИКРОСЕЙСМИЧЕСКИХ СОБЫТИЙ НА ОСНОВЕ МОДЕЛИРОВАНИЯ ПОЛНЫХ ВОЛНОВЫХ ПОЛЕЙ В ГОРИЗОНТАЛЬНО-СЛОИСТЫХ СРЕДАХ*
Рассмотрена проблема сейсмического мониторинга геодинамических процессов, протекающих в массиве горных пород входе разработки месторождений полезных ископаемых. Решается задача восстановления механизма очага точечного микросейсмического источника по данным площадных поверхностных систем наблюдения. Предполагается, что координаты гипоцентра уже найдены с использованием известных подходов локации микросейсмических событий. Предложенная методика основана на сопоставлении зарегистрированного поля с раннее насчитанными полями от источников с различными возможными механизмами. Для моделирования волновых полей использован метод «wavenumber integration» позволяющий эффективно рассчитывать синтетические сейсмограммы в случае горизонтально-слоистых сред. Приведены результаты восстановления механизма точечного источника в синтетической скоростной модели углепородного массива. Исследована устойчивость восстановления механизма к ошибкам в задании скоростной модели среды. Показано, что в случае неточного задания скоростной модели вмещающей среды необходимо учитывать сдвиги во временах пробега сейсмических волн. Ключевые слова: микросейсмический источник, механизм источника, тензор сейсмического момента, уравнения упругости, численное моделирование, обратная задача.
Введение
Микросейсмический мониторинг (МСМ) является методом дистанционного контроля состояния массива горных пород. Технологии МСМ основаны на регистрации и обработке сейсмических полей, генерируемых источниками, кото-
* Работа выполнена при поддержке РФФИ, проект № 16-35-00513 мол_а.
ISSN 0236-1493. Горный информационно-аналитический бюллетень. 2016. № 10. С. 131-143. © 2016. А.В. Азаров, А.С. Сердюков, А.В. Яблоков.
рые могут иметь различные механизмы. Вопросы определения и исследования механизмов очагов микросейсмических событий занимают в обработке данных микросейсмического мониторинга особое положение, так как они несут в себе важную информацию о геодинамических процессах, которую возможно использовать при контроле разработки твердых полезных ископаемых.
В большинстве существующих работ рассматривается задача локации гипоцентров событий. Предложен целый ряд различных алгоритмов. В тоже время, восстановление механизмов очагов событий так же важно, так как позволяет судить о природе геодинамических процессов следствием которых является наблюдаемая сейсмичность. В ряде ситуаций для уточнения результатов локации гипоцентров событий так же важно учитывать механизмы источников. Например, в случае мониторинга большой области массива горных пород наземными системами наблюдения обычно используют алгоритмы получившие название сейсмо-эмиссионной томографии (СЭТ) [1, 2]. Они основаны на суммировании сейсмограмм, и дают объемную плотность интенсивности излучения во всей контролируемой области породного массива и предназначены для обработки данных от множества микросейсмических событий различной интенсивности. В большинстве существующих алгоритмов СЭТ не учитываются различные механизмы микросейсмических источников. Обычно предполагается, что они представляют собой центры расширения. Это может приводить к неверным результатам из-за не учета фазовых характеристик диаграммы направленности источника. Так, при деформации горных пород могут возникать источники типа сдвиговой дислокации, волновое поле от которых значительно отличается от поля источников типа центр расширения [3].
В данной работе рассматривается задача восстановления механизма источника в горизонтально слоистой среде. Предлагается методика определения механизма источника с известной локацией и исследуется устойчивость алгоритма к ошибкам в скоростной модели среды. Описан подход к моделированию упругих волновых полей, позволяющий эффективно рассчитывать синтетические сейсмограммы в горизонтально слоистых средах.
Восстановление механизма источника
Запишем систему уравнений упругости:
д
р — и = V • с,
д
-a = XV- uI + |~(Vu) + (Vu)T], (1)
dt L J
где u — скорости смещения среды, a — напряжения.
В работе рассматриваются точечные модели источников, для описания которых можно использовать понятие тензора сейсмического момента M [4]. В трехмерном случае M имеет шесть компонент, однозначно определяющие механизм источника. Задание компоненты M.. данного тензора означает, что к точке прикладывается пара сил в направлении i с моментом M.. и с плечом в направлении j. Введение источника заданного тензором M в систему (1) осуществляется добавлением во второе уравнение члена Mb(t - t0)d(x - x0) , где х0 определяет координаты источника, а t0 время срабатывания. В работе решается обратная задача по восстановлению механизма очага, т.е. определении тензора сейсмического момента M по известному полю u. В рамках работы полагается, что локация источника х0 известна.
Методика восстановления механизма состоит из двух основных этапов. Первый этап заключается в вычисление шести полей от источников с тензором сейсмического момента имеющего одну ненулевую компоненту M... Второй этап заключается в поиске комбинации шести насчитанных полей дающих в результате зарегистрированное поле. Ниже приведен алгоритм реализации предлагаемой схемы:
1. Расчет набора волновых полей uM . uM — поле от источника, имеющего только одну ненулевую компоненту M.. в тензоре сейсмического момента M. В трехмерном случае необходимо рассчитать 6 полей: uMu, uM^, um33 , uMa, uM^, uM3 .
2. Выполнение сдвигов по времени сеймотрасс от uM^ , чтобы учесть изменения времен прихода волн при наличии ошибки в модели среды или в локации источника.
3. Построение usum = ^ wM,, uMi, , где wM,, — вес с которым
м S
ij
uM входит в u , u — комбинация насчитанных полей. п
Mij sum sum ц ц 2
4. Построение и минимизация функционала \\uregSsstr - u I , где uregistr — зарегистрированное поле от источника, механизм которого необходимо определить.
II II2
Функционал \\uregSstr - u I представляет собой многочлен второй степени от переменных wM,, , число которых в данном случае равно шести. Для нахождения точки минимума полученного многочлена продифференцируем его по каждой пере-
менной-\игеаЫг - и\|2 = 0 . Отсюда получаем шесть линей-
дтм " у 112
ных уравнений на нахождение шести неизвестных . После решения системы линейных уравнений получаем веса ^ , с которыми поля от источников с ненулевой компонентой М в тензоре сейсмического момента входят в комбинацию и8ит. Таким образом, коэффициенты шм,, определяют тензор сейсмического момента, то есть механизм источника.
Также в работе ставится задача исследования устойчивости алгоритма к ошибкам в модели среды. С этой целью проводится восстановление источника в моделях с измененными упругими параметрами. Это может привести к изменению времен прихода волн на приемники, из-за чего минимизация функцио-
II 1|2
нала ||игедЫг - иит|| не даст корректного результата. Поэтому
в предлагаемом алгоритме восстановления механизма реализован этап, на котором сейсмотрассы, насчитанные в разных моделях, сдвигались относительно друг друга для нахождения максимум корреляции сигналов. Коэффициент корреляции k двух трасс а и Ь находится как я
ЕаЬ
к = , ы (2)
Е ь2
i=i V i=i
где N — число отчетов по времени. После чего применялся соответствующий сдвиг по времени и выполнялись шаги 1—3. Аналогично можно восстанавливать механизм, когда координаты источника содержат ошибки. В этом случае также возникают изменения времен прихода волн, что может сделать невозможным определения механизма.
Прямое моделирование волновых полей
В работах [5, 6] моделирование сейсмических полей от источников с различными механизмами в сложно-построенных средах проводилось конечно-разностными методами. В представленной работе рассматривается горизонтально слоистая среда, в которой моделирование волновых полей целесообразно проводить с использованием методики «wavenumber integration» [7, 8]. При использовании данной методики считается, что каждый слой является однородным и изотропным. Введение цилиндри-
ческой системы г, 0} координат и рассмотрение сейсмического источников распределенных на оси цилиндра, позволяет понизить размерность задачи, так как решение в таком случае становится независимо от координаты 0, что приводит к снижению вычислительных затрат и повышению скорости расчета синтетических сейсмограмм. Решение в каждом слое представляется в виде суперпозиции поля генерируемого источником внутри слоя при отсутствии границ и неизвестных полей, которые необходимы, чтобы удовлетворить граничным условиям на поверхностях раздела сред. В каждом слое решается задача упругости, решение которой можно выразить в терминах скалярных потенциалов ф, Л, у удовлетворяющие волновым уравнениям:
(V2 + к2)ф = 0,
(V2 + ¿2)(Л, у) = 0, где h и k волновые числа для продольных и поперечных волн соответственно. Скалярный потенциал ф характеризует поле продольных волн, скалярные потенциалы Л и у характеризуют поле поперечных волн. В цилиндрической системе координат решения уравнений (2) можно записать в следующем виде:
Ф(г, в, г) = £У (г,г) {""<» йеП,
т=0 [ $,т(т в)
т, ч Гсоз(т в)
Л(г, в, г) = £Лт (г, г) \ к I, (4)
т=0 [ 81п(т в)
у(г,в,г) = ±Г(г,г)| -п<твП,
т=0 [- со8(т в)
Подставив представление (4) в (3) и применив преобразование Ханкеля по переменной г получим выражения для фт, Лт, ут:
да
Фт (г, г) = |[< (з)в -га(з) + < (з)вга(з) (гз)йз,
0
да
Лт (г, г) = |[Ь1т (з)в~гр(з) + Ь2т (5)вгр(з) ]з/т (гз)йз, (5)
0
да
ут (г, г) = {[< (з)в-гр(з) + С2т (з)вгв(з) ]з^ (гз)йз,
0
где Jm — Бессель функции порядка т, 5 — горизонтальное волновое число,
a(s) = PCs) =
[ (s2 - й2)1/2,при s2 > Re(h2), [i(h2 - s2)1/2,npè s2 < Re(h2),
f Cs2 - ¿2)1/2,при s2 > ReCk2),
[¿Ck2 - 52)1/2,при s2 < ReCk2),
Чтобы решить рассматриваемую задачу упругости в горизонтально-слоистой среде необходимо определить коэффициенты
m m 7 m т m m m
aj , a2 , bj , b2 , cj , c2 , удовлетворяющие граничным условиям, на отыскание которых составляется система линейных уравнений.
Несмотря на то, что в рамках методики «wavenumber integration» вычисляются решения, не зависящие от угла ф, с помощью нее можно рассчитывать поля от источников с любыми механизмами. Для этого, как показано в работе [9] необходимо вычислить 10 функций Грина (ZSS, ZDS, ZDD, RSS, RDS, RDD, TSS, TDS, ZEP, REP) не зависящие от ф, из которых можно собрать любой механизм. Каждая функция Грина соответствует одной из компонент смещения от определенного механизма источника. В табл. 1 приведено описание всех 10 функций.
Используя приведенные функции Грина, поле от источника с произвольным тензором сейсмического момента записывается следующим образом:
uz (r, z = 0, ю, ф) = My
ZSS , ZDD ZEX -cos(2 ф)---1--
Таблица 1 Функции Грина
Название функции Грина Компонента смещения Механизм источника
ZSS вертикальная сдвиг (скольжение по простиранию)
ZDS вертикальная сброс под углом 90°
ZDD вертикальная сброс под углом 45°
RSS радиальная сдвиг (скольжение по простиранию)
RDS радиальная сброс под углом 90°
RDD радиальная сброс под углом 45°
TSS тангенсальная сдвиг (скольжение по простиранию)
TDS тангенсальная сброс под углом 90°
ZEP вертикальная центр расширения
REP радиальная центр расширения
+м.
22
-ZSS ч ZDD ZEX соз(2ф) - -
М
33
ZDD ZEX
2 '63
+М12 [ZSS 81п(2ф)] + М13 [ZDScos(ф)] + М23 [ZSSs1п(ф)]
ЯЕХ
щ (г, г = 0, ю, ф) =
RSS ч RDD
-С08(2ф) -
+м.
22
-RSS ч RDD REX
-С08(2ф)---1--
м
6 3 RDD REX
33
+М12 [RSSs1п(2ф)] + М13 [RDScos(ф)] + М23 [RSSs1п(ф)]
(6)
и (г, г = 0, ю, ф) = м1
81п(2ф)
■м.
22
81п(2ф)
+М12 [-ГSScos(2ф)] + М13 [ГDSs1п(ф)] + М23 [-ГDScos(ф)]
где М{. компоненты тензора сейсмического момента, угол ф определяет азимут на приемник.
Чтобы получить синтетические сейсмограммы необходимо на первом шаге вычислить 10 функций Грина, на втором шаге построить поля в каждом приемнике для конкретного механизма по формулам (5), и на третьем шаге свернуть поле с заданной формулой импульса.
Исследование алгоритма восстановления механизма
источника
В данном разделе представлены результаты реализации и исследования предлагаемого алгоритма восстановления меха-
источник А - приемник
Рис. 1. Горизонтально-слоистая модель среды (а), график скоростей продольных, поперечных волн и плотности (б)
Таблица 2
Точная модель среды
Номер слоя Мощность, м Ур, м/с Уs, м/с р, кг/м3
1 5 600 300 1600
2 5 1000 500 1750
3 10 1600 900 1800
4 140 2300 1350 2000
5 30 2700 1600 2200
6 10 1700 1000 1360
7 0 2300 1350 2000
низма источника на синтетических данных. Для расчетов синтетических сейсмограмм использовалась трехмерная область, моделирующая углепородный массив, которая представляла собой куб со стороной 300 м (рис. 1).
На рис. 1, б графики упругих волн и плотности от глубины. Верхняя горизонтальная прямая — шкала скорости, нижняя — шкала плотности. В табл. 2 приведены упругие параметры горизонтально-слоистой среды. Два верхних низкоскоростных слоя моделируют земную поверхность, слой под номером 6, условно считался угольным.
В качестве источника был выбран чистый сдвиг с тензором сейсмического момента
Г0 1
м = м0
1 0 0 0 0 0
где М0 — скалярный сейсмический момент, определяемый для сдвига как М0 = ц и Б , ц — модуль сдвига, и — смещение на разрыве, S — площадь поверхности разрыва. Для расчетов М0 принимался равным 1. Данный источник устанавливался на глубину 195 м, то есть непосредственно в угольный пласт. Форма сигнала бралась в виде импульса Риккера с центральной частотой 100 герц. Регистрация сейсмического поля производилась наземной площадной системой наблюдения. Система наблюдения состояла из 100 приемников равномерно разбросанных по площади поверхности. Далее производилось моделирование сейсмического поля для рассматриваемого источника и модели среды с использованием методики «wavenumber
Рис. 2. Смещение компоненты z на одном из приемников
integration». На рис. 2 представлено смещение компоненты z зарегистрированное одним из приемников расположенным на поверхности:
Также на рис. 3 представлены сейсмограмма, полученная с поверхностной системы наблюдения для рассматриваемого случая.
После проведения прямого моделирования производилось восстановление данного источника по предложенному алгоритму. В табл. 3 показаны результаты восстановления. В данном случае компоненты тензора сейсмического момента восстановились точно.
025 03
время, с
Рис. 3. Сейсмограмма, полученная с поверхностной системы наблюдения
Таблица 3
Восстановление механизма источника
Компоненты тензора сейсмического момента Значения компонент для сдвига Восстановленные компоненты в «реальной» модели Восстановленные компоненты в «испорченной» модель (без учета изменения времени прихода волны) Восстановленные компоненты в «испорченной» модель (с учетом изменения времени прихода волны)
Мп 0,00 0,00 -0,1049 0,0081
м12 1,00 1,00 -0,4890 0,6954
Мв 0,00 0,00 0,0000 0,0000
М22 0,00 0,00 0,1049 -0,0081
М23 0,00 0,00 0,0000 0,0000
М33 0,00 0,00 0,0000 0,0000
Далее в работе производилась проверка устойчивости алгоритма восстановления механизма к ошибкам модели среды. С этой целью была рассмотрена модель среда с ошибками, где изменены упругие параметры двух верхних низкоскоростных слоев, в которых были понижены скорости продольных и поперечных волн (табл. 4).
Используя модель с ошибками, насчитывались поля от источников имеющих одну ненулевую компоненту в тензоре сейсмического момента. В точной модели вычислялось поле от рассматриваемого источника типа сдвиговой дислокации. Таким образом, моделировалась реальная практическая ситуация, ког-
Таблица 4
Модель среды с ошибками
Номер слоя Мощность, м Vp, м/с \&, м/с , кг/м3
1 5 500 260 1500
2 5 900 450 1750
3 10 1600 900 1800
4 140 2300 1350 2000
5 30 2700 1600 2200
6 10 1700 1000 1360
7 0 2300 1350 2000
да точная модель среды не известна. В табл. 3 в четвертой колонке показан результат восстановления механизма источника в данном случае. Как видно, в этом случае механизм источника не смог восстановится. Главной причиной такого результата является изменение времен прихода волн в приемники. В этом случае на первом этапе проводился поиск необходимого сдвига сейсмотрасс по времени путем поиска максимума корреляции между наблюдаемыми и рассчитанными сейсмограммами, после чего выполнялись шаги 1—3 предлагаемого алгоритма. В табл. 3 приведен результат восстановления рассматриваемого источника с учетом изменения времени прихода волн на приемники. В этом случае механизм восстановился достаточно точно, но значение скалярного момента стало более низким, по сравнению с исходным источником.
Выводы
В работе был предложен и протестирован на синтетических данных алгоритм восстановления механизма очага сейсмического источника в горизонтально слоистой среде с использованием наземной системы наблюдения. Описан способ эффективного расчета синтетических сейсмограмм для слоистых сред и продемонстрированы результаты прямого моделирования сейсмического источника с механизмом чистый сдвиг. Показана важность учета изменения времен приходы волн, при не точном задании скоростной модели вмещающей среды. Полученные результаты могут быть использованы при решении задач сейсмического мониторинга геодинамических процессов, протекающих в массиве горных пород.
СПИСОК ЛИТЕРАТУРЫ
1. Хогоев Е. А., Колесников Ю.И. Применение сейсмоэмиссионной томографии для изучения геодинамически активных зон // Seismic Technology. - 2011. - Т. 8. - № 1. - С. 59-65.
2. Chebotareva I. Y. New algorithms of emission tomography for passive seismic monitoring of a producing hydrocarbon deposit: Part I. Algorithms of processing and numerical simulation // Izvestiya, Physics of the Solid Earth. - 2010. - Т. 46. - № 3. - Pp. 187-198.
3. Маловичко Д. А. Изучение механизмов событий в рудниках Верхнекамского месторождений калийных солей: дис. ... на соискание степени к.ф.-м.н. - M.: ОИФЗ, 2004. - 178 с.
4. Aki K, Richards P. Quantitative Seismology: Theory and Methods. San Francisco, 1980. - 520 p.
5. Азаров А. В., Сердюков А. С., Никитин А. А. Численное моделирование микросейсмического шума в блочно-иерархических неоднород-
ных средах // Горный информационно-аналитический бюллетень. — 2015. - № 9. - С. 328-335.
6. Курленя М. В., Сердюков А. С., Азаров А. В., Никитин А. А. Численное моделирование волновых полей от микросейсмических событий при подземной добыче полезных ископаемых // Физико-технические проблемы разработки полезных ископаемых. — 2015. — № 4. — С. 61-69.
7. Lucifredi I., Ishii M. A coupled wavenumber integration approach for calculating the wavefield in large-scale laterally varying structures // Geophysical Journal International. - 2012. - Т. 188. - № 1. - Pp. 195-210.
8. Schmidt H., Glattetre J. A fast field model for three dimensional wave propagation in stratified environments based on the global matrix method // The Journal of the Acoustical Society of America. - 1985. - Т. 78. - № 6. -Pp. 2105-2114.
9. Jost M. L., Herrmann R. B. A student's guide to and review of moment tensors // Seismological Research Letters. - 1989. - Т. 60. - № 2. -Pp. 37-57. ЕШ
КОРОТКО ОБ АВТОРАХ
Азаров Антон Витальевич - младший научный сотрудник, аспирант, Институт горного дела им. Чинакала СО РАН, e-mail: antonazv@mail.ru,
Сердюков Александр Сергеевич - кандидат физико-математических наук, младший научный сотрудник, Институт нефтегазовой геологии и геофизики им. А.А.Трофимука СО РАН, e-mail: aleksanderserdyukov@ya.ru, Яблоков Александр Викторович - магистрант, Новосибирский государственный университет, e-mail: sasha.yablokov@mail.ru.
UDC 550.34
Gornyy informatsionno-analiticheskiy byulleten'. 2016. No. 10, pp. 131-143. A.V. Azarov, A.S. Serdyukov, A.V. Yablokov METHODS OF THE FOCAL MECHANISM DETERMINATION OF MICROSEISMIC EVENTS BASED ON MODELING FULL WAVE FIELDS IN A HORIZONTALLY STRATIFIED MEDIA
This paper is dedicated to the seismic monitoring of geodynamic processes occurring in the rock mass in mining. One of the work's aims is to recover mechanism of a point micro-seismic source using the data from areal surface observation systems. We suppose that the hypocenter's coordinates have been found with the help of the well-known methods determining the locations of microseismic events. The method suggested in the paper is based on the comparison of the registered field and the fields calculated by the sources with various mechanisms. In order to design wave fields we have used «wavenumber integration» method allowing us to calculate synthetic seismograms for a horizontally stratified media. The paper presents the results of the point source mechanism recovery in the coal bed's synthetic ve-
locity model. Moreover, we have surveyed the mechanism recovery stability to faults in determination of the environment velocity model. Thus, in case when the velocity model is not determined accurately, it is necessary to take into account seismic waves' travel time shifts. To calculate the shifts one can find the maximum correlation between the observed and calculated seismograms.
Key words: microseismic source, source mechanism, seismic moment tensor, elastic equations, numerical modeling, inverse problem.
AUTHORS
Azarov A.V., Junior Researcher, Graduate Student, Chinakal Institute of Mining of Siberian Branch of Russian Academy of Sciences, 630091, Novosibirsk, Russia, e-mail: antonazv@mail.ru,
Serdyukov A.S., Candidate of Physical and Mathematical Sciences,
Junior Researcher, e-mail: aleksanderserdyukov@ya.ru,
Trofimuk Institute of Petroleum Geology and Geophysics of Siberian,
Siberian Branch of Russian Academy of Sciences,
630090, Novosibirsk, Russia,
YablokovA.V., Master's Degree Student,
Novosibirsk State University,
630090, Novosibirsk, Russia,
e-mail: sasha.yablokov@mail.ru.
ACKNOWLEDGEMENTS
The study has been supported by the Russian Foundation for Basic Research, Project No. 163500513mol_a.
REFERENCES
1. Khogoev E. A., Kolesnikov Yu. I. Seismic Technology. 2011, vol. 8, no 1, pp. 59—65.
2. Chebotareva I. Y. Izvestiya, Physics of the Solid Earth. 2010, vol. 46, no 3, pp. 187— 198.
3. Malovichko D. A. Izuchenie mekhanizmov sobytiy v rudnikakh Verkhnekamskogo mestorozhdeniy kaliynykh soley (Investigation of mechanisms of events in mines of the Upper Kama Potash Salt Deposit), Candidate's thesis, Moscow, OIFZ, 2004, 178 p.
4. Aki K., Richards P. Quantitative Seismology: Theory and Methods. San Francisco, 1980.520 p.
5. Azarov A. V, Serdyukov A. S., Nikitin A. A. Gornyy informatsionno-analiticheskiy byul-leten'. 2015, no 9, pp. 328-335.
6. Kurlenya M. V., Serdyukov A. S., Azarov A. V, Nikitin A. A. Fiziko-tekhnicheskie prob-lemy razrabotki poleznykh iskopaemykh. 2015, no 4, pp. 61-69.
7. Lucifredi I., Ishii M. Geophysical Journal International. 2012, vol. 188, no 1, pp. 195-210.
8. Schmidt H., Glattetre J. The Journal of the Acoustical Society of America. 1985, vol. 78, no 6, pp. 2105-2114.
9. Jost M. L., Herrmann R. B. Seismological Research Letters. 1989, vol. 60, no 2, pp. 37-57.
A_
