Научная статья на тему 'Исследование сейсмических волновых полей в зоне залегания коллекторов при помощи численного моделирования полей методами непродольного вертикального сейсмического профилирования и общей точки возбуждения'

Исследование сейсмических волновых полей в зоне залегания коллекторов при помощи численного моделирования полей методами непродольного вертикального сейсмического профилирования и общей точки возбуждения Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

CC BY
53
12
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ / СЕЙСМИЧЕСКИЕ ПОЛЯ / КОЛЛЕКТОР / СЕЙСМИЧЕСКОЕ ЭКРАНИРОВАНИЕ / СЛОИ БИО / NUMERICAL MODELING / SEISMIC FIELDS / COLLECTOR / SCREENING EFFECT / BIOT LAYERS

Аннотация научной статьи по наукам о Земле и смежным экологическим наукам, автор научной работы — Ковтун Александр Андреевич, Лыскова Евгения Леонидовна, Санников Константин Юрьевич

В первой части работы рассмотрена частная модель двухслойной высокоскоростной толщи, состоящей из покрышки и коллектора, залегающей в низкоскоростной вмещающей среде. На основе результатов интерпретации модельных сейсмических волновых полей методами непродольного вертикального сейсмического профилирования (НВСП) и общей точки возбуждения (ОТВ) установлено, что целевая граница кровля коллектора может быть выделена в случае достаточно большой мощности высокоскоростной покрышки по полю отраженных обменных волн типа P1S2S2P1 и P1S2S21, порождаемых проходящей обменной поперечной волной, образующейся на кровле покрышки и распространяющейся вглубь среды. Волны P1S2S2P1 и P1S2S21 существуют только в определенном интервале эпицентральных удалений от источника, протяженность которого зависит от контрастности высокоскоростной покрышки, а также ее мощности и глубины залегания. Надежное определение границы между высокоскоростной покрышкой и коллектором может быть затруднено при использовании только наземных сейсмических данных вследствие проявления экранирования. Дополнительную информацию можно получить при привлечении материалов НВСП и корректном использовании данных отраженных обменных волн типа P1S2S2P1 и P1S2S21. Во второй части работы результаты численных модельных исследований полноволнового сейсмического отклика от тонкослоистой пачки насыщенных упруго-пористых слоев Био, моделирующей резервуарную зону, показывают, что низкочастотные амплитудные аномалии, наблюдаемые на сейсмических разрезах ниже резервуарной зоны, могут быть вызваны образованием поля отраженных обменных поперечных волн, запаздывающих относительно отклика продольной волны и имеющих на малых и средних офсетах значимые амплитуды на вертикальной компоненте поля. Амплитуды и временные задержки суммарной обменной поперечной волны могут варьироваться в широких пределах в зависимости от конкретных условий: глубины залегания и мощности резервуара, частоты падающей продольной волны, наличия в резервуарной зоне границ с контактом проскальзывания или протяженных тонких трещин, заполненных флюидом, и других факторов.This paper presents a study of seismic wave fields in reservoir zones by numerical simulation. In the first part of the study, numerical modeling of full wavefield was performed to synthesize the multi-offset vertical profiling data and data of OPV (Common Source Point Technique) for the model of the medium in which high-velocity layers (seal and collector) are placed in a low-velocity containing medium. It has been established that the principal boundary, the roof of the collector, can be distinguished by data only in the case of great thickness of high-velocity seal from the field of converted P1S2S2P1 and P1S2S21 waves. These waves are formed on the roof of the seal during the passage of a converted transverse wave and exist only in a certain interval of epicentral distances from the source, the length of which depends on the contrast of the high-velocity seal, as well as on its thickness and depth. Reliable determination of the boundary between the high-velocity seal and the collector can be difficult when using only ground-based seismic data due to the screening effect. However, it becomes possible with the use of multi-offset vertical profiling materials and the correct use of the data of reflected converted P1S2S2P1 and P1S2S21 waves. In the second part of the work, we analyze the results of numerical model studies of the full-wave seismic response from a thin-layered pack of saturated elastic porous Biot layers simulating a reservoir zone. The results of numerical model studies show that the low-frequency amplitude anomalies observed in seismic sections below the reservoir zone may be caused by the formation of a field of reflected converted transverse waves that are delayed relative to the response of the longitudinal wave and having significant amplitudes on the vertical field component for the small and middle offsets. The amplitudes and time delays of the total converted transverse wave can vary widely depending on specific conditions: depth and reservoir thickness, incident longitudinal wave frequency, the presence of boundaries in the reservoir area with sliding contact or extended thin cracks filled with fluid, and other factors.

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по наукам о Земле и смежным экологическим наукам , автор научной работы — Ковтун Александр Андреевич, Лыскова Евгения Леонидовна, Санников Константин Юрьевич

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Исследование сейсмических волновых полей в зоне залегания коллекторов при помощи численного моделирования полей методами непродольного вертикального сейсмического профилирования и общей точки возбуждения»

УДК 550.834

Вестник СПбГУ. Науки о Земле. 2020. Т. 65. Вып. 1

Исследование сейсмических волновых полей в зоне залегания коллекторов при помощи численного моделирования полей методами непродольного вертикального сейсмического профилирования и общей точки возбуждения А. А. Ковтун, Е. Л. Лыскова, К. Ю. Санников

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7-9

Для цитирования: Ковтун, А. А., Лыскова, Е. Л., Санников, К. Ю. (2020). Исследование сейсмических волновых полей в зоне залегания коллекторов при помощи численного моделирования полей методами непродольного вертикального сейсмического профилирования и общей точки возбуждения. Вестник Санкт-Петербургского университета. Науки о Земле, 65 (1), 33-50. https://doi.org/10.21638/spbu07.2020.103

В первой части работы рассмотрена частная модель двухслойной высокоскоростной толщи, состоящей из покрышки и коллектора, залегающей в низкоскоростной вмещающей среде. На основе результатов интерпретации модельных сейсмических волновых полей методами непродольного вертикального сейсмического профилирования (НВСП) и общей точки возбуждения (ОТВ) установлено, что целевая граница — кровля коллектора — может быть выделена в случае достаточно большой мощности высокоскоростной покрышки по полю отраженных обменных волн типа Р18282Р1 и Р182Б21, порождаемых проходящей обменной поперечной волной, образующейся на кровле покрышки и распространяющейся вглубь среды. Волны Р182Б2Р1 и Р^2821 существуют только в определенном интервале эпицентральных удалений от источника, протяженность которого зависит от контрастности высокоскоростной покрышки, а также ее мощности и глубины залегания. Надежное определение границы между высокоскоростной покрышкой и коллектором может быть затруднено при использовании только наземных сейсмических данных вследствие проявления экранирования. Дополнительную информацию можно получить при привлечении материалов НВСП и корректном использовании данных отраженных обменных волн типа Р18282Р1 и Р182Б21. Во второй части работы результаты численных модельных исследований полноволнового сейсмического отклика от тонкослоистой пачки насыщенных упруго-пористых слоев Био, моделирующей резервуарную зону, показывают, что низкочастотные амплитудные аномалии, наблюдаемые на сейсмических разрезах ниже резервуарной зоны, могут быть вызваны образованием поля отраженных обменных поперечных волн, запаздывающих относительно отклика продольной волны и имеющих на малых и средних офсетах значимые амплитуды на вертикальной компоненте поля. Амплитуды и временные задержки суммарной обменной поперечной волны могут варьироваться в широких пределах в зависимости от конкретных условий: глубины залегания и мощности резервуара, частоты падающей продольной волны, наличия в резервуарной зоне границ с контактом проскальзывания или протяженных тонких трещин, заполненных флюидом, и других факторов.

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

© Санкт-Петербургский государственный университет, 2020

1. Введение

Задачи, связанные с выделением и изучением коллекторов по сейсмическим данным, не перестают быть актуальными. В ряде работ изучаются особенности полей отраженных волн для среды, содержащей карбонатный слой-коллектор и высокоскоростную покрышку. Так, в работе (Левянт и др., 2010) отмечены сложности в определении границы между высокоскоростной покрышкой и основным телом карбонатного резервуара по данным наземной сейсмики. На большом количестве экспериментального материала показано, что кровля карбонатного резервуара однозначно не прослеживается на сейсмических разрезах общей глубинной точки и может быть определена в основном только по данным акустического каротажа, по уменьшению акустического импеданса. Одной из возможных причин отсутствия отражений от кровли резервуара может быть не только шероховатость этой границы (Левянт и др., 2010), но и явление сейсмического экранирования. В работе (Голикова и др., 2016а; 2016б) при помощи численного моделирования исследовано влияние высокоскоростного слоя на процесс волнообразования. Показано, что состав поля проходящих и отраженных волн и их свойства могут существенно изменяться в зависимости от эпицентрального удаления и мощности слоя.

Следующая особенность волновых полей коллекторов связана с тем, что ниже пласта-коллектора могут возникать амплитудные аномалии продольных волн, имеющие значительные временные задержки (до 200 мс) относительно отражений от кровли резервуарной зоны (Chabyshova and Goloshubin, 2014; Ahmad et al., 2017). В работе (Chabyshova and Goloshubin, 2014) проверяется гипотеза, согласно которой амплитудные аномалии и их временные задержки могут быть вызваны конверсионными переходами типа «быстрая-медленная-быстрая» Р-волна, генерируемыми в тонкослоистых пористых проницаемых средах, насыщенных флюидами. В рамках асимптотической теории пороупругости (Silin and Goloshubin, 2010; Голошубин и др., 2015) авторы приводят расчеты откликов от модельного тонкослоистого многослойного резервуара в случае нормального плосковолнового распространения продольных волн. Однако полученные результаты не дают однозначного ответа на поставленные вопросы.

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

Численное моделирование волновых полей в обеих задачах проводилось при помощи программы OASP из комплекса OASES (Shmidt, 2004). Эта программа реализует двумерное вычисление полного волнового поля, возбуждаемого нестацио-

нарным сосредоточенным источником, в слоисто-однородной осесимметрической упругой среде с плоско-параллельными границами раздела на основе использования специальных методик численного интегрирования в частотной области и применения метода глобальной матрицы (S^midt and Jensen, 1985). Используемый метод вычислений позволяет учесть все возможные обмены между быстрой и медленной продольными волнами, а также обмены между продольными и поперечными волнами на всех границах раздела модельной среды.

2. Модель коллектора с высокоскоростной покрышкой 2.1. Описание модели

Многочисленные исследования нефтегазоносных провинций указывают на присутствие во вмещающих толщах (осадочных, низкоскоростных породах) высокоскоростных геологических образований (соляных куполов, пластов плотных карбонатов, интрузивных и эффузивных тел). В этой связи возникает необходимость в правильной интерпретации сейсмических волновых полей, проходящих через высокоскоростные слои. В частности, подобная проблема возникает при интерпретации сейсморазведочных данных в некоторых районах северных морей и Ледовитого океана, связанная с залегающими на разных глубинах высокоскоростными эффузивными покровами, которые оказывают влияние на поля проходящих и отраженных волн (Bruvoll et al., 2012).

Рассмотрим модель коллектора с высокоскоростной покрышкой. Удобно использовать следующую номенклатуру волн. Прямую падающую на первый слой (покрышку) продольную волну обозначим P1, отраженную от кровли слоя — как P1P1. Монотипной волне, проходящей через слои 1-2-3, присвоим обозначение P123. Таким образом, нижний индекс у волны обозначает номера слоев, через которые она последовательно проходит. Прописная буква в обозначении поля проходящих и отраженных волн будет обозначать смену типа волны (проходящее поле P1S234) или изменение направления распространения — отраженные поля P1P1, P123P3 (рис. 1).

R

Источник pi<3,soi p,!iMq„, РюэР321

Vp1=2.0 км/с VS1 =1.1 км/С Рт =1.9 г/см3

X

Vp2=5.0 км/с Vs2=2.6 км/с Р2=2.5 г/см3

Vp3 =4.3 км/с Vs3=2.2 км/с р3=2.4г/см3

Vp4=3.0 км/с Vs4=1 -5 км/с Р4=2.3 г/см3

Z

Рис. 1. Схема наблюдений, модель среды и номенклатура волн

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

Параметры базовой контрастной модели М-1 имели следующие значения:

1) покрывающее полупространство (г < 0): у^ = 2.0 км/с, у^ = 1.1 км/с, р: = 1.9 г/см3;

2) слой-покрышка (0 < г < ^ = 100 м: у^ = 5.0 км/с, у^ = 2.6 км/с, р2 = 2.5 г/см3;

3) слой-коллектор ф2 < г < ^ + ^ =50 м: у^=4.3 км/с, у^ = 2.2 км/с, р3 = 2.4 г/см3;

4) нижнее полупространство (г > ^+у^ = 3.0 км/с, у^ = 1.5 км/с, р4 = 2.3 г/ см3.

Дневная поверхность в модели отсутствует, поскольку ее введение приводит

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

На всех границах задавались граничные условия типа жесткого контакта. Сосредоточенный источник продольных волн типа «центр расширения» располагался в верхнем полупространстве на удалении 1000 м от кровли покрышки. В качестве временной функции в источнике принимались широкополосный аналитический однофазный импульс или его производная (двухфазный импульс) с максимумом модуля частотного спектра на 40 Гц.

Отметим, что помимо базовой модели рассматривались ее варианты, в которых значения мощности покрышки и коллектора были уменьшены: ^ = 50 м, ^ = 25 м и ^ = 25 м, ^ = 25 м. Кроме того, изучалась модель М-2 с пониженной контрастностью вмещающей среды.

2.2. Результаты численного моделирования волновых полей

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

На эпицентральном удалении И = 0.5 км, при котором наблюдается «лучевое» до-критическое прохождение падающей продольной волны (Р1) через высокоскоростную толщу, по точкам образования отраженных продольных (Р1Р1, Р12Р21, Р123Р321) и обменных волн (Р&8Л, Р&38321) на сейсмограмме НВСП надежно определяется кровля покрышки и менее надежно — кровля и подошва коллектора (рис. 2).

В случае эпицентрального удаления И = 1 км наиболее информативными оказываются поля, которые формируются обменной поперечной волной Р18234. Эта волна образуется на кровле высокоскоростной покрышки, «лучевым» образом распространяется через покрышку и коллектор и возбуждает отражения от их границ (рис. 3). Для наглядности на рис. 4 представлено разложение волнового поля НВСП на проходящие и восходящие поля продольных и поперечных волн. Отраженные волны Р18282, Р1823832 образуют на кровле покрышки помимо восходящих поперечных волн Р182821 и Р18238321 (рис. 4, г) также и обменные продольные волны Р18282Р1 и Р1823832Р1 (рис. 4, б).

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 I, С 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 !, с

Рис. 2. Волновое поле НВСП при Я = 0.5 км: а — горизонтальная компонента, б — вертикальная компонента. Надежно определяется кровля покрышки (Р1Р1) и менее надежно — кровля (Р12Р21, Р^!) и подошва коллектора (Р123Р321, Р1§23§321)

0.6 0.7 0.8 0.9 (, С 0.6 0.7 0.8 0.9 I, С

Рис. 3. Волновое поле НВСП при Я = 1 км: а — горизонтальная компонента, б — вертикальная компонента. Наиболее информативными оказываются поля Р182821 и Р18238321, Р18282Р1 и Р1823 832Р1, которые формируются обменной поперечной волной Р18234 (см. рис. 4). Р122 — слабое интерференционное поле продольной волны внутри покрышки

Рис. 4. Разложение волнового поля НВСП при Я = 1 км, вертикальные компоненты поля проходящих продольных волн (а), восходящих продольных волн (б), проходящих поперечных волн (в) и восходящих поперечных волн (г)

Рис. 5. Вертикальная компонента волнового поля НВСП: а — при Я = 1.5 км, б — при Я = 3 км

Поле продольной волны преломляется в высокоскоростную толщу (0 < г < 112) только в области докритических углов падения. На удалении Я = 1 км оно распространяется внутри высокоскоростной покрышки в виде слабого интерференционного поля Р122 (рис. 3, а), сформировавшегося за счет конструктивного наложения полей кратных волн. Преломленное в слой-коллектор поле Р1223 (рис. 3, а, б) вследствие малой интенсивности не дает значимых отражений от границ резервуара.

По мере увеличения эпицентрального расстояния Я > 1.5 км за счет экранирования происходит затухание (до полного исчезновения) поля Р1223 продольных волн (рис. 5, а), а при Я > 2 км начинается ослабление проходящего поля обменной поперечной волны Р1823 (рис. 5, б).

В случае «наземных» наблюдений ОТВ отраженные «лучевые» волны Р12Р21 и Р12ЗР321 наблюдаются только на малых горизонтальных удалениях (Я < 0.5 км), тогда как обменные волны Р18282Р1 и Р1823832Р1 выделяются на сейсмограмме ОТВ (рис. 6) в интервале удалений 1 - 2.5 км. На больших офсетах эти волны перестают существовать. При этом неверная интерпретация таких волн как чисто монотипных продольных может привести к ошибкам в определении границ резервуара. Поперечные обменные волны Р182821, Р18238321 также могут дать полезную информацию о коллекторе, особенно в случае регистрации горизонтальной компоненты волнового поля, где интенсивности этих волн будут максимальны.

В модели М-1 с уменьшенными значениями мощности покрышки и коллектора (^ = 50 м, ^ = 25 м) основные свойства полей проходящих и отраженных

Рис. 6. Модель М-1. Вертикальная компонента поля ОТВ на базе наблюдения 0 т 3 км. Отраженные волны Р12Р21 и РшР321 наблюдаются только на малых горизонтальных удалениях Я < 0.5 км, обменные волны Р^282Р1 и Р1823832Р1 выделяются в интервале удалений 1 т 2.5 км

волн в целом сохраняются такими, как в описанном выше случае при Ь2 = 100 м, Ь3 = 50 м. При этом уменьшаются временные задержки между отражениями от покрышки и коллектора и сокращается протяженность областей существования обменных продольных волн Р18282Р1 и Р1823832Р1 и поперечных волн Р182821 и Р18238321 на сейсмограммах ОТВ. При еще большем уменьшении мощности высокоскоростной покрышки поле продольной волны в области закритического падения проникает в коллектор в виде экранированной Р-волны, которая распространяется через высокоскоростной слой как неоднородное возмущение с коэффициентом затухания, зависящим от толщины слоя. При этом очевидно, что отклики от покрышки и коллектора будут находиться в тесной интерференции. Начиная с удаления Я = 2.5 км, формы импульсов отраженных волн Р1Р1 и Р181 стабилизируются и характеризуют суммарный отклик от всей высокоскоростной толщи «покрышка + коллектор».

Вычисление волновых полей выполнялось и для модели М-2 с пониженной контрастностью, в которой параметры покрышки и коллектора остаются такими, как в базовой модели, но увеличены скорости вмещающей среды (параметры покрывающего полупространства: ур1 = 3.0 км/с, у81 = 1.6 км/с, р1 = 12.2 г/см3; параметры подстилающего полупространства: ур4 = 3.2 км/с, у84 = 1.65 км/с, р4 = 2.3 г/см3).

0.6 0.7 0.8 0.9 1.0 1, с 1.0 1.1 1.2 1.3 1, С

Рис. 7. Модель М-2. Вертикальная компонента поля НВСП: а — при Я = 2 км, б — при Я = 3 км. На удалениях 2 и 3 км отклик от кровли и подошвы коллектора определяется отражен-

ными обменными волнами P1S2S2P1, P1S23S32P1 и P1S2S21, P1S23S321

В модели М-2 область докритического прохождения падающей продольной волны расширяется, вследствие чего при R = 1 км сохраняется «лучевое» прохождение Р-волной высокоскоростной толщи («покрышка + коллектор»), подобное тому, что наблюдается при R = 0.5 км в случае модели М-1. На еще больших удалениях R = 2 - 3 км (рис. 7) уже не наблюдается «лучевого» прохождения падающей Р-волны, и отклик от кровли и подошвы коллектора определяется отраженными обменными волнами PiS2S2Pi, PiS23S32P1 и P1S2S21, P1S23S321, как в случае модели М-1 при удалениях R = 1 - 2 км. На сейсмограммах ОТВ (рис. 8) для модели М-2 видно, что область «лучевой» регистрации отраженных продольных волн от коллектора ограничена эпицентральным удалением ~2 км. При этом монотипные продольные волны P12P21 и P123P321 прослеживаются на удалениях до 1.5 км, а обменные продольные волны P1S2S2P1 и P1S23S32P1 — на интервале удалений 2 - 5 км.

3. Модель тонкослоистого песчанистого резервуара

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

0 J-.-LJ-,-,-,-,-,-,-,-,-,-,-—r-

0.4 0.6 0.8 1.0 1.2 1.4 1 6 1.8 2.0 2.2 2.4 2.6 t, С

Рис. 8. Модель М-2. Вертикальная компонента поля ОТВ на базе наблюдений 0 т 5.5 км. Монотипные продольные волны Р12Р21 и Р123Рз21 прослеживаются на удалениях до 1.5 км, а обменные продольные волны Р18282Р1 и Р1823832Р1 — на интервале удалений 2 т 5 км

и времен вступлений обычных и конвертированных Р-волновых отражений на низких сейсмических частотах (10 Гц) оценивались в (Chabyshova and Goloshubin, 2014) при помощи асимптотического анализа коэффициентов отражения и обменов в случае нормального плосковолнового падения на тонкослоистый резервуар (Silin and Goloshubin, 2010). Расчеты в (Chabyshova and Goloshubin, 2014) показали, что амплитудные аномалии вследствие конверсии Р-волн могут иметь заметные временные задержки (до 200 мс) относительно отражения быстрой Р-волны, однако амплитуды модельных конвертированных Р-волн были гораздо ниже, чем амплитудные аномалии, наблюдаемые в разведочных случаях.

В данном исследовании изучается полное отраженное волновое поле от тонкослоистой пачки насыщенных упруго-пористых слоев Био, моделирующей ре-зервуарную зону. Программные коды OASP (Schmidt, 2004) допускают вычисление волновых полей в тонкослоистых средах, составленных из изотропных упруго-пористых слоев, с учетом всех возможных конверсионных переходов между быстрой и медленной продольными волнами, а также поперечной волной на всех границах модели среды. Описание упруго-пористой среды в комплексе OASES дается на основе известного формализма теории Био (Biot, 1962) с учетом механизма диссипации Био — Столл (Biot, 1962; Schmidt, 2004). На границах пористых слоев задаются стандартные граничные условия для случая проницаемых границ, описанные,

к примеру, в (Молотков, 2001; Schmidt, 2004). Кроме того, дополнительно использовались граничные условия контакта с проскальзыванием (Молотков, 2001). Условия такого контакта моделировались посредством включения на границе между слоями очень тонкого (~1 мм) флюидного пропластка.

3.1. Описание модели

Для описания модели среды были использованы данные из (Chabyshova and Goloshubin, 2014), отвечающие насыщенному резервуару, сложенному из тонкослоистого песчаника. Песчаник залегает между толщами глинистых сланцев. Полная мощность слоистого резервуара из 50 слоев составляет 20 м. Толщины пористых слоев варьируются от 0.1 до 0.6 м. Значения пористости в слоях изменялись по аналогии с (Chabyshova and Goloshubin, 2014) в пределах от 6 до 30 % Внутри резер-вуарной зоны скорости быстрой продольной волны изменялись от 3.2 до 4.8 км/с, а скорости S-волн — от 1.7 до 2.9 км/с. Верхняя часть резервуарной толщи насыщена газом, а нижняя — водой (рапой).

Для характеристики пористых слоев Био задавались значения следующих величин: коэффициента пористости ф, плотности флюида pf, модуля всестороннего сжатия флюида Kf, вязкости флюида п (в случае водного насыщения: pw = 1040 кг/м3, Kw = 2.8 ГПа, nw = 1.6 мПа-с; в случае газового насыщения: pgas = 0.34 кг/м3, Kgas = 0.29 ГПа, ngas = 0.01 мПа-с), плотности материала каркаса pg = 2650 кг/м3 и модуля всестороннего сжатия зерен каркаса Kr = 38 ГПа. Для каждого слоя вычислялись значения модуля сдвига ц и модуля всестороннего сжатия пористого каркаса K на основе заданного в (Chabyshova and Goloshubin, 2014) распределения скоростей внутри резервуарной зоны. Остальные параметры выражались через коэффициент пористости на основе эмпирических формул, приведенных в (Shmidt, 2004; Chabyshova and Goloshubin, 2014): для проницаемости пористой среды — к ~ 4ф2, для коэффициента извилистости порового пространства — а = 1 + 0.227(1-ф)/ф.

Фоновая среда выше и ниже резервуарной зоны принималась чисто упругой и характеризовалась значениями скоростей P и S волн и плотности (для верхнего полупространства: vp01 = 3.29 км/c, vs01 = 1.7 км/c, p01 = 2.4 г/см3, для нижнего полупространства: Vp02 = 3.6 км/c, vs02 = 1.8 км/c, Р02 = 2.4 г/см3), типичными для глинистых сланцев, окружающих такой резервуар.

Сосредоточенный источник возбуждения продольных волн (типа «центр расширения») расположен в верхнем полупространстве на расстоянии 1000 м от кровли резервуара. В качестве временной функции в источнике использовались двухфазные импульсы с максимумом амплитудного спектра на 10, 20, 30 и 40 Гц.

Для рассматриваемой модели среды рассчитывались сейсмограммы поля отраженных волн в точках (виртуального) горизонтального профиля наблюдений, проходящего внутри верхнего полупространства на уровне 500 или 900 м от кровли резервуара (так, чтобы можно было разделить по времени падающую и отраженные волны), в интервале горизонтальных удалений R = 25 ^ 500 м. Важно отметить, что с помощью программы OASP технически невозможно вычислять волновое поле в случае нормального падения волны, поскольку численное интегрирование ведется по волновому числу.

3.2. Результаты численного моделирования в случае тонкослоистого коллектора

Результаты вычислений по программе OASP полного волнового поля, отраженного от тонкослоистого насыщенного пористого резервуара, при разных частотах импульса падающей волны (10, 20, 30, 40 Гц) свидетельствуют об отсутствии какого-либо значимого по амплитуде конверсионного поля, связанного с обменами между быстрой и медленной продольными волнами на внутренних границах резервуара. На сейсмограммах ОТВ (где поле падающей прямой волны отсечено по времени) выделяются только поля отраженной продольной PP-волны и обменной поперечной PS-волны (рис. 9).

Отметим, что интенсивность вертикальной компоненты поля отраженной обменной PS-волны на малых офсетах (R ~ 100 м) составляет ~ 0.1 от интенсивности отраженной PP-волны (рис. 9). При этом временные запаздывания между PP- и PS-волнами могут иметь значения, сопоставимые с получаемыми в работе (Chabyshova and Goloshubin, 2014) для задержек между монотипной быстрой продольной волной и волнами, претерпевшими конверсионные переходы по схеме «быстрая-медленная-быстрая» внутри пачки тонких слоев.

0.5 0.6 0.7 0.8 0.9 1.0 t, с 0.5 ' о!б ' 07 ' 08 ' CL9 ' L0 .,с

Рис. 9. Вертикальная компонента волнового поля, отраженного от пачки пористых слоев, регистрация на уровне 900 м над кровлей резервуара: а — при 30 Гц, б — при 10 Гц

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

Рис. 10. Вертикальная компонента волнового поля, отраженного от пачки пористых слоев, в случае скользящего контакта на границе газового и водного насыщения; регистрация на уровне 900 м над кровлей резервуара: а — при 30 Гц, б — при 10 Гц

скользящем контакте на границе между газонасыщенной и водонасыщенной частями резервуара относительная амплитуда вертикальной компоненты поля отраженной Р8-волны увеличивается в несколько раз (рис. 10). При этом наибольшее усиление амплитуды РБ-волны отмечается на низких частотах (10 Гц), где амплитуды РР- и Р8-волн оказываются примерно одинаковыми на вертикальной компоненте поля, а на горизонтальной компоненте амплитуда Р8-волны становится значительно больше, чем у РР-волны. Подобное соотношение для амплитуд РР- и Р8-волн получается и в случае представления модели резервуара чисто упругими слоями, разделенными флюидными прослоями толщиной 1 мм. В частности, на рис. 11 приведено поле отраженных волн от пачки, состоящей из 27 непористых упругих слоев общей мощностью 15 м, на границах которых введены миллиметровые прослои, заполненные в верхней части резервуара газом, а в нижней части — водой.

Рис. 11. Вертикальная компонента волнового поля при 30 Гц, отраженного от пачки из 27 упругих (непористых) слоев общей мощностью 15 м с условиями проскальзывания на всех внутренних границах; регистрация на уровне 500 м над кровлей пачки

3.3. Возможные механизмы образования низкочастотных амплитудных аномалий

В работах (Ahmad et al., 2017; Chabyshova and Goloshubin, 2014) обсуждаются возможные причины образования низкочастотных амплитудных аномалий ниже залежей углеводорода и предлагаются два возможных физических механизма. Первый механизм — это вышеупомянутая гипотеза Чабышовой и Голошубина, предполагающая, что низкочастотные медленные Р-волны и конвертированные «быстрая-медленная-быстрая» продольные волны могут служить объяснением этого явления (Chabyshova and Goloshubin, 2014). В этой гипотезе низкочастотные аномалии связаны со степенью неоднородности в пористости и проницаемости резервуара. В исследовании (Ahmad et al., 2017) было обнаружено, что изучаемые два резервуара имеют характер двойной пористости, в которых интервалы с относительно высокой эффективной пористостью (высокая проницаемость) способствуют течению, в то время как материал с низкой пористостью сохраняет флюид (Barenblatt et al., 1960). В работе (Ahmad et al., 2017) отмечается также, что распространение волн через резервуар такого типа (Био — Баренблатт) может привести к генерации медленных Р-волн на каждой проницаемой/непроницаемой границе. Общая сумма модовых конверсий может привести к большему затуханию высокочастотного содержимого и задержке низкочастотных отражений более чем на 100 мс, как описано Чабышовой и Голошубиным (Chabyshova and Goloshubin, 2014). Более того, конверсионные переходы типа «быстрая-медленная-быстрая» Р-волны в тонкослоистых резервуарах и их интерференция могут способствовать высоким амплитудам резонансного рассеяния и резонансного усиления. Авторы цитированных работ полагают, что эти возможные механизмы требуют дальнейшего изучения, особенно интерференционных эффектов, возникающих за счет конверсионных переходов на границах коллектора.

В качестве второго механизма авторы в (Ahmad et al., 2017) рассматривают явление возбуждения в упругих средах, содержащих трещины, заполненные жидким флюидом, слабо затухающих низкочастотных (приблизительно 10 Гц), низкоскоростных (приблизительно 0.3 км/с) и высокоамплитудных волн Крауклиса (Krauklis et al., 1996). Однако в реальных примерах сейсмических данных, показанных в (Ahmad et al., 2017), не наблюдались подобные низкочастотные высокоамплитудные волновые поля, которые имели бы сходные с описанными в (Krauklis et al., 1996) временные задержки и интервальные скорости между кровлей резервуара и ниже его подошвы.

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

4. Заключение

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

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

Подводя итоги результатов численных исследований полноволнового сейсмического отклика от тонкослоистой пачки насыщенных упруго-пористых слоев Био, моделирующей резервуарную зону, можно сделать вывод, что обсуждаемые в работах (Chabyshova and Goloshubin, 2014; Ahmad et al., 2017) низкочастотные амплитудные аномалии, наблюдаемые на сейсмических разрезах ниже резервуарной зоны, не удается объяснить влиянием каких-либо типов модовых обменов между быстрой и медленной продольными волнами, возникающих в модели тонкослоистого пористого насыщенного резервуара, оставаясь в рамках обычной теории по-роупругости Био (Biot, 1962). Отсутствие конверсионного продольного волнового поля в наших вычислениях, по-видимому, вызвано тем, что, согласно исследованиям известных частотно-зависимых механизмов диссипации в рамках модели Био (Muller et al., 2010), в области сейсмических частот в пористой среде Био распространяются только быстрая продольная и поперечная волны, а вторая продольная мода является диффузионной (она становится распространяющейся волной только на достаточно высоких частотах >1 КГц) и успевает затухнуть даже в очень тонких слоях. Возможно, использование иных, более сложных теоретических моделей пористых сред, рассматриваемых в (Muller et al., 2010), позволит объяснить данный механизм.

Результаты наших численных модельных исследований показывают, что амплитудные аномалии, подобные описанным в (Chabyshova and Goloshubin, 2014; Ahmad et al., 2017), могут быть вызваны образованием поля отраженных обменных поперечных волн, запаздывающих относительно отклика продольной волны и имеющих на малых и средних офсетах значимые амплитуды на вертикальной компоненте поля. Амплитуды и временные задержки суммарной обменной поперечной волны могут варьироваться в широких пределах в зависимости от конкретных условий: глубины залегания и мощности резервуара, частоты падающей продольной волны, наличия в резервуарной зоне границ с контактом проскальзыва-

ния или протяженных тонких трещин, заполненных флюидом, и других факторов. Отметим также, что наличие в резервуарной зоне нескольких границ с контактом проскальзывания может приводить к формированию внутри таких интервалов среды интенсивных интерференционных полей кратных поперечных волн, переизлучающих на границах часть своей энергии в виде слабого поля обменных продольных волн типа SP (Голикова и др., 2006). Такое вторичное поле SP-волн будет иметь некоторое запаздывание относительно обычных отраженных Р-волн, а его частота будет зависеть от толщины указанных интервалов. Отметим, что в случае модели тонкослоистого резервуара с контактами проскальзывания на границах будет формироваться высокочастотное поле вторичных SP-волн, поэтому такая тонкослоистая модель не дает объяснения образования низкочастотных амплитудных аномалий.

Литература

Голикова, Г. В., Дараган-Сущова, Л. А., Ковтун, Ал. А., Лыскова, Е. Л., Санников, К. Ю. (2016а). Особенности сейсмического волнового поля в северной части Баренцева моря по результатам численного моделирования. Региональная геология и металлогения, 67, 70-78. Голикова, Г. В., Ковтун, А. А., Лыскова, Е. Л., Санников, К. Ю. (2016б). Влияние слоев высокой скорости на формирование сейсмических волновых полей. Технологии сейсморазведки, 4, 33-44. https://doi.org/10.18303/1813-4254-2016-4-33-44 Голикова, Г. В., Ковтун А. А., Решетников, В. В. (2006). Численные исследования интерференционных волновых полей в слоистых средах, содержащих границы разделов с контактом проскальзывания. Вопросы геофизики. Выпуск 39. (Ученые записки СПбГУ; № 439), 20-37. Голошубин, Г. М., Силин, Д. Б., Баюк, И. О. (2015). Продольные сейсмические волны в проницаемых средах: асимптотическое представление. Технологии сейсморазведки, 4, 15-22. https://doi. org/10.18303/1813-4254-2015-4-15-22 Левянт, В. Б., Хромова, И. Ю., Козлов, Е. А., Керусов, И. Н., Кащеев, Д. Е., Колесов, В. В., Мармалев-ский, Н. Я. (2010). Методические рекомендации по использованию данных сейсморазведки для подсчета запасов углеводородов в условиях карбонатных пород с пористостью трещинно-ка-вернозного типа. Москва: ЦГЭ. Молотков, Л. А. (2001). Исследование распространения волн в пористых и трещиноватых средах на

основе эффективных моделей Био и слоистых сред. Санкт-Петербург: Наука. Ahmad, S. S., Brown, R. J., Escalona, A., Rosland, B. O. (2017). Frequency-dependent velocity analysis and offset-dependent low-frequency amplitude anomalies from hydrocarbon-bearing reservoirs in the southern North Sea, Norwegian sector. Geophysics, 82(6), N51-N60. https://doi.org/10.1190/ GEO2016-0555.1

Barenblatt, G. I., Zheltov, I. P., Kochina, I. N. (1960). Basic concepts in the theory of seepage of homogeneous

liquids in fissured rocks [strata]. Journal of Applied Mathematics and Mechanics, 24, 1286-1303. Biot, M. A. (1962). Generalized theory of acoustic propagation in porous dissipative media. J. Acoust. Soc. Amer., 34(5), 1254-1264.

Bruvoll, V., Kristoffersen, Y., Coakley, B. J., Hopper, J. R., Planke, S., Kandilarov, A. (2012). The nature of acoustic basement on Mendeleev and northwestern Alpha ridges, Arctic Ocean. Tectonophysics, 514517, 123-145. https://doi.org/10.1016/j.tecto.2011.10.015 Chabyshova, E., Goloshubin, G. (2014). Seismic modeling of low-frequency "shadows" beneath gas reservoirs. Geophysics, 79(6), D417-D423. https://doi.org/10.1190/geo2013-0379.! Krauklis, P., Goloshubin, G., Krauklis L. (1996). A slow wave in a fluid layer simulating an oil-saturated

seam. Journal of Mathematical Sciences, 79(4), 1224-1230. https://doi.org/10.1007/BF02362888 Muller, T., Gurevich, B., Lebedev, M. (2010). Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks — review. Geophysics, 75(5), 75A147-75A164. https://doi. org/10.1190/1.3463417

Schmidt, H. (2004). OASES Version 3.1. User Guide and Reference Manual. Department of Ocean Engineering Massachusetts Institute of Technology.

Schmidt, H., Jensen, F. (1985). A full wave solution for propagation in multilayered viscoelastic media with application to Gaussian beam reflection at fluid-solid interfaces. J. Acoust. Soc. Amer., 77, 813-825. https://doi.org/10.1121/1.392050 Silin, D., Goloshubin, G. (2010). An asymptotic model of seismic reflection from a permeable layer. Transport in Porous Media, 83, 233-256. https://doi.org/10.1007/s11242-010-9533-8

Статья поступила в редакцию 30 января 2019 г.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Статья рекомендована в печать 27 ноября 2019 г.

Контактная информация:

Ковтун Александр Андреевич — al_an_kovtun_46@mail.ru Лыскова Евгения Леонидовна — e.lyskova@spbu.ru Санников Константин Юрьевич — k.sannikov@spbu.ru

The study of seismic wave fields in reservoir zones by numerical simulation

of the multi-offset vertical profiling data and data of common source point technique

A. A. Kovtun, E. L. Lyskova, K. Yu. Sannikov St. Petersburg State University,

7-9, Universitetskaya nab., St. Petersburg, 199034, Russian Federation

For citation: Kovtun, A. A., Lyskova, E. L., Sannikov, K. Yu. (2020). The study of seismic wave fields in reservoir zones by numerical simulation of the multi-offset vertical profiling data and data of common source point technique. Vestnik of Saint Petersburg University. Earth Sciences, 65 (1), 33-50. https://doi.org/10.21638/spbu07.2020.103 (In Russian)

This paper presents a study of seismic wave fields in reservoir zones by numerical simulation. In the first part of the study, numerical modeling of full wavefield was performed to synthesize the multi-offset vertical profiling data and data of OPV (Common Source Point Technique) for the model of the medium in which high-velocity layers (seal and collector) are placed in a low-velocity containing medium. It has been established that the principal boundary, the roof of the collector, can be distinguished by data only in the case of great thickness of high-velocity seal from the field of converted P1S2S2P1 and P1S2S21 waves. These waves are formed on the roof of the seal during the passage of a converted transverse wave and exist only in a certain interval of epicentral distances from the source, the length of which depends on the contrast of the high-velocity seal, as well as on its thickness and depth. Reliable determination of the boundary between the high-velocity seal and the collector can be difficult when using only ground-based seismic data due to the screening effect. However, it becomes possible with the use of multi-offset vertical profiling materials and the correct use of the data of reflected converted P1S2S2P1 and P1S2S21 waves. In the second part of the work, we analyze the results of numerical model studies of the full-wave seismic response from a thin-layered pack of saturated elastic porous Biot layers simulating a reservoir zone. The results of numerical model studies show that the low-frequency amplitude anomalies observed in seismic sections below the reservoir zone may be caused by the formation of a field of reflected converted transverse waves that are delayed relative to the response of the longitudinal wave and having significant amplitudes on the vertical field component for the small and middle offsets. The amplitudes and time delays of the total converted transverse wave can vary widely depending on specific conditions: depth and reservoir thickness, incident longitudinal wave frequency, the presence of boundaries in the reservoir area with sliding contact or extended thin cracks filled with fluid, and other factors. Keywords: numerical modeling, seismic fields, collector, screening effect, Biot layers.

References

Ahmad, S. S., Brown, R. J., Escalona, A., Rosland, B. O. (2017). Frequency-dependent velocity analysis and offset-dependent low-frequency amplitude anomalies from hydrocarbon-bearing reservoirs in the southern North Sea, Norwegian sector. Geophysics, 82(6), N51-N60. https://doi.org/10.1190/ GEO2016-0555.1

Barenblatt, G. I., Zheltov, I. P., Kochina, I. N. (1960). Basic concepts in the theory of seepage of homogeneous

liquids in fissured rocks [strata]. Journal of Applied Mathematics and Mechanics, 24, 1286-1303. Biot, M. A. (1962). Generalized theory of acoustic propagation in porous dissipative media. J. Acoust. Soc. Amer., 34(5), 1254-1264.

Bruvoll, V., Kristoffersen, Y., Coakley, B. J., Hopper, J. R., Planke, S., Kandilarov, A. (2012). The nature of acoustic basement on Mendeleev and northwestern Alpha ridges, Arctic Ocean. Tectonophysics, 514517, 123-145. https://doi.org/10.1016Zj.tecto.2011.10.015 Chabyshova, E., Goloshubin, G. (2014). Seismic modeling of low-frequency "shadows" beneath gas reservoirs. Geophysics, 79(6), D417-D423. https://doi.org/10.1190/geo2013-0379.! Golikova, G. V., Daragan-Sushchova, L. A., Kovtun A. A., Lyskova, E. L., Sannikov, K. Yu. (2016a). Features of seismic wave field in the northern part of the Barents Sea from numerical simulation. Regional'naia geologiia i metallogeniia, 67, 70-78. (In Russian) Golikova, G. V., Kovtun A. A., Lyskova, E. L., Sannikov, K. Yu. (2016b). The effect of high-velocity layers on seismic wave field formation. Tekhnologii seismorazvedki, 4, 33-44. (In Russian) https://doi. org/10.18303/1813-4254-2016-4-33-44 Golikova, G. V., Kovtun A. A., Reshetnikov, V. V. (2006). Numerical investigation of interference wave field are affected by sliding contacts in stratified media. Voprosygeofiziki. Vypusk 39. (Uchenye zapiski SPb-GU; N 439), 20-37. (In Russian) Goloshubin, G. M., Silin, D. B., Bayuk, I. O. (2015). Compressional seismic waves in permeable media: an asymptotic representation. Tekhnologii seismorazvedki, 4, 15-22. (In Russian) https://doi. org/10.18303/1813-4254-2015-4-15-22. Krauklis, P., Goloshubin, G., Krauklis L. (1996). A slow wave in a fluid layer simulating an oil-saturated

seam. Journal of Mathematical Sciences, 79(4), 1224-1230. https://doi.org/10.1007/BF02362888 Leviant, V. B., Khromova, I. Iu., Kozlov, E. A., Kerusov, I. N., Kashcheev, D. E., Kolesov, V. V., Marma-levskii, N. Ia. (2010). Methodical recommendations on the use of seismic data for the calculation of hydrocarbon reserves in the conditions of carbonate rocks with a fractured-cavern-type porosity. Moscow: TsGe Publ. (In Russian)

Molotkov, L. A. (2001). Investigation of wave propagation in porous and fractured media based on effective

Biot and layered media models. St. Petersburg: Nauka Publ. (In Russian) Muller, T., Gurevich, B., Lebedev, M. (2010). Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks — review. Geophysics, 75(5), 75A147-75A164. https://doi. org/10.1190/1.3463417

Schmidt, H. (2004). OASES Version 3.1. User Guide and Reference Manual. Department of Ocean Engineering Massachusetts Institute of Technology. Schmidt, H., Jensen, F. (1985). A full wave solution for propagation in multilayered viscoelastic media with application to Gaussian beam reflection at fluid-solid interfaces. J. Acoust. Soc. Amer., 77, 813-825. https://doi.org/10.1121/1.392050 Silin, D., Goloshubin, G. (2010). An asymptotic model of seismic reflection from a permeable layer. Transport in Porous Media, 83, 233-256. https://doi.org/10.1007/s11242-010-9533-8

Received: January 30, 2019 Accepted: November 27, 2019

Contact information:

Alexander A. Kovtun — al_an_kovtun_46@mail.ru Eugenia L. Lyskova — e.lyskova@spbu.ru Konstantin Yu. Sannikov — k.sannikov@spbu.ru

i Надоели баннеры? Вы всегда можете отключить рекламу.