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

И.И. Кубышта, А.Ю. Колцун, Научно-Технический Центр «Газпром нефти» (ООО «Газпромнефть НТЦ») г. Санкт-Петербург, В.С. Дакиров, Научно-Технический Центр «Газпром нефти» (ООО «Газпромнефть НТЦ») г. Тюмень, А.В. Малюгин, И.В. Желудков, ООО «Газпромнефть-Ямал».

Источник: Журнал «Геофизика»

ВВЕДЕНИЕ

 Инверсионные технологии преобразования сейсмических данных 3D в кубы упругих параметров широко используются при динамической интерпретации и уже давно вошли в граф стандартных интерпретационных процедур. В ряде случаев выполнение этого пункта технического задания сводится к формальному получению некондиционных кубов, которые затем не используются при построении сейсмогеологических моделей. Существует множество возможных причин, объясняющих некорректность полученного результата. Наряду с отсутствием необходимых каротажных данных и некорректной их подготовкой, отсутствием специализированного программного обеспечения и опыта работы с ним, объективным ограничением является недостаточное качество сейсмических данных. Основой для выполнения акустической инверсии может послужить практически любой полнократный сейсмокуб, независимо от примененных агрессивных постсуммарных процедур. Современные сейсмоданные 3D высокой кратности с большими удалениями источник – приемник, позволяющие сформировать представительный набор частично кратных сумм, не часто оказываются в руках у интерпретатора. Обладая таким уникальным набором, можно выполнить целый ряд тестирований и исследовательских работ. В нашем проекте рассматривается один из наиболее часто возникающих вопросов: «Как влияет количество частично кратных сумм (или, другими словами, соответствующие им максимальные диапазоны углов/удалений) на результаты синхронной инверсии?» Для этого выполнена серия расчетов синхронной инверсии с различным набором офсетных сумм – 3, 5, 7, 10 и 12 сумм. 

Исходные сейсмические данные зарегистрированы на одной из площадей Пуровского района ЯНАО Тюменской области с применением вибрационного источника колебаний. Основные параметры системы полевых наблюдений: расстояние между линиями возбуждения и приема – 300 м, шаг пунктов возбуждения и пунктов приема – 50 м, в активной расстановке количество линий – 32 по 240 каналов, общее количество пунктов приема в активной расстановке – 7680, максимальное удаление взрыв – прием – 7648,6 м. Номинальная кратность наблюдений составляет 320. При таких параметрах система наблюдений обеспечивает регистрацию отраженных волн от целевого ачимовского интервала на углах до 45–50°. Вопросы проектирования, выполнения полевых работ и обработки такой уникальной для Западной Сибири съемки могут послужить темой самостоятельного исследования и в настоящей статье не рассматриваются. Необходимо отметить, что система наблюдений спроектирована с применением технологии лучевого моделирования. Обработка выполнена с сохранением истинного соотношения амплитуд в полном диапазоне удалений, использована временная анизотропная миграция до суммирования. Все ключевые процедуры обработки выполнены под интерпретационным контролем с калибровкой на скважинные данные.

Объект исследования

Разрез рассматриваемой площади характеризуется широким интервалом нефте-, газоносности, первоочередным целевым объектом являются пласты ачимовской толщи – Ач5 1 и Ач5 2-3. Строение пласта Ач5 очень сложное, в его составе выделяются несколько песчаных пластов, к которым приурочены нефтяные залежи, залегающие на различных гипсометрических уровнях. 

В пределах рассматриваемой площади расположены пять скважин, четыре (A, B, D, E) из которых вскрыли разрез ачимовской толщи и одна (скважина С) – разрез аномального строения баженовской свиты. Необходимые для инверсионных преобразований каротажные диаграммы акустического (АК), широкополосного акустического (АКш) и плотностного (ГГК-п) каротажа зарегистрированы только в скважине В, для остальных скважин необходимые кривые были синтезированы в ходе соответствующих работ по петрофизической интерпретации. Суммарная мощность пластов Ач5 1 и Ач5 2-3 по скважинам составляет 40–83 м, суммарные эффективные толщины изменяются от 31,7 до 46,6 м. Глубина залегания пластов ∼ 4000 м.

Обработанные сейсмические данные 3D характеризуются высоким качеством (соотношение сигнал – шум SNR в среднем составляет 8) и ограниченным частотным диапазоном на уровне целевых пластов – не более 30–35 Гц (рис. 1). В сейсмогеологических условиях площади на уровне целевых пластов длина волны λ составляет 115–135 м. Вертикальная разрешенность сейсморазведки, оцениваемая как ¼ λ, не превышает 28–33 м. Целевые пласты Ач5 1 и Ач5 2-3 разделены глинистой перемычкой мощностью 3–6 м, и их границы не формируют самостоятельных отраженных волн. Коэффициент отражения, оцениваемый для кровли пласта Ач5 , не превышает -0,08, что отвечает слабоконтрастным средам. В условиях слабоконтрастных сред, когда коэффициент отражения от кровли перспективного резервуара < 0,1, трехчленная аппроксимация Аки – Ричардса полных уравнений Цеппритца дает незначительную погрешность для распределения коэффициентов отражения продольных волн (Р-волн) в зависимости от угла падения на границу [2]. В такой ситуации отклонения коэффициента отражения от его значения для нормального падения определяются величиной соотношения скоростей продольных и поперечных волн Vp/Vs на границе контактирующих сред [1]. В этом несложно убедиться, построив график изменения коэффициентов отражения Р-волн в зависимости от угла падения по формуле трехчленной аппроксимации Аки – Ричардса (1).


где θ – угол падения, Vp = (Vp1 +Vp2 ) / 2, Vs = (Vs1 +Vs2 ) / 2), ρ = (ρ1 + ρ2 ) / 2, Ip = Vp * ρ, Is = Vs * ρ, ΔIp = Ip2 – Ip1 , ΔIs = Is1 – Is2 , Δρ = ρ1 – ρ2 .


Рис. 1 Фрагмент временного сейсмического разреза в целевом интервале с положением коллекторов по данным ГИС 

Для этого модель целевого пласта Ач5 по данным скважины В зададим следующими упругими параметрами: для покрывающего слоя скорость продольной волны Vp1 = 4250 м/с, скорость поперечной волны Vs1 = 2450 м/с, плотность ρ1 = 2,52 г/см3 ; для подстилающего слоя (целевого пласта Ач5 ) скорость продольной волны Vp2 = 3800 м/с, скорость поперечной волны Vs2 = 2050 м/с, плотность ρ2 = 2,4 г/см3 . Полученный график AVO-характеристики показан на рис. 2 голубым цветом (R): наблюдается незначительное уменьшение общего коэффициента отражения (Rpp) с увеличением удаления. Кривыми а, b и c показано поведение каждого из многочленов уравнения (1). Традиционно многочлен a приблизительно равен коэффициенту отражения продольной волны для нормального падения. Начиная с углов 8–10° влияние многочлена b, зависящего от скоростей поперечных волн, начинает расти и в условиях нашего проекта достигает значительной доли от общего Rpp на углах падения 40–45°. Влияние третьего многочлена c, наиболее зависимого от плотностных характеристик среды, остается незначительным во всем диапазоне удалений. Другими словами, в сейсмогеологических условиях изучаемого месторождения в результате синхронной инверсии существуют предпосылки ожидать хорошее восстановление таких упругих параметров, как поперечный импеданс и отношение скоростей Vp/Vs, предпосылки к восстановлению плотности ограниченны.

Методологические подходы и результаты синхронной инверсии

В алгоритме синхронной инверсии частично кратных сумм одновременно используются несколько угловых или ограниченных по выносу (офсетных) суммарных сейсмических кубов для восстановления упругого продольного, поперечного импедансов, отношения Vp/Vs и плотности [1, 3]. 

Выбор количества частично кратных кубов и соответствующих им диапазонов углов/удалений – задача, напрямую зависящая от номинальной кратности съемки и качества обработанных сейсмограмм. В то же время разделенные суммы должны характеризоваться высоким SNR и обеспечивать хорошую сходимость фактических AVO-характеристик с теоретическими. Способ формирования сумм – в углах или в офсетах – по большей части определяется используемым для инверсии программным обеспечением (ПО). Традиционно предпочтение отдаетсяугловым диапазонам. Используемое при выполнении проекта ПО позволяет работать с офсетными суммами, для которых дополнительно необходима глубинно-скоростная модель, удовлетворяющая требованиям привязки скважин и отражающая геологическое строение площади. Выбор в пользу офсетов отдан по причине недостаточной сходимости скоростной модели, полученной после временной миграции, со скоростными характеристиками среды по данным скважин. Шаг офсетных коридоров выбран переменным: исключены сейсмограммы из ближней зоны до 700 м, сумма № 1 включает интервал 700–1500 м, что приближенно отвечает диапазону отражений для нормального падения, диапазон до 6250 м разбит с шагом 400–500 м на 10 сумм и заключительная сумма собрана в диапазоне 6250–7680 м. Такой широкий диапазон суммирования объясняется низкой кратностью на дальних удалениях и непротяженным по времени интервалом наличия отраженных волн. Все частично кратные суммы были подвергнуты необходимому «интерпретационному постпроцессингу», включающему выравнивание и нормировку амплитуд. Важным шагом при выполнении синхронной инверсии, так же как и акустической, является оценка формы и характеристик сейсмического импульса по всем частично кратным суммам. Оценка сигналов выполнена детерминистическим способом по данным четырех скважин (A, B, D, E). Форма всех импульсов близка к симметричной (нуль-фазовой), фазовый спектр стабилен в сейсмической полосе частот. Несмотря на ограниченный временной интервал наличия сейсмоданных, для последней суммы также удалось определить свой импульс. Сейсмические данные имеют ограниченный частотный диапазон (8–12 Гц – 35–40 Гц), поэтому для восстановления отсутствующей низкочастотной компоненты требуется построение низкочастотных или начальных моделей упругих свойств. Данные модели играют роль геологического ограничения в решении инверсии. Выполнение синхронной инверсии возможно в нескольких параметризациях и является зависимым от исходных скважинных данных. Для текущего проекта была выбрана параметризация продольный импеданс (P-impedance)/отношение скоростей Vp/Vs/плотность (Density). Низкочастотные модели указанных упругих свойств были построены по данным всех пяти скважин на стратиграфическом каркасе, алгоритм интерполяции – inverse distance weighted. Подбор параметров инверсии для каждого расчета не имеет принципиальных отличий, низкочастотные модели упругих свойств использованы одинаковые. 

Рис. 2 Коэффициенты отражения Р-волн в зависимости от угла падения на кровлю пласта Ач5 . (Комментарии в тексте.) 

Для решения поставленной задачи выполнены расчеты синхронной инверсии с различными максимальными удалениями сейсмоданных: инверсия с 3 суммами (максимальные удаления 3050 м) далее в тексте условно названа 3ST, с 5 суммами (максимальные удаления 3900 м) – 5ST, с 7 суммами (максимальные удаления 4600 м) – 7ST, с 10 суммами (максимальные удаления 7450 м) – 10ST, с 12 суммами (максимальные удаления 7680 м) – 12ST.

Качество решения всех расчетов проконтролировано традиционными способами: сопоставлением извлеченных в точках скважин псевдокривых с записанными; анализом полученных разрезов петрофизических свойств; сравнением реального сейсмического поля, поля, восстановленного в результате инверсии, и их разницы; анализом рассчитанных автоматически карт коэффициента корреляции сейсмики и синтетики, карт инвертированного отношения сигнал – шум для каждой частичной суммы. Все варианты расчета показывают одинаково высокое качество восстановления сейсмических данных. Наибольшее отличие, как и ожидалось, наблюдается в восстановлении таких упругих параметров, как отношение скоростей Vp/Vs и плотность (рис. 3). С увеличением количества частично кратных сумм повышается когерентность разрезов упругих свойств, упрощая геологическую интерпретацию полученных материалов. Несмотря на ограниченные предпосылки, в результате инверсии 10ST и 12ST удалось получить кубы плотности, отражающие основные особенности строения Ачимовской толщи. Рассмотренный ранее расчет AVO-характеристики для отражающего горизонта Ач5 (рис. 2) можно назвать «худшим случаем», остальные геологические границы обладают более высокими Rpp, и для них изменения, вызванные влиянием плотности, начинают фиксироваться на меньших углах падения.

Рис. 3 Разрезы отношения Vp/Vs (a – e) и плотности (f – j) по составной линии через скважины Е – D – B Разрезы а и f – результат инверсии с тремя суммами, b и g – с пятью, c и h – с семью, d и i – с десятью, e и j – с двенадцатью. Разрезы упругих параметров показаны в едином динамическом диапазоне цветовой шкалы. Сплошными линиями черного цвета вынесена корреляция отражающих горизонтов 

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

По данным скважины В в ходе работ по петрофизической интерпретации для ачимовских пластов выявлены взаимосвязи упругих параметров с литологическим составом, пористостью и насыщением. Установлено, что в поле значений акустического импеданса наблюдается разделение коллектор – неколлектор. В терминах математической статистики акустический импеданс – случайная величина, распределенная по нормальному закону, и для нее выполняется правило трех сигм: модуль ее отклонения от среднего значения не превосходит трех стандартных отклонений (σ). Если рассмотреть значения акустического импеданса в интервале пластов Ач5 1 и Ач5 2-3 для «литотипов» коллектор и неколлектор раздельно и построить графики плотности вероятности (PDF – Probability Density Function), можно установить, что PDF «литотипов» коллектор и неколлектор пересекаются в точке +/- 1,2σ. Другимисловами, 75% площади PDF коллектора не пересекается со значениями акустического импеданса «литотипа» неколлектор. Возможен площадной прогноз зон коллекторов на основе куба акустического импеданса. В качестве граничной отсечки примем более жесткое значение, составляющее среднее + 1σ и численно равное 9750 г/см3 *м/с. Для прогноза зон коллекторов в интервале пластов Ач5 1 и Ач5 2-3 выполнена объемная интерпретация по кубу акустического импеданса, полученного в результате синхронной инверсии для каждого из вариантов расчета (3ST, 5ST, 7ST, 10ST и 12ST). Не совсем корректно называть результат синхронной инверсии именно «акустическим импедансом» (более подходящее название – упругий импеданс продольных волн), но для упрощения терминологии далее в тексте будет использовано такое название. Как уже отмечалось, куб акустического импеданса можно получить и в результате более простой и менее ресурсоемкой акустической инверсии полнократного куба, что и было сделано для сравнения. Акустическая инверсия выполнена с той же низкочастотной моделью, которая была использована для серии расчетов упругой инверсии. Сейсмический импульс также определен детерминистическим способом по четырем скважинам. 

В результате объемной интерпретации по кубам акустического импеданса по граничной отсечке (7,8–9,75) * 1000 г/м3 *м/с получен объемный прогноз временных толщин коллекторов пластов Ач5 1 и Ач5 2-3. Через интервальную скорость прогноз переведен в глубины и верифицирован на суммарные эффективные толщины (Нэф) в скважинах. При таком подходе все скважины можно считать проверочными для полученного прогноза Нэф. Стоит акцентировать внимание, что прогноз Нэф получен без использования каких-либо регрессий «значение Нэф в скважине – значение сейсмического атрибута»: фактически прогнозная карта Нэф является двумерным представлением объема среды, обладающего определенными значениями акустического импеданса. В условиях изучаемой площади, когда целевые пласты вскрыты всего четырьмя скважинами, такой подход является предпочтительным. Результаты верификации сведены в таблицу 1. 

Как видно, результат синхронной инверсии дает более точный прогноз Нэф, чем акустической инверсии. В прогнозах по результатам различных вариантов синхронной инверсии более точный получен для инверсии 12ST. Средняя ошибка прогноза Нэф по результату инверсии 12ST составила минус 1,13 м, 10ST – минус 1,77 м, 7ST – минус 2,65 м, 5ST – минус 4,23 м, 3ST – минус 5,91 м. Знак «минус» означает, что прогноз занижен относительно скважинных значений. Средняя ошибка прогноза Нэф по результату акустической инверсии составляет 8,06 м, прогноз завышен относительно скважинных значений, причем для двух скважин (В и Е), более чем на 40%.

Интересно сравнение прогнозных схем Нэф, полученных по результату акустической и синхронной инверсии (рис. 4). 

Рис. 4 Прогнозные схемы эффективных толщин пластов Ач5 1 и Ач5 2-3, полученные по результатам акустической инверсии (а) и синхронной инверсии 12ST (b) 

Прогнозная схема – результат акустической инверсии показывает область развития повышенных Нэф к западу от скважины D и большую ошибку в скважине Е. Прогнозная схема – результат синхронной инверсии показывает тренд на уменьшение Нэф в западном направлении, согласный со скважинным трендом, и занижает прогноз в скважине D. Аналогичная область значительных отличий располагается к северу от скважины D. Это – области неопределенности, но в то же время прогнозная схема – результат синхронной инверсии, на наш взгляд, дает более скромный и достоверный прогноз. В ходе петрофизического моделирования установлено, что на отношение скоростей Vp/Vs оказывают комплексное влияние и изменение глинистости (повышая значения), и наличие флюида (понижая значения). По гистограммам значений Vp/Vs «литотипов» коллектор и неколлектор для построения прогнозов выбрана граничная отсечка, равная 1,77. Аналогичным описанному ранее способом выполнен объемный прогноз временных толщин улучшенных коллекторов по граничной отсечке Vp/Vs 1,5–1,77. 


Рис. 5 Прогнозные схемы зон улучшенных коллекторских свойств, построенные по результатам пяти вариантов синхронной инверсии 

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

Варианты геологической интерпретации

Перспективные пласты Ач5 1 и Ач5 2-3 по данным скважин представлены средне- и грубозернистыми песчаниками и обладают хорошими фильтрационноемкостными свойствами (ФЕС). Предполагается, что для пластов Ач5 1 , Ач5 2-3 характерно развитие нескольких налегающих друг на друга линзовидных тел (лопастей), контролирующих формирование залежей с разными водонефтяными контактами (ВНК). В скважине В был получен фонтан водонефтяного притока, в скважине А, занимающей более низкое гипсометрическое положение, безводный приток нефти. Подтверждение гипотезы о наличии разделенных литологически экранированных тел (конусов выноса) можно увидеть на прогнозной схеме зон улучшенных коллекторских свойств (рис. 6, b). Предполагаемые конусы выноса показаны заливкой серого цвета, выделены по результатам синхронной инверсии 12ST. Как видно, предполагаемая модель состоит из трех уровней проградирующих в северозападном направлении конусов выноса. Наиболее крупный верхний конус, условно названный «К1», расположен в районе скважины В – на прогнозной схеме видна сама лопасть конуса и распределительные каналы. По периферии этого песчаного тела к западу, северо-западу и северу выделяются конусы «К2», «К3» и «К4». Конусы «К2» и «К3» полностью располагаются в пределах сейсмической съемки, и вскрывшие их скважины А и D характеризуются различным типом насыщения (А – приток безводной нефти, D – приток нефти с водой). 

Скважина D расположена гипсометрически ниже скважины А и не противоречит гипотезе о гидродинамической связи конусов «К2» и «К3» или их объединению в один с несколькими подводящими каналами. К северу в непосредственной близости от границы сейсмической съемки расположена скважина, в которой из пласта Ач5 2-3 получен приток безводной нефти. 

Предположительно, она характеризует обособляющийся северный конус «К4», часть лопасти которого и мощный подводящий канал видны на прогнозной схеме на рис. 6, b. Скважина Е, расположенная в северо-западной части площади, скорее всего, также вскрывает отложения обособленного песчаного тела. Из-за значительно уменьшенных мощностей отложений вдали от источников сноса оконтуривание этого предполагаемого конуса третьего уровня («К5») затруднено.

Прогнозная схема эффективных толщин – результат акустической инверсии (рис. 6, а) не позволяет районировать территорию и выделять обособленные литологические объекты.

Рис. 6 Выделение линзовидных тел (конусов выноса) в интервале пласта Ач5 Цветом показаны схемы прогнозных эффективных толщин пластов Ач5 1 и Ач5 2-3. На прогнозную схему b) сплошной заливкой серого цвета вынесены зоны
улучшенных коллекторских свойств

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

ЗАКЛЮЧЕНИЕ

В статье подробно представлены результаты детерминистических инверсий – акустической полнократного сейсмокуба 3D высокой кратности и синхронной, выполненной с различным набором частично кратных сумм. На основе петрофизической интерпретации скважинных данных обоснована возможность прогноза эффективных толщин по абсолютным значениям акустического импеданса и возможность прогноза зон улучшенных коллекторских свойств по абсолютным значениям второго упругого параметра – отношения скоростей Vp/Vs. Количественно обоснованные граничные отсечки упругих параметров легли в основу объемной интерпретации временных кубов – результатов инверсии. Через интервальную скорость прогноз переведен в глубины и верифицирован на суммарные Нэф пласта Ач5 в скважинах. Ошибка прогноза Нэф в точках скважин приравнена к оценке достоверности количественного прогноза по результатам инверсии. Установлено, что результат синхронной инверсии дает более точный прогноз Нэф, чем акустической инверсии, в среднем на 15%. В прогнозах по результатам различных вариантов синхронной инверсии более точный получен для инверсии с максимальным количеством частично кратных кубов и максимальными удалениями. Абсолютная ошибка прогноза Нэф для различных вариантов синхронной инверсии изменяется от 6 до 15%. 

Привлечение к интерпретации второго упругого параметра – отношения скоростей Vp/Vs – позволяет оконтурить области улучшенных коллекторских свойств, геометрический облик которых схож с конусами выноса, их лопастями и подводящими/ распределительными (?) каналами. Выделенные литологические тела хорошо увязываются с геологической концепцией строения целевого интервала, сформированной на основе скважинных данных. 

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

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

ЛИТЕРАТУРА

1. Воскресенский Ю.Н. Изучение изменений амплитуд сейсмических отражений для поисков и разведки залежей углеводородов. М.: РГУ нефти и газа, 2001. 68 с.
2. Haase A.B., 2004. Modelling of linearized Zoeppritz approximations: CREWES Research Report, 16.
3. Rutherford S.R., Williams R.H. Amplitude-versus-offset variations in gas sands. Geophysics, 1989. V. 54. № 6. Pp. 680–688.

REFERENCES

1. Voskresenskij YuN. Izuchenie izmenenij amplitud sejsmicheskih otrazhenij dlya poiskov i razvedki zalezhej uglevodorodov. Moscow: RGU nefti i gaza, 2001. 68 p. (in Russian).
2. Haase A.B. 2004. Modelling of linearized Zoeppritz approximations: CREWES Research Report, 16.
3. Rutherford SR, Williams RH. Amplitude-versus-offset variations in gas sands. Geophysics, 1989; 54 (6): 680–688.

Положительная рецензия от 01.08.2018
Решение редколлегии о публикации от 31.08.2018

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