2.5D МОДЕЛИРОВАНИЕ СЕЙСМИЧЕСКИХ КОЛЕБАНИЙ В УГЛЕПОРОДНОМ МАССИВЕ Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

ISSN 0136-4545

^Курнал теоретической и прикладной механики.

№1(58) / 2017.

ГЕОМЕХАНИКА, РАЗРУШЕНИЕ ГОРНЫХ ПОРОД, РУДНИЧНАЯ АЭРОГАЗОДИНАМИКА И ГОРНАЯ ТЕПЛОФИЗИКА

УДК 550.834:622.12 ©2017. А.А. Глухов

2.5Б МОДЕЛИРОВАНИЕ СЕЙСМИЧЕСКИХ КОЛЕБАНИЙ В УГЛЕПОРОДНОМ МАССИВЕ

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

Ключевые слова: метод пластовой сейсморазведки, каналовые волны, 2,5Б модели.

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

Полное трехмерное моделирование не представляет математической сложности, однако требует огромных ресурсов машинного времени и объема оперативной памяти. Поэтому его можно практически реализовать только на суперкомпьютерах. Компромиссным решением является использование модели среды, параметры которой не изменяются вдоль некоторого направления (2.5D - моделирование).

В настоящей статье кратко изложен алгоритм 2,5D моделирования, разработанный в РАНИМИ специально для задач шахтной пластовой сейсморазведки.

2. Математическая модель. Алгоритм 2.5D моделирования процесса распространения сейсмических колебаний в среде, параметры которой не изменяются

Автор благодарит Кутепова Д. В. и Ребенка Е. В. за оформлении статьи к публикации

вдоль заданного направления базируются на использовании системы уравнений в виде

д2и д

рй = (\ + 2у)—2 + (\ + у) —

ду д—

ду дх

д2у д I ди д—

ру = (Л + 2у) + (Л + у) — ( — +

д2— д рт = (А + 2у) + (Л + у) —

дг

ди ду

дх ду

+ у

+ у

+ у

'д2 и д2и^

<9у2 у

д2у д2уА

д2— д2—

дх2 ду2

(1)

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

=

Л + 2у

р

^ =

(2)

Источник колебаний моделируется последовательным набором значений п(Ь), ю(Ь) и «(^.описывающие х, у и г компоненты смещений среды в заданной точке.

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

Рис. 1. Модель среды, структура и характеристики пород которой не изменяются вдоль заданного

направления.

К уравнениям системы применяется преобразование Фурье по схеме

1 оо

й = —= / ие~гку(1у = и, л/2п-о 1 оо

ди

— = гк — ду \/2тт -

/ ие гкуdy = гки,

(3)

д2 и

1 оо

_= -к2

ду'2 У27г-

/ ие-гкуdy = -к2и,

где к = ку - волновое число вдоль направления у. Вещественную часть получаемой в результате преобразований системы уравнений можно записать в виде

д2и д дШ ( 0 д2и\

ри =(Л + ^ +(л + ^ аГаГ + " ~к и +

рУ = -к2 (Л + 2р) V + р

д2У д2УЛ ■ +

дх2 дх2

(4)

д2Ш д ди (д2Ш 0 ^

= <Л + 2"> + (Л + ''> & Л + '' Ь? " к №'

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

А

Ат п —

(Лтп + 2р

рт п Ь>2

2

рт п Т

13т п —

рт пМ

Ст п — 4

1 (Лт п + рт п) Т

От

Е =

Ет п —

рт пЬ2 рт пк Т

рт п

(Лтп + 2р

рт п

Рт п — 2(1 2Вт п) Ет п) &т п — 2(1 Ат п Вт п От п) )

уравнения системы (4) можно переписать в конечно-разностном виде

(5)

ттР+1 = Г1 ТТР - ТТР-1 + А (ттр + ТТР \ +

иmn — ^тпитп итп + Ат п \и т+1 п + ит-1 п) +

+Втп (иРп+1 + и^пп-!) + Стп № +1 п+1 ^т — 1 п+1 ^т +1 п—1 + ^т-1 п—1) ,

= От п - + Ат п № п+1 + ^ п—1) + (6)

+Вт п (№Р+1 п + —1 п) + Ст п (иР+1 п+1 — Рт —1 п+1 — Рт+1 п—1 + Рт —1 п—1) ,

Т/Р+1 — И (УР + У/Р + УР + УР ) + Р УР — Л/Р—1

Ут п — Вт п т+1 п + Ут—1 п + Ут п+1 + Ут п—1) + Рт п Ут п Ут п ■

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

и — ие—ах у — уе—ау — — -в—аг (7)

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

иР+п — От пит п - Рит п + Ат п (иР+1 п + Рт —1 п) + Вт п {РР п+1 + Рт п— 1) + +Ст п (№т+1 п+1 — —1 п+1 — №т+1 п—1 + —1 п—1) ,

— О'т п№Рп - тт—п + Ат п № п+1 + №Рп—1) + (8)

+Вт п (^^+1 п + 1 п) + Ст п (иР+1 п+1 — Рт—1 п+1 — Рт+1 п—1 + Рт—1 п—1) ,

т/Р+1 _ д (лгР I ЛТР I ЛТР I т/Р \ I р/ т/Р _ Я\гР— 1

ут п — Вт п т+1 п + У т—1 п + Ут п+1 + Ут п—1) + Р тпут п Рут п

где в — 1 + ат,

п — 2 (1 - 2Втп) - Етп + ат - (ат)2,

От п — 2 (1 - Ат п - Вт п - Вт п) + ат - (ат)2 ■ (9)

Нетрудно заметить, что система соотношений (8) для к=0 представляет собой конечно разностную схему решения 2D задачи в плоскости у=0. Далее для ограниченной заданным периодом времени последовательности моментов времени р вычисляем соотношения (8) для дискретного набора равномерно распределённых ^ =]А.к и получаем решение при у=0 по формуле

1 гж Ак

и(х, О, г, = — / (1ки{х, к, г, ¿)егА''° рз — ^ и(х, 2п ¿-со п Г^

_ - (Ю)

п 1=0

Поскольку и (х,к, г, Ь) это чётная функция к, поле и(х,0,г,Ь) можно получить, суммируя вклады для всех положительных к1. При практической практических целей мы должны использовать конечное число волновых чисел, ограниченное предельным значением ктах, которое определяется по условию Найквиста (ктах < п/у). Величина дискретизации для волнового числа определяется из свойства дискретного преобразования Фурье: Ак = п/ (утах — утш). Программные модули реализующие 2.5D моделирование процесса распространения сейсмических колебаний в среде, параметры которой не изменяются вдоль заданного направления базируются на разработанном в РАНИМИ программном комплексе моделирования сейсмических колебаний в углепородном массиве [3].

На первом этапе для каждой ячейки счетной решетки согласно (5) и (9) вычисляются коэффициенты соотношений (8).

На втором этапе итерационно производится расчет последовательных стадий процесса распространения волнового поля:

- для каждого к1 по соотношениям (8) вычисляются значения ,УР+1 ,ЖР+1 для каждого узла расчетной решетки;

- по формуле (10) путем суммирования решений для каждого к1 вычисляются поля смещений волнового поля вдоль направлений х, у и г в плоскости у=0 (См. рис. 2).

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

Максимальное время расчета 1000 итераций для модели 500 х 500 х 100 - 4 часа 15 минут.

3. Результаты расчетов. На рис.3 для сравнения представлены результаты 2D и 2.5D расчета процесса распространения сейсмоакустических колебаний для модели, представляющей собой упругий слой, расположенный параллельно плоскости (ху) в более плотных упругих полупространствах. Ось г выбрана перпендикулярно слою, инвариантное направление выбрано вдоль оси у. Толщина слоя равна 1.5м. Параметры слоя: У^=1900 м/с, У51=900 м/с, р1=1800 кг/м3. Параметры для полупространств: Ур2=3500 м/с, У^2=1900 м/с, р1=2700 кг/м3. Параметры слоев соответствуют условиям залегания угольного пласта в алевролите, что типично для углей марки Ж в условиях Донецкого угольного бассейна [1].

На рисунках четко различимы несколько типов волн. Две волны имеют фронты в виде концентрических окружностей. Это классические волны сжатия и сдвига, распространяющиеся со скоростями Ур и У3 волн во вмещающей породе, соответ-

Рис. 2. Поля смещений волнового поля вдоль направлений х, у и г в плоскости у=0 вычисляются по формуле (10) путем суммирования решений для каждого к^.

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

Сравнение результатов расчета для 2D и 2.5D модели (см. рис. 3) показывает, что перечисленные типы волн на качественном уровне практически не различаются.

Оба результата содержат и цуги нормальных волн. Они образуются в слое-волноводе в результате интерференции колебаний, многократно отраженных от границ раздела сред. В шахтной сейсморазведке такие волны называют каналовыми и широко используются. В случае 2D модели цуг имеет скорость распространения около 1700 м/с. Классических колебаний, которые должны в пласте такой мощности распространяться со скоростью близкой к значению У81 не наблюдается. Строго говоря, каналовая волна присутствует, но амплитуда составляет около 10% от максимальной по области расчета, что мало для проявления в этом режиме показа. В случае 2.5D модели начальная часть цуга располагается там же, однако четко про-

2D 2.5D

Рис. 3. Результаты 2В и 2.5В расчета процесса распространения сейсмоакустических колебаний для упругого слоя в более плотных упругих полупространствах.

является и теоретически ожидаемая низкоскоростная составляющая. Для данных параметров угольного пласта преимущество использования 2.5 D модели очевидно.

Следует отметить, что в практике моделирования было отмечено, что 2D решение достаточно точно передает структуру поля в угольных пластах в плане числа волновых пакетов, скоростей их распространения. В то же время соотношения амплитуд различных волновых пакетов всегда вызывали вопросы [1,2]. Возможно, использование 2.5D моделирования позволит решить эту проблему.

На данном этапе работы расчеты были проведены на 12 различных моделях. Варьировались параметры угля и вмещающих пород в диапазоне значений, характерных для условий Донбасса. Во всех случаях, когда мощность пласта составляла 1.5м и более, при 2.5D моделировании проявлялись высокоамплитудные каналовые волны . Этот факт соответствует реалиям шахтной сейсморазведки [1]. В то же время, амплитула каналовой волны при 2D моделирование не превышала 50% амплитуды сигнала в целом.

4. Заключение. Таким образом, применение 2.5D подхода позволяет в сравнении с 2D решением более точно моделировать характеристики волновых полей, в частности - амплитуды каналовых волн, используемых в шахтной пластовой сейсморазведке для прогноза геологических нарушений. Это важный в методическом плане результат, поскольку на анализе аплитудных характеристик базируется ряд прогнозных методик, в частности, сейсмическая томография [1].

1. Анциферов А.В. Теория и практика шахтной сейсморазведки / А.В Анциферов. - Донецк: ООО "Алан", 2002. - 312 с.

2. Анциферов А.В. Математическое моделирование в шахтной сейсморазведке / А.В Анциферов., А.А. Глухов - Киев: Наук. думка, 2012. - 255 с.

3. Глухов А.А. О программном комплексе моделирования распространения сейсмических колебаний в угленосной толще / А.А. Глухов // Науков1 пращ Донецького нащонального техшчно-

го ушверситету. Сер. Обчислювальна технжа та автоматизащя. - Донецьк: ДонНТУ, 2005. -Вип. 88. - С. 106-113.

Glukhov A.A.

2.5-D Seismic Velocity Modeling in Coal-Rock Mass.

The article deals with 2.5-D modeling scheme for seismic-wave propagation in coal-rock mass adapted for the problems of in-seam seismic survey in mines. It is shown that in comparison with 2-D solution the employment of 2.5-D approach enables to simulate wave field attributes more precisely, including channel wave amplitudes used to predict geologic faults in coal seams. Keywords: in-seam seismic method, channel waves, 2,5D models..

Республиканский академический научно-исследовательский и Получено 26.12.16

проектно-конструкторский институт горной геологии, геомеханики, геофизики и маркшейдерского дела (РАНИМИ), Донецк

Glukhov1964@yandex.ru