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

УДК 519.876.5:550.8

А.Г. Черников, М.Б. Матушкин, А.Ю. Дегтерёв

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

Ключевые слова:

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

Keywords:

three-dimensional,

modelling,

geology,

petrophysics,

interpolation,

fuzzy,

Markov sequence,

recognition,

forecast,

filtration-capacitive properties.

В настоящее время результативные решения по всем аспектам разведки и эксплуатации объектов геотехнического комплекса, включая действующие и проектируемые горные предприятия, принимаются с учетом данных их трехмерного моделирования. В этой связи корректность 3D-моделей и их адекватность реальным геологическим объектам имеют большое практическое значение. Как правило, фактического материала, полученного в результате непосредственного изучения геологических объектов, недостаточно для того, чтобы составить полное представление об их внутреннем строении, состоянии и свойствах. Поэтому адекватность геологических моделей существенно зависит не только от исходных данных, но и от численного метода, использованного при построении модели. Для решения задачи трехмерного геологического моделирования в современных программных средствах используются в основном численные методы, основанные на предположениях статистического характера об изменчивости переменных, например метод Криге (кригинг), последовательная гауссова симуляция. Данные алгоритмы имеют ряд ограничений, не во всех случаях позволяя автоматически получать геологически достоверное распределение свойств. В связи с этим производители программных средств геологического моделирования, в частности таких распространенных, как Schlumberger Petrel, Roxar RMS и др., пополнили свои продукты специальными инструментами моделирования, дающими исследователю больше возможностей для внесения в модель субъективных параметров, основанных на субъективном понимании закономерностей строения моделируемого объекта. При условии работы высококвалифицированного геолога-моделиста это позволяет скорректировать работу автоматизированных методов, устранив присущие им недостатки, однако принятие ошибочных решений, например из-за усталости или недостаточной квалификации оператора, может привести к еще большим ошибкам моделирования.

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

Одним из возможных подходов при решении задач, связанных с моделированием систем, характеризующихся нелинейной изменчивостью свойств, является

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

При изучении последовательности геотехнических событий можно рассматривать любые интервалы в качестве «прошлого» (/ - 1) и «настоящего» (/) и тем самым воспользоваться марковским свойством для предсказания состояния системы на основе изменения состояний во времени (в разрезе) ее элементов. Аналогичный подход можно использовать для любых чередований событий, определенных на какой-либо шкале последовательностей. В качестве такой последовательности событий можно рассматривать множество состояний каких-либо переменных, например последовательность чередования типов пород в терри-генном разрезе (Ь1, Ь2, Ь3, Ь4 ... Ы) [2].

С позиций системологии [3] марковская последовательность представляет собой реализацию объекта моделирования в виде списка наблюденных переменных на базе, выраженной через абстрактное «время». Переменные выражены через их определенные (формализованные) значения. Последовательность переменных представляет собой вектор, идентифицированный меткой - параметром. Формирование векторов переменных целесообразно выполнять в виде графов. Графы подразделяются на простые и составные. Простые графы бывают двух видов: открытые и замкнутые. Открытые графы представляют собой нормальную временную последовательность переменных. Примером такой последовательности служит порядок залегания литотипов в стратиграфическом разрезе.

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

графа служит граф, условно названный «транс-позитным» (от лат. йашроБШо - перестановка). Переменные, составляющие вектор этого графа, образуют полный комплекс сочетаний с перестановкой. Используется он в случае включения в вектор переменных, не упорядоченных на какой-либо шкале (переменная с номинальной шкалой). Например, координаты пространства. Число подобных сочетаний определяется по формуле N = п2 - п, где N - суммарное число переменных в транспозитном векторе, а п -число исходных переменных, включенных в вектор. Использование последовательностей этого типа обеспечивает однозначное позиционирование для списка переменных с номинальной шкалой [2].

Существенным моментом формализации (абстрагирования) для параметров и переменных с целью использования их математических свойств и отношений является операция, вводящая переменную (параметр) как образ свойства. Канал наблюдения о,, с помощью которого свойство Л1 представляется переменной V, реализуется функцией о,: Л1 ^ V,. Эта функция гомоморфна относительно предполагаемых свойств множеств Л1 и V. Для некоторых параметров и переменных каналы наблюдения могут представлять собой явно заданные функции о,, для других функция задается исследователем. Например, можно на упорядоченном множестве, представленном координатной шкалой, выраженной в географических градусах или километрах, выделить некоторое количество классов, задав их граничными значениями классификации. Возможно разделение, когда множество значений будет разбито на группы, соответствующие подмножествам, лимитированным границами при равновероятном распределении. Существуют и другие приемы введения образа свойств. Особенностью разделения на классы являются ограничения, накладываемые на количество классов величиной погрешности наблюдения, объемом и плотностью распределения на шкале состояний переменной.

Далее переменной присваивается код - целочисленное значение в интервале (I... ?), обозначающее как ее метку, так и ее состояние. Следующей переменной присваивается значение в интервале + I... п) и т.д. Таким образом, меткой становится интервал, а состоянием переменной - число в этом интервале.

Приведенные рассуждения касались в основном простых переменных. Однако

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

Для сжатого описания поведения марковской цепи используются матрицы условных вероятностей перехода. Каждый элемент этой матрицы представляет собой вероятность перехода из заданного состояния (этому состоянию в матрице соответствует строка) в следующее состояние (которому соответствует столбец). Таким образом, относительно данного множества состояний все возможные переходы предусмотрены. Матрицу вероятностей переходов ||Р:||, соответствующую п состояниям, можно записать следующим образом:

Р, =

Р,

Р.,

Р.

Р..

где р - вероятность перехода из состояний Б: 51, 52, ..., Бп, соответствующих строкам матрицы, в состояния Б,: Бх, 52, ..., Бп, соответствующие столбцам. Таким образом, вероятности перехода из состояния Б1 в состояния Б1, 52, ..., Бп записаны в верхней строке матрицы Р., из состояния Б2 - во второй строке и т.д. То есть г-я строка матрицы вероятностей переходов показывает вероятности переходов из г -го состояния в любое возможное состояние, а ,-й столбец - вероятности переходов из любого возможного состояния в,-е. При г = ] вероятность показывает продолжительность нахождения переменной в состоянии г. Сумма вероятностей

в строке равна 1. Операционно матрица носит название «матрицы переходных вероятностей», и ей присваивается обозначение О.

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

Математическая постановка задачи распознавания сводится к следующему. Пусть задано множество М состояний системы, на этом множестве существует разбиение на конечное число подмножеств (классов) й.: М = и(й,), г = 1 ... т. Разбиение определено не полностью. Задана лишь некоторая информация /о о классах й,.. Объекты ю задаются значениями признаков Х,} = 1, ..., N (этот набор всегда один и тот же для всех объектов, рассматриваемых при решении определенной задачи). Совокупность значений признаков Х] определяет описание /(ю) объекта. В качестве образа признаков Х, принимается матрица переходных вероятностей О, (таблица). Последовательность формирования «учителя»: описание объекта /(ю) = А,(ю) ^ У(ю) ^ О, является стандартным и принимает значение из множества допустимых значений.

Задача распознавания со стандартной информацией состоит в том, чтобы для данного объекта ю' и набора классов йь ..., йт по обучающей информации /„(й^ ..., йт) о классах и описанию /(ю) вычислить значения вероятности отнесения объекта к одному

Концептуальная модель обучения системы распознавания с «учителем» при марковском прогнозировании

Объект Матрица признаков Класс

ю1 Юг О.

Юг+1 юк [О,}к йк

юк+1 Ют {О,}к+1 {О,}т Йт

ю' О й?

Рис. 1. Блок-схема марковской системы прогнозирования УесРг^п

из классов ^ т. Полученные значения Р(Ц^) являются гипотезами. В качестве решения принимаются гипотезы с максимальными значениями вероятности: ю = тах Р(йк).

Приведенный алгоритм программно реализован в системе марковского прогнозирования VecProgn (рис. 1).

С использованием в качестве алгоритмической базы представленной ранее системы марковского прогнозирования разработана принципиально новая система 2Б-, 3Б-, 4Б-моделирования сложнопостроенных геологических объектов. В основу алгоритма марковского 2-, 3-, 4-мерного моделирования положен подход, предполагающий удержание в памяти не координатных точек пространства со значениями соответствующих им геолого-геофизических показателей, а вероятностного образа свойств геологических тел на Ы-мерной решетке. Создание вероятностного образа осуществляется на совокупности данных, в которой координаты точек рассматриваются как координаты фазового пространства, последовательность числовых значений которых представлена в виде марковской последовательности со строго определенными количеством и чередованием шагов. В отличие от представленного ранее алгоритма система использует нечеткие каналы наблюдения с ветвящимися графами марковских состояний: переход из нечеткого состояния (/' - 1) в нечеткое состояние (?) представляется в виде ветвящегося графа (рис. 2).

Суммарная вероятность возможных переходов равна 1 (см. ранее). Принципиальное

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

Последовательным накоплением векторов (последовательностей значений) координат получаются матрицы накопленных частот, отнесенные к конкретному значению или интервалу значений (классу) моделируемого геологического признака. Дальнейшее преобразование матриц частот в матрицы переходных вероятностей и создает модель распределения признака в координатах Ы-мерного пространства. Распознавание ведется по точкам фазового (декартового) пространства. Генерация координат точек осуществляется по выбранным

12,6

10 | Р = 0,74:

20 | Р = 0,26:

„90 | Р = 0,80 '100 | Р = 0,20 ,90 | Р = 0,80 '100 | Р = 0,20

Рис. 2. Пример ветвящегося графа нечетких марковских состояний

исследователем плоскостям сечения. Задаются граничные точки координат и шаг дискретизации. На выходе модели получается вероятностная оценка распределения геолого-геофизических характеристик модели на плоскости сечения массива (3Б-модель) или последовательность сечений в выбранном временном интервале (4Б-модель). Визуализация прогноза осуществляется на экране дисплея в двухмерных координатах в виде разрезов моделируемого геологического объекта.

Таким образом, прогнозирование подразумевает два этапа:

1) моделирование, заключающееся в создании вероятностного образа свойств объекта, выраженного в условных единицах нечеткой привязки к границам классов, определенных на координатных шкалах;

2) распознавание, заключающееся в том, что по структурированной совокупности пространственных координат дается вероятностная оценка параметра в любой точке модели, принимаемая в качестве вероятностного

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

Описанный алгоритм использован в программе марковского 4Б-моделирования и опробован при решении широкого спектра задач геологического моделирования [4-9], в том числе успешно применялся для решения задач ЗБ-моделирования геологического строения подземных хранилищ газа (ПХГ) [10] (рис. 3).

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

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

0 2000 4000 6000 8000 10000 12000 14000 16000 18000 20000 22000 24000 26000

Расстояние по профилю, м

Рис. 3. Вертикальная томограмма нечеткой марковской ЭБ-модели ПХГ

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

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

а это в свою очередь приводит к снижению качества прогноза. Для устранения влияния перечисленных факторов вводится технология нормирования данных на геофизическую нормальную модель исследуемого разреза [6]. Получение нормальной геофизической модели разреза состоит из последовательно выполняемых в автоматическом или интерактивном режимах операций. Программно технология реализована следующим образом: значения глу -бин по всем скважинам, объединенным в пакет, нормируются на стандартный интервал глубин (АИнорм) и ранжируются по возрастанию нормализованной глубины. В результате получают осложненный шумами интегральный образ пласта по ГИС (рис. 4).

Для получения нормального образа исследуемого объекта проводится фильтрация интегральной кривой с подавлением высокочастотной составляющей по интервалам АИнорм. Результатом является модельная диаграмма исследуемого объекта. На рис. 4 она для наглядности наложена на интегральную кривую.

Грамши ПфЫ

. < 7<2 733 .'. :> 778 311 7*0 760 73Й 750 73* 7*9 335 797 733 757 775 751 775 305 77* Е52 800 737 317 737 ЯШ 73* 332 Е85 813 303 793 30* 783 3*3 305 339 813 332 333 327 833 333 3*1 3*3 3*3

[1ЕРТ

Рис. 4. Просуммированная (интегральная) кривая потенциалов самопроизвольной поляризации (ПС), нормализованных по глубине и амплитуде (синяя), и ее сопоставление с отфильтрованной модельной кривой ПС (желтая) в программе УесРгс^

— исходная — модельная — нормированная на модельную кривую

Рис. 5. Сопоставление кривых ПС в пакете до и после приведения к нормализованному виду в программе УесРго^

— Кш — стратиграфические разбивки — нормированная — модельная

Рис. 6. Фрагменты пакетной петрофизической модели Коп, полученной по результатам

многофакторного прогнозирования

реальные скважины + виртуальные скважины — площадь подсчета

Рис. 7. Схема расположения реальных и виртуальных скважин, использовавшихся для построения ЭБ-модели ПХГ и расчета порового объема пласта-коллектора

K — стратиграфические разбивки

Рис. 8. Прогнозные кривые пористости, рассчитанные для виртуальных скважин

На первом этапе моделирования выполняется приведение всех кривых геофизических исследований в скважинах на ПХГ к стандартным значениям (рис. 5). По нормированным геофизическим данным с привлечением комплекса геолого-геофизических характеристик в системе УесР1^п делается многофакторный прогноз фильтрационно-емкостных свойств (ФЕС) моделируемого геологического объекта. Результат в виде прогнозных значений открытой пористости Коп приведен на рис. 6. Аналогичным образом устанавливаются зависимости для других показателей.

Для расчета 3Б-модели распределения ФЕС необходимо компенсировать неоднородность в распределении скважин по площади ПХГ. Для этой цели был применен способ многофакторного прогнозирования свойств с помощью системы нелинейной марковской статистики. Для получения 3Б-модели распределения пористости сгенерирована равномерная сеть виртуальных скважин (рис. 7).

В виртуальных скважинах с шагом 250 м по горизонтали и 0,5 м по вертикали на основании описанной ранее петрофизической модели и с учетом латеральной изменчивости рассчитаны значения Коп (рис. 8). Расчет пористости в виртуальных скважинах осуществлялся системой распознавания на основании степени близости к классам обучающей выборки в многофазовом пространстве признаков. Помимо координат в пространство признаков входили показания геофизических методов,

3 5 7 9 11 15 19 23 27

I K ,%

Рис. 9. Фронтальный срез (томограмма) петрофизической ЭБ-модели продуктивного объекта ПХГ, сопоставленный с рельефом кровли пласта-коллектора

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

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

7 = 7100м

N -200

269000 X, м Г= 7000 м

263000

264000

265000

266000

267000

268000

-100

-150

-200

-250

269000 X, м Г = 6900 м

I

263000

264000 265000 266000 267000

268000

N -100

-150

-200

-250

269000 X, м 7 = 6800м

28 -26 -24 -22 -20 - 18 - 16

- 14

- 12 - 10 - 8 -6

4

К , %

263000

264000

265000

266000

267000

268000

269000 X, м Г= 6700 м

^ -150

-200

-250

-300

263000 264000 265000 266000 267000 268000 269000 Хм Рис. 10. Вертикальные разрезы строения порового пространства продуктивной толщи

взяты рельеф кровли верхнего контрольного горизонта (рис. 9) и результаты прогноза значений пористости в реальных и виртуальных скважинах (см. рис. 6, 8).

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

Список литературы

1. Вистелиус А. Б. Основы математической геологии / А.Б. Вистелиус. - Л.: Наука, 1980.

2. Гриб Н.Н. Методологические основы системного исследования массива горных пород / Н.Н. Гриб, А.В. Самохин,

А.Г. Черников. - Якутск: Изд-во ЯНЦ СО РАН, 2000. - 104 с.

3. Клир Дж. Системология. Автоматизация решения системных задач: пер. с англ. / Дж. Клир. - М.: Радио и связь, 1990. -544 с., ил.

4. Гриб Н.Н. Технология прогнозирования пространственной изменчивости показателей качества углей / Н.Н. Гриб, А.В. Самохин, А.Г. Черников // Наука и образование. - Якутск: АН Республики Саха (Якутия), 2001. - № 3. -С. 131-139.

5. Гриб Н.Н. Пространственное моделирование свойств угольных пластов сложного строения / Н.Н. Гриб, А.Г. Черников // Горный информационно-аналитический бюллетень. -2003. - № 8. - С. 45-46.

6. Матушкин М.Б. Пакетная обработка геолого-геофизических данных для создания геолого-технологической модели / М.Б. Матушкин, А.Г. Черников // Горный информационно-аналитический бюллетень. - 2011. - № 12. -С. 115-124.

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

7. Черников А.Г. Методика прогнозирования литогенетических свойств углей

по комплексу геолого-геофизических данных / А.Г. Черников // Сб. материалов IX Всесоюзного угольного совещания. Кн. 2. -Ростов н/Д, 1991. - С. 72-74.

8. Черников А.Г. Использование марковской гипсотомографии при геологических исследованиях в океанологии / А.Г. Черников, Н.В. Либина // Океанология. - 2011. - Т. 51. -№ 3. - С. 561-565.

9. Черников А.Г. Оценка возможностей марковской гипсотомографии

при моделировании глубинного строения нефтегазоносных акваторий / А.Г. Черников, Н.В. Либина // Вестник ЦКР Роснедра. -2011. - № 6.

10. Исхаков А.Я. Моделирование изменчивости свойств породного массива на основе нечетких марковских последовательностей / А.Я. Исхаков, М.Б. Матушкин, Р.Г. Темиргалеев и др. // Подземное хранение газа. Полвека

в России: опыт и перспективы: сб. науч. тр. -М.: ВНИИГАЗ, 2008. - С. 255-265.

11. Дегтерёв А. Ю. Перспективы использования интерполяционного предсказательного моделирования при решении задач эксплуатации подземных хранилищ газа / А.Ю. Дегтерёв, В.А. Алькин, А.Г. Черников и др. // Газовая промышленность. - 2015. -№ 8. - С. 75-79.