CC BY
5УДК 519.6 Научная статья
https://doi.org/10.23947/2587-8999-2023-6-1-27-33
Прямое численное моделирование: топография дневной поверхности и анизотропия верхней части разреза
В. И. Голубев 1,2 О И, А. В. Шевченко 1,2 , А. В. Екименко3 , В. Ю. Петрухин4
'Московский физико-технический институт (национальный исследовательский университет), Российская Федерация, г. Долгопрудный, Институтский переулок, 9
^Институт автоматизации проектирования Российской Академии Наук, Российская Федерация, г. Москва, ул. 2-я Брестская, 19/18
3ООО «Газпромнефть НТЦ», Российская Федерация, г. Санкт-Петербург, наб. реки Мойки, 75-79д 4ПАО Сбербанк, г. Москва, ул. Вавилова, 19
И w. golubev@mail.ru
Аннотация
Введение. Одной из проблем разработки нефтегазовых месторождений является построение корректных геологических моделей подповерхностного пространства. Исследователями из различных научных групп во всем мире предложены различные способы повышения точности компьютерного моделирования, используемого в этом процессе. Цель настоящего исследования — оценка степени влияния рельефа дневной поверхности и анизотропии верхней части геологического разреза на регистрируемый сейсмический сигнал на примере реалистичной модели Оренбургского месторождения.
Материалы и методы. Рассмотрена сейсмогеологическая модель, описывающая нижнепермский интервал Оренбургского геологического разреза. По скважинным данным оценены упругие свойства геологических пластов: плотность и скорости распространения продольных и поперечных волн. Отмечается высокий контраст скоростей продольных волн, оцененных по диаграммам акустического каротажа. Резервуар в этой модели приурочен к нижним пластам. Он сложен сульфатно-карбонатными средами, однородными по плотности и акустическим свойствам. С использованием сеточно-характеристического метода были рассчитаны синтетические сейсмограммы нулевых удалений. Выбор структурных криволинейных расчетных сеток позволил корректно провести учет рельефа дневной поверхности.
Результаты исследования. В работе было проведено сопоставление двух различных моделей. В первую модель была включена анизотропия верхней части разреза и топография дневной поверхности. Во второй модели верхняя граница расчетной области была плоской, и вся среда рассматривалась в рамках изотропной линейно упругой модели. Анализ синтетических сейсмограмм показал, что анизотропия, присущая этой модели, не оказывает существенного влияния на регистрируемое сейсмическое волновое поле. Однако учет рельефа дневной поверхности заметно сдвигает времена прихода отраженных волн.
Обсуждение и заключение. Представленный в работе алгоритм может быть использован для верификации графа обработки полевых данных, поскольку оценка анизотропии среды является стандартным шагом при построении скоростной модели. Представленный подход может быть расширен на трехмерные модели реалистичных размеров.
Ключевые слова: математическое моделирование, сейсмическая разведка, сеточно-характеристический численный метод, анизотропные среды, топография дневной поверхности
Финансирование. Исследование выполнено за счет гранта Российского научного фонда № 19-71-10060, https://rscf.ru/project/19-71-10060/.
Для цитирования. Прямое численное моделирование: топография дневной поверхности и анизотропия верхней части разреза / В. И. Голубев, А. В. Шевченко, А. В. Екименко, В. Ю. Петрухин // Computational Mathematics and Information Technologies. — 2023. — Т. 6, № 1. — С. 27-33. https://doi.org/10.23947/2587-8999-2023-6-1-27-33
© Голубев В. И., Шевченко А. В., Екименко А. В., Петрухин В. Ю., 2023
27
Original article
Direct seismic modeling: day surface topography and shallow subsurface anisotropy V. I. Golubev12 , A. V. Shevchenko12 , A. V. Ekimenko3 , V. Yu. Petrukhin4
1Moscow Institute of Physics and Technology (National Research University), 9, Institutskii lane, Dolgoprudny, Russian Federation institute of Computer Aided Design of RAS, 19/18, 2-nd Brestskaya St., Moscow, Russian Federation 3Gazprom Neft Science and Technology Centre, 75-79d, Moika emb. river, St. Petersburg, Russian Federation 4PAO Sberbank, 19, Vavilova St., Moscow, Russian Federation
H w.golubcv@mail.ru Abstract
Introduction. The article is devoted to one of the problems in the development of oil and gas fields—the construction of correct geological models of the subsurface space. Researchers from various scientific groups around the world have proposed various ways to improve the accuracy of the computer simulations used in this process. The purpose of this study is to assess the degree of influence of the day surface relief and the anisotropy of the upper part of the geological section on the recorded seismic signal using a realistic model of the Orenburg field as an example.
Materials and methods. A seismogeological model describing the Lower Permian interval of the Orenburg geological section is considered. According to well data, the elastic properties of geological formations were estimated: density and propagation velocities of longitudinal and transverse waves. There is a high contrast of P-wave velocities estimated from sonic logs. The reservoir in this model is confined to the lower layers. It is composed of sulfate-carbonate media, uniform in density and acoustic properties. Using the grid-characteristic method, zero-offset synthetic seismograms were calculated. The choice of structural curvilinear computational grids made it possible to correctly consider the relief of the day surface. Research results. In this work, two different models were compared. The first model included the anisotropy of the upper part of the section and the topography of the day surface. In the second model, the upper boundary of the computational domain was flat, and the entire medium was considered within the framework of an isotropic linear elastic model. The analysis of synthetic seismograms showed that the anisotropy inherent in this model does not significantly affect the recorded seismic wave field. However, considering the relief of the day surface significantly shifts the times of arrival of reflected waves. Discussion and Conclusion. The algorithm presented in the paper can be used to verify the field data processing graph, since the assessment of the anisotropy of the medium is a standard step in building a velocity model. The presented approach can be extended to 3D models of realistic dimensions.
Keywords: mathematical modeling, seismic exploration, grid-characteristic numerical method, anisotropic media, day surface topography.
Funding information. The research was carried out at the expense of the grant of the Russian Science Foundation No. 1971-10060, https://rscf.ru/project/19-71-10060/.
For citation. Direct seismic modeling: day surface topography and shallow subsurface anisotropy / V. I. Golubev, A. V. Shevchenko, A. V Ekimenko, V Yu. Petrukhin // Computational Mathematics and Information Technologies. — 2023. — Vol. 6, no. 1. — P. 27-33. https://doi.org/10.23947/2587-8999-2023-6-1-27-33
Введение. Проблема учета в вычислительной модели топографии дневной поверхности при проведении сейсмического моделирования исследуется многими учеными. Простейший подход основан на ступенчатом описании границы и конечно-разностных схемах [1] на прямоугольных сетках. Однако он приводит к появлению дополнительных артефактов-источников из-за дифракции. Метод повернутых сдвинутых сеток был предложен для сопоставления сетки с границами области [2]. Использование конечно-разностных схем на деформированных сетках не позволяет описывать очень крутые наклоны в рельефе [3]. В работе [4] был предложен псевдоспектральный метод с отображением прямоугольных сеток на криволинейную поверхность, который также оказался неустойчив вблизи сильно наклонных границ. По-видимому, приемлемым решением является использование конечно-элементных методов на неструктурных сетках [5]. Недостатком данного подхода является значительное возрастание вычислительной сложности задачи и требуемой оперативной памяти.
В работе [6] был предложен подход по комбинированию разрывного метода Галёркина с методом конечных разностей. Он позволяет использовать неструктурные сетки только вблизи поверхности, покрывая основную часть геологического массива структурной сеткой.
Для расчета волновых полей в упругих средах с резким изменением упругих параметров может быть использован сеточно-характеристический метод [7, 8]. С ним был объединен метод наложенных сеток, позволяющий корректно описывать криволинейные границы [9]. В настоящем исследовании используется другой подход. Отдельная криволинейная сетка используется для покрытия каждого геологического слоя, совместно с явной постановкой контактных условий между отдельными слоями. Данный метод сохраняет преимущества сеточно-характеристического метода, не приводя к чрезмерному возрастанию вычислительной сложности задачи. В работе учитывается тонкослоистость верхней части геологического разреза с применением широко используемой УТ1 модели среды [10, 11].
Целью данной работы являлась оценка степени влияния рельефа дневной поверхности и анизотропии верхней части геологического разреза на регистрируемый сейсмический сигнал на примере реалистичной модели Оренбургского месторождения.
Материалы и методы. Значительный объем сейсморазведки связан с изучением нефтегазоносных бассейнов. Эти территории, сложенные осадочными породами, представляют собой мощные толщи тонкослойного переслаивания. Условия образования осадочных пород определяют различие физических свойств каждого слоя. Геологическая история каждого региона определяет современное положение сформировавшихся пластов горных пород. Поэтому характерным признаком является наличие пластов, вытянутых в поперечном направлении, но маломощных по отношению к используемым длинам сейсмических волн. Это дает возможность описать данную геологическую среду моделью УТ1. Она была сформулирована в связи с регистрацией разных скоростей для вертикальных и наклонных сейсмических лучей в 1930-х годах. В работе рассмотрена сейсмогеологическая модель, описывающая нижнепермский интервал Оренбургского геологического разреза. Физические свойства оценивались по скважинным данным Оренбургского нефтегазоконденсатного месторождения (рис. 1 а). Выделены две основные части модели: 1) тонкослоистый верхний интервал, соответствующий кунгурскому ярусу, 2) нижний интервал, описывающий филипповский горизонт и артинский ярус. Первая часть представлена чистыми солями и прослоями каменной соли, ангидритов и доломитов. Отмечается высокий контраст скоростей продольных волн, оцененных по диаграммам акустического каротажа. Типичное значение для каменной соли составило 4500 м/с, для ангидритов — 6300 м/с, а для доломитов — в диапазоне от 4500 м/с до 6500 м/с (рис. 1 б). Этот интервал влияет на кинематические характеристики волнового поля и вызывает анизотропию.
На рис. 2 представлены полевые сейсмические данные. Чередование доломитовых, солевых и ангидритовых слоев приводит к наличию динамически выраженных отражений. Разрез показывает сложную геометрию отражающих границ верхнего интервала и процессы галокинеза, сопровождающиеся проникновением солей в вышележащие терригенные толщи. Можно предположить, что в областях расположения контрастного субгоризонтального слоя параметры анизотропии будут иными, чем у соляных диапиров. Резервуар в этой модели приурочен к нижним пластам. Он сложен сульфатно-карбонатными средами, однородными по плотности и акустическим свойствам. Отражающие границы в этом интервале малоконтрастны. Этот факт объясняется отсутствием акустически жестких границ. По данным геофизических исследований скважин интервал относительно однороден (рис. 1 а). На основе приведенного выше описания была разработана упругая двумерная модель. Для достижения требуемой точности моделью УТ1 был описан тонкослоистый интервал.
Артинский ярус-Р, аг
Кунгурский ярус-Р, к
Уфа + Казан ь-Р2 u + Р2 kz
Филипповский горизонт-Р,fl
Иренский горизонт-Р, ir
Горизонт
Глубина, м
кш
iOLf
J4 С
Изотропная среда
4-1---*■*-
Анизотропная тонкослойная среда VTI
а)
1200 1000 800 600 400 200 0
4000 4500 5000 5500 6000 6500 4000 4500 5000 5500 6000 6500 4000 4500 5000 5500 6000 6500 Cp, м/с Cp, м/с Cp, м/с
б)
Рис. 1. Результаты натурных измерений: а — скважинные данные с Оренбургского газоконденсатного месторождения; б—гистограммы распределения скорости продольных волн по акустическим данным
Следы 380 420 460 500 540 580 620 660
Рис. 2. Сейсмический разрез Оренбургского месторождения
Определяющая система уравнений в частных производных гиперболического типа имеет вид:
Я, + АЧх + АтЯу = ?•
где вектор о, состоит из всех неизвестных задачи (компоненты вектора скорости и тензора напряжений). Значения элементов матриц А , А2 определяются параметрами геологической среды. Для решения исходной системы используется процедура расщепления по пространственным направлениям. Сначала решается одномерная система О, + А1ох = 0, потом решается о, + А2оу = 0 и после — О, = /■ Таким образом, многомерная задача сводится к набору одномерных задач. Учет правой части может быть эффективно реализован, поскольку он выполняется независимо в каждом узле расчетной сетки. Каждая решаемая одномерная система имеет вид:
0 + АО = 0
Гиперболичность системы означает, что матрица А имеет полный набор собственных векторов и может быть представлена в виде:
А = ОГ'ЛО.,
где матрица Л диагональная и состоит из собственных значений матрицы А. Вводя в рассмотрение переход в инварианты м> = ООО, получаем набор независимых линейных уравнений переноса вида: + Лй^ = 0 . Для каждого из них справедливо характеристическое соотношение:
й{х, + т) = й{х Т,)
где нижний индекс обозначает соответствующую компоненты вектора. Обратный переход к переменным задачи осуществляется невырожденным преобразованием: 0 = ^
Критически важным отличием при реализации данного метода на прямоугольных и криволинейных сетках является тот факт, что в первом случае существуют лишь два выделенных направления, вдоль которых матрицы О могут быть предвычислены автоматически. В случае криволинейных сеток в каждой точке строится локальная система координат (4, п), вдоль осей которой вычисляются матрица О и обратная к ней с помощью численной итерационной процедуры.
Результаты исследования. Было проведено распространение упругих волн в описанной геологической модели. Использовалась структурная сетка, построенная с использованием программного обеспечения иМЛМЛЬЬЛ [12] и покрывающая область размерами 10 км х 5 км. Отдельные одномерные уравнения переноса решались по схеме Русанова. Пространственный шаг сетки составлял порядка 2 м. Шаг по времени выбирался из условия устойчивости Куранта. Рассматривалось две постановки задачи. Первая — включающая анизотропию и топографию, вторая — не включающая их. Результаты расчетов представлены на рис. 3. Тщательный анализ синтетических данных показал, что анизотропия, присущая этой модели, не оказывает существенного влияния на сейсмическое волновое поле. Однако рельеф заметно сдвигает времена прихода отраженных волн.
Время, с 0
Рис. 3. Синтетические сейсмограммы: а — изотропной модели с плоской дневной поверхностью;
б—анизотропной модели с топографией
Обсуждение и заключения. Проведено исследование влияния топографии дневной поверхности и анизотропии верхней части разреза на данные сейсмической разведки. Рассматриваемая геологическая модель основывалась на реальном месторождении нефти и газа (Оренбургское месторождение). Сеточно-характеристический метод на структурированных сетках был успешно применен для учета анизотропии и криволинейных границ в модели. Рассчитанное волновое поле содержит основные особенности исходного сейсмического разреза, что подтверждает адекватность построенной модели. Результаты проведенных расчетов свидетельствуют о значительном влиянии рассматриваемых факторов на регистрируемый на дневной поверхности сигнал.
Представленный в работе алгоритм может быть использован для верификации графа обработки полевых данных, поскольку оценка VTI анизотропии является стандартным шагом при построении скоростной модели [13]. Представленный подход может быть расширен на трехмерные модели реалистичных размеров.
Список литературы
1. Virieux J. P-SV wave propagation in heterogeneous media: Velocity — stress finite-difference method / J. Virieux // Geophysics. — 1986. — Vol. 51, no. 4. — P. 889-901.
2. Saenger, E. H. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid / E. H. and T. B. Saenger // Geophysics. — 2004. — Vol. 69, no. 2. — P. 583-591.
3. Tarrass, I. New curvilinear scheme for elastic wave propagation in presence of curved topography / I. Tarrass, L. Giraud and P. Thore // Geophys. — 2011. — Vol. 59, no. 5. — P. 889-906.
4. Tessmer, E. Elastic wave propagation simulation in the presence of surface topography / E. Tessmer, D. and A. B. Kosloff // Geophys. — 1992. — Vol. 108, no. 2. — P. 621-632.
5. An hpadaptive discontinuous Galerkin finite-element method for 3D elastic wave modelling / V. Etienne, E. Chaljub, J. Virieux and N. Glinsky // Geophysical Journal International. — 2010. — Vol. 183.— P. 941-962.
6. Lisitsa, V. Combination of the discontinuous Galerkin method with finite differences for simulation of seismic wave propagation / V Lisitsa, V Tcheverda and C. Botter // Journal of Computational Physics. — 2016. — Vol. 311.— 142-157.
7. Compact Grid-Characteristic Scheme for the Acoustic System with the Piece-Wise Constant Coefficients / V. Golubev, A. Shevchenko, N. Khokhlov, I. Petrov and M. Malovichko // International Journal of Applied Mechanics. — 2022. — Vol. 14, no. 2.
8. Golubev, V. Raising convergence order of grid-characteristic schemes for 2D linear elasticity problems using operator splitting / V. Golubev, A. Shevchenko and I. Petrov // Computer Research and Modeling — 2022. — Vol. 14, no. 4. — P. 899-910.
9. Khokhlov, N. I. Overset grids approach for topography modeling in elastic-wave modeling using the grid-characteristic method. / N. I. Khokhlov, V. O. Stetsyuk and I. A. Mitskovets // Computer Research and Modeling. — 2019. — Vol. 11, no. 6. — P. 1049-1059.
10. Petrov, I. B. Simulation of Seismic Waves in Anisotropic Media / I. B.Petrov, V. I. Golubev, V. Y. Petrukhin and I. S. Nikitin // Mathematics. — 2021. — Vol. 103, no. 3. — P. 146-150.
11. Golubev, V. Explicit simulation of seismic waves in fractured VTI media. / V. Golubev and A. Shevchenko // 82nd EAGE Annual Conference and Exhibition. — 2021.
12. Generating Quality Structured Convex Grids on Irregular Regions / P. Barrera, F. Dominguez, G. F. Gonzalez G. and Tinoco // Electronic Transactions on Numerical Analysis. — 2009. — Vol. 34.— P. 76-89.
13. Deep Convolutional Neural Networks in Seismic Exploration Problems / A. Vasyukov, I. Nikitin, A. Stankevich and V. Golubev // Interfacial Phenomena and Heat Transfer. — 2022. — Vol. 10, no. 3. — P. 61-74.
Поступила в редакцию 07.02.2023. Поступила после рецензирования 02.03.2023. Принята к публикации 03.03.2023.
Об авторах:
Голубев Василий Иванович, ведущий научный сотрудник лаборатории прикладной вычислительной геофизики Московского физико-технического института (национального исследовательского университета) (РФ, 141701, Московская область, г. Долгопрудный, Институтский переулок, 9), доктор физико-математических наук, доцент, ORCID, w.golubev@ mail.ru
Шевченко Алексей Владимирович, аспирант Московского физико-технического института (национального исследовательского университета) (РФ, 141701, Московская область, г. Долгопрудный, Институтский переулок, 9), ORCID, alexshevchenko@phystech.edu
Екименко Антон Валерьевич, эксперт блока интегрированных решений ООО «Газпромнефть НТЦ» (РФ, 190000, г. Санкт-Петербург, наб. реки Мойки, 75-79д), кандидат геолого-минералогических наук, ORCID, ekimenko.av@gazpromneft-ntc.ru
Петрухин Вячеслав Юрьевич, руководитель направления ПАО Сбербанк (РФ, 117312, г. Москва, ул. Вавилова, 19), ORCГО, v.y.petrukhin@gmail.com
Заявленный вклад соавторов:
В. И. Голубев — постановка цели и задачи исследования, подготовка текста, формулирование выводов. А. В. Шевченко — проведение численных расчетов по изотропной модели без рельефа. А. В. Екименко — построение геологической модели, анализ результатов численных расчетов. В. Ю. Петрухин — проведение расчетов по анизотропной модели с рельефом.
Конфликт интересов.
Авторы заявляют об отсутствии конфликта интересов.
Все авторы прочитали и одобрили окончательный вариант рукописи.
