CC BY
10УДК550.34:531.2:519.6
ЧИСЛЕННОЕ РЕШЕНИЕ ДВУХ СЕЙСМИЧЕСКИХ ЗАДАЧ СЕТОЧНО-ХАРАКТЕРИСТИЧЕСКИМ МЕТОДОМ В ПОЛНОЙ ТРЕХМЕРНОЙ ПОСТАНОВКЕ
В. И. Голубев 1 2, И. Б. Петров 1 2
1 Московский физико-технический институт (государственный университет), 2 Федеральный научный центр Научно-исследовательский институт системных исследований РАН, w.golubev@mail. ru, petrov@mipt. ru
Сейсмическая разведка является наиболее распространенным методом поиска месторождений нефти и газа. Ввиду открытия нетрадиционных коллекторов важной задачей является разработка новых процедур обработки данных полевых работ. Развитие высокопроизводительных вычислительных систем позволяет с высокой точностью проводить моделирование процессов распространения сейсмических волн в гетерогенных геологических средах для формирования синтетических сейсмограмм. В настоящей работе использован сеточно-характеристический метод на гексаэдральных расчетных сетках, позволяющий корректно описать криволинейные границы между геологическими слоями и в явной форме выделить трещиноватые включения. Решена задача о распространении волн из очага землетрясения, происходящего на шельфе, получена полная трехмерная картина процесса. Проведено исследование прохождения сейсмических волн через трещиноватую среду с вариацией параметров трещиноватости. Оценена кинематическая и динамическая анизотропия регистрируемого на дневной поверхности сигнала.
Ключевые слова: сейсмические волны, математическое моделирование, сеточно-характеристический численный метод.
NUMERICAL SOLUTION OF TWO SEISMIC PROBLEMS WITH GRID-CHARACTERISTIC METHOD IN A FULL THREE-DIMENSIONAL SETTING
V. I. Golubev 1 2, I. B. Petrov 1 2
1 Moscow Institute of Physics and Technology, 2 System Research Institute, Russian Academy of Sciences, w.golubev@mail. ru, petrov@mipt. ru
Seismic prospecting is the routine method of locating oil and gas fields. In view of the discovery of non-traditional reservoirs, an important task is the development of new procedures for processing field data. The development of high-performance computing systems allows simulating seismic wave propagation processes in heterogeneous geological environments to prepare synthetic-seismograms with high precision. The article presents the use of the grid-characteristic method on hexahedral grids, which describes correctly the curved boundaries between geological layers and explicitly identifies fractured inclusions. The problem of waves' propagation from an earthquake source occurring in the shelf area is solved; a complete three-dimensional picture of the process is obtained. A study of the passage of seismic waves through a fractured medium with a variation in the fracture parameters is undertaken. The kinematic and dynamic anisotropy of the signal recorded on the day surface is estimated.
Keywords: seismic waves, mathematical modeling, grid-characteristic numerical method.
Введение. Решение сейсмических задач является актуальной проблемой ввиду их широкой распространенности. Они возникают в процессе проведения сейсмической разведки месторождений нефти и газа, при инициации сейсмической активности (землетрясений) в
сеисмоопасных регионах, проведении неразрушающего контроля элементов железнодорожной инфраструктуры, исследовании качества изготовления деталей из композитных материалов. Общим свойством всех перечисленных процессов является то, что они описываются в рамках динамического поведения гетерогенных сред. При этом могут быть использованы акустическое, упругое, упруго-пластическое, вязко-упругое приближения. К сожалению, актуальные постановки задач, включающие сложную геометрию расчетной области, существенную пространственную неоднородность физической среды, необходимость получения полноволнового решения не могут быть решены аналитически. При этом, вне зависимости от выбранной реологии определяющая система уравнений является гиперболической и должна быть решена численно с использованием адекватных численных методов.
За последние десятилетия было разработано множество разнообразных подходов к моделированию сейсмических процессов [1-2]. Наиболее популярными в вычислительной геофизике стали разновидности конечно-разностного (FDM) и спектрального методов (SEM), а также разрывный метода Галёркина (DGM) [6-8]. Параллельно с ними активно развивается сеточно-характеристический подход на структурированных [9-10] и неструктурированных сетках [11], который позволяет проводить явное выделение контактных поверхностей раздела сред и различных неоднородностей (в том числе трещиноватых кластеров).
В настоящей работе с использованием сеточно-характеристического метода на гекса-эдральных сетках решаются две исследовательские задачи. Первая посвящена распространению сейсмических волн из очага землетрясения на шельфе, модель которого включает в себя многослойный геологический массив и верхний пласт воды. Вторая заключается в исследовании процесса распространения сейсмических волн при проведении сейсмической разведки в трещиноватых средах, оценки влияния параметров трещиноватости на анизотропию регистрируемого на дневной поверхности сигнала.
Определяющие уравнения и численный метод. Рассмотрим подробнее подход, позволяющий проводить компьютерное моделирование процесса распространения сейсмических волн в неоднородных геологических средах. Динамическое состояние линейно-упругой среды описывается определяющей системой уравнений в частных производных вида
°ij = qijkiski •
P V <JH s,,
где н - плотность среды; i - компоненты вектора скорости смещения; 1 и kl - компоненты тензоров напряжений Коши и деформации; ^j - ковариантная производная по j -й
координате. Вид компонент тензора 4-го порядка qijkl определяется реологией среды. Для линейно-упругого изотропного случая они имеют вид
j = +M(ölköß + ôaô]k ).
g
В этом соотношении, которое обобщает закон Гука, ^ и ^ - параметры Ламе, а 1 -символ Кронекера. В настоящей работе численно решается данная определяющая система уравнений на гексаэдральных расчетных сетках сеточно-характеристическим методом. При этом на этапе построения сетки производится явное выделение границ неоднородностей -наборов узлов с каждого из бортов трещин. Поскольку решение уравнений акустики внутри трещины чрезвычайно ресурсоемко, используется подход с постановкой корректных граничных условий, обеспечивающий совпадение сейсмических откликов от трещиноватых объектов [12]. В узлах на противоположных бортах «флюидонасыщенной трещины» совпадают нормальные к плоскости трещины компоненты скорости, нормальные к плоскости трещины компоненты силы компенсируют друг друга, а тангенциальные компоненты силы равны нулю.
Распространение сейсмических волн из очага землетрясения. В работе исследовался процесс распространения упругих волн в геологическом массиве при инициации сей-
смической активности на шельфе. Модель среды была плоско-слоистой, включающей верхний слой жидкости, в полной трехмерной постановке [13]. Она состояла из пяти слоев с различными параметрами, соответствующими слагающим их породам. Первый слой описывался в рамках теории акустики и имел плотность 1 000 кг/м3, скорость распространения продольных волн 1 500 м/с, толщину 250 м. Все остальные слои считались упругими и описывались уравнениями линейной теории упругости. Второй слой (слой осадков) имел плотность 1 500 кг/м3, скорость распространения продольных волн 1 600 м/с, скорость распространения поперечных волн 60 м/с, толщину 50 м. Третий слой имел плотность 2 100 кг/м3, скорость распространения продольных волн 2 500 м/с, скорость распространения поперечных волн 1 000 м/с, толщину 300 м. Четвертый слой имел плотность 2 500 кг/м3, скорость распространения продольных волн 3 500 м/с, скорость распространения поперечных волн 1 300 м/с, толщину 400 м. Пятый слой (кристаллический фундамент) имел плотность 2 500 кг/м3, скорость распространения продольных волн 4 000 м/с, скорость распространения поперечных волн 2 500 м/с, толщину 500 м. В качестве источника сейсмического сигнала использовалась модель «подвижки по разлому» [14] с геометрическими размерами источника -50 х 150 х 150 м, углами наклона плоскости разрыва (rake=dip=slip=45 градусов) и глубиной залегания - 1 150 м.
Определяющая система уравнений в частных производных решалась численно сеточ-но-характеристическим методом с третьим порядком точности по пространству. Шаг расчетной сетки составлял 5 м, физическое время моделирования - 1,25 с. Использовался метод явного выделения контактных границ геологических слоев, что обеспечивает корректный расчет не только кинематических характеристик сейсмического сигнала, но и динамических характеристик (амплитуд) [15]. Необходимо отметить, что использованное контактное условие между первым слоем (жидкостью) и вторым слоем (слоем осадков) позволяет с высокой точностью воспроизвести наблюдения, осуществляемые при шельфовой сейсморазведке.
На рис. 1 приведены волновые картины, полученные в последовательные моменты времени. Продольные и поперечные волны от источника начинают распространяться к дневной поверхности и, проходя дно, порождают набор акустических волн разных частот и амплитуд.
Рис. 1. Распределение модуля скорости в шельфовой модели. Слева - 0,375 с (волны в третьем слое); справа - 0,875 с (волны в воде)
Распространение сейсмических волн в процессе сейсмической разведки.
Авторами была построена модель трещиноватой геологической среды. Она представляла собой трехмерный слоистый параллелепипед, содержащий флюидонасыщенные трещины. Размер модели составлял 10 500 х 10 500 х 3 900 м. Мощности слоев сверху вниз были 550 м, 200 м, 3 150 м. Первый и третий слой характе-ризовались параметрами, близкими к
параметрам геологических пород: скорость продольных волн - 4 500 м/с; скорость поперечных волн - 2 250 м/с; плотность - 2 500 кг/м3. Второй слой был контрастным и обладал другими характеристиками: скорость продольных волн - 6 000 м/с; скорость поперечных волн - 3 000 м/с; плотность - 2 500 кг/м3. На глубине 600 м во втором слое располагались вертикальные флюидонасыщенные трещины, задаваемые в виде прямоугольников с нормалью вдоль оси OX и размерами 100 м (вдоль OY и OZ). Расстояние вдоль OX и OY между ближайшими трещинами составляло 50 м. Трещиноватая среда покрывала в горизонтальной плоскости прямоугольник размерами 2 400 на 2 400 м. Описанная модель получила название verticaloriented. В дальнейшем, на основе описанной базовой модели была построена еще одна модель. В ней все трещины вертикальные, ориентированы случайно с равномерным распределением по азимуту, расстояние между ближайшими 150 м, названа vertical_random360. Также использовалась слоистая модель, полученная из базовой исключением трещиноватого блока, названная layered.
При проведении всех расчетов в работе использовалась площадная система наблюдений с расстоянием 10 м между ближайшими приемниками. Шаг по времени был выбран 0,8 мс, а всего моделировалось 1 440 шагов. В качестве источника использовалось точечное возбуждение, заглубленное на 50 м со временной зависимостью в виде импульса Рикера с частотой 35 Гц.
На рис. 2 представлены времена регистрации первого отклика для двух моделей: layered и verticaloriented. Как видно, они совпадают с высокой точностью. Для модели vertical_random360 времена прихода также не отличаются от приведенных, что свидетельствует о том, что структура трещиноватой среды не влияет на кинематические характеристики отклика.
150 200
Azimuth angle, qrad
Рис. 2. Времена регистрации сигнала в точках дневной поверхности. Зависимость от расстояния до источника и от азимута. Представлены результаты для моделей layered и vertical_oriented
В работе была изучена зависимость амплитуды зарегистрированного сигнала (максимум первого отклика) от азимута на приемник для всех описанных моделей. На рис. 3 представлены результаты для моделей vertical_oriented и vertical_random360. Для модели vertical_oriented наблюдается значимая зависимость амплитуды сигнала от азимута на уровне десятков процентов. В том случае, если ориентация трещин по азимуту хаотична, можно ожидать существенное снижение анизотропии сигнала, хотя некоторая ее часть и может остаться, что обусловлено наличием выделенного направления OZ даже в такой рандомизированной среде. На рис. 3 кривые для vertical_random360 значительно ближе к горизонтальной для всех рассматриваемых удалений, что подтверждает высказанные предположения.
Azimutti anqlc, qrad
Рис. 3. Амплитуды сигнала, зарегистрированного на дневной поверхности на различных удалениях от источника. Сравнение результатов для моделей vertical_oriented и vertical_random360
Заключение. В работе проведен численный расчет распространения сейсмических волн из очага землетрясения в геологической модели шельфа России. Показана возможность определения волновых полей как в многослойной упругой среде, так и в водном слое при использовании корректных контактных условий на границе акустическое-упругое тела. Проведена серия расчетов сейсмического отклика от трещиноватой геологической среды в полной трехмерной постановке с изменяющейся внутренней структурой. Характерные размеры трещиноватых моделей 10 х 10 х 4 км, частота сигнала источника 35 Гц, скоростные характеристики слоев 4 500-6 000 м/с, протяженность трещин 100 м с интервалом между ближайшими 50 м. Система наблюдений выбиралась так, чтобы обеспечить регистрацию отраженных волн с углами отражения вплоть до критических. Анализ результатов расчетов подтвердил, что изученные варианты структуры трещиноватой среды не влияют на время прихода сейсмического отклика. Получена значимая анизотропия регистрируемого сигнала при наличии ориентированной структуры в среде, которая значительно снижается при хаотизации среды по азимуту. Во всех проведенных расчетах использовался сеточно-характеристический метод на гексаэдральных расчетных сетках.
Работа была выполнена с использованием оборудования центра коллективного пользования «Комплекс моделирования и обработки данных исследовательских установок мега-класса» НИЦ «Курчатовский институт», http://ckp.nrcki.ru/.
Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 16-29-02018 офи_м.
Литература
1. Virieux J., Calandra H., Plessix R. E. A review of the spectral, pseudo-spectral, finite-difference and finite-element modelling techniques for geophysical imaging // Geophys Prosp. 2011. V. 59. № 5. P. 794-813.
2. Carcione J. M., Herman G. C., Kroode P. E. Seismic modeling // Review Literature And Arts Of The Americas. 2002. V. 67. № 4. P. 1304-1325.
3. Virieux J., Etienne V. et al. Modeling seismic wave propagation for geophysical imaging // Seismic Waves: Research and Analysis. InTech. 2012. Р. 253-304.
4. Moczo P., Kristek J., Galis M. The Finite-Difference Modeling of Earthquake Motions: Waves and Ruptures. Cambridge : Cambridge Univ. Press, 2014.
5. LeVeque R. J. Finite Volume Methods for Hyperbolic Problems. New York : Cambridge University Press, 2002.
6. Etienne V., Chaljub E., Virieux J. et al. An hp-adaptive discontinuous Galerkin finite-element method for 3D elastic wave modelling // Geophys J Int. 2010. Т. 183. № 2. P. 941-962.
7. Hermann V., Käser M., Castro C. E. Non-conforming hybrid meshes for efficient 2D wave propagation using the Discontinuous Galerkin Method // Geophysical Journal International. 2011. Т. 184. № 2. С. 746-758.
8. Mercerat E. D., Glinsky N. A nodal high order discontinuous Galerkin method for elastic wave propagation in arbitrary heterogeneous media // Geophys J Int. 2015. № 201 (2). Р. 10991116.
9. Голубев В. И., Петров И. Б., Хохлов Н. И. Моделирование волновых процессов внутри планеты с помощью гибридного сеточно-характеристического метода // Матем. моделирование. 2015. Т. 27. № 2. С. 139-148.
10. Голубев В. И., Петров И. Б., Хохлов Н. И., Шульц К. И. Численный расчет волновых процессов в трещиноватых средах на гексаэдральных сетках сеточно-характеристическим методом // Журн. вычисл. математики и матем. Физики. 2015. Т. 55. № 3. С. 512-522.
11. Chelnokov F. B. Explicit representation of grid-characteristic schemes for elasticity equations in 2D and 3D spaces // Mat Model. 2016. № 18 (6). Р. 96-108.
12. Левянт В. Б., Миряха В. А., Муратов М. В., Петров И. Б. Оценка влияния на сейсмический отклик степени раскрытости трещины и доли площади локальных контактов к ее поверхности // Технологии сейсморазведки. 2015. № 3. С. 16-30.
13. Заславский Ю. М., Кержаков Б. В., Кулинич В. В. Вертикальное сейсмическое профилирование на морском шельфе // Акустич. журн. 2008. Т. 54. № 3. С. 483-490.
14. Golubev V. I., Kvasov I. E., Petrov I. B. Influence of natural disasters on ground facilities // Mathematical Models and Computer Simulations. 2012. № 4 (2). Р. 129-134.
15. Голубев В. И., Петров И. Б. Опыт расчета сейсмических откликов от криволинейных геологических границ на основе их явного выделения // Технологии сейсморазведки. 2016. № 4. С. 45-50.
