ИССЛЕДОВАНИЕ ОСОБЕННОСТЕЙ ГЕОЛОГИЧЕСКИ СОГЛАСОВАННОЙ АДАПТАЦИИ ГИДРОДИНАМИЧЕСКИХ МОДЕЛЕЙ. ЧАСТЬ 1. БЕЗГРАДИЕНТНЫЕ МЕТОДЫ ОПТИМИЗАЦИИ Текст научной статьи по специальности «Науки о Земле и смежные экологические науки»

Актуальные проблемы нефти и газа ■ Вып. 2(21) 2018 ■ http://oilgasjournal.ru

ИССЛЕДОВАНИЕ ОСОБЕННОСТЕЙ ГЕОЛОГИЧЕСКИ СОГЛАСОВАННОЙ АДАПТАЦИИ ГИДРОДИНАМИЧЕСКИХ МОДЕЛЕЙ.

ЧАСТЬ 1. БЕЗГРАДИЕНТНЫЕ МЕТОДЫ ОПТИМИЗАЦИИ

И.М. Ширяев, И.М. Индрупский, Э.С. Закиров ИПНГ РАН, e-mail: shiryaev.ivn@gmail.com

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

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

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

В следующем разделе будет введено понятие геологически согласованной адаптации и дано определение основного инструмента геостатистики - вариограммы. Далее будет представлена постановка обратной задачи и методы оптимизации для ее

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

Геологическая согласованность

В предлагаемом подходе сохранение заложенных в модель геологических принципов при адаптации достигается за счет уточнения не значений ФЕС в ячейках пласта, а параметров вариограмм, используемых для построения распределений ФЕС в модели. Вариограмма является количественной характеристикой неопределенности пространственно распределенных данных и представляет собой дисперсию разницы значений параметра в двух точках как функцию расстояния и направления между ними [5, 6]. Значения вариограммы возрастают от минимального в нуле, обусловленного только дисперсией ошибки измерений, до максимального на расстоянии, на котором взаимная корреляция значений параметра в двух точках перестает прослеживаться. Для описания различия в пространственном распределении данных вдоль разных направлений используют анизотропную вариограмму. Основные (и обычно наиболее неопределенные) параметры анизотропной вариограммы - два перпендикулярных радиуса в горизонтальной плоскости модели (в главном направлении и во второстепенном направлении й2), радиус вдоль вертикальной оси й3 и угол поворота главного направления в плоскости относительно одной из координатных осей. При этом главный радиус обычно задается вдоль основного тренда распространения ФЕС в пласте.

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

Вариограмма используется во всех методах так называемой двухточечной геостатистики [5]. Существуют как детерминированные алгоритмы, которые представляют собой интерполяцию статических скважинных данных и позволяют получить одну наиболее вероятную модель распределения ФЕС, так и стохастические методы, которые позволяют получить множество равновероятных моделей (реализаций).

В данной работе использован наиболее известный из детерминированных методов -точечный кригинг. В кригинге неизвестное значение параметра (ФЕС) рассчитывается следующим образом:

N

ЯЫ = .

1=1

(1)

Здесь А^ - веса кригинга, х0 - точка, где оценивается значение свойства, а г(х1),г(х2), ...,г(хм) - известные значения свойства в других точках (статические данные на скважинах). Веса выбираются таким образом, чтобы оценка в статистическом смысле была несмещенной, а дисперсия минимальной. Это приводит к известной системе из N + 1 линейного уравнения - системе обычного кригинга:

(2)

Ун У12 . ■ Уш 1 -У10-

У21 У22 ■ У2М 1 ¿2 У20

Уы1 Ум2 . ■ Уыы 1 Уыо

[ 1 1. 1 0- -V- [ 1-

а дисперсию получаемой оценки (1 ) можно выразить в терминах вариограммы:

а2

= - ^^ ЩУ1] А1Ую - У 00,

(3)

I )

где - значение вариограммы соответственно между I и ] точками замеренных данных, Ум - между точкой данных I и точкой, в которой оценивается значение поля, у00 -величина эффекта самородка [5], д - множитель Лагранжа для условия несмещенности оценки.

Обратная задача

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

п к

1=1 ]=1

и л

(4)

Здесь У[ - замеренные показатели работы скважин (дебиты фаз, обводненность, забойные давления и т.д.), у(и)\ - рассчитанные на модели показатели работы скважин, Wi - веса, регулирующие вклад отдельных показателей работы скважин в функционал; индексы / и

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

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

Для минимизации функционала (4) могут использоваться итерационные градиентные и безградиентные методы. Из градиентных методов наибольшее распространение получили квазиньютоновские, метод сопряженных градиентов, а также специализированные методы для задач наименьших квадратов (Гаусса-Ньютона и Левенберга-Марквардта). Градиентные методы отличаются высокой вычислительной эффективностью, однако требуют доступа к исходному коду решения прямой задачи, что не позволяет использовать их с готовыми гидродинамическими симуляторами. Из безградиентных методов наиболее известны методы Нелдера-Мида, Хука-Дживса, координатная стратегия и прямой поиск, а также группа генетических алгоритмов и метод роя частиц. Метод Нелдера-Мида прост в реализации, при этом показывает неплохие результаты как в плане скорости вычислений, так и в эффективности поиска минимума в сравнении с другими безградиентными процедурами. При этом, как и другие безградиентные методы, его можно использовать без доступа к коду решения прямой задачи, что сделало его популярным для включения в модули автоадаптации ряда коммерческих программных комплексов гидродинамического моделирования. По этой причине в данной работе использован именно метод Нелдера-Мида в сочетании со свободно распространяемым гидродинамическим симулятором OPMFlow [7].

Описание тестовой модели и постановка задачи

Для проведения исследований была подготовлена синтетическая неоднородная модель элемента пятиточечной системы заводнения (рис. 1) с четырьмя добывающими скважинами по углам и одной нагнетательной в центре. Размерность модели 25^25x10 = 6250 ячеек. Размер каждой ячейки 100x100x10 м. Кровля пласта задана на глубине 2400 м, подошва соответственно - на 2500 м. Исходные значения пористости заданы постоянными вдоль стволов скважин на всю толщину пласта. Проницаемость рассчитывается по зависимости от пористости. В табл. 1 приведены заданные значения

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

L,

/

щшшшшшшшшш.

■■■■I мят

INJ1313

ШШШШШШШШШШШШШШШШ

Рис. 1. Элемент пятиточечной системы заводнения

Таблица 1

Значения ФЕС на скважинах «реальной» модели. PROD - добывающие скважины, INJ - нагнетательная

Скважина PROD11 PROD125 PROD2525 PROD251 INJ1313

Положение скважины Юго- западный угол Северозападный угол Северовосточный угол Юго-восточный угол Центр модели

Значение пористости, ед. 0,25 0,18 0,25 0,1 0,25

Значение проницаемости, мД 900 255 900 60 900

На рис. 2 приведены поля пористости и проницаемости «реальной» модели, построенные по данным на скважинах. Поле пористости построено кригингом со следующими параметрами вариограммы: Яг = 1800, Я2 = 500, Я3 = 50 м, ^ = 135о. Проницаемость (в мДарси) задается по зависимости от пористости Кпр = а* ехр(Ь * Кп), где а = 10, Ь = 18.

Рис. 2. «Реальные» поля пористости (слева) и проницаемости (справа)

Элемент расположен в чисто нефтяной зоне. Начальное пластовое давление - 250 бар (1 бар = 10-1 МПа). Другие параметры пласта и флюидов заданы по аналогии с одним из реальных объектов. Рассматриваемый период эксплуатации модели - 15 лет с ежемесячными измерениями показателей, всего 180 моментов замеров. С начала расчета реализуется заводнение. На нагнетательной скважине задан расход 14 400 м3/сут, при забойном давлении не более 400 бар. На добывающих скважинах задан дебит по жидкости 2 700, 5 500, 900 и 5 500м3/сут, соответственно считая от северо-западной к юго-западной скважине по часовой стрелке. Забойное давление для всех добывающих скважин - не ниже 130 бар. На рис. 3 приведены начальное и конечное состояние модели.

Рис. 3. Распределение нефтенасыщенности «реальной» модели на начальном

и конечном временном шаге

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

Для решения обратной задачи рассмотрено два различных начальных приближения. В качестве уточняемых параметров выбраны следующие: радиусы вариограммы в главном и второстепенном направлениях в плоскости напластования и й2, угол поворота главного направления вариограммы ^ и параметры а и Ъ зависимости проницаемости от пористости Кпр = а * ехр(Ь * Кп). В первом варианте начального приближения все параметры, кроме й2, были занижены относительно истинных; во-втором, наоборот, завышены, а й2 занижен (табл. 2).

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

Результаты

Результаты решения обратной задачи по рассматриваемым вариантам иллюстрируют поля распределения ФЕС: «реальное» поле, начальное приближение и результат адаптации (рис. 4, а-в). Также рассчитаны поля дисперсии значений ФЕС, которые отражают неопределенность оцененных значений поля в каждой ячейке модели. Чем дальше от точек с измеренными статическими данными (скважин), тем выше дисперсия [5]. Поля дисперсии позволяют наглядно оценить влияние каждого из параметров вариограммы на итоговое распределение ФЕС (см. рис. 4, в).

Сводные результаты решения обратной задачи для рассмотренных вариантов представлены в табл. 2. На рис. 4, б в качестве наиболее информативных приведены поля проницаемости для всех моделей. Из сравнения рис. 4, б с «реальным» полем

проницаемости (см. рис. 2) хорошо видно, что начальное приближение № 1 (см. рис. 4, а) позволяет намного точнее решить обратную задачу. Здесь в целом корректная конфигурация поля ФЕС восстанавливается во всех 5 случаях, с закономерным улучшением при включении в функционал обоих видов замеряемых показателей. При этом достаточно хорошее приближение к искомому полю проницаемости обеспечивается, даже несмотря на не вполне точное восстановление значений отдельных управляющих параметров (см. табл. 2).

Таблица 2

Результаты адаптации моделей с различными начальными приближениями и распределением весов между показателями в функционале

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

1 0,75 0,5 0,25 0 Относительные веса при дебитах по добывающим скважинам и расходу по нагнетательной скважине

0 0,25 0,5 0,75 1 Относительные веса при забойных давлениях на добывающих скважинах и нагнетательной скважине

Начальное приближение 1 Номер модели: 1-1 1-2 1-3 1-4 1-5 Параметры реальной модели: R1 = 1800 R2 = 500 Phi= 135 a = 10 b = 18

R1 1600 1500 1650 1785 1800 1875

R2 700 930 800 416 404 364

Phi 110 124 127 128 129 117

a 8 10,3 10,2 10,1 9,9 9

b 16 17,8 17,8 18 18 18,2

errortotai: 2,47 1,78 0,17 0,12 1,5

errordebt: 2,47 1,66 0,14 0,08 0

errorpressure: 0 0,12 0,03 0,04 1,5

iters: 75 59 64 56 27

Начальное приближение 2 Номер модели: 2-1 2-2 2-3 2-4 2-5

R1 2000 1860 2015 1175 2100 2070

R2 300 784 795 1000 785 782

Phi 160 157 145 211 188 159

a 12 14,4 11 8,9 9,7 10,5

b 20 16,3 17,4 18,6 18,2 17,7

errortotai: 5,26 2,33 3,59 1,96 0,46

errordebt: 5,26 2,26 2,92 1,34 0

errorpressure: 0 0,07 0,67 0,62 0,46

iters: 56 16 22 62 35

Примечание: R1 - радиус вариограммы в главном направлении, м; R2 - радиус вариограммы в направлении, перпендикулярном главному, м; Phi - угол поворота главного направления, м; а и b - параметры зависимости «пористость-проницаемость», Кпр = а*ехр(Ь*Кп); errordebt - суммарная невязка по дебитам и расходу; errorpressure - суммарная невязка по забойным давлениям; iters -количество итераций алгоритма.

Cell Results: VAR

-0.8

Рис. 4, а. Поля дисперсии (вверху) и проницаемости (внизу) для начального приближения № 1 (слева) и № 2 (справа)

А А Ж шЯ /а

■Им А !2j ■Яш

Cell Results: PERMX

900 800

6Б0

Б00

ЗБО

200

60

Рис. 4, б. Итоговые поля проницаемости для начального приближения № 1 (вверху) и № 2 (внизу) при различном распределении весов в соответствии с табл. 2

Cell Results: VAR

-0.8

Рис. 4, в. Итоговые поля дисперсии для начального приближения № 1 (вверху) и № 2 (внизу) при различном распределении весов в соответствии с табл. 2

При втором начальном приближении алгоритм Нелдера-Мида значительно хуже справляется с решением обратной задачи как по общей конфигурации итогового поля ФЕС (см. рис. 4, б и 4, в), так и по значениям отдельных параметров (см. табл. 2). Неожиданным является тот факт, что учет в функционале невязок по обоим видам показателей не улучшает, а ухудшает получаемое решение. Наилучший результат достигается в случае, когда учитывается невязка только по забойным давлениям на скважинах (вес при дебитах/расходе равен нулю).

На рис. 5 приведены зависимости итоговой суммарной невязки в функционале (общей и по показателям) от весов для двух видов замеряемых показателей работы скважины. Картина для двух начальных приближений существенно различается. Если для начального приближения 1 ожидаемо лучшие результаты дает комбинация обоих видов измерений, то во втором случае - чем больше вес забойных давлений, тем лучше результаты. С одной стороны, это показывает большую чувствительность задачи к невязке забойных давлений. С другой стороны, найденная точка явно далека от оптимальной. Итоговые невязки кратно выше, чем для начального приближения 1, а полученные оценки управляющих параметров и общая конфигурация поля распределения ФЕС далеки от истинных (см. табл. 2 и рис. 4).

Приведенные выводы подтверждаются на рис. 6, а и 6, б, показывающих графики дебитов и давлений для моделей 1-3 и 2-3 (начальные приближения 1 и 2, веса равно распределены между показателями). Здесь сопоставлены динамики показателей работы скважин для «реальной» модели, модели начального приближения и адаптированной модели. Видно, что численные значения невязок хорошо отражают соответствие между

«реальными» замерами и динамиками показателей согласно адаптированной модели. Также подтверждается, что начальная невязка по забойным давлениям в случае начального приближения № 2 намного выше, чем в случае начального приближения № 1.

Рис. 5. Зависимость невязки от распределения весов для начального приближения № 1 (вверху)

и № 2 (внизу)

Рис. 6, a. Графики дебитов нефтяных скважин PROD 11 и PROD 125 для начального приближения

№ 1 (вверху) и № 2 (внизу)

^ Забойное давление на скважине PROD 125

Рис. 6, б. Графики забойных давлений нефтяных скважин PROD11 и PROD125 для начального

приближения № 1 (вверху) и № 2 (внизу)

Для анализа эффективности поиска решения алгоритмом Нелдера-Мида на рис. 7 и 8 приведены зависимости соответственно невязок и изменения уточняемых параметров от номера итерации для моделей 1-3 и 2-3. На рис. 8 параметры нормализованы таким образом, что значение 1 соответствует величине параметра в «реальной» модели. Хорошо видно, что в случае начального приближения № 1 большинство параметров уже к 20-й итерации приближаются к своему искомому значению и далее корректируются незначительно. Иная ситуация наблюдается в случае с начальным приближением № 2. Алгоритм попадает в локальный минимум и выйти из него не может, хотя и невязки, и значения параметров далеки от оптимальных. Данная особенность характеризует проблемы поиска минимума методом Нелдера-Мида: несмотря на то, что он является безградиентным, наблюдается возможность стагнации в неоптимальной точке и сильная чувствительность к начальному приближению.

Модель 1-3

0 20 40 60 80

Количество итераций

Рис. 7. Зависимость невязок от номера итерации алгоритма. Модель с начальным приближением № 1 (вверху) и начальным приближением № 2 (внизу). Веса равно распределены, по 0,5 для дебитов/расхода и забойных давлений

3 2,0 о.

I-

ш

га о. га га с ш 2

1,5 1,0 0,5

ш к

X т о

> о,о

d

-V

20

40

60

Количество итераций

R2 РЫ а Ь

80

Рис. 8. Зависимость уточняемых параметров от номера итерации алгоритма. Модель с начальным приближением № 1 (слева) и начальным приближением № 2 (справа). Веса равно распределены,

по 0,5 для дебитов/расхода и забойных давлений

Отдельно стоит отметить, что найденное значение радиуса вариограммы для всех рассчитанных вариантов значительно отличается от «реального», что говорит о его малом влиянии на показатели работы скважин. Это подтверждают рис. 2 и 4: изменение й2 в рассматриваемой модели сказывается практически только на ширине высокопроницаемой «диагонали», и при его изменении в достаточно широком диапазоне не происходит существенного изменения ФЕС в окрестностях добывающих скважин, а также темпов прорыва воды вдоль первоочередного ЮЗ-СВ направления.

0

Выводы

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

В частности, подтверждены и наглядно показаны:

- возможность корректного воспроизведения главного направления и в целом полей распределения ФЕС с существенно анизотропной конфигурацией;

- влияние относительных весов различных типов замеров в функционале на результат решения обратной задачи;

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

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

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

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

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

С другой стороны, данная особенность лишает безградиентные методы соответствующего преимущества над более эффективными в других аспектах методами гладкой оптимизации. Их применению для решения обратных задач в геологически согласованной постановке посвящены исследования [3, 4, 8].

Статья написана в рамках выполнения государственного задания (тема «Научное обоснование новых экологически чистых технологий разработки месторождений углеводородов в сложных горно-геологических условиях на основе 3D-компьютерных экспериментов», № АААА-А16-116022510270-1).

ЛИТЕРАТУРА

1. Oliver D.S., Reynolds A.C., Liu N. Inverse Theory for Petroleum Reservoir Characterization and History Matching. New York: Cambridge University Press, 2008. 394 p.

2. ЗакировЭ.С. Трехмерные многофазные задачи прогнозирования, анализа и регулирования разработки месторождений нефти и газа. М.: Изд. дом «Грааль», 2001. 302 с.

3. Zakirov E.S., Indrupskiy I.M., Lubimova O.V., Shiriaev I.M. Geostatistically-Consistent History Matching // Proceedings of the 14th European Conference on the Mathematics of Oil Recovery. Catania, Sicily, Italy, 8-11 September 2014. 13 р.

4. Закиров Э.С., Индрупский И.М., Любимова О.В., Ширяев И.М., Аникеев Д.П. Согласованная адаптация геостатистических моделей залежей нефти и газа // Доклады Академии наук. 2017. Т. 476. № 4. С. 421-425.

5. Deutsch C.V. Geostatistical reservoir modeling. New York: Oxford University Press, 2002. 376 p.

6. Дюбрул О. Использование геостатистики для включения в геологическую модель сейсмических данных. EAGE, 2002. 296 с.

7. The Open Porous Media Initiative web site [Электронный ресурс]. - Режим доступа: https://opm-project.org (Дата обращения 13.07.2018).

8. Закиров Э.С., Любимова О.В., Индрупский И.М., Аникеев Д.П., Архипова Е.Ю. Автоматизированная геологически-согласованная адаптация распределения фаций и свойств пласта в межскважинном пространстве на основе сопряженных методов // SPE Russian Petroleum Technology Conference, 16-18 October 2017, Moscow, Russia. SPE-187803-RU. 23 с.