Сейсмогеологическое моделирование с целью определения влияния полноты исходной информации и геологических условий на результат прогноза емкостных свойств коллекторов по сейсмическим данным

Егоров С.В., ООО «Газпромнефть НТЦ» Egorov.SVa@gazpromneft-ntc.ru

Приезжев И.И., РГУ нефти и газа (НИУ) имени И.М. Губкина ivanpriez@gmail.com

Журнал «Нефтегазовая геология. Теория и практика»

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

Ключевые слова: геологическое моделирование, сейсмическое моделирование, алгоритм машинного обучения, фильтрационно-емкостные свойства коллектора.

Введение

Последние несколько лет все больше научных работ в нефтегазовой отрасли посвящены вопросам внедрения современных цифровых технологий с целью увеличения эффективности работы с большими объемами информации (Big Data). Особую роль среди этих технологий занимают алгоритмы на основе машинного обучения, которые уже на протяжении последнего десятилетия эффективно используются в различных отраслях науки и техники [Флах, 2015]. Данные алгоритмы позволяют с высокой точностью решать задачи кластеризации, распознавания и прогнозирования.

Одним из наиболее обширных и не до конца изученных источников информации в нефтегазовой отрасли является сейсмическое волновое поле. Данные сейсморазведочных работ - основной источник пространственной информации на геологоразведочной стадии изучения недр. Преобразования позволяют получить информацию в двухмерном или трехмерном виде о физических свойствах изучаемого свойства в межскважинном пространстве. Последнее время для извлечения этой информации стали активно применять алгоритмы на основе машинного обучения.

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

Данная статья посвящена изучению возможности использования алгоритмов машинного обучения взамен линейной функции. Важно заметить, что изучаемый вопрос освещается в работах многих отечественных [Ампилов, 2008; Ампилов, Барков, Яковлев, 2008; Приезжев, Егоров, 2017, 2018; Приезжев, Егоров, Гладков, 2018] и зарубежных ученых [Chopra, Marfurt, 2007; Herron, Latimer, 2011]. Методика использования нейронных сетей описана в инструкции Центральной геофизической экспедиции [Методические рекомендации..., 2006]. Ученые отмечают высокую точность алгоритмов при работе с обучающей выборкой (входные данные) и низкую при проверке результатов на контрольной выборке. Проблема в терминологии анализа данных именуются «переобучением» и до сих пор не решена полностью [Флах, 2015]. Однако, одной из причин недавнего подъема интереса изучаемых алгоритмов является разработка процедур, позволяющих существенно сократить влияние описанной проблемы. Вопрос использования современных решений машинного обучения для прогноза фильтрационно-емкостных свойств (ФЕС) коллекторов с помощью сейсмических данных рассматривался преимущество зарубежными исследователями, в их работах рассматриваемые алгоритмы используются для решения задачи инверсии [Schuster, Chen, 2019].

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

Методология

Целью выполненного моделирования являлась оценка влияния следующих параметров на результат прогнозных построений:

1) обстановка осадконакопления;

2) характер связи между прогнозируемым свойством и волновым полем;

3) количество скважин и освещенность ими геологического разреза;

4) тип предсказываемого емкостного параметра коллектора (эффективная толщина и средняя пористость);

5) сейсмические атрибуты, не имеющиеся связи с прогнозируемым свойством;

6) преобразующая функция.

Синтетическое моделирование выполнялось в несколько этапов (рис. 1).

1.1.PNG

Первый этап. Построены трехмерные геологические модели по 2 пластам (одного из месторождений Западной Сибири) с условными названиями RES1 и RES2.

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

Третий этап. Расчёт сейсмических атрибутов по каждому синтетическому кубу.

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

Четвертый этап. Моделирование различных вариантов расстановки скважин для определения влияния количества скважин и их размещения на достоверность результата прогнозирования.

Пятый этап. Расчет карт средней пористости и эффективных толщин по пластам RES1 и RES2. Карты использовались для двух целей. Во-первых, для снятия значений в точках скважин, размещенных в предыдущем этапе. Во-вторых, для контроля качества прогнозных построений.

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

Обсуждение

Поэтапное описание технических деталей выполненных процедур

Первый этап. Созданные геологические модели различаются пространственным распределением свойств. Для пласта RES1 характерно относительно плавное изменение свойств по латерали и более широкий диапазон значений эффективной толщины (0-24 м). В тоже время, характерной чертой пласта RES2 являются высокая латеральная изменчивость и более узкий диапазон значений эффективной толщины (0-16 м). Кроме того, зависимость между картами эффективной толщины и средней пористости в пластах RES1 и RES2 существенно различается. Причиной этому служат различия в вертикальном строении пластов: отложения RES1 сложены мощными и выдержанными толщами «коллектора» и «не коллектора», а пласт RES2 представлен тонким переслаиванием пород (рис. 2). Две модели созданы с целью количественной оценки влияния различных исходных геологических условий (обстановок осадконакопления) на точность прогноза.

Второй этап. Чтобы оценить, как на точность прогноза влияет характер связи между волновым полем и прогнозируемыми свойствами, сейсмические данные синтезировались 3 путями (рис. 3).

1 способ

1) Куб плотности получен путем интерполяции данных в скважинах, с использованием алгоритма кокрикинг, где в качестве тренда выступал куб пористости (тесная корреляционная связь подтверждена результатами сопоставления показаний кривой ГГКП и керна).

2) Куб скоростей продольных волн получен путем трансформации куба плотности через линейную функцию.

2.PNG

3.PNG

3) Затем путем перемножения двух созданных объемных распределений скорости и плотности получен куб акустического импеданса.

4) Трансформация во временную область проведена с использованием куба интервальных скоростей, рассчитанного с применением карт интервальных скоростей, полученных по каждому пласту.

5) Затем выполнена процедура свертки, где в качестве сигнала использовался теоретический вейвлет Рикера с центральной частотой в 30 Гц.

2 способ

1) На втором этапе куб скорости получен через кокрикинг, где в качестве тренда выступал куб плотности. Заданный коэффициент корреляции между основной переменой и вспомогательной составил 0,8 (оценен по данным ГИС реального месторождения). Все остальные этапы аналогичны первому способу.

3 способ

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

3) В рамках каждого способа генерации сейсмических данных, для приближения синтетического материала к реальным сейсмическим данным, к результирующим кубам добавлялся случайный шум (15%), а скоростная модель рассчитывалась таким образом, чтобы добавить определенные искажения (в пределах 10 мс) во временные поверхности по отношению к глубинным.

Третий этап. По каждому синтетическому кубу рассчитан следующий набор сейсмических атрибутов: карта среднеквадратичных амплитуд, карты максимальных и минимальных значений амплитуд, стратиграфическая нарезка слайсов куба амплитуд от кровли до подошвы изучаемого интервала. Кроме того, для оценки влияния «физически не обоснованных» атрибутов на результат прогноза, во входные данные добавлен параметр двойного времени пробега продольной волны до горизонта. Важно заметить, что ни один из используемых атрибутов не имеет высокой корреляционной связи (коэффициент корреляции не превышает значения 0,4) с картами пористости и эффективной толщины, что не позволяет выполнить прогноз, привлекая лишь один параметр.

Четвертый этап. Смоделировано несколько возможных случаев разведки месторождения (рис. 4):

1. пробурено 10 скважин, бурением не охвачена зона максимальных ФЕС;

2. пробурено 10 скважин, бурением охвачены все зоны;

3. пробурено 30 скважин, бурением не охвачена зона максимальных ФЕС;

4. пробурено 30 скважин, бурением охвачены все зоны.

Пятый этап. Чтобы оценить влияние прогнозируемого свойства на результат, для прогноза выбраны два параметра ФЕС: эффективная толщина и средняя открытая пористость. Выбор обусловлен тем, что данные параметры используются для подсчета запасов и оценки ресурсов объемным методом и, как правило, могут оказывать ощутимое влияние на волновое поле. Стоит так же отметить, что использованные методики синтеза сейсмических данных обеспечивают изменение динамики волновой картины в зависимости от значений пористости и эффективной толщины, что, в свою очередь, является необходимым фактором для решения поставленных задач.

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

4.PNG

Многомерная линейная функция представляет собой вариацию линейной функции для задачи, где предсказывание происходит не по одной переменной, а по нескольким, каждой из которых присваивается свой весовой коэффициент. Проблемой данного алгоритма, как и многих методов, где используется множество переменных, является мультиколлинеарность или наличие корреляционной связи между предсказывающими переменными, которая влечет за собой некорректную оценку весовых коэффициентов [Флах, 2015]. Для решения данной проблемы целевая функция уравнения рассчитывалась с учетом коэффициента Тихонова [Tikhonov, Arsenin, 1977].

Нейронные сети – огромная группа алгоритмов, которых объединяет аппаратная реализация, построенная по образу и подобию сетей нервных клеток живого организма. В данной работе использовался алгоритм на основе гибридного метода обучения [Priezzhev, Kobrunov, 2016; Приезжев, Егоров, 2017].

ACE - метод итерационного вычисления условного математического ожидания. Данный алгоритм подбирает отдельную преобразующую функцию для каждой предсказывающей переменной [Breiman, Friedman, 1985].

«Случайный лес» – алгоритм, основанный на теории решающих деревьев и бутстрэпа. Решающее дерево и деревья решений – это способ представления правил в иерархической, последовательной структуре, где каждому объекту соответствует единственный узел, дающий решение. Другими словами, это метод находит решение, обучаясь на исходных данных путем их последовательной классификации по имеющимся признакам [Приезжев, Егоров, 2018].

Метод ближайшего соседа - согласно этому методу в пространстве признаков находятся несколько объектов обучения (в нашем случае это скважины) с максимально похожими свойствами (набору атрибутов) на текущую точку изучаемого пласта. Следующим шагом прогнозируемое свойство интерполирует в текущую точку пласта на основе параметров близости значений используемых атрибутов в этой точке и значений атрибутов найденных скважин. При этом используется параметр количества ближайших соседей (скважин). Если он равен единице, то значение прогнозируемого параметра сносится в текущую точку. Если используется несколько ближайших соседей (скважин), то можно найти весовое среднее параметра в зависимости от степени его близости к скважине [Приезжев, Егоров, Гладков, 2018].

Результаты

Характер связи между прогнозируемым свойством и упругим параметром. На рис. 5 показано сравнение прогнозных карт пористости по трем кубам, полученных при помощи алгоритма на основе нейронных сетей по данным пласта RES1, где 30 скважин расположены равномерно. Видно, что различия между прогнозными картами, рассчитанными из куба, полученного способом 1 и 2, несущественны, что является следствием фактического отсутствия различий между типами зависимости пористости от акустического импеданса - в обоих случаях одна линейная зависимость, не зависимо от литологии. В свою очередь, прогноз, полученный к кубу 3, характеризуется худшей корреляцией с картой пористости из геологической модели. Таким образом, результаты демонстрируют, что с усложнением литологического строения пласта (увеличение количества литотипов, обособленных в поле упругих параметров) ухудшается качество прогноза.

Тип предсказываемого параметра коллектора. Согласно результатам проведённых тестовых работ, точность прогноза пористости ощутимо выше предсказывания эффективной толщины (максимальный коэффициент корреляции с соответствующими картами из геологической модели порядка 0,8 для пористости и 0,65 для толщины). Одной из основных причин данного результата является наличие корреляционной связи между значениями объёмной плотности и пористостью пород (обоснованной по скважинным исследованиям), которая заложена в синтетические модели. Кроме того, важно отметить тот факт, что само понятие коллектор и методы его определения являются довольно относительными [Бакирова, 1990].

Обстановка осадконакопления. Прогнозные карты, полученные по пласту RES1, характеризуются более высокой достоверностью нежели по пласту RES2. Объяснить результат можно ключевым различием в строении моделируемых пластов: мощностью и распределением коллекторов в толще пород. В отложениях RES1 толщина пропластков коллектора, как правило, составляет 3-10 м, а сосредоточены пропластки вблизи друг друга. В толще RES2 мощности коллекторов изменяются в пределах первых метров, которые распределены равномерно по всему пласту. Вследствие вышеописанного, часть песчаных прослоев нивелируется в процессе преобразования геологической модели в синтетический сейсмический куб из-за интерференции.

Количество скважин и освещенность ими геологического разреза. Карты, полученные с использованием 10 и 30 скважин, не охватывающих зону максимальной пористости, демонстрируют, что точность прогноза существенно зависит от расположения скважин. В тоже время, однозначно лучшая сходимость прогноза по 30 скважинам, чем по 10, охватывающим все зоны, показывает, что количество данных так же не мало важно. Таким образом, прогноз в условиях неполной охарактеризованности бурением изучаемого разреза маловероятен даже на качественном уровне.

5.PNG

Сейсмические атрибуты, не имеющие связи с прогнозируемым свойством. Карты, полученные с учетом и без «некорректного» атрибута, практически не различаются по качественным и количественным критериям оценки. Причиной слабого влияния на результат, несомненно, является наличие большого количества атрибутов (14 шт.) и низкой корреляции данного параметра со скважинной информацией (0,1–0,2).

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

Выводы

Исследованиями показано, что при условии расстановки скважин в пределах одной геолого-геофизической зоны и сложного литологического строения коллектора, прогноз на основе поиска зависимости между картой сейсмического атрибута и РИГИС (скважины) практически бесполезен. В тоже время, если бурением охарактеризованы зоны пониженных и повышенных ФЕС пласта, а петроупругое моделирование показывает хотя бы одну четкую зависимость между свойствами пласта и упругими методами ГИС, то регрессионные алгоритмы на основе карт атрибутов могут дать достоверный прогноз. Среди всех протестированных алгоритмов наиболее точным оказались нейронные сети (с регуляризацией), что свидетельствует о том, что методы машинного обучения позволяют получить более достоверный прогноз в сравнении с линейной регрессией в идентичных условиях.

Литература

Ампилов Ю.П. От сейсмической интерпретации к моделированию и оценке месторождений нефти и газа. - Москва: Изд-во Спектр, 2008. - 384 с.

Ампилов Ю.П., Барков А.Ю., Яковлев А.В. Применение нейронных сетей для прогноза фильтрационно-емкостных свойств по сейсмическим атрибутам и результатам акустической инверсии: стендовый доклад Р105 // Геонауки – от новых идей к новым открытиям: третья Международная выставка-конференция (г. Санкт-Петербург, 7-10 апреля 2008 г.). - СПб.: SEG-EAGE, 2008. – С. 105.

Бакирова Э.А. Геология нефти и газа. - М.: Недра - 1990. - 240 с.

Методические рекомендации по использованию данных сейсморазведки (2D, 3D) для подсчета запасов нефти и газа /   Ю.П.   Ампилов,   В.М.   Глоговский,   В.В. Колесов, М.Б. Коростышевский, В.Б. Левянт, С.Н. Птецов. - М.: ОАО «ЦГЭ», 2006. - 39 с.

Приезжев И.И., Егоров С.В. Гибридное обучение нейронных сетей с целью прогноза параметров нефтегазовой продуктивности горных пород // Сейсмические технологии. - 2017. - С. 205 - 208.

Приезжев И.И., Егоров С.В. Прогноз кубов упругих свойств по данным сейсморазведки 3D и ГИС при помощи алгоритма «случайного леса» // Геофизика. - 2018. - №2. - С. 10-16.

Приезжев И.И., Егоров С.В., Гладков Е.А. Применение алгоритмов машинного обучения для решения задач количественного прогноза ФЕС по сейсмическим и скважинным данным // Геофизика. - 2018. - №3. - С.33-38.

Флах П. Машинное обучение. Наука и искусство построения алгоритмов, которые извлекают знания из данных. - М.: ДМК Пресс, 2015. - 400 с.

Breiman L., Friedman J.H. Estimating Optimal Transformations for Multiple Regression and Correlation // Journal of American Statistical Association. - 1985. - Vol. 80. - No. 391. - P. 580-598.

Chopra S., Marfurt K.J. Seismic attributes for prospect identification and reservoir characterization. - Tulsa: Society of Exploration Geophysicists, 2007. - 481 p.

Herron D.A., Latimer R.B. First Steps in Seismic Interpretation // Geophysical monograph series, no. 16. - Tulsa: Society of Exploration Geophysicists, 2011. - 203 p.

Priezzhev I., Kobrunov А. Hybrid combination genetic algorithm and controlled gradient method   to   train   a   neural   network   //   Geophysics.   –   2016.   -   Vol   81.   –   P.   35-43. DOI: https://doi.org/10.1190/geo2015-0297.1

Schuster G.T., Chen Yu. Seismic Inversion by Newtonian Machine Learning // Research Gate GEO-2019-0434.R2: Preprint (PDF Available) April 2019. - 57 р. - https://www.researchgate.net/publication/332630844_Seismic_Inversion_by_Newtonian_Machine_ Learning

Tikhonov A.N., Arsenin V.Y. Solutions of ill-posed problems. - Washington: V.H. Winston and Sons. – 1977. – P. 259.



Egorov S.V.

LLC Gazpromneft STC, St. Petersburg, Russia, Egorov.SVa@gazpromneft-ntc.ru

Priezzhev I.I.

National University of Oil and Gas "Gubkin University" (Gubkin University), Moscow, Russia, ivanpriez@gmail.com

SEISMOGEOLOGICAL MODELING IN ORDER TO DETERMINE THE INFLUENCE OF THE COMPLETENESS OF INITIAL INFORMATION

AND GEOLOGICAL CONDITIONS ON THE RESULT OF FORECAST OF PORO-PERM RESERVOIRS PROPERTIES

The article is devoted to the study of the predictive capabilities of several of the most common machine learning algorithms in various geological and geophysical conditions. To solve this problem, several synthetic models have been calculated that simulate the possible conditions for reservoir properties. It is shown which algorithms allow to achieve the better positive effect.

Keywords: geological modeling, seismic modeling, machine learning algorithm, poro-perm reservoir properties.

References

Ampilov Yu.P. Ot seysmicheskoy interpretatsii k modelirovaniyu i otsenke mestorozhdeniy nefti i gaza [From seismic interpretation to modeling and evaluation of oil and gas fields]. Moscow: Izd-vo Spektr, 2008, 384 p.

Ampilov Yu.P., Barkov A.Yu., Yakovlev A.V. Primenenie neyronnykh setey dlya prognoza fil'tratsionno-emkostnykh svoystv po seysmicheskim atributam i rezul'tatam akusticheskoy inversii: stendovyy doklad R105 [Тhe implementation of neural network to рredict reservoir characteristics from seismic attributes and inversion results]. Geonauki – ot novykh idey k novym otkrytiyam: tret'ya Mezhdunarodnaya vystavka-konferentsiya (St. Petersburg, 7-10 Apr 2008). St. Petersburg: SEG- EAGE, 2008, p. 105.

Bakirova E.A. Geologiya nefti i gaza [Petroleum geology]. Moscow: Nedra, 1990, 240 p. Breiman L., Friedman J.H. Estimating Optimal Transformations for Multiple Regression and

Correlation // Journal of American Statistical Association, 1985, vol. 80, no. 391, pp. 580-598.

Chopra S., Marfurt K.J. Seismic attributes for prospect identification and reservoir characterization. - Tulsa: Society of Exploration Geophysicists, 2007, 481 p.

Flakh P. Mashinnoe obuchenie. Nauka i iskusstvo postroeniya algoritmov, kotorye izvlekayut znaniya iz dannykh [Machine learning. The science and art of building algorithms that extract knowledge from data]. Moscow: DMK Press, 2015, 400 p.

Herron D.A., Latimer R.B. First Steps in Seismic Interpretation // Geophysical monograph series, no. 16. - Tulsa: Society of Exploration Geophysicists, 2011, 203 p.

Metodicheskie rekomendatsii po ispol'zovaniyu dannykh seysmorazvedki (2D, 3D) dlya podscheta zapasov nefti i gaza [Guidelines for using seismic data (2D, 3D) to calculate oil and gas reserves]. Yu.P. Ampilov, V.M. Glogovskiy, V.V. Kolesov, M.B. Korostyshevskiy, V.B. Levyant,

S.N. Ptetsov. Moscow: OAO «TsGE», 2006, 39 p.

Priezzhev I., Kobrunov А. Hybrid combination genetic algorithm and controlled gradient method    to    train    a    neural    network    //    Geophysics,    2016,    vol.    81,    pp.    35-43. DOI: https://doi.org/10.1190/geo2015-0297.1

Priezzhev I.I., Egorov S.V. Gibridnoe obuchenie neyronnykh setey s tsel'yu prognoza parametrov neftegazovoy produktivnosti gornykh porod [Hybrid training of neural networks to predict the parameters of petroleum productivity]. Seysmicheskie tekhnologii, 2017, pp. 205-208.

Priezzhev I.I., Egorov S.V. Prognoz kubov uprugikh svoystv po dannym seysmorazvedki 3D i GIS pri pomoshchi algoritma «sluchaynogo lesa» [Spatial rock elastic characteristics prediction using 3D seismic and well logs data based on random forest algorithm]. Geofizika, 2018, no.2, pp. 10-16.

Priezzhev I.I., Egorov S.V., Gladkov E.A. Primenenie algoritmov mashinnogo obucheniya dlya resheniya zadach kolichestvennogo prognoza FES po seysmicheskim i skvazhinnym dannym [Quantitative prediction of reservoir properties using 3D seismic and well data based on machine learning algorithms]. Geofizika, 2018, no.3, pp.33-38.

Schuster G.T., Chen Yu. Seismic Inversion by Newtonian Machine Learning // Research Gate GEO-2019-0434.R2: Preprint (PDF Available) April 2019, 57 р. - https://www.researchgate.net/publication/332630844_Seismic_Inversion_by_Newtonian_Machine_ Learning

Tikhonov A.N., Arsenin V.Y. Solutions of ill-posed problems. - Washington: V.H. Winston and Sons, 1977, 259 p.

© Егоров С.В., Приезжев И.И., 2020


Возврат к списку